Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-09-07 04:53:47

0001 //  TSelector-based code for getting the HCAL resp. correction
0002 //  from physics events. Works for DiJet and IsoTrack calibration.
0003 //
0004 //  Anton Anastassov (Northwestern)
0005 //  Email: aa@fnal.gov
0006 //
0007 //
0008 
0009 #include "Calibration/HcalCalibAlgos/interface/hcalCalib.h"
0010 #include "Calibration/HcalCalibAlgos/interface/hcalCalibUtils.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 
0013 #include <iostream>
0014 #include <fstream>
0015 #include <sstream>
0016 #include <string>
0017 
0018 #include <map>
0019 #include <numeric>
0020 #include <algorithm>
0021 #include <set>
0022 
0023 #include "Calibration/Tools/interface/MinL3AlgoUniv.h"
0024 
0025 #include "DataFormats/DetId/interface/DetId.h"
0026 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0027 
0028 void hcalCalib::Begin(TTree* /*tree*/) {
0029   TString option = GetOption();
0030 
0031   nEvents = 0;
0032 
0033   if (APPLY_PHI_SYM_COR_FLAG && !ReadPhiSymCor()) {
0034     edm::LogError("HcalCalib") << "\nERROR: Failed to read the phi symmetry corrections.\n"
0035                                << "Check if the filename is correct. If the corrections are not needed, set the "
0036                                   "corresponding flag to \"false\"\n"
0037                                << "\nThe program will be terminated\n";
0038 
0039     exit(1);
0040   }
0041 
0042   //  cellEnergies.reserve(1000000);
0043   //  cellIds.reserve(1000000);
0044   //  targetEnergies.reserve(1000000);
0045 
0046   histoFile = new TFile(HISTO_FILENAME.Data(), "RECREATE");
0047 
0048   h1_trkP = new TH1F("h1_trkP", "Track momenta; p_{trk} (GeV); Number of tracks", 100, 0, 200);
0049   h1_allTrkP = new TH1F("h1_allTrkP", "Track momenta - all tracks; p_{trk} (GeV); Number of tracks", 100, 0, 200);
0050 
0051   h1_selTrkP_iEta10 = new TH1F(
0052       "h1_selTrkP_iEta10", "Track momenta - tracks with |iEta|<10; p_{trk} (GeV); Number of tracks", 100, 0, 200);
0053 
0054   if (CALIB_TYPE == "ISO_TRACK")
0055     h1_rawSumE = new TH1F("h1_rawSumE", "Cluster Energy; E_{cl} (GeV); Number of tracks", 100, 0, 200);
0056   else
0057     h1_rawSumE = new TH1F("h1_rawSumE", "Cluster Energy; E_{cl} (GeV); Number of tracks", 1000, 0, 2000);
0058 
0059   h1_rawResp = new TH1F("h1_rawResp", "Uncorrected response: |iEta|<24;  E_{had}/p; Number of tracks", 300, 0, 3);
0060   h1_corResp = new TH1F("h1_corResp", "Corrected response: |iEta|<24; E_{had}/p; Number of tracks", 300, 0, 3);
0061 
0062   h1_rawRespBarrel =
0063       new TH1F("h1_rawRespBarrel", "Uncorrected response: |iEta|<15;  E_{had}/p; Number of tracks", 300, 0, 3);
0064   h1_corRespBarrel =
0065       new TH1F("h1_corRespBarrel", "Corrected response: |iEta|<15; E_{had}/p; Number of tracks", 300, 0, 3);
0066 
0067   h1_rawRespEndcap =
0068       new TH1F("h1_rawRespEndcap", "Uncorrected response:  17<|iEta|<24;  E_{had}/p; Number of tracks", 300, 0, 3);
0069   h1_corRespEndcap =
0070       new TH1F("h1_corRespEndcap", "Corrected response: 17<|iEta|<24; E_{had}/p; Number of tracks", 300, 0, 3);
0071 
0072   h1_numEventsTwrIEta = new TH1F("h1_numEventsTwrIEta", "h1_numEventsTwrIEta", 80, -40, 40);
0073 
0074   h2_dHitRefBarrel = new TH2F("h2_dHitRefBarrel",
0075                               "{#Delta}i{#phi} vs {#Delta}i{#eta} of hit and most energetic "
0076                               "tower(|i{#eta}|<16);{#Delta}i{#eta}; {#Delta}i{#phi}",
0077                               10,
0078                               -5,
0079                               5,
0080                               10,
0081                               -5,
0082                               5);
0083   h2_dHitRefEndcap = new TH2F("h2_dHitRefEndcap",
0084                               "{#Delta}i{#phi} vs {#Delta}i{#eta} of hit and most energetic tower (16<|i{#eta}|<25) "
0085                               ";{#Delta}i{#eta}; {#Delta}i{#phi}",
0086                               10,
0087                               -5,
0088                               5,
0089                               10,
0090                               -5,
0091                               5);
0092 
0093   TString histoName = "isoTrack_";
0094 
0095   for (Int_t i = 0; i < 48; ++i) {
0096     Long_t iEta;
0097     if (i < 24)
0098       iEta = i - 24;
0099     else
0100       iEta = i - 23;
0101     TString hn = histoName + iEta;
0102     h1_corRespIEta[i] = new TH1F(hn, hn, 300, 0, 3.0);
0103   }
0104 
0105 }  // end of Begin()
0106 
0107 //void hcalCalib::SlaveBegin(TTree * /*tree*/) {
0108 //  TString option = GetOption();
0109 //}
0110 
0111 Bool_t hcalCalib::Process(Long64_t entry) {
0112   //  fChain->GetTree()->GetEntry(entry);
0113   GetEntry(entry);
0114 
0115   std::set<UInt_t> uniqueIds;  // for testing: check if there are duplicate cells   (AA)
0116 
0117   Bool_t acceptEvent = kTRUE;
0118 
0119   ++nEvents;
0120 
0121   if (!(nEvents % 100000))
0122     edm::LogVerbatim("HcalCalib") << "event: " << nEvents;
0123 
0124   h1_allTrkP->Fill(targetE);
0125 
0126   if (targetE < MIN_TARGET_E || targetE > MAX_TARGET_E)
0127     return kFALSE;
0128   ;
0129 
0130   // make local copy as the cells may be modified  due to phi/depth sum, phi corrections etc
0131   std::vector<TCell> selectCells;
0132 
0133   if (cells->GetSize() == 0)
0134     return kFALSE;
0135 
0136   if (CALIB_TYPE == "DI_JET" && probeJetEmFrac > 0.999)
0137     return kTRUE;
0138 
0139   for (Int_t i = 0; i < cells->GetSize(); ++i) {
0140     TCell* thisCell = (TCell*)cells->At(i);
0141 
0142     if (HcalDetId(thisCell->id()).subdet() == HcalOuter)
0143       continue;  // reject HO, make a switch!
0144 
0145     if (HcalDetId(thisCell->id()).subdet() != HcalBarrel && HcalDetId(thisCell->id()).subdet() != HcalEndcap &&
0146         HcalDetId(thisCell->id()).subdet() != HcalForward) {
0147       edm::LogWarning("HcalCalib") << "Unknown or wrong hcal subdetector: " << HcalDetId(thisCell->id()).subdet();
0148     }
0149 
0150     // Apply phi symmetry corrections if the flag is set
0151     if (APPLY_PHI_SYM_COR_FLAG)
0152       thisCell->SetE(phiSymCor[thisCell->id()] * thisCell->e());
0153 
0154     if (thisCell->e() > MIN_CELL_E)
0155       selectCells.push_back(*thisCell);
0156   }
0157 
0158   if (selectCells.empty()) {
0159     edm::LogWarning("HcalCalib") << "NO CELLS ABOVE THRESHOLD FOUND FOR TARGET!!!";
0160   }
0161 
0162   if (SUM_DEPTHS)
0163     sumDepths(selectCells);
0164   else if (SUM_SMALL_DEPTHS)
0165     sumSmallDepths(selectCells);  // depth 1,2 in twrs 15,16
0166 
0167   // most energetic tower (IsoTracks) or centroid of probe jet (DiJets)
0168   std::pair<Int_t, UInt_t> refPos;
0169 
0170   Int_t dEtaHitRef = 999;
0171   Int_t dPhiHitRef = 999;
0172 
0173   if (CALIB_TYPE == "ISO_TRACK") {
0174     Int_t iEtaMaxE;   // filled by reference in getIEtaIPhiForHighestE
0175     UInt_t iPhiMaxE;  //
0176 
0177     getIEtaIPhiForHighestE(selectCells, iEtaMaxE, iPhiMaxE);
0178 
0179     dEtaHitRef = iEtaMaxE - iEtaHit;
0180     dPhiHitRef = iPhiMaxE - iPhiHit;
0181 
0182     if (dPhiHitRef < -36)
0183       dPhiHitRef += 72;
0184     if (dPhiHitRef > 36)
0185       dPhiHitRef -= 72;
0186 
0187     if (iEtaHit * iEtaMaxE < 0) {
0188       if (dEtaHitRef < 0)
0189         dEtaHitRef += 1;
0190       if (dEtaHitRef > 0)
0191         dEtaHitRef -= 1;
0192     }
0193 
0194     if (abs(iEtaHit) < 16)
0195       h2_dHitRefBarrel->Fill(dEtaHitRef, dPhiHitRef);
0196     if (abs(iEtaHit) > 16 && abs(iEtaHit) < 25)
0197       h2_dHitRefEndcap->Fill(dEtaHitRef, dPhiHitRef);
0198 
0199     // --------------------------------------------------
0200     // Choice of cluster definition
0201     //
0202     // fixed size NxN clusters as specified in to config file
0203     if (!USE_CONE_CLUSTERING) {
0204       if (abs(iEtaMaxE) < 16 && HB_CLUSTER_SIZE == 3)
0205         filterCells3x3(selectCells, iEtaMaxE, iPhiMaxE);
0206       if (abs(iEtaMaxE) > 15 && HE_CLUSTER_SIZE == 3)
0207         filterCells3x3(selectCells, iEtaMaxE, iPhiMaxE);
0208 
0209       if (abs(iEtaMaxE) < 16 && HB_CLUSTER_SIZE == 5)
0210         filterCells5x5(selectCells, iEtaMaxE, iPhiMaxE);
0211       if (abs(iEtaMaxE) > 15 && HE_CLUSTER_SIZE == 5)
0212         filterCells5x5(selectCells, iEtaMaxE, iPhiMaxE);
0213     } else {
0214       //  calculate distance at hcal surface
0215       const GlobalPoint hitPositionHcal(xTrkHcal, yTrkHcal, zTrkHcal);
0216       filterCellsInCone(selectCells, hitPositionHcal, MAX_CONE_DIST, theCaloGeometry);
0217     }
0218 
0219     refPos.first = iEtaMaxE;
0220     refPos.second = iPhiMaxE;
0221 
0222   } else if (CALIB_TYPE == "DI_JET") {  // Apply selection cuts on DiJet events here
0223 
0224     if (etVetoJet > MAX_ET_THIRD_JET)
0225       acceptEvent = kFALSE;
0226 
0227     Float_t jetsDPhi = probeJetP4->DeltaPhi(*tagJetP4);
0228     if (fabs(jetsDPhi * 180.0 / M_PI) < MIN_DPHI_DIJETS)
0229       acceptEvent = kFALSE;
0230 
0231     if (probeJetEmFrac > MAX_PROBEJET_EMFRAC)
0232       acceptEvent = kFALSE;
0233     if (fabs(probeJetP4->Eta()) < MIN_PROBEJET_ABSETA)
0234       acceptEvent = kFALSE;
0235     if (fabs(tagJetP4->Eta()) > MAX_TAGJET_ABSETA)
0236       acceptEvent = kFALSE;
0237     if (fabs(tagJetP4->Et()) < MIN_TAGJET_ET)
0238       acceptEvent = kFALSE;
0239 
0240     if (acceptEvent) {
0241       Int_t iEtaMaxE;   // filled by reference in getIEtaIPhiForHighestE
0242       UInt_t iPhiMaxE;  //
0243 
0244       getIEtaIPhiForHighestE(selectCells, iEtaMaxE, iPhiMaxE);
0245 
0246       // The ref position for the jet is not used in the minimization at this time.
0247       // It will be needed if we attempt to do matrix inversion: then the question is
0248       // which value is better suited: the centroid of the jet or the hottest tower...
0249 
0250       //    refPos.first  = iEtaHit;
0251       //    refPos.second = iPhiHit;
0252 
0253       refPos.first = iEtaMaxE;
0254       refPos.second = iPhiMaxE;
0255 
0256       if (abs(iEtaMaxE) > 40)
0257         acceptEvent = kFALSE;  // for testing :  set as parameter (AA)
0258     }
0259   }
0260 
0261   if (COMBINE_PHI)
0262     combinePhi(selectCells);
0263 
0264   // fill the containers for the minimization prcedures
0265   std::vector<Float_t> energies;
0266   std::vector<UInt_t> ids;
0267 
0268   for (std::vector<TCell>::iterator i_it = selectCells.begin(); i_it != selectCells.end(); ++i_it) {
0269     // for testing : fill only unique id's
0270 
0271     if (uniqueIds.insert(i_it->id()).second) {
0272       energies.push_back(i_it->e());
0273       ids.push_back(i_it->id());
0274     }
0275   }
0276 
0277   if (CALIB_TYPE == "ISO_TRACK") {
0278     if (accumulate(energies.begin(), energies.end(), 0.0) / targetE < MIN_EOVERP)
0279       acceptEvent = kFALSE;
0280     if (accumulate(energies.begin(), energies.end(), 0.0) / targetE > MAX_EOVERP)
0281       acceptEvent = kFALSE;
0282 
0283     if (emEnergy > MAX_TRK_EME)
0284       acceptEvent = kFALSE;
0285 
0286     if (abs(dEtaHitRef) > 1 || abs(dPhiHitRef) > 1)
0287       acceptEvent = kFALSE;
0288 
0289     // Have to check if for |iEta|>20 (and neighboring region) the dPhiHitRef
0290     // should be relaxed to 2. The neighboring towers have dPhi=2...
0291   }
0292 
0293   h1_rawSumE->Fill(accumulate(energies.begin(), energies.end(), 0.0));
0294 
0295   // here we fill the information for the minimization procedure
0296   if (acceptEvent) {
0297     cellEnergies.push_back(energies);
0298     cellIds.push_back(ids);
0299     targetEnergies.push_back(targetE);
0300     refIEtaIPhi.push_back(refPos);
0301 
0302     if (abs(refPos.first) <= 10)
0303       h1_selTrkP_iEta10->Fill(targetE);
0304   }
0305 
0306   // Clean up
0307   energies.clear();
0308   ids.clear();
0309   selectCells.clear();
0310 
0311   return kTRUE;
0312 }
0313 
0314 //void hcalCalib::SlaveTerminate() {}
0315 
0316 void hcalCalib::Terminate() {
0317   edm::LogVerbatim("HcalCalib") << "\n\nFinished reading the events.\n"
0318                                 << "Number of input objects: " << cellIds.size()
0319                                 << "\nPerforming minimization: depending on selected method can take some time...\n\n";
0320 
0321   for (std::vector<std::pair<Int_t, UInt_t> >::iterator it_rp = refIEtaIPhi.begin(); it_rp != refIEtaIPhi.end();
0322        ++it_rp) {
0323     Float_t weight = (abs(it_rp->first) < 21) ? 1.0 / 72.0 : 1.0 / 36.0;
0324     h1_numEventsTwrIEta->Fill(it_rp->first, weight);
0325   }
0326 
0327   if (CALIB_METHOD == "MATRIX_INV_OF_ETA_AVE") {
0328     GetCoefFromMtrxInvOfAve();
0329   } else if (CALIB_METHOD == "L3" || CALIB_METHOD == "L3_AND_MTRX_INV") {
0330     int eventWeight = 2;  // 2 is the default (try at some point 0,1,2,3)
0331     MinL3AlgoUniv<UInt_t>* thisL3Algo = new MinL3AlgoUniv<UInt_t>(eventWeight);
0332     int numIterations = 10;  // default 10
0333 
0334     solution = thisL3Algo->iterate(cellEnergies, cellIds, targetEnergies, numIterations);
0335 
0336     // in order to handle the case where sumDepths="false", but the flag to sum depths 1,2 in HB towers 15, 16
0337     // is set (sumSmallDepths) we create entries in "solution" to create equal correction here
0338     // for each encountered coef in depth one.
0339 
0340     if (!SUM_DEPTHS && SUM_SMALL_DEPTHS) {
0341       std::vector<UInt_t> idForSummedCells;
0342 
0343       for (std::map<UInt_t, Float_t>::iterator m_it = solution.begin(); m_it != solution.end(); ++m_it) {
0344         if (HcalDetId(m_it->first).ietaAbs() != 15 && HcalDetId(m_it->first).ietaAbs() != 16)
0345           continue;
0346         if (HcalDetId(m_it->first).subdet() != HcalBarrel)
0347           continue;
0348         if (HcalDetId(m_it->first).depth() == 1)
0349           idForSummedCells.push_back(HcalDetId(m_it->first));
0350       }
0351 
0352       for (std::vector<UInt_t>::iterator v_it = idForSummedCells.begin(); v_it != idForSummedCells.end(); ++v_it) {
0353         UInt_t addCoefId = HcalDetId(HcalBarrel, HcalDetId(*v_it).ieta(), HcalDetId(*v_it).iphi(), 2);
0354         solution[addCoefId] = solution[*v_it];
0355       }
0356 
0357     }  // end of special treatment for "sumSmallDepths" mode
0358 
0359     if (CALIB_METHOD == "L3_AND_MTRX_INV") {
0360       GetCoefFromMtrxInvOfAve();
0361 
0362       // loop over the solution from L3 and multiply by the additional factor from
0363       // the matrix inversion. Set coef outside of the valid calibration region =1.
0364       for (std::map<UInt_t, Float_t>::iterator it_s = solution.begin(); it_s != solution.end(); ++it_s) {
0365         Int_t iEtaSol = HcalDetId(it_s->first).ieta();
0366         if (abs(iEtaSol) < CALIB_ABS_IETA_MIN || abs(iEtaSol) > CALIB_ABS_IETA_MAX)
0367           it_s->second = 1.0;
0368         else
0369           it_s->second *= iEtaCoefMap[iEtaSol];
0370       }
0371 
0372     }  // if CALIB_METHOD=="L3_AND_MTRX_INV"
0373 
0374   }  // end of L3 or L3 + mtrx inversion minimization
0375 
0376   // done getting the constants -> write the formatted file
0377   makeTextFile();
0378 
0379   // fill some histograms
0380   Float_t rawResp = 0;
0381   Float_t corResp = 0;
0382   Int_t maxIEta = 0;
0383   Int_t minIEta = 999;
0384 
0385   for (unsigned int i = 0; i < cellEnergies.size(); ++i) {
0386     Int_t iEta;
0387     for (unsigned int j = 0; j < (cellEnergies[i]).size(); ++j) {
0388       iEta = HcalDetId(cellIds[i][j]).ieta();
0389       rawResp += (cellEnergies[i])[j];
0390 
0391       if (CALIB_METHOD == "L3_AND_MTRX_INV") {
0392         corResp += (solution[cellIds[i][j]] * cellEnergies[i][j]);
0393       } else if (CALIB_METHOD == "L3") {
0394         corResp += (solution[cellIds[i][j]] * (cellEnergies[i][j]));
0395       } else if (CALIB_METHOD == "MATRIX_INV_OF_ETA_AVE") {
0396         corResp += (iEtaCoefMap[iEta] * cellEnergies[i][j]);
0397       }
0398 
0399       if (maxIEta < abs(iEta))
0400         maxIEta = abs(iEta);
0401       if (minIEta > abs(iEta))
0402         minIEta = abs(iEta);
0403     }
0404 
0405     rawResp /= targetEnergies[i];
0406     corResp /= targetEnergies[i];
0407 
0408     // fill histograms based on iEta on ref point of the cluster (for now: hottest tower)
0409     // expect range |iEta|<=24 (to do: add flexibility for arbitrary range)
0410 
0411     if (CALIB_TYPE == "ISO_TRACK") {
0412       Int_t ind = refIEtaIPhi[i].first;
0413       (ind < 0) ? (ind += 24) : (ind += 23);
0414       if (ind >= 0 && ind < 48) {
0415         h1_corRespIEta[ind]->Fill(corResp);
0416       }
0417 
0418       // fill histograms for cases where all towers are in the barrel or endcap
0419       if (maxIEta < 25) {
0420         h1_rawResp->Fill(rawResp);
0421         h1_corResp->Fill(corResp);
0422       }
0423       if (maxIEta < 15) {
0424         h1_rawRespBarrel->Fill(rawResp);
0425         h1_corRespBarrel->Fill(corResp);
0426       } else if (maxIEta < 25 && minIEta > 16) {
0427         h1_rawRespEndcap->Fill(rawResp);
0428         h1_corRespEndcap->Fill(corResp);
0429       }
0430     }  // histograms for isotrack calibration
0431 
0432     else {
0433       // put jet plots here
0434     }
0435 
0436     rawResp = 0;
0437     corResp = 0;
0438     maxIEta = 0;
0439     minIEta = 999;
0440   }
0441 
0442   // save the histograms
0443   h1_trkP->Write();
0444   h1_allTrkP->Write();
0445 
0446   h1_selTrkP_iEta10->Write();
0447 
0448   h1_rawSumE->Write();
0449   h1_rawResp->Write();
0450   h1_corResp->Write();
0451   h1_rawRespBarrel->Write();
0452   h1_corRespBarrel->Write();
0453   h1_rawRespEndcap->Write();
0454   h1_corRespEndcap->Write();
0455   h1_numEventsTwrIEta->Write();
0456   h2_dHitRefBarrel->Write();
0457   h2_dHitRefEndcap->Write();
0458   for (Int_t i = 0; i < 48; ++i) {
0459     h1_corRespIEta[i]->Write();
0460   }
0461 
0462   histoFile->Write();
0463   histoFile->Close();
0464 
0465   edm::LogVerbatim("HcalCalib") << "\n Finished calibration.\n ";
0466 
0467 }  // end of Terminate()
0468 
0469 //**************************************************************************************************************
0470 
0471 void hcalCalib::GetCoefFromMtrxInvOfAve() {
0472   // these maps are keyed by iEta
0473   std::map<Int_t, Float_t> aveTargetE;
0474   std::map<Int_t, Int_t> nEntries;  // count hits
0475 
0476   //  iEtaRef  iEtaCell, energy
0477   std::map<Int_t, std::map<Int_t, Float_t> > aveHitE;  // add energies in the loop, normalize after that
0478 
0479   for (unsigned int i = 0; i < cellEnergies.size(); ++i) {
0480     Int_t iEtaRef = refIEtaIPhi[i].first;
0481     aveTargetE[iEtaRef] += targetEnergies[i];
0482     nEntries[iEtaRef]++;
0483 
0484     // for hybrid method: matrix inv of averages preceeded by L3
0485     if (CALIB_METHOD == "L3_AND_MTRX_INV") {
0486       for (unsigned int j = 0; j < (cellEnergies[i]).size(); ++j) {
0487         aveHitE[iEtaRef][HcalDetId(cellIds[i][j]).ieta()] += (solution[cellIds[i][j]] * cellEnergies[i][j]);
0488       }
0489     } else if (CALIB_METHOD == "MATRIX_INV_OF_ETA_AVE") {
0490       for (unsigned int j = 0; j < (cellEnergies[i]).size(); ++j) {
0491         aveHitE[iEtaRef][HcalDetId(cellIds[i][j]).ieta()] += cellEnergies[i][j];
0492       }
0493     }
0494 
0495   }  // end of loop of entries
0496 
0497   // scale by number of entries to get the averages
0498   Float_t norm = 1.0;
0499   for (std::map<Int_t, Float_t>::iterator m_it = aveTargetE.begin(); m_it != aveTargetE.end(); ++m_it) {
0500     Int_t iEta = m_it->first;
0501     norm = (nEntries[iEta] > 0) ? 1.0 / (nEntries[iEta]) : 1.0;
0502     aveTargetE[iEta] *= norm;
0503 
0504     std::map<Int_t, Float_t>::iterator n_it = (aveHitE[iEta]).begin();
0505 
0506     Float_t sumRawE = 0;
0507     for (; n_it != (aveHitE[iEta]).end(); ++n_it) {
0508       (n_it->second) *= norm;
0509       sumRawE += (n_it->second);
0510     }
0511 
0512   }  // end of scaling by number of entries
0513 
0514   Int_t ONE_SIDE_IETA_RANGE = CALIB_ABS_IETA_MAX - CALIB_ABS_IETA_MIN + 1;
0515 
0516   // conversion from iEta to index for the linear system
0517   // contains elements only in the valid range for *matrix inversion*
0518   std::vector<Int_t> iEtaList;
0519 
0520   for (Int_t i = -CALIB_ABS_IETA_MAX; i <= CALIB_ABS_IETA_MAX; ++i) {
0521     if (abs(i) < CALIB_ABS_IETA_MIN)
0522       continue;
0523     iEtaList.push_back(i);
0524   }
0525 
0526   TMatrixD A(2 * ONE_SIDE_IETA_RANGE, 2 * ONE_SIDE_IETA_RANGE);
0527   TMatrixD b(2 * ONE_SIDE_IETA_RANGE, 1);
0528   TMatrixD x(2 * ONE_SIDE_IETA_RANGE, 1);
0529 
0530   for (Int_t i = 0; i < 2 * ONE_SIDE_IETA_RANGE; ++i) {
0531     for (Int_t j = 0; j < 2 * ONE_SIDE_IETA_RANGE; ++j) {
0532       A(i, j) = 0.0;
0533     }
0534   }
0535 
0536   for (UInt_t i = 0; i < iEtaList.size(); ++i) {
0537     b(i, 0) = aveTargetE[iEtaList[i]];
0538 
0539     std::map<Int_t, Float_t>::iterator n_it = aveHitE[iEtaList[i]].begin();
0540     for (; n_it != aveHitE[iEtaList[i]].end(); ++n_it) {
0541       if (fabs(n_it->first) > CALIB_ABS_IETA_MAX || fabs(n_it->first) < CALIB_ABS_IETA_MIN)
0542         continue;
0543       Int_t j = Int_t(find(iEtaList.begin(), iEtaList.end(), n_it->first) - iEtaList.begin());
0544       A(i, j) = n_it->second;
0545     }
0546   }
0547 
0548   TMatrixD coef = b;
0549   TDecompQRH qrh(A);
0550   Bool_t hasSolution = qrh.MultiSolve(coef);
0551 
0552   if (hasSolution && (CALIB_METHOD == "L3_AND_MTRX_INV" || CALIB_METHOD == "MATRIX_INV_OF_ETA_AVE")) {
0553     for (UInt_t i = 0; i < iEtaList.size(); ++i) {
0554       iEtaCoefMap[iEtaList[i]] = coef(i, 0);
0555     }
0556   }
0557 }
0558 
0559 Bool_t hcalCalib::ReadPhiSymCor() {
0560   std::ifstream phiSymFile(PHI_SYM_COR_FILENAME.Data());
0561 
0562   if (!phiSymFile) {
0563     edm::LogWarning("HcalCalib") << "\nERROR: Can not find file with phi symmetry constants \""
0564                                  << PHI_SYM_COR_FILENAME.Data() << "\"";
0565     return kFALSE;
0566   }
0567 
0568   // assumes the format used in CSA08, first line is a comment
0569 
0570   Int_t iEta;
0571   UInt_t iPhi;
0572   UInt_t depth;
0573   TString sdName;
0574   UInt_t detId;
0575 
0576   Float_t value;
0577   HcalSubdetector sd;
0578 
0579   std::string line;
0580 
0581   while (getline(phiSymFile, line)) {
0582     if (line.empty() || line[0] == '#')
0583       continue;
0584 
0585     std::istringstream linestream(line);
0586     linestream >> iEta >> iPhi >> depth >> sdName >> value >> std::hex >> detId;
0587 
0588     if (sdName == "HB")
0589       sd = HcalBarrel;
0590     else if (sdName == "HE")
0591       sd = HcalEndcap;
0592     else if (sdName == "HO")
0593       sd = HcalOuter;
0594     else if (sdName == "HF")
0595       sd = HcalForward;
0596     else {
0597       edm::LogWarning("HcalCalib") << "\nInvalid detector name in phi symmetry constants file: " << sdName.Data()
0598                                    << "\nCheck file and rerun!\n";
0599       return kFALSE;
0600     }
0601 
0602     // check if the data is consistent
0603 
0604     if (HcalDetId(sd, iEta, iPhi, depth) != HcalDetId(detId)) {
0605       edm::LogWarning("HcalCalib")
0606           << "\nInconsistent info in phi symmetry file: subdet, iEta, iPhi, depth do not match rawId!\n";
0607       return kFALSE;
0608     }
0609     HcalDetId hId(detId);
0610     if (!topo_->valid(hId)) {
0611       edm::LogWarning("HcalCalib") << "\nInvalid DetId from: iEta=" << iEta << " iPhi=" << iPhi << " depth=" << depth
0612                                    << " subdet=" << sdName.Data() << " detId=" << detId << "\n";
0613       return kFALSE;
0614     }
0615 
0616     phiSymCor[HcalDetId(sd, iEta, iPhi, depth)] = value;
0617   }
0618 
0619   return kTRUE;
0620 }
0621 
0622 void hcalCalib::makeTextFile() {
0623   //{ HcalEmpty=0, HcalBarrel=1, HcalEndcap=2, HcalOuter=3, HcalForward=4, HcalTriggerTower=5, HcalOther=7 };
0624 
0625   TString sdName[8] = {"EMPTY", "HB", "HE", "HO", "HF", "TRITWR", "UNDEF", "OTHER"};
0626 
0627   FILE* constFile = fopen(OUTPUT_COR_COEF_FILENAME.Data(), "w+");
0628 
0629   // header of the constants file
0630   fprintf(constFile, "%1s%16s%16s%16s%16s%9s%11s\n", "#", "eta", "phi", "depth", "det", "value", "DetId");
0631 
0632   // Order loops to produce sequence of constants as for phi symmetry
0633 
0634   for (Int_t sd = 1; sd <= 4; sd++) {
0635     for (Int_t e = 1; e <= 50; e++) {
0636       Int_t eta;
0637 
0638       for (Int_t side = 0; side < 2; side++) {
0639         eta = (side == 0) ? -e : e;  //ta *= (-1);
0640 
0641         for (Int_t phi = 1; phi <= 72; phi++) {
0642           for (Int_t d = 1; d < 5; d++) {
0643             if (!topo_->valid(HcalDetId(HcalSubdetector(sd), eta, phi, d)))
0644               continue;
0645             HcalDetId id(HcalSubdetector(sd), eta, phi, d);
0646             Float_t corrFactor = 1.0;
0647 
0648             if (abs(eta) >= CALIB_ABS_IETA_MIN && abs(eta) <= CALIB_ABS_IETA_MAX && HcalSubdetector(sd) != HcalOuter) {
0649               //        if (abs(eta)>=CALIB_ABS_IETA_MIN && abs(eta)<=22 && HcalSubdetector(sd)!=HcalOuter) {
0650 
0651               // need some care when depths were summed for iEta=16 =>
0652               // the coeficients are saved in depth 1 of HB: affects
0653               Int_t subdetInd = sd;
0654 
0655               if (abs(eta) == 16 && HcalSubdetector(sd) == HcalEndcap && SUM_DEPTHS) {
0656                 subdetInd = 1;
0657               }
0658 
0659               if (CALIB_METHOD == "L3" || CALIB_METHOD == "L3_AND_MTRX_INV") {
0660                 if (SUM_DEPTHS && COMBINE_PHI)
0661                   corrFactor = solution[HcalDetId(HcalSubdetector(subdetInd), eta, 1, 1)];
0662                 else if (SUM_DEPTHS)
0663                   corrFactor = solution[HcalDetId(HcalSubdetector(subdetInd), eta, phi, 1)];
0664                 else if (COMBINE_PHI)
0665                   corrFactor = solution[HcalDetId(HcalSubdetector(sd), eta, 1, d)];
0666                 else
0667                   corrFactor = solution[HcalDetId(HcalSubdetector(sd), eta, phi, d)];
0668 
0669                 // Remark: a new case was added (sumSmallDepths) where the first two depths in towers 15,16
0670                 // are summed and stored in  depth 1.
0671                 // For now we create the correction coef for depth 2 (set equal to depth 1)
0672                 // after the call to the L3 minimizer so that this case is also handled without modifying the
0673                 // logic above. Probably it is better to move it here?
0674 
0675               }  // L3
0676 
0677               else if (CALIB_METHOD == "MATRIX_INV_OF_ETA_AVE") {
0678                 corrFactor = iEtaCoefMap[eta];
0679               }
0680 
0681               else if (CALIB_METHOD == "ISO_TRK_PHI_SYM") {
0682                 corrFactor = solution[HcalDetId(HcalSubdetector(sd), eta, phi, d)];
0683               }
0684 
0685             }  // if within the calibration range
0686 
0687             fprintf(constFile, "%17i%16i%16i%16s%9.5f%11X\n", eta, phi, d, sdName[sd].Data(), corrFactor, id.rawId());
0688           }
0689         }
0690       }
0691     }
0692   }
0693 
0694   return;
0695 }
0696 
0697 inline void hcalCalib::Init(TTree* tree) {
0698   // The Init() function is called when the selector needs to initialize
0699   // a new tree or chain. Typically here the branch addresses and branch
0700   // pointers of the tree will be set.
0701   // It is normaly not necessary to make changes to the generated
0702   // code, but the routine can be extended by the user if needed.
0703   // Init() will be called many times when running on PROOF
0704   // (once per file to be processed).
0705 
0706   // Set object pointer
0707   cells = nullptr;
0708   tagJetP4 = nullptr;
0709   probeJetP4 = nullptr;
0710 
0711   // Set branch addresses and branch pointers
0712   if (!tree)
0713     return;
0714   fChain = tree;
0715 
0716   //      fChain->SetMakeClass(1);
0717 
0718   fChain->SetBranchAddress("eventNumber", &eventNumber, &b_eventNumber);
0719   fChain->SetBranchAddress("runNumber", &runNumber, &b_runNumber);
0720   fChain->SetBranchAddress("iEtaHit", &iEtaHit, &b_iEtaHit);
0721   fChain->SetBranchAddress("iPhiHit", &iPhiHit, &b_iPhiHit);
0722   fChain->SetBranchAddress("cells", &cells, &b_cells);
0723   fChain->SetBranchAddress("emEnergy", &emEnergy, &b_emEnergy);
0724   fChain->SetBranchAddress("targetE", &targetE, &b_targetE);
0725   fChain->SetBranchAddress("etVetoJet", &etVetoJet, &b_etVetoJet);
0726 
0727   fChain->SetBranchAddress("xTrkHcal", &xTrkHcal, &b_xTrkHcal);
0728   fChain->SetBranchAddress("yTrkHcal", &yTrkHcal, &b_yTrkHcal);
0729   fChain->SetBranchAddress("zTrkHcal", &zTrkHcal, &b_zTrkHcal);
0730   fChain->SetBranchAddress("xTrkEcal", &xTrkEcal, &b_xTrkEcal);
0731   fChain->SetBranchAddress("yTrkEcal", &yTrkEcal, &b_yTrkEcal);
0732   fChain->SetBranchAddress("zTrkEcal", &zTrkEcal, &b_zTrkEcal);
0733 
0734   fChain->SetBranchAddress("tagJetEmFrac", &tagJetEmFrac, &b_tagJetEmFrac);
0735   fChain->SetBranchAddress("probeJetEmFrac", &probeJetEmFrac, &b_probeJetEmFrac);
0736 
0737   fChain->SetBranchAddress("tagJetP4", &tagJetP4, &b_tagJetP4);
0738   fChain->SetBranchAddress("probeJetP4", &probeJetP4, &b_probeJetP4);
0739 }
0740 
0741 inline Bool_t hcalCalib::Notify() {
0742   // The Notify() function is called when a new file is opened. This
0743   // can be either for a new TTree in a TChain or when when a new TTree
0744   // is started when using PROOF. It is normaly not necessary to make changes
0745   // to the generated code, but the routine can be extended by the
0746   // user if needed. The return value is currently not used.
0747 
0748   return kTRUE;
0749 }