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_2DProjectiveSpacal(PHG4Reco *g4Reco, double radius, const int crossings);
31 
32 R__LOAD_LIBRARY(libcalo_reco.so)
33 R__LOAD_LIBRARY(libg4calo.so)
34 R__LOAD_LIBRARY(libg4detectors.so)
35 R__LOAD_LIBRARY(libg4eval.so)
36 R__LOAD_LIBRARY(libqa_modules.so)
37 
38 namespace Enable
39 {
40  bool CEMC = false;
41  bool CEMC_ABSORBER = false;
42  bool CEMC_OVERLAPCHECK = false;
43  bool CEMC_CELL = false;
44  bool CEMC_TOWER = false;
45  bool CEMC_CLUSTER = false;
46  bool CEMC_EVAL = false;
47  bool CEMC_QA = false;
48  bool CEMC_G4Hit = true;
49  int CEMC_VERBOSITY = 0;
50 } // namespace Enable
51 
52 namespace G4CEMC
53 {
54  int Min_cemc_layer = 1;
55  int Max_cemc_layer = 1;
56 
57  // Digitization (default photon digi):
59  // directly pass the energy of sim tower to digitized tower
60  // kNo_digitization
61  // simple digitization with photon statistics, single amplitude ADC conversion and pedestal
62  // kSimple_photon_digitization
63  // digitization with photon statistics on SiPM with an effective pixel N, ADC conversion and pedestal
64  // kSiPM_photon_digitization
65 
66  // 2D azimuthal projective SPACAL (slow)
68 
70  {
72 
74  };
75 
79  // enu_Cemc_clusterizer Cemc_clusterizer = kCemcGraphClusterizer;
80 
81 } // namespace G4CEMC
82 
83 // black hole parameters are set in CEmc function
84 // needs a dummy argument to play with current G4Setup_sPHENIX.C
85 void CEmcInit(const int i = 0)
86 {
87 }
88 
90 double
91 CEmc(PHG4Reco *g4Reco, double radius, const int crossings)
92 {
93  return CEmc_2DProjectiveSpacal(/*PHG4Reco**/ g4Reco, /*double*/ radius, /*const int */ crossings);
94 }
95 
97 double
98 CEmc_2DProjectiveSpacal(PHG4Reco *g4Reco, double radius, const int crossings)
99 {
100  bool AbsorberActive = Enable::ABSORBER || Enable::CEMC_ABSORBER;
101  bool OverlapCheck = Enable::OVERLAPCHECK || Enable::CEMC_OVERLAPCHECK;
102 
103  double emc_inner_radius = 92; // emc inner radius from engineering drawing
104  double cemcthickness = 22.50000 - no_overlapp;
105 
106  // max radius is 116 cm;
107  double emc_outer_radius = emc_inner_radius + cemcthickness; // outer radius
108  assert(emc_outer_radius < 116);
109 
110  if (radius > emc_inner_radius)
111  {
112  cout << "inconsistency: preshower radius+thickness: " << radius
113  << " larger than emc inner radius: " << emc_inner_radius << endl;
114  gSystem->Exit(-1);
115  }
116 
117  // the radii are only to determined the thickness of the cemc
118  radius = emc_inner_radius;
119 
120  // 1.5cm thick teflon as an approximation for EMCAl light collection + electronics (10% X0 total estimated)
121  PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("CEMC_ELECTRONICS", 0);
122  cyl->set_double_param("radius", radius);
123  cyl->set_string_param("material", "G4_TEFLON");
124  cyl->set_double_param("thickness", 1.5 - no_overlapp);
125  cyl->SuperDetector("CEMC_ELECTRONICS");
126  cyl->OverlapCheck(OverlapCheck);
127  if (AbsorberActive) cyl->SetActive();
128  g4Reco->registerSubsystem(cyl);
129 
130  radius += 1.5;
131  cemcthickness -= 1.5 + no_overlapp;
132 
133  // 0.5cm thick Stainless Steel as an approximation for EMCAl support system
134  cyl = new PHG4CylinderSubsystem("CEMC_SPT", 0);
135  cyl->SuperDetector("CEMC_SPT");
136  cyl->set_double_param("radius", radius + cemcthickness - 0.5);
137  cyl->set_string_param("material", "SS310"); // SS310 Stainless Steel
138  cyl->set_double_param("thickness", 0.5 - no_overlapp);
139  cyl->OverlapCheck(OverlapCheck);
140  if (AbsorberActive) cyl->SetActive();
141  g4Reco->registerSubsystem(cyl);
142 
143  // this is the z extend and outer radius of the support structure and therefore the z extend
144  // and radius of the surrounding black holes
145  double sptlen = PHG4Utils::GetLengthForRapidityCoverage(radius + cemcthickness);
148  BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, radius + cemcthickness);
149 
150  cemcthickness -= 0.5 + no_overlapp;
151 
152  int ilayer = 0;
154 
155  cemc = new PHG4SpacalSubsystem("CEMC", ilayer);
156 
157  cemc->set_int_param("virualize_fiber", 0);
158  cemc->set_int_param("azimuthal_seg_visible", 1);
159  cemc->set_int_param("construction_verbose", 0);
160  cemc->Verbosity(0);
162  cemc->SetCalibrationFileDir(string(getenv("CALIBRATIONROOT")) + string("/CEMC/Geometry_2023ProjTilted/"));
163  cemc->set_double_param("radius", radius); // overwrite minimal radius
164  cemc->set_double_param("thickness", cemcthickness); // overwrite thickness
166 
167  cemc->SetActive();
168  cemc->SuperDetector("CEMC");
169  if (AbsorberActive) cemc->SetAbsorberActive();
170  cemc->OverlapCheck(OverlapCheck);
171 
172  g4Reco->registerSubsystem(cemc);
173 
174  if (ilayer > G4CEMC::Max_cemc_layer)
175  {
176  cout << "layer discrepancy, current layer " << ilayer
177  << " max cemc layer: " << G4CEMC::Max_cemc_layer << endl;
178  }
179 
180  radius += cemcthickness;
181  radius += no_overlapp;
182 
183  return radius;
184 }
185 
187 {
189 
191 
193  {
194  PHG4CylinderCellReco *cemc_cells = new PHG4CylinderCellReco("CEMCCYLCELLRECO");
195  cemc_cells->Detector("CEMC");
196  cemc_cells->Verbosity(verbosity);
198  {
199  // cemc_cells->etaphisize(i, 0.024, 0.024);
200  const double radius = 95;
201  cemc_cells->cellsize(i, 2 * M_PI / 256. * radius, 2 * M_PI / 256. * radius);
202  }
203  se->registerSubsystem(cemc_cells);
204  }
206  {
207  if (!Enable::CEMC_G4Hit) return;
208  PHG4FullProjSpacalCellReco *cemc_cells = new PHG4FullProjSpacalCellReco("CEMCCYLCELLRECO");
209  cemc_cells->Detector("CEMC");
210  cemc_cells->Verbosity(verbosity);
212  string(getenv("CALIBRATIONROOT")) + string("/CEMC/LightCollection/Prototype3Module.xml"),
213  "data_grid_light_guide_efficiency", "data_grid_fiber_trans");
214  se->registerSubsystem(cemc_cells);
215  }
216  else
217  {
218  cout << "G4_CEmc_Spacal.C::CEmc - Fatal Error - unrecognized SPACAL configuration #"
219  << G4CEMC::Cemc_spacal_configuration << ". Force exiting..." << endl;
220  gSystem->Exit(-1);
221  return;
222  }
223 
224  return;
225 }
226 
228 {
230 
232  if (Enable::CEMC_G4Hit)
233  {
234  RawTowerBuilder *TowerBuilder = new RawTowerBuilder("EmcRawTowerBuilder");
235  TowerBuilder->Detector("CEMC");
236  TowerBuilder->set_sim_tower_node_prefix("SIM");
237  TowerBuilder->Verbosity(verbosity);
238  se->registerSubsystem(TowerBuilder);
239  }
240  double sampling_fraction = 1;
241  // 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
242  // 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
243  // sampling_fraction = 1.90951e-02; // 2017 Tilt porjective SPACAL, 8 GeV photon, eta = 0.3 - 0.4
244  sampling_fraction = 2e-02; // 2017 Tilt porjective SPACAL, tower-by-tower calibration
245  const double photoelectron_per_GeV = 500; // 500 photon per total GeV deposition
246 
247  RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("EmcRawTowerDigitizer");
248  TowerDigitizer->Detector("CEMC");
249  TowerDigitizer->Verbosity(verbosity);
250  TowerDigitizer->set_digi_algorithm(G4CEMC::TowerDigi);
251  TowerDigitizer->set_variable_pedestal(true); // read ped central and width from calibrations file comment next 2 lines if true
252  // TowerDigitizer->set_pedstal_central_ADC(0);
253  // TowerDigitizer->set_pedstal_width_ADC(8); // eRD1 test beam setting
254  TowerDigitizer->set_photonelec_ADC(1); // not simulating ADC discretization error
255  TowerDigitizer->set_photonelec_yield_visible_GeV(photoelectron_per_GeV / sampling_fraction);
256  TowerDigitizer->set_variable_zero_suppression(true); // read zs values from calibrations file comment next line if true
257  // TowerDigitizer->set_zero_suppression_ADC(16); // eRD1 test beam setting
258  if (!Enable::CEMC_G4Hit) TowerDigitizer->set_towerinfo(RawTowerDigitizer::ProcessTowerType::kTowerInfoOnly); // just use towerinfo
259  if (Enable::CDB)
260  {
261  TowerDigitizer->GetParameters().ReadFromCDB("EMCTOWERCALIB");
262  }
263  else
264  {
265  TowerDigitizer->GetParameters().ReadFromFile("CEMC", "xml", 0, 0,
266  string(getenv("CALIBRATIONROOT")) + string("/CEMC/TowerCalibCombinedParams_2020/")); // calibration database
267  }
268  se->registerSubsystem(TowerDigitizer);
269 
270  RawTowerCalibration *TowerCalibration = new RawTowerCalibration("EmcRawTowerCalibration");
271  TowerCalibration->Detector("CEMC");
272  TowerCalibration->Verbosity(verbosity);
273  if (!Enable::CEMC_G4Hit) TowerCalibration->set_towerinfo(RawTowerCalibration::ProcessTowerType::kTowerInfoOnly); // just use towerinfo
275  {
276  // just use sampling fraction set previously
277  TowerCalibration->set_calib_const_GeV_ADC(1.0 / sampling_fraction);
278  }
279  else
280  {
282  if (Enable::CDB)
283  {
284  TowerCalibration->GetCalibrationParameters().ReadFromCDB("EMCTOWERCALIB");
285  }
286  else
287  {
288  TowerCalibration->GetCalibrationParameters().ReadFromFile("CEMC", "xml", 0, 0,
289  string(getenv("CALIBRATIONROOT")) + string("/CEMC/TowerCalibCombinedParams_2020/")); // calibration database
290  }
291  TowerCalibration->set_variable_GeV_ADC(true); // read GeV per ADC from calibrations file comment next line if true
292  // TowerCalibration->set_calib_const_GeV_ADC(1. / photoelectron_per_GeV / 0.9715); // overall energy scale based on 4-GeV photon simulations
293  TowerCalibration->set_variable_pedestal(true); // read pedestals from calibrations file comment next line if true
294  // TowerCalibration->set_pedstal_ADC(0);
295  }
296  se->registerSubsystem(TowerCalibration);
297 
298  return;
299 }
300 
302 {
304 
306 
308  {
309  RawClusterBuilderTemplate *ClusterBuilder = new RawClusterBuilderTemplate("EmcRawClusterBuilderTemplate");
310  ClusterBuilder->Detector("CEMC");
311  ClusterBuilder->Verbosity(verbosity);
312  ClusterBuilder->set_threshold_energy(0.030); // This threshold should be the same as in CEMCprof_Thresh**.root file below
313  std::string emc_prof = getenv("CALIBRATIONROOT");
314  emc_prof += "/EmcProfile/CEMCprof_Thresh30MeV.root";
315  ClusterBuilder->LoadProfile(emc_prof);
316  if (!Enable::CEMC_G4Hit) ClusterBuilder->set_UseTowerInfo(1); // just use towerinfo
317  // ClusterBuilder->set_UseTowerInfo(1); // to use towerinfo objects rather than old RawTower
318  se->registerSubsystem(ClusterBuilder);
319  }
321  {
322  RawClusterBuilderGraph *ClusterBuilder = new RawClusterBuilderGraph("EmcRawClusterBuilderGraph");
323  ClusterBuilder->Detector("CEMC");
324  ClusterBuilder->Verbosity(verbosity);
325  se->registerSubsystem(ClusterBuilder);
326  }
327  else
328  {
329  cout << "CEMC_Clusters - unknown clusterizer setting!" << endl;
330  exit(1);
331  }
332 
333  RawClusterPositionCorrection *clusterCorrection = new RawClusterPositionCorrection("CEMC");
334  if (!Enable::CEMC_G4Hit) clusterCorrection->set_UseTowerInfo(1); // just use towerinfo
335  // clusterCorrection->set_UseTowerInfo(1); // to use towerinfo objects rather than old RawTower
336 
337 // if (Enable::CDB)
338 // {
339 // clusterCorrection->Get_eclus_CalibrationParameters().ReadFromCDB("CEMCRECALIB");
340 // clusterCorrection->Get_ecore_CalibrationParameters().ReadFromCDB("CEMC_ECORE_RECALIB");
341 // }
342 // else
343 // {
344 // clusterCorrection->Get_eclus_CalibrationParameters().ReadFromFile("CEMC_RECALIB", "xml", 0, 0,
345 // // raw location
346 // string(getenv("CALIBRATIONROOT")) + string("/CEMC/PositionRecalibration_EMCal_9deg_tilt/"));
347 // clusterCorrection->Get_ecore_CalibrationParameters().ReadFromFile("CEMC_ECORE_RECALIB", "xml", 0, 0,
348 // // raw location
349 // string(getenv("CALIBRATIONROOT")) + string("/CEMC/PositionRecalibration_EMCal_9deg_tilt/"));
350 // }
351 
352  clusterCorrection->Verbosity(verbosity);
353  se->registerSubsystem(clusterCorrection);
354 
355  return;
356 }
357 void CEMC_Eval(const std::string &outputfile)
358 {
360 
362 
363  CaloEvaluator *eval = new CaloEvaluator("CEMCEVALUATOR", "CEMC", outputfile);
364  eval->Verbosity(verbosity);
365  se->registerSubsystem(eval);
366 
367  return;
368 }
369 
370 void CEMC_QA()
371 {
373 
376  qa->Verbosity(verbosity);
377  se->registerSubsystem(qa);
378 
379  return;
380 }
381 
382 #endif