HighDCATrackCounting/macros/Jin_BJet.C

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