VandySkimmer/macros/Fun4All_VandySkimmerTruth.C

0001 #include <string>
0002 #include <sstream>
0003 #include <iostream>
0004 #include <format>
0005
0006 #include <fun4all/Fun4AllBase.h>
0007 #include <fun4all/Fun4AllUtils.h>
0008 #include <fun4all/Fun4AllServer.h>
0009 #include <fun4all/Fun4AllInputManager.h>
0010 #include <fun4all/Fun4AllDstInputManager.h>
0011 #include <fun4all/SubsysReco.h>
0012 #include <ffamodules/CDBInterface.h>
0013 #include "Calo_CalibLocal.C"
0014
0015 #include <caloreco/RawClusterBuilderTopo.h>
0016
0017 #include <jetbase/JetReco.h>
0018 #include <jetbase/TowerJetInput.h>
0019 #include <jetbase/JetCalib.h>
0020 #include <jetbackground/FastJetAlgoSub.h>
0021
0022 #include <g4jets/TruthJetInput.h>
0023 //#include <globalvertex/GlobalVertexReco.h>
0024 //#include <GlobalVertex.h>
0025 //#include <MbdDigitization.h>
0026 //#include <MbdReco.h>
0027
0028 #include <jetbackground/TimingCut.h>
0029
0030 #include <vandyskimmer/VandyJetDSTSkimmer.h>
0031 #include <jetbackground/RetowerCEMC.h>
0032
0033
0034 R__LOAD_LIBRARY(libfun4all.so)
0035 R__LOAD_LIBRARY(libfun4allraw.so)
0036 R__LOAD_LIBRARY(libcalo_io.so)
0037 R__LOAD_LIBRARY(libffamodules.so)
0038 R__LOAD_LIBRARY(libVandySkimmer.so)
0039 R__LOAD_LIBRARY(libjetbase.so)
0040 R__LOAD_LIBRARY(libjetbackground.so)
0041 R__LOAD_LIBRARY(libg4dst.so)
0042
0043 void Fun4All_VandySkimmerTruth(const std::string caloDSTlist, const std::string jetDSTlist, const std::string g4HitsDSTlist, const std::string globalDSTlist, const std::string outDir = "/sphenix/tg/tg01/jets/bkimelman/wEEC/", const std::string nfiles="25")
0044 {
0045
0046   bool doSim = true;
0047
0048   Fun4AllServer* se=Fun4AllServer::instance();
0049   se->Verbosity(0);
0050
0051
0052   int runnumber = 0;
0053   int seg = 0;
0054   int n=0;
0055   n = std::stoi(nfiles);
0056
0057   std::ifstream ifs(caloDSTlist);
0058   std::string filepath;
0059   std::getline(ifs,filepath);
0060   std::pair<int,int> runseg = Fun4AllUtils::GetRunSegment(filepath);
0061   runnumber = runseg.first;
0062   seg = runseg.second;
0063
0064   std::string sample_name {"Jet20"};
0065   std::stringstream fullfilename (caloDSTlist);
0066   std::string temp1, temp2;
0067   while(std::getline(fullfilename, temp1, '/'))
0068   {
0069       if(temp1.find("Jet") == std::string::npos) continue;
0070       std::stringstream filetag (temp1);
0071       while(std::getline(filetag, temp2, '_'))
0072       {
0073           if(temp2.find("data") == std::string::npos){
0074               sample_name = temp2;
0075               break;
0076           }
0077       }
0078       break;
0079   }
0080   Fun4AllInputManager *inCalo = new Fun4AllDstInputManager("InputManagerCalo");
0081   inCalo->AddListFile(caloDSTlist);
0082   se->registerInputManager(inCalo);
0083
0084   Fun4AllInputManager *inJet = new Fun4AllDstInputManager("InputManagerJet");
0085   inJet->AddListFile(jetDSTlist);
0086   se->registerInputManager(inJet);
0087
0088   Fun4AllInputManager *inTruth = new Fun4AllDstInputManager("InputManagerG4Hits");
0089   inTruth->AddListFile(g4HitsDSTlist);
0090   se->registerInputManager(inTruth);
0091
0092   Fun4AllInputManager *inGlobal = new Fun4AllDstInputManager("InputManagerGlobal");
0093   inGlobal->AddListFile(globalDSTlist);
0094   se->registerInputManager(inGlobal);
0095
0096   auto rc = recoConsts::instance();
0097   rc->set_StringFlag("CDB_GLOBALTAG", "MDC2");
0098   rc->set_uint64Flag("TIMESTAMP", 28);
0099   CDBInterface::instance()->Verbosity(0);
0100   Process_Calo_Calib();
0101
0102   RawClusterBuilderTopo* ClusterBuilder = new RawClusterBuilderTopo("HcalRawClusterBuilderTopo");
0103   ClusterBuilder->Verbosity(0);
0104   ClusterBuilder->set_nodename("TOPOCLUSTER_ALLCALO");
0105   ClusterBuilder->set_enable_HCal(true);
0106   ClusterBuilder->set_enable_EMCal(true);
0107   ClusterBuilder->set_noise(0.0053, 0.0351, 0.0684); // 3sigma of pedestal noise
0108   ClusterBuilder->set_significance(4.0, 2.0, 1.0);
0109   ClusterBuilder->allow_corner_neighbor(true);
0110   ClusterBuilder->set_do_split(true);
0111   ClusterBuilder->set_minE_local_max(1.0, 2.0, 0.5);
0112   ClusterBuilder->set_R_shower(0.025);
0113   ClusterBuilder->set_use_only_good_towers(true);
0114   ClusterBuilder->set_absE(true);
0115   se->registerSubsystem(ClusterBuilder);
0116
0117   RetowerCEMC* rtcemc = new RetowerCEMC("RetowerCEMV");
0118   rtcemc->Verbosity(0);
0119   rtcemc->set_towerinfo(true);
0120   rtcemc->set_frac_cut(0.5);
0121   se->registerSubsystem(rtcemc);
0122
0123   JetReco *towerjetreco = new JetReco("towerJetReco");
0124   TowerJetInput* emtji = new TowerJetInput(Jet::CEMC_TOWERINFO_RETOWER,"TOWERINFO_CALIB_RETOWER");
0125   TowerJetInput* ohtji = new TowerJetInput(Jet::HCALIN_TOWERINFO,"TOWERINFO_CALIB");
0126   TowerJetInput* ihtji = new TowerJetInput(Jet::HCALOUT_TOWERINFO,"TOWERINFO_CALIB");
0127   emtji->set_GlobalVertexType(GlobalVertex::VTXTYPE::MBD);
0128   ohtji->set_GlobalVertexType(GlobalVertex::VTXTYPE::MBD);
0129   ihtji->set_GlobalVertexType(GlobalVertex::VTXTYPE::MBD);
0130   towerjetreco->add_input(emtji);
0131   towerjetreco->add_input(ohtji);
0132   towerjetreco->add_input(ihtji);
0133   towerjetreco->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.2), "AntiKt_r02");
0134   towerjetreco->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.3), "AntiKt_r03");
0135   towerjetreco->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.4), "AntiKt_r04");
0136   towerjetreco->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.5), "AntiKt_r05");
0137   towerjetreco->set_algo_node("ANTIKT");
0138   towerjetreco->set_input_node("TOWER");
0139   towerjetreco->Verbosity(0);
0140   se->registerSubsystem(towerjetreco);
0141
0142   JetCalib *jetCalib02 = new JetCalib("JetCalib02");
0143   jetCalib02->set_InputNode("AntiKt_r02");
0144   jetCalib02->set_OutputNode("AntiKt_r02_calib");
0145   jetCalib02->set_JetRadius(0.2);
0146   jetCalib02->set_ZvrtxNode("GlobalVertexMap");
0147   jetCalib02->set_ApplyZvrtxDependentCalib(true);
0148   jetCalib02->set_ApplyEtaDependentCalib(true);
0149   se->registerSubsystem(jetCalib02);
0150
0151   JetCalib *jetCalib03 = new JetCalib("JetCalib03");
0152   jetCalib03->set_InputNode("AntiKt_r03");
0153   jetCalib03->set_OutputNode("AntiKt_r03_calib");
0154   jetCalib03->set_JetRadius(0.3);
0155   jetCalib03->set_ZvrtxNode("GlobalVertexMap");
0156   jetCalib03->set_ApplyZvrtxDependentCalib(true);
0157   jetCalib03->set_ApplyEtaDependentCalib(true);
0158   se->registerSubsystem(jetCalib03);
0159
0160   JetCalib *jetCalib04 = new JetCalib("JetCalib04");
0161   jetCalib04->set_InputNode("AntiKt_r04");
0162   jetCalib04->set_OutputNode("AntiKt_r04_calib");
0163   jetCalib04->set_JetRadius(0.4);
0164   jetCalib04->set_ZvrtxNode("GlobalVertexMap");
0165   jetCalib04->set_ApplyZvrtxDependentCalib(true);
0166   jetCalib04->set_ApplyEtaDependentCalib(true);
0167   se->registerSubsystem(jetCalib04);
0168
0169   JetCalib *jetCalib05 = new JetCalib("JetCalib05");
0170   jetCalib05->set_InputNode("AntiKt_r05");
0171   jetCalib05->set_OutputNode("AntiKt_r05_calib");
0172   jetCalib05->set_JetRadius(0.5);
0173   jetCalib05->set_ZvrtxNode("GlobalVertexMap");
0174   jetCalib05->set_ApplyZvrtxDependentCalib(true);
0175   jetCalib05->set_ApplyEtaDependentCalib(true);
0176   se->registerSubsystem(jetCalib05);
0177
0178   if(doSim)
0179   {
0180     JetReco *truthjetreco = new JetReco("truthJetReco");
0181     truthjetreco->add_input(new TruthJetInput(Jet::PARTICLE));
0182     truthjetreco->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.2), "AntiKt_Truth_r02");
0183     truthjetreco->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.3), "AntiKt_Truth_r03");
0184     truthjetreco->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.4), "AntiKt_Truth_r04");
0185     truthjetreco->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.5), "AntiKt_Truth_r05");
0186     truthjetreco->set_algo_node("ANTIKT");
0187     truthjetreco->set_input_node("TRUTH");
0188     truthjetreco->Verbosity(0);
0189     se->registerSubsystem(truthjetreco);
0190   }
0191
0192   if(!doSim)
0193   {
0194     TimingCut* tc = new TimingCut("AntiKt_r04","TimingCutModule",false); //arguments: 1. jet node name to use, 2. name for fun4all module, 3. doAbort. The last of these tells the module whether or not to abort events that fail the cuts.
0195     tc->Verbosity(0);
0196     se->registerSubsystem(tc);
0197   }
0198   VandyJetDSTSkimmer *vs = new VandyJetDSTSkimmer("VandyJetDSTSkimmer");
0199   vs->SetRunnumber(runnumber);
0200   vs->SetOutfileName(std::format("{}/VandyDSTs_run2pp_ana521_2025p007_v001_{}-{}-{:06d}_to-{:06d}.root", outDir, sample_name, runnumber,seg, seg+n).c_str());
0201   vs->SetDoSim(doSim);
0202   std::cout<<"Sample name: " <<sample_name<<std::endl;
0203   //if(doSim) vs->SetDoCalib(false);
0204   if(doSim) vs->SetSimSample(sample_name);
0205   std::cout<<vs->GetSimSample() <<std::endl;
0206   se->registerSubsystem(vs);
0207   se->run(0);
0208   se->End();
0209   se->PrintTimer();
0210
0211   delete se;
0212   gSystem->Exit(0);
0213 }