Changeset 6491 for trunk/MagicSoft/Mars


Ignore:
Timestamp:
02/15/05 17:20:40 (20 years ago)
Author:
moralejo
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r6489 r6491  
    2121
    2222                                                 -*-*- END OF LINE -*-*-
     23
     24 2005/02/15 Abelardo Moralejo
     25
     26   * mmontecarlo/MMcCollectionAreaCalc.[h,cc]
     27   * mhistmc/MHMcCollectionArea.[h,cc]
     28     - Changed the way of calculating final effective area for data
     29       analysis. The new approach requires the use of MC files produced
     30       with the current CVS version of camera. We now make use of the
     31       true total number of produced MC events, and allow for the
     32       setting of a "tentative" differential gamma spectrum to be used
     33       in the calculation of effective areas.
     34
     35   * macros/collarea.C
     36     - Adapted to the new way of calculating effective areas.
     37
    2338
    2439 2005/02/15 Thomas Bretz

  • trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.cc

    r3848 r6491  
    11/* ======================================================================== *\
    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): Thomas Bretz  12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
    19 !   Author(s): Harald Kornmayer 1/2001
    20 !
    21 !   Copyright: MAGIC Software Development, 2000-2002
    22 !
    23 !
    24 \* ======================================================================== */
     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): Thomas Bretz  12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
     19   !   Author(s): Harald Kornmayer 1/2001
     20   !   Author(s): Abelardo Moralejo 2/2005 <mailto:moralejo@pd.infn.it>
     21   !
     22   !   Copyright: MAGIC Software Development, 2000-2005
     23   !
     24   !
     25   \* ======================================================================== */
    2526
    2627//////////////////////////////////////////////////////////////////////////////
    … …  
    3334
    3435#include <TH2.h>
     36#include <TH3.h>
    3537#include <TCanvas.h>
     38#include <THStack.h>
     39#include <TLegend.h>
     40#include <TArrayD.h>
    3641
    3742#include "MH.h"
    3843#include "MBinning.h"
    39 #include "MHMcEfficiency.h"
    40 #include "MHMcEnergyImpact.h"
    4144
    4245#include "MLog.h"
    4346#include "MLogManip.h"
    4447
     48
    4549ClassImp(MHMcCollectionArea);
    4650
    4751using namespace std;
    4852
    49 // --------------------------------------------------------------------------
    50 //
    51 //  Creates the three necessary histograms:
     53////////////////////////////////////////////////////////////////////////////////
     54//
     55//  Constructor. Creates the three necessary histograms:
    5256//   - selected showers (input)
    5357//   - all showers (input)
    5458//   - collection area (result)
    5559//
    56 MHMcCollectionArea::MHMcCollectionArea(const char *name, const char *title)
    57   : fMCAreaRadius(300.)
     60MHMcCollectionArea::MHMcCollectionArea(const char *name, const char *title):
     61  fImpactBins(50), fImpactMax(500.), fMinEvents(10)
    5862{
    59     //   initialize the histogram for the distribution r vs E
    60     //
    61     //   we set the energy range from 2 Gev to 20000 GeV (in log 4 orders
    62     //   of magnitude) and for each order we take 25 subdivision --> 100 xbins
    63     //
    64     //   we set the radius range from 0 m to 500 m with 10 m bin --> 50 ybins
    65     //
    66     fName  = name  ? name  : "MHMcCollectionArea";
    67     fTitle = title ? title : "Collection Area vs. Energy";
    68 
    69     MBinning binsx;
    70     MBinning binsy;
    71     binsx.SetEdgesLog(100, 2., 20000);
    72     binsy.SetEdges   (50,  0,   500);
    73 
    74     fHistAll = new TH2D;
    75     fHistSel = new TH2D;
    76     fHistCol = new TH1D;
     63  fName  = name  ? name  : "CollectionArea";
     64  fTitle = title ? title : "Collection Area vs. Theta vs. Energy";
     65
     66  //
     67  //   Initialize the histogram for the distribution
     68  //   Theta vs impact parameter vs E (z, y, x)
     69  //
     70  //   As default we set the energy range from 2 Gev to 20000 GeV (in log 4
     71  //   orders of magnitude) and for each order we take 25 subdivisions -->
     72  //   100 xbins. We set the radius range from 0 m to 500 m with 10 m bin -->
     73  //   50 ybins. We make bins equally spaced in cos(theta)
     74  //
     75  //   The coarse binning (of fHistColCoarse) is not set by default, the
     76  //   PreProcess of mmc/MMcCollectionAreaCalc will do it with the binnings
     77  //   found in the parameter list.
     78  //
     79
     80  MBinning binsx;
     81  MBinning binsy;
     82  MBinning binsz;
     83
     84  Int_t nbins = 32;
     85  TArrayD edges(nbins+1);
     86   
     87  edges[0] = 0;
     88   
     89  for(int i = 0; i < nbins; i++)
     90    {
     91      Double_t x = 1 - i*0.01;  // x = cos(theta)
     92      edges[i+1] = acos(x-0.005)*kRad2Deg;
     93    }
     94
     95  binsx.SetEdgesLog(100, 2., 20000);               // Energy [GeV]
     96  binsy.SetEdges   (fImpactBins, 0, fImpactMax);   // Impact parameter [m]
     97  binsz.SetEdges   (edges);                        // Theta [deg]
     98
     99  fHistAll = new TH3D();
     100  fHistSel = new TH3D();
     101  fHistCol = new TH2D();
     102  fHistColCoarse = new TH2D();
    77103 
    78     MH::SetBinning(fHistAll, &binsx, &binsy);
    79     MH::SetBinning(fHistSel, &binsx, &binsy);
    80 
    81     fHistCol->SetName(fName);
    82     fHistAll->SetName("AllEvents");
    83     fHistSel->SetName("SelectedEvents");
    84 
    85     fHistCol->SetTitle(fTitle);
    86     fHistAll->SetTitle("All showers - Radius vs Energy distribution");
    87     fHistSel->SetTitle("Selected showers - Radius vs Energy distribution");
    88 
    89     fHistAll->SetDirectory(NULL);
    90     fHistSel->SetDirectory(NULL);
    91     fHistCol->SetDirectory(NULL);
    92 
    93     fHistAll->UseCurrentStyle();
    94     fHistSel->UseCurrentStyle();
    95     fHistCol->UseCurrentStyle();
    96 
    97     fHistAll->SetXTitle("E [GeV]");
    98     fHistAll->SetYTitle("r [m]");
    99     fHistAll->SetZTitle("N");
    100 
    101     fHistSel->SetXTitle("E [GeV]");
    102     fHistSel->SetYTitle("r [m]");
    103     fHistSel->SetYTitle("N");
    104 
    105     fHistCol->SetXTitle("E [GeV]");
    106     fHistCol->SetYTitle("A [m^{2}]");
     104  MH::SetBinning(fHistAll, &binsx, &binsy, &binsz);
     105  MH::SetBinning(fHistSel, &binsx, &binsy, &binsz);
     106
     107  fHistColCoarse->SetName(fName);
     108  fHistCol->SetName("CollAreaFineBins");
     109  fHistAll->SetName("AllEvents");
     110  fHistSel->SetName("SelectedEvents");
     111
     112  fHistAll->Sumw2();
     113  fHistSel->Sumw2();
     114
     115  fHistColCoarse->SetTitle(fTitle);
     116  fHistCol->SetTitle(fTitle);
     117  fHistAll->SetTitle("All showers - Theta vs Radius vs Energy");
     118  fHistSel->SetTitle("Selected showers - Theta vs Radius vs Energy");
     119
     120  fHistAll->SetDirectory(NULL);
     121  fHistSel->SetDirectory(NULL);
     122  fHistCol->SetDirectory(NULL);
     123  fHistColCoarse->SetDirectory(NULL);
     124
     125  fHistAll->UseCurrentStyle();
     126  fHistSel->UseCurrentStyle();
     127  fHistCol->UseCurrentStyle();
     128  fHistColCoarse->UseCurrentStyle();
     129
     130  fHistAll->SetXTitle("E [GeV]");
     131  fHistAll->SetYTitle("r [m]");
     132  fHistAll->SetZTitle("\\theta [\\circ]");
     133
     134  fHistSel->SetXTitle("E [GeV]");
     135  fHistSel->SetYTitle("r [m]");
     136  fHistSel->SetZTitle("\\theta [\\circ]");
     137
     138  fHistCol->SetXTitle("E [GeV]");
     139  fHistCol->SetYTitle("\\theta [\\circ]");
     140  fHistCol->SetZTitle("A [m^{2}]");
     141
     142  fHistColCoarse->SetXTitle("E [GeV]");
     143  fHistColCoarse->SetYTitle("\\theta [\\circ]");
     144  fHistColCoarse->SetZTitle("A [m^{2}]");
    107145}
    108146
    … …  
    113151MHMcCollectionArea::~MHMcCollectionArea()
    114152{
    115     delete fHistAll;
    116     delete fHistSel;
    117     delete fHistCol;
     153  delete fHistAll;
     154  delete fHistSel;
     155  delete fHistCol;
     156}
     157
     158// --------------------------------------------------------------------------
     159//
     160// Set the (fine) binnings of histograms fHistAll, fHistSel used in the
     161// calculations. We do not need to change impact parameter binning.
     162//
     163void MHMcCollectionArea::SetBinnings(const MBinning &binsEnergy,
     164                                     const MBinning &binsTheta)
     165{
     166  MBinning binsImpact;
     167  binsImpact.SetEdges(fImpactBins, 0., fImpactMax);
     168
     169  MH::SetBinning(fHistAll, &binsEnergy, &binsImpact, &binsTheta);
     170  MH::SetBinning(fHistSel, &binsEnergy, &binsImpact, &binsTheta);
     171
     172  fHistAll->Sumw2();
     173  fHistSel->Sumw2();
     174}
     175
     176
     177// --------------------------------------------------------------------------
     178//
     179// Set the binnings of the histogram fHistColCoarse, the effective areas
     180// in the coarse bins used in the analysis.
     181//
     182//
     183void MHMcCollectionArea::SetCoarseBinnings(const MBinning &binsEnergy,
     184                                           const MBinning &binsTheta)
     185{
     186  MH::SetBinning(fHistColCoarse, &binsEnergy, &binsTheta);
    118187}
    119188
    … …  
    122191// Fill data into the histogram which contains all showers
    123192//
    124 void MHMcCollectionArea::FillAll(Double_t energy, Double_t radius)
    125 {
    126     fHistAll->Fill(energy, radius);
     193void MHMcCollectionArea::FillAll(Double_t energy, Double_t radius, Double_t theta)
     194{
     195  fHistAll->Fill(energy, radius, theta);
    127196}
    128197
    … …  
    131200// Fill data into the histogram which contains the selected showers
    132201//
    133 void MHMcCollectionArea::FillSel(Double_t energy, Double_t radius)
    134 {
    135     fHistSel->Fill(energy, radius);
    136 }
    137 
    138 // --------------------------------------------------------------------------
    139 //
    140 // Draw the histogram with all showers
    141 //
    142 void MHMcCollectionArea::DrawAll(Option_t* option)
    143 {
    144     if (!gPad)
    145         MH::MakeDefCanvas(fHistAll);
    146 
    147     fHistAll->Draw(option);
    148 
    149     gPad->SetLogx();
    150 
    151     gPad->Modified();
    152     gPad->Update();
    153 }
    154 
    155 // --------------------------------------------------------------------------
    156 //
    157 // Draw the histogram with the selected showers only.
    158 //
    159 void MHMcCollectionArea::DrawSel(Option_t* option)
    160 {
    161     if (!gPad)
    162         MH::MakeDefCanvas(fHistSel);
    163 
    164     fHistSel->Draw(option);
    165 
    166     gPad->SetLogx();
    167 
    168     gPad->Modified();
    169     gPad->Update();
    170 }
    171 
     202void MHMcCollectionArea::FillSel(Double_t energy, Double_t radius, Double_t theta)
     203{
     204  fHistSel->Fill(energy, radius, theta);
     205}
     206
     207// --------------------------------------------------------------------------
     208//
     209// Draw
     210//
    172211void MHMcCollectionArea::Draw(Option_t* option)
    173212{
    174     TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
    175     pad->SetBorderMode(0);
    176 
    177     fHistCol->Draw(option);
    178 
    179     pad->SetLogx();
    180 
    181     pad->Modified();
    182     pad->Update();
    183 }
    184 
    185 // --------------------------------------------------------------------------
    186 //
    187 //  Calculate the Efficiency (collection area) and set the 'ReadyToSave'
    188 //  flag
    189 //
    190 void MHMcCollectionArea::CalcEfficiency()
    191 {
    192     Calc(*fHistSel, *fHistAll);
    193 }
    194 
    195 // --------------------------------------------------------------------------
    196 //
    197 //  Calculate the Efficiency (collection area) and set the 'ReadyToSave'
    198 //  flag
    199 //
    200 void MHMcCollectionArea::Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall)
    201 {
    202     Calc((TH2D&)*mcsel.GetHist(), (TH2D&)*mcall.GetHist());
    203 }
    204 
    205 // --------------------------------------------------------------------------
    206 //
    207 //  Calculate the Efficiency (collection area) and set the 'ReadyToSave'
    208 //  flag
    209 //
    210 void MHMcCollectionArea::Calc(TH2D &hsel, TH2D &hall)
    211 {
    212     MH::SetBinning(fHistCol, hsel.GetXaxis());
    213 
    214     fHistCol->Sumw2();
    215 
    216     TH1D &psel = *hsel.ProjectionX();
    217     TH1D &pall = *hall.ProjectionX();
    218 
    219     TH1D &pally = *hall.ProjectionY();
    220 
    221     Double_t max = 0.;
    222     for (Int_t i = hall.GetNbinsY(); i > 0; i--)
    223       if (pally.GetBinContent(i) > 0)
    224         {
    225           max = pally.GetBinLowEdge(i+1);
    226           break;
    227         }
    228 
    229     fHistCol->Divide(&psel, &pall, TMath::Pi()*max*max, 1);
    230     fHistCol->SetEntries(hsel.GetEntries());
    231 
    232     delete &psel;
    233     delete &pall;
    234     delete &pally;
    235 
    236     SetReadyToSave();
    237 }
    238 
    239 // --------------------------------------------------------------------------
    240 //
    241 //  Calculate the Efficiency (collection area) and set the 'ReadyToSave'
    242 //  flag
    243 //
    244 void MHMcCollectionArea::CalcEfficiency2()
    245 {
    246     TH1D &histsel = *fHistSel->ProjectionX();
    247     TH1D &histall = *fHistAll->ProjectionX();
    248 
    249     TAxis &xaxis = *histsel.GetXaxis();
    250     MH::SetBinning(fHistCol, &xaxis);
    251     const Int_t nbinx = xaxis.GetNbins();
    252 
    253     // -----------------------------------------------------------
    254     //
    255     // Impact parameter range:  TO BE FIXED! Impact parameter shoud be
    256     // read from run header, but it is not yet in!!
    257     //
    258     const Float_t r1 = 0;
    259     const Float_t r2 = fMCAreaRadius;
    260 
    261     *fLog << warn << endl << dbginf << "WARNING! I will assume a maximum impact parameter of " << fMCAreaRadius << " meters for the MC events. Check that this is the true one!" <<endl<<endl;
    262     const Float_t total_area = TMath::Pi() * (r2*r2 - r1*r1);
    263 
    264     for (Int_t ix=1; ix<=nbinx; ix++)
     213  //
     214  // Lego plot
     215  //
     216  TCanvas *c1 = new TCanvas(); 
     217  c1->SetLogx();
     218  c1->SetLogz();
     219  c1->SetGridx();
     220  c1->SetGridy();
     221
     222  fHistCol->Draw("lego2");
     223   
     224  //
     225  // Averagye Aeff
     226  //
     227  TCanvas *c2 = new TCanvas();
     228  c2->SetLogx();
     229  c2->SetLogy();
     230  c2->SetGridx();
     231  c2->SetGridy();
     232
     233  TH1D* harea = fHistCol->ProjectionX();
     234  harea->Draw("e1");
     235 
     236  //
     237  // Plot the Aeff for the different theta
     238  //
     239  TCanvas *c3 = new TCanvas();
     240  c3->SetLogx();
     241  c3->SetLogy();
     242  c3->SetGridx();
     243  c3->SetGridy();
     244
     245  TLegend * leg = new TLegend(0.73,0.65,0.89,0.89);
     246   
     247  TAxis* yaxis = fHistCol->GetYaxis();
     248  const Int_t nbiny = fHistCol->GetYaxis()->GetNbins();
     249   
     250  THStack* hs = new THStack("aa","aa");
     251
     252  hs->Add(harea,"e1");
     253  leg->AddEntry(harea,"All","l");
     254
     255  for(Int_t iy=1; iy<=nbiny; iy++)
    265256    {
    266         const Float_t Na = histall.GetBinContent(ix);
    267 
    268         if (Na <= 0)
    269             continue;
    270 
    271         const Float_t Ns = histsel.GetBinContent(ix);
    272 
    273         // Since Na is an estimate of the total number of showers generated
    274         // in the energy bin, it may happen that Ns (triggered showers) is
    275         // larger than Na. In that case, the bin is skipped:
    276 
    277         if (Na < Ns)
    278             continue;
    279 
    280         const Double_t eff = Ns/Na;
    281 
    282         const Double_t efferr = sqrt((1.-eff)*Ns)/Na;
    283        
    284         const Float_t col_area  =  eff * total_area;
    285         const Float_t col_area_error =  efferr * total_area;
    286 
    287         fHistCol->SetBinContent(ix, col_area);
    288         fHistCol->SetBinError(ix, col_area_error);
     257
     258      TH1D* h1=  fHistCol->ProjectionX(Form("%d",iy),iy,iy);
     259
     260      if(h1->GetEntries()==0)
     261        continue;
     262
     263      cout <<h1->GetEntries() << endl;
     264
     265      leg->AddEntry(h1,Form("\\theta = %.0f",yaxis->GetBinCenter(iy)),"l");
     266      h1->SetLineColor(iy);
     267      hs->Add(h1,"e1");
    289268    }
    290 
    291     delete &histsel;
    292 
    293     SetReadyToSave();
    294 }
    295 
    296 // --------------------------------------------------------------------------
    297 //
    298 //  Calculate the Efficiency (collection area) and set the 'ReadyToSave'
    299 //  flag
    300 //
    301 void MHMcCollectionArea::Calc(const MHMcEfficiency &heff)
    302 {
    303     //
    304     //  now calculate the Collection area for different
    305     //  energies
    306     //
    307     TH2D &h = (TH2D&)*heff.GetHist();
    308 
    309     MH::SetBinning(fHistCol, h.GetXaxis());
    310 
    311     const Int_t nbinx = h.GetXaxis()->GetNbins();
    312     const Int_t nbiny = h.GetYaxis()->GetNbins();
    313 
    314     for (Int_t ix=1; ix<=nbinx; ix++)
    315     {
    316         Double_t errA = 0;
    317         Double_t colA = 0;
    318 
    319         for (Int_t iy=1; iy<=nbiny; iy++)
    320         {
    321             TAxis *yaxis = h.GetYaxis();
    322 
    323             const Double_t r1  = yaxis->GetBinLowEdge(iy);
    324             const Double_t r2  = yaxis->GetBinLowEdge(iy+1);
    325 
    326             const Double_t A   = TMath::Pi() * (r2*r2 - r1*r1);
    327 
    328             const Double_t eff = h.GetCellContent(ix, iy);
    329             const Double_t efferr = h.GetCellError(ix, iy);
    330 
    331             colA += eff*A;
    332             errA += A*A * efferr*efferr;
    333         }
    334 
    335         fHistCol->SetBinContent(ix, colA);
    336         fHistCol->SetBinError(ix, sqrt(errA));
    337     }
    338 
    339     SetReadyToSave();
    340 }
     269  hs->SetMinimum(1);
     270 
     271  hs->Draw("nostack");
     272  leg->Draw();
     273
     274}
     275
     276// --------------------------------------------------------------------------
     277//
     278// Calculate the collection area and set the 'ReadyToSave' flag
     279// We first calculate the area in fine energy bins, and then do a
     280// weighted mean to obtain the area in coarse bins. The weights in
     281// the coarse bins are intended to account for the effect of the
     282// energy spectrum in the effective area itself. The weights
     283// are taken from the tentative differential spectrum dN_gam/dE given
     284// through the function "spectrum". If no such function is supplied,
     285// then no weights are applied (and hence the spectrum will be as a
     286// flat spectrum in dN_gam/dE). Of course we have a "generated" MC
     287// spectrum, but within each fine bin the differences in spectrum
     288// should not change the result (if bins are fine enough). With no
     289// supplied tentative spectrum, each fine bin is weighted equally in
     290// calculating the area in the coarse bin, and so it is like having a
     291// flat spectrum.
     292//
     293// You can run this Calc procedure on an already existing
     294// MHMcCollectionArea object, as long as it is filled.
     295//
     296void MHMcCollectionArea::Calc(TF1 *spectrum)
     297{
     298  // Search last impact parameter bin containing events
     299  // FIXME: this should be done independently for each theta angle.
     300  //
     301  TH1D &himpact = *(TH1D*)fHistAll->Project3D("y");
     302
     303  Int_t impbin;
     304  for (impbin = himpact.GetNbinsX(); impbin > 0; impbin--)
     305    if (himpact.GetBinContent(impbin)>0)
     306      break;
     307
     308  Float_t max_radius = himpact.GetBinLowEdge(impbin);
     309
     310  Float_t total_area = TMath::Pi()*max_radius*max_radius;
     311
     312  for (Int_t ix = 1; ix <= fHistAll->GetNbinsX(); ix++)
     313    for (Int_t iz = 1; iz <= fHistAll->GetNbinsZ(); iz++)
     314      {
     315        fHistAll->SetBinContent(ix, impbin, iz, 0.);
     316        fHistSel->SetBinContent(ix, impbin, iz, 0.);
     317        fHistAll->SetBinError(ix, impbin, iz, 0.);
     318        fHistSel->SetBinError(ix, impbin, iz, 0.);
     319      }
     320
     321  TH2D &histsel = *(TH2D*)fHistSel->Project3D("zx,e");
     322  TH2D &histall = *(TH2D*)fHistAll->Project3D("zx,e");
     323  // "e" option means that errors are computed!
     324
     325
     326  TAxis &xaxis = *histsel.GetXaxis();
     327  TAxis &yaxis = *histsel.GetYaxis();
     328  MH::SetBinning(fHistCol, &xaxis, &yaxis);
     329
     330  cout << "Total considered MC area = pi * " << max_radius
     331       << "^2 square meters" << endl;
     332
     333  fHistCol->Sumw2();
     334  fHistCol->Divide(&histsel, &histall, total_area, 1., "b");
     335
     336  //
     337  // Now get the effective area in the selected coarse bins. Weight
     338  // the values in the small bins according the supplied tentative
     339  // spectrum, if it has been supplied as argument of Calc.
     340  //
     341
     342  for (Int_t ibin = 1; ibin <= fHistColCoarse->GetNbinsX(); ibin++)
     343    for (Int_t jbin = 1; jbin <= fHistColCoarse->GetNbinsY(); jbin++)
     344      {
     345        Float_t maxenergy = fHistColCoarse->GetXaxis()->GetBinUpEdge(ibin);
     346        Float_t minenergy = fHistColCoarse->GetXaxis()->GetBinLowEdge(ibin);
     347
     348        Float_t maxtheta = fHistColCoarse->GetYaxis()->GetBinUpEdge(jbin);
     349        Float_t mintheta = fHistColCoarse->GetYaxis()->GetBinLowEdge(jbin);
     350
     351        // Fine bins ranges covered by the coarse bin ibin, jbin:
     352        Int_t ibin2max = fHistCol->GetXaxis()->FindBin(maxenergy);
     353        Int_t ibin2min = fHistCol->GetXaxis()->FindBin(minenergy);
     354
     355        Int_t jbin2max = fHistCol->GetYaxis()->FindBin(maxtheta);
     356        Int_t jbin2min = fHistCol->GetYaxis()->FindBin(mintheta);
     357
     358        Float_t area = 0.;
     359        Float_t errarea = 0.;
     360        Float_t norm = 0;
     361
     362        for (Int_t ibin2 = ibin2min; ibin2 <= ibin2max; ibin2++)
     363          {
     364            Float_t weight = spectrum? spectrum->
     365              Eval(fHistCol->GetXaxis()->GetBinCenter(ibin2)) : 1.;
     366
     367            for (Int_t jbin2 = jbin2min; jbin2 <= jbin2max; jbin2++)
     368              {
     369                // Skip bins with too few produced MC events
     370                if (histall.GetBinContent(ibin2,jbin2) < fMinEvents)
     371                  continue;
     372
     373                area += weight * fHistCol->GetBinContent(ibin2,jbin2);
     374                norm += weight;
     375                errarea += pow(weight * fHistCol->
     376                               GetBinError(ibin2,jbin2), 2.);
     377              }
     378          }
     379        if (norm > 0.)
     380          {
     381            area /= norm;
     382            errarea = sqrt(errarea)/norm;
     383          }
     384
     385        fHistColCoarse->SetBinContent(ibin, jbin, area);
     386        fHistColCoarse->SetBinError(ibin, jbin, errarea);
     387      }
     388
     389  SetReadyToSave();
     390}
     391     
     392

  • trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.h

    r3848 r6491  
    66#endif
    77
     8#include "MBinning.h"
     9#include "TF1.h"
     10
    811class TH1D;
    912class TH2D;
     13class TH3D;
    1014
    11 class MHMcEfficiency;
    12 class MHMcEnergyImpact;
     15//
     16// Version 2:
     17//
     18// Now all histograms are saved with the object. Added members
     19// fImpactMax, fImpactBins, fMinEvents
     20//
    1321
    1422class MHMcCollectionArea : public MH
    1523{
    1624private:
    17     TH2D *fHistAll; //! all simulated showers
    18     TH2D *fHistSel; //! the selected showers
     25  TH3D *fHistAll;        // all simulated showers
     26  TH3D *fHistSel;        // the selected showers
     27  TH2D *fHistCol;        // the collection area in fine bins
     28  TH2D *fHistColCoarse;  // the coll. area in coarse bins
    1929
    20     TH1D *fHistCol; //  the collection area
     30  Int_t   fImpactBins;   // number of bins in i.p. for histograms
     31  Float_t fImpactMax;    // [m] Maximum impact parameter for histograms
    2132
    22     Float_t fMCAreaRadius; // [m] Radius of circle within which the cores
    23                            // of Corsika events have been uniformly
    24                            // distributed.
     33  Int_t   fMinEvents;   
     34  // minimum number of events in each of the "fine bins" of fHistAll
     35  // so that the bin is used in the averaging for the coarse bins.
    2536
    26     void Calc(TH2D &hsel, TH2D &hall);
     37  void Calc(TH3D &hsel, TH3D &hall);
    2738
    2839public:
    29     MHMcCollectionArea(const char *name=NULL, const char *title=NULL);
    30     ~MHMcCollectionArea();
     40  MHMcCollectionArea(const char *name=NULL, const char *title=NULL);
     41  ~MHMcCollectionArea();
    3142
    32     void FillAll(Double_t energy, Double_t radius);
    33     void FillSel(Double_t energy, Double_t radius);
     43  void FillAll(Double_t energy, Double_t radius, Double_t theta);
     44  void FillSel(Double_t energy, Double_t radius, Double_t theta);
    3445
    35     void DrawAll(Option_t *option="");
    36     void DrawSel(Option_t *option="");
     46  const TH2D *GetHist()       const { return fHistCol; }
     47  const TH2D *GetHistCoarse() const { return fHistColCoarse; }
     48  const TH3D *GetHistAll()    const { return fHistAll; }
     49  const TH3D *GetHistSel()    const { return fHistSel; }
    3750
    38     const TH1D *GetHist()       { return fHistCol; }
    39     const TH1D *GetHist() const { return fHistCol; }
     51  void SetBinnings(const MBinning &binsEnergy, const MBinning &binsTheta);
     52  void SetCoarseBinnings(const MBinning &binsEnergy, const MBinning &binsTheta);
    4053
    41     TH2D *GetHistAll()    { return fHistAll; }
     54  void Draw(Option_t *option="");
    4255
    43     void Draw(Option_t *option="");
     56  void Calc(TF1 *spectrum);
    4457
    45     void CalcEfficiency();
    46     void CalcEfficiency2();
    47 
    48     void Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall);
    49     void Calc(const MHMcEfficiency &heff);
    50 
    51     void SetMCAreaRadius(Float_t x) { fMCAreaRadius = x; }
    52 
    53     ClassDef(MHMcCollectionArea, 1)  // Data Container to calculate Collection Area
    54 };
     58  ClassDef(MHMcCollectionArea, 2)  // Data Container to store Collection Area
     59    };
    5560
    5661#endif

  • trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc

    r5340 r6491  
    1717!
    1818!   Author(s): Thomas Bretz  12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
    19 !   Author(s): Harald Kornmayer 1/2001
     19!   Author(s): Harald Kornmayer  1/2001
     20!   Author(s): Abelardo Moralejo 2/2005 <mailto:moralejo@pd.infn.it>
    2021!
    2122!   Copyright: MAGIC Software Development, 2000-2002
    … …  
    2627//////////////////////////////////////////////////////////////////////////////
    2728//
    28 //  MHMcCollectionAreaCalc
    29 //
    30 //  Remark: The initialization is mainly done in the ReInit function.
    31 //          Please make sure, that you don't use MReadTree when processing
    32 //          a file. Use a 'ReInit'-calling task like MReadMarsFile
     29//  MMcCollectionAreaCalc
    3330//
    3431//////////////////////////////////////////////////////////////////////////////
    … …  
    4239
    4340#include "MMcEvt.hxx"
    44 #include "MMcTrig.hxx"
     41#include "MMcEvtBasic.h"
     42
    4543#include "MMcRunHeader.hxx"
    4644#include "MMcCorsikaRunHeader.h"
    … …  
    5250using namespace std;
    5351
    54 MMcCollectionAreaCalc::MMcCollectionAreaCalc(const char *input,
    55                                              const char *name, const char *title)
     52////////////////////////////////////////////////////////////////////////////////
     53//
     54//  Constructor
     55//
     56MMcCollectionAreaCalc::MMcCollectionAreaCalc(const char *input, const char *name,
     57                                             const char *title):
     58  fBinsTheta(0), fBinsEnergy(0), fSpectrum(0)
    5659{
    5760    fName  = name  ? name  : "MMcCollectionAreaCalc";
    5861    fTitle = title ? title : "Task to calculate the collection area";
    59 
    60     fObjName = input ? input : "MMcTrig";
    61     AddToBranchList(Form("%s.fNumFirstLevel", input));
    62 
    63     AddToBranchList("MMcEvt.fEnergy");
    64     AddToBranchList("MMcEvt.fImpact");
    65 
    66     fCorsikaVersion           =  0;
    67     fAllEvtsTriggered         =  kFALSE;
    68     fTotalNumSimulatedShowers =  0;
    69     fThetaMin                 = -1.;
    70     fThetaMax                 = -1.;
    71 
    7262}
    7363
     64////////////////////////////////////////////////////////////////////////////////
     65//
     66// PreProcess.
     67// Create MHMcCollectionArea object if necessary.
     68// Search in parameter list for MBinning objects with binning in theta and E.
     69// These contain the coarse binning to be used in the analysis. Then search
     70// for other necessary input containers:
     71// if MMcEvt is found, it means we are in the loop over the Events tree,
     72// and so we must fill the histogram MHMcCollectionArea::fHistSel (using
     73// MHMcCollectionArea::FillSel). If MMcEvt is not found, it means that we are in
     74// the loop over the "OriginalMC" tree, containing all the original Corsika events
     75// produced, and hence we must fill the histogram  fHistAll through
     76// MHMcCollectionArea::FillAll.
     77//
    7478Int_t MMcCollectionAreaCalc::PreProcess (MParList *pList)
    7579{
    76     // connect the raw data with this task
     80  fCollArea = (MHMcCollectionArea*)pList->FindCreateObj("MHMcCollectionArea");
    7781
    78     fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
    79     if (!fMcEvt)
     82  // Look for the binnings of the histograms if they have not been already
     83  // found in a previous loop.
     84
     85  if (!fBinsTheta)
    8086    {
    81         *fLog << err << "MMcEvt not found... abort." << endl;
    82         return kFALSE;
     87      fBinsTheta = (MBinning*) pList->FindObject("binsTheta");
     88      if (!fBinsTheta)
     89        {
     90          *fLog << err << "Coarse Theta binning not found... Aborting."
     91            << endl;
     92          return kFALSE;
     93        }
     94    }
     95
     96  if (!fBinsEnergy)
     97    {
     98      fBinsEnergy = (MBinning*) pList->FindObject("binsEnergy");
     99      if (!fBinsEnergy)
     100        {
     101          *fLog << err << "Coarse Energy binning not found... Aborting."
     102                << endl;
     103          return kFALSE;
     104        }
     105    }
     106
     107  fCollArea->SetCoarseBinnings(*fBinsEnergy, *fBinsTheta);
     108
     109
     110  // Look for the input containers
     111
     112  fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
     113  if (fMcEvt)
     114    {
     115        *fLog << inf << "MMcEvt found... I will fill MHMcCollectionArea.fHistSel..." << endl;
     116        return kTRUE;
     117    }
     118
     119  *fLog << inf << "MMcEvt not found... looking for MMcEvtBasic..." << endl;
     120
     121
     122  fMcEvtBasic = (MMcEvtBasic*) pList->FindObject("MMcEvtBasic");
     123   
     124  if (fMcEvtBasic)
     125    {
     126      *fLog << inf << "MMcEvtBasic found... I will fill MHMcCollectionArea.fHistAll..." << endl;
     127      return kTRUE;
     128    }
     129
     130  *fLog << err << "MMcEvtBasic not found. Aborting..." << endl;
     131
     132  return kFALSE;
     133}
     134
     135////////////////////////////////////////////////////////////////////////////////
     136//
     137// Process. Depending on whether fMcEvt or fMcEvtBasic are available, fill
     138// MHMcCollectionArea::fHistAll, else, if fMcEvt is available, fill
     139// MHMcCollectionArea::fHistSel
     140//
     141Int_t MMcCollectionAreaCalc::Process()
     142{
     143    Double_t energy = fMcEvt? fMcEvt->GetEnergy() : fMcEvtBasic->GetEnergy();
     144    Double_t impact = fMcEvt? fMcEvt->GetImpact()/100. : fMcEvtBasic->GetImpact()/100.; // in m
     145    Double_t theta  = fMcEvt? fMcEvt->GetTelescopeTheta()*TMath::RadToDeg() :
     146      fMcEvtBasic->GetTelescopeTheta()*TMath::RadToDeg(); // in deg
     147
     148    if (fMcEvt)
     149      fCollArea->FillSel(energy, impact, theta);
     150    else
     151      fCollArea->FillAll(energy, impact, theta);
     152
     153    return kTRUE;
     154}
     155
     156///////////////////////////////////////////////////////////////////////////////
     157//
     158// Postprocess. If both fHistAll and fHistSel are already filled, calculate
     159// effective areas. Else it means we still have to run one more loop.
     160//
     161Int_t MMcCollectionAreaCalc::PostProcess()
     162{
     163  if ( ((TH2D*)fCollArea->GetHistAll())->GetEntries() > 0 &&
     164       ((TH2D*)fCollArea->GetHistSel())->GetEntries() > 0)
     165    {
     166      *fLog << inf << "Calculation Collection Area..." << endl;
     167      fCollArea->Calc(fSpectrum);
    83168    }
    84169
    85170    return kTRUE;
    86171}
    87 
    88 Bool_t MMcCollectionAreaCalc::ReInit(MParList *plist)
    89 {
    90     fCollArea=0;
    91 
    92     MMcRunHeader *runheader = (MMcRunHeader*)plist->FindObject("MMcRunHeader");
    93     if (!runheader)
    94     {
    95         *fLog << err << "Error - MMcRunHeader not found... abort." << endl;
    96         return kFALSE;
    97     }
    98 
    99     fTotalNumSimulatedShowers += runheader->GetNumSimulatedShowers();
    100     *fLog << inf << "Total Number of Simulated showers: " << fTotalNumSimulatedShowers << endl;
    101 
    102 
    103     if ( (fThetaMin >= 0 && fThetaMin != runheader->GetShowerThetaMin()) ||
    104          (fThetaMax >= 0 && fThetaMax != runheader->GetShowerThetaMax()) )
    105     {
    106         *fLog << warn << "Warning - Read files have different Theta ranges (";
    107         *fLog << "(" << fThetaMin << ", " << fThetaMax << ")   vs   (" <<
    108           runheader->GetShowerThetaMin() << ", " <<
    109           runheader->GetShowerThetaMax() << ")..." << endl;
    110     }
    111     fThetaMin = runheader->GetShowerThetaMin();
    112     fThetaMax = runheader->GetShowerThetaMax();
    113 
    114 
    115     if (fCorsikaVersion!=0 && fCorsikaVersion!=runheader->GetCorsikaVersion())
    116         *fLog << warn << "Warning - Read files have different Corsika versions..." << endl;
    117     fCorsikaVersion = runheader->GetCorsikaVersion();
    118 
    119 
    120     fAllEvtsTriggered |= runheader->GetAllEvtsTriggered();
    121     *fLog << inf << "Only triggered events avail: " << (fAllEvtsTriggered?"yes":"no") << endl;
    122 
    123 
    124     fCollArea = (MHMcCollectionArea*)plist->FindCreateObj("MHMcCollectionArea");
    125     if (!fCollArea)
    126         return kFALSE;
    127 
    128     if (!fAllEvtsTriggered)
    129     {
    130         fMcTrig = (MMcTrig*)plist->FindObject(fObjName);
    131         if (!fMcTrig)
    132         {
    133             *fLog << err << fObjName << " [MMcTrig] not found... anort." << endl;
    134             return kFALSE;
    135         }
    136         return kTRUE;
    137     }
    138 
    139     MMcCorsikaRunHeader *corrunheader  = (MMcCorsikaRunHeader*)plist->FindObject("MMcCorsikaRunHeader");
    140     if (!corrunheader)
    141     {
    142         *fLog << err << "MMcCorsikaRunHeader not found... abort." << endl;
    143         return kFALSE;
    144     }
    145 
    146     //
    147     // Calculate approximately the original number of events in each
    148     // energy bin:
    149     //
    150 
    151     const Float_t emin = corrunheader->GetELowLim();
    152     const Float_t emax = corrunheader->GetEUppLim();
    153     const Float_t expo = 1 + corrunheader->GetSlopeSpec();
    154     const Float_t k = runheader->GetNumSimulatedShowers() /
    155         (pow(emax,expo) - pow(emin,expo)) ;
    156 
    157     TH2 &hall = *fCollArea->GetHistAll();
    158 
    159     const Int_t nbinx = hall.GetNbinsX();
    160 
    161     TAxis &axe = *hall.GetXaxis();
    162     for (Int_t i = 1; i <= nbinx; i++)
    163     {
    164         const Float_t e1 = axe.GetBinLowEdge(i);
    165         const Float_t e2 = axe.GetBinLowEdge(i+1);
    166 
    167         if (e1 < emin || e2 > emax)
    168             continue;
    169 
    170         const Float_t events = k * (pow(e2, expo) - pow(e1, expo));
    171         //
    172         // We fill the i-th energy bin, with the total number of events
    173         // Second argument of Fill would be impact parameter of each
    174         // event, but we don't really need it for the collection area,
    175         // so we just put a dummy value (1.)
    176         //
    177 
    178         const Float_t energy = (e1+e2)/2.;
    179         hall.Fill(energy, 1., events);
    180     }
    181 
    182     return kTRUE;
    183 }
    184 
    185 Int_t MMcCollectionAreaCalc::Process()
    186 {
    187     const Double_t energy = fMcEvt->GetEnergy();
    188     const Double_t impact = fMcEvt->GetImpact()/100.;
    189 
    190     //
    191     // This happens for camera files created with Camera 0.5
    192     //
    193     if (TMath::IsNaN(impact))
    194     {
    195         *fLog << err << dbginf << "ERROR - Impact=NaN (Not a number)... abort." << endl;
    196         return kERROR;
    197     }
    198 
    199     if (!fAllEvtsTriggered)
    200     {
    201         fCollArea->FillAll(energy, impact);
    202 
    203         if (fMcTrig->GetFirstLevel() <= 0)
    204             return kTRUE;
    205     }
    206 
    207     fCollArea->FillSel(energy, impact);
    208 
    209     return kTRUE;
    210 }
    211 
    212 Int_t MMcCollectionAreaCalc::PostProcess()
    213 {
    214     if (!fCollArea)
    215         return kTRUE;
    216 
    217     //
    218     //   do the calculation of the effectiv area
    219     //
    220     *fLog << inf << "Calculation Collection Area..." << endl;
    221 
    222     if (!fAllEvtsTriggered)
    223         fCollArea->CalcEfficiency();
    224     else
    225     {
    226         *fLog << inf << "Total number of showers: " << fTotalNumSimulatedShowers << endl;
    227         fCollArea->CalcEfficiency2();
    228     }
    229 
    230     return kTRUE;
    231 }
    232 

  • trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.h

    r5340 r6491  
    77
    88#include <TH2.h>
     9#include <TF1.h>
    910
    1011class MParList;
    1112class MMcEvt;
     13class MMcEvtBasic;
    1214class MMcTrig;
    1315class MHMcCollectionArea;
     16class MBinning;
    1417
    1518class MMcCollectionAreaCalc : public MTask
    1619{
    1720private:
    18     const MMcEvt  *fMcEvt;
    19     const MMcTrig *fMcTrig;
     21    const MMcEvt       *fMcEvt;
     22    const MMcEvtBasic  *fMcEvtBasic;
     23    const MMcTrig      *fMcTrig;
     24
     25    MBinning           *fBinsTheta;
     26    MBinning           *fBinsEnergy;
     27    // Coarse zenith angle and energy bins used in the analysis
     28
     29    TF1                *fSpectrum;
     30    // Tentative energy spectrum. This modifies slightly the calculation
     31    // of the effective area (see MHMcCollectionArea::Calc)
     32
    2033
    2134    MHMcCollectionArea *fCollArea;
    … …  
    2336    TString fObjName;
    2437
    25     UInt_t fTotalNumSimulatedShowers;
    26     UInt_t fCorsikaVersion;
    27     Bool_t fAllEvtsTriggered;
    28     Float_t fThetaMin;
    29     Float_t fThetaMax;
    30 
    31     Bool_t ReInit(MParList *plist);
    32 
    33     Int_t PreProcess(MParList *pList);
    34     Int_t Process();
    35     Int_t PostProcess();
     38    Int_t  PreProcess(MParList *pList);
     39    Int_t  Process();
     40    Int_t  PostProcess();
    3641
    3742public:
    38     MMcCollectionAreaCalc(const char *input=NULL,
    39                           const char *name=NULL, const char *title=NULL);
     43    MMcCollectionAreaCalc(const char *input = NULL,
     44                          const char *name = NULL, const char *title = NULL);
     45
     46    void SetSpectrum(TF1 *f) { fSpectrum = f; }
    4047
    4148    ClassDef(MMcCollectionAreaCalc, 0) // Task to calculate the collection area histogram
Note: See TracChangeset for help on using the changeset viewer.