Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HFMLTriggerHepMCTrigger.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HFMLTriggerHepMCTrigger.C
2 
6 #include <phool/PHTimeServer.h>
7 #include <phool/PHTimer.h>
8 #include <phool/getClass.h>
9 
10 #include <phool/PHCompositeNode.h>
11 
12 #include <pdbcalbase/PdbParameterMap.h>
13 
14 #include <g4main/PHG4Hit.h>
15 #include <g4main/PHG4Particle.h>
17 
18 
19 #include <g4detectors/PHG4Cell.h>
25 
26 //#include <trackbase_historic/SvtxCluster.h>
27 //#include <trackbase_historic/SvtxClusterMap.h>
28 //#include <trackbase_historic/SvtxHit.h>
29 //#include <trackbase_historic/SvtxHitMap.h>
30 //#include <trackbase_historic/SvtxTrack.h>
31 //#include <trackbase_historic/SvtxTrackMap.h>
32 //#include <trackbase_historic/SvtxVertex.h>
33 //#include <trackbase_historic/SvtxVertexMap.h>
34 
35 #include <g4eval/SvtxEvalStack.h>
36 
37 #include <HepMC/GenEvent.h>
38 #include <HepMC/GenRanges.h>
39 #include <HepMC/GenVertex.h>
40 #include <ffaobjects/FlagSavev1.h>
43 
44 #include <TDatabasePDG.h>
45 #include <TFile.h>
46 #include <TH2F.h>
47 #include <TH3F.h>
48 #include <TLorentzVector.h>
49 #include <TTree.h>
50 
51 #include <rapidjson/document.h>
52 #include <rapidjson/ostreamwrapper.h>
53 #include <rapidjson/prettywriter.h>
54 
55 #include <boost/format.hpp>
56 #include <boost/property_tree/json_parser.hpp>
57 #include <boost/property_tree/ptree.hpp>
58 #include <boost/property_tree/xml_parser.hpp>
59 
60 #include <cassert>
61 #include <cmath>
62 #include <cstddef>
63 #include <iostream>
64 
65 using namespace std;
66 
68  const std::string& filename)
69  : SubsysReco(moduleName)
70  , _ievent(0)
71  , m_RejectReturnCode(Fun4AllReturnCodes::ABORTEVENT)
72  , _f(nullptr)
73  , _eta_min(-1)
74  , _eta_max(1)
75  , _embedding_id(1)
76  , m_Geneventmap(nullptr)
77  , m_Flags(nullptr)
78  , m_hNorm(nullptr)
79  , m_DRapidity(nullptr)
80 {
82 }
83 
85 {
86  _ievent = 0;
87 
88  _f = new TFile((_foutname + string(".root")).c_str(), "RECREATE");
89 
90  m_hNorm = new TH1D("hNormalization", //
91  "Normalization;Items;Summed quantity", 10, .5, 10.5);
92  int i = 1;
93  m_hNorm->GetXaxis()->SetBinLabel(i++, "IntegratedLumi");
94  m_hNorm->GetXaxis()->SetBinLabel(i++, "Event");
95  m_hNorm->GetXaxis()->SetBinLabel(i++, "D0");
96  m_hNorm->GetXaxis()->SetBinLabel(i++, "D0->PiK");
97  m_hNorm->GetXaxis()->SetBinLabel(i++, "D0-Pair");
98  m_hNorm->GetXaxis()->SetBinLabel(i++, "D0->PiK-Pair");
99  m_hNorm->GetXaxis()->SetBinLabel(i++, "Accepted");
100 
101  m_hNorm->GetXaxis()->LabelsOption("v");
102 
103  m_DRapidity = new TH2F("hDRapidity", //
104  "hDRapidity;Rapidity of D0 meson;Accepted", 1000, -5, 5, 2, -.5, 1.5);
105 
107 }
108 
110 {
111  m_Geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
112  if (!m_Geneventmap)
113  {
114  std::cout << PHWHERE << " - Fatal error - missing node PHHepMCGenEventMap" << std::endl;
116  }
117 
118  PHNodeIterator iter(topNode);
119 
120  PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
121  if (!dstNode)
122  {
123  cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
124  throw std::runtime_error(
125  "Failed to find DST node in RawTowerBuilder::CreateNodes");
126  }
127 
128  m_Flags = findNode::getClass<PdbParameterMap>(dstNode, "HFMLTrigger_HepMCTriggerFlags");
129  if (m_Flags == nullptr)
130  {
131  m_Flags = new PdbParameterMap();
132 
133  dstNode->addNode(new PHDataNode<PHObject>(m_Flags, "HFMLTrigger_HepMCTriggerFlags", "PHObject"));
134  }
135 
137 }
138 
140 {
142 
144  if (!genevt)
145  {
146  std::cout << PHWHERE << " - Fatal error - node PHHepMCGenEventMap missing subevent with embedding ID " << _embedding_id;
147  std::cout << ". Print PHHepMCGenEventMap:";
150  }
151 
152  HepMC::GenEvent* theEvent = genevt->getEvent();
153  assert(theEvent);
154  if (Verbosity() >= VERBOSITY_MORE)
155  {
156  cout << "HFMLTriggerHepMCTrigger::process_event - process HepMC::GenEvent with signal_process_id = "
157  << theEvent->signal_process_id();
158  if (theEvent->signal_process_vertex())
159  {
160  cout << " and signal_process_vertex : ";
161  theEvent->signal_process_vertex()->print();
162  }
163  cout << " - Event record:" << endl;
164  theEvent->print();
165  }
166 
167  TDatabasePDG* pdg = TDatabasePDG::Instance();
168 
169  int targetPID = std::abs(pdg->GetParticle("D0")->PdgCode());
170  int daughter1PID = std::abs(pdg->GetParticle("pi+")->PdgCode());
171  int daughter2PID = std::abs(pdg->GetParticle("K+")->PdgCode());
172 
173  bool acceptEvent = false;
174 
175  assert(m_hNorm);
176  m_hNorm->Fill("Event", 1);
177 
178  unsigned int nD0(0);
179  unsigned int nD0PiK(0);
180 
181  auto range = theEvent->particle_range();
182  for (HepMC::GenEvent::particle_const_iterator piter = range.begin(); piter != range.end(); ++piter)
183  {
184  const HepMC::GenParticle* p = *piter;
185  assert(p);
186 
187  if (std::abs(p->pdg_id()) == targetPID)
188  {
189  if (Verbosity())
190  {
191  cout << "HFMLTriggerHepMCTrigger::process_event - Accept signal particle : ";
192  p->print();
193  cout << endl;
194  }
195 
196  m_hNorm->Fill("D0", 1);
197  ++nD0;
198 
200  const double rapidity = 0.5 * log((p->momentum().e() + p->momentum().z()) /
201  (p->momentum().e() - p->momentum().z()));
202 
203  m_DRapidity->Fill(rapidity, 0);
204 
205  const HepMC::GenVertex* decayVertex = p->end_vertex();
206 
207  int hasDecay1(0);
208  int hasDecay2(0);
209  int hasDecayOther(0);
210 
211  if (decayVertex)
212  {
213  for (auto diter = decayVertex->particles_out_const_begin();
214  diter != decayVertex->particles_out_const_end();
215  ++diter)
216 
217  {
218  const HepMC::GenParticle* pd = *diter;
219  assert(pd);
220 
221  if (Verbosity())
222  {
223  cout << "HFMLTriggerHepMCTrigger::process_event - Testing daughter particle: ";
224  pd->print();
225  cout << endl;
226  }
227 
228  if (std::abs(pd->pdg_id()) == daughter1PID)
229  {
230  const double eta = pd->momentum().eta();
231 
232  if (eta > _eta_min and eta < _eta_max)
233  ++hasDecay1;
234  }
235  else if (std::abs(pd->pdg_id()) == daughter2PID)
236  {
237  const double eta = pd->momentum().eta();
238 
239  if (eta > _eta_min and eta < _eta_max)
240  ++hasDecay2;
241  }
242  else
243  ++hasDecayOther;
244  }
245 
246  if (hasDecay1 == 1 and hasDecay2 == 1 and hasDecayOther == 0)
247  {
248  m_hNorm->Fill("D0->PiK", 1);
249  ++nD0PiK;
250 
251  acceptEvent = true;
252 
253  m_DRapidity->Fill(rapidity, 1);
254  }
255 
256  } // if (decayVertex)
257  else
258  {
259  cout << "HFMLTriggerHepMCTrigger::process_event - Warning - target particle did not have decay vertex : ";
260  p->print();
261  cout << endl;
262  }
263 
264  } // if (std::abs(p-> pdg_id()) == targetPID)
265  } // for (HepMC::GenEvent::particle_const_iterator piter = range.begin(); piter != range.end(); ++piter)
266 
267  if (nD0 >= 2)
268  {
269  cout <<"HFMLTriggerHepMCTrigger::process_event - D0-Pair with nD0 = "<<nD0<<endl;
270  m_hNorm->Fill("D0-Pair", nD0 * (nD0 - 1) / 2);
271  }
272  if (nD0PiK >= 2)
273  {
274  m_hNorm->Fill("D0->PiK-Pair", nD0PiK * (nD0PiK - 1) / 2);
275  }
276 
277  ++_ievent;
278 
279  if (Verbosity())
280  {
281  cout << "HFMLTriggerHepMCTrigger::process_event - acceptEvent = " << acceptEvent;
282  cout << endl;
283 
284  if (acceptEvent)
285  {
286  cout << "HFMLTriggerHepMCTrigger::process_event - processed HepMC::GenEvent with signal_process_id = "
287  << theEvent->signal_process_id();
288  if (theEvent->signal_process_vertex())
289  {
290  cout << " and signal_process_vertex : ";
291  theEvent->signal_process_vertex()->print();
292  }
293  cout << " - Event record:" << endl;
294  theEvent->print();
295  }
296  }
297 
298  assert(m_Flags);
299  m_Flags->set_int_param(Name(), acceptEvent);
300 
301  if (acceptEvent)
302  {
303  m_hNorm->Fill("Accepted", 1);
305  }
306  else
307  return m_RejectReturnCode;
308 }
309 
311 {
312  PHGenIntegral* integral_node = findNode::getClass<PHGenIntegral>(topNode, "PHGenIntegral");
313  if (integral_node)
314  {
315  assert(m_hNorm);
316  m_hNorm->Fill("IntegratedLumi", integral_node->get_Integrated_Lumi());
317  }
318 
319  if (_f)
320  {
321  _f->cd();
322  _f->Write();
323  //_f->Close();
324 
325  // m_hitLayerMap->Write();
326  }
327 
328  cout << "HFMLTriggerHepMCTrigger::End - output to " << _foutname << ".*" << endl;
329 
331 }