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