analysis/FrawleyTrackEvalMacros/quarkonia_reconstruction_embedded.C

0001 #include <TChain.h>
0002 #include <TCanvas.h>
0003 #include <TH1.h>
0004 #include <TH1F.h>
0005 #include <TH2D.h>
0006 #include <TF1.h>
0007 #include <TFile.h>
0008 #include <TGraph.h>
0009 #include <TStyle.h>
0010 #include <TROOT.h>
0011 #include <TLegend.h>
0012 #include <TLine.h>
0013 #include <TLatex.h>
0014 #include <TRandom1.h>
0015 #include <TPolyLine.h>
0016 #include <iostream>
0017 #include <fstream>
0018 #include <TMath.h>
0019 //#include <iomanip.h>
0020 #include <TLorentzVector.h>
0021
0022 //#define TEST
0023
0024 void quarkonia_reconstruction_embedded()
0025 {
0026   gROOT->SetStyle("Plain");
0027   gStyle->SetOptStat(1);
0028   gStyle->SetOptTitle(1);
0029
0030   int verbose = 0;
0031   int embed_flag = 3;  // embedding flag used during Upsilon generation
0032
0033   // track cuts
0034   double quality_cut = 3.0;
0035   double dca_cut = 0.1;
0036
0037   char lepton[100];
0038   sprintf(lepton,"electron");
0039
0040   double decaymass=0.000511;
0041   cout << "Assuming decay particle mass is " << decaymass << endl;
0042
0043   // Define some histograms
0044
0045   int nbpt = 20;
0046   double ptmax = 10.0;
0047
0048   TH1F* hrquality = new TH1F("hrquality", "Reconstructed track Quality", 1000, 0, 10);
0049   TH1F* hrdca2d = new TH1F("hrdca2d", "Reconstructed track dca2d", 1000, 0, 0.05);
0050   TH1F* hrpt = new TH1F("hrpt"," pT", nbpt, 0.0, ptmax);
0051   TH1F* hgpt = new TH1F("hgpt","g4 pT", nbpt, 0.0, ptmax);
0052
0053   TH1D *g4mass = new TH1D("g4mass","G4 input invariant mass",200,7.0,11.0);
0054   g4mass->GetXaxis()->SetTitle("invariant mass (GeV/c^{2})");
0055   TH1D *g4mass_primary = new TH1D("g4mass_primary","G4 input invariant mass",200,7.0,11.0);
0056   g4mass_primary->GetXaxis()->SetTitle("invariant mass (GeV/c^{2})");
0057
0058   TH1D *recomass = new TH1D("recomass","Reconstructed invariant mass",200,7.0,11.0);
0059   recomass->GetXaxis()->SetTitle("invariant mass (GeV/c^{2})");
0060   TH1D *recomass_primary = new TH1D("recomass_primary","Reconstructed invariant mass",200,7.0,11.0);
0061   recomass_primary->GetXaxis()->SetTitle("invariant mass (GeV/c^{2})");
0062
0063   int ups_state = 1;   // used in naming of output files
0064   //int ups_state = 2;   // used in naming of output files
0065   //int ups_state = 3;   // used in naming of output files
0066
0067   // Number of upsilons to process - quit after this number is reached
0068   int nups_requested = 0;
0069   double pair_acc = 0.60 * 0.9 * 0.98;  // fraction into 1 unit of rapidity * pair eID fraction * trigger fraction
0070   // These values are for 1 year of pp running
0071   if(ups_state == 1)  nups_requested = pair_acc * (8769 * 1.126) / (0.38 * 0.9 * 0.98);  // pair_acc * yield for 197 /pb (= 15590)
0072   if(ups_state == 2)  nups_requested = pair_acc * (2205 * 1.126) / (0.38 * 0.9 * 0.98); //  pair_acc * yield for 197 /pb (= 3920)
0073   if(ups_state == 3)  nups_requested = pair_acc * (1156 * 1.126) / (0.38 * 0.9 * 0.98); //  pair_acc * yield for 197 /pb (=2055)
0074
0075   nups_requested = 100000;  // take them all
0076
0077   cout << "ups_state = " << ups_state << " Upsilons requested = " << nups_requested << endl;
0078
0079   int nrecog4mass = 0;
0080   int nrecormass = 0;
0081
0082   // Open the g4 evaluator output file
0083   cout << "Reading electron ntuples " << endl;
0084   int nevents = 0;
0085   double nhittpcin_cum = 0;
0086   double nhittpcin_wt = 0;
0087
0088     // The condor job output files -  we open them one at a time and process them
0089   for(int i=0;i<2000;i++)
0090     {
0091       if(nrecormass > nups_requested)
0092     {
0093       cout << "Reached requested number of reco Upsilons, quit out of file loop" << endl;
0094       break;
0095     }
0096
0097       char name[500];
0098
0099       //sprintf(name,"/sphenix/user/frawley/macros_newTPC_june6/macros/macros/g4simulations/intt_6layers_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",i);
0100       //sprintf(name,"/sphenix/user/frawley/macros_newTPC_june6/macros/macros/g4simulations/intt_4layers_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",i);
0101       //sprintf(name,"/sphenix/user/frawley/macros_newTPC_june6/macros/macros/g4simulations/intt_8layers_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",i);
0102       sprintf(name,"/sphenix/user/frawley/macros_newTPC_june6/macros/macros/g4simulations/intt_0layers_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",i);
0103
0104
0105       /*
0106       // ups states, 100 poins, 80 ns
0107       if(ups_state == 1)
0108     {
0109       //sprintf(name,"/sphenix/user/frawley/fresh_mar8_testing/macros/macros/g4simulations/eval_output_ups1s_100pions_80ns/g4svtx_eval_%i.root_g4svtx_eval.root_primary_eval.root",i);
0110       //sprintf(name,"/sphenix/user/frawley/fresh_mar8_testing/macros/macros/g4simulations/eval_output_ups1s_100pions_massres_121_30kevts/g4svtx_eval_%i.root_g4svtx_eval.root_primary_eval.root",i);
0111     }
0112       else if(ups_state == 2)
0113     {
0114       //sprintf(name,"/sphenix/user/frawley/fresh_mar8_testing/macros/macros/g4simulations/eval_output_ups2s_100pions_80ns/g4svtx_eval_%i.root_g4svtx_eval.root_primary_eval.root",i);
0115       sprintf(name,"/sphenix/user/frawley/fresh_mar8_testing/macros/macros/g4simulations/eval_output_ups2s_100pions_massres_121_20kevts/g4svtx_eval_%i.root_g4svtx_eval.root_primary_eval.root",i);
0116     }
0117       else if(ups_state == 3)
0118     {
0119       //sprintf(name,"/sphenix/user/frawley/fresh_mar8_testing/macros/macros/g4simulations/eval_output_ups3s_100pions_80ns/g4svtx_eval_%i.root_g4svtx_eval.root_primary_eval.root",i);
0120       sprintf(name,"/sphenix/user/frawley/fresh_mar8_testing/macros/macros/g4simulations/eval_output_ups3s_100pions_massres_121_20kevts/g4svtx_eval_%i.root_g4svtx_eval.root_primary_eval.root",i);
0121     }
0122       else
0123     {
0124       cout << "Oops!" << endl;
0125       exit(1);
0126     }
0127       */
0128
0129
0130       cout << "Adding file " << name << endl;
0131
0132       TChain* ntp_track = new TChain("ntp_track","reco tracks");
0133       TChain* ntp_gtrack = new TChain("ntp_gtrack","g4 tracks");
0134       TChain* ntp_vertex = new TChain("ntp_vertex","events");
0135       TChain *ntp_cluster = new TChain("ntp_cluster","clusters");
0136       TChain *ntp_g4hit = new TChain("ntp_g4hit","G4 hits");
0137
0138       ntp_vertex->Add(name);
0139       ntp_track->Add(name);
0140       ntp_gtrack->Add(name);
0141
0142
0143       cout << "The ntuples contain " << ntp_vertex->GetEntries() << " events " << endl;
0144
0145
0146       // Ntuple access variables
0147       // This include file contains the definitions of the ntuple variables
0148 #include "ntuple_variables.C"
0149
0150       //=======================
0151       // Loop over events
0152       //=======================
0153
0154       ntracks = ntp_track->GetEntries();
0155       ngtracks = ntp_gtrack->GetEntries();
0156
0157       int nr = 0;
0158       int ng = 0;
0159       //int nev = 1;
0160       int nev = ntp_vertex->GetEntries();
0161
0162       for(int iev=0;iev<nev;iev++)
0163     {
0164       // drop out when the requested number of reco'd Upsilons has been reached
0165       if(nrecormass > nups_requested)
0166         {
0167           cout << "Reached requested number of reco Upsilons, quit" << endl;
0168
0169           break;
0170         }
0171
0172       nevents++;
0173
0174       int recoget = ntp_vertex->GetEntry(iev);
0175
0176       nhittpcin_cum += nhittpcin;
0177       nhittpcin_wt += 1.0;
0178
0179       if(verbose)
0180         cout
0181           << "iev " << iev
0182           << " event " << event
0183           << " ntracks  (reco) " << ntracks
0184           << " ngtracks (g4) " << ngtracks
0185           //<< " gvz " << egvz
0186           //     << " vz " << evz
0187           << endl;
0188
0189
0190       //============================
0191       // process G4 tracks
0192       // for this event
0193       //============================
0194
0195       int ng4trevt_elec = -1;
0196       int ng4trevt_pos = -1;
0197       int g4trnum_elec[1000];
0198       int g4trnum_pos[1000];
0199
0200
0201
0202       for(int ig=ng;ig<ng+ngtracks;ig++)
0203         {
0204           int recoget1 = ntp_gtrack->GetEntry(ig);
0205           if(!recoget1)
0206         {
0207           if(verbose > 0) cout << "Did not get entry for ig = " << ig << endl;
0208           break;
0209         }
0210
0211
0212           // This bookkeeping is needed because the evaluator records for each event the total track count in ntp_vertex, even
0213           // when it writes out only embedded tracks
0214           if(tevent != iev)
0215         {
0216           if(verbose > 0) cout << " reached new event for ig = " << ig << " tevent = " << tevent << endl;
0217           ng = ig;
0218           break;
0219         }
0220           if(ig == ng+ngtracks - 1)
0221         {
0222           if(verbose > 0) cout << " last time through loop for ig = " << ig << " revent = " << tevent << endl;
0223           ng = ig+1;
0224         }
0225
0226
0227           if(tembed != embed_flag)
0228         //if( tgtrackid != 1 && tgtrackid != 2)
0229         continue;
0230
0231
0232           // we want only electrons or positrons
0233           if(tflavor != 11 && tflavor != -11)
0234         continue;
0235
0236           if(tflavor == 11)
0237         {
0238           // electron
0239           ng4trevt_elec++;
0240           if(ng4trevt_elec > 999)
0241             continue;
0242
0243           if(verbose)
0244             cout << " Found electron:" << endl
0245              << "  ig " << ig
0246              << " ng4trevt_elec " << ng4trevt_elec
0247              << " gtrackID " << tgtrackid
0248              << " gflavor " << tflavor
0249              << " tpx " << tpx
0250              << " tpy " << tpy
0251              << " tpz " << tpz
0252              << endl;
0253
0254           g4trnum_elec[ng4trevt_elec] = ig;
0255         }
0256           else
0257         {
0258           // positron
0259           ng4trevt_pos++;
0260           if(ng4trevt_pos > 999)
0261             continue;
0262
0263           if(verbose)
0264             cout << " Found positron:" << endl
0265              << "  ig " << ig
0266              << " ng4trevt_pos " << ng4trevt_pos
0267              << " gtrackID " << tgtrackid
0268              << " gflavor " << tflavor
0269              << " tpx " << tpx
0270              << " tpy " << tpy
0271              << " tpz " << tpz
0272              << endl;
0273
0274           g4trnum_pos[ng4trevt_pos] = ig;
0275         }
0276         }
0277       ng4trevt_elec++;
0278       ng4trevt_pos++;
0279
0280       if(verbose)
0281         cout << "For this event found " << ng4trevt_elec << " g4 electrons and " << ng4trevt_pos << " g4 positrons"  << endl;
0282
0283       // make all pairs of g4 (truth) electrons and positrons
0284       for(int ielec=0;ielec<ng4trevt_elec;ielec++)
0285         {
0286           int recoelec = ntp_gtrack->GetEntry(g4trnum_elec[ielec]);
0287
0288
0289           if(tembed != embed_flag)
0290         continue;
0291
0292           double elec_pT = sqrt(tpx*tpx+tpy*tpy);
0293           double elec_eta = asinh(tpz/sqrt(tpx*tpx+tpy*tpy));
0294
0295           int gtrid1 = tgtrackid;
0296
0297           TLorentzVector t1;
0298           double E1 = sqrt( pow(tpx,2) + pow(tpy,2) + pow(tpz,2) + pow(decaymass,2));
0299           t1.SetPxPyPzE(tpx,tpy,tpz,E1);
0300
0301           // print out track details
0302           if(verbose > 1)
0303         cout << "  Pair electron:  iev " << iev << " ielec " << ielec
0304              << " g4trnum_elec " << g4trnum_elec[ielec]
0305              << " tgtrackid " << tgtrackid
0306              << " tflavor " << tflavor
0307              << " tpx " << tpx
0308              << " tpy " << tpy
0309              << " tpz " << tpz
0310              << " elec_eta " << elec_eta
0311              << " elec_gpT " << elec_pT
0312              << endl;
0313
0314           for(int ipos =0;ipos<ng4trevt_pos;ipos++)
0315         {
0316           int recopos = ntp_gtrack->GetEntry(g4trnum_pos[ipos]);
0317
0318           int gtrid2 = tgtrackid;
0319
0320           double pos_pT = sqrt(tpx*tpx+tpy*tpy);
0321           double pos_eta = asinh(tpz/sqrt(tpx*tpx+tpy*tpy));
0322
0323           // print out track details
0324           if(verbose > 1)
0325             cout << "  Pair positron: iev " << iev << " ipos " << ipos
0326              << " g4trnum_pos " << g4trnum_pos[ipos]
0327              << " tgtrackid " << tgtrackid
0328              << " tflavor " << tflavor
0329              << " tpx " << tpx
0330              << " tpy " << tpy
0331              << " tpz " << tpz
0332              << " pos_eta " << pos_eta
0333              << " pos_gpT " << pos_pT
0334              << endl;
0335
0336           // Make G4 invariant mass
0337
0338           TLorentzVector t2;
0339           double E2 = sqrt( pow(tpx,2) + pow(tpy,2) + pow(tpz,2) + pow(decaymass,2));
0340           t2.SetPxPyPzE(tpx,tpy,tpz,E2);
0341
0342           TLorentzVector t = t1+t2;
0343
0344           if(verbose)
0345             cout << "                       reco'd g4 mass = " << t.M() << endl << endl;
0346
0347           if(t.M() > 7.0 && t.M() < 11.0)
0348             {
0349               nrecog4mass++;
0350               g4mass->Fill(t.M());
0351               hgpt->Fill(t.Pt());
0352
0353               // Capture the mass spectrum where both tracks are the primary Upsilon decay electrons
0354               g4mass_primary->Fill(t.M());
0355
0356             }
0357         }  // end of ipos loop
0358         } // end of ielec loop
0359
0360
0361       if(verbose)
0362         {
0363           cout << " # of g4 electron tracks = " << ng4trevt_elec
0364            << " # of g4 positron tracks = " << ng4trevt_pos << endl;
0365         }
0366
0367
0368       //=============================
0369       // process reconstructed tracks
0370       // for this event
0371       //=============================
0372
0373       int nrtrevt = 0;
0374       int nr_elec = -1;
0375       int nr_pos = -1;
0376       int rectrnum_elec[1000];
0377       int rectrnum_pos[1000];
0378
0379       //cout << "Number of ntp_track entries = " << recoget << endl;
0380
0381       for(int ir=nr;ir<nr+ntracks;ir++)
0382         {
0383           int recoget = ntp_track->GetEntry(ir);
0384           if(!recoget)
0385         {
0386           if(verbose > 0) cout << "Did not get entry for ir = " << ir << endl;
0387           break;
0388         }
0389
0390
0391           // This bookkeeping is needed because the evaluator records for each event the total track count in ntp_vertex, even
0392           // when it writes out only embedded tracks
0393           if(revent != iev)
0394         {
0395           if(verbose > 1) cout << " reached new event for ir = " << ir << " revent = " << revent << endl;
0396           nr = ir;
0397           break;
0398         }
0399           if(ir == nr+ntracks - 1)
0400         {
0401           if(verbose > 0)  cout << " last time through loop for ir = " << ir << " revent = " << revent << endl;
0402           nr = ir+1;
0403         }
0404
0405           if(rgembed != embed_flag)
0406         continue;
0407
0408           hrquality->Fill(rquality);
0409           hrdca2d->Fill(rdca2d);
0410
0411           // track quality cuts
0412           if(rquality > quality_cut || fabs(rdca2d) > dca_cut)
0413         continue;
0414
0415           // make a list of electrons and positrons
0416           if(rcharge == -1)
0417         {
0418           nr_elec++;
0419           rectrnum_elec[nr_elec] = ir;
0420         }
0421
0422           if(rcharge == 1)
0423         {
0424           nr_pos++;
0425           rectrnum_pos[nr_pos] = ir;
0426         }
0427
0428
0429           double rpT = sqrt(rpx*rpx+rpy*rpy);
0430           double reta = asinh(rpz/sqrt(rpx*rpx+rpy*rpy));
0431
0432           // print out track details
0433           if(verbose)
0434         cout << "     revent " << revent << " ir " << ir
0435              << " rgtrackid " << rgtrackid
0436              << " rgflavor " << rgflavor
0437              << " rvz " << rvz
0438              << " reta " << reta
0439              << " rpT " << rpT
0440              << endl;
0441         }
0442
0443       nr_elec++;
0444       nr_pos++;
0445
0446       for(int ielec = 0;ielec<nr_elec;ielec++)
0447         {
0448
0449           TLorentzVector t1;
0450
0451           int recoget1 = ntp_track->GetEntry(rectrnum_elec[ielec]);
0452
0453           int trid1 = rgtrackid;
0454
0455           double E1 = sqrt( pow(rpx,2) + pow(rpy,2) + pow(rpz,2) + pow(decaymass,2));
0456           t1.SetPxPyPzE(rpx,rpy,rpz,E1);
0457
0458           for(int ipos = 0;ipos<nr_pos;ipos++)
0459         {
0460           int recoget2 = ntp_track->GetEntry(rectrnum_pos[ipos]);
0461
0462           int trid2 = rgtrackid;
0463
0464           TLorentzVector t2;
0465           double E2 = sqrt( pow(rpx,2) + pow(rpy,2) + pow(rpz,2) + pow(decaymass,2));
0466
0467           t2.SetPxPyPzE(rpx,rpy,rpz,E2);
0468
0469           TLorentzVector t = t1+t2;
0470
0471           if(verbose)
0472             cout << " reco'd track mass = " << t.M() << endl;
0473
0474           if(t.M() > 7.0 && t.M() < 11.0)
0475             {
0476               nrecormass++;
0477               recomass->Fill(t.M());
0478               hrpt->Fill(t.Pt());
0479
0480               // Capture the mass spectrum where both tracks are the primary Upsilon decay electrons
0481               if( (trid1 == 1 || trid1 == 2) && (trid2 == 1 || trid2 == 2) )
0482             recomass_primary->Fill(t.M());
0483             }
0484         }
0485         }
0486         }
0487
0488       delete ntp_gtrack;
0489       delete ntp_track;
0490       delete ntp_cluster;
0491       delete ntp_vertex;
0492     }
0493
0494   cout << "nevents = " << nevents << endl;
0495
0496   cout << "nrecog4mass = " << nrecog4mass << endl;
0497
0498   cout << "nrecormass = " << nrecormass << endl;
0499
0500   cout << "nhittpcin per event = " << nhittpcin_cum/nhittpcin_wt << endl;
0501
0502
0503   //======================================================
0504   // End of loop over events and generation of mass histos
0505   //=====================================================
0506
0507   TCanvas *cq = new TCanvas("cq","cq",5,5,600,600 );
0508   cq->Divide(1,2);
0509   cq->cd(1);
0510   hrquality->Draw();
0511   cq->cd(2);
0512   gPad->SetLogy(1);
0513   hrdca2d->Draw();
0514
0515   // Mass histos
0516
0517   TCanvas *cmass = new TCanvas("cmass","cmass",10,10,800,600);
0518
0519   recomass_primary->SetLineColor(kRed);
0520   recomass->SetLineColor(kBlack);
0521   recomass->DrawCopy();
0522   recomass_primary->DrawCopy("same");
0523
0524   TCanvas *cm_comp = new TCanvas("cm_comp","cm_comp",10,10,800,800);
0525   cm_comp->Divide(2,1);
0526   cm_comp->cd(1);
0527   recomass->Draw();
0528   recomass_primary->Draw("same");
0529
0530   // we want from 7 to 11 GeV/c^2 - the whole range
0531   double yreco = recomass->Integral();
0532  double yreco_primary = recomass_primary->Integral();
0533  cout << "Reconstructed mass spectrum has " << yreco_primary << " entries from primary tracks and " << yreco << " entries total " << endl;
0534
0535   cm_comp->cd(2);
0536
0537   g4mass_primary->SetLineColor(kRed);
0538   g4mass->Draw();
0539   g4mass_primary->Draw("same");
0540
0541   double yg4_primary = g4mass_primary->Integral();
0542   double yg4 = g4mass->Integral();
0543   cout << "G4 mass spectrum has " << yg4_primary << " entries from primary tracks and  " << yg4 << " entries total" << endl;
0544
0545   cout << "Reconstruction efficiency is " << yreco_primary/yg4_primary << endl;
0546
0547   // Output histos for reconstructed Upsilons
0548
0549   char fname[500];
0550
0551   sprintf(fname,"root_files/ups%is_qual%.2f_dca2d%.2f.root",ups_state,quality_cut,dca_cut);
0552
0553   cout << "Create output file " << fname << endl;
0554
0555   TFile *fout1 = new TFile(fname,"recreate");
0556   recomass->Write();
0557   recomass_primary->Write();
0558   g4mass->Write();
0559   g4mass_primary->Write();
0560   hrpt->Write();
0561   hrquality->Write();
0562   hrdca2d->Write();
0563   fout1->Close();
0564
0565   cout << "Finished write to file " << fname << endl;
0566
0567 }