EIC-IR/macros_fun4all/Fun4All_G4_EICIR.C

0001
0002 int Fun4All_G4_EICIR(
0003              const int nEvents = 1,
0004              const char * inputFile = "/gpfs02/phenix/prod/sPHENIX/preCDR/pro.1-beta.5/single_particle/spacal1d/fieldmap/G4Hits_sPHENIX_e-_eta0_16GeV.root",
0005              const char * outputFile = "G4EICIR.root"
0006              )
0007 {
0008   //===============
0009   // Input options
0010   //===============
0011
0012   // Use particle gun
0013   const bool pgun = true;
0014
0015   //======================
0016   // What to run
0017   //======================
0018
0019   bool do_pipe = false;
0020
0021   bool do_magnet = true;
0022
0023   // Extended IR
0024   bool do_ExtendedIR = true;
0025
0026   //---------------
0027   // Load libraries
0028   //---------------
0029
0030   gSystem->Load("libfun4all.so");
0031   gSystem->Load("libg4detectors.so");
0032   gSystem->Load("libphhepmc.so");
0033   gSystem->Load("libg4testbench.so");
0034   gSystem->Load("libg4hough.so");
0035   gSystem->Load("libg4eval.so");
0036
0037   // establish the geometry and reconstruction setup
0038   gROOT->LoadMacro("G4Setup_EICIR.C");
0039   G4Init(do_magnet,do_pipe,do_ExtendedIR);
0040
0041   int absorberactive = 0; // set to 1 to make all absorbers active volumes
0042   //  const string magfield = "1.5"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
0043   const string magfield = "/phenix/upgrades/decadal/fieldmaps/sPHENIX.2d.root"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
0044   const float magfield_rescale = 1.4/1.5; // scale the map to a 1.4 T field
0045
0046   //---------------
0047   // Fun4All server
0048   //---------------
0049
0050   Fun4AllServer *se = Fun4AllServer::instance();
0051   //se->Verbosity(100); // uncomment for batch production running with minimal output messages
0052   se->Verbosity(Fun4AllServer::VERBOSITY_SOME); // uncomment for some info for interactive running
0053   // just if we set some flags somewhere in this macro
0054   recoConsts *rc = recoConsts::instance();
0055   // By default every random number generator uses
0056   // PHRandomSeed() which reads /dev/urandom to get its seed
0057   // if the RANDOMSEED flag is set its value is taken as seed
0058   // You can either set this to a random value using PHRandomSeed()
0059   // which will make all seeds identical (not sure what the point of
0060   // this would be:
0061   //  rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
0062   // or set it to a fixed value so you can debug your code
0063   // rc->set_IntFlag("RANDOMSEED", 12345);
0064
0065   //-----------------
0066   // Event generation
0067   //-----------------
0068
0069   if (pgun)
0070     {
0071       /* angle of particle phi:
0072      pz = p * cos(psi)
0073      px = p * sin(psi) */
0074       double psi_mrad = 0;
0075
0076       double ptot = 250*1;
0077
0078       double vx = 0;
0079       double vy = 0;
0080       double vz = 0;
0081
0082       double px = ptot * sin(psi_mrad / 1000.);
0083       double py = 0;
0084       double pz = ptot * cos(psi_mrad / 1000.);
0085
0086       PHG4ParticleGun*gun = new PHG4ParticleGun();
0087       gun->set_name("proton");
0088       gun->set_vtx(vx,vy,vz);
0089       gun->set_mom(px,py,pz);
0090       se->registerSubsystem(gun);
0091     }
0092   else
0093     {
0094       cout << "WARNING: No events being generated!" << endl;
0095     }
0096
0097   //---------------------
0098   // Detector description
0099   //---------------------
0100
0101   G4Setup(absorberactive, magfield, TPythia6Decayer::kAll,do_magnet,do_pipe,do_ExtendedIR,magfield_rescale);
0102
0103
0104
0105   //--------------
0106   // IO management
0107   //--------------
0108
0109   // for single particle generators we just need something which drives
0110   // the event loop, the Dummy Input Mgr does just that
0111   Fun4AllInputManager *in = new Fun4AllDummyInputManager( "JADE");
0112   se->registerInputManager( in );
0113
0114   //Convert DST to human command readable TTree for quick poke around the outputs
0115   gROOT->LoadMacro("G4_DSTReader_EICIR.C");
0116
0117   G4DSTreader_EICIR( outputFile, //
0118              /*int*/ absorberactive );
0119
0120   Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputFile);
0121   //if (do_dst_compress) DstCompress(out);
0122   se->registerOutputManager(out);
0123
0124   if (nEvents == 0 && !readhits && !readhepmc)
0125     {
0126       cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
0127       cout << "it will run forever, so I just return without running anything" << endl;
0128       return;
0129     }
0130
0131   if (nEvents < 0)
0132     {
0133       PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco("PHG4RECO");
0134       g4->ApplyCommand("/control/execute vis.mac");
0135       //g4->StartGui();
0136       se->run(1);
0137
0138       se->End();
0139       std::cout << "All done" << std::endl;
0140
0141       std::cout << "==== Useful display commands ==" << std::endl;
0142       cout << "draw axis: " << endl;
0143       cout << " G4Cmd(\"/vis/scene/add/axes 0 0 0 50 cm\")" << endl;
0144       cout << "zoom" << endl;
0145       cout << " G4Cmd(\"/vis/viewer/zoom 1\")" << endl;
0146       cout << "viewpoint:" << endl;
0147       cout << " G4Cmd(\"/vis/viewer/set/viewpointThetaPhi 0 0\")" << endl;
0148       cout << "panTo:" << endl;
0149       cout << " G4Cmd(\"/vis/viewer/panTo 0 0 cm\")" << endl;
0150       cout << "print to eps:" << endl;
0151       cout << " G4Cmd(\"/vis/ogl/printEPS\")" << endl;
0152       cout << "set background color:" << endl;
0153       cout << " G4Cmd(\"/vis/viewer/set/background white\")" << endl;
0154       std::cout << "===============================" << std::endl;
0155     }
0156   else
0157     {
0158       se->run(nEvents);
0159
0160       se->End();
0161       std::cout << "All done" << std::endl;
0162       delete se;
0163       gSystem->Exit(0);
0164     }
0165
0166 }
0167
0168
0169 void
0170 G4Cmd(const char * cmd)
0171 {
0172   Fun4AllServer *se = Fun4AllServer::instance();
0173   PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco("PHG4RECO");
0174   g4->ApplyCommand(cmd);
0175 }