HarpoAnalysePmm2Monitor.cxx

Go to the documentation of this file.

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

#include "MakeNiceHisto.h"

#include "TFile.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TGraph.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TArrayD.h"
#include "TMath.h"
#include "TLine.h"
#include "TGaxis.h"
#include "TROOT.h"

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

ClassImp(HarpoAnalysePmm2Monitor)

void   HarpoAnalysePmm2Monitor::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   HarpoAnalysePmm2Monitor::process()
{
  //  Bool_t drawEvent = kFALSE;
  nEvents++;
  //  std::cout << " Processing Event " << nEvents<< std::endl;
  if(nEvents%10000==0)
    Info("process"," Processing Event %ld",nEvents);
  
  // TCanvas *ct = (TCanvas *)gROOT->GetListOfCanvases()->FindObject("cPmm2");
  // if(!ct){
  //   ct = new TCanvas("cPmm2","Monitor PMm2",0,0,1200,800);
  //   ct->Clear();
  //   ct->Divide(2);//,0.001,0.001);
  // }
  
  if (!gHDetSet->isExist(PMM2)) return;
  Pmm2Event *anaEvt = static_cast<Pmm2Event *>(fEvt->GetDetEvent(PMM2));

  fHistInfo->Reset();
  fHistInfo->SetBinContent(1,gHConfig->GetRunNo());
  fHistInfo->SetBinContent(2,fEvt->GetHeader()->GetEvtNo());
  Double_t timestamp = fEvt->GetTimeStamp(PMM2)*10;  
  if(gHarpoDebug>1) Info("process","%g",timestamp);
  if(timestamp==0) return;  
  //  fHistInfo->SetBinContent(3,fEvt->GetTimeStamp(PMM2)*10);
  fHistInfo->SetBinContent(3,timestamp);
  int imes;
  int esize = anaEvt->GetHeader()->eventSize;
  //  int ievnum = anaEvt->GetHeader()->eventNumb;
  //  Info("process","ievnum = %d, esize = %d",ievnum,esize);
  Pmm2MesVect * pm = anaEvt->GetMesurements();
  
    
  fHistSize->Fill(esize);
  //int itsold = 0;
  int hit[16];
  for(int ch = 0; ch<16; ch++) hit[ch] = 0;
  for(imes=0;imes<esize;imes++) {
    //    Pmm2MesVect &m = *pm;
    
    int ch = pm->at(imes).getChNum(); // channel
    hit[ch] = 1;
    rate[ch] += 1;
    //fHistRatePmm2->Fill(ch);
    fHistRatePmm2Tot->Fill(ch);
  }
  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]);
  
  fHistNscint->Fill(nScint);
  
  bool trigger = nScint > 1;
  if(trigger){
    fHistSizeTrig->Fill(esize);
  }

  for(int ch = 0; ch<16; ch++){
    int ch2 = ch%2 ? ch+1 : ch-1;
    if(trigger){
      if(hit[ch2]){
        if(!hit[ch])
          fHistMissed->Fill(ch);
        fHistTriggered->Fill(ch);
      }
    }
  }

  TString str("");
  str += nEvents;
  str += " : ";
  for(imes=0;imes<esize;imes++) {
    // i cann't use operator[] on pointer
    Pmm2MesVect &m = *pm;
    
    int ch      = pm->at(imes).getChNum(); // channel
    int iadc    = m[imes].getCharge(); // Charge in adc counts
    int trigtype= m[imes].GetNU0(); // Fine Time

    str += ch;
    str += " ";
    str += iadc;
    str += " ";
    str += trigtype;
    str += " | ";

    if(ch<0 || ch>15){
      Info("process","Unknown PMM2 channel %d",ch);
      continue;
    }
    fHistSpectrum[ch]->Fill(iadc);
    fHistSpectrum2D->Fill(ch,iadc);
    if(trigger){
      fHistSpectrumTrig[ch]->Fill(iadc);
      fHistSpectrumTrig2D->Fill(ch,iadc);
    }
  } // Digi Loop
  if(gHarpoDebug>1)
    Info("process"," %s",str.Data());

  // if(timestamp<itsold-1){
  //   Info("process","%d - %d = %d",timestamp, itsold, timestamp - itsold);
  //   Info("","************************************************");
  //   startRun = kTRUE;
  // }

  for(int ch1 = 0; ch1<16; ch1++){
    for(int ch2 = 0; ch2<16; ch2++){
      if(hit[ch1] && hit[ch2]){
        fHistCoincidence->Fill(ch1,ch2);
      }
      if(hit[ch1])
        fHistCoincidenceNorm->Fill(ch1,ch2);
    }
  }


  if(timestamp>timestampold+fRateSample*1e8){ // every fRateSample seconds
    Double_t freq = 1e8/(timestamp-timestampold); // Hz
    timestampold = timestamp;
    Double_t time = timestamp/1e8; // seconds
    //    fTimeTot += time;
    Double_t ratetot = 0;
    for(Int_t i = 0; i<16; i++){
      // if(fHistRatePmm2->GetBinContent(i+1) == 0) continue;
      // gRatePmm2[i]->SetPoint(gRatePmm2[i]->GetN(),time,fHistRatePmm2->GetBinContent(i+1)*freq);
      // gRatePmm2[i]->SetPointError(gRatePmm2[i]->GetN()-1,fRateSample/2,TMath::Sqrt(fHistRatePmm2->GetBinContent(i+1))*freq);
      // ratetot += fHistRatePmm2->GetBinContent(i+1);
      //      if(rate[i] <= 0) continue;
      if(rate[i]){
        gRatePmm2[i]->SetPoint(gRatePmm2[i]->GetN(),time,rate[i]*freq);
        gRatePmm2[i]->SetPointError(gRatePmm2[i]->GetN()-1,fRateSample/2,TMath::Sqrt(rate[i])*freq);
        ratetot += rate[i];
      }
      fHistRatePmm2->SetBinContent(i+1,rate[i]*freq);
      //      fHistRatePmm2->SetBinContent(i+1,rate[i]+ fHistRatePmm2->GetBinContent(i+1));
      fHistRatePmm2->SetBinError(i+1,TMath::Sqrt(rate[i])*freq);
      rate[i] = 0;
    }
    //Info("process","Ratetot %g",ratetot);
    //    fHistRatePmm2->Scale(freq);
    //    fHistRatePmm2->Reset();
    gRatePmm2Tot->SetPoint(gRatePmm2Tot->GetN(),time,ratetot*freq);
    gRatePmm2Tot->SetPointError(gRatePmm2Tot->GetN()-1,fRateSample/2,TMath::Sqrt(ratetot)*freq);
  }

  fTimeTot = timestamp;
  gTimeVsEvent->SetPoint(gTimeVsEvent->GetN(),nEvents,timestamp/1e8);

  if(fDisplay%2 == 0 && (nEvents-1)%10 == 0)
    mfile->Update();


}

void   HarpoAnalysePmm2Monitor::Init()
{
  
  // Initialise histograms here
  
  fTimeTot = 0;

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


  fDisplay = 0;
  Long64_t display = 0;
  if ( ! gHConfig->Lookup("monitor.displayPmm2",display) )
    Info("Init","Use default address %d",fDisplay);
  else
    fDisplay = display;

  if(fDisplay%2 == 0){
    Long64_t fAddress = 0x7ff3ab9f9000;
    Long64_t address = 0;
    if ( ! gHConfig->Lookup("monitor.addressPmm2",address) )
      Info("Init","Use default address 0x%llx",fAddress);
    else{
      fAddress = address;
      TMapFile::SetMapAddress(fAddress);  
    }
    //  TMapFile::SetMapAddress(0x7ff3ab9f9000);
    mfile = TMapFile::Create("/tmp/harpoPmm2.map","RECREATE", 10000000,
                             "Memory mapped file with monitoring histograms");
  }

  for(Int_t ch = 0; ch<16; ch++){
    fHistSpectrum[ch] = new TH1F(Form("fHistSpectrum%d",ch),"; signal [ADC]",1024,0,1024);
    fHistSpectrum[ch]->SetLineColor(51 + ch*49./15);
    
    fHistSpectrumTrig[ch] = new TH1F(Form("fHistSpectrumTrig%d",ch),"",1024,0,1024);
    fHistSpectrumTrig[ch]->SetLineColor(51 + ch*49./15);
  }
  fHistSpectrum2D = new TH2F("fHistSpectrum2D","; channel; signal [ADC]",16,0,16,1024,0,1024);
  fHistSize = new TH1F("fHistSize","; Number of PM hit",16,0,16);
  
  fHistSpectrumTrig2D = new TH2F("fHistSpectrumTrig2D","",16,0,16,1024,0,1024);
  fHistSizeTrig = new TH1F("fHistSizeTrig","",16,0,16);
  
  
  fHistInfo = new TH1F("fHistInfo","",10,0,10);
  fHistNscint = new TH1F("fHistNscint","",6,0,6);
  fHistMissed = new TH1F("fHistMissed","; channel; efficiency",16,0,16);
  fHistTriggered = new TH1F("fHistTriggered","; channel; efficiency",16,0,16);
  fHistCoincidence = new TH2F("fHistCoincidence","",16,0,16,16,0,16);
  fHistCoincidenceNorm = new TH2F("fHistCoincidenceNorm","",16,0,16,16,0,16);
  fHistRatePmm2 = new TH1F("fHistRatePmm2","; channel; rate",16,0,16);
  fHistRatePmm2Tot = new TH1F("fHistRatePmm2Tot","; channel; rate",16,0,16);
  

  for(Int_t i = 0; i<16; i++){
    gRatePmm2[i] = new TGraphErrors();
    if(fDisplay%2 == 0)
      mfile->Add(gRatePmm2[i],Form("gRatePmm2_%d",i));
    rate[i] = 0;
  }
  timestampold = 0;
  gRatePmm2Tot = new TGraphErrors();
  if(fDisplay%2 == 0)
    mfile->Add(gRatePmm2Tot,"gRatePmm2Tot");
 
  fRateSample = 200;
  Double_t rateSample = 0;
  if ( ! gHConfig->Lookup("monitor.RateSample",rateSample) )
    Info("InitCfg","Use default fRateSample %g",fRateSample);
  else
    fRateSample = rateSample;

  gTimeVsEvent = new TGraphErrors();

  itscorr = 0;
  itsold = 0;
  itsold2 = 0;
  nEventsTrig = 0;
  startRun = kFALSE;
  
  if(fDisplay%2==0)
    mfile->Print();
  

  
}



void   HarpoAnalysePmm2Monitor::Save(char * /* mode */)
 {


   if(fDisplay==0) return;


   if(gHarpoDebug>0) Info("Save","Displaying histograms");

   gStyle->SetOptStat(0);
   
  TLegend* lCh = new TLegend(0.75,0.4,0.9,0.9);
  Int_t nSignal = 0;
  for(Int_t ch = 0; ch<16; ch++){
    fHistSpectrum[ch]->Rebin(8);
    fHistSpectrumTrig[ch]->Rebin(8);
    lCh->AddEntry(fHistSpectrum[ch],Form("Ch %d",ch));
    nSignal += fHistSize->GetBinContent(ch+1)*ch;
  }

  TLegend* l = new TLegend(0.6,0.6,0.9,0.9);


  // TCanvas* ct0 = new TCanvas("cPmm2","Monitor PMm2",0,0,1200,800);
  // ct0->Clear();
  // ct0->Divide(2,2);//,0.001,0.001);
  
  // TH1F* hDummyEvt = new TH1F("hDummyEvt",";event number; time stamp",nEvents/10,0,nEvents); 
  // hDummyEvt->SetMaximum(16777216);
  // ct0->cd(1);
  // hDummyEvt->DrawCopy();
  // fGraphTimeVsEvent->Draw("Psame");

  // ct0->cd(2);
  // fHistTime2D->Draw("colz");
  // ct0->GetPad(4)->SetLogz();

  // ct0->cd(3);
  // hDummyEvt->DrawCopy();
  // fGraphTimeVsEventTrig->Draw("Psame");

  // ct0->cd(4);
  // fHistTimeTrig2D->Draw("colz");
  // ct0->GetPad(4)->SetLogz();
  
  // ct0->Update();
  TCanvas* ct1 = new TCanvas("cPmm2_1","Monitor PMm2",0,0,1200,800);
  ct1->Clear();
  ct1->Divide(2,2);//,0.001,0.001);
  
  TH1F* hThr = new TH1F("hThr","",16,0,16);
  TH1F* hMean = new TH1F("hMean","",16,0,16);
  TH1F* hTruncMax = new TH1F("hTruncMax","",16,0,16);
  TH1F* hTruncMin = new TH1F("hTruncMin","",16,0,16);
  TH1F* hSat = new TH1F("hSat","",16,0,16);

  ct1->cd(1);
  Double_t max = 1;
  Double_t meanThr = 0, normThr = 0;
  for(Int_t ch = 0; ch<16; ch++){
    if(fHistSpectrum[ch]->GetMaximum()>max)
      max = fHistSpectrum[ch]->GetMaximum();
    if(ch)
       fHistSpectrum[ch]->Draw("same");
    else
      fHistSpectrum[ch]->Draw();

    Double_t thr = GetThr(fHistSpectrum[ch]);
    if(!thr) continue;
    meanThr += thr;
    normThr++;
    // TLine* line = new TLine(thr,0.5,thr,1.1*max);
    // line->SetLineStyle(2);
    // line->SetLineColor(fHistSpectrum[ch]->GetLineColor());
    // line->Draw();
    hThr->SetBinContent(ch+1,thr);
  }
  //   fHistSpectrum[0]->GetYaxis()->SetRangeUser(0,1.1*max);
  fHistSpectrum[0]->SetMaximum(1.1*max);
  fHistSpectrum[0]->SetMinimum(0.5);
  ct1->GetPad(1)->SetLogy();
  lCh->Draw();  

  ct1->cd(2);
  fHistSpectrum2D->Draw("colz");
  ct1->GetPad(2)->SetLogz();
  hThr->SetLineWidth(3);
  hThr->SetLineStyle(2);
  hThr->Draw("same");  
  TLatex* latex = new TLatex(); 
  latex->SetTextFont(132);
  latex->SetTextAlign(12);
  latex->SetTextSize(0.07);
  latex->DrawLatex(1,900,Form("mean Thr %g",meanThr/normThr));

  ct1->cd(3);
  max = 1;
  Double_t meanMean = 0, normMean = 0;
  for(Int_t ch = 0; ch<16; ch++){
    if(fHistSpectrumTrig[ch]->GetMaximum()>max)
      max = fHistSpectrumTrig[ch]->GetMaximum();
    if(ch)
       fHistSpectrumTrig[ch]->Draw("same");
    else
      fHistSpectrumTrig[ch]->Draw();


    Double_t thr = GetThr(fHistSpectrum[ch]);
    if(!thr) continue;
    // TLine* line = new TLine(thr,0.5,thr,1.1*max);
    // line->SetLineWidth(2);
    // line->SetLineStyle(2);
    // line->SetLineColor(fHistSpectrum[ch]->GetLineColor());
    // line->Draw();

    Double_t truncMax = 0, truncMin = 0;
    Double_t mean = GetTruncMean(fHistSpectrumTrig[ch],truncMin,truncMax,0.1,0.6);
    Info("Save","mean ch %d: %g",ch,mean);
    meanMean += mean;
    normMean++;
    hMean->SetBinContent(ch+1,mean);
    hTruncMin->SetBinContent(ch+1,truncMin);
    hTruncMax->SetBinContent(ch+1,truncMax);
    // TLine* lineM = new TLine(truncMax,0.5,truncMax,1.1*max);
    // lineM->SetLineWidth(2);
    // lineM->SetLineColor(fHistSpectrum[ch]->GetLineColor());
    // lineM->Draw();


    fHistSpectrumTrig[ch]->GetXaxis()->SetRangeUser(150,1024);
    Double_t nEvTrig = fHistSpectrumTrig[ch]->Integral();
    fHistSpectrumTrig[ch]->GetXaxis()->SetRangeUser(600,1024);
    Double_t nEvSat = fHistSpectrumTrig[ch]->Integral();
    fHistSpectrumTrig[ch]->GetXaxis()->SetRangeUser(0,1024);

    if(nEvTrig)
      hSat->SetBinContent(ch+1,nEvSat/nEvTrig);


    // if(nEvTrig)
    //   Info("Save","Taux saturation ch %d: %g", ch, nEvSat/nEvTrig);

  }

  //   fHistSpectrum[0]->GetYaxis()->SetRangeUser(0,1.1*max);
  fHistSpectrumTrig[0]->SetMaximum(1.1*max);
  fHistSpectrumTrig[0]->SetMinimum(0.5);
  ct1->GetPad(3)->SetLogy();
  lCh->Draw();  
  
  ct1->cd(4);
  fHistSpectrumTrig2D->Draw("colz");
  ct1->GetPad(4)->SetLogz();
  hMean->SetLineWidth(3);
  hMean->Draw("same");
  hTruncMin->SetLineStyle(2);
  hTruncMin->SetLineWidth(3);
  hTruncMin->Draw("same");
  hTruncMax->SetLineStyle(2);
  hTruncMax->SetLineWidth(3);
  hTruncMax->Draw("same");
  latex->DrawLatex(1,900,Form("mean Mean %g",meanMean/normMean));

  TCanvas* ct2 = new TCanvas("cPmm2_2","Monitor time PMm2",0,0,1200,800);
  ct2->Divide(2,2);
  ct2->cd(1);
  Info("Save","fTimeTot %g, timestampold %g",fTimeTot,timestampold);
  Double_t freq = 1e8/(fTimeTot-timestampold); // Hz
  fTimeTot /= 1e8;
  Double_t ratetot = 0;
  for(Int_t i = 0; i<16; i++){
    gRatePmm2[i]->SetPoint(gRatePmm2[i]->GetN(),fTimeTot,rate[i]*freq);
    gRatePmm2[i]->SetPointError(gRatePmm2[i]->GetN()-1,fRateSample/2,TMath::Sqrt(rate[i])*freq);
    ratetot += rate[i];
    fHistRatePmm2->SetBinContent(i+1,rate[i]*freq);
    fHistRatePmm2->SetBinError(i+1,TMath::Sqrt(rate[i])*freq);
    rate[i] = 0;
  }
  gRatePmm2Tot->SetPoint(gRatePmm2Tot->GetN(),fTimeTot,ratetot*freq);
  gRatePmm2Tot->SetPointError(gRatePmm2Tot->GetN()-1,fRateSample/2,TMath::Sqrt(ratetot)*freq);
  TH1F* hDummyTime = new TH1F("hDummyTime",";time [s]; rate [Hz]",1000,0,fTimeTot);
  hDummyTime->SetMaximum(gRatePmm2Tot->GetYaxis()->GetXmax());
  // fGraphTimeVsEventCorr->GetHistogram()->SetTitle(";timestamp;event");
  hDummyTime->DrawCopy();
  gRatePmm2Tot->Draw("Psame");
  for(Int_t i = 0; i<16; i++){
    if(!gRatePmm2[i]->GetN()) continue;
    Int_t color = 51 + i*49./17;
    gRatePmm2[i]->SetLineColor(color);
    gRatePmm2[i]->Draw("PSAME");
  }

  ct2->cd(3);
  fHistSize->Draw();
  fHistSizeTrig->SetLineColor(kRed);
  fHistSizeTrig->Draw("same");
  l->AddEntry(fHistSize,"All events","L");
  l->AddEntry(fHistSizeTrig,"Triggered events","L");
  l->Draw();


  Double_t N3 = fHistSizeTrig->GetBinContent(4);
  Double_t N4 = fHistSizeTrig->GetBinContent(5);
  Double_t N5 = fHistSizeTrig->GetBinContent(6);

  Double_t y = 0.5*fHistSize->GetMaximum(), x = 8;
  Double_t dy = 0.1*fHistSize->GetMaximum();
  latex->DrawLatex(x,y,Form("N_{4}/(N_{3}+N_{4}) = %.2g",N4/(N3+N4)));
  y -= dy;
  latex->DrawLatex(x,y,Form("N_{3}/(N_{3}+N_{4}) = %.2g",N3/(N3+N4)));
  y -= dy;
  latex->DrawLatex(x,y,Form("N_{5}/(N_{5}+N_{4}) = %.2g",N5/(N5+N4)));
  y -= dy;
  latex->DrawLatex(x,y,Form("N_{trig}/N_{tot} = %.2g",fHistSizeTrig->GetEntries()/fHistSize->GetEntries()));
  y -= dy;

  ct2->cd(4);
  fHistCoincidenceNorm->Scale(1./fHistCoincidenceNorm->GetMaximum());
  for(int ch1 = 0; ch1<16; ch1++){
    for(int ch2 = 0; ch2<16; ch2++){
      if(fHistCoincidenceNorm->GetBinContent(ch1+1,ch2+1) == 0)
        fHistCoincidenceNorm->SetBinContent(ch1+1,ch2+1,1);
    }
  }
  fHistCoincidence->Divide(fHistCoincidenceNorm);
  fHistCoincidence->Draw("colz");
  ct2->GetPad(4)->SetLogz();
  ct2->cd(2);
  // fHistCoincidenceNorm->Draw("colz");
  // fHistNscint->Draw();
  // hSat->SetLineColor(kRed);
  // hSat->Draw("same");
  // fHistRatePmm2Tot->DrawCopy();
  // fHistRatePmm2->DrawCopy("same"); 
  
  TH1F* fHistEfficiency = (TH1F*)fHistTriggered->Clone();
  fHistEfficiency->SetTitle("PM efficiency;Channel;Efficiency [%]");
  fHistEfficiency->Add(fHistMissed,-1);
  fHistEfficiency->Divide(fHistTriggered);
  fHistEfficiency->Scale(100);
  fHistEfficiency->SetMinimum(0);
  fHistEfficiency->SetMaximum(100);
  MakeNice1dHisto(fHistEfficiency,ct2->GetPad(2));
  
  

  Double_t scaletmp = 80./fHistRatePmm2Tot->GetMaximum();
  fHistRatePmm2Tot->Scale(scaletmp);   // Nevents -> arbitrary
  fHistRatePmm2Tot->SetLineColor(kRed);
  fHistRatePmm2->SetLineColor(kRed+2);
  fHistRatePmm2->SetLineStyle(2);
  fHistRatePmm2->Scale(fTimeTot); // Hz -> Nevents
  fHistRatePmm2->Scale(scaletmp);   // Nevents -> arbitrary
  // draw an axis on the right side
  TGaxis *axis = new TGaxis(16,0,
                            16,100,
                            0,100./scaletmp/fTimeTot,510,"+L");
  axis->SetLineColor(kRed);
  axis->SetTextColor(kRed);
  axis->SetLabelColor(kRed);
  axis->SetTitle("Rate [Hz]");
  axis->SetTitleSize(0.07);
  axis->SetTitleOffset(0.7);
  axis->SetTitleFont(132);
  axis->SetLabelFont(132);
  axis->SetLabelSize(0.06);
  axis->Draw();
  MakeNiceHisto(fHistRatePmm2Tot,ct2->GetPad(2),"same");
  MakeNiceHisto(fHistRatePmm2,ct2->GetPad(2),"same");
  ct2->GetPad(2)->SetRightMargin(0.12);
  

  ct2->cd(3);
  gTimeVsEvent->Draw("A*");
  
   TString * hstFile = gHConfig->GetHistFile();
   if ( hstFile == NULL ) {
     std::cout << gHConfig->GetProgramName() << " " <<
       "No Hist File name given, use default" << std::endl;
     hstFile = new TString(Form("outputs/pmm2%lli.root",gHConfig->GetRunNo()));
   }


   TFile *hf = new TFile(hstFile->Data(),"RECREATE");
   // Save histograms here

   for(Int_t ch = 0; ch<16; ch++){
     fHistSpectrum[ch]->Write();
     fHistSpectrumTrig[ch]->Write();
   }
   fHistSpectrum2D->Write();
   //   fHistTime2D->Write();
   fHistSize->Write();

   fHistSpectrumTrig2D->Write();
   //  fHistTimeTrig2D->Write();
   fHistSizeTrig->Write();


   fHistNscint->Write();
   fHistCoincidence->Write();
   fHistCoincidenceNorm->Write();



   hf->Close();
   printf("fNewFile %s closed\n", hstFile->Data() );
 
 }


Double_t HarpoAnalysePmm2Monitor::GetThr(TH1F* h)
{
  if(h->GetEntries()<10) return 0;

  h->GetXaxis()->SetRangeUser(150,1000);
  Int_t maxbin = h->GetMaximumBin();
  Double_t max = h->GetBinContent(maxbin);
  h->GetXaxis()->SetRangeUser(0,1024);

  Int_t nbins = h->GetXaxis()->GetNbins();

  Int_t i = 0;
  while(i<nbins){
    i++;
    if(h->GetXaxis()->GetBinCenter(i)<150) continue;
    Double_t val = h->GetBinContent(i);
    if(val>max/2) break;
  }

  Double_t thr = 0, W = 0;
  for(Int_t k = -1; k<2; k++){
    Double_t w = h->GetBinContent(i+k) - h->GetBinContent(i+k-1);
    thr += w*h->GetXaxis()->GetBinCenter(i+k);
    W += w;
  }
  thr = h->GetXaxis()->GetBinCenter(i);
  W = 1;

  // if(W)
  //   Info("GetThr","%g / %g = %g",thr, W, thr/W);

  if(W)
    return thr/W;
  else
    return -1;

}

Double_t HarpoAnalysePmm2Monitor::GetTruncMean(TH1F* h, Double_t &min, Double_t &max, Double_t tl, Double_t th)
{

  Int_t nbins = h->GetXaxis()->GetNbins();
  Int_t ntot = h->GetEntries();
  Int_t nt = 0, norm =0;
  Double_t truncmean = 0;

  // TArrayD* arr = new TArrayD(ntot);
  //  Int_t l = 0;
  min = 150;
  for(Int_t i = 1; i<=nbins; i++){
    Double_t q = h->GetXaxis()->GetBinCenter(i);
    Int_t n = h->GetBinContent(i);
    nt += n;
    if(nt>th*ntot) break;
    if(nt<tl*ntot){
      min = q;
      continue;
    }
    if(h->GetXaxis()->GetBinCenter(i)<150) {}; // ??? empty if ?????
    truncmean += q*n;
    norm += n;
    max = q;
    // for(Int_t k = 0; k<n; k++)
    //   arr->AddAt(q,l++);
  }

  return truncmean/norm;
  // return TruncMean(arr,tl,th);
}


Double_t HarpoAnalysePmm2Monitor::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]);
  delete[] index;
  if(tHigh-tLow) {
    return truncMean/(tHigh-tLow);
  }
  else return 0;
}