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