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 "AliESDpid.h"
  50 #include "AliESDEvent.h"
  51 #include "TFile.h"
  52 #include "TH1F.h"
  53 #include "AliTOFcalib.h"
  54 #include "AliTOFRunParams.h"
  55 #include "TRandom.h"
  56
  57 ClassImp(AliTOFT0maker)
  58
  59 //____________________________________________________________________________
  60 AliTOFT0maker::AliTOFT0maker():
  61   TObject(),
  62   fT0TOF(NULL),
  63   fPIDesd(NULL),
  64   fExternalPIDFlag(kFALSE),
  65   fTOFcalib(NULL),
  66   fNoTOFT0(0),
  67   fTimeResolution(100),
  68   fT0sigma(1000),
  69   fHmapChannel(0),
  70   fKmask(0),
  71   fT0width(150.)
  72 {
  73   // ctr
  74   fCalculated[0] = 0;
  75   fCalculated[1] = 0;
  76   fCalculated[2] = 0;
  77   fCalculated[3] = 0;
  78
  79   fT0TOF = new AliTOFT0v1();
  80   if(AliPID::ParticleMass(0) == 0) new AliPID();
  81
  82   fPIDesd = new AliESDpid();
  83
  84   fPtCutMin[0] = 0.3;
  85   fPtCutMin[1] = 0.5;
  86   fPtCutMin[2] = 0.6;
  87   fPtCutMin[3] = 0.7;
  88   fPtCutMin[4] = 0.8;
  89   fPtCutMin[5] = 0.9;
  90   fPtCutMin[6] = 1;
  91   fPtCutMin[7] = 1.2;
  92   fPtCutMin[8] = 1.5;
  93   fPtCutMin[9] = 2;
  94
  95   fPtCutMax[0] = 0.5;
  96   fPtCutMax[1] = 0.6;
  97   fPtCutMax[2] = 0.7;
  98   fPtCutMax[3] = 0.8;
  99   fPtCutMax[4] = 0.9;
 100   fPtCutMax[5] = 1;
 101   fPtCutMax[6] = 1.2;
 102   fPtCutMax[7] = 1.5;
 103   fPtCutMax[8] = 2;
 104   fPtCutMax[9] = 3;
 105
 106   /* init arrays */
 107   for (Int_t i = 0; i < 10; i++) {
 108     fT0pt[i] = 0.;
 109     fT0ptSigma[i] = 0.;
 110   }
 111 }
 112 //____________________________________________________________________________
 113 AliTOFT0maker::AliTOFT0maker(AliESDpid *externalPID, AliTOFcalib *tofCalib):
 114     TObject(),
 115     fT0TOF(NULL),
 116     fPIDesd(externalPID),
 117     fExternalPIDFlag(kTRUE),
 118     fTOFcalib(tofCalib),
 119     fNoTOFT0(0),
 120     fTimeResolution(100),
 121     fT0sigma(1000),
 122     fHmapChannel(0),
 123     fKmask(0),
 124     fT0width(150.)
 125 {
 126   // ctr
 127   fCalculated[0] = 0;
 128   fCalculated[1] = 0;
 129   fCalculated[2] = 0;
 130   fCalculated[3] = 0;
 131
 132   fT0TOF = new AliTOFT0v1();
 133   if(AliPID::ParticleMass(0) == 0) new AliPID();
 134
 135   if(!fPIDesd){
 136     fPIDesd = new AliESDpid();
 137     fExternalPIDFlag = kFALSE;
 138   }
 139
 140   fPtCutMin[0] = 0.3;
 141   fPtCutMin[1] = 0.5;
 142   fPtCutMin[2] = 0.6;
 143   fPtCutMin[3] = 0.7;
 144   fPtCutMin[4] = 0.8;
 145   fPtCutMin[5] = 0.9;
 146   fPtCutMin[6] = 1;
 147   fPtCutMin[7] = 1.2;
 148   fPtCutMin[8] = 1.5;
 149   fPtCutMin[9] = 2;
 150
 151   fPtCutMax[0] = 0.5;
 152   fPtCutMax[1] = 0.6;
 153   fPtCutMax[2] = 0.7;
 154   fPtCutMax[3] = 0.8;
 155   fPtCutMax[4] = 0.9;
 156   fPtCutMax[5] = 1;
 157   fPtCutMax[6] = 1.2;
 158   fPtCutMax[7] = 1.5;
 159   fPtCutMax[8] = 2;
 160   fPtCutMax[9] = 3;
 161
 162   /* init arrays */
 163   for (Int_t i = 0; i < 10; i++) {
 164     fT0pt[i] = 0.;
 165     fT0ptSigma[i] = 0.;
 166   }
 167 }
 168
 169 /* copy-constructor and operator= suppressed
 170
 171
 172 //____________________________________________________________________________
 173 AliTOFT0maker::AliTOFT0maker(const AliTOFT0maker & t) :
 174   TObject(t),
 175   fT0TOF(t.fT0TOF),
 176   fPIDESD(t.fPIDESD),
 177   fNoTOFT0(t.fNoTOFT0),
 178   fTimeResolution(t.fTimeResolution),
 179   fT0sigma(t.fT0sigma),
 180   fHmapChannel(t.fHmapChannel),
 181   fKmask(t.fKmask)
 182   {
 183   // copy ctr
 184 }
 185
 186 //____________________________________________________________________________
 187 AliTOFT0maker& AliTOFT0maker::operator=(const AliTOFT0maker &t)
 188 {
 189   //
 190   // assign. operator
 191   //
 192
 193   if (this == &t)
 194     return *this;
 195
 196   TObject::operator=(t);
 197   fTimeResolution = t.fTimeResolution;
 198   fT0sigma = t.fT0sigma;
 199
 200   return *this;
 201 }
 202
 203 */
 204
 205 //____________________________________________________________________________
 206 AliTOFT0maker::~AliTOFT0maker()
 207 {
 208   // dtor
 209
 210   delete fT0TOF;
 211   if (!fExternalPIDFlag) delete fPIDesd;
 212 }
 213 //____________________________________________________________________________
 214 Double_t* AliTOFT0maker::ComputeT0TOF(AliESDEvent *esd,Double_t t0time,Double_t t0sigma){
 215   //
 216   // Remake TOF PID probabilities
 217   //
 218
 219   Double_t t0tof[4];
 220
 221   if(fKmask) ApplyMask(esd);
 222
 223   /* get T0 spread from TOFcalib if available otherwise use default value */
 224   if (fTOFcalib && esd) {
 225     AliTOFRunParams *runParams = fTOFcalib->GetRunParams();
 226     if (runParams && runParams->GetTimestamp(0) != 0) {
 227       Float_t t0spread = runParams->EvalT0Spread(esd->GetTimeStamp());
 228       SetT0FillWidth(t0spread);
 229     }
 230   }
 231
 232   fT0TOF->Init(esd);
 233   AliTOFT0v1* t0maker= fT0TOF;
 234   t0maker->SetTimeResolution(fTimeResolution*1e-12);
 235
 236   t0maker->DefineT0("all",1.5,3.0);
 237   t0tof[0] = t0maker->GetResult(0);
 238   t0tof[1] = t0maker->GetResult(1);
 239   t0tof[2] = t0maker->GetResult(2);
 240   t0tof[3] = t0maker->GetResult(3);
 241
 242   Float_t lT0Current=0.;
 243   fT0sigma=1000;
 244
 245 //   Int_t nrun = esd->GetRunNumber();
 246   Double_t t0fill = 0.;
 247
 248   t0time += t0fill;
 249
 250   Float_t sigmaFill = fT0width;
 251
 252   if(sigmaFill < 20) sigmaFill = 140;
 253
 254   fCalculated[0]=-1000*t0tof[0]; // best t0
 255   fCalculated[1]=1000*t0tof[1]; // sigma best t0
 256   fCalculated[2] = t0fill;    //t0 fill
 257   fCalculated[3] = t0tof[2];  // n TOF tracks
 258   fCalculated[4]=-1000*t0tof[0]; // TOF t0
 259   fCalculated[5]=1000*t0tof[1]; // TOF t0 sigma
 260   fCalculated[6]=sigmaFill; // sigma t0 fill
 261   fCalculated[7] = t0tof[3];  // n TOF tracks used for T0
 262
 263   if(fCalculated[1] < sigmaFill && TMath::Abs(fCalculated[0] - t0fill) < 500 && fCalculated[1] < fTimeResolution*1.2){
 264     fT0sigma=fCalculated[1];
 265     lT0Current=fCalculated[0];
 266   }
 267   else{
 268     fCalculated[4] = t0fill;
 269     fCalculated[5] = sigmaFill;
 270   }
 271
 272   if(fCalculated[1] < 1 || fT0sigma > sigmaFill || fCalculated[1] > fTimeResolution* 1.2){
 273     fT0sigma =1000;
 274     fCalculated[4] = t0fill;
 275     fCalculated[5] = sigmaFill;
 276   }
 277
 278   if(t0sigma < 1000){
 279     if(fT0sigma < 1000){
 280       Double_t w1 = 1./t0sigma/t0sigma;
 281       Double_t w2 = 1./fCalculated[1]/fCalculated[1];
 282
 283       Double_t wtot = w1+w2;
 284
 285       lT0Current = (w1*t0time + w2*fCalculated[0]) / wtot;
 286       fT0sigma = TMath::Sqrt(1./wtot);
 287     }
 288     else{
 289       lT0Current=t0time;
 290       fT0sigma=t0sigma;
 291     }
 292   }
 293
 294   if(fT0sigma < sigmaFill && TMath::Abs(lT0Current - t0fill) < 500){
 295     fCalculated[1]=fT0sigma;
 296     fCalculated[0]=lT0Current;
 297   }
 298
 299   if(fT0sigma >= 1000 || fNoTOFT0){
 300     lT0Current = t0fill;
 301     fT0sigma = sigmaFill;
 302
 303     fCalculated[0] = t0fill;
 304     fCalculated[1] = sigmaFill;
 305   }
 306
 307   // T0 pt bin
 308   for(Int_t i=0;i<10;i++){
 309    t0maker->DefineT0("all",fPtCutMin[i],fPtCutMax[i]);
 310     t0tof[0] = t0maker->GetResult(0);
 311     t0tof[1] = t0maker->GetResult(1);
 312     t0tof[2] = t0maker->GetResult(2);
 313     t0tof[3] = t0maker->GetResult(3);
 314     fT0pt[i] =-1000*t0tof[0]; // best t0
 315     fT0ptSigma[i] =1000*t0tof[1]; // sigma best t0
 316
 317     if(fT0ptSigma[i] < sigmaFill  && fT0ptSigma[i] < fTimeResolution * 1.2 && TMath::Abs(fT0pt[i] - t0fill) < 500){
 318       // Ok T0
 319     }
 320     else{
 321       fT0pt[i] = t0fill;
 322       fT0ptSigma[i] = sigmaFill;
 323     }
 324   }
 325   //----
 326   SetTOFResponse();
 327
 328   return fCalculated;
 329 }
 330 //____________________________________________________________________________
 331 Double_t  *AliTOFT0maker::GetT0p(Float_t p){// [0]=to -- [1] = sigma T0
 332   Int_t i=0;
 333   while(p > fPtCutMin[i] && i < 10) i++;
 334   if(i > 0) i--;
 335
 336   fT0cur[0] = fT0pt[i];
 337   fT0cur[1] = fT0ptSigma[i];
 338   return fT0cur;
 339 }
 340 //____________________________________________________________________________
 341 void AliTOFT0maker::SetTOFResponse(){
 342     fPIDesd->GetTOFResponse().SetTimeResolution(TMath::Sqrt(fT0sigma*fT0sigma + fTimeResolution*fTimeResolution));
 343 }
 344 //____________________________________________________________________________
 345 Float_t AliTOFT0maker::GetExpectedSigma(Float_t mom, Float_t tof, Float_t mass){
 346   Double_t *sigmaT0 = GetT0p(mom);
 347   fPIDesd->GetTOFResponse().SetTimeResolution(TMath::Sqrt(sigmaT0[1]*sigmaT0[1] + fTimeResolution*fTimeResolution));
 348   Float_t sigma = fPIDesd->GetTOFResponse().GetExpectedSigma(mom,tof,mass);
 349   fPIDesd->GetTOFResponse().SetTimeResolution(TMath::Sqrt(fT0sigma*fT0sigma + fTimeResolution*fTimeResolution));
 350
 351   return sigma;
 352 }
 353 //____________________________________________________________________________
 354 void AliTOFT0maker::ApplyT0TOF(AliESDEvent *esd){
 355   //
 356   // Recalculate TOF PID probabilities
 357   //
 358
 359   // subtruct t0 for each track
 360   Int_t ntracks = esd->GetNumberOfTracks();
 361
 362   while (ntracks--) {
 363     AliESDtrack *t=esd->GetTrack(ntracks);
 364
 365     if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
 366
 367     Double_t time=t->GetTOFsignal();
 368     Float_t p = t->GetP();
 369
 370     Double_t *t0=GetT0p(p);
 371     time -= t0[0];
 372     t->SetTOFsignal(time);
 373   }
 374   //
 375 }
 376 //____________________________________________________________________________
 377 void  AliTOFT0maker::LoadChannelMap(char *filename){
 378   // Load the histo with the channel off map
 379   TFile *f= new TFile(filename);
 380   if(!f){
 381     printf("Cannot open the channel map file (%s)\n",filename);
 382     return;
 383   }
 384
 385   fHmapChannel = (TH1F *) f->Get("hChEnabled");
 386
 387   if(!fHmapChannel){
 388     printf("Cannot laod the channel map histo (from %s)\n",filename);
 389     return;
 390   }
 391
 392 }
 393 //____________________________________________________________________________
 394 void AliTOFT0maker::ApplyMask(AliESDEvent * const esd){
 395   // Switch off the disable channel
 396   if(!fHmapChannel){
 397     printf("Channel Map is not available\n");
 398     return;
 399   }
 400
 401   Int_t ntracks = esd->GetNumberOfTracks();
 402
 403   while (ntracks--) {
 404     AliESDtrack *t=esd->GetTrack(ntracks);
 405
 406     if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
 407
 408     Int_t chan = t->GetTOFCalChannel();
 409
 410     if(fHmapChannel->GetBinContent(chan) < 0.01){
 411       t->ResetStatus(AliESDtrack::kTOFout);
 412     }
 413   }
 414 }
 415
 416 Float_t
 417 AliTOFT0maker::TuneForMC(AliESDEvent *esd){ // return true T0 event
 418   //
 419   // tune for MC data
 420   //
 421
 422   Float_t TOFtimeResolutionDefault=80;
 423
 424   Float_t t0 = gRandom->Gaus(0.,fT0width);
 425
 426   Float_t extraSmearing = 0;
 427
 428   if(fTimeResolution > TOFtimeResolutionDefault){
 429     extraSmearing = TMath::Sqrt(fTimeResolution*fTimeResolution - TOFtimeResolutionDefault*TOFtimeResolutionDefault);
 430   }
 431
 432   // subtruct t0 for each track
 433   Int_t ntracks = esd->GetNumberOfTracks();
 434
 435   while (ntracks--) {
 436     AliESDtrack *t=esd->GetTrack(ntracks);
 437
 438     if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
 439
 440     /* check if channel is enabled */
 441     if (fTOFcalib && !fTOFcalib->IsChannelEnabled(t->GetTOFCalChannel())) {
 442       /* reset TOF status */
 443       t->ResetStatus(AliESDtrack::kTOFin);
 444       t->ResetStatus(AliESDtrack::kTOFout);
 445       t->ResetStatus(AliESDtrack::kTOFrefit);
 446       t->ResetStatus(AliESDtrack::kTOFpid);
 447     }
 448
 449     Double_t time=t->GetTOFsignal();
 450
 451     time += t0;
 452
 453     if(extraSmearing>0){
 454       Float_t smearing = gRandom->Gaus(0.,extraSmearing);
 455       time += smearing;
 456     }
 457
 458     t->SetTOFsignal(time);
 459   }
 460   //
 461   return t0;
 462 }