LargeRLCaloENC/src/LargeRLENC.h

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