beamtime/cern-oct12/go4/AnalysisMacros/readtree_spadic.cxx (r4864/r4843)

Go to the documentation of this file.

// Test script to retrieve complete composite event from Go4 written ROOT Tree
// May be used as template for user analysis in root
// v0.1 16-Jun-2011 J.Adamczewski-Musch, GSI - initial version
// v0.2 22-Jun-2011 S.Linev, GSI - reading of complete event
// v0.3 7-Nov-2011 JAM - show how to access single spadic data from beamtime event

// Test file can be easily produced by analysis of old lmd files (from CERN beamtime in Nov 2010)
// [shell]  go4analysis -file nov20_run44_0000.lmd -store test.root

// NOTE; you need to run this script with ACLIC and set include path!
// The following does not work?

//#ifndef __CINT__

#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 "TH2.h"
#include "TH1.h"
#include "TF1.h"

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

#include "TFile.h"
#include "TTree.h"
#include "TBranch.h"
#include "TMath.h"
#include "TLegend.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 "langaus.C"

//#endif

// this is workaround for interpreter mode:
#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;
TPrincipal* principal = new TPrincipal(NUM_SPADIC_CHA,"ND");
Int_t pidcuts[18] = {0};
using namespace std;

class MFSpadicEvent
{
public:
  Float_t rawPulse[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE];
  Float_t cleanPulse[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE];

  Bool_t signalCH[NUM_SPADIC_CHA];

  Bool_t underFlowEv;
  Bool_t underFlowCH[NUM_SPADIC_CHA];
  Bool_t underFlowTB[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE];

  Bool_t overFlowEv;
  Bool_t overFlowCH[NUM_SPADIC_CHA];
  Bool_t overFlowTB[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE];


  MFSpadicEvent() {
    for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {
      for (Int_t tb = 0; tb < SPADIC_TRACE_SIZE; tb++) {
        rawPulse[ch][tb] = 0.0;
        cleanPulse[ch][tb] = 0.0;
        underFlowTB[ch][tb] = kFALSE;
        overFlowTB[ch][tb] = kFALSE;
      }
      signalCH[ch] = kFALSE;
      underFlowCH[ch] = kFALSE;
      overFlowCH[ch] = kFALSE;
    }
    underFlowEv = kFALSE;
    overFlowEv = kFALSE;
    /*
      rawPulse[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE] = {0.0};
      cleanPulse[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE] = {0.0};
      signalCH[NUM_SPADIC_CHA] = {kFALSE};
      underFlowEv = kFALSE;
      underFlowCH[NUM_SPADIC_CHA] = {kFALSE};
      underFlowTB[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE] = {kFALSE};
      overFlowEv = kFALSE;
      overFlowCH[NUM_SPADIC_CHA] = {kFALSE};
      overFlowTB[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE] = {kFALSE};
    */
  }

private:

};

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 Clusterizer(Double_t pulse[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE],Int_t nSpadicSignalCh, Bool_t SignalChannel[NUM_SPADIC_CHA], Int_t maxch, Int_t maxT[NUM_SPADIC_CHA], Double_t maxA[NUM_SPADIC_CHA]) {
  Bool_t Run_TimeClustering(0), Run_AmplitudeClustering(0);
  nSpadicSignalCh = 0;
  /*
  if (Run_TimeClustering) {
    for (UInt_t ch = 0; ch < NUM_SPADIC_CHA; ch++){
      if (maxT[ch] >= maxT[maxch] - 5 && maxT[ch] <= maxT[maxch] + 5){ // if signal is correlated in time domain
        //spdata->fSpadicSignalCh[ch] = true;
        SignalChannel[ch] = true;
        nSpadicSignalCh++;
      }
      else {
        SignalChannel[ch] = false;
      }
    }
  }
  if (Run_AmplitudeClustering) {
    for (UInt_t ch = 0; ch < NUM_SPADIC_CHA; ch++){
      if (maxA[ch] > 15 + pulse[ch][0]){                               //   if signal is above a minimum threshold with respect to the first time bin of this channel
        //spdata->fSpadicSignalCh[ch] = true;
        SignalChannel[ch] = true;
        nSpadicSignalCh++;        
      }
      else {
        SignalChannel[ch] = false;
      }
    }
  }
  fSpadic_ClusterWidth[sid]->Fill(nSpadicSignalCh);  
  */
}

void NoiseReduction(TH1D* inPulse[NUM_SPADIC_CHA], TH1D* outPulse[NUM_SPADIC_CHA], Bool_t SignalChannel[NUM_SPADIC_CHA]) {
  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;
    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 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)
        {
          SignalChannel[ch] = false;
          noise_ch_counter++;
          //thisisnoise->Add(input[ch]);
          for(Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++){
            thisisnoise[bin] += inputarray[bin][ch];
          }
        }
      else { 
        SignalChannel[ch] = true;
      } 
    }
  //thisisnoise->Scale(1./(double)noise_ch_counter);
  for(Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++){
    thisisnoise[bin] /= (double)noise_ch_counter;
  }
 
  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();

  /*
    Double_t AverTrace[SPADIC_TRACE_SIZE] = {0.0};
    for (UInt_t bin = 0; bin < SPADIC_TRACE_SIZE-1; bin++){ 
    AverTrace[bin] = 0.0;  
    for (UInt_t ch = 0; ch < NUM_SPADIC_CHA; ch++){
    if (!SignalChannel[ch]){
    AverTrace[bin] += rawPulse[ch]bin];
    }
  
    if(nSpadicSignalCh<NUM_SPADIC_CHA) // Debug check validity of array? JAM
    {
    AverTrace[bin] /= (NUM_SPADIC_CHA - nSpadicSignalCh);
    }
    else
    {
    //cout <<"WWWWWWWWWW arning: nSpadicSignalCh="<<nSpadicSignalCh<<" exceeds NUM_SPADIC_CHA "<<NUM_SPADIC_CHA << endl;
    }

    for (UInt_t ch = 0; ch < NUM_SPADIC_CHA; ch++){ 
    cleanPulse[ch][bin] = rawPulse[ch]bin] - AverTrace[bin];
    }
  */
}
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++){
    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 BaselineCorrection(Double_t pulse[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE]) {
  for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {
    //cout << endl;
    Float_t baseline(0);        
    for (Int_t bin = 0; bin < noBaselineTB; bin++) {    
        baseline += pulse[ch][bin];
    }
    baseline /= Float_t(noBaselineTB);
    for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {       
      pulse[ch][bin] -= baseline;
    }
  }
}
/*
void CopyData(Double_t pulse[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE], TSpadicData* theSpadicA) {
  for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {       
    for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {       
      pulse[ch][bin] = theSpadicA->fSpadicCompensated[ch][bin];
    }
  }
  BaselineCorrection(pulse);
}
void MergeData(Double_t pulse[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE], TSpadicData* theSpadicA, TSpadicData* theSpadicB) {
  Double_t pulseA[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE] = {{0.0}};
  Double_t pulseB[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE] = {{0.0}};
  CopyData(pulseA, theSpadicA);
  CopyData(pulseB, theSpadicB);
  for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {       
    for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {
      pulse[ch][bin] = pulseA[ch][bin] + pulseB[7-ch][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 GetPidCuts(TString runId)
{
  Bool_t MS_Pid = true;
  Bool_t FFM_Pid = false;
  //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] =  500; // // El Ch1 min
  pidcuts[1] = 9999; // // El Ch1 max
  pidcuts[2] =  450; // // El Ch2 min
  pidcuts[3] = 9999; // // El Ch2 max
  pidcuts[4] = 1050; // // El Pb min
  pidcuts[5] = 9999; // // El Pb max
  pidcuts[6] =    0; // // Pi Ch1 min
  pidcuts[7] =  375; // // Pi Ch1 max
  pidcuts[8] =    0; // // Pi Ch2 min
  pidcuts[9] =  340; // // Pi Ch2 max
  pidcuts[10]=  500;//350; // // Pi Pb min
  pidcuts[11]= 1000; // // Pi Pb max
  pidcuts[12]=  400; // // Mu Ch1 min
  pidcuts[13]=  700; // // Mu Ch1 max
  pidcuts[14]=  360; // // Mu Ch2 min
  pidcuts[15]=  450; // // Mu Ch2 max
  pidcuts[16]=  500; // // Mu Pb min
  pidcuts[17]= 1000; // // Mu Pb max
  /*
  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 (MS_Pid) {
    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 (FFM_Pid){
  
    if (runnumber < 1910000 )       cout << "[WARN] No PID cuts set! Invalid runnumber " << runnumber << endl;
    else if (runnumber >= 2010000 && runnumber <= 2010008) {//3GeV
      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 > 2110008 && runnumber <= 2110015) {//3GeV
      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) {//3GeV
      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[12]=  1; pidcuts[13]=   1; pidcuts[14]=  1; pidcuts[15]=   1; pidcuts[16]=   1; pidcuts[17]=   1;
    }
    else if (runnumber > 2110020 && runnumber <= 2210000) {//3GeV
      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) {//3GeV
      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) {//3GeV
      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) {//3GeV
      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 <= 2310000) {//3GeV
      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) {//3GeV
      pidcuts[0] =600; pidcuts[1] =1500; pidcuts[2] =550; pidcuts[3] =1500; 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 > 2310001 && runnumber <= 2310004) {//3GeV
      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) {//3GeV
      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 <= 2410010) {//3GeV
      pidcuts[0] = 500; pidcuts[1] =9999; pidcuts[2] =550; pidcuts[3] =9999; pidcuts[4] =1000; pidcuts[5] =9999; 
      pidcuts[6] = 300; pidcuts[7] = 340; pidcuts[8] =315; pidcuts[9] = 375; 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 > 2410010 && runnumber <= 2410020) {//2GeV
      pidcuts[0] = 500; pidcuts[1] =9999; pidcuts[2] =550; pidcuts[3] =9999; pidcuts[4] =1000; 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 > 2410020 && runnumber <= 2510003) {//4GeV
      pidcuts[0] = 600; 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 > 2510003 && runnumber <= 2510012) {//6GeV
      pidcuts[0] = 500; 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 == 2510013) {//10GeV
      pidcuts[0] = 400; pidcuts[1] =9999; pidcuts[2] =500; 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 > 2510014 && runnumber <= 2610005) {//8GeV
      pidcuts[0] = 500; 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 > 2610005 && runnumber <= 2710019) {//4GeV
      pidcuts[0] = 500; 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 > 2710019 && runnumber <= 2810003) {//2GeV
      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 > 2810003 && runnumber <= 2910013) {//4GeV
      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 > 2910016 && runnumber <= 3010015) {//2GeV
      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 > 2910013 && 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;}
    
  
  
    //  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 > 2110000 && runnumber <= 2110015) {
      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 > 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[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] =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 > 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] =3500; 
      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;
    }  
}
/*
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;
}
*/
/*
void Distance2Track(Double_t xPos[4], Double_t distance[4]) {
  Double_t zMean(0.0), xMean(0.0), m(0.0), b(0.0), SUMzx(0.0), SUMzz(0.0);
  //Least squares: x = m * z + b
  Double_t zPos[4] = {0, 320, 638, 958};
  for (Int_t j = 0; j < 4; j++) {
    for (Int_t i = 0; i < 4; i++) {
      if (i != j) {
        zMean += zPos[i];
        xMean += xPos[i];
      }
    }
    zMean /= 3.;
    xMean /= 3.;
    for (Int_t i = 0; i < 4; i++) {
      if (i != j) {
        SUMzx += (zPos[i] - zMean) * (xPos[i] - xMean);
        SUMzz += (zPos[i] - zMean) * (zPos[i] - zMean);
      }
    }
    m = SUMzx / SUMzz;
    b = xMean - m * zMean;
    distance[j] = xPos[j] - (m * zPos[j] + b);
  }
}
*/
Bool_t Distance2TrackFFM(Double_t x[8], Double_t distance[8]/*, Double_t mAlig[8][8], Double_t bAlig[8][8] /*, Bool_t MS = true*/) {
  //Alignment correction
  // use suid = 0 as fixed reference
  Double_t xPos[8] = {0.0};
  Int_t nPos = 0;
  for (Int_t i = 4; i< 8; i++) {
    xPos[i] = x[i];
    if (xPos[i] != 0.0) nPos++;
  }
  if (nPos < 4) return false;
  
  Double_t zMean(0.0), xMean(0.0), m(0.0), b(0.0), SUMzx(0.0), SUMzz(0.0);
  //Least squares fit: x = m * z + b
  //                  MS  Ms   Ms   MS  F    F    F    F
  Double_t zPos[8] = {0, 320, 638, 958, 0, 320, 638, 958};
  for (Int_t j = 4; j < 8; j++) {
    Int_t count = 0;
    for (Int_t i = 4; i < 8; i++) {
      if (i != j && x[i] != 0.0) {
        zMean += zPos[i];
        xMean += xPos[i];
        count++;
      }
    }
    //if (count < 2) return false;
    zMean /= count;
    xMean /= count;
    for (Int_t i = 4; i < 8; i++) {
      if (i != j && x[i] != 0.0) {
        SUMzx += (zPos[i] - zMean) * (xPos[i] - xMean);
        SUMzz += (zPos[i] - zMean) * (zPos[i] - zMean);
      }
    }
    m = SUMzx / SUMzz;
    b = xMean - m * zMean;
    distance[j] = xPos[j] - (m * zPos[j] + b);
  }
  return true;
}
Bool_t Distance2TrackMS(Double_t x[8], Double_t distance[8]/*, Double_t mAlig[8][8], Double_t bAlig[8][8] /*, Bool_t MS = true*/) {
  //Alignment correction
  // use suid = 0 as fixed reference
  Double_t xPos[4] = {0.0};
  Int_t nPos = 0;
  for (Int_t i = 0; i< 4; i++) {
    xPos[i] = x[i];
    if (xPos[i] != 0.0) nPos++;
  }
  if (nPos < 4) return false;
 
  Double_t zMean(0.0), xMean(0.0), m(0.0), b(0.0), SUMzx(0.0), SUMzz(0.0);
  //Least squares fit: x = m * z + b

  Double_t zPos[4] = {0, 320, 638, 958};
  for (Int_t j = 0; j < 4; j++) {
    Int_t count = 0;
    for (Int_t i = 0; i < 4; i++) {
      if (i != j && x[i] != 0.0) {
        zMean += zPos[i];
        xMean += xPos[i];
        count++;
      }
    }
    //if (count < 2) return false;
    zMean /= count;
    xMean /= count;
    for (Int_t i = 0; i < 4; i++) {
      if (i != j && x[i] != 0.0) {
        SUMzx += (zPos[i] - zMean) * (xPos[i] - xMean);
        SUMzz += (zPos[i] - zMean) * (zPos[i] - zMean);
      }
    }
    m = SUMzx / SUMzz;
    b = xMean - m * zMean;
    distance[j] = xPos[j] - (m * zPos[j] + b);
  }
  return true;
}


void readtree_spadic(TString nameFileSet, Bool_t debug=0) {

  Bool_t raw = true;
  Bool_t msMerge = false;

  ifstream fileSet;
  fileSet.open(nameFileSet.Data());

  if (fileSet.is_open()) {
    TString ofname = "data2012/Analysis/" + TString(nameFileSet(0,nameFileSet.Sizeof()-5)).Append(".root");
    TFile outputFile(ofname,"RECREATE","output of readtree_spadic.cxx");
    TString mergeName[2] = {"NotMerged","Merged"};
    TString filterName[4] = {"wOwU","woOwU","woOwoU","wOwoU"}; // w == with; wo == without; O == overflow; U == undershoot
    TString folder = "data2012/TTrees/";
    TString picEx = ".pdf";
    const Int_t usedSusibos = 3;
    const Int_t recoSets = 9;
    const Int_t nMerge = 2;
    const Int_t nFilter = 4;
    //const Int_t noBaselineTB = 5;
    Int_t topSuId[3]            = {        11,         10,       18};
    Int_t usedSuId[usedSusibos] = {        11,         10,       18};//{         2,         12,        6,        18,         4,         5,        16,        17  };
    Float_t padWidth[usedSusibos] = {7.125, 7.125, 7.125};//{       5,          5,        8,         8,         5,         5,         5,         5  };
    TString recoSetName[recoSets] = {"ClusterInteg","3PadInteg","1PadInteg","ClusterIntegWindow","3PadIntegWindow","1PadIntegWindow","ClusterAmpl","3PadAmpl","1PadAmpl"};
    Int_t allPion[usedSusibos][nMerge] = {{0}};
    Int_t allElectron[usedSusibos][nMerge] = {{0}};
    Int_t allGoodEvent[usedSusibos][nMerge] = {{0}};
    Int_t allEvents = 0;
    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,W[a],par,h[b],K3,K3,W[a],par,h[b]);

    TF1 *mathiesonMS336 = new TF1("mathieson",
                                  " -1. / (2. * atan(sqrt([0]))) * (atan(sqrt([0]) *tanh(3.14159265 * (-2. + sqrt([0]) ) * (5.0 + 2.* x * 1.0) / (8.* 3.0) )) +  atan(sqrt([0]) *  tanh(3.14159265 * (-2. + sqrt([0]) ) * (5.0 - 2.* x * 1.0) / (8.* 3.0) )) )"
                                  ,-10,10);
    mathiesonMS336->SetLineStyle(2);
    mathiesonMS336->SetLineColor(3);
    TF1 *mathiesonMS444 = new TF1("mathieson",
                                  " -1. / (2. * atan(sqrt([0]))) * (atan(sqrt([0]) *tanh(3.14159265 * (-2. + sqrt([0]) ) * (8.0 + 2.* x * 1.0) / (8.* 4.0) )) +  atan(sqrt([0]) *  tanh(3.14159265 * (-2. + sqrt([0]) ) * (8.0 - 2.* x * 1.0) / (8.* 4.0) )) )"
                                  ,-10,10);
    mathiesonMS444->SetLineStyle(2);
    mathiesonMS444->SetLineColor(3);
    TF1 *gausFit = new TF1("gausFit","gaus",20,50);
    TF1 *sigmaFit = new TF1("sigmaFit","gaus",-10,10);
    TF1 *alignmentFit= new TF1("alignmentFit","[0]*x+[1]",0,60);
    // Histogramms
    TString name, title;

    //TH1D* hCompensatedAll[usedSusibos];
    TH2F* 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);

    TH2F* 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);

    TH2F* 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);

    TH2F* 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);

    TH2F* 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);

    TH2F* 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);

    TH1D* preCompSpectraEl[usedSusibos][recoSets][nMerge][nFilter];
    TH1D* preCompSpectraPi[usedSusibos][recoSets][nMerge][nFilter];
    TH1I* hitTimeEl[usedSusibos];
    TH1I* hitTimePi[usedSusibos];

    TH2I* amplitudeCorr[usedSusibos][usedSusibos];

    TH1I* amplitudeSpectra[usedSusibos];

    TH1F* chaGain[usedSusibos];

    TProfile* pulseShapeEl[usedSusibos];
    TProfile* pulseShapePi[usedSusibos];

    TH1I* beamProfile[usedSusibos];
    TProfile* chamberAlignment[usedSusibos][usedSusibos]; 

    TH2D* prf[usedSusibos];
    TProfile* prfProfile[usedSusibos];

    TH1I *chaGainCalibSpectra[usedSusibos][NUM_SPADIC_CHA];

    for (Int_t suid = 0; suid < usedSusibos; suid++) {
      for (Int_t rs = 0; rs < recoSets; rs++) {
        for (Int_t imerge = 0; imerge < nMerge; imerge++) {
          for (Int_t ifilter = 0; ifilter < nFilter; ifilter++) {
            title.Form("preCompSpectra%sEl_%i_%s_%s",recoSetName[rs].Data(),usedSuId[suid],mergeName[imerge].Data(),filterName[ifilter].Data());
            name.Form("pre cleaned go4 output (%s) electron Sus%i (%s %s)",recoSetName[rs].Data(),usedSuId[suid],mergeName[imerge].Data(),filterName[ifilter].Data());
            preCompSpectraEl[suid][rs][imerge][ifilter] = new TH1D(title, name, NUM_SPADIC_CHA*SPADIC_TRACE_SIZE*SPADIC_ADC_MAX,0,NUM_SPADIC_CHA*SPADIC_TRACE_SIZE*SPADIC_ADC_MAX);
            preCompSpectraEl[suid][rs][imerge][ifilter]->SetXTitle("charge [a.u.]");
            preCompSpectraEl[suid][rs][imerge][ifilter]->SetYTitle("normalized counts");
            preCompSpectraEl[suid][rs][imerge][ifilter]->SetLineColor(2);

            title.Form("preCompSpectra%sPi_%i_%s_%s",recoSetName[rs].Data(),usedSuId[suid],mergeName[imerge].Data(),filterName[ifilter].Data());
            name.Form("pre cleaned go4 output (%s) pion Sus%i (%s %s)",recoSetName[rs].Data(),usedSuId[suid],mergeName[imerge].Data(),filterName[ifilter].Data());
            preCompSpectraPi[suid][rs][imerge][ifilter] = new TH1D(title, name, NUM_SPADIC_CHA*SPADIC_TRACE_SIZE*SPADIC_ADC_MAX,0,NUM_SPADIC_CHA*SPADIC_TRACE_SIZE*SPADIC_ADC_MAX);
            preCompSpectraPi[suid][rs][imerge][ifilter]->SetXTitle("charge [a.u.]");
            preCompSpectraPi[suid][rs][imerge][ifilter]->SetYTitle("normalized counts");
          }
        }
      }
      for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {
        title.Form("chaGainCalibSpectraCh%i_%i",ch,usedSuId[suid]);
        name.Form("Channel Spectra Ch%i Sus%i",ch,usedSuId[suid]);
        chaGainCalibSpectra[suid][ch] = new TH1I(title, name, SPADIC_TRACE_SIZE*SPADIC_ADC_MAX+50,-50,SPADIC_TRACE_SIZE*SPADIC_ADC_MAX);
        chaGainCalibSpectra[suid][ch]->SetLineColor(ch+1);
      }
      for (Int_t suid2 = 0; suid2 < usedSusibos; suid2++) {
        title.Form("amplitudeCorr_%i_%i",usedSuId[suid],usedSuId[suid2]);
        name.Form("max amplitude correlation (Sus%i,Sus%i)",usedSuId[suid],usedSuId[suid2]);
        amplitudeCorr[suid][suid2] = new TH2I(title, name, SPADIC_ADC_MAX+50,-50,SPADIC_ADC_MAX, SPADIC_ADC_MAX+50,-50,SPADIC_ADC_MAX);
        amplitudeCorr[suid][suid2]->SetContour(99); 

        title.Form("chamberAlignment_%i_%i",usedSuId[suid],usedSuId[suid2]);
        name.Form("chamber alignment Sus%i Sus%i",usedSuId[suid],usedSuId[suid2]);
        chamberAlignment[suid][suid2] = new TProfile(title, name, NUM_SPADIC_CHA*padWidth[suid], 0-0.5, NUM_SPADIC_CHA*padWidth[suid]-0.5*padWidth[suid]);
        chamberAlignment[suid][suid2]->GetYaxis()->SetRangeUser(-(NUM_SPADIC_CHA-0.5)*padWidth[suid2],(NUM_SPADIC_CHA-0.5)*padWidth[suid2]);
      }
      title.Form("beamProfile_%i",usedSuId[suid]);
      name.Form("pad max distribution (beam profile) Sus%i",usedSuId[suid]);
      beamProfile[suid] = new TH1I(title, name, NUM_SPADIC_CHA, 0-0.5, NUM_SPADIC_CHA-0.5);
      beamProfile[suid]->SetXTitle("pad with maximum Q");
      beamProfile[suid]->SetYTitle("#");
      /*
        if (suid == usedSusibos-1) {
        title.Form("chamberAlignment_%i_%i",usedSuId[suid],usedSuId[0]);
        name.Form("chamber alignment Sus%i Sus%i",usedSuId[suid],usedSuId[0]);
        }
        else {
        title.Form("chamberAlignment_%i_%i",usedSuId[suid],usedSuId[suid+1]);
        name.Form("chamber alignment Sus%i Sus%i",usedSuId[suid],usedSuId[suid+1]);
        }
        chamberAlignment[suid] = new TProfile(title, name, NUM_SPADIC_CHA*padWidth[suid], 0-0.5, NUM_SPADIC_CHA*padWidth[suid]-0.5*padWidth[suid]);
        //printf("%f  %f  %f\n", NUM_SPADIC_CHA*padWidth[suid], 0-0.5, NUM_SPADIC_CHA*padWidth[suid]-0.5*padWidth[suid]);
        */
      title.Form("chaGain_%i",usedSuId[suid]);
      name.Form("channel gain pad max Sus%i",usedSuId[suid]);
      chaGain[suid] = new TH1F(title, name, NUM_SPADIC_CHA, 0-0.5, NUM_SPADIC_CHA-0.5);

      title.Form("amplitudeSpectra_%i",usedSuId[suid]);
      name.Form("amplitude spectra pad max Sus%i",usedSuId[suid]);
      amplitudeSpectra[suid] = new TH1I(title, name, SPADIC_ADC_MAX+50,-50,SPADIC_ADC_MAX);

      title.Form("hitTimeEl_%i",usedSuId[suid]);
      name.Form("hit time electron Sus%i",usedSuId[suid]);
      hitTimeEl[suid] = new TH1I(title, name, SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      hitTimeEl[suid]->SetXTitle("time-bin of maximum amplitude");
      hitTimeEl[suid]->SetYTitle("#");
      hitTimeEl[suid]->SetLineColor(2);

      title.Form("hitTimePi_%i",usedSuId[suid]);
      name.Form("hit time pion Sus%i",usedSuId[suid]);
      hitTimePi[suid] = new TH1I(title, name, SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      hitTimePi[suid]->SetXTitle("time-bin of maximum amplitude");
      hitTimePi[suid]->SetYTitle("#");

      title.Form("pulseShapeEl_%i",usedSuId[suid]);
      name.Form("pulse shape electron Sus%i",usedSuId[suid]);
      pulseShapeEl[suid] = new TProfile(title, name, SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      pulseShapeEl[suid]->SetLineColor(2);
      pulseShapeEl[suid]->GetYaxis()->SetRangeUser(-20,256);
      pulseShapeEl[suid]->SetXTitle("time-bin");
      pulseShapeEl[suid]->SetYTitle("average ADC-value");

      title.Form("pulseShapePi_%i",usedSuId[suid]);
      name.Form("pulse shape pion Sus%i",usedSuId[suid]);
      pulseShapePi[suid] = new TProfile(title, name, SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      pulseShapePi[suid]->GetYaxis()->SetRangeUser(-20,256);
      pulseShapePi[suid]->SetXTitle("time-bin");
      pulseShapePi[suid]->SetYTitle("average ADC-value");

      title.Form("PRF_%i",usedSuId[suid]);
      name.Form("PRF Sus%i",usedSuId[suid]);
      prf[suid] = new TH2D(title, name, 300, -1.5, 1.5, 100, 0, 1);
      prf[suid]->SetXTitle("relative pad charge");
      prf[suid]->SetYTitle("reco. hit position [pad units]");
      prf[suid]->SetContour(99); 

      title.Form("PRF_Profile%i",usedSuId[suid]);
      name.Form("PRF Profile Sus%i",usedSuId[suid]);
      prfProfile[suid] = new TProfile(title, name, 30 * padWidth[suid], -1.5 * padWidth[suid], 1.5 * padWidth[suid]);
      prfProfile[suid]->SetXTitle("relative pad charge");
      prfProfile[suid]->SetYTitle("reco. hit position [mm]");
    }
    TH1D *rawPulseA[NUM_SPADIC_CHA];
    TH1D *rawPulseB[NUM_SPADIC_CHA];
    TH1D *rawPulse[NUM_SPADIC_CHA];
    TH1D *baselinePulse[NUM_SPADIC_CHA];
    TH1D *baselinePulseA[NUM_SPADIC_CHA];
    TH1D *baselinePulseB[NUM_SPADIC_CHA];
    TH1D *baselineNoisePulse[NUM_SPADIC_CHA];
    TH1D *baselineNoisePulseA[NUM_SPADIC_CHA];
    TH1D *baselineNoisePulseB[NUM_SPADIC_CHA];
    //TString name;
    for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {     
      name.Form("Ch%02i",ch);
      rawPulseA[ch] = new TH1D("rawPulseA"+name,"rawPulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE); 
      rawPulseB[ch] = new TH1D("rawPulseB"+name,"rawPulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      rawPulse[ch] = new TH1D("rawPulse"+name,"rawPulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      baselinePulseA[ch] = new TH1D("baselinePulseA"+name,"baselinePulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      baselinePulseB[ch] = new TH1D("baselinePulseB"+name,"baselinePulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      baselinePulse[ch] = new TH1D("baselinePulse"+name,"baselinePulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      baselineNoisePulseA[ch] = new TH1D("baselineNoisePulseA"+name,"baselineNoisePulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      baselineNoisePulseB[ch] = new TH1D("baselineNoisePulseB"+name,"baselineNoisePulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
      baselineNoisePulse[ch] = new TH1D("baselineNoisePulse"+name,"baselineNoisePulse"+name,SPADIC_TRACE_SIZE,0,SPADIC_TRACE_SIZE);
    }
    while (fileSet.good()) {
      TString filename;
      fileSet >> filename;
      printf("  %s\n",filename.Data());
      if (filename.Sizeof() > 5) {
        TString ifname = filename;
        if(!ifname.EndsWith(".root")) ifname.Append(".root");
        //get_Pb_Ch_Cut(ifname);
 
        /*
          hCompensatedAll->GetYaxis()->SetRangeUser(-20,256);
          myc->cd();
          hCompensatedAll->Draw();
        */
        ifname = folder + filename;
        if(!ifname.EndsWith(".root")) ifname.Append(".root");
        TFile inputFile(ifname,"READ");
   
        TTree* theTree=0;
        TKey* kee=0;
        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 "<< ifname.Data() << endl;
          return;
        }
  


        // Thresholds
        Double_t minAmplitude[usedSusibos] =      {        15,        15,        15,        15,        15,        15,        15,       15   };
        Int_t HitTimeWindow[usedSusibos][2] =     {   {13,20},   {13,20},   {13,20},   {13,20},   {13,20},   {13,20},   {13,20},   {13,20}  };
        Int_t IntegrationWindow[usedSusibos][2] = {   {14,22},   {14,22},   {14,22},   {14,22},   {14,22},   {14,22},   {14,22},   {14,22}  };

        Double_t dxPosmm[usedSusibos] = {0.0};
 

        TCernOct11UnpackEvent* evnt = new TCernOct11UnpackEvent;
        TGo4EventElement* theBase=evnt;
        evnt->synchronizeWithTree(theTree, &theBase);
        //evnt->synchronizeWithTree(theTree, &evnt);

        // with go4 4.4.x:

        // evnt->synchronizeWithTree(theTree);
        // theTree->SetBranchAddress(theTree->GetListOfBranches()->First()->GetName(), &evnt);

        Int_t entries = theTree->GetEntries();
        allEvents += entries;
        Int_t nPion[usedSusibos][nMerge] = {{0}};
        Int_t nElectron[usedSusibos][nMerge] = {{0}};
        Int_t nGoodEvents[usedSusibos][nMerge] = {{0}};

        TMbsCrateEvent* fCrateInputEvent;
        TSpadicEvent* SpadicInputEvent;
        TSpadicData* theSpadic;
        TSpadicData* theSpadicA;
        TSpadicData* theSpadicB;
        //MFSpadicEvent* theMFSpadic = new MFSpadicEvent();
        Float_t Ch1(0), Ch2(0), Pb(0);
        cout << "   " << entries << " Entries" << endl;
        if (debug) 
          entries = 10000;
        GetPidCuts(filename);

        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"));
          //int sid=2; // example: get spadic of id 1
  
          isPion = false;
          isElectron = false;
          isMyon = 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);
          }
          Double_t maxAmplitudeCorr[usedSusibos] = {0.0};
          for (Int_t suid = 0; suid < usedSusibos; suid++){
   
            theSpadic=dynamic_cast<TSpadicData*>(SpadicInputEvent->getEventElement(usedSuId[suid]));
            Double_t pulse[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE] = {{0.0}};
            Bool_t signalChannel[NUM_SPADIC_CHA] = {false};
            Bool_t overFlow = theSpadic->fSpadicOverFlowEvent;//false;
            Bool_t underShoot = theSpadic->fSpadicUndershootsEvent;
            maxAmplitudeCorr[suid] = -40.0;
            dxPosmm[suid] = 0.0;

            for (Int_t imerge = 0; imerge < nMerge; imerge++) {
              msMerge = Bool_t(imerge);
              if (raw) {
                if((suid < 4) && (msMerge)) {
                  theSpadicA = dynamic_cast<TSpadicData*>(SpadicInputEvent->getEventElement(usedSuId[suid]));
                  if(theSpadicA==0) continue; // skip not configured ones
                  if(!theSpadicA->IsValid()) continue; // required to suppress empty spadic data!

                  theSpadicB = dynamic_cast<TSpadicData*>(SpadicInputEvent->getEventElement(topSuId[suid]));
                  if(theSpadicB==0) continue; // skip not configured ones
                  if(!theSpadicB->IsValid()) continue; // required to suppress empty spadic data!

                  if (theSpadicA->fSpadicOverFlowEvent || theSpadicB->fSpadicOverFlowEvent) overFlow = true;
                  if (theSpadicA->fSpadicUndershootsEvent || theSpadicB->fSpadicUndershootsEvent) underShoot = true;
                  //MergeData(pulse, theSpadicA, theSpadicB);

                  for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {
                    rawPulseA[ch]->Reset();
                    rawPulseB[ch]->Reset();
                    baselinePulseA[ch]->Reset();
                    baselinePulseB[ch]->Reset();
                    baselineNoisePulseA[ch]->Reset();
                    baselineNoisePulseB[ch]->Reset();
                    for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {       
                      rawPulseA[ch]->Fill(bin, (Double_t)theSpadicA->fSpadicPulse[ch][bin]);
                      rawPulseB[ch]->Fill(bin, (Double_t)theSpadicB->fSpadicPulse[ch][bin]);
                      //pulseA[ch][bin] = theSpadicA->fSpadicPulse[ch][bin];
                      //pulseB[ch][bin] = theSpadicB->fSpadicPulse[ch][bin];
                    }
                  }

                  BaselineCompensation(rawPulseA, baselinePulseA);
                  NoiseReduction(baselinePulseA, baselineNoisePulseA, signalChannel);

                  BaselineCompensation(rawPulseB, baselinePulseB);
                  NoiseReduction(baselinePulseB, baselineNoisePulseB, signalChannel);

                  for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {       
                    for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {
                      if (baselineNoisePulseA[ch]->GetBinContent(baselineNoisePulseA[ch]->FindBin(bin)) > 0)
                        pulse[ch][bin] = baselineNoisePulseA[ch]->GetBinContent(baselineNoisePulseA[ch]->FindBin(bin));//pulseA[ch][bin];
                      if (baselineNoisePulseB[7-ch]->GetBinContent(baselineNoisePulseB[7-ch]->FindBin(bin))) // channel remaping for MS top spadics
                        pulse[ch][bin] += baselineNoisePulseB[7-ch]->GetBinContent(baselineNoisePulseB[7-ch]->FindBin(bin));//pulseB[7-ch][bin];
                    }
                  }
                }
                else {
                  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)theSpadicA->fSpadicPulse[ch][bin]);// = theSpadic->fSpadicPulse[ch][bin];
                    }
                  }
                  BaselineCompensation(rawPulse, baselinePulse);
                  NoiseReduction(baselinePulse, baselineNoisePulse, signalChannel);
                  for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {       
                    for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {       
                      pulse[ch][bin] = baselineNoisePulse[ch]->GetBinContent(baselineNoisePulse[ch]->FindBin(bin));
                    }
                  }
                }
              } 
              else {
                if((suid < 4) && (msMerge)) {
                  theSpadicA = dynamic_cast<TSpadicData*>(SpadicInputEvent->getEventElement(usedSuId[suid]));
                  if(theSpadicA==0) continue; // skip not configured ones
                  if(!theSpadicA->IsValid()) continue; // required to suppress empty spadic data!

                  theSpadicB = dynamic_cast<TSpadicData*>(SpadicInputEvent->getEventElement(topSuId[suid]));
                  if(theSpadicB==0) continue; // skip not configured ones
                  if(!theSpadicB->IsValid()) continue; // required to suppress empty spadic data!

                  if (theSpadicA->fSpadicOverFlowEvent || theSpadicB->fSpadicOverFlowEvent) overFlow = true;
                  if (theSpadicA->fSpadicUndershootsEvent || theSpadicB->fSpadicUndershootsEvent) underShoot = true;
                  //MergeData(pulse, theSpadicA, theSpadicB);
                  Double_t pulseA[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE] = {{0.0}};
                  Double_t pulseB[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE] = {{0.0}};
        
                  for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {       
                    for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {       
                      pulseA[ch][bin] = theSpadicA->fSpadicCompensated[ch][bin];
                      pulseB[ch][bin] = theSpadicB->fSpadicCompensated[ch][bin];
                      if (theSpadicA->fSpadicSignalCh[ch] || theSpadicB->fSpadicSignalCh[ch])
                        signalChannel[ch] = true;;
                    }
                  }
                  BaselineCorrection(pulseA);
                  BaselineCorrection(pulseB);
                  for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {       
                    for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {
                      if (pulseA[ch][bin] > 0)
                        pulse[ch][bin] = pulseA[ch][bin];
                      if (pulseB[7-ch][bin] > 0) // channel remaping for MS top spadics
                        pulse[ch][bin] += pulseB[7-ch][bin];
                    }
                  }          
                }
                else{
                  //CopyData(pulse, theSpadic);
                  for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {       
                    for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {       
                      pulse[ch][bin] = theSpadic->fSpadicCompensated[ch][bin];
                      signalChannel[ch] = theSpadic->fSpadicSignalCh[ch];
                    }
                  }
                  BaselineCorrection(pulse);
                }
              }

              if(theSpadic==0) continue; // skip not configured ones
              if(!theSpadic->IsValid()) continue; // required to suppress empty spadic data!


              //if (overFlow || underShoot) continue; // skip fragmented signals

              // Find pad with maximum signal
              Double_t maxAmplitude = 0.0;
              Double_t maxAmplitudeCh[NUM_SPADIC_CHA] = {0.0};
              Int_t maxCh = -1;
              Int_t maxTB[NUM_SPADIC_CHA] = {-1};

              //cout << " - - - - - - - - " << endl;
              /*
                for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {
                //cout << endl;
                Float_t baseline(0);    
                for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {   
                pulse[ch][bin] = theSpadic->fSpadicCompensated[ch][bin];
                //cout << pulse[ch][bin] << " ";
                if (bin < noBaselineTB)
                baseline += pulse[ch][bin];
                }
                baseline /= Float_t(noBaselineTB);
                for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {   
                pulse[ch][bin] -= baseline;
                }
                }
              */

              for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {   
                for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {
                  if (theSpadic->fSpadicOverflows[ch][bin])
                    overFlow = true;
                  if (pulse[ch][bin] > maxAmplitudeCh[ch]){
                    maxAmplitudeCh[ch] = pulse[ch][bin];
                    maxTB[ch] = bin;
                  }
                }
                if(maxAmplitudeCh[ch] > maxAmplitude){
                  maxAmplitude = maxAmplitudeCh[ch];
                  maxCh = ch;
                }
              }
              if (!msMerge){
                if (maxCh != -1) {
                  maxAmplitudeCorr[suid] = pulse[maxCh][maxTB[maxCh]];
                  amplitudeSpectra[suid]->Fill(pulse[maxCh][maxTB[maxCh]]);
                  chaGainCalibSpectra[suid][maxCh]->Fill(pulse[maxCh][maxTB[maxCh]]);
                  //for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {     
                  //  chaGainCalibSpectra[suid][ch]->Fill(pulse[ch][maxTB[maxCh]]);
                  //}
                }
                else {
                  if (debug) printf("\n    <<<<< found no pad max\n");
                }
              }
              if (isElectron){
                hitTimeEl[suid]->Fill(maxTB[maxCh]);
              }
              if (isPion){
                hitTimePi[suid]->Fill(maxTB[maxCh]);
              }
              if (debug) {
                if (overFlow != theSpadic->fSpadicOverFlowEvent)
                  cout << "wrong overflow" << endl;
            
                if (maxCh != theSpadic->fSpadicHighestChannel)
                  cout << "online: " << theSpadic->fSpadicHighestChannel <<  "  maxA: " << maxAmplitudeCh[theSpadic->fSpadicHighestChannel]<< "  offline: " << maxCh <<  "  maxA: " << maxAmplitudeCh[maxCh] << endl;
              }
              if (maxCh < 1 || maxCh > 6)  continue; // skip fragmented signals

              // calc signal integral
              Double_t amplCluster(0), ampl3Pad(0), ampl1Pad(0);
              Double_t intensCluster(0), intens3Pad(0), intens1Pad(0), intensClusterTBW(0), intens3PadTBW(0), intens1PadTBW(0);
              Double_t intensCh[NUM_SPADIC_CHA] = {0.0};
              for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {
                for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {
                  if (pulse[ch][bin] > 0.0) {
                    intensCh[ch] +=  pulse[ch][bin];
                    if(signalChannel[ch]/*theSpadic->fSpadicSignalCh[ch]*/){
                      intensCluster += pulse[ch][bin];
                      if (bin > IntegrationWindow[suid][0] && bin < IntegrationWindow[suid][1])
                        intensClusterTBW += pulse[ch][bin];
                    }
                    if (ch == maxCh){
                      intens1Pad += pulse[ch][bin];
                      if (bin > IntegrationWindow[suid][0] && bin < IntegrationWindow[suid][1])
                        intens3PadTBW += pulse[ch][bin];
                    }
                    if (ch <= maxCh+1 && ch >= maxCh-1){
                      intens3Pad += pulse[ch][bin];
                      if (bin > IntegrationWindow[suid][0] && bin < IntegrationWindow[suid][1])
                        intens1PadTBW += pulse[ch][bin];
                    }
                  }
                }
                if (pulse[ch][maxTB[maxCh]] > 0.0) {
                  if(signalChannel[ch]/*theSpadic->fSpadicSignalCh[ch]*/)
                    amplCluster += pulse[ch][maxTB[maxCh]];
                  if (ch == maxCh)
                    ampl1Pad += pulse[ch][maxTB[maxCh]];
                  if (ch <= maxCh+1 && ch >= maxCh-1)
                    ampl3Pad += pulse[ch][maxTB[maxCh]];
                  /*
                    if(theSpadic->fSpadicSignalCh[ch])
                    amplCluster += maxAmplitudeCh[ch];
                    if (ch == maxCh)
                    ampl1Pad += maxAmplitudeCh[ch];
                    if (ch <= maxCh+1 && ch >= maxCh-1)
                    ampl3Pad += maxAmplitudeCh[ch];
                  */
                }               
              }
              // fill histos after cuts
              if (maxTB[maxCh] > HitTimeWindow[suid][0] && maxTB[maxCh] < HitTimeWindow[suid][1]){
                if (maxAmplitude > minAmplitude[suid]) {

                  if (imerge == 0) {
                    beamProfile[suid]->Fill(maxCh);
                    if (maxCh >= 1 && maxCh <= NUM_SPADIC_CHA - 1){
                      if (intensCh[maxCh-1] > 0.0 && intensCh[maxCh] > 0.0 && intensCh[maxCh+1] > 0.0) {
                        Double_t dxPos = 0.5 * log((intensCh[maxCh+1] / intensCh[maxCh-1])) / log((pow(intensCh[maxCh],2) / (intensCh[maxCh+1] * intensCh[maxCh-1])));
                        dxPosmm[suid] = dxPos * padWidth[suid] + maxCh * padWidth[suid];
                        Double_t chargePercent = intensCh[maxCh-1] / (intensCh[maxCh-1] + intensCh[maxCh] + intensCh[maxCh+1]);
                        prf[suid]->Fill(-dxPos-1,chargePercent);
                        prfProfile[suid]->Fill((-dxPos-1) * padWidth[suid],chargePercent);
                        chargePercent = intensCh[maxCh] / (intensCh[maxCh-1] + intensCh[maxCh] + intensCh[maxCh+1]);
                        prf[suid]->Fill(dxPos,chargePercent);
                        prfProfile[suid]->Fill(dxPos * padWidth[suid],chargePercent);
                        chargePercent = intensCh[maxCh+1] / (intensCh[maxCh-1] + intensCh[maxCh] + intensCh[maxCh+1]);
                        prf[suid]->Fill(-dxPos+1,chargePercent);
                        prfProfile[suid]->Fill((-dxPos+1) * padWidth[suid],chargePercent);
                      }
                    } 
                  }

                  for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {
                    if (isElectron){
                      pulseShapeEl[suid]->Fill(bin,pulse[maxCh][bin]);
                    }
                    if (isPion){
                      pulseShapePi[suid]->Fill(bin,pulse[maxCh][bin]);
                    }
                  }
                  for (Int_t ifilter = 0; ifilter < nFilter; ifilter++) { //{"wOwU","woOwU","woOwoU","wOwoU"}
                    if ((ifilter == 0 ) ||
                        (ifilter == 1 && !overFlow) ||
                        (ifilter == 2 && !overFlow && !underShoot) ||
                        (ifilter == 3 && !underShoot)
                        ) {
                      if (ifilter == 0 )
                        nGoodEvents[suid][imerge]++;
                      if (isElectron){
                        preCompSpectraEl[suid][0][imerge][ifilter]->Fill(intensCluster);
                        preCompSpectraEl[suid][1][imerge][ifilter]->Fill(intens3Pad);
                        preCompSpectraEl[suid][2][imerge][ifilter]->Fill(intens1Pad);
                        preCompSpectraEl[suid][3][imerge][ifilter]->Fill(intensClusterTBW);
                        preCompSpectraEl[suid][4][imerge][ifilter]->Fill(intens3PadTBW);
                        preCompSpectraEl[suid][5][imerge][ifilter]->Fill(intens1PadTBW);                  
                        preCompSpectraEl[suid][6][imerge][ifilter]->Fill(amplCluster);
                        preCompSpectraEl[suid][7][imerge][ifilter]->Fill(ampl3Pad);
                        preCompSpectraEl[suid][8][imerge][ifilter]->Fill(ampl1Pad);
                        if (ifilter == 0 )
                          nElectron[suid][imerge]++;
                      }
                      if (isPion){
                        preCompSpectraPi[suid][0][imerge][ifilter]->Fill(intensCluster);
                        preCompSpectraPi[suid][1][imerge][ifilter]->Fill(intens3Pad);
                        preCompSpectraPi[suid][2][imerge][ifilter]->Fill(intens1Pad);
                        preCompSpectraPi[suid][3][imerge][ifilter]->Fill(intensClusterTBW);
                        preCompSpectraPi[suid][4][imerge][ifilter]->Fill(intens3PadTBW);
                        preCompSpectraPi[suid][5][imerge][ifilter]->Fill(intens1PadTBW);
                        preCompSpectraPi[suid][6][imerge][ifilter]->Fill(amplCluster);
                        preCompSpectraPi[suid][7][imerge][ifilter]->Fill(ampl3Pad);
                        preCompSpectraPi[suid][8][imerge][ifilter]->Fill(ampl1Pad);
                        if (ifilter == 0 )
                          nPion[suid][imerge]++;
                      }
                    }
                  }
                } // if minimum amplitude threshold
              } // if hit time window
            } // for merge
          } // for suid

          Double_t distance[usedSusibos] = {0.0};
          //Distance2Track(dxPosmm, distance);

          Bool_t trackfoundMS = Distance2TrackMS(dxPosmm, distance/*, m, b*/);
          Bool_t trackfoundFFM = Distance2TrackFFM(dxPosmm, distance/*, m, b*/);

          Int_t MShitCount(0), FFMhitCount(0);
          for (Int_t suid = 0; suid < usedSusibos; suid++) {
            if (dxPosmm[suid] != 0.0){
              if (suid < 4)
                MShitCount++;
              else
                FFMhitCount++;
            }
          }
          for (Int_t suid = 0; suid < usedSusibos; suid++) {
            /*
              if (suid == usedSusibos-1) {
              chamberAlignment[suid]->Fill(dxPosmm[suid],dxPosmm[suid]-dxPosmm[suid]);
              //chamberAlignment2D[suid]->Fill(dxPosmm[suid],dxPosmm[suid]-dxPosmm[0]);
              }
              else {
              chamberAlignment[suid]->Fill(dxPosmm[suid],dxPosmm[suid]-dxPosmm[suid+1]);
              //chamberAlignment2D[suid]->Fill(dxPosmm[suid],dxPosmm[suid]-dxPosmm[suid+1]);
              }
            */
            for (Int_t suid2 = 0; suid2 < usedSusibos; suid2++) {
              if (dxPosmm[suid] > 0.0 && dxPosmm[suid2] > 0.0) {
                if (suid == suid2) {    
                  //if (debug) printf("%i   %.4f \n",suid,distance[suid]);
                  if(suid < 4){
                    if (trackfoundMS) {
                      chamberAlignment[suid][suid2]->Fill(dxPosmm[suid],distance[suid]);
                    }
                  }
                  else {
                    if (trackfoundFFM) {
                      chamberAlignment[suid][suid2]->Fill(dxPosmm[suid],distance[suid]);
                    }
                  }
                }
                else {
                  //if (suid < 4 && suid2 < 4){
                  //if (hitCount == 4){
                  chamberAlignment[suid][suid2]->Fill(dxPosmm[suid],dxPosmm[suid]-dxPosmm[suid2]);
                  //}
                  //}
                  //else {
                  //chamberAlignment[suid][suid2]->Fill(dxPosmm[suid],dxPosmm[suid]-dxPosmm[suid2]);
                  //}
                  if (debug) printf("  %.2f  %.2f  %.2f\n",dxPosmm[suid],dxPosmm[suid2],dxPosmm[suid]-dxPosmm[suid2]);
                }
              }
              amplitudeCorr[suid][suid2]->Fill(maxAmplitudeCorr[suid],maxAmplitudeCorr[suid2]);
            }
          }


        } // for tree entries

        for (Int_t suid = 0; suid < usedSusibos; suid++){
          printf("   Susibo%3i    %6i [%6i] electrons (%4.1f%%) [(%4.1f%%)]   %6i [%6i] pions (%4.1f%%) [(%4.1f%%)]   %6i [%6i] good events (%4.1f%%) [(%4.1f%%)]\n"
                 ,usedSuId[suid]
                 ,nElectron[suid][0],nElectron[suid][1],Float_t(nElectron[suid][0]*100)/entries,Float_t(nElectron[suid][1]*100)/entries
                 ,nPion[suid][0],nPion[suid][1],Float_t(nPion[suid][0]*100)/entries,Float_t(nPion[suid][1]*100)/entries
                 ,nElectron[suid][0]+nPion[suid][0],nElectron[suid][1]+nPion[suid][1],Float_t((nGoodEvents[suid][0])*100)/entries,Float_t((nGoodEvents[suid][1])*100)/entries
                 );
          for (Int_t imerge = 0; imerge < nMerge; imerge++) {
            allElectron[suid][imerge] += nElectron[suid][imerge];
            allPion[suid][imerge] += nPion[suid][imerge];
            allGoodEvent[suid][imerge] += nGoodEvents[suid][imerge];
          }
        }
        cout <<"*** DONE."<< endl;
 
        inputFile.Close();
      }
    } // while dataset is open
    
    Int_t c_x(1680*4), c_y(1050*4);
    if (debug){
      c_x /= 2;
      c_y /= 2;
    }

    TCanvas* cSpectra = new TCanvas("Spectra","Spectra",c_x,c_y);
    TCanvas* cRest = new TCanvas("Rest","Rest",c_x,c_x);
    TCanvas* cBeamMonitor = new TCanvas("BeamMonitor","Beam Monitor",c_x/8,c_y/8);
    TCanvas* cACorr = new TCanvas("cACorr","cACorr",c_x,c_y);
    TCanvas* cPCorr = new TCanvas("cPCorr","cPCorr",c_x,c_y);

    cSpectra->Divide(usedSusibos,recoSets);
    cRest->Divide(usedSusibos,12);
    cBeamMonitor->Divide(2,2);
    cACorr->Divide(usedSusibos,usedSusibos);
    cPCorr->Divide(usedSusibos,usedSusibos);

    cBeamMonitor->cd(1);
    Ch1_Pb->DrawCopy("colz");
    cBeamMonitor->cd(3);
    Ch2_Pb->DrawCopy("colz");
    cBeamMonitor->cd(2);
    noCh1_Pb->DrawCopy("colz");
    cBeamMonitor->cd(4);
    noCh2_Pb->DrawCopy("colz");  
    cBeamMonitor->Update();

    for (Int_t suid = 0; suid < usedSusibos; suid++) {
      for (Int_t suid2 = 0; suid2 < usedSusibos; suid2++) {
        cACorr->cd(suid2*usedSusibos+suid+1)->SetLogz(1);
        amplitudeCorr[suid][suid2]->DrawCopy("colz");
        outputFile.cd();
        amplitudeCorr[suid][suid2]->Write("", TObject::kOverwrite);

        cPCorr->cd(suid2*usedSusibos+suid+1);
        chamberAlignment[suid][suid2]->DrawCopy();
        chamberAlignment[suid][suid2]->Fit("alignmentFit","Q","0",0.25*NUM_SPADIC_CHA*padWidth[suid],0.75*NUM_SPADIC_CHA*padWidth[suid]);
        alignmentFit->SetLineColor(2);
        alignmentFit->DrawCopy("same");
        chamberAlignment[suid][suid2]->DrawCopy("same");
        chamberAlignment[suid][suid2]->Write("", TObject::kOverwrite);
      }
    }

    outputFile.cd();
    Ch1_Pb->Write("", TObject::kOverwrite);
    noCh1_Pb->Write("", TObject::kOverwrite);
    Ch2_Pb->Write("", TObject::kOverwrite);
    noCh2_Pb->Write("", TObject::kOverwrite);
    cBeamMonitor->Write("", TObject::kOverwrite);
    cACorr->Write("", TObject::kOverwrite);
    cPCorr->Write("", TObject::kOverwrite);

    Double_t MeanElectronSpectra(0), MeanPionSpectra(0);
    TLatex EP_Ratio;
    TString EP_temp;
    EP_Ratio.SetTextSize(0.04);
    TLegend* leg[usedSusibos][recoSets][nMerge][nFilter];
    Int_t x_max = 15000;
    Int_t rBin = 1;
    //Float_t x_pos(0.5), y_pos(0.75);
    for (Int_t imerge = 0; imerge < nMerge; imerge++) {
      for (Int_t ifilter = 0; ifilter < nFilter; ifilter++) {
        for (Int_t suid = 0; suid < usedSusibos; suid++){  
          for (Int_t rs = 0; rs < recoSets; rs++) {  
            if (rs > 5 && rs < recoSets){
              x_max = 1000;
              //rBin = 2;
            }
            else
              {
                x_max = 15000;
                //rBin = 17;
              }
            leg[suid][rs][imerge][ifilter] = new TLegend(0.5,0.6,0.85,0.85);
            leg[suid][rs][imerge][ifilter]->SetLineColor(0);
            leg[suid][rs][imerge][ifilter]->SetLineStyle(0);
            leg[suid][rs][imerge][ifilter]->SetFillColor(0);
            leg[suid][rs][imerge][ifilter]->SetFillStyle(0);
            leg[suid][rs][imerge][ifilter]->SetTextSize(0.035);

            preCompSpectraEl[suid][rs][imerge][ifilter]->Scale(1./Float_t(preCompSpectraEl[suid][rs][imerge][ifilter]->GetEntries()));
            preCompSpectraEl[suid][rs][imerge][ifilter]->Rebin(rBin);
            if (preCompSpectraEl[suid][rs][imerge][ifilter]->GetEntries() > 0)
              MeanElectronSpectra = preCompSpectraEl[suid][rs][imerge][ifilter]->GetMean(1);

            preCompSpectraPi[suid][rs][imerge][ifilter]->Scale(1./Float_t(preCompSpectraPi[suid][rs][imerge][ifilter]->GetEntries()));
            preCompSpectraPi[suid][rs][imerge][ifilter]->Rebin(rBin);
            if (preCompSpectraPi[suid][rs][imerge][ifilter]->GetEntries() > 0)
              MeanPionSpectra = preCompSpectraPi[suid][rs][imerge][ifilter]->GetMean(1);

            preCompSpectraEl[suid][rs][imerge][ifilter]->GetXaxis()->SetRangeUser(1,x_max);
            preCompSpectraPi[suid][rs][imerge][ifilter]->GetXaxis()->SetRangeUser(1,x_max);
            EP_temp.Form("#LTe^{-}#GT = %.3f",MeanElectronSpectra);
            leg[suid][rs][imerge][ifilter]->AddEntry(preCompSpectraEl[suid][rs][imerge][ifilter],EP_temp,"l");
            EP_temp.Form("#LT#pi^{-}#GT = %.3f",MeanPionSpectra);
            leg[suid][rs][imerge][ifilter]->AddEntry(preCompSpectraPi[suid][rs][imerge][ifilter],EP_temp,"l");
            if (MeanPionSpectra != 0)
              EP_temp.Form("#LTe^{-}#GT/#LT#pi^{-}#GT = %.3f",MeanElectronSpectra/MeanPionSpectra);
            else
              EP_temp.Form("#LTe^{-}#GT/#LT#pi^{-}#GT = %.3f",0.0);
            leg[suid][rs][imerge][ifilter]->AddEntry((TObject*)0,EP_temp,"");
            if (debug) printf("        %s  %s\n",recoSetName[rs].Data(),EP_temp.Data());
            if (preCompSpectraEl[suid][rs][imerge][ifilter]->GetEntries() > 0 &&  preCompSpectraPi[suid][rs][imerge][ifilter]->GetEntries() > 0) {
              cSpectra->cd(rs * usedSusibos + suid + 1);
              preCompSpectraPi[suid][rs][imerge][ifilter]->DrawCopy();
              preCompSpectraEl[suid][rs][imerge][ifilter]->DrawCopy("same");
              leg[suid][rs][imerge][ifilter]->Draw("same");
            }
          }
        }
        cSpectra->Update();
        cSpectra->SaveAs("data2012/Analysis/Pics/Spectra_" + TString(nameFileSet(0,nameFileSet.Sizeof()-5)) + "_" + mergeName[imerge] + "_" + filterName[ifilter] + /*picEx*/".png");
        cSpectra->Write(TString("Spectra_" + TString(nameFileSet(0,nameFileSet.Sizeof()-5)) + "_" + mergeName[imerge] + "_" + filterName[ifilter]), TObject::kOverwrite);
      }
    }
    TLegend *lPRF[usedSusibos];
    for (Int_t suid = 0; suid < usedSusibos; suid++){  
      lPRF[suid] = new TLegend(0.5,0.7,0.95,0.95);
      cRest->cd(1+suid);
      pulseShapeEl[suid]->DrawCopy();
      pulseShapePi[suid]->DrawCopy("same");
      cRest->cd(9+suid);
      hitTimePi[suid]->DrawCopy();
      hitTimeEl[suid]->DrawCopy("same");
      cRest->cd(17+suid)->SetLogy(1);
      amplitudeSpectra[suid]->DrawCopy();
      cRest->cd(25+suid);
      prf[suid]->DrawCopy("colz");
      cRest->cd(33+suid);
      prfProfile[suid]->GetYaxis()->SetRangeUser(0,1);
      prfProfile[suid]->DrawCopy();
      prfProfile[suid]->Fit("sigmaFit","RQ0");
      TString sl;
      sl.Form("#sigma = %f",sigmaFit->GetParameter(2));
      lPRF[suid]->AddEntry(sigmaFit,sl,"lep");
      if (suid >= 0 && suid < 2) {
        prfProfile[suid]->Fit("mathiesonMS336","RQ0");
        sl.Form("K_{3} = %f",mathiesonMS336->GetParameter(0));
      }
      if (suid >= 2 && suid < 4) {
        prfProfile[suid]->Fit("mathiesonMS444","RQ0");
        sl.Form("K_{3} = %f",mathiesonMS444->GetParameter(0));
      }      
      lPRF[suid]->AddEntry(prfProfile[suid],sl,"lep");
      sigmaFit->SetLineColor(2);
      sigmaFit->DrawCopy("same");
      if (suid < 4)
        prfProfile[suid]->DrawCopy("same");
      prfProfile[suid]->DrawCopy("same");
      lPRF[suid]->Draw("same");
      cRest->cd(41+suid);
      Int_t suid2min;     
      if (suid < 4) {
        suid2min = 0;
      }
      else {
        suid2min = 4;
      }
      Int_t suid2max = suid2min + 4;
      for (Int_t suid2 = suid2min; suid2 < suid2max; suid2++){
        preCompSpectraPi[suid2][1][0][0]->SetLineColor(suid2+1);
        preCompSpectraPi[suid2][1][0][0]->GetXaxis()->SetRangeUser(0,6000);
        if (suid == suid2)
          preCompSpectraPi[suid2][1][0][0]->SetLineStyle(2);
        if (suid2 == 0)
          preCompSpectraPi[suid2][1][0][0]->DrawCopy();
        else 
          preCompSpectraPi[suid2][1][0][0]->DrawCopy("same");
        preCompSpectraPi[suid2][1][0][0]->SetLineColor(1);
        preCompSpectraPi[suid2][1][0][0]->SetLineStyle(1);
        preCompSpectraPi[suid2][1][0][0]->GetXaxis()->SetRangeUser(1,x_max);
      }
      cRest->cd(49+suid)->SetLogy(1);
      for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {
        chaGainCalibSpectra[suid][ch]->GetXaxis()->SetRangeUser(-50,255);
        if (ch == 0)
          chaGainCalibSpectra[suid][ch]->DrawCopy();
        else
          chaGainCalibSpectra[suid][ch]->DrawCopy("same");

        chaGainCalibSpectra[suid][ch]->Fit("gausFit","RQ0");
        gausFit->DrawCopy("same");
        chaGain[suid]->Fill(ch, gausFit->GetParameter(1));
      }
      cRest->cd(57+suid);
      chaGain[suid]->DrawCopy();
  
      cRest->cd(65+suid);
      beamProfile[suid]->DrawCopy();

      cRest->cd(73+suid);
      /*
        chamberAlignment[suid]->DrawCopy();
        chamberAlignment[suid]->Fit("alignmentFit","RQ0");
        alignmentFit->DrawCopy("same");
      */
      cRest->Update();      
    

      //-------------------------------------------------------------------------------  write to file
      outputFile.cd();

      for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {
        chaGainCalibSpectra[suid][ch]->Write("", TObject::kOverwrite);
      }
   
      for (Int_t rs = 0; rs < recoSets; rs++) {
        for (Int_t imerge = 0; imerge < nMerge; imerge++) {
          for (Int_t ifilter = 0; ifilter < nFilter; ifilter++) {
            preCompSpectraEl[suid][rs][imerge][ifilter]->Write("", TObject::kOverwrite);
            preCompSpectraPi[suid][rs][imerge][ifilter]->Write("", TObject::kOverwrite);   
          }
        }
      }
    

      beamProfile[suid]->Write("", TObject::kOverwrite);

      chaGain[suid]->Write("", TObject::kOverwrite);

      amplitudeSpectra[suid]->Write("", TObject::kOverwrite);

      pulseShapeEl[suid]->Write("", TObject::kOverwrite);
      pulseShapePi[suid]->Write("", TObject::kOverwrite);

      hitTimeEl[suid]->Write("", TObject::kOverwrite);
      hitTimePi[suid]->Write("", TObject::kOverwrite);

      prf[suid]->Write("", TObject::kOverwrite);
      prfProfile[suid]->Write("", TObject::kOverwrite);
      cRest->Write("", TObject::kOverwrite);
      
      printf("         Susibo%3i    %6i [%6i] electrons (%4.1f%%) [(%4.1f%%)]   %6i [%6i] pions (%4.1f%%) [(%4.1f%%)]   %6i [%6i] good events (%4.1f%%) [(%4.1f%%)]\n"
             ,usedSuId[suid]
             ,allElectron[suid][0],allElectron[suid][1],Float_t(allElectron[suid][0]*100)/allEvents,Float_t(allElectron[suid][1]*100)/allEvents
             ,allPion[suid][0],allPion[suid][1],Float_t(allPion[suid][0]*100)/allEvents,Float_t(allPion[suid][1]*100)/allEvents
             ,allPion[suid][0]+allElectron[suid][0],allPion[suid][1]+allElectron[suid][1],Float_t((allGoodEvent[suid][0])*100)/allEvents,Float_t((allGoodEvent[suid][1])*100)/allEvents
             );

     
      
    } // for suid
    
      //------------------------------------------------------------------------------- save figures
    cBeamMonitor->SaveAs("data2012/Analysis/Pics/BeamMonitor_" + TString(nameFileSet(0,nameFileSet.Sizeof()-5)) + picEx);
  
    cRest->SaveAs("data2012/Analysis/Pics/Rest_" + TString(nameFileSet(0,nameFileSet.Sizeof()-5)) + /*picEx*/".png");

    cACorr->SaveAs("data2012/Analysis/Pics/ACorr_" + TString(nameFileSet(0,nameFileSet.Sizeof()-5)) + /*picEx*/".png");

    cPCorr->SaveAs("data2012/Analysis/Pics/PR/PCorr_" + TString(nameFileSet(0,nameFileSet.Sizeof()-5)) + /*picEx*/".png");

    //cChaGainCal->SaveAs("data2012/Analysis/Pics/Gain_" + TString(nameFileSet(0,nameFileSet.Sizeof()-5)) + picEx);
    outputFile.Close();
    /*
      for (Int_t suid = 0; suid < usedSusibos; suid++){
      for (Int_t rs = 0; rs < recoSets; rs++) {
      for (Int_t imerge = 0; imerge < nMerge; imerge++) {
      for (Int_t ifilter = 0; ifilter < nFilter; ifilter++) {
      delete preCompSpectraEl[suid][rs][imerge][ifilter];
      delete preCompSpectraPi[suid][rs][imerge][ifilter];
      delete leg[suid][rs][imerge][ifilter];
      }
      }
      }
    
      delete hitTimeEl[suid];
      delete hitTimePi[suid];
      delete pulseShapeEl[suid];
      delete pulseShapePi[suid];
      delete prf[suid];
      }
    
      delete cBeamMonitor;
      delete cSpectra;
      delete cRest;
      delete Ch1_Pb;
      delete noCh1_Pb;
      delete Ch2_Pb;
      delete noCh2_Pb;
    */
  } // if dataset is open
  else {
    printf("no input file list in %s",nameFileSet.Data());
    return;
  }
}