Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
quarkonia_reconstruction_embedded.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file quarkonia_reconstruction_embedded.C
1 #include <TChain.h>
2 #include <TCanvas.h>
3 #include <TH1.h>
4 #include <TH1F.h>
5 #include <TH2D.h>
6 #include <TF1.h>
7 #include <TFile.h>
8 #include <TGraph.h>
9 #include <TStyle.h>
10 #include <TROOT.h>
11 #include <TLegend.h>
12 #include <TLine.h>
13 #include <TLatex.h>
14 #include <TRandom1.h>
15 #include <TPolyLine.h>
16 #include <iostream>
17 #include <fstream>
18 #include <TMath.h>
19 //#include <iomanip.h>
20 #include <TLorentzVector.h>
21 
22 //#define TEST
23 
24 void quarkonia_reconstruction_embedded()
25 {
26  gROOT->SetStyle("Plain");
27  gStyle->SetOptStat(1);
28  gStyle->SetOptTitle(1);
29 
30  int verbose = 0;
31  int embed_flag = 3; // embedding flag used during Upsilon generation
32 
33  // track cuts
34  double quality_cut = 3.0;
35  double dca_cut = 0.1;
36 
37  char lepton[100];
38  sprintf(lepton,"electron");
39 
40  double decaymass=0.000511;
41  cout << "Assuming decay particle mass is " << decaymass << endl;
42 
43  // Define some histograms
44 
45  int nbpt = 20;
46  double ptmax = 10.0;
47 
48  TH1F* hrquality = new TH1F("hrquality", "Reconstructed track Quality", 1000, 0, 10);
49  TH1F* hrdca2d = new TH1F("hrdca2d", "Reconstructed track dca2d", 1000, 0, 0.05);
50  TH1F* hrpt = new TH1F("hrpt"," pT", nbpt, 0.0, ptmax);
51  TH1F* hgpt = new TH1F("hgpt","g4 pT", nbpt, 0.0, ptmax);
52 
53  TH1D *g4mass = new TH1D("g4mass","G4 input invariant mass",200,7.0,11.0);
54  g4mass->GetXaxis()->SetTitle("invariant mass (GeV/c^{2})");
55  TH1D *g4mass_primary = new TH1D("g4mass_primary","G4 input invariant mass",200,7.0,11.0);
56  g4mass_primary->GetXaxis()->SetTitle("invariant mass (GeV/c^{2})");
57 
58  TH1D *recomass = new TH1D("recomass","Reconstructed invariant mass",200,7.0,11.0);
59  recomass->GetXaxis()->SetTitle("invariant mass (GeV/c^{2})");
60  TH1D *recomass_primary = new TH1D("recomass_primary","Reconstructed invariant mass",200,7.0,11.0);
61  recomass_primary->GetXaxis()->SetTitle("invariant mass (GeV/c^{2})");
62 
63  int ups_state = 1; // used in naming of output files
64  //int ups_state = 2; // used in naming of output files
65  //int ups_state = 3; // used in naming of output files
66 
67  // Number of upsilons to process - quit after this number is reached
68  int nups_requested = 0;
69  double pair_acc = 0.60 * 0.9 * 0.98; // fraction into 1 unit of rapidity * pair eID fraction * trigger fraction
70  // These values are for 1 year of pp running
71  if(ups_state == 1) nups_requested = pair_acc * (8769 * 1.126) / (0.38 * 0.9 * 0.98); // pair_acc * yield for 197 /pb (= 15590)
72  if(ups_state == 2) nups_requested = pair_acc * (2205 * 1.126) / (0.38 * 0.9 * 0.98); // pair_acc * yield for 197 /pb (= 3920)
73  if(ups_state == 3) nups_requested = pair_acc * (1156 * 1.126) / (0.38 * 0.9 * 0.98); // pair_acc * yield for 197 /pb (=2055)
74 
75  nups_requested = 100000; // take them all
76 
77  cout << "ups_state = " << ups_state << " Upsilons requested = " << nups_requested << endl;
78 
79  int nrecog4mass = 0;
80  int nrecormass = 0;
81 
82  // Open the g4 evaluator output file
83  cout << "Reading electron ntuples " << endl;
84  int nevents = 0;
85  double nhittpcin_cum = 0;
86  double nhittpcin_wt = 0;
87 
88  // The condor job output files - we open them one at a time and process them
89  for(int i=0;i<2000;i++)
90  {
91  if(nrecormass > nups_requested)
92  {
93  cout << "Reached requested number of reco Upsilons, quit out of file loop" << endl;
94  break;
95  }
96 
97  char name[500];
98 
99  //sprintf(name,"/sphenix/user/frawley/macros_newTPC_june6/macros/macros/g4simulations/intt_6layers_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",i);
100  //sprintf(name,"/sphenix/user/frawley/macros_newTPC_june6/macros/macros/g4simulations/intt_4layers_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",i);
101  //sprintf(name,"/sphenix/user/frawley/macros_newTPC_june6/macros/macros/g4simulations/intt_8layers_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",i);
102  sprintf(name,"/sphenix/user/frawley/macros_newTPC_june6/macros/macros/g4simulations/intt_0layers_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",i);
103 
104 
105  /*
106  // ups states, 100 poins, 80 ns
107  if(ups_state == 1)
108  {
109  //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);
110  //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);
111  }
112  else if(ups_state == 2)
113  {
114  //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);
115  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);
116  }
117  else if(ups_state == 3)
118  {
119  //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);
120  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);
121  }
122  else
123  {
124  cout << "Oops!" << endl;
125  exit(1);
126  }
127  */
128 
129 
130  cout << "Adding file " << name << endl;
131 
132  TChain* ntp_track = new TChain("ntp_track","reco tracks");
133  TChain* ntp_gtrack = new TChain("ntp_gtrack","g4 tracks");
134  TChain* ntp_vertex = new TChain("ntp_vertex","events");
135  TChain *ntp_cluster = new TChain("ntp_cluster","clusters");
136  TChain *ntp_g4hit = new TChain("ntp_g4hit","G4 hits");
137 
138  ntp_vertex->Add(name);
139  ntp_track->Add(name);
140  ntp_gtrack->Add(name);
141 
142 
143  cout << "The ntuples contain " << ntp_vertex->GetEntries() << " events " << endl;
144 
145 
146  // Ntuple access variables
147  // This include file contains the definitions of the ntuple variables
148 #include "ntuple_variables.C"
149 
150  //=======================
151  // Loop over events
152  //=======================
153 
154  ntracks = ntp_track->GetEntries();
155  ngtracks = ntp_gtrack->GetEntries();
156 
157  int nr = 0;
158  int ng = 0;
159  //int nev = 1;
160  int nev = ntp_vertex->GetEntries();
161 
162  for(int iev=0;iev<nev;iev++)
163  {
164  // drop out when the requested number of reco'd Upsilons has been reached
165  if(nrecormass > nups_requested)
166  {
167  cout << "Reached requested number of reco Upsilons, quit" << endl;
168 
169  break;
170  }
171 
172  nevents++;
173 
174  int recoget = ntp_vertex->GetEntry(iev);
175 
176  nhittpcin_cum += nhittpcin;
177  nhittpcin_wt += 1.0;
178 
179  if(verbose)
180  cout
181  << "iev " << iev
182  << " event " << event
183  << " ntracks (reco) " << ntracks
184  << " ngtracks (g4) " << ngtracks
185  //<< " gvz " << egvz
186  // << " vz " << evz
187  << endl;
188 
189 
190  //============================
191  // process G4 tracks
192  // for this event
193  //============================
194 
195  int ng4trevt_elec = -1;
196  int ng4trevt_pos = -1;
197  int g4trnum_elec[1000];
198  int g4trnum_pos[1000];
199 
200 
201 
202  for(int ig=ng;ig<ng+ngtracks;ig++)
203  {
204  int recoget1 = ntp_gtrack->GetEntry(ig);
205  if(!recoget1)
206  {
207  if(verbose > 0) cout << "Did not get entry for ig = " << ig << endl;
208  break;
209  }
210 
211 
212  // This bookkeeping is needed because the evaluator records for each event the total track count in ntp_vertex, even
213  // when it writes out only embedded tracks
214  if(tevent != iev)
215  {
216  if(verbose > 0) cout << " reached new event for ig = " << ig << " tevent = " << tevent << endl;
217  ng = ig;
218  break;
219  }
220  if(ig == ng+ngtracks - 1)
221  {
222  if(verbose > 0) cout << " last time through loop for ig = " << ig << " revent = " << tevent << endl;
223  ng = ig+1;
224  }
225 
226 
227  if(tembed != embed_flag)
228  //if( tgtrackid != 1 && tgtrackid != 2)
229  continue;
230 
231 
232  // we want only electrons or positrons
233  if(tflavor != 11 && tflavor != -11)
234  continue;
235 
236  if(tflavor == 11)
237  {
238  // electron
239  ng4trevt_elec++;
240  if(ng4trevt_elec > 999)
241  continue;
242 
243  if(verbose)
244  cout << " Found electron:" << endl
245  << " ig " << ig
246  << " ng4trevt_elec " << ng4trevt_elec
247  << " gtrackID " << tgtrackid
248  << " gflavor " << tflavor
249  << " tpx " << tpx
250  << " tpy " << tpy
251  << " tpz " << tpz
252  << endl;
253 
254  g4trnum_elec[ng4trevt_elec] = ig;
255  }
256  else
257  {
258  // positron
259  ng4trevt_pos++;
260  if(ng4trevt_pos > 999)
261  continue;
262 
263  if(verbose)
264  cout << " Found positron:" << endl
265  << " ig " << ig
266  << " ng4trevt_pos " << ng4trevt_pos
267  << " gtrackID " << tgtrackid
268  << " gflavor " << tflavor
269  << " tpx " << tpx
270  << " tpy " << tpy
271  << " tpz " << tpz
272  << endl;
273 
274  g4trnum_pos[ng4trevt_pos] = ig;
275  }
276  }
277  ng4trevt_elec++;
278  ng4trevt_pos++;
279 
280  if(verbose)
281  cout << "For this event found " << ng4trevt_elec << " g4 electrons and " << ng4trevt_pos << " g4 positrons" << endl;
282 
283  // make all pairs of g4 (truth) electrons and positrons
284  for(int ielec=0;ielec<ng4trevt_elec;ielec++)
285  {
286  int recoelec = ntp_gtrack->GetEntry(g4trnum_elec[ielec]);
287 
288 
289  if(tembed != embed_flag)
290  continue;
291 
292  double elec_pT = sqrt(tpx*tpx+tpy*tpy);
293  double elec_eta = asinh(tpz/sqrt(tpx*tpx+tpy*tpy));
294 
295  int gtrid1 = tgtrackid;
296 
297  TLorentzVector t1;
298  double E1 = sqrt( pow(tpx,2) + pow(tpy,2) + pow(tpz,2) + pow(decaymass,2));
299  t1.SetPxPyPzE(tpx,tpy,tpz,E1);
300 
301  // print out track details
302  if(verbose > 1)
303  cout << " Pair electron: iev " << iev << " ielec " << ielec
304  << " g4trnum_elec " << g4trnum_elec[ielec]
305  << " tgtrackid " << tgtrackid
306  << " tflavor " << tflavor
307  << " tpx " << tpx
308  << " tpy " << tpy
309  << " tpz " << tpz
310  << " elec_eta " << elec_eta
311  << " elec_gpT " << elec_pT
312  << endl;
313 
314  for(int ipos =0;ipos<ng4trevt_pos;ipos++)
315  {
316  int recopos = ntp_gtrack->GetEntry(g4trnum_pos[ipos]);
317 
318  int gtrid2 = tgtrackid;
319 
320  double pos_pT = sqrt(tpx*tpx+tpy*tpy);
321  double pos_eta = asinh(tpz/sqrt(tpx*tpx+tpy*tpy));
322 
323  // print out track details
324  if(verbose > 1)
325  cout << " Pair positron: iev " << iev << " ipos " << ipos
326  << " g4trnum_pos " << g4trnum_pos[ipos]
327  << " tgtrackid " << tgtrackid
328  << " tflavor " << tflavor
329  << " tpx " << tpx
330  << " tpy " << tpy
331  << " tpz " << tpz
332  << " pos_eta " << pos_eta
333  << " pos_gpT " << pos_pT
334  << endl;
335 
336  // Make G4 invariant mass
337 
338  TLorentzVector t2;
339  double E2 = sqrt( pow(tpx,2) + pow(tpy,2) + pow(tpz,2) + pow(decaymass,2));
340  t2.SetPxPyPzE(tpx,tpy,tpz,E2);
341 
342  TLorentzVector t = t1+t2;
343 
344  if(verbose)
345  cout << " reco'd g4 mass = " << t.M() << endl << endl;
346 
347  if(t.M() > 7.0 && t.M() < 11.0)
348  {
349  nrecog4mass++;
350  g4mass->Fill(t.M());
351  hgpt->Fill(t.Pt());
352 
353  // Capture the mass spectrum where both tracks are the primary Upsilon decay electrons
354  g4mass_primary->Fill(t.M());
355 
356  }
357  } // end of ipos loop
358  } // end of ielec loop
359 
360 
361  if(verbose)
362  {
363  cout << " # of g4 electron tracks = " << ng4trevt_elec
364  << " # of g4 positron tracks = " << ng4trevt_pos << endl;
365  }
366 
367 
368  //=============================
369  // process reconstructed tracks
370  // for this event
371  //=============================
372 
373  int nrtrevt = 0;
374  int nr_elec = -1;
375  int nr_pos = -1;
376  int rectrnum_elec[1000];
377  int rectrnum_pos[1000];
378 
379  //cout << "Number of ntp_track entries = " << recoget << endl;
380 
381  for(int ir=nr;ir<nr+ntracks;ir++)
382  {
383  int recoget = ntp_track->GetEntry(ir);
384  if(!recoget)
385  {
386  if(verbose > 0) cout << "Did not get entry for ir = " << ir << endl;
387  break;
388  }
389 
390 
391  // This bookkeeping is needed because the evaluator records for each event the total track count in ntp_vertex, even
392  // when it writes out only embedded tracks
393  if(revent != iev)
394  {
395  if(verbose > 1) cout << " reached new event for ir = " << ir << " revent = " << revent << endl;
396  nr = ir;
397  break;
398  }
399  if(ir == nr+ntracks - 1)
400  {
401  if(verbose > 0) cout << " last time through loop for ir = " << ir << " revent = " << revent << endl;
402  nr = ir+1;
403  }
404 
405  if(rgembed != embed_flag)
406  continue;
407 
408  hrquality->Fill(rquality);
409  hrdca2d->Fill(rdca2d);
410 
411  // track quality cuts
412  if(rquality > quality_cut || fabs(rdca2d) > dca_cut)
413  continue;
414 
415  // make a list of electrons and positrons
416  if(rcharge == -1)
417  {
418  nr_elec++;
419  rectrnum_elec[nr_elec] = ir;
420  }
421 
422  if(rcharge == 1)
423  {
424  nr_pos++;
425  rectrnum_pos[nr_pos] = ir;
426  }
427 
428 
429  double rpT = sqrt(rpx*rpx+rpy*rpy);
430  double reta = asinh(rpz/sqrt(rpx*rpx+rpy*rpy));
431 
432  // print out track details
433  if(verbose)
434  cout << " revent " << revent << " ir " << ir
435  << " rgtrackid " << rgtrackid
436  << " rgflavor " << rgflavor
437  << " rvz " << rvz
438  << " reta " << reta
439  << " rpT " << rpT
440  << endl;
441  }
442 
443  nr_elec++;
444  nr_pos++;
445 
446  for(int ielec = 0;ielec<nr_elec;ielec++)
447  {
448 
449  TLorentzVector t1;
450 
451  int recoget1 = ntp_track->GetEntry(rectrnum_elec[ielec]);
452 
453  int trid1 = rgtrackid;
454 
455  double E1 = sqrt( pow(rpx,2) + pow(rpy,2) + pow(rpz,2) + pow(decaymass,2));
456  t1.SetPxPyPzE(rpx,rpy,rpz,E1);
457 
458  for(int ipos = 0;ipos<nr_pos;ipos++)
459  {
460  int recoget2 = ntp_track->GetEntry(rectrnum_pos[ipos]);
461 
462  int trid2 = rgtrackid;
463 
464  TLorentzVector t2;
465  double E2 = sqrt( pow(rpx,2) + pow(rpy,2) + pow(rpz,2) + pow(decaymass,2));
466 
467  t2.SetPxPyPzE(rpx,rpy,rpz,E2);
468 
469  TLorentzVector t = t1+t2;
470 
471  if(verbose)
472  cout << " reco'd track mass = " << t.M() << endl;
473 
474  if(t.M() > 7.0 && t.M() < 11.0)
475  {
476  nrecormass++;
477  recomass->Fill(t.M());
478  hrpt->Fill(t.Pt());
479 
480  // Capture the mass spectrum where both tracks are the primary Upsilon decay electrons
481  if( (trid1 == 1 || trid1 == 2) && (trid2 == 1 || trid2 == 2) )
482  recomass_primary->Fill(t.M());
483  }
484  }
485  }
486  }
487 
488  delete ntp_gtrack;
489  delete ntp_track;
490  delete ntp_cluster;
491  delete ntp_vertex;
492  }
493 
494  cout << "nevents = " << nevents << endl;
495 
496  cout << "nrecog4mass = " << nrecog4mass << endl;
497 
498  cout << "nrecormass = " << nrecormass << endl;
499 
500  cout << "nhittpcin per event = " << nhittpcin_cum/nhittpcin_wt << endl;
501 
502 
503  //======================================================
504  // End of loop over events and generation of mass histos
505  //=====================================================
506 
507  TCanvas *cq = new TCanvas("cq","cq",5,5,600,600 );
508  cq->Divide(1,2);
509  cq->cd(1);
510  hrquality->Draw();
511  cq->cd(2);
512  gPad->SetLogy(1);
513  hrdca2d->Draw();
514 
515  // Mass histos
516 
517  TCanvas *cmass = new TCanvas("cmass","cmass",10,10,800,600);
518 
519  recomass_primary->SetLineColor(kRed);
520  recomass->SetLineColor(kBlack);
521  recomass->DrawCopy();
522  recomass_primary->DrawCopy("same");
523 
524  TCanvas *cm_comp = new TCanvas("cm_comp","cm_comp",10,10,800,800);
525  cm_comp->Divide(2,1);
526  cm_comp->cd(1);
527  recomass->Draw();
528  recomass_primary->Draw("same");
529 
530  // we want from 7 to 11 GeV/c^2 - the whole range
531  double yreco = recomass->Integral();
532  double yreco_primary = recomass_primary->Integral();
533  cout << "Reconstructed mass spectrum has " << yreco_primary << " entries from primary tracks and " << yreco << " entries total " << endl;
534 
535  cm_comp->cd(2);
536 
537  g4mass_primary->SetLineColor(kRed);
538  g4mass->Draw();
539  g4mass_primary->Draw("same");
540 
541  double yg4_primary = g4mass_primary->Integral();
542  double yg4 = g4mass->Integral();
543  cout << "G4 mass spectrum has " << yg4_primary << " entries from primary tracks and " << yg4 << " entries total" << endl;
544 
545  cout << "Reconstruction efficiency is " << yreco_primary/yg4_primary << endl;
546 
547  // Output histos for reconstructed Upsilons
548 
549  char fname[500];
550 
551  sprintf(fname,"root_files/ups%is_qual%.2f_dca2d%.2f.root",ups_state,quality_cut,dca_cut);
552 
553  cout << "Create output file " << fname << endl;
554 
555  TFile *fout1 = new TFile(fname,"recreate");
556  recomass->Write();
557  recomass_primary->Write();
558  g4mass->Write();
559  g4mass_primary->Write();
560  hrpt->Write();
561  hrquality->Write();
562  hrdca2d->Write();
563  fout1->Close();
564 
565  cout << "Finished write to file " << fname << endl;
566 
567 }