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 //-----------------------------------------------------------------
  26
  27 #include <TGraph.h>
  28 #include <TObjArray.h>
  29 #include <TSpline.h>
  30 #include <TBits.h>
  31 #include <TMath.h>
  32
  33 #include <AliLog.h>
  34 #include "AliExternalTrackParam.h"
  35 #include "AliVTrack.h"
  36 #include "AliTPCPIDResponse.h"
  37 #include "AliTPCdEdxInfo.h"
  38
  39 ClassImp(AliTPCPIDResponse);
  40
  41 const char* AliTPCPIDResponse::fgkGainScenarioName[fgkNumberOfGainScenarios+1]=
  42 {
  43   "", //default - no name
  44   "1", //all high
  45   "2", //OROC only
  46   "unknown"
  47 };
  48
  49 //_________________________________________________________________________
  50 AliTPCPIDResponse::AliTPCPIDResponse():
  51   TNamed(),
  52   fMIP(50.),
  53   fRes0(),
  54   fResN2(),
  55   fKp1(0.0283086),
  56   fKp2(2.63394e+01),
  57   fKp3(5.04114e-11),
  58   fKp4(2.12543),
  59   fKp5(4.88663),
  60   fUseDatabase(kFALSE),
  61   fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
  62   fVoltageMap(72),
  63   fLowGainIROCthreshold(-40),
  64   fBadIROCthreshhold(-70),
  65   fLowGainOROCthreshold(-40),
  66   fBadOROCthreshhold(-40),
  67   fMaxBadLengthFraction(0.5),
  68   fCurrentResponseFunction(NULL),
  69   fCurrentdEdx(0.0),
  70   fCurrentNPoints(0),
  71   fCurrentGainScenario(kGainScenarioInvalid),
  72   fMagField(0.)
  73 {
  74   //
  75   //  The default constructor
  76   //
  77   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=0.07;fResN2[i]=0.0;}
  78 }
  79
  80 //_________________________________________________________________________
  81 AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param):
  82   TNamed(),
  83   fMIP(param[0]),
  84   fRes0(),
  85   fResN2(),
  86   fKp1(0.0283086),
  87   fKp2(2.63394e+01),
  88   fKp3(5.04114e-11),
  89   fKp4(2.12543),
  90   fKp5(4.88663),
  91   fUseDatabase(kFALSE),
  92   fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
  93   fVoltageMap(72),
  94   fLowGainIROCthreshold(-40),
  95   fBadIROCthreshhold(-70),
  96   fLowGainOROCthreshold(-40),
  97   fBadOROCthreshhold(-40),
  98   fMaxBadLengthFraction(0.5),
  99   fCurrentResponseFunction(NULL),
 100   fCurrentdEdx(0.0),
 101   fCurrentNPoints(0),
 102   fCurrentGainScenario(kGainScenarioInvalid),
 103   fMagField(0.)
 104 {
 105   //
 106   //  The main constructor
 107   //
 108   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=param[1];fResN2[i]=param[2];}
 109 }
 110
 111 //_________________________________________________________________________
 112 AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that):
 113   TNamed(that),
 114   fMIP(that.fMIP),
 115   fRes0(),
 116   fResN2(),
 117   fKp1(that.fKp1),
 118   fKp2(that.fKp2),
 119   fKp3(that.fKp3),
 120   fKp4(that.fKp4),
 121   fKp5(that.fKp5),
 122   fUseDatabase(that.fUseDatabase),
 123   fResponseFunctions(that.fResponseFunctions),
 124   fVoltageMap(that.fVoltageMap),
 125   fLowGainIROCthreshold(that.fLowGainIROCthreshold),
 126   fBadIROCthreshhold(that.fBadIROCthreshhold),
 127   fLowGainOROCthreshold(that.fLowGainOROCthreshold),
 128   fBadOROCthreshhold(that.fBadOROCthreshhold),
 129   fMaxBadLengthFraction(that.fMaxBadLengthFraction),
 130   fCurrentResponseFunction(NULL),
 131   fCurrentdEdx(0.0),
 132   fCurrentNPoints(0),
 133   fCurrentGainScenario(kGainScenarioInvalid),
 134   fMagField(that.fMagField)
 135 {
 136   //copy ctor
 137   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
 138 }
 139
 140 //_________________________________________________________________________
 141 AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that)
 142 {
 143   //assignment
 144   if (&that==this) return *this;
 145   TNamed::operator=(that);
 146   fMIP=that.fMIP;
 147   fKp1=that.fKp1;
 148   fKp2=that.fKp2;
 149   fKp3=that.fKp3;
 150   fKp4=that.fKp4;
 151   fKp5=that.fKp5;
 152   fUseDatabase=that.fUseDatabase;
 153   fResponseFunctions=that.fResponseFunctions;
 154   fVoltageMap=that.fVoltageMap;
 155   fLowGainIROCthreshold=that.fLowGainIROCthreshold;
 156   fBadIROCthreshhold=that.fBadIROCthreshhold;
 157   fLowGainOROCthreshold=that.fLowGainOROCthreshold;
 158   fBadOROCthreshhold=that.fBadOROCthreshhold;
 159   fMaxBadLengthFraction=that.fMaxBadLengthFraction;
 160   fMagField=that.fMagField;
 161   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
 162   return *this;
 163 }
 164
 165 //_________________________________________________________________________
 166 Double_t AliTPCPIDResponse::Bethe(Double_t betaGamma) const {
 167   //
 168   // This is the Bethe-Bloch function normalised to 1 at the minimum
 169   // WARNING
 170   // Simulated and reconstructed Bethe-Bloch differs
 171   //           Simulated  curve is the dNprim/dx
 172   //           Reconstructed is proportianal dNtot/dx
 173   // Temporary fix for production -  Simple linear correction function
 174   // Future    2 Bethe Bloch formulas needed
 175   //           1. for simulation
 176   //           2. for reconstructed PID
 177   //
 178
 179 //   const Float_t kmeanCorrection =0.1;
 180   Double_t bb=
 181     AliExternalTrackParam::BetheBlochAleph(betaGamma,fKp1,fKp2,fKp3,fKp4,fKp5);
 182   return bb*fMIP;
 183 }
 184
 185 //_________________________________________________________________________
 186 void AliTPCPIDResponse::SetBetheBlochParameters(Double_t kp1,
 187                              Double_t kp2,
 188                              Double_t kp3,
 189                              Double_t kp4,
 190                              Double_t kp5) {
 191   //
 192   // Set the parameters of the ALEPH Bethe-Bloch formula
 193   //
 194   fKp1=kp1;
 195   fKp2=kp2;
 196   fKp3=kp3;
 197   fKp4=kp4;
 198   fKp5=kp5;
 199 }
 200
 201 //_________________________________________________________________________
 202 void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2) {
 203   //
 204   // Set the relative resolution  sigma_rel = res0 * sqrt(1+resN2/npoint)
 205   //
 206   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=res0;fResN2[i]=resN2;}
 207 }
 208
 209 //_________________________________________________________________________
 210 Double_t AliTPCPIDResponse::GetExpectedSignal(const Float_t mom,
 211                                               AliPID::EParticleType n) const {
 212   //
 213   // Calculates the expected PID signal as the function of
 214   // the information stored in the track, for the specified particle type
 215   //
 216   // At the moment, these signals are just the results of calling the
 217   // Bethe-Bloch formula.
 218   // This can be improved. By taking into account the number of
 219   // assigned clusters and/or the track dip angle, for example.
 220   //
 221
 222   Double_t mass=AliPID::ParticleMassZ(n);
 223   if (!fUseDatabase) return Bethe(mom/mass);
 224   //
 225   const TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n);
 226
 227   if (!responseFunction) return Bethe(mom/mass);
 228   //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
 229   // !!! Splines for light nuclei need to be normalised to this factor !!!
 230   const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(n),2.3);
 231   return fMIP*responseFunction->Eval(mom/mass)*chargeFactor;
 232
 233 }
 234
 235 //_________________________________________________________________________
 236 Double_t AliTPCPIDResponse::GetExpectedSigma(const Float_t mom,
 237                                              const Int_t nPoints,
 238                                              AliPID::EParticleType n) const {
 239   //
 240   // Calculates the expected sigma of the PID signal as the function of
 241   // the information stored in the track, for the specified particle type
 242   //
 243
 244   if (nPoints != 0)
 245     return GetExpectedSignal(mom,n)*fRes0[0]*sqrt(1. + fResN2[0]/nPoints);
 246   else
 247     return GetExpectedSignal(mom,n)*fRes0[0];
 248 }
 249
 250 ////////////////////////////////////////////////////NEW//////////////////////////////
 251
 252 //_________________________________________________________________________
 253 void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario) {
 254   //
 255   // Set the relative resolution  sigma_rel = res0 * sqrt(1+resN2/npoint)
 256   //
 257   if ((Int_t)gainScenario>(Int_t)fgkNumberOfGainScenarios) return; //TODO: better handling!
 258   fRes0[gainScenario]=res0;
 259   fResN2[gainScenario]=resN2;
 260 }
 261
 262 //_________________________________________________________________________
 263 Double_t AliTPCPIDResponse::GetExpectedSignal(Double_t mom,
 264                                               AliPID::EParticleType species,
 265                                               const TSpline3* responseFunction) const
 266 {
 267   // Calculates the expected PID signal as the function of
 268   // the information stored in the track, for the specified particle type
 269   //
 270   // At the moment, these signals are just the results of calling the
 271   // Bethe-Bloch formula.
 272   // This can be improved. By taking into account the number of
 273   // assigned clusters and/or the track dip angle, for example.
 274   //
 275
 276
 277   Double_t mass=AliPID::ParticleMass(species);
 278
 279   if (!fUseDatabase) return Bethe(mom/mass);
 280
 281   if (!responseFunction) return Bethe(mom/mass);
 282   return fMIP*responseFunction->Eval(mom/mass);
 283 }
 284
 285 //_________________________________________________________________________
 286 Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
 287                                               AliPID::EParticleType species,
 288                                               ETPCdEdxSource dedxSource)
 289 {
 290   // Calculates the expected PID signal as the function of
 291   // the information stored in the track, for the specified particle type
 292   //
 293   // At the moment, these signals are just the results of calling the
 294   // Bethe-Bloch formula.
 295   // This can be improved. By taking into account the number of
 296   // assigned clusters and/or the track dip angle, for example.
 297   //
 298
 299   Double_t mom = track->GetTPCmomentum();
 300   Double_t mass=AliPID::ParticleMass(species);
 301
 302   if (!fUseDatabase) return Bethe(mom/mass);
 303
 304   const TSpline3* responseFunction = GetResponseFunction(track,species,dedxSource);
 305   if (!responseFunction) return Bethe(mom/mass);
 306   return fMIP*responseFunction->Eval(mom/mass);
 307
 308 }
 309
 310 //_________________________________________________________________________
 311 TSpline3* AliTPCPIDResponse::GetResponseFunction( AliPID::EParticleType type,
 312                                                   AliTPCPIDResponse::ETPCgainScenario gainScenario ) const
 313 {
 314   //get response function
 315   return dynamic_cast<TSpline3*>(fResponseFunctions.At(ResponseFunctionIndex(type,gainScenario)));
 316 }
 317
 318 //_________________________________________________________________________
 319 TSpline3* AliTPCPIDResponse::GetResponseFunction( const AliVTrack* track,
 320                                AliPID::EParticleType species,
 321                                ETPCdEdxSource dedxSource )
 322 {
 323   //the splines are stored in an array, different scenarios
 324
 325   if (ResponseFunctiondEdxN(track,
 326                             species,
 327                             dedxSource))
 328     return fCurrentResponseFunction;
 329
 330   return NULL;
 331 }
 332
 333 //_________________________________________________________________________
 334 void AliTPCPIDResponse::ResetSplines()
 335 {
 336   //reset the array with splines
 337   for (Int_t i=0;i<fResponseFunctions.GetEntriesFast();i++)
 338   {
 339     fResponseFunctions.AddAt(NULL,i);
 340   }
 341 }
 342 //_________________________________________________________________________
 343 Int_t AliTPCPIDResponse::ResponseFunctionIndex( AliPID::EParticleType species,
 344                                                 ETPCgainScenario gainScenario ) const
 345 {
 346   //get the index in fResponseFunctions given type and scenario
 347   return Int_t(species)+Int_t(gainScenario)*fgkNumberOfParticleSpecies;
 348 }
 349
 350 //_________________________________________________________________________
 351 void AliTPCPIDResponse::SetResponseFunction( TObject* o,
 352                                              AliPID::EParticleType species,
 353                                              ETPCgainScenario gainScenario )
 354 {
 355   fResponseFunctions.AddAtAndExpand(o,ResponseFunctionIndex(species,gainScenario));
 356 }
 357
 358 //_________________________________________________________________________
 359 Double_t AliTPCPIDResponse::GetExpectedSigma( Double_t mom,
 360                                               Int_t nPoints,
 361                                               AliPID::EParticleType species,
 362                                               ETPCgainScenario gainScenario,
 363                                               const TSpline3* responseFunction) const
 364 {
 365   // Calculates the expected sigma of the PID signal as the function of
 366   // the information stored in the track, for the specified particle type
 367   // and dedx scenario
 368   //
 369
 370   if (!responseFunction) return 999;
 371   if (nPoints != 0)
 372     return GetExpectedSignal(mom,species,responseFunction) *
 373            fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints);
 374   else
 375     return GetExpectedSignal(mom,species,responseFunction)*fRes0[gainScenario];
 376 }
 377
 378 //_________________________________________________________________________
 379 Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
 380                                              AliPID::EParticleType species,
 381                                              ETPCdEdxSource dedxSource)
 382 {
 383   // Calculates the expected sigma of the PID signal as the function of
 384   // the information stored in the track, for the specified particle type
 385   // and dedx scenario
 386   //
 387
 388   if (!ResponseFunctiondEdxN(track,
 389                              species,
 390                              dedxSource)) return 998; //TODO: better handling!
 391
 392   return GetExpectedSigma( track->GetTPCmomentum(),
 393                            fCurrentNPoints,
 394                            species,
 395                            fCurrentGainScenario,
 396                            fCurrentResponseFunction );
 397 }
 398
 399 //_________________________________________________________________________
 400 Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track,
 401                              AliPID::EParticleType species,
 402                              ETPCdEdxSource dedxSource)
 403 {
 404   //calculates the number of sigmas of the PID signal from the expected value
 405   //for a gicven particle species in the presence of multiple gain scenarios
 406   //inside the TPC
 407
 408   if (!ResponseFunctiondEdxN(track,
 409                              species,
 410                              dedxSource)) return 997; //TODO: better handling!
 411
 412   Double_t mom = track->GetTPCmomentum();
 413   Double_t bethe=GetExpectedSignal(mom,species,fCurrentResponseFunction);
 414   Double_t sigma=GetExpectedSigma( mom,
 415                                    fCurrentNPoints,
 416                                    species,
 417                                    fCurrentGainScenario,
 418                                    fCurrentResponseFunction );
 419   return (fCurrentdEdx-bethe)/sigma;
 420 }
 421
 422 //_________________________________________________________________________
 423 Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track,
 424                                                  AliPID::EParticleType species,
 425                                                  ETPCdEdxSource dedxSource )
 426 {
 427   // Calculates the right parameters for PID
 428   //   dEdx parametrization for the proper gain scenario, dEdx
 429   //   and NPoints used for dEdx
 430   // based on the track geometry (which chambers it crosses) for the specified particle type
 431   // and preferred source of dedx.
 432   // returns true on success
 433
 434   Double32_t signal[4]; //0: IROC, 1: OROC medium, 2:OROC long, 3: OROC all (def. truncation used)
 435   Char_t ncl[3];        //same
 436   Char_t nrows[3];      //same
 437   AliTPCdEdxInfo* dEdxInfo = track->GetTPCdEdxInfo();
 438
 439   if (!dEdxInfo && dedxSource!=kdEdxDefault)  //in one case its ok if we dont have the info
 440   {
 441     AliError("AliTPCdEdxInfo not available");
 442     InvalidateCurrentValues();
 443     return kFALSE;
 444   }
 445
 446   if (dEdxInfo) dEdxInfo->GetTPCSignalRegionInfo(signal,ncl,nrows);
 447
 448   //check if we cross a bad OROC in which case we reject
 449   EChamberStatus trackStatus = TrackStatus(track,2);
 450   if (trackStatus==kChamberOff || trackStatus==kChamberLowGain)
 451   {
 452     InvalidateCurrentValues();
 453     return kFALSE;
 454   }
 455
 456   switch (dedxSource)
 457   {
 458     case kdEdxDefault:
 459       {
 460         fCurrentdEdx = track->GetTPCsignal();
 461         fCurrentNPoints = track->GetTPCsignalN();
 462         fCurrentGainScenario = kDefault;
 463         break;
 464       }
 465     case kdEdxOROC:
 466       {
 467         fCurrentdEdx = signal[3];
 468         fCurrentNPoints = ncl[2]+ncl[1];
 469         fCurrentGainScenario = kOROChigh;
 470         break;
 471       }
 472     case kdEdxHybrid:
 473       {
 474         //if we cross a bad IROC we use OROC dedx, if we dont we use combined
 475         EChamberStatus status = TrackStatus(track,1);
 476         if (status!=kChamberHighGain)
 477         {
 478           fCurrentdEdx = signal[3];
 479           fCurrentNPoints = ncl[2]+ncl[1];
 480           fCurrentGainScenario = kOROChigh;
 481         }
 482         else
 483         {
 484           fCurrentdEdx = track->GetTPCsignal();
 485           fCurrentNPoints = track->GetTPCsignalN();
 486           fCurrentGainScenario = kALLhigh;
 487         }
 488         break;
 489       }
 490     default:
 491       {
 492          fCurrentdEdx = 0.;
 493          fCurrentNPoints = 0;
 494          fCurrentGainScenario = kGainScenarioInvalid;
 495          return kFALSE;
 496       }
 497   }
 498   TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,fCurrentGainScenario));
 499   fCurrentResponseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast?
 500   return kTRUE;
 501 }
 502
 503 //_________________________________________________________________________
 504 Bool_t AliTPCPIDResponse::sectorNumbersInOut(const AliVTrack* track,
 505                                              Double_t innerRadius,
 506                                              Double_t outerRadius,
 507                                              Float_t& inphi,
 508                                              Float_t& outphi,
 509                                              Int_t& in,
 510                                              Int_t& out ) const
 511 {
 512   //calculate the sector numbers (equivalent to IROC chamber numbers) a track crosses
 513   //for OROC chamber numbers add 36
 514   //returned angles are between (0,2pi)
 515
 516   Double_t trackPositionInner[3];
 517   Double_t trackPositionOuter[3];
 518
 519   Bool_t trackAtInner=kTRUE;
 520   Bool_t trackAtOuter=kTRUE;
 521   const AliExternalTrackParam* ip = track->GetInnerParam();
 522
 523   //if there is no inner param this could mean we're using the AOD track,
 524   //we still can extrapolate from the vertex - so use those params.
 525   if (ip) track=ip;
 526
 527   if (!track->GetXYZAt(innerRadius, fMagField, trackPositionInner))
 528     trackAtInner=kFALSE;
 529
 530   if (!track->GetXYZAt(outerRadius, fMagField, trackPositionOuter))
 531     trackAtOuter=kFALSE;
 532
 533   if (!trackAtInner)
 534   {
 535     //if we dont even enter inner radius we do nothing
 536     inphi=0.0;
 537     outphi=0.0;
 538     in=0;
 539     out=0;
 540     return kFALSE;
 541   }
 542
 543   if (!trackAtOuter)
 544   {
 545     //if we don't reach the outer radius check that the apex is indeed within the outer radius and use apex position
 546     Bool_t haveApex = TrackApex(track, fMagField, trackPositionOuter);
 547     Float_t apexRadius = TMath::Sqrt(trackPositionOuter[0]*trackPositionOuter[0]+trackPositionOuter[1]*trackPositionOuter[1]);
 548     if ( haveApex && apexRadius<=outerRadius && apexRadius>innerRadius)
 549     {
 550       //printf("pt: %.2f, apexRadius: %.2f(%s), x: %.2f, y: %.2f\n",track->Pt(),apexRadius,(haveApex)?"OK":"BAD",trackPositionOuter[0],trackPositionOuter[1]);
 551     }
 552     else
 553     {
 554       inphi=0.0;
 555       outphi=0.0;
 556       in=0;
 557       out=0;
 558       return kFALSE;
 559     }
 560   }
 561
 562   inphi = TMath::ATan2(trackPositionInner[1],trackPositionInner[0]);
 563   outphi = TMath::ATan2(trackPositionOuter[1], trackPositionOuter[0]);
 564
 565   if (inphi<0) {inphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi
 566   if (outphi<0) {outphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi
 567
 568   in = sectorNumber(inphi);
 569   out = sectorNumber(outphi);
 570
 571   //for the C side (positive z) offset by 18
 572   if (trackPositionInner[2]>0.0)
 573   {
 574     in+=18;
 575     out+=18;
 576   }
 577   return kTRUE;
 578 }
 579
 580 //_____________________________________________________________________________
 581 Int_t AliTPCPIDResponse::sectorNumber(Double_t phi) const
 582 {
 583   //calculate sector number
 584   const Float_t width=TMath::TwoPi()/18.;
 585   return TMath::Floor(phi/width);
 586 }
 587
 588 //_____________________________________________________________________________
 589 void AliTPCPIDResponse::Print(Option_t* /*option*/) const
 590 {
 591   //Print info
 592   fResponseFunctions.Print();
 593 }
 594
 595 //_____________________________________________________________________________
 596 AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack* track, Int_t layer) const
 597 {
 598   //status of the track: if it crosses any bad chambers return kChamberOff;
 599   //IROC:layer=1, OROC:layer=2
 600   if (layer<1 || layer>2) layer=1;
 601   Int_t in=0;
 602   Int_t out=0;
 603   Float_t inphi=0.;
 604   Float_t outphi=0.;
 605   Float_t innerRadius = (layer==1)?83.0:133.7;
 606   Float_t outerRadius = (layer==1)?133.5:247.7;
 607   if (!sectorNumbersInOut(track,
 608                           innerRadius,
 609                           outerRadius,
 610                           inphi,
 611                           outphi,
 612                           in,
 613                           out)) return kChamberOff;
 614
 615   Bool_t sideA = kTRUE;
 616
 617   if (((in/18)%2==1) && ((out/18)%2==1)) sideA=kFALSE;
 618
 619   in=in%18;
 620   out=out%18;
 621
 622   if (TMath::Abs(in-out)>9)
 623   {
 624     if (TMath::Max(in,out)==out)
 625     {
 626       Int_t tmp=out;
 627       out = in;
 628       in = tmp;
 629       Float_t tmpphi=outphi;
 630       outphi=inphi;
 631       inphi=tmpphi;
 632     }
 633     in-=18;
 634     inphi-=TMath::TwoPi();
 635   }
 636   else
 637   {
 638     if (TMath::Max(in,out)==in)
 639     {
 640       Int_t tmp=out;
 641       out=in;
 642       in=tmp;
 643       Float_t tmpphi=outphi;
 644       outphi=inphi;
 645       inphi=tmpphi;
 646     }
 647   }
 648
 649   Float_t trackLengthInBad = 0.;
 650   Float_t trackLengthInLowGain = 0.;
 651   Float_t trackLengthTotal = TMath::Abs(outphi-inphi);
 652
 653   const Float_t sectorWidth = TMath::TwoPi()/18.;
 654
 655   for (Int_t i=in; i<=out; i++)
 656   {
 657     int j=i;
 658     if (i<0) j+=18;    //correct for the negative values
 659     if (!sideA) j+=18; //move to the correct side
 660
 661     Float_t deltaPhi = 0.;
 662     Float_t phiEdge=sectorWidth*i;
 663     if (inphi>phiEdge) {deltaPhi=phiEdge+sectorWidth-inphi;}
 664     else if ((outphi>=phiEdge) && (outphi<(phiEdge+sectorWidth))) {deltaPhi=outphi-phiEdge;}
 665     else {deltaPhi=sectorWidth;}
 666
 667     Float_t v = fVoltageMap[(layer==1)?(j):(j+36)];
 668     if (v<=fBadIROCthreshhold) trackLengthInBad+=deltaPhi;
 669     if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold) trackLengthInLowGain+=deltaPhi;
 670   }
 671
 672   //for now low gain and bad (off) chambers are treated equally
 673   Float_t lengthFractionInBadSectors = (trackLengthInLowGain+trackLengthInBad)/trackLengthTotal;
 674
 675   //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);
 676
 677   if (lengthFractionInBadSectors>fMaxBadLengthFraction)
 678   {
 679     //printf("%%%%%%%% %s kChamberLowGain\n",(layer==1)?"IROC":"OROC");
 680     return kChamberLowGain;
 681   }
 682
 683   return kChamberHighGain;
 684 }
 685
 686 //_____________________________________________________________________________
 687 Float_t AliTPCPIDResponse::MaxClusterRadius(const AliVTrack* track) const
 688 {
 689   //return the radius of the outermost padrow containing a cluster in TPC
 690   //for the track
 691   const TBits* clusterMap=track->GetTPCClusterMapPtr();
 692   if (!clusterMap) return 0.;
 693
 694   //from AliTPCParam, radius of first IROC row
 695   const Float_t rfirstIROC = 8.52249984741210938e+01;
 696   const Float_t padrowHeightIROC = 0.75;
 697   const Float_t rfirstOROC0 = 1.35100006103515625e+02;
 698   const Float_t padrowHeightOROC0 = 1.0;
 699   const Float_t rfirstOROC1 = 1.99350006103515625e+02;
 700   const Float_t padrowHeightOROC1 = 1.5;
 701
 702   Int_t maxPadRow=160;
 703   while ((--maxPadRow)>0 && !clusterMap->TestBitNumber(maxPadRow)){}
 704   if (maxPadRow>126) return rfirstOROC1+(maxPadRow-126-1)*padrowHeightOROC1;
 705   if (maxPadRow>62) return rfirstOROC0+(maxPadRow-62-1)*padrowHeightOROC0;
 706   if (maxPadRow>0) return rfirstIROC+(maxPadRow-1)*padrowHeightIROC;
 707   return 0.0;
 708 }
 709
 710 //_____________________________________________________________________________
 711 void AliTPCPIDResponse::InvalidateCurrentValues()
 712 {
 713   //make the "current" stored values from the last processed track invalid
 714   fCurrentResponseFunction=NULL;
 715   fCurrentdEdx=0.;
 716   fCurrentNPoints=0;
 717   fCurrentGainScenario=kGainScenarioInvalid;
 718 }
 719
 720 //_____________________________________________________________________________
 721 Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const
 722 {
 723   //calculate the coordinates of the apex of the track
 724   Double_t x[3];
 725   track->GetXYZ(x);
 726   Double_t p[3];
 727   track->GetPxPyPz(p);
 728   Double_t r = 1./track->OneOverPt()/0.0299792458/magField; //signed - will determine the direction of b
 729   //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);
 730   //find orthogonal vector (Gram-Schmidt)
 731   Double_t alpha = (p[0]*x[0] + p[1]*x[1])/(p[0]*p[0] + p[1]*p[1]);
 732   Double_t b[2];
 733   b[0]=x[0]-alpha*p[0];
 734   b[1]=x[1]-alpha*p[1];
 735
 736   Double_t norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
 737   if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
 738   b[0]/=norm;
 739   b[1]/=norm;
 740   b[0]*=r;
 741   b[1]*=r;
 742   b[0]+=x[0];
 743   b[1]+=x[1];
 744   //printf("center: x:%.2f, y:%.2f\n",b[0],b[1]);
 745
 746   norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
 747   if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
 748
 749   position[0]=b[0]+b[0]*TMath::Abs(r)/norm;
 750   position[1]=b[1]+b[1]*TMath::Abs(r)/norm;
 751   position[2]=0.;
 752   return kTRUE;
 753 }
 754