Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pi0ClusterAna.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pi0ClusterAna.cc
1 #include "pi0ClusterAna.h"
2 
3 //Fun4all stuff
4 #include <fun4all/Fun4AllReturnCodes.h>
5 #include <fun4all/Fun4AllServer.h>
6 #include <fun4all/Fun4AllHistoManager.h>
7 #include <phool/PHCompositeNode.h>
8 #include <phool/getClass.h>
9 #include <phool/phool.h>
10 #include <ffaobjects/EventHeader.h>
11 
12 //ROOT stuff
13 #include <TH1F.h>
14 #include <TH2F.h>
15 #include <TH3F.h>
16 #include <TFile.h>
17 #include <TLorentzVector.h>
18 #include <TTree.h>
19 
20 //for emc clusters
21 #include <calobase/RawCluster.h>
22 #include <calobase/RawClusterContainer.h>
23 #include <calobase/RawClusterUtility.h>
24 #include <calobase/RawTowerGeomContainer.h>
25 #include <calobase/RawTower.h>
26 #include <calobase/RawTowerContainer.h>
27 
28 //for vetex information
29 #include <g4vertex/GlobalVertex.h>
30 #include <g4vertex/GlobalVertexMap.h>
31 
32 //tracking
33 #include <trackbase_historic/SvtxTrack.h>
34 #include <trackbase_historic/SvtxTrackMap.h>
35 #include <trackbase_historic/SvtxVertex.h>
36 #include <trackbase_historic/SvtxVertexMap.h>
37 
38 //truth information
39 #include <g4main/PHG4TruthInfoContainer.h>
40 #include <g4main/PHG4Particle.h>
41 #include <g4main/PHG4VtxPoint.h>
42 #pragma GCC diagnostic push
43 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
44 #include <HepMC/GenEvent.h>
45 #include <HepMC/GenParticle.h>
46 #include <HepMC/GenVertex.h>
47 #include <HepMC/IteratorRange.h>
48 #include <HepMC/SimpleVector.h>
49 #include <HepMC/GenParticle.h>
50 #pragma GCC diagnostic pop
51 #include <CLHEP/Geometry/Point3D.h>
52 
53 
54 //____________________________________________________________________________..
55 pi0ClusterAna::pi0ClusterAna(const std::string &name, const std::string &outName = "pi0ClusterAnaOut"):
56 SubsysReco(name)
57  , clusters_Towers(nullptr)
58  , truth_photon(nullptr)
59  , truth_pi0(nullptr)
60  //, caloevalstack(NULL)
61  , m_eta_center()
62  , m_phi_center()
63  , m_tower_energy()
64 , m_cluster_eta()
65  , m_cluster_phi()
66 , m_cluster_e()
67  , m_cluster_chi2()
68  , m_cluster_prob()
69 , m_cluster_nTowers()
70  , m_asym()
71  , m_deltaR()
72  , m_lead_E()
73 , m_sublead_E()
74 , m_lead_phi()
75 , m_lead_eta()
76  , m_sublead_phi()
77  , m_sublead_eta()
78  , m_pi0_E()
79 , m_pi0_eta()
80 , m_pi0_phi()
81  , Outfile(outName)
82 
83 {
84 
85 
86  std::cout << "pi0ClusterAna::pi0ClusterAna(const std::string &name) Calling ctor" << std::endl;
87 }
88 
89 //____________________________________________________________________________..
90 pi0ClusterAna::~pi0ClusterAna()
91 {
92  std::cout << "pi0ClusterAna::~pi0ClusterAna() Calling dtor" << std::endl;
93 }
94 
95 //____________________________________________________________________________..
96 int pi0ClusterAna::Init(PHCompositeNode *topNode)
97 {
98 
99  clusters_Towers = new TTree("TreeClusterTower","Tree for cluster and tower information");
100  clusters_Towers -> Branch("towerEtaCEnter",&m_eta_center);
101  clusters_Towers -> Branch("towerPhiCenter",& m_phi_center);
102  clusters_Towers -> Branch("towerEnergy",&m_tower_energy);
103  clusters_Towers -> Branch("clusterEta",&m_cluster_eta);
104  clusters_Towers -> Branch("clusterPhi", &m_cluster_phi);
105  clusters_Towers -> Branch("clusterE",&m_cluster_e);
106  clusters_Towers -> Branch("clusterChi2",&m_cluster_chi2);
107  clusters_Towers -> Branch("clusterProb",&m_cluster_prob);
108  clusters_Towers -> Branch("clusterNTowers",&m_cluster_nTowers);
109 
110  truth_photon = new TTree("TreeTruthPhoton", "Tree for truth photons");
111  truth_photon -> Branch("pairAsym",&m_asym);
112  truth_photon -> Branch("pairDeltaR",&m_deltaR);
113  truth_photon -> Branch("leadPhotonPhi",&m_lead_phi);
114  truth_photon -> Branch("leadPhotonEta",&m_lead_eta);
115  truth_photon -> Branch("subleadPhotonPhi",&m_sublead_phi);
116  truth_photon -> Branch("subleadPhotonEta",&m_sublead_eta);
117  truth_photon -> Branch("leadPhotonE", &m_lead_E);
118  truth_photon -> Branch("subLeadPhotonE", &m_sublead_E);
119 
120  truth_pi0 = new TTree("TreeTruthPi0", "Tree for Truth pi0");
121  truth_pi0 -> Branch("pi0E",&m_pi0_E);
122  truth_pi0 -> Branch("pi0_phi",&m_pi0_phi);
123  truth_pi0 -> Branch("pi0_eta",&m_pi0_eta);
124 
125  //so that the histos actually get written out
126  Fun4AllServer *se = Fun4AllServer::instance();
127  se -> Print("NODETREE");
128  hm = new Fun4AllHistoManager("MYHISTOS");
129 
130  se -> registerHistoManager(hm);
131 
132  se -> registerHisto(truth_pi0 -> GetName(),truth_pi0);
133  se -> registerHisto(truth_photon -> GetName(), truth_photon);
134  se -> registerHisto(clusters_Towers -> GetName(), clusters_Towers);
135 
136  out = new TFile(Outfile.c_str(),"RECREATE");
137 
138  std::cout << "pi0ClusterAna::Init(PHCompositeNode *topNode) Initializing" << std::endl;
139  return Fun4AllReturnCodes::EVENT_OK;
140 }
141 
142 //____________________________________________________________________________..
143 int pi0ClusterAna::InitRun(PHCompositeNode *topNode)
144 {
145  std::cout << "pi0ClusterAna::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
146 
147  return Fun4AllReturnCodes::EVENT_OK;
148 }
149 
150 //____________________________________________________________________________..
151 int pi0ClusterAna::process_event(PHCompositeNode *topNode)
152 {
153 
154  //Information on clusters
155  RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_POS_COR_CEMC");
156  //RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
157  if(!clusterContainer)
158  {
159  std::cout << PHWHERE << "pi0Efficiency::process_event - Fatal Error - CLUSTER_CEMC node is missing. " << std::endl;
160 
161  return 0;
162  }
163 
164 
165  //Vertex information
166  GlobalVertexMap *vtxContainer = findNode::getClass<GlobalVertexMap>(topNode,"GlobalVertexMap");
167  if (!vtxContainer)
168  {
169  std::cout << PHWHERE << "pi0Efficiency::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." << std::endl;
170  assert(vtxContainer); // force quit
171 
172  return 0;
173  }
174 
175  if (vtxContainer->empty())
176  {
177  std::cout << PHWHERE << "pi0Efficiency::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." << std::endl;
178  return 0;
179  }
180 
181  //More vertex information
182  GlobalVertex *vtx = vtxContainer->begin()->second;
183  if(!vtx)
184  {
185 
186  std::cout << PHWHERE << "pi0Efficiency::process_event Could not find vtx from vtxContainer" << std::endl;
187  return Fun4AllReturnCodes::ABORTEVENT;
188  }
189 
190  //Tower geometry node for location information
191  RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
192  if (!towergeom)
193  {
194  std::cout << PHWHERE << "pi0Efficiency::process_event Could not find node TOWERGEOM_CEMC" << std::endl;
195  return Fun4AllReturnCodes::ABORTEVENT;
196  }
197 
198  //Raw tower information
199  RawTowerContainer *towerContainer = findNode::getClass<RawTowerContainer>(topNode,"TOWER_CALIB_CEMC");
200  if(!towerContainer)
201  {
202  std::cout << PHWHERE << "pi0Efficiency::process_event Could not find node TOWER_CALIB_CEMC" << std::endl;
203  return Fun4AllReturnCodes::ABORTEVENT;
204  }
205 
206  //truth particle information
207  PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
208  if(!truthinfo)
209  {
210  std::cout << PHWHERE << "pi0Efficiency::process_event Could not find node G4TruthInfo" << std::endl;
211  return Fun4AllReturnCodes::ABORTEVENT;
212  }
213 
214 
215  float pi0Eta = -99999;
216  float photonEtaMax = 1.1;
217  float mesonEtaMax = 0.3;
218  //Now we go and find our truth pi0s and just take its eta coordinate for now
219  PHG4TruthInfoContainer::Range truthRange = truthinfo->GetPrimaryParticleRange();
220  PHG4TruthInfoContainer::ConstIterator truthIter;
221 
222  PHG4Particle *truthPar = NULL;
223 
224  for(truthIter = truthRange.first; truthIter != truthRange.second; truthIter++)
225  {
226  truthPar = truthIter->second;
227 
228  if(truthPar -> get_pid() == 111 && truthPar -> get_parent_id() == 0)
229  {
230  pi0Eta = getEta(truthPar);
231  if(abs(pi0Eta) >= mesonEtaMax)
232  {
233  return 0;
234  }
235  }
236  }
237 
238  int firstphotonflag = 0;
239  int firstfirstphotonflag = 0;
240  int secondphotonflag = 0;
241 
242  //int secondsecondphotonflag = 0;
243 
244  PHG4TruthInfoContainer::Range truthRangeDecay1 = truthinfo->GetSecondaryParticleRange();
245  PHG4TruthInfoContainer::ConstIterator truthIterDecay1;
246 
247 
248  TLorentzVector photon1;
249  TLorentzVector photon2;
250  int nParticles = 0;
251 
252  //loop over all our decay photons.
253  //Make sure they fall within the desired acceptance
254  //Toss Dalitz decays
255  //Retain photon kinematics to determine lead photon for eta binning
256  for(truthIterDecay1 = truthRangeDecay1.first; truthIterDecay1 != truthRangeDecay1.second; truthIterDecay1++)
257  {
258  PHG4Particle *decay = truthIterDecay1 -> second;
259 
260  int dumtruthpid = decay -> get_pid();
261  int dumparentid = decay -> get_parent_id();
262 
263  //if the parent is the pi0 and it's a photon and we haven't marked one yet
264  if(dumparentid == 1 && dumtruthpid == 22 && !firstphotonflag)
265  {
266  if(abs(getEta(decay)) > photonEtaMax)
267  {
268  return 0;
269  }
270  photon1.SetPxPyPzE(decay -> get_px(), decay -> get_py(), decay -> get_pz(), decay -> get_e());
271 
272  firstphotonflag = 1;
273  }
274 
275  if(dumparentid == 1 && dumtruthpid == 22 && firstphotonflag && firstfirstphotonflag)
276  {
277  if(abs(getEta(decay)) > photonEtaMax)
278  {
279  return 0;
280  }
281  photon2.SetPxPyPzE(decay -> get_px(), decay -> get_py(), decay -> get_pz(), decay -> get_e()) ;
282 
283  secondphotonflag = 1;
284  }
285 
286  //Need this flag to make it skip the first photon slot
287  if(firstphotonflag) firstfirstphotonflag = 1;
288  if(dumparentid == 1)nParticles ++;
289  }
290 
291  if((!firstphotonflag || !secondphotonflag) && nParticles > 1) //Dalitz
292  {
293  return 0;
294  }
295  else if((!firstphotonflag || !secondphotonflag) && nParticles == 1) //One photon falls outside simulation acceptance
296  {
297  return 0;
298  }
299  else if((!firstphotonflag || !secondphotonflag) && nParticles == 0) //Both photons fall outside simulation acceptance
300  {
301  return 0;
302  }
303 
304 
305  TLorentzVector pi0;
306  pi0.SetPxPyPzE(truthPar -> get_px(), truthPar -> get_py(), truthPar -> get_pz(), truthPar -> get_e());
307  TLorentzVector leadPho, subLeadPho;
308  if(photon1.Energy() >= photon2.Energy())
309  {
310  leadPho = photon1;
311  subLeadPho = photon2;
312  }
313  else
314  {
315  leadPho = photon2;
316  subLeadPho = photon1;
317  }
318 
319  float asym = abs(photon1.Energy() - photon2.Energy())/(photon1.Energy() + photon2.Energy());
320 
321  m_asym.push_back(asym);
322 
323  float deltaR = photon1.DeltaR(photon2);
324 
325  m_deltaR.push_back(deltaR);
326 
327  m_lead_phi.push_back(leadPho.Phi());
328  m_sublead_phi.push_back(subLeadPho.Phi());
329 
330  m_lead_eta.push_back(leadPho.PseudoRapidity());
331  m_sublead_eta.push_back(subLeadPho.PseudoRapidity());
332 
333  m_lead_E.push_back(leadPho.Energy());
334  m_sublead_E.push_back(subLeadPho.Energy());
335 
336  m_pi0_E.push_back(pi0.Energy());
337  m_pi0_phi.push_back(pi0.Phi());
338  m_pi0_eta.push_back(pi0.PseudoRapidity());
339 
340  //grab all the towers and fill their energies.
341  RawTowerContainer::ConstRange tower_range = towerContainer -> getTowers();
342  for(RawTowerContainer::ConstIterator tower_iter = tower_range.first; tower_iter!= tower_range.second; tower_iter++)
343  {
344  int phibin = tower_iter -> second -> get_binphi();
345  int etabin = tower_iter -> second -> get_bineta();
346  double phi = towergeom -> get_phicenter(phibin);
347  double eta = towergeom -> get_etacenter(etabin);
348  double energy = tower_iter -> second -> get_energy();
349 
350  m_phi_center.push_back(phi);
351  m_eta_center.push_back(eta);
352  m_tower_energy.push_back(energy);
353  }
354 
355  //This is how we iterate over items in the container.
356  RawClusterContainer::ConstRange clusterEnd = clusterContainer -> getClusters();
357  RawClusterContainer::ConstIterator clusterIter;
358 
359  for(clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
360  {
361  RawCluster *recoCluster = clusterIter -> second;
362 
363  CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
364  CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
365  float clusE = E_vec_cluster.mag();
366  float clus_eta = E_vec_cluster.pseudoRapidity();
367  float clus_phi = E_vec_cluster.phi();
368 
369  m_cluster_eta.push_back(clus_eta);
370  m_cluster_phi.push_back(clus_phi);
371  m_cluster_e.push_back(clusE);
372 
373  m_cluster_chi2.push_back(recoCluster -> get_chi2());
374  m_cluster_prob.push_back(recoCluster -> get_prob());
375  m_cluster_nTowers.push_back(recoCluster -> getNTowers());
376  }
377 
378 
379  clusters_Towers -> Fill();
380  truth_photon -> Fill();
381  truth_pi0 -> Fill();
382 
383 
384  return Fun4AllReturnCodes::EVENT_OK;
385 }
386 
387 //____________________________________________________________________________..
388 int pi0ClusterAna::ResetEvent(PHCompositeNode *topNode)
389 {
390  //std::cout << "pi0ClusterAna::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
391 
392  m_eta_center.clear();
393  m_phi_center.clear();
394  m_tower_energy.clear();
395  m_cluster_eta.clear();
396  m_cluster_phi.clear();
397  m_cluster_e.clear();
398  m_cluster_chi2.clear();
399  m_cluster_prob.clear();
400  m_cluster_nTowers.clear();
401  m_asym.clear();
402  m_deltaR.clear();
403  m_lead_E.clear();
404  m_sublead_E.clear();
405  m_lead_phi.clear();
406  m_lead_eta.clear();
407  m_sublead_phi.clear();
408  m_sublead_eta.clear();
409  m_pi0_E.clear();
410  m_pi0_eta.clear();
411  m_pi0_phi.clear();
412  return Fun4AllReturnCodes::EVENT_OK;
413 }
414 
415 //____________________________________________________________________________..
416 int pi0ClusterAna::EndRun(const int runnumber)
417 {
418  std::cout << "pi0ClusterAna::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
419  return Fun4AllReturnCodes::EVENT_OK;
420 }
421 
422 //____________________________________________________________________________..
423 int pi0ClusterAna::End(PHCompositeNode *topNode)
424 {
425  std::cout << "pi0ClusterAna::End(PHCompositeNode *topNode) This is the End..." << std::endl;
426  out -> cd();
427 
428  truth_pi0 -> Write();
429  truth_photon -> Write();
430  clusters_Towers -> Write();
431 
432  out -> Close();
433  delete out;
434 
435  hm -> dumpHistos(Outfile.c_str(),"UPDATE");
436 
437 
438  return Fun4AllReturnCodes::EVENT_OK;
439 }
440 
441 //____________________________________________________________________________..
442 int pi0ClusterAna::Reset(PHCompositeNode *topNode)
443 {
444  std::cout << "pi0ClusterAna::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
445  return Fun4AllReturnCodes::EVENT_OK;
446 }
447 
448 //____________________________________________________________________________..
449 void pi0ClusterAna::Print(const std::string &what) const
450 {
451  std::cout << "pi0ClusterAna::Print(const std::string &what) const Printing info for " << what << std::endl;
452 }
453 //____________________________________________________________________________..
454 float pi0ClusterAna::getEta(PHG4Particle *particle)
455 {
456  float px = particle -> get_px();
457  float py = particle -> get_py();
458  float pz = particle -> get_pz();
459  float p = sqrt(pow(px,2) + pow(py,2) + pow(pz,2));
460 
461  return 0.5*log((p+pz)/(p-pz));
462 }