HarpoAnalyseMonitorNew.cxx

Go to the documentation of this file.

// 
// File HarpoAnalyseMonitorNew.cxx
//
#include "HarpoAnalyseMonitorNew.h"
#include "HarpoConfig.h"
#include "HarpoDetSet.h"
//#include "HarpoDetSet.h"
#include "HarpoDccEvent.h"
#include "Pmm2Event.h"
#include "Pmm2Header.h"
#include "Pmm2Status.h"
#include "HarpoEvent.h"

#include "TROOT.h"
#include "TStyle.h"
#include "TApplication.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TGraph.h"
#include "TMath.h"
#include "TBox.h"
#include "TMapFile.h"
#include "TLinearFitter.h"
#include "TRandom.h"
#include "TSystem.h"

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

ClassImp(HarpoAnalyseMonitorNew)



enum HarpoMonitorEventType {
  E_CREATED,
  E_TRAVERSING,
  E_EMPTY,
  E_SATURATED,
  E_OTHER
};


void   HarpoAnalyseMonitorNew::print()
{
  UInt_t 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 HarpoAnalyseMonitorNew::process()
{
  nEvents++;
  if(nEvents%100 == 0)
    //    std::cout << " Processing Event " << nEvents<< std::endl;
    Info("process","Event %ld",nEvents);
  if(nEvents%100 == 0)
    Save();
  hTev[0]->Reset();
  hTev[1]->Reset();
  hTevZs[0]->Reset();
  hTevZs[1]->Reset();
  hQev[0]->Reset();
  hQev[1]->Reset();
  hQevZs[0]->Reset();
  hQevZs[1]->Reset();
  //hInfo->Reset();
  hPmm2->Reset();
  
  Double_t qtot[2] = {-1,-1};
  //  Double_t qtb[NADC][2];
  // Int_t npix = 0, nch = 0, ntb = 0;
  //  for(Int_t ndet=0;ndet<gkNDetectors;ndet++) {
  // Int_t qmin = 10000, qmax = 300000;
  // Double_t tmin = 150, lTpc = 200;
  Double_t timestamp = -1;
  //const Int_t TRIG_TT = 0;

  Bool_t isTT = false;
  
  hInfo->SetBinContent(1,gHConfig->GetRunNo());
  hInfo->SetBinContent(2,fEvt->GetHeader()->GetEvtNo());
  hInfo->SetBinContent(6,gHConfig->GetSkipEvents()+1);
  hInfo->SetBinContent(7,gHConfig->GetFromEvent());
  Int_t triggerType = -1;
  if(gHDetSet->isExist(PMM2)) {
    hInfo->SetBinContent(5,fEvt->GetTimeStamp(PMM2));
    timestamp = fEvt->GetTimeStamp(PMM2)*10;



    if(nEvents==1){
      
      Pmm2Header* pmm2Hdr = (Pmm2Header*)fRunHeader->GetHeader(PMM2);
      fTriggerTimeTot = 0;
      for(Int_t i = 0; i<kNtriggerTypes; i++){
        fTriggerTime[i] = 0;
        fTriggerDownscaling[i] = 1;
        if(pmm2Hdr == NULL) continue;
        Pmm2Status* pmm2Stat = pmm2Hdr->getStatus();
        if(pmm2Stat == NULL) continue;
        unsigned int triggerinfo = pmm2Stat->getTrigger(i); // 32bits: 
        // time (16bits * XXXms), 
        // downscaling (16bits)
        fTriggerTime[i] = triggerinfo/65536;
        fTriggerDownscaling[i] = triggerinfo%65536;
      }
      for(Int_t i = 0; i<kNtriggerTypes; i++)
        fTriggerTimeTot += fTriggerTime[i];
      for(Int_t i = 0; i<kNtriggerTypes; i++){
        Double_t trigCorr  = 1.;
        if(fTriggerTime[i])
          trigCorr /= fTriggerTime[i]/fTriggerTimeTot;
        if(fTriggerDownscaling[i])
          trigCorr *= fTriggerDownscaling[i];
        hTriggerCorrection->SetBinContent(i+2,trigCorr);
        Info("Init","Trigger param %d %d => %g",fTriggerTime[i],fTriggerDownscaling[i],trigCorr);
      }
      hTriggerCorrection->SetBinContent(1,1);
    }










    Pmm2Event *anaEvt = static_cast<Pmm2Event *>(fEvt->GetDetEvent(PMM2));
    if(gHarpoDebug>0)
      std::cout << "pmm2Event: " << anaEvt << std::endl;    

    triggerType = anaEvt->GetTriggerType();
    //Info("process","** trigger type: %d",triggerType);

    int imes;
    int esize = anaEvt->GetHeader()->eventSize;
    //int ievnum = anaEvt->GetHeader()->eventNumb;
    //Info("process","ievnum = %d, esize = %d",ievnum,esize);
    //Info("process","esize = %d",esize);
    hPmm2Size->Fill(esize);
    Pmm2MesVect * pm = anaEvt->GetMesurements();
    int hit[16];
    for(int ch = 0; ch<16; ch++) hit[ch] = 0;
    for(imes=0;imes<esize;imes++) {
      int ch    = pm->at(imes).getChNum(); // channel
      int q     = pm->at(imes).getCharge(); // charge
      //Info("process","ch = %d, q = %d",ch,q);
      hPmm2->SetBinContent(ch+1,q);
      hPmm2Spectrum2D->Fill(ch,q);
      hit[ch] = 1;
      if(gHarpoDebug>1)
        std::cout << ch << "  ";
    }
    if(gHarpoDebug>1)
      std::cout << std::endl;

// #if 0
//     // FIXME: fake trigger types
//   triggerType = 7+gRandom->Rndm()*6;

//   int nScint = (hit[1] && hit[2])
//     + (hit[3] && hit[4])
//     + (hit[5] && hit[6])
//     + (hit[7] && hit[8])
//     + (hit[9] && hit[10])
//     + (hit[11] && hit[12]);
//   if(nScint>1) triggerType = 14;

    if(hit[3] && hit[4] && hit[5] && hit[6]){
      if(gHarpoDebug>1) 
        Info("process","TT");
      isTT = true;
      hInfo->SetBinContent(9,hInfo->GetBinContent(9)+1);
    }
//   if(((hit[1] && hit[2]) || 
//       (hit[3] && hit[4]) || 
//       (hit[7] && hit[8]) || 
//       (hit[9] && hit[10]) || 
//       (hit[11] && hit[12])) &&
//      !(hit[5] && hit[6]))
//     triggerType = 5;

// #endif


    for(int ch = 0; ch<16; ch++){
      int ch2 = ch%2 ? ch+1 : ch-1;
        if(hit[ch2]){
          if(!hit[ch])
            hPmm2Missed->Fill(ch);
          hPmm2Triggered->Fill(ch);
      }
    }
    
    for(int ch1 = 0; ch1<16; ch1++){
      for(int ch2 = 0; ch2<16; ch2++){
        //       if(hit[ch1] && hit[ch2] && ch1 != ch2){
        if(hit[ch1] && hit[ch2]){
          hPmm2Coincidence->Fill(ch1,ch2);
        }
      }
    }
  }else{
    for(Int_t ch = 0; ch<16; ch++) hPmm2->SetBinContent(ch+1,-1);
  }
  


  if(gHarpoDebug>0)
    Info("process","isTT = %d",isTT);


  HarpoMonitorEventType eventType = E_OTHER;
  for(Int_t ndet=0;ndet<2;ndet++) {
    if (gHDetSet->isExist(ndet)) {
      hInfo->SetBinContent(3+ndet,fEvt->GetTimeStamp(ndet));
      timestamp = fEvt->GetTimeStamp(ndet);
      HarpoDccMap * m = fEvt->GetDccMap(ndet);
      //      std::cout << "Event " << nEvents << ", detector " << ndet << "  " << m << std::endl;
      if ( m != NULL ) {

        hZs[ndet]->Reset();
        h[ndet]->Reset();

        if(nEvents<10){
          
          Int_t Qnoise = 0;
          for(Int_t i=0;i<NADC;i++){ //Time bins
            for(Int_t j=NCHAN;j<NALL;j++){ // Noise Channels
              Double_t q = m->GetData(j,i);
              Qnoise += q;
            }
          }
          fIsZs = Qnoise<100*NADC;
        }
        h[ndet]->SetTitle(Form("X direction, evt: %li",nEvents));

        Double_t qTot = 0, qTotZs = 0, Tmin = -10, Xmin = -500, Tmax = -10;
        for(Int_t i=0;i<NADC;i++){ //Time bins
          //      qtb[i][ndet] = 0;
          Double_t qTotT = 0, qTotTZs = 0, meanCh = 0;
          Int_t nCh = 0;
          for(Int_t j=0;j<NALL;j++){ // Channels
            Double_t q = m->GetData(j,i);
            if(q<0) continue;
            h[ndet]->SetBinContent(i+1,j+1,q);
            qTotT += q;
            nCh++;
            meanCh += q*j;
            qTot += q;
            hCh[ndet]->Fill(j+0.5,q);
            hT[ndet]->Fill(i+0.5,q);
            hTev[ndet]->Fill(i+0.5,q);
            hQ[ndet]->Fill(q);
            hQev[ndet]->Fill(q);
            //}         
          }
          if(gHarpoDebug>0)
            Info("process","TT %d %f",i,qTotT);
          if(qTotT>0){
            hQtotT[ndet]->Fill(qTotT);
            hQtotTev[ndet]->Fill(qTotT);
            //      if(triggerType == TRIG_TT){
            if(isTT){
              if(gHarpoDebug>1){
                Info("process","TT %d %f",i,qTotT);
                std::cout << i << " " << qTotT << std::endl;
              }
              hQvsT_TT[ndet]->Fill(i,qTotT);
              hQvsT_TT2D[ndet]->Fill(i,qTotT);
            }
            meanCh /= qTotT;
          }
          
          // if(!fIsZs){

          //   Int_t Qnoise = 0;
          //   for(Int_t ii=0;i<NADC;i++){ //Time bins
          //     for(Int_t j=NCHAN;j<NALL;j++){ // Noise Channels
          //    Double_t q = m->GetData(j,ii);
          //   if(q<0) continue;
          //    Qnoise += q;
          //     }
          //   }
          //   //           Bool_t isZs = Qnoise<100*NADC;

          //   TLinearFitter* fFitCh = new TLinearFitter(1,"pol2","D");
          //   //           Double_t qTotT = 0;
          //   TArrayD* valuesTB = new TArrayD(NALL);
          //   for(Int_t j=0;j<NCHAN;j++){ // Channels
          //     Double_t q = m->GetData(j,i);
          //   if(q<0) continue;
          //     valuesTB->AddAt(q,j);
          //   }
          //   Double_t noiseTB = TruncMean(valuesTB,0.16,0.84);
          //   //if(i>490) Info("process","t = %d, noiseTB = %g",i,noiseTB);
          //   TGraph* gMinusCh = new TGraph();
          //   for(Int_t j=0;j<NCHAN;j++){ // Channels
          //     Double_t q = m->GetData(j,i);
          //     q -= noiseTB;
          //     m->SetData(j,i,q);
          //     gMinusCh->SetPoint(gMinusCh->GetN(),j,q);
          //   }
          //   fFitCh->AssignData(gMinusCh->GetN(), 1, gMinusCh->GetX(), gMinusCh->GetY());
          //   fFitCh->Eval();
          //   delete gMinusCh;
          //   Double_t mean = 0, rms = 0;
          //   Int_t nnoise = 0;
          //   for(Int_t j=0;j<NCHAN;j++){ // Channels
          //     Double_t q = m->GetData(j,i);
          //   if(q<0) continue;
          //     q -= fFitCh->GetParameter(0) + j*fFitCh->GetParameter(1) + j*j*fFitCh->GetParameter(2);
          //     m->SetData(j,i,q);
          //     if(q<100){
          //    mean += q;
          //    rms += q*q;
          //    nnoise++;
          //     }
          //   }
          //   rms = TMath::Sqrt(rms/nnoise - mean*mean/nnoise/nnoise);
          //   for(Int_t j=0;j<NALL;j++){ // Channels
          //     Double_t q = m->GetData(j,i);
          //   if(q<0) continue;
          //     if(q>2*rms){
          //    m->SetData(j,i,q);
          //     }else{
          //    m->SetData(j,i,0);
          //     }
          //   }
          //   fFitCh->Delete();
          // }
          
          
          
          Double_t Qcl = 0, Xcl = 0;
          for(Int_t j=0;j<NCHAN;j++){ // Channels
            Double_t q = m->GetData(j,i);
            if(q<0) continue;
            if(q>40){
              Qcl += q;
              Xcl += q*j;
            }
            //if(i>250 && j==200) Info("process","t = %d, ch = %d, Q = %g",i,j,q);
            hCumul[ndet]->Fill(i+0.5,j+0.5,q);
            hZs[ndet]->SetBinContent(i+1,j+1,q);
            qTotTZs += q;
            qTotZs += q;
            hChZs[ndet]->Fill(j+0.5,q);
            hTZs[ndet]->Fill(i+0.5,q);
            hTevZs[ndet]->Fill(i+0.5,q);
            hQZs[ndet]->Fill(q);
            hQevZs[ndet]->Fill(q);
          }
          if(Qcl>250 && Tmin == -10 && i>5){
            Tmin = i;
            Xmin = Xcl/Qcl;
          }
          if(Qcl>250 && Tmin > 0) Tmax = i;
          hQtotTZs[ndet]->Fill(qTotTZs);
          //hQtotTZs[ndet]->Fill(Qcl);
          hQtotTevZs[ndet]->Fill(qTotTZs);
        }               
        hTmin[ndet]->Fill(Tmin);
        hTmax[ndet]->Fill(Tmax);
        gXTmin[ndet]->SetPoint(gXTmin[ndet]->GetN(),Tmin,Xmin);

        hQtot[ndet]->Fill(qTot);
        hQtotZs[ndet]->Fill(qTotZs);
        qtot[ndet] = qTot;
        if(Tmin>fTmin) eventType = E_CREATED;
        if(Tmax-Tmin>fLTpc) eventType = E_TRAVERSING; 
        //      if(triggerType == TRIG_TT)
        if(isTT)
          hTdrift_TT[ndet]->Fill(Tmax-Tmin);
      } else {
        std::cout << "No Data Map for plane "<< ndet << std::endl;
      }
    }
  }
  
  if(qtot[0]>fQmax || qtot[1]>fQmax) eventType = E_SATURATED;
  if(qtot[0]<fQmin && qtot[1]<fQmin) eventType = E_EMPTY;

  hInfo->SetBinContent(8,triggerType);  
  //Info("process","trigger type: %d",triggerType);
  switch(eventType){
  case E_CREATED:
    hTriggerCreated->Fill(triggerType);
    fNcreated[triggerType+1]++;
    break;
  case E_TRAVERSING:
    hTriggerTraversing->Fill(triggerType);
    fNtraversing[triggerType+1]++;
    break;
  case E_EMPTY:
    hTriggerEmpty->Fill(triggerType);
    fNempty[triggerType+1]++;
    break;
  case E_SATURATED:
    hTriggerSaturated->Fill(triggerType);
    fNsaturated[triggerType+1]++;
    break;
  case E_OTHER:
    hTriggerOther->Fill(triggerType);
    fNother[triggerType+1]++;
    break;
  }
  

  hDeadTime[triggerType+1]->Fill((timestamp-timestampprev)*1e-5);
  timestampprev = timestamp;

  nSamples++;
  qtotX += qtot[0];
  qtotY += qtot[1];
  qtot2X += qtot[0]*qtot[0];
  qtot2Y += qtot[1]*qtot[1];
  if(timestamp>timestampold+fRateSample*1e8){ // every second
    Double_t freq = 1e8/(timestamp-timestampold); // Hz
    freq *= gHConfig->GetSkipEvents()+1;
    timestampold = timestamp;
    Double_t time = (timestamp+timestampold)/2e8; // seconds
    //      Info("process","freq = 1e7/(%g - %g) = %g",timestamp,timestampold,freq);
    for(Int_t i = 0; i<kNtriggerTypes; i++){
      Double_t trigCorr  = hTriggerCorrection->GetBinContent(i+1);
      // Double_t trigCorr  = 1.;
      // if(fTriggerTime[i])
      //        trigCorr = fTriggerTimeTot/fTriggerTime[i];
      // if(fTriggerDownscaling[i])
      //        trigCorr *= fTriggerDownscaling[i];
      if(gHarpoDebug>1)
        Info("process","trigCorr = %g",trigCorr);
      gRateCreated[i]->SetPoint(gRateCreated[i]->GetN(),time,fNcreated[i]*freq*trigCorr);
      //gRateCreated[i]->SetPointError(gRateCreated[i]->GetN()-1,fRateSample/2,TMath::Sqrt(fNcreated[i])*freq*trigCorr);

      gRateTraversing[i]->SetPoint(gRateTraversing[i]->GetN(),time,fNtraversing[i]*freq*trigCorr);
      //gRateTraversing[i]->SetPointError(gRateTraversing[i]->GetN()-1,fRateSample/2,TMath::Sqrt(fNtraversing[i])*freq*trigCorr);

      gRateEmpty[i]->SetPoint(gRateEmpty[i]->GetN(),time,fNempty[i]*freq*trigCorr);
      //gRateEmpty[i]->SetPointError(gRateEmpty[i]->GetN()-1,fRateSample/2,TMath::Sqrt(fNempty[i])*freq*trigCorr);

      gRateSaturated[i]->SetPoint(gRateSaturated[i]->GetN(),time,fNsaturated[i]*freq*trigCorr);
      //gRateSaturated[i]->SetPointError(gRateSaturated[i]->GetN()-1,fRateSample/2,TMath::Sqrt(fNsaturated[i])*freq*trigCorr);

      gRateOther[i]->SetPoint(gRateOther[i]->GetN(),time,fNother[i]*freq*trigCorr);
      //      gRateOther[i]->SetPointError(gRateOther[i]->GetN()-1,fRateSample/2,TMath::Sqrt(fNother[i])*freq*trigCorr);

      Int_t nTot = fNother[i]+fNtraversing[i]+fNempty[i]+fNcreated[i]+fNsaturated[i];
      if(nTot){
        gRateTotal[i]->SetPoint(gRateTotal[i]->GetN(),time,nTot*freq*trigCorr);
        //      gRateTotal[i]->SetPointError(gRateTotal[i]->GetN()-1,fRateSample/2,TMath::Sqrt(nTot)*freq*trigCorr);
        gRateTotal[i]->SetPointError(gRateTotal[i]->GetN()-1,0,0);
      }
      fNcreated[i] = 0;
      fNtraversing[i] = 0;
      fNempty[i] = 0;
      fNsaturated[i] = 0;
      fNother[i] = 0;
    }

    if(nSamples){
      qtotX /= nSamples;
      qtot2X /= nSamples;
      qtotY /= nSamples;
      qtot2Y /= nSamples;
        
      gQtotTT_X->SetPoint(gQtotTT_X->GetN(),time,qtotX);
      gQtotTT_X->SetPointError(gQtotTT_X->GetN()-1,fRateSample/2,TMath::Sqrt(qtot2X-qtotX*qtotX)/2);
      gQtotTT_Y->SetPoint(gQtotTT_Y->GetN(),time,qtotY);
      gQtotTT_Y->SetPointError(gQtotTT_Y->GetN()-1,fRateSample/2,TMath::Sqrt(qtot2Y-qtotY*qtotY)/2);
    }      
    nSamples = 0;
    qtotX = 0;
    qtotY = 0;
    qtot2X = 0;
    qtot2Y = 0;
  }

  // for(Int_t i = 0; i<gkHarpoMonNhist; i++){
  //   HarpoHistoIdentifier id = HarpoHistoIdentifier(i);
  //   GetHist(id,0)->Write(Form("h%d_0",i));
  //   GetHist(id,1)->Write(Form("h%d_1",i));
  // }
  if((nEvents-1)%fSaveRate == 0)
    mfile->Update(); 
}

void   HarpoAnalyseMonitorNew::Init()
{
  // m_hmap[0] = new HarpoMap(0);
  // m_hmap[1] = new HarpoMap(1);





  Long64_t fAddress = 0x7fef1fc55000;
  Long64_t address = 0;
  if ( ! gHConfig->Lookup("monitor.address",address) )
    Info("Init","Use default address %lld",fAddress);
  else{
    fAddress = address;
    TMapFile::SetMapAddress(fAddress);
  }

  //   TMapFile::SetMapAddress(0xb46a5000);
  //  TMapFile::SetMapAddress(0x7fef1fc55000);
  mfile = TMapFile::Create("/tmp/harpo.map","RECREATE", 50000000,
                           "Memory mapped file with monitoring histograms");
  
  // for(Int_t ndet = 0; ndet<1; ndet++){
  //   if(ndet == 0) color = kBlue;
  //   else color = kRed;
  //   MakeHist(H_,ndet,";",
  fIsZs = 1;
  Long64_t isZs = 1;
  if ( ! gHConfig->Lookup("monitor.IsZs",isZs) )
    Info("Init","Use default fIsZs %d",fIsZs);
  else
    fIsZs = isZs;

  fQmin = 10000;
  Long64_t qmin = 0;
  if ( ! gHConfig->Lookup("monitor.Qmin",qmin) )
    Info("Init","Use default fQmin %d",fQmin);
  else
    fQmin = qmin;
  
  fQmax = 300000;
  Long64_t qmax = 0;
  if ( ! gHConfig->Lookup("monitor.Qmax",qmax) )
    Info("Init","Use default fQmax %d",fQmax);
  else
    fQmax = qmax;

  fTmin = 150;
  Double_t tmin = 0;
  if ( ! gHConfig->Lookup("monitor.Tmin",tmin) )
    Info("InitCfg","Use default fTmin %g",fTmin);
  else
    fTmin = tmin;

  fLTpc = 200;
  Double_t lTpc = 0;
  if ( ! gHConfig->Lookup("monitor.LTpc",lTpc) )
    Info("InitCfg","Use default fLTpc %g",fLTpc);
  else
    fLTpc = lTpc;

  fRateSample = 200;
  Double_t rateSample = 0;
  if ( ! gHConfig->Lookup("monitor.RateSample",rateSample) )
    Info("InitCfg","Use default fRateSample %g",fRateSample);
  else
    fRateSample = rateSample;

  fSaveRate = 1;
  Long64_t saveRate = 0;
  if ( ! gHConfig->Lookup("monitor.SaveRate",saveRate) )
    Info("InitCfg","Use default fSaveRate %d",fSaveRate);
  else
    fSaveRate = saveRate;
  if(fSaveRate<1) fSaveRate = 1;

  

  hQvsT_TT[0] = new TH1F("hQvsT_TT_X","Q vs T (traversing tracks);T [tb]; Q [ADC]",511,0,511);
  hQvsT_TT[1] = new TH1F("hQvsT_TT_Y","Q vs T (traversing tracks);T [tb]; Q [ADC]",511,0,511);
  hQvsT_TT2D[0] = new TH2F("hQvsT_TT2D_X","Q vs T (traversing tracks);T [tb]; Q [ADC]",511,0,511,500,0,50000);
  hQvsT_TT2D[1] = new TH2F("hQvsT_TT2D_Y","Q vs T (traversing tracks);T [tb]; Q [ADC]",511,0,511,500,0,50000);
  hTdrift_TT[0] = new TH1F("hTdrift_TT_X",";T_{max}-T_{min} [tb]",511,0,511);
  hTdrift_TT[1] = new TH1F("hTdrift_TT_Y",";T_{max}-T_{min} [tb]",511,0,511);


  hTmin[0] = new TH1F("hTmin_X","Tmin;T_{min} [tb]",511,0,511);
  hTmin[1] = new TH1F("hTmin_Y","Tmin;T_{min} [tb]",511,0,511);

  hTmax[0] = new TH1F("hTmax_X",";T_{max} [tb]",511,0,511);
  hTmax[1] = new TH1F("hTmax_Y",";T_{max} [tb]",511,0,511);

  hT[0] = new TH1F("hT_X",";time [tb]",511,0,511);
  hT[1] = new TH1F("hT_Y",";time [tb]",511,0,511);

  hCh[0] = new TH1F("hCh_X","Channels (cumulated);channel",304,0,304);
  hCh[1] = new TH1F("hCh_Y","Channels (cumulated);channel",304,0,304);

  hQ[0] = new TH1F("hQ_X","Pixel Charge;Q_{pixel} [ADC]",4096,0,4096);
  hQ[1] = new TH1F("hQ_Y","Pixel Charge;Q_{pixel} [ADC]",4096,0,4096);

  hQev[0] = new TH1F("hQev_X",";Q_{pixel} [ADC]",4096,0,4096);
  hQev[1] = new TH1F("hQev_Y",";Q_{pixel} [ADC]",4096,0,4096);


  hTZs[0] = new TH1F("hTZs_X",";time [tb]",511,0,511);
  hTZs[1] = new TH1F("hTZs_Y",";time [tb]",511,0,511);

  hChZs[0] = new TH1F("hChZs_X","Channels (cumulated);channel",304,0,304);
  hChZs[1] = new TH1F("hChZs_Y","Channels (cumulated);channel",304,0,304);

  hQZs[0] = new TH1F("hQZs_X",";Q_{pixel} [ADC]",4096,0,4096);
  hQZs[1] = new TH1F("hQZs_Y",";Q_{pixel} [ADC]",4096,0,4096);

  hQevZs[0] = new TH1F("hQevZs_X",";Q_{pixel} [ADC]",4096,0,4096);
  hQevZs[1] = new TH1F("hQevZs_Y",";Q_{pixel} [ADC]",4096,0,4096);


  const Int_t nbins1 = 500;
  const Int_t nbins2 = 500;
  Double_t xmin = 10;
  Double_t xmin2 = 10;
  Double_t xmax = 2e5;
  Double_t xmax2 = 1e8;
  Double_t logxmin = TMath::Log(xmin);
  Double_t logxmin2 = TMath::Log(xmin2);
  Double_t logxmax = TMath::Log(xmax);
  Double_t logxmax2 = TMath::Log(xmax2);
  Double_t binwidth = (logxmax-logxmin)/nbins1;
  Double_t binwidth2 = (logxmax2-logxmin2)/nbins2;
  Double_t xbins1[nbins1+1];
  Double_t xbins2[nbins2+1];
  xbins1[0] = xmin;
  xbins2[0] = xmin2;
  for (Int_t i=1;i<=nbins1;i++)
    xbins1[i] = TMath::Exp(logxmin+i*binwidth);
  for (Int_t i=1;i<=nbins2;i++)
    xbins2[i] = TMath::Exp(logxmin2+i*binwidth2);

  hQtotT[0] = new TH1F("hQtotT_X","Cluster Charge;Q_{cluster} [ADC]",nbins1,xbins1);
  hQtotT[1] = new TH1F("hQtotT_Y","Cluster Charge;Q_{cluster} [ADC]",nbins1,xbins1);

  hQtot[0] = new TH1F("hQtot_X","Event Charge;Q_{event} [ADC]",nbins2,xbins2);
  hQtot[1] = new TH1F("hQtot_Y","Event Charge;Q_{event} [ADC]",nbins2,xbins2);

  hQtotTZs[0] = new TH1F("hQtotTZs_X",";Q_{cluster} [ADC]",nbins1,xbins1);
  hQtotTZs[1] = new TH1F("hQtotTZs_Y",";Q_{cluster} [ADC]",nbins1,xbins1);

  hQtotZs[0] = new TH1F("hQtotZs_X","Event Charge;Q_{event} [ADC]",nbins2,xbins2);
  hQtotZs[1] = new TH1F("hQtotZs_Y","Event Charge;Q_{event} [ADC]",nbins2,xbins2);

  // hQtotT[0] = new TH1F("hQtotT_X",";Q_{cluster} [ADC]",1000,0,10000);
  // hQtotT[1] = new TH1F("hQtotT_Y",";Q_{cluster} [ADC]",1000,0,10000);

  // hQtot[0] = new TH1F("hQtot_X","Q_{event} [ADC]",1000,0,1e7);
  // hQtot[1] = new TH1F("hQtot_Y","Q_{event} [ADC]",1000,0,1e7);

  // hQtotTZs[0] = new TH1F("hQtotTZs_X",";Q_{cluster} [ADC]",1000,0,10000);
  // hQtotTZs[1] = new TH1F("hQtotTZs_Y",";Q_{cluster} [ADC]",1000,0,10000);

  // hQtotZs[0] = new TH1F("hQtotZs_X","Q_{event} [ADC]",1000,0,1e7);
  // hQtotZs[1] = new TH1F("hQtotZs_Y","Q_{event} [ADC]",1000,0,1e7);

  hInfo = new TH1F("hInfo","",15,0.5,15.5);
  hPmm2 = new TH1F("hPmm2","",16,0,16);

  hTmin[0]->SetLineColor(kBlue);
  hTmax[0]->SetLineColor(kBlue+3);
  hTmin[1]->SetLineColor(kRed);
  hTmax[1]->SetLineColor(kRed+3);


  hCh[0]->SetLineColor(kBlue);
  hT[0]->SetLineColor(kBlue);
  hQtotT[0]->SetLineColor(kBlue);
  hQtot[0]->SetLineColor(kBlue);
  hQ[0]->SetLineColor(kBlue);
  hQev[0]->SetLineColor(kBlue);
  hCh[1]->SetLineColor(kRed);
  hT[1]->SetLineColor(kRed);
  hQtotT[1]->SetLineColor(kRed);
  hQtot[1]->SetLineColor(kRed);
  hQ[1]->SetLineColor(kRed);
  hQev[1]->SetLineColor(kRed);

  hQvsT_TT[0]->SetLineColor(kBlue);
  hQvsT_TT[1]->SetLineColor(kRed);
  hTdrift_TT[0]->SetLineColor(kBlue);
  hTdrift_TT[1]->SetLineColor(kRed);

  hChZs[0]->SetLineColor(kBlue);
  hTZs[0]->SetLineColor(kBlue);
  hQtotTZs[0]->SetLineColor(kBlue);
  hQtotZs[0]->SetLineColor(kBlue);
  hQevZs[0]->SetLineColor(kBlue);
  hQZs[0]->SetLineColor(kBlue);
  hChZs[1]->SetLineColor(kRed);
  hTZs[1]->SetLineColor(kRed);
  hQtotTZs[1]->SetLineColor(kRed);
  hQtotZs[1]->SetLineColor(kRed);
  hQevZs[1]->SetLineColor(kRed);
  hQZs[1]->SetLineColor(kRed);

  hChZs[0]->SetLineStyle(2);
  hTZs[0]->SetLineStyle(2);
  hQtotTZs[0]->SetLineStyle(2);
  hQtotZs[0]->SetLineStyle(2);
  hQevZs[0]->SetLineStyle(2);
  hQZs[0]->SetLineStyle(2);
  hChZs[1]->SetLineStyle(2);
  hTZs[1]->SetLineStyle(2);
  hQtotTZs[1]->SetLineStyle(2);
  hQtotZs[1]->SetLineStyle(2);
  hQevZs[1]->SetLineStyle(2);
  hQZs[1]->SetLineStyle(2);


   hTev[0] = new TH1F("hTev_X","",511,0,511);
   hTev[1] = new TH1F("hTev_Y","",511,0,511);
   hTevZs[0] = new TH1F("hTevZs_X","",511,0,511);
   hTevZs[1] = new TH1F("hTevZs_Y","",511,0,511);
   
   Int_t nbins = 1000;
   Int_t Qmax = 10000;
   hQtotTev[0] = new TH1F("hQtotTev_X","",nbins,0,Qmax);
   hQtotTev[1] = new TH1F("hQtotTev_Y","",nbins,0,Qmax);
   hQtotTevZs[0] = new TH1F("hQtotTevZs_X","",nbins,0,Qmax);
   hQtotTevZs[1] = new TH1F("hQtotTevZs_Y","",nbins,0,Qmax);

   hTev[0]->SetLineColor(kBlue);
   hTev[1]->SetLineColor(kRed);
   hQtotTev[0]->SetLineColor(kBlue);
   hQtotTev[1]->SetLineColor(kRed);
   hTevZs[0]->SetLineColor(kBlue);
   hTevZs[1]->SetLineColor(kRed);
   hQtotTevZs[0]->SetLineColor(kBlue);
   hQtotTevZs[1]->SetLineColor(kRed);
   hTevZs[0]->SetLineStyle(2);
   hTevZs[1]->SetLineStyle(2);
   hQtotTevZs[0]->SetLineStyle(2);
   hQtotTevZs[1]->SetLineStyle(2);
   

   h[0] = new TH2F("hX","X direction;time [t.b.];channel [1mm]",511,0,511,304,0,304);
   hZs[0] = new TH2F("hXZs","X direction;time [t.b.];channel [1mm]",511,0,511,304,0,304);
   h[1] = new TH2F("hY","Y direction;time [t.b.];channel [1mm]",511,0,511,304,0,304);
   hZs[1] = new TH2F("hYZs","Y direction;time [t.b.];channel [1mm]",511,0,511,304,0,304);
   hCumul[0] = new TH2F("hCumul_X","X direction;time [t.b.];channel [1mm]",511,0,511,304,0,304);
   hCumul[1] = new TH2F("hCumul_Y","Y direction;time [t.b.];channel [1mm]",511,0,511,304,0,304);

   hTriggerCorrection = new TH1F("hTriggerCorrection",";Trigger Type",kNtriggerTypes+1,-1,kNtriggerTypes);
   hTriggerCreated = new TH1F("hTriggerCreated",";Trigger Type",kNtriggerTypes+1,-1,kNtriggerTypes);
   hTriggerTraversing = new TH1F("hTriggerTraversing",";Trigger Type",kNtriggerTypes+1,-1,kNtriggerTypes);
   hTriggerEmpty = new TH1F("hTriggerEmpty",";Trigger Type",kNtriggerTypes+1,-1,kNtriggerTypes);
   hTriggerSaturated = new TH1F("hTriggerSaturated",";Trigger Type",kNtriggerTypes+1,-1,kNtriggerTypes);
   hTriggerOther = new TH1F("hTriggerOther",";Trigger Type",kNtriggerTypes+1,-1,kNtriggerTypes);


   hPmm2Size = new TH1F("hPmm2Size",";#PM",12,1,13);
   hPmm2Triggered = new TH1F("hPmm2Triggered",";#PM",12,1,13);
   hPmm2Missed = new TH1F("hPmm2Missed",";#PM",12,1,13);
   hPmm2Spectrum2D = new TH2F("hPmm2Spectrum2D",";Channel;Charge [ADC]",12,1,13,1024,0,1024);
   hPmm2Coincidence = new TH2F("hPmm2Coincidence",";Channel;Channel",12,1,13,12,1,13);

   gXTmin[0] = new TGraph();
   mfile->Add(gXTmin[0],"gXTmin");
   gXTmin[1] = new TGraph();
   mfile->Add(gXTmin[1],"gYTmin");
   for(Int_t i = 0; i<kNtriggerTypes+1; i++){
     gRateCreated[i] = new TGraphErrors();
     mfile->Add(gRateCreated[i],Form("gRateCreated_%d",i));
     gRateTraversing[i] = new TGraphErrors();
     mfile->Add(gRateTraversing[i],Form("gRateTraversing_%d",i));
     gRateEmpty[i] = new TGraphErrors();
     mfile->Add(gRateEmpty[i],Form("gRateEmpty_%d",i));
     gRateSaturated[i] = new TGraphErrors();
     mfile->Add(gRateSaturated[i],Form("gRateSaturated_%d",i));
     gRateOther[i] = new TGraphErrors();
     mfile->Add(gRateOther[i],Form("gRateOther_%d",i));
     gRateTotal[i] = new TGraphErrors();
     mfile->Add(gRateTotal[i],Form("gRateTotal_%d",i));
     fNcreated[i] = 0;
     fNtraversing[i] = 0;
     fNempty[i] = 0;
     fNsaturated[i] = 0;
     fNother[i] = 0;
     hDeadTime[i] = new TH1F(Form("hDeadTime_%d",i),"; #deltaT [ms]; counts",500,0,500);
   }

   gQtotTT_X = new TGraphErrors();
   mfile->Add(gQtotTT_X,"gQtotTT_X");
   gQtotTT_Y = new TGraphErrors();
   mfile->Add(gQtotTT_Y,"gQtotTT_Y");
   qtotX = 0;
   qtotY = 0;
   qtot2X = 0;
   qtot2Y = 0;
   nSamples = 0;

   mfile->Print();


   timestampold = 0;
}


TH1* HarpoAnalyseMonitorNew::GetHist(HarpoHistoIdentifiers ihist, Int_t ndet)
{

  // Int_t ndet = i%2;
  // Int_t ihist = i/2;

  switch(ihist){
  case H_MAP: return h[ndet];
  case H_MAP_ZS: return hZs[ndet];
  case H_Q: return hQ[ndet];
  case H_Q_EV: return hQev[ndet];
  case H_T: return hT[ndet];
  case H_CH: return hCh[ndet];
  case H_QTOT_T: return hQtotT[ndet];
  case H_QTOT: return hQtot[ndet];
  case H_Q_ZS: return hQZs[ndet];
  case H_Q_EV_ZS: return hQevZs[ndet];
  case H_T_ZS: return hTZs[ndet];
  case H_CH_ZS: return hChZs[ndet];
  case H_QTOTT_ZS: return hQtotTZs[ndet];
  case H_QTOT_ZS: return hQtotZs[ndet];
  case H_T_EV: return hTev[ndet];
  case H_T_EV_ZS: return hTevZs[ndet];
  case H_QTOTT_EV: return hQtotTev[ndet];
  case H_QTOTT_EV_ZS: return hQtotTevZs[ndet];
  default: return 0;
  }

  return 0;
}



void   HarpoAnalyseMonitorNew::Save(char * /* mode */)
{
  TString * hstFile = gHConfig->GetHistFile();
  // if ( hstFile == NULL ) {
  //   std::cout << gHConfig->GetProgramName() << " " <<
  //     "No Hist File name gived, use default" << std::endl;
  //   hstFile = new TString("onlinemonitor.root");
  // }
  if ( hstFile == NULL ) {
    const char* dir = gSystem->ExpandPathName("$HARPO_ANA_DIR/monitor");
    //const char* dir = "$HARPO_ANA_DIR/monitor";
    if(gSystem->AccessPathName(dir))
      hstFile = new TString(Form("onlinemonitor%lli.root",gHConfig->GetRunNo()));
    else    
      hstFile = new TString(Form("%s/onlinemonitor%lli.root",dir,gHConfig->GetRunNo()));
    Info("Save","No Hist File name given, use default %s",hstFile->Data());
  }
  
  TFile *hf = new TFile(hstFile->Data(),"RECREATE");
  hQvsT_TT[0]->Write();
  hQvsT_TT[1]->Write();
  hQvsT_TT2D[0]->Write();
  hQvsT_TT2D[1]->Write();
  hTdrift_TT[0]->Write();
  hTdrift_TT[1]->Write();
  hTmin[0]->Write();
  hTmin[1]->Write();
  hTmax[0]->Write();
  hTmax[1]->Write();
  hT[0]->Write();
  hT[1]->Write();
  hCh[0]->Write();
  hCh[1]->Write();
  hQ[0]->Write();
  hQ[1]->Write();
  hQev[0]->Write();
  hQev[1]->Write();
  hTZs[0]->Write();
  hTZs[1]->Write();
  hChZs[0]->Write();
  hChZs[1]->Write();
  hQZs[0]->Write();
  hQZs[1]->Write();
  hQevZs[0]->Write();
  hQevZs[1]->Write();
  hQtotT[0]->Write();
  hQtotT[1]->Write();
  hQtot[0]->Write();
  hQtot[1]->Write();
  hQtotTZs[0]->Write();
  hQtotTZs[1]->Write();
  hQtotZs[0]->Write();
  hQtotZs[1]->Write();
  hInfo->Write();
  hPmm2->Write();
  hTev[0]->Write();
  hTev[1]->Write();
  hTevZs[0]->Write();
  hTevZs[1]->Write();
  hQtotTev[0]->Write();
  hQtotTev[1]->Write();
  hQtotTevZs[0]->Write();
  hQtotTevZs[1]->Write();
  h[0]->Write();
  hZs[0]->Write();
  h[1]->Write();
  hZs[1]->Write();
  hCumul[0]->Write();
  hCumul[1]->Write();
  hTriggerCorrection->Write();
  hTriggerCreated->Write();
  hTriggerTraversing->Write();
  hTriggerEmpty->Write();
  hTriggerSaturated->Write();
  hTriggerOther->Write();
  hPmm2Size->Write();
  hPmm2Triggered->Write();
  hPmm2Missed->Write();
  hPmm2Spectrum2D->Write();
  hPmm2Coincidence->Write();




   hf->Close();
}




Double_t HarpoAnalyseMonitorNew::TruncMean(TArrayD* vect, 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;
  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]);
  if(tHigh-tLow) {
    return truncMean/(tHigh-tLow);
  }
  else return 0;
}