source: trunk/MagicSoft/Mars/mtemp/mifae/macros/lightcurve.C@ 4393

Last change on this file since 4393 was 4387, checked in by jlopez, 21 years ago
*** empty log message ***
File size: 49.8 KB
Line 
1Double_t ChiSquareNDof(TH1D *h1, TH1D *h2)
2{
3 Double_t chiq = 0.;
4 Double_t chi;
5 Double_t error;
6 Int_t nbinsnozero = 0;
7
8 Int_t nbins = h1->GetNbinsX();
9 if (nbins != h2->GetNbinsX() || nbins == 0)
10 return -1;
11
12 for (UInt_t bin=1; bin<=nbins; bin++)
13 {
14 error = sqrt(h1->GetBinError(bin)*h1->GetBinError(bin) +
15 h2->GetBinError(bin)*h2->GetBinError(bin));
16 if (error != 0)
17 {
18 chi = (h1->GetBinContent(bin)-h2->GetBinContent(bin))/error;
19 chiq += chi*chi;
20 nbinsnozero++;
21 }
22 }
23
24 return chiq/nbinsnozero;
25}
26
27Int_t GetBin(Double_t size, Int_t numberSizeBins, Double_t sizeBins[])
28{
29
30 Int_t result = -1;
31
32 Int_t lowerbin = 0;
33 Int_t upperbin = numberSizeBins;
34 Int_t bin;
35
36 Int_t count = 0;
37
38 if (size >= sizeBins[0])
39 {
40 while (upperbin - lowerbin > 1 && count++<=numberSizeBins)
41 {
42 bin = (upperbin+lowerbin)/2;
43 if (size >= sizeBins[bin])
44 lowerbin = bin;
45 else
46 upperbin = bin;
47 }
48 result = count<=numberSizeBins?lowerbin:-1;
49 }
50
51 return result;
52
53}
54
55
56void lightcurve(TString f_on_name = "../HillasFiles/Mrk421/*_H.root",
57 TString f_off_name = "../HillasFiles/OffMrk421/*_H.root")
58{
59
60 gROOT->Reset();
61 gStyle->SetTimeOffset(-3600);
62
63 // Configuration
64 const Bool_t debug = kFALSE;
65 const Double_t timebin = 480.; //[sec]
66 TString psname = "./20040422_Mrk421." + timebin + "s.ps";
67 cout << psname << endl;
68 return;
69
70 //Constanst
71 const Double_t kConvDegToRad = TMath::ACos(-1)/180.;
72 const Double_t kSec = 1e6; //[sec/microsec]
73 const Double_t kDay = 24.*60.*60.; //[Day/sec]
74
75 UInt_t numberTimeBins = 1;
76 TArrayI numberOnEvents(1);
77 TArrayD numberBkgEventsToNormOn(1);
78 TArrayD timeOn(1);
79
80 TArrayD meanTimeBinOn(1);
81 TArrayD widthTimeBinOn(1);
82
83 TArrayD zenithMinimumOn(1);
84 TArrayD zenithMaximumOn(1);
85
86 TObjArray coszenithHistoOn;
87 TObjArray alphaHistoOn;
88 TObjArray srcposHistoOn;
89
90
91 //none cuts
92// Double_t sizemin = 2000.; //[Photons]
93// Double_t sizemax = 10000000000.; //[Photons]
94// Double_t widthmin = 0.;
95// Double_t widthmax = 2.0;
96// Double_t lengthmin = 0.;
97// Double_t lengthmax = 2.0;
98// Double_t distmin = 0.;
99// Double_t distmax = 4.0;
100// Double_t alphamin = 0.;
101// Double_t alphamax = 90.;
102
103
104 const Int_t numberSizeBins = 4;
105 Double_t sizeBins[numberSizeBins] = {1292., 1897., 2785., 4087.}; //[Photons]
106
107 //cuts
108 Double_t sizemin = 2000.; //[Photons]
109 Double_t sizemax = 10E9.; //[Photons]
110
111 Double_t widthmin[numberSizeBins] = { 0.06, 0.06, 0.06, 0.06 };
112 Double_t widthmax[numberSizeBins] = { 0.10, 0.10, 0.12, 0.12 };
113 Double_t lengthmin[numberSizeBins] = { 0.10, 0.10, 0.10, 0.10 };
114 Double_t lengthmax[numberSizeBins] = { 0.20, 0.25, 0.26, 0.36 };
115 Double_t distmin[numberSizeBins] = { 0.30, 0.30, 0.30, 0.30,};
116 Double_t distmax[numberSizeBins] = { 1.20, 1.20, 1.20, 1.20,};
117
118 //alpha plot integration for gammas
119 Double_t sigexccmin = 0.;
120 Double_t sigexccmax = 20.;
121 Double_t bkgnormmin = 40.;
122 Double_t bkgnormmax = 80.;
123
124 gStyle->SetOptStat(111111);
125 gStyle->SetOptFit();
126
127 //
128 // Make a loop only for the ON data:
129 //
130
131 MParList plist_on;
132 MTaskList tlist_on;
133 plist_on.AddToList(&tlist_on);
134
135 // ON containers
136 MGeomCamMagic geomcam;
137 MSrcPosCam source_on;
138 MHillas hillas;
139 MHillasSrc hillasscr;
140
141 plist_on.AddToList(&geomcam);
142 plist_on.AddToList(&source_on);
143 plist_on.AddToList(&hillas);
144 plist_on.AddToList(&hillasscr);
145
146 //create some 1-dim histo to test only for the ON distribution of dist, width , length, size...
147 MH3 hDist_on("MHillasSrc.fDist/MGeomCam.fConvMm2Deg");
148 hDist_on.SetName("Dist_on");
149 plist_on.AddToList(&hDist_on);
150 MBinning binsDist_on("BinningDist_on");
151 Int_t nbins_Dist = 20;
152 Double_t min_Dist = 0.;
153 Double_t max_Dist = distmax[0]*1.2;
154 binsDist_on.SetEdges(nbins_Dist, min_Dist, max_Dist);
155 plist_on.AddToList(&binsDist_on);
156
157 MH3 hWidth_on("MHillas.fWidth/MGeomCam.fConvMm2Deg");
158 hWidth_on.SetName("Width_on");
159 plist_on.AddToList(&hWidth_on);
160 MBinning binsWidth_on("BinningWidth_on");
161 Int_t nbins_Width = 20;
162 Double_t min_Width = 0.;
163 Double_t max_Width = widthmax[0]*1.2;
164 binsWidth_on.SetEdges(nbins_Width, min_Width, max_Width);
165 plist_on.AddToList(&binsWidth_on);
166
167 MH3 hLength_on("MHillas.fLength/MGeomCam.fConvMm2Deg");
168 hLength_on.SetName("Length_on");
169 plist_on.AddToList(&hLength_on);
170 MBinning binsLength_on("BinningLength_on");
171 Int_t nbins_Length = 20;
172 Double_t min_Length = 0.;
173 Double_t max_Length = lengthmax[0]*1.2;
174 binsLength_on.SetEdges(nbins_Length, min_Length, max_Length);
175 plist_on.AddToList(&binsLength_on);
176
177 MH3 hSize_on("log10(MHillas.fSize)");
178 hSize_on.SetName("Size_on");
179 plist_on.AddToList(&hSize_on);
180 MBinning binsSize_on("BinningSize_on");
181 Int_t nbins_Size = 60;
182 Double_t min_Size = log10(sizemin)*0.8;
183 Double_t max_Size = log10(1000000)*1.2;
184 binsSize_on.SetEdges(nbins_Size, min_Size, max_Size);
185 plist_on.AddToList(&binsSize_on);
186
187 //create a histo to fill the alpha values: one alpha plot form 0 to +90 deg in abs value
188 MH3 hAlpha_on_abs("abs(MHillasSrc.fAlpha)");
189 hAlpha_on_abs.SetName("Alpha_on_abs");
190 plist_on.AddToList(&hAlpha_on_abs);
191 MBinning binsAlpha_on_abs("BinningAlpha_on_abs");
192 Int_t nbins_abs = 18;
193 Double_t minalpha_abs = 0.;
194 Double_t maxalpha_abs =90.;
195 binsAlpha_on_abs.SetEdges(nbins_abs, minalpha_abs, maxalpha_abs);
196 plist_on.AddToList(&binsAlpha_on_abs);
197
198 //create a histo to fill the alpha values: one alpha plot form -90 to +90 deg.
199 MH3 hAlpha_on("MHillasSrc.fAlpha");
200 hAlpha_on.SetName("Alpha_on");
201 plist_on.AddToList(&hAlpha_on);
202 MBinning binsAlpha_on("BinningAlpha_on");
203 Int_t nbins = nbins_abs*2;
204 Double_t minalpha = -90.;
205 Double_t maxalpha = 90.;
206 binsAlpha_on.SetEdges(nbins, minalpha, maxalpha);
207 plist_on.AddToList(&binsAlpha_on);
208
209
210 //create a histo to fill the source position
211 MH3 hSrcPos_on("MSrcPosCam.fX","MSrcPosCam.fY");
212 hSrcPos_on.SetName("SrcPos_on");
213 plist_on.AddToList(&hSrcPos_on);
214 MBinning binsSrcPos_onX("BinningSrcPos_onX");
215 MBinning binsSrcPos_onY("BinningSrcPos_onY");
216 Int_t nbins_srcpos = 400;
217 Double_t minsrcpos = -600.;
218 Double_t maxsrcpos = 600.;
219 binsSrcPos_onX.SetEdges(nbins_srcpos, minsrcpos, maxsrcpos);
220 binsSrcPos_onY.SetEdges(nbins_srcpos, minsrcpos, maxsrcpos);
221 plist_on.AddToList(&binsSrcPos_onX);
222 plist_on.AddToList(&binsSrcPos_onY);
223
224 MH3 hDAQEvtNumber_on("MRawEvtHeader.fDAQEvtNumber");
225 hDAQEvtNumber_on.SetName("DAQEvtNumber_on");
226 plist_on.AddToList(&hDAQEvtNumber_on);
227 MBinning binsDAQEvtNumber_onX("BinningDAQEvtNumber_onX");
228 Int_t nbins_evtnum = 1000;
229 Double_t minevtnum = 0.;
230 Double_t maxevtnum = 1000.;
231 binsDAQEvtNumber_onX.SetEdges(nbins_evtnum,minevtnum,maxevtnum);
232 plist_on.AddToList(&binsDAQEvtNumber_onX);
233
234
235 //
236 //tasks
237 //
238
239 MReadTree read_on("Parameters", f_on_name);
240 read_on.DisableAutoScheme();
241
242 MSrcPlace srcplace;
243
244 //cuts
245 TString sizestr = "(MHillas.fSize < ";
246 sizestr += sizemin;
247 sizestr += ") || (";
248 sizestr += "MHillas.fSize > ";
249 sizestr += sizemax;
250 sizestr += ")";
251 MF sizefilter(sizestr);
252
253 TString widthstr = "({MHillas.fWidth/MGeomCam.fConvMm2Deg} < ";
254 widthstr += widthmin[0];
255 widthstr += ") || (";
256 widthstr += "{MHillas.fWidth/MGeomCam.fConvMm2Deg} > ";
257 widthstr += widthmax[0];
258 widthstr += ")";
259 MF widthfilter(widthstr);
260
261 TString lengthstr = "({MHillas.fLength/MGeomCam.fConvMm2Deg} < ";
262 lengthstr += lengthmin[0];
263 lengthstr += ") || (";
264 lengthstr += "{MHillas.fLength/MGeomCam.fConvMm2Deg} > ";
265 lengthstr += lengthmax[0];
266 lengthstr += ")";
267 MF lengthfilter(lengthstr);
268
269 TString diststr = "({MHillasSrc.fDist/MGeomCam.fConvMm2Deg} < ";
270 diststr += distmin[0];
271 diststr += ") || (";
272 diststr += "{MHillasSrc.fDist/MGeomCam.fConvMm2Deg} > ";
273 diststr += distmax[0];
274 diststr += ")";
275 MF distfilter(diststr);
276
277 MF evenfilter("{MRawEvtHeader.fDAQEvtNumber%3}<0.5");
278 MF oddfilter("{MRawEvtHeader.fDAQEvtNumber%3}>0.5");
279
280 MContinue cont_size(&sizefilter);
281 MContinue cont_width(&widthfilter);
282 MContinue cont_length(&lengthfilter);
283 MContinue cont_dist(&distfilter);
284 MContinue cont_even(&evenfilter);
285 MContinue cont_odd(&oddfilter);
286
287 MHillasSrcCalc csrc_on;
288
289 // fill all histograms
290 MFillH falpha_on_abs(&hAlpha_on_abs);
291 MFillH falpha_on(&hAlpha_on);
292 MFillH fdist_on(&hDist_on);
293 MFillH fwidth_on(&hWidth_on);
294 MFillH flength_on(&hLength_on);
295 MFillH fsize_on(&hSize_on);
296 MFillH fsrcpos_on(&hSrcPos_on);
297 MFillH fevtnum_on(&hDAQEvtNumber_on);
298
299 // prints
300 MPrint pevent("MRawEvtHeader");
301 MPrint phillas("MHillas");
302 MPrint phillassrc("MHillasSrc");
303 MPrint psrcpos("MSrcPosCam");
304
305 //tasklist
306 tlist_on.AddToList(&read_on);
307 tlist_on.AddToList(&srcplace);
308 tlist_on.AddToList(&csrc_on);
309 tlist_on.AddToList(&fsrcpos_on);
310 // tlist_on.AddToList(&cont_odd);
311 tlist_on.AddToList(&cont_size);
312 tlist_on.AddToList(&cont_width);
313 tlist_on.AddToList(&cont_length);
314 tlist_on.AddToList(&cont_dist);
315 tlist_on.AddToList(&falpha_on_abs);
316 tlist_on.AddToList(&falpha_on);
317 tlist_on.AddToList(&fdist_on);
318 tlist_on.AddToList(&fwidth_on);
319 tlist_on.AddToList(&flength_on);
320 tlist_on.AddToList(&fsize_on);
321 // tlist_on.AddToList(&fevtnum_on);
322
323 // Create and setup the eventloop
324 MEvtLoop loop_on;
325 loop_on.SetParList(&plist_on);
326 //loop_on.SetDisplay(display);
327
328// MProgressBar bar;
329// loop_on.SetProgressBar(&bar);
330
331// if (!loop_on.Eventloop())
332// return;
333
334 if (!loop_on.PreProcess())
335 return;
336
337 MRawRunHeader* runheader_on = plist_on.FindObject("MRawRunHeader");
338 MRawEvtHeader* evtheader_on = plist_on.FindObject("MRawEvtHeader");
339 MTime* time_on = plist_on.FindObject("MTime");
340 MHillas* hillas_on = plist_on.FindObject("MHillas");
341 MHillasSrc* hillassrc_on = plist_on.FindObject("MHillasSrc");
342 MReportDrive* drive_on = plist_on.FindObject("MReportDrive");
343 MSrcPosCam* srcpos_on = plist_on.FindObject("MSrcPosCam");
344
345
346 Double_t mjdFirstEventofBin=0;
347 Double_t mjdFirstEvent=0;
348
349 MTime mtimeFirstEventofBin;
350 MTime mtimeFirstEvent;
351
352 Double_t mjdLastEvent=0;
353 MTime mtimeLastEvent=0;
354 UInt_t evtLastEvent;
355 UInt_t runLastEvent=0;
356 MTime mtimeLastEvent;
357
358 Bool_t flag = kFALSE;
359
360 zenithMinimumOn[numberTimeBins-1] = 100.;
361 zenithMaximumOn[numberTimeBins-1] = 0.;
362 timeOn[numberTimeBins-1] = 0;
363
364 //create histos needed in the time bins
365 TString alphaTitle = Form("%s%02i","hAlphaOn",numberTimeBins-1);
366 TH1F *hAlpha_on_abs_timebin = new TH1F (alphaTitle,"",nbins_abs,minalpha_abs,maxalpha_abs);
367 Int_t nbins_srcpos = 200;
368 Double_t minsrcpos = -600.; //[mm]
369 Double_t maxsrcpos = 600.; //[mm] //!!! position precition 3mm ~ 0.01 deg
370 TString srcposTitle = Form("%s%02i","hSrcPosOn",numberTimeBins-1);
371 TH2F *hSrcPos_on_timebin = new TH2F (srcposTitle,"",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos);
372 Int_t nbins_coszenith = 200;
373 Double_t mincoszenith = 0.;
374 Double_t maxcoszenith = 1.;
375 TString coszenithTitle = Form("%s%02i","hCosZenithOn",numberTimeBins-1);
376 TH1F *hCosZenith_on_timebin = new TH1F (coszenithTitle,"",nbins_coszenith,mincoszenith,maxcoszenith);
377
378 while(tlist_on.Process())
379 {
380 // Compute live time
381 UInt_t run = runheader_on->GetRunNumber();
382 UInt_t evt = evtheader_on->GetDAQEvtNumber();
383 Double_t mjd = time_on->GetMjd();
384
385 if (mjd == 0)
386 {
387
388 if (debug)
389 {
390 cout << "MTime::GetTime() == 0" << endl;
391 cout.precision(15);
392 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
393 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl;
394 mtimeLastEvent.Print();
395 time_on->Print();
396 }
397
398 flag = kTRUE;
399 }
400 else if (mjd-mjdLastEvent < 0 && flag)
401 {
402
403 if (debug)
404 {
405 cout << "mjd-mjdLastEvent < 0" << endl;
406 cout.precision(15);
407 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
408 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl;
409 mtimeLastEvent.Print();
410 time_on->Print();
411 }
412
413 flag = kTRUE;
414 }
415 else if (mjd-mjdLastEvent == 0 && flag)
416 {
417
418 if (debug)
419 {
420 cout << "mjd-mjdLastEvent == 0" << endl;
421 cout.precision(15);
422 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
423 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl;
424 mtimeLastEvent.Print();
425 time_on->Print();
426 }
427
428 flag = kTRUE;
429 }
430 else if (evt-evtLastEvent <= 0)
431 {
432
433 if (debug)
434 {
435 cout << "evt-evtLastEvent <= 0" << endl;
436 cout.precision(15);
437 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
438 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl;
439 mtimeLastEvent.Print();
440 time_on->Print();
441 }
442
443 flag = kTRUE;
444 }
445
446 else if ((Int_t)mjd == mjd)
447 {
448
449 if (debug)
450 {
451 cout << "(Int_t)mjd == mjd" << endl;
452 cout.precision(15);
453 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
454 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl;
455 mtimeLastEvent.Print();
456 time_on->Print();
457 }
458
459 flag = kTRUE;
460 }
461 else
462 {
463
464 if (flag && debug)
465 {
466
467 cout << "flag = TRUE" << endl;
468 cout.precision(15);
469 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
470 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl;
471 mtimeLastEvent.Print();
472 time_on->Print();
473
474 flag = kFALSE;
475
476 }
477
478 if ( run != runLastEvent )
479 {
480 if ( mjdLastEvent != 0 )
481 {
482 timeOn[numberTimeBins-1] += (mjdLastEvent-mjdFirstEvent)*kDay;
483
484 cout.precision(10);
485 cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl;
486 if (flag && debug)
487 {
488 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
489 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl;
490 mtimeLastEvent.Print();
491 time_on->Print();
492 }
493
494 }
495
496 mjdFirstEvent = mjd;
497 }
498
499 if (mjdFirstEventofBin == 0)
500 {
501 mjdFirstEvent = mjd;
502 mjdFirstEventofBin = mjd;
503 }
504
505 evtLastEvent = evt;
506 runLastEvent = run;
507 mjdLastEvent = mjd;
508 mtimeLastEvent = *time_on;
509
510
511
512 Double_t size = hillas_on->GetSize();
513 Double_t width = hillas_on->GetWidth()*geomcam.GetConvMm2Deg();
514 Double_t length = hillas_on->GetLength()*geomcam.GetConvMm2Deg();
515 Double_t dist = hillassrc_on->GetDist()*geomcam.GetConvMm2Deg();
516 Double_t alpha = hillassrc_on->GetAlpha();
517 Double_t srcposx = srcpos_on->GetX();
518 Double_t srcposy = srcpos_on->GetY();
519 Double_t zenith = drive_on->GetNominalZd();
520
521
522 Int_t sizebin = GetBin(size,numberSizeBins,sizeBins);
523
524 if (sizebin >= 0)
525 {
526 if (width > widthmin[sizebin] && width < widthmax[sizebin])
527 {
528 if (length > lengthmin[sizebin] && length < lengthmax[sizebin])
529 {
530 if (dist > distmin[sizebin] && dist < distmax[sizebin])
531 {
532 hAlpha_on_abs_timebin->Fill(TMath::Abs(alpha));
533 hSrcPos_on_timebin->Fill(srcposx,srcposy);
534 hCosZenith_on_timebin->Fill(TMath::Cos(zenith*kConvDegToRad));
535
536 if (zenith != 0 && zenith < zenithMinimumOn[numberTimeBins-1])
537 zenithMinimumOn[numberTimeBins-1] = zenith;
538 if (zenith>zenithMaximumOn[numberTimeBins-1])
539 zenithMaximumOn[numberTimeBins-1] = zenith;
540
541 }
542 }
543 }
544 }
545
546
547 // Check if you are overload the maxim time bin
548 if ((mjd-mjdFirstEventofBin)*kDay >= timebin)
549 {
550 //Compute the time on
551 timeOn[numberTimeBins-1] += (mjdLastEvent-mjdFirstEvent)*kDay;
552 widthTimeBinOn[numberTimeBins-1] = (mjd-mjdFirstEventofBin)/2;
553 meanTimeBinOn[numberTimeBins-1] = mjdFirstEventofBin + widthTimeBinOn[numberTimeBins-1];
554
555 cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl;
556 cout << "mjd " << mjd << " mjdFirstEventofBin " << mjdFirstEventofBin << " widthTimeBinOn " << widthTimeBinOn[numberTimeBins-1] << " meanTimeBinOn " << meanTimeBinOn[numberTimeBins-1] << ' ' << mjdFirstEventofBin + widthTimeBinOn[numberTimeBins-1] << endl;
557
558 //Compute the number of on events
559 numberOnEvents[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
560 numberBkgEventsToNormOn[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
561 //Initialize variables
562 mjdFirstEvent = mjd;
563 mjdFirstEventofBin = mjd;
564
565 alphaHistoOn.AddLast(hAlpha_on_abs_timebin);
566 srcposHistoOn.AddLast(hSrcPos_on_timebin);
567 coszenithHistoOn.AddLast(hCosZenith_on_timebin);
568
569// cout << "hAlpha_on_abs_timebin " << hAlpha_on_abs_timebin << " alphaHistoOn [" << numberTimeBins-1 << "] " << alphaHistoOn[numberTimeBins-1] << endl;
570
571 //Increase the number of time bins and all needed arrays
572 numberTimeBins++;
573
574 timeOn.Set(numberTimeBins);
575 numberOnEvents.Set(numberTimeBins);
576 numberBkgEventsToNormOn.Set(numberTimeBins);
577 widthTimeBinOn.Set(numberTimeBins);
578 meanTimeBinOn.Set(numberTimeBins);
579 zenithMinimumOn.Set(numberTimeBins);
580 zenithMaximumOn.Set(numberTimeBins);
581
582 timeOn[numberTimeBins-1] = 0.;
583 zenithMinimumOn[numberTimeBins-1] = 100.;
584 zenithMaximumOn[numberTimeBins-1] = 0.;
585
586 alphaTitle = Form("%s%02i","hAlphaOn",numberTimeBins-1);
587 hAlpha_on_abs_timebin = new TH1F (alphaTitle,"",nbins_abs,minalpha_abs,maxalpha_abs);
588
589 srcposTitle = Form("%s%02i","hSrcPosOn",numberTimeBins-1);
590 hSrcPos_on_timebin = new TH2F (srcposTitle,"",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos);
591
592 coszenithTitle = Form("%s%02i","hCosZenithOn",numberTimeBins-1);
593 hCosZenith_on_timebin = new TH1F (coszenithTitle,"",nbins_coszenith,mincoszenith,maxcoszenith);
594
595 flag = kTRUE;
596 }
597 }
598 }
599
600
601 //Compute the time on for last time bin
602 timeOn[numberTimeBins-1] += (mjdLastEvent-mjdFirstEvent)*kDay;
603 widthTimeBinOn[numberTimeBins-1] = (mjd-mjdFirstEventofBin)/2;
604 meanTimeBinOn[numberTimeBins-1] = mjdFirstEventofBin + widthTimeBinOn[numberTimeBins-1];
605 //Compute the number of on events for the last time bin
606 numberOnEvents[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
607 numberBkgEventsToNormOn[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
608
609 cout.precision(10);
610 cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl;
611 cout << "mjd " << mjd << " mjdFirstEventofBin " << mjdFirstEventofBin << " meanTimeBinOn " << meanTimeBinOn[numberTimeBins-1] << " widthTimeBinOn " << widthTimeBinOn[numberTimeBins-1] << endl;
612
613 alphaHistoOn.AddLast(hAlpha_on_abs_timebin);
614 srcposHistoOn.AddLast(hSrcPos_on_timebin);
615 coszenithHistoOn.AddLast(hCosZenith_on_timebin);
616
617 // cout << "hAlpha_on_abs_timebin " << hAlpha_on_abs_timebin << " alphaHistoOn [" << numberTimeBins-1 << "] " << alphaHistoOn[numberTimeBins-1] << endl;
618
619 loop_on.PostProcess();
620
621 tlist_on.PrintStatistics();
622
623 for (UInt_t bin=0; bin<numberTimeBins; bin++)
624 cout << bin << " timeOn " << timeOn[bin] << " min-max zenithOn " << zenithMinimumOn[bin] << "-" << zenithMaximumOn[bin] << " min-max cos(zenithOn) " << TMath::Cos(zenithMinimumOn[bin]*kConvDegToRad) << "-" << TMath::Cos(zenithMaximumOn[bin]*kConvDegToRad) << " numberOnEvents " << numberOnEvents[bin] << endl;
625
626
627 //-----------------------OFF------------------------
628
629 TObjArray alphaHistoOff;
630 TObjArray srcposHistoOff;
631
632 TArrayD timeOff(numberTimeBins);
633
634 TArrayI numberOffEvents(numberTimeBins);
635
636 TArrayD numberBkgEventsToNormOff(numberTimeBins);
637 TArrayD onOffNormFactor(numberTimeBins);
638
639 TArrayD numberExcessEvents(numberTimeBins);
640 TArrayD errorNumberExcessEvents(numberTimeBins);
641
642 TArrayD numberExcessEventsPerSec(numberTimeBins);
643 TArrayD errorNumberExcessEventsPerSec(numberTimeBins);
644
645 TArrayD numberExcessEventsPerMin(numberTimeBins);
646 TArrayD errorNumberExcessEventsPerMin(numberTimeBins);
647
648 timeOff.Set(numberTimeBins);
649
650 TH1F* hAlpha_off_abs_timebin;
651 TH2F* hSrcPos_off_timebin;
652 for (UInt_t bin=0; bin<numberTimeBins; bin++)
653 {
654 alphaTitle = Form("%s%02i","hAlphaOff",bin);
655 hAlpha_off_abs_timebin = new TH1F (alphaTitle,"",nbins_abs,minalpha_abs,maxalpha_abs);;
656 alphaHistoOff.AddLast(hAlpha_off_abs_timebin);
657 // cout << "hAlpha_off_abs_timebin " << hAlpha_off_abs_timebin << " alphaHistoOff [" << bin << "] " << alphaHistoOff[bin] << endl;
658
659 srcposTitle = Form("%s%02i","hSrcPosOff",bin);
660 hSrcPos_off_timebin = new TH2F (srcposTitle,"",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos);
661 srcposHistoOff.AddLast(hSrcPos_off_timebin);
662 // cout << "hSrcPos_off_timebin " << hSrcPos_off_timebin << " srcposHistoOff [" << bin << "] " << srcposHistoOff[bin] << endl;
663 }
664
665 //
666 // Make a loop only for the OFF data:
667 //
668
669 MParList plist_off;
670 MTaskList tlist_off;
671 plist_off.AddToList(&tlist_off);
672
673 MSrcPosCam source_off;
674
675 plist_off.AddToList(&geomcam);
676 plist_off.AddToList(&source_off);
677 plist_off.AddToList(&hillas);
678 plist_off.AddToList(&hillasscr);
679
680 //create some 1-dim histo to test only for the OFF distribution of dist, width , length, size...
681 MH3 hDist_off("MHillasSrc.fDist/MGeomCam.fConvMm2Deg");
682 hDist_off.SetName("Dist_off");
683 plist_off.AddToList(&hDist_off);
684 MBinning binsDist_off("BinningDist_off");
685 binsDist_off.SetEdges(nbins_Dist, min_Dist, max_Dist);
686 plist_off.AddToList(&binsDist_off);
687
688 MH3 hWidth_off("MHillas.fWidth/MGeomCam.fConvMm2Deg");
689 hWidth_off.SetName("Width_off");
690 plist_off.AddToList(&hWidth_off);
691 MBinning binsWidth_off("BinningWidth_off");
692 binsWidth_off.SetEdges(nbins_Width, min_Width, max_Width);
693 plist_off.AddToList(&binsWidth_off);
694
695 MH3 hLength_off("MHillas.fLength/MGeomCam.fConvMm2Deg");
696 hLength_off.SetName("Length_off");
697 plist_off.AddToList(&hLength_off);
698 MBinning binsLength_off("BinningLength_off");
699 binsLength_off.SetEdges(nbins_Length, min_Length, max_Length);
700 plist_off.AddToList(&binsLength_off);
701
702 MH3 hSize_off("log10(MHillas.fSize)");
703 hSize_off.SetName("Size_off");
704 plist_off.AddToList(&hSize_off);
705 MBinning binsSize_off("BinningSize_off");
706 binsSize_off.SetEdges(nbins_Size, min_Size, max_Size);
707 plist_off.AddToList(&binsSize_off);
708
709 //create a histo to fill the alpha values: from 0 to 90 deg -> abs value
710 MH3 hAlpha_off_abs("abs(MHillasSrc.fAlpha)");
711 hAlpha_off_abs.SetName("Alpha_off_abs");
712 plist_off.AddToList(&hAlpha_off_abs);
713 MBinning binsAlpha_off_abs("BinningAlpha_off_abs");
714 binsAlpha_off_abs.SetEdges(nbins_abs, minalpha_abs, maxalpha_abs);
715 plist_off.AddToList(&binsAlpha_off_abs);
716
717 //create a histo to fill the alpha values: from -90 to 90 deg
718 MH3 hAlpha_off("MHillasSrc.fAlpha");
719 hAlpha_off.SetName("Alpha_off");
720 plist_off.AddToList(&hAlpha_off);
721 MBinning binsAlpha_off("BinningAlpha_off");
722 binsAlpha_off.SetEdges(nbins, minalpha, maxalpha);
723 plist_off.AddToList(&binsAlpha_off);
724
725
726 MH3 hSrcPos_off("MSrcPosCam.fX","MSrcPosCam.fY");
727 hSrcPos_off.SetName("SrcPos_off");
728 plist_off.AddToList(&hSrcPos_off);
729 MBinning binsSrcPos_offX("BinningSrcPos_offX");
730 MBinning binsSrcPos_offY("BinningSrcPos_offY");
731 binsSrcPos_offX.SetEdges(nbins_srcpos, minsrcpos, maxsrcpos);
732 binsSrcPos_offY.SetEdges(nbins_srcpos, minsrcpos, maxsrcpos);
733 plist_off.AddToList(&binsSrcPos_offX);
734 plist_off.AddToList(&binsSrcPos_offY);
735
736 MH3 hDAQEvtNumber_off("MRawEvtHeader.fDAQEvtNumber");
737 hDAQEvtNumber_off.SetName("DAQEvtNumber_off");
738 plist_off.AddToList(&hDAQEvtNumber_off);
739 MBinning binsDAQEvtNumber_offX("BinningDAQEvtNumber_offX");
740 Int_t nbins_evtnum = 100;
741 Double_t minevtnum = 0.;
742 Double_t maxevtnum = 100.;
743 binsDAQEvtNumber_offX.SetEdges(nbins_evtnum,minevtnum,maxevtnum);
744 plist_off.AddToList(&binsDAQEvtNumber_offX);
745
746 //tasks
747 MReadTree read_off("Parameters", f_off_name);
748 read_off.DisableAutoScheme();
749
750 // taks to set the src position in the off data for the time bins
751 // --------------------------------------------------------------
752 MSrcPlace srcplace_timebin;
753 srcplace_timebin.SetCreateHisto(kFALSE);
754 srcplace_timebin.SetMode(MSrcPlace::kOff);
755 // --------------------------------------------------------------
756
757 srcplace.SetMode(MSrcPlace::kOff);
758
759 MHillasSrcCalc csrc_off;
760
761 // fill all histograms
762 MFillH falpha_off_abs(&hAlpha_off_abs);
763 MFillH falpha_off(&hAlpha_off);
764 MFillH fdist_off(&hDist_off);
765 MFillH fwidth_off(&hWidth_off);
766 MFillH flength_off(&hLength_off);
767 MFillH fsize_off(&hSize_off);
768 MFillH fsrcpos_off(&hSrcPos_off);
769 MFillH fevtnum_off(&hDAQEvtNumber_off);
770
771 //tasklist
772 tlist_off.AddToList(&read_off);
773 tlist_off.AddToList(&srcplace);
774 tlist_off.AddToList(&csrc_off);
775 tlist_off.AddToList(&fsrcpos_off);
776// tlist_off.AddToList(&cont_even);
777// tlist_off.AddToList(&cont_size);
778// tlist_off.AddToList(&cont_width);
779// tlist_off.AddToList(&cont_length);
780// tlist_off.AddToList(&cont_dist);
781// tlist_off.AddToList(&falpha_off_abs);
782// tlist_off.AddToList(&falpha_off);
783// tlist_off.AddToList(&fdist_off);
784// tlist_off.AddToList(&fwidth_off);
785// tlist_off.AddToList(&flength_off);
786// tlist_off.AddToList(&fsize_off);
787// tlist_off.AddToList(&fevtnum_off);
788 tlist_off.AddToList(&srcplace_timebin); // This is just to run the preprocess correctely
789
790 // Create and setup the eventloop
791 MEvtLoop loop_off;
792 loop_off.SetParList(&plist_off);
793 //loop_off.SetDisplay(display);
794
795// MProgressBar bar_off;
796// loop_off.SetProgressBar(&bar_off);
797
798// if (!loop_off.Eventloop(numEntries))
799// return;
800
801 if (!loop_off.PreProcess())
802 return;
803
804
805 MHillas* hillas_off = plist_off.FindObject("MHillas");
806 MHillasSrc* hillassrc_off = plist_off.FindObject("MHillasSrc");
807 MReportDrive* drive_off = plist_off.FindObject("MReportDrive");
808 MSrcPosCam* srcpos_off = plist_off.FindObject("MSrcPosCam");
809
810 TArrayI sizecut(numberTimeBins);
811 TArrayI widthcut(numberTimeBins);
812 TArrayI lengthcut(numberTimeBins);
813 TArrayI distcut(numberTimeBins);
814 TArrayI alphacut(numberTimeBins);
815
816 sizecut.Reset(0);
817 widthcut.Reset(0);
818 lengthcut.Reset(0);
819 distcut.Reset(0);
820 alphacut.Reset(0);
821
822 while(tlist_off.Process())
823 {
824
825 //Now fill just the alpha plot in the time/zenith bins defined in on
826 Double_t zenith = drive_off->GetNominalZd();
827
828 for (UInt_t bin=0; bin<numberTimeBins; bin++)
829 {
830 srcplace_timebin.SetPositionHisto((TH2F*)srcposHistoOn[bin]);
831 srcplace_timebin.CallProcess();
832 csrc_off.CallProcess();
833
834 Double_t size = hillas_off->GetSize();
835 Double_t width = hillas_off->GetWidth()*geomcam.GetConvMm2Deg();
836 Double_t length = hillas_off->GetLength()*geomcam.GetConvMm2Deg();
837 Double_t dist = hillassrc_off->GetDist()*geomcam.GetConvMm2Deg();
838 Double_t alpha = hillassrc_off->GetAlpha();
839 Double_t srcposx = srcpos_off->GetX();
840 Double_t srcposy = srcpos_off->GetY();
841
842 Int_t sizebin = GetBin(size,numberSizeBins,sizeBins);
843
844 if (sizebin >= 0)
845 {
846 if (width > widthmin[sizebin] && width < widthmax[sizebin])
847 {
848 if (length > lengthmin[sizebin] && length < lengthmax[sizebin])
849 {
850 if (dist > distmin[sizebin] && dist < distmax[sizebin])
851 {
852 ((TH1F*)alphaHistoOff[bin])->Fill(TMath::Abs(alpha));
853 ((TH2F*)srcposHistoOff[bin])->Fill(srcposx,srcposy);
854 }
855 else
856 distcut[bin]++;
857 }
858 else
859 lengthcut[bin]++;
860 }
861 else
862 widthcut[bin]++;
863 }
864 }
865
866// if (zenith > zenithMaximumOn[numberTimeBins-1])
867// {
868// cout << "Exit loop because we arrive to zenith = " << zenith << endl;
869// break;
870// }
871 }
872
873 loop_off.PostProcess();
874
875 tlist_off.PrintStatistics();
876
877// for (UInt_t bin=0; bin<numberTimeBins; bin++)
878// {
879// cout << bin << ' ' << ((TH1F*)alphaHistoOff[bin])->GetEntries() << ' ' << ((TH2F*)srcposHistoOff[bin])->GetEntries() << endl;
880// cout << "sizecut " << sizecut[bin] << endl;
881// cout << "width " << widthcut[bin] << endl;
882// cout << "length " << lengthcut[bin] << endl;
883// cout << "dist " << distcut[bin] << endl;
884// cout << "alpha " << alphacut[bin] << endl;
885// }
886
887
888 TArrayD meanTimeBinOnInSec(numberTimeBins);
889 TArrayD widthTimeBinOnInSec(numberTimeBins);
890
891 for (UInt_t bin=0; bin<numberTimeBins; bin++)
892 {
893
894 meanTimeBinOnInSec[bin] = (meanTimeBinOn[bin]-(Int_t)meanTimeBinOn[bin]);
895 if (meanTimeBinOnInSec[bin] > 0.5)
896 meanTimeBinOnInSec[bin] = meanTimeBinOnInSec[bin] - 1 ;
897 meanTimeBinOnInSec[bin] *= kDay;
898
899 widthTimeBinOnInSec[bin] = widthTimeBinOn[bin]*kDay;
900
901 numberOffEvents[bin] = (Double_t) ((TH1F*)alphaHistoOff[bin])->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
902 numberBkgEventsToNormOff[bin] = (Double_t) ((TH1F*)alphaHistoOff[bin])->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
903 if (numberOffEvents[bin] != 0 && numberBkgEventsToNormOff[bin] != 0)
904 {
905 onOffNormFactor[bin] = numberBkgEventsToNormOn[bin]/numberBkgEventsToNormOff[bin];
906 numberExcessEvents[bin] = numberOnEvents[bin] - onOffNormFactor[bin]*numberOffEvents[bin];
907 errorNumberExcessEvents[bin] = TMath::Sqrt(numberOnEvents[bin] + onOffNormFactor[bin]*onOffNormFactor[bin]*numberOffEvents[bin]);
908 numberExcessEventsPerSec[bin] = numberExcessEvents[bin]/timeOn[bin];
909 errorNumberExcessEventsPerSec[bin] = errorNumberExcessEvents[bin]/timeOn[bin];
910 numberExcessEventsPerMin[bin] = 60.*numberExcessEvents[bin]/timeOn[bin];
911 errorNumberExcessEventsPerMin[bin] = 60.*errorNumberExcessEvents[bin]/timeOn[bin];
912 }
913 }
914
915 for (UInt_t bin=0; bin<numberTimeBins; bin++)
916 {
917 cout.precision(5);
918 cout << bin << " timeOn " << timeOn[bin] << " mean-width time " << meanTimeBinOn[bin] << "-" << widthTimeBinOn[bin] << endl;
919 cout << " numberOnEvents\t " << numberOnEvents[bin] << endl;
920 cout << " numberOffEvents\t " << numberOffEvents[bin] << endl;
921 cout << " numberBkgEventsToNormOn\t " << numberBkgEventsToNormOn[bin] << endl;
922 cout << " numberBkgEventsToNormOff\t " << numberBkgEventsToNormOff[bin] << endl;
923 cout << " onOffNormFactor\t " << onOffNormFactor[bin] << endl;
924 cout << " numberExcessEvents\t " << numberExcessEvents[bin] << " +- " << errorNumberExcessEvents[bin] << endl;
925 cout << " Excess/Sec\t " << numberExcessEventsPerSec[bin] << " +- " << errorNumberExcessEventsPerSec[bin] << endl;
926 }
927
928
929
930 TString openpsname = psname + "(";
931 TString closepsname = psname + ")";
932
933 TCanvas *c1 = new TCanvas;
934 c1->cd(1);
935 TGraphErrors* lightcurvegraph = new TGraphErrors(numberTimeBins,meanTimeBinOnInSec.GetArray(),numberExcessEventsPerMin.GetArray(),widthTimeBinOnInSec.GetArray(),errorNumberExcessEventsPerMin.GetArray());
936 lightcurvegraph->SetTitle("LightCurve");
937 lightcurvegraph->SetMinimum(0.);
938 lightcurvegraph->SetMarkerStyle(21);
939 lightcurvegraph->SetMarkerSize(0.03);
940 lightcurvegraph->Draw("AP");
941 lightcurvegraph->GetXaxis()->SetTimeDisplay(1);
942 c1->Print(openpsname);
943
944 TCanvas** c = new TCanvas*[numberTimeBins];
945
946 for (UInt_t bin=0; bin<numberTimeBins-1; bin++)
947 {
948 c[bin] = new TCanvas;
949 c[bin]->cd(1);
950
951 ((TH1F*)alphaHistoOn[bin])->Sumw2();
952 ((TH1F*)alphaHistoOff[bin])->SetStats(0);
953 ((TH1F*)alphaHistoOff[bin])->Sumw2();
954 ((TH1F*)alphaHistoOff[bin])->Scale(onOffNormFactor[bin]);
955 ((TH1F*)alphaHistoOn[bin])->SetStats(0); //-> Do NOT show the legend with statistics
956 ((TH1F*)alphaHistoOn[bin])->SetLineColor(kBlack);
957 ((TH1F*)alphaHistoOn[bin])->SetMarkerStyle(21);
958 ((TH1F*)alphaHistoOn[bin])->SetMarkerColor(kBlack);
959 ((TH1F*)alphaHistoOn[bin])->SetMarkerSize(0.7);
960 ((TH1F*)alphaHistoOff[bin])->SetFillColor(46);
961 ((TH1F*)alphaHistoOff[bin])->SetLineColor(46);
962 ((TH1F*)alphaHistoOff[bin])->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
963 ((TH1F*)alphaHistoOff[bin])->SetMinimum(0.);
964 alphaTitle = Form("%s%02i","hAlphaOnOff",bin);
965 ((TH1F*)alphaHistoOn[bin])->SetTitle(alphaTitle);
966
967
968 ((TH1F*)alphaHistoOn[bin])->DrawCopy("E1P");
969 ((TH1F*)alphaHistoOff[bin])->DrawCopy("HISTSAME");
970 ((TH1F*)alphaHistoOff[bin])->DrawCopy("ESAME");
971 ((TH1F*)alphaHistoOn[bin])->DrawCopy("E1PSAME");
972
973 c[bin]->Print(psname);
974 }
975
976 c[numberTimeBins-1] = new TCanvas;
977 c[numberTimeBins-1]->cd(1);
978
979 ((TH1F*)alphaHistoOn[numberTimeBins-1])->Sumw2();
980 ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetStats(0);
981 ((TH1F*)alphaHistoOff[numberTimeBins-1])->Sumw2();
982 ((TH1F*)alphaHistoOff[numberTimeBins-1])->Scale(onOffNormFactor[numberTimeBins-1]);
983 ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetStats(0); //-> Do NOT show the legend with statistics
984 ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetLineColor(kBlack);
985 ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetMarkerStyle(21);
986 ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetMarkerColor(kBlack);
987 ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetMarkerSize(0.7);
988 ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetFillColor(46);
989 ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetLineColor(46);
990 ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
991 ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetMinimum(0.);
992 alphaTitle = Form("%s%02i","hAlphaOnOff",numberTimeBins-1);
993 ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetTitle(alphaTitle);
994
995
996 ((TH1F*)alphaHistoOn[numberTimeBins-1])->DrawCopy("E1P");
997 ((TH1F*)alphaHistoOff[numberTimeBins-1])->DrawCopy("HISTSAME");
998 ((TH1F*)alphaHistoOff[numberTimeBins-1])->DrawCopy("ESAME");
999 ((TH1F*)alphaHistoOn[numberTimeBins-1])->DrawCopy("E1PSAME");
1000
1001 c[numberTimeBins-1]->Print(closepsname);
1002
1003 TFile input("./prueba.root", "RECREATE");
1004
1005 for (UInt_t bin=0; bin<numberTimeBins; bin++)
1006 {
1007 ((TH1F*)alphaHistoOn[bin])->Write();
1008 ((TH1F*)alphaHistoOff[bin])->Write();
1009 ((TH2F*)srcposHistoOn[bin])->Write();
1010 ((TH2F*)srcposHistoOff[bin])->Write();
1011 ((TH1F*)coszenithHistoOn[bin])->Write();
1012 }
1013
1014 input.Close();
1015
1016 if(kFALSE)
1017 {
1018
1019 // ############################################################################
1020 // look for the histograms
1021 // ############################################################################
1022
1023 TH1F *hist_size_on = (TH1F*)hSize_on.GetHist();
1024 TH1F *hist_size_off = (TH1F*)hSize_off.GetHist();
1025
1026 TH1F *hist_dist_on = (TH1F*)hDist_on.GetHist();
1027 TH1F *hist_dist_off = (TH1F*)hDist_off.GetHist();
1028
1029 TH1F *hist_width_on = (TH1F*)hWidth_on.GetHist();
1030 TH1F *hist_width_off = (TH1F*)hWidth_off.GetHist();
1031
1032 TH1F *hist_length_on = (TH1F*)hLength_on.GetHist();
1033 TH1F *hist_length_off = (TH1F*)hLength_off.GetHist();
1034
1035 TH1F *hist_on_abs = (TH1F*)hAlpha_on_abs.GetHist();
1036 TH1F *hist_off_abs = (TH1F*)hAlpha_off_abs.GetHist();
1037
1038 TH1F *hist_on = (TH1F*)hAlpha_on.GetHist();
1039 TH1F *hist_off = (TH1F*)hAlpha_off.GetHist();
1040
1041
1042 // ############################################################################
1043 // Calculate significance and excess:
1044 // ############################################################################
1045
1046 Double_t norm_on_abs = (Double_t) hist_on_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
1047 Double_t exces_on_abs = (Double_t) hist_on_abs->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
1048 Double_t norm_off_abs = (Double_t) hist_off_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
1049 Double_t exces_off_abs = (Double_t) hist_off_abs->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
1050 Double_t norm = norm_on_abs/norm_off_abs;
1051
1052 char text_tit_alpha[256];
1053 sprintf(text_tit_alpha, " Alpha Plot On and Off ");
1054 hist_off_abs->SetTitle(text_tit_alpha);
1055 hist_on_abs->SetTitle(text_tit_alpha);
1056
1057 Double_t excess = exces_on_abs - exces_off_abs*norm;
1058 Double_t sign = excess / sqrt( exces_on_abs + norm*norm*exces_off_abs );
1059 Double_t int_off = (Double_t) hist_off_abs->Integral(1, 18);
1060 int hist_on_entries = (int) hist_on_abs->GetEntries();
1061 int hist_off_entries = (int) hist_off_abs->GetEntries();
1062
1063 cout << "---> Normalization F factor =\t" << norm <<endl;
1064 cout << "---> Excess =\t\t\t" << excess <<endl;
1065 cout << "---> Significancia =\t\t" << sign <<endl;
1066 cout << "---> entries on =\t\t" << hist_on_entries <<endl;
1067 cout << "---> entries off =\t\t" << hist_off_entries <<endl;
1068 cout << "---> integral off =\t\t" << int_off <<endl;
1069
1070 Double_t shiftx;
1071
1072 //
1073 //Create the display -> from now on, all histos are plotted
1074 MStatusDisplay *display = new MStatusDisplay;
1075 display->SetUpdateTime(3000);
1076 display->Resize(850,700);
1077
1078 // ############################################################################
1079 // Draw SIZE
1080 // ############################################################################
1081 display->AddTab("SIZE");
1082
1083 gPad->cd();
1084
1085 gPad->SetLogy();
1086 hist_size_on->Sumw2();
1087 hist_size_off->Sumw2();
1088 hist_size_off->Scale(norm);
1089 hist_size_on->SetLineColor(kBlack);
1090 hist_size_on->SetMarkerStyle(21);
1091 hist_size_on->SetMarkerSize(0.7);
1092 hist_size_on->SetMarkerColor(kBlack);
1093 hist_size_off->SetFillColor(46);
1094 hist_size_off->SetLineColor(46);
1095 hist_size_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
1096 hist_size_off->SetMinimum(0.1);
1097 hist_size_on->SetMinimum(0.1);
1098 hist_size_on->SetTitle("SIZE distribution");
1099 hist_size_off->SetTitle("SIZE distribution");
1100
1101 hist_size_on->DrawCopy("E1P");
1102
1103 // move stat box to make them all visible
1104 gPad->Update();
1105 TPaveStats* pavs_on_size = (TPaveStats*) hist_size_on->GetListOfFunctions()->FindObject("stats");
1106 if(pavs_on_size){
1107 shiftx = pavs_on_size->GetX2NDC() - pavs_on_size->GetX1NDC();
1108 pavs_on_size->SetX1NDC(pavs_on_size->GetX1NDC() - shiftx);
1109 pavs_on_size->SetX2NDC(pavs_on_size->GetX2NDC() - shiftx);
1110 }
1111 gPad->Modified();
1112 gPad->Update();
1113
1114 hist_size_off->DrawCopy("HISTSAME");
1115 hist_size_off->DrawCopy("ESAME");
1116
1117 gPad->Modified();
1118 gPad->Update();
1119
1120 Double_t chisize = ChiSquareNDof((TH1D*)hist_size_on,(TH1D*)hist_size_off);
1121
1122 Double_t x_label_pos = log10(1000000)*0.7;
1123 Double_t y_label_pos = log10((hist_size_on->GetBinContent(hist_size_on->GetMaximumBin()))/2.);
1124 Double_t textsize = 0.03;
1125
1126 char text_size[256];
1127 sprintf(text_size,"ChiSquare/NDof = %4.2f",chisize);
1128
1129 TLatex *tsize = new TLatex(x_label_pos, y_label_pos, text_size);
1130 tsize->SetTextSize(textsize);
1131// tsize->Draw();
1132
1133 gPad->Modified();
1134 gPad->Update();
1135
1136 // ############################################################################
1137 // DrawCopy DIST
1138 // ############################################################################
1139 display->AddTab("DIST");
1140
1141 gPad->cd();
1142
1143 hist_dist_on->Sumw2();
1144 hist_dist_off->Sumw2();
1145 hist_dist_off->Scale(norm);
1146 hist_dist_on->SetLineColor(kBlack);
1147 hist_dist_on->SetMarkerStyle(21);
1148 hist_dist_on->SetMarkerSize(0.7);
1149 hist_dist_on->SetMarkerColor(kBlack);
1150 hist_dist_off->SetFillColor(46);
1151 hist_dist_off->SetLineColor(46);
1152 hist_dist_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
1153 hist_dist_off->SetMinimum(0.);
1154 hist_dist_on->SetTitle("DIST distribution");
1155 hist_dist_off->SetTitle("DIST distribution");
1156
1157 hist_dist_on->DrawCopy("E1P");
1158
1159 // move stat box to make them all visible
1160 gPad->Update();
1161 TPaveStats* pavs_on_dist = (TPaveStats*) hist_dist_on->GetListOfFunctions()->FindObject("stats");
1162 if(pavs_on_dist){
1163 shiftx = pavs_on_dist->GetX2NDC() - pavs_on_dist->GetX1NDC();
1164 pavs_on_dist->SetX1NDC(pavs_on_dist->GetX1NDC() - shiftx);
1165 pavs_on_dist->SetX2NDC(pavs_on_dist->GetX2NDC() - shiftx);
1166 }
1167 gPad->Modified();
1168 gPad->Update();
1169
1170 hist_dist_off->DrawCopy("HISTSAME");
1171 hist_dist_off->DrawCopy("ESAME");
1172 hist_dist_on->DrawCopy("E1PSAME");
1173
1174 Double_t chidist = ChiSquareNDof((TH1D*)hist_dist_on,(TH1D*)hist_dist_off);
1175
1176 x_label_pos = distmax[0]*0.7;
1177 y_label_pos = hist_dist_on->GetBinContent(hist_dist_on->GetMaximumBin())/2.;
1178
1179 char text_dist[256];
1180 sprintf(text_size,"ChiSquare/NDof = %4.2f",chidist);
1181
1182 TLatex *tdist = new TLatex(x_label_pos, y_label_pos, text_dist);
1183 tdist->SetTextSize(textsize);
1184// tdist->Draw();
1185
1186 gPad->Modified();
1187 gPad->Update();
1188
1189 // ############################################################################
1190 // DrawCopy WIDTH
1191 // ############################################################################
1192 display->AddTab("WIDTH");
1193
1194 gPad->cd();
1195
1196 hist_width_off->Sumw2();
1197 hist_width_off->Scale(norm);
1198 hist_width_on->SetLineColor(kBlack);
1199 hist_width_on->SetMarkerStyle(21);
1200 hist_width_on->SetMarkerSize(0.7);
1201 hist_width_on->SetMarkerColor(kBlack);
1202 hist_width_off->SetFillColor(46);
1203 hist_width_off->SetLineColor(46);
1204 hist_width_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
1205 hist_width_off->SetMinimum(0.);
1206 hist_width_on->SetTitle("WIDTH distribution");
1207 hist_width_off->SetTitle("WIDTH distribution");
1208
1209 hist_width_on->DrawCopy("E1P");
1210
1211 // move stat box to make them all visible
1212 gPad->Update();
1213 TPaveStats* pavs_on_width = (TPaveStats*) hist_width_on->GetListOfFunctions()->FindObject("stats");
1214 if(pavs_on_width){
1215 shiftx = pavs_on_width->GetX2NDC() - pavs_on_width->GetX1NDC();
1216 pavs_on_width->SetX1NDC(pavs_on_width->GetX1NDC() - shiftx);
1217 pavs_on_width->SetX2NDC(pavs_on_width->GetX2NDC() - shiftx);
1218 }
1219 gPad->Modified();
1220 gPad->Update();
1221
1222 hist_width_off->DrawCopy("HISTSAME");
1223 hist_width_off->DrawCopy("ESAME");
1224 hist_width_on->DrawCopy("E1PSAME");
1225
1226 Double_t chiwidth = ChiSquareNDof((TH1D*)hist_width_on,(TH1D*)hist_width_off);
1227
1228 x_label_pos = widthmax[0]*0.7;
1229 y_label_pos = hist_width_on->GetBinContent(hist_width_on->GetMaximumBin())/2.;
1230
1231 char text_width[256];
1232 sprintf(text_size,"ChiSquare/NDof = %4.2f",chiwidth);
1233
1234 TLatex *twidth = new TLatex(x_label_pos, y_label_pos, text_width);
1235 twidth->SetTextSize(textsize);
1236// twidth->Draw();
1237
1238 gPad->Modified();
1239 gPad->Update();
1240
1241 // ############################################################################
1242 // DrawCopy LENGTH
1243 // ############################################################################
1244 display->AddTab("LENGTH");
1245
1246 gPad->cd();
1247
1248 hist_length_on->Sumw2();
1249 hist_length_off->Sumw2();
1250 hist_length_off->Scale(norm);
1251 hist_length_on->SetLineColor(kBlack);
1252 hist_length_on->SetMarkerStyle(21);
1253 hist_length_on->SetMarkerSize(0.7);
1254 hist_length_on->SetMarkerColor(kBlack);
1255 hist_length_off->SetFillColor(46);
1256 hist_length_off->SetLineColor(46);
1257 hist_length_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
1258 hist_length_off->SetMinimum(0.);
1259 hist_length_on->SetTitle("LENGTH distribution");
1260 hist_length_off->SetTitle("LENGTH distribution");
1261
1262 hist_length_on->DrawCopy("E1P");
1263
1264 // move stat box to make them all visible
1265 gPad->Update();
1266 TPaveStats* pavs_on_length = (TPaveStats*) hist_length_on->GetListOfFunctions()->FindObject("stats");
1267 if(pavs_on_length){
1268 shiftx = pavs_on_length->GetX2NDC() - pavs_on_length->GetX1NDC();
1269 pavs_on_length->SetX1NDC(pavs_on_length->GetX1NDC() - shiftx);
1270 pavs_on_length->SetX2NDC(pavs_on_length->GetX2NDC() - shiftx);
1271 }
1272 gPad->Modified();
1273 gPad->Update();
1274
1275 hist_length_off->DrawCopy("HISTSAME");
1276 hist_length_off->DrawCopy("ESAME");
1277 hist_length_on->DrawCopy("E1PSAME");
1278
1279 Double_t chilength = ChiSquareNDof((TH1D*)hist_length_on,(TH1D*)hist_length_off);
1280
1281 x_label_pos = lengthmax[0]*0.7;
1282 y_label_pos = hist_length_on->GetBinContent(hist_length_on->GetMaximumBin())/2.;
1283
1284 char text_length[256];
1285 sprintf(text_size,"ChiSquare/NDof = %4.2f",chilength);
1286
1287 TLatex *tlength = new TLatex(x_label_pos, y_label_pos, text_length);
1288 tlength->SetTextSize(textsize);
1289// tlength->Draw();
1290
1291 gPad->Modified();
1292 gPad->Update();
1293
1294 // ############################################################################
1295 // DrawCopy normalized ALPHA plot
1296 // ############################################################################
1297 display->AddTab("ALPHA");
1298
1299 gPad->cd();
1300
1301 hist_on_abs->Sumw2();
1302 hist_off_abs->SetStats(0);
1303 hist_off_abs->Sumw2();
1304 hist_off_abs->Scale(norm);
1305 hist_on_abs->SetStats(0); //-> Do NOT show the legend with statistics
1306 hist_on_abs->SetLineColor(kBlack);
1307 hist_on_abs->SetMarkerStyle(21);
1308 //hist_on_abs->SetMarkerSize();
1309 hist_on_abs->SetMarkerColor(kBlack);
1310 hist_on_abs->SetMarkerSize(0.7);
1311 hist_off_abs->SetFillColor(46);
1312 hist_off_abs->SetLineColor(46);
1313 hist_off_abs->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
1314 hist_off_abs->SetMinimum(0.);
1315 hist_on_abs->SetTitle("Alpha plot");
1316 hist_off_abs->SetTitle("Alpha plot");
1317
1318
1319 hist_on_abs->DrawCopy("E1P");
1320 hist_off_abs->DrawCopy("HISTSAME");
1321 hist_off_abs->DrawCopy("ESAME");
1322 hist_on_abs->DrawCopy("E1PSAME");
1323
1324
1325 //draw the LEGEND with excess and significance values in the alpha plot:
1326 char text_Fnorm[256], text_excess[256], text_sign[256];
1327 char text_entries_on[256], text_entries_off[256], text_integral_off[256];
1328 int hist_on_entries = (int) hist_on_abs->GetEntries();
1329 int hist_off_entries = (int) hist_off_abs->GetEntries();
1330 sprintf(text_Fnorm, " F norm = %.3f", norm);
1331 sprintf(text_excess, " Excess = %.3f", excess);
1332 sprintf(text_sign, " Significance = %.3f", sign);
1333 sprintf(text_entries_on, " Entries ON = %d", hist_on_entries);
1334 sprintf(text_entries_off, " Entries OFF = %d", hist_off_entries);
1335 sprintf(text_integral_off," Integral OFF = %d", int_off);
1336
1337 x_label_pos = alphamax[0]*0.7;
1338 y_label_pos = (hist_on_abs->GetBinContent(hist_on_abs->GetMaximumBin()))/1.6; //2.;
1339 Double_t y_label_step = y_label_pos / 8.;
1340
1341 TLatex *t0 = new TLatex(x_label_pos, y_label_pos - y_label_step*0, text_Fnorm);
1342 t0->SetTextSize(textsize);
1343 t0->Draw();
1344 TLatex *t1 = new TLatex(x_label_pos, y_label_pos - y_label_step*1, text_excess);
1345 t1->SetTextSize(textsize);
1346 t1->Draw();
1347 TLatex *t2 = new TLatex(x_label_pos, y_label_pos - y_label_step*2, text_sign);
1348 t2->SetTextSize(textsize);
1349 t2->Draw();
1350 TLatex *t3 = new TLatex(x_label_pos, y_label_pos - y_label_step*3, text_entries_on);
1351 t3->SetTextSize(textsize);
1352 t3->Draw();
1353 TLatex *t4 = new TLatex(x_label_pos, y_label_pos - y_label_step*4, text_entries_off);
1354 t4->SetTextSize(textsize);
1355 t4->Draw();
1356 TLatex *t5 = new TLatex(x_label_pos, y_label_pos - y_label_step*5, text_integral_off);
1357 t5->SetTextSize(textsize);
1358 t5->Draw();
1359
1360
1361 Double_t chialpha = ChiSquareNDof((TH1D*)hist_on_abs,(TH1D*)hist_off_abs);
1362
1363 y_label_pos = (hist_on_abs->GetBinContent(hist_on_abs->GetMaximumBin()))/2.;
1364
1365 char text_alpha[256];
1366 sprintf(text_size,"ChiSquare/NDof = %4.2f",chialpha);
1367
1368 TLatex *talpha = new TLatex(x_label_pos, y_label_pos, text_alpha);
1369 talpha->SetTextSize(textsize);
1370// talpha->Draw();
1371
1372 gPad->Modified();
1373 gPad->Update();
1374
1375 // ############################################################################
1376 // DrawCopy normalized alpha histos for alpha form -90 to 90 deg.
1377 // ############################################################################
1378 display->AddTab("ALPHA +-90");
1379
1380 gPad->cd();
1381
1382 hist_on->Sumw2();
1383 hist_off->SetStats(0);
1384 hist_off->Sumw2();
1385 hist_off->Scale(norm);
1386 hist_off->SetFillColor(46);
1387 hist_off->SetLineColor(46);
1388 hist_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
1389 hist_off->SetMinimum(0.);
1390 hist_on->SetStats(0); //-> Do NOT show the legend with statistics
1391 hist_on->SetLineColor(kBlack);
1392 hist_on->SetMarkerStyle(21);
1393 hist_on->SetMarkerSize(0.7);
1394 hist_on->SetMarkerColor(kBlack);
1395 hist_on->SetTitle("Alpha plot form -90 to 90 deg");
1396 hist_off->SetTitle("Alpha plot form -90 to 90 deg");
1397
1398 hist_on->DrawCopy("E1P");
1399 hist_off->DrawCopy("HISTSAME");
1400 hist_off->DrawCopy("ESAME");
1401 hist_on->DrawCopy("E1PSAME");
1402
1403 Double_t chialpha90 = ChiSquareNDof((TH1D*)hist_on,(TH1D*)hist_off);
1404
1405 x_label_pos = alphamax[0]*0.5;
1406 y_label_pos = hist_on->GetBinContent(hist_on->GetMaximumBin())/2.;
1407
1408 char text_alpha90[256];
1409 sprintf(text_alpha90,"ChiSquare/NDof = %4.2f",chialpha90);
1410
1411 TLatex *talpha90 = new TLatex(x_label_pos, y_label_pos, text_alpha90);
1412 talpha90->SetTextSize(textsize);
1413// talpha90->Draw();
1414
1415 gPad->Update();
1416 gPad->Modified();
1417
1418
1419 cout << "---> ChiSquare/NDof [Size] =\t\t" << chisize << endl;
1420 cout << "---> ChiSquare/NDof [Dist] =\t\t" << chidist << endl;
1421 cout << "---> ChiSquare/NDof [Width] =\t\t" << chiwidth << endl;
1422 cout << "---> ChiSquare/NDof [Length] =\t\t" << chilength << endl;
1423 cout << "---> ChiSquare/NDof [Abs(Alpha)] =\t" << chialpha << endl;
1424 cout << "---> ChiSquare/NDof [Alpha] =\t\t" << chialpha90 << endl;
1425
1426
1427 display->AddTab("SRCPOS ON");
1428 TH2F *hist_srcpos_on = (TH2F*)hSrcPos_on.GetHist();
1429 hist_srcpos_on->DrawCopy("BOX");
1430
1431 display->AddTab("SRCPOS OFF");
1432 TH2F *hist_srcpos_off = (TH2F*)hSrcPos_off.GetHist();
1433 hist_srcpos_off->DrawCopy("BOX");
1434
1435// display->AddTab("EVTNUM ON");
1436// TH1F *hist_evtnum_on = (TH1F*)hDAQEvtNumber_on.GetHist();
1437// hist_evtnum_on->DrawCopy();
1438
1439// display->AddTab("EVTNUM OFF");
1440// TH1F *hist_evtnum_off = (TH1F*)hDAQEvtNumber_off.GetHist();
1441// hist_evtnum_off->DrawCopy();
1442
1443 }
1444
1445 cout << "Done!!" <<endl;
1446
1447}
1448
1449
1450
Note: See TracBrowser for help on using the repository browser.