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

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