epic-analysis/src/Analysis.h

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 Christopher Dilks, Connor Pecar, Duane Byer, Sanghwa Park, Matthew McEneaney, Brian Page
0003
0004 #pragma once
0005
0006 #include <stdlib.h>
0007 #include <stdio.h>
0008 #include <string>
0009 #include <iostream>
0010 #include <sstream>
0011 #include <fstream>
0012 #include <vector>
0013 #include <map>
0014 #include <set>
0015 #include <stdexcept>
0016 #include <functional>
0017 #include <fmt/format.h>
0018
0019 // root
0020 #include <TChain.h>
0021 #include <TObjArray.h>
0022 #include <TFile.h>
0023 #include <TRegexp.h>
0024
0025 // adage
0026 #include <adage/BinSet.h>
0027
0028 // epic-analysis
0029 #include "DataModel.h"
0030 #include "Histos.h"
0031 #include "HistosDAG.h"
0032 #include "Kinematics.h"
0033 #ifndef EXCLUDE_DELPHES
0034 #include "KinematicsJets.h"
0035 #endif
0036 #include "SidisTree.h"
0037 #include "HFSTree.h"
0038 #include "ParticleTree.h"
0039 #include "Weights.h"
0040 #include "CommonConstants.h"
0041
0042 class Analysis
0043 {
0044 public:
0045   Analysis(
0046        TString configFileName_="",
0047        TString outfilePrefix_=""
0048        );
0049   ~Analysis();
0050
0051   // number of bins for histograms
0052   const Int_t NBINS = 50;
0053   const Int_t NBINS_FULL = 10;
0054
0055   // bin schemes
0056     void AddBinScheme(TString varname); // add a new bin scheme
0057   BinSet *BinScheme(TString varname); // access bin scheme by name
0058   std::map<TString,BinSet*> GetBinSchemes() { return binSchemes; }; // get full set of bin schemes
0059
0060   // add a new final state bin
0061   void AddFinalState(TString finalStateN);
0062
0063   // common settings
0064   Bool_t verbose; // if true, print a lot more information
0065   Bool_t writeSidisTree;   // if true, write SidisTree (not binned)
0066   Bool_t writeHFSTree;      // if true, write HFSTree (not binned)
0067   Bool_t writeParticleTree; // if true, write ParticleTree (not binned)
0068   Long64_t maxEvents; /* default=0, which runs all events;
0069                * if > 0, run a maximum number of `maxEvents` events (useful for quick tests)
0070                          */
0071   Bool_t useBreitJets; // if true, use Breit jets, if using finalState `jets` (requires centauro)
0072   // set kinematics reconstruction method; see constructor for available methods
0073   void SetReconMethod(TString reconMethod_) { reconMethod=reconMethod_; };
0074   // choose which output sets to include
0075   std::map<TString,Bool_t> includeOutputSet;
0076   // maximum number of errors to print
0077   Long64_t errorCntMax;
0078
0079   // Jet Definition Quantities
0080   int jetAlg;
0081   double jetRad, jetMin;
0082   double jetMatchDR;
0083
0084   // add a group of files to the analysis, where all of these files have a
0085   // common cross section `xs`, and Q2 range `Q2min` to `Q2max`
0086   void AddFileGroup(
0087         std::vector<std::string> fileNames,
0088         Long64_t totalEntries,
0089         Double_t xs,
0090         Double_t Q2min,
0091         Double_t Q2max,
0092         Double_t manualWeight
0093             );
0094
0095   // access HistosDAG
0096   std::shared_ptr<HistosDAG> GetHistosDAG() { return HD; };
0097
0098
0099   Double_t GetEventQ2Weight(Double_t Q2, Int_t guess=0);
0100   // after adding all files, estimate the weights for events in each Q2 range
0101     void CalculateEventQ2Weights();
0102
0103     Int_t GetEventQ2Idx(Double_t Q2, Int_t guess=0);
0104
0105     // run the analysis
0106     virtual void Execute() = 0;
0107
0108   protected:
0109
0110     // prepare to perform the analysis; in derived classes, define a method `Execute()`, which
0111     // will run the event loop; the first line of `Execute()` should call `Analysis::Prepare()`,
0112     // which set up common things like output files, `HistosDAG`, etc.
0113     void Prepare();
0114
0115     // finish the analysis; call `Analysis::Finish()` at the end of derived `Execute()` methods
0116     void Finish();
0117
0118    // print an error; if more than `errorCntMax` errors are printed, printing is suppressed
0119     void ErrorPrint(std::string message);
0120
0121     // `FillHistos(weight)` methods: fill histograms
0122     void FillHistosInclusive(Double_t wgt); // inclusive kinematics
0123     void FillHistos1h(Double_t wgt);        // single-hadron kinematics
0124     void FillHistosJets(Double_t wgt);      // jet kinematics
0125
0126     // shared objects
0127     std::unique_ptr<SidisTree>   ST;
0128     std::unique_ptr<HFSTree>      HFST;
0129     std::unique_ptr<ParticleTree> PT;
0130     std::shared_ptr<Kinematics>   kin, kinTrue;
0131 #ifndef EXCLUDE_DELPHES
0132     std::shared_ptr<KinematicsJets> kinJet, kinJetTrue;
0133 #endif
0134     std::shared_ptr<HistosDAG>    HD;
0135     std::unique_ptr<Weights> weightInclusive, weightTrack, weightJet;
0136     Double_t wInclusiveTotal, wTrackTotal, wJetTotal;
0137     Long64_t entriesTot;
0138     Long64_t errorCnt;
0139     const TString sep = "--------------------------------------------";
0140
0141     // setup / common settings
0142     std::vector<std::vector<std::string> > infiles;
0143     // A lookup index for guessing which Q2 range an event belongs to.
0144     std::vector<std::size_t> inLookup;
0145     std::vector<Double_t> Q2xsecs;
0146     std::vector<Double_t> Q2mins;
0147     std::vector<Double_t> Q2maxs;
0148     std::vector<Long64_t> Q2entries;
0149     std::vector<Double_t> Q2weights;
0150     std::vector<int> total_events = {0};
0151     TString configFileName,outfileName,outfilePrefix;
0152     TFile *outFile;
0153     Double_t eleBeamEn; // GeV
0154     Double_t ionBeamEn; // GeV
0155     Double_t crossingAngle; // mrad
0156     Double_t totalCrossSection;
0157     TString reconMethod;
0158
0159     // event loop objects
0160     Long64_t ENT;
0161     Double_t eleP,maxEleP;
0162     Double_t elePtrue, maxElePtrue;
0163     int pid;
0164     TString finalStateID;
0165
0166 #ifndef EXCLUDE_DELPHES
0167   fastjet::PseudoJet jet;
0168 #endif
0169
0170   // binning names / titles / etc.
0171   std::map<TString,TString> availableBinSchemes;
0172   std::map<TString,BinSet*> binSchemes;
0173   std::map<TString,TString> reconMethodToTitle;
0174   std::map<TString,TString> finalStateToTitle;
0175   std::map<int,TString> PIDtoFinalState;
0176   std::set<TString> activeFinalStates;
0177
0178   // check if Q2 `val` is between `min` and `max`; if `max==0`, only `val>=min` is checked
0179   template<class T> bool InQ2Range(T val, T min, T max, bool ignoreZero=false) {
0180     if (ignoreZero && !(val>0)) return true;
0181     if (max>0.0) return val>=min && val<=max;
0182     else         return val>=min;
0183   }
0184
0185   // container printing
0186   // mostly for debugging; if we need more than this, switch to using a common pretty printer library
0187   template<class O> void PrintStdVector(std::vector<O> vec, std::string name="") {
0188     if(name!="") fmt::print("{}: ",name);
0189     fmt::print("[ ");
0190     for(const auto elem : vec) fmt::print("{}, ",elem);
0191     fmt::print("]\n");
0192   }
0193   template<class O> void PrintStdVector2D(std::vector<std::vector<O>> vec, std::string name="") {
0194     if(name!="") fmt::print("{} = ",name);
0195     fmt::print("[\n");
0196     for(const auto elem : vec) PrintStdVector(elem,"  ");
0197     fmt::print("]\n");
0198   }
0199   template<class K, class V> void PrintStdMap(std::map<K,V> hash, std::string name="") {
0200     if(name!="") fmt::print("{} = ",name);
0201     fmt::print("{{\n");
0202     for(const auto [key,val] : hash) fmt::print("  {} => {},\n",key,val);
0203     fmt::print("}}\n");
0204   }
0205
0206 private:
0207
0208   // fill histograms, according to `fill_payload`
0209   void FillHistos(std::function<void(Histos*)> fill_payload);
0210
0211   ClassDef(Analysis,1);
0212 };