CalorimeterJets/src/MyJetAnalysis.h

0001 #ifndef MYJETANALYSIS_H
0002 #define MYJETANALYSIS_H
0003
0004 #include <fun4all/SubsysReco.h>
0005
0006 #include <memory>
0007 #include <string>
0008 #include <vector>
0009 #include <TLorentzVector.h>
0010 #include <utility>  // std::pair, std::make_pair
0011
0012 #include <array>
0013 #include <g4jets/Jet.h>
0014
0015 #include <HepMC/GenEvent.h>
0016 #include <HepMC/GenParticle.h>
0017
0018 #include <fastjet/PseudoJet.hh>
0019
0020 class PHCompositeNode;
0021 class JetEvalStack;
0022 class CaloEvalStack;
0023 class TTree;
0024 class TH1;
0025 class FastJetAlgo;
0026
0027 /// \class MyJetAnalysis
0028 class MyJetAnalysis : public SubsysReco
0029 {
0030  public:
0031   MyJetAnalysis(
0032       const std::string &recojetname = "AntiKt_Tower_r02",
0033       const std::string &truthjetname = "AntiKt_Truth_r02",
0034       const std::string &outputfilename = "myjetanalysis.root");
0035
0036   virtual ~MyJetAnalysis();
0037
0038   //! set eta range
0039   void
0040   setEtaRange(double low, double high)
0041   {
0042     m_etaRange.first = low;
0043     m_etaRange.second = high;
0044   }
0045   //! set eta range
0046   void
0047   setPtRange(double low, double high)
0048   {
0049     m_ptRange.first = low;
0050     m_ptRange.second = high;
0051   }
0052
0053   void
0054   setMindR(double jetradius)
0055   {
0056       m_trackJetMatchingRadius = jetradius;
0057   }
0058   void use_initial_vertex(const bool b = true) {initial_vertex = b;}
0059   int Init(PHCompositeNode *topNode);
0060   int InitRun(PHCompositeNode *topNode);
0061   int process_event(PHCompositeNode *topNode);
0062   int process_event_bimp(PHCompositeNode *topNode);
0063   int End(PHCompositeNode *topNode);
0064   void initializeTrees();
0065   void Clean();
0066
0067  private:
0068   //! cache the jet evaluation modules
0069   std::shared_ptr<JetEvalStack> m_jetEvalStack;
0070   FastJetAlgo * m_fastjet;
0071   std::vector<Jet*> m_inputs_smallR;
0072   std::vector<Jet*> m_inputs_reco_smallR;
0073   std::vector<Jet*> m_inputs_recoMatch_smallR;
0074   std::vector<Jet*> m_output_largeR;
0075
0076   std::string m_recoJetName;
0077   std::string m_truthJetName;
0078   std::string m_outputFileName;
0079
0080   std::vector<fastjet::PseudoJet> m_truth_pseudojet;
0081
0082   //! eta range
0083   std::pair<double, double> m_etaRange;
0084
0085   //! pT range
0086   std::pair<double, double> m_ptRange;
0087
0088   //! flag to use initial vertex in track evaluator
0089   bool initial_vertex = false;
0090
0091   //! max track-jet matching radius
0092   double m_trackJetMatchingRadius;
0093
0094   //! Output histograms
0095   TH1 *m_hInclusiveE;
0096   TH1 *m_hInclusiveEta;
0097   TH1 *m_hInclusivePhi;
0098
0099   //! Output Tree variables
0100   TTree *m_T;
0101   TTree *m_Treclus;
0102   //TTree *m_TreEcoor;
0103   std::vector<double> m_recotruthE;
0104   std::vector<double> m_recotruthEta;
0105   std::vector<double> m_recotruthPhi;
0106   std::vector<double> m_recotruthPt;
0107   std::vector<double> m_recotruthPx;
0108   std::vector<double> m_recotruthPy;
0109   std::vector<double> m_recotruthPz;
0110   std::vector<double> m_recotruthnComponent;
0111
0112   int m_event;
0113   std::vector<double> m_id;
0114   std::vector<double> m_nComponent;
0115   std::vector<double> m_eta;
0116   std::vector<double> m_phi;
0117   std::vector<double> m_e;
0118   std::vector<double> m_new_e;
0119   std::vector<double> m_pt;
0120   std::vector<double> m_px;
0121   std::vector<double> m_py;
0122   std::vector<double> m_pz;
0123   std::vector<double> m_dR;
0124   std::vector<double> m_cent;
0125   std::vector<double> m_CAL_ID;
0126   std::vector<double> m_Constit_E;
0127   std::vector<double> m_Constit_Cent;
0128   //std::vector<double> m_Embedded_Count;
0129   //std::vector<double> m_Embedded_IHCAL;
0130   //std::vector<double> m_Embedded_CEMC;
0131   //std::vector<double> m_Embedded_OHCAL;
0132   //std::vector<double> m_Background_IHCAL;
0133   //std::vector<double> m_Background_CEMC;
0134   //std::vector<double> m_Background_OHCAL;
0135   //std::vector<double> m_Unidentified_OHCAL;
0136   //std::vector<double> m_Unidentified_IHCAL;
0137   //std::vector<double> m_Total_Count;
0138   std::vector<double> towers_id;
0139   std::vector<double> towers_primary;
0140   std::vector<double> m_IHCAL_Tower_Energy;
0141   std::vector<double> m_IHCAL_Cent;
0142   std::vector<double> m_EMCAL_Tower_Energy;
0143   std::vector<double> m_EMCAL_Cent;
0144   std::vector<double> m_OHCAL_Tower_Energy;
0145   std::vector<double> m_OHCAL_Cent;
0146
0147   double_t dPhi;
0148   double_t dPhi_temp;
0149
0150   CaloEvalStack* _caloevalstackHCALIN;
0151   CaloEvalStack* _caloevalstackHCALOUT;
0152   CaloEvalStack* _caloevalstackCEMC;
0153
0154   float temp_eta;
0155   float temp_phi;
0156   float temp_e;
0157   float temp_dR;
0158
0159   float temp_E_Matched;
0160   float temp_Phi_Matched;
0161   float temp_Eta_Matched;
0162
0163   float temp_TE_Matched;
0164   float temp_TPhi_Matched;
0165   float temp_TEta_Matched;
0166
0167   std::vector<double> m_truthNComponent;
0168   std::vector<double> m_truthEta;
0169   std::vector<double> m_truthPhi;
0170   std::vector<double> m_truthE;
0171   std::vector<double> m_truthPt;
0172   std::vector<double> m_truthdR;
0173   std::vector<double> m_truthPx;
0174   std::vector<double> m_truthPy;
0175   std::vector<double> m_truthPz;
0176   std::vector<double> m_truthConstitID;
0177   std::vector<double> m_truthConstitE;
0178   std::vector<double> m_truthConstitPt;
0179
0180   std::vector<double> m_E_Matched;
0181   std::vector<double> m_Phi_Matched;
0182   std::vector<double> m_Eta_Matched;
0183
0184   std::vector<double> m_TE_Matched;
0185   std::vector<double> m_TPhi_Matched;
0186   std::vector<double> m_TEta_Matched;
0187
0188   //std::vector<double> m_PairedID_Truth;
0189   //std::vector<double> m_PairedID_Reco;
0190   //std::vector<double> m_impactparam;
0191
0192   //! number of matched tracks
0193   int m_nMatchedTrack;
0194
0195   enum
0196   {
0197     //! max number of tracks
0198     kMaxMatchedTrack = 1000
0199   };
0200   std::array<float, kMaxMatchedTrack> m_trackdR;
0201   std::array<float, kMaxMatchedTrack> m_trackpT;
0202   std::array<float, kMaxMatchedTrack> m_trackPID;
0203 };
0204
0205 #endif  // MYJETANALYSIS_H