HF-Particle/KFParticle_sPHENIX/Fun4All_MDC2reco.C

0001 #include "HFReco.C"
0002
0003 #include <g4main/Fun4AllDstPileupInputManager.h>
0004
0005 #include <G4_Global.C>
0006 #include <G4_Magnet.C>
0007 #include <G4_ActsGeom.C>
0008 #include <G4_TrkrSimulation.C>
0009
0010 #include <Trkr_Clustering.C>
0011 #include <Trkr_Reco.C>
0012 #include <Trkr_RecoInit.C>
0013
0014 #include <FROG.h>
0015 #include <decayfinder/DecayFinder.h>
0016 #include <hftrackefficiency/HFTrackEfficiency.h>
0017
0018 #include <fun4all/Fun4AllDstInputManager.h>
0019
0020 #include <g4eval/SvtxEvaluator.h>
0021 #include <g4eval/SvtxTruthRecoTableEval.h>
0022
0023 #include <ffamodules/CDBInterface.h>
0024 #include <phool/recoConsts.h>
0025 #include <globalvertex/GlobalVertexReco.h>
0026
0027 R__LOAD_LIBRARY(libfun4all.so)
0028 R__LOAD_LIBRARY(libdecayfinder.so)
0029 R__LOAD_LIBRARY(libhftrackefficiency.so)
0030 R__LOAD_LIBRARY(libg4vertex.so)
0031 R__LOAD_LIBRARY(libglobalvertex.so)
0032
0033 using namespace std;
0034 using namespace HeavyFlavorReco;
0035
0036 /****************************/
0037 /*     MDC2 Reco for MDC2   */
0038 /* Cameron Dean, LANL, 2021 */
0039 /*      cdean@bnl.gov       */
0040 /****************************/
0041
0042 void Fun4All_MDC2reco(vector<string> myInputLists = {"condorJob/fileLists/productionFiles-CHARM-dst_tracks-00000.list"}, const int nEvents = 0)
0043 {
0044   bool runDemo = true; //Run the pure Bs2JpsiKs0 reco
0045
0046   gSystem->Load("libg4dst.so");
0047   gSystem->Load("libFROG.so");
0048   FROG *fr = new FROG();
0049
0050   string outDir = "./";
0051   if (!runDemo)
0052   {
0053     //The next set of lines figures out folder revisions, file numbers etc
0054     if (outDir.substr(outDir.size() - 1, 1) != "/") outDir += "/";
0055     outDir += reconstructionName + "/";
0056
0057     string fileNumber = myInputLists[0];
0058     size_t findLastDash = fileNumber.find_last_of("-");
0059     if (findLastDash != string::npos) fileNumber.erase(0, findLastDash + 1);
0060     string remove_this = ".list";
0061     size_t pos = fileNumber.find(remove_this);
0062     if (pos != string::npos) fileNumber.erase(pos, remove_this.length());
0063     string outputFileName = "outputFile_" + reconstructionName + "_" + fileNumber + ".root";
0064     string outputHFEffFileName = "outputHFEff_" + reconstructionName + "_" + fileNumber + ".root";
0065
0066     string outputRecoDir = outDir + "inReconstruction/";
0067     string makeDirectory = "mkdir -p " + outputRecoDir;
0068     system(makeDirectory.c_str());
0069     outputRecoFile = outputRecoDir + outputFileName;
0070     if (runTrackEff) outputHFEffFile = outputRecoDir + outputHFEffFileName;
0071   }
0072
0073   recoConsts *rc = recoConsts::instance();
0074   // global tag
0075   rc->set_StringFlag("CDB_GLOBALTAG", "MDC2");
0076   // 64 bit timestamp
0077   rc->set_uint64Flag("TIMESTAMP", 6);
0078
0079   //Create the server
0080   Fun4AllServer *se = Fun4AllServer::instance();
0081   se->Verbosity(1);
0082   if (runDemo)
0083   {
0084     Fun4AllInputManager *infile = new Fun4AllDstInputManager("DSTin");
0085     infile->AddListFile("demo.list");
0086     se->registerInputManager(infile);
0087
0088     myDemoReco();
0089   }
0090   else
0091   {
0092     //Add all required input files
0093     for (unsigned int i = 0; i < myInputLists.size(); ++i)
0094     {
0095       Fun4AllInputManager *infile = new Fun4AllDstInputManager("DSTin_" + to_string(i));
0096       infile->AddListFile(myInputLists[i]);
0097       se->registerInputManager(infile);
0098     }
0099
0100     // Runs decay finder to trigger on your decay. Useful for signal cleaning
0101     if (runTruthTrigger)
0102     {
0103       DecayFinder *myFinder = new DecayFinder("myFinder");
0104       myFinder->Verbosity(INT_MAX);
0105       myFinder->setDecayDescriptor(decayDescriptor);
0106       myFinder->saveDST(1);
0107       myFinder->allowPi0(0);
0108       myFinder->allowPhotons(1);
0109       myFinder->triggerOnDecay(1);
0110       myFinder->setPTmin(0.2); //Note: sPHENIX min pT is 0.2 GeV for tracking
0111       myFinder->setEtaRange(-1.1, 1.1); //Note: sPHENIX acceptance is |eta| <= 1.1
0112       se->registerSubsystem(myFinder);
0113     }
0114
0115     //Run the tracking if not already done
0116     if (runTracking)
0117     {
0118       Enable::MICROMEGAS=true;
0119
0120       G4MAGNET::magfield_rescale = 1.;
0121       MagnetInit();
0122       MagnetFieldInit();
0123
0124       Mvtx_Cells();
0125       Intt_Cells();
0126       TPC_Cells();
0127       Micromegas_Cells();
0128
0129       TrackingInit();
0130
0131       Mvtx_Clustering();
0132       Intt_Clustering();
0133       TPC_Clustering();
0134       Micromegas_Clustering();
0135
0136       Tracking_Reco();
0137
0138       Global_Reco();
0139     }
0140
0141     if (buildTruthTable)
0142     {
0143       SvtxTruthRecoTableEval *tables = new SvtxTruthRecoTableEval();
0144       tables->Verbosity(VERBOSITY);
0145       se->registerSubsystem(tables);
0146     }
0147
0148     if (runTrackEff)
0149     {
0150       HFTrackEfficiency *myTrackEff = new HFTrackEfficiency("myTrackEff");
0151       myTrackEff->Verbosity(1);
0152       myTrackEff->setDFNodeName("myFinder");
0153       myTrackEff->triggerOnDecay(1);
0154       myTrackEff->writeSelectedTrackMap(true);
0155       myTrackEff->writeOutputFile(true);
0156       myTrackEff->setOutputFileName(outputHFEffFile);
0157       se->registerSubsystem(myTrackEff);
0158     }
0159
0160     //Now run the actual reconstruction
0161     myHeavyFlavorReco();
0162   }
0163
0164   se->run(nEvents);
0165   se->End();
0166
0167   if (!runDemo)
0168   {
0169     ifstream file(outputRecoFile.c_str());
0170     if (file.good())
0171     {
0172       string moveOutput = "mv " + outputRecoFile + " " + outDir;
0173       system(moveOutput.c_str());
0174     }
0175
0176     if (runTrackEff)
0177     {
0178       ifstream evalfile(outputHFEffFile.c_str());
0179       if (evalfile.good())
0180       {
0181         string moveOutput = "mv " + outputHFEffFile + " " + outDir;
0182         system(moveOutput.c_str());
0183       }
0184     }
0185   }
0186
0187   std::cout << "All done" << std::endl;
0188   delete se;
0189   gSystem->Exit(0);
0190
0191   return;
0192 }