ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4_CEmc_EIC.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4_CEmc_EIC.C
1 #ifndef MACRO_G4CEMCEIC_C
2 #define MACRO_G4CEMCEIC_C
3 
4 #include <GlobalVariables.C>
5 
8 
11 
12 #include <g4eval/CaloEvaluator.h>
13 
14 #include <g4main/PHG4Reco.h>
15 
16 #include <caloreco/RawClusterBuilderGraph.h>
17 #include <caloreco/RawClusterBuilderTemplate.h>
18 #include <caloreco/RawTowerCalibration.h>
19 
20 #include <fun4all/Fun4AllServer.h>
21 
22 #include <cmath>
23 
24 R__LOAD_LIBRARY(libcalo_reco.so)
25 R__LOAD_LIBRARY(libg4calo.so)
26 R__LOAD_LIBRARY(libg4detectors.so)
27 R__LOAD_LIBRARY(libg4eval.so)
28 
29 namespace Enable
30 {
31  bool CEMC = false;
32  bool CEMC_ABSORBER = false;
33  bool CEMC_OVERLAPCHECK = false;
34  bool CEMC_CELL = false;
35  bool CEMC_TOWER = false;
36  bool CEMC_CLUSTER = false;
37  bool CEMC_EVAL = false;
38  int CEMC_VERBOSITY = 0;
39 } // namespace Enable
40 
41 namespace G4CEMC
42 {
43  double cemcdepth = 9;
44  // tungs to scint width ratio of ~10:1
45  // corresponds to approx 2% sampling fraction
46 
47  // 18 radiation lengths for 40 layers
48  double scint_width = 0.05;
49  double tungs_width = 0.245;
50  double electronics_width = 0.5;
51 
52  int min_cemc_layer = 1;
53  int max_cemc_layer = 41;
54 
55  double topradius = 106.8; // cm
56  double bottomradius = 95; // cm
57  double negrapidity = -1.5;
58  double posrapidity = 1.24;
59  // this is default set to -1.5<eta<1.24 for 2018 Letter of Intent
60  // if the user changes these, the z position of the
61  // calorimeter must be changed in the function CEmc(...)
62 
63  // Digitization (default photon digi):
65  // directly pass the energy of sim tower to digitized tower
66  // kNo_digitization
67  // simple digitization with photon statistics, single amplitude ADC conversion and pedestal
68  // kSimple_photon_digitization
69  // digitization with photon statistics on SiPM with an effective pixel N, ADC conversion and pedestal
70  // kSiPM_photon_digitization
71 
73  {
75 
77  };
78 
79  // default: template clusterizer, RawClusterBuilderTemplate, as developed by Sasha Bazilevsky
81  // graph clusterizer, RawClusterBuilderGraph
82  // enu_Cemc_clusterizer Cemc_clusterizer = kCemcGraphClusterizer;
83 } // namespace G4CEMC
84 
85 namespace CEMC_TOWER
86 {
87  double emin = NAN;
88 }
89 
90 // Black hole and size parameters set in CEmc function
91 void CEmcInit(const int nslats = 1)
92 {
93 }
94 
95 double CEmc(PHG4Reco *g4Reco, double radius)
96 {
97  bool AbsorberActive = Enable::ABSORBER || Enable::CEMC_ABSORBER;
98  bool OverlapCheck = Enable::OVERLAPCHECK || Enable::CEMC_OVERLAPCHECK;
99 
100  if (radius > 95)
101  {
102  cout << "inconsistency, radius: " << radius
103  << " larger than allowed inner radius for CEMC = 95 cm" << endl;
104  gSystem->Exit(-1);
105  }
106 
107  radius = 95;
108 
109  PHG4CylinderSubsystem *cemc;
110 
111  // determine the length of the calorimeter
112  // can adjust length coverage by just adjusting these values
113  // rapidity coverage will be determined by z shift of EMCAl
114  // as indicated in the loop below
115 
116  // eta = -ln(tan(theta/2))
117  double theta1 = 2. * TMath::ATan(TMath::Exp(-1 * G4CEMC::posrapidity));
118  double theta2 = 2. * TMath::ATan(TMath::Exp(-1 * G4CEMC::negrapidity));
119  // get the angle between the beam pipe and negative pseudorapidity axis
120  theta2 = M_PI - theta2;
121 
122  double z1 = G4CEMC::topradius / TMath::Tan(theta1);
123  double z2 = G4CEMC::topradius / TMath::Tan(theta2);
124 
125  double z3 = G4CEMC::bottomradius / TMath::Tan(theta1);
126  double z4 = G4CEMC::bottomradius / TMath::Tan(theta2);
127 
128  // this is the top layer length
129  double totaltoplength = z1 + z2;
130  // this is the bottom layer length
131  double totalbottomlength = z3 + z4;
132 
133  //Added by Barak, 12/12/19
134  double ztemp = 0;
135  double layer_shift = 0;
136 
137  double height = 0;
138  for (int thislayer = G4CEMC::min_cemc_layer; thislayer <= G4CEMC::max_cemc_layer;
139  thislayer++)
140  {
141  // the length for a particular layer is determined from the bottom length
142  double thislength = totalbottomlength + (height / TMath::Tan(theta1)) + (height / TMath::Tan(theta2));
143 
144  cemc = new PHG4CylinderSubsystem("ABSORBER_CEMC", thislayer);
145  cemc->set_double_param("radius", radius);
146  cemc->set_string_param("material", "Spacal_W_Epoxy");
147  cemc->set_double_param("thickness", G4CEMC::tungs_width);
148  cemc->set_double_param("length", thislength);
149  cemc->set_int_param("lengthviarapidity", 0);
150 
151  // starts centered around IP
152  // shift backwards 30 cm for total 370 cm length to cover -1.5<eta<1.24
153  //cemc->set_double_param("place_z", -30);
154 
155  //Modified by Barak, 12/12/19
156  ztemp = radius / TMath::Tan(theta2);
157  layer_shift = -1. * (ztemp - (thislength / 2.));
158  cemc->set_double_param("place_z", layer_shift);
159 
160  cemc->SuperDetector("ABSORBER_CEMC");
161  if (AbsorberActive) cemc->SetActive();
162  cemc->OverlapCheck(OverlapCheck);
163 
164  g4Reco->registerSubsystem(cemc);
165 
166  radius += G4CEMC::tungs_width;
167  radius += no_overlapp;
168 
169  height += G4CEMC::tungs_width;
170  height += no_overlapp; //Added by Barak, 12/13/19
171 
172  //Added by Barak, 12/13/19
173  thislength = totalbottomlength + (height / TMath::Tan(theta1)) + (height / TMath::Tan(theta2));
174 
175  cemc = new PHG4CylinderSubsystem("CEMC", thislayer);
176  cemc->set_double_param("radius", radius);
177  cemc->set_string_param("material", "PMMA");
178  cemc->set_double_param("thickness", G4CEMC::scint_width);
179  cemc->set_int_param("lightyield", 1);
180  cemc->set_int_param("lengthviarapidity", 0);
181  cemc->set_double_param("length", thislength);
182 
183  // shift back -30 cm to cover -1.4<eta<1.1
184  //cemc->set_double_param("place_z", -30);
185 
186  //Modified by Barak, 12/12/19
187  cemc->set_double_param("place_z", layer_shift);
188 
189  cemc->SuperDetector("CEMC");
190 
191  cemc->SetActive();
192  cemc->OverlapCheck(OverlapCheck);
193  g4Reco->registerSubsystem(cemc);
194 
195  radius += G4CEMC::scint_width;
196  radius += no_overlapp;
197 
198  height += G4CEMC::scint_width;
199  height += no_overlapp; //Added by Barak, 12/13/19
200  }
201 
202  PHG4CylinderSubsystem *cemc_cyl = new PHG4CylinderSubsystem("CEMC_ELECTRONICS", 0);
203  cemc_cyl->set_double_param("radius", radius);
204  cemc_cyl->set_string_param("material", "G4_TEFLON");
205  cemc_cyl->set_double_param("thickness", G4CEMC::electronics_width);
206 
207  double l1 = (radius + 0.5) / TMath::Tan(theta1);
208  double l2 = (radius + 0.5) / TMath::Tan(theta2);
209  cemc_cyl->set_int_param("lengthviarapidity", 0);
210  cemc_cyl->set_double_param("length", l1 + l2);
211  // shift back -30 cm to cover -1.4<eta<1.1
212  //cemc_cyl->set_double_param("place_z", -30);
213 
214  //Modified by Barak, 12/12/19
215  layer_shift = -1. * ((l2 - l1) / 2.);
216  cemc_cyl->set_double_param("place_z", layer_shift);
217 
218  if (AbsorberActive) cemc_cyl->SetActive();
219  cemc_cyl->OverlapCheck(OverlapCheck);
220  g4Reco->registerSubsystem(cemc_cyl);
221  // update black hole settings since we have the values here
223  BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, layer_shift + (l1 + l2) / 2.);
224  BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, layer_shift - (l1 + l2) / 2.);
225 
226  return radius;
227 }
228 
230 {
232 
234 
235  PHG4CylinderCellReco *cemc_cells = new PHG4CylinderCellReco("CEMCCYLCELLRECO");
236  cemc_cells->Detector("CEMC");
237  cemc_cells->Verbosity(verbosity);
238  double radius = 95;
239  for (int i = G4CEMC::min_cemc_layer; i <= G4CEMC::max_cemc_layer; i++)
240  {
241  //Added by Barak, 12/13/19
242  radius += (G4CEMC::tungs_width + no_overlapp);
243  if (i > 1) radius += (G4CEMC::scint_width + no_overlapp);
244 
245  cemc_cells->cellsize(i, 2. * M_PI / 256. * radius, 2. * M_PI / 256. * radius);
246  }
247  se->registerSubsystem(cemc_cells);
248 
249  return;
250 }
251 
253 {
255 
257 
258  RawTowerBuilder *CemcTowerBuilder = new RawTowerBuilder("EmcRawTowerBuilder");
259  CemcTowerBuilder->Detector("CEMC");
260  CemcTowerBuilder->set_sim_tower_node_prefix("SIM");
262  {
263  CemcTowerBuilder->EminCut(CEMC_TOWER::emin);
264  }
265  CemcTowerBuilder->Verbosity(verbosity);
266  se->registerSubsystem(CemcTowerBuilder);
267 
268  const double photoelectron_per_GeV = 500; // 500 photon per total GeV deposition
269  // just set a 4% sampling fraction - already tuned by tungs/scint width ratio
270  double sampling_fraction = 4e-02;
271  RawTowerDigitizer *CemcTowerDigitizer = new RawTowerDigitizer("EmcRawTowerDigitizer");
272  CemcTowerDigitizer->Detector("CEMC");
273  CemcTowerDigitizer->Verbosity(verbosity);
274  CemcTowerDigitizer->set_digi_algorithm(G4CEMC::TowerDigi);
275  CemcTowerDigitizer->set_pedstal_central_ADC(0);
276  CemcTowerDigitizer->set_pedstal_width_ADC(8); // eRD1 test beam setting
277  CemcTowerDigitizer->set_photonelec_ADC(1); // not simulating ADC discretization error
278  CemcTowerDigitizer->set_photonelec_yield_visible_GeV(photoelectron_per_GeV / sampling_fraction);
279  CemcTowerDigitizer->set_zero_suppression_ADC(16); // eRD1 test beam setting
280  se->registerSubsystem(CemcTowerDigitizer);
281 
282  RawTowerCalibration *CemcTowerCalibration = new RawTowerCalibration("EmcRawTowerCalibration");
283  CemcTowerCalibration->Detector("CEMC");
284  CemcTowerCalibration->Verbosity(verbosity);
287  {
288  CemcTowerCalibration->set_calib_const_GeV_ADC(1.0 / 0.023); // 2.3% sampling fraction from test beam
289  }
290  else
291  {
292  CemcTowerCalibration->set_calib_const_GeV_ADC(1. / photoelectron_per_GeV / 0.9715);
293  }
294  CemcTowerCalibration->set_pedstal_ADC(0);
295 
296  se->registerSubsystem(CemcTowerCalibration);
297 
298  return;
299 }
300 
302 {
304 
306 
308  {
309  RawClusterBuilderTemplate *cemc_clusterbuilder = new RawClusterBuilderTemplate("EmcRawClusterBuilderTemplate");
310  cemc_clusterbuilder->Detector("CEMC");
311  cemc_clusterbuilder->Verbosity(verbosity);
312 
313  cemc_clusterbuilder->set_threshold_energy(0.030); // This threshold should be the same as in CEMCprof_Thresh**.root file below
314  std::string femc_prof = getenv("CALIBRATIONROOT");
315  femc_prof += "/EmcProfile/CEMCprof_Thresh30MeV.root";
316  cemc_clusterbuilder->LoadProfile(femc_prof);
317  se->registerSubsystem(cemc_clusterbuilder);
318  }
320  {
321  RawClusterBuilderGraph *cemc_clusterbuilder = new RawClusterBuilderGraph("EmcRawClusterBuilderGraph");
322  cemc_clusterbuilder->Detector("CEMC");
323  cemc_clusterbuilder->Verbosity(verbosity);
324  se->registerSubsystem(cemc_clusterbuilder);
325  }
326  else
327  {
328  cout << "CEMC_Clusters - unknown clusterizer setting!! " << endl;
329  exit(1);
330  }
331  return;
332 }
333 void CEMC_Eval(const std::string &outputfile)
334 {
337  CaloEvaluator *eval = new CaloEvaluator("CEMCEVALUATOR", "CEMC", outputfile.c_str());
338  eval->Verbosity(verbosity);
339  se->registerSubsystem(eval);
340  return;
341 }
342 #endif