EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FarForwardEvaluator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FarForwardEvaluator.cc
1 #include "FarForwardEvaluator.h"
2 
3 #include <g4main/PHG4Particle.h>
4 #include <g4main/PHG4Shower.h>
5 #include <g4main/PHG4TruthInfoContainer.h>
6 #include <g4main/PHG4VtxPoint.h>
7 
8 #include <g4vertex/GlobalVertex.h>
9 #include <g4vertex/GlobalVertexMap.h>
10 
11 #include <calobase/RawCluster.h>
12 #include <calobase/RawClusterContainer.h>
13 #include <calobase/RawClusterUtility.h>
14 #include <calobase/RawTower.h>
15 #include <calobase/RawTowerContainer.h>
16 #include <calobase/RawTowerGeom.h>
17 #include <calobase/RawTowerGeomContainer.h>
18 
19 #include <fun4all/Fun4AllReturnCodes.h>
20 #include <fun4all/SubsysReco.h>
21 
22 #include <phool/getClass.h>
23 #include <phool/phool.h>
24 
25 #include <TFile.h>
26 #include <TNtuple.h>
27 #include <TTree.h>
28 
29 #include <CLHEP/Vector/ThreeVector.h>
30 
31 #include <cmath>
32 #include <cstdlib>
33 #include <iostream>
34 #include <set>
35 #include <utility>
36 
37 #include <sstream>
38 
39 #include <TSystem.h>
40 
41 // G4Hits includes
42 #include <g4main/PHG4Hit.h>
43 #include <g4main/PHG4HitContainer.h>
44 
45 FarForwardEvaluator::FarForwardEvaluator(const std::string& name, const std::string& ffrname, const std::string& filename, const std::string& ip_str)
46  : SubsysReco(name)
47  , _ffrname(ffrname)
48  , _ip_str(ip_str)
49  , g4hitntuple(nullptr)
50  , g4b0hitntuple(nullptr)
51  , g4rphitntuple(nullptr)
52  , h2_ZDC_XY(nullptr)
53  , h2_ZDC_XY_double(nullptr)
54  , h2_B0_XY(nullptr)
55  , h2_RP_XY(nullptr)
56  , hm(nullptr)
57  , _ievent(0)
58  , _towerID_debug(0)
59  , _ieta_debug(0)
60  , _iphi_debug(0)
61  , _eta_debug(0)
62  , _phi_debug(0)
63  , _e_debug(0)
64  , _x_debug(0)
65  , _y_debug(0)
66  , _z_debug(0)
67  , _truth_trace_embed_flags()
68  //, _truth_e_threshold(0.0)
69  //, 0 GeV before reco is traced
70  //,_reco_e_threshold(0.0)
71  //, 0 GeV before reco is traced
72  //,_caloevalstack(nullptr)
73  //, _strict(false)
74  , _do_gpoint_eval(true)
75  , _do_gshower_eval(true)
76  , _do_tower_eval(true)
77  , _do_cluster_eval(true)
78  , _ntp_gpoint(nullptr)
79  , _ntp_gshower(nullptr)
80  , _ntp_tower(nullptr)
81  , _tower_debug(nullptr)
82  , _ntp_cluster(nullptr)
83  , ZDC_hit(0)
84  , event_itt(0)
85  , _filename(filename)
86  , _tfile(nullptr)
87  , h1_E_dep_smeared(nullptr)
88  , h1_E_dep(nullptr)
89  , h1_B0_E_dep(nullptr)
90  , h1_B0_E(nullptr)
91  , h1_B0_E_abs(nullptr)
92  , h1_RP_E_dep(nullptr)
93  , h1_RP_E_abs(nullptr)
94 {
95 }
96 
97 int FarForwardEvaluator::Init(PHCompositeNode* topNode)
98 {
99  _ievent = 0;
100 
101  _tfile = new TFile(_filename.c_str(), "RECREATE");
102 
103  if (_do_gpoint_eval) _ntp_gpoint = new TNtuple("ntp_gpoint", "primary vertex => best (first) vertex",
104  "event:gvx:gvy:gvz:"
105  "vx:vy:vz");
106 
107  if (_do_gshower_eval) _ntp_gshower = new TNtuple("ntp_gshower", "truth shower => best cluster",
108  "event:gparticleID:gflavor:gnhits:"
109  "geta:gphi:ge:gpt:gvx:gvy:gvz:gembed:gedep:"
110  "clusterID:ntowers:eta:x:y:z:phi:e:efromtruth");
111 
112  //Barak: Added TTree to will allow the TowerID to be set correctly as integer
113  if (_do_tower_eval)
114  {
115  _ntp_tower = new TNtuple("ntp_tower", "tower => max truth primary",
116  "event:towerID:ieta:iphi:eta:phi:e:x:y:z:"
117  "gparticleID:gflavor:gnhits:"
118  "geta:gphi:ge:gpt:gvx:gvy:gvz:"
119  "gembed:gedep:"
120  "efromtruth");
121 
122  //Make Tree
123  _tower_debug = new TTree("tower_debug", "tower => max truth primary");
124 
125  _tower_debug->Branch("event", &_ievent, "event/I");
126  _tower_debug->Branch("towerID", &_towerID_debug, "towerID/I");
127  _tower_debug->Branch("ieta", &_ieta_debug, "ieta/I");
128  _tower_debug->Branch("iphi", &_iphi_debug, "iphi/I");
129  _tower_debug->Branch("eta", &_eta_debug, "eta/F");
130  _tower_debug->Branch("phi", &_phi_debug, "phi/F");
131  _tower_debug->Branch("e", &_e_debug, "e/F");
132  _tower_debug->Branch("x", &_x_debug, "x/F");
133  _tower_debug->Branch("y", &_y_debug, "y/F");
134  _tower_debug->Branch("z", &_z_debug, "z/F");
135  }
136 
137  if (_do_cluster_eval) _ntp_cluster = new TNtuple("ntp_cluster", "cluster => max truth primary",
138  "event:clusterID:ntowers:eta:x:y:z:phi:e:"
139  "gparticleID:gflavor:gnhits:"
140  "geta:gphi:ge:gpt:gvx:gvy:gvz:"
141  "gembed:gedep:"
142  "efromtruth");
143 
144  hm = new Fun4AllHistoManager(Name());
145 
146  g4hitntuple = new TNtuple("hitntup", "G4Hits", "x0:y0:z0:x1:y1:z1:edep");
147  g4b0hitntuple = new TNtuple("b0hit", "B0 Hits", "layer:type:x0:y0:z0:x1:y1:z1:t0:t1:edep");
148  g4rphitntuple = new TNtuple("rphit", "RP Hits", "layer:type:x0:y0:z0:x1:y1:z1:t0:t1:edep");
149 
150  std::cout << "diff_tagg_ana::Init(PHCompositeNode *topNode) Initializing" << std::endl;
151 
152  event_itt = 0;
153 
154  //----------------------------
155  // ZDC Occupancy
156 
157  gDirectory->mkdir("ZDC");
158  gDirectory->cd("ZDC");
159 
160  h2_ZDC_XY = new TH2F("ZDC_XY", "ZDC XY", 200, -50, 50, 200, -50, 50);
161 
162  h2_ZDC_XY_double = new TH2F("ZDC_XY_double", "ZDC XY Double gamma", 200, -50, 50, 200, -50, 50);
163 
164  h1_E_dep = new TH1F("E_dep", "E Dependence", 120, 0.0, 60.0);
165 
166  h1_E_dep_smeared = new TH1F("E_dep_smeared", "E Dependence Smeared", 120, 0.0, 60.0);
167 
168  gDirectory->cd("/");
169 
170  //----------------------------
171  // B0 Occupancy
172 
173  gDirectory->mkdir("B0");
174  gDirectory->cd("B0");
175 
176  h2_B0_XY = new TH2F("B0_XY", "B0 XY", 400, -200, 200, 200, -50, 50);
177 
178  h1_B0_E_dep = new TH1F("B0_E_dep", "E Dependence", 120, 0.0, 60.0);
179  h1_B0_E_abs = new TH1F("B0_E_abs", "E Absorbed", 120, 0.0, 60.0);
180  h1_B0_E = new TH1F("B0_E", "E layer", 10, 0.0, 10.0);
181 
182  gDirectory->cd("/");
183 
184  //----------------------------
185  // RP Occupancy
186 
187  gDirectory->mkdir("RP");
188  gDirectory->cd("RP");
189 
190  h2_RP_XY = new TH2F("RP_XY", "RP XY;X (cm);Y (cm)", 400, -200, 200, 200, -50, 50);
191  h1_RP_E_dep = new TH1F("RP_E_dep", "E Dependence", 120, 0.0, 60.0);
192  h1_RP_E_abs = new TH1F("RP_E_abs", "E Absorbed", 120, 0.0, 60.0);
193 
194  std::string paramFile = std::string(getenv("CALIBRATIONROOT")) + "/RomanPots/RP_parameters_IP6.dat";
195  if (_ip_str != "IP6")
196  {
197  paramFile = std::string(getenv("CALIBRATIONROOT")) + "/RomanPots/RP_parameters_IP8.dat";
198  }
199 
200  int Nlayers = GetParameterFromFile <int> (paramFile, "Number_layers");
201 
202  for( int layer = 0; layer < Nlayers; layer++ ) {
203  h2_RP_layers_XY.push_back( new TH2F(Form("RP_XY_layer%i",layer),"RP XY;X (cm);Y (cm)", 400, -200, 200, 200, -50, 50) );
204  h2_RP_virtlayers_XY.push_back( new TH2F(Form("RP_XY_virtlayer%i",layer),"RP XY;X (cm);Y (cm)", 4000, -200, 200, 200, -10, 10) );
205  }
206 
207  gDirectory->cd("/");
208 
209  return Fun4AllReturnCodes::EVENT_OK;
210 }
211 //
212 int FarForwardEvaluator::process_event(PHCompositeNode* topNode)
213 {
214  ZDC_hit = 0;
215 
216  event_itt++;
217 
218  if (event_itt % 100 == 0)
219  std::cout << "Event Processing Counter: " << event_itt << std::endl;
220 
221  process_g4hits_ZDC(topNode);
222 
223  process_g4hits_RomanPots(topNode);
224 
225  process_g4hits_B0(topNode);
226 
227  return Fun4AllReturnCodes::EVENT_OK;
228 }
229 
230 //***************************************************
231 
232 int FarForwardEvaluator::process_g4hits_ZDC(PHCompositeNode* topNode)
233 {
234  std::ostringstream nodename;
235 
236  nodename.str("");
237  nodename << "G4HIT_"
238  << "ZDCsurrogate";
239 
240  PHG4HitContainer* hits = findNode::getClass<PHG4HitContainer>(topNode, nodename.str().c_str());
241 
242  if (hits)
243  {
244  PHG4HitContainer::ConstRange hit_range = hits->getHits();
245  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
246  {
247  ZDC_hit++;
248  }
249 
250  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
251  {
252  float smeared_E = 0;
253  g4hitntuple->Fill(hit_iter->second->get_x(0),
254  hit_iter->second->get_y(0),
255  hit_iter->second->get_z(0),
256  hit_iter->second->get_x(1),
257  hit_iter->second->get_y(1),
258  hit_iter->second->get_z(1),
259  hit_iter->second->get_edep());
260 
261  float x_offset;
262 
263  if (_ip_str == "IP6")
264  {
265  x_offset = 90;
266  }
267  else
268  {
269  x_offset = -120;
270  }
271 
272  h2_ZDC_XY->Fill(hit_iter->second->get_x(0) + x_offset, hit_iter->second->get_y(0));
273 
274  //
275  // smeared_E = EMCAL_Smear(hit_iter->second->get_edep());
276  smeared_E = hit_iter->second->get_edep();
277  //
278  if (ZDC_hit == 2)
279  {
280  // cout << hit_iter->second->get_x(0)-90 << " " << hit_iter->second->get_y(0) << endl;
281  // h2_ZDC_XY_double->Fill(hit_iter->second->get_x(0)-90, hit_iter->second->get_y(0));
282  h2_ZDC_XY_double->Fill(hit_iter->second->get_x(0) + x_offset, hit_iter->second->get_y(0));
283  // h1_E_dep->Fill(hit_iter->second->get_edep());
284 
285  h1_E_dep->Fill(hit_iter->second->get_edep());
286  h1_E_dep_smeared->Fill(smeared_E);
287  //
288  }
289  }
290  }
291 
292  return Fun4AllReturnCodes::EVENT_OK;
293 }
294 
295 //***************************************************
296 // Getting the RomanPots hits
297 
298 int FarForwardEvaluator::process_g4hits_RomanPots(PHCompositeNode* topNode)
299 {
300  std::string nodename = "G4HIT_rpTruth";
301  PHG4HitContainer* hits = findNode::getClass<PHG4HitContainer>(topNode, nodename);
302 
303  std::string nodenameVirt = nodename + "_VirtSheet";
304  PHG4HitContainer* hits_virt = findNode::getClass<PHG4HitContainer>(topNode, nodenameVirt);
305 
306  if (hits)
307  {
308  // this returns an iterator to the beginning and the end of our G4Hits
309  PHG4HitContainer::ConstRange hit_range = hits->getHits();
310 
311  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
312  {
313  h2_RP_XY->Fill(hit_iter->second->get_x(0), hit_iter->second->get_y(0));
314 
315  if( h2_RP_layers_XY[ hit_iter->second->get_layer() ] ) {
316  h2_RP_layers_XY[ hit_iter->second->get_layer() ]->Fill(hit_iter->second->get_x(0), hit_iter->second->get_y(0));
317  }
318 
319  g4rphitntuple->Fill(hit_iter->second->get_layer(),
320  hit_iter->second->get_hit_type(),
321  hit_iter->second->get_x(0),
322  hit_iter->second->get_y(0),
323  hit_iter->second->get_z(0),
324  hit_iter->second->get_x(1),
325  hit_iter->second->get_y(1),
326  hit_iter->second->get_z(1),
327  hit_iter->second->get_t(0),
328  hit_iter->second->get_t(1),
329  hit_iter->second->get_edep());
330  if (hit_iter->second->get_hit_type())
331  h1_RP_E_dep->Fill(hit_iter->second->get_edep());
332  else
333  h1_RP_E_abs->Fill(hit_iter->second->get_edep());
334  }
335  }
336 
337  // hits on virtual layer without beam hole
338  if( hits_virt ) {
339 
340  PHG4HitContainer::ConstRange hit_range = hits_virt->getHits();
341 
342  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
343  {
344  if( h2_RP_virtlayers_XY[ hit_iter->second->get_layer() ] ) {
345  h2_RP_virtlayers_XY[ hit_iter->second->get_layer() ]->Fill(hit_iter->second->get_x(0), hit_iter->second->get_y(0));
346  }
347  }
348 
349  }
350 
351  return Fun4AllReturnCodes::EVENT_OK;
352 }
353 
354 //-----------------------------------
355 
356 //***************************************************
357 // Getting the B0 hits
358 
359 int FarForwardEvaluator::process_g4hits_B0(PHCompositeNode* topNode)
360 {
361  std::ostringstream nodename;
362 
363  nodename.str("");
364  nodename << "G4HIT_"
365  << "b0Truth";
366 
367  PHG4HitContainer* hits = findNode::getClass<PHG4HitContainer>(topNode, nodename.str().c_str());
368 
369  if (hits)
370  {
371  // // this returns an iterator to the beginning and the end of our G4Hits
372  PHG4HitContainer::ConstRange hit_range = hits->getHits();
373 
374  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
375  {
376  // cout << "B0 hits? " << endl;
377  // cout << "This is where you can fill your loop " << endl;
378 
379  h2_B0_XY->Fill(hit_iter->second->get_x(0), hit_iter->second->get_y(0));
380  g4b0hitntuple->Fill(hit_iter->second->get_layer(),
381  hit_iter->second->get_hit_type(),
382  hit_iter->second->get_x(0),
383  hit_iter->second->get_y(0),
384  hit_iter->second->get_z(0),
385  hit_iter->second->get_x(1),
386  hit_iter->second->get_y(1),
387  hit_iter->second->get_z(1),
388  hit_iter->second->get_t(0),
389  hit_iter->second->get_t(1),
390  hit_iter->second->get_edep());
391  if (hit_iter->second->get_hit_type())
392  h1_B0_E_dep->Fill(hit_iter->second->get_edep());
393  else
394  h1_B0_E_abs->Fill(hit_iter->second->get_edep());
395  h1_B0_E->Fill(hit_iter->second->get_layer(), hit_iter->second->get_edep());
396  }
397  }
398 
399  return Fun4AllReturnCodes::EVENT_OK;
400 }
401 
402 //***************************************************
403 
404 int FarForwardEvaluator::End(PHCompositeNode* topNode)
405 {
406  _tfile->cd();
407 
408  g4hitntuple->Write();
409  g4b0hitntuple->Write();
410  g4rphitntuple->Write();
411  _tfile->Write();
412  _tfile->Close();
413  delete _tfile;
414 
415  return Fun4AllReturnCodes::EVENT_OK;
416 }
417 
418 //--------------------------------------------------------
419 template <class T>
420 T GetParameterFromFile(std::string filename, std::string param)
421 {
422  std::ifstream infile;
423  std::string line;
424 
425  infile.open( filename );
426 
427  if( ! infile.is_open() )
428  {
429  std::cout << "ERROR in FarForwardEvaluator: Failed to open parameter file " << filename << std::endl;
430  gSystem->Exit(1);
431  }
432 
433  while( std::getline(infile, line) ) {
434 
435  std::string name;
436  double value;
437 
438  std::istringstream iss( line );
439 
440  // skip comment lines
441  if( line.find("#") != std::string::npos ) { continue; }
442 
443  if( !(iss >> name >> value) ) {
444  std::cout << "Could not decode " << line << std::endl;
445  gSystem->Exit(1);
446  }
447 
448  if( name.compare(param) == 0 ) {
449  return value;
450  }
451  }
452 
453  infile.close();
454  return 0;
455 }
456