source: trunk/MagicSoft/Mars/macros/starmcstereo.C@ 3093

/* ======================================================================== *\
!
! *
! * 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
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  STARMCSTEREO - STandard Analysis and Reconstruction (for MC stereo files)
//
//  This macro is the standard converter to convert raw data from stereo 
//  camera simulation into image parameters
//
/////////////////////////////////////////////////////////////////////////////

void starmcstereo(Int_t ct1 = 2, Int_t ct2 = 3)
{
  Int_t CT[2] = {ct1, ct2};

  Int_t NCTs = sizeof(CT)/sizeof(CT[0]);


  // ------------- user change -----------------

  Char_t* AnalysisFilename = "G_XX_*_w0.root";  // File to be analyzed
  Char_t* OutFileTag      = "gammas";           // Output file tag

  Float_t CleanLev[2] = {4., 3.}; // Tail cuts for image analysis

  Int_t BinsHigh[2] = {0, 5}; // First and last FADC bin of the range to be integrated,
  Int_t BinsLow[2]  = {0, 5}; // for high and low gain respectively.

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

  //
  // This is a demonstration program which calculates the image 
  // parameters from a Magic Monte Carlo root file. Units of Size are
  // for the moment, FADC counts.
  //
  //
  // 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);

  TString s = " MSrcPosCam;";
  s += CT[0];
  MSrcPosCam src1(s);
  src1.SetReadyToSave();
  plist.AddToList(&src1);

  if (NCTs==2)
    {
      s = " MSrcPosCam;";
      s += CT[1];
      MSrcPosCam src2(s);
      src2.SetReadyToSave();
      plist.AddToList(&src2);
    }
  //
  // Now setup the tasks and tasklist:
  // ---------------------------------
  //
  MReadMarsFile read("Events");
  read.DisableAutoScheme();

  read.AddFile(AnalysisFilename);

  MGeomApply*        apply = new MGeomApply[NCTs];

  MMcPedestalCopy*   pcopy = new MMcPedestalCopy[NCTs];

  MExtractSignal*      sigextract = new MExtractSignal[NCTs];

  MMcCalibrationUpdate*  mccalibupdate = new MMcCalibrationUpdate[NCTs];
  MCalibrate* calib = new MCalibrate[NCTs];

  MImgCleanStd**     clean = new MImgCleanStd*[NCTs];

  MHillasCalc*       hcalc = new MHillasCalc[NCTs];
  MHillasSrcCalc*    scalc = new MHillasSrcCalc[NCTs];

  TString outfile = "star_";
  outfile += CT[0];
  if (NCTs==2)
    {
      outfile += "_";
      outfile += CT[1];
    }

  //
  // We have two output files (will be later train and test sampls for random forest)
  //
  outfile += "_";
  outfile += OutFileTag;
  outfile += "_train.root";
  MWriteRootFile    write1(outfile);

  outfile = "star_";
  outfile += CT[0];
  if (NCTs==2)
    {
      outfile += "_";
      outfile += CT[1];
    }

  outfile += "_";
  outfile += OutFileTag; 
  outfile += "_test.root";

  MWriteRootFile write2(outfile);

  for (Int_t i = 0; i < NCTs; i++)
    {
      apply[i]->SetSerialNumber(CT[i]);

      pcopy[i]->SetSerialNumber(CT[i]);

      sigextract[i]->SetSerialNumber(CT[i]);
      sigextract[i].SetRange(BinsHigh[0], BinsHigh[1], BinsLow[0], BinsLow[1]);

      mccalibupdate[i]->SetSerialNumber(CT[i]);
      calib[i]->SetSerialNumber(CT[i]);

      clean[i] = new MImgCleanStd(CleanLev[0], CleanLev[1]);
      clean[i]->SetSerialNumber(CT[i]);

      hcalc[i]->SetSerialNumber(CT[i]);
      scalc[i]->SetSerialNumber(CT[i]);

      write1.SetSerialNumber(CT[i]);
      write2.SetSerialNumber(CT[i]);

      write1.AddContainer(write1.AddSerialNumber("MMcEvt"),       "Events");
      write1.AddContainer(write1.AddSerialNumber("MHillas"),      "Events");
      write1.AddContainer(write1.AddSerialNumber("MHillasExt"),   "Events");
      write1.AddContainer(write1.AddSerialNumber("MHillasSrc"),   "Events");
      write1.AddContainer(write1.AddSerialNumber("MNewImagePar"), "Events");
      write1.AddContainer(write1.AddSerialNumber("MSrcPosCam"),   "RunHeaders");
      write2.AddContainer(write2.AddSerialNumber("MMcEvt"),       "Events");
      write2.AddContainer(write2.AddSerialNumber("MHillas"),      "Events");
      write2.AddContainer(write2.AddSerialNumber("MHillasExt"),   "Events");
      write2.AddContainer(write2.AddSerialNumber("MHillasSrc"),   "Events");
      write2.AddContainer(write2.AddSerialNumber("MNewImagePar"), "Events");
      write2.AddContainer(write2.AddSerialNumber("MSrcPosCam"),   "RunHeaders");
    }

  if (NCTs==2)
    {
      write1.AddContainer("MStereoPar", "Events");
      write2.AddContainer("MStereoPar", "Events");
    }

  write1.AddContainer("MRawRunHeader", "RunHeaders");
  write1.AddContainer("MMcRunHeader",  "RunHeaders");

  write2.AddContainer("MRawRunHeader", "RunHeaders");
  write2.AddContainer("MMcRunHeader",  "RunHeaders");

  tlist.AddToList(&read);

  for (i = 0; i < NCTs; i++)
    {
      tlist.AddToList(&(apply[i]));
      tlist.AddToList(&(pcopy[i]));
      tlist.AddToList(&(sigextract[i]));
      tlist.AddToList(&(mccalibupdate[i]));
      tlist.AddToList(&(calib[i]));
      tlist.AddToList(clean[i]);
      tlist.AddToList(&(hcalc[i]));
      tlist.AddToList(&(scalc[i]));
    }

  MStereoCalc stereocalc;
  stereocalc.SetCTids(CT[0],CT[1]);

  //
  // FIXME: telescope coordinates must be introduced by the user, since
  // they are not available nor in the camera file, neither in the reflector 
  // output.
  //

  if (CT[0] == 1)
    stereocalc.SetCT1coor(0.,0.);
  else if (CT[0] == 2)
    stereocalc.SetCT1coor(0.,0.);
  else if (CT[0] == 3)
    stereocalc.SetCT1coor(60.,60.);  // in meters
  else if (CT[0] == 4)
    stereocalc.SetCT1coor(60.,-60.);
  else if (CT[0] == 5)
    stereocalc.SetCT1coor(-60.,60.);
  else if (CT[0] == 6)
    stereocalc.SetCT1coor(-60.,-60.);
  else
    {
      cout << "Unknown CT id!" << endl;
      exit;
    }

  if (NCTs==2)
    {
      if (CT[1] == 1)
        stereocalc.SetCT2coor(0.,0.);
      else if (CT[1] == 2)
        stereocalc.SetCT2coor(0.,0.);
      else if (CT[1] == 3)
        stereocalc.SetCT2coor(60.,60.);  // in meters
      else if (CT[1] == 4)
        stereocalc.SetCT2coor(60.,-60.);
      else if (CT[1] == 5)
        stereocalc.SetCT2coor(-60.,60.);
      else if (CT[1] == 6)
        stereocalc.SetCT2coor(-60.,-60.);
      else
        {
          cout << "Unknown CT id!" << endl;
          exit;
        }

      tlist.AddToList(&stereocalc);
    }

  MF filter1("MMcEvt;1.fEvtNumber<2501");
  MF filter2("MMcEvt;1.fEvtNumber>2500");
  //
  // ^^^^ Filters to divide output in two: test and train samples.
  // FIXME: It is better to separate the events in odd and even 
  // event numbers (it is independent of the number of events in 
  // the file), but it is not yet possible because we cannot use
  // the modulus operator(%) in the filter yet.
  //

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

  tlist.AddToList(&filter1);
  tlist.AddToList(&write1);
  tlist.AddToList(&filter2);
  tlist.AddToList(&write2);

  //
  // Create and set up the eventloop
  //
  MProgressBar bar;

  MEvtLoop evtloop;
  evtloop.SetProgressBar(&bar);
  evtloop.SetParList(&plist);

  //
  // Execute your analysis
  //
  if (!evtloop.Eventloop())
    return;

  for (Int_t i= 0; i < NCTs; i++ )
    delete clean[i];

  tlist.PrintStatistics();
}