Reference Manual

Go to the documentation of this file.
00001 #include "GeneratorInterface/GenFilters/interface/EMEnrichingFilterAlgo.h"
00002 
00003 #include "FWCore/Framework/interface/ESHandle.h"
00004 
00005 #include "DataFormats/Math/interface/deltaR.h"
00006 #include "DataFormats/Math/interface/deltaPhi.h"
00007 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
00008 #include "DataFormats/GeometrySurface/interface/Plane.h"
00009 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00010 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
00011 
00012 #include "MagneticField/Engine/interface/MagneticField.h"
00013 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00014 
00015 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
00016 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00017 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00018 
00019 #include "CLHEP/Vector/LorentzVector.h"
00020 
00021 
00022 using namespace edm;
00023 using namespace std;
00024 
00025 
00026 EMEnrichingFilterAlgo::EMEnrichingFilterAlgo(const edm::ParameterSet& iConfig) { 
00027 
00028   //set constants
00029   FILTER_TKISOCUT_=4;
00030   FILTER_CALOISOCUT_=7;
00031   FILTER_ETA_MIN_=0;
00032   FILTER_ETA_MAX_=2.4;
00033   ECALBARRELMAXETA_=1.479;
00034   ECALBARRELRADIUS_=129.0;
00035   ECALENDCAPZ_=304.5;
00036 
00037   
00038   
00039   isoGenParETMin_=(float) iConfig.getParameter<double>("isoGenParETMin");
00040   isoGenParConeSize_=(float) iConfig.getParameter<double>("isoGenParConeSize");
00041   clusterThreshold_=(float) iConfig.getParameter<double>("clusterThreshold");
00042   isoConeSize_=(float) iConfig.getParameter<double>("isoConeSize");
00043   hOverEMax_=(float) iConfig.getParameter<double>("hOverEMax");
00044   tkIsoMax_=(float) iConfig.getParameter<double>("tkIsoMax");
00045   caloIsoMax_=(float) iConfig.getParameter<double>("caloIsoMax");
00046   requireTrackMatch_=iConfig.getParameter<bool>("requireTrackMatch");
00047   genParSource_=iConfig.getParameter<edm::InputTag>("genParSource");
00048 
00049 }
00050 
00051 EMEnrichingFilterAlgo::~EMEnrichingFilterAlgo() {
00052 }
00053 
00054 
00055 bool EMEnrichingFilterAlgo::filter(const edm::Event& iEvent, const edm::EventSetup& iSetup)  {
00056 
00057 
00058   Handle<reco::GenParticleCollection> genParsHandle;
00059   iEvent.getByLabel(genParSource_,genParsHandle);
00060   reco::GenParticleCollection genPars=*genParsHandle;
00061 
00062   //bending of traj. of charged particles under influence of B-field
00063   std::vector<reco::GenParticle> genParsCurved=applyBFieldCurv(genPars,iSetup);
00064 
00065   bool result1=filterPhotonElectronSeed(clusterThreshold_,
00066                                         isoConeSize_,
00067                                         hOverEMax_,
00068                                         tkIsoMax_,
00069                                         caloIsoMax_,
00070                                         requireTrackMatch_,
00071                                         genPars,
00072                                         genParsCurved);
00073     
00074   bool result2=filterIsoGenPar(isoGenParETMin_,isoGenParConeSize_,genPars,genParsCurved);
00075 
00076 
00077   bool result=result1 || result2;
00078   
00079   return result;
00080 
00081 }
00082 
00083 
00084 
00085 //filter that uses clustering around photon and electron seeds
00086 //only electrons, photons, charged pions, and charged kaons are clustered
00087 //additional requirements:
00088 
00089 
00090 //seed threshold, total threshold, and cone size/shape are specified separately for the barrel and the endcap
00091 bool EMEnrichingFilterAlgo::filterPhotonElectronSeed(float clusterthreshold,
00092                                                  float isoConeSize,
00093                                                  float hOverEMax,
00094                                                  float tkIsoMax,
00095                                                  float caloIsoMax,
00096                                                  bool requiretrackmatch,
00097                                                  const std::vector<reco::GenParticle> &genPars,
00098                                                  const std::vector<reco::GenParticle> &genParsCurved) {
00099 
00100   float seedthreshold=5;
00101   float conesizeendcap=15;
00102 
00103   bool retval=false;
00104   
00105   vector<reco::GenParticle> seeds;
00106   //find electron and photon seeds - must have E>seedthreshold GeV
00107   for (uint32_t is=0;is<genParsCurved.size();is++) {
00108     reco::GenParticle gp=genParsCurved.at(is);
00109     if (gp.status()!=1 || fabs(gp.eta()) > FILTER_ETA_MAX_ || fabs(gp.eta()) < FILTER_ETA_MIN_) continue;
00110     int absid=abs(gp.pdgId());
00111     if (absid!=11 && absid!=22) continue;
00112     if (gp.et()>seedthreshold) seeds.push_back(gp);
00113   }
00114   
00115   bool matchtrack=false;
00116   
00117   //for every seed, try to cluster stable particles about it in cone
00118   for (uint32_t is=0;is<seeds.size();is++) {
00119     float eTInCone=0;//eT associated to the electron cluster
00120     float tkIsoET=0;//tracker isolation energy
00121     float caloIsoET=0;//calorimeter isolation energy
00122     float hadET=0;//isolation energy from heavy hadrons that goes in the same area as the "electron" - so contributes to H/E
00123     bool isBarrel=fabs(seeds.at(is).eta())<ECALBARRELMAXETA_;
00124     for (uint32_t ig=0;ig<genParsCurved.size();ig++) {
00125       reco::GenParticle gp=genParsCurved.at(ig);
00126       reco::GenParticle gpUnCurv=genPars.at(ig);//for tk isolation, p at vertex
00127       if (gp.status()!=1) continue;
00128       int gpabsid=abs(gp.pdgId());
00129       if (gp.et()<1) continue;//ignore very soft particles
00130       //BARREL
00131       if (isBarrel) {
00132         float dr=deltaR(seeds.at(is),gp);
00133         float dphi=deltaPhi(seeds.at(is).phi(),gp.phi());
00134         float deta=fabs(seeds.at(is).eta()-gp.eta());
00135         if (deta<0.03 && dphi<0.2) {
00136           if (gpabsid==22 || gpabsid==11 || gpabsid==211 || gpabsid==321) {
00137             //contributes to electron
00138             eTInCone+=gp.et();
00139             //check for a matched track with at least 5 GeV
00140             if ((gpabsid==11 || gpabsid==211 || gpabsid==321) && gp.et()>5) matchtrack=true;
00141           } else {
00142             //contributes to H/E
00143             hadET+=gp.et();
00144           }
00145         } else {
00146           float drUnCurv=deltaR(seeds.at(is),gpUnCurv);
00147           if ((gp.charge()==0 && dr<isoConeSize && gpabsid!=22) ||
00148               (gp.charge()!=0 && drUnCurv<isoConeSize)) {      
00149             //contributes to calo isolation energy
00150             caloIsoET+=gp.et();
00151           }
00152           if (gp.charge()!=0 && drUnCurv<isoConeSize) {
00153             //contributes to track isolation energy
00154             tkIsoET+=gp.et();
00155           }
00156         }
00157         //ENDCAP
00158       } else {
00159         float drxy=deltaRxyAtEE(seeds.at(is),gp);
00160         float dr=deltaR(seeds.at(is),gp);//the isolation is done in dR
00161         if (drxy<conesizeendcap) {
00162           if (gpabsid==22 || gpabsid==11 || gpabsid==211 || gpabsid==321) {
00163             //contributes to electron
00164             eTInCone+=gp.et();
00165             //check for a matched track with at least 5 GeV
00166             if ((gpabsid==11 || gpabsid==211 || gpabsid==321) && gp.et()>5) matchtrack=true;
00167           } else {
00168             //contributes to H/E
00169             hadET+=gp.et();
00170           }
00171         } else {
00172           float drUnCurv=deltaR(seeds.at(is),gpUnCurv);
00173           if ((gp.charge()==0 && dr<isoConeSize && gpabsid!=22) ||
00174               (gp.charge()!=0 && drUnCurv<isoConeSize)) {
00175             //contributes to calo isolation energy
00176             caloIsoET+=gp.et();
00177           }
00178           if (gp.charge()!=0 && drUnCurv<isoConeSize) {
00179             //contributes to track isolation energy
00180             tkIsoET+=gp.et();
00181           }
00182         }
00183       }
00184     }
00185 
00186     if (eTInCone>clusterthreshold && (!requiretrackmatch || matchtrack)) {
00187       //       cout <<"isoET: "<<isoET<<endl;
00188       if (hadET/eTInCone<hOverEMax && tkIsoET<tkIsoMax && caloIsoET<caloIsoMax) retval=true;
00189       break;
00190     }
00191   }
00192   
00193   return retval;
00194 }
00195 
00196 
00197 
00198 
00199 //make new genparticles vector taking into account the bending of charged particles in the b field
00200 //only stable-final-state (status==1) particles, with ET>=1 GeV, have their trajectories bent
00201 std::vector<reco::GenParticle> EMEnrichingFilterAlgo::applyBFieldCurv(const std::vector<reco::GenParticle> &genPars, const edm::EventSetup& iSetup) {
00202   
00203   
00204   vector<reco::GenParticle> curvedPars;
00205 
00206   edm::ESHandle<MagneticField> magField;
00207   iSetup.get<IdealMagneticFieldRecord>().get(magField);
00208 
00209   Cylinder::CylinderPointer theBarrel=Cylinder::build(Surface::PositionType(0,0,0),Surface::RotationType(),ECALBARRELRADIUS_);
00210   Plane::PlanePointer endCapPlus=Plane::build(Surface::PositionType(0,0,ECALENDCAPZ_),Surface::RotationType());
00211   Plane::PlanePointer endCapMinus=Plane::build(Surface::PositionType(0,0,-1*ECALENDCAPZ_),Surface::RotationType());
00212 
00213   AnalyticalPropagator propagator(&(*magField),alongMomentum);  
00214 
00215   for (uint32_t ig=0;ig<genPars.size();ig++) {
00216     reco::GenParticle gp=genPars.at(ig);
00217     //don't bend trajectories of neutral particles, unstable particles, particles with < 1 GeV
00218     //particles with < ~0.9 GeV don't reach the barrel
00219     //so just put them as-is into the new vector
00220     if (gp.charge()==0 || gp.status()!=1 || gp.et()<1) {
00221       curvedPars.push_back(gp);
00222       continue;
00223     }
00224     GlobalPoint vertex(gp.vx(),gp.vy(),gp.vz());
00225     GlobalVector gvect(gp.px(),gp.py(),gp.pz());
00226     FreeTrajectoryState fts(vertex,gvect,gp.charge(),&(*magField));
00227     TrajectoryStateOnSurface impact;
00228     //choose to propagate to barrel, +Z endcap, or -Z endcap, according to eta
00229     if (fabs(gp.eta())<ECALBARRELMAXETA_) {
00230       impact=propagator.propagate(fts,*theBarrel);
00231     } else if (gp.eta()>0) {
00232       impact=propagator.propagate(fts,*endCapPlus);
00233     } else {
00234       impact=propagator.propagate(fts,*endCapMinus);
00235     }
00236     //in case the particle doesn't reach the barrel/endcap, just put it as-is into the new vector
00237     //it should reach though.
00238     if (!impact.isValid()) {
00239       curvedPars.push_back(gp);
00240       continue;
00241     }
00242     math::XYZTLorentzVector newp4;
00243 
00244     //the magnitude of p doesn't change, only the direction
00245     //NB I do get some small change in magnitude by the following,
00246     //I think it is a numerical precision issue
00247     float et=gp.et();
00248     float phinew=impact.globalPosition().phi();
00249     float pxnew=et*cos(phinew);
00250     float pynew=et*sin(phinew);
00251     newp4.SetPx(pxnew);
00252     newp4.SetPy(pynew);
00253     newp4.SetPz(gp.pz());
00254     newp4.SetE(gp.energy());
00255     reco::GenParticle gpnew=gp;
00256     gpnew.setP4(newp4);
00257     curvedPars.push_back(gpnew);
00258   }
00259   return curvedPars;
00260 
00261 
00262 }
00263 
00264 //calculate the difference in the xy-plane positions of gp1 and gp1 at the ECAL endcap
00265 //if they go in different z directions returns 9999.
00266 float EMEnrichingFilterAlgo::deltaRxyAtEE(const reco::GenParticle &gp1, const reco::GenParticle &gp2) {
00267   
00268   if (gp1.pz()*gp2.pz() < 0) return 9999.;
00269   
00270   float rxy1=ECALENDCAPZ_*tan(gp1.theta());
00271   float x1=cos(gp1.phi())*rxy1;
00272   float y1=sin(gp1.phi())*rxy1;
00273 
00274   float rxy2=ECALENDCAPZ_*tan(gp2.theta());
00275   float x2=cos(gp2.phi())*rxy2;
00276   float y2=sin(gp2.phi())*rxy2;
00277   
00278   float dxy=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
00279   return dxy;
00280 
00281 
00282 }
00283 
00284 
00285 //filter looking for isolated charged pions, charged kaons, and electrons.
00286 //isolation done in cone of given size, looking at charged particles and neutral hadrons
00287 //photons aren't counted in the isolation requirements
00288 
00289 //need to have both the the curved and uncurved genpar collections
00290 //because tracker iso has to be treated differently than calo iso
00291 bool EMEnrichingFilterAlgo::filterIsoGenPar(float etmin, float conesize,const reco::GenParticleCollection &gph,
00292                                         const reco::GenParticleCollection &gphCurved
00293                                         ) {
00294   
00295   for (uint32_t ip=0;ip<gph.size();ip++) {
00296 
00297     reco::GenParticle gp=gph.at(ip);
00298     reco::GenParticle gpCurved=gphCurved.at(ip);    
00299     int gpabsid=abs(gp.pdgId());
00300     //find potential candidates
00301     if (gp.et()<=etmin || gp.status()!=1) continue;
00302     if (gpabsid!=11 && gpabsid != 211 && gpabsid!= 321) continue;
00303     if (fabs(gp.eta()) < FILTER_ETA_MIN_) continue;
00304     if (fabs(gp.eta()) > FILTER_ETA_MAX_) continue;
00305     
00306     //check isolation
00307     float tkiso=0;
00308     float caloiso=0;
00309     for (uint32_t jp=0;jp<gph.size();jp++) {
00310       if (jp==ip) continue;
00311       reco::GenParticle pp=gph.at(jp);
00312       reco::GenParticle ppCurved=gphCurved.at(jp);
00313       if (pp.status() != 1) continue;
00314       float dr=deltaR(gp,pp);
00315       float drCurved=deltaR(gpCurved,ppCurved);
00316       if (abs(pp.charge())==1 && pp.et()>2 && dr<conesize) {
00317         tkiso+=pp.et();
00318       }
00319       if (pp.et()>2 && abs(pp.pdgId())!=22 && drCurved<conesize) {
00320         caloiso+=pp.energy();
00321       }
00322     }
00323     if (tkiso<FILTER_TKISOCUT_ && caloiso<FILTER_CALOISOCUT_) return true;
00324   }
00325   return false;
00326 }
00327
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/GeneratorInterface/GenFilters/src/EMEnrichingFilterAlgo.cc
