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

/////////////////////////////////////////////////////////////////////////////
//
//  mmcCleaning - calibrate/clean mc events
//
/////////////////////////////////////////////////////////////////////////////


void mmcCleaning(TString NoiseFilename, TString NoiseOutFilename)
{
  //
  // This is a demonstration program which reads in MC camera files
  // and produces and output with calibrated events (signal in photons
  // for all pixels, in MCerPhotEvt containers).
  //

  // ------------- user change -----------------
  Char_t*  NonoiseFilename = "/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0.root";
  //  Char_t*  NoiseFilename   = "/mnt/wdflix/root_0.73mirror/wuerzburg/gammas/Gamma_zbin3_0_7_1160to1169_w0.root"; // File to be analyzed

  //  Char_t* NoiseOutFilename   = "/mnt/users/jrico/prueba.root";   // Output file name

  const Int_t wsize=2;
  // Signal sigextract
  // (other extraction methods can be used)
  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);

  MPedestalCam  pedcam;
  MGeomCamMagic geomcam;
  MCerPhotEvt   nphot;

  plist.AddToList(&pedcam);
  plist.AddToList(&geomcam);
  plist.AddToList(&nphot);

  //
  // Now setup the tasks and tasklist:
  // ---------------------------------
  //
  MReadMarsFile read("Events");
  read.AddFile(NonoiseFilename);
  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.

  MDisplay       display(&nphot,&geomcam);
  MHillasDisplay display2(&nphot,&geomcam);

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

  MMcCalibrationUpdate  mccalibupdate;

  MCalibrate calib; 
  // MCalibrate 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(MCalibrate::kFfactor);

  MImgCleanStd clean(4.0,3.5);
  clean.SetCleanRings(1);
  //  clean.SetRemoveSingle(kFALSE);
  //
  // Applies tail cuts to image. Since the calibration is performed on 
  // noiseless camera files, the precise values of the cleaning levels 
  // are unimportant (in any case, only pixels without any C-photon will
  // be rejected).
  //

  MHillasCalc hcalc; // Calculates Hillas parameters not dependent on source position.
  MHillasSrcCalc hsrccalc;

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

  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(&hcalc);
  tlist.AddToList(&hsrccalc); 
  tlist.AddToList(&mccalibcalc);
  //tlist.AddToList(&display2);

  //
  // First loop: No noise
  //
  MProgressBar bar;
  bar.SetWindowName("No noise...");

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

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


  //
  // Second loop: process file with noise
  //  
  MIslands      isl;
  MArrivalTimeCam timecam;      
  MTopology     topology;
  plist.AddToList(&isl);
  plist.AddToList(&timecam);
  plist.AddToList(&topology);

  MArrivalTimeCalc2 timecalc;
  MIslandsCalc      islandcalc;
  islandcalc.SetOutputName("MIslands");
  islandcalc.SetAlgorithm(1);
  MIslandsClean     islclean(40);
  islclean.SetInputName("MIslands");
  islclean.SetMethod(1);
  MTopologyCalc     topcalc;

  MReadMarsFile read2("Events");
  read2.AddFile(NoiseFilename);
  read2.DisableAutoScheme();
  tlist.AddToListBefore(&read2, &read, "All");
  tlist.RemoveFromList(&read);

  tlist.AddToListBefore(&timecalc,&mccalibupdate,"All");
  tlist.AddToListBefore(&islandcalc,&hcalc,"All");
  tlist.AddToListBefore(&topcalc,&hcalc,"All");

  MWriteRootFile write(NoiseOutFilename); // Writes output
  write.AddContainer("MRawRunHeader", "RunHeaders");
  write.AddContainer("MMcRunHeader",  "RunHeaders");
  write.AddContainer("MSrcPosCam",    "RunHeaders");
  write.AddContainer("MGeomCam",      "RunHeaders");
  write.AddContainer("MMcConfigRunHeader",  "RunHeaders");
  write.AddContainer("MMcCorsikaRunHeader", "RunHeaders");
  write.AddContainer("MMcFadcHeader",  "RunHeaders");
  write.AddContainer("MMcTrigHeader",  "RunHeaders");

  write.AddContainer("MMcEvt",        "Events");
  write.AddContainer("MRawEvtHeader", "Events");
  write.AddContainer("MHillas",       "Events");
  write.AddContainer("MHillasSrc",    "Events");
  write.AddContainer("MHillasExt",    "Events");
  write.AddContainer("MConcentration","Events");
  write.AddContainer("MImagePar",     "Events");
  write.AddContainer("MNewImagePar",  "Events");
  write.AddContainer("MIslands",      "Events");
  write.AddContainer("MPointingPos",  "Events");
  write.AddContainer("MTopology",      "Events");

  tlist.RemoveFromList(&mccalibcalc);
  tlist.AddToList(&write);
  //  tlist.AddToListBefore(&display2,&write, "All");

  bar.SetWindowName("Calibrating/Cleaning...");
  //  clean.SetRemoveSingle();
  //  clean.SetMethod(MImgCleanStd::kDemocratic);

  if (!evtloop.Eventloop())
    return;

  tlist.PrintStatistics();
}