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