Skip to content

Commit

Permalink
Merge pull request #26209 from bsunanda/Run2-alca152
Browse files Browse the repository at this point in the history 
Run2-alca152 Update the macros used for calibration of HCAL
  • Loading branch information
@cmsbuild
cmsbuild committed Mar 25, 2019
2 parents c1e72b4 + bbc7e9f commit 17b243b
Show file tree
Hide file tree
Showing 4 changed files with 86 additions and 40 deletions.
74 changes: 59 additions & 15 deletions Calibration/HcalCalibAlgos/macros/CalibCorr.C
View file
Expand Up @@ -67,25 +67,69 @@ unsigned int truncateId(unsigned int detId, int truncateFlag, bool debug=false){
double puFactor(int type, int ieta, double pmom, double eHcal, double ediff) {

double fac(1.0);
double frac = (type == 1) ? 0.02 : 0.03;
if (pmom > 0 && ediff > frac*pmom) {
double a1(0), a2(0);
if (type == 1) {
a1 = -0.35; a2 = -0.65;
if (std::abs(ieta) == 25) {
a2 = -0.30;
} else if (std::abs(ieta) > 25) {
a1 = -0.45; a2 = -0.10;
if (type <=2) {
double frac = (type == 1) ? 0.02 : 0.03;
if (pmom > 0 && ediff > frac*pmom) {
double a1(0), a2(0);
if (type == 1) {
a1 = -0.35; a2 = -0.65;
if (std::abs(ieta) == 25) {
a2 = -0.30;
} else if (std::abs(ieta) > 25) {
a1 = -0.45; a2 = -0.10;
}
} else {
a1 = -0.39; a2 = -0.59;
if (std::abs(ieta) >= 25) {
a1 = -0.283; a2 = -0.272;
} else if (std::abs(ieta) > 22) {
a1 = -0.238; a2 = -0.241;
}
}
fac = (1.0+a1*(eHcal/pmom)*(ediff/pmom)*(1+a2*(ediff/pmom)));
} else {
a1 = -0.39; a2 = -0.59;
if (std::abs(ieta) >= 25) {
a1 = -0.283; a2 = -0.272;
} else if (std::abs(ieta) > 22) {
a1 = -0.238; a2 = -0.241;
int jeta = std::abs(ieta);
double d2p = (ediff/pmom);
const double DELTA_CUT = 0.03;
const int PU_IETA_3 = 25;
if (type == 3) { // 16pu
const double CONST_COR_COEF[4] = { 0.971, 1.008, 0.985, 1.086 };
const double LINEAR_COR_COEF[4] = { 0, -0.359, -0.251, -0.535 };
const double SQUARE_COR_COEF[4] = { 0, 0, 0.048, 0.143 };
const int PU_IETA_1 = 9;
const int PU_IETA_2 = 16;
unsigned icor = (unsigned(jeta >= PU_IETA_1) +
unsigned(jeta >= PU_IETA_2) +
unsigned(jeta >= PU_IETA_3));
if (d2p > DELTA_CUT) fac = (CONST_COR_COEF[icor] +
LINEAR_COR_COEF[icor]*d2p +
SQUARE_COR_COEF[icor]*d2p*d2p);
} else if (type == 4) { // 17pu
const double CONST_COR_COEF[4] = { 0.974, 1.023, 0.989, 1.077 };
const double LINEAR_COR_COEF[4] = { 0, -0.524, -0.268, -0.584 };
const double SQUARE_COR_COEF[4] = { 0, 0, 0.053, 0.170 };
const int PU_IETA_1 = 9;
const int PU_IETA_2 = 18;
unsigned icor = (unsigned(jeta >= PU_IETA_1) +
unsigned(jeta >= PU_IETA_2) +
unsigned(jeta >= PU_IETA_3));
if (d2p > DELTA_CUT) fac = (CONST_COR_COEF[icor] +
LINEAR_COR_COEF[icor]*d2p +
SQUARE_COR_COEF[icor]*d2p*d2p);
} else { // 18pu
const double CONST_COR_COEF[4] = { 0.973, 0.998, 0.992, 0.965 };
const double LINEAR_COR_COEF[4] = { 0, -0.318, -0.261, -0.406 };
const double SQUARE_COR_COEF[4] = { 0, 0, 0.047, 0.089 };
const int PU_IETA_1 = 7;
const int PU_IETA_2 = 16;
unsigned icor = (unsigned(jeta >= PU_IETA_1) +
unsigned(jeta >= PU_IETA_2) +
unsigned(jeta >= PU_IETA_3));
if (d2p > DELTA_CUT) fac = (CONST_COR_COEF[icor] +
LINEAR_COR_COEF[icor]*d2p +
SQUARE_COR_COEF[icor]*d2p*d2p);
}
}
fac = (1.0+a1*(eHcal/pmom)*(ediff/pmom)*(1+a2*(ediff/pmom)));
}
return fac;
}
Expand Down
27 changes: 5 additions & 22 deletions Calibration/HcalCalibAlgos/macros/CalibFitPlots.C
View file
Expand Up @@ -154,6 +154,8 @@
#include <fstream>
#include <cstdlib>

#include "CalibCorr.C"

struct cfactors {
int ieta, depth;
double corrf, dcorr;
Expand Down Expand Up @@ -214,25 +216,6 @@ std::pair<double,double> GetWidth(TH1D* hist, double xmin, double xmax) {
return std::pair<double,double>(rms,err);
}

std::vector <std::string> splitString (const std::string& fLine) {
std::vector <std::string> result;
int start = 0;
bool empty = true;
for (unsigned i = 0; i <= fLine.size (); i++) {
if (fLine [i] == ' ' || i == fLine.size ()) {
if (!empty) {
std::string item (fLine, start, i-start);
result.push_back (item);
empty = true;
}
start = i+1;
} else {
if (empty) empty = false;
}
}
return result;
}

Double_t langaufun(Double_t *x, Double_t *par) {

//Fit parameters:
Expand Down Expand Up @@ -371,7 +354,7 @@ results fitTwoGauss (TH1D* hist, bool debug) {
startvalues[4] = Fit->Value(1); lowValue[4] = 0.5*startvalues[4]; highValue[4] = 2.*startvalues[4];
startvalues[5] = 2.0*Fit->Value(2); lowValue[5] = 0.5*startvalues[5]; highValue[5] = 100.*startvalues[5];
//fitrange[0] = mean - 2.0*rms; fitrange[1] = mean + 2.0*rms;
fitrange[0] = Fit->Value(1) - 2.0*Fit->Value(2); fitrange[1] = Fit->Value(1) + 2.0*Fit->Value(2);
fitrange[0] = Fit->Value(1) - Fit->Value(2); fitrange[1] = Fit->Value(1) + Fit->Value(2);
TFitResultPtr Fitfun = functionFit(hist, fitrange, startvalues, lowValue, highValue);
double wt1 = (Fitfun->Value(0))*(Fitfun->Value(2));
double value1 = Fitfun->Value(1);
Expand Down Expand Up @@ -660,7 +643,7 @@ void FitHistExtended(const char* infile, const char* outfile,std::string prefix,
}
if (hist0->GetEntries() > 10) {
double rms;
results meaner0 = fitTwoGauss(hist0, debug);
results meaner0 = fitOneGauss(hist0,true,debug);
std::pair<double,double> meaner1 = GetMean(hist0,0.2,2.0,rms);
std::pair<double,double> meaner2 = GetWidth(hist0,0.2,2.0);
if (debug) {
Expand Down Expand Up @@ -704,7 +687,7 @@ void FitHistExtended(const char* infile, const char* outfile,std::string prefix,
value = meaner.mean; error = meaner.errmean;
width = meaner.width; werror= meaner.errwidth;
} else {
results meaner = fitOneGauss(hist,false,debug);
results meaner = fitOneGauss(hist,true,debug);
value = meaner.mean; error = meaner.errmean;
width = meaner.width; werror= meaner.errwidth;
}
Expand Down
23 changes: 21 additions & 2 deletions Calibration/HcalCalibAlgos/macros/CalibMonitor.C
View file
Expand Up @@ -796,12 +796,16 @@ void CalibMonitor::Loop() {
unsigned int kp1 = ps_.size() - 1;
unsigned int kv1 = 0;
std::vector<int> kounts(kp1,0);
std::vector<int> kount50(20,0);
for (Long64_t jentry=0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
if (jentry%100000 == 0) std::cout << "Entry " << jentry << " Run " << t_Run
<< " Event " << t_Event << std::endl;
double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
bool p4060= ((pmom >= 40.0) && (pmom <= 60.0));
if (p4060) ++kount50[0];
bool select = (std::find(entries_.begin(),entries_.end(),jentry) == entries_.end());
if (!select) {
++duplicate;
Expand All @@ -810,11 +814,14 @@ void CalibMonitor::Loop() {
<< " " << t_p << std::endl;
continue;
}
if (p4060) ++kount50[1];
bool selRun = (includeRun_ ? ((t_Run >= runlo_) && (t_Run <= runhi_)) :
((t_Run < runlo_) || (t_Run > runhi_)));
if (select && p4060) ++kount50[2];
select = (selRun && (fabs(t_ieta) >= etalo_) &&
(fabs(t_ieta) <= etahi_) && (t_nVtx >= nvxlo_) &&
(t_nVtx <= nvxhi_));
if (select && p4060) ++kount50[3];
if (!select) {
if (debug)
std::cout << "Run # " << t_Run << " out of range of " << runlo_ << ":"
Expand All @@ -825,18 +832,20 @@ void CalibMonitor::Loop() {
}
if (cSelect_ != nullptr) {
if (exclude_) {
if (cSelect_->isItRBX(t_DetIds)) continue;
if (cSelect_->isItRBX(t_DetIds)) continue;
} else {
if (!(cSelect_->isItRBX(t_ieta,t_iphi))) continue;
}
}
if (p4060) ++kount50[4];
select = (!cutL1T_ || (t_mindR1 >= 0.5));
if (!select) {
if (debug)
std::cout << "Reject Run # " << t_Run << " Event # " << t_Event
<< " too close to L1 trigger " << t_mindR1 << std::endl;
continue;
}
if (p4060) ++kount50[5];
select = ((events_.size() == 0) ||
(std::find(events_.begin(),events_.end(),
std::pair<int,int>(t_Run,t_Event)) != events_.end()));
Expand All @@ -846,10 +855,10 @@ void CalibMonitor::Loop() {
<< " not in the selection list" << std::endl;
continue;
}
if (p4060) ++kount50[6];

// if (Cut(ientry) < 0) continue;
int kp(-1), jp(-1), jp1(-1);
double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
for (unsigned int k=1; k<ps_.size(); ++k ) {
if (pmom >= ps_[k-1] && pmom < ps_[k]) {
kp = k - 1; break;
Expand Down Expand Up @@ -911,7 +920,9 @@ void CalibMonitor::Loop() {
}
}
bool goodTk = goodTrack(eHcal, cut, debug);
if (goodTk && p4060) ++kount50[7];
bool selPhi = selectPhi(debug);
if (goodTk && selPhi && p4060) ++kount50[8];
if (pmom > 0) rat = (eHcal/(pmom-t_eMipDR));
if (debug) {
std::cout << "Entry " << jentry << " p|eHcal|ratio " << pmom << "|"
Expand All @@ -924,6 +935,7 @@ void CalibMonitor::Loop() {
<< ":" << kd << ":" << kd1 << ":" << jp << std::endl;
}
if (goodTk && kp >=0 && selPhi) {
if (p4060) ++kount50[9];
if (t_eHcal < 0.01) {
std::map<int,counter>::const_iterator itr = runEn1.find(t_Run);
if (itr == runEn1.end()) {
Expand All @@ -939,6 +951,7 @@ void CalibMonitor::Loop() {
}
}
if (t_eMipDR < 0.01 && t_eHcal < 0.01) {
if (p4060) ++kount50[10];
std::map<int,counter>::const_iterator itr = runEn2.find(t_Run);
if (itr == runEn2.end()) {
counter knt;
Expand All @@ -953,6 +966,7 @@ void CalibMonitor::Loop() {
}
}
if (rat > rcut) {
if (p4060) ++kount50[11];
if (plotType_ <= 1) {
h_etaX[kp][kv]->Fill(eta,rat,t_EventWeight);
h_etaX[kp][kv1]->Fill(eta,rat,t_EventWeight);
Expand Down Expand Up @@ -983,6 +997,7 @@ void CalibMonitor::Loop() {
}
}
if ((!dataMC_) || (t_mindR1 > 0.5) || (t_DataType == 1)) {
if (p4060) ++kount50[12];
++kounts[kp];
if (plotType_ <= 1) {
if (jp > 0) h_etaF[kp][jp]->Fill(rat,t_EventWeight);
Expand Down Expand Up @@ -1017,6 +1032,7 @@ void CalibMonitor::Loop() {
h_etaB[kp][jp1]->Fill(rat,t_EventWeight);
h_etaB[kp][jp2]->Fill(rat,t_EventWeight);
}
if (p4060) ++kount50[13];
}
}
}
Expand Down Expand Up @@ -1066,6 +1082,9 @@ void CalibMonitor::Loop() {
for (unsigned int k=1; k<ps_.size(); ++k)
if (ps_[k] > 21) std::cout << ps_[k-1] <<":"<< ps_[k] << " "
<< kounts[k-1] << std::endl;
std::cout << "Number in each step: ";
for (unsigned int k=0; k < 14; ++k) std::cout << " [" << k << "] " << kount50[k];
std::cout <<std::endl;
}

bool CalibMonitor::goodTrack(double& eHcal, double &cuti, bool debug) {
Expand Down
2 changes: 1 addition & 1 deletion Calibration/HcalCalibAlgos/macros/CalibTree.C
View file
Expand Up @@ -551,7 +551,7 @@ Double_t CalibTree::Loop(int loop, TFile *fout, bool useweight, int nMin,
else isItRBX = !(temp);
}
++kprint;
if (cSelect_ && !(isItRBX)) {
if (!(isItRBX)) {
for (unsigned int idet=0; idet<(*t_DetIds).size(); idet++) {
if (selectPhi((*t_DetIds)[idet])) {
unsigned int detid = truncateId((*t_DetIds)[idet],
Expand Down

0 comments on commit 17b243b

Please sign in to comment.