Demonstrator_RMS1.C

Example of a ROOT macro that computes the average samples and the RMS for all channels, from aTPC demonstrator run file.

*  [ROOT]  .L Demonstrator_RMS1.C
*  [ROOT]  PlotDemonstrator_RMS1 ( "../../data/Demonstrator/run_baseline.dat", 0, 0, 1, 17 );
*  

The processing is done in the macro itself, from channels output content.

//======================================================================
//  ROOT macros for TPC demonstrator system
//    - read a test file (2 CoBo modules / 8 AsAd)
//    - compute mean sample for all channels
//    - compute RMS sample for all channels
//    - display result for 1 channel
//
//======================================================================
/*
These macros require to load the GET libraries
  gROOT->Macro ( (getenv("GET_BASE")+string("/root/GETRootLibs.C")).c_str() );

To execute the macro
  .L Demonstrator_RMS1.C
  PlotDemonstrator_RMS1 ( "../../data/Demonstrator/run_baseline.dat", 0, 0, 1, 17 );
*/

//----------------------------------------------------------------------
//  This function reads all events from a TPC demonstrator run file
//  (a baseline run).
//  For all channels, it creates the average signal.
//  The the file is read a second time to compute the RMS of the
//  signal.
//
//  The average and the RMS are computed directly in the macro.
//
//  The result is plotted for one channel.
void PlotDemonstrator_RMS1 ( const string & data_file,
                             u_short cobo_id,
                             u_short asad_id,
                             u_short aget_id,
                             u_short chan_id )
{
  GSetVerboseLevel(3);

  // create the GET system
  u_short ncobo = 2;        // 2 CoBo modules (TPC demonstrator)
  u_int   n     = 512;      // 512 samples / channel
  double  dt    = 0.01;     // 0.01 us sampling time
  GETSystem gDem(ncobo,n,dt);

  u_int nch_tot = gDem.GetTotalChannelCount();


  // create GET systems to store the mean and RMS samples
  GETSystem gDemMean(ncobo,n,dt);
  GETSystem gDemRMS (ncobo,n,dt);

  //------------------------------------------------------------
  // Compute the average for all channels
  cout << endl
       << "Computing averages" << endl;

  u_int nr = 0;   // number of events read

  gDem.OpenRunFile ( data_file );
  while ( gDem.ReadEvent ( false ) == 0 )
  {
    nr++;   // increment the number of events

    for ( u_short ic = 0; ic < ncobo; ++ic )    // loop on CoBo modules
    {
      for ( u_short ia = 0; ia < 4; ++ia )        // loop on AsAd boards
      {
        for ( u_short ig = 0; ig < 4; ++ig )        // loop on AGet chips
        {
          for ( u_short ich = 0; ich < 68; ++ich )     // loop on channels
          {
            gDemMean[ic][ia][ig][ich].OutSample() += gDem[ic][ia][ig][ich].OutSample();
          } // end channels loop
        } // end AsAd loop
      } // end AsAd loop
    } // end CoBo loop
  }
  gDem.CloseRunFile ( );

  cout << "  - Number of events read: " << nr << endl
       << endl;

  if (nr > 0)
  {
    // divide by the number of events
    for ( u_int i = 0; i < nch_tot; ++i )
    {
      gDemMean.GetChannel(i)->OutSample() *= (1./nr);
    }

  //------------------------------------------------------------
  // Compute the RMS for all channels
    cout << "Computing RMS"  << endl
         << endl;

    nr = 0;
    gDem.OpenRunFile ( data_file );
    while ( gDem.ReadEvent ( false ) == 0 )
    {
      nr++;
      for ( u_int i = 0; i < nch_tot; ++i )
      {
        GETSample dif = gDemMean.GetChannel(i)->OutSample() - gDem.GetChannel(i)->OutSample();

        gDemRMS.GetChannel(i)->OutSample() += dif*dif;
      }
    }
    gDem.CloseRunFile ( );

    // square root and event number normalization
    for ( u_int i = 0; i < nch_tot; ++i )
    {
      GETSample & rms  = gDemRMS.GetChannel(i)->OutSample();
      double    * data = rms.GetFunctionData(false);

      for ( u_short j = 0; j < n; ++j )
      {
        double val = data[j];
        rms.SetFunctionValue ( j, sqrt(val/(nr-1.)) );
      }

      // see RRealSampleFFT (GETSample) documentation
      //      (not needed here since we used SetFunctionValue)
      rms.FunctionSet ( true );
    }


  //------------------------------------------------------------
  // Plot the results for one channel
    GETSample sample_mean = gDemMean[cobo_id][asad_id][aget_id][chan_id].OutSample();;
    GETSample sample_rms  = gDemRMS [cobo_id][asad_id][aget_id][chan_id].OutSample();;

    // compute the mean sample +/- the RMS
    GETSample sample_inf = sample_mean - sample_rms;
    GETSample sample_sup = sample_mean + sample_rms;

    // the GETSample is a RRealSampleFFT object (see GFFT documentation)
    TGraph * gr_mean = sample_mean.CreateFunctionGraph ( );
    gr_mean->SetTitle ( "Average + RMS channel output" );
    gr_mean->SetLineColor ( kRed+1 );
    gr_mean->SetLineWidth ( 2 );
    gr_mean->SetLineStyle ( 1 );
    gr_mean->GetXaxis()->SetRangeUser ( 0., (n-1)*dt );

    TGraph * gr_inf = sample_inf.CreateFunctionGraph ( );
    gr_inf->SetLineColor ( kGreen+1 );
    gr_inf->SetLineWidth ( 1 );
    gr_inf->SetLineStyle ( 1 );

    TGraph * gr_sup = sample_sup.CreateFunctionGraph ( );
    gr_sup->SetLineColor ( kGreen+1 );
    gr_sup->SetLineWidth ( 1 );
    gr_sup->SetLineStyle ( 1 );

    TCanvas * canvas = new TCanvas ( "PlotRMS", "PLOTRMS", 800, 400 );
    gr_mean->Draw ( "AL" );
    gr_inf->Draw ( "L" );
    gr_sup->Draw ( "L" );

    // create an EPS file
    gSystem->mkdir ( "plots", kTRUE );
    canvas->SaveAs ( "plots/Demonstrator_RMS1.eps" );
  }
}