mcr/calib/tofino/StabilityFromPeakTable.C

/*
======================================================================
  StabilityFromPeakTable(const TString &runtype,const Int_t &run_num_start,const Int_t &run_num_stop, const TString &directory, const TString &calib_filename, const TString &root_filename, Bool_t estimate_threshold)
    Read Peak Tables and make histograms about trends of pedestal, Bhabha
    and Kaon peaks.
======================================================================
*/
#include "TH2.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TApplication.h"
#include "TTree.h"
#include "TF1.h"
#include "TNamed.h"
#include "TGraph.h"
#include "TProfile.h"

#include "TFndHdt.h"
//#include "FROOT.h"

using namespace FROOT;

TH1F *p_ped_e[12];
TH1F *p_ped_p[12];
TH1F *p_bha_e[12];
TH1F *p_bha_p[12];
TH1F *p_kaon_e[12];
TH1F *p_kaon_p[12];
TH1F *p_eq_e[12];
TH1F *p_eq_p[12];

TH1F *p_EQbha_e[12];
TH1F *p_EQbha_p[12];
//TH1F *p_EQbha_sqr[12];
TH1F *p_EQkaon_e[12];
TH1F *p_EQkaon_p[12];
//TH1F *p_EQkaon_sqr[12];

TH1F *p_EQ_e[12];
TH1F *p_EQ_p[12];

TH1F *p_OverKaonThreshold[12];

Int_t CurrentRunNumber;

Int_t KaonThreshold[K_N_TOFINO_SLABS];
Int_t NevOverThreshold[K_N_TOFINO_SLABS];
Int_t NevAllKaon[K_N_TOFINO_SLABS];

Float_t ped_table_e[K_N_TOFINO_SLABS];
Float_t ped_table_p[K_N_TOFINO_SLABS];
Float_t eq_table_e[K_N_TOFINO_SLABS];
Float_t eq_table_p[K_N_TOFINO_SLABS];

TFile *fHdtFile = 0;

TFndHdt *fCurHdtEv = 0;
TTree *fHdtTree = 0;

void DeclareProfiles(Int_t const &run_num_start, Int_t const &run_num_stop){

  TString hname;
  TString htitle;

  Int_t N=run_num_stop-run_num_start+2;

  for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
    hname.Form("PedE_%d",sl+1);
    htitle.Form("Tofino Pedestal slab: %d e",sl+1);
    p_ped_e[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("PedP_%d",sl+1);
    htitle.Form("Tofino Pedestal slab: %d p",sl+1);
    p_ped_p[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("BhaE_%d",sl+1);
    htitle.Form("Tofino Bhabha Peak slab: %d e",sl+1);
    p_bha_e[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("BhaP_%d",sl+1);
    htitle.Form("Tofino Bhabha Peak slab: %d p",sl+1);
    p_bha_p[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("EqBhaE_%d",sl+1);
    htitle.Form("Tofino Bhabha Peak (eq.) slab: %d e",sl+1);
    p_EQbha_e[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("EqBhaP_%d",sl+1);
    htitle.Form("Tofino Bhabha Peak (eq.) slab: %d p",sl+1);
    p_EQbha_p[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

//      hname.Form("EqBhaSqr_%d",sl+1);
//      htitle.Form("Tofino Bhabha Peak slab sqrt(e*p) (eq.): %d",sl+1);
//      p_EQbha_sqr[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("KaonE_%d",sl+1);
    htitle.Form("Tofino Kaon Peak slab: %d e",sl+1);
    p_kaon_e[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("KaonP_%d",sl+1);
    htitle.Form("Tofino Kaon Peak slab: %d p",sl+1);
    p_kaon_p[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("EqKaonE_%d",sl+1);
    htitle.Form("Tofino Kaon Peak (eq.) slab: %d e",sl+1);
    p_EQkaon_e[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("EqKaonP_%d",sl+1);
    htitle.Form("Tofino Kaon Peak slab (eq.): %d p",sl+1);
    p_EQkaon_p[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

//      hname.Form("EqKaonSqr_%d",sl+1);
//      htitle.Form("Tofino Kaon Peak slab sqrt(e*p) (eq.): %d",sl+1);
//      p_EQkaon_sqr[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("EqE_%d",sl+1);
    htitle.Form("Tofino Equalization Factor slab: %d e",sl+1);
    p_eq_e[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("EqP_%d",sl+1);
    htitle.Form("Tofino Equalization Factor slab: %d p",sl+1);
    p_eq_p[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("EqFactorE_%d",sl+1);
    htitle.Form("Tofino Equalization Factor slab: %d e",sl+1);
    p_EQ_e[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("EqFactorP_%d",sl+1);
    htitle.Form("Tofino Equalization Factor slab: %d p",sl+1);
    p_EQ_p[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);

    hname.Form("OverKaonThreshold_%d",sl+1);
    htitle.Form("Tofino Ratio over kaon threshold slab: %d p",sl+1);
    p_OverKaonThreshold[sl] = new TH1F(hname, htitle, N, run_num_start-0.5, run_num_stop+1.5);
  }
}

void SaveProfiles(const TString &fname){
  TFile *fOutFile = new TFile(fname,"RECREATE");

  for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
    p_ped_e[sl]->Write();
    p_ped_p[sl]->Write();
    p_bha_e[sl]->Write();
    p_bha_p[sl]->Write();
    p_kaon_e[sl]->Write();
    p_kaon_p[sl]->Write(); 
    p_eq_e[sl]->Write();
    p_eq_p[sl]->Write();
    p_EQbha_e[sl]->Write();
    p_EQbha_p[sl]->Write();
//      p_EQbha_sqr[sl]->Write();
    p_EQkaon_e[sl]->Write();
    p_EQkaon_p[sl]->Write();
//      p_EQkaon_sqr[sl]->Write();
    p_EQ_e[sl]->Write();
    p_EQ_p[sl]->Write();
    p_OverKaonThreshold[sl]->Write();
  }

  fOutFile->Close();
  delete fOutFile;
  fOutFile = 0;

}

void CloseHdtFile(){
  
  if (fHdtTree){
    delete fHdtTree;
    // fHdtTree->Delete();
    fHdtTree = 0;
  }
  
  if(fHdtFile){
    if(fHdtFile->IsOpen()) fHdtFile->Close();
    delete fHdtFile;
    fHdtFile = 0;
  }
  
}

Int_t OpenHdtFile(const TString &runtype,const Int_t &run_num){
  // return value:
  //              0: ok
  //             -1: file not found

  CloseHdtFile();

  TString runname;
  TString fnam;
  //  fnam.Form("%s/%s.hdt.root",ExpandPathName("$FND_HDT").Data(),runname.Data());
  //  if (run_num>=7000){
    runname = BuildRunName(runtype,run_num, kTRUE);
    //  fnam.Form("%s/run_group_%d000/%s.hdt.root",ExpandPathName("$FND_HDT").Data(),
    //    Int_t(run_num/1000), runname.Data());
    //} else {
  //runname = BuildRunName(runtype,run_num, kFALSE);
    //fnam.Form("%s_v9.0/%s.hdt.root",ExpandPathName("$FND_HDT").Data(),
    //runname.Data());
  //}
  fnam.Form("%s/%s.hdt.root",ExpandPathName("$FND_HDT").Data(),runname.Data());
  
  if(gSystem->AccessPathName(fnam)){
    gROOT->Warning("OpenHdtFile","file \"%s\" not found",fnam.Data());
    return -1;
  }

  gROOT->Info("OpenHdtFile","opening file \"%s\"",fnam.Data());
  fHdtFile = new TFile(fnam,"OPEN");

  return 0;
}


Int_t LoadHdtTree(){
  // return value:
  //              0: ok
  //             -1: Problems while getting HDT tree from file
  
  if (fHdtTree){
    delete fHdtTree;
    // fHdtTree->Delete();
    fHdtTree = 0;
  }
  
  fHdtTree = (TTree *) ( fHdtFile->Get("FIN_HDT_TREE") );
  if(!fHdtTree){
    gROOT->Warning("LoadHdtTree","Problems while getting HDT tree from file");
    return -1;
  }
  //  fHdtTree->Print();

  fHdtTree->SetBranchAddress("fndhdt",&fCurHdtEv);

  return 0;
}

void GetMultiplicity(Int_t *mult){
  // mult[0] ... low threshold
  // mult[1] ... high threshold
  mult[0]=0; mult[1]=0;

  TClonesArray *TofHits = fCurHdtEv->GetTofHits();
  Int_t lay, slab;
  TFndHTof *cur_hit = 0;

  for(Int_t ic=0;ic< TofHits->GetEntries(); ic++) {
    cur_hit =  (TFndHTof *) ( TofHits->At(ic) );
    cur_hit->GetLayer_Slab(lay, slab);
    if(lay==E_FIN_OUTER_LAYER) continue;

    Bool_t both_adc = kFALSE;
    Bool_t both_tdcl = kFALSE;
    Bool_t both_tdch = kFALSE;
    if(cur_hit->GetBothAdc() == TFndHit::E_HIT_BOTH_CH ) both_adc = kTRUE;
    if(cur_hit->GetBothLowTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdcl = kTRUE;
    if(cur_hit->GetBothHighTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdch = kTRUE;

    Double_t adce=cur_hit->GetAdcE();
    Double_t adcp=cur_hit->GetAdcP();

    if(both_tdch&&adce>0&&adcp>0){ // kaons
      mult[1]++;
    }   
    if(both_tdcl&&adce>0&&adcp>0){ // bhabha+kaons
      mult[0]++;
    }
  }
}
  
Int_t GetTrigType(){
  TFndTrig *cur_trig = (TFndTrig *) (fCurHdtEv->GetGts());
  Short_t hyp_bit=cur_trig->GetPU_cont_bit(0, 0);
  Short_t bha_bit=cur_trig->GetPU_cont_bit(0, 1);
  Short_t cosm_bit=cur_trig->GetPU_cont_bit(0, 5);
  if ((CurrentRunNumber>=5872&&CurrentRunNumber<=5880)||
      (CurrentRunNumber>=5887&&CurrentRunNumber<=5905)||
      (CurrentRunNumber>=5912&&CurrentRunNumber<=5935)||
      (CurrentRunNumber>=5941&&CurrentRunNumber<=5957)||
      (CurrentRunNumber>=5998&&CurrentRunNumber<=6014)||
      (CurrentRunNumber>=6020&&CurrentRunNumber<=6014)||
      (CurrentRunNumber>=7290&&CurrentRunNumber<=7316)||
      (CurrentRunNumber>=7319&&CurrentRunNumber<=7329)||
      (CurrentRunNumber>=7332&&CurrentRunNumber<=7346)||
      (CurrentRunNumber>=7349&&CurrentRunNumber<=7364)||
      (CurrentRunNumber>=7367&&CurrentRunNumber<=7382)||
      (CurrentRunNumber>=7385&&CurrentRunNumber<=7400)||
      (CurrentRunNumber>=7403&&CurrentRunNumber<=7405)||
      (CurrentRunNumber>=7403&&CurrentRunNumber<=7405)||
      (CurrentRunNumber>=7407&&CurrentRunNumber<=7421)||
      (CurrentRunNumber>=7423&&CurrentRunNumber<=7433)||
      (CurrentRunNumber>=7436&&CurrentRunNumber<=7445)||
      (CurrentRunNumber>=7447&&CurrentRunNumber<=7453)||
      (CurrentRunNumber>=7458&&CurrentRunNumber<=7463)||
      (CurrentRunNumber>=7465&&CurrentRunNumber<=7475)||
      (CurrentRunNumber>=7477&&CurrentRunNumber<=7486)||
      (CurrentRunNumber>=7488&&CurrentRunNumber<=7497)||
      (CurrentRunNumber>=7499&&CurrentRunNumber<=7517)||
      (CurrentRunNumber>=7519&&CurrentRunNumber<=7528)||
      (CurrentRunNumber>=7530&&CurrentRunNumber<=7541)||
      (CurrentRunNumber>=7543&&CurrentRunNumber<=7552)||
      (CurrentRunNumber>=7554&&CurrentRunNumber<=7563)||
      (CurrentRunNumber>=7565&&CurrentRunNumber<=7574)||
      (CurrentRunNumber>=7576&&CurrentRunNumber<=7585))      
    { hyp_bit=1; bha_bit=0; cosm_bit=0;}
  if ((CurrentRunNumber>=5881&&CurrentRunNumber<=5886)||
      (CurrentRunNumber>=5906&&CurrentRunNumber<=5911)||
      (CurrentRunNumber>=5936&&CurrentRunNumber<=5940)||
      (CurrentRunNumber>=5969&&CurrentRunNumber<=5973)||
      (CurrentRunNumber>=5986&&CurrentRunNumber<=5991)||
      (CurrentRunNumber>=7289&&CurrentRunNumber<=7289)||
      (CurrentRunNumber>=7317&&CurrentRunNumber<=7318)||
      (CurrentRunNumber>=7330&&CurrentRunNumber<=7331)||
      (CurrentRunNumber>=7347&&CurrentRunNumber<=7348)||
      (CurrentRunNumber>=7365&&CurrentRunNumber<=7366)||
      (CurrentRunNumber>=7383&&CurrentRunNumber<=7384)||
      (CurrentRunNumber>=7401&&CurrentRunNumber<=7402)||
      (CurrentRunNumber>=7406&&CurrentRunNumber<=7406)||
      (CurrentRunNumber>=7422&&CurrentRunNumber<=7422)||
      (CurrentRunNumber>=7434&&CurrentRunNumber<=7435)||
      (CurrentRunNumber>=7446&&CurrentRunNumber<=7446)||
      (CurrentRunNumber>=7454&&CurrentRunNumber<=7457)||
      (CurrentRunNumber>=7464&&CurrentRunNumber<=7464)||
      CurrentRunNumber==7476||CurrentRunNumber==7487||
      CurrentRunNumber==7498||CurrentRunNumber==7518||
      CurrentRunNumber==7529||CurrentRunNumber==7542||
      CurrentRunNumber==7553||CurrentRunNumber==7564||
      CurrentRunNumber==7575||CurrentRunNumber==7586)
    { hyp_bit=0; bha_bit=1; cosm_bit=0;}
  if ((CurrentRunNumber>=7587&&CurrentRunNumber<=8335)||
      (CurrentRunNumber>=9027&&CurrentRunNumber<=9084)||
      (CurrentRunNumber>=9394&&CurrentRunNumber<=9585))
  {
    if (CurrentRunNumber==7597||
        CurrentRunNumber==7608||
        CurrentRunNumber==7619||
        CurrentRunNumber==7630||
        CurrentRunNumber==7639||
        CurrentRunNumber==7650||
        CurrentRunNumber==7661||
        CurrentRunNumber==7670||
        CurrentRunNumber==7690||
        CurrentRunNumber==7691||
        CurrentRunNumber==7702||
        CurrentRunNumber==7713||
        CurrentRunNumber==7719||
        CurrentRunNumber==7733||
        CurrentRunNumber==7744||
        CurrentRunNumber==7756||
        CurrentRunNumber==7777||
        CurrentRunNumber==7788||
        CurrentRunNumber==7792||
        CurrentRunNumber==7793||
        (CurrentRunNumber>=7799&&CurrentRunNumber<=7802)||
        CurrentRunNumber==7814||
        CurrentRunNumber==7825||
        CurrentRunNumber==7836||
        CurrentRunNumber==7847||
        CurrentRunNumber==7858||
        CurrentRunNumber==7869||
        CurrentRunNumber==7873||
        CurrentRunNumber==7874||
        CurrentRunNumber==7885||
        CurrentRunNumber==7896||
        CurrentRunNumber==7907||
        CurrentRunNumber==7919||
        CurrentRunNumber==7921||
        CurrentRunNumber==7932||
        CurrentRunNumber==7943||
        CurrentRunNumber==7954||
        CurrentRunNumber==7965||
        CurrentRunNumber==7976||
        CurrentRunNumber==7987||
        CurrentRunNumber==7991||
        CurrentRunNumber==7992||
        CurrentRunNumber==7993||
        CurrentRunNumber==8004||
        CurrentRunNumber==8005||
        CurrentRunNumber==8023||
        CurrentRunNumber==8064||
        CurrentRunNumber==8065||
        CurrentRunNumber==8074||
        CurrentRunNumber==8085||
        CurrentRunNumber==8098||
        CurrentRunNumber==8109||
        CurrentRunNumber==8120||
        CurrentRunNumber==8131||
        CurrentRunNumber==8142||
        CurrentRunNumber==8153||
        CurrentRunNumber==8163||
        CurrentRunNumber==8174||
        CurrentRunNumber==8185||
        CurrentRunNumber==8208||
        CurrentRunNumber==8222||
        CurrentRunNumber==8233||
        CurrentRunNumber==8244||
        CurrentRunNumber==8255||
        CurrentRunNumber==8266||
        CurrentRunNumber==8277||
        CurrentRunNumber==8288||
        CurrentRunNumber==8299||
        CurrentRunNumber==8303||
        CurrentRunNumber==8308||
        CurrentRunNumber==8319||
        CurrentRunNumber==8334||
        CurrentRunNumber==8335||
        CurrentRunNumber==9029||
        CurrentRunNumber==9040||
        CurrentRunNumber==9051||
        CurrentRunNumber==9062||
        CurrentRunNumber==9073||
        CurrentRunNumber==9401||
        CurrentRunNumber==9412||
        CurrentRunNumber==9416||
        CurrentRunNumber==9428||
        CurrentRunNumber==9439||
        CurrentRunNumber==9450||
        CurrentRunNumber==9461||
        CurrentRunNumber==9472||
        CurrentRunNumber==9483||
        CurrentRunNumber==9494||
        CurrentRunNumber==9505||
        CurrentRunNumber==9516||
        CurrentRunNumber==9527||
        CurrentRunNumber==9538||
        CurrentRunNumber==9549||
        CurrentRunNumber==9560||
        CurrentRunNumber==9571||
        CurrentRunNumber==9582)
      { hyp_bit=0; bha_bit=1; cosm_bit=0;}
    else { hyp_bit=1; bha_bit=0; cosm_bit=0;}
  }
  return (Int_t) hyp_bit+2*bha_bit+32*cosm_bit;
}

void FillHistos(){
  
  TClonesArray *TofHits = fCurHdtEv->GetTofHits();
  
  Int_t lay, slab;
  TFndHTof *cur_hit = 0;
  Int_t mult[2];
  GetMultiplicity(mult);
  if (GetTrigType()!=1||mult[1]!=2) return;

  for(Int_t ic=0;ic< TofHits->GetEntries(); ic++) {
    cur_hit =  (TFndHTof *) ( TofHits->At(ic) );
    //    cur_hit->PrintHit();
    cur_hit->GetLayer_Slab(lay, slab);
    //     TThread::Printf("%d %d %d",lay, cham, wire);
    //     TThread::Printf("");
    if(lay==E_FIN_OUTER_LAYER) continue;

    Bool_t both_adc = kFALSE;
    Bool_t both_tdcl = kFALSE;
    Bool_t both_tdch = kFALSE;
    if(cur_hit->GetBothAdc() == TFndHit::E_HIT_BOTH_CH ) both_adc = kTRUE;
    if(cur_hit->GetBothLowTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdcl = kTRUE;
    if(cur_hit->GetBothHighTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdch = kTRUE;

    //     cout << "a: " << cur_hit->GetBothAdc() << " ; l: " << cur_hit->GetBothLowTdc() << " ; h: " << cur_hit->GetBothHighTdc() << endl;
    //     if(both_tdcl) continue; // pedestal
    //     if(!both_tdcl || both_tdch) continue; // bhabha
    //     if(!both_tdcl) continue; // bhabha + kaons
    //    if(!both_tdch) continue; // kaons

    Double_t adce=cur_hit->GetAdcE();
    Double_t adcp=cur_hit->GetAdcP();
    Double_t sqradc;
    sqradc=sqrt((adce-ped_table_e[slab-1])*eq_table_e[slab-1]
                *(adcp-ped_table_p[slab-1])*eq_table_p[slab-1]);

    if(both_tdch){ // kaons
      if (adce>1&&adcp>1){
        NevAllKaon[slab-1]++;
        if (sqradc>KaonThreshold[slab-1]) NevOverThreshold[slab-1]++;
      }   
    }
  }
}

void EventLoop(const Int_t &nevs){
  
  Int_t entries = fHdtTree->GetEntries();
  
  Int_t cur_ev = 0;
  while(1){
    if(cur_ev >= entries){      
      gROOT->Info("EventLoop","run completed");
      break;
    }
    if(nevs != -1 && cur_ev >= nevs){      
      gROOT->Info("EventLoop","requested number of events reached");
      break;
    }
    //
    fHdtTree->GetEntry(cur_ev);
    //
    FillHistos();
    if(cur_ev%3000 == 0 ) cout << "cur_ev: " << cur_ev << endl;
    //     fCurHdtEv->PrintTofHits();
    //
    cur_ev++;
  }
}

void FillTofFromHdt(const TString &runtype,const Int_t &run_num_start,const Int_t &run_num_stop,const TString &root_filename,const Int_t &nevs, Int_t B){
  
  gStyle->SetOptFit(1);
//   TofCalibration = new TofiCalib();

  // --- loop in runs
  for(Int_t cur_run=run_num_start; cur_run <=run_num_stop ; cur_run++){
    cout << "Filling from run \"" << runtype.Data() << "\" " << cur_run << endl;
    if(OpenHdtFile(runtype,cur_run) !=0) continue; 
    if(LoadHdtTree() != 0) continue;
    
    CurrentRunNumber = cur_run;
    
    for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
     NevOverThreshold[sl]=0;
     NevAllKaon[sl]=0;
    }
    EventLoop(nevs);
    for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
      if (NevAllKaon[sl]>0){
        p_OverKaonThreshold[sl]->Fill(cur_run, (Double_t)NevOverThreshold[sl]/NevAllKaon[sl]);
      }
    }
    CloseHdtFile();
  }

  // --- loop in runs

  // Store information  
  //   TofCalibration->SetPedE(4,998);
  //   cout << "================> " << TofCalibration->GetPedE(4) << endl;
  // --

}

void FillTofFromHdt_TS(const TString &runtype,const Int_t &run_num){
  FillTofFromHdt(runtype, run_num, run_num, "", -1, 1);
}
void StabilityFromPeakTable(const TString &runtype,const Int_t &run_num_start,const Int_t &run_num_stop, const TString &directory, const TString &calib_filename, const TString &root_filename, Bool_t estimate_threshold){
   
  // Threshold For Kaon Selection (very tentative ver.20070126)
  KaonThreshold[0]=350;
  KaonThreshold[1]=400;
  KaonThreshold[2]=400;
  KaonThreshold[3]=450;
  KaonThreshold[4]=470;
  KaonThreshold[5]=590;
  KaonThreshold[6]=500;
  KaonThreshold[7]=500;
  KaonThreshold[8]=500;
  KaonThreshold[9]=300;
  KaonThreshold[10]=360;
  KaonThreshold[11]=400;

  Int_t start_time, stop_time;
  TString temp1, temp2, temp3, temp4, temp5;
  TString temp6, temp7, temp8, temp9, temp10;

  std::fstream cFstrm (calib_filename.Data(), std::ios::in);
  if (cFstrm.fail()) return;

  Int_t line=-1;
  Char_t buff[1000];
  while (!cFstrm.eof()){
    cFstrm.getline(buff, 1000);
    line++;
  }
  cout << "LINE=" << line << endl;
  cFstrm.clear();
  cFstrm.seekg(0, std::ios::beg);

  Int_t start_run[line], stop_run[line];
  Float_t PED_e[line+1][K_N_TOFINO_SLABS], PED_p[line+1][K_N_TOFINO_SLABS];
  Float_t EQ_e[line+1][K_N_TOFINO_SLABS], EQ_p[line+1][K_N_TOFINO_SLABS];
    
    for (Int_t i=0; i<line; i++){ 
      cFstrm >> start_run[i] >> start_time >> temp1 >> temp2 >> temp3 >> temp4 >> temp5 
             >>  stop_run[i] >> stop_time >> temp6 >> temp7 >> temp8 >> temp9 >> temp10;
      cout << start_run[i] << " " << stop_run[i] << endl;
      for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
        cFstrm >> PED_e[i][sl];
      }
      for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
        cFstrm >> PED_p[i][sl];
      }
      for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
        cFstrm >> EQ_e[i][sl];
      }
      for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
        cFstrm >> EQ_p[i][sl];
      }
    }

    for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
      PED_e[line][sl]=0;
      PED_p[line][sl]=0;
      EQ_e[line][sl]=0;
      EQ_p[line][sl]=0;
    }

  DeclareProfiles(run_num_start, run_num_stop);

  Int_t i0=0;

  for (Int_t run_num=run_num_start; run_num<=run_num_stop; run_num++){

    Int_t i;
    // get calibration constants obtained already
    for (i=i0; i<line; i++){
      if (run_num>=start_run[i]&&run_num<=stop_run[i]) break;
    }
      
    for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
      ped_table_e[sl]=PED_e[i][sl];
      ped_table_p[sl]=PED_p[i][sl];
      eq_table_e[sl]=EQ_e[i][sl];
      eq_table_p[sl]=EQ_p[i][sl];
    }

    if (i!=line){
      cout << "Use the calibration table adjusted between run " <<
        start_run[i] << " and " << stop_run[i] << endl;
      i0=i;
      if (estimate_threshold) FillTofFromHdt_TS(runtype, run_num);
    }

    // read the peak table
 
    TString table_filename="./";
    table_filename += directory;
    table_filename += "/";
    table_filename += runtype;
    table_filename += "_";
    table_filename += run_num;
    table_filename += ".txt";
    cout << table_filename << endl;
    std::fstream cFstrm_table (table_filename.Data(), std::ios::in);
    if (cFstrm_table.fail()) continue;

    Double_t ped_e[12], ped_p[12];
    Int_t bhaentry_e[12], bhaentry_p[12];
    Double_t bha_e[12], bha_p[12];
    Int_t kaonentry_e[12], kaonentry_p[12];
    Double_t kaon_e[12], kaon_p[12];
    Double_t eq_e[12], eq_p[12];
    
    for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
      cFstrm_table >> ped_e[sl] >> ped_p[sl]
                   >> bhaentry_e[sl] >> bha_e[sl] 
                   >> bhaentry_p[sl] >> bha_p[sl] 
                   >> kaonentry_e[sl] >> kaon_e[sl] 
                   >> kaonentry_p[sl] >> kaon_p[sl] 
                   >> eq_e[sl] >> eq_p[sl];

      p_ped_e[sl]->Fill(run_num, ped_e[sl]);
      p_ped_p[sl]->Fill(run_num, ped_p[sl]);

      if (ped_e[sl]>0&&bhaentry_e[sl]>500){
        p_bha_e[sl]->Fill(run_num, bha_e[sl]);
        p_EQbha_e[sl]->Fill(run_num, (bha_e[sl]-PED_e[i][sl])*EQ_e[i][sl]);
        p_eq_e[sl]->Fill(run_num, eq_e[sl]);
      }

      if (ped_p[sl]>0&&bhaentry_p[sl]>500){
        p_bha_p[sl]->Fill(run_num, bha_p[sl]);
        p_EQbha_p[sl]->Fill(run_num, (bha_p[sl]-PED_p[i][sl])*EQ_p[i][sl]);
        p_eq_p[sl]->Fill(run_num, eq_p[sl]);
      }

      if (ped_e[sl]>0&&kaon_e[sl]>ped_e[sl]){
        p_kaon_e[sl]->Fill(run_num, kaon_e[sl]);
        p_EQkaon_e[sl]->Fill(run_num, (kaon_e[sl]-PED_e[i][sl])*EQ_e[i][sl]);
      }

      if (ped_p[sl]>0&&kaon_p[sl]>ped_p[sl]){
        p_kaon_p[sl]->Fill(run_num, kaon_p[sl]);
        p_EQkaon_p[sl]->Fill(run_num, (kaon_p[sl]-PED_p[i][sl])*EQ_p[i][sl]);
      }

      p_EQ_e[sl]->Fill(run_num, EQ_e[i][sl]);
      p_EQ_p[sl]->Fill(run_num, EQ_p[i][sl]);
    }
  }
   SaveProfiles(root_filename);
}

void StabilityFromPeakTable(const TString &runtype,const Int_t &run_num_start,const Int_t &run_num_stop, const TString &directory, const TString &calib_filename, const TString &root_filename){
  StabilityFromPeakTable(runtype, run_num_start, run_num_stop, directory,
                         calib_filename, root_filename, false);
}

---