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
00014
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
00072 findSeeds( rechits );
00073
00074
00075 buildTopoClusters( rechits );
00076
00077
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:
00103
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:
00122
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
00152 #ifdef PFLOW_DEBUG
00153 if(debug_)
00154 cout<<"PFClusterAlgo::findSeeds : start"<<endl;
00155 #endif
00156
00157
00158
00159 for(EH ih = eRecHits_.begin(); ih != eRecHits_.end(); ih++ ) {
00160
00161 unsigned rhi = ih->second;
00162
00163 if(! masked(rhi) ) continue;
00164
00165
00166 double rhenergy = ih->first;
00167 const reco::PFRecHit& wannaBeSeed = rechit(rhi, rechits);
00168
00169 if( seedStates_[rhi] == NO ) continue;
00170
00171
00172
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
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
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
00233 if( neighbour.energy() > wannaBeSeed.energy() ) {
00234 seedStates_[rhi] = NO;
00235 break;
00236 }
00237 }
00238
00239 if ( seedStates_[rhi] == YES ) {
00240
00241
00242 seeds_.push_back( rhi );
00243
00244
00245 paint(rhi, SEED);
00246
00247
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
00278
00279
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
00332
00333 cluster.push_back( rhi );
00334
00335
00336 usedInTopo_[ rhi ] = true;
00337
00338
00339
00340
00341 const std::vector< unsigned >& nbs = rh.neighbours4();
00342
00343 for(unsigned i=0; i<nbs.size(); i++) {
00344
00345
00346
00347
00348
00349
00350
00351
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
00377
00378
00379
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
00401
00402 calculateClusterPosition( cluster,
00403 true );
00404
00405
00406
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
00418 seedsintopocluster.push_back( rhi );
00419
00420 }
00421 }
00422
00423
00424
00425
00426 int posCalcNCrystal = seedsintopocluster.size()>1 ? posCalcNCrystal_:-1;
00427
00428
00429
00430 unsigned iter = 0;
00431 unsigned niter = 50;
00432 double diff = 1.;
00433
00434
00435 while ( iter++ < niter && diff > 1E-8 ) {
00436
00437
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
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
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
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
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
00511
00512
00513 math::XYZVector deltav = cposxyzclust;
00514 deltav -= cposxyzcell;
00515 double d = deltav.R() / showerSigma_;
00516
00517
00518
00519
00520
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
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
00546
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
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
00587 }
00588 #endif
00589
00590 continue;
00591 }
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604 if ( dist[ic] < 6. || frac[ic] > 0.99999 ) {
00605
00606 reco::PFRecHitRef recHitRef = createRecHitRef( rechits, rhindex );
00607 reco::PFRecHitFraction rhf( recHitRef,frac[ic] );
00608 curpfclusters[ic].addRecHitFraction( rhf );
00609 }
00610 }
00611
00612
00613 dist.clear();
00614 frac.clear();
00615 }
00616
00617
00618
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
00635
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
00644
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
00673
00674
00675 map <PFLayer::Layer, double> layers;
00676 unsigned seedIndex = 0;
00677 bool seedIndexFound = false;
00678
00679
00680
00681
00682 for (unsigned ic=0; ic<cluster.rechits_.size(); ic++ ) {
00683
00684 unsigned rhi = cluster.rechits_[ic].recHitRef().index();
00685
00686
00687 const reco::PFRecHit& rh = *(cluster.rechits_[ic].recHitRef());
00688 double fraction = cluster.rechits_[ic].fraction();
00689
00690
00691
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
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
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
00725 cluster.setLayer( layer );
00726
00727
00728
00729
00730 double p1 = posCalcP1_;
00731 if( p1 < 0 ) {
00732
00733
00734
00735
00736
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
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
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 ) {
00774 p1 = 1e-9;
00775 }
00776
00777
00778
00779 reco::PFCluster::REPPoint clusterpos;
00780 math::XYZPoint clusterposxyz;
00781 math::XYZPoint firstrechitposxyz;
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
00792
00793 const reco::PFRecHit& rh = *(cluster.rechits_[ic].recHitRef());
00794
00795 if(rhi != seedIndex) {
00796 if( posCalcNCrystal == 5 ) {
00797 if(!rh.isNeighbour4(seedIndex) ) {
00798 continue;
00799 }
00800 }
00801 if( posCalcNCrystal == 9 ) {
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
00825 normalize += norm;
00826 }
00827
00828
00829
00830 if( normalize < 1e-9 ) {
00831
00832
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
00851
00852
00853
00854 if( depcor &&
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
00864
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:
00873 depth = corra * ( corrb + log(cluster.energy_) );
00874 break;
00875 case 2:
00876 depth = corra;
00877 break;
00878 default:
00879 cerr<<"PFClusterAlgo::calculateClusterPosition : unknown function for depth correction! "<<endl;
00880 assert(0);
00881 }
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
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
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
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
00941 math::XYZVector rechitaxis = rh.getAxisXYZ();
00942
00943
00944 math::XYZVector rechitaxisu( rechitaxis );
00945 rechitaxisu /= sqrt( rechitaxis.Mag2() );
00946
00947 math::XYZVector displacement( rechitaxisu );
00948
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
00966 normalize += norm;
00967 }
00968
00969
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 }