CMSSW Reference Manual

00001 // -*- C++ -*-
00002 //
00003 // Package:    EgammaElectronAlgos
00004 // Class:      ElectronPixelSeedGenerator.
00005 // 
00013 //
00014 // Original Author:  Ursula Berthon, Claude Charlot
00015 //         Created:  Mon Mar 27 13:22:06 CEST 2006
00016 // $Id: ElectronPixelSeedGenerator.cc,v 1.54 2008/05/27 11:37:22 elmer Exp $
00017 //
00018 //
00019 
00020 #include "RecoEgamma/EgammaElectronAlgos/interface/PixelHitMatcher.h" 
00021 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronPixelSeedGenerator.h" 
00022 
00023 #include "RecoTracker/TransientTrackingRecHit/interface/TSiPixelRecHit.h" 
00024 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h" 
00025 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
00026 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
00027 #include "RecoTracker/TkNavigation/interface/SimpleNavigationSchool.h"
00028 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
00029 #include "RecoTracker/Record/interface/NavigationSchoolRecord.h"
00030 
00031 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00032 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
00033 //#include "DataFormats/EgammaReco/interface/EcalCluster.h"
00034 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00035 
00036 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
00037 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
00038 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
00039 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00040 #include "DataFormats/Common/interface/Handle.h"
00041 
00042 #include "MagneticField/Engine/interface/MagneticField.h"
00043 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00044 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00045 
00046 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
00047 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00048 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
00049 
00050 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00051 #include "FWCore/Framework/interface/ESHandle.h"
00052 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00053 
00054 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
00055 #include <vector>
00056 #include <utility>
00057 
00058 ElectronPixelSeedGenerator::ElectronPixelSeedGenerator(const edm::ParameterSet &pset)
00059   :   dynamicphiroad_(pset.getParameter<bool>("dynamicPhiRoad")),
00060       fromTrackerSeeds_(pset.getParameter<bool>("fromTrackerSeeds")),
00061       //      initialSeeds_(pset.getParameter<edm::InputTag>("initialSeeds")),
00062       lowPtThreshold_(pset.getParameter<double>("LowPtThreshold")),
00063       highPtThreshold_(pset.getParameter<double>("HighPtThreshold")),
00064       sizeWindowENeg_(pset.getParameter<double>("SizeWindowENeg")),
00065       phimin2_(pset.getParameter<double>("PhiMin2")),      
00066       phimax2_(pset.getParameter<double>("PhiMax2")),
00067       deltaPhi1Low_(pset.getParameter<double>("DeltaPhi1Low")),
00068       deltaPhi1High_(pset.getParameter<double>("DeltaPhi1High")),
00069       deltaPhi2_(pset.getParameter<double>("DeltaPhi2")),
00070       myMatchEle(0), myMatchPos(0),
00071       thePropagator(0), theMeasurementTracker(0), 
00072       theSetup(0), pts_(0),
00073       cacheIDMagField_(0),cacheIDGeom_(0),cacheIDNavSchool_(0),cacheIDCkfComp_(0),cacheIDTrkGeom_(0)
00074 { 
00075      // Instantiate the pixel hit matchers
00076 
00077   myMatchEle = new PixelHitMatcher( pset.getParameter<double>("ePhiMin1"), 
00078                                     pset.getParameter<double>("ePhiMax1"),
00079                                     pset.getParameter<double>("PhiMin2"),
00080                                     pset.getParameter<double>("PhiMax2"),
00081                                     pset.getParameter<double>("z2MinB"),
00082                                     pset.getParameter<double>("z2MaxB"),
00083                                     pset.getParameter<double>("r2MinF"),
00084                                     pset.getParameter<double>("r2MaxF"),
00085                                     pset.getParameter<double>("rMinI"),
00086                                     pset.getParameter<double>("rMaxI"),
00087                                     pset.getParameter<bool>("searchInTIDTEC"));
00088 
00089   myMatchPos = new PixelHitMatcher( pset.getParameter<double>("pPhiMin1"),
00090                                     pset.getParameter<double>("pPhiMax1"),
00091                                     pset.getParameter<double>("PhiMin2"),
00092                                     pset.getParameter<double>("PhiMax2"),
00093                                     pset.getParameter<double>("z2MinB"),
00094                                     pset.getParameter<double>("z2MaxB"),
00095                                     pset.getParameter<double>("r2MinF"),
00096                                     pset.getParameter<double>("r2MaxF"),
00097                                     pset.getParameter<double>("rMinI"),
00098                                     pset.getParameter<double>("rMaxI"),
00099                                     pset.getParameter<bool>("searchInTIDTEC"));
00100 
00101   theUpdator = new KFUpdator();
00102 }
00103 
00104 ElectronPixelSeedGenerator::~ElectronPixelSeedGenerator() {
00105 
00106   delete myMatchEle;
00107   delete myMatchPos;
00108   delete thePropagator;
00109   delete theUpdator;
00110 }
00111 
00112 void ElectronPixelSeedGenerator::setupES(const edm::EventSetup& setup) {
00113 
00114   // get records if necessary (called once per event)
00115   bool tochange=false;
00116 
00117   if (cacheIDMagField_!=setup.get<IdealMagneticFieldRecord>().cacheIdentifier()) {
00118     setup.get<IdealMagneticFieldRecord>().get(theMagField);
00119     cacheIDMagField_=setup.get<IdealMagneticFieldRecord>().cacheIdentifier();
00120     if (thePropagator) delete thePropagator;
00121     thePropagator = new PropagatorWithMaterial(alongMomentum,.000511,&(*theMagField)); 
00122     tochange=true;
00123   }
00124   if (cacheIDGeom_!=setup.get<TrackerRecoGeometryRecord>().cacheIdentifier()) {
00125     setup.get<TrackerRecoGeometryRecord>().get( theGeomSearchTracker );
00126     cacheIDGeom_=setup.get<TrackerRecoGeometryRecord>().cacheIdentifier();
00127   }
00128 
00129   if (cacheIDNavSchool_!=setup.get<NavigationSchoolRecord>().cacheIdentifier()) {
00130     edm::ESHandle<NavigationSchool> nav;
00131     setup.get<NavigationSchoolRecord>().get("SimpleNavigationSchool", nav);
00132     cacheIDNavSchool_=setup.get<NavigationSchoolRecord>().cacheIdentifier();
00133   
00134     theNavigationSchool = nav.product();
00135   }                                
00136 
00137   if (cacheIDCkfComp_!=setup.get<CkfComponentsRecord>().cacheIdentifier()) {
00138     edm::ESHandle<MeasurementTracker>    measurementTrackerHandle;
00139     setup.get<CkfComponentsRecord>().get(measurementTrackerHandle);
00140     cacheIDCkfComp_=setup.get<CkfComponentsRecord>().cacheIdentifier();
00141     theMeasurementTracker = measurementTrackerHandle.product();
00142     tochange=true;
00143   }
00144  
00145   edm::ESHandle<TrackerGeometry> trackerGeometryHandle;
00146   if (cacheIDTrkGeom_!=setup.get<TrackerDigiGeometryRecord>().cacheIdentifier()) {
00147     cacheIDTrkGeom_=setup.get<TrackerDigiGeometryRecord>().cacheIdentifier();
00148     setup.get<TrackerDigiGeometryRecord>().get(trackerGeometryHandle);
00149     tochange=true; //FIXME
00150   }
00151 
00152   if (tochange) {
00153     myMatchEle->setES(&(*theMagField),theMeasurementTracker,trackerGeometryHandle.product());
00154     myMatchPos->setES(&(*theMagField),theMeasurementTracker,trackerGeometryHandle.product());
00155   }
00156 
00157 }
00158 
00159 void  ElectronPixelSeedGenerator::run(edm::Event& e, const edm::EventSetup& setup, const reco::SuperClusterRefVector &sclRefs, TrajectorySeedCollection *seeds, reco::ElectronPixelSeedCollection & out){
00160 
00161   theInitialSeedColl=seeds;
00162 
00163   theSetup= &setup; 
00164   NavigationSetter theSetter(*theNavigationSchool);
00165 
00166   // get initial TrajectorySeeds if necessary
00167   //  if (fromTrackerSeeds_) e.getByLabel(initialSeeds_, theInitialSeedColl);
00168  
00169   // get the beamspot from the Event:
00170   e.getByType(theBeamSpot);
00171 
00172   // get its position
00173   double sigmaZ=theBeamSpot->sigmaZ();
00174   double sigmaZ0Error=theBeamSpot->sigmaZ0Error();
00175   double sq=sqrt(sigmaZ*sigmaZ+sigmaZ0Error*sigmaZ0Error);
00176   zmin1_=theBeamSpot->position().z()-3*sq;
00177   zmax1_=theBeamSpot->position().z()+3*sq;
00178 
00179   theMeasurementTracker->update(e); 
00180   
00181   for  (unsigned int i=0;i<sclRefs.size();++i) {
00182     // Find the seeds
00183     recHits_.clear();
00184 
00185     LogDebug ("run") << "new cluster, calling seedsFromThisCluster";
00186     seedsFromThisCluster(sclRefs[i],out);
00187   }
00188   
00189   LogDebug ("run") << ": For event "<<e.id();
00190   LogDebug ("run") <<"Nr of superclusters after filter: "<<sclRefs.size()
00191    <<", no. of ElectronPixelSeeds found  = " << out.size();
00192 }
00193 
00194 void ElectronPixelSeedGenerator::seedsFromThisCluster( edm::Ref<reco::SuperClusterCollection> seedCluster, reco::ElectronPixelSeedCollection& result)
00195 {
00196 
00197   float clusterEnergy = seedCluster->energy();
00198   GlobalPoint clusterPos(seedCluster->position().x(),
00199                          seedCluster->position().y(), 
00200                          seedCluster->position().z());
00201 
00202   const GlobalPoint 
00203    vertexPos(theBeamSpot->position().x(),theBeamSpot->position().y(),theBeamSpot->position().z());
00204   LogDebug("") << "[ElectronPixelSeedGenerator::seedsFromThisCluster] new supercluster with energy: " << clusterEnergy;
00205   LogDebug("") << "[ElectronPixelSeedGenerator::seedsFromThisCluster] and position: " << clusterPos;
00206 
00207   myMatchEle->set1stLayerZRange(zmin1_,zmax1_);
00208   myMatchPos->set1stLayerZRange(zmin1_,zmax1_);
00209   
00210   if (dynamicphiroad_)
00211     {
00212 
00213       float clusterEnergyT = clusterEnergy*sin(seedCluster->position().theta()) ;
00214 
00215       float deltaPhi1 = 0.875/clusterEnergyT + 0.055; 
00216       if (clusterEnergyT < lowPtThreshold_) deltaPhi1= deltaPhi1Low_;
00217       if (clusterEnergyT > highPtThreshold_) deltaPhi1= deltaPhi1High_;
00218 
00219       float ephimin1 = -deltaPhi1*sizeWindowENeg_ ;
00220       float ephimax1 =  deltaPhi1*(1.-sizeWindowENeg_);
00221       float pphimin1 = -deltaPhi1*(1.-sizeWindowENeg_);
00222       float pphimax1 =  deltaPhi1*sizeWindowENeg_;
00223 
00224       float phimin2  = -deltaPhi2_/2. ;
00225       float phimax2  =  deltaPhi2_/2. ;
00226 
00227       myMatchEle->set1stLayer(ephimin1,ephimax1);
00228       myMatchPos->set1stLayer(pphimin1,pphimax1);
00229       myMatchEle->set2ndLayer(phimin2,phimax2);
00230       myMatchPos->set2ndLayer(phimin2,phimax2);
00231 
00232     }
00233 
00234   PropagationDirection dir = alongMomentum;
00235   
00236    // try electron
00237   double aCharge=-1.;
00238  
00239   if (!fromTrackerSeeds_) {
00240     std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> > elePixelHits = 
00241       myMatchEle->compatibleHits(clusterPos,vertexPos, clusterEnergy, aCharge);
00242  
00243     float vertexZ = myMatchEle->getVertex();
00244     GlobalPoint eleVertex(theBeamSpot->position().x(),theBeamSpot->position().y(),vertexZ);
00245 
00246     if (!elePixelHits.empty() ) {
00247       LogDebug("ElectronPixelSeedGenerator") << "seedsFromThisCluster: electron compatible hits found ";
00248 
00249       std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> >::iterator v;
00250      
00251       for (v = elePixelHits.begin(); v != elePixelHits.end(); v++) {
00252         (*v).first.invert();
00253         bool valid = prepareElTrackSeed((*v).first.recHit(),(*v).second,eleVertex);
00254         if (valid) {
00255           reco::ElectronPixelSeed s(seedCluster,*pts_,recHits_,dir);
00256           result.push_back(s);
00257           delete pts_;
00258           pts_=0;
00259         }
00260       }
00261     } 
00262   } else {
00263     std::vector<TrajectorySeed> elePixelSeeds=
00264       myMatchEle->compatibleSeeds(theInitialSeedColl,clusterPos,vertexPos, clusterEnergy, aCharge);
00265       std::vector<TrajectorySeed>::iterator s;
00266       for (s = elePixelSeeds.begin(); s != elePixelSeeds.end(); s++) {
00267           reco::ElectronPixelSeed seed(seedCluster,*s);
00268           result.push_back(seed);
00269       }
00270   }
00271 
00272   // try positron
00273   aCharge=1.;  
00274   
00275   if (!fromTrackerSeeds_) {
00276     std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> > posPixelHits = 
00277       myMatchPos->compatibleHits(clusterPos,vertexPos, clusterEnergy, aCharge);
00278  
00279     float vertexZ = myMatchPos->getVertex();
00280     GlobalPoint posVertex(theBeamSpot->position().x(),theBeamSpot->position().y(),vertexZ);
00281 
00282     if (!posPixelHits.empty() ) {
00283       LogDebug("ElectronPixelSeedGenerator") << "seedsFromThisCluster: positron compatible hits found ";
00284 
00285       std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> >::iterator v;
00286       for (v = posPixelHits.begin(); v != posPixelHits.end(); v++) {
00287         bool valid = prepareElTrackSeed((*v).first.recHit(),(*v).second,posVertex);
00288         if (valid) {
00289           reco::ElectronPixelSeed s(seedCluster,*pts_,recHits_,dir);    
00290           result.push_back(s);
00291           delete pts_;
00292           pts_=0;
00293         }
00294       }
00295     }
00296   } else {
00297     std::vector<TrajectorySeed> posPixelSeeds=
00298       myMatchPos->compatibleSeeds(theInitialSeedColl,clusterPos,vertexPos, clusterEnergy, aCharge);
00299     std::vector<TrajectorySeed>::iterator s;
00300     for (s = posPixelSeeds.begin(); s != posPixelSeeds.end(); s++) {
00301       reco::ElectronPixelSeed seed(seedCluster,*s);
00302       result.push_back(seed);
00303     }
00304   }
00305   return ;
00306 }
00307 
00308 bool ElectronPixelSeedGenerator::prepareElTrackSeed(ConstRecHitPointer innerhit,
00309                                                     ConstRecHitPointer outerhit,
00310                                                     const GlobalPoint& vertexPos)
00311 {
00312   
00313   // debug prints
00314   LogDebug("") <<"[ElectronPixelSeedGenerator::prepareElTrackSeed] inner PixelHit   x,y,z "<<innerhit->globalPosition();
00315   LogDebug("") <<"[ElectronPixelSeedGenerator::prepareElTrackSeed] outer PixelHit   x,y,z "<<outerhit->globalPosition();
00316 
00317   pts_=0;
00318   recHits_.clear();
00319   
00320   SiPixelRecHit *pixhit=0;
00321   SiStripMatchedRecHit2D *striphit=0;
00322   const SiPixelRecHit* constpixhit = dynamic_cast <const SiPixelRecHit*> (innerhit->hit());
00323   if (constpixhit) {
00324     pixhit=new SiPixelRecHit(*constpixhit);
00325     recHits_.push_back(pixhit); 
00326   } else  return false;
00327   constpixhit =  dynamic_cast <const SiPixelRecHit *> (outerhit->hit());
00328   if (constpixhit) {
00329     pixhit=new SiPixelRecHit(*constpixhit);
00330     recHits_.push_back(pixhit); 
00331   } else {
00332     const SiStripMatchedRecHit2D * conststriphit=dynamic_cast <const SiStripMatchedRecHit2D *> (outerhit->hit());
00333     if (conststriphit) {
00334       striphit = new SiStripMatchedRecHit2D(*conststriphit);
00335       recHits_.push_back(striphit);   
00336     } else return false;
00337   }
00338 
00339   typedef TrajectoryStateOnSurface     TSOS;
00340   
00341   // make a spiral
00342   FastHelix helix(outerhit->globalPosition(),innerhit->globalPosition(),vertexPos,*theSetup);
00343   if ( !helix.isValid()) {
00344     return false;
00345   }
00346   FreeTrajectoryState fts = helix.stateAtVertex();
00347   TSOS propagatedState = thePropagator->propagate(fts,innerhit->det()->surface()) ;
00348   if (!propagatedState.isValid()) 
00349     return false;
00350   TSOS updatedState = theUpdator->update(propagatedState, *innerhit);
00351   
00352   TSOS propagatedState_out = thePropagator->propagate(updatedState,outerhit->det()->surface()) ;
00353   if (!propagatedState_out.isValid()) 
00354     return false;
00355   TSOS updatedState_out = theUpdator->update(propagatedState_out, *outerhit);
00356   // debug prints
00357   LogDebug("") <<"[ElectronPixelSeedGenerator::prepareElTrackSeed] final TSOS, position: "<<updatedState_out.globalPosition()<<" momentum: "<<updatedState_out.globalMomentum();
00358   LogDebug("") <<"[ElectronPixelSeedGenerator::prepareElTrackSeed] final TSOS Pt: "<<updatedState_out.globalMomentum().perp();
00359   pts_ =  transformer_.persistentState(updatedState_out, outerhit->geographicalId().rawId());
00360 
00361   return true;
00362 }