source: trunk/MagicSoft/Mars/macros/starmc2.C@ 6373

/* ======================================================================== *\
   !
   ! *
   ! * This file is part of MARS, the MAGIC Analysis and Reconstruction
   ! * Software. It is distributed to you in the hope that it can be a useful
   ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
   ! * It is distributed WITHOUT ANY WARRANTY.
   ! *
   ! * Permission to use, copy, modify and distribute this software and its
   ! * documentation for any purpose is hereby granted without fee,
   ! * provided that the above copyright notice appear in all copies and
   ! * that both that copyright notice and this permission notice appear
   ! * in supporting documentation. It is provided "as is" without express
   ! * or implied warranty.
   ! *
   !
   !
   !   Author(s): Abelardo Moralejo 1/2004 <mailto:moralejo@pd.infn.it>
   !              Thomas Bretz  5/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
   !
   !   Copyright: MAGIC Software Development, 2000-2004
   !
   !
   \* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  STARMC2 - STandard Analysis and Reconstruction (MC example)
//
//  This macro converts into image  parameters an input file of MC data 
//  previously calibrated (see mccalibrate.C).
//  
//
/////////////////////////////////////////////////////////////////////////////

#include "MImgCleanStd.h"

void starmc2()
{
  Char_t* AnalysisFilename = "calibrated_gamma.root"; // File to be analyzed

  TString* OutFilename1;
  TString* OutFilename2;

  // Change output file names as desired. If you want only one output, comment
  // out the initialization of OutFilename2:
  OutFilename1 = new TString("star_gamma_train.root");   // Output file name 1 (test)
  OutFilename2 = new TString("star_gamma_test.root");    // Output file name 2 (train)

  MImgCleanStd      clean(4.5, 3.); // Applies tail cuts to image.
  // WARNING: the tightness of the tail cuts depends on the signal extraction method
  // used in mccalibrate.C!! (some methods result in positively biased signals)

  // ------------------------------------------------------------------

  //
  // Create a empty Parameter List and an empty Task List
  // The tasklist is identified in the eventloop by its name
  //
  MParList  plist;

  MTaskList tlist;

  plist.AddToList(&tlist);

  MSrcPosCam src;
  //
  // FOR WOBBLE MODE!! Set source position on camera here.
  //  src.SetX(120.);  // units: mm

  src.SetReadyToSave();

  plist.AddToList(&src);

  //
  // Now setup the tasks and tasklist:
  // ---------------------------------
  //
  MReadMarsFile read("Events");

  read.AddFile(AnalysisFilename);

  read.DisableAutoScheme();


  MHillasCalc       hcalc; // Calculates Hillas parameters not dependent on source position.
  hcalc.Enable(MHillasCalc::kCalcHillasSrc);

  tlist.AddToList(&read);
  tlist.AddToList(&clean);
  tlist.AddToList(&hcalc);

  //
  // Open output file(s):
  //
  MWriteRootFile write1(OutFilename1->Data()); // Writes output   
  //
  // Add MC containers only if they exist. In this way you can also run on real calibrated data.
  //
  write1.AddContainer("MRawRunHeader", "RunHeaders");
  write1.AddContainer("MMcRunHeader",  "RunHeaders", kFALSE);
  write1.AddContainer("MGeomCam",      "RunHeaders", kFALSE);
  write1.AddContainer("MMcConfigRunHeader",  "RunHeaders", kFALSE);
  write1.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
  write1.AddContainer("MMcFadcHeader",  "RunHeaders", kFALSE);
  write1.AddContainer("MMcTrigHeader",  "RunHeaders", kFALSE);

  write1.AddContainer("MMcEvt",        "Events", kFALSE);
  write1.AddContainer("MPointingPos",  "Events", kFALSE);
  write1.AddContainer("MMcTrig",       "Events", kFALSE);
  write1.AddContainer("MSrcPosCam",    "Events", kFALSE);
  write1.AddContainer("MRawEvtHeader", "Events");
  write1.AddContainer("MHillas",       "Events");
  write1.AddContainer("MHillasExt",    "Events");
  write1.AddContainer("MHillasSrc",    "Events");
  write1.AddContainer("MImagePar",     "Events");
  write1.AddContainer("MNewImagePar",  "Events");
  write1.AddContainer("MConcentration","Events");

  if (OutFilename2) // Second output file, in case we want a split output
    {
      MWriteRootFile write2(OutFilename2->Data()); // Writes output
      write2.AddContainer("MRawRunHeader", "RunHeaders");
      write2.AddContainer("MMcRunHeader",  "RunHeaders", kFALSE);
      write2.AddContainer("MGeomCam",      "RunHeaders", kFALSE);
      write2.AddContainer("MMcConfigRunHeader",  "RunHeaders", kFALSE);
      write2.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
      write2.AddContainer("MMcFadcHeader",  "RunHeaders", kFALSE);
      write2.AddContainer("MMcTrigHeader",  "RunHeaders", kFALSE);

      write2.AddContainer("MMcEvt",        "Events", kFALSE);
      write2.AddContainer("MPointingPos",  "Events", kFALSE);
      write2.AddContainer("MMcTrig",       "Events", kFALSE);
      write2.AddContainer("MSrcPosCam",    "Events", kFALSE);
      write2.AddContainer("MRawEvtHeader", "Events");
      write2.AddContainer("MHillas",       "Events");
      write2.AddContainer("MHillasExt",    "Events");
      write2.AddContainer("MHillasSrc",    "Events");
      write2.AddContainer("MImagePar",     "Events");
      write2.AddContainer("MNewImagePar",  "Events");
      write2.AddContainer("MConcentration","Events");

      //
      // Divide output in train and test samples, using the event number
      // (odd/even) to achieve otherwise unbiased event samples:
      //
      
      MF filter1("{MMcEvt.fEvtNumber%2}>0.5");
      MF filter2("{MMcEvt.fEvtNumber%2}<0.5");

      write1.SetFilter(&filter1);
      write2.SetFilter(&filter2);

      tlist.AddToList(&filter1);
      tlist.AddToList(&filter2);
    }


  tlist.AddToList(&write1);

  if (OutFilename2)
    tlist.AddToList(&write2);   
  

  //
  // analysis loop
  //

  MEvtLoop evtloop;
  MProgressBar bar;
  bar.SetWindowName("Analyzing...");
  evtloop.SetProgressBar(&bar);
  evtloop.SetParList(&plist);

  if (!evtloop.Eventloop())
    return;

  tlist.PrintStatistics();

  return;
}