[u/mrichter/AliRoot.git] / STEER / STEERBase / 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 <TAxis.h>
#include <TClass.h>
#include <TDirectory.h>
#include <TFile.h>
#include <TH1.h>
#include <TKey.h>
#include <TMath.h>
#include <TObjArray.h>
#include <TROOT.h> 
#include <TString.h>
#include <TSystem.h>

#include "AliLog.h"
 
#include "AliVTrack.h"

#include "AliTRDPIDResponseObject.h"
#include "AliTRDPIDResponse.h"
//#include "AliTRDTKDInterpolator.h"
#include "AliTRDNDFast.h"

ClassImp(AliTRDPIDResponse)

//____________________________________________________________
AliTRDPIDResponse::AliTRDPIDResponse():
  TObject()
  ,fkPIDResponseObject(NULL)
  ,fGainNormalisationFactor(1.)
{
  //
  // Default constructor
  //
}

//____________________________________________________________
AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
  TObject(ref)
  ,fkPIDResponseObject(NULL)
  ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
{
  //
  // Copy constructor
  //
}

//____________________________________________________________
AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
  //
  // Assignment operator
  //
  if(this == &ref) return *this;
  
  // Make copy
  TObject::operator=(ref);
  fGainNormalisationFactor = ref.fGainNormalisationFactor;
  fkPIDResponseObject = ref.fkPIDResponseObject;
  
  return *this;
}

//____________________________________________________________
AliTRDPIDResponse::~AliTRDPIDResponse(){
  //
  // Destructor
  //
    if(IsOwner()) delete fkPIDResponseObject;
}

//____________________________________________________________
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();
  fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(in->Get("TRDPIDResponse")->Clone());
  in->Close(); delete in;
  owd->cd();
  SetBit(kIsOwner, kTRUE);
  AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDResponseObject->GetNumberOfMomentumBins()));
  return kTRUE;
}


//____________________________________________________________
Double_t AliTRDPIDResponse::GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType type) const
{
  //
  //calculate the TRD nSigma
  //

  const Double_t badval = -9999;
  Double_t info[5]; for(int i=0; i<5; i++){info[i]=badval;}

  const Double_t delta = GetSignalDelta(track, type, kFALSE, info);

  const Double_t mean = info[0];
  const Double_t res = info[1];
  if(res<0){
    return badval;
  }

  const Double_t sigma = mean*res;
  const Double_t eps = 1e-12;
  return delta/(sigma + eps);
}

//____________________________________________________________
Double_t AliTRDPIDResponse::GetSignalDelta( const AliVTrack* track, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/, Double_t *info/*=0x0*/) const
{
  //
  //calculate the TRD signal difference w.r.t. the expected
  //output other information in info[]
  //

  //============================================================================================>>>>>>>>>>>>
  //Data production dependent, in the end should go to OADB
  //============================================================================================
  //~/.task/TRDdEdx/2013output/resolution$ r bgmean.root
  // hfit;1  1.521e+00 , 7.294e-01 , 8.066e+00 , 2.146e-02 , 3.137e+01 , 7.160e-15 , 2.838e+00 , 1.100e+01 , (14290.45/67) conv_runlist_LHC13b.pass3.root trdncls/trdnch>=18 && trdnch>=6 && tpcncls>120 && pidV0>=1          
  // hscn0_wid;1     -4.122e-01 , 1.019e+00 , -1.319e-01 , (15071.70/120) conv_runlist_LHC13b.pass3.root tpcncls>120 && pidV0>=1

  //mean signal parametrization
  Double_t meanpars[]={1.521e+00 , 7.294e-01 , 8.066e+00 , 2.146e-02 , 3.137e+01 , 7.160e-15 , 2.838e+00 , 1.100e+01};
  const Int_t nmp = sizeof(meanpars)/sizeof(Double_t);

  if(type==AliPID::kProton){
    const Double_t tmpproton[]={2.026e+00 , -1.462e-04 , 1.202e+00 , 1.162e-01 , 2.092e+00 , -3.018e-02 , 3.665e+00 , 3.487e-07};
    for(Int_t ii=0; ii<nmp; ii++){
      meanpars[ii] = tmpproton[ii];
    }
  }
  else if(type==AliPID::kPion || type==AliPID::kElectron){
    const Double_t tmppe[]={1.508e+00 , 7.356e-01 , 8.002e+00 , 1.932e-02 , 3.058e+01 , 5.114e-18 , 3.561e+00 , 1.408e+01};
    for(Int_t ii=0; ii<nmp; ii++){
      meanpars[ii] = tmppe[ii];
    }
  }

  //resolution parametrization
  const Double_t respars[]={-4.122e-01 , 1.019e+00 , -1.319e-01};

  //cut on number of chambers
  const Double_t cutch = 6;

  //cut on mean number of clusters per chamber
  const Double_t cutclsperch = 18;
  //============================================================================================<<<<<<<<<<<<<


  const Double_t badval = -9999;

  const Double_t nch = track->GetTRDNchamberdEdx();
  if(nch < cutch){
    return badval;
  }

  const Double_t ncls = track->GetTRDNclusterdEdx();
  if( ncls/nch < cutclsperch){
    return badval;
  }

  Double_t pTRD = -999;
  for(Int_t ich=0; ich<6; ich++){
    pTRD = track->GetTRDmomentum(ich);
    if(pTRD>0)
      break;
  }
  if(pTRD<0){
    return badval;
  }

  const Double_t bg = pTRD/AliPID::ParticleMass(type);
  const Double_t expsig = MeandEdxTR(&bg, meanpars);

  if(info){
    info[0]= expsig;
    info[1]= ResolutiondEdxTR(&ncls, respars);
  }

  const Double_t eps = 1e-10;

  if(ratio){
    return track->GetTRDsignal()/(expsig +  eps);
  }
  else{
    return track->GetTRDsignal() - expsig;
  }
}


Double_t AliTRDPIDResponse::ResolutiondEdxTR(const Double_t * xx,  const Double_t * par)
{
  //
  //resolution 
  //npar=3
  //

  const Double_t ncls = xx[0];

  return par[0]+par[1]*TMath::Power(ncls, par[2]);
}

Double_t AliTRDPIDResponse::MeandEdxTR(const Double_t * xx,  const Double_t * pin)
{
  //
  //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
  //npar = 8 = 3+5
  //

  Double_t ptr[4]={0};
  for(int ii=0; ii<3; ii++){
    ptr[ii+1]=pin[ii];
  }
  return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
}

Double_t AliTRDPIDResponse::MeanTR(const Double_t * xx,  const Double_t * par)
{
  //
  //LOGISTIC parametrization for TR, in unit of MIP
  //npar = 4
  //
                                                                                                                                                                                                                                                        
  const Double_t bg = xx[0];
  const Double_t gamma = sqrt(1+bg*bg);

  const Double_t p0 = TMath::Abs(par[1]);
  const Double_t p1 = TMath::Abs(par[2]);
  const Double_t p2 = TMath::Abs(par[3]);

  const Double_t zz = TMath::Log(gamma);
  const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );

  return par[0]+tryield;
}

Double_t AliTRDPIDResponse::MeandEdx(const Double_t * xx,  const Double_t * par)
{
  //
  //ALEPH parametrization for dEdx
  //npar = 5
  //

  const Double_t bg = xx[0];
  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);

  const Double_t p0 = TMath::Abs(par[0]);
  const Double_t p1 = TMath::Abs(par[1]);

  const Double_t p2 = TMath::Abs(par[2]);

  const Double_t p3 = TMath::Abs(par[3]);
  const Double_t p4 = TMath::Abs(par[4]);

  const Double_t aa = TMath::Power(beta, p3);

  const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );

  return (p1-aa-bb)*p0/aa;

}


//____________________________________________________________
Int_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES],ETRDPIDMethod PIDmethod,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
  //   number of tracklets used for PID, 0 if no PID
    //
    AliDebug(3,Form(" Response for PID method: %d",PIDmethod));


    if(!fkPIDResponseObject){
	AliDebug(3,"Missing reference data. PID calculation not possible.");
	return 0;
    }

    for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
    Double_t prLayer[AliPID::kSPECIES];
    Double_t dE[10], s(0.);
    Int_t ntrackletsPID=0;
    for(Int_t il(kNlayer); il--;){
	memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
	if(!CookdEdx(n, &dedx[il*n], &dE[0],PIDmethod)) continue;
        s=0.;
        Bool_t filled=kTRUE;
	for(Int_t is(AliPID::kSPECIES); is--;){
	    //if((PIDmethod==kLQ2D)&&(!(is==0||is==2)))continue;
	    if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], &dE[0],PIDmethod);
	    AliDebug(3, Form("Probability for Species %d in Layer %d: %e", is, il, prLayer[is]));
	    if(prLayer[is]<1.e-30){
		AliDebug(2, Form("Null for species %d species prob layer[%d].",is,il));
		filled=kFALSE;
		break;
	    }
	    s+=prLayer[is];
	}
	if(!filled){
	    continue;
	}
	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];
	}
	ntrackletsPID++;
    }
    if(!kNorm) return ntrackletsPID;

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

//____________________________________________________________
Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx,ETRDPIDMethod PIDmethod) 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));

  Double_t probLayer = 0.;

  Float_t pLower, pUpper;
	
  AliTRDNDFast *refUpper = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,PIDmethod)),
      *refLower = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,PIDmethod));

  // Do Interpolation exept for underflow and overflow
  if(refLower && refUpper){
      Double_t probLower = refLower->Eval(dEdx);
      Double_t probUpper = refUpper->Eval(dEdx);

      probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
  } else if(refLower){
     // underflow
      probLayer = refLower->Eval(dEdx);
  } else if(refUpper){
      // overflow
      probLayer = refUpper->Eval(dEdx);
  } else {
      AliError("No references available");
  }
  AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));

  return probLayer;

/* old implementation

switch(PIDmethod){
case kNN: // NN
      break;
  case kLQ2D: // 2D LQ
      {
	  if(species==0||species==2){ // references only for electrons and pions
	      Double_t error = 0.;
	      Double_t point[kNslicesLQ2D];
	      for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}

	      AliTRDTKDInterpolator *refUpper = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ2D)),
	      *refLower = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));
	      // Do Interpolation exept for underflow and overflow
	      if(refLower && refUpper){
                  Double_t probLower=0,probUpper=0;
		  refLower->Eval(point,probLower,error);
                  refUpper->Eval(point,probUpper,error);
		  probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
	      } else if(refLower){
		  // underflow
		  refLower->Eval(point,probLayer,error);
		  } else if(refUpper){
		  // overflow
		  refUpper->Eval(point,probLayer,error);
	      } else {
		  AliError("No references available");
	      }
	      AliDebug(2,Form("Eval 2D Q0 %f Q1 %f P %e Err %e",point[0],point[1],probLayer,error));
	  }
      }
      break;
  case kLQ1D: // 1D LQ
      {
	  TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ1D)),
	      *refLower = dynamic_cast<TH1 *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ1D));
	  // Do Interpolation exept for underflow and overflow
	  if(refLower && refUpper){
	      Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
	      Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));

	      probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
	  } else if(refLower){
	      // underflow
	      probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
	  } else if(refUpper){
	      // overflow
	      probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
	  } else {
	      AliError("No references available");
	  }
	  AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
      }
      break;
  default:
      break;
      }
      return probLayer;
      */

}

//____________________________________________________________
void AliTRDPIDResponse::SetOwner(){
  //
  // Make Deep Copy of the Reference Histograms
  //
    if(!fkPIDResponseObject || IsOwner()) return;
    const AliTRDPIDResponseObject *tmp = fkPIDResponseObject;
    fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(tmp->Clone());
    SetBit(kIsOwner, kTRUE);
}

//____________________________________________________________
Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out,ETRDPIDMethod PIDmethod) const
{
    //
    // Recalculate dE/dx
    //
  switch(PIDmethod){
  case kNN: // NN 
      break;
  case kLQ2D: // 2D LQ
      out[0]=0;
      out[1]=0;
      for(Int_t islice = 0; islice < nSlice; islice++){
	  if(in[islice]<=0){out[0]=0;out[1]=0;return kFALSE;}  // Require that all slices are filled
	  if(islice<fkPIDResponseObject->GetNSlicesQ0())out[0]+= in[islice];
	  else out[1]+= in[islice];
      }
      // normalize signal to number of slices
      out[0]*=1./Double_t(fkPIDResponseObject->GetNSlicesQ0());
      out[1]*=1./Double_t(nSlice-fkPIDResponseObject->GetNSlicesQ0());
      if(out[0] < 1e-6) return kFALSE;
      AliDebug(3,Form("CookdEdx Q0 %f Q1 %f",out[0],out[1]));
      break;
  case kLQ1D: // 1D LQ
      out[0]= 0.;
      for(Int_t islice = 0; islice < nSlice; islice++)
	  if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor;   // Protect against negative values for slices having no dE/dx information
      out[0]*=1./Double_t(nSlice);
      if(out[0] < 1e-6) return kFALSE;
      AliDebug(3,Form("CookdEdx dEdx %f",out[0]));
      break;

  default:
      return kFALSE;
  }
  return kTRUE;
}

//____________________________________________________________
Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level,Double_t centrality,ETRDPIDMethod PIDmethod) const {
    //
    // Check whether particle is identified as electron assuming a certain electron efficiency level
    // Only electron and pion hypothesis is taken into account
    //
    // Inputs:
    //         Number of tracklets
    //         Likelihood values
    //         Momentum
    //         Electron efficiency level
    //
    // If the function fails when the params are not accessible, the function returns true
    //
  if(!fkPIDResponseObject){
    AliDebug(3,"No PID Param object available");
    return kTRUE;
  } 
  Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
  Double_t params[4];
  if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,centrality,PIDmethod)){
    AliError("No Params found for the given configuration");
    return kTRUE;
  }
  Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
  if(probEle > TMath::Max(TMath::Min(threshold, 0.99), 0.2)) return kTRUE; // truncate the threshold upperwards to 0.999 and lowerwards to 0.2 and exclude unphysical values
  return kFALSE;
}

//____________________________________________________________
Bool_t AliTRDPIDResponse::SetPIDResponseObject(const AliTRDPIDResponseObject * obj){

    fkPIDResponseObject = obj;
    if((AliLog::GetDebugLevel("",IsA()->GetName()))>0)fkPIDResponseObject->Print("");
    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	//
ffb1ee30	25	#include <TAxis.h>
ea235c90	26	#include <TClass.h>
ea235c90	27	#include <TDirectory.h>
ffb1ee30	28	#include <TFile.h>
ffb1ee30	29	#include <TH1.h>
ce5d6908	30	#include <TKey.h>
ea235c90	31	#include <TMath.h>
ffb1ee30	32	#include <TObjArray.h>
ce5d6908	33	#include <TROOT.h>
ea235c90	34	#include <TString.h>
ffb1ee30	35	#include <TSystem.h>
	36
	37	#include "AliLog.h"
5cd0300d	38
5cd0300d	39	#include "AliVTrack.h"
ffb1ee30	40
db0e2c5f	41	#include "AliTRDPIDResponseObject.h"
ffb1ee30	42	#include "AliTRDPIDResponse.h"
239fe91c	43	//#include "AliTRDTKDInterpolator.h"
239fe91c	44	#include "AliTRDNDFast.h"
ffb1ee30	45
	46	ClassImp(AliTRDPIDResponse)
	47
ffb1ee30	48	//____________________________________________________________
e0de37e9	49	AliTRDPIDResponse::AliTRDPIDResponse():
e0de37e9	50	TObject()
db0e2c5f	51	,fkPIDResponseObject(NULL)
e0de37e9	52	,fGainNormalisationFactor(1.)
ffb1ee30	53	{
e0de37e9	54	//
	55	// Default constructor
	56	//
ffb1ee30	57	}
	58
	59	//____________________________________________________________
	60	AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
e0de37e9	61	TObject(ref)
db0e2c5f	62	,fkPIDResponseObject(NULL)
e0de37e9	63	,fGainNormalisationFactor(ref.fGainNormalisationFactor)
ffb1ee30	64	{
e0de37e9	65	//
e0de37e9	66	// Copy constructor
e0de37e9	67	//
ffb1ee30	68	}
	69
	70	//____________________________________________________________
	71	AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
e0de37e9	72	//
e0de37e9	73	// Assignment operator
e0de37e9	74	//
	75	if(this == &ref) return *this;
	76
	77	// Make copy
	78	TObject::operator=(ref);
51a0ce25	79	fGainNormalisationFactor = ref.fGainNormalisationFactor;
db0e2c5f	80	fkPIDResponseObject = ref.fkPIDResponseObject;
e0de37e9	81
e0de37e9	82	return *this;
ffb1ee30	83	}
	84
	85	//____________________________________________________________
	86	AliTRDPIDResponse::~AliTRDPIDResponse(){
e0de37e9	87	//
	88	// Destructor
	89	//
db0e2c5f	90	if(IsOwner()) delete fkPIDResponseObject;
ffb1ee30	91	}
	92
	93	//____________________________________________________________
db0e2c5f	94	Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
e0de37e9	95	//
	96	// Load References into the toolkit
	97	//
	98	AliDebug(1, "Loading reference histos from root file");
2ad33025	99	TDirectory *owd = gDirectory;// store old working directory
e0de37e9	100
	101	if(!filename)
	102	filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
	103	TFile *in = TFile::Open(filename);
	104	if(!in){
	105	AliError("Ref file not available.");
	106	return kFALSE;
	107	}
	108
	109	gROOT->cd();
db0e2c5f	110	fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(in->Get("TRDPIDResponse")->Clone());
51a0ce25	111	in->Close(); delete in;
e0de37e9	112	owd->cd();
e0de37e9	113	SetBit(kIsOwner, kTRUE);
db0e2c5f	114	AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDResponseObject->GetNumberOfMomentumBins()));
e0de37e9	115	return kTRUE;
ffb1ee30	116	}
ffb1ee30	117
5cd0300d	118
	119
	120	//____________________________________________________________
	121	Double_t AliTRDPIDResponse::GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType type) const
	122	{
	123	//
	124	//calculate the TRD nSigma
	125	//
	126
	127	const Double_t badval = -9999;
	128	Double_t info[5]; for(int i=0; i<5; i++){info[i]=badval;}
	129
	130	const Double_t delta = GetSignalDelta(track, type, kFALSE, info);
	131
	132	const Double_t mean = info[0];
	133	const Double_t res = info[1];
	134	if(res<0){
	135	return badval;
	136	}
	137
	138	const Double_t sigma = mean*res;
bebe75fc	139	const Double_t eps = 1e-12;
bebe75fc	140	return delta/(sigma + eps);
5cd0300d	141	}
	142
	143	//____________________________________________________________
	144	Double_t AliTRDPIDResponse::GetSignalDelta( const AliVTrack* track, AliPID::EParticleType type, Bool_t ratio/=kFALSE/, Double_t info/=0x0*/) const
	145	{
	146	//
	147	//calculate the TRD signal difference w.r.t. the expected
	148	//output other information in info[]
	149	//
	150
	151	//============================================================================================>>>>>>>>>>>>
	152	//Data production dependent, in the end should go to OADB
	153	//============================================================================================
	154	//~/.task/TRDdEdx/2013output/resolution$ r bgmean.root
	155	// hfit;1 1.521e+00 , 7.294e-01 , 8.066e+00 , 2.146e-02 , 3.137e+01 , 7.160e-15 , 2.838e+00 , 1.100e+01 , (14290.45/67) conv_runlist_LHC13b.pass3.root trdncls/trdnch>=18 && trdnch>=6 && tpcncls>120 && pidV0>=1
	156	// hscn0_wid;1 -4.122e-01 , 1.019e+00 , -1.319e-01 , (15071.70/120) conv_runlist_LHC13b.pass3.root tpcncls>120 && pidV0>=1
	157
	158	//mean signal parametrization
	159	Double_t meanpars[]={1.521e+00 , 7.294e-01 , 8.066e+00 , 2.146e-02 , 3.137e+01 , 7.160e-15 , 2.838e+00 , 1.100e+01};
	160	const Int_t nmp = sizeof(meanpars)/sizeof(Double_t);
	161
	162	if(type==AliPID::kProton){
	163	const Double_t tmpproton[]={2.026e+00 , -1.462e-04 , 1.202e+00 , 1.162e-01 , 2.092e+00 , -3.018e-02 , 3.665e+00 , 3.487e-07};
	164	for(Int_t ii=0; ii<nmp; ii++){
	165	meanpars[ii] = tmpproton[ii];
	166	}
	167	}
	168	else if(type==AliPID::kPion \|\| type==AliPID::kElectron){
	169	const Double_t tmppe[]={1.508e+00 , 7.356e-01 , 8.002e+00 , 1.932e-02 , 3.058e+01 , 5.114e-18 , 3.561e+00 , 1.408e+01};
	170	for(Int_t ii=0; ii<nmp; ii++){
	171	meanpars[ii] = tmppe[ii];
	172	}
	173	}
	174
	175	//resolution parametrization
	176	const Double_t respars[]={-4.122e-01 , 1.019e+00 , -1.319e-01};
	177
	178	//cut on number of chambers
	179	const Double_t cutch = 6;
	180
	181	//cut on mean number of clusters per chamber
	182	const Double_t cutclsperch = 18;
	183	//============================================================================================<<<<<<<<<<<<<
	184
	185
	186	const Double_t badval = -9999;
	187
	188	const Double_t nch = track->GetTRDNchamberdEdx();
	189	if(nch < cutch){
	190	return badval;
	191	}
	192
	193	const Double_t ncls = track->GetTRDNclusterdEdx();
	194	if( ncls/nch < cutclsperch){
	195	return badval;
	196	}
	197
	198	Double_t pTRD = -999;
	199	for(Int_t ich=0; ich<6; ich++){
	200	pTRD = track->GetTRDmomentum(ich);
	201	if(pTRD>0)
	202	break;
	203	}
	204	if(pTRD<0){
205	return badval;
206	}
207
208	const Double_t bg = pTRD/AliPID::ParticleMass(type);
209	const Double_t expsig = MeandEdxTR(&bg, meanpars);
210
211	if(info){
212	info[0]= expsig;
213	info[1]= ResolutiondEdxTR(&ncls, respars);
214	}
215
bebe75fc	216	const Double_t eps = 1e-10;
bebe75fc	217
5cd0300d	218	if(ratio){
bebe75fc	219	return track->GetTRDsignal()/(expsig + eps);
5cd0300d	220	}
	221	else{
	222	return track->GetTRDsignal() - expsig;
	223	}
	224	}
	225
	226
	227	Double_t AliTRDPIDResponse::ResolutiondEdxTR(const Double_t * xx, const Double_t * par)
	228	{
	229	//
	230	//resolution
	231	//npar=3
	232	//
	233
	234	const Double_t ncls = xx[0];
	235
	236	return par[0]+par[1]*TMath::Power(ncls, par[2]);
	237	}
	238
	239	Double_t AliTRDPIDResponse::MeandEdxTR(const Double_t * xx, const Double_t * pin)
	240	{
	241	//
	242	//ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
	243	//npar = 8 = 3+5
	244	//
	245
	246	Double_t ptr[4]={0};
	247	for(int ii=0; ii<3; ii++){
	248	ptr[ii+1]=pin[ii];
	249	}
	250	return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
	251	}
	252
	253	Double_t AliTRDPIDResponse::MeanTR(const Double_t * xx, const Double_t * par)
	254	{
	255	//
	256	//LOGISTIC parametrization for TR, in unit of MIP
	257	//npar = 4
	258	//
	259
	260	const Double_t bg = xx[0];
	261	const Double_t gamma = sqrt(1+bg*bg);
	262
	263	const Double_t p0 = TMath::Abs(par[1]);
	264	const Double_t p1 = TMath::Abs(par[2]);
	265	const Double_t p2 = TMath::Abs(par[3]);
	266
	267	const Double_t zz = TMath::Log(gamma);
	268	const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
	269
	270	return par[0]+tryield;
	271	}
	272
	273	Double_t AliTRDPIDResponse::MeandEdx(const Double_t * xx, const Double_t * par)
	274	{
	275	//
	276	//ALEPH parametrization for dEdx
	277	//npar = 5
	278	//
	279
	280	const Double_t bg = xx[0];
	281	const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
	282
	283	const Double_t p0 = TMath::Abs(par[0]);
284	const Double_t p1 = TMath::Abs(par[1]);
285
286	const Double_t p2 = TMath::Abs(par[2]);
287
288	const Double_t p3 = TMath::Abs(par[3]);
289	const Double_t p4 = TMath::Abs(par[4]);
290
291	const Double_t aa = TMath::Power(beta, p3);
292
293	const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
294
295	return (p1-aa-bb)*p0/aa;
296
297	}
298
299
ffb1ee30	300	//____________________________________________________________
239fe91c	301	Int_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES],ETRDPIDMethod PIDmethod,Bool_t kNorm) const
ffb1ee30	302	{
e0de37e9	303	//
	304	// Calculate TRD likelihood values for the track based on dedx and
	305	// momentum values. The likelihoods are calculated by query the
	306	// reference data depending on the PID method selected
	307	//
	308	// Input parameter :
	309	// n - number of dedx slices/chamber
	310	// dedx - array of dedx slices organized layer wise
	311	// p - array of momentum measurements organized layer wise
	312	//
	313	// Return parameters
	314	// prob - probabilities allocated by TRD for particle specis
	315	// kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
	316	//
	317	// Return value
239fe91c	318	// number of tracklets used for PID, 0 if no PID
bd58d4b9	319	//
	320	AliDebug(3,Form(" Response for PID method: %d",PIDmethod));
	321
ffb1ee30	322
db0e2c5f	323	if(!fkPIDResponseObject){
31746593	324	AliDebug(3,"Missing reference data. PID calculation not possible.");
239fe91c	325	return 0;
db0e2c5f	326	}
ffb1ee30	327
db0e2c5f	328	for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
	329	Double_t prLayer[AliPID::kSPECIES];
	330	Double_t dE[10], s(0.);
239fe91c	331	Int_t ntrackletsPID=0;
db0e2c5f	332	for(Int_t il(kNlayer); il--;){
db0e2c5f	333	memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
bd58d4b9	334	if(!CookdEdx(n, &dedx[il*n], &dE[0],PIDmethod)) continue;
239fe91c	335	s=0.;
db0e2c5f	336	Bool_t filled=kTRUE;
db0e2c5f	337	for(Int_t is(AliPID::kSPECIES); is--;){
239fe91c	338	//if((PIDmethod==kLQ2D)&&(!(is==0\|\|is==2)))continue;
bd58d4b9	339	if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], &dE[0],PIDmethod);
f2762b1c	340	AliDebug(3, Form("Probability for Species %d in Layer %d: %e", is, il, prLayer[is]));
db0e2c5f	341	if(prLayer[is]<1.e-30){
	342	AliDebug(2, Form("Null for species %d species prob layer[%d].",is,il));
	343	filled=kFALSE;
	344	break;
	345	}
	346	s+=prLayer[is];
	347	}
	348	if(!filled){
	349	continue;
	350	}
	351	if(s<1.e-30){
	352	AliDebug(2, Form("Null all species prob layer[%d].", il));
	353	continue;
	354	}
	355	for(Int_t is(AliPID::kSPECIES); is--;){
	356	if(kNorm) prLayer[is] /= s;
	357	prob[is] *= prLayer[is];
	358	}
239fe91c	359	ntrackletsPID++;
e0de37e9	360	}
239fe91c	361	if(!kNorm) return ntrackletsPID;
db0e2c5f	362
	363	s=0.;
	364	for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
e0de37e9	365	if(s<1.e-30){
db0e2c5f	366	AliDebug(2, "Null total prob.");
239fe91c	367	return 0;
e0de37e9	368	}
db0e2c5f	369	for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
239fe91c	370	return ntrackletsPID;
ffb1ee30	371	}
ffb1ee30	372
ffb1ee30	373	//____________________________________________________________
bd58d4b9	374	Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx,ETRDPIDMethod PIDmethod) const {
e0de37e9	375	//
	376	// Get the non-normalized probability for a certain particle species as coming
	377	// from the reference histogram
	378	// Interpolation between momentum bins
	379	//
	380	AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
db0e2c5f	381
51a0ce25	382	Double_t probLayer = 0.;
239fe91c	383
db0e2c5f	384	Float_t pLower, pUpper;
db0e2c5f	385
239fe91c	386	AliTRDNDFast refUpper = dynamic_cast<AliTRDNDFast >(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,PIDmethod)),
	387	refLower = dynamic_cast<AliTRDNDFast >(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,PIDmethod));
	388
	389	// Do Interpolation exept for underflow and overflow
	390	if(refLower && refUpper){
	391	Double_t probLower = refLower->Eval(dEdx);
	392	Double_t probUpper = refUpper->Eval(dEdx);
	393
	394	probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
	395	} else if(refLower){
	396	// underflow
	397	probLayer = refLower->Eval(dEdx);
	398	} else if(refUpper){
	399	// overflow
	400	probLayer = refUpper->Eval(dEdx);
	401	} else {
	402	AliError("No references available");
	403	}
	404	AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
	405
	406	return probLayer;
	407
	408	/* old implementation
	409
	410	switch(PIDmethod){
	411	case kNN: // NN
db0e2c5f	412	break;
	413	case kLQ2D: // 2D LQ
	414	{
f2762b1c	415	if(species==0\|\|species==2){ // references only for electrons and pions
2285168d	416	Double_t error = 0.;
f2762b1c	417	Double_t point[kNslicesLQ2D];
f2762b1c	418	for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}
db0e2c5f	419
239fe91c	420	AliTRDTKDInterpolator refUpper = dynamic_cast<AliTRDTKDInterpolator >(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ2D)),
	421	refLower = dynamic_cast<AliTRDTKDInterpolator >(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));
	422	// Do Interpolation exept for underflow and overflow
	423	if(refLower && refUpper){
	424	Double_t probLower=0,probUpper=0;
	425	refLower->Eval(point,probLower,error);
	426	refUpper->Eval(point,probUpper,error);
	427	probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
	428	} else if(refLower){
	429	// underflow
f2762b1c	430	refLower->Eval(point,probLayer,error);
239fe91c	431	} else if(refUpper){
	432	// overflow
	433	refUpper->Eval(point,probLayer,error);
	434	} else {
f2762b1c	435	AliError("No references available");
	436	}
	437	AliDebug(2,Form("Eval 2D Q0 %f Q1 %f P %e Err %e",point[0],point[1],probLayer,error));
db0e2c5f	438	}
	439	}
	440	break;
	441	case kLQ1D: // 1D LQ
	442	{
	443	TH1 refUpper = dynamic_cast<TH1 >(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ1D)),
	444	refLower = dynamic_cast<TH1 >(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ1D));
	445	// Do Interpolation exept for underflow and overflow
	446	if(refLower && refUpper){
	447	Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
	448	Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
	449
	450	probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
	451	} else if(refLower){
	452	// underflow
	453	probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
	454	} else if(refUpper){
	455	// overflow
	456	probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
	457	} else {
	458	AliError("No references available");
	459	}
f2762b1c	460	AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
db0e2c5f	461	}
	462	break;
	463	default:
	464	break;
239fe91c	465	}
	466	return probLayer;
	467	*/
	468
ffb1ee30	469	}
	470
	471	//____________________________________________________________
	472	void AliTRDPIDResponse::SetOwner(){
e0de37e9	473	//
	474	// Make Deep Copy of the Reference Histograms
	475	//
db0e2c5f	476	if(!fkPIDResponseObject \|\| IsOwner()) return;
	477	const AliTRDPIDResponseObject *tmp = fkPIDResponseObject;
	478	fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(tmp->Clone());
	479	SetBit(kIsOwner, kTRUE);
ffb1ee30	480	}
	481
	482	//____________________________________________________________
bd58d4b9	483	Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out,ETRDPIDMethod PIDmethod) const
ffb1ee30	484	{
f2762b1c	485	//
	486	// Recalculate dE/dx
	487	//
bd58d4b9	488	switch(PIDmethod){
e0de37e9	489	case kNN: // NN
db0e2c5f	490	break;
	491	case kLQ2D: // 2D LQ
	492	out[0]=0;
	493	out[1]=0;
	494	for(Int_t islice = 0; islice < nSlice; islice++){
239fe91c	495	if(in[islice]<=0){out[0]=0;out[1]=0;return kFALSE;} // Require that all slices are filled
db0e2c5f	496	if(islice<fkPIDResponseObject->GetNSlicesQ0())out[0]+= in[islice];
	497	else out[1]+= in[islice];
	498	}
239fe91c	499	// normalize signal to number of slices
	500	out[0]*=1./Double_t(fkPIDResponseObject->GetNSlicesQ0());
	501	out[1]*=1./Double_t(nSlice-fkPIDResponseObject->GetNSlicesQ0());
f2762b1c	502	if(out[0] < 1e-6) return kFALSE;
f2762b1c	503	AliDebug(3,Form("CookdEdx Q0 %f Q1 %f",out[0],out[1]));
db0e2c5f	504	break;
	505	case kLQ1D: // 1D LQ
	506	out[0]= 0.;
	507	for(Int_t islice = 0; islice < nSlice; islice++)
	508	if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor; // Protect against negative values for slices having no dE/dx information
239fe91c	509	out[0]*=1./Double_t(nSlice);
db0e2c5f	510	if(out[0] < 1e-6) return kFALSE;
f2762b1c	511	AliDebug(3,Form("CookdEdx dEdx %f",out[0]));
db0e2c5f	512	break;
db0e2c5f	513
e0de37e9	514	default:
db0e2c5f	515	return kFALSE;
e0de37e9	516	}
e0de37e9	517	return kTRUE;
ffb1ee30	518	}
ffb1ee30	519
ea235c90	520	//____________________________________________________________
bd58d4b9	521	Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level,Double_t centrality,ETRDPIDMethod PIDmethod) const {
db0e2c5f	522	//
	523	// Check whether particle is identified as electron assuming a certain electron efficiency level
	524	// Only electron and pion hypothesis is taken into account
	525	//
	526	// Inputs:
	527	// Number of tracklets
	528	// Likelihood values
	529	// Momentum
	530	// Electron efficiency level
	531	//
	532	// If the function fails when the params are not accessible, the function returns true
	533	//
	534	if(!fkPIDResponseObject){
31746593	535	AliDebug(3,"No PID Param object available");
ea235c90	536	return kTRUE;
	537	}
	538	Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
ce487a7f	539	Double_t params[4];
bd58d4b9	540	if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,centrality,PIDmethod)){
ea235c90	541	AliError("No Params found for the given configuration");
	542	return kTRUE;
	543	}
ce487a7f	544	Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
ea235c90	545	if(probEle > TMath::Max(TMath::Min(threshold, 0.99), 0.2)) return kTRUE; // truncate the threshold upperwards to 0.999 and lowerwards to 0.2 and exclude unphysical values
	546	return kFALSE;
	547	}
	548
f2762b1c	549	//____________________________________________________________
	550	Bool_t AliTRDPIDResponse::SetPIDResponseObject(const AliTRDPIDResponseObject * obj){
	551
	552	fkPIDResponseObject = obj;
	553	if((AliLog::GetDebugLevel("",IsA()->GetName()))>0)fkPIDResponseObject->Print("");
	554	return kTRUE;
	555	}