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