source: trunk/MagicSoft/Mars/mtemp/mpisa/macros/timedist3.C@ 5138

/* ======================================================================== *\
!
! *
! * 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 timedist3(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);

  MArrivalTimeCam 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);
  MRawEvtData *evtData;
  //
  // Create and setup the eventloop
  //
  MEvtLoop evtloop4;
  evtloop4.SetParList(&plist4);
  histTime = new TH1F("histTime","Arrival Time Dist",15,0.0,15.0);
  Int_t i=0;
  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())
    {
      i++;
      //disp1.SetCamContent(nphot);
      
      //gPad->Modified();
      //gPad->Update();
      //      cout << i << endl;
      //if((i == 16)||(i == 41)||(i == 111)||(i == 173)){
      if((i == 41)){
        cout << "Merda" << endl;
        evtData = (MRawEvtData *)plist4->FindObject("MRawEvtData");
        cout << "Merda1" << tlist4->FindObject("MRawEvtData") << endl;
        MRawEvtPixelIter *pixIter = new MRawEvtPixelIter(evtData);
        cout << "Merda2" << endl;
        for (Int_t pix=1;pix<576;pix++){          
          pixIter->Next();
          cout << "Merda3" << endl;
          if(pixIter->GetSumHiGainSamples()>200){
            cout << "Merda5" << endl;
            cout << "pix #" << pixIter->GetPixelId() << "  FADC Charge = " << pixIter->GetSumHiGainSamples() << endl;
            histTime->Fill((times[pix])->GetArrivalTimeHiGain());
          }
          cout << i << endl;
        }
      }
      
      /*    
      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();
  histTime->DrawCopy();
}

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