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[5];
 186   track->GetIntegratedTimes(expt);
 187   if (type<=AliPID::kProton) return expt[type];
 188   else {
 189     Double_t p = track->P();
 190     Double_t massZ = AliPID::ParticleMassZ(type);
 191     return expt[0]/p*massZ*TMath::Sqrt(1.+p*p/massZ/massZ);
 192   }
 193 }
 194 //_________________________________________________________________________
 195 Int_t AliTOFPIDResponse::GetMomBin(Float_t p) const{
 196   //
 197   // Returns the momentum bin index
 198   //
 199
 200   Int_t i=0;
 201   while(p > fPCutMin[i] && i < fNmomBins) i++;
 202   if(i > 0) i--;
 203
 204   return i;
 205 }
 206 //_________________________________________________________________________
 207 void AliTOFPIDResponse::SetMomBoundary(){
 208   //
 209   // Set boundaries for momentum bins
 210   //
 211
 212   fPCutMin[0] = 0.3;
 213   fPCutMin[1] = 0.5;
 214   fPCutMin[2] = 0.6;
 215   fPCutMin[3] = 0.7;
 216   fPCutMin[4] = 0.8;
 217   fPCutMin[5] = 0.9;
 218   fPCutMin[6] = 1;
 219   fPCutMin[7] = 1.2;
 220   fPCutMin[8] = 1.5;
 221   fPCutMin[9] = 2;
 222   fPCutMin[10] = 3;
 223 }
 224 //_________________________________________________________________________
 225 Float_t AliTOFPIDResponse::GetStartTime(Float_t mom) const {
 226   //
 227   // Returns event_time value as estimated by TOF combinatorial algorithm
 228   //
 229
 230   Int_t ibin = GetMomBin(mom);
 231   return GetT0bin(ibin);
 232
 233 }
 234 //_________________________________________________________________________
 235 Float_t AliTOFPIDResponse::GetStartTimeRes(Float_t mom) const {
 236   //
 237   // Returns event_time resolution as estimated by TOF combinatorial algorithm
 238   //
 239
 240   Int_t ibin = GetMomBin(mom);
 241   return GetT0binRes(ibin);
 242
 243 }
 244 //_________________________________________________________________________
 245 Int_t AliTOFPIDResponse::GetStartTimeMask(Float_t mom) const {
 246   //
 247   // Returns event_time mask
 248   //
 249
 250   Int_t ibin = GetMomBin(mom);
 251   return GetT0binMask(ibin);
 252
 253 }
 254 //_________________________________________________________________________
 255 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)
 256 {
 257
 258   // To add mismatch
 259   Float_t mismAdd = addmism*0.01;
 260   if(pt>1.0) mismAdd /= pt;
 261
 262   if(mismAdd > 0.01){ // apply additional mismatch
 263     if(gRandom->Rndm() < mismAdd){
 264       return GetMismatchRandomValue(eta)-time;
 265     }
 266   }
 267
 268   if(fTOFtailResponse)
 269     return fTOFtailResponse->GetRandom();
 270   else
 271     return 0.0;
 272 }
 273 //_________________________________________________________________________
 274 Double_t AliTOFPIDResponse::GetMismatchRandomValue(Float_t eta) // generate a random value for mismatched tracks (for MC analyses)
 275 {
 276   if(!fHmismTOF){
 277     TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
 278     if(fmism) fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
 279     if(!fHmismTOF){
 280       printf("I cannot retrive TOF mismatch histos... skipped!");
 281       return -10000.;
 282     }
 283     fHmismTOF->Scale(TMath::Sqrt(2*TMath::Pi())/(fHmismTOF->Integral(1,fHmismTOF->GetNbinsX()) * fHmismTOF->GetBinWidth(1)));
 284
 285     TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
 286     if(fchDist) fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist");
 287     if(!fHchannelTOFdistr){
 288       printf("I cannot retrive TOF channel distance distribution... skipped!");
 289       return -10000.;
 290     }
 291   }
 292
 293   Float_t etaAbs = TMath::Abs(eta);
 294   Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 295   if(channel < 1 || etaAbs > 1) channel = 1;
 296   Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
 297
 298   return fHmismTOF->GetRandom() + distIP*3.35655419905265973e+01;
 299 }
 300 //_________________________________________________________________________
 301 Int_t AliTOFPIDResponse::GetTOFchannel(AliVParticle *trk) const{
 302   Float_t etaAbs = TMath::Abs(trk->Eta());
 303   Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
 304   if(channel < 1 || etaAbs > 1) channel = 1;
 305
 306   return channel;
 307 }