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

Last change on this file since 4416 was 4408, checked in by jlopez, 21 years ago
*** empty log message ***
File size: 46.3 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 (nbinsnozero>0?chiq/nbinsnozero:0);
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->SetCanvasColor(0);
62 gStyle->SetCanvasBorderMode(0);
63 gStyle->SetPadBorderMode(0);
64 gStyle->SetFrameBorderMode(0);
65 gStyle->SetTimeOffset(-3600);
66
67 // Configuration
68 const Bool_t debug = kFALSE;
69 const Double_t timebin = 600.; //[sec]
70 TString psname = "./20040420_Mrk421.600s.ps";
71
72 //Constanst
73 const Double_t kConvDegToRad = TMath::Pi()/180.;
74 const Double_t kSec = 1e3; //[sec/milisec]
75 const Double_t kDay = 24.*60.*60.; //[Day/sec]
76
77 UInt_t numberTimeBins = 1;
78 TArrayI numberOnEvents(1);
79 TArrayD numberOnEventsAfterCleaning(1);
80 TArrayD numberBkgEventsToNormOn(1);
81 TArrayD timeOn(1);
82
83 TArrayD meanTimeBinOn(1);
84 TArrayD widthTimeBinOn(1);
85
86 TArrayD zenithMinimumOn(1);
87 TArrayD zenithMaximumOn(1);
88
89 TArrayD deadFractionOn(1);
90
91 TObjArray coszenithHistoOn;
92 TObjArray alphaHistoOn;
93 TObjArray srcposHistoOn;
94 TObjArray timediffHistoOn;
95
96 const Int_t numberSizeBins = 4;
97 Double_t sizeBins[numberSizeBins] = {1292., 1897., 2785., 4087.}; //[Photons]
98
99 //cuts
100
101 Double_t widthmin[numberSizeBins] = { 0.06, 0.06, 0.06, 0.06 };
102 Double_t widthmax[numberSizeBins] = { 0.10, 0.10, 0.12, 0.12 };
103 Double_t lengthmin[numberSizeBins] = { 0.10, 0.10, 0.10, 0.10 };
104 Double_t lengthmax[numberSizeBins] = { 0.20, 0.25, 0.26, 0.36 };
105 Double_t distmin[numberSizeBins] = { 0.30, 0.30, 0.30, 0.30 };
106 Double_t distmax[numberSizeBins] = { 1.20, 1.20, 1.20, 1.20 };
107
108 //alpha plot integration for gammas
109 Double_t sigexccmin = 0.;
110 Double_t sigexccmax = 20.;
111 Double_t bkgnormmin = 40.;
112 Double_t bkgnormmax = 80.;
113
114 gStyle->SetOptStat(111111);
115 gStyle->SetOptFit();
116
117 //
118 // Make a loop only for the ON data:
119 //
120
121 MParList plist_on;
122 MTaskList tlist_on;
123 plist_on.AddToList(&tlist_on);
124
125 // ON containers
126 MGeomCamMagic geomcam;
127 MRawRunHeader runheader_on;
128 MRawEvtHeader evtheader_on;
129 MTime time_on;
130 MHillas hillas_on;
131 MHillasSrc hillassrc_on;
132 MSrcPosCam srcpos_on;
133 MReportDrive drive_on;
134
135 plist_on.AddToList(&geomcam);
136 plist_on.AddToList(&runheader_on);
137 plist_on.AddToList(&evtheader_on);
138 plist_on.AddToList(&time_on);
139 plist_on.AddToList(&hillas_on);
140 plist_on.AddToList(&hillassrc_on);
141 plist_on.AddToList(&srcpos_on);
142 plist_on.AddToList(&drive_on);
143
144 Int_t nbins_Size = 60;
145 Double_t min_Size = log10(sizeBins[0])*0.8;
146 Double_t max_Size = log10(1000000)*1.2;
147 TH1F *hSize_on = new TH1F ("hSizeOn","",nbins_Size,min_Size,max_Size);
148
149 Int_t nbins_Width = 20;
150 Double_t min_Width = 0.;
151 Double_t max_Width = widthmax[numberSizeBins-1]*1.2;
152 TH1F *hWidth_on = new TH1F ("hWidthOn","",nbins_Width,min_Width,max_Width);
153
154 Int_t nbins_Length = 20;
155 Double_t min_Length = 0.;
156 Double_t max_Length = lengthmax[numberSizeBins-1]*1.2;
157 TH1F *hLength_on = new TH1F ("hLengthOn","",nbins_Length,min_Length,max_Length);
158
159 Int_t nbins_Dist = 20;
160 Double_t min_Dist = 0.;
161 Double_t max_Dist = distmax[numberSizeBins-1]*1.2;
162 TH1F *hDist_on = new TH1F ("hDistOn","",nbins_Dist,min_Dist,max_Dist);
163
164 Int_t nbins_abs = 18;
165 Double_t minalpha_abs = 0.;
166 Double_t maxalpha_abs =90.;
167 TH1F *hAlpha_on_abs = new TH1F ("hAbsAlphaOn","",nbins_abs,minalpha_abs,maxalpha_abs);
168
169 Int_t nbins = nbins_abs*2;
170 Double_t minalpha = -90.;
171 Double_t maxalpha = 90.;
172 TH1F *hAlpha_on = new TH1F ("hAlphaOn","",nbins,minalpha,maxalpha);
173
174 Int_t nbins_srcpos = 200;
175 Double_t minsrcpos = -200.;
176 Double_t maxsrcpos = 200.; //[mm] //!!! position precition 3mm ~ 0.01 deg
177 TH2F *hSrcPos_on = new TH2F ("hSrcPosOn","",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos);
178
179 //
180 //tasks
181 //
182
183 MReadTree read_on("Parameters", f_on_name);
184 read_on.DisableAutoScheme();
185
186 MSrcPlace srcplace;
187 MHillasSrcCalc csrc_on;
188
189// // prints
190// MPrint pevent("MRawEvtHeader");
191// MPrint phillas("MHillas");
192// MPrint phillassrc("MHillasSrc");
193// MPrint psrcpos("MSrcPosCam");
194
195 //tasklist
196 tlist_on.AddToList(&read_on);
197 tlist_on.AddToList(&srcplace);
198 tlist_on.AddToList(&csrc_on);
199
200 // Create and setup the eventloop
201 MEvtLoop loop_on;
202 loop_on.SetParList(&plist_on);
203 //loop_on.SetDisplay(display);
204
205// MProgressBar bar;
206// loop_on.SetProgressBar(&bar);
207
208// if (!loop_on.Eventloop())
209// return;
210
211 if (!loop_on.PreProcess())
212 return;
213
214 Double_t mjdFirstEventofBin=0;
215 Double_t mjdFirstEvent=0;
216
217 MTime mtimeFirstEventofBin;
218 MTime mtimeFirstEvent;
219
220 Double_t mjdLastEvent=0;
221 MTime mtimeLastEvent=0;
222 UInt_t evtLastEvent;
223 UInt_t runLastEvent=0;
224 MTime mtimeLastEvent;
225
226 Bool_t flag = kFALSE;
227
228 numberOnEventsAfterCleaning[numberTimeBins-1] = 0;
229 zenithMinimumOn[numberTimeBins-1] = 100.;
230 zenithMaximumOn[numberTimeBins-1] = 0.;
231 timeOn[numberTimeBins-1] = 0;
232
233 //create histos needed in the time bins
234
235 TString alphaTitle = Form("%s%02i","hAlphaOn",numberTimeBins-1);
236 TH1F *hAlpha_on_abs_timebin = new TH1F (alphaTitle,"",nbins_abs,minalpha_abs,maxalpha_abs);
237
238 TString srcposTitle = Form("%s%02i","hSrcPosOn",numberTimeBins-1);
239 TH2F *hSrcPos_on_timebin = new TH2F (srcposTitle,"",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos);
240
241 Int_t nbins_coszenith = 200;
242 Double_t mincoszenith = 0.;
243 Double_t maxcoszenith = 1.;
244 TString coszenithTitle = Form("%s%02i","hCosZenithOn",numberTimeBins-1);
245 TH1F *hCosZenith_on_timebin = new TH1F (coszenithTitle,"",nbins_coszenith,mincoszenith,maxcoszenith);
246
247 Int_t nbins_timediff = 2000;
248 Double_t mintimediff = 0.;
249 Double_t maxtimediff = 40.;
250 TString timediffTitle = Form("%s%02i","hTimeDiffOn",numberTimeBins-1);
251 TH1F *hTimeDiff_on_timebin = new TH1F (timediffTitle,"",nbins_timediff,mintimediff,maxtimediff);
252 TF1 *f1 = new TF1("exp","expo",mintimediff,maxtimediff);
253
254 while(tlist_on.Process())
255 {
256 // Compute live time
257 UInt_t run = runheader_on.GetRunNumber();
258 UInt_t evt = evtheader_on.GetDAQEvtNumber();
259 Double_t mjd = time_on.GetMjd();
260
261 numberOnEventsAfterCleaning[numberTimeBins-1]++;
262
263 if (mjd == 0)
264 {
265
266 if (debug)
267 {
268 cout << "MTime::GetTime() == 0" << endl;
269 cout.precision(15);
270 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
271 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl;
272 mtimeLastEvent.Print();
273 time_on.Print();
274 }
275
276 flag = kTRUE;
277 }
278 else if (mjd-mjdLastEvent < 0 && flag)
279 {
280
281 if (debug)
282 {
283 cout << "mjd-mjdLastEvent < 0" << endl;
284 cout.precision(15);
285 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
286 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl;
287 mtimeLastEvent.Print();
288 time_on.Print();
289 }
290
291 flag = kTRUE;
292 }
293 else if (mjd-mjdLastEvent == 0 && flag)
294 {
295
296 if (debug)
297 {
298 cout << "mjd-mjdLastEvent == 0" << endl;
299 cout.precision(15);
300 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
301 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl;
302 mtimeLastEvent.Print();
303 time_on.Print();
304 }
305
306 flag = kTRUE;
307 }
308 else if (evt-evtLastEvent <= 0)
309 {
310
311 if (debug)
312 {
313 cout << "evt-evtLastEvent <= 0" << endl;
314 cout.precision(15);
315 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
316 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl;
317 mtimeLastEvent.Print();
318 time_on.Print();
319 }
320
321 flag = kTRUE;
322 }
323
324 else if ((Int_t)mjd == mjd)
325 {
326
327 if (debug)
328 {
329 cout << "(Int_t)mjd == mjd" << endl;
330 cout.precision(15);
331 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
332 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl;
333 mtimeLastEvent.Print();
334 time_on.Print();
335 }
336
337 flag = kTRUE;
338 }
339 else
340 {
341
342 if (flag && debug)
343 {
344
345 cout << "flag = TRUE" << endl;
346 cout.precision(15);
347 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
348 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl;
349 mtimeLastEvent.Print();
350 time_on.Print();
351
352 flag = kFALSE;
353
354 }
355
356 if ( run != runLastEvent )
357 {
358 if ( mjdLastEvent != 0 )
359 {
360 timeOn[numberTimeBins-1] += (mjdLastEvent-mjdFirstEvent)*kDay;
361
362 cout.precision(10);
363 cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl;
364 if (flag && debug)
365 {
366 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
367 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl;
368 mtimeLastEvent.Print();
369 time_on.Print();
370 }
371
372 }
373
374 mjdFirstEvent = mjd;
375 }
376 else if (mjdFirstEventofBin != 0)
377 hTimeDiff_on_timebin->Fill((mjd-mjdLastEvent)*kDay*kSec);
378
379 if (mjdFirstEventofBin == 0)
380 {
381 mjdFirstEvent = mjd;
382 mjdFirstEventofBin = mjd;
383 }
384
385 evtLastEvent = evt;
386 runLastEvent = run;
387 mjdLastEvent = mjd;
388 mtimeLastEvent = time_on;
389
390 Double_t size = hillas_on.GetSize();
391 Double_t width = hillas_on.GetWidth()*geomcam.GetConvMm2Deg();
392 Double_t length = hillas_on.GetLength()*geomcam.GetConvMm2Deg();
393 Double_t meanx = hillas_on.GetMeanX()*geomcam.GetConvMm2Deg();
394 Double_t meany = hillas_on.GetMeanY()*geomcam.GetConvMm2Deg();
395 Double_t centerdist = TMath::Sqrt(meanx*meanx + meany*meany);
396 Double_t dist = hillassrc_on.GetDist()*geomcam.GetConvMm2Deg();
397 Double_t alpha = hillassrc_on.GetAlpha();
398 Double_t srcposx = srcpos_on.GetX();
399 Double_t srcposy = srcpos_on.GetY();
400 Double_t zenith = drive_on.GetNominalZd();
401
402
403 Int_t sizebin = GetBin(size,numberSizeBins,sizeBins);
404
405 if (sizebin >= 0)
406 {
407 if (width > widthmin[sizebin] && width < widthmax[sizebin])
408 {
409 if (length > lengthmin[sizebin] && length < lengthmax[sizebin])
410 {
411 if (dist > distmin[sizebin] && centerdist < distmax[sizebin])
412 {
413
414 //General histos
415 hSize_on->Fill(log10(size));
416 hWidth_on->Fill(width);
417 hLength_on->Fill(length);
418 hDist_on->Fill(dist);
419 hAlpha_on_abs->Fill(TMath::Abs(alpha));
420 hAlpha_on->Fill(alpha);
421 hSrcPos_on->Fill(srcposx,srcposy);
422
423 // Time bin histos
424 hAlpha_on_abs_timebin->Fill(TMath::Abs(alpha));
425 hSrcPos_on_timebin->Fill(srcposx,srcposy);
426 hCosZenith_on_timebin->Fill(TMath::Cos(zenith*kConvDegToRad));
427
428 if (zenith != 0 && zenith < zenithMinimumOn[numberTimeBins-1])
429 zenithMinimumOn[numberTimeBins-1] = zenith;
430 if (zenith>zenithMaximumOn[numberTimeBins-1])
431 zenithMaximumOn[numberTimeBins-1] = zenith;
432
433 }
434 }
435 }
436 }
437
438
439 // Check if you are overload the maxim time bin
440 if ((mjd-mjdFirstEventofBin)*kDay >= timebin)
441 {
442 //Compute the time on
443 timeOn[numberTimeBins-1] += (mjdLastEvent-mjdFirstEvent)*kDay;
444 widthTimeBinOn[numberTimeBins-1] = (mjdLastEvent-mjdFirstEventofBin)/2;
445 meanTimeBinOn[numberTimeBins-1] = mjdFirstEventofBin + widthTimeBinOn[numberTimeBins-1];
446
447 cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl;
448 cout << "mjd " << mjd << " mjdFirstEventofBin " << mjdFirstEventofBin << " widthTimeBinOn " << widthTimeBinOn[numberTimeBins-1] << " meanTimeBinOn " << meanTimeBinOn[numberTimeBins-1] << ' ' << mjdFirstEventofBin + widthTimeBinOn[numberTimeBins-1] << endl;
449
450 //Compute the number of on events
451 numberOnEvents[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
452 numberBkgEventsToNormOn[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
453 //Initialize variables
454 mjdFirstEvent = mjd;
455 mjdFirstEventofBin = mjd;
456
457 hTimeDiff_on_timebin->Fit("exp","RQ0");
458 deadFractionOn[numberTimeBins-1] = (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/(TMath::Abs(f1->GetParameter(1))*hTimeDiff_on_timebin->GetEntries());
459 cout << "1-> Exp(" << f1->GetParameter(0) << " + " << f1->GetParameter(1) << "*x) +- [" << f1->GetParError(0) << ' ' << f1->GetParError(1) << "]" << endl;
460 cout << "Calc entries " << (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/(TMath::Abs(f1->GetParameter(1))) << " +- Nt*(" << TMath::Sqrt( f1->GetParameter(0)*f1->GetParError(0)*f1->GetParameter(0)*f1->GetParError(0) + (f1->GetParError(1)*f1->GetParError(1))/(f1->GetParameter(1)*f1->GetParameter(1))) << ") meas entries " << hTimeDiff_on_timebin->GetEntries() << " dead fraction " << deadFractionOn[numberTimeBins-1] << endl;
461
462 alphaHistoOn.AddLast(hAlpha_on_abs_timebin);
463 srcposHistoOn.AddLast(hSrcPos_on_timebin);
464 coszenithHistoOn.AddLast(hCosZenith_on_timebin);
465 timediffHistoOn.AddLast(hTimeDiff_on_timebin);
466
467// cout << "hAlpha_on_abs_timebin " << hAlpha_on_abs_timebin << " alphaHistoOn [" << numberTimeBins-1 << "] " << alphaHistoOn[numberTimeBins-1] << endl;
468
469 //Increase the number of time bins and all needed arrays
470 numberTimeBins++;
471
472 timeOn.Set(numberTimeBins);
473 numberOnEvents.Set(numberTimeBins);
474 numberBkgEventsToNormOn.Set(numberTimeBins);
475 widthTimeBinOn.Set(numberTimeBins);
476 meanTimeBinOn.Set(numberTimeBins);
477 zenithMinimumOn.Set(numberTimeBins);
478 zenithMaximumOn.Set(numberTimeBins);
479 deadFractionOn.Set(numberTimeBins);
480 numberOnEventsAfterCleaning.Set(numberTimeBins);
481
482 timeOn[numberTimeBins-1] = 0.;
483 zenithMinimumOn[numberTimeBins-1] = 100.;
484 zenithMaximumOn[numberTimeBins-1] = 0.;
485 numberOnEventsAfterCleaning[numberTimeBins-1] = 0;
486
487 alphaTitle = Form("%s%02i","hAlphaOn",numberTimeBins-1);
488 hAlpha_on_abs_timebin = new TH1F (alphaTitle,"",nbins_abs,minalpha_abs,maxalpha_abs);
489
490 srcposTitle = Form("%s%02i","hSrcPosOn",numberTimeBins-1);
491 hSrcPos_on_timebin = new TH2F (srcposTitle,"",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos);
492
493 coszenithTitle = Form("%s%02i","hCosZenithOn",numberTimeBins-1);
494 hCosZenith_on_timebin = new TH1F (coszenithTitle,"",nbins_coszenith,mincoszenith,maxcoszenith);
495
496 timediffTitle = Form("%s%02i","hTimeDiffOn",numberTimeBins-1);
497 hTimeDiff_on_timebin = new TH1F (timediffTitle,"",nbins_timediff,mintimediff,maxtimediff);
498
499 flag = kTRUE;
500 }
501
502 }
503
504 }
505
506
507 //Compute the time on for last time bin
508 timeOn[numberTimeBins-1] += (mjdLastEvent-mjdFirstEvent)*kDay;
509 widthTimeBinOn[numberTimeBins-1] = (mjd-mjdFirstEventofBin)/2;
510 meanTimeBinOn[numberTimeBins-1] = mjdFirstEventofBin + widthTimeBinOn[numberTimeBins-1];
511 //Compute the number of on events for the last time bin
512 numberOnEvents[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
513 numberBkgEventsToNormOn[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
514
515 hTimeDiff_on_timebin->Fit("exp","RQ0");
516 deadFractionOn[numberTimeBins-1] = (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/(TMath::Abs(f1->GetParameter(1))*hTimeDiff_on_timebin->GetEntries());
517
518 cout.precision(5);
519 cout << "2-> Exp(" << f1->GetParameter(0) << " + " << f1->GetParameter(1) << "*x) +- [" << f1->GetParError(0) << ' ' << f1->GetParError(1) << "]" << endl;
520 cout << "Calc entries " << (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/TMath::Abs(f1->GetParameter(1)) << " +- Nt*(" << TMath::Sqrt( f1->GetParameter(0)*f1->GetParError(0)*f1->GetParameter(0)*f1->GetParError(0) + (f1->GetParError(1)*f1->GetParError(1))/(f1->GetParameter(1)*f1->GetParameter(1))) << ") meas entries " << hTimeDiff_on_timebin->GetEntries() << " dead fraction " << deadFractionOn[numberTimeBins-1] << endl;
521
522 cout.precision(10);
523 cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl;
524 cout << "mjd " << mjd << " mjdFirstEventofBin " << mjdFirstEventofBin << " meanTimeBinOn " << meanTimeBinOn[numberTimeBins-1] << " widthTimeBinOn " << widthTimeBinOn[numberTimeBins-1] << endl;
525
526 alphaHistoOn.AddLast(hAlpha_on_abs_timebin);
527 srcposHistoOn.AddLast(hSrcPos_on_timebin);
528 coszenithHistoOn.AddLast(hCosZenith_on_timebin);
529 timediffHistoOn.AddLast(hTimeDiff_on_timebin);
530
531 // cout << "hAlpha_on_abs_timebin " << hAlpha_on_abs_timebin << " alphaHistoOn [" << numberTimeBins-1 << "] " << alphaHistoOn[numberTimeBins-1] << endl;
532
533 loop_on.PostProcess();
534
535 tlist_on.PrintStatistics();
536
537 for (UInt_t bin=0; bin<numberTimeBins; bin++)
538 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] << " numberOnEventsAfterCleaning " << numberOnEventsAfterCleaning[bin] << " deadFractionOn " << deadFractionOn[bin] << endl;
539
540
541
542
543 //-----------------------OFF------------------------
544
545 TObjArray alphaHistoOff;
546 TObjArray srcposHistoOff;
547
548 TArrayD timeOff(numberTimeBins);
549
550 TArrayI numberOffEvents(numberTimeBins);
551
552 TArrayD numberBkgEventsToNormOff(numberTimeBins);
553 TArrayD onOffNormFactor(numberTimeBins);
554
555 TArrayD numberExcessEvents(numberTimeBins);
556 TArrayD errorNumberExcessEvents(numberTimeBins);
557
558 TArrayD numberExcessEventsPerSec(numberTimeBins);
559 TArrayD errorNumberExcessEventsPerSec(numberTimeBins);
560
561 TArrayD numberExcessEventsPerMin(numberTimeBins);
562 TArrayD errorNumberExcessEventsPerMin(numberTimeBins);
563
564 timeOff.Set(numberTimeBins);
565
566 TH1F* hAlpha_off_abs_timebin;
567 TH2F* hSrcPos_off_timebin;
568 for (UInt_t bin=0; bin<numberTimeBins; bin++)
569 {
570 alphaTitle = Form("%s%02i","hAlphaOff",bin);
571 hAlpha_off_abs_timebin = new TH1F (alphaTitle,"",nbins_abs,minalpha_abs,maxalpha_abs);;
572 alphaHistoOff.AddLast(hAlpha_off_abs_timebin);
573 // cout << "hAlpha_off_abs_timebin " << hAlpha_off_abs_timebin << " alphaHistoOff [" << bin << "] " << alphaHistoOff[bin] << endl;
574
575 srcposTitle = Form("%s%02i","hSrcPosOff",bin);
576 hSrcPos_off_timebin = new TH2F (srcposTitle,"",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos);
577 srcposHistoOff.AddLast(hSrcPos_off_timebin);
578 // cout << "hSrcPos_off_timebin " << hSrcPos_off_timebin << " srcposHistoOff [" << bin << "] " << srcposHistoOff[bin] << endl;
579 }
580
581 //
582 // Make a loop only for the OFF data:
583 //
584
585 MParList plist_off;
586 MTaskList tlist_off;
587 plist_off.AddToList(&tlist_off);
588
589 MRawRunHeader runheader_off;
590 MRawEvtHeader evtheader_off;
591 MTime time_off;
592 MHillas hillas_off;
593 MHillasSrc hillassrc_off;
594 MSrcPosCam srcpos_off;
595 MReportDrive drive_off;
596
597 plist_off.AddToList(&geomcam);
598 plist_off.AddToList(&runheader_off);
599 plist_off.AddToList(&evtheader_off);
600 plist_off.AddToList(&time_off);
601 plist_off.AddToList(&hillas_off);
602 plist_off.AddToList(&hillassrc_off);
603 plist_off.AddToList(&srcpos_off);
604 plist_off.AddToList(&drive_off);
605
606 TH1F *hSize_off = new TH1F ("hSizeOff","",nbins_Size,min_Size,max_Size);
607 TH1F *hWidth_off = new TH1F ("hWidthOff","",nbins_Width,min_Width,max_Width);
608 TH1F *hLength_off = new TH1F ("hLengthOff","",nbins_Length,min_Length,max_Length);
609 TH1F *hDist_off = new TH1F ("hDistOff","",nbins_Dist,min_Dist,max_Dist);
610 TH1F *hAlpha_off_abs = new TH1F ("hAbsAlphaOff","",nbins_abs,minalpha_abs,maxalpha_abs);
611 TH1F *hAlpha_off = new TH1F ("hAlphaOff","",nbins,minalpha,maxalpha);
612 TH2F *hSrcPos_off = new TH2F ("hSrcPosOff","",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos);
613
614 //tasks
615 MReadTree read_off("Parameters", f_off_name);
616 read_off.DisableAutoScheme();
617
618 // taks to set the src position in the off data for the time bins
619 // --------------------------------------------------------------
620 MSrcPlace srcplace_timebin;
621 srcplace_timebin.SetCreateHisto(kFALSE);
622 srcplace_timebin.SetMode(MSrcPlace::kOff);
623 // --------------------------------------------------------------
624
625 // srcplace.SetMode(MSrcPlace::kOff);
626
627 MHillasSrcCalc csrc_off;
628
629 //tasklist
630 tlist_off.AddToList(&read_off);
631 tlist_off.AddToList(&srcplace_timebin); // This is just to run the preprocess correctely
632 tlist_off.AddToList(&csrc_off); // This is just to run the preprocess correctely
633
634 // Create and setup the eventloop
635 MEvtLoop loop_off;
636 loop_off.SetParList(&plist_off);
637 //loop_off.SetDisplay(display);
638
639// MProgressBar bar_off;
640// loop_off.SetProgressBar(&bar_off);
641
642// if (!loop_off.Eventloop(numEntries))
643// return;
644
645 if (!loop_off.PreProcess())
646 return;
647
648 while(tlist_off.Process())
649 {
650
651 //First read the variables source(i.e. time bin) independent
652 Double_t zenith = drive_off.GetNominalZd();
653
654 Double_t size = hillas_off.GetSize();
655 Double_t width = hillas_off.GetWidth()*geomcam.GetConvMm2Deg();
656 Double_t length = hillas_off.GetLength()*geomcam.GetConvMm2Deg();
657 Double_t meanx = hillas_off.GetMeanX()*geomcam.GetConvMm2Deg();
658 Double_t meany = hillas_off.GetMeanY()*geomcam.GetConvMm2Deg();
659 Double_t centerdist = TMath::Sqrt(meanx*meanx + meany*meany);
660
661 Int_t sizebin = GetBin(size,numberSizeBins,sizeBins);
662
663 if (sizebin >= 0)
664 {
665 if (width > widthmin[sizebin] && width < widthmax[sizebin])
666 {
667 if (length > lengthmin[sizebin] && length < lengthmax[sizebin])
668 {
669 if (centerdist < distmax[sizebin])
670 {
671 //General histos
672
673 srcplace_timebin.SetPositionHisto(hSrcPos_on);
674 srcplace_timebin.CallProcess();
675 csrc_off.CallProcess();
676
677 Double_t dist = hillassrc_off.GetDist()*geomcam.GetConvMm2Deg();
678 Double_t alpha = hillassrc_off.GetAlpha();
679 Double_t srcposx = srcpos_off.GetX();
680 Double_t srcposy = srcpos_off.GetY();
681
682 if (dist > distmin[sizebin] )
683 {
684 hSize_off->Fill(log10(size));
685 hWidth_off->Fill(width);
686 hLength_off->Fill(length);
687
688 hDist_off->Fill(dist);
689 hAlpha_off_abs->Fill(TMath::Abs(alpha));
690 hAlpha_off->Fill(alpha);
691 hSrcPos_off->Fill(srcposx,srcposy);
692 }
693
694 // Time bin histos
695 for (UInt_t bin=0; bin<numberTimeBins; bin++)
696 {
697 srcplace_timebin.SetPositionHisto((TH2F*)srcposHistoOn[bin]);
698 srcplace_timebin.CallProcess();
699 csrc_off.CallProcess();
700
701 dist = hillassrc_off.GetDist()*geomcam.GetConvMm2Deg();
702 alpha = hillassrc_off.GetAlpha();
703 srcposx = srcpos_off.GetX();
704 srcposy = srcpos_off.GetY();
705
706 if (dist > distmin[sizebin])
707 {
708 ((TH1F*)alphaHistoOff[bin])->Fill(TMath::Abs(alpha));
709 ((TH2F*)srcposHistoOff[bin])->Fill(srcposx,srcposy);
710 }
711 }
712 }
713 }
714 }
715 }
716
717 // if (zenith > zenithMaximumOn[numberTimeBins-1])
718 // {
719 // cout << "Exit loop because we arrive to zenith = " << zenith << endl;
720 // break;
721 // }
722 }
723
724 loop_off.PostProcess();
725
726 tlist_off.PrintStatistics();
727
728 TArrayD meanTimeBinOnInSec(numberTimeBins);
729 TArrayD widthTimeBinOnInSec(numberTimeBins);
730
731 TArrayD meanTriggerRateOn(numberTimeBins);
732 TArrayD errorMeanTriggerRateOn(numberTimeBins);
733
734 for (UInt_t bin=0; bin<numberTimeBins; bin++)
735 {
736
737 meanTriggerRateOn[bin] = numberOnEventsAfterCleaning[bin]/timeOn[bin];
738 errorMeanTriggerRateOn[bin] = TMath::Sqrt(numberOnEventsAfterCleaning[bin])/timeOn[bin];
739
740 meanTimeBinOnInSec[bin] = (meanTimeBinOn[bin]-(Int_t)meanTimeBinOn[bin]);
741 if (meanTimeBinOnInSec[bin] > 0.5)
742 meanTimeBinOnInSec[bin] = meanTimeBinOnInSec[bin] - 1 ;
743 meanTimeBinOnInSec[bin] *= kDay;
744
745 widthTimeBinOnInSec[bin] = widthTimeBinOn[bin]*kDay;
746
747 numberOffEvents[bin] = (Double_t) ((TH1F*)alphaHistoOff[bin])->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
748 numberBkgEventsToNormOff[bin] = (Double_t) ((TH1F*)alphaHistoOff[bin])->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
749 if (numberOffEvents[bin] != 0 && numberBkgEventsToNormOff[bin] != 0)
750 {
751 onOffNormFactor[bin] = numberBkgEventsToNormOn[bin]/numberBkgEventsToNormOff[bin];
752 numberExcessEvents[bin] = numberOnEvents[bin] - onOffNormFactor[bin]*numberOffEvents[bin];
753 errorNumberExcessEvents[bin] = TMath::Sqrt(numberOnEvents[bin] + onOffNormFactor[bin]*onOffNormFactor[bin]*numberOffEvents[bin]);
754 numberExcessEventsPerSec[bin] = numberExcessEvents[bin]/timeOn[bin]*(deadFractionOn[bin]>1.?deadFractionOn[bin]:1.);
755 errorNumberExcessEventsPerSec[bin] = errorNumberExcessEvents[bin]/timeOn[bin];
756 numberExcessEventsPerMin[bin] = 60.*numberExcessEvents[bin]/timeOn[bin]*(deadFractionOn[bin]>1.?deadFractionOn[bin]:1.);
757 errorNumberExcessEventsPerMin[bin] = 60.*errorNumberExcessEvents[bin]/timeOn[bin];
758 }
759 }
760
761 for (UInt_t bin=0; bin<numberTimeBins; bin++)
762 {
763 cout.precision(5);
764 cout << bin << " timeOn " << timeOn[bin] << " mean-width time " << meanTimeBinOn[bin] << "-" << widthTimeBinOn[bin] << endl;
765 cout << " numberOnEvents\t " << numberOnEvents[bin] << endl;
766 cout << " numberOffEvents\t " << numberOffEvents[bin] << endl;
767 cout << " numberBkgEventsToNormOn\t " << numberBkgEventsToNormOn[bin] << endl;
768 cout << " numberBkgEventsToNormOff\t " << numberBkgEventsToNormOff[bin] << endl;
769 cout << " onOffNormFactor\t " << onOffNormFactor[bin] << endl;
770 cout << " numberExcessEvents\t " << numberExcessEvents[bin] << " +- " << errorNumberExcessEvents[bin] << endl;
771 cout << " deadFraction\t" << deadFractionOn[bin] << endl;
772 cout << " Excess/Sec\t " << numberExcessEventsPerSec[bin] << " +- " << errorNumberExcessEventsPerSec[bin] << endl;
773 cout << " Trigger Rate\t " << meanTriggerRateOn[bin] << " +- " << errorMeanTriggerRateOn[bin] << endl;
774 }
775
776
777
778 TString openpsname = psname + "(";
779 TString closepsname = psname + ")";
780
781 TCanvas *c0 = new TCanvas;
782 c0->cd(1);
783 TGraphErrors* lightcurvegraph = new TGraphErrors(numberTimeBins,meanTimeBinOnInSec.GetArray(),numberExcessEventsPerMin.GetArray(),widthTimeBinOnInSec.GetArray(),errorNumberExcessEventsPerMin.GetArray());
784 lightcurvegraph->SetTitle("LightCurve");
785 lightcurvegraph->SetMinimum(0.);
786 lightcurvegraph->SetMarkerStyle(21);
787 lightcurvegraph->SetMarkerSize(0.03);
788 lightcurvegraph->Draw("AP");
789 lightcurvegraph->GetYaxis()->SetTitle("Excess/min");
790 lightcurvegraph->GetXaxis()->SetTitle("UTC Time");
791 lightcurvegraph->GetXaxis()->SetTimeDisplay(1);
792 c0->Print(openpsname);
793
794 TCanvas *c00 = new TCanvas;
795 c00->cd(1);
796 TGraphErrors* cosmicrategraph = new TGraphErrors(numberTimeBins,meanTimeBinOnInSec.GetArray(),meanTriggerRateOn.GetArray(),widthTimeBinOnInSec.GetArray(),errorMeanTriggerRateOn.GetArray());
797 cosmicrategraph->SetTitle("Cosmic Rate");
798 cosmicrategraph->SetMinimum(0.);
799 cosmicrategraph->SetMarkerStyle(21);
800 cosmicrategraph->SetMarkerSize(0.03);
801 cosmicrategraph->Draw("AP");
802 cosmicrategraph->GetYaxis()->SetTitle("[Hz]");
803 cosmicrategraph->GetXaxis()->SetTitle("UTC Time");
804 cosmicrategraph->GetXaxis()->SetTimeDisplay(1);
805 c00->Print(psname);
806
807 TCanvas** c = new TCanvas*[numberTimeBins];
808
809 //Compute the maximum of all hAlphaOn histograms
810 Float_t maxAlphaHistoHeight = 0;
811
812 for (UInt_t bin=0; bin<numberTimeBins-1; bin++)
813 {
814 for (UInt_t i=1; i<=nbins_abs; i++)
815 if (((TH1F*)alphaHistoOn[bin])->GetBinContent(i) > maxAlphaHistoHeight)
816 maxAlphaHistoHeight = ((TH1F*)alphaHistoOn[bin])->GetBinContent(i);
817 }
818
819
820 for (UInt_t bin=0; bin<numberTimeBins-1; bin++)
821 {
822 c[bin] = new TCanvas;
823 c[bin]->cd(1);
824
825 ((TH1F*)alphaHistoOn[bin])->Sumw2();
826 ((TH1F*)alphaHistoOff[bin])->SetStats(0);
827 ((TH1F*)alphaHistoOff[bin])->Sumw2();
828 ((TH1F*)alphaHistoOff[bin])->Scale(onOffNormFactor[bin]);
829 ((TH1F*)alphaHistoOn[bin])->SetStats(0); //-> Do NOT show the legend with statistics
830 ((TH1F*)alphaHistoOn[bin])->SetLineColor(kBlack);
831 ((TH1F*)alphaHistoOn[bin])->SetMarkerStyle(21);
832 ((TH1F*)alphaHistoOn[bin])->SetMarkerColor(kBlack);
833 ((TH1F*)alphaHistoOn[bin])->SetMarkerSize(0.7);
834 ((TH1F*)alphaHistoOff[bin])->SetFillColor(46);
835 ((TH1F*)alphaHistoOff[bin])->SetLineColor(46);
836 ((TH1F*)alphaHistoOff[bin])->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
837 ((TH1F*)alphaHistoOn[bin])->SetMaximum(maxAlphaHistoHeight*1.2);
838 ((TH1F*)alphaHistoOff[bin])->SetMinimum(0.);
839 alphaTitle = Form("%s%02i","hAlphaOnOff",bin);
840 ((TH1F*)alphaHistoOn[bin])->SetTitle(alphaTitle);
841
842
843 ((TH1F*)alphaHistoOn[bin])->DrawCopy("E1P");
844 ((TH1F*)alphaHistoOff[bin])->DrawCopy("HISTSAME");
845 ((TH1F*)alphaHistoOff[bin])->DrawCopy("ESAME");
846 ((TH1F*)alphaHistoOn[bin])->DrawCopy("E1PSAME");
847
848 c[bin]->Print(psname);
849 }
850
851 c[numberTimeBins-1] = new TCanvas;
852 c[numberTimeBins-1]->cd(1);
853
854 ((TH1F*)alphaHistoOn[numberTimeBins-1])->Sumw2();
855 ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetStats(0);
856 ((TH1F*)alphaHistoOff[numberTimeBins-1])->Sumw2();
857 ((TH1F*)alphaHistoOff[numberTimeBins-1])->Scale(onOffNormFactor[numberTimeBins-1]);
858 ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetStats(0); //-> Do NOT show the legend with statistics
859 ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetLineColor(kBlack);
860 ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetMarkerStyle(21);
861 ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetMarkerColor(kBlack);
862 ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetMarkerSize(0.7);
863 ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetFillColor(46);
864 ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetLineColor(46);
865 ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
866 ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetMaximum(maxAlphaHistoHeight*1.2);
867 ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetMinimum(0.);
868 alphaTitle = Form("%s%02i","hAlphaOnOff",numberTimeBins-1);
869 ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetTitle(alphaTitle);
870
871
872 ((TH1F*)alphaHistoOn[numberTimeBins-1])->DrawCopy("E1P");
873 ((TH1F*)alphaHistoOff[numberTimeBins-1])->DrawCopy("HISTSAME");
874 ((TH1F*)alphaHistoOff[numberTimeBins-1])->DrawCopy("ESAME");
875 ((TH1F*)alphaHistoOn[numberTimeBins-1])->DrawCopy("E1PSAME");
876
877 c[numberTimeBins-1]->Print(psname);
878
879 // TString rootname = psname.ReplaceAll(".ps",".root");
880 TString rootname = "./prueba.root";
881 TFile input(rootname, "RECREATE");
882
883 for (UInt_t bin=0; bin<numberTimeBins; bin++)
884 {
885 ((TH1F*)alphaHistoOn[bin])->Write();
886 ((TH2F*)srcposHistoOn[bin])->Write();
887 ((TH1F*)coszenithHistoOn[bin])->Write();
888 ((TH1F*)timediffHistoOn[bin])->Write();
889 ((TH1F*)alphaHistoOff[bin])->Write();
890 ((TH2F*)srcposHistoOff[bin])->Write();
891 }
892
893 input.Close();
894
895 // ############################################################################
896 // Calculate significance and excess:
897 // ############################################################################
898
899 Double_t norm_on_abs = (Double_t) hAlpha_on_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
900 Double_t exces_on_abs = (Double_t) hAlpha_on_abs->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
901 Double_t norm_off_abs = (Double_t) hAlpha_off_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
902 Double_t exces_off_abs = (Double_t) hAlpha_off_abs->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
903 Double_t norm = norm_on_abs/norm_off_abs;
904
905 char text_tit_alpha[256];
906 sprintf(text_tit_alpha, " Alpha Plot On and Off ");
907 hAlpha_off_abs->SetTitle(text_tit_alpha);
908 hAlpha_on_abs->SetTitle(text_tit_alpha);
909
910 Double_t excess = exces_on_abs - exces_off_abs*norm;
911 Double_t sign = excess / sqrt( exces_on_abs + norm*norm*exces_off_abs );
912 Double_t int_off = (Double_t) hAlpha_off_abs->Integral(1, 18);
913 int hAlpha_on_entries = (int) hAlpha_on_abs->GetEntries();
914 int hAlpha_off_entries = (int) hAlpha_off_abs->GetEntries();
915
916 cout << "---> Normalization F factor =\t" << norm <<endl;
917 cout << "---> Excess =\t\t\t" << excess <<endl;
918 cout << "---> Significancia =\t\t" << sign <<endl;
919 cout << "---> entries on =\t\t" << hAlpha_on_entries <<endl;
920 cout << "---> entries off =\t\t" << hAlpha_off_entries <<endl;
921 cout << "---> integral off =\t\t" << int_off <<endl;
922
923 Double_t shiftx;
924
925 // ############################################################################
926 // Draw SIZE
927 // ############################################################################
928 TCanvas *c1 = new TCanvas;
929
930 gPad->cd();
931
932 gPad->SetLogy();
933 hSize_on->Sumw2();
934 hSize_off->Sumw2();
935 hSize_off->Scale(norm);
936 hSize_on->SetLineColor(kBlack);
937 hSize_on->SetMarkerStyle(21);
938 hSize_on->SetMarkerSize(0.7);
939 hSize_on->SetMarkerColor(kBlack);
940 hSize_off->SetFillColor(46);
941 hSize_off->SetLineColor(46);
942 hSize_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
943 hSize_off->SetMinimum(0.1);
944 hSize_on->SetMinimum(0.1);
945 hSize_on->SetTitle("SIZE distribution");
946 hSize_off->SetTitle("SIZE distribution");
947
948 hSize_on->DrawCopy("E1P");
949
950 // move stat box to make them all visible
951 gPad->Update();
952 TPaveStats* pavs_on_size = (TPaveStats*) hSize_on->GetListOfFunctions()->FindObject("stats");
953 if(pavs_on_size){
954 shiftx = pavs_on_size->GetX2NDC() - pavs_on_size->GetX1NDC();
955 pavs_on_size->SetX1NDC(pavs_on_size->GetX1NDC() - shiftx);
956 pavs_on_size->SetX2NDC(pavs_on_size->GetX2NDC() - shiftx);
957 }
958 gPad->Modified();
959 gPad->Update();
960
961 hSize_off->DrawCopy("HISTSAME");
962 hSize_off->DrawCopy("ESAME");
963
964 gPad->Modified();
965 gPad->Update();
966
967 Double_t chisize = ChiSquareNDof((TH1D*)hSize_on,(TH1D*)hSize_off);
968
969 Double_t x_label_pos = log10(1000000)*0.7;
970 Double_t y_label_pos = log10((hSize_on->GetBinContent(hSize_on->GetMaximumBin()))/2.);
971 Double_t textsize = 0.03;
972
973 char text_size[256];
974 sprintf(text_size,"ChiSquare/NDof = %4.2f",chisize);
975
976 TLatex *tsize = new TLatex(x_label_pos, y_label_pos, text_size);
977 tsize->SetTextSize(textsize);
978// tsize->Draw();
979
980 gPad->Modified();
981 gPad->Update();
982
983 c1->Print(psname);
984
985 // ############################################################################
986 // DrawCopy DIST
987 // ############################################################################
988 TCanvas *c2 = new TCanvas;
989
990 gPad->cd();
991
992 hDist_on->Sumw2();
993 hDist_off->Sumw2();
994 hDist_off->Scale(norm);
995 hDist_on->SetLineColor(kBlack);
996 hDist_on->SetMarkerStyle(21);
997 hDist_on->SetMarkerSize(0.7);
998 hDist_on->SetMarkerColor(kBlack);
999 hDist_off->SetFillColor(46);
1000 hDist_off->SetLineColor(46);
1001 hDist_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
1002 hDist_off->SetMinimum(0.);
1003 hDist_on->SetTitle("DIST distribution");
1004 hDist_off->SetTitle("DIST distribution");
1005
1006 hDist_on->DrawCopy("E1P");
1007
1008 // move stat box to make them all visible
1009 gPad->Update();
1010 TPaveStats* pavs_on_dist = (TPaveStats*) hDist_on->GetListOfFunctions()->FindObject("stats");
1011 if(pavs_on_dist){
1012 shiftx = pavs_on_dist->GetX2NDC() - pavs_on_dist->GetX1NDC();
1013 pavs_on_dist->SetX1NDC(pavs_on_dist->GetX1NDC() - shiftx);
1014 pavs_on_dist->SetX2NDC(pavs_on_dist->GetX2NDC() - shiftx);
1015 }
1016 gPad->Modified();
1017 gPad->Update();
1018
1019 hDist_off->DrawCopy("HISTSAME");
1020 hDist_off->DrawCopy("ESAME");
1021 hDist_on->DrawCopy("E1PSAME");
1022
1023 Double_t chidist = ChiSquareNDof((TH1D*)hDist_on,(TH1D*)hDist_off);
1024
1025 x_label_pos = distmax[numberSizeBins-1]*0.7;
1026 y_label_pos = hDist_on->GetBinContent(hDist_on->GetMaximumBin())/2.;
1027
1028 char text_dist[256];
1029 sprintf(text_size,"ChiSquare/NDof = %4.2f",chidist);
1030
1031 TLatex *tdist = new TLatex(x_label_pos, y_label_pos, text_dist);
1032 tdist->SetTextSize(textsize);
1033// tdist->Draw();
1034
1035 gPad->Modified();
1036 gPad->Update();
1037
1038 c2->Print(psname);
1039
1040 // ############################################################################
1041 // DrawCopy WIDTH
1042 // ############################################################################
1043 TCanvas *c3 = new TCanvas;
1044
1045 gPad->cd();
1046
1047 hWidth_off->Sumw2();
1048 hWidth_off->Scale(norm);
1049 hWidth_on->SetLineColor(kBlack);
1050 hWidth_on->SetMarkerStyle(21);
1051 hWidth_on->SetMarkerSize(0.7);
1052 hWidth_on->SetMarkerColor(kBlack);
1053 hWidth_off->SetFillColor(46);
1054 hWidth_off->SetLineColor(46);
1055 hWidth_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
1056 hWidth_off->SetMinimum(0.);
1057 hWidth_on->SetTitle("WIDTH distribution");
1058 hWidth_off->SetTitle("WIDTH distribution");
1059
1060 hWidth_on->DrawCopy("E1P");
1061
1062 // move stat box to make them all visible
1063 gPad->Update();
1064 TPaveStats* pavs_on_width = (TPaveStats*) hWidth_on->GetListOfFunctions()->FindObject("stats");
1065 if(pavs_on_width){
1066 shiftx = pavs_on_width->GetX2NDC() - pavs_on_width->GetX1NDC();
1067 pavs_on_width->SetX1NDC(pavs_on_width->GetX1NDC() - shiftx);
1068 pavs_on_width->SetX2NDC(pavs_on_width->GetX2NDC() - shiftx);
1069 }
1070 gPad->Modified();
1071 gPad->Update();
1072
1073 hWidth_off->DrawCopy("HISTSAME");
1074 hWidth_off->DrawCopy("ESAME");
1075 hWidth_on->DrawCopy("E1PSAME");
1076
1077 Double_t chiwidth = ChiSquareNDof((TH1D*)hWidth_on,(TH1D*)hWidth_off);
1078
1079 x_label_pos = widthmax[numberSizeBins-1]*0.7;
1080 y_label_pos = hWidth_on->GetBinContent(hWidth_on->GetMaximumBin())/2.;
1081
1082 char text_width[256];
1083 sprintf(text_size,"ChiSquare/NDof = %4.2f",chiwidth);
1084
1085 TLatex *twidth = new TLatex(x_label_pos, y_label_pos, text_width);
1086 twidth->SetTextSize(textsize);
1087// twidth->Draw();
1088
1089 gPad->Modified();
1090 gPad->Update();
1091
1092 c3->Print(psname);
1093
1094 // ############################################################################
1095 // DrawCopy LENGTH
1096 // ############################################################################
1097 TCanvas *c4 = new TCanvas;
1098
1099 gPad->cd();
1100
1101 hLength_on->Sumw2();
1102 hLength_off->Sumw2();
1103 hLength_off->Scale(norm);
1104 hLength_on->SetLineColor(kBlack);
1105 hLength_on->SetMarkerStyle(21);
1106 hLength_on->SetMarkerSize(0.7);
1107 hLength_on->SetMarkerColor(kBlack);
1108 hLength_off->SetFillColor(46);
1109 hLength_off->SetLineColor(46);
1110 hLength_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
1111 hLength_off->SetMinimum(0.);
1112 hLength_on->SetTitle("LENGTH distribution");
1113 hLength_off->SetTitle("LENGTH distribution");
1114
1115 hLength_on->DrawCopy("E1P");
1116
1117 // move stat box to make them all visible
1118 gPad->Update();
1119 TPaveStats* pavs_on_length = (TPaveStats*) hLength_on->GetListOfFunctions()->FindObject("stats");
1120 if(pavs_on_length){
1121 shiftx = pavs_on_length->GetX2NDC() - pavs_on_length->GetX1NDC();
1122 pavs_on_length->SetX1NDC(pavs_on_length->GetX1NDC() - shiftx);
1123 pavs_on_length->SetX2NDC(pavs_on_length->GetX2NDC() - shiftx);
1124 }
1125 gPad->Modified();
1126 gPad->Update();
1127
1128 hLength_off->DrawCopy("HISTSAME");
1129 hLength_off->DrawCopy("ESAME");
1130 hLength_on->DrawCopy("E1PSAME");
1131
1132 Double_t chilength = ChiSquareNDof((TH1D*)hLength_on,(TH1D*)hLength_off);
1133
1134 x_label_pos = lengthmax[numberSizeBins-1]*0.7;
1135 y_label_pos = hLength_on->GetBinContent(hLength_on->GetMaximumBin())/2.;
1136
1137 char text_length[256];
1138 sprintf(text_size,"ChiSquare/NDof = %4.2f",chilength);
1139
1140 TLatex *tlength = new TLatex(x_label_pos, y_label_pos, text_length);
1141 tlength->SetTextSize(textsize);
1142// tlength->Draw();
1143
1144 gPad->Modified();
1145 gPad->Update();
1146
1147 c4->Print(psname);
1148
1149 // ############################################################################
1150 // DrawCopy normalized ALPHA plot
1151 // ############################################################################
1152 TCanvas *c5 = new TCanvas;
1153
1154 gPad->cd();
1155
1156 hAlpha_on_abs->Sumw2();
1157 hAlpha_off_abs->SetStats(0);
1158 hAlpha_off_abs->Sumw2();
1159 hAlpha_off_abs->Scale(norm);
1160 hAlpha_on_abs->SetStats(0); //-> Do NOT show the legend with statistics
1161 hAlpha_on_abs->SetLineColor(kBlack);
1162 hAlpha_on_abs->SetMarkerStyle(21);
1163 //hAlpha_on_abs->SetMarkerSize();
1164 hAlpha_on_abs->SetMarkerColor(kBlack);
1165 hAlpha_on_abs->SetMarkerSize(0.7);
1166 hAlpha_off_abs->SetFillColor(46);
1167 hAlpha_off_abs->SetLineColor(46);
1168 hAlpha_off_abs->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
1169 hAlpha_off_abs->SetMinimum(0.);
1170 hAlpha_on_abs->SetTitle("Alpha plot");
1171 hAlpha_off_abs->SetTitle("Alpha plot");
1172
1173
1174 hAlpha_on_abs->DrawCopy("E1P");
1175 hAlpha_off_abs->DrawCopy("HISTSAME");
1176 hAlpha_off_abs->DrawCopy("ESAME");
1177 hAlpha_on_abs->DrawCopy("E1PSAME");
1178
1179
1180 //draw the LEGEND with excess and significance values in the alpha plot:
1181 char text_Fnorm[256], text_excess[256], text_sign[256];
1182 char text_entries_on[256], text_entries_off[256], text_integral_off[256];
1183 int hAlpha_on_entries = (int) hAlpha_on_abs->GetEntries();
1184 int hAlpha_off_entries = (int) hAlpha_off_abs->GetEntries();
1185 sprintf(text_Fnorm, " F norm = %.3f", norm);
1186 sprintf(text_excess, " Excess = %.3f", excess);
1187 sprintf(text_sign, " Significance = %.3f", sign);
1188 sprintf(text_entries_on, " Entries ON = %d", hAlpha_on_entries);
1189 sprintf(text_entries_off, " Entries OFF = %d", hAlpha_off_entries);
1190 sprintf(text_integral_off," Integral OFF = %d", int_off);
1191
1192 x_label_pos = 90.*0.7;
1193 y_label_pos = (hAlpha_on_abs->GetBinContent(hAlpha_on_abs->GetMaximumBin())); //2.;
1194 Double_t y_label_step = y_label_pos / 8.;
1195
1196 TLatex *t0 = new TLatex(x_label_pos, y_label_pos - y_label_step*0, text_Fnorm);
1197 t0->SetTextSize(textsize);
1198 t0->Draw();
1199 TLatex *t1 = new TLatex(x_label_pos, y_label_pos - y_label_step*1, text_excess);
1200 t1->SetTextSize(textsize);
1201 t1->Draw();
1202 TLatex *t2 = new TLatex(x_label_pos, y_label_pos - y_label_step*2, text_sign);
1203 t2->SetTextSize(textsize);
1204 t2->Draw();
1205 TLatex *t3 = new TLatex(x_label_pos, y_label_pos - y_label_step*3, text_entries_on);
1206 t3->SetTextSize(textsize);
1207 t3->Draw();
1208 TLatex *t4 = new TLatex(x_label_pos, y_label_pos - y_label_step*4, text_entries_off);
1209 t4->SetTextSize(textsize);
1210 t4->Draw();
1211 TLatex *t5 = new TLatex(x_label_pos, y_label_pos - y_label_step*5, text_integral_off);
1212 t5->SetTextSize(textsize);
1213 t5->Draw();
1214
1215
1216 Double_t chialpha = ChiSquareNDof((TH1D*)hAlpha_on_abs,(TH1D*)hAlpha_off_abs);
1217
1218 y_label_pos = (hAlpha_on_abs->GetBinContent(hAlpha_on_abs->GetMaximumBin()))/2.;
1219
1220 char text_alpha[256];
1221 sprintf(text_size,"ChiSquare/NDof = %4.2f",chialpha);
1222
1223 TLatex *talpha = new TLatex(x_label_pos, y_label_pos, text_alpha);
1224 talpha->SetTextSize(textsize);
1225// talpha->Draw();
1226
1227 gPad->Modified();
1228 gPad->Update();
1229
1230 c5->Print(psname);
1231
1232 // ############################################################################
1233 // DrawCopy normalized alpha histos for alpha form -90 to 90 deg.
1234 // ############################################################################
1235 TCanvas *c6 = new TCanvas;
1236
1237 gPad->cd();
1238
1239 hAlpha_on->Sumw2();
1240 hAlpha_off->SetStats(0);
1241 hAlpha_off->Sumw2();
1242 hAlpha_off->Scale(norm);
1243 hAlpha_off->SetFillColor(46);
1244 hAlpha_off->SetLineColor(46);
1245 hAlpha_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
1246 hAlpha_off->SetMinimum(0.);
1247 hAlpha_on->SetStats(0); //-> Do NOT show the legend with statistics
1248 hAlpha_on->SetLineColor(kBlack);
1249 hAlpha_on->SetMarkerStyle(21);
1250 hAlpha_on->SetMarkerSize(0.7);
1251 hAlpha_on->SetMarkerColor(kBlack);
1252 hAlpha_on->SetTitle("Alpha plot form -90 to 90 deg");
1253 hAlpha_off->SetTitle("Alpha plot form -90 to 90 deg");
1254
1255 hAlpha_on->DrawCopy("E1P");
1256 hAlpha_off->DrawCopy("HISTSAME");
1257 hAlpha_off->DrawCopy("ESAME");
1258 hAlpha_on->DrawCopy("E1PSAME");
1259
1260 Double_t chialpha90 = ChiSquareNDof((TH1D*)hAlpha_on,(TH1D*)hAlpha_off);
1261
1262 x_label_pos = 90.*0.5;
1263 y_label_pos = hAlpha_on->GetBinContent(hAlpha_on->GetMaximumBin())/2.;
1264
1265 char text_alpha90[256];
1266 sprintf(text_alpha90,"ChiSquare/NDof = %4.2f",chialpha90);
1267
1268 TLatex *talpha90 = new TLatex(x_label_pos, y_label_pos, text_alpha90);
1269 talpha90->SetTextSize(textsize);
1270// talpha90->Draw();
1271
1272 gPad->Update();
1273 gPad->Modified();
1274
1275 c6->Print(psname);
1276
1277 cout << "---> ChiSquare/NDof [Size] =\t\t" << chisize << endl;
1278 cout << "---> ChiSquare/NDof [Dist] =\t\t" << chidist << endl;
1279 cout << "---> ChiSquare/NDof [Width] =\t\t" << chiwidth << endl;
1280 cout << "---> ChiSquare/NDof [Length] =\t\t" << chilength << endl;
1281 cout << "---> ChiSquare/NDof [Abs(Alpha)] =\t" << chialpha << endl;
1282 cout << "---> ChiSquare/NDof [Alpha] =\t\t" << chialpha90 << endl;
1283
1284
1285 TCanvas *c7 = new TCanvas;
1286 hSrcPos_on->DrawCopy("BOX");
1287 c7->Print(psname);
1288
1289 TCanvas *c8 = new TCanvas;
1290 hSrcPos_off->DrawCopy("BOX");
1291 c8->Print(closepsname);
1292
1293 cout << "Done!!" <<endl;
1294
1295}
1296
1297
1298
Note: See TracBrowser for help on using the repository browser.