Here the comment:
[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
42 ClassImp(AliTPCPIDResponse);
43
44 const char* AliTPCPIDResponse::fgkGainScenarioName[fgkNumberOfGainScenarios+1]=
45 {
46   "", //default - no name
47   "1", //all high
48   "2", //OROC only
49   "unknown"
50 };
51
52 //_________________________________________________________________________
53 AliTPCPIDResponse::AliTPCPIDResponse():
54   TNamed(),
55   fMIP(50.),
56   fRes0(),
57   fResN2(),
58   fKp1(0.0283086),
59   fKp2(2.63394e+01),
60   fKp3(5.04114e-11),
61   fKp4(2.12543),
62   fKp5(4.88663),
63   fUseDatabase(kFALSE),
64   fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
65   fVoltageMap(72),
66   fLowGainIROCthreshold(-40),
67   fBadIROCthreshhold(-70),
68   fLowGainOROCthreshold(-40),
69   fBadOROCthreshhold(-40),
70   fMaxBadLengthFraction(0.5),
71   fMagField(0.),
72   fhEtaCorr(0x0),
73   fhEtaSigmaPar1(0x0),
74   fSigmaPar0(0.0)
75 {
76   //
77   //  The default constructor
78   //
79   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=0.07;fResN2[i]=0.0;}
80 }
81 /*TODO remove?
82 //_________________________________________________________________________
83 AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param):
84   TNamed(),
85   fMIP(param[0]),
86   fRes0(),
87   fResN2(),
88   fKp1(0.0283086),
89   fKp2(2.63394e+01),
90   fKp3(5.04114e-11),
91   fKp4(2.12543),
92   fKp5(4.88663),
93   fUseDatabase(kFALSE),
94   fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
95   fVoltageMap(72),
96   fLowGainIROCthreshold(-40),
97   fBadIROCthreshhold(-70),
98   fLowGainOROCthreshold(-40),
99   fBadOROCthreshhold(-40),
100   fMaxBadLengthFraction(0.5),
101   fMagField(0.),
102   fhEtaCorr(0x0),
103   fhEtaSigmaPar1(0x0),
104   fSigmaPar0(0.0)
105 {
106   //
107   //  The main constructor
108   //
109   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=param[1];fResN2[i]=param[2];}
110 }
111 */
112
113 //_________________________________________________________________________
114 AliTPCPIDResponse::~AliTPCPIDResponse()
115 {
116   //
117   // Destructor
118   //
119   
120   delete fhEtaCorr;
121   fhEtaCorr = 0x0;
122   
123   delete fhEtaSigmaPar1;
124   fhEtaSigmaPar1 = 0x0;
125 }
126
127
128 //_________________________________________________________________________
129 AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that):
130   TNamed(that),
131   fMIP(that.fMIP),
132   fRes0(),
133   fResN2(),
134   fKp1(that.fKp1),
135   fKp2(that.fKp2),
136   fKp3(that.fKp3),
137   fKp4(that.fKp4),
138   fKp5(that.fKp5),
139   fUseDatabase(that.fUseDatabase),
140   fResponseFunctions(that.fResponseFunctions),
141   fVoltageMap(that.fVoltageMap),
142   fLowGainIROCthreshold(that.fLowGainIROCthreshold),
143   fBadIROCthreshhold(that.fBadIROCthreshhold),
144   fLowGainOROCthreshold(that.fLowGainOROCthreshold),
145   fBadOROCthreshhold(that.fBadOROCthreshhold),
146   fMaxBadLengthFraction(that.fMaxBadLengthFraction),
147   fMagField(that.fMagField),
148   fhEtaCorr(0x0),
149   fhEtaSigmaPar1(0x0),
150   fSigmaPar0(that.fSigmaPar0)
151 {
152   //copy ctor
153   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
154  
155   // Copy eta maps
156   if (that.fhEtaCorr) {
157     fhEtaCorr = new TH2D(*(that.fhEtaCorr));
158     fhEtaCorr->SetDirectory(0);
159   }
160   
161   if (that.fhEtaSigmaPar1) {
162     fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1));
163     fhEtaSigmaPar1->SetDirectory(0);
164   }
165 }
166
167 //_________________________________________________________________________
168 AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that)
169 {
170   //assignment
171   if (&that==this) return *this;
172   TNamed::operator=(that);
173   fMIP=that.fMIP;
174   fKp1=that.fKp1;
175   fKp2=that.fKp2;
176   fKp3=that.fKp3;
177   fKp4=that.fKp4;
178   fKp5=that.fKp5;
179   fUseDatabase=that.fUseDatabase;
180   fResponseFunctions=that.fResponseFunctions;
181   fVoltageMap=that.fVoltageMap;
182   fLowGainIROCthreshold=that.fLowGainIROCthreshold;
183   fBadIROCthreshhold=that.fBadIROCthreshhold;
184   fLowGainOROCthreshold=that.fLowGainOROCthreshold;
185   fBadOROCthreshhold=that.fBadOROCthreshhold;
186   fMaxBadLengthFraction=that.fMaxBadLengthFraction;
187   fMagField=that.fMagField;
188   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
189
190   delete fhEtaCorr;
191   fhEtaCorr=0x0;
192   if (that.fhEtaCorr) {
193     fhEtaCorr = new TH2D(*(that.fhEtaCorr));
194     fhEtaCorr->SetDirectory(0);
195   }
196   
197   delete fhEtaSigmaPar1;
198   fhEtaSigmaPar1=0x0;
199   if (that.fhEtaSigmaPar1) {
200     fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1));
201     fhEtaSigmaPar1->SetDirectory(0);
202   }
203   
204   fSigmaPar0 = that.fSigmaPar0;
205
206   return *this;
207 }
208
209 //_________________________________________________________________________
210 Double_t AliTPCPIDResponse::Bethe(Double_t betaGamma) const {
211   //
212   // This is the Bethe-Bloch function normalised to 1 at the minimum
213   // WARNING
214   // Simulated and reconstructed Bethe-Bloch differs
215   //           Simulated  curve is the dNprim/dx
216   //           Reconstructed is proportianal dNtot/dx
217   // Temporary fix for production -  Simple linear correction function
218   // Future    2 Bethe Bloch formulas needed
219   //           1. for simulation
220   //           2. for reconstructed PID
221   //
222   
223 //   const Float_t kmeanCorrection =0.1;
224   Double_t bb=
225     AliExternalTrackParam::BetheBlochAleph(betaGamma,fKp1,fKp2,fKp3,fKp4,fKp5);
226   return bb*fMIP;
227 }
228
229 //_________________________________________________________________________
230 void AliTPCPIDResponse::SetBetheBlochParameters(Double_t kp1,
231                              Double_t kp2,
232                              Double_t kp3,
233                              Double_t kp4,
234                              Double_t kp5) {
235   //
236   // Set the parameters of the ALEPH Bethe-Bloch formula
237   //
238   fKp1=kp1;
239   fKp2=kp2;
240   fKp3=kp3;
241   fKp4=kp4;
242   fKp5=kp5;
243 }
244
245 //_________________________________________________________________________
246 void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2) {
247   //
248   // Set the relative resolution  sigma_rel = res0 * sqrt(1+resN2/npoint)
249   //
250   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=res0;fResN2[i]=resN2;}
251 }
252
253 //_________________________________________________________________________
254 Double_t AliTPCPIDResponse::GetExpectedSignal(const Float_t mom,
255                                               AliPID::EParticleType n) const {
256   //
257   // Deprecated function (for backward compatibility). Please use 
258   // GetExpectedSignal(const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource,
259   //                   Bool_t correctEta = kTRUE);
260   // instead!
261   //
262   //
263   // Calculates the expected PID signal as the function of 
264   // the information stored in the track, for the specified particle type 
265   //  
266   // At the moment, these signals are just the results of calling the 
267   // Bethe-Bloch formula. 
268   // This can be improved. By taking into account the number of
269   // assigned clusters and/or the track dip angle, for example.  
270   //
271   
272   //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
273   // !!! Splines for light nuclei need to be normalised to this factor !!!
274   const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(n),2.3);
275   
276   Double_t mass=AliPID::ParticleMassZ(n);
277   if (!fUseDatabase) return Bethe(mom/mass) * chargeFactor;
278   //
279   const TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n);
280
281   if (!responseFunction) return Bethe(mom/mass) * chargeFactor;
282   
283   return fMIP*responseFunction->Eval(mom/mass)*chargeFactor;
284
285 }
286
287 //_________________________________________________________________________
288 Double_t AliTPCPIDResponse::GetExpectedSigma(const Float_t mom, 
289                                              const Int_t nPoints,
290                                              AliPID::EParticleType n) const {
291   //
292   // Deprecated function (for backward compatibility). Please use 
293   // GetExpectedSigma(onst AliVTrack* track, AliPID::EParticleType species, 
294   // ETPCdEdxSource dedxSource, Bool_t correctEta) instead!
295   //
296   //
297   // Calculates the expected sigma of the PID signal as the function of 
298   // the information stored in the track, for the specified particle type 
299   //  
300   
301   if (nPoints != 0) 
302     return GetExpectedSignal(mom,n)*fRes0[0]*sqrt(1. + fResN2[0]/nPoints);
303   else
304     return GetExpectedSignal(mom,n)*fRes0[0];
305 }
306
307 ////////////////////////////////////////////////////NEW//////////////////////////////
308
309 //_________________________________________________________________________
310 void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario) {
311   //
312   // Set the relative resolution  sigma_rel = res0 * sqrt(1+resN2/npoint)
313   //
314   if ((Int_t)gainScenario>(Int_t)fgkNumberOfGainScenarios) return; //TODO: better handling!
315   fRes0[gainScenario]=res0;
316   fResN2[gainScenario]=resN2;
317 }
318
319
320 //_________________________________________________________________________
321 Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
322                                               AliPID::EParticleType species,
323                                               Double_t /*dEdx*/,
324                                               const TSpline3* responseFunction,
325                                               Bool_t correctEta) const 
326 {
327   // Calculates the expected PID signal as the function of 
328   // the information stored in the track and the given parameters,
329   // for the specified particle type 
330   //  
331   // At the moment, these signals are just the results of calling the 
332   // Bethe-Bloch formula plus, if desired, taking into account the eta dependence. 
333   // This can be improved. By taking into account the number of
334   // assigned clusters and/or the track dip angle, for example.  
335   //
336   
337   Double_t mom=track->GetTPCmomentum();
338   Double_t mass=AliPID::ParticleMass(species);
339   
340   //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
341   // !!! Splines for light nuclei need to be normalised to this factor !!!
342   const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3);
343   
344   if (!responseFunction)
345     return Bethe(mom/mass) * chargeFactor;
346   
347   Double_t dEdxSplines = fMIP*responseFunction->Eval(mom/mass) * chargeFactor;
348   
349   if (!correctEta)
350     return dEdxSplines;
351   
352   //TODO Alternatively take current track dEdx
353   //return dEdxSplines * GetEtaCorrection(track, dEdx);
354   return dEdxSplines * GetEtaCorrection(track, dEdxSplines);
355 }
356
357
358 //_________________________________________________________________________
359 Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
360                                               AliPID::EParticleType species,
361                                               ETPCdEdxSource dedxSource,
362                                               Bool_t correctEta) const
363 {
364   // Calculates the expected PID signal as the function of 
365   // the information stored in the track, for the specified particle type 
366   //  
367   // At the moment, these signals are just the results of calling the 
368   // Bethe-Bloch formula plus, if desired, taking into account the eta dependence. 
369   // This can be improved. By taking into account the number of
370   // assigned clusters and/or the track dip angle, for example.  
371   //
372   
373   if (!fUseDatabase) {
374     //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
375     // !!! Splines for light nuclei need to be normalised to this factor !!!
376     const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3);
377   
378     return Bethe(track->GetTPCmomentum() / AliPID::ParticleMass(species)) * chargeFactor;
379   }
380   
381   Double_t dEdx = -1;
382   Int_t nPoints = -1;
383   ETPCgainScenario gainScenario = kGainScenarioInvalid;
384   TSpline3* responseFunction = 0x0;
385     
386   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) {
387     // Something is wrong with the track -> Return obviously invalid value
388     return -999;
389   }
390   
391   // Charge factor already taken into account inside the following function call
392   return GetExpectedSignal(track, species, dEdx, responseFunction, correctEta);
393 }
394   
395 //_________________________________________________________________________
396 TSpline3* AliTPCPIDResponse::GetResponseFunction( AliPID::EParticleType type,
397                                                   AliTPCPIDResponse::ETPCgainScenario gainScenario ) const
398 {
399   //get response function
400   return dynamic_cast<TSpline3*>(fResponseFunctions.At(ResponseFunctionIndex(type,gainScenario)));
401 }
402
403 //_________________________________________________________________________
404 TSpline3* AliTPCPIDResponse::GetResponseFunction( const AliVTrack* track,
405                                AliPID::EParticleType species,
406                                ETPCdEdxSource dedxSource) const 
407 {
408   //the splines are stored in an array, different scenarios
409
410   Double_t dEdx = -1;
411   Int_t nPoints = -1;
412   ETPCgainScenario gainScenario = kGainScenarioInvalid;
413   TSpline3* responseFunction = 0x0;
414   
415   if (ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
416     return responseFunction;
417   
418   return NULL;
419 }
420
421 //_________________________________________________________________________
422 void AliTPCPIDResponse::ResetSplines()
423 {
424   //reset the array with splines
425   for (Int_t i=0;i<fResponseFunctions.GetEntriesFast();i++)
426   {
427     fResponseFunctions.AddAt(NULL,i);
428   }
429 }
430 //_________________________________________________________________________
431 Int_t AliTPCPIDResponse::ResponseFunctionIndex( AliPID::EParticleType species,
432                                                 ETPCgainScenario gainScenario ) const
433 {
434   //get the index in fResponseFunctions given type and scenario
435   return Int_t(species)+Int_t(gainScenario)*fgkNumberOfParticleSpecies;
436 }
437
438 //_________________________________________________________________________
439 void AliTPCPIDResponse::SetResponseFunction( TObject* o,
440                                              AliPID::EParticleType species,
441                                              ETPCgainScenario gainScenario )
442 {
443   fResponseFunctions.AddAtAndExpand(o,ResponseFunctionIndex(species,gainScenario));
444 }
445
446
447 //_________________________________________________________________________
448 Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, 
449                                              AliPID::EParticleType species,
450                                              ETPCgainScenario gainScenario,
451                                              Double_t dEdx,
452                                              Int_t nPoints,
453                                              const TSpline3* responseFunction,
454                                              Bool_t correctEta) const 
455 {
456   // Calculates the expected sigma of the PID signal as the function of 
457   // the information stored in the track and the given parameters,
458   // for the specified particle type 
459   //
460   
461   if (!responseFunction)
462     return 999;
463   
464   // If no sigma map is available or if no eta correction is requested (sigma maps only for corrected eta!), use the old parametrisation
465   if (!fhEtaSigmaPar1 || !correctEta) {  
466     if (nPoints != 0) 
467       return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE) *
468                fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints);
469     else
470       return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE)*fRes0[gainScenario];
471   }
472     
473   if (nPoints > 0) {
474     Double_t sigmaPar1 = GetSigmaPar1(track, species, dEdx, responseFunction);
475     
476     return GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE)*sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints);
477   }
478   else { 
479     // One should never have/take tracks with 0 dEdx clusters!
480     return 999;
481   }
482 }
483
484
485 //_________________________________________________________________________
486 Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, 
487                                              AliPID::EParticleType species,
488                                              ETPCdEdxSource dedxSource,
489                                              Bool_t correctEta) const 
490 {
491   // Calculates the expected sigma of the PID signal as the function of 
492   // the information stored in the track, for the specified particle type 
493   // and dedx scenario
494   //
495   
496   Double_t dEdx = -1;
497   Int_t nPoints = -1;
498   ETPCgainScenario gainScenario = kGainScenarioInvalid;
499   TSpline3* responseFunction = 0x0;
500   
501   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
502     return 999; //TODO: better handling!
503   
504   return GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta);
505 }
506
507
508 //_________________________________________________________________________
509 Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track, 
510                              AliPID::EParticleType species,
511                              ETPCdEdxSource dedxSource,
512                              Bool_t correctEta) const
513 {
514   //Calculates the number of sigmas of the PID signal from the expected value
515   //for a given particle species in the presence of multiple gain scenarios
516   //inside the TPC
517                              
518   Double_t dEdx = -1;
519   Int_t nPoints = -1;
520   ETPCgainScenario gainScenario = kGainScenarioInvalid;
521   TSpline3* responseFunction = 0x0;
522   
523   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
524     return -999; //TODO: Better handling!
525                              
526   Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta);
527   Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta);
528   // 999 will be returned by GetExpectedSigma e.g. in case of 0 dEdx clusters
529   if (sigma >= 998) 
530     return -999;
531   else
532     return (dEdx-bethe)/sigma;
533 }
534
535 //_________________________________________________________________________
536 Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track, 
537                                                  AliPID::EParticleType species,
538                                                  ETPCdEdxSource dedxSource,
539                                                  Double_t& dEdx,
540                                                  Int_t& nPoints,
541                                                  ETPCgainScenario& gainScenario,
542                                                  TSpline3** responseFunction) const 
543 {
544   // Calculates the right parameters for PID
545   //   dEdx parametrization for the proper gain scenario, dEdx 
546   //   and NPoints used for dEdx
547   // based on the track geometry (which chambers it crosses) for the specified particle type 
548   // and preferred source of dedx.
549   // returns true on success
550   
551   
552   if (dedxSource == kdEdxDefault) {
553     // Fast handling for default case. In addition: Keep it simple (don't call additional functions) to
554     // avoid possible bugs
555     
556     // GetTPCsignalTunedOnData will be non-positive, if it has not been set (i.e. in case of MC NOT tuned to data).
557     // If this is the case, just take the normal signal
558     dEdx = track->GetTPCsignalTunedOnData();
559     if (dEdx <= 0) {
560       dEdx = track->GetTPCsignal();
561     }
562     
563     nPoints = track->GetTPCsignalN();
564     gainScenario = kDefault;
565     
566     TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario));
567     *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast?
568   
569     return kTRUE;
570   }
571   
572   //TODO Proper handle of tuneMConData for other dEdx sources
573   
574   Double32_t signal[4]; //0: IROC, 1: OROC medium, 2:OROC long, 3: OROC all (def. truncation used)
575   Char_t ncl[3];        //same
576   Char_t nrows[3];      //same
577   const AliTPCdEdxInfo* dEdxInfo = track->GetTPCdEdxInfo();
578   
579   if (!dEdxInfo && dedxSource!=kdEdxDefault)  //in one case its ok if we dont have the info
580   {
581     AliError("AliTPCdEdxInfo not available");
582     return kFALSE;
583   }
584
585   if (dEdxInfo) dEdxInfo->GetTPCSignalRegionInfo(signal,ncl,nrows);
586
587   //check if we cross a bad OROC in which case we reject
588   EChamberStatus trackOROCStatus = TrackStatus(track,2);
589   if (trackOROCStatus==kChamberOff || trackOROCStatus==kChamberLowGain)
590   {
591     return kFALSE;
592   }
593
594   switch (dedxSource)
595   {
596     case kdEdxOROC:
597       {
598         if (trackOROCStatus==kChamberInvalid) return kFALSE; //never reached OROC
599         dEdx = signal[3];
600         nPoints = ncl[2]+ncl[1];
601         gainScenario = kOROChigh;
602         break;
603       }
604     case kdEdxHybrid:
605       {
606         //if we cross a bad IROC we use OROC dedx, if we dont we use combined
607         EChamberStatus status = TrackStatus(track,1);
608         if (status!=kChamberHighGain)
609         {
610           dEdx = signal[3];
611           nPoints = ncl[2]+ncl[1];
612           gainScenario = kOROChigh;
613         }
614         else
615         {
616           dEdx = track->GetTPCsignal();
617           nPoints = track->GetTPCsignalN();
618           gainScenario = kALLhigh;
619         }
620         break;
621       }
622     default:
623       {
624          dEdx = 0.;
625          nPoints = 0;
626          gainScenario = kGainScenarioInvalid;
627          return kFALSE;
628       }
629   }
630   TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario));
631   *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast?
632   
633   return kTRUE;
634 }
635
636
637 //_________________________________________________________________________
638 Double_t AliTPCPIDResponse::GetEtaCorrection(const AliVTrack *track, Double_t dEdxSplines) const
639 {
640   //
641   // Get eta correction for the given parameters.
642   //
643   
644   if (!fhEtaCorr)
645     return 1.;
646   
647   Double_t tpcSignal = dEdxSplines;
648   
649   if (tpcSignal < 1.)
650     return 1.;
651   
652   // For ESD tracks, the local tanTheta could be used (esdTrack->GetInnerParam()->GetTgl()).
653   // However, this value is not available for AODs and, thus, not for AliVTrack.
654   // Fortunately, the following formula allows to approximate the local tanTheta with the 
655   // global theta angle -> This is for by far most of the tracks the same, but gives at
656   // maybe the percent level differences within +- 0.2 in tanTheta -> Which is still ok.
657   Double_t tanTheta = TMath::Tan(-track->Theta() + TMath::Pi() / 2.0);
658   Int_t binX = fhEtaCorr->GetXaxis()->FindBin(tanTheta);
659   Int_t binY = fhEtaCorr->GetYaxis()->FindBin(1. / tpcSignal);
660   
661   if (binX == 0) 
662     binX = 1;
663   if (binX > fhEtaCorr->GetXaxis()->GetNbins())
664     binX = fhEtaCorr->GetXaxis()->GetNbins();
665   
666   if (binY == 0)
667     binY = 1;
668   if (binY > fhEtaCorr->GetYaxis()->GetNbins())
669     binY = fhEtaCorr->GetYaxis()->GetNbins();
670   
671   return fhEtaCorr->GetBinContent(binX, binY);
672 }
673
674
675 //_________________________________________________________________________
676 Double_t AliTPCPIDResponse::GetEtaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
677 {
678   //
679   // Get eta correction for the given track.
680   //
681   
682   if (!fhEtaCorr)
683     return 1.;
684   
685   Double_t dEdx = -1;
686   Int_t nPoints = -1;
687   ETPCgainScenario gainScenario = kGainScenarioInvalid;
688   TSpline3* responseFunction = 0x0;
689   
690   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
691     return 1.; 
692   
693   Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE);
694   
695   //TODO Alternatively take current track dEdx
696   //return GetEtaCorrection(track, dEdx);
697   
698   return GetEtaCorrection(track, dEdxSplines);
699 }
700
701
702 //_________________________________________________________________________
703 Double_t AliTPCPIDResponse::GetEtaCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
704 {
705   //
706   // Get eta corrected dEdx for the given track. For the correction, the expected dEdx of
707   // the specified species will be used. If the species is set to AliPID::kUnknown, the
708   // dEdx of the track is used instead.
709   // WARNING: In the latter case, the eta correction might not be as good as if the
710   // expected dEdx is used, which is the way the correction factor is designed
711   // for.
712   // In any case, one has to decide carefully to which expected signal one wants to
713   // compare the corrected value - to the corrected or uncorrected.
714   // Anyhow, a safer way of looking e.g. at the n-sigma is to call
715   // the corresponding function GetNumberOfSigmas!
716   //
717   
718   Double_t dEdx = -1;
719   Int_t nPoints = -1;
720   ETPCgainScenario gainScenario = kGainScenarioInvalid;
721   TSpline3* responseFunction = 0x0;
722   
723   // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case
724   // it is not used anyway, so this causes no trouble.
725   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
726     return -1.;
727   
728   Double_t etaCorr = 0.;
729   
730   if (species < AliPID::kUnknown) {
731     Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE);
732     etaCorr = GetEtaCorrection(track, dEdxSplines);
733   }
734   else {
735     etaCorr = GetEtaCorrection(track, dEdx);
736   }
737     
738   if (etaCorr <= 0)
739     return -1.;
740   
741   return dEdx / etaCorr; 
742 }
743
744
745
746 //_________________________________________________________________________
747 Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, Double_t dEdx, const TSpline3* responseFunction) const
748 {
749   //
750   // Get parameter 1 of sigma parametrisation of TPC dEdx from the histogram for the given track.
751   //
752   
753   if (!fhEtaSigmaPar1)
754     return 999;
755   
756   // The sigma maps are created with data sets that are already eta corrected and for which the 
757   // splines have been re-created. Therefore, the value for the lookup needs to be the value of
758   // the splines without any additional eta correction.
759   // NOTE: This is due to the method the maps are created. The track dEdx (not the expected one!)
760   // is corrected to uniquely related a momemtum bin with an expected dEdx, where the expected dEdx
761   // equals the track dEdx for all eta (thanks to the correction and the re-fit of the splines).
762   // Consequently, looking up the uncorrected expected dEdx at a given tanTheta yields the correct
763   // sigma parameter!
764   // Also: It has to be the spline dEdx, since one wants to get the sigma for the assumption(!)
765   // of such a particle, which by assumption then has this dEdx value
766     
767   Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE);
768   
769   if (dEdxExpected < 1.)
770     return 999;
771   
772   // For ESD tracks, the local tanTheta could be used (esdTrack->GetInnerParam()->GetTgl()).
773   // However, this value is not available for AODs and, thus, not or AliVTrack.
774   // Fortunately, the following formula allows to approximate the local tanTheta with the 
775   // global theta angle -> This is for by far most of the tracks the same, but gives at
776   // maybe the percent level differences within +- 0.2 in tanTheta -> Which is still ok.
777   Double_t tanTheta = TMath::Tan(-track->Theta() + TMath::Pi() / 2.0);
778   Int_t binX = fhEtaSigmaPar1->GetXaxis()->FindBin(tanTheta);
779   Int_t binY = fhEtaSigmaPar1->GetYaxis()->FindBin(1. / dEdxExpected);
780     
781   if (binX == 0) 
782     binX = 1;
783   if (binX > fhEtaSigmaPar1->GetXaxis()->GetNbins())
784     binX = fhEtaSigmaPar1->GetXaxis()->GetNbins();
785     
786   if (binY == 0)
787     binY = 1;
788   if (binY > fhEtaSigmaPar1->GetYaxis()->GetNbins())
789     binY = fhEtaSigmaPar1->GetYaxis()->GetNbins();
790     
791   return fhEtaSigmaPar1->GetBinContent(binX, binY);
792 }
793
794
795 //_________________________________________________________________________
796 Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
797 {
798   //
799   // Get eta correction for the given track.
800   //
801   
802   if (!fhEtaSigmaPar1)
803     return 999;
804   
805   Double_t dEdx = -1;
806   Int_t nPoints = -1;
807   ETPCgainScenario gainScenario = kGainScenarioInvalid;
808   TSpline3* responseFunction = 0x0;
809   
810   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
811     return 999; 
812   
813   return GetSigmaPar1(track, species, dEdx, responseFunction);
814 }
815
816
817 //_________________________________________________________________________
818 Bool_t AliTPCPIDResponse::SetEtaCorrMap(TH2D* hMap)
819 {
820   //
821   // Load map for TPC eta correction (a copy is stored and will be deleted automatically).
822   // If hMap is 0x0,the eta correction will be disabled and kFALSE is returned. 
823   // If the map can be set, kTRUE is returned.
824   //
825   
826   delete fhEtaCorr;
827   
828   if (!hMap) {
829     fhEtaCorr = 0x0;
830     
831     return kFALSE;
832   }
833   
834   fhEtaCorr = (TH2D*)(hMap->Clone());
835   fhEtaCorr->SetDirectory(0);
836       
837   return kTRUE;
838 }
839
840 //_________________________________________________________________________
841 Bool_t AliTPCPIDResponse::SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0)
842 {
843   //
844   // Load map for TPC sigma map (a copy is stored and will be deleted automatically):
845   // Parameter 1 is stored as a 2D map (1/dEdx vs. tanTheta_local) and parameter 0 is
846   // a constant. If hSigmaPar1Map is 0x0, the old sigma parametrisation will be used
847   // (and sigmaPar0 is ignored!) and kFALSE is returned. 
848   // If the map can be set, sigmaPar0 is also set and kTRUE will be returned.
849   //
850   
851   delete fhEtaSigmaPar1;
852   
853   if (!hSigmaPar1Map) {
854     fhEtaSigmaPar1 = 0x0;
855     fSigmaPar0 = 0.0;
856     
857     return kFALSE;
858   }
859   
860   fhEtaSigmaPar1 = (TH2D*)(hSigmaPar1Map->Clone());
861   fhEtaSigmaPar1->SetDirectory(0);
862   fSigmaPar0 = sigmaPar0;
863   
864   return kTRUE;
865 }
866
867
868 //_________________________________________________________________________
869 Bool_t AliTPCPIDResponse::sectorNumbersInOut(Double_t* trackPositionInner,
870                                              Double_t* trackPositionOuter,
871                                              Float_t& inphi,
872                                              Float_t& outphi,
873                                              Int_t& in,
874                                              Int_t& out ) const
875 {
876   //calculate the sector numbers (equivalent to IROC chamber numbers) a track crosses
877   //for OROC chamber numbers add 36
878   //returned angles are between (0,2pi)
879
880   inphi = TMath::ATan2(trackPositionInner[1],trackPositionInner[0]);
881   outphi = TMath::ATan2(trackPositionOuter[1], trackPositionOuter[0]);
882
883   if (inphi<0) {inphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi
884   if (outphi<0) {outphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi
885
886   in = sectorNumber(inphi);
887   out = sectorNumber(outphi);
888
889   //for the C side (positive z) offset by 18
890   if (trackPositionInner[2]>0.0)
891   {
892     in+=18;
893     out+=18;
894   }
895   return kTRUE;
896 }
897
898
899 //_____________________________________________________________________________
900 Int_t AliTPCPIDResponse::sectorNumber(Double_t phi) const
901 {
902   //calculate sector number
903   const Float_t width=TMath::TwoPi()/18.;
904   return TMath::Floor(phi/width);
905 }
906
907
908 //_____________________________________________________________________________
909 void AliTPCPIDResponse::Print(Option_t* /*option*/) const
910 {
911   //Print info
912   fResponseFunctions.Print();
913 }
914
915
916 //_____________________________________________________________________________
917 AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack* track, Int_t layer) const
918 {
919   //status of the track: if it crosses any bad chambers return kChamberOff;
920   //IROC:layer=1, OROC:layer=2
921   if (layer<1 || layer>2) layer=1;
922   Int_t in=0;
923   Int_t out=0;
924   Float_t inphi=0.;
925   Float_t outphi=0.;
926   Float_t innerRadius = (layer==1)?83.0:133.7;
927   Float_t outerRadius = (layer==1)?133.5:247.7;
928
929   /////////////////////////////////////////////////////////////////////////////
930   //find out where track enters and leaves the layer.
931   //
932   Double_t trackPositionInner[3];
933   Double_t trackPositionOuter[3];
934  
935   //if there is no inner param this could mean we're using the AOD track,
936   //we still can extrapolate from the vertex - so use those params.
937   const AliExternalTrackParam* ip = track->GetInnerParam();
938   if (ip) track=ip;
939
940   Bool_t trackAtInner = track->GetXYZAt(innerRadius, fMagField, trackPositionInner);
941   Bool_t trackAtOuter = track->GetXYZAt(outerRadius, fMagField, trackPositionOuter);
942
943   if (!trackAtInner)
944   {
945     //if we dont even enter inner radius we do nothing and return invalid
946     inphi=0.0;
947     outphi=0.0;
948     in=0;
949     out=0;
950     return kChamberInvalid;
951   }
952
953   if (!trackAtOuter)
954   {
955     //if we don't reach the outer radius check that the apex is indeed within the outer radius and use apex position
956     Bool_t haveApex = TrackApex(track, fMagField, trackPositionOuter);
957     Float_t apexRadius = TMath::Sqrt(trackPositionOuter[0]*trackPositionOuter[0]+trackPositionOuter[1]*trackPositionOuter[1]);
958     if ( haveApex && apexRadius<=outerRadius && apexRadius>innerRadius)
959     {
960       //printf("pt: %.2f, apexRadius: %.2f(%s), x: %.2f, y: %.2f\n",track->Pt(),apexRadius,(haveApex)?"OK":"BAD",trackPositionOuter[0],trackPositionOuter[1]);
961     }
962     else
963     {
964       inphi=0.0;
965       outphi=0.0;
966       in=0;
967       out=0;
968       return kChamberInvalid;
969     }
970   }
971
972
973   if (!sectorNumbersInOut(trackPositionInner,
974                           trackPositionOuter,
975                           inphi,
976                           outphi,
977                           in,
978                           out)) return kChamberInvalid;
979
980   /////////////////////////////////////////////////////////////////////////////
981   //now we have the location of the track we can check
982   //if it is in a good/bad chamber
983   //
984   Bool_t sideA = kTRUE;
985  
986   if (((in/18)%2==1) && ((out/18)%2==1)) sideA=kFALSE;
987
988   in=in%18;
989   out=out%18;
990
991   if (TMath::Abs(in-out)>9)
992   {
993     if (TMath::Max(in,out)==out)
994     {
995       Int_t tmp=out;
996       out = in;
997       in = tmp;
998       Float_t tmpphi=outphi;
999       outphi=inphi;
1000       inphi=tmpphi;
1001     }
1002     in-=18;
1003     inphi-=TMath::TwoPi();
1004   }
1005   else
1006   {
1007     if (TMath::Max(in,out)==in)
1008     {
1009       Int_t tmp=out;
1010       out=in;
1011       in=tmp;
1012       Float_t tmpphi=outphi;
1013       outphi=inphi;
1014       inphi=tmpphi;
1015     }
1016   }
1017
1018   Float_t trackLengthInBad = 0.;
1019   Float_t trackLengthInLowGain = 0.;
1020   Float_t trackLengthTotal = TMath::Abs(outphi-inphi);
1021   Float_t lengthFractionInBadSectors = 0.;
1022
1023   const Float_t sectorWidth = TMath::TwoPi()/18.;  
1024  
1025   for (Int_t i=in; i<=out; i++)
1026   {
1027     int j=i;
1028     if (i<0) j+=18;    //correct for the negative values
1029     if (!sideA) j+=18; //move to the correct side
1030    
1031     Float_t deltaPhi = 0.;
1032     Float_t phiEdge=sectorWidth*i;
1033     if (inphi>phiEdge) {deltaPhi=phiEdge+sectorWidth-inphi;}
1034     else if ((outphi>=phiEdge) && (outphi<(phiEdge+sectorWidth))) {deltaPhi=outphi-phiEdge;}
1035     else {deltaPhi=sectorWidth;}
1036    
1037     Float_t v = fVoltageMap[(layer==1)?(j):(j+36)];
1038     if (v<=fBadIROCthreshhold)
1039     {
1040       trackLengthInBad+=deltaPhi;
1041       lengthFractionInBadSectors=1.;
1042     }
1043     if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold)
1044     {
1045       trackLengthInLowGain+=deltaPhi;
1046       lengthFractionInBadSectors=1.;
1047     }
1048   }
1049
1050   //for now low gain and bad (off) chambers are treated equally
1051   if (trackLengthTotal>0)
1052     lengthFractionInBadSectors = (trackLengthInLowGain+trackLengthInBad)/trackLengthTotal;
1053
1054   //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);
1055  
1056   if (lengthFractionInBadSectors>fMaxBadLengthFraction)
1057   {
1058     //printf("%%%%%%%% %s kChamberLowGain\n",(layer==1)?"IROC":"OROC");
1059     return kChamberLowGain;
1060   }
1061  
1062   return kChamberHighGain;
1063 }
1064
1065
1066 //_____________________________________________________________________________
1067 Float_t AliTPCPIDResponse::MaxClusterRadius(const AliVTrack* track) const
1068 {
1069   //return the radius of the outermost padrow containing a cluster in TPC
1070   //for the track
1071   const TBits* clusterMap=track->GetTPCClusterMapPtr();
1072   if (!clusterMap) return 0.;
1073
1074   //from AliTPCParam, radius of first IROC row
1075   const Float_t rfirstIROC = 8.52249984741210938e+01;
1076   const Float_t padrowHeightIROC = 0.75;
1077   const Float_t rfirstOROC0 = 1.35100006103515625e+02;
1078   const Float_t padrowHeightOROC0 = 1.0;
1079   const Float_t rfirstOROC1 = 1.99350006103515625e+02;
1080   const Float_t padrowHeightOROC1 = 1.5;
1081
1082   Int_t maxPadRow=160;
1083   while ((--maxPadRow)>0 && !clusterMap->TestBitNumber(maxPadRow)){}
1084   if (maxPadRow>126) return rfirstOROC1+(maxPadRow-126-1)*padrowHeightOROC1;
1085   if (maxPadRow>62) return rfirstOROC0+(maxPadRow-62-1)*padrowHeightOROC0;
1086   if (maxPadRow>0) return rfirstIROC+(maxPadRow-1)*padrowHeightIROC;
1087   return 0.0;
1088 }
1089
1090
1091 //_____________________________________________________________________________
1092 Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const
1093 {
1094   //calculate the coordinates of the apex of the track
1095   Double_t x[3];
1096   track->GetXYZ(x);
1097   Double_t p[3];
1098   track->GetPxPyPz(p);
1099   Double_t r = 1./track->OneOverPt()/0.0299792458/magField; //signed - will determine the direction of b
1100   //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);
1101   //find orthogonal vector (Gram-Schmidt)
1102   Double_t alpha = (p[0]*x[0] + p[1]*x[1])/(p[0]*p[0] + p[1]*p[1]);
1103   Double_t b[2];
1104   b[0]=x[0]-alpha*p[0];
1105   b[1]=x[1]-alpha*p[1];
1106  
1107   Double_t norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
1108   if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
1109   b[0]/=norm;
1110   b[1]/=norm;
1111   b[0]*=r;
1112   b[1]*=r;
1113   b[0]+=x[0];
1114   b[1]+=x[1];
1115   //printf("center: x:%.2f, y:%.2f\n",b[0],b[1]);
1116  
1117   norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
1118   if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
1119  
1120   position[0]=b[0]+b[0]*TMath::Abs(r)/norm;
1121   position[1]=b[1]+b[1]*TMath::Abs(r)/norm;
1122   position[2]=0.;
1123   return kTRUE;
1124 }