Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:20:35 

0001 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,00,0)
0002 #include <caloreco/RawTowerCalibration.h>
0003 
0004 #include <fun4all/SubsysReco.h>
0005 #include <fun4all/Fun4AllDstOutputManager.h>
0006 #include <fun4all/Fun4AllDummyInputManager.h>
0007 #include <fun4all/Fun4AllInputManager.h>
0008 #include <fun4all/Fun4AllOutputManager.h>
0009 #include <fun4all/Fun4AllServer.h>
0010 
0011 #include <g4calo/RawTowerBuilder.h>
0012 #include <g4calo/RawTowerDigitizer.h>
0013 
0014 #include <g4caloprototype/Prototype2RawTowerBuilder.h>
0015 #include <g4detectors/PHG4BlockSubsystem.h>
0016 #include <g4detectors/PHG4FullProjSpacalCellReco.h>
0017 #include <g4caloprototype/PHG4Prototype2HcalCellReco.h>
0018 #include <g4caloprototype/PHG4Prototype2InnerHcalSubsystem.h>
0019 #include <g4caloprototype/PHG4Prototype2OuterHcalSubsystem.h>
0020 #include <g4caloprototype/PHG4SpacalPrototypeSubsystem.h>
0021 
0022 #include <g4eval/PHG4DSTReader.h>
0023 
0024 #include <g4histos/G4HitNtuple.h>
0025 
0026 #include <g4main/PHG4ParticleGun.h>
0027 #include <g4main/PHG4Reco.h>
0028 #include <g4main/PHG4ParticleGenerator.h>
0029 #include <g4main/PHG4SimpleEventGenerator.h>
0030 #include <g4main/PHG4TruthSubsystem.h>
0031 
0032 #include <phool/recoConsts.h>
0033 
0034 #include <qa_modules/QAHistManagerDef.h>
0035 #include <qa_modules/QAG4SimulationCalorimeter.h>
0036 
0037 R__LOAD_LIBRARY(libcalo_reco.so)
0038 R__LOAD_LIBRARY(libfun4all.so)
0039 R__LOAD_LIBRARY(libg4calo.so)
0040 R__LOAD_LIBRARY(libg4caloprototype.so)
0041 R__LOAD_LIBRARY(libg4eval.so)
0042 R__LOAD_LIBRARY(libg4histos.so)
0043 R__LOAD_LIBRARY(libg4testbench.so)
0044 R__LOAD_LIBRARY(libqa_modules.so)
0045 
0046 #endif
0047 
0048 void Fun4All_G4_Prototype3(int nEvents = 1)
0049 {
0050 
0051   gSystem->Load("libfun4all");
0052   gSystem->Load("libg4caloprototype");
0053   gSystem->Load("libg4testbench");
0054   gSystem->Load("libg4histos");
0055   gSystem->Load("libg4eval.so");
0056   gSystem->Load("libqa_modules");
0057   gSystem->Load("libg4calo");
0058   gSystem->Load("libcalo_reco");
0059 
0060   bool cemc_on = true;
0061   bool cemc_cell = cemc_on && true;
0062   bool cemc_twr = cemc_cell && true;
0063   bool cemc_digi = cemc_twr && true;
0064   bool cemc_twrcal = cemc_digi && true;
0065   bool ihcal_on = true;
0066   bool ihcal_cell = ihcal_on && false;
0067   bool ihcal_twr = ihcal_cell && false;
0068   bool ihcal_digi = ihcal_twr && false;
0069   bool ihcal_twrcal = ihcal_digi && false;
0070   bool ohcal_on = true;
0071   bool ohcal_cell = ohcal_on && false;
0072   bool ohcal_twr = ohcal_cell && false;
0073   bool ohcal_digi = ohcal_twr && false;
0074   bool ohcal_twrcal =  ohcal_digi && false;
0075   bool cryo_on = true;
0076   bool bh_on = false; // the surrounding boxes need some further thinking
0077   bool dstreader = true;
0078   bool hit_ntuple = false;
0079   bool dstoutput = false;
0080 
0081   ///////////////////////////////////////////
0082   // Make the Server
0083   //////////////////////////////////////////
0084   Fun4AllServer *se = Fun4AllServer::instance();
0085   se->Verbosity(1);
0086   recoConsts *rc = recoConsts::instance();
0087   // only set this if you want a fixed random seed to make
0088   // results reproducible for testing
0089 //    rc->set_IntFlag("RANDOMSEED",12345678);
0090 
0091   // simulated setup sits at eta=1, theta=40.395 degrees
0092   double theta = 90-46.4;
0093   // shift in x with respect to midrapidity setup
0094   double add_place_x = 183.-173.93+2.54/2.;
0095   // Test beam generator
0096   PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
0097   gen->add_particles("e-", 1); // mu-,e-,anti_proton,pi-
0098   gen->set_vertex_distribution_mean(0.0, 0.0, 0);
0099   gen->set_vertex_distribution_width(0.0, .7, .7); // Rough beam profile size @ 16 GeV measured by Abhisek
0100   gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Gaus,
0101                     PHG4SimpleEventGenerator::Gaus, 
0102                                         PHG4SimpleEventGenerator::Gaus); // Gauss beam profile
0103   double angle = theta*TMath::Pi()/180.;
0104   double eta = -1.*TMath::Log(TMath::Tan(angle/2.));
0105   gen->set_eta_range(eta-0.001,eta+0.001); // 1mrad angular divergence
0106   gen->set_phi_range(-0.001, 0.001); // 1mrad angular divergence
0107   const double momentum = 32;
0108   gen->set_p_range(momentum,momentum, momentum*2e-2); // 2% momentum smearing
0109   se->registerSubsystem(gen);
0110 
0111   PHG4ParticleGenerator *pgen = new PHG4ParticleGenerator();
0112   pgen->set_name("geantino");
0113   //pgen->set_name(particle);
0114   pgen->set_vtx(0, 0, 0);
0115   //pgen->set_vtx(0, ypos, 0);
0116   angle = theta*TMath::Pi()/180.;
0117   eta = -1.*TMath::Log(TMath::Tan(angle/2.));
0118   pgen->set_eta_range(0.2*eta, 1.8*eta);
0119   //pgen->set_phi_range(-0.001, 0.001); // 1mrad angular diverpgence
0120   //pgen->set_phi_range(-0.5/180.*TMath::Pi(), 0.5/180.*TMath::Pi());
0121   //pgen->set_eta_range(-1., 1.);
0122   //pgen->set_phi_range(-0./180.*TMath::Pi(), 0./180.*TMath::Pi());
0123   pgen->set_phi_range(-20/180.*TMath::Pi(), 20/180.*TMath::Pi());
0124   pgen->set_mom_range(1, 1);
0125   //  se->registerSubsystem(pgen);
0126 
0127   // Simple single particle generator
0128   PHG4ParticleGun *gun = new PHG4ParticleGun();
0129   gun->set_name("geantino");
0130   //  gun->set_name("proton");
0131   gun->set_vtx(0, 0, 0);
0132   angle = theta*TMath::Pi()/180.;
0133   gun->set_mom(sin(angle),0.,cos(angle));
0134 //  se->registerSubsystem(gun);
0135 
0136 
0137   PHG4Reco* g4Reco = new PHG4Reco();
0138   g4Reco->set_field(0);
0139   //  g4Reco->SetPhysicsList("QGSP_BERT_HP"); // uncomment this line to enable the high-precision neutron simulation physics list, QGSP_BERT_HP
0140 
0141   //----------------------------------------
0142   // EMCal G4
0143   //----------------------------------------
0144   if (cemc_on)
0145     {
0146       PHG4SpacalPrototypeSubsystem *cemc;
0147       cemc = new PHG4SpacalPrototypeSubsystem("CEMC");
0148       cemc->SetActive();
0149       cemc->SuperDetector("CEMC");
0150       cemc->SetAbsorberActive();
0151       cemc->OverlapCheck(true);
0152 //      cemc->Verbosity(2);
0153 //      cemc->set_int_param("construction_verbose",2);
0154       cemc->UseCalibFiles(PHG4DetectorSubsystem::xml);
0155       cemc->SetCalibrationFileDir(string(getenv("CALIBRATIONROOT")) + string("/Prototype3/Geometry/") );
0156 //      cemc->SetCalibrationFileDir("./test_geom/" );
0157       //  cemc->set_double_param("z_rotation_degree", 15); // rotation around CG
0158 //      cemc->set_double_param("xpos", (116.77 + 137.0)*.5 - 26.5 - 10.2); // location in cm of EMCal CG. Updated with final positioning of EMCal
0159 //      cemc->set_double_param("ypos", 4); // put it some where in UIUC blocks
0160 //      cemc->set_double_param("zpos", 4); // put it some where in UIUC blocks
0161       g4Reco->registerSubsystem(cemc);
0162     }
0163   //----------------------------------------
0164   // HCal G4
0165   //----------------------------------------
0166   if (ihcal_on)
0167     {
0168       PHG4Prototype2InnerHcalSubsystem *innerhcal = new PHG4Prototype2InnerHcalSubsystem("HCalIn");
0169       innerhcal->set_int_param("hi_eta",1);
0170       innerhcal->set_double_param("place_x",add_place_x);
0171       innerhcal->set_double_param("place_z",144);
0172       innerhcal->SetActive();
0173       innerhcal->SetAbsorberActive();
0174       innerhcal->SetAbsorberTruth(1);
0175       innerhcal->OverlapCheck(true);
0176       innerhcal->SuperDetector("HCALIN");
0177       g4Reco->registerSubsystem(innerhcal);
0178     }
0179   if (ohcal_on)
0180     {
0181       PHG4Prototype2OuterHcalSubsystem *outerhcal = new PHG4Prototype2OuterHcalSubsystem("HCalOut");
0182       outerhcal->set_int_param("hi_eta",1);
0183       outerhcal->set_double_param("place_x",add_place_x);
0184       outerhcal->set_double_param("place_z",229.5);
0185       outerhcal->SetActive();
0186       outerhcal->SetAbsorberActive();
0187       outerhcal->SetAbsorberTruth(1);
0188       outerhcal->OverlapCheck(true);
0189       outerhcal->SuperDetector("HCALOUT");
0190       g4Reco->registerSubsystem(outerhcal);
0191     }
0192   if (cryo_on)
0193     {
0194       double place_z = 175.;
0195       // Cryostat from engineering drawing
0196       PHG4BlockSubsystem *cryo1 = new PHG4BlockSubsystem("cryo1",1);
0197       cryo1->set_double_param("size_x",0.95);
0198       cryo1->set_double_param("size_y",60.96);
0199       cryo1->set_double_param("size_z",60.96);
0200       cryo1->set_double_param("place_x",141.96+0.95/2.+add_place_x);
0201       cryo1->set_double_param("place_z",place_z);
0202       cryo1->set_string_param("material","G4_Al");
0203       cryo1->SetActive(); // it is an active volume - save G4Hits
0204       cryo1->SuperDetector("CRYO");
0205       g4Reco->registerSubsystem(cryo1);
0206 
0207       PHG4BlockSubsystem *cryo2 = new PHG4BlockSubsystem("cryo2",2);
0208       cryo2->set_double_param("size_x",8.89);
0209       cryo2->set_double_param("size_y",60.96);
0210       cryo2->set_double_param("size_z",60.96);
0211       cryo2->set_double_param("place_x",150.72+8.89/2.+add_place_x);
0212       cryo2->set_double_param("place_z",place_z);
0213       cryo2->set_string_param("material","G4_Al");
0214       cryo2->SetActive(); // it is an active volume - save G4Hits
0215       cryo2->SuperDetector("CRYO");
0216       g4Reco->registerSubsystem(cryo2);
0217 
0218       PHG4BlockSubsystem *cryo3 = new PHG4BlockSubsystem("cryo3",3);
0219       cryo3->set_double_param("size_x",2.54);
0220       cryo3->set_double_param("size_y",60.96);
0221       cryo3->set_double_param("size_z",60.96);
0222       cryo3->set_double_param("place_x",173.93+2.54/2.+add_place_x);
0223       cryo3->set_double_param("place_z",place_z);
0224       cryo3->set_string_param("material","G4_Al");
0225       cryo3->SetActive(); // it is an active volume - save G4Hits
0226       cryo3->SuperDetector("CRYO");
0227       g4Reco->registerSubsystem(cryo3);
0228     }
0229   if (bh_on)
0230     { 
0231       // BLACKHOLE, box surrounding the prototype to check for leakage
0232       PHG4BlockSubsystem *bh[5];
0233       // surrounding outer hcal
0234       // top
0235       bh[0] = new PHG4BlockSubsystem("bh1",1);
0236       bh[0]->set_double_param("size_x",270.);
0237       bh[0]->set_double_param("size_y",0.01);
0238       bh[0]->set_double_param("size_z",165.);
0239       bh[0]->set_double_param("place_x",270./2.);
0240       bh[0]->set_double_param("place_y",125./2.);
0241       // bottom
0242       bh[1] = new PHG4BlockSubsystem("bh2",2);
0243       bh[1]->set_double_param("size_x",270.);
0244       bh[1]->set_double_param("size_y",0.01);
0245       bh[1]->set_double_param("size_z",165.);
0246       bh[1]->set_double_param("place_x",270./2.);
0247       bh[1]->set_double_param("place_y",-125./2.);
0248       // right side
0249       bh[2] = new PHG4BlockSubsystem("bh3",3);
0250       bh[2]->set_double_param("size_x",200.);
0251       bh[2]->set_double_param("size_y",125.);
0252       bh[2]->set_double_param("size_z",0.01);
0253       bh[2]->set_double_param("place_x",200./2.);
0254       bh[2]->set_double_param("place_z",165./2.);
0255       // left side
0256       bh[3] = new PHG4BlockSubsystem("bh4",4);
0257       bh[3]->set_double_param("size_x",270.);
0258       bh[3]->set_double_param("size_y",125.);
0259       bh[3]->set_double_param("size_z",0.01);
0260       bh[3]->set_double_param("place_x",270./2.);
0261       bh[3]->set_double_param("place_z",-165./2.);
0262       // back
0263       bh[4] = new PHG4BlockSubsystem("bh5",5);
0264       bh[4]->set_double_param("size_x",0.01);
0265       bh[4]->set_double_param("size_y",125.);
0266       bh[4]->set_double_param("size_z",165.);
0267       bh[4]->set_double_param("place_x",270.);
0268       for (int i=0; i<5; i++)
0269     {
0270       bh[i]->BlackHole();
0271       bh[i]->SetActive();
0272       bh[i]->SuperDetector("BlackHole");
0273       bh[i]->OverlapCheck(true);
0274       g4Reco->registerSubsystem(bh[i]);
0275     }
0276     }
0277   PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
0278   g4Reco->registerSubsystem(truth);
0279 
0280   se->registerSubsystem( g4Reco );
0281   //----------------------------------------
0282   // EMCal digitization
0283   //----------------------------------------
0284   if (cemc_cell)
0285     {
0286       PHG4FullProjSpacalCellReco *cemc_cells = new PHG4FullProjSpacalCellReco("CEMCCYLCELLRECO");
0287       cemc_cells->Detector("CEMC");
0288       cemc_cells->set_timing_window(0.,60.);
0289       cemc_cells->get_light_collection_model().load_data_file(string(getenv("CALIBRATIONROOT")) + string("/CEMC/LightCollection/Prototype2Module.xml"),"data_grid_light_guide_efficiency","data_grid_fiber_trans");
0290 
0291       se->registerSubsystem(cemc_cells);
0292     }
0293   if (cemc_twr)
0294     {
0295       RawTowerBuilder *TowerBuilder = new RawTowerBuilder("EmcRawTowerBuilder");
0296       TowerBuilder->Detector("CEMC");
0297       TowerBuilder->set_sim_tower_node_prefix("SIM");
0298       se->registerSubsystem(TowerBuilder);
0299     }
0300   const double sampling_fraction = 0.0190134; //  +/-   0.000224984  from 0 Degree indenting 32 GeV electron showers
0301   const double photoelectron_per_GeV = 500; //500 photon per total GeV deposition
0302   const double ADC_per_photoelectron_HG = 3.8; // From Sean Stoll, Mar 29
0303   const double ADC_per_photoelectron_LG = 0.24; // From Sean Stoll, Mar 29
0304 
0305   // low gains
0306   if (cemc_digi)
0307     {
0308       RawTowerDigitizer *TowerDigitizer = new RawTowerDigitizer("EmcRawTowerDigitizerLG");
0309       TowerDigitizer->Detector("CEMC");
0310       TowerDigitizer->set_raw_tower_node_prefix("RAW_LG");
0311       TowerDigitizer->set_digi_algorithm(RawTowerDigitizer::kSimple_photon_digitalization);
0312       TowerDigitizer->set_pedstal_central_ADC(0);
0313       TowerDigitizer->set_pedstal_width_ADC(1); // From Jin's guess. No EMCal High Gain data yet! TODO: update
0314       TowerDigitizer->set_photonelec_ADC(1. / ADC_per_photoelectron_LG);
0315       TowerDigitizer->set_photonelec_yield_visible_GeV(photoelectron_per_GeV / sampling_fraction);
0316       TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
0317       se->registerSubsystem(TowerDigitizer);
0318       // high gains
0319       TowerDigitizer = new RawTowerDigitizer("EmcRawTowerDigitizerHG");
0320       TowerDigitizer->Detector("CEMC");
0321       TowerDigitizer->set_raw_tower_node_prefix("RAW_HG");
0322       TowerDigitizer->set_digi_algorithm(
0323                      RawTowerDigitizer::kSimple_photon_digitalization);
0324       TowerDigitizer->set_pedstal_central_ADC(0);
0325       TowerDigitizer->set_pedstal_width_ADC(15); // From John Haggerty, Mar 29
0326       TowerDigitizer->set_photonelec_ADC(1. / ADC_per_photoelectron_HG);
0327       TowerDigitizer->set_photonelec_yield_visible_GeV(photoelectron_per_GeV / sampling_fraction);
0328       TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
0329       se->registerSubsystem(TowerDigitizer);
0330     }
0331   if (cemc_twrcal)
0332     {
0333       RawTowerCalibration *TowerCalibration = new RawTowerCalibration("EmcRawTowerCalibrationLG");
0334       TowerCalibration->Detector("CEMC");
0335       TowerCalibration->set_raw_tower_node_prefix("RAW_LG");
0336       TowerCalibration->set_calib_tower_node_prefix("CALIB_LG");
0337       TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0338       TowerCalibration->set_calib_const_GeV_ADC(1. / ADC_per_photoelectron_LG / photoelectron_per_GeV);
0339       TowerCalibration->set_pedstal_ADC(0);
0340       TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
0341       se->registerSubsystem(TowerCalibration);
0342 
0343 
0344       TowerCalibration = new RawTowerCalibration("EmcRawTowerCalibrationHG");
0345       TowerCalibration->Detector("CEMC");
0346       TowerCalibration->set_raw_tower_node_prefix("RAW_HG");
0347       TowerCalibration->set_calib_tower_node_prefix("CALIB_HG");
0348       TowerCalibration->set_calib_algorithm(
0349                         RawTowerCalibration::kSimple_linear_calibration);
0350       TowerCalibration->set_calib_const_GeV_ADC(1. / ADC_per_photoelectron_HG / photoelectron_per_GeV);
0351       TowerCalibration->set_pedstal_ADC(0);
0352       TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
0353       se->registerSubsystem(TowerCalibration);
0354     }
0355 
0356   //----------------------------------------
0357   // HCal towering
0358   //----------------------------------------
0359   if (ihcal_cell)
0360     {
0361       PHG4Prototype2HcalCellReco *hccell = new PHG4Prototype2HcalCellReco("HCALinCellReco");
0362       hccell->Detector("HCALIN");
0363       se->registerSubsystem(hccell);
0364     }
0365   if (ihcal_twr)
0366     {
0367       Prototype2RawTowerBuilder *hcaltwr = new Prototype2RawTowerBuilder("HCALinRawTowerBuilder");
0368       hcaltwr->Detector("HCALIN");
0369       hcaltwr->set_sim_tower_node_prefix("SIM");
0370       se->registerSubsystem(hcaltwr);
0371     }
0372 
0373 
0374   if (ohcal_cell)
0375     {
0376       PHG4Prototype2HcalCellReco *hccell = new PHG4Prototype2HcalCellReco("HCALoutCellReco");
0377       hccell->Detector("HCALOUT");
0378       se->registerSubsystem(hccell);
0379     }
0380   if (ohcal_twr)
0381     {
0382       Prototype2RawTowerBuilder *hcaltwr = new Prototype2RawTowerBuilder("HCALoutRawTowerBuilder");
0383       hcaltwr->Detector("HCALOUT");
0384       hcaltwr->set_sim_tower_node_prefix("SIM");
0385       se->registerSubsystem(hcaltwr);
0386     }
0387 
0388   //----------------------------------------
0389   // HCal digitization
0390   //----------------------------------------
0391   //  From: Abhisek Sen [mailto:sen.abhisek@gmail.com]
0392   //  Sent: Tuesday, April 19, 2016 10:55 PM
0393   //  To: Huang, Jin <jhuang@bnl.gov>; Haggerty, John <haggerty@bnl.gov>
0394 
0395   //  HCALIN:
0396   //     1/5 pixel / HG ADC channel
0397   //     32/5 pixel / LG ADC channel
0398   //     0.4 MeV/ LG ADC
0399   //     0.4/32 MeV/ HG ADC
0400 
0401   //  HCALOUT:
0402   //     1/5 pixel / HG ADC channel
0403   //     16/5 pixel / LG ADC channel
0404   //     0.2 MeV/ LG ADC
0405   //     0.2/16 MeV/ HG ADC
0406   RawTowerDigitizer *TowerDigitizer = NULL;
0407   if (ihcal_digi)
0408     {
0409       TowerDigitizer = new RawTowerDigitizer("HCALinTowerDigitizerLG");
0410       TowerDigitizer->Detector("HCALIN");
0411       TowerDigitizer->set_raw_tower_node_prefix("RAW_LG");
0412       TowerDigitizer->set_digi_algorithm(RawTowerDigitizer::kSimple_photon_digitalization);
0413       TowerDigitizer->set_pedstal_central_ADC(0);
0414       TowerDigitizer->set_pedstal_width_ADC(1); // From Jin's guess. No EMCal High Gain data yet! TODO: update
0415       TowerDigitizer->set_photonelec_ADC(32. / 5.);
0416       TowerDigitizer->set_photonelec_yield_visible_GeV(32. / 5 / (0.4e-3));
0417       TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
0418       se->registerSubsystem(TowerDigitizer);
0419 
0420       TowerDigitizer = new RawTowerDigitizer("HCALinTowerDigitizerHG");
0421       TowerDigitizer->Detector("HCALIN");
0422       TowerDigitizer->set_raw_tower_node_prefix("RAW_HG");
0423       TowerDigitizer->set_digi_algorithm(RawTowerDigitizer::kSimple_photon_digitalization);
0424       TowerDigitizer->set_pedstal_central_ADC(0);
0425       TowerDigitizer->set_pedstal_width_ADC(1); // From Jin's guess. No EMCal High Gain data yet! TODO: update
0426       TowerDigitizer->set_photonelec_ADC(1. / 5.);
0427       TowerDigitizer->set_photonelec_yield_visible_GeV(1. / 5 / (0.4e-3 / 32));
0428       TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
0429       se->registerSubsystem(TowerDigitizer);
0430     }
0431   if (ohcal_digi)
0432     {
0433       TowerDigitizer = new RawTowerDigitizer("HCALoutTowerDigitizerLG");
0434       TowerDigitizer->Detector("HCALOUT");
0435       TowerDigitizer->set_raw_tower_node_prefix("RAW_LG");
0436       TowerDigitizer->set_digi_algorithm(RawTowerDigitizer::kSimple_photon_digitalization);
0437       TowerDigitizer->set_pedstal_central_ADC(0);
0438       TowerDigitizer->set_pedstal_width_ADC(1); // From Jin's guess. No EMCal High Gain data yet! TODO: update
0439       TowerDigitizer->set_photonelec_ADC(16. / 5.);
0440       TowerDigitizer->set_photonelec_yield_visible_GeV(16. / 5 / (0.2e-3));
0441       TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
0442       se->registerSubsystem(TowerDigitizer);
0443 
0444       TowerDigitizer = new RawTowerDigitizer("HCALoutTowerDigitizerHG");
0445       TowerDigitizer->Detector("HCALOUT");
0446       TowerDigitizer->set_raw_tower_node_prefix("RAW_HG");
0447       TowerDigitizer->set_digi_algorithm(RawTowerDigitizer::kSimple_photon_digitalization);
0448       TowerDigitizer->set_pedstal_central_ADC(0);
0449       TowerDigitizer->set_pedstal_width_ADC(1); // From Jin's guess. No EMCal High Gain data yet! TODO: update
0450       TowerDigitizer->set_photonelec_ADC(1. / 5.);
0451       TowerDigitizer->set_photonelec_yield_visible_GeV(1. / 5 / (0.2e-3 / 16));
0452       TowerDigitizer->set_zero_suppression_ADC(-1000); // no-zero suppression
0453       se->registerSubsystem(TowerDigitizer);
0454     }
0455   //----------------------------------------
0456   // HCal calibration
0457   //----------------------------------------
0458   // 32 GeV Pi+ scan
0459   const double visible_sample_fraction_HCALIN = 7.19505e-02 ; // 1.34152e-02
0460   const double visible_sample_fraction_HCALOUT = 0.0313466 ; //  +/-   0.0067744
0461   RawTowerCalibration *TowerCalibration = NULL;
0462   if (ihcal_twrcal)
0463     {
0464       TowerCalibration = new RawTowerCalibration("HCALinRawTowerCalibrationLG");
0465       TowerCalibration->Detector("HCALIN");
0466       TowerCalibration->set_raw_tower_node_prefix("RAW_LG");
0467       TowerCalibration->set_calib_tower_node_prefix("CALIB_LG");
0468       TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0469       TowerCalibration->set_calib_const_GeV_ADC(0.4e-3 / visible_sample_fraction_HCALIN);
0470       TowerCalibration->set_pedstal_ADC(0);
0471       TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
0472       se->registerSubsystem(TowerCalibration);
0473 
0474       TowerCalibration = new RawTowerCalibration("HCALinRawTowerCalibrationHG");
0475       TowerCalibration->Detector("HCALIN");
0476       TowerCalibration->set_raw_tower_node_prefix("RAW_HG");
0477       TowerCalibration->set_calib_tower_node_prefix("CALIB_HG");
0478       TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0479       TowerCalibration->set_calib_const_GeV_ADC(0.4e-3 / 32 / visible_sample_fraction_HCALIN);
0480       TowerCalibration->set_pedstal_ADC(0);
0481       TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
0482       se->registerSubsystem(TowerCalibration);
0483     }
0484   if (ohcal_twrcal)
0485     {
0486       TowerCalibration = new RawTowerCalibration("HCALoutRawTowerCalibrationLG");
0487       TowerCalibration->Detector("HCALOUT");
0488       TowerCalibration->set_raw_tower_node_prefix("RAW_LG");
0489       TowerCalibration->set_calib_tower_node_prefix("CALIB_LG");
0490       TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0491       TowerCalibration->set_calib_const_GeV_ADC(0.2e-3 / visible_sample_fraction_HCALOUT);
0492       TowerCalibration->set_pedstal_ADC(0);
0493       TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
0494       se->registerSubsystem(TowerCalibration);
0495 
0496       TowerCalibration = new RawTowerCalibration("HCALoutRawTowerCalibrationHG");
0497       TowerCalibration->Detector("HCALOUT");
0498       TowerCalibration->set_raw_tower_node_prefix("RAW_HG");
0499       TowerCalibration->set_calib_tower_node_prefix("CALIB_HG");
0500       TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0501       TowerCalibration->set_calib_const_GeV_ADC(0.2e-3 / 16 / visible_sample_fraction_HCALOUT);
0502       TowerCalibration->set_pedstal_ADC(0);
0503       TowerCalibration->set_zero_suppression_GeV(-1); // no-zero suppression
0504       se->registerSubsystem(TowerCalibration);
0505     }
0506   //----------------------
0507   // QA Histograms
0508   //----------------------
0509   if (cemc_on)
0510     {
0511       se->registerSubsystem(new QAG4SimulationCalorimeter("CEMC",QAG4SimulationCalorimeter::kProcessG4Hit));
0512     }
0513   if (ihcal_on)
0514     {
0515       // TODO: disable QA for HCal right now as there is a hit->particle truth association error at the moment
0516       // se->registerSubsystem(new QAG4SimulationCalorimeter("HCALIN",QAG4SimulationCalorimeter::kProcessG4Hit));
0517     }
0518   if (ohcal_on)
0519     {
0520       // se->registerSubsystem(new QAG4SimulationCalorimeter("HCALOUT",QAG4SimulationCalorimeter::kProcessG4Hit));
0521     }
0522   //----------------------
0523   // G4HitNtuple
0524   //----------------------
0525   if (hit_ntuple)
0526     {
0527       G4HitNtuple *hit = new G4HitNtuple("G4HitNtuple","g4hitntuple.root");
0528       hit->AddNode("HCALIN", 0);
0529       hit->AddNode("HCALOUT", 1);
0530       hit->AddNode("CRYO", 2);
0531       hit->AddNode("BlackHole", 3);
0532       hit->AddNode("ABSORBER_HCALIN", 10);
0533       hit->AddNode("ABSORBER_HCALOUT", 11);
0534       se->registerSubsystem(hit);
0535     }
0536   // G4ScintillatorSlatTTree *scintcell = new G4ScintillatorSlatTTree("inslat");
0537   // scintcell->Detector("HCALIN");
0538   // se->registerSubsystem(scintcell);
0539 
0540   // scintcell = new G4ScintillatorSlatTTree("outslat");
0541   // scintcell->Detector("HCALOUT");
0542   // se->registerSubsystem(scintcell);
0543 
0544 
0545   //----------------------
0546   // save a comprehensive  evaluation file
0547   //----------------------
0548   if (dstreader)
0549     {
0550       PHG4DSTReader* ana = new PHG4DSTReader(string("DSTReader.root"));
0551       ana->set_save_particle(true);
0552       ana->set_load_all_particle(false);
0553       ana->set_load_active_particle(false);
0554       ana->set_save_vertex(true);
0555       ana->set_tower_zero_sup(-1000); // no zero suppression
0556 
0557       //  ana->AddNode("CEMC");
0558       //  if (absorberactive)
0559       //    {
0560       //      ana->AddNode("ABSORBER_CEMC");
0561       //    }
0562 
0563       if (cemc_twr)
0564         ana->AddTower("SIM_CEMC");
0565       if (cemc_digi)
0566         ana->AddTower("RAW_LG_CEMC");
0567       if (cemc_twrcal)
0568         ana->AddTower("CALIB_LG_CEMC"); // Low gain CEMC
0569       if (cemc_digi)
0570         ana->AddTower("RAW_HG_CEMC");
0571       if (cemc_twrcal)
0572         ana->AddTower("CALIB_HG_CEMC"); // High gain CEMC
0573 
0574       if (ohcal_twr)
0575         ana->AddTower("SIM_HCALOUT");
0576       if (ihcal_twr)
0577         ana->AddTower("SIM_HCALIN");
0578 
0579       if (ihcal_digi)
0580         ana->AddTower("RAW_LG_HCALIN");
0581       if (ihcal_digi)
0582         ana->AddTower("RAW_HG_HCALIN");
0583       if (ohcal_digi)
0584         ana->AddTower("RAW_LG_HCALOUT");
0585       if (ohcal_digi)
0586         ana->AddTower("RAW_HG_HCALOUT");
0587 
0588       if (ihcal_twrcal)
0589         ana->AddTower("CALIB_LG_HCALIN");
0590       if (ihcal_twrcal)
0591         ana->AddTower("CALIB_HG_HCALIN");
0592       if (ohcal_twrcal)
0593         ana->AddTower("CALIB_LG_HCALOUT");
0594       if (ohcal_twrcal)
0595         ana->AddTower("CALIB_HG_HCALOUT");
0596 
0597       if (bh_on)
0598         ana->AddNode("BlackHole"); // add a G4Hit node
0599 
0600       se->registerSubsystem(ana);
0601     }
0602 
0603   // Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT","/phenix/scratch/pinkenbu/G4Prototype2Hcalin.root");
0604   // out->AddNode("G4RootScintillatorSlat_HCALIN");
0605   // se->registerOutputManager(out);
0606 
0607   // out = new Fun4AllDstOutputManager("DSTHCOUT","/phenix/scratch/pinkenbu/G4Prototype2Hcalout.root");
0608   // out->AddNode("G4RootScintillatorSlat_HCALOUT");
0609   // se->registerOutputManager(out);
0610 
0611   if (dstoutput)
0612     {
0613       Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT","G4Prototype3New.root");
0614       se->registerOutputManager(out);
0615     }
0616 
0617   Fun4AllInputManager *in = new Fun4AllDummyInputManager( "JADE");
0618   se->registerInputManager( in );
0619   if (nEvents <= 0)
0620     {
0621       return;
0622     }
0623   se->run(nEvents);
0624 
0625   se->End();
0626 
0627   QAHistManagerDef::saveQARootFile("G4Prototype2_qa.root");
0628 
0629 
0630   //   std::cout << "All done" << std::endl;
0631   delete se;
0632   gSystem->Exit(0);
0633 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!