00001
00002 #include "DataFormats/TrackReco/interface/TrackBase.h"
00003 #include "Calibration/IsolatedParticles/plugins/IsolatedTracksHcalScale.h"
00004
00005 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
00006 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
00007 #include "Calibration/IsolatedParticles/interface/eCone.h"
00008 #include "Calibration/IsolatedParticles/interface/eHCALMatrixExtra.h"
00009 #include "MagneticField/Engine/interface/MagneticField.h"
00010 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00011 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
00012 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
00013 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
00014 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
00015
00016 IsolatedTracksHcalScale::IsolatedTracksHcalScale(const edm::ParameterSet& iConfig) {
00017
00018
00019 doMC = iConfig.getUntrackedParameter<bool>("DoMC", false);
00020 myverbose = iConfig.getUntrackedParameter<int>("Verbosity", 5 );
00021 theTrackQuality = iConfig.getUntrackedParameter<std::string>("TrackQuality","highPurity");
00022 reco::TrackBase::TrackQuality trackQuality_=reco::TrackBase::qualityByName(theTrackQuality);
00023 selectionParameters.minPt = iConfig.getUntrackedParameter<double>("MinTrackPt", 10.0);
00024 selectionParameters.minQuality = trackQuality_;
00025 selectionParameters.maxDxyPV = iConfig.getUntrackedParameter<double>("MaxDxyPV", 0.2);
00026 selectionParameters.maxDzPV = iConfig.getUntrackedParameter<double>("MaxDzPV", 5.0);
00027 selectionParameters.maxChi2 = iConfig.getUntrackedParameter<double>("MaxChi2", 5.0);
00028 selectionParameters.maxDpOverP = iConfig.getUntrackedParameter<double>("MaxDpOverP", 0.1);
00029 selectionParameters.minOuterHit = iConfig.getUntrackedParameter<int>("MinOuterHit", 4);
00030 selectionParameters.minLayerCrossed = iConfig.getUntrackedParameter<int>("MinLayerCrossed", 8);
00031 selectionParameters.maxInMiss = iConfig.getUntrackedParameter<int>("MaxInMiss", 0);
00032 selectionParameters.maxOutMiss = iConfig.getUntrackedParameter<int>("MaxOutMiss", 0);
00033 a_coneR = iConfig.getUntrackedParameter<double>("ConeRadius",34.98);
00034 a_charIsoR = a_coneR + 28.9;
00035 a_neutIsoR = a_charIsoR*0.726;
00036 a_mipR = iConfig.getUntrackedParameter<double>("ConeRadiusMIP",14.0);
00037 tMinE_ = iConfig.getUntrackedParameter<double>("TimeMinCutECAL", -500.);
00038 tMaxE_ = iConfig.getUntrackedParameter<double>("TimeMaxCutECAL", 500.);
00039
00040 if (myverbose>=0) {
00041 std::cout <<"Parameters read from config file \n"
00042 <<" doMC " << doMC
00043 <<"\t myverbose " << myverbose
00044 <<"\t minPt " << selectionParameters.minPt
00045 <<"\t theTrackQuality " << theTrackQuality
00046 <<"\t minQuality " << selectionParameters.minQuality
00047 <<"\t maxDxyPV " << selectionParameters.maxDxyPV
00048 <<"\t maxDzPV " << selectionParameters.maxDzPV
00049 <<"\t maxChi2 " << selectionParameters.maxChi2
00050 <<"\t maxDpOverP " << selectionParameters.maxDpOverP
00051 <<"\t minOuterHit " << selectionParameters.minOuterHit
00052 <<"\t minLayerCrossed " << selectionParameters.minLayerCrossed
00053 <<"\t maxInMiss " << selectionParameters.maxInMiss
00054 <<"\t maxOutMiss " << selectionParameters.maxOutMiss
00055 <<"\t a_coneR " << a_coneR
00056 <<"\t a_charIsoR " << a_charIsoR
00057 <<"\t a_neutIsoR " << a_neutIsoR
00058 <<"\t a_mipR " << a_mipR
00059 <<"\t time Range (" << tMinE_ << ":" << tMaxE_ << ")"
00060 << std::endl;
00061 }
00062 initL1 = false;
00063
00064 }
00065
00066 IsolatedTracksHcalScale::~IsolatedTracksHcalScale() {}
00067
00068 void IsolatedTracksHcalScale::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
00069
00070 edm::ESHandle<MagneticField> bFieldH;
00071 iSetup.get<IdealMagneticFieldRecord>().get(bFieldH);
00072 bField = bFieldH.product();
00073
00074
00075 edm::ESHandle<CaloGeometry> pG;
00076 iSetup.get<CaloGeometryRecord>().get(pG);
00077 const CaloGeometry* geo = pG.product();
00078
00079 edm::ESHandle<CaloTopology> theCaloTopology;
00080 iSetup.get<CaloTopologyRecord>().get(theCaloTopology);
00081 const CaloTopology *caloTopology = theCaloTopology.product();
00082
00083
00084
00085
00086
00087
00088
00089
00090 edm::ESHandle<EcalChannelStatus> ecalChStatus;
00091 iSetup.get<EcalChannelStatusRcd>().get(ecalChStatus);
00092 const EcalChannelStatus* theEcalChStatus = ecalChStatus.product();
00093
00094
00095
00096
00097
00098
00099
00100
00101 clearTreeVectors();
00102
00103 nEventProc++;
00104
00105 t_RunNo = iEvent.id().run();
00106 t_EvtNo = iEvent.id().event();
00107 t_Lumi = iEvent.luminosityBlock();
00108 t_Bunch = iEvent.bunchCrossing();
00109 if (myverbose>0) std::cout << nEventProc << " Run " << t_RunNo << " Event " << t_EvtNo << " Lumi " << t_Lumi << " Bunch " << t_Bunch << std::endl;
00110
00111 edm::Handle<reco::TrackCollection> trkCollection;
00112 iEvent.getByLabel("generalTracks", trkCollection);
00113
00114 edm::Handle<reco::VertexCollection> recVtxs;
00115 iEvent.getByLabel("offlinePrimaryVertices",recVtxs);
00116
00117
00118 edm::Handle<reco::BeamSpot> beamSpotH;
00119 iEvent.getByLabel("offlineBeamSpot", beamSpotH);
00120
00121 math::XYZPoint leadPV(0,0,0);
00122 if (recVtxs->size()>0 && !((*recVtxs)[0].isFake())) {
00123 leadPV = math::XYZPoint( (*recVtxs)[0].x(),(*recVtxs)[0].y(), (*recVtxs)[0].z() );
00124 } else if (beamSpotH.isValid()) {
00125 leadPV = beamSpotH->position();
00126 }
00127
00128 if (myverbose>0) {
00129 std::cout << "Primary Vertex " << leadPV;
00130 if (beamSpotH.isValid()) std::cout << " Beam Spot " << beamSpotH->position();
00131 std::cout << std::endl;
00132 }
00133
00134 std::vector<spr::propagatedTrackDirection> trkCaloDirections;
00135 spr::propagateCALO(trkCollection, geo, bField, theTrackQuality, trkCaloDirections, (myverbose>2));
00136 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
00137
00138 edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
00139 edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
00140 iEvent.getByLabel("ecalRecHit","EcalRecHitsEB",barrelRecHitsHandle);
00141 iEvent.getByLabel("ecalRecHit","EcalRecHitsEE",endcapRecHitsHandle);
00142
00143 edm::Handle<HBHERecHitCollection> hbhe;
00144 iEvent.getByLabel("hbhereco",hbhe);
00145 const HBHERecHitCollection Hithbhe = *(hbhe.product());
00146
00147
00148 edm::Handle<edm::SimTrackContainer> SimTk;
00149 edm::SimTrackContainer::const_iterator simTrkItr;
00150 edm::Handle<edm::SimVertexContainer> SimVtx;
00151 edm::SimVertexContainer::const_iterator vtxItr = SimVtx->begin();
00152
00153
00154 edm::Handle<edm::PCaloHitContainer> pcaloeb;
00155 edm::Handle<edm::PCaloHitContainer> pcaloee;
00156 edm::Handle<edm::PCaloHitContainer> pcalohh;
00157
00158
00159 TrackerHitAssociator* associate=0;
00160
00161 if (doMC) {
00162 iEvent.getByLabel("g4SimHits",SimTk);
00163 iEvent.getByLabel("g4SimHits",SimVtx);
00164 iEvent.getByLabel("g4SimHits", "EcalHitsEB", pcaloeb);
00165 iEvent.getByLabel("g4SimHits", "EcalHitsEE", pcaloee);
00166 iEvent.getByLabel("g4SimHits", "HcalHits", pcalohh);
00167 associate = new TrackerHitAssociator(iEvent);
00168 }
00169
00170 unsigned int nTracks=0;
00171 for (trkDetItr = trkCaloDirections.begin(),nTracks=0; trkDetItr != trkCaloDirections.end(); trkDetItr++,nTracks++){
00172 const reco::Track* pTrack = &(*(trkDetItr->trkItr));
00173 if (spr::goodTrack(pTrack,leadPV,selectionParameters,(myverbose>2)) && trkDetItr->okECAL && trkDetItr->okHCAL) {
00174 int nRH_eMipDR=0, nRH_eDR=0, nNearTRKs=0, nRecHitsCone=-99;
00175 double distFromHotCell=-99.0, distFromHotCell2=-99.0;
00176 int ietaHotCell=-99, iphiHotCell=-99;
00177 int ietaHotCell2=-99, iphiHotCell2=-99;
00178 GlobalPoint gposHotCell(0.,0.,0.), gposHotCell2(0.,0.,0.);
00179 std::vector<DetId> coneRecHitDetIds, coneRecHitDetIds2;
00180 std::pair<double, bool> e11x11_20SigP, e15x15_20SigP;
00181 double hCone = spr::eCone_hcal(geo, hbhe, trkDetItr->pointHCAL, trkDetItr->pointECAL,
00182 a_coneR, trkDetItr->directionHCAL, nRecHitsCone,
00183 coneRecHitDetIds, distFromHotCell,
00184 ietaHotCell, iphiHotCell, gposHotCell);
00185 double eHCALDR = spr::eCone_hcal(geo, hbhe, trkDetItr->pointHCAL, trkDetItr->pointECAL,
00186 a_charIsoR, trkDetItr->directionHCAL, nRecHitsCone,
00187 coneRecHitDetIds2, distFromHotCell2,
00188 ietaHotCell2, iphiHotCell2, gposHotCell2);
00189
00190 double conehmaxNearP = spr::chargeIsolationCone(nTracks, trkCaloDirections, a_charIsoR, nNearTRKs, (myverbose>3));
00191
00192 double eMipDR = spr::eCone_ecal(geo, barrelRecHitsHandle, endcapRecHitsHandle,
00193 trkDetItr->pointHCAL, trkDetItr->pointECAL,
00194 a_mipR, trkDetItr->directionECAL, nRH_eMipDR);
00195 double eECALDR = spr::eCone_ecal(geo, barrelRecHitsHandle, endcapRecHitsHandle,
00196 trkDetItr->pointHCAL, trkDetItr->pointECAL,
00197 a_neutIsoR, trkDetItr->directionECAL, nRH_eDR);
00198
00199 HcalDetId closestCell = (HcalDetId)(trkDetItr->detIdHCAL);
00200
00201 edm::ESHandle<EcalSeverityLevelAlgo> sevlv;
00202 iSetup.get<EcalSeverityLevelAlgoRcd>().get(sevlv);
00203
00204 e11x11_20SigP = spr::eECALmatrix(trkDetItr->detIdECAL,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),5,5, 0.060, 0.300, tMinE_,tMaxE_);
00205 e15x15_20SigP = spr::eECALmatrix(trkDetItr->detIdECAL,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),7,7, 0.060, 0.300, tMinE_,tMaxE_);
00206
00207
00208 t_trackP ->push_back( pTrack->p() );
00209 t_trackPt ->push_back( pTrack->pt() );
00210 t_trackEta ->push_back( pTrack->momentum().eta() );
00211 t_trackPhi ->push_back( pTrack->momentum().phi() );
00212 t_trackHcalEta ->push_back( closestCell.ieta() );
00213 t_trackHcalPhi ->push_back( closestCell.iphi() );
00214 t_hCone ->push_back( hCone);
00215 t_conehmaxNearP ->push_back( conehmaxNearP);
00216 t_eMipDR ->push_back( eMipDR);
00217 t_eECALDR ->push_back( eECALDR);
00218 t_eHCALDR ->push_back( eHCALDR);
00219 t_e11x11_20Sig ->push_back( e11x11_20SigP.first );
00220 t_e15x15_20Sig ->push_back( e15x15_20SigP.first );
00221
00222 if (myverbose > 0) {
00223 std::cout << "Track p " << pTrack->p() << " pt " << pTrack->pt()
00224 << " eta " << pTrack->momentum().eta() << " phi "
00225 << pTrack->momentum().phi() << " ieta/iphi ("
00226 << closestCell.ieta() << ", " << closestCell.iphi()
00227 << ") Energy in cone " << hCone << " Charge Isolation "
00228 << conehmaxNearP << " eMIP " << eMipDR
00229 << " Neutral isolation (ECAL) " << eECALDR-eMipDR
00230 << " (ECAL NxN) " << e15x15_20SigP.first-e11x11_20SigP.first
00231 << " (HCAL) " << eHCALDR-hCone << std::endl;
00232 }
00233
00234 if (doMC) {
00235 int nSimHits = -999;
00236 double hsim;
00237 std::map<std::string, double> hsimInfo;
00238 std::vector<int> multiplicity;
00239 hsim = spr::eCone_hcal(geo, pcalohh, trkDetItr->pointHCAL,
00240 trkDetItr->pointECAL, a_coneR,
00241 trkDetItr->directionHCAL, nSimHits);
00242 hsimInfo = spr::eHCALSimInfoCone(iEvent, pcalohh, SimTk, SimVtx,
00243 pTrack, *associate, geo,
00244 trkDetItr->pointHCAL,
00245 trkDetItr->pointECAL, a_coneR,
00246 trkDetItr->directionHCAL,
00247 multiplicity);
00248
00249 t_hsimInfoMatched ->push_back(hsimInfo["eMatched" ]);
00250 t_hsimInfoRest ->push_back(hsimInfo["eRest" ]);
00251 t_hsimInfoPhoton ->push_back(hsimInfo["eGamma" ]);
00252 t_hsimInfoNeutHad ->push_back(hsimInfo["eNeutralHad"]);
00253 t_hsimInfoCharHad ->push_back(hsimInfo["eChargedHad"]);
00254 t_hsimInfoPdgMatched->push_back(hsimInfo["pdgMatched" ]);
00255 t_hsimInfoTotal ->push_back(hsimInfo["eTotal" ]);
00256
00257 t_hsimInfoNMatched ->push_back(multiplicity.at(0));
00258 t_hsimInfoNTotal ->push_back(multiplicity.at(1));
00259 t_hsimInfoNNeutHad ->push_back(multiplicity.at(2));
00260 t_hsimInfoNCharHad ->push_back(multiplicity.at(3));
00261 t_hsimInfoNPhoton ->push_back(multiplicity.at(4));
00262 t_hsimInfoNRest ->push_back(multiplicity.at(5));
00263
00264 t_hsim ->push_back(hsim );
00265 t_nSimHits ->push_back(nSimHits );
00266
00267 if (myverbose > 0) {
00268 std::cout << "Matched (E) " << hsimInfo["eMatched"] << " (N) "
00269 << multiplicity.at(0) << " Rest (E) " << hsimInfo["eRest"]
00270 << " (N) " << multiplicity.at(1) << " Gamma (E) "
00271 << hsimInfo["eGamma"] << " (N) " << multiplicity.at(2)
00272 << " Neutral Had (E) " << hsimInfo["eNeutralHad"]
00273 << " (N) " << multiplicity.at(3) << " Charged Had (E) "
00274 << hsimInfo["eChargedHad"] << " (N) " << multiplicity.at(4)
00275 << " Total (E) " << hsimInfo["eTotal"] << " (N) "
00276 << multiplicity.at(5) << " PDG " << hsimInfo["pdgMatched"]
00277 << " Total E " << hsim << " NHit " << nSimHits <<std::endl;
00278 }
00279 }
00280 }
00281 }
00282
00283
00284 if (associate) delete associate;
00285 tree->Fill();
00286 }
00287
00288 void IsolatedTracksHcalScale::beginJob() {
00289
00290 nEventProc=0;
00291
00292
00294 tree = fs->make<TTree>("tree", "tree");
00295 tree->SetAutoSave(10000);
00296
00297 tree->Branch("t_RunNo" ,&t_RunNo ,"t_RunNo/I");
00298 tree->Branch("t_Lumi" ,&t_Lumi ,"t_Lumi/I");
00299 tree->Branch("t_Bunch" ,&t_Bunch ,"t_Bunch/I");
00300
00301 t_trackP = new std::vector<double>();
00302 t_trackPt = new std::vector<double>();
00303 t_trackEta = new std::vector<double>();
00304 t_trackPhi = new std::vector<double>();
00305 t_trackHcalEta = new std::vector<double>();
00306 t_trackHcalPhi = new std::vector<double>();
00307 t_hCone = new std::vector<double>();
00308 t_conehmaxNearP = new std::vector<double>();
00309 t_eMipDR = new std::vector<double>();
00310 t_eECALDR = new std::vector<double>();
00311 t_eHCALDR = new std::vector<double>();
00312 t_e11x11_20Sig = new std::vector<double>();
00313 t_e15x15_20Sig = new std::vector<double>();
00314
00315 tree->Branch("t_trackP", "vector<double>", &t_trackP );
00316 tree->Branch("t_trackPt", "vector<double>", &t_trackPt );
00317 tree->Branch("t_trackEta", "vector<double>", &t_trackEta );
00318 tree->Branch("t_trackPhi", "vector<double>", &t_trackPhi );
00319 tree->Branch("t_trackHcalEta", "vector<double>", &t_trackHcalEta );
00320 tree->Branch("t_trackHcalPhi", "vector<double>", &t_trackHcalPhi );
00321 tree->Branch("t_hCone", "vector<double>", &t_hCone );
00322 tree->Branch("t_conehmaxNearP", "vector<double>", &t_conehmaxNearP );
00323 tree->Branch("t_eMipDR", "vector<double>", &t_eMipDR );
00324 tree->Branch("t_eECALDR", "vector<double>", &t_eECALDR );
00325 tree->Branch("t_eHCALDR", "vector<double>", &t_eHCALDR );
00326 tree->Branch("t_e11x11_20Sig", "vector<double>", &t_e11x11_20Sig );
00327 tree->Branch("t_e15x15_20Sig", "vector<double>", &t_e15x15_20Sig );
00328
00329 if (doMC) {
00330 t_hsimInfoMatched = new std::vector<double>();
00331 t_hsimInfoRest = new std::vector<double>();
00332 t_hsimInfoPhoton = new std::vector<double>();
00333 t_hsimInfoNeutHad = new std::vector<double>();
00334 t_hsimInfoCharHad = new std::vector<double>();
00335 t_hsimInfoPdgMatched = new std::vector<double>();
00336 t_hsimInfoTotal = new std::vector<double>();
00337 t_hsimInfoNMatched = new std::vector<int>();
00338 t_hsimInfoNTotal = new std::vector<int>();
00339 t_hsimInfoNNeutHad = new std::vector<int>();
00340 t_hsimInfoNCharHad = new std::vector<int>();
00341 t_hsimInfoNPhoton = new std::vector<int>();
00342 t_hsimInfoNRest = new std::vector<int>();
00343 t_hsim = new std::vector<double>();
00344 t_nSimHits = new std::vector<int>();
00345
00346 tree->Branch("t_hsimInfoMatched", "vector<double>", &t_hsimInfoMatched );
00347 tree->Branch("t_hsimInfoRest", "vector<double>", &t_hsimInfoRest );
00348 tree->Branch("t_hsimInfoPhoton", "vector<double>", &t_hsimInfoPhoton );
00349 tree->Branch("t_hsimInfoNeutHad", "vector<double>", &t_hsimInfoNeutHad );
00350 tree->Branch("t_hsimInfoCharHad", "vector<double>", &t_hsimInfoCharHad );
00351 tree->Branch("t_hsimInfoPdgMatched", "vector<double>", &t_hsimInfoPdgMatched );
00352 tree->Branch("t_hsimInfoTotal", "vector<double>", &t_hsimInfoTotal );
00353 tree->Branch("t_hsimInfoNMatched", "vector<int>", &t_hsimInfoNMatched );
00354 tree->Branch("t_hsimInfoNTotal", "vector<int>", &t_hsimInfoNTotal );
00355 tree->Branch("t_hsimInfoNNeutHad", "vector<int>", &t_hsimInfoNNeutHad );
00356 tree->Branch("t_hsimInfoNCharHad", "vector<int>", &t_hsimInfoNCharHad );
00357 tree->Branch("t_hsimInfoNPhoton", "vector<int>", &t_hsimInfoNPhoton );
00358 tree->Branch("t_hsimInfoNRest", "vector<int>", &t_hsimInfoNRest );
00359 tree->Branch("t_hsim", "vector<double>", &t_hsim );
00360 tree->Branch("t_nSimHits", "vector<int>", &t_nSimHits );
00361 }
00362 }
00363
00364 void IsolatedTracksHcalScale::endJob() {
00365
00366 std::cout << "Number of Events Processed " << nEventProc << std::endl;
00367 }
00368
00369 void IsolatedTracksHcalScale::clearTreeVectors() {
00370
00371 t_trackP ->clear();
00372 t_trackPt ->clear();
00373 t_trackEta ->clear();
00374 t_trackPhi ->clear();
00375 t_trackHcalEta ->clear();
00376 t_trackHcalPhi ->clear();
00377 t_hCone ->clear();
00378 t_conehmaxNearP ->clear();
00379 t_eMipDR ->clear();
00380 t_eECALDR ->clear();
00381 t_eHCALDR ->clear();
00382 t_e11x11_20Sig ->clear();
00383 t_e15x15_20Sig ->clear();
00384
00385 if (doMC) {
00386 t_hsimInfoMatched ->clear();
00387 t_hsimInfoRest ->clear();
00388 t_hsimInfoPhoton ->clear();
00389 t_hsimInfoNeutHad ->clear();
00390 t_hsimInfoCharHad ->clear();
00391 t_hsimInfoPdgMatched ->clear();
00392 t_hsimInfoTotal ->clear();
00393 t_hsimInfoNMatched ->clear();
00394 t_hsimInfoNTotal ->clear();
00395 t_hsimInfoNNeutHad ->clear();
00396 t_hsimInfoNCharHad ->clear();
00397 t_hsimInfoNPhoton ->clear();
00398 t_hsimInfoNRest ->clear();
00399 t_hsim ->clear();
00400 t_nSimHits ->clear();
00401 }
00402 }
00403
00404
00405
00406 DEFINE_FWK_MODULE(IsolatedTracksHcalScale);
1.7.3