//////////////////////////////////////////////////////////////
//                                                          //
// This macro produces all final results and QA plots for   //
// the DEUTERON nuclei analysis.                            //
//                                                          //
//                                                          //
//                                                          //
//////////////////////////////////////////////////////////////

#include "TFile.h"
#include "TList.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TH3D.h"
#include "THnSparse.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TFractionFitter.h"
#include "TColor.h"



using namespace std;

//
// global variables
//

Char_t * fileNameMcEnhancedNew = "MC/McEnhancedNew.root";
Char_t * fileNameMcEnhancedOld = "MC/McEnhancedOld2.root";
Char_t * fileNameMcUnEnhanced  = "MC/McUnEnhanced.root";
Char_t * fileNameData          = "data/dataFinal.root";

//
// definition of functions
//
void       MakeEfficiency();
TH1D    *  ExtractEfficiency(TList * listMC, Int_t sign = -1, Int_t centrality = 0, Bool_t isTOF = kTRUE);
TF1     *  MakeEtaRapidityGenCorrection(Bool_t drawQA = kFALSE);
//
void       MakeFinalSpectra();
TH1D    *  MakeRawSpectraTPC(TList * list, Int_t sign = -1, Int_t centralityBin = 0);
TH1D    *  MakeRawSpectra(TList * list, Int_t sign = -1, Int_t centralityBin = 0);
void       NormalizeSpectrum(TH1D * spectrum, Float_t dy, Float_t numberOfEvents);
void       PlotSpectraComparison();
void       CombinationAndSystematics();
TH1D    *  GetSpectrumWithSyst(TH1D * histStatError);
//
void       MakeQAplots();
TCanvas *  PlotTpcQA(Char_t * fileName = "", Bool_t isMC = kTRUE,  Bool_t isMCold = kFALSE);
void       PlotSpectraQA(TList * list, Int_t particleType = 0, Int_t sign = -1);
//
void       MakeMaterialCorrection();
Float_t    GetMaterialCorrection(Float_t ptBin, Int_t centralityBin, Bool_t pureMC = kFALSE);



//
// implementation
//
//_______________________________________________________________________
void CombinationAndSystematics() {
  //
  // produce the final spectra for the paper
  //
  Int_t kMaxCentrality = 5;
  //
  TFile fileOut("output/spectraDeuteronCombined.root", "RECREATE");
  fileOut.Close();
  //
  TFile * fileIn = TFile::Open("output/spectraDeuteron.root");
  //
  TCanvas * canvDeuteron = new TCanvas("canvDeuteron","deuteron spectra");
  //
  for(Int_t iCentr = 0; iCentr < kMaxCentrality; iCentr++) {
    //
    TH1D * tpc = (TH1D *) fileIn->Get(Form("cent%i_sign%s_TPC",iCentr,"Pos"));
    TH1D * tof = (TH1D *) fileIn->Get(Form("cent%i_sign%s",iCentr,"Pos"));
    //
    TH1D * deuteron = (TH1D *) tof->Clone();
    deuteron->SetNameTitle(Form("deuteron_centrality%i",iCentr), Form("deuteron_centrality%i",iCentr));
    deuteron->GetYaxis()->SetTitle("d#it{N}/(d#it{y}d#it{p_{T}})");
    deuteron->GetXaxis()->SetTitle("p_{T} (GeV/#it{c})");
    //
    // merge the two
    //
    deuteron->SetBinContent(4, tpc->GetBinContent(4));  deuteron->SetBinError(4, tpc->GetBinError(4));
    deuteron->SetBinContent(5, tpc->GetBinContent(5));  deuteron->SetBinError(5, tpc->GetBinError(4));
    //
    // get systematics
    //
    TH1D * deuteronSyst = GetSpectrumWithSyst(deuteron);
    //
    // drawing and saving
    //
    if (iCentr == 0) {
      deuteronSyst->DrawCopy("E2");
      deuteron->DrawCopy("SAME");
    } else {
      deuteronSyst->DrawCopy("E2,SAME");
      deuteron->DrawCopy("SAME");
    }
    //

    //
    TFile fileIn2("output/fitsBW.root");
    TF1 * fit = (TF1 * ) fileIn2.Get(Form("fit%i_signNeg", iCentr));
    if (fit) {
      fit->Draw("SAME");
    }
    fileIn2.Close();
    //
    TFile fileOut("output/spectraDeuteronCombined.root", "UPDATE");
    deuteron->Write();
    deuteronSyst->Write();
    if (fit) fit->Write();
    fileOut.Close();

  }



}

//_______________________________________________________________________
TH1D * GetSpectrumWithSyst(TH1D * histStatError) {
  //
  // add systematic error to the points
  //
  TH1D * histSyst = (TH1D *) histStatError->Clone();
  histSyst->SetNameTitle(Form("%s_SYST", histStatError->GetName()), 
			 Form("%s_SYST", histStatError->GetName()));
  histSyst->SetFillColor(histStatError->GetMarkerColor());
  histSyst->SetFillStyle(0);
  histSyst->SetDrawOption("E2");
  //
  for(Int_t i=0; i <  histSyst->GetXaxis()->GetNbins(); i++) { // begin loop over pt-bins
      Float_t syst = 0.;
      Float_t pt = histSyst->GetXaxis()->GetBinCenter(i);
      //
      // (1.) tracking and matching normally 4%, but we add 2% for the pt-correction
      //
      syst += 0.06*0.06;
      //
      // (2.) PID -- difference between bin counting and fit in TOF
      //
      if (pt > 1.4) syst += 0.05*0.05;
      //
      // (3.) material knock-out -- 20% of correction
      //
      Float_t corr = 1.5*TMath::Exp(1.27259 - 2.527*pt);
      syst += (0.2*corr)*(0.2*corr);
      //
      // (4.) absorption in material (shifted proton exponential, difference between Eulogio's and geant production)
      // ==> for anti-deuterons correspondingly more. ==> see plot by Natasha (backup slide APW)
      //
      syst += 0.03*0.03;
      //
      // (5.) matching efficiency
      //
      //if (pt > 1.) syst += 0.05*0.05;
      //
      // put it to the spectrum
      //
      Float_t  sum = TMath::Sqrt(syst)*histSyst->GetBinContent(i);
      histSyst->SetBinError(i,sum);
    } // end loop over pt-bins

  //
  return histSyst;


}


//_______________________________________________________________________
void PlotSpectraComparison() {
  //
  // plot TPC & TOF spectra for comparison
  //
  Int_t kMaxCentrality = 5;
  //
  //
  TFile * fileIn = TFile::Open("output/spectraDeuteron.root");
  //
  TCanvas * canvRatio = new TCanvas("canvRatio","canvRatio");
  //
  for(Int_t iCentr = 0; iCentr < kMaxCentrality; iCentr++) {
    //
    TCanvas * canvComparison = new TCanvas(Form("canvComparison_%i",iCentr),Form("spectra in centrality %i", iCentr));
    //
    TH1D * tpc = (TH1D *) fileIn->Get(Form("cent%i_sign%s_TPC",iCentr,"Pos"));
    TH1D * tof = (TH1D *) fileIn->Get(Form("cent%i_sign%s",iCentr,"Pos"));
    //
    TH1D * tpcNeg = (TH1D *) fileIn->Get(Form("cent%i_sign%s_TPC",iCentr,"Neg"));
    TH1D * tofNeg = (TH1D *) fileIn->Get(Form("cent%i_sign%s",iCentr,"Neg"));
    tpcNeg->SetMarkerStyle(21);
    tofNeg->SetMarkerStyle(25);
    //
    tof->DrawCopy();
    tpc->DrawCopy("SAME");
    tpcNeg->DrawCopy("SAME");
    tofNeg->DrawCopy("SAME");
    //
    TFile fileIn2("output/fitsBW.root");
    TF1 * fit = (TF1 * ) fileIn2.Get(Form("fit%i_signNeg", iCentr));
    fit->Draw("SAME");
    fileIn2.Close();
    //
    //TCanvas * canvRatio = new TCanvas(Form("canvRatio_%i",iCentr),Form("Ratio in centrality %i", iCentr));
    canvRatio->cd();
    tofNeg->Divide(tof);
    if (iCentr == 0) {
      tofNeg->DrawCopy();
    } else {
      tofNeg->DrawCopy("SAME");
    }
  }


}


//_______________________________________________________________________
void MakeFinalSpectra() {
  //
  // make final spectra for all centrality bins
  //
  Int_t kMaxCentrality = 5;
  TFile fileOut("output/spectraDeuteron.root", "RECREATE");
  fileOut.Close();
  //
  // get efficiencies and data
  //
  TFile * inFileEff = TFile::Open("output/efficiencies.root");
  TH1D  * efficiencyTrackingNegNew = (TH1D*) inFileEff->Get("efficiencyTrackingNegNew");
  TH1D  * efficiencyTrackingPosNew = (TH1D*) inFileEff->Get("efficiencyTrackingPosNew");
  TH1D  * efficiencyTofPosNew      = (TH1D*) inFileEff->Get("efficiencyTofPosNew");
  TH1D  * efficiencyTofNegNew      = (TH1D*) inFileEff->Get("efficiencyTofNegNew");
  //
  TFile inFileData(fileNameData);
  TList * list = (TList * ) inFileData.Get("akalweit_Nuclei");
  //
  // positive particles
  //
  TCanvas * canvAllPos = new TCanvas("canvAllPos","deuteron spectra");
  for(Int_t iCentr = 0; iCentr < kMaxCentrality; iCentr++) {
    //
    TH1D * rawSpecTpc = MakeRawSpectraTPC(list, +1, iCentr);
    TH1D * rawSpecTof = MakeRawSpectra(list, +1, iCentr);
    rawSpecTpc->Divide(efficiencyTrackingPosNew);
    rawSpecTof->Divide(efficiencyTofPosNew);
    //
    canvAllPos->cd();
    if (iCentr == 0) {
      rawSpecTpc->DrawCopy("EP");
      rawSpecTof->DrawCopy("EPSAME");
    } else {
      rawSpecTpc->DrawCopy("EPSAME");
      rawSpecTof->DrawCopy("EPSAME");
    }
    //
    // apply material correction
    //
    TFile matFile("output/materialCorrection.root");
    TH1D * materialCorrection = (TH1D *) matFile.Get(Form("matCorr_%i",iCentr));
    for(Int_t iBin = 0; iBin < rawSpecTpc->GetXaxis()->GetNbins(); iBin++) {
      Double_t pT = rawSpecTpc->GetXaxis()->GetBinCenter(iBin);
      if (pT > 0.4 && pT < 2.) {
	Float_t primFactor = 1. - materialCorrection->GetBinContent(iBin);
	rawSpecTof->SetBinContent(iBin, rawSpecTof->GetBinContent(iBin)*primFactor);
	rawSpecTof->SetBinError(iBin, rawSpecTof->GetBinError(iBin)*primFactor);
	rawSpecTpc->SetBinContent(iBin, rawSpecTpc->GetBinContent(iBin)*primFactor);
	rawSpecTpc->SetBinError(iBin, rawSpecTpc->GetBinError(iBin)*primFactor);
      }
    }
    //
    matFile.Close();
    //
    TFile fileOutput("output/spectraDeuteron.root", "UPDATE");
    rawSpecTpc->Write();
    rawSpecTof->Write();
    fileOutput.Close();
  }
  canvAllPos->Print("QAplots/posSpectra.pdf");
  //
  // negative particles
  //
  TCanvas * canvAllNeg = new TCanvas("canvAllNeg","anti-deuteron spectra");
  for(Int_t iCentr = 0; iCentr < kMaxCentrality; iCentr++) {
    //
    TH1D * rawSpecTpc = MakeRawSpectraTPC(list, -1, iCentr);
    TH1D * rawSpecTof = MakeRawSpectra(list, -1, iCentr);
    rawSpecTpc->Divide(efficiencyTrackingNegNew);
    rawSpecTof->Divide(efficiencyTofNegNew);
    //
    canvAllPos->cd();
    if (iCentr == 0) {
      rawSpecTpc->DrawCopy("EP");
      rawSpecTof->DrawCopy("EPSAME");
    } else {
      rawSpecTpc->DrawCopy("EPSAME");
      rawSpecTof->DrawCopy("EPSAME");
    }
    TFile fileOutput("output/spectraDeuteron.root", "UPDATE");
    rawSpecTpc->Write();
    rawSpecTof->Write();
    fileOutput.Close();
  }
  canvAllNeg->Print("QAplots/negSpectra.pdf");
  


}


//_______________________________________________________________________
TH1D * MakeRawSpectraTPC(TList * list, Int_t sign, Int_t centralityBin) {
  //
  // Make raw spectra for given sign and centrality bin
  //
  // (0.) assumed particle: 0. deuteron, 1. triton, 2. He-3
  // (1.) multiplicity or centrality -- number of accepted ESD tracks per events (deprecated), but now classes from 1 to 10, 0: Min. Bias
  // (2.) pT
  // (3.) sign
  // (4.) rapidity --> filled 4xW
  // (5.)  pull TPC dEx --> filled 4x
  // (6.) has valid TOF pid signal
  // (7.) nsigma TOF --> filled 4x XXXXXXXXX no mass*mass
  // (8..) dca_xy
  // (9.) CODE -- only MC 0-generated, 1-true rec. primaries, 2-misident prim, 3-second weak, 4-second material, 5-misident sec
  //
  TH1D * histCentr = (TH1D *)  list->FindObject("fHistCentrality");
  //
  THnSparseF * hist = (THnSparseF *) list->FindObject("fHistRealTracks");
  hist->GetAxis(0)->SetRangeUser(0,0);                     // select deuterons
  //
  // centrality 0: 0-5% / 1: 5-10% / 2: 10-20% / 3: 20-30% / 4: 30-40% / 5: 40-50% / 6: 50-60% / 7: 60-70% / 8: 70-80%
  //
  Int_t centrality = 0;
  if (centralityBin == 0) centrality = 0; //  0-10%
  if (centralityBin == 1) centrality = 2; // 10-20%
  if (centralityBin == 2) centrality = 3; // 20-40%
  if (centralityBin == 3) centrality = 5; // 40-60%
  if (centralityBin == 4) centrality = 7; // 60-80%
  //
  hist->GetAxis(1)->SetRangeUser(centrality,centrality);   // select centrality bin
  Float_t norm = histCentr->GetBinContent(centrality+2);
  if (centralityBin != 1) {
    hist->GetAxis(1)->SetRangeUser(centrality,centrality+1);   // select centrality bin
    norm += histCentr->GetBinContent(centrality+3);
  }
  //
  // apply cuts -- Has to be consistent with efficiency
  //
  hist->GetAxis(4)->SetRangeUser(-0.49, 0.49);             // select rapidity range
  hist->GetAxis(3)->SetRangeUser(sign,sign);               // select anti-deuterons
  hist->GetAxis(8)->SetRangeUser(-0.5,0.5);  // DCA-cut / TODO: do proper unfolding
  hist->GetAxis(5)->SetRangeUser(-2.5,2.5);  // TPC-PID cut  
  //
  TH1D * spec    = hist->Projection(2);
  spec->SetNameTitle("spec","SPEC");
  TH1D * rawSpec = hist->Projection(2);
  rawSpec->Reset();
  for(Int_t iBin = 0; iBin < spec->GetXaxis()->GetNbins(); iBin++) {
    Double_t pT = spec->GetXaxis()->GetBinCenter(iBin);
    if (pT > 0.6 && pT < 1.2) {
      rawSpec->SetBinContent(iBin, spec->GetBinContent(iBin));
      rawSpec->SetBinError(iBin, spec->GetBinError(iBin));
    }
  }
  //
  //
  NormalizeSpectrum(rawSpec,1,norm);
  //
  if (sign < 0) rawSpec->SetNameTitle(Form("cent%i_sign%s_TPC",centralityBin,"Neg"),
				      Form("cent%i_sign%s_TPC",centralityBin,"Neg"));
  if (sign > 0) rawSpec->SetNameTitle(Form("cent%i_sign%s_TPC",centralityBin,"Pos"),
				      Form("cent%i_sign%s_TPC",centralityBin,"Pos"));
  //
  Int_t colorList[9] = {600+4, 600+3, 600+2, 600+1, 600, 600-4, 600-7, 600-9, 600-10};
  rawSpec->SetMarkerStyle(20);
  rawSpec->SetMarkerSize(1.4);
  rawSpec->SetMarkerColor(colorList[centrality]);
  rawSpec->SetLineColor(colorList[centrality]);
  //
  hist->GetAxis(5)->SetRangeUser(-1000.,1000);  // TPC-PID cut  
  hist->GetAxis(6)->SetRangeUser(-1000.,1000);  // TOF-PID cut  
  hist->GetAxis(7)->SetRangeUser(-1000.,1000);  // TOF-PID cut    
  hist->GetAxis(8)->SetRangeUser(-1000.,1000);  // TOF-PID cut    

  //
  return rawSpec;

}



//_______________________________________________________________________
TH1D * MakeRawSpectra(TList * list, Int_t sign, Int_t centralityBin) {
  //
  // Make raw spectra for given sign and centrality bin
  //
  //
  // (0.) assumed particle: 0. deuteron, 1. triton, 2. He-3
  // (1.) multiplicity or centrality -- number of accepted ESD tracks per events (deprecated), but now classes from 1 to 10, 0: Min. Bias
  // (2.) pT
  // (3.) sign
  // (4.) rapidity --> filled 4xW
  // (5.)  pull TPC dEx --> filled 4x
  // (6.) has valid TOF pid signal
  // (7.) nsigma TOF --> filled 4x XXXXXXXXX no mass*mass
  // (8..) dca_xy
  // (9.) CODE -- only MC 0-generated, 1-true rec. primaries, 2-misident prim, 3-second weak, 4-second material, 5-misident sec
  //
  TH1D * histCentr = (TH1D *)  list->FindObject("fHistCentrality");
  //
  THnSparseF * hist = (THnSparseF *) list->FindObject("fHistRealTracks");
  hist->GetAxis(0)->SetRangeUser(0,0);                     // select deuterons
  //
  // centrality 0: 0-5% / 1: 5-10% / 2: 10-20% / 3: 20-30% / 4: 30-40% / 5: 40-50% / 6: 50-60% / 7: 60-70% / 8: 70-80%
  //
  Int_t centrality = 0;
  if (centralityBin == 0) centrality = 0; //  0-10%
  if (centralityBin == 1) centrality = 2; // 10-20%
  if (centralityBin == 2) centrality = 3; // 20-40%
  if (centralityBin == 3) centrality = 5; // 40-60%
  if (centralityBin == 4) centrality = 7; // 60-80%
  //
  hist->GetAxis(1)->SetRangeUser(centrality,centrality);   // select centrality bin
  Float_t norm = histCentr->GetBinContent(centrality+2);
  if (centralityBin != 1) {
    hist->GetAxis(1)->SetRangeUser(centrality,centrality+1);   // select centrality bin
    norm += histCentr->GetBinContent(centrality+3);
  }
  //
  //
  //
  hist->GetAxis(4)->SetRangeUser(-0.49, 0.49);             // select rapidity range
  hist->GetAxis(3)->SetRangeUser(sign,sign);               // select anti-deuterons
  //
  hist->GetAxis(8)->SetRangeUser(-0.5,0.5);  // DCA-cut / TODO: do proper unfolding
  //
  hist->GetAxis(5)->SetRangeUser(-2.5,2.5);  // TPC-PID cut  
  //hist->GetAxis(7)->SetRangeUser(-2.,2.5);
  //
  // PID in the TOF
  //
  hist->GetAxis(6)->SetRangeUser(1,1);        // TOF-PID cut: require hasTOF
  //
  TH1D * rawSpec = hist->Projection(2);
  rawSpec->Reset();
  if (sign < 0) rawSpec->SetNameTitle(Form("cent%i_sign%s",centralityBin,"Neg"),
				      Form("cent%i_sign%s",centralityBin,"Neg"));
  if (sign > 0) rawSpec->SetNameTitle(Form("cent%i_sign%s",centralityBin,"Pos"),
				      Form("cent%i_sign%s",centralityBin,"Pos"));
  //
  // --> simple PID: hist->GetAxis(7)->SetRangeUser(-0.8.,0.8);  // TOF-PID cut
  //TF1 fitFunc("fitFunc","gaus(0) + [3] + [4]*x",-1.5,1.5);
  TF1 fTOFsignal("fTOFsignal", "(x <= ([3] + [1])) * [0] * TMath::Gaus(x, [1], [2]) + (x > ([3] + [1])) * [0] * TMath::Gaus([3] + [1], [1], [2]) * TMath::Exp(-([3]) * (x - [3] - [1]) / ([2] * [2]))    +    [4] + [5]*x   +   [6]*TMath::Exp(-[7]*x) ", -2.2,2.2); //0:norm, 1:mean, 2:sigma, 3:tail
  //
  //
  TH2D * m2fit = (TH2D *) hist->Projection(7,2);
  m2fit->SetNameTitle(Form("m2fits_%i_%i",sign,centralityBin),
		      Form("m**2 fits for sign %i and centrality %i",sign,centralityBin));
  hist->GetAxis(5)->SetRangeUser(-1000.,1000);  // TPC-PID cut  
  hist->GetAxis(6)->SetRangeUser(-1000.,1000);  // TOF-PID cut  
  hist->GetAxis(7)->SetRangeUser(-1000.,1000);  // TOF-PID cut    
  hist->GetAxis(8)->SetRangeUser(-1000.,1000);  // TOF-PID cut    
  //
  // rebin because of statistics
  //
  m2fit->RebinY(2);
  //
  TCanvas * canvM2fits = new TCanvas(Form("canvM2fits_%i_%i",sign,centralityBin),
				     Form("m**2 fits for sign %i and centrality %i",sign,centralityBin));
  canvM2fits->Print(Form("QA/M2fits_%i_%i.pdf(",sign,centralityBin));
  //
  Double_t oldParStartValue[8] = {664.78 , 0.0178744 , 0.212441 , 0.229438 , 171.889, 0. , 0.0203292};
  //
  for(Int_t iBin = 1; iBin < m2fit->GetXaxis()->GetNbins(); iBin++) {
    //
    Float_t ptBin = m2fit->GetXaxis()->GetBinCenter(iBin);
    if (ptBin < 0.6) continue;
    //
    TH1D * m2distr = m2fit->ProjectionY(Form("m2_%f",ptBin),iBin,iBin);
    //    m2distr->RebinX(2);
    //
    m2distr->SetMarkerColor(kBlue);
    m2distr->SetMarkerStyle(20);
    m2distr->SetMarkerSize(1.4);
    //
    // initialize signal with background parameters
    //
    fTOFsignal.SetParameters(oldParStartValue);
    //
    // improve fit convergence
    //
    fTOFsignal.SetParLimits(1,-0.2,0.2);
    fTOFsignal.SetParLimits(2,0.05,0.8);
    fTOFsignal.SetParLimits(3,0.,1.5);
    //
    fTOFsignal.SetParameter(4, m2distr->GetBinContent(m2distr->GetXaxis()->FindBin(-1.)));
    fTOFsignal.SetParLimits(4, 0.,  m2distr->GetBinContent(m2distr->GetXaxis()->FindBin(0.)));
    //
    fTOFsignal.SetParLimits(7,0.,10.);
    fTOFsignal.SetParLimits(6,0., 5.);
    //
    if (centralityBin < 3 && ptBin < 2.4) {
      fTOFsignal.SetRange(-1.3,2.2);
    } else{
      fTOFsignal.SetRange(-2.2,2.2);
    }
    //
    Int_t fitStatus = m2distr->Fit(&fTOFsignal, "QNR");
    cout << fTOFsignal.GetParameter(3) << endl;
    //
    for(Int_t ival =0; ival < 7; ival++) oldParStartValue[ival] = fTOFsignal.GetParameter(ival);
    //
    canvM2fits->cd();
    m2distr->GetYaxis()->SetRangeUser(0.1, m2distr->GetBinContent(m2distr->GetXaxis()->FindBin(0.))*1.5);
    m2distr->DrawCopy("EP");
    fTOFsignal.DrawCopy("same");
    //
    //
    // bin counting at low pt and subtract the background at high pt
    //
    Float_t yield =1;//  m2distr->Integral(m2distr->GetXaxis()->FindBin(-0.5), m2distr->GetXaxis()->FindBin(+0.5)) - 
      (m2distr->Integral(m2distr->GetXaxis()->FindBin(-1.), m2distr->GetXaxis()->FindBin(-0.5)) + m2distr->Integral(m2distr->GetXaxis()->FindBin(0.5), m2distr->GetXaxis()->FindBin(1.)));
    //
    //
    //
    cout << "pT: " << ptBin 
	 << " , " << fTOFsignal.GetParameter(0)
	 << " , " << fTOFsignal.GetParameter(1)
	 << " , " << fTOFsignal.GetParameter(2)
	 << " , " << fTOFsignal.GetParameter(3)
	 << " , " << fTOFsignal.GetParameter(4)
	 << " , " << fTOFsignal.GetParameter(5)
	 << " , " << fTOFsignal.GetParameter(6)
	 << " , " << fTOFsignal.GetParameter(7) << endl;
    //
    //
    fTOFsignal.SetParameter(4,0.);
    fTOFsignal.SetParameter(5,0.);
    fTOFsignal.SetParameter(6,0.);
    //
    //if (ptBin > 2.3) 
    yield = fTOFsignal.Integral(-1.,1.5)/m2distr->GetBinWidth(10);
    //     yield = m2distr->Integral(m2distr->GetXaxis()->FindBin(-1.), m2distr->GetXaxis()->FindBin(+1.));
    //
    Float_t yieldErr = TMath::Sqrt(yield);
    //    if (fTOFsignal.GetParameter(0) != 0 ) yieldErr = yield*fTOFsignal.GetParError(0)/fTOFsignal.GetParameter(0);
    //
    Float_t maxPt = 4.2;
    if (centralityBin > 2) maxPt = 3.2;
    if (ptBin > 0.5 && ptBin < maxPt && yield > 0) {
      rawSpec->SetBinContent(iBin, yield);
      rawSpec->SetBinError(iBin, yieldErr);
    }
    //
    //
    canvM2fits->Print(Form("QA/M2fits_%i_%i.pdf",sign,centralityBin));
  }
  canvM2fits->Print(Form("QA/M2fits_%i_%i.pdf)",sign,centralityBin));
  //
  rawSpec->Sumw2();
  //
  NormalizeSpectrum(rawSpec,1,norm);
  //
  Int_t colorList[9] = {600+4, 600+3, 600+2, 600+1, 600, 600-4, 600-7, 600-9, 600-10};
  rawSpec->SetMarkerStyle(24);
  rawSpec->SetMarkerSize(1.4);
  rawSpec->SetMarkerColor(colorList[centrality]);
  rawSpec->SetLineColor(colorList[centrality]);
  //
  return rawSpec;

}

//_______________________________________________________________________
void NormalizeSpectrum(TH1D * spectrum, Float_t dy, Float_t numberOfEvents) {
  //
  // make (1 / Nev)*  dN/(dy*dpt)
  //
  spectrum->Scale(1./dy);
  spectrum->Scale(1./numberOfEvents);
  //    
  Int_t nBins = spectrum->GetNbinsX(); 
  for(Int_t i=0; i < nBins+1; i++) {
    Double_t Content =spectrum->GetBinContent(i);
    Double_t error =spectrum->GetBinError(i);
    if (spectrum->GetBinWidth(i) == 0) continue;
    spectrum->SetBinContent(i, Content/spectrum->GetBinWidth(i));
    spectrum->SetBinError(i, error/spectrum->GetBinWidth(i));
  }

}



//_______________________________________________________________________
void MakeEfficiency() {
  //
  // make efficiencies for TPC and TOF part.
  //
  TFile inFileMCold(fileNameMcEnhancedOld);
  TList * listMCold = (TList * ) inFileMCold.Get("akalweit_Nuclei");
  //
  //
  // TPC and tracking as well as matching efficiencies
  //
  //
  TH1D * efficiencyTrackingNegOld = ExtractEfficiency(listMCold, -1, 0, kFALSE);
  TH1D * efficiencyTrackingPosOld = ExtractEfficiency(listMCold, +1, 0, kFALSE);
  efficiencyTrackingNegOld->SetNameTitle("efficiencyTrackingNegOld","efficiencyTrackingNegOld");
  efficiencyTrackingPosOld->SetNameTitle("efficiencyTrackingPosOld","efficiencyTrackingPosOld");
  //
  TH1D * efficiencyTofNegOld = ExtractEfficiency(listMCold, -1, 0, kTRUE);
  TH1D * efficiencyTofPosOld = ExtractEfficiency(listMCold, +1, 0, kTRUE);
  efficiencyTofNegOld->SetNameTitle("efficiencyTofNegOld","efficiencyTofNegOld");
  efficiencyTofPosOld->SetNameTitle("efficiencyTofPosOld","efficiencyTofPosOld");
  //
  TF1 * etaGenCorrection = MakeEtaRapidityGenCorrection(kFALSE);
  efficiencyTrackingPosOld->Divide(etaGenCorrection);
  efficiencyTrackingNegOld->Divide(etaGenCorrection);
  efficiencyTofPosOld->Divide(etaGenCorrection);
  efficiencyTofNegOld->Divide(etaGenCorrection);
  //
  // new simulation
  //
  TFile inFileMCnew(fileNameMcEnhancedNew);
  TList * listMCnew = (TList * ) inFileMCnew.Get("akalweit_Nuclei");
  //
  TH1D * efficiencyTrackingNegNew = ExtractEfficiency(listMCnew, -1, 0, kFALSE);
  TH1D * efficiencyTrackingPosNew = ExtractEfficiency(listMCnew, +1, 0, kFALSE);
  efficiencyTrackingNegNew->SetNameTitle("efficiencyTrackingNegNew","efficiencyTrackingNegNew");
  efficiencyTrackingPosNew->SetNameTitle("efficiencyTrackingPosNew","efficiencyTrackingPosNew");
  efficiencyTrackingNegNew->SetMarkerStyle(21);
  efficiencyTrackingPosNew->SetMarkerStyle(21);
  //
  TH1D * efficiencyTofNegNew = ExtractEfficiency(listMCnew, -1, 0, kTRUE);
  TH1D * efficiencyTofPosNew = ExtractEfficiency(listMCnew, +1, 0, kTRUE);
  efficiencyTofNegNew->SetNameTitle("efficiencyTofNegNew","efficiencyTofNegNew");
  efficiencyTofPosNew->SetNameTitle("efficiencyTofPosNew","efficiencyTofPosNew");
  efficiencyTofNegNew->SetMarkerStyle(24);
  efficiencyTofPosNew->SetMarkerStyle(24);  
  //
  TCanvas * canvEfficiencyTracking = new TCanvas("canvEfficiencyTracking","Efficiencies Tracking");
  efficiencyTrackingNegOld->GetYaxis()->SetRangeUser(0.,1.3);
  efficiencyTrackingNegOld->DrawCopy("E");
  efficiencyTrackingPosOld->DrawCopy("ESame");
  //
  efficiencyTrackingNegNew->DrawCopy("ESame");
  efficiencyTrackingPosNew->DrawCopy("ESame");
  //
  TCanvas * canvEfficiencyTof = new TCanvas("canvEfficiencyTof","Efficiencies Tracking");
  efficiencyTofNegNew->GetYaxis()->SetRangeUser(0.,1.3);
  efficiencyTofNegNew->DrawCopy("E");
  efficiencyTofPosNew->DrawCopy("ESame");
  efficiencyTofNegOld->DrawCopy("ESAME");
  efficiencyTofPosOld->DrawCopy("ESame");
  //
  canvEfficiencyTracking->Print("QAplots/efficiencyTracking.pdf");
  canvEfficiencyTof->Print("QAplots/efficiencyTof.pdf");
  //
  // store output
  //
  TFile * outFile = new TFile("output/efficiencies.root","RECREATE");
  efficiencyTrackingNegNew->Write();
  efficiencyTrackingPosNew->Write();
  efficiencyTofNegNew->Write();
  efficiencyTofPosNew->Write();
  outFile->Close();
  //
  //
  TCanvas * canvEfficiencyMatching = new TCanvas("canvEfficiencyMatching","matching efficiency");
  efficiencyTofNegNew->Divide(efficiencyTrackingNegNew);
  efficiencyTofPosNew->Divide(efficiencyTrackingPosNew);
  efficiencyTofPosNew->DrawCopy();
  efficiencyTofNegNew->DrawCopy("SAME");
  //
  TCanvas * canvEfficiencyRatio = new TCanvas("canvEfficiencyRatio","Efficiencies ratio");
  efficiencyTrackingPosNew->Divide(efficiencyTrackingPosOld);
  efficiencyTrackingPosNew->DrawCopy();
  efficiencyTrackingNegNew->Divide(efficiencyTrackingNegOld);
  efficiencyTrackingNegNew->DrawCopy("SAME");

  

}



//_______________________________________________________________________
TH1D * ExtractEfficiency(TList * listMC, Int_t sign, Int_t centrality, Bool_t isTOF) {
  //
  // Extract the efficiency
  //
  //
  // (0.) assumed particle: 0. deuteron, 1. triton, 2. He-3
  // (1.) multiplicity or centrality -- number of accepted ESD tracks per events (deprecated), but now classes from 1 to 10, 0: Min. Bias
  // (2.) pT
  // (3.) sign
  // (4.) rapidity --> filled 4x
  // (5.)  pull TPC dEx --> filled 4x
  // (6.) has valid TOF pid signal
  // (7.) nsigma TOF --> filled 4x XXXXXXXXX no mass*mass
  // (8..) dca_xy
  // (9.) CODE -- only MC 0-generated, 1-true rec. primaries, 2-misident prim, 3-second weak, 4-second material, 5-misident sec
  //
  THnSparseF * hist = (THnSparseF *) listMC->FindObject("fHistMCparticles");
  //
  // common selections for generated and reconstructed particles
  // 
  hist->GetAxis(0)->SetRangeUser(0,0);                     // select deuterons
  hist->GetAxis(4)->SetRangeUser(-0.49, 0.49);             // select rapidity range
  hist->GetAxis(3)->SetRangeUser(sign,sign);               // select anti-deuterons
  //assume eff. indepedent of centrality  ---> hist->GetAxis(1)->SetRangeUser(centrality,centrality);   // select centrality bin
  //
  hist->GetAxis(9)->SetRangeUser(0,0);
  TH1D * generated = hist->Projection(2);
  generated->SetNameTitle(Form("generated_cent%i_sign%i",centrality,sign),
			  Form("generated_cent%i_sign%i",centrality,sign));
  //
  // further cuts on reconstructed particles
  //
  hist->GetAxis(8)->SetRangeUser(-0.5,0.5);  // DCA-cut / TODO: do proper unfolding
  //
  hist->GetAxis(5)->SetRangeUser(-4.5,4.5);  // TPC-PID cut  
  //
  if (isTOF) {
    //    hist->GetAxis(5)->SetRangeUser(-2.,2.);  // TPC-PID cut  
    hist->GetAxis(6)->SetRangeUser(1.,1.);  // TOF-PID cut  
    hist->GetAxis(7)->SetRangeUser(-1.2,1.2);  // TOF-PID cut  
  }
  //
  hist->GetAxis(9)->SetRangeUser(1,1);
  TH1D * rawSpec = hist->Projection(2);
  rawSpec->Sumw2();
  rawSpec->SetNameTitle(Form("eff_cent%i_sign%i",centrality,sign),
			Form("eff_cent%i_sign%i",centrality,sign));
  //
  rawSpec->Divide(generated);
  //
  if (sign > 0) {
    rawSpec->SetMarkerColor(kRed);
    rawSpec->SetLineColor(kRed);
  } else{
    rawSpec->SetMarkerColor(kBlue);
    rawSpec->SetLineColor(kBlue);
  }
  //
  // reset ranges
  //
  hist->GetAxis(5)->SetRangeUser(-1000.,1000);  // TPC-PID cut  
  hist->GetAxis(6)->SetRangeUser(-1000.,1000);  // TOF-PID cut  
  hist->GetAxis(7)->SetRangeUser(-1000.,1000);  // TOF-PID cut    
  hist->GetAxis(8)->SetRangeUser(-1000.,1000);  // TOF-PID cut    
  //
  //
  //
  return rawSpec;

}

//__________________________________________________________________
TF1 * MakeEtaRapidityGenCorrection(Bool_t drawQA) {
  //
  // the enhanced sample LHC11b9_1 has a stupid cut on the generated level
  //
  TF1 * funcDeut = new TF1("funcDeut","TMath::ASinH(TMath::Sqrt(1 + (1.876*1.876)/(x*x))*TMath::SinH(0.5))",0.5,5);
  TF1 * etaCut = new TF1("etaCut","0.9",0.5,5);
  //
  funcDeut->SetLineColor(kBlue);
  //
  if (drawQA) {
    TCanvas * canv1 = new TCanvas("canv1","original functions");
    funcDeut->DrawCopy();
    funcDeut->GetHistogram()->GetYaxis()->SetRange(0,1.8);
    etaCut->DrawCopy("same");
  }
  //
  TF1 * funcCorr = new TF1("funcCorr","funcDeut/etaCut*(x - 1.1 < 0) + 1*(x - 1.1 >= 0)",0.5,2);
  if (drawQA) {
    TCanvas * canv2 = new TCanvas("canv2","correction function");
    funcCorr->Draw();
  }
  //
  return funcCorr;

}


//__________________________________________________________________
void MakeQAplots() {
  //
  // make the QA plots
  //
  TCanvas * canvData          =  PlotTpcQA(fileNameData,          kFALSE);
  TCanvas * canvMcEnhancedNew =  PlotTpcQA(fileNameMcEnhancedNew, kTRUE);
  TCanvas * canvMcUnEnhanced  =  PlotTpcQA(fileNameMcUnEnhanced,  kTRUE, kTRUE);
  TCanvas * canvMcEnhancedOld =  PlotTpcQA(fileNameMcEnhancedOld, kTRUE, kTRUE);
  //
  TCanvas * canvQA = new TCanvas("canvQA","canvQA");
  canvQA->Print("QAplots/QAplots.pdf(");
  //
  // add all QA plots to file
  //
  canvData->Print("QAplots/QAplots.pdf");
  canvMcEnhancedNew->Print("QAplots/QAplots.pdf");
  canvMcUnEnhanced->Print("QAplots/QAplots.pdf");
  canvMcEnhancedOld->Print("QAplots/QAplots.pdf");
  //
  canvQA->Print("QAplots/QAplots.pdf)");
  //
  // spectra QA
  //
  TFile inFile(fileNameData);
  TList * list = (TList * ) inFile.Get("akalweit_Nuclei");
  PlotSpectraQA(list);


}


//__________________________________________________________________
TCanvas * PlotTpcQA(Char_t * fileName, Bool_t isMC, Bool_t isMCold) {
  //
  // plot QA of TPC-PID
  //
  TFile inFileMC(fileName);
  TList * list = (TList * ) inFileMC.Get("akalweit_Nuclei");
  //
  TF1 * fBBlines[6];
  const Double_t masses[6] = {0.000511, 0.138, 0.4936, 0.938, 2*0.938, 3*0.938};
  //
  // set BB parameters
  //
  Double_t parData[5] = {1.45802, 27.4992, 4.00313e-15, 2.48485, 8.31768};
  Double_t parMC[5] = {1.17329, 27.4992, 4.00313e-15, 2.1204316, 4.1373729}; 
  Double_t parMCold[5] = {1.17329, 27.4992, 4.00313e-15, 2.35563, 9.47569}; // OLD FOR LHC11b9_1 !!
  //
  for(Int_t ii = 0; ii < 6; ii++) {
    fBBlines[ii] = new TF1(Form("fBB_%i",ii),"AliExternalTrackParam::BetheBlochAleph(x/([0]),[1],[2],[3], [4], [5])",0.3,10);
    
    fBBlines[ii]->SetParameters(masses[ii], parData[0], parData[1], parData[2], parData[3], parData[4]);
    if (isMC) {
      fBBlines[ii]->SetParameters(masses[ii], parMC[0], parMC[1], parMC[2], parMC[3], parMC[4]);
      if (isMCold) fBBlines[ii]->SetParameters(masses[ii], parMCold[0], parMCold[1], parMCold[2], parMCold[3], parMCold[4]);
    }
    fBBlines[ii]->SetLineColor(kRed);
    fBBlines[ii]->SetLineWidth(2);
  }
  //
  TCanvas * canvDeDx = new TCanvas(Form("canvDeDxTpcQA_%s",fileName),"control histogram for dE/dx");
  canvDeDx->Divide(1,2);
  canvDeDx->cd(1);
  TH3D * histTPC = (TH3D *) list->FindObject("fHistPidQA");
  histTPC->GetZaxis()->SetRangeUser(-1,-1);
  TH2D * histNegTPC = (TH2D *) histTPC->Project3D("yx");
  histNegTPC->SetNameTitle("histNegTPC", Form("histNegTPC_%s", fileName));
  gPad->SetLogx();
  histNegTPC->DrawCopy("colZ");
  for(Int_t ii = 0; ii < 6; ii++) fBBlines[ii]->DrawCopy("same");
  //
  canvDeDx->cd(2);
  histTPC->GetZaxis()->SetRangeUser(+1,+1);
  TH2D * histPosTPC = (TH2D*) histTPC->Project3D("yx");
  histPosTPC->SetNameTitle("histPosTPC", Form("histPosTPC_%s", fileName));
  gPad->SetLogx();
  histPosTPC->DrawCopy("colZ");
  //  for(Int_t ii = 0; ii < 6; ii++) fBBlines[ii]->DrawCopy("same");
  fBBlines[4]->DrawCopy("same"); // Draw deuteron line
  //
  // He-3 goes extra
  //
  TF1 * funcHel3 = new TF1("funcHel3","4*AliExternalTrackParam::BetheBlochAleph(2*x/(0.938*3),1.74962,27.4992,4.00313e-15,2.42485,8.31768)",0.2,6.);
  if (isMC) funcHel3 = new TF1("funcHel3MC","4*AliExternalTrackParam::BetheBlochAleph(2*x/(0.938*3), 1.17329, 27.4992, 4.00313e-15, 2.1204316, 4.1373729)",0.2,6.);
  //
  fBBlines[5]->DrawCopy("same"); // Draw deuteron line
  funcHel3->DrawCopy("same");

  return canvDeDx;

}



//_______________________________________________________________________
void PlotSpectraQA(TList * list, Int_t particleType, Int_t sign) {
  //
  // Make some basic QA plots for TOF related PID
  //
  THnSparseF * hist = (THnSparseF *) list->FindObject("fHistRealTracks");
  hist->GetAxis(4)->SetRangeUser(-0.49, 0.49);             // select Rapidity range
  hist->GetAxis(0)->SetRangeUser(particleType,particleType); // select deuterons
  hist->GetAxis(3)->SetRangeUser(sign,sign); // select deuterons
  //
  TCanvas * canvSpectraQA1 = new TCanvas("canvSpectraQA1","QA for nsigma TPC");
  TH2F * histTpcNsigma = (TH2F *) hist->Projection(5,2);
  histTpcNsigma->SetNameTitle("histTpcNsigma","histTpcNsigma");
  histTpcNsigma->DrawCopy("colZ");
  //
  TCanvas * canvSpectraQA2 = new TCanvas("canvSpectraQA2","QA for nsigma TOF");
  hist->GetAxis(6)->SetRangeUser(1,1);   // require proper match to TOF
  hist->GetAxis(5)->SetRangeUser(-4.,4); // pre-select tracks with TPC-PID
  TH2F * histTofNsigma = (TH2F *) hist->Projection(7,2);
  histTofNsigma->RebinY(4);
  histTofNsigma->SetNameTitle("histTofNsigma","histTofNsigma");
  histTofNsigma->DrawCopy("colZ");

}



//_______________________________________________________________________
void MakeMaterialCorrection(){
  //
  // compare fitted values with MC pt-dependence
  //
  TFile outFile("output/materialCorrection.root","RECREATE");
  outFile.Close();
  //
  Int_t kMaxCentrality = 5;
  TFile * fileIn = TFile::Open("output/efficiencies.root");
  //
  //
  for(Int_t iCentr = 0; iCentr < kMaxCentrality; iCentr++) {
    //
    TH1D * matCorrMC = (TH1D *) fileIn->Get("efficiencyTrackingNegNew");
    matCorrMC->Reset();
    matCorrMC->SetNameTitle(Form("matCorrMC_%i",iCentr), Form("matCorrMC_i",iCentr));
    //
    TH1D * matCorr = (TH1D *) fileIn->Get("efficiencyTrackingNegNew");
    matCorr->SetLineColor(kRed);
    matCorr->Reset();
    matCorr->SetNameTitle(Form("matCorr_%i",iCentr), Form("matCorr_i",iCentr));    
    //
    for(Int_t iBin = 0; iBin < matCorrMC->GetXaxis()->GetNbins(); iBin++) {
      Double_t pT = matCorrMC->GetXaxis()->GetBinCenter(iBin);
      if (pT > 0.6 && pT < 2.8) {
	Float_t corr = GetMaterialCorrection(pT,iCentr,kTRUE);
	matCorrMC->SetBinContent(iBin, corr);
	if (pT < 2.) {
	  Float_t corrReal = GetMaterialCorrection(pT,iCentr,kFALSE);
	  matCorr->SetBinContent(iBin, corrReal);
	}
      }
    }
    //
    //
    TCanvas * canvMaterial = new TCanvas("canvMaterial","canvMaterial");
    canvMaterial->cd();
    if (iCentr == 0) {
      matCorr->DrawCopy();
      matCorrMC->DrawCopy("SAME");
    } else {
      matCorrMC->DrawCopy("SAME");
    }
    TFile * outputFile = new TFile("output/materialCorrection.root","UPDATE");
    matCorr->Write();
    outputFile->Close();
    delete outputFile;
  }


}



//_______________________________________________________________________
Float_t GetMaterialCorrection(Float_t ptBin, Int_t centralityBin, Bool_t pureMC) {
  //
  // subtract the material contamination from positive particles
  //
  //
  Float_t dcaCutInAnalysis = 0.5;
  //
  TFile * inFileData = TFile::Open("data/dataFinal.root");
  TList * listData = (TList *) inFileData->Get("akalweit_Nuclei");
  THnSparse * histData = (THnSparse *) listData->FindObject("fHistRealTracks");
  //
  TFile * inFileMC = TFile::Open("MC/McCombined.root");
  TList * listMC = (TList *) inFileMC->Get("akalweit_Nuclei");
  THnSparse * histMC = (THnSparse *) listMC->FindObject("fHistMCparticles");
  //
  //
  // (1.) select deuterons
  //
  histData->GetAxis(0)->SetRangeUser(0,0);
  histMC->GetAxis(0)->SetRangeUser(0,0);
  histData->GetAxis(3)->SetRangeUser(+1,+1); // sign
  histMC->GetAxis(3)->SetRangeUser(+1,+1);   // sign
  //
  histData->GetAxis(4)->SetRangeUser(-0.49, 0.49);             // select rapidity range
  histMC->GetAxis(4)->SetRangeUser(-0.49, 0.49);             // select rapidity range
  //
  //  hist->GetAxis(8)->SetRangeUser(-1.,1.);  // DCA-cut / TODO: do proper unfolding
  histData->GetAxis(5)->SetRangeUser(-2.5,2.5);  // TPC-PID cut  
  histData->GetAxis(6)->SetRangeUser(1.,1.);  // TOF-PID cut  
  histData->GetAxis(7)->SetRangeUser(-0.5,0.5);  // TOF-PID cut  
  //
  // (2.) select pt-range
  //
  histData->GetAxis(2)->SetRangeUser(ptBin,ptBin);
  histMC->GetAxis(2)->SetRangeUser(ptBin,ptBin);
  //
  // (4.a) get different MC templates
  //
  histMC->GetAxis(9)->SetRangeUser(1,1);
  TH1D * prim = (TH1D*) histMC->Projection(8);
  prim->SetNameTitle("prim","prim");
  prim->SetLineColor(kRed);
  //
  histMC->GetAxis(9)->SetRangeUser(4,4);
  //
  // centrality can oonly be selected after the primary yield
  //
  Int_t centrality = 0;
  if (centralityBin == 0) centrality = 0; //  0-10%
  if (centralityBin == 1) centrality = 2; // 10-20%
  if (centralityBin == 2) centrality = 3; // 20-40%
  if (centralityBin == 3) centrality = 5; // 40-60%
  if (centralityBin == 4) centrality = 7; // 60-80%
  //
  histData->GetAxis(1)->SetRangeUser(centrality,centrality);       // select centrality bin
  histMC->GetAxis(1)->SetRangeUser(centrality,centrality);       // select centrality bin
  if (centralityBin != 1) {
    histData->GetAxis(1)->SetRangeUser(centrality,centrality+1);   // select centrality bin
    histMC->GetAxis(1)->SetRangeUser(centrality,centrality+1);     // select centrality bin
  }
  //
  TH1D * material = (TH1D*) histMC->Projection(8);
  histMC->GetAxis(2)->SetRangeUser(ptBin,ptBin);
  material->SetNameTitle("material","material");  
  material->SetLineColor(kGreen);
  //
  TH1D * dcaDistr = (TH1D*) histData->Projection(8);
  //
  // if only MC, we just return the values
  //
  if (pureMC) {
    Int_t fractionBinLow = dcaDistr->FindBin(-dcaCutInAnalysis);
    Int_t fractionBinUp  = dcaDistr->FindBin(+dcaCutInAnalysis);
    Float_t materialFraction = 0.;
    prim->Add(material); // has to be scaled with number of events somehow
    materialFraction = material->Integral(fractionBinLow,fractionBinUp)/prim->Integral(fractionBinLow,fractionBinUp);
    //
    delete listData;
    delete listMC;
    inFileData->Close();
    inFileMC->Close();
    delete inFileData;
    delete inFileMC;
    //
    return materialFraction;
  }
  //
  // (6.) do the fit
  //
  Float_t dcaUp = 0.5;//0.0182 + 0.035/TMath::Power(ptBin,1.01);
  Float_t dcaLow  = -1.*dcaUp;
  Int_t binLow = dcaDistr->GetXaxis()->FindBin(dcaLow);
  Int_t binUp  = dcaDistr->GetXaxis()->FindBin(dcaUp);
  //
  TObjArray *mcShapes = new TObjArray(3);        // MC histograms are put in this array
  mcShapes->Add(prim);
  mcShapes->Add(material);
  //
  TFractionFitter* fit = new TFractionFitter(dcaDistr, mcShapes); // initialise
  fit->SetRangeX(binLow,binUp);
  fit->Constrain(1,0.,1.);
  fit->Constrain(2,0.,1.);
  //
  //
  //
  Int_t status = fit->Fit();
  if (status==0) fit->GetPlot()->Draw("same");
  //
  // (7.) Rescaling of histograms and final plots
  //
  TH1D * result = 0x0;
  Double_t yieldPrim, yieldSec, yieldMat, error;
  if (status == 0) {
    fit->GetResult(0,yieldPrim,error);
    fit->GetResult(1,yieldMat,error);
    result = (TH1D*) fit->GetPlot();
    result->SetLineWidth(3);
  }
  //
  prim->Scale(1./prim->Integral(binLow,binUp));
  if (status == 0) prim->Scale(yieldPrim*result->Integral(binLow,binUp));
  //
  material->Scale(1./material->Integral(binLow,binUp));
  if (status == 0) material->Scale(yieldMat*result->Integral(binLow,binUp));
  //
  NormalizeSpectrum(dcaDistr,1,1);
  NormalizeSpectrum(prim,1,1);
  NormalizeSpectrum(material,1,1);
  if (status == 0) NormalizeSpectrum(result,1,1);
  //
  TH1D * sumHist = (TH1D*) prim->Clone(); // result is not equal to weighted sum ---> see FractionFitter doku and paper (?)
  sumHist->SetNameTitle("sumHist","sumHist");
  sumHist->Add(material);
  sumHist->SetLineColor(kOrange);
  //
  TCanvas * canvFeed = new TCanvas("canvFeed", "feed down and material correction");
  dcaDistr->SetMarkerStyle(24);  dcaDistr->SetMarkerSize(1.4);
  dcaDistr->DrawCopy("Ep");
  prim->DrawCopy("same");
  material->DrawCopy("same");
  sumHist->DrawCopy("same");
  //
  // Final-Results
  //
  Int_t fractionBinLow = dcaDistr->FindBin(-dcaCutInAnalysis);
  Int_t fractionBinUp  = dcaDistr->FindBin(+dcaCutInAnalysis);
  //
  Float_t materialFraction = 0.;
  if (status == 0) {
    result->DrawCopy("same");
    materialFraction = material->Integral(fractionBinLow,fractionBinUp)/dcaDistr->Integral(fractionBinLow,fractionBinUp);
  }
  Printf("RESULTING MATERIAL FRACTION: %f", materialFraction);
  Printf("TEMPLATE  MATERIAL FRACTION: %f", yieldMat);
  //
  delete sumHist;
  delete result;
  delete mcShapes;
  delete fit;
  //  delete histData;
  //delete histMC;
  delete listData;
  delete listMC;
  inFileData->Close();
  inFileMC->Close();
  delete inFileData;
  delete inFileMC;
  //
  return materialFraction;

}