Reference Manual

Go to the documentation of this file.
00001 #include <vector>
00002 #include <memory>
00003 #include <algorithm>
00004 
00005 // Class header file
00006 #include "Calibration/HcalIsolatedTrackReco/interface/IsolatedPixelTrackCandidateProducer.h"
00007 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidateFwd.h"
00008 // Framework
00009 #include "DataFormats/Common/interface/Handle.h"
00010 #include "FWCore/Framework/interface/ESHandle.h"
00011 #include "FWCore/Framework/interface/ESTransientHandle.h"
00012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00013 #include "FWCore/Utilities/interface/Exception.h"
00014 //
00015 #include "DataFormats/Common/interface/TriggerResults.h"
00016 // L1Extra
00017 #include "DataFormats/L1Trigger/interface/L1EmParticle.h"
00018 #include "DataFormats/L1Trigger/interface/L1JetParticleFwd.h"
00020 
00021 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
00022 
00023 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
00024 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
00025 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
00026 //#include "DataFormats/L1GlobalTrigger/interface/L1GtLogicParser.h"
00027 
00028 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapFwd.h"
00029 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMap.h"
00030 
00031 // Math
00032 #include "Math/GenVector/VectorUtil.h"
00033 #include "Math/GenVector/PxPyPzE4D.h"
00034 #include "DataFormats/Math/interface/deltaR.h"
00035 
00036 //vertices
00037 #include "DataFormats/VertexReco/interface/Vertex.h"
00038 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00039 
00040 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00041 #include "DetectorDescription/Core/interface/DDLogicalPart.h"
00042 #include "DetectorDescription/Core/interface/DDSolid.h"
00043 
00044 //magF
00045 #include "MagneticField/Engine/interface/MagneticField.h"
00046 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00047 #include "MagneticField/VolumeBasedEngine/interface/VolumeBasedMagneticField.h"
00048 
00049 //for ECAL geometry
00050 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
00051 #include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
00052 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00053 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00054 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00055 
00056 
00057 IsolatedPixelTrackCandidateProducer::IsolatedPixelTrackCandidateProducer(const edm::ParameterSet& config){
00058    
00059   l1eTauJetsSource_           = config.getParameter<edm::InputTag>("L1eTauJetsSource");
00060   tauAssocCone_               = config.getParameter<double>("tauAssociationCone"); 
00061   tauUnbiasCone_              = config.getParameter<double>("tauUnbiasCone");
00062   pixelTracksSources_         = config.getParameter<std::vector<edm::InputTag> >("PixelTracksSources");
00063   prelimCone_                 = config.getParameter<double>("ExtrapolationConeSize");
00064   pixelIsolationConeSizeAtEC_ = config.getParameter<double>("PixelIsolationConeSizeAtEC");
00065   hltGTseedlabel_             = config.getParameter<edm::InputTag>("L1GTSeedLabel");
00066   vtxCutSeed_                 = config.getParameter<double>("MaxVtxDXYSeed");
00067   vtxCutIsol_                 = config.getParameter<double>("MaxVtxDXYIsol");
00068   vertexLabel_                = config.getParameter<edm::InputTag>("VertexLabel");
00069   bfield_                     = config.getParameter<std::string>("MagFieldRecordName");
00070   minPTrackValue_             = config.getParameter<double>("minPTrack");
00071   maxPForIsolationValue_      = config.getParameter<double>("maxPTrackForIsolation");
00072   ebEtaBoundary_              = config.getParameter<double>("EBEtaBoundary");
00073   rEB_ = zEE_ = -1;
00074 
00075   // Register the product
00076   produces< reco::IsolatedPixelTrackCandidateCollection >();
00077 }
00078 
00079 IsolatedPixelTrackCandidateProducer::~IsolatedPixelTrackCandidateProducer() {
00080 
00081 }
00082 
00083 void IsolatedPixelTrackCandidateProducer::beginJob () {}
00084 
00085 void IsolatedPixelTrackCandidateProducer::beginRun(edm::Run &run, const edm::EventSetup &theEventSetup)
00086 {
00087 
00088   edm::ESHandle<CaloGeometry> pG;
00089   theEventSetup.get<CaloGeometryRecord>().get(pG);   
00090   
00091   const double rad (dynamic_cast<const EcalBarrelGeometry*>( pG->getSubdetectorGeometry(DetId::Ecal, EcalBarrel ))->avgRadiusXYFrontFaceCenter() ) ;
00092   
00093   const double zz (dynamic_cast<const EcalEndcapGeometry*>( pG->getSubdetectorGeometry(DetId::Ecal, EcalEndcap ))->avgAbsZFrontFaceCenter() ) ;
00094 
00095   rEB_=rad;
00096   zEE_=zz;
00097 
00098 }
00099 
00100 void IsolatedPixelTrackCandidateProducer::produce(edm::Event& theEvent, const edm::EventSetup& theEventSetup) {
00101 
00102   reco::IsolatedPixelTrackCandidateCollection* trackCollection=new reco::IsolatedPixelTrackCandidateCollection;
00103 
00104   //create vector of refs from input collections
00105   std::vector<reco::TrackRef> pixelTrackRefs;
00106 
00107   for (unsigned int iPix=0; iPix<pixelTracksSources_.size(); iPix++)
00108     {
00109       edm::Handle<reco::TrackCollection> iPixCol;
00110       theEvent.getByLabel(pixelTracksSources_[iPix],iPixCol);
00111       for (reco::TrackCollection::const_iterator pit=iPixCol->begin(); pit!=iPixCol->end(); pit++)
00112         {
00113           pixelTrackRefs.push_back(reco::TrackRef(iPixCol,pit-iPixCol->begin()));
00114         }
00115     }
00116 
00117   edm::Handle<l1extra::L1JetParticleCollection> l1eTauJets;
00118   theEvent.getByLabel(l1eTauJetsSource_,l1eTauJets);
00119 
00120   edm::Handle<reco::VertexCollection> pVert;
00121   theEvent.getByLabel(vertexLabel_,pVert);
00122 
00123   double ptTriggered  = -10;
00124   double etaTriggered = -100;
00125   double phiTriggered = -100;
00126   
00127   edm::Handle<trigger::TriggerFilterObjectWithRefs> l1trigobj;
00128   theEvent.getByLabel(hltGTseedlabel_, l1trigobj);
00129   
00130   std::vector< edm::Ref<l1extra::L1JetParticleCollection> > l1tauobjref;
00131   std::vector< edm::Ref<l1extra::L1JetParticleCollection> > l1jetobjref;
00132   std::vector< edm::Ref<l1extra::L1JetParticleCollection> > l1forjetobjref;
00133   
00134   l1trigobj->getObjects(trigger::TriggerL1TauJet, l1tauobjref);
00135   l1trigobj->getObjects(trigger::TriggerL1CenJet, l1jetobjref);
00136   l1trigobj->getObjects(trigger::TriggerL1ForJet, l1forjetobjref);
00137   
00138   for (unsigned int p=0; p<l1tauobjref.size(); p++)
00139     {
00140       if (l1tauobjref[p]->pt()>ptTriggered)
00141         {
00142           ptTriggered  = l1tauobjref[p]->pt(); 
00143           phiTriggered = l1tauobjref[p]->phi();
00144           etaTriggered = l1tauobjref[p]->eta();
00145         }
00146     }
00147   for (unsigned int p=0; p<l1jetobjref.size(); p++)
00148     {
00149       if (l1jetobjref[p]->pt()>ptTriggered)
00150         {
00151           ptTriggered  = l1jetobjref[p]->pt();
00152           phiTriggered = l1jetobjref[p]->phi();
00153           etaTriggered = l1jetobjref[p]->eta();
00154         }
00155     }
00156   for (unsigned int p=0; p<l1forjetobjref.size(); p++)
00157     {
00158       if (l1forjetobjref[p]->pt()>ptTriggered)
00159         {
00160           ptTriggered=l1forjetobjref[p]->pt();
00161           phiTriggered=l1forjetobjref[p]->phi();
00162           etaTriggered=l1forjetobjref[p]->eta();
00163         }
00164     }
00165 
00166   double minPTrack_    = minPTrackValue_;
00167   double drMaxL1Track_ = tauAssocCone_;
00168   
00169   int ntr = 0;
00170   
00171   //loop to select isolated tracks
00172   for (unsigned iSeed=0; iSeed<pixelTrackRefs.size(); iSeed++)
00173     {
00174       if(pixelTrackRefs[iSeed]->p()<minPTrack_) continue;
00175 
00176       bool good     = false;
00177       bool vtxMatch = false;
00178 
00179       //associate to vertex (in Z) 
00180       reco::VertexCollection::const_iterator vitSel;
00181       double minDZ = 100;
00182       for (reco::VertexCollection::const_iterator vit=pVert->begin(); vit!=pVert->end(); vit++)
00183         {
00184           if (fabs(pixelTrackRefs[iSeed]->dz(vit->position()))<minDZ)
00185             {
00186               minDZ  = fabs(pixelTrackRefs[iSeed]->dz(vit->position()));
00187               vitSel = vit;
00188             }
00189         }
00190       //cut on dYX:
00191       if (minDZ!=100&&fabs(pixelTrackRefs[iSeed]->dxy(vitSel->position()))<vtxCutSeed_) vtxMatch=true;
00192       if (minDZ==100) vtxMatch=true;
00193 
00194       //select tracks not matched to triggered L1 jet
00195       double R=deltaR(etaTriggered, phiTriggered, pixelTrackRefs[iSeed]->eta(), pixelTrackRefs[iSeed]->phi());
00196       if (R<tauUnbiasCone_) continue;
00197       
00198       //check taujet matching
00199       bool tmatch=false;
00200       l1extra::L1JetParticleCollection::const_iterator selj;
00201       for (l1extra::L1JetParticleCollection::const_iterator tj=l1eTauJets->begin(); tj!=l1eTauJets->end(); tj++) 
00202         {
00203           if(ROOT::Math::VectorUtil::DeltaR(pixelTrackRefs[iSeed]->momentum(),tj->momentum()) > drMaxL1Track_) continue;
00204           selj   = tj;
00205           tmatch = true;
00206         } //loop over L1 tau
00207       
00208       //propagate seed track to ECAL surface:
00209       std::pair<double,double> seedCooAtEC;
00210       // in case vertex is found:
00211       if (minDZ!=100) seedCooAtEC=GetEtaPhiAtEcal(theEventSetup, pixelTrackRefs[iSeed]->eta(), pixelTrackRefs[iSeed]->phi(), pixelTrackRefs[iSeed]->pt(), pixelTrackRefs[iSeed]->charge(), vitSel->z());
00212       //in case vertex is not found:
00213       else seedCooAtEC=GetEtaPhiAtEcal(theEventSetup, pixelTrackRefs[iSeed]->eta(), pixelTrackRefs[iSeed]->phi(), pixelTrackRefs[iSeed]->pt(), pixelTrackRefs[iSeed]->charge(), 0);
00214 
00215       //calculate isolation
00216       double maxP = 0;
00217       double sumP = 0;
00218       for (unsigned iSurr=0; iSurr<pixelTrackRefs.size(); iSurr++)
00219         {
00220           if(iSeed==iSurr) continue;
00221           //define preliminary cone around seed track impact point from which tracks will be extrapolated:
00222           if (deltaR(seedCooAtEC.first, seedCooAtEC.second, pixelTrackRefs[iSurr]->eta(), pixelTrackRefs[iSurr]->phi())>prelimCone_) continue;
00223           //associate to vertex (in Z):
00224           double minDZ2=100;
00225           reco::VertexCollection::const_iterator vitSel2;
00226           for (reco::VertexCollection::const_iterator vit=pVert->begin(); vit!=pVert->end(); vit++)
00227             {
00228               if (fabs(pixelTrackRefs[iSurr]->dz(vit->position()))<minDZ2)
00229                 {
00230                   minDZ2  = fabs(pixelTrackRefs[iSurr]->dz(vit->position()));
00231                   vitSel2 = vit;
00232                 }
00233             }
00234           //cut ot dXY:
00235           if (minDZ2!=100&&fabs(pixelTrackRefs[iSurr]->dxy(vitSel2->position()))>vtxCutIsol_) continue;
00236           //propagate to ECAL surface:
00237           std::pair<double,double> cooAtEC;
00238           // in case vertex is found:
00239           if (minDZ2!=100) cooAtEC=GetEtaPhiAtEcal(theEventSetup, pixelTrackRefs[iSurr]->eta(), pixelTrackRefs[iSurr]->phi(), pixelTrackRefs[iSurr]->pt(), pixelTrackRefs[iSurr]->charge(), vitSel2->z());
00240           // in case vertex is not found:
00241           else cooAtEC=GetEtaPhiAtEcal(theEventSetup, pixelTrackRefs[iSurr]->eta(), pixelTrackRefs[iSurr]->phi(), pixelTrackRefs[iSurr]->pt(), pixelTrackRefs[iSurr]->charge(), 0);
00242           
00243           //calculate distance at ECAL surface and update isolation: 
00244           if (getDistInCM(seedCooAtEC.first, seedCooAtEC.second, cooAtEC.first, cooAtEC.second)<pixelIsolationConeSizeAtEC_)
00245             {
00246               sumP+=pixelTrackRefs[iSurr]->p();
00247               if(pixelTrackRefs[iSurr]->p()>maxP) maxP=pixelTrackRefs[iSurr]->p();
00248             }
00249         }
00250 
00251       if (tmatch||vtxMatch) good=true;
00252 
00253       if (good&&maxP<maxPForIsolationValue_)
00254         {
00255           reco::IsolatedPixelTrackCandidate newCandidate(pixelTrackRefs[iSeed], l1extra::L1JetParticleRef(l1eTauJets,selj-l1eTauJets->begin()), maxP, sumP);
00256           trackCollection->push_back(newCandidate);
00257           ntr++;
00258         }
00259     }//loop over pixel tracks
00260 
00261   // put the product in the event
00262   std::auto_ptr< reco::IsolatedPixelTrackCandidateCollection > outCollection(trackCollection);
00263   theEvent.put(outCollection);
00264 
00265 }
00266 
00267 
00268 double IsolatedPixelTrackCandidateProducer::getDistInCM(double eta1, double phi1, double eta2, double phi2)
00269 {
00270   double Rec;
00271   double theta1=2*atan(exp(-eta1));
00272   double theta2=2*atan(exp(-eta2));
00273   if (fabs(eta1)<1.479) Rec=129; //radius of ECAL barrel
00274   else if (fabs(eta1)>1.479&&fabs(eta1)<7.0) Rec=tan(theta1)*317; //distance from IP to ECAL endcap
00275   else return 1000;
00276 
00277   //|vect| times tg of acos(scalar product)
00278   double angle=acos((sin(theta1)*sin(theta2)*(sin(phi1)*sin(phi2)+cos(phi1)*cos(phi2))+cos(theta1)*cos(theta2)));
00279   if (angle<acos(-1)/2) return fabs((Rec/sin(theta1))*tan(angle));
00280   else return 1000;
00281 }
00282 
00283 
00284 std::pair<double,double>
00285 IsolatedPixelTrackCandidateProducer::GetEtaPhiAtEcal(const edm::EventSetup& iSetup, double etaIP, double phiIP, double pT, int charge, double vtxZ)
00286 {
00287   edm::ESHandle<MagneticField> vbfField;
00288   iSetup.get<IdealMagneticFieldRecord>().get(vbfField);
00289   const VolumeBasedMagneticField* vbfCPtr = dynamic_cast<const VolumeBasedMagneticField*>(&(*vbfField));
00290   GlobalVector BField=vbfCPtr->inTesla(GlobalPoint(0,0,0));
00291  //test
00292  //int curvSgn=int(BField.z()/fabs(BField.z())); 
00293 
00294   double bfVal=BField.mag();
00295 
00296   double deltaPhi=0;
00297   double etaEC = 100;
00298   double phiEC = 100;
00299 
00300   double Rcurv = 9999999;
00301   if (bfVal!=0) Rcurv=pT*33.3*100/(bfVal*10); //r(m)=pT(GeV)*33.3/B(kG)
00302 
00303   double ecDist = zEE_;  //distance to ECAL andcap from IP (cm), 317 - ecal (not preshower), preshower -300
00304   double ecRad  = rEB_;  //radius of ECAL barrel (cm)
00305   double theta=2*atan(exp(-etaIP));
00306   double zNew=0;
00307   if (theta>0.5*acos(-1)) theta=acos(-1)-theta;
00308   if (fabs(etaIP)<ebEtaBoundary_)
00309     {
00310       if ((0.5*ecRad/Rcurv)>1)
00311         {
00312           etaEC=10000;
00313           deltaPhi=0;
00314         }
00315       else
00316         {
00317           deltaPhi      =-charge*asin(0.5*ecRad/Rcurv);
00318           double alpha1 = 2*asin(0.5*ecRad/Rcurv);
00319           double z      = ecRad/tan(theta);
00320           if (etaIP>0) zNew = z*(Rcurv*alpha1)/ecRad+vtxZ; //new z-coordinate of track
00321           else         zNew =-z*(Rcurv*alpha1)/ecRad+vtxZ; //new z-coordinate of track
00322           double zAbs=fabs(zNew);
00323           if (zAbs<ecDist)
00324             {
00325               etaEC    = -log(tan(0.5*atan(ecRad/zAbs)));
00326               deltaPhi = -charge*asin(0.5*ecRad/Rcurv);
00327             }
00328           if (zAbs>ecDist)
00329             {
00330               zAbs           = (fabs(etaIP)/etaIP)*ecDist;
00331               double Zflight = fabs(zAbs-vtxZ);
00332               double alpha   = (Zflight*ecRad)/(z*Rcurv);
00333               double Rec     = 2*Rcurv*sin(alpha/2);
00334               deltaPhi       =-charge*alpha/2;
00335               etaEC          =-log(tan(0.5*atan(Rec/ecDist)));
00336             }
00337         }
00338     }
00339   else
00340     {
00341       zNew           = (fabs(etaIP)/etaIP)*ecDist;
00342       double Zflight = fabs(zNew-vtxZ);
00343       double Rvirt   = fabs(Zflight*tan(theta));
00344       double Rec     = 2*Rcurv*sin(Rvirt/(2*Rcurv));
00345       deltaPhi       =-(charge)*(Rvirt/(2*Rcurv));
00346       etaEC          =-log(tan(0.5*atan(Rec/ecDist)));
00347     }
00348 
00349   if (zNew<0) etaEC=-etaEC;
00350   phiEC            = phiIP+deltaPhi;
00351 
00352   if (phiEC<-acos(-1)) phiEC = 2*acos(-1)+phiEC;
00353   if (phiEC>acos(-1))  phiEC =-2*acos(-1)+phiEC;
00354 
00355   std::pair<double,double> retVal(etaEC,phiEC);
00356   return retVal;
00357 }
00358 
00359
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/Calibration/HcalIsolatedTrackReco/src/IsolatedPixelTrackCandidateProducer.cc
