46 virtual void endJob() ;
93 Int_t numValidTrkHits[50], numValidTrkStrips[50], numLayers[50];
119 parameters_.loadParameters( parameters );
120 trackAssociator_.useDefaultPropagator();
122 associationConeSize_=iConfig.
getParameter<
double>(
"associationConeSize");
123 clusterConeSize_=iConfig.
getParameter<
double>(
"clusterConeSize");
125 trackIsolationCone_ = iConfig.
getParameter<
double>(
"trackIsolationCone");
126 neutralIsolationCone_ = iConfig.
getParameter<
double>(
"neutralIsolationCone");
135 Float_t resprecal = 1.;
136 Float_t pfrecal = 1.;
138 if(Respcorr_) resprecal = respRecalib -> getValues(
id)->getValue();
139 if(PFcorr_) pfrecal = pfRecalib -> getValues(
id)->getValue();
141 Float_t factor = resprecal*pfrecal;
155 respRecalib = recalibCorrs.
product();
165 LogWarning(
"CalibConstants")<<
" Not Found!! ";
181 ev.
getByLabel(
"ecalRecHit",
"EcalRecHitsEE",ecalEE);
185 ev.
getByLabel(
"ecalRecHit",
"EcalRecHitsEB",ecalEB);
191 {tmpEcalRecHitCollection->push_back(*recHit);}
193 {tmpEcalRecHitCollection->push_back(*recHit);}
198 ev.
getByLabel(
"generalTracks", generalTracks);
212 parameters_.useEcal =
true;
213 parameters_.useHcal =
true;
214 parameters_.useCalo =
false;
215 parameters_.useMuon =
false;
216 parameters_.dREcal = 0.5;
217 parameters_.dRHcal = 0.6;
226 stepPropF->setMaterialMode(
false);
227 stepPropF->applyRadX0Correction(
true);
248 std::cout <<
"no HepMCProduct found" << std::endl;
257 double maxPt = -99999.;
262 HepMC::GenEvent * myGenEvent =
new HepMC::GenEvent(*(evtMC->GetEvent()));
263 for ( HepMC::GenEvent::particle_iterator
p = myGenEvent->particles_begin();
264 p != myGenEvent->particles_end(); ++
p )
266 phiParticle = (*p)->momentum().phi();
267 etaParticle = (*p)->momentum().eta();
268 double pt = (*p)->momentum().perp();
269 mom_MC = (*p)->momentum().rho();
270 if(pt > maxPt) { npart++; maxPt = pt; }
271 GlobalVector mom ((*p)->momentum().x(),(*p)->momentum().y(),(*p)->momentum().z());
274 if(
abs((*p)->pdg_id())==211) charge = (*p)->pdg_id()/
abs((*p)->pdg_id());
287 GlobalPoint barrelMC(0,0,0), endcapMC(0,0,0), forwardMC1(0,0,0), forwardMC2(0,0,0);
303 if(fabs(etaParticle)<1.392) {
308 if(steppingHelixstateinfo_.
isValid() )
315 if(fabs(etaParticle)>=1.392 && fabs(etaParticle)<5.191) {
331 if (
abs(etaParticle)>5.191)
continue;
333 if(
abs(etaParticle)>2.99) {
344 if(steppingHelixstateinfo_.
isValid() )
347 steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*aPlane2));
348 if(steppingHelixstateinfo_.
isValid() )
357 Int_t iphitrue = -10;
358 Int_t ietatrue = 100;
362 if (
abs(etaParticle)<1.392)
364 gPointHcal = barrelMC;
367 if (
abs(etaParticle)>=1.392 &&
abs(etaParticle)<3.0)
369 gPointHcal = endcapMC;
372 if (
abs(etaParticle)>=3.0 &&
abs(etaParticle)<5.191)
382 gPointHcal = forwardMC1;
391 ietatrue = tempId.
ieta();
392 iphitrue = tempId.
iphi();
394 etaGPoint = gPointHcal.
eta();
395 phiGPoint = gPointHcal.
phi();
403 if (gPointHcal.
x()==0 && gPointHcal.
y()==0 && gPointHcal.
z()==0)
407 float etahcal=gPointHcal.
eta();
409 if (
abs(etahcal)>5.192)
continue;
427 GlobalPoint gPointEcal(xTrkEcal,yTrkEcal,zTrkEcal);
457 recal = RecalibFactor(hhit->detid());
460 int iphihit = (hhit->id()).iphi();
461 int ietahit = (hhit->id()).ieta();
462 int depthhit = (hhit->id()).depth();
463 float enehit = hhit->energy()* recal;
465 if (depthhit!=1)
continue;
469 if(distAtHcal < clusterConeSize_)
474 int iphihit2 = (hhit2->id()).iphi();
475 int ietahit2 = (hhit2->id()).ieta();
476 int depthhit2 = (hhit2->id()).depth();
477 float enehit2 = hhit2->energy() * recal;
479 if (iphihit==iphihit2 && ietahit==ietahit2 && depthhit!=depthhit2) enehit = enehit+enehit2;
487 MaxHit.
ietahitm = (hhit->id()).ieta();
488 MaxHit.
iphihitm = (hhit->id()).iphi();
489 MaxHit.
dr = distAtHcal;
500 recal = RecalibFactor(hhit->detid());
504 int iphihit = (hhit->id()).iphi();
505 int ietahit = (hhit->id()).ieta();
506 int depthhit = (hhit->id()).depth();
507 float enehit = hhit->energy()* recal;
511 if(distAtHcal < associationConeSize_)
515 int iphihit2 = (hhit2->id()).iphi();
516 int ietahit2 = (hhit2->id()).ieta();
517 int depthhit2 = (hhit2->id()).depth();
518 float enehit2 = hhit2->energy() * recal;
520 if (iphihit==iphihit2 && ietahit==ietahit2 && depthhit!=depthhit2) enehit = enehit+enehit2;
527 MaxHit.
ietahitm = (hhit->id()).ieta();
528 MaxHit.
iphihitm = (hhit->id()).iphi();
529 MaxHit.
dr = distAtHcal;
544 recal = RecalibFactor(hhit->detid());
549 int iphihit = (hhit->id()).iphi();
550 int ietahit = (hhit->id()).ieta();
551 int depthhit = (hhit->id()).depth();
552 float enehit = hhit->energy()* recal;
557 int iphihitNoise = iphihit >36 ? iphihit-36 : iphihit+36;
558 int ietahitNoise = -ietahit;
559 int depthhitNoise = depthhit;
562 if(distAtHcal < clusterConeSize_ && MaxHit.hitenergy > 0.)
569 if(MaxHit.
ietahitm*(hhit->id()).ieta()>0)
570 { DIETA = MaxHit.
ietahitm - (hhit->id()).ieta();}
571 if(MaxHit.
ietahitm*(hhit->id()).ieta()<0)
572 { DIETA = MaxHit.
ietahitm - (hhit->id()).ieta(); DIETA = DIETA>0 ? DIETA-1 : DIETA+1;}
573 int DIPHI =
abs(MaxHit.
iphihitm - (hhit->id()).iphi());
574 DIPHI = DIPHI>36 ? 72-DIPHI : DIPHI;
578 {e5x5 += hhit->energy();}
580 {e3x3 += hhit->energy();}
586 recal = RecalibFactor(hhit2->detid());
587 int iphihit2 = (hhit2->id()).iphi();
588 int ietahit2 = (hhit2->id()).ieta();
589 int depthhit2 = (hhit2->id()).depth();
590 float enehit2 = hhit2->energy()* recal;
592 if (iphihitNoise == iphihit2 && ietahitNoise == ietahit2 && depthhitNoise == depthhit2 && enehit2>0.)
594 eHcalConeNoise += hhit2->energy()*recal;
606 recal = RecalibFactor(hhit->detid());
612 int iphihit = (hhit->id()).iphi();
613 int ietahit = (hhit->id()).ieta();
614 int depthhit = (hhit->id()).depth();
615 float enehit = hhit->energy()* recal;
618 int iphihitNoise = iphihit >36 ? iphihit-36 : iphihit+36;
619 int ietahitNoise = -ietahit;
620 int depthhitNoise = depthhit;
624 if(distAtHcal < clusterConeSize_ && MaxHit.hitenergy > 0.)
632 if(MaxHit.
ietahitm*(hhit->id()).ieta()>0)
633 { DIETA = MaxHit.
ietahitm - (hhit->id()).ieta();}
634 if(MaxHit.
ietahitm*(hhit->id()).ieta()<0)
635 { DIETA = MaxHit.
ietahitm - (hhit->id()).ieta(); DIETA = DIETA>0 ? DIETA-1 : DIETA+1;}
636 int DIPHI =
abs(MaxHit.
iphihitm - (hhit->id()).iphi());
637 DIPHI = DIPHI>36 ? 72-DIPHI : DIPHI;
641 {e5x5 += hhit->energy();}
643 {e3x3 += hhit->energy();}
647 recal = RecalibFactor(hhit2->detid());
649 int iphihit2 = (hhit2->id()).iphi();
650 int ietahit2 = (hhit2->id()).ieta();
651 int depthhit2 = (hhit2->id()).depth();
652 float enehit2 = hhit2->energy()* recal;
654 if (iphihitNoise == iphihit2 && ietahitNoise == ietahit2 && depthhitNoise == depthhit2 && enehit2>0.01)
656 eHcalConeNoise += hhit2->energy()*recal;
673 float delR_track_particle = 100;
675 for (reco::TrackCollection::const_iterator track1=generalTracks->begin(); track1!=generalTracks->end(); track1++)
677 delR_track_particle =
deltaR(etaParticle, phiParticle, track1->eta(), track1->phi());
681 trackP[nTracks] =
sqrt(track1->px()*track1->px() + track1->py()*track1->py() + track1->pz()*track1->pz());
683 delRmc[nTracks] = delR_track_particle;
684 numValidTrkHits[nTracks] = track1->hitPattern().numberOfValidHits();
685 numValidTrkStrips[nTracks] = track1->hitPattern().numberOfValidStripTECHits();
686 numLayers[nTracks] = track1->hitPattern().trackerLayersWithMeasurement();
702 diphi_M_P = diphi_M_P>36 ? 72-diphi_M_P : diphi_M_P;
703 iDr =
sqrt(diphi_M_P*diphi_M_P+dieta_M_P*dieta_M_P);
713 Bool_t passCuts = kFALSE;
754 pfTree = fs -> make<TTree>(
"pfTree",
"Tree for pf info");
757 pfTree->Branch(
"nTracks", &nTracks,
"nTracks/I");
758 pfTree->Branch(
"trackEta",
trackEta,
"trackEta[nTracks]/F");
759 pfTree->Branch(
"trackPhi",
trackPhi,
"trackPhi[nTracks]/F");
760 pfTree->Branch(
"trackP", trackP,
"trackP[nTracks]/F");
762 pfTree->Branch(
"delRmc", delRmc,
"delRmc[nTracks]/F");
763 pfTree->Branch(
"numValidTrkHits", numValidTrkHits,
"numValidTrkHits[nTracks]/I");
764 pfTree->Branch(
"numValidTrkStrips", numValidTrkStrips,
"numValidTrkStrips[nTracks]/I");
765 pfTree->Branch(
"numLayers", numLayers,
"numLayers[nTracks]/I");
766 pfTree->Branch(
"trkQual", trkQual,
"trkQual[nTracks]/O");
769 pfTree->Branch(
"eEcalCone", &eEcalCone,
"eEcalCone/F");
770 pfTree->Branch(
"eHcalCone", &eHcalCone,
"eHcalCone/F");
771 pfTree->Branch(
"eHcalConeNoise", &eHcalConeNoise,
"eHcalConeNoise/F");
773 pfTree->Branch(
"UsedCellsNoise", &UsedCellsNoise,
"UsedCellsNoise/I");
774 pfTree->Branch(
"UsedCells", &UsedCells,
"UsedCells/I");
776 pfTree->Branch(
"eCentHit", &eCentHit ,
"eCentHit/F");
778 pfTree->Branch(
"eParticle", &eParticle,
"eParticle/F");
779 pfTree->Branch(
"etaParticle", &etaParticle,
"etaParticle/F");
780 pfTree->Branch(
"phiParticle", &phiParticle,
"phiParticle/F");
782 pfTree->Branch(
"etaGPoint", &etaGPoint,
"etaGPoint/F");
783 pfTree->Branch(
"phiGPoint", &phiGPoint,
"phiGPoint/F");
785 pfTree->Branch(
"xAtHcal", &xAtHcal,
"xAtHcal/F");
786 pfTree->Branch(
"yAtHcal", &yAtHcal,
"yAtHcal/F");
787 pfTree->Branch(
"zAtHcal", &zAtHcal,
"zAtHcal/F");
789 pfTree->Branch(
"eECAL09cm", &eECAL09cm,
"eECAL09cm/F");
790 pfTree->Branch(
"eECAL40cm", &eECAL40cm,
"eECAL40cm/F");
791 pfTree->Branch(
"eECAL", &eECAL,
"eECAL/F");
793 pfTree->Branch(
"e3x3 ", &e3x3 ,
"e3x3/F");
794 pfTree->Branch(
"e5x5", &e5x5 ,
"e5x5/F");
796 pfTree->Branch(
"iDr", &iDr,
"iDr/F");
797 pfTree->Branch(
"delR", &delR,
"delR/F");
799 pfTree->Branch(
"iEta", &iEta,
"iEta/I");
800 pfTree->Branch(
"iPhi", &iPhi,
"iPhi/I");
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
const HcalGeometry * geoHcal
const HcalPFCorrs * pfRecalib
double trackIsolationCone_
double ecalEnergyInCone(const GlobalPoint center, double radius, const EcalRecHitCollection ecalCol, const CaloGeometry *geo)
#define DEFINE_FWK_MODULE(type)
Geom::Phi< T > phi() const
ReturnType plane(const PositionType &pos, const RotationType &rot) const
TrackAssociatorParameters parameters_
std::vector< EcalRecHit >::const_iterator const_iterator
bool getByType(Handle< PROD > &result) const
tuple SteppingHelixPropagator
FreeTrajectoryState * freeState(bool withErrors=true) const
int ieta() const
get the cell ieta
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
const HcalRespCorrs * respRecalib
const_iterator end() const
double deltaR(double eta1, double eta2, double phi1, double phi2)
virtual DetId getClosestCell(const GlobalPoint &r) const
int iphi() const
get the cell iphi
virtual void analyze(edm::Event const &ev, edm::EventSetup const &c)
GlobalPoint position() const
SteppingHelixPropagator * stepPropF
T const * product() const
T const * product() const
MagneticField * theMagField
TrackDetectorAssociator trackAssociator_
double RecalibFactor(HcalDetId id)
math::XYZPoint trkGlobPosAtEcal
Track position at different parts of the calorimeter.
double getDistInPlaneSimple(const GlobalPoint caloPoint, const GlobalPoint rechitPoint)
edm::Service< TFileService > fs
const_iterator begin() const
Global3DVector GlobalVector
HcalCorrPFCalculation(edm::ParameterSet const &conf)