16f334a5bc9f1df3bcdd9e03110978cfdc54b098
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliTPCPIDResponse.cxx
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 //           Implementation of the TPC PID class
18 // Very naive one... Should be made better by the detector experts...
19 //      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
20 // With many additions and modifications suggested by
21 //      Alexander Kalweit, GSI, alexander.philipp.kalweit@cern.ch
22 //      Dariusz Miskowiec, GSI, D.Miskowiec@gsi.de
23 // ...and some modifications by
24 //      Mikolaj Krzewicki, GSI, mikolaj.krzewicki@cern.ch
25 // ...and some modifications plus eta correction functions by
26 //      Benjamin Hess, University of Tuebingen, bhess@cern.ch
27 //-----------------------------------------------------------------
28
29 #include <TGraph.h>
30 #include <TObjArray.h>
31 #include <TSpline.h>
32 #include <TBits.h>
33 #include <TMath.h>
34 #include <TH2D.h>
35
36 #include <AliLog.h>
37 #include "AliExternalTrackParam.h"
38 #include "AliVTrack.h"
39 #include "AliTPCPIDResponse.h"
40 #include "AliTPCdEdxInfo.h"
41 #include "TFile.h"
42 #include "TSpline.h"
43
44 ClassImp(AliTPCPIDResponse);
45
46
47 AliTPCPIDResponse *AliTPCPIDResponse::fgInstance =0;
48
49 const char* AliTPCPIDResponse::fgkGainScenarioName[fgkNumberOfGainScenarios+1]=
50 {
51   "", //default - no name
52   "1", //all high
53   "2", //OROC only
54   "unknown"
55 };
56
57 //_________________________________________________________________________
58 AliTPCPIDResponse::AliTPCPIDResponse():
59   TNamed(),
60   fMIP(50.),
61   fRes0(),
62   fResN2(),
63   fKp1(0.0283086),
64   fKp2(2.63394e+01),
65   fKp3(5.04114e-11),
66   fKp4(2.12543),
67   fKp5(4.88663),
68   fUseDatabase(kFALSE),
69   fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
70   fVoltageMap(72),
71   fLowGainIROCthreshold(-40),
72   fBadIROCthreshhold(-70),
73   fLowGainOROCthreshold(-40),
74   fBadOROCthreshhold(-40),
75   fMaxBadLengthFraction(0.5),
76   fMagField(0.),
77   fhEtaCorr(0x0),
78   fhEtaSigmaPar1(0x0),
79   fSigmaPar0(0.0),
80   fCurrentEventMultiplicity(0),
81   fCorrFuncMultiplicity(0x0),
82   fCorrFuncMultiplicityTanTheta(0x0),
83   fCorrFuncSigmaMultiplicity(0x0)
84 {
85   //
86   //  The default constructor
87   //
88   
89   AliLog::SetClassDebugLevel("AliTPCPIDResponse", AliLog::kInfo); 
90   
91   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=0.07;fResN2[i]=0.0;}
92   
93   fCorrFuncMultiplicity = new TF1("fCorrFuncMultiplicity", 
94                                   "[0] + [1]*TMath::Max([4], TMath::Min(x, [3])) + [2] * TMath::Power(TMath::Max([4], TMath::Min(x, [3])), 2)",
95                                   0., 0.2);
96   fCorrFuncMultiplicityTanTheta = new TF1("fCorrFuncMultiplicityTanTheta", "[0] * (x - [2]) + [1] * (x * x - [2] * [2])", -1.5, 1.5);
97   fCorrFuncSigmaMultiplicity = new TF1("fCorrFuncSigmaMultiplicity",
98                                        "TMath::Max(0, [0] + [1]*TMath::Min(x, [3]) + [2] * TMath::Power(TMath::Min(x, [3]), 2))", 0., 0.2);
99
100   
101   ResetMultiplicityCorrectionFunctions();
102   fgInstance=this;
103 }
104 /*TODO remove?
105 //_________________________________________________________________________
106 AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param):
107   TNamed(),
108   fMIP(param[0]),
109   fRes0(),
110   fResN2(),
111   fKp1(0.0283086),
112   fKp2(2.63394e+01),
113   fKp3(5.04114e-11),
114   fKp4(2.12543),
115   fKp5(4.88663),
116   fUseDatabase(kFALSE),
117   fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
118   fVoltageMap(72),
119   fLowGainIROCthreshold(-40),
120   fBadIROCthreshhold(-70),
121   fLowGainOROCthreshold(-40),
122   fBadOROCthreshhold(-40),
123   fMaxBadLengthFraction(0.5),
124   fMagField(0.),
125   fhEtaCorr(0x0),
126   fhEtaSigmaPar1(0x0),
127   fSigmaPar0(0.0)
128 {
129   //
130   //  The main constructor
131   //
132   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=param[1];fResN2[i]=param[2];}
133 }
134 */
135
136 //_________________________________________________________________________
137 AliTPCPIDResponse::~AliTPCPIDResponse()
138 {
139   //
140   // Destructor
141   //
142   
143   delete fhEtaCorr;
144   fhEtaCorr = 0x0;
145   
146   delete fhEtaSigmaPar1;
147   fhEtaSigmaPar1 = 0x0;
148   
149   delete fCorrFuncMultiplicity;
150   fCorrFuncMultiplicity = 0x0;
151   
152   delete fCorrFuncMultiplicityTanTheta;
153   fCorrFuncMultiplicityTanTheta = 0x0;
154   
155   delete fCorrFuncSigmaMultiplicity;
156   fCorrFuncSigmaMultiplicity = 0x0;
157   if (fgInstance==this) fgInstance=0;
158 }
159
160
161 //_________________________________________________________________________
162 AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that):
163   TNamed(that),
164   fMIP(that.fMIP),
165   fRes0(),
166   fResN2(),
167   fKp1(that.fKp1),
168   fKp2(that.fKp2),
169   fKp3(that.fKp3),
170   fKp4(that.fKp4),
171   fKp5(that.fKp5),
172   fUseDatabase(that.fUseDatabase),
173   fResponseFunctions(that.fResponseFunctions),
174   fVoltageMap(that.fVoltageMap),
175   fLowGainIROCthreshold(that.fLowGainIROCthreshold),
176   fBadIROCthreshhold(that.fBadIROCthreshhold),
177   fLowGainOROCthreshold(that.fLowGainOROCthreshold),
178   fBadOROCthreshhold(that.fBadOROCthreshhold),
179   fMaxBadLengthFraction(that.fMaxBadLengthFraction),
180   fMagField(that.fMagField),
181   fhEtaCorr(0x0),
182   fhEtaSigmaPar1(0x0),
183   fSigmaPar0(that.fSigmaPar0),
184   fCurrentEventMultiplicity(that.fCurrentEventMultiplicity),
185   fCorrFuncMultiplicity(0x0),
186   fCorrFuncMultiplicityTanTheta(0x0),
187   fCorrFuncSigmaMultiplicity(0x0)
188 {
189   //copy ctor
190   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
191  
192   // Copy eta maps
193   if (that.fhEtaCorr) {
194     fhEtaCorr = new TH2D(*(that.fhEtaCorr));
195     fhEtaCorr->SetDirectory(0);
196   }
197   
198   if (that.fhEtaSigmaPar1) {
199     fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1));
200     fhEtaSigmaPar1->SetDirectory(0);
201   }
202   
203   // Copy multiplicity correction functions
204   if (that.fCorrFuncMultiplicity) {
205     fCorrFuncMultiplicity = new TF1(*(that.fCorrFuncMultiplicity));
206   }
207   
208   if (that.fCorrFuncMultiplicityTanTheta) {
209     fCorrFuncMultiplicityTanTheta = new TF1(*(that.fCorrFuncMultiplicityTanTheta));
210   }
211   
212   if (that.fCorrFuncSigmaMultiplicity) {
213     fCorrFuncSigmaMultiplicity = new TF1(*(that.fCorrFuncSigmaMultiplicity));
214   }
215 }
216
217 //_________________________________________________________________________
218 AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that)
219 {
220   //assignment
221   if (&that==this) return *this;
222   TNamed::operator=(that);
223   fMIP=that.fMIP;
224   fKp1=that.fKp1;
225   fKp2=that.fKp2;
226   fKp3=that.fKp3;
227   fKp4=that.fKp4;
228   fKp5=that.fKp5;
229   fUseDatabase=that.fUseDatabase;
230   fResponseFunctions=that.fResponseFunctions;
231   fVoltageMap=that.fVoltageMap;
232   fLowGainIROCthreshold=that.fLowGainIROCthreshold;
233   fBadIROCthreshhold=that.fBadIROCthreshhold;
234   fLowGainOROCthreshold=that.fLowGainOROCthreshold;
235   fBadOROCthreshhold=that.fBadOROCthreshhold;
236   fMaxBadLengthFraction=that.fMaxBadLengthFraction;
237   fMagField=that.fMagField;
238   fCurrentEventMultiplicity=that.fCurrentEventMultiplicity;
239   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
240
241   delete fhEtaCorr;
242   fhEtaCorr=0x0;
243   if (that.fhEtaCorr) {
244     fhEtaCorr = new TH2D(*(that.fhEtaCorr));
245     fhEtaCorr->SetDirectory(0);
246   }
247   
248   delete fhEtaSigmaPar1;
249   fhEtaSigmaPar1=0x0;
250   if (that.fhEtaSigmaPar1) {
251     fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1));
252     fhEtaSigmaPar1->SetDirectory(0);
253   }
254   
255   fSigmaPar0 = that.fSigmaPar0;
256   
257   delete fCorrFuncMultiplicity;
258   fCorrFuncMultiplicity = 0x0;
259   if (that.fCorrFuncMultiplicity) {
260     fCorrFuncMultiplicity = new TF1(*(that.fCorrFuncMultiplicity));
261   }
262   
263   delete fCorrFuncMultiplicityTanTheta;
264   fCorrFuncMultiplicityTanTheta = 0x0;
265   if (that.fCorrFuncMultiplicityTanTheta) {
266     fCorrFuncMultiplicityTanTheta = new TF1(*(that.fCorrFuncMultiplicityTanTheta));
267   }
268   
269   delete fCorrFuncSigmaMultiplicity;
270   fCorrFuncSigmaMultiplicity = 0x0;
271   if (that.fCorrFuncSigmaMultiplicity) {
272     fCorrFuncSigmaMultiplicity = new TF1(*(that.fCorrFuncSigmaMultiplicity));
273   }
274
275   return *this;
276 }
277
278 //_________________________________________________________________________
279 Double_t AliTPCPIDResponse::Bethe(Double_t betaGamma) const {
280   //
281   // This is the Bethe-Bloch function normalised to 1 at the minimum
282   // WARNING
283   // Simulated and reconstructed Bethe-Bloch differs
284   //           Simulated  curve is the dNprim/dx
285   //           Reconstructed is proportianal dNtot/dx
286   // Temporary fix for production -  Simple linear correction function
287   // Future    2 Bethe Bloch formulas needed
288   //           1. for simulation
289   //           2. for reconstructed PID
290   //
291   
292 //   const Float_t kmeanCorrection =0.1;
293   Double_t bb=
294     AliExternalTrackParam::BetheBlochAleph(betaGamma,fKp1,fKp2,fKp3,fKp4,fKp5);
295   return bb*fMIP;
296 }
297
298 //_________________________________________________________________________
299 void AliTPCPIDResponse::SetBetheBlochParameters(Double_t kp1,
300                              Double_t kp2,
301                              Double_t kp3,
302                              Double_t kp4,
303                              Double_t kp5) {
304   //
305   // Set the parameters of the ALEPH Bethe-Bloch formula
306   //
307   fKp1=kp1;
308   fKp2=kp2;
309   fKp3=kp3;
310   fKp4=kp4;
311   fKp5=kp5;
312 }
313
314 //_________________________________________________________________________
315 void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2) {
316   //
317   // Set the relative resolution  sigma_rel = res0 * sqrt(1+resN2/npoint)
318   //
319   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=res0;fResN2[i]=resN2;}
320 }
321
322 //_________________________________________________________________________
323 Double_t AliTPCPIDResponse::GetExpectedSignal(Float_t mom,
324                                               AliPID::EParticleType n) const {
325   //
326   // Deprecated function (for backward compatibility). Please use 
327   // GetExpectedSignal(const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource,
328   //                   Bool_t correctEta, Bool_t correctMultiplicity);
329   // instead!
330   //
331   //
332   // Calculates the expected PID signal as the function of 
333   // the information stored in the track, for the specified particle type 
334   //  
335   // At the moment, these signals are just the results of calling the 
336   // Bethe-Bloch formula. 
337   // This can be improved. By taking into account the number of
338   // assigned clusters and/or the track dip angle, for example.  
339   //
340   
341   //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
342   // !!! Splines for light nuclei need to be normalised to this factor !!!
343   const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(n),2.3);
344   
345   Double_t mass=AliPID::ParticleMassZ(n);
346   if (!fUseDatabase) return Bethe(mom/mass) * chargeFactor;
347   //
348   const TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n);
349
350   if (!responseFunction) return Bethe(mom/mass) * chargeFactor;
351   
352   return fMIP*responseFunction->Eval(mom/mass)*chargeFactor;
353
354 }
355
356 //_________________________________________________________________________
357 Double_t AliTPCPIDResponse::GetExpectedSigma(Float_t mom, 
358                                              Int_t nPoints,
359                                              AliPID::EParticleType n) const {
360   //
361   // Deprecated function (for backward compatibility). Please use 
362   // GetExpectedSigma(onst AliVTrack* track, AliPID::EParticleType species, 
363   // ETPCdEdxSource dedxSource, Bool_t correctEta) instead!
364   //
365   //
366   // Calculates the expected sigma of the PID signal as the function of 
367   // the information stored in the track, for the specified particle type 
368   //  
369   
370   if (nPoints != 0) 
371     return GetExpectedSignal(mom,n)*fRes0[0]*sqrt(1. + fResN2[0]/nPoints);
372   else
373     return GetExpectedSignal(mom,n)*fRes0[0];
374 }
375
376 ////////////////////////////////////////////////////NEW//////////////////////////////
377
378 //_________________________________________________________________________
379 void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario) {
380   //
381   // Set the relative resolution  sigma_rel = res0 * sqrt(1+resN2/npoint)
382   //
383   if ((Int_t)gainScenario>=(Int_t)fgkNumberOfGainScenarios) return; //TODO: better handling!
384   fRes0[gainScenario]=res0;
385   fResN2[gainScenario]=resN2;
386 }
387
388
389 //_________________________________________________________________________
390 Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
391                                               AliPID::EParticleType species,
392                                               Double_t /*dEdx*/,
393                                               const TSpline3* responseFunction,
394                                               Bool_t correctEta,
395                                               Bool_t correctMultiplicity) const 
396 {
397   // Calculates the expected PID signal as the function of 
398   // the information stored in the track and the given parameters,
399   // for the specified particle type 
400   //  
401   // At the moment, these signals are just the results of calling the 
402   // Bethe-Bloch formula plus, if desired, taking into account the eta dependence
403   // and the multiplicity dependence (for PbPb). 
404   // This can be improved. By taking into account the number of
405   // assigned clusters and/or the track dip angle, for example.  
406   //
407   
408   Double_t mom=track->GetTPCmomentum();
409   Double_t mass=AliPID::ParticleMassZ(species);
410   
411   //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
412   // !!! Splines for light nuclei need to be normalised to this factor !!!
413   const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3);
414   
415   if (!responseFunction)
416     return Bethe(mom/mass) * chargeFactor;
417   
418   Double_t dEdxSplines = fMIP*responseFunction->Eval(mom/mass) * chargeFactor;
419     
420   if (!correctEta && !correctMultiplicity)
421     return dEdxSplines;
422   
423   Double_t corrFactorEta = 1.0;
424   Double_t corrFactorMultiplicity = 1.0;
425   
426   if (correctEta) {
427     corrFactorEta = GetEtaCorrectionFast(track, dEdxSplines);
428     //TODO Alternatively take current track dEdx
429     //corrFactorEta = GetEtaCorrectionFast(track, dEdx);
430   }
431   
432   if (correctMultiplicity)
433     corrFactorMultiplicity = GetMultiplicityCorrectionFast(track, dEdxSplines * corrFactorEta, fCurrentEventMultiplicity);
434
435   return dEdxSplines * corrFactorEta * corrFactorMultiplicity;
436 }
437
438
439 //_________________________________________________________________________
440 Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
441                                               AliPID::EParticleType species,
442                                               ETPCdEdxSource dedxSource,
443                                               Bool_t correctEta,
444                                               Bool_t correctMultiplicity) const
445 {
446   // Calculates the expected PID signal as the function of 
447   // the information stored in the track, for the specified particle type 
448   //  
449   // At the moment, these signals are just the results of calling the 
450   // Bethe-Bloch formula plus, if desired, taking into account the eta dependence
451   // and the multiplicity dependence (for PbPb). 
452   // This can be improved. By taking into account the number of
453   // assigned clusters and/or the track dip angle, for example.  
454   //
455   
456   if (!fUseDatabase) {
457     //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
458     // !!! Splines for light nuclei need to be normalised to this factor !!!
459     const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3);
460   
461     return Bethe(track->GetTPCmomentum() / AliPID::ParticleMassZ(species)) * chargeFactor;
462   }
463   
464   Double_t dEdx = -1;
465   Int_t nPoints = -1;
466   ETPCgainScenario gainScenario = kGainScenarioInvalid;
467   TSpline3* responseFunction = 0x0;
468     
469   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) {
470     // Something is wrong with the track -> Return obviously invalid value
471     return -999;
472   }
473   
474   // Charge factor already taken into account inside the following function call
475   return GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity);
476 }
477   
478 //_________________________________________________________________________
479 TSpline3* AliTPCPIDResponse::GetResponseFunction( AliPID::EParticleType type,
480                                                   AliTPCPIDResponse::ETPCgainScenario gainScenario ) const
481 {
482   //get response function
483   return dynamic_cast<TSpline3*>(fResponseFunctions.At(ResponseFunctionIndex(type,gainScenario)));
484 }
485
486 //_________________________________________________________________________
487 TSpline3* AliTPCPIDResponse::GetResponseFunction( const AliVTrack* track,
488                                AliPID::EParticleType species,
489                                ETPCdEdxSource dedxSource) const 
490 {
491   //the splines are stored in an array, different scenarios
492
493   Double_t dEdx = -1;
494   Int_t nPoints = -1;
495   ETPCgainScenario gainScenario = kGainScenarioInvalid;
496   TSpline3* responseFunction = 0x0;
497   
498   if (ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
499     return responseFunction;
500   
501   return NULL;
502 }
503
504 //_________________________________________________________________________
505 void AliTPCPIDResponse::ResetSplines()
506 {
507   //reset the array with splines
508   for (Int_t i=0;i<fResponseFunctions.GetEntriesFast();i++)
509   {
510     fResponseFunctions.AddAt(NULL,i);
511   }
512 }
513 //_________________________________________________________________________
514 Int_t AliTPCPIDResponse::ResponseFunctionIndex( AliPID::EParticleType species,
515                                                 ETPCgainScenario gainScenario ) const
516 {
517   //get the index in fResponseFunctions given type and scenario
518   return Int_t(species)+Int_t(gainScenario)*fgkNumberOfParticleSpecies;
519 }
520
521 //_________________________________________________________________________
522 void AliTPCPIDResponse::SetResponseFunction( TObject* o,
523                                              AliPID::EParticleType species,
524                                              ETPCgainScenario gainScenario )
525 {
526   fResponseFunctions.AddAtAndExpand(o,ResponseFunctionIndex(species,gainScenario));
527 }
528
529
530 //_________________________________________________________________________
531 Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, 
532                                              AliPID::EParticleType species,
533                                              ETPCgainScenario gainScenario,
534                                              Double_t dEdx,
535                                              Int_t nPoints,
536                                              const TSpline3* responseFunction,
537                                              Bool_t correctEta,
538                                              Bool_t correctMultiplicity) const 
539 {
540   // Calculates the expected sigma of the PID signal as the function of 
541   // the information stored in the track and the given parameters,
542   // for the specified particle type 
543   //
544   
545   if (!responseFunction)
546     return 999;
547     
548   //TODO Check whether it makes sense to set correctMultiplicity to kTRUE while correctEta might be kFALSE
549   
550   // If eta correction (=> new sigma parametrisation) is requested, but no sigma map is available, print error message
551   if (correctEta && !fhEtaSigmaPar1) {
552     AliError("New sigma parametrisation requested, but sigma map not initialised (usually via AliPIDResponse). Old sigma parametrisation will be used!");
553   }
554   
555   // If no sigma map is available or if no eta correction is requested (sigma maps only for corrected eta!), use the old parametrisation
556   if (!fhEtaSigmaPar1 || !correctEta) {  
557     if (nPoints != 0) 
558       return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity) *
559                fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints);
560     else
561       return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity)*fRes0[gainScenario];
562   }
563     
564   if (nPoints > 0) {
565     // Use eta correction (+ eta-dependent sigma)
566     Double_t sigmaPar1 = GetSigmaPar1Fast(track, species, dEdx, responseFunction);
567     
568     if (correctMultiplicity) {
569       // In addition, take into account multiplicity dependence of mean and sigma of dEdx
570       Double_t dEdxExpectedEtaCorrected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE);
571       
572       // GetMultiplicityCorrection and GetMultiplicitySigmaCorrection both need the eta corrected dEdxExpected
573       Double_t multiplicityCorrFactor = GetMultiplicityCorrectionFast(track, dEdxExpectedEtaCorrected, fCurrentEventMultiplicity);
574       Double_t multiplicitySigmaCorrFactor = GetMultiplicitySigmaCorrectionFast(dEdxExpectedEtaCorrected, fCurrentEventMultiplicity);
575       
576       // multiplicityCorrFactor to correct dEdxExpected for multiplicity. In addition: Correction factor for sigma
577       return (dEdxExpectedEtaCorrected * multiplicityCorrFactor) 
578               * (sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints) * multiplicitySigmaCorrFactor);
579     }
580     else {
581       return GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE)*
582              sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints);
583     }
584   }
585   else { 
586     // One should never have/take tracks with 0 dEdx clusters!
587     return 999;
588   }
589 }
590
591
592 //_________________________________________________________________________
593 Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, 
594                                              AliPID::EParticleType species,
595                                              ETPCdEdxSource dedxSource,
596                                              Bool_t correctEta,
597                                              Bool_t correctMultiplicity) const 
598 {
599   // Calculates the expected sigma of the PID signal as the function of 
600   // the information stored in the track, for the specified particle type 
601   // and dedx scenario
602   //
603   
604   Double_t dEdx = -1;
605   Int_t nPoints = -1;
606   ETPCgainScenario gainScenario = kGainScenarioInvalid;
607   TSpline3* responseFunction = 0x0;
608   
609   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
610     return 999; //TODO: better handling!
611   
612   return GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity);
613 }
614
615
616 //_________________________________________________________________________
617 Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track, 
618                              AliPID::EParticleType species,
619                              ETPCdEdxSource dedxSource,
620                              Bool_t correctEta,
621                              Bool_t correctMultiplicity) const
622 {
623   //Calculates the number of sigmas of the PID signal from the expected value
624   //for a given particle species in the presence of multiple gain scenarios
625   //inside the TPC
626   
627   Double_t dEdx = -1;
628   Int_t nPoints = -1;
629   ETPCgainScenario gainScenario = kGainScenarioInvalid;
630   TSpline3* responseFunction = 0x0;
631   
632   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
633     return -999; //TODO: Better handling!
634     
635   Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity);
636   Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity);
637   // 999 will be returned by GetExpectedSigma e.g. in case of 0 dEdx clusters
638   if (sigma >= 998) 
639     return -999;
640   else
641     return (dEdx-bethe)/sigma;
642 }
643
644 //_________________________________________________________________________
645 Float_t AliTPCPIDResponse::GetSignalDelta(const AliVTrack* track,
646                                           AliPID::EParticleType species,
647                                           ETPCdEdxSource dedxSource,
648                                           Bool_t correctEta,
649                                           Bool_t correctMultiplicity,
650                                           Bool_t ratio/*=kFALSE*/)const
651 {
652   //Calculates the number of sigmas of the PID signal from the expected value
653   //for a given particle species in the presence of multiple gain scenarios
654   //inside the TPC
655
656   Double_t dEdx = -1;
657   Int_t nPoints = -1;
658   ETPCgainScenario gainScenario = kGainScenarioInvalid;
659   TSpline3* responseFunction = 0x0;
660
661   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
662     return -9999.; //TODO: Better handling!
663
664   const Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity);
665
666   Double_t delta=-9999.;
667   if (!ratio) delta=dEdx-bethe;
668   else if (bethe>1.e-20) delta=dEdx/bethe;
669   
670   return delta;
671 }
672
673 //_________________________________________________________________________
674 Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track, 
675                                                  AliPID::EParticleType species,
676                                                  ETPCdEdxSource dedxSource,
677                                                  Double_t& dEdx,
678                                                  Int_t& nPoints,
679                                                  ETPCgainScenario& gainScenario,
680                                                  TSpline3** responseFunction) const 
681 {
682   // Calculates the right parameters for PID
683   //   dEdx parametrization for the proper gain scenario, dEdx 
684   //   and NPoints used for dEdx
685   // based on the track geometry (which chambers it crosses) for the specified particle type 
686   // and preferred source of dedx.
687   // returns true on success
688   
689   
690   if (dedxSource == kdEdxDefault) {
691     // Fast handling for default case. In addition: Keep it simple (don't call additional functions) to
692     // avoid possible bugs
693     
694     // GetTPCsignalTunedOnData will be non-positive, if it has not been set (i.e. in case of MC NOT tuned to data).
695     // If this is the case, just take the normal signal
696     dEdx = track->GetTPCsignalTunedOnData();
697     if (dEdx <= 0) {
698       dEdx = track->GetTPCsignal();
699     }
700     
701     nPoints = track->GetTPCsignalN();
702     gainScenario = kDefault;
703     
704     TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario));
705     *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast?
706   
707     return kTRUE;
708   }
709   
710   //TODO Proper handle of tuneMConData for other dEdx sources
711   
712   Double32_t signal[4]; //0: IROC, 1: OROC medium, 2:OROC long, 3: OROC all (def. truncation used)
713   Char_t ncl[3];        //same
714   Char_t nrows[3];      //same
715   const AliTPCdEdxInfo* dEdxInfo = track->GetTPCdEdxInfo();
716   
717   if (!dEdxInfo && dedxSource!=kdEdxDefault)  //in one case its ok if we dont have the info
718   {
719     AliError("AliTPCdEdxInfo not available");
720     return kFALSE;
721   }
722
723   if (dEdxInfo) dEdxInfo->GetTPCSignalRegionInfo(signal,ncl,nrows);
724
725   //check if we cross a bad OROC in which case we reject
726   EChamberStatus trackOROCStatus = TrackStatus(track,2);
727   if (trackOROCStatus==kChamberOff || trackOROCStatus==kChamberLowGain)
728   {
729     return kFALSE;
730   }
731
732   switch (dedxSource)
733   {
734     case kdEdxOROC:
735       {
736         if (trackOROCStatus==kChamberInvalid) return kFALSE; //never reached OROC
737         dEdx = signal[3];
738         nPoints = ncl[2]+ncl[1];
739         gainScenario = kOROChigh;
740         break;
741       }
742     case kdEdxHybrid:
743       {
744         //if we cross a bad IROC we use OROC dedx, if we dont we use combined
745         EChamberStatus status = TrackStatus(track,1);
746         if (status!=kChamberHighGain)
747         {
748           dEdx = signal[3];
749           nPoints = ncl[2]+ncl[1];
750           gainScenario = kOROChigh;
751         }
752         else
753         {
754           dEdx = track->GetTPCsignal();
755           nPoints = track->GetTPCsignalN();
756           gainScenario = kALLhigh;
757         }
758         break;
759       }
760     default:
761       {
762          dEdx = 0.;
763          nPoints = 0;
764          gainScenario = kGainScenarioInvalid;
765          return kFALSE;
766       }
767   }
768   TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario));
769   *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast?
770   
771   return kTRUE;
772 }
773
774
775 //_________________________________________________________________________
776 Double_t AliTPCPIDResponse::GetEtaCorrectionFast(const AliVTrack *track, Double_t dEdxSplines) const
777 {
778   // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly.
779   //
780   // Get eta correction for the given parameters.
781   //
782   
783   if (!fhEtaCorr) {
784     // Calling this function means to request eta correction in some way. Print error message, if no map is available!
785     AliError("Eta correction requested, but map not initialised (usually via AliPIDResponse). Returning eta correction factor 1!");
786     return 1.;
787   }
788   
789   Double_t tpcSignal = dEdxSplines;
790   
791   if (tpcSignal < 1.)
792     return 1.;
793   
794   Double_t tanTheta = GetTrackTanTheta(track); 
795   Int_t binX = fhEtaCorr->GetXaxis()->FindFixBin(tanTheta);
796   Int_t binY = fhEtaCorr->GetYaxis()->FindFixBin(1. / tpcSignal);
797   
798   if (binX == 0) 
799     binX = 1;
800   if (binX > fhEtaCorr->GetXaxis()->GetNbins())
801     binX = fhEtaCorr->GetXaxis()->GetNbins();
802   
803   if (binY == 0)
804     binY = 1;
805   if (binY > fhEtaCorr->GetYaxis()->GetNbins())
806     binY = fhEtaCorr->GetYaxis()->GetNbins();
807   
808   return fhEtaCorr->GetBinContent(binX, binY);
809 }
810
811
812 //_________________________________________________________________________
813 Double_t AliTPCPIDResponse::GetEtaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
814 {
815   //
816   // Get eta correction for the given track.
817   //
818   
819   if (!fhEtaCorr)
820     return 1.;
821   
822   Double_t dEdx = -1;
823   Int_t nPoints = -1;
824   ETPCgainScenario gainScenario = kGainScenarioInvalid;
825   TSpline3* responseFunction = 0x0;
826   
827   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
828     return 1.; 
829   
830   // For the eta correction, do NOT take the multiplicity corrected value of dEdx
831   Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE);
832   
833   //TODO Alternatively take current track dEdx
834   //return GetEtaCorrectionFast(track, dEdx);
835   
836   return GetEtaCorrectionFast(track, dEdxSplines);
837 }
838
839
840 //_________________________________________________________________________
841 Double_t AliTPCPIDResponse::GetEtaCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
842 {
843   //
844   // Get eta corrected dEdx for the given track. For the correction, the expected dEdx of
845   // the specified species will be used. If the species is set to AliPID::kUnknown, the
846   // dEdx of the track is used instead.
847   // WARNING: In the latter case, the eta correction might not be as good as if the
848   // expected dEdx is used, which is the way the correction factor is designed
849   // for.
850   // In any case, one has to decide carefully to which expected signal one wants to
851   // compare the corrected value - to the corrected or uncorrected.
852   // Anyhow, a safer way of looking e.g. at the n-sigma is to call
853   // the corresponding function GetNumberOfSigmas!
854   //
855   
856   Double_t dEdx = -1;
857   Int_t nPoints = -1;
858   ETPCgainScenario gainScenario = kGainScenarioInvalid;
859   TSpline3* responseFunction = 0x0;
860   
861   // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case
862   // it is not used anyway, so this causes no trouble.
863   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
864     return -1.;
865   
866   Double_t etaCorr = 0.;
867   
868   if (species < AliPID::kUnknown) {
869     // For the eta correction, do NOT take the multiplicity corrected value of dEdx
870     Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE);
871     etaCorr = GetEtaCorrectionFast(track, dEdxSplines);
872   }
873   else {
874     etaCorr = GetEtaCorrectionFast(track, dEdx);
875   }
876     
877   if (etaCorr <= 0)
878     return -1.;
879   
880   return dEdx / etaCorr; 
881 }
882
883
884 //_________________________________________________________________________
885 Double_t AliTPCPIDResponse::GetSigmaPar1Fast(const AliVTrack *track, AliPID::EParticleType species, Double_t dEdx,
886                                              const TSpline3* responseFunction) const
887 {
888   // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly.
889   //
890   // Get parameter 1 of sigma parametrisation of TPC dEdx from the histogram for the given track.
891   //
892   
893   if (!fhEtaSigmaPar1) {
894     // Calling this function means to request new sigma parametrisation in some way. Print error message, if no map is available!
895     AliError("New sigma parametrisation requested, but sigma map not initialised (usually via AliPIDResponse). Returning error value for sigma parameter1 = 999!");
896     return 999;
897   }
898   
899   // The sigma maps are created with data sets that are already eta corrected and for which the 
900   // splines have been re-created. Therefore, the value for the lookup needs to be the value of
901   // the splines without any additional eta correction.
902   // NOTE: This is due to the method the maps are created. The track dEdx (not the expected one!)
903   // is corrected to uniquely related a momemtum bin with an expected dEdx, where the expected dEdx
904   // equals the track dEdx for all eta (thanks to the correction and the re-fit of the splines).
905   // Consequently, looking up the uncorrected expected dEdx at a given tanTheta yields the correct
906   // sigma parameter!
907   // Also: It has to be the spline dEdx, since one wants to get the sigma for the assumption(!)
908   // of such a particle, which by assumption then has this dEdx value
909   
910   // For the eta correction, do NOT take the multiplicity corrected value of dEdx
911   Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE);
912   
913   if (dEdxExpected < 1.)
914     return 999;
915   
916   Double_t tanTheta = GetTrackTanTheta(track);
917   Int_t binX = fhEtaSigmaPar1->GetXaxis()->FindFixBin(tanTheta);
918   Int_t binY = fhEtaSigmaPar1->GetYaxis()->FindFixBin(1. / dEdxExpected);
919     
920   if (binX == 0) 
921     binX = 1;
922   if (binX > fhEtaSigmaPar1->GetXaxis()->GetNbins())
923     binX = fhEtaSigmaPar1->GetXaxis()->GetNbins();
924     
925   if (binY == 0)
926     binY = 1;
927   if (binY > fhEtaSigmaPar1->GetYaxis()->GetNbins())
928     binY = fhEtaSigmaPar1->GetYaxis()->GetNbins();
929     
930   return fhEtaSigmaPar1->GetBinContent(binX, binY);
931 }
932
933
934 //_________________________________________________________________________
935 Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
936 {
937   //
938   // Get parameter 1 of sigma parametrisation of TPC dEdx from the histogram for the given track.
939   //
940   
941   if (!fhEtaSigmaPar1)
942     return 999;
943   
944   Double_t dEdx = -1;
945   Int_t nPoints = -1;
946   ETPCgainScenario gainScenario = kGainScenarioInvalid;
947   TSpline3* responseFunction = 0x0;
948   
949   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
950     return 999; 
951   
952   return GetSigmaPar1Fast(track, species, dEdx, responseFunction);
953 }
954
955
956 //_________________________________________________________________________
957 Bool_t AliTPCPIDResponse::SetEtaCorrMap(TH2D* hMap)
958 {
959   //
960   // Load map for TPC eta correction (a copy is stored and will be deleted automatically).
961   // If hMap is 0x0,the eta correction will be disabled and kFALSE is returned. 
962   // If the map can be set, kTRUE is returned.
963   //
964   
965   delete fhEtaCorr;
966   
967   if (!hMap) {
968     fhEtaCorr = 0x0;
969     
970     return kFALSE;
971   }
972   
973   fhEtaCorr = (TH2D*)(hMap->Clone());
974   fhEtaCorr->SetDirectory(0);
975       
976   return kTRUE;
977 }
978
979
980 //_________________________________________________________________________
981 Bool_t AliTPCPIDResponse::SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0)
982 {
983   //
984   // Load map for TPC sigma map (a copy is stored and will be deleted automatically):
985   // Parameter 1 is stored as a 2D map (1/dEdx vs. tanTheta_local) and parameter 0 is
986   // a constant. If hSigmaPar1Map is 0x0, the old sigma parametrisation will be used
987   // (and sigmaPar0 is ignored!) and kFALSE is returned. 
988   // If the map can be set, sigmaPar0 is also set and kTRUE will be returned.
989   //
990   
991   delete fhEtaSigmaPar1;
992   
993   if (!hSigmaPar1Map) {
994     fhEtaSigmaPar1 = 0x0;
995     fSigmaPar0 = 0.0;
996     
997     return kFALSE;
998   }
999   
1000   fhEtaSigmaPar1 = (TH2D*)(hSigmaPar1Map->Clone());
1001   fhEtaSigmaPar1->SetDirectory(0);
1002   fSigmaPar0 = sigmaPar0;
1003   
1004   return kTRUE;
1005 }
1006
1007
1008 //_________________________________________________________________________
1009 Double_t AliTPCPIDResponse::GetTrackTanTheta(const AliVTrack *track) const
1010 {
1011   // Extract the tanTheta from the information available in the AliVTrack
1012   
1013   // For ESD tracks, the local tanTheta could be used (esdTrack->GetInnerParam()->GetTgl()).
1014   // However, this value is not available for AODs and, thus, not for AliVTrack.
1015   // Fortunately, the following formula allows to approximate the local tanTheta with the 
1016   // global theta angle -> This is for by far most of the tracks the same, but gives at
1017   // maybe the percent level differences within +- 0.2 in tanTheta -> Which is still ok.
1018   
1019   /*
1020   const AliExternalTrackParam* innerParam = track->GetInnerParam();
1021   Double_t tanTheta = 0;
1022   if (innerParam) 
1023     tanTheta = innerParam->GetTgl();
1024   else
1025     tanTheta = TMath::Tan(-track->Theta() + TMath::Pi() / 2.0);
1026   
1027   // Constant in formula for B in kGauss (factor 0.1 to convert B from Tesla to kGauss),
1028   // pT in GeV/c (factor c*1e-9 to arrive at GeV/c) and curvature in 1/cm (factor 0.01 to get from m to cm)
1029   const Double_t constant = TMath::C()* 1e-9 * 0.1 * 0.01; 
1030   const Double_t curvature = fMagField * constant / track->Pt(); // in 1./cm
1031   
1032   Double_t averageddzdr = 0.;
1033   Int_t nParts = 0;
1034
1035   for (Double_t r = 85; r < 245; r++) {
1036     Double_t sinPhiLocal = TMath::Abs(r*curvature*0.5);
1037     
1038     // Cut on |sin(phi)| as during reco
1039     if (TMath::Abs(sinPhiLocal) <= 0.95) {
1040       const Double_t phiLocal = TMath::ASin(sinPhiLocal);
1041       const Double_t tanPhiLocal = TMath::Tan(phiLocal);
1042       
1043       averageddzdr += tanTheta * TMath::Sqrt(1. + tanPhiLocal * tanPhiLocal); 
1044       nParts++;
1045     }
1046   }
1047   
1048   if (nParts > 0)
1049     averageddzdr /= nParts; 
1050   else {
1051     AliError("Problems during determination of dz/dr. Returning pure tanTheta as best estimate!");
1052     return tanTheta;
1053   }
1054   
1055   //printf("pT: %f\nFactor/magField(kGs)/curvature^-1: %f / %f /%f\ntanThetaGlobalFromTheta/tanTheta/Averageddzdr: %f / %f / %f\n\n",
1056   //          track->Pt(), constant, fMagField, 1./curvature, TMath::Tan(-track->Theta() + TMath::Pi() / 2.0), tanTheta, averageddzdr);
1057   
1058   return averageddzdr;
1059   */
1060   
1061   
1062   // Alternatively (in average, the effect is found to be negligable!):
1063   // Take local tanTheta from TPC inner wall, if available (currently only for ESDs available)
1064   //const AliExternalTrackParam* innerParam = track->GetInnerParam();
1065   //if (innerParam) {
1066   //  return innerParam->GetTgl();
1067   //}
1068   
1069   return TMath::Tan(-track->Theta() + TMath::Pi() / 2.0);
1070 }
1071
1072
1073 //_________________________________________________________________________
1074 Double_t AliTPCPIDResponse::GetMultiplicityCorrectionFast(const AliVTrack *track, Double_t dEdxExpected, Int_t multiplicity) const
1075 {
1076   // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly.
1077   //
1078   // Calculate the multiplicity correction factor for this track for the given multiplicity.
1079   // The parameter dEdxExpected should take into account the eta correction already!
1080   
1081   // Multiplicity depends on pure dEdx. Therefore, correction factor depends indirectly on eta
1082   // => Use eta corrected value for dEdexExpected.
1083   
1084   if (dEdxExpected <= 0 || multiplicity <= 0)
1085     return 1.0;
1086   
1087   const Double_t dEdxExpectedInv = 1. / dEdxExpected;
1088   Double_t relSlope = fCorrFuncMultiplicity->Eval(dEdxExpectedInv);
1089   
1090   const Double_t tanTheta = GetTrackTanTheta(track);
1091   relSlope += fCorrFuncMultiplicityTanTheta->Eval(tanTheta);
1092
1093   return (1. + relSlope * multiplicity);
1094 }
1095
1096
1097 //_________________________________________________________________________
1098 Double_t AliTPCPIDResponse::GetMultiplicityCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
1099 {
1100   //
1101   // Get multiplicity correction for the given track (w.r.t. the mulitplicity of the current event)
1102   //
1103   
1104   //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse
1105   
1106   // No correction in case of multiplicity <= 0
1107   if (fCurrentEventMultiplicity <= 0)
1108     return 1.;
1109   
1110   Double_t dEdx = -1;
1111   Int_t nPoints = -1;
1112   ETPCgainScenario gainScenario = kGainScenarioInvalid;
1113   TSpline3* responseFunction = 0x0;
1114   
1115   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
1116     return 1.; 
1117   
1118   //TODO Does it make sense to use the multiplicity correction WITHOUT eta correction?! Are the fit parameters still valid?
1119   // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course)
1120   Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE);
1121   
1122   return GetMultiplicityCorrectionFast(track, dEdxExpected, fCurrentEventMultiplicity);
1123 }
1124
1125
1126 //_________________________________________________________________________
1127 Double_t AliTPCPIDResponse::GetMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
1128 {
1129   //
1130   // Get multiplicity corrected dEdx for the given track. For the correction, the expected dEdx of
1131   // the specified species will be used. If the species is set to AliPID::kUnknown, the
1132   // dEdx of the track is used instead.
1133   // WARNING: In the latter case, the correction might not be as good as if the
1134   // expected dEdx is used, which is the way the correction factor is designed
1135   // for.
1136   // In any case, one has to decide carefully to which expected signal one wants to
1137   // compare the corrected value - to the corrected or uncorrected.
1138   // Anyhow, a safer way of looking e.g. at the n-sigma is to call
1139   // the corresponding function GetNumberOfSigmas!
1140   //
1141   
1142   //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse
1143   
1144   Double_t dEdx = -1;
1145   Int_t nPoints = -1;
1146   ETPCgainScenario gainScenario = kGainScenarioInvalid;
1147   TSpline3* responseFunction = 0x0;
1148   
1149   // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case
1150   // it is not used anyway, so this causes no trouble.
1151   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
1152     return -1.;
1153   
1154   
1155   // No correction in case of multiplicity <= 0
1156   if (fCurrentEventMultiplicity <= 0)
1157     return dEdx;
1158   
1159   Double_t multiplicityCorr = 0.;
1160   
1161   // TODO Normally, one should use the eta corrected values in BOTH of the following cases. Does it make sense to use the uncorrected dEdx values?
1162   if (species < AliPID::kUnknown) {
1163     // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course).
1164     // However, one needs the eta corrected value!
1165     Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE);
1166     multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdxSplines, fCurrentEventMultiplicity);
1167   }
1168   else {
1169     // One needs the eta corrected value to determine the multiplicity correction factor!
1170     Double_t etaCorr = GetEtaCorrectionFast(track, dEdx);
1171     multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdx * etaCorr, fCurrentEventMultiplicity);
1172   }
1173     
1174   if (multiplicityCorr <= 0)
1175     return -1.;
1176   
1177   return dEdx / multiplicityCorr; 
1178 }
1179
1180
1181 //_________________________________________________________________________
1182 Double_t AliTPCPIDResponse::GetEtaAndMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species,
1183                                                                     ETPCdEdxSource dedxSource) const
1184 {
1185   //
1186   // Get multiplicity and eta corrected dEdx for the given track. For the correction,
1187   // the expected dEdx of the specified species will be used. If the species is set 
1188   // to AliPID::kUnknown, the dEdx of the track is used instead.
1189   // WARNING: In the latter case, the correction might not be as good as if the
1190   // expected dEdx is used, which is the way the correction factor is designed
1191   // for.
1192   // In any case, one has to decide carefully to which expected signal one wants to
1193   // compare the corrected value - to the corrected or uncorrected.
1194   // Anyhow, a safer way of looking e.g. at the n-sigma is to call
1195   // the corresponding function GetNumberOfSigmas!
1196   //
1197   
1198   Double_t dEdx = -1;
1199   Int_t nPoints = -1;
1200   ETPCgainScenario gainScenario = kGainScenarioInvalid;
1201   TSpline3* responseFunction = 0x0;
1202   
1203   // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case
1204   // it is not used anyway, so this causes no trouble.
1205   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
1206     return -1.;
1207   
1208   Double_t multiplicityCorr = 0.;
1209   Double_t etaCorr = 0.;
1210     
1211   if (species < AliPID::kUnknown) {
1212     // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course)
1213     Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE);
1214     etaCorr = GetEtaCorrectionFast(track, dEdxSplines);
1215     multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdxSplines * etaCorr, fCurrentEventMultiplicity);
1216   }
1217   else {
1218     etaCorr = GetEtaCorrectionFast(track, dEdx);
1219     multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdx * etaCorr, fCurrentEventMultiplicity);
1220   }
1221     
1222   if (multiplicityCorr <= 0 || etaCorr <= 0)
1223     return -1.;
1224   
1225   return dEdx / multiplicityCorr / etaCorr; 
1226 }
1227
1228
1229 //_________________________________________________________________________
1230 Double_t AliTPCPIDResponse::GetMultiplicitySigmaCorrectionFast(Double_t dEdxExpected, Int_t multiplicity) const
1231 {
1232   // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly.
1233   //
1234   // Calculate the multiplicity sigma correction factor for the corresponding expected dEdx and for the given multiplicity.
1235   // The parameter dEdxExpected should take into account the eta correction already!
1236   
1237   // Multiplicity dependence of sigma depends on the real dEdx at zero multiplicity,
1238   // i.e. the eta (only) corrected dEdxExpected value has to be used
1239   // since all maps etc. have been created for ~zero multiplicity
1240   
1241   if (dEdxExpected <= 0 || multiplicity <= 0)
1242     return 1.0;
1243   
1244   Double_t relSigmaSlope = fCorrFuncSigmaMultiplicity->Eval(1. / dEdxExpected);
1245  
1246   return (1. + relSigmaSlope * multiplicity);
1247 }
1248
1249
1250 //_________________________________________________________________________
1251 Double_t AliTPCPIDResponse::GetMultiplicitySigmaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
1252 {
1253   //
1254   // Get multiplicity sigma correction for the given track (w.r.t. the mulitplicity of the current event)
1255   //
1256   
1257   //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse
1258   
1259   // No correction in case of multiplicity <= 0
1260   if (fCurrentEventMultiplicity <= 0)
1261     return 1.;
1262   
1263   Double_t dEdx = -1;
1264   Int_t nPoints = -1;
1265   ETPCgainScenario gainScenario = kGainScenarioInvalid;
1266   TSpline3* responseFunction = 0x0;
1267   
1268   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
1269     return 1.; 
1270   
1271   //TODO Does it make sense to use the multiplicity correction WITHOUT eta correction?! Are the fit parameters still valid?
1272   // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course)
1273   Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE);
1274   
1275   return GetMultiplicitySigmaCorrectionFast(dEdxExpected, fCurrentEventMultiplicity);
1276 }
1277
1278
1279 //_________________________________________________________________________
1280 void AliTPCPIDResponse::ResetMultiplicityCorrectionFunctions()
1281 {
1282   // Default values: No correction, i.e. overall correction factor should be one
1283   
1284   // This function should always return 0.0, if multcorr disabled
1285   fCorrFuncMultiplicity->SetParameter(0, 0.);
1286   fCorrFuncMultiplicity->SetParameter(1, 0.);
1287   fCorrFuncMultiplicity->SetParameter(2, 0.);
1288   fCorrFuncMultiplicity->SetParameter(3, 0.);
1289   fCorrFuncMultiplicity->SetParameter(4, 0.);
1290     
1291   // This function should always return 0., if multCorr disabled
1292   fCorrFuncMultiplicityTanTheta->SetParameter(0, 0.);
1293   fCorrFuncMultiplicityTanTheta->SetParameter(1, 0.);
1294   fCorrFuncMultiplicityTanTheta->SetParameter(2, 0.);
1295   
1296   // This function should always return 0.0, if mutlcorr disabled
1297   fCorrFuncSigmaMultiplicity->SetParameter(0, 0.);
1298   fCorrFuncSigmaMultiplicity->SetParameter(1, 0.);
1299   fCorrFuncSigmaMultiplicity->SetParameter(2, 0.);
1300   fCorrFuncSigmaMultiplicity->SetParameter(3, 0.);
1301 }
1302
1303
1304 //_________________________________________________________________________
1305 Bool_t AliTPCPIDResponse::sectorNumbersInOut(Double_t* trackPositionInner,
1306                                              Double_t* trackPositionOuter,
1307                                              Float_t& inphi,
1308                                              Float_t& outphi,
1309                                              Int_t& in,
1310                                              Int_t& out ) const
1311 {
1312   //calculate the sector numbers (equivalent to IROC chamber numbers) a track crosses
1313   //for OROC chamber numbers add 36
1314   //returned angles are between (0,2pi)
1315
1316   inphi = TMath::ATan2(trackPositionInner[1],trackPositionInner[0]);
1317   outphi = TMath::ATan2(trackPositionOuter[1], trackPositionOuter[0]);
1318
1319   if (inphi<0) {inphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi
1320   if (outphi<0) {outphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi
1321
1322   in = sectorNumber(inphi);
1323   out = sectorNumber(outphi);
1324
1325   //for the C side (positive z) offset by 18
1326   if (trackPositionInner[2]>0.0)
1327   {
1328     in+=18;
1329     out+=18;
1330   }
1331   return kTRUE;
1332 }
1333
1334
1335 //_____________________________________________________________________________
1336 Int_t AliTPCPIDResponse::sectorNumber(Double_t phi) const
1337 {
1338   //calculate sector number
1339   const Float_t width=TMath::TwoPi()/18.;
1340   return TMath::Floor(phi/width);
1341 }
1342
1343
1344 //_____________________________________________________________________________
1345 void AliTPCPIDResponse::Print(Option_t* /*option*/) const
1346 {
1347   //Print info
1348   fResponseFunctions.Print();
1349 }
1350
1351
1352 //_____________________________________________________________________________
1353 AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack* track, Int_t layer) const
1354 {
1355   //status of the track: if it crosses any bad chambers return kChamberOff;
1356   //IROC:layer=1, OROC:layer=2
1357   if (layer<1 || layer>2) layer=1;
1358   Int_t in=0;
1359   Int_t out=0;
1360   Float_t inphi=0.;
1361   Float_t outphi=0.;
1362   Float_t innerRadius = (layer==1)?83.0:133.7;
1363   Float_t outerRadius = (layer==1)?133.5:247.7;
1364
1365   /////////////////////////////////////////////////////////////////////////////
1366   //find out where track enters and leaves the layer.
1367   //
1368   Double_t trackPositionInner[3];
1369   Double_t trackPositionOuter[3];
1370  
1371   //if there is no inner param this could mean we're using the AOD track,
1372   //we still can extrapolate from the vertex - so use those params.
1373   const AliExternalTrackParam* ip = track->GetInnerParam();
1374   if (ip) track=ip;
1375
1376   Bool_t trackAtInner = track->GetXYZAt(innerRadius, fMagField, trackPositionInner);
1377   Bool_t trackAtOuter = track->GetXYZAt(outerRadius, fMagField, trackPositionOuter);
1378
1379   if (!trackAtInner)
1380   {
1381     //if we dont even enter inner radius we do nothing and return invalid
1382     inphi=0.0;
1383     outphi=0.0;
1384     in=0;
1385     out=0;
1386     return kChamberInvalid;
1387   }
1388
1389   if (!trackAtOuter)
1390   {
1391     //if we don't reach the outer radius check that the apex is indeed within the outer radius and use apex position
1392     Bool_t haveApex = TrackApex(track, fMagField, trackPositionOuter);
1393     Float_t apexRadius = TMath::Sqrt(trackPositionOuter[0]*trackPositionOuter[0]+trackPositionOuter[1]*trackPositionOuter[1]);
1394     if ( haveApex && apexRadius<=outerRadius && apexRadius>innerRadius)
1395     {
1396       //printf("pt: %.2f, apexRadius: %.2f(%s), x: %.2f, y: %.2f\n",track->Pt(),apexRadius,(haveApex)?"OK":"BAD",trackPositionOuter[0],trackPositionOuter[1]);
1397     }
1398     else
1399     {
1400       inphi=0.0;
1401       outphi=0.0;
1402       in=0;
1403       out=0;
1404       return kChamberInvalid;
1405     }
1406   }
1407
1408
1409   if (!sectorNumbersInOut(trackPositionInner,
1410                           trackPositionOuter,
1411                           inphi,
1412                           outphi,
1413                           in,
1414                           out)) return kChamberInvalid;
1415
1416   /////////////////////////////////////////////////////////////////////////////
1417   //now we have the location of the track we can check
1418   //if it is in a good/bad chamber
1419   //
1420   Bool_t sideA = kTRUE;
1421  
1422   if (((in/18)%2==1) && ((out/18)%2==1)) sideA=kFALSE;
1423
1424   in=in%18;
1425   out=out%18;
1426
1427   if (TMath::Abs(in-out)>9)
1428   {
1429     if (TMath::Max(in,out)==out)
1430     {
1431       Int_t tmp=out;
1432       out = in;
1433       in = tmp;
1434       Float_t tmpphi=outphi;
1435       outphi=inphi;
1436       inphi=tmpphi;
1437     }
1438     in-=18;
1439     inphi-=TMath::TwoPi();
1440   }
1441   else
1442   {
1443     if (TMath::Max(in,out)==in)
1444     {
1445       Int_t tmp=out;
1446       out=in;
1447       in=tmp;
1448       Float_t tmpphi=outphi;
1449       outphi=inphi;
1450       inphi=tmpphi;
1451     }
1452   }
1453
1454   Float_t trackLengthInBad = 0.;
1455   Float_t trackLengthInLowGain = 0.;
1456   Float_t trackLengthTotal = TMath::Abs(outphi-inphi);
1457   Float_t lengthFractionInBadSectors = 0.;
1458
1459   const Float_t sectorWidth = TMath::TwoPi()/18.;  
1460  
1461   for (Int_t i=in; i<=out; i++)
1462   {
1463     int j=i;
1464     if (i<0) j+=18;    //correct for the negative values
1465     if (!sideA) j+=18; //move to the correct side
1466    
1467     Float_t deltaPhi = 0.;
1468     Float_t phiEdge=sectorWidth*i;
1469     if (inphi>phiEdge) {deltaPhi=phiEdge+sectorWidth-inphi;}
1470     else if ((outphi>=phiEdge) && (outphi<(phiEdge+sectorWidth))) {deltaPhi=outphi-phiEdge;}
1471     else {deltaPhi=sectorWidth;}
1472    
1473     Float_t v = fVoltageMap[(layer==1)?(j):(j+36)];
1474     if (v<=fBadIROCthreshhold)
1475     {
1476       trackLengthInBad+=deltaPhi;
1477       lengthFractionInBadSectors=1.;
1478     }
1479     if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold)
1480     {
1481       trackLengthInLowGain+=deltaPhi;
1482       lengthFractionInBadSectors=1.;
1483     }
1484   }
1485
1486   //for now low gain and bad (off) chambers are treated equally
1487   if (trackLengthTotal>0)
1488     lengthFractionInBadSectors = (trackLengthInLowGain+trackLengthInBad)/trackLengthTotal;
1489
1490   //printf("### side: %s, pt: %.2f, pz: %.2f, in: %i, out: %i, phiIN: %.2f, phiOUT: %.2f, rIN: %.2f, rOUT: %.2f\n",(sideA)?"A":"C",track->Pt(),track->Pz(),in,out,inphi,outphi,innerRadius,outerRadius);
1491  
1492   if (lengthFractionInBadSectors>fMaxBadLengthFraction)
1493   {
1494     //printf("%%%%%%%% %s kChamberLowGain\n",(layer==1)?"IROC":"OROC");
1495     return kChamberLowGain;
1496   }
1497  
1498   return kChamberHighGain;
1499 }
1500
1501
1502 //_____________________________________________________________________________
1503 Float_t AliTPCPIDResponse::MaxClusterRadius(const AliVTrack* track) const
1504 {
1505   //return the radius of the outermost padrow containing a cluster in TPC
1506   //for the track
1507   const TBits* clusterMap=track->GetTPCClusterMapPtr();
1508   if (!clusterMap) return 0.;
1509
1510   //from AliTPCParam, radius of first IROC row
1511   const Float_t rfirstIROC = 8.52249984741210938e+01;
1512   const Float_t padrowHeightIROC = 0.75;
1513   const Float_t rfirstOROC0 = 1.35100006103515625e+02;
1514   const Float_t padrowHeightOROC0 = 1.0;
1515   const Float_t rfirstOROC1 = 1.99350006103515625e+02;
1516   const Float_t padrowHeightOROC1 = 1.5;
1517
1518   Int_t maxPadRow=160;
1519   while ((--maxPadRow)>0 && !clusterMap->TestBitNumber(maxPadRow)){}
1520   if (maxPadRow>126) return rfirstOROC1+(maxPadRow-126-1)*padrowHeightOROC1;
1521   if (maxPadRow>62) return rfirstOROC0+(maxPadRow-62-1)*padrowHeightOROC0;
1522   if (maxPadRow>0) return rfirstIROC+(maxPadRow-1)*padrowHeightIROC;
1523   return 0.0;
1524 }
1525
1526
1527 //_____________________________________________________________________________
1528 Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const
1529 {
1530   //calculate the coordinates of the apex of the track
1531   Double_t x[3];
1532   track->GetXYZ(x);
1533   Double_t p[3];
1534   track->GetPxPyPz(p);
1535   Double_t r = 1./track->OneOverPt()/0.0299792458/magField; //signed - will determine the direction of b
1536   //printf("b: %.2f, x:%.2f, y:%.2f, pt: %.2f, px:%.2f, py%.2f, r: %.2f\n",magField, x[0],x[1],track->Pt(), p[0],p[1],r);
1537   //find orthogonal vector (Gram-Schmidt)
1538   Double_t alpha = (p[0]*x[0] + p[1]*x[1])/(p[0]*p[0] + p[1]*p[1]);
1539   Double_t b[2];
1540   b[0]=x[0]-alpha*p[0];
1541   b[1]=x[1]-alpha*p[1];
1542  
1543   Double_t norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
1544   if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
1545   b[0]/=norm;
1546   b[1]/=norm;
1547   b[0]*=r;
1548   b[1]*=r;
1549   b[0]+=x[0];
1550   b[1]+=x[1];
1551   //printf("center: x:%.2f, y:%.2f\n",b[0],b[1]);
1552  
1553   norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
1554   if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
1555  
1556   position[0]=b[0]+b[0]*TMath::Abs(r)/norm;
1557   position[1]=b[1]+b[1]*TMath::Abs(r)/norm;
1558   position[2]=0.;
1559   return kTRUE;
1560 }
1561
1562 Double_t AliTPCPIDResponse::EvaldEdxSpline(Double_t bg,Int_t entry){
1563   //
1564   // Evaluate the dEdx response for given entry
1565   //
1566   TSpline * spline = (TSpline*)fSplineArray.At(entry);
1567   if (spline) return spline->Eval(bg);
1568   return 0;
1569 }
1570
1571
1572 Bool_t   AliTPCPIDResponse::RegisterSpline(const char * name, Int_t index){
1573   //
1574   // register spline to be used for drawing comparisons
1575   // 
1576   TFile * fTPCBB = TFile::Open("$ALICE_ROOT/OADB/COMMON/PID/data/TPCPIDResponse.root");
1577   TObjArray  *arrayTPCPID= (TObjArray*)  fTPCBB->Get("TPCPIDResponse");
1578   if (fSplineArray.GetEntriesFast()<index) fSplineArray.Expand(index*2);
1579   TSpline3 *spline=0;
1580   if (arrayTPCPID){
1581     spline = (TSpline3*)arrayTPCPID->FindObject(name);
1582     if (spline) fSplineArray.AddAt(spline->Clone(),index);    
1583   }
1584   delete arrayTPCPID;
1585   delete fTPCBB;
1586   return (spline!=0);
1587 }