Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
eicsmear_dvmp_tree.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file eicsmear_dvmp_tree.C
1 #include <vector>
2 int eicsmear_dvmp_tree(TString filename_mc = "/sphenix/user/gregtom3/data/Fall2018/JPsi_reco_studies/reconstructedQ2/18x275_DVMP_1M_ascii_converted.root",
3  TString filename_mc_smeared = "/sphenix/user/gregtom3/data/Fall2018/JPsi_reco_studies/reconstructedQ2/18x275_DVMP_1M_ascii_converted.smear.root")
4 {
5 
6  /* PRELIMINARY ROOT STUFF */
7 
8  /* Loading libraries and setting sphenix style */
9  gSystem->Load("libeicsmear");
10  gROOT->LoadMacro("./sPHENIXStyle/sPhenixStyle.C");
11  SetsPhenixStyle();
12 
13 
14 
15  /* INPUT */
16 
17 
18 
19  /* Input files */
20  TFile *file_mc = new TFile(filename_mc, "OPEN");
21  TFile *file_mc_smeared = new TFile(filename_mc_smeared, "OPEN");
22 
23  /* Trees */
24  TTree *tree = (TTree*)file_mc->Get("EICTree");
25  Double32_t Q2fromEICTree;
26  tree->SetBranchAddress("trueQ2",&Q2fromEICTree);
27  TTree *tree_smeared= (TTree*)file_mc_smeared->Get("Smeared");
28 
29  /* Output File */
30  TFile *myfile = new TFile("/sphenix/user/gregtom3/data/Fall2018/JPsi_reco_studies/reconstructedQ2/analysisTree.root","RECREATE");
31  /* Add friend */
32  tree->AddFriend(tree_smeared);
33 
34 
35  /* Create data tree because this is getting hectic */
36  TTree *theTree = new TTree("DVMP_Tree","");
37  Float_t de_phi_truth, de_eta_truth, de_pt_truth, de_p_truth;
38  Float_t dp_phi_truth, dp_eta_truth, dp_pt_truth, dp_p_truth;
39  Float_t sp_phi_truth, sp_eta_truth, sp_pt_truth, sp_p_truth;
40  Float_t se_phi_truth, se_eta_truth, se_pt_truth, se_p_truth;
41  Float_t jpsi_phi_truth, jpsi_eta_truth, jpsi_pt_truth, jpsi_p_truth;
42  Float_t de_phi_reco, de_eta_reco, de_pt_reco, de_p_reco;
43  Float_t dp_phi_reco, dp_eta_reco, dp_pt_reco, dp_p_reco;
44  Float_t sp_phi_reco, sp_eta_reco, sp_pt_reco, sp_p_reco;
45  Float_t se_phi_reco, se_eta_reco, se_pt_reco, se_p_reco;
46  Float_t jpsi_phi_reco, jpsi_eta_reco, jpsi_pt_reco, jpsi_p_reco;
47  Float_t invariant_reco_decay;
48  Float_t invariant_reco_scatter;
49  Float_t jpsi_px_truth, jpsi_py_truth, jpsi_pz_truth;
50  Float_t jpsi_px_reco, jpsi_py_reco, jpsi_pz_reco;
51  Double_t Q2;
52 
53  Float_t t;
54  theTree->Branch("Q2",&Q2,"Event QSquared");
55  theTree->Branch("t",&t,"Event t");
56  theTree->Branch("de_phi_truth",&de_phi_truth,"Decay Electron Truth Phi");
57  theTree->Branch("de_eta_truth",&de_eta_truth,"Decay Electron Truth Eta");
58  theTree->Branch("de_pt_truth",&de_pt_truth,"Decay Electron Truth Pt");
59  theTree->Branch("de_p_truth",&de_p_truth,"Decay Electron Truth P");
60  theTree->Branch("dp_phi_truth",&dp_phi_truth,"Decay Positron Truth Phi");
61  theTree->Branch("dp_eta_truth",&dp_eta_truth,"Decay Positron Truth Eta");
62  theTree->Branch("dp_pt_truth",&dp_pt_truth,"Decay Positron Truth Pt");
63  theTree->Branch("dp_p_truth",&dp_p_truth,"Decay Positron Truth P");
64  theTree->Branch("sp_phi_truth",&sp_phi_truth,"Scattered Proton Truth Phi");
65  theTree->Branch("sp_eta_truth",&sp_eta_truth,"Scattered Proton Truth Eta");
66  theTree->Branch("sp_pt_truth",&sp_pt_truth,"Scattered Proton Truth Pt");
67  theTree->Branch("sp_p_truth",&sp_p_truth,"Scattered Proton Truth P");
68  theTree->Branch("se_phi_truth",&se_phi_truth,"Scattered Electron Truth Phi");
69  theTree->Branch("se_eta_truth",&se_eta_truth,"Scattered Electron Truth Eta");
70  theTree->Branch("se_pt_truth",&se_pt_truth,"Scattered Electron Truth Pt");
71  theTree->Branch("se_p_truth",&se_p_truth,"Scattered Electron Truth P");
72  theTree->Branch("jpsi_phi_truth",&jpsi_phi_truth,"Jpsi Truth Phi");
73  theTree->Branch("jpsi_eta_truth",&jpsi_eta_truth,"Jpsi Truth Eta");
74  theTree->Branch("jpsi_pt_truth",&jpsi_pt_truth,"Jpsi Truth Pt");
75  theTree->Branch("jpsi_p_truth",&jpsi_p_truth,"Jpsi Truth P");
76  theTree->Branch("de_phi_reco",&de_phi_reco,"Decay Electron Reco Phi");
77  theTree->Branch("de_eta_reco",&de_eta_reco,"Decay Electron Reco Eta");
78  theTree->Branch("de_pt_reco",&de_pt_reco,"Decay Electron Reco Pt");
79  theTree->Branch("de_p_reco",&de_p_reco,"Decay Electron Reco P");
80  theTree->Branch("dp_phi_reco",&dp_phi_reco,"Decay Positron Reco Phi");
81  theTree->Branch("dp_eta_reco",&dp_eta_reco,"Decay Positron Reco Eta");
82  theTree->Branch("dp_pt_reco",&dp_pt_reco,"Decay Positron Reco Pt");
83  theTree->Branch("dp_p_reco",&dp_p_reco,"Decay Positron Reco P");
84  theTree->Branch("sp_phi_reco",&sp_phi_reco,"Scattered Proton Reco Phi");
85  theTree->Branch("sp_eta_reco",&sp_eta_reco,"Scattered Proton Reco Eta");
86  theTree->Branch("sp_pt_reco",&sp_pt_reco,"Scattered Proton Reco Pt");
87  theTree->Branch("sp_p_reco",&sp_p_reco,"Scattered Proton Reco P");
88  theTree->Branch("se_phi_reco",&se_phi_reco,"Scattered Electron Reco Phi");
89  theTree->Branch("se_eta_reco",&se_eta_reco,"Scattered Electron Reco Eta");
90  theTree->Branch("se_pt_reco",&se_pt_reco,"Scattered Electron Reco Pt");
91  theTree->Branch("se_p_reco",&se_p_reco,"Scattered Electron Reco P");
92  theTree->Branch("jpsi_phi_reco",&jpsi_phi_reco,"Jpsi Reco Phi");
93  theTree->Branch("jpsi_eta_reco",&jpsi_eta_reco,"Jpsi Reco Eta");
94  theTree->Branch("jpsi_pt_reco",&jpsi_pt_reco,"Jpsi Reco Pt");
95  theTree->Branch("jpsi_p_reco",&jpsi_p_reco,"Jpsi Reco P");
96  theTree->Branch("jpsi_px_truth",&jpsi_px_truth,"Jpsi Truth Px");
97  theTree->Branch("jpsi_py_truth",&jpsi_py_truth,"Jpsi Truth Py");
98  theTree->Branch("jpsi_pz_truth",&jpsi_pz_truth,"Jpsi Truth Pz");
99  theTree->Branch("jpsi_px_reco",&jpsi_px_reco,"Jpsi Reco Px");
100  theTree->Branch("jpsi_py_reco",&jpsi_py_reco,"Jpsi Reco Py");
101  theTree->Branch("jpsi_pz_reco",&jpsi_pz_reco,"Jpsi Reco Pz");
102  theTree->Branch("invariant_reco_decay",&invariant_reco_decay,"Invariant Mass Decay Particles");
103  theTree->Branch("invariant_reco_scatter",&invariant_reco_scatter,"Invariant Mass Scattered Electron");
104 
105  /* Momentum cut for plotting */
106  TCut pcut("Smeared.particles.p<=50&&Smeared.particles.p>1");
107 
108  /* Eta cut for plotting */
109  TCut ecut("-log(tan(Smeared.particles.theta/2))<=4&&-log(tan(Smeared.particles.theta/2))>=-4");
110 
111  /* Scattered lepton cut */
112  TCut cut_scattered_lepton("particles.id==11&&particles.KS==2&&Smeared.particles.p<=50&&Smeared.particles.p>1");
113 
114  /* Decay electron cut */
115  TCut cut_decay_electron("particles.id==11&&particles.KS==1&&Smeared.particles.p<=50&&Smeared.particles.p>1");
116 
117  /* Decay positron cut */
118  TCut cut_decay_positron("particles.id==-11&&particles.KS==1&&Smeared.particles.p<=50&&Smeared.particles.p>1");
119 
120 
121 
122  /* PLOTTING */
123 
124  erhic::EventMilou *event = NULL;
125  Smear::Event *eventS = NULL;
126 
127  tree->SetBranchAddress("event",&event);
128  tree_smeared->SetBranchAddress("eventS", &eventS);
129 
130  unsigned max_event = tree_smeared->GetEntries();
131  std::vector<float> particle_eta;
132  particle_eta.push_back(0);
133  particle_eta.push_back(0);
134  particle_eta.push_back(0);
135  particle_eta.push_back(0);
136  particle_eta.push_back(0);
137  std::vector<float> particle_theta;
138  particle_theta.push_back(0);
139  particle_theta.push_back(0);
140  particle_theta.push_back(0);
141  particle_theta.push_back(0);
142  particle_theta.push_back(0);
143  std::vector<float> particle_phi;
144  particle_phi.push_back(0);
145  particle_phi.push_back(0);
146  particle_phi.push_back(0);
147  particle_phi.push_back(0);
148  particle_phi.push_back(0);
149  std::vector<float> particle_pt;
150  particle_pt.push_back(0);
151  particle_pt.push_back(0);
152  particle_pt.push_back(0);
153  particle_pt.push_back(0);
154  particle_pt.push_back(0);
155  std::vector<float> particle_p;
156  particle_p.push_back(0);
157  particle_p.push_back(0);
158  particle_p.push_back(0);
159  particle_p.push_back(0);
160  particle_p.push_back(0);
161 
162  std::vector<float> true_particle_eta;
163  true_particle_eta.push_back(0);
164  true_particle_eta.push_back(0);
165  true_particle_eta.push_back(0);
166  true_particle_eta.push_back(0);
167  true_particle_eta.push_back(0);
168  std::vector<float> true_particle_theta;
169  true_particle_theta.push_back(0);
170  true_particle_theta.push_back(0);
171  true_particle_theta.push_back(0);
172  true_particle_theta.push_back(0);
173  true_particle_theta.push_back(0);
174  std::vector<float> true_particle_phi;
175  true_particle_phi.push_back(0);
176  true_particle_phi.push_back(0);
177  true_particle_phi.push_back(0);
178  true_particle_phi.push_back(0);
179  true_particle_phi.push_back(0);
180  std::vector<float> true_particle_p;
181  true_particle_p.push_back(0);
182  true_particle_p.push_back(0);
183  true_particle_p.push_back(0);
184  true_particle_p.push_back(0);
185  true_particle_p.push_back(0);
186  std::vector<float> true_particle_pt;
187  true_particle_pt.push_back(0);
188  true_particle_pt.push_back(0);
189  true_particle_pt.push_back(0);
190  true_particle_pt.push_back(0);
191  true_particle_pt.push_back(0);
192  for ( unsigned ievent = 0; ievent < max_event/1000; ievent++ )
193  {
194  if ( ievent%1000 == 0 )
195  cout << "Processing event " << ievent << endl;
196 
197  /* load event */
198  tree->GetEntry(ievent);
199  tree_smeared->GetEntry(ievent);
200  cout << Q2fromEICTree << endl;
201 
202  unsigned ntracks = eventS->GetNTracks();
203 
204 
205 
206 
207  for ( unsigned itrack = 0; itrack < ntracks; itrack++ )
208  {
209 
210  Smear::ParticleMCS * sparticle = eventS->GetTrack( itrack );
211  erhic::ParticleMC * tparticle = event->GetTrack(itrack);
212  int type_of_particle = -1;
213 
214  // Decay Electron
215  if(tparticle->Id().Code()==11&&tparticle->GetStatus()==1)
216  {
217  type_of_particle = 0;
218  }
219 
220  // Scattered Electron
221  if(tparticle->Id().Code()==11&&tparticle->GetStatus()==2)
222  {
223  type_of_particle = 1;
224  }
225 
226  // Decay Positron
227  if(tparticle->Id().Code()==-11&&tparticle->GetStatus()==1)
228  {
229  type_of_particle = 2;
230  }
231 
232  // Scattered Proton
233  if(tparticle->Id().Code()==2212&&tparticle->GetStatus()==1)
234  {
235  type_of_particle = 3;
236  }
237  // JPsi
238  if(tparticle->Id().Code()==443)
239  {
240  type_of_particle = 4;
241  }
242 
243  // Test that the smear particle was found
244  if(sparticle==NULL&&type_of_particle!=4)
245  {
246  particle_eta.at(type_of_particle)=NULL;
247  particle_theta.at(type_of_particle)=NULL;
248  particle_phi.at(type_of_particle)=NULL;
249  particle_pt.at(type_of_particle)=NULL;
250  particle_p.at(type_of_particle)=NULL;
251 
252  true_particle_eta.at(type_of_particle)=tparticle->GetEta();
253  true_particle_theta.at(type_of_particle)=tparticle->GetTheta();
254  true_particle_phi.at(type_of_particle)=tparticle->GetPhi();
255  true_particle_p.at(type_of_particle)=tparticle->GetP();
256  true_particle_pt.at(type_of_particle)=tparticle->GetPt();
257 
258  //true_particle_eta.at(type_of_particle)=NULL;
259  //true_particle_phi.at(type_of_particle)=NULL;
260  //true_particle_p.at(type_of_particle)=NULL;
261  //true_particle_pt.at(type_of_particle)=NULL;
262  }
263  else if(type_of_particle!=4)
264  {
265 
266  particle_eta.at(type_of_particle)=sparticle->GetEta();
267  particle_theta.at(type_of_particle)=2*atan(exp(-sparticle->GetEta()));
268  particle_phi.at(type_of_particle)=sparticle->GetPhi();
269  particle_pt.at(type_of_particle)=sparticle->GetPt();
270  particle_p.at(type_of_particle)=sparticle->GetP();
271 
272  true_particle_eta.at(type_of_particle)=tparticle->GetEta();
273  true_particle_theta.at(type_of_particle)=tparticle->GetTheta();
274  true_particle_phi.at(type_of_particle)=tparticle->GetPhi();
275  true_particle_p.at(type_of_particle)=tparticle->GetP();
276  true_particle_pt.at(type_of_particle)=tparticle->GetPt();
277  }
278  else if(type_of_particle==4)
279  {
280  true_particle_eta.at(type_of_particle)=tparticle->GetEta();
281  true_particle_theta.at(type_of_particle)=tparticle->GetTheta();
282  true_particle_phi.at(type_of_particle)=tparticle->GetPhi();
283  true_particle_p.at(type_of_particle)=tparticle->GetP();
284  true_particle_pt.at(type_of_particle)=tparticle->GetPt();
285  }
286 
287  } // end loop over particles
288 
289 
290  bool can_we_reconstruct_3_final_state_leptons = true;
291  if((particle_p.at(0)<1||particle_p.at(1)<1)||particle_p.at(2)<1)
292  {
293  can_we_reconstruct_3_final_state_leptons = false;
294  }
295 
296  bool can_we_reconstruct_scattered_lepton = true;
297  if(particle_p.at(1)<1)
298  {
299  can_we_reconstruct_scattered_lepton = false;
300  }
301 
302  bool can_we_reconstruct_decay_particles = true;
303  if(particle_p.at(0)<1||particle_p.at(2)<1)
304  {
305  can_we_reconstruct_decay_particles = false;
306  }
307 
308  // From the eta, phi and p arrays, calculate the invariant mass
309  // First, try the decay electron and decay positron
310  if(can_we_reconstruct_3_final_state_leptons)
311  {
312  float m = sqrt(2*particle_pt.at(0)*particle_pt.at(2)*(cosh(particle_eta.at(0)-particle_eta.at(2))-cos(particle_phi.at(0)-particle_phi.at(2))));
313  invariant_reco_decay=m;
314  }
315 
316  // Now try the scattered electron and decay positron
317 
318  if(can_we_reconstruct_3_final_state_leptons)
319  {
320  float m = sqrt(2*particle_pt.at(1)*particle_pt.at(2)*(cosh(particle_eta.at(1)-particle_eta.at(2))-cos(particle_phi.at(1)-particle_phi.at(2))));
321  invariant_reco_scatter=m;
322  }
323  // cout << particle_theta.at(0) << endl;
324  // Next, fill in reco J/Psi stuff
325  if(can_we_reconstruct_3_final_state_leptons)
326  {
327  float sum_px = particle_pt.at(0)*cos(particle_phi.at(0))+particle_pt.at(2)*cos(particle_phi.at(2));
328  float sum_py = particle_pt.at(0)*sin(particle_phi.at(0))+particle_pt.at(2)*sin(particle_phi.at(2));
329  float sum_pz = particle_p.at(0)*cos(particle_theta.at(0))+particle_p.at(2)*cos(particle_theta.at(2));
330  float sum_p = sqrt(sum_px*sum_px+sum_py*sum_py+sum_pz*sum_pz);
331  jpsi_px_reco=sum_px;
332  jpsi_py_reco=sum_py;
333  jpsi_pz_reco=sum_pz;
334 
335 
336  float sum_px_true = true_particle_pt.at(0)*cos(true_particle_phi.at(0))+true_particle_pt.at(2)*cos(true_particle_phi.at(2));
337  float sum_py_true = true_particle_pt.at(0)*sin(true_particle_phi.at(0))+true_particle_pt.at(2)*sin(true_particle_phi.at(2));
338  float sum_pz_true = true_particle_p.at(0)*cos(true_particle_theta.at(0))+true_particle_p.at(2)*cos(true_particle_theta.at(2));
339  jpsi_px_truth=sum_px_true;
340  jpsi_py_truth=sum_py_true;
341  jpsi_pz_truth=sum_pz_true;
342 
343 
344 
345  // if(sum_p>4.8&&sum_p<5)
346  //cout << ievent << " - " << true_particle_p.at(4) << " : " << sum_p << " : " << sum_pz << " : " << sum_px << " : " << sum_py << " : " << endl;
347  //float jpsi_eta = -log(tan( (acos( (particle_p.at(0)*cos(2*atan(exp(-particle_eta.at(0))))+particle_p.at(2)*cos(2*atan(exp(-particle_eta.at(2))))/(sum_p) )))/2));
348  //float jpsi_eta = -log(tan((acos(particle_p.at(0)*cos(particle_theta.at(0))+particle_p.at(2)*cos(particle_theta.at(2))))
349  // float jpsi_eta = -log(tan(atan( (particle_pt.at(0)+particle_pt.at(2))/(particle_p.at(0)*cos(particle_theta.at(0))+particle_p.at(2)*cos(particle_theta.at(2))))/2));
350  // float jpsi_phi = acos( (particle_pt.at(0)*cos((particle_phi.at(0)))+particle_pt.at(2)*cos((particle_phi.at(2))))/(sqrt(sum_px*sum_px+sum_py*sum_py)) );
351  float jpsi_eta = -log(tan( (acos(jpsi_pz_reco/sum_p))/2 ) );
352  float jpsi_phi = atan(jpsi_py_reco/jpsi_px_reco);
353  if(jpsi_py_reco>0&&jpsi_px_reco<0)
354  jpsi_phi=3.14159265+jpsi_phi;
355  else if(jpsi_py_reco<0&&jpsi_px_reco<0)
356  jpsi_phi=3.14159265+jpsi_phi;
357  else if(jpsi_py_reco<0&&jpsi_px_reco>0)
358  jpsi_phi=2*3.14159265+jpsi_phi;
359 
360  //cout << true_particle_phi.at(4) << " : " << jpsi_px_truth << " : " << jpsi_py_truth << " : " << atan(jpsi_py_truth/jpsi_px_truth)<< endl;
361  //cout << true_particle_phi.at(4) << " : " << jpsi_phi<< endl;
362  //cout << jpsi_phi << " : " << true_particle_phi.at(4) << endl;
363  particle_p.at(4)=sum_p;
364  particle_pt.at(4)=sum_p*sin(jpsi_phi);
365  particle_eta.at(4)=jpsi_eta;
366  particle_phi.at(4)=jpsi_phi;
367 
368 
369  // Event Kinematics
370 
371 
372  t = 2*275*true_particle_p.at(3)*(1-cos(2*atan(exp(-true_particle_eta.at(3)))));
373  }
374 
375  de_phi_truth=true_particle_phi.at(0);
376  de_eta_truth=true_particle_eta.at(0);
377  de_pt_truth=true_particle_pt.at(0);
378  de_p_truth=true_particle_p.at(0);
379  dp_phi_truth=true_particle_phi.at(2);
380  dp_eta_truth=true_particle_eta.at(2);
381  dp_pt_truth=true_particle_pt.at(2);
382  dp_p_truth=true_particle_p.at(2);
383  sp_phi_truth=true_particle_phi.at(3);
384  sp_eta_truth=true_particle_eta.at(3);
385  sp_pt_truth=true_particle_pt.at(3);
386  sp_p_truth=true_particle_p.at(3);
387  se_phi_truth=true_particle_phi.at(1);
388  se_eta_truth=true_particle_eta.at(1);
389  se_pt_truth=true_particle_pt.at(1);
390  se_p_truth=true_particle_p.at(1);
391  jpsi_phi_truth=true_particle_phi.at(4);
392  jpsi_eta_truth=true_particle_eta.at(4);
393  jpsi_pt_truth=true_particle_pt.at(4);
394  jpsi_p_truth=true_particle_p.at(4);
395  de_phi_reco=particle_phi.at(0);
396  de_eta_reco=particle_eta.at(0);
397  de_pt_reco=particle_pt.at(0);
398  de_p_reco=particle_p.at(0);
399  dp_phi_reco=particle_phi.at(2);
400  dp_eta_reco=particle_eta.at(2);
401  dp_pt_reco=particle_pt.at(2);
402  dp_p_reco=particle_p.at(2);
403  sp_phi_reco=particle_phi.at(3);
404  sp_eta_reco=particle_eta.at(3);
405  sp_pt_reco=particle_pt.at(3);
406  sp_p_reco=particle_p.at(3);
407  se_phi_reco=particle_phi.at(1);
408  se_eta_reco=particle_eta.at(1);
409  se_pt_reco=particle_pt.at(1);
410  se_p_reco=particle_p.at(1);
411  jpsi_phi_reco=particle_phi.at(4);
412  jpsi_eta_reco=particle_eta.at(4);
413  jpsi_pt_reco=particle_pt.at(4);
414  jpsi_p_reco=particle_p.at(4);
415  // if(ievent==770)
416  // cout << jpsi_p_truth << " : " << jpsi_p_reco << endl;
417  theTree->Fill();
418  } // end loop over events
419 
421 
422 
423 
424  myfile->cd();
425  theTree->Write("Tree");
426  myfile->Write();
427  myfile->Close();
428  return 0;
429 }
430 
431 float eta_to_theta(float eta)
432 {
433  return 2*atan(exp(-1*eta));
434 }
435 
436 float theta_to_eta(float theta)
437 {
438  return -log(tan(theta/2));
439 }