beamtime/cern-oct12/go4/AnalysisMacros/spadic_noise_Analysis.cxx (r4864/r4859)

Go to the documentation of this file.

#include <iostream>
#include <string>
#include <typeinfo>
#include <iostream>
#include <vector>
#include <list>


#include "Riostream.h"
#include "TRint.h"
#include "TROOT.h"
#include "TStyle.h"
#include "Riostream.h"
#include "TRandom2.h"

#include "TH2.h"
#include "TH1.h"
#include "TF1.h"

#include "TGraph.h"
#include "TGraphErrors.h"
#include "TMultiGraph.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TProfile.h"
#include "TStopwatch.h"

#include "TFile.h"
#include "TTree.h"
#include "TBranch.h"
#include "TMath.h"
#include "TLegend.h"
#include "TLine.h"
#include "TColor.h"
#include "TText.h"
#include "TKey.h"
#include "TLatex.h"
#include "TPaveText.h"
#include "TCanvas.h"

#include "TPrincipal.h"
#include "TMatrix.h"

#include "TRocEvent.h"
#include "roc/Message.h"
#include "TSpadicEvent.h"
#include "TEpicsEvent.h"
#include "TFiberHodEvent.h"
#include "TMbsCrateEvent.h"
#include "TBeamMonitorEvent.h"
#include "TCernOct12UnpackEvent.h"
#include "TCernOct12DetectorEvent.h"

#include "PR_tools.cxx"

//#include "CBMsimTR.C"

#define NUM_SPADIC_CHA             8
#define SPADIC_TRACE_SIZE         45
#define SPADIC_ADC_MAX           255
#define noBaselineTB               5 
Double_t fHistoScale = 25;//default to get good calibration fit
//Double_t fHistoScale = 700/3.674/3;//800/3.674;//769.599/3.674;//830.609/3.674;
//Double_t fHistoScale = SPADIC_TRACE_SIZE * SPADIC_ADC_MAX / 100.;// uncalibrated data analysis
Double_t dEdx_min =   0.; //range of the spectrum to calc. pion efficiency full range [0,100]
Double_t dEdx_max = 100.; //range of the spectrum to calc. pion efficiency full range [0,100]
Int_t Pi_Ch1_Min, Pi_Ch1_Max, Pi_Ch2_Min, Pi_Ch2_Max, Pi_Pb_Min, Pi_Pb_Max, El_Ch1_Min, El_Ch1_Max, El_Ch2_Min, El_Ch2_Max, El_Pb_Min, El_Pb_Max;
Bool_t isElectron, isPion, isMyon, isOverFlowEvent, isUnderFlowEvent, isPulser;
TPrincipal* principal = new TPrincipal(NUM_SPADIC_CHA,"ND");
//TMatrixD* cova = new TMatrixD(NUM_SPADIC_CHA,NUM_SPADIC_CHA);
/*
TH1I *clusterSize = new TH1I("clusterSize","clusterSize",NUM_SPADIC_CHA+1,-0.5,NUM_SPADIC_CHA+0.5);
TH1I *maxAmplitudeValue = new TH1I("maxAmplitudeValue","maxAmplitudeValue",256*2,0,256);
TH2I *maxAmplitudeHitTime = new TH2I("maxAmplitudeHitTime","maxAmplitudeHitTime",256,0,256,45,0,45);
TH1I *covaMatixValue = new TH1I("covaMatixValue","covaMatixValue",1000,-0.13,0.13);
TH2I *covaMatixValueMaxAmplitude = new TH2I("covaMatixValueMaxAmplitude","covaMatixValueMaxAmplitude",1000,-0.13,0.13,256,0,256);
TH2I *covaMatixValueHitTime = new TH2I("covaMatixValueHitTime","covaMatixValueHitTime",1000,-0.13,0.13,45,0,45);
TH1I *signalChDistance = new TH1I("signalChDistance","signalChDistance",2*NUM_SPADIC_CHA+1,-NUM_SPADIC_CHA-0.5,NUM_SPADIC_CHA+0.5);
TH2I *averageSignal_2D = new TH2I("averageSignal_2D","averageSignal_2D",SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE,356,-100,256);
TH1I *noiseDistribution = new TH1I("noiseDistribution","noiseDistribution",2*SPADIC_ADC_MAX,-SPADIC_ADC_MAX,SPADIC_ADC_MAX);
TH2I *baselineDistribution = new TH2I("baselineDistribution","baselineDistribution",NUM_SPADIC_CHA,-0.5,NUM_SPADIC_CHA-0.5,SPADIC_ADC_MAX,0,SPADIC_ADC_MAX);
*/
Int_t momentum(0), AnodeV(0), DriftV(0);
Int_t susid = 0;
Int_t pidcuts[18] = {0};
Double_t calibrationParameter[4] = {0.0};

void SetSusid(Int_t sid=0){
  susid = sid;
}

Double_t MPV_kev(0.0), MPV_adc(0.0), sigma_kev(0.0), sigma_adc(0.0), ChiSqr_kev(0.0), ChiSqr_adc(0.0);

void Statusbar(Int_t i, Int_t n) {
 
  if (Int_t(i * 100 / Float_t(n)) - Int_t((i-1) * 100 / Float_t(n)) >= 1){
    //cout << "     " << Int_t(i * 100 / Float_t(n)) << "%" << flush << "\r";
    printf("     %3i%%",Int_t(i * 100 / Float_t(n)));
    cout << flush << "\r";
  }
  //if (int(i * 100 / float(n)) >=99 && int((i-1) * 100 / float(n)) < 99) 
  //cout << endl;
  /*
    if (int(i * 100 / float(n)) - int((i-1) * 100 / float(n)) >= 1) {
    if (int(i * 100 / float(n)) == 1 || i == 1 || i == 0) 
    cout << "[" << flush;
    cout << "-" << flush;
    if (int(i * 10 / float(n)) - int((i-1) * 10 / float(n)) >= 1) 
    cout << "|";
    if (int(i * 100 / float(n)) >=99) 
    cout << "]" <<endl;
    }
  */
}

Double_t scaling2keV(Double_t ADCvalue) {
  if (fHistoScale == SPADIC_TRACE_SIZE * SPADIC_ADC_MAX / 100.){
    return ADCvalue;
  }
  //calibrationParameter[1] = 0.0;
  //return (ADCvalue - calibrationParameter[1]) * calibrationParameter[2]/*3.60898e+00*/ / ((calibrationParameter[0] + calibrationParameter[1])/* * 0.85*/);
  //return (ADCvalue / (calibrationParameter[0] * 0.75)) * calibrationParameter[2]; // based on the maximum with correction factor 2011 standart
  //return (ADCvalue / calibrationParameter[0]) * calibrationParameter[2]; // based on landau MPV
  //return (ADCvalue - calibrationParameter[1]) / (calibrationParameter[0] - calibrationParameter[1]) * calibrationParameter[2]; // based on landau MPV including offset
  //return (ADCvalue) * calibrationParameter[3];// based on 3 points of the slope 
  /*
    if (calibrationParameter[3] < 0.005)
    return ADCvalue * (0.005283 * 1.22);
    else
  */
  //return (MPV_kev - sigma_kev) / (MPV_adc - sigma_adc) * ADCvalue - (MPV_kev - sigma_kev) / (MPV_adc - sigma_adc) * MPV_adc + MPV_kev; //1)

  //return (MPV_kev - sigma_kev) / (MPV_adc - sigma_adc) * ADCvalue - (MPV_kev - sigma_kev) / (MPV_adc - sigma_adc) * MPV_adc + MPV_kev; //1)
  //return (sigma_kev / sigma_adc * ADCvalue + MPV_kev - sigma_kev/ sigma_adc * MPV_adc);// * 1.25; //2)
  return (MPV_kev / MPV_adc * ADCvalue); // 10.07.13 bugfix
  
  //return (ADCvalue - calibrationParameter[1]) / (calibrationParameter[0] - calibrationParameter[1]) * calibrationParameter[2]; // based on landau MPV including offset
  //return (ADCvalue / (calibrationParameter[0]/* * 0.75*/)) * calibrationParameter[2]; // based on the maximum with correction factor 2011 standart
  //return (ADCvalue - calibrationParameter[1]) * (calibrationParameter[3]); // based on 3 points of the slope including offset
}

void CalcPionEfficiency(TH1D* hElectronSpectra, TH1D* hPionSpectra, TGraphErrors* hEfficiency, TGraphErrors* fitEfficiency, const Int_t nChambers, Bool_t debug=false){  
  hEfficiency->SetNameTitle("hEfficiency","hEfficiency");
  fitEfficiency->SetNameTitle("fitEfficiency","fitEfficiency");

  TString cHist, cTitle;
  TH1D *sLikelihoodPi[nChambers+1];
  TH1D *sLikelihoodEl[nChambers+1];
  Int_t nLikeliBins[] = {220,220,220,220,220,220,220,220,220,220,220,220,220};
  for (Int_t Chamberloop = 1; Chamberloop <= nChambers; Chamberloop++){
    cHist.Form("L_Pi_Prob_temp_%i",Chamberloop);
    //sLikelihoodPi[Chamberloop] = new TH1D(cHist,"e and pi Prob",220,-0.05,1.05);
    sLikelihoodPi[Chamberloop] = new TH1D(cHist,"e and pi Prob",nLikeliBins[Chamberloop-1],-0.05,1.05);
    cHist.Form("L_El_Prob_temp_%i",Chamberloop);
    //sLikelihoodEl[Chamberloop] = new TH1D(cHist,"e and pi Prob",220,-0.05,1.05);
    sLikelihoodEl[Chamberloop] = new TH1D(cHist,"e and pi Prob",nLikeliBins[Chamberloop-1],-0.05,1.05);
  }
  Int_t nEvents = 5e7;//1000000;//200000
  TH1D* hElectronProbability = (TH1D*)hElectronSpectra->Clone("ElectronProbability");
  hElectronProbability->Scale(1./hElectronProbability->Integral());
  TH1D* hPionProbability = (TH1D*)hPionSpectra->Clone("hPionProbability");
  hPionProbability->Scale(1./hPionProbability->Integral());
  hEfficiency->SetMarkerStyle(20);
  hEfficiency->GetXaxis()->SetTitle("Number of detector hits");
  hEfficiency->GetYaxis()->SetTitle("Pion efficiency [%]");
  fitEfficiency->SetMarkerStyle(24);
  fitEfficiency->GetXaxis()->SetTitle("Number of detector hits");
  fitEfficiency->GetYaxis()->SetTitle("Pion efficiency [%]");

  for (Int_t Chamberloop = 1; Chamberloop <= nChambers; Chamberloop++){
    cHist.Form("sLikelihoodPi%02d",Chamberloop);
    sLikelihoodPi[Chamberloop]->Reset();
    sLikelihoodPi[Chamberloop]->SetNameTitle(cHist,cHist);
    sLikelihoodPi[Chamberloop]->SetXTitle("Likelihood");
    sLikelihoodPi[Chamberloop]->SetYTitle("#");
    cHist.Form("sLikelihoodEl%02d",Chamberloop);
    sLikelihoodEl[Chamberloop]->Reset();
    sLikelihoodEl[Chamberloop]->SetNameTitle(cHist,cHist);
    sLikelihoodEl[Chamberloop]->SetXTitle("Likelihood");
    sLikelihoodEl[Chamberloop]->SetYTitle("#");
if (hElectronSpectra->Integral() > 0 && hPionSpectra->Integral() > 0){
  Bool_t lowStatistic = false;
  Int_t nElTries(0), nPiTries(0);
  Double_t P_e(1.0), P_pi(1.0);
  Double_t fEl(0.0), fPi(0.0);
  Double_t L_e(0.0), L_p(0.0);
  Double_t chEl(0.0), chPi(0.0);
  for(Int_t iE = 0;iE < nEvents;iE++){
    if (debug) Statusbar(iE,nEvents);
    P_e = 1;
    P_pi = 1;
         
    for(int iCh = 1; iCh <= Chamberloop; iCh++){
      if (lowStatistic)
        break;
      // ladung würfeln
      fEl = 0;
      fPi = 0;
      while (fEl == 0 || fPi == 0) {
        chEl = hElectronSpectra->GetRandom(/*dEdx_min,dEdx_max*//*0.0,100.0*/); 
        //while (hElectronProbability->FindBin(chEl) > dEdx_max) {
        //chEl = hElectronSpectra->GetRandom(/*dEdx_min,dEdx_max*//*0.0,100.0*/); 
        //}     
        fEl = hElectronProbability->GetBinContent( hElectronProbability->FindBin(chEl));
        fPi = hPionProbability->GetBinContent( hPionProbability->FindBin(chEl));
        nElTries++;
        if (nElTries/Float_t(Chamberloop*nEvents) > 15) {
          printf("too low electron statistic\n");
          lowStatistic = true;
          break;
        }
        
      }
      P_e = P_e * (fEl);
      P_pi = P_pi * (fPi);
        
    }      
    L_e = 0;
    if((P_e+P_pi)>0){
      L_e = P_e /(P_e+P_pi);
    }
    sLikelihoodEl[Chamberloop]->Fill(L_e);
      
    P_e = 1;
    P_pi = 1;
          
    L_p = 0;


    for(Int_t iCh = 1;iCh<= Chamberloop;iCh++){
      if (lowStatistic)
        break;
      // ladung würfeln
      fEl = 0;
      fPi = 0;
      while (fEl == 0 || fPi == 0) {
        chPi = hPionSpectra->GetRandom(/*dEdx_min,dEdx_max*//*0.0,100.0*/);
        // while (hElectronProbability->FindBin(chPi) > dEdx_max) {
        //chPi = hElectronSpectra->GetRandom(/*dEdx_min,dEdx_max*//*0.0,100.0*/); 
        //}
        fPi = hPionProbability->GetBinContent( hPionProbability->FindBin(chPi));
        fEl = hElectronProbability->GetBinContent( hElectronProbability->FindBin(chPi));
        nPiTries++;
        if (nPiTries/Float_t(Chamberloop*nEvents) > 15) {
          printf("too low pion statistic\n");
          lowStatistic = true;
          break;
        }
        
      }
      P_e = P_e * (fEl);
      P_pi = P_pi * (fPi);
              
    }    
    L_p = 0; 
    if((P_e+P_pi)>0){
      L_p = 1. - (P_pi /(P_e+P_pi));
    }   
        
    if (lowStatistic){
      printf(" ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! too low statistic\n");
      break;
    }
    else
      sLikelihoodPi[Chamberloop]->Fill(L_p);
  }
  Double_t ElectronIntegral =  0.0;
  Double_t PionIntegral = 0.0;
  Double_t lastPionIntegral = 0.0;
  Int_t ElectronBin = nLikeliBins[Chamberloop-1];         
  Double_t ElectronIntAll = sLikelihoodEl[Chamberloop]->Integral(0,nLikeliBins[Chamberloop-1]);
  Double_t PionIntAll = sLikelihoodPi[Chamberloop]->Integral(0,nLikeliBins[Chamberloop-1]); // misidentified pions of all pions
  if (ElectronIntAll > 0 && PionIntAll > 0) {
    while (ElectronIntegral < 0.90000 && ElectronBin > 0){
      ElectronBin -= 1;
      ElectronIntegral = (sLikelihoodEl[Chamberloop]->Integral(ElectronBin,nLikeliBins[Chamberloop-1]))/ElectronIntAll;
      //PionIntAll = (sLikelihoodPi[Chamberloop]->Integral(ElectronBin,nLikeliBins[Chamberloop-1]) + sLikelihoodEl[Chamberloop]->Integral(ElectronBin,nLikeliBins[Chamberloop-1]));//new remaining pions within 90% electrons
      PionIntegral     = (sLikelihoodPi[Chamberloop]->Integral(ElectronBin,nLikeliBins[Chamberloop-1]))/PionIntAll;
      if (sLikelihoodPi[Chamberloop]->GetBinContent(ElectronBin+1) > 0.0){
        lastPionIntegral = (sLikelihoodPi[Chamberloop]->Integral(ElectronBin+1,nLikeliBins[Chamberloop-1]))/PionIntAll;
      } 
    }
  }
  if (lastPionIntegral == 0.0)
    lastPionIntegral = (sLikelihoodPi[Chamberloop]->Integral(ElectronBin-1,nLikeliBins[Chamberloop-1]))/PionIntAll;
  /*
    if (ElectronBin / nLikeliBins[Chamberloop-1] > 0.955){
    fitLimit--;
    }
  */
  if (debug) {
    printf("  layer %d\n",Chamberloop);
    printf("    Pioneneffizienz     : %8.4f Bin %.3f\n",PionIntegral*100,(ElectronBin)/(Float_t)nLikeliBins[Chamberloop-1]);
    printf("    remaining Pieff.    : %8.4f Bin %.3f\n",PionIntegral*100*PionIntAll/(PionIntegral*PionIntAll+ElectronIntegral*ElectronIntAll),(ElectronBin)/(Float_t)nLikeliBins[Chamberloop-1]);
    printf("    Electroneneffizienz : %8.4f Bin %.3f\n\n",ElectronIntegral*100,(ElectronBin)/(Float_t)nLikeliBins[Chamberloop-1]);
  }
  //hEfficiency->Fill(Chamberloop,PionIntegral*100.);
  hEfficiency->SetPoint(Chamberloop-1, Chamberloop, PionIntegral*100.);
  hEfficiency->SetPointError(Chamberloop-1,0, fabs(PionIntegral - lastPionIntegral) *100.);// Delta pion integral between electron integral < 0.9 and >= 0.9
  
  printf("   %2i %7.4f +- %7.4f  (%.4f EE) (Ebin:%3i, %.4f) %3i bins\n", Chamberloop, PionIntegral*100.,  fabs(PionIntegral - lastPionIntegral) *100., ElectronIntegral, nLikeliBins[Chamberloop-1] - ElectronBin, ElectronBin/Float_t(nLikeliBins[Chamberloop-1]), nLikeliBins[Chamberloop-1]);
 
 }
    else {
      printf("  empty spectra: no pion efficiency calculated\n");
    }
  }
  TF1 *simPionEff = new TF1("simPionEff","exp([0]+[1]*x)",0.5,6.5/*nChambers*/);
  hEfficiency->Fit("simPionEff","RQ0");
  printf("           fit PE:%.6f   (ChiSqr/NDF:%.6f)\n",(exp(simPionEff->GetParameter(0) + simPionEff->GetParameter(1) * nChambers)),simPionEff->GetChisquare() / Float_t(simPionEff->GetNDF()));
  for (Int_t iChamber = 1; iChamber <= nChambers; iChamber++){
    fitEfficiency->SetPoint(iChamber-1, iChamber, (exp(simPionEff->GetParameter(0) + simPionEff->GetParameter(1) * iChamber)));
    fitEfficiency->SetPointError(iChamber-1, 0, fabs((exp(simPionEff->GetParameter(0) + simPionEff->GetParameter(1) * iChamber)) 
                                                     - (exp((simPionEff->GetParameter(0) + simPionEff->GetParError(0)) + (simPionEff->GetParameter(1) + simPionEff->GetParError(1)) * iChamber))));
  }
  delete simPionEff;
  //TFile *outfile = new TFile("Likelihoods.root","RECREATE");
  //outfile->cd();
  for (Int_t Chamberloop = 1; Chamberloop <= nChambers; Chamberloop++){
    //sLikelihoodPi[Chamberloop]->Write(TString(sLikelihoodPi[Chamberloop]->GetName())+TString(hElectronSpectra->GetName()), TObject::kOverwrite);
    //sLikelihoodEl[Chamberloop]->Write(TString(sLikelihoodEl[Chamberloop]->GetName())+TString(hElectronSpectra->GetName()), TObject::kOverwrite);
    delete sLikelihoodPi[Chamberloop];   
    delete sLikelihoodEl[Chamberloop];
  }
  //outfile->Close();
 
}



void CalcPionEfficiency(TH1D* hElectronSpectra, TH1D* hPionSpectra, TProfile* hEfficiency, const Int_t nChambers, Bool_t debug){    
  Bool_t useAntiProbability = false;
  TString cHist, cTitle;
  TH1D *sLikelihoodPi[nChambers+1];
  TH1D *sLikelihoodEl[nChambers+1];
  for (Int_t Chamberloop = 1; Chamberloop <= nChambers; Chamberloop++){
    cHist.Form("L_Pi_Prob_temp_%i",Chamberloop);
    sLikelihoodPi[Chamberloop] = new TH1D(cHist,"e and pi Prob",220,-0.05,1.05);
    cHist.Form("L_El_Prob_temp_%i",Chamberloop);
    sLikelihoodEl[Chamberloop] = new TH1D(cHist,"e and pi Prob",220,-0.05,1.05);
  }
  Int_t nEvents = 5e7;//25000000;//100000;
  TH1D* hElectronProbability = (TH1D*)hElectronSpectra->Clone("ElectronProbability");
  hElectronProbability->Scale(1./hElectronProbability->Integral());
  TH1D* hPionProbability = (TH1D*)hPionSpectra->Clone("hPionProbability");
  hPionProbability->Scale(1./hPionProbability->Integral());
  TCanvas *cTemp2 = NULL;
  TCanvas *cTemp = NULL;
  if (debug) {
    cTemp2 = new TCanvas("temp2","temp2",800,600);
    cTemp2->Divide(1,1);
    cTemp = new TCanvas("temp","temp",800,600);
    cTemp->Divide(2,2);
    cTemp->cd(1);
    hElectronSpectra->DrawCopy();
    cTemp->cd(2);
    hElectronProbability->DrawCopy();
    cTemp->cd(3);
    hPionSpectra->DrawCopy();
    cTemp->cd(4);
    hPionProbability->DrawCopy();
    cTemp->Update();
  }
  hEfficiency->Reset();
  //hEfficiency->SetNameTitle(title,title);
  hEfficiency->GetYaxis()->SetRangeUser(0.01,100.5);
  hEfficiency->SetMarkerStyle(20);
  hEfficiency->SetXTitle("Number of detector layer");
  hEfficiency->SetYTitle("Pion efficiency [%]");

  //if (hElectronSpectra->Integral() > 0 && hPionSpectra->Integral() > 0){
  for (Int_t Chamberloop = 1; Chamberloop <= nChambers; Chamberloop++){
    cHist.Form("sLikelihoodPi%02d",Chamberloop);
    sLikelihoodPi[Chamberloop]->Reset();
    sLikelihoodPi[Chamberloop]->SetNameTitle(cHist,cHist);
    sLikelihoodPi[Chamberloop]->SetXTitle("Likelihood");
    sLikelihoodPi[Chamberloop]->SetYTitle("#");
    cHist.Form("sLikelihoodEl%02d",Chamberloop);
    sLikelihoodEl[Chamberloop]->Reset();
    sLikelihoodEl[Chamberloop]->SetNameTitle(cHist,cHist);
    sLikelihoodEl[Chamberloop]->SetXTitle("Likelihood");
    sLikelihoodEl[Chamberloop]->SetYTitle("#");
    if (hElectronSpectra->Integral() > 0 && hPionSpectra->Integral() > 0){
      Bool_t lowStatistic = false;
      Int_t nElTries(0), nPiTries(0);

      for(Int_t iE = 0;iE < nEvents;iE++){
        if (debug) Statusbar(iE,nEvents);
        Double_t P_e = 1;
        Double_t P_pi = 1;
         
        for(int iCh = 1; iCh <= Chamberloop; iCh++){
          if (lowStatistic)
            break;
          // ladung würfeln
          Double_t fEl = 0;
          Double_t fPi = 0;
          while (fEl == 0 || fPi == 0) {
            Double_t chEl = hElectronSpectra->GetRandom(/*dEdx_min,dEdx_max*//*0.0,100.0*/); 
            //while (hElectronProbability->FindBin(chEl) > dEdx_max) {
            //chEl = hElectronSpectra->GetRandom(/*dEdx_min,dEdx_max*//*0.0,100.0*/); 
            //} 
            fEl = hElectronProbability->GetBinContent( hElectronProbability->FindBin(chEl));
            fPi = hPionProbability->GetBinContent( hPionProbability->FindBin(chEl));
            nElTries++;
            if (nElTries/Float_t(Chamberloop*nEvents) > 15) {
              printf("too low electron statistic\n");
              lowStatistic = true;
              break;
            }
        
          }
          P_e = P_e * (fEl * 1e9);
          P_pi = P_pi * (fPi * 1e9);
        
        }      
        Double_t L_e = 0;
        if((P_e+P_pi)>0){
          L_e = P_e /(P_e+P_pi);
        }
        sLikelihoodEl[Chamberloop]->Fill(L_e);
      
        P_e = 1;
        P_pi = 1;
          
        Double_t L_p = 0;
        if (useAntiProbability) {
          for(Int_t iCh = 1;iCh<= Chamberloop;iCh++){
            if (lowStatistic)
              break;
            // ladung würfeln
            Double_t fEl = 0;
            Double_t fPi = 0;
            while (fEl == 0 || fPi == 0) {
              Double_t chPi = hPionSpectra->GetRandom(/*dEdx_min,dEdx_max*//*0.0,100.0*/);
              //while (hElectronProbability->FindBin(chPi) > dEdx_max) {
              //chPi = hElectronSpectra->GetRandom(/*dEdx_min,dEdx_max*//*0.0,100.0*/); 
              //}
              fEl = hPionProbability->GetBinContent( hPionProbability->FindBin(chPi));
              fPi = hElectronProbability->GetBinContent( hElectronProbability->FindBin(chPi));
              nPiTries++;
              if (nPiTries/Float_t(Chamberloop*nEvents) > 15) {
                printf("too low pion statistic\n");
                lowStatistic = true;
                break;
              }
        
            }
            P_e = P_e * (fEl * 1e9);
            P_pi = P_pi * (fPi * 1e9);
        
          }      
          L_p = 0;
          if((P_e+P_pi)>0){
            L_p = P_pi /(P_e+P_pi);
          }  
        } 
        else {
          for(Int_t iCh = 1;iCh<= Chamberloop;iCh++){
            if (lowStatistic)
              break;
            // ladung würfeln
            Double_t fEl = 0;
            Double_t fPi = 0;
            while (fEl == 0 || fPi == 0) {
              Double_t chPi = hPionSpectra->GetRandom(/*dEdx_min,dEdx_max*//*0.0,100.0*/);
              // while (hElectronProbability->FindBin(chPi) > dEdx_max) {
              //chPi = hElectronSpectra->GetRandom(/*dEdx_min,dEdx_max*//*0.0,100.0*/); 
              //}
              fPi = hPionProbability->GetBinContent( hPionProbability->FindBin(chPi));
              fEl = hElectronProbability->GetBinContent( hElectronProbability->FindBin(chPi));
              nPiTries++;
              if (nPiTries/Float_t(Chamberloop*nEvents) > 15) {
                printf("too low pion statistic\n");
                lowStatistic = true;
                break;
              }
        
            }
            P_e = P_e * (fEl * 1e9);
            P_pi = P_pi * (fPi * 1e9);
              
          }    
          L_p = 0; 
          if((P_e+P_pi)>0){
            L_p = 1. - (P_pi /(P_e+P_pi));
          }   
        }
        if (lowStatistic){
          printf(" ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! too low statistic\n");
          break;
        }
        else
          sLikelihoodPi[Chamberloop]->Fill(L_p);
      }
      Double_t ElectronIntegral =  0.0;
      Double_t PionIntegral = 0.0;
      Int_t ElectronBin = 220;    
      Double_t ElectronIntAll = sLikelihoodEl[Chamberloop]->Integral(0,220);
      Double_t PionIntAll = sLikelihoodPi[Chamberloop]->Integral(0,220);
      if (ElectronIntAll > 0 && PionIntAll > 0) {
        while (ElectronIntegral < 0.90000 && ElectronBin > 0){
          ElectronBin -= 1;
          ElectronIntegral = (sLikelihoodEl[Chamberloop]->Integral(ElectronBin,220))/ElectronIntAll;
          PionIntegral     = (sLikelihoodPi[Chamberloop]->Integral(ElectronBin,220))/PionIntAll;
        }
      }
      if (debug) {
        printf("  layer %d\n",Chamberloop);
        printf("    Pioneneffizienz     : %8.4f Bin %.3f\n",PionIntegral*100,(ElectronBin)/220.);
        printf("    remaining Pieff.    : %8.4f Bin %.3f\n",PionIntegral*100*PionIntAll/(PionIntegral*PionIntAll+ElectronIntegral*ElectronIntAll),(ElectronBin)/220.);
        printf("    Electroneneffizienz : %8.4f Bin %.3f\n\n",ElectronIntegral*100,(ElectronBin)/220.);
      }
      hEfficiency->Fill(Chamberloop,PionIntegral*100.);
      if (debug) {
        cTemp2->cd(1)->SetLogy(1);
        hEfficiency->DrawCopy();
        cTemp2->Update();
      }
    }
    else {
      printf("  empty spectra: no pion efficiency calculated\n");
    }
  }
  //TFile *output = new TFile(file,"UPDATE");
  //output->cd();
  for (Int_t Chamberloop = 1; Chamberloop <= nChambers; Chamberloop++){
    sLikelihoodPi[Chamberloop]->Write("", TObject::kOverwrite);
    sLikelihoodEl[Chamberloop]->Write("", TObject::kOverwrite);
    delete sLikelihoodPi[Chamberloop];   
    delete sLikelihoodEl[Chamberloop];
  }
  // file->cd();
  // hEfficiency->Write("", TObject::kOverwrite);
  if (debug) {
    cTemp->Close();
    delete cTemp;
    cTemp2->Close();
    delete cTemp2;
  }
}


void GetPidCuts(TString runId)
{
  //cout << runId << endl;
  runId.ReplaceAll("Be_run","");
  //cout << runId << endl;
  runId.ReplaceAll("Te_run","");
  //cout << runId << endl;
  runId.ReplaceAll(".root","");
  //cout << runId << endl;
  //int runnumber=runId.Atoi();
  //Int_t pidcuts[18];
  runId.ReplaceAll("data2012/TTrees/merged/","");
  Int_t Momentum = TString(runId(0,1)).Atoi();
  printf("     Momentum %iGeV/c\n",Momentum);
  if (Momentum == 2){
    pidcuts[0] = 1000; // // El Ch1 min
    pidcuts[1] = 4000; // // El Ch1 max
    pidcuts[2] =  800; // // El Ch2 min
    pidcuts[3] = 4000; // // El Ch2 max
    pidcuts[4] =  850;//300; // // El Pb min
    pidcuts[5] = 4000; // // El Pb max
    pidcuts[6] =    0; // // Pi Ch1 min
    pidcuts[7] =  375; // // Pi Ch1 max
    pidcuts[8] =    0; // // Pi Ch2 min
    pidcuts[9] =  350; // // Pi Ch2 max
    pidcuts[10]=  350;//500; // // Pi Pb min
    pidcuts[11]= 4000; // // Pi Pb max
    pidcuts[12]=  375; // // Mu Ch1 min
    pidcuts[13]= 1000; // // Mu Ch1 max
    pidcuts[14]=  350; // // Mu Ch2 min
    pidcuts[15]=  800; // // Mu Ch2 max
    pidcuts[16]=  550; // // Mu Pb min
    pidcuts[17]=  850; // // Mu Pb max
  } else if (Momentum == 3) {
    if (runId == "3GeV_1775V_500V_FFM0.5x100_FFM0.5x200_0" ||
        runId == "3GeV_1775V_500V_FFM0.7x100_FFM0.7x100_0" ||
        runId == "3GeV_1775V_500V_FFM0.7x350_FFM1.2x200_0" ||
        runId == "3GeV_1775V_500V_FFM1.2x100_FFM1.2x100_0" ||
        runId == "3GeV_1775V_500V_H25_H25_0" ||
        runId == "3GeV_1775V_500V_H_FFM0.7x200_0"
        ){
      pidcuts[0] =  800; // // El Ch1 min
      pidcuts[1] = 4000; // // El Ch1 max
      pidcuts[2] =  475; // // El Ch2 min
      pidcuts[3] = 4000; // // El Ch2 max
      pidcuts[4] =  850; // // El Pb min
      pidcuts[5] = 4000; // // El Pb max
      pidcuts[6] =    0; // // Pi Ch1 min
      pidcuts[7] =  375; // // Pi Ch1 max
      pidcuts[8] =    0; // // Pi Ch2 min
      pidcuts[9] =  350; // // Pi Ch2 max
      pidcuts[10]=  350;//500; // // Pi Pb min
      pidcuts[11]= 4000; // // Pi Pb max
      pidcuts[12]=  375; // // Mu Ch1 min
      pidcuts[13]= 1000; // // Mu Ch1 max
      pidcuts[14]=  350; // // Mu Ch2 min
      pidcuts[15]=  600; // // Mu Ch2 max
      pidcuts[16]=  600; // // Mu Pb min
      pidcuts[17]=  850; // // Mu Pb max
    } else if (runId == "3GeV_1800V_500V_B_K_A" ||
               runId == "3GeV_B_K_A"
               ){
      pidcuts[0] =  800; // // El Ch1 min
      pidcuts[1] = 4000; // // El Ch1 max
      pidcuts[2] =  475; // // El Ch2 min
      pidcuts[3] = 4000; // // El Ch2 max
      pidcuts[4] =  850; // // El Pb min
      pidcuts[5] = 4000; // // El Pb max
      pidcuts[6] =    0; // // Pi Ch1 min
      pidcuts[7] =  375; // // Pi Ch1 max
      pidcuts[8] =    0; // // Pi Ch2 min
      pidcuts[9] =  350; // // Pi Ch2 max
      pidcuts[10]=  350;//500; // // Pi Pb min
      pidcuts[11]= 4000; // // Pi Pb max
      pidcuts[12]=  375; // // Mu Ch1 min
      pidcuts[13]= 1000; // // Mu Ch1 max
      pidcuts[14]=  350; // // Mu Ch2 min
      pidcuts[15]=  700; // // Mu Ch2 max
      pidcuts[16]=  350; // // Mu Pb min
      pidcuts[17]=  850; // // Mu Pb max
    } else {
      pidcuts[0] =  650; // // El Ch1 min
      pidcuts[1] = 4000; // // El Ch1 max
      pidcuts[2] =  475; // // El Ch2 min
      pidcuts[3] = 4000; // // El Ch2 max
      pidcuts[4] =  850; // // El Pb min
      pidcuts[5] = 4000; // // El Pb max
      pidcuts[6] =    0; // // Pi Ch1 min
      pidcuts[7] =  375; // // Pi Ch1 max
      pidcuts[8] =    0; // // Pi Ch2 min
      pidcuts[9] =  350; // // Pi Ch2 max
      pidcuts[10]=  350;//500; // // Pi Pb min
      pidcuts[11]= 4000; // // Pi Pb max
      pidcuts[12]=  375; // // Mu Ch1 min
      pidcuts[13]=  650; // // Mu Ch1 max
      pidcuts[14]=  350; // // Mu Ch2 min
      pidcuts[15]=  475; // // Mu Ch2 max
      pidcuts[16]=  600; // // Mu Pb min
      pidcuts[17]=  850; // // Mu Pb max
    }
  } else if (Momentum == 4) {
    pidcuts[0] =  600; // // El Ch1 min
    pidcuts[1] = 4000; // // El Ch1 max
    pidcuts[2] =  400; // // El Ch2 min
    pidcuts[3] = 4000; // // El Ch2 max
    pidcuts[4] =  900; // // El Pb min
    pidcuts[5] = 4000; // // El Pb max
    pidcuts[6] =    0; // // Pi Ch1 min
    pidcuts[7] =  375; // // Pi Ch1 max
    pidcuts[8] =    0; // // Pi Ch2 min
    pidcuts[9] =  350; // // Pi Ch2 max
    pidcuts[10]=  350; // // Pi Pb min
    pidcuts[11]= 4000; // // Pi Pb max
    pidcuts[12]=  375; // // Mu Ch1 min
    pidcuts[13]=  900; // // Mu Ch1 max
    pidcuts[14]=  350; // // Mu Ch2 min
    pidcuts[15]=  550; // // Mu Ch2 max
    pidcuts[16]=  600; // // Mu Pb min
    pidcuts[17]=  900; // // Mu Pb max
  } else if (Momentum == 6) {
    pidcuts[0] =  500; // // El Ch1 min
    pidcuts[1] = 4000; // // El Ch1 max
    pidcuts[2] =  375; // // El Ch2 min
    pidcuts[3] = 4000; // // El Ch2 max
    pidcuts[4] = 1000; // // El Pb min
    pidcuts[5] = 4000; // // El Pb max
    pidcuts[6] =    0; // // Pi Ch1 min
    pidcuts[7] =  375; // // Pi Ch1 max
    pidcuts[8] =    0; // // Pi Ch2 min
    pidcuts[9] =  350; // // Pi Ch2 max
    pidcuts[10]=  350; // // Pi Pb min
    pidcuts[11]= 4000; // // Pi Pb max
    pidcuts[12]=  375; // // Mu Ch1 min
    pidcuts[13]= 1000; // // Mu Ch1 max
    pidcuts[14]=  350; // // Mu Ch2 min
    pidcuts[15]= 1000; // // Mu Ch2 max
    pidcuts[16]=  400; // // Mu Pb min
    pidcuts[17]=  600; // // Mu Pb max
  } else if (Momentum == 8) { 
    pidcuts[0] =  450; // // El Ch1 min
    pidcuts[1] = 4000; // // El Ch1 max
    pidcuts[2] =  375; // // El Ch2 min
    pidcuts[3] = 4000; // // El Ch2 max
    pidcuts[4] =  800; // // El Pb min
    pidcuts[5] = 4000; // // El Pb max
    pidcuts[6] =    0; // // Pi Ch1 min
    pidcuts[7] =  375; // // Pi Ch1 max
    pidcuts[8] =    0; // // Pi Ch2 min
    pidcuts[9] =  350; // // Pi Ch2 max
    pidcuts[10]=  350; // // Pi Pb min
    pidcuts[11]= 4000; // // Pi Pb max
    pidcuts[12]=  375;//400; // // Mu Ch1 min
    pidcuts[13]= 1000;//700; // // Mu Ch1 max
    pidcuts[14]=  350; // // Mu Ch2 min
    pidcuts[15]= 4000;//700;//450; // // Mu Ch2 max
    pidcuts[16]=  425;//500; // // Mu Pb min
    pidcuts[17]=  550;//1000; // // Mu Pb max
  }
}
void getPID(/*TMbsCrateEvent* input, TString runId*/ Float_t ch1, Float_t ch2, Float_t pb)
{
  //int *pidcuts = GetPidCuts(runId);
  
  //cout << "test: " << pidcuts[7] << endl;
  isPion = false;
  isElectron = false;
  isMyon = false;
  if ((ch1 < 350 && pb < 350) || (ch2 < 350 && pb < 350)) return; // noise trigger or event without pid at allgo 
  if (ch1 > pidcuts[0] && // ELECTRON
      ch1 < pidcuts[1] &&
      ch2 > pidcuts[2] &&
      ch2 < pidcuts[3] &&
      pb  > pidcuts[4] &&
      pb  < pidcuts[5]){
    isElectron=true;
  }  
  else if (ch1 > pidcuts[6]  && // PION
           ch1 < pidcuts[7]  &&
           ch2 > pidcuts[8]  &&  // was 340
           ch2 < pidcuts[9]  &&  // was 348
           pb  > pidcuts[10] &&
           pb  < pidcuts[11]){
    isPion=true;
  }     
  else if (((ch1 > pidcuts[12] && ch1 < pidcuts[13]) 
            ||
            (ch2 > pidcuts[14] && ch2 < pidcuts[15]))
           &&
           (pb  > pidcuts[16] && pb  < pidcuts[17])){
    isMyon=true;
  }  
}

Int_t GetPadMax(TH1D* inPulse[NUM_SPADIC_CHA]) {
  Int_t maxCh = -1;
  Double_t maxIntegral = 0.0;
  Double_t tempIntegral = 0.0;
  //Double_t Integral[NUM_SPADIC_CHA] = {0.0};
  for(Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++){
    tempIntegral = 0.0;
    for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {
      if (inPulse[ch]->GetBinContent(inPulse[ch]->FindBin(bin)) > 0.0)
        tempIntegral += inPulse[ch]->GetBinContent(inPulse[ch]->FindBin(bin));
    }
    //Integral[ch] = tempIntegral;
    if (tempIntegral > maxIntegral) {
      maxCh = ch;
      maxIntegral = tempIntegral;
    }
  }

  /*
    Double_t maximum = 0.0;
    Int_t maxCh = -1;
    for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++){
    if (inPulse[ch]->GetBincontent(inPulse[ch]->GetMaximumBin()) > maximum){
    maximum = inPulse[ch]->GetBincontent(inPulse[ch]->GetMaximumBin());
    maxCh = ch;
    }
    }
  */
  return maxCh;
}

Double_t Get3PadIntegral(TH1D* inPulse[NUM_SPADIC_CHA])
{
  Int_t maxCh = GetPadMax(inPulse);
  Double_t signal = inPulse[maxCh]->Integral();
  if (maxCh-1 >= 0)
    signal += inPulse[maxCh-1]->Integral();
  if (maxCh+1 < NUM_SPADIC_CHA)
    signal += inPulse[maxCh+1]->Integral();
  return signal;
}
Double_t Get3PadPosAmplitude(TH1D* inPulse[NUM_SPADIC_CHA])
{
  Int_t maxCh = GetPadMax(inPulse);
  Int_t maxBin = inPulse[maxCh]->GetMaximumBin();
  Double_t signal = 0.0;
  if (inPulse[maxCh]->GetBinContent(maxBin) > 0)
    signal = inPulse[maxCh]->GetBinContent(maxBin);
  if (maxCh-1 >= 0 && inPulse[maxCh-1]->GetBinContent(maxBin) > 0)
    signal += inPulse[maxCh-1]->GetBinContent(maxBin);
  if (maxCh+1 < NUM_SPADIC_CHA && inPulse[maxCh+1]->GetBinContent(maxBin) > 0)
    signal += inPulse[maxCh+1]->GetBinContent(maxBin);
  return signal;
}
Double_t Get3PadPosIntegral(TH1D* inPulse[NUM_SPADIC_CHA])
{
  const Int_t minBin = 0;
  const Int_t maxBin = SPADIC_TRACE_SIZE-1;
  Int_t maxCh = -1;
  Double_t maxIntegral = 0.0;
  Double_t tempIntegral = 0.0;
  Double_t Integral[NUM_SPADIC_CHA] = {0.0};
  for(Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++){
    tempIntegral = 0.0;
    for (Int_t bin = minBin; bin < maxBin; bin++) {
      if (inPulse[ch]->GetBinContent(inPulse[ch]->FindBin(bin)) > 0.0)
        tempIntegral += inPulse[ch]->GetBinContent(inPulse[ch]->FindBin(bin));
    }
    Integral[ch] = tempIntegral;
    if (tempIntegral > maxIntegral) {
      maxCh = ch;
      maxIntegral = tempIntegral;
    }
  }
  Double_t signal = Integral[maxCh];
  if (maxCh-1 >= 0)
    signal += Integral[maxCh-1];
  if (maxCh+1 < NUM_SPADIC_CHA)
    signal += Integral[maxCh+1];
  return signal;
}
Double_t Get3PadPosIntegralCrossTalk(TH1D* inPulse1[NUM_SPADIC_CHA], TH1D* inPulse2[NUM_SPADIC_CHA], Int_t *deltaChMax)
{
  const Int_t minBin = 0;
  const Int_t maxBin = SPADIC_TRACE_SIZE-1;
  Int_t maxCh1(-1), maxCh2(-1);
  Double_t maxIntegral1(0.0), maxIntegral2(0.0);
  Double_t tempIntegral1(0.0), tempIntegral2(0.0);
  Double_t Integral1[NUM_SPADIC_CHA] = {0.0};
  Double_t Integral2[NUM_SPADIC_CHA] = {0.0};
  for(Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++){
    tempIntegral1 = 0.0;
    tempIntegral2 = 0.0;
    for (Int_t bin = minBin; bin < maxBin; bin++) {
      if (inPulse1[ch]->GetBinContent(inPulse1[ch]->FindBin(bin)) > 0.0)
        tempIntegral1 += inPulse1[ch]->GetBinContent(inPulse1[ch]->FindBin(bin));
      if (inPulse2[ch]->GetBinContent(inPulse2[ch]->FindBin(bin)) > 0.0)
        tempIntegral2 += inPulse2[ch]->GetBinContent(inPulse2[ch]->FindBin(bin));
    }
    Integral1[ch] = tempIntegral1;
    Integral2[ch] = tempIntegral2;
    if (tempIntegral1 > maxIntegral1) {
      maxCh1 = ch;
      maxIntegral1 = tempIntegral1;
    }
    if (tempIntegral2 > maxIntegral2) {
      maxCh2 = ch;
      maxIntegral2 = tempIntegral2;
    }
  }
  *deltaChMax = fabs(maxCh1 - maxCh2) ;
   for(Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++){
       Integral1[ch] += Integral2[ch];
   }
   Double_t signal = Integral1[maxCh1];
   if (maxCh1-1 >= 0)
     signal += Integral1[maxCh1-1];
   if (maxCh1+1 < NUM_SPADIC_CHA)
     signal += Integral1[maxCh1+1];
   return signal;

 }
/*
  Double_t GetPosition(TH1D* inPulse[NUM_SPADIC_CHA], Double_t chargePercentage[3]) 
  {
  const Int_t minBin = 0;
  const Int_t maxBin = SPADIC_TRACE_SIZE-1;
  Double_t dxPos = 0;
  Int_t maxCh = GetPadMax(inPulse);
  Double_t intensCh[NUM_SPADIC_CHA] = {0.0};
  if (maxCh >= 1 && maxCh <= NUM_SPADIC_CHA - 1){
  for (Int_t ch = maxCh-1; ch <= maxCh+1; ch++){
  for (Int_t bin = minBin; bin < maxBin; bin++) {
  if (inPulse[ch]->GetBinContent(inPulse[ch]->FindBin(bin)) > 0.0)
  intensCh[ch] += inPulse[ch]->GetBinContent(inPulse[ch]->FindBin(bin));
  }
  }
  chargePercentage[0] = intensCh[maxCh-1] / (intensCh[maxCh-1] + intensCh[maxCh] + intensCh[maxCh+1]);
  chargePercentage[1] = intensCh[maxCh]   / (intensCh[maxCh-1] + intensCh[maxCh] + intensCh[maxCh+1]);
  chargePercentage[2] = intensCh[maxCh+1] / (intensCh[maxCh-1] + intensCh[maxCh] + intensCh[maxCh+1]);
  dxPos = 0.5 * log((intensCh[maxCh+1] / intensCh[maxCh-1])) / log((pow(intensCh[maxCh],2) / (intensCh[maxCh+1] * intensCh[maxCh-1])));
  }
  return dxPos;
  }
*/
Double_t calcK3at0hs(Double_t ra=0.02/*[mm]*/, Double_t s=2.5/*[mm]*/){
  // valid within ra/s [1.5E-3,7.5E-3]
  return -20.19*(ra/s)+0.85; //lin fit
  //return 0.68*TMath::Exp(-26.70*(ra/s)+0.22);// exp fit
  //return 0.86-25.34*(ra/s)+569.86*(ra/s)*(ra/s);// poly 2 fit
}

Double_t calcK3(Float_t h=4.0){
  // valid within h/s [0.5,1.6]
  Double_t s=2.5;//[mm]
  Double_t ra=0.02;//[mm]
  return -0.245 * (h/s) + calcK3at0hs(ra,s);
  //return -0.245 * (h/2.5) + 0.70; //PhD status 08.05.2013
}
void PrfAna(TH1D* inPulse[NUM_SPADIC_CHA], TH2D* clusterwCharge, Int_t maxCh = -1, Int_t clusterWidth = 0){
  const Double_t W = 7.125;//[mm]


  const Double_t par = 1.0;
  const Double_t h = 3.5; //[mm]
  const Double_t K3 = calcK3(h);//0.525; 
  const Double_t SqrtK3 = sqrt(K3);
  const Int_t minBin = 0;
  const Int_t maxBin = SPADIC_TRACE_SIZE-1;
  if(maxCh == -1) maxCh = GetPadMax(inPulse);
  Double_t intensCh[NUM_SPADIC_CHA] = {0.0};
  
  //Int_t clusterWidth = 0;
  if (maxCh >= 1 && maxCh <= NUM_SPADIC_CHA - 1){
    for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++){
      for (Int_t bin = minBin; bin < maxBin; bin++) {
        if (inPulse[ch]->GetBinContent(inPulse[ch]->FindBin(bin)) > 0.0)
          intensCh[ch] += inPulse[ch]->GetBinContent(inPulse[ch]->FindBin(bin));
      }
    }
    /*
    for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++){
      if(intensCh[ch] >= 0.1 * intensCh[maxCh]){
        clusterWidth++;
      }
    }
    */
    Double_t d = 0.5 * W * log((intensCh[maxCh+1] / intensCh[maxCh-1])) / log((pow(intensCh[maxCh],2) / (intensCh[maxCh+1] * intensCh[maxCh-1])));
    Double_t chargeTheo = -1. / (2. * atan(SqrtK3)) * (atan(SqrtK3 * tanh(TMath::Pi() * (-2. + SqrtK3 ) * (W + 2.* d * par) / (8.* h) )) +  atan(SqrtK3 *  tanh(TMath::Pi() * (-2. + SqrtK3 ) * (W - 2.* d * par) / (8.* h) )) );       
    Double_t chargePercent = intensCh[maxCh] / (intensCh[maxCh-1] + intensCh[maxCh] + intensCh[maxCh+1]);
    clusterwCharge->Fill(clusterWidth, chargeTheo - chargePercent);
  }

}

void CalcPrf(TH1D* inPulse[NUM_SPADIC_CHA], TH2D* prf, TProfile* prfProfile, Double_t padWidth, Int_t maxCh = -1)
{
  const Int_t minBin = 0;
  const Int_t maxBin = SPADIC_TRACE_SIZE-1;
  if(maxCh == -1) maxCh = GetPadMax(inPulse);
  //Double_t intensCh[NUM_SPADIC_CHA] = {0.0};
  Double_t intensCh[3] = {0.0};
  if (maxCh >= 1 && maxCh <= NUM_SPADIC_CHA - 1){
    //for (Int_t ch = maxCh-1; ch <= maxCh+1; ch++){
    for (Int_t ch = 0; ch < 3; ch++){
      for (Int_t bin = minBin; bin < maxBin; bin++) {
        //if (inPulse[ch]->GetBinContent(inPulse[ch]->FindBin(bin)) > 0.0)
        if (inPulse[maxCh-1+ch]->GetBinContent(inPulse[maxCh-1+ch]->FindBin(bin)) > 0.0)
          //intensCh[ch] += inPulse[ch]->GetBinContent(inPulse[ch]->FindBin(bin));
          intensCh[ch] += inPulse[maxCh-1+ch]->GetBinContent(inPulse[maxCh-1+ch]->FindBin(bin));
      }
    }
    //Double_t dxPos = 0.5 * log((intensCh[maxCh+1] / intensCh[maxCh-1])) / log((pow(intensCh[maxCh],2) / (intensCh[maxCh+1] * intensCh[maxCh-1])));
    Double_t dxPos = calcSimpleDisplacement( padWidth, intensCh);//0.5 * padWidth * log((intensCh[2] / intensCh[0])) / log((pow(intensCh[1],2) / (intensCh[2] * intensCh[0])));
    //Double_t dxPos = calcDisplacementHyperbolicSecantSquared( padWidth, intensCh);
    //Double_t dxPos = calcDisplacementPRFFitting( padWidth, intensCh);

    Double_t QSum = (intensCh[0] + intensCh[1] + intensCh[2]);

    //Double_t chargePercent = intensCh[maxCh-1] / (intensCh[maxCh-1] + intensCh[maxCh] + intensCh[maxCh+1]);
    Double_t chargePercent = intensCh[0] / QSum;
    prf->Fill(-dxPos-padWidth, chargePercent);
    prfProfile->Fill(-dxPos-padWidth ,chargePercent);
    //chargePercent = intensCh[maxCh] / (intensCh[maxCh-1] + intensCh[maxCh] + intensCh[maxCh+1]);
    chargePercent = intensCh[1] / QSum;
    prf->Fill(dxPos,chargePercent);
    prfProfile->Fill(dxPos ,chargePercent);
    //chargePercent = intensCh[maxCh+1] / (intensCh[maxCh-1] + intensCh[maxCh] + intensCh[maxCh+1]);
    chargePercent = intensCh[2] / QSum;
    prf->Fill(-dxPos+padWidth,chargePercent);
    prfProfile->Fill(-dxPos+padWidth ,chargePercent);
  }
}
void FillAverageSignal(TH1D* inPulse[NUM_SPADIC_CHA], TProfile* averagePulse, TH2I* averagePulse2D)
{
  const Int_t minBin = 0;
  const Int_t maxBin = SPADIC_TRACE_SIZE-1;
  Int_t maxCh = GetPadMax(inPulse);
  for (Int_t bin = minBin; bin < maxBin; bin++) {
    averagePulse->Fill(bin,inPulse[maxCh]->GetBinContent(inPulse[maxCh]->FindBin(bin)));
    averagePulse2D->Fill(bin,inPulse[maxCh]->GetBinContent(inPulse[maxCh]->FindBin(bin)));
  }
}

  /*
    void get_Pb_Ch_Cut(TString fname) {
    ifstream cutLib;
    cutLib.open("cutConditions.txt");
    TString run;
    string line;
    Bool_t found = false;
    if (cutLib.is_open()) {
    getline (cutLib, line);
    while (run != fname && cutLib.good()){
    cutLib >> run >> Pi_Ch1_Min >> Pi_Ch1_Max >> Pi_Ch2_Min >> Pi_Ch2_Max >> Pi_Pb_Min >> Pi_Pb_Max >> El_Ch1_Min >> El_Ch1_Max >> El_Ch2_Min >> El_Ch2_Max >> El_Pb_Min >> El_Pb_Max;
    //printf("%s\n",run.Data());
    if(run == fname){
    printf("          run                  Pi_Ch1_Min, Pi_Ch1_Max, Pi_Ch2_Min, Pi_Ch2_Max,  Pi_Pb_Min,  Pi_Pb_Max, El_Ch1_Min, El_Ch1_Max, El_Ch2_Min, El_Ch2_Max,  El_Pb_Min,  El_Pb_Max\n          %s %11i %11i %11i %11i %11i %11i %11i %11i %11i %11i %11i %11i\n",run.Data(), Pi_Ch1_Min, Pi_Ch1_Max, Pi_Ch2_Min, Pi_Ch2_Max, Pi_Pb_Min, Pi_Pb_Max, El_Ch1_Min, El_Ch1_Max, El_Ch2_Min, El_Ch2_Max, El_Pb_Min, El_Pb_Max);
    found = true;
    }
    }
    if (!found) printf("did not found pid cut conditions for %s\n",fname.Data());
    }
    cutLib.close();
    }


    void getPID(Float_t Ch1, Float_t Ch2, Float_t Pb) {
    if (Ch1 < Pi_Ch1_Max && Ch1 > Pi_Ch1_Min)
    if (Ch2 < Pi_Ch2_Max && Ch2 > Pi_Ch2_Min)
    if (Pb < Pi_Pb_Max && Pb > Pi_Pb_Min)
    isPion = true;
    if (Ch1 < El_Ch1_Max && Ch1 > El_Ch1_Min)
    if (Ch2 < El_Ch2_Max && Ch2 > El_Ch2_Min)
    if (Pb < El_Pb_Max && Pb > El_Pb_Min)
    isElectron = true;
    }
  */

Double_t GetMaxAplitude(TH1D* inPulse[NUM_SPADIC_CHA], Int_t maxCh = -1) {
  if (maxCh == -1) maxCh = GetPadMax(inPulse);
  return inPulse[maxCh]->GetBinContent(inPulse[maxCh]->GetMaximumBin());
}


  Bool_t HitTimeWindow(TH1D* inPulse[NUM_SPADIC_CHA], Int_t maxCh = -1) {
    if (maxCh == -1) maxCh = GetPadMax(inPulse);
    if (inPulse[maxCh]->GetMaximumBin() > 10 && inPulse[maxCh]->GetMaximumBin() < 20) // 13, 20
      return true;

    return false;
  }

  Bool_t BorderPadMax(TH1D* inPulse[NUM_SPADIC_CHA], Int_t maxCh = -1) {
    if (maxCh == -1) maxCh = GetPadMax(inPulse);
    if (maxCh == 0 || maxCh == NUM_SPADIC_CHA-1)
      return true;
    return false;
  }
Bool_t minimumIntegralThreshold(TH1D* inPulse[NUM_SPADIC_CHA]) {
  if (Get3PadPosIntegral(inPulse) > 25
      //25 
      //50
      )//2011: 25
    return true;
  return false;
}

Bool_t minimumAmplitudeThreshold(TH1D* inPulse[NUM_SPADIC_CHA], Int_t maxCh = -1) {
  /*
    Int_t maxCh = GetPadMax(inPulse);
    if (maxCh > -1){
    if (inPulse[maxCh]->GetBinContent(inPulse[maxCh]->GetMaximumBin()) > 15){
    return true;
    }
    }
  */
  /*
    for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++)
    maxAmplitudeValue->Fill(inPulse[ch]->GetBinContent(inPulse[ch]->GetMaximumBin()));
  */
  Double_t minAmpl[3] = {12,15,15};//{10,15,15};
  if (maxCh == -1){ 
    for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++){
      if (inPulse[ch]->GetBinContent(inPulse[ch]->GetMaximumBin()) > minAmpl[susid]){ //2011: 15
        return true;
      }    
    }
    return false;
  }
  else {
    if (inPulse[maxCh]->GetBinContent(inPulse[maxCh]->GetMaximumBin()) > minAmpl[susid]){
      return true;
    }
    return false;
  }
  cout << ">>>>>>>>>>>>>>>>  BIG SHIT HAS HAPPEND!!!!!!!!!!!!!!" << endl;
  return false;
}

  void BaselineCompensation(TH1D* inPulse[NUM_SPADIC_CHA], TH1D* outPulse[NUM_SPADIC_CHA], TH2I* baselineDistribution=NULL)
  {
    Double_t baseline[NUM_SPADIC_CHA] = {0.0};
    for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++){
      outPulse[ch]->Reset();
      for (Int_t bin = 0; bin < noBaselineTB; bin++){
        baseline[ch] += inPulse[ch]->GetBinContent(inPulse[ch]->FindBin(bin));
      }
      baseline[ch] /= (Double_t)noBaselineTB;
      if (baseline[ch] > 1)
        if(baselineDistribution)
          baselineDistribution->Fill(ch,baseline[ch]);
    }
   
    for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++){
      for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++){
        outPulse[ch]->Fill(bin, inPulse[ch]->GetBinContent(inPulse[ch]->FindBin(bin)) - baseline[ch]);
      }
    }
  }
  /*
    void NoiseCompensation(TH1D* inPulse[NUM_SPADIC_CHA],TH1D* outPulse[NUM_SPADIC_CHA])
    {  
    Double_t minIntegral = 265*45;
    Double_t tempIntegral = 0.0;
    int minch = -1;//GetMinimumChannel(input);
    Int_t max2minCh[NUM_SPADIC_CHA] = {-1};
    Double_t Integral[NUM_SPADIC_CHA] = {0.0};
    for(Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++){
    tempIntegral = 0.0;   
    Integral[ch] = inPulse[ch]->Integral();
    if (tempIntegral < minIntegral){
    minIntegral = tempIntegral;
    minch = ch;
    }
    }
    }
  */
Int_t GetClusterSize(){
  return -1;
}

Int_t CancelNoise_Cova(TH1D* inPulse[NUM_SPADIC_CHA],
                       TH1D* outPulse[NUM_SPADIC_CHA], 
                       TH1I*covaMatixValue=NULL,
                       TH2I*covaMatixValueMaxAmplitude=NULL,
                       TH2I*covaMatixValueHitTime=NULL,
                       TH2I*maxAmplitudeHitTime=NULL,
                       TH1I*noiseDistribution=NULL,
                       TH1I*clusterSize=NULL,
                       TH1I*signalChDistance=NULL,
                       TH2I*averageNoise_2D=NULL,
                       TH2I*averageSignal_2D=NULL,
                       TH2I*covaMatixValueClusterSize=NULL,
                       Bool_t debug=false){
  Bool_t rawStatistic = false;
  Double_t lowercorthreshold[3] = {0.120, 0.112, 0.111};//{0.121, 0.112, 0.111};
  Double_t minAmplitude[3] = {30, 16, 20};
  //const double uppercorthreshold = 0.112;
  //const double nosignalcorthreshold = 0.1;
  Double_t cortest[NUM_SPADIC_CHA];
  Double_t inputarray[SPADIC_TRACE_SIZE][NUM_SPADIC_CHA] = {{0.0}};
  Double_t minIntegral = 265*45;
  Double_t tempIntegral = 0.0;
  Int_t minCh = -1;//GetMinimumChannel(input);
  Int_t maxCh = -1;//GetPadMax(inPulse);
  Int_t noise_ch_counter=0;
  Int_t signal_ch_counter=0;

  for(Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++){
    tempIntegral = 0.0;
    outPulse[ch]->Reset();
    for(Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++){  
      inputarray[bin][ch] = inPulse[ch]->GetBinContent(inPulse[ch]->FindBin(bin));
      tempIntegral += inputarray[bin][ch];
    }
    if (tempIntegral < minIntegral){
      minIntegral = tempIntegral;
      minCh = ch;
    }
  }

  Double_t thisisnoise[SPADIC_TRACE_SIZE] = {0.0};

  for(int t=0;t<SPADIC_TRACE_SIZE;t++){
    principal->AddRow(inputarray[t]);
  }
  
  principal->MakePrincipals();
  TMatrixD* cova = (TMatrixD*)principal->GetCovarianceMatrix();
  
  Int_t lastSigCh = -1;
  std::list<Double_t> covaList;
  std::list<Double_t>::iterator it;
  for(int ch=0;ch<NUM_SPADIC_CHA;ch++){
    cortest[ch] = TMatrixDRow (*cova,minCh)(ch);
    if(covaMatixValueClusterSize){
      covaList.push_back(cortest[ch]);
    }
    if((cortest[ch] > lowercorthreshold[susid]) 
       && (
           (inPulse[ch]->GetBinContent(inPulse[ch]->GetMaximumBin()) < minAmplitude[susid])
           || ((inPulse[ch]->GetMaximumBin() < 5 || inPulse[ch]->GetMaximumBin() > 20))
           )
       // && (fabs(inPulse[ch]->GetBinContent(inPulse[ch]->GetMinimumBin())) < minAmplitude[susid])
       ){
      noise_ch_counter++;
      for(Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++){
        thisisnoise[bin] += inputarray[bin][ch];
      }
    }
    else { //debug
      if (lastSigCh != -1){
        if(signalChDistance){
          signalChDistance->Fill(ch - lastSigCh/*maxCh-ch*/); //debug
        }
      }
      lastSigCh = ch;
    } //debug
  }
  if(covaMatixValueClusterSize) {
    covaList.sort();
    Int_t counter = 0;//NUM_SPADIC_CHA-1;
    for (it = covaList.begin(); it != covaList.end(); it++){
      if(debug) cout << *it << " (" << counter << ") | ";
      covaMatixValueClusterSize->Fill(*it,counter);
      counter++;
    }
    if(debug) cout << endl;
    covaList.clear();
  }
  //if (noise_ch_counter > 0){ // have found a noise channel to subtract
  //if (noise_ch_counter < NUM_SPADIC_CHA && noise_ch_counter > 0){
  for(Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++){
    thisisnoise[bin] /= (double)noise_ch_counter;
    if (noise_ch_counter < NUM_SPADIC_CHA && noise_ch_counter > 0){
      if(noiseDistribution){
        noiseDistribution->Fill(thisisnoise[bin]);
      }
    }    
  } /*//else {cout << "very bad thing happened!!" << endl;}
    //} else 
  if (noise_ch_counter == 0){ // if no noise channels could be found just subtract the lowest channel to avoid underflows
    for(Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++){
      thisisnoise[bin] = inputarray[bin][minCh];
    }
  }
    */

  for (int ch=0; ch<NUM_SPADIC_CHA; ch++){     
    for (int k=0; k<SPADIC_TRACE_SIZE; k++){  
      outPulse[ch]->Fill(k, inPulse[ch]->GetBinContent(inPulse[ch]->FindBin(k)) - thisisnoise[k]);
      if ( TString(inPulse[0]->GetTitle()).BeginsWith("baselinePulse")){
        if(averageSignal_2D && averageNoise_2D){
          if (rawStatistic){
            averageSignal_2D->Fill(k,inPulse[ch]->GetBinContent(inPulse[ch]->FindBin(k)));
          } else {
            averageSignal_2D->Fill(k,outPulse[ch]->GetBinContent(outPulse[ch]->FindBin(k)));
          }
          averageNoise_2D->Fill(k,thisisnoise[k]);
        }
      }
    }
    if ( TString(inPulse[ch]->GetTitle()).BeginsWith("baselinePulse")){
      if (rawStatistic){
        if (ch != minCh)
          if (covaMatixValue)
            covaMatixValue->Fill(cortest[ch]);
        if (ch != minCh)
          if(covaMatixValueMaxAmplitude)
            covaMatixValueMaxAmplitude->Fill(cortest[ch], inPulse[ch]->GetBinContent(inPulse[ch]->GetMaximumBin()));
        if (ch != minCh)
          if(covaMatixValueHitTime)
            covaMatixValueHitTime->Fill(cortest[ch], inPulse[ch]->GetMaximumBin());
        if(maxAmplitudeHitTime)
          if (inPulse[ch]->GetBinContent(inPulse[ch]->GetMaximumBin()) > 0.0){
            maxAmplitudeHitTime->Fill(inPulse[ch]->GetBinContent(inPulse[ch]->GetMaximumBin()), inPulse[ch]->GetMaximumBin());
          }
      } else {
        if (ch != minCh){
          if (covaMatixValue){
            covaMatixValue->Fill(cortest[ch]);
          }     
          if(covaMatixValueMaxAmplitude){
            covaMatixValueMaxAmplitude->Fill(cortest[ch], outPulse[ch]->GetBinContent(outPulse[ch]->GetMaximumBin()));
          }     
          if(covaMatixValueHitTime){
            covaMatixValueHitTime->Fill(cortest[ch], outPulse[ch]->GetMaximumBin());
          }
        }
        if(maxAmplitudeHitTime){
          if (outPulse[ch]->GetBinContent(outPulse[ch]->GetMaximumBin()) > 0.0){
            maxAmplitudeHitTime->Fill(outPulse[ch]->GetBinContent(outPulse[ch]->GetMaximumBin()), outPulse[ch]->GetMaximumBin());
          }
        }       
      }
    }
  }
  if ( TString(inPulse[0]->GetTitle()).BeginsWith("baselinePulse")){
    if (clusterSize){      
      if (rawStatistic){
        maxCh = GetPadMax(inPulse);
        for (Int_t iCh = maxCh; iCh < NUM_SPADIC_CHA; iCh++){
          if (
              (cortest[iCh] > lowercorthreshold[susid]) 
              && (
                  (inPulse[iCh]->GetBinContent(inPulse[iCh]->GetMaximumBin()) < minAmplitude[susid])
                  || ((inPulse[iCh]->GetMaximumBin() < 5 || inPulse[iCh]->GetMaximumBin() > 20))        
                  )             
              ) break;
          signal_ch_counter++;
        }
        for (Int_t iCh = maxCh-1; iCh >= 0; iCh--){
          if (
              (cortest[iCh] > lowercorthreshold[susid]) 
              && (
                  (inPulse[iCh]->GetBinContent(inPulse[iCh]->GetMaximumBin()) < minAmplitude[susid])
                  || ((inPulse[iCh]->GetMaximumBin() < 5 || inPulse[iCh]->GetMaximumBin() > 20))
                  )             
              ) break;
          signal_ch_counter++;
        } 
      } else {
        maxCh = GetPadMax(outPulse);
        for (Int_t iCh = maxCh; iCh < NUM_SPADIC_CHA; iCh++){
          if (
              (cortest[iCh] > lowercorthreshold[susid]) 
              && (
                  (inPulse[iCh]->GetBinContent(inPulse[iCh]->GetMaximumBin()) < minAmplitude[susid])
                  || ((inPulse[iCh]->GetMaximumBin() < 5 || inPulse[iCh]->GetMaximumBin() > 20))        
                  )             
              ) break;
          signal_ch_counter++;
        }
        for (Int_t iCh = maxCh-1; iCh >= 0; iCh--){
          if (
              (cortest[iCh] > lowercorthreshold[susid]) 
              && (
                  (inPulse[iCh]->GetBinContent(inPulse[iCh]->GetMaximumBin()) < minAmplitude[susid])
                  || ((inPulse[iCh]->GetMaximumBin() < 5 || inPulse[iCh]->GetMaximumBin() > 20))
                  )             
              ) break;
          signal_ch_counter++;
        } 
      }      
      //clusterSize->Fill(NUM_SPADIC_CHA-noise_ch_counter); //debug
      clusterSize->Fill(signal_ch_counter); //debug
    }
  }
  principal->Clear();
  return signal_ch_counter;
}
void CalcDeltaSim(TH1D *delta, TH1D *sim/*, TH1D *measurement*/) {
  delta->SetYTitle(/*"#Delta [%]"*/"#frac{measured}{simulated}");
  delta->GetXaxis()->SetLabelSize(0.08);
  delta->GetYaxis()->SetLabelSize(0.08);
  delta->GetXaxis()->SetTitleSize(0.09);
  delta->GetYaxis()->SetTitleSize(0.09);
  delta->GetYaxis()->SetTitleOffset(0.5);
  delta->GetXaxis()->SetTitleOffset(1.0);
  //delta->Sumw2();
  delta->Divide(sim);
 
  /*
    for (Int_t bin = 1; bin < delta->GetNbinsX(); bin++) {
    Float_t dEdx = delta->GetBinCenter(bin);
    Float_t simValue = sim->GetBinContent(sim->FindBin(dEdx));
    Float_t measuredValue = measurement->GetBinContent(measurement->FindBin(dEdx));
    if (simValue > 0 && measuredValue > 0 && measurement->FindBin(dEdx) < measurement->GetNbinsX()+1)     
    delta->SetBinContent(bin, measuredValue/simValue);

    }
  */
}

  Double_t Calibration(TString file, Int_t suid) {

    TFile *output = new TFile(file,"READ");
    if (output->IsOpen()) {
      TF1 *uLandau = NULL;
      TF1 *sLandau = NULL;
      TH1D* uncalibrated = (TH1D*)output->Get("Uncalibrated_Pions"); 
      //TH1D* uncalibratedElectron = (TH1D*)output->Get("Uncalibrated_Electrons"); 
      /*
      // offset
      Int_t xBin = 1;
      Float_t threshold = uncalibrated->GetEntries() * 5.0e-4;
      while (uncalibrated->GetBinContent(xBin) <= threshold)
      xBin++;
      if (xBin > 1)
      calibrationParameter[1] = uncalibrated->GetBinCenter(xBin);//0.5 * (uncalibrated->GetBinCenter(xBin) + uncalibrated->GetBinCenter(xBin+1));
      else
      calibrationParameter[1] = 0.0;

      if (suid == 11)
      calibrationParameter[1] *= 0.75;
      else if (suid == 10)
      calibrationParameter[1] *= 0.90;
      else if (suid == 18)
      calibrationParameter[1] *= 0.65;
    
      // with offset compensation
      Double_t uncMax           = uncalibrated->GetBinContent(uncalibrated->GetMaximumBin());
      Double_t uncFullHeight    = uncalibrated->GetBinCenter(uncalibrated->FindLastBinAbove(0.50 * uncMax)) - calibrationParameter[1];
      Double_t uncPercentHeight = uncalibrated->GetBinCenter(uncalibrated->FindLastBinAbove(0.10 * uncMax)) - calibrationParameter[1];
      Double_t uncPermillHeight = uncalibrated->GetBinCenter(uncalibrated->FindLastBinAbove(0.05 * uncMax)) - calibrationParameter[1];
      */
      /*
      // without offset compensation
      Double_t uncMax           = uncalibrated->GetBinContent(uncalibrated->GetMaximumBin());
      Double_t uncFullHeight    = uncalibrated->GetBinCenter(uncalibrated->FindLastBinAbove(0.50 * uncMax));
      Double_t uncPercentHeight = uncalibrated->GetBinCenter(uncalibrated->FindLastBinAbove(0.10 * uncMax));
      Double_t uncPermillHeight = uncalibrated->GetBinCenter(uncalibrated->FindLastBinAbove(0.05 * uncMax));
      */
      const Double_t yFitRange = 0.50;
      //landau = new TF1("landau","landau",0,uncalibrated->GetBinCenter(uncalibrated->FindLastBinAbove(yFitRange * uncalibrated->GetBinContent(uncalibrated->GetMaximumBin()))));//3000);//2000
      uLandau = new TF1("uLandau","[0]*TMath::Landau(x,[1],[2],0)+[3]*TMath::Gaus(x,[1],[4],0)",200,uncalibrated->GetBinCenter(uncalibrated->FindLastBinAbove(yFitRange * uncalibrated->GetBinContent(uncalibrated->GetMaximumBin()))));//3000);//2000
      uLandau->SetParName(0,"lConstant");
      uLandau->SetParameter("lConstant",uncalibrated->GetBinContent(uncalibrated->GetMaximumBin()));

      uLandau->SetParName(1,"MPV");
      uLandau->SetParameter("MPV",700);
      uLandau->SetParLimits(1,500,810);//(1,514,644);

      uLandau->SetParName(2,"lSigma");
      uLandau->SetParameter("lSigma",190);
      uLandau->SetParLimits(2,100.0,220.0);

      uLandau->SetParName(3,"gConstant");
      uLandau->SetParameter("gConstant",0.025 * uncalibrated->GetBinContent(uncalibrated->GetMaximumBin()));
      uLandau->SetParLimits(3,0.0,uncalibrated->GetBinContent(uncalibrated->GetMaximumBin()));

      uLandau->SetParName(4,"gSigma");
      uLandau->SetParameter("gSigma",160);
      uLandau->SetParLimits(4,0.0,200);
      //landau->SetParLimits(2,114.7,175.0);//160.13); 
  
      uncalibrated->Fit(uLandau,"0QRB");//LL"); 
      printf("Unc:         Landau    Gaus\nConst.  %9.4f     %9.4f\nMPV     %9.4f     %9.4f\nSigma   %9.4f     %9.4f\n\n",uLandau->GetParameter("lConstant"),uLandau->GetParameter("gConstant"),uLandau->GetParameter("MPV"),uLandau->GetParameter("MPV"),uLandau->GetParameter("lSigma"),uLandau->GetParameter("gSigma"));    
      MPV_adc = uLandau->GetParameter("MPV");
      sigma_adc = uLandau->GetParameter("lSigma");
      ChiSqr_adc = uLandau->GetChisquare()/uLandau->GetNDF();
      TH2I *paras = NULL;
      TProfile *momentumMPV = NULL;
      TProfile *momentumSigma = NULL;
      TProfile *HVMPV = NULL;
      TProfile *HVSigma = NULL;
      TFile *fitFile = new TFile("CalibrationFitParameter.root","UPDATE");
      if (fitFile->IsOpen()) {
        TString title;
        title.Form("sigma_MPV_%i",suid);
        paras = (TH2I*)fitFile->Get(title);
        if (NULL == paras)
          paras = new TH2I(title,title,500,0,500,1500,0,1500);
        paras->Fill(sigma_adc, MPV_adc);

        title.Form("momentum_MPV_%i",suid);
        momentumMPV = (TProfile*)fitFile->Get(title);
        if (NULL == momentumMPV)
          momentumMPV = new TProfile(title,title,10,0.5,10.5);
        if (AnodeV == 1775)
          momentumMPV->Fill(momentum, MPV_adc);

        title.Form("momentum_Sigma_%i",suid);
        momentumSigma = (TProfile*)fitFile->Get(title);
        if (NULL == momentumSigma)
          momentumSigma = new TProfile(title,title,10,0.5,10.5);
        if (AnodeV == 1775)
          momentumSigma->Fill(momentum, sigma_adc);

        title.Form("HV_MPV_%i",suid);
        HVMPV = (TProfile*)fitFile->Get(title);
        if (NULL == HVMPV)
          HVMPV = new TProfile(title,title,3,1737.5,1812.5);
        if (momentum == 3)
          HVMPV->Fill(AnodeV, MPV_adc);

        title.Form("HV_Sigma_%i",suid);
        HVSigma = (TProfile*)fitFile->Get(title);
        if (NULL == HVSigma)
          HVSigma = new TProfile(title,title,3,1737.5,1812.5);
        if (momentum == 3)
          HVSigma->Fill(AnodeV, sigma_adc);

        fitFile->cd();
        paras->Write("", TObject::kOverwrite);
        momentumMPV->Write("", TObject::kOverwrite);
        momentumSigma->Write("", TObject::kOverwrite);
        HVMPV->Write("", TObject::kOverwrite);
        HVSigma->Write("", TObject::kOverwrite);
        fitFile->Close();
      }
      calibrationParameter[0] = uLandau->GetParameter("MPV");//uncalibrated->GetMean();
      //printf("\nuncalibrated:\n  MPV = %7.3f  offset = %7.3f (%7.3f) \n",calibrationParameter[0], calibrationParameter[1], threshold);
   
      //calibrationParameter[1] = landau->GetParameter("Sigma");
  
      TString simfile = "data2012/root/sim/Sim_xeco20_12mm_03GeV.root";
      simfile.Form("data2012/root/sim/Sim_xeco20_12mm_%02iGeV.root",momentum);
      //cout << simfile << endl;
      TFile *sim = new TFile(simfile.Data(),"READ");
      //TFile *sim = new TFile(/*"Pions_dEdx_3GeV_12mm_XeCo2.root"*/"data2012/root/sim/Sim_xeco20_12mm_03GeV.root","READ");
      TH1D* simulated= (TH1D*)sim->Get("hPi"); 
      /*
        Double_t simMax           = simulated->GetBinContent(simulated->GetMaximumBin());
        Double_t simFullHeight    = simulated->GetBinCenter(simulated->FindLastBinAbove(0.50 * simMax));
        Double_t simPercentHeight = simulated->GetBinCenter(simulated->FindLastBinAbove(0.10 * simMax));
        Double_t simPermillHeight = simulated->GetBinCenter(simulated->FindLastBinAbove(0.05 * simMax));
      */
      //landau = new TF1("landau","landau",0,simulated->GetBinCenter(simulated->FindLastBinAbove(yFitRange * simulated->GetBinContent(simulated->GetMaximumBin()))));//30);//10
      sLandau = new TF1("sLandau","[0]*TMath::Landau(x,[1],[2],0)+[3]*TMath::Gaus(x,[1],[4],0)",1,simulated->GetBinCenter(simulated->FindLastBinAbove(yFitRange * simulated->GetBinContent(simulated->GetMaximumBin()))));//30);//10
    
      sLandau->SetParName(0,"lConstant");
      sLandau->SetParameter("lConstant",simulated->GetBinContent(simulated->GetMaximumBin()));

      sLandau->SetParName(1,"MPV");
      sLandau->SetParLimits(1,0.0,5);
      sLandau->SetParameter("MPV",3.75);

      sLandau->SetParName(2,"lSigma");
      sLandau->SetParLimits(2,0.0,2);
      sLandau->SetParameter("lSigma",1.13);

      sLandau->SetParName(3,"gConstant");
      sLandau->SetParLimits(3,0.0, simulated->GetBinContent(simulated->GetMaximumBin()));
      sLandau->SetParameter("gConstant",0.025 * simulated->GetBinContent(simulated->GetMaximumBin()));

      sLandau->SetParName(4,"gSigma");
      sLandau->SetParLimits(4,0.0,10);
      sLandau->SetParameter("gSigma",1.5);

      simulated->Fit( sLandau,"0QRB");//LL");
      printf("Sim:         Landau    Gaus\nConst.  %9.4f     %9.4f\nMPV     %9.4f     %9.4f\nSigma   %9.4f     %9.4f\n\n",sLandau->GetParameter("lConstant"),sLandau->GetParameter("gConstant"),sLandau->GetParameter("MPV"),sLandau->GetParameter("MPV"),sLandau->GetParameter("lSigma"),sLandau->GetParameter("gSigma"));
      MPV_kev = sLandau->GetParameter("MPV");
      sigma_kev = sLandau->GetParameter("lSigma");
      ChiSqr_kev = sLandau->GetChisquare()/ sLandau->GetNDF();
      calibrationParameter[2] =  sLandau->GetParameter("MPV");//simulated->GetMean();
      printf("\nfit range: [0, %f]  [0, %f]\n",uncalibrated->GetBinCenter(uncalibrated->FindLastBinAbove(yFitRange * uncalibrated->GetBinContent(uncalibrated->GetMaximumBin()))),simulated->GetBinCenter(simulated->FindLastBinAbove(yFitRange * simulated->GetBinContent(simulated->GetMaximumBin()))));
      printf("\n     MPV      sigma     ChiSqr/NDF\nADC: %7.3f  %7.3f  %7.3f\nkev: %7.3f  %7.3f  %7.3f\n  1) m = %7.5f, b = %7.5f\n  2) m = %7.5f, b = %7.5f\n",MPV_adc,sigma_adc,ChiSqr_adc,MPV_kev,sigma_kev,ChiSqr_kev,
             (MPV_kev-sigma_kev)/(MPV_adc-sigma_adc),MPV_kev-(MPV_kev-sigma_kev)/(MPV_adc-sigma_adc)*MPV_adc,
             sigma_kev/sigma_adc,MPV_kev-sigma_kev/sigma_adc*MPV_adc
             );
      /*
        if (susid == 0){
        sigma_adc *= 0.85;
        MPV_adc *= 0.90;
        }
        if (susid == 2){
        sigma_adc *= 0.85;
        MPV_adc *= 0.85;
        }
      */
      if (susid == 0){
        if (momentum == 2){
          sigma_adc = 128;
          MPV_adc = 510;
        } else if (momentum == 3){
          if (AnodeV == 1750){
            sigma_adc = 150;
            MPV_adc = 515;//529.84
          } else if (AnodeV == 1775){
            sigma_adc = 175;
            MPV_adc = 629.78;//595.72//657.05
          } else if (AnodeV == 1800){
            sigma_adc = 190;
            MPV_adc = 696.98;
          }
        } else if (momentum == 4){
          sigma_adc = 170;
          MPV_adc = 750.00;
        } else if (momentum == 6){
          sigma_adc = 180.;
          MPV_adc = 650.;//630.12
        } else if (momentum == 8){
          sigma_adc = 190;//218.74
          MPV_adc = 650.00;//651.86
        }   
      }
      if (susid == 1){
        if (momentum == 2){
          sigma_adc = 140;//163.24
          MPV_adc = 600;
        } else if (momentum == 3){
          if (AnodeV == 1750){
            sigma_adc = 124.48;
            MPV_adc = 534.72;
          } else if (AnodeV == 1775){
            sigma_adc = 146.69;//136.99//144.65
            MPV_adc = 604.33;//580.00//608.40
          } else if (AnodeV == 1800){
            sigma_adc = 153;
            MPV_adc = 629.1;
          }
        } else if (momentum == 4){
          sigma_adc = 140;
          MPV_adc = 668.04;
        } else if (momentum == 6){
          sigma_adc = 160;
          MPV_adc = 654.10;
        } else if (momentum == 8){
          sigma_adc = 173.48;
          MPV_adc = 670.05;
        }   
      }
      if (susid == 2){
        if (momentum == 2){

        } else if (momentum == 3){
          if (AnodeV == 1750){
            sigma_adc = 150;
            MPV_adc = 520;
          } else if (AnodeV == 1775){
            sigma_adc = 170;
            MPV_adc = 560;//850
          } else if (AnodeV == 1800){
            sigma_adc = 180;
            MPV_adc = 550;
          }
        } else if (momentum == 4){

        } else if (momentum == 6){

        } else if (momentum == 8){

        }
      }
      
  
   
      if (suid == 11){
        //MPV_adc = 595.506;
        //sigma_adc = 136.429;
        //calibrationParameter[3] = 1.22 * (simFullHeight / uncFullHeight + simPercentHeight / uncPercentHeight + simPermillHeight / uncPermillHeight) / 3.0;
        //calibrationParameter[2] *= 1.25;//without offset
        //calibrationParameter[2] *= 1.25;//with offset
      }
      else if (suid == 10){
        //calibrationParameter[3] = 1.00 * (simFullHeight / uncFullHeight + simPercentHeight / uncPercentHeight + simPermillHeight / uncPermillHeight) / 3.0;
        //calibrationParameter[2] *= 1.22;//without offset
        //calibrationParameter[2] *= 1.10;//with offset
      }
      else if (suid == 18){
        //calibrationParameter[3] = 2.00 * (simFullHeight / uncFullHeight + simPercentHeight / uncPercentHeight + simPermillHeight / uncPermillHeight) / 3.0;
        //calibrationParameter[2] *= 2.25;//without offset
        //calibrationParameter[2] *= 2.35;//with offset
      }
      /*
        printf("simulated:\n  MPV = %7.3f  offset = %7.3f \n",landau->GetParameter("MPV"), 0.0);
        printf("\n  (ADCvalue / %7.3f ) * %7.3f     -> ADCvalue * %7.6f\n",calibrationParameter[0],calibrationParameter[2],calibrationParameter[2]/calibrationParameter[0]);
        printf("\n  ADCvalue * %7.6f\n\n",calibrationParameter[3]);

        printf("  maxA:%.2f  %.2f  %.2f  %.2f\n  maxA:%.2f  %.2f  %.2f  %.2f\n",uncMax,uncFullHeight,uncPercentHeight,uncPermillHeight,simMax,simFullHeight,simPercentHeight,simPermillHeight);
      */
      printf("\nfit range: [0, %f]  [0, %f]\n",uncalibrated->GetBinCenter(uncalibrated->FindLastBinAbove(yFitRange * uncalibrated->GetBinContent(uncalibrated->GetMaximumBin()))),simulated->GetBinCenter(simulated->FindLastBinAbove(yFitRange * simulated->GetBinContent(simulated->GetMaximumBin()))));
      printf("\n     MPV      sigma     ChiSqr/NDF\nADC: %7.3f  %7.3f  %7.3f\nkev: %7.3f  %7.3f  %7.3f\n  1) m = %7.5f, b = %7.5f\n  2) m = %7.5f, b = %7.5f\n",MPV_adc,sigma_adc,ChiSqr_adc,MPV_kev,sigma_kev,ChiSqr_kev,
             (MPV_kev-sigma_kev)/(MPV_adc-sigma_adc),MPV_kev-(MPV_kev-sigma_kev)/(MPV_adc-sigma_adc)*MPV_adc,
             sigma_kev/sigma_adc,MPV_kev-sigma_kev/sigma_adc*MPV_adc
             );
      output->Close();
      return  MPV_adc;//calibrationParameter[0];//landau->GetParameter("MPV");//*0.85;
    }
    else
      return 1;
  }

void FillUnCalibratedPionSpectrum(TString filename, TString file,  Int_t usedSuId, Bool_t amplitude, Bool_t debug) {
  printf("---FillUnCalibratedPionSpectrum---\n");
  TFile *output = new TFile(file,"UPDATE");
  if (output->IsOpen()) {
    TH1D* uncalibrated = new TH1D("Uncalibrated_Pions","Uncalibrated Pion Spectrum",200,0,100*fHistoScale);
    TH1D* uncalibratedElectron = new TH1D("Uncalibrated_Electrons","Uncalibrated Electron Spectrum",200,0,100*fHistoScale);
    //TH1I* hitTime = new TH1I();
    //TH1I *amplitude = new TH1I();

    TH1I *clusterSize = new TH1I("UC_clusterSize","clusterSize",NUM_SPADIC_CHA+1,-0.5,NUM_SPADIC_CHA+0.5);
    TH1I *maxAmplitudeValue = new TH1I("UC_maxAmplitudeValue","maxAmplitudeValue",256*2,0,256);
    TH2I *maxAmplitudeHitTime = new TH2I("UC_maxAmplitudeHitTime","maxAmplitudeHitTime",256,0,256,45,0,45);
    TH1I *covaMatixValue = new TH1I("UC_covaMatixValue","covaMatixValue",1000,-0.13,0.13);
    TH2I *covaMatixValueMaxAmplitude = new TH2I("UC_covaMatixValueMaxAmplitude","covaMatixValueMaxAmplitude",1000,-0.13,0.13,256,0,256);
    TH2I *covaMatixValueHitTime = new TH2I("UC_covaMatixValueHitTime","covaMatixValueHitTime",1000,-0.13,0.13,45,0,45);
    TH1I *signalChDistance = new TH1I("UC_signalChDistance","signalChDistance",2*NUM_SPADIC_CHA+1,-NUM_SPADIC_CHA-0.5,NUM_SPADIC_CHA+0.5);
    TH2I *averageSignal_2D = new TH2I("UC_averageSignal_2D","averageSignal_2D",SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE,356,-100,256);
    TH2I *averageNoise_2D = new TH2I("UC_averageNoise_2D","averageNoise_2D",SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE,356,-100,256);
    TH1I *noiseDistribution = new TH1I("UC_noiseDistribution","noiseDistribution",2*SPADIC_ADC_MAX,-SPADIC_ADC_MAX,SPADIC_ADC_MAX);
    TH2I *baselineDistribution = new TH2I("UC_baselineDistribution","baselineDistribution",NUM_SPADIC_CHA,-0.5,NUM_SPADIC_CHA-0.5,SPADIC_ADC_MAX,0,SPADIC_ADC_MAX);

    TH1D *rawPulse[NUM_SPADIC_CHA];
    TH1D *baselinePulse[NUM_SPADIC_CHA];
    TH1D *baselineNoisePulse[NUM_SPADIC_CHA];
    TString name;
    for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {
      name.Form("Ch%02i",ch);
      rawPulse[ch] = new TH1D("rawPulse"+name,"rawPulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      baselinePulse[ch] = new TH1D("baselinePulse"+name,"baselinePulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      baselineNoisePulse[ch] = new TH1D("baselineNoisePulse"+name,"baselineNoisePulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
    }
    name = filename;
    name.ReplaceAll(".root","");
    GetPidCuts(name.ReplaceAll("data2012/TTrees/merged/",""));
    TFile inputFile(filename,"READ");
    if (inputFile.IsOpen()) {
      TTree* theTree=NULL;
      TKey* kee=NULL;
      TIter iter(inputFile.GetListOfKeys());
      while ( ( kee=dynamic_cast<TKey*>(iter()) ) !=0 ) {
        theTree = dynamic_cast<TTree*>(kee->ReadObj());
        if (theTree)
          break; // we take first Tree in file, disregarding its name...
      }
      if(theTree==0) {
        cout <<"Error: no Tree in file "<< filename.Data() << endl;
        return;
      }
      TCernOct12UnpackEvent* evnt = new TCernOct12UnpackEvent;
      TGo4EventElement* theBase=evnt;
      evnt->synchronizeWithTree(theTree, &theBase);
      Int_t entries = theTree->GetEntries();
      Int_t nPion(0),nElectron(0);
      printf("%i Events found in TTree\n",entries);
      TMbsCrateEvent* fCrateInputEvent=NULL;
      TSpadicEvent* SpadicInputEvent=NULL;
      TSpadicData* theSpadic=NULL;
      Float_t Ch1(0), Ch2(0), Pb(0);
      //Bool_t isPulser = false;
      Int_t clusterCh = 0;
      if (debug) entries /= 100;
      for(Int_t i=0;i<entries;++i) {
        Statusbar(i,entries);
        theTree->GetEntry(i);
        fCrateInputEvent=dynamic_cast<TMbsCrateEvent*>(evnt->GetSubEvent("MBSCRATE"));    
        SpadicInputEvent=dynamic_cast<TSpadicEvent*>(evnt->GetSubEvent("SPADIC"));
        isPion = false;
        isElectron = false;
        isOverFlowEvent = false;
        isUnderFlowEvent = false;
        isPulser = fCrateInputEvent->fIsPulser;
        if (isPulser) continue;
        if ((Ch1 < 350 && Pb < 350) || (Ch2 < 350 && Pb < 350)) continue;
    
        /*
        //2011
        Ch1 = fCrateInputEvent->fData1182[1][0];
        Ch2 = fCrateInputEvent->fData1182[0][1];
        Pb  = fCrateInputEvent->fData1182[1][5];      
        */
        //2012
        Ch1 = fCrateInputEvent->fData1182[1][0];
        Ch2 = fCrateInputEvent->fData1182[0][0];
        Pb  = fCrateInputEvent->fData1182[1][1];
        getPID(Ch1, Ch2, Pb);
        if (isPion)
          nPion++;
        if (isElectron)
          nElectron++;
      
        //for (Int_t suid = 0; suid < /*usedSusibos*/1; suid++){
        theSpadic=dynamic_cast<TSpadicData*>(SpadicInputEvent->getEventElement(usedSuId));

        if(NULL == theSpadic || !theSpadic) {cout << "no data found for Susibo " << usedSuId << endl; return;}
        if (isPion || isElectron){
          for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {       
            rawPulse[ch]->Reset();
            baselinePulse[ch]->Reset();
            baselineNoisePulse[ch]->Reset();
     
            for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {       
              rawPulse[ch]->Fill(bin, (Double_t)theSpadic->fSpadicPulse[ch][bin]);
            }     
          }
   
          BaselineCompensation(rawPulse, baselinePulse, baselineDistribution);
          clusterCh = CancelNoise_Cova(baselinePulse, baselineNoisePulse, covaMatixValue,covaMatixValueMaxAmplitude,covaMatixValueHitTime,maxAmplitudeHitTime,noiseDistribution,clusterSize,signalChDistance,averageNoise_2D,averageSignal_2D,NULL,false);


          //inPulse[ch]->GetBinContent(inPulse[ch]->GetMaximumBin())

          if(amplitude){
            if (isElectron || isPion)
              for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++)
                maxAmplitudeValue->Fill(baselineNoisePulse[ch]->GetBinContent(baselineNoisePulse[ch]->GetMaximumBin()));
            if(minimumAmplitudeThreshold(baselineNoisePulse) && !BorderPadMax(baselineNoisePulse) && HitTimeWindow(baselineNoisePulse) && minimumIntegralThreshold(baselineNoisePulse)){
              if (isPion){
                uncalibrated->Fill(Get3PadPosAmplitude(baselineNoisePulse));
              }
              else if (isElectron){
                uncalibratedElectron->Fill(Get3PadPosAmplitude(baselineNoisePulse));
              }
            }
          }     
          else{

            for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++)
              maxAmplitudeValue->Fill(baselineNoisePulse[ch]->GetBinContent(baselineNoisePulse[ch]->GetMaximumBin()));
            if(minimumAmplitudeThreshold(baselineNoisePulse) && !BorderPadMax(baselineNoisePulse) && HitTimeWindow(baselineNoisePulse) && minimumIntegralThreshold(baselineNoisePulse)){
              if (isPion){
                uncalibrated->Fill(Get3PadPosIntegral(baselineNoisePulse));     
              }
              else if (isElectron){
                uncalibratedElectron->Fill(Get3PadPosIntegral(baselineNoisePulse));
              }
            }
          }
        }
      }
      printf("%.2f%% pions, %.2f%% electrons and %.2f%% no PID or noise\n",Float_t(nPion)/entries*100,Float_t(nElectron)/entries*100, Float_t(entries-(nElectron+nPion))/entries*100);
      output->cd();
      uncalibrated->Write("", TObject::kOverwrite);
      uncalibratedElectron->Write("", TObject::kOverwrite);

      maxAmplitudeValue->Write("", TObject::kOverwrite);
      covaMatixValue->Write("", TObject::kOverwrite);
      covaMatixValueMaxAmplitude->Write("", TObject::kOverwrite);
      covaMatixValueHitTime->Write("", TObject::kOverwrite);
      clusterSize->Write("", TObject::kOverwrite);
      maxAmplitudeHitTime->Write("", TObject::kOverwrite);
  
      noiseDistribution->SetXTitle("noise [ADC]");
      noiseDistribution->SetYTitle("#");
      noiseDistribution->Write("", TObject::kOverwrite);
      baselineDistribution->SetXTitle("channel");
      baselineDistribution->SetYTitle("ADC");
      baselineDistribution->Write("", TObject::kOverwrite);

      output->Close();
      inputFile.Close();
     
      /*
        delete theTree;
        delete kee;
        delete evnt;
        delete theBase;
        delete fCrateInputEvent;
        delete SpadicInputEvent;
        delete theSpadic;
      
        delete uncalibrated; 
      
        for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {
        delete rawPulse[ch];
        delete baselinePulse[ch];
        delete baselineNoisePulse[ch];
        }
      */
    }
    else {
      printf("ERROR::FillUnCalibratedPionSpectrum: %s not found",filename.Data());
    }
  }
  else {
    printf("ERROR::FillUnCalibratedPionSpectrum: %s not found",file.Data());

  }
}
  //void spadic_noise_Analysis(){

  //}

void spadic_noise_Analysis(TString filename=/*"data2011/TTrees/Be_run2110020.root"*/"data2012/TTrees/Be_run2310004.root",TString radiator = "D", const Int_t suid = 1, Bool_t debug=false, Bool_t simple = false, Bool_t plotFromFile = true)
{
  TStopwatch timer;
  timer.Start();
  Bool_t fast = true;
  gROOT->Reset(); 
  gROOT->SetStyle("Plain");
  gStyle->SetPalette(1,0);
  gStyle->SetOptStat(kFALSE);
  gStyle->SetOptTitle(kFALSE);
  gStyle->SetPadTickX(1);                         //<-- tic marks on all axes
  gStyle->SetPadTickY(1);  
  Bool_t amplitude = false;
 
  if (suid > 4){
    printf("ERROR:: did not find susibo %i\n",suid);
    return;
  }
  susid = suid;
  TString tempcast = filename;
  TString replaceString;
  momentum = TString((tempcast.ReplaceAll("data2012/TTrees/merged/",""))(0,1)).Atoi();
  replaceString.Form("%iGeV_",momentum);
  AnodeV = TString((tempcast.ReplaceAll(replaceString,""))(0,4)).Atoi();
  replaceString.Form("%iV_",AnodeV);
  DriftV = TString((tempcast.ReplaceAll(replaceString,""))(0,3)).Atoi();
  if (AnodeV == 0)
    AnodeV = 1800;
  if (DriftV == 0)
    DriftV = 500;
  printf("\n  momentum: %iGeV  anode: %iV  drift: %iV\n", momentum, AnodeV, DriftV);
  TString simfile = "data2012/root/sim/Sim_xeco20_12mm_03GeV.root";
  simfile.Form("data2012/root/sim/Sim_xeco20_12mm_%02iGeV.root",momentum);
  cout << "  sim. file path: " << simfile << endl;
  const Int_t usedSusibos = 5;
  //Int_t usedSuId[usedSusibos] = {         2,         12,        6,        18,         4,         5,        16,        17  };
  //TString usedSuName[usedSusibos] = {"MS336/10","MS336/11","MS444/10","MS444/11","FFM01","FFM02","FFM03","FFM04"};
  //Double_t padWidth[usedSusibos] = {5,5,8,8,0,0,0,0};
  //Double_t h[usedSusibos] = {3,3,4,4,0,0,0,0};
  Int_t usedSuId[usedSusibos] = {            11,      10,        18,        5,       19};
  TString usedSuName[usedSusibos] = {"MS_R_I","MS_N_II","MS_J_III","FFM_5+5","FFM_4+4"};
  Double_t padWidth[usedSusibos] = {     7.125,   7.125,     7.125,    7.125,    7.125};
  Double_t h[usedSusibos] = {              3.5,     3.5,       3.5,        5,        4};
  Double_t K3 = calcK3(h[0]);//0.525;
  Double_t par = 1.0;
  TString outpic;
  TString outfilename;
  outpic = filename;
  outpic.ReplaceAll("data2012/TTrees/merged/","");
  outpic.ReplaceAll(".root","");
  outfilename.Form("data2012/root/Spectra_%s_Rad_%s_Sus%02i.root",outpic.Data(),radiator.Data(),usedSuId[suid]);
  outpic = outfilename;
  outpic.ReplaceAll("root/","pics/");
  outpic.ReplaceAll(".root","");
  Double_t MPV_pion = Calibration(outfilename, usedSuId[suid]);
  if (MPV_pion == 1) {
    FillUnCalibratedPionSpectrum(filename,outfilename, usedSuId[suid], amplitude, debug);
    MPV_pion = Calibration(outfilename, usedSuId[suid]);
  }
  //printf("MPV_pion: %.3e\n",MPV_pion);
  //Double_t scaling2keV = 3.60898e+00/MPV_pion;//3.65338e+00/MPV_pion;//3.67179/MPV_pion;//7.87153e+02;//3.67179/6.67190e+02;
  Int_t nbins(200),lbin(0/*-1000*/),hbin(100);
  
  if(amplitude) {
    lbin /= 20;
    hbin /= 20;
  }
  else {
    nbins = 200;
    lbin = 0;
    if (fHistoScale == SPADIC_TRACE_SIZE * SPADIC_ADC_MAX / 100.)
      hbin = 100*fHistoScale;
  }
  TFile *output = new TFile(outfilename,"UPDATE");
  TH2F* Ch1_Pb = NULL;
  TH2F* allCh1_Pb = NULL;
  TH2F* noCh1_Pb = NULL;
  TH2F* Ch2_Pb = NULL;
  TH2F* allCh2_Pb = NULL;
  TH2F* noCh2_Pb = NULL;
  TH2F* Ch1_Ch2 = NULL;
  TH2F* allCh1_Ch2 = NULL;
  TH2F* noCh1_Ch2 = NULL;
  TH1I* hCh_1 = NULL;
  TH1I *hCh_2 = NULL;
  TH1I *hPb = NULL;
  TH2D* clusterwCharge = NULL;
  TH2D* prf = NULL;
  TProfile* prfProfile = NULL;

  TH1I *clusterSize = NULL;
  TH1I *maxAmplitudeValue = NULL;
  TH2I *maxAmplitudeHitTime = NULL;
  TH1I *covaMatixValue = NULL;
  TH2I *covaMatixValueMaxAmplitude = NULL;
  TH2I *covaMatixValueHitTime = NULL;
  TH1I *signalChDistance = NULL;
  TH2I *averageSignal_2D = NULL;
  TH2I *averageNoise_2D = NULL;
  TH1I *noiseDistribution = NULL;
  TH2I *baselineDistribution = NULL;


  TH1D *preCompPulse[NUM_SPADIC_CHA] = {NULL};
  TH1D *preCompBaselinePulse[NUM_SPADIC_CHA] = {NULL};
  TH1D *rawPulse[NUM_SPADIC_CHA] = {NULL};
  TH1D *baselinePulse[NUM_SPADIC_CHA] = {NULL};
  TH1D *noisePulse[NUM_SPADIC_CHA] = {NULL};
  TH1D *baselineNoisePulse[NUM_SPADIC_CHA] = {NULL};
  TH1D *noiseBaselinePulse[NUM_SPADIC_CHA] = {NULL};
  TH1D *lastPulse[NUM_SPADIC_CHA] = {NULL};
  TH2I *offSpillPulse[NUM_SPADIC_CHA] = {NULL};

  TProfile *averagePulse_e = NULL;
  TProfile *averagePulse_p = NULL;


  TH2I *averagePulse_2D_e = NULL;
  TH2I *averagePulse_2D_p = NULL;

  TH1D *preCompSpectrum = NULL;
  TH1D *preCompBaselineSpectrum = NULL;
  TH1D *preCompSpectrumElectron = NULL;
  TH1D *preCompBaselineSpectrumElectron = NULL;
  TH1D *rawSpectrum = NULL;
  TH1D *baselineSpectrum = NULL;
  TH1D *noiseSpectrum = NULL;
  TH1D *baselineNoiseSpectrum = NULL;
  TH1D *baselineNoiseSpectrumWoOWoU = NULL;
  TH1D *baselineNoiseSpectrumElectron = NULL;
  TH1D *baselineNoiseSpectrumWoOWoUElectron = NULL;
  TH1D *noiseBaselineSpectrum = NULL;
  TH1D *probabilityOverflow = NULL;
  TH1D *probabilityUnderflow = NULL;
  TH1D *probabilityOverflowE = NULL;
  TH1D *probabilityUnderflowE = NULL;
  TH2I* keV2Amplitude = new TH2I("keV2Amplitude","keV2Amplitude",100,0,100,256,0,256);
  keV2Amplitude->SetContour(99);
  TProfile* keV2AmplitudeProfile = new TProfile("keV2AmplitudeProfile","keV2AmplitudeProfile",100,0,100);


  TH2I* IntegralvsAmplitude = new TH2I("IntegralvsAmplitude","IntegralvsAmplitude",NUM_SPADIC_CHA*SPADIC_TRACE_SIZE*256,0,NUM_SPADIC_CHA*SPADIC_TRACE_SIZE*256,NUM_SPADIC_CHA*SPADIC_TRACE_SIZE,0,NUM_SPADIC_CHA*SPADIC_TRACE_SIZE);
  IntegralvsAmplitude->SetXTitle("Integrated ADC values 3 pad cluster");
  IntegralvsAmplitude->SetYTitle("Maximum ADC amplitude on pad max");
  IntegralvsAmplitude->SetContour(99);

  TProfile *crossTalk = NULL;//new TProfile("CrossTalk","CrossTalk",6,-0.5,5.5);

  TH1D *baselineNoiseSpectrumCrossTalkE[6] = {NULL};
  TH1D *baselineNoiseSpectrumCrossTalkP[6] = {NULL};

  TString name;
  Int_t color[NUM_SPADIC_CHA] = {1,2,3,4,800,6,7,8};

  if (plotFromFile) {
    output->cd();
    clusterSize = (TH1I*)output->Get("clusterSize");
    maxAmplitudeValue = (TH1I*)output->Get("maxAmplitudeValue");
    maxAmplitudeHitTime = (TH2I*)output->Get("maxAmplitudeHitTime");
    covaMatixValue = (TH1I*)output->Get("covaMatixValue");
    covaMatixValueMaxAmplitude = (TH2I*)output->Get("covaMatixValueMaxAmplitude");
    covaMatixValueHitTime = (TH2I*)output->Get("covaMatixValueHitTime");
    signalChDistance = (TH1I*)output->Get("signalChDistance");
    averageSignal_2D = (TH2I*)output->Get("averageSignal_2D");
    averageNoise_2D = (TH2I*)output->Get("averageNoise_2D");
    noiseDistribution = (TH1I*)output->Get("noiseDistribution");
    baselineDistribution = (TH2I*)output->Get("baselineDistribution");
    Ch1_Pb = (TH2F*)output->Get("Ch1_Pb");
    allCh1_Pb = (TH2F*)output->Get("allCh1_Pb");
    if (!Ch1_Pb)
      cout << "Ch1_Pb" << endl;
    noCh1_Pb = (TH2F*)output->Get("noCh1_Pb");
    if (!noCh1_Pb)
      cout << "noCh1_Pb" << endl;
    Ch2_Pb = (TH2F*)output->Get("Ch2_Pb");
    allCh2_Pb = (TH2F*)output->Get("allCh2_Pb");
    if (!Ch2_Pb)
      cout << "Ch2_Pb" << endl;
    noCh2_Pb = (TH2F*)output->Get("noCh2_Pb");
    if (!noCh2_Pb)
      cout << "noCh2_Pb" << endl;
    Ch1_Ch2 = (TH2F*)output->Get("Ch1_Ch2");
    allCh1_Ch2 = (TH2F*)output->Get("allCh1_Ch2");
    if (!Ch1_Ch2)
      cout << "Ch1_Ch2" << endl;
    noCh1_Ch2 = (TH2F*)output->Get("noCh1_Ch2");
    if (!noCh1_Ch2)
      cout << "noCh1_Ch2" << endl;

    clusterwCharge = (TH2D*)output->Get("clusterwCharge");

    prf = (TH2D*)output->Get("PRF_2D");
    if (!prf)
      cout << "PRF_2D" << endl;
    prfProfile = (TProfile*)output->Get("PRF");
    if (!prfProfile)
      cout << "PRF" << endl;
  
    averagePulse_e = (TProfile*)output->Get("averagePulse_e");
    if (!averagePulse_e)
      cout << "averagePulse_e" << endl;
    averagePulse_p = (TProfile*)output->Get("averagePulse_p");
    if (!averagePulse_p)
      cout << "averagePulse_p" << endl;

    averagePulse_2D_e = (TH2I*)output->Get("averagePulse_2D_e");
    if (!averagePulse_2D_e)
      cout << "averagePulse_2D_e" << endl;
    averagePulse_2D_p = (TH2I*)output->Get("averagePulse_2D_p");
    if (!averagePulse_2D_p)
      cout << "averagePulse_2D_p" << endl;


    preCompSpectrum = (TH1D*)output->Get("preCompSpectrum");
    if (!preCompSpectrum)
      cout << "preCompSpectrum" << endl;
    preCompBaselineSpectrum = (TH1D*)output->Get("preCompBaselineSpectrum");
    if (!preCompBaselineSpectrum)
      cout << "preCompBaselineSpectrum" << endl;
    preCompSpectrumElectron = (TH1D*)output->Get("preCompSpectrumElectron");
    if (!preCompSpectrumElectron)
      cout << "preCompSpectrumElectron" << endl;
    preCompBaselineSpectrumElectron = (TH1D*)output->Get("preCompBaselineSpectrumElectron");
    if (!preCompBaselineSpectrumElectron)
      cout << "preCompBaselineSpectrumElectron" << endl;
    rawSpectrum = (TH1D*)output->Get("rawSpectrum");
    if (!rawSpectrum)
      cout << "rawSpectrum" << endl;
    baselineSpectrum = (TH1D*)output->Get("baselineSpectrum");
    if (!baselineSpectrum)
      cout << "baselineSpectrum" << endl;
    noiseSpectrum = (TH1D*)output->Get("noiseSpectrum");
    if (!noiseSpectrum)
      cout << "noiseSpectrum" << endl;
    baselineNoiseSpectrum = (TH1D*)output->Get("baselineNoiseSpectrum");
    if (!baselineNoiseSpectrum)
      cout << "baselineNoiseSpectrum" << endl;
    baselineNoiseSpectrumWoOWoU = (TH1D*)output->Get("baselineNoiseSpectrumWoOWoU");
    if (!baselineNoiseSpectrumWoOWoU)
      cout << "baselineNoiseSpectrumWoOWoU" << endl;
    baselineNoiseSpectrumElectron = (TH1D*)output->Get("baselineNoiseSpectrumElectron");
    if (!baselineNoiseSpectrumElectron)
      cout << "baselineNoiseSpectrumElectron" << endl;
    baselineNoiseSpectrumWoOWoUElectron = (TH1D*)output->Get("baselineNoiseSpectrumWoOWoUElectron");
    if (!baselineNoiseSpectrumWoOWoUElectron)
      cout << "baselineNoiseSpectrumWoOWoUElectron" << endl;
    noiseBaselineSpectrum = (TH1D*)output->Get("noiseBaselineSpectrum");
    if (!noiseBaselineSpectrum)
      cout << "noiseBaselineSpectrum" << endl;
    probabilityOverflow = (TH1D*)output->Get("probabilityOverflow");
    if (!probabilityOverflow)
      cout << "probabilityOverflow" << endl;
    probabilityUnderflow = (TH1D*)output->Get("probabilityUnderflow");
    if (!probabilityUnderflow)
      cout << "probabilityUnderflow" << endl;
    probabilityOverflowE = (TH1D*)output->Get("probabilityOverflowE");
    if (!probabilityOverflowE)
      cout << "probabilityOverflowE" << endl;
    probabilityUnderflowE = (TH1D*)output->Get("probabilityUnderflowE");
    if (!probabilityUnderflowE)
      cout << "probabilityUnderflowE" << endl;

    crossTalk = (TProfile *)output->Get("CrossTalk");
    if (!crossTalk)
      cout << "crossTalk" << endl;

    for(Int_t delta = 0; delta < 6; delta++) {
      name.Form("_delta%02i",delta);
      baselineNoiseSpectrumCrossTalkE[delta] = (TH1D*)output->Get("baselineNoiseSpectrumCrossTalkE"+name);
      if (!baselineNoiseSpectrumCrossTalkE[delta])
        cout << TString("baselineNoiseSpectrumCrossTalkE"+name) << endl;
      baselineNoiseSpectrumCrossTalkP[delta] = (TH1D*)output->Get("baselineNoiseSpectrumCrossTalkP"+name);
      if (!baselineNoiseSpectrumCrossTalkP[delta])
        cout << TString("baselineNoiseSpectrumCrossTalkP"+name) << endl;
    }
  }
  else {


    clusterSize = new TH1I("clusterSize","clusterSize",NUM_SPADIC_CHA+1,-0.5,NUM_SPADIC_CHA+0.5);
    maxAmplitudeValue = new TH1I("maxAmplitudeValue","maxAmplitudeValue",256*2,0,256);
    maxAmplitudeHitTime = new TH2I("maxAmplitudeHitTime","maxAmplitudeHitTime",256,0,256,45,0,45);
    covaMatixValue = new TH1I("covaMatixValue","covaMatixValue",1000,-0.13,0.13);
    covaMatixValueMaxAmplitude = new TH2I("covaMatixValueMaxAmplitude","covaMatixValueMaxAmplitude",1000,-0.13,0.13,256,0,256);
    covaMatixValueHitTime = new TH2I("covaMatixValueHitTime","covaMatixValueHitTime",1000,-0.13,0.13,45,0,45);
    signalChDistance = new TH1I("signalChDistance","signalChDistance",2*NUM_SPADIC_CHA+1,-NUM_SPADIC_CHA-0.5,NUM_SPADIC_CHA+0.5);
    averageSignal_2D = new TH2I("averageSignal_2D","averageSignal_2D",SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE,356,-100,256);
    averageNoise_2D = new TH2I("averageNoise_2D","averageNoise_2D",SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE,356,-100,256);
    noiseDistribution = new TH1I("noiseDistribution","noiseDistribution",2*SPADIC_ADC_MAX,-SPADIC_ADC_MAX,SPADIC_ADC_MAX);
    baselineDistribution = new TH2I("baselineDistribution","baselineDistribution",NUM_SPADIC_CHA,-0.5,NUM_SPADIC_CHA-0.5,SPADIC_ADC_MAX,0,SPADIC_ADC_MAX);
    clusterSize->SetXTitle("cluster width [channels]");
    clusterSize->SetYTitle("normalized counts");
    signalChDistance->SetXTitle("distance to next signal channel [channels]");
    signalChDistance->SetYTitle("normalized counts");
    maxAmplitudeHitTime->SetXTitle("max. amplitude [ADC]");
    maxAmplitudeHitTime->SetYTitle("max. time bin");
    maxAmplitudeHitTime->SetContour(99);
    baselineDistribution->SetContour(99);
    covaMatixValueMaxAmplitude->SetContour(99);
    covaMatixValueHitTime->SetContour(99);

    maxAmplitudeValue->SetXTitle("max. amplitude [ADC]");
    covaMatixValue->SetXTitle("cova. value");

    covaMatixValueMaxAmplitude->SetYTitle("max. amplitude [ADC]");
    covaMatixValueMaxAmplitude->SetXTitle("cova. value");

    covaMatixValueHitTime->SetXTitle("cova. value");
    covaMatixValueHitTime->SetYTitle("max. time bin");

    hCh_1 = new TH1I("Ch_1","Ch_1",1000,0,4000);
    hCh_1->SetXTitle("Cherenkov 1 [a.u.]");

    hCh_2 = new TH1I("Ch_2","Ch_2",1000,0,4000);
    hCh_2->SetXTitle("Cherenkov 2 [a.u.]");

    hPb = new TH1I("Pb","Pb",1000,0,4000);
    hPb->SetYTitle("Pb-glass calorimeter [a.u.]");

    Ch1_Pb = new TH2F("Ch1_Pb","good PID Ch1_Pb",400,0,4000,400,0,4000);
    Ch1_Pb->SetXTitle("Cherenkov 1 [a.u.]");
    Ch1_Pb->SetYTitle("Pb-glass calorimeter [a.u.]");
    Ch1_Pb->SetContour(99);

    allCh1_Pb = new TH2F("allCh1_Pb","all PID Ch1_Pb",400,0,4000,400,0,4000);
    allCh1_Pb->SetXTitle("Cherenkov 1 [a.u.]");
    allCh1_Pb->SetYTitle("Pb-glass calorimeter [a.u.]");
    allCh1_Pb->SetContour(99);

    noCh1_Pb = new TH2F("noCh1_Pb","no good PID Ch1_Pb",400,0,4000,400,0,4000);
    noCh1_Pb->SetXTitle("Cherenkov 1 [a.u.]");
    noCh1_Pb->SetYTitle("Pb-glass calorimeter [a.u.]");
    noCh1_Pb->SetContour(99);

    Ch2_Pb = new TH2F("Ch2_Pb","good PID Ch2_Pb",400,0,4000,400,0,4000);
    Ch2_Pb->SetXTitle("Cherenkov 2 [a.u.]");
    Ch2_Pb->SetYTitle("Pb-glass calorimeter [a.u.]");
    Ch2_Pb->SetContour(99);

    allCh2_Pb = new TH2F("allCh2_Pb","all PID Ch2_Pb",400,0,4000,400,0,4000);
    allCh2_Pb->SetXTitle("Cherenkov 2 [a.u.]");
    allCh2_Pb->SetYTitle("Pb-glass calorimeter [a.u.]");
    allCh2_Pb->SetContour(99);

    noCh2_Pb = new TH2F("noCh2_Pb","no good PID Ch2_Pb",400,0,4000,400,0,4000);
    noCh2_Pb->SetXTitle("Cherenkov 2 [a.u.]");
    noCh2_Pb->SetYTitle("Pb-glass calorimeter [a.u.]");
    noCh2_Pb->SetContour(99);

    Ch1_Ch2 = new TH2F("Ch1_Ch2","good PID Ch1_Ch2",400,0,4000,400,0,4000);
    Ch1_Ch2->SetXTitle("Cherenkov 1 [a.u.]");
    Ch1_Ch2->SetYTitle("Cherenkov 2 [a.u.]");
    Ch1_Ch2->SetContour(99);

    allCh1_Ch2 = new TH2F("allCh1_Ch2","all PID Ch1_Ch2",400,0,4000,400,0,4000);
    allCh1_Ch2->SetXTitle("Cherenkov 1 [a.u.]");
    allCh1_Ch2->SetYTitle("Cherenkov 2 [a.u.]");
    allCh1_Ch2->SetContour(99);

    noCh1_Ch2 = new TH2F("noCh1_Ch2","no good PID Ch1_Ch2",400,0,4000,400,0,4000);
    noCh1_Ch2->SetXTitle("Cherenkov 1 [a.u.]");
    noCh1_Ch2->SetYTitle("Cherenkov 2 [a.u.]");
    noCh1_Ch2->SetContour(99);

    clusterwCharge = new TH2D("clusterwCharge", "clusterwCharge", NUM_SPADIC_CHA,0,NUM_SPADIC_CHA,2000,-1,1);
    clusterwCharge->SetXTitle("cluster size");
    clusterwCharge->SetYTitle("relative charge - theo.");
    clusterwCharge->SetContour(99); 

    prf = new TH2D("PRF_2D", "PRF_2D", 30 * padWidth[suid], -1.5 * padWidth[suid] , 1.5 * padWidth[suid], 100, 0, 1);
    prf->SetYTitle("relative pad charge");
    prf->SetXTitle("reco. hit position [mm]");
    prf->SetContour(99); 
    prfProfile = new TProfile("PRF", "PRF", 6 * padWidth[suid], -1.5 * padWidth[suid], 1.5 * padWidth[suid]);
    prfProfile->SetMarkerStyle(24);
    prfProfile->SetMarkerColor(1);
    prfProfile->SetLineColor(1);
    prfProfile->SetYTitle("relative pad charge");
    prfProfile->SetXTitle("reco. hit position [mm]");

    for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {
      name.Form("Ch%02i",ch);
      preCompPulse[ch] = new TH1D("preCompPulse"+name,"preCompPulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      preCompPulse[ch]->SetLineColor(color[ch]);
      preCompPulse[ch]->SetXTitle("time bin []");
      preCompPulse[ch]->SetYTitle("ADC value");
      preCompBaselinePulse[ch] = new TH1D("preCompBaselinePulse"+name,"preCompBaselinePulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      preCompBaselinePulse[ch]->SetLineColor(color[ch]);
      preCompBaselinePulse[ch]->SetXTitle("time bin []");
      preCompBaselinePulse[ch]->SetYTitle("ADC value");
      rawPulse[ch] = new TH1D("rawPulse"+name,"rawPulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      rawPulse[ch]->SetLineColor(color[ch]);
      rawPulse[ch]->SetXTitle("time bin []");
      rawPulse[ch]->SetYTitle("ADC value");
      baselinePulse[ch] = new TH1D("baselinePulse"+name,"baselinePulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      baselinePulse[ch]->SetLineColor(color[ch]);
      baselinePulse[ch]->SetXTitle("time bin []");
      baselinePulse[ch]->SetYTitle("ADC value");
      noisePulse[ch] = new TH1D("noisePulse"+name,"noisePulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      noisePulse[ch]->SetLineColor(color[ch]);
      noisePulse[ch]->SetXTitle("time bin []");
      noisePulse[ch]->SetYTitle("ADC value");
      baselineNoisePulse[ch] = new TH1D("baselineNoisePulse"+name,"baselineNoisePulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      baselineNoisePulse[ch]->SetLineColor(color[ch]);
      baselineNoisePulse[ch]->SetXTitle("time bin []");
      baselineNoisePulse[ch]->SetYTitle("ADC value");
      noiseBaselinePulse[ch] = new TH1D("noiseBaselinePulse"+name,"noiseBaselinePulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      noiseBaselinePulse[ch]->SetLineColor(color[ch]);
      noiseBaselinePulse[ch]->SetXTitle("time bin []");
      noiseBaselinePulse[ch]->SetYTitle("ADC value");
      offSpillPulse[ch] = new TH2I("offSpillPulse"+name,"offSpillPulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE,256,0,256);
      offSpillPulse[ch]->SetLineColor(color[ch]);
      name.Form("Ch%02i",ch);
      lastPulse[ch] = new TH1D("lastPulse"+name,"lastPulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
    }

    averagePulse_e = new TProfile("averagePulse_e","averagePulse_e",SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE); 
    averagePulse_e->SetLineColor(kRed);
    averagePulse_e->SetMarkerStyle(20);
    averagePulse_e->SetMarkerColor(kRed);
    averagePulse_e->SetXTitle("TimeBin");
    averagePulse_e->SetYTitle("ADC value");

    averagePulse_p = new TProfile("averagePulse_p","averagePulse_p",SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
    averagePulse_p->SetLineColor(kBlue);
    averagePulse_p->SetMarkerStyle(22);
    averagePulse_p->SetMarkerColor(kBlue);
    averagePulse_p->SetXTitle("TimeBin");
    averagePulse_p->SetYTitle("ADC value");

    averagePulse_2D_e = new TH2I("averagePulse_2D_e","averagePulse_2D_e",SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE,276,-20,256); 
    //averagePulse_2D_e->SetLineColor(kRed);
    //averagePulse_2D_e->SetMarkerStyle(20);
    //averagePulse_2D_e->SetMarkerColor(kRed);
    averagePulse_2D_e->SetXTitle("TimeBin");
    averagePulse_2D_e->SetYTitle("ADC value");
    averagePulse_2D_e->SetContour(99); 
    averagePulse_2D_p = new TH2I("averagePulse_2D_p","averagePulse_2D_p",SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE,276,-20,256);
    //averagePulse_2D_p->SetLineColor(kBlue);
    //averagePulse_2D_p->SetMarkerStyle(22);
    //averagePulse_2D_p->SetMarkerColor(kBlue);
    averagePulse_2D_p->SetXTitle("TimeBin");
    averagePulse_2D_p->SetYTitle("ADC value");
    averagePulse_2D_p->SetContour(99); 
    averageSignal_2D->SetContour(99);
    averageNoise_2D->SetContour(99);
    preCompSpectrum = new TH1D("preCompSpectrum","preCompSpectrum",nbins,lbin,hbin);
    preCompSpectrum->SetLineColor(2);
    preCompSpectrum->SetXTitle("dE/dx [keV]");
    preCompSpectrum->SetYTitle("normalized counts");
    preCompBaselineSpectrum = new TH1D("preCompBaselineSpectrum","preCompBaselineSpectrum",nbins,lbin,hbin);
    preCompBaselineSpectrum->SetLineColor(7);
    preCompBaselineSpectrum->SetXTitle("dE/dx [keV]");
    preCompBaselineSpectrum->SetYTitle("normalized counts");
    preCompSpectrumElectron = new TH1D("preCompSpectrumElectron","preCompSpectrumElectron",nbins,lbin,hbin);
    preCompSpectrumElectron->SetLineColor(2);
    preCompSpectrumElectron->SetXTitle("dE/dx [keV]");
    preCompSpectrumElectron->SetYTitle("normalized counts");
    preCompBaselineSpectrumElectron = new TH1D("preCompBaselineSpectrumElectron","preCompBaselineSpectrumElectron",nbins,lbin,hbin);
    preCompBaselineSpectrumElectron->SetLineColor(3);
    preCompBaselineSpectrumElectron->SetXTitle("dE/dx [keV]");
    preCompBaselineSpectrumElectron->SetYTitle("normalized counts");
    rawSpectrum = new TH1D("rawSpectrum","rawSpectrum",nbins,lbin,hbin);
    rawSpectrum->SetLineColor(1);
    rawSpectrum->SetXTitle("dE/dx [keV]");
    rawSpectrum->SetYTitle("normalized counts");
    baselineSpectrum = new TH1D("baselineSpectrum","baselineSpectrum",nbins,lbin,hbin);
    baselineSpectrum->SetLineColor(2);
    baselineSpectrum->SetXTitle("dE/dx [keV]");
    baselineSpectrum->SetYTitle("normalized counts");
    noiseSpectrum = new TH1D("noiseSpectrum","noiseSpectrum",nbins,lbin,hbin);
    noiseSpectrum->SetLineColor(3);
    noiseSpectrum->SetXTitle("dE/dx [keV]");
    noiseSpectrum->SetYTitle("normalized counts");
    baselineNoiseSpectrum = new TH1D("baselineNoiseSpectrum","baselineNoiseSpectrum",nbins,lbin,hbin);
    baselineNoiseSpectrum->SetLineColor(4);
    baselineNoiseSpectrum->SetXTitle("dE/dx [keV]");
    baselineNoiseSpectrum->SetYTitle("normalized counts");
    baselineNoiseSpectrumWoOWoU = new TH1D("baselineNoiseSpectrumWoOWoU","baselineNoiseSpectrumWoOWoU",nbins,lbin,hbin);
    baselineNoiseSpectrumWoOWoU->SetLineColor(6);
    baselineNoiseSpectrumWoOWoU->SetXTitle("dE/dx [keV]");
    baselineNoiseSpectrumWoOWoU->SetYTitle("normalized counts");
    baselineNoiseSpectrumElectron = new TH1D("baselineNoiseSpectrumElectron","baselineNoiseSpectrumElectron",nbins,lbin,hbin);
    baselineNoiseSpectrumElectron->SetLineColor(4);
    baselineNoiseSpectrumElectron->SetXTitle("dE/dx [keV]");
    baselineNoiseSpectrumElectron->SetYTitle("normalized counts");
    baselineNoiseSpectrumWoOWoUElectron = new TH1D("baselineNoiseSpectrumWoOWoUElectron","baselineNoiseSpectrumWoOWoUElectron",nbins,lbin,hbin);
    baselineNoiseSpectrumWoOWoUElectron->SetLineColor(6);
    baselineNoiseSpectrumWoOWoUElectron->SetLineStyle(1);
    baselineNoiseSpectrumWoOWoUElectron->SetXTitle("dE/dx [keV]");
    baselineNoiseSpectrumWoOWoUElectron->SetYTitle("normalized counts");
    noiseBaselineSpectrum = new TH1D("noiseBaselineSpectrum","noiseBaselineSpectrum",nbins,lbin,hbin);
    noiseBaselineSpectrum->SetLineColor(800);
    noiseBaselineSpectrum->SetXTitle("dE/dx [keV]");
    noiseBaselineSpectrum->SetYTitle("normalized counts");
    probabilityOverflow = new TH1D("probabilityOverflow","probabilityOverflow",nbins,lbin,hbin);
    probabilityOverflow->SetLineColor(800);
    probabilityOverflow->SetFillColor(800);
    probabilityOverflow->SetFillStyle(3002);
    probabilityUnderflow = new TH1D("probabilityUnderflow","probabilityUnderflow",nbins,lbin,hbin); 
    probabilityUnderflow->SetLineColor(8);
    probabilityUnderflow->SetFillColor(8);
    probabilityUnderflow->SetFillStyle(3005);
    probabilityOverflowE = new TH1D("probabilityOverflowE","probabilityOverflowE",nbins,lbin,hbin);
    probabilityOverflowE->SetLineColor(800);
    probabilityOverflowE->SetFillColor(800);
    probabilityOverflowE->SetFillStyle(3002);
    probabilityUnderflowE = new TH1D("probabilityUnderflowE","probabilityUnderflowE",nbins,lbin,hbin); 
    probabilityUnderflowE->SetLineColor(8);
    probabilityUnderflowE->SetFillColor(8);
    probabilityUnderflowE->SetFillStyle(3005);

    crossTalk = new TProfile("CrossTalk","CrossTalk",6,-0.5,5.5);

    for(Int_t delta = 0; delta < 6; delta++) {
      name.Form("_delta%02i",delta);
      baselineNoiseSpectrumCrossTalkE[delta] = new TH1D("baselineNoiseSpectrumCrossTalkE"+name,"baselineNoiseSpectrumCrossTalkE"+name,nbins,lbin,hbin);
      baselineNoiseSpectrumCrossTalkE[delta]->SetLineColor(12+delta);
      baselineNoiseSpectrumCrossTalkE[delta]->SetXTitle("dE/dx [keV]");
      baselineNoiseSpectrumCrossTalkE[delta]->SetYTitle("normalized counts");
      baselineNoiseSpectrumCrossTalkP[delta] = new TH1D("baselineNoiseSpectrumCrossTalkP"+name,"baselineNoiseSpectrumCrossTalkP"+name,nbins,lbin,hbin);
      baselineNoiseSpectrumCrossTalkP[delta]->SetLineColor(12+delta);
      baselineNoiseSpectrumCrossTalkP[delta]->SetLineStyle(2);
      baselineNoiseSpectrumCrossTalkP[delta]->SetXTitle("dE/dx [keV]");
      baselineNoiseSpectrumCrossTalkP[delta]->SetYTitle("normalized counts");
    }
  }



  TString formula(""), K3formula("-0.245*([0]/[2])+(-20.19*([1]/[2])+0.85)");//K3formula("0.35");;
 
  
  formula.Form(" -1. / (2. * atan(sqrt(%f))) * (atan(sqrt(%f) *tanh(3.14159265 * (-2. + sqrt(%f) ) * (%f + 2.* x * %f) / (8.* %f) )) +  atan(sqrt(%f) *  tanh(3.14159265 * (-2. + sqrt(%f) ) * (%f - 2.* x * %f) / (8.* %f) )) )",calcK3(h[suid]),calcK3(h[suid]),calcK3(h[suid]),padWidth[suid],par,h[suid],calcK3(h[suid]),calcK3(h[suid]),padWidth[suid],par,h[suid]);
  TF1 *mathieson = new TF1("mathieson", formula, -1.5 * padWidth[suid], 1.5 * padWidth[suid]);
  /*
    formula.Form(" -1. / (2. * atan(sqrt([0]))) * (atan(sqrt([0]) *tanh(3.14159265 * (-2. + sqrt([0]) ) * (%f + 2.* x * %f) / (8.* %f) )) +  atan(sqrt([0]) *  tanh(3.14159265 * (-2. + sqrt([0]) ) * (%f - 2.* x * %f) / (8.* %f) )) )",padWidth[suid],par,h[suid],padWidth[suid],par,h[suid]);
    formula.Form(" -1. / (2. * atan(sqrt([0]))) * (atan(sqrt([0]) *tanh(3.14159265 * (-2. + sqrt([0]) ) * (%f + 2.* x * %f) / (8.* [1]) )) +  atan(sqrt([0]) *  tanh(3.14159265 * (-2. + sqrt([0]) ) * (%f - 2.* x * %f) / (8.* [1]) )) )",padWidth[suid],par,padWidth[suid],par);
  */
  formula.Form(" -1. / (2. * atan(sqrt(%s))) * (atan(sqrt(%s) *tanh(3.14159265 * (-2. + sqrt(%s) ) * (%f + 2.* x * %f) / (8.* [0]) )) +  atan(sqrt(%s) *  tanh(3.14159265 * (-2. + sqrt(%s) ) * (%f - 2.* x * %f) / (8.* [0]) )) )",  K3formula.Data(),K3formula.Data(),K3formula.Data(),padWidth[suid],par,K3formula.Data(),K3formula.Data(),padWidth[suid],par);// [0] = h; [1] = r_a; [2] = s
  TF1 *mathiesonFit = new TF1("mathiesonFit", formula, -1.5 * padWidth[suid], 1.5 * padWidth[suid]);
  mathiesonFit->SetLineColor(1);
  mathiesonFit->SetLineStyle(2);
  mathiesonFit->SetLineWidth(1);
  mathieson->SetLineColor(1);
  mathieson->SetLineStyle(1);
  mathieson->SetLineWidth(1);
  const Double_t r_a = 0.002;
  const Double_t s_a = 2.5;
  mathiesonFit->SetParName(0,"h");
  mathiesonFit->SetParameter(0,h[suid]);
  mathiesonFit->SetParLimits(0,1.25,4.0);
  mathiesonFit->SetParName(1,"r_a");
  mathiesonFit->SetParameter(1,r_a);
  mathiesonFit->SetParLimits(1,0.0018,0.0022);
  mathiesonFit->SetParName(2,"s");
  mathiesonFit->SetParameter(2,s_a);
  mathiesonFit->SetParLimits(2,2.4,2.6);


  Int_t nPions = 0;
  Int_t nElectrons = 0;
  Int_t entries = 0;
  TCanvas *c = NULL;
  if (!plotFromFile) {
    name = filename;
    //get_Pb_Ch_Cut(name.ReplaceAll("data2011/TTrees/",""));
    GetPidCuts(name.ReplaceAll("data2012/TTrees/merged",""));

    TFile inputFile(filename,"READ");
    if (!inputFile.IsOpen())
      cout << "file not found: " << filename << endl;
    else
      cout << "file found: " << filename << endl;
    TTree* theTree=NULL;
    TKey* kee=NULL;
    TIter iter(inputFile.GetListOfKeys());
    while ( ( kee=dynamic_cast<TKey*>(iter()) ) !=0 ) {
      theTree = dynamic_cast<TTree*>(kee->ReadObj());
      if (theTree)
        break; // we take first Tree in file, disregarding its name...
    }
    if(theTree==NULL) {
      cout <<"Error: no Tree in file "<< filename.Data() << endl;
      return;
    }
    TCernOct12UnpackEvent* evnt = new TCernOct12UnpackEvent;
    TGo4EventElement* theBase=evnt;
    evnt->synchronizeWithTree(theTree, &theBase);
    entries = theTree->GetEntries();
    printf("%i Events found in TTree\n",entries);
    TMbsCrateEvent* fCrateInputEvent = NULL;
    TSpadicEvent* SpadicInputEvent = NULL;
    TSpadicData* theSpadic = NULL;
    Float_t Ch1(0), Ch2(0), Pb(0);
    Int_t clusterCh = 0;
    if (debug) {
      c = new TCanvas("c","c",1750,700);
      c->Divide(7,2);
    }
    TCanvas *cDebug = NULL;
    if (debug) {
      cDebug = new TCanvas("cDebug","cDebug",800,600);
      cDebug->Divide(2,1);
    }

    //Int_t nAfterMinA = 0;
    //Int_t nAfterBorderPads = 0;
    //Int_t nAfterOverUnder = 0;
    if (debug) entries = 200;
  
    for(Int_t i=0;i<entries;++i) {
      Statusbar(i,entries);
      theTree->GetEntry(i);
      fCrateInputEvent=dynamic_cast<TMbsCrateEvent*>(evnt->GetSubEvent("MBSCRATE"));    
      SpadicInputEvent=dynamic_cast<TSpadicEvent*>(evnt->GetSubEvent("SPADIC"));
      if (fCrateInputEvent->fIsPulser) continue;
      isPion = false;
      isElectron = false;
      isOverFlowEvent = false;
      isUnderFlowEvent = false;
      /* 
      //2011
      Ch1 = fCrateInputEvent->fData1182[1][0];
      Ch2 = fCrateInputEvent->fData1182[0][1];
      Pb  = fCrateInputEvent->fData1182[1][5];
      */
      //2012
      Ch1 = fCrateInputEvent->fData1182[1][0];
      Ch2 = fCrateInputEvent->fData1182[0][0];
      Pb  = fCrateInputEvent->fData1182[1][1];
      getPID(Ch1, Ch2, Pb);

      hCh_1->Fill(Ch1);
      hCh_2->Fill(Ch2);
      hPb->Fill(Pb);

      if ((Ch1 > 350 && Pb > 350) || (Ch2 > 350 && Pb > 350)){

        allCh1_Pb->Fill(Ch1,Pb);
        allCh2_Pb->Fill(Ch2,Pb);
        allCh1_Ch2->Fill(Ch1,Ch2);

        if (isPion || isElectron) {
          Ch1_Pb->Fill(Ch1,Pb);
          Ch2_Pb->Fill(Ch2,Pb);
          Ch1_Ch2->Fill(Ch1,Ch2);
        }
        if (!isPion && !isElectron) {
          noCh1_Pb->Fill(Ch1,Pb);
          noCh2_Pb->Fill(Ch2,Pb);
          noCh1_Ch2->Fill(Ch1,Ch2);
        }
      }
      //for (Int_t suid = 0; suid < /*usedSusibos*/1; suid++){
      theSpadic=dynamic_cast<TSpadicData*>(SpadicInputEvent->getEventElement(usedSuId[suid]));
      if (!theSpadic) { return;}
      if (NULL == theSpadic || !theSpadic) { return;}
      if (isPion || isElectron) {
        for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) { 
          rawPulse[ch]->Reset();
          rawPulse[ch]->GetYaxis()->SetRangeUser(-50,256);
          baselinePulse[ch]->Reset();
          baselinePulse[ch]->GetYaxis()->SetRangeUser(-50,256);
          noisePulse[ch]->Reset();
          noisePulse[ch]->GetYaxis()->SetRangeUser(-50,256);
          baselineNoisePulse[ch]->Reset();
          baselineNoisePulse[ch]->GetYaxis()->SetRangeUser(-50,256);
          noiseBaselinePulse[ch]->Reset();
          noiseBaselinePulse[ch]->GetYaxis()->SetRangeUser(-50,256);
          preCompPulse[ch]->Reset();
          preCompPulse[ch]->GetYaxis()->SetRangeUser(-50,256);
          preCompBaselinePulse[ch]->Reset();
          preCompBaselinePulse[ch]->GetYaxis()->SetRangeUser(-50,256);
          for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) { 
            rawPulse[ch]->Fill(bin, (Double_t)theSpadic->fSpadicPulse[ch][bin]);
            preCompPulse[ch]->Fill(bin, (Double_t)theSpadic->fSpadicCompensated[ch][bin]);
            if ((Int_t)theSpadic->fSpadicPulse[ch][bin] == 0 && bin < SPADIC_TRACE_SIZE-1) 
              isUnderFlowEvent = true;
          
            if ((Int_t)theSpadic->fSpadicPulse[ch][bin] == 255)
              isOverFlowEvent = true;
          
          }      
        }
        if (fast) {
          BaselineCompensation(rawPulse, baselinePulse, baselineDistribution);
          clusterCh = CancelNoise_Cova(baselinePulse, baselineNoisePulse, covaMatixValue,covaMatixValueMaxAmplitude,covaMatixValueHitTime,maxAmplitudeHitTime,noiseDistribution,clusterSize,signalChDistance,averageNoise_2D,averageSignal_2D,NULL,false);
        } else {
          BaselineCompensation(preCompPulse, preCompBaselinePulse, baselineDistribution);
          BaselineCompensation(rawPulse, baselinePulse, baselineDistribution);
          clusterCh = CancelNoise_Cova(rawPulse, noisePulse, covaMatixValue,covaMatixValueMaxAmplitude,covaMatixValueHitTime,maxAmplitudeHitTime,noiseDistribution,clusterSize,signalChDistance,averageNoise_2D,averageSignal_2D,NULL,false);
          clusterCh = CancelNoise_Cova(baselinePulse, baselineNoisePulse, covaMatixValue,covaMatixValueMaxAmplitude,covaMatixValueHitTime,maxAmplitudeHitTime,noiseDistribution,clusterSize,signalChDistance,averageNoise_2D,averageSignal_2D,NULL,false);
          BaselineCompensation(noisePulse, noiseBaselinePulse, baselineDistribution);
        }

        if (isElectron || isPion) 
          for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++)
            maxAmplitudeValue->Fill(baselineNoisePulse[ch]->GetBinContent(baselineNoisePulse[ch]->GetMaximumBin()));

        if(minimumAmplitudeThreshold(baselineNoisePulse) && !BorderPadMax(baselineNoisePulse) && HitTimeWindow(baselineNoisePulse) && minimumIntegralThreshold(baselineNoisePulse)){
          CalcPrf(baselineNoisePulse, prf, prfProfile, padWidth[suid]);
          PrfAna(baselineNoisePulse, clusterwCharge, -1, clusterCh);
        }
        Int_t deltaChMax = -1;     
        Double_t temp = scaling2keV(Get3PadPosIntegralCrossTalk(baselineNoisePulse, lastPulse, &deltaChMax));     
     
        if (isElectron) {
          nElectrons++;
          if(amplitude){         
            if(minimumAmplitudeThreshold(preCompPulse) && !BorderPadMax(preCompPulse) && HitTimeWindow(preCompPulse) && minimumIntegralThreshold(preCompPulse))
              preCompSpectrumElectron->Fill(scaling2keV(Get3PadPosAmplitude(preCompPulse)));
            if(minimumAmplitudeThreshold(preCompBaselinePulse) && !BorderPadMax(preCompBaselinePulse) && HitTimeWindow(preCompBaselinePulse) && minimumIntegralThreshold(preCompBaselinePulse))
              preCompBaselineSpectrumElectron->Fill(scaling2keV(Get3PadPosAmplitude(preCompBaselinePulse)));
            if(minimumAmplitudeThreshold(baselineNoisePulse) && !BorderPadMax(baselineNoisePulse) && HitTimeWindow(baselineNoisePulse) && minimumIntegralThreshold(baselineNoisePulse)){
              baselineNoiseSpectrumElectron->Fill(scaling2keV(Get3PadPosAmplitude(baselineNoisePulse)));
              if (isOverFlowEvent)
                probabilityOverflowE->Fill(Get3PadPosAmplitude(baselineNoisePulse));
              if(isUnderFlowEvent)
                probabilityUnderflowE->Fill(Get3PadPosAmplitude(baselineNoisePulse));
            }
          }
          else{
            if (!fast){
              if(minimumAmplitudeThreshold(preCompPulse) && !BorderPadMax(preCompPulse) && HitTimeWindow(preCompPulse) && minimumIntegralThreshold(preCompPulse))
                preCompSpectrumElectron->Fill(scaling2keV(Get3PadPosIntegral(preCompPulse)));
              if(minimumAmplitudeThreshold(preCompBaselinePulse) && !BorderPadMax(preCompBaselinePulse) && HitTimeWindow(preCompBaselinePulse) && minimumIntegralThreshold(preCompBaselinePulse))
                preCompBaselineSpectrumElectron->Fill(scaling2keV(Get3PadPosIntegral(preCompBaselinePulse)));
            }
            if(minimumAmplitudeThreshold(baselineNoisePulse) && !BorderPadMax(baselineNoisePulse) && HitTimeWindow(baselineNoisePulse) && minimumIntegralThreshold(baselineNoisePulse)){
              baselineNoiseSpectrumElectron->Fill(scaling2keV(Get3PadPosIntegral(baselineNoisePulse)));
              FillAverageSignal(baselineNoisePulse, averagePulse_e, averagePulse_2D_e);
              if (deltaChMax >= 0 && deltaChMax < 6){
                if(minimumAmplitudeThreshold(lastPulse) && !BorderPadMax(lastPulse) && HitTimeWindow(lastPulse) && minimumIntegralThreshold(lastPulse)){
                  baselineNoiseSpectrumCrossTalkE[deltaChMax]->Fill(temp);
                  //crossTalk->Fill(deltaChMax,(temp*scaling2keV - Get3PadPosIntegral(baselineNoisePulse)*scaling2keV) / (temp*scaling2keV) * 100);
                }
              }
              if (isOverFlowEvent)
                probabilityOverflowE->Fill(scaling2keV(Get3PadPosIntegral(baselineNoisePulse)));
              if(isUnderFlowEvent)
                probabilityUnderflowE->Fill(scaling2keV(Get3PadPosIntegral(baselineNoisePulse)));      
              if (!isOverFlowEvent && !isUnderFlowEvent){
                if(amplitude){
                  baselineNoiseSpectrumWoOWoUElectron->Fill(Get3PadPosAmplitude(baselineNoisePulse));     
                }
                else{
                  baselineNoiseSpectrumWoOWoUElectron->Fill(scaling2keV(Get3PadPosIntegral(baselineNoisePulse)));          
                }
              }
            }
          }
        }
        if (isPion){
          nPions++;
          IntegralvsAmplitude->Fill(Get3PadPosIntegral(noiseBaselinePulse), GetMaxAplitude(noiseBaselinePulse));//Get3PadPosAmplitude(noiseBaselinePulse));
          if(minimumAmplitudeThreshold(rawPulse) && !BorderPadMax(rawPulse) && HitTimeWindow(rawPulse) && minimumIntegralThreshold(rawPulse)){
            if(amplitude)
              rawSpectrum->Fill(Get3PadPosAmplitude(rawPulse));
            else
              rawSpectrum->Fill(scaling2keV(Get3PadPosIntegral(rawPulse)));
          }
          if(amplitude){
            if(minimumAmplitudeThreshold(baselinePulse) && !BorderPadMax(baselinePulse) && HitTimeWindow(baselinePulse) && minimumIntegralThreshold(baselinePulse))
              baselineSpectrum->Fill(scaling2keV(Get3PadPosAmplitude(baselinePulse)));
            if(minimumAmplitudeThreshold(noiseBaselinePulse) && !BorderPadMax(noiseBaselinePulse) && HitTimeWindow(noiseBaselinePulse) && minimumIntegralThreshold(noiseBaselinePulse))
              noiseBaselineSpectrum->Fill(scaling2keV(Get3PadPosAmplitude(noiseBaselinePulse)));
            if(minimumAmplitudeThreshold(noisePulse) && !BorderPadMax(noisePulse) && HitTimeWindow(noisePulse) && minimumIntegralThreshold(noisePulse))
              noiseSpectrum->Fill(scaling2keV(Get3PadPosAmplitude(noisePulse)));        
            if(minimumAmplitudeThreshold(preCompPulse) && !BorderPadMax(preCompPulse) && HitTimeWindow(preCompPulse) && minimumIntegralThreshold(preCompPulse))
              preCompSpectrum->Fill(scaling2keV(Get3PadPosAmplitude(preCompPulse)));
            if(minimumAmplitudeThreshold(preCompBaselinePulse) && !BorderPadMax(preCompBaselinePulse) && HitTimeWindow(preCompBaselinePulse) && minimumIntegralThreshold(preCompBaselinePulse))
              preCompBaselineSpectrum->Fill(scaling2keV(Get3PadPosAmplitude(preCompBaselinePulse)));
            if(minimumAmplitudeThreshold(baselineNoisePulse) && !BorderPadMax(baselineNoisePulse) && HitTimeWindow(baselineNoisePulse) && minimumIntegralThreshold(baselineNoisePulse)){
              baselineNoiseSpectrum->Fill(scaling2keV(Get3PadPosAmplitude(baselineNoisePulse)));
              if (!isOverFlowEvent && !isUnderFlowEvent)
                baselineNoiseSpectrumWoOWoU->Fill(Get3PadPosAmplitude(baselineNoisePulse));
              if (isOverFlowEvent)
                probabilityOverflow->Fill(Get3PadPosAmplitude(baselineNoisePulse));
              if(isUnderFlowEvent)
                probabilityUnderflow->Fill(Get3PadPosAmplitude(baselineNoisePulse));
            }
          }
          else{
            if (!fast){
              if(minimumAmplitudeThreshold(baselinePulse) && !BorderPadMax(baselinePulse) && HitTimeWindow(baselinePulse) && minimumIntegralThreshold(baselinePulse))
                baselineSpectrum->Fill(scaling2keV(Get3PadPosIntegral(baselinePulse)));
              if(minimumAmplitudeThreshold(noisePulse) && !BorderPadMax(noisePulse) && HitTimeWindow(noisePulse) && minimumIntegralThreshold(noisePulse))
                noiseSpectrum->Fill(scaling2keV(Get3PadPosIntegral(noisePulse)));
              if(minimumAmplitudeThreshold(noiseBaselinePulse) && !BorderPadMax(noiseBaselinePulse) && HitTimeWindow(noiseBaselinePulse) && minimumIntegralThreshold(noiseBaselinePulse)){
                noiseBaselineSpectrum->Fill(scaling2keV(Get3PadPosIntegral(noiseBaselinePulse)));
              }
              if(minimumAmplitudeThreshold(preCompPulse) && !BorderPadMax(preCompPulse) && HitTimeWindow(preCompPulse) && minimumIntegralThreshold(preCompPulse))
                preCompSpectrum->Fill(scaling2keV(Get3PadPosIntegral(preCompPulse)));
              if(minimumAmplitudeThreshold(preCompBaselinePulse) && !BorderPadMax(preCompBaselinePulse) && HitTimeWindow(preCompBaselinePulse) && minimumIntegralThreshold(preCompBaselinePulse))
                preCompBaselineSpectrum->Fill(scaling2keV(Get3PadPosIntegral(preCompBaselinePulse)));
            }
            if(minimumAmplitudeThreshold(baselineNoisePulse) && !BorderPadMax(baselineNoisePulse) && HitTimeWindow(baselineNoisePulse) && minimumIntegralThreshold(baselineNoisePulse)){
              baselineNoiseSpectrum->Fill(scaling2keV(Get3PadPosIntegral(baselineNoisePulse)));
              keV2AmplitudeProfile->Fill(scaling2keV(Get3PadPosIntegral(noiseBaselinePulse)), GetMaxAplitude(noiseBaselinePulse));
              keV2Amplitude->Fill(scaling2keV(Get3PadPosIntegral(noiseBaselinePulse)), GetMaxAplitude(noiseBaselinePulse));
              //IntegralvsAmplitude->Fill(Get3PadPosIntegral(noiseBaselinePulse), GetMaxAplitude(noiseBaselinePulse));//Get3PadPosAmplitude(noiseBaselinePulse));
              FillAverageSignal(baselineNoisePulse, averagePulse_p, averagePulse_2D_p);
              if (deltaChMax >= 0 && deltaChMax < 6){
                if(minimumAmplitudeThreshold(lastPulse) && !BorderPadMax(lastPulse) && HitTimeWindow(lastPulse) && minimumIntegralThreshold(lastPulse)){
                  baselineNoiseSpectrumCrossTalkP[deltaChMax]->Fill(temp);
                  crossTalk->Fill(deltaChMax,(temp - scaling2keV(Get3PadPosIntegral(baselineNoisePulse))) / (temp) * 100);
                }
              }
              if (!isOverFlowEvent && !isUnderFlowEvent)
                baselineNoiseSpectrumWoOWoU->Fill(scaling2keV(Get3PadPosIntegral(baselineNoisePulse)));
              if (isOverFlowEvent)
                probabilityOverflow->Fill(scaling2keV(Get3PadPosIntegral(baselineNoisePulse)));
              if(isUnderFlowEvent)
                probabilityUnderflow->Fill(scaling2keV(Get3PadPosIntegral(baselineNoisePulse)));
            }
          
          }
        }
        if (debug) {
          if (i%100000 == 0) {
            for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {     
              c->cd(1);
              if(ch==0)
                rawPulse[ch]->DrawCopy();
              else
                rawPulse[ch]->DrawCopy("same");
              c->cd(2);
              if(ch==0)
                baselinePulse[ch]->DrawCopy();
              else
                baselinePulse[ch]->DrawCopy("same");
              c->cd(3);
              if(ch==0)
                baselineNoisePulse[ch]->DrawCopy();
              else
                baselineNoisePulse[ch]->DrawCopy("same");
              c->cd(4);
              if(ch==0)
                noisePulse[ch]->DrawCopy();
              else
                noisePulse[ch]->DrawCopy("same");
              c->cd(5);
              if(ch==0)
                noiseBaselinePulse[ch]->DrawCopy();
              else
                noiseBaselinePulse[ch]->DrawCopy("same");
              c->cd(6);
              if(ch==0)
                preCompPulse[ch]->DrawCopy();
              else
                preCompPulse[ch]->DrawCopy("same");
              c->cd(7);
              if(ch==0)
                preCompBaselinePulse[ch]->DrawCopy();
              else
                preCompBaselinePulse[ch]->DrawCopy("same");

              rawPulse[ch]->Reset();
              baselinePulse[ch]->Reset();
              baselineNoisePulse[ch]->Reset();
            }
          }
          c->cd(8)->SetLogy(1);
          rawSpectrum->DrawCopy();
          c->cd(9)->SetLogy(1);
          baselineSpectrum->DrawCopy();
          c->cd(10)->SetLogy(1);
          baselineNoiseSpectrum->DrawCopy();
          c->cd(11)->SetLogy(1);
          noiseSpectrum->DrawCopy();
          c->cd(12)->SetLogy(1);
          noiseBaselineSpectrum->DrawCopy();
          c->cd(13)->SetLogy(1);
          preCompSpectrum->DrawCopy();
          c->cd(14)->SetLogy(1);
          preCompBaselineSpectrum->DrawCopy();
          c->Update();
          cDebug->cd(1);
          clusterSize->DrawCopy();
          cDebug->cd(2);
          signalChDistance->DrawCopy();
          cDebug->Update();
          cDebug->SaveAs("pics/Debug.pdf");
        }
    
        if(minimumAmplitudeThreshold(baselineNoisePulse) && !BorderPadMax(baselineNoisePulse) && HitTimeWindow(baselineNoisePulse) && minimumIntegralThreshold(baselineNoisePulse)){
          for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++){
            for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {
              lastPulse[ch]->SetBinContent(lastPulse[ch]->FindBin(bin),baselineNoisePulse[ch]->GetBinContent(baselineNoisePulse[ch]->FindBin(bin)));
            }
          }
        }
      }
    }
  } // if !plotFromFile

  TFile *sim = new TFile(simfile.Data(),"READ");
  TH1D* simSpectrum = (TH1D*)sim->Get("hPi");
 
  simSpectrum->SetXTitle("dE/dx [keV]");
  simSpectrum->SetYTitle("normalized counts");
  simSpectrum->SetStats(false);
  simSpectrum->SetLineStyle(1);
  simSpectrum->SetLineColor(1);
  simSpectrum->Scale(/*baselineNoiseSpectrum->GetEntries()*/1./(Double_t)simSpectrum->Integral());
  simSpectrum->GetXaxis()->SetRangeUser(-0.5,70);
  if (debug) {
    c->SaveAs("pics/Spectra_overview.pdf");
    c->Clear();
    c->cd(1)->SetLogy(1);
    rawSpectrum->GetYaxis()->SetRangeUser(0.5,1e4);
    rawSpectrum->DrawCopy();
    c->cd(2)->SetLogy(1);
    baselineSpectrum->DrawCopy();
    c->cd(3)->SetLogy(1);
    baselineNoiseSpectrum->DrawCopy();
    c->cd(4)->SetLogy(1);
    noiseSpectrum->DrawCopy();
    c->cd(5)->SetLogy(1);
    noiseBaselineSpectrum->DrawCopy();
    c->cd(6)->SetLogy(1);
  }
  TCanvas *cResults = new TCanvas("cResults","cResults",1600,600);
  //cResults->Divide(1,1);
  TPad *pad1 = new TPad("pad1","top pad",0,0.35,0.5,1);
  //ad *pad1 = new TPad("pad1","top pad",0,0,1,1);
  pad1->SetBottomMargin(0);
  pad1->Draw();
  
  TPad *pad2 = new TPad("pad2","bottom pad",0,0,0.5,0.35);
  pad2->SetTopMargin(0);
  pad2->SetBottomMargin(0.2);
  pad2->Draw();
  TPad *pad3 = new TPad("pad3","top pad",0.5,0.35,1,1);
  //ad *pad1 = new TPad("pad1","top pad",0,0,1,1);
  pad3->SetBottomMargin(0);
  pad3->Draw();
  
  TPad *pad4 = new TPad("pad4","bottom pad",0.5,0,1,0.35);
  pad4->SetTopMargin(0);
  pad4->SetBottomMargin(0.2);
  pad4->Draw();
  
  pad1->cd()->SetLogy(1);
  TLegend *leg = new TLegend(0.15,0.5,0.85,0.85);
  leg->SetNColumns(2);
  leg->SetFillStyle(0);
  leg->SetTextSize(0.035);
  leg->SetLineStyle(0);
  leg->SetLineColor(0);
  TString calib = "";
  calib.Form(" (#sigma:%.2f MPV:%.2f #Chi^{2}/NDF:%.2f m:%.3fE-3 b:%.3f)",sigma_adc,MPV_adc,ChiSqr_adc,
             //(MPV_kev-sigma_kev)/(MPV_adc-sigma_adc),MPV_kev-(MPV_kev-sigma_kev)/(MPV_adc-sigma_adc)*MPV_adc, // 1)
             sigma_kev/sigma_adc*1000,MPV_kev-sigma_kev/sigma_adc*MPV_adc // 2)
             );//calibrationParameter[3],calibrationParameter[2]/calibrationParameter[0],calibrationParameter[1]);
  leg->SetHeader("#pi^{-} spectra, " + usedSuName[suid] + ", p=3GeV/c, 12mm XeCO_{2} " + calib);
  leg->AddEntry(simSpectrum,"simulated","l");
  leg->AddEntry((TObject*)0, "","");
  if (simple) {
    leg->AddEntry(baselineNoiseSpectrum,"measured: 1.baseline 2.noise","l");
    leg->AddEntry(baselineNoiseSpectrumWoOWoU,"measured: 1.baseline 2.noise(woO && woU)","l");
    //leg->AddEntry(preCompSpectrum,"measured: MS noise","l");
    //leg->AddEntry(preCompBaselineSpectrum,"measured: MS noise + baseline","l");
    //leg->AddEntry((TObject*)0, "               woO && woU", "");
    //leg->AddEntry(noiseBaselineSpectrum,"measured: 1. noise 2. baseline","l");
    leg->AddEntry(probabilityOverflow,"overflow probability","f");
    leg->AddEntry(probabilityUnderflow,"underflow probability","f");
  } else {
    leg->AddEntry(baselineNoiseSpectrum,"measured: FFM noise + baseline","l");
    leg->AddEntry(baselineNoiseSpectrumWoOWoU,"measured: FFM noise + baseline (woO && woU)","l");
    leg->AddEntry(preCompSpectrum,"measured: MS noise","l");
    leg->AddEntry(preCompBaselineSpectrum,"measured: MS noise + baseline","l");
    //leg->AddEntry((TObject*)0, "               woO && woU", "");
    //leg->AddEntry(noiseBaselineSpectrum,"measured: 1. noise 2. baseline","l");
    leg->AddEntry(probabilityOverflow,"overflow probability","f");
    leg->AddEntry(probabilityUnderflow,"underflow probability","f");
  }

  simSpectrum->GetYaxis()->SetRangeUser(5e-6, 10);
  simSpectrum->GetXaxis()->SetLabelSize(0.04);
  simSpectrum->GetYaxis()->SetLabelSize(0.04);
  simSpectrum->GetXaxis()->SetTitleSize(0.045);
  simSpectrum->GetYaxis()->SetTitleSize(0.045);
  simSpectrum->GetYaxis()->SetTitleOffset(0.9);
  simSpectrum->GetXaxis()->SetTitleOffset(0.5);

  simSpectrum->DrawCopy("");
  if (!plotFromFile){
    probabilityOverflow->Scale(1./(Double_t)/*probabilityOverflow*/baselineNoiseSpectrum->Integral());
    probabilityUnderflow->Scale(1./(Double_t)/*probabilityUnderflow*/baselineNoiseSpectrum->Integral());
  }
  probabilityUnderflow->DrawCopy("same");
  probabilityOverflow->DrawCopy("same");
  //rawSpectrum->DrawCopy();
  //baselineSpectrum->DrawCopy("same");
  if (!plotFromFile)
    baselineNoiseSpectrum->Scale(1./(Double_t)baselineNoiseSpectrum->Integral());
  baselineNoiseSpectrum->DrawCopy("same");
  if (!plotFromFile)
    baselineNoiseSpectrumWoOWoU->Scale(1./(Double_t)baselineNoiseSpectrumWoOWoU->Integral());
  baselineNoiseSpectrumWoOWoU->DrawCopy("same");
  //noiseSpectrum->DrawCopy("same");
  //noiseBaselineSpectrum->Scale(1./(Float_t)noiseBaselineSpectrum->GetEntries());
  //noiseBaselineSpectrum->DrawCopy("same");
  if (!simple) {
    if (!plotFromFile)
      preCompSpectrum->Scale(1./(Double_t)preCompSpectrum->Integral());
    preCompSpectrum->DrawCopy("same");
    if (!plotFromFile)
      preCompBaselineSpectrum->Scale(1./(Double_t)preCompBaselineSpectrum->Integral());
    preCompBaselineSpectrum->DrawCopy("same");
  }
  simSpectrum->DrawCopy("same");
  leg->Draw("same");
  
  //TH1D *deltaSim = (TH1D*)simSpectrum->Clone("deltaSim");
  //deltaSim->Reset();
  TH1D *deltaSim = (TH1D*)baselineNoiseSpectrum->Clone("deltaSim");
  CalcDeltaSim(deltaSim, simSpectrum/*, baselineNoiseSpectrum*/);
  deltaSim->SetLineColor(4);
  deltaSim->SetLineStyle(1);
  deltaSim->SetMarkerStyle(20);
  deltaSim->SetMarkerColor(4);

  //TH1D *deltaSim2 = (TH1D*)simSpectrum->Clone("deltaSim2");
  //deltaSim2->Reset();
  TH1D *deltaSim2 = (TH1D*)baselineNoiseSpectrumWoOWoU->Clone("deltaSim2"); 
  CalcDeltaSim(deltaSim2, simSpectrum/*, baselineNoiseSpectrumWoOWoU*/);
  deltaSim2->SetLineColor(6);
  deltaSim2->SetLineStyle(1);
  deltaSim2->SetMarkerStyle(20);
  deltaSim2->SetMarkerColor(6);

  //TH1D *deltaSimMS = (TH1D*)simSpectrum->Clone("deltaSimMS");
  //deltaSimMS->Reset();
  if (!simple) {
    TH1D *deltaSimMS = (TH1D*)preCompSpectrum->Clone("deltaSimMS"); 
    CalcDeltaSim(deltaSimMS, simSpectrum/*, preCompSpectrum*/);
    deltaSimMS->SetLineColor(2);
    deltaSimMS->SetLineStyle(1);
    deltaSimMS->SetMarkerStyle(20);
    deltaSimMS->SetMarkerColor(2);

    //TH1D *deltaSimMS2 = (TH1D*)simSpectrum->Clone("deltaSimMS2");
    //deltaSimMS2->Reset();
    TH1D *deltaSimMS2 = (TH1D*)preCompBaselineSpectrum->Clone("deltaSimMS2");
    CalcDeltaSim(deltaSimMS2, simSpectrum/*, preCompBaselineSpectrum*/);

    deltaSimMS2->SetLineColor(7);
    deltaSimMS2->SetLineStyle(1);
    deltaSimMS2->SetMarkerStyle(20);
    deltaSimMS2->SetMarkerColor(7);
    pad2->cd();
    deltaSimMS->GetYaxis()->SetRangeUser(0.1, 2.1);
    deltaSimMS->DrawCopy("same");
    deltaSimMS2->GetYaxis()->SetRangeUser(0.1, 2.1);
    deltaSimMS2->DrawCopy("same");
  }

  pad2->cd();
  deltaSim->GetXaxis()->SetRangeUser(-0.5,70);
  deltaSim->GetYaxis()->SetRangeUser(0.1, 2.1);
  deltaSim->DrawCopy("");
  deltaSim2->GetYaxis()->SetRangeUser(0.1, 2.1);
  deltaSim2->DrawCopy("same");





  pad3->cd()->SetLogy(1);
  //TFile dEdxFile("Electrons_dEdx_3GeV_12mm_XeCo2.root","READ");
  TFile dEdxFile(/*"Electrons_3GeV_12mm_XeCo2.root"*//*"data2012/root/sim/Sim_xeco20_12mm_03GeV.root"*/simfile.Data(),"READ");
  TString Tr_simfile = "data2012/root/sim/Sim_xeco20_12mm_03GeV_TR.root";
  Tr_simfile.Form("data2012/root/sim/Sim_xeco20_12mm_%02iGeV_TR.root",momentum);
  TFile TR_simFile(Tr_simfile.Data(),"READ");
  TH1D *hElDedx = (TH1D*)dEdxFile.Get("hElDedx");
  TH1D *hElTr;
  TH1D *hElTrAbs;
  TString radiator2;
  
  TH1D *spectrum; 
  TH1D* spectrumAbs;
  Float_t attenuation = 1.0;
  if (radiator == "") { // No radiator
    hElTrAbs = (TH1D*)dEdxFile.Get("hElDedx");
  }
  else {
    if (radiator == "Bshort" || radiator == "B" || radiator == "B++" || radiator == "C" || radiator == "D" || 
        radiator == "E" || radiator == "F" || radiator == "K" ||  radiator == "Kshort" || radiator == "K++" || 
        radiator == "FFM0.5x100" || radiator == "FFM0.5x200" || radiator == "FFM0.5x350" || 
        radiator == "FFM0.7x100" || radiator == "FFM0.7x200" || radiator == "FFM0.7x350" || 
        radiator == "FFM1.2x100" || radiator == "FFM1.2x200" || radiator == "FFM1.2x350"){
      spectrumAbs = (TH1D*)TR_simFile.Get("Absorption_("+radiator+")");
      spectrum = (TH1D*)TR_simFile.Get("Production_("+radiator+")");
      if (radiator == "B")
        attenuation = 0.65;//ok
      if (radiator == "B++")
        attenuation = 0.65;//
      if (radiator == "Bshort")
        attenuation = 0.65;//
      if (radiator == "K")
        attenuation = 0.65;//
      if (radiator == "K++")
        attenuation = 0.65;//
      if (radiator == "Kshort")
        attenuation = 0.65;//
      if (radiator == "C")
        attenuation = 0.30;//ok
      if (radiator == "D")
        attenuation = 0.45;//ok
      if (radiator == "E")
        attenuation = 0.60;//ok
      if (radiator == "F")
        attenuation = 0.65;//ok
      if (radiator == "FFM0.5x100")
        attenuation = 0.75;//0.65;//
      if (radiator == "FFM0.5x200")
        attenuation = 0.75;//0.75;//
      if (radiator == "FFM0.5x350")
        attenuation = 0.75;//0.65;//
      if (radiator == "FFM0.7x100")
        attenuation = 0.75;//0.75;//
      if (radiator == "FFM0.7x200")
        attenuation = 0.75;//0.80;//
      if (radiator == "FFM0.7x350")
        attenuation = 0.75;//0.65;//
      if (radiator == "FFM1.2x100")
        attenuation = 0.75;//0.75;//
      if (radiator == "FFM1.2x200")
        attenuation = 0.75;//0.75;//
      if (radiator == "FFM1.2x350")
        attenuation = 0.75;//0.65;//
    }
    else if (radiator == "A" || radiator == "G10" || radiator == "G20" || radiator == "G30" || 
             radiator == "H25" || radiator == "H" || radiator == "H++" || radiator == "Iu"  || 
             radiator == "Ip" || radiator == "J" || radiator == "L" || radiator == "M" || 
             radiator == "N25" || radiator == "O5x5" || radiator == "HF110" || radiator == "3xEF700"){
      TString irradiator;
      if (radiator == "A"){
        irradiator = "C";
        attenuation = 0.40;//ok
      }
      if (radiator == "G10"){
        irradiator = "C";
        attenuation = 0.35;//ok
      }
      if (radiator == "G20"){
        irradiator = "C";
        attenuation = 0.50;//ok
      }
      if (radiator == "G30"){
        irradiator = "C";
        attenuation = 0.78;//0.65;//ok
      }
      if (radiator == "H25"){
        irradiator = "C";
        attenuation = 0.35;//
      }
      if (radiator == "H"){
        irradiator = "C";
        attenuation = 0.78;//0.75;//0.65;//ok
      }
      if (radiator == "H++"){
        irradiator = "C";
        attenuation = 0.60;//ok
      }
      if (radiator == "Ip"){
        irradiator = "C";
        attenuation = 0.60;//ok
      }
      if (radiator == "Iu"){
        irradiator = "C";
        attenuation = 0.55;//ok
      }
      if (radiator == "J"){
        irradiator = "C";
        attenuation = 0.20;//ok
      }
      if (radiator == "L"){
        irradiator = "C";
        attenuation = 0.60;//0.50;//0.45;//
      }
      if (radiator == "M"){
        irradiator = "C";
        attenuation = 0.50;//
      }
      if (radiator == "N25"){
        irradiator = "C";
        attenuation = 0.42;//0.40;//0.55;//
      }
      if (radiator == "O5x5"){
        irradiator = "C";
        attenuation = 0.39;//0.38;//0.35;//
      }
      if (radiator == "HF110"){
        irradiator = "C";
        attenuation = 0.65;//0.35;//0.40;//
      }
      if (radiator == "3xEF700"){
        irradiator = "C";
        attenuation = 0.40;//0.35;//0.30;//
      }
      spectrumAbs = (TH1D*)TR_simFile.Get("Absorption_("+/*ir*/radiator+")");
      spectrum = (TH1D*)TR_simFile.Get("Production_("+/*ir*/radiator+")");
    }
    else {
      printf("ERROR::Main: no or wrong radiator type");
      return;
    }
  
    hElTr = (TH1D*)hElDedx->Clone("hElTr_"+TString(spectrum->GetTitle()));
    hElTr->Reset();
    hElTrAbs = (TH1D*)hElDedx->Clone("hElTrAbs_"+TString(spectrumAbs->GetTitle()));
    hElTrAbs->Reset();
    TRandom2 *rEfficiency = new TRandom2();
 

    Int_t nElec = 0.2E6;
    radiator2.Form("%s (attenuated to %.1f%%)",radiator.Data(),attenuation*100);
    Double_t meanPhot = spectrum->Integral() * spectrum->GetBinWidth(1) * attenuation;
    for (Int_t iElec = 0; iElec < nElec; iElec++) {
      Int_t nPhot = rEfficiency->Poisson(meanPhot);
      Float_t dEdx = hElDedx->GetRandom();
      Float_t Tr = 0;
      for (Int_t iPhot = 0; iPhot < nPhot; iPhot++) {
        Tr += spectrum->GetRandom();
      }
      hElTr->Fill(dEdx+Tr);
    }
    meanPhot = spectrumAbs->Integral() * spectrumAbs->GetBinWidth(1) * attenuation;
    for (Int_t iElec = 0; iElec < nElec; iElec++) {
      Int_t nPhot = rEfficiency->Poisson(meanPhot);
      Float_t dEdx = hElDedx->GetRandom();
      Float_t Tr = 0;
      for (Int_t iPhot = 0; iPhot < nPhot; iPhot++) {
        Tr += spectrumAbs->GetRandom();
      }
      hElTrAbs->Fill(dEdx+Tr);
    }
  }
  //TFile *sim2 = new TFile("Electrons_3GeV_12mm_XeCo2.root","READ");
  //TH1D* simSpectrumE = (TH1D*)sim2->Get("ElectronSpectrum_Abs_(D)");
  TH1D* simSpectrumE = hElTrAbs;
  simSpectrumE->SetXTitle("dE/dx [keV]");
  simSpectrumE->SetYTitle("normalized counts");
  simSpectrumE->SetStats(false);
  simSpectrumE->SetLineStyle(1);
  simSpectrumE->SetLineColor(1);
  simSpectrumE->Scale(/*baselineNoiseSpectrum->GetEntries()*/1./(Double_t)simSpectrumE->Integral());
  simSpectrumE->GetXaxis()->SetRangeUser(-0.5,70);
  TLegend *legE = new TLegend(0.15,0.5,0.85,0.85);
  legE->SetNColumns(2);
  legE->SetFillStyle(0);
  legE->SetTextSize(0.035);
  legE->SetLineStyle(0);
  legE->SetLineColor(0);
  legE->SetHeader("e^{-} spectra, " + usedSuName[suid] + ", p=3GeV/c, 12mm XeCO_{2}"); 
  legE->AddEntry(simSpectrumE,"simulated radiator: "+radiator2,"l");
  legE->AddEntry((TObject*)0, "","");
  if (simple) {
    legE->AddEntry(baselineNoiseSpectrumElectron,"measured: 1.baseline 2.noise","l");
    legE->AddEntry(baselineNoiseSpectrumWoOWoUElectron,"measured: 1.baseline 2.noise(woO && woU)","l");
    //legE->AddEntry(preCompSpectrumElectron,"measured: MS noise","l");
    //legE->AddEntry(preCompBaselineSpectrumElectron,"measured: MS noise + baseline","l");
    //legE->AddEntry((TObject*)0, "               woO && woU", "");
    //legE->AddEntry(noiseBaselineSpectrum,"measured: 1. noise 2. baseline","l");
    legE->AddEntry(probabilityOverflowE,"overflow probability","f");
    legE->AddEntry(probabilityUnderflowE,"underflow probability","f");
  } else {
    legE->AddEntry(baselineNoiseSpectrumElectron,"measured: baseline + FFM noise","l");
    legE->AddEntry(baselineNoiseSpectrumWoOWoUElectron,"measured: baseline + FFM noise (woO && woU)","l");
    legE->AddEntry(preCompSpectrumElectron,"measured: MS noise","l");
    legE->AddEntry(preCompBaselineSpectrumElectron,"measured: MS noise + baseline","l");
    //legE->AddEntry((TObject*)0, "               woO && woU", "");
    //legE->AddEntry(noiseBaselineSpectrum,"measured: 1. noise 2. baseline","l");
    legE->AddEntry(probabilityOverflowE,"overflow probability","f");
    legE->AddEntry(probabilityUnderflowE,"underflow probability","f");
  }
  simSpectrumE->GetYaxis()->SetRangeUser(5e-6, 10);
  simSpectrumE->GetXaxis()->SetLabelSize(0.04);
  simSpectrumE->GetYaxis()->SetLabelSize(0.04);
  simSpectrumE->GetXaxis()->SetTitleSize(0.045);
  simSpectrumE->GetYaxis()->SetTitleSize(0.045);
  simSpectrumE->GetYaxis()->SetTitleOffset(0.9);
  simSpectrumE->GetXaxis()->SetTitleOffset(0.5);
  simSpectrumE->DrawCopy();
  if (!plotFromFile){
    probabilityOverflowE->Scale(1./(Double_t)/*probabilityOverflow*/baselineNoiseSpectrumElectron->Integral());
    probabilityUnderflowE->Scale(1./(Double_t)/*probabilityUnderflow*/baselineNoiseSpectrumElectron->Integral());
  }
  probabilityUnderflowE->DrawCopy("same");
  probabilityOverflowE->DrawCopy("same");
  if (!plotFromFile)
    baselineNoiseSpectrumElectron->Scale(1./(Double_t)baselineNoiseSpectrumElectron->Integral());
  baselineNoiseSpectrumElectron->DrawCopy("same");
  if (!plotFromFile)
    baselineNoiseSpectrumWoOWoUElectron->Scale(1./(Double_t)baselineNoiseSpectrumWoOWoUElectron->Integral());
  baselineNoiseSpectrumWoOWoUElectron->DrawCopy("same");
  if (!simple) {
    if (!plotFromFile)
      preCompSpectrumElectron->Scale(1./(Double_t)preCompBaselineSpectrumElectron->Integral());
    preCompSpectrumElectron->DrawCopy("same");
    if (!plotFromFile)
      preCompBaselineSpectrumElectron->Scale(1./(Double_t)preCompBaselineSpectrumElectron->Integral());
    preCompBaselineSpectrumElectron->DrawCopy("same");
  }
  simSpectrumE->DrawCopy("same");
  legE->Draw("same");


  //TH1D *deltaSimE = (TH1D*)simSpectrumE->Clone("deltaSimE");
  //deltaSimE->Reset();
  TH1D *deltaSimE = (TH1D*)baselineNoiseSpectrumElectron->Clone("deltaSimE");
  CalcDeltaSim(deltaSimE, simSpectrumE/*, baselineNoiseSpectrumElectron*/);
  deltaSimE->SetLineColor(4);
  deltaSimE->SetLineStyle(1);
  deltaSimE->SetMarkerStyle(20);
  deltaSimE->SetMarkerColor(4);

  //TH1D *deltaSimE2 = (TH1D*)simSpectrumE->Clone("deltaSimE2");
  //deltaSimE2->Reset();
  TH1D *deltaSimE2 = (TH1D*)baselineNoiseSpectrumWoOWoUElectron->Clone("deltaSimE2");
  CalcDeltaSim(deltaSimE2, simSpectrumE/*, baselineNoiseSpectrumWoOWoUElectron*/);
  deltaSimE2->SetLineColor(6);
  deltaSimE2->SetLineStyle(1);
  deltaSimE2->SetMarkerStyle(20);
  deltaSimE2->SetMarkerColor(6);

  if (!simple) {
    //TH1D *deltaSimMSE = (TH1D*)simSpectrumE->Clone("deltaSimMSE");
    //deltaSimMSE->Reset();
    TH1D *deltaSimMSE = (TH1D*)preCompSpectrumElectron->Clone("deltaSimMSE");
    CalcDeltaSim(deltaSimMSE, simSpectrumE/*, preCompSpectrumElectron*/);
    deltaSimMSE->SetLineColor(2);
    deltaSimMSE->SetLineStyle(1);
    deltaSimMSE->SetMarkerStyle(20);
    deltaSimMSE->SetMarkerColor(2);
 
    //TH1D *deltaSimMSE2 = (TH1D*)simSpectrumE->Clone("deltaSimMSE2");
    //deltaSimMSE2->Reset();
    TH1D *deltaSimMSE2 = (TH1D*)preCompBaselineSpectrumElectron->Clone("deltaSimMSE2");
    CalcDeltaSim(deltaSimMSE2, simSpectrumE/*, preCompBaselineSpectrumElectron*/);
    deltaSimMSE2->SetLineColor(7);
    deltaSimMSE2->SetLineStyle(1);
    deltaSimMSE2->SetMarkerStyle(20);
    deltaSimMSE2->SetMarkerColor(7);
    pad4->cd();
    deltaSimMSE->GetYaxis()->SetRangeUser(0.1, 2.1);
    deltaSimMSE->DrawCopy("same");
    deltaSimMSE2->GetYaxis()->SetRangeUser(0.1, 2.1);
    deltaSimMSE2->DrawCopy("same");
  }

  pad4->cd();
  deltaSimE->GetYaxis()->SetRangeUser(0.1, 2.1);
  deltaSimE->GetXaxis()->SetRangeUser(-0.5,70);
  deltaSimE->DrawCopy("");
  deltaSimE2->GetYaxis()->SetRangeUser(0.1, 2.1);
  deltaSimE2->DrawCopy("same");

  TLine *l2 = new TLine(0,1,70,1);
  l2->SetLineStyle(2);
  pad4->cd();
  l2->Draw("same");
  pad2->cd();
  l2->Draw("same");


 
  //cResults->Update();
  cResults->SaveAs(outpic+".pdf");
  cResults->SaveAs(outpic+".png");



  TCanvas *cPRF = new TCanvas("cPRF","cPRF",800,600);
  TLegend* lPRF = new TLegend(0.65,0.65,0.85,0.85);
  lPRF->SetFillStyle(0);
  lPRF->SetTextSize(0.025);
  lPRF->SetLineStyle(0);
  lPRF->SetLineColor(0);
  TString values;
  lPRF->AddEntry(prfProfile,"average PRF","lp");
  lPRF->AddEntry(mathieson,"mathieson formula","l");
  lPRF->AddEntry((TObject*)0," K_{3}=0.35","");
  prfProfile->Fit("mathiesonFit","QR0");
  values.Form(" K_{3}=%.3f",-0.245 * (mathiesonFit->GetParameter(0)/mathiesonFit->GetParameter(2)) + (-20.19*(mathiesonFit->GetParameter(1)/mathiesonFit->GetParameter(2))+0.85));//mathiesonFit->GetParameter(0));
  lPRF->AddEntry(mathiesonFit,"mathieson fit","l");
  lPRF->AddEntry((TObject*)0,values,"");
  values.Form(" h=%.3f",mathiesonFit->GetParameter(0));
  lPRF->AddEntry((TObject*)0,values,"");
  cPRF->cd()->SetLogz(1);
  prf->DrawCopy("colz");
  prfProfile->Fit("mathiesonFit","QR0");
  mathiesonFit->DrawCopy("same");
  mathieson->DrawCopy("same");
  prfProfile->DrawCopy("same");

  lPRF->Draw("same");
  cPRF->SaveAs(outpic+"_PRF.pdf");
  cPRF->SaveAs(outpic+"_PRF.png");

  const Int_t nChambers = 12;
  TCanvas *cEfficiency = new TCanvas("cEfficiency","cEfficiency",800,600);
  TProfile* pEfficiency = new TProfile("pEfficiency","pEfficiency",12,0.5,12.5);
  pEfficiency->GetYaxis()->SetRangeUser(0.0001,100.5);
  TGraphErrors* hEfficiency = new TGraphErrors(nChambers);
  TGraphErrors* fitEfficiency = new TGraphErrors(nChambers);
  TLegend* leg3 = new TLegend(0.15,0.15,0.35,0.35);
  leg3->SetFillStyle(0);
  leg3->SetTextSize(0.025);
  leg3->SetLineStyle(0);
  leg3->SetLineColor(0);
  leg3->SetHeader(usedSuName[suid] + ", radiator: "+radiator2+", p=3GeV/c, 12mm XeCO_{2}");
  TMultiGraph *g1 = new TMultiGraph();    
  //TMultiGraph *g2 = new TMultiGraph();
  //TFile *output = new TFile(outfilename,"UPDATE");
  output->cd();
  if (true){
    if (simple) {
      TString names;
      for(Int_t delta = 0; delta < 6; delta++) {
        hEfficiency->GetYaxis()->SetRangeUser(1e-2,1e2);
        hEfficiency->SetMarkerColor(baselineNoiseSpectrumCrossTalkP[delta]->GetLineColor());
        hEfficiency->SetLineColor(baselineNoiseSpectrumCrossTalkP[delta]->GetLineColor());
        CalcPionEfficiency(baselineNoiseSpectrumCrossTalkE[delta], baselineNoiseSpectrumCrossTalkP[delta], hEfficiency, fitEfficiency, nChambers, debug);
        names.Form("Efficiency_Crosstalk_delta%02i",delta);
        hEfficiency->Write(names, TObject::kOverwrite);
        fitEfficiency->Write(names + TString("_fit"), TObject::kOverwrite);
      }
     
      CalcPionEfficiency(simSpectrumE, simSpectrum, hEfficiency, fitEfficiency, nChambers, debug);
      hEfficiency->SetMarkerColor(simSpectrumE->GetLineColor());
      hEfficiency->SetLineColor(simSpectrumE->GetLineColor());
      fitEfficiency->SetMarkerColor(simSpectrumE->GetLineColor());
      fitEfficiency->SetLineColor(simSpectrumE->GetLineColor());
      fitEfficiency->SetFillColor(simSpectrumE->GetLineColor());
      fitEfficiency->SetFillStyle(3003);
      output->cd();
      hEfficiency->Write("Efficiency_simulated", TObject::kOverwrite);
      fitEfficiency->Write("Efficiency_simulated_fit", TObject::kOverwrite);
      leg3->AddEntry(simSpectrum,"simulated","l");
      cEfficiency->cd()->SetLogy(1);
      pEfficiency->DrawCopy();
      g1->Add(hEfficiency,"P");
      g1->Add(fitEfficiency,"a3");
      //hEfficiency->DrawClone("P,text0");
      CalcPionEfficiency(preCompBaselineSpectrumElectron, preCompBaselineSpectrum, hEfficiency, fitEfficiency, nChambers, debug);
      hEfficiency->SetMarkerColor(preCompBaselineSpectrumElectron->GetLineColor());
      hEfficiency->SetLineColor(preCompBaselineSpectrumElectron->GetLineColor());
      fitEfficiency->SetMarkerColor(preCompBaselineSpectrumElectron->GetLineColor());
      fitEfficiency->SetLineColor(preCompBaselineSpectrumElectron->GetLineColor());
      fitEfficiency->SetFillColor(preCompBaselineSpectrumElectron->GetLineColor());
      fitEfficiency->SetFillStyle(3003);
      output->cd();
      hEfficiency->Write("Efficiency_measured_MS", TObject::kOverwrite);
      fitEfficiency->Write("Efficiency_measured_MS_fit", TObject::kOverwrite);
      //if (!simple) leg3->AddEntry(preCompBaselineSpectrumElectron,"measured: MS noise + baseline","l");
      cEfficiency->cd();
      //if (!simple) hEfficiency->DrawClone("P,text0,same");
      CalcPionEfficiency(baselineNoiseSpectrumElectron, baselineNoiseSpectrum, hEfficiency, fitEfficiency, nChambers, debug);
      hEfficiency->SetMarkerColor(baselineNoiseSpectrumElectron->GetLineColor());
      hEfficiency->SetLineColor(baselineNoiseSpectrumElectron->GetLineColor());
      fitEfficiency->SetMarkerColor(baselineNoiseSpectrumElectron->GetLineColor());
      fitEfficiency->SetLineColor(baselineNoiseSpectrumElectron->GetLineColor());
      fitEfficiency->SetFillColor(baselineNoiseSpectrumElectron->GetLineColor());
      fitEfficiency->SetFillStyle(3003);
      output->cd();
      hEfficiency->Write("Efficiency_measured_FFM", TObject::kOverwrite);
      fitEfficiency->Write("Efficiency_measured_FFM_fit", TObject::kOverwrite);
      leg3->AddEntry(baselineNoiseSpectrumElectron,"measured: 1.baseline 2.noise","l");
      cEfficiency->cd();
      g1->Add(hEfficiency,"P");
      g1->Add(fitEfficiency,"a3");
      //hEfficiency->DrawClone("P,text0,same");
      CalcPionEfficiency(baselineNoiseSpectrumWoOWoUElectron, baselineNoiseSpectrumWoOWoU, hEfficiency, fitEfficiency, nChambers, debug);
      hEfficiency->SetMarkerColor(baselineNoiseSpectrumWoOWoUElectron->GetLineColor());
      hEfficiency->SetLineColor(baselineNoiseSpectrumWoOWoUElectron->GetLineColor());
      fitEfficiency->SetMarkerColor(baselineNoiseSpectrumWoOWoUElectron->GetLineColor());
      fitEfficiency->SetLineColor(baselineNoiseSpectrumWoOWoUElectron->GetLineColor());
      fitEfficiency->SetFillColor(baselineNoiseSpectrumWoOWoUElectron->GetLineColor());
      fitEfficiency->SetFillStyle(3003);
      output->cd();
      hEfficiency->Write("Efficiency_measured_FFM_WoOWo", TObject::kOverwrite);
      fitEfficiency->Write("Efficiency_measured_FFM_WoOWo_fit", TObject::kOverwrite);
      leg3->AddEntry(baselineNoiseSpectrumWoOWoUElectron,"measured: 1.baseline 2.noise (woO && woU)","l");
      cEfficiency->cd();
      g1->Add(hEfficiency,"P");
      g1->Add(fitEfficiency,"a3");
      g1->Draw("same");
      //hEfficiency->DrawClone("P,text0,same");
      leg3->Draw("same");
      TLine *l1 = new TLine(0.5,1,12.5,1);
      l1->SetLineColor(2);
      l1->SetLineStyle(2);
      l1->Draw("same");
      outpic.ReplaceAll("Spectra","Efficiency");
      cEfficiency->SaveAs(outpic+".pdf");
      cEfficiency->SaveAs(outpic+".png");
    }
    if (!simple) {
      hEfficiency->GetXaxis()->SetRangeUser(1e-2,1e2);
      TString names;
      for(Int_t delta = 0; delta < 6; delta++) {
        hEfficiency->GetYaxis()->SetRangeUser(1e-2,1e2);
        hEfficiency->SetMarkerColor(baselineNoiseSpectrumCrossTalkP[delta]->GetLineColor());
        hEfficiency->SetLineColor(baselineNoiseSpectrumCrossTalkP[delta]->GetLineColor());
        CalcPionEfficiency(baselineNoiseSpectrumCrossTalkE[delta], baselineNoiseSpectrumCrossTalkP[delta], hEfficiency, fitEfficiency, nChambers, debug);
        names.Form("Efficiency_Crosstalk_delta%02i",delta);
        hEfficiency->Write(names, TObject::kOverwrite);
        fitEfficiency->Write(names + TString("_fit"), TObject::kOverwrite);
      }

      CalcPionEfficiency(simSpectrumE, simSpectrum, hEfficiency, fitEfficiency, nChambers, debug);
      hEfficiency->SetMarkerColor(simSpectrumE->GetLineColor());
      hEfficiency->SetLineColor(simSpectrumE->GetLineColor());
      fitEfficiency->SetMarkerColor(simSpectrumE->GetLineColor());
      fitEfficiency->SetLineColor(simSpectrumE->GetLineColor());
      fitEfficiency->SetFillColor(simSpectrumE->GetLineColor());
      fitEfficiency->SetFillStyle(3003);
      output->cd();
      hEfficiency->Write("Efficiency_simulated", TObject::kOverwrite);
      fitEfficiency->Write("Efficiency_simulated_fit", TObject::kOverwrite);
      leg3->AddEntry(simSpectrum,"simulated","l");
      cEfficiency->cd()->SetLogy(1);
      //hEfficiency->DrawClone("P,text0");
      pEfficiency->DrawCopy("P");
      g1->Add(hEfficiency,"P");
      g1->Add(fitEfficiency,"P");
      CalcPionEfficiency(preCompBaselineSpectrumElectron, preCompBaselineSpectrum, hEfficiency, fitEfficiency, nChambers, debug);
      hEfficiency->SetMarkerColor(preCompBaselineSpectrumElectron->GetLineColor());
      hEfficiency->SetLineColor(preCompBaselineSpectrumElectron->GetLineColor());
      fitEfficiency->SetMarkerColor(preCompBaselineSpectrumElectron->GetLineColor());
      fitEfficiency->SetLineColor(preCompBaselineSpectrumElectron->GetLineColor());
      fitEfficiency->SetFillColor(preCompBaselineSpectrumElectron->GetLineColor());
      fitEfficiency->SetFillStyle(3003);
      output->cd();
      hEfficiency->Write("Efficiency_measured_MS", TObject::kOverwrite);
      fitEfficiency->Write("Efficiency_measured_MS_fit", TObject::kOverwrite);
      leg3->AddEntry(preCompBaselineSpectrumElectron,"measured: MS 1.baseline 2.noise","l");
      cEfficiency->cd();
      g1->Add(hEfficiency,"P");
      g1->Add(fitEfficiency,"P");
      //hEfficiency->DrawClone("P,text0,same");
      CalcPionEfficiency(baselineNoiseSpectrumElectron, baselineNoiseSpectrum, hEfficiency, fitEfficiency, nChambers, debug);
      hEfficiency->SetMarkerColor(baselineNoiseSpectrumElectron->GetLineColor());
      hEfficiency->SetLineColor(baselineNoiseSpectrumElectron->GetLineColor());
      fitEfficiency->SetMarkerColor(baselineNoiseSpectrumElectron->GetLineColor());
      fitEfficiency->SetLineColor(baselineNoiseSpectrumElectron->GetLineColor());
      fitEfficiency->SetFillColor(baselineNoiseSpectrumElectron->GetLineColor());
      fitEfficiency->SetFillStyle(3003);
      output->cd();
      hEfficiency->Write("Efficiency_measured_FFM", TObject::kOverwrite);
      fitEfficiency->Write("Efficiency_measured_FFM_fit", TObject::kOverwrite);
      leg3->AddEntry(baselineNoiseSpectrumElectron,"measured: FFM 1.baseline 2.noise","l");
      cEfficiency->cd();
      //hEfficiency->DrawClone("P,text0,same");
      g1->Add(hEfficiency,"P");
      g1->Add(fitEfficiency,"P");
      CalcPionEfficiency(baselineNoiseSpectrumWoOWoUElectron, baselineNoiseSpectrumWoOWoU, hEfficiency, fitEfficiency, nChambers, debug);
      hEfficiency->SetMarkerColor(baselineNoiseSpectrumWoOWoUElectron->GetLineColor());
      hEfficiency->SetLineColor(baselineNoiseSpectrumWoOWoUElectron->GetLineColor());
      fitEfficiency->SetMarkerColor(baselineNoiseSpectrumWoOWoUElectron->GetLineColor());
      fitEfficiency->SetLineColor(baselineNoiseSpectrumWoOWoUElectron->GetLineColor());
      fitEfficiency->SetFillColor(baselineNoiseSpectrumWoOWoUElectron->GetLineColor());
      fitEfficiency->SetFillStyle(3003);
      output->cd();
      hEfficiency->Write("Efficiency_measured_FFM_WoOWo", TObject::kOverwrite);
      fitEfficiency->Write("Efficiency_measured_FFM_WoOWo_fit", TObject::kOverwrite);
      leg3->AddEntry(baselineNoiseSpectrumWoOWoUElectron,"measured: FFM 1.baseline 2.noise (woO && woU)","l");
      cEfficiency->cd();
      g1->Add(hEfficiency,"P");
      g1->Add(fitEfficiency,"P");
      //hEfficiency->DrawClone("P,text0,same");
      g1->Draw("same");
      leg3->Draw("same");
      TLine *l1 = new TLine(0.5,1,12.5,1);
      l1->SetLineColor(2);
      l1->SetLineStyle(2);
      l1->Draw("same");
      outpic.ReplaceAll("Spectra","Efficiency");
      cEfficiency->SaveAs(outpic+".pdf");
      cEfficiency->SaveAs(outpic+".png");
    }
  }
  TLegend *lShape = new TLegend(0.5,0.75,0.85,0.85);
  lShape->SetFillStyle(0);
  lShape->SetTextSize(0.025);
  lShape->SetLineStyle(0);
  lShape->SetLineColor(0);
  TCanvas *r = new TCanvas("r","r",2*800,2*600); 
  r->Divide(2,2);
  r->cd(1); 
  averagePulse_e->GetYaxis()->SetRangeUser(-1,256);
  averagePulse_e->DrawCopy();
  lShape->AddEntry(averagePulse_e,"electrons","lp");
  averagePulse_p->DrawCopy("same");
  lShape->AddEntry(averagePulse_p,"pions","lp");
  lShape->Draw("same");  
  r->cd(3)->SetLogz(1); 
  averagePulse_2D_e->DrawCopy("colz");
  r->cd(4)->SetLogz(1);
  averagePulse_2D_p->DrawCopy("colz");
  r->cd(2)->SetLogz(1);
  averageSignal_2D->DrawCopy("colz");
  r->SaveAs(outpic+"PulseShape.pdf");
  r->SaveAs(outpic+"PulseShape.png");
  TCanvas *s = new TCanvas("s","s",800,600); 
  s->cd(1)->SetLogz(1);   
  keV2Amplitude->SetXTitle("dE/dx [keV]");
  keV2Amplitude->SetYTitle("Amplitude [ADCs]");
  keV2Amplitude->DrawCopy("colz");
  keV2AmplitudeProfile->SetMarkerStyle(24);
  keV2AmplitudeProfile->DrawCopy("same");
  TF1 *linear = new TF1("linear","pol1",2,30);
  keV2Amplitude->Fit("linear","RQ0");
  linear->SetLineColor(2);
  linear->SetLineStyle(1);
  linear->Draw("same");
  keV2AmplitudeProfile->DrawCopy("same");
  TString fitResult;
  fitResult.Form("f(x) = (%.3f#pm%.3f) x + (%.3f#pm%.3f)",linear->GetParameter(1), linear->GetParError(1), linear->GetParameter(0), linear->GetParError(0));
  TPaveText *pt = new TPaveText(30,20,95,50);
  //pt->SetFillStyle(0);
  pt->SetLineColor(0);
  pt->SetTextSize(0.035);
  pt->AddText(fitResult.Data());
  pt->Draw("same");
  s->SaveAs(outpic+"_Pions_ADC_vs_Amplitude.pdf");
  s->SaveAs(outpic+"_Pions_ADC_vs_Amplitude.png");
  IntegralvsAmplitude->Draw("colz");
  s->SaveAs(outpic+"_Pions_Integral_vs_Amplitude.pdf");
  s->SaveAs(outpic+"_Pions_Integral_vs_Amplitude.png");
  TCanvas* cBeamMonitor = new TCanvas("BeamMonitor","Beam Monitor",2*800,600);
  cBeamMonitor->Divide(2,1);
  cBeamMonitor->cd(1)->SetLogz(1);
  allCh1_Pb->GetYaxis()->SetTitleOffset(1.35);
  allCh1_Pb->DrawCopy("colz");
  cBeamMonitor->cd(2)->SetLogz(1);
  allCh2_Pb->GetYaxis()->SetTitleOffset(1.35);
  allCh2_Pb->DrawCopy("colz");
  cBeamMonitor->SaveAs(outpic+"Ch_Pb.pdf");
  cBeamMonitor->SaveAs(outpic+"Ch_Pb.png");
  TCanvas* cCutMonitor = new TCanvas("CutMonitor","Cut Value Monitor",3*800,3*600);
  cCutMonitor->Divide(3,3);
  //cCutMonitor->cd(1)->SetLogx(1);
  cCutMonitor->cd(1)->SetLogy(1);
  maxAmplitudeValue->DrawCopy();
  cCutMonitor->cd(2)->SetLogy(1);   
  covaMatixValue->DrawCopy();
  cCutMonitor->cd(3)->SetLogy(0); 
  clusterSize->DrawCopy();
  cCutMonitor->cd(4)->SetLogz(1); 
  covaMatixValueMaxAmplitude->DrawCopy("colz");
  cCutMonitor->cd(5)->SetLogz(1);
  covaMatixValueHitTime->DrawCopy("colz");
  cCutMonitor->cd(6)->SetLogz(1);
  maxAmplitudeHitTime->DrawCopy("colz");
  cCutMonitor->cd(7)->SetLogz(1);
  clusterwCharge->DrawCopy("colz");
  cCutMonitor->SaveAs(outpic+"Cuts.pdf");
  cCutMonitor->SaveAs(outpic+"Cuts.png");
  TCanvas *cSimple = new TCanvas("cSimple","cSimple",800,600);
  cSimple->cd()->SetLogy(1);
  if (fHistoScale == SPADIC_TRACE_SIZE * SPADIC_ADC_MAX / 100.)
    baselineNoiseSpectrum->SetXTitle("ADC values");
  else
    baselineNoiseSpectrum->SetXTitle("dE/dx [keV]");
  baselineNoiseSpectrum->SetYTitle("normalized counts");
  baselineNoiseSpectrum->SetStats(false);
  baselineNoiseSpectrum->SetLineStyle(1);
  baselineNoiseSpectrum->SetLineColor(1);
  //simSpectrumE->Scale(/*baselineNoiseSpectrum->GetEntries()*/1./(Float_t)simSpectrumE->GetEntries());
  if (fHistoScale < SPADIC_TRACE_SIZE * SPADIC_ADC_MAX / 100.)
    baselineNoiseSpectrum->GetXaxis()->SetRangeUser(-0.5,70);
  baselineNoiseSpectrum->DrawCopy();
  baselineNoiseSpectrumElectron->SetLineColor(2);
  baselineNoiseSpectrumElectron->DrawCopy("same");
  outpic.ReplaceAll("Efficiency","Spectra");
  cSimple->SaveAs(outpic+"_Simple.pdf");
  cSimple->SaveAs(outpic+"_Simple.png");
  if (!plotFromFile) {
    printf("%i Pion (%.2f%%) & %i Electron (%.2f%%) & %i no PID or noise (%.2f%%) events found\n",nPions,nPions*100./entries,nElectrons,nElectrons*100./entries,entries-(nPions+nElectrons),(entries-(nPions+nElectrons))*100./entries);
    hCh_1->Write("", TObject::kOverwrite);
    hCh_2->Write("", TObject::kOverwrite);
    hPb->Write("", TObject::kOverwrite);
    Ch1_Pb->Write("", TObject::kOverwrite);
    allCh1_Pb->Write("", TObject::kOverwrite);
    noCh1_Pb->Write("", TObject::kOverwrite);
    Ch2_Pb->Write("", TObject::kOverwrite);
    allCh2_Pb->Write("", TObject::kOverwrite);
    noCh2_Pb->Write("", TObject::kOverwrite);
    Ch1_Ch2->Write("", TObject::kOverwrite);
    allCh1_Ch2->Write("", TObject::kOverwrite);
    noCh1_Ch2->Write("", TObject::kOverwrite);
    prfProfile->Write("", TObject::kOverwrite);
    prf->Write("", TObject::kOverwrite);
    clusterwCharge->Write("", TObject::kOverwrite);
    mathieson->Write("", TObject::kOverwrite);
    mathiesonFit->Write("", TObject::kOverwrite);
    cPRF->Write("", TObject::kOverwrite);
    cEfficiency->Write("", TObject::kOverwrite);
    simSpectrum->Write("", TObject::kOverwrite);
    simSpectrumE->Write("", TObject::kOverwrite);
    rawSpectrum->Write("", TObject::kOverwrite);
    baselineSpectrum->Write("", TObject::kOverwrite);
    baselineNoiseSpectrum->Write("", TObject::kOverwrite);
    noiseSpectrum->Write("", TObject::kOverwrite);
    noiseBaselineSpectrum->Write("", TObject::kOverwrite);
    baselineNoiseSpectrumElectron->Write("", TObject::kOverwrite);
    baselineNoiseSpectrumWoOWoUElectron->Write("", TObject::kOverwrite);
    baselineNoiseSpectrumWoOWoU->Write("", TObject::kOverwrite);
    preCompSpectrumElectron->Write("", TObject::kOverwrite);
    preCompBaselineSpectrumElectron->Write("", TObject::kOverwrite);
    preCompSpectrum->Write("", TObject::kOverwrite);
    preCompBaselineSpectrum->Write("", TObject::kOverwrite);
    averagePulse_e->Write("", TObject::kOverwrite);
    averagePulse_p->Write("", TObject::kOverwrite);
    averagePulse_2D_e->Write("", TObject::kOverwrite);
    averagePulse_2D_p->Write("", TObject::kOverwrite);
    probabilityOverflow->Write("", TObject::kOverwrite);
    probabilityUnderflow->Write("", TObject::kOverwrite);
    probabilityOverflowE->Write("", TObject::kOverwrite);
    probabilityUnderflowE->Write("", TObject::kOverwrite);
    keV2AmplitudeProfile->Write("", TObject::kOverwrite);
    keV2Amplitude->Write("", TObject::kOverwrite);
    IntegralvsAmplitude->Write("", TObject::kOverwrite);
    for(Int_t delta = 0; delta < 6; delta++) {
      baselineNoiseSpectrumCrossTalkE[delta]->Scale(1./(Double_t)baselineNoiseSpectrumCrossTalkE[delta]->Integral());
      baselineNoiseSpectrumCrossTalkE[delta]->Write("", TObject::kOverwrite);
      baselineNoiseSpectrumCrossTalkP[delta]->Scale(1./(Double_t)baselineNoiseSpectrumCrossTalkP[delta]->Integral());
      baselineNoiseSpectrumCrossTalkP[delta]->Write("", TObject::kOverwrite);  
    }
    //deltaSim->Write("", TObject::kOverwrite);
    //deltaSim2->Write("", TObject::kOverwrite);
    maxAmplitudeValue->Write("", TObject::kOverwrite);
    covaMatixValue->Write("", TObject::kOverwrite);
    covaMatixValueMaxAmplitude->Write("", TObject::kOverwrite);
    covaMatixValueHitTime->Write("", TObject::kOverwrite);
    clusterSize->Write("", TObject::kOverwrite);
    maxAmplitudeHitTime->Write("", TObject::kOverwrite);
    crossTalk->Write("", TObject::kOverwrite);
    cResults->Write("", TObject::kOverwrite);
    noiseDistribution->SetXTitle("noise [ADC]");
    noiseDistribution->SetYTitle("#");
    noiseDistribution->Write("", TObject::kOverwrite);
    baselineDistribution->SetXTitle("channel");
    baselineDistribution->SetYTitle("ADC");
    baselineDistribution->Write("", TObject::kOverwrite);
  }
  output->Close();
  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();

  printf("\n\n******************** Reading Test  **********************\n");
  printf("   RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
  printf("*********************************************************\n\n");

}