beamtime/cern-oct11/go4/AnalysisMacros/spadic_noise_Analysis.cxx (r4864/r4851)

Go to the documentation of this file.

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


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

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

#include "TGraph.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 "TCanvas.h"

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

#include "TRocEvent.h"
#include "roc/Message.h"
#include "TSpadicEvent.h"


#include "TMbsCrateEvent.h"
#include "TBeamMonitorEvent.h"
#include "TCernOct11UnpackEvent.h"

//#include "CBMsimTR.C"

#define NUM_SPADIC_CHA             8
#define SPADIC_TRACE_SIZE         45
#define SPADIC_ADC_MAX           255
#define noBaselineTB               5 
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;
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 *signalChDistance = new TH1I("signalChDistance","signalChDistance",2*NUM_SPADIC_CHA+1,-NUM_SPADIC_CHA-0.5,NUM_SPADIC_CHA+0.5);

Int_t pidcuts[18] = {0};

void Statusbar(Int_t i, Int_t n) {
  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;
  }
}

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 = 100000;
  TH1D* hElectronProbability = (TH1D*)hElectronSpectra->Clone("ElectronProbability");
  //hElectronProbability->Scale(1./hElectronProbability->GetEntries());
  TH1D* hPionProbability = (TH1D*)hPionSpectra->Clone("hPionProbability");
  //hPionProbability->Scale(1./hPionProbability->GetEntries());
  TCanvas *cTemp2;
  TCanvas *cTemp;
  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 [%]");

  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("#");

    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();        
          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();  
            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();  
            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);
          
    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();
    }
  }

  for (Int_t Chamberloop = 1; Chamberloop <= nChambers; Chamberloop++){
    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];
  pidcuts[0] =1; // // El Ch1 min
  pidcuts[1] =1; // // El Ch1 max
  pidcuts[2] =1; // // El Ch2 min
  pidcuts[3] =1; // // El Ch2 max
  pidcuts[4] =1; // // El Pb min
  pidcuts[5] =1; // // El Pb max
  pidcuts[6] =1; // // Pi Ch1 min
  pidcuts[7] =1; // // Pi Ch1 max
  pidcuts[8] =1; // // Pi Ch2 min
  pidcuts[9] =1; // // Pi Ch2 max
  pidcuts[10]=1; // // Pi Pb min
  pidcuts[11]=1; // // Pi Pb max
  pidcuts[12]=1; // // Mu Ch1 min
  pidcuts[13]=1; // // Mu Ch1 max
  pidcuts[14]=1; // // Mu Ch2 min
  pidcuts[15]=1; // // Mu Ch2 max
  pidcuts[16]=1; // // Mu Pb min
  pidcuts[17]=1; // // Mu Pb max
  if (runnumber < 1910000 )       cout << "[WARN] No PID cuts set! Invalid runnumber " << runnumber << endl;
  else if (runnumber >= 2010010 && runnumber <= 2410021) {//3GeV
    pidcuts[0] =500; pidcuts[1] =9999; pidcuts[2] =500; pidcuts[3] =9999; pidcuts[4] =1000; pidcuts[5] =9999; 
    pidcuts[6] =  0; pidcuts[7] = 400; pidcuts[8] =  0; pidcuts[9] = 400; pidcuts[10]=   0; pidcuts[11]= 900; 
    pidcuts[12]=  1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
  }
  else if (runnumber > 2410021 && runnumber <= 2510014) {//3GeV
    pidcuts[0] =500; pidcuts[1] =9999; pidcuts[2] =500; pidcuts[3] =9999; pidcuts[4] =2000; pidcuts[5] =9999; 
    pidcuts[6] =  0; pidcuts[7] = 400; pidcuts[8] =  0; pidcuts[9] = 400; pidcuts[10]=   0; pidcuts[11]=1500; 
    pidcuts[12]=  1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
  }
  else if (runnumber > 2510014 && runnumber <= 2610002) {//3GeV
    pidcuts[0] =500; pidcuts[1] =9999; pidcuts[2] =500; pidcuts[3] =9999; pidcuts[4] =2000; pidcuts[5] =9999; 
    pidcuts[6] =  0; pidcuts[7] = 400; pidcuts[8] =  0; pidcuts[9] = 400; pidcuts[10]=   0; pidcuts[11]=1500; 
    pidcuts[12]=  1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
  }
  else if (runnumber > 2610002 && runnumber <= 2610005) {//3GeV
    pidcuts[0] =500; pidcuts[1] =9999; pidcuts[2] =500; pidcuts[3] =9999; pidcuts[4] =2500; pidcuts[5] =9999; 
    pidcuts[6] =  0; pidcuts[7] = 400; pidcuts[8] =  0; pidcuts[9] = 400; pidcuts[10]=   0; pidcuts[11]=1500; 
    pidcuts[12]=  1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
  }
  else if (runnumber > 2610005 && runnumber <= 3010011) {//3GeV
    pidcuts[0] =500; pidcuts[1] =9999; pidcuts[2] =500; pidcuts[3] =9999; pidcuts[4] =2000; pidcuts[5] =9999; 
    pidcuts[6] =  0; pidcuts[7] = 400; pidcuts[8] =  0; pidcuts[9] = 400; pidcuts[10]=   0; pidcuts[11]=1500; 
    pidcuts[12]=  1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
  }
 
  else if (runnumber > 3010015 )  cout << "[WARN] No PID cuts set! Invalid runnumber " << runnumber << endl;
  else                           {cout << "[FATAL] Runnumber totaly f**ked up... exiting root NOW! " << endl;}
  /*
    if (runnumber < 1910000 )       cout << "[WARN] No PID cuts set! Invalid runnumber " << runnumber << endl;
  */
  //  else if (runnumber == 2210011) {pidcuts[0]=600; pidcuts[1]=1050; pidcuts[2]=550; pidcuts[3]=850; pidcuts[4]=1100; pidcuts[5]=1500; pidcuts[6]=305; pidcuts[7]=318; pidcuts[8]=335; pidcuts[9]=350; pidcuts[10]=400; pidcuts[11]=1000; pidcuts[12]=1; pidcuts[13]=1; pidcuts[14]=1; pidcuts[15]=1; pidcuts[16]=1; pidcuts[17]=1;}
  //else if (runnumber >= 2010000 && runnumber <= 2010008) {
  //pidcuts[0] =600; pidcuts[1] =9999; pidcuts[2] =550; pidcuts[3] =9999; pidcuts[4] =1100; pidcuts[5] =9999; 
  //pidcuts[6] =  0; pidcuts[7] = 400; pidcuts[8] =  0; pidcuts[9] = 400; pidcuts[10]=   0; pidcuts[11]=1000; 
  //pidcuts[12]=  1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
  /*
    }
    else if (runnumber > 2110000 && runnumber <= 2110015) {
    pidcuts[0] =600; pidcuts[1] =9999; pidcuts[2] =550; pidcuts[3] =9999; pidcuts[4] =1100; pidcuts[5] =9999; 
    pidcuts[6] =305; pidcuts[7] = 318; pidcuts[8] =335; pidcuts[9] = 350; pidcuts[10]= 400; pidcuts[11]=9999; 
    pidcuts[12]=  1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
    }
    else if (runnumber > 2110015 && runnumber <= 2110020) {
    pidcuts[0] =500; pidcuts[1] =9999; pidcuts[2] =500; pidcuts[3] =9999; pidcuts[4] =1000; pidcuts[5] =9999; 
    // pidcuts[6] =290; pidcuts[7] = 320; pidcuts[8] =330; pidcuts[9] = 350; pidcuts[10]=   1; pidcuts[11]=9999;
    pidcuts[6] =  0; pidcuts[7] = 400; pidcuts[8] =  0; pidcuts[9] = 400; pidcuts[10]=   1; pidcuts[11]=9999;
    pidcuts[12]=  1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
    }
    else if (runnumber > 2110020 && runnumber <= 2210000) {
    pidcuts[0] =500; pidcuts[1] =9999; pidcuts[2] =520; pidcuts[3] =9999; pidcuts[4] =1000; pidcuts[5] =9999; 
    pidcuts[6] =303; pidcuts[7] = 323; pidcuts[8] =328; pidcuts[9] = 346; pidcuts[10]=   1; pidcuts[11]=9999; 
    pidcuts[12]=  1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
    }
    else if (runnumber > 2210000 && runnumber <= 2210009) {
    pidcuts[0] =500; pidcuts[1] =9999; pidcuts[2] =520; pidcuts[3] =9999; pidcuts[4] =1000; pidcuts[5] =9999; 
    pidcuts[6] =305; pidcuts[7] = 325; pidcuts[8] =325; pidcuts[9] = 349; pidcuts[10]=   1; pidcuts[11]=9999; 
    pidcuts[12]=  1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
    }
    else if (runnumber > 2210009 && runnumber <= 2210011) {
    pidcuts[0] =500; pidcuts[1] =9999; pidcuts[2] =520; pidcuts[3] =9999; pidcuts[4] =1000; pidcuts[5] =9999; 
    pidcuts[6] =303; pidcuts[7] = 323; pidcuts[8] =325; pidcuts[9] = 348; pidcuts[10]=   1; pidcuts[11]=9999; 
    pidcuts[12]=  1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
    }
    else if (runnumber > 2210011 && runnumber <= 2210012) {
    pidcuts[0] =500; pidcuts[1] =9999; pidcuts[2] =520; pidcuts[3] =9999; pidcuts[4] =1000; pidcuts[5] =9999; 
    pidcuts[6] =303; pidcuts[7] = 323; pidcuts[8] =325; pidcuts[9] = 348; pidcuts[10]=   1; pidcuts[11]=9999; 
    pidcuts[12]=  1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
    } 
    else if (runnumber > 2210012 && runnumber <= 2310001) {
    pidcuts[0] =500; pidcuts[1] =9999; pidcuts[2] =520; pidcuts[3] =9999; pidcuts[4] =1000; pidcuts[5] =9999; 
    pidcuts[6] =300; pidcuts[7] = 325; pidcuts[8] =325; pidcuts[9] = 348; pidcuts[10]=   1; pidcuts[11]=9999; 
    pidcuts[12]=  1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
    }
    else if (runnumber > 2310001 && runnumber <= 2310001) {
    pidcuts[0] =600; pidcuts[1] =9999; pidcuts[2] =550; pidcuts[3] =9999; pidcuts[4] =1100; pidcuts[5] =9999; 
    pidcuts[6] =305; pidcuts[7] = 318; pidcuts[8] =335; pidcuts[9] = 350; pidcuts[10]= 400; pidcuts[11]=1000; 
    pidcuts[12]=  1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
    }
    else if (runnumber > 2310001 && runnumber <= 2310004) {
    pidcuts[0] =500; pidcuts[1] =9999; pidcuts[2] =520; pidcuts[3] =9999; pidcuts[4] =1000; pidcuts[5] =9999; 
    pidcuts[6] =310; pidcuts[7] = 330; pidcuts[8] =325; pidcuts[9] = 355; pidcuts[10]=   1; pidcuts[11]=9999; 
    pidcuts[12]=  1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
    }  
    else if (runnumber > 2310004 && runnumber <= 2310008) {
    pidcuts[0] =500; pidcuts[1] =9999; pidcuts[2] =520; pidcuts[3] =9999; pidcuts[4] =1000; pidcuts[5] =9999; 
    pidcuts[6] =310; pidcuts[7] = 330; pidcuts[8] =325; pidcuts[9] = 355; pidcuts[10]=   1; pidcuts[11]=9999; 
    pidcuts[12]=  1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
    }
    else if (runnumber > 2310008 && runnumber <= 2710001) {
    pidcuts[0] =1000; pidcuts[1] =9999; pidcuts[2] =550; pidcuts[3] =9999; pidcuts[4] =2000; pidcuts[5] =9999; 
    pidcuts[6] = 308; pidcuts[7] = 328; pidcuts[8] =350; pidcuts[9] = 374; pidcuts[10]=   1; pidcuts[11]=9999; 
    pidcuts[12]=   1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
    }
    else if (runnumber > 2710001 && runnumber <= 2910006) {
    pidcuts[0] =1000; pidcuts[1] =9999; pidcuts[2] =600; pidcuts[3] =9999; pidcuts[4] =2000; pidcuts[5] =9999; 
    pidcuts[6] = 315; pidcuts[7] = 330; pidcuts[8] =355; pidcuts[9] = 380; pidcuts[10]=   1; pidcuts[11]=9999; 
    pidcuts[12]=   1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
    }
    else if (runnumber > 2710006 && runnumber <= 2910010) {
    pidcuts[0] =1000; pidcuts[1] =9999; pidcuts[2] =600; pidcuts[3] =9999; pidcuts[4] =2000; pidcuts[5] =9999; 
    pidcuts[6] = 300; pidcuts[7] = 320; pidcuts[8] =330; pidcuts[9] = 360; pidcuts[10]=   1; pidcuts[11]=9999; 
    pidcuts[12]=   1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
    }

    else if (runnumber > 2910010 && runnumber <= 3010015) {
    pidcuts[0] =1000; pidcuts[1] =2000; pidcuts[2] =650; pidcuts[3] =1500; pidcuts[4] =2000; pidcuts[5] =9999; 
    pidcuts[6] = 290; pidcuts[7] = 320; pidcuts[8] =325; pidcuts[9] = 360; pidcuts[10]= 400; pidcuts[11]= 800; 
    pidcuts[12]=   1; pidcuts[13]=9999; pidcuts[14]=  1; pidcuts[15]=9999; pidcuts[16]= 850; pidcuts[17]=1100;
    }
    else if (runnumber > 3010015 )  cout << "[WARN] No PID cuts set! Invalid runnumber " << runnumber << endl;
    else                           {cout << "[FATAL] Runnumber totaly f**ked up... exiting root NOW! " << endl;}
  */
 
}
void getPID(/*TMbsCrateEvent* input, TString runId*/ Float_t ch1, Float_t ch2, Float_t pb)
{
  //int *pidcuts = GetPidCuts(runId);
  
  //cout << "test: " << pidcuts[7] << endl;

  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] && // MUON
           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;
  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;
  }
*/
void CalcPrf(TH1D* inPulse[NUM_SPADIC_CHA], TH2D* prf, TProfile* prfProfile, Double_t padWidth)
{
  const Int_t minBin = 0;
  const Int_t maxBin = SPADIC_TRACE_SIZE-1;
  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));
      }
    }
    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 chargePercent = intensCh[maxCh-1] / (intensCh[maxCh-1] + intensCh[maxCh] + intensCh[maxCh+1]);
    prf->Fill((-dxPos-1)*padWidth, chargePercent);
    prfProfile->Fill((-dxPos-1)*padWidth ,chargePercent);
    chargePercent = intensCh[maxCh] / (intensCh[maxCh-1] + intensCh[maxCh] + intensCh[maxCh+1]);
    prf->Fill((dxPos)*padWidth,chargePercent);
    prfProfile->Fill((dxPos)*padWidth ,chargePercent);
    chargePercent = intensCh[maxCh+1] / (intensCh[maxCh-1] + intensCh[maxCh] + intensCh[maxCh+1]);
    prf->Fill((-dxPos+1)*padWidth,chargePercent);
    prfProfile->Fill((-dxPos+1)*padWidth ,chargePercent);
  }
}
void FillAverageSignal(TH1D* inPulse[NUM_SPADIC_CHA], TProfile* averagePulse)
{
  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)));
  }
}

  /*
    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;
    }
  */




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

    return false;
  }

  Bool_t BorderPadMax(TH1D* inPulse[NUM_SPADIC_CHA]) {
    Int_t 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) > 50)
      return true;
    return false;
  }

  Bool_t minimumAmplitudeThreshold(TH1D* inPulse[NUM_SPADIC_CHA]) {
    for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++)
      if (inPulse[ch]->GetBinContent(inPulse[ch]->GetMaximumBin()) > 20/*17/*15*/)
        return true;
    return false;
  }


  void BaselineCompensation(TH1D* inPulse[NUM_SPADIC_CHA], TH1D* outPulse[NUM_SPADIC_CHA])
  {
    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] /= noBaselineTB;
    }

    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;
    }
    }
    }
  */
  void CancelNoise_Cova(TH1D* inPulse[NUM_SPADIC_CHA],TH1D* outPulse[NUM_SPADIC_CHA], bool debug=false)
  {
    const double lowercorthreshold = 0.112;
    //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 minch = -1;//GetMinimumChannel(input);
    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;
      }
    }

    //Int_t maxCh = GetPadMax(inPulse); //debug

    int noise_ch_counter=0;

    Double_t thisisnoise[SPADIC_TRACE_SIZE] = {0.0};

    //double val_in=0;
    //double val_noise=0;


    for(int t=0;t<SPADIC_TRACE_SIZE;t++)
      {
        principal->AddRow(inputarray[t]);
      }
  
    principal->MakePrincipals();
    //TMatrixD* cova = new TMatrixD(NUM_SPADIC_CHA,NUM_SPADIC_CHA);
    TMatrixD* cova = (TMatrixD*)principal->GetCovarianceMatrix();
  
    Int_t lastSigCh = -1;

    for(int ch=0;ch<NUM_SPADIC_CHA;ch++)
      {
        cortest[ch] = TMatrixDRow (*cova,minch)(ch);
        if(cortest[ch]>lowercorthreshold)
          {
            noise_ch_counter++;
            //thisisnoise->Add(input[ch]);
            for(Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++){
              thisisnoise[bin] += inputarray[bin][ch];
            }
          }
        else { //debug
          if (lastSigCh != -1)
            signalChDistance->Fill(ch - lastSigCh/*maxCh-ch*/); //debug
          lastSigCh = ch;
        } //debug
      }
    //thisisnoise->Scale(1./(double)noise_ch_counter);
    for(Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++){
      thisisnoise[bin] /= (double)noise_ch_counter;
    }
  
    clusterSize->Fill(NUM_SPADIC_CHA-noise_ch_counter); //debug



    for (int ch=0; ch<NUM_SPADIC_CHA; ch++)
      {
        for (int k=0; k<SPADIC_TRACE_SIZE; k++)
          {
            /*
              val_in=input[ch]->GetBinContent(k);
              val_noise=thisisnoise->GetBinContent(k);
              input[ch]->SetBinContent(k, val_in-val_noise);
            */
            //outPulse[ch][k] = inPulse[ch][k] - thisisnoise[k];
            outPulse[ch]->Fill(k, inPulse[ch]->GetBinContent(inPulse[ch]->FindBin(k)) - thisisnoise[k]);
          }
      }

    principal->Clear();
  }
  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) {
    TFile *output = new TFile(file,"READ");
    if (output->IsOpen()) {
      TH1D* uncalibrated = (TH1D*)output->Get("Uncalibrated_Pions"); 
      TF1 *landau = new TF1("landau","landau",-1000,15000);
      uncalibrated->Fit(landau,"0QR");
      return 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,-1000,15000); 
      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("data2011/TTrees/",""));
      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;
        }
        TCernOct11UnpackEvent* evnt = new TCernOct11UnpackEvent;
        TGo4EventElement* theBase=evnt;
        evnt->synchronizeWithTree(theTree, &theBase);
        Int_t entries = theTree->GetEntries();
        printf("%i Events found in TTree\n",entries);
        TMbsCrateEvent* fCrateInputEvent;
        TSpadicEvent* SpadicInputEvent;
        TSpadicData* theSpadic;
        Float_t Ch1(0), Ch2(0), Pb(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;
          Ch1 = fCrateInputEvent->fData1182[1][0];
          Ch2 = fCrateInputEvent->fData1182[0][1];
          Pb  = fCrateInputEvent->fData1182[1][5];
          getPID(Ch1, Ch2, Pb);
      
          //for (Int_t suid = 0; suid < /*usedSusibos*/1; suid++){
          theSpadic=dynamic_cast<TSpadicData*>(SpadicInputEvent->getEventElement(usedSuId));
          if (isPion){
            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);
            CancelNoise_Cova(baselinePulse, baselineNoisePulse, false);

            if(amplitude){
              if(minimumAmplitudeThreshold(baselineNoisePulse) && !BorderPadMax(baselineNoisePulse) && HitTimeWindow(baselineNoisePulse) && minimumIntegralThreshold(baselineNoisePulse))
                uncalibrated->Fill(Get3PadPosAmplitude(baselinePulse));
            }   
            else{
              if(minimumAmplitudeThreshold(baselineNoisePulse) && !BorderPadMax(baselineNoisePulse) && HitTimeWindow(baselineNoisePulse) && minimumIntegralThreshold(baselineNoisePulse))
                uncalibrated->Fill(Get3PadPosIntegral(baselinePulse));    
            }
          }
        }
        output->cd();
        uncalibrated->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(TString filename=/*"data2011/TTrees/Be_run2110020.root"*/"data2011/TTrees/Be_run2310004.root",TString radiator = "D", const Int_t suid = 1, Bool_t debug=false, Bool_t simple = false, Bool_t plotFromFile = true)
{
  gROOT->Reset(); 
  gROOT->SetStyle("Plain");
  gStyle->SetPalette(1,0);
  gStyle->SetOptStat(kFALSE);
  gStyle->SetOptTitle(kFALSE);
  Bool_t amplitude = false;
  clusterSize->SetXTitle("cluster width [channels]");
  clusterSize->SetYTitle("normalized counts");
  signalChDistance->SetXTitle("distance to next signal channel [channels]");
  signalChDistance->SetYTitle("normalized counts");
  if (suid > 7){
    printf("ERROR:: did not find susibo %i\n",suid);
    return;
  }
  const Int_t usedSusibos = 8;
  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};
  Double_t K3 = 0.525;
  Double_t par = 1.0;
  TString outpic;
  TString outfilename;
  outpic = filename;
  outpic.ReplaceAll("data2011/TTrees/","");
  outpic.ReplaceAll(".root","");
  outfilename.Form("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);
  if (MPV_pion == 1) {
    FillUnCalibratedPionSpectrum(filename,outfilename, usedSuId[suid], amplitude, debug);
    MPV_pion = Calibration(outfilename);
  }
  printf("MPV_pion: %.3e\n",MPV_pion);
  Double_t scaling2keV = 3.65338e+00/MPV_pion;//3.67179/MPV_pion;//7.87153e+02;//3.67179/6.67190e+02;
  Int_t nbins(200),lbin(-1000),hbin(15000);
  
  if(amplitude) {
    lbin /= 20;
    hbin /= 20;
  }
  else {
    nbins = 200;
    lbin = 0;
    hbin = 100;
  }
  TFile *output = new TFile(outfilename,"UPDATE");
  TH2F* Ch1_Pb;
  TH2F* noCh1_Pb;
  TH2F* Ch2_Pb;
  TH2F* noCh2_Pb;
  TH2F* Ch1_Ch2;
  TH2F* noCh1_Ch2;

  TH2D* prf;
  TProfile* prfProfile;

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

  TProfile *averagePulse_e;
  TProfile *averagePulse_p;

  TH1D *preCompSpectrum;
  TH1D *preCompBaselineSpectrum;
  TH1D *preCompSpectrumElectron;
  TH1D *preCompBaselineSpectrumElectron;
  TH1D *rawSpectrum;
  TH1D *baselineSpectrum;
  TH1D *noiseSpectrum;
  TH1D *baselineNoiseSpectrum;
  TH1D *baselineNoiseSpectrumWoOWoU;
  TH1D *baselineNoiseSpectrumElectron;
  TH1D *baselineNoiseSpectrumWoOWoUElectron;
  TH1D *noiseBaselineSpectrum;
  TH1D *probabilityOverflow;
  TH1D *probabilityUnderflow;
  TH1D *probabilityOverflowE;
  TH1D *probabilityUnderflowE;

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

  TH1D *baselineNoiseSpectrumCrossTalkE[6];
  TH1D *baselineNoiseSpectrumCrossTalkP[6];

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

  if (plotFromFile) {
    output->cd();
    Ch1_Pb = (TH2F*)output->Get("Ch1_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");
    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");
    if (!Ch1_Ch2)
      cout << "Ch1_Ch2" << endl;
    noCh1_Ch2 = (TH2F*)output->Get("noCh1_Ch2");
    if (!noCh1_Ch2)
      cout << "noCh1_Ch2" << endl;

    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;

    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 {
    Ch1_Pb = new TH2F("Ch1_Pb","good PID Ch1_Pb",2000,0,4000,2000,0,4000);
    Ch1_Pb->SetXTitle("Cherenkov 1 [a.u.]");
    Ch1_Pb->SetYTitle("Pb-glass calorimeter [a.u.]");
    Ch1_Pb->SetContour(99);

    noCh1_Pb = new TH2F("noCh1_Pb","no good PID Ch1_Pb",2000,0,4000,2000,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",2000,0,4000,2000,0,4000);
    Ch2_Pb->SetXTitle("Cherenkov 2 [a.u.]");
    Ch2_Pb->SetYTitle("Pb-glass calorimeter [a.u.]");
    Ch2_Pb->SetContour(99);

    noCh2_Pb = new TH2F("noCh2_Pb","no good PID Ch2_Pb",2000,0,4000,2000,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",2000,0,4000,2000,0,4000);
    Ch1_Ch2->SetXTitle("Cherenkov 1 [a.u.]");
    Ch1_Ch2->SetYTitle("Cherenkov 2 [a.u.]");
    Ch1_Ch2->SetContour(99);

    noCh1_Ch2 = new TH2F("noCh1_Ch2","no good PID Ch1_Ch2",2000,0,4000,2000,0,4000);
    noCh1_Ch2->SetXTitle("Cherenkov 1 [a.u.]");
    noCh1_Ch2->SetYTitle("Cherenkov 2 [a.u.]");
    noCh1_Ch2->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");

    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;
  
  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) )) )",K3,K3,K3,padWidth[suid],par,h[suid],K3,K3,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]);
  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);
 



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

    TFile inputFile(filename,"READ");
    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;
    }
    TCernOct11UnpackEvent* evnt = new TCernOct11UnpackEvent;
    TGo4EventElement* theBase=evnt;
    evnt->synchronizeWithTree(theTree, &theBase);
    Int_t entries = theTree->GetEntries();
    printf("%i Events found in TTree\n",entries);
    TMbsCrateEvent* fCrateInputEvent;
    TSpadicEvent* SpadicInputEvent;
    TSpadicData* theSpadic;
    Float_t Ch1(0), Ch2(0), Pb(0);
  
    if (debug) {
      c = new TCanvas("c","c",1750,700);
      c->Divide(7,2);
    }
    TCanvas *cDebug;
    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"));
      isPion = false;
      isElectron = false;
      isOverFlowEvent = false;
      isUnderFlowEvent = false;
      Ch1 = fCrateInputEvent->fData1182[1][0];
      Ch2 = fCrateInputEvent->fData1182[0][1];
      Pb  = fCrateInputEvent->fData1182[1][5];
      getPID(Ch1, Ch2, Pb);

      if (isPion || isElectron) {
        Ch1_Pb->Fill(Ch1,Pb);
        Ch2_Pb->Fill(Ch2,Pb);
      }
      if (!isPion && !isElectron) {
        noCh1_Pb->Fill(Ch1,Pb);
        noCh2_Pb->Fill(Ch2,Pb);
      }
      
      //for (Int_t suid = 0; suid < /*usedSusibos*/1; suid++){
      theSpadic=dynamic_cast<TSpadicData*>(SpadicInputEvent->getEventElement(usedSuId[suid]));
      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;
          
          }      
        }
        BaselineCompensation(preCompPulse, preCompBaselinePulse);
        BaselineCompensation(rawPulse, baselinePulse);
        CancelNoise_Cova(rawPulse, noisePulse, false);
        CancelNoise_Cova(baselinePulse, baselineNoisePulse, false);
        BaselineCompensation(noisePulse, noiseBaselinePulse);

        if(minimumAmplitudeThreshold(baselineNoisePulse) && !BorderPadMax(baselineNoisePulse) && HitTimeWindow(baselineNoisePulse) && minimumIntegralThreshold(baselineNoisePulse))
          CalcPrf(baselineNoisePulse, prf, prfProfile, padWidth[suid]);

        Int_t deltaChMax = -1;     
        Double_t temp = Get3PadPosIntegralCrossTalk(baselineNoisePulse, lastPulse, &deltaChMax);     
     
        if (isElectron) {
          nElectrons++;
          if(amplitude){         
            if(minimumAmplitudeThreshold(preCompPulse) && !BorderPadMax(preCompPulse) && HitTimeWindow(preCompPulse) && minimumIntegralThreshold(preCompPulse))
              preCompSpectrumElectron->Fill(Get3PadPosAmplitude(preCompPulse)*scaling2keV);
            if(minimumAmplitudeThreshold(preCompBaselinePulse) && !BorderPadMax(preCompBaselinePulse) && HitTimeWindow(preCompBaselinePulse) && minimumIntegralThreshold(preCompBaselinePulse))
              preCompBaselineSpectrumElectron->Fill(Get3PadPosAmplitude(preCompBaselinePulse)*scaling2keV);
            if(minimumAmplitudeThreshold(baselineNoisePulse) && !BorderPadMax(baselineNoisePulse) && HitTimeWindow(baselineNoisePulse) && minimumIntegralThreshold(baselineNoisePulse)){
              baselineNoiseSpectrumElectron->Fill(Get3PadPosAmplitude(baselineNoisePulse)*scaling2keV);
              if (isOverFlowEvent)
                probabilityOverflowE->Fill(Get3PadPosAmplitude(baselineNoisePulse));
              if(isUnderFlowEvent)
                probabilityUnderflowE->Fill(Get3PadPosAmplitude(baselineNoisePulse));
            }
          }
          else{
            if(minimumAmplitudeThreshold(preCompPulse) && !BorderPadMax(preCompPulse) && HitTimeWindow(preCompPulse) && minimumIntegralThreshold(preCompPulse))
              preCompSpectrumElectron->Fill(Get3PadPosIntegral(preCompPulse)*scaling2keV);
            if(minimumAmplitudeThreshold(preCompBaselinePulse) && !BorderPadMax(preCompBaselinePulse) && HitTimeWindow(preCompBaselinePulse) && minimumIntegralThreshold(preCompBaselinePulse))
              preCompBaselineSpectrumElectron->Fill(Get3PadPosIntegral(preCompBaselinePulse)*scaling2keV);
            if(minimumAmplitudeThreshold(baselineNoisePulse) && !BorderPadMax(baselineNoisePulse) && HitTimeWindow(baselineNoisePulse) && minimumIntegralThreshold(baselineNoisePulse)){
              baselineNoiseSpectrumElectron->Fill(Get3PadPosIntegral(baselineNoisePulse)*scaling2keV);
              FillAverageSignal(baselineNoisePulse, averagePulse_e);
              if (deltaChMax >= 0 && deltaChMax < 6){
                if(minimumAmplitudeThreshold(lastPulse) && !BorderPadMax(lastPulse) && HitTimeWindow(lastPulse) && minimumIntegralThreshold(lastPulse)){
                  baselineNoiseSpectrumCrossTalkE[deltaChMax]->Fill(temp*scaling2keV);
                  //crossTalk->Fill(deltaChMax,(temp*scaling2keV - Get3PadPosIntegral(baselineNoisePulse)*scaling2keV) / (temp*scaling2keV) * 100);
                }
              }
              if (isOverFlowEvent)
                probabilityOverflowE->Fill(Get3PadPosIntegral(baselineNoisePulse)*scaling2keV);
              if(isUnderFlowEvent)
                probabilityUnderflowE->Fill(Get3PadPosIntegral(baselineNoisePulse)*scaling2keV);      
              if (!isOverFlowEvent && !isUnderFlowEvent){
                if(amplitude){
                  baselineNoiseSpectrumWoOWoUElectron->Fill(Get3PadPosAmplitude(baselineNoisePulse));     
                }
                else{
                  baselineNoiseSpectrumWoOWoUElectron->Fill(Get3PadPosIntegral(baselineNoisePulse)*scaling2keV);           
                }
              }
            }
          }
        }
        if (isPion){
          nPions++;
          if(minimumAmplitudeThreshold(rawPulse) && !BorderPadMax(rawPulse) && HitTimeWindow(rawPulse) && minimumIntegralThreshold(rawPulse)){
            if(amplitude)
              rawSpectrum->Fill(Get3PadPosAmplitude(rawPulse));
            else
              rawSpectrum->Fill(Get3PadPosIntegral(rawPulse)*scaling2keV);
          }
          if(amplitude){
            if(minimumAmplitudeThreshold(baselinePulse) && !BorderPadMax(baselinePulse) && HitTimeWindow(baselinePulse) && minimumIntegralThreshold(baselinePulse))
              baselineSpectrum->Fill(Get3PadPosAmplitude(baselinePulse)*scaling2keV);
            if(minimumAmplitudeThreshold(noiseBaselinePulse) && !BorderPadMax(noiseBaselinePulse) && HitTimeWindow(noiseBaselinePulse) && minimumIntegralThreshold(noiseBaselinePulse))
              noiseBaselineSpectrum->Fill(Get3PadPosAmplitude(noiseBaselinePulse)*scaling2keV);
            if(minimumAmplitudeThreshold(noisePulse) && !BorderPadMax(noisePulse) && HitTimeWindow(noisePulse) && minimumIntegralThreshold(noisePulse))
              noiseSpectrum->Fill(Get3PadPosAmplitude(noisePulse)*scaling2keV); 
            if(minimumAmplitudeThreshold(preCompPulse) && !BorderPadMax(preCompPulse) && HitTimeWindow(preCompPulse) && minimumIntegralThreshold(preCompPulse))
              preCompSpectrum->Fill(Get3PadPosAmplitude(preCompPulse)*scaling2keV);
            if(minimumAmplitudeThreshold(preCompBaselinePulse) && !BorderPadMax(preCompBaselinePulse) && HitTimeWindow(preCompBaselinePulse) && minimumIntegralThreshold(preCompBaselinePulse))
              preCompBaselineSpectrum->Fill(Get3PadPosAmplitude(preCompBaselinePulse)*scaling2keV);
            if(minimumAmplitudeThreshold(baselineNoisePulse) && !BorderPadMax(baselineNoisePulse) && HitTimeWindow(baselineNoisePulse) && minimumIntegralThreshold(baselineNoisePulse)){
              baselineNoiseSpectrum->Fill(Get3PadPosAmplitude(baselineNoisePulse)*scaling2keV);
              if (!isOverFlowEvent && !isUnderFlowEvent)
                baselineNoiseSpectrumWoOWoU->Fill(Get3PadPosAmplitude(baselineNoisePulse));
              if (isOverFlowEvent)
                probabilityOverflow->Fill(Get3PadPosAmplitude(baselineNoisePulse));
              if(isUnderFlowEvent)
                probabilityUnderflow->Fill(Get3PadPosAmplitude(baselineNoisePulse));
            }
          }
          else{
            if(minimumAmplitudeThreshold(baselinePulse) && !BorderPadMax(baselinePulse) && HitTimeWindow(baselinePulse) && minimumIntegralThreshold(baselinePulse))
              baselineSpectrum->Fill(Get3PadPosIntegral(baselinePulse)*scaling2keV);
            if(minimumAmplitudeThreshold(noisePulse) && !BorderPadMax(noisePulse) && HitTimeWindow(noisePulse) && minimumIntegralThreshold(noisePulse))
              noiseSpectrum->Fill(Get3PadPosIntegral(noisePulse)*scaling2keV);
            if(minimumAmplitudeThreshold(noiseBaselinePulse) && !BorderPadMax(noiseBaselinePulse) && HitTimeWindow(noiseBaselinePulse) && minimumIntegralThreshold(noiseBaselinePulse))
              noiseBaselineSpectrum->Fill(Get3PadPosIntegral(noiseBaselinePulse)*scaling2keV);     
            if(minimumAmplitudeThreshold(preCompPulse) && !BorderPadMax(preCompPulse) && HitTimeWindow(preCompPulse) && minimumIntegralThreshold(preCompPulse))
              preCompSpectrum->Fill(Get3PadPosIntegral(preCompPulse)*scaling2keV);
            if(minimumAmplitudeThreshold(preCompBaselinePulse) && !BorderPadMax(preCompBaselinePulse) && HitTimeWindow(preCompBaselinePulse) && minimumIntegralThreshold(preCompBaselinePulse))
              preCompBaselineSpectrum->Fill(Get3PadPosIntegral(preCompBaselinePulse)*scaling2keV);
            if(minimumAmplitudeThreshold(baselineNoisePulse) && !BorderPadMax(baselineNoisePulse) && HitTimeWindow(baselineNoisePulse) && minimumIntegralThreshold(baselineNoisePulse)){
              baselineNoiseSpectrum->Fill(Get3PadPosIntegral(baselineNoisePulse)*scaling2keV);
              FillAverageSignal(baselineNoisePulse, averagePulse_p);
              if (deltaChMax >= 0 && deltaChMax < 6){
                if(minimumAmplitudeThreshold(lastPulse) && !BorderPadMax(lastPulse) && HitTimeWindow(lastPulse) && minimumIntegralThreshold(lastPulse)){
                  baselineNoiseSpectrumCrossTalkP[deltaChMax]->Fill(temp*scaling2keV);
                  crossTalk->Fill(deltaChMax,(temp*scaling2keV - Get3PadPosIntegral(baselineNoisePulse)*scaling2keV) / (temp*scaling2keV) * 100);
                }
              }
              if (!isOverFlowEvent && !isUnderFlowEvent)
                baselineNoiseSpectrumWoOWoU->Fill(Get3PadPosIntegral(baselineNoisePulse)*scaling2keV);
              if (isOverFlowEvent)
                probabilityOverflow->Fill(Get3PadPosIntegral(baselineNoisePulse)*scaling2keV);
              if(isUnderFlowEvent)
                probabilityUnderflow->Fill(Get3PadPosIntegral(baselineNoisePulse)*scaling2keV);
            }
          
          }
        }
        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(/*"Pions_dEdx_3GeV_12mm_XeCo2.root"*/"root/Sim_xeco20_12mm_03GeV.root","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./(Float_t)simSpectrum->GetEntries());
  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);
  leg->SetHeader("#pi^{-} spectra, " + usedSuName[suid] + ", p=3GeV/c, 12mm XeCO_{2}");
  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./(Float_t)/*probabilityOverflow*/baselineNoiseSpectrum->GetEntries());
    probabilityUnderflow->Scale(1./(Float_t)/*probabilityUnderflow*/baselineNoiseSpectrum->GetEntries());
  }
  probabilityUnderflow->DrawCopy("same");
  probabilityOverflow->DrawCopy("same");
  //rawSpectrum->DrawCopy();
  //baselineSpectrum->DrawCopy("same");
  if (!plotFromFile)
    baselineNoiseSpectrum->Scale(1./(Float_t)baselineNoiseSpectrum->GetEntries());
  baselineNoiseSpectrum->DrawCopy("same");
  if (!plotFromFile)
    baselineNoiseSpectrumWoOWoU->Scale(1./(Float_t)baselineNoiseSpectrumWoOWoU->GetEntries());
  baselineNoiseSpectrumWoOWoU->DrawCopy("same");
  //noiseSpectrum->DrawCopy("same");
  //noiseBaselineSpectrum->Scale(1./(Float_t)noiseBaselineSpectrum->GetEntries());
  //noiseBaselineSpectrum->DrawCopy("same");
  if (!simple) {
    if (!plotFromFile)
      preCompSpectrum->Scale(1./(Float_t)preCompSpectrum->GetEntries());
    preCompSpectrum->DrawCopy("same");
    if (!plotFromFile)
      preCompBaselineSpectrum->Scale(1./(Float_t)preCompBaselineSpectrum->GetEntries());
    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"*/"root/Sim_xeco20_12mm_03GeV.root","READ");
  TFile TR_simFile("root/Sim_xeco20_12mm_03GeV_TR.root","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 == "B" || radiator == "C" || radiator == "D" || radiator == "E" || radiator == "F"){
      spectrumAbs = (TH1D*)TR_simFile.Get("Absorption_("+radiator+")");
      spectrum = (TH1D*)TR_simFile.Get("Production_("+radiator+")");
      if (radiator == "B")
        attenuation = 0.55;//ok
      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.55;//ok
    }
    else if (radiator == "A" || radiator == "G10" || radiator == "G20" || radiator == "G30" || radiator == "H" || radiator == "H++" || radiator == "Iu"  || radiator == "Ip" || radiator == "J"){
      TString irradiator;
      if (radiator == "A"){
        irradiator = "C";
        attenuation = 0.35;//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.65;//ok
      }
      if (radiator == "H"){
        irradiator = "C";
        attenuation = 0.65;//ok
      }
      if (radiator == "H++"){
        irradiator = "C";
        attenuation = 0.65;//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
      }
      spectrumAbs = (TH1D*)TR_simFile.Get("Absorption_("+irradiator+")");
      spectrum = (TH1D*)TR_simFile.Get("Production_("+irradiator+")");
    }
    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();
    TRandom *rEfficiency = new TRandom();
 

    
    radiator2.Form("%s (%.1f%% attenuation)",radiator.Data(),attenuation*100);
    for (Int_t iElec = 0; iElec < 200000; iElec++) {
      Int_t nPhot = rEfficiency->Poisson(spectrum->Integral() * spectrum->GetBinWidth(1) * attenuation);
      Float_t dEdx = hElDedx->GetRandom();
      Float_t Tr = 0;
      for (Int_t iPhot = 0; iPhot < nPhot; iPhot++) {
        Tr += spectrum->GetRandom();
      }
      hElTr->Fill(dEdx+Tr);
    }
    for (Int_t iElec = 0; iElec < 200000; iElec++) {
      Int_t nPhot = rEfficiency->Poisson(spectrumAbs->Integral() * spectrumAbs->GetBinWidth(1) * attenuation);
      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./(Float_t)simSpectrumE->GetEntries());
  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./(Float_t)/*probabilityOverflow*/baselineNoiseSpectrumElectron->GetEntries());
    probabilityUnderflowE->Scale(1./(Float_t)/*probabilityUnderflow*/baselineNoiseSpectrumElectron->GetEntries());
  }
  probabilityUnderflowE->DrawCopy("same");
  probabilityOverflowE->DrawCopy("same");
  if (!plotFromFile)
    baselineNoiseSpectrumElectron->Scale(1./(Float_t)baselineNoiseSpectrumElectron->GetEntries());
  baselineNoiseSpectrumElectron->DrawCopy("same");
  if (!plotFromFile)
    baselineNoiseSpectrumWoOWoUElectron->Scale(1./(Float_t)baselineNoiseSpectrumWoOWoUElectron->GetEntries());
  baselineNoiseSpectrumWoOWoUElectron->DrawCopy("same");
  if (!simple) {
    if (!plotFromFile)
      preCompSpectrumElectron->Scale(1./(Float_t)preCompBaselineSpectrumElectron->GetEntries());
    preCompSpectrumElectron->DrawCopy("same");
    if (!plotFromFile)
      preCompBaselineSpectrumElectron->Scale(1./(Float_t)preCompBaselineSpectrumElectron->GetEntries());
    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.525","");
  prfProfile->Fit("mathiesonFit","QR0");
  values.Form(" K_{3}=%.3f",mathiesonFit->GetParameter(0));
  lPRF->AddEntry(mathiesonFit,"mathieson fit","l");
  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* hEfficiency = new TProfile("hEfficiency","hEfficiency",12,0.5,12.5);
  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}");

  //TFile *output = new TFile(outfilename,"UPDATE");
  output->cd();

  hEfficiency->GetXaxis()->SetRangeUser(1e-2,1e2);
  CalcPionEfficiency(simSpectrumE, simSpectrum, hEfficiency, nChambers, debug);
  hEfficiency->SetMarkerColor(simSpectrumE->GetLineColor());
  hEfficiency->SetLineColor(simSpectrumE->GetLineColor());
  output->cd();
  hEfficiency->Write("Efficiency_simulated", TObject::kOverwrite);
  leg3->AddEntry(simSpectrum,"simulated","l");
  cEfficiency->cd()->SetLogy(1);
  hEfficiency->DrawCopy("P,text0");
  CalcPionEfficiency(preCompBaselineSpectrumElectron, preCompBaselineSpectrum, hEfficiency, nChambers, debug);
  hEfficiency->SetMarkerColor(preCompBaselineSpectrumElectron->GetLineColor());
  hEfficiency->SetLineColor(preCompBaselineSpectrumElectron->GetLineColor());
  output->cd();
  hEfficiency->Write("Efficiency_measured_MS", TObject::kOverwrite);
  leg3->AddEntry(preCompBaselineSpectrumElectron,"measured: MS noise + baseline","l");
  cEfficiency->cd();
  hEfficiency->DrawCopy("P,text0,same");
  CalcPionEfficiency(baselineNoiseSpectrumElectron, baselineNoiseSpectrum, hEfficiency, nChambers, debug);
  hEfficiency->SetMarkerColor(baselineNoiseSpectrumElectron->GetLineColor());
  hEfficiency->SetLineColor(baselineNoiseSpectrumElectron->GetLineColor());
  output->cd();
  hEfficiency->Write("Efficiency_measured_FFM", TObject::kOverwrite);
  leg3->AddEntry(baselineNoiseSpectrumElectron,"measured: FFM noise + baseline","l");
  cEfficiency->cd();
  hEfficiency->DrawCopy("P,text0,same");
  CalcPionEfficiency(baselineNoiseSpectrumWoOWoUElectron, baselineNoiseSpectrumWoOWoU, hEfficiency, nChambers, debug);
  hEfficiency->SetMarkerColor(baselineNoiseSpectrumWoOWoUElectron->GetLineColor());
  hEfficiency->SetLineColor(baselineNoiseSpectrumWoOWoUElectron->GetLineColor());
  output->cd();
  hEfficiency->Write("Efficiency_measured_FFM_WoOWo", TObject::kOverwrite);
  leg3->AddEntry(baselineNoiseSpectrumWoOWoUElectron,"measured: FFM noise + baseline (woO && woU)","l");
  cEfficiency->cd();
  hEfficiency->DrawCopy("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");

  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",800,600); 
  r->Divide(1,1);
  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->SaveAs(outpic+"PulseShape.pdf");
  r->SaveAs(outpic+"PulseShape.png");

  TCanvas* cBeamMonitor = new TCanvas("BeamMonitor","Beam Monitor",2*800,600);
  cBeamMonitor->Divide(2,1);
  cBeamMonitor->cd(1)->SetLogz(1);
  Ch1_Pb->DrawCopy("colz");
  cBeamMonitor->cd(2)->SetLogz(1);
  Ch2_Pb->DrawCopy("colz");
  /*//Add cut lines and statistic legend here

   */

  cBeamMonitor->SaveAs(outpic+"Ch_Pb.pdf");
  cBeamMonitor->SaveAs(outpic+"Ch_Pb.png");
  if (!plotFromFile) {
    printf("%i Pion & %i Electron events found\n",nPions,nElectrons);
    Ch1_Pb->Write("", TObject::kOverwrite);
    noCh1_Pb->Write("", TObject::kOverwrite);
    Ch2_Pb->Write("", TObject::kOverwrite);
    noCh2_Pb->Write("", TObject::kOverwrite);
    Ch1_Ch2->Write("", TObject::kOverwrite);
    noCh1_Ch2->Write("", TObject::kOverwrite);
    prfProfile->Write("", TObject::kOverwrite);
    prf->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);
    probabilityOverflow->Write("", TObject::kOverwrite);
    probabilityUnderflow->Write("", TObject::kOverwrite);
    probabilityOverflowE->Write("", TObject::kOverwrite);
    probabilityUnderflowE->Write("", TObject::kOverwrite);
    for(Int_t delta = 0; delta < 6; delta++) {
      baselineNoiseSpectrumCrossTalkE[delta]->Scale(1./(Float_t)baselineNoiseSpectrumCrossTalkE[delta]->GetEntries());
      baselineNoiseSpectrumCrossTalkE[delta]->Write("", TObject::kOverwrite);
      baselineNoiseSpectrumCrossTalkP[delta]->Scale(1./(Float_t)baselineNoiseSpectrumCrossTalkP[delta]->GetEntries());
      baselineNoiseSpectrumCrossTalkP[delta]->Write("", TObject::kOverwrite);  
    }
    //deltaSim->Write("", TObject::kOverwrite);
    //deltaSim2->Write("", TObject::kOverwrite);
    crossTalk->Write("", TObject::kOverwrite);
    cResults->Write("", TObject::kOverwrite);
  }
  output->Close();
}