Reference Manual

Go to the documentation of this file.
00001 // ROOT includes
00002 #include <Math/VectorUtil.h>
00003 
00004 #include "RecoHI/HiEgammaAlgos/interface/CxCalculator.h"
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006 
00007 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00008 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00009 
00010 #include "DataFormats/Common/interface/Handle.h"
00011 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00012 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00013 #include "DataFormats/Math/interface/Vector3D.h"
00014 
00015 
00016 using namespace edm;
00017 using namespace reco;
00018 using namespace std;
00019 using namespace ROOT::Math::VectorUtil; 
00020 
00021 #define PI 3.141592653589793238462643383279502884197169399375105820974945
00022 
00023 CxCalculator::CxCalculator (const edm::Event &iEvent, const edm::EventSetup &iSetup, edm::InputTag barrelLabel, edm::InputTag endcapLabel)
00024 {
00025 //InputTag("islandBasicClusters:islandBarrelBasicClusters")
00026 //InputTag("islandBasicClusters:islandEndcapBasicClusters")
00027    Handle<BasicClusterCollection> pEBclusters;
00028    iEvent.getByLabel(barrelLabel, pEBclusters);
00029    fEBclusters_ = pEBclusters.product(); 
00030 
00031    Handle<BasicClusterCollection> pEEclusters;
00032    iEvent.getByLabel(endcapLabel, pEEclusters);
00033    fEEclusters_ = pEEclusters.product(); 
00034 
00035    ESHandle<CaloGeometry> geometryHandle;
00036    iSetup.get<CaloGeometryRecord>().get(geometryHandle);
00037 
00038    geometry_ = geometryHandle.product();
00039 } 
00040 
00041 double CxCalculator::getCx(const reco::SuperClusterRef cluster, double x, double threshold)
00042 {
00043    using namespace edm;
00044    using namespace reco;
00045 
00046    if(!fEBclusters_) {       
00047       LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
00048       return -100;
00049    }
00050 
00051    if(!fEEclusters_) {       
00052       LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
00053       return -100;
00054    }
00055 
00056    math::XYZVector SClusPoint(cluster->position().x(),
00057                           cluster->position().y(),
00058                           cluster->position().z());
00059 
00060    double TotalEt = 0;
00061 
00062    TotalEt = - cluster->rawEnergy()/cosh(cluster->eta());
00063 
00064    // Loop over barrel basic clusters 
00065    for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
00066        iclu != fEBclusters_->end(); ++iclu) {
00067       const BasicCluster *clu = &(*iclu);
00068       math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
00069       double eta = ClusPoint.eta();
00070 
00071       double dR = ROOT::Math::VectorUtil::DeltaR(ClusPoint,SClusPoint);
00072       
00073       if (dR<x*0.1) {
00074          double et = clu->energy()/cosh(eta);
00075          if (et<threshold) et=0;
00076          TotalEt += et;
00077       } 
00078    }
00079 
00080    for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
00081        iclu != fEEclusters_->end(); ++iclu) {
00082       const BasicCluster *clu = &(*iclu);
00083       const GlobalPoint clusPoint(clu->x(),clu->y(),clu->z());
00084       math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
00085       double eta = ClusPoint.eta();
00086 
00087       double dR = ROOT::Math::VectorUtil::DeltaR(ClusPoint,SClusPoint);
00088 
00089       if (dR<x*0.1) {
00090          double et = clu->energy()/cosh(eta);
00091          if (et<threshold) et=0;
00092          TotalEt += et;
00093       } 
00094    }
00095 
00096    return TotalEt;
00097 }
00098 
00099 double CxCalculator::getCxRemoveSC(const reco::SuperClusterRef cluster, double x, double threshold)
00100 {
00101    // Calculate Cx and remove the basicClusters used by superCluster
00102 
00103    using namespace edm;
00104    using namespace reco;
00105 
00106    if(!fEBclusters_) {       
00107       LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
00108       return -100;
00109    }
00110 
00111    if(!fEEclusters_) {       
00112       LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
00113       return -100;
00114    }
00115 
00116    math::XYZVector SClusPoint(cluster->position().x(),
00117                           cluster->position().y(),
00118                           cluster->position().z());
00119 
00120    double TotalEt = 0;
00121 
00122    TotalEt = 0;
00123 
00124    // Loop over barrel basic clusters 
00125    for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
00126        iclu != fEBclusters_->end(); ++iclu) {
00127       const BasicCluster *clu = &(*iclu);
00128       math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
00129       double eta = ClusPoint.eta();
00130 
00131       double dR = ROOT::Math::VectorUtil::DeltaR(ClusPoint,SClusPoint);
00132       
00133       // check if this basic cluster is used in the target supercluster
00134       bool inSuperCluster = checkUsed(cluster,clu);
00135 
00136       if (dR<x*0.1&&inSuperCluster==false) {
00137          double et = clu->energy()/cosh(eta);
00138          if (et<threshold) et=0;
00139          TotalEt += et;
00140       } 
00141      
00142    }
00143 
00144    for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
00145        iclu != fEEclusters_->end(); ++iclu) {
00146       const BasicCluster *clu = &(*iclu);
00147       const GlobalPoint clusPoint(clu->x(),clu->y(),clu->z());
00148       math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
00149       double eta = ClusPoint.eta();
00150 
00151       double dR = ROOT::Math::VectorUtil::DeltaR(ClusPoint,SClusPoint);
00152 
00153       // check if this basic cluster is used in the target supercluster
00154       bool inSuperCluster = checkUsed(cluster,clu);
00155 
00156       if (dR<x*0.1&&inSuperCluster==false) {
00157          double et = clu->energy()/cosh(eta);
00158          if (et<threshold) et=0;
00159          TotalEt += et;
00160       } 
00161    }
00162 
00163    return TotalEt;
00164 }
00165 
00166 double CxCalculator::getCCx(const reco::SuperClusterRef cluster, double x, double threshold)
00167 {
00168    using namespace edm;
00169    using namespace reco;
00170 
00171 
00172    if(!fEBclusters_) {       
00173       LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
00174       return -100;
00175    }
00176 
00177    if(!fEEclusters_) {       
00178       LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
00179       return -100;
00180    }
00181 
00182    double SClusterEta = cluster->eta();
00183    double SClusterPhi = cluster->phi();
00184    double TotalEt = 0;
00185 
00186    TotalEt = - cluster->rawEnergy()/cosh(cluster->eta());
00187 
00188    for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
00189        iclu != fEBclusters_->end(); ++iclu) {
00190       const BasicCluster *clu = &(*iclu);
00191       math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
00192       double eta = ClusPoint.eta();
00193 
00194       double dEta = fabs(eta-SClusterEta);
00195  
00196      if (dEta<x*0.1) {
00197          double et = clu->energy()/cosh(eta);
00198          if (et<threshold) et=0;
00199          TotalEt += et;
00200       } 
00201    }
00202 
00203    for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
00204        iclu != fEEclusters_->end(); ++iclu) {
00205       const BasicCluster *clu = &(*iclu);
00206       math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
00207       double eta = ClusPoint.eta();
00208       double phi = ClusPoint.phi();
00209 
00210       double dEta = fabs(eta-SClusterEta);
00211       double dPhi = fabs(phi-SClusterPhi);
00212       while (dPhi>2*PI) dPhi-=2*PI;
00213       if (dPhi>PI) dPhi=2*PI-dPhi;
00214 
00215       if (dEta<x*0.1) {
00216          double et = clu->energy()/cosh(eta);
00217          if (et<threshold) et=0;
00218          TotalEt += et;
00219       } 
00220    }
00221 
00222    double Cx = getCx(cluster,x,threshold);
00223    double CCx = Cx - TotalEt / 40.0 * x;
00224 
00225    return CCx;
00226 }
00227 
00228 
00229 double CxCalculator::getCCxRemoveSC(const reco::SuperClusterRef cluster, double x, double threshold)
00230 {
00231    using namespace edm;
00232    using namespace reco;
00233 
00234 
00235    if(!fEBclusters_) {       
00236       LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
00237       return -100;
00238    }
00239 
00240    if(!fEEclusters_) {       
00241       LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
00242       return -100;
00243    }
00244 
00245    double SClusterEta = cluster->eta();
00246    double SClusterPhi = cluster->phi();
00247    double TotalEt = 0;
00248 
00249    TotalEt = 0;
00250 
00251    for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
00252        iclu != fEBclusters_->end(); ++iclu) {
00253       const BasicCluster *clu = &(*iclu);
00254       math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
00255       double eta = ClusPoint.eta();
00256 
00257       double dEta = fabs(eta-SClusterEta);
00258 
00259       // check if this basic cluster is used in the target supercluster
00260       bool inSuperCluster = checkUsed(cluster,clu);
00261  
00262      if (dEta<x*0.1&&inSuperCluster==false) {
00263          double et = clu->energy()/cosh(eta);
00264          if (et<threshold) et=0;
00265          TotalEt += et;
00266       } 
00267    }
00268 
00269    for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
00270        iclu != fEEclusters_->end(); ++iclu) {
00271       const BasicCluster *clu = &(*iclu);
00272       math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
00273       double eta = ClusPoint.eta();
00274       double phi = ClusPoint.phi();
00275 
00276       double dEta = fabs(eta-SClusterEta);
00277       double dPhi = fabs(phi-SClusterPhi);
00278       while (dPhi>2*PI) dPhi-=2*PI;
00279       if (dPhi>PI) dPhi=2*PI-dPhi;
00280 
00281       // check if this basic cluster is used in the target supercluster
00282       bool inSuperCluster = checkUsed(cluster,clu);
00283 
00284       if (dEta<x*0.1&&inSuperCluster==false) {
00285          double et = clu->energy()/cosh(eta);
00286          if (et<threshold) et=0;
00287          TotalEt += et;
00288       } 
00289    }
00290 
00291    double Cx = getCxRemoveSC(cluster,x,threshold);
00292    double CCx = Cx - TotalEt / 40.0 * x;
00293 
00294    return CCx;
00295 }
00296 
00297 
00298 bool CxCalculator::checkUsed(const reco::SuperClusterRef sc, const reco::BasicCluster* bc)
00299 {
00300    reco::CaloCluster_iterator theEclust = sc->clustersBegin();
00301 
00302    // Loop over the basicClusters inside the target superCluster
00303    for(;theEclust != sc->clustersEnd(); theEclust++) {
00304      if ((**theEclust) == (*bc) ) return  true; //matched, so it's used.
00305    }
00306    return false;
00307 }
00308 
00309 double CxCalculator::getBCMax(const reco::SuperClusterRef cluster,int i)
00310 {
00311    reco::CaloCluster_iterator theEclust = cluster->clustersBegin();
00312 
00313    double energyMax=0,energySecond=0;
00314    // Loop over the basicClusters inside the target superCluster
00315    for(;theEclust != cluster->clustersEnd(); theEclust++) {
00316      if ((*theEclust)->energy()>energyMax ) {
00317         energySecond=energyMax;
00318         energyMax=(*theEclust)->energy();
00319      } else if ((*theEclust)->energy()>energySecond) {
00320         energySecond=(*theEclust)->energy();
00321      }
00322    }
00323    if (i==1) return energyMax;
00324    return energySecond;
00325 }
00326 
00327 
00328 double CxCalculator::getCorrection(const reco::SuperClusterRef cluster, double x, double y,double threshold)
00329 {
00330    using namespace edm;
00331    using namespace reco;
00332 
00333    // doesn't really work now ^^; (Yen-Jie)
00334    if(!fEBclusters_) {       
00335       LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
00336       return -100;
00337    }
00338 
00339    if(!fEEclusters_) {       
00340       LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
00341       return -100;
00342    }
00343 
00344    double SClusterEta = cluster->eta();
00345    double SClusterPhi = cluster->phi();
00346    double TotalEnergy = 0;
00347    double TotalBC = 0;
00348 
00349    TotalEnergy = 0;
00350 
00351    double Area = PI * (-x*x+y*y) / 100.0;
00352    double nCrystal = Area / 0.0174 / 0.0174; // ignore the difference between endcap and barrel for the moment....
00353 
00354    for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
00355        iclu != fEBclusters_->end(); ++iclu) {
00356       const BasicCluster *clu = &(*iclu);
00357       const GlobalPoint clusPoint(clu->x(),clu->y(),clu->z());
00358       double eta = clusPoint.eta();
00359       double phi = clusPoint.phi();
00360       double dEta = fabs(eta-SClusterEta);
00361       double dPhi = fabs(phi-SClusterPhi);
00362       while (dPhi>2*PI) dPhi-=2*PI;
00363       if (dPhi>PI) dPhi=2*PI-dPhi;
00364       double dR = sqrt(dEta*dEta+dPhi*dPhi);
00365  
00366      if (dR>x*0.1&&dR<y*0.1) {
00367          double e = clu->energy();
00368          if (e<threshold) e=0;
00369          TotalEnergy += e;
00370          if (e!=0) TotalBC+=clu->size();  // number of crystals
00371    
00372       } 
00373    }
00374 
00375    for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
00376        iclu != fEEclusters_->end(); ++iclu) {
00377       const BasicCluster *clu = &(*iclu);
00378       const GlobalPoint clusPoint(clu->x(),clu->y(),clu->z());
00379       double eta = clusPoint.eta();
00380       double phi = clusPoint.phi();
00381       double dEta = fabs(eta-SClusterEta);
00382       double dPhi = fabs(phi-SClusterPhi);
00383       while (dPhi>2*PI) dPhi-=2*PI;
00384       if (dPhi>PI) dPhi=2*PI-dPhi;
00385       double dR = sqrt(dEta*dEta+dPhi*dPhi);
00386  
00387      if (dR>x*0.1&&dR<y*0.1) {
00388          double e = clu->energy();
00389          if (e<threshold) e=0;
00390          TotalEnergy += e;
00391          if (e!=0) TotalBC += clu->size(); // number of crystals
00392       } 
00393    }
00394 
00395 
00396   if (TotalBC==0) return 0;
00397   return TotalEnergy/nCrystal;
00398 }
00399 
00400 double CxCalculator::getAvgBCEt(const reco::SuperClusterRef cluster, double x,double phi1,double phi2, double threshold)
00401 // x: eta cut, phi1: deltaPhiMin cut, phi2: deltaPhiMax
00402 {
00403    using namespace edm;
00404    using namespace reco;
00405 
00406 
00407    if(!fEBclusters_) {       
00408       LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
00409       return -100;
00410    }
00411 
00412    if(!fEEclusters_) {       
00413       LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
00414       return -100;
00415    }
00416 
00417    double SClusterEta = cluster->eta();
00418    double SClusterPhi = cluster->phi();
00419 
00420    double TotalEt = 0;    // Total E
00421    double TotalN = 0;     // Total N
00422 
00423    TotalEt = - cluster->rawEnergy()/cosh(cluster->eta());
00424 
00425    if (fabs(SClusterEta) < 1.479) {
00426       //Barrel    
00427       for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
00428           iclu != fEBclusters_->end(); ++iclu) {
00429          const BasicCluster *clu = &(*iclu);
00430          math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
00431          double eta = ClusPoint.eta();
00432          double phi = ClusPoint.phi();  
00433 
00434          double dEta = fabs(eta-SClusterEta);
00435          double dPhi = fabs(phi-SClusterPhi);
00436          while (dPhi>2*PI) dPhi-=2*PI;
00437 
00438          bool inSuperCluster = checkUsed(cluster,clu);
00439 
00440          if (dEta<x*0.1&&inSuperCluster==false&&dPhi>phi1*0.1&&dPhi<phi2*0.1) {
00441             double et = clu->energy()/cosh(eta);
00442             if (et<threshold) et=0;
00443             TotalEt += et;
00444             TotalN ++;
00445          } 
00446       }   
00447    } else {
00448       //Endcap
00449       for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
00450           iclu != fEEclusters_->end(); ++iclu) {
00451          const BasicCluster *clu = &(*iclu);
00452          math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
00453          double eta = ClusPoint.eta();
00454          double phi = ClusPoint.phi();  
00455    
00456          double dEta = fabs(eta-SClusterEta);
00457          double dPhi = fabs(phi-SClusterPhi);
00458          while (dPhi>2*PI) dPhi-=2*PI;
00459 
00460          bool inSuperCluster = checkUsed(cluster,clu);
00461 
00462          if (dEta<x*0.1&&inSuperCluster==false&&dPhi>phi1*0.1&&dPhi<phi2*0.1) {
00463             double et = clu->energy()/cosh(eta);
00464             if (et<threshold) et=0;
00465             TotalEt += et;
00466             TotalN ++;
00467          } 
00468       }
00469    }
00470    return TotalEt / TotalN;
00471 }
00472 
00473 double CxCalculator::getNBC(const reco::SuperClusterRef cluster, double x,double phi1,double phi2, double threshold)
00474 // x: eta cut, phi1: deltaPhiMin cut, phi2: deltaPhiMax
00475 {
00476    using namespace edm;
00477    using namespace reco;
00478 
00479 
00480    if(!fEBclusters_) {       
00481       LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
00482       return -100;
00483    }
00484 
00485    if(!fEEclusters_) {       
00486       LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
00487       return -100;
00488    }
00489 
00490    double SClusterEta = cluster->eta();
00491    double SClusterPhi = cluster->phi();
00492 
00493    double TotalEt = 0;    // Total E
00494    double TotalN = 0;     // Total N
00495 
00496    TotalEt = - cluster->rawEnergy()/cosh(cluster->eta());
00497 
00498    
00499 
00500    if (fabs(SClusterEta) < 1.479) {
00501       //Barrel    
00502       for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
00503           iclu != fEBclusters_->end(); ++iclu) {
00504          const BasicCluster *clu = &(*iclu);
00505          math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
00506          double eta = ClusPoint.eta();
00507          double phi = ClusPoint.phi();  
00508 
00509          double dEta = fabs(eta-SClusterEta);
00510          double dPhi = fabs(phi-SClusterPhi);
00511          while (dPhi>2*PI) dPhi-=2*PI;
00512 
00513          bool inSuperCluster = checkUsed(cluster,clu);
00514 
00515          if (dEta<x*0.1&&inSuperCluster==false&&dPhi>phi1*0.1&&dPhi<phi2*0.1) {
00516             double et = clu->energy()/cosh(eta);
00517             if (et<threshold) et=0;
00518             TotalEt += et;
00519             TotalN ++;
00520          } 
00521       }   
00522    } else {
00523       //Endcap
00524       for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
00525           iclu != fEEclusters_->end(); ++iclu) {
00526          const BasicCluster *clu = &(*iclu);
00527          math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
00528          double eta = ClusPoint.eta();
00529          double phi = ClusPoint.phi();  
00530    
00531          double dEta = fabs(eta-SClusterEta);
00532          double dPhi = fabs(phi-SClusterPhi);
00533          while (dPhi>2*PI) dPhi-=2*PI;
00534 
00535          bool inSuperCluster = checkUsed(cluster,clu);
00536 
00537          if (dEta<x*0.1&&inSuperCluster==false&&dPhi>phi1*0.1&&dPhi<phi2*0.1) {
00538             double et = clu->energy()/cosh(eta);
00539             if (et<threshold) et=0;
00540             TotalEt += et;
00541             TotalN ++;
00542          } 
00543       }
00544    }
00545    return TotalN;
00546 }
/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoHI/HiEgammaAlgos/src/CxCalculator.cc
