Reference Manual

Go to the documentation of this file.
00001 #ifndef RecoParticleFlow_PFProducer_PFBlockAlgo_h
00002 #define RecoParticleFlow_PFProducer_PFBlockAlgo_h 
00003 
00004 #include <set>
00005 #include <vector>
00006 #include <iostream>
00007 
00008 // #include "FWCore/Framework/interface/Handle.h"
00009 #include "DataFormats/Common/interface/Handle.h"
00010 // #include "FWCore/Framework/interface/OrphanHandle.h"
00011 #include "DataFormats/Common/interface/OrphanHandle.h"
00012 
00013 
00014 #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
00015 #include "DataFormats/ParticleFlowReco/interface/PFRecTrackFwd.h"
00016 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
00017 #include "DataFormats/ParticleFlowReco/interface/PFRecTrack.h"
00018 #include "DataFormats/ParticleFlowReco/interface/PFDisplacedTrackerVertex.h" // gouzevitch
00019 #include "DataFormats/ParticleFlowReco/interface/PFConversionFwd.h"
00020 #include "DataFormats/ParticleFlowReco/interface/PFConversion.h"
00021 #include "DataFormats/ParticleFlowReco/interface/PFV0Fwd.h"
00022 #include "DataFormats/ParticleFlowReco/interface/PFV0.h"
00023 #include "DataFormats/ParticleFlowReco/interface/GsfPFRecTrack.h"
00024 #include "DataFormats/ParticleFlowReco/interface/GsfPFRecTrackFwd.h"
00025 #include "DataFormats/ParticleFlowReco/interface/PFBrem.h"
00026 #include "DataFormats/ParticleFlowReco/interface/PFTrajectoryPoint.h"
00027 
00028 #include "DataFormats/ParticleFlowReco/interface/PFBlockElement.h"
00029 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
00030 #include "DataFormats/ParticleFlowReco/interface/PFBlockFwd.h"
00031 
00032 // Glowinski & Gouzevitch
00033 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"             
00034 #include "RecoParticleFlow/PFProducer/interface/KDTreeLinkerTrackEcal.h" 
00035 #include "RecoParticleFlow/PFProducer/interface/KDTreeLinkerTrackHcal.h" 
00036 #include "RecoParticleFlow/PFProducer/interface/KDTreeLinkerPSEcal.h" 
00037 // !Glowinski & Gouzevitch
00038 
00039 // #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
00040 
00041 #include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h"
00042 #include "RecoParticleFlow/PFProducer/interface/PhotonSelectorAlgo.h"
00043 #include "RecoParticleFlow/PFProducer/interface/PFBlockElementSCEqual.h"
00044 #include "RecoParticleFlow/PFClusterTools/interface/PFResolutionMap.h"
00045 
00046 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00047 #include "DataFormats/MuonReco/interface/Muon.h"
00048 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
00049 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00050 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
00051 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
00052 #include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
00053 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
00054 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
00055 #include "DataFormats/EgammaCandidates/interface/Photon.h"
00056 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementSuperCluster.h"
00057 
00058 #include "RecoParticleFlow/PFClusterTools/interface/ClusterClusterMapping.h"
00059 
00060 #include "RecoParticleFlow/PFProducer/interface/PFBlockLink.h"
00061 
00062 #include <map>
00063 
00064 
00066 
00071 class PFBlockAlgo {
00072 
00073  public:
00074 
00075   PFBlockAlgo();
00076 
00077   ~PFBlockAlgo();
00078   
00079 
00080   void setParameters( std::vector<double>& DPtovPtCut, 
00081                       std::vector<unsigned>& NHitCut,
00082                       bool useConvBremPFRecTracks,
00083                       bool useIterTracking,
00084                       int nuclearInteractionsPurity,
00085                       bool useEGPhotons,
00086                       std::vector<double> & photonSelectionCuts
00087                       );
00088   
00089   // Glowinski & Gouzevitch
00090   void setUseOptimization(bool useKDTreeTrackEcalLinker);
00091   // ! Glowinski & Gouzevitch
00092 
00093   typedef std::vector<bool> Mask;
00094 
00096   template< template<typename> class T>
00097     void  setInput(const T<reco::PFRecTrackCollection>&    trackh,
00098                    const T<reco::GsfPFRecTrackCollection>&    gsftrackh,
00099                    const T<reco::GsfPFRecTrackCollection>&    convbremgsftrackh,
00100                    const T<reco::MuonCollection>&    muonh,
00101                    const T<reco::PFDisplacedTrackerVertexCollection>&  nuclearh,
00102                    const T<reco::PFRecTrackCollection>&    nucleartrackh,
00103                    const T<reco::PFConversionCollection>&  conv,
00104                    const T<reco::PFV0Collection>&  v0,
00105                    const T<reco::PFClusterCollection>&  ecalh,
00106                    const T<reco::PFClusterCollection>&  hcalh,
00107                    const T<reco::PFClusterCollection>&  hfemh,
00108                    const T<reco::PFClusterCollection>&  hfhadh,
00109                    const T<reco::PFClusterCollection>&  psh,
00110                    const T<reco::PhotonCollection>&  egphh,
00111                    const Mask& trackMask = dummyMask_,
00112                    const Mask& gsftrackMask = dummyMask_,
00113                    const Mask& ecalMask = dummyMask_,
00114                    const Mask& hcalMask = dummyMask_,
00115                    const Mask& hfemMask = dummyMask_,              
00116                    const Mask& hfhadMask = dummyMask_,
00117                    const Mask& psMask = dummyMask_,
00118                    const Mask& phMask = dummyMask_); 
00119   
00121   template< template<typename> class T >
00122     void setInput(const T<reco::PFRecTrackCollection>&    trackh,
00123                   const T<reco::MuonCollection>&    muonh,
00124                   const T<reco::PFClusterCollection>&  ecalh,
00125                   const T<reco::PFClusterCollection>&  hcalh,
00126                   const T<reco::PFClusterCollection>&  hfemh,
00127                   const T<reco::PFClusterCollection>&  hfhadh,
00128                   const T<reco::PFClusterCollection>&  psh,
00129                   const Mask& trackMask = dummyMask_,
00130                   const Mask& ecalMask = dummyMask_,
00131                   const Mask& hcalMask = dummyMask_,
00132                   const Mask& psMask = dummyMask_ ) {
00133     T<reco::GsfPFRecTrackCollection> gsftrackh;
00134     T<reco::GsfPFRecTrackCollection> convbremgsftrackh;
00135     //T<reco::MuonCollection> muonh;
00136     T<reco::PFDisplacedTrackerVertexCollection> nuclearh;
00137     T<reco::PFRecTrackCollection>    nucleartrackh;
00138     T<reco::PFConversionCollection> convh;
00139     T<reco::PFV0Collection> v0;
00140     T<reco::PhotonCollection> phh;
00141     setInput<T>( trackh, gsftrackh, convbremgsftrackh, muonh, nuclearh, nucleartrackh, convh, v0, 
00142                  ecalh, hcalh, hfemh, hfhadh, psh, phh,
00143                  trackMask, ecalMask, hcalMask, psMask); 
00144   }
00145   
00147   template< template<typename> class T >
00148     void setInput(const T<reco::PFRecTrackCollection>&    trackh,
00149                   const T<reco::GsfPFRecTrackCollection>&    gsftrackh,
00150                   const T<reco::PFClusterCollection>&  ecalh,
00151                   const T<reco::PFClusterCollection>&  hcalh,
00152                   const T<reco::PFClusterCollection>&  psh,
00153                   const Mask& trackMask = dummyMask_,
00154                   const Mask& gsftrackMask = dummyMask_,
00155                   const Mask& ecalMask = dummyMask_,
00156                   const Mask& hcalMask = dummyMask_,
00157                   const Mask& psMask = dummyMask_ ) {
00158     T<reco::GsfPFRecTrackCollection> convbremgsftrackh;
00159     T<reco::MuonCollection> muonh;
00160     T<reco::PFDisplacedTrackerVertexCollection>  nuclearh;
00161     T<reco::PFRecTrackCollection>    nucleartrackh;
00162     T<reco::PFConversionCollection> convh;
00163     T<reco::PFV0Collection> v0;
00164     T<reco::PhotonCollection> egphh;
00165     setInput<T>( trackh, gsftrackh, convbremgsftrackh, muonh, nuclearh, nucleartrackh, convh, v0, ecalh, hcalh, psh, egphh,
00166                  trackMask, gsftrackMask,ecalMask, hcalMask, psMask); 
00167   }
00168   
00169   
00171   void setDebug( bool debug ) {debug_ = debug;}
00172   
00174   void findBlocks();
00175   
00176 
00178   /*   const  reco::PFBlockCollection& blocks() const {return *blocks_;} */
00179   const std::auto_ptr< reco::PFBlockCollection >& blocks() const 
00180     {return blocks_;}
00181   
00183   std::auto_ptr< reco::PFBlockCollection > transferBlocks() {return blocks_;}
00184 
00186   typedef std::list< reco::PFBlockElement* >::iterator IE;
00187   typedef std::list< reco::PFBlockElement* >::const_iterator IEC;  
00188   typedef reco::PFBlockCollection::const_iterator IBC;
00189   
00190  private:
00191   
00198   IE associate(IE next, IE last, std::vector<PFBlockLink>& links);
00199 
00202   void packLinks(reco::PFBlock& block, 
00203                  const std::vector<PFBlockLink>& links) const; 
00204 
00206   void checkDisplacedVertexLinks( reco::PFBlock& block ) const;
00207   
00211   void buildGraph(); 
00212 
00214   inline bool linkPrefilter(const reco::PFBlockElement* last, const reco::PFBlockElement* next) const;
00215 
00217   void link( const reco::PFBlockElement* el1, 
00218              const reco::PFBlockElement* el2, 
00219              PFBlockLink::Type& linktype, 
00220              reco::PFBlock::LinkTest& linktest,
00221              double& dist) const;
00222                  
00225   double testTrackAndPS(const reco::PFRecTrack& track,
00226                         const reco::PFCluster& ps) const;
00227   
00230   double testECALAndHCAL(const reco::PFCluster& ecal, 
00231                          const reco::PFCluster& hcal) const;
00232                          
00235   double testPS1AndPS2(const reco::PFCluster& ps1,
00236                        const reco::PFCluster& ps2) const;
00237 
00239   double testLinkBySuperCluster(const reco::PFClusterRef & elt1,
00240                                 const reco::PFClusterRef & elt2) const;   
00241 
00243   double testSuperClusterPFCluster(const reco::SuperClusterRef & sct1,
00244                                    const reco::PFClusterRef & elt2) const;
00245 
00249   void checkMaskSize( const reco::PFRecTrackCollection& tracks,
00250                       const reco::GsfPFRecTrackCollection& gsftracks,
00251                       const reco::PFClusterCollection&  ecals,
00252                       const reco::PFClusterCollection&  hcals,
00253                       const reco::PFClusterCollection&  hfems,
00254                       const reco::PFClusterCollection&  hfhads,
00255                       const reco::PFClusterCollection&  pss, 
00256                       const reco::PhotonCollection&  egphh, 
00257                       const Mask& trackMask,
00258                       const Mask& gsftrackMask, 
00259                       const Mask& ecalMask, 
00260                       const Mask& hcalMask,
00261                       const Mask& hfemMask,
00262                       const Mask& hfhadMask,                  
00263                       const Mask& psMask,
00264                       const Mask& phMask) const;
00265 
00267   // PFResolutionMap* openResolutionMap(const char* resMapName);
00268 
00270   bool goodPtResolution( const reco::TrackRef& trackref);
00271 
00272   double testLinkByVertex(const reco::PFBlockElement* elt1,
00273                           const reco::PFBlockElement* elt2) const;
00274 
00275   int muAssocToTrack( const reco::TrackRef& trackref,
00276                       const edm::Handle<reco::MuonCollection>& muonh) const;
00277   int muAssocToTrack( const reco::TrackRef& trackref,
00278                       const edm::OrphanHandle<reco::MuonCollection>& muonh) const;
00279 
00280   template< template<typename> class T>
00281     void fillFromPhoton(const T<reco::PhotonCollection> &, unsigned isc, reco::PFBlockElementSuperCluster * pfbe);
00282 
00283   std::auto_ptr< reco::PFBlockCollection >    blocks_;
00284   
00286   // std::vector< reco::PFCandidate >   particles_;
00287 
00288   // the test elements will be transferred to the blocks
00289   std::list< reco::PFBlockElement* >     elements_;
00290 
00291   // Glowinski & Gouzevitch
00292   bool useKDTreeTrackEcalLinker_;
00293   KDTreeLinkerTrackEcal TELinker_;
00294   KDTreeLinkerTrackHcal THLinker_;
00295   KDTreeLinkerPSEcal    PSELinker_;
00296   // !Glowinski & Gouzevitch
00297 
00298   static const Mask                      dummyMask_;
00299 
00301   std::vector<double> DPtovPtCut_;
00302   
00304   std::vector<unsigned> NHitCut_;
00305   
00307   bool useIterTracking_;
00308 
00310   bool useEGPhotons_;
00311 
00312   // This parameters defines the level of purity of
00313   // nuclear interactions choosen.
00314   // Level 1 is only high Purity sample labeled as isNucl.
00315   // Level 2 isNucl + isNucl_Loose (2 secondary tracks vertices)
00316   // Level 3 isNucl + isNucl_Loose + isNucl_Kink
00317   //         (low purity sample made of 1 primary and 1 secondary track)
00318   // By default the level 1 is teh safest one.
00319   int nuclearInteractionsPurity_;
00320 
00322   bool  useConvBremPFRecTracks_;
00323 
00325   const PhotonSelectorAlgo * photonSelector_;
00327   std::vector<reco::SuperClusterRef > superClusters_;
00328 
00330   //  std::map<reco::PFClusterRef,int>  pfcRefSCMap_;
00331   std::vector<int> pfcSCVec_;
00332 
00333   // A boolean to avoid to compare ECAL and ECAl if there i no superclusters in the event
00334   bool bNoSuperclus_;
00335 
00337   std::vector<std::vector<reco::PFClusterRef> > scpfcRefs_;
00339   bool   debug_;
00340   
00341   friend std::ostream& operator<<(std::ostream&, const PFBlockAlgo&);
00342 
00343 };
00344 
00345 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementGsfTrack.h"
00346 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementBrem.h"
00347 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementTrack.h"
00348 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementCluster.h"
00349 #include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
00350 
00351 template< template<typename> class T >
00352 void
00353 PFBlockAlgo::setInput(const T<reco::PFRecTrackCollection>&    trackh,
00354                       const T<reco::GsfPFRecTrackCollection>&    gsftrackh, 
00355                       const T<reco::GsfPFRecTrackCollection>&    convbremgsftrackh,
00356                       const T<reco::MuonCollection>&    muonh,
00357                       const T<reco::PFDisplacedTrackerVertexCollection>&  nuclearh,
00358                       const T<reco::PFRecTrackCollection>&    nucleartrackh,
00359                       const T<reco::PFConversionCollection>&  convh,
00360                       const T<reco::PFV0Collection>&  v0,
00361                       const T<reco::PFClusterCollection>&  ecalh,
00362                       const T<reco::PFClusterCollection>&  hcalh,
00363                       const T<reco::PFClusterCollection>&  hfemh,
00364                       const T<reco::PFClusterCollection>&  hfhadh,
00365                       const T<reco::PFClusterCollection>&  psh,
00366                       const T<reco::PhotonCollection>& egphh,
00367                       const Mask& trackMask,
00368                       const Mask& gsftrackMask,
00369                       const Mask& ecalMask,
00370                       const Mask& hcalMask,
00371                       const Mask& hfemMask,
00372                       const Mask& hfhadMask,
00373                       const Mask& psMask,
00374                       const Mask& phMask) {
00375 
00376 
00377   checkMaskSize( *trackh,
00378                  *gsftrackh,
00379                  *ecalh,
00380                  *hcalh,
00381                  *hfemh,
00382                  *hfhadh,
00383                  *psh,
00384                  *egphh,
00385                  trackMask,
00386                  gsftrackMask,
00387                  ecalMask,
00388                  hcalMask,
00389                  hfemMask,
00390                  hfhadMask,
00391                  psMask,
00392                  phMask);
00393 
00394   /*
00395   if (nucleartrackh.isValid()){
00396     for(unsigned i=0;i<nucleartrackh->size(); i++) {
00397       reco::PFRecTrackRef trackRef(nucleartrackh,i);
00398       std::cout << *trackRef << std::endl;
00399     }
00400   }
00401   */
00402 
00404   std::vector<reco::PFRecTrackRef> convBremPFRecTracks;
00405   convBremPFRecTracks.clear();
00406   // Super cluster mapping
00407   superClusters_.clear();
00408   scpfcRefs_.clear();
00409   pfcSCVec_.clear();
00410 
00411   if(gsftrackh.isValid() ) {
00412     const  reco::GsfPFRecTrackCollection PFGsfProd = *(gsftrackh.product());
00413     for(unsigned i=0;i<gsftrackh->size(); ++i) {
00414       if( !gsftrackMask.empty() &&
00415           !gsftrackMask[i] ) continue;
00416       reco::GsfPFRecTrackRef refgsf(gsftrackh,i );   
00417    
00418       if((refgsf).isNull()) continue;
00419       reco::GsfTrackRef gsf=refgsf->gsfTrackRef();
00420 
00421       // retrieve and save the SC if ECAL-driven - Florian
00422       if(gsf->extra().isAvailable() && gsf->extra()->seedRef().isAvailable()) {
00423         reco::ElectronSeedRef seedRef=  gsf->extra()->seedRef().castTo<reco::ElectronSeedRef>();
00424         // check that the seed is valid
00425         if(seedRef.isAvailable() && seedRef->isEcalDriven()) {
00426           reco::SuperClusterRef scRef = seedRef->caloCluster().castTo<reco::SuperClusterRef>();
00427           if(scRef.isNonnull())   {          
00428             std::vector<reco::SuperClusterRef>::iterator itcheck=find(superClusters_.begin(),superClusters_.end(),scRef);
00429             // not present, add it
00430             if(itcheck==superClusters_.end())
00431               {
00432                 superClusters_.push_back(scRef);
00433                 reco::PFBlockElementSuperCluster * sce =
00434                   new reco::PFBlockElementSuperCluster(scRef);
00435                 sce->setFromGsfElectron(true);
00436                 elements_.push_back(sce);
00437               }
00438             else // it is already present, update the PFBE
00439               {
00440                 PFBlockElementSCEqual myEqual(scRef);
00441                 std::list<reco::PFBlockElement*>::iterator itcheck=find_if(elements_.begin(),elements_.end(),myEqual);
00442                 if(itcheck!=elements_.end())
00443                   {
00444                     reco::PFBlockElementSuperCluster* thePFBE=dynamic_cast<reco::PFBlockElementSuperCluster*>(*itcheck);
00445                     thePFBE->setFromGsfElectron(true);
00446                     //              std::cout << " Updating element to add electron information" << std::endl;
00447                   }
00448 //              else
00449 //                {
00450 //                  std::cout << " Missing element " << std::endl;
00451 //                }
00452               }
00453           }
00454         }
00455       }
00456 
00457       reco::PFBlockElement* gsfEl;
00458       
00459       const  std::vector<reco::PFTrajectoryPoint> 
00460         PfGsfPoint =  PFGsfProd[i].trajectoryPoints();
00461   
00462       unsigned int c_gsf=0;
00463       bool PassTracker = false;
00464       bool GetPout = false;
00465       unsigned int IndexPout = 0;
00466       
00467       typedef std::vector<reco::PFTrajectoryPoint>::const_iterator IP;
00468       for(IP itPfGsfPoint =  PfGsfPoint.begin();  
00469           itPfGsfPoint!= PfGsfPoint.end();++itPfGsfPoint) {
00470         
00471         if (itPfGsfPoint->isValid()){
00472           int layGsfP = itPfGsfPoint->layer();
00473           if (layGsfP == -1) PassTracker = true;
00474           if (PassTracker && layGsfP > 0 && GetPout == false) {
00475             IndexPout = c_gsf-1;
00476             GetPout = true;
00477           }
00478           //const math::XYZTLorentzVector GsfMoment = itPfGsfPoint->momentum();
00479           ++c_gsf;
00480         }
00481       }
00482       math::XYZTLorentzVector pin = PfGsfPoint[0].momentum();      
00483       math::XYZTLorentzVector pout = PfGsfPoint[IndexPout].momentum();
00484 
00486       if(useConvBremPFRecTracks_) {
00487         const std::vector<reco::PFRecTrackRef>& temp_convBremPFRecTracks(refgsf->convBremPFRecTrackRef());
00488         if(temp_convBremPFRecTracks.size() > 0) {
00489           for(unsigned int iconv = 0; iconv <temp_convBremPFRecTracks.size(); ++iconv) {
00490             convBremPFRecTracks.push_back(temp_convBremPFRecTracks[iconv]);
00491           }
00492         }
00493       }
00494       
00495       gsfEl = new reco::PFBlockElementGsfTrack(refgsf, pin, pout);
00496       
00497       elements_.push_back( gsfEl);
00498 
00499       std::vector<reco::PFBrem> pfbrem = refgsf->PFRecBrem();
00500       
00501       for (unsigned i2=0;i2<pfbrem.size(); ++i2) {
00502         const double DP = pfbrem[i2].DeltaP();
00503         const double SigmaDP =  pfbrem[i2].SigmaDeltaP(); 
00504         const unsigned int TrajP = pfbrem[i2].indTrajPoint();
00505         if(TrajP == 99) continue;
00506 
00507         reco::PFBlockElement* bremEl;
00508         bremEl = new reco::PFBlockElementBrem(refgsf,DP,SigmaDP,TrajP);
00509         elements_.push_back(bremEl);
00510         
00511       }
00512     }
00513 
00514   }
00515 
00517   if(useEGPhotons_ && egphh.isValid()) {
00518     unsigned size=egphh->size();
00519     for(unsigned isc=0; isc<size; ++isc) {
00520       if(!phMask.empty() && !(phMask)[isc] ) continue;
00521       if(!photonSelector_->passPhotonSelection((*egphh)[isc])) continue;
00522       //      std::cout << " Selected a supercluster" << std::endl;
00523       // Add only the super clusters not already included 
00524       reco::SuperClusterRef scRef((*egphh)[isc].superCluster());
00525       std::vector<reco::SuperClusterRef>::iterator itcheck=find(superClusters_.begin(),superClusters_.end(),(*egphh)[isc].superCluster());
00526       if(itcheck==superClusters_.end())
00527         {
00528           superClusters_.push_back(scRef);
00529           reco::PFBlockElementSuperCluster * sce =
00530             new reco::PFBlockElementSuperCluster((*egphh)[isc].superCluster());
00531           fillFromPhoton(egphh,isc,sce);
00532           elements_.push_back(sce);
00533         }
00534       else
00535         {
00536           PFBlockElementSCEqual myEqual(scRef);
00537           std::list<reco::PFBlockElement*>::iterator itcheck=find_if(elements_.begin(),elements_.end(),myEqual);
00538           if(itcheck!=elements_.end())
00539             {
00540               reco::PFBlockElementSuperCluster* thePFBE=dynamic_cast<reco::PFBlockElementSuperCluster*>(*itcheck);
00541               fillFromPhoton(egphh,isc,thePFBE);
00542               thePFBE->setFromPhoton(true);
00543               //              std::cout << " Updating element to add Photon information " << photonSelector_->passPhotonSelection((*egphh)[isc]) << std::endl;
00544 
00545             }
00546 //        else
00547 //          {
00548 //            std::cout << " Missing element " << std::endl;
00549 //          }
00550         }
00551     }
00552   }
00553   
00554   // set the vector to the right size so to allow random access
00555   scpfcRefs_.resize(superClusters_.size());
00556 
00558 
00560 
00561   if(convh.isValid() ) {
00562     reco::PFBlockElement* trkFromConversionElement;
00563     for(unsigned i=0;i<convh->size(); ++i) {
00564       reco::PFConversionRef convRef(convh,i);
00565 
00566       unsigned int trackSize=(convRef->pfTracks()).size();
00567       if ( convRef->pfTracks().size() < 2) continue;
00568       for(unsigned iTk=0;iTk<trackSize; ++iTk) {
00569         
00570         reco::PFRecTrackRef compPFTkRef = convRef->pfTracks()[iTk];     
00571         trkFromConversionElement = new reco::PFBlockElementTrack(convRef->pfTracks()[iTk]);
00572         trkFromConversionElement->setConversionRef( convRef->originalConversion(), reco::PFBlockElement::T_FROM_GAMMACONV);
00573 
00574         elements_.push_back( trkFromConversionElement );
00575 
00576 
00577         if (debug_){
00578           std::cout << "PF Block Element from Conversion electron " << 
00579             (*trkFromConversionElement).trackRef().key() << std::endl;
00580           std::cout << *trkFromConversionElement << std::endl;
00581         }
00582        
00583       }     
00584     }  
00585   }
00586   
00587   
00589   
00591   
00592   if(v0.isValid() ) {
00593     reco::PFBlockElement* trkFromV0Element = 0;
00594     for(unsigned i=0;i<v0->size(); ++i) {
00595       reco::PFV0Ref v0Ref( v0, i );
00596       unsigned int trackSize=(v0Ref->pfTracks()).size();
00597       for(unsigned iTk=0;iTk<trackSize; ++iTk) {
00598 
00599         reco::PFRecTrackRef newPFRecTrackRef = (v0Ref->pfTracks())[iTk]; 
00600         reco::TrackBaseRef newTrackBaseRef(newPFRecTrackRef->trackRef());
00601         bool bNew = true;
00602         
00605         for(IE iel = elements_.begin(); iel != elements_.end(); ++iel){
00606           reco::TrackBaseRef elemTrackBaseRef((*iel)->trackRef());
00607           if (newTrackBaseRef == elemTrackBaseRef){         
00608             trkFromV0Element = *iel;
00609             bNew = false;
00610             continue;
00611           }
00612         } 
00613 
00615         if (bNew) {
00616           trkFromV0Element = new reco::PFBlockElementTrack(v0Ref->pfTracks()[iTk]);
00617           elements_.push_back( trkFromV0Element );
00618         }
00619 
00620         trkFromV0Element->setV0Ref( v0Ref->originalV0(),
00621                                     reco::PFBlockElement::T_FROM_V0 );
00622           
00623         if (debug_){
00624           std::cout << "PF Block Element from V0 track New = " << bNew 
00625                     << (*trkFromV0Element).trackRef().key() << std::endl;
00626           std::cout << *trkFromV0Element << std::endl;
00627         }
00628 
00629         
00630       }
00631     }
00632   }
00633   
00635 
00637 
00638   if(nuclearh.isValid()) {
00639     reco::PFBlockElement* trkFromDisplacedVertexElement = 0;
00640     for(unsigned i=0;i<nuclearh->size(); ++i) {
00641 
00642       const reco::PFDisplacedTrackerVertexRef dispacedVertexRef( nuclearh, i );
00643 
00644       //      std::cout << "Nuclear Interactions Purity " <<  nuclearInteractionsPurity_ << std::endl;
00645       //     dispacedVertexRef->displacedVertexRef()->Dump();
00646       //bool bIncludeVertices = true;
00647       // We add a cut at rho > 2.7 since this corresponds to the lower edge of the beam pipe
00648       // This cut have to be changer when a new beam pipe would be installed
00649       
00650       bool bIncludeVertices = false; 
00651       bool bNucl = dispacedVertexRef->displacedVertexRef()->isNucl()
00652         && dispacedVertexRef->displacedVertexRef()->position().rho()>  2.7;
00653       bool bNucl_Loose = dispacedVertexRef->displacedVertexRef()->isNucl_Loose();
00654       bool bNucl_Kink = dispacedVertexRef->displacedVertexRef()->isNucl_Kink();
00655 
00656       if (nuclearInteractionsPurity_ >= 1) bIncludeVertices = bNucl;
00657       if (nuclearInteractionsPurity_ >= 2) bIncludeVertices = bIncludeVertices || bNucl_Loose;
00658       if (nuclearInteractionsPurity_ >= 3) bIncludeVertices = bIncludeVertices || bNucl_Kink;
00659 
00660       if (bIncludeVertices){
00661 
00662         unsigned int trackSize= dispacedVertexRef->pfRecTracks().size();
00663         if (debug_){
00664           std::cout << "" << std::endl;
00665           std::cout << "Displaced Vertex " << i << std::endl;
00666           dispacedVertexRef->displacedVertexRef()->Dump();
00667         }
00668         for(unsigned iTk=0;iTk < trackSize; ++iTk) {
00669 
00670 
00671           // This peace of code looks weired at first but it seems to be necessary to let 
00672           // PFRooTEvent work properly. Since the track position called REPPoint is transient 
00673           // it has to be calculated and the PFRecTrack collection associted to Displaced Vertex is not
00674           // anymore the original collection. So here we match both collections if possible
00675           reco::PFRecTrackRef newPFRecTrackRef = dispacedVertexRef->pfRecTracks()[iTk]; 
00676           reco::TrackBaseRef constTrackBaseRef(newPFRecTrackRef->trackRef());
00677 
00678           
00679           if (nucleartrackh.isValid()){
00680             for(unsigned i=0;i<nucleartrackh->size(); ++i) {
00681               reco::PFRecTrackRef transientPFRecTrackRef(nucleartrackh,i);
00682               reco::TrackBaseRef transientTrackBaseRef(transientPFRecTrackRef->trackRef());
00683               if (constTrackBaseRef==transientTrackBaseRef){
00684                 newPFRecTrackRef = transientPFRecTrackRef;
00685                 break;
00686               }
00687             }
00688           }
00689           reco::TrackBaseRef newTrackBaseRef(newPFRecTrackRef->trackRef());
00690           
00691 
00692 
00693           bool bNew = true;
00694           reco::PFBlockElement::TrackType blockType;
00695 
00698           for(IE iel = elements_.begin(); iel != elements_.end(); ++iel){
00699             reco::TrackBaseRef elemTrackBaseRef((*iel)->trackRef());
00700             if (newTrackBaseRef == elemTrackBaseRef){
00701               trkFromDisplacedVertexElement = *iel;
00702               bNew = false;
00703               continue;
00704             }
00705           }
00706 
00707 
00709           if (bNew) { 
00710             
00711 
00712 
00713             trkFromDisplacedVertexElement = new reco::PFBlockElementTrack(newPFRecTrackRef);
00714             elements_.push_back( trkFromDisplacedVertexElement );
00715           }
00716 
00717           if (dispacedVertexRef->isIncomingTrack(newPFRecTrackRef)) 
00718             blockType = reco::PFBlockElement::T_TO_DISP;
00719           else if (dispacedVertexRef->isOutgoingTrack(newPFRecTrackRef)) 
00720             blockType = reco::PFBlockElement::T_FROM_DISP;
00721           else 
00722             blockType = reco::PFBlockElement::DEFAULT;
00723 
00725           trkFromDisplacedVertexElement->setDisplacedVertexRef( dispacedVertexRef, blockType );
00726 
00727 
00728           if (debug_){
00729             std::cout << "PF Block Element from DisplacedTrackingVertex track New = " << bNew
00730                       << (*trkFromDisplacedVertexElement).trackRef().key() << std::endl;
00731             std::cout << *trkFromDisplacedVertexElement << std::endl;
00732           }
00733         
00734         
00735         }
00736       }
00737     }  
00738 
00739     if (debug_) std::cout << "" << std::endl;
00740 
00741   }
00742 
00744 
00748 
00749   if(trackh.isValid() ) {
00750 
00751     if (debug_) std::cout << "Tracks already in from Displaced Vertices " << std::endl;
00752 
00753     Mask trackMaskVertex;
00754 
00755     for(unsigned i=0;i<trackh->size(); ++i) {
00756       reco::PFRecTrackRef pfRefTrack( trackh,i );
00757       reco::TrackRef trackRef = pfRefTrack->trackRef();
00758 
00759       bool bMask = true;
00760       for(IE iel = elements_.begin(); iel != elements_.end(); ++iel){
00761         reco::TrackRef elemTrackRef = (*iel)->trackRef();
00762         if( trackRef == elemTrackRef ) {
00763           if (debug_) std::cout << " " << trackRef.key();
00764           bMask = false; continue;
00765         }
00766       }
00767     
00768       trackMaskVertex.push_back(bMask);
00769     }
00770 
00771     if (debug_) std::cout << "" << std::endl;
00772 
00773     if (debug_) std::cout << "Additionnal tracks from main collection " << std::endl;
00774 
00775     for(unsigned i=0;i<trackh->size(); ++i) {
00776 
00777 
00778       // this track has been disabled
00779       if( !trackMask.empty() && !trackMask[i] ) continue;
00780       
00781       reco::PFRecTrackRef ref( trackh,i );
00782 
00783       if (debug_) std::cout << " " << ref->trackRef().key();
00784 
00785       // Get the eventual muon associated to this track
00786       int muId_ = muAssocToTrack( ref->trackRef(), muonh );
00787       bool thisIsAPotentialMuon = false;
00788       if( muId_ != -1 ) {
00789         reco::MuonRef muonref( muonh, muId_ );
00790         thisIsAPotentialMuon = 
00791           PFMuonAlgo::isLooseMuon(muonref) || 
00792           PFMuonAlgo::isMuon(muonref);
00793       }
00794       // Reject bad tracks (except if identified as muon
00795       if( !thisIsAPotentialMuon && !goodPtResolution( ref->trackRef() ) ) continue;
00796 
00797       if (thisIsAPotentialMuon && debug_) std::cout << "Potential Muon P " <<  ref->trackRef()->p() 
00798                                                     << " pt " << ref->trackRef()->p() << std::endl; 
00799 
00800 
00801 
00802       reco::PFBlockElement* primaryElement = new reco::PFBlockElementTrack( ref );
00803 
00804       if( muId_ != -1 ) {
00805         // if a muon has been found
00806         reco::MuonRef muonref( muonh, muId_ );
00807 
00808         // If this track was already added to the collection, we just need to find the associated element and 
00809         // attach to it the reference
00810         if (!trackMaskVertex.empty() && !trackMaskVertex[i]){
00811           reco::TrackRef primaryTrackRef = ref->trackRef();
00812           for(IE iel = elements_.begin(); iel != elements_.end(); ++iel){
00813             reco::TrackRef elemTrackRef = (*iel)->trackRef();
00814             if( primaryTrackRef == elemTrackRef ) {
00815               (*iel)->setMuonRef( muonref );
00816               if (debug_) std::cout << "One of the tracks identified in displaced vertices collections was spotted as muon" <<std:: endl;
00817             }
00818           }
00819         } else primaryElement->setMuonRef( muonref );
00820 
00821       }
00822 
00823       if (!trackMaskVertex.empty() && !trackMaskVertex[i]) continue;
00824 
00825       
00826       // set track type T_FROM_GAMMA for pfrectracks associated to conv brems
00827       if(useConvBremPFRecTracks_) {
00828         if(convBremPFRecTracks.size() > 0.) {
00829           for(unsigned int iconv = 0; iconv < convBremPFRecTracks.size(); ++iconv) {
00830             if((*ref).trackRef() == (*convBremPFRecTracks[iconv]).trackRef()) {
00831               bool value = true;
00832               primaryElement->setTrackType(reco::PFBlockElement::T_FROM_GAMMACONV, value);
00833             }
00834           }
00835         }
00836       }
00837       elements_.push_back( primaryElement );
00838     }
00839 
00840     if (debug_) std::cout << " " << std::endl;
00841 
00842   }
00843 
00844  
00845   // -------------- GSF tracks and brems for Conversion Recovery ----------
00846    
00847   if(convbremgsftrackh.isValid() ) {
00848     
00849  
00850     const  reco::GsfPFRecTrackCollection ConvPFGsfProd = *(convbremgsftrackh.product());
00851     for(unsigned i=0;i<convbremgsftrackh->size(); ++i) {
00852 
00853       reco::GsfPFRecTrackRef refgsf(convbremgsftrackh,i );   
00854       
00855       if((refgsf).isNull()) continue;
00856       
00857       reco::PFBlockElement* gsfEl;
00858       
00859       const  std::vector<reco::PFTrajectoryPoint> 
00860         PfGsfPoint =  ConvPFGsfProd[i].trajectoryPoints();
00861       
00862       unsigned int c_gsf=0;
00863       bool PassTracker = false;
00864       bool GetPout = false;
00865       unsigned int IndexPout = -1;
00866       
00867       typedef std::vector<reco::PFTrajectoryPoint>::const_iterator IP;
00868       for(IP itPfGsfPoint =  PfGsfPoint.begin();  
00869           itPfGsfPoint!= PfGsfPoint.end();++itPfGsfPoint) {
00870         
00871         if (itPfGsfPoint->isValid()){
00872           int layGsfP = itPfGsfPoint->layer();
00873           if (layGsfP == -1) PassTracker = true;
00874           if (PassTracker && layGsfP > 0 && GetPout == false) {
00875             IndexPout = c_gsf-1;
00876             GetPout = true;
00877           }
00878           //const math::XYZTLorentzVector GsfMoment = itPfGsfPoint->momentum();
00879           ++c_gsf;
00880         }
00881       }
00882       math::XYZTLorentzVector pin = PfGsfPoint[0].momentum();      
00883       math::XYZTLorentzVector pout = PfGsfPoint[IndexPout].momentum();
00884       
00885       
00886     
00887       gsfEl = new reco::PFBlockElementGsfTrack(refgsf, pin, pout);
00888       
00889       bool valuegsf = true;
00890       // IMPORTANT SET T_FROM_GAMMACONV trackType() FOR CONVERSIONS
00891       gsfEl->setTrackType(reco::PFBlockElement::T_FROM_GAMMACONV, valuegsf);
00892 
00893       
00894 
00895       elements_.push_back( gsfEl);
00896       std::vector<reco::PFBrem> pfbrem = refgsf->PFRecBrem();
00897       
00898       for (unsigned i2=0;i2<pfbrem.size(); ++i2) {
00899         const double DP = pfbrem[i2].DeltaP();
00900         const double SigmaDP =  pfbrem[i2].SigmaDeltaP(); 
00901         const unsigned int TrajP = pfbrem[i2].indTrajPoint();
00902         if(TrajP == 99) continue;
00903 
00904         reco::PFBlockElement* bremEl;
00905         bremEl = new reco::PFBlockElementBrem(refgsf,DP,SigmaDP,TrajP);
00906         elements_.push_back(bremEl);
00907         
00908       }
00909     }
00910   }
00911 
00912   
00913   // -------------- ECAL clusters ---------------------
00914 
00915 
00916   if(ecalh.isValid() ) {
00917     //  pfcSCVec_.resize(ecalh->size(),-1);
00918  
00919     bNoSuperclus_  = (superClusters_.size() == 0);
00920     if (!bNoSuperclus_) pfcSCVec_.resize(ecalh->size(),-1);
00921 
00922 
00923     for(unsigned i=0;i<ecalh->size(); ++i)  {
00924 
00925       // this ecal cluster has been disabled
00926       if( !ecalMask.empty() &&
00927           !ecalMask[i] ) continue;
00928 
00929       reco::PFClusterRef ref( ecalh,i );
00930       reco::PFBlockElement* te
00931         = new reco::PFBlockElementCluster( ref,
00932                                            reco::PFBlockElement::ECAL);
00933       elements_.push_back( te );
00934 
00935       if (!bNoSuperclus_) {
00936 
00937         // Now mapping with Superclusters
00938         int scindex= ClusterClusterMapping::checkOverlap(*ref,superClusters_);
00939 
00940         if(scindex>=0)  {
00941           pfcSCVec_[ref.key()]=scindex;
00942           scpfcRefs_[scindex].push_back(ref);
00943         }
00944       }
00945 
00946     }
00947 
00948     bNoSuperclus_ = (scpfcRefs_.size() == 0);
00949 
00950   }
00951 
00952   // -------------- HCAL clusters ---------------------
00953 
00954   if(hcalh.isValid() ) {
00955     
00956     for(unsigned i=0;i<hcalh->size(); ++i)  {
00957       
00958       // this hcal cluster has been disabled
00959       if( !hcalMask.empty() &&
00960           !hcalMask[i] ) continue;
00961       
00962       reco::PFClusterRef ref( hcalh,i );
00963       reco::PFBlockElement* th
00964         = new reco::PFBlockElementCluster( ref,
00965                                            reco::PFBlockElement::HCAL );
00966       elements_.push_back( th );
00967     }
00968   }
00969 
00970 
00971   // -------------- HFEM clusters ---------------------
00972 
00973   if(hfemh.isValid() ) {
00974     
00975     for(unsigned i=0;i<hfemh->size(); ++i)  {
00976       
00977       // this hfem cluster has been disabled
00978       if( !hfemMask.empty() &&
00979           !hfemMask[i] ) continue;
00980       
00981       reco::PFClusterRef ref( hfemh,i );
00982       reco::PFBlockElement* th
00983         = new reco::PFBlockElementCluster( ref,
00984                                            reco::PFBlockElement::HFEM );
00985       elements_.push_back( th );
00986     }
00987   }
00988 
00989 
00990   // -------------- HFHAD clusters ---------------------
00991 
00992   if(hfhadh.isValid() ) {
00993     
00994     for(unsigned i=0;i<hfhadh->size(); ++i)  {
00995       
00996       // this hfhad cluster has been disabled
00997       if( !hfhadMask.empty() &&
00998           !hfhadMask[i] ) continue;
00999       
01000       reco::PFClusterRef ref( hfhadh,i );
01001       reco::PFBlockElement* th
01002         = new reco::PFBlockElementCluster( ref,
01003                                            reco::PFBlockElement::HFHAD );
01004       elements_.push_back( th );
01005     }
01006   }
01007 
01008 
01009 
01010 
01011   // -------------- PS clusters ---------------------
01012 
01013   if(psh.isValid() ) {
01014     for(unsigned i=0;i<psh->size(); ++i)  {
01015 
01016       // this ps cluster has been disabled
01017       if( !psMask.empty() &&
01018           !psMask[i] ) continue;
01019       reco::PFBlockElement::Type type = reco::PFBlockElement::NONE;
01020       reco::PFClusterRef ref( psh,i );
01021       // two types of elements:  PS1 (V) and PS2 (H) 
01022       // depending on layer:  PS1 or PS2
01023       switch(ref->layer()){
01024       case PFLayer::PS1:
01025         type = reco::PFBlockElement::PS1;
01026         break;
01027       case PFLayer::PS2:
01028         type = reco::PFBlockElement::PS2;
01029         break;
01030       default:
01031         break;
01032       }
01033       reco::PFBlockElement* tp
01034         = new reco::PFBlockElementCluster( ref,
01035                                            type );
01036       elements_.push_back( tp );
01037     }
01038   }
01039 
01040 
01041   // -------------- Loop over block elements ---------------------
01042 
01043   // Here we provide to all KDTree linkers the collections to link.
01044   // Glowinski & Gouzevitch
01045   
01046   for (std::list< reco::PFBlockElement* >::iterator it = elements_.begin();
01047        it != elements_.end(); ++it) {
01048     switch ((*it)->type()){
01049         
01050     case reco::PFBlockElement::TRACK:
01051       if (useKDTreeTrackEcalLinker_) {
01052         if ( (*it)->trackRefPF()->extrapolatedPoint( reco::PFTrajectoryPoint::ECALShowerMax ).isValid() )
01053           TELinker_.insertTargetElt(*it);
01054         if ( (*it)->trackRefPF()->extrapolatedPoint( reco::PFTrajectoryPoint::HCALEntrance ).isValid() )
01055           THLinker_.insertTargetElt(*it);
01056       }
01057       
01058       break;
01059 
01060     case reco::PFBlockElement::PS1:
01061     case reco::PFBlockElement::PS2:
01062       if (useKDTreeTrackEcalLinker_)
01063         PSELinker_.insertTargetElt(*it);
01064       break;
01065 
01066     case reco::PFBlockElement::HCAL:
01067       if (useKDTreeTrackEcalLinker_)
01068         THLinker_.insertFieldClusterElt(*it);
01069       break;
01070         
01071     case reco::PFBlockElement::ECAL:
01072       if (useKDTreeTrackEcalLinker_) {
01073         TELinker_.insertFieldClusterElt(*it);
01074         PSELinker_.insertFieldClusterElt(*it);
01075       }
01076       break;
01077 
01078     default:
01079       break;
01080     }
01081   }
01082 }
01083 
01084 
01085 template< template<typename> class T>
01086   void PFBlockAlgo::fillFromPhoton(const T<reco::PhotonCollection> & egh, unsigned isc, reco::PFBlockElementSuperCluster * pfbe) {
01087   reco::PhotonRef photonRef(egh,isc);
01088     pfbe->setTrackIso(photonRef->trkSumPtHollowConeDR04());
01089     pfbe->setEcalIso(photonRef->ecalRecHitSumEtConeDR04());
01090     pfbe->setHcalIso(photonRef->hcalTowerSumEtConeDR04());
01091     pfbe->setHoE(photonRef->hadronicOverEm());
01092     pfbe->setPhotonRef(photonRef);
01093     pfbe->setFromPhoton(true);
01094   }
01095 
01096 #endif
01097 
01098
/data/refman/pasoursint/CMSSW_4_4_5_patch3/src/RecoParticleFlow/PFProducer/interface/PFBlockAlgo.h
