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