[u/mrichter/AliRoot.git] / TRD / AliTRDprobdist.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//-----------------------------------------------------------------
// Class for dE/dx and Time Bin of Max. Cluster for Electrons and 
// pions in TRD. 
// It is instantiated in class AliTRDpidESD for particle identification
// in TRD
// Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>
//-----------------------------------------------------------------

#include "AliTRDprobdist.h"

ClassImp(AliTRDprobdist)

//_________________________________________________________________________
AliTRDprobdist::AliTRDprobdist(Char_t *responseFile)
{
  //
  //  The main constructor
  //
  Char_t *partName[5] = {"electron", "muon", "pion", "kaon", "proton"};
  for(Int_t i=0; i<5; i++) {
    fpartName[i] = partName[i]; 
  }
  // ADC Gain normalization
  fADCNorm=1.0;
  fNMom=kNMom;
  // Set track momenta for which response functions are available
  Double_t trackMomentum[kNMom]= {0.6, 0.8, 1, 1.5, 2, 3, 4, 5, 6, 8, 10};
  for(Int_t imom=0; imom<fNMom; imom++) {
    fTrackMomentum[imom] = trackMomentum[imom];
  }
  // Read the histogram file
  ReadData(responseFile);
  //
}

Bool_t AliTRDprobdist::ReadData(Char_t *responseFile)
{
  //
  // Read the TRD dEdx histograms.
  //
  // Read histogram Root file  
  TFile *histFile = new TFile(responseFile);
  if (!histFile || !histFile->IsOpen()) {
    Error("AliTRDprobdist", "opening TRD histgram file %s failed", responseFile);
    return kFALSE;
  }

  // Read histograms
  Char_t text[200];
  for (Int_t imom = 0; imom < kNMom; imom++) {
    sprintf(text,"h1dEdxEL%01d",imom+1);
    fh1dEdxEL[imom] = (TH1F*)histFile->Get(text);
    fh1dEdxEL[imom]->Scale(1.0/fh1dEdxEL[imom]->Integral());

    sprintf(text,"h1dEdxPI%01d",imom+1);
    fh1dEdxPI[imom] = (TH1F*)histFile->Get(text);
    fh1dEdxPI[imom]->Scale(1.0/fh1dEdxPI[imom]->Integral());

    sprintf(text,"h1dEdxMU%01d",imom+1);
    fh1dEdxMU[imom] = (TH1F*)histFile->Get(text);
    fh1dEdxMU[imom]->Scale(1.0/fh1dEdxMU[imom]->Integral());

    sprintf(text,"h1dEdxKA%01d",imom+1);
    fh1dEdxKA[imom] = (TH1F*)histFile->Get(text);
    fh1dEdxKA[imom]->Scale(1.0/fh1dEdxKA[imom]->Integral());

    sprintf(text,"h1dEdxPR%01d",imom+1);
    fh1dEdxPR[imom] = (TH1F*)histFile->Get(text);
    fh1dEdxPR[imom]->Scale(1.0/fh1dEdxPR[imom]->Integral());

    sprintf(text,"h1MaxTimBinEL%01d",imom+1);
    fh1MaxTimBinEL[imom] = (TH1F*)histFile->Get(text);
    fh1MaxTimBinEL[imom]->Scale(1.0/fh1MaxTimBinEL[imom]->Integral());

    sprintf(text,"h1MaxTimBinPI%01d",imom+1);
    fh1MaxTimBinPI[imom] = (TH1F*)histFile->Get(text);
    fh1MaxTimBinPI[imom]->Scale(1.0/fh1MaxTimBinPI[imom]->Integral());
  }
  // Number of bins and bin size
  fNbins = fh1dEdxPI[1]->GetNbinsX();
  fBinSize = fh1dEdxPI[1]->GetBinWidth(1);
  return kTRUE;
}

AliTRDprobdist::AliTRDprobdist(const AliTRDprobdist& pd):TNamed()
{
  //
  // Copy constructor.
  //
  for(Int_t i=0; i<5; i++) {
    fpartName[i] = pd.fpartName[i]; 
  }
  // ADC Gain normalization
  fADCNorm=pd.fADCNorm;
  fNMom=pd.fNMom;
  // Set track momenta for which response functions are available
  for(Int_t imom=0; imom<fNMom; imom++) {
    fTrackMomentum[imom] = pd.fTrackMomentum[imom];
  }

  fNbins=pd.fNbins;
  fBinSize=pd.fBinSize;  

  for (Int_t imom = 0; imom < kNMom; imom++) {
    fh1dEdxEL[imom] = pd.fh1dEdxEL[imom];
    fh1dEdxPI[imom] = pd.fh1dEdxPI[imom];
    fh1dEdxMU[imom] = pd.fh1dEdxMU[imom];
    fh1dEdxKA[imom] = pd.fh1dEdxKA[imom];
    fh1dEdxPR[imom] = pd.fh1dEdxPR[imom];
    fh1MaxTimBinEL[imom] = pd.fh1MaxTimBinEL[imom];
    fh1MaxTimBinPI[imom] = pd.fh1MaxTimBinPI[imom];
  }
}


AliTRDprobdist::~AliTRDprobdist()
{
  // Destructor

}


//_________________________________________________________________________
Double_t  AliTRDprobdist::GetMean(Int_t k, Int_t ip) const
{
  //
  // Gets mean of de/dx dist. of e
  printf("Mean for particle = %s and momentum = %.2f is:\n", fpartName[k], fTrackMomentum[ip]);
  if(k==0) return fh1dEdxEL[ip]->GetMean();
  if(k==1) return fh1dEdxMU[ip]->GetMean();
  if(k==2) return fh1dEdxPI[ip]->GetMean();
  if(k==3) return fh1dEdxKA[ip]->GetMean();
  if(k==4) return fh1dEdxPR[ip]->GetMean();
  return fh1dEdxPR[ip]->GetMean();
}

//_________________________________________________________________________
Double_t  AliTRDprobdist::GetNormalization(Int_t k, Int_t ip) const
{
  //
  // Gets Normalization of de/dx dist. of e

  printf("Normalization for particle = %s and momentum = %.2f is:\n",fpartName[k], fTrackMomentum[ip]);
  if(k==0) return fh1dEdxEL[ip]->Integral();
  if(k==1) return fh1dEdxMU[ip]->Integral();
  if(k==2) return fh1dEdxPI[ip]->Integral();
  if(k==3) return fh1dEdxKA[ip]->Integral();
  if(k==4) return fh1dEdxPR[ip]->Integral();
  return fh1dEdxPR[ip]->Integral();
}

TH1F* AliTRDprobdist::GetHistogram(Int_t k, Int_t ip) const
{
  //
  //
  printf("Histogram for particle = %s and momentum = %.2f is:\n", fpartName[k], fTrackMomentum[ip]);
  if(k==0) return fh1dEdxEL[ip];
  if(k==1) return fh1dEdxMU[ip];
  if(k==2) return fh1dEdxPI[ip];
  if(k==3) return fh1dEdxKA[ip];
  if(k==4) return fh1dEdxPR[ip];
  return fh1dEdxPR[ip];
}

//_________________________________________________________________________
Double_t AliTRDprobdist::GetProbability(Int_t k, Double_t mom, Double_t dedx1) const
{
  //
  // Gets the Probability of having dedx at a given momentum (mom)
  // and particle type k (0 for e) and (2 for pi)
  // from the precalculated de/dx distributions 
  Double_t probability = 1.0;
  Double_t dedx = dedx1/fADCNorm;
  Int_t iEnBin= ((Int_t) (dedx/fBinSize+1));
  if(iEnBin > fNbins) iEnBin = fNbins;

  ////Electron//////////////////////////
  if(k==0){    // electron 
    Double_t slop;
    // Lower limit
    if(mom<=fTrackMomentum[0]) {
      slop=(fh1dEdxEL[1]->GetBinContent(iEnBin)-fh1dEdxEL[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
      probability= fh1dEdxEL[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
      return probability;
    }
    // Upper Limit
    if(mom>=fTrackMomentum[fNMom-1]) {
      slop=(fh1dEdxEL[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxEL[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
      probability= fh1dEdxEL[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
      return probability;
    }
    // In the range
    for(Int_t ip=1; ip<fNMom; ip++) 
      if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
	slop=(fh1dEdxEL[ip]->GetBinContent(iEnBin)-fh1dEdxEL[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
	// Linear Interpolation
	probability= fh1dEdxEL[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
	return probability;
      }
  }

  ////Pion//////////////////////////
  if(k==2){    // Pion
    Double_t slop;
    // Lower limit
    if(mom<=fTrackMomentum[0]) {
      slop=(fh1dEdxPI[1]->GetBinContent(iEnBin)-fh1dEdxPI[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
      probability= fh1dEdxPI[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
      return probability;
    }
    // Upper Limit
    if(mom>=fTrackMomentum[fNMom-1]) {
      slop=(fh1dEdxPI[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxPI[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
      probability= fh1dEdxPI[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
      return probability;
    }
    // In the range
    for(Int_t ip=1; ip<fNMom; ip++) 
      if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
	slop=(fh1dEdxPI[ip]->GetBinContent(iEnBin)-fh1dEdxPI[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
	// Linear Interpolation
	probability= fh1dEdxPI[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
	return probability;
      }
  }

  ////Muon//////////////////////////
  if(k==1){    // Muon
    Double_t slop;
    // Lower limit
    if(mom<=fTrackMomentum[0]) {
      slop=(fh1dEdxMU[1]->GetBinContent(iEnBin)-fh1dEdxMU[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
      probability= fh1dEdxMU[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
      return probability;
    }
    // Upper Limit
    if(mom>=fTrackMomentum[fNMom-1]) {
      slop=(fh1dEdxMU[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxMU[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
      probability= fh1dEdxMU[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
      return probability;
    }
    // In the range
    for(Int_t ip=1; ip<fNMom; ip++) 
      if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
	slop=(fh1dEdxMU[ip]->GetBinContent(iEnBin)-fh1dEdxMU[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
	// Linear Interpolation
	probability= fh1dEdxMU[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
	return probability;
      }
  }

  ////Kaon//////////////////////////
  if(k==3){    // Kaon
    Double_t slop;
    // Lower limit
    if(mom<=fTrackMomentum[0]) {
      slop=(fh1dEdxKA[1]->GetBinContent(iEnBin)-fh1dEdxKA[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
      probability= fh1dEdxKA[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
      return probability;
    }
    // Upper Limit
    if(mom>=fTrackMomentum[fNMom-1]) {
      slop=(fh1dEdxKA[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxKA[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
      probability= fh1dEdxKA[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
      return probability;
    }
    // In the range
    for(Int_t ip=1; ip<fNMom; ip++) 
      if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
	slop=(fh1dEdxKA[ip]->GetBinContent(iEnBin)-fh1dEdxKA[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
	// Linear Interpolation
	probability= fh1dEdxKA[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
	return probability;
      }
  }

  ////Proton//////////////////////////
  if(k==4){    // Proton
    Double_t slop;
    // Lower limit
    if(mom<=fTrackMomentum[0]) {
      slop=(fh1dEdxPR[1]->GetBinContent(iEnBin)-fh1dEdxPR[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
      probability= fh1dEdxPR[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
      return probability;
    }
    // Upper Limit
    if(mom>=fTrackMomentum[fNMom-1]) {
      slop=(fh1dEdxPR[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxPR[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
      probability= fh1dEdxPR[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
      return probability;
    }
    // In the range
    for(Int_t ip=1; ip<fNMom; ip++) 
      if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
	slop=(fh1dEdxPR[ip]->GetBinContent(iEnBin)-fh1dEdxPR[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
	// Linear Interpolation
	probability= fh1dEdxPR[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
	return probability;
      }
  }

  return probability;
}


//_________________________________________________________________________
Double_t AliTRDprobdist::GetProbabilityT(Int_t k, Double_t mom, Int_t timbin) const
{
  //
  // Gets the Probability of having timbin at a given momentum (mom)
  // and particle type k (0 for e) and (2 for pi)
  // from the precalculated timbin distributions 
  Double_t probabilityT = 1.0;
  if(timbin<=0) return 0.;
  Int_t iTBin=timbin+1;

  if(k==0){    // electron
    if(mom<=fTrackMomentum[0]) probabilityT = fh1MaxTimBinEL[0]->GetBinContent(iTBin);
    if(mom>=fTrackMomentum[fNMom-1]) probabilityT = fh1MaxTimBinEL[fNMom-1]->GetBinContent(iTBin);
  }
  if(k==1||k==2||k==3||k==4){    // pion
    if(mom<=fTrackMomentum[0]) probabilityT = fh1MaxTimBinPI[0]->GetBinContent(iTBin);
    if(mom>=fTrackMomentum[fNMom-1]) probabilityT = fh1MaxTimBinPI[fNMom-1]->GetBinContent(iTBin);
  }

  if(k==0)    // electron
  for(Int_t ip=1; ip<fNMom; ip++)
  if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
    Double_t slop=(fh1MaxTimBinEL[ip]->GetBinContent(iTBin)-fh1MaxTimBinEL[ip-1]->GetBinContent(iTBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
    // Linear Interpolation
    probabilityT= fh1MaxTimBinEL[ip-1]->GetBinContent(iTBin) + slop*(mom-fTrackMomentum[ip-1]);
    return probabilityT;
  }

  if(k==1||k==2||k==3||k==4)   // pion and other particles
  for(Int_t ip=1; ip<fNMom; ip++)
  if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
    Double_t slop=(fh1MaxTimBinPI[ip]->GetBinContent(iTBin)-fh1MaxTimBinPI[ip-1]->GetBinContent(iTBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
    // Linear Interpolation
    probabilityT= fh1MaxTimBinPI[ip-1]->GetBinContent(iTBin) + slop*(mom-fTrackMomentum[ip-1]);
    return probabilityT;
  }
  return probabilityT;
}
Commit	Line	Data
33a848f6	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//-----------------------------------------------------------------
	17	// Class for dE/dx and Time Bin of Max. Cluster for Electrons and
	18	// pions in TRD.
	19	// It is instantiated in class AliTRDpidESD for particle identification
	20	// in TRD
	21	// Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>
	22	//-----------------------------------------------------------------
	23
	24	#include "AliTRDprobdist.h"
	25
	26	ClassImp(AliTRDprobdist)
	27
	28	//_________________________________________________________________________
defc9040	29	AliTRDprobdist::AliTRDprobdist(Char_t *responseFile)
33a848f6	30	{
	31	//
	32	// The main constructor
	33	//
defc9040	34	Char_t *partName[5] = {"electron", "muon", "pion", "kaon", "proton"};
	35	for(Int_t i=0; i<5; i++) {
	36	fpartName[i] = partName[i];
	37	}
	38	// ADC Gain normalization
	39	fADCNorm=1.0;
	40	fNMom=kNMom;
	41	// Set track momenta for which response functions are available
	42	Double_t trackMomentum[kNMom]= {0.6, 0.8, 1, 1.5, 2, 3, 4, 5, 6, 8, 10};
	43	for(Int_t imom=0; imom<fNMom; imom++) {
	44	fTrackMomentum[imom] = trackMomentum[imom];
	45	}
	46	// Read the histogram file
	47	ReadData(responseFile);
	48	//
33a848f6	49	}
33a848f6	50
defc9040	51	Bool_t AliTRDprobdist::ReadData(Char_t *responseFile)
33a848f6	52	{
defc9040	53	//
	54	// Read the TRD dEdx histograms.
	55	//
	56	// Read histogram Root file
	57	TFile *histFile = new TFile(responseFile);
	58	if (!histFile \|\| !histFile->IsOpen()) {
	59	Error("AliTRDprobdist", "opening TRD histgram file %s failed", responseFile);
	60	return kFALSE;
33a848f6	61	}
defc9040	62
	63	// Read histograms
	64	Char_t text[200];
	65	for (Int_t imom = 0; imom < kNMom; imom++) {
	66	sprintf(text,"h1dEdxEL%01d",imom+1);
	67	fh1dEdxEL[imom] = (TH1F*)histFile->Get(text);
	68	fh1dEdxEL[imom]->Scale(1.0/fh1dEdxEL[imom]->Integral());
	69
	70	sprintf(text,"h1dEdxPI%01d",imom+1);
	71	fh1dEdxPI[imom] = (TH1F*)histFile->Get(text);
	72	fh1dEdxPI[imom]->Scale(1.0/fh1dEdxPI[imom]->Integral());
	73
	74	sprintf(text,"h1dEdxMU%01d",imom+1);
	75	fh1dEdxMU[imom] = (TH1F*)histFile->Get(text);
	76	fh1dEdxMU[imom]->Scale(1.0/fh1dEdxMU[imom]->Integral());
	77
	78	sprintf(text,"h1dEdxKA%01d",imom+1);
	79	fh1dEdxKA[imom] = (TH1F*)histFile->Get(text);
	80	fh1dEdxKA[imom]->Scale(1.0/fh1dEdxKA[imom]->Integral());
	81
	82	sprintf(text,"h1dEdxPR%01d",imom+1);
	83	fh1dEdxPR[imom] = (TH1F*)histFile->Get(text);
	84	fh1dEdxPR[imom]->Scale(1.0/fh1dEdxPR[imom]->Integral());
	85
	86	sprintf(text,"h1MaxTimBinEL%01d",imom+1);
	87	fh1MaxTimBinEL[imom] = (TH1F*)histFile->Get(text);
	88	fh1MaxTimBinEL[imom]->Scale(1.0/fh1MaxTimBinEL[imom]->Integral());
	89
	90	sprintf(text,"h1MaxTimBinPI%01d",imom+1);
	91	fh1MaxTimBinPI[imom] = (TH1F*)histFile->Get(text);
	92	fh1MaxTimBinPI[imom]->Scale(1.0/fh1MaxTimBinPI[imom]->Integral());
	93	}
	94	// Number of bins and bin size
	95	fNbins = fh1dEdxPI[1]->GetNbinsX();
	96	fBinSize = fh1dEdxPI[1]->GetBinWidth(1);
	97	return kTRUE;
33a848f6	98	}
33a848f6	99
defc9040	100	AliTRDprobdist::AliTRDprobdist(const AliTRDprobdist& pd):TNamed()
33a848f6	101	{
33a848f6	102	//
defc9040	103	// Copy constructor.
	104	//
	105	for(Int_t i=0; i<5; i++) {
	106	fpartName[i] = pd.fpartName[i];
	107	}
	108	// ADC Gain normalization
	109	fADCNorm=pd.fADCNorm;
	110	fNMom=pd.fNMom;
	111	// Set track momenta for which response functions are available
	112	for(Int_t imom=0; imom<fNMom; imom++) {
	113	fTrackMomentum[imom] = pd.fTrackMomentum[imom];
	114	}
	115
	116	fNbins=pd.fNbins;
	117	fBinSize=pd.fBinSize;
	118
	119	for (Int_t imom = 0; imom < kNMom; imom++) {
	120	fh1dEdxEL[imom] = pd.fh1dEdxEL[imom];
	121	fh1dEdxPI[imom] = pd.fh1dEdxPI[imom];
	122	fh1dEdxMU[imom] = pd.fh1dEdxMU[imom];
	123	fh1dEdxKA[imom] = pd.fh1dEdxKA[imom];
	124	fh1dEdxPR[imom] = pd.fh1dEdxPR[imom];
	125	fh1MaxTimBinEL[imom] = pd.fh1MaxTimBinEL[imom];
	126	fh1MaxTimBinPI[imom] = pd.fh1MaxTimBinPI[imom];
33a848f6	127	}
33a848f6	128	}
33a848f6	129
defc9040	130
	131	AliTRDprobdist::~AliTRDprobdist()
	132	{
	133	// Destructor
	134
	135	}
	136
	137
33a848f6	138	//_________________________________________________________________________
defc9040	139	Double_t AliTRDprobdist::GetMean(Int_t k, Int_t ip) const
33a848f6	140	{
33a848f6	141	//
defc9040	142	// Gets mean of de/dx dist. of e
	143	printf("Mean for particle = %s and momentum = %.2f is:\n", fpartName[k], fTrackMomentum[ip]);
	144	if(k==0) return fh1dEdxEL[ip]->GetMean();
	145	if(k==1) return fh1dEdxMU[ip]->GetMean();
	146	if(k==2) return fh1dEdxPI[ip]->GetMean();
	147	if(k==3) return fh1dEdxKA[ip]->GetMean();
	148	if(k==4) return fh1dEdxPR[ip]->GetMean();
	149	return fh1dEdxPR[ip]->GetMean();
33a848f6	150	}
	151
	152	//_________________________________________________________________________
defc9040	153	Double_t AliTRDprobdist::GetNormalization(Int_t k, Int_t ip) const
33a848f6	154	{
	155	//
	156	// Gets Normalization of de/dx dist. of e
defc9040	157
	158	printf("Normalization for particle = %s and momentum = %.2f is:\n",fpartName[k], fTrackMomentum[ip]);
	159	if(k==0) return fh1dEdxEL[ip]->Integral();
	160	if(k==1) return fh1dEdxMU[ip]->Integral();
	161	if(k==2) return fh1dEdxPI[ip]->Integral();
	162	if(k==3) return fh1dEdxKA[ip]->Integral();
	163	if(k==4) return fh1dEdxPR[ip]->Integral();
	164	return fh1dEdxPR[ip]->Integral();
	165	}
	166
	167	TH1F* AliTRDprobdist::GetHistogram(Int_t k, Int_t ip) const
	168	{
	169	//
	170	//
	171	printf("Histogram for particle = %s and momentum = %.2f is:\n", fpartName[k], fTrackMomentum[ip]);
	172	if(k==0) return fh1dEdxEL[ip];
	173	if(k==1) return fh1dEdxMU[ip];
	174	if(k==2) return fh1dEdxPI[ip];
	175	if(k==3) return fh1dEdxKA[ip];
	176	if(k==4) return fh1dEdxPR[ip];
	177	return fh1dEdxPR[ip];
33a848f6	178	}
	179
	180	//_________________________________________________________________________
defc9040	181	Double_t AliTRDprobdist::GetProbability(Int_t k, Double_t mom, Double_t dedx1) const
33a848f6	182	{
	183	//
	184	// Gets the Probability of having dedx at a given momentum (mom)
	185	// and particle type k (0 for e) and (2 for pi)
	186	// from the precalculated de/dx distributions
1642e650	187	Double_t probability = 1.0;
defc9040	188	Double_t dedx = dedx1/fADCNorm;
	189	Int_t iEnBin= ((Int_t) (dedx/fBinSize+1));
	190	if(iEnBin > fNbins) iEnBin = fNbins;
33a848f6	191
defc9040	192	////Electron//////////////////////////
defc9040	193	if(k==0){ // electron
eee41b89	194	Double_t slop;
defc9040	195	// Lower limit
eee41b89	196	if(mom<=fTrackMomentum[0]) {
defc9040	197	slop=(fh1dEdxEL[1]->GetBinContent(iEnBin)-fh1dEdxEL[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
defc9040	198	probability= fh1dEdxEL[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
eee41b89	199	return probability;
eee41b89	200	}
defc9040	201	// Upper Limit
eee41b89	202	if(mom>=fTrackMomentum[fNMom-1]) {
defc9040	203	slop=(fh1dEdxEL[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxEL[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
defc9040	204	probability= fh1dEdxEL[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
eee41b89	205	return probability;
eee41b89	206	}
defc9040	207	// In the range
	208	for(Int_t ip=1; ip<fNMom; ip++)
	209	if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
	210	slop=(fh1dEdxEL[ip]->GetBinContent(iEnBin)-fh1dEdxEL[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
	211	// Linear Interpolation
	212	probability= fh1dEdxEL[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
	213	return probability;
	214	}
33a848f6	215	}
eee41b89	216
defc9040	217	////Pion//////////////////////////
defc9040	218	if(k==2){ // Pion
eee41b89	219	Double_t slop;
defc9040	220	// Lower limit
eee41b89	221	if(mom<=fTrackMomentum[0]) {
defc9040	222	slop=(fh1dEdxPI[1]->GetBinContent(iEnBin)-fh1dEdxPI[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
defc9040	223	probability= fh1dEdxPI[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
eee41b89	224	return probability;
eee41b89	225	}
defc9040	226	// Upper Limit
eee41b89	227	if(mom>=fTrackMomentum[fNMom-1]) {
defc9040	228	slop=(fh1dEdxPI[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxPI[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
defc9040	229	probability= fh1dEdxPI[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
eee41b89	230	return probability;
eee41b89	231	}
defc9040	232	// In the range
	233	for(Int_t ip=1; ip<fNMom; ip++)
	234	if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
	235	slop=(fh1dEdxPI[ip]->GetBinContent(iEnBin)-fh1dEdxPI[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
	236	// Linear Interpolation
	237	probability= fh1dEdxPI[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
	238	return probability;
	239	}
33a848f6	240	}
33a848f6	241
defc9040	242	////Muon//////////////////////////
defc9040	243	if(k==1){ // Muon
eee41b89	244	Double_t slop;
defc9040	245	// Lower limit
eee41b89	246	if(mom<=fTrackMomentum[0]) {
defc9040	247	slop=(fh1dEdxMU[1]->GetBinContent(iEnBin)-fh1dEdxMU[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
defc9040	248	probability= fh1dEdxMU[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
eee41b89	249	return probability;
eee41b89	250	}
defc9040	251	// Upper Limit
eee41b89	252	if(mom>=fTrackMomentum[fNMom-1]) {
defc9040	253	slop=(fh1dEdxMU[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxMU[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
defc9040	254	probability= fh1dEdxMU[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
eee41b89	255	return probability;
eee41b89	256	}
defc9040	257	// In the range
	258	for(Int_t ip=1; ip<fNMom; ip++)
	259	if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
	260	slop=(fh1dEdxMU[ip]->GetBinContent(iEnBin)-fh1dEdxMU[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
	261	// Linear Interpolation
	262	probability= fh1dEdxMU[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
	263	return probability;
	264	}
eee41b89	265	}
eee41b89	266
defc9040	267	////Kaon//////////////////////////
defc9040	268	if(k==3){ // Kaon
eee41b89	269	Double_t slop;
defc9040	270	// Lower limit
eee41b89	271	if(mom<=fTrackMomentum[0]) {
defc9040	272	slop=(fh1dEdxKA[1]->GetBinContent(iEnBin)-fh1dEdxKA[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
defc9040	273	probability= fh1dEdxKA[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
eee41b89	274	return probability;
eee41b89	275	}
defc9040	276	// Upper Limit
eee41b89	277	if(mom>=fTrackMomentum[fNMom-1]) {
defc9040	278	slop=(fh1dEdxKA[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxKA[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
defc9040	279	probability= fh1dEdxKA[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
eee41b89	280	return probability;
eee41b89	281	}
defc9040	282	// In the range
	283	for(Int_t ip=1; ip<fNMom; ip++)
	284	if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
	285	slop=(fh1dEdxKA[ip]->GetBinContent(iEnBin)-fh1dEdxKA[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
	286	// Linear Interpolation
	287	probability= fh1dEdxKA[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
	288	return probability;
	289	}
eee41b89	290	}
eee41b89	291
defc9040	292	////Proton//////////////////////////
defc9040	293	if(k==4){ // Proton
eee41b89	294	Double_t slop;
defc9040	295	// Lower limit
eee41b89	296	if(mom<=fTrackMomentum[0]) {
defc9040	297	slop=(fh1dEdxPR[1]->GetBinContent(iEnBin)-fh1dEdxPR[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
defc9040	298	probability= fh1dEdxPR[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
eee41b89	299	return probability;
eee41b89	300	}
defc9040	301	// Upper Limit
eee41b89	302	if(mom>=fTrackMomentum[fNMom-1]) {
defc9040	303	slop=(fh1dEdxPR[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxPR[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
defc9040	304	probability= fh1dEdxPR[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
eee41b89	305	return probability;
eee41b89	306	}
defc9040	307	// In the range
	308	for(Int_t ip=1; ip<fNMom; ip++)
	309	if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
	310	slop=(fh1dEdxPR[ip]->GetBinContent(iEnBin)-fh1dEdxPR[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
	311	// Linear Interpolation
	312	probability= fh1dEdxPR[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
	313	return probability;
	314	}
eee41b89	315	}
eee41b89	316
1642e650	317	return probability;
33a848f6	318	}
	319
	320
	321	//_________________________________________________________________________
1642e650	322	Double_t AliTRDprobdist::GetProbabilityT(Int_t k, Double_t mom, Int_t timbin) const
33a848f6	323	{
	324	//
	325	// Gets the Probability of having timbin at a given momentum (mom)
	326	// and particle type k (0 for e) and (2 for pi)
	327	// from the precalculated timbin distributions
1642e650	328	Double_t probabilityT = 1.0;
33a848f6	329	if(timbin<=0) return 0.;
defc9040	330	Int_t iTBin=timbin+1;
33a848f6	331
33a848f6	332	if(k==0){ // electron
defc9040	333	if(mom<=fTrackMomentum[0]) probabilityT = fh1MaxTimBinEL[0]->GetBinContent(iTBin);
defc9040	334	if(mom>=fTrackMomentum[fNMom-1]) probabilityT = fh1MaxTimBinEL[fNMom-1]->GetBinContent(iTBin);
33a848f6	335	}
eee41b89	336	if(k==1\|\|k==2\|\|k==3\|\|k==4){ // pion
defc9040	337	if(mom<=fTrackMomentum[0]) probabilityT = fh1MaxTimBinPI[0]->GetBinContent(iTBin);
defc9040	338	if(mom>=fTrackMomentum[fNMom-1]) probabilityT = fh1MaxTimBinPI[fNMom-1]->GetBinContent(iTBin);
33a848f6	339	}
	340
	341	if(k==0) // electron
	342	for(Int_t ip=1; ip<fNMom; ip++)
	343	if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
defc9040	344	Double_t slop=(fh1MaxTimBinEL[ip]->GetBinContent(iTBin)-fh1MaxTimBinEL[ip-1]->GetBinContent(iTBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
33a848f6	345	// Linear Interpolation
defc9040	346	probabilityT= fh1MaxTimBinEL[ip-1]->GetBinContent(iTBin) + slop*(mom-fTrackMomentum[ip-1]);
1642e650	347	return probabilityT;
33a848f6	348	}
33a848f6	349
eee41b89	350	if(k==1\|\|k==2\|\|k==3\|\|k==4) // pion and other particles
33a848f6	351	for(Int_t ip=1; ip<fNMom; ip++)
33a848f6	352	if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
defc9040	353	Double_t slop=(fh1MaxTimBinPI[ip]->GetBinContent(iTBin)-fh1MaxTimBinPI[ip-1]->GetBinContent(iTBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
33a848f6	354	// Linear Interpolation
defc9040	355	probabilityT= fh1MaxTimBinPI[ip-1]->GetBinContent(iTBin) + slop*(mom-fTrackMomentum[ip-1]);
1642e650	356	return probabilityT;
33a848f6	357	}
1642e650	358	return probabilityT;
33a848f6	359	}
	360
	361