#include <PythiaFilterIsolatedTrack.h>

Inheritance diagram for PythiaFilterIsolatedTrack:

Public Member Functions
virtual bool	filter (edm::Event &, const edm::EventSetup &)

	PythiaFilterIsolatedTrack (const edm::ParameterSet &)

	~PythiaFilterIsolatedTrack ()

Public Member Functions inherited from edm::EDFilter
	EDFilter ()

virtual	~EDFilter ()

Public Member Functions inherited from edm::ProducerBase
	ProducerBase ()

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

boost::function< void(const BranchDescription &)>	registrationCallback () const
	used by the fwk to register list of products More...

virtual	~ProducerBase ()

Static Public Member Functions
static double	getDistInCM (double eta1, double phi1, double eta2, double phi2)

static std::pair< double, double >	GetEtaPhiAtEcal (double etaIP, double phiIP, double pT, int charge, double vtxZ)

Static Public Member Functions inherited from edm::EDFilter
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &)

Private Attributes
CLHEP::RandFlat *	flatDistribution_

double	IsolCone_

double	MaxSeedEta_

double	MinIsolTrackMom_

double	MinSeedMom_

std::string	ModuleLabel_

double	PixelEfficiency_

edm::Service < edm::RandomNumberGenerator >	rng_

Additional Inherited Members
Public Types inherited from edm::EDFilter
typedef EDFilter	ModuleType

typedef WorkerT< EDFilter >	WorkerType

Public Types inherited from edm::ProducerBase
typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

Protected Member Functions inherited from edm::EDFilter
CurrentProcessingContext const *	currentContext () const

Protected Member Functions inherited from edm::ProducerBase
template<class TProducer , class TMethod >
void	callWhenNewProductsRegistered (TProducer *iProd, TMethod iMethod)

Detailed Description

PythiaFilterGammaJet filter implements generator-level preselections for ChargedHadron+jet like events to be used in jet energy calibration. Ported from fortran code written by V.Konoplianikov.

Author: O.Kodolova, SINP

Cleaned up code, and added pixel efficiency functionality

Author: J.P. Chou, Brown University

Definition at line 33 of file PythiaFilterIsolatedTrack.h.

Constructor & Destructor Documentation

PythiaFilterIsolatedTrack::PythiaFilterIsolatedTrack ( const edm::ParameterSet & iConfig )

explicit

Definition at line 79 of file PythiaFilterIsolatedTrack.cc.

References edm::hlt::Exception, flatDistribution_, edm::RandomNumberGenerator::getEngine(), edm::Service< T >::isAvailable(), and rng_.

                                                                                    :
   ModuleLabel_(iConfig.getUntrackedParameter("ModuleLabel",std::string("generator"))),
   MaxSeedEta_(iConfig.getUntrackedParameter<double>("MaxSeedEta", 2.3)),
   MinSeedMom_(iConfig.getUntrackedParameter<double>("MinSeedMom", 20.)),
   MinIsolTrackMom_(iConfig.getUntrackedParameter<double>("MinIsolTrackMom",2.0)),
   IsolCone_(iConfig.getUntrackedParameter<double>("IsolCone", 40.0)),
   PixelEfficiency_(iConfig.getUntrackedParameter<double>("PixelEfficiency", 0.8))
 { 
 
   // initialize the random number generator service
   if(!rng_.isAvailable()) {
     throw cms::Exception("Configuration") << "PythiaFilterIsolatedTrack requires the RandomNumberGeneratorService\n";
   }
   CLHEP::HepRandomEngine& engine = rng_->getEngine();
   flatDistribution_ = new CLHEP::RandFlat(engine, 0.0, 1.0);
 }

PythiaFilterIsolatedTrack::~PythiaFilterIsolatedTrack ( )

Definition at line 97 of file PythiaFilterIsolatedTrack.cc.

References flatDistribution_.

 {
   delete flatDistribution_;
 }

Member Function Documentation

bool PythiaFilterIsolatedTrack::filter	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

virtual

Implements edm::EDFilter.

Definition at line 104 of file PythiaFilterIsolatedTrack.cc.

References abs, MBUEandQCDValidation_cff::chargedParticles, flatDistribution_, configurableAnalysis::GenParticle, edm::Event::getByLabel(), edm::EventSetup::getData(), getDistInCM(), GetEtaPhiAtEcal(), IsolCone_, MaxSeedEta_, MinIsolTrackMom_, MinSeedMom_, ModuleLabel_, AlCaHLTBitMon_ParallelJobs::p, p1, p2, PixelEfficiency_, and ntuplemaker::status.

                                                                                    {
 
   edm::ESHandle<ParticleDataTable> pdt;
   iSetup.getData( pdt );
 
   edm::Handle<edm::HepMCProduct> evt;
   iEvent.getByLabel(ModuleLabel_, evt);
 
   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
 
   // all of the stable, charged particles with momentum>MinIsolTrackMom_ and |eta|<MaxSeedEta_+0.5
   std::vector<const HepMC::GenParticle *> chargedParticles;
 
   // all of the stable, charged particles with momentum>MinSeedMom_ and |eta|<MaxSeedEta_
   std::vector<const HepMC::GenParticle *> seeds;
 
   // fill the vector of charged particles and seeds in the event
   for(HepMC::GenEvent::particle_const_iterator iter=myGenEvent->particles_begin(); iter!=myGenEvent->particles_end(); ++iter) {
     const HepMC::GenParticle *p=*iter;
     int charge3 = pdt->particle(p->pdg_id())->ID().threeCharge();
     int status = p->status();
     double momentum = p->momentum().rho();
     double abseta = fabs(p->momentum().eta());
 
     // only consider stable, charged particles
     if(abs(charge3)==3 && status==1 && momentum>MinIsolTrackMom_ && abseta<MaxSeedEta_+0.5) {
       chargedParticles.push_back(p);
       if(momentum>MinSeedMom_ && abseta<MaxSeedEta_) {
     seeds.push_back(p);
       }
     }
   }
 
   // loop over all the seeds and see if any of them are isolated
   for(std::vector<const HepMC::GenParticle *>::const_iterator it1=seeds.begin(); it1!=seeds.end(); ++it1) {
     const HepMC::GenParticle *p1=*it1;
 
     std::pair<double,double> EtaPhi1=GetEtaPhiAtEcal(p1->momentum().eta(),
                              p1->momentum().phi(),
                              p1->momentum().perp(),
                              (pdt->particle(p1->pdg_id()))->ID().threeCharge()/3,
                              0.0);
 
     // loop over all of the other charged particles in the event, and see if any are close by
     bool failsIso=false;
     for(std::vector<const HepMC::GenParticle *>::const_iterator it2=chargedParticles.begin(); it2!=chargedParticles.end(); ++it2) {
       const HepMC::GenParticle *p2=*it2;
 
       // don't consider the seed particle among the other charge particles
       if(p1==p2) continue;
 
       std::pair<double,double> EtaPhi2=GetEtaPhiAtEcal(p2->momentum().eta(),
                                p2->momentum().phi(),
                                p2->momentum().perp(),
                                (pdt->particle(p2->pdg_id()))->ID().threeCharge()/3,
                                0.0);
 
       // find out how far apart the particles are
       // if the seed fails the isolation requirement, try a different seed
       // occasionally allow a seed to pass to isolation requirement
       if(getDistInCM(EtaPhi1.first, EtaPhi1.second, EtaPhi2.first, EtaPhi2.second) < IsolCone_ &&
      flatDistribution_->fire() < PixelEfficiency_) {
     failsIso=true;
     break;
       }
     }
 
     if(!failsIso) return true;
 
   } //loop over seeds
 
   return false;
 }

double PythiaFilterIsolatedTrack::getDistInCM	(	double	eta1,
		double	phi1,
		double	eta2,
		double	phi2
	)

static

Definition at line 63 of file PythiaFilterIsolatedTrack.cc.

References funct::cos(), PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi::dR, mathSSE::sqrt(), and detailsBasic3DVector::z.

Referenced by filter().

 {
   double dR, Rec;
   if (fabs(eta1)<1.479) Rec=129;
   else Rec=317;
   double ce1=cosh(eta1);
   double ce2=cosh(eta2);
   double te1=tanh(eta1);
   double te2=tanh(eta2);
   
   double z=cos(phi1-phi2)/ce1/ce2+te1*te2;
   if(z!=0) dR=fabs(Rec*ce1*sqrt(1./z/z-1.));
   else     dR=999999.;
   return dR;
 }

std::pair< double, double > PythiaFilterIsolatedTrack::GetEtaPhiAtEcal	(	double	etaIP,
		double	phiIP,
		double	pT,
		int	charge,
		double	vtxZ
	)

static

Definition at line 14 of file PythiaFilterIsolatedTrack.cc.

References alpha, DeDxDiscriminatorTools::charge(), SiPixelRawToDigiRegional_cfi::deltaPhi, create_public_lumi_plots::exp, create_public_lumi_plots::log, funct::sin(), funct::tan(), theta(), and detailsBasic3DVector::z.

Referenced by filter().

 {
   double deltaPhi;
   double etaEC=100;
   double phiEC=100;
   double Rcurv=pT*33.3*100/38; //r(m)=pT(GeV)*33.3/B(kG)
   double ecDist=317;  //distance to ECAL andcap from IP (cm), 317 - ecal (not preshower), preshower -300
   double ecRad=129;  //radius of ECAL barrel (cm)
   double theta=2*atan(exp(-etaIP));
   double zNew;
   if (theta>0.5*acos(-1)) theta=acos(-1)-theta;
   if (fabs(etaIP)<1.479) {
     deltaPhi=-charge*asin(0.5*ecRad/Rcurv);
     double alpha1=2*asin(0.5*ecRad/Rcurv);
     double z=ecRad/tan(theta);
     if (etaIP>0) zNew=z*(Rcurv*alpha1)/ecRad+vtxZ; //new z-coordinate of track
     else  zNew=-z*(Rcurv*alpha1)/ecRad+vtxZ; //new z-coordinate of track
     double zAbs=fabs(zNew);
     if (zAbs<ecDist) {
       etaEC=-log(tan(0.5*atan(ecRad/zAbs)));
       deltaPhi=-charge*asin(0.5*ecRad/Rcurv);
     }
     if (zAbs>ecDist) {
       zAbs=(fabs(etaIP)/etaIP)*ecDist;
       double Zflight=fabs(zAbs-vtxZ);
       double alpha=(Zflight*ecRad)/(z*Rcurv);
       double Rec=2*Rcurv*sin(alpha/2);
       deltaPhi=-charge*alpha/2;
       etaEC=-log(tan(0.5*atan(Rec/ecDist)));
     }
   } else {
     zNew=(fabs(etaIP)/etaIP)*ecDist;
     double Zflight=fabs(zNew-vtxZ);
     double Rvirt=fabs(Zflight*tan(theta));
     double Rec=2*Rcurv*sin(Rvirt/(2*Rcurv));
     deltaPhi=-(charge)*(Rvirt/(2*Rcurv));
     etaEC=-log(tan(0.5*atan(Rec/ecDist)));
   }
   
   if (zNew<0) etaEC=-etaEC;
   phiEC=phiIP+deltaPhi;
   
   if (phiEC<-acos(-1)) phiEC=2*acos(-1)+phiEC;
   if (phiEC>acos(-1)) phiEC=-2*acos(-1)+phiEC;
   
   std::pair<double,double> retVal(etaEC,phiEC);
   return retVal;
 }