00001 #include "RecoEgamma/EgammaTools/interface/ggPFClusters.h"
00002 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00003 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00004 #include "DataFormats/Math/interface/deltaR.h"
00005 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00006 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
00007 #include "TMath.h"
00008
00009
00010
00011
00012
00013 ggPFClusters::ggPFClusters(
00014
00015 edm::Handle<EcalRecHitCollection>& EBReducedRecHits,
00016 edm::Handle<EcalRecHitCollection>& EEReducedRecHits,
00017 const CaloSubdetectorGeometry* geomBar,
00018 const CaloSubdetectorGeometry* geomEnd
00019
00020 ):
00021 EBReducedRecHits_(EBReducedRecHits),
00022 EEReducedRecHits_(EEReducedRecHits),
00023 geomBar_(geomBar),
00024 geomEnd_(geomEnd)
00025 {
00026
00027 }
00028
00029
00030 ggPFClusters::~ggPFClusters(){
00031
00032 }
00033
00034
00035 vector<reco::CaloCluster> ggPFClusters::getPFClusters(reco::SuperCluster sc){
00036 vector<reco::CaloCluster> PFClust;
00037 reco::CaloCluster_iterator it=sc.clustersBegin();
00038
00039 for(;it!=sc.clustersEnd();++it){
00040 std::vector< std::pair<DetId, float> >bcCells=(*it)->hitsAndFractions();
00041 DetId seedXtalId = bcCells[0].first ;
00042 int detector = seedXtalId.subdetId();
00043 bool isEb;
00044 if(detector==1)isEb=true;
00045 else isEb=false;
00046
00047 float ClusSum=SumPFRecHits(bcCells, isEb);
00048 CaloCluster calo(ClusSum, (*it)->position());
00049 for(unsigned int i=0; i<bcCells.size(); ++i){
00050 calo.addHitAndFraction(bcCells[i].first, bcCells[i].second);
00051 }
00052 PFClust.push_back(calo);
00053 }
00054 return PFClust;
00055 }
00056
00057 float ggPFClusters::SumPFRecHits(std::vector< std::pair<DetId, float> >& bcCells, bool isEB){
00058 float ClustSum=0;
00059 for(unsigned int i=0; i<bcCells.size(); ++i){
00060 EcalRecHitCollection::const_iterator eb;
00061 EcalRecHitCollection::const_iterator ee;
00062
00063 if(isEB){
00064 for(eb=EBReducedRecHits_->begin();eb!= EBReducedRecHits_->end();++eb){
00065 if(eb->id().rawId()==bcCells[i].first.rawId()){
00066 float cellE=eb->energy()* bcCells[i].second;
00067 ClustSum=ClustSum+cellE;
00068 }
00069 }
00070 }
00071 else{
00072 for(ee=EEReducedRecHits_->begin();ee!= EEReducedRecHits_->end();++ee){
00073 if(ee->id().rawId()==bcCells[i].first.rawId()){
00074 float cellE=ee->energy()* bcCells[i].second;
00075 ClustSum=ClustSum+cellE;
00076 break;
00077 }
00078 }
00079 }
00080 }
00081 return ClustSum;
00082 }
00083
00084 float ggPFClusters::getPFSuperclusterOverlap(reco::CaloCluster PFClust, reco::SuperCluster sc){
00085 float overlap=0;
00086 reco::CaloCluster_iterator pit=sc.clustersBegin();
00087 std::vector< std::pair<DetId, float> >bcCellsPF=PFClust.hitsAndFractions();
00088
00089 DetId seedXtalId = bcCellsPF[0].first ;
00090 int detector = seedXtalId.subdetId() ;
00091 bool isEb;
00092 std::vector< std::pair<DetId, float> >bcCellsreco;
00093 if(detector==1)isEb=true;
00094 else isEb=false;
00095 for(;pit!=sc.clustersEnd();++pit){
00096 std::vector< std::pair<DetId, float> >bcCells2=(*pit)->hitsAndFractions();
00097 for(unsigned int h=0; h<bcCells2.size(); ++h)bcCellsreco.push_back(bcCells2[h]);
00098 }
00099 float clustOverlap=0;
00100 clustOverlap=PFRecHitsSCOverlap(
00101 bcCellsPF,
00102 bcCellsreco,
00103 isEb);
00104 overlap=clustOverlap;
00105 return overlap;
00106
00107
00108 }
00109
00110
00111 float ggPFClusters::getPFSuperclusterOverlap(reco::CaloCluster PFClust, reco::Photon phot){
00112 float overlap=0;
00113 SuperClusterRef recoSC=phot.superCluster();
00114 reco::CaloCluster_iterator pit=recoSC->clustersBegin();
00115
00116 std::vector< std::pair<DetId, float> >bcCellsPF=PFClust.hitsAndFractions();
00117
00118 DetId seedXtalId = bcCellsPF[0].first ;
00119 int detector = seedXtalId.subdetId() ;
00120 bool isEb;
00121 std::vector< std::pair<DetId, float> >bcCellsreco;
00122 if(detector==1)isEb=true;
00123 else isEb=false;
00124 for(;pit!=recoSC->clustersEnd();++pit){
00125 std::vector< std::pair<DetId, float> >bcCells2=(*pit)->hitsAndFractions();
00126 for(unsigned int h=0; h<bcCells2.size(); ++h)bcCellsreco.push_back(bcCells2[h]);
00127 }
00128 float clustOverlap=0;
00129 clustOverlap=PFRecHitsSCOverlap(
00130 bcCellsPF,
00131 bcCellsreco,
00132 isEb);
00133 overlap=clustOverlap;
00134 return overlap;
00135 }
00136
00137 float ggPFClusters::PFRecHitsSCOverlap(
00138 std::vector< std::pair<DetId, float> >& bcCells1,
00139 std::vector< std::pair<DetId, float> >& bcCells2,
00140 bool isEB){
00141 float OverlapSum=0;
00142 multimap <DetId, float> pfDID;
00143 multimap <DetId, float> recoDID;
00144 multimap<DetId, float>::iterator pfit;
00145 multimap<DetId, float>::iterator pit;
00146 vector<DetId>matched;
00147 vector<float>matchedfrac;
00148
00149 for(unsigned int i=0; i<bcCells1.size(); ++i){
00150 pfDID.insert(make_pair(bcCells1[i].first, bcCells1[i].second));
00151 }
00152
00153 for(unsigned int i=0; i<bcCells2.size(); ++i){
00154 recoDID.insert(make_pair(bcCells2[i].first, bcCells2[i].second));
00155 }
00156 pit=recoDID.begin();
00157 pfit=pfDID.begin();
00158
00159 for(; pfit!=pfDID.end(); ++pfit){
00160 pit=recoDID.find(pfit->first);
00161 if(pit!=recoDID.end()){
00162
00163 matched.push_back(pfit->first);
00164 matchedfrac.push_back(pfit->second);
00165 }
00166 }
00167
00168 for(unsigned int m=0; m<matched.size(); ++m){
00169 DetId det=matched[m];
00170 EcalRecHitCollection::const_iterator eb;
00171 EcalRecHitCollection::const_iterator ee;
00172 if(isEB){
00173 for(eb=EBReducedRecHits_->begin();eb!= EBReducedRecHits_->end();++eb){
00174 if(eb->id().rawId()==det.rawId()){
00175 float cellE=eb->energy()* matchedfrac[m];
00176 OverlapSum=OverlapSum+cellE;
00177 break;
00178 }
00179 }
00180 }
00181 else{
00182 for(ee=EEReducedRecHits_->begin();ee!= EEReducedRecHits_->end();++ee){
00183 if(ee->id().rawId()==det.rawId()){
00184 float cellE=ee->energy() * matchedfrac[m];
00185 OverlapSum=OverlapSum+cellE;
00186 break;
00187 }
00188 }
00189 }
00190
00191 }
00192 return OverlapSum;
00193 }
00194
00195 void ggPFClusters::localCoordsEB( reco::CaloCluster clus, float &etacry, float &phicry, int &ieta, int &iphi, float &thetatilt, float &phitilt){
00196 const math::XYZPoint position_ = clus.position();
00197 double Theta = -position_.theta()+0.5*TMath::Pi();
00198 double Eta = position_.eta();
00199 double Phi = TVector2::Phi_mpi_pi(position_.phi());
00200
00201
00202
00203 const float X0 = 0.89; const float T0 = 7.4;
00204 double depth = X0 * (T0 + log(clus.energy()));
00205
00206
00207 std::vector< std::pair<DetId, float> > crystals_vector = clus.hitsAndFractions();
00208 float drmin = 999.;
00209 EBDetId crystalseed(crystals_vector[0].first);
00210
00211 for (unsigned int icry=0; icry!=crystals_vector.size(); ++icry) {
00212
00213 EBDetId crystal(crystals_vector[icry].first);
00214
00215 const CaloCellGeometry* cell=geomBar_->getGeometry(crystal);
00216 GlobalPoint center_pos = (dynamic_cast<const TruncatedPyramid*>(cell))->getPosition(depth);
00217 double EtaCentr = center_pos.eta();
00218 double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
00219
00220 float dr = reco::deltaR(Eta,Phi,EtaCentr,PhiCentr);
00221 if (dr<drmin) {
00222 drmin = dr;
00223 crystalseed = crystal;
00224 }
00225
00226 }
00227
00228 ieta = crystalseed.ieta();
00229 iphi = crystalseed.iphi();
00230
00231
00232 const CaloCellGeometry* cell=geomBar_->getGeometry(crystalseed);
00233 const TruncatedPyramid *cpyr = dynamic_cast<const TruncatedPyramid*>(cell);
00234
00235 thetatilt = cpyr->getThetaAxis();
00236 phitilt = cpyr->getPhiAxis();
00237
00238 GlobalPoint center_pos = cpyr->getPosition(depth);
00239
00240 double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
00241 double PhiWidth = (TMath::Pi()/180.);
00242 phicry = (TVector2::Phi_mpi_pi(Phi-PhiCentr))/PhiWidth;
00243
00244 if (ieta<0) phicry *= -1.;
00245
00246 double ThetaCentr = -center_pos.theta()+0.5*TMath::Pi();
00247 double ThetaWidth = (TMath::Pi()/180.)*TMath::Cos(ThetaCentr);
00248 etacry = (Theta-ThetaCentr)/ThetaWidth;
00249
00250
00251 if (ieta<0) etacry *= -1.;
00252 return;
00253
00254 }
00255
00256 void ggPFClusters::localCoordsEE(reco::CaloCluster clus, float &xcry, float &ycry, int &ix, int &iy, float &thetatilt, float &phitilt){
00257 const math::XYZPoint position_ = clus.position();
00258
00259 double Eta = position_.eta();
00260 double Phi = TVector2::Phi_mpi_pi(position_.phi());
00261 double X = position_.x();
00262 double Y = position_.y();
00263
00264
00265
00266 const float X0 = 0.89; float T0 = 1.2;
00267
00268 if (TMath::Abs(clus.eta())<1.653) T0 = 3.1;
00269
00270 double depth = X0 * (T0 + log(clus.energy()));
00271
00272
00273 std::vector< std::pair<DetId, float> > crystals_vector = clus.hitsAndFractions();
00274 float drmin = 999.;
00275 EEDetId crystalseed(crystals_vector[0].first);
00276
00277 for (unsigned int icry=0; icry!=crystals_vector.size(); ++icry) {
00278
00279 EEDetId crystal(crystals_vector[icry].first);
00280
00281 const CaloCellGeometry* cell=geomEnd_->getGeometry(crystal);
00282 GlobalPoint center_pos = (dynamic_cast<const TruncatedPyramid*>(cell))->getPosition(depth);
00283 double EtaCentr = center_pos.eta();
00284 double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
00285
00286 float dr = reco::deltaR(Eta,Phi,EtaCentr,PhiCentr);
00287 if (dr<drmin) {
00288 drmin = dr;
00289 crystalseed = crystal;
00290 }
00291
00292 }
00293
00294 ix = crystalseed.ix();
00295 iy = crystalseed.iy();
00296
00297
00298 const CaloCellGeometry* cell=geomEnd_->getGeometry(crystalseed);
00299 const TruncatedPyramid *cpyr = dynamic_cast<const TruncatedPyramid*>(cell);
00300
00301 thetatilt = cpyr->getThetaAxis();
00302 phitilt = cpyr->getPhiAxis();
00303
00304 GlobalPoint center_pos = cpyr->getPosition(depth);
00305
00306 double XCentr = center_pos.x();
00307 double XWidth = 2.59;
00308 xcry = (X-XCentr)/XWidth;
00309
00310
00311 double YCentr = center_pos.y();
00312 double YWidth = 2.59;
00313 ycry = (Y-YCentr)/YWidth;
00314
00315 }
00316
00317 float ggPFClusters::get5x5Element(int i, int j,
00318 std::vector< std::pair<DetId, float> >&
00319 bcCells,
00320 bool isEB){
00321
00322 Fill5x5Map(bcCells,isEB);
00323 if(abs(i)>2 ||abs(j)>2)return 0;
00324 int ind1=i+2;
00325 int ind2=j+2;
00326 float E=e5x5_[ind1][ind2];
00327 return E;
00328 }
00329
00330 void ggPFClusters::Fill5x5Map(std::vector< std::pair<DetId, float> >& bcCells,
00331 bool isEB){
00332
00333 for(int i=0; i<5; ++i)
00334 for(int j=0; j<5; ++j)e5x5_[i][j]=0;
00335
00336 EcalRecHitCollection::const_iterator eb;
00337 EcalRecHitCollection::const_iterator ee;
00338
00339 DetId idseed=FindSeed(bcCells, isEB);
00340
00341 for(unsigned int i=0; i<bcCells.size(); ++i){
00342 DetId id=bcCells[i].first;
00343 if(isEB){
00344 EBDetId EBidSeed=EBDetId(idseed.rawId());
00345 int deta=EBDetId::distanceEta(id,idseed);
00346 int dphi=EBDetId::distancePhi(id,idseed);
00347 if(abs(dphi)<=2 && abs(deta)<=2){
00348 for(eb=EBReducedRecHits_->begin();eb!= EBReducedRecHits_->end();++eb){
00349 if(eb->id().rawId()==bcCells[i].first.rawId()){
00350 EBDetId EBid=EBDetId(id.rawId());
00351 int ind1=EBidSeed.ieta()-EBid.ieta();
00352 int ind2=EBid.iphi()-EBidSeed.iphi();
00353 if(EBidSeed.ieta() * EBid.ieta() > 0){
00354 ind1=EBid.ieta()-EBidSeed.ieta();
00355 }
00356 else{
00357
00358
00359 ind1 = (EBid.ieta()-EBidSeed.ieta())*(abs(EBid.ieta()-EBidSeed.ieta())-1)/abs(EBid.ieta()-EBidSeed.ieta());
00360 }
00361 int iEta=ind1+2;
00362 int iPhi=ind2+2;
00363 e5x5_[iEta][iPhi]=eb->energy()* bcCells[i].second;
00364 }
00365 }
00366 }
00367 }
00368 else{
00369 EEDetId EEidSeed=EEDetId(idseed.rawId());
00370 int dx=EEDetId::distanceX(id,idseed);
00371 int dy=EEDetId::distanceY(id,idseed);
00372 if(abs(dx)<=2 && abs(dy)<=2){
00373 for(ee=EEReducedRecHits_->begin();ee!= EEReducedRecHits_->end();++ee){
00374 if(ee->id().rawId()==bcCells[i].first.rawId()){
00375 EEDetId EEid=EEDetId(id.rawId());
00376 int ind1=EEid.ix()-EEidSeed.ix();
00377 int ind2=EEid.iy()-EEidSeed.iy();
00378 int ix=ind1+2;
00379 int iy=ind2+2;
00380 e5x5_[ix][iy]=ee->energy()* bcCells[i].second;
00381 }
00382 }
00383 }
00384 }
00385 }
00386 }
00387
00388 DetId ggPFClusters::FindSeed(std::vector< std::pair<DetId, float> >& bcCells, bool isEB){
00389
00390 EcalRecHitCollection::const_iterator eb;
00391 EcalRecHitCollection::const_iterator ee;
00392 DetId idseed=bcCells[0].first;
00393 float PFSeedE=0;
00394
00395 for(unsigned int i=0; i<bcCells.size(); ++i){
00396 if(isEB){
00397 for(eb=EBReducedRecHits_->begin();eb!= EBReducedRecHits_->end();++eb){
00398 if(eb->id().rawId()==bcCells[i].first.rawId()){
00399 float cellE=eb->energy()* bcCells[i].second;
00400 if(PFSeedE<cellE){
00401 PFSeedE=cellE;
00402 idseed=bcCells[i].first;
00403 }
00404 break;
00405 }
00406 }
00407 }
00408 else{
00409 for(ee=EEReducedRecHits_->begin();ee!= EEReducedRecHits_->end();++ee){
00410
00411 if(ee->id().rawId()==bcCells[i].first.rawId()){
00412 float cellE=ee->energy()* bcCells[i].second;
00413 if(PFSeedE<cellE){
00414 PFSeedE=cellE;
00415 idseed=bcCells[i].first;
00416 }
00417 break;
00418 }
00419 }
00420 }
00421 }
00422 return idseed;
00423 }
00424
00425 std::pair<float, float> ggPFClusters::ClusterWidth(
00426 vector<reco::CaloCluster>& PFClust){
00427 pair<float, float> widths(0,0);
00428 multimap<float, unsigned int>ClusterMap;
00429
00430 for(unsigned int i=0; i<PFClust.size();++i)ClusterMap.insert(make_pair(PFClust[i].energy(), i));
00431
00432 std::multimap<float, unsigned int>::reverse_iterator it;
00433 it=ClusterMap.rbegin();
00434 unsigned int max_c=(*it).second;
00435 std::vector< std::pair<DetId, float> >bcCells=PFClust[max_c].hitsAndFractions();
00436 EcalRecHitCollection::const_iterator eb;
00437 EcalRecHitCollection::const_iterator ee;
00438 float numeratorEtaWidth=0;
00439 float numeratorPhiWidth=0;
00440 float denominator=0;
00441
00442 DetId seedXtalId = bcCells[0].first ;
00443 int detector = seedXtalId.subdetId() ;
00444 bool isEb;
00445 if(detector==1)isEb=true;
00446 else isEb=false;
00447
00448
00449 if(isEb){
00450 for(unsigned int i=0; i<bcCells.size(); ++i){
00451 for(eb=EBReducedRecHits_->begin();eb!= EBReducedRecHits_->end();++eb){
00452 if(eb->id().rawId()==bcCells[i].first.rawId()){
00453
00454 double energyHit = eb->energy()*bcCells[i].second;
00455 DetId id=bcCells[i].first;
00456 float eta=geomBar_->getGeometry(id)->getPosition().eta();
00457 float phi=geomBar_->getGeometry(id)->getPosition().phi();
00458 float dEta = eta - PFClust[max_c].eta();
00459 float dPhi = phi - PFClust[max_c].phi();
00460 if (dPhi > + TMath::Pi()) { dPhi = TMath::TwoPi() - dPhi; }
00461 if (dPhi < - TMath::Pi()) { dPhi = TMath::TwoPi() + dPhi; }
00462 numeratorEtaWidth=dEta*dEta*energyHit+numeratorEtaWidth;
00463 numeratorPhiWidth=dPhi*dPhi*energyHit+numeratorPhiWidth;
00464 denominator=energyHit+denominator;
00465 break;
00466 }
00467 }
00468
00469 }
00470 }
00471 else{
00472 for(unsigned int i=0; i<bcCells.size(); ++i){
00473 for(ee=EEReducedRecHits_->begin();ee!= EEReducedRecHits_->end();++ee){
00474 if(ee->id().rawId()==bcCells[i].first.rawId()){
00475 double energyHit = ee->energy()*bcCells[i].second;
00476 DetId id=bcCells[i].first;
00477 float eta=geomEnd_->getGeometry(id)->getPosition().eta();
00478 float phi=geomEnd_->getGeometry(id)->getPosition().phi();
00479 float dEta = eta - PFClust[max_c].eta();
00480 float dPhi = phi - PFClust[max_c].phi();
00481 if (dPhi > + TMath::Pi()) { dPhi = TMath::TwoPi() - dPhi; }
00482 if (dPhi < - TMath::Pi()) { dPhi = TMath::TwoPi() + dPhi; }
00483 numeratorEtaWidth=dEta*dEta*energyHit+numeratorEtaWidth;
00484 numeratorPhiWidth=dPhi*dPhi*energyHit+numeratorPhiWidth;
00485 denominator=energyHit+denominator;
00486 break;
00487 }
00488 }
00489 }
00490 }
00491 float etaWidth=sqrt(numeratorEtaWidth/denominator);
00492 float phiWidth=sqrt(numeratorPhiWidth/denominator);
00493 widths=make_pair(etaWidth, phiWidth);
00494 return widths;
00495 }
00496
00497 double ggPFClusters::LocalEnergyCorrection(const GBRForest *ReaderLCEB, const GBRForest *ReaderLCEE,reco::CaloCluster PFClust,float beamspotZ ){
00498
00499
00500 double PFClustCorr=PFClust.energy();
00501 float ClusEta=PFClust.eta();
00502 float ClusE=PFClust.energy();
00503 float ClusPhi=PFClust.phi();
00504
00505 std::vector<float>inputs;
00506 std::vector< std::pair<DetId, float> >bcCells=PFClust.hitsAndFractions();
00507 bool isEb;
00508 DetId seedXtalId = bcCells[0].first ;
00509 int detector = seedXtalId.subdetId();
00510 if(detector==1)isEb=true;
00511 else isEb=false;
00512
00513 Fill5x5Map(bcCells, isEb);
00514 float Eseed=get5x5Element(0, 0, bcCells, isEb);
00515
00516 float Etop=get5x5Element(0, 1, bcCells, isEb);
00517 float Ebottom=get5x5Element(0, -1, bcCells, isEb);
00518 float Eleft=get5x5Element(-1, 0, bcCells, isEb);
00519 float Eright=get5x5Element(1, 0, bcCells, isEb);
00520 float E5x5=0; float E3x3=0;
00521 for(int i=-2; i<3; ++i)
00522 for(int j=-2; j<3; ++j){
00523 float e=get5x5Element(i, j, bcCells, isEb);
00524 if(abs(i)<2)E3x3=E3x3+e;
00525 E5x5=e+E5x5;
00526 }
00527
00528 float E3x1=e5x5_[2][2]+e5x5_[1][2]+e5x5_[3][2];
00529 float E1x3=e5x5_[2][2]+e5x5_[2][1]+e5x5_[2][3];
00530 float E1x5=e5x5_[2][2]+e5x5_[2][0]+e5x5_[2][1]+e5x5_[2][3]+e5x5_[2][4];
00531
00532 float E2x5Top=e5x5_[0][4]+e5x5_[1][4]+e5x5_[2][4]+e5x5_[3][4]+e5x5_[4][4]
00533 +e5x5_[0][3]+e5x5_[1][3]+e5x5_[2][3]+e5x5_[3][3]+e5x5_[4][3];
00534
00535 float E2x5Bottom=e5x5_[0][0]+e5x5_[1][0]+e5x5_[2][0]+e5x5_[3][0]+e5x5_[4][0]
00536 +e5x5_[0][1]+e5x5_[1][1]+e5x5_[2][1]+e5x5_[3][1]+e5x5_[4][1];
00537
00538 float E2x5Left=e5x5_[0][1]+e5x5_[0][1]+e5x5_[0][2]+e5x5_[0][3]+e5x5_[0][4]
00539 +e5x5_[1][0]+e5x5_[1][1]+e5x5_[1][2]+e5x5_[1][3]+e5x5_[1][4];
00540
00541 float E2x5Right=e5x5_[4][0]+e5x5_[4][1]+e5x5_[4][2]+e5x5_[4][3]+e5x5_[4][4]
00542 +e5x5_[3][0]+e5x5_[3][1]+e5x5_[3][2]+e5x5_[3][3]+e5x5_[3][4];
00543
00544 float centerstrip=e5x5_[2][2]+e5x5_[2][0]+e5x5_[2][1]+e5x5_[2][3]+e5x5_[2][4];
00545 float rightstrip=e5x5_[3][2]+e5x5_[3][0]+e5x5_[3][1]+e5x5_[3][3]+e5x5_[3][4];
00546 float leftstrip=e5x5_[1][2]+e5x5_[1][0]+e5x5_[1][1]+e5x5_[1][3]+e5x5_[1][4];
00547 float E2x5Max=0;
00548 if(rightstrip>leftstrip)E2x5Max=rightstrip+centerstrip;
00549 else E2x5Max=leftstrip+centerstrip;
00550
00551 if(isEb){
00552 float etacry; float phicry; int ieta; int iphi; float thetatilt; float phitilt;
00553 int iEtaCrack=3;
00554 if(abs(ieta)==1 || abs(ieta)==2 )
00555 iEtaCrack=abs(ieta);
00556 if(abs(ieta)>2 && abs(ieta)<24)
00557 iEtaCrack=3;
00558 if(abs(ieta)==24)
00559 iEtaCrack=4;
00560 if(abs(ieta)==25)
00561 iEtaCrack=5;
00562 if(abs(ieta)==26)
00563 iEtaCrack=6;
00564 if(abs(ieta)==27)
00565 iEtaCrack=7;
00566 if(abs(ieta)>27 && abs(ieta)<44)
00567 iEtaCrack=8;
00568 if(abs(ieta)==44)
00569 iEtaCrack=9;
00570 if(abs(ieta)==45)
00571 iEtaCrack=10;
00572 if(abs(ieta)==46)
00573 iEtaCrack=11;
00574 if(abs(ieta)==47)
00575 iEtaCrack=12;
00576 if(abs(ieta)>47 && abs(ieta)<64)
00577 iEtaCrack=13;
00578 if(abs(ieta)==64)
00579 iEtaCrack=14;
00580 if(abs(ieta)==65)
00581 iEtaCrack=15;
00582 if(abs(ieta)==66)
00583 iEtaCrack=16;
00584 if(abs(ieta)==67)
00585 iEtaCrack=17;
00586 if(abs(ieta)>67 && abs(ieta)<84)
00587 iEtaCrack=18;
00588 if(abs(ieta)==84)
00589 iEtaCrack=19;
00590 if(abs(ieta)==85)
00591 iEtaCrack=20;
00592 localCoordsEB(PFClust, etacry, phicry, ieta, iphi, thetatilt, phitilt);
00593 inputs.push_back(beamspotZ);
00594 inputs.push_back(ClusEta/fabs(ClusEta));
00595 inputs.push_back(fabs(ClusEta));
00596 inputs.push_back(fabs(ClusPhi));
00597 inputs.push_back(log(ClusE));
00598 inputs.push_back((Eseed/ClusE));
00599 inputs.push_back((Etop/ClusE));
00600 inputs.push_back((Ebottom/ClusE));
00601 inputs.push_back((Eleft/ClusE));
00602 inputs.push_back((Eright/ClusE));
00603 inputs.push_back(E3x3/ClusE);
00604 inputs.push_back(E1x3/ClusE);
00605 inputs.push_back(E3x1/ClusE);
00606 inputs.push_back(E5x5/ClusE);
00607 inputs.push_back(E1x5/ClusE);
00608 inputs.push_back(E2x5Max/ClusE);
00609 inputs.push_back(E2x5Top/ClusE);
00610 inputs.push_back(E2x5Bottom/ClusE);
00611 inputs.push_back(E2x5Left/ClusE);
00612 inputs.push_back(E2x5Right/ClusE);
00613 inputs.push_back(etacry);
00614 inputs.push_back(phicry);
00615 inputs.push_back(iphi%2);
00616 inputs.push_back(ieta%5);
00617 inputs.push_back(iphi%20);
00618 inputs.push_back(iEtaCrack);
00619 int size=inputs.size();
00620 float PFInputs[26];
00621 for(int i=0; i<size; ++i)PFInputs[i]=inputs[i];
00622 PFClustCorr= ReaderLCEB->GetResponse(PFInputs)*ClusE;
00623 }
00624 else{
00625 float xcry; float ycry; int ix; int iy; float thetatilt; float phitilt;
00626 localCoordsEE(PFClust, xcry, ycry, ix, iy, thetatilt, phitilt);
00627 inputs.push_back(beamspotZ);
00628 inputs.push_back(ClusEta/fabs(ClusEta));
00629 inputs.push_back(fabs(ClusEta));
00630 inputs.push_back(fabs(ClusPhi));
00631 inputs.push_back(log(ClusE));
00632 inputs.push_back((Eseed/ClusE));
00633 inputs.push_back((Etop/ClusE));
00634 inputs.push_back((Ebottom/ClusE));
00635 inputs.push_back((Eleft/ClusE));
00636 inputs.push_back((Eright/ClusE));
00637 inputs.push_back(E3x3/ClusE);
00638 inputs.push_back(E1x3/ClusE);
00639 inputs.push_back(E3x1/ClusE);
00640 inputs.push_back(E5x5/ClusE);
00641 inputs.push_back(E1x5/ClusE);
00642 inputs.push_back(E2x5Max/ClusE);
00643 inputs.push_back(E2x5Top/ClusE);
00644 inputs.push_back(E2x5Bottom/ClusE);
00645 inputs.push_back(E2x5Left/ClusE);
00646 inputs.push_back(E2x5Right/ClusE);
00647 inputs.push_back(xcry);
00648 inputs.push_back(ycry);
00649 int size=inputs.size();
00650 float PFInputs[22];
00651 for(int i=0; i<size; ++i)PFInputs[i]=inputs[i];
00652 PFClustCorr= ReaderLCEE->GetResponse(PFInputs) *ClusE;
00653 }
00654
00655 return PFClustCorr;
00656 }
00657 void ggPFClusters::BasicClusterPFCandLink(
00658 reco::SuperCluster sc,
00659 std::vector<reco::PFCandidatePtr>&insideBox,
00660 std::vector<DetId>& MatchedRH
00661 ){
00662 std::vector<reco::PFCandidatePtr>Linked;
00663 for(unsigned int p=0; p<insideBox.size(); ++p){
00664 math::XYZPointF position_ = insideBox[p]->positionAtECALEntrance();
00665
00666
00667
00668
00669
00670 if(insideBox[p]->pdgId()==22){
00671 double Theta = -position_.theta()+0.5*TMath::Pi();
00672 double Eta = position_.eta();
00673 double Phi = TVector2::Phi_mpi_pi(position_.phi());
00674 double X = position_.x();
00675 double Y = position_.y();
00676 reco::CaloCluster_iterator cit=sc.clustersBegin();
00677 std::vector< std::pair<DetId, float> > crystals_vector = (*cit)->hitsAndFractions();
00678 DetId seedXtalId = crystals_vector[0].first;
00679 int detector = seedXtalId.subdetId();
00680 bool isEb=false;
00681 float X0 = 0.89; float T0 = 7.4;
00682 double depth = X0 * (T0 + log((*cit)->energy()));
00683 if(detector==1){
00684 X0 = 0.89; T0 = 7.4;
00685 depth = X0 * (T0 + log((*cit)->energy()));
00686 isEb=true;
00687 }
00688 else{
00689 X0 = 0.89; T0 = 1.2;
00690 if(fabs(Eta)<1.653)T0=3.1;
00691 depth = X0 * (T0 + log((*cit)->energy()));
00692
00693 isEb=false;
00694 }
00695 crystals_vector.clear();
00696 for(; cit!=sc.clustersEnd(); ++cit){
00697 bool matchBC=false;
00698 crystals_vector = (*cit)->hitsAndFractions();
00699
00700
00701 for (unsigned int icry=0; icry<crystals_vector.size(); ++icry){
00702
00703 if(isEb){
00704
00705 EBDetId crystal(crystals_vector[icry].first);
00706 const CaloCellGeometry* cell=geomBar_->getGeometry(crystal);
00707
00708 GlobalPoint center_pos = (dynamic_cast<const TruncatedPyramid*>(cell))->getPosition(depth);
00709
00710
00711
00712
00713
00714 double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
00715 double PhiWidth = (TMath::Pi()/180.);
00716 double ThetaCentr = -center_pos.theta()+0.5*TMath::Pi();
00717 double ThetaWidth = (TMath::Pi()/180.)*TMath::Cos(ThetaCentr);
00718 double phicry = (TVector2::Phi_mpi_pi(Phi-PhiCentr))/PhiWidth;
00719 double etacry=(Theta-ThetaCentr)/ThetaWidth;
00720 if(fabs(etacry)<0.6 && fabs(phicry)<0.6){
00721 Linked.push_back(insideBox[p]);
00722 matchBC=true;
00723 MatchedRH.push_back(crystals_vector[icry].first);
00724 break;
00725
00726 }
00727
00728 }
00729 else{
00730
00731 EEDetId crystal(crystals_vector[icry].first);
00732 const CaloCellGeometry* cell=geomEnd_->getGeometry(crystal);
00733
00734 GlobalPoint center_pos = (dynamic_cast<const TruncatedPyramid*>(cell))->getPosition(depth);
00735
00736
00737 double XCentr = center_pos.x();
00738 double XWidth = 2.59;
00739 double xcry = (X-XCentr)/XWidth;
00740 double YCentr = center_pos.y();
00741 double YWidth = 2.59;
00742 double ycry = (Y-YCentr)/YWidth;
00743 if(fabs(xcry)<0.6 && fabs(ycry)<0.6){
00744 Linked.push_back(insideBox[p]);
00745 MatchedRH.push_back(crystals_vector[icry].first);
00746 matchBC=true;
00747 break;
00748 }
00749
00750 }
00751
00752 }
00753 if(matchBC)break;
00754 }
00755 }
00756 if(abs(insideBox[p]->pdgId())==211){
00757 float drmin = 999.;
00758 float deta=999;
00759 float dphi=999;
00760 reco::CaloCluster_iterator cit=sc.clustersBegin();
00761 DetId seedCrystal=(*cit)->hitsAndFractions()[0].first;
00762 for(; cit!=sc.clustersEnd(); ++cit){
00763 DetId crystals_vector_seed=(*cit)->hitsAndFractions()[0].first;
00764 math::XYZVector photon_directionWrtVtx((*cit)->position().x()-insideBox[p]->vx(),
00765 (*cit)->position().y()-insideBox[p]->vy(),
00766 (*cit)->position().z()-insideBox[p]->vz()
00767 );
00768 float dR=deltaR(photon_directionWrtVtx.eta(), photon_directionWrtVtx.phi(), position_.eta(), position_.phi());
00769
00770 if(dR<drmin){
00771 drmin=dR;
00772 deta=fabs(photon_directionWrtVtx.eta()-position_.eta());
00773 dphi=acos(cos(photon_directionWrtVtx.phi()- position_.phi()));
00774 seedCrystal=crystals_vector_seed;
00775 }
00776 }
00777 if(deta<0.05 && dphi<0.07){
00778 Linked.push_back(insideBox[p]);
00779 MatchedRH.push_back(seedCrystal);
00780 }
00781 }
00782
00783 }
00784 insideBox.clear();
00785 insideBox=Linked;
00786 }