Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Run2-alca117 Modify all macros for IsoTrack Calibration #22133

Merged
merged 2 commits into from
Feb 8, 2018
Merged

Run2-alca117 Modify all macros for IsoTrack Calibration #22133

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6fc38e6
Modify all macros for IsoTrack Calibration
Feb 6, 2018
94ae008
change 0 to nullptr
Feb 7, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
38 changes: 32 additions & 6 deletions Calibration/HcalCalibAlgos/macros/CalibCorr.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,31 @@ void unpackDetId(unsigned int detId, int& subdet, int& zside, int& ieta,
}
}

unsigned int truncateId(unsigned int detId, int truncateFlag, bool debug=false){
//Truncate depth information of DetId's
unsigned int id(detId);
if (debug) {
std::cout << "Truncate 1 " << std::hex << detId << " " << id
<< std::dec << std::endl;
}
int subdet, depth, zside, ieta, iphi;
unpackDetId(detId, subdet, zside, ieta, iphi, depth);
if (truncateFlag == 1) {
//Ignore depth index of ieta values of 15 and 16 of HB
if ((subdet == 1) && (ieta > 14)) depth = 1;
} else if (truncateFlag == 2) {
//Ignore depth index of all ieta values
depth = 1;
}
id = (subdet<<25) | (0x1000000) | ((depth&0xF)<<20) | ((zside>0)?(0x80000|(ieta<<10)):(ieta<<10));
if (debug) {
std::cout << "Truncate 2: " << subdet << " " << zside*ieta << " "
<< depth << " " << std::hex << id << " input " << detId
<< std::dec << std::endl;
}
return id;
}

class CalibCorr {
public :
CalibCorr(const char* infile, bool debug=false);
Expand Down Expand Up @@ -72,7 +97,9 @@ float CalibCorr::getCorr(int run, unsigned int id) {
ip = (int)(i); break;
}
}
if (debug_) std::cout << "Run " << run << " Perdiod " << ip << std::endl;
if (debug_) {
std::cout << "Run " << run << " Perdiod " << ip << std::endl;
}
unsigned idx = correctDetId(id);
if (ip >= 0) {
std::map<unsigned int,float>::iterator itr = corrFac_[ip].find(idx);
Expand Down Expand Up @@ -107,7 +134,8 @@ void CalibCorr::readCorr(const char* infile) {
int run = std::atoi (items[n].c_str());
runlow_.push_back(run);
}
std::cout << ncorr << ":" << runlow_.size() << " Run ranges\n";
std::cout << ncorr << ":" << runlow_.size() << " Run ranges"
<< std::endl;
for (unsigned int n=0; n<runlow_.size(); ++n)
std::cout << " [" << n << "] " << runlow_[n];
std::cout << std::endl;
Expand Down Expand Up @@ -244,10 +272,8 @@ bool CalibSelectRBX::isItRBX(const std::vector<unsigned int> * detId) {
if (ok) break;
}
}
if (debug_) {
std::cout << "isItRBX: size " << detId->size() << " OK " << ok
<< std::endl;
}
if (debug_) std::cout << "isItRBX: size " << detId->size() << " OK " << ok
<< std::endl;
return ok;
}

Expand Down
167 changes: 105 additions & 62 deletions Calibration/HcalCalibAlgos/macros/CalibFitPlots.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -196,10 +196,11 @@ std::pair<double,double> fitTwoGauss (TH1D* hist, bool debug) {
g1->SetLineColor(kGreen);
TFitResultPtr Fit = hist->Fit(g1,option.c_str(),"");

if (debug)
if (debug) {
for (int k=0; k<3; ++k)
std::cout << "Initial Parameter[" << k << "] = " << Fit->Value(k)
<< " +- " << Fit->FitResult::Error(k) << std::endl;
}
double startvalues[6], fitrange[2], lowValue[6], highValue[6];
startvalues[0] = Fit->Value(0); lowValue[0] = 0.5*startvalues[0]; highValue[0] = 2.*startvalues[0];
startvalues[1] = Fit->Value(1); lowValue[1] = 0.5*startvalues[1]; highValue[1] = 2.*startvalues[1];
Expand Down Expand Up @@ -240,14 +241,17 @@ std::pair<double,double> fitOneGauss (TH1D* hist, bool debug) {
double value = Fit->Value(1);
double error = Fit->FitResult::Error(1);
std::pair<double,double> meaner = GetMean(hist,0.2,2.0);
if (debug) std::cout << "Fit " << value << ":" << error << ":"
<< hist->GetMeanError() << " Mean "
<< meaner.first << ":" << meaner.second;
if (debug) {
std::cout << "Fit " << value << ":" << error << ":" << hist->GetMeanError()
<< " Mean " << meaner.first << ":" << meaner.second;
}
double minvalue(0.30);
if (value < minvalue || value > 2.0 || error > 0.5) {
value = meaner.first; error = meaner.second;
}
if (debug) std::cout << " Final " << value << ":" << error << std::endl;
if (debug) {
std::cout << " Final " << value << ":" << error << std::endl;
}
return std::pair<double,double>(value,error);
}

Expand Down Expand Up @@ -362,8 +366,10 @@ void FitHistStandard(std::string infile,std::string outfile,std::string prefix,
double LowEdge = histo->GetBinLowEdge(jmin);
double HighEdge= histo->GetBinLowEdge(jmax)+histo->GetBinWidth(jmax);
TFitResultPtr Fit = histo->Fit("pol0","+QRWLS","",LowEdge,HighEdge);
if (debug) std::cout << "Fit to Pol0: " << Fit->Value(0) << " +- "
<< Fit->FitResult::Error(0) << std::endl;
if (debug) {
std::cout << "Fit to Pol0: " << Fit->Value(0) << " +- "
<< Fit->FitResult::Error(0) << std::endl;
}
histo->GetXaxis()->SetTitle(xname[m1].c_str());
histo->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
histo->GetYaxis()->SetRangeUser(0.4,1.6);
Expand Down Expand Up @@ -411,13 +417,16 @@ void FitHistExtended(const char* infile, const char* outfile,std::string prefix,
TFile *file = new TFile(infile);
std::vector<TH1D*> hists;
char name[200];
if (debug) std::cout << infile << " " << file << std::endl;
if (debug) {
std::cout << infile << " " << file << std::endl;
}
if (file != 0) {
sprintf (name, "%s%s%d0", prefix.c_str(), sname.c_str(), iname);
TH1D* hist0 = (TH1D*)file->FindObjectAny(name);
bool ok = (hist0 != 0);
if (debug) std::cout << name << " Pointer " << hist0 << " " << ok
<< std::endl;
if (debug) {
std::cout << name << " Pointer " << hist0 << " " << ok << std::endl;
}
if (ok) {
TH1D* histo(0);
if (numb > 0) {
Expand All @@ -427,63 +436,78 @@ void FitHistExtended(const char* infile, const char* outfile,std::string prefix,
if (hist0->GetEntries() > 10) {
std::pair<double,double> meaner0 = fitTwoGauss(hist0, debug);
std::pair<double,double> meaner1 = GetMean(hist0,0.2,2.0);
if (debug) std::cout << "Fit " << meaner0.first << ":"
<< meaner0.second << " Mean1 "
<< hist0->GetMean() << ":"
<< hist0->GetMeanError() << " Mean2 "
<< meaner1.first << ":" << meaner1.second;
if (debug) {
std::cout << "Fit " << meaner0.first << ":" << meaner0.second
<< " Mean1 " << hist0->GetMean() << ":"
<< hist0->GetMeanError() << " Mean2 " << meaner1.first
<< ":" << meaner1.second << std::endl;
}
}
int nv1(100), nv2(0);
int jmin(numb), jmax(0);
for (int j=0; j<=numb; ++j) {
sprintf (name, "%s%s%d%d", prefix.c_str(), sname.c_str(), iname, j);
TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
if (debug) std::cout << "Get Histogram for " << name << " at " << hist1
<< std::endl;
TH1D* hist = (TH1D*)hist1->Clone();
double value(0), error(0), total(0);
if (hist->GetEntries() > 0) {
value = hist->GetMean(); error = hist->GetRMS();
for (int i=1; i<=hist->GetNbinsX(); ++i)
total += hist->GetBinContent(i);
if (debug) {
std::cout << "Get Histogram for " << name << " at " << hist1
<< std::endl;
}
if (total > 4) {
if (nv1 > j) nv1 = j;
if (nv2 < j) nv2 = j;
if (j == 0) {
sprintf (name, "%sOne", hist1->GetName());
TH1D* hist2 = (TH1D*)hist1->Clone(name);
fitOneGauss(hist2,debug);
hists.push_back(hist2);
std::pair<double,double> meaner0 = fitTwoGauss(hist,debug);
value = meaner0.first;
error = meaner0.second;
} else {
std::pair<double,double> meaner = fitOneGauss(hist,debug);
value = meaner.first; error = meaner.second;
double value(0), error(0), total(0);
if (hist1 == 0) {
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please use nullptr (and also below where appropriate)

value = 1.0;
} else {
TH1D* hist = (TH1D*)hist1->Clone();
if (hist->GetEntries() > 0) {
value = hist->GetMean(); error = hist->GetRMS();
for (int i=1; i<=hist->GetNbinsX(); ++i)
total += hist->GetBinContent(i);
}
if (j != 0) {
if (j < jmin) jmin = j;
if (j > jmax) jmax = j;
if (total > 4) {
if (nv1 > j) nv1 = j;
if (nv2 < j) nv2 = j;
if (j == 0) {
sprintf (name, "%sOne", hist1->GetName());
TH1D* hist2 = (TH1D*)hist1->Clone(name);
fitOneGauss(hist2,debug);
hists.push_back(hist2);
std::pair<double,double> meaner0 = fitTwoGauss(hist,debug);
value = meaner0.first;
error = meaner0.second;
} else {
std::pair<double,double> meaner = fitOneGauss(hist,debug);
value = meaner.first; error = meaner.second;
}
if (j != 0) {
if (j < jmin) jmin = j;
if (j > jmax) jmax = j;
}
}
hists.push_back(hist);
}
if (debug) {
std::cout << "Hist " << j << " Value " << value << " +- " << error
<< std::endl;
}
hists.push_back(hist);
if (j != 0) {
histo->SetBinContent(j, value);
histo->SetBinError(j, error);
}
}
if (histo != 0) {
if (histo->GetEntries() > 2 && fiteta) {
if (debug) std::cout << "Jmin/max " << jmin << ":" << jmax
<< ":" << histo->GetNbinsX() << std::endl;
if (debug) {
std::cout << "Jmin/max " << jmin << ":" << jmax
<< ":" << histo->GetNbinsX() << std::endl;
}
double LowEdge = histo->GetBinLowEdge(jmin);
double HighEdge= histo->GetBinLowEdge(jmax)+histo->GetBinWidth(jmax);
TFitResultPtr Fit = histo->Fit("pol0","+QRWLS","",LowEdge,HighEdge);
if (debug) std::cout << "Fit to Pol0: " << Fit->Value(0) << " +- "
<< Fit->FitResult::Error(0) << " in range "
<< nv1 << ":" << xbins[nv1] << ":" << nv2 << ":"
<< xbins[nv2] << std::endl;
if (debug) {
std::cout << "Fit to Pol0: " << Fit->Value(0) << " +- "
<< Fit->FitResult::Error(0) << " in range "
<< nv1 << ":" << xbins[nv1] << ":" << nv2 << ":"
<< xbins[nv2] << std::endl;
}
histo->GetXaxis()->SetTitle("i#eta");
histo->GetYaxis()->SetTitle("MPV(E_{HCAL}/(p-E_{ECAL}))");
histo->GetYaxis()->SetRangeUser(0.4,1.6);
Expand All @@ -497,8 +521,10 @@ void FitHistExtended(const char* infile, const char* outfile,std::string prefix,
for (int j=1; j<=3; ++j) {
sprintf (name, "%s%s%d%d", prefix.c_str(), ename.c_str(), iname, j);
TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
if (debug) std::cout << "Get Histogram for " << name << " at " << hist1
<< std::endl;
if (debug) {
std::cout << "Get Histogram for " << name << " at " << hist1
<< std::endl;
}
if (hist1!=0) {
TH1D* hist = (TH1D*)hist1->Clone();
double value(0), error(0), total(0);
Expand All @@ -516,28 +542,35 @@ void FitHistExtended(const char* infile, const char* outfile,std::string prefix,
value = meaner0.first;
error = meaner0.second;
std::pair<double,double> meaner = GetMean(hist,0.2,2.0);
if (debug) std::cout << "Fit " << value << ":" << error << ":"
<< hist->GetMeanError() << " Mean "
<< meaner.first << ":" << meaner.second
<< std::endl;
if (debug) {
std::cout << "Fit " << value << ":" << error << ":"
<< hist->GetMeanError() << " Mean "
<< meaner.first << ":" << meaner.second << std::endl;
}
}
hists.push_back(hist);
}
}
}
TFile* theFile(0);
if (append) {
if (debug) std::cout << "Open file " << outfile << " in append mode\n";
if (debug) {
std::cout << "Open file " << outfile << " in append mode" << std::endl;
}
theFile = new TFile(outfile, "UPDATE");
} else {
if (debug) std::cout << "Open file " << outfile << " in recreate mode\n";
if (debug) {
std::cout << "Open file " << outfile << " in recreate mode" <<std::endl;
}
theFile = new TFile(outfile, "RECREATE");
}

theFile->cd();
for (unsigned int i=0; i<hists.size(); ++i) {
TH1D* hnew = (TH1D*)hists[i]->Clone();
if (debug) std::cout << "Write Histogram " << hnew->GetTitle() << std::endl;
if (debug) {
std::cout << "Write Histogram " << hnew->GetTitle() << std::endl;
}
hnew->Write();
}
theFile->Close();
Expand All @@ -556,13 +589,17 @@ void FitHistRBX(const char* infile, const char* outfile,std::string prefix,
std::vector<TH1D*> hists;
sprintf (name, "%s%s%d", prefix.c_str(), lname.c_str(), iname);
TH1D *histo = new TH1D(name, "", numb, 0, numb);
if (debug) std::cout << infile << " " << file << std::endl;
if (debug) {
std::cout << infile << " " << file << std::endl;
}
if (file != 0) {
for (int j=0; j<numb; ++j) {
sprintf (name, "%s%s%d%d", prefix.c_str(), sname.c_str(), iname, j+1);
TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
if (debug) std::cout << "Get Histogram for " << name << " at " << hist1
<< std::endl;
if (debug) {
std::cout << "Get Histogram for " << name << " at " << hist1
<< std::endl;
}
TH1D* hist = (TH1D*)hist1->Clone();
double value(0), error(0), total(0);
if (hist->GetEntries() > 0) {
Expand All @@ -585,17 +622,23 @@ void FitHistRBX(const char* infile, const char* outfile,std::string prefix,

TFile* theFile(0);
if (append) {
if (debug) std::cout << "Open file " << outfile << " in append mode\n";
if (debug) {
std::cout << "Open file " << outfile << " in append mode" << std::endl;
}
theFile = new TFile(outfile, "UPDATE");
} else {
if (debug) std::cout << "Open file " << outfile << " in recreate mode\n";
if (debug) {
std::cout << "Open file " << outfile << " in recreate mode" <<std::endl;
}
theFile = new TFile(outfile, "RECREATE");
}

theFile->cd();
for (unsigned int i=0; i<hists.size(); ++i) {
TH1D* hnew = (TH1D*)hists[i]->Clone();
if (debug) std::cout << "Write Histogram " << hnew->GetTitle() << std::endl;
if (debug) {
std::cout << "Write Histogram " << hnew->GetTitle() << std::endl;
}
hnew->Write();
}
theFile->Close();
Expand Down