18 trackerHitAssociatorConfig_(consumesCollector()) {
55 std::cout <<
"Parameters read from config file \n"
57 <<
"\t myverbose " << myverbose
59 <<
"\t theTrackQuality " << theTrackQuality
69 <<
"\t a_coneR " << a_coneR
123 if (recVtxs->size()>0 && !((*recVtxs)[0].isFake())) {
124 leadPV =
math::XYZPoint( (*recVtxs)[0].
x(),(*recVtxs)[0].
y(), (*recVtxs)[0].
z() );
125 }
else if (beamSpotH.
isValid()) {
126 leadPV = beamSpotH->position();
130 std::cout <<
"Primary Vertex " << leadPV;
131 if (beamSpotH.
isValid())
std::cout <<
" Beam Spot " << beamSpotH->position();
135 std::vector<spr::propagatedTrackDirection> trkCaloDirections;
137 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
150 edm::SimTrackContainer::const_iterator simTrkItr;
159 std::unique_ptr<TrackerHitAssociator> associate;
170 unsigned int nTracks=0;
171 for (trkDetItr = trkCaloDirections.begin(),nTracks=0; trkDetItr != trkCaloDirections.end(); trkDetItr++,nTracks++){
172 const reco::Track* pTrack = &(*(trkDetItr->trkItr));
174 int nRH_eMipDR=0, nRH_eDR=0, nNearTRKs=0, nRecHitsCone=-99;
175 double distFromHotCell=-99.0, distFromHotCell2=-99.0;
176 int ietaHotCell=-99, iphiHotCell=-99;
177 int ietaHotCell2=-99, iphiHotCell2=-99;
178 GlobalPoint gposHotCell(0.,0.,0.), gposHotCell2(0.,0.,0.);
179 std::vector<DetId> coneRecHitDetIds, coneRecHitDetIds2;
180 std::pair<double, bool> e11x11_20SigP, e15x15_20SigP;
182 trkDetItr->pointECAL,
183 a_coneR, trkDetItr->directionHCAL,
184 nRecHitsCone, coneRecHitDetIds,
185 distFromHotCell, ietaHotCell, iphiHotCell,
188 trkDetItr->pointECAL,
189 a_coneR, trkDetItr->directionHCAL,
190 nRecHitsCone, coneRecHitDetIds,
191 distFromHotCell, ietaHotCell,
192 iphiHotCell, gposHotCell,
196 trkDetItr->directionHCAL, nRecHitsCone,
197 coneRecHitDetIds2, distFromHotCell2,
198 ietaHotCell2, iphiHotCell2, gposHotCell2,
202 trkDetItr->directionHCAL, nRecHitsCone,
203 coneRecHitDetIds2, distFromHotCell2,
204 ietaHotCell2, iphiHotCell2,
210 endcapRecHitsHandle,trkDetItr->pointHCAL,
211 trkDetItr->pointECAL,
a_mipR,
212 trkDetItr->directionECAL, nRH_eMipDR);
214 endcapRecHitsHandle,trkDetItr->pointHCAL,
216 trkDetItr->directionECAL, nRH_eDR);
218 endcapRecHitsHandle,trkDetItr->pointHCAL,
219 trkDetItr->pointECAL,
a_mipR,
220 trkDetItr->directionECAL, nRH_eMipDR,
223 endcapRecHitsHandle,trkDetItr->pointHCAL,
225 trkDetItr->directionECAL, nRH_eDR,
228 endcapRecHitsHandle,trkDetItr->pointHCAL,
229 trkDetItr->pointECAL,
a_mipR,
230 trkDetItr->directionECAL, nRH_eMipDR,
233 endcapRecHitsHandle,trkDetItr->pointHCAL,
235 trkDetItr->directionECAL, nRH_eDR,
243 e11x11_20SigP =
spr::eECALmatrix(trkDetItr->detIdECAL,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.
product(),5,5, 0.060, 0.300,
tMinE_,
tMaxE_);
244 e15x15_20SigP =
spr::eECALmatrix(trkDetItr->detIdECAL,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.
product(),7,7, 0.060, 0.300,
tMinE_,
tMaxE_);
268 std::cout <<
"Track p " << pTrack->
p() <<
" pt " << pTrack->
pt()
269 <<
" eta " << pTrack->
momentum().eta() <<
" phi "
270 << pTrack->
momentum().phi() <<
" ieta/iphi ("
271 << closestCell.
ieta() <<
", " << closestCell.
iphi()
272 <<
") Energy in cone " << hCone <<
" Charge Isolation "
273 << conehmaxNearP <<
" eMIP (" << eMipDR <<
", "
274 << eMipDR_1 <<
", " << eMipDR_2 <<
")"
275 <<
" Neutral isolation (ECAL) (" << eECALDR-eMipDR <<
", "
276 << eECALDR_1-eMipDR_1 <<
", " << eECALDR_2-eMipDR_2 <<
")"
277 <<
" (ECAL NxN) " << e15x15_20SigP.first-e11x11_20SigP.first
278 <<
" (HCAL) " << eHCALDR-hCone << std::endl;
284 std::map<std::string, double> hsimInfo;
285 std::vector<int> multiplicity;
288 trkDetItr->directionHCAL, nSimHits);
289 hsimInfo = spr::eHCALSimInfoCone(iEvent, pcalohh, SimTk, SimVtx,
290 pTrack, *associate, geo,
291 trkDetItr->pointHCAL,
293 trkDetItr->directionHCAL,
311 t_hsim ->push_back(hsim );
315 std::cout <<
"Matched (E) " << hsimInfo[
"eMatched"] <<
" (N) "
316 << multiplicity.at(0) <<
" Rest (E) " << hsimInfo[
"eRest"]
317 <<
" (N) " << multiplicity.at(1) <<
" Gamma (E) "
318 << hsimInfo[
"eGamma"] <<
" (N) " << multiplicity.at(2)
319 <<
" Neutral Had (E) " << hsimInfo[
"eNeutralHad"]
320 <<
" (N) " << multiplicity.at(3) <<
" Charged Had (E) "
321 << hsimInfo[
"eChargedHad"] <<
" (N) " << multiplicity.at(4)
322 <<
" Total (E) " << hsimInfo[
"eTotal"] <<
" (N) "
323 << multiplicity.at(5) <<
" PDG " << hsimInfo[
"pdgMatched"]
324 <<
" Total E " << hsim <<
" NHit " << nSimHits <<std::endl;
340 tree->SetAutoSave(10000);
346 t_trackP =
new std::vector<double>();
352 t_hCone =
new std::vector<double>();
354 t_eMipDR =
new std::vector<double>();
367 tree->Branch(
"t_trackPt",
"vector<double>", &t_trackPt );
368 tree->Branch(
"t_trackEta",
"vector<double>", &t_trackEta );
369 tree->Branch(
"t_trackPhi",
"vector<double>", &t_trackPhi );
370 tree->Branch(
"t_trackHcalEta",
"vector<double>", &t_trackHcalEta );
371 tree->Branch(
"t_trackHcalPhi",
"vector<double>", &t_trackHcalPhi );
372 tree->Branch(
"t_hCone",
"vector<double>", &t_hCone );
373 tree->Branch(
"t_conehmaxNearP",
"vector<double>", &t_conehmaxNearP );
374 tree->Branch(
"t_eMipDR",
"vector<double>", &t_eMipDR );
375 tree->Branch(
"t_eECALDR",
"vector<double>", &t_eECALDR );
376 tree->Branch(
"t_eHCALDR",
"vector<double>", &t_eHCALDR );
377 tree->Branch(
"t_e11x11_20Sig",
"vector<double>", &t_e11x11_20Sig );
378 tree->Branch(
"t_e15x15_20Sig",
"vector<double>", &t_e15x15_20Sig );
379 tree->Branch(
"t_eMipDR_1",
"vector<double>", &t_eMipDR_1 );
380 tree->Branch(
"t_eECALDR_1",
"vector<double>", &t_eECALDR_1 );
381 tree->Branch(
"t_eMipDR_2",
"vector<double>", &t_eMipDR_2 );
382 tree->Branch(
"t_eECALDR_2",
"vector<double>", &t_eECALDR_2 );
383 tree->Branch(
"t_hConeHB",
"vector<double>", &t_hConeHB );
384 tree->Branch(
"t_eHCALDRHB",
"vector<double>", &t_eHCALDRHB );
400 t_hsim =
new std::vector<double>();
404 tree->Branch(
"t_hsimInfoRest",
"vector<double>", &t_hsimInfoRest );
405 tree->Branch(
"t_hsimInfoPhoton",
"vector<double>", &t_hsimInfoPhoton );
406 tree->Branch(
"t_hsimInfoNeutHad",
"vector<double>", &t_hsimInfoNeutHad );
407 tree->Branch(
"t_hsimInfoCharHad",
"vector<double>", &t_hsimInfoCharHad );
408 tree->Branch(
"t_hsimInfoPdgMatched",
"vector<double>", &t_hsimInfoPdgMatched );
409 tree->Branch(
"t_hsimInfoTotal",
"vector<double>", &t_hsimInfoTotal );
410 tree->Branch(
"t_hsimInfoNMatched",
"vector<int>", &t_hsimInfoNMatched );
411 tree->Branch(
"t_hsimInfoNTotal",
"vector<int>", &t_hsimInfoNTotal );
412 tree->Branch(
"t_hsimInfoNNeutHad",
"vector<int>", &t_hsimInfoNNeutHad );
413 tree->Branch(
"t_hsimInfoNCharHad",
"vector<int>", &t_hsimInfoNCharHad );
414 tree->Branch(
"t_hsimInfoNPhoton",
"vector<int>", &t_hsimInfoNPhoton );
415 tree->Branch(
"t_hsimInfoNRest",
"vector<int>", &t_hsimInfoNRest );
416 tree->Branch(
"t_hsim",
"vector<double>", &t_hsim );
417 tree->Branch(
"t_nSimHits",
"vector<int>", &t_nSimHits );
double p() const
momentum vector magnitude
std::vector< double > * t_hsimInfoMatched
edm::Service< TFileService > fs
EventNumber_t event() const
T getUntrackedParameter(std::string const &, T const &) const
std::vector< double > * t_eHCALDR
std::vector< int > * t_hsimInfoNRest
std::string theTrackQuality
std::vector< double > * t_hCone
std::vector< spr::propagatedTrackID > propagateCALO(edm::Handle< reco::TrackCollection > &trkCollection, const CaloGeometry *geo, const MagneticField *bField, std::string &theTrackQuality, bool debug=false)
IsolatedTracksHcalScale(const edm::ParameterSet &)
edm::EDGetTokenT< HBHERecHitCollection > tok_hbhe_
std::vector< int > * t_hsimInfoNTotal
std::vector< int > * t_hsimInfoNCharHad
edm::EDGetTokenT< edm::PCaloHitContainer > tok_caloEE_
TrackerHitAssociator::Config trackerHitAssociatorConfig_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
std::vector< int > * t_hsimInfoNMatched
TrackQuality
track quality
#define DEFINE_FWK_MODULE(type)
edm::EDGetTokenT< EcalRecHitCollection > tok_EE_
std::vector< double > * t_hsim
int bunchCrossing() const
std::vector< double > * t_hsimInfoNeutHad
edm::LuminosityBlockNumber_t luminosityBlock() const
T * make(const Args &...args) const
make new ROOT object
std::vector< double > * t_eMipDR_2
std::vector< int > * t_hsimInfoNPhoton
double chargeIsolationCone(unsigned int trkIndex, std::vector< spr::propagatedTrackDirection > &trkDirs, double dR, int &nNearTRKs, bool debug=false)
const Vector & momentum() const
track momentum vector
~IsolatedTracksHcalScale()
edm::EDGetTokenT< EcalRecHitCollection > tok_EB_
std::vector< double > * t_e11x11_20Sig
std::vector< int > * t_hsimInfoNNeutHad
bool goodTrack(const reco::Track *pTrack, math::XYZPoint leadPV, trackSelectionParameters parameters, bool debug=false)
edm::EDGetTokenT< edm::PCaloHitContainer > tok_caloEB_
double pt() const
track transverse momentum
void analyze(const edm::Event &, const edm::EventSetup &)
edm::EDGetTokenT< reco::VertexCollection > tok_recVtx_
std::vector< double > * t_trackEta
int ieta() const
get the cell ieta
std::vector< double > * t_eECALDR_2
edm::EDGetTokenT< edm::SimTrackContainer > tok_simTk_
std::vector< double > * t_eECALDR
spr::trackSelectionParameters selectionParameters
double eCone_ecal(const CaloGeometry *geo, edm::Handle< T > &barrelhits, edm::Handle< T > &endcaphits, const GlobalPoint &hpoint1, const GlobalPoint &point1, double dR, const GlobalVector &trackMom, int &nRecHits, double ebThr=-100, double eeThr=-100, double tMin=-500, double tMax=500, bool debug=false)
std::vector< double > * t_eMipDR_1
int iphi() const
get the cell iphi
std::vector< double > * t_trackPt
static TrackQuality qualityByName(const std::string &name)
const MagneticField * bField
T const * product() const
std::vector< double > * t_trackHcalEta
std::vector< double > * t_trackP
double eCone_hcal(const CaloGeometry *geo, edm::Handle< T > &hits, const GlobalPoint &hpoint1, const GlobalPoint &point1, double dR, const GlobalVector &trackMom, int &nRecHits, double hbThr=-100, double heThr=-100, double hfThr=-100, double hoThr=-100, double tMin=-500, double tMax=500, int detOnly=-1, bool debug=false)
std::vector< double > * t_eECALDR_1
XYZPointD XYZPoint
point in space with cartesian internal representation
std::vector< double > * t_trackHcalPhi
T const * product() const
std::vector< double > * t_hsimInfoPhoton
std::vector< double > * t_hsimInfoTotal
edm::EDGetTokenT< edm::PCaloHitContainer > tok_caloHH_
std::vector< double > * t_hsimInfoCharHad
reco::TrackBase::TrackQuality minQuality
edm::EDGetTokenT< reco::TrackCollection > tok_genTrack_
std::vector< double > * t_trackPhi
std::vector< double > * t_eHCALDRHB
std::vector< double > * t_hsimInfoPdgMatched
edm::EDGetTokenT< reco::BeamSpot > tok_bs_
edm::EDGetTokenT< edm::SimVertexContainer > tok_simVtx_
std::vector< int > * t_nSimHits
std::vector< double > * t_hConeHB
std::vector< double > * t_hsimInfoRest
std::vector< double > * t_conehmaxNearP
std::vector< double > * t_e15x15_20Sig
std::vector< double > * t_eMipDR
double eECALmatrix(const DetId &detId, edm::Handle< T > &hitsEB, edm::Handle< T > &hitsEE, const CaloGeometry *geo, const CaloTopology *caloTopology, int ieta, int iphi, double ebThr=-100, double eeThr=-100, double tMin=-500, double tMax=500, bool debug=false)