HCalJetPhiShift/macros/Fun4All_HCalJetPhiShift.C

0001 #ifndef MACRO_FUN4ALLG4SPHENIX_C
0002 #define MACRO_FUN4ALLG4SPHENIX_C
0003
0004 #include <GlobalVariables.C>
0005
0006 #include <DisplayOn.C>
0007 #include <G4Setup_sPHENIX.C>
0008 #include <G4_Bbc.C>
0009 #include <G4_CaloTrigger.C>
0010 #include <G4_Centrality.C>
0011 #include <G4_DSTReader.C>
0012 #include <G4_Global.C>
0013 #include <G4_HIJetReco.C>
0014 #include <G4_Input.C>
0015 #include <G4_Jets.C>
0016 #include <G4_KFParticle.C>
0017 #include <G4_ParticleFlow.C>
0018 #include <G4_Production.C>
0019 #include <G4_TopoClusterReco.C>
0020
0021 #include <G4_User.C>
0022 #include <QA.C>
0023
0024 #include <ffamodules/FlagHandler.h>
0025 #include <ffamodules/HeadReco.h>
0026 #include <ffamodules/SyncReco.h>
0027 #include <ffamodules/CDBInterface.h>
0028
0029 #include <fun4all/Fun4AllDstOutputManager.h>
0030 #include <fun4all/Fun4AllOutputManager.h>
0031 #include <fun4all/Fun4AllServer.h>
0032
0033 #include <phool/PHRandomSeed.h>
0034 #include <phool/recoConsts.h>
0035
0036 #include <hcaljetphishift/HCalJetPhiShift.h>
0037
0038 R__LOAD_LIBRARY(libfun4all.so)
0039 R__LOAD_LIBRARY(libffamodules.so)
0040 R__LOAD_LIBRARY(libHCalJetPhiShift.so)
0041
0042 // For HepMC Hijing
0043 // try inputFile = /sphenix/sim/sim01/sphnxpro/sHijing_HepMC/sHijing_0-12fm.dat
0044
0045 int Fun4All_HCalJetPhiShift(
0046     const int nEvents = 1,
0047     const int index = 0,
0048     const string &outdir = ".")
0049 {
0050
0051   int event_number = (int)index*nEvents;
0052
0053   //  string outputFile = "HCalJetPhiShift"
0054   string outputroot = outdir + "/HCalJetPhiShift";
0055   if (index>-1) {
0056     outputroot += "_" + to_string(index);
0057   }
0058
0059   string outputFile = outputroot + ".root";
0060
0061   Fun4AllServer *se = Fun4AllServer::instance();
0062   se->Verbosity(0);
0063
0064   //Opt to print all random seed used for debugging reproducibility. Comment out to reduce stdout prints.
0065 //  PHRandomSeed::Verbosity(1);
0066
0067   // just if we set some flags somewhere in this macro
0068   recoConsts *rc = recoConsts::instance();
0069   // By default every random number generator uses
0070   // PHRandomSeed() which reads /dev/urandom to get its seed
0071   // if the RANDOMSEED flag is set its value is taken as seed
0072   // You can either set this to a random value using PHRandomSeed()
0073   // which will make all seeds identical (not sure what the point of
0074   // this would be:
0075   //  rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
0076   // or set it to a fixed value so you can debug your code
0077   //  rc->set_IntFlag("RANDOMSEED", 12345);
0078
0079
0080   //===============
0081   // Input options
0082   //===============
0083   // verbosity setting (applies to all input managers)
0084   Input::VERBOSITY = 0;
0085
0086   Input::SIMPLE = true;
0087   // Input::SIMPLE_NUMBER = 2; // if you need 2 of them
0088   Input::SIMPLE_VERBOSITY = false;
0089
0090
0091   //-----------------
0092   // Initialize the selected Input/Event generation
0093   //-----------------
0094   // This creates the input generator(s)
0095   InputInit();
0096
0097   //--------------
0098   // Set generator specific options
0099   //--------------
0100   // can only be set after InputInit() is called
0101
0102   // Simple Input generator:
0103   // if you run more than one of these Input::SIMPLE_NUMBER > 1
0104   // add the settings for other with [1], next with [2]...
0105   if (Input::SIMPLE)
0106   {
0107     INPUTGENERATOR::SimpleEventGenerator[0]->add_particles("pi-", 1);
0108     INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_function(PHG4SimpleEventGenerator::Gaus,
0109                                                                             PHG4SimpleEventGenerator::Gaus,
0110                                                                             PHG4SimpleEventGenerator::Gaus);
0111     INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_mean(0., 0., 0.);
0112     INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_width(0.0, 0.0, 0.0);
0113     INPUTGENERATOR::SimpleEventGenerator[0]->set_eta_range(-1, 1);
0114     INPUTGENERATOR::SimpleEventGenerator[0]->set_phi_range(-M_PI, M_PI);
0115     INPUTGENERATOR::SimpleEventGenerator[0]->set_pt_range(1., 30.);
0116   }
0117
0118   //--------------
0119   // Set Input Manager specific options
0120   //--------------
0121
0122   // register all input generators with Fun4All
0123   InputRegister();
0124
0125   rc->set_IntFlag("RUNNUMBER",1);
0126
0127   SyncReco *sync = new SyncReco();
0128   se->registerSubsystem(sync);
0129
0130   HeadReco *head = new HeadReco();
0131   se->registerSubsystem(head);
0132
0133 // Flag Handler is always needed to read flags from input (if used)
0134 // and update our rc flags with them. At the end it saves all flags
0135 // again on the DST in the Flags node under the RUN node
0136   FlagHandler *flag = new FlagHandler();
0137   se->registerSubsystem(flag);
0138
0139   //======================
0140   // Write the DST
0141   //======================
0142
0143   //Enable::DSTOUT = true;
0144   Enable::DSTOUT_COMPRESS = false;
0145   DstOut::OutputDir = outdir;
0146   DstOut::OutputFile = outputFile;
0147
0148   //Option to convert DST to human command readable TTree for quick poke around the outputs
0149   //  Enable::DSTREADER = true;
0150
0151   // turn the display on (default off)
0152    //Enable::DISPLAY = true;
0153
0154   //======================
0155   // What to run
0156   //======================
0157
0158   // QA, main switch
0159   Enable::QA = false;
0160
0161   // Global options (enabled for all enables subsystems - if implemented)
0162   //  Enable::ABSORBER = true;
0163   //  Enable::OVERLAPCHECK = true;
0164   //  Enable::VERBOSITY = 1;
0165
0166   // Enable::BBC = true;
0167   // Enable::BBC_SUPPORT = true; // save hist in bbc support structure
0168   // Enable::BBCRECO = Enable::BBC && true
0169   Enable::BBCFAKE = false;  // Smeared vtx and t0, use if you don't want real BBC in simulation
0170
0171   Enable::PIPE = false;
0172   Enable::PIPE_ABSORBER = false;
0173
0174   // central tracking
0175   Enable::MVTX = false;
0176   Enable::MVTX_CELL = Enable::MVTX && true;
0177   Enable::MVTX_CLUSTER = Enable::MVTX_CELL && true;
0178   Enable::MVTX_QA = Enable::MVTX_CLUSTER && Enable::QA && true;
0179
0180   Enable::INTT = false;
0181 //  Enable::INTT_ABSORBER = true; // enables layerwise support structure readout
0182 //  Enable::INTT_SUPPORT = true; // enable global support structure readout
0183   Enable::INTT_CELL = Enable::INTT && true;
0184   Enable::INTT_CLUSTER = Enable::INTT_CELL && true;
0185   Enable::INTT_QA = Enable::INTT_CLUSTER && Enable::QA && true;
0186
0187   Enable::TPC = false;
0188   Enable::TPC_ABSORBER = false;
0189   Enable::TPC_CELL = Enable::TPC && true;
0190   Enable::TPC_CLUSTER = Enable::TPC_CELL && true;
0191   Enable::TPC_QA = Enable::TPC_CLUSTER && Enable::QA && true;
0192
0193   Enable::MICROMEGAS = false;
0194   Enable::MICROMEGAS_CELL = Enable::MICROMEGAS && true;
0195   Enable::MICROMEGAS_CLUSTER = Enable::MICROMEGAS_CELL && true;
0196   Enable::MICROMEGAS_QA = Enable::MICROMEGAS_CLUSTER && Enable::QA && true;
0197
0198   Enable::TRACKING_TRACK = (Enable::MICROMEGAS_CLUSTER && Enable::TPC_CLUSTER && Enable::INTT_CLUSTER && Enable::MVTX_CLUSTER) && true;
0199   Enable::TRACKING_EVAL = Enable::TRACKING_TRACK && true;
0200   Enable::TRACKING_QA = Enable::TRACKING_TRACK && Enable::QA && true;
0201
0202   //G4CEMC::TowerDigi = RawTowerDigitizer::kNo_digitization;
0203
0204   Enable::CEMC = true;
0205   Enable::CEMC_ABSORBER = true;
0206
0207   Enable::CEMC_CELL = Enable::CEMC && true;
0208   Enable::CEMC_TOWER = Enable::CEMC_CELL && true;
0209   Enable::CEMC_CLUSTER = Enable::CEMC_TOWER && true;
0210   Enable::CEMC_EVAL = Enable::CEMC_CLUSTER && true;
0211   Enable::CEMC_QA = Enable::CEMC_CLUSTER && Enable::QA && true;
0212
0213   Enable::HCALIN = true;
0214   Enable::HCALIN_ABSORBER = true;
0215   Enable::HCALIN_CELL = Enable::HCALIN && true;
0216   Enable::HCALIN_TOWER = Enable::HCALIN_CELL && true;
0217   Enable::HCALIN_CLUSTER = Enable::HCALIN_TOWER && true;
0218   Enable::HCALIN_EVAL = Enable::HCALIN_CLUSTER && true;
0219   Enable::HCALIN_QA = Enable::HCALIN_CLUSTER && Enable::QA && true;
0220
0221   Enable::MAGNET = false;
0222   Enable::MAGNET_ABSORBER = false;
0223
0224   Enable::HCALOUT = true;
0225   Enable::HCALOUT_ABSORBER = true;
0226   Enable::HCALOUT_CELL = Enable::HCALOUT && true;
0227   Enable::HCALOUT_TOWER = Enable::HCALOUT_CELL && true;
0228   Enable::HCALOUT_CLUSTER = Enable::HCALOUT_TOWER && true;
0229   Enable::HCALOUT_EVAL = Enable::HCALOUT_CLUSTER && true;
0230   Enable::HCALOUT_QA = Enable::HCALOUT_CLUSTER && Enable::QA && true;
0231
0232   Enable::EPD = false;
0233   Enable::EPD_TILE = Enable::EPD && true;
0234
0235   Enable::BEAMLINE = false;
0236 //  Enable::BEAMLINE_ABSORBER = true;  // makes the beam line magnets sensitive volumes
0237 //  Enable::BEAMLINE_BLACKHOLE = true; // turns the beamline magnets into black holes
0238   Enable::ZDC = false;
0239 //  Enable::ZDC_ABSORBER = true;
0240 //  Enable::ZDC_SUPPORT = true;
0241 //  Enable::ZDC_TOWER = Enable::ZDC && true;
0242 //  Enable::ZDC_EVAL = Enable::ZDC_TOWER && true;
0243
0244   //! forward flux return plug door. Out of acceptance and off by default.
0245   //Enable::PLUGDOOR = true;
0246   Enable::PLUGDOOR_ABSORBER = false;
0247
0248 //  Enable::GLOBAL_RECO = (Enable::BBCFAKE || Enable::TRACKING_TRACK) && true;
0249   //Enable::GLOBAL_FASTSIM = true;
0250
0251   //Enable::KFPARTICLE = true;
0252   //Enable::KFPARTICLE_VERBOSITY = 1;
0253   //Enable::KFPARTICLE_TRUTH_MATCH = true;
0254   //Enable::KFPARTICLE_SAVE_NTUPLE = true;
0255
0256   Enable::CALOTRIGGER = Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
0257
0258   Enable::JETS = (Enable::GLOBAL_RECO || Enable::GLOBAL_FASTSIM) && false;
0259   Enable::JETS_EVAL = Enable::JETS && true;
0260   Enable::JETS_QA = Enable::JETS && Enable::QA && true;
0261
0262   // HI Jet Reco for p+Au / Au+Au collisions (default is false for
0263   // single particle / p+p-only simulations, or for p+Au / Au+Au
0264   // simulations which don't particularly care about jets)
0265   Enable::HIJETS = Enable::JETS && Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
0266
0267   // 3-D topoCluster reconstruction, potentially in all calorimeter layers
0268   Enable::TOPOCLUSTER = Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
0269   // particle flow jet reconstruction - needs topoClusters!
0270   Enable::PARTICLEFLOW = Enable::TOPOCLUSTER && true;
0271   // centrality reconstruction
0272   Enable::CENTRALITY = false;
0273
0274   // new settings using Enable namespace in GlobalVariables.C
0275   Enable::BLACKHOLE = false;
0276   //Enable::BLACKHOLE_SAVEHITS = false; // turn off saving of bh hits
0277   //Enable::BLACKHOLE_FORWARD_SAVEHITS = false; // disable forward/backward hits
0278   //BlackHoleGeometry::visible = true;
0279
0280   // run user provided code (from local G4_User.C)
0281   //Enable::USER = true;
0282
0283   //===============
0284   // conditions DB flags
0285   //===============
0286   Enable::CDB = false;
0287   // global tag
0288   rc->set_StringFlag("CDB_GLOBALTAG",CDB::global_tag);
0289   // 64 bit timestamp
0290   rc->set_uint64Flag("TIMESTAMP",CDB::timestamp);
0291   //---------------
0292   // World Settings
0293   //---------------
0294   G4WORLD::PhysicsList = "FTFP_BERT"; //FTFP_BERT_HP best for calo
0295   //  G4WORLD::WorldMaterial = "G4_AIR"; // set to G4_GALACTIC for material scans
0296
0297   //---------------
0298   // Magnet Settings
0299   //---------------
0300
0301   //G4MAGNET::magfield =  string(getenv("CALIBRATIONROOT"))+ string("/Field/Map/sphenix3dbigmapxyz.root");  // default map from the calibration database
0302   G4MAGNET::magfield = "0."; // alternatively to specify a constant magnetic field, give a float number, which will be translated to solenoidal field in T, if string use as fieldmap name (including path)
0303 //  G4MAGNET::magfield_rescale = 1.;  // make consistent with expected Babar field strength of 1.4T
0304
0305   //---------------
0306   // Pythia Decayer
0307   //---------------
0308   // list of decay types in
0309   // $OFFLINE_MAIN/include/g4decayer/EDecayType.hh
0310   // default is All:
0311   // G4P6DECAYER::decayType = EDecayType::kAll;
0312
0313   // Initialize the selected subsystems
0314   G4Init();
0315
0316   //---------------------
0317   // GEANT4 Detector description
0318   //---------------------
0319   if (!Input::READHITS)
0320   {
0321     G4Setup();
0322   }
0323
0324   //------------------
0325   // Detector Division
0326   //------------------
0327
0328 //  if ((Enable::BBC && Enable::BBCRECO) || Enable::BBCFAKE) Bbc_Reco();
0329 //
0330 //  if (Enable::MVTX_CELL) Mvtx_Cells();
0331 //  if (Enable::INTT_CELL) Intt_Cells();
0332 //  if (Enable::TPC_CELL) TPC_Cells();
0333 //  if (Enable::MICROMEGAS_CELL) Micromegas_Cells();
0334
0335   if (Enable::CEMC_CELL) CEMC_Cells();
0336
0337   if (Enable::HCALIN_CELL) HCALInner_Cells();
0338
0339   if (Enable::HCALOUT_CELL) HCALOuter_Cells();
0340
0341   //-----------------------------
0342   // CEMC towering and clustering
0343   //-----------------------------
0344
0345   if (Enable::CEMC_TOWER) CEMC_Towers();
0346   if (Enable::CEMC_CLUSTER) CEMC_Clusters();
0347
0348   //--------------
0349   // EPD tile reconstruction
0350   //--------------
0351
0352 //  if (Enable::EPD_TILE) EPD_Tiles();
0353
0354   //-----------------------------
0355   // HCAL towering and clustering
0356   //-----------------------------
0357
0358   if (Enable::HCALIN_TOWER) HCALInner_Towers();
0359   if (Enable::HCALIN_CLUSTER) HCALInner_Clusters();
0360
0361   if (Enable::HCALOUT_TOWER) HCALOuter_Towers();
0362   if (Enable::HCALOUT_CLUSTER) HCALOuter_Clusters();
0363
0364   // if enabled, do topoClustering early, upstream of any possible jet reconstruction
0365 //  if (Enable::TOPOCLUSTER) TopoClusterReco();
0366
0367   //--------------
0368   // SVTX tracking
0369   //--------------
0370 //  if(Enable::TRACKING_TRACK)
0371 //    {
0372 //      TrackingInit();
0373 //    }
0374 //  if (Enable::MVTX_CLUSTER) Mvtx_Clustering();
0375 //  if (Enable::INTT_CLUSTER) Intt_Clustering();
0376 //  if (Enable::TPC_CLUSTER)
0377 //    {
0378 //      if(G4TPC::ENABLE_DIRECT_LASER_HITS || G4TPC::ENABLE_CENTRAL_MEMBRANE_HITS)
0379 //  {
0380 //    TPC_LaserClustering();
0381 //  }
0382 //      else
0383 //  {
0384 //    TPC_Clustering();
0385 //  }
0386 //    }
0387 //  if (Enable::MICROMEGAS_CLUSTER) Micromegas_Clustering();
0388 //
0389 //  if (Enable::TRACKING_TRACK)
0390 //  {
0391 //    Tracking_Reco();
0392 //  }
0393 //
0394 //  if(Enable::TRACKING_DIAGNOSTICS)
0395 //    {
0396 //      const std::string kshortFile = "./kshort_" + outputFile;
0397 //      const std::string residualsFile = "./residuals_" + outputFile;
0398 //
0399 //      G4KshortReconstruction(kshortFile);
0400 //      seedResiduals(residualsFile);
0401 //    }
0402
0403   //-----------------
0404   // Global Vertexing
0405   //-----------------
0406
0407 //  if (Enable::GLOBAL_RECO && Enable::GLOBAL_FASTSIM)
0408 //  {
0409 //    cout << "You can only enable Enable::GLOBAL_RECO or Enable::GLOBAL_FASTSIM, not both" << endl;
0410 //    gSystem->Exit(1);
0411 //  }
0412 //  if (Enable::GLOBAL_RECO)
0413 //  {
0414 //    Global_Reco();
0415 //  }
0416 //  else if (Enable::GLOBAL_FASTSIM)
0417 //  {
0418 //    Global_FastSim();
0419 //  }
0420
0421   //-----------------
0422   // Centrality Determination
0423   //-----------------
0424
0425   if (Enable::CENTRALITY)
0426   {
0427       Centrality();
0428   }
0429
0430   //-----------------
0431   // Calo Trigger Simulation
0432   //-----------------
0433
0434   if (Enable::CALOTRIGGER)
0435   {
0436     CaloTrigger_Sim();
0437   }
0438
0439   //---------
0440   // Jet reco
0441   //---------
0442
0443 //  if (Enable::JETS) Jet_Reco();
0444 //  if (Enable::HIJETS) HIJetReco();
0445 //
0446 //  if (Enable::PARTICLEFLOW) ParticleFlow();
0447
0448   //----------------------
0449   // Simulation evaluation
0450   //----------------------
0451 //  string outputroot = outputFile;
0452 //  string remove_this = ".root";
0453 //  size_t pos = outputroot.find(remove_this);
0454 //  if (pos != string::npos)
0455 //  {
0456 //    outputroot.erase(pos, remove_this.length());
0457 //  }
0458 //  int start_event = (int) index*10;
0459 //  if (Enable::HCALOUT_EVAL) HCALOuter_Eval(outputroot + "_g4hcalout_eval.root", start_event);
0460 //  if (Enable::HCALIN_EVAL) HCALInner_Eval(outputroot + "_g4hcalin_eval.root", start_event);
0461
0462   //--------------
0463   // Set up Input Managers
0464   //--------------
0465
0466   InputManagers();
0467
0468   if (Enable::PRODUCTION)
0469   {
0470     Production_CreateOutputDir();
0471   }
0472
0473   if (Enable::DSTOUT)
0474   {
0475     string FullOutFile = DstOut::OutputDir + "/" + DstOut::OutputFile;
0476     Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
0477     if (Enable::DSTOUT_COMPRESS)
0478     {
0479       ShowerCompress();
0480       DstCompress(out);
0481     }
0482     se->registerOutputManager(out);
0483   }
0484   //-----------------
0485   // Event processing
0486   //-----------------
0487
0488   // if we use a negative number of events we go back to the command line here
0489   if (nEvents < 0)
0490   {
0491     return 0;
0492   }
0493   // if we run the particle generator and use 0 it'll run forever
0494   // for embedding it runs forever if the repeat flag is set
0495   if (nEvents == 0 && !Input::HEPMC && !Input::READHITS && INPUTEMBED::REPEAT)
0496   {
0497     cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
0498     cout << "it will run forever, so I just return without running anything" << endl;
0499     return 0;
0500   }
0501
0502   RawTowerCalibration *TowerCalibration = new RawTowerCalibration("EmcRawTowerCalibration");
0503   TowerCalibration->Detector("CEMC");
0504   TowerCalibration->set_calib_algorithm(RawTowerCalibration::kSimple_linear_calibration);
0505   TowerCalibration->set_calib_const_GeV_ADC(1. / 0.023);
0506   TowerCalibration->set_pedstal_ADC(0);
0507   se->registerSubsystem(TowerCalibration);
0508
0509   HCalJetPhiShift *caloPhiShift = new HCalJetPhiShift("caloPhiShift",outputFile);
0510   caloPhiShift->SetEventNumber(event_number);
0511   se->registerSubsystem(caloPhiShift);
0512   se->run(nEvents);
0513
0514   //-----
0515   // QA output
0516   //-----
0517
0518   if (Enable::QA) QA_Output(outputroot + "_qa.root");
0519
0520   //-----
0521   // Exit
0522   //-----
0523
0524 //  CDBInterface::instance()->Print(); // print used DB files
0525   se->End();
0526   std::cout << "All done" << std::endl;
0527   delete se;
0528
0529   gSystem->Exit(0);
0530   return 0;
0531 }
0532 #endif
0533
0534 // initial val setter in evaluator, set ievent --> call function to set ievent in macro
0535 // then modify g4 macro or put line in my macro