CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_2_SLHC4_patch1/src/RecoParticleFlow/PFClusterTools/src/PFPhotonClusters.cc

Go to the documentation of this file.
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 }