analysis/sPhenixAj/Wrapper.cxx

0001 int Min_cemc_layer = -1;
0002 int Max_cemc_layer = -1;
0003 int Min_hcal_in_layer = -1;
0004 int Max_hcal_in_layer = -1;
0005 int Min_hcal_out_layer = -1;
0006 int Max_hcal_out_layer = -1;
0007 int Min_preshower_layer = -1;
0008 int Max_preshower_layer = -1;
0009 int Min_si_layer = -1;
0010 int Max_si_layer = -1;
0011 int Cemc_slats_per_cell = 72; // make it 2*2*2*3*3 so we can try other combinations
0012
0013 // #include<src/PJTranslator.h>
0014
0015 int Wrapper(
0016         int nEvents = 1,
0017         const char * fileList = "/star/u/kkauder/sphenix/analysis/sPhenixAj/SingleEventPythia.list",
0018         // const char * fileList = "/star/u/kkauder/sphenix/analysis/sPhenixAj/singleparticle.list",
0019         const char * outputFile = "G4sPHENIXCells.root"
0020         )
0021 {
0022   //===============
0023   // Input options
0024   //===============
0025
0026   // Either:
0027   // read previously generated g4-hits files, in this case it opens a DST and skips
0028   // the simulations step completely. The G4Setup macro is only loaded to get information
0029   // about the number of layers used for the cell reco code
0030   const bool readhits = false;
0031   // Or:
0032   // read files in HepMC format (typically output from event generators like hijing or pythia)
0033   const bool readhepmc = true; // read HepMC files
0034   // Or:
0035   // Use particle generator
0036
0037   //======================
0038   // What to run
0039   //======================
0040
0041   bool do_svtx = true;
0042   bool do_svtx_cell = true;
0043   bool do_svtx_track = true;
0044   bool do_svtx_eval = true;
0045
0046   bool do_preshower = false;
0047
0048   bool do_cemc = true;
0049   bool do_cemc_cell = true;
0050   bool do_cemc_twr = true;
0051   bool do_cemc_cluster = true;
0052   bool do_cemc_eval = true;
0053
0054   bool do_hcalin = true;
0055   bool do_hcalin_cell = true;
0056   bool do_hcalin_twr = true;
0057   bool do_hcalin_cluster = true;
0058   bool do_hcalin_eval = true;
0059
0060   bool do_magnet = true;
0061
0062   bool do_hcalout = true;
0063   bool do_hcalout_cell = true;
0064   bool do_hcalout_twr = true;
0065   bool do_hcalout_cluster = true;
0066   bool do_hcalout_eval = true;
0067
0068   //---------------
0069   // Load libraries
0070   //---------------
0071
0072   gSystem->Load("libfun4all.so");
0073   gSystem->Load("libg4detectors.so");
0074   gSystem->Load("libphhepmc.so");
0075   gSystem->Load("libg4testbench.so");
0076   gSystem->Load("libg4hough.so");
0077   gSystem->Load("libg4calo.so");
0078   gSystem->Load("libg4eval.so");
0079   gSystem->Load("libPHAJMaker.so");
0080
0081   // establish the geometry and reconstruction setup
0082   gROOT->LoadMacro("G4Setup_sPHENIX.C");
0083   G4Init(do_svtx,do_preshower,do_cemc,do_hcalin,do_magnet,do_hcalout);
0084
0085   int absorberactive = 1; // set to 1 to make all absorbers active volumes
0086   //  const string magfield = "1.5"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
0087   const string magfield = "/phenix/upgrades/decadal/fieldmaps/bPHENIX.dp.root"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
0088
0089
0090   //---------------
0091   // Fun4All server
0092   //---------------
0093   Fun4AllServer *se = Fun4AllServer::instance();
0094   se->Verbosity(0);
0095
0096   // unique seed
0097   int uniqueseed = TRandom3(0).GetSeed();
0098   recoConsts *rc = recoConsts::instance();
0099   rc->set_IntFlag("RANDOMSEED", uniqueseed);
0100
0101   //-----------------
0102   // Event generation
0103   //-----------------
0104
0105   if (readhits)
0106     {
0107       // Get the hits from a file
0108       // The input manager is declared later
0109     }
0110   else if (readhepmc)
0111     {
0112       // this module is needed to read the HepMC records into our G4 sims
0113       // but only if you read HepMC input files
0114       HepMCNodeReader *hr = new HepMCNodeReader();
0115       se->registerSubsystem(hr);
0116     }
0117   else
0118     {
0119       // toss low multiplicity dummy events
0120       PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
0121       gen->add_particles("e-",5); // mu+,e+,proton,pi+,Upsilon
0122       gen->add_particles("e+",5); // mu-,e-,anti_proton,pi-
0123       if (readhepmc) {
0124     gen->set_reuse_existing_vertex(true);
0125     gen->set_existing_vertex_offset_vector(0.0,0.0,0.0);
0126       } else {
0127     gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
0128                           PHG4SimpleEventGenerator::Uniform,
0129                           PHG4SimpleEventGenerator::Uniform);
0130     gen->set_vertex_distribution_mean(0.0,0.0,0.0);
0131     gen->set_vertex_distribution_width(0.0,0.0,5.0);
0132       }
0133       gen->set_vertex_size_function(PHG4SimpleEventGenerator::Uniform);
0134       gen->set_vertex_size_parameters(0.0,0.0);
0135       gen->set_eta_range(-0.5, 0.5);
0136       gen->set_phi_range(-1.0*TMath::Pi(), 1.0*TMath::Pi());
0137       gen->set_pt_range(0.1, 10.0);
0138       gen->set_embedflag(1);
0139       gen->set_seed(uniqueseed);
0140       gen->set_verbosity(0);
0141       se->registerSubsystem(gen);
0142     }
0143
0144   if (!readhits)
0145     {
0146       //---------------------
0147       // Detector description
0148       //---------------------
0149
0150       G4Setup(absorberactive, magfield, TPythia6Decayer::kAll,
0151           do_svtx, do_preshower, do_cemc, do_hcalin, do_magnet, do_hcalout);
0152     }
0153
0154   //------------------
0155   // Detector Division
0156   //------------------
0157
0158
0159   if (do_svtx_cell) Svtx_Cells();
0160
0161   if (do_cemc_cell) CEMC_Cells();
0162
0163   if (do_hcalin_cell) HCALInner_Cells();
0164
0165   if (do_hcalout_cell) HCALOuter_Cells();
0166
0167   //-----------------------------
0168   // CEMC towering and clustering
0169   //-----------------------------
0170
0171   if (do_cemc_twr) CEMC_Towers();
0172   if (do_cemc_cluster) CEMC_Clusters();
0173
0174   //-----------------------------
0175   // HCAL towering and clustering
0176   //-----------------------------
0177
0178   if (do_hcalin_cell) HCALInner_Towers();
0179   if (do_hcalin_cell) HCALInner_Clusters();
0180
0181   if (do_hcalout_cell) HCALOuter_Towers();
0182   if (do_hcalout_cell) HCALOuter_Clusters();
0183
0184   //--------------
0185   // SVTX tracking
0186   //--------------
0187
0188   if (do_svtx_track) Svtx_Reco();
0189
0190   //----------------------
0191   // Simulation evaluation
0192   //----------------------
0193
0194   if (do_svtx_eval) Svtx_Eval("g4svtx_eval.root");
0195
0196   if (do_cemc_eval) CEMC_Eval("g4cemc_eval.root");
0197
0198   if (do_hcalin_cell) HCALInner_Eval("g4hcalin_eval.root");
0199
0200   if (do_hcalout_cell) HCALOuter_Eval("g4hcalout_eval.root");
0201
0202   // PHAJMaker *ana = new PHAJMaker(outputFile);
0203   // se->registerSubsystem(ana);
0204
0205   PJTranslator *trans = new PJTranslator( "PJTranslator" );
0206   se->registerSubsystem(trans);
0207
0208   // PHAJMaker *ana = new PHAJMaker(outputFile);
0209   PHAJMaker *AJ = new PHAJMaker("AJ");
0210   se->registerSubsystem(AJ);
0211
0212   //--------------
0213   // IO management
0214   //--------------
0215
0216   if (readhits)
0217     {
0218       // Hits file
0219       Fun4AllInputManager *hitsin = new Fun4AllDstInputManager("DSTin");
0220       throw (-1);
0221       //      hitsin->fileopen(inputFile);
0222       // hitsin->AddListFile(fileList);
0223       se->registerInputManager(hitsin);
0224     }
0225   if (readhepmc)
0226     {
0227       Fun4AllInputManager *in = new Fun4AllHepMCInputManager( "DSTIN");
0228       in->AddListFile(fileList);
0229       se->registerInputManager( in );
0230       //se->fileopen( in->Name(), inputFile );
0231     }
0232   else
0233     {
0234       // for single particle generators we just need something which drives
0235       // the event loop, the Dummy Input Mgr does just that
0236       throw(-1);
0237       Fun4AllInputManager *in = new Fun4AllDummyInputManager( "JADE");
0238       se->registerInputManager( in );
0239     }
0240
0241
0242   // Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputFile);
0243   // se->registerOutputManager(out);
0244
0245   //-----------------
0246   // Event processing
0247   //-----------------
0248   se->run(nEvents);
0249
0250   //-----
0251   // Exit
0252   //-----
0253
0254   se->End();
0255   se->dumpHistos();
0256   std::cout << "All done" << std::endl;
0257   delete se;
0258   gSystem->Exit(0);
0259 }
0260
0261 int Get_Min_cemc_layer()
0262 {
0263   return Min_cemc_layer;
0264 }
0265
0266 int Get_Max_cemc_layer()
0267 {
0268   return Max_cemc_layer;
0269 }
0270
0271 int Get_Min_hcal_in_layer()
0272 {
0273   return Min_hcal_in_layer;
0274 }
0275
0276 int Get_Max_hcal_in_layer()
0277 {
0278   return Max_hcal_in_layer;
0279 }
0280
0281 int Get_Min_hcal_out_layer()
0282 {
0283   return Min_hcal_out_layer;
0284 }
0285
0286 int Get_Max_hcal_out_layer()
0287 {
0288   return Max_hcal_out_layer;
0289 }
0290
0291 int Get_Min_si_layer()
0292 {
0293   return Min_si_layer;
0294 }
0295
0296 int Get_Max_si_layer()
0297 {
0298   return Max_si_layer;
0299 }
0300
0301 int Get_Min_preshower_layer()
0302 {
0303   return Min_preshower_layer;
0304 }
0305
0306 int Get_Max_preshower_layer()
0307 {
0308   return Max_preshower_layer;
0309 }