source: branches/Corsika7500Compatibility/mhistmc/MHMcCT1CollectionArea.cc@ 19272

Last change on this file since 19272 was 2314, checked in by wittek, 21 years ago
*** empty log message ***
File size: 16.4 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): A. Moralejo 3/2003 <mailto:moralejo@pd.infn.it>
19!
20! Copyright: MAGIC Software Development, 2000-2003
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26// //
27// MHMcCT1CollectionArea //
28// //
29//////////////////////////////////////////////////////////////////////////////
30
31#include "MHMcCT1CollectionArea.h"
32
33#include <TH2.h>
34#include <TCanvas.h>
35
36#include "MMcEvt.hxx"
37#include "MH.h"
38#include "MBinning.h"
39#include "MParList.h"
40#include "MLog.h"
41#include "MLogManip.h"
42
43ClassImp(MHMcCT1CollectionArea);
44
45using namespace std;
46
47// --------------------------------------------------------------------------
48//
49// Creates the three necessary histograms:
50// - selected showers (input)
51// - all showers (input)
52// - collection area (result)
53//
54MHMcCT1CollectionArea::MHMcCT1CollectionArea(const char *name, const char *title) :
55 fEaxis(kLog10)
56{
57 //
58 // nbins, minEnergy, maxEnergy defaults:
59 // we set the energy range from 100 Gev to 30000 GeV (in log, 3.5 orders
60 // of magnitude) and for each order we take 10 subdivisions --> 35 xbins
61 // we set the theta range from 12.5 to 48 deg, with 6 bins (the latter
62 // choice has been done to make the bin centers as close as possible to
63 // the actual zenith angles in the CT1 MC sample).
64 //
65
66 fName = name ? name : "MHMcCT1CollectionArea";
67 fTitle = title ? title : "Collection Area vs. log10 Energy";
68
69 fHistAll = new TH2D;
70 fHistSel = new TH2D;
71 fHistCol = new TH2D;
72
73 fHistCol->SetName(fName);
74 fHistAll->SetName("AllEvents");
75 fHistSel->SetName("SelectedEvents");
76
77 fHistCol->SetTitle(fTitle);
78 fHistAll->SetTitle("All showers - Theta vs log10 Energy distribution");
79 fHistSel->SetTitle("Selected showers - Theta vs log10 Energy distribution");
80
81 fHistAll->SetDirectory(NULL);
82 fHistSel->SetDirectory(NULL);
83 fHistCol->SetDirectory(NULL);
84
85 fHistAll->SetXTitle("log10(E [GeV])");
86 fHistAll->SetYTitle("\\Theta [\\circ]");
87 fHistAll->SetZTitle("Counts");
88
89 fHistSel->SetXTitle("log10(E [GeV])");
90 fHistSel->SetYTitle("\\Theta [\\circ]");
91 fHistSel->SetZTitle("Counts");
92
93 fHistCol->SetXTitle("log10(E [GeV])");
94 fHistCol->SetYTitle("theta [deg]");
95 fHistCol->SetZTitle("A [m^{2}]");
96
97}
98
99// --------------------------------------------------------------------------
100//
101// Delete the three histograms
102//
103MHMcCT1CollectionArea::~MHMcCT1CollectionArea()
104{
105 delete fHistAll;
106 delete fHistSel;
107 delete fHistCol;
108}
109
110// --------------------------------------------------------------------------
111//
112// Set the binnings and prepare the filling of the histograms
113//
114Bool_t MHMcCT1CollectionArea::SetupFill(const MParList *plist)
115{
116 const MBinning* binsenergy = (MBinning*)plist->FindObject("BinningE");
117 const MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta");
118
119 if (!binsenergy || !binstheta)
120 {
121 *fLog << err << dbginf << "At least one MBinning not found... aborting.";
122 *fLog << endl;
123 return kFALSE;
124 }
125
126 SetBinning(fHistAll, binsenergy, binstheta);
127 SetBinning(fHistSel, binsenergy, binstheta);
128 SetBinning(fHistCol, binsenergy, binstheta);
129
130 fHistAll->Sumw2();
131 fHistSel->Sumw2();
132 fHistCol->Sumw2();
133
134 if (fEaxis == kLinear)
135 {
136 fTitle = "Collection Area vs. Energy";
137 fHistCol->SetTitle(fTitle);
138 fHistAll->SetTitle("All showers - Theta vs Energy distribution");
139 fHistSel->SetTitle("Selected showers - Theta vs Energy distribution");
140 fHistCol->SetXTitle("E [GeV]");
141 fHistAll->SetXTitle("E [GeV]");
142 fHistSel->SetXTitle("E [GeV]");
143 }
144
145 return kTRUE;
146}
147
148
149// --------------------------------------------------------------------------
150//
151// Fill data into the histogram which contains the selected showers
152//
153Bool_t MHMcCT1CollectionArea::Fill(const MParContainer *par, const Stat_t w)
154{
155 MMcEvt &mcevt = *(MMcEvt*)par;
156
157 const Double_t E = fEaxis==kLinear ? mcevt.GetEnergy() : log10(mcevt.GetEnergy());
158
159 fHistSel->Fill(E, kRad2Deg*mcevt.GetTelescopeTheta(), w);
160
161 return kTRUE;
162}
163
164// --------------------------------------------------------------------------
165//
166// Draw the histogram with all showers
167//
168void MHMcCT1CollectionArea::DrawAll(Option_t* option)
169{
170 if (!gPad)
171 MH::MakeDefCanvas(fHistAll);
172
173 fHistAll->Draw(option);
174
175 gPad->Modified();
176 gPad->Update();
177}
178
179// --------------------------------------------------------------------------
180//
181// Draw the histogram with the selected showers only.
182//
183void MHMcCT1CollectionArea::DrawSel(Option_t* option)
184{
185 if (!gPad)
186 MH::MakeDefCanvas(fHistSel);
187
188 fHistSel->Draw(option);
189
190 gPad->Modified();
191 gPad->Update();
192}
193
194// --------------------------------------------------------------------------
195//
196// Creates a new canvas and draws the histogram into it.
197// Be careful: The histogram belongs to this object and won't get deleted
198// together with the canvas.
199//
200TObject *MHMcCT1CollectionArea::DrawClone(Option_t* option) const
201{
202 TCanvas &c = *MakeDefCanvas("CollArea", "Collection area plots", 600, 600);
203 c.Divide(2,2);
204
205 //
206 // This is necessary to get the expected behaviour of DrawClone
207 //
208 gROOT->SetSelectedPad(NULL);
209
210 c.cd(1);
211 if (fEaxis == kLinear)
212 gPad->SetLogx();
213 fHistCol->SetDirectory(NULL);
214 fHistCol->DrawCopy(option);
215
216 c.cd(2);
217 if (fEaxis == kLinear)
218 gPad->SetLogx();
219 fHistSel->SetDirectory(NULL);
220 fHistSel->DrawCopy(option);
221
222 c.cd(3);
223 if (fEaxis == kLinear)
224 gPad->SetLogx();
225 fHistAll->SetDirectory(NULL);
226 fHistAll->DrawCopy(option);
227
228
229 c.Modified();
230 c.Update();
231
232 return &c;
233}
234
235void MHMcCT1CollectionArea::Draw(Option_t* option)
236{
237 if (!gPad)
238 MH::MakeDefCanvas(fHistCol);
239
240 fHistCol->Draw(option);
241
242 gPad->Modified();
243 gPad->Update();
244}
245
246//
247// Calculate the Efficiency (collection area) for the CT1 MC sample
248// and set the 'ReadyToSave' flag
249//
250void MHMcCT1CollectionArea::CalcEfficiency()
251{
252 //
253 // Here we estimate the total number of showers in each energy bin
254 // from the known the energy range and spectral index of the generated
255 // showers. This procedure is intended for the CT1 MC files. The total
256 // number of generated events, collection area, spectral index etc will be
257 // set here by hand, so make sure that the MC sample you are using is the
258 // right one (check all these quantities in your files and compare with
259 // what is written below. In some theta bins, there are two different
260 // productions, with different energy limits but with the same spectral
261 // slope. We account for this when calculating the original number of
262 // events in each energy bin.
263 //
264 // The theta angle with which the MC data (from D. Kranich) were produced
265 // is not exactly the center of the theta bins we are using (the bin limits
266 // should be 0.0, 17.5, 23.5, 29.5, 35.5, 42., 50.). The theta variable in
267 // the MC root file has nevertheless been changed (W.Wittek) to correspond
268 // to the centers of these bins. Only in the first bin is the difference big:
269 // the data were produced at theta = 15 degrees, whreas the bin center is at
270 // 8.75 degrees. Howeverm at such low z.a. the shower characteristics change
271 // very slowly with theta.
272 //
273 //
274 //
275 // Only for the binning taken from D. Kranich :
276 //
277
278 for (Int_t thetabin = 1; thetabin <= fHistAll->GetNbinsY(); thetabin++)
279 {
280 // This theta is not exactly the one of the MC events, just about
281 // the same (bins have been selected so):
282
283 Float_t theta = fHistAll->GetYaxis()->GetBinCenter(thetabin);
284 Float_t thetalo = fHistAll->GetYaxis()->GetBinLowEdge(thetabin);
285 Float_t thetahi = fHistAll->GetYaxis()->GetBinLowEdge(thetabin+1);
286
287 Float_t emin[4]; // Minimum energy in MC sample
288 Float_t emax[4]; // Maximum energy in MC sample
289 Float_t index[4]; // Spectral index
290 Float_t numevts[4]; // Number of events
291 Float_t multfactor[4]; // Factor by which the original number of events in an MC
292 // sample has been multiplied to account for the differences
293 // in the generation areas of the various samples.
294 Float_t rmax; // Maximum impact parameter range (on ground up to 45 degrees,
295 // on a plane perpendicular to Shower axis for 55 and 65 deg).
296
297 memset(emin, 0, 4*sizeof(Float_t));
298 memset(emax, 0, 4*sizeof(Float_t));
299 memset(index, 0, 4*sizeof(Float_t));
300 memset(numevts, 0, 4*sizeof(Float_t));
301 rmax = 0.;
302
303 multfactor[0] = 1.;
304 multfactor[1] = 1.;
305 multfactor[2] = 1.;
306 multfactor[3] = 1.;
307
308 //
309 // rmin and rmax are the minimum and maximum values of the impact
310 // parameter of the shower on the ground (horizontal plane).
311 //
312
313 Int_t num_MC_samples = 0;
314
315 //if (theta > 8 && theta < 9) // 8.75 deg
316 if ( thetalo<8.75 && 8.75<thetahi) // 8.75 deg
317 {
318 emin[0] = 300.;
319 emax[0] = 400.; // Energies in GeV.
320 index[0] = 1.5;
321 numevts[0] = 4000.;
322
323 emin[1] = 400.;
324 emax[1] = 30000.;
325 index[1] = 1.5;
326 numevts[1] = 25740.;
327
328 rmax = 250.; //meters
329 num_MC_samples = 2;
330 }
331 //else if (theta > 20 && theta < 21) // 20.5 deg
332 else if ( thetalo<20.5 && 20.5<thetahi) // 20.5 deg
333 {
334 emin[0] = 300.;
335 emax[0] = 400.; // Energies in GeV.
336 index[0] = 1.5;
337 numevts[0] = 6611.;
338
339 emin[1] = 400.;
340 emax[1] = 30000.;
341 index[1] = 1.5;
342 numevts[1] = 24448.;
343
344 rmax = 263.;
345 num_MC_samples = 2;
346 }
347 //else if (theta > 26 && theta < 27) // 26.5 degrees
348 else if ( thetalo<26.5 && 26.5<thetahi) // 26.5 deg
349 {
350 emin[0] = 300.;
351 emax[0] = 400.; // Energies in GeV.
352 index[0] = 1.5;
353 numevts[0] = 4000.;
354
355 emin[1] = 400.;
356 emax[1] = 30000.;
357 index[1] = 1.5;
358 numevts[1] = 26316.;
359
360 rmax = 290.; //meters
361 num_MC_samples = 2;
362 }
363 //else if (theta > 32 && theta < 33) // 32.5 degrees
364 else if ( thetalo<32.5 && 32.5<thetahi) // 32.5 deg
365 {
366 emin[0] = 300.;
367 emax[0] = 30000.; // Energies in GeV.
368 index[0] = 1.5;
369 numevts[0] = 33646.;
370
371 rmax = 350.; //meters
372 num_MC_samples = 1;
373 }
374 //else if (theta > 38 && theta < 39) // 38.75 degrees
375 else if ( thetalo<38.75 && 38.75<thetahi) // 38.75 deg
376 {
377 emin[0] = 300.;
378 emax[0] = 30000.; // Energies in GeV.
379 index[0] = 1.5;
380 numevts[0] = 38415.;
381
382 rmax = 380.; //meters
383 num_MC_samples = 1;
384 }
385 //else if (theta > 45 && theta < 47) // 46 degrees
386 else if ( thetalo<46 && 46<thetahi) // 46 deg
387 {
388 emin[0] = 300.;
389 emax[0] = 50000.; // Energies in GeV.
390 index[0] = 1.5;
391 numevts[0] = 30197.;
392
393 rmax = 565.; //meters
394 num_MC_samples = 1;
395 }
396 //else if (theta > 54 && theta < 56) // 55 degrees
397 else if ( thetalo<55 && 55<thetahi) // 55 deg
398 {
399 //
400 // The value of numevts in the first sample (below) has been
401 // changed to simplify calculations. We have multiplied it
402 // times 1.2808997 to convert it to the number it would be if
403 // the generation area was equal to that of the other samples
404 // at 55 degrees (pi*600**2 m2). This has to be taken into account
405 // in the error in the number of events.
406 //
407
408 emin[0] = 500.;
409 emax[0] = 50000.; // Energies in GeV.
410 index[0] = 1.5;
411 numevts[0] = 3298.;
412 multfactor[0] = 1.2808997;
413
414 emin[1] = 1500.;
415 emax[1] = 50000.; // Energies in GeV.
416 index[1] = 1.5;
417 numevts[1] = 22229.;
418
419 emin[2] = 1500.;
420 emax[2] = 50000.; // Energies in GeV.
421 index[2] = 1.7;
422 numevts[2] = 7553.;
423
424 rmax = 600; //meters
425 num_MC_samples = 3;
426 }
427
428 //else if (theta > 64 && theta < 66) // 65 degrees
429 else if ( thetalo<65 && 65<thetahi) // 65 deg
430 {
431 emin[0] = 2000.;
432 emax[0] = 50000.; // Energies in GeV.
433 index[0] = 1.5;
434 numevts[0] = 16310.;
435
436 emin[1] = 2000.;
437 emax[1] = 50000.; // Energies in GeV.
438 index[1] = 1.7;
439 numevts[1] = 3000.;
440
441 //
442 // The value of numevts in the next two samples (below) has been
443 // changed to simplify calculations. We have converted them to the
444 // number it would be if the generation area was equal to that of
445 // the first two samples at 65 degrees (pi*800**2 m2) (four times
446 // as many, since the original maximum impact parameter was 400
447 // instead of 800. This is taken into account in the error too.
448 //
449
450 emin[2] = 5000.;
451 emax[2] = 50000.; // Energies in GeV.
452 index[2] = 1.5;
453 numevts[2] = 56584.;
454 multfactor[2] = 4;
455
456 emin[3] = 5000.;
457 emax[3] = 50000.; // Energies in GeV.
458 index[3] = 1.7;
459 numevts[3] = 11464;
460 multfactor[3] = 4;
461
462 rmax = 800; // meters
463 num_MC_samples = 4;
464 }
465
466
467 for (Int_t i=1; i <= fHistAll->GetNbinsX(); i++)
468 {
469 Float_t e1;
470 Float_t e2;
471
472 if (fEaxis == kLog10)
473 {
474 e1 = pow(10.,fHistAll->GetXaxis()->GetBinLowEdge(i));
475 e2 = pow(10.,fHistAll->GetXaxis()->GetBinLowEdge(i+1));
476 }
477 else
478 {
479 e1 = fHistAll->GetXaxis()->GetBinLowEdge(i);
480 e2 = fHistAll->GetXaxis()->GetBinLowEdge(i+1);
481 }
482
483 Float_t events = 0.;
484 Float_t errevents = 0.;
485
486 for (Int_t sample = 0; sample < num_MC_samples; sample++)
487 {
488 Float_t expo = 1.-index[sample];
489 Float_t k = numevts[sample] / (pow(emax[sample],expo) - pow(emin[sample],expo));
490
491 if (e2 < emin[sample] || e1 > emax[sample])
492 continue;
493
494 if (emin[sample] > e1)
495 e1 = emin[sample];
496
497 if (emax[sample] < e2)
498 e2 = emax[sample];
499
500 events += k * (pow(e2, expo) - pow(e1, expo));
501 errevents += multfactor[sample] * events;
502 }
503
504 errevents= sqrt(errevents);
505
506 fHistAll->SetBinContent(i, thetabin, events);
507 fHistAll->SetBinError(i, thetabin, errevents);
508 }
509
510 // -----------------------------------------------------------
511
512 const Float_t dr = TMath::Pi() * rmax * rmax;
513
514 for (Int_t ix = 1; ix <= fHistAll->GetNbinsX(); ix++)
515 {
516 const Float_t Na = fHistAll->GetBinContent(ix,thetabin);
517
518 if (Na <= 0)
519 {
520 //
521 // If energy is large, this case means that no or very few events
522 // were generated at this energy bin. In this case we assign it
523 // the effective area of the bin below it in energy. If energy is
524 // below 1E4, it means that no events triggered -> eff area = 0
525 //
526 // NOW DISABLED: because collection area after analysis does not
527 // saturate at high E!
528 //
529
530 /*
531 if (fHistSel->GetXaxis()->GetBinLowEdge(ix) > 4.)
532 {
533 fHistCol->SetBinContent(ix, thetabin, fHistCol->GetBinContent(ix-1, thetabin));
534 fHistCol->SetBinError(ix, thetabin, fHistCol->GetBinError(ix-1, thetabin));
535 }
536 */
537 continue;
538 }
539
540 const Float_t Ns = fHistSel->GetBinContent(ix,thetabin);
541
542 // Since Na is an estimate of the total number of showers generated
543 // in the energy bin, it may happen that Ns (triggered showers) is
544 // larger than Na. In that case, the bin is skipped:
545
546 if (Na < Ns)
547 continue;
548
549 const Double_t eff = Ns/Na;
550 const Double_t efferr = sqrt((1.-eff)*Ns)/Na;
551
552 //
553 // Now we get the total area, perpendicular to the observation direction
554 // in which the events were generated (correct for cos theta):
555 //
556
557 Float_t area = dr;
558
559 if (theta < 50)
560 area *= cos(theta*TMath::Pi()/180.);
561
562 // Above 50 degrees MC was generated with Corsika 6.xx, and the cores
563 // were distributed on a circle perpendicular to the observation direction,
564 // and not on ground, hence the correction for cos(theta) is not necessary.
565 //
566
567
568 fHistCol->SetBinContent(ix, thetabin, eff*area);
569 fHistCol->SetBinError(ix, thetabin, efferr*area);
570
571 }
572 }
573
574 SetReadyToSave();
575}
Note: See TracBrowser for help on using the repository browser.