Skip to content

Commit 44dc539

Browse files
author
Sunanda
committed
Change FitPlot and PlotPropertiies to chck differences between EGamma and AlCaRaw
1 parent 0dd653a commit 44dc539

File tree

2 files changed

+424
-24
lines changed

2 files changed

+424
-24
lines changed

‎Calibration/HcalCalibAlgos/macros/CalibFitPlots.C‎

Lines changed: 103 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -568,6 +568,93 @@ results fitOneGauss(TH1D* hist, bool fitTwice, bool debug) {
568568
return results(value, error, width, werror);
569569
}
570570

571+
double DoubleSidedCrystalballFunction(double *x, double *par) {
572+
double alpha_l = par[0];
573+
double alpha_h = par[1];
574+
double n_l = par[2];
575+
double n_h = par[3];
576+
double mean = par[4];
577+
double sigma = par[5];
578+
double N = par[6];
579+
float t = (x[0]-mean)/sigma;
580+
double result;
581+
double fact1TLessMinosAlphaL = alpha_l/n_l;
582+
double fact2TLessMinosAlphaL = (n_l/alpha_l) - alpha_l -t;
583+
double fact1THihgerAlphaH = alpha_h/n_h;
584+
// double fact2THigherAlphaH = (n_h/alpha_h) - alpha_h + t;
585+
586+
if (-alpha_l <= t && alpha_h >= t) {
587+
result = exp(-0.5 * t * t);
588+
} else if (t < -alpha_l) {
589+
result = exp(-0.5 * alpha_l * alpha_l) * pow(fact1TLessMinosAlphaL * fact2TLessMinosAlphaL, -n_l);
590+
} else if (t > alpha_h) {
591+
result = exp(-0.5 * alpha_l * alpha_l) * pow(fact1THihgerAlphaH * fact1THihgerAlphaH, -n_h);
592+
}
593+
return N * result;
594+
}
595+
596+
results fitDoubleSidedCrystalball(TH1D* hist, bool /* fitTwice */, bool debug) {
597+
double rms;
598+
std::pair<double, double> mrms = GetMean(hist, 0.2, 2.0, rms);
599+
double mean = mrms.first;
600+
double LowEdge = ((mean - fitrangeFactor * rms) < 0.5) ? 0.5 : (mean - fitrangeFactor * rms);
601+
double diff = mean - LowEdge;
602+
double HighEdge = (diff > fitrangeFactor1 * rms) ? (mean + fitrangeFactor1 * rms) : (mean + diff);
603+
if (debug)
604+
std::cout << hist->GetName() << " Mean " << mean << " RMS " << rms << " Range " << LowEdge << ":" << HighEdge
605+
<< "\n";
606+
std::string option = (hist->GetEntries() > 100) ? "QRS" : "QRWLS";
607+
TObject* ob = gROOT->FindObject("g1");
608+
if (ob)
609+
ob->Delete();
610+
TF1* g1 = new TF1("g1", "gaus", LowEdge, HighEdge);
611+
g1->SetLineColor(kGreen);
612+
TFitResultPtr Fit1 = hist->Fit(g1, option.c_str(), "");
613+
double value = Fit1->Value(1);
614+
double error = Fit1->FitResult::Error(1);
615+
double width = Fit1->Value(2);
616+
double werror = Fit1->FitResult::Error(2);
617+
618+
LowEdge = Fit1->Value(1) - fitrangeFactor * Fit1->Value(2);
619+
if (LowEdge < 0.5)
620+
LowEdge = 0.5;
621+
diff = Fit1->Value(1) - LowEdge;
622+
HighEdge = (diff > fitrangeFactor1 * rms) ? (Fit1->Value(1) + fitrangeFactor1 * Fit1->Value(2)) : (Fit1->Value(1) + 1.5 * diff);
623+
if (HighEdge > 5.0)
624+
HighEdge = 5.0;
625+
if (debug)
626+
std::cout << " Range for second Fit " << LowEdge << ":" << HighEdge << std::endl;
627+
TObject* ob2 = gROOT->FindObject("fitDSCB");
628+
if (ob2 != nullptr)
629+
ob2->Delete();
630+
TF1 *fitDSCB = new TF1("fitDSCB", DoubleSidedCrystalballFunction, LowEdge, HighEdge, 7);
631+
fitDSCB->SetParameters(1, 2, 2, 1, value, width, hist->Integral(LowEdge, HighEdge));
632+
fitDSCB->SetParNames ("alpha_{low}","alpha_{high}","n_{low}", "n_{high}", "mean", "sigma", "Norm");
633+
TFitResultPtr Fit = hist->Fit(fitDSCB, option.c_str(), "", LowEdge, HighEdge);
634+
value = Fit->Value(4);
635+
error = Fit->FitResult::Error(4);
636+
width = Fit->Value(5);
637+
werror = Fit->FitResult::Error(5);
638+
639+
std::pair<double, double> meaner = GetMean(hist, 0.2, 2.0, rms);
640+
std::pair<double, double> rmserr = GetWidth(hist, 0.2, 2.0);
641+
if (debug) {
642+
std::cout << "FitDSCB " << value << ":" << error << ":" << hist->GetMeanError() << " Mean " << meaner.first << ":"
643+
<< meaner.second << " Width " << rmserr.first << ":" << rmserr.second;
644+
}
645+
double minvalue(0.30);
646+
if (value < minvalue || value > 2.0 || error > 0.5) {
647+
value = meaner.first;
648+
error = meaner.second;
649+
width = rmserr.first;
650+
werror = rmserr.second;
651+
}
652+
if (debug) {
653+
std::cout << " Final " << value << "+-" << error << ":" << width << "+-" << werror << std::endl;
654+
}
655+
return results(value, error, width, werror);
656+
}
657+
571658
void readCorrFactors(char* infile,
572659
double scale,
573660
std::map<int, cfactors>& cfacs,
@@ -640,7 +727,7 @@ void FitHistStandard(std::string infile,
640727
std::string xname[4] = {"i#eta", "i#eta", "d_{L1}", "# Vertex"};
641728
int numb[4] = {10, 8, 8, 5};
642729

643-
if (type == 0) {
730+
if ((type % 10) == 0) {
644731
numb[0] = 8;
645732
for (int i = 0; i < 9; ++i)
646733
xbins[i] = xbin0[i];
@@ -683,7 +770,7 @@ void FitHistStandard(std::string infile,
683770
}
684771
if (hist->GetEntries() > 4) {
685772
bool flag = (j == 0) ? true : false;
686-
results meaner = fitOneGauss(hist, flag, debug);
773+
results meaner = (((type / 10) % 10) == 0) ? fitOneGauss(hist, flag, debug) : fitDoubleSidedCrystalball(hist, flag, debug);
687774
value = meaner.mean;
688775
error = meaner.errmean;
689776
width = meaner.width;
@@ -756,11 +843,11 @@ void FitHistExtended(const char* infile,
756843
double xbins[99];
757844
double xbin[23] = {-23.0, -21.0, -19.0, -17.0, -15.0, -13.0, -11.0, -9.0, -7.0, -5.0, -3.0, 0.0,
758845
3.0, 5.0, 7.0, 9.0, 11.0, 13.0, 15.0, 17.0, 19.0, 21.0, 23.0};
759-
if (type == 2) {
846+
if ((type % 10) == 2) {
760847
numb = 22;
761848
for (int k = 0; k <= numb; ++k)
762849
xbins[k] = xbin[k];
763-
} else if (type == 1) {
850+
} else if ((type % 10) == 1) {
764851
numb = 1;
765852
xbins[0] = -25;
766853
xbins[1] = 25;
@@ -802,7 +889,7 @@ void FitHistExtended(const char* infile,
802889
}
803890
if (hist0->GetEntries() > 10) {
804891
double rms;
805-
results meaner0 = fitOneGauss(hist0, true, debug);
892+
results meaner0 = (((type / 10) % 10) == 0) ? fitOneGauss(hist0, true, debug) : fitDoubleSidedCrystalball(hist0, true, debug);
806893
std::pair<double, double> meaner1 = GetMean(hist0, 0.2, 2.0, rms);
807894
std::pair<double, double> meaner2 = GetWidth(hist0, 0.2, 2.0);
808895
if (debug) {
@@ -844,15 +931,15 @@ void FitHistExtended(const char* infile,
844931
if (j == 0) {
845932
sprintf(name, "%sOne", hist1->GetName());
846933
TH1D* hist2 = (TH1D*)hist1->Clone(name);
847-
fitOneGauss(hist2, true, debug);
934+
fitOneGauss(hist2, true, debug);
848935
hists.push_back(hist2);
849-
results meaner = fitOneGauss(hist, true, debug);
936+
results meaner = (((type / 10) % 10) == 0) ? fitOneGauss(hist, true, debug) : fitDoubleSidedCrystalball(hist, true, debug);
850937
value = meaner.mean;
851938
error = meaner.errmean;
852939
width = meaner.width;
853940
werror = meaner.errwidth;
854941
} else {
855-
results meaner = fitOneGauss(hist, true, debug);
942+
results meaner = (((type / 10) % 10) == 0) ? fitOneGauss(hist, true, debug) : fitDoubleSidedCrystalball(hist, true, debug);
856943
value = meaner.mean;
857944
error = meaner.errmean;
858945
width = meaner.width;
@@ -923,7 +1010,7 @@ void FitHistExtended(const char* infile,
9231010
if (total > 4) {
9241011
sprintf(name, "%sOne", hist1->GetName());
9251012
TH1D* hist2 = (TH1D*)hist1->Clone(name);
926-
results meanerr = fitOneGauss(hist2, true, debug);
1013+
results meanerr = (((type / 10) % 10) == 0) ? fitOneGauss(hist2, true, debug) : fitDoubleSidedCrystalball(hist2, true, debug);
9271014
value = meanerr.mean;
9281015
error = meanerr.errmean;
9291016
width = meanerr.width;
@@ -940,7 +1027,7 @@ void FitHistExtended(const char* infile,
9401027
hists.push_back(hist3);
9411028
}
9421029
// results meaner0 = fitTwoGauss(hist, debug);
943-
results meaner0 = fitOneGauss(hist, true, debug);
1030+
results meaner0 = (((type / 10) % 10) == 0) ? fitOneGauss(hist, true, debug) : fitDoubleSidedCrystalball(hist, true, debug);
9441031
value = meaner0.mean;
9451032
error = meaner0.errmean;
9461033
double rms;
@@ -984,6 +1071,7 @@ void FitHistExtended2(const char* infile,
9841071
const char* outfile,
9851072
std::string prefix,
9861073
int numb = 50,
1074+
int type = 0,
9871075
bool append = true,
9881076
int iname = 3,
9891077
bool debug = false) {
@@ -1055,15 +1143,15 @@ void FitHistExtended2(const char* infile,
10551143
if (j == 0) {
10561144
sprintf(name, "%sOne", hist1->GetName());
10571145
TH1D* hist2 = (TH1D*)hist1->Clone(name);
1058-
fitOneGauss(hist2, true, debug);
1146+
fitOneGauss(hist2, true, debug);
10591147
hists.push_back(hist2);
1060-
results meaner = fitOneGauss(hist, true, debug);
1148+
results meaner = (((type / 10) % 10) == 0) ? fitOneGauss(hist, true, debug) : fitDoubleSidedCrystalball(hist, true, debug);
10611149
value = meaner.mean;
10621150
error = meaner.errmean;
10631151
width = meaner.width;
10641152
werror = meaner.errwidth;
10651153
} else {
1066-
results meaner = fitOneGauss(hist, true, debug);
1154+
results meaner = (((type / 10) % 10) == 0) ? fitOneGauss(hist, true, debug) : fitDoubleSidedCrystalball(hist, true, debug);
10671155
value = meaner.mean;
10681156
error = meaner.errmean;
10691157
width = meaner.width;
@@ -1139,7 +1227,7 @@ void FitHistExtended2(const char* infile,
11391227
if (total > 4) {
11401228
sprintf(name, "%sOne", hist1->GetName());
11411229
TH1D* hist2 = (TH1D*)hist1->Clone(name);
1142-
results meanerr = fitOneGauss(hist2, true, debug);
1230+
results meanerr = (((type / 10) % 10) == 0) ? fitOneGauss(hist2, true, debug) : fitDoubleSidedCrystalball(hist2, true, debug);
11431231
value = meanerr.mean;
11441232
error = meanerr.errmean;
11451233
width = meanerr.width;
@@ -1149,7 +1237,7 @@ void FitHistExtended2(const char* infile,
11491237
std::cout << hist2->GetName() << " MPV " << value << " +- " << error << " Width " << width << " +- "
11501238
<< werror << " W/M " << wbyv << " +- " << wverr << std::endl;
11511239
hists.push_back(hist2);
1152-
results meaner0 = fitOneGauss(hist, true, debug);
1240+
results meaner0 = (((type / 10) % 10) == 0) ? fitOneGauss(hist, true, debug) : fitDoubleSidedCrystalball(hist, true, debug);
11531241
value = meaner0.mean;
11541242
error = meaner0.errmean;
11551243
double rms;

0 commit comments

Comments
 (0)