source: trunk/MagicSoft/Mars/mtemp/mifae/macros/makeHillasMC.C@ 6723

/* ======================================================================== *\
   !
   ! *
   ! * 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): Thomas Bretz      5/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
   !              Abelardo Moralejo 1/2004 <mailto:moralejo@pd.infn.it>
   !              Eva Domingo       1/2005 <mailto:domingo@ifae.es>
   !             
   !   Copyright: MAGIC Software Development, 2000-2004
   !
   !
   \* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  STARMC - STandard Analysis and Reconstruction (MC example)
//
//  This macro is a version of the standard converter to convert raw data 
//  into image parameters, made to show how to run analysis on MC files. 
//
//  (1/2005):: Spline extractor added and additional containers written out
//
/////////////////////////////////////////////////////////////////////////////

#include "MImgCleanStd.h"

void makeHillasMC()
{
  //
  // This is a demonstration program which calculates the image 
  // parameters from Magic Monte Carlo files (output of camera).


  TString* CalibrationFilename;
  TString* AnalysisFilename;
  TString* OutFilename1;
  TString* OutFilename2;

  TString  rffilename;
  TString  dispfilename;

  // ------------- USER CHANGE -----------------
  // Comment line starting "CalibrationFileName" to disable calibration. 
  // In that case the units of the MHillas.fSize parameter will be ADC counts 
  // (rather, equivalent ADC counts in inner pixels, since we correct for the 
  // possible differences in gain of outer pixels)

  // File to be used in the calibration (must be a camera file without added noise)
  CalibrationFilename = new TString("/mnt/local_wdjrico/jrico/mc/Gammas/Gamma_zbin3_90_7_1290to1299_w0_nonoise.root");
  // File to be analyzed
  AnalysisFilename = new TString("/mnt/users/blanch/magic/TestSample/file*.root");

  // Change output file names as desired. 
  // If you want only one output (not dividing the events in Train and Test samples),
  // comment the initialization of OutFilename2.
  OutFilename1 = new TString("/mnt/users/blanch/Temp/Prova.root");
  //  OutFilename2 = new TString("star_test.root");    // Output file name 2 (train)

  // File to read RandomForest
  rffilename="/mnt/users/blanch/magic/Mars-All/Mars_Standard06/mtemp/mifae/programs/RFstd.root";
  // File to read Disp
  dispfilename="/mnt/users/blanch/magic/Mars-All/Mars_Standard06/mtemp/mifae/programs/DISPstd.root";

  // Fraction of events (taken at random) which one wants to process from the 
  // file to be analyzed (useful to make smaller files if starting sample is too large).
  Float_t accepted_fraction = 1.;


  // ------------- USER CHANGE -----------------
  // Cleaning levels
  Float_t CleanLev[2] = {4.0, 3.5}; // Tail cuts for image analysis

  // Signal extractor
  // (default=Spline, other extraction methods can be used)
  const Int_t wsize=2;
  MExtractor* sigextract = new MExtractTimeAndChargeSpline();
  ((MExtractTimeAndChargeSpline*)sigextract)->SetTimeType(MExtractTimeAndChargeSpline::kHalfMaximum);
  ((MExtractTimeAndChargeSpline*)sigextract)->SetChargeType(MExtractTimeAndChargeSpline::kIntegral);      
  ((MExtractTimeAndChargeSpline*)sigextract)->SetRiseTime((Float_t)wsize*0.25); 
  ((MExtractTimeAndChargeSpline*)sigextract)->SetFallTime((Float_t)wsize*0.75);  
  // Define FADC slices to be integrated in high and low gain:
  sigextract->SetRange(1, 11, 2, 12);
  // Defines when to switch to low gain
  sigextract->SetSaturationLimit(240);


  // -------------------------------------------
  //
  // 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);

  //  MGeomCamMagic geomcam;

  //  MSrcPosCam src;
  // WOBBLE MODE!!
  //  src.SetX(120.);  // units: mm
  //  src.SetXY(0.,0.4/geomcam.GetConvMm2Deg());
  //  src.SetXY(0.,0.); 
  //  src.SetReadyToSave();

  MBadPixelsCam badpix;

  //  plist.AddToList(&geomcam);
  //  plist.AddToList(&src);
  plist.AddToList(&badpix);

  MCerPhotEvt     nphot;
  // Stores number of phe
  plist.AddToList(&nphot);
  
  // Now setup the tasks and tasklist:
  // ---------------------------------
  //
  MReadMarsFile read("Events");
  if (CalibrationFilename)
    read.AddFile(CalibrationFilename->Data());
  read.DisableAutoScheme();

  MGeomApply geom; 
  // Reads in geometry from MC file and sets the right sizes for
  // several parameter containers.

  MMcPedestalCopy   pcopy; 
  // Copies pedestal data from the MC file run fadc header to the MPedestalCam container.

  MPointingPosCalc pointcalc;
  // Creates MPointingPos object and fill it with the telescope orientation
  // information taken from MMcEvt.

  MMcCalibrationUpdate  mccalibupdate;

  MCalibrateData calib; 
  // Transforms signals from ADC counts into photons. In the first
  // loop it applies a "dummy" calibration supplied by MMcCalibrationUpdate, just 
  // to equalize inner and outer pixels. At the end of the first loop, in the
  // PostProcess of MMcCalibrationCalc (see below) the true calibration constants
  // are calculated.
  calib.SetCalibrationMode(MCalibrateData::kFfactor);

  //    MBlindPixelCalc   blind;
  //    blind.SetUseInterpolation();

  MMcCalibrationCalc mccalibcalc; 
  // Calculates calibration constants to convert from ADC counts to photons.

  MAddGainFluctuation gainfluc;
  //gainfluc.FillHist(1,0.5);
  gainfluc.FillHist(0,0.1); // defaul line to not add fluctuations
  // Adds Gain fluctuations

  MImgCleanStd  clean(CleanLev[0], CleanLev[1]); 
  // Applies tail cuts to image. Since the calibration is performed on 
  // noiseless camera files, the precise values of the cleaning levels 
  // are unimportant for the calibration file processing (in any case, 
  // only pixels without any C-photon will be rejected).
  clean.SetCleanRings(1);
  //  clean.SetRemoveSingle(kFALSE);
  //  clean.SetMethod(MImgCleanStd::kDemocratic);

  MHillasCalc   hcalc; 
  // Calculates Hillas parameters.

  tlist.AddToList(&read);
  tlist.AddToList(&geom);
  tlist.AddToList(&pcopy);
  tlist.AddToList(&pointcalc);
  tlist.AddToList(sigextract);
  tlist.AddToList(&mccalibupdate);
  tlist.AddToList(&calib);
  tlist.AddToList(&clean);
  //  tlist.AddToList(&blind);
  tlist.AddToList(&hcalc);
  tlist.AddToList(&mccalibcalc);


  //
  // Open output files:
  //

  MWriteRootFile write1(OutFilename1->Data()); // Writes output1
  write1.AddContainer("MRawRunHeader",       "RunHeaders");
  write1.AddContainer("MMcRunHeader",        "RunHeaders");
  write1.AddContainer("MSrcPosCam",          "RunHeaders");
  write1.AddContainer("MGeomCam",            "RunHeaders");
  write1.AddContainer("MMcConfigRunHeader",  "RunHeaders");
  write1.AddContainer("MMcCorsikaRunHeader", "RunHeaders");
  write1.AddContainer("MMcFadcHeader",       "RunHeaders");
  write1.AddContainer("MMcTrigHeader",       "RunHeaders");
  write1.AddContainer("MGainFluctuationCam", "RunHeaders");

  write1.AddContainer("MRawEvtHeader",   "Events");
  write1.AddContainer("MMcEvt",          "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");
  write1.AddContainer("MIslands",        "Events");
  write1.AddContainer("MTopology",       "Events");
  write1.AddContainer("MPointingPos",    "Events");
  write1.AddContainer("MHadronness",     "Events");
  write1.AddContainer("MImageParDisp"  , "Events");
  //write1.AddContainer("MArrivalTimeCam", "Events");
  //write1.AddContainer("MCerPhotEvt",     "Events");

  if (OutFilename2)
    {
      MWriteRootFile write2(OutFilename2->Data()); // Writes output2
      write2.AddContainer("MRawRunHeader",       "RunHeaders");
      write2.AddContainer("MMcRunHeader",        "RunHeaders");
      write2.AddContainer("MSrcPosCam",          "RunHeaders");
      write2.AddContainer("MGeomCam",            "RunHeaders");
      write2.AddContainer("MMcConfigRunHeader",  "RunHeaders");
      write2.AddContainer("MMcCorsikaRunHeader", "RunHeaders");
      write2.AddContainer("MMcFadcHeader",       "RunHeaders");
      write2.AddContainer("MMcTrigHeader",       "RunHeaders");
      write2.AddContainer("MGainFluctuationCam", "RunHeaders");

      write2.AddContainer("MRawEvtHeader",   "Events");
      write2.AddContainer("MMcEvt",          "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");
      write2.AddContainer("MIslands",        "Events");
      write2.AddContainer("MTopology",       "Events");
      write2.AddContainer("MPointingPos",    "Events");
      write2.AddContainer("MHadronness",     "Events");
      write2.AddContainer("MImageParDisp"  , "Events");
      //write2.AddContainer("MArrivalTimeCam", "Events");
      //write2.AddContainer("MCerPhotEvt",     "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);
    }


  //
  // FIRST LOOP: Calibration loop (with no noise file)
  //
  MProgressBar bar;
  bar.SetWindowName("Calibrating...");

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

  if (CalibrationFilename)
    {
      if (!evtloop.Eventloop())
        return;
      mccalibcalc.GetHistADC2PhotEl()->Write();
      mccalibcalc.GetHistPhot2PhotEl()->Write();
      // Writes out the histograms used for calibration.
    }

  //
  // SECOND LOOP: TO READ THE RANDOM FOREST FILE(S) 
  //

  MParList   plistrf;
  MTaskList  tlistrf;
  MRanForest ranforest;
  MRanTree   rantree;

  if(rffilename.Length())
    {
      plistrf.AddToList(&tlistrf);
      plistrf.AddToList(&ranforest);
      plistrf.AddToList(&rantree);

      MReadTree readrf("Tree",rffilename);
      readrf.DisableAutoScheme();
      
      MRanForestFill rffill;
      rffill.SetNumTrees(100);
      
      tlistrf.AddToList(&readrf);
      tlistrf.AddToList(&rffill);
      
      MEvtLoop evtlooprf;
      evtlooprf.SetParList(&plistrf);
      if (!evtlooprf.Eventloop())
        return;
      
      tlistrf.PrintStatistics();
    }
  

  //
  // THIRD LOOP: Analysis loop (process file with noise)
  //
  bar.SetWindowName("Analyzing...");

  MIslands         isl;
  MArrivalTimeCam  timecam;      
  MTopology        topology;
  plist.AddToList(&isl);
  plist.AddToList(&timecam);
  plist.AddToList(&topology);
  plist.AddToList(&rantree);
  plist.AddToList(&ranforest);
 
  MArrivalTimeCalc2  timecalc;
  // Calculates the arrival time as the mean time of the fWindowSize 
  // time slices which have the highest integral content.
  //  MArrivalTimeCalc timecalc;            
  // Calculates the arrival times of photons. Returns the absolute 
  // maximum of the spline that interpolates the FADC slices.

  MIslandsCalc       islcalc;
  islcalc.SetOutputName("MIslands");
  islcalc.SetAlgorithm(1);
  //  MIslandsClean      islclean(40);
  //  islclean.SetInputName("MIslands");
  //  islclean.SetMethod(1);

  MTopologyCalc      topcalc;

  // Change the read task by another one which reads the file we want to analyze:
  MReadMarsFile read2("Events","/mnt/users/blanch/magic/TestSample/file*.root");
  //  MReadMarsFile read2("Events","/mnt/users/blanch/magic/TestSample/file53.root");
  //  read2.AddFile(AnalysisFilename->Data());
  read2.DisableAutoScheme();
  tlist.AddToListBefore(&read2, &read);
  tlist.RemoveFromList(&read);


  // Add new tasks (Islands and Topology calculations)
  tlist.AddToListBefore(&timecalc,&mccalibupdate);
  tlist.AddToListBefore(&islcalc,&hcalc);
  //  tlist.AddToListBefore(&islclean,&hcalc);
  tlist.AddToListBefore(&topcalc,&hcalc);

  // Add Task to Add Gain Fluctuations
  tlist.AddToListBefore(&gainfluc, &hcalc);

  MGeomCamMagic geomcam;// = new MGeomCam();
  plist.AddToList(&geomcam);


  MHillasDisplay*  disphillas=NULL;
  disphillas = new MHillasDisplay(&nphot,&geomcam);

  // Analyze only the desired fraction of events, skip the rest:
  MFEventSelector eventselector;
  Float_t rejected_fraction = 1. - accepted_fraction; 
  eventselector.SetSelectionRatio(rejected_fraction);
  MContinue skip(&eventselector);
  tlist.AddToListBefore(&skip, sigextract);

  // Remove calibration task from list.
  tlist.RemoveFromList(&mccalibcalc); 

  // disp
  // (read in optimum Disp parameter values)  
  MDispParameters dispin;
  TArrayD dispPar;
  if(dispfilename.Length())
    {
      TFile inparam(dispfilename);
      dispin.Read("MDispParameters");
      cout << "Optimum parameter values taken for calculating Disp : " << endl;

      dispPar =  dispin.GetParameters();
      for (Int_t i=0; i<dispPar.GetSize(); i++)
        cout << dispPar[i] << ",  ";
      cout << endl;
      
      inparam.Close();
    }
  // Disp results container
  MImageParDisp imagepardisp;
  // Disp histograms
  MHDisp hdisp;
  hdisp.SetName("MHDispAsym");
  hdisp.SetSelectedPos(4);
  MFillH filldisp("MHDispAsym[MHDisp]", "");
  
  // Add RandomForest and Disp task computation
  MRanForestCalc hadrcalc;
  MDispCalc         dispcalc;  

  tlist.AddToList(&dispcalc);
  tlist.AddToList(&filldisp);
  if(rffilename.Length())
    tlist.AddToList(&hadrcalc);
 
  // Add tasks to write output:
  if (OutFilename2)
    {
      tlist.AddToList(&filter1);
      tlist.AddToList(&filter2);
      tlist.AddToList(&write2); 
    }
  tlist.AddToList(&write1); 

  // Add display image with hillas
  disphillas->SetPSFile();
  disphillas->SetPSFileName("provaGF.ps");
  disphillas->SetPause(kFALSE); 
  //tlist.AddToList(disphillas);
  

  if (!evtloop.Eventloop())
    return;

  tlist.PrintStatistics();
}