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