STEER/STEERBase/AliTPCPIDResponse.cxx

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