]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TOF/AliTOFT0maker.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
bugfix: corrected calculation of slice and partition from track point Id
[u/mrichter/AliRoot.git] / TOF / AliTOFT0maker.cxx
... / ...
CommitLineData
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/* $Id: AliTOFT0maker.cxx,v 1.8 2010/01/19 16:32:20 noferini Exp $ */
16
17/////////////////////////////////////////////////////////////////////////////
18// //
19// This class contains the basic functions for the time zero //
20// evaluation with TOF detector informations. //
21// Use case in an analysis task: //
22// //
23// Create the object in the task constructor (fTOFmaker is a private var) //
24// AliESDpid *extPID=new AliESDpid(); //
25// fTOFmaker = new AliTOFT0maker(extPID); //
26// fTOFmaker->SetTimeResolution(100.0); // if you want set the TOF res //
27// 115 ps is the TOF default resolution value //
28// //
29// Use the RemakePID method in the task::Exec //
30// Double_t* calcolot0; //
31// calcolot0=fTOFmaker->RemakePID(fESD); //
32// //calcolot0[0] = calculated event time //
33// //calcolot0[1] = event time time resolution //
34// //calcolot0[2] = average event time for the current fill //
35// //calcolot0[3] = tracks at TOF //
36// //calcolot0[4] = calculated event time (only TOF) //
37// //calcolot0[5] = event time time resolution (only TOF) //
38// //calcolot0[6] = sigma t0 fill //
39// //calcolot0[7] = tracks at TOF really used in tht algorithm //
40// //
41// Let consider that: //
42// - the PIF is automatically recalculated with the event time subtrction //
43// //
44/////////////////////////////////////////////////////////////////////////////
45
46#include "AliTOFT0v1.h"
47#include "AliTOFT0maker.h"
48#include "AliPID.h"
49#include "AliLog.h"
50#include "AliESDpid.h"
51#include "AliESDEvent.h"
52#include "TFile.h"
53#include "TH1F.h"
54#include "AliTOFcalib.h"
55#include "AliTOFRunParams.h"
56#include "TRandom.h"
57#include "AliTOFHeader.h"
58
59ClassImp(AliTOFT0maker)
60
61//____________________________________________________________________________
62AliTOFT0maker::AliTOFT0maker():
63 TObject(),
64 fT0TOF(NULL),
65 fPIDesd(NULL),
66 fExternalPIDFlag(kFALSE),
67 fTOFcalib(NULL),
68 fNoTOFT0(0),
69 fNmomBins(0),
70 fTimeResolution(100),
71 fT0sigma(1000),
72 fHmapChannel(0),
73 fKmask(0),
74 fT0width(150.),
75 fT0spreadExt(-1.),
76 fT0fillExt(0)
77{
78 // ctr
79 fCalculated[0] = 0;
80 fCalculated[1] = 0;
81 fCalculated[2] = 0;
82 fCalculated[3] = 0;
83
84 fT0cur[0]=0.;
85 fT0cur[1]=0.;
86
87 if(AliPID::ParticleMass(0) == 0) new AliPID();
88
89 fPIDesd = new AliESDpid();
90
91 fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
92 SetTOFResponse();
93
94 fT0TOF = new AliTOFT0v1(fPIDesd);
95
96}
97//____________________________________________________________________________
98AliTOFT0maker::AliTOFT0maker(AliESDpid *externalPID, AliTOFcalib *tofCalib):
99 TObject(),
100 fT0TOF(NULL),
101 fPIDesd(externalPID),
102 fExternalPIDFlag(kTRUE),
103 fTOFcalib(tofCalib),
104 fNoTOFT0(0),
105 fNmomBins(0),
106 fTimeResolution(100),
107 fT0sigma(1000),
108 fHmapChannel(0),
109 fKmask(0),
110 fT0width(150.),
111 fT0spreadExt(-1.),
112 fT0fillExt(0)
113{
114 // ctr
115 fCalculated[0] = 0;
116 fCalculated[1] = 0;
117 fCalculated[2] = 0;
118 fCalculated[3] = 0;
119
120 fT0cur[0]=0.;
121 fT0cur[1]=0.;
122
123 if(AliPID::ParticleMass(0) == 0) new AliPID();
124
125 if(!fPIDesd){
126 fPIDesd = new AliESDpid();
127 fExternalPIDFlag = kFALSE;
128 }
129
130 fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
131 SetTOFResponse();
132
133 fT0TOF = new AliTOFT0v1(fPIDesd);
134
135}
136
137//____________________________________________________________________________
138AliTOFT0maker::~AliTOFT0maker()
139{
140 // dtor
141
142 delete fT0TOF;
143 if (!fExternalPIDFlag) delete fPIDesd;
144}
145//____________________________________________________________________________
146Double_t* AliTOFT0maker::ComputeT0TOF(AliESDEvent *esd,Double_t t0time,Double_t t0sigma){
147 //
148 // Remake TOF PID probabilities
149 //
150 Double_t t0tof[6];
151
152 if(fKmask) ApplyMask(esd);
153
154 Double_t t0fill = 0.;
155
156 fPIDesd->GetTOFResponse().ResetT0info();
157
158 /* get T0 spread from TOFcalib if available otherwise use default value */
159 if (fTOFcalib && esd) {
160 AliTOFRunParams *runParams = fTOFcalib->GetRunParams();
161 if (runParams && runParams->GetTimestamp(0) != 0) {
162 Float_t t0spread = runParams->EvalT0Spread(esd->GetTimeStamp());
163 if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
164 else{
165 SetT0FillWidth(t0spread);
166 t0fill = fT0fillExt;
167 }
168 }
169 else{
170 if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
171 t0fill = fT0fillExt;
172 }
173 }
174 else if(esd){
175 Float_t t0spread = esd->GetSigma2DiamondZ(); // vertex pread ^2
176 if(t0spread > 0) t0spread = TMath::Sqrt(t0spread)/0.0299792458;
177
178 if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
179 else{
180 SetT0FillWidth(t0spread);
181 t0fill = fT0fillExt;
182 }
183 }
184
185 Float_t thrGood = TMath::Max(Float_t(500.),fT0width*3);
186
187 fT0TOF->Init(esd);
188 AliTOFT0v1* t0maker= fT0TOF;
189
190 t0maker->DefineT0("all",1.5,3.0);
191 t0tof[0] = t0maker->GetResult(0);
192 t0tof[1] = t0maker->GetResult(1);
193 t0tof[2] = t0maker->GetResult(2);
194 t0tof[3] = t0maker->GetResult(3);
195 t0tof[4] = t0maker->GetResult(4);
196 t0tof[5] = t0maker->GetResult(5);
197
198 Float_t lT0Current=0.;
199 fT0sigma=1000;
200
201// Int_t nrun = esd->GetRunNumber();
202
203 t0time += t0fill;
204
205 Float_t sigmaFill = fT0width;
206
207 if(sigmaFill < 20) sigmaFill = 140;
208
209 fCalculated[0]=-1000*t0tof[0]; // best t0
210 fCalculated[1]=1000*t0tof[1]; // sigma best t0
211 fCalculated[2] = t0fill; //t0 fill
212 fCalculated[3] = t0tof[2]; // n TOF tracks
213 fCalculated[4]=-1000*t0tof[0]; // TOF t0
214 fCalculated[5]=1000*t0tof[1]; // TOF t0 sigma
215 fCalculated[6]=sigmaFill; // sigma t0 fill
216 fCalculated[7] = t0tof[3]; // n TOF tracks used for T0
217
218 if(fCalculated[7] > 30) thrGood = 10000000;
219
220 //statistics
221 fCalculated[8] = t0tof[4]; // real time in s
222 fCalculated[9] = t0tof[5]; // cpu time in s
223
224 if(fCalculated[1] < sigmaFill && TMath::Abs(fCalculated[0] - t0fill) < thrGood && fCalculated[1] < fTimeResolution*1.2){
225 fT0sigma=fCalculated[1];
226 lT0Current=fCalculated[0];
227 }
228 else{
229 fCalculated[4] = t0fill;
230 fCalculated[5] = sigmaFill;
231 }
232
233 if(fCalculated[1] < 1 || fT0sigma > sigmaFill || fCalculated[1] > fTimeResolution* 1.2){
234 fT0sigma =1000;
235 fCalculated[4] = t0fill;
236 fCalculated[5] = sigmaFill;
237 }
238
239 if(t0sigma < 1000){
240 if(fT0sigma < 1000){
241 Double_t w1 = 1./t0sigma/t0sigma;
242 Double_t w2 = 1./fCalculated[1]/fCalculated[1];
243
244 Double_t wtot = w1+w2;
245
246 lT0Current = (w1*t0time + w2*fCalculated[0]) / wtot;
247 fT0sigma = TMath::Sqrt(1./wtot);
248 }
249 else{
250 lT0Current=t0time;
251 fT0sigma=t0sigma;
252 }
253 }
254
255 if(fT0sigma < sigmaFill && TMath::Abs(lT0Current - t0fill) < thrGood){
256 fCalculated[1]=fT0sigma;
257 fCalculated[0]=lT0Current;
258 }
259
260 if(fT0sigma >= 1000 || fNoTOFT0){
261 lT0Current = t0fill;
262 fT0sigma = sigmaFill;
263
264 fCalculated[0] = t0fill;
265 fCalculated[1] = sigmaFill;
266 }
267
268 // T0 pt bin
269 Float_t *t0values = new Float_t[fNmomBins];
270 Float_t *t0resolution = new Float_t[fNmomBins];
271 if(fCalculated[7] < 100){
272 for(Int_t i=0;i<fNmomBins;i++){
273 t0maker->DefineT0("all",fPIDesd->GetTOFResponse().GetMinMom(i),fPIDesd->GetTOFResponse().GetMaxMom(i));
274 t0tof[0] = t0maker->GetResult(0);
275 t0tof[1] = t0maker->GetResult(1);
276 t0tof[2] = t0maker->GetResult(2);
277 t0tof[3] = t0maker->GetResult(3);
278
279 Float_t t0bin =-1000*t0tof[0]; // best t0
280 Float_t t0binRes =1000*t0tof[1]; // sigma best t0
281
282 if(t0binRes < sigmaFill && t0binRes < fTimeResolution * 1.2 && TMath::Abs(t0bin - t0fill) < thrGood){
283 // Ok T0
284 if(t0sigma < 1000){
285 Double_t w1 = 1./t0sigma/t0sigma;
286 Double_t w2 = 1./t0binRes/t0binRes;
287
288 Double_t wtot = w1+w2;
289
290 t0bin = (w1*t0time + w2*t0bin) / wtot;
291 t0binRes = TMath::Sqrt(1./wtot);
292 }
293 }
294 else{
295 t0bin = t0fill;
296 t0binRes = sigmaFill;
297 if(t0sigma < 1000){
298 t0bin = t0time;
299 t0binRes= t0sigma;
300 }
301 }
302 t0values[i] = t0bin;
303 t0resolution[i] = t0binRes;
304 }
305 }
306 else{
307 for(Int_t i=0;i<fNmomBins;i++){
308 t0values[i] = lT0Current;
309 t0resolution[i] = fT0sigma;
310 }
311 }
312 for(Int_t i=0;i<fNmomBins;i++){
313 fPIDesd->GetTOFResponse().SetT0bin(i,t0values[i]);
314 fPIDesd->GetTOFResponse().SetT0binRes(i,t0resolution[i]);
315 }
316
317 delete[] t0values;
318 delete[] t0resolution;
319
320 return fCalculated;
321}
322//____________________________________________________________________________
323Double_t *AliTOFT0maker::GetT0p(Float_t p){// [0]=to -- [1] = sigma T0
324 Int_t i=fPIDesd->GetTOFResponse().GetMomBin(p);
325
326 fT0cur[0] = fPIDesd->GetTOFResponse().GetT0bin(i);
327 fT0cur[1] = fPIDesd->GetTOFResponse().GetT0binRes(i);
328
329 return fT0cur;
330}
331//____________________________________________________________________________
332void AliTOFT0maker::SetTOFResponse(){
333 fPIDesd->GetTOFResponse().SetTimeResolution(fTimeResolution);
334}
335//____________________________________________________________________________
336Float_t AliTOFT0maker::GetExpectedSigma(Float_t mom, Float_t tof, Float_t mass){
337 Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,tof,mass);
338
339 return sigma;
340}
341//____________________________________________________________________________
342void AliTOFT0maker::ApplyT0TOF(AliESDEvent *esd){
343 //
344 // Recalculate TOF PID probabilities
345 //
346
347// subtruct t0 for each track
348 Int_t ntracks = esd->GetNumberOfTracks();
349
350 while (ntracks--) {
351 AliESDtrack *t=esd->GetTrack(ntracks);
352
353 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
354
355 Double_t time=t->GetTOFsignal();
356 Float_t p = t->GetP();
357
358 Double_t *t0=GetT0p(p);
359 time -= t0[0];
360 t->SetTOFsignal(time);
361 }
362
363 for(Int_t i=0;i<fNmomBins;i++){
364 fPIDesd->GetTOFResponse().SetT0bin(i,0.0);
365 }
366
367 //
368}
369//____________________________________________________________________________
370void AliTOFT0maker::LoadChannelMap(char *filename){
371 // Load the histo with the channel off map
372 TFile *f= new TFile(filename);
373 if(!f){
374 printf("Cannot open the channel map file (%s)\n",filename);
375 return;
376 }
377
378 fHmapChannel = (TH1F *) f->Get("hChEnabled");
379
380 if(!fHmapChannel){
381 printf("Cannot laod the channel map histo (from %s)\n",filename);
382 return;
383 }
384
385}
386//____________________________________________________________________________
387void AliTOFT0maker::ApplyMask(AliESDEvent * const esd){
388 // Switch off the disable channel
389 if(!fHmapChannel){
390 printf("Channel Map is not available\n");
391 return;
392 }
393
394 Int_t ntracks = esd->GetNumberOfTracks();
395
396 while (ntracks--) {
397 AliESDtrack *t=esd->GetTrack(ntracks);
398
399 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
400
401 Int_t chan = t->GetTOFCalChannel();
402
403 if(fHmapChannel->GetBinContent(chan) < 0.01){
404 t->ResetStatus(AliESDtrack::kTOFout);
405 }
406 }
407}
408
409Float_t
410AliTOFT0maker::TuneForMC(AliESDEvent *esd){ // return true T0 event
411 //
412 // tune for MC data
413 //
414
415 Float_t TOFtimeResolutionDefault=80;
416
417 Float_t t0 = gRandom->Gaus(0.,fT0width);
418
419 Float_t extraSmearing = 0;
420
421 if(fTimeResolution > TOFtimeResolutionDefault){
422 extraSmearing = TMath::Sqrt(fTimeResolution*fTimeResolution - TOFtimeResolutionDefault*TOFtimeResolutionDefault);
423 }
424
425 // subtruct t0 for each track
426 Int_t ntracks = esd->GetNumberOfTracks();
427
428 while (ntracks--) {
429 AliESDtrack *t=esd->GetTrack(ntracks);
430
431 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
432
433 /* check if channel is enabled */
434 if (fTOFcalib){
435 if(!fTOFcalib->IsChannelEnabled(t->GetTOFCalChannel())) {
436 /* reset TOF status */
437 t->ResetStatus(AliESDtrack::kTOFin);
438 t->ResetStatus(AliESDtrack::kTOFout);
439 t->ResetStatus(AliESDtrack::kTOFmismatch);
440 t->ResetStatus(AliESDtrack::kTOFpid);
441 }
442 }
443
444 Double_t time=t->GetTOFsignal();
445
446 time += t0;
447
448 if(extraSmearing>0){
449 Float_t smearing = gRandom->Gaus(0.,extraSmearing);
450 time += smearing;
451 }
452
453 t->SetTOFsignal(time);
454 }
455 //
456 return t0;
457}
458//_________________________________________________________________________
459void AliTOFT0maker::WriteInESD(AliESDEvent *esd){
460 //
461 // Write t0Gen, t0ResGen, nt0;
462 // t0resESD[0:nt0], it0ESD[0:nt0]
463 // in the AliESDEvent
464 //
465 Int_t nMomBins = fPIDesd->GetTOFResponse().GetNmomBins();
466
467 Int_t nt0=0;
468 Float_t *t0 = new Float_t[nMomBins];
469 Float_t *t0res = new Float_t[nMomBins];
470 Int_t *multT0 = new Int_t[nMomBins];
471
472 for(Int_t i=0;i<nMomBins;i++){
473// printf("START %i) %f %f\n",i,fT0event[i],fT0resolution[i]);
474 multT0[i]=0;
475 Bool_t kNewT0=kTRUE;
476 for(Int_t j=0;j < nt0;j++){
477 if(TMath::Abs(fPIDesd->GetTOFResponse().GetT0bin(i) - t0[j])<0.1){
478 kNewT0=kFALSE;
479 multT0[j]++;
480 j=nMomBins*10;
481 }
482 }
483 if(kNewT0){
484 t0[nt0]=fPIDesd->GetTOFResponse().GetT0bin(i);
485 t0res[nt0]=fPIDesd->GetTOFResponse().GetT0binRes(i);
486 nt0++;
487 }
488 }
489
490 Int_t iMultT0=0,nmult=0;
491 for(Int_t j=0;j < nt0;j++){ // find the most frequent
492 if(multT0[j] > nmult){
493 iMultT0 = j;
494 nmult = multT0[j];
495 }
496 }
497
498 Float_t *t0ESD = new Float_t[nMomBins];
499 Float_t *t0resESD = new Float_t[nMomBins];
500 Int_t *it0ESD = new Int_t[nMomBins];
501
502 Float_t t0Gen,t0ResGen;
503 t0Gen = t0[iMultT0];
504 t0ResGen = t0res[iMultT0];
505 nt0=0;
506 // printf("T0 to ESD\n%f %f\n",t0Gen,t0ResGen);
507 for(Int_t i=0;i<nMomBins;i++){
508 if(TMath::Abs(fPIDesd->GetTOFResponse().GetT0bin(i) - t0Gen)>0.1){
509 t0ESD[nt0]=fPIDesd->GetTOFResponse().GetT0bin(i);
510 t0resESD[nt0]=fPIDesd->GetTOFResponse().GetT0binRes(i);
511 it0ESD[nt0]=i;
512// printf("%i) %f %f %i\n",nt0,t0ESD[nt0],t0resESD[nt0],it0ESD[nt0]);
513 nt0++;
514 }
515 }
516
517 // Write t0Gen,t0ResGen; nt0; t0resESD[0:nt0],it0ESD[0:nt0] in the AliESDEvent
518
519 AliTOFHeader *tofHeader =
520 new AliTOFHeader(t0Gen,t0ResGen,nt0,
521 t0ESD,t0resESD,it0ESD,fTimeResolution,fT0width);
522
523 esd->SetTOFHeader(tofHeader);
524
525 AliDebug(1,Form("resTOF=%f T0spread=%f t0Gen=%f t0resGen=%f",fTimeResolution,fT0width,t0Gen,t0ResGen));
526 AliDebug(1,Form("%d ",nt0));
527 for (Int_t ii=0; ii<nt0; ii++)
528 AliDebug(1,Form("pBin=%d t0val=%f t0res=%f",it0ESD[ii],t0ESD[ii],t0resESD[ii]));
529
530 delete[] t0;
531 t0 = NULL;
532 delete[] t0res;
533 t0res = NULL;
534 delete[] multT0;
535 multT0 = NULL;
536 delete[] t0ESD;
537 t0ESD = NULL;
538 delete[] t0resESD;
539 t0resESD = NULL;
540 delete[] it0ESD;
541 it0ESD = NULL;
542
543}