LightFlavorRatios/single_track_abundance/get_singletrack_abundance.C

0001 #include "TFile.h"
0002 #include "TTree.h"
0003
0004 #include "../util/HistogramContainer.h"
0005
0006 void get_singletrack_abundance(std::string infile = "/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/data/nTP/nTuplizer_tracks_00053877_00001_00000.root", std::string outfile="output/test_out.root")
0007 {
0008   const float dedx_relative_diff_cutoff = 0.4;
0009
0010   TFile* f = TFile::Open(infile.c_str());
0011   TTree* t = (TTree*)f->Get("ntp_track");
0012
0013   HistogramContainer h_pi("pi","#pi^{+}+#pi^{-}");
0014   HistogramContainer h_piplus("piplus","#pi^{+}");
0015   HistogramContainer h_piminus("piminus","#pi^{-}");
0016   HistogramContainer h_K("K","K^{+} + K^{-}");
0017   HistogramContainer h_Kplus("Kplus","K^{+}");
0018   HistogramContainer h_Kminus("Kminus","K^{-}");
0019   HistogramContainer h_p("p","p^{+} + p^{-}");
0020   HistogramContainer h_pplus("pplus","p^{+}");
0021   HistogramContainer h_pminus("pminus","p^{-}");
0022   HistogramContainer h_track("track","all tracks");
0023   HistogramContainer h_trackplus("trackplus","positively-charged track");
0024   HistogramContainer h_trackminus("trackminus","negatively-charged track");
0025
0026   float pt;
0027   float eta;
0028   float phi;
0029   float ntracks;
0030   float charge;
0031   float measured_dedx;
0032   float expected_pi_dedx;
0033   float expected_K_dedx;
0034   float expected_p_dedx;
0035
0036   t->SetBranchAddress("pt",&pt);
0037   t->SetBranchAddress("eta",&eta);
0038   t->SetBranchAddress("phi",&phi);
0039   t->SetBranchAddress("ntrk",&ntracks);
0040   t->SetBranchAddress("charge",&charge);
0041   t->SetBranchAddress("dedx",&measured_dedx);
0042   t->SetBranchAddress("pidedx",&expected_pi_dedx);
0043   t->SetBranchAddress("kdedx",&expected_K_dedx);
0044   t->SetBranchAddress("prdedx",&expected_p_dedx);
0045
0046   for(size_t entry=0; entry<t->GetEntries(); ++entry)
0047   {
0048     if(entry % 10000 == 0) std::cout << entry << std::endl;
0049     t->GetEntry(entry);
0050
0051     h_track.Fill(pt,eta,phi,ntracks);
0052     if(charge>0.)
0053     {
0054       h_trackplus.Fill(pt,eta,phi,ntracks);
0055     }
0056     else if(charge<0.)
0057     {
0058       h_trackminus.Fill(pt,eta,phi,ntracks);
0059     }
0060
0061     float pi_diff = fabs(measured_dedx-expected_pi_dedx)/expected_pi_dedx;
0062     float K_diff = fabs(measured_dedx-expected_K_dedx)/expected_K_dedx;
0063     float p_diff = fabs(measured_dedx-expected_p_dedx)/expected_p_dedx;
0064
0065     if(pi_diff<dedx_relative_diff_cutoff && pi_diff<K_diff && pi_diff<p_diff)
0066     {
0067       h_pi.Fill(pt,eta,phi,ntracks);
0068       if(charge>0.)
0069       {
0070         h_piplus.Fill(pt,eta,phi,ntracks);
0071       }
0072       else if(charge<0.)
0073       {
0074         h_piminus.Fill(pt,eta,phi,ntracks);
0075       }
0076     }
0077     else if(K_diff<dedx_relative_diff_cutoff && K_diff<pi_diff && K_diff<p_diff)
0078     {
0079       h_K.Fill(pt,eta,phi,ntracks);
0080       if(charge>0.)
0081       {
0082         h_Kplus.Fill(pt,eta,phi,ntracks);
0083       }
0084       else if(charge<0.)
0085       {
0086         h_Kminus.Fill(pt,eta,phi,ntracks);
0087       }
0088     }
0089     else if(p_diff<dedx_relative_diff_cutoff && p_diff<pi_diff && p_diff<K_diff)
0090     {
0091       h_p.Fill(pt,eta,phi,ntracks);
0092       if(charge>0.)
0093       {
0094         h_pplus.Fill(pt,eta,phi,ntracks);
0095       }
0096       else if(charge<0.)
0097       {
0098         h_pminus.Fill(pt,eta,phi,ntracks);
0099       }
0100     }
0101   }
0102
0103   TFile* outf = new TFile(outfile.c_str(),"RECREATE");
0104   h_pi.Write();
0105   h_piplus.Write();
0106   h_piminus.Write();
0107   h_K.Write();
0108   h_Kplus.Write();
0109   h_Kminus.Write();
0110   h_p.Write();
0111   h_pplus.Write();
0112   h_pminus.Write();
0113   h_track.Write();
0114   h_trackplus.Write();
0115   h_trackminus.Write();
0116   std::cout << "Finished" << std::endl;
0117 }