00001 #ifndef CalibrationIsolatedParticlesIsolatedGenParticles_h
00002 #define CalibrationIsolatedParticlesIsolatedGenParticles_h
00003
00004
00005 #include <memory>
00006
00007
00008 #include "FWCore/Framework/interface/Frameworkfwd.h"
00009 #include "FWCore/Framework/interface/EDAnalyzer.h"
00010
00011 #include "FWCore/Framework/interface/Event.h"
00012 #include "FWCore/Framework/interface/MakerMacros.h"
00013
00014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00015
00016 #include "FWCore/ServiceRegistry/interface/Service.h"
00017 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00018
00019
00020 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00021 #include "SimTracker/Records/interface/TrackAssociatorRecord.h"
00022 #include "DataFormats/GeometrySurface/interface/GloballyPositioned.h"
00023
00024
00025 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00026 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00027
00028 #include "DataFormats/DetId/interface/DetId.h"
00029 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00030 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00031 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
00032
00033 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00034 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
00035
00036 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
00037
00038
00039 #include "DataFormats/L1Trigger/interface/L1JetParticle.h"
00040 #include "DataFormats/L1Trigger/interface/L1JetParticleFwd.h"
00041 #include "DataFormats/L1Trigger/interface/L1EmParticle.h"
00042 #include "DataFormats/L1Trigger/interface/L1EmParticleFwd.h"
00043 #include "DataFormats/L1Trigger/interface/L1MuonParticle.h"
00044 #include "DataFormats/L1Trigger/interface/L1MuonParticleFwd.h"
00045
00046 #include "CondFormats/L1TObjects/interface/L1GtTriggerMenu.h"
00047 #include "CondFormats/DataRecord/interface/L1GtTriggerMenuRcd.h"
00048 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
00049 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
00050 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
00051 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
00052 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapFwd.h"
00053 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMap.h"
00054
00055 #include "DataFormats/Math/interface/LorentzVector.h"
00056 #include "Calibration/IsolatedParticles/interface/GenSimInfo.h"
00057
00058
00059 #include "TROOT.h"
00060 #include "TSystem.h"
00061 #include "TFile.h"
00062 #include "TH1I.h"
00063 #include "TH2D.h"
00064 #include "TProfile.h"
00065 #include "TDirectory.h"
00066 #include "TTree.h"
00067
00068
00069 namespace{
00070 class ParticlePtGreater{
00071 public:
00072 int operator()(const HepMC::GenParticle * p1,
00073 const HepMC::GenParticle * p2) const{
00074 return p1->momentum().perp() > p2->momentum().perp();
00075 }
00076 };
00077
00078 class ParticlePGreater{
00079 public:
00080 int operator()(const HepMC::GenParticle * p1,
00081 const HepMC::GenParticle * p2) const{
00082 return p1->momentum().rho() > p2->momentum().rho();
00083 }
00084 };
00085 }
00086
00087
00088 class IsolatedGenParticles : public edm::EDAnalyzer {
00089
00090 public:
00091 explicit IsolatedGenParticles(const edm::ParameterSet&);
00092 ~IsolatedGenParticles();
00093
00094 private:
00095 virtual void beginJob() ;
00096 virtual void analyze(const edm::Event&, const edm::EventSetup&);
00097 virtual void endJob() ;
00098
00099 double DeltaPhi(double v1, double v2);
00100 double DeltaR(double eta1, double phi1, double eta2, double phi2);
00101 double DeltaR2(double eta1, double phi1, double eta2, double phi2);
00102
00103 void fillTrack (GlobalPoint& posVec, math::XYZTLorentzVector& momVec, GlobalPoint& posECAL, int pdgId, bool okECAL, bool accpet);
00104 void fillIsolatedTrack(math::XYZTLorentzVector& momVec, GlobalPoint& posECAL, int pdgId);
00105 void BookHistograms();
00106 void clearTreeVectors();
00107 int particleCode(int);
00108
00109
00110 static const int NPBins = 3;
00111 static const int NEtaBins = 4;
00112 static const int Particles=12;
00113 int nEventProc;
00114 double genPartPBins[NPBins+1], genPartEtaBins[NEtaBins+1];
00115 double pSeed, ptMin, etaMax, pCutIsolate;
00116 bool a_Isolation;
00117
00118 std::string genSrc_;
00119 const MagneticField *bField;
00120
00121 bool initL1, useHepMC;
00122 static const size_t nL1BitsMax=128;
00123 std::string algoBitToName[nL1BitsMax];
00124 double a_coneR, a_charIsoR, a_neutIsoR, a_mipR;
00125
00126 bool debugL1Info_;
00127 int verbosity, debugTrks_;
00128 bool printTrkHitPattern_;
00129 int myverbose_;
00130 bool useJetTrigger_;
00131 double drLeadJetVeto_, ptMinLeadJet_;
00132 edm::InputTag L1extraTauJetSource_, L1extraCenJetSource_, L1extraFwdJetSource_;
00133 edm::InputTag L1extraMuonSource_, L1extraIsoEmSource_, L1extraNonIsoEmSource_;
00134 edm::InputTag L1GTReadoutRcdSource_, L1GTObjectMapRcdSource_;
00135
00136
00137 edm::Service<TFileService> fs;
00138
00139 TH1I *h_L1AlgoNames;
00140 TH1I *h_NEventProc;
00141 TH2D *h_pEta[Particles];
00142
00143 TTree *tree;
00144
00145 std::vector<double> *t_isoTrkPAll;
00146 std::vector<double> *t_isoTrkPtAll;
00147 std::vector<double> *t_isoTrkPhiAll;
00148 std::vector<double> *t_isoTrkEtaAll;
00149 std::vector<double> *t_isoTrkPdgIdAll;
00150 std::vector<double> *t_isoTrkDEtaAll;
00151 std::vector<double> *t_isoTrkDPhiAll;
00152
00153 std::vector<double> *t_isoTrkP;
00154 std::vector<double> *t_isoTrkPt;
00155 std::vector<double> *t_isoTrkEne;
00156 std::vector<double> *t_isoTrkEta;
00157 std::vector<double> *t_isoTrkPhi;
00158 std::vector<double> *t_isoTrkEtaEC;
00159 std::vector<double> *t_isoTrkPhiEC;
00160 std::vector<double> *t_isoTrkPdgId;
00161
00162 std::vector<double> *t_maxNearP31x31;
00163 std::vector<double> *t_cHadronEne31x31, *t_cHadronEne31x31_1, *t_cHadronEne31x31_2, *t_cHadronEne31x31_3;
00164 std::vector<double> *t_nHadronEne31x31;
00165 std::vector<double> *t_photonEne31x31;
00166 std::vector<double> *t_eleEne31x31;
00167 std::vector<double> *t_muEne31x31;
00168
00169 std::vector<double> *t_maxNearP25x25;
00170 std::vector<double> *t_cHadronEne25x25, *t_cHadronEne25x25_1, *t_cHadronEne25x25_2, *t_cHadronEne25x25_3;
00171 std::vector<double> *t_nHadronEne25x25;
00172 std::vector<double> *t_photonEne25x25;
00173 std::vector<double> *t_eleEne25x25;
00174 std::vector<double> *t_muEne25x25;
00175
00176 std::vector<double> *t_maxNearP21x21;
00177 std::vector<double> *t_cHadronEne21x21, *t_cHadronEne21x21_1, *t_cHadronEne21x21_2, *t_cHadronEne21x21_3;
00178 std::vector<double> *t_nHadronEne21x21;
00179 std::vector<double> *t_photonEne21x21;
00180 std::vector<double> *t_eleEne21x21;
00181 std::vector<double> *t_muEne21x21;
00182
00183 std::vector<double> *t_maxNearP15x15;
00184 std::vector<double> *t_cHadronEne15x15, *t_cHadronEne15x15_1, *t_cHadronEne15x15_2, *t_cHadronEne15x15_3;
00185 std::vector<double> *t_nHadronEne15x15;
00186 std::vector<double> *t_photonEne15x15;
00187 std::vector<double> *t_eleEne15x15;
00188 std::vector<double> *t_muEne15x15;
00189
00190 std::vector<double> *t_maxNearP11x11;
00191 std::vector<double> *t_cHadronEne11x11, *t_cHadronEne11x11_1, *t_cHadronEne11x11_2, *t_cHadronEne11x11_3;
00192 std::vector<double> *t_nHadronEne11x11;
00193 std::vector<double> *t_photonEne11x11;
00194 std::vector<double> *t_eleEne11x11;
00195 std::vector<double> *t_muEne11x11;
00196
00197 std::vector<double> *t_maxNearP9x9;
00198 std::vector<double> *t_cHadronEne9x9, *t_cHadronEne9x9_1, *t_cHadronEne9x9_2, *t_cHadronEne9x9_3;
00199 std::vector<double> *t_nHadronEne9x9;
00200 std::vector<double> *t_photonEne9x9;
00201 std::vector<double> *t_eleEne9x9;
00202 std::vector<double> *t_muEne9x9;
00203
00204 std::vector<double> *t_maxNearP7x7;
00205 std::vector<double> *t_cHadronEne7x7, *t_cHadronEne7x7_1, *t_cHadronEne7x7_2, *t_cHadronEne7x7_3;
00206 std::vector<double> *t_nHadronEne7x7;
00207 std::vector<double> *t_photonEne7x7;
00208 std::vector<double> *t_eleEne7x7;
00209 std::vector<double> *t_muEne7x7;
00210
00211 std::vector<double> *t_maxNearP3x3;
00212 std::vector<double> *t_cHadronEne3x3, *t_cHadronEne3x3_1, *t_cHadronEne3x3_2, *t_cHadronEne3x3_3;
00213 std::vector<double> *t_nHadronEne3x3;
00214 std::vector<double> *t_photonEne3x3;
00215 std::vector<double> *t_eleEne3x3;
00216 std::vector<double> *t_muEne3x3;
00217
00218 std::vector<double> *t_maxNearP1x1;
00219 std::vector<double> *t_cHadronEne1x1, *t_cHadronEne1x1_1, *t_cHadronEne1x1_2, *t_cHadronEne1x1_3;
00220 std::vector<double> *t_nHadronEne1x1;
00221 std::vector<double> *t_photonEne1x1;
00222 std::vector<double> *t_eleEne1x1;
00223 std::vector<double> *t_muEne1x1;
00224
00225 std::vector<double> *t_maxNearPHC1x1;
00226 std::vector<double> *t_cHadronEneHC1x1, *t_cHadronEneHC1x1_1, *t_cHadronEneHC1x1_2, *t_cHadronEneHC1x1_3;
00227 std::vector<double> *t_nHadronEneHC1x1;
00228 std::vector<double> *t_photonEneHC1x1;
00229 std::vector<double> *t_eleEneHC1x1;
00230 std::vector<double> *t_muEneHC1x1;
00231
00232 std::vector<double> *t_maxNearPHC3x3;
00233 std::vector<double> *t_cHadronEneHC3x3, *t_cHadronEneHC3x3_1, *t_cHadronEneHC3x3_2, *t_cHadronEneHC3x3_3;
00234 std::vector<double> *t_nHadronEneHC3x3;
00235 std::vector<double> *t_photonEneHC3x3;
00236 std::vector<double> *t_eleEneHC3x3;
00237 std::vector<double> *t_muEneHC3x3;
00238
00239 std::vector<double> *t_maxNearPHC5x5;
00240 std::vector<double> *t_cHadronEneHC5x5, *t_cHadronEneHC5x5_1, *t_cHadronEneHC5x5_2, *t_cHadronEneHC5x5_3;
00241 std::vector<double> *t_nHadronEneHC5x5;
00242 std::vector<double> *t_photonEneHC5x5;
00243 std::vector<double> *t_eleEneHC5x5;
00244 std::vector<double> *t_muEneHC5x5;
00245
00246 std::vector<double> *t_maxNearPHC7x7;
00247 std::vector<double> *t_cHadronEneHC7x7, *t_cHadronEneHC7x7_1, *t_cHadronEneHC7x7_2, *t_cHadronEneHC7x7_3;
00248 std::vector<double> *t_nHadronEneHC7x7;
00249 std::vector<double> *t_photonEneHC7x7;
00250 std::vector<double> *t_eleEneHC7x7;
00251 std::vector<double> *t_muEneHC7x7;
00252
00253 std::vector<double> *t_maxNearPR;
00254 std::vector<double> *t_cHadronEneR, *t_cHadronEneR_1, *t_cHadronEneR_2, *t_cHadronEneR_3;
00255 std::vector<double> *t_nHadronEneR;
00256 std::vector<double> *t_photonEneR;
00257 std::vector<double> *t_eleEneR;
00258 std::vector<double> *t_muEneR;
00259
00260 std::vector<double> *t_maxNearPIsoR;
00261 std::vector<double> *t_cHadronEneIsoR, *t_cHadronEneIsoR_1, *t_cHadronEneIsoR_2, *t_cHadronEneIsoR_3;
00262 std::vector<double> *t_nHadronEneIsoR;
00263 std::vector<double> *t_photonEneIsoR;
00264 std::vector<double> *t_eleEneIsoR;
00265 std::vector<double> *t_muEneIsoR;
00266
00267 std::vector<double> *t_maxNearPHCR;
00268 std::vector<double> *t_cHadronEneHCR, *t_cHadronEneHCR_1, *t_cHadronEneHCR_2, *t_cHadronEneHCR_3;
00269 std::vector<double> *t_nHadronEneHCR;
00270 std::vector<double> *t_photonEneHCR;
00271 std::vector<double> *t_eleEneHCR;
00272 std::vector<double> *t_muEneHCR;
00273
00274 std::vector<double> *t_maxNearPIsoHCR;
00275 std::vector<double> *t_cHadronEneIsoHCR, *t_cHadronEneIsoHCR_1, *t_cHadronEneIsoHCR_2, *t_cHadronEneIsoHCR_3;
00276 std::vector<double> *t_nHadronEneIsoHCR;
00277 std::vector<double> *t_photonEneIsoHCR;
00278 std::vector<double> *t_eleEneIsoHCR;
00279 std::vector<double> *t_muEneIsoHCR;
00280
00281 std::vector<int> *t_L1Decision;
00282 std::vector<double> *t_L1CenJetPt, *t_L1CenJetEta, *t_L1CenJetPhi;
00283 std::vector<double> *t_L1FwdJetPt, *t_L1FwdJetEta, *t_L1FwdJetPhi;
00284 std::vector<double> *t_L1TauJetPt, *t_L1TauJetEta, *t_L1TauJetPhi;
00285 std::vector<double> *t_L1MuonPt, *t_L1MuonEta, *t_L1MuonPhi;
00286 std::vector<double> *t_L1IsoEMPt, *t_L1IsoEMEta, *t_L1IsoEMPhi;
00287 std::vector<double> *t_L1NonIsoEMPt, *t_L1NonIsoEMEta, *t_L1NonIsoEMPhi;
00288 std::vector<double> *t_L1METPt, *t_L1METEta, *t_L1METPhi;
00289
00290 spr::genSimInfo isoinfo1x1, isoinfo3x3, isoinfo7x7, isoinfo9x9, isoinfo11x11;
00291 spr::genSimInfo isoinfo15x15, isoinfo21x21, isoinfo25x25, isoinfo31x31;
00292 spr::genSimInfo isoinfoHC1x1, isoinfoHC3x3, isoinfoHC5x5, isoinfoHC7x7;
00293 spr::genSimInfo isoinfoR, isoinfoIsoR, isoinfoHCR, isoinfoIsoHCR;
00294
00295 };
00296
00297 #endif