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

Legend:

Unmodified
Added
Removed

  • 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
Note: See TracChangeset for help on using the changeset viewer.