GenFitTrackProp/macros/Fun4All_single_particle.C

0001 int Fun4All_single_particle (
0002         const int nEvents = 10,
0003         const char * inputFile = NULL,
0004         const char * outputFile = "SvtxClusters.root",
0005         const char * embed_input_file = NULL, //"Hijing_G4Hits.root",
0006         const int which_tracking = 15,
0007         const bool do_embedding = false,
0008         const int  particle_pid = 211
0009         )
0010 {
0011     //===============
0012     // Input options
0013     //===============
0014
0015     // Either:
0016     // read previously generated g4-hits files, in this case it opens a DST and skips
0017     // the simulations step completely. The G4Setup macro is only loaded to get information
0018     // about the number of layers used for the cell reco code
0019     //
0020     // In case reading production output, please double check your G4Setup_sPHENIX.C and G4_*.C consistent with those in the production macro folder
0021     // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
0022     const bool readhits = false;
0023     // Or:
0024     // read files in HepMC format (typically output from event generators like hijing or pythia)
0025     const bool readhepmc = false; // read HepMC files
0026     // Or:
0027     // Use pythia
0028     const bool runpythia8 = false;
0029     const bool runpythia6 = false;
0030     // else
0031     // Use particle generator (default simple generator)
0032     // or gun/ very simple generator
0033     const bool usegun = false;
0034     // And
0035     // Further choose to embed newly simulated events to a previous simulation. Not compatible with `readhits = true`
0036     // In case embedding into a production output, please double check your G4Setup_sPHENIX.C and G4_*.C consistent with those in the production macro folder
0037     // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
0038     //const bool do_embedding = true;
0039
0040     //======================
0041     // What to run
0042     //======================
0043
0044     bool do_bbc = true;
0045
0046     bool do_pipe = true;
0047
0048     bool do_svtx = true;
0049     bool do_svtx_cell = true;
0050     bool do_svtx_cluster = true;
0051     bool do_svtx_track = false;
0052     bool do_svtx_eval = false;
0053
0054 //  if(which_tracking == 14 || which_tracking == 15) {
0055 //      do_svtx_track = false;
0056 //      do_svtx_eval = false;
0057 //  }
0058
0059     bool do_preshower = false;
0060
0061     bool do_cemc = false;
0062     bool do_cemc_cell = do_cemc && true;
0063     bool do_cemc_twr = do_cemc_cell && true;
0064     bool do_cemc_cluster = do_cemc_twr && true;
0065     bool do_cemc_eval = do_cemc_cluster && true;
0066
0067     bool do_hcalin = false;
0068     bool do_hcalin_cell = do_hcalin && true;
0069     bool do_hcalin_twr = do_hcalin_cell && true;
0070     bool do_hcalin_cluster = do_hcalin_twr && true;
0071     bool do_hcalin_eval = do_hcalin_cluster && true;
0072
0073     bool do_magnet = false;
0074
0075     bool do_hcalout = false;
0076     bool do_hcalout_cell = do_hcalout && true;
0077     bool do_hcalout_twr = do_hcalout_cell && true;
0078     bool do_hcalout_cluster = do_hcalout_twr && true;
0079     bool do_hcalout_eval = do_hcalout_cluster && true;
0080
0081     bool do_global = false;
0082     bool do_global_fastsim = false;
0083
0084     bool do_jet_reco = false;
0085     bool do_jet_eval = false;
0086
0087     bool do_dst_compress = false;
0088
0089     //Option to convert DST to human command readable TTree for quick poke around the outputs
0090     bool do_DSTReader = false;
0091     //---------------
0092     // Load libraries
0093     //---------------
0094
0095     gSystem->Load("libfun4all.so");
0096     gSystem->Load("libg4detectors.so");
0097     gSystem->Load("libphhepmc.so");
0098     gSystem->Load("libg4testbench.so");
0099     gSystem->Load("libg4hough.so");
0100   gSystem->Load("libg4calo.so");
0101     gSystem->Load("libg4eval.so");
0102
0103     // establish the geometry and reconstruction setup
0104     gROOT->LoadMacro("G4Setup_sPHENIX.C");
0105     G4Init(do_svtx,do_preshower,do_cemc,do_hcalin,do_magnet,do_hcalout,do_pipe,which_tracking);
0106
0107     int absorberactive = 1; // set to 1 to make all absorbers active volumes
0108     //  const string magfield = "1.5"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
0109     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)
0110     const float magfield_rescale = 1.4/1.5; // scale the map to a 1.4 T field
0111
0112     //---------------
0113     // Fun4All server
0114     //---------------
0115
0116     Fun4AllServer *se = Fun4AllServer::instance();
0117     //se->Verbosity(100);
0118     // just if we set some flags somewhere in this macro
0119     recoConsts *rc = recoConsts::instance();
0120     // By default every random number generator uses
0121     // PHRandomSeed() which reads /dev/urandom to get its seed
0122     // if the RANDOMSEED flag is set its value is taken as seed
0123     // You ca neither set this to a random value using PHRandomSeed()
0124     // which will make all seeds identical (not sure what the point of
0125     // this would be:
0126     //  rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
0127     // or set it to a fixed value so you can debug your code
0128     //  rc->set_IntFlag("RANDOMSEED", 12345);
0129
0130     //-----------------
0131     // Event generation
0132     //-----------------
0133
0134     if (readhits)
0135     {
0136         // Get the hits from a file
0137         // The input manager is declared later
0138
0139         if (do_embedding)
0140         {
0141             cout <<"Do not support read hits and embed background at the same time."<<endl;
0142             exit(1);
0143         }
0144
0145     }
0146     else if (readhepmc)
0147     {
0148         // this module is needed to read the HepMC records into our G4 sims
0149         // but only if you read HepMC input files
0150         HepMCNodeReader *hr = new HepMCNodeReader();
0151         se->registerSubsystem(hr);
0152     }
0153     else if (runpythia8)
0154     {
0155         gSystem->Load("libPHPythia8.so");
0156
0157         PHPy8JetTrigger *theTrigger = new PHPy8JetTrigger();
0158         //theTrigger->Verbosity(10);
0159         theTrigger->SetEtaHighLow(-0.6, 0.6);
0160         theTrigger->SetJetR(.4);
0161         theTrigger->SetMinJetPt(20);
0162
0163         PHPythia8* pythia8 = new PHPythia8();
0164         // see coresoftware/generators/PHPythia8 for example config
0165         pythia8->set_config_file("phpythia8.cfg");
0166         pythia8->register_trigger(theTrigger);
0167         se->registerSubsystem(pythia8);
0168
0169         HepMCNodeReader *hr = new HepMCNodeReader();
0170         hr->Embed(10);
0171         se->registerSubsystem(hr);
0172     }
0173     else if (runpythia6)
0174     {
0175         gSystem->Load("libPHPythia6.so");
0176
0177         PHPythia6 *pythia6 = new PHPythia6();
0178         pythia6->set_config_file("phpythia6.cfg");
0179         se->registerSubsystem(pythia6);
0180
0181         HepMCNodeReader *hr = new HepMCNodeReader();
0182         se->registerSubsystem(hr);
0183     }
0184     else
0185     {
0186         // toss low multiplicity dummy events
0187         PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
0188         //gen->set_seed(1710374143);
0189         // mu+,e+,proton,pi+,kaon+,Upsilon
0190         gen->add_particles(particle_pid,1);
0191 //      gen->add_particles("pi+",5);
0192 //      gen->add_particles("pi-",5);
0193         if (readhepmc || do_embedding)
0194         {
0195             gen->set_reuse_existing_vertex(true);
0196             gen->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
0197         }
0198         else
0199         {
0200             gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
0201                     PHG4SimpleEventGenerator::Uniform,
0202                     PHG4SimpleEventGenerator::Uniform);
0203             gen->set_vertex_distribution_mean(0.0, 0.0, 0.0);
0204             gen->set_vertex_distribution_width(0.0, 0.0, 0.0);
0205         }
0206         gen->set_vertex_size_function(PHG4SimpleEventGenerator::Uniform);
0207         gen->set_vertex_size_parameters(0.0, 0.0);
0208         gen->set_eta_range(0,0);
0209         gen->set_phi_range(-1 * TMath::Pi(), 1 * TMath::Pi());
0210 //      gen->set_eta_range(0, 0);
0211 //      gen->set_phi_range(0, 0);
0212         gen->set_pt_range(0.1, 5);
0213         gen->Embed(10);
0214         gen->Verbosity(0);
0215         if (! usegun)
0216         {
0217             se->registerSubsystem(gen);
0218         }
0219         else
0220         {
0221             PHG4ParticleGun *gun = new PHG4ParticleGun();
0222             //  gun->set_name("anti_proton");
0223             gun->set_name("geantino");
0224             gun->set_vtx(0, 0, 0);
0225             gun->set_mom(10, 0, 0.01);
0226             // gun->AddParticle("geantino",1.7776,-0.4335,0.);
0227             // gun->AddParticle("geantino",1.7709,-0.4598,0.);
0228             // gun->AddParticle("geantino",2.5621,0.60964,0.);
0229             // gun->AddParticle("geantino",1.8121,0.253,0.);
0230             //    se->registerSubsystem(gun);
0231             PHG4ParticleGenerator *pgen = new PHG4ParticleGenerator();
0232             pgen->set_name("geantino");
0233             pgen->set_z_range(0,0);
0234             pgen->set_eta_range(0.01,0.01);
0235             pgen->set_mom_range(10,10);
0236             pgen->set_phi_range(5.3/180.*TMath::Pi(),5.3/180.*TMath::Pi());
0237             se->registerSubsystem(pgen);
0238             pgen = new PHG4ParticleGenerator();
0239             pgen->set_name("geantino");
0240             pgen->set_z_range(0,0);
0241             pgen->set_eta_range(0.01,0.01);
0242             pgen->set_mom_range(10,10);
0243             pgen->set_phi_range(-0.2/180.*TMath::Pi(),0.2/180.*TMath::Pi());
0244             se->registerSubsystem(pgen);
0245         }
0246     }
0247
0248     if (!readhits)
0249     {
0250         //---------------------
0251         // Detector description
0252         //---------------------
0253
0254         G4Setup(absorberactive, magfield, TPythia6Decayer::kAll,
0255                 do_svtx, do_preshower, do_cemc, do_hcalin, do_magnet, do_hcalout, do_pipe, magfield_rescale);
0256     }
0257
0258     //---------
0259     // BBC Reco
0260     //---------
0261
0262     if (do_bbc)
0263     {
0264         gROOT->LoadMacro("G4_Bbc.C");
0265         BbcInit();
0266         Bbc_Reco();
0267     }
0268     //------------------
0269     // Detector Division
0270     //------------------
0271
0272     if (do_svtx_cell) Svtx_Cells();
0273
0274     if (do_cemc_cell) CEMC_Cells();
0275
0276     if (do_hcalin_cell) HCALInner_Cells();
0277
0278     if (do_hcalout_cell) HCALOuter_Cells();
0279
0280     //-----------------------------
0281     // CEMC towering and clustering
0282     //-----------------------------
0283
0284     if (do_cemc_twr) CEMC_Towers();
0285     if (do_cemc_cluster) CEMC_Clusters();
0286
0287     //-----------------------------
0288     // HCAL towering and clustering
0289     //-----------------------------
0290
0291     if (do_hcalin_twr) HCALInner_Towers();
0292     if (do_hcalin_cluster) HCALInner_Clusters();
0293
0294     if (do_hcalout_twr) HCALOuter_Towers();
0295     if (do_hcalout_cluster) HCALOuter_Clusters();
0296
0297     if (do_dst_compress) ShowerCompress();
0298
0299     //--------------
0300     // SVTX tracking
0301     //--------------
0302
0303     if(which_tracking == 14 || which_tracking == 15) {
0304         if (do_svtx_cluster) Svtx_Clustering();
0305         if (do_svtx_track) Svtx_Tracking();
0306     } else {
0307         if (do_svtx_cluster || do_svtx_track) Svtx_Reco();
0308     }
0309
0310     //-----------------
0311     // Global Vertexing
0312     //-----------------
0313
0314     if (do_global)
0315     {
0316         gROOT->LoadMacro("G4_Global.C");
0317         Global_Reco();
0318     }
0319
0320     else if (do_global_fastsim)
0321     {
0322         gROOT->LoadMacro("G4_Global.C");
0323         Global_FastSim();
0324     }
0325
0326     //---------
0327     // Jet reco
0328     //---------
0329
0330     if (do_jet_reco)
0331     {
0332         gROOT->LoadMacro("G4_Jets.C");
0333         Jet_Reco();
0334     }
0335
0336     // HF jet trigger moudle
0337     //assert (gSystem->Load("libHFJetTruthGeneration") == 0);
0338     //{
0339     //  if (do_jet_reco)
0340     //  {
0341     //      HFJetTruthTrigger * jt = new HFJetTruthTrigger(
0342     //              "HFJetTruthTrigger.root", 5 , "AntiKt_Truth_r04");
0343     //      //jt->Verbosity(HFJetTruthTrigger::VERBOSITY_MORE);
0344     //      jt->set_pt_min(20);
0345     //      se->registerSubsystem(jt);
0346     //  }
0347     //}
0348
0349     //----------------------
0350     // Simulation evaluation
0351     //----------------------
0352
0353     if (do_svtx_eval) Svtx_Eval("g4svtx_eval.root");
0354
0355     if (do_cemc_eval) CEMC_Eval("g4cemc_eval.root");
0356
0357     if (do_hcalin_eval) HCALInner_Eval("g4hcalin_eval.root");
0358
0359     if (do_hcalout_eval) HCALOuter_Eval("g4hcalout_eval.root");
0360
0361     if (do_jet_eval) Jet_Eval("g4jet_eval.root");
0362
0363     //--------------
0364     // IO management
0365     //--------------
0366
0367     if (readhits)
0368     {
0369         // Hits file
0370         Fun4AllInputManager *hitsin = new Fun4AllDstInputManager("DSTin");
0371         hitsin->fileopen(inputFile);
0372         se->registerInputManager(hitsin);
0373     }
0374     if (do_embedding)
0375     {
0376         if (embed_input_file == NULL)
0377         {
0378             cout << "Missing embed_input_file! Exit";
0379             exit(3);
0380         }
0381
0382         Fun4AllDstInputManager *in1 = new Fun4AllNoSyncDstInputManager("DSTinEmbed");
0383         in1->AddFile(embed_input_file); // if one use a single input file
0384         //in1->AddListFile(embed_input_file); // RecommendedL: if one use a text list of many input files
0385         se->registerInputManager(in1);
0386     }
0387     if (readhepmc)
0388     {
0389         Fun4AllInputManager *in = new Fun4AllHepMCInputManager( "DSTIN");
0390         se->registerInputManager( in );
0391         se->fileopen( in->Name().c_str(), inputFile );
0392     }
0393     else
0394     {
0395         // for single particle generators we just need something which drives
0396         // the event loop, the Dummy Input Mgr does just that
0397         Fun4AllInputManager *in = new Fun4AllDummyInputManager( "JADE");
0398         se->registerInputManager( in );
0399     }
0400
0401     if (do_DSTReader)
0402     {
0403         //Convert DST to human command readable TTree for quick poke around the outputs
0404         gROOT->LoadMacro("G4_DSTReader.C");
0405
0406         G4DSTreader( outputFile, //
0407                 /*int*/ absorberactive ,
0408                 /*bool*/ do_svtx ,
0409                 /*bool*/ do_preshower ,
0410                 /*bool*/ do_cemc ,
0411                 /*bool*/ do_hcalin ,
0412                 /*bool*/ do_magnet ,
0413                 /*bool*/ do_hcalout ,
0414                 /*bool*/ do_cemc_twr ,
0415                 /*bool*/ do_hcalin_twr ,
0416                 /*bool*/ do_magnet  ,
0417                 /*bool*/ do_hcalout_twr
0418                 );
0419     }
0420
0421     Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputFile);
0422     if (do_dst_compress) DstCompress(out);
0423     se->registerOutputManager(out);
0424
0425     //-----------------
0426     // Event processing
0427     //-----------------
0428     if (nEvents < 0)
0429     {
0430         return;
0431     }
0432     // if we run the particle generator and use 0 it'll run forever
0433     if (nEvents == 0 && !readhits && !readhepmc)
0434     {
0435         cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
0436         cout << "it will run forever, so I just return without running anything" << endl;
0437         return;
0438     }
0439
0440     se->run(nEvents);
0441
0442     //-----
0443     // Exit
0444     //-----
0445
0446     se->End();
0447
0448     //getchar();
0449
0450     std::cout << "All done" << std::endl;
0451     //delete se;
0452     //gSystem->Exit(0);
0453 }