Improved QAchecker and plots for DQM
[u/mrichter/AliRoot.git] / TOF / AliTOFT0maker.cxx
CommitLineData
8f589502 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 **************************************************************************/
6a4e212e 15/* $Id: AliTOFT0maker.cxx,v 1.8 2010/01/19 16:32:20 noferini Exp $ */
8f589502 16
17/////////////////////////////////////////////////////////////////////////////
18// //
19// This class contains the basic functions for the time zero //
6a4e212e 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) //
5b4ed716 24// AliESDpid *extPID=new AliESDpid(); //
25// fTOFmaker = new AliTOFT0maker(extPID); //
26// fTOFmaker->SetTimeResolution(100.0); // if you want set the TOF res //
6a4e212e 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 //
03bd764d 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 //
6a4e212e 40// //
41// Let consider that: //
42// - the PIF is automatically recalculated with the event time subtrction //
8f589502 43// //
44/////////////////////////////////////////////////////////////////////////////
45
69ccfdf6 46#include "AliTOFT0v1.h"
536031f2 47#include "AliTOFT0maker.h"
536031f2 48#include "AliPID.h"
f858b00e 49#include "AliLog.h"
536031f2 50#include "AliESDpid.h"
03bd764d 51#include "AliESDEvent.h"
52#include "TFile.h"
53#include "TH1F.h"
5b4ed716 54#include "AliTOFcalib.h"
55#include "AliTOFRunParams.h"
56#include "TRandom.h"
f858b00e 57#include "AliTOFHeader.h"
536031f2 58
59ClassImp(AliTOFT0maker)
60
61//____________________________________________________________________________
5b4ed716 62AliTOFT0maker::AliTOFT0maker():
63 TObject(),
64 fT0TOF(NULL),
65 fPIDesd(NULL),
66 fExternalPIDFlag(kFALSE),
67 fTOFcalib(NULL),
68 fNoTOFT0(0),
2a258f40 69 fNmomBins(0),
5b4ed716 70 fTimeResolution(100),
71 fT0sigma(1000),
72 fHmapChannel(0),
73 fKmask(0),
2a258f40 74 fT0width(150.),
75 fT0spreadExt(-1.),
a558aa69 76 fT0fillExt(0),
77 fTOFT0algorithm(1)
5b4ed716 78{
79 // ctr
80 fCalculated[0] = 0;
81 fCalculated[1] = 0;
82 fCalculated[2] = 0;
83 fCalculated[3] = 0;
84
020dbfc2 85 fT0cur[0]=0.;
86 fT0cur[1]=0.;
87
5b4ed716 88 if(AliPID::ParticleMass(0) == 0) new AliPID();
89
90 fPIDesd = new AliESDpid();
91
2a258f40 92 fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
93 SetTOFResponse();
94
95 fT0TOF = new AliTOFT0v1(fPIDesd);
96
5b4ed716 97}
98//____________________________________________________________________________
99AliTOFT0maker::AliTOFT0maker(AliESDpid *externalPID, AliTOFcalib *tofCalib):
100 TObject(),
101 fT0TOF(NULL),
102 fPIDesd(externalPID),
103 fExternalPIDFlag(kTRUE),
104 fTOFcalib(tofCalib),
03bd764d 105 fNoTOFT0(0),
2a258f40 106 fNmomBins(0),
5b4ed716 107 fTimeResolution(100),
03bd764d 108 fT0sigma(1000),
109 fHmapChannel(0),
5b4ed716 110 fKmask(0),
2a258f40 111 fT0width(150.),
112 fT0spreadExt(-1.),
a558aa69 113 fT0fillExt(0),
114 fTOFT0algorithm(1)
536031f2 115{
6a4e212e 116 // ctr
6a4e212e 117 fCalculated[0] = 0;
118 fCalculated[1] = 0;
119 fCalculated[2] = 0;
03bd764d 120 fCalculated[3] = 0;
6a4e212e 121
020dbfc2 122 fT0cur[0]=0.;
123 fT0cur[1]=0.;
124
536031f2 125 if(AliPID::ParticleMass(0) == 0) new AliPID();
03bd764d 126
5b4ed716 127 if(!fPIDesd){
128 fPIDesd = new AliESDpid();
129 fExternalPIDFlag = kFALSE;
130 }
131
2a258f40 132 fNmomBins = fPIDesd->GetTOFResponse().GetNmomBins();
133 SetTOFResponse();
5b4ed716 134
2a258f40 135 fT0TOF = new AliTOFT0v1(fPIDesd);
8f589502 136
8f589502 137}
5b4ed716 138
8f589502 139//____________________________________________________________________________
536031f2 140AliTOFT0maker::~AliTOFT0maker()
141{
142 // dtor
5b4ed716 143 delete fT0TOF;
144 if (!fExternalPIDFlag) delete fPIDesd;
536031f2 145}
146//____________________________________________________________________________
5b4ed716 147Double_t* AliTOFT0maker::ComputeT0TOF(AliESDEvent *esd,Double_t t0time,Double_t t0sigma){
8f589502 148 //
149 // Remake TOF PID probabilities
150 //
2a258f40 151 Double_t t0tof[6];
536031f2 152
03bd764d 153 if(fKmask) ApplyMask(esd);
154
2a258f40 155 Double_t t0fill = 0.;
156
157 fPIDesd->GetTOFResponse().ResetT0info();
158
5b4ed716 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());
2a258f40 164 if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
165 else{
166 SetT0FillWidth(t0spread);
167 t0fill = fT0fillExt;
168 }
169 }
170 else{
90beec49 171 if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
2a258f40 172 t0fill = fT0fillExt;
5b4ed716 173 }
536031f2 174 }
bc4fdef6 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 }
03bd764d 185
bf40b7fa 186 Float_t thrGood = TMath::Max(Float_t(500.),fT0width*3);
187
5b4ed716 188
a558aa69 189 fT0TOF->Init(esd);
190 AliTOFT0v1* t0maker = fT0TOF;
191 if (fTOFT0algorithm==2) t0maker->SetOptimization(kTRUE);
5b4ed716 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);
2a258f40 197 t0tof[4] = t0maker->GetResult(4);
198 t0tof[5] = t0maker->GetResult(5);
536031f2 199
8f589502 200 Float_t lT0Current=0.;
536031f2 201 fT0sigma=1000;
202
03bd764d 203// Int_t nrun = esd->GetRunNumber();
5b4ed716 204
03bd764d 205 t0time += t0fill;
8f589502 206
5b4ed716 207 Float_t sigmaFill = fT0width;
208
209 if(sigmaFill < 20) sigmaFill = 140;
8f589502 210
03bd764d 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
bf40b7fa 220 if(fCalculated[7] > 30) thrGood = 10000000;
221
2a258f40 222 //statistics
223 fCalculated[8] = t0tof[4]; // real time in s
224 fCalculated[9] = t0tof[5]; // cpu time in s
225
bf40b7fa 226 if(fCalculated[1] < sigmaFill && TMath::Abs(fCalculated[0] - t0fill) < thrGood && fCalculated[1] < fTimeResolution*1.2){
6a4e212e 227 fT0sigma=fCalculated[1];
228 lT0Current=fCalculated[0];
536031f2 229 }
03bd764d 230 else{
231 fCalculated[4] = t0fill;
232 fCalculated[5] = sigmaFill;
233 }
234
5b4ed716 235 if(fCalculated[1] < 1 || fT0sigma > sigmaFill || fCalculated[1] > fTimeResolution* 1.2){
03bd764d 236 fT0sigma =1000;
237 fCalculated[4] = t0fill;
238 fCalculated[5] = sigmaFill;
239 }
536031f2 240
241 if(t0sigma < 1000){
242 if(fT0sigma < 1000){
243 Double_t w1 = 1./t0sigma/t0sigma;
6a4e212e 244 Double_t w2 = 1./fCalculated[1]/fCalculated[1];
536031f2 245
246 Double_t wtot = w1+w2;
247
6a4e212e 248 lT0Current = (w1*t0time + w2*fCalculated[0]) / wtot;
536031f2 249 fT0sigma = TMath::Sqrt(1./wtot);
250 }
251 else{
8f589502 252 lT0Current=t0time;
536031f2 253 fT0sigma=t0sigma;
254 }
255 }
256
bf40b7fa 257 if(fT0sigma < sigmaFill && TMath::Abs(lT0Current - t0fill) < thrGood){
03bd764d 258 fCalculated[1]=fT0sigma;
259 fCalculated[0]=lT0Current;
260 }
261
262 if(fT0sigma >= 1000 || fNoTOFT0){
8f589502 263 lT0Current = t0fill;
03bd764d 264 fT0sigma = sigmaFill;
8f589502 265
6a4e212e 266 fCalculated[0] = t0fill;
03bd764d 267 fCalculated[1] = sigmaFill;
536031f2 268 }
269
5b4ed716 270 // T0 pt bin
721ed687 271 Float_t *t0values = new Float_t[fNmomBins];
272 Float_t *t0resolution = new Float_t[fNmomBins];
2a258f40 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);
2a258f40 280
281 Float_t t0bin =-1000*t0tof[0]; // best t0
282 Float_t t0binRes =1000*t0tof[1]; // sigma best t0
283
bf40b7fa 284 if(t0binRes < sigmaFill && t0binRes < fTimeResolution * 1.2 && TMath::Abs(t0bin - t0fill) < thrGood){
2a258f40 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 }
721ed687 304 t0values[i] = t0bin;
305 t0resolution[i] = t0binRes;
5b4ed716 306 }
2a258f40 307 }
308 else{
309 for(Int_t i=0;i<fNmomBins;i++){
721ed687 310 t0values[i] = lT0Current;
311 t0resolution[i] = fT0sigma;
5b4ed716 312 }
313 }
721ed687 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;
536031f2 321
6a4e212e 322 return fCalculated;
536031f2 323}
324//____________________________________________________________________________
5b4ed716 325Double_t *AliTOFT0maker::GetT0p(Float_t p){// [0]=to -- [1] = sigma T0
2a258f40 326 Int_t i=fPIDesd->GetTOFResponse().GetMomBin(p);
03bd764d 327
2a258f40 328 fT0cur[0] = fPIDesd->GetTOFResponse().GetT0bin(i);
329 fT0cur[1] = fPIDesd->GetTOFResponse().GetT0binRes(i);
330
5b4ed716 331 return fT0cur;
536031f2 332}
333//____________________________________________________________________________
5b4ed716 334void AliTOFT0maker::SetTOFResponse(){
2a258f40 335 fPIDesd->GetTOFResponse().SetTimeResolution(fTimeResolution);
5b4ed716 336}
337//____________________________________________________________________________
338Float_t AliTOFT0maker::GetExpectedSigma(Float_t mom, Float_t tof, Float_t mass){
5b4ed716 339 Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,tof,mass);
8f589502 340
5b4ed716 341 return sigma;
536031f2 342}
343//____________________________________________________________________________
5b4ed716 344void AliTOFT0maker::ApplyT0TOF(AliESDEvent *esd){
8f589502 345 //
5b4ed716 346 // Recalculate TOF PID probabilities
8f589502 347 //
348
90beec49 349// subtruct t0 for each track
350 Int_t ntracks = esd->GetNumberOfTracks();
536031f2 351
90beec49 352 while (ntracks--) {
353 AliESDtrack *t=esd->GetTrack(ntracks);
03bd764d 354
90beec49 355 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
5b4ed716 356
90beec49 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
5b4ed716 369 //
536031f2 370}
03bd764d 371//____________________________________________________________________________
777fcfde 372void AliTOFT0maker::LoadChannelMap(const char *filename){
03bd764d 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
5b4ed716 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 */
bc4fdef6 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 }
5b4ed716 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}
f858b00e 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,
c85a1b24 523 t0ESD,t0resESD,it0ESD,fTimeResolution,fT0width);
f858b00e 524
525 esd->SetTOFHeader(tofHeader);
526
634de978 527 delete tofHeader;
528
c85a1b24 529 AliDebug(1,Form("resTOF=%f T0spread=%f t0Gen=%f t0resGen=%f",fTimeResolution,fT0width,t0Gen,t0ResGen));
530 AliDebug(1,Form("%d ",nt0));
f858b00e 531 for (Int_t ii=0; ii<nt0; ii++)
c85a1b24 532 AliDebug(1,Form("pBin=%d t0val=%f t0res=%f",it0ESD[ii],t0ESD[ii],t0resESD[ii]));
f858b00e 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}