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