diff --git a/JetProduction/Fun4All_JetSkimmedProductionYear3.C b/JetProduction/Fun4All_JetSkimmedProductionYear3.C index 69bca39e2..116748c7d 100644 --- a/JetProduction/Fun4All_JetSkimmedProductionYear3.C +++ b/JetProduction/Fun4All_JetSkimmedProductionYear3.C @@ -35,6 +35,45 @@ #include #include +#include +#include + +std::string ResolveInputFile(const std::string &baseDir, const std::string &fname) +{ + if (fname.find('/') != std::string::npos) + { + return fname; + } + + std::string candidate = baseDir + "/" + fname; + if (!gSystem->AccessPathName(candidate.c_str())) // returns kFALSE if file exists + { + return candidate; + } + std::string pattern = baseDir + "/run_*/" + fname; + + // build the shell command + std::string cmd = "ls " + pattern + " 2>/dev/null | head -n1"; + + // GetFromPipe returns a TString + TString t = gSystem->GetFromPipe(cmd.c_str()); + std::string result = std::string(t.Data()); + + // strip trailing newline(s) + while (!result.empty() && (result.back() == '\n' || result.back() == '\r')) + { + result.pop_back(); + } + if (!result.empty()) + { + return result; // found a match + } + std::cerr << "ResolveInputFile: could not find " << fname + << " under " << baseDir << std::endl; + return fname; +} + + R__LOAD_LIBRARY(libfun4all.so) R__LOAD_LIBRARY(libfun4allraw.so) R__LOAD_LIBRARY(libcalo_reco.so) @@ -47,6 +86,19 @@ R__LOAD_LIBRARY(libcalovalid.so) R__LOAD_LIBRARY(libjetbackground.so) R__LOAD_LIBRARY(libJetDSTSkimmer.so) +void Fun4All_JetSkimmedProductionYear3( + int nEvents = 1000, + const std::string &fname = + //"DST_CALOFITTING_run3auau_new_newcdbtag_v008-00075142-00000.root", + "DST_CALOFITTING_run3pp_new_newcdbtag_v008-00079151-00010.root", + const std::string &outfile_low = "DST_JETCALO-00000000-000000.root", + const std::string &outfile_high = "DST_Jet-00000000-000000.root", + const std::string &outfile_hist = "HIST_JETQA-00000000-000000.root", + const std::string &dbtag = "newcdbtag", + const std::string &basedir = + //"/sphenix/lustre01/sphnxpro/production2/run3auau/physics/calofitting/new_newcdbtag_v008") + "/sphenix/lustre01/sphnxpro/production2/run3pp/physics/calofitting/new_newcdbtag_v008") +/* void Fun4All_JetSkimmedProductionYear3( int nEvents = 1000, const std::string &fname = @@ -55,6 +107,7 @@ void Fun4All_JetSkimmedProductionYear3( const std::string &outfile_high = "DST_Jet-00000000-000000.root", const std::string &outfile_hist = "HIST_JETQA-00000000-000000.root", const std::string &dbtag = "newcdbtag") +*/ { Fun4AllServer *se = Fun4AllServer::instance(); recoConsts::instance()->set_IntFlag("PHOOL_VERBOSITY", 0); @@ -90,7 +143,7 @@ void Fun4All_JetSkimmedProductionYear3( Enable::QA = true; // Au+Au mode - HIJETS::is_pp = false; + HIJETS::is_pp = true; // QA options JetQA::HasTracks = false; @@ -106,12 +159,11 @@ void Fun4All_JetSkimmedProductionYear3( { std::cout << ">>> Running Centrality()" << std::endl; Centrality(); + } - // HIJetReco macro - std::cout << ">>> Running HIJetReco()" << std::endl; - HIJetReco(); - + HIJetReco(); + // Jet DST skimmer with map-based threshold std::cout << ">>> Configuring JetDSTSkimmer" << std::endl; JetDSTSkimmer *jetDSTSkimmer = new JetDSTSkimmer(); @@ -119,15 +171,15 @@ void Fun4All_JetSkimmedProductionYear3( // Jet thresholds: map from node name [set min jet pT] std::map jetNodePts; - jetNodePts["AntiKt_Tower_r04_Sub1"] = 7.0F; - jetDSTSkimmer->SetJetNodeThresholds(jetNodePts); + //jetNodePts["AntiKt_Tower_r04_Sub1"] = 7.0F; + //jetDSTSkimmer->SetJetNodeThresholds(jetNodePts); // Cluster thresholds: map from cluster node [set min cluster pT] std::map clusterNodePts; - clusterNodePts["CLUSTERINFO_CEMC"] = 5.0F; - jetDSTSkimmer->SetClusterNodeThresholds(clusterNodePts); + // clusterNodePts["CLUSTERINFO_CEMC"] = 5.0F; + //jetDSTSkimmer->SetClusterNodeThresholds(clusterNodePts); - jetDSTSkimmer->Verbosity(0); + jetDSTSkimmer->Verbosity(10); se->registerSubsystem(jetDSTSkimmer); std::cout << ">>> JetDSTSkimmer registered" << std::endl; @@ -156,7 +208,17 @@ void Fun4All_JetSkimmedProductionYear3( Fun4AllInputManager *In = new Fun4AllDstInputManager("in"); In->Verbosity(1); - In->AddFile(fname); + + const std::string inFile = ResolveInputFile(basedir, fname); + std::cout << ">>> Input file resolved to: " << inFile << std::endl; + + if (gSystem->AccessPathName(inFile.c_str())) + { + std::cerr << "FATAL: cannot access input file: " << inFile << std::endl; + gSystem->Exit(1); + } + + In->AddFile(inFile); se->registerInputManager(In); Fun4AllDstOutputManager *outlower = diff --git a/JetProduction/Fun4All_JetSkimmedProductionYear3_pp.C b/JetProduction/Fun4All_JetSkimmedProductionYear3_pp.C new file mode 100644 index 000000000..e3ca16fcb --- /dev/null +++ b/JetProduction/Fun4All_JetSkimmedProductionYear3_pp.C @@ -0,0 +1,217 @@ +#ifndef FUN4ALL_JETSKIMMEDPRODUCTIONYEAR3_C +#define FUN4ALL_JETSKIMMEDPRODUCTIONYEAR3_C + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include + + +R__LOAD_LIBRARY(libfun4all.so) +R__LOAD_LIBRARY(libfun4allraw.so) +R__LOAD_LIBRARY(libcalo_reco.so) +R__LOAD_LIBRARY(libcalotrigger.so) +R__LOAD_LIBRARY(libcentrality.so) +R__LOAD_LIBRARY(libffamodules.so) +R__LOAD_LIBRARY(libmbd.so) +R__LOAD_LIBRARY(libglobalvertex.so) +R__LOAD_LIBRARY(libcalovalid.so) +R__LOAD_LIBRARY(libjetbackground.so) +R__LOAD_LIBRARY(libJetDSTSkimmer.so) + +void Fun4All_JetSkimmedProductionYear3_pp( + int nEvents = 1000, + const std::string &fname = + "DST_CALOFITTING_run3pp_new_newcdbtag_v008-00081000-00000.root", + const std::string &outfile_low = "DST_JETCALO-00000000-000000.root", + const std::string &outfile_high = "DST_Jet-00000000-000000.root", + const std::string &outfile_hist = "HIST_JETQA-00000000-000000.root", + const std::string &dbtag = "newcdbtag") +{ + Fun4AllServer *se = Fun4AllServer::instance(); + recoConsts::instance()->set_IntFlag("PHOOL_VERBOSITY", 0); + se->Verbosity(0); + + recoConsts *rc = recoConsts::instance(); + + std::pair runseg = Fun4AllUtils::GetRunSegment(fname); + int runnumber = runseg.first; + + // conditions DB flags and timestamp + rc->set_StringFlag("CDB_GLOBALTAG", dbtag); + rc->set_uint64Flag("TIMESTAMP", runnumber); + CDBInterface::instance()->Verbosity(0); + + FlagHandler *flag = new FlagHandler(); + se->registerSubsystem(flag); + + // MBD/BBC Reconstruction + MbdReco *mbdreco = new MbdReco(); + se->registerSubsystem(mbdreco); + + // Official vertex storage + GlobalVertexReco *gvertex = new GlobalVertexReco(); + se->registerSubsystem(gvertex); + + ///////////////////////////////////////////////////// + // Set status of CALO towers, Calibrate towers, Cluster + std::cout << ">>> Processing Calo Calibration" << std::endl; + Process_Calo_Calib(); + + // Jet reco related flags + Enable::QA = true; + + // is Au+Au mode? + HIJETS::is_pp = true; + + //is Au+Au mode or o+o (true)? + HIJETS::is_oo = false; + // QA options + JetQA::HasTracks = false; + JetQA::DoInclusive = true; + JetQA::DoTriggered = true; + JetQA::RestrictPtToTrig = false; + JetQA::RestrictEtaByR = true; + JetQA::DoPP = HIJETS::is_pp; + JetQA::DoOO = HIJETS::is_oo; + JetQA::UseBkgdSub = true; + JetQA::HasCalos = true; + + if (!HIJETS::is_pp) + { + std::cout << ">>> Running Centrality()" << std::endl; + Centrality(); + + } + + HIJetReco(); + + // Jet DST skimmer with map-based threshold + std::cout << ">>> Configuring JetDSTSkimmer" << std::endl; + JetDSTSkimmer *jetDSTSkimmer = new JetDSTSkimmer(); + + // Jet thresholds: map from node name [set min jet pT] + std::map jetNodePts; + + //jetNodePts["AntiKt_Tower_r04_Sub1"] = 7.0F; + //jetDSTSkimmer->SetJetNodeThresholds(jetNodePts); + + // Cluster thresholds: map from cluster node [set min cluster pT] + std::map clusterNodePts; + // clusterNodePts["CLUSTERINFO_CEMC"] = 5.0F; + //jetDSTSkimmer->SetClusterNodeThresholds(clusterNodePts); + + jetDSTSkimmer->Verbosity(10); + se->registerSubsystem(jetDSTSkimmer); + std::cout << ">>> JetDSTSkimmer registered" << std::endl; + + // Filter out beam-background events + BeamBackgroundFilterAndQA *filter = + new BeamBackgroundFilterAndQA("BeamBackgroundFilterAndQA"); + filter->Verbosity(std::max(Enable::QA_VERBOSITY, Enable::JETQA_VERBOSITY)); + filter->SetConfig( + { + .debug = false, + .doQA = Enable::QA, + .doEvtAbort = false, + .filtersToApply = {"StreakSideband"} + }); + se->registerSubsystem(filter); + + // Register modules necessary for QA + if (Enable::QA) + { + std::cout << ">>> Enabling QA: DoRhoCalculation + Jet_QA" << std::endl; + DoRhoCalculation(); + Jet_QA(); + } + + std::cout << ">>> Setting up input and output managers" << std::endl; + + Fun4AllInputManager *In = new Fun4AllDstInputManager("in"); + In->Verbosity(1); + + + In->AddFile(fname); + se->registerInputManager(In); + + Fun4AllDstOutputManager *outlower = + new Fun4AllDstOutputManager("DSTOUTLOW", outfile_low); + outlower->AddNode("Sync"); + outlower->AddNode("EventHeader"); + outlower->AddNode("GL1Packet"); + outlower->AddNode("TOWERS_HCALIN"); + outlower->AddNode("TOWERS_HCALOUT"); + outlower->AddNode("TOWERS_SEPD"); + outlower->AddNode("TOWERS_ZDC"); + outlower->AddNode("MbdOut"); + outlower->AddNode("MbdPmtContainer"); + outlower->AddNode("MBDPackets"); + outlower->AddNode("TriggerRunInfo"); + se->registerOutputManager(outlower); + + Fun4AllDstOutputManager *outhigher = + new Fun4AllDstOutputManager("DSTOUTHIGH", outfile_high); + outhigher->StripNode("TOWERINFO_CALIB_HCALIN"); + outhigher->StripNode("TOWERS_HCALIN"); + outhigher->StripNode("TOWERINFO_CALIB_HCALOUT"); + outhigher->StripNode("TOWERS_HCALOUT"); + outhigher->StripNode("TOWERINFO_CALIB_CEMC"); + outhigher->StripNode("TOWERS_CEMC"); + outhigher->StripNode("TOWERS_SEPD"); + outhigher->StripNode("TOWERINFO_CALIB_ZDC"); + outhigher->StripNode("TOWERS_ZDC"); + outhigher->StripNode("MBDPackets"); + outhigher->StripNode("MbdOut"); + outhigher->StripNode("MbdPmtContainer"); + se->registerOutputManager(outhigher); + + se->run(nEvents); + se->End(); + + if (Enable::QA) + { + std::cout << ">>> Saving QA histograms to " << outfile_hist << std::endl; + QAHistManagerDef::saveQARootFile(outfile_hist); + } + + CDBInterface::instance()->Print(); // print used DB files + se->PrintTimer(); + delete se; + std::cout << "All done!" << std::endl; + gSystem->Exit(0); +} + +#endif diff --git a/JetProduction/HIJetReco.C b/JetProduction/HIJetReco.C new file mode 100644 index 000000000..e17ca3840 --- /dev/null +++ b/JetProduction/HIJetReco.C @@ -0,0 +1,420 @@ +#ifndef MACRO_HIJETRECO_C +#define MACRO_HIJETRECO_C + +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include + +#include + +R__LOAD_LIBRARY(libg4jets.so) +R__LOAD_LIBRARY(libglobalvertex.so) +R__LOAD_LIBRARY(libjetbackground.so) +R__LOAD_LIBRARY(libjetbase.so) +R__LOAD_LIBRARY(libparticleflow.so) +R__LOAD_LIBRARY(libeventplaneinfo.so) + +// ---------------------------------------------------------------------------- +//! General options for background-subtracted jet reconstruction +// ---------------------------------------------------------------------------- +namespace Enable +{ + bool HIJETS = false; ///< do HI jet reconstruction + int HIJETS_VERBOSITY = 0; ///< verbosity + bool HIJETS_MC = false; ///< is simulation + bool HIJETS_TRUTH = false; ///< make truth jets + bool HIJETS_TOWER = true; ///< make tower jets + bool HIJETS_TRACK = false; ///< make track jets + bool HIJETS_PFLOW = false; ///< make particle flow jets +} // namespace Enable + +// ---------------------------------------------------------------------------- +//! Options specific to background-subtracted jet reconstruction +// ---------------------------------------------------------------------------- +namespace HIJETS +{ + // do_flow = 0 --noflow + // do_flow = 1 --psi2 derived from calo + // do_flow = 2 --psi2 derived from HIJING + // do_flow = 3 --psi2 derived from sEPD + int do_flow = 0; + + ///! do constituent subtraction + bool do_CS = false; + + ///! turn on/off functionality only relevant for + ///! nucleon collisions + bool is_pp = true; + + bool is_oo = true; + ///! sets prefix of nodes to use as tower jet + ///! input + std::string tower_prefix = "TOWERINFO_CALIB"; + + ///! turn on/off vertex specification + bool do_vertex_type = true; + + ///! if specifying vertex, set vertex + ///! to use to this one + GlobalVertex::VTXTYPE vertex_type = GlobalVertex::MBD; + + ///! Base fastjet options to use. Note that the + ///! resolution parameter will be overwritten + ///! to R = 0.2, 0.3, 0.4, and 0.5 + FastJetOptions fj_opts({Jet::ANTIKT, JET_R, 0.4, VERBOSITY, static_cast(Enable::HIJETS_VERBOSITY)}); + + ///! sets jet node name + std::string jet_node = "ANTIKT"; + + ///! sets prefix of nodes to store jets + std::string algo_prefix = "AntiKt"; + + ///! sets the embedding g4 flag used to determine which particles are classified as 'truth' in MakeHITruthJets + ///! if you don't know what this is, you probably don't need to change it + ///! 0 = all particles (use for all g4particles !not useful for HIJING+Pythia samples) + ///! 1 = pythia/herwig particles only (use for pythia/herwig jets in HIJING+Pythia samples) + ///! 2 = pythia particles from the HIJING+Pythia samples (use for pythia jets in HIJING+Pythia samples) + ///! negative values are typically background particles (HIJING particles) + int embedding_flag = 1; + + ///! enumerates reconstructed resolution + ///! parameters + enum Res + { + R02 = 0, + R03 = 1, + R04 = 2, + R05 = 3 + }; + + // -------------------------------------------------------------------------- + //! Helper method to generate releveant FastJet algorithms + // -------------------------------------------------------------------------- + FastJetAlgoSub *GetFJAlgo(const float reso) + { + // grab current options & update + // reso parameter + FastJetOptions opts = fj_opts; + opts.jet_R = reso; + + // create new algorithm + return new FastJetAlgoSub(opts); + + } // end 'GetFJAlgo()' +} // namespace HIJETS + +// ---------------------------------------------------------------------------- +//! Make jets out of appropriate truth particles +// ---------------------------------------------------------------------------- +void MakeHITruthJets() +{ + // set verbosity + int verbosity = std::max(Enable::VERBOSITY, Enable::HIJETS_VERBOSITY); + + //--------------- + // Fun4All server + //--------------- + Fun4AllServer *se = Fun4AllServer::instance(); + + // if making track jets, make truth jets out of only charged particles + if (Enable::HIJETS_TRACK) + { + // configure truth jet input for charged particles + TruthJetInput *ctji = new TruthJetInput(Jet::SRC::CHARGED_PARTICLE); + ctji->add_embedding_flag(HIJETS::embedding_flag); // changes depending on signal vs. embedded + + // book jet reconstruction on chargedparticles + JetReco *chargedtruthjetreco = new JetReco(); + chargedtruthjetreco->add_input(ctji); + chargedtruthjetreco->add_algo(HIJETS::GetFJAlgo(0.2), HIJETS::algo_prefix + "_ChargedTruth_r02"); + chargedtruthjetreco->add_algo(HIJETS::GetFJAlgo(0.3), HIJETS::algo_prefix + "_ChargedTruth_r03"); + chargedtruthjetreco->add_algo(HIJETS::GetFJAlgo(0.4), HIJETS::algo_prefix + "_ChargedTruth_r04"); + chargedtruthjetreco->add_algo(HIJETS::GetFJAlgo(0.5), HIJETS::algo_prefix + "_ChargedTruth_r05"); + chargedtruthjetreco->set_algo_node(HIJETS::jet_node); + chargedtruthjetreco->set_input_node("TRUTH"); + chargedtruthjetreco->Verbosity(verbosity); + se->registerSubsystem(chargedtruthjetreco); + } + + // if making tower or pflow jets, make truth jets out of all particles + if (Enable::HIJETS_TOWER || Enable::HIJETS_PFLOW) + { + // configure truth jet input for all particles + TruthJetInput *tji = new TruthJetInput(Jet::PARTICLE); + tji->add_embedding_flag(HIJETS::embedding_flag); // changes depending on signal vs. embedded + + // book jet reconstruction on all particles + JetReco *truthjetreco = new JetReco(); + truthjetreco->add_input(tji); + truthjetreco->add_algo(HIJETS::GetFJAlgo(0.2), HIJETS::algo_prefix + "_Truth_r02"); + truthjetreco->add_algo(HIJETS::GetFJAlgo(0.3), HIJETS::algo_prefix + "_Truth_r03"); + truthjetreco->add_algo(HIJETS::GetFJAlgo(0.4), HIJETS::algo_prefix + "_Truth_r04"); + truthjetreco->add_algo(HIJETS::GetFJAlgo(0.5), HIJETS::algo_prefix + "_Truth_r05"); + truthjetreco->set_algo_node(HIJETS::jet_node); + truthjetreco->set_input_node("TRUTH"); + truthjetreco->Verbosity(verbosity); + se->registerSubsystem(truthjetreco); + } + + // exit back to HIJetReco() + return; +} + +// ---------------------------------------------------------------------------- +//! Make jets out of subtracted towers +// ---------------------------------------------------------------------------- +void MakeHITowerJets() +{ + int verbosity = std::max(Enable::VERBOSITY, Enable::HIJETS_VERBOSITY); + + //--------------- + // Fun4All server + //--------------- + Fun4AllServer *se = Fun4AllServer::instance(); + + if (HIJETS::do_flow == 3) + { + EventPlaneReco *epreco = new EventPlaneReco(); + se->registerSubsystem(epreco); + } + + RetowerCEMC *rcemc = new RetowerCEMC(); + rcemc->Verbosity(verbosity); + rcemc->set_towerinfo(true); + rcemc->set_frac_cut(0.5); // fraction of retower that must be masked to mask the full retower + rcemc->set_towerNodePrefix(HIJETS::tower_prefix); + se->registerSubsystem(rcemc); + + JetReco *towerjetreco = new JetReco(); + TowerJetInput *incemc = new TowerJetInput(Jet::CEMC_TOWERINFO_RETOWER, HIJETS::tower_prefix); + TowerJetInput *inihcal = new TowerJetInput(Jet::HCALIN_TOWERINFO, HIJETS::tower_prefix); + TowerJetInput *inohcal = new TowerJetInput(Jet::HCALOUT_TOWERINFO, HIJETS::tower_prefix); + if (HIJETS::do_vertex_type) + { + incemc->set_GlobalVertexType(HIJETS::vertex_type); + inihcal->set_GlobalVertexType(HIJETS::vertex_type); + inohcal->set_GlobalVertexType(HIJETS::vertex_type); + } + towerjetreco->add_input(incemc); + towerjetreco->add_input(inihcal); + towerjetreco->add_input(inohcal); + towerjetreco->add_algo(HIJETS::GetFJAlgo(0.2), HIJETS::algo_prefix + "_TowerInfo_HIRecoSeedsRaw_r02"); + towerjetreco->set_algo_node(HIJETS::jet_node); + towerjetreco->set_input_node("TOWER"); + towerjetreco->Verbosity(verbosity); + se->registerSubsystem(towerjetreco); + + DetermineTowerBackground *dtb = new DetermineTowerBackground(); + dtb->SetBackgroundOutputName("TowerInfoBackground_Sub1"); + dtb->SetFlow(HIJETS::do_flow); + dtb->SetSeedType(0); + dtb->SetSeedJetD(3); + // dtb->set_towerinfo(true); + dtb->Verbosity(verbosity); + dtb->set_towerNodePrefix(HIJETS::tower_prefix); + se->registerSubsystem(dtb); + + CopyAndSubtractJets *casj = new CopyAndSubtractJets(); + casj->SetFlowModulation(HIJETS::do_flow); + casj->Verbosity(verbosity); + casj->set_towerinfo(true); + casj->set_towerNodePrefix(HIJETS::tower_prefix); + se->registerSubsystem(casj); + + DetermineTowerBackground *dtb2 = new DetermineTowerBackground(); + dtb2->SetBackgroundOutputName("TowerInfoBackground_Sub2"); + dtb2->SetFlow(HIJETS::do_flow); + dtb2->SetSeedType(1); + dtb2->SetSeedJetPt(7); + dtb2->Verbosity(verbosity); + // dtb2->set_towerinfo(true); + dtb2->set_towerNodePrefix(HIJETS::tower_prefix); + se->registerSubsystem(dtb2); + + SubtractTowers *st = new SubtractTowers(); + st->SetFlowModulation(HIJETS::do_flow); + st->Verbosity(verbosity); + st->set_towerinfo(true); + st->set_towerNodePrefix(HIJETS::tower_prefix); + se->registerSubsystem(st); + + towerjetreco = new JetReco(); + incemc = new TowerJetInput(Jet::CEMC_TOWERINFO_SUB1, HIJETS::tower_prefix); + inihcal = new TowerJetInput(Jet::HCALIN_TOWERINFO_SUB1, HIJETS::tower_prefix); + inohcal = new TowerJetInput(Jet::HCALOUT_TOWERINFO_SUB1, HIJETS::tower_prefix); + if (HIJETS::do_vertex_type) + { + incemc->set_GlobalVertexType(HIJETS::vertex_type); + inihcal->set_GlobalVertexType(HIJETS::vertex_type); + inohcal->set_GlobalVertexType(HIJETS::vertex_type); + } + towerjetreco->add_input(incemc); + towerjetreco->add_input(inihcal); + towerjetreco->add_input(inohcal); + towerjetreco->add_algo(HIJETS::GetFJAlgo(0.2), HIJETS::algo_prefix + "_Tower_r02_Sub1"); + towerjetreco->add_algo(HIJETS::GetFJAlgo(0.3), HIJETS::algo_prefix + "_Tower_r03_Sub1"); + towerjetreco->add_algo(HIJETS::GetFJAlgo(0.4), HIJETS::algo_prefix + "_Tower_r04_Sub1"); + towerjetreco->add_algo(HIJETS::GetFJAlgo(0.5), HIJETS::algo_prefix + "_Tower_r05_Sub1"); + towerjetreco->set_algo_node(HIJETS::jet_node); + towerjetreco->set_input_node("TOWER"); + towerjetreco->Verbosity(verbosity); + se->registerSubsystem(towerjetreco); + + return; +} + +// ---------------------------------------------------------------------------- +//! Make jets out of tracks with background subtraction +// ---------------------------------------------------------------------------- +void MakeHITrackJets() +{ + // set verbosity + int verbosity = std::max(Enable::VERBOSITY, Enable::HIJETS_VERBOSITY); + + //--------------- + // Fun4All server + //--------------- + Fun4AllServer *se = Fun4AllServer::instance(); + + // emit warning: background sub will be added later + std::cerr << "WARNING: Background subtraction for track jets is still in development!\n" + << " If you want to do jet reco without background subtraction, please\n" + << " use NoBkgdSubJetReco()" + << std::endl; + + // book jet reconstruction routines on tracks + JetReco *trackjetreco = new JetReco(); + trackjetreco->add_input(new TrackJetInput(Jet::SRC::TRACK)); + trackjetreco->add_algo(HIJETS::GetFJAlgo(0.2), HIJETS::algo_prefix + "_Track_r02"); + trackjetreco->add_algo(HIJETS::GetFJAlgo(0.3), HIJETS::algo_prefix + "_Track_r03"); + trackjetreco->add_algo(HIJETS::GetFJAlgo(0.4), HIJETS::algo_prefix + "_Track_r04"); + trackjetreco->add_algo(HIJETS::GetFJAlgo(0.5), HIJETS::algo_prefix + "_Track_r05"); + trackjetreco->set_algo_node(HIJETS::jet_node); + trackjetreco->set_input_node("TRACK"); + trackjetreco->Verbosity(verbosity); + se->registerSubsystem(trackjetreco); + + // exit back to HIJetReco() + return; +} + +// ---------------------------------------------------------------------------- +//! Make jets out of particle-flow elements with background subtraction +// ---------------------------------------------------------------------------- +void MakeHIPFlowJets() +{ + // set verbosity + int verbosity = std::max(Enable::VERBOSITY, Enable::HIJETS_VERBOSITY); + + //--------------- + // Fun4All server + //--------------- + Fun4AllServer *se = Fun4AllServer::instance(); + + // emit warning: background sub will be added later + std::cerr << "WARNING: Background subtraction for particle-flow jets is still in development!\n" + << " If you want to do jet reco without background subtraction, please\n" + << " use NoBkgdSubJetReco.C" + << std::endl; + + // book jet reconstruction routines on pflow elements + JetReco *pflowjetreco = new JetReco(); + pflowjetreco->add_input(new ParticleFlowJetInput()); + pflowjetreco->add_algo(HIJETS::GetFJAlgo(0.2), HIJETS::algo_prefix + "_ParticleFlow_r02"); + pflowjetreco->add_algo(HIJETS::GetFJAlgo(0.3), HIJETS::algo_prefix + "_ParticleFlow_r03"); + pflowjetreco->add_algo(HIJETS::GetFJAlgo(0.4), HIJETS::algo_prefix + "_ParticleFlow_r04"); + pflowjetreco->add_algo(HIJETS::GetFJAlgo(0.5), HIJETS::algo_prefix + "_ParticleFlow_r05"); + pflowjetreco->set_algo_node(HIJETS::jet_node); + pflowjetreco->set_input_node("ELEMENT"); + pflowjetreco->Verbosity(verbosity); + se->registerSubsystem(pflowjetreco); + + // exit back to HIJetReco() + return; +} + +// ---------------------------------------------------------------------------- +//! Run background-subtracted jet reconstruction +// ---------------------------------------------------------------------------- +void HIJetReco() +{ + // if simulation, make appropriate truth jets + if (Enable::HIJETS_MC && Enable::HIJETS_TRUTH) MakeHITruthJets(); + + // run approriate jet reconstruction routines + if (Enable::HIJETS_TOWER) MakeHITowerJets(); + if (Enable::HIJETS_TRACK) MakeHITrackJets(); + if (Enable::HIJETS_PFLOW) MakeHIPFlowJets(); +} + +// ---------------------------------------------------------------------------- +//! Determine rho from tower input to jet reco (necessary for jet QA) +// ---------------------------------------------------------------------------- +void DoRhoCalculation() +{ + // set verbosity + // int verbosity = std::max(Enable::VERBOSITY, Enable::HIJETS_VERBOSITY); + + //--------------- + // Fun4All server + //--------------- + Fun4AllServer *se = Fun4AllServer::instance(); + + // run rho calculations w/ towers + if (Enable::HIJETS_TOWER) + { + DetermineTowerRho *towRhoCalc = new DetermineTowerRho(); + towRhoCalc->add_method(TowerRho::Method::AREA); + towRhoCalc->add_method(TowerRho::Method::MULT); + TowerJetInput *rho_incemc = new TowerJetInput(Jet::CEMC_TOWERINFO_RETOWER, HIJETS::tower_prefix); + TowerJetInput *rho_inihcal = new TowerJetInput(Jet::HCALIN_TOWERINFO, HIJETS::tower_prefix); + TowerJetInput *rho_inohcal = new TowerJetInput(Jet::HCALOUT_TOWERINFO, HIJETS::tower_prefix); + if (HIJETS::do_vertex_type) + { + rho_incemc->set_GlobalVertexType(HIJETS::vertex_type); + rho_inihcal->set_GlobalVertexType(HIJETS::vertex_type); + rho_inohcal->set_GlobalVertexType(HIJETS::vertex_type); + } + towRhoCalc->add_tower_input(rho_incemc); + towRhoCalc->add_tower_input(rho_inihcal); + towRhoCalc->add_tower_input(rho_inohcal); + se->registerSubsystem(towRhoCalc); + } + + // run rho calculations w/ tracks + if (Enable::HIJETS_TRACK) + { + DetermineEventRho *trkRhoCalc = new DetermineEventRho(); + trkRhoCalc->add_method(EventRho::Method::AREA, "EventRho_AREA"); + trkRhoCalc->add_method(EventRho::Method::MULT, "EventRho_MULT"); + trkRhoCalc->add_input(new TrackJetInput(Jet::SRC::TRACK)); + se->registerSubsystem(trkRhoCalc); + } + + // exit back to main macro + return; +} + +#endif diff --git a/JetProduction/Jet_QA.C b/JetProduction/Jet_QA.C new file mode 100644 index 000000000..2526174b1 --- /dev/null +++ b/JetProduction/Jet_QA.C @@ -0,0 +1,806 @@ +#ifndef MACRO_JETQA_C +#define MACRO_JETQA_C + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +R__LOAD_LIBRARY(libjetqa.so) + +// ---------------------------------------------------------------------------- +//! General options for jet QA +// ---------------------------------------------------------------------------- +namespace Enable +{ + int JETQA_VERBOSITY = 0; ///< verbosity of jet qa +} + +// ---------------------------------------------------------------------------- +//! Namespace to collect various enums, default arguments, etc. +// ---------------------------------------------------------------------------- +namespace JetQA +{ + + // flags -------------------------------------------------------------------- + + ///! Set to true if input jets utilize tracks (e.g. via particle flow) + bool HasTracks = false; + + ///! Set to true if input jets utilize calorimeters(e.g. via particle flow) + bool HasCalos = false; + + ///! Set to true if jets have background subtracted + ///! (e.g. HIJetReco.C output) + bool UseBkgdSub = false; + + ///! Set to true to generate histograms for no trigger selection + bool DoInclusive = true; + + ///! Set to true to generate histograms for a specified set of triggers + bool DoTriggered = true; + + ///! Set to true to use pp-specific options (otherwise assumes AuAU) + bool DoPP = false; + + ///! Set to true to use oo-specific options (otherwise assumes AuAU) + bool DoOO = false; + + ///! Set to true to restrict minimum jet pt to trigger threshold + bool RestrictPtToTrig = false; + + ///! Set to true to restrict jet eta acceptance by resolution parameter + bool RestrictEtaByR = true; + + // enums -------------------------------------------------------------------- + + ///! enumerates possible jet types + enum Type + { + AntiKtTowerSubR02, + AntiKtTowerSubR03, + AntiKtTowerSubR04, + AntiKtTowerSubR05, + AntiKtTowerR02, + AntiKtTowerR03, + AntiKtTowerR04, + AntiKtTowerR05, + AntiKtTrackR02, + AntiKtTrackR03, + AntiKtTrackR04, + AntiKtTrackR05 + }; + + ///! enumerates possible jet sources + enum class Source : uint32_t + { + Calos = 0, + Tracks = 1 + }; + + // constants ---------------------------------------------------------------- + + ///! Min jet pt in GeV/c + double MinJetPt = 6.; + + ///! Max jet pt in GeV/c + double MaxJetPt = 100.; + + ///! Min eta acceptance + double MinAcceptEta = -1.1; + + ///! Max eta acceptance + double MaxAcceptEta = 1.1; + + // maps --------------------------------------------------------------------- + + ///! Map from trigger to histogram tag + std::map GL1Tag = { + {JetQADefs::GL1::Clock, "clock"}, + {JetQADefs::GL1::ZDCS, "zdcs"}, + {JetQADefs::GL1::ZDCN, "zdcn"}, + {JetQADefs::GL1::ZDCNS, "zdcns"}, + {JetQADefs::GL1::HCalSingle, "hcalsingle"}, + {JetQADefs::GL1::MBDS, "mbds"}, + {JetQADefs::GL1::MBDN, "mbdn"}, + {JetQADefs::GL1::MBDNS1, "mbdns1"}, + {JetQADefs::GL1::MBDNS2, "mbdns2"}, + {JetQADefs::GL1::MBDNSVtx10, "mbdnsvtx10"}, + {JetQADefs::GL1::MBDNSVtx30, "mbdnsvtx30"}, + {JetQADefs::GL1::MBDNSVtx60, "mbdnsvtx60"}, + {JetQADefs::GL1::MBDNSHCalSingle, "mbdnshcalsingle"}, + {JetQADefs::GL1::MBDNSJet1, "mbdnsjet1"}, + {JetQADefs::GL1::MBDNSJet2, "mbdnsjet2"}, + {JetQADefs::GL1::MBDNSJet3, "mbdnsjet3"}, + {JetQADefs::GL1::MBDNSJet4, "mbdnsjet4"}, + {JetQADefs::GL1::Jet1, "jet1"}, + {JetQADefs::GL1::Jet2, "jet2"}, + {JetQADefs::GL1::Jet3, "jet3"}, + {JetQADefs::GL1::Jet4, "jet4"}, + {JetQADefs::GL1::MBDNSPhoton1, "mbdnsphoton1"}, + {JetQADefs::GL1::MBDNSPhoton2, "mbdnsphoton2"}, + {JetQADefs::GL1::MBDNSPhoton3, "mbdnsphoton3"}, + {JetQADefs::GL1::MBDNSPhoton4, "mbdnsphoton4"}, + {JetQADefs::GL1::Photon1, "photon1"}, + {JetQADefs::GL1::Photon2, "photon2"}, + {JetQADefs::GL1::Photon3, "photon3"}, + {JetQADefs::GL1::Photon4, "photon4"}, + {JetQADefs::GL1::MBDNS2Vtx10, "mbdns2vtx10"}, + {JetQADefs::GL1::MBDNS2Vtx30, "mbdns2vtx30"}, + {JetQADefs::GL1::MBDNS2Vtx60, "mbdns2vtx60"}, + {JetQADefs::GL1::MBDNS2Vtx150, "mbdns2vtx150"}, + {JetQADefs::GL1::MBDNS2Photon6Vtx10, "mbdns2photon6vtx10"}, + {JetQADefs::GL1::MBDNS2Photon8Vtx10, "mbdns2photon8vtx10"}, + {JetQADefs::GL1::MBDNS2Photon10Vtx10, "mbdns2photon10vtx10"}, + {JetQADefs::GL1::MBDNS2Photon12Vtx10, "mbdns2photon12vtx10"}, + {JetQADefs::GL1::MBDNS2Photon6Vtx150, "mbdns2photon6vtx150"}, + {JetQADefs::GL1::MBDNS2Photon8Vtx150, "mbdns2photon8vtx150"}, + {JetQADefs::GL1::MBDNS2Photon10Vtx150, "mbdns2photon10vtx150"}, + {JetQADefs::GL1::MBDNS2Photon12Vtx150, "mbdns2photon12vtx150"}, + {JetQADefs::GL1::Inclusive, "inclusive"}}; + + ///! Map from jet type to input node + std::map JetInput = { + {Type::AntiKtTowerSubR02, "AntiKt_Tower_r02_Sub1"}, + {Type::AntiKtTowerSubR03, "AntiKt_Tower_r03_Sub1"}, + {Type::AntiKtTowerSubR04, "AntiKt_Tower_r04_Sub1"}, + {Type::AntiKtTowerSubR05, "AntiKt_Tower_r05_Sub1"}, + {Type::AntiKtTowerR02, "AntiKt_Tower_r02"}, + {Type::AntiKtTowerR03, "AntiKt_Tower_r03"}, + {Type::AntiKtTowerR04, "AntiKt_Tower_r04"}, + {Type::AntiKtTowerR05, "AntiKt_Tower_r05"}, + {Type::AntiKtTrackR02, "AntiKt_Track_r02"}, + {Type::AntiKtTrackR03, "AntiKt_Track_r03"}, + {Type::AntiKtTrackR04, "AntiKt_Track_r04"}, + {Type::AntiKtTrackR05, "AntiKt_Track_r05"}}; + + ///! Map from jet type to histogram tag + std::map JetTag = { + {Type::AntiKtTowerSubR02, "towersub1_antikt_r02"}, + {Type::AntiKtTowerSubR03, "towersub1_antikt_r03"}, + {Type::AntiKtTowerSubR04, "towersub1_antikt_r04"}, + {Type::AntiKtTowerSubR05, "towersub1_antikt_r05"}, + {Type::AntiKtTowerR02, "tower_antikt_r02"}, + {Type::AntiKtTowerR03, "tower_antikt_r03"}, + {Type::AntiKtTowerR04, "tower_antikt_r04"}, + {Type::AntiKtTowerR05, "tower_antikt_r05"}, + {Type::AntiKtTrackR02, "track_antikt_r02"}, + {Type::AntiKtTrackR03, "track_antikt_r03"}, + {Type::AntiKtTrackR04, "track_antikt_r04"}, + {Type::AntiKtTrackR05, "track_antikt_r05"}}; + + ///! Map from jet type to resolution parameter + std::map JetRes = { + {Type::AntiKtTowerSubR02, 0.2}, + {Type::AntiKtTowerSubR03, 0.3}, + {Type::AntiKtTowerSubR04, 0.4}, + {Type::AntiKtTowerSubR05, 0.5}, + {Type::AntiKtTowerR02, 0.2}, + {Type::AntiKtTowerR03, 0.3}, + {Type::AntiKtTowerR04, 0.4}, + {Type::AntiKtTowerR05, 0.5}, + {Type::AntiKtTrackR02, 0.2}, + {Type::AntiKtTrackR03, 0.3}, + {Type::AntiKtTrackR04, 0.4}, + {Type::AntiKtTrackR05, 0.5}}; + + // methods ------------------------------------------------------------------ + + // -------------------------------------------------------------------------- + //! Get trigger tag + // -------------------------------------------------------------------------- + /*! If a trigger index is provided (i.e. a trigger is being selected), + * function returns the correpsonding tag. Otherwise (i.e. NO trigger + * is being selected), returns the inclusive tag. + */ + inline std::string GetTriggerTag(std::optional trg = std::nullopt) + { + std::string tag; + if (trg.has_value()) + { + tag.append("_" + GL1Tag[trg.value()]); + } + else + { + tag.append("_" + GL1Tag[JetQADefs::GL1::Inclusive]); + } + return tag; + + } // end 'GetTriggerTag(std::optional)' + + // -------------------------------------------------------------------------- + //! Get minimum jet pt based on which trigger fired + // -------------------------------------------------------------------------- + // FIXME it might be prudent to allow for thresholds to change vs. run + // number... Before run 46038, the Jet1, 2 thresholds were different + inline double GetMinJetPt(const uint32_t trg = JetQADefs::GL1::MBDNSJet1) + { + // by defult, set min to global constant + double ptJetMin = MinJetPt; + + // if restricting pt to trigger threhsold, pick out relevant threshold + if (RestrictPtToTrig) + { + switch (trg) + { + // Jet4 threshold + case JetQADefs::GL1::MBDNSJet4: + [[fallthrough]]; + case JetQADefs::GL1::Jet4: + ptJetMin = 11.; + break; + + // Jet3 threshold + case JetQADefs::GL1::MBDNSJet3: + [[fallthrough]]; + case JetQADefs::GL1::Jet3: + ptJetMin = 10.; + break; + + // Jet2 threshold + case JetQADefs::GL1::MBDNSJet2: + [[fallthrough]]; + case JetQADefs::GL1::Jet2: + ptJetMin = 9.; + break; + + // Jet1 threshold (default value) + case JetQADefs::GL1::MBDNSJet1: + [[fallthrough]]; + case JetQADefs::GL1::Jet1: + [[fallthrough]]; + default: + ptJetMin = 6.; + break; + } + } + return ptJetMin; + + } // end 'GetMinJetPt(uint32_t)' + + // -------------------------------------------------------------------------- + //! Get default jet pt range + // -------------------------------------------------------------------------- + inline std::pair GetJetPtRange(std::optional trg = std::nullopt) + { + std::pair ptJetRange; + if (trg.has_value()) + { + ptJetRange = {GetMinJetPt(trg.value()), MaxJetPt}; + } + else + { + ptJetRange = {GetMinJetPt(), MaxJetPt}; + } + return ptJetRange; + + } // end 'GetJetPtRange(std::optional)' + + // -------------------------------------------------------------------------- + //! Get default jet eta range + // -------------------------------------------------------------------------- + inline std::pair GetJetEtaRange(const double res = 0.) + { + // determine relevant min/max + const double etaMin = RestrictEtaByR ? MinAcceptEta + res : MinAcceptEta; + const double etaMax = RestrictEtaByR ? MaxAcceptEta - res : MaxAcceptEta; + + // return range + std::pair etaJetRange = {etaMin, etaMax}; + return etaJetRange; + + } // end 'GetJetEtaRange(double)' + + // -------------------------------------------------------------------------- + //! Get default list of triggers to analyze + // -------------------------------------------------------------------------- + inline std::vector GetDefaultTriggerList() + { + // default jets for pp + std::vector vecDefaultTrigsPP = { + JetQADefs::GL1::MBDNS1, + JetQADefs::GL1::MBDNSJet1, + JetQADefs::GL1::MBDNSJet2, + JetQADefs::GL1::MBDNSJet3, + JetQADefs::GL1::MBDNSJet4, + JetQADefs::GL1::Jet1, + JetQADefs::GL1::Jet2, + JetQADefs::GL1::Jet3, + JetQADefs::GL1::Jet4}; + + // default triggers for 2025 AuAu + std::vector vecDefaultTrigsAA = { + JetQADefs::GL1::MBDNS2Vtx10, + JetQADefs::GL1::MBDNS2Vtx150}; + + std::vector vecDefaultTrigsOO = { + JetQADefs::GL1::MBDNS1, + JetQADefs::GL1::MBDNSJet1, + JetQADefs::GL1::MBDNSJet2, + JetQADefs::GL1::MBDNSJet3, + JetQADefs::GL1::MBDNSJet4, + JetQADefs::GL1::Jet1, + JetQADefs::GL1::Jet2, + JetQADefs::GL1::Jet3, + JetQADefs::GL1::Jet4}; + + if (JetQA::DoPP) + { + return vecDefaultTrigsPP; + } + if (JetQA::DoOO) + { + return vecDefaultTrigsOO; + } + + return vecDefaultTrigsAA; + + } // end 'GetDefaultTriggerList()' + // -------------------------------------------------------------------------- + //! Get list of jets to analyze + // -------------------------------------------------------------------------- + inline std::vector GetJetsToQA(Source source) + { + std::vector vecJetsToQA; + switch (source) + { + case Source::Calos: + if (UseBkgdSub) + { + vecJetsToQA.insert(vecJetsToQA.end(), + {Type::AntiKtTowerSubR02, + Type::AntiKtTowerSubR03, + Type::AntiKtTowerSubR04, + Type::AntiKtTowerSubR05}); + } + else + { + vecJetsToQA.insert(vecJetsToQA.end(), + {Type::AntiKtTowerR02, + Type::AntiKtTowerR03, + Type::AntiKtTowerR04, + Type::AntiKtTowerR05}); + } + break; + + case Source::Tracks: + vecJetsToQA.insert(vecJetsToQA.end(), + {Type::AntiKtTrackR02, + Type::AntiKtTrackR03, + Type::AntiKtTrackR04, + Type::AntiKtTrackR05}); + break; + } + return vecJetsToQA; + + } // end 'GetJetsToQA(Source)' + +} // namespace JetQA + +// ---------------------------------------------------------------------------- +//! Create QA modules for all jets, regardless of constituents +// ---------------------------------------------------------------------------- +void CommonJetQA(std::optional trg = std::nullopt) +{ + // set verbosity + int verbosity = std::max(Enable::JETQA_VERBOSITY, Enable::QA_VERBOSITY); + if (verbosity > 1) + { + std::cout << ">>> Entering CommonJetQA() with trg=" + << (trg.has_value() ? std::to_string(trg.value()) : "none") + << std::endl; + } + + // grab appropriate trigger tag + std::string trig_tag = JetQA::GetTriggerTag(trg); + + // grab default pt, eta ranges + std::pair ptJetRange = JetQA::GetJetPtRange(trg); + std::pair etaJetMaxRange = JetQA::GetJetEtaRange(); + + // connect to f4a server + Fun4AllServer* se = Fun4AllServer::instance(); + + // create modules that are independent of/take multiple R values ------------ + + // initialize and register mass, eta, and pt qa module + // for calos if need be + if (JetQA::HasCalos) + { + auto vecCalos = JetQA::GetJetsToQA(JetQA::Source::Calos); + if (!vecCalos.empty()) + { + JetKinematicCheck* kinematicQA = new JetKinematicCheck( + "JetKinematicCheck" + trig_tag + (JetQA::UseBkgdSub ? "_towersub1_antikt" : "_tower_antikt"), + JetQA::JetInput[vecCalos[0]], + JetQA::JetInput[vecCalos[1]], + JetQA::JetInput[vecCalos[2]], + JetQA::JetInput[vecCalos[3]]); + kinematicQA->Verbosity(verbosity); + kinematicQA->setHistTag(""); + kinematicQA->setRestrictEtaRange(JetQA::RestrictEtaByR); + kinematicQA->setPtRange(ptJetRange.first, ptJetRange.second); + kinematicQA->setEtaRange(etaJetMaxRange.first, etaJetMaxRange.second); + if (trg.has_value()) + { + kinematicQA->setTrgToSelect(trg.value()); + } + se->registerSubsystem(kinematicQA); + } + else + { + std::cerr << ">>> Warning: trying to register JetKinematicCheck for tower jets, but jet list is empty!" + << std::endl; + } + } + + // initialize and register mass, eta, and pt qa module + // for tracks if need be + if (JetQA::HasTracks) + { + auto vecTrk = JetQA::GetJetsToQA(JetQA::Source::Tracks); + if (!vecTrk.empty()) + { + JetKinematicCheck* kinematicQA_Tracks = new JetKinematicCheck( + "JetKinematicCheck" + trig_tag + "_track_antikt", + JetQA::JetInput[vecTrk[0]], + JetQA::JetInput[vecTrk[1]], + JetQA::JetInput[vecTrk[2]], + JetQA::JetInput[vecTrk[3]]); + kinematicQA_Tracks->Verbosity(verbosity); + kinematicQA_Tracks->setHistTag(""); + kinematicQA_Tracks->setRestrictEtaRange(JetQA::RestrictEtaByR); + kinematicQA_Tracks->setPtRange(ptJetRange.first, ptJetRange.second); + kinematicQA_Tracks->setEtaRange(etaJetMaxRange.first, etaJetMaxRange.second); + if (trg.has_value()) + { + kinematicQA_Tracks->setTrgToSelect(trg.value()); + } + se->registerSubsystem(kinematicQA_Tracks); + } + else + { + std::cerr << ">>> Warning: trying to register JetKinematicCheck for track jets, but jet list is empty!" + << std::endl; + } + } + + // create modules that take single R values --------------------------------- + + // get relevant jets to QA + std::vector vecJetsToQA; + if (JetQA::HasCalos) + { + auto vecCalJets = JetQA::GetJetsToQA(JetQA::Source::Calos); + vecJetsToQA.insert(vecJetsToQA.end(), vecCalJets.begin(), vecCalJets.end()); + } + if (JetQA::HasTracks) + { + auto vecTrkJets = JetQA::GetJetsToQA(JetQA::Source::Tracks); + vecJetsToQA.insert(vecJetsToQA.end(), vecTrkJets.begin(), vecTrkJets.end()); + } + + // loop over resolution parameters + for (uint32_t jet : vecJetsToQA) + { + // get R-dependent eta range + std::pair etaJetRange = JetQA::GetJetEtaRange(JetQA::JetRes[jet]); + + // intialize and register dijetQA + DijetQA* dijetQA = new DijetQA( + "DijetQA" + trig_tag + "_" + JetQA::JetTag[jet], + JetQA::JetInput[jet]); + dijetQA->Verbosity(verbosity); + dijetQA->setPtLeadRange(ptJetRange.first, ptJetRange.second); + dijetQA->setPtSubRange(1.0, ptJetRange.second); + dijetQA->setEtaRange(etaJetRange.first, etaJetRange.second); + if (trg.has_value()) + { + dijetQA->setTrgToSelect(trg.value()); + } + se->registerSubsystem(dijetQA); + + } // end jet loop + return; + +} // end 'CommonJetQA(std::optional)' + +// ---------------------------------------------------------------------------- +//! Create QA modules for jets with tracks +// ---------------------------------------------------------------------------- +void JetsWithTracksQA(std::optional trg = std::nullopt) +{ + // set verbosity + int verbosity = std::max(Enable::JETQA_VERBOSITY, Enable::QA_VERBOSITY); + if (verbosity > 1) + { + std::cout << ">>> Entering JetsWithTracksQA() with trg=" + << (trg.has_value() ? std::to_string(trg.value()) : "none") + << std::endl; + } + + // grab appropriate trigger tag + std::string trig_tag = JetQA::GetTriggerTag(trg); + + // grab default pt, eta ranges + std::pair ptJetRange = JetQA::GetJetPtRange(trg); + std::pair etaJetRange = JetQA::GetJetEtaRange(); + + // get list of jet nodes to analyze + std::vector vecJetsToQA = JetQA::GetJetsToQA(JetQA::Source::Tracks); + + // connect to f4a server + Fun4AllServer* se = Fun4AllServer::instance(); + + // initialize and register event-wise rho check + RhosinEvent* evtRhoQA = new RhosinEvent("EventWiseTrackRho" + trig_tag, ""); + evtRhoQA->Verbosity(verbosity); + evtRhoQA->add_area_rho_node("EventRho_AREA"); + evtRhoQA->add_mult_rho_node("EventRho_MULT"); + if (trg.has_value()) + { + evtRhoQA->setTrgToSelect(trg.value()); + } + se->registerSubsystem(evtRhoQA); + + // create modules that take single R values --------------------------------- + + // loop over resolution parameters + for (uint32_t jet : vecJetsToQA) + { + // intialize and register sum track vs. jet pt qa module + StructureinJets* jetStructQA = new StructureinJets( + "StructureInJets" + trig_tag + "_" + JetQA::JetTag[jet], + JetQA::JetInput[jet], + "SvtxTrackMap", + ""); + jetStructQA->Verbosity(verbosity); + jetStructQA->writeToOutputFile(false); + jetStructQA->setTrkPtCut(0.1); + jetStructQA->setTrkQualityCut(100.0); // set to be big for now, should refine later + jetStructQA->setTrkNSilCut(4); + jetStructQA->setTrkNTPCCut(25); + jetStructQA->setJetRadius(JetQA::JetRes[jet]); + jetStructQA->setJetPtRange(ptJetRange.first, ptJetRange.second); + jetStructQA->setJetEtaRange(etaJetRange.first, etaJetRange.second); + if (trg.has_value()) + { + jetStructQA->setTrgToSelect(trg.value()); + } + se->registerSubsystem(jetStructQA); + + // initialize and register track jet qa module + TrksInJetQA* trksInJetQA = new TrksInJetQA("TrksInJet" + trig_tag + "_" + JetQA::JetTag[jet]); + trksInJetQA->SetHistSuffix(""); + trksInJetQA->Configure( + {.outMode = TrksInJetQA::OutMode::QA, + .verbose = verbosity, + .doDebug = false, + .doInclusive = false, // n.b. inclusive case covered by track QA + .doInJet = true, + .doHitQA = false, + .doClustQA = true, + .doTrackQA = true, + .doJetQA = true, + .doSubsysHist = false, + .doOptHist = false, // turn off optional histograms + .rJet = JetQA::JetRes[jet], + .jetInNode = JetQA::JetInput[jet]}); + if (trg.has_value()) + { + trksInJetQA->SetTrgToSelect(trg.value()); + } + se->registerSubsystem(trksInJetQA); + + } // end jet loop + return; + +} // end 'JetWithTracksQA(std::optional)' + +// ---------------------------------------------------------------------------- +//! Create QA modules for jets with clusters (Calo) +// ---------------------------------------------------------------------------- +void JetsWithCaloQA(std::optional trg = std::nullopt) +{ + // set verbosity + int verbosity = std::max(Enable::JETQA_VERBOSITY, Enable::QA_VERBOSITY); + if (verbosity > 1) + { + std::cout << ">>> Entering JetsWithCaloQA() with trg=" + << (trg.has_value() ? std::to_string(trg.value()) : "none") + << std::endl; + } + + // grab appropriate trigger tag + std::string trig_tag = JetQA::GetTriggerTag(trg); + + // grab default pt, eta ranges + std::pair ptJetRange = JetQA::GetJetPtRange(trg); + // std::pair etaJetRange = JetQA::GetJetEtaRange(); + + // connect to f4a server + Fun4AllServer* se = Fun4AllServer::instance(); + + // get list of jet nodes to analyze + std::vector vecJetsToQA = JetQA::GetJetsToQA(JetQA::Source::Calos); + + // initialize and register emcluster jet kinematic QA + EMClusterKinematics* emclusterJetsQA = new EMClusterKinematics( + "EMClusterKinematics" + trig_tag, + "CLUSTERINFO_CEMC", + ""); + emclusterJetsQA->SetDoOptHist(true); + if (trg.has_value()) + { + emclusterJetsQA->SetTrgToSelect(trg.value()); + } + emclusterJetsQA->Verbosity(0); + se->registerSubsystem(emclusterJetsQA); + + // initialize and register event-wise rho check + RhosinEvent* evtRhoQA = new RhosinEvent("EventWiseCaloRho" + trig_tag, ""); + evtRhoQA->Verbosity(verbosity); + evtRhoQA->add_area_rho_node("TowerRho_AREA"); + evtRhoQA->add_mult_rho_node("TowerRho_MULT"); + if (trg.has_value()) + { + evtRhoQA->setTrgToSelect(trg.value()); + } + se->registerSubsystem(evtRhoQA); + + // initialize and register jet seed counter qa module + if (!JetQA::DoPP) + { + JetSeedCount* jetSeedQA = new JetSeedCount( + "JetSeedCount" + trig_tag + "_towersub1_antikt_r02", + "AntiKt_Tower_r04", // n.b. unused in module + "AntiKt_TowerInfo_HIRecoSeedsRaw_r02", + "AntiKt_TowerInfo_HIRecoSeedsSub_r02"); + jetSeedQA->Verbosity(verbosity); + jetSeedQA->setHistTag(""); + jetSeedQA->setPtRange(ptJetRange.first, ptJetRange.second); + jetSeedQA->setEtaRange(JetQA::MinAcceptEta, JetQA::MaxAcceptEta); + jetSeedQA->setWriteToOutputFile(false); + jetSeedQA->setPPMode(JetQA::DoPP); + if (trg.has_value()) + { + jetSeedQA->setTrgToSelect(trg.value()); + } + se->registerSubsystem(jetSeedQA); + } + + // register underlying event module + if (!JetQA::DoPP && JetQA::DoOO) + { + UEvsEtaCentrality* UEvsEtaCentralityQA = new UEvsEtaCentrality("UEvsEtaCent" + trig_tag); + UEvsEtaCentralityQA->SetConfig( + { + .debug = false, + .histTag = "", + }); + if (trg.has_value()) + { + UEvsEtaCentralityQA->SetConfig( + { + .doTrgSelect = true, + .trgToSelect = trg.value(), + }); + } + se->registerSubsystem(UEvsEtaCentralityQA); + } + // create modules that take single R values --------------------------------- + + // loop over resolution parameters + for (uint32_t jet : vecJetsToQA) + { + // get R-dependent eta range + std::pair etaJetRange = JetQA::GetJetEtaRange(JetQA::JetRes[jet]); + + // initialize and register constituent checks + ConstituentsinJets* jetCstQA = new ConstituentsinJets( + "ConstituentsInJets" + trig_tag + "_" + JetQA::JetTag[jet], + JetQA::JetInput[jet], + "TowerInfoBackground_Sub1", + "", + "TOWERINFO_CALIB_CEMC_RETOWER", + "TOWERINFO_CALIB_HCALIN", + "TOWERINFO_CALIB_HCALOUT"); + jetCstQA->Verbosity(verbosity); + jetCstQA->setPtRange(ptJetRange.first, ptJetRange.second); + jetCstQA->setEtaRange(etaJetRange.first, etaJetRange.second); + jetCstQA->setPPMode(JetQA::DoPP); + if (trg.has_value()) + { + jetCstQA->setTrgToSelect(trg.value()); + } + se->registerSubsystem(jetCstQA); + + } // end jet loop + return; + +} // end 'JetsWithCaloQA(std::optional)' + +// ---------------------------------------------------------------------------- +//! Run jet QA +// ---------------------------------------------------------------------------- +void Jet_QA(const std::vector& vecTrigsToUse = JetQA::GetDefaultTriggerList()) +{ + // if running with calos, instantiate & register the + // CaloStatusMapper QA module + if (JetQA::HasCalos) + { + // set verbosity + int verbosity = std::max(Enable::JETQA_VERBOSITY, Enable::QA_VERBOSITY); + + // connect to f4a server + Fun4AllServer* se = Fun4AllServer::instance(); + + // register module + CaloStatusMapper* caloStatusQA = new CaloStatusMapper("CaloStatusMapper"); + caloStatusQA->SetConfig( + { + .debug = false, + .doNorm = false, // do NOT try to normalize histograms + .doOptHist = false, // turn off extra histograms + .histTag = "", + .doTrgSelect = false // n.b. differential in trigger not useful here + }); + + se->Verbosity(verbosity); + se->registerSubsystem(caloStatusQA); + } + + // run in inclusive mode if needed + if (JetQA::DoInclusive) + { + CommonJetQA(); + if (JetQA::HasTracks) + { + JetsWithTracksQA(); + } + if (JetQA::HasCalos) + { + JetsWithCaloQA(); + } + } + + // run in triggered mode if needed + if (JetQA::DoTriggered) + { + for (uint32_t trg : vecTrigsToUse) + { + CommonJetQA(trg); + if (JetQA::HasTracks) + { + JetsWithTracksQA(trg); + } + if (JetQA::HasCalos) + { + JetsWithCaloQA(trg); + } + } // end trigger loop + } + return; + +} // end 'Jet_QA(std::vector)' + +#endif + +// end ------------------------------------------------------------------------ diff --git a/common/Jet_QA.C b/common/Jet_QA.C index 42dc6971d..84bdd46c2 100644 --- a/common/Jet_QA.C +++ b/common/Jet_QA.C @@ -625,6 +625,7 @@ void JetsWithCaloQA(std::optional trg = std::nullopt) // get list of jet nodes to analyze std::vector vecJetsToQA = JetQA::GetJetsToQA(JetQA::Source::Calos); + /* // initialize and register emcluster jet kinematic QA EMClusterKinematics* emclusterJetsQA = new EMClusterKinematics( "EMClusterKinematics" + trig_tag, @@ -637,7 +638,21 @@ void JetsWithCaloQA(std::optional trg = std::nullopt) } emclusterJetsQA->Verbosity(verbosity); se->registerSubsystem(emclusterJetsQA); - + */ + EMClusterKinematics* em = new EMClusterKinematics( + "EMClusterKinematics" + trig_tag, + "CLUSTERINFO_CEMC", + ""); + + // make the output name explicit + em->SetHistTag(JetQA::GetTriggerTag(trg)); // or just JetQA::GL1Tag[trg.value()] when trg has value + + // for debugging: disable trigger selection entirely + if (trg.has_value()) em->SetTrgToSelect(trg.value()); + + em->SetDoOptHist(false); + se->registerSubsystem(em); + // initialize and register event-wise rho check RhosinEvent* evtRhoQA = new RhosinEvent("EventWiseCaloRho" + trig_tag, ""); evtRhoQA->Verbosity(verbosity);