Skip to content

Commit

Permalink
Modify most of the scripts needed to get IsoTrack calibration of HCAL
Browse files Browse the repository at this point in the history 
Code check
  • Loading branch information
Sunanda committed Mar 2, 2023
1 parent 35e58d0 commit b4cd7b3
Show file tree
Hide file tree
Showing 5 changed files with 175 additions and 57 deletions.
158 changes: 131 additions & 27 deletions Calibration/HcalCalibAlgos/macros/CalibCorr.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,30 @@ void unpackDetId(unsigned int detId, int& subdet, int& zside, int& ieta, int& ip
}
}

unsigned int packDetId(int subdet, int ieta, int iphi, int depth) {
// The maskings are defined in DataFormats/DetId/interface/DetId.h
// and in DataFormats/HcalDetId/interface/HcalDetId.h
// The macro does not invoke the classes there and use them
const int det(4);
unsigned int id = (((det & 0xF) << 28) | ((subdet & 0x7) << 25));
id |= ((0x1000000) | ((depth & 0xF) << 20) | ((ieta > 0) ? (0x80000 | (ieta << 10)) : ((-ieta) << 10)) |
(iphi & 0x3FF));
return id;
}

unsigned int matchDetId(unsigned int detId) {
if ((detId & 0x1000000) == 0) {
int subdet = ((detId >> 25) & (0x7));
int ieta = ((detId >> 7) & 0x3F);
int zside = (detId & 0x2000) ? (1) : (-1);
int depth = ((detId >> 14) & 0x1F);
int iphi = (detId & 0x3F);
return packDetId(subdet, (ieta * zside), iphi, depth);
} else {
return detId;
}
}

unsigned int truncateId(unsigned int detId, int truncateFlag, bool debug = false) {
//Truncate depth information of DetId's
unsigned int id(detId);
Expand Down Expand Up @@ -453,24 +477,29 @@ public:
float getCorr(int run, unsigned int id);
double getCorr(const Long64_t& entry);
double getTrueCorr(const Long64_t& entry);
double getPhiCorr(unsigned int id);
bool absent(const Long64_t& entry);
bool absent() { return (good_ == 0); }
bool present(const Long64_t& entry);

private:
void readCorrRun(const char* infile);
void readCorrDepth(const char* infile);
void readCorrResp(const char* infile);
void readCorrPU(const char* infile);
unsigned int readCorrRun(const char* infile);
unsigned int readCorrDepth(const char* infile);
unsigned int readCorrResp(const char* infile);
unsigned int readCorrPU(const char* infile);
unsigned int readCorrPhi(const char* infile);
unsigned int getDetIdHE(int ieta, int iphi, int depth);
unsigned int getDetId(int subdet, int ieta, int iphi, int depth);
unsigned int correctDetId(const unsigned int& detId);

static const unsigned int nmax_ = 10;
int flag_;
bool debug_;
unsigned int good_;
std::map<unsigned int, float> corrFac_[nmax_], corrFacDepth_, corrFacResp_;
std::map<Long64_t, double> cfactors_;
std::vector<int> runlow_;
std::map<unsigned int, double> corrPhiSym_;
};

class CalibSelectRBX {
Expand Down Expand Up @@ -510,17 +539,19 @@ CalibCorrFactor::CalibCorrFactor(const char* infile, int useScale, double scale,
double CalibCorrFactor::getCorr(unsigned int id) {
double cfac(1.0);
if (corrE_) {
int subdet, zside, ieta, iphi, depth;
unpackDetId(id, subdet, zside, ieta, iphi, depth);
std::map<std::pair<int, int>, double>::const_iterator itr =
cfactors_.find(std::pair<int, int>(zside * ieta, depth));
if (itr != cfactors_.end()) {
cfac = itr->second;
} else if (etaMax_) {
if (zside > 0 && ieta > etamp_)
cfac = (depth < depMax_) ? cfacmp_[depth] : cfacmp_[depMax_ - 1];
if (zside < 0 && ieta > -etamn_)
cfac = (depth < depMax_) ? cfacmn_[depth] : cfacmn_[depMax_ - 1];
if (cfactors_.size() > 0) {
int subdet, zside, ieta, iphi, depth;
unpackDetId(id, subdet, zside, ieta, iphi, depth);
std::map<std::pair<int, int>, double>::const_iterator itr =
cfactors_.find(std::pair<int, int>(zside * ieta, depth));
if (itr != cfactors_.end()) {
cfac = itr->second;
} else if (etaMax_) {
if (zside > 0 && ieta > etamp_)
cfac = (depth < depMax_) ? cfacmp_[depth] : cfacmp_[depMax_ - 1];
if (zside < 0 && ieta > -etamn_)
cfac = (depth < depMax_) ? cfacmn_[depth] : cfacmn_[depMax_ - 1];
}
}
} else if (useScale_ != 0) {
int subdet, zside, ieta, iphi, depth;
Expand Down Expand Up @@ -590,17 +621,21 @@ double CalibCorrFactor::getFactor(const int& ieta) {
CalibCorr::CalibCorr(const char* infile, int flag, bool debug) : flag_(flag), debug_(debug) {
std::cout << "CalibCorr is created with flag " << flag << ":" << flag_ << " for i/p file " << infile << std::endl;
if (flag == 1)
readCorrDepth(infile);
good_ = readCorrDepth(infile);
else if (flag == 2)
readCorrResp(infile);
good_ = readCorrResp(infile);
else if (flag == 3)
readCorrPU(infile);
good_ = readCorrPU(infile);
else if (flag == 3)
good_ = readCorrPhi(infile);
else
readCorrRun(infile);
good_ = readCorrRun(infile);
}

float CalibCorr::getCorr(int run, unsigned int id) {
float cfac(1.0);
if (good_ == 0)
return cfac;
unsigned idx = correctDetId(id);
if (flag_ == 1) {
std::map<unsigned int, float>::iterator itr = corrFacDepth_.find(idx);
Expand All @@ -610,6 +645,10 @@ float CalibCorr::getCorr(int run, unsigned int id) {
std::map<unsigned int, float>::iterator itr = corrFacResp_.find(idx);
if (itr != corrFacResp_.end())
cfac = itr->second;
} else if (flag_ == 4) {
std::map<unsigned int, double>::iterator itr = corrPhiSym_.find(idx);
if (itr != corrPhiSym_.end())
cfac = itr->second;
} else {
int ip(-1);
for (unsigned int k = 0; k < runlow_.size(); ++k) {
Expand Down Expand Up @@ -638,6 +677,8 @@ float CalibCorr::getCorr(int run, unsigned int id) {
}

double CalibCorr::getCorr(const Long64_t& entry) {
if (good_ == 0)
return 1.0;
double cfac(0.0);
std::map<Long64_t, double>::iterator itr = cfactors_.find(entry);
if (itr != cfactors_.end())
Expand All @@ -646,26 +687,43 @@ double CalibCorr::getCorr(const Long64_t& entry) {
}

double CalibCorr::getTrueCorr(const Long64_t& entry) {
if (good_ == 0)
return 1.0;
double cfac(0.0);
std::map<Long64_t, double>::iterator itr = cfactors_.find(entry);
if (itr != cfactors_.end())
cfac = itr->second;
return cfac;
}

double CalibCorr::getPhiCorr(unsigned int idx) {
double cfac(1.0);
if (good_ == 0)
return cfac;
std::map<unsigned int, double>::iterator itr = corrPhiSym_.find(idx);
if (itr != corrPhiSym_.end())
cfac = itr->second;
if (debug_) {
int subdet, zside, ieta, iphi, depth;
unpackDetId(idx, subdet, zside, ieta, iphi, depth);
std::cout << "ID " << std::hex << idx << std::dec << " (Sub " << subdet << " eta " << zside * ieta << " phi "
<< iphi << " depth " << depth << ") Factor " << cfac << std::endl;
}
return cfac;
}

bool CalibCorr::absent(const Long64_t& entry) { return (cfactors_.find(entry) == cfactors_.end()); }

bool CalibCorr::present(const Long64_t& entry) { return (cfactors_.find(entry) != cfactors_.end()); }

void CalibCorr::readCorrRun(const char* infile) {
unsigned int CalibCorr::readCorrRun(const char* infile) {
std::cout << "Enters readCorrRun for " << infile << std::endl;
std::ifstream fInput(infile);
unsigned int ncorr(0);
unsigned int all(0), good(0), ncorr(0);
if (!fInput.good()) {
std::cout << "Cannot open file " << infile << std::endl;
} else {
char buffer[1024];
unsigned int all(0), good(0);
while (fInput.getline(buffer, 1024)) {
++all;
std::string bufferString(buffer);
Expand Down Expand Up @@ -711,16 +769,17 @@ void CalibCorr::readCorrRun(const char* infile) {
std::cout << "Reads total of " << all << " and " << good << " good records of run dependent corrections from "
<< infile << std::endl;
}
return good;
}

void CalibCorr::readCorrDepth(const char* infile) {
unsigned int CalibCorr::readCorrDepth(const char* infile) {
std::cout << "Enters readCorrDepth for " << infile << std::endl;
unsigned int all(0), good(0);
std::ifstream fInput(infile);
if (!fInput.good()) {
std::cout << "Cannot open file " << infile << std::endl;
} else {
char buffer[1024];
unsigned int all(0), good(0);
while (fInput.getline(buffer, 1024)) {
++all;
std::string bufferString(buffer);
Expand Down Expand Up @@ -755,16 +814,17 @@ void CalibCorr::readCorrDepth(const char* infile) {
std::cout << "Reads total of " << all << " and " << good << " good records of depth dependent factors from "
<< infile << std::endl;
}
return good;
}

void CalibCorr::readCorrResp(const char* infile) {
unsigned int CalibCorr::readCorrResp(const char* infile) {
std::cout << "Enters readCorrResp for " << infile << std::endl;
unsigned int all(0), good(0), other(0);
std::ifstream fInput(infile);
if (!fInput.good()) {
std::cout << "Cannot open file " << infile << std::endl;
} else {
char buffer[1024];
unsigned int all(0), good(0), other(0);
while (fInput.getline(buffer, 1024)) {
++all;
std::string bufferString(buffer);
Expand Down Expand Up @@ -800,9 +860,10 @@ void CalibCorr::readCorrResp(const char* infile) {
std::cout << "Reads total of " << all << " and " << good << " good and " << other
<< " detector records of depth dependent factors from " << infile << std::endl;
}
return good;
}

void CalibCorr::readCorrPU(const char* infile) {
unsigned int CalibCorr::readCorrPU(const char* infile) {
if (std::string(infile) != "") {
std::ifstream fInput(infile);
if (!fInput.good()) {
Expand All @@ -821,6 +882,49 @@ void CalibCorr::readCorrPU(const char* infile) {
}
}
std::cout << "Reads " << cfactors_.size() << " PU correction factors from " << infile << std::endl;
return cfactors_.size();
}

unsigned int CalibCorr::readCorrPhi(const char* infile) {
std::cout << "Enters readCorrPhi for " << infile << std::endl;
unsigned int all(0), good(0);
std::ifstream fInput(infile);
if (!fInput.good()) {
std::cout << "Cannot open file " << infile << std::endl;
} else {
char buffer[1024];
while (fInput.getline(buffer, 1024)) {
++all;
std::string bufferString(buffer);
if (bufferString.substr(0, 1) == "#") {
continue; //ignore other comments
} else {
std::vector<std::string> items = splitString(bufferString);
if (items.size() < 5) {
std::cout << "Ignore line: " << buffer << " Size " << items.size();
for (unsigned int k = 0; k < items.size(); ++k)
std::cout << " [" << k << "] : " << items[k];
std::cout << std::endl;
} else {
++good;
int subdet = std::atoi(items[0].c_str());
int ieta = std::atoi(items[1].c_str());
int iphi = std::atoi(items[2].c_str());
int depth = std::atoi(items[3].c_str());
double corrf = std::atof(items[4].c_str());
unsigned int id = packDetId(subdet, ieta, iphi, depth);
corrPhiSym_[id] = corrf;
if (debug_)
std::cout << "ID " << std::hex << id << std::dec << ":" << id << " (subdet " << subdet << " eta " << ieta
<< " phi " << iphi << " depth " << depth << ") " << corrPhiSym_[id] << std::endl;
}
}
}
fInput.close();
std::cout << "Reads total of " << all << " and " << good << " good records of phi-symmetry factors from " << infile
<< std::endl;
}
return good;
}

unsigned int CalibCorr::getDetIdHE(int ieta, int iphi, int depth) { return getDetId(2, ieta, iphi, depth); }
Expand Down
17 changes: 10 additions & 7 deletions Calibration/HcalCalibAlgos/macros/CalibMonitor.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,16 +39,17 @@
// rcorFileName (char*) = name of the text file having the correction
// factors as a function of run numbers or depth
// or entry number to be used for raddam/depth/
// pileup dependent correction (default="",
// no corr.)
// pileup/phisym dependent correction
// (default="", no correction)
// puCorr (int) = PU correction to be applied or not: 0 no
// correction; < 0 use eDelta; > 0 rho dependent
// correction (-8)
// flag (int) = 7 digit integer (xmlthdo) with control
// information (x=3/2/1/0 for having 1000/500/50/
// 100 bins for response distribution in (0:5);
// m=1/0 for (not) making plots for each RBX;
// l=3/2/1/0 for type of rcorFileName (3 for
// l=4/3/2/1/0 for type of rcorFileName (4 for
// using results from phi-symmetry; 3 for
// pileup correction using machine learning
// method; 2 for overall response corrections;
// 1 for depth dependence corrections;
Expand Down Expand Up @@ -396,8 +397,10 @@ CalibMonitor::CalibMonitor(const char *fname,
Init(chain, dupFileName, comFileName, outFName);
if (std::string(rcorFileName) != "") {
cFactor_ = new CalibCorr(rcorFileName, ifDepth_, false);
if (cFactor_->absent())
ifDepth_ = -1;
} else {
ifDepth_ = 0;
ifDepth_ = -1;
}
if (rbx != 0)
cSelect_ = new CalibSelectRBX(rbx, false);
Expand Down Expand Up @@ -1118,7 +1121,7 @@ void CalibMonitor::Loop(Long64_t nmax) {
unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
cfac = corrFactor_->getCorr(id);
}
if ((cFactor_ != nullptr) && (ifDepth_ != 3))
if ((cFactor_ != nullptr) && (ifDepth_ != 3) && (ifDepth_ > 0))
cfac *= cFactor_->getCorr(t_Run, (*t_DetIds)[k]);
eHcal += (cfac * ((*t_HitEnergies)[k]));
if (debug) {
Expand Down Expand Up @@ -1698,15 +1701,15 @@ void CalibMonitor::correctEnergy(double &eHcal, const Long64_t &entry) {
for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
unsigned int id = truncateId((*t_DetIds1)[idet], truncateFlag_, false);
double cfac = corrFactor_->getCorr(id);
if (cFactor_ != 0)
if ((cFactor_ != 0) && (ifDepth_ != 3) && (ifDepth_ > 0))
cfac *= cFactor_->getCorr(t_Run, (*t_DetIds1)[idet]);
double hitEn = cfac * (*t_HitEnergies1)[idet];
Etot1 += hitEn;
}
for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
unsigned int id = truncateId((*t_DetIds3)[idet], truncateFlag_, false);
double cfac = corrFactor_->getCorr(id);
if (cFactor_ != 0)
if ((cFactor_ != 0) && (ifDepth_ != 3) && (ifDepth_ > 0))
cfac *= cFactor_->getCorr(t_Run, (*t_DetIds3)[idet]);
double hitEn = cfac * (*t_HitEnergies3)[idet];
Etot3 += hitEn;
Expand Down
Loading

0 comments on commit b4cd7b3

Please sign in to comment.