Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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/InvRingCalib.h"
0007 #include "Calibration/EcalCalibAlgos/interface/IMACalibBlock.h"
0008 #include "Calibration/EcalCalibAlgos/interface/L3CalibBlock.h"
0009 #include "DataFormats/Common/interface/Handle.h"
0010 #include "Calibration/Tools/interface/calibXMLwriter.h"
0011 #include "CalibCalorimetry/CaloMiscalibTools/interface/MiscalibReaderFromXMLEcalBarrel.h"
0012 #include "CalibCalorimetry/CaloMiscalibTools/interface/MiscalibReaderFromXMLEcalEndcap.h"
0013 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0014 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0015 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0016 
0017 #include "Calibration/EcalCalibAlgos/interface/MatrixFillMap.h"
0018 #include "Calibration/EcalCalibAlgos/interface/ClusterFillMap.h"
0019 
0020 //Not to remain in the final version
0021 #include "TH2.h"
0022 #include "TFile.h"
0023 //----------------------------------------------------------------
0024 //ctor
0025 
0026 InvRingCalib::InvRingCalib(const edm::ParameterSet& iConfig)
0027     : m_barrelAlCa(iConfig.getParameter<edm::InputTag>("barrelAlca")),
0028       m_endcapAlCa(iConfig.getParameter<edm::InputTag>("endcapAlca")),
0029       m_ElectronLabel(iConfig.getParameter<edm::InputTag>("ElectronLabel")),
0030       m_recoWindowSidex(iConfig.getParameter<int>("recoWindowSidex")),
0031       m_recoWindowSidey(iConfig.getParameter<int>("recoWindowSidey")),
0032       m_minEnergyPerCrystal(iConfig.getParameter<double>("minEnergyPerCrystal")),
0033       m_maxEnergyPerCrystal(iConfig.getParameter<double>("maxEnergyPerCrystal")),
0034       m_etaStart(iConfig.getParameter<int>("etaStart")),
0035       m_etaEnd(iConfig.getParameter<int>("etaEnd")),
0036       m_etaWidth(iConfig.getParameter<int>("etaWidth")),
0037       m_maxSelectedNumPerRing(iConfig.getParameter<int>("maxNumPerRing")),
0038       m_minCoeff(iConfig.getParameter<double>("minCoeff")),
0039       m_maxCoeff(iConfig.getParameter<double>("maxCoeff")),
0040       m_usingBlockSolver(iConfig.getParameter<int>("usingBlockSolver")),
0041       m_startRing(iConfig.getParameter<int>("startRing")),
0042       m_endRing(iConfig.getParameter<int>("endRing")),
0043       m_EBcoeffFile(iConfig.getParameter<std::string>("EBcoeffs")),
0044       m_EEcoeffFile(iConfig.getParameter<std::string>("EEcoeffs")),
0045       m_EEZone(iConfig.getParameter<int>("EEZone")),
0046       m_ebRecHitToken(consumes<EBRecHitCollection>(m_barrelAlCa)),
0047       m_eeRecHitToken(consumes<EERecHitCollection>(m_endcapAlCa)),
0048       m_gsfElectronToken(consumes<reco::GsfElectronCollection>(m_ElectronLabel)),
0049       m_geometryToken(esConsumes()) {
0050   //controls if the parameters inputed are correct
0051   if ((m_etaEnd * m_etaStart) > 0)
0052     assert(!((m_etaEnd - m_etaStart) % m_etaWidth));
0053   if ((m_etaEnd * m_etaStart) < 0)
0054     assert(!((m_etaEnd - m_etaStart - 1) % m_etaWidth));
0055 
0056   assert(m_etaStart >= -85 && m_etaStart <= 86);
0057   assert(m_etaEnd >= m_etaStart && m_etaEnd <= 86);
0058   assert(m_startRing > -1 && m_startRing <= 40);
0059   assert(m_endRing >= m_startRing && m_endRing <= 40);
0060 
0061   assert(!((m_endRing - m_startRing) % m_etaWidth));
0062   assert((abs(m_EEZone) <= 1));
0063 
0064   m_loops = (unsigned int)iConfig.getParameter<int>("loops") - 1;
0065   //LP CalibBlock vector instantiation
0066   edm::LogInfo("IML") << "[InvRingCalib][ctor] Calib Block";
0067   std::string algorithm = iConfig.getParameter<std::string>("algorithm");
0068   m_mapFillerType = iConfig.getParameter<std::string>("FillType");
0069   int eventWeight = iConfig.getUntrackedParameter<int>("L3EventWeight", 1);
0070 
0071   for (int i = 0; i < EBRegionNum(); ++i) {
0072     if (algorithm == "IMA")
0073       m_IMACalibBlocks.push_back(new IMACalibBlock(m_etaWidth));
0074     else if (algorithm == "L3")
0075       m_IMACalibBlocks.push_back(new L3CalibBlock(m_etaWidth, eventWeight));
0076     else {
0077       edm::LogError("building") << algorithm << " is not a valid calibration algorithm";
0078       exit(1);
0079     }
0080   }
0081   int EEBlocks = 0;
0082   if (m_EEZone == 0)
0083     EEBlocks = 2 * EERegionNum();
0084   if (m_EEZone == 1 || m_EEZone == -1)
0085     EEBlocks = EERegionNum();
0086 
0087   for (int i = 0; i < EEBlocks; ++i) {
0088     if (algorithm == "IMA")
0089       m_IMACalibBlocks.push_back(new IMACalibBlock(m_etaWidth));
0090     else if (algorithm == "L3")
0091       m_IMACalibBlocks.push_back(new L3CalibBlock(m_etaWidth, eventWeight));
0092     else {
0093       edm::LogError("building") << algorithm << " is not a valid calibration algorithm";
0094       exit(1);
0095     }
0096   }
0097   edm::LogInfo("IML") << " [InvRingCalib][ctor] end of creator";
0098 }
0099 
0100 //-------------------------------------------------------------- end ctor
0101 //!destructor
0102 
0103 InvRingCalib::~InvRingCalib() {}
0104 
0105 //---------------------------------------------------
0106 
0107 //!BeginOfJob
0108 void InvRingCalib::beginOfJob() { isfirstcall_ = true; }
0109 
0110 //--------------------------------------------------------
0111 
0112 //!startingNewLoop
0113 void InvRingCalib::startingNewLoop(unsigned int ciclo) {
0114   edm::LogInfo("IML") << "[InvMatrixCalibLooper][Start] entering loop " << ciclo;
0115   for (std::vector<VEcalCalibBlock*>::iterator calibBlock = m_IMACalibBlocks.begin();
0116        calibBlock != m_IMACalibBlocks.end();
0117        ++calibBlock) {
0118     //LP empties the energies vector, to fill DuringLoop.
0119     (*calibBlock)->reset();
0120   }
0121   for (std::map<int, int>::const_iterator ring = m_xtalRing.begin(); ring != m_xtalRing.end(); ++ring)
0122     m_RingNumOfHits[ring->second] = 0;
0123   return;
0124 }
0125 
0126 //--------------------------------------------------------
0127 
0128 //!duringLoop
0129 edm::EDLooper::Status InvRingCalib::duringLoop(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0130   if (isfirstcall_) {
0131     edm::LogInfo("IML") << "[InvRingCalib][beginOfJob]";
0132     //gets the geometry from the event setup
0133     const auto& geometry = iSetup.getData(m_geometryToken);
0134     edm::LogInfo("IML") << "[InvRingCalib] Event Setup read";
0135     //fills a vector with all the cells
0136     m_barrelCells = geometry.getValidDetIds(DetId::Ecal, EcalBarrel);
0137     m_endcapCells = geometry.getValidDetIds(DetId::Ecal, EcalEndcap);
0138     //Defines the EB regions
0139     edm::LogInfo("IML") << "[InvRingCalib] Defining Barrel Regions";
0140     EBRegionDef();
0141     //Defines what is a ring in the EE
0142     edm::LogInfo("IML") << "[InvRingCalib] Defining endcap Rings";
0143     EERingDef(iSetup);
0144     //Defines the regions in the EE
0145     edm::LogInfo("IML") << "[InvRingCalib] Defining endcap Regions";
0146     EERegionDef();
0147     if (m_mapFillerType == "Cluster")
0148       m_MapFiller = new ClusterFillMap(m_recoWindowSidex,
0149                                        m_recoWindowSidey,
0150                                        m_xtalRegionId,
0151                                        m_minEnergyPerCrystal,
0152                                        m_maxEnergyPerCrystal,
0153                                        m_RinginRegion,
0154                                        &m_barrelMap,
0155                                        &m_endcapMap);
0156     if (m_mapFillerType == "Matrix")
0157       m_MapFiller = new MatrixFillMap(m_recoWindowSidex,
0158                                       m_recoWindowSidey,
0159                                       m_xtalRegionId,
0160                                       m_minEnergyPerCrystal,
0161                                       m_maxEnergyPerCrystal,
0162                                       m_RinginRegion,
0163                                       &m_barrelMap,
0164                                       &m_endcapMap);
0165     edm::LogInfo("IML") << "[InvRingCalib] Initializing the coeffs";
0166     //Sets the initial coefficients to 1.
0167     //Graphs to check ring, regions and so on, not needed in the final version
0168     TH2F EBRegion("EBRegion", "EBRegion", 171, -85, 86, 360, 1, 361);
0169     TH2F EBRing("EBRing", "EBRing", 171, -85, 86, 360, 1, 361);
0170     for (std::vector<DetId>::const_iterator it = m_barrelCells.begin(); it != m_barrelCells.end(); ++it) {
0171       EBDetId eb(*it);
0172       EBRing.Fill(eb.ieta(), eb.iphi(), m_RinginRegion[it->rawId()]);
0173       EBRegion.Fill(eb.ieta(), eb.iphi(), m_xtalRegionId[it->rawId()]);
0174     }
0175 
0176     TH2F EEPRegion("EEPRegion", "EEPRegion", 100, 1, 101, 100, 1, 101);
0177     TH2F EEPRing("EEPRing", "EEPRing", 100, 1, 101, 100, 1, 101);
0178     TH2F EEPRingReg("EEPRingReg", "EEPRingReg", 100, 1, 101, 100, 1, 101);
0179     TH2F EEMRegion("EEMRegion", "EEMRegion", 100, 1, 101, 100, 1, 101);
0180     TH2F EEMRing("EEMRing", "EEMRing", 100, 1, 101, 100, 1, 101);
0181     TH2F EEMRingReg("EEMRingReg", "EEMRingReg", 100, 1, 101, 100, 1, 101);
0182     //  TH1F eta ("eta","eta",250,-85,165);
0183     for (std::vector<DetId>::const_iterator it = m_endcapCells.begin(); it != m_endcapCells.end(); ++it) {
0184       EEDetId ee(*it);
0185       if (ee.zside() > 0) {
0186         EEPRegion.Fill(ee.ix(), ee.iy(), m_xtalRegionId[ee.rawId()]);
0187         EEPRing.Fill(ee.ix(), ee.iy(), m_xtalRing[ee.rawId()]);
0188         EEPRingReg.Fill(ee.ix(), ee.iy(), m_RinginRegion[ee.rawId()]);
0189       }
0190       if (ee.zside() < 0) {
0191         EEMRegion.Fill(ee.ix(), ee.iy(), m_xtalRegionId[ee.rawId()]);
0192         EEMRing.Fill(ee.ix(), ee.iy(), m_xtalRing[ee.rawId()]);
0193         EEMRingReg.Fill(ee.ix(), ee.iy(), m_RinginRegion[ee.rawId()]);
0194       }
0195     }
0196 
0197     //  for (std::map<int,float>::iterator it=m_eta.begin();
0198     //        it!=m_eta.end();++it)
0199     //         eta.Fill(it->first,it->second);
0200     TFile out("EBZone.root", "recreate");
0201     EBRegion.Write();
0202     EBRing.Write();
0203     EEPRegion.Write();
0204     EEPRing.Write();
0205     EEPRingReg.Write();
0206     EEMRegion.Write();
0207     EEMRing.Write();
0208     //  eta.Write();
0209     EEMRingReg.Write();
0210     out.Close();
0211     edm::LogInfo("IML") << "[InvRingCalib] Start to acquire the coeffs";
0212     CaloMiscalibMapEcal EBmap;
0213     EBmap.prefillMap();
0214     MiscalibReaderFromXMLEcalBarrel barrelreader(EBmap);
0215     if (!m_EBcoeffFile.empty())
0216       barrelreader.parseXMLMiscalibFile(m_EBcoeffFile);
0217     EcalIntercalibConstants costants(EBmap.get());
0218     m_barrelMap = costants.getMap();
0219     CaloMiscalibMapEcal EEmap;
0220     EEmap.prefillMap();
0221     MiscalibReaderFromXMLEcalEndcap endcapreader(EEmap);
0222     if (!m_EEcoeffFile.empty())
0223       endcapreader.parseXMLMiscalibFile(m_EEcoeffFile);
0224     EcalIntercalibConstants EEcostants(EEmap.get());
0225     m_endcapMap = EEcostants.getMap();
0226 
0227     isfirstcall_ = false;
0228   }  // if isfirstcall
0229 
0230   //gets the barrel recHits
0231   double pSubtract = 0.;
0232   double pTk = 0.;
0233   const EcalRecHitCollection* barrelHitsCollection = nullptr;
0234   edm::Handle<EBRecHitCollection> barrelRecHitsHandle;
0235   iEvent.getByToken(m_ebRecHitToken, barrelRecHitsHandle);
0236   barrelHitsCollection = barrelRecHitsHandle.product();
0237 
0238   if (!barrelRecHitsHandle.isValid()) {
0239     edm::LogError("IML") << "[EcalEleCalibLooper] barrel rec hits not found";
0240     return kContinue;
0241   }
0242 
0243   //gets the endcap recHits
0244   const EcalRecHitCollection* endcapHitsCollection = nullptr;
0245   edm::Handle<EERecHitCollection> endcapRecHitsHandle;
0246   iEvent.getByToken(m_eeRecHitToken, endcapRecHitsHandle);
0247   endcapHitsCollection = endcapRecHitsHandle.product();
0248 
0249   if (!endcapRecHitsHandle.isValid()) {
0250     edm::LogError("IML") << "[EcalEleCalibLooper] endcap rec hits not found";
0251     return kContinue;
0252   }
0253 
0254   //gets the electrons
0255   edm::Handle<reco::GsfElectronCollection> pElectrons;
0256   iEvent.getByToken(m_gsfElectronToken, pElectrons);
0257 
0258   if (!pElectrons.isValid()) {
0259     edm::LogError("IML") << "[EcalEleCalibLooper] electrons not found";
0260     return kContinue;
0261   }
0262 
0263   //loops over the electrons in the event
0264   for (reco::GsfElectronCollection::const_iterator eleIt = pElectrons->begin(); eleIt != pElectrons->end(); ++eleIt) {
0265     pSubtract = 0;
0266     pTk = eleIt->trackMomentumAtVtx().R();
0267     std::map<int, double> xtlMap;
0268     DetId Max = 0;
0269     if (std::abs(eleIt->eta()) < 1.49)
0270       Max = EcalClusterTools::getMaximum(eleIt->superCluster()->hitsAndFractions(), barrelHitsCollection).first;
0271     else
0272       Max = EcalClusterTools::getMaximum(eleIt->superCluster()->hitsAndFractions(), endcapHitsCollection).first;
0273     if (Max.det() == 0)
0274       continue;
0275     m_MapFiller->fillMap(
0276         eleIt->superCluster()->hitsAndFractions(), Max, barrelHitsCollection, endcapHitsCollection, xtlMap, pSubtract);
0277     if (m_xtalRegionId[Max.rawId()] == -1)
0278       continue;
0279     pSubtract += eleIt->superCluster()->preshowerEnergy();
0280     ++m_RingNumOfHits[m_xtalRing[Max.rawId()]];
0281     //fills the calibBlock
0282     m_IMACalibBlocks.at(m_xtalRegionId[Max.rawId()])->Fill(xtlMap.begin(), xtlMap.end(), pTk, pSubtract);
0283   }
0284   return kContinue;
0285 }  //end of duringLoop
0286 
0287 //-------------------------------------
0288 
0289 //EndOfLoop
0290 edm::EDLooper::Status InvRingCalib::endOfLoop(const edm::EventSetup& dumb, unsigned int iCounter) {
0291   std::map<int, double> InterRings;
0292   edm::LogInfo("IML") << "[InvMatrixCalibLooper][endOfLoop] Start to invert the matrixes";
0293   //call the autoexplaining "solve" method for every calibBlock
0294   for (std::vector<VEcalCalibBlock*>::iterator calibBlock = m_IMACalibBlocks.begin();
0295        calibBlock != m_IMACalibBlocks.end();
0296        ++calibBlock)
0297     (*calibBlock)->solve(m_usingBlockSolver, m_minCoeff, m_maxCoeff);
0298 
0299   edm::LogInfo("IML") << "[InvRingLooper][endOfLoop] Starting to write the coeffs";
0300   TH1F* coeffDistr = new TH1F("coeffdistr", "coeffdistr", 100, 0.7, 1.4);
0301   TH1F* coeffMap = new TH1F("coeffRingMap", "coeffRingMap", 250, -85, 165);
0302   TH1F* ringDistr = new TH1F("ringDistr", "ringDistr", 250, -85, 165);
0303   TH1F* RingFill = new TH1F("RingFill", "RingFill", 250, -85, 165);
0304   for (std::map<int, int>::const_iterator it = m_xtalRing.begin(); it != m_xtalRing.end(); ++it)
0305     ringDistr->Fill(it->second + 0.1);
0306 
0307   int ID;
0308   std::map<int, int> flag;
0309   for (std::map<int, int>::const_iterator it = m_xtalRing.begin(); it != m_xtalRing.end(); ++it)
0310     flag[it->second] = 0;
0311 
0312   for (std::vector<DetId>::const_iterator it = m_barrelCells.begin(); it != m_barrelCells.end(); ++it) {
0313     ID = it->rawId();
0314     if (m_xtalRegionId[ID] == -1)
0315       continue;
0316     if (flag[m_xtalRing[ID]])
0317       continue;
0318     flag[m_xtalRing[ID]] = 1;
0319     RingFill->Fill(m_xtalRing[ID], m_RingNumOfHits[m_xtalRing[ID]]);
0320     InterRings[m_xtalRing[ID]] = m_IMACalibBlocks.at(m_xtalRegionId[ID])->at(m_RinginRegion[ID]);
0321     coeffMap->Fill(m_xtalRing[ID] + 0.1, InterRings[m_xtalRing[ID]]);
0322     coeffDistr->Fill(InterRings[m_xtalRing[ID]]);
0323   }
0324 
0325   for (std::vector<DetId>::const_iterator it = m_endcapCells.begin(); it != m_endcapCells.end(); ++it) {
0326     ID = it->rawId();
0327     if (m_xtalRegionId[ID] == -1)
0328       continue;
0329     if (flag[m_xtalRing[ID]])
0330       continue;
0331     flag[m_xtalRing[ID]] = 1;
0332     InterRings[m_xtalRing[ID]] = m_IMACalibBlocks.at(m_xtalRegionId[ID])->at(m_RinginRegion[ID]);
0333     RingFill->Fill(m_xtalRing[ID], m_RingNumOfHits[m_xtalRing[ID]]);
0334     coeffMap->Fill(m_xtalRing[ID], InterRings[m_xtalRing[ID]]);
0335     coeffDistr->Fill(InterRings[m_xtalRing[ID]]);
0336   }
0337 
0338   char filename[80];
0339   sprintf(filename, "coeff%d.root", iCounter);
0340   TFile out(filename, "recreate");
0341   coeffDistr->Write();
0342   coeffMap->Write();
0343   ringDistr->Write();
0344   RingFill->Write();
0345   out.Close();
0346   for (std::vector<DetId>::const_iterator it = m_barrelCells.begin(); it != m_barrelCells.end(); ++it) {
0347     m_barrelMap[*it] *= InterRings[m_xtalRing[it->rawId()]];
0348   }
0349   for (std::vector<DetId>::const_iterator it = m_endcapCells.begin(); it != m_endcapCells.end(); ++it)
0350     m_endcapMap[*it] *= InterRings[m_xtalRing[it->rawId()]];
0351   if (iCounter < m_loops - 1)
0352     return kContinue;
0353   else
0354     return kStop;
0355 }
0356 
0357 //---------------------------------------
0358 
0359 //LP endOfJob
0360 void InvRingCalib::endOfJob() {
0361   edm::LogInfo("IML") << "[InvMatrixCalibLooper][endOfJob] saving calib coeffs";
0362   calibXMLwriter barrelWriter(EcalBarrel);
0363   calibXMLwriter endcapWriter(EcalEndcap);
0364   for (std::vector<DetId>::const_iterator barrelIt = m_barrelCells.begin(); barrelIt != m_barrelCells.end();
0365        ++barrelIt) {
0366     EBDetId eb(*barrelIt);
0367     barrelWriter.writeLine(eb, m_barrelMap[eb]);
0368   }
0369   for (std::vector<DetId>::const_iterator endcapIt = m_endcapCells.begin(); endcapIt != m_endcapCells.end();
0370        ++endcapIt) {
0371     EEDetId ee(*endcapIt);
0372     endcapWriter.writeLine(ee, m_endcapMap[ee]);
0373   }
0374 }
0375 
0376 //------------------------------------//
0377 //      definition of functions       //
0378 //------------------------------------//
0379 
0380 //------------------------------------------------------------
0381 
0382 //!EE ring definition
0383 void InvRingCalib::EERingDef(const edm::EventSetup& iSetup) {
0384   //Gets the geometry of the endcap
0385   const auto& geometry = iSetup.getData(m_geometryToken);
0386   const CaloSubdetectorGeometry* endcapGeometry = geometry.getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0387   //for every xtal gets the position Vector and the phi position
0388 
0389   // for (std::vector<DetId>::const_iterator barrelIt = m_barrelCells.begin();
0390   //    barrelIt!=m_barrelCells.end();
0391   //    ++barrelIt) {
0392   //     const CaloCellGeometry *cellGeometry = barrelGeometry->getGeometry(*barrelIt);
0393   //     GlobalPoint point;
0394   //     EBDetId eb (*barrelIt);
0395   //     point=cellGeometry->getPosition();
0396   //     m_eta[eb.ieta()]=point.eta() ;    //cellGeometry->getPosition().eta();
0397   //    }
0398   for (std::vector<DetId>::const_iterator endcapIt = m_endcapCells.begin(); endcapIt != m_endcapCells.end();
0399        ++endcapIt) {
0400     auto cellGeometry = endcapGeometry->getGeometry(*endcapIt);
0401     m_cellPos[endcapIt->rawId()] = cellGeometry->getPosition();
0402     m_cellPhi[endcapIt->rawId()] = cellGeometry->getPosition().phi();
0403   }
0404   //takes the first 39 xtals at a fixed y varying the x coordinate and saves their eta coordinate
0405   float eta_ring[39];
0406   for (int ring = 0; ring < 39; ring++)
0407     if (EEDetId::validDetId(ring, 50, 1)) {
0408       EEDetId det = EEDetId(ring, 50, 1, EEDetId::XYMODE);
0409       eta_ring[ring] = m_cellPos[det.rawId()].eta();
0410     }
0411   //defines the bonduary of the rings as the average eta of a xtal
0412   double etaBonduary[40];
0413   etaBonduary[0] = 1.49;
0414   etaBonduary[39] = 4.0;
0415   for (int ring = 1; ring < 39; ++ring)
0416     etaBonduary[ring] = (eta_ring[ring] + eta_ring[ring - 1]) / 2.;
0417   //assign to each xtal a ring number
0418   int CRing;
0419   for (int ring = 0; ring < 39; ring++)
0420     for (std::vector<DetId>::const_iterator endcapIt = m_endcapCells.begin(); endcapIt != m_endcapCells.end();
0421          ++endcapIt) {
0422       if (fabs(m_cellPos[endcapIt->rawId()].eta()) > etaBonduary[ring] &&
0423           fabs(m_cellPos[endcapIt->rawId()].eta()) < etaBonduary[ring + 1]) {
0424         EEDetId ee(*endcapIt);
0425         if (ee.zside() > 0)
0426           CRing = ring + 86;
0427         else
0428           CRing = ring + 125;
0429         m_xtalRing[endcapIt->rawId()] = CRing;
0430         //              m_eta[CRing]=m_cellPos[endcapIt->rawId()].eta();
0431       }
0432     }
0433   return;
0434 }
0435 
0436 //------------------------------------------------------------
0437 
0438 //!Gives the Id of the region given the id of the xtal
0439 int InvRingCalib::EERegId(int id) {
0440   int reg;
0441   int ring;
0442   EEDetId ee(id);
0443   //sets the reg to -1 if the ring doesn't exist or is outside the region of interest
0444   if (m_xtalRing[id] == -1)
0445     return -1;
0446   //avoid the calibration in the wrong zside
0447   if (m_EEZone == 1) {
0448     if (ee.zside() < 0)
0449       return -1;
0450     ring = m_xtalRing[id] - 86;
0451     if (ring >= m_endRing)
0452       return -1;
0453     if (ring < m_startRing)
0454       return -1;
0455     reg = (ring - m_startRing) / m_etaWidth;
0456     m_RinginRegion[id] = (ring - m_startRing) % m_etaWidth;
0457     return reg;
0458   }
0459   if (m_EEZone == -1) {
0460     if (ee.zside() > 0)
0461       return -1;
0462     ring = m_xtalRing[id] - 125;
0463     if (ring >= m_endRing)
0464       return -1;
0465     if (ring < m_startRing)
0466       return -1;
0467     reg = (ring - m_startRing) / m_etaWidth;
0468     m_RinginRegion[id] = (ring - m_startRing) % m_etaWidth;
0469     return reg;
0470   }
0471   if (ee.zside() > 0)
0472     ring = m_xtalRing[id] - 86;
0473   else
0474     ring = m_xtalRing[id] - 125;
0475   if (ring >= m_endRing)
0476     return -1;
0477   if (ring < m_startRing)
0478     return -1;
0479   reg = (ring - m_startRing) / m_etaWidth;
0480   m_RinginRegion[id] = (ring - m_startRing) % m_etaWidth;
0481   return reg;
0482 }
0483 //----------------------------------------
0484 //!Loops over all the endcap xtals and sets for each xtal the value of the region
0485 //!the xtal is in, and the ringNumber inside the region
0486 void InvRingCalib::EERegionDef() {
0487   int reg;
0488   for (std::vector<DetId>::const_iterator endcapIt = m_endcapCells.begin(); endcapIt != m_endcapCells.end();
0489        ++endcapIt) {
0490     EEDetId ee(*endcapIt);
0491     reg = EERegId(endcapIt->rawId());
0492     //If the ring is not of interest saves only the region Id(-1)
0493     if (reg == -1)
0494       m_xtalRegionId[endcapIt->rawId()] = reg;
0495     //sums the number of region in EB or EB+EE to have different regionsId in different regions
0496     else {
0497       if (ee.zside() > 0)
0498         reg += EBRegionNum();
0499       else
0500         reg += EBRegionNum() + EERegionNum();
0501       m_xtalRegionId[endcapIt->rawId()] = reg;
0502     }
0503   }
0504 }
0505 
0506 //------------------------------------------------------------
0507 
0508 //!Number of Regions in EE
0509 inline int InvRingCalib::EERegionNum() const { return ((m_endRing - m_startRing) / m_etaWidth); }
0510 
0511 //! number of Ring in EB
0512 int InvRingCalib::EBRegionNum() const {
0513   if ((m_etaEnd * m_etaStart) > 0)
0514     return ((m_etaEnd - m_etaStart) / m_etaWidth);
0515 
0516   if ((m_etaEnd * m_etaStart) < 0)
0517     return ((m_etaEnd - m_etaStart - 1) / m_etaWidth);
0518 
0519   return 0;
0520 }
0521 //!Divides the barrel in region, necessary to take into
0522 //! account the missing 0 xtal
0523 void InvRingCalib::RegPrepare() {
0524   int k = 0;
0525   for (int i = m_etaStart; i < m_etaEnd; ++i) {
0526     if (i == 0)
0527       continue;
0528     m_Reg[i] = k / m_etaWidth;
0529     ++k;
0530   }
0531 }
0532 //! gives the region Id given ieta
0533 int InvRingCalib::EBRegId(const int ieta) {
0534   if (ieta < m_etaStart || ieta >= m_etaEnd)
0535     return -1;
0536   else
0537     return (m_Reg[ieta]);
0538 }
0539 
0540 //------------------------------------------------------------
0541 
0542 //EB Region Definition
0543 void InvRingCalib::EBRegionDef() {
0544   RegPrepare();
0545   for (std::vector<DetId>::const_iterator it = m_barrelCells.begin(); it != m_barrelCells.end(); ++it) {
0546     EBDetId eb(it->rawId());
0547     m_xtalRing[eb.rawId()] = eb.ieta();
0548     m_xtalRegionId[eb.rawId()] = EBRegId(eb.ieta());
0549     if (m_xtalRegionId[eb.rawId()] == -1)
0550       continue;
0551     m_RinginRegion[eb.rawId()] = (eb.ieta() - m_etaStart) % m_etaWidth;
0552   }
0553 }
0554 //------------------------------------------------------------
[ 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!