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