HcalCalibAlgos/src/hcalCalib.cc

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     for (; n_it != (aveHitE[iEta]).end(); ++n_it) {
0507       (n_it->second) *= norm;
0508     }
0509
0510   }  // end of scaling by number of entries
0511
0512   Int_t ONE_SIDE_IETA_RANGE = CALIB_ABS_IETA_MAX - CALIB_ABS_IETA_MIN + 1;
0513
0514   // conversion from iEta to index for the linear system
0515   // contains elements only in the valid range for *matrix inversion*
0516   std::vector<Int_t> iEtaList;
0517
0518   for (Int_t i = -CALIB_ABS_IETA_MAX; i <= CALIB_ABS_IETA_MAX; ++i) {
0519     if (abs(i) < CALIB_ABS_IETA_MIN)
0520       continue;
0521     iEtaList.push_back(i);
0522   }
0523
0524   TMatrixD A(2 * ONE_SIDE_IETA_RANGE, 2 * ONE_SIDE_IETA_RANGE);
0525   TMatrixD b(2 * ONE_SIDE_IETA_RANGE, 1);
0526   TMatrixD x(2 * ONE_SIDE_IETA_RANGE, 1);
0527
0528   for (Int_t i = 0; i < 2 * ONE_SIDE_IETA_RANGE; ++i) {
0529     for (Int_t j = 0; j < 2 * ONE_SIDE_IETA_RANGE; ++j) {
0530       A(i, j) = 0.0;
0531     }
0532   }
0533
0534   for (UInt_t i = 0; i < iEtaList.size(); ++i) {
0535     b(i, 0) = aveTargetE[iEtaList[i]];
0536
0537     std::map<Int_t, Float_t>::iterator n_it = aveHitE[iEtaList[i]].begin();
0538     for (; n_it != aveHitE[iEtaList[i]].end(); ++n_it) {
0539       if (fabs(n_it->first) > CALIB_ABS_IETA_MAX || fabs(n_it->first) < CALIB_ABS_IETA_MIN)
0540         continue;
0541       Int_t j = Int_t(find(iEtaList.begin(), iEtaList.end(), n_it->first) - iEtaList.begin());
0542       A(i, j) = n_it->second;
0543     }
0544   }
0545
0546   TMatrixD coef = b;
0547   TDecompQRH qrh(A);
0548   Bool_t hasSolution = qrh.MultiSolve(coef);
0549
0550   if (hasSolution && (CALIB_METHOD == "L3_AND_MTRX_INV" || CALIB_METHOD == "MATRIX_INV_OF_ETA_AVE")) {
0551     for (UInt_t i = 0; i < iEtaList.size(); ++i) {
0552       iEtaCoefMap[iEtaList[i]] = coef(i, 0);
0553     }
0554   }
0555 }
0556
0557 Bool_t hcalCalib::ReadPhiSymCor() {
0558   std::ifstream phiSymFile(PHI_SYM_COR_FILENAME.Data());
0559
0560   if (!phiSymFile) {
0561     edm::LogWarning("HcalCalib") << "\nERROR: Can not find file with phi symmetry constants \""
0562                                  << PHI_SYM_COR_FILENAME.Data() << "\"";
0563     return kFALSE;
0564   }
0565
0566   // assumes the format used in CSA08, first line is a comment
0567
0568   Int_t iEta;
0569   UInt_t iPhi;
0570   UInt_t depth;
0571   TString sdName;
0572   UInt_t detId;
0573
0574   Float_t value;
0575   HcalSubdetector sd;
0576
0577   std::string line;
0578
0579   while (getline(phiSymFile, line)) {
0580     if (line.empty() || line[0] == '#')
0581       continue;
0582
0583     std::istringstream linestream(line);
0584     linestream >> iEta >> iPhi >> depth >> sdName >> value >> std::hex >> detId;
0585
0586     if (sdName == "HB")
0587       sd = HcalBarrel;
0588     else if (sdName == "HE")
0589       sd = HcalEndcap;
0590     else if (sdName == "HO")
0591       sd = HcalOuter;
0592     else if (sdName == "HF")
0593       sd = HcalForward;
0594     else {
0595       edm::LogWarning("HcalCalib") << "\nInvalid detector name in phi symmetry constants file: " << sdName.Data()
0596                                    << "\nCheck file and rerun!\n";
0597       return kFALSE;
0598     }
0599
0600     // check if the data is consistent
0601
0602     if (HcalDetId(sd, iEta, iPhi, depth) != HcalDetId(detId)) {
0603       edm::LogWarning("HcalCalib")
0604           << "\nInconsistent info in phi symmetry file: subdet, iEta, iPhi, depth do not match rawId!\n";
0605       return kFALSE;
0606     }
0607     HcalDetId hId(detId);
0608     if (!topo_->valid(hId)) {
0609       edm::LogWarning("HcalCalib") << "\nInvalid DetId from: iEta=" << iEta << " iPhi=" << iPhi << " depth=" << depth
0610                                    << " subdet=" << sdName.Data() << " detId=" << detId << "\n";
0611       return kFALSE;
0612     }
0613
0614     phiSymCor[HcalDetId(sd, iEta, iPhi, depth)] = value;
0615   }
0616
0617   return kTRUE;
0618 }
0619
0620 void hcalCalib::makeTextFile() {
0621   //{ HcalEmpty=0, HcalBarrel=1, HcalEndcap=2, HcalOuter=3, HcalForward=4, HcalTriggerTower=5, HcalOther=7 };
0622
0623   TString sdName[8] = {"EMPTY", "HB", "HE", "HO", "HF", "TRITWR", "UNDEF", "OTHER"};
0624
0625   FILE* constFile = fopen(OUTPUT_COR_COEF_FILENAME.Data(), "w+");
0626
0627   // header of the constants file
0628   fprintf(constFile, "%1s%16s%16s%16s%16s%9s%11s\n", "#", "eta", "phi", "depth", "det", "value", "DetId");
0629
0630   // Order loops to produce sequence of constants as for phi symmetry
0631
0632   for (Int_t sd = 1; sd <= 4; sd++) {
0633     for (Int_t e = 1; e <= 50; e++) {
0634       Int_t eta;
0635
0636       for (Int_t side = 0; side < 2; side++) {
0637         eta = (side == 0) ? -e : e;  //ta *= (-1);
0638
0639         for (Int_t phi = 1; phi <= 72; phi++) {
0640           for (Int_t d = 1; d < 5; d++) {
0641             if (!topo_->valid(HcalDetId(HcalSubdetector(sd), eta, phi, d)))
0642               continue;
0643             HcalDetId id(HcalSubdetector(sd), eta, phi, d);
0644             Float_t corrFactor = 1.0;
0645
0646             if (abs(eta) >= CALIB_ABS_IETA_MIN && abs(eta) <= CALIB_ABS_IETA_MAX && HcalSubdetector(sd) != HcalOuter) {
0647               //        if (abs(eta)>=CALIB_ABS_IETA_MIN && abs(eta)<=22 && HcalSubdetector(sd)!=HcalOuter) {
0648
0649               // need some care when depths were summed for iEta=16 =>
0650               // the coeficients are saved in depth 1 of HB: affects
0651               Int_t subdetInd = sd;
0652
0653               if (abs(eta) == 16 && HcalSubdetector(sd) == HcalEndcap && SUM_DEPTHS) {
0654                 subdetInd = 1;
0655               }
0656
0657               if (CALIB_METHOD == "L3" || CALIB_METHOD == "L3_AND_MTRX_INV") {
0658                 if (SUM_DEPTHS && COMBINE_PHI)
0659                   corrFactor = solution[HcalDetId(HcalSubdetector(subdetInd), eta, 1, 1)];
0660                 else if (SUM_DEPTHS)
0661                   corrFactor = solution[HcalDetId(HcalSubdetector(subdetInd), eta, phi, 1)];
0662                 else if (COMBINE_PHI)
0663                   corrFactor = solution[HcalDetId(HcalSubdetector(sd), eta, 1, d)];
0664                 else
0665                   corrFactor = solution[HcalDetId(HcalSubdetector(sd), eta, phi, d)];
0666
0667                 // Remark: a new case was added (sumSmallDepths) where the first two depths in towers 15,16
0668                 // are summed and stored in  depth 1.
0669                 // For now we create the correction coef for depth 2 (set equal to depth 1)
0670                 // after the call to the L3 minimizer so that this case is also handled without modifying the
0671                 // logic above. Probably it is better to move it here?
0672
0673               }  // L3
0674
0675               else if (CALIB_METHOD == "MATRIX_INV_OF_ETA_AVE") {
0676                 corrFactor = iEtaCoefMap[eta];
0677               }
0678
0679               else if (CALIB_METHOD == "ISO_TRK_PHI_SYM") {
0680                 corrFactor = solution[HcalDetId(HcalSubdetector(sd), eta, phi, d)];
0681               }
0682
0683             }  // if within the calibration range
0684
0685             fprintf(constFile, "%17i%16i%16i%16s%9.5f%11X\n", eta, phi, d, sdName[sd].Data(), corrFactor, id.rawId());
0686           }
0687         }
0688       }
0689     }
0690   }
0691
0692   return;
0693 }
0694
0695 inline void hcalCalib::Init(TTree* tree) {
0696   // The Init() function is called when the selector needs to initialize
0697   // a new tree or chain. Typically here the branch addresses and branch
0698   // pointers of the tree will be set.
0699   // It is normaly not necessary to make changes to the generated
0700   // code, but the routine can be extended by the user if needed.
0701   // Init() will be called many times when running on PROOF
0702   // (once per file to be processed).
0703
0704   // Set object pointer
0705   cells = nullptr;
0706   tagJetP4 = nullptr;
0707   probeJetP4 = nullptr;
0708
0709   // Set branch addresses and branch pointers
0710   if (!tree)
0711     return;
0712   fChain = tree;
0713
0714   //      fChain->SetMakeClass(1);
0715
0716   fChain->SetBranchAddress("eventNumber", &eventNumber, &b_eventNumber);
0717   fChain->SetBranchAddress("runNumber", &runNumber, &b_runNumber);
0718   fChain->SetBranchAddress("iEtaHit", &iEtaHit, &b_iEtaHit);
0719   fChain->SetBranchAddress("iPhiHit", &iPhiHit, &b_iPhiHit);
0720   fChain->SetBranchAddress("cells", &cells, &b_cells);
0721   fChain->SetBranchAddress("emEnergy", &emEnergy, &b_emEnergy);
0722   fChain->SetBranchAddress("targetE", &targetE, &b_targetE);
0723   fChain->SetBranchAddress("etVetoJet", &etVetoJet, &b_etVetoJet);
0724
0725   fChain->SetBranchAddress("xTrkHcal", &xTrkHcal, &b_xTrkHcal);
0726   fChain->SetBranchAddress("yTrkHcal", &yTrkHcal, &b_yTrkHcal);
0727   fChain->SetBranchAddress("zTrkHcal", &zTrkHcal, &b_zTrkHcal);
0728   fChain->SetBranchAddress("xTrkEcal", &xTrkEcal, &b_xTrkEcal);
0729   fChain->SetBranchAddress("yTrkEcal", &yTrkEcal, &b_yTrkEcal);
0730   fChain->SetBranchAddress("zTrkEcal", &zTrkEcal, &b_zTrkEcal);
0731
0732   fChain->SetBranchAddress("tagJetEmFrac", &tagJetEmFrac, &b_tagJetEmFrac);
0733   fChain->SetBranchAddress("probeJetEmFrac", &probeJetEmFrac, &b_probeJetEmFrac);
0734
0735   fChain->SetBranchAddress("tagJetP4", &tagJetP4, &b_tagJetP4);
0736   fChain->SetBranchAddress("probeJetP4", &probeJetP4, &b_probeJetP4);
0737 }
0738
0739 inline Bool_t hcalCalib::Notify() {
0740   // The Notify() function is called when a new file is opened. This
0741   // can be either for a new TTree in a TChain or when when a new TTree
0742   // is started when using PROOF. It is normaly not necessary to make changes
0743   // to the generated code, but the routine can be extended by the
0744   // user if needed. The return value is currently not used.
0745
0746   return kTRUE;
0747 }