source: branches/Corsika7500Compatibility/macros/datatrigcheck.C@ 19921

/* ======================================================================== *\
!
! *
! * 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): Nicola Galante, 3/2004 <mailto:nicola.galante@pi.infn.it>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  datatrigcheck 
//
//  This macro is the standard trigger's parameters checking stuff
//  
//
/////////////////////////////////////////////////////////////////////////////
Double_t myPoisson(Double_t *x, Double_t *par)
{
  Double_t a=0.0, b=0.0, c=0.0;
  Double_t fitval = 0.0;

  Int_t floorX = (Int_t)TMath::Floor(x[0]);
  a = TMath::Power(par[0],floorX);
  b = TMath::Exp(-par[0]);
  c = TMath::Factorial(floorX);
 
  if (c!=0)
    fitval = par[1]*a*b/c;

  return fitval;
}

Double_t myExp(Double_t *x, Double_t *par)
{
  Double_t a=0.0;
  Double_t fitval = 0.0;

  a = TMath::Exp(-x[0]);
  fitval = par[0]*a;

  return fitval;
}

// The macro
//  Input:
//   1. The directory path starting from the Period (the complete path is
//       written in the basedir variable)
//   2. The data filename (without extension)
//
void datatrigcheck(TString dirin="Period014/rootdata/2004_02_17/", TString filein="20040117_03695_D_Crab-On_E")
{
 
  // Set file name and INPUT and OUTPUT directories
  
  TString basedir="/data1/earth/magic/data/";
  
  TString filenamein(filein);
  TString dirnamein(dirin);
  
  dirnamein = basedir+dirnamein;
  filenamein = dirnamein + filein + ".root";
  
  // Here the output plots are saved
  TString webdatadir = basedir + "webdata/" + filein + "/";   
  
  //
  // 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);
  
  MReadMarsFile read("Events");
  read.DisableAutoScheme();
  read.AddFile(filenamein);
  
  tlist.AddToList(&read);
  
  //
  // Extract Number of events
  //
  const Int_t numEvt = read.GetEntries();

  cout << "Opened file " << filenamein << " with "<< numEvt << " events." << endl;
  // ==================
  // Create histos
  // FIXME: should be created using MH or similar?
  
  // HISTcom-hTrigPatternEvt
  // The trigger pattern is shown as a decimal number
  // The trigger pattern consists of 16 bits (8+8 bits) generated by the trigger system. 
  // The first 8 bits correspond to the trigger configuration before the prescaling, 
  // the others after prescaling.
  //    Bit structure: 
  //          not prscd | prscaled
  //          xxxx xxxx xxxx xxxx    <-- pattern (x=0,1)
  //    bit   7654 3210 7654 3210
  //
    TH1I *hTrigPatternEvt = new TH1I("hTrigPatternEvt","Trigger Pattern vs. Event Number",numEvt,0,(Axis_t)numEvt);

    // HISTcom-hTimeEvt
    // Absolute time vs event number 
    // The Absolute time has been calculated using the MTime->GetTime() (millisec) and the MTime->GetMjd() (nanosec)
    TH1D *hTimeEvt = new TH1D("hTimeEvt","DAQ Time vs. Event Number",numEvt,0,(Axis_t)numEvt);
    
    // HISTcom-hTimeDiffEvt
    TH1D *hTimeDiffEvt = new TH1D("hTimeDiffEvt","Differences of time vs. Event Number",numEvt-1,0,(Axis_t)numEvt-1);
 
    // HISTcom-hTimeDiff
    TH1D *hTimeDiff = new TH1D("hTimeDiff","Distribution of differences of time (ns)",200,0,1000000);
    

    TF1 *func1 = new TF1("myPoisson",myPoisson,0,10,2);
    TF1 *func2 = new TF1("myExp",myExp,0,10,1);

    cout << "*func defined" << endl;
    func1->SetParameters(1.0,1.0);
    func2->SetParameter(0,400.0);
    cout << "parameters setted" << endl;
    func1->SetParNames("MU","CONST");
    func2->SetParNames("A");

    //
    // Create and set up the eventloop
    //
    MProgressBar bar;

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

    Double_t tOld;  // buffer variable
    Double_t tMin = 1E+10;
    Double_t tMax = -1E+10;
    Double_t DAQNumMin = 1E+10;
    //
    // Execute your analysis
    //
    //if (!evtloop.Eventloop())
    //    return;
    if (!evtloop.PreProcess())
    return;
    
    while (tlist.Process())
      {
        // Create the container MRawEvtHeader where time and trigger 
        //  info are read and stored.
        MRawEvtHeader *evtHeader = (MRawEvtHeader *)plist->FindObject("MRawEvtHeader");
        MTime *evtTime = (MTime *)plist->FindObject("MTime");

        Long_t ms = evtTime->GetTime();
        Double_t ns = ((evtTime->GetMjd()-TMath::Floor(evtTime->GetMjd()))*(Double_t)evtTime->kDay-(Double_t)ms)*1E+6;
        Double_t t = (Double_t)ms+ns/1.E+6;
        if (t<tMin) tMin=t;

        // Fill the histos
        hTrigPatternEvt->Fill(evtHeader->GetDAQEvtNumber(),evtHeader->GetTriggerID());
        hTimeEvt->Fill(evtHeader->GetDAQEvtNumber(),t);
        if (evtHeader->GetDAQEvtNumber() != 0){
          hTimeDiffEvt->Fill(evtHeader->GetDAQEvtNumber(),t-tOld);
          hTimeDiff->Fill((t-tOld)*1E+6);
        }
        if (t-tOld > tMax) tMax = t-tOld;
        tOld = t;
      }

    //TString command = "./macros/create_web_folder.sh " + webdatadir + filein ;

    // Create new directory with the same name of the run.
    // Delete old content.
    //
    TString command = "mkdir -p " + webdatadir ;
    system ( command.Data() ) ;
    command = "rm -f " + webdatadir + "*.gif" ;
    system ( command.Data() ) ;
    command = "rm -f " + webdatadir + "*.txt" ;
    system ( command.Data() ) ;

    tlist.PrintStatistics();

    // Draw Plots and save them as gif images
    //
    c1 = new TCanvas( "c1" );
    hTrigPatternEvt->SetXTitle("Event");
    hTrigPatternEvt->SetYTitle("Trigger Pattern");
    hTrigPatternEvt->Draw();
    c1->Print(webdatadir+hTrigPatternEvt->GetName()+".gif");

    c2 = new TCanvas( "c2" );
    hTimeEvt->SetXTitle("Event");
    hTimeEvt->SetYTitle("Abs. time (ms)");
    hTimeEvt->SetMinimum(tMin);
    hTimeEvt->Draw();
    c2->Print(webdatadir+hTimeEvt->GetName()+".gif");

    c3 = new TCanvas( "c3" );
    hTimeDiffEvt->SetXTitle("Event");
    hTimeDiffEvt->SetYTitle("Time diff (ms)");
    hTimeDiffEvt->SetMaximum(tMax*1.1);
    hTimeDiffEvt->Draw();
    c3->Print(webdatadir+hTimeDiffEvt->GetName()+".gif");

    c4 = new TCanvas( "c4" );
    hTimeDiff->SetXTitle("Time diff. (ms)");
    hTimeDiff->SetYTitle("Events");
    hTimeDiff->Draw();
    c4->Print(webdatadir+hTimeDiff->GetName()+".gif");

}