source: trunk/MagicSoft/Mars/mtemp/mpisa/macros/production3.C@ 6724

/* ======================================================================== *\
!
! *
! * 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): Javier López, 12/2003 <mailto:jlopez@ifae.es>
!              Markus Gaug , 04/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */
/////////////////////////////////////////////////////////////////////////////
//
//  bootcampstandardanalysis.C
//
//  Updated version of the macro designed at the Wuerzburg bootcamp and 
//  compatible with the latest changes in Mars for general usage at the
//  Udine bootcamp. 
//
//  Needs as arguments the run number of a pedestal file ("*_P_*.root"), 
//  one of a calibration file ("*_C_*.root") and one of a data file 
//  ("*_D_*.root"). Performs the pedestal calculation, the calibration 
/// constants calculation and the calibration of the data. 
//
//  The TString inpath has to be set correctly.
//
//  The macro searches for the pulser colour which corresponds to the calibration
//  run number. If the run number is smaller than 20000, pulser colour "CT1" 
//  is assumed, otherwise, it searches for the strings "green", "blue", "uv" or 
//  "ct1" in the filenames. If no colour or multiple colours are found, the 
//  execution is aborted.  
//
//////////////////////////////////////////////////////////////////////////////////
const Bool_t usedisplay = kTRUE;
const TString defpath = "/data1/earth/magic/data/Period015/rootdata/2004_03_19/";
const TString outpath = "./";
const TString defrout = "output_Mrk421_20040319.root";
// Mrk (13)
const Int_t defpedr[] = {20823};
//const Int_t defpedr [] = {20590,20593,20596,20599,20601,20602,20605,
//                        20608,20611,20614,20649,20652,20654};
// OffMrk (14)
//const Int_t defpedr [] = {20617,20619,20621,20623,20625,20628,20630,
//                        20632,20635,20637,20640,20643,20645,20647};


//const Int_t defcalr [] = {20651};
const Int_t defcalr [] = {20822};

// Mrk (20)
const Int_t defdatar[] = {20824};
//const Int_t defdatar[] = {20591,20592,20594,20595,20597,20598,20600,
//                        20603,20604,20606,20607,20609,20610,20612,
//                        20613,20615,20616,20650,20653,20655};
// OffMrk (18)
// const Int_t defdatar[] = {20618,20620,20622,20624,20626,20627,20629,
//                        20631,20633,20634,20636,20638,20639,20641,
//                        20642,20644,20646,20648};


void production3(const TString inpath=defpath, 
                    const Int_t psize=1, const Int_t pedruns[]=defpedr, 
                    const Int_t csize=1, const Int_t calruns[]=defcalr, 
                    const Int_t dsize=1, const Int_t dataruns[]=defdatar, 
                    const TString resname=defrout)

{

  //
  // Choose the signal Extractor:
  //
  //MExtractSlidingWindow extractor;
  MExtractFixedWindowPeakSearch extractor;
  //MExtractFixedWindow    extractor;   

  //
  // Set Ranges or Windows
  //
  extractor.SetRange(3,14,3,14);
  //extractor.SetWindows(8,8);

  //
  // Choose the arrival time Extractor:
  //
  MExtractTimeFastSpline       timeext;
  //  MExtractTimeHighestIntegral timeext;
  //  MExtractTimeSpline          timeext;
  //
  // Set Ranges or Windows
  //
  timeext.SetRange(2,12,4,14);


  MRunIter pruns;
  MRunIter cruns;
  MRunIter druns;

  for (Int_t i=0;i<psize;i++) {
    cout << "Adding pedestal run: " << pedruns[i] << endl;
    pruns.AddRun(pedruns[i],inpath);
  }
  for (Int_t i=0;i<csize;i++) {
    cout << "Adding calibration run: " << calruns[i] << endl;
    cruns.AddRun(calruns[i],inpath);
  }
  for (Int_t i=0;i<dsize;i++) {
    cout << "Adding data run: " << dataruns[i] << endl;
    druns.AddRun(dataruns[i],inpath);
  }

  
  //
  // Now setup the tasks and tasklist for the pedestals:
  // ---------------------------------------------------
  //
  MBadPixelsCam     badcam;
  MGeomCamMagic     geomcam;
  MGeomApply        geomapl;
  //
  // If you want to exclude pixels from the beginning, read 
  // an ascii-file with the corr. pixel numbers (see MBadPixelsCam)
  //
  //  badcam.AsciiRead("badpixels.dat");
  
  /************************************/
  /* FIRST LOOP: PEDESTAL COMPUTATION */
  /************************************/

  MJPedestal pedloop;
  pedloop.SetInput(&pruns);
  if (usedisplay)
    {
      MStatusDisplay *display = new MStatusDisplay;
      display->SetUpdateTime(3000);
      display->Resize(850,700);
      pedloop.SetDisplay(display);
    }
  pedloop.SetBadPixels(badcam);
  pedloop.SetExtractor(&extractor);

  if (!pedloop.Process())
    return;

  /****************************************/
  /* SECOND LOOP: CALIBRATION COMPUTATION */
  /****************************************/

  //
  // Now setup the new tasks for the calibration:
  // ---------------------------------------------------
  //
  MCalibrationQECam qecam;
  MJCalibration     calloop;
  calloop.SetInput(&cruns);
  calloop.SetExtractor(&extractor);
  //
  // Apply rel. time calibration:
  //
  calloop.SetRelTimeCalibration();
  calloop.SetTimeExtractor(&timeext);
  //
  // Set the corr. cams:
  //
  calloop.SetQECam(qecam);
  calloop.SetBadPixels(pedloop.GetBadPixels());
  //
  // The next two commands are for the display:
  //
  if (usedisplay)
      calloop.SetDisplay(display);
  cout << "Minchia 0" << endl;
  //
  // Do the event-loop:
  //
  if (!calloop.Process(pedloop.GetPedestalCam()))
    return;
  
  cout << "Minchia 1" << endl;
  /*************************************/
  /* THIRD LOOP: PEDESTAL CALIBRATION  */
  /*************************************/

  //
  // Create a empty Parameter List and an empty Task List
  //
  MParList  plist3;
  MTaskList tlist3;
  plist3.AddToList(&tlist3);
  cout << "Minchia 2" << endl;
  //
  // Now setup the tasks and tasklist to calculate the pedestal rms in number of photons
  // -----------------------------------------------------------------------------------
  //
  
  MCerPhotEvt    nphot;
  MPedPhotCam    nphotrms;
  
  plist3.AddToList(&geomcam);
  cout << "Minchia 3" << endl;
  //
  // Retrieve the cameras from the previous runs:
  //
  plist3.AddToList(&pedloop.GetPedestalCam());
  plist3.AddToList(&calloop.GetCalibrationCam());
  plist3.AddToList(&calloop.GetQECam());
  plist3.AddToList(&calloop.GetRelTimeCam());
  plist3.AddToList(&calloop.GetBadPixels());
  plist3.AddToList(&nphot);
  plist3.AddToList(&nphotrms);
    cout << "Minchia 4" << endl;
  //tasks
  MReadMarsFile read3("Events");
  read3.DisableAutoScheme();
  static_cast<MRead&>(read3).AddFiles(pruns);  

  //MCalibrate       photcalc;
  MCalibrateData       photcalc;
  photcalc.SetCalibrationMode(MCalibrate::kFfactor);
  MPedPhotCalc  photrmscalc;  //It doesn't exist yet
  
  tlist3.AddToList(&read3);
  tlist3.AddToList(&geomapl);
  tlist3.AddToList(&extractor);
  tlist3.AddToList(&photrmscalc);
  //tlist3.AddToList(&photcalc);
  
  
  //
  // Create and setup the eventloop
  //
  MEvtLoop evtloop3;
  evtloop3.SetParList(&plist3);
  
  //
  // Execute first analysis
  //
  if (!evtloop3.Eventloop())
    return;
  
  tlist3.PrintStatistics();
  
  /*************************************/
  /* FOURTH LOOP: DATA CALIBRATION     */
  /*************************************/

  //
  // Create a empty Parameter List and an empty Task List
  //
  MParList  plist4;
  MTaskList tlist4;
  plist4.AddToList(&tlist4);
  
  //
  // Now setup the tasks and tasklist to analize the data
  // -----------------------------------------------------
  //
  
  plist4.AddToList(&geomcam);
  //
  // Retrieve the cameras from the previous runs:
  //
  plist4.AddToList(&pedloop.GetPedestalCam());
  plist4.AddToList(&calloop.GetCalibrationCam());
  plist4.AddToList(&calloop.GetQECam());
  plist4.AddToList(&calloop.GetRelTimeCam());
  plist4.AddToList(&calloop.GetBadPixels());
  plist4.AddToList(&nphot);
  plist4.AddToList(&nphotrms);

  MArrivalTime times;
  plist4.AddToList(&times);
  
  //tasks
  MReadMarsFile read4("Events");
  read4.DisableAutoScheme();
  static_cast<MRead&>(read4).AddFiles(druns);  
  
  MArrivalTimeCalc2   timecalc;
  MCalibrateRelTimes  timecal;

  tlist4.AddToList(&read4);
  tlist4.AddToList(&geomapl);
  tlist4.AddToList(&extractor);
  tlist4.AddToList(&timecalc);
  tlist4.AddToList(&photcalc);
  tlist4.AddToList(&timecal);
  
  //
  // Create and setup the eventloop
  //
  MEvtLoop evtloop4;
  evtloop4.SetParList(&plist4);

  if (!evtloop4.PreProcess())
    return;
  
  TCanvas *c1 = new TCanvas;
  MHCamera disp1(geomcam);
  disp1.SetPrettyPalette();
  //disp1.SetInvDeepBlueSeaPalette()
  disp1.Draw();
  gPad->SetLogy();
  gPad->cd(1);
  
  /*
  TCanvas *c2 = new TCanvas;
  MHCamera disp2(geomcam);
  disp2.SetPrettyPalette();
  //disp2.SetInvDeepBlueSeaPalette()
  disp2.Draw();
  gPad->SetLogy();
  gPad->cd(1);
  */
  while (tlist4.Process())
    {
      disp1.SetCamContent(nphot);
      
      gPad->Modified();
      gPad->Update();

      /*    
      disp2.SetCamContent(times);
      
      gPad->Modified();
      gPad->Update();
      */

      // Remove the comments if you want to go through the file
      // event-by-event:
      if (!HandleInput())
        break;
    }

  evtloop4.PostProcess();
  
}

Bool_t HandleInput()
{
    TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);

    while (1)
    {
        //
        // While reading the input process gui events asynchronously
        //
        timer.TurnOn();
        TString input = Getline("Type 'q' to exit, <return> to go on: ");
        timer.TurnOff();

        if (input=="q\n")
            return kFALSE;

        if (input=="\n")
            return kTRUE;
    };

    return kFALSE;
}