Test/W-Shashlik/Fun4All_G4_W.C

0001 int
0002 Fun4All_G4_W(const int nEvents = 100, const double momentum = 5, const char *outfile = "test.root")
0003 {
0004
0005   gSystem->Load("libfun4all");
0006   gSystem->Load("libg4detectors");
0007   gSystem->Load("libg4testbench");
0008   gSystem->Load("libg4eval");
0009   gSystem->Load("libg4histos");
0010
0011   ///////////////////////////////////////////
0012   // Make the Server
0013   //////////////////////////////////////////
0014   Fun4AllServer *se = Fun4AllServer::instance();
0015   se->Verbosity(1);
0016
0017 //  gSystem->Load("libPHPythia8.so");
0018 //
0019 //  PHPythia8* pythia8 = new PHPythia8();
0020 //  // see coresoftware/generators/PHPythia8 for example config
0021 //  pythia8->set_config_file("./phpythia8_510.cfg");
0022 //  se->registerSubsystem(pythia8);
0023 //
0024 //  HepMCNodeReader *hr = new HepMCNodeReader();
0025 //  se->registerSubsystem(hr);
0026
0027 //  // particle gun
0028 //  PHG4ParticleGun *gun = new PHG4ParticleGun("PGUN");
0029 //  //  gun->set_name("anti_proton");
0030 //  gun->set_name("mu-");
0031 //  //  gun->set_name("proton");
0032 //  gun->set_vtx(-0,0,0);
0033 //  gun->set_mom(10,0,0.);
0034 //  se->registerSubsystem(gun);
0035
0036   // Fun4All G4 module
0037   PHG4Reco* g4Reco = new PHG4Reco();
0038   // no magnetic field
0039   g4Reco->set_field(0);
0040   g4Reco->set_field_rescale(1);
0041   // size of the world - every detector has to fit in here
0042   g4Reco->SetWorldSizeX(500);
0043   g4Reco->SetWorldSizeY(500);
0044   g4Reco->SetWorldSizeZ(500);
0045   // shape of our world - it is a box
0046   g4Reco->SetWorldShape("G4BOX");
0047   // this is what our world is filled with
0048   g4Reco->SetWorldMaterial("G4_AIR");
0049   // Geant4 Physics list to use
0050   g4Reco->SetPhysicsList("QGSP_BERT");
0051
0052   ///////////////////////////////////////////
0053   // Event generator
0054   //////////////////////////////////////////
0055   // toss low multiplicity dummy events
0056   PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
0057   gen->add_particles("e-", 1); // mu+,e+,proton,pi+,Upsilon
0058   // gen->add_particles("e+",5); // mu-,e-,anti_proton,pi-
0059
0060   gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
0061       PHG4SimpleEventGenerator::Uniform, PHG4SimpleEventGenerator::Uniform);
0062   gen->set_vertex_distribution_mean(0.0, 0.0, 0.0);
0063   gen->set_vertex_distribution_width(0.0, 0.0, 5.0);
0064   gen->set_vertex_size_function(PHG4SimpleEventGenerator::Uniform);
0065   gen->set_vertex_size_parameters(0.0, 0.0);
0066   gen->set_eta_range(-0.1, 0.1);
0067   gen->set_phi_range(-1.0 * TMath::Pi(), 1.0 * TMath::Pi());
0068   gen->set_p_range(momentum, momentum);
0069   gen->Embed(1);
0070   gen->Verbosity(0);
0071   se->registerSubsystem(gen);
0072
0073   ///////////////////////////////////////////
0074   // W-shashlyk
0075   //////////////////////////////////////////
0076   // radiation length/layer = 0.12/0.3504 + 0.03/1.436 + 0.15/41.31 = 0.37
0077   // 20 Radiation length = 54 layers
0078   // Total length = 16 cm
0079   double scintiwidth = 0.15;
0080   double tungstenwidth = 0.12;
0081   double copperwidth = 0.015;
0082   double no_overlapp = 0.0001;
0083   double radius = 95;
0084   int Max_cemc_layer = 54;
0085   int absorberactive = 1;
0086   int overlapcheck = 1;
0087
0088   PHG4CylinderSubsystem *cemc;
0089
0090   for (int ilayer = 0; ilayer <= Max_cemc_layer; ilayer++)
0091     {
0092
0093       cemc = new PHG4CylinderSubsystem("ABSORBER_CEMC", 3 * ilayer + 0);
0094       cemc->set_double_param("radius", radius);
0095       cemc->set_string_param("material", "G4_Cu");
0096       cemc->set_double_param("thickness", copperwidth);
0097       if (absorberactive)
0098         cemc->SetActive();
0099       cemc->OverlapCheck(overlapcheck);
0100       g4Reco->registerSubsystem(cemc);
0101       cemc->SuperDetector("ABSORBER_CEMC");
0102
0103       radius += copperwidth;
0104       radius += no_overlapp;
0105
0106       cemc = new PHG4CylinderSubsystem("ABSORBER_CEMC", 3 * ilayer + 1);
0107       cemc->set_double_param("radius", radius);
0108       cemc->set_string_param("material", "G4_W");
0109       cemc->set_double_param("thickness", tungstenwidth);
0110       if (absorberactive)
0111         cemc->SetActive();
0112       cemc->OverlapCheck(overlapcheck);
0113       g4Reco->registerSubsystem(cemc);
0114       cemc->SuperDetector("ABSORBER_CEMC");
0115
0116       radius += tungstenwidth;
0117       radius += no_overlapp;
0118
0119       cemc = new PHG4CylinderSubsystem("ABSORBER_CEMC", 3 * ilayer + 2);
0120       cemc->set_double_param("radius", radius);
0121       cemc->set_string_param("material", "G4_Cu");
0122       cemc->set_double_param("thickness", copperwidth);
0123       if (absorberactive)
0124         cemc->SetActive();
0125       cemc->OverlapCheck(overlapcheck);
0126       g4Reco->registerSubsystem(cemc);
0127       cemc->SuperDetector("ABSORBER_CEMC");
0128
0129       radius += copperwidth;
0130       radius += no_overlapp;
0131
0132       cemc = new PHG4CylinderSubsystem("CEMC", 3 * ilayer + 0);
0133       cemc->set_double_param("radius", radius);
0134       cemc->set_string_param("material", "G4_POLYSTYRENE");
0135       cemc->set_double_param("thickness", scintiwidth);
0136       if (absorberactive)
0137         cemc->SetActive();
0138       cemc->OverlapCheck(overlapcheck);
0139       g4Reco->registerSubsystem(cemc);
0140       cemc->SuperDetector("CEMC");
0141
0142       radius += scintiwidth;
0143       radius += no_overlapp;
0144     }
0145
0146   //----------------------------------------
0147   // BLACKHOLE
0148
0149   // swallow all particles coming out of the backend of sPHENIX
0150   PHG4CylinderSubsystem *blackhole = new PHG4CylinderSubsystem("BH", 1);
0151   blackhole->set_double_param("radius", radius + 10); // add 10 cm
0152
0153   blackhole->set_int_param("lengthviarapidity", 0);
0154   blackhole->set_double_param("length", g4Reco->GetWorldSizeZ() - no_overlapp); // make it cover the world in length
0155   blackhole->BlackHole();
0156   blackhole->set_double_param("thickness", 0.1); // it needs some thickness
0157   blackhole->SetActive(); // always see what leaks out
0158   blackhole->OverlapCheck(overlapcheck);
0159   g4Reco->registerSubsystem(blackhole);
0160
0161   //----------------------------------------
0162   // FORWARD BLACKHOLEs
0163   // +Z
0164   blackhole = new PHG4CylinderSubsystem("BH_FORWARD_PLUS", 1);
0165   blackhole->SuperDetector("BH_FORWARD_PLUS");
0166   blackhole->set_double_param("radius", 0); // add 10 cm
0167   blackhole->set_int_param("lengthviarapidity", 0);
0168   blackhole->set_double_param("length", 0.1); // make it cover the world in length
0169   blackhole->set_double_param("place_z",
0170       g4Reco->GetWorldSizeZ() / 2. - 0.1 - no_overlapp);
0171   blackhole->BlackHole();
0172   blackhole->set_double_param("thickness", radius - no_overlapp); // it needs some thickness
0173   blackhole->SetActive(); // always see what leaks out
0174   blackhole->OverlapCheck(overlapcheck);
0175   g4Reco->registerSubsystem(blackhole);
0176
0177   blackhole = new PHG4CylinderSubsystem("BH_FORWARD_NEG", 1);
0178   blackhole->SuperDetector("BH_FORWARD_NEG");
0179   blackhole->set_double_param("radius", 0); // add 10 cm
0180   blackhole->set_int_param("lengthviarapidity", 0);
0181   blackhole->set_double_param("length", 0.1); // make it cover the world in length
0182   blackhole->set_double_param("place_z",
0183       -g4Reco->GetWorldSizeZ() / 2. + 0.1 + no_overlapp);
0184   blackhole->BlackHole();
0185   blackhole->set_double_param("thickness", radius - no_overlapp); // it needs some thickness
0186   blackhole->SetActive(); // always see what leaks out
0187   blackhole->OverlapCheck(overlapcheck);
0188   g4Reco->registerSubsystem(blackhole);
0189
0190   //----------------------------------------
0191   // PHG4TruthSubsystem
0192
0193   PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
0194   g4Reco->registerSubsystem(truth);
0195
0196   se->registerSubsystem(g4Reco);
0197
0198   ///////////////////////////////////////////
0199   // Output
0200   ///////////////////////////////////////////
0201
0202 //  if (outfile)
0203 //    {
0204 //      Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT",outfile);
0205 //      se->registerOutputManager(out);
0206 //
0207 //    }
0208   if (outfile)
0209     {
0210       // save a comprehensive  evaluation file
0211       PHG4DSTReader* ana = new PHG4DSTReader(
0212           string(outfile) + string("_DSTReader.root"));
0213       ana->set_save_particle(true);
0214       ana->set_load_all_particle(false);
0215       ana->set_load_active_particle(false);
0216       ana->set_save_vertex(true);
0217       if (nEvents > 0 && nEvents < 2)
0218         {
0219           ana->Verbosity(2);
0220         }
0221       ana->AddNode("CEMC");
0222       se->registerSubsystem(ana);
0223
0224       gSystem->Load("libqa_modules");
0225
0226       QAG4SimulationCalorimeter * qa = new QAG4SimulationCalorimeter("CEMC",
0227           QAG4SimulationCalorimeter::kProcessG4Hit);
0228       se->registerSubsystem(qa);
0229
0230     }
0231
0232   // input - we need a dummy to drive the event loop
0233   Fun4AllInputManager *in = new Fun4AllDummyInputManager("JADE");
0234   se->registerInputManager(in);
0235
0236   //-----------------
0237   // Event processing
0238   //-----------------
0239   if (nEvents < 0)
0240     {
0241       return;
0242     }
0243   // if we run the particle generator and use 0 it'll run forever
0244   if (nEvents == 0)
0245     {
0246       cout
0247           << "using 0 for number of events is a bad idea when using particle generators"
0248           << endl;
0249       cout << "it will run forever, so I just return without running anything"
0250           << endl;
0251       return;
0252     }
0253
0254   se->run(nEvents);
0255
0256 //  if (!readhits)
0257 //    {
0258 //      // Geometry export:
0259 //      PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco("PHG4RECO");
0260 //      g4->Dump_GDML("sPHENIX.gdml");
0261 //    }
0262   //-----
0263   // Exit
0264   //-----
0265   if (outfile)
0266   {
0267     gSystem->Load("libqa_modules");
0268     QAHistManagerDef::saveQARootFile(string(outfile) + "_qa.root");
0269
0270   }
0271
0272   se->End();
0273   std::cout << "All done" << std::endl;
0274   delete se;
0275   gSystem->Exit(0);
0276
0277   return;
0278
0279 }
0280
0281 void
0282 G4Cmd(const char * cmd)
0283 {
0284   Fun4AllServer *se = Fun4AllServer::instance();
0285   PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco("PHG4RECO");
0286   g4->ApplyCommand(cmd);
0287 }
0288
0289 PHG4ParticleGun *
0290 get_gun(const char *name = "PGUN")
0291 {
0292   Fun4AllServer *se = Fun4AllServer::instance();
0293   PHG4ParticleGun *pgun = (PHG4ParticleGun *) se->getSubsysReco(name);
0294   return pgun;
0295 }
0296