Small update of the calibration code for the case of runs without TRD (Raphaelle)
[u/mrichter/AliRoot.git] / STEER / AliTRDPIDResponse.cxx
CommitLineData
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>
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
51a0ce25 38#include "AliTRDPIDReference.h"
ffb1ee30 39#include "AliTRDPIDResponse.h"
40
41ClassImp(AliTRDPIDResponse)
42
ffb1ee30 43//____________________________________________________________
e0de37e9 44AliTRDPIDResponse::AliTRDPIDResponse():
45 TObject()
51a0ce25 46 ,fkPIDReference(NULL)
e0de37e9 47 ,fGainNormalisationFactor(1.)
48 ,fPIDmethod(kLQ1D)
ffb1ee30 49{
e0de37e9 50 //
51 // Default constructor
52 //
ffb1ee30 53}
54
55//____________________________________________________________
56AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
e0de37e9 57 TObject(ref)
51a0ce25 58 ,fkPIDReference(ref.fkPIDReference)
e0de37e9 59 ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
60 ,fPIDmethod(ref.fPIDmethod)
ffb1ee30 61{
e0de37e9 62 //
63 // Copy constructor
e0de37e9 64 //
ffb1ee30 65}
66
67//____________________________________________________________
68AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
e0de37e9 69 //
70 // Assignment operator
e0de37e9 71 //
72 if(this == &ref) return *this;
73
74 // Make copy
75 TObject::operator=(ref);
51a0ce25 76 fGainNormalisationFactor = ref.fGainNormalisationFactor;
77 fkPIDReference = ref.fkPIDReference;
e0de37e9 78 fPIDmethod = ref.fPIDmethod;
79
80 return *this;
ffb1ee30 81}
82
83//____________________________________________________________
84AliTRDPIDResponse::~AliTRDPIDResponse(){
e0de37e9 85 //
86 // Destructor
87 //
51a0ce25 88 if(IsOwner()) delete fkPIDReference;
ffb1ee30 89}
90
91//____________________________________________________________
51a0ce25 92Bool_t AliTRDPIDResponse::Load(const Char_t * filename, const Char_t *refName){
e0de37e9 93 //
94 // Load References into the toolkit
95 //
96 AliDebug(1, "Loading reference histos from root file");
2ad33025 97 TDirectory *owd = gDirectory;// store old working directory
e0de37e9 98
99 if(!filename)
100 filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
101 TFile *in = TFile::Open(filename);
102 if(!in){
103 AliError("Ref file not available.");
104 return kFALSE;
105 }
106
107 gROOT->cd();
51a0ce25 108 fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(in->Get(refName)->Clone());
109 in->Close(); delete in;
e0de37e9 110 owd->cd();
e0de37e9 111 SetBit(kIsOwner, kTRUE);
51a0ce25 112 AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDReference->GetNumberOfMomentumBins()));
e0de37e9 113 return kTRUE;
ffb1ee30 114}
115
116//____________________________________________________________
51a0ce25 117Bool_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES], Bool_t kNorm) const
ffb1ee30 118{
e0de37e9 119 //
120 // Calculate TRD likelihood values for the track based on dedx and
121 // momentum values. The likelihoods are calculated by query the
122 // reference data depending on the PID method selected
123 //
124 // Input parameter :
125 // n - number of dedx slices/chamber
126 // dedx - array of dedx slices organized layer wise
127 // p - array of momentum measurements organized layer wise
128 //
129 // Return parameters
130 // prob - probabilities allocated by TRD for particle specis
131 // kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
132 //
133 // Return value
134 // true if calculation success
135 //
ffb1ee30 136
51a0ce25 137 if(!fkPIDReference){
e0de37e9 138 AliWarning("Missing reference data. PID calculation not possible.");
139 return kFALSE;
140 }
ffb1ee30 141
e0de37e9 142 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
143 Double_t prLayer[AliPID::kSPECIES];
51a0ce25 144 Double_t dE[10], s(0.);
e0de37e9 145 for(Int_t il(kNlayer); il--;){
146 memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
51a0ce25 147 if(!CookdEdx(n, &dedx[il*n], &dE[0])) continue;
ffb1ee30 148
e0de37e9 149 s=0.;
150 for(Int_t is(AliPID::kSPECIES); is--;){
51a0ce25 151 if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], dE[0]);
e0de37e9 152 AliDebug(3, Form("Probability for Species %d in Layer %d: %f", is, il, prLayer[is]));
153 s+=prLayer[is];
154 }
155 if(s<1.e-30){
156 AliDebug(2, Form("Null all species prob layer[%d].", il));
157 continue;
158 }
159 for(Int_t is(AliPID::kSPECIES); is--;){
160 if(kNorm) prLayer[is] /= s;
161 prob[is] *= prLayer[is];
162 }
163 }
164 if(!kNorm) return kTRUE;
ffb1ee30 165
e0de37e9 166 s=0.;
167 for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
168 if(s<1.e-30){
169 AliDebug(2, "Null total prob.");
170 return kFALSE;
171 }
172 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
173 return kTRUE;
ffb1ee30 174}
175
ffb1ee30 176//____________________________________________________________
b439f460 177Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const {
e0de37e9 178 //
179 // Get the non-normalized probability for a certain particle species as coming
180 // from the reference histogram
181 // Interpolation between momentum bins
182 //
183 AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
51a0ce25 184 Float_t pLower, pUpper;
185 Double_t probLayer = 0.;
186 TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDReference->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper)),
187 *refLower = dynamic_cast<TH1 *>(fkPIDReference->GetLowerReference((AliPID::EParticleType)species, plocal, pLower));
e0de37e9 188 // Do Interpolation exept for underflow and overflow
51a0ce25 189 if(refLower && refUpper){
190 Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx));
191 Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx));
e0de37e9 192
51a0ce25 193 probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
194 } else if(refLower){
195 // underflow
196 probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx));
197 } else if(refUpper){
198 // overflow
199 probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx));
200 } else {
201 AliError("No references available");
e0de37e9 202 }
51a0ce25 203 AliDebug(1, Form("Probability %f", probLayer));
204 return probLayer;
ffb1ee30 205}
206
207//____________________________________________________________
208void AliTRDPIDResponse::SetOwner(){
e0de37e9 209 //
210 // Make Deep Copy of the Reference Histograms
211 //
51a0ce25 212 if(!fkPIDReference || IsOwner()) return;
213 const AliTRDPIDReference *tmp = fkPIDReference;
214 fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(tmp->Clone());
e0de37e9 215 SetBit(kIsOwner, kTRUE);
ffb1ee30 216}
217
218//____________________________________________________________
51a0ce25 219Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out) const
ffb1ee30 220{
51a0ce25 221 //
222 // Recalculate dE/dx
223 //
e0de37e9 224 switch(fPIDmethod){
225 case kNN: // NN
226 break;
227 case kLQ2D: // 2D LQ
228 break;
229 case kLQ1D: // 1D LQ
230 out[0]= 0.;
231 for(Int_t islice = 0; islice < nSlice; islice++) out[0] += in[islice] * fGainNormalisationFactor;
232 if(out[0] < 1e-6) return kFALSE;
233 break;
234 default:
235 return kFALSE;
236 }
237 return kTRUE;
ffb1ee30 238}
239