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

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

//
// User change. 
//

Float_t ctx[7] = {0., 0., 0., 0., 0., 0., 0.}; 
Float_t cty[7] = {-70., -40., -30., 30., 50., 60., 70.}; // in meters
//
// FIXME: unfortunately present version of reflector was not prepared for
// stereo configurations and keeps no track of CT position. So the positions
// must be set above by the user, making sure that they correspond to the
// files one is analysing.
//

void starmcstereo(Int_t ct1 = 2, Int_t ct2 = 5)
{
  if (ct1 > sizeof(ctx)/sizeof(ctx[0]) || 
      ct2 > sizeof(ctx)/sizeof(ctx[0]) )
    {
      cout << endl << "Wrong CT id number!" << endl;
      return;
    }

  Int_t CT[2] = {ct1, ct2};  // Only 2-telescope analysis for the moment
  Int_t NCTs = sizeof(CT)/sizeof(CT[0]);


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

  Char_t* AnalysisFilename = "gam-yXX-00001.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] = {5, 9}; // First and last FADC bin of the range to be integrated,
  Int_t BinsLow[2]  = {5, 9}; // 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);

  MSrcPosCam src[NCTs];
  MBadPixelsCam badpix[NCTs];

  for (Int_t ict = 0; ict < NCTs; ict++)
    {
      TString s = "MSrcPosCam;";
      s += CT[ict];
      src[ict].SetName(s);
      src[ict].SetReadyToSave();
      plist.AddToList(&(src[ict]));

      TString b = "MBadPixelsCam;";
      b += CT[ict];
      badpix[ict].SetName(b);
      badpix[ict].SetReadyToSave();
      plist.AddToList(&(badpix[ict]));
    }
  //
  // 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.
  //

  stereocalc.SetCT1coor(ctx[CT[0]-1],cty[CT[0]-1]);
  stereocalc.SetCT2coor(ctx[CT[1]-1],cty[CT[1]-1]);

  tlist.AddToList(&stereocalc);


  MF filter1("{MMcEvt;1.fEvtNumber%2}<0.5");
  MF filter2("{MMcEvt;1.fEvtNumber%2}>0.5");
  //
  // ^^^^ Filters to divide output in two: test and train samples.
  //

  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();
}