STEER/STEERBase/AliTOFPIDResponse.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 //                                                                 //
  18 //           Implementation of the TOF PID class                   //
  19 //      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch         //
  20 //                                                                 //
  21 //-----------------------------------------------------------------//
  22
  23 #include "TMath.h"
  24 #include "AliLog.h"
  25 #include "TF1.h"
  26 #include "TH1F.h"
  27 #include "TH1D.h"
  28 #include "TFile.h"
  29 #include "TRandom.h"
  30
  31 #include "AliTOFPIDResponse.h"
  32
  33 ClassImp(AliTOFPIDResponse)
  34
  35 TF1 *AliTOFPIDResponse::fTOFtailResponse = NULL; // function to generate a TOF tail
  36 TH1F *AliTOFPIDResponse::fHmismTOF = NULL; // TOF mismatch distribution
  37 TH1D *AliTOFPIDResponse::fHchannelTOFdistr=NULL;  // TOF channel distance distribution
  38
  39 //_________________________________________________________________________
  40 AliTOFPIDResponse::AliTOFPIDResponse():
  41   fSigma(0),
  42   fPmax(0),         // zero at 0.5 GeV/c for pp
  43   fTime0(0)
  44 {
  45   fPar[0] = 0.008;
  46   fPar[1] = 0.008;
  47   fPar[2] = 0.002;
  48   fPar[3] = 40.0;
  49
  50   if(!fTOFtailResponse){
  51     fTOFtailResponse = new TF1("fTOFtail","[0]*TMath::Exp(-(x-[1])*(x-[1])/2/[2]/[2])* (x < [1]+[3]*[2]) + (x > [1]+[3]*[2])*[0]*TMath::Exp(-(x-[1]-[3]*[2]*0.5)*[3]/[2] * 0.0111)*0.018",-1000,1000);
  52     fTOFtailResponse->SetParameter(0,1);
  53     fTOFtailResponse->SetParameter(1,-26);
  54     fTOFtailResponse->SetParameter(2,1);
  55     fTOFtailResponse->SetParameter(3,0.89);
  56     fTOFtailResponse->SetNpx(10000);
  57   }
  58
  59
  60   // Reset T0 info
  61   ResetT0info();
  62   SetMomBoundary();
  63 }
  64 //_________________________________________________________________________
  65 AliTOFPIDResponse::AliTOFPIDResponse(Double_t *param):
  66   fSigma(param[0]),
  67   fPmax(0),          // zero at 0.5 GeV/c for pp
  68   fTime0(0)
  69 {
  70   //
  71   //  The main constructor
  72   //
  73   //
  74
  75   //fPmax=TMath::Exp(-0.5*3*3)/fSigma; // ~3 sigma at 0.5 GeV/c for PbPb
  76
  77   fPar[0] = 0.008;
  78   fPar[1] = 0.008;
  79   fPar[2] = 0.002;
  80   fPar[3] = 40.0;
  81
  82   if(!fTOFtailResponse){
  83     fTOFtailResponse = new TF1("fTOFtail","[0]*TMath::Exp(-(x-[1])*(x-[1])/2/[2]/[2])* (x < [1]+[3]*[2]) + (x > [1]+[3]*[2])*[0]*TMath::Exp(-(x-[1]-[3]*[2]*0.5)*[3]/[2] * 0.0111)*0.018",-1000,1000);
  84     fTOFtailResponse->SetParameter(0,1);
  85     fTOFtailResponse->SetParameter(1,-26);
  86     fTOFtailResponse->SetParameter(2,1);
  87     fTOFtailResponse->SetParameter(3,0.89);
  88     fTOFtailResponse->SetNpx(10000);
  89   }
  90
  91   // Reset T0 info
  92   ResetT0info();
  93   SetMomBoundary();
  94 }
  95 //_________________________________________________________________________
  96 void AliTOFPIDResponse::SetTOFtail(Float_t tail){
  97   if(!fTOFtailResponse){
  98     fTOFtailResponse = new TF1("fTOFtail","[0]*TMath::Exp(-(x-[1])*(x-[1])/2/[2]/[2])* (x < [1]+[3]*[2]) + (x > [1]+[3]*[2])*[0]*TMath::Exp(-(x-[1]-[3]*[2]*0.5)*[3]/[2] * 0.0111)*0.018",-1000,1000);
  99     fTOFtailResponse->SetParameter(0,1);
 100     fTOFtailResponse->SetParameter(1,-26);
 101     fTOFtailResponse->SetParameter(2,1);
 102     fTOFtailResponse->SetParameter(3,tail);
 103     fTOFtailResponse->SetNpx(10000);
 104   }
 105   else{
 106     fTOFtailResponse->SetParameter(3,tail);
 107   }
 108 }
 109 //_________________________________________________________________________
 110 Double_t
 111 AliTOFPIDResponse::GetMismatchProbability(Double_t time,Double_t eta) const {
 112   if(!fHmismTOF){
 113     TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
 114     if(fmism) fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
 115     if(!fHmismTOF){
 116       printf("I cannot retrive TOF mismatch histos... skipped!");
 117       return 1E-4;
 118     }
 119     fHmismTOF->Scale(TMath::Sqrt(2*TMath::Pi())/(fHmismTOF->Integral(1,fHmismTOF->GetNbinsX()) * fHmismTOF->GetBinWidth(1)));
 120
 121     TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
 122     if(fchDist) fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist");
 123     if(!fHchannelTOFdistr){
 124       printf("I cannot retrive TOF channel distance distribution... skipped!");
 125       return 1E-4;
 126     }
 127   }
 128
 129   Float_t etaAbs = TMath::Abs(eta);
 130   Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 131   if(channel < 1 || etaAbs > 1) channel = 1;
 132   Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 133
 134   Double_t mismWeight = fHmismTOF->Interpolate(time - distIP*3.35655419905265973e+01);
 135
 136   return mismWeight;
 137 }
 138 //_________________________________________________________________________
 139 Double_t AliTOFPIDResponse::GetExpectedSigma(Float_t mom, Float_t time, Float_t mass) const {
 140   //
 141   // Return the expected sigma of the PID signal for the specified
 142   // particle mass/Z.
 143   // If the operation is not possible, return a negative value.
 144   //
 145
 146   Double_t dpp=fPar[0] + fPar[1]*mom + fPar[2]*mass/mom;      //mean relative pt resolution;
 147
 148
 149   Double_t sigma = dpp*time/(1.+ mom*mom/(mass*mass));
 150
 151   Int_t index = GetMomBin(mom);
 152
 153   Double_t t0res = fT0resolution[index];
 154
 155   return TMath::Sqrt(sigma*sigma + fPar[3]*fPar[3]/mom/mom + fSigma*fSigma + t0res*t0res);
 156
 157 }
 158 //_________________________________________________________________________
 159 Double_t AliTOFPIDResponse::GetExpectedSigma(Float_t mom, Float_t time, AliPID::EParticleType  type) const {
 160   //
 161   // Return the expected sigma of the PID signal for the specified
 162   // particle type.
 163   // If the operation is not possible, return a negative value.
 164   //
 165
 166   Double_t mass = AliPID::ParticleMassZ(type);
 167   Double_t dpp=fPar[0] + fPar[1]*mom + fPar[2]*mass/mom;      //mean relative pt resolution;
 168
 169
 170   Double_t sigma = dpp*time/(1.+ mom*mom/(mass*mass));
 171
 172   Int_t index = GetMomBin(mom);
 173
 174   Double_t t0res = fT0resolution[index];
 175
 176   return TMath::Sqrt(sigma*sigma + fPar[3]*fPar[3]/mom/mom + fSigma*fSigma + t0res*t0res);
 177
 178 }
 179 //_________________________________________________________________________
 180 Double_t AliTOFPIDResponse::GetExpectedSignal(const AliVTrack* track,AliPID::EParticleType type) const {
 181   //
 182   // Return the expected signal of the PID signal for the particle type
 183   // If the operation is not possible, return a negative value.
 184   //
 185   Double_t expt[AliPID::kSPECIESC];
 186   track->GetIntegratedTimes(expt,AliPID::kSPECIESC);
 187   if (type<=AliPID::kProton) return expt[type];
 188   else {
 189     if (expt[type]<1.E-1) {
 190       Double_t p = track->P();
 191       Double_t massZ = AliPID::ParticleMassZ(type);
 192       return expt[0]/p*massZ*TMath::Sqrt(1.+p*p/massZ/massZ);
 193     } else return expt[type];
 194   }
 195 }
 196 //_________________________________________________________________________
 197 Int_t AliTOFPIDResponse::GetMomBin(Float_t p) const{
 198   //
 199   // Returns the momentum bin index
 200   //
 201
 202   Int_t i=0;
 203   while(p > fPCutMin[i] && i < fNmomBins) i++;
 204   if(i > 0) i--;
 205
 206   return i;
 207 }
 208 //_________________________________________________________________________
 209 void AliTOFPIDResponse::SetMomBoundary(){
 210   //
 211   // Set boundaries for momentum bins
 212   //
 213
 214   fPCutMin[0] = 0.3;
 215   fPCutMin[1] = 0.5;
 216   fPCutMin[2] = 0.6;
 217   fPCutMin[3] = 0.7;
 218   fPCutMin[4] = 0.8;
 219   fPCutMin[5] = 0.9;
 220   fPCutMin[6] = 1;
 221   fPCutMin[7] = 1.2;
 222   fPCutMin[8] = 1.5;
 223   fPCutMin[9] = 2;
 224   fPCutMin[10] = 3;
 225 }
 226 //_________________________________________________________________________
 227 Float_t AliTOFPIDResponse::GetStartTime(Float_t mom) const {
 228   //
 229   // Returns event_time value as estimated by TOF combinatorial algorithm
 230   //
 231
 232   Int_t ibin = GetMomBin(mom);
 233   return GetT0bin(ibin);
 234
 235 }
 236 //_________________________________________________________________________
 237 Float_t AliTOFPIDResponse::GetStartTimeRes(Float_t mom) const {
 238   //
 239   // Returns event_time resolution as estimated by TOF combinatorial algorithm
 240   //
 241
 242   Int_t ibin = GetMomBin(mom);
 243   return GetT0binRes(ibin);
 244
 245 }
 246 //_________________________________________________________________________
 247 Int_t AliTOFPIDResponse::GetStartTimeMask(Float_t mom) const {
 248   //
 249   // Returns event_time mask
 250   //
 251
 252   Int_t ibin = GetMomBin(mom);
 253   return GetT0binMask(ibin);
 254
 255 }
 256 //_________________________________________________________________________
 257 Double_t AliTOFPIDResponse::GetTailRandomValue(Float_t pt,Float_t eta,Float_t time,Float_t addmism) // generate a random value to add a tail to TOF time (for MC analyses)
 258 {
 259
 260   // To add mismatch
 261   Float_t mismAdd = addmism*0.01;
 262   if(pt>1.0) mismAdd /= pt;
 263
 264   if(mismAdd > 0.01){ // apply additional mismatch
 265     if(gRandom->Rndm() < mismAdd){
 266       return GetMismatchRandomValue(eta)-time;
 267     }
 268   }
 269
 270   if(fTOFtailResponse)
 271     return fTOFtailResponse->GetRandom();
 272   else
 273     return 0.0;
 274 }
 275 //_________________________________________________________________________
 276 Double_t AliTOFPIDResponse::GetMismatchRandomValue(Float_t eta) // generate a random value for mismatched tracks (for MC analyses)
 277 {
 278   if(!fHmismTOF){
 279     TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
 280     if(fmism) fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
 281     if(!fHmismTOF){
 282       printf("I cannot retrive TOF mismatch histos... skipped!");
 283       return -10000.;
 284     }
 285     fHmismTOF->Scale(TMath::Sqrt(2*TMath::Pi())/(fHmismTOF->Integral(1,fHmismTOF->GetNbinsX()) * fHmismTOF->GetBinWidth(1)));
 286
 287     TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
 288     if(fchDist) fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist");
 289     if(!fHchannelTOFdistr){
 290       printf("I cannot retrive TOF channel distance distribution... skipped!");
 291       return -10000.;
 292     }
 293   }
 294
 295   Float_t etaAbs = TMath::Abs(eta);
 296   Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 297   if(channel < 1 || etaAbs > 1) channel = 1;
 298   Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 299
 300   return fHmismTOF->GetRandom() + distIP*3.35655419905265973e+01;
 301 }
 302 //_________________________________________________________________________
 303 Int_t AliTOFPIDResponse::GetTOFchannel(AliVParticle *trk) const{
 304   Float_t etaAbs = TMath::Abs(trk->Eta());
 305   Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 306   if(channel < 1 || etaAbs > 1) channel = 1;
 307
 308   return channel;
 309 }