source: trunk/MagicSoft/Mars/macros/flux.C@ 5675

Last change on this file since 5675 was 1543, checked in by tbretz, 23 years ago
*** empty log message ***
File size: 22.5 KB
Line 
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 Cerenkov 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!
18! Author(s): Wolfgang Wittek 4/2002 <mailto:wittek@mppmu.mpg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2002
21!
22!
23\* ======================================================================== */
24//////////////////////////////////////////////////////////////////////////////
25// //
26// This macro calculates photon fluxes as a function of //
27// - energy and time //
28// - energy and Theta //
29// //
30// It is assumed that the data have been taken in the wobble mode. //
31// This means that there is only one set of data, from which both //
32// 'on' and 'off' data are constructed. //
33// //
34// Necessary input from MC : //
35// - migration matrix (E_est, E_true) as a functioin of Theta //
36// - effective collection areas as a function of energy and Theta //
37// //
38// //
39// The input from MC has to refer to the wobble mode too. //
40// //
41// The macro includes : //
42// - the calculation of Hillas parameters //
43// - the calculation of the effective on time //
44// - the unfolding of the distributions in the estimated energy E_est //
45// //
46// For testing purposes (as long as no real data, no energy estimator, //
47// no migration matrices and no effective collection areas are available) //
48// - event times are generated (according to a Poissonian with dead time)//
49// - a dummy energy estimator is provided //
50// - migration matrices are generated //
51// - dummy effective collection areas are used //
52// //
53//////////////////////////////////////////////////////////////////////////////
54void flux()
55{
56 //--------------------------------------------------------------
57 // empty lists of parameter containers and tasks
58 //
59 MParList parlist;
60 MTaskList tasklist;
61
62 //--------------------------------------------------------------
63 // Geometry information of the MAGIC camera
64 //
65 MGeomCamMagic geomcam;
66
67 //--------------------------------------------------------------
68 // Define the two source positions
69 //
70 MSrcPosCam source("Source");
71 MSrcPosCam antisrc("AntiSource");
72 source.SetXY(0, 0);
73 antisrc.SetXY(+240, 0);
74
75 //--------------------------------------------------------------
76 // Setup the binning for the histograms
77 //
78
79 //...............................................................
80 // These are NOT the binnings for the flux determination
81
82 MBinning binswidth("BinningWidth");
83 binswidth.SetEdges(100, 0, 0.5); // 100 bins from 0 to 0.5 deg
84
85 MBinning binslength("BinningLength");
86 binslength.SetEdges(100, 0, 1); // 100 bins from 0 to 1 deg
87
88 MBinning binscamera("BinningCamera");
89 binscamera.SetEdges(90, -2.25, 2.25); // 90 bins from -2.25 to 2.25 deg
90
91 MBinning binsdist("BinningDist");
92 binsdist.SetEdges(90, 0, 2.25); // 90 bins from 0 to 2.25 deg
93
94 MBinning binsalpha("BinningAlpha");
95 binsalpha.SetEdges(100, -100, 100);
96
97 MBinning binsasym("BinningAsym");
98 binsasym.SetEdges(100, -0.5, 0.5); // 100 bins from -0.5 to 0.5 deg
99
100 MBinning binsasymn("BinningAsymn");
101 binsasymn.SetEdges(100, -0.005, 0.005); // 100 bins from -0.005 to 0.005 deg
102
103 MBinning binsm3long("BinningM3Long");
104 binsm3long.SetEdges(100, -0.5, 0.5); // 100 bins from -0.5 to 0.5 deg
105
106 MBinning binsm3trans("BinningM3Trans");
107 binsm3trans.SetEdges(100, -0.5, 0.5); // 100 bins from -0.5 to 0.5 deg
108
109 //...............................................................
110 // These ARE the binnings for the flux determination
111
112 MBinning binsalphaflux("BinningAlphaFlux");
113 binsalphaflux.SetEdges(100, 0, 100);
114
115 MBinning binse("BinningE");
116 binse.SetEdgesLog(10, 10, 1e3);
117
118 MBinning binstheta("BinningTheta");
119 binstheta.SetEdges(7, -2.5, 32.5);
120
121 MBinning binstime("BinningTime");
122 binstime.SetEdges(5, 0, 100);
123
124 MBinning binsdifftime("BinningTimeDiff");
125 binsdifftime.SetEdges(50, 0, 0.1);
126
127
128 //--------------------------------------------------------------
129 // Fill list of parameter containers
130 //
131
132 parlist.AddToList(&tasklist);
133 parlist.AddToList(&geomcam);
134
135 parlist.AddToList(&source);
136 parlist.AddToList(&antisrc);
137
138 parlist.AddToList(&binswidth);
139 parlist.AddToList(&binslength);
140 parlist.AddToList(&binscamera);
141
142 parlist.AddToList(&binsdist);
143
144 parlist.AddToList(&binsalpha);
145
146 parlist.AddToList(&binsasym);
147 parlist.AddToList(&binsasymn);
148 parlist.AddToList(&binsm3long);
149 parlist.AddToList(&binsm3trans);
150
151 parlist.AddToList(&binsalphaflux);
152 parlist.AddToList(&binse);
153 parlist.AddToList(&binstheta);
154 parlist.AddToList(&binstime);
155 parlist.AddToList(&binsdifftime);
156
157 //--------------------------------------------------------------
158 // Setup the tasks
159 //
160
161 //.....................................
162 // input file
163 //
164 // next statement commented out
165 // MReadMarsFile reader("Events", "~/data/Gamma*.root");
166
167 MReadMarsFile reader("Events");
168
169 //reader.AddFile("/hd31/rkb/Camera/mpi/Gamma*.root");
170 //reader.AddFile("/hd31/rkb/Camera/mpi/prot_N_*.root");
171 //reader.AddFile("/hd31/rkb/Camera/mpi/prot_NS_*.root");
172 //reader.AddFile("protondata/gamma_10_n*.root");
173 reader.AddFile("protondata/gamma_10_cn*.root");
174
175 reader.EnableBranch("MMcEvt.fTheta");
176
177 //.....................................
178 // filters
179 MFTriggerLvl1 lvl1;
180
181 //.....................................
182 // generation of event time
183 MMcTimeGenerate rand;
184
185 //.....................................
186 MMcPedestalCopy pcopy;
187 MMcPedestalNSBAdd pnsb;
188 MCerPhotCalc ncalc;
189 MImgCleanStd clean;
190 MBlindPixelCalc blind;
191
192 //.....................................
193 // source independent image parameters
194 MHillasCalc hcalc;
195 MHillasExt hext;
196 parlist.AddToList(&hext);
197
198 //.....................................
199 // source dependent image parameters
200 MHillasSrcCalc hsrc1("Source", "HillasSrc");
201 MHillasSrcCalc hsrc2("AntiSource", "HillasAntiSrc");
202
203 //.....................................
204 // energy estimation
205 MEnergyEstimate estim;
206
207 //.....................................
208 // migration matrix (MC)
209 MFillH eestetrue("MHMcEnergyMigration", "HillasSrc");
210 eestetrue.SetTitle("Task to determine the migration matrix for the energy");
211
212 //.....................................
213 // plots of source independent parameters (length, width, delta, size, center)
214 MFillH hfill1h("MHHillas", "MHillas");
215 hfill1h.SetTitle("Task to plot the source independent image parameters");
216 hfill1h.SetFilter(&lvl1);
217
218 MFillH hfill1x("MHHillasExt", "MHillas");
219 hfill1x.SetTitle("Task to plot the extended image parameters");
220 hfill1x.SetFilter(&lvl1);
221
222 //.....................................
223 // plots of source dependent parameters (alpha, dist)
224 MFillH hfill2s("HSource [MHHillasSrc]", "HillasSrc");
225 hfill2s.SetTitle("Task to plot the source dependent image parameters (Source)");
226 hfill2s.SetFilter(&lvl1);
227 MFillH hfill2a("HAntiSource [MHHillasSrc]", "HillasAntiSrc");
228 hfill2a.SetTitle("Task to plot the source dependent image parameters (AntiSource)");
229 hfill2a.SetFilter(&lvl1);
230
231 //.....................................
232 // star map
233 MFillH hfill1m("MHStarMap", "MHillas");
234 hfill1m.SetTitle("Task to plot the Star map");
235
236 //.....................................
237 const Float_t alpha0 = 15; // [deg]
238 MFAlpha fsrc ("HillasSrc", '>', alpha0);
239 fsrc.SetTitle("Task to check alphaSrc > alpha0");
240
241 MFAlpha fasrc("HillasAntiSrc", '>', alpha0);
242 fasrc.SetTitle("Task to check alphaAntiSrc > alpha0");
243
244 //.....................................
245 // 1-D profile plots (<Theta>, time)
246 // (<Theta>, Theta)
247 MFillH fillthetabartime ("ThetabarTime [MHThetabarTime]", "MMcEvt");
248 fillthetabartime.SetTitle("Task to plot <Theta> vs time");
249 // fillthetabartime.SetFilter(&lvl1);
250
251 MFillH fillthetabartheta("ThetabarTheta [MHThetabarTheta]", "MMcEvt");
252 fillthetabartheta.SetTitle("Task to plot <Theta> vs theta");
253 // fillthetabartheta.SetFilter(&lvl1);
254
255 //.....................................
256 // 2-D plots (Delta(time), time)
257 // (Delta(time), Theta)
258 MFillH filldtimetime ("EffOnTime [MHTimeDiffTime]", "MMcEvt");
259 filldtimetime.SetTitle("Task to plot Delta(time) vs time");
260 // filldtime.SetFilter(&lvl1);
261
262 MFillH filldtimetheta("EffOnTheta [MHTimeDiffTheta]", "MMcEvt");
263 filldtimetheta.SetTitle("Task to plot Delta(time) vs theta");
264 // fillontheta.SetFilter(&lvl1);
265
266 //.....................................
267 // 3-D plots (alpha, E_est, time)
268 MFillH fillsptime ("SrcTime [MHAlphaEnergyTime]", "HillasSrc");
269 fillsptime.SetTitle("Task for 3D plots (alpha, E_est, time) (Source)");
270 fillsptime.SetFilter(&fasrc);
271
272 MFillH fillasptime ("ASrcTime [MHAlphaEnergyTime]", "HillasAntiSrc");
273 fillasptime.SetTitle("Task for 3D plots (alpha, E_est, time) (AntiSource)");
274 fillasptime.SetFilter(&fsrc);
275
276 //.....................................
277 // 3-D plots (alpha, E_est, Theta)
278 MFillH fillsptheta ("SrcTheta [MHAlphaEnergyTheta]", "HillasSrc");
279 fillsptheta.SetTitle("Task for 3D plots (alpha, E_est, Theta) (Source)");
280 fillsptheta.SetFilter(&fasrc);
281
282 MFillH fillasptheta("ASrcTheta [MHAlphaEnergyTheta]", "HillasAntiSrc");
283 fillasptheta.SetTitle("Task for 3D plots (alpha, E_est, Theta) (AntiSource)");
284 fillasptheta.SetFilter(&fsrc);
285
286 //.....................................
287 // MFillH fethetaall("AllTheta [MHEnergyTheta]", "MMcEvt");
288 // MFillH fethetasel("SelTheta [MHEnergyTheta]", "MMcEvt");
289 // fethetasel.SetFilter(&lvl1);
290
291 // MFillH fetimeall ("AllTime [MHEnergyTime]", "MMcEvt");
292 // MFillH fetimesel ("SelTime [MHEnergyTime]", "MMcEvt");
293 // fetimesel.SetFilter(&lvl1);
294
295
296 //---------------------------------------------------------------------
297 // Setup the task list
298 //
299 tasklist.AddToList(&reader);
300
301 tasklist.AddToList(&lvl1);
302 tasklist.AddToList(&rand);
303
304 tasklist.AddToList(&pcopy);
305 tasklist.AddToList(&pnsb);
306 tasklist.AddToList(&ncalc);
307 tasklist.AddToList(&clean);
308 tasklist.AddToList(&blind);
309
310 tasklist.AddToList(&hcalc);
311
312 tasklist.AddToList(&hsrc1);
313 tasklist.AddToList(&hsrc2);
314 tasklist.AddToList(&estim);
315 tasklist.AddToList(&eestetrue);
316
317 tasklist.AddToList(&hfill1h);
318 tasklist.AddToList(&hfill1x);
319
320 tasklist.AddToList(&hfill1m);
321 tasklist.AddToList(&hfill2s);
322 tasklist.AddToList(&hfill2a);
323
324 tasklist.AddToList(&fillthetabartime);
325 tasklist.AddToList(&fillthetabartheta);
326
327 tasklist.AddToList(&filldtimetime);
328 tasklist.AddToList(&filldtimetheta);
329
330
331
332 // tasklist.AddToList(&fethetaall);
333 // tasklist.AddToList(&fethetasel);
334 // tasklist.AddToList(&fetimeall);
335 // tasklist.AddToList(&fetimesel);
336
337 tasklist.AddToList(&fsrc);
338 tasklist.AddToList(&fasrc);
339
340 tasklist.AddToList(&fillsptime);
341 tasklist.AddToList(&fillasptime);
342
343 tasklist.AddToList(&fillsptheta);
344 tasklist.AddToList(&fillasptheta);
345
346
347 //----------------------------------------------------------------------
348 // Event loop
349 //
350 MEvtLoop magic;
351 magic.SetParList(&parlist);
352
353 //
354 // loop over all events
355 //
356 if (!magic.Eventloop())
357 return;
358 //----------------------------------------------------------------------
359
360
361 tasklist.PrintStatistics(10);
362
363
364 MHHillas *mhhillas = (MHHillas*)parlist.FindObject("MHHillas");
365 mhhillas->SetTitle("Source indep. image parameters");
366 mhhillas->DrawClone();
367
368
369
370 MHHillasExt *mhhillasext = (MHHillasExt*)parlist.FindObject("MHHillasExt");
371 mhhillasext->SetTitle("Extended image parameters");
372 mhhillasext->DrawClone();
373
374
375 MHHillasSrc *hsource = (MHHillasSrc*)parlist.FindObject("HSource");
376 hsource->SetTitle("Source dependent image parameters (Source)");
377 hsource->DrawClone();
378
379 MHHillasSrc *hantisource = (MHHillasSrc*)parlist.FindObject("HAntiSource");
380 hantisource->SetTitle("Source dependent image parameters (AntiSource)");
381 hantisource->DrawClone();
382
383
384 /*
385 parlist.FindObject("MHStarMap")->DrawClone();
386 */
387
388 //----------------------------------------------------------------------
389 // average Theta versus time
390 // and versus Theta
391 // This is needed for selecting later the right collection area
392 // (at the right Theta) for each bin in time or Theta
393 //
394 MHThetabarTime *bartime = (MHThetabarTime*)parlist.FindObject("ThetabarTime");
395 MHThetabarTheta *bartheta = (MHThetabarTheta*)parlist.FindObject("ThetabarTheta");
396
397 /*
398 bartime->DrawClone();
399 bartheta->DrawClone();
400 */
401
402 //----------------------------------------------------------------------
403 // Effective on time
404
405 //....................................
406 // get plots of time differences for different intervals in time
407 // and plots of time differences for different intervals in Theta
408
409 MHTimeDiffTime *dtimetime = (MHTimeDiffTime*)parlist.FindObject("EffOnTime");
410 MHTimeDiffTheta *dtimetheta = (MHTimeDiffTheta*)parlist.FindObject("EffOnTheta");
411
412 /*
413 dtimetime->DrawClone();
414 dtimetheta->DrawClone();
415 */
416
417
418 //....................................
419 // calculate effective on time for different intervals in Time
420 // and for different intervals in Theta
421
422 MHEffOnTime effontime ("Time", "[s]");
423 MHEffOnTime effontheta("Theta", "[\\circ]");
424
425 effontime.SetupFill(&parlist);
426 effontheta.SetupFill(&parlist);
427
428
429 // Draw == 0 don't draw the individual distributions of time differences
430 // != 0 draw them
431 Bool_t draw=kFALSE;
432 effontime.Calc (dtimetime->GetHist(), draw);
433 effontheta.Calc(dtimetheta->GetHist(),draw);
434
435
436 // plot effective on time versus Time
437 // and versus Theta
438 /*
439 effontime.DrawClone();
440 effontheta.DrawClone();
441 */
442
443
444 //======================================================================
445 //
446 // A : fully differential plots (Alpha, E-est, Var)
447 //
448 //----------------------------------------------------------------------
449 // 3D-plots of alpha, E_est and time
450 // and of alpha, E_est and Theta
451 // with alpha calculated relative to the source position
452 // and relative to the anti-source position
453
454 MHAlphaEnergyTime *evtsptime = (MHAlphaEnergyTime*)parlist.FindObject("SrcTime");
455 MHAlphaEnergyTime *evtasptime = (MHAlphaEnergyTime*)parlist.FindObject("ASrcTime");
456 MHAlphaEnergyTheta *evtsptheta = (MHAlphaEnergyTheta*)parlist.FindObject("SrcTheta");
457 MHAlphaEnergyTheta *evtasptheta = (MHAlphaEnergyTheta*)parlist.FindObject("ASrcTheta");
458
459
460 /*
461 evtsptime->SetTitle("Source : 3D-plot Al-En-time, alpha\_{asp} .gt. alpha\_0");
462 evtsptime->DrawClone();
463 evtasptime->SetTitle("AntiSource : 3D-plot Al-En-time, alpha\_{sp} .gt. alpha\_0");
464 evtasptime->DrawClone();
465
466 evtsptheta->SetTitle("Source : 3D-plot Al-En-Theta, alpha\_{asp} .gt. alpha\_0");
467 evtsptheta->DrawClone();
468 evtasptheta->SetTitle("AntiSource : 3D-plot Al-En-Theta, alpha\_{sp} .gt. alpha\_0");
469 evtasptheta->DrawClone();
470 */
471
472 //----------------------------------------------------------------------
473 // Difference between source and antisource (= 'gamma' sample)
474 // for 3D-plots of alpha, E_est and time
475 // and of alpha, E_est and Theta
476 // 3rd and 4th argument are name and title of the resulting histogram
477
478 MHGamma gamma;
479
480 TH3D *hsubtime = gamma.Subtract( evtsptime->GetHist(),
481 evtasptime->GetHist(),
482 "Al-En-time", "3D-plot of Alpha,E-est,time", draw);
483
484 TH3D *hsubtheta = gamma.Subtract( evtsptheta->GetHist(),
485 evtasptheta->GetHist(),
486 "Al-En-time", "3D-plot of Alpha,E-est,Theta", draw);
487
488
489 //----------------------------------------------------------------------
490 // get number of 'gammas' in the alpha range (lo, up) as a function of
491 // E_est and time
492 // or E_est and Theta
493
494 Axis_t lo = 0; // [deg]
495 Axis_t up = 10; // [deg]
496 const TH2D &evttime = *(gamma.GetAlphaProjection(hsubtime,
497 lo, up, draw));
498
499 const TH2D &evttheta = *(gamma.GetAlphaProjection(hsubtheta,
500 lo, up, draw));
501
502
503 //----------------------------------------------------------------------
504 // plot migration matrix: E-true E-est Theta
505
506
507 /*
508 MHMcEnergyMigration *migrationmatrix =
509 (MHMcEnergyMigration*)parlist.FindObject("MHMcEnergyMigration");
510 migrationmatrix->SetTitle("Migration matrix");;
511 migrationmatrix->DrawClone();
512 */
513
514
515 //----------------------------------------------------------------------
516 // determine gamma flux from distributions of E-est :
517 // - do unfolding to get the distributions in E-unfold
518 // - divide by bin width in energy
519 // effective ontime and
520 // effective collection area
521 // to get the absolute gamma flux
522
523 //.............................................
524 // get flux spectrum for different bins in Time
525 //
526 // use dummy histogram *aeff which eventually will contain
527 // the effective collection area as a function of Etrue and Theta
528 TH2D aeff;
529
530 // arguments of MHFlux constructor are :
531 // - 2D histogram containing the no.of gammas as a function of
532 // E-est and Var
533 // - name of variable Var ("Time" or "Theta")
534 // - units of variable Var ( "[s]" or "[\\circ]")
535
536 // Draw == kTRUE draw the no.of photons vs. E-est
537 // for the individual bins of the variable Var
538 draw=kTRUE;
539 MHFlux fluxtime(evttime, draw, "Time", "[s]");
540 fluxtime.Unfold(draw);
541 fluxtime.CalcFlux(effontime.GetHist(), bartime.GetHist(),
542 &aeff, draw);
543
544 fluxtime.DrawClone();
545
546
547
548 //..............................................
549 // get flux spectrum for different bins in Theta
550 MHFlux fluxtheta(evttheta, draw, "Theta", "[\\circ]");
551 fluxtheta.Unfold(draw);
552 fluxtheta.CalcFlux(effontheta.GetHist(), bartheta.GetHist(),
553 &aeff, draw);
554
555 fluxtheta.DrawClone();
556
557
558return;
559 //----------------------------------------------------------------------
560
561
562 //======================================================================
563 //
564 // B : plots integrated over E-est
565 //
566 //----------------------------------------------------------------------
567 // integrate 3D-plot of (alpha, E_est and Time) over E-est
568 // to get 2D-plot of (alpha, Time)
569
570 Draw = kFALSE;
571 TH2D *evttimesp = evtsptime.IntegrateEest ("AlphaSP vs. Time", Draw);
572 TH2D *evttimeasp = evtasptime.IntegrateEest("AlphaASP vs. Time",Draw);
573
574
575 //======================================================================
576 //
577 // C : plots integrated over Time
578 //
579 //----------------------------------------------------------------------
580 // integrate 3D-plot of (alpha, E_est and Time) over the Time
581 // to get 2D-plot of (alpha, E_est)
582
583 Draw = kFALSE;
584 TH2D *evteestsp = evtsptime.IntegrateTime ("AlphaSP vs. E-est", Draw);
585 TH2D *evteestasp = evtasptime.IntegrateTime("AlphaASP vs. E-est",Draw);
586
587
588
589 //======================================================================
590 //
591 // D : plots integrated over E-est and Time
592 //
593
594 //----------------------------------------------------------------------
595 // integrate 3D-plot of (alpha, E_est and Time) over E-est and Time
596 // to get 1D-plot of (alpha)
597
598 Draw = kFALSE;
599 TH1D *evtalphasp = evtsptime.IntegrateEestTime ("AlphaSP", Draw);
600 TH1D *evtalphaasp = evtasptime.IntegrateEestTime("AlphaASP",Draw);
601
602 //----------------------------------------------------------------------
603 // list of tasks for processing
604
605 cout << " " << endl;
606 cout << "List of tasks for processing :" << endl;
607 cout << "==============================" << endl;
608 cout << " " << endl;
609
610 // was folgt geht nicht :
611
612 // Error: Function Next() is not defined in current scope FILE:macros/flux.C LINE:636
613
614 //while ( (task=(MTask*)Next()) )
615 //{
616 // cout << tasklist.GetName() << " : " << tasklist.GetTitle() << endl;
617 //}
618
619 return;
620
621 //======================================================================
622}
623//===========================================================================
624
625
Note: See TracBrowser for help on using the repository browser.