Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4_CEmc_Spacal.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4_CEmc_Spacal.C
1 #ifndef MACRO_G4CEMCSPACAL_C
2 #define MACRO_G4CEMCSPACAL_C
3 
4 #include <GlobalVariables.C>
5 #include <QA.C>
6 
12 
13 #include <g4calo/RawTowerBuilder.h>
15 
16 #include <g4eval/CaloEvaluator.h>
17 
18 #include <g4main/PHG4Reco.h>
19 #include <g4main/PHG4Utils.h>
20 
21 #include <caloreco/RawClusterBuilderGraph.h>
22 #include <caloreco/RawClusterBuilderTemplate.h>
23 #include <caloreco/RawClusterPositionCorrection.h>
24 #include <caloreco/RawTowerCalibration.h>
25 #include <qa_modules/QAG4SimulationCalorimeter.h>
26 
27 #include <fun4all/Fun4AllServer.h>
28 
29 double
30 CEmc_1DProjectiveSpacal(PHG4Reco *g4Reco, double radius, const int crossings);
31 
32 double
33 CEmc_2DProjectiveSpacal(PHG4Reco *g4Reco, double radius, const int crossings);
34 
35 R__LOAD_LIBRARY(libcalo_reco.so)
36 R__LOAD_LIBRARY(libg4calo.so)
37 R__LOAD_LIBRARY(libg4detectors.so)
38 R__LOAD_LIBRARY(libg4eval.so)
39 R__LOAD_LIBRARY(libqa_modules.so)
40 
41 namespace Enable
42 {
43  bool CEMC = false;
44  bool CEMC_ABSORBER = false;
45  bool CEMC_OVERLAPCHECK = false;
46  bool CEMC_CELL = false;
47  bool CEMC_TOWER = false;
48  bool CEMC_CLUSTER = false;
49  bool CEMC_EVAL = false;
50  bool CEMC_QA = false;
51  int CEMC_VERBOSITY = 0;
52 } // namespace Enable
53 
54 namespace G4CEMC
55 {
56  int Min_cemc_layer = 1;
57  int Max_cemc_layer = 1;
58 
59  // Digitization (default photon digi):
61  // directly pass the energy of sim tower to digitized tower
62  // kNo_digitization
63  // simple digitization with photon statistics, single amplitude ADC conversion and pedestal
64  // kSimple_photon_digitization
65  // digitization with photon statistics on SiPM with an effective pixel N, ADC conversion and pedestal
66  // kSiPM_photon_digitization
67 
68  // set a default value for SPACAL configuration
69  // // 1D azimuthal projective SPACAL (fast)
70  //int Cemc_spacal_configuration = PHG4CylinderGeom_Spacalv1::k1DProjectiveSpacal;
71  // 2D azimuthal projective SPACAL (slow)
73 
75  {
77 
79  };
80 
84  //enu_Cemc_clusterizer Cemc_clusterizer = kCemcGraphClusterizer;
85 
86 } // namespace G4CEMC
87 
88 // black hole parameters are set in CEmc function
89 // needs a dummy argument to play with current G4Setup_sPHENIX.C
90 void CEmcInit(const int i = 0)
91 {
92 }
93 
95 double
96 CEmc(PHG4Reco *g4Reco, double radius, const int crossings)
97 {
99  {
100  return CEmc_1DProjectiveSpacal(/*PHG4Reco**/ g4Reco, /*double*/ radius, /*const int */ crossings);
101  }
103  {
104  return CEmc_2DProjectiveSpacal(/*PHG4Reco**/ g4Reco, /*double*/ radius, /*const int */ crossings);
105  }
106  else
107  {
108  std::cout
109  << "G4_CEmc_Spacal.C::CEmc - Fatal Error - unrecognized SPACAL configuration #"
110  << G4CEMC::Cemc_spacal_configuration << ". Force exiting..." << std::endl;
111  exit(-1);
112  return 0;
113  }
114 }
115 
117 double
118 CEmc_1DProjectiveSpacal(PHG4Reco *g4Reco, double radius, const int crossings)
119 {
120  bool AbsorberActive = Enable::ABSORBER || Enable::CEMC_ABSORBER;
121  bool OverlapCheck = Enable::OVERLAPCHECK || Enable::CEMC_OVERLAPCHECK;
122 
123  double emc_inner_radius = 95.; // emc inner radius from engineering drawing
124  double cemcthickness = 12.7;
125  double emc_outer_radius = emc_inner_radius + cemcthickness; // outer radius
126 
127  if (radius > emc_inner_radius)
128  {
129  cout << "inconsistency: pstof outer radius: " << radius
130  << " larger than emc inner radius: " << emc_inner_radius
131  << endl;
132  gSystem->Exit(-1);
133  }
134 
135  // boundary check
136  if (radius > emc_inner_radius - 1.5 - no_overlapp)
137  {
138  cout << "G4_CEmc_Spacal.C::CEmc() - expect radius < " << emc_inner_radius - 1.5 - no_overlapp << " to install SPACAL" << endl;
139  exit(1);
140  }
141  radius = emc_inner_radius - 1.5 - no_overlapp;
142 
143  // 1.5cm thick teflon as an approximation for EMCAl light collection + electronics (10% X0 total estimated)
144  PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("CEMC_ELECTRONICS", 0);
145  cyl->SuperDetector("CEMC_ELECTRONICS");
146  cyl->set_double_param("radius", radius);
147  cyl->set_string_param("material", "G4_TEFLON");
148  cyl->set_double_param("thickness", 1.5);
149  if (AbsorberActive) cyl->SetActive();
150  g4Reco->registerSubsystem(cyl);
151 
152  radius += 1.5;
153  radius += no_overlapp;
154 
155  int ilayer = G4CEMC::Min_cemc_layer;
156  PHG4SpacalSubsystem *cemc = new PHG4SpacalSubsystem("CEMC", ilayer);
157  cemc->set_double_param("radius", emc_inner_radius);
158  cemc->set_double_param("thickness", cemcthickness);
159 
160  cemc->SetActive();
161  cemc->SuperDetector("CEMC");
162  if (AbsorberActive) cemc->SetAbsorberActive();
163  cemc->OverlapCheck(OverlapCheck);
164 
165  g4Reco->registerSubsystem(cemc);
166 
167  if (ilayer > G4CEMC::Max_cemc_layer)
168  {
169  cout << "layer discrepancy, current layer " << ilayer
170  << " max cemc layer: " << G4CEMC::Max_cemc_layer << endl;
171  }
172 
173  radius += cemcthickness;
174  radius += no_overlapp;
175 
176  // 0.5cm thick Stainless Steel as an approximation for EMCAl support system
177  cyl = new PHG4CylinderSubsystem("CEMC_SPT", 0);
178  cyl->SuperDetector("CEMC_SPT");
179  cyl->set_double_param("radius", radius);
180  cyl->set_string_param("material", "SS310"); // SS310 Stainless Steel
181  cyl->set_double_param("thickness", 0.5);
182  if (AbsorberActive) cyl->SetActive();
183  g4Reco->registerSubsystem(cyl);
184  radius += 0.5;
185  // this is the z extend and outer radius of the support structure and therefore the z extend
186  // and radius of the surrounding black holes
190  radius += no_overlapp;
191 
192  return radius;
193 }
194 
196 double
197 CEmc_2DProjectiveSpacal(PHG4Reco *g4Reco, double radius, const int crossings)
198 {
199  bool AbsorberActive = Enable::ABSORBER || Enable::CEMC_ABSORBER;
200  bool OverlapCheck = Enable::OVERLAPCHECK || Enable::CEMC_OVERLAPCHECK;
201 
202  double emc_inner_radius = 92; // emc inner radius from engineering drawing
203  double cemcthickness = 24.00000 - no_overlapp;
204 
205  //max radius is 116 cm;
206  double emc_outer_radius = emc_inner_radius + cemcthickness; // outer radius
207  assert(emc_outer_radius < 116);
208 
209  if (radius > emc_inner_radius)
210  {
211  cout << "inconsistency: preshower radius+thickness: " << radius
212  << " larger than emc inner radius: " << emc_inner_radius << endl;
213  gSystem->Exit(-1);
214  }
215 
216  // the radii are only to determined the thickness of the cemc
217  radius = emc_inner_radius;
218 
219  // 1.5cm thick teflon as an approximation for EMCAl light collection + electronics (10% X0 total estimated)
220  PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("CEMC_ELECTRONICS", 0);
221  cyl->set_double_param("radius", radius);
222  cyl->set_string_param("material", "G4_TEFLON");
223  cyl->set_double_param("thickness", 1.5 - no_overlapp);
224  cyl->SuperDetector("CEMC_ELECTRONICS");
225  cyl->OverlapCheck(OverlapCheck);
226  if (AbsorberActive) cyl->SetActive();
227  g4Reco->registerSubsystem(cyl);
228 
229  radius += 1.5;
230  cemcthickness -= 1.5 + no_overlapp;
231 
232  // 0.5cm thick Stainless Steel as an approximation for EMCAl support system
233  cyl = new PHG4CylinderSubsystem("CEMC_SPT", 0);
234  cyl->SuperDetector("CEMC_SPT");
235  cyl->set_double_param("radius", radius + cemcthickness - 0.5);
236  cyl->set_string_param("material", "SS310"); // SS310 Stainless Steel
237  cyl->set_double_param("thickness", 0.5 - no_overlapp);
238  cyl->OverlapCheck(OverlapCheck);
239  if (AbsorberActive) cyl->SetActive();
240  g4Reco->registerSubsystem(cyl);
241 
242  // this is the z extend and outer radius of the support structure and therefore the z extend
243  // and radius of the surrounding black holes
244  double sptlen = PHG4Utils::GetLengthForRapidityCoverage(radius + cemcthickness);
247  BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, radius + cemcthickness);
248 
249  cemcthickness -= 0.5 + no_overlapp;
250 
251  int ilayer = 0;
253 
254  cemc = new PHG4SpacalSubsystem("CEMC", ilayer);
255 
256  cemc->set_int_param("virualize_fiber", 0);
257  cemc->set_int_param("azimuthal_seg_visible", 1);
258  cemc->set_int_param("construction_verbose", 0);
259  cemc->Verbosity(0);
260 
262  cemc->SetCalibrationFileDir(string(getenv("CALIBRATIONROOT")) + string("/CEMC/Geometry_2018ProjTilted/"));
263  cemc->set_double_param("radius", radius); // overwrite minimal radius
264  cemc->set_double_param("thickness", cemcthickness); // overwrite thickness
265 
266  cemc->SetActive();
267  cemc->SuperDetector("CEMC");
268  if (AbsorberActive) cemc->SetAbsorberActive();
269  cemc->OverlapCheck(OverlapCheck);
270 
271  g4Reco->registerSubsystem(cemc);
272 
273  if (ilayer > G4CEMC::Max_cemc_layer)
274  {
275  cout << "layer discrepancy, current layer " << ilayer
276  << " max cemc layer: " << G4CEMC::Max_cemc_layer << endl;
277  }
278 
279  radius += cemcthickness;
280  radius += no_overlapp;
281 
282  return radius;
283 }
284 
286 {
288 
290 
292  {
293  PHG4CylinderCellReco *cemc_cells = new PHG4CylinderCellReco("CEMCCYLCELLRECO");
294  cemc_cells->Detector("CEMC");
295  cemc_cells->Verbosity(verbosity);
297  {
298  // cemc_cells->etaphisize(i, 0.024, 0.024);
299  const double radius = 95;
300  cemc_cells->cellsize(i, 2 * M_PI / 256. * radius, 2 * M_PI / 256. * radius);
301  }
302  se->registerSubsystem(cemc_cells);
303  }
305  {
306  PHG4FullProjSpacalCellReco *cemc_cells = new PHG4FullProjSpacalCellReco("CEMCCYLCELLRECO");
307  cemc_cells->Detector("CEMC");
308  cemc_cells->Verbosity(verbosity);
310  string(getenv("CALIBRATIONROOT")) + string("/CEMC/LightCollection/Prototype3Module.xml"),
311  "data_grid_light_guide_efficiency", "data_grid_fiber_trans");
312  se->registerSubsystem(cemc_cells);
313  }
314  else
315  {
316  cout << "G4_CEmc_Spacal.C::CEmc - Fatal Error - unrecognized SPACAL configuration #"
317  << G4CEMC::Cemc_spacal_configuration << ". Force exiting..." << endl;
318  gSystem->Exit(-1);
319  return;
320  }
321 
322  return;
323 }
324 
326 {
328 
330 
331  RawTowerBuilder *TowerBuilder = new RawTowerBuilder("EmcRawTowerBuilder");
332  TowerBuilder->Detector("CEMC");
333  TowerBuilder->set_sim_tower_node_prefix("SIM");
334  TowerBuilder->Verbosity(verbosity);
335  se->registerSubsystem(TowerBuilder);
336 
337  double sampling_fraction = 1;
339  {
340  sampling_fraction = 0.0234335; //from production:/gpfs02/phenix/prod/sPHENIX/preCDR/pro.1-beta.3/single_particle/spacal1d/zerofield/G4Hits_sPHENIX_e-_eta0_8GeV.root
341  }
343  {
344  // sampling_fraction = 0.02244; //from production: /gpfs02/phenix/prod/sPHENIX/preCDR/pro.1-beta.3/single_particle/spacal2d/zerofield/G4Hits_sPHENIX_e-_eta0_8GeV.root
345  // sampling_fraction = 2.36081e-02; //from production: /gpfs02/phenix/prod/sPHENIX/preCDR/pro.1-beta.5/single_particle/spacal2d/zerofield/G4Hits_sPHENIX_e-_eta0_8GeV.root
346  // sampling_fraction = 1.90951e-02; // 2017 Tilt porjective SPACAL, 8 GeV photon, eta = 0.3 - 0.4
347  sampling_fraction = 2e-02; // 2017 Tilt porjective SPACAL, tower-by-tower calibration
348  }
349  else
350  {
351  std::cout
352  << "G4_CEmc_Spacal.C::CEMC_Towers - Fatal Error - unrecognized SPACAL configuration #"
353  << G4CEMC::Cemc_spacal_configuration << ". Force exiting..." << std::endl;
354  exit(-1);
355  return;
356  }
357 
358  const double photoelectron_per_GeV = 500; //500 photon per total GeV deposition
359 
360  RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("EmcRawTowerDigitizer");
361  TowerDigitizer->Detector("CEMC");
362  TowerDigitizer->Verbosity(verbosity);
363  TowerDigitizer->set_digi_algorithm(G4CEMC::TowerDigi);
364  TowerDigitizer->set_variable_pedestal(true); //read ped central and width from calibrations file comment next 2 lines if true
365  // TowerDigitizer->set_pedstal_central_ADC(0);
366  // TowerDigitizer->set_pedstal_width_ADC(8); // eRD1 test beam setting
367  TowerDigitizer->set_photonelec_ADC(1); //not simulating ADC discretization error
368  TowerDigitizer->set_photonelec_yield_visible_GeV(photoelectron_per_GeV / sampling_fraction);
369  TowerDigitizer->set_variable_zero_suppression(true); //read zs values from calibrations file comment next line if true
370  // TowerDigitizer->set_zero_suppression_ADC(16); // eRD1 test beam setting
371  TowerDigitizer->GetParameters().ReadFromFile("CEMC", "xml", 0, 0,
372  string(getenv("CALIBRATIONROOT")) + string("/CEMC/TowerCalibCombinedParams_2020/")); // calibration database
373  se->registerSubsystem(TowerDigitizer);
374 
375  RawTowerCalibration *TowerCalibration = new RawTowerCalibration("EmcRawTowerCalibration");
376  TowerCalibration->Detector("CEMC");
377  TowerCalibration->Verbosity(verbosity);
378 
380  {
382  {
383  // just use sampling fraction set previously
384  TowerCalibration->set_calib_const_GeV_ADC(1.0 / sampling_fraction);
385  }
386  else
387  {
389  TowerCalibration->set_calib_const_GeV_ADC(1. / photoelectron_per_GeV);
390  TowerCalibration->set_pedstal_ADC(0);
391  }
392  }
394  {
396  {
397  // just use sampling fraction set previously
398  TowerCalibration->set_calib_const_GeV_ADC(1.0 / sampling_fraction);
399  }
400  else
401  {
403  TowerCalibration->GetCalibrationParameters().ReadFromFile("CEMC", "xml", 0, 0,
404  string(getenv("CALIBRATIONROOT")) + string("/CEMC/TowerCalibCombinedParams_2020/")); // calibration database
405  TowerCalibration->set_variable_GeV_ADC(true); //read GeV per ADC from calibrations file comment next line if true
406  // TowerCalibration->set_calib_const_GeV_ADC(1. / photoelectron_per_GeV / 0.9715); // overall energy scale based on 4-GeV photon simulations
407  TowerCalibration->set_variable_pedestal(true); //read pedestals from calibrations file comment next line if true
408  // TowerCalibration->set_pedstal_ADC(0);
409  }
410  }
411  else
412  {
413  cout << "G4_CEmc_Spacal.C::CEMC_Towers - Fatal Error - unrecognized SPACAL configuration #"
414  << G4CEMC::Cemc_spacal_configuration << ". Force exiting..." << endl;
415  gSystem->Exit(-1);
416  return;
417  }
418  se->registerSubsystem(TowerCalibration);
419 
420  return;
421 }
422 
424 {
426 
428 
430  {
431  RawClusterBuilderTemplate *ClusterBuilder = new RawClusterBuilderTemplate("EmcRawClusterBuilderTemplate");
432  ClusterBuilder->Detector("CEMC");
433  ClusterBuilder->Verbosity(verbosity);
434  ClusterBuilder->set_threshold_energy(0.030); // This threshold should be the same as in CEMCprof_Thresh**.root file below
435  std::string emc_prof = getenv("CALIBRATIONROOT");
436  emc_prof += "/EmcProfile/CEMCprof_Thresh30MeV.root";
437  ClusterBuilder->LoadProfile(emc_prof);
438  se->registerSubsystem(ClusterBuilder);
439  }
441  {
442  RawClusterBuilderGraph *ClusterBuilder = new RawClusterBuilderGraph("EmcRawClusterBuilderGraph");
443  ClusterBuilder->Detector("CEMC");
444  ClusterBuilder->Verbosity(verbosity);
445  se->registerSubsystem(ClusterBuilder);
446  }
447  else
448  {
449  cout << "CEMC_Clusters - unknown clusterizer setting!" << endl;
450  exit(1);
451  }
452 
453  RawClusterPositionCorrection *clusterCorrection = new RawClusterPositionCorrection("CEMC");
454 
455  clusterCorrection->Get_eclus_CalibrationParameters().ReadFromFile("CEMC_RECALIB", "xml", 0, 0,
456  //raw location
457  string(getenv("CALIBRATIONROOT")) + string("/CEMC/PositionRecalibration_EMCal_9deg_tilt/"));
458 
459  clusterCorrection->Get_ecore_CalibrationParameters().ReadFromFile("CEMC_ECORE_RECALIB", "xml", 0, 0,
460  //raw location
461  string(getenv("CALIBRATIONROOT")) + string("/CEMC/PositionRecalibration_EMCal_9deg_tilt/"));
462 
463  clusterCorrection->Verbosity(verbosity);
464  se->registerSubsystem(clusterCorrection);
465 
466  return;
467 }
468 void CEMC_Eval(const std::string &outputfile)
469 {
471 
473 
474  CaloEvaluator *eval = new CaloEvaluator("CEMCEVALUATOR", "CEMC", outputfile);
475  eval->Verbosity(3);
476  se->registerSubsystem(eval);
477 
478  return;
479 }
480 
481 void CEMC_QA()
482 {
484 
487  qa->Verbosity(verbosity);
488  se->registerSubsystem(qa);
489 
490  return;
491 }
492 
493 #endif