ResonanceJetTaggingOutputs/macro/Fun4All_MDC2reco.C

0001 #include "HFReco.C"
0002
0003 #include <resonancejettaggingoutputs/BuildResonanceJetTaggingTree.h>
0004
0005 #include <resonancejettagging/ResonanceJetTagging.h>
0006
0007 #include <g4main/Fun4AllDstPileupInputManager.h>
0008
0009 #include <G4_Magnet.C>
0010 #include <QA.C>
0011 #include <G4_Global.C>
0012 #include <FROG.h>
0013 #include <decayfinder/DecayFinder.h>
0014 #include <fun4all/Fun4AllDstInputManager.h>
0015
0016 #include <g4eval/SvtxEvaluator.h>
0017 #include <g4eval/SvtxTruthRecoTableEval.h>
0018
0019 #include <caloreco/RawClusterBuilderTopo.h>
0020 #include <particleflowreco/ParticleFlowReco.h>
0021 #include <phool/recoConsts.h>
0022
0023 R__LOAD_LIBRARY(libdecayfinder.so)
0024 R__LOAD_LIBRARY(libcalo_reco.so)
0025 R__LOAD_LIBRARY(libfun4all.so)
0026 R__LOAD_LIBRARY(libparticleflow.so)
0027 R__LOAD_LIBRARY(libresonancejettaggingoutputs.so)
0028
0029 using namespace std;
0030
0031 /****************************/
0032 /*     MDC2 Reco for MDC2   */
0033 /* Cameron Dean, LANL, 2021 */
0034 /*      cdean@bnl.gov       */
0035 /****************************/
0036 /****************************/
0037 /*   Modified for D0Jets    */
0038 /* Antonio Silva, ISU, 2022 */
0039 /*antonio.sphenix@gmail.com */
0040 /****************************/
0041 /****************************/
0042 /*        Contributor       */
0043 /* Jakub Kvapil, LANL, 2023 */
0044 /*   jakub.kvapil@cern.ch   */
0045 /****************************/
0046
0047 void Fun4All_MDC2reco(vector<string> myInputLists = {"condorJob/fileLists/productionFiles-CHARM-dst_tracks-00000.LIST"}, const int nEvents = 10, ResonanceJetTagging::TAG tag = ResonanceJetTagging::TAG::D0)
0048 {
0049   int verbosity = 0;
0050
0051   gSystem->Load("libg4dst.so");
0052   gSystem->Load("libFROG.so");
0053   FROG *fr = new FROG();
0054   recoConsts *rc = recoConsts::instance();
0055   Enable::CDB = true;
0056   // global tag
0057   rc->set_StringFlag("CDB_GLOBALTAG",CDB::global_tag);
0058   // 64 bit timestamp
0059   rc->set_uint64Flag("TIMESTAMP",CDB::timestamp);
0060
0061   std::string particle_name;
0062   switch (tag) {
0063     case ResonanceJetTagging::TAG::D0:
0064       particle_name = HFjets::KFPARTICLE::D0Name;
0065       break;
0066     case ResonanceJetTagging::TAG::D0TOK3PI:
0067       particle_name = HFjets::KFPARTICLE::D0toK3piName;
0068       break;
0069     case ResonanceJetTagging::TAG::DPLUS:
0070       particle_name = HFjets::KFPARTICLE::DplusName;
0071       break;
0072     case ResonanceJetTagging::TAG::LAMBDAC:
0073       particle_name = HFjets::KFPARTICLE::LambdacName;
0074       break;
0075     default:
0076       std::cout<<"ERROR:decay not implemented, ABORTING!";
0077       return Fun4AllReturnCodes::ABORTRUN;
0078       break;
0079   }
0080
0081   //The next set of lines figures out folder revisions, file numbers etc
0082   string outDir = "./";
0083   if (outDir.substr(outDir.size() - 1, 1) != "/") outDir += "/";
0084   outDir += HFjets::Enable::reconstructionName + "_" + particle_name + "/";
0085
0086   string fileNumber = myInputLists[0];
0087   size_t findLastDash = fileNumber.find_last_of("-");
0088   if (findLastDash != string::npos) fileNumber.erase(0, findLastDash + 1);
0089   string remove_this = ".list";
0090   size_t pos = fileNumber.find(remove_this);
0091   if (pos != string::npos) fileNumber.erase(pos, remove_this.length());
0092   string outputFileName = "outputData_" + HFjets::Enable::reconstructionName + "_" + fileNumber + ".root";
0093
0094   string outputRecoDir = outDir + "/inReconstruction/";
0095   string makeDirectory = "mkdir -p " + outputRecoDir;
0096   system(makeDirectory.c_str());
0097   string outputRecoFile = outputRecoDir + outputFileName;
0098
0099   //Create the server
0100   Fun4AllServer *se = Fun4AllServer::instance();
0101   se->Verbosity(verbosity);
0102
0103   //Add all required input files
0104   for (unsigned int i = 0; i < myInputLists.size(); ++i)
0105   {
0106     Fun4AllInputManager *infile = new Fun4AllDstInputManager("DSTin_" + to_string(i));
0107     infile->AddListFile(myInputLists[i]);
0108     se->registerInputManager(infile);
0109   }
0110
0111   // Runs decay finder to trigger on your decay. Useful for signal cleaning
0112   if (HFjets::Enable::runTruthTrigger)
0113   {
0114     DecayFinder *myFinder = new DecayFinder("myFinder");
0115     myFinder->Verbosity(verbosity);
0116     switch (tag) {
0117     case ResonanceJetTagging::TAG::D0:
0118       myFinder->setDecayDescriptor(HFjets::KFPARTICLE::D0toKpiDecayDescriptor);
0119       break;
0120     case ResonanceJetTagging::TAG::D0TOK3PI:
0121       myFinder->setDecayDescriptor(HFjets::KFPARTICLE::D0toK3piDecayDescriptor);
0122       break;
0123     case ResonanceJetTagging::TAG::DPLUS:
0124       myFinder->setDecayDescriptor(HFjets::KFPARTICLE::DplustoK2piDecayDescriptor);
0125       break;
0126     case ResonanceJetTagging::TAG::LAMBDAC:
0127       myFinder->setDecayDescriptor(HFjets::KFPARTICLE::LambdacDecayDescriptor);
0128       break;
0129     default:
0130       std::cout<<"ERROR:decay not implemented, ABORTING!";
0131       return Fun4AllReturnCodes::ABORTRUN;
0132       break;
0133    }
0134     myFinder->saveDST(true);
0135     myFinder->allowPi0(false);
0136     myFinder->allowPhotons(false);
0137     myFinder->triggerOnDecay(true);
0138     myFinder->setPTmin(0.); //Note: sPHENIX min pT is 0.2 GeV for tracking
0139     myFinder->setEtaRange(-10., 10.); //Note: sPHENIX acceptance is |eta| <= 1.1
0140     se->registerSubsystem(myFinder);
0141   }
0142
0143   Global_Reco();
0144
0145   //Now run the actual reconstruction
0146   if (HFjets::KFParticle_Set_Reco(tag) == Fun4AllReturnCodes::ABORTRUN) return;
0147
0148   //Set Calo towers
0149   RawClusterBuilderTopo* ClusterBuilder1 = new RawClusterBuilderTopo("HcalRawClusterBuilderTopo1");
0150   ClusterBuilder1->Verbosity(verbosity);
0151   ClusterBuilder1->set_nodename("TOPOCLUSTER_EMCAL");
0152   ClusterBuilder1->set_enable_HCal(false);
0153   ClusterBuilder1->set_enable_EMCal(true);
0154   ClusterBuilder1->set_noise(0.0025, 0.006, 0.03);
0155   ClusterBuilder1->set_significance(4.0, 2.0, 0.0);
0156   ClusterBuilder1->allow_corner_neighbor(true);
0157   ClusterBuilder1->set_do_split(true);
0158   ClusterBuilder1->set_minE_local_max(1.0, 2.0, 0.5);
0159   ClusterBuilder1->set_R_shower(0.025);
0160   se->registerSubsystem(ClusterBuilder1);
0161
0162   RawClusterBuilderTopo* ClusterBuilder2 = new RawClusterBuilderTopo("HcalRawClusterBuilderTopo2");
0163   ClusterBuilder2->Verbosity(verbosity);
0164   ClusterBuilder2->set_nodename("TOPOCLUSTER_HCAL");
0165   ClusterBuilder2->set_enable_HCal(true);
0166   ClusterBuilder2->set_enable_EMCal(false);
0167   ClusterBuilder2->set_noise(0.0025, 0.006, 0.03);
0168   ClusterBuilder2->set_significance(4.0, 2.0, 0.0);
0169   ClusterBuilder2->allow_corner_neighbor(true);
0170   ClusterBuilder2->set_do_split(true);
0171   ClusterBuilder2->set_minE_local_max(1.0, 2.0, 0.5);
0172   ClusterBuilder2->set_R_shower(0.025);
0173   se->registerSubsystem(ClusterBuilder2);
0174
0175   ParticleFlowReco *pfr = new ParticleFlowReco();
0176   pfr->set_energy_match_Nsigma(1.5);
0177   pfr->Verbosity(verbosity);
0178   se->registerSubsystem(pfr);
0179
0180
0181   std::string jetTagRecoFile = outputRecoDir + particle_name + "JetTree_" + fileNumber + ".root";
0182
0183   ResonanceJetTagging *jetTag = new ResonanceJetTagging(particle_name + "Tagging", tag, particle_name + "_KFParticle_Container");
0184   jetTag->Verbosity(verbosity);
0185   /*
0186   jetTag->setAddTracks(false);
0187   jetTag->setAddEMCalClusters(false);
0188   jetTag->setTrackPtAcc(0.2, 9999.);
0189   jetTag->setTrackEtaAcc(-1.1, 1.1);
0190   jetTag->setEMCalClusterPtAcc(0.3, 9999.);
0191   jetTag->setEMCalClusterEtaAcc(-1.1, 1.1);
0192   jetTag->setHCalClusterPtAcc(0.3, 9999.);
0193   jetTag->setHCalClusterEtaAcc(-1.1, 1.1);
0194   */
0195   jetTag->setParticleFlowEtaAcc(-1.7, 1.7); // -1.1 1.1
0196   jetTag->setJetParameters(0.4, ResonanceJetTagging::ALGO::ANTIKT, ResonanceJetTagging::RECOMB::E_SCHEME);
0197   jetTag->setJetContainerName(particle_name + "Jets");
0198   jetTag->setDoTruth(true);
0199   se->registerSubsystem(jetTag);
0200
0201   BuildResonanceJetTaggingTree *buildTree = new BuildResonanceJetTaggingTree(particle_name + "JetTree", jetTagRecoFile, tag);
0202   buildTree->setDoTruth(true);
0203   buildTree->setTagContainerName(particle_name + "_KFParticle_Container");
0204   buildTree->setJetContainerName(particle_name + "Jets_Jet_Container");
0205   buildTree->setTruthJetContainerName(particle_name + "Jets_Truth_Jet_Container");
0206   se->registerSubsystem(buildTree);
0207
0208   se->run(nEvents);
0209   se->End();
0210
0211   ifstream file(outputRecoFile.c_str());
0212   if (file.good())
0213   {
0214     std::string moveOutput = "mv " + outputRecoFile + " " + outDir;
0215     system(moveOutput.c_str());
0216   }
0217
0218   ifstream jetTagfile(jetTagRecoFile.c_str());
0219   if(jetTagfile.good())
0220   {
0221     std::string moveOutput = "mv " + jetTagRecoFile + " " + outDir;
0222     system(moveOutput.c_str());
0223   }
0224
0225   std::cout << "All done" << std::endl;
0226   delete se;
0227   gSystem->Exit(0);
0228
0229   return;
0230 }