ParticleGen-Jingyu/macro/Fun4All_physiTuto.C

0001 // Fun4All headers
0002 #include <fun4all/Fun4AllServer.h>
0003 #include <fun4all/Fun4AllDstInputManager.h>
0004
0005 #include <GlobalVariables.C>
0006
0007 //#include <G4Setup_sPHENIX_Bbc.C>
0008 #include <G4Setup_sPHENIX_macro.C>
0009 // #include <G4_Bbc.C>
0010 // #include <G4_CaloTrigger.C>
0011 // #include <G4_Centrality.C>
0012 // #include <G4_DSTReader.C>
0013 // #include <G4_Global.C>
0014 // #include <G4_HIJetReco.C>
0015 #include <G4_Input.C>
0016 // #include <G4_Jets.C>
0017 // #include <G4_KFParticle.C>
0018 // #include <G4_ParticleFlow.C>
0019 // #include <G4_Production.C>
0020 // #include <G4_TopoClusterReco.C>
0021
0022 #include <Trkr_RecoInit.C>
0023 #include <Trkr_Clustering.C>
0024 #include <Trkr_LaserClustering.C>
0025 #include <Trkr_Reco.C>
0026 #include <Trkr_Eval.C>
0027 #include <G4_ActsGeom.C>
0028 // #include <Trkr_QA.C>
0029
0030 // #include <Trkr_Diagnostics.C>
0031 // #include <G4_User.C>
0032 // #include <QA.C>
0033
0034 #include <ffamodules/FlagHandler.h>
0035 #include <ffamodules/HeadReco.h>
0036 #include <ffamodules/SyncReco.h>
0037 #include <ffamodules/CDBInterface.h>
0038
0039 #include <fun4all/Fun4AllDstOutputManager.h>
0040 #include <fun4all/Fun4AllOutputManager.h>
0041 #include <fun4all/Fun4AllServer.h>
0042
0043 #include <phool/PHRandomSeed.h>
0044 #include <phool/recoConsts.h>
0045
0046 #include <iostream>
0047 #include <typeinfo>
0048
0049 R__LOAD_LIBRARY(libfun4all.so)
0050
0051 // #include <physiTuto/tutorial.h>
0052 #include <tutorial/tutorial.h>
0053 #include <tutorial/TruthToSvtxTrack.h>
0054 R__LOAD_LIBRARY( libtutorial.so )
0055
0056 #include <caloreco/CaloGeomMappingv2.h>
0057 #include <calogeomtest/CaloGeomTest.h>
0058 R__LOAD_LIBRARY(libcalo_reco.so)
0059 R__LOAD_LIBRARY(libcalogeomtest.so)
0060
0061 void Fun4All_physiTuto(
0062     int nEvents = 3,
0063     const double particle_pT = 5.0,
0064     const double energy_range_up = 0.0,
0065     const double energy_range_down = 0.0,
0066     const string particle_species = "Electron",
0067     const string &output_directory = "/sphenix/user/jzhang1/INTT-EMCAL/InttSeedingTrackDev/ParticleGen/output/",
0068     const string &output_filename = "results_test.root"
0069 )
0070 {
0071     std::cout << "CEMC_Clusters function is located at: " << (void*) CEMC_Clusters << std::endl;
0072
0073   const int skip = 0;
0074     const bool is_pythia = false;
0075   const int Nparticle = 1;
0076
0077   // const bool is_pythia = true;
0078   map<string, string> particle_map  = {
0079     {"MuonPlus", "mu+"},
0080     {"Muon", "mu-"},
0081     {"PionMinus", "pi-"},
0082     {"PionPlus", "pi+"},
0083     {"Electron", "e-"},
0084     {"Positron", "e+"}
0085     // {"Proton", "proton"},
0086     // {"Kaon", "kaon-"}
0087   };
0088
0089   if (particle_map.find(particle_species) == particle_map.end())
0090   {
0091     cout << "Invalid particle species: " << particle_species << endl;
0092     return;
0093   }
0094
0095   Fun4AllServer *se = Fun4AllServer::instance();
0096   //  se->Print("NODETREE"); // useless
0097   se->Verbosity(0);
0098
0099   //Opt to print all random seed used for debugging reproducibility. Comment out to reduce stdout prints.
0100   //PHRandomSeed::Verbosity(1);
0101
0102   // just if we set some flags somewhere in this macro
0103   recoConsts *rc = recoConsts::instance();
0104   // By default every random number generator uses
0105   // PHRandomSeed() which reads /dev/urandom to get its seed
0106   // if the RANDOMSEED flag is set its value is taken as seed
0107   // You can either set this to a random value using PHRandomSeed()
0108   // which will make all seeds identical (not sure what the point of
0109   // this would be:
0110   //  rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
0111   // or set it to a fixed value so you can debug your code
0112   //  rc->set_IntFlag("RANDOMSEED", 12345);
0113
0114   // INPUTREADHITS::filename[0] = inputFile;
0115
0116   // Force to use Pythia8 or GUN
0117   Input::PYTHIA8 = is_pythia;
0118   Input::GUN = !( is_pythia );
0119   Input::SIMPLE = true;
0120   if( Input::PYTHIA8 == false )
0121     {
0122       Input::GUN_NUMBER = 1; // if you need 3 of them
0123       Input::GUN_VERBOSITY = 1;
0124     }
0125
0126   // Initialize the selected Input/Event generation
0127   InputInit();
0128
0129   ///////////////
0130   // Set generator specific options
0131   ////////////////
0132   // can only be set after InputInit() is called
0133
0134   // pythia8
0135   if (Input::PYTHIA8)
0136     {
0137       //! Nominal collision geometry is selected by Input::BEAM_CONFIGURATION
0138       Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia8);
0139     }
0140   else // for GUN
0141     {
0142       // particle gun
0143       // if you run more than one of these Input::GUN_NUMBER > 1
0144       // add the settings for other with [1], next with [2]...
0145       // INPUTGENERATOR::Gun[0]->AddParticle("pi-", 0, 1, 0.1 );
0146       // INPUTGENERATOR::Gun[0]->set_vtx(0, 0, 0);
0147
0148       INPUTGENERATOR::SimpleEventGenerator[0] -> add_particles(particle_map[particle_species.c_str()], Nparticle );
0149       // INPUTGENERATOR::SimpleEventGenerator[0] -> set_name("mu-");
0150       INPUTGENERATOR::SimpleEventGenerator[0] -> set_vtx(0, 0, 0);
0151       INPUTGENERATOR::SimpleEventGenerator[0] -> set_pt_range(particle_pT - energy_range_down, particle_pT + energy_range_up);
0152       // INPUTGENERATOR::SimpleEventGenerator[0] -> set_mom(120, 0, 0);
0153       INPUTGENERATOR::SimpleEventGenerator[0] -> set_eta_range(-1.1, 1.1); // +-1
0154       // INPUTGENERATOR::SimpleEventGenerator[0] -> set_eta_range(-0.01, 0.01);
0155       INPUTGENERATOR::SimpleEventGenerator[0] -> set_phi_range(-M_PI, M_PI);
0156     }
0157
0158
0159
0160   InputRegister();
0161
0162   // Flag Handler is always needed to read flags from input (if used)
0163   // and update our rc flags with them. At the end it saves all flags
0164   // again on the DST in the Flags node under the RUN node
0165   FlagHandler *flag = new FlagHandler();
0166   se->registerSubsystem(flag);
0167
0168   //======================
0169   // Write the DST
0170   //======================
0171
0172   //Enable::DSTOUT = true;
0173   Enable::DSTOUT_COMPRESS = false;
0174   // DstOut::OutputDir = outdir;
0175   // DstOut::OutputFile = outputFile;
0176
0177   //Option to convert DST to human command readable TTree for quick poke around the outputs
0178   //  Enable::DSTREADER = true;
0179
0180   // turn the display on (default off)
0181   //Enable::DISPLAY = true;
0182
0183   //======================
0184   // What to run
0185   //======================
0186   // QA, main switch
0187   Enable::QA                    = false;
0188
0189   Enable::PIPE                  = true;
0190   Enable::PIPE_ABSORBER = true;
0191
0192   // central tracking
0193   Enable::MVTX                  = true;
0194   Enable::MVTX_CELL         = Enable::MVTX && true;
0195   Enable::MVTX_CLUSTER  = Enable::MVTX_CELL && true;
0196   Enable::MVTX_QA               = Enable::MVTX_CLUSTER && Enable::QA && true;
0197
0198   Enable::INTT                  = true;
0199   //  Enable::INTT_ABSORBER         = true; // enables layerwise support structure readout
0200   //  Enable::INTT_SUPPORT          = true; // enable global support structure readout
0201   Enable::INTT_CELL             = Enable::INTT && true;
0202   Enable::INTT_CLUSTER      = Enable::INTT_CELL && true;
0203   Enable::INTT_QA                 = Enable::INTT_CLUSTER && Enable::QA && true;
0204
0205   Enable::TPC                     = true;
0206   Enable::TPC_ABSORBER  = true;
0207   Enable::TPC_CELL          = Enable::TPC && true;
0208   Enable::TPC_CLUSTER       = Enable::TPC_CELL &&true;
0209   Enable::TPC_QA                = Enable::TPC_CLUSTER && Enable::QA && false;
0210
0211   Enable::MICROMEGAS                = true;
0212   Enable::MICROMEGAS_CELL           = Enable::MICROMEGAS && true;
0213   Enable::MICROMEGAS_CLUSTER    = Enable::MICROMEGAS_CELL && true;
0214   Enable::MICROMEGAS_QA             = Enable::MICROMEGAS_CLUSTER && Enable::QA && false;
0215
0216   Enable::TRACKING_TRACK    =  true;
0217
0218   Enable::TRACKING_EVAL     = Enable::TRACKING_TRACK && false;
0219   Enable::TRACKING_QA           = Enable::TRACKING_TRACK && Enable::QA && true; // <-- this one
0220
0221
0222   //! forward flux return plug door. Out of acceptance and off by default.
0223   //Enable::PLUGDOOR                = true;
0224   Enable::PLUGDOOR_ABSORBER         = true;
0225
0226
0227   Enable::CEMC = true;
0228   Enable::CEMC_ABSORBER = true;
0229   Enable::CEMC_CELL = Enable::CEMC && true;
0230   Enable::CEMC_TOWER = Enable::CEMC_CELL && true;
0231   Enable::CEMC_CLUSTER = Enable::CEMC_TOWER && true;
0232   Enable::CEMC_EVAL = Enable::CEMC_CLUSTER && true;
0233   Enable::CEMC_QA = Enable::CEMC_CLUSTER && Enable::QA && false;
0234
0235   Enable::HCALIN = true;
0236   Enable::HCALIN_ABSORBER = true;
0237   Enable::HCALIN_CELL = Enable::HCALIN && true;
0238   Enable::HCALIN_TOWER = Enable::HCALIN_CELL && true;
0239   Enable::HCALIN_CLUSTER = Enable::HCALIN_TOWER && true;
0240   Enable::HCALIN_EVAL = Enable::HCALIN_CLUSTER && true;
0241   Enable::HCALIN_QA = Enable::HCALIN_CLUSTER && Enable::QA && false;
0242
0243
0244   Enable::MAGNET                = true;
0245   Enable::MAGNET_ABSORBER           = true;
0246
0247
0248   Enable::HCALOUT = true;
0249   Enable::HCALOUT_ABSORBER = true;
0250   Enable::HCALOUT_CELL = Enable::HCALOUT && true;
0251   Enable::HCALOUT_TOWER = Enable::HCALOUT_CELL && true;
0252   Enable::HCALOUT_CLUSTER = Enable::HCALOUT_TOWER && true;
0253   Enable::HCALOUT_EVAL = Enable::HCALOUT_CLUSTER && true;
0254   Enable::HCALOUT_QA = Enable::HCALOUT_CLUSTER && Enable::QA && false;
0255
0256
0257   // new settings using Enable namespace in GlobalVariables.C
0258   Enable::BLACKHOLE             = true;
0259   //Enable::BLACKHOLE_SAVEHITS          = false; // turn off saving of bh hits
0260   //Enable::BLACKHOLE_FORWARD_SAVEHITS      = false; // disable forward/backward hits
0261   //BlackHoleGeometry::visible          = true;
0262
0263   //////////////////////////////
0264   // conditions DB flags      //
0265   //////////////////////////////
0266   Enable::CDB                   = true;
0267   // global tag
0268   rc->set_StringFlag("CDB_GLOBALTAG",CDB::global_tag);
0269   // 64 bit timestamp
0270   rc->set_uint64Flag("TIMESTAMP",CDB::timestamp);
0271
0272   /////////////////
0273   // Magnet Settings
0274   /////////////////
0275   G4MAGNET::magfield_rescale = 1.;  // for zero field
0276
0277   // Initialize the selected subsystems
0278   G4Init();
0279
0280   ///////////////////////
0281   // GEANT4 Detector description
0282   ///////////////////////
0283   if (!Input::READHITS) G4Setup();
0284
0285   ///////////////////////
0286   // Detector Division //
0287   ///////////////////////
0288
0289   //  if ((Enable::MBD && Enable::MBDRECO) || Enable::MBDFAKE) Mbd_Reco();
0290   if (Enable::MVTX_CELL) Mvtx_Cells();
0291   if (Enable::INTT_CELL) Intt_Cells();
0292
0293   //////////////////////////////////
0294   // CEMC towering and clustering //
0295   //////////////////////////////////
0296   if (Enable::CEMC_CELL) CEMC_Cells();
0297   if (Enable::CEMC_TOWER) CEMC_Towers();
0298   if (Enable::CEMC_CLUSTER) CEMC_Clusters();
0299
0300   if (Enable::HCALIN_CELL) HCALInner_Cells();
0301   if (Enable::HCALOUT_CELL) HCALOuter_Cells();
0302   if (Enable::HCALIN_TOWER) HCALInner_Towers();
0303   if (Enable::HCALIN_CLUSTER) HCALInner_Clusters();
0304   if (Enable::HCALOUT_TOWER) HCALOuter_Towers();
0305   if (Enable::HCALOUT_CLUSTER) HCALOuter_Clusters();
0306
0307
0308   /////////////////
0309   // SVTX tracking
0310   ////////////////
0311   if(Enable::TRACKING_TRACK) TrackingInit();
0312   if (Enable::MVTX_CLUSTER) Mvtx_Clustering();
0313   if (Enable::INTT_CLUSTER) Intt_Clustering();
0314   if (Enable::TRACKING_TRACK) Tracking_Reco();
0315
0316   //////////////////////
0317   // Simulation evaluation
0318   //////////////////////
0319   // string outputroot = outputFile;
0320   // string remove_this = ".root";
0321   // size_t pos = outputroot.find(remove_this);
0322   // if (pos != string::npos) outputroot.erase(pos, remove_this.length());
0323
0324   // if (Enable::TRACKING_EVAL) Tracking_Eval(outputroot + "_g4svtx_eval.root");
0325
0326     // Load the modified geometry
0327     CaloGeomMappingv2 *cgm = new CaloGeomMappingv2();
0328     cgm->set_detector_name("CEMC");
0329     cgm->setTowerGeomNodeName("TOWERGEOM_CEMCv3");
0330     se->registerSubsystem(cgm);
0331
0332     TruthToSvtxTrack *t2t = new TruthToSvtxTrack();
0333     se->registerSubsystem(t2t);
0334
0335     auto truthProjection = new PHActsTrackProjection("TruthTrackProjection");
0336     bool doEMcalRadiusCorr = true; // Set to true to apply the radius correction
0337     float new_cemc_rad = 99;
0338     if (doEMcalRadiusCorr)
0339     {
0340         truthProjection->setLayerRadius(SvtxTrack::CEMC, new_cemc_rad);
0341     }
0342     se->registerSubsystem(truthProjection);
0343
0344     // // tutorial info get
0345     tutorial* analysis_module = new tutorial( "name", output_directory, output_filename);
0346     analysis_module->setTowerGeomNodeName("TOWERGEOM_CEMCv3");
0347     // string output_path = "tutorial_results_";
0348     // output_path += "pythia8_MC.root";
0349     // analysis_module->SetOutputPath(output_path);
0350     se->registerSubsystem(analysis_module);
0351
0352     // Flag Handler is always needed to read flags from input (if used)
0353     // and update our rc flags with them. At the end it saves all flags
0354     // again on the DST in the Flags node under the RUN node
0355     // FlagHandler *flag = new FlagHandler();
0356     // se->registerSubsystem(flag);
0357     se->skip(skip);
0358     se->run(nEvents);
0359     se->End();
0360
0361     std::cout << "All done!" << std::endl;
0362
0363     gSystem->Exit(0);
0364     return;
0365 }
0366
0367