CalorimeterJets/src/MyJetAnalysis.cc

0001 #include "MyJetAnalysis.h"
0002 #include <vector>
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004 #include <fun4all/Fun4AllServer.h>
0005 #include <fun4all/PHTFileServer.h>
0006 //#include <sPHAnalysis.h>
0007
0008 #include <HepMC/GenVertex.h>
0009 #include <HepMC/IteratorRange.h>
0010 #include <HepMC/SimpleVector.h>
0011
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/getClass.h>
0014
0015 //#include <g4eval/JetEvalStack.h>
0016
0017 #include <trackbase_historic/SvtxTrackMap.h>
0018 #include <g4jets/FastJetAlgo.h>
0019 #include <jetbackground/FastJetAlgoSub.h>
0020 #include <g4jets/JetMap.h>
0021 #include <g4jets/Jetv1.h>
0022
0023 #include <g4vertex/GlobalVertex.h>
0024 #include <g4vertex/GlobalVertexMap.h>
0025
0026 #include <calobase/RawTowerContainer.h>
0027 #include <calobase/RawTower.h>
0028 #include <calobase/RawTowerGeomContainer.h>
0029 #include <calobase/RawTowerGeom.h>
0030 #include <calobase/TowerInfoContainer.h>
0031 #include <calobase/TowerInfo.h>
0032
0033 #include <g4main/PHG4Particle.h>
0034 #include <g4main/PHG4TruthInfoContainer.h>
0035
0036 #include <g4eval/CaloRawTowerEval.h>
0037 #include <g4eval/CaloEvalStack.h>
0038
0039 #include <ffaobjects/EventHeader.h>
0040 #include <ffaobjects/EventHeaderv1.h>
0041 #include <centrality/CentralityInfo.h>
0042 // fastjet includes
0043 #include <fastjet/ClusterSequence.hh>
0044 #include <fastjet/JetDefinition.hh>
0045 #include <fastjet/PseudoJet.hh>
0046 #include <jetbackground/TowerBackground.h>
0047 #include <phhepmc/PHHepMCGenEvent.h>
0048 #include <phhepmc/PHHepMCGenEventMap.h>
0049
0050 #include <TFile.h>
0051 #include <TH1F.h>
0052 #include <TH2F.h>
0053 #include <TString.h>
0054 #include <TTree.h>
0055 #include <TVector3.h>
0056 #include <TLorentzVector.h>
0057 #include <TArrayF.h>
0058
0059 #include <fstream>
0060 #include <algorithm>
0061 #include <cassert>
0062 #include <cmath>
0063 #include <iostream>
0064 #include <limits>
0065 #include <stdexcept>
0066 //#include <memory>
0067
0068 double DeltaPhi(double phi1, double phi2){
0069    double_t dPhi_temp = phi1 - phi2;
0070    if(dPhi_temp<-TMath::Pi()) dPhi_temp += TMath::TwoPi();
0071    if(dPhi_temp>=TMath::Pi()) dPhi_temp -= TMath::TwoPi();
0072    return dPhi_temp;
0073    }
0074
0075 //#include "/sphenix/user/stapiaar/jet_correlator/EnergyEnergyCorrelators/eec/include/EEC.hh"
0076
0077 using namespace std;
0078
0079 MyJetAnalysis::MyJetAnalysis(const std::string& recojetname, const std::string& truthjetname, const std::string& outputfilename)
0080   : SubsysReco("MyJetAnalysis_" + recojetname + "_" + truthjetname)
0081   , m_recoJetName(recojetname)
0082   , m_truthJetName(truthjetname)
0083   , m_outputFileName(outputfilename)
0084   , m_etaRange(-1.1, 1.1)
0085   , m_ptRange(5, 100)
0086   , m_trackJetMatchingRadius(.7)
0087   , m_hInclusiveE(nullptr)
0088   , m_hInclusiveEta(nullptr)
0089   , m_hInclusivePhi(nullptr)
0090   , m_T(nullptr)
0091   , m_event(-1)
0092   , _caloevalstackHCALIN(nullptr)
0093   , _caloevalstackHCALOUT(nullptr)
0094   , _caloevalstackCEMC(nullptr)
0095   /* , m_id(-1)
0096   , m_nComponent(-1)
0097   , m_eta(numeric_limits<float>::signaling_NaN())
0098   , m_phi(numeric_limits<float>::signaling_NaN())
0099  // , m_e(numeric_limits<vector>::signaling_NaN())
0100   , m_pt(numeric_limits<float>::signaling_NaN())
0101  // , m_px(numeric_limits<vector>::signaling_NaN())
0102  // , m_py(numeric_limits<vector>::signaling_NaN())
0103  // , m_pz(numeric_limits<vector>::signaling_NaN())
0104     //m_v4(numeric_limits<float>::signaling_NaN())
0105   , //m_truthID(-1)
0106   , m_truthNComponent(-1)
0107   , m_truthEta(numeric_limits<float>::signaling_NaN())
0108   , m_truthPhi(numeric_limits<float>::signaling_NaN())
0109  // , m_truthE(numeric_limits<vector>::signaling_NaN())
0110   , m_truthPt(numeric_limits<float>::signaling_NaN())
0111  // , m_truthPx(numeric_limits<vector>::signaling_NaN())
0112  // , m_truthPy(numeric_limits<vector>::signaling_NaN())
0113  /  , m_truthPz(numeric_limits<vecto>::signaling_NaN())
0114   */
0115   , m_nMatchedTrack(-1)
0116
0117 {
0118   m_trackdR.fill(numeric_limits<float>::signaling_NaN());
0119   m_trackpT.fill(numeric_limits<float>::signaling_NaN());
0120   m_trackPID.fill(numeric_limits<float>::signaling_NaN());
0121 }
0122
0123 MyJetAnalysis::~MyJetAnalysis()
0124 {
0125 }
0126
0127 int MyJetAnalysis::Init(PHCompositeNode* topNode)
0128 {
0129   if (Verbosity() >= MyJetAnalysis::VERBOSITY_SOME)
0130     cout << "MyJetAnalysis::Init - Outoput to " << m_outputFileName << endl;
0131
0132   PHTFileServer::get().open(m_outputFileName, "RECREATE");
0133
0134   cout << "MyJetAnalysis::Init - Outoput to " << m_outputFileName << endl;
0135
0136
0137   // Histograms
0138   m_hInclusiveE = new TH1F(
0139       "hInclusive_E",  //
0140       TString(m_recoJetName) + " inclusive jet E;Total jet energy (GeV)", 100, 0, 100);
0141
0142   m_hInclusiveEta =
0143       new TH1F(
0144           "hInclusive_eta",  //
0145           TString(m_recoJetName) + " inclusive jet #eta;#eta;Jet energy density", 50, -1, 1);
0146   m_hInclusivePhi =
0147       new TH1F(
0148           "hInclusive_phi",  //
0149           TString(m_recoJetName) + " inclusive jet #phi;#phi;Jet energy density", 50, -M_PI, M_PI);
0150
0151   initializeTrees();
0152   return Fun4AllReturnCodes::EVENT_OK;
0153 }
0154
0155 int MyJetAnalysis::End(PHCompositeNode* topNode)
0156 {
0157   PHTFileServer::get().cd(m_outputFileName);
0158   m_T->Write();
0159   cout << "MyJetAnalysis::End - Output to " << m_outputFileName << endl;
0160   return Fun4AllReturnCodes::EVENT_OK;
0161 }
0162
0163 int MyJetAnalysis::InitRun(PHCompositeNode* topNode)
0164 {
0165 return Fun4AllReturnCodes::EVENT_OK;
0166 }
0167
0168 int MyJetAnalysis::process_event(PHCompositeNode* topNode)
0169 {
0170   //cout << "MyJetAnalysis::process_event" << endl;
0171   if (Verbosity() >= MyJetAnalysis::VERBOSITY_SOME)
0172     cout << "MyJetAnalysis::process_event() entered" << endl;
0173
0174   ++m_event;
0175   if( (m_event % 10) == 0 ) cout << "Event number = "<< m_event << endl;
0176
0177   // Making Jets (Needed if not using JetMap)
0178   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0179
0180   // interface to jets (JetMap)
0181   //JetMap* jets = findNode::getClass<JetMap>(topNode, m_recoJetName);
0182   JetMap* jetsMC = findNode::getClass<JetMap>(topNode, m_truthJetName);
0183
0184   // Making Jets (Non JetMap)
0185
0186   RawTowerContainer *ihcal_towers_sub = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN_SUB1");
0187   RawTowerContainer *ihcal_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
0188   RawTowerContainer *cemc_towers_sub = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER_SUB1");
0189   RawTowerContainer *cemc_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
0190   RawTowerContainer *ohcal_towers_sub = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT_SUB1");
0191   RawTowerContainer *ohcal_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
0192
0193   //TowerInfoContainer *ihcal_towers_sub = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN_SUB1");
0194   //TowerInfoContainer *ihcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0195   //TowerInfoContainer *cemc_towers_sub = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER_SUB1");
0196   //TowerInfoContainer *cemc_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0197   //TowerInfoContainer *ohcal_towers_sub = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT_SUB1");
0198   //TowerInfoContainer *ohcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0199
0200   RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0201   RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0202   RawTowerGeomContainer *geomEC = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0203
0204   if (!geomIH || !geomOH ){
0205      cout << "Cannot find Tower Geom Nodes" << endl;
0206      exit(-1);
0207   }
0208
0209   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0210   if (!truthinfo){
0211      cout << "Error cant find G4 Truth Info" << endl;
0212      exit(-1);
0213   }
0214
0215    CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0216   if (!cent_node){
0217      cout << "Error can't find CentralityInfo" << endl;
0218      exit(-1);
0219   }
0220
0221   float cent_mbd = cent_node->get_centile(CentralityInfo::PROP::mbd_NS);
0222   m_cent.push_back(cent_mbd);
0223
0224   /*
0225   std::vector<fastjet::PseudoJet> pseudojet_truth;
0226   PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
0227      for(PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter){
0228         PHG4Particle *part = iter->second;
0229         if ((abs(part->get_pid()) >= 12) && (abs(part->get_pid()) <= 16)) continue;
0230         if ((part->get_px() == 0.0) && (part->get_py() == 0.0)) continue;
0231         float eta = asinh(part->get_pz()/sqrt(pow(part->get_px(),2) + pow(part->get_py(),2)));
0232         if (eta < -1.1) continue;
0233         if (eta > 1.1) continue;
0234         float Part_ID = truthinfo->isEmbeded(part->get_track_id());
0235         if (Part_ID != 2) continue;
0236         m_truthConstitID.push_back(Part_ID);
0237         m_truthConstitE.push_back(part->get_e());
0238         //m_truthConstitPt.push_back(part->get_pt());
0239         if (Part_ID != 2){
0240         cout << "TruthJet has a particle with ID = " << Part_ID << endl;
0241         }
0242         fastjet::PseudoJet jet_truth(part->get_px(), part->get_py(), part->get_pz(), part->get_e());
0243         //jet_truth.set_user_index(Part_ID);
0244         pseudojet_truth.push_back(jet_truth);
0245      }
0246   */
0247   GlobalVertex *vtx = vertexmap->begin()->second;
0248   float vtxz = NAN;
0249   if (vtx) vtxz = vtx->get_z();
0250
0251   std::vector<fastjet::PseudoJet> pseudojet_reco;
0252   std::vector<fastjet::PseudoJet> pseudojet_reco_truth;
0253   CaloRawTowerEval* towerevalHCALIN = new CaloRawTowerEval(topNode, "HCALIN");
0254   CaloRawTowerEval* towerevalHCALOUT = new CaloRawTowerEval(topNode, "HCALOUT");
0255   CaloRawTowerEval* towerevalCEMC = new CaloRawTowerEval(topNode, "CEMC");
0256
0257   RawTowerContainer::ConstRange ICB_begin_end = ihcal_towers->getTowers();
0258   RawTowerContainer::ConstIterator ICB_rtiter;
0259   for (ICB_rtiter = ICB_begin_end.first; ICB_rtiter != ICB_begin_end.second; ++ICB_rtiter){
0260      float Tower_Embedded = 0;
0261      RawTower *towerB = ICB_rtiter->second;
0262      if (towerB->get_energy() <= 0) continue;
0263      PHG4Particle* primary = towerevalHCALIN->max_truth_primary_particle_by_energy(towerB);
0264      if (primary){
0265         if (truthinfo->isEmbeded(primary->get_track_id()) == 2){
0266            Tower_Embedded = 0;
0267            }
0268         else{
0269            Tower_Embedded = 1;
0270            }
0271         }
0272      else Tower_Embedded = -10;
0273      towers_id.push_back(towerB->get_id());
0274      towers_primary.push_back(Tower_Embedded);
0275   }
0276
0277   RawTowerContainer::ConstRange IC_begin_end = ihcal_towers_sub->getTowers();
0278   RawTowerContainer::ConstIterator IC_rtiter;
0279   for (IC_rtiter = IC_begin_end.first; IC_rtiter != IC_begin_end.second; ++IC_rtiter){
0280      RawTower *tower = IC_rtiter->second;
0281      RawTowerGeom *tower_geom = geomIH->get_tower_geometry(tower->get_key());
0282      assert(tower_geom);
0283      if(tower->get_energy() <= 0) continue;
0284      float HCALIN_Embedded = -2;
0285      //for (ICB_rtiter = ICB_begin_end.first; ICB_rtiter != ICB_begin_end.second; ++ICB_rtiter){
0286         //RawTower *towerB = ICB_rtiter->second;
0287         //if (towerB->get_id() != tower->get_id()) continue;
0288      for (int j = 0; j < (int)towers_id.size(); ++j){
0289         if (tower->get_id() != towers_id.at(j)) continue;
0290         HCALIN_Embedded = towers_primary.at(j);
0291      }
0292
0293      double r = tower_geom->get_center_radius();
0294      double phi = atan2(tower_geom->get_center_y(), tower_geom->get_center_x());
0295      double z0 = tower_geom->get_center_z();
0296      double z = z0 - vtxz;
0297
0298      double eta = asinh(z/r);
0299      double pt = tower->get_energy()/cosh(eta);
0300      double px = pt*cos(phi);
0301      double py = pt*sin(phi);
0302      double pz = pt*sinh(eta);
0303
0304      m_IHCAL_Tower_Energy.push_back(tower->get_energy());
0305      m_IHCAL_Cent.push_back(cent_mbd);
0306      //if (HCALIN_Embedded == -10 || HCALIN_Embedded == 1) continue;
0307      fastjet::PseudoJet pseudojet(px, py, pz, tower->get_energy());
0308      pseudojet.set_user_index(HCALIN_Embedded);
0309      pseudojet_reco.push_back(pseudojet);
0310      if (HCALIN_Embedded == -10 || HCALIN_Embedded == 1) continue;
0311      pseudojet_reco_truth.push_back(pseudojet);
0312   }
0313
0314   towers_id.clear();
0315   towers_primary.clear();
0316
0317   RawTowerContainer::ConstRange EC_begin_end = cemc_towers_sub->getTowers();
0318   RawTowerContainer::ConstRange ECB_begin_end = cemc_towers->getTowers();
0319   RawTowerContainer::ConstIterator EC_rtiter;
0320   RawTowerContainer::ConstIterator ECB_rtiter;
0321   for (EC_rtiter = EC_begin_end.first; EC_rtiter != EC_begin_end.second; ++EC_rtiter){
0322      RawTower *tower = EC_rtiter->second;
0323      RawTowerGeom *tower_geom = geomIH->get_tower_geometry(tower->get_key());
0324      assert(tower_geom);
0325      if(tower->get_energy() <= 0) continue;
0326      int this_ECetabin = geomIH->get_etabin(tower_geom->get_eta());
0327      int this_ECphibin = geomIH->get_phibin(tower_geom->get_phi());
0328      float HCEMC_Embedded = -2;
0329      float Tower_Energy = 0;
0330      for (ECB_rtiter = ECB_begin_end.first; ECB_rtiter != ECB_begin_end.second; ++ECB_rtiter){
0331         RawTower *towerB = ECB_rtiter->second;
0332         RawTowerGeom *towerB_geom = geomEC->get_tower_geometry(towerB->get_key());
0333         int this_ECBetabin = geomIH->get_etabin(towerB_geom->get_eta());
0334         int this_ECBphibin = geomIH->get_phibin(towerB_geom->get_phi());
0335         if (this_ECBetabin != this_ECetabin && this_ECBphibin != this_ECphibin) continue;
0336         if(towerB->get_energy() < Tower_Energy) continue;
0337         PHG4Particle* primaryEC = towerevalCEMC->max_truth_primary_particle_by_energy(towerB);
0338         if (primaryEC){
0339            if (truthinfo->isEmbeded(primaryEC->get_track_id()) == 2){
0340            HCEMC_Embedded = 2;
0341            }
0342            else{
0343            HCEMC_Embedded = 3;
0344            }
0345         }
0346         else{
0347            HCEMC_Embedded = -11;
0348         }
0349         Tower_Energy = towerB->get_energy();
0350      }
0351
0352
0353      double r = tower_geom->get_center_radius();
0354      double phi = atan2(tower_geom->get_center_y(), tower_geom->get_center_x());
0355      double z0 = tower_geom->get_center_z();
0356      double z = z0 - vtxz;
0357
0358      double eta = asinh(z/r);
0359      double pt = tower->get_energy()/cosh(eta);
0360      double px = pt*cos(phi);
0361      double py = pt*sin(phi);
0362      double pz = pt*sinh(eta);
0363
0364      m_EMCAL_Tower_Energy.push_back(tower->get_energy());
0365      m_EMCAL_Cent.push_back(cent_mbd);
0366      //if(HCEMC_Embedded == -11  || HCEMC_Embedded == 3) continue;
0367      fastjet::PseudoJet pseudojet(px, py, pz, tower->get_energy());
0368      pseudojet.set_user_index(HCEMC_Embedded);
0369      pseudojet_reco.push_back(pseudojet);
0370      if(HCEMC_Embedded == -11  || HCEMC_Embedded == 3) continue;
0371      pseudojet_reco_truth.push_back(pseudojet);
0372   }
0373
0374   RawTowerContainer::ConstRange OCB_begin_end = ohcal_towers->getTowers();
0375   RawTowerContainer::ConstIterator OCB_rtiter;
0376   for (OCB_rtiter = OCB_begin_end.first; OCB_rtiter != OCB_begin_end.second; ++OCB_rtiter){
0377      float Tower_Embedded = 0;
0378      RawTower *towerB = OCB_rtiter->second;
0379      if (towerB->get_energy() <= 0) continue;
0380      PHG4Particle* primary = towerevalHCALOUT->max_truth_primary_particle_by_energy(towerB);
0381      if (primary){
0382         if (truthinfo->isEmbeded(primary->get_track_id()) == 2){
0383            Tower_Embedded = 4;
0384            }
0385         else{
0386            Tower_Embedded = 5;
0387            }
0388         }
0389      else Tower_Embedded = -12;
0390      towers_id.push_back(towerB->get_id());
0391      towers_primary.push_back(Tower_Embedded);
0392   }
0393
0394   RawTowerContainer::ConstRange OC_begin_end = ohcal_towers_sub->getTowers();
0395   RawTowerContainer::ConstIterator OC_rtiter;
0396   for (OC_rtiter = OC_begin_end.first; OC_rtiter != OC_begin_end.second; ++OC_rtiter){
0397      RawTower *tower = OC_rtiter->second;
0398      RawTowerGeom *tower_geom = geomOH->get_tower_geometry(tower->get_key());
0399      assert(tower_geom);
0400      if(tower->get_energy() <= 0) continue;
0401      float HCALOUT_Embedded = -1;
0402      for (int j = 0; j < (int)towers_id.size(); ++j){
0403         if (tower->get_id() != towers_id.at(j)) continue;
0404         HCALOUT_Embedded = towers_primary.at(j);
0405      }
0406      double r = tower_geom->get_center_radius();
0407      double phi = atan2(tower_geom->get_center_y(), tower_geom->get_center_x());
0408      double z0 = tower_geom->get_center_z();
0409      double z = z0 - vtxz;
0410
0411      double eta = asinh(z/r);
0412      double pt = tower->get_energy()/cosh(eta);
0413      double px = pt*cos(phi);
0414      double py = pt*sin(phi);
0415      double pz = pt*sinh(eta);
0416
0417      m_OHCAL_Tower_Energy.push_back(tower->get_energy());
0418      m_OHCAL_Cent.push_back(cent_mbd);
0419      //if (HCALOUT_Embedded == -12 || HCALOUT_Embedded == 5) continue;
0420      fastjet::PseudoJet pseudojet(px, py, pz, tower->get_energy());
0421      pseudojet.set_user_index(HCALOUT_Embedded);
0422      pseudojet_reco.push_back(pseudojet);
0423      if (HCALOUT_Embedded == -12 || HCALOUT_Embedded == 5) continue;
0424      pseudojet_reco_truth.push_back(pseudojet);
0425   }
0426
0427   towers_id.clear();
0428   towers_primary.clear();
0429
0430   fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, 0.2, fastjet::E_scheme,fastjet::Best);
0431   fastjet::ClusterSequence jetFinder_reco_truth(pseudojet_reco_truth, jetDef);
0432   fastjet::ClusterSequence jetFinder_reco(pseudojet_reco, jetDef);
0433   std::vector<fastjet::PseudoJet> fastjets_reco_truth = jetFinder_reco_truth.inclusive_jets();
0434   std::vector<fastjet::PseudoJet> fastjets_reco = jetFinder_reco.inclusive_jets();
0435
0436   pseudojet_reco_truth.clear();
0437   pseudojet_reco.clear();
0438
0439   //Saves Truth Jet Info
0440   for (JetMap::Iter iter = jetsMC->begin(); iter != jetsMC->end(); ++iter){
0441   //for (int j = 0; j < (int)fastjets.size(); ++j){
0442     Jet* truthjet = iter->second;
0443     assert(truthjet);
0444     bool eta_cut = (truthjet->get_eta() >= m_etaRange.first) and (truthjet->get_eta() <= m_etaRange.second);
0445     if (not eta_cut) continue;
0446     //cout << "Jet Truth = " << j << endl;
0447     /*
0448     vector<fastjet::PseudoJet> consts = jetFinder.constituents(fastjets[j]);
0449     //cout << "Size of Truth Jet" << (int)consts.size() << endl;
0450     for (int i = 0; i < (int)consts.size(); i++){
0451        if (consts[i].user_index() != -10){
0452           //cout << consts[i].user_index() << endl;
0453        }
0454     }
0455
0456     for(Jet::ConstIter comp = truthjet->begin_comp(); comp != truthjet->end_comp(); ++comp){
0457        const int this_embed_id = truthinfo->isEmbeded((*comp).second);
0458        if (this_embed_id != 2){
0459           cout << "Truth jet has non embedded particle with id = " << this_embed_id << endl;
0460        }
0461     }
0462     */
0463     m_truthEta.push_back(truthjet->get_eta());
0464     m_truthPhi.push_back(truthjet->get_phi());
0465     m_truthE.push_back(truthjet->get_e());
0466     m_truthPt.push_back(truthjet->get_pt());
0467     m_truthPx.push_back(truthjet->get_px());
0468     m_truthPy.push_back(truthjet->get_py());
0469     m_truthPz.push_back(truthjet->get_pz());
0470     m_truthNComponent.push_back(truthjet->size_comp());
0471
0472   }
0473     //Saving Reco Truth Jets
0474     for (int j = 0; j < (int)fastjets_reco_truth.size(); ++j){
0475
0476        vector<fastjet::PseudoJet> consts_reco_truth = jetFinder_reco_truth.constituents(fastjets_reco_truth[j]);
0477
0478        m_recotruthnComponent.push_back((int)consts_reco_truth.size());
0479        m_recotruthEta.push_back(fastjets_reco_truth[j].eta());
0480        m_recotruthPhi.push_back(fastjets_reco_truth[j].phi_std());
0481        m_recotruthPt.push_back(fastjets_reco_truth[j].perp());
0482        m_recotruthPx.push_back(fastjets_reco_truth[j].px());
0483        m_recotruthPy.push_back(fastjets_reco_truth[j].py());
0484        m_recotruthPz.push_back(fastjets_reco_truth[j].pz());
0485        m_recotruthE.push_back(fastjets_reco_truth[j].e());
0486
0487     }
0488     //saving closest recojets
0489     //for (JetMap::Iter iter = jets->begin(); iter != jets->end(); ++iter)
0490     for (int j = 0; j < (int)fastjets_reco.size(); ++j)
0491     {
0492         vector<fastjet::PseudoJet> consts = jetFinder_reco.constituents(fastjets_reco[j]);
0493
0494         //Jet* jet = iter->second;
0495
0496         //Needed For Making Jets
0497         /*
0498         float nEmbedded = 0;
0499         float nTotal = 0;
0500         float nUnidentified_IHCAL = 0;
0501         float nUnidentified_OHCAL = 0;
0502         float nIHCAL_EM = 0;
0503         float nIHCAL_BC = 0;
0504         float nCEMC_EM = 0;
0505         float nCEMC_BC = 0;
0506         float nOHCAL_EM = 0;
0507         float nOHCAL_BC = 0;
0508         */
0509         //ID ordered by ID = 0 embedded Inner, ID = 1 nonembedded Inner, continues like this for EMCal then Outer
0510         for (int i = 0; i < (int)consts.size(); i++){
0511             m_CAL_ID.push_back(consts[i].user_index());
0512             m_Constit_E.push_back(consts[i].e());
0513             m_Constit_Cent.push_back(cent_mbd);
0514             /*
0515             if (consts[i].user_index() == 0 || consts[i].user_index() == 2 || consts[i].user_index() == 4){
0516                nEmbedded = nEmbedded + 1;
0517                if (consts[i].user_index() == 0){
0518                   nIHCAL_EM = nIHCAL_EM + 1;
0519                }
0520                if (consts[i].user_index() == 2){
0521                   nCEMC_EM = nCEMC_EM + 1;
0522                }
0523                if (consts[i].user_index() == 4){
0524                   nOHCAL_EM = nOHCAL_EM + 1;
0525                }
0526             }
0527             else{
0528                if (consts[i].user_index() == 1){
0529                   nIHCAL_BC = nIHCAL_BC + 1;
0530                }
0531                if (consts[i].user_index() == 3){
0532                   nCEMC_BC = nCEMC_BC + 1;
0533                }
0534                if (consts[i].user_index() == 5){
0535                   nOHCAL_BC = nOHCAL_BC + 1;
0536                }
0537             }
0538
0539             if (consts[i].user_index() == -11){
0540                nUnidentified_IHCAL = nUnidentified_IHCAL + 1;
0541             }
0542             if (consts[i].user_index() == -10){
0543                nUnidentified_OHCAL = nUnidentified_OHCAL + 1;
0544             }
0545             nTotal = nTotal + 1;
0546         */
0547         }
0548         /*
0549         m_Embedded_IHCAL.push_back(nIHCAL_EM);
0550         m_Embedded_CEMC.push_back(nCEMC_EM);
0551         m_Embedded_OHCAL.push_back(nOHCAL_EM);
0552         m_Background_IHCAL.push_back(nIHCAL_BC);
0553         m_Background_CEMC.push_back(nCEMC_BC);
0554         m_Background_OHCAL.push_back(nOHCAL_BC);
0555         m_Embedded_Count.push_back(nEmbedded);
0556         m_Unidentified_IHCAL.push_back(nUnidentified_IHCAL);
0557         m_Unidentified_OHCAL.push_back(nUnidentified_OHCAL);
0558         m_Total_Count.push_back(nTotal);
0559         */
0560         // Save reco jet information (In code jet reconstruction)
0561         m_eta.push_back(fastjets_reco[j].eta());
0562         m_phi.push_back(fastjets_reco[j].phi_std());
0563         m_pt.push_back(fastjets_reco[j].perp());
0564         m_px.push_back(fastjets_reco[j].px());
0565         m_py.push_back(fastjets_reco[j].py());
0566         m_pz.push_back(fastjets_reco[j].pz());
0567         m_e.push_back(fastjets_reco[j].e());
0568
0569         //Save the reco jet Information (Jet Map)
0570
0571         //m_id.push_back(jet->get_id());
0572         //m_nComponent.push_back(jet->size_comp());
0573         //m_eta.push_back(jet->get_eta());
0574         //m_phi.push_back(jet->get_phi());
0575         //m_pt.push_back(jet->get_pt());
0576         //m_px.push_back(jet->get_px());
0577         //m_py.push_back(jet->get_py());
0578         //m_pz.push_back(jet->get_pz());
0579         //m_e.push_back(jet->get_e());
0580
0581         float m_Loop_dR = 100000;// dR used for Loop to save closer truth jets
0582
0583         //for (int i = 0; i < (int)fastjets.size(); ++i){
0584         for(std::size_t i = 0, max = (int)m_truthPhi.size(); i != max; ++i){
0585          //if(m_truthE.at(i) < 0) continue;
0586
0587              double dEta_temp = fastjets_reco[j].eta() - m_truthEta.at(i);
0588              double dPhi_temp = DeltaPhi(fastjets_reco[j].phi_std(), m_truthPhi.at(i));
0589
0590              //double dEta_temp = jet->get_eta() - m_truthEta.at(i);
0591              //double dPhi_temp = DeltaPhi(jet->get_phi(), m_truthPhi.at(i));
0592
0593              float dR_temp = sqrt(dEta_temp * dEta_temp + dPhi_temp * dPhi_temp);
0594          if(dR_temp > m_Loop_dR) continue;
0595
0596              // Saving truth R = 0.2 matched Info
0597              temp_TE_Matched = m_truthE.at(i);
0598              temp_TPhi_Matched = m_truthPhi.at(i);
0599              temp_TEta_Matched = m_truthEta.at(i);
0600
0601
0602              m_Loop_dR = dR_temp;//Resets cutoff dR to mae sure to save only closer truth jets
0603           }
0604
0605       m_E_Matched.push_back(fastjets_reco[j].e());
0606       m_Phi_Matched.push_back(fastjets_reco[j].phi_std());
0607       m_Eta_Matched.push_back(fastjets_reco[j].eta());
0608
0609       m_TE_Matched.push_back(temp_TE_Matched);
0610       m_TPhi_Matched.push_back(temp_TPhi_Matched);
0611       m_TEta_Matched.push_back(temp_TEta_Matched);
0612       m_dR.push_back(m_Loop_dR);
0613     }
0614
0615    m_T->Fill();
0616    m_IHCAL_Tower_Energy.clear();
0617    m_IHCAL_Cent.clear();
0618    m_EMCAL_Tower_Energy.clear();
0619    m_EMCAL_Cent.clear();
0620    m_OHCAL_Tower_Energy.clear();
0621    m_OHCAL_Cent.clear();
0622    //clearing truth info
0623    //m_truthID.clear();
0624    //m_truthNComponent.clear();
0625    m_truthEta.clear();
0626    m_truthPhi.clear();
0627    m_truthE.clear();
0628    m_truthPt.clear();
0629    m_truthPx.clear();
0630    m_truthPy.clear();
0631    m_truthPz.clear();
0632    //m_truthConstitE.clear();
0633    //m_truthConstitID.clear();
0634    //m_truthConstitPt.clear();
0635    m_cent.clear();
0636
0637    //clearing reco info
0638    //m_id.clear();
0639    m_nComponent.clear();
0640    m_eta.clear();
0641    m_phi.clear();
0642    m_pt.clear();
0643    m_e.clear();
0644    m_px.clear();
0645    m_py.clear();
0646    m_pz.clear();
0647
0648    m_recotruthnComponent.clear();
0649    m_recotruthEta.clear();
0650    m_recotruthPhi.clear();
0651    m_recotruthPt.clear();
0652    m_recotruthE.clear();
0653    m_recotruthPx.clear();
0654    m_recotruthPy.clear();
0655    m_recotruthPz.clear();
0656
0657    //MAtched infor clear
0658    m_TE_Matched.clear();
0659    m_TPhi_Matched.clear();
0660    m_TEta_Matched.clear();
0661
0662    m_E_Matched.clear();
0663    m_Phi_Matched.clear();
0664    m_Eta_Matched.clear();
0665
0666    m_dR.clear();
0667    m_cent.clear();
0668    m_CAL_ID.clear();
0669    m_Constit_E.clear();
0670    m_Constit_Cent.clear();
0671    /*
0672    m_Embedded_IHCAL.clear();
0673    m_Embedded_CEMC.clear();
0674    m_Embedded_OHCAL.clear();
0675    m_Background_IHCAL.clear();
0676    m_Background_CEMC.clear();
0677    m_Background_OHCAL.clear();
0678    m_Embedded_Count.clear();
0679    m_Total_Count.clear();
0680    m_Unidentified_IHCAL.clear();
0681    m_Unidentified_OHCAL.clear();
0682    */
0683    Clean();
0684
0685    return Fun4AllReturnCodes::EVENT_OK;
0686 }
0687
0688 void MyJetAnalysis::initializeTrees()
0689 {
0690
0691 m_T = new TTree("T", "MyJetAnalysis Tree");
0692
0693 //Storing All Reco Info
0694 m_T->Branch("m_event", &m_event, "event/I");
0695 m_T->Branch("id", &m_id);
0696 m_T->Branch("nComponent", &m_nComponent);
0697 m_T->Branch("eta", &m_eta);
0698 m_T->Branch("phi", &m_phi);
0699 m_T->Branch("e", &m_e);
0700 m_T->Branch("pt", &m_pt);
0701 m_T->Branch("px", &m_px);
0702 m_T->Branch("py", &m_py);
0703 m_T->Branch("pz", &m_pz);
0704
0705 m_T->Branch("IHCAL_Tower_Energy",&m_IHCAL_Tower_Energy);
0706 m_T->Branch("IHCAL_Cent",&m_IHCAL_Cent);
0707 m_T->Branch("EMCAL_Tower_Energy",&m_EMCAL_Tower_Energy);
0708 m_T->Branch("EMCAL_Cent",&m_EMCAL_Cent);
0709 m_T->Branch("OHCAL_Tower_Energy",&m_OHCAL_Tower_Energy);
0710 m_T->Branch("OHCAL_Cent",&m_OHCAL_Cent);
0711
0712
0713 //Storing All Reco Truth Info
0714 m_T->Branch("recotruthnComponent", &m_recotruthnComponent);
0715 m_T->Branch("recotruthEta", &m_recotruthEta);
0716 m_T->Branch("recotruthPhi", &m_recotruthPhi);
0717 m_T->Branch("recotruthPt", &m_recotruthPt);
0718 m_T->Branch("recotruthPx", &m_recotruthPx);
0719 m_T->Branch("recotruthPy", &m_recotruthPy);
0720 m_T->Branch("recotruthPz", &m_recotruthPz);
0721 m_T->Branch("recotruthE", &m_recotruthE);
0722
0723 //Storing Matched Jets Info
0724 m_T->Branch("dR", &m_dR);
0725
0726 //Storing Reco Jet Info
0727 m_T->Branch("cent", &m_cent);
0728 m_T->Branch("CAL_ID", &m_CAL_ID);
0729 //m_T->Branch("CAL_ID", &m_CAL_ID_RT);
0730 m_T->Branch("Constit_Energy", &m_Constit_E);
0731 //m_T->Branch("Constit_Energy_RT", &m_Constit_RT);
0732 m_T->Branch("Constit_Cent", &m_Constit_Cent);
0733 //m_T->Branch("Constit_Cent_RT", &m_Constit_Cent_RT);
0734 //m_T->Branch("Embedded_Count", &m_Embedded_Count);
0735 //m_T->Branch("Embedded_nIHCAL", &m_Embedded_IHCAL);
0736 //m_T->Branch("Embedded_nCEMC", &m_Embedded_CEMC);
0737 //m_T->Branch("Embedded_nOHCAL", &m_Embedded_OHCAL);
0738 //m_T->Branch("Background_nIHCAL", &m_Background_IHCAL);
0739 //m_T->Branch("Background_nCEMC", &m_Background_CEMC);
0740 //m_T->Branch("Background_nOHCAL", &m_Background_OHCAL);
0741 //m_T->Branch("Unidentified_IHCAL", &m_Unidentified_IHCAL);
0742 //m_T->Branch("Unidentified_OHCAL", &m_Unidentified_OHCAL);
0743 //m_T->Branch("Total_Count", &m_Total_Count);
0744 //m_T->Branch("truthID", &m_truthID);
0745 //m_T->Branch("truthNComponent", &m_truthNComponent);
0746 //m_T->Branch("truthConstitE", &m_truthConstitE);
0747 //m_T->Branch("truthConstitID", &m_truthConstitID);
0748 m_T->Branch("truthEta", &m_truthEta);
0749 m_T->Branch("truthPhi", &m_truthPhi);
0750 m_T->Branch("truthE", &m_truthE);
0751 m_T->Branch("truthPt", &m_truthPt);
0752 m_T->Branch("truthPx", &m_truthPx);
0753 m_T->Branch("truthPy", &m_truthPy);
0754 m_T->Branch("truthPz", &m_truthPz);
0755 //m_T->Branch("truthdR", &m_truthdR, "truthdR/F");
0756
0757 }
0758
0759  void MyJetAnalysis::Clean()
0760 {
0761 //m_cent.clear();
0762 }
0763
0764
0765
0766
0767
0768
0769