TOF/AliTOFT0maker.cxx

   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 ClassImp(AliTOFT0maker)
  60
  61 //____________________________________________________________________________
  62 AliTOFT0maker::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 //____________________________________________________________________________
  99 AliTOFT0maker::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 //____________________________________________________________________________
 140 AliTOFT0maker::~AliTOFT0maker()
 141 {
 142   // dtor
 143   delete fT0TOF;
 144   if (!fExternalPIDFlag) delete fPIDesd;
 145 }
 146 //____________________________________________________________________________
 147 Double_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 //____________________________________________________________________________
 325 Double_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 //____________________________________________________________________________
 334 void AliTOFT0maker::SetTOFResponse(){
 335   fPIDesd->GetTOFResponse().SetTimeResolution(fTimeResolution);
 336 }
 337 //____________________________________________________________________________
 338 Float_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 //____________________________________________________________________________
 344 void 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 //____________________________________________________________________________
 372 void  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 //____________________________________________________________________________
 389 void 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
 411 Float_t
 412 AliTOFT0maker::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 //_________________________________________________________________________
 461 void     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 }