LargeRLCaloENC/src/LargeRLENC_ThresholdScan.h

0001 // Tell emacs that this is a C++ source
0002 //  -*- C++ -*-.
0003 #ifndef LARGERLENC_H
0004 #define LARGERLENC_H
0005 #pragma GCC diagnostic push
0006 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0007 //fun4all
0008 #include <fun4all/SubsysReco.h>
0009 #include <fun4all/Fun4AllBase.h>
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011
0012 //phool
0013 #include <phool/PHCompositeNode.h>
0014 #include <phool/getClass.h>
0015
0016 //Calo towers
0017 #include <calobase/TowerInfoContainer.h>
0018 #include <calobase/TowerInfoContainerv1.h>
0019 #include <calobase/TowerInfoContainerv2.h>
0020 #include <calobase/TowerInfov2.h>
0021 #include <calobase/TowerInfov1.h>
0022 #include <calobase/TowerInfo.h>
0023 #include <calobase/RawTowerDefs.h>
0024 #include <calobase/RawTowerContainer.h>
0025 #include <calobase/RawTowerGeomContainer.h>
0026 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
0027
0028 //G4 objects
0029 #include <g4main/PHG4Particle.h>
0030 #include <g4main/PHG4Hit.h>
0031 #include <g4main/PHG4TruthInfoContainer.h>
0032
0033 #include <phhepmc/PHHepMCGenEvent.h>
0034 #include <phhepmc/PHHepMCGenEventMap.h>
0035 #include <HepMC/GenEvent.h>
0036
0037 //jetbase objects
0038 #include <jetbase/JetContainer.h>
0039 #include <jetbase/JetContainerv1.h>
0040 #include <jetbase/JetMapv1.h>
0041 #include <jetbase/JetMap.h>
0042 #include <jetbase/Jetv1.h>
0043 #include <jetbase/Jetv2.h>
0044 #include <jetbase/Jet.h>
0045 #include <jetbase/TowerJetInput.h>
0046 #include <jetbase/ClusterJetInput.h>
0047
0048 //fastjet
0049 #include <fastjet/PseudoJet.hh>
0050 #include <fastjet/ClusterSequence.hh>
0051
0052 //vertex stuff
0053 #include <globalvertex/GlobalVertex.h>
0054 #include <globalvertex/GlobalVertexMap.h>
0055
0056 //trigger
0057 #include <ffarawobjects/Gl1Packetv2.h>
0058 #include <ffarawobjects/Gl1Packetv1.h>
0059 #include <calotrigger/TriggerAnalyzer.h>
0060 #include <calotrigger/TriggerRunInfov1.h>
0061
0062 //c++
0063 #include <thread>
0064 #include <map>
0065 #include <utility>
0066 #include <vector>
0067 #include <string>
0068 #include <math.h>
0069 #include <array>
0070 #include <iostream>
0071 #include <fstream>
0072 #include <set>
0073 //root
0074 #include <TH1.h>
0075 #include <TH2.h>
0076 #include <TFile.h>
0077 #include <TTree.h>
0078 #include <TInterpreter.h>
0079 #include <TDirectory.h>
0080 //Homebrews
0081 #include <calorimetertowerenc/MethodHistograms.h>
0082 #include "DijetEventCuts.h"
0083 #include "HelperStructs.h"
0084
0085 #define PI 3.14159265358979323464
0086 class PHCompositeNode;
0087 class Jet;
0088 class JetContainer;
0089 class PHG4Hit;
0090 class PHG4Particle;
0091 class PHG4TruthInfoContainer;
0092 class TriggerAnalyzer;
0093
0094
0095 class LargeRLENC : public SubsysReco
0096 {
0097
0098  public:
0099     enum Regions{
0100         Full,
0101         Towards,
0102         Away,
0103         TransverseMax,
0104         TransverseMin
0105     };
0106     enum Calorimeter{
0107         All,
0108         EMCAL,
0109         IHCAL,
0110         OHCAL,
0111         TRUTH
0112     };
0113     LargeRLENC(const int n_run=0, const int n_segment=0, const float jet_min_pT=1.0, const bool data=false, const bool pedestal=false, std::fstream* ofs=nullptr, const std::string vari="E", const std::string &name = "LargeRLENC");
0114
0115     ~LargeRLENC() override {};
0116
0117   /** Called during initialization.
0118       Typically this is where you can book histograms, and e.g.
0119       register them to Fun4AllServer (so they can be output to file
0120       using Fun4AllServer::dumpHistos() method).
0121    */
0122     int Init(PHCompositeNode *topNode/*, bool* has_retower, bool *has_jets*/) override {
0123     // this first section is a backup way to get the status of the towers or jets
0124     /*  *has_retower=true;
0125         *has_jets=true;
0126         if(data){
0127             auto jets = findNode::getClass<JetContainer>(topNode, "AntiKt_Towers_r04");
0128             auto retower = findNode::getClass<TowerInfoContainer(topNode, "TOWERINFO_CALIB_CEMC_RETOWER");
0129             if(!jets) (*has_jets)=false;
0130             if(!retower) (*has_retower)=false;
0131         }*/
0132         return Fun4AllReturnCodes::EVENT_OK;
0133     }
0134
0135     /** Called for first event when run number is known.
0136         Typically this is where you may want to fetch data from
0137         database, because you know the run number. A place
0138         to book histograms which have to know the run number.
0139     */
0140     int InitRun(PHCompositeNode *topNode) override {return Fun4AllReturnCodes::EVENT_OK;}
0141
0142     /** Called for each event.
0143         This is where you do the real work.
0144     */
0145     int process_event(PHCompositeNode *topNode) override;
0146
0147     /// Clean up internals after each event.
0148     int ResetEvent(PHCompositeNode *topNode) override {return Fun4AllReturnCodes::EVENT_OK;}
0149
0150     /// Called at the end of each run.
0151     int EndRun(const int runnumber) override {return Fun4AllReturnCodes::EVENT_OK;}
0152
0153     /// Called at the end of all processing.
0154     int End(PHCompositeNode *topNode) override {
0155 //      if(write_to_file) this->text_out_file->close();
0156         return Fun4AllReturnCodes::EVENT_OK;
0157     };
0158
0159     /// Reset
0160     int Reset(PHCompositeNode * /*topNode*/) override {return Fun4AllReturnCodes::EVENT_OK;};
0161
0162     void Print(const std::string &what = "ALL") const override;
0163
0164     //and now for the unique stuff
0165     void addTower(int, TowerInfoContainer*, RawTowerGeomContainer_Cylinderv1*, std::map<std::array<float, 3>, float>*, RawTowerDefs::CalorimeterId);
0166
0167     std::array<float,3> HadronicEnergyBalence(Jet*, float, PHCompositeNode*);
0168     std::vector<std::array<float,3>> getJetEnergyRatios(JetContainerv1*, float, PHCompositeNode*);
0169     JetContainerv1* getJets(std::string, std::string, std::array<float, 3>, float ohcal_rat, PHCompositeNode*);
0170     void TruthRegion (std::map<std::array<float, 3>, float> , float,  std::string, std::array<float, 3>, float);
0171
0172     void CaloRegion(std::map<std::array<float, 3>, float>, std::map<std::array<float, 3>, float>, std::map<std::array<float, 3>, float>, float, std::string, std::array<float, 3>, float);
0173
0174     void SingleCaloENC( std::map<std::array<float, 3>, float>, float, std::array<float, 3>, bool, bool, std::map<int, std::pair<float, float>>, LargeRLENC::Calorimeter);
0175
0176     void CalculateENC(StrippedDownTowerWithThreshold*, std::vector<StrippedDownTowerWithThreshold>, bool, bool);
0177
0178     void JetEventObservablesBuilding(std::array<float, 3>, std::map<std::array<float, 3>, float>, std::map<float, float>*);
0179     float getR(float, float, float, float, bool print=false);
0180     bool triggerCut(bool, PHCompositeNode*);
0181     void Merger(TowerOutput*, std::vector<TowerOutput*>, std::set<float>, std::set<std::array<float, 3>>);
0182     DijetEventCuts* eventCut;
0183     void MakeEMCALRetowerMap(RawTowerGeomContainer_Cylinderv1* em_geom, TowerInfoContainer* emcal, RawTowerGeomContainer_Cylinderv1* h_geom, TowerInfoContainer* hcal );
0184     std::vector<std::array<float, 4>> Thresholds;
0185     std::map<int, std::pair<float, float>> emcal_lookup_table;
0186  private:
0187     std::string algo, radius, output_file_name;
0188     std::string ohcal_energy_towers="TOWERINFO_CALIB_HCALOUT", ihcal_energy_towers="TOWERINFO_CALIB_HCALIN", emcal_energy_towers="TOWERINFO_CALIB_CEMC";
0189     bool isRealData, pedestalData;
0190     int nRun, nSegment, m_Njets, n_evts, n_with_jets=0;
0191     float jetMinpT, MinpTComp;
0192     float ptoE=1.; //need to actually do some studies into this in order to get a meaningful conversion factor
0193     int m_region, m_calo;
0194     std::vector<std::pair<int, std::vector< std::pair<float, float> > > > t_e2c, t_e3c; //make this an indexable object
0195     float m_rl, m_rm, m_ri, m_e2c, m_e3c, m_Et, m_vtx, m_vty, m_vtz, m_philead, m_etalead;
0196     std::vector<std::pair< int, std::vector< std::pair< std::array<float, 3>, float> > > > t_e3c_full;
0197     std::map< std::string, std::array< std::map< std::pair< float, float >, float >, 3 > > t_pt_evt;
0198     float m_phi, m_eta;
0199     std::string which_variable; //Which varaible are we caluclating the EEC over (E, E_T, p, p_T)
0200     TTree* DijetQA, *EEC/*, *JetEvtObs*/;
0201     std::vector<std::vector<std::vector<MethodHistograms*>>> Region_vector, Tr_Region_vector;
0202     float m_etotal, m_eemcal, m_eihcal, m_eohcal;
0203     std::array<float, 3> m_vertex;
0204     std::vector<std::array<float, 3>> m_dijets;
0205     float m_xjl, m_Ajjl;
0206     TH1F* emcal_occup, *ihcal_occup, *ohcal_occup, *ohcal_rat_h;
0207     TH2F* ohcal_rat_occup, *ohcal_bad_hits, *bad_occ_em_oh, *bad_occ_em_h;
0208     TH1F* IE, *badIE, *goodIE;
0209     TH2F* E_IE, *badE_IE, *goodE_IE;
0210     std::vector<std::array<float, 3>> m_emcal, m_ihcal, m_ohcal; //3 points, eta, phi, et
0211     std::array<std::array<TH1F*, 3>, 3> Et_miss_hists;
0212     std::array<float, 5> thresh_mins;
0213     float ohcal_min=0.01; //7.5 MeV from jet 30 and 10 study
0214     float emcal_min=0.01; //50 MeV
0215     float ihcal_min=0.005; //7.5 MeV
0216     float all_min=0.065; //65 MeV
0217     //all these are conservative vals
0218     int n_steps=10;
0219 };
0220 #endif // LARGERLENC_H