[u/mrichter/AliRoot.git] / PWGPP / TPC / macros / PIDCalib / correctShapeEtaTree.C

#include "TTree.h"
#include "TH2D.h"
#include "TList.h"
#include "TString.h"
#include "TFile.h"
#include "TGraph.h"
#include "TMath.h"
#include "TDatime.h"
#include "TSpline.h"
#include "TF1.h"

#include "AliPID.h"

#include <iostream>

#include "ProgressBar.h"

Double_t getEtaCorrection(TH2D* hMap, Double_t tanTheta, Double_t dEdx)
{
  if (!hMap) 
    return 1.;
    
  if (dEdx < 1.)
    return 1.;
  
  
  //TODO WARNING: Here we can use tanTheta directly but in data we have to extract tanTheta from thetaGlobal!!!
    
  /*
  // For ESD tracks, the local tanTheta could be used (esdTrack->GetInnerParam()->GetTgl()).
  // However, this value is not available for AODs and, thus, not for AliVTrack.
  // Fortunately, the following formula allows to approximate the local tanTheta with the 
  // global theta angle -> This is for by far most of the tracks the same, but gives at
  // maybe the percent level differences within +- 0.2 in tanTheta -> Which is still ok.
  Double_t tanThetaLocal = TMath::Tan(-thetaGlobal + TMath::Pi() / 2.0);
  */
  Int_t binX = hMap->GetXaxis()->FindBin(tanTheta);
  Int_t binY = hMap->GetYaxis()->FindBin(1. / dEdx);
  
  if (binX == 0) 
    binX = 1;
  if (binX > hMap->GetXaxis()->GetNbins())
    binX = hMap->GetXaxis()->GetNbins();
  
  if (binY == 0)
    binY = 1;
  if (binY > hMap->GetYaxis()->GetNbins())
    binY = hMap->GetYaxis()->GetNbins();
  
  return hMap->GetBinContent(binX, binY);
}

Int_t correctShapeEtaTree(Bool_t correctData, TString pathMap, TString fileNameMap, TString mapSuffix,
                          Bool_t recalculateExpecteddEdx, TString pathNameSplinesFile, TString prSplinesName,
                          TString pathTree,
                          Bool_t hasMultiplicity,
                          Bool_t correctMult,
                          TString fileNameTree = "bhess_PIDetaTree.root", TString treeName = "fTree") 
{ 
  const Double_t massProton = AliPID::ParticleMass(AliPID::kProton);
  
  if (!correctData && !recalculateExpecteddEdx) {
    std::cout << "Nothing to be done: Correction AND recalculation of expected dEdx are disabled!" << std::endl;
    return -1;
  }
  // Extract the splines, if desired
  TSpline3* splPr = 0x0;
  if (recalculateExpecteddEdx) {
    std::cout << "Loading splines to recalculate expected dEdx!" << std::endl << std::endl;
    
    TFile* fSpl = TFile::Open(pathNameSplinesFile.Data());
    if (!fSpl) {
      std::cout << "Failed to open spline file \"" << pathNameSplinesFile.Data() << "\"!" << std::endl;
      return 0x0;
    }
    
    TObjArray* TPCPIDResponse = (TObjArray*)fSpl->Get("TPCPIDResponse");
    if (!TPCPIDResponse) {
      splPr = (TSpline3*)fSpl->Get(prSplinesName.Data());
      
      // If splines are in file directly, without TPCPIDResponse object, try to load them
      if (!splPr) {
        std::cout << "Failed to load object array from spline file \"" << pathNameSplinesFile.Data() << "\"!" << std::endl;
        return 0x0;
      }
    }
    else {
      splPr = (TSpline3*)TPCPIDResponse->FindObject(prSplinesName.Data());
      
      if (!splPr) {
        std::cout << "Failed to load splines from file \"" << pathNameSplinesFile.Data() << "\"!" << std::endl;
        return 0x0;
      }
    }
  }
  else
    std::cout << "Taking dEdxExpected from Tree..." << std::endl << std::endl;
  
  if (correctMult)
    std::cout << "Correcting multiplicity dependence..." << std::endl << std::endl;
  else
    std::cout << "NOT correcting multiplicity dependence..." << std::endl << std::endl;
  
  // Extract the correction map, if desired
  TH2D* hMap = 0x0;
  if (correctData) {
    std::cout << "Loading map to correct dEdx (NOT dEdx_expected in this special case!)!" << std::endl << std::endl;
    
    TFile* fMap = TFile::Open(Form("%s/%s", pathMap.Data(), fileNameMap.Data()));
    if (!fMap)  {
      std::cout << "Failed to open file \"" << Form("%s/%s", pathMap.Data(), fileNameMap.Data()) << "\"!" << std::endl;
      return -1;
    }
    
    hMap = dynamic_cast<TH2D*>(fMap->Get(Form("hRefined%s", mapSuffix.Data())));
    if (!hMap) {
      std::cout << "Failed to load correction map!" << std::endl;
      return -1;
    }
  }
  else
    std::cout << "NOT correcting eta dependence..." << std::endl << std::endl;
  
  // Extract the data Tree		       
  TFile* fTree = TFile::Open(Form("%s/%s", pathTree.Data(), fileNameTree.Data()));
  if (!fTree)  {
    std::cout << "Failed to open file \"" << Form("%s/%s", pathTree.Data(), fileNameTree.Data()) << "\"!" << std::endl;
    return -1;
  }
  
  TTree* tree = dynamic_cast<TTree*>(fTree->Get(treeName.Data()));
  if (!tree) {
    std::cout << "Failed to load data tree!" << std::endl << std::endl;
    return -1;
  }

  Long64_t nTreeEntries = tree->GetEntriesFast();
  
  Double_t pTPC = 0.; // Momentum at TPC inner wall
  //Double_t pT = 0.;
  Double_t dEdx = 0.; // Measured dE/dx
  Double_t dEdxExpected = 0.; // Expected dE/dx according to parametrisation
  Double_t tanTheta = 0.; // Tangens of (local) theta at TPC inner wall
  //Double_t sinAlpha = 0.; // Sine of (local) phi at TPC inner wall
  //Double_t y = 0.; // Local Y
  Double_t phiPrime = 0; // Phi prime
  UShort_t tpcSignalN = 0; // Number of TPC clusters for PID
  Int_t multiplicity = 0;
  UChar_t  pidType = 0; // Type of identification (TPC dEdx, V0, ...)
  
  tree->SetBranchAddress("pTPC", &pTPC);
  //tree->SetBranchAddress("pT", &Pt);
  tree->SetBranchAddress("dEdx", &dEdx);
  if (!recalculateExpecteddEdx)
    tree->SetBranchAddress("dEdxExpected", &dEdxExpected);
  tree->SetBranchAddress("tanTheta", &tanTheta);
  //tree->SetBranchAddress("sinAlpha", &sinAlpha);
  //tree->SetBranchAddress("y", &y);
  //TODO not needed for the moment tree->SetBranchAddress("phiPrime", &phiPrime);
  tree->SetBranchAddress("tpcSignalN", &tpcSignalN);
  if (hasMultiplicity)
    tree->SetBranchAddress("fMultiplicity", &multiplicity);
  tree->SetBranchAddress("pidType", &pidType);
  
  // Output file
  TDatime daTime;  
  TString savefileName = "";
  if (correctData) {
    savefileName = Form("%s_%sCorrectedWithMap_%s___%04d_%02d_%02d__%02d_%02d.root", fileNameTree.ReplaceAll(".root", "").Data(),
                        correctMult ? "multiplicityCorrected_" : "",
                        fileNameMap.ReplaceAll(".root", "").Data(), daTime.GetYear(), 
                        daTime.GetMonth(), daTime.GetDay(), daTime.GetHour(), daTime.GetMinute());
  }
  else if (recalculateExpecteddEdx) {
    savefileName = Form("%s_%sNewSplines___%04d_%02d_%02d__%02d_%02d.root", fileNameTree.ReplaceAll(".root", "").Data(),
                        correctMult ? "multiplicityCorrected_" : "",
                        daTime.GetYear(), daTime.GetMonth(), daTime.GetDay(), daTime.GetHour(), daTime.GetMinute());
  }
  else
    return -1;
    
    
  TFile* fSave = TFile::Open(Form("%s/%s", pathTree.Data(), savefileName.Data()), "recreate");
  if (!fSave) {
    std::cout << "Failed to open save file \"" << Form("%s/%s", pathTree.Data(), savefileName.Data()) << "\"!" << std::endl;
    return -1;
  }
  
  
  fSave->cd();
  TTree* treeCorrected = new TTree("fTree", "Tree for analysis of #eta dependence of TPC signal; #eta corrected");
  treeCorrected->Write(0, TObject::kOverwrite);
  
  treeCorrected->Branch("pTPC", &pTPC);
  //treeCorrected->Branch("pT", &Pt);
  treeCorrected->Branch("dEdx", &dEdx);
  treeCorrected->Branch("dEdxExpected", &dEdxExpected);
  treeCorrected->Branch("tpcSignalN", &tpcSignalN);
  
  treeCorrected->Branch("tanTheta", &tanTheta);
  //treeCorrected->Branch("sinAlpha", &sinAlpha);
  //treeCorrected->Branch("y", &y);
  treeCorrected->Branch("phiPrime", &phiPrime);
  if (hasMultiplicity)
    treeCorrected->Branch("fMultiplicity", &multiplicity);
  treeCorrected->Branch("pidType", &pidType);
  
  Double_t corrFactor = 0;
  
  
  TF1 corrFuncMult("corrFuncMult", "[0] + [1]*TMath::Max([4], TMath::Min(x, [3])) + [2] * TMath::Power(TMath::Max([4], TMath::Min(x, [3])), 2)",
                   0., 0.2);
  TF1 corrFuncMultTanTheta("corrFuncMultTanTheta", "[0] * (x -[2]) + [1] * (x * x - [2] * [2])", -1.5, 1.5);
  //OLD TF1 corrFuncMult("corrFuncMult", "[0] + [1]*TMath::Min(x, [3]) + [2] * TMath::Power(TMath::Min(x, [3]), 2)", 0., 0.1);
  //OLD acceptable try, but with coarser tanTheta binning and absTanTheta TF1 corrFuncMultTanTheta("corrFuncMultTanTheta", "[0] * (TMath::Abs(x) -TMath::Abs([2])) + [1] * (x * x - [2] * [2])", -1.5, 1.5);
  
  // LHC13b.pass2
  if (correctMult)
    printf("Using corr Parameters for 13b.pass2\n!");
  
  corrFuncMult.SetParameter(0, -6.27187e-06);
  corrFuncMult.SetParameter(1, -4.60649e-04);
  corrFuncMult.SetParameter(2, -4.26450e-02);
  corrFuncMult.SetParameter(3, 2.40590e-02);
  corrFuncMult.SetParameter(4, 0);
  
  corrFuncMultTanTheta.SetParameter(0, -5.338e-06);
  corrFuncMultTanTheta.SetParameter(1,  1.220e-05);
  corrFuncMultTanTheta.SetParameter(2, -0.5);
  
  
  /*
  // LHC11a10a
  if (correctMult)
    printf("Using corr Parameters for 11a10a\n!");
  
  corrFuncMult.SetParameter(0, 6.90133e-06);
  corrFuncMult.SetParameter(1, -1.22123e-03);
  corrFuncMult.SetParameter(2, 1.80220e-02);
  corrFuncMult.SetParameter(3, 0.1);
  corrFuncMult.SetParameter(4, 6.45306e-03);
  
  corrFuncMultTanTheta.SetParameter(0, -2.85505e-07);
  corrFuncMultTanTheta.SetParameter(1, -1.31911e-06);
  corrFuncMultTanTheta.SetParameter(2, -0.5);
  
  /* OLD very good try, but with fewer pBins for the fitting
  corrFuncMult.SetParameter(0, 6.88365e-06);
  corrFuncMult.SetParameter(1, -1.22324e-03);
  corrFuncMult.SetParameter(2, 1.81625e-02);
  corrFuncMult.SetParameter(3, 0.1);
  corrFuncMult.SetParameter(4, 6.36890e-03);
  
  corrFuncMultTanTheta.SetParameter(0, -2.85505e-07);
  corrFuncMultTanTheta.SetParameter(1, -1.31911e-06);
  corrFuncMultTanTheta.SetParameter(2, -0.5);
  */
  /*OLD good try
  corrFuncMult.SetParameter(0, 7.50321e-06);
  corrFuncMult.SetParameter(1, -1.25250e-03);
  corrFuncMult.SetParameter(2, 1.85437e-02);
  corrFuncMult.SetParameter(3, 0.1);
  corrFuncMult.SetParameter(4, 6.21192e-03);
  
  corrFuncMultTanTheta.SetParameter(0, -1.43112e-07);
  corrFuncMultTanTheta.SetParameter(1, -1.53e-06);
  corrFuncMultTanTheta.SetParameter(2, 0.3);
  */
  /* OLD acceptable try, but with coarser tanTheta binning and absTanTheta
  corrFuncMult.SetParameter(0, 6.78255e-6);
  corrFuncMult.SetParameter(1, -0.00117312);
  corrFuncMult.SetParameter(2, 0.0162423);
  corrFuncMult.SetParameter(3, 0.0563968);
  corrFuncMult.SetParameter(4, 0.00663576);
  
  corrFuncMultTanTheta.SetParameter(0, -1.85779e-6);
  corrFuncMultTanTheta.SetParameter(1, 5.40642e-7);
  corrFuncMultTanTheta.SetParameter(2, 0.35);
  */
  
  /*OLD
  corrFuncMult.SetParameter(0, 6.798e-6);
  corrFuncMult.SetParameter(1, -0.001176);
  corrFuncMult.SetParameter(2, 0.01603);
  corrFuncMult.SetParameter(3, 0.1955);
  */
  
  
  /*
  // LHC10h.pass2
  if (correctMult)
    printf("Using corr Parameters for 10h.pass2\n!");
  
  corrFuncMult.SetParameter(0, 3.21636e-07);
  corrFuncMult.SetParameter(1, -6.65876e-04);
  corrFuncMult.SetParameter(2, 1.28786e-03);
  corrFuncMult.SetParameter(3, 1.47677e-02);
  corrFuncMult.SetParameter(4, 0.);
  
  corrFuncMultTanTheta.SetParameter(0, 7.23591e-08);
  corrFuncMultTanTheta.SetParameter(1, 2.7469e-06);
  corrFuncMultTanTheta.SetParameter(2, -0.5);
  */
  /*OLD bad try
  corrFuncMult.SetParameter(0, 2.71514e-07);
  corrFuncMult.SetParameter(1, -6.92031e-04);
  corrFuncMult.SetParameter(2, 3.56042e-03);
  corrFuncMult.SetParameter(3, 1.47497e-02);
  corrFuncMult.SetParameter(4, 0.);
  
  corrFuncMultTanTheta.SetParameter(0, 8.53204e-08);
  corrFuncMultTanTheta.SetParameter(1, 2.85591e-06);
  corrFuncMultTanTheta.SetParameter(2, -0.5);
  */
  
  /*OLD OLD acceptable try, but with coarser tanTheta binning and absTanTheta
  corrFuncMult.SetParameter(0, 1.167e-6);
  corrFuncMult.SetParameter(1, -0.0009747);
  corrFuncMult.SetParameter(2, 0.02117);
  corrFuncMult.SetParameter(3, 0.01778);
  corrFuncMult.SetParameter(4, 0.);
  
  corrFuncMultTanTheta.SetParameter(0, 7.036e-7);
  corrFuncMultTanTheta.SetParameter(1, 1.868e-6);
  corrFuncMultTanTheta.SetParameter(2, 0.5);
  */
  
  
  /* OLD definition of multiplicity with nContriutorsToPrimVertex
  TF1 corrFuncMult("corrFuncMult", "pol2", 0, 0.1);
  corrFuncMult.SetParameter(0, 5.5e-6);
  corrFuncMult.SetParameter(1, -0.00436);
  corrFuncMult.SetParameter(2, 0.103);
  */

  progressbar(0.);
  for (Long64_t i = 0; i < nTreeEntries; i++) {
    tree->GetEntry(i);

    if (recalculateExpecteddEdx) {
      //Double_t old = dEdxExpected;
      dEdxExpected = 50. * splPr->Eval(pTPC / massProton); //WARNING: What, if MIP is different from 50.? Seems not to be used (tested for pp, MC_pp, PbPb and MC_PbPb), but can in principle happen
      /*if (TMath::Abs(dEdxExpected - old) > 1e-4) {
       p *rintedSomething = kTRUE;
       printf("%.2f - ", dEdxExpected);
       printf("%.2f\n", old);
      }*/
    }
    
    if (correctData) {
      corrFactor = getEtaCorrection(hMap, tanTheta, dEdxExpected);
      
      if (corrFactor <= 0)  {
        printf("Error: Bad correction factor (%f)\n", corrFactor);
        printedSomething = kTRUE;
        continue;
      }
      
      // Correct the track dEdx and leave the expected dEdx as it is, when creating the sigma map!
      // NOTE: This is due to the method the maps are created. The track dEdx (not the expected one!)
      // is corrected to uniquely relate a momemtum bin with an expected dEdx, where the expected dEdx
      // equals the track dEdx for all eta (thanks to the correction and the re-fit of the splines).
      // Consequently, looking up the uncorrected expected dEdx at a given tanTheta yields the correct
      // sigma parameter!
      
      // In principle, during map creation, one could also calculate for each tanTheta-p bin pair the expected dEdx
      // and then fill the map at the corresponding dEdxExpected-tanTheta bin. Then, one would also need to correct
      // the dEdxExpected to create the sigma maps.
      
      // In summary, both correction methods are valid in this case, but one has to be consistent! For the different methods,
      // the maps will be slightly distorted, but overall give the same results.
      
      
      
      //Scale eta dependence -> Doesn't seem to work!!!
      //corrFactor = (corrFactor - 1.)*(2e-6*30*multiplicity + 1.0) + 1.0; 
      
      dEdx /= corrFactor;
    }
    
    if (correctMult) {
      // Multiplicity depends on pure dEdx. Therefore, correction factor depends indirectly on eta
      // => Use eta correction factor to scale dEdxExpected accordingly
      Double_t dEdxExpectedInv = 1. / (dEdxExpected * (correctData ? corrFactor : 1.));
      Double_t relSlope = corrFuncMult.Eval(dEdxExpectedInv);
      
      //Correct eta dependence of slope
      relSlope += corrFuncMultTanTheta(tanTheta);
      
      Double_t corrFactorMult = relSlope * multiplicity;
      dEdx /= 1. + corrFactorMult;
    }
    
    treeCorrected->Fill();
    
    if (i % 10000 == 0)
      progressbar(100. * (((Double_t)i) / nTreeEntries));
  }
  
  progressbar(100.);
  
  fSave->cd();
  treeCorrected->Write(0, TObject::kOverwrite);
  fSave->Close();
  
  return 0;
}
Commit	Line	Data
88b71b9f	1	#include "TTree.h"
	2	#include "TH2D.h"
	3	#include "TList.h"
	4	#include "TString.h"
	5	#include "TFile.h"
	6	#include "TGraph.h"
	7	#include "TMath.h"
	8	#include "TDatime.h"
	9	#include "TSpline.h"
	10	#include "TF1.h"
	11
	12	#include "AliPID.h"
	13
	14	#include <iostream>
	15
	16	#include "ProgressBar.h"
	17
	18	Double_t getEtaCorrection(TH2D* hMap, Double_t tanTheta, Double_t dEdx)
	19	{
	20	if (!hMap)
	21	return 1.;
	22
	23	if (dEdx < 1.)
	24	return 1.;
	25
	26
	27	//TODO WARNING: Here we can use tanTheta directly but in data we have to extract tanTheta from thetaGlobal!!!
	28
	29	/*
	30	// For ESD tracks, the local tanTheta could be used (esdTrack->GetInnerParam()->GetTgl()).
	31	// However, this value is not available for AODs and, thus, not for AliVTrack.
	32	// Fortunately, the following formula allows to approximate the local tanTheta with the
	33	// global theta angle -> This is for by far most of the tracks the same, but gives at
	34	// maybe the percent level differences within +- 0.2 in tanTheta -> Which is still ok.
	35	Double_t tanThetaLocal = TMath::Tan(-thetaGlobal + TMath::Pi() / 2.0);
	36	*/
	37	Int_t binX = hMap->GetXaxis()->FindBin(tanTheta);
	38	Int_t binY = hMap->GetYaxis()->FindBin(1. / dEdx);
	39
	40	if (binX == 0)
	41	binX = 1;
	42	if (binX > hMap->GetXaxis()->GetNbins())
	43	binX = hMap->GetXaxis()->GetNbins();
	44
	45	if (binY == 0)
	46	binY = 1;
	47	if (binY > hMap->GetYaxis()->GetNbins())
	48	binY = hMap->GetYaxis()->GetNbins();
	49
	50	return hMap->GetBinContent(binX, binY);
	51	}
	52
	53	Int_t correctShapeEtaTree(Bool_t correctData, TString pathMap, TString fileNameMap, TString mapSuffix,
	54	Bool_t recalculateExpecteddEdx, TString pathNameSplinesFile, TString prSplinesName,
	55	TString pathTree,
	56	Bool_t hasMultiplicity,
	57	Bool_t correctMult,
	58	TString fileNameTree = "bhess_PIDetaTree.root", TString treeName = "fTree")
	59	{
	60	const Double_t massProton = AliPID::ParticleMass(AliPID::kProton);
	61
	62	if (!correctData && !recalculateExpecteddEdx) {
	63	std::cout << "Nothing to be done: Correction AND recalculation of expected dEdx are disabled!" << std::endl;
	64	return -1;
65	}
66	// Extract the splines, if desired
67	TSpline3* splPr = 0x0;
68	if (recalculateExpecteddEdx) {
69	std::cout << "Loading splines to recalculate expected dEdx!" << std::endl << std::endl;
70
71	TFile* fSpl = TFile::Open(pathNameSplinesFile.Data());
72	if (!fSpl) {
73	std::cout << "Failed to open spline file \"" << pathNameSplinesFile.Data() << "\"!" << std::endl;
74	return 0x0;
75	}
76
77	TObjArray* TPCPIDResponse = (TObjArray*)fSpl->Get("TPCPIDResponse");
78	if (!TPCPIDResponse) {
79	splPr = (TSpline3*)fSpl->Get(prSplinesName.Data());
80
81	// If splines are in file directly, without TPCPIDResponse object, try to load them
82	if (!splPr) {
83	std::cout << "Failed to load object array from spline file \"" << pathNameSplinesFile.Data() << "\"!" << std::endl;
84	return 0x0;
85	}
86	}
87	else {
88	splPr = (TSpline3*)TPCPIDResponse->FindObject(prSplinesName.Data());
89
90	if (!splPr) {
91	std::cout << "Failed to load splines from file \"" << pathNameSplinesFile.Data() << "\"!" << std::endl;
92	return 0x0;
93	}
94	}
95	}
96	else
97	std::cout << "Taking dEdxExpected from Tree..." << std::endl << std::endl;
98
99	if (correctMult)
100	std::cout << "Correcting multiplicity dependence..." << std::endl << std::endl;
101	else
102	std::cout << "NOT correcting multiplicity dependence..." << std::endl << std::endl;
103
104	// Extract the correction map, if desired
105	TH2D* hMap = 0x0;
106	if (correctData) {
107	std::cout << "Loading map to correct dEdx (NOT dEdx_expected in this special case!)!" << std::endl << std::endl;
108
109	TFile* fMap = TFile::Open(Form("%s/%s", pathMap.Data(), fileNameMap.Data()));
110	if (!fMap) {
111	std::cout << "Failed to open file \"" << Form("%s/%s", pathMap.Data(), fileNameMap.Data()) << "\"!" << std::endl;
112	return -1;
113	}
114
115	hMap = dynamic_cast<TH2D*>(fMap->Get(Form("hRefined%s", mapSuffix.Data())));
116	if (!hMap) {
117	std::cout << "Failed to load correction map!" << std::endl;
118	return -1;
119	}
120	}
121	else
122	std::cout << "NOT correcting eta dependence..." << std::endl << std::endl;
123
124	// Extract the data Tree
125	TFile* fTree = TFile::Open(Form("%s/%s", pathTree.Data(), fileNameTree.Data()));
126	if (!fTree) {
127	std::cout << "Failed to open file \"" << Form("%s/%s", pathTree.Data(), fileNameTree.Data()) << "\"!" << std::endl;
128	return -1;
129	}
130
131	TTree* tree = dynamic_cast<TTree*>(fTree->Get(treeName.Data()));
132	if (!tree) {
133	std::cout << "Failed to load data tree!" << std::endl << std::endl;
134	return -1;
135	}
136
137	Long64_t nTreeEntries = tree->GetEntriesFast();
138
139	Double_t pTPC = 0.; // Momentum at TPC inner wall
140	//Double_t pT = 0.;
141	Double_t dEdx = 0.; // Measured dE/dx
142	Double_t dEdxExpected = 0.; // Expected dE/dx according to parametrisation
143	Double_t tanTheta = 0.; // Tangens of (local) theta at TPC inner wall
144	//Double_t sinAlpha = 0.; // Sine of (local) phi at TPC inner wall
145	//Double_t y = 0.; // Local Y
146	Double_t phiPrime = 0; // Phi prime
147	UShort_t tpcSignalN = 0; // Number of TPC clusters for PID
148	Int_t multiplicity = 0;
149	UChar_t pidType = 0; // Type of identification (TPC dEdx, V0, ...)
150
151	tree->SetBranchAddress("pTPC", &pTPC);
152	//tree->SetBranchAddress("pT", &Pt);
153	tree->SetBranchAddress("dEdx", &dEdx);
154	if (!recalculateExpecteddEdx)
155	tree->SetBranchAddress("dEdxExpected", &dEdxExpected);
156	tree->SetBranchAddress("tanTheta", &tanTheta);
157	//tree->SetBranchAddress("sinAlpha", &sinAlpha);
158	//tree->SetBranchAddress("y", &y);
159	//TODO not needed for the moment tree->SetBranchAddress("phiPrime", &phiPrime);
160	tree->SetBranchAddress("tpcSignalN", &tpcSignalN);
161	if (hasMultiplicity)
162	tree->SetBranchAddress("fMultiplicity", &multiplicity);
163	tree->SetBranchAddress("pidType", &pidType);
164
165	// Output file
166	TDatime daTime;
167	TString savefileName = "";
168	if (correctData) {
169	savefileName = Form("%s_%sCorrectedWithMap_%s___%04d_%02d_%02d__%02d_%02d.root", fileNameTree.ReplaceAll(".root", "").Data(),
170	correctMult ? "multiplicityCorrected_" : "",
171	fileNameMap.ReplaceAll(".root", "").Data(), daTime.GetYear(),
172	daTime.GetMonth(), daTime.GetDay(), daTime.GetHour(), daTime.GetMinute());
173	}
174	else if (recalculateExpecteddEdx) {
175	savefileName = Form("%s_%sNewSplines___%04d_%02d_%02d__%02d_%02d.root", fileNameTree.ReplaceAll(".root", "").Data(),
176	correctMult ? "multiplicityCorrected_" : "",
177	daTime.GetYear(), daTime.GetMonth(), daTime.GetDay(), daTime.GetHour(), daTime.GetMinute());
178	}
179	else
180	return -1;
181
182
183	TFile* fSave = TFile::Open(Form("%s/%s", pathTree.Data(), savefileName.Data()), "recreate");
184	if (!fSave) {
185	std::cout << "Failed to open save file \"" << Form("%s/%s", pathTree.Data(), savefileName.Data()) << "\"!" << std::endl;
186	return -1;
187	}
188
189
190	fSave->cd();
191	TTree* treeCorrected = new TTree("fTree", "Tree for analysis of #eta dependence of TPC signal; #eta corrected");
192	treeCorrected->Write(0, TObject::kOverwrite);
193
194	treeCorrected->Branch("pTPC", &pTPC);
195	//treeCorrected->Branch("pT", &Pt);
196	treeCorrected->Branch("dEdx", &dEdx);
197	treeCorrected->Branch("dEdxExpected", &dEdxExpected);
198	treeCorrected->Branch("tpcSignalN", &tpcSignalN);
199
200	treeCorrected->Branch("tanTheta", &tanTheta);
201	//treeCorrected->Branch("sinAlpha", &sinAlpha);
202	//treeCorrected->Branch("y", &y);
203	treeCorrected->Branch("phiPrime", &phiPrime);
204	if (hasMultiplicity)
205	treeCorrected->Branch("fMultiplicity", &multiplicity);
206	treeCorrected->Branch("pidType", &pidType);
207
208	Double_t corrFactor = 0;
209
210
211	TF1 corrFuncMult("corrFuncMult", "[0] + [1]TMath::Max([4], TMath::Min(x, [3])) + [2] TMath::Power(TMath::Max([4], TMath::Min(x, [3])), 2)",
212	0., 0.2);
213	TF1 corrFuncMultTanTheta("corrFuncMultTanTheta", "[0] * (x -[2]) + [1] * (x * x - [2] * [2])", -1.5, 1.5);
214	//OLD TF1 corrFuncMult("corrFuncMult", "[0] + [1]TMath::Min(x, [3]) + [2] TMath::Power(TMath::Min(x, [3]), 2)", 0., 0.1);
215	//OLD acceptable try, but with coarser tanTheta binning and absTanTheta TF1 corrFuncMultTanTheta("corrFuncMultTanTheta", "[0] * (TMath::Abs(x) -TMath::Abs([2])) + [1] * (x * x - [2] * [2])", -1.5, 1.5);
216
217	// LHC13b.pass2
218	if (correctMult)
219	printf("Using corr Parameters for 13b.pass2\n!");
220
221	corrFuncMult.SetParameter(0, -6.27187e-06);
222	corrFuncMult.SetParameter(1, -4.60649e-04);
223	corrFuncMult.SetParameter(2, -4.26450e-02);
224	corrFuncMult.SetParameter(3, 2.40590e-02);
225	corrFuncMult.SetParameter(4, 0);
226
227	corrFuncMultTanTheta.SetParameter(0, -5.338e-06);
228	corrFuncMultTanTheta.SetParameter(1, 1.220e-05);
229	corrFuncMultTanTheta.SetParameter(2, -0.5);
230
231
232	/*
233	// LHC11a10a
234	if (correctMult)
235	printf("Using corr Parameters for 11a10a\n!");
236
237	corrFuncMult.SetParameter(0, 6.90133e-06);
238	corrFuncMult.SetParameter(1, -1.22123e-03);
239	corrFuncMult.SetParameter(2, 1.80220e-02);
240	corrFuncMult.SetParameter(3, 0.1);
241	corrFuncMult.SetParameter(4, 6.45306e-03);
242
243	corrFuncMultTanTheta.SetParameter(0, -2.85505e-07);
244	corrFuncMultTanTheta.SetParameter(1, -1.31911e-06);
245	corrFuncMultTanTheta.SetParameter(2, -0.5);
246
247	/* OLD very good try, but with fewer pBins for the fitting
248	corrFuncMult.SetParameter(0, 6.88365e-06);
249	corrFuncMult.SetParameter(1, -1.22324e-03);
250	corrFuncMult.SetParameter(2, 1.81625e-02);
251	corrFuncMult.SetParameter(3, 0.1);
252	corrFuncMult.SetParameter(4, 6.36890e-03);
253
254	corrFuncMultTanTheta.SetParameter(0, -2.85505e-07);
255	corrFuncMultTanTheta.SetParameter(1, -1.31911e-06);
256	corrFuncMultTanTheta.SetParameter(2, -0.5);
257	*/
258	/*OLD good try
259	corrFuncMult.SetParameter(0, 7.50321e-06);
260	corrFuncMult.SetParameter(1, -1.25250e-03);
261	corrFuncMult.SetParameter(2, 1.85437e-02);
262	corrFuncMult.SetParameter(3, 0.1);
263	corrFuncMult.SetParameter(4, 6.21192e-03);
264
265	corrFuncMultTanTheta.SetParameter(0, -1.43112e-07);
266	corrFuncMultTanTheta.SetParameter(1, -1.53e-06);
267	corrFuncMultTanTheta.SetParameter(2, 0.3);
268	*/
269	/* OLD acceptable try, but with coarser tanTheta binning and absTanTheta
270	corrFuncMult.SetParameter(0, 6.78255e-6);
271	corrFuncMult.SetParameter(1, -0.00117312);
272	corrFuncMult.SetParameter(2, 0.0162423);
273	corrFuncMult.SetParameter(3, 0.0563968);
274	corrFuncMult.SetParameter(4, 0.00663576);
275
276	corrFuncMultTanTheta.SetParameter(0, -1.85779e-6);
277	corrFuncMultTanTheta.SetParameter(1, 5.40642e-7);
278	corrFuncMultTanTheta.SetParameter(2, 0.35);
279	*/
280
281	/*OLD
282	corrFuncMult.SetParameter(0, 6.798e-6);
283	corrFuncMult.SetParameter(1, -0.001176);
284	corrFuncMult.SetParameter(2, 0.01603);
285	corrFuncMult.SetParameter(3, 0.1955);
286	*/
287
288
289	/*
290	// LHC10h.pass2
291	if (correctMult)
292	printf("Using corr Parameters for 10h.pass2\n!");
293
294	corrFuncMult.SetParameter(0, 3.21636e-07);
295	corrFuncMult.SetParameter(1, -6.65876e-04);
296	corrFuncMult.SetParameter(2, 1.28786e-03);
297	corrFuncMult.SetParameter(3, 1.47677e-02);
298	corrFuncMult.SetParameter(4, 0.);
299
300	corrFuncMultTanTheta.SetParameter(0, 7.23591e-08);
301	corrFuncMultTanTheta.SetParameter(1, 2.7469e-06);
302	corrFuncMultTanTheta.SetParameter(2, -0.5);
303	*/
304	/*OLD bad try
305	corrFuncMult.SetParameter(0, 2.71514e-07);
306	corrFuncMult.SetParameter(1, -6.92031e-04);
307	corrFuncMult.SetParameter(2, 3.56042e-03);
308	corrFuncMult.SetParameter(3, 1.47497e-02);
309	corrFuncMult.SetParameter(4, 0.);
310
311	corrFuncMultTanTheta.SetParameter(0, 8.53204e-08);
312	corrFuncMultTanTheta.SetParameter(1, 2.85591e-06);
313	corrFuncMultTanTheta.SetParameter(2, -0.5);
314	*/
315
316	/*OLD OLD acceptable try, but with coarser tanTheta binning and absTanTheta
317	corrFuncMult.SetParameter(0, 1.167e-6);
318	corrFuncMult.SetParameter(1, -0.0009747);
319	corrFuncMult.SetParameter(2, 0.02117);
320	corrFuncMult.SetParameter(3, 0.01778);
321	corrFuncMult.SetParameter(4, 0.);
322
323	corrFuncMultTanTheta.SetParameter(0, 7.036e-7);
324	corrFuncMultTanTheta.SetParameter(1, 1.868e-6);
325	corrFuncMultTanTheta.SetParameter(2, 0.5);
326	*/
327
328
329	/* OLD definition of multiplicity with nContriutorsToPrimVertex
330	TF1 corrFuncMult("corrFuncMult", "pol2", 0, 0.1);
331	corrFuncMult.SetParameter(0, 5.5e-6);
332	corrFuncMult.SetParameter(1, -0.00436);
333	corrFuncMult.SetParameter(2, 0.103);
334	*/
335
336	progressbar(0.);
337	for (Long64_t i = 0; i < nTreeEntries; i++) {
338	tree->GetEntry(i);
339
340	if (recalculateExpecteddEdx) {
341	//Double_t old = dEdxExpected;
342	dEdxExpected = 50. * splPr->Eval(pTPC / massProton); //WARNING: What, if MIP is different from 50.? Seems not to be used (tested for pp, MC_pp, PbPb and MC_PbPb), but can in principle happen
343	/*if (TMath::Abs(dEdxExpected - old) > 1e-4) {
344	p *rintedSomething = kTRUE;
345	printf("%.2f - ", dEdxExpected);
346	printf("%.2f\n", old);
347	}*/
348	}
349
350	if (correctData) {
351	corrFactor = getEtaCorrection(hMap, tanTheta, dEdxExpected);
352
353	if (corrFactor <= 0) {
354	printf("Error: Bad correction factor (%f)\n", corrFactor);
355	printedSomething = kTRUE;
356	continue;
357	}
358
359	// Correct the track dEdx and leave the expected dEdx as it is, when creating the sigma map!
360	// NOTE: This is due to the method the maps are created. The track dEdx (not the expected one!)
361	// is corrected to uniquely relate a momemtum bin with an expected dEdx, where the expected dEdx
362	// equals the track dEdx for all eta (thanks to the correction and the re-fit of the splines).
363	// Consequently, looking up the uncorrected expected dEdx at a given tanTheta yields the correct
364	// sigma parameter!
365
366	// In principle, during map creation, one could also calculate for each tanTheta-p bin pair the expected dEdx
367	// and then fill the map at the corresponding dEdxExpected-tanTheta bin. Then, one would also need to correct
368	// the dEdxExpected to create the sigma maps.
369
370	// In summary, both correction methods are valid in this case, but one has to be consistent! For the different methods,
371	// the maps will be slightly distorted, but overall give the same results.
372
373
374
375	//Scale eta dependence -> Doesn't seem to work!!!
376	//corrFactor = (corrFactor - 1.)(2e-630*multiplicity + 1.0) + 1.0;
377
378	dEdx /= corrFactor;
379	}
380
381	if (correctMult) {
382	// Multiplicity depends on pure dEdx. Therefore, correction factor depends indirectly on eta
383	// => Use eta correction factor to scale dEdxExpected accordingly
384	Double_t dEdxExpectedInv = 1. / (dEdxExpected * (correctData ? corrFactor : 1.));
385	Double_t relSlope = corrFuncMult.Eval(dEdxExpectedInv);
386
387	//Correct eta dependence of slope
388	relSlope += corrFuncMultTanTheta(tanTheta);
389
390	Double_t corrFactorMult = relSlope * multiplicity;
391	dEdx /= 1. + corrFactorMult;
392	}
393
394	treeCorrected->Fill();
395
396	if (i % 10000 == 0)
397	progressbar(100. * (((Double_t)i) / nTreeEntries));
398	}
399
400	progressbar(100.);
401
402	fSave->cd();
403	treeCorrected->Write(0, TObject::kOverwrite);
404	fSave->Close();
405
406	return 0;
407	}