STEER/STEERBase/AliTPCPIDResponse.cxx

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