HarpoAnalyseGainStudy.cxx

Go to the documentation of this file.

// 
// File HarpoAnalyseGainStudy.cxx
//
#include "HarpoAnalyseGainStudy.h"
#include "HarpoConfig.h"
#include "HarpoDetSet.h"
#include "HarpoDebug.h"
#include "HarpoDccEvent.h"
#include "HarpoEvent.h"
#include "Pmm2Event.h"
#include "HarpoRecoEvent.h"
#include "HarpoSimEvent.h"

#include "TFile.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TGraph.h"
#include "TF1.h"
#include "TMath.h"
#include "TSystem.h"

#include <cstdlib>
#include <cstring>
#include <cassert>
#include <fstream>
#include <iostream>

ClassImp(HarpoAnalyseGainStudy);

void   HarpoAnalyseGainStudy::print()
 {
   unsigned int nd; // number of detectors
   HarpoEventHeader *hdr = fEvt->GetHeader();

   assert(hdr != NULL);
   hdr->print();
   
   for (nd = 0; nd < gkNDetectors; nd++) {
     //     if (fEvt->isdataExist(nd)) {
       HarpoDccMap *plane = fEvt->GetDccMap(nd);
       if (plane != NULL )plane->print();
   }
 }

void   HarpoAnalyseGainStudy::process()
{
  //   Bool_t drawEvent = kFALSE;
  nEvents++;
  
  if(nEvents%100 == 0)
    std::cout << "Event " << nEvents << std::endl;
  
  
  
  
  
  
  if (gHDetSet->isExist(SIMDET)){
    
    // Process RAW data
    Double_t qTotEvt = 0;
    for(UInt_t ndet=0;ndet<2;ndet++) {
      //std::cout << "NDET: " << ndet << " => " << gHDetSet->isExist(ndet) << std::endl;
      if (!gHDetSet->isExist(ndet)) continue;
      HarpoDccMap* m = fEvt->GetDccMap(ndet);
      if ( m == NULL ) continue;
      for(Int_t i=0;i<NADC;i++){ //Time bins
        for(Int_t j=0;j<NCHAN;j++){ // Channels
          Double_t q = m->GetData(j,i);// - 250.;
          if(q <= -1000) continue;
          qTotEvt += q;
        }
      }
    }
    
    
    
    HarpoSimEvent *simEvt = (HarpoSimEvent *)(fEvt->GetDetEvent(SIMDET));
    Int_t nElectronsTot = 0;
    Int_t nIonTr = simEvt->GetNtracks();
    for(Int_t i = 0; i<nIonTr; i++){
      TpcSimIonisationTrack* tr = simEvt->GetIonisationTrack(i);
      Int_t nPoints = tr->GetNpoints();
      for(Int_t j = 0; j<nPoints; j++){
        TpcSimIonisationPoint* point = tr->GetPoint(j);
        nElectronsTot += point->GetNe();
      }
    }

    if(gHarpoDebug>1)
      Info("process","%g / %d = %g",qTotEvt,nElectronsTot,qTotEvt/nElectronsTot);
    hQoverNe->Fill(qTotEvt/nElectronsTot);

  }
  



   // if (gHDetSet->isExist(PMM2)){
   //   Pmm2Event *anaEvt = static_cast<Pmm2Event *>(fEvt->GetDetEvent(PMM2));
   //   Int_t esize = anaEvt->GetHeader()->eventSize;
   //   Pmm2MesVect * pm = anaEvt->GetMesurements();
   //   Int_t hit[16];
   //   for(Int_t ch = 0; ch<16; ch++) hit[ch] = 0;
   //   for(Int_t imes=0;imes<esize;imes++) {
   //     int ch    = pm->at(imes).getChNum(); // channel
   //     int q     = pm->at(imes).getCharge(); // charge
   //     hit[ch] = 1;
   //   }
   //   if(!(hit[3] && hit[4] && hit[5] && hit[6]))
   //     return;
   // }

   // Process Reconstructed data (clusters, tracks, matched tracks)
   HarpoRecoEvent* reco = fEvt->GetRecoEvent();
   TClonesArray* clArray = reco->GetClustersArray();
   TClonesArray* trArray = reco->GetTracksArray();
   // Int_t NtrX=reco->GetNtracksXEvt();
   // Int_t NtrY=reco->GetNtracksYEvt();


   // if(reco->GetNclXEvt()<20)
   //   return;
   // if(reco->GetQtotXEvt()/reco->GetNclXEvt()<2000)
   //   return;

   // if(NtrX!=1) return;
   // if(NtrY!=1) return;
   const Int_t nTr = trArray->GetEntries();
   if(gHarpoDebug>1)
     Info("process","nTr = %d",nTr);
   if(nTr != 2) return;
   Int_t nCl = clArray->GetEntries();

   Double_t angle[2];
   //   Double_t angleCut = 100;//0.5;
   for(Int_t itr = 0; itr<nTr; itr++){
     HarpoRecoTracks* track = (HarpoRecoTracks*)trArray->At(itr);
     Int_t plane = track->GetPlane();
     angle[plane] = track->GetAngleTrack();
     hAngleTrack[plane]->Fill(angle[plane]);
   }

   if(gHarpoDebug>1)
     Info("process","angles = %g %g",angle[0],angle[1]);
   if(angle[0]==0) return;
   if(angle[1]==0) return;

   // Double_t corr = 1./TMath::Cos(angle[0])/TMath::Cos(angle[1]);
   // Double_t corr2 = TMath::Abs(1./TMath::Sin(angle[0])/TMath::Sin(angle[1]));
   Double_t alpha = TMath::Abs(TMath::Cos(angle[0])*TMath::Cos(angle[1]));
   Double_t alpha2 = TMath::Abs(TMath::Sin(angle[0])*TMath::Sin(angle[1]));
   hAlpha->Fill(alpha);

   Double_t tmin[2], tmax[2];
   tmin[0] = 512;
   tmin[1] = 512;
   tmax[0] = 0;
   tmax[1] = 0;
   const Int_t maxN = 20;
   Double_t Qx[511][maxN], Qy[511][maxN];
   // Double_t Wx[511][maxN], Wy[511][maxN];
   Double_t X[511][maxN], Y[511][maxN];
   Int_t Nx[511], Ny[511];
   Double_t Qxtot[511], Qytot[511];
   for(Int_t i = 0; i<511; i++) {
     for(Int_t j = 0; j<maxN; j++) {
       Qx[i][j] = 0;
       Qy[i][j] = 0;
       // Wx[i][j] = 0;
       // Wy[i][j] = 0;
       X[i][j] = 0;
       Y[i][j] = 0;
     }
     Nx[i] = 0;
     Ny[i] = 0;
     Qxtot[i] = 0;
     Qytot[i] = 0;
   }
   Double_t QtotX = 0, QtotY = 0;
   for(Int_t icl = 0; icl<nCl; icl++){
     HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
     //     Int_t size = cluster->GetWidth();
     //     if(size<4) continue;
     Int_t plane = cluster->GetPlane();
     Int_t index = cluster->GetIndex();
     //     Int_t width = cluster->GetWidth();
     Double_t mean = cluster->GetMean();
     Double_t q = cluster->GetQ();

     if(cluster->GetType()==CCLUSTER) {
       if(!(TMath::Abs(fmod(mean+16.5,33)-11.5)<5))
         hQclVsAlpha2[plane]->Fill(alpha,q*alpha2);
       if(alpha<0.05) continue;
       if(alpha>0.4) continue;
       if(!(TMath::Abs(fmod(index+16.5,33)-11.5)<5))
         hQclVsT2[plane]->Fill(mean,q*alpha2);
       hQclVsCh2[plane]->Fill(index,q*alpha2);
     }
     
     if(cluster->GetType()!=TCLUSTER) continue;
     if(index<tmin[plane]) tmin[plane] = index;
     if(index>tmax[plane]) tmax[plane] = index;
     hMean[plane]->Fill(mean);
     if(plane==XPLANE){
       QtotX += q;
       Qxtot[index] += q;
       if(Nx[index]<maxN){
         Qx[index][Nx[index]] = q;
         //      Wx[index][Nx[index]] = width;
         X[index][Nx[index]] = mean;
         Nx[index]++;
       }
     }
     if(plane==YPLANE){
       QtotY += q;
       Qytot[index] += q;
       if(Ny[index]<maxN){
         Qy[index][Ny[index]] = q;
         //      Wy[index][Ny[index]] = width;
         Y[index][Ny[index]] = mean;
         Ny[index]++;
       }
     }
     if(!(TMath::Abs(fmod(mean+16.5,33)-11.5)<5))
       hMeanCut[plane]->Fill(mean);
     //     hQclVsAngle[plane]->Fill(index,alpha,q*corr);
     // hQcl2D[plane]->Fill(index,q);
     // hQcl[plane]->Fill(index,q);
     if(!(TMath::Abs(fmod(mean+16.5,33)-11.5)<5))
       hQclVsAlpha[plane]->Fill(alpha,q*alpha);
     if(alpha<0.5) continue;
     if(alpha>0.95) continue;
     if(!(TMath::Abs(fmod(mean+16.5,33)-11.5)<5))
       hQclVsT[plane]->Fill(index,q*alpha);
     hQclVsCh[plane]->Fill(mean,q*alpha);
     
   }

   for(Int_t i = 0; i<511; i++){
     for(Int_t j = 0; j<Nx[i]; j++){
       if(Qx[i][j] == 0) continue;
       hRatioClVsX->Fill(X[i][j],Qx[i][j]/Qytot[i]);
       //       Qxtot += Qx[i][j];
       for(Int_t k = 0; k<Ny[i]; k++){
         if(Qy[i][k] == 0) continue;
         //fNtuple->Fill(i,Qx[i][j],Qy[i][k],Wx[i][j],Wy[i][k],Nx[i],Ny[i],X[i][j],Y[i][k],alpha,QtotX,QtotY);
         hRatioClVsT->Fill(i,Qx[i][k]/Qy[i][j]);
         if(j!=0) continue;
         hRatioClVsY->Fill(Y[i][k],Qy[i][k]/Qxtot[i]);
       }
     }
     if(Qxtot[i] == 0) continue;
     if(Qytot[i] == 0) continue;
     hQVsT[0]->Fill(i,Qxtot[i]*alpha);
     hQVsT[1]->Fill(i,Qytot[i]*alpha);
     if(i>tmin[0]+10 && i<tmax[0]-10)
       hQVsTcut[0]->Fill(i,Qxtot[i]*alpha);
     if(i>tmin[1]+10 && i<tmax[1]-10)
       hQVsTcut[1]->Fill(i,Qytot[i]*alpha);
     
     hQVsAlpha[0]->Fill(alpha,Qxtot[i]*alpha);
     hQVsAlpha[1]->Fill(alpha,Qytot[i]*alpha);
     //     fNtuple->Fill(i,Qxtot[i],Qytot[i],Wx[i][0],Wy[i][0],Nx[i],Ny[i],X[i][0],Y[i][0],alpha,QtotX,QtotY);
     hRatioVsAngle->Fill(i,alpha,Qxtot[i]/Qytot[i]);
   }
   //std::cout << "PROCESSING RAW DATA" << std::endl;








   //   Info("process","Check Sim");
   if (gHDetSet->isExist(SIMDET)){
     //     Info("process","Sim Data exists");
     HarpoSimEvent *simEvt = (HarpoSimEvent *)(fEvt->GetDetEvent(SIMDET));

     TpcSimMCEvent* mcEvent = (TpcSimMCEvent*)simEvt->GetMCEvent();
     TpcSimMCTrack* mcTrack = mcEvent->GetTrack(0);
     Double_t pz = mcTrack->GetPz();
     //     Info("process","%p",simEvt);

     Int_t nIonTr = simEvt->GetNtracks();
     //     Info("process","%p   %d",simEvt, nIonTr);
     TH1F* hTmp = new TH1F("hTmp","",512,0,512);
     for(Int_t i = 0; i<nIonTr; i++){
       //       Info("process","%d", i);
       TpcSimIonisationTrack* tr = simEvt->GetIonisationTrack(i);
       Int_t nPoints = tr->GetNpoints();
       for(Int_t j = 0; j<nPoints; j++){
         TpcSimIonisationPoint* point = tr->GetPoint(j);
         hTmp->Fill(point->GetZ(),point->GetNe()*pz);
         hNeVsLambda->Fill(point->GetZ()/pz,point->GetNe());
       }
     }
     for(Int_t i = 0; i<512; i++)
       hNeVsZSim1->Fill(hTmp->GetXaxis()->GetBinCenter(i+1),400*hTmp->GetBinContent(i+1)/1.);
     hTmp->Rebin(2);
     for(Int_t i = 0; i<256; i++)
       hNeVsZSim2->Fill(hTmp->GetXaxis()->GetBinCenter(i+1),400*hTmp->GetBinContent(i+1)/2.);
     hTmp->Rebin(2);
     for(Int_t i = 0; i<128; i++)
       hNeVsZSim3->Fill(hTmp->GetXaxis()->GetBinCenter(i+1),400*hTmp->GetBinContent(i+1)/4.);
     hTmp->Rebin(2);
     for(Int_t i = 0; i<64; i++)
       hNeVsZSim4->Fill(hTmp->GetXaxis()->GetBinCenter(i+1),400*hTmp->GetBinContent(i+1)/8.);
     
     hTmp->Delete();
   }

   // Process RAW data
   //   Double_t sigQ[NADC][2];
   //   for(Int_t ndet=0;ndet<gkNDetectors;ndet++) {
   // for(UInt_t ndet=0;ndet<2;ndet++) {
   //   //std::cout << "NDET: " << ndet << " => " << gHDetSet->isExist(ndet) << std::endl;
   //   if (!gHDetSet->isExist(ndet)) continue;
   //   ULong_t timestamp = fEvt->GetDetEvent(ndet)->GetTimeStamp();
   //   HarpoDccMap* m = fEvt->GetDccMap(ndet);
   //   if ( m == NULL ) continue;
   //   Double_t qTotEvt = 0;
   //   Int_t nCl = 0;
   //   for(Int_t i=0;i<NADC;i++){ //Time bins
   //     if(tmin[ndet]>100 && i<tmin[ndet]+20) continue;
   //     if(tmax[ndet]<400 && i>tmax[ndet]-20) continue;
   //     Double_t qTot = 0;
   //     Double_t sigmaQ = 0;
   //     for(Int_t j=0;j<NCHAN;j++){ // Channels
   //    Double_t q = m->GetData(j,i);// - 250.;
   //    if(q<2) continue;
   //    qTot += q;
   //    sigmaQ += q*q;
   //    hPosVsAngle[ndet]->Fill(i,j,alpha);
   //    hMap[ndet]->Fill(i,j,q);
   //     }
   //     if(qTot>0){
   //    nCl++;
   //    hQraw[ndet]->Fill(i,qTot);
   //    //      hQrawVsAngle[ndet]->Fill(i,alpha,qTot*corr);
   //    qTot *= alpha;
   //    hQrawCut[ndet]->Fill(i,qTot);
   //     }
   //     qTotEvt += qTot;
   //     //       if(angle[ndet]>angleCut && angle[ndet]<TMath::Pi()-angleCut) continue;
   //     // if(alpha>angleCut && alpha<TMath::Pi()-angleCut) continue;
   //     // if(qTot>0)
   //     //     hQrawCut[ndet]->Fill(i,qTot);
   //     // // sigmaQ /= NCHAN;
   //     // // qTot /= NCHAN;
   //     // sigmaQ -= qTot*qTot;
   //     // sigQ[i][ndet] = sigmaQ;
   //     // if(sigmaQ) 
   //     //     hSigmaQ[ndet]->Fill(sigmaQ);
   //   }
   //   if(nCl)
   //     qTotEvt /= nCl;
   //   hQ[ndet]->Fill(qTotEvt);
   // }


   // for(Int_t i=0;i<NADC;i++)
   //   hSigmaQXY->Fill(sigQ[i][0],sigQ[i][1]);



   //   return;



   // hTminX->Fill(tmin[0]);
   // hTminY->Fill(tmin[1]);
   // hTmaxX->Fill(tmax[0]);
   // hTmaxY->Fill(tmax[1]);
   // hTminmaxX->Fill(tmin[0],tmax[0]);
   // hTminmaxY->Fill(tmin[1],tmax[1]);
   // hTdriftX->Fill(tmax[0]-tmin[0]);
   // hTdriftY->Fill(tmax[1]-tmin[1]);



   if(nEvents%100 == 0)
     std::cout << "Event " << nEvents << " done" << std::endl;
 }

 void   HarpoAnalyseGainStudy::Init()
 {

  hAlpha = new TH1F("hAlpha",";track angle",100,0,1);

  // Initialise histograms here
  hAngleTrack[0] = new TH1F("hAngleTrackX",";track angle",100,0,2*TMath::Pi());
  hAngleTrack[1] = new TH1F("hAngleTrackY",";track angle",100,0,2*TMath::Pi());

   hMean[0] = new TH1F("hMeanX",";X",300,0,300);
   hMean[1] = new TH1F("hMeanY",";Y",300,0,300);
   hMeanCut[0] = new TH1F("hMeanCutX",";X",300,0,300);
   hMeanCut[1] = new TH1F("hMeanCutY",";Y",300,0,300);

   hRatioVsAngle = new TH3F("hRatioVsAngle",";t;#alpha;Q_{x}/Q_{y}",512,0,512,100,0,1,500,0,5);

   fNeventCut[0] = 0;
   fNeventCut[1] = 0;

   hQoverNe = new TH1F("hQoverNe","",1000,0,1000);

 
   fNtuple = new TNtuple("fNtuple","ntuple","t:Qx:Qy:Wx:Wy:Nx:Ny:X:Y:alpha:QtotX:QtotY");
   fNtuple->SetDirectory(0);

  const Int_t nbinsQ = 500;
  Double_t xminQ = 10;
  Double_t xmaxQ = 1e5;
  Double_t logxminQ = TMath::Log(xminQ);
  Double_t logxmaxQ = TMath::Log(xmaxQ);
  Double_t binwidthQ = (logxmaxQ-logxminQ)/nbinsQ;
  Double_t xbinsQ[nbinsQ+1];
  xbinsQ[0] = xminQ;
  for (Int_t i=1;i<=nbinsQ;i++)
    xbinsQ[i] = TMath::Exp(logxminQ+i*binwidthQ);
  hQ[0] = new TH1F("hQX",";Q_{cl}",nbinsQ,xbinsQ);
  hQ[1] = new TH1F("hQY",";Q_{cl}",nbinsQ,xbinsQ);
  hQclVsT[0] = new TH2F("hQclVsTX",";t;Q_{cl}",512,0,512,nbinsQ,xbinsQ);
  hQclVsT[1] = new TH2F("hQclVsTY",";t;Q_{cl}",512,0,512,nbinsQ,xbinsQ);
  hQclVsCh[0] = new TH2F("hQclVsChX",";channel;Q_{cl}",304,0,304,nbinsQ,xbinsQ);
  hQclVsCh[1] = new TH2F("hQclVsChY",";channel;Q_{cl}",304,0,304,nbinsQ,xbinsQ);
  hQclVsAlpha[0] = new TH2F("hQclVsAlphaX",";#alpha;Q",100,0,1,nbinsQ,xbinsQ);
  hQclVsAlpha[1] = new TH2F("hQclVsAlphaY",";#alpha;Q",100,0,1,nbinsQ,xbinsQ);
  hQclVsT2[0] = new TH2F("hQclVsT2X",";t;Q_{cl}",512,0,512,nbinsQ,xbinsQ);
  hQclVsT2[1] = new TH2F("hQclVsT2Y",";t;Q_{cl}",512,0,512,nbinsQ,xbinsQ);
  hQclVsCh2[0] = new TH2F("hQclVsCh2X",";channel;Q_{cl}",304,0,304,nbinsQ,xbinsQ);
  hQclVsCh2[1] = new TH2F("hQclVsCh2Y",";channel;Q_{cl}",304,0,304,nbinsQ,xbinsQ);
  hQclVsAlpha2[0] = new TH2F("hQclVsAlpha2X",";#alpha;Q",100,0,1,nbinsQ,xbinsQ);
  hQclVsAlpha2[1] = new TH2F("hQclVsAlpha2Y",";#alpha;Q",100,0,1,nbinsQ,xbinsQ);

  hNeVsZSim1 = new TH2F("hNeVsZSim1",";Z [mm];N_{electrons}",512,0,512,nbinsQ,xbinsQ);
  hNeVsZSim2 = new TH2F("hNeVsZSim2",";Z [mm];N_{electrons}",256,0,512,nbinsQ,xbinsQ);
  hNeVsZSim3 = new TH2F("hNeVsZSim3",";Z [mm];N_{electrons}",128,0,512,nbinsQ,xbinsQ);
  hNeVsZSim4 = new TH2F("hNeVsZSim4",";Z [mm];N_{electrons}",64,0,512,nbinsQ,xbinsQ);
  hNeVsLambda = new TH2F("hNeVsLambda",";Z [mm];N_{electrons}",1000,0,1000,nbinsQ,xbinsQ);

  hQVsT[0] = new TH2F("hQVsTX",";channel;Q_{cl}",512,0,512,nbinsQ,xbinsQ);
  hQVsT[1] = new TH2F("hQVsTY",";channel;Q_{cl}",512,0,512,nbinsQ,xbinsQ);
  hQVsTcut[0] = new TH2F("hQVsTcutX",";channel;Q_{cl}",512,0,512,nbinsQ,xbinsQ);
  hQVsTcut[1] = new TH2F("hQVsTcutY",";channel;Q_{cl}",512,0,512,nbinsQ,xbinsQ);
  hQVsAlpha[0] = new TH2F("hQVsAlphaX",";#alpha;Q",100,0,1,nbinsQ,xbinsQ);
  hQVsAlpha[1] = new TH2F("hQVsAlphaY",";#alpha;Q",100,0,1,nbinsQ,xbinsQ);

  hRatioClVsX = new TH2F("hRatioClVsX",";X;Q",304,0,304,100,0,5);
  hRatioClVsY = new TH2F("hRatioClVsY",";Y;Q",304,0,304,100,0,5);
  hRatioClVsT = new TH2F("hRatioClVsT",";t;Q",512,0,512,100,0,5);


  // const Int_t nbins2 = 500;
  // Double_t xmin2 = 1;
  // Double_t xmax2 = 1e6;
  // Double_t logxmin2 = TMath::Log(xmin2);
  // Double_t logxmax2 = TMath::Log(xmax2);
  // Double_t binwidth2 = (logxmax2-logxmin2)/nbins2;
  // Double_t xbins2[nbins2+1];
  // xbins2[0] = xmin2;
  // for (Int_t i=1;i<=nbins2;i++)
  //   xbins2[i] = TMath::Exp(logxmin2+i*binwidth2);
  // hSigmaQ[0] = new TH1F("hSigmaQX",";#sigma_{Q}",nbins2,xbins2);
  // hSigmaQ[1] = new TH1F("hSigmaQY",";#sigma_{Q}",nbins2,xbins2);
  // hSigmaQXY = new TH2F("hSigmaQXY",";#sigma_{Q}",nbins2,xbins2,nbins2,xbins2);

  // const Int_t nbinsPrf = 100;
  // Double_t xminPrf = 1e-3;
  // Double_t xmaxPrf = 2;
  // Double_t logxminPrf = TMath::Log(xminPrf);
  // Double_t logxmaxPrf = TMath::Log(xmaxPrf);
  // Double_t binwidthPrf = (logxmaxPrf-logxminPrf)/nbinsPrf;
  // Double_t xbinsPrf[nbinsPrf+1];
  // xbinsPrf[0] = xminPrf;
  // for (Int_t i=1;i<=nbinsPrf;i++)
  //   xbinsPrf[i] = xminPrf + TMath::Exp(logxminPrf+i*binwidthPrf);
  // hPrf[0] = new TH2F("hPrf_X",";",1000,-10,10,nbinsPrf,xbinsPrf);
  // hPrf[1] = new TH2F("hPrf_Y",";",1000,-10,10,nbinsPrf,xbinsPrf);
  // hPrf2[0] = new TH2F("hPrf2_X",";",1000,-10,10,nbinsPrf,xbinsPrf);
  // hPrf2[1] = new TH2F("hPrf2_Y",";",1000,-10,10,nbinsPrf,xbinsPrf);
  // hPrfT[0] = new TH2F("hPrfT_X",";",200,-20,30,nbinsPrf,xbinsPrf);
  // hPrfT[1] = new TH2F("hPrfT_Y",";",200,-20,30,nbinsPrf,xbinsPrf);
  // hPrfTcut[0] = new TH2F("hPrfTcut_X",";",200,-20,30,nbinsPrf,xbinsPrf);
  // hPrfTcut[1] = new TH2F("hPrfTcut_Y",";",200,-20,30,nbinsPrf,xbinsPrf);



   // gTimeStampVsEvent[0] = new TGraph();
   // gTimeStampVsEvent[1] = new TGraph();

 }

void   HarpoAnalyseGainStudy::Save(char * /* mode */)
 {
   OpenHistFile("gain");

   // Save histograms here
   hAlpha->Write();
   hAngleTrack[0]->Write();
   hAngleTrack[1]->Write();
   hRatioVsAngle->Write();
   hQ[0]->Write();
   hQ[1]->Write();

   hMean[0]->Write();
   hMean[1]->Write();
   hMeanCut[0]->Write();
   hMeanCut[1]->Write();

  hQVsT[0]->Write();
  hQVsT[1]->Write();
  hQVsTcut[0]->Write();
  hQVsTcut[1]->Write();
  hQVsAlpha[0]->Write();
  hQVsAlpha[1]->Write();
  hQclVsT[0]->Write();
  hQclVsT[1]->Write();
  hQclVsCh[0]->Write();
  hQclVsCh[1]->Write();
  hQclVsAlpha[0]->Write();
  hQclVsAlpha[1]->Write();
  hQclVsT2[0]->Write();
  hQclVsT2[1]->Write();
  hQclVsCh2[0]->Write();
  hQclVsCh2[1]->Write();
  hQclVsAlpha2[0]->Write();
  hQclVsAlpha2[1]->Write();
  hRatioClVsX->Write();
  hRatioClVsY->Write();
  hRatioClVsT->Write();

  hQoverNe->Write();
  hNeVsZSim1->Write();
  hNeVsZSim2->Write();
  hNeVsZSim3->Write();
  hNeVsZSim4->Write();
  hNeVsLambda->Write();
 
 }



Double_t HarpoAnalyseGainStudy::TruncSigma(TArrayD* vect, Double_t &truncS, Double_t &truncM, Double_t tl, Double_t th)
{
  //
  // Calculates the truncated mean mean of the non zero value in vect
  // The truncation is done on the 100*tl% lowest and 100*(1-th) highest value
  //
  //  debug(2,"truncated mean");


  Int_t size = vect->GetSize();
  Double_t truncMean = 0;
  Double_t truncSigma = 0;
  Int_t* index = new Int_t[size];
  TMath::Sort(size,vect->GetArray(),index,kFALSE);
  Int_t t = 0, tLow, tHigh;


  while(vect->At(index[t])<10&&t<size-1) t++;
  tLow  = TMath::FloorNint(t + (size - t)*tl);
  tHigh = TMath::FloorNint(t + (size - t)*th);

  for(Int_t tind = tLow; tind<tHigh; tind++){
    truncMean += vect->At(index[tind]);
    truncSigma += vect->At(index[tind])*vect->At(index[tind]);
  }

  //  if(tHigh-tLow == 0) return 0;
  if(tHigh-tLow < 15) return 0;

  truncMean /= (tHigh-tLow);
  truncSigma /= (tHigh-tLow);
  truncSigma -= truncMean*truncMean;

  truncM = truncMean;
  truncS = TMath::Sqrt(truncSigma);

  return truncSigma;
}