Line data Source code
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 176 : 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 : TH1D *AliTOFPIDResponse::fHTOFtailResponse=NULL; // histogram to generate a TOF tail
39 :
40 : //_________________________________________________________________________
41 14 : AliTOFPIDResponse::AliTOFPIDResponse():
42 14 : fSigma(0),
43 14 : fPmax(0), // zero at 0.5 GeV/c for pp
44 14 : fTime0(0)
45 70 : {
46 14 : AliLog::SetClassDebugLevel("AliTOFPIDResponse",0);
47 14 : fPar[0] = 0.008;
48 14 : fPar[1] = 0.008;
49 14 : fPar[2] = 0.002;
50 14 : fPar[3] = 40.0;
51 :
52 14 : if(!fTOFtailResponse){
53 24 : 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);
54 8 : fTOFtailResponse->SetParameter(0,1);
55 8 : fTOFtailResponse->SetParameter(1,-26);
56 8 : fTOFtailResponse->SetParameter(2,1);
57 8 : fTOFtailResponse->SetParameter(3,0.89);
58 8 : fTOFtailResponse->SetNpx(10000);
59 : }
60 :
61 14 : LoadTOFtailHisto();
62 :
63 :
64 : // Reset T0 info
65 14 : ResetT0info();
66 14 : SetMomBoundary();
67 28 : }
68 : //_________________________________________________________________________
69 0 : AliTOFPIDResponse::AliTOFPIDResponse(Double_t *param):
70 0 : fSigma(param[0]),
71 0 : fPmax(0), // zero at 0.5 GeV/c for pp
72 0 : fTime0(0)
73 0 : {
74 : //
75 : // The main constructor
76 : //
77 : //
78 :
79 : //fPmax=TMath::Exp(-0.5*3*3)/fSigma; // ~3 sigma at 0.5 GeV/c for PbPb
80 :
81 0 : fPar[0] = 0.008;
82 0 : fPar[1] = 0.008;
83 0 : fPar[2] = 0.002;
84 0 : fPar[3] = 40.0;
85 :
86 0 : if(!fTOFtailResponse){
87 0 : 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);
88 0 : fTOFtailResponse->SetParameter(0,1);
89 0 : fTOFtailResponse->SetParameter(1,-26);
90 0 : fTOFtailResponse->SetParameter(2,1);
91 0 : fTOFtailResponse->SetParameter(3,0.89);
92 0 : fTOFtailResponse->SetNpx(10000);
93 : }
94 :
95 0 : LoadTOFtailHisto();
96 :
97 : // Reset T0 info
98 0 : ResetT0info();
99 0 : SetMomBoundary();
100 0 : }
101 : //_________________________________________________________________________
102 : void AliTOFPIDResponse::SetTOFtail(Float_t tail){
103 0 : LoadTOFtailHisto();
104 :
105 0 : if(!fTOFtailResponse){
106 0 : 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);
107 0 : fTOFtailResponse->SetParameter(0,1);
108 0 : fTOFtailResponse->SetParameter(1,-26);
109 0 : fTOFtailResponse->SetParameter(2,1);
110 0 : fTOFtailResponse->SetParameter(3,tail);
111 0 : fTOFtailResponse->SetNpx(10000);
112 0 : }
113 : else{
114 0 : fTOFtailResponse->SetParameter(3,tail);
115 :
116 0 : if(fHTOFtailResponse){ // adjust the TOF tail histo
117 0 : fHTOFtailResponse->Reset();
118 0 : for(Int_t i=1;i<=200;i++){
119 0 : Float_t x = fHTOFtailResponse->GetBinCenter(i);
120 0 : Float_t wx = fHTOFtailResponse->GetBinWidth(i)*0.5;
121 0 : fHTOFtailResponse->SetBinContent(i,fTOFtailResponse->Integral(x-wx,x+wx));
122 : }
123 0 : }
124 : }
125 0 : }
126 : void AliTOFPIDResponse::SetTOFtailAllPara(Float_t mean,Float_t tail){
127 0 : LoadTOFtailHisto();
128 :
129 0 : if(!fTOFtailResponse){
130 0 : 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);
131 0 : fTOFtailResponse->SetParameter(0,1);
132 0 : fTOFtailResponse->SetParameter(1,mean);
133 0 : fTOFtailResponse->SetParameter(2,1);
134 0 : fTOFtailResponse->SetParameter(3,tail);
135 0 : fTOFtailResponse->SetNpx(10000);
136 0 : }
137 : else{
138 0 : fTOFtailResponse->SetParameter(1,mean);
139 0 : fTOFtailResponse->SetParameter(3,tail);
140 :
141 :
142 0 : if(fHTOFtailResponse){ // adjust the TOF tail histo
143 0 : fHTOFtailResponse->Reset();
144 0 : for(Int_t i=1;i<=200;i++){
145 0 : Float_t x = fHTOFtailResponse->GetBinCenter(i);
146 0 : Float_t wx = fHTOFtailResponse->GetBinWidth(i)*0.5;
147 0 : fHTOFtailResponse->SetBinContent(i,fTOFtailResponse->Integral(x-wx,x+wx));
148 : }
149 0 : }
150 : }
151 0 : }
152 :
153 : //_________________________________________________________________________
154 : Double_t
155 : AliTOFPIDResponse::GetMismatchProbability(Double_t time,Double_t eta) const {
156 0 : if(!fHmismTOF){
157 0 : TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
158 0 : if(fmism) fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
159 0 : if(!fHmismTOF){
160 0 : printf("I cannot retrive TOF mismatch histos... skipped!");
161 0 : return 1E-4;
162 : }
163 0 : fHmismTOF->Scale(TMath::Sqrt(2*TMath::Pi())/(fHmismTOF->Integral(1,fHmismTOF->GetNbinsX()) * fHmismTOF->GetBinWidth(1)));
164 :
165 0 : TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
166 0 : if(fchDist) fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist");
167 0 : if(!fHchannelTOFdistr){
168 0 : printf("I cannot retrive TOF channel distance distribution... skipped!");
169 0 : return 1E-4;
170 : }
171 0 : }
172 :
173 0 : Float_t etaAbs = TMath::Abs(eta);
174 0 : Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
175 0 : if(channel < 1 || etaAbs > 1) channel = 1;
176 0 : Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
177 :
178 0 : Double_t mismWeight = fHmismTOF->Interpolate(time - distIP*3.35655419905265973e+01);
179 :
180 : return mismWeight;
181 0 : }
182 : //_________________________________________________________________________
183 : Double_t AliTOFPIDResponse::GetExpectedSigma(Float_t mom, Float_t time, Float_t mass) const {
184 : //
185 : // Return the expected sigma of the PID signal for the specified
186 : // particle mass/Z.
187 : // If the operation is not possible, return a negative value.
188 : //
189 :
190 20880990 : Double_t dpp=fPar[0] + fPar[1]*mom + fPar[2]*mass/mom; //mean relative pt resolution;
191 :
192 :
193 10440495 : Double_t sigma = dpp*time/(1.+ mom*mom/(mass*mass));
194 :
195 10440495 : Int_t index = GetMomBin(mom);
196 :
197 10440495 : Double_t t0res = fT0resolution[index];
198 :
199 10440495 : return TMath::Sqrt(sigma*sigma + fPar[3]*fPar[3]/mom/mom + fSigma*fSigma + t0res*t0res);
200 :
201 : }
202 : //_________________________________________________________________________
203 : Double_t AliTOFPIDResponse::GetExpectedSigma(Float_t mom, Float_t time, AliPID::EParticleType type) const {
204 : //
205 : // Return the expected sigma of the PID signal for the specified
206 : // particle type.
207 : // If the operation is not possible, return a negative value.
208 : //
209 :
210 0 : Double_t mass = AliPID::ParticleMassZ(type);
211 0 : Double_t dpp=fPar[0] + fPar[1]*mom + fPar[2]*mass/mom; //mean relative pt resolution;
212 :
213 :
214 0 : Double_t sigma = dpp*time/(1.+ mom*mom/(mass*mass));
215 :
216 0 : Int_t index = GetMomBin(mom);
217 :
218 0 : Double_t t0res = fT0resolution[index];
219 :
220 0 : return TMath::Sqrt(sigma*sigma + fPar[3]*fPar[3]/mom/mom + fSigma*fSigma + t0res*t0res);
221 :
222 : }
223 : //_________________________________________________________________________
224 : Double_t AliTOFPIDResponse::GetExpectedSignal(const AliVTrack* track,AliPID::EParticleType type) const {
225 : //
226 : // Return the expected signal of the PID signal for the particle type
227 : // If the operation is not possible, return a negative value.
228 : //
229 0 : Double_t expt[AliPID::kSPECIESC];
230 0 : track->GetIntegratedTimes(expt,AliPID::kSPECIESC);
231 0 : if (type<=AliPID::kProton) return expt[type];
232 : else {
233 0 : if (expt[type]<1.E-1) {
234 0 : Double_t p = track->P();
235 0 : Double_t massZ = AliPID::ParticleMassZ(type);
236 0 : return expt[0]/p*massZ*TMath::Sqrt(1.+p*p/massZ/massZ);
237 0 : } else return expt[type];
238 : }
239 0 : }
240 : //_________________________________________________________________________
241 : Int_t AliTOFPIDResponse::GetMomBin(Float_t p) const{
242 : //
243 : // Returns the momentum bin index
244 : //
245 :
246 : Int_t i=0;
247 178067727 : while(p > fPCutMin[i] && i < fNmomBins) i++;
248 20880990 : if(i > 0) i--;
249 :
250 10440495 : return i;
251 : }
252 : //_________________________________________________________________________
253 : void AliTOFPIDResponse::SetMomBoundary(){
254 : //
255 : // Set boundaries for momentum bins
256 : //
257 :
258 28 : fPCutMin[0] = 0.3;
259 14 : fPCutMin[1] = 0.5;
260 14 : fPCutMin[2] = 0.6;
261 14 : fPCutMin[3] = 0.7;
262 14 : fPCutMin[4] = 0.8;
263 14 : fPCutMin[5] = 0.9;
264 14 : fPCutMin[6] = 1;
265 14 : fPCutMin[7] = 1.2;
266 14 : fPCutMin[8] = 1.5;
267 14 : fPCutMin[9] = 2;
268 14 : fPCutMin[10] = 3;
269 14 : }
270 : //_________________________________________________________________________
271 : Float_t AliTOFPIDResponse::GetStartTime(Float_t mom) const {
272 : //
273 : // Returns event_time value as estimated by TOF combinatorial algorithm
274 : //
275 :
276 0 : Int_t ibin = GetMomBin(mom);
277 0 : return GetT0bin(ibin);
278 :
279 : }
280 : //_________________________________________________________________________
281 : Float_t AliTOFPIDResponse::GetStartTimeRes(Float_t mom) const {
282 : //
283 : // Returns event_time resolution as estimated by TOF combinatorial algorithm
284 : //
285 :
286 0 : Int_t ibin = GetMomBin(mom);
287 0 : return GetT0binRes(ibin);
288 :
289 : }
290 : //_________________________________________________________________________
291 : Int_t AliTOFPIDResponse::GetStartTimeMask(Float_t mom) const {
292 : //
293 : // Returns event_time mask
294 : //
295 :
296 0 : Int_t ibin = GetMomBin(mom);
297 0 : return GetT0binMask(ibin);
298 :
299 : }
300 : //_________________________________________________________________________
301 : 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)
302 : {
303 :
304 : // To add mismatch
305 0 : Float_t mismAdd = addmism*0.01;
306 0 : if(pt>1.0) mismAdd /= pt;
307 :
308 0 : if(mismAdd > 0.01){ // apply additional mismatch
309 0 : if(gRandom->Rndm() < mismAdd){
310 0 : return GetMismatchRandomValue(eta)-time;
311 : }
312 : }
313 :
314 0 : if(fHTOFtailResponse)
315 0 : return fHTOFtailResponse->GetRandom();
316 : else
317 0 : return 0.0;
318 0 : }
319 : //_________________________________________________________________________
320 : Double_t AliTOFPIDResponse::GetMismatchRandomValue(Float_t eta) // generate a random value for mismatched tracks (for MC analyses)
321 : {
322 0 : if(!fHmismTOF){
323 0 : TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
324 0 : if(fmism) fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
325 0 : if(!fHmismTOF){
326 0 : printf("I cannot retrive TOF mismatch histos... skipped!");
327 0 : return -10000.;
328 : }
329 0 : fHmismTOF->Scale(TMath::Sqrt(2*TMath::Pi())/(fHmismTOF->Integral(1,fHmismTOF->GetNbinsX()) * fHmismTOF->GetBinWidth(1)));
330 :
331 0 : TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
332 0 : if(fchDist) fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist");
333 0 : if(!fHchannelTOFdistr){
334 0 : printf("I cannot retrive TOF channel distance distribution... skipped!");
335 0 : return -10000.;
336 : }
337 0 : }
338 :
339 0 : Float_t etaAbs = TMath::Abs(eta);
340 0 : Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
341 0 : if(channel < 1 || etaAbs > 1) channel = 1;
342 0 : Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
343 :
344 0 : return fHmismTOF->GetRandom() + distIP*3.35655419905265973e+01;
345 0 : }
346 : //_________________________________________________________________________
347 : Int_t AliTOFPIDResponse::GetTOFchannel(AliVParticle *trk) const{
348 0 : Float_t etaAbs = TMath::Abs(trk->Eta());
349 0 : Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
350 0 : if(channel < 1 || etaAbs > 1) channel = 1;
351 :
352 0 : return channel;
353 : }
354 : //_________________________________________________________________________
355 : Int_t AliTOFPIDResponse::LoadTOFtailHisto(){
356 28 : if(! fHTOFtailResponse){
357 8 : TFile *fAddTail = new TFile("$ALICE_ROOT/TOF/data/addTOFtail.root");
358 16 : if(fAddTail) fHTOFtailResponse = (TH1D *) fAddTail->Get("hTOFTail");
359 16 : if(! fHTOFtailResponse){
360 0 : AliError("Cannot retrive TOF tail histogram from file $ALICE_ROOT/TOF/data/addTOFtail.root ... skipped!");
361 0 : delete fAddTail;
362 0 : return 2;
363 : }
364 : else{
365 16 : AliInfo("Loaded TOF tail histogram from file $ALICE_ROOT/TOF/data/addTOFtail.root");
366 8 : fHTOFtailResponse->SetDirectory(0x0);
367 : }
368 16 : delete fAddTail;
369 8 : return 0;
370 : }
371 :
372 6 : return 1;
373 14 : }
|
