BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
BmnResultsTOF700.C
Go to the documentation of this file.
1#include "../../../../bmndata/BmnEnums.h"
2// #include "../../../../decoder/BmnTof2Raw2DigitNew.h"
3#include "../../../run/bmnloadlibs.C"
4#include "TCanvas.h"
5#include "TChain.h"
6#include "TClonesArray.h"
7#include "TFile.h"
8#include "TGeoManager.h"
9#include "TH1F.h"
10#include "TH2F.h"
11#include "TMath.h"
12#include "TProfile.h"
13#include "TStyle.h"
14#include "TVector3.h"
15
16#include <vector>
17
18#define NUMBER_CHAMBERS 59
19#define IN_WIDTH_LIMITS 1
20#define WIDTH_MAX 5000
21
22//-------------- Time difference between chambers -------------------------
23
24// No stack chambers
25// #define DIFF_CHAMB_1 -1
26// #define DIFF_CHAMB_2 -1
27
28// Stack chambers for run periods 6-7
29#define DIFF_CHAMB_1 9
30#define DIFF_CHAMB_2 18
31
32// all strips
33#define STRIP_1 0
34#define STRIP_2 31
35#define STRIP_1_2 0
36#define STRIP_2_2 31
37
38// strip 21
39// #define STRIP_1 21
40// #define STRIP_2 21
41// #define STRIP_1_2 21
42// #define STRIP_2_2 21
43
44//-------------------------------------------------------------------------
45
46int champos[NUMBER_CHAMBERS] = {0};
47int wmi[NUMBER_CHAMBERS] = {2800};
48int wma[NUMBER_CHAMBERS] = {4000};
49#define NDX 8
50#define NDY 8
51
52float idchambers[59] = {27.1, 28.1, 3.1, 1.1, 29.1, 4.1, 33.1, 30.1, 5.1, 19.3, 31.1, 6.1,
53 2.1, 32.1, 15.2, 16.2, 17.2, 18.2, 19.2, 20.2, 7.1, 115.2, 113.1, 117.1,
54 35.1, 9.1, 37.1, 11.1, 39.1, 13.1, 34.1, 8.1, 36.1, 10.1, 38.1, 12.1,
55 21.2, 23.2, 25.2, 22.2, 24.2, 26.2, 107.2, 108.2, 109.2, 110.2, 111.2, 112.2,
56 114.1, 116.2, 118.1, 14.1, 40.1, 119.2, 120.2, 121.2, 122.2, 123.2, 124.2};
57
58int chamnum_left[30];
59
60int chamnum_right[NUMBER_CHAMBERS - 30];
61
62int champos_left[30] = {25, 26, 11, 9, 27, 12, 31, 28, 13, 23, 29, 14, 10, 30, 17,
63 18, 19, 20, 21, 22, 15, 16, 8, 24, 1, 3, 5, 33, 35, 37};
64
65int champos_right[NUMBER_CHAMBERS - 30] = {25, 9, 27, 11, 29, 13, 24, 8, 26, 10, 28, 12, 17, 19, 21,
66 18, 20, 22, 1, 3, 5, 15, 23, 31, 14, 30, 33, 35, 37};
67
68using namespace TMath;
69
70void draw_vs_cham(char*, TH1F**);
71void draw_vs_cham_prof(char*, TProfile**);
72double FitIn(TH1F* h, Double_t xmin, Double_t xmax);
73
74class BmnTof2Raw2DigitNew;
75
76BmnTof2Raw2DigitNew* TOF2 = NULL;
77
78void BmnResultsTOF700(const char* fname = "./5070/bmn_run5070_digi.root")
79{
80
81#if ROOT_VERSION_CODE < ROOT_VERSION(5, 99, 99)
82 gROOT->LoadMacro("$VMCWORKDIR/macro/run/bmnloadlibs.C");
83#endif
84 bmnloadlibs(); // load BmnRoot libraries
85
86 const char* bname = gSystem->BaseName(fname);
87 char runfile[128];
88 strcpy(runfile, bname);
89 char* cdigi = strstr(runfile, "digi");
90 *cdigi = '\0';
91 strcat(runfile, "raw.root");
92 printf("Decoded raw ROOT file name is %s\n", runfile);
93
94 TOF2 = new BmnTof2Raw2DigitNew("TOF700_map_period_7.txt", fname, 0, 0, "TOF700_geometry_run7.txt");
95
96 if (NUMBER_CHAMBERS > 24) {
97 for (int i = 0; i < NUMBER_CHAMBERS; i++) {
98 champos[i] = i;
99 wmi[i] = 3000;
100 wma[i] = 8000;
101 }
102 }
103 for (int i = 0; i < 30; i++) {
104 if (i < 24)
105 chamnum_left[i] = i;
106 else if (i < 27) // i = 24,25,26 => 42,43,44
107 chamnum_left[i] = i + 18;
108 else // i = 27,28,29 => 53,54,55
109 chamnum_left[i] = i + 26;
110 }
111 for (int i = 24; i < 53; i++) // 24,25,26,...,50,51,52 (29 chambers)
112 {
113 if (i < 42)
114 chamnum_right[i - 24] = i;
115 else if (i < 50) // i = 42-49 => 45-52
116 chamnum_right[i - 24] = i + 3;
117 else // i = 50,51,52 => 56,57,58
118 chamnum_right[i - 24] = i + 6;
119 }
120
121 TClonesArray* tof700Digits = NULL;
122 TChain* bmnTree = new TChain("cbmsim");
123 if (bmnTree->Add(fname) == 0) {
124 printf("Can't find BMN digits tree in file %s!\n", fname);
125 return;
126 } else {
127 bmnTree->SetBranchAddress("TOF700", &tof700Digits);
128 }
129
130 char name[128], title[128];
131 Int_t startEvent = 0;
132 Int_t nEvents = bmnTree->GetEntries();
133
134 printf("\nNumber of events in tree is %d\n", nEvents);
135
136 gStyle->SetOptFit(111);
137
138 TH1F* htime[NUMBER_CHAMBERS] = {0};
139 TH1F* hwidth[NUMBER_CHAMBERS] = {0};
140 TH1F* hlrdiff[NUMBER_CHAMBERS] = {0};
141 TH1F* hstrips[NUMBER_CHAMBERS] = {0};
142 TH1F* hstripsmax[NUMBER_CHAMBERS] = {0};
143 TH1F* htimemax[NUMBER_CHAMBERS] = {0};
144 TProfile* timeprof[NUMBER_CHAMBERS] = {0};
145 TH1F* hdiff = 0;
146 TH2F* htvst = 0;
147 TH2F *htvss1 = 0, *htvss2 = 0, *hsvss = 0;
148 TH1F *htimemax1 = 0, *htimemax2 = 0;
149 TH1F* hxhit = 0;
150
151 TFile* f = NULL;
152 if (strlen(fname)) {
153 char fname_root[128], *delim = 0;
154 strcpy(fname_root, gSystem->BaseName(fname));
155 if ((delim = strrchr(fname_root, (int)'.')))
156 *delim = '\0';
157 strcat(fname_root, "_results.root");
158 f = new TFile(fname_root, "RECREATE", "Results of BmnTOF700");
159 }
160
161 for (int i = 0; i < NUMBER_CHAMBERS; i++) {
162 sprintf(name, "Time_chamber_%.1f", idchambers[i]);
163 sprintf(title, "Time, chamber %.1f", idchambers[i]);
164 htime[i] = new TH1F(name, title, 4000, -50, 50);
165 sprintf(name, "Width_chamber_%.1f", idchambers[i]);
166 sprintf(title, "Width, chamber %.1f", idchambers[i]);
167 hwidth[i] = new TH1F(name, title, 2500, 0, 5000);
168 sprintf(name, "LR-Diff_chamber_%.1f", idchambers[i]);
169 sprintf(title, "LR-Diff, chamber %.1f", idchambers[i]);
170 hlrdiff[i] = new TH1F(name, title, 200, -5, 5);
171 sprintf(name, "Strips_rate_chamber_%.1f", idchambers[i]);
172 sprintf(title, "Strips rate, chamber %.1f", idchambers[i]);
173 hstrips[i] = new TH1F(name, title, 32, 0, 32);
174 sprintf(name, "Time_chamber_%.1f_max_strip", idchambers[i]);
175 sprintf(title, "Time, chamber %.1f, max strip", idchambers[i]);
176 htimemax[i] = new TH1F(name, title, 4000, -50, 50);
177 sprintf(name, "Strips_rate_chamber_%.1f_max_strip", idchambers[i]);
178 sprintf(title, "Strips rate, chamber %.1f, max strip", idchambers[i]);
179 hstripsmax[i] = new TH1F(name, title, 32, 0, 32);
180 sprintf(name, "Time_profile_chamber_%.1f", idchambers[i]);
181 sprintf(title, "Time profile, chamber %.1f", idchambers[i]);
182 timeprof[i] = new TProfile(name, title, 1000, 0, 5000, -20, +20, "e");
183 }
184
185 sprintf(name, "X_coordinate");
186 sprintf(title, "X coordinate of TOF700 hit");
187 hxhit = new TH1F(name, title, 300, -150, +150);
188
189 sprintf(name, "Time_diff_chambers_%.1f_%.1f", idchambers[DIFF_CHAMB_1], idchambers[DIFF_CHAMB_2]);
190 sprintf(title, "Time difference, chambers %.1f and %.1f", idchambers[DIFF_CHAMB_1], idchambers[DIFF_CHAMB_2]);
191 hdiff = new TH1F(name, title, 10000, -200, 200);
192
193 sprintf(name, "Time_vs_time_chambers_%.1f_%.1f", idchambers[DIFF_CHAMB_1], idchambers[DIFF_CHAMB_2]);
194 sprintf(title, "Time vs Time, chambers %.1f and %.1f", idchambers[DIFF_CHAMB_1], idchambers[DIFF_CHAMB_2]);
195 htvst = new TH2F(name, title, 1000, -20, 20, 1000, -20, 20);
196
197 sprintf(name, "Time_chamber_%.1f_max_strip_both", idchambers[DIFF_CHAMB_1]);
198 sprintf(title, "Time, chamber %.1f, max strip, both chamber", idchambers[DIFF_CHAMB_1]);
199 htimemax1 = new TH1F(name, title, 2000, -50, 50);
200
201 sprintf(name, "Time_chamber_%.1f_max_strip_both", idchambers[DIFF_CHAMB_2]);
202 sprintf(title, "Time, chamber %.1f, max strip, both chamber", idchambers[DIFF_CHAMB_2]);
203 htimemax2 = new TH1F(name, title, 2000, -50, 50);
204
205 sprintf(name, "Time_vs_strip_chamber_%.1f", idchambers[DIFF_CHAMB_1]);
206 sprintf(title, "Time vs strip, chamber %.1f", idchambers[DIFF_CHAMB_1]);
207 htvss1 = new TH2F(name, title, 32, 0, 32, 2000, -50, 50);
208
209 sprintf(name, "Time_vs_strip_chamber_%.1f", idchambers[DIFF_CHAMB_2]);
210 sprintf(title, "Time vs strip, chamber %.1f", idchambers[DIFF_CHAMB_2]);
211 htvss2 = new TH2F(name, title, 32, 0, 32, 2000, -50, 50);
212
213 sprintf(name, "Strip_vs_strip_chambers_%.1f_%.1f _", idchambers[DIFF_CHAMB_1], idchambers[DIFF_CHAMB_2]);
214 sprintf(title, "Strip vs strip, chambers %.1f and %.1f", idchambers[DIFF_CHAMB_1], idchambers[DIFF_CHAMB_2]);
215 hsvss = new TH2F(name, title, 32, 0, 32, 32, 0, 32);
216
217 Int_t nchamb1 = 0, nchamb2 = 0, nchamb12 = 0, nchamb21 = 0, nchamb12n = 0, nchamb12nn = 0;
218
219 for (Int_t iEv = startEvent; iEv < startEvent + nEvents; iEv++) {
220 bmnTree->GetEntry(iEv);
221
222 if (iEv % 10000 == 0) {
223 cout << "Event: " << iEv << "/" << startEvent + nEvents << endl;
224 }
225 float wmax[NUMBER_CHAMBERS] = {0.};
226 float tmax[NUMBER_CHAMBERS] = {-10000000.};
227 int smax[NUMBER_CHAMBERS] = {-1};
228 int smax1[NUMBER_CHAMBERS] = {-1};
229
230 for (int i = 0; i < NUMBER_CHAMBERS; i++) {
231 wmax[i] = 0.;
232 tmax[i] = -1000000.;
233 smax[i] = -1;
234 smax1[i] = -1;
235 }
236
237 for (Int_t iDig = 0; iDig < tof700Digits->GetEntriesFast(); ++iDig) {
238 BmnTof2Digit* digit = (BmnTof2Digit*)tof700Digits->At(iDig);
239 if (digit == NULL)
240 continue;
241 Short_t plane = digit->GetPlane();
242 Short_t strip = digit->GetStrip();
243 Float_t time = digit->GetTime();
244 Float_t width = digit->GetAmplitude();
245 Float_t lrdiff = digit->GetDiff();
246 // if (iEv % 10 == 0) cout << "plane = " << plane << " strip = " << strip << " time = " << time <<
247 // " width = " << width << endl; if (iEv % 10 == 0) cout << "wmi = " << wmi[plane] << " wma = "
248 // << wma[plane] << endl;
249 if (plane >= NUMBER_CHAMBERS)
250 continue;
251 // if (lrdiff < -0.5 || lrdiff > 0.5) continue;
252 hwidth[plane]->Fill(width);
253 smax1[plane] = strip;
254 if (width > wmax[plane]) {
255 wmax[plane] = width;
256 tmax[plane] = time;
257 smax[plane] = strip;
258 }
259 if (width > WIDTH_MAX)
260 continue;
261 if (IN_WIDTH_LIMITS)
262 if (width < wmi[plane] || width > wma[plane])
263 continue;
264 htime[plane]->Fill(time);
265 hstrips[plane]->Fill(strip);
266 hlrdiff[plane]->Fill(lrdiff);
267 hxhit->Fill(TOF2->get_hit_x(plane, strip, lrdiff));
268 } // tof700Digits loop
269
270 if (DIFF_CHAMB_1 >= 0)
271 if ((smax1[DIFF_CHAMB_1] >= STRIP_1) && (smax1[DIFF_CHAMB_1] >= STRIP_2)) {
272 nchamb1++;
273 if ((smax1[DIFF_CHAMB_2] >= 0))
274 nchamb12++;
275 if ((smax1[DIFF_CHAMB_2 - 1] >= 0) || (smax1[DIFF_CHAMB_2] >= 0) || (smax1[DIFF_CHAMB_2 + 1] >= 0))
276 nchamb12n++;
277 if ((smax1[DIFF_CHAMB_2 - 1] >= 0) || (smax1[DIFF_CHAMB_2] >= 0) || (smax1[DIFF_CHAMB_2 + 1] >= 0)
278 || (smax1[8] >= 0) || (smax1[10] >= 0))
279 nchamb12nn++;
280 }
281 if (DIFF_CHAMB_1 >= 0)
282 if ((smax1[DIFF_CHAMB_2] >= STRIP_1_2) && (smax1[DIFF_CHAMB_2] >= STRIP_2_2)) {
283 nchamb2++;
284 if ((smax1[DIFF_CHAMB_1] >= 0))
285 nchamb21++;
286 }
287 for (int i = 0; i < NUMBER_CHAMBERS; i++) {
288 if (smax[i] > -1) {
289 timeprof[i]->Fill(wmax[i], tmax[i]);
290 if (wmax[i] > WIDTH_MAX)
291 continue;
292 if (IN_WIDTH_LIMITS)
293 if (wmax[i] < wmi[i] || wmax[i] > wma[i])
294 continue;
295 htimemax[i]->Fill(tmax[i]);
296 hstripsmax[i]->Fill(smax[i]);
297 }
298 }
299 if (DIFF_CHAMB_1 < 0)
300 continue;
301 if (IN_WIDTH_LIMITS)
302 if (wmax[DIFF_CHAMB_1] < wmi[DIFF_CHAMB_1] || wmax[DIFF_CHAMB_1] > wma[DIFF_CHAMB_1])
303 continue;
304 if (IN_WIDTH_LIMITS)
305 if (wmax[DIFF_CHAMB_2] < wmi[DIFF_CHAMB_2] || wmax[DIFF_CHAMB_2] > wma[DIFF_CHAMB_2])
306 continue;
307 if (smax[DIFF_CHAMB_1] >= 0)
308 htvss1->Fill(smax[DIFF_CHAMB_1], tmax[DIFF_CHAMB_1]);
309 if (smax[DIFF_CHAMB_2] >= 0)
310 htvss2->Fill(smax[DIFF_CHAMB_2], tmax[DIFF_CHAMB_2]);
311 if (smax[DIFF_CHAMB_1] >= 0 && smax[DIFF_CHAMB_2] >= 0)
312 hsvss->Fill(smax[DIFF_CHAMB_2], smax[DIFF_CHAMB_1]);
313 if ((smax[DIFF_CHAMB_1] >= STRIP_1) && (smax[DIFF_CHAMB_2] >= STRIP_1_2) && (smax[DIFF_CHAMB_1] <= STRIP_2)
314 && (smax[DIFF_CHAMB_2] <= STRIP_2_2))
315 {
316 htvst->Fill(tmax[DIFF_CHAMB_2], tmax[DIFF_CHAMB_1]);
317 htimemax1->Fill(tmax[DIFF_CHAMB_1]);
318 htimemax2->Fill(tmax[DIFF_CHAMB_2]);
319 // if ((tmax[DIFF_CHAMB_1] > 0.3) && (tmax[DIFF_CHAMB_2] < 0.3) )
320 hdiff->Fill(tmax[DIFF_CHAMB_1] - tmax[DIFF_CHAMB_2]);
321 }
322 } // event loop
323
324 if (nchamb1 > 0)
325 printf("Efficiency (normalize to chamber %3d) = %.1f %%\n", DIFF_CHAMB_1, (float)nchamb12 / nchamb1 * 100.);
326 if (nchamb1 > 0)
327 printf("Efficiency (normalize to chamber %3d) = %.1f %% (3 chambers)\n", DIFF_CHAMB_1,
328 (float)nchamb12n / nchamb1 * 100.);
329 if (nchamb1 > 0)
330 printf("Efficiency (normalize to chamber %3d) = %.1f %% (5 chambers)\n", DIFF_CHAMB_1,
331 (float)nchamb12nn / nchamb1 * 100.);
332 if (nchamb2 > 0)
333 printf("Efficiency (normalize to chamber %3d) = %.1f %%\n", DIFF_CHAMB_2, (float)nchamb21 / nchamb1 * 100.);
334 if (DIFF_CHAMB_1 >= 0)
335 FitIn(hdiff, -0.2, +0.2);
336 for (int i = 0; i < NUMBER_CHAMBERS; i++)
337 FitIn(htime[i], -0.2, +0.2);
338 for (int i = 0; i < NUMBER_CHAMBERS; i++)
339 FitIn(htimemax[i], -0.2, +0.2);
340
341 // FitIn2Sigma(hdiff);
342 // FitIn2Sigma(htimemax1);
343 // FitIn2Sigma(htimemax2);
344 // for (int i=0; i<NUMBER_CHAMBERS; i++) FitIn2Sigma(htime[i]);
345 // for (int i=0; i<NUMBER_CHAMBERS; i++) FitIn2Sigma(htimemax[i]);
346
347 draw_vs_cham((char*)"Time", htime);
348 draw_vs_cham((char*)"TimeMax", htimemax);
349 draw_vs_cham((char*)"Rate", hstrips);
350 draw_vs_cham((char*)"RateMax", hstripsmax);
351 draw_vs_cham((char*)"Width", hwidth);
352 draw_vs_cham((char*)"Diff", hlrdiff);
353 draw_vs_cham_prof((char*)"Profile", timeprof);
354
355 TCanvas* cx = new TCanvas("cx", "X-coordinate", 900, 900);
356 cx->Divide(1, 1);
357 cx->cd();
358 cx->cd(1);
359 hxhit->Draw();
360
361 /*
362 TCanvas *ct = new TCanvas("ct", "RPC time - T0 time", 900, 900);
363 ct->Divide(NDX,NDY);
364 ct->cd();
365 for (int i=0; i<NUMBER_CHAMBERS; i++)
366 {
367 ct->cd(champos[i]+1);
368 htime[i]->Draw();
369 gPad->AddExec("exselt","select_hist()");
370 }
371
372 TCanvas *c1m = new TCanvas("c1m", "RPC time max strip - T0 time", 900, 900);
373 c1m->Divide(NDX,NDY);
374 c1m->cd();
375 for (int i=0; i<NUMBER_CHAMBERS; i++)
376 {
377 c1m->cd(champos[i]+1);
378 htimemax[i]->Draw();
379 gPad->AddExec("exselt","select_hist()");
380 }
381
382
383 TCanvas *c2 = new TCanvas("c2", "RPC strip rate", 900, 900);
384 c2->Divide(NDX,NDY);
385 c2->cd();
386 for (int i=0; i<NUMBER_CHAMBERS; i++)
387 {
388 c2->cd(champos[i]+1);
389 hstrips[i]->Draw();
390 gPad->AddExec("exselt","select_hist()");
391 }
392
393 TCanvas *c2m = new TCanvas("c2m", "RPC maximal strip rate", 900, 900);
394 c2m->Divide(NDX,NDY);
395 c2m->cd();
396 for (int i=0; i<NUMBER_CHAMBERS; i++)
397 {
398 c2m->cd(champos[i]+1);
399 hstripsmax[i]->Draw();
400 gPad->AddExec("exselt","select_hist()");
401 }
402
403 TCanvas *cw = new TCanvas("cw", "RPC pulse width", 900, 900);
404 cw->Divide(NDX,NDY);
405 cw->cd();
406 for (int i=0; i<NUMBER_CHAMBERS; i++)
407 {
408 cw->cd(champos[i]+1);
409 hwidth[i]->Draw();
410 gPad->AddExec("exselt","select_hist()");
411 }
412
413 TCanvas *clr = new TCanvas("clr", "LR Time Difference", 900, 900);
414 clr->Divide(NDX,NDY);
415 clr->cd();
416 for (int i=0; i<NUMBER_CHAMBERS; i++)
417 {
418 clr->cd(champos[i]+1);
419 hlrdiff[i]->Draw();
420 gPad->AddExec("exselt","select_hist()");
421 }
422
423 TCanvas *ctim = new TCanvas("ctim", "Time Profiles", 900, 900);
424 ctim->Divide(NDX,NDY);
425 ctim->cd();
426 for (int i=0; i<NUMBER_CHAMBERS; i++)
427 {
428 ctim->cd(champos[i]+1);
429 timeprof[i]->Draw();
430 gPad->AddExec("exselt","select_hist()");
431 }
432 */
433
434 if (DIFF_CHAMB_1 < 0)
435 return;
436
437 TCanvas* cdif = new TCanvas("cdif", "Time Difference", 900, 900);
438 cdif->Divide(2, 2);
439 cdif->cd();
440 cdif->cd(1);
441 hdiff->Draw();
442 gPad->AddExec("exselt", "select_hist()");
443 cdif->cd(2);
444 htimemax1->Draw();
445 gPad->AddExec("exselt", "select_hist()");
446 cdif->cd(3);
447 htvst->Draw();
448 gPad->AddExec("exselt", "select_hist()");
449 cdif->cd(4);
450 htimemax2->Draw();
451 gPad->AddExec("exselt", "select_hist()");
452
453 TCanvas* cstr = new TCanvas("cstr", "Time vs Strip", 900, 900);
454 cstr->Divide(1, 2);
455 cstr->cd();
456 cstr->cd(1);
457 htvss1->Draw();
458 gPad->AddExec("exselt", "select_hist()");
459 cstr->cd(2);
460 htvss2->Draw();
461 gPad->AddExec("exselt", "select_hist()");
462
463 TCanvas* cstrstr = new TCanvas("cstrstr", "Strip vs strip", 900, 900);
464 cstrstr->Divide(1, 1);
465 cstrstr->cd();
466 cstrstr->cd(1);
467 hsvss->Draw();
468 gPad->AddExec("exselt", "select_hist()");
469
470 if (f)
471 f->Write();
472}
473
474double FitIn2Sigma(TH1F* h)
475{
476 double mean, sigm, sigmw = 0., lev1 = 1.75, lev2 = 3.;
477 double total = 0, peak = 0;
478 int mib, mab, nbi;
479
480 if ((total = h->GetEntries()) < 5.) {
481 printf("Number of histogram '%s' entries is too few!\n", h->GetTitle());
482 return 0.;
483 }
484 int nonzero = 0;
485 for (int i = 1; i < h->GetNbinsX(); i++)
486 if (h->GetBinContent(i) > 0)
487 nonzero++;
488 if (nonzero < 3) {
489 if (nonzero < 2) {
490 printf("Number of non-empty bins is too few in histogram '%s'!\n", h->GetTitle());
491 return 0.;
492 }
493 lev1 = 5.;
494 mean = h->GetMean();
495 sigm = h->GetRMS();
496 sigmw = sigm;
497 if (sigm < h->GetBinWidth(1))
498 sigmw = h->GetBinWidth(1);
499 TF1* fg = new TF1("fg", "gaus", mean - lev1 * sigmw, mean + lev1 * sigmw);
500 fg->SetParameter(0, h->GetMaximum());
501 fg->SetParameter(1, mean);
502 fg->FixParameter(1, mean);
503 fg->SetParameter(2, sigm);
504 fg->FixParameter(2, sigm);
505 h->Fit(fg, "Q0IB", "", mean - lev1 * sigmw, mean + lev1 * sigmw);
506 fg = h->GetFunction("fg");
507 if (fg)
508 fg->ResetBit(TF1::kNotDraw);
509 h->SetAxisRange(mean - lev1 * sigmw, mean + lev1 * sigmw);
510 fg->ResetBit(TF1::kNotDraw);
511 if (fg->GetParameter(0) > 1.)
512 return 100.;
513 }
514 if (nonzero == 3) {
515 lev1 = 5.;
516 mean = h->GetMean();
517 sigm = h->GetRMS();
518 sigmw = sigm;
519 if (sigm < h->GetBinWidth(1))
520 sigmw = h->GetBinWidth(1);
521 TF1* fg = new TF1("fg", "gaus", mean - lev1 * sigmw, mean + lev1 * sigmw);
522 fg->SetParameter(0, h->GetMaximum());
523 fg->SetParameter(1, mean);
524 fg->FixParameter(1, mean);
525 fg->SetParameter(2, sigm);
526 h->Fit(fg, "Q0IB", "", mean - lev1 * sigmw, mean + lev1 * sigmw);
527 fg = h->GetFunction("fg");
528 if (fg)
529 fg->ResetBit(TF1::kNotDraw);
530 h->SetAxisRange(mean - lev1 * sigmw, mean + lev1 * sigmw);
531 fg->ResetBit(TF1::kNotDraw);
532 if (fg->GetParameter(0) > 1.)
533 return 100.;
534 }
535 // mean = h->GetMean();
536 mean = h->GetBinCenter(h->GetMaximumBin());
537 // if (fabs(mean) > 1.) mean = 0.;
538 // mean = 0.;
539 sigm = h->GetRMS();
540 // if (sigm > 1.) sigm = 1.;
541 // sigm = 0.5;
542 mib = h->FindBin(mean - lev1 * sigm);
543 mab = h->FindBin(mean + lev1 * sigm);
544 if ((mab - mib + 1) < 4) {
545 lev1 = 2.;
546 mib = h->FindBin(mean - lev1 * sigm);
547 mab = h->FindBin(mean + lev1 * sigm);
548 if ((mab - mib + 1) < 4) {
549 lev1 = 2.25;
550 }
551 }
552 double min0 = mean - lev1 * sigm;
553 double max0 = mean + lev1 * sigm;
554 h->Fit("gaus", "Q0I", "", mean - lev1 * sigm, mean + lev1 * sigm);
555 TF1* ff = h->GetFunction("gaus");
556 if (ff)
557 mean = ff->GetParameter(1);
558 if (ff)
559 sigm = ff->GetParameter(2);
560 if (mean < min0 || mean > max0) {
561 h->Fit("gaus", "Q0I", "", mean - lev1 * sigm, mean + lev1 * sigm);
562 ff = h->GetFunction("gaus");
563 mean = ff->GetParameter(1);
564 sigm = ff->GetParameter(2);
565 }
566 h->Fit("gaus", "Q0I", "", mean - lev1 * sigm, mean + lev1 * sigm);
567 ff = h->GetFunction("gaus");
568 if (ff)
569 mean = ff->GetParameter(1);
570 if (ff)
571 sigm = ff->GetParameter(2);
572 h->Fit("gaus", "Q0I", "", mean - lev1 * sigm, mean + lev1 * sigm);
573 ff = h->GetFunction("gaus");
574 if (ff)
575 mean = ff->GetParameter(1);
576 if (ff)
577 sigm = ff->GetParameter(2);
578 if (ff)
579 ff->ResetBit(TF1::kNotDraw);
580 mib = h->FindBin(mean - lev2 * sigm);
581 mab = h->FindBin(mean + lev2 * sigm);
582 nbi = h->GetNbinsX();
583 peak = h->Integral(mib < 1 ? 1 : mib, mab > nbi ? nbi : mab);
584 h->SetAxisRange(mean - 10. * sigm, mean + 10. * sigm);
585 return total > 0. ? peak / total * 100. : 0.;
586}
587
588double FitIn(TH1F* h, Double_t xmin, Double_t xmax)
589{
590 double mean, sigm;
591 double total = 0, peak = 0;
592 int mib, mab, nbi;
593
594 if ((total = h->GetEntries()) < 5.) {
595 printf("Number of histogram '%s' entries is too few (%d)!\n", h->GetTitle(), (int)h->GetEntries());
596 return 0.;
597 }
598 int nonzero = 0;
599 for (int i = 1; i < h->GetNbinsX(); i++)
600 if (h->GetBinContent(i) > 0)
601 nonzero++;
602 if (nonzero <= 3) {
603 printf("Number of histogram '%s' non-zero bins is too few (%d)!\n", h->GetTitle(), nonzero);
604 return 0.;
605 }
606 TF1* ff = new TF1("ff", "gaus", xmin, xmax);
607 h->Fit("ff", "QIRM0");
608 // h->Fit("gaus", "Q0I");
609 // ff = h->GetFunction("ff");
610 mean = ff->GetParameter(1);
611 sigm = ff->GetParameter(2);
612 ff = h->GetFunction("ff");
613 if (ff)
614 ff->ResetBit(TF1::kNotDraw);
615 mib = h->FindBin(mean - sigm);
616 mab = h->FindBin(mean + sigm);
617 nbi = h->GetNbinsX();
618 peak = h->Integral(mib < 1 ? 1 : mib, mab > nbi ? nbi : mab);
619 h->SetAxisRange(mean - 5. * sigm, mean + 5. * sigm);
620 return total > 0. ? peak / total * 100. : 0.;
621}
622
623void select_hist()
624{
625 // example of macro called when a mouse event occurs in a pad.
626 // Example:
627 // Root > TFile f("hsimple.root");
628 // Root > hpxpy.Draw();
629 // Root > c1.AddExec("exsel","select_hist()");
630 // When moving the mouse in the canvas, a second canvas shows the
631 // projection along X of the bin corresponding to the Y position
632 // of the mouse. The resulting histogram is fitted with a gaussian.
633 // A "dynamic" line shows the current bin position in Y.
634 // This more elaborated example can be used as a starting point
635 // to develop more powerful interactive applications exploiting CINT
636 // as a development engine.
637 // Author: Rene Brun
638
639 int event = gPad->GetEvent();
640 if (event != 11)
641 return;
642
643 TObject* select = gPad->GetSelected();
644 if (!select)
645 return;
646 if (!(select->InheritsFrom("TH1")))
647 return;
648
649 // create or set the new canvas c2
650 TCanvas* csingle = 0;
651 if ((csingle = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("csingle")))
652 ;
653 else
654 csingle = new TCanvas("csingle");
655 csingle->cd();
656
657 // draw slice corresponding to mouse position
658 if (!(select->InheritsFrom("TH1F"))) {
659 TH1F* h = (TH1F*)select;
660 h->Draw();
661 } else if (!(select->InheritsFrom("TH2F"))) {
662 TH2F* h = (TH2F*)select;
663 h->Draw();
664 } else if (!(select->InheritsFrom("TProfile"))) {
665 TProfile* h = (TProfile*)select;
666 h->Draw();
667 } else if (!(select->InheritsFrom("TH1D"))) {
668 TH1D* h = (TH1D*)select;
669 h->Draw();
670 } else if (!(select->InheritsFrom("TH2D"))) {
671 TH2D* h = (TH2D*)select;
672 h->Draw();
673 } else if (!(select->InheritsFrom("TProfile2D"))) {
674 TProfile2D* h = (TProfile2D*)select;
675 h->Draw();
676 } else
677 return;
678 csingle->Update();
679}
680
681void draw_vs_cham(char* name, TH1F** h)
682{
683 char namet[128], namen[128];
684 strcpy(namet, name);
685 strcpy(namen, name);
686 TCanvas* callr1 = new TCanvas(strcat(namen, "l"), strcat(namet, " left side"), 900, 700);
687 callr1->cd();
688 callr1->Divide(8, 5);
689 for (int i = 0; i < 30; i++) {
690 callr1->cd(champos_left[i] + 1);
691 h[chamnum_left[i]]->Draw();
692 gPad->AddExec("exselt", "select_hist()");
693 }
694 strcpy(namet, name);
695 strcpy(namen, name);
696 TCanvas* callr2 = new TCanvas(strcat(namen, "r"), strcat(namet, " right side"), 900, 700);
697 callr2->cd();
698 callr2->Divide(8, 5);
699 // chambers 24-41, 45-52, 56-58
700 for (int i = 24; i < 53; i++) // 24,25,26,...,50,51,52 (29 chambers)
701 {
702 callr2->cd(champos_right[i - 24] + 1);
703 h[chamnum_right[i - 24]]->Draw();
704 gPad->AddExec("exselt", "select_hist()");
705 }
706 return;
707}
708
709void draw_vs_cham_prof(char* name, TProfile** h)
710{
711 char namet[128], namen[128];
712 strcpy(namet, name);
713 strcpy(namen, name);
714 TCanvas* callr1 = new TCanvas(strcat(namen, "l"), strcat(namet, " left side"), 900, 700);
715 callr1->cd();
716 callr1->Divide(8, 5);
717 for (int i = 0; i < 30; i++) {
718 callr1->cd(champos_left[i] + 1);
719 h[chamnum_left[i]]->Draw();
720 gPad->AddExec("exselt", "select_hist()");
721 }
722 strcpy(namet, name);
723 strcpy(namen, name);
724 TCanvas* callr2 = new TCanvas(strcat(namen, "r"), strcat(namet, " right side"), 900, 700);
725 callr2->cd();
726 callr2->Divide(8, 5);
727 // chambers 24-41, 45-52, 56-58
728 for (int i = 24; i < 53; i++) // 24,25,26,...,50,51,52 (29 chambers)
729 {
730 callr2->cd(champos_right[i - 24] + 1);
731 h[chamnum_right[i - 24]]->Draw();
732 gPad->AddExec("exselt", "select_hist()");
733 }
734 return;
735}
i
int i
Definition P4_F32vec4.h:22
Int_t
int Int_t
Float_t
float Float_t
Short_t
short Short_t
Double_t
double Double_t
gROOT
#define gROOT
gStyle
R__EXTERN TStyle * gStyle
gSystem
R__EXTERN TSystem * gSystem
gPad
#define gPad
BmnTof1Digit::GetPlane
Short_t GetPlane() const
Definition BmnTof1Digit.h:23
BmnTof1Digit::GetStrip
Short_t GetStrip() const
Definition BmnTof1Digit.h:24
BmnTof1Digit::GetTime
Float_t GetTime() const
Definition BmnTof1Digit.h:27
BmnTof1Digit::GetAmplitude
Float_t GetAmplitude() const
Definition BmnTof1Digit.h:26
TCanvas::cd
TVirtualPad * cd(Int_t subpadnumber=0) override
TCanvas::Update
void Update() override
TChain::SetBranchAddress
Int_t SetBranchAddress(const char *bname, T **add, TBranch **ptr=nullptr)
TChain::Add
virtual Int_t Add(const char *name, Long64_t nentries=TTree::kMaxEntries)
TChain::GetEntries
Long64_t GetEntries() const override
TChain::GetEntry
Int_t GetEntry(Long64_t entry=0, Int_t getall=0) override
TF1::GetParameter
virtual Double_t GetParameter(const TString &name) const
TF1::SetParameter
virtual void SetParameter(const TString &name, Double_t value)
TF1::FixParameter
virtual void FixParameter(Int_t ipar, Double_t value)
TH1::GetBinCenter
virtual Double_t GetBinCenter(Int_t bin) const
TH1::GetMean
virtual Double_t GetMean(Int_t axis=1) const
TH1::Fit
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Double_t xmin=0, Double_t xmax=0)
TH1::GetMaximum
virtual Double_t GetMaximum(Double_t maxval=FLT_MAX) const
TH1::GetNbinsX
virtual Int_t GetNbinsX() const
TH1::Draw
void Draw(Option_t *option="") override
TH1::Fill
virtual Int_t Fill(const char *name, Double_t w)
TH1::GetRMS
Double_t GetRMS(Int_t axis=1) const
TH1::GetEntries
virtual Double_t GetEntries() const
TH1::GetFunction
virtual TF1 * GetFunction(const char *name) const
TH1::Integral
virtual Double_t Integral(Int_t binx1, Int_t binx2, Option_t *option="") const
TH1::GetMaximumBin
virtual Int_t GetMaximumBin() const
TH1::GetBinContent
virtual Double_t GetBinContent(Int_t bin) const
TH1::GetBinWidth
virtual Double_t GetBinWidth(Int_t bin) const
TH1::SetAxisRange
virtual void SetAxisRange(Double_t xmin, Double_t xmax, Option_t *axis="X")
TH1::FindBin
virtual Int_t FindBin(Double_t x, Double_t y=0, Double_t z=0)
TH2::Fill
virtual Int_t Fill(const char *namex, const char *namey, Double_t w)
TNamed::GetTitle
const char * GetTitle() const override
TObjArray::GetEntriesFast
Int_t GetEntriesFast() const
TObjArray::At
TObject * At(Int_t idx) const override
TObject::ResetBit
void ResetBit(UInt_t f)
TPad::Divide
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
TProfile::Fill
Int_t Fill(const char *namex, Double_t y) override
TVirtualPad::Draw
void Draw(Option_t *option="") override=0
IN_WIDTH_LIMITS
#define IN_WIDTH_LIMITS
Definition BmnResultsTOF700.C:19
STRIP_1
#define STRIP_1
Definition BmnResultsTOF700.C:33
STRIP_2
#define STRIP_2
Definition BmnResultsTOF700.C:34
WIDTH_MAX
#define WIDTH_MAX
Definition BmnResultsTOF700.C:20
draw_vs_cham
void draw_vs_cham(char *, TH1F **)
Definition BmnResultsTOF700.C:681
wma
int wma[NUMBER_CHAMBERS]
Definition BmnResultsTOF700.C:48
champos_left
int champos_left[30]
Definition BmnResultsTOF700.C:62
FitIn2Sigma
double FitIn2Sigma(TH1F *h)
Definition BmnResultsTOF700.C:474
NUMBER_CHAMBERS
#define NUMBER_CHAMBERS
Definition BmnResultsTOF700.C:18
draw_vs_cham_prof
void draw_vs_cham_prof(char *, TProfile **)
Definition BmnResultsTOF700.C:709
idchambers
float idchambers[59]
Definition BmnResultsTOF700.C:52
STRIP_1_2
#define STRIP_1_2
Definition BmnResultsTOF700.C:35
DIFF_CHAMB_1
#define DIFF_CHAMB_1
Definition BmnResultsTOF700.C:29
BmnResultsTOF700
void BmnResultsTOF700(const char *fname="./5070/bmn_run5070_digi.root")
Definition BmnResultsTOF700.C:78
DIFF_CHAMB_2
#define DIFF_CHAMB_2
Definition BmnResultsTOF700.C:30
champos
int champos[NUMBER_CHAMBERS]
Definition BmnResultsTOF700.C:46
chamnum_right
int chamnum_right[NUMBER_CHAMBERS - 30]
Definition BmnResultsTOF700.C:60
wmi
int wmi[NUMBER_CHAMBERS]
Definition BmnResultsTOF700.C:47
STRIP_2_2
#define STRIP_2_2
Definition BmnResultsTOF700.C:36
chamnum_left
int chamnum_left[30]
Definition BmnResultsTOF700.C:58
FitIn
double FitIn(TH1F *h, Double_t xmin, Double_t xmax)
Definition BmnResultsTOF700.C:588
champos_right
int champos_right[NUMBER_CHAMBERS - 30]
Definition BmnResultsTOF700.C:65
select_hist
void select_hist()
Definition BmnResultsTOF700.C:623