tpc/dEdx/GlobaldEdxFitter.cc

0001 #include "GlobaldEdxFitter.h"
0002
0003 #include "bethe_bloch.h"
0004 #include "TF1.h"
0005 #include "TF2.h"
0006 #include "TF3.h"
0007 #include "TChain.h"
0008 #include "TGraph.h"
0009 #include "Math/Minimizer.h"
0010 #include "Math/Functor.h"
0011 #include "Math/Factory.h"
0012
0013 void GlobaldEdxFitter::processResidualData(const std::string& infile, size_t ntracks, size_t skip)
0014 {
0015   std::unique_ptr<TChain> t = std::make_unique<TChain>();
0016   t->Add((infile+"?#residualtree").c_str());
0017 //  TFile* f = TFile::Open(infile.c_str());
0018 //  TTree* t = (TTree*)f->Get("residualtree");
0019
0020   float px;
0021   float py;
0022   float pz;
0023   float dedx;
0024   float eta;
0025   int nmaps;
0026   int nintt;
0027   int ntpc;
0028   float dcaxy;
0029
0030   t->SetBranchAddress("px",&px);
0031   t->SetBranchAddress("py",&py);
0032   t->SetBranchAddress("pz",&pz);
0033   t->SetBranchAddress("dedx",&dedx);
0034   t->SetBranchAddress("eta",&eta);
0035   t->SetBranchAddress("nmaps",&nmaps);
0036   t->SetBranchAddress("nintt",&nintt);
0037   t->SetBranchAddress("ntpc",&ntpc);
0038   t->SetBranchAddress("dcaxy",&dcaxy);
0039
0040   size_t total_entries = t->GetEntriesFast();
0041
0042   for(size_t entry=skip; entry<(skip+ntracks); entry++)
0043   {
0044     if(entry==total_entries)
0045     {
0046       break;
0047     }
0048     if(entry % 1000 == 0)
0049     {
0050       std::cout << entry << std::endl;
0051     }
0052     t->GetEntry(entry);
0053     if(nmaps>0 && nintt>0 && std::fabs(eta)<1. && dcaxy<0.5 && ntpc>30)
0054     {
0055       p.push_back(std::sqrt(px*px+py*py+pz*pz));
0056       dEdx.push_back(dedx);
0057     }
0058   }
0059   std::cout << "number of good tracks: " << p.size() << std::endl;
0060   //f->Close();
0061 }
0062
0063 void GlobaldEdxFitter::addTrack(double trk_dEdx, double trk_p)
0064 {
0065   dEdx.push_back(trk_dEdx);
0066   p.push_back(trk_p);
0067 }
0068
0069 double GlobaldEdxFitter::get_fitquality_new(double A)
0070 {
0071   //double chi2 = 0.;
0072   //double ndf = -1.;
0073
0074   double pi_chi2 = 0.;
0075   double K_chi2 = 0.;
0076   double p_chi2 = 0.;
0077   double d_chi2 = 0.;
0078   double pi_ndf = -1.;
0079   double K_ndf = -1.;
0080   double p_ndf = -1.;
0081   double d_ndf = -1.;
0082
0083   for(size_t i=0; i<dEdx.size(); i++)
0084   {
0085     const double dedx_pi = bethe_bloch_new_1D(p[i]/dedx_constants::m_pi,A);
0086     const double dedx_K = bethe_bloch_new_1D(p[i]/dedx_constants::m_K,A);
0087     const double dedx_p = bethe_bloch_new_1D(p[i]/dedx_constants::m_p,A);
0088     const double dedx_d = bethe_bloch_new_1D(p[i]/dedx_constants::m_d,A);
0089
0090     //std::cout << "dedx: (" << dedx_pi << ", " << dedx_K << ", " << dedx_p << ", " << dedx_d << ")" << std::endl;
0091     //std::cout << "measured: " << dEdx[i] << std::endl;
0092
0093     const double pi_dist = fabs(dEdx[i]-dedx_pi);
0094     const double K_dist = fabs(dEdx[i]-dedx_K);
0095     const double p_dist = fabs(dEdx[i]-dedx_p);
0096     const double d_dist = fabs(dEdx[i]-dedx_d);
0097
0098     //std::cout << "dist: (" << pi_dist << ", " << K_dist << ", " << p_dist << ", " << d_dist << ")" << std::endl;
0099
0100     if(pi_dist<K_dist && pi_dist<p_dist && pi_dist<d_dist)
0101     {
0102       //std::cout << "pion" << std::endl;
0103       //chi2 += pi_dist*pi_dist;
0104       pi_chi2 += pi_dist*pi_dist;
0105       pi_ndf += 1.;
0106     }
0107     else if(K_dist<pi_dist && K_dist<p_dist && K_dist<d_dist)
0108     {
0109       //std::cout << "K" << std::endl;
0110       //chi2 += K_dist*K_dist;
0111       K_chi2 += K_dist*K_dist;
0112       K_ndf += 1.;
0113     }
0114     else if(p_dist<pi_dist && p_dist<K_dist && p_dist<d_dist)
0115     {
0116       //std::cout << "p" << std::endl;
0117       //chi2 += p_dist*p_dist;
0118       p_chi2 += p_dist*p_dist;
0119       p_ndf += 1.;
0120     }
0121
0122     else if(d_dist<pi_dist && d_dist<K_dist && d_dist<p_dist)
0123     {
0124       //std::cout << "d" << std::endl;
0125       //chi2 += d_dist*d_dist;
0126       d_chi2 += d_dist*d_dist;
0127       d_ndf += 1.;
0128     }
0129
0130     //ndf += 1.;
0131     //std::cout << "chi2: " << chi2 << std::endl;
0132     //std::cout << "ndf: " << ndf << std::endl;
0133   }
0134
0135   if(pi_ndf<1.) {pi_chi2 = 0.; pi_ndf = 1.;}
0136   if(K_ndf<1.) {K_chi2 = 0.; K_ndf = 1.;}
0137   if(p_ndf<1.) {p_chi2 = 0.; p_ndf = 1.;}
0138   if(d_ndf<1.) {d_chi2 = 0.; d_ndf = 1.;}
0139
0140   const double quality = sqrt((pi_chi2*pi_chi2)/(pi_ndf*pi_ndf)+(K_chi2*K_chi2)/(K_ndf*K_ndf)+(p_chi2*p_chi2)/(p_ndf*p_ndf)+(d_chi2*d_chi2)/(d_ndf*d_ndf));
0141   //const double quality = chi2/ndf;
0142
0143   //std::cout << "A: " << A << " quality: " << quality << std::endl;
0144
0145   return quality;
0146 }
0147
0148 double GlobaldEdxFitter::get_fitquality(double norm, double ZS_loss)
0149 {
0150   float pi_chi2 = 0.;
0151   float K_chi2 = 0.;
0152   float p_chi2 = 0.;
0153   float d_chi2 = 0.;
0154   float pi_ndf = -1.;
0155   float K_ndf = -1.;
0156   float p_ndf = -1.;
0157   float d_ndf = -1.;
0158   //float ndf = -1;
0159
0160   for(size_t i=0; i<dEdx.size(); i++)
0161   {
0162
0163     const float dedx_pi  = norm*bethe_bloch_total(p[i]/dedx_constants::m_pi) - ZS_loss;
0164     const float dedx_K = norm*bethe_bloch_total(p[i]/dedx_constants::m_K) - ZS_loss;
0165     const float dedx_p = norm*bethe_bloch_total(p[i]/dedx_constants::m_p) - ZS_loss;
0166     const float dedx_d = norm*bethe_bloch_total(p[i]/dedx_constants::m_d) - ZS_loss;
0167
0168     const float pi_dist = fabs(dedx_pi-dEdx[i]);
0169     const float K_dist = fabs(dedx_K-dEdx[i]);
0170     const float p_dist = fabs(dedx_p-dEdx[i]);
0171     const float d_dist = fabs(dedx_d-dEdx[i]);
0172
0173     if(pi_dist<K_dist && pi_dist<p_dist && pi_dist<d_dist)
0174     {
0175       pi_chi2 += pi_dist*pi_dist/std::fabs(dedx_pi);
0176       pi_ndf += 2.;
0177     }
0178     else if(K_dist<pi_dist && K_dist<p_dist && K_dist<d_dist)
0179     {
0180       K_chi2 += K_dist*K_dist/std::fabs(dedx_K);
0181       K_ndf += 2.;
0182     }
0183     else if(p_dist<pi_dist && p_dist<K_dist && p_dist<d_dist)
0184     {
0185       p_chi2 += p_dist*p_dist/std::fabs(dedx_p);
0186       p_ndf += 2.;
0187     }
0188     else if(d_dist<pi_dist && d_dist<K_dist && d_dist<p_dist)
0189     {
0190       d_chi2 += d_dist*d_dist/std::fabs(dedx_d);
0191       d_ndf += 2.;
0192     }
0193     //ndf += 2;
0194   }
0195
0196   const double quality = std::sqrt((pi_chi2*pi_chi2)/(pi_ndf*pi_ndf) +
0197                         (K_chi2*K_chi2)/(K_ndf*K_ndf) +
0198                         (p_chi2*p_chi2)/(p_ndf*p_ndf) +
0199                         (d_chi2*d_chi2)/(d_ndf*d_ndf));
0200
0201 /*
0202   const double quality = sqrt((pi_chi2*pi_chi2 +
0203                         K_chi2*K_chi2 +
0204                         p_chi2*p_chi2)
0205                         /(ndf*ndf));
0206 */
0207   //std::cout << "norm: " << norm << " ZS: " << ZS_loss << std::endl;
0208   //std::cout << "quality: " << quality << std::endl;
0209
0210   return quality;
0211 }
0212
0213 std::vector<double> GlobaldEdxFitter::get_betagamma(double A)
0214 {
0215   std::vector<double> betagamma;
0216   for(size_t i=0; i<dEdx.size(); i++)
0217   {
0218     const double dedx_pi = bethe_bloch_new_1D(p[i]/dedx_constants::m_pi,A);
0219     const double dedx_K = bethe_bloch_new_1D(p[i]/dedx_constants::m_K,A);
0220     const double dedx_p = bethe_bloch_new_1D(p[i]/dedx_constants::m_p,A);
0221     const double dedx_d = bethe_bloch_new_1D(p[i]/dedx_constants::m_d,A);
0222
0223     //std::cout << "dedx: (" << dedx_pi << ", " << dedx_K << ", " << dedx_p << ", " << dedx_d << ")" << std::endl;
0224     //std::cout << "measured: " << dEdx[i] << std::endl;
0225
0226     const double pi_dist = fabs(dEdx[i]-dedx_pi);
0227     const double K_dist = fabs(dEdx[i]-dedx_K);
0228     const double p_dist = fabs(dEdx[i]-dedx_p);
0229     const double d_dist = fabs(dEdx[i]-dedx_d);
0230
0231     if(pi_dist<K_dist && pi_dist<p_dist && pi_dist<d_dist)
0232     {
0233       betagamma.push_back(p[i]/dedx_constants::m_pi);
0234     }
0235     else if(K_dist<pi_dist && K_dist<p_dist && K_dist<d_dist)
0236     {
0237       betagamma.push_back(p[i]/dedx_constants::m_K);
0238     }
0239     else if(p_dist<pi_dist && p_dist<K_dist && p_dist<d_dist)
0240     {
0241       betagamma.push_back(p[i]/dedx_constants::m_p);
0242     }
0243     else if(d_dist<pi_dist && d_dist<K_dist && d_dist<p_dist)
0244     {
0245       betagamma.push_back(p[i]/dedx_constants::m_d);
0246     }
0247   }
0248   return betagamma;
0249 }
0250
0251 double GlobaldEdxFitter::get_fitquality_functor(const double* x)
0252 {
0253   return get_fitquality_new(x[0]);
0254 }
0255
0256 double GlobaldEdxFitter::get_fitquality_wrapper(double* x, double* /* par */)
0257 {
0258   return get_fitquality(x[0]);
0259 }
0260
0261 double GlobaldEdxFitter::get_fitquality_wrapper_ZS(double* x, double* /* par */)
0262 {
0263   return get_fitquality(x[0],x[1]);
0264 }
0265
0266 double GlobaldEdxFitter::get_fitquality_wrapper_new(double* x, double* /* par */)
0267 {
0268   return get_fitquality_new(x[0]);
0269 }
0270
0271 TF3* GlobaldEdxFitter::create_TF3_new(const std::string& name)
0272 {
0273   TF3* f = new TF3(name.c_str(),this,&GlobaldEdxFitter::get_fitquality_wrapper_new,min_norm,max_norm,min_B,max_B,min_ZS,max_ZS,0,"GlobaldEdxFitter","get_fitquality_wrapper_new");
0274   return f;
0275 }
0276
0277 TF1* GlobaldEdxFitter::create_TF1(const std::string& name)
0278 {
0279   TF1* f = new TF1(name.c_str(),this,&GlobaldEdxFitter::get_fitquality_wrapper_new,min_norm,max_norm,0,"GlobaldEdxFitter","get_fitquality_wrapper");
0280   return f;
0281 }
0282
0283 TF2* GlobaldEdxFitter::create_TF2(const std::string& name)
0284 {
0285   TF2* f = new TF2(name.c_str(),this,&GlobaldEdxFitter::get_fitquality_wrapper_ZS,min_norm,max_norm,min_ZS,max_ZS,0,"GlobaldEdxFitter","get_fitquality_wrapper_ZS");
0286   return f;
0287 }
0288
0289 double GlobaldEdxFitter::get_minimum_new()
0290 {
0291 /*
0292   TF3* f = create_TF3_new("temp");
0293   double minA;
0294   double minB;
0295   double minC;
0296   f->GetMinimumXYZ(minA,minB,minC);
0297   delete f;
0298   return std::make_tuple(minA,minB,minC);
0299 */
0300   ROOT::Math::Minimizer* minimizer = ROOT::Math::Factory::CreateMinimizer("Minuit2");
0301   minimizer->SetMaxFunctionCalls(1000000);
0302   minimizer->SetMaxIterations(10000);
0303   minimizer->SetTolerance(0.1);
0304   minimizer->SetPrintLevel(1);
0305   ROOT::Math::Functor f(this,&GlobaldEdxFitter::get_fitquality_functor,1);
0306   double step[1] = {.01};
0307   double variable[1] = {20.};
0308   minimizer->SetFunction(f);
0309   minimizer->SetVariable(0,"A",variable[0],step[0]);
0310   minimizer->Minimize();
0311   const double *xs = minimizer->X();
0312   delete minimizer;
0313   return xs[0];
0314 }
0315
0316 double GlobaldEdxFitter::get_minimum()
0317 {
0318   TF1* f = create_TF1("temp");
0319   f->SetNpx(1000);
0320   double minX = f->GetMinimumX();
0321   delete f;
0322   return minX;
0323 }
0324
0325 std::pair<double,double> GlobaldEdxFitter::get_minimum_ZS()
0326 {
0327   TF2* f = create_TF2("temp");
0328   double minX;
0329   double minY;
0330   f->GetMinimumXY(minX,minY);
0331   delete f;
0332   return std::make_pair(minX,minY);
0333 }
0334
0335 TGraph* GlobaldEdxFitter::graph_vsbetagamma(double A)
0336 {
0337   std::vector<double> betagamma = get_betagamma(A);
0338   TGraph* g = new TGraph(dEdx.size(),betagamma.data(),dEdx.data());
0339   return g;
0340 }
0341
0342 TGraph* GlobaldEdxFitter::graph_vsp()
0343 {
0344   TGraph* g = new TGraph(dEdx.size(),p.data(),dEdx.data());
0345   return g;
0346 }