[u/mrichter/AliRoot.git] / STEER / AliTRDPIDResponse.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.                  *
**************************************************************************/
//
// PID Response class for the TRD detector
// Based on 1D Likelihood approach
// Calculation of probabilities using Bayesian approach
// Attention: This method is only used to separate electrons from pions
//
// Authors:
//  Markus Fasel <M.Fasel@gsi.de>
//  Anton Andronic <A.Andronic@gsi.de>
//
#include <TClass.h>
#include <TAxis.h>
#include <TFile.h>
#include <TH1.h>
#include <TKey.h>
#include <TObjArray.h>
#include <TString.h>
#include <TROOT.h> 
#include <TSystem.h>
#include <TDirectory.h>

#include "AliLog.h"

#include "AliTRDPIDResponse.h"

ClassImp(AliTRDPIDResponse)

const Double_t AliTRDPIDResponse::fgkPBins[kNPBins] = {1, 2, 3, 4, 5, 6};

//____________________________________________________________
AliTRDPIDResponse::AliTRDPIDResponse():
  TObject()
  ,fReferences(NULL)
  ,fGainNormalisationFactor(1.)
  ,fPIDmethod(kLQ1D)
{
  //
  // Default constructor
  //
  for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++)
    for(Int_t ipbin = 0; ipbin < kNPBins; ipbin++)
      fMapRefHists[ispec][ipbin] = -1;
}

//____________________________________________________________
AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
  TObject(ref)
  ,fReferences(ref.fReferences)
  ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
  ,fPIDmethod(ref.fPIDmethod)
{
  //
  // Copy constructor
  // Flat copy of the reference histos
  // For creating a deep copy call SetOwner
  //
  Int_t size  = (AliPID::kSPECIES)*(kNPBins);
  memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * size);
}

//____________________________________________________________
AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
  //
  // Assignment operator
  // Performs a flat copy of the reference histos
  //
  if(this == &ref) return *this;
  
  // Make copy
  TObject::operator=(ref);
  if(IsOwner()){
    if(fReferences){
      fReferences->Delete();
      delete fReferences;
    }
    SetBit(kIsOwner, kFALSE);
  } else if(fReferences) delete fReferences;
  fReferences = ref.fReferences;
  Int_t size  = (AliPID::kSPECIES)*(kNPBins);
  memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * size);
  fPIDmethod = ref.fPIDmethod;
  
  return *this;
}

//____________________________________________________________
AliTRDPIDResponse::~AliTRDPIDResponse(){
  //
  // Destructor
  //
  if(IsOwner()){
    // Destroy histos
    if(fReferences) fReferences->Delete();
  }
  if(fReferences) delete fReferences;
}

//____________________________________________________________
Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
  //
  // Load References into the toolkit
  //
  AliDebug(1, "Loading reference histos from root file");
  TDirectory *owd = gDirectory;// store old working directory
  
  if(!filename)
    filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
  TFile *in = TFile::Open(filename);
  if(!in){
    AliError("Ref file not available.");
    return kFALSE;
  }
  
  gROOT->cd();
  fReferences = new TObjArray(AliPID::kSPECIES*kNPBins);
  fReferences->SetOwner();
  TIter iter(in->GetListOfKeys());
  TKey *histkey = NULL;
  TObject *tmp = NULL;
  Int_t arrayPos = 0, pbin = 0; 
  AliPID::EParticleType species;
  TString histname;
  TH1 *hnew = NULL;
  while((histkey = dynamic_cast<TKey *>(iter()))){
    tmp = histkey->ReadObj();
    histname = tmp->GetName();
    if(histname.BeginsWith("fHQel")){ // Electron histogram
      histname.ReplaceAll("fHQel_p","");
      species = AliPID::kElectron;
    } else if(histname.BeginsWith("fHQmu")){ // Muon histogram
      histname.ReplaceAll("fHQmu_p","");
      species = AliPID::kMuon;
    } else if(histname.BeginsWith("fHQpi")){ // Pion histogram
      histname.ReplaceAll("fHQpi_p","");
      species = AliPID::kPion;
    }else if(histname.BeginsWith("fHQka")){ // Kaon histogram
      histname.ReplaceAll("fHQka_p","");
      species = AliPID::kKaon;
    }else if(histname.BeginsWith("fHQpr")){ // Proton histogram
      histname.ReplaceAll("fHQpr_p","");
      species = AliPID::kProton;
    } else continue;
    pbin = histname.Atoi() - 1;
    AliDebug(1, Form("Species %d, Histname %s, Pbin %d, Position in container %d", species, histname.Data(), pbin, arrayPos));
    fMapRefHists[species][pbin] = arrayPos;
    fReferences->AddAt((hnew =new TH1F(*dynamic_cast<TH1F *>(tmp))), arrayPos);
    hnew->SetDirectory(gROOT);
    arrayPos++;
  } 
  in->Close();
  owd->cd();
  AliDebug(2, Form("Successfully loaded %d Reference Histograms", arrayPos));
  SetBit(kIsOwner, kTRUE);
  delete in;
  return kTRUE;
}

//____________________________________________________________
Bool_t AliTRDPIDResponse::Load(const TObjArray *histos){
  //
  // Load Reference into the PID Response (from OCDB/OADB)
  //
  AliDebug(1, "Setting references (from the database)");
  if(fReferences) fReferences->Clear();
  else
    fReferences = new TObjArray(AliPID::kSPECIES*kNPBins);
  fReferences->SetOwner();
  TIter iter(histos);
  Int_t arrayPos = 0, pbin = 0; 
  AliPID::EParticleType species;
  TObject *tmp = NULL;
  TString histname;
  TH1 *hnew = NULL;
  while((tmp = iter())){
    histname = tmp->GetName();
    if(histname.BeginsWith("fHQel")){ // Electron histogram
      histname.ReplaceAll("fHQel_p","");
      species = AliPID::kElectron;
    } else if(histname.BeginsWith("fHQmu")){ // Muon histogram
      histname.ReplaceAll("fHQmu_p","");
      species = AliPID::kMuon;
    } else if(histname.BeginsWith("fHQpi")){ // Pion histogram
      histname.ReplaceAll("fHQpi_p","");
      species = AliPID::kPion;
    }else if(histname.BeginsWith("fHQka")){ // Kaon histogram
      histname.ReplaceAll("fHQka_p","");
      species = AliPID::kKaon;
    }else if(histname.BeginsWith("fHQpr")){ // Proton histogram
      histname.ReplaceAll("fHQpr_p","");
      species = AliPID::kProton;
    } else continue;
    pbin = histname.Atoi() - 1;
    AliDebug(1, Form("Species %d, Histname %s, Pbin %d, Position in container %d", species, histname.Data(), pbin, arrayPos));
    fMapRefHists[species][pbin] = arrayPos;
    fReferences->AddAt((hnew = new TH1F(*dynamic_cast<TH1F *>(tmp))), arrayPos);
    hnew->SetDirectory(gROOT);
    arrayPos++;
  }
  if(!arrayPos){
    AliError("Failed loading References");
    return kFALSE;
  }
  AliDebug(2, Form("Successfully loaded %d Reference Histograms", arrayPos));
  SetBit(kIsOwner, kTRUE);
  return kTRUE;
}

//____________________________________________________________
Bool_t AliTRDPIDResponse::GetResponse(Int_t n, Double_t *dedx, Float_t *p, Double_t prob[AliPID::kSPECIES], Bool_t kNorm) const
{
  //
  // Calculate TRD likelihood values for the track based on dedx and 
  // momentum values. The likelihoods are calculated by query the 
  // reference data depending on the PID method selected
  //
  // Input parameter :
  //   n - number of dedx slices/chamber
  //   dedx - array of dedx slices organized layer wise
  //   p - array of momentum measurements organized layer wise
  // 
  // Return parameters
  //   prob - probabilities allocated by TRD for particle specis
  //   kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
  // 
  // Return value
  //   true if calculation success
  // 

  if(!fReferences){
    AliWarning("Missing reference data. PID calculation not possible.");
    return kFALSE;
  }

  for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
  Double_t prLayer[AliPID::kSPECIES];
  Double_t DE[10], s(0.);
  for(Int_t il(kNlayer); il--;){
    memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
    if(!CookdEdx(n, &dedx[il*n], &DE[0])) continue;

    s=0.;
    for(Int_t is(AliPID::kSPECIES); is--;){
      if((DE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], DE[0]);
      AliDebug(3, Form("Probability for Species %d in Layer %d: %f", is, il, prLayer[is]));
      s+=prLayer[is];
    }
    if(s<1.e-30){
      AliDebug(2, Form("Null all species prob layer[%d].", il));
      continue;
    }
    for(Int_t is(AliPID::kSPECIES); is--;){
      if(kNorm) prLayer[is] /= s;
      prob[is] *= prLayer[is];
    }
  }
  if(!kNorm) return kTRUE;

  s=0.;
  for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
  if(s<1.e-30){
    AliDebug(2, "Null total prob.");
    return kFALSE;
  }
  for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
  return kTRUE;
}

//____________________________________________________________
Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const {
  //
  // Get the non-normalized probability for a certain particle species as coming
  // from the reference histogram
  // Interpolation between momentum bins
  //
  AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
  Int_t pbin = GetLowerMomentumBin(plocal);  
  AliDebug(1, Form("Bin %d", pbin));
  Double_t pLayer = 0.;
  // Do Interpolation exept for underflow and overflow
  if(pbin >= 0 && pbin < kNPBins-1){
    TH1 *refHistUpper = NULL, *refHistLower = NULL;
    if(fMapRefHists[species][pbin] >= 0)
      refHistLower = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin]));
    if(fMapRefHists[species][pbin+1] >= 0)
      refHistUpper = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin+1]));
    AliDebug(1, Form("Reference Histos (Upper/Lower): [%p|%p]", refHistUpper, refHistLower));
  
    if (refHistLower && refHistUpper ) {
      Double_t pLower = refHistLower->GetBinContent(refHistLower->GetXaxis()->FindBin(dEdx));
      Double_t pUpper = refHistUpper->GetBinContent(refHistUpper->GetXaxis()->FindBin(dEdx));
  
      pLayer = pLower + (pUpper - pLower)/(fgkPBins[pbin+1]-fgkPBins[pbin]) * (plocal - fgkPBins[pbin]); 
    }
  }
  else{
    TH1 *refHist = NULL;
    if(pbin < 0){
      // underflow
      if(fMapRefHists[species][0] >= 0){
        refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][0])); 
      }
    } else {
      // overflow
      if(fMapRefHists[species][kNPBins-1] > 0)
        refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][kNPBins-1])); 
    }
    if (refHist)
      pLayer = refHist->GetBinContent(refHist->GetXaxis()->FindBin(dEdx));
  }
  AliDebug(1, Form("Probability %f", pLayer));
  return pLayer;
}

//____________________________________________________________
Int_t AliTRDPIDResponse::GetLowerMomentumBin(Double_t p) const {
  //
  // Get the momentum bin for a given momentum value
  //
  Int_t bin = -1;
  if(p > fgkPBins[kNPBins-1]) return kNPBins-1;
  for(Int_t ibin = 0; ibin < kNPBins - 1; ibin++){
    if(p >= fgkPBins[ibin] && p < fgkPBins[ibin+1]){
      bin = ibin;
      break;
    }
  }
  return bin;
}

//____________________________________________________________
void AliTRDPIDResponse::SetOwner(){
  //
  // Make Deep Copy of the Reference Histograms
  //
  if(!fReferences || IsOwner()) return;
  TObjArray *ctmp = new TObjArray();
  for(Int_t ien = 0; ien < fReferences->GetEntriesFast(); ien++)
    ctmp->AddAt(fReferences->UncheckedAt(ien)->Clone(), ien);
  delete fReferences;
  fReferences = ctmp;
  SetBit(kIsOwner, kTRUE);
}

//____________________________________________________________
Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, Double_t *in, Double_t *out) const
{
  switch(fPIDmethod){
  case kNN: // NN 
    break;
  case kLQ2D: // 2D LQ 
    break;
  case kLQ1D: // 1D LQ 
    out[0]= 0.;
    for(Int_t islice = 0; islice < nSlice; islice++) out[0] += in[islice] * fGainNormalisationFactor;
    if(out[0] < 1e-6) return kFALSE;
    break;
  default:
    return kFALSE;
  }
  return kTRUE;
}
Commit	Line	Data
ffb1ee30	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	// PID Response class for the TRD detector
	17	// Based on 1D Likelihood approach
	18	// Calculation of probabilities using Bayesian approach
	19	// Attention: This method is only used to separate electrons from pions
	20	//
	21	// Authors:
	22	// Markus Fasel <M.Fasel@gsi.de>
	23	// Anton Andronic <A.Andronic@gsi.de>
	24	//
	25	#include <TClass.h>
	26	#include <TAxis.h>
	27	#include <TFile.h>
	28	#include <TH1.h>
ce5d6908	29	#include <TKey.h>
ffb1ee30	30	#include <TObjArray.h>
ffb1ee30	31	#include <TString.h>
ce5d6908	32	#include <TROOT.h>
ffb1ee30	33	#include <TSystem.h>
2ad33025	34	#include <TDirectory.h>
ffb1ee30	35
	36	#include "AliLog.h"
	37
	38	#include "AliTRDPIDResponse.h"
	39
	40	ClassImp(AliTRDPIDResponse)
	41
	42	const Double_t AliTRDPIDResponse::fgkPBins[kNPBins] = {1, 2, 3, 4, 5, 6};
	43
	44	//____________________________________________________________
e0de37e9	45	AliTRDPIDResponse::AliTRDPIDResponse():
	46	TObject()
	47	,fReferences(NULL)
	48	,fGainNormalisationFactor(1.)
	49	,fPIDmethod(kLQ1D)
ffb1ee30	50	{
e0de37e9	51	//
	52	// Default constructor
	53	//
ce5d6908	54	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++)
	55	for(Int_t ipbin = 0; ipbin < kNPBins; ipbin++)
	56	fMapRefHists[ispec][ipbin] = -1;
ffb1ee30	57	}
	58
	59	//____________________________________________________________
	60	AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
e0de37e9	61	TObject(ref)
	62	,fReferences(ref.fReferences)
	63	,fGainNormalisationFactor(ref.fGainNormalisationFactor)
	64	,fPIDmethod(ref.fPIDmethod)
ffb1ee30	65	{
e0de37e9	66	//
	67	// Copy constructor
	68	// Flat copy of the reference histos
	69	// For creating a deep copy call SetOwner
	70	//
ce5d6908	71	Int_t size = (AliPID::kSPECIES)*(kNPBins);
e0de37e9	72	memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * size);
ffb1ee30	73	}
	74
	75	//____________________________________________________________
	76	AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
e0de37e9	77	//
	78	// Assignment operator
	79	// Performs a flat copy of the reference histos
	80	//
	81	if(this == &ref) return *this;
	82
	83	// Make copy
	84	TObject::operator=(ref);
	85	if(IsOwner()){
	86	if(fReferences){
	87	fReferences->Delete();
	88	delete fReferences;
	89	}
	90	SetBit(kIsOwner, kFALSE);
	91	} else if(fReferences) delete fReferences;
	92	fReferences = ref.fReferences;
	93	Int_t size = (AliPID::kSPECIES)*(kNPBins);
	94	memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * size);
	95	fPIDmethod = ref.fPIDmethod;
	96
	97	return *this;
ffb1ee30	98	}
	99
	100	//____________________________________________________________
	101	AliTRDPIDResponse::~AliTRDPIDResponse(){
e0de37e9	102	//
	103	// Destructor
	104	//
	105	if(IsOwner()){
	106	// Destroy histos
	107	if(fReferences) fReferences->Delete();
	108	}
	109	if(fReferences) delete fReferences;
ffb1ee30	110	}
	111
	112	//____________________________________________________________
	113	Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
e0de37e9	114	//
	115	// Load References into the toolkit
	116	//
	117	AliDebug(1, "Loading reference histos from root file");
2ad33025	118	TDirectory *owd = gDirectory;// store old working directory
e0de37e9	119
	120	if(!filename)
	121	filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
	122	TFile *in = TFile::Open(filename);
	123	if(!in){
	124	AliError("Ref file not available.");
	125	return kFALSE;
	126	}
	127
	128	gROOT->cd();
	129	fReferences = new TObjArray(AliPID::kSPECIES*kNPBins);
	130	fReferences->SetOwner();
	131	TIter iter(in->GetListOfKeys());
	132	TKey *histkey = NULL;
	133	TObject *tmp = NULL;
	134	Int_t arrayPos = 0, pbin = 0;
	135	AliPID::EParticleType species;
	136	TString histname;
	137	TH1 *hnew = NULL;
	138	while((histkey = dynamic_cast<TKey *>(iter()))){
	139	tmp = histkey->ReadObj();
	140	histname = tmp->GetName();
	141	if(histname.BeginsWith("fHQel")){ // Electron histogram
	142	histname.ReplaceAll("fHQel_p","");
	143	species = AliPID::kElectron;
	144	} else if(histname.BeginsWith("fHQmu")){ // Muon histogram
	145	histname.ReplaceAll("fHQmu_p","");
	146	species = AliPID::kMuon;
	147	} else if(histname.BeginsWith("fHQpi")){ // Pion histogram
	148	histname.ReplaceAll("fHQpi_p","");
	149	species = AliPID::kPion;
	150	}else if(histname.BeginsWith("fHQka")){ // Kaon histogram
	151	histname.ReplaceAll("fHQka_p","");
	152	species = AliPID::kKaon;
	153	}else if(histname.BeginsWith("fHQpr")){ // Proton histogram
	154	histname.ReplaceAll("fHQpr_p","");
	155	species = AliPID::kProton;
	156	} else continue;
	157	pbin = histname.Atoi() - 1;
	158	AliDebug(1, Form("Species %d, Histname %s, Pbin %d, Position in container %d", species, histname.Data(), pbin, arrayPos));
	159	fMapRefHists[species][pbin] = arrayPos;
	160	fReferences->AddAt((hnew =new TH1F(dynamic_cast<TH1F >(tmp))), arrayPos);
	161	hnew->SetDirectory(gROOT);
	162	arrayPos++;
	163	}
	164	in->Close();
	165	owd->cd();
	166	AliDebug(2, Form("Successfully loaded %d Reference Histograms", arrayPos));
	167	SetBit(kIsOwner, kTRUE);
	168	delete in;
	169	return kTRUE;
	170	}
2ad33025	171
e0de37e9	172	//____________________________________________________________
	173	Bool_t AliTRDPIDResponse::Load(const TObjArray *histos){
	174	//
	175	// Load Reference into the PID Response (from OCDB/OADB)
	176	//
	177	AliDebug(1, "Setting references (from the database)");
	178	if(fReferences) fReferences->Clear();
	179	else
	180	fReferences = new TObjArray(AliPID::kSPECIES*kNPBins);
	181	fReferences->SetOwner();
	182	TIter iter(histos);
	183	Int_t arrayPos = 0, pbin = 0;
	184	AliPID::EParticleType species;
	185	TObject *tmp = NULL;
	186	TString histname;
	187	TH1 *hnew = NULL;
	188	while((tmp = iter())){
	189	histname = tmp->GetName();
	190	if(histname.BeginsWith("fHQel")){ // Electron histogram
	191	histname.ReplaceAll("fHQel_p","");
	192	species = AliPID::kElectron;
	193	} else if(histname.BeginsWith("fHQmu")){ // Muon histogram
	194	histname.ReplaceAll("fHQmu_p","");
	195	species = AliPID::kMuon;
	196	} else if(histname.BeginsWith("fHQpi")){ // Pion histogram
	197	histname.ReplaceAll("fHQpi_p","");
	198	species = AliPID::kPion;
	199	}else if(histname.BeginsWith("fHQka")){ // Kaon histogram
	200	histname.ReplaceAll("fHQka_p","");
	201	species = AliPID::kKaon;
	202	}else if(histname.BeginsWith("fHQpr")){ // Proton histogram
	203	histname.ReplaceAll("fHQpr_p","");
	204	species = AliPID::kProton;
	205	} else continue;
	206	pbin = histname.Atoi() - 1;
	207	AliDebug(1, Form("Species %d, Histname %s, Pbin %d, Position in container %d", species, histname.Data(), pbin, arrayPos));
	208	fMapRefHists[species][pbin] = arrayPos;
	209	fReferences->AddAt((hnew = new TH1F(dynamic_cast<TH1F >(tmp))), arrayPos);
	210	hnew->SetDirectory(gROOT);
	211	arrayPos++;
	212	}
	213	if(!arrayPos){
	214	AliError("Failed loading References");
	215	return kFALSE;
	216	}
	217	AliDebug(2, Form("Successfully loaded %d Reference Histograms", arrayPos));
	218	SetBit(kIsOwner, kTRUE);
	219	return kTRUE;
ffb1ee30	220	}
	221
	222	//____________________________________________________________
b439f460	223	Bool_t AliTRDPIDResponse::GetResponse(Int_t n, Double_t dedx, Float_t p, Double_t prob[AliPID::kSPECIES], Bool_t kNorm) const
ffb1ee30	224	{
e0de37e9	225	//
	226	// Calculate TRD likelihood values for the track based on dedx and
	227	// momentum values. The likelihoods are calculated by query the
	228	// reference data depending on the PID method selected
	229	//
	230	// Input parameter :
	231	// n - number of dedx slices/chamber
	232	// dedx - array of dedx slices organized layer wise
	233	// p - array of momentum measurements organized layer wise
	234	//
	235	// Return parameters
	236	// prob - probabilities allocated by TRD for particle specis
	237	// kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
	238	//
	239	// Return value
	240	// true if calculation success
	241	//
ffb1ee30	242
e0de37e9	243	if(!fReferences){
	244	AliWarning("Missing reference data. PID calculation not possible.");
	245	return kFALSE;
	246	}
ffb1ee30	247
e0de37e9	248	for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
	249	Double_t prLayer[AliPID::kSPECIES];
	250	Double_t DE[10], s(0.);
	251	for(Int_t il(kNlayer); il--;){
	252	memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
	253	if(!CookdEdx(n, &dedx[il*n], &DE[0])) continue;
ffb1ee30	254
e0de37e9	255	s=0.;
	256	for(Int_t is(AliPID::kSPECIES); is--;){
	257	if((DE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], DE[0]);
	258	AliDebug(3, Form("Probability for Species %d in Layer %d: %f", is, il, prLayer[is]));
	259	s+=prLayer[is];
	260	}
	261	if(s<1.e-30){
	262	AliDebug(2, Form("Null all species prob layer[%d].", il));
	263	continue;
	264	}
	265	for(Int_t is(AliPID::kSPECIES); is--;){
	266	if(kNorm) prLayer[is] /= s;
	267	prob[is] *= prLayer[is];
	268	}
	269	}
	270	if(!kNorm) return kTRUE;
ffb1ee30	271
e0de37e9	272	s=0.;
	273	for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
	274	if(s<1.e-30){
	275	AliDebug(2, "Null total prob.");
	276	return kFALSE;
	277	}
	278	for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
	279	return kTRUE;
ffb1ee30	280	}
ffb1ee30	281
ffb1ee30	282	//____________________________________________________________
b439f460	283	Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const {
e0de37e9	284	//
	285	// Get the non-normalized probability for a certain particle species as coming
	286	// from the reference histogram
	287	// Interpolation between momentum bins
	288	//
	289	AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
	290	Int_t pbin = GetLowerMomentumBin(plocal);
	291	AliDebug(1, Form("Bin %d", pbin));
	292	Double_t pLayer = 0.;
	293	// Do Interpolation exept for underflow and overflow
	294	if(pbin >= 0 && pbin < kNPBins-1){
	295	TH1 refHistUpper = NULL, refHistLower = NULL;
	296	if(fMapRefHists[species][pbin] >= 0)
	297	refHistLower = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin]));
	298	if(fMapRefHists[species][pbin+1] >= 0)
	299	refHistUpper = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin+1]));
	300	AliDebug(1, Form("Reference Histos (Upper/Lower): [%p\|%p]", refHistUpper, refHistLower));
	301
	302	if (refHistLower && refHistUpper ) {
	303	Double_t pLower = refHistLower->GetBinContent(refHistLower->GetXaxis()->FindBin(dEdx));
	304	Double_t pUpper = refHistUpper->GetBinContent(refHistUpper->GetXaxis()->FindBin(dEdx));
	305
	306	pLayer = pLower + (pUpper - pLower)/(fgkPBins[pbin+1]-fgkPBins[pbin]) * (plocal - fgkPBins[pbin]);
	307	}
	308	}
	309	else{
	310	TH1 *refHist = NULL;
	311	if(pbin < 0){
	312	// underflow
	313	if(fMapRefHists[species][0] >= 0){
	314	refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][0]));
	315	}
	316	} else {
	317	// overflow
	318	if(fMapRefHists[species][kNPBins-1] > 0)
	319	refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][kNPBins-1]));
	320	}
	321	if (refHist)
	322	pLayer = refHist->GetBinContent(refHist->GetXaxis()->FindBin(dEdx));
	323	}
	324	AliDebug(1, Form("Probability %f", pLayer));
	325	return pLayer;
ffb1ee30	326	}
	327
	328	//____________________________________________________________
b439f460	329	Int_t AliTRDPIDResponse::GetLowerMomentumBin(Double_t p) const {
e0de37e9	330	//
	331	// Get the momentum bin for a given momentum value
	332	//
	333	Int_t bin = -1;
	334	if(p > fgkPBins[kNPBins-1]) return kNPBins-1;
	335	for(Int_t ibin = 0; ibin < kNPBins - 1; ibin++){
	336	if(p >= fgkPBins[ibin] && p < fgkPBins[ibin+1]){
	337	bin = ibin;
	338	break;
	339	}
	340	}
	341	return bin;
ffb1ee30	342	}
	343
	344	//____________________________________________________________
	345	void AliTRDPIDResponse::SetOwner(){
e0de37e9	346	//
	347	// Make Deep Copy of the Reference Histograms
	348	//
	349	if(!fReferences \|\| IsOwner()) return;
	350	TObjArray *ctmp = new TObjArray();
	351	for(Int_t ien = 0; ien < fReferences->GetEntriesFast(); ien++)
	352	ctmp->AddAt(fReferences->UncheckedAt(ien)->Clone(), ien);
	353	delete fReferences;
	354	fReferences = ctmp;
	355	SetBit(kIsOwner, kTRUE);
ffb1ee30	356	}
	357
	358	//____________________________________________________________
e0de37e9	359	Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, Double_t in, Double_t out) const
ffb1ee30	360	{
e0de37e9	361	switch(fPIDmethod){
	362	case kNN: // NN
	363	break;
	364	case kLQ2D: // 2D LQ
	365	break;
	366	case kLQ1D: // 1D LQ
	367	out[0]= 0.;
	368	for(Int_t islice = 0; islice < nSlice; islice++) out[0] += in[islice] * fGainNormalisationFactor;
	369	if(out[0] < 1e-6) return kFALSE;
	370	break;
	371	default:
	372	return kFALSE;
	373	}
	374	return kTRUE;
ffb1ee30	375	}
ffb1ee30	376