test/stubs/EcalShowerContainmentAnalyzer.cc

0001 /**
0002  * \file EcalShowerContainmentAnalyzer
0003  *
0004  * Analyzer to test Shower Containment Corrections
0005  *
0006  * \author S. Argiro'
0007  *
0008 */
0009
0010 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013
0014 #include "CondFormats/EcalObjects/interface/EcalADCToGeVConstant.h"
0015 #include "CondFormats/DataRecord/interface/EcalADCToGeVConstantRcd.h"
0016 #include "CondFormats/EcalCorrections/interface/EcalShowerContainmentCorrections.h"
0017 #include "CondFormats/DataRecord/interface/EcalShowerContainmentCorrectionsRcd.h"
0018 #include "FWCore/ParameterSet/interface/FileInPath.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020
0021 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0022 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0023 #include "TBDataFormats/EcalTBObjects/interface/EcalTBEventHeader.h"
0024 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0025 #include "TBDataFormats/EcalTBObjects/interface/EcalTBHodoscopeRecInfo.h"
0026 #include "TBDataFormats/EcalTBObjects/interface/EcalTBTDCRecInfo.h"
0027 #include "RecoTBCalo/EcalTBAnalysisCoreTools/interface/TBPositionCalc.h"
0028
0029 // geometry
0030 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0031 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0032 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0033 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0034
0035 #include "TFile.h"
0036 #include "TTree.h"
0037
0038 #include <vector>
0039 #include <map>
0040
0041 class EcalShowerContainmentAnalyzer : public edm::one::EDAnalyzer<> {
0042 public:
0043   explicit EcalShowerContainmentAnalyzer(const edm::ParameterSet& ps);
0044   ~EcalShowerContainmentAnalyzer() override;
0045
0046   void analyze(edm::Event const& iEvent, const edm::EventSetup& iSetup) override;
0047   void endJob() override;
0048
0049 protected:
0050   void readIntercalibrationConstants();
0051
0052   std::vector<EBDetId> Xtals3x3(const edm::Event& iEvent, EBDetId centerXtal);
0053   std::vector<EBDetId> Xtals5x5(const edm::Event& iEvent, EBDetId centerXtal);
0054
0055   double energy3x3(const edm::Event& iEvent, EBDetId centerXtal);
0056   double energy5x5(const edm::Event& iEvent, EBDetId centerXtal);
0057   std::pair<double, double> contCorrection(const edm::Event& iEvent,
0058                                            const edm::EventSetup& iSetup,
0059                                            const EBDetId& centerXtal);
0060
0061   /// fill a map of (detid,amplitude)
0062   void fillxtalmap(edm::Event const& iEvent);
0063
0064   double getAdcToGevConstant(const edm::EventSetup& eSetup);
0065
0066   TFile* file_;
0067   TTree* tree_;
0068
0069   // map of < xtal, (coefficient, goodness)>
0070   typedef std::map<EBDetId, std::pair<double, double> > CalibMap;
0071   CalibMap calibmap_;
0072
0073   // map of <xtal, amplitude>
0074   std::map<EBDetId, double> xtalmap_;
0075
0076   const edm::InputTag inputTag_;
0077   const std::string intercalibrationCoeffFile_;
0078
0079   const std::string recHitProducer_;
0080   const std::string ebRecHitCollection_;
0081
0082   const edm::EDGetTokenT<EcalTBEventHeader> tbEvtHeaderToken_;
0083   const edm::EDGetTokenT<EcalTBHodoscopeRecInfo> tbHodoRecInfoToken_;
0084   const edm::EDGetTokenT<EcalUncalibratedRecHitCollection> uncalibRecHitsToken_;
0085   const edm::EDGetTokenT<EBRecHitCollection> ebRecHitsToken_;
0086
0087   const edm::ESGetToken<EcalShowerContainmentCorrections, EcalShowerContainmentCorrectionsRcd> showerContCorrToken_;
0088   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geometryToken_;
0089   const edm::ESGetToken<EcalADCToGeVConstant, EcalADCToGeVConstantRcd> adcToGeVConstantToken_;
0090
0091   struct Ntuple {
0092     int run;
0093     int event;
0094     double energy1x1;
0095     double energy3x3;
0096     double energy5x5;
0097     double corrected3x3;
0098     double corrected5x5;
0099     double posx;
0100     double posy;
0101     double hodo_posx;
0102     double hodo_posy;
0103
0104     void zero() {
0105       run = event = 0;
0106       energy1x1 = energy3x3 = energy5x5 = corrected3x3 = corrected5x5 = posx = posy = hodo_posx = hodo_posy = 0;
0107     }
0108
0109   } ntuple_;
0110 };
0111
0112 DEFINE_FWK_MODULE(EcalShowerContainmentAnalyzer);
0113
0114 EcalShowerContainmentAnalyzer::EcalShowerContainmentAnalyzer(const edm::ParameterSet& ps)
0115     : inputTag_(ps.getUntrackedParameter<edm::InputTag>("EcalUnCalibratedRecHitsLabel")),
0116       intercalibrationCoeffFile_(ps.getUntrackedParameter<std::string>("IntercalibFile")),
0117       recHitProducer_("ecal2006TBRecHit"),
0118       ebRecHitCollection_("EcalRecHitsEB"),
0119       tbEvtHeaderToken_(consumes<EcalTBEventHeader>(edm::InputTag("ecalTBunpack"))),
0120       tbHodoRecInfoToken_(consumes<EcalTBHodoscopeRecInfo>(
0121           edm::InputTag("ecal2006TBHodoscopeReconstructor", "EcalTBHodoscopeRecInfo"))),
0122       uncalibRecHitsToken_(consumes<EcalUncalibratedRecHitCollection>(inputTag_)),
0123       ebRecHitsToken_(consumes<EBRecHitCollection>(edm::InputTag(recHitProducer_, ebRecHitCollection_))),
0124       geometryToken_(esConsumes()),
0125       adcToGeVConstantToken_(esConsumes()) {
0126   std::string outfilename = ps.getUntrackedParameter<std::string>("OutputFile", "testContCorrect.root");
0127
0128   file_ = new TFile(outfilename.c_str(), "recreate");
0129   tree_ = new TTree("test", "test");
0130
0131   tree_->Branch("tb",
0132                 &ntuple_,
0133                 "run/I:"
0134                 "event/I:"
0135                 "energy1x1/D:"
0136                 "energy3x3/D:"
0137                 "energy5x5/D:"
0138                 "corrected3x3/D:"
0139                 "corrected5x5/D:"
0140                 "posx/D:"
0141                 "posy/D:"
0142                 "hodo_posx/D:"
0143                 "hodo_posy/D");
0144
0145   readIntercalibrationConstants();
0146 }
0147
0148 EcalShowerContainmentAnalyzer::~EcalShowerContainmentAnalyzer() {
0149   delete file_;
0150   delete tree_;
0151 }
0152
0153 void EcalShowerContainmentAnalyzer::endJob() {
0154   file_->cd();
0155   tree_->Write();
0156 }
0157
0158 void EcalShowerContainmentAnalyzer::analyze(edm::Event const& iEvent, const edm::EventSetup& iSetup) {
0159   using namespace edm;
0160   using namespace std;
0161
0162   fillxtalmap(iEvent);
0163
0164   Handle<EcalTBEventHeader> pHeader;
0165   iEvent.getByToken(tbEvtHeaderToken_, pHeader);
0166
0167   ntuple_.run = pHeader->runNumber();
0168   ntuple_.event = pHeader->eventNumber();
0169
0170   Handle<EcalTBHodoscopeRecInfo> pHodoscope;
0171   iEvent.getByToken(tbHodoRecInfoToken_, pHodoscope);
0172
0173   ntuple_.hodo_posx = pHodoscope->posX();
0174   ntuple_.hodo_posy = pHodoscope->posY();
0175
0176   Handle<EcalUncalibratedRecHitCollection> pIn;
0177   iEvent.getByToken(uncalibRecHitsToken_, pIn);
0178
0179   // pick the hottest xtal
0180
0181   EcalUncalibratedRecHitCollection::const_iterator iter;
0182   double maxampl = 0;
0183   EBDetId maxid;
0184
0185   for (iter = pIn->begin(); iter != pIn->end(); ++iter) {
0186     EBDetId detid = iter->id();
0187
0188     if (iter->amplitude() > maxampl) {
0189       maxampl = iter->amplitude() * calibmap_[detid].first;
0190       maxid = iter->id();
0191     }
0192   }  // for
0193
0194   double adctogev = getAdcToGevConstant(iSetup);
0195
0196   ntuple_.energy1x1 = maxampl * calibmap_[maxid].first * adctogev;
0197   ntuple_.energy3x3 = energy3x3(iEvent, maxid) * adctogev;
0198   ntuple_.energy5x5 = energy5x5(iEvent, maxid) * adctogev;
0199
0200   std::pair<double, double> correction = contCorrection(iEvent, iSetup, maxid);
0201   ntuple_.corrected3x3 = ntuple_.energy3x3 / correction.first;
0202   ntuple_.corrected5x5 = ntuple_.energy5x5 / correction.second;
0203
0204   tree_->Fill();
0205
0206   ntuple_.zero();
0207 }
0208
0209 /// return vector of detids of the xtals around centerXtal (3x3)
0210 std::vector<EBDetId> EcalShowerContainmentAnalyzer::Xtals3x3(const edm::Event& iEvent, EBDetId centerXtal) {
0211   using namespace edm;
0212   using namespace std;
0213
0214   Handle<EcalUncalibratedRecHitCollection> pIn;
0215   iEvent.getByToken(uncalibRecHitsToken_, pIn);
0216
0217   // find the ids of the 3x3 matrix
0218   vector<EBDetId> Xtals3x3;
0219
0220   for (unsigned int icry = 0; icry < 9; icry++) {
0221     unsigned int row = icry / 3;
0222     unsigned int column = icry % 3;
0223
0224     try {
0225       Xtals3x3.push_back(EBDetId(centerXtal.ieta() + column - 1, centerXtal.iphi() + row - 1, EBDetId::ETAPHIMODE));
0226     } catch (cms::Exception& e) {
0227       Xtals3x3.clear();
0228       return Xtals3x3;
0229     }  //catch
0230   }  // for
0231   return Xtals3x3;
0232 }
0233
0234 /// return vector of detids of the xtals around centerXtal (5x5)
0235 std::vector<EBDetId> EcalShowerContainmentAnalyzer::Xtals5x5(const edm::Event& iEvent, EBDetId centerXtal) {
0236   using namespace edm;
0237   using namespace std;
0238
0239   Handle<EcalUncalibratedRecHitCollection> pIn;
0240   iEvent.getByToken(uncalibRecHitsToken_, pIn);
0241
0242   // find the ids of the 3x3 matrix
0243   vector<EBDetId> Xtals5x5;
0244
0245   for (unsigned int icry = 0; icry < 25; icry++) {
0246     unsigned int row = icry / 5;
0247     unsigned int column = icry % 5;
0248
0249     try {
0250       Xtals5x5.push_back(EBDetId(centerXtal.ieta() + column - 2, centerXtal.iphi() + row - 2, EBDetId::ETAPHIMODE));
0251     } catch (cms::Exception& e) {
0252       Xtals5x5.clear();
0253       return Xtals5x5;
0254     }  //catch
0255   }  // for
0256
0257   return Xtals5x5;
0258 }
0259
0260 void EcalShowerContainmentAnalyzer::readIntercalibrationConstants() {
0261   using namespace std;
0262
0263   ifstream calibrationfile(intercalibrationCoeffFile_.c_str());
0264
0265   // skip header lines
0266   int kHeaderLines = 6;
0267   for (int i = 0; i < kHeaderLines; i++) {
0268     char line[64];
0269     calibrationfile.getline(line, 64);
0270   }
0271
0272   int xtal;
0273   double coeff, sigma;
0274   int nev, good;
0275
0276   while (calibrationfile >> xtal >> coeff >> sigma >> nev >> good) {
0277     EBDetId detid(1, xtal, EBDetId::SMCRYSTALMODE);
0278     calibmap_[detid] = pair<double, double>(coeff, good);
0279   }
0280 }
0281
0282 void EcalShowerContainmentAnalyzer::fillxtalmap(const edm::Event& iEvent) {
0283   using namespace edm;
0284
0285   Handle<EcalUncalibratedRecHitCollection> pIn;
0286   iEvent.getByToken(uncalibRecHitsToken_, pIn);
0287
0288   EcalUncalibratedRecHitCollection::const_iterator iter;
0289
0290   for (iter = pIn->begin(); iter != pIn->end(); ++iter) {
0291     EBDetId detid = iter->id();
0292     double amplitude = iter->amplitude();
0293     xtalmap_[detid] = amplitude;
0294   }  // for
0295 }
0296
0297 /** return a pair with correction coeff for 3x3 in first place
0298     and correction 5x5 in second   */
0299 std::pair<double, double> EcalShowerContainmentAnalyzer::contCorrection(const edm::Event& iEvent,
0300                                                                         const edm::EventSetup& iESetup,
0301                                                                         const EBDetId& centerXtal) {
0302   using namespace std;
0303   using namespace edm;
0304
0305   const auto& gapCorr = iESetup.getData(showerContCorrToken_);
0306
0307   Handle<EBRecHitCollection> pEBRecHits;
0308   iEvent.getByToken(ebRecHitsToken_, pEBRecHits);
0309   const EBRecHitCollection* EBRecHits = pEBRecHits.product();
0310
0311   map<string, double> posparam;
0312   posparam["LogWeighted"] = 1.0;
0313   posparam["X0"] = 0.89;
0314   posparam["T0"] = 6.2;
0315   posparam["W0"] = 4.0;
0316
0317   const auto& geometry = iESetup.getData(geometryToken_);
0318   const CaloSubdetectorGeometry* geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0319
0320   edm::FileInPath mapfile("Geometry/EcalTestBeam/data/BarrelSM1CrystalCenterElectron120GeV.dat");
0321   TBPositionCalc pos(posparam, mapfile.fullPath(), geometry_p);
0322
0323   vector<EBDetId> myDets = Xtals3x3(iEvent, centerXtal);
0324
0325   CLHEP::Hep3Vector clusterPos;
0326   try {
0327     clusterPos = pos.CalculateTBPos(myDets, centerXtal.ic(), EBRecHits);
0328   } catch (cms::Exception& e) {
0329     return make_pair(-99., -99.);
0330   }
0331
0332   ntuple_.posx = clusterPos.x() * 10.;
0333   ntuple_.posy = clusterPos.y() * 10.;
0334
0335   math::XYZPoint mathpoint(clusterPos.x(), clusterPos.y(), clusterPos.z());
0336
0337   double correction3x3 = gapCorr.correction3x3(centerXtal, mathpoint);
0338   double correction5x5 = gapCorr.correction5x5(centerXtal, mathpoint);
0339
0340   return std::make_pair(correction3x3, correction5x5);
0341 }
0342
0343 // retrieve adctogev constant for xtalid from the database
0344 double EcalShowerContainmentAnalyzer::getAdcToGevConstant(const edm::EventSetup& eSetup) {
0345   const auto& ical = eSetup.getData(adcToGeVConstantToken_);
0346   return ical.getEBValue();
0347 }
0348
0349 double EcalShowerContainmentAnalyzer::energy3x3(const edm::Event& iEvent, EBDetId centerXtal) {
0350   using namespace edm;
0351   using namespace std;
0352
0353   vector<EBDetId> xtals3x3 = Xtals3x3(iEvent, centerXtal);
0354
0355   // get the energy
0356   double energy3x3 = 0;
0357   vector<EBDetId>::iterator detIt;
0358   for (detIt = xtals3x3.begin(); detIt != xtals3x3.end(); ++detIt) {
0359     energy3x3 += xtalmap_[*detIt] * calibmap_[*detIt].first;
0360   }
0361
0362   return energy3x3;
0363 }
0364
0365 double EcalShowerContainmentAnalyzer::energy5x5(const edm::Event& iEvent, EBDetId centerXtal) {
0366   using namespace edm;
0367   using namespace std;
0368
0369   vector<EBDetId> xtals5x5 = Xtals5x5(iEvent, centerXtal);
0370
0371   // get the energy
0372   double energy5x5 = 0;
0373   vector<EBDetId>::iterator detIt;
0374   for (detIt = xtals5x5.begin(); detIt != xtals5x5.end(); ++detIt) {
0375     energy5x5 += xtalmap_[*detIt] * calibmap_[*detIt].first;
0376   }
0377
0378   return energy5x5;
0379 }