macros/leptoquarks_fun4all/G4_CEmc_Spacal.C

0001
0002 int Min_cemc_layer = 1;
0003 int Max_cemc_layer = 1;
0004
0005 // set a default value for SPACAL configuration
0006 //  // 1D azimuthal projective SPACAL (fast)
0007 //int Cemc_spacal_configuration = PHG4CylinderGeom_Spacalv1::k1DProjectiveSpacal;
0008 //   2D azimuthal projective SPACAL (slow)
0009 int Cemc_spacal_configuration = PHG4CylinderGeom_Spacalv1::k2DProjectiveSpacal;
0010
0011 enum enu_Cemc_clusterizer
0012 {
0013   kCemcGraphClusterizer,
0014
0015   kCemcTemplateClusterizer
0016 };
0017
0018 //! template clusterizer, RawClusterBuilderTemplate, as developed by Sasha Bazilevsky
0019 enu_Cemc_clusterizer Cemc_clusterizer = kCemcTemplateClusterizer;
0020 //! graph clusterizer, RawClusterBuilderGraph
0021 //enu_Cemc_clusterizer Cemc_clusterizer = kCemcGraphClusterizer;
0022
0023 #include <iostream>
0024
0025 // just a dummy parameter used by the tilted plate geom
0026 void CEmcInit(const int nslats = 1)
0027 {
0028   Min_cemc_layer = 1;
0029   Max_cemc_layer = 1;
0030 }
0031
0032 //! EMCal main setup macro
0033 double
0034 CEmc(PHG4Reco *g4Reco, double radius, const int crossings,
0035      const int absorberactive = 0)
0036 {
0037   if (Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k1DProjectiveSpacal)
0038   {
0039     return CEmc_1DProjectiveSpacal(/*PHG4Reco**/ g4Reco, /*double*/ radius, /*const int */
0040                                    crossings, /*const int*/ absorberactive);
0041   }
0042   else if (Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k2DProjectiveSpacal)
0043   {
0044     return CEmc_2DProjectiveSpacal(/*PHG4Reco**/ g4Reco, /*double*/ radius, /*const int */
0045                                    crossings, /*const int*/ absorberactive);
0046   }
0047   else
0048   {
0049     std::cout
0050         << "G4_CEmc_Spacal.C::CEmc - Fatal Error - unrecognized SPACAL configuration #"
0051         << Cemc_spacal_configuration << ". Force exiting..." << std::endl;
0052     exit(-1);
0053     return 0;
0054   }
0055 }
0056
0057 //! EMCal setup macro - 1D azimuthal projective SPACAL
0058 double
0059 CEmc_1DProjectiveSpacal(PHG4Reco *g4Reco, double radius, const int crossings, const int absorberactive = 0)
0060 {
0061   double emc_inner_radius = 95.;  // emc inner radius from engineering drawing
0062   double cemcthickness = 12.7;
0063   double emc_outer_radius = emc_inner_radius + cemcthickness;  // outer radius
0064
0065   if (radius > emc_inner_radius)
0066   {
0067     cout << "inconsistency: pstof outer radius: " << radius
0068          << " larger than emc inner radius: " << emc_inner_radius
0069          << endl;
0070     gSystem->Exit(-1);
0071   }
0072
0073   //---------------
0074   // Load libraries
0075   //---------------
0076
0077   gSystem->Load("libg4detectors.so");
0078   gSystem->Load("libg4testbench.so");
0079
0080   //  boundary check
0081   if (radius > emc_inner_radius - 1.5 - no_overlapp)
0082   {
0083     cout << "G4_CEmc_Spacal.C::CEmc() - expect radius < " << emc_inner_radius - 1.5 - no_overlapp << " to install SPACAL" << endl;
0084     exit(1);
0085   }
0086   radius = emc_inner_radius - 1.5 - no_overlapp;
0087
0088   // 1.5cm thick teflon as an approximation for EMCAl light collection + electronics (10% X0 total estimated)
0089   PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("CEMC_ELECTRONICS", 0);
0090   cyl->SuperDetector("CEMC_ELECTRONICS");
0091   cyl->set_double_param("radius", radius);
0092   cyl->set_string_param("material", "G4_TEFLON");
0093   cyl->set_double_param("thickness", 1.5);
0094   if (absorberactive) cyl->SetActive();
0095   g4Reco->registerSubsystem(cyl);
0096
0097   radius += 1.5;
0098   radius += no_overlapp;
0099
0100   int ilayer = Min_cemc_layer;
0101   PHG4SpacalSubsystem *cemc;
0102   cemc = new PHG4SpacalSubsystem("CEMC", ilayer);
0103   cemc->set_double_param("radius",emc_inner_radius);
0104   cemc->set_double_param("thickness", cemcthickness);
0105
0106   cemc->SetActive();
0107   cemc->SuperDetector("CEMC");
0108   if (absorberactive) cemc->SetAbsorberActive();
0109   cemc->OverlapCheck(overlapcheck);
0110
0111   g4Reco->registerSubsystem(cemc);
0112
0113   if (ilayer > Max_cemc_layer)
0114   {
0115     cout << "layer discrepancy, current layer " << ilayer
0116          << " max cemc layer: " << Max_cemc_layer << endl;
0117   }
0118
0119   radius += cemcthickness;
0120   radius += no_overlapp;
0121
0122   // 0.5cm thick Stainless Steel as an approximation for EMCAl support system
0123   cyl = new PHG4CylinderSubsystem("CEMC_SPT", 0);
0124   cyl->SuperDetector("CEMC_SPT");
0125   cyl->set_double_param("radius", radius);
0126   cyl->set_string_param("material", "SS310");  // SS310 Stainless Steel
0127   cyl->set_double_param("thickness", 0.5);
0128   if (absorberactive)
0129     cyl->SetActive();
0130   g4Reco->registerSubsystem(cyl);
0131
0132   radius += 0.5;
0133   radius += no_overlapp;
0134
0135   return radius;
0136 }
0137
0138 //! 2D full projective SPACAL
0139 double
0140 CEmc_2DProjectiveSpacal(PHG4Reco *g4Reco, double radius, const int crossings,
0141                         const int absorberactive = 0)
0142 {
0143   double emc_inner_radius = 92;  // emc inner radius from engineering drawing
0144   double cemcthickness = 24.00000 - no_overlapp;
0145
0146   //max radius is 116 cm;
0147   double emc_outer_radius = emc_inner_radius + cemcthickness;  // outer radius
0148   assert(emc_outer_radius < 116);
0149
0150   if (radius > emc_inner_radius)
0151   {
0152     cout << "inconsistency: preshower radius+thickness: " << radius
0153          << " larger than emc inner radius: " << emc_inner_radius << endl;
0154     gSystem->Exit(-1);
0155   }
0156
0157   //---------------
0158   // Load libraries
0159   //---------------
0160
0161   gSystem->Load("libg4detectors.so");
0162
0163   // the radii are only to determined the thickness of the cemc
0164   radius = emc_inner_radius;
0165
0166   //---------------
0167   // Load libraries
0168   //---------------
0169
0170   // 1.5cm thick teflon as an approximation for EMCAl light collection + electronics (10% X0 total estimated)
0171   PHG4CylinderSubsystem *cyl = new PHG4CylinderSubsystem("CEMC_ELECTRONICS", 0);
0172   cyl->set_double_param("radius", radius);
0173   cyl->set_string_param("material", "G4_TEFLON");
0174   cyl->set_double_param("thickness", 1.5 - no_overlapp);
0175   cyl->SuperDetector("CEMC_ELECTRONICS");
0176   cyl->OverlapCheck(overlapcheck);
0177   if (absorberactive) cyl->SetActive();
0178   g4Reco->registerSubsystem(cyl);
0179
0180   radius += 1.5;
0181   cemcthickness -= 1.5 + no_overlapp;
0182
0183   // 0.5cm thick Stainless Steel as an approximation for EMCAl support system
0184   cyl = new PHG4CylinderSubsystem("CEMC_SPT", 0);
0185   cyl->SuperDetector("CEMC_SPT");
0186   cyl->set_double_param("radius", radius + cemcthickness - 0.5);
0187   cyl->set_string_param("material", "SS310");  // SS310 Stainless Steel
0188   cyl->set_double_param("thickness", 0.5 - no_overlapp);
0189   cyl->OverlapCheck(overlapcheck);
0190   if (absorberactive)
0191     cyl->SetActive();
0192   g4Reco->registerSubsystem(cyl);
0193
0194   cemcthickness -= 0.5 + no_overlapp;
0195
0196   //---------------
0197   // Load libraries
0198   //---------------
0199
0200   int ilayer = 0;
0201   PHG4SpacalSubsystem *cemc;
0202
0203   const bool use_2015_design = false;
0204   if (use_2015_design)
0205   {
0206     cemc = new PHG4SpacalSubsystem("CEMC", ilayer);
0207
0208     cemc->set_int_param("config", PHG4CylinderGeom_Spacalv1::kFullProjective_2DTaper_SameLengthFiberPerTower);
0209     cemc->set_double_param("radius", radius);            // overwrite minimal radius
0210     cemc->set_double_param("thickness", cemcthickness);  // overwrite thickness
0211     cemc->set_int_param("azimuthal_n_sec", 32);
0212     //    cemc->set_int_param("construction_verbose", 2);
0213
0214     cemc->SetActive();
0215     cemc->SuperDetector("CEMC");
0216     if (absorberactive)
0217       cemc->SetAbsorberActive();
0218     cemc->OverlapCheck(overlapcheck);
0219   }
0220
0221   else
0222   {
0223     cemc = new PHG4SpacalSubsystem("CEMC", ilayer);
0224
0225     cemc->set_int_param("virualize_fiber", 0);
0226     cemc->set_int_param("azimuthal_seg_visible", 1);
0227     cemc->set_int_param("construction_verbose", 0);
0228     cemc->Verbosity(0);
0229
0230     cemc->UseCalibFiles(PHG4DetectorSubsystem::xml);
0231     cemc->SetCalibrationFileDir(string(getenv("CALIBRATIONROOT")) + string("/CEMC/Geometry_2017ProjTilted/"));
0232     cemc->set_double_param("radius", radius);            // overwrite minimal radius
0233     cemc->set_double_param("thickness", cemcthickness);  // overwrite thickness
0234
0235     cemc->SetActive();
0236     cemc->SuperDetector("CEMC");
0237     if (absorberactive)
0238       cemc->SetAbsorberActive();
0239     cemc->OverlapCheck(overlapcheck);
0240   }
0241
0242   g4Reco->registerSubsystem(cemc);
0243
0244   if (ilayer > Max_cemc_layer)
0245   {
0246     cout << "layer discrepancy, current layer " << ilayer
0247          << " max cemc layer: " << Max_cemc_layer << endl;
0248   }
0249
0250   radius += cemcthickness;
0251   radius += no_overlapp;
0252
0253   return radius;
0254 }
0255
0256 void CEMC_Cells(int verbosity = 0)
0257 {
0258   gSystem->Load("libfun4all.so");
0259   gSystem->Load("libg4detectors.so");
0260   Fun4AllServer *se = Fun4AllServer::instance();
0261
0262   if (Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k1DProjectiveSpacal)
0263   {
0264     PHG4CylinderCellReco *cemc_cells = new PHG4CylinderCellReco("CEMCCYLCELLRECO");
0265     cemc_cells->Detector("CEMC");
0266     cemc_cells->Verbosity(verbosity);
0267     for (int i = Min_cemc_layer; i <= Max_cemc_layer; i++)
0268     {
0269       //          cemc_cells->etaphisize(i, 0.024, 0.024);
0270       const double radius = 95;
0271       cemc_cells->cellsize(i, 2 * TMath::Pi() / 256. * radius, 2 * TMath::Pi() / 256. * radius);
0272
0273     }
0274     se->registerSubsystem(cemc_cells);
0275   }
0276   else if (Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k2DProjectiveSpacal)
0277   {
0278     PHG4FullProjSpacalCellReco *cemc_cells = new PHG4FullProjSpacalCellReco("CEMCCYLCELLRECO");
0279     cemc_cells->Detector("CEMC");
0280     cemc_cells->Verbosity(verbosity);
0281     cemc_cells->get_light_collection_model().load_data_file(
0282         string(getenv("CALIBRATIONROOT")) + string("/CEMC/LightCollection/Prototype3Module.xml"),
0283         "data_grid_light_guide_efficiency", "data_grid_fiber_trans");
0284     se->registerSubsystem(cemc_cells);
0285   }
0286   else
0287   {
0288     std::cout
0289         << "G4_CEmc_Spacal.C::CEmc - Fatal Error - unrecognized SPACAL configuration #"
0290         << Cemc_spacal_configuration << ". Force exiting..." << std::endl;
0291     exit(-1);
0292     return;
0293   }
0294
0295   return;
0296 }
0297
0298 void CEMC_Towers(int verbosity = 0)
0299 {
0300   gSystem->Load("libg4calo.so");
0301   gSystem->Load("libcalo_reco.so");
0302   Fun4AllServer *se = Fun4AllServer::instance();
0303
0304   RawTowerBuilder *TowerBuilder = new RawTowerBuilder("EmcRawTowerBuilder");
0305   TowerBuilder->Detector("CEMC");
0306   TowerBuilder->set_sim_tower_node_prefix("SIM");
0307   TowerBuilder->Verbosity(verbosity);
0308   se->registerSubsystem(TowerBuilder);
0309
0310   double sampling_fraction = 1;
0311   if (Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k1DProjectiveSpacal)
0312   {
0313     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
0314   }
0315   else if (Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k2DProjectiveSpacal)
0316   {
0317     //      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
0318 //    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
0319 //    sampling_fraction = 1.90951e-02; // 2017 Tilt porjective SPACAL, 8 GeV photon, eta = 0.3 - 0.4
0320     sampling_fraction = 2e-02; // 2017 Tilt porjective SPACAL, tower-by-tower calibration
0321   }
0322   else
0323   {
0324     std::cout
0325         << "G4_CEmc_Spacal.C::CEMC_Towers - Fatal Error - unrecognized SPACAL configuration #"
0326         << Cemc_spacal_configuration << ". Force exiting..." << std::endl;
0327     exit(-1);
0328     return;
0329   }
0330
0331   const double photoelectron_per_GeV = 500;  //500 photon per total GeV deposition
0332
0333   RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("EmcRawTowerDigitizer");
0334   TowerDigitizer->Detector("CEMC");
0335   TowerDigitizer->Verbosity(verbosity);
0336   TowerDigitizer->set_digi_algorithm(RawTowerDigitizer::kSimple_photon_digitization);
0337   TowerDigitizer->set_pedstal_central_ADC(0);
0338   TowerDigitizer->set_pedstal_width_ADC(8);  // eRD1 test beam setting
0339   TowerDigitizer->set_photonelec_ADC(1);     //not simulating ADC discretization error
0340   TowerDigitizer->set_photonelec_yield_visible_GeV(photoelectron_per_GeV / sampling_fraction);
0341   TowerDigitizer->set_zero_suppression_ADC(16);  // eRD1 test beam setting
0342   se->registerSubsystem(TowerDigitizer);
0343
0344   if (Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k1DProjectiveSpacal)
0345   {
0346     RawTowerCalibration *TowerCalibration = new RawTowerCalibration("EmcRawTowerCalibration");
0347     TowerCalibration->Detector("CEMC");
0348     TowerCalibration->Verbosity(verbosity);
0349     TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0350     TowerCalibration->set_calib_const_GeV_ADC(1. / photoelectron_per_GeV);
0351     TowerCalibration->set_pedstal_ADC(0);
0352     se->registerSubsystem(TowerCalibration);
0353   }
0354   else if (Cemc_spacal_configuration == PHG4CylinderGeom_Spacalv1::k2DProjectiveSpacal)
0355   {
0356     RawTowerCalibration *TowerCalibration = new RawTowerCalibration("EmcRawTowerCalibration");
0357     TowerCalibration->Detector("CEMC");
0358     TowerCalibration->Verbosity(verbosity);
0359     TowerCalibration->set_calib_algorithm(RawTowerCalibration::kTower_by_tower_calibration);
0360     TowerCalibration->GetCalibrationParameters().ReadFromFile("CEMC","xml",0,0,
0361         string(getenv("CALIBRATIONROOT")) + string("/CEMC/TowerCalib_2017ProjTilted/")); // calibration database
0362     TowerCalibration->set_calib_const_GeV_ADC(1. / photoelectron_per_GeV / 0.9715 ); // overall energy scale based on 4-GeV photon simulations
0363     TowerCalibration->set_pedstal_ADC(0);
0364     se->registerSubsystem(TowerCalibration);
0365   }
0366   else
0367   {
0368     std::cout
0369         << "G4_CEmc_Spacal.C::CEMC_Towers - Fatal Error - unrecognized SPACAL configuration #"
0370         << Cemc_spacal_configuration << ". Force exiting..." << std::endl;
0371     exit(-1);
0372     return;
0373   }
0374
0375   return;
0376 }
0377
0378 void CEMC_Clusters(int verbosity = 0)
0379 {
0380   gSystem->Load("libcalo_reco.so");
0381   Fun4AllServer *se = Fun4AllServer::instance();
0382
0383   if (Cemc_clusterizer == kCemcTemplateClusterizer)
0384   {
0385     RawClusterBuilderTemplate *ClusterBuilder = new RawClusterBuilderTemplate("EmcRawClusterBuilderTemplate");
0386     ClusterBuilder->Detector("CEMC");
0387     ClusterBuilder->Verbosity(verbosity);
0388     se->registerSubsystem(ClusterBuilder);
0389   }
0390   else if (Cemc_clusterizer == kCemcGraphClusterizer)
0391   {
0392     RawClusterBuilderGraph *ClusterBuilder = new RawClusterBuilderGraph("EmcRawClusterBuilderGraph");
0393     ClusterBuilder->Detector("CEMC");
0394     ClusterBuilder->Verbosity(verbosity);
0395     se->registerSubsystem(ClusterBuilder);
0396   }
0397   else
0398   {
0399     cout <<"CEMC_Clusters - unknown clusterizer setting!"<<endl;
0400     exit(1);
0401   }
0402
0403
0404   RawClusterPositionCorrection *clusterCorrection = new RawClusterPositionCorrection("CEMC");
0405
0406     clusterCorrection->Get_eclus_CalibrationParameters().ReadFromFile("CEMC_RECALIB","xml",0,0,
0407                             //raw location
0408                             string(getenv("CALIBRATIONROOT"))+string("/CEMC/PositionRecalibration/"));
0409   clusterCorrection->Get_ecore_CalibrationParameters().ReadFromFile("CEMC_ECORE_RECALIB","xml",0,0,
0410                                //raw location
0411                                string(getenv("CALIBRATIONROOT"))+string("/CEMC/PositionRecalibration"));
0412
0413   clusterCorrection->Verbosity(verbosity);
0414   se->registerSubsystem(clusterCorrection);
0415
0416   return;
0417 }
0418 void CEMC_Eval(std::string outputfile, int verbosity = 0)
0419 {
0420   gSystem->Load("libfun4all.so");
0421   gSystem->Load("libg4eval.so");
0422   Fun4AllServer *se = Fun4AllServer::instance();
0423
0424   CaloEvaluator *eval = new CaloEvaluator("CEMCEVALUATOR", "CEMC", outputfile.c_str());
0425   eval->Verbosity(verbosity);
0426   se->registerSubsystem(eval);
0427
0428   return;
0429 }