Ignore:
Timestamp:
02/15/05 17:33:43 (20 years ago)
Author:
moralejo
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/macros/collarea.C

    r2098 r6493  
    1717!
    1818!   Author(s): Thomas Bretz, 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
     19!   Author(s): Abelardo Moralejo, 2/2005 <mailto:moralejo@pd.infn.it>
    1920!
    20 !   Copyright: MAGIC Software Development, 2000-2003
     21!   Copyright: MAGIC Software Development, 2000-2005
    2122!
    2223!
    2324\* ======================================================================== */
    2425
     26//
     27// Example macro on how to calculate effective areas. The input is a file
     28// containing Hillas parameters of the MC events surviving a certain kind
     29// of analysis. It must also contain a tree called "OriginalMC" with a branch
     30// "MMcEvtBasic" and one entry per original simulated shower used to produce
     31// the events in that input file. The current (2/2005) camera simulation in
     32// the cvs already writes out all events with the MMcEvtBasic container to
     33// the Events tree. In Mars/macros/mccalibrate.C and starmc2.C you can see
     34// how this branch is put to a separate tree "OriginalMC" and written out
     35// to the "star" file. This will not work with star files not containing this
     36// extra tree. Please contact moralejo@pd.infn.it in case of questions
     37// about this new approach to the calculation of effective areas.
     38//
    2539
    26 void collarea(TString filename="camera.root", TString outname="")
     40void collarea(TString filename="star_gamma_test.root", TString outname="area.root")
    2741{
    28     MStatusDisplay *d = new MStatusDisplay;
     42  MStatusDisplay *d = new MStatusDisplay;
    2943    // redirect logging output to GUI, kFALSE to disable stream to stdout
    30     d->SetLogStream(&gLog, kTRUE);
     44  d->SetLogStream(&gLog, kTRUE);
    3145
    32     //
    33     // first we have to create our empty lists
    34     //
    35     MParList  parlist;
    36     MTaskList tasklist;
     46  //
     47  // First loop: fill the MHMcCollectionArea::fHistAll histogram containing
     48  // all the events produced at the Corsika level:
     49  //
     50  MParList  parlist;
     51  MTaskList tasklist;
    3752
    38     parlist.AddToList(&tasklist);
     53  parlist.AddToList(&tasklist);
    3954
    40     //
    41     // Setup out tasks:
    42     //  - First we have to read the events
    43     //  - Then we can fill the efficiency histograms
    44     //
    45     MReadMarsFile reader("Events", filename);
    46     MMcCollectionAreaCalc effi;
     55  //
     56  // Setup out tasks:
     57  //  - First we have to read the events
     58  //
     59  MReadMarsFile reader("OriginalMC", filename);
     60  reader.DisableAutoScheme();
    4761
    48     tasklist.AddToList(&reader);
    49     tasklist.AddToList(&effi);
     62  //
     63  // Create collection area object and add it to the task list
     64  //
     65  MHMcCollectionArea collarea;
     66  parlist.AddToList(&collarea);
    5067
    51     //
    52     // set up the loop for the processing
    53     //
    54     MEvtLoop magic;
    55     magic.SetParList(&parlist);
     68  //
     69  // Set the COARSE binnings in which you will divide your data, to get
     70  // the collection area in these same binnings. You can set also this
     71  // binnings on an existing and filled MHMcCollectionAreaObject, and
     72  // call Calc again to get the effective area in other coarse bins.
     73  //
     74  MBinning binsTheta("binsTheta");
     75  MBinning binsEnergy("binsEnergy"); // name must be standardized!
     76  Int_t nbins = 8;
     77  TArrayD edges(nbins+1);
     78  edges[0] = 0;
     79  edges[1] = 10.;
     80  edges[2] = 20.;
     81  edges[3] = 30.;
     82  edges[4] = 40.;
     83  edges[5] = 50.;
     84  edges[6] = 60.;
     85  edges[7] = 70.;
     86  edges[8] = 80.;
     87  binsTheta.SetEdges(edges);              // Theta [deg]
     88  binsEnergy.SetEdgesLog(20, 2., 20000);  // Energy [GeV]
     89  parlist.AddToList(&binsTheta);
     90  parlist.AddToList(&binsEnergy);
    5691
    57     //
    58     // Start to loop over all events
    59     //
    60     magic.SetDisplay(d);
    61     if (!magic.Eventloop())
    62         return;
     92  // Task which fills the necessary histograms in MHMcCollectionArea
     93  MMcCollectionAreaCalc effi;
    6394
    64     tasklist.PrintStatistics();
     95  // Tentative energy spectrum. This may slightly modify the calculation
     96  // of effective areas. If not added to the parameter list, a flat spectrum
     97  // in dN/dE will be assumed. "x" stands for the energy. The name of the
     98  // function must be standardized!
    6599
    66     //
    67     // Now the histogram we wanted to get out of the data is
    68     // filled and can be displayed
    69     //
    70     if ((d = magic.GetDisplay()))
    71         d->AddTab("Collection Area");
    72     else
    73         new TCanvas("CollArea", "Result of the collection area calculation");
     100  TF1 Spectrum("Spectrum", "pow(x,-5.)");
     101  effi.SetSpectrum(&Spectrum);
    74102
    75     parlist.FindObject("MHMcCollectionArea")->DrawClone("nonew");
     103  tasklist.AddToList(&reader);
     104  tasklist.AddToList(&effi);
     105
     106  //
     107  // set up the loop for the processing
     108  //
     109  MEvtLoop magic;
     110  magic.SetParList(&parlist);
     111
     112  //
     113  // Start to loop over all events
     114  //
     115  magic.SetDisplay(d);
     116
     117  if (!magic.Eventloop())
     118    return;
     119
     120  tasklist.PrintStatistics();
    76121
    77122
    78     //
    79     // Write histogram to a file in case an output
    80     // filename has been supplied
    81     //
    82     if (outname.IsNull())
    83         return;
     123  //
     124  // Second loop: now fill the histogram MHMcCollectionArea::fHistSel
     125  //
    84126
    85     TFile f(outname, "recreate");
    86     if (!f)
    87         return;
     127  MParList parlist2;
     128  MTaskList tasklist2;
     129  parlist2.AddToList(&tasklist2);
     130  parlist2.AddToList((MHMcCollectionArea*)parlist->FindObject("MHMcCollectionArea"));
    88131
    89     parlist.FindObject("MHMcCollectionArea")->Write();
     132  MReadMarsFile reader2("Events", filename);
     133  reader2.DisableAutoScheme();
     134
     135  tasklist2.AddToList(&reader2);
     136
     137  // The same task as before, now will fill MHMcCollectionArea::fHistSel
     138  // and call MHMcCollectionArea::Calc in the PostProcess
     139  tasklist2.AddToList(&effi);
     140
     141  //
     142  // set up the loop for the processing
     143  //
     144  MEvtLoop magic2;
     145  magic2.SetParList(&parlist2);
     146
     147  //
     148  // Start to loop over all events
     149  //
     150  magic2.SetDisplay(d);
     151
     152  if (!magic2.Eventloop())
     153    return;
     154
     155  tasklist2.PrintStatistics();
     156
     157  //
     158  // Now the histogram we wanted to get out of the data is
     159  // filled and can be displayed
     160  //
     161  if ((d = magic2.GetDisplay()))
     162    d->AddTab("Collection Area");
     163  else
     164    new TCanvas("CollArea", "Result of the collection area calculation");
     165
     166  gStyle->SetPalette(1,0);
     167  gPad->SetLogx();
     168  TH2D* areaplot = ((MHMcCollectionArea*)parlist2.FindObject
     169                    ("MHMcCollectionArea"))->GetHistCoarse();
     170
     171  areaplot->SetStats(0);
     172  areaplot->SetTitle("Effective area (m^{2})");
     173  areaplot->GetXaxis()->SetTitle("Energy (GeV)");
     174  areaplot->GetYaxis()->SetTitle("\\theta (deg)");
     175  areaplot->DrawCopy("zcol");
     176
     177  //
     178  // Write histogram to a file in case an output
     179  // filename has been supplied
     180  //
     181
     182  if (outname.IsNull())
     183    return;
     184
     185  TFile f(outname, "recreate");
     186  if (!f)
     187    return;
     188
     189  // Write to the output file:
     190  parlist2.FindObject("MHMcCollectionArea")->Write();
    90191}
Note: See TracChangeset for help on using the changeset viewer.