Changeset 6493

Timestamp:

02/15/05 17:33:43 (20 years ago)

Author:

moralejo

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/macros/collarea.C (modified) (1 diff)

trunk/MagicSoft/Mars/macros/collarea.C

@@ -17 +17 @@
 !
 !   Author(s): Thomas Bretz, 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
+!   Author(s): Abelardo Moralejo, 2/2005 <mailto:moralejo@pd.infn.it>
 !
-!   Copyright: MAGIC Software Development, 2000-2003
+!   Copyright: MAGIC Software Development, 2000-2005
 !
 !
 \* ======================================================================== */
 
+// 
+// Example macro on how to calculate effective areas. The input is a file 
+// containing Hillas parameters of the MC events surviving a certain kind 
+// of analysis. It must also contain a tree called "OriginalMC" with a branch
+// "MMcEvtBasic" and one entry per original simulated shower used to produce
+// the events in that input file. The current (2/2005) camera simulation in 
+// the cvs already writes out all events with the MMcEvtBasic container to
+// the Events tree. In Mars/macros/mccalibrate.C and starmc2.C you can see
+// how this branch is put to a separate tree "OriginalMC" and written out
+// to the "star" file. This will not work with star files not containing this
+// extra tree. Please contact moralejo@pd.infn.it in case of questions
+// about this new approach to the calculation of effective areas.
+//
 
-void collarea(TString filename="camera.root", TString outname="")
+void collarea(TString filename="star_gamma_test.root", TString outname="area.root")
 {
-    MStatusDisplay *d = new MStatusDisplay;
+  MStatusDisplay *d = new MStatusDisplay;
     // redirect logging output to GUI, kFALSE to disable stream to stdout
-    d->SetLogStream(&gLog, kTRUE);
+  d->SetLogStream(&gLog, kTRUE);
 
-    //
-    // first we have to create our empty lists
-    //
-    MParList  parlist;
-    MTaskList tasklist;
+  //
+  // First loop: fill the MHMcCollectionArea::fHistAll histogram containing
+  // all the events produced at the Corsika level:
+  //
+  MParList  parlist;
+  MTaskList tasklist;
 
-    parlist.AddToList(&tasklist);
+  parlist.AddToList(&tasklist);
 
-    //
-    // Setup out tasks:
-    //  - First we have to read the events
-    //  - Then we can fill the efficiency histograms
-    //
-    MReadMarsFile reader("Events", filename);
-    MMcCollectionAreaCalc effi;
+  //
+  // Setup out tasks:
+  //  - First we have to read the events
+  //
+  MReadMarsFile reader("OriginalMC", filename);
+  reader.DisableAutoScheme();
 
-    tasklist.AddToList(&reader);
-    tasklist.AddToList(&effi);
+  //
+  // Create collection area object and add it to the task list
+  //
+  MHMcCollectionArea collarea;
+  parlist.AddToList(&collarea);
 
-    //
-    // set up the loop for the processing
-    //
-    MEvtLoop magic;
-    magic.SetParList(&parlist);
+  //
+  // Set the COARSE binnings in which you will divide your data, to get
+  // the collection area in these same binnings. You can set also this
+  // binnings on an existing and filled MHMcCollectionAreaObject, and
+  // call Calc again to get the effective area in other coarse bins.
+  //
+  MBinning binsTheta("binsTheta");
+  MBinning binsEnergy("binsEnergy"); // name must be standardized!
+  Int_t nbins = 8;
+  TArrayD edges(nbins+1);
+  edges[0] = 0;
+  edges[1] = 10.;
+  edges[2] = 20.;
+  edges[3] = 30.;
+  edges[4] = 40.;
+  edges[5] = 50.;
+  edges[6] = 60.;
+  edges[7] = 70.;
+  edges[8] = 80.;
+  binsTheta.SetEdges(edges);              // Theta [deg]
+  binsEnergy.SetEdgesLog(20, 2., 20000);  // Energy [GeV]
+  parlist.AddToList(&binsTheta);
+  parlist.AddToList(&binsEnergy);
 
-    //
-    // Start to loop over all events
-    //
-    magic.SetDisplay(d);
-    if (!magic.Eventloop())
-        return;
+  // Task which fills the necessary histograms in MHMcCollectionArea
+  MMcCollectionAreaCalc effi;
 
-    tasklist.PrintStatistics();
+  // Tentative energy spectrum. This may slightly modify the calculation
+  // of effective areas. If not added to the parameter list, a flat spectrum
+  // in dN/dE will be assumed. "x" stands for the energy. The name of the 
+  // function must be standardized!
 
-    //
-    // Now the histogram we wanted to get out of the data is
-    // filled and can be displayed
-    //
-    if ((d = magic.GetDisplay()))
-        d->AddTab("Collection Area");
-    else
-        new TCanvas("CollArea", "Result of the collection area calculation");
+  TF1 Spectrum("Spectrum", "pow(x,-5.)");
+  effi.SetSpectrum(&Spectrum);
 
-    parlist.FindObject("MHMcCollectionArea")->DrawClone("nonew");
+  tasklist.AddToList(&reader);
+  tasklist.AddToList(&effi);
+
+  //
+  // set up the loop for the processing
+  //
+  MEvtLoop magic;
+  magic.SetParList(&parlist);
+
+  //
+  // Start to loop over all events
+  //
+  magic.SetDisplay(d);
+
+  if (!magic.Eventloop())
+    return;
+
+  tasklist.PrintStatistics();
 
 
-    //
-    // Write histogram to a file in case an output
-    // filename has been supplied
-    //
-    if (outname.IsNull())
-        return;
+  //
+  // Second loop: now fill the histogram MHMcCollectionArea::fHistSel
+  //
 
-    TFile f(outname, "recreate");
-    if (!f)
-        return;
+  MParList parlist2;
+  MTaskList tasklist2;
+  parlist2.AddToList(&tasklist2);
+  parlist2.AddToList((MHMcCollectionArea*)parlist->FindObject("MHMcCollectionArea"));
 
-    parlist.FindObject("MHMcCollectionArea")->Write();
+  MReadMarsFile reader2("Events", filename);
+  reader2.DisableAutoScheme();
+
+  tasklist2.AddToList(&reader2);
+
+  // The same task as before, now will fill MHMcCollectionArea::fHistSel
+  // and call MHMcCollectionArea::Calc in the PostProcess
+  tasklist2.AddToList(&effi);
+
+  //
+  // set up the loop for the processing
+  //
+  MEvtLoop magic2;
+  magic2.SetParList(&parlist2);
+
+  //
+  // Start to loop over all events
+  //
+  magic2.SetDisplay(d);
+
+  if (!magic2.Eventloop())
+    return;
+
+  tasklist2.PrintStatistics();
+
+  //
+  // Now the histogram we wanted to get out of the data is
+  // filled and can be displayed
+  //
+  if ((d = magic2.GetDisplay()))
+    d->AddTab("Collection Area");
+  else
+    new TCanvas("CollArea", "Result of the collection area calculation");
+
+  gStyle->SetPalette(1,0);
+  gPad->SetLogx();
+  TH2D* areaplot = ((MHMcCollectionArea*)parlist2.FindObject
+                    ("MHMcCollectionArea"))->GetHistCoarse();
+
+  areaplot->SetStats(0);
+  areaplot->SetTitle("Effective area (m^{2})");
+  areaplot->GetXaxis()->SetTitle("Energy (GeV)");
+  areaplot->GetYaxis()->SetTitle("\\theta (deg)");
+  areaplot->DrawCopy("zcol");
+
+  //
+  // Write histogram to a file in case an output
+  // filename has been supplied
+  //
+
+  if (outname.IsNull())
+    return;
+
+  TFile f(outname, "recreate");
+  if (!f)
+    return;
+
+  // Write to the output file:
+  parlist2.FindObject("MHMcCollectionArea")->Write();
 }