![]() |
![]() |
00001 #include "DataFormats/EcalDetId/interface/EBDetId.h" 00002 #include "DataFormats/EcalDetId/interface/EEDetId.h" 00003 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h" 00004 #include "RecoParticleFlow/PFClusterTools/interface/PFPhotonClusters.h" 00005 00006 #include <TMath.h> 00007 #include <TVector2.h> 00008 using namespace reco; 00009 PFPhotonClusters::PFPhotonClusters(PFClusterRef PFClusterRef): 00010 PFClusterRef_(PFClusterRef) 00011 { 00012 if(PFClusterRef_->layer()==PFLayer:: ECAL_BARREL )isEB_=true; 00013 else isEB_=false; 00014 SetSeed(); 00015 PFCrystalCoor(); 00016 for(int i=0; i<5; ++i) 00017 for(int j=0; j<5; ++j)e5x5_[i][j]=0; 00018 FillClusterShape(); 00019 FillClusterWidth(); 00020 } 00021 00022 void PFPhotonClusters::SetSeed(){ 00023 double PFSeedE=0; 00024 math::XYZVector axis; 00025 math::XYZVector position; 00026 DetId idseed; 00027 const std::vector< reco::PFRecHitFraction >& PFRecHits= 00028 PFClusterRef_->recHitFractions(); 00029 00030 for(std::vector< reco::PFRecHitFraction >::const_iterator it = PFRecHits.begin(); 00031 it != PFRecHits.end(); ++it){ 00032 const PFRecHitRef& RefPFRecHit = it->recHitRef(); 00033 double frac=it->fraction(); 00034 float E= RefPFRecHit->energy()* frac; 00035 if(E>PFSeedE){ 00036 PFSeedE=E; 00037 axis=RefPFRecHit->getAxisXYZ(); 00038 position=RefPFRecHit->position(); 00039 idseed = RefPFRecHit->detId(); 00040 } 00041 } 00042 idseed_=idseed; 00043 seedPosition_=position; 00044 seedAxis_=axis; 00045 } 00046 00047 void PFPhotonClusters::PFCrystalCoor(){ 00048 if(PFClusterRef_->layer()==PFLayer:: ECAL_BARREL ){//is Barrel 00049 isEB_=true; 00050 EBDetId EBidSeed=EBDetId(idseed_.rawId()); 00051 CrysIEta_=EBidSeed.ieta(); 00052 CrysIPhi_=EBidSeed.iphi(); 00053 double depth = PFClusterRef_->getDepthCorrection(PFClusterRef_->energy(), false, false); 00054 math::XYZVector center_pos = seedPosition_+depth*seedAxis_; 00055 //Crystal Coordinates: 00056 double Pi=TMath::Pi(); 00057 float Phi=PFClusterRef_->position().phi(); 00058 double Theta = -(PFClusterRef_->position().theta())+0.5* Pi; 00059 double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi()); 00060 double PhiWidth = (Pi/180.); 00061 double PhiCry = (TVector2::Phi_mpi_pi(Phi-PhiCentr))/PhiWidth; 00062 double ThetaCentr = -center_pos.theta()+0.5*Pi; 00063 double ThetaWidth = (Pi/180.)*cos(ThetaCentr); 00064 00065 double EtaCry = (Theta-ThetaCentr)/ThetaWidth; 00066 CrysEta_=EtaCry; 00067 CrysPhi_=PhiCry; 00068 00069 if(abs(CrysIEta_)==1 || abs(CrysIEta_)==2 ) 00070 CrysIEtaCrack_=abs(CrysIEta_); 00071 if(abs(CrysIEta_)>2 && abs(CrysIEta_)<24) 00072 CrysIEtaCrack_=3; 00073 if(abs(CrysIEta_)==24) 00074 CrysIEtaCrack_=4; 00075 if(abs(CrysIEta_)==25) 00076 CrysIEtaCrack_=5; 00077 if(abs(CrysIEta_)==26) 00078 CrysIEtaCrack_=6; 00079 if(abs(CrysIEta_)==27) 00080 CrysIEtaCrack_=7; 00081 if(abs(CrysIEta_)>27 && abs(CrysIEta_)<44) 00082 CrysIEtaCrack_=8; 00083 if(abs(CrysIEta_)==44) 00084 CrysIEtaCrack_=9; 00085 if(abs(CrysIEta_)==45) 00086 CrysIEtaCrack_=10; 00087 if(abs(CrysIEta_)==46) 00088 CrysIEtaCrack_=11; 00089 if(abs(CrysIEta_)==47) 00090 CrysIEtaCrack_=12; 00091 if(abs(CrysIEta_)>47 && abs(CrysIEta_)<64) 00092 CrysIEtaCrack_=13; 00093 if(abs(CrysIEta_)==64) 00094 CrysIEtaCrack_=14; 00095 if(abs(CrysIEta_)==65) 00096 CrysIEtaCrack_=15; 00097 if(abs(CrysIEta_)==66) 00098 CrysIEtaCrack_=16; 00099 if(abs(CrysIEta_)==67) 00100 CrysIEtaCrack_=17; 00101 if(abs(CrysIEta_)>67 && abs(CrysIEta_)<84) 00102 CrysIEtaCrack_=18; 00103 if(abs(CrysIEta_)==84) 00104 CrysIEtaCrack_=19; 00105 if(abs(CrysIEta_)==85) 00106 CrysIEtaCrack_=20; 00107 } 00108 else{ 00109 isEB_=false; 00110 EEDetId EEidSeed=EEDetId(idseed_.rawId()); 00111 CrysIX_=EEidSeed.ix(); 00112 CrysIY_=EEidSeed.iy(); 00113 float X0 = 0.89; float T0 = 1.2; 00114 if(fabs(PFClusterRef_->eta())>1.653)T0=3.1; 00115 double depth = X0 * (T0 + log(PFClusterRef_->energy())); 00116 math::XYZVector center_pos=(seedPosition_)+depth*seedAxis_; 00117 double XCentr = center_pos.x(); 00118 double YCentr = center_pos.y(); 00119 double XWidth = 2.59; 00120 double YWidth = 2.59; 00121 00122 CrysX_=(PFClusterRef_->x()-XCentr)/XWidth; 00123 CrysY_=(PFClusterRef_->y()-YCentr)/YWidth; 00124 } 00125 } 00126 00127 void PFPhotonClusters::FillClusterShape(){ 00128 const std::vector< reco::PFRecHitFraction >& PFRecHits=PFClusterRef_->recHitFractions(); 00129 for(std::vector< reco::PFRecHitFraction >::const_iterator it = PFRecHits.begin(); it != PFRecHits.end(); ++it){ 00130 const PFRecHitRef& RefPFRecHit = it->recHitRef(); 00131 DetId id=RefPFRecHit->detId(); 00132 double frac=it->fraction(); 00133 float E=RefPFRecHit->energy()*frac; 00134 if(isEB_){ 00135 int deta=EBDetId::distanceEta(id,idseed_); 00136 int dphi=EBDetId::distancePhi(id,idseed_); 00137 00138 if(abs(deta)>2 ||abs(dphi)>2)continue; 00139 00140 //f(abs(dphi)<=2 && abs(deta)<=2){ 00141 EBDetId EBidSeed=EBDetId(idseed_.rawId()); 00142 EBDetId EBid=EBDetId(id.rawId()); 00143 int ind1=EBidSeed.ieta()-EBid.ieta(); 00144 int ind2=EBidSeed.iphi()-EBid.iphi(); 00145 if(EBidSeed.ieta() * EBid.ieta() > 0){ 00146 ind1=EBid.ieta()-EBidSeed.ieta(); 00147 } 00148 else{ //near EB+ EB- 00149 int shift(EBidSeed.ieta()>0 ? -1 : 1); 00150 ind1=EBidSeed.ieta()-EBid.ieta()+shift; 00151 } 00152 00153 // more tricky edges in phi. Note that distance is already <2 at this point 00154 if( (EBidSeed.iphi()<5&&EBid.iphi()>355) || (EBidSeed.iphi()>355&&EBid.iphi()<5)) { 00155 int shift(EBidSeed.iphi()<180 ? EBDetId::MAX_IPHI:-EBDetId::MAX_IPHI) ; 00156 ind2 = shift + EBidSeed.iphi() - EBid.iphi(); 00157 // std::cout << " Phi check " << EBidSeed.iphi() << " " << EBid.iphi() << " " << ind2 << std::endl; 00158 } 00159 00160 int iEta=ind1+2; 00161 int iPhi=ind2+2; 00162 //std::cout<<"IEta, IPhi "<<iEta<<", "<<iPhi<<std::endl; 00163 // if(iPhi >= 5 || iPhi <0) { 00164 // std::cout << "dphi "<< EBDetId::distancePhi(id,idseed_) << " iphi " << EBid.iphi() << " iphiseed " << EBidSeed.iphi() << " iPhi " << iPhi << std::endl;} 00165 // if(iEta >= 5 || iEta <0) { 00166 // std::cout << "deta "<< EBDetId::distanceEta(id,idseed_) << " ieta " << EBid.ieta() << " ietaseed " << EBidSeed.ieta() << "ind1 " << ind1 << " iEta " << iEta << " " ; 00167 // ind1=ind1prime; 00168 // iEta=ind1+2; 00169 // std::cout << " new iEta " << iEta << std::endl; 00170 // } 00171 // assert(iEta < 5); 00172 // assert(iEta >= 0); 00173 // assert(iPhi < 5); 00174 // assert(iPhi >= 0); 00175 e5x5_[iEta][iPhi]=E; 00176 } 00177 else{ 00178 int dx=EEDetId::distanceX(id,idseed_); 00179 int dy=EEDetId::distanceY(id,idseed_); 00180 if(abs(dx)>2 ||abs(dy>2))continue; 00181 EEDetId EEidSeed=EEDetId(idseed_.rawId()); 00182 EEDetId EEid=EEDetId(id.rawId()); 00183 int ind1=EEid.ix()-EEidSeed.ix(); 00184 int ind2=EEid.iy()-EEidSeed.iy(); 00185 int ix=ind1+2; 00186 int iy=ind2+2; 00187 //std::cout<<"IX, IY "<<ix<<", "<<iy<<std::endl; 00188 // assert(ix < 5); 00189 // assert(ix >= 0); 00190 // assert(iy < 5); 00191 // assert(iy >= 0); 00192 e5x5_[ix][iy]=E; 00193 } 00194 } 00195 } 00196 00197 void PFPhotonClusters::FillClusterWidth(){ 00198 double numeratorEtaWidth = 0.; 00199 double numeratorPhiWidth = 0.; 00200 double numeratorEtaPhiWidth = 0.; 00201 double ClustEta=PFClusterRef_->eta(); 00202 double ClustPhi=PFClusterRef_->phi(); 00203 const std::vector< reco::PFRecHitFraction >& PFRecHits=PFClusterRef_->recHitFractions(); 00204 for(std::vector< reco::PFRecHitFraction >::const_iterator it = PFRecHits.begin(); it != PFRecHits.end(); ++it){ 00205 const PFRecHitRef& RefPFRecHit = it->recHitRef(); 00206 float E=RefPFRecHit->energy() * it->fraction(); 00207 double dEta = RefPFRecHit->position().eta() - ClustEta; 00208 double dPhi = RefPFRecHit->position().phi() - ClustPhi; 00209 if (dPhi > + TMath::Pi()) { dPhi = TMath::TwoPi() - dPhi; } 00210 if (dPhi < - TMath::Pi()) { dPhi = TMath::TwoPi() + dPhi; } 00211 numeratorEtaWidth += E * dEta * dEta; 00212 numeratorPhiWidth += E * dPhi * dPhi; 00213 numeratorEtaPhiWidth += E * fabs(dPhi) * fabs(dEta); 00214 } 00215 double denominator=PFClusterRef_->energy(); 00216 sigetaeta_ = sqrt(numeratorEtaWidth / denominator); 00217 sigphiphi_ = sqrt(numeratorPhiWidth / denominator); 00218 sigetaphi_ = sqrt(numeratorEtaPhiWidth / denominator); 00219 }