Tools/src/CalibElectron.cc

0001 #include "Calibration/Tools/interface/CalibElectron.h"
0002 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0003 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0004 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0005 #include "DataFormats/EgammaCandidates/interface/Electron.h"
0006 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0007 #include "Calibration/Tools/interface/EcalRingCalibrationTools.h"
0008 #include "Calibration/Tools/interface/EcalIndexingTools.h"
0009 #include <iostream>
0010
0011 using namespace calib;
0012 using namespace std;
0013
0014 CalibElectron::CalibElectron() : theElectron_(nullptr), theHits_(nullptr), theEEHits_(nullptr) {}
0015
0016 std::vector<std::pair<int, float> > CalibElectron::getCalibModulesWeights(TString calibtype) {
0017   std::vector<std::pair<int, float> > theWeights;
0018
0019   if (calibtype == "RING") {
0020     float w_ring[EcalRingCalibrationTools::N_RING_TOTAL];
0021
0022     for (int i = 0; i < EcalRingCalibrationTools::N_RING_TOTAL; ++i)
0023       w_ring[i] = 0.;
0024
0025     std::vector<std::pair<DetId, float> > scDetIds = theElectron_->superCluster()->hitsAndFractions();
0026
0027     for (std::vector<std::pair<DetId, float> >::const_iterator idIt = scDetIds.begin(); idIt != scDetIds.end();
0028          ++idIt) {
0029       const EcalRecHit* rh = nullptr;
0030       if ((*idIt).first.subdetId() == EcalBarrel)
0031         rh = &*(theHits_->find((*idIt).first));
0032       else if ((*idIt).first.subdetId() == EcalEndcap)
0033         rh = &*(theEEHits_->find((*idIt).first));
0034       if (!rh)
0035         std::cout << "CalibElectron::BIG ERROR::RecHit NOT FOUND" << std::endl;
0036       w_ring[EcalRingCalibrationTools::getRingIndex((*idIt).first)] += rh->energy();
0037     }
0038
0039     for (int i = 0; i < EcalRingCalibrationTools::N_RING_TOTAL; ++i)
0040       if (w_ring[i] != 0.)
0041         theWeights.push_back(std::pair<int, float>(i, w_ring[i]));
0042     // std::cout << " ring " << i << " - energy sum " << w_ring[i] << std::endl;
0043   }
0044
0045   else if (calibtype == "MODULE") {
0046     float w_ring[EcalRingCalibrationTools::N_MODULES_BARREL];
0047
0048     for (int i = 0; i < EcalRingCalibrationTools::N_MODULES_BARREL; ++i)
0049       w_ring[i] = 0.;
0050
0051     std::vector<std::pair<DetId, float> > scDetIds = theElectron_->superCluster()->hitsAndFractions();
0052
0053     for (std::vector<std::pair<DetId, float> >::const_iterator idIt = scDetIds.begin(); idIt != scDetIds.end();
0054          ++idIt) {
0055       const EcalRecHit* rh = nullptr;
0056       if ((*idIt).first.subdetId() == EcalBarrel)
0057         rh = &*(theHits_->find((*idIt).first));
0058       else if ((*idIt).first.subdetId() == EcalEndcap)
0059         rh = &*(theEEHits_->find((*idIt).first));
0060       if (!rh)
0061         std::cout << "CalibElectron::BIG ERROR::RecHit NOT FOUND" << std::endl;
0062       w_ring[EcalRingCalibrationTools::getModuleIndex((*idIt).first)] += rh->energy();
0063     }
0064
0065     for (int i = 0; i < EcalRingCalibrationTools::N_MODULES_BARREL; ++i)
0066       if (w_ring[i] != 0.)
0067         theWeights.push_back(std::pair<int, float>(i, w_ring[i]));
0068     // std::cout << " ring " << i << " - energy sum " << w_ring[i] << std::endl;
0069
0070   }
0071
0072   else if (calibtype == "ABS_SCALE") {
0073     std::cout << "ENTERING CalibElectron, ABS SCALE mode" << std::endl;
0074
0075     float w_ring(0.);
0076
0077     std::vector<std::pair<DetId, float> > scDetIds = theElectron_->superCluster()->hitsAndFractions();
0078
0079     for (std::vector<std::pair<DetId, float> >::const_iterator idIt = scDetIds.begin(); idIt != scDetIds.end();
0080          ++idIt) {
0081       const EcalRecHit* rh = nullptr;
0082       if ((*idIt).first.subdetId() == EcalBarrel)
0083         rh = &*(theHits_->find((*idIt).first));
0084       else if ((*idIt).first.subdetId() == EcalEndcap)
0085         rh = &*(theEEHits_->find((*idIt).first));
0086       if (!rh)
0087         std::cout << "CalibElectron::BIG ERROR::RecHit NOT FOUND" << std::endl;
0088
0089       w_ring += rh->energy();
0090     }
0091
0092     if (w_ring != 0.)
0093       theWeights.push_back(std::pair<int, float>(0, w_ring));
0094     std::cout << " ABS SCALE  - energy sum " << w_ring << std::endl;
0095
0096   }
0097
0098   else if (calibtype == "ETA_ET_MODE") {
0099     float w_ring[200];
0100
0101     for (int i = 0; i < EcalIndexingTools::getInstance()->getNumberOfChannels(); ++i)
0102       w_ring[i] = 0.;
0103
0104     std::vector<std::pair<DetId, float> > scDetIds = theElectron_->superCluster()->hitsAndFractions();
0105
0106     for (std::vector<std::pair<DetId, float> >::const_iterator idIt = scDetIds.begin(); idIt != scDetIds.end();
0107          ++idIt) {
0108       const EcalRecHit* rh = nullptr;
0109       if ((*idIt).first.subdetId() == EcalBarrel)
0110         rh = &*(theHits_->find((*idIt).first));
0111       else if ((*idIt).first.subdetId() == EcalEndcap)
0112         rh = &*(theEEHits_->find((*idIt).first));
0113       if (!rh)
0114         std::cout << "CalibElectron::BIG ERROR::RecHit NOT FOUND" << std::endl;
0115
0116       float eta = fabs(theElectron_->eta());
0117       float theta = 2. * atan(exp(-eta));
0118       float et = theElectron_->superCluster()->energy() * sin(theta);
0119
0120       int in = EcalIndexingTools::getInstance()->getProgressiveIndex(eta, et);
0121
0122       w_ring[in] += rh->energy();
0123       //w_ring[in]+=theElectron_->superCluster()->energy();
0124
0125       std::cout << "CalibElectron::filling channel " << in << " with value " << theElectron_->superCluster()->energy()
0126                 << std::endl;
0127     }
0128
0129     for (int i = 0; i < EcalIndexingTools::getInstance()->getNumberOfChannels(); ++i) {
0130       if (w_ring[i] != 0.) {
0131         theWeights.push_back(std::pair<int, float>(i, w_ring[i]));
0132         std::cout << " ring " << i << " - energy sum " << w_ring[i] << std::endl;
0133       }
0134     }
0135
0136   }
0137
0138   else {
0139     std::cout << "CalibType not yet implemented" << std::endl;
0140   }
0141
0142   return theWeights;
0143 }