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