Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
isoCluster.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file isoCluster.cc
1 
2 //____________________________________________________________________________..
3 
4 //our base code
5 #include "isoCluster.h"
6 
7 #include <fun4all/Fun4AllReturnCodes.h>
8 
9 #include <phool/PHCompositeNode.h>
10 
11 //Fun4all stuff
12 #include <fun4all/Fun4AllReturnCodes.h>
13 #include <fun4all/Fun4AllServer.h>
14 #include <fun4all/Fun4AllHistoManager.h>
15 #include <phool/PHCompositeNode.h>
16 #include <phool/getClass.h>
17 #include <phool/phool.h>
18 #include <ffaobjects/EventHeader.h>
19 
20 //ROOT stuff
21 #include <TH1F.h>
22 #include <TH2F.h>
23 #include <TH3F.h>
24 #include <TFile.h>
25 #include <TLorentzVector.h>
26 #include <TTree.h>
27 
28 //for emc clusters
29 #include <calobase/RawCluster.h>
30 #include <calobase/RawClusterContainer.h>
31 #include <calobase/RawClusterUtility.h>
32 #include <calobase/RawTowerGeomContainer.h>
33 #include <calobase/RawTower.h>
34 #include <calobase/RawTowerContainer.h>
35 
36 //caloEvalStack for cluster to truth matching
37 #include <g4eval/CaloEvalStack.h>
38 #include <g4eval/CaloRawClusterEval.h>
39 
40 //for vetex information
41 #include <g4vertex/GlobalVertex.h>
42 #include <g4vertex/GlobalVertexMap.h>
43 
44 //tracking
45 #include <trackbase_historic/SvtxTrack.h>
46 #include <trackbase_historic/SvtxTrackMap.h>
47 #include <trackbase_historic/SvtxVertex.h>
48 #include <trackbase_historic/SvtxVertexMap.h>
49 
50 //truth information
51 #include <g4main/PHG4TruthInfoContainer.h>
52 #include <g4main/PHG4Particle.h>
53 #include <g4main/PHG4VtxPoint.h>
54 #include <phhepmc/PHHepMCGenEvent.h>
55 #include <phhepmc/PHHepMCGenEventMap.h>
56 #pragma GCC diagnostic push
57 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
58 #include <HepMC/GenEvent.h>
59 #include <HepMC/GenParticle.h>
60 #include <HepMC/GenVertex.h>
61 #include <HepMC/IteratorRange.h>
62 #include <HepMC/SimpleVector.h>
63 #include <HepMC/GenParticle.h>
64 #pragma GCC diagnostic pop
65 
66 //____________________________________________________________________________..
67 isoCluster::isoCluster(const std::string &name):
68 SubsysReco(name),
69  m_clusterPhi(0),
70 m_clusterEta(0),
71  m_clusterE(0),
72 m_clusterPt(0),
73  m_clusterEtIso(0),
74  m_clusterChisq(0),
75  m_photonPhi(0),
76  m_photonEta(0),
77  m_photonE(0),
78 m_photonPt(0),
79 fout(NULL),
80  outname(name),
81  getEvent(-9999)
82 {
83 std::cout << "isoCluster::isoCluster(const std::string &name) Calling ctor" << std::endl;
84 }
85 //____________________________________________________________________________..
86 isoCluster::~isoCluster()
87 {
88  std::cout << "isoCluster::~isoCluster() Calling dtor" << std::endl;
89 }
90 
91 //____________________________________________________________________________..
92 int isoCluster::Init(PHCompositeNode *topNode)
93 {
94  std::cout << "isoCluster::Init(PHCompositeNode *topNode) Initializing" << std::endl;
95 
96  fout = new TFile(outname.c_str(),"RECREATE");
97 
98  T = new TTree("T","T");
99 
100  T -> Branch("clusterPhi",&m_clusterPhi);
101  T -> Branch("clusterEta",&m_clusterEta);
102  T -> Branch("clusterE",&m_clusterE);
103  T -> Branch("clusterPt",&m_clusterPt);
104  T -> Branch("clusterEtIso",&m_clusterEtIso);
105  T -> Branch("clusterchisq",&m_clusterChisq);
106 
107  T -> Branch("photonPhi",&m_photonPhi);
108  T -> Branch("photonEta",&m_photonEta);
109  T -> Branch("photonE",&m_photonE);
110  T -> Branch("photonPt",&m_photonPt);
111 
112  return Fun4AllReturnCodes::EVENT_OK;
113 }
114 
115 //____________________________________________________________________________..
116 int isoCluster::InitRun(PHCompositeNode *topNode)
117 {
118  std::cout << "isoCluster::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
119  return Fun4AllReturnCodes::EVENT_OK;
120 }
121 
122 //____________________________________________________________________________..
123 int isoCluster::process_event(PHCompositeNode *topNode)
124 {
125  //Information on clusters
126  // RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_POS_COR_CEMC");
127  RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
128  if(!clusterContainer)
129  {
130  std::cout << PHWHERE << "isoCluster::process_event - Fatal Error - CLUSTER_POS_COR_CEMC node is missing. " << std::endl;
131 
132  return 0;
133  }
134 
135  //Vertex information
136  GlobalVertexMap *vtxContainer = findNode::getClass<GlobalVertexMap>(topNode,"GlobalVertexMap");
137  if (!vtxContainer)
138  {
139  std::cout << PHWHERE << "isoCluster::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;
140  assert(vtxContainer); // force quit
141 
142  return 0;
143  }
144 
145  if (vtxContainer->empty())
146  {
147  std::cout << PHWHERE << "isoCluster::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;
148  return 0;
149  }
150 
151  //More vertex information
152  GlobalVertex *vtx = vtxContainer->begin()->second;
153  if(!vtx)
154  {
155 
156  std::cout << PHWHERE << "isoCluster::process_event Could not find vtx from vtxContainer" << std::endl;
157  return Fun4AllReturnCodes::ABORTEVENT;
158  }
159 
160  //truth particle information
161  PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
162  if(!truthinfo)
163  {
164  std::cout << PHWHERE << "isoCluster::process_event Could not find node G4TruthInfo" << std::endl;
165  return Fun4AllReturnCodes::ABORTEVENT;
166  }
167 
168  //For eventgen ancestory information
169  PHHepMCGenEventMap *genEventMap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
170  if(!genEventMap)
171  {
172  std::cout << PHWHERE << "isoCluster::process_event Could not find PHHepMCGenEventMap" << std::endl;
173  return Fun4AllReturnCodes::ABORTEVENT;
174  }
175 
176  //event level information of the above
177  PHHepMCGenEvent *genEvent = genEventMap -> get(getEvent);
178  if(!genEvent)
179  {
180  std::cout << PHWHERE << "isoCluster::process_event Could not find PHHepMCGenEvent" << std::endl;
181  return Fun4AllReturnCodes::ABORTEVENT;
182  }
183 
184  HepMC::GenEvent *theEvent = genEvent -> getEvent();
185 
186  CaloEvalStack caloevalstack(topNode, "CEMC");
187  CaloRawClusterEval *clustereval = caloevalstack.get_rawcluster_eval();
188  if(!clustereval)
189  {
190  std::cout << PHWHERE << "isoCluster::process_event cluster eval not available" << std::endl;
191  return Fun4AllReturnCodes::ABORTEVENT;
192  }
193 
194  //grab all the clusters in the event
195  RawClusterContainer::ConstRange clusterEnd = clusterContainer -> getClusters();
196  RawClusterContainer::ConstIterator clusterIter;
197  for(clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
198  {
199  RawCluster *recoCluster = clusterIter -> second;
200  CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
201  CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
202  if(E_vec_cluster.mag() < 3) continue;
203  m_clusterE.push_back(E_vec_cluster.mag());
204  m_clusterEta.push_back(E_vec_cluster.pseudoRapidity());
205  m_clusterPhi.push_back(E_vec_cluster.phi());
206  m_clusterPt.push_back(E_vec_cluster.perp());
207  m_clusterChisq.push_back(recoCluster -> get_chi2());
208  m_clusterEtIso.push_back(recoCluster -> get_et_iso(3,1,1));
209  }
210 
211 
212  //grab all the direct photons in the event
213  PHG4TruthInfoContainer::Range truthRange = truthinfo->GetPrimaryParticleRange();
214  PHG4TruthInfoContainer::ConstIterator truthIter;
215  //from the HepMC event log
216  for(truthIter = truthRange.first; truthIter != truthRange.second; truthIter++)
217  {
218  //PHG4Particle* part = truthinfo->GetParticle(truthIter->second->get_trkid())
219 
220  PHG4Particle *truthPar = truthIter->second;
221  int previousparentid = -999;
222  if(truthPar -> get_pid() != 22) continue;//end with a photon
223 
224 
225  int parentid = truthPar->get_parent_id();
226  previousparentid = truthPar->get_track_id();
227  while (parentid != 0)
228  {
229  previousparentid = parentid;
230  PHG4Particle* mommypart = truthinfo->GetParticle(parentid);
231  parentid = mommypart->get_parent_id();
232  }
233 
234  if(previousparentid > 0)truthPar = truthinfo -> GetParticle(previousparentid);
235  if(truthPar -> get_pid() != 22) continue;//start with a photon
236 
237  for(HepMC::GenEvent::particle_const_iterator p = theEvent -> particles_begin(); p != theEvent -> particles_end(); ++p)
238  {
239  assert(*p);
240 
241  if((*p) -> barcode() == truthPar -> get_barcode())
242  {
243  for(HepMC::GenVertex::particle_iterator mother = (*p)->production_vertex()->particles_begin(HepMC::parents);
244  mother != (*p)->production_vertex()->particles_end(HepMC::parents); ++mother)
245  {
246  //std::cout << "Mother pid: " << (*mother) -> pdg_id() << std::endl;
247  //if((*mother) -> pdg_id() != 21 || (*mother) -> pdg_id() != 1 || (*mother) -> pdg_id() != 2) continue;
248  if((*mother) -> pdg_id() != 22) continue;
249 
250  HepMC::FourVector moMomentum = (*mother) -> momentum();
251 
252  m_photonE.push_back(moMomentum.e());
253  m_photonEta.push_back(moMomentum.pseudoRapidity());
254  m_photonPhi.push_back(moMomentum.phi());
255  m_photonPt.push_back(sqrt(moMomentum.px()*moMomentum.px()+moMomentum.py()*moMomentum.py()));
256 
257 
258  }
259 
260  }
261  }
262  }
263 
264 
265  T -> Fill();
266 
267 return Fun4AllReturnCodes::EVENT_OK;
268 }
269 
270 //____________________________________________________________________________..
271 int isoCluster::ResetEvent(PHCompositeNode *topNode)
272 {
273 
274  m_clusterPhi.clear();
275  m_clusterEta.clear();
276  m_clusterE.clear();
277  m_clusterPt.clear();
278  m_clusterEtIso.clear();
279  m_clusterChisq.clear();
280  m_photonPhi.clear();
281  m_photonEta.clear();
282  m_photonE.clear();
283  m_photonPt.clear();
284 
285 
286  return Fun4AllReturnCodes::EVENT_OK;
287 }
288 
289 //____________________________________________________________________________..
290 int isoCluster::EndRun(const int runnumber)
291 {
292  std::cout << "isoCluster::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
293  return Fun4AllReturnCodes::EVENT_OK;
294 }
295 
296 //____________________________________________________________________________..
297 int isoCluster::End(PHCompositeNode *topNode)
298 {
299  std::cout << "isoCluster::End(PHCompositeNode *topNode) This is the End..." << std::endl;
300  fout -> cd();
301  T -> Write();
302  fout -> Close();
303  return Fun4AllReturnCodes::EVENT_OK;
304 }
305 
306 //____________________________________________________________________________..
307 int isoCluster::Reset(PHCompositeNode *topNode)
308 {
309  std::cout << "isoCluster::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
310  return Fun4AllReturnCodes::EVENT_OK;
311 }
312 
313 //____________________________________________________________________________..
314 void isoCluster::Print(const std::string &what) const
315 {
316  std::cout << "isoCluster::Print(const std::string &what) const Printing info for " << what << std::endl;
317 }