Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:58:40 

0001 #include <memory>
0002 #include <cmath>
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "Calibration/EcalCalibAlgos/interface/IMACalibBlock.h"
0007 #include "Calibration/EcalCalibAlgos/interface/L3CalibBlock.h"
0008 #include "Calibration/EcalCalibAlgos/interface/EcalEleCalibLooper.h"
0009 #include "DataFormats/Common/interface/Handle.h"
0010 #include "DataFormats/EgammaReco/interface/ClusterShape.h"
0011 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0012 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0013 #include "Calibration/Tools/interface/calibXMLwriter.h"
0014 #include "CalibCalorimetry/CaloMiscalibTools/interface/CaloMiscalibTools.h"
0015 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0016 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0017 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0018 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0019 #include "DataFormats/EgammaReco/interface/ClusterShapeFwd.h"
0020 #include "Calibration/EcalCalibAlgos/interface/MatrixFillMap.h"
0021 #include "Calibration/EcalCalibAlgos/interface/ClusterFillMap.h"
0022 #include "TH1.h"
0023 #include "TH2.h"
0024 #include "TFile.h"
0025 
0026 //----------------------------------------------------------------
0027 
0028 //!LP ctor
0029 EcalEleCalibLooper::EcalEleCalibLooper(const edm::ParameterSet& iConfig)
0030     : m_barrelAlCa(iConfig.getParameter<edm::InputTag>("alcaBarrelHitCollection")),
0031       m_endcapAlCa(iConfig.getParameter<edm::InputTag>("alcaEndcapHitCollection")),
0032       m_ElectronLabel(iConfig.getParameter<edm::InputTag>("electronLabel")),
0033       m_recoWindowSidex(iConfig.getParameter<int>("recoWindowSidex")),
0034       m_recoWindowSidey(iConfig.getParameter<int>("recoWindowSidey")),
0035       m_etaWidth(iConfig.getParameter<int>("etaWidth")),
0036       m_phiWidthEB(iConfig.getParameter<int>("phiWidthEB")),
0037       m_etaStart(etaShifter(iConfig.getParameter<int>("etaStart"))),
0038       m_etaEnd(etaShifter(iConfig.getParameter<int>("etaEnd"))),
0039       m_phiStartEB(iConfig.getParameter<int>("phiStartEB")),
0040       m_phiEndEB(iConfig.getParameter<int>("phiEndEB")),
0041       m_radStart(iConfig.getParameter<int>("radStart")),
0042       m_radEnd(iConfig.getParameter<int>("radEnd")),
0043       m_radWidth(iConfig.getParameter<int>("radWidth")),
0044       m_phiStartEE(iConfig.getParameter<int>("phiStartEE")),
0045       m_phiEndEE(iConfig.getParameter<int>("phiEndEE")),
0046       m_phiWidthEE(iConfig.getParameter<int>("phiWidthEE")),
0047       m_maxSelectedNumPerXtal(iConfig.getParameter<int>("maxSelectedNumPerCrystal")),
0048       m_minEnergyPerCrystal(iConfig.getParameter<double>("minEnergyPerCrystal")),
0049       m_maxEnergyPerCrystal(iConfig.getParameter<double>("maxEnergyPerCrystal")),
0050       m_minCoeff(iConfig.getParameter<double>("minCoeff")),
0051       m_maxCoeff(iConfig.getParameter<double>("maxCoeff")),
0052       m_usingBlockSolver(iConfig.getParameter<int>("usingBlockSolver")),
0053       m_loops(iConfig.getParameter<int>("loops")),
0054       m_ebRecHitToken(consumes<EBRecHitCollection>(m_barrelAlCa)),
0055       m_eeRecHitToken(consumes<EERecHitCollection>(m_endcapAlCa)),
0056       m_gsfElectronToken(consumes<reco::GsfElectronCollection>(m_ElectronLabel)),
0057       m_geometryToken(esConsumes()) {
0058   edm::LogInfo("IML") << "[EcalEleCalibLooper][ctor] asserts";
0059   assert(!((m_etaEnd - m_etaStart) % m_etaWidth));
0060 
0061   assert(m_etaStart >= 0 && m_etaStart <= 171);
0062   assert(m_etaEnd >= m_etaStart && m_etaEnd <= 171);
0063   assert((m_radEnd - m_radStart) % m_radWidth == 0);
0064   assert(m_radStart >= 0 && m_radStart <= 50);
0065   assert(m_radEnd >= m_radStart && m_radEnd <= 50);
0066   edm::LogInfo("IML") << "[EcalEleCalibLooper][ctor] entering ";
0067   edm::LogInfo("IML") << "[EcalEleCalibLooper][ctor] region definition";
0068   EBRegionDefinition();
0069   EERegionDefinition();
0070   ///Graphs to ckeck the region definition
0071   TH2F* EBRegion = new TH2F("EBRegion", "EBRegion", 170, 0, 170, 360, 0, 360);
0072   for (int eta = 0; eta < 170; ++eta)
0073     for (int phi = 0; phi < 360; ++phi) {
0074       EBRegion->Fill(eta, phi, m_xtalRegionId[EBDetId::unhashIndex(eta * 360 + phi).rawId()]);
0075     }
0076   TH2F* EERegion = new TH2F("EERegion", "EERegion", 100, 0, 100, 100, 0, 100);
0077   for (int x = 0; x < 100; ++x)
0078     for (int y = 0; y < 100; ++y) {
0079       if (EEDetId::validDetId(x + 1, y + 1, 1))
0080         EERegion->Fill(x, y, m_xtalRegionId[EEDetId(x + 1, y + 1, -1).rawId()]);
0081     }
0082 
0083   TFile out("EBZone.root", "recreate");
0084   EBRegion->Write();
0085   EERegion->Write();
0086   out.Close();
0087   delete EERegion;
0088   delete EBRegion;
0089   ///End of Graphs
0090 
0091   //PG build the calibration algorithms for the regions
0092   //PG ------------------------------------------------
0093 
0094   edm::LogInfo("IML") << "[EcalEleCalibLooper][ctor] Calib Block";
0095   std::string algorithm = iConfig.getParameter<std::string>("algorithm");
0096   int eventWeight = iConfig.getUntrackedParameter<int>("L3EventWeight", 1);
0097 
0098   //PG loop over the regions set
0099   for (int region = 0; region < EBregionsNum() + 2 * EEregionsNum(); ++region) {
0100     if (algorithm == "IMA")
0101       m_EcalCalibBlocks.push_back(new IMACalibBlock(m_regions.at(region)));
0102     else if (algorithm == "L3")
0103       m_EcalCalibBlocks.push_back(new L3CalibBlock(m_regions.at(region), eventWeight));
0104     else {
0105       edm::LogError("building") << algorithm << " is not a valid calibration algorithm";
0106       exit(1);
0107     }
0108   }  //PG loop over the regions set
0109   std::string mapFiller = iConfig.getParameter<std::string>("FillType");
0110   if (mapFiller == "Cluster")
0111     m_MapFiller = new ClusterFillMap(m_recoWindowSidex,
0112                                      m_recoWindowSidey,
0113                                      m_xtalRegionId,
0114                                      m_minEnergyPerCrystal,
0115                                      m_maxEnergyPerCrystal,
0116                                      m_xtalPositionInRegion,
0117                                      &m_barrelMap,
0118                                      &m_endcapMap);
0119   if (mapFiller == "Matrix")
0120     m_MapFiller = new MatrixFillMap(m_recoWindowSidex,
0121                                     m_recoWindowSidey,
0122                                     m_xtalRegionId,
0123                                     m_minEnergyPerCrystal,
0124                                     m_maxEnergyPerCrystal,
0125                                     m_xtalPositionInRegion,
0126                                     &m_barrelMap,
0127                                     &m_endcapMap);
0128 }  //end ctor
0129 
0130 //---------------------------------------------------------------------------
0131 
0132 //!LP destructor
0133 EcalEleCalibLooper::~EcalEleCalibLooper() {
0134   edm::LogInfo("IML") << "[EcalEleCalibLooper][dtor]";
0135   for (std::vector<VEcalCalibBlock*>::iterator calibBlock = m_EcalCalibBlocks.begin();
0136        calibBlock != m_EcalCalibBlocks.end();
0137        ++calibBlock)
0138     delete (*calibBlock);
0139 }
0140 
0141 //---------------------------------------------------------------------------
0142 
0143 //!BeginOfJob
0144 void EcalEleCalibLooper::beginOfJob() { isfirstcall_ = true; }
0145 
0146 //----------------------------------------------------------------------
0147 
0148 //!startingNewLoop
0149 //!empties the map of the calibBlock so that it can be filled
0150 void EcalEleCalibLooper::startingNewLoop(unsigned int ciclo) {
0151   edm::LogInfo("IML") << "[InvMatrixCalibLooper][Start] entering loop " << ciclo;
0152 
0153   for (std::vector<VEcalCalibBlock*>::iterator calibBlock = m_EcalCalibBlocks.begin();
0154        calibBlock != m_EcalCalibBlocks.end();
0155        ++calibBlock)
0156     (*calibBlock)->reset();
0157   for (std::map<int, int>::iterator it = m_xtalNumOfHits.begin(); it != m_xtalNumOfHits.end(); ++it)
0158     it->second = 0;
0159   return;
0160 }
0161 
0162 //----------------------------------------------------------------------
0163 
0164 //!duringLoop
0165 //!return the status Kcontinue, fills the calibBlock with the recHits
0166 edm::EDLooper::Status EcalEleCalibLooper::duringLoop(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0167   // this chunk used to belong to beginJob(isetup). Moved here
0168   // with the beginJob without arguments migration
0169 
0170   if (isfirstcall_) {
0171     const auto& geometry = iSetup.getData(m_geometryToken);
0172     m_barrelCells = geometry.getValidDetIds(DetId::Ecal, EcalBarrel);
0173     m_endcapCells = geometry.getValidDetIds(DetId::Ecal, EcalEndcap);
0174     for (std::vector<DetId>::const_iterator barrelIt = m_barrelCells.begin(); barrelIt != m_barrelCells.end();
0175          ++barrelIt) {
0176       m_barrelMap[*barrelIt] = 1;
0177       m_xtalNumOfHits[barrelIt->rawId()] = 0;
0178     }
0179     for (std::vector<DetId>::const_iterator endcapIt = m_endcapCells.begin(); endcapIt != m_endcapCells.end();
0180          ++endcapIt) {
0181       m_endcapMap[*endcapIt] = 1;
0182       m_xtalNumOfHits[endcapIt->rawId()] = 0;
0183     }
0184 
0185     isfirstcall_ = false;
0186   }
0187 
0188   //take the collection of recHits in the barrel
0189   const EBRecHitCollection* barrelHitsCollection = nullptr;
0190   edm::Handle<EBRecHitCollection> barrelRecHitsHandle;
0191   iEvent.getByToken(m_ebRecHitToken, barrelRecHitsHandle);
0192   barrelHitsCollection = barrelRecHitsHandle.product();
0193   if (!barrelRecHitsHandle.isValid()) {
0194     edm::LogError("reading") << "[EcalEleCalibLooper] barrel rec hits not found";
0195     return kContinue;
0196   }
0197 
0198   //take the collection of rechis in the endcap
0199   const EERecHitCollection* endcapHitsCollection = nullptr;
0200   edm::Handle<EERecHitCollection> endcapRecHitsHandle;
0201   iEvent.getByToken(m_eeRecHitToken, endcapRecHitsHandle);
0202   endcapHitsCollection = endcapRecHitsHandle.product();
0203   if (!endcapRecHitsHandle.isValid()) {
0204     edm::LogError("reading") << "[EcalEleCalibLooper] endcap rec hits not found";
0205     return kContinue;
0206   }
0207 
0208   //Takes the electron collection of the pixel detector
0209   edm::Handle<reco::GsfElectronCollection> pElectrons;
0210   iEvent.getByToken(m_gsfElectronToken, pElectrons);
0211   if (!pElectrons.isValid()) {
0212     edm::LogError("reading") << "[EcalEleCalibLooper] electrons not found";
0213     return kContinue;
0214   }
0215 
0216   //Start the loop over the electrons
0217   for (reco::GsfElectronCollection::const_iterator eleIt = pElectrons->begin(); eleIt != pElectrons->end(); ++eleIt) {
0218     double pSubtract = 0;
0219     double pTk = 0;
0220     std::map<int, double> xtlMap;
0221     DetId Max = 0;
0222     if (std::abs(eleIt->eta()) < 1.49)
0223       Max = EcalClusterTools::getMaximum(eleIt->superCluster()->hitsAndFractions(), barrelHitsCollection).first;
0224     else
0225       Max = EcalClusterTools::getMaximum(eleIt->superCluster()->hitsAndFractions(), endcapHitsCollection).first;
0226     if (Max.det() == 0)
0227       continue;
0228     m_MapFiller->fillMap(
0229         eleIt->superCluster()->hitsAndFractions(), Max, barrelHitsCollection, endcapHitsCollection, xtlMap, pSubtract);
0230     if (m_maxSelectedNumPerXtal > 0 && m_xtalNumOfHits[Max.rawId()] > m_maxSelectedNumPerXtal)
0231       continue;
0232     ++m_xtalNumOfHits[Max.rawId()];
0233     if (m_xtalRegionId[Max.rawId()] == -1)
0234       continue;
0235     pTk = eleIt->trackMomentumAtVtx().R();
0236     m_EcalCalibBlocks.at(m_xtalRegionId[Max.rawId()])->Fill(xtlMap.begin(), xtlMap.end(), pTk, pSubtract);
0237   }  //End of the loop over the electron collection
0238 
0239   return kContinue;
0240 }  //end of duringLoop
0241 
0242 //----------------------------------------------------------------------
0243 
0244 //! EndOfLoop
0245 //!Return kContinue if there's still another loop to be done; otherwise stops returnig kStop;
0246 //!Takes the coefficients solving the calibBlock;
0247 edm::EDLooper::Status EcalEleCalibLooper::endOfLoop(const edm::EventSetup& dumb, unsigned int iCounter) {
0248   edm::LogInfo("IML") << "[InvMatrixCalibLooper][endOfLoop] entering...";
0249   for (std::vector<VEcalCalibBlock*>::iterator calibBlock = m_EcalCalibBlocks.begin();
0250        calibBlock != m_EcalCalibBlocks.end();
0251        ++calibBlock)
0252     (*calibBlock)->solve(m_usingBlockSolver, m_minCoeff, m_maxCoeff);
0253 
0254   TH1F* EBcoeffEnd = new TH1F("EBRegion", "EBRegion", 100, 0.5, 2.1);
0255   TH2F* EBcoeffMap = new TH2F("EBcoeff", "EBcoeff", 171, -85, 85, 360, 1, 361);
0256   TH1F* EEPcoeffEnd = new TH1F("EEPRegion", "EEPRegion", 100, 0.5, 2.1);
0257   TH1F* EEMcoeffEnd = new TH1F("EEMRegion", "EEMRegion", 100, 0.5, 2.1);
0258   TH2F* EEPcoeffMap = new TH2F("EEPcoeffMap", "EEPcoeffMap", 101, 1, 101, 101, 0, 101);
0259   TH2F* EEMcoeffMap = new TH2F("EEMcoeffMap", "EEMcoeffMap", 101, 1, 101, 101, 0, 101);
0260   //loop over the barrel xtals to get the coeffs
0261   for (std::vector<DetId>::const_iterator barrelIt = m_barrelCells.begin(); barrelIt != m_barrelCells.end();
0262        ++barrelIt) {
0263     EBDetId ee(*barrelIt);
0264     int index = barrelIt->rawId();
0265     if (m_xtalRegionId[index] == -1)
0266       continue;
0267     m_barrelMap[*barrelIt] *= m_EcalCalibBlocks.at(m_xtalRegionId[index])->at(m_xtalPositionInRegion[index]);
0268     EBcoeffEnd->Fill(m_barrelMap[*barrelIt]);
0269     EBcoeffMap->Fill(ee.ieta(), ee.iphi(), m_barrelMap[*barrelIt]);
0270   }  //PG loop over phi
0271 
0272   // loop over the EndCap to get the recalib coefficients
0273   for (std::vector<DetId>::const_iterator endcapIt = m_endcapCells.begin(); endcapIt != m_endcapCells.end();
0274        ++endcapIt) {
0275     EEDetId ee(*endcapIt);
0276     int index = endcapIt->rawId();
0277     if (ee.zside() > 0) {
0278       if (m_xtalRegionId[index] == -1)
0279         continue;
0280       m_endcapMap[*endcapIt] *= m_EcalCalibBlocks.at(m_xtalRegionId[index])->at(m_xtalPositionInRegion[index]);
0281       EEPcoeffEnd->Fill(m_endcapMap[*endcapIt]);
0282       EEPcoeffMap->Fill(ee.ix(), ee.iy(), m_endcapMap[*endcapIt]);
0283     } else {
0284       m_endcapMap[*endcapIt] *= m_EcalCalibBlocks.at(m_xtalRegionId[index])->at(m_xtalPositionInRegion[index]);
0285       EEMcoeffEnd->Fill(m_endcapMap[*endcapIt]);
0286       EEMcoeffMap->Fill(ee.ix(), ee.iy(), m_endcapMap[*endcapIt]);
0287     }
0288   }  // loop over the EndCap to get the recalib coefficients
0289 
0290   edm::LogInfo("IML") << "[InvMatrixCalibLooper][endOfLoop] End of endOfLoop";
0291 
0292   char filename[80];
0293   sprintf(filename, "coeffs%d.root", iCounter);
0294   TFile zout(filename, "recreate");
0295   EBcoeffEnd->Write();
0296   EBcoeffMap->Write();
0297   EEPcoeffEnd->Write();
0298   EEPcoeffMap->Write();
0299   EEMcoeffEnd->Write();
0300   EEMcoeffMap->Write();
0301   zout.Close();
0302   delete EBcoeffEnd;
0303   delete EBcoeffMap;
0304   delete EEPcoeffEnd;
0305   delete EEMcoeffEnd;
0306   delete EEPcoeffMap;
0307   delete EEMcoeffMap;
0308   if (iCounter < m_loops - 1)
0309     return kContinue;
0310   else
0311     return kStop;
0312 }
0313 
0314 //-------------------------------------------------------------------
0315 
0316 //!LP endOfJob
0317 //!writes the coefficients in the xml format and exits
0318 void EcalEleCalibLooper::endOfJob() {
0319   edm::LogInfo("IML") << "[InvMatrixCalibLooper][endOfJob] saving calib coeffs";
0320 
0321   //Writes the coeffs
0322   calibXMLwriter barrelWriter(EcalBarrel);
0323   calibXMLwriter endcapWriter(EcalEndcap);
0324   for (std::vector<DetId>::const_iterator barrelIt = m_barrelCells.begin(); barrelIt != m_barrelCells.end();
0325        ++barrelIt) {
0326     EBDetId eb(*barrelIt);
0327     barrelWriter.writeLine(eb, m_barrelMap[*barrelIt]);
0328   }
0329   for (std::vector<DetId>::const_iterator endcapIt = m_endcapCells.begin(); endcapIt != m_endcapCells.end();
0330        ++endcapIt) {
0331     EEDetId ee(*endcapIt);
0332     endcapWriter.writeLine(ee, m_endcapMap[*endcapIt]);
0333   }
0334   edm::LogInfo("IML") << "[InvMatrixCalibLooper][endOfJob] Exiting";
0335 }
0336 
0337 //------------------------------------//
0338 //      definition of functions       //
0339 //------------------------------------//
0340 
0341 //! Tells if you are in the region to be calibrated
0342 int EcalEleCalibLooper::EBregionCheck(const int eta, const int phi) const {
0343   if (eta < m_etaStart)
0344     return 1;
0345   if (eta >= m_etaEnd)
0346     return 2;
0347   if (phi < m_phiStartEB)
0348     return 3;
0349   if (phi >= m_phiEndEB)
0350     return 4;
0351   return 0;
0352 }
0353 
0354 //--------------------------------------------
0355 
0356 //! def degrees
0357 inline double degrees(double radiants) { return radiants * 180 * (1 / M_PI); }
0358 
0359 //--------------------------------------------
0360 
0361 //!DS def radiants
0362 inline double radiants(int degrees) { return degrees * M_PI * (1 / 180); }
0363 
0364 //--------------------------------------------
0365 
0366 //! copes with the infinitives of the tangent
0367 double EcalEleCalibLooper::giveLimit(int degrees) {
0368   //PG 200 > atan (50/0.5)
0369   if (degrees == 90)
0370     return 90;
0371   return tan(radiants(degrees));
0372 }
0373 
0374 //--------------------------------------------
0375 
0376 //! autoexplaining
0377 inline double Mod(double phi) {
0378   if (phi >= 360 && phi < 720)
0379     return phi - 360;
0380   if (phi >= 720)
0381     return phi - 720;
0382   return phi;
0383 }
0384 
0385 //----------------------------------------
0386 
0387 //!Reg Id generator EB ----- for the barrel
0388 int EcalEleCalibLooper::EBRegionId(const int etaXtl, const int phiXtl) const {
0389   if (EBregionCheck(etaXtl, phiXtl))
0390     return -1;
0391   int phifake = m_phiStartEB;
0392   if (m_phiStartEB > m_phiEndEB)
0393     phifake = m_phiStartEB - 360;
0394   int Nphi = (m_phiEndEB - phifake) / m_phiWidthEB;
0395   int etaI = (etaXtl - m_etaStart) / m_etaWidth;
0396   int phiI = (phiXtl - m_phiStartEB) / m_phiWidthEB;
0397   int regionNumEB = phiI + Nphi * etaI;
0398   return (int)regionNumEB;
0399 }
0400 
0401 //----------------------------------------
0402 
0403 //! Gives the id of the region
0404 int EcalEleCalibLooper::EERegionId(const int ics, const int ips) const {
0405   if (EEregionCheck(ics, ips))
0406     return -1;
0407   int phifake = m_phiStartEE;
0408   if (m_phiStartEE > m_phiEndEE)
0409     phifake = m_phiStartEE - 360;
0410   double radius = (ics - 50) * (ics - 50) + (ips - 50) * (ips - 50);
0411   radius = sqrt(radius);
0412   int Nphi = (m_phiEndEE - phifake) / m_phiWidthEE;
0413   double phi = atan2(static_cast<double>(ips - 50), static_cast<double>(ics - 50));
0414   phi = degrees(phi);
0415   if (phi < 0)
0416     phi += 360;
0417   int radI = static_cast<int>((radius - m_radStart) / m_radWidth);
0418   int phiI = static_cast<int>((m_phiEndEE - phi) / m_phiWidthEE);
0419   int regionNumEE = phiI + Nphi * radI;
0420   return regionNumEE;
0421 }
0422 
0423 //----------------------------------------
0424 
0425 //!DS Number of regions in EE
0426 inline int EcalEleCalibLooper::EEregionsNum() const {
0427   int phifake = m_phiStartEE;
0428   if (m_phiStartEE > m_phiEndEE)
0429     phifake = m_phiStartEE - 360;
0430   return ((m_radEnd - m_radStart) / m_radWidth) * ((m_phiEndEE - phifake) / m_phiWidthEE);
0431 }
0432 
0433 //----------------------------------------
0434 
0435 //!DS number of regions in EB
0436 inline int EcalEleCalibLooper::EBregionsNum() const {
0437   int phi = m_phiStartEB;
0438   if (m_phiStartEB > m_phiEndEB)
0439     phi = m_phiStartEB - 360;
0440   return ((m_etaEnd - m_etaStart) / m_etaWidth) * ((m_phiEndEB - phi) / m_phiWidthEB);
0441 }
0442 
0443 //----------------------------------------
0444 
0445 //!DS EB Region Definition
0446 void EcalEleCalibLooper::EBRegionDefinition() {
0447   int reg = -1;
0448   for (int it = 0; it < EBregionsNum(); ++it)
0449     m_regions.push_back(0);
0450   for (int eta = 0; eta < 170; ++eta)
0451     for (int phi = 0; phi < 360; ++phi) {
0452       reg = EBRegionId(eta, phi);
0453       m_xtalRegionId[EBDetId::unhashIndex(eta * 360 + phi).rawId()] = reg;
0454       if (reg == -1)
0455         continue;
0456       m_xtalPositionInRegion[EBDetId::unhashIndex(eta * 360 + phi).rawId()] = m_regions.at(reg);
0457       ++m_regions.at(reg);
0458     }
0459 }
0460 
0461 //----------------------------------------
0462 
0463 //DS EE Region Definition
0464 void EcalEleCalibLooper::EERegionDefinition() {
0465   // reset
0466   int EBnum = EBregionsNum();
0467   int EEnum = EEregionsNum();
0468   for (int it = 0; it < 2 * EEnum; ++it)
0469     m_regions.push_back(0);
0470   // loop sui xtl
0471   int reg = -1;
0472   for (int ics = 0; ics < 100; ++ics)
0473     for (int ips = 0; ips < 100; ++ips) {
0474       int ireg = EERegionId(ics, ips);
0475       if (ireg == -1)
0476         reg = -1;
0477       else
0478         reg = EBnum + ireg;
0479       if (EEDetId::validDetId(ics + 1, ips + 1, 1)) {
0480         m_xtalRegionId[EEDetId(ics + 1, ips + 1, 1).rawId()] = reg;
0481         if (reg == -1)
0482           continue;
0483         m_xtalPositionInRegion[EEDetId(ics + 1, ips + 1, 1).rawId()] = m_regions.at(reg);
0484         ++m_regions.at(reg);
0485       }
0486       if (reg != -1)
0487         reg += EEnum;
0488       if (EEDetId::validDetId(ics + 1, ips + 1, -1)) {
0489         m_xtalRegionId[EEDetId(ics + 1, ips + 1, -1).rawId()] = reg;
0490         if (reg == -1)
0491           continue;
0492         m_xtalPositionInRegion[EEDetId(ics + 1, ips + 1, -1).rawId()] = m_regions.at(reg);
0493         ++m_regions.at(reg);
0494       }
0495     }
0496 }
0497 
0498 //-----------------------------------------
0499 
0500 //!returns zero if the coordinates are in the right place.
0501 int EcalEleCalibLooper::EEregionCheck(const int ics, const int ips) const {
0502   int x = ics - 50;
0503   int y = ips - 50;
0504   double radius2 = x * x + y * y;
0505   if (radius2 < 10 * 10)
0506     return 1;  //center of the donut
0507   if (radius2 > 50 * 50)
0508     return 1;  //outer part of the donut
0509   if (radius2 < m_radStart * m_radStart)
0510     return 2;
0511   if (radius2 >= m_radEnd * m_radEnd)
0512     return 2;
0513   double phi = atan2(static_cast<double>(y), static_cast<double>(x));
0514   phi = degrees(phi);
0515   if (phi < 0)
0516     phi += 360;
0517   if (m_phiStartEE < m_phiEndEE && phi > m_phiStartEE && phi < m_phiEndEE)
0518     return 0;
0519   if (m_phiStartEE > m_phiEndEE && (phi > m_phiStartEE || phi < m_phiEndEE))
0520     return 0;
0521   return 3;
0522 }
0523 
0524 //--------------------------------------------
0525 
0526 //Shifts eta in other coordinates (from 0 to 170)
0527 int EcalEleCalibLooper::etaShifter(const int etaOld) const {
0528   if (etaOld < 0)
0529     return etaOld + 85;
0530   else if (etaOld > 0)
0531     return etaOld + 84;
0532   assert(0 != etaOld);  // etaOld = 0, apparently not a foreseen value, so fail
0533   return 999;           // dummy statement to silence compiler warning
0534 }
0535 
0536 //--------------------------------------------
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!