CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_2_SLHC4_patch1/src/RecoParticleFlow/PFProducer/src/KDTreeLinkerPSEcal.cc

Go to the documentation of this file.
00001 #include "RecoParticleFlow/PFProducer/interface/KDTreeLinkerPSEcal.h"
00002 
00003 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
00004 #include "TMath.h"
00005 
00006 KDTreeLinkerPSEcal::KDTreeLinkerPSEcal()
00007   : KDTreeLinkerBase(),
00008     resPSpitch_ (0.19),
00009     resPSlength_ (6.1),
00010     ps1ToEcal_ (1.072),
00011     ps2ToEcal_ (1.057)
00012 {}
00013 
00014 KDTreeLinkerPSEcal::~KDTreeLinkerPSEcal()
00015 {
00016   clear();
00017 }
00018 
00019 
00020 void
00021 KDTreeLinkerPSEcal::insertTargetElt(reco::PFBlockElement        *psCluster)
00022 {
00023   reco::PFClusterRef clusterref = psCluster->clusterRef();
00024 
00025   // This test is more or less done in PFBlockAlgo.h. In others cases, it should be switch on.
00026   //   if (!((clusterref->layer() == PFLayer::PS1) || 
00027   //    (clusterref->layer() == PFLayer::PS2)))
00028   //     return;
00029 
00030   targetSet_.insert(psCluster);
00031 }
00032 
00033 
00034 void
00035 KDTreeLinkerPSEcal::insertFieldClusterElt(reco::PFBlockElement  *ecalCluster)
00036 {
00037   reco::PFClusterRef clusterref = ecalCluster->clusterRef();
00038 
00039   if (clusterref->layer() != PFLayer::ECAL_ENDCAP)
00040     return;
00041           
00042   const std::vector<reco::PFRecHitFraction> &fraction = clusterref->recHitFractions();
00043 
00044   // We create a list of cluster
00045   fieldClusterSet_.insert(ecalCluster);
00046 
00047   double clusterz = clusterref->position().Z();
00048   RecHitSet& rechitsSet = (clusterz < 0) ? rechitsNegSet_ : rechitsPosSet_;
00049 
00050   for(size_t rhit = 0; rhit < fraction.size(); ++rhit) {
00051     const reco::PFRecHitRef& rh = fraction[rhit].recHitRef();
00052     double fract = fraction[rhit].fraction();
00053 
00054     if ((rh.isNull()) || (fract < 1E-4))
00055       continue;
00056       
00057     const reco::PFRecHit& rechit = *rh;
00058       
00059     // We save the links rechit to Clusters
00060     rechit2ClusterLinks_[&rechit].insert(ecalCluster);
00061     
00062     // We create a liste of rechits
00063     rechitsSet.insert(&rechit);
00064   }
00065 }
00066 
00067 void 
00068 KDTreeLinkerPSEcal::buildTree()
00069 {
00070   buildTree(rechitsNegSet_, treeNeg_);
00071   buildTree(rechitsPosSet_, treePos_);
00072 }
00073 
00074 void 
00075 KDTreeLinkerPSEcal::buildTree(const RecHitSet   &rechitsSet,
00076                               KDTreeLinkerAlgo  &tree)
00077 {
00078   // List of pseudo-rechits that will be used to create the KDTree
00079   std::vector<KDTreeNodeInfo> eltList;
00080 
00081   // Filling of this eltList
00082   for(RecHitSet::const_iterator it = rechitsSet.begin(); 
00083       it != rechitsSet.end(); it++) {
00084 
00085     const reco::PFRecHit* rh = *it;
00086     const math::XYZPoint& posxyz = rh->position();
00087         
00088     KDTreeNodeInfo rhinfo (rh, posxyz.X(), posxyz.Y());
00089     eltList.push_back(rhinfo);
00090   }
00091 
00092   // xmin-xmax, ymain-ymax
00093   KDTreeBox region(-150., 150., -150., 150.);
00094 
00095   // We may now build the KDTree
00096   tree.build(eltList, region);
00097 }
00098 
00099 void
00100 KDTreeLinkerPSEcal::searchLinks()
00101 {
00102   // Must of the code has been taken from LinkByRecHit.cc
00103 
00104   // We iterate over the PS clusters.
00105   for(BlockEltSet::iterator it = targetSet_.begin(); 
00106       it != targetSet_.end(); it++) {
00107 
00108     (*it)->setIsValidMultilinks(true);
00109         
00110     reco::PFClusterRef clusterPSRef = (*it)->clusterRef();
00111     const reco::PFCluster& clusterPS = *clusterPSRef;
00112 
00113     // PS cluster position, extrapolated to ECAL
00114     double zPS = clusterPS.position().Z();
00115     double xPS = clusterPS.position().X();
00116     double yPS = clusterPS.position().Y();
00117 
00118     double etaPS = fabs(clusterPS.positionREP().eta());
00119     double deltaX = 0.;
00120     double deltaY = 0.;
00121     double xPSonEcal = xPS;
00122     double yPSonEcal = yPS;
00123 
00124     if (clusterPS.layer() == PFLayer::PS1) { // PS1
00125 
00126       // vertical strips, measure x with pitch precision
00127       deltaX = resPSpitch_;
00128       deltaY = resPSlength_;
00129       xPSonEcal *= ps1ToEcal_;
00130       yPSonEcal *= ps1ToEcal_;
00131 
00132     } else { // PS2
00133 
00134       // horizontal strips, measure y with pitch precision
00135       deltaY = resPSpitch_;
00136       deltaX = resPSlength_;
00137       xPSonEcal *= ps2ToEcal_;
00138       yPSonEcal *= ps2ToEcal_;
00139 
00140     }
00141  
00142     
00143     // Estimate the maximal envelope in phi/eta that will be used to find rechit candidates.
00144     // Same envelope for cap et barrel rechits.
00145     
00146     
00147     double maxEcalRadius = getCristalXYMaxSize() / 2.;
00148 
00149     // The inflation factor includes the approximate projection from Preshower to ECAL
00150     double inflation = 2.4 - (etaPS-1.6);
00151     double rangeX = maxEcalRadius * (1 + (0.05 + 1.0 / maxEcalRadius * deltaX / 2.)) * inflation; 
00152     double rangeY = maxEcalRadius * (1 + (0.05 + 1.0 / maxEcalRadius * deltaY / 2.)) * inflation; 
00153     
00154     // We search for all candidate recHits, ie all recHits contained in the maximal size envelope.
00155     std::vector<KDTreeNodeInfo> recHits;
00156     KDTreeBox trackBox(xPSonEcal - rangeX, xPSonEcal + rangeX, 
00157                   yPSonEcal - rangeY, yPSonEcal + rangeY);
00158 
00159     if (zPS < 0)
00160       treeNeg_.search(trackBox, recHits);
00161     else
00162       treePos_.search(trackBox, recHits);
00163 
00164 
00165     for(std::vector<KDTreeNodeInfo>::const_iterator rhit = recHits.begin(); 
00166         rhit != recHits.end(); ++rhit) {
00167            
00168       const std::vector< math::XYZPoint >& corners = rhit->ptr->getCornersXYZ();
00169       if(corners.size() != 4) continue;
00170 
00171       // Find all clusters associated to given rechit
00172       RecHit2BlockEltMap::iterator ret = rechit2ClusterLinks_.find(rhit->ptr);
00173       
00174       for(BlockEltSet::const_iterator clusterIt = ret->second.begin(); 
00175           clusterIt != ret->second.end(); clusterIt++) {
00176         
00177         reco::PFClusterRef clusterref = (*clusterIt)->clusterRef();
00178         double clusterz = clusterref->position().Z();
00179 
00180         const math::XYZPoint& posxyz = rhit->ptr->position() * zPS / clusterz;
00181 
00182         double x[5];
00183         double y[5];
00184         for ( unsigned jc=0; jc<4; ++jc ) {
00185           math::XYZPoint cornerpos = corners[jc] * zPS / clusterz;
00186           x[jc] = cornerpos.X() + (cornerpos.X()-posxyz.X()) * (0.05 +1.0/fabs((cornerpos.X()-posxyz.X()))*deltaX/2.);
00187           y[jc] = cornerpos.Y() + (cornerpos.Y()-posxyz.Y()) * (0.05 +1.0/fabs((cornerpos.Y()-posxyz.Y()))*deltaY/2.);
00188         }
00189 
00190         x[4] = x[0];
00191         y[4] = y[0];
00192         
00193         bool isinside = TMath::IsInside(xPS,
00194                                         yPS,
00195                                         5,x,y);
00196         
00197         // Check if the track and the cluster are linked
00198         if( isinside )
00199           target2ClusterLinks_[*it].insert(*clusterIt); 
00200       }
00201     }
00202     
00203   }
00204 }
00205 
00206 void
00207 KDTreeLinkerPSEcal::updatePFBlockEltWithLinks()
00208 {
00209   //TODO YG : Check if cluster positionREP() is valid ?
00210 
00211   // Here we save in each track the list of phi/eta values of linked clusters.
00212   for (BlockElt2BlockEltMap::iterator it = target2ClusterLinks_.begin();
00213        it != target2ClusterLinks_.end(); ++it) {
00214     reco::PFMultiLinksTC multitracks(true);
00215 
00216     for (BlockEltSet::iterator jt = it->second.begin();
00217          jt != it->second.end(); ++jt) {
00218 
00219       double clusterPhi = (*jt)->clusterRef()->positionREP().Phi();
00220       double clusterEta = (*jt)->clusterRef()->positionREP().Eta();
00221 
00222       multitracks.linkedClusters.push_back(std::make_pair(clusterPhi, clusterEta));
00223     }
00224 
00225     it->first->setMultilinks(multitracks);
00226   }
00227 }
00228 
00229 void
00230 KDTreeLinkerPSEcal::clear()
00231 {
00232   targetSet_.clear();
00233   fieldClusterSet_.clear();
00234 
00235   rechitsNegSet_.clear();
00236   rechitsPosSet_.clear();
00237 
00238   rechit2ClusterLinks_.clear();
00239   target2ClusterLinks_.clear();
00240 
00241   treeNeg_.clear();
00242   treePos_.clear();
00243 }