Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Forward_pi0s.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Forward_pi0s.C
1 #include "Forward_pi0s.h"
2 
3 #include <fun4all/Fun4AllServer.h>
4 #include <g4main/PHG4Hit.h>
5 #include <g4main/PHG4Particle.h>
6 #include <g4main/PHG4TruthInfoContainer.h>
7 #include <phool/PHCompositeNode.h>
8 #include <phool/getClass.h>
9 
10 #include <TLorentzVector.h>
11 #include <calobase/RawTower.h>
12 #include <calobase/RawTowerContainer.h>
13 #include <calobase/RawTowerGeom.h>
14 #include <calobase/RawTowerGeomContainer.h>
15 #include <g4jets/Jet.h>
16 #include <g4jets/JetMap.h>
17 #include <iostream>
18 
19 #include <calobase/RawCluster.h>
20 #include <calobase/RawClusterContainer.h>
21 #include <calobase/RawClusterUtility.h>
22 #include <g4detectors/PHG4ScintillatorSlatContainer.h>
23 #include <g4eval/JetEvalStack.h>
24 #include <g4eval/SvtxEvalStack.h>
25 
26 #include <g4vertex/GlobalVertex.h>
27 #include <g4vertex/GlobalVertexMap.h>
28 #include <sstream>
29 
30 #include <HepMC/GenEvent.h>
31 #include <HepMC/GenVertex.h>
32 #include <phhepmc/PHHepMCGenEvent.h>
33 using namespace std;
34 
35 #include <cassert>
36 #include <iostream>
37 
38 Forward_pi0s::Forward_pi0s(const std::string &name)
39  : SubsysReco("FORWARD_PI0S")
40 {
41  outfilename = name;
42 
43  //add other initializers here
44  //default use isocone algorithm
45  use_isocone = 1;
46 
47  //default central arm
48  _etalow = -1;
49  _etahi = 1;
50 
51  //default to only central arm
52  _useforwardarm = 0;
53 
54  mincluspt = 0.3;
55 
56  nevents = 0;
57  //default use 0.3 jet cone
58  jet_cone_size = 0.3;
59 
60 }
61 
62 int Forward_pi0s::Init(PHCompositeNode *topnode)
63 {
64  file = new TFile(outfilename.c_str(), "RECREATE");
65 
66  histo = new TH1F("histo", "histo", 100, -3, 3);
67 
68  tree = new TTree("tree", "a tree");
69  tree->Branch("nevents", &nevents, "nevents/I");
70 
71  Set_Tree_Branches();
72  initialize_things();
73 
74  return 0;
75 }
76 
77 int Forward_pi0s::process_event(PHCompositeNode *topnode)
78 {
79  nevents++;
80  cout << "at event number " << nevents << endl;
81 
82  //reset the truth variables
83  tphote2 = -999;
84  tphotpx2 = -999;
85  tphotpy2 = -999;
86  tphotpz2 = -999;
87  tphotpid2 = -999;
88  tphotparentid2 = -999;
89  tphotpt2 = -999;
90  tphotphi2 = -999;
91  tphoteta2 = -999;
92 
93  tphote1 = -999;
94  tphotpx1 = -999;
95  tphotpy1 = -999;
96  tphotpz1 = -999;
97  tphotpid1 = -999;
98  tphotparentid1 = -999;
99  tphotpt1 = -999;
100  tphotphi1 = -999;
101  tphoteta1 = -999;
102 
103  tpi0e = -999;
104  tpi0px = -999;
105  tpi0py = -999;
106  tpi0pz = -999;
107  tpi0pid = -999;
108  tpi0pt = -999;
109  tpi0phi = -999;
110  tpi0eta = -999;
111 
112  //reset cluster variables
113  fclusenergy = -999;
114  fclus_eta = -999;
115  fclus_phi = -999;
116  fclus_theta = -999;
117  fclus_pt = -999;
118  fclus_px = -999;
119  fclus_py = -999;
120  fclus_pz = -999;
121 
122 
123  clus_energy = -999;
124  clus_eta = -999;
125  clus_theta = -999;
126  clus_pt = -999;
127  clus_phi = -999;
128  clus_px = -999;
129  clus_py = -999;
130  clus_pz = -999;
131 
132 
133 
134 
135 
136 
137  //get the nodes from the NodeTree
138 
139  PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topnode, "G4TruthInfo");
140  RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topnode, "CLUSTER_CEMC");
141  RawClusterContainer *fclusters = 0;
142  if (_useforwardarm)
143  {
144  fclusters = findNode::getClass<RawClusterContainer>(topnode, "CLUSTER_FEMC");
145  }
146  PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
147 
148  GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topnode, "GlobalVertexMap");
149  if (!vertexmap)
150  {
151  cout << "Forward_pi0s::process_event - Fatal Error - GlobalVertexMap node is missing. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << endl;
152  assert(vertexmap); // force quit
153 
154  return 0;
155  }
156 
157  if (vertexmap->empty())
158  {
159  cout << "Forward_pi0s::process_event - Fatal Error - GlobalVertexMap node is empty. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << endl;
160  return 0;
161  }
162 
163  GlobalVertex *vtx = vertexmap->begin()->second;
164  if (vtx == nullptr) return 0;
165 
166 
167  if (!fclusters && _useforwardarm)
168  {
169  cout << "not forward cluster info" << endl;
170  return 0;
171  }
172  if (!truthinfo)
173  {
174  cout << "no truth track info" << endl;
175  return 0;
176  }
177 
178  if (!clusters)
179  {
180  cout << "no cluster info" << endl;
181  return 0;
182  }
183 
184 
185  for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
186  {
187  PHG4Particle *truth = iter->second;
188 
189  truthpx = truth->get_px();
190  truthpy = truth->get_py();
191  truthpz = truth->get_pz();
192  truthp = sqrt(truthpx * truthpx + truthpy * truthpy + truthpz * truthpz);
193  truthenergy = truth->get_e();
194 
195  TLorentzVector vec;
196  vec.SetPxPyPzE(truthpx, truthpy, truthpz, truthenergy);
197 
198  truthpt = sqrt(truthpx * truthpx + truthpy * truthpy);
199 
200  truthphi = vec.Phi();
201  trutheta = vec.Eta();
202  truthpid = truth->get_pid();
203 
204  //get the pi0 information
205  if (truthpid == 111)
206  {
207  tpi0e = truth->get_e();
208  tpi0px = truth->get_px();
209  tpi0py = truth->get_py();
210  tpi0pz = truth->get_pz();
211  tpi0pid = truth->get_pid();
212  tpi0pt = sqrt(tpi0px * tpi0px + tpi0py * tpi0py);
213  TLorentzVector pi0;
214  pi0.SetPxPyPzE(tpi0px, tpi0py, tpi0pz, tpi0e);
215  tpi0phi = pi0.Phi();
216  tpi0eta = pi0.Eta();
217  }
218 
219  truth_g4particles->Fill();
220  }
221 
222  if (Verbosity() > 1)
223  cout << "Truth pi0 energy is :" << tpi0e << endl;
224 
225  int firstphotonflag = 0;
226  int firstfirstphotonflag = 0;
227  int secondphotonflag = 0;
228  int secondsecondphotonflag = 0;
229  PHG4TruthInfoContainer::Range ranger = truthinfo->GetSecondaryParticleRange();
230  //This sets the global variables in the if statements to the pion first decay photon truth values and second decay photon truth values
231  for (PHG4TruthInfoContainer::ConstIterator iter = ranger.first; iter != ranger.second; ++iter)
232  {
233  PHG4Particle *truth = iter->second;
234 
235  int dumtruthpid = truth->get_pid();
236  int dumparentid = truth->get_parent_id();
237 
238 
239  //if the parent is the pi0 and its a photon and we haven't marked one yet
240  if (dumparentid == 1 && dumtruthpid == 22 && !firstphotonflag)
241  {
242  tphote1 = truth->get_e();
243  tphotpx1 = truth->get_px();
244  tphotpy1 = truth->get_py();
245  tphotpz1 = truth->get_pz();
246  tphotpid1 = truth->get_pid();
247  tphotparentid1 = truth->get_parent_id();
248  tphotpt1 = sqrt(tphotpx1 * tphotpx1 + tphotpy1 * tphotpy1);
249 
250  TLorentzVector vec;
251  vec.SetPxPyPzE(tphotpx1, tphotpy1, tphotpz1, tphote1);
252 
253  tphotphi1 = vec.Phi();
254  tphoteta1 = vec.Eta();
255  firstphotonflag = 1;
256  }
257 
258  //if the parent is the pi0 and its a photon and we have
259  //marked the first photon but not the second photon
260  if (dumparentid == 1 &&
261  dumtruthpid == 22 &&
262  firstphotonflag && firstfirstphotonflag)
263  {
264  tphote2 = truth->get_e();
265  tphotpx2 = truth->get_px();
266  tphotpy2 = truth->get_py();
267  tphotpz2 = truth->get_pz();
268  tphotpid2 = truth->get_pid();
269  tphotparentid2 = truth->get_parent_id();
270  tphotpt2 = sqrt(tphotpx2 * tphotpx2 + tphotpy2 * tphotpy2);
271 
272  TLorentzVector vec;
273  vec.SetPxPyPzE(tphotpx2, tphotpy2, tphotpz2, tphote2);
274 
275  tphotphi2 = vec.Phi();
276  tphoteta2 = vec.Eta();
277  secondphotonflag = 1;
278  }
279 
280  //this is a decay to photon + ee pair i.e. dalitz decay, lets just skip it
281  if (dumparentid == 1 &&
282  firstphotonflag && secondphotonflag && secondsecondphotonflag)
283  {
284  cout << "I AM GOING TO SKIP THIS EVENT BECAUSE IT APPEARS TO BE A DALITZ DECAY" << endl;
285  cout << "THERE WERE 3 PARTICLES MARKED WITH PARENT ID == 1" << endl;
286 
287  return 0;
288  }
289 
290  //we need these extra flags because otherwise the dalitz check
291  //will always have first and secondphoton flags marked as true
292  //so this allows the dalitz check to occur after marking the
293  //second truth photon as a decay component
294  if (firstphotonflag)
295  firstfirstphotonflag = 1;
296  if (secondphotonflag)
297  secondsecondphotonflag = 1;
298  }
299  if (Verbosity() > 1)
300  {
301  cout << "truth decay photon 1 energy | phi | eta " << tphote1
302  << " | " << tphotphi1 << " | " << tphoteta1 << endl;
303  cout << "truth decay photon 2 energy | phi | eta " << tphote2
304  << " | " << tphotphi2 << " | " << tphoteta2 << endl;
305  }
306 
307 
308 
309 
310 
311 
312  /***********************************************
313 
314  GET THE FORWARD EMCAL CLUSTERS
315 
316  ************************************************/
317 
318  if (_useforwardarm)
319  {
320  RawClusterContainer::ConstRange fclus = fclusters->getClusters();
321  RawClusterContainer::ConstIterator fclusiter;
322 
323  for (fclusiter = fclus.first; fclusiter != fclus.second; ++fclusiter)
324  {
325  RawCluster *cluster = fclusiter->second;
326 
327  CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
328  CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*cluster, vertex);
329  fclusenergy = E_vec_cluster.mag();
330  fclus_eta = E_vec_cluster.pseudoRapidity();
331  fclus_theta = E_vec_cluster.getTheta();
332  fclus_pt = E_vec_cluster.perp();
333  fclus_phi = E_vec_cluster.getPhi();
334 
335 
336  if (fclusenergy < mincluspt)
337  continue;
338 
339 
340  fclus_px = fclus_pt * TMath::Cos(fclus_phi);
341  fclus_py = fclus_pt * TMath::Sin(fclus_phi);
342  fclus_pz = fclusenergy * TMath::Cos(fclus_theta);
343 
344  //reset second cluster variables
345  fclusenergy2 = -999;
346  fclus_eta2 = -999;
347  fclus_phi2 = -999;
348  fclus_theta2 = -999;
349  fclus_pt2 = -999;
350  fclus_px2 = -999;
351  fclus_py2 = -999;
352  fclus_pz2 = -999;
353 
354 
355  if(Verbosity() > 1)
356  {
357  cout<<"Found one good cluster"<<endl;
358  cout<<"energy | phi | eta: "<<fclusenergy<<" | "<<fclus_phi<<" | "<<fclus_eta<<endl;
359  }
360 
361  //found a first good cluster, lets look for a second
362  RawClusterContainer::ConstRange fclus2 = fclusters->getClusters();
363  RawClusterContainer::ConstIterator fclusiter2;
364  for (fclusiter2 = fclus2.first; fclusiter2 != fclus2.second; ++fclusiter2)
365  {
366  RawCluster *cluster2 = fclusiter2->second;
367 
368  CLHEP::Hep3Vector E_vec_cluster2 = RawClusterUtility::GetEVec(*cluster2, vertex);
369  fclusenergy2 = E_vec_cluster2.mag();
370  //just got the same photon, skip
371  if (fclusenergy == fclusenergy2)
372  {
373  //reset it
374  fclusenergy2 = -999;
375  continue;
376  }
377 
378  fclus_eta2 = E_vec_cluster2.pseudoRapidity();
379  fclus_theta2 = E_vec_cluster2.getTheta();
380  fclus_pt2 = E_vec_cluster2.perp();
381  fclus_phi2 = E_vec_cluster2.getPhi();
382 
383 
384 
385  if (fclusenergy2 < mincluspt)
386  continue;
387 
388  //avoid double counting, leading photon should always be highest energy
389  if (fclusenergy < fclusenergy2)
390  continue;
391 
392  if (Verbosity() > 1)
393  {
394  cout << "identified reco photon 1 energy | phi | eta: "
395  << fclusenergy << " | " << fclus_phi << " | " << fclus_eta << endl;
396  cout << "identified reco photon 2 energy | phi | eta: "
397  << fclusenergy2 << " | " << fclus_phi2 << " | " << fclus_eta2 << endl;
398  }
399 
400  fclus_px2 = fclus_pt2 * TMath::Cos(fclus_phi2);
401  fclus_py2 = fclus_pt2 * TMath::Sin(fclus_phi2);
402  fclus_pz2 = fclusenergy2 * TMath::Cos(fclus_theta2);
403 
404  TLorentzVector phot1, phot2;
405  phot1.SetPtEtaPhiE(fclus_pt, fclus_eta, fclus_phi, fclusenergy);
406  phot2.SetPtEtaPhiE(fclus_pt2, fclus_eta2, fclus_phi2, fclusenergy2);
407 
408  TLorentzVector pi0;
409  pi0 = phot1 + phot2;
410 
411  //get the pi0 invmass
412  invmass = pi0.M();
413 
414  if (Verbosity() > 1)
415  cout << "pi0 reco invmass is " << invmass << endl;
416  }
417 
418 
419  if(Verbosity() > 1)
420  {
421  cout<<"Final clusters found and writing to tree"<<endl;
422  cout<<"clus 1 - energy | phi | eta: "<<fclusenergy<<" | "<<fclus_phi
423  <<" | "<<fclus_eta<<endl;
424  cout<<"clus 2 - energy | phi | eta: "<<fclusenergy2<<" | "<<fclus_phi2
425  <<" | "<<fclus_eta2<<endl;
426 
427  }
428  fcluster_tree->Fill();
429  }
430  }
431 
432 
433 
434 
435 
436 
437  /***********************************************
438 
439  GET THE CENTRAL ARM EMCAL CLUSTERS
440 
441  ************************************************/
442 
443  RawClusterContainer::ConstRange begin_end = clusters->getClusters();
444  RawClusterContainer::ConstIterator clusiter;
445 
446  for (clusiter = begin_end.first; clusiter != begin_end.second; ++clusiter)
447  {
448  RawCluster *cluster = clusiter->second;
449 
450  CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
451  CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*cluster, vertex);
452  clus_energy = E_vec_cluster.mag();
453  clus_eta = E_vec_cluster.pseudoRapidity();
454  clus_theta = E_vec_cluster.getTheta();
455  clus_pt = E_vec_cluster.perp();
456  clus_phi = E_vec_cluster.getPhi();
457 
458  if (clus_pt < mincluspt)
459  continue;
460 
461  TLorentzVector *clus = new TLorentzVector();
462  clus->SetPtEtaPhiE(clus_pt, clus_eta, clus_phi, clus_energy);
463 
464  float dumarray[4];
465  clus->GetXYZT(dumarray);
466  clus_x = dumarray[0];
467  clus_y = dumarray[1];
468  clus_z = dumarray[2];
469  clus_t = dumarray[3];
470 
471  clus_px = clus_pt * TMath::Cos(clus_phi);
472  clus_py = clus_pt * TMath::Sin(clus_phi);
473  clus_pz = sqrt(clus_energy * clus_energy - clus_px * clus_px - clus_py * clus_py);
474 
475 
476 
477 
478  //reset second cluster variables
479  clus_energy2 = -999;
480  clus_eta2 = -999;
481  clus_theta2 = -999;
482  clus_pt2 = -999;
483  clus_phi2 = -999;
484  clus_px2 = -999;
485  clus_py2 = -999;
486  clus_pz2 = -999;
487 
488 
489 
490 
491 
492  //found a first good cluster, lets look for a second
493  RawClusterContainer::ConstRange begin_end2 = clusters->getClusters();
494  RawClusterContainer::ConstIterator clusiter2;
495 
496  for (clusiter2 = begin_end2.first; clusiter2 != begin_end2.second; ++clusiter2)
497  {
498  RawCluster *cluster2 = clusiter2->second;
499  CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
500  CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*cluster2, vertex);
501  clus_energy2 = E_vec_cluster.mag();
502  if(clus_energy == clus_energy2)
503  {
504  //reset it
505  clus_energy2 = -999;
506  continue;
507  }
508  clus_eta2 = E_vec_cluster.pseudoRapidity();
509  clus_theta2 = E_vec_cluster.getTheta();
510  clus_pt2 = E_vec_cluster.perp();
511  clus_phi2 = E_vec_cluster.getPhi();
512 
513  if(clus_energy == clus_energy2 and clus_eta == clus_eta2 and clus_phi == clus_phi2)
514  continue;
515 
516  if(clus_energy2 < mincluspt)
517  continue;
518 
519  //avoid double counting
520  if(clus_energy < clus_energy2)
521  continue;
522 
523  if (Verbosity() > 1)
524  {
525  cout << "CENTRAL ARM ID: "<<endl;
526  cout << "identified reco photon 1 energy | phi | eta: "
527  << clus_energy << " | " << clus_phi << " | " << clus_eta << endl;
528  cout << "identified reco photon 2 energy | phi | eta: "
529  << clus_energy2 << " | " << clus_phi2 << " | " << clus_eta2 << endl;
530  }
531 
532  clus_px2 = clus_pt2 * TMath::Cos(clus_phi2);
533  clus_py2 = clus_pt2 * TMath::Sin(clus_phi2);
534  clus_pz2 = sqrt(clus_energy2 * clus_energy2 - clus_px2 * clus_px2 - clus_py2 * clus_py2);
535 
536 
537  TLorentzVector phot1, phot2;
538  phot1.SetPtEtaPhiE(clus_pt, clus_eta, clus_phi, clus_energy);
539  phot2.SetPtEtaPhiE(clus_pt2, clus_eta2, clus_phi2, clus_energy2);
540 
541  TLorentzVector pi0;
542  pi0 = phot1 + phot2;
543 
544  //get the pi0 invmass
545  cent_invmass = pi0.M();
546 
547  if (Verbosity() > 1)
548  cout << "Central pi0 reco invmass is " << cent_invmass << endl;
549 
550 
551  cluster_tree->Fill();
552  }
553 
554  }
555 
556  tree->Fill();
557 
558 
559 
560  return 0;
561 }
562 
563 int Forward_pi0s::End(PHCompositeNode *topnode)
564 {
565  std::cout << " DONE PROCESSING " << endl;
566 
567  file->Write();
568  file->Close();
569  return 0;
570 }
571 
572 void Forward_pi0s::Set_Tree_Branches()
573 {
574  inhcal_tree = new TTree("hcalclustree", "a tree with inner hcal cluster information");
575  inhcal_tree->Branch("hcal_energy", &hcal_energy, "hcal_energy/F");
576  inhcal_tree->Branch("hcal_eta", &hcal_eta, "hcal_eta/F");
577  inhcal_tree->Branch("hcal_phi", &hcal_phi, "hcal_phi/F");
578  inhcal_tree->Branch("hcal_pt", &hcal_pt, "hcal_pt/F");
579  inhcal_tree->Branch("hcal_theta", &hcal_theta, "hcal_theta/F");
580  inhcal_tree->Branch("hcal_x", &hcal_x, "hcal_x/F");
581  inhcal_tree->Branch("hcal_y", &hcal_y, "hcal_y/F");
582  inhcal_tree->Branch("hcal_z", &hcal_z, "hcal_z/F");
583  inhcal_tree->Branch("hcal_t", &hcal_t, "hcal_t/F");
584  inhcal_tree->Branch("hcal_px", &hcal_px, "hcal_px/F");
585  inhcal_tree->Branch("hcal_py", &hcal_py, "hcal_py/F");
586  inhcal_tree->Branch("hcal_pz", &hcal_pz, "hcal_pz/F");
587  inhcal_tree->Branch("nevents", &nevents, "nevents/I");
588  inhcal_tree->Branch("hclustruthpx", &hclustruthpx, "hclustruthpx/F");
589  inhcal_tree->Branch("hclustruthpy", &hclustruthpy, "hclustruthpy/F");
590  inhcal_tree->Branch("hclustruthpz", &hclustruthpz, "hclustruthpz/F");
591  inhcal_tree->Branch("hclustruthenergy", &hclustruthenergy, "hclustruthenergy/F");
592  inhcal_tree->Branch("hclustruthpt", &hclustruthpt, "hclustruthpt/F");
593  inhcal_tree->Branch("hclustruthphi", &hclustruthphi, "hclustruthphi/F");
594  inhcal_tree->Branch("hclustrutheta", &hclustrutheta, "hclustrutheta/F");
595 
596  cluster_tree = new TTree("clustertree", "A tree with EMCal cluster information");
597  cluster_tree->Branch("clus_energy", &clus_energy, "clus_energy/F");
598  cluster_tree->Branch("clus_eta", &clus_eta, "clus_eta/F");
599  cluster_tree->Branch("clus_phi", &clus_phi, "clus_phi/F");
600  cluster_tree->Branch("clus_pt", &clus_pt, "clus_pt/F");
601  cluster_tree->Branch("clus_theta", &clus_theta, "clus_theta/F");
602  cluster_tree->Branch("clus_x", &clus_x, "clus_x/F");
603  cluster_tree->Branch("clus_y", &clus_y, "clus_y/F");
604  cluster_tree->Branch("clus_z", &clus_z, "clus_z/F");
605  cluster_tree->Branch("clus_t", &clus_t, "clus_t/F");
606  cluster_tree->Branch("clus_px", &clus_px, "clus_px/F");
607  cluster_tree->Branch("clus_py", &clus_py, "clus_py/F");
608  cluster_tree->Branch("clus_pz", &clus_pz, "clus_pz/F");
609  cluster_tree->Branch("nevents", &nevents, "nevents/I");
610 
611  cluster_tree->Branch("clus_px2", &clus_px2, "clus_px2/F");
612  cluster_tree->Branch("clus_py2", &clus_py2, "clus_py2/F");
613  cluster_tree->Branch("clus_pz2", &clus_pz2, "clus_pz2/F");
614  cluster_tree->Branch("clus_energy2", &clus_energy2, "clus_energy2/F");
615  cluster_tree->Branch("clus_eta2", &clus_eta2, "clus_eta2/F");
616  cluster_tree->Branch("clus_phi2", &clus_phi2, "clus_phi2/F");
617  cluster_tree->Branch("clus_pt2", &clus_pt2, "clus_pt2/F");
618  cluster_tree->Branch("clus_theta2", &clus_theta2, "clus_theta2/F");
619  cluster_tree->Branch("Cent_invmass", &cent_invmass,"cent_invmass/F");
620 
621  cluster_tree->Branch("tpi0e", &tpi0e, "tpi0e/F");
622  cluster_tree->Branch("tpi0px", &tpi0px, "tpi0px/F");
623  cluster_tree->Branch("tpi0py", &tpi0py, "tpi0py/F");
624  cluster_tree->Branch("tpi0pz", &tpi0pz, "tpi0pz/F");
625  cluster_tree->Branch("tpi0pid", &tpi0pid, "tpi0pid/F");
626  cluster_tree->Branch("tpi0pt", &tpi0pt, "tpi0pt/F");
627  cluster_tree->Branch("tpi0phi", &tpi0phi, "tpi0phi/F");
628  cluster_tree->Branch("tpi0eta", &tpi0eta, "tpi0eta/F");
629  cluster_tree->Branch("tphote1", &tphote1, "tphote1/F");
630  cluster_tree->Branch("tphotpx1", &tphotpx1, "tphotpx1/F");
631  cluster_tree->Branch("tphotpy1", &tphotpy1, "tphotpy1/F");
632  cluster_tree->Branch("tphotpz1", &tphotpz1, "tphotpz1/F");
633  cluster_tree->Branch("tphotpt1", &tphotpt1, "tphotpt1/F");
634  cluster_tree->Branch("tphotpid1", &tphotpid1, "tphotpid1/I");
635  cluster_tree->Branch("tphotparentid1", &tphotparentid1, "tphotparentid1/I");
636  cluster_tree->Branch("tphotphi1", &tphotphi1, "tphotphi1/F");
637  cluster_tree->Branch("tphoteta1", &tphoteta1, "tphoteta1/F");
638 
639  cluster_tree->Branch("tphote2", &tphote2, "tphote2/F");
640  cluster_tree->Branch("tphotpx2", &tphotpx2, "tphotpx2/F");
641  cluster_tree->Branch("tphotpy2", &tphotpy2, "tphotpy2/F");
642  cluster_tree->Branch("tphotpz2", &tphotpz2, "tphotpz2/F");
643  cluster_tree->Branch("tphotpt2", &tphotpt2, "tphotpt2/F");
644  cluster_tree->Branch("tphotpid2", &tphotpid2, "tphotpid2/I");
645  cluster_tree->Branch("tphotparentid2", &tphotparentid2, "tphotparentid2/I");
646  cluster_tree->Branch("tphotphi2", &tphotphi2, "tphotphi2/F");
647  cluster_tree->Branch("tphoteta2", &tphoteta2, "tphoteta2/F");
648 
649 
650 
651  fcluster_tree = new TTree("fclustertree", "A tree with FEMCal cluster information");
652  fcluster_tree->Branch("invmass", &invmass, "invmass/F");
653  fcluster_tree->Branch("fclusenergy", &fclusenergy, "fclusenergy/F");
654  fcluster_tree->Branch("fclus_eta", &fclus_eta, "fclus_eta/F");
655  fcluster_tree->Branch("fclus_phi", &fclus_phi, "fclus_phi/F");
656  fcluster_tree->Branch("fclus_pt", &fclus_pt, "fclus_pt/F");
657  fcluster_tree->Branch("fclus_theta", &fclus_theta, "fclus_theta/F");
658  fcluster_tree->Branch("fclus_px", &fclus_px, "fclus_px/F");
659  fcluster_tree->Branch("fclus_py", &fclus_py, "fclus_py/F");
660  fcluster_tree->Branch("fclus_pz", &fclus_pz, "fclus_pz/F");
661  fcluster_tree->Branch("nevents", &nevents, "nevents/I");
662  fcluster_tree->Branch("fclusenergy2", &fclusenergy2, "fclusenergy2/F");
663  fcluster_tree->Branch("fclus_eta2", &fclus_eta2, "fclus_eta2/F");
664  fcluster_tree->Branch("fclus_phi2", &fclus_phi2, "fclus_phi2/F");
665  fcluster_tree->Branch("fclus_pt2", &fclus_pt2, "fclus_pt2/F");
666  fcluster_tree->Branch("fclus_theta2", &fclus_theta2, "fclus_theta2/F");
667  fcluster_tree->Branch("fclus_px2", &fclus_px2, "fclus_px2/F");
668  fcluster_tree->Branch("fclus_py2", &fclus_py2, "fclus_py2/F");
669  fcluster_tree->Branch("fclus_pz2", &fclus_pz2, "fclus_pz2/F");
670 
671  fcluster_tree->Branch("tpi0e", &tpi0e, "tpi0e/F");
672  fcluster_tree->Branch("tpi0px", &tpi0px, "tpi0px/F");
673  fcluster_tree->Branch("tpi0py", &tpi0py, "tpi0py/F");
674  fcluster_tree->Branch("tpi0pz", &tpi0pz, "tpi0pz/F");
675  fcluster_tree->Branch("tpi0pid", &tpi0pid, "tpi0pid/F");
676  fcluster_tree->Branch("tpi0pt", &tpi0pt, "tpi0pt/F");
677  fcluster_tree->Branch("tpi0phi", &tpi0phi, "tpi0phi/F");
678  fcluster_tree->Branch("tpi0eta", &tpi0eta, "tpi0eta/F");
679  fcluster_tree->Branch("tphote1", &tphote1, "tphote1/F");
680  fcluster_tree->Branch("tphotpx1", &tphotpx1, "tphotpx1/F");
681  fcluster_tree->Branch("tphotpy1", &tphotpy1, "tphotpy1/F");
682  fcluster_tree->Branch("tphotpz1", &tphotpz1, "tphotpz1/F");
683  fcluster_tree->Branch("tphotpt1", &tphotpt1, "tphotpt1/F");
684  fcluster_tree->Branch("tphotpid1", &tphotpid1, "tphotpid1/I");
685  fcluster_tree->Branch("tphotparentid1", &tphotparentid1, "tphotparentid1/I");
686  fcluster_tree->Branch("tphotphi1", &tphotphi1, "tphotphi1/F");
687  fcluster_tree->Branch("tphoteta1", &tphoteta1, "tphoteta1/F");
688 
689  fcluster_tree->Branch("tphote2", &tphote2, "tphote2/F");
690  fcluster_tree->Branch("tphotpx2", &tphotpx2, "tphotpx2/F");
691  fcluster_tree->Branch("tphotpy2", &tphotpy2, "tphotpy2/F");
692  fcluster_tree->Branch("tphotpz2", &tphotpz2, "tphotpz2/F");
693  fcluster_tree->Branch("tphotpt2", &tphotpt2, "tphotpt2/F");
694  fcluster_tree->Branch("tphotpid2", &tphotpid2, "tphotpid2/I");
695  fcluster_tree->Branch("tphotparentid2", &tphotparentid2, "tphotparentid2/I");
696  fcluster_tree->Branch("tphotphi2", &tphotphi2, "tphotphi2/F");
697  fcluster_tree->Branch("tphoteta2", &tphoteta2, "tphoteta2/F");
698 
699  truth_g4particles = new TTree("truthtree_g4", "a tree with all truth g4 particles");
700  truth_g4particles->Branch("truthpx", &truthpx, "truthpx/F");
701  truth_g4particles->Branch("truthpy", &truthpy, "truthpy/F");
702  truth_g4particles->Branch("truthpz", &truthpz, "truthpz/F");
703  truth_g4particles->Branch("truthp", &truthp, "truthp/F");
704  truth_g4particles->Branch("truthenergy", &truthenergy, "truthenergy/F");
705  truth_g4particles->Branch("truthphi", &truthphi, "truthphi/F");
706  truth_g4particles->Branch("trutheta", &trutheta, "trutheta/F");
707  truth_g4particles->Branch("truthpt", &truthpt, "truthpt/F");
708  truth_g4particles->Branch("truthpid", &truthpid, "truthpid/I");
709  truth_g4particles->Branch("nevents", &nevents, "nevents/I");
710 
711  truth_g4particles->Branch("tphote1", &tphote1, "tphote1/F");
712  truth_g4particles->Branch("tphotpx1", &tphotpx1, "tphotpx1/F");
713  truth_g4particles->Branch("tphotpy1", &tphotpy1, "tphotpy1/F");
714  truth_g4particles->Branch("tphotpz1", &tphotpz1, "tphotpz1/F");
715  truth_g4particles->Branch("tphotpt1", &tphotpt1, "tphotpt1/F");
716  truth_g4particles->Branch("tphotpid1", &tphotpid1, "tphotpid1/I");
717  truth_g4particles->Branch("tphotparentid1", &tphotparentid1, "tphotparentid1/I");
718  truth_g4particles->Branch("tphotphi1", &tphotphi1, "tphotphi1/F");
719  truth_g4particles->Branch("tphoteta1", &tphoteta1, "tphoteta1/F");
720 
721  truth_g4particles->Branch("tphote2", &tphote2, "tphote2/F");
722  truth_g4particles->Branch("tphotpx2", &tphotpx2, "tphotpx2/F");
723  truth_g4particles->Branch("tphotpy2", &tphotpy2, "tphotpy2/F");
724  truth_g4particles->Branch("tphotpz2", &tphotpz2, "tphotpz2/F");
725  truth_g4particles->Branch("tphotpt2", &tphotpt2, "tphotpt2/F");
726  truth_g4particles->Branch("tphotpid2", &tphotpid2, "tphotpid2/I");
727  truth_g4particles->Branch("tphotparentid2", &tphotparentid2, "tphotparentid2/I");
728  truth_g4particles->Branch("tphotphi2", &tphotphi2, "tphotphi2/F");
729  truth_g4particles->Branch("tphoteta2", &tphoteta2, "tphoteta2/F");
730 }
731 
732 void Forward_pi0s::initialize_things()
733 {
734  hcal_energy = -999;
735  hcal_eta = -999;
736  hcal_phi = -999;
737  hcal_pt = -999;
738  hcal_px = -999;
739  hcal_py = -999;
740  hcal_pz = -999;
741  hcal_theta = -999;
742  hcal_x = -999;
743  hcal_y = -999;
744  hcal_z = -999;
745  hcal_t = -999;
746 
747  clus_energy = -999;
748  clus_eta = -999;
749  clus_phi = -999;
750  clus_pt = -999;
751  clus_px = -999;
752  clus_py = -999;
753  clus_pz = -999;
754  clus_theta = -999;
755  clus_x = -999;
756  clus_y = -999;
757  clus_z = -999;
758  clus_t = -999;
759 
760  fclusenergy2 = -999;
761  fclus_eta2 = -999;
762  fclus_phi2 = -999;
763  fclus_theta2 = -999;
764  fclus_pt2 = -999;
765  fclus_px2 = -999;
766  fclus_py2 = -999;
767  fclus_pz2 = -999;
768  invmass = -999;
769 
770  tphote1 = -999;
771  tphotpx1 = -999;
772  tphotpy1 = -999;
773  tphotpz1 = -999;
774  tphotpid1 = -999;
775  tphotparentid1 = -999;
776  tphotpt1 = -999;
777  tphotphi1 = -999;
778  tphoteta1 = -999;
779  tphote2 = -999;
780  tphotpx2 = -999;
781  tphotpy2 = -999;
782  tphotpz2 = -999;
783  tphotpid2 = -999;
784  tphotparentid2 = -999;
785  tphotpt2 = -999;
786  tphotphi2 = -999;
787  tphoteta2 = -999;
788  fclusenergy = -999;
789  fclus_eta = -999;
790  fclus_phi = -999;
791  fclus_theta = -999;
792  fclus_pt = -999;
793  fclustruthpid = -999;
794  fclustruthpx = -999;
795  fclustruthpy = -999;
796  fclustruthpz = -999;
797  fclustruthenergy = -999;
798  fclustruthpt = -999;
799  fclustruthphi = -999;
800  fclustrutheta = -999;
801  fclus_px = -999;
802  fclus_py = -999;
803  fclus_pz = -999;
804 
805  tpi0e = -999;
806  tpi0px = -999;
807  tpi0py = -999;
808  tpi0pz = -999;
809  tpi0pid = -999;
810  tpi0pt = -999;
811  tpi0phi = -999;
812  tpi0eta = -999;
813 }