CMS 3D CMS Logo

PFClusterAlgo.cc

Go to the documentation of this file.
00001 #include "RecoParticleFlow/PFClusterAlgo/interface/PFClusterAlgo.h"
00002 #include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
00003 #include "Math/GenVector/VectorUtil.h"
00004 
00005 
00006 #include <stdexcept>
00007 #include <string>
00008 
00009 using namespace std;
00010 
00011 unsigned PFClusterAlgo::prodNum_ = 1;
00012 
00013 //for debug only 
00014 //#define PFLOW_DEBUG
00015 
00016 PFClusterAlgo::PFClusterAlgo() :
00017   pfClusters_( new vector<reco::PFCluster> ),
00018   threshBarrel_(0.),
00019   threshSeedBarrel_(0.2),
00020   threshEndcap_(0.),
00021   threshSeedEndcap_(0.6),
00022   nNeighbours_(4),
00023   posCalcNCrystal_(-1),
00024   posCalcP1_(-1),
00025   showerSigma_(5),
00026   debug_(false) {}
00027 
00028 
00029 
00030 void PFClusterAlgo::doClustering( const PFRecHitHandle& rechitsHandle ) {
00031   rechitsHandle_ = rechitsHandle;
00032   doClustering( *rechitsHandle );
00033 }
00034 
00035 void PFClusterAlgo::doClustering( const reco::PFRecHitCollection& rechits ) {
00036 
00037 
00038   if(pfClusters_.get() ) pfClusters_->clear();
00039   else 
00040     pfClusters_.reset( new std::vector<reco::PFCluster> );
00041 
00042 
00043   eRecHits_.clear();
00044 
00045   bool initMask = false;
00046   if( mask_.size() != rechits.size() ) {
00047     initMask = true;
00048     mask_.clear();
00049     mask_.reserve( rechits.size() );
00050 
00051     if( ! mask_.empty() ) 
00052       cerr<<"PClusterAlgo::doClustering: map size should be "<<mask_.size()
00053           <<". Will be reinitialized."<<endl;    
00054   }
00055   
00056   color_.clear(); 
00057   color_.reserve( rechits.size() );
00058   seedStates_.clear();
00059   seedStates_.reserve( rechits.size() );
00060   usedInTopo_.clear();
00061   usedInTopo_.reserve( rechits.size() );
00062   
00063   for ( unsigned i = 0; i < rechits.size(); i++ ) {
00064     eRecHits_.insert( make_pair( rechit(i, rechits).energy(), i) );
00065     if(initMask) mask_.push_back( true );
00066     color_.push_back( 0 );     
00067     seedStates_.push_back( UNKNOWN ); 
00068     usedInTopo_.push_back( false ); 
00069   }  
00070 
00071   // look for seeds.
00072   findSeeds( rechits );
00073 
00074   // build topological clusters around seeds
00075   buildTopoClusters( rechits );
00076   
00077   // look for PFClusters inside each topological cluster (one per seed)
00078   for(unsigned i=0; i<topoClusters_.size(); i++) {
00079 
00080     const std::vector< unsigned >& topocluster = topoClusters_[i];
00081     buildPFClusters( topocluster, rechits ); 
00082   }
00083 }
00084 
00085 
00086 void PFClusterAlgo::setMask( const std::vector<bool>& mask ) {
00087   mask_ = mask;
00088 }
00089 
00090 
00091 
00092 
00093 double PFClusterAlgo::parameter( Parameter paramtype, 
00094                                  PFLayer::Layer layer)  const {
00095   
00096 
00097   double value = 0;
00098 
00099   switch( layer ) {
00100   case PFLayer::ECAL_BARREL:
00101   case PFLayer::HCAL_BARREL1:
00102   case PFLayer::HCAL_BARREL2: // I think this is HO. 
00103                               // should not do anything for HO !
00104     switch(paramtype) {
00105     case THRESH:
00106       value = threshBarrel_;
00107       break;
00108     case SEED_THRESH:
00109       value = threshSeedBarrel_;
00110       break;
00111     default:
00112       cerr<<"PFClusterAlgo::parameter : unknown parameter type "
00113           <<paramtype<<endl;
00114       assert(0);
00115     }
00116     break;
00117   case PFLayer::ECAL_ENDCAP:
00118   case PFLayer::HCAL_ENDCAP:
00119   case PFLayer::PS1:
00120   case PFLayer::PS2:
00121   case PFLayer::HCAL_HF: // should also be removed. probably no clusters 
00122     // and no particle flow in VFCAL
00123     switch(paramtype) {
00124     case THRESH:
00125       value = threshEndcap_;
00126       break;
00127     case SEED_THRESH:
00128       value = threshSeedEndcap_;
00129       break;
00130     default:
00131       cerr<<"PFClusterAlgo::parameter : unknown parameter type "
00132           <<paramtype<<endl;
00133       assert(0);
00134     }
00135     break;
00136   default:
00137     cerr<<"PFClusterAlgo::parameter : unknown layer "<<layer<<endl;
00138     assert(0);
00139     break;
00140   }
00141 
00142   return value;
00143 }
00144 
00145 
00146 
00147 void PFClusterAlgo::findSeeds( const reco::PFRecHitCollection& rechits ) {
00148 
00149   seeds_.clear();
00150 
00151   // should replace this by the message logger.
00152 #ifdef PFLOW_DEBUG
00153   if(debug_) 
00154     cout<<"PFClusterAlgo::findSeeds : start"<<endl;
00155 #endif
00156 
00157 
00158   // loop on rechits (sorted by decreasing energy - not E_T)
00159   for(EH ih = eRecHits_.begin(); ih != eRecHits_.end(); ih++ ) {
00160 
00161     unsigned  rhi      = ih->second; 
00162 
00163     if(! masked(rhi) ) continue;
00164     // rechit was asked to be processed
00165 
00166     double    rhenergy = ih->first;   
00167     const reco::PFRecHit& wannaBeSeed = rechit(rhi, rechits);
00168      
00169     if( seedStates_[rhi] == NO ) continue;
00170     // this hit was already tested, and is not a seed
00171  
00172     // determine seed energy threshold depending on the detector
00173     int layer = wannaBeSeed.layer();
00174     double seedThresh = parameter( SEED_THRESH, 
00175                                    static_cast<PFLayer::Layer>(layer) );
00176 
00177     
00178 #ifdef PFLOW_DEBUG
00179     if(debug_) 
00180       cout<<"layer:"<<layer<<" seedThresh:"<<seedThresh<<endl;
00181 #endif
00182 
00183 
00184     if( rhenergy < seedThresh ) {
00185       seedStates_[rhi] = NO; 
00186       continue;
00187     } 
00188 
00189       
00190     // Find the cell unused neighbours
00191     const vector<unsigned>* nbp;
00192 
00193     switch ( layer ) { 
00194     case PFLayer::ECAL_BARREL:         
00195     case PFLayer::ECAL_ENDCAP:       
00196     case PFLayer::HCAL_BARREL1:
00197     case PFLayer::HCAL_BARREL2:
00198     case PFLayer::HCAL_ENDCAP:
00199     case PFLayer::HCAL_HF:
00200       if( nNeighbours_ == 4 ) {
00201         nbp = & wannaBeSeed.neighbours4();
00202       }
00203       else if( nNeighbours_ == 8 )
00204         nbp = & wannaBeSeed.neighbours8();
00205       else {
00206         cerr<<"you're not allowed to set n neighbours to "
00207             <<nNeighbours_<<endl;
00208         assert(0);
00209       }
00210       break;
00211     case PFLayer::PS1:       
00212     case PFLayer::PS2:     
00213       nbp = & wannaBeSeed.neighbours4();
00214       break;
00215 
00216     default:
00217       cerr<<"CellsEF::PhotonSeeds : unknown layer "<<layer<<endl;
00218       assert(0);
00219     }
00220 
00221     const vector<unsigned>& neighbours = *nbp;
00222 
00223       
00224     // Select as a seed if all neighbours have a smaller energy
00225 
00226     seedStates_[rhi] = YES;
00227     for(unsigned in=0; in<neighbours.size(); in++) {
00228         
00229       const reco::PFRecHit& neighbour = rechit( neighbours[in], 
00230                                                 rechits ); 
00231         
00232       // one neighbour has a higher energy -> the tested rechit is not a seed
00233       if( neighbour.energy() > wannaBeSeed.energy() ) {
00234         seedStates_[rhi] = NO;
00235         break;
00236       }
00237     }
00238       
00239     if ( seedStates_[rhi] == YES ) {
00240 
00241       // seeds_ contains the indices of all seeds. 
00242       seeds_.push_back( rhi );
00243       
00244       // marking the rechit
00245       paint(rhi, SEED);
00246         
00247       // then all neighbours cannot be seeds and are flagged as such
00248       for(unsigned in=0; in<neighbours.size(); in++) {
00249         seedStates_[ neighbours[in] ] = NO;
00250       }
00251     }
00252   }  
00253 
00254 #ifdef PFLOW_DEBUG
00255   if(debug_) 
00256     cout<<"PFClusterAlgo::findSeeds : done"<<endl;
00257 #endif
00258 }
00259 
00260 
00261 
00262   
00263 void PFClusterAlgo::buildTopoClusters( const reco::PFRecHitCollection& rechits ){
00264 
00265   topoClusters_.clear(); 
00266   
00267 #ifdef PFLOW_DEBUG
00268   if(debug_) 
00269     cout<<"PFClusterAlgo::buildTopoClusters start"<<endl;
00270 #endif
00271   
00272   for(unsigned is = 0; is<seeds_.size(); is++) {
00273     
00274     unsigned rhi = seeds_[is];
00275 
00276     if( !masked(rhi) ) continue;
00277     // rechit was masked to be processed
00278 
00279     // already used in a topological cluster
00280     if( usedInTopo_[rhi] ) {
00281 #ifdef PFLOW_DEBUG
00282       if(debug_) 
00283         cout<<rhi<<" used"<<endl; 
00284 #endif
00285       continue;
00286     }
00287     
00288     vector< unsigned > topocluster;
00289     buildTopoCluster( topocluster, rhi, rechits );
00290    
00291     if(topocluster.empty() ) continue;
00292     
00293     topoClusters_.push_back( topocluster );
00294   } 
00295 
00296 #ifdef PFLOW_DEBUG
00297   if(debug_) 
00298     cout<<"PFClusterAlgo::buildTopoClusters done"<<endl;
00299 #endif
00300   
00301   return;
00302 }
00303 
00304 
00305 void PFClusterAlgo::buildTopoCluster( vector< unsigned >& cluster, 
00306                                       unsigned rhi, 
00307                                       const reco::PFRecHitCollection& rechits ){
00308 
00309 
00310 #ifdef PFLOW_DEBUG
00311   if(debug_)
00312     cout<<"PFClusterAlgo::buildTopoCluster in"<<endl;
00313 #endif
00314 
00315   const reco::PFRecHit& rh = rechit( rhi, rechits); 
00316 
00317   double e = rh.energy();
00318   int layer = rh.layer();
00319   
00320   double thresh = parameter( THRESH, 
00321                              static_cast<PFLayer::Layer>(layer) );
00322 
00323   if( e < thresh ) {
00324 #ifdef PFLOW_DEBUG
00325     if(debug_)
00326       cout<<"return : "<<e<<"<"<<thresh<<endl; 
00327 #endif
00328     return;
00329   }
00330 
00331   // add hit to cluster
00332 
00333   cluster.push_back( rhi );
00334   // idUsedRecHits_.insert( rh.detId() );
00335 
00336   usedInTopo_[ rhi ] = true;
00337 
00338   //   cout<<" hit ptr "<<hit<<endl;
00339 
00340   // get neighbours (1 common side)
00341   const std::vector< unsigned >& nbs = rh.neighbours4();
00342   
00343   for(unsigned i=0; i<nbs.size(); i++) {
00344 
00345 //     const reco::PFRecHit& neighbour = rechit( nbs[i], rechits );
00346 
00347 //     set<unsigned>::iterator used 
00348 //       = idUsedRecHits_.find( neighbour.detId() );
00349 //     if(used != idUsedRecHits_.end() ) continue;
00350     
00351     // already used
00352     if( usedInTopo_[ nbs[i] ] ) {
00353 #ifdef PFLOW_DEBUG
00354       if(debug_) 
00355         cout<<rhi<<" used"<<endl; 
00356 #endif
00357       continue;
00358     }
00359                              
00360     buildTopoCluster( cluster, nbs[i], rechits );
00361   }
00362 #ifdef PFLOW_DEBUG
00363   if(debug_)
00364     cout<<"PFClusterAlgo::buildTopoCluster out"<<endl;
00365 #endif
00366 
00367 }
00368 
00369 
00370 void 
00371 PFClusterAlgo::buildPFClusters( const std::vector< unsigned >& topocluster,
00372                                 const reco::PFRecHitCollection& rechits ) 
00373 {
00374 
00375 
00376   //  bool debug = false;
00377 
00378 
00379   // several rechits may be seeds. initialize PFClusters on these seeds. 
00380   
00381   vector<reco::PFCluster> curpfclusters;
00382   vector< unsigned > seedsintopocluster;
00383 
00384 
00385   for(unsigned i=0; i<topocluster.size(); i++ ) {
00386 
00387     unsigned rhi = topocluster[i];
00388 
00389     if( seedStates_[rhi] == YES ) {
00390 
00391       reco::PFCluster cluster;
00392 
00393       double fraction = 1.0; 
00394       
00395       reco::PFRecHitRef  recHitRef = createRecHitRef( rechits, rhi ); 
00396         
00397       cluster.addRecHitFraction( reco::PFRecHitFraction( recHitRef, 
00398                                                          fraction ) );
00399 
00400     // cluster.addRecHit( rhi, fraction );
00401       
00402       calculateClusterPosition( cluster, 
00403                                 true );    
00404 
00405 
00406 //       cout<<"PFClusterAlgo: 2"<<endl;
00407       curpfclusters.push_back( cluster );
00408 #ifdef PFLOW_DEBUG
00409       if(debug_) {
00410         cout << "PFClusterAlgo::buildPFClusters: seed "
00411              << rechit( rhi, rechits) <<endl;
00412         cout << "PFClusterAlgo::buildPFClusters: pfcluster initialized : "
00413              << cluster <<endl;
00414       }
00415 #endif
00416 
00417       // keep track of the seed of each topocluster
00418       seedsintopocluster.push_back( rhi );
00419       
00420     }
00421   }
00422 
00423   // if only one seed in the topocluster, use all crystals
00424   // in the position calculation (posCalcNCrystal = -1)
00425   // otherwise, use the user specified value
00426   int posCalcNCrystal = seedsintopocluster.size()>1 ? posCalcNCrystal_:-1;
00427     
00428   // Find iteratively the energy and position
00429   // of each pfcluster in the topological cluster
00430   unsigned iter = 0;
00431   unsigned niter = 50;
00432   double diff = 1.;
00433 
00434   // if(debug_) niter=2;
00435   while ( iter++ < niter && diff > 1E-8 ) {
00436 
00437     // Store previous iteration's result and reset pfclusters     
00438     vector<double> ener;
00439     vector<math::XYZVector> tmp;
00440 
00441     for ( unsigned ic=0; ic<curpfclusters.size(); ic++ ) {
00442       ener.push_back( curpfclusters[ic].energy() );
00443       
00444       math::XYZVector v;
00445       v = curpfclusters[ic].position();
00446 
00447       tmp.push_back( v );
00448 
00449 #ifdef PFLOW_DEBUG
00450       if(debug_)  {
00451         cout<<"saving photon pos "<<ic<<" "<<curpfclusters[ic]<<endl;
00452         cout<<tmp[ic].X()<<" "<<tmp[ic].Y()<<" "<<tmp[ic].Z()<<endl;
00453       }
00454 #endif
00455 
00456       curpfclusters[ic].reset();
00457     }
00458 
00459 
00460     // Loop over topocluster cells
00461     for( unsigned irh=0; irh<topocluster.size(); irh++ ) {
00462       
00463       unsigned rhindex = topocluster[irh];
00464       
00465       const reco::PFRecHit& rh = rechit( rhindex, rechits);
00466       
00467       // int layer = rh.layer();
00468              
00469       vector<double> dist;
00470       vector<double> frac;
00471       double fractot = 0.;
00472 
00473       bool isaseed = isSeed(rhindex);
00474 
00475       math::XYZVector cposxyzcell;
00476       cposxyzcell = rh.position();
00477 
00478 #ifdef PFLOW_DEBUG
00479       if(debug_) { 
00480         cout<<rh<<endl;
00481         cout<<"start loop on curpfclusters"<<endl;
00482       }
00483 #endif
00484 
00485       // Loop over pfclusters
00486       for ( unsigned ic=0; ic<tmp.size(); ic++) {
00487         
00488 #ifdef PFLOW_DEBUG
00489         if(debug_) cout<<"pfcluster "<<ic<<endl;
00490 #endif
00491         
00492         double frc=0.;
00493         bool seedexclusion=true;
00494 
00495         // convert cluster coordinates to xyz
00496         math::XYZVector cposxyzclust( tmp[ic].X(), tmp[ic].Y(), tmp[ic].Z() );
00497         
00498 #ifdef PFLOW_DEBUG
00499         if(debug_) {
00500           
00501           cout<<"CLUSTER "<<cposxyzclust.X()<<","
00502               <<cposxyzclust.Y()<<","
00503               <<cposxyzclust.Z()<<"\t\t"
00504               <<"CELL "<<cposxyzcell.X()<<","
00505               <<cposxyzcell.Y()<<","
00506               <<cposxyzcell.Z()<<endl;
00507         }  
00508 #endif
00509         
00510         // Compute the distance between the current cell 
00511         // and the current PF cluster, normalized to a 
00512         // number of "sigma"
00513         math::XYZVector deltav = cposxyzclust;
00514         deltav -= cposxyzcell;
00515         double d = deltav.R() / showerSigma_;
00516         
00517         // if distance cell-cluster is too large, it means that 
00518         // we're right on the junction between 2 subdetectors (HCAL/VFCAL...)
00519         // in this case, distance is calculated in the xy plane
00520         // could also be a large supercluster... 
00521 #ifdef PFLOW_DEBUG
00522         if( d > 6. && debug_ ) { 
00523           for( unsigned jrh=0; jrh<topocluster.size(); jrh++ ) {
00524             paint(jrh, SPECIAL);
00525           }
00526           cout<<"PFClusterAlgo Warning: distance too large"<<d<<endl;
00527         }
00528 #endif
00529         dist.push_back( d );
00530 
00531         // the current cell is the seed from the current photon.
00532         if( rhindex == seedsintopocluster[ic] && seedexclusion ) {
00533           frc = 1.;
00534 #ifdef PFLOW_DEBUG
00535           if(debug_) cout<<"this cell is a seed for the current photon"<<endl;
00536 #endif
00537         }
00538         else if( isaseed && seedexclusion ) {
00539           frc = 0.;
00540 #ifdef PFLOW_DEBUG
00541           if(debug_) cout<<"this cell is a seed for another photon"<<endl;
00542 #endif
00543         }
00544         else {
00545           // Compute the fractions of the cell energy to be assigned to 
00546           // each curpfclusters in the cluster.
00547           frc = ener[ic] * exp ( - dist[ic]*dist[ic] / 2. );
00548 
00549 #ifdef PFLOW_DEBUG
00550           if(debug_) {
00551             cout<<"dist["<<ic<<"] "<<dist[ic]
00552                 <<", sigma="<<sigma
00553                 <<", frc="<<frc<<endl;
00554           }  
00555 #endif
00556         
00557         }
00558         fractot += frc;
00559         frac.push_back(frc);
00560       }      
00561 
00562       // Add the relevant fraction of the cell to the curpfclusters
00563 #ifdef PFLOW_DEBUG
00564       if(debug_) cout<<"start add cell"<<endl;
00565 #endif
00566       for ( unsigned ic=0; ic<tmp.size(); ++ic ) {
00567 #ifdef PFLOW_DEBUG
00568         if(debug_) 
00569           cout<<" frac["<<ic<<"] "<<frac[ic]<<" "<<fractot<<" "<<rh<<endl;
00570 #endif
00571 
00572         if( fractot ) 
00573           frac[ic] /= fractot;
00574         else { 
00575 #ifdef PFLOW_DEBUG
00576           if( debug_ ) {
00577             int layer = rh.layer();
00578             cerr<<"fractot = 0 ! "<<layer<<endl;
00579             
00580             for( unsigned trh=0; trh<topocluster.size(); trh++ ) {
00581               unsigned tindex = topocluster[trh];
00582               const reco::PFRecHit& rh = rechit( tindex, rechits);
00583               cout<<rh<<endl;
00584             }
00585 
00586             // assert(0)
00587           }
00588 #endif
00589 
00590           continue;
00591         }
00592 
00593         // if the fraction has been set to 0, the cell 
00594         // is now added to the cluster - careful ! (PJ, 19/07/08)
00595         // BUT KEEP ONLY CLOSE CELLS OTHERWISE MEMORY JUST EXPLOSES
00596         // (PJ, 15/09/08 <- similar to what existed before the 
00597         // previous bug fix, but keeps the close seeds inside, 
00598         // even if their fraction was set to zero.)
00599         // Also add a protection to keep the seed in the cluster 
00600         // when the latter gets far from the former. These cases
00601         // (about 1% of the clusters) need to be studied, as 
00602         // they create fake photons, in general.
00603         // (PJ, 16/09/08) 
00604         if ( dist[ic] < 6. || frac[ic] > 0.99999 ) { 
00605           // if ( dist[ic] > 6. ) cout << "Warning : PCluster is getting very far from its seeding cell" << endl;
00606           reco::PFRecHitRef  recHitRef = createRecHitRef( rechits, rhindex ); 
00607           reco::PFRecHitFraction rhf( recHitRef,frac[ic] );
00608           curpfclusters[ic].addRecHitFraction( rhf );
00609         }
00610       }
00611       // if(debug_) cout<<" end add cell"<<endl;
00612       
00613       dist.clear();
00614       frac.clear();
00615     }
00616     
00617     // Determine the new cluster position and check 
00618     // the distance with the previous iteration
00619     diff = 0.;
00620     for (  unsigned ic=0; ic<tmp.size(); ++ic ) {
00621       calculateClusterPosition(curpfclusters[ic], true, posCalcNCrystal);
00622 #ifdef PFLOW_DEBUG
00623       if(debug_) cout<<"new iter "<<ic<<endl;
00624       if(debug_) cout<<curpfclusters[ic]<<endl;
00625 #endif
00626 
00627       double delta = ROOT::Math::VectorUtil::DeltaR(curpfclusters[ic].position(),tmp[ic]);
00628       if ( delta > diff ) diff = delta;
00629     }
00630     ener.clear();
00631     tmp.clear();
00632   }
00633   
00634   // Issue a warning message if the number of iterations 
00635   // exceeds 50
00636 #ifdef PFLOW_DEBUG
00637   if ( iter >= 50 && debug_ ) 
00638     cout << "PFClusterAlgo Warning: "
00639          << "more than "<<niter<<" iterations in pfcluster finding: " 
00640          <<  setprecision(10) << diff << endl;
00641 #endif
00642   
00643   // There we go
00644   // add all clusters to the list of pfClusters.
00645   for(unsigned ic=0; ic<curpfclusters.size(); ic++) {
00646     calculateClusterPosition(curpfclusters[ic], true, posCalcNCrystal);
00647     pfClusters_->push_back(curpfclusters[ic]); 
00648   }
00649 }
00650 
00651 
00652 
00653 void 
00654 PFClusterAlgo::calculateClusterPosition(reco::PFCluster& cluster,
00655                                         bool depcor, 
00656                                         int posCalcNCrystal) {
00657 
00658   if( posCalcNCrystal_ != -1 && 
00659       posCalcNCrystal_ != 5 && 
00660       posCalcNCrystal_ != 9 ) {
00661     throw "PFCluster::calculatePosition : posCalcNCrystal_ must be -1, 5, or 9.";
00662   }  
00663 
00664   if(!posCalcNCrystal) posCalcNCrystal = posCalcNCrystal_; 
00665 
00666   cluster.position_.SetXYZ(0,0,0);
00667 
00668   cluster.energy_ = 0;
00669   
00670   double normalize = 0;
00671 
00672   // calculate total energy, average layer, and look for seed  ---------- //
00673 
00674   // double layer = 0;
00675   map <PFLayer::Layer, double> layers; 
00676   unsigned seedIndex = 0;
00677   bool     seedIndexFound = false;
00678 
00679   //Colin: the following can be simplified!
00680 
00681   // loop on rechit fractions
00682   for (unsigned ic=0; ic<cluster.rechits_.size(); ic++ ) {
00683 
00684     unsigned rhi = cluster.rechits_[ic].recHitRef().index();
00685     // const reco::PFRecHit& rh = rechit( rhi, rechits );
00686 
00687     const reco::PFRecHit& rh = *(cluster.rechits_[ic].recHitRef());
00688     double fraction =  cluster.rechits_[ic].fraction();
00689 
00690     // Find the seed of this sub-cluster (excluding other seeds found in the topological
00691     // cluster, the energy fraction of which were set to 0 fpr the position determination.
00692     if( isSeed(rhi) && fraction > 1e-9 ) {
00693       seedIndex = rhi;
00694       seedIndexFound = true;
00695     }
00696 
00697 
00698     double recHitEnergy = rh.energy() * fraction;
00699     cluster.energy_ += recHitEnergy;
00700 
00701     // sum energy in each layer
00702     PFLayer::Layer layer = rh.layer();                              
00703     map <PFLayer::Layer, double>:: iterator it = layers.find(layer);
00704     if (it != layers.end()) 
00705       it->second += recHitEnergy;
00706     else 
00707       layers.insert(make_pair(layer, recHitEnergy));
00708   }  
00709 
00710   assert(seedIndexFound);
00711 
00712   // loop over pairs to find layer with max energy          
00713   double Emax = 0.;
00714   PFLayer::Layer layer = PFLayer::NONE;
00715   for (map<PFLayer::Layer, double>::iterator it = layers.begin();
00716        it != layers.end(); ++it) {
00717     double e = it->second;
00718     if(e > Emax){ 
00719       Emax = e; 
00720       layer = it->first;
00721     }
00722   }
00723   
00724   //setlayer here
00725   cluster.setLayer( layer ); // take layer with max energy
00726 
00727   // layer /= cluster.energy_;
00728   // cluster.layer_ = lrintf(layer); // nearest integer
00729 
00730   double p1 =  posCalcP1_;
00731   if( p1 < 0 ) { 
00732     // automatic (and hopefully best !) determination of the parameter
00733     // for position determination.
00734     
00735     // Remove the ad-hoc determination of p1, and set it to the 
00736     // seed threshold.
00737     switch(cluster.layer() ) {
00738     case PFLayer::ECAL_BARREL:
00739     case PFLayer::HCAL_BARREL1:
00740     case PFLayer::HCAL_BARREL2:
00741       p1 = threshBarrel_;
00742       break;
00743     case PFLayer::ECAL_ENDCAP:
00744     case PFLayer::HCAL_ENDCAP:
00745     case PFLayer::HCAL_HF:
00746       p1 = threshEndcap_;
00747       break;
00748 
00749     /*
00750     switch(cluster.layer() ) {
00751     case PFLayer::ECAL_BARREL:
00752       p1 = 0.004 + 0.022*cluster.energy_; // 27 feb 2006 
00753       break;
00754     case PFLayer::ECAL_ENDCAP:
00755       p1 = 0.014 + 0.009*cluster.energy_; // 27 feb 2006 
00756       break;
00757     case PFLayer::HCAL_BARREL1:
00758     case PFLayer::HCAL_BARREL2:
00759     case PFLayer::HCAL_ENDCAP:
00760     case PFLayer::HCAL_HF:
00761       p1 = 5.41215e-01 * log( cluster.energy_ / 1.29803e+01 );
00762       if(p1<0.01) p1 = 0.01;
00763       break;
00764     */
00765 
00766     default:
00767       cerr<<"Clusters weight_p1 -1 not yet allowed for layer "<<layer
00768           <<". Chose a better value in the opt file"<<endl;
00769       assert(0);
00770       break;
00771     }
00772   } 
00773   else if( p1< 1e-9 ) { // will divide by p1 later on
00774     p1 = 1e-9;
00775   }
00776 
00777   // calculate uncorrected cluster position --------------------------------
00778 
00779   reco::PFCluster::REPPoint clusterpos;   // uncorrected cluster pos 
00780   math::XYZPoint clusterposxyz;           // idem, xyz coord 
00781   math::XYZPoint firstrechitposxyz;       // pos of the rechit with highest E
00782 
00783   double maxe = -9999;
00784   double x = 0;
00785   double y = 0;
00786   double z = 0;
00787   
00788   for (unsigned ic=0; ic<cluster.rechits_.size(); ic++ ) {
00789     
00790     unsigned rhi = cluster.rechits_[ic].recHitRef().index();
00791 //     const reco::PFRecHit& rh = rechit( rhi, rechits );
00792 
00793     const reco::PFRecHit& rh = *(cluster.rechits_[ic].recHitRef());
00794 
00795     if(rhi != seedIndex) { // not the seed
00796       if( posCalcNCrystal == 5 ) { // pos calculated from the 5 neighbours only
00797         if(!rh.isNeighbour4(seedIndex) ) {
00798           continue;
00799         }
00800       }
00801       if( posCalcNCrystal == 9 ) { // pos calculated from the 9 neighbours only
00802         if(!rh.isNeighbour8(seedIndex) ) {
00803           continue;
00804         }
00805       }
00806     }
00807     
00808     double fraction =  cluster.rechits_[ic].fraction();
00809     double recHitEnergy = rh.energy() * fraction;
00810 
00811     double norm = max(0., log(recHitEnergy/p1 ));
00812     
00813     const math::XYZPoint& rechitposxyz = rh.position();
00814     
00815     if( recHitEnergy > maxe ) {
00816       firstrechitposxyz = rechitposxyz;
00817       maxe = recHitEnergy;
00818     }
00819 
00820     x += rechitposxyz.X() * norm;
00821     y += rechitposxyz.Y() * norm;
00822     z += rechitposxyz.Z() * norm;
00823     
00824     // clusterposxyz += rechitposxyz * norm;
00825     normalize += norm;
00826   }
00827   
00828   // normalize uncorrected position
00829   // assert(normalize);
00830   if( normalize < 1e-9 ) {
00831     //    cerr<<"--------------------"<<endl;
00832     //    cerr<<(*this)<<endl;
00833     cout << "Watch out : cluster too far from its seeding cell, set to 0,0,0" << endl;
00834     clusterposxyz.SetXYZ(0,0,0);
00835     clusterpos.SetCoordinates(0,0,0);
00836     return;
00837   }
00838   else {
00839     x /= normalize;
00840     y /= normalize; 
00841     z /= normalize; 
00842     clusterposxyz.SetCoordinates( x, y, z);
00843     clusterpos.SetCoordinates( clusterposxyz.Rho(), clusterposxyz.Eta(), clusterposxyz.Phi() );
00844   }  
00845 
00846   cluster.posrep_ = clusterpos;
00847   cluster.position_ = clusterposxyz;
00848 
00849 
00850   // correction of the rechit position, 
00851   // according to the depth, only for ECAL 
00852 
00853 
00854   if( depcor &&   // correction requested and ECAL
00855       ( cluster.layer() == PFLayer::ECAL_BARREL ||       
00856         cluster.layer() == PFLayer::ECAL_ENDCAP ) ) {
00857 
00858     
00859     double corra = reco::PFCluster::depthCorA_;
00860     double corrb = reco::PFCluster::depthCorB_;
00861     if( abs(clusterpos.Eta() ) < 2.6 && 
00862         abs(clusterpos.Eta() ) > 1.65   ) { 
00863       // if crystals under preshower, correction is not the same  
00864       // (shower depth smaller)
00865       corra = reco::PFCluster::depthCorAp_;
00866       corrb = reco::PFCluster::depthCorBp_;
00867     }
00868 
00869     double depth = 0;
00870 
00871     switch( reco::PFCluster::depthCorMode_ ) {
00872     case 1: // for e/gamma 
00873       depth = corra * ( corrb + log(cluster.energy_) ); 
00874       break;
00875     case 2: // for hadrons
00876       depth = corra;
00877       break;
00878     default:
00879       cerr<<"PFClusterAlgo::calculateClusterPosition : unknown function for depth correction! "<<endl;
00880       assert(0);
00881     }
00882 
00883     // calculate depth vector:
00884     // its mag is depth
00885     // its direction is the cluster direction (uncorrected)
00886 
00887 //     double xdepthv = clusterposxyz.X();
00888 //     double ydepthv = clusterposxyz.Y();
00889 //     double zdepthv = clusterposxyz.Z();
00890 //     double mag = sqrt( xdepthv*xdepthv + 
00891 //                     ydepthv*ydepthv + 
00892 //                     zdepthv*zdepthv );
00893     
00894 
00895 //     math::XYZPoint depthv(clusterposxyz); 
00896 //     depthv.SetMag(depth);
00897     
00898     
00899     math::XYZVector depthv( clusterposxyz.X(), 
00900                             clusterposxyz.Y(),
00901                             clusterposxyz.Z() );
00902     depthv /= sqrt(depthv.Mag2() );
00903     depthv *= depth;
00904 
00905 
00906     // now calculate corrected cluster position:    
00907     math::XYZPoint clusterposxyzcor;
00908 
00909     maxe = -9999;
00910     x = 0;
00911     y = 0;
00912     z = 0;
00913     cluster.posrep_.SetXYZ(0,0,0);
00914     normalize = 0;
00915     for (unsigned ic=0; ic<cluster.rechits_.size(); ic++ ) {
00916 
00917       unsigned rhi = cluster.rechits_[ic].recHitRef().index();
00918 //       const reco::PFRecHit& rh = rechit( rhi, rechits );
00919       
00920       const reco::PFRecHit& rh = *(cluster.rechits_[ic].recHitRef());
00921 
00922       if(rhi != seedIndex) {
00923         if( posCalcNCrystal == 5 ) {
00924           if(!rh.isNeighbour4(seedIndex) ) {
00925             continue;
00926           }
00927         }
00928         if( posCalcNCrystal == 9 ) {
00929           if(!rh.isNeighbour8(seedIndex) ) {
00930             continue;
00931           }
00932         }
00933       }
00934     
00935       double fraction =  cluster.rechits_[ic].fraction();
00936       double recHitEnergy = rh.energy() * fraction;
00937       
00938       const math::XYZPoint&  rechitposxyz = rh.position();
00939 
00940       // rechit axis not correct ! 
00941       math::XYZVector rechitaxis = rh.getAxisXYZ();
00942       // rechitaxis -= math::XYZVector( rechitposxyz.X(), rechitposxyz.Y(), rechitposxyz.Z() );
00943       
00944       math::XYZVector rechitaxisu( rechitaxis );
00945       rechitaxisu /= sqrt( rechitaxis.Mag2() );
00946 
00947       math::XYZVector displacement( rechitaxisu );
00948       // displacement /= sqrt( displacement.Mag2() );    
00949       displacement *= rechitaxisu.Dot( depthv );
00950       
00951       math::XYZPoint rechitposxyzcor( rechitposxyz );
00952       rechitposxyzcor += displacement;
00953 
00954       if( recHitEnergy > maxe ) {
00955         firstrechitposxyz = rechitposxyzcor;
00956         maxe = recHitEnergy;
00957       }
00958       
00959       double norm = max(0., log(recHitEnergy/p1 ));
00960       
00961       x += rechitposxyzcor.X() * norm;
00962       y += rechitposxyzcor.Y() * norm;
00963       z += rechitposxyzcor.Z() * norm;
00964       
00965       // clusterposxyzcor += rechitposxyzcor * norm;
00966       normalize += norm;
00967     }
00968     
00969     // normalize
00970     if(normalize < 1e-9) {
00971       cerr<<"--------------------"<<endl;
00972       cerr<< cluster <<endl;
00973       assert(0);
00974     }
00975     else {
00976       x /= normalize;
00977       y /= normalize;
00978       z /= normalize;
00979       
00980 
00981       clusterposxyzcor.SetCoordinates(x,y,z);
00982       cluster.posrep_.SetCoordinates( clusterposxyzcor.Rho(), 
00983                                       clusterposxyzcor.Eta(), 
00984                                       clusterposxyzcor.Phi() );
00985       cluster.position_  = clusterposxyzcor;
00986       clusterposxyz = clusterposxyzcor;
00987     }
00988   }
00989 }
00990 
00991 
00992 
00993 const reco::PFRecHit& 
00994 PFClusterAlgo::rechit(unsigned i, 
00995                       const reco::PFRecHitCollection& rechits ) {
00996 
00997   if( i < 0 || i >= rechits.size() ) {
00998     string err = "PFClusterAlgo::rechit : out of range";
00999     throw std::out_of_range(err);
01000   }
01001   
01002   return rechits[i];
01003 }
01004 
01005 
01006 
01007 bool PFClusterAlgo::masked(unsigned rhi) const {
01008 
01009   if(rhi<0 || rhi>=mask_.size() ) {
01010     string err = "PFClusterAlgo::masked : out of range";
01011     throw std::out_of_range(err);
01012   }
01013   
01014   return mask_[rhi];
01015 }
01016 
01017 
01018 unsigned PFClusterAlgo::color(unsigned rhi) const {
01019 
01020   if(rhi<0 || rhi>=color_.size() ) {
01021     string err = "PFClusterAlgo::color : out of range";
01022     throw std::out_of_range(err);
01023   }
01024   
01025   return color_[rhi];
01026 }
01027 
01028 
01029 
01030 bool PFClusterAlgo::isSeed(unsigned rhi) const {
01031 
01032   if(rhi<0 || rhi>=seedStates_.size() ) {
01033     string err = "PFClusterAlgo::isSeed : out of range";
01034     throw std::out_of_range(err);
01035   }
01036   
01037   return seedStates_[rhi]>0 ? true : false;
01038 }
01039 
01040 
01041 void PFClusterAlgo::paint(unsigned rhi, unsigned color ) {
01042 
01043   if(rhi<0 || rhi>=color_.size() ) {
01044     string err = "PFClusterAlgo::color : out of range";
01045     throw std::out_of_range(err);
01046   }
01047   
01048   color_[rhi] = color;
01049 }
01050 
01051 
01052 reco::PFRecHitRef 
01053 PFClusterAlgo::createRecHitRef( const reco::PFRecHitCollection& rechits, 
01054                                 unsigned rhi ) {
01055 
01056   if( rechitsHandle_.isValid() ) {
01057     return reco::PFRecHitRef(  rechitsHandle_, rhi );
01058   } 
01059   else {
01060     return reco::PFRecHitRef(  &rechits, rhi );
01061   }
01062 }
01063 
01064 
01065 ostream& operator<<(ostream& out,const PFClusterAlgo& algo) {
01066   if(!out) return out;
01067   out<<"PFClusterAlgo parameters : "<<endl;
01068   out<<"-----------------------------------------------------"<<endl;
01069   out<<"threshBarrel     : "<<algo.threshBarrel_     <<endl;
01070   out<<"threshSeedBarrel : "<<algo.threshSeedBarrel_ <<endl;
01071   out<<"threshEndcap     : "<<algo.threshEndcap_     <<endl;
01072   out<<"threshSeedEndcap : "<<algo.threshSeedEndcap_ <<endl;
01073   out<<"nNeighbours      : "<<algo.nNeighbours_      <<endl;
01074   out<<"posCalcNCrystal  : "<<algo.posCalcNCrystal_  <<endl;
01075   out<<"posCalcP1        : "<<algo.posCalcP1_        <<endl;
01076   out<<"showerSigma      : "<<algo.showerSigma_      <<endl;
01077 
01078   out<<endl;
01079   out<<algo.pfClusters_->size()<<" clusters:"<<endl;
01080 
01081   for(unsigned i=0; i<algo.pfClusters_->size(); i++) {
01082     out<<(*algo.pfClusters_)[i]<<endl;
01083     
01084     if(!out) return out;
01085   }
01086   
01087   return out;
01088 }

Generated on Tue Jun 9 17:44:40 2009 for CMSSW by  doxygen 1.5.4