Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HFMLTriggerOccupancy.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HFMLTriggerOccupancy.C
1 #include "HFMLTriggerOccupancy.h"
2 
3 #include <fun4all/Fun4AllReturnCodes.h>
4 #include <pdbcalbase/PdbParameterMap.h>
5 #include <phool/PHCompositeNode.h>
6 #include <phool/PHTimeServer.h>
7 #include <phool/PHTimer.h>
8 #include <phool/getClass.h>
9 
10 #include <phool/PHCompositeNode.h>
11 
12 #include <g4main/PHG4Hit.h>
13 #include <g4main/PHG4Particle.h>
14 #include <g4main/PHG4VtxPoint.h>
15 #include <g4main/PHG4TruthInfoContainer.h>
16 
17 #include <g4detectors/PHG4Cell.h>
18 #include <g4detectors/PHG4CellContainer.h>
19 #include <g4detectors/PHG4CylinderCellGeom.h>
20 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
21 #include <g4detectors/PHG4CylinderGeom.h>
22 #include <g4detectors/PHG4CylinderGeomContainer.h>
23 
24 #include <g4eval/SvtxEvalStack.h>
25 
26 #include <trackbase/TrkrHitSetContainerv1.h>
27 #include <trackbase/TrkrHitSetv1.h>
28 #include <trackbase/TrkrHitv1.h>
29 
30 #include <mvtx/CylinderGeom_Mvtx.h>
31 #include <mvtx/MvtxDefs.h>
32 
33 #include <HepMC/GenEvent.h>
34 #include <HepMC/GenVertex.h>
35 #include <phhepmc/PHHepMCGenEvent.h>
36 #include <phhepmc/PHHepMCGenEventMap.h>
37 
38 #include <TDatabasePDG.h>
39 #include <TFile.h>
40 #include <TH1D.h>
41 #include <TH1F.h>
42 #include <TH2F.h>
43 #include <TH3F.h>
44 #include <TLorentzVector.h>
45 #include <TString.h>
46 #include <TTree.h>
47 
48 #include <rapidjson/document.h>
49 #include <rapidjson/ostreamwrapper.h>
50 #include <rapidjson/prettywriter.h>
51 
52 #include <boost/format.hpp>
53 #include <boost/property_tree/json_parser.hpp>
54 #include <boost/property_tree/ptree.hpp>
55 #include <boost/property_tree/xml_parser.hpp>
56 
57 #include <cassert>
58 #include <cmath>
59 #include <cstddef>
60 #include <iostream>
61 
62 using namespace std;
63 
64 HFMLTriggerOccupancy::HFMLTriggerOccupancy(std::string filename)
65  : SubsysReco("HFMLTriggerOccupancy")
66  , _ievent(0)
67  , _f(nullptr)
68  , _eta_min(-1)
69  , _eta_max(1)
70  , _embedding_id(1)
71  , _nlayers_maps(3)
72  , _svtxevalstack(nullptr)
73  , m_hitsetcont(nullptr)
74  , m_Geoms(nullptr)
75  , m_truthInfo(nullptr)
76  , m_Flags(nullptr)
77  , m_hNorm(nullptr)
78  , m_hNChEta(nullptr)
79  , m_hVertexZ(nullptr)
80  , m_hitStaveLayer(nullptr)
81  , m_hitModuleHalfStave(nullptr)
82  , m_hitChipModule(nullptr)
83  , m_hitLayerMap(nullptr)
84  , m_hitPixelPhiMap(nullptr)
85  , m_hitPixelPhiMapHL(nullptr)
86  , m_hitPixelZMap(nullptr)
87  , m_Multiplicity(nullptr)
88  , m_LayerMultiplicity(nullptr)
89  , m_LayerChipMultiplicity(nullptr)
90 {
91  _foutname = filename;
92 }
93 
94 int HFMLTriggerOccupancy::Init(PHCompositeNode* topNode)
95 {
96  _ievent = 0;
97 
98  _f = new TFile((_foutname + string(".root")).c_str(), "RECREATE");
99 
100  // m_jsonOut.open((_foutname + string(".json")).c_str(), fstream::out);
101  //
102  // m_jsonOut << "{" << endl;
103  // m_jsonOut << "\"Events\" : [" << endl;
104 
105  // _h2 = new TH2D("h2", "", 100, 0, 100.0, 40, -2, +2);
106  // _h2_b = new TH2D("h2_b", "", 100, 0, 100.0, 40, -2, +2);
107  // _h2_c = new TH2D("h2_c", "", 100, 0, 100.0, 40, -2, +2);
108  // _h2all = new TH2D("h2all", "", 100, 0, 100.0, 40, -2, +2);
109 
110  m_hNorm = new TH1F("hNormalization", //
111  "Normalization;Items;Summed quantity", 10, .5, 10.5);
112  int i = 1;
113  m_hNorm->GetXaxis()->SetBinLabel(i++, "IntegratedLumi");
114  m_hNorm->GetXaxis()->SetBinLabel(i++, "Event");
115  m_hNorm->GetXaxis()->LabelsOption("v");
116 
117  m_hNChEta = new TH1F("hNChEta", //
118  "Charged particle #eta distribution;#eta;Count",
119  1000, -5, 5);
120  m_hVertexZ = new TH1F("hVertexZ", //
121  "Vertex z distribution;z [cm];Count",
122  1000, -200, 200);
123 
124  m_hitLayerMap = new TH3F("hitLayerMap", "hitLayerMap", 600, -10, 10, 600, -10, 10, 10, -.5, 9.5);
125  m_hitLayerMap->SetTitle("hitLayerMap;x [mm];y [mm];Half Layers");
126 
127 // m_hitPixelPhiMap = new TH3F("hitPixelPhiMap", "hitPixelPhiMap", 16000, -.5, 16000 - .5, 600, -M_PI, M_PI, 10, -.5, 9.5);
128 // m_hitPixelPhiMap->SetTitle("hitPixelPhiMap;PixelPhiIndex in layer;phi [rad];Half Layers Index");
129 // m_hitPixelPhiMapHL = new TH3F("hitPixelPhiMapHL", "hitPixelPhiMapHL", 16000, -.5, 16000 - .5, 600, -M_PI, M_PI, 10, -.5, 9.5);
130 // m_hitPixelPhiMapHL->SetTitle("hitPixelPhiMap;PixelPhiIndex in half layer;phi [rad];Half Layers Index");
131 // m_hitPixelZMap = new TH3F("hitPixelZMap", "hitPixelZMap", 16000, -.5, 16000 - .5, 600, 15, 15, 10, -.5, 9.5);
132 // m_hitPixelZMap->SetTitle("hitPixelZMap;hitPixelZMap;z [cm];Half Layers");
133 
134  m_hitStaveLayer = new TH2F("hitStaveLayer", "hitStaveLayer", 100, -.5, 100 - .5, 10, -.5, 9.5);
135  m_hitStaveLayer->SetTitle("hitStaveLayer;Stave index;Half Layers");
136  m_hitModuleHalfStave = new TH2F("hitModuleHalfStave", "hitModuleHalfStave", 100, -.5, 100 - .5, 10, -.5, 9.5);
137  m_hitModuleHalfStave->SetTitle("hitModuleHalfStave;Module index;Half Stave");
138  m_hitChipModule = new TH2F("hitChipModule", "hitChipModule", 100, -.5, 100 - .5, 10, -.5, 9.5);
139  m_hitChipModule->SetTitle("hitChipModule;Chip;Module");
140 
141  m_Multiplicity = new TH1F("Multiplicity", "Multiplicity", 10000, -.5, 10000 - .5);
142  m_Multiplicity->SetTitle("Multiplicity;Chip Multiplicity;Count");
143 
144  m_LayerMultiplicity = new TH2F("LayerMultiplicity", "LayerMultiplicity", 3, -.5, 2.5, 1000, -.5, 1000 - .5);
145  m_LayerMultiplicity->SetTitle("LayerMultiplicity;Layer;Chip Multiplicity");
146 
147  m_LayerChipMultiplicity = new TH3F("LayerChipMultiplicity", "LayerChipMultiplicity", 3, -.5, 2.5, kNCHip, -.5, kNCHip - .5, 1000, -.5, 1000 - .5);
148  m_LayerChipMultiplicity->SetTitle("LayerChipMultiplicity;Layer;Chip");
149 
150  return Fun4AllReturnCodes::EVENT_OK;
151 }
152 
153 int HFMLTriggerOccupancy::InitRun(PHCompositeNode* topNode)
154 {
155  m_hitsetcont = findNode::getClass<TrkrHitSetContainerv1>(topNode, "TRKR_HITSET");
156  if (!m_hitsetcont)
157  {
158  std::cout << PHWHERE << "ERROR: Can't find node TrkrHitSetContainer" << std::endl;
159  return Fun4AllReturnCodes::ABORTRUN;
160  }
161  m_Geoms =
162  findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
163  if (!m_Geoms)
164  {
165  std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_MVTX" << std::endl;
166  return Fun4AllReturnCodes::ABORTRUN;
167  }
168 
169  m_truthInfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
170  if (!m_truthInfo)
171  {
172  std::cout << PHWHERE << "ERROR: Can't find node G4TruthInfo" << std::endl;
173  return Fun4AllReturnCodes::ABORTRUN;
174  }
175 
176  m_Flags = findNode::getClass<PdbParameterMap>(topNode, "HFMLTrigger_HepMCTriggerFlags");
177  if (!m_Flags)
178  {
179  cout << "HFMLTriggerOccupancy::InitRun - WARNING - missing HFMLTrigger_HepMCTriggerFlags" << endl;
180  }
181 
182  return Fun4AllReturnCodes::EVENT_OK;
183 }
184 
185 int HFMLTriggerOccupancy::process_event(PHCompositeNode* topNode)
186 {
187  if (!_svtxevalstack)
188  {
189  _svtxevalstack = new SvtxEvalStack(topNode);
190  _svtxevalstack->set_strict(0);
191  _svtxevalstack->set_verbosity(Verbosity() + 1);
192  }
193  else
194  {
195  _svtxevalstack->next_event(topNode);
196  }
197 
198  // SvtxVertexEval* vertexeval = _svtxevalstack->get_vertex_eval();
199  // SvtxTrackEval* trackeval = _svtxevalstack->get_track_eval();
200  // SvtxClusterEval* clustereval = _svtxevalstack->get_cluster_eval();
201  // SvtxHitEval* hiteval = _svtxevalstack->get_hit_eval();
202  // SvtxTruthEval* trutheval = _svtxevalstack->get_truth_eval();
203 
204  // PHHepMCGenEventMap* geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
205  // if (!geneventmap)
206  // {
207  // std::cout << PHWHERE << " - Fatal error - missing node PHHepMCGenEventMap" << std::endl;
208  // return Fun4AllReturnCodes::ABORTRUN;
209  // }
210  //
211  // PHHepMCGenEvent* genevt = geneventmap->get(_embedding_id);
212  // if (!genevt)
213  // {
214  // std::cout << PHWHERE << " - Fatal error - node PHHepMCGenEventMap missing subevent with embedding ID " << _embedding_id;
215  // std::cout << ". Print PHHepMCGenEventMap:";
216  // geneventmap->identify();
217  // return Fun4AllReturnCodes::ABORTRUN;
218  // }
219  //
220  // HepMC::GenEvent* theEvent = genevt->getEvent();
221  // assert(theEvent);
222  // if (Verbosity())
223  // {
224  // cout << "HFMLTriggerOccupancy::process_event - process HepMC::GenEvent with signal_process_id = "
225  // << theEvent->signal_process_id();
226  // if (theEvent->signal_process_vertex())
227  // {
228  // cout << " and signal_process_vertex : ";
229  // theEvent->signal_process_vertex()->print();
230  // }
231  // cout << " - Event record:" << endl;
232  // theEvent->print();
233  // }
234 
235  assert(m_hNorm);
236  m_hNorm->Fill("Event", 1.);
237 
238  assert(m_truthInfo);
239  assert(m_hNChEta);
240  assert(m_hVertexZ);
241 
242  const PHG4VtxPoint* vtx =
243  m_truthInfo->GetPrimaryVtx(m_truthInfo->GetPrimaryVertexIndex());
244  if (vtx)
245  {
246  m_hVertexZ -> Fill(vtx->get_z());
247  }
248 
249  const auto primary_range =
250  m_truthInfo->GetPrimaryParticleRange();
251  for (auto particle_iter =
252  primary_range.first;
253  particle_iter != primary_range.second;
254  ++particle_iter)
255  {
256  const PHG4Particle *p = particle_iter->second;
257  assert(p);
258  TParticlePDG *pdg_p = TDatabasePDG::Instance()->GetParticle(
259  p->get_pid());
260  assert(pdg_p);
261  if (fabs(pdg_p->Charge()) > 0)
262  {
263  TVector3 pvec(p->get_px(), p->get_py(), p->get_pz());
264  if (pvec.Perp2() > 0)
265  {
266  assert(m_hNChEta);
267  m_hNChEta->Fill(pvec.PseudoRapidity());
268  }
269  }
270  } // if (_load_all_particle) else
271 
272 
273  // chip counting
274  vector<vector<vector<int> > >
275  multiplicity_vec(_nlayers_maps);
276 
277  //
278  for (unsigned int layer = 0; layer < _nlayers_maps; ++layer)
279  {
280  CylinderGeom_Mvtx* geom = dynamic_cast<CylinderGeom_Mvtx*>(m_Geoms->GetLayerGeom(layer));
281  assert(geom);
282 
283  const int n_stave = geom->get_N_staves();
284 
285  assert(multiplicity_vec[layer].empty());
286  multiplicity_vec[layer].resize(n_stave, vector<int>(kNCHip, 0));
287  }
288 
289  //set<unsigned int> mapsHits; //internal consistency check for later stages of truth tracks
290  assert(m_hitsetcont);
291 
292  TrkrHitSetContainer::ConstRange hitset_range = m_hitsetcont->getHitSets(TrkrDefs::TrkrId::mvtxId);
293  for (TrkrHitSetContainer::ConstIterator hitset_iter = hitset_range.first;
294  hitset_iter != hitset_range.second;
295  ++hitset_iter)
296  {
297  TrkrDefs::hitsetkey hitSetKey = hitset_iter->first;
298  TrkrHitSet::ConstRange hit_range = hitset_iter->second->getHits();
299  for (TrkrHitSet::ConstIterator hit_iter = hit_range.first;
300  hit_iter != hit_range.second;
301  ++hit_iter)
302  {
303  TrkrDefs::hitkey hitKey = hit_iter->first;
304  TrkrHit* hit = hit_iter->second;
305  assert(hit);
306  unsigned int layer = TrkrDefs::getLayer(hitset_iter->first);
307  if (layer < _nlayers_maps)
308  {
309  // unsigned int hitID = hit->get_id();
310  //mapsHits.insert(hitID);
311 
312  //PHG4Cell* cell = hiteval->get_cell(hit);
313  //assert(cell);
314 
315  CylinderGeom_Mvtx* geom = dynamic_cast<CylinderGeom_Mvtx*>(m_Geoms->GetLayerGeom(layer));
316  assert(geom);
317  unsigned int pixel_x = MvtxDefs::getRow(hitKey);
318  unsigned int pixel_z = MvtxDefs::getCol(hitKey);
319  assert((int)pixel_x < geom->get_NX());
320  assert((int)pixel_z < geom->get_NZ());
321  unsigned int chip = MvtxDefs::getChipId(hitSetKey);
322  unsigned int stave = MvtxDefs::getStaveId(hitSetKey);
323  TVector3 local_coords = geom->get_local_coords_from_pixel(pixel_x, pixel_z);
324  TVector3 world_coords = geom->get_world_from_local_coords(stave,
325  0,
326  0,
327  chip,
328  local_coords);
329 
330  //unsigned int pixel_x(cell->get_pixel_index() % geom->get_NX());
331  //unsigned int pixel_z(cell->get_pixel_index() / geom->get_NX());
332  //unsigned int halflayer = (int) pixel_x >= geom->get_NX() / 2 ? 0 : 1;
333 
334  //unsigned int halfLayerIndex(layer * 2 + halflayer);
335 // unsigned int pixelPhiIndex(
336 // cell->get_stave_index() * geom->get_NX() + pixel_x);
337 // unsigned int pixelPhiIndexHL(
338 // cell->get_stave_index() * geom->get_NX() / 2 + pixel_x % (geom->get_NX() / 2));
339 // unsigned int pixelZIndex(cell->get_chip_index() * geom->get_NZ() + pixel_z);
340 
341  // // ptree hitTree;
342  // rapidjson::Value hitTree(rapidjson::kObjectType);
343  //
344  // // ptree hitIDTree;
345  // rapidjson::Value hitIDTree(rapidjson::kObjectType);
346  // hitIDTree.AddMember("HitSequenceInEvent", hitID, alloc);
347  //
348  // hitIDTree.AddMember("PixelHalfLayerIndex", halfLayerIndex, alloc);
349  // hitIDTree.AddMember("PixelPhiIndexInLayer", pixelPhiIndex, alloc);
350  // hitIDTree.AddMember("PixelPhiIndexInHalfLayer", pixelPhiIndexHL, alloc);
351  // hitIDTree.AddMember("PixelZIndex", pixelZIndex, alloc);
352  //
353  // hitIDTree.AddMember("Layer", layer, alloc);
354  // hitIDTree.AddMember("HalfLayer", halflayer, alloc);
355  // hitIDTree.AddMember("Stave", cell->get_stave_index(), alloc);
356  // // hitIDTree.put("HalfStave", cell->get_half_stave_index());
357  // // hitIDTree.put("Module", cell->get_module_index());
358  // hitIDTree.AddMember("Chip", cell->get_chip_index(), alloc);
359  // hitIDTree.AddMember("Pixel", cell->get_pixel_index(), alloc);
360  // hitTree.AddMember("ID", hitIDTree, alloc);
361  //
362  // hitTree.AddMember("Coordinate",
363  // loadCoordinate(world_coords.x(),
364  // world_coords.y(),
365  // world_coords.z()),
366  // alloc);
367 
368  if (Verbosity() >= 2)
369  cout << "HFMLTriggerOccupancy::process_event - MVTX hit "
370  << "layer " << layer << "\t"
371  << "Stave " << stave << "\t"
372  << "Chip " << chip << "\t"
373  << "Pixel x " << pixel_x << "\t"
374  << "Pixel z " << pixel_z << "\t"
375  << "Coordinate"
376  << "(" << world_coords.x() << "," << world_coords.y() << "," << world_coords.z() << ")" << endl;
377 
378  // rawHitsTree.add_child("MVTXHit", hitTree);
379  // rawHitsTree.PushBack(hitTree, alloc);
380 
381  m_hitStaveLayer->Fill(stave, layer);
382  m_hitModuleHalfStave->Fill(stave, stave);
383  m_hitChipModule->Fill(chip, stave);
384 
385  m_hitLayerMap->Fill(world_coords.x(), world_coords.y(), layer);
386 // m_hitPixelPhiMap->Fill(pixelPhiIndex, atan2(world_coords.y(), world_coords.x()), halfLayerIndex);
387 // m_hitPixelPhiMapHL->Fill(pixelPhiIndexHL, atan2(world_coords.y(), world_coords.x()), halfLayerIndex);
388 // m_hitPixelZMap->Fill(pixelZIndex, world_coords.z(), halfLayerIndex);
389 
390  assert(layer < multiplicity_vec.size());
391  auto& multiplicity_layer = multiplicity_vec[layer];
392  assert(static_cast<size_t>(stave) < multiplicity_layer.size());
393  auto& multiplicity_stave = multiplicity_layer[stave];
394  assert(static_cast<size_t>(chip) < multiplicity_stave.size());
395  multiplicity_stave[chip]++;
396 
397  } // if (layer < _nlayers_maps)
398 
399  } // hit_iter loop
400 
401  } // hitset_iter loop
402  // rawHitTree.AddMember("MVTXHits", rawHitsTree, alloc);
403 
404  // m_Multiplicity = new TH1F("Multiplicity", "Multiplicity", 10000, -.5, 10000 - .5);
405  // m_Multiplicity->SetTitle("Multiplicity;Chip Multiplicity;Count");
406  //
407  // m_LayerMultiplicity = new TH2F("LayerMultiplicity", "LayerMultiplicity", 3, -.5, 2.5, 1000, -.5, 1000 - .5);
408  // m_LayerMultiplicity->SetTitle("LayerMultiplicity;Layer;Chip Multiplicity");
409  //
410  // m_LayerChipMultiplicity = new TH3F("LayerChipMultiplicity", "LayerChipMultiplicity", 3, -.5, 2.5, 9, -.5, 8.5, 1000, -.5, 1000 - .5);
411  // m_LayerChipMultiplicity->SetTitle("LayerChipMultiplicity;Layer;Chip");
412 
413  int layer = 0;
414  for (auto& multiplicity_layer : multiplicity_vec)
415  {
416  int stave = 0;
417  for (auto& multiplicity_stave : multiplicity_layer)
418  {
419  int chip = 0;
420  for (auto& multiplicity_chip : multiplicity_stave)
421  {
422  m_Multiplicity->Fill(multiplicity_chip);
423  m_LayerMultiplicity->Fill(layer, multiplicity_chip);
424  m_LayerChipMultiplicity->Fill(layer, chip, multiplicity_chip);
425 
426  chip++;
427  }
428 
429  if (Verbosity() >= 2)
430  cout << "HFMLTriggerOccupancy::process_event - fill layer " << layer << " stave. " << stave << "Accumulated chips = "
431  << m_Multiplicity->GetSum()
432  << endl;
433 
434  stave++;
435  }
436  layer++;
437  }
438 
439  ++_ievent;
440 
441  return Fun4AllReturnCodes::EVENT_OK;
442 }
443 
444 int HFMLTriggerOccupancy::End(PHCompositeNode* topNode)
445 {
446  if (_f)
447  {
448  _f->cd();
449  _f->Write();
450  //_f->Close();
451 
452  // m_hitLayerMap->Write();
453  }
454 
455  // if (m_jsonOut.is_open())
456  // {
457  // m_jsonOut << "]" << endl;
458  // m_jsonOut << "}" << endl;
459  //
460  // m_jsonOut.close();
461  // }
462 
463  cout << "HFMLTriggerOccupancy::End - output to " << _foutname << ".*" << endl;
464 
465  return Fun4AllReturnCodes::EVENT_OK;
466 }