tof/src/TofScan.cpp

/**
* TOFScan
* Author  Nagni
* Version 1.2
* Modified by G.De Rosa
* Date 27 Apr 2006
* Modified by G.De Rosa
* Date 03 Jul 2006
* Modified by W. Menn to select helium particles for PMT gain check
* Date 09 Aug 2007
* Last version 08 Oct 2007
*
* Description:
* Describe the performance of the TOF.
*
* Parameters:
*       TString base            - the path to the root directory for the specific Pamela unpack session 
*       TString outDirectory    - the path where to save the output image (Default = base) 
*       TString format          - the format which will be used fo rsave the produced images (Default = "gif") 
*/

#include <TROOT.h>
#include <TH1.h>
#include <TFile.h>
#include <TObjArray.h>
#include <TString.h>
#include <TObjString.h>
#include <TTree.h>
#include <TBranch.h>
#include <TGraph.h>
#include <TStyle.h>
#include <TH2S.h>
#include <TPaveText.h>
#include <TCanvas.h>
#include <physics/tof/TofEvent.h>
#include <cstdlib>
#include <sys/stat.h>
#include <math.h>
#include <iostream>
#include <fstream> 

using namespace std;

void TofScan(TString base, TString outDirectory = "", TString format = ""){
   
   std::stringstream sst;
   if (outDirectory == "") outDirectory = base.Data();
   TString filename = ((TObjString*)base.Tokenize('/')->Last())->GetString();

   TFile *file  =new TFile(base.Data()) ;
   if (!file){
    printf("file not Found \n");
    return;
   }

   TTree *PhysicsTr = (TTree*)file->Get("Physics");
   TBranch *TofBr = PhysicsTr->GetBranch("Tof");
   pamela::tof::TofEvent *tofEvent = 0;
   PhysicsTr->SetBranchAddress("Tof", &tofEvent);
   
   Long64_t nevents    = TofBr->GetEntries();
   if (nevents <= 0) {
     printf("nevents = %llu \n", nevents);
     file->Close();
     return;
   }

   /*
    * Array to convert hdc/adc to the real Photomultiplier
    * The array rows definitions are:
    * tof[0][] = chxxA (strip or channel xxA)
    * tof[1][] = hbxxA (halfboard xxA)
    * tof[2][] = chxxB (strip or channel xxB)
    * tof[3][] = hbxxB (halfboard xxB)
    *
    * Each single row is a sequence of photomultipliers in this shape
    * - The elements from  0 to  7  correspond to S11_1->S11_8
    * - The elements from  8 to 13  correspond to S12_1->S12_6
    * - The elements from 14 to 15  correspond to S21_1->S21_2
    * - The elements from 16 to 17  correspond to S22_1->S22_2
    * - The elements from 18 to 20  correspond to S31_1->S31_3
    * - The elements from 21 to 23  correspond to S32_1->S32_3
    * 
    * Example: 
    * -------> the tdc of the S12_3B photomultiplier correspond to tdc[(tof[2][10])][(tof[3][10])] 
    * -------> the tdc of the S31_3A photomultiplier correspond to tdc[(tof[0][20])][(tof[1][20])] 
    */
   short tof[4][24] = {
     {4, 4,  4,  4,  1,  1, 2, 2,  3,  3, 3, 3,  3,  3, 1, 1,  1,  1, 2, 3,  3, 3, 3,  4},
     {1, 3,  5,  7, 10, 12, 2, 4,  2,  4, 6, 8, 10, 12, 1, 5,  3,  9, 7, 9, 11, 1, 5,  9},
     {2, 2,  2,  2,  1,  1, 1, 1,  4,  4, 4, 4,  4,  4, 2, 1,  2,  1, 2, 2,  2, 3, 3,  4},
     {6, 8, 12, 10,  8,  6, 4, 2, 12, 10, 8, 6,  4,  2, 9, 7, 11, 11, 5, 3,  1, 3, 7, 11}
   };
   
   TString photoS[48] = {
     "S11_1A", "S11_1B", "S11_2A", "S11_2B", "S11_3A", "S11_3B", "S11_4A", "S11_4B", 
     "S11_5A", "S11_5B", "S11_6A", "S11_6B", "S11_7A", "S11_7B", "S11_8A", "S11_8B", 
     "S12_1A", "S12_1B", "S12_2A", "S12_2B", "S12_3A", "S12_3B", "S12_4A", "S12_4B", "S12_5A",  "S12_5B", "S12_6A", "S12_6B",
     "S21_1A", "S21_1B", "S21_2A", "S21_2B", 
     "S22_1A", "S22_1B", "S22_2A", "S22_2B",
     "S31_1A", "S31_1B", "S31_2A", "S31_2B", "S31_3A", "S31_3B", 
     "S32_1A", "S32_1B", "S32_2A", "S32_2B", "S32_3A", "S32_3B"
   };
   
   const Int_t nh = 48;
   TH1F *htdc[nh];
   TH1F *hadc[nh];
   
   TObjArray *hhtdc = new TObjArray(nh);
   TObjArray *hhadc = new TObjArray(nh);
   char tdcname[48]="";
   char adcname[48]="";

   char htitle[50];
   TH1F *adche[48];
   for(int i=0;i<48;i++) {
   sprintf(htitle, "adche_%d",(i+1));
   adche[i] = new TH1F(htitle,htitle,100,0.,1500.);
                         }


   Float_t adca[48];  // vector with adc values according to "ind"=pmt_id 
   Float_t tdca[48];  // the same for tdc
   
   int j  = 0;
   int k  = 0;
   int z  = 0;
   int ch = 0;
   int hb = 0;
   int ind =0;

   int heevent =0;

// upper and lower limits  for the helium selection
Float_t A_l[24]={200,190,300,210,220,200,210,60, 60, 120,220,120,160,50, 300,200,120,250,350,300,350,250,280,300}; 
Float_t A_h[24]={550,490,800,600,650,600,600,260,200,380,620,380,550,200,850,560,400,750,900,800,880,800,750,800};

// The k1 constants for the beta calculation, only for S1-S3
// k2 constant is taken to be the standard 2D/c
Float_t k1[72] = {50,59.3296,28.4328,-26.0818,5.91253,-19.588,-9.26316,24.7544,2.32465,
-50.5058,-15.3195,-39.1443,-91.2546,-58.6243,-84.5641,-63.1516,-32.2091,
-58.3358,13.8084,45.5322,33.2416,-11.5313,51.3271,75,-14.1141,
42.8466,15.1794,-63.6672,-6.07739,-32.164,-41.771,10.5274,-9.46096,
-81.7404,-28.783,-52.7167,-127.394,-69.6166,-93.4655,-98.9543,-42.863,
-67.8244,-19.3238,31.1221,8.7319,-43.1627,5.55573,-14.4078,-83.4466,
-47.4647,-77.8379,-108.222,-75.986,-101.297,-96.0205,-63.1881,-90.1372,
-22.7347,8.31409,-19.6912,-7.49008,23.6979,-1.66677,1.81556,34.4668,
6.23693,-100,-59.5861,-90.9159,-141.639,-89.2521,-112.881} ;

//-------------------------------------------------------------------


   for (int i=0; i < nevents; i++){

     TofBr->GetEntry(i);

     k = 0;
     while (k < 24){
       j = 0;
       while (j < 2){
         ch = tof[2*j][k]     - 1;
         hb = tof[2*j + 1][k] - 1;
         ind = 2*k + j;

         if(i==0){
           sprintf(tdcname,"TDChist%4.4d",ind);
           sprintf(adcname,"ADChist%4.4d",ind);
       
           htdc[ind] = new TH1F(tdcname,tdcname,409,0,4096);
           hadc[ind] = new TH1F(adcname,adcname,409,0,4096);
       
           hhtdc->Add(htdc[ind]);
           hhadc->Add(hadc[ind]);
         }

         htdc[ind]->Fill(tofEvent->tdc[ch][hb]);
         hadc[ind]->Fill(tofEvent->adc[ch][hb]);
         tdca[ind]=tofEvent->tdc[ch][hb];
         adca[ind]=tofEvent->adc[ch][hb];
         j++;
       }
       k++;
     }

//============   calculate beta and select helium ====================

// find hitted paddle by looking for ADC values on both sides
// since we looking for helium this gives decent results

//Int_t tof11_i,tof12_i,tof21_i,tof22_i,tof31_i,tof32_i;
Int_t tof11_i,tof12_i,tof31_i,tof32_i;
Float_t  a1,a2;
Int_t jj;

//  reset values
tof11_i = -1;
tof12_i = -1;
//tof21_i = -1;
//tof22_i = -1;
tof31_i = -1;
tof32_i = -1;

    for(jj=0; jj<8; jj++){
    a1 = adca[2*jj];
    a2 = adca[2*jj+1];
    if ((a1 < 3000) && (a2 < 3000)) tof11_i = jj;
                           }
    for(jj=0; jj<6; jj++){
    a1 = adca[16+2*jj];
    a2 = adca[16+2*jj+1];
    if ((a1 < 3000) && (a2 < 3000)) tof12_i = jj;
                           }
    /*    for(jj=0; jj<2; jj++){
    a1 = adca[28+2*jj];
    a2 = adca[28+2*jj+1];
    //    if ((a1 < 3000) && (a2 < 3000)) tof21_i = jj;
                           }
    for(jj=0; jj<2; jj++){
    a1 = adca[32+2*jj];
    a2 = adca[32+2*jj+1];
    //    if ((a1 < 3000) && (a2 < 3000)) tof22_i = jj;
    }*/
    for(jj=0; jj<3; jj++){
    a1 = adca[36+2*jj];
    a2 = adca[36+2*jj+1];
    if ((a1 < 3000) && (a2 < 3000)) tof31_i = jj;
                           }
    for(jj=0; jj<3; jj++){
    a1 = adca[42+2*jj];
    a2 = adca[42+2*jj+1];
    if ((a1 < 3000) && (a2 < 3000)) tof32_i = jj;
                           }


//----------------------------------------------------------------

Float_t zin[6] = {53.74, 53.04, 23.94, 23.44, -23.49, -24.34};
Float_t  c1,c2,xhelp,xhelp1,xhelp2,ds,dist,F;
Float_t sw,sxw,beta_mean_tof,w_i;
Float_t  theta,x1,x2,y1,y2,dx,dy,dr;
Int_t  ihelp;
Int_t ipmt[4];
Float_t time[4];
Float_t beta1[4];

// Only use events with:  S11 and S12 and S31 and S32

if ( (tof11_i>-1) && (tof12_i>-1) && (tof31_i>-1) && (tof32_i>-1) ) {

//  calculate zenith angle theta using the locations of the hitted paddles


        Float_t tof11_x[8] = {-17.85,-12.75,-7.65,-2.55,2.55,7.65,12.75,17.85};
        Float_t tof12_y[6] = { -13.75,-8.25,-2.75,2.75,8.25,13.75};
//        Float_t tof21_y[2] = { 3.75,-3.75};
//        Float_t tof22_x[2] = { -4.5,4.5};
        Float_t tof31_x[3] = { -6.0,0.,6.0};
        Float_t tof32_y[3] = { -5.0,0.0,5.0};

        //  S11 8 paddles  33.0 x 5.1 cm
        //  S12 6 paddles  40.8 x 5.5 cm
        //  S21 2 paddles  18.0 x 7.5 cm
        //  S22 2 paddles  15.0 x 9.0 cm
        //  S31 3 paddles  15.0 x 6.0 cm
        //  S32 3 paddles  18.0 x 5.0 cm

      x1 = 0.;
      x2 = 0.;
      y1 = 0.;
      y2 = 0.;

      x1 = tof11_x[tof11_i] ;
      y1 = tof12_y[tof12_i] ;
      x2 = tof31_x[tof31_i] ;
      y2 = tof32_y[tof32_i] ;

      theta=0.;
      dx=0.;
      dy=0.;
      dr=0.;

      dx = x1-x2;
      dy = y1-y2;
      dr = sqrt(dx*dx+dy*dy);
      theta = atan(dr/77.5);


      beta_mean_tof=100.;

      for (Int_t jj=0; jj< 4; jj++) beta1[jj] = 100. ;


//----------------------------------------------------------------
//---------  S1 - S3 ---------------------------------------------
//----------------------------------------------------------------

//---------  S11 - S31 -------------------------------------------

if ((tof11_i>-1)&&(tof31_i>-1)) {

dist = zin[0] - zin[4];
c2 = (2.*0.01*dist)/(3.E08*50.E-12);
F = 1./cos(theta);

           ipmt[0]   = (tof11_i)*2;
           ipmt[1]   = (tof11_i)*2+1;
           ipmt[2]   = 36+(tof31_i)*2;
           ipmt[3]   = 36+(tof31_i)*2+1;

   for (Int_t jj=0; jj< 4; jj++) time[jj] = tdca[(ipmt[jj])] ;

   if ((time[0]<4095)&&(time[1]<4095)&&(time[2]<4095)&&(time[3]<4095)) {
         xhelp1 = time[0] + time[1] ;
         xhelp2 = time[2] + time[3] ;
         ds = xhelp1-xhelp2;
         ihelp=0+(tof11_i)*3+tof31_i ;
         c1 = k1[ihelp] ;
         beta1[0] = c2*F/(ds-c1);
}
}

//---------  S11 - S32 -------------------------------------------

if ((tof11_i>-1)&&(tof32_i>-1)) {

dist = zin[0] - zin[5];
F = 1./cos(theta);
c2 = (2.*0.01*dist)/(3.E08*50.E-12);

           ipmt[0]   = (tof11_i)*2;
           ipmt[1]   = (tof11_i)*2+1;
           ipmt[2]   = 42+(tof32_i)*2;
           ipmt[3]   = 42+(tof32_i)*2+1;

   for (Int_t jj=0; jj< 4; jj++) time[jj] = tdca[(ipmt[jj])] ;

   if ((time[0]<4095)&&(time[1]<4095)&&(time[2]<4095)&&(time[3]<4095)) {
         xhelp1 = time[0] + time[1] ;
         xhelp2 = time[2] + time[3] ;
         ds = xhelp1-xhelp2;
         ihelp=24+(tof11_i)*3+tof32_i ;
         c1 = k1[ihelp] ;
         beta1[1] = c2*F/(ds-c1);
}
}

//---------  S12 - S31 -------------------------------------------

if ((tof12_i>-1)&&(tof31_i>-1)) {

dist = zin[1] - zin[4];
F = 1./cos(theta);
c2 = (2.*0.01*dist)/(3.E08*50.E-12);

           ipmt[0]   = 16+(tof12_i)*2;
           ipmt[1]   = 16+(tof12_i)*2+1;
           ipmt[2]   = 36+(tof31_i)*2;
           ipmt[3]   = 36+(tof31_i)*2+1;

   for (Int_t jj=0; jj< 4; jj++) time[jj] = tdca[(ipmt[jj])] ;

   if ((time[0]<4095)&&(time[1]<4095)&&(time[2]<4095)&&(time[3]<4095)) {
         xhelp1 = time[0] + time[1] ;
         xhelp2 = time[2] + time[3] ;
         ds = xhelp1-xhelp2;
         ihelp=48+(tof12_i)*3+tof31_i ;
         c1 = k1[ihelp] ;
         beta1[2] = c2*F/(ds-c1);
}
}

//---------  S12 - S32 -------------------------------------------

if ((tof12_i>-1)&&(tof32_i>-1)) {

dist = zin[1] - zin[5];
F = 1./cos(theta);
c2 = (2.*0.01*dist)/(3.E08*50.E-12);

           ipmt[0]   = 16+(tof12_i)*2;
           ipmt[1]   = 16+(tof12_i)*2+1;
           ipmt[2]   = 42+(tof32_i)*2;
           ipmt[3]   = 42+(tof32_i)*2+1;

    for (Int_t jj=0; jj< 4; jj++) time[jj] = tdca[(ipmt[jj])] ;

   if ((time[0]<4095)&&(time[1]<4095)&&(time[2]<4095)&&(time[3]<4095)) {
         xhelp1 = time[0] + time[1] ;
         xhelp2 = time[2] + time[3] ;
         ds = xhelp1-xhelp2;
         ihelp=66+(tof12_i)*3+tof32_i ;
         c1 = k1[ihelp] ;
         beta1[3] = c2*F/(ds-c1);
}
}

//----------------------   calculate  beta mean -----------------

      sw=0.;
      sxw=0.;
      beta_mean_tof=100.;

        for (Int_t jj=0; jj<4;jj++){
        if ((beta1[jj]>0.1) && (beta1[jj]<1.5)) {
            w_i=1./(0.13*0.13);
            sxw=sxw + beta1[jj]*w_i ;
            sw =sw + w_i ;
         }
         }

      if (sw>0) beta_mean_tof=sxw/sw;

} //  if tof11_i > -1 && ...... beta calculation


   Float_t   beta_help = beta_mean_tof ;  // pow(beta_mean_tof,1.0) gave best results

//-----------------------  Select helium   --------------------------

   Int_t  icount=0;

   for (jj=0; jj<24; jj++){
   a1 = adca[2*jj]*cos(theta);
   a2 = adca[2*jj+1]*cos(theta);

   xhelp = 100000.;
   if ((a1 < 3000) && (a2 < 3000)) xhelp = sqrt(a1*a2);  // geometric mean
// if geometric mean multiplied by beta_help  is inside helium limits, increase counter
   if ((beta_mean_tof>0.6) && (beta_mean_tof<1.1) && 
        ((beta_help*xhelp)>A_l[jj]) && ((beta_help*xhelp)<A_h[jj])) icount++ ;
                     }

  Int_t   iz=0;
//  if (icount > 3) iz=2;   // if more than three paddles see helium, then set Z=2
  if (icount > 4) iz=2;   

//----------------------  Z=2 fill histograms  -----------------------------

    if (iz==2) {

    heevent++;
    for (jj=0; jj<48; jj++)  adche[jj]->Fill(adca[jj]);
                           
                           }  // iz0==2


//=====================  end beta and helium part ===========================

   }  //  i < nevents
 

   float *X       = new float[48];
   float *means   = new float[48];
   float *entries = new float[48];
   int *entriestdc = new int[48];
   int *entriesadc = new int[48];
   
   const char *saveas = format;
   
   int i=0;
   
   gStyle->SetStatW(0.4);
   gStyle->SetStatH(0.4);
   gStyle->SetOptStat("nmri");
   gStyle->SetTitleH(0.10);
   gStyle->SetTitleW(0.96);
   
   TCanvas *SCanvas = new TCanvas("SCanvas","SCanvas", 1280, 1024);
   SCanvas->Divide(4,2);

   j = 0;
   while (j < 12){
     k = 0;
     z = 0;
     if (gROOT->IsBatch()) {
       SCanvas = new TCanvas("SCanvas","SCanvas", 1280, 1024);
       SCanvas->Divide(4,2);
     } else {
       if (j > 0) SCanvas->DrawClone();
     }
     
     
     while(k < 4){
       if (k > 1) z = 2;
       i = j*4 + k;
       X[i] = i;
       
       SCanvas->cd(k+3+z);
       htdc[i] = (TH1F*)hhtdc->At(i);  
       entriestdc[i] = (Int_t)htdc[i]->Integral();
       sst.str("");
       sst << "TDC - " << photoS[i].Data() << " (Nev < 4096 = " << entriestdc[i] << ")";
       htdc[i]->SetTitle(sst.str().c_str());
       htdc[i]->SetTitleSize(10);
       htdc[i]->SetAxisRange(690,1510);
       htdc[i]->DrawCopy();
       htdc[i]->ComputeIntegral();
       entries[i] = htdc[i]->Integral();
       
       SCanvas->cd(k+1+z);       
       hadc[i] = (TH1F*)hhadc->At(i);
       entriesadc[i] = (Int_t)hadc[i]->Integral();
       sst.str("");
       sst << "ADC - " << photoS[i].Data() << " (Nev < 4096 = " << entriesadc[i] << ")";
       hadc[i]->SetTitle(sst.str().c_str());
       hadc[i]->SetAxisRange(-10,710); 
       hadc[i]->DrawCopy();
       means[i]   = hadc[i]->GetMean();
       
       k++;
     }
     

     if ( !strcmp(saveas,"ps") ) {
       sst.str("");
       sst << outDirectory.Data() << filename.Data() << "TOFScan.ps(";
       SCanvas->Print(sst.str().c_str());
     } else {
       sst.str("");
       sst << outDirectory.Data() << filename.Data() << "TOFScan" << j+1 << "." << saveas;
       SCanvas->SaveAs(sst.str().c_str());
       
     }
     j++;
   }

   if (gROOT->IsBatch()) SCanvas->Close();
   
   /*
    * This Canvas will represent a summary of the performances for TOF TDC/ADC channels
    */
//   TCanvas *performanceCanvas = new TCanvas("performanceCanvas","performanceCanvas", 1280, 1024);
   TCanvas *performanceCanvas = new TCanvas("performanceCanvas","performanceCanvas", 1024, 1280);
   performanceCanvas->Divide(1,3);
   
   gStyle->SetTitleW(.9);

   performanceCanvas->cd(1);
   TGraph *adcMeans   = new TGraph(48, X, means);
   sst.str("");
   sst << "ADCMean"  << " - Data in " << base.Data() << " - Nevents in the run = " << nevents;
   adcMeans->SetTitle(sst.str().c_str());
   adcMeans->SetFillColor(35);
   adcMeans->GetXaxis()->SetTitle("Photomultipliers");
   adcMeans->GetXaxis()->CenterTitle();
   adcMeans->GetXaxis()->SetLimits(-0.5, 47.5);
   adcMeans->GetYaxis()->SetTitle("ADCMean");
   adcMeans->GetYaxis()->CenterTitle(); 
   adcMeans->Draw("AB");
   
   performanceCanvas->cd(2);
   TGraph *tdcEntries = new TGraph(48, X, entries);
   sst.str("");
   sst << "TDCEntries"  << " - Data in " << base.Data() << " - Nevents in the run = " << nevents;
   tdcEntries->SetTitle(sst.str().c_str());
   tdcEntries->SetFillColor(35);
   tdcEntries->GetXaxis()->SetTitle("Photomultipliers");
   tdcEntries->GetXaxis()->CenterTitle();
   tdcEntries->GetXaxis()->SetLimits(-0.5, 47.5);
   tdcEntries->GetYaxis()->SetTitle("TDCIntegral");
   tdcEntries->GetYaxis()->CenterTitle();       
   tdcEntries->Draw("AB");
 
//--------- new part PMT gain check  -----------------------------

   performanceCanvas->cd(3);

Float_t  xc[48],xmean1[48],xmeana[48];
Float_t xmean_arr[12][48];

//  xmean values from 2-3 april 2007

char date_info[]="Reference Data: apr-2007";

Float_t  xmean[48] = {
491.609,509.241,400.786,530.122,699.674,555.747,521.04,486.363,
470.173,227.752,611.038,455.889,553.601,520.54,403.527,382.099,
349.697,365.113,447.653,377.667,517.815,572.932,338.501,436.681,
485.696,450.491,395.375,329.631,751.258,626.681,385.561,578.476,
374.454,356.733,641.888,562.767,582.849,521.748,527.043,505.89,
489.828,628.408,532.924,506.511,482.872,532.236,554.554,498.849 };

// new 01-oct-2007
int channelmap[] = {0,7,3,6,2,8,1,5,3,7,3,6,1,7,2,10,
                    10,10,10,5,0,7,0,5,0,6,1,5,
                    2,8,3,8,2,6,1,8,
                    11,9,11,11,9,11,4,4,4,9,9,4};


int colormap[] = {46,2,29,4,5,6,7,8,9,11,28,34};
//int colormap[] = {417,400,632,617,603,600,434,419,591,625,403,424};


for (Int_t j=0; j<48; j++) xmeana[j]=0.;
for (Int_t j=0; j<24; j++) xmeana[2*j]=xmean[2*j];

for (Int_t i=0; i<12; i++) {
for (Int_t j=0; j<48; j++) {
xmean_arr[i][j]=0.;
}
}

for (Int_t j=0; j<48; j++) {
Int_t ichan = channelmap[j];
xmean_arr[ichan][j]=xmean[j];
}

// get results from ADC histogram
for (Int_t j=0; j<48; j++) {
xc[j]=j;
xmean1[j]=adche[j]->GetMean();
}


   gStyle->SetTitleW(.5);
   gStyle->SetTitleH(.05);

TH2F *hr = new TH2F("frame","2-Dim",2,-0.5,47.5,2,-300.,100.);
hr->SetStats(kFALSE);
hr->GetXaxis()->CenterTitle();
hr->GetXaxis()->SetTitle("Photomultipliers");
hr->GetYaxis()->CenterTitle();  
hr->GetYaxis()->SetTitle("Mean ADC Difference");
hr->SetTitle("Difference between Reference and Actual Values");
hr->Draw();

Int_t npoint=48;

for (Int_t j=0; j<12; j++) {
for (Int_t i=0; i<48; i++) xmeana[i] = 0.;
for (Int_t i=0; i<48; i++) {
if (xmean_arr[j][i] != 0)  xmeana[i] = xmean1[i] - xmean_arr[j][i];
                           }


TGraph *graph1 = new TGraph(npoint,xc,xmeana);
graph1->SetFillColor(colormap[j]);
graph1->GetXaxis()->SetLimits(-0.5, 47.5);
graph1->Draw("BP");
}

Float_t tp[10];
tp[0] = 15.5;
tp[1] = 27.5;
tp[2] = 31.5;
tp[3] = 35.5;
tp[4] = 41.5;

for (Int_t ii=0; ii<5; ii++) {
TLine *l1=new TLine(tp[ii],-300,tp[ii],100);
l1->SetLineColor(38);
l1->Draw("same");
                        }

for (Int_t j=0; j<12; j++) {
sprintf(htitle, "HV_%d",j);
TText *text1 = new TText(0+j*4,80,htitle);
text1->SetTextColor(colormap[j]);
//text1->SetTextSize(0.03);
text1->SetTextSize(0.05);
text1->Draw();
}


TText *text1 = new TText(0,-185,date_info);
text1->SetTextColor(kBlack);
text1->SetTextSize(0.023);
text1->Draw();


sprintf(htitle, "Helium Events: %d",heevent);
TText *text2 = new TText(20,-185,htitle);
text2->SetTextColor(kBlack);
text2->SetTextSize(0.023);
text2->Draw();


for (Int_t i=0; i<6; i++) {
for (Int_t j=0; j<8; j++) {
Int_t ihelp = i*8+j;
sprintf(htitle, "%d: %.0f/%.0f",(ihelp+1),xmean[ihelp],xmean1[ihelp]);
TText *text1 = new TText(0+j*6,-200-i*15,htitle);
text1->SetTextColor(kBlack);
text1->SetTextSize(0.023);
text1->Draw();
}
}

//-------- end new part -------------------------

  
   //------print the ps

   if ( !strcmp(saveas,"ps") ) {
     sst.str("");
     sst << outDirectory.Data() << filename.Data() << "TOFScan.ps)";
     performanceCanvas->Print(sst.str().c_str());
     
   } else {
     sst.str("");
     sst << outDirectory.Data() << filename.Data() << "TOFScan13." << saveas;
     performanceCanvas->SaveAs(sst.str().c_str());
   }
   if (gROOT->IsBatch()) {
     SCanvas->Close();
     performanceCanvas->Close();
   }


}

int main(int argc, char* argv[]){
  TString path;
  TString outDir        ="./";
  TString format        ="ps";
    
 if (argc < 2){
    printf("You have to insert at least the file to analyze \n");
    printf("Try '--help' for more information. \n");
    exit(1);
  }  

  if (!strcmp(argv[1], "--help")){
        printf( "Usage: TofScan FILE [OPTION] \n");
        printf( "\t --help                  Print this help and exit \n");      
        printf( "\t -outDir[path]           Path where to put the output [default ./] \n");
        printf( "\t -format[ps]             Format for output files [default 'ps'] \n");
        exit(1);
  }


  path=argv[1];

   for (int i = 2; i < argc; i++){

     if (!strcmp(argv[i], "-outDir")){
       if (++i >= argc){
         printf( "-outDir needs arguments. \n");
         printf( "Try '--help' for more information. \n");
         exit(1);
        } 
       else{
         outDir = argv[i];
         continue;
       }
     }
     
     
     if (!strcmp(argv[i], "-format")){
       if (++i >= argc){
         printf( "-format needs arguments. \n");
         printf( "Try '--help' for more information. \n");
         exit(1);
       } 
       else{
         format = argv[i];
         continue;
       }
     }
  }

   TofScan(argv[1], outDir, format);

}
