1 | /* ======================================================================== *\
|
---|
2 | !
|
---|
3 | ! *
|
---|
4 | ! * This file is part of MARS, the MAGIC Analysis and Reconstruction
|
---|
5 | ! * Software. It is distributed to you in the hope that it can be a useful
|
---|
6 | ! * and timesaving tool in analysing Data of imaging Cherenkov telescopes.
|
---|
7 | ! * It is distributed WITHOUT ANY WARRANTY.
|
---|
8 | ! *
|
---|
9 | ! * Permission to use, copy, modify and distribute this software and its
|
---|
10 | ! * documentation for any purpose is hereby granted without fee,
|
---|
11 | ! * provided that the above copyright notice appear in all copies and
|
---|
12 | ! * that both that copyright notice and this permission notice appear
|
---|
13 | ! * in supporting documentation. It is provided "as is" without express
|
---|
14 | ! * or implied warranty.
|
---|
15 | ! *
|
---|
16 | !
|
---|
17 | ! Author(s): Daniel Mazin, 05/2004 <mailto:mazin@imppmu.mpg.de>
|
---|
18 | !
|
---|
19 | ! Copyright: MAGIC Software Development, 2000-2004
|
---|
20 | !
|
---|
21 | !
|
---|
22 | \* ======================================================================== */
|
---|
23 |
|
---|
24 |
|
---|
25 | // Daniel Mazin 18.05.2004 mazin@mppmu.mpg.de
|
---|
26 | // **********************************************************************************
|
---|
27 | // this macro is used to produce an alpha plot/plots for a chosen position in the sky map.
|
---|
28 | // In case of several subsamples, the overall alpha plot is also produced.
|
---|
29 | // The derotation is possible if Zd and Az are available for each event.
|
---|
30 | // Using of OFF data to estimate Nex and Significance is optional (if such data available)
|
---|
31 | //
|
---|
32 | // input: hillas parameter file/files
|
---|
33 | // optional: a) ascii file with the estimated position of the source
|
---|
34 | // b) ascii file with fit parameters from the OFF sample
|
---|
35 | // (needed for estimation of the background from OFF dada)
|
---|
36 | // output: alpha plot for each subsample, overall alpha plot
|
---|
37 | // **********************************************************************************
|
---|
38 |
|
---|
39 |
|
---|
40 | #define XSOUR 0.05 // [deg]
|
---|
41 | #define YSOUR 0.6 // [deg]
|
---|
42 | #define NUM 0
|
---|
43 |
|
---|
44 |
|
---|
45 | TString sourcename = "M87On";
|
---|
46 |
|
---|
47 | const Bool_t ROTOPTION = kFALSE; // kFALSE: do not derotate, use camera coordinates
|
---|
48 | // kTRUE: derotate into "quasi" sky coordinates
|
---|
49 |
|
---|
50 | const Bool_t USEOFFDATA = kFALSE; // kFALSE: do not use OFF data to estimate the significance
|
---|
51 | // kTRUE: use in addition OFF data to estimate the significance
|
---|
52 |
|
---|
53 | const char * offfile = "paramOffCrab2701.dat"; // file with parameters of the OFF sample
|
---|
54 | // needed only if USEOFFDATA = kTRUE
|
---|
55 |
|
---|
56 | const Bool_t USEFILE = kFALSE; // kFALSE : use XSOUR, YSOUR
|
---|
57 | // kTRUE: use "posfile"
|
---|
58 | // to specify the position for which alpha plot has to be produced
|
---|
59 |
|
---|
60 | const Bool_t SAMPLES = kFALSE; // kFALSE: one sample only
|
---|
61 | // kTRUE: several samples, usage of posfile is preferred!
|
---|
62 |
|
---|
63 | const char * posfile = "data/trackCrab1502_1305_Berlin_5deg.dat"; // file with best position
|
---|
64 | // if USEFILE == kFALSE : not needed
|
---|
65 |
|
---|
66 | TString dirname = "~/data/M87/2004_02_17/";
|
---|
67 | TString filename = "M87HillasON.root"; // NOTE in case of several subsamples (SAMPLES=kTRUE)
|
---|
68 | // folllowing name construction is assumed:
|
---|
69 | // file = dirname + sample + <number> + / + filename
|
---|
70 | // e.g: ~/data/Crab/2004_01_27/sample2/CrabHillasON.root
|
---|
71 | #define tolerance 1e-3
|
---|
72 | /* ****************************************************** */
|
---|
73 | /* dynamical cuts Crab 27th Jan 2004 */
|
---|
74 | /*
|
---|
75 | #define LENGTHMINParA 0.136 // deg
|
---|
76 | #define LENGTHMINParB 0.036 //
|
---|
77 | #define LENGTHMINParC -0.0038 //
|
---|
78 | #define LENGTHMAXParA 0.332 // deg
|
---|
79 | #define LENGTHMAXParB 0.037 //
|
---|
80 | #define LENGTHMAXParC 0.0261 //
|
---|
81 | #define WIDTHMINParA 0.063 // deg
|
---|
82 | #define WIDTHMINParB 0.013 //
|
---|
83 | #define WIDTHMINParC 0.0003 //
|
---|
84 | #define WIDTHMAXParA 0.123 // deg
|
---|
85 | #define WIDTHMAXParB 0.019 //
|
---|
86 | #define WIDTHMAXParC 0.0005 //
|
---|
87 | #define DISTMINParA 0.6 // deg
|
---|
88 | #define DISTMINParB 0.059 //
|
---|
89 | #define DISTMINParC 0. //
|
---|
90 | #define DISTMAXParA 1.25 // deg
|
---|
91 | #define DISTMAXParB 0.059 //
|
---|
92 | #define DISTMAXParC 0. //
|
---|
93 | */
|
---|
94 |
|
---|
95 | /* dynamical cuts Mrk 421 and Crab 15th Feb 2004 */
|
---|
96 |
|
---|
97 | #define LENGTHMINParA 0.12 // deg
|
---|
98 | #define LENGTHMINParB 0.034 //
|
---|
99 | #define LENGTHMINParC 0. //
|
---|
100 | #define LENGTHMAXParA 0.32 // deg
|
---|
101 | #define LENGTHMAXParB 0.034 //
|
---|
102 | #define LENGTHMAXParC 0. //
|
---|
103 | #define WIDTHMINParA 0.055 // deg
|
---|
104 | #define WIDTHMINParB 0.013 //
|
---|
105 | #define WIDTHMINParC 0.0 //
|
---|
106 | #define WIDTHMAXParA 0.12 // deg
|
---|
107 | #define WIDTHMAXParB 0.013 //
|
---|
108 | #define WIDTHMAXParC 0.0 //
|
---|
109 | #define DISTMINParA 0.6 // deg
|
---|
110 | #define DISTMINParB 0.059 //
|
---|
111 | #define DISTMINParC 0. //
|
---|
112 | #define DISTMAXParA 1.25 // deg
|
---|
113 | #define DISTMAXParB 0.059 //
|
---|
114 | #define DISTMAXParC 0. //
|
---|
115 |
|
---|
116 |
|
---|
117 | #define SIZEMIN 2000
|
---|
118 | #define LEAKMAX 0.25
|
---|
119 | #define ASYMMIN -0.1
|
---|
120 |
|
---|
121 | #include "mtools.C"
|
---|
122 |
|
---|
123 | #define ALPHAMAX 15.
|
---|
124 |
|
---|
125 | #define ALOFFMAX 90.
|
---|
126 | #define ALOFFMIN 30.
|
---|
127 |
|
---|
128 | // histogram to store the sum of alpha plots
|
---|
129 | TH1F histalphaAll("alpha plot", "alpha plot, size cut 2000, 15 deg", 36, 0., 90.);
|
---|
130 | histalphaAll.SetXTitle("alpha (deg)");
|
---|
131 | histalphaAll.SetYTitle("Counts");
|
---|
132 | histalphaAll.SetDirectory(NULL);
|
---|
133 | histalphaAll.SetFillStyle(4000);
|
---|
134 | histalphaAll.UseCurrentStyle();
|
---|
135 |
|
---|
136 |
|
---|
137 | void alphaplots(Int_t num, Double_t XSOURCE, Double_t YSOURCE)
|
---|
138 | {
|
---|
139 |
|
---|
140 | // for USEOFFDATA = kTRUE the ascii file with fit parameters of the OFF sample is read.
|
---|
141 |
|
---|
142 | if (USEOFFDATA == kTRUE)
|
---|
143 | {
|
---|
144 | Int_t event, dummy, numbinsoff;
|
---|
145 | Float_t binwidthOff, xpos, ypos, aparOff, bparOff, chi2par;
|
---|
146 |
|
---|
147 | FILE *fp;
|
---|
148 | fp = fopen(offfile,"r");
|
---|
149 | while(fscanf(fp,"%d %f %f %f %f %f %f %f %f %f",
|
---|
150 | &event, &xpos, &ypos, &aparOff, &bparOff, &chi2par,
|
---|
151 | &binwidthOff, &dummy, &dummy, &dummy) != EOF)
|
---|
152 | {
|
---|
153 | //cout << event << " " << xpos << " " << ypos << endl;
|
---|
154 | if(TMath::Abs(xpos-XSOURCE) < tolerance && TMath::Abs(ypos-YSOURCE) < tolerance) break;
|
---|
155 | }
|
---|
156 | fclose(fp);
|
---|
157 | }
|
---|
158 |
|
---|
159 |
|
---|
160 | gStyle->SetCanvasBorderMode(0);
|
---|
161 | gStyle->SetCanvasBorderSize(0);
|
---|
162 | gStyle->SetCanvasColor(10);
|
---|
163 | gStyle->SetPadBorderMode(0);
|
---|
164 | gStyle->SetPadBorderSize(0);
|
---|
165 | gStyle->SetPadColor(10);
|
---|
166 | gStyle->SetOptFit(1);
|
---|
167 | gStyle->SetStatColor(10);
|
---|
168 | // gStyle->SetOptStat(0);
|
---|
169 | gStyle->SetPalette(1,0);
|
---|
170 |
|
---|
171 |
|
---|
172 |
|
---|
173 | // the name of the Hillas parameter file, which has to be read in
|
---|
174 | TString ddumy;
|
---|
175 |
|
---|
176 | if (SAMPLES == kFALSE)
|
---|
177 | filename = dirname + filename;
|
---|
178 | else
|
---|
179 | {
|
---|
180 | ddumy = dirname;
|
---|
181 | ddumy += "sample";
|
---|
182 | ddumy += num;
|
---|
183 | filename = ddumy + filename;
|
---|
184 | }
|
---|
185 |
|
---|
186 | cout << "file to read :" << filename << endl;
|
---|
187 |
|
---|
188 |
|
---|
189 | // create histograms
|
---|
190 | MBinning bins;
|
---|
191 |
|
---|
192 | TH1F histlength;
|
---|
193 | histlength.SetName("Length");
|
---|
194 | histlength.SetTitle("Length");
|
---|
195 | histlength.SetXTitle("Length [deg]");
|
---|
196 | histlength.SetYTitle("Counts");
|
---|
197 | histlength.SetDirectory(NULL);
|
---|
198 | histlength.SetFillStyle(4000);
|
---|
199 | histlength.UseCurrentStyle();
|
---|
200 |
|
---|
201 |
|
---|
202 | bins.SetEdges(100, 0., 1.);
|
---|
203 | bins.Apply(histlength);
|
---|
204 |
|
---|
205 |
|
---|
206 | TH1F histwidth;
|
---|
207 | histwidth.SetName("Width");
|
---|
208 | histwidth.SetTitle("Width");
|
---|
209 | histwidth.SetXTitle("Width [deg]");
|
---|
210 | histwidth.SetYTitle("Counts");
|
---|
211 | histwidth.SetDirectory(NULL);
|
---|
212 | histwidth.SetFillStyle(4000);
|
---|
213 | histwidth.UseCurrentStyle();
|
---|
214 |
|
---|
215 | bins.SetEdges(100, 0., 0.5);
|
---|
216 | bins.Apply(histwidth);
|
---|
217 |
|
---|
218 |
|
---|
219 | TH1F histsize;
|
---|
220 | histsize.SetName("Size");
|
---|
221 | histsize.SetTitle("Size");
|
---|
222 | histsize.SetXTitle("Size");
|
---|
223 | histsize.SetYTitle("Counts");
|
---|
224 | histsize.SetDirectory(NULL);
|
---|
225 | histsize.SetFillStyle(4000);
|
---|
226 | histsize.UseCurrentStyle();
|
---|
227 |
|
---|
228 | // bins.SetEdges(100, 100., 2e5);
|
---|
229 | // bins.Apply(histsize);
|
---|
230 |
|
---|
231 | bins.SetEdgesLog(100, 100., 10e5);
|
---|
232 | bins.Apply(histsize);
|
---|
233 |
|
---|
234 | TH1F histalpha;
|
---|
235 | histalpha.SetName("Alpha");
|
---|
236 | histalpha.SetTitle("Alpha");
|
---|
237 | histalpha.SetXTitle("Alpha [deg]");
|
---|
238 | histalpha.SetYTitle("Counts");
|
---|
239 | histalpha.SetDirectory(NULL);
|
---|
240 | histalpha.SetFillStyle(4000);
|
---|
241 | histalpha.UseCurrentStyle();
|
---|
242 |
|
---|
243 | bins.SetEdges(100, -100., 100.);
|
---|
244 | bins.Apply(histalpha);
|
---|
245 |
|
---|
246 | TH1F histdist;
|
---|
247 | histdist.SetName("Dist");
|
---|
248 | histdist.SetTitle("Dist");
|
---|
249 | histdist.SetXTitle("Dist [deg]");
|
---|
250 | histdist.SetYTitle("Counts");
|
---|
251 | histdist.SetDirectory(NULL);
|
---|
252 | histdist.SetFillStyle(4000);
|
---|
253 | histdist.UseCurrentStyle();
|
---|
254 |
|
---|
255 | bins.SetEdges(100, 0., 2.);
|
---|
256 | bins.Apply(histdist);
|
---|
257 |
|
---|
258 | TH1F histmeanx;
|
---|
259 | histmeanx.SetName("MeanX");
|
---|
260 | histmeanx.SetTitle("MeanX");
|
---|
261 | histmeanx.SetXTitle("MeanX [deg]");
|
---|
262 | histmeanx.SetYTitle("Counts");
|
---|
263 | histmeanx.SetDirectory(NULL);
|
---|
264 | histmeanx.SetFillStyle(4000);
|
---|
265 | histmeanx.UseCurrentStyle();
|
---|
266 |
|
---|
267 | bins.SetEdges(100, -1.8, 1.8);
|
---|
268 | bins.Apply(histmeanx);
|
---|
269 |
|
---|
270 | TH1F histmeany;
|
---|
271 | histmeany.SetName("MeanY");
|
---|
272 | histmeany.SetTitle("MeanY");
|
---|
273 | histmeany.SetXTitle("MeanY [deg]");
|
---|
274 | histmeany.SetYTitle("Counts");
|
---|
275 | histmeany.SetDirectory(NULL);
|
---|
276 | histmeany.SetFillStyle(4000);
|
---|
277 | histmeany.UseCurrentStyle();
|
---|
278 |
|
---|
279 |
|
---|
280 | bins.SetEdges(100, -1.8, 1.8);
|
---|
281 | bins.Apply(histmeany);
|
---|
282 |
|
---|
283 | TH1F histalphafinal;
|
---|
284 | histalphafinal.SetName("ALPHA");
|
---|
285 | histalphafinal.SetTitle("ALPHA");
|
---|
286 | histalphafinal.SetXTitle("alpha [deg]");
|
---|
287 | histalphafinal.SetYTitle("Counts");
|
---|
288 | histalphafinal.SetDirectory(NULL);
|
---|
289 | histalphafinal.SetFillStyle(4000);
|
---|
290 | histalphafinal.UseCurrentStyle();
|
---|
291 |
|
---|
292 | bins.SetEdges(36, 0.0, 90.);
|
---|
293 | bins.Apply(histalphafinal);
|
---|
294 |
|
---|
295 | TH1F histAssym;
|
---|
296 | histAssym.SetName("Assymetry");
|
---|
297 | histAssym.SetTitle("Assymetry");
|
---|
298 | histAssym.SetXTitle("Assymetry");
|
---|
299 | histAssym.SetYTitle("Counts");
|
---|
300 | histAssym.SetDirectory(NULL);
|
---|
301 | histAssym.SetFillStyle(4000);
|
---|
302 | histAssym.UseCurrentStyle();
|
---|
303 |
|
---|
304 | bins.SetEdges(100, -1, 1.);
|
---|
305 | bins.Apply(histAssym);
|
---|
306 |
|
---|
307 | TH1F histAssymM3;
|
---|
308 | histAssymM3.SetName("Assymetry 3M");
|
---|
309 | histAssymM3.SetTitle("Assymetry 3rd moment");
|
---|
310 | histAssymM3.SetXTitle("Assymetry 3rd moment");
|
---|
311 | histAssymM3.SetYTitle("Counts");
|
---|
312 | histAssymM3.SetDirectory(NULL);
|
---|
313 | histAssymM3.SetFillStyle(4000);
|
---|
314 | histAssymM3.UseCurrentStyle();
|
---|
315 |
|
---|
316 | bins.SetEdges(100, -1., 1.);
|
---|
317 | bins.Apply(histAssymM3);
|
---|
318 |
|
---|
319 | TH2F hist2xy("CoG","Center of Gravity", 100, -1.8, 1.8, 100, -1.8, 1.8);
|
---|
320 | hist2xy.SetXTitle("MeanX [deg]");
|
---|
321 | hist2xy.SetYTitle("MeanY [deg]");
|
---|
322 | hist2xy.SetDirectory(NULL);
|
---|
323 | hist2xy.SetFillStyle(4000);
|
---|
324 | hist2xy.UseCurrentStyle();
|
---|
325 |
|
---|
326 | TH1F histLoverS;
|
---|
327 | histLoverS.SetName("LoverS");
|
---|
328 | histLoverS.SetTitle("LoverS");
|
---|
329 | histLoverS.SetXTitle("LoverS");
|
---|
330 | histLoverS.SetYTitle("Counts");
|
---|
331 | histLoverS.SetDirectory(NULL);
|
---|
332 | histLoverS.SetFillStyle(4000);
|
---|
333 | histLoverS.UseCurrentStyle();
|
---|
334 |
|
---|
335 | bins.SetEdges(100, -0., 0.0006);
|
---|
336 | bins.Apply(histLoverS);
|
---|
337 |
|
---|
338 | TH2F hist2lw("Length-Width", "correlation Length-Width", 100, 0.0, 1.0, 100, 0.0, 0.5);
|
---|
339 | hist2lw.SetXTitle("Length [deg]");
|
---|
340 | hist2lw.SetYTitle("Width [deg]");
|
---|
341 | hist2lw.SetDirectory(NULL);
|
---|
342 | hist2lw.SetFillStyle(4000);
|
---|
343 | hist2lw.UseCurrentStyle();
|
---|
344 |
|
---|
345 | TH2F hist2lalpha("Length-Alpha", "correlation Length-Alpha", 100, 0.0, 1.0, 100, -100., 100.);
|
---|
346 | hist2lalpha.SetXTitle("Length [deg]");
|
---|
347 | hist2lalpha.SetYTitle("Alpha [deg]");
|
---|
348 | hist2lalpha.SetDirectory(NULL);
|
---|
349 | hist2lalpha.SetFillStyle(4000);
|
---|
350 | hist2lalpha.UseCurrentStyle();
|
---|
351 |
|
---|
352 | TH2F hist2ldist("Length-Dist","correlation Length-Dist", 100, 0.0, 1.0, 100, 0.0, 1.7);
|
---|
353 | hist2ldist.SetXTitle("Length [deg]");
|
---|
354 | hist2ldist.SetYTitle("Dist [deg]");
|
---|
355 | hist2ldist.SetDirectory(NULL);
|
---|
356 | hist2ldist.SetFillStyle(4000);
|
---|
357 | hist2ldist.UseCurrentStyle();
|
---|
358 |
|
---|
359 | TH2F hist2walpha("Width-Alpha","correlation Width-Alpha", 100, 0.0, 0.5, 100, -100., 100.);
|
---|
360 | hist2walpha.SetXTitle("Width [deg]");
|
---|
361 | hist2walpha.SetYTitle("Alpha [deg]");
|
---|
362 | hist2walpha.SetDirectory(NULL);
|
---|
363 | hist2walpha.SetFillStyle(4000);
|
---|
364 | hist2walpha.UseCurrentStyle();
|
---|
365 |
|
---|
366 | TH2F hist2wdist("Width-Dist","correlation Width-Dist", 100, 0.0, 0.5, 100, 0.0, 1.7);
|
---|
367 | hist2wdist.SetXTitle("Width [deg]");
|
---|
368 | hist2wdist.SetYTitle("Dist [deg]");
|
---|
369 | hist2wdist.SetDirectory(NULL);
|
---|
370 | hist2wdist.SetFillStyle(4000);
|
---|
371 | hist2wdist.UseCurrentStyle();
|
---|
372 |
|
---|
373 | TH2F hist2alphadist("Alpha-Dist","correlation Alpha-Dist", 100, -100., 100, 100, 0.0, 1.7);
|
---|
374 | hist2alphadist.SetXTitle("Alpha [deg]");
|
---|
375 | hist2alphadist.SetYTitle("Dist [deg]");
|
---|
376 | hist2alphadist.SetDirectory(NULL);
|
---|
377 | hist2alphadist.SetFillStyle(4000);
|
---|
378 | hist2alphadist.UseCurrentStyle();
|
---|
379 |
|
---|
380 | TH1F histphi;
|
---|
381 | histphi.SetName("TelPhia");
|
---|
382 | histphi.SetTitle("Telescope Phi");
|
---|
383 | histphi.SetXTitle("Phi [rad]");
|
---|
384 | histphi.SetYTitle("Counts");
|
---|
385 | histphi.SetDirectory(NULL);
|
---|
386 | histphi.SetFillStyle(4000);
|
---|
387 | histphi.UseCurrentStyle();
|
---|
388 |
|
---|
389 | bins.SetEdges(100, -10, 10);
|
---|
390 | bins.Apply(histphi);
|
---|
391 |
|
---|
392 | TH1F histtheta;
|
---|
393 | histtheta.SetName("TelTheta");
|
---|
394 | histtheta.SetTitle("Telescope Theta");
|
---|
395 | histtheta.SetXTitle("Theta [rad]");
|
---|
396 | histtheta.SetYTitle("Counts");
|
---|
397 | histtheta.SetDirectory(NULL);
|
---|
398 | histtheta.SetFillStyle(4000);
|
---|
399 | histtheta.UseCurrentStyle();
|
---|
400 |
|
---|
401 | bins.SetEdges(100, -2, 2.);
|
---|
402 | bins.Apply(histtheta);
|
---|
403 |
|
---|
404 | TH1F aftercuthistlength;
|
---|
405 | aftercuthistlength.SetName("Length");
|
---|
406 | aftercuthistlength.SetTitle("Length");
|
---|
407 | aftercuthistlength.SetXTitle("Length [deg]");
|
---|
408 | aftercuthistlength.SetYTitle("Counts");
|
---|
409 | aftercuthistlength.SetDirectory(NULL);
|
---|
410 | aftercuthistlength.SetFillStyle(4000);
|
---|
411 | aftercuthistlength.UseCurrentStyle();
|
---|
412 |
|
---|
413 |
|
---|
414 | bins.SetEdges(100, 0., 1.);
|
---|
415 | bins.Apply(aftercuthistlength);
|
---|
416 |
|
---|
417 | TH1F aftercuthistwidth;
|
---|
418 | aftercuthistwidth.SetName("Width");
|
---|
419 | aftercuthistwidth.SetTitle("Width");
|
---|
420 | aftercuthistwidth.SetXTitle("Width [deg]");
|
---|
421 | aftercuthistwidth.SetYTitle("Counts");
|
---|
422 | aftercuthistwidth.SetDirectory(NULL);
|
---|
423 | aftercuthistwidth.SetFillStyle(4000);
|
---|
424 | aftercuthistwidth.UseCurrentStyle();
|
---|
425 |
|
---|
426 | bins.SetEdges(100, 0., 0.5);
|
---|
427 | bins.Apply(aftercuthistwidth);
|
---|
428 |
|
---|
429 | TH1F aftercuthistsize;
|
---|
430 | aftercuthistsize.SetName("Size");
|
---|
431 | aftercuthistsize.SetTitle("Size");
|
---|
432 | aftercuthistsize.SetXTitle("Size [photons]");
|
---|
433 | aftercuthistsize.SetYTitle("Counts");
|
---|
434 | aftercuthistsize.SetDirectory(NULL);
|
---|
435 | aftercuthistsize.SetFillStyle(4000);
|
---|
436 | aftercuthistsize.UseCurrentStyle();
|
---|
437 |
|
---|
438 | bins.SetEdgesLog(100, 100., 10e5);
|
---|
439 | bins.Apply(aftercuthistsize);
|
---|
440 |
|
---|
441 | TH1F aftercuthistalpha;
|
---|
442 | aftercuthistalpha.SetName("Alpha");
|
---|
443 | aftercuthistalpha.SetTitle("Alpha");
|
---|
444 | aftercuthistalpha.SetXTitle("Alpha [deg]");
|
---|
445 | aftercuthistalpha.SetYTitle("Counts");
|
---|
446 | aftercuthistalpha.SetDirectory(NULL);
|
---|
447 | aftercuthistalpha.SetFillStyle(4000);
|
---|
448 | aftercuthistalpha.UseCurrentStyle();
|
---|
449 |
|
---|
450 | bins.SetEdges(20, 0., 100.);
|
---|
451 | bins.Apply(aftercuthistalpha);
|
---|
452 |
|
---|
453 | TH1F aftercuthistdist;
|
---|
454 | aftercuthistdist.SetName("Dist");
|
---|
455 | aftercuthistdist.SetTitle("Dist");
|
---|
456 | aftercuthistdist.SetXTitle("Dist [deg]");
|
---|
457 | aftercuthistdist.SetYTitle("Counts");
|
---|
458 | aftercuthistdist.SetDirectory(NULL);
|
---|
459 | aftercuthistdist.SetFillStyle(4000);
|
---|
460 | aftercuthistdist.UseCurrentStyle();
|
---|
461 |
|
---|
462 | bins.SetEdges(100, 0., 2.);
|
---|
463 | bins.Apply(aftercuthistdist);
|
---|
464 |
|
---|
465 | TH1F aftercuthistmeanx;
|
---|
466 | aftercuthistmeanx.SetName("MeanX");
|
---|
467 | aftercuthistmeanx.SetTitle("MeanX");
|
---|
468 | aftercuthistmeanx.SetXTitle("MeanX [deg]");
|
---|
469 | aftercuthistmeanx.SetYTitle("Counts");
|
---|
470 | aftercuthistmeanx.SetDirectory(NULL);
|
---|
471 | aftercuthistmeanx.SetFillStyle(4000);
|
---|
472 | aftercuthistmeanx.UseCurrentStyle();
|
---|
473 |
|
---|
474 | bins.SetEdges(100, -1.8, 1.8);
|
---|
475 | bins.Apply(aftercuthistmeanx);
|
---|
476 |
|
---|
477 | TH1F aftercuthistmeany;
|
---|
478 | aftercuthistmeany.SetName("MeanY");
|
---|
479 | aftercuthistmeany.SetTitle("MeanY");
|
---|
480 | aftercuthistmeany.SetXTitle("MeanY [deg]");
|
---|
481 | aftercuthistmeany.SetYTitle("Counts");
|
---|
482 | aftercuthistmeany.SetDirectory(NULL);
|
---|
483 | aftercuthistmeany.SetFillStyle(4000);
|
---|
484 | aftercuthistmeany.UseCurrentStyle();
|
---|
485 |
|
---|
486 |
|
---|
487 | bins.SetEdges(100, -1.8, 1.8);
|
---|
488 | bins.Apply(aftercuthistmeany);
|
---|
489 |
|
---|
490 | TH2F aftercuthist2xy("CoG","Center of Gravity", 100, -1.8, 1.8, 100, -1.8, 1.8);
|
---|
491 | aftercuthist2xy.SetXTitle("MeanX [deg]");
|
---|
492 | aftercuthist2xy.SetYTitle("MeanY [deg]");
|
---|
493 | aftercuthist2xy.SetDirectory(NULL);
|
---|
494 | aftercuthist2xy.SetFillStyle(4000);
|
---|
495 | aftercuthist2xy.UseCurrentStyle();
|
---|
496 |
|
---|
497 |
|
---|
498 |
|
---|
499 | const Int_t n = 100;
|
---|
500 | Double_t binsize[n];
|
---|
501 |
|
---|
502 | Float_t nmin = 100.;
|
---|
503 | Float_t nmax = 1e7;
|
---|
504 |
|
---|
505 | for(Int_t i=0; i<n; i++)
|
---|
506 | {
|
---|
507 | binsize[i] = pow(10., log10(nmin) + i * (log10(nmax) - log10(nmin)) / (n-1.));
|
---|
508 | }
|
---|
509 |
|
---|
510 | TH2F hist2wsize("Width-Size", "correlation Width-Size", 100, 0.0, 0.5, n-1, binsize);
|
---|
511 | hist2wsize.SetXTitle("Width [deg]");
|
---|
512 | hist2wsize.SetYTitle("Size");
|
---|
513 | hist2wsize.SetDirectory(NULL);
|
---|
514 | hist2wsize.SetFillStyle(4000);
|
---|
515 | hist2wsize.UseCurrentStyle();
|
---|
516 |
|
---|
517 |
|
---|
518 |
|
---|
519 | TH2F hist2alphasize("Alpha-Size","correlation Alpha-Size", 100, -100., 100., n-1, binsize);
|
---|
520 | hist2alphasize.SetXTitle("Alpha [deg]");
|
---|
521 | hist2alphasize.SetYTitle("Size");
|
---|
522 | hist2alphasize.SetDirectory(NULL);
|
---|
523 | hist2alphasize.SetFillStyle(4000);
|
---|
524 | hist2alphasize.UseCurrentStyle();
|
---|
525 |
|
---|
526 | TH2F hist2distsize("Dist-Size","correlation Dist-Size", 100, 0.0, 1.7, n-1, binsize);
|
---|
527 | hist2distsize.SetXTitle("Dist [deg]");
|
---|
528 | hist2distsize.SetYTitle("Size");
|
---|
529 | hist2distsize.SetDirectory(NULL);
|
---|
530 | hist2distsize.SetFillStyle(4000);
|
---|
531 | hist2distsize.UseCurrentStyle();
|
---|
532 | // end create histograms
|
---|
533 |
|
---|
534 |
|
---|
535 |
|
---|
536 | //
|
---|
537 | // Now setup the tasks and tasklist:
|
---|
538 | // ---------------------------------
|
---|
539 | //
|
---|
540 |
|
---|
541 | MParList plist;
|
---|
542 |
|
---|
543 | MTaskList tlist;
|
---|
544 | plist.AddToList(&tlist);
|
---|
545 |
|
---|
546 |
|
---|
547 | MReadMarsFile read("Events");
|
---|
548 | read.DisableAutoScheme();
|
---|
549 |
|
---|
550 | read.AddFile(filename);
|
---|
551 |
|
---|
552 | MHillas mhillas;
|
---|
553 | plist.AddToList(&mhillas);
|
---|
554 |
|
---|
555 | MHillasSrc mhillassrc;
|
---|
556 | plist.AddToList(&mhillassrc);
|
---|
557 |
|
---|
558 | MHillasExt mhillasext;
|
---|
559 | plist.AddToList(&mhillasext);
|
---|
560 |
|
---|
561 | MNewImagePar mnewimpar;
|
---|
562 | plist.AddToList(&mnewimpar);
|
---|
563 |
|
---|
564 | MGeomCamMagic cam;
|
---|
565 | plist.AddToList(&cam);
|
---|
566 |
|
---|
567 | MMcEvt mcevt;
|
---|
568 | plist.AddToList(&mcevt);
|
---|
569 |
|
---|
570 | MPointingPos mpoint;
|
---|
571 | plist.AddToList(&mpoint);
|
---|
572 |
|
---|
573 | MObservatory observ;
|
---|
574 | plist.AddToList(&observ);
|
---|
575 |
|
---|
576 | // MRawRunHeader header;
|
---|
577 | // plist.AddToList(&header);
|
---|
578 |
|
---|
579 | tlist.AddToList(&read);
|
---|
580 |
|
---|
581 | MEvtLoop evtloop;
|
---|
582 | evtloop.SetParList(&plist);
|
---|
583 |
|
---|
584 | if (!tlist.PreProcess(&plist))
|
---|
585 | return;
|
---|
586 |
|
---|
587 | Float_t fMm2Deg = cam->GetConvMm2Deg();
|
---|
588 | Int_t event = 0;
|
---|
589 | Int_t filenumber = 0;
|
---|
590 |
|
---|
591 | Float_t ftheta, fphi, flength, fwidth, fsize, fmeanx, fmeany, falpha, fdist;
|
---|
592 | Float_t fsingam, fcosgam;
|
---|
593 | Double_t xsournew, ysournew;
|
---|
594 | Float_t fdelta, fleak, fconc1, fcosda, fassym, fassymM3;
|
---|
595 | Int_t AsGrNull=0, AsLessNull=0;
|
---|
596 | Int_t AsGrNullAfter=0, AsLessNullAfter=0;
|
---|
597 | Float_t logsize, lgsize, lgsize2, tanbeta, beta;
|
---|
598 | const Float_t LOG3000 = log(3000.);
|
---|
599 | Char_t stringtriv1[80], stringlima[80], stringNex[80], stringsig[80];
|
---|
600 | Char_t stringNexOnOff[80], stringLiMaOnOff[80];
|
---|
601 |
|
---|
602 | // initial values:
|
---|
603 | Float_t xsource = XSOURCE;
|
---|
604 | Float_t ysource = YSOURCE;
|
---|
605 |
|
---|
606 |
|
---|
607 | while (tlist.Process())
|
---|
608 | {
|
---|
609 | event++;
|
---|
610 |
|
---|
611 | if (mhillas->GetLength() != -1.)
|
---|
612 | {
|
---|
613 | // parameters:
|
---|
614 | flength = (mhillas->GetLength()) * fMm2Deg;
|
---|
615 | fwidth = (mhillas->GetWidth())*fMm2Deg;
|
---|
616 | fsize = mhillas->GetSize();
|
---|
617 | fmeanx = (mhillas->GetMeanX())*fMm2Deg;
|
---|
618 | fmeany = (mhillas->GetMeanY())*fMm2Deg;
|
---|
619 | falpha = mhillassrc->GetAlpha();
|
---|
620 | fdist = (mhillassrc->GetDist())*fMm2Deg;
|
---|
621 | fdelta = mhillas->GetDelta();
|
---|
622 | fconc1 = (mnewimpar->GetConc1());
|
---|
623 | fleak = mnewimpar->GetLeakage1();
|
---|
624 |
|
---|
625 |
|
---|
626 | // ftheta = mcevt->GetTelescopeTheta();
|
---|
627 | ftheta = mpoint->GetZd();
|
---|
628 | // fphi = mcevt->GetTelescopePhi();
|
---|
629 | fphi = mpoint->GetAz();
|
---|
630 | // cout << " phi : " << fphi << " theta : " << ftheta << endl;
|
---|
631 | observ.RotationAngle(ftheta, fphi, fsingam, fcosgam);
|
---|
632 |
|
---|
633 | fassym = (mhillasext->GetAsym()) * fMm2Deg;
|
---|
634 | fassymM3 = (mhillasext->GetM3Long()) * fMm2Deg;
|
---|
635 | fcosda = mhillassrc->GetCosDeltaAlpha();
|
---|
636 |
|
---|
637 | if ((fassymM3*TMath::Sign(1.,fcosda)) > 0.) AsGrNull++;
|
---|
638 | else AsLessNull++;
|
---|
639 |
|
---|
640 | if (ROTOPTION == kTRUE) // derotate into sky coordinates
|
---|
641 | {
|
---|
642 | /* derotation : correct sky coordinates into camera coordinates */
|
---|
643 | xsournew = fcosgam * xsource - fsingam * ysource;
|
---|
644 | ysournew = fsingam * xsource + fcosgam * ysource;
|
---|
645 | /* end derotatiom */
|
---|
646 | }
|
---|
647 | else // do not derotate, plot into camera coordinates
|
---|
648 | {
|
---|
649 | xsournew = xsource;
|
---|
650 | ysournew = ysource;
|
---|
651 | }
|
---|
652 |
|
---|
653 | // basic plots:
|
---|
654 |
|
---|
655 | // if (fsize > 3000.)
|
---|
656 | if (fsize > 0.)
|
---|
657 | {
|
---|
658 | histphi.Fill(fphi,1.);
|
---|
659 | histtheta.Fill(ftheta,1.);
|
---|
660 |
|
---|
661 | histlength.Fill(flength,1.);
|
---|
662 | histwidth.Fill(fwidth,1.);
|
---|
663 | histsize.Fill(fsize,1.);
|
---|
664 | histLoverS.Fill(flength/fsize,1.);
|
---|
665 | histmeanx.Fill(fmeanx,1.);
|
---|
666 | histmeany.Fill(fmeany,1.);
|
---|
667 | histalpha.Fill(falpha,1.);
|
---|
668 | histdist.Fill(fdist,1.);
|
---|
669 | hist2xy.Fill(fmeanx, fmeany, 1.);
|
---|
670 | }
|
---|
671 |
|
---|
672 | // some cuts:
|
---|
673 | if (flength > 0.1 && flength < 0.32)
|
---|
674 | if (fwidth > 0.06 && fwidth < 0.15)
|
---|
675 | if (fdist > 0.6 && fdist < 1.3)
|
---|
676 | if (fsize > 3000.)
|
---|
677 | // if(sqrt(fmeanx*fmeanx + fmeany*fmeany) < 1.1)
|
---|
678 | // if((fassymM3*fcosda < 0.3 && fassymM3*fcosda > 0.02) ||
|
---|
679 | // (fassymM3*fcosda < -0.02 && fassymM3*fcosda > -0.2) )
|
---|
680 | // if(fassym*fcosda > 0.)
|
---|
681 | {
|
---|
682 | histAssymM3.Fill(fassymM3*TMath::Sign(1.,fcosda), 1.);
|
---|
683 | histAssym.Fill(fassym*TMath::Sign(1.,fcosda), 1.);
|
---|
684 | if ((fassymM3*TMath::Sign(1.,fcosda)) > 0.) AsGrNullAfter++;
|
---|
685 | else AsLessNullAfter++;
|
---|
686 | }
|
---|
687 |
|
---|
688 | // ********************************************************************** //
|
---|
689 | // calculate alpha and dist according to the source location:
|
---|
690 | tanbeta = (fmeany - ysournew) / (fmeanx - xsournew);
|
---|
691 | beta = TMath::ATan(tanbeta);
|
---|
692 | falpha = (fdelta - beta) * 180./ TMath::Pi();
|
---|
693 | fdist = sqrt((fmeany - ysournew) * (fmeany - ysournew) +
|
---|
694 | (fmeanx - xsournew) * (fmeanx - xsournew));
|
---|
695 |
|
---|
696 | if(falpha > 90.) falpha -= 180.;
|
---|
697 | if(falpha < -90.) falpha += 180.;
|
---|
698 |
|
---|
699 | // ********************************************************************** //
|
---|
700 |
|
---|
701 |
|
---|
702 |
|
---|
703 | if (fsize > 3000.)
|
---|
704 | {
|
---|
705 | // correlations:
|
---|
706 | hist2lw.Fill(flength, fwidth, 1.);
|
---|
707 | hist2lalpha.Fill(flength, falpha, 1.);
|
---|
708 | hist2ldist.Fill(flength, fdist, 1.);
|
---|
709 | hist2walpha.Fill(fwidth, falpha, 1.);
|
---|
710 | hist2wdist.Fill(fwidth, fdist, 1.);
|
---|
711 | hist2alphadist.Fill(falpha, fdist, 1.);
|
---|
712 | hist2wsize.Fill(fwidth, fsize, 1.);
|
---|
713 | hist2alphasize.Fill(falpha, fsize, 1.);
|
---|
714 | hist2distsize.Fill(fdist, fsize, 1.);
|
---|
715 | }
|
---|
716 |
|
---|
717 | // cuts:
|
---|
718 | //cout << " before the cuts" << "size :" << fsize << endl;
|
---|
719 | logsize = log(fsize);
|
---|
720 | lgsize = logsize-LOG3000;
|
---|
721 | lgsize2 = lgsize*lgsize;
|
---|
722 | if ( fsize > SIZEMIN )
|
---|
723 | if ( fleak < LEAKMAX )
|
---|
724 | if ( flength > (LENGTHMINParA + LENGTHMINParB*lgsize + LENGTHMINParC*lgsize2) &&
|
---|
725 | flength < (LENGTHMAXParA + LENGTHMAXParB*lgsize + LENGTHMAXParC*lgsize2))
|
---|
726 | if ( fwidth > (WIDTHMINParA + WIDTHMINParB*lgsize + WIDTHMINParC*lgsize2) &&
|
---|
727 | fwidth < (WIDTHMAXParA + WIDTHMAXParB*lgsize + WIDTHMAXParC*lgsize2) )
|
---|
728 | if ( fdist > (DISTMINParA + DISTMINParB*lgsize + DISTMINParC*lgsize2) &&
|
---|
729 | fdist < (DISTMAXParA + DISTMAXParB*lgsize + DISTMAXParC*lgsize2) )
|
---|
730 | // if ((fassym*TMath::Sign(1.,fcosda)) > ASYMMIN) // asymmcut
|
---|
731 | {
|
---|
732 | falpha = TMath::Abs(falpha);
|
---|
733 | histalphafinal.Fill(falpha,1.);
|
---|
734 | histalphaAll.Fill(falpha,1.);
|
---|
735 |
|
---|
736 | aftercuthistlength.Fill(flength,1.);
|
---|
737 | aftercuthistwidth.Fill(fwidth,1.);
|
---|
738 | aftercuthistsize.Fill(fsize,1.);
|
---|
739 | aftercuthistmeanx.Fill(fmeanx,1.);
|
---|
740 | aftercuthistmeany.Fill(fmeany,1.);
|
---|
741 | aftercuthistalpha.Fill(falpha,1.);
|
---|
742 | aftercuthistdist.Fill(fdist,1.);
|
---|
743 | aftercuthist2xy.Fill(fmeanx, fmeany, 1.);
|
---|
744 |
|
---|
745 | }
|
---|
746 |
|
---|
747 |
|
---|
748 | }
|
---|
749 | else filenumber++;
|
---|
750 | }
|
---|
751 |
|
---|
752 | // cout << " conversion factor is: " << fMm2Deg << endl;
|
---|
753 | cout << " events read in from file : " << event << endl;
|
---|
754 | cout << " runs found in the file : " << filenumber << endl;
|
---|
755 |
|
---|
756 | Int_t startbinoff;
|
---|
757 | Float_t Nex, Non, Noff, Sign, SignLiMa;
|
---|
758 | Float_t normf, integon, integoff, NexOnOff, NoffOFF, SignOnOff, SignLiMaOnOff;
|
---|
759 | Float_t binwidth = histalphafinal.GetBinWidth(1);
|
---|
760 | Float_t numbinMax = ALPHAMAX/binwidth;
|
---|
761 |
|
---|
762 | // ********************************************************************** //
|
---|
763 | /* fit parabel from 30 to 90 degrees */
|
---|
764 | TF1 * fitbgpar = new TF1("fbgpar", "[0]*x*x + [1]", ALOFFMIN, ALOFFMAX);
|
---|
765 | fitbgpar->SetLineColor(2);
|
---|
766 |
|
---|
767 | histalphafinal.Fit("fbgpar","WNR");
|
---|
768 |
|
---|
769 | Double_t apar = fitbgpar->GetParameter(0);
|
---|
770 | Double_t bpar = fitbgpar->GetParameter(1);
|
---|
771 |
|
---|
772 | TF1 * bgoff = new TF1("bgoffON", parabfunc, 0., 90., 3);
|
---|
773 | bgoff->SetParameters(apar, bpar, 1.);
|
---|
774 | bgoff->FixParameter(0, apar);
|
---|
775 | bgoff->FixParameter(1, bpar);
|
---|
776 | bgoff->FixParameter(2, 1.);
|
---|
777 | bgoff->SetLineColor(9);
|
---|
778 |
|
---|
779 | /* end of the fit parabel from 30 to 90 degrees*/
|
---|
780 | // ********************************************************************** //
|
---|
781 |
|
---|
782 |
|
---|
783 |
|
---|
784 | if (!tlist.PostProcess())
|
---|
785 | return;
|
---|
786 |
|
---|
787 |
|
---|
788 | gStyle->SetOptStat(11);
|
---|
789 |
|
---|
790 | // calculate significance: DO NOT USE FIT FOR Non!!!
|
---|
791 | Non = 0.;
|
---|
792 | for(Int_t i=1; i<=numbinMax;i++) Non += histalphafinal.GetBinContent(i);
|
---|
793 |
|
---|
794 | Noff = (1./3. * (fitbgpar->GetParameter(0)) * pow(ALPHAMAX,3.) +
|
---|
795 | (fitbgpar->GetParameter(1)) * ALPHAMAX) / binwidth;
|
---|
796 | Nex = Non - Noff;
|
---|
797 |
|
---|
798 | Sign = Nex / sqrt(Nex + 2.* Noff);
|
---|
799 |
|
---|
800 | cout << " Non : " << Non << " Noff : " << Noff << " Nex : " << Nex << endl;
|
---|
801 | cout << " significance : " << Sign << " sigma" << endl;
|
---|
802 |
|
---|
803 | SignLiMa = LiMa17(Non,Noff,1.);
|
---|
804 | cout << " significance Li and Ma (17): " << SignLiMa << " sigma" << endl;
|
---|
805 |
|
---|
806 | Char_t stringsig[80];
|
---|
807 |
|
---|
808 | sprintf(stringsig,"S = %.2f sigma", Sign);
|
---|
809 | sprintf(stringtriv1,"Signif: S = %.2f sigma", Sign);
|
---|
810 | sprintf(stringlima,"Li&Ma 17: S = %.2f sigma", SignLiMa);
|
---|
811 | sprintf(stringNex,"N excess: Nex = %.d ", Nex);
|
---|
812 |
|
---|
813 | // ********************************************************************** //
|
---|
814 | // use OFF data to estimate background ******************************* //
|
---|
815 | TF1 * bgoff2 = new TF1("bgoffOFF", parabfunc, 0., 90., 3);
|
---|
816 | if (USEOFFDATA == kTRUE)
|
---|
817 | {
|
---|
818 | // ON:
|
---|
819 |
|
---|
820 | integon = 0.; // number of events between 30 and 90 degrees
|
---|
821 | numbinsoff = TMath::Nint((ALOFFMAX - ALOFFMIN)/binwidth);
|
---|
822 | startbinoff = TMath::Nint(ALOFFMIN/binwidth) + 1;
|
---|
823 |
|
---|
824 | for (Int_t ik = 0; ik < numbinsoff; ik++)
|
---|
825 | {
|
---|
826 | integon += histalphafinal.GetBinContent(startbinoff+ik);
|
---|
827 | }
|
---|
828 | // OFF:
|
---|
829 |
|
---|
830 | integoff = ((1./3. * aparOff * pow(90.,3.) + bparOff * 90.) -
|
---|
831 | (1./3. * aparOff * pow(30.,3.) + bparOff * 30.)) / binwidthOff;
|
---|
832 |
|
---|
833 | normf = integoff / integon;
|
---|
834 |
|
---|
835 | NoffOFF = (1./3. * aparOff * pow(ALPHAMAX,3.) +
|
---|
836 | (bparOff * ALPHAMAX)) / binwidthOff / normf;
|
---|
837 |
|
---|
838 | NexOnOff = Non - NoffOFF;
|
---|
839 |
|
---|
840 | SignOnOff = NexOnOff / sqrt(NexOnOff + 2.* NoffOFF);
|
---|
841 |
|
---|
842 | // calculate according to Li Ma:
|
---|
843 | SignLiMaOnOff = LiMa17(Non,NoffOFF*normf,1./normf);
|
---|
844 |
|
---|
845 | cout << " integon: " << integon << ", integoff : " << integoff << ", normf: " << normf
|
---|
846 | << ", NoffOFF : " << NoffOFF << ", Non : " << Non
|
---|
847 | << ", NexOnOff : " << NexOnOff << ", SignOnOff : " << SignOnOff << endl;
|
---|
848 | cout << " significance (LiMa 17): " << SignLiMaOnOff << " sigma" << endl;
|
---|
849 |
|
---|
850 | sprintf(stringNexOnOff,"N excess (ON - OFF) = %.d ", NexOnOff);
|
---|
851 | sprintf(stringLiMaOnOff,"Signif (ON - OFF) = %.2f ", SignLiMaOnOff);
|
---|
852 |
|
---|
853 |
|
---|
854 | bgoff2->SetParameters(aparOff, bparOff, normf/binwidth*binwidthOff);
|
---|
855 | bgoff2->FixParameter(0, aparOff);
|
---|
856 | bgoff2->FixParameter(1, bparOff);
|
---|
857 | bgoff2->FixParameter(2, normf/binwidth*binwidthOff);
|
---|
858 | bgoff2->SetLineColor(2);
|
---|
859 |
|
---|
860 | }
|
---|
861 |
|
---|
862 | /*
|
---|
863 | TCanvas canv("c1", "basic histograms", 600, 500);
|
---|
864 | canv.SetBorderMode(0);
|
---|
865 | canv.Divide(3,3);
|
---|
866 |
|
---|
867 | canv.cd(1);
|
---|
868 | gPad->SetBorderMode(0);
|
---|
869 | histlength.Draw();
|
---|
870 |
|
---|
871 | canv.cd(2);
|
---|
872 | gPad->SetBorderMode(0);
|
---|
873 | histwidth.Draw();
|
---|
874 |
|
---|
875 | canv.cd(3);
|
---|
876 | gPad->SetBorderMode(0);
|
---|
877 | gPad->SetLogx();
|
---|
878 | gPad->SetLogy();
|
---|
879 | histsize.Draw();
|
---|
880 |
|
---|
881 | canv.cd(4);
|
---|
882 | gPad->SetBorderMode(0);
|
---|
883 | histalpha.Draw();
|
---|
884 |
|
---|
885 | canv.cd(5);
|
---|
886 | gPad->SetBorderMode(0);
|
---|
887 | histdist.Draw();
|
---|
888 |
|
---|
889 | canv.cd(6);
|
---|
890 | gPad->SetBorderMode(0);
|
---|
891 | histmeanx.Draw();
|
---|
892 |
|
---|
893 | canv.cd(7);
|
---|
894 | gPad->SetBorderMode(0);
|
---|
895 | histmeany.Draw();
|
---|
896 |
|
---|
897 | canv.cd(8);
|
---|
898 | gPad->SetBorderMode(0);
|
---|
899 | hist2xy.Draw();
|
---|
900 |
|
---|
901 | canv.cd(9);
|
---|
902 | gPad->SetBorderMode(0);
|
---|
903 | histLoverS.Draw();
|
---|
904 |
|
---|
905 | canv.Modified();
|
---|
906 | canv.Update();
|
---|
907 |
|
---|
908 | canv.DrawClone();
|
---|
909 |
|
---|
910 |
|
---|
911 |
|
---|
912 |
|
---|
913 | TCanvas canvcor("c2", "correlation histograms", 600, 500);
|
---|
914 | canvcor.SetBorderMode(0);
|
---|
915 | canvcor.Divide(3,3);
|
---|
916 |
|
---|
917 | canvcor.cd(1);
|
---|
918 | gPad->SetBorderMode(0);
|
---|
919 | hist2lw.Draw();
|
---|
920 |
|
---|
921 | canvcor.cd(2);
|
---|
922 | gPad->SetBorderMode(0);
|
---|
923 | hist2lalpha.Draw();
|
---|
924 |
|
---|
925 | canvcor.cd(3);
|
---|
926 | gPad->SetBorderMode(0);
|
---|
927 | hist2ldist.Draw();
|
---|
928 |
|
---|
929 | canvcor.cd(4);
|
---|
930 | gPad->SetBorderMode(0);
|
---|
931 | hist2walpha.Draw();
|
---|
932 |
|
---|
933 | canvcor.cd(5);
|
---|
934 | gPad->SetBorderMode(0);
|
---|
935 | hist2wdist.Draw();
|
---|
936 |
|
---|
937 | canvcor.cd(6);
|
---|
938 | gPad->SetBorderMode(0);
|
---|
939 | hist2alphadist.Draw();
|
---|
940 |
|
---|
941 | canvcor.cd(7);
|
---|
942 | gPad->SetBorderMode(0);
|
---|
943 | gPad->SetLogy();
|
---|
944 | hist2wsize.Draw();
|
---|
945 |
|
---|
946 |
|
---|
947 | canvcor.cd(8);
|
---|
948 | gPad->SetBorderMode(0);
|
---|
949 | gPad->SetLogy();
|
---|
950 | hist2alphasize.Draw();
|
---|
951 |
|
---|
952 | canvcor.cd(9);
|
---|
953 | gPad->SetBorderMode(0);
|
---|
954 | gPad->SetLogy();
|
---|
955 | hist2distsize.Draw();
|
---|
956 |
|
---|
957 |
|
---|
958 | canvcor.Modified();
|
---|
959 | canvcor.Update();
|
---|
960 |
|
---|
961 | canvcor.DrawClone();
|
---|
962 | */
|
---|
963 |
|
---|
964 | /**********************************************************/
|
---|
965 | /* plot the alpha plot for the current sample */
|
---|
966 |
|
---|
967 | Char_t titelname[80];
|
---|
968 | sprintf(titelname,"alpha plot. sample %d assumed source position: x = %.2f y = %.2f", num, XSOURCE, YSOURCE);
|
---|
969 |
|
---|
970 |
|
---|
971 | TCanvas canval("c3", "canvas for alpha", 600, 500);
|
---|
972 | canval.SetBorderMode(0);
|
---|
973 |
|
---|
974 | gPad->SetBorderMode(0);
|
---|
975 | histalphafinal.SetMarkerStyle(20);
|
---|
976 | histalphafinal.SetTitle(titelname);
|
---|
977 | histalphafinal.SetFillColor(8);
|
---|
978 | histalphafinal.Draw();
|
---|
979 | bgoff->Draw("same");
|
---|
980 | if (USEOFFDATA == kTRUE) bgoff2->Draw("same");
|
---|
981 |
|
---|
982 | leg = new TLegend(0.1,0.15,0.52,0.35);
|
---|
983 | // leg->Draw();
|
---|
984 | leg->AddEntry(fitbgpar,"fit for OFF region (30-90 deg)","l");
|
---|
985 | leg->SetHeader("Legend");
|
---|
986 | leg->SetFillColor(19);
|
---|
987 | // leg->Draw();
|
---|
988 |
|
---|
989 | if (USEOFFDATA == kFALSE)
|
---|
990 | {
|
---|
991 | text = new TPaveText(0.53,0.45,0.9,0.65,"NDC");
|
---|
992 | text->AddText(0.4, 0.6, stringNex);
|
---|
993 | text->AddText(0.5, 0.3, stringlima);
|
---|
994 | text->SetTextSize(0.032);
|
---|
995 | text->Draw();
|
---|
996 | }
|
---|
997 | else
|
---|
998 | {
|
---|
999 | text = new TPaveText(0.53,0.25,0.9,0.55,"NDC");
|
---|
1000 | text->AddText(0.4, 0.8, stringNex);
|
---|
1001 | text->AddText(0.45, 0.6, stringlima);
|
---|
1002 | text->AddText(0.5, 0.4, stringNexOnOff);
|
---|
1003 | text->AddText(0.4, 0.2, stringLiMaOnOff);
|
---|
1004 | text->SetTextSize(0.032);
|
---|
1005 | text->Draw();
|
---|
1006 | }
|
---|
1007 | canval.Modified();
|
---|
1008 | canval.Update();
|
---|
1009 |
|
---|
1010 | // canval.DrawClone();
|
---|
1011 |
|
---|
1012 | TString strin = dirname + sourcename;
|
---|
1013 | // TString strin = "data/plots/alpha/dummySize";
|
---|
1014 | strin += SIZEMIN;
|
---|
1015 | strin += "Sample";
|
---|
1016 | if (num<10) strin +="0";
|
---|
1017 | strin += num;
|
---|
1018 | strin += ".root";
|
---|
1019 |
|
---|
1020 |
|
---|
1021 | canval.SaveAs(strin); // please enable if you want to save alphaplots
|
---|
1022 | cout << " alpha plot for the sample " << num << " has been saved into " << strin << endl;
|
---|
1023 |
|
---|
1024 | if(bgoff2) delete bgoff2;
|
---|
1025 | delete bgoff;
|
---|
1026 |
|
---|
1027 | /**********************************************************/
|
---|
1028 |
|
---|
1029 | /*
|
---|
1030 | TCanvas canvt("c4", "telescope", 600, 500);
|
---|
1031 | canvt.SetBorderMode(0);
|
---|
1032 | canvt.Divide(2,1);
|
---|
1033 |
|
---|
1034 | canvt.cd(1);
|
---|
1035 | gPad->SetBorderMode(0);
|
---|
1036 | histphi.Draw();
|
---|
1037 |
|
---|
1038 | canvt.cd(2);
|
---|
1039 | gPad->SetBorderMode(0);
|
---|
1040 | histtheta.Draw();
|
---|
1041 |
|
---|
1042 | canvt.DrawClone();
|
---|
1043 |
|
---|
1044 | TCanvas canvass("c5", "assymetry", 600, 500);
|
---|
1045 | canvass.SetBorderMode(0);
|
---|
1046 | canvass.Divide(2,1);
|
---|
1047 |
|
---|
1048 | canvass.cd(1);
|
---|
1049 | gPad->SetGridx();
|
---|
1050 | gPad->SetGridy();
|
---|
1051 | gPad->SetBorderMode(0);
|
---|
1052 | histAssym.Draw();
|
---|
1053 |
|
---|
1054 | canvass.cd(2);
|
---|
1055 | gPad->SetGridx();
|
---|
1056 | gPad->SetGridy();
|
---|
1057 | gPad->SetBorderMode(0);
|
---|
1058 | histAssymM3.Draw();
|
---|
1059 |
|
---|
1060 | canvass.DrawClone();
|
---|
1061 |
|
---|
1062 | TCanvas aftercutcanv("c1a", "basic histograms", 600, 500);
|
---|
1063 | aftercutcanv.SetBorderMode(0);
|
---|
1064 | aftercutcanv.Divide(3,3);
|
---|
1065 |
|
---|
1066 | aftercutcanv.cd(1);
|
---|
1067 | gPad->SetBorderMode(0);
|
---|
1068 | gPad->SetGridx();
|
---|
1069 | gPad->SetGridy();
|
---|
1070 | aftercuthistlength.Draw();
|
---|
1071 |
|
---|
1072 | aftercutcanv.cd(2);
|
---|
1073 | gPad->SetBorderMode(0);
|
---|
1074 | gPad->SetGridx();
|
---|
1075 | gPad->SetGridy();
|
---|
1076 | aftercuthistwidth.Draw();
|
---|
1077 |
|
---|
1078 | aftercutcanv.cd(3);
|
---|
1079 | gPad->SetBorderMode(0);
|
---|
1080 | gPad->SetLogx();
|
---|
1081 | gPad->SetLogy();
|
---|
1082 | gPad->SetGridx();
|
---|
1083 | gPad->SetGridy();
|
---|
1084 | aftercuthistsize.Draw();
|
---|
1085 |
|
---|
1086 | aftercutcanv.cd(4);
|
---|
1087 | gPad->SetBorderMode(0);
|
---|
1088 | gPad->SetGridx();
|
---|
1089 | gPad->SetGridy();
|
---|
1090 | aftercuthistalpha.Draw();
|
---|
1091 |
|
---|
1092 | aftercutcanv.cd(5);
|
---|
1093 | gPad->SetBorderMode(0);
|
---|
1094 | gPad->SetGridx();
|
---|
1095 | gPad->SetGridy();
|
---|
1096 | aftercuthistdist.Draw();
|
---|
1097 |
|
---|
1098 | aftercutcanv.cd(6);
|
---|
1099 | gPad->SetBorderMode(0);
|
---|
1100 | gPad->SetGridx();
|
---|
1101 | gPad->SetGridy();
|
---|
1102 | aftercuthistmeanx.Draw();
|
---|
1103 |
|
---|
1104 | aftercutcanv.cd(7);
|
---|
1105 | gPad->SetBorderMode(0);
|
---|
1106 | gPad->SetGridx();
|
---|
1107 | gPad->SetGridy();
|
---|
1108 | aftercuthistmeany.Draw();
|
---|
1109 |
|
---|
1110 |
|
---|
1111 | aftercutcanv.cd(8);
|
---|
1112 | gPad->SetBorderMode(0);
|
---|
1113 | gPad->SetGridx();
|
---|
1114 | gPad->SetGridy();
|
---|
1115 | aftercuthist2xy.Draw();
|
---|
1116 |
|
---|
1117 | aftercutcanv.Modified();
|
---|
1118 | aftercutcanv.Update();
|
---|
1119 |
|
---|
1120 | aftercutcanv.DrawClone();
|
---|
1121 | */
|
---|
1122 | Double_t rat1, rat2;
|
---|
1123 | rat1 = (double)AsGrNull / (double)AsLessNull;
|
---|
1124 | rat2 = (double)AsGrNullAfter / (double)AsLessNullAfter;
|
---|
1125 |
|
---|
1126 | cout << " Asymmetry M3 > 0 : " << AsGrNull << endl;
|
---|
1127 | cout << " Asymmetry M3 < 0 : " << AsLessNull << endl;
|
---|
1128 | cout << " Ratio (before cuts) : " << rat1 << endl;
|
---|
1129 | cout << " Asymmetry M3 > 0 (after): " << AsGrNullAfter << endl;
|
---|
1130 | cout << " Asymmetry M3 < 0 (after): " << AsLessNullAfter << endl;
|
---|
1131 | cout << " Ratio (after cuts) : " << rat2 << endl;
|
---|
1132 |
|
---|
1133 |
|
---|
1134 | }
|
---|
1135 |
|
---|
1136 | void callalphaplot()
|
---|
1137 | {
|
---|
1138 |
|
---|
1139 |
|
---|
1140 | Int_t num;
|
---|
1141 | Double_t xpeakM, ypeakM, xpeakB, ypeakB;
|
---|
1142 |
|
---|
1143 | FILE *fp;
|
---|
1144 | // fp = fopen("data/trackMrk421_0505_2000.dat", "r");
|
---|
1145 | if (USEFILE == kTRUE)
|
---|
1146 | {
|
---|
1147 | fp = fopen("data/trackCrab1502_1305_Berlin_5deg.dat", "r");
|
---|
1148 |
|
---|
1149 | for(Int_t i = 0; i < 2; i++)
|
---|
1150 | {
|
---|
1151 | fscanf(fp,"%d %lf %lf %lf %lf", &num, &xpeakM, &ypeakM, &xpeakB, &ypeakB);
|
---|
1152 | cout << endl << " SUBS NUMBER " << num << ", xpeakM = " << xpeakM <<
|
---|
1153 | ", ypeakM = " << ypeakM << endl;
|
---|
1154 | cout << " xpeakB = " << xpeakB << ", ypeakB = " << ypeakB << endl;
|
---|
1155 | if (num > 0 ) alphaplots(num, xpeakM, ypeakM);
|
---|
1156 | }
|
---|
1157 |
|
---|
1158 | fclose(fp);
|
---|
1159 | cout << "FERTIG" << endl;
|
---|
1160 | }
|
---|
1161 | else
|
---|
1162 | alphaplots(NUM, XSOUR, YSOUR);
|
---|
1163 |
|
---|
1164 |
|
---|
1165 | /**********************************************************/
|
---|
1166 | /* now calculate Nex and S for the overall alpha plot */:
|
---|
1167 |
|
---|
1168 |
|
---|
1169 |
|
---|
1170 | TF1 * fitbgparAll = new TF1("fbgparA", "[0]*x*x + [1]", 30., 90.);
|
---|
1171 | fitbgparAll->SetLineColor(2);
|
---|
1172 | fitbgparAll->SetLineWidth(3);
|
---|
1173 |
|
---|
1174 | histalphaAll.Fit("fbgparA","WR");
|
---|
1175 |
|
---|
1176 | Double_t apar = fitbgparAll->GetParameter(0);
|
---|
1177 | Double_t bpar = fitbgparAll->GetParameter(1);
|
---|
1178 | Double_t normf = 1.;
|
---|
1179 |
|
---|
1180 | TF1 * bgoff = new TF1("bgoff", parabfunc, 0., 90., 3);
|
---|
1181 | bgoff->SetParameters(apar, bpar, normf);
|
---|
1182 | bgoff->FixParameter(0, apar);
|
---|
1183 | bgoff->FixParameter(1, bpar);
|
---|
1184 | bgoff->FixParameter(2, normf);
|
---|
1185 | bgoff->SetLineColor(8);
|
---|
1186 |
|
---|
1187 | // calc significance:
|
---|
1188 | Double_t Sign, Non, Noff, Nex;
|
---|
1189 | Double_t binwidth = histalphaAll.GetBinWidth(1);
|
---|
1190 | Double_t numbinMax = ALPHAMAX / binwidth;
|
---|
1191 |
|
---|
1192 | Non = 0.;
|
---|
1193 | for(Int_t i=1; i<=numbinMax;i++) Non += histalphaAll.GetBinContent(i);
|
---|
1194 |
|
---|
1195 | //cout << histalphaAll.GetBinContent(1) + histalphaAll.GetBinContent(2) +
|
---|
1196 | // histalphaAll.GetBinContent(3) + histalphaAll.GetBinContent(4) << endl;
|
---|
1197 | Noff = (1./3. * apar * pow(ALPHAMAX,3.) + bpar * ALPHAMAX) / binwidth;
|
---|
1198 | Nex = Non - Noff;
|
---|
1199 |
|
---|
1200 | Sign = LiMa17(Non,Noff,1.);
|
---|
1201 |
|
---|
1202 | cout << " Non : " << Non << " Noff : " << Noff << " Nex : " << Nex << endl;
|
---|
1203 | cout << " significance : " << Sign << " sigma" << endl;
|
---|
1204 |
|
---|
1205 |
|
---|
1206 | Char_t stringsig[80], stringNex[80];
|
---|
1207 |
|
---|
1208 | sprintf(stringsig,"Signif: S = %.2f sigma", Sign);
|
---|
1209 | sprintf(stringNex,"N excess: Nex = %.d ", Nex);
|
---|
1210 | /**********************************************************/
|
---|
1211 |
|
---|
1212 |
|
---|
1213 | // plot all alpha plots together
|
---|
1214 |
|
---|
1215 | TCanvas canvA("cA", "alphacanvas", 600, 500);
|
---|
1216 | canvA.SetBorderMode(0);
|
---|
1217 |
|
---|
1218 | gPad->SetBorderMode(0);
|
---|
1219 | histalphaAll.SetXTitle("alpha [deg]");
|
---|
1220 | histalphaAll.SetYTitle("Counts");
|
---|
1221 | histalphaAll.SetMarkerStyle(20);
|
---|
1222 | histalphaAll.SetFillColor(17);
|
---|
1223 | histalphaAll.Draw();
|
---|
1224 | bgoff->Draw("same");
|
---|
1225 |
|
---|
1226 | // leg = new TLegend(0.1,0.15,0.52,0.35);
|
---|
1227 | // leg->Draw();
|
---|
1228 | // leg->AddEntry(fitbgpar,"fit for OFF region (30-90 deg)","l");
|
---|
1229 | // leg->SetHeader("Legend");
|
---|
1230 | // leg->SetFillColor(19);
|
---|
1231 | // leg->Draw();
|
---|
1232 |
|
---|
1233 | text = new TPaveText(0.53,0.45,0.9,0.65,"NDC");
|
---|
1234 | text->AddText(0.4, 0.6, stringNex);
|
---|
1235 | text->AddText(0.45, 0.3, stringsig);
|
---|
1236 | text->SetTextSize(0.032);
|
---|
1237 | text->Draw();
|
---|
1238 |
|
---|
1239 | canvA.Modified();
|
---|
1240 | canvA.Update();
|
---|
1241 |
|
---|
1242 |
|
---|
1243 | Char_t xstr[20], ystr[20];
|
---|
1244 | if (USEFILE == kTRUE)
|
---|
1245 | {
|
---|
1246 | sprintf(xstr,"%3.2f",0.);
|
---|
1247 | sprintf(ystr,"%3.2f",0.);
|
---|
1248 | }
|
---|
1249 | else
|
---|
1250 | {
|
---|
1251 | sprintf(xstr,"%3.2f",XSOUR);
|
---|
1252 | sprintf(ystr,"%3.2f",YSOUR);
|
---|
1253 | }
|
---|
1254 |
|
---|
1255 |
|
---|
1256 | TString string = dirname + sourcename;
|
---|
1257 | string += "X";
|
---|
1258 | string += xstr;
|
---|
1259 | string += "Y";
|
---|
1260 | string += ystr;
|
---|
1261 | string += "Size";
|
---|
1262 | string += SIZEMIN;
|
---|
1263 | string += "alpha";
|
---|
1264 | string += TMath::Nint(ALPHAMAX);
|
---|
1265 | string += ".root";
|
---|
1266 |
|
---|
1267 | canvA.SaveAs(string);
|
---|
1268 | cout << " alpha plot has been saved into " << string << endl;
|
---|
1269 |
|
---|
1270 | }
|
---|
1271 |
|
---|