7 #include "CLHEP/Random/RandomEngine.h"
22 double Rcurv=pT*33.3*100/38;
27 if (theta>0.5*acos(-1)) theta=acos(-1)-
theta;
28 if (fabs(etaIP)<1.479) {
29 deltaPhi=-charge*asin(0.5*ecRad/Rcurv);
30 double alpha1=2*asin(0.5*ecRad/Rcurv);
31 double z=ecRad/
tan(theta);
32 if (etaIP>0) zNew=z*(Rcurv*alpha1)/ecRad+vtxZ;
33 else zNew=-z*(Rcurv*alpha1)/ecRad+vtxZ;
34 double zAbs=fabs(zNew);
36 etaEC=-
log(
tan(0.5*atan(ecRad/zAbs)));
37 deltaPhi=-charge*asin(0.5*ecRad/Rcurv);
40 zAbs=(fabs(etaIP)/etaIP)*ecDist;
41 double Zflight=fabs(zAbs-vtxZ);
42 double alpha=(Zflight*ecRad)/(z*Rcurv);
43 double Rec=2*Rcurv*
sin(alpha/2);
44 deltaPhi=-charge*alpha/2;
45 etaEC=-
log(
tan(0.5*atan(Rec/ecDist)));
48 zNew=(fabs(etaIP)/etaIP)*ecDist;
49 double Zflight=fabs(zNew-vtxZ);
50 double Rvirt=fabs(Zflight*
tan(theta));
51 double Rec=2*Rcurv*
sin(Rvirt/(2*Rcurv));
52 deltaPhi=-(
charge)*(Rvirt/(2*Rcurv));
53 etaEC=-
log(
tan(0.5*atan(Rec/ecDist)));
56 if (zNew<0) etaEC=-etaEC;
59 if (phiEC<-acos(-1)) phiEC=2*acos(-1)+phiEC;
60 if (phiEC>acos(-1)) phiEC=-2*acos(-1)+phiEC;
62 std::pair<double,double> retVal(etaEC,phiEC);
69 if (fabs(eta1)<1.479) Rec=129;
71 double ce1=cosh(eta1);
72 double ce2=cosh(eta2);
73 double te1=tanh(eta1);
74 double te2=tanh(eta2);
76 double z=
cos(phi1-phi2)/ce1/ce2+te1*te2;
77 if(z!=0) dR=fabs(Rec*ce1*
sqrt(1./z/z-1.));
83 token_(consumes<edm::
HepMCProduct>(iConfig.getUntrackedParameter(
"ModuleLabel",edm::
InputTag(
"generator",
"unsmeared")))),
84 MaxSeedEta_(iConfig.getUntrackedParameter<double>(
"MaxSeedEta", 2.3)),
85 MinSeedMom_(iConfig.getUntrackedParameter<double>(
"MinSeedMom", 20.)),
86 MinIsolTrackMom_(iConfig.getUntrackedParameter<double>(
"MinIsolTrackMom",2.0)),
87 IsolCone_(iConfig.getUntrackedParameter<double>(
"IsolCone", 40.0)),
88 PixelEfficiency_(iConfig.getUntrackedParameter<double>(
"PixelEfficiency", 0.8))
93 throw cms::Exception(
"Configuration") <<
"PythiaFilterIsolatedTrack requires the RandomNumberGeneratorService\n";
114 const HepMC::GenEvent* myGenEvent = evt->GetEvent();
120 std::vector<const HepMC::GenParticle *> seeds;
123 for(HepMC::GenEvent::particle_const_iterator iter=myGenEvent->particles_begin(); iter!=myGenEvent->particles_end(); ++iter) {
125 int charge3 = pdt->particle(p->pdg_id())->
ID().threeCharge();
127 double momentum = p->momentum().rho();
128 double abseta = fabs(p->momentum().eta());
132 chargedParticles.push_back(p);
140 for(std::vector<const HepMC::GenParticle *>::const_iterator it1=seeds.begin(); it1!=seeds.end(); ++it1) {
144 p1->momentum().phi(),
145 p1->momentum().perp(),
146 (pdt->particle(p1->pdg_id()))->ID().threeCharge()/3,
151 for(std::vector<const HepMC::GenParticle *>::const_iterator it2=chargedParticles.begin(); it2!=chargedParticles.end(); ++it2) {
158 p2->momentum().phi(),
159 p2->momentum().perp(),
160 (pdt->particle(p2->pdg_id()))->ID().threeCharge()/3,
173 if(!failsIso)
return true;
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Sin< T >::type sin(const T &t)
Geom::Theta< T > theta() const
void getData(T &iHolder) const
Cos< T >::type cos(const T &t)
virtual bool filter(edm::Event &, const edm::EventSetup &)
edm::EDGetTokenT< edm::HepMCProduct > token_
Tan< T >::type tan(const T &t)
Abs< T >::type abs(const T &t)
static double getDistInCM(double eta1, double phi1, double eta2, double phi2)
PythiaFilterIsolatedTrack(const edm::ParameterSet &)
virtual CLHEP::HepRandomEngine & getEngine(StreamID const &) const =0
Use this engine in event methods.
edm::Service< edm::RandomNumberGenerator > rng_
~PythiaFilterIsolatedTrack()
StreamID streamID() const
static std::pair< double, double > GetEtaPhiAtEcal(double etaIP, double phiIP, double pT, int charge, double vtxZ)