8 #include <TFitResult.h>
9 #include <TFitResultPtr.h>
19 #define CORRECTIONS false // when true the correction parameters will be applied
20 #define TAGM true // include TAGM when true
21 #define TAGH false // include TAGH when true
26 int bad_channels_tm = 0;
29 double fit_tm[3] = {0};
34 int bad_channels_th = 0;
37 double fit_th[3] = {0};
46 TCanvas *
c1 = (TCanvas*)gROOT->FindObject(
"c1");
48 c1 =
new TCanvas(
"c1",
"c1",0,20,600,500);
51 std::cout <<
"Getting ROOT file" << std::endl;
52 TFile *
infile =
new TFile(inputFile);
53 TFile *
outfile =
new TFile(
"profiles.root",
"recreate");
56 std::cout <<
"Creating output file directories" << std::endl;
57 outfile->mkdir(
"TAGM");
58 outfile->mkdir(
"TAGH");
63 std::ofstream fout_tm(
"Eparms-TAGM.out");
64 #endif // !CORRECTIONS
67 std::fstream parmsFile(
"Eparms-TAGM.out", ios_base::in);
68 string parm0, parm1, parm2;
73 while (parmsFile >> std::ws && std::getline(parmsFile, line))
75 std::istringstream input(line);
76 input >> parm0 >> parm1 >> parm2;
77 p0 = atof(parm0.c_str());
78 p1 = atof(parm1.c_str());
79 p2 = atof(parm2.c_str());
83 std::cout <<
"Getting original TAGM histograms" << std::endl;
90 h_psE_vs_psEl_tm[col] = (TH2F*)infile->Get(Form(
"TAGM/h_psE_vs_psEl_tm_%i",col+1));
91 psE_vs_psEl_tm_corr[col] =
new TProfile(Form(
"psE_vs_psEl_tm_corr_%i",col+1),
92 Form(
"PS E vs PS left, TAGM col %i;\
93 Energy asymmetry;PS energy (GeV)",col+1),
95 psE_vs_psEl_tm_corr[col]->SetMaximum(4.9);
96 psE_vs_psEl_tm_corr[col]->SetMinimum(2.3);
97 h_tmp = h_psE_vs_psEl_tm[col]->ProjectionY();
98 int bMax = h_tmp->GetMaximumBin();
101 for (
int i = 0; i < (
NEb_PS-bMax); ++i) {
102 if (h_tmp->GetBinContent(bMax+i) < 2) {
107 for (
int i = 0; i < bMax; ++i) {
108 if (h_tmp->GetBinContent(bMax-i) < 2) {
113 psE_vs_psEl_tm[col] = h_psE_vs_psEl_tm[col]->ProfileX(Form(
"col_%i",col+1),yMin,yMax);
115 for (
int i = 1; i <= 50; ++i) {
116 double bEntries = psE_vs_psEl_tm[col]->GetBinEntries(i);
118 psE_vs_psEl_tm[col]->SetBinEntries(i,0);
122 int nentries = psE_vs_psEl_tm[col]->GetEntries();
123 if (nentries == 0 || bad_bins > 45) {
124 for (
int i = 0; i < 3; ++i) {
125 fitResults_tm[col][i] = 0;
131 for (
int i = 1; i <= 50; ++i) {
132 double content = psE_vs_psEl_tm[col]->GetBinContent(i);
133 double asym = (i+0.5)*0.02;
135 psE_vs_psEl_tm_corr[col]->Fill(asym,content);
137 #endif // CORRECTIONS
140 TFitResultPtr ptr = psE_vs_psEl_tm[col]->Fit(
"pol2",
"sq");
141 fitResults_tm[col][0] = ptr->Parameter(0);
142 fitResults_tm[col][1] = ptr->Parameter(1);
143 fitResults_tm[col][2] = ptr->Parameter(2);
146 max_E_tm[col] = -1*fitResults_tm[col][1]*fitResults_tm[col][1];
147 max_E_tm[col] /= 4*fitResults_tm[col][2];
148 max_E_tm[col] += fitResults_tm[col][0];
150 for (
int i = 0; i < 3; ++i) {
151 fitResults_tm[col][i] /= max_E_tm[col];
152 fit_tm[i] += fitResults_tm[col][i];
154 #endif // !CORRECTIONS
159 psE_vs_psEl_tm_corr[col]->Write();
160 #endif // CORRECTIONS
164 psE_vs_psEl_tm[col]->Write();
165 #endif // !CORRECTIONS
168 for (
int i = 0; i < 3; ++i) {
169 fit_tm[i] /= (MAX_COLUMNS - bad_channels_tm);
172 fout_tm << fit_tm[0] << std::setw(15)
173 << fit_tm[1] << std::setw(15)
174 << fit_tm[2] << std::endl;
176 #endif // !CORRECTIONS
182 std::ofstream fout_th(
"Eparms-TAGH.out");
183 #endif // !CORRECTIONS
186 std::fstream parmsFile(
"Eparms-TAGH.out", ios_base::in);
187 string parm0, parm1, parm2;
192 while (parmsFile >> std::ws && std::getline(parmsFile, line))
194 std::istringstream input(line);
195 input >> parm0 >> parm1 >> parm2;
196 p0 = atof(parm0.c_str());
197 p1 = atof(parm1.c_str());
198 p2 = atof(parm2.c_str());
201 #endif // CORRECTIONS
203 std::cout <<
"Getting original TAGH histograms" << std::endl;
210 h_psE_vs_psEl_th[hodo] = (TH2F*)infile->Get(Form(
"TAGH/psE_vs_psEl_th_%i",hodo+1));
211 h_tmp = h_psE_vs_psEl_th[hodo]->ProjectionY();
212 int bMax = h_tmp->GetMaximumBin();
215 for (
int i = 0; i < (
NEb_PS-bMax); ++i) {
216 if (h_tmp->GetBinContent(bMax+i) < 2) {
221 for (
int i = 0; i < bMax; ++i) {
222 if (h_tmp->GetBinContent(bMax-i) < 2) {
227 psE_vs_psEl_th[hodo] = h_psE_vs_psEl_th[hodo]->ProfileX(Form(
"hodo_%i",hodo+1),yMin,yMax);
229 for (
int i = 1; i <= 50; ++i) {
230 double bEntries = psE_vs_psEl_th[hodo]->GetBinEntries(i);
232 psE_vs_psEl_th[hodo]->SetBinEntries(i,0);
236 int nentries = psE_vs_psEl_th[hodo]->GetEntries();
237 if (nentries == 0 || bad_bins > 45) {
238 for (
int i = 0; i < 3; ++i) {
239 fitResults_th[hodo][i] = 0;
246 psE_vs_psEl_th_corr[hodo] =
new TProfile(Form(
"psE_vs_psEl_th_corr_%i",hodo+1),
247 Form(
"PS E vs PS left, TAGH counter %i;\
248 Energy asymmetry;PS energy (GeV)",hodo+1),
250 psE_vs_psEl_th_corr[hodo]->SetMaximum(4.9);
251 psE_vs_psEl_th_corr[hodo]->SetMinimum(2.3);
252 for (
int i = 1; i <= 50; ++i) {
253 double content = psE_vs_psEl_th[hodo]->GetBinContent(i);
254 double asym = (i+0.5)*0.02;
255 content /= (p0 + p1*asym + p2*asym*asym);
256 psE_vs_psEl_th_corr[hodo]->Fill(asym,content);
258 #endif // CORRECTIONS
261 TFitResultPtr ptr = psE_vs_psEl_th[hodo]->Fit(
"pol2",
"sq");
262 fitResults_th[hodo][0] = ptr->Parameter(0);
263 fitResults_th[hodo][1] = ptr->Parameter(1);
264 fitResults_th[hodo][2] = ptr->Parameter(2);
267 max_E_th[hodo] = -1*fitResults_th[hodo][1]*fitResults_th[hodo][1];
268 max_E_th[hodo] /= 4*fitResults_th[hodo][2];
269 max_E_th[hodo] += fitResults_th[hodo][0];
271 for (
int i = 0; i < 3; ++i) {
272 fitResults_th[hodo][i] /= max_E_th[hodo];
273 fit_th[i] += fitResults_th[hodo][i];
275 #endif // !CORRECTIONS
280 psE_vs_psEl_th_corr[hodo]->Write();
281 #endif // CORRECTIONS
285 h_psE_vs_psEl_th[hodo]->Write();
286 psE_vs_psEl_th[hodo]->Write();
287 #endif // !CORRECTIONS
291 std::cout <<
"bad channels " << bad_channels_th << std::endl;
293 for (
int i = 0; i < 3; ++i) {
294 fit_th[i] /= (MAX_COUNTERS - bad_channels_th);
297 fout_th << fit_th[0] << std::setw(15)
298 << fit_th[1] << std::setw(15)
299 << fit_th[2] << std::endl;
301 #endif // !CORRECTIONS
static char index(char c)
static TH2F * h_psE_vs_psEl_th[MAX_COUNTERS]
void PSEcorr(char const *inputFile)
static TH2F * h_psE_vs_psEl_tm[MAX_COLUMNS]