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