CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/Calibration/IsolatedParticles/src/eECALMatrix.cc

Go to the documentation of this file.
00001 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
00002 #include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
00003 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
00004 
00005 #include<iostream>
00006 
00007 namespace spr{
00008 
00009   std::pair<double,bool> eECALmatrix(const DetId& detId, edm::Handle<EcalRecHitCollection>& hitsEB, edm::Handle<EcalRecHitCollection>& hitsEE, const EcalChannelStatus& chStatus, const CaloGeometry* geo, const CaloTopology* caloTopology, int ieta, int iphi, double ebThr, double eeThr, double tMin, double tMax, bool debug) {
00010 
00011     std::vector<DetId> vdets;
00012     spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);
00013 
00014     const EcalRecHitCollection * recHitsEB = 0;
00015     if (hitsEB.isValid())  recHitsEB = hitsEB.product();
00016     bool flag = true;
00017     if (debug) {
00018       std::cout << "Inside eECALmatrix " << 2*ieta+1 << "X" << 2*iphi+1
00019                 << " nXtals " << vdets.size() << std::endl;
00020    }
00021 
00022     double energySum = 0.0;
00023     for (unsigned int i1=0; i1<vdets.size(); i1++) {
00024       if (vdets[i1] != DetId(0)) {
00025         bool ok = true;
00026         std::vector<EcalRecHitCollection::const_iterator> hit;
00027         if (vdets[i1].subdetId()==EcalBarrel) {
00028           spr::findHit(hitsEB,vdets[i1],hit,debug);
00029           ok  = (EcalSeverityLevelAlgo::severityLevel(vdets[i1], (*recHitsEB), chStatus) != EcalSeverityLevelAlgo::kWeird);
00030         } else if (vdets[i1].subdetId()==EcalEndcap) {
00031           spr::findHit(hitsEE,vdets[i1],hit,debug);
00032         }
00033         if (debug) std::cout << "Xtal 0x" <<std::hex << vdets[i1]() <<std::dec;
00034         double ener=0, ethr=ebThr;
00035         if (vdets[i1].subdetId() !=EcalBarrel) ethr = eeThr;
00036         for (unsigned int ihit=0; ihit<hit.size(); ihit++) {
00037           double en=0, tt=0;
00038           if (vdets[i1].subdetId()==EcalBarrel) {
00039             if (hit[ihit] != hitsEB->end()) {
00040               en = hit[ihit]->energy();
00041               tt = hit[ihit]->time();
00042             }
00043           } else if (vdets[i1].subdetId()==EcalEndcap) {
00044             if (hit[ihit] != hitsEE->end()) {
00045               en = hit[ihit]->energy();
00046               tt = hit[ihit]->time();
00047             }
00048           }
00049           if (debug) std::cout << " " << ihit << " " << en;
00050           if (tt > tMin && tt < tMax) ener += en;
00051         }
00052         if (!ok) {
00053           flag = false;
00054           if (debug) std::cout << " detected to be a spike";
00055         }
00056         if (debug) std::cout << "\n";
00057         if (ener > ethr) energySum += ener;
00058       }
00059     }
00060     if (debug) std::cout << "energyECAL: energySum = " << energySum << " flag = " << flag << std::endl;
00061     return std::pair<double,bool>(energySum,flag);
00062   }
00063 
00064   std::pair<double,bool> eECALmatrix(const DetId& detId, edm::Handle<EcalRecHitCollection>& hitsEB, edm::Handle<EcalRecHitCollection>& hitsEE, const EcalChannelStatus& chStatus, const CaloGeometry* geo, const CaloTopology* caloTopology, const EcalTrigTowerConstituentsMap& ttMap, int ieta, int iphi, double ebThr, double eeThr, double tMin, double tMax, bool debug) {
00065 
00066     std::vector<DetId> vdets;
00067     spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, debug);
00068 
00069     const EcalRecHitCollection * recHitsEB = 0;
00070     if (hitsEB.isValid())  recHitsEB = hitsEB.product();
00071     bool flag = true;
00072     if (debug) {
00073       std::cout << "Inside eECALmatrix " << 2*ieta+1 << "X" << 2*iphi+1
00074                 << " nXtals " << vdets.size() << std::endl;
00075    }
00076 
00077     double energySum = 0.0;
00078     for (unsigned int i1=0; i1<vdets.size(); i1++) {
00079       if (vdets[i1] != DetId(0)) {
00080         double eTower = spr::energyECALTower(vdets[i1], hitsEB, hitsEE, ttMap, debug);
00081         bool ok = true;
00082         if      (vdets[i1].subdetId()==EcalBarrel) ok = (eTower > ebThr);
00083         else if (vdets[i1].subdetId()==EcalEndcap) ok = (eTower > eeThr);
00084         if (debug) std::cout << "Crystal 0x" <<std::hex << vdets[i1]() 
00085                              <<std::dec << " Flag " << ok;
00086         if (ok) {
00087           std::vector<EcalRecHitCollection::const_iterator> hit;
00088           if (vdets[i1].subdetId()==EcalBarrel) {
00089             spr::findHit(hitsEB,vdets[i1],hit,debug);
00090             ok  = (EcalSeverityLevelAlgo::severityLevel(vdets[i1], (*recHitsEB), chStatus) != EcalSeverityLevelAlgo::kWeird);
00091           } else if (vdets[i1].subdetId()==EcalEndcap) {
00092             spr::findHit(hitsEE,vdets[i1],hit,debug);
00093           }
00094           double ener=0;
00095           for (unsigned int ihit=0; ihit<hit.size(); ihit++) {
00096             double en=0, tt=0;
00097             if (vdets[i1].subdetId()==EcalBarrel) {
00098               if (hit[ihit] != hitsEB->end()) {
00099                 en = hit[ihit]->energy();
00100                 tt = hit[ihit]->time();
00101               }
00102             } else if (vdets[i1].subdetId()==EcalEndcap) {
00103               if (hit[ihit] != hitsEE->end()) {
00104                 en = hit[ihit]->energy();
00105                 tt = hit[ihit]->time();
00106               }
00107             }
00108             if (debug) std::cout << " " << ihit << " E " << en << " T " << tt;
00109             if (tt > tMin && tt < tMax) ener += en;
00110           }
00111           if (!ok) {
00112             flag = false;
00113             if (debug) std::cout << " detected to be a spike";
00114           }
00115           energySum += ener;
00116         }
00117         if (debug) std::cout << "\n";
00118       }
00119     }
00120     if (debug) std::cout << "energyECAL: energySum = " << energySum << " flag = " << flag << std::endl;
00121     return std::pair<double,bool>(energySum,flag);
00122   }
00123 }