TENDER/TenderSupplies/AliTOFTenderSupply.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
  16
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 // TOF tender: for a description of what this tender does see
  20 //             https://twiki.cern.ch/twiki/bin/viewauth/ALICE/TOF
  21 //             It has to be used on LHC2010 data (and old MC productions)
  22 //
  23 // Contacts: Pietro.Antonioli@bo.infn.it                                     //
  24 //           Francesco.Noferini@bo.infn.it                                   //
  25 ///////////////////////////////////////////////////////////////////////////////
  26 #include <TFile.h>
  27 #include <TChain.h>
  28
  29 #include <TMath.h>
  30 #include <TRandom.h>
  31 #include <AliLog.h>
  32 #include <AliESDEvent.h>
  33 #include <AliESDtrack.h>
  34 #include <AliESDInputHandler.h>
  35 #include <AliAnalysisManager.h>
  36 #include <AliESDpid.h>
  37 #include <AliTender.h>
  38
  39 #include <AliTOFcalib.h>
  40 #include <AliTOFT0maker.h>
  41
  42 #include <AliGeomManager.h>
  43 #include <AliCDBManager.h>
  44 #include <AliCDBEntry.h>
  45
  46 #include <AliOADBContainer.h>
  47 #include <AliTOFPIDParams.h>
  48
  49 #include <AliT0CalibSeasonTimeShift.h>
  50
  51 #include <AliMultiplicity.h>
  52
  53 #include "AliTOFTenderSupply.h"
  54
  55 ClassImp(AliTOFTenderSupply)
  56
  57 Float_t AliTOFTenderSupply::fgT0Aresolution = 75.;
  58 Float_t AliTOFTenderSupply::fgT0Cresolution = 65.;
  59
  60 AliTOFTenderSupply::AliTOFTenderSupply() :
  61   AliTenderSupply(),
  62   fESDpid(0x0),
  63   fTenderNoAction(kTRUE),
  64   fIsMC(kFALSE),
  65   fCorrectExpTimes(kTRUE),
  66   fCorrectTRDBug(kFALSE),
  67   fLHC10dPatch(kFALSE),
  68   fT0DetectorAdjust(kFALSE),
  69   fDebugLevel(0),
  70   fAutomaticSettings(kTRUE),
  71   fRecoPass(0),
  72   fUserRecoPass(0),
  73   fForceCorrectTRDBug(kFALSE),
  74   fT0Simulate(kFALSE),
  75   fTOFPIDParams(0x0),
  76   fTOFCalib(0x0),
  77   fTOFT0maker(0x0),
  78   fT0IntercalibrationShift(0),
  79   fGeomSet(kFALSE),
  80   fIsEnteringInTRD(kFALSE),
  81   fInTRD(kFALSE),
  82   fIsComingOutTRD(kFALSE),
  83   fOutTRD(kFALSE),
  84   fRhoTRDin(288.38), // cm
  85   fRhoTRDout(366.38), // cm
  86   fStep(0.5),
  87   fMagField(0.),
  88   fCDBkey(0)
  89
  90
  91
  92 {
  93   //
  94   // default ctor
  95   //
  96   fT0shift[0] = 0;
  97   fT0shift[1] = 0;
  98   fT0shift[2] = 0;
  99   fT0shift[3] = 0;
 100 }
 101
 102 //_____________________________________________________
 103 AliTOFTenderSupply::AliTOFTenderSupply(const char *name, const AliTender *tender) :
 104   AliTenderSupply(name,tender),
 105   fESDpid(0x0),
 106   fTenderNoAction(kTRUE),
 107   fIsMC(kFALSE),
 108   fCorrectExpTimes(kTRUE),
 109   fCorrectTRDBug(kFALSE),
 110   fLHC10dPatch(kFALSE),
 111   fT0DetectorAdjust(kFALSE),
 112   fDebugLevel(0),
 113   fAutomaticSettings(kTRUE),
 114   fRecoPass(0),
 115   fUserRecoPass(0),
 116   fForceCorrectTRDBug(kFALSE),
 117   fT0Simulate(kFALSE),
 118   fTOFPIDParams(0x0),
 119   fTOFCalib(0x0),
 120   fTOFT0maker(0x0),
 121   fT0IntercalibrationShift(0),
 122   fGeomSet(kFALSE),
 123   fIsEnteringInTRD(kFALSE),
 124   fInTRD(kFALSE),
 125   fIsComingOutTRD(kFALSE),
 126   fOutTRD(kFALSE),
 127   fRhoTRDin(288.38), // cm
 128   fRhoTRDout(366.38), // cm
 129   fStep(0.5),
 130   fMagField(0.),
 131   fCDBkey(0)
 132
 133 {
 134   //
 135   // named ctor
 136   //
 137
 138   fT0shift[0] = 0;
 139   fT0shift[1] = 0;
 140   fT0shift[2] = 0;
 141   fT0shift[3] = 0;
 142 }
 143
 144 //_____________________________________________________
 145 void AliTOFTenderSupply::Init()
 146 {
 147
 148   fTenderNoAction = kFALSE;
 149   // Initialise TOF tender (this is called at each detected run change)
 150   AliLog::SetClassDebugLevel("AliTOFTenderSupply",10);
 151
 152   // Setup PID object, check for MC, set AliTOFcalib and TOFT0 maker conf
 153   Int_t run = fTender->GetRun();
 154   if (run == 0) return;                // to skip first init, when we don't have yet a run number
 155
 156   fGeomSet=kFALSE;
 157   // Even if the user didn't set fIsMC, we force it on if we find the MC handler
 158   AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
 159   if (mgr->GetMCtruthEventHandler() && !(fIsMC) ) {
 160     AliWarning("This ESD is MC, fIsMC found OFF: fIsMC turned ON");
 161     fIsMC=kTRUE;
 162   }
 163
 164   if (fAutomaticSettings) {
 165     if (!fIsMC) {
 166       if (fUserRecoPass == 0) DetectRecoPass();
 167       else fRecoPass = fUserRecoPass;
 168     }
 169     if (run<114737) {
 170       fTenderNoAction = kTRUE;
 171     }
 172     else if (run>=114737&&run<=117223) {      //period="LHC10B";
 173       if (!fIsMC) {
 174         if (fRecoPass == 2) {
 175           fCorrectExpTimes=kTRUE;
 176           fCorrectTRDBug=kFALSE;}
 177         else if (fRecoPass == 3) {
 178           fCorrectExpTimes=kFALSE;
 179           fCorrectTRDBug=kTRUE;
 180         }
 181         fLHC10dPatch=kFALSE;
 182         fT0IntercalibrationShift = 0;
 183         fT0DetectorAdjust=kFALSE;  // previously was true (we acted as a T0 tender)
 184       } else if (fIsMC) {
 185         fT0DetectorAdjust=kFALSE;
 186         fT0Simulate=kTRUE;
 187         fCorrectExpTimes=kTRUE;
 188       }
 189     }
 190     else if (run>=118503&&run<=121040) { //period="LHC10C";
 191       if (!fIsMC) {
 192         if (fRecoPass == 2) {fCorrectExpTimes=kTRUE; fCorrectTRDBug=kFALSE;}
 193         else if (fRecoPass == 3) {fCorrectExpTimes=kFALSE; fCorrectTRDBug=kTRUE;}
 194         fLHC10dPatch=kFALSE;
 195         fT0IntercalibrationShift = 0;
 196         fT0DetectorAdjust=kFALSE;  // previously was true (we acted as a T0 tender)
 197       } else if (fIsMC) {
 198         fT0DetectorAdjust=kFALSE;
 199         fT0Simulate=kTRUE;
 200         fCorrectExpTimes=kTRUE;
 201       }
 202     }
 203     else if (run>=122195&&run<=126437) { //period="LHC10D";
 204       if (!fIsMC) {
 205         fCorrectExpTimes=kFALSE;
 206         fLHC10dPatch=kTRUE;
 207         fT0DetectorAdjust=kFALSE;  // previously was true (we acted as a T0 tender)
 208         fT0IntercalibrationShift = 0;
 209       } else if (fIsMC) {
 210         fCorrectExpTimes=kTRUE;    // for old MC the expected times bug was there
 211         fLHC10dPatch=kFALSE;       // but not the fake geometry
 212         fT0DetectorAdjust=kFALSE;
 213         fT0Simulate=kTRUE;
 214       }
 215     }
 216     else if (run>=127719&&run<=130850) { //period="LHC10E";
 217       if (!fIsMC) {
 218         fCorrectExpTimes=kFALSE;
 219         fLHC10dPatch=kFALSE;
 220         fT0DetectorAdjust=kFALSE;            // previously was true (we acted as a T0 tender)
 221         fT0IntercalibrationShift = 0.;       // this was 30 before, but it is now handled via TOFPIDResponse inside OADB
 222       } else if (fIsMC) {
 223         fCorrectExpTimes=kTRUE;    // this is not fully correct for newer productions like LHC11b2 but we live with this
 224         fT0DetectorAdjust=kFALSE;
 225         fT0Simulate=kTRUE;
 226       }
 227     }
 228     else if (run>=133004&&run<=135029) { //period="LHC10F";
 229       fTenderNoAction=kTRUE;
 230     }
 231     else if (run>=135654&&run<=136377) { //period="LHC10G";
 232       fTenderNoAction=kTRUE;
 233     }
 234     else if (run>=136851&&run<=139517) { //period="LHC10H" - pass2;
 235       fCorrectExpTimes=kFALSE;
 236       fLHC10dPatch=kFALSE;
 237       fT0IntercalibrationShift = 0.;
 238       fT0DetectorAdjust=kFALSE;          // it was kTRUE
 239     }
 240     else if (run>=139699) {              //period="LHC11A";
 241       fTenderNoAction=kTRUE;
 242     }
 243   }
 244
 245   if (fTenderNoAction) {
 246     AliInfo(" |---------------------------------------------------------------------------|");
 247     AliInfo(" |                                                                           |");
 248     AliInfo(Form(" |  TOF tender is not supported for run %d                           |",run));
 249     AliInfo(" |  TOF tender will do nothing.                                              |");
 250     AliInfo(" |  Check TOF tender usage for run/periods at:                               |");
 251     AliInfo(" |  https://twiki.cern.ch/twiki/bin/view/ALICE/TOF.                          |");
 252     AliInfo(" |---------------------------------------------------------------------------|");
 253     AliInfo(" ");
 254     return;
 255   }
 256
 257
 258   // Load from OADB TOF resolution
 259   LoadTOFPIDParams(run);
 260
 261   // Check if another tender wagon already created the esd pid object
 262   // if not we create it and set it to the ESD input handler
 263   fESDpid=fTender->GetESDhandler()->GetESDpid();
 264   if (!fESDpid) {
 265     fESDpid=new AliESDpid;
 266     fTender->GetESDhandler()->SetESDpid(fESDpid);
 267   }
 268
 269
 270   // Configure TOF calibration class
 271   if (!fTOFCalib)fTOFCalib=new AliTOFcalib();  // create if needed
 272   fTOFCalib->SetRemoveMeanT0(!(fIsMC));        // must be kFALSE on MC (default is kTRUE)
 273   fTOFCalib->SetCalibrateTOFsignal(!(fIsMC));  // must be kFALSE on MC (no new calibration) (default is kTRUE)
 274   fTOFCalib->SetCorrectTExp(fCorrectExpTimes); // apply a fine tuning on the expected times at low momenta
 275                                                // (this is done for LHC10b/LHC10c pass2)
 276
 277   // Configure TOFT0 maker class
 278   //  if (!fTOFT0maker) fTOFT0maker = new AliTOFT0maker(fESDpid,fTOFCalib); // create if needed
 279   if (!fTOFT0maker) fTOFT0maker = new AliTOFT0maker(fESDpid); // without passing AliTOFCalib it uses the diamond
 280   fTOFT0maker->SetTimeResolution(fTOFPIDParams->GetTOFresolution());     // set TOF resolution for the PID
 281   fTOFT0maker->SetTOFT0algorithm(2);
 282
 283   AliInfo("|******************************************************|");
 284   AliInfo(Form("|    Alice TOF Tender Initialisation (Run %d)  |",fTender->GetRun()));
 285   AliInfo("|    Settings:                                         |");
 286   AliInfo(Form("|    Correct Exp Times              :  %d               |",fCorrectExpTimes));
 287   AliInfo(Form("|    Correct TRD Bug                :  %d               |",fCorrectTRDBug));
 288   AliInfo(Form("|    LHC10d patch                   :  %d               |",fLHC10dPatch));
 289   AliInfo(Form("|    TOF resolution for TOFT0 maker :  %5.2f (ps)     |",fTOFPIDParams->GetTOFresolution()));
 290   AliInfo(Form("|    MC flag (start time added)     :  %d               |",fIsMC));
 291   AliInfo(Form("|    T0 detector offsets applied    :  %d               |",fT0DetectorAdjust));
 292   AliInfo(Form("|    T0 signal re-sampled           :  %d               |",fT0Simulate));
 293   AliInfo(Form("|    TOF/T0 intercalibration shift   :  %5.2f (ps)     |",fT0IntercalibrationShift));
 294   AliInfo("|******************************************************|");
 295
 296
 297 }
 298
 299 //_____________________________________________________
 300 void AliTOFTenderSupply::ProcessEvent()
 301 {
 302   //
 303   // Use updated calibrations for TOF and T0, reapply PID information
 304   // For MC: timeZero sampling and additional smearing for T0
 305
 306   if (fDebugLevel > 1) AliInfo("process event");
 307
 308   AliESDEvent *event=fTender->GetEvent();
 309   if (!event) return;
 310   if (fDebugLevel > 1) AliInfo("event read");
 311
 312
 313
 314   if (fTender->RunChanged()){
 315
 316     Init();
 317
 318     if (fTenderNoAction) return;
 319     Int_t versionNumber = GetOCDBVersion(fTender->GetRun());
 320     fTOFCalib->SetRunParamsSpecificVersion(versionNumber);
 321     fTOFCalib->Init(fTender->GetRun());
 322
 323     if(event->GetT0TOF()){ // read T0 detector correction from OCDB
 324         // OCDB instance
 325         if (fT0DetectorAdjust) {
 326           AliCDBManager* ocdbMan = AliCDBManager::Instance();
 327           ocdbMan->SetRun(fTender->GetRun());
 328           AliCDBEntry *entry = ocdbMan->Get("T0/Calib/TimeAdjust/");
 329           if(entry) {
 330             AliT0CalibSeasonTimeShift *clb = (AliT0CalibSeasonTimeShift*) entry->GetObject();
 331             Float_t *t0means= clb->GetT0Means();
 332             //      Float_t *t0sigmas = clb->GetT0Sigmas();
 333             fT0shift[0] = t0means[0] + fT0IntercalibrationShift;
 334             fT0shift[1] = t0means[1] + fT0IntercalibrationShift;
 335             fT0shift[2] = t0means[2] + fT0IntercalibrationShift;
 336             fT0shift[3] = t0means[3] + fT0IntercalibrationShift;
 337           } else {
 338             for (Int_t i=0;i<4;i++) fT0shift[i]=0;
 339             AliWarning("TofTender no T0 entry found T0shift set to 0");
 340           }
 341         } else {
 342           for (Int_t i=0;i<4;i++) fT0shift[i]=0;
 343         }
 344     }
 345   }
 346
 347   if (fTenderNoAction) return;
 348
 349   fTOFCalib->CalibrateESD(event);   //recalculate TOF signal (no harm for MC, see settings inside init)
 350
 351
 352   // patches for various reconstruction bugs
 353   if (fLHC10dPatch && !(fIsMC)) RecomputeTExp(event);   // LHC10d pass2: fake full TRD geometry
 354   if ( (fCorrectTRDBug && !(fIsMC)) || (fForceCorrectTRDBug)) FixTRDBug(event);     // LHC10b,c pass3: wrong TRD dE/dx
 355
 356   Double_t startTime = 0.;
 357   if (fIsMC) startTime = fTOFCalib->TuneForMC(event,fTOFPIDParams->GetTOFresolution());   // this is for old MC when we didn't jitter startTime in MC
 358
 359   if (fDebugLevel > 1) Printf(" TofTender: startTime %f",startTime);
 360   if (fDebugLevel > 1) Printf(" TofTender: T0 time (orig) %f %f %f",event->GetT0TOF(0),event->GetT0TOF(1),event->GetT0TOF(2));
 361
 362   // event by event TO detector treatment
 363   if(event->GetT0TOF()){   // protection: we adjust T0 only if it is there....
 364
 365     if (event->GetT0TOF(0) == 0) event->SetT0TOF(0, 9999999.); // in case no information we set to unknown
 366     if (event->GetT0TOF(1) == 0) event->SetT0TOF(1, 99999.);
 367     if (event->GetT0TOF(2) == 0) event->SetT0TOF(2, 99999.);
 368
 369     if ( (fT0DetectorAdjust) && !(fIsMC) ) {  // DATA: apply shifts to align around T0: this is like a T0 tender!!
 370         event->SetT0TOF(0,event->GetT0TOF(0) - fT0shift[0]);
 371         event->SetT0TOF(1,event->GetT0TOF(1) - fT0shift[1]);
 372         event->SetT0TOF(2,event->GetT0TOF(2) - fT0shift[2]);
 373     }
 374     if (fIsMC) {
 375       if (fT0DetectorAdjust)  { // MC case 1: add an additional contribution to resolution
 376         // MC: add smearing for realistic T0A and T0C resolution
 377         Double_t defResolutionT0A = 33.;   // in future we will get this from ESDrun data structure or via OCDB
 378         Double_t defResolutionT0C = 30.;   // for the moment we don't trust them
 379         if ( (fgT0Aresolution > defResolutionT0A) && (event->GetT0TOF(1)<90000.) ) { // add smearing only if signal is there
 380           Double_t addedSmearingT0A = TMath::Sqrt(fgT0Aresolution*fgT0Aresolution - defResolutionT0A*defResolutionT0A);
 381           Double_t smearingT0A = gRandom->Gaus(0.,addedSmearingT0A);
 382           event->SetT0TOF(1,event->GetT0TOF(1) + smearingT0A);
 383         }
 384         if ( (fgT0Cresolution > defResolutionT0C) && (event->GetT0TOF(2)<90000.) ) { // add smearing only if signal is there
 385         Double_t addedSmearingT0C = TMath::Sqrt(fgT0Cresolution*fgT0Cresolution - defResolutionT0C*defResolutionT0C);
 386         Double_t smearingT0C = gRandom->Gaus(0.,addedSmearingT0C);
 387         event->SetT0TOF(2,event->GetT0TOF(2) + smearingT0C);
 388         }
 389         if (event->GetT0TOF(0)<90000.) { // we recompute the AND only if it is already there...
 390           Double_t smearedT0AC = (event->GetT0TOF(1)+event->GetT0TOF(2))/2.;
 391           event->SetT0TOF(0,smearedT0AC);
 392         }
 393         if (fDebugLevel > 1) Printf(" TofTender: T0 time (postSmear) %f %f %f",event->GetT0TOF(0),event->GetT0TOF(1),event->GetT0TOF(2));
 394       }
 395       else if (fT0Simulate) {  // MC case 2: we completely simulate signal in T0 based on multiplicity in ITS and vtx position
 396         event->SetT0TOF(0, 9999999.); // we wipe-out whatever is there
 397         event->SetT0TOF(1, 99999.);
 398         event->SetT0TOF(2, 99999.);
 399         if (fDebugLevel > 1) Printf(" TofTender: T0 time (after wipe-out) %f %f %f",event->GetT0TOF(0),event->GetT0TOF(1),event->GetT0TOF(2));
 400         AliESDVertex *fvtx = (AliESDVertex*)event->GetPrimaryVertex();
 401         Double_t zvtx = fvtx->GetZ();
 402         Double_t tracklets[2] = {0.,0.};
 403         GetTrackletsForT0(event,&tracklets[0],&tracklets[1]);
 404         if (fDebugLevel > 1) Printf(" TofTender: T0 simul (z vtx tracklets A/C) %f %f %f",zvtx,tracklets[0],tracklets[1]);
 405         for (Int_t side = 0; side < 2; side ++) {    // side 0 = T0A - side 1 = T0C
 406           Double_t signal = SampleT0Signal(side,zvtx,tracklets[side]);  // if not fired we return 99999.
 407           event->SetT0TOF(side+1,signal);            // but for the T0 structure we need to add 1...
 408         }
 409         if ( (event->GetT0TOF(1) < 1000.) && (event->GetT0TOF(2) < 1000.) ) { // both signals are there
 410           Double_t meanT0AC=(event->GetT0TOF(1)+event->GetT0TOF(2))/2.;
 411           event->SetT0TOF(0,meanT0AC);
 412         }
 413         if (fDebugLevel > 1) Printf(" TofTender: T0 simul (AC A C) %f %f %f",event->GetT0TOF(0),event->GetT0TOF(1),event->GetT0TOF(2));
 414       }
 415       // add the startTime offset also to the T0 detector information
 416       event->SetT0TOF(0,event->GetT0TOF(0) + startTime);
 417       event->SetT0TOF(1,event->GetT0TOF(1) + startTime);
 418       event->SetT0TOF(2,event->GetT0TOF(2) + startTime);
 419       if (fDebugLevel > 1) Printf(" TofTender: T0 time (postStart AC A C) %f %f %f",event->GetT0TOF(0),event->GetT0TOF(1),event->GetT0TOF(2));
 420     }
 421     // after shifts adjust (data) or smearing+offset (MC) we 'clean' to default if signals not there
 422     if(event->GetT0TOF(0) > 900000) event->SetT0TOF(0, 999999.);
 423     if(event->GetT0TOF(1) > 90000)  event->SetT0TOF(1, 99999.);
 424     if(event->GetT0TOF(2) > 90000)  event->SetT0TOF(2, 99999.);
 425   }
 426   if (fDebugLevel > 1) Printf(" TofTender: T0 time (FINAL) %f %f %f",event->GetT0TOF(0),event->GetT0TOF(1),event->GetT0TOF(2));
 427
 428   //compute timeZero of the event via TOF-TO
 429   fTOFT0maker->ComputeT0TOF(event);
 430   fTOFT0maker->WriteInESD(event);
 431
 432   //  set preferred startTime: this is now done via AliPIDResponseTask
 433   fESDpid->SetTOFResponse(event, (AliESDpid::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
 434
 435   // recalculate PID probabilities
 436   // this is for safety, especially if the user doesn't attach a PID tender after TOF tender
 437   Int_t ntracks=event->GetNumberOfTracks();
 438   //  AliESDtrack *track = NULL;
 439   //  Float_t tzeroTrack = 0;
 440   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 441     //    track = event->GetTrack(itrack);
 442     //    tzeroTrack = fESDpid->GetTOFResponse().GetStartTime(track->P());
 443     //    Printf("================> Track # %d mom: %f tzeroTrack %f",itrack,track->P(),tzeroTrack);
 444     fESDpid->MakeTOFPID(event->GetTrack(itrack),0);
 445   }
 446
 447
 448 }
 449
 450
 451 //_____________________________________________________
 452 void AliTOFTenderSupply::RecomputeTExp(AliESDEvent *event) const
 453 {
 454   /*
 455    * calibrate TExp
 456    */
 457
 458
 459   /* loop over tracks */
 460   AliESDtrack *track = NULL;
 461   for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) {
 462     /* get track and calibrate */
 463     track = event->GetTrack(itrk);
 464     RecomputeTExp(track);
 465   }
 466
 467 }
 468
 469 //_____________________________________________________
 470 void AliTOFTenderSupply::RecomputeTExp(AliESDtrack *track) const
 471 {
 472   /***
 473        THIS METHOD IS BASED ON THEORETICAL EXPECTED TIME COMPUTED
 474        USING AVERAGE MOMENTUM BETWEEN INNER/OUTER TRACK PARAMS
 475        IT IS A ROUGH APPROXIMATION APPLIED TO FIX LHC10d-pass2 DATA
 476        WHERE A WRONG GEOMETRY (FULL TRD) WAS INSERTED
 477   ***/
 478
 479   Double_t texp[AliPID::kSPECIESC];
 480   if (!track || !(track->GetStatus() & AliESDtrack::kTOFout)) return;
 481
 482
 483   /* get track params */
 484   Float_t l = track->GetIntegratedLength();
 485   Float_t p = track->P();
 486   if (track->GetInnerParam() && track->GetOuterParam()) {
 487     Float_t pin = track->GetInnerParam()->P();
 488     Float_t pout = track->GetOuterParam()->P();
 489     p = 0.5 * (pin + pout);
 490   }
 491   /* loop over particle types and compute expected time */
 492   for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++)
 493     texp[ipart] = GetExpTimeTh(AliPID::ParticleMass(ipart), p, l) - 37.;
 494   // 37 is a final semiempirical offset to further adjust (calibrations were
 495   // done with "standard" integratedTimes)
 496
 497   // in old ESDs like this horridous LHC10d pass2 light nuclei were not supported...
 498   for (Int_t ipart = AliPID::kProton+1; ipart < AliPID::kSPECIESC; ipart++) texp[ipart]=0.;
 499
 500   /* set integrated times */
 501   track->SetIntegratedTimes(texp);
 502
 503 }
 504
 505
 506 //______________________________________________________________________________
 507 void AliTOFTenderSupply::DetectRecoPass()
 508 {
 509   //
 510   // Detect reconstruction information
 511   //
 512
 513   //reset information
 514   fRecoPass=0;
 515
 516   //Get the current file to check the reconstruction pass (UGLY, but not stored in ESD... )
 517   AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
 518   AliVEventHandler *inputHandler=mgr->GetInputEventHandler();
 519   if (!inputHandler) return;
 520
 521   TTree *tree= (TTree*)inputHandler->GetTree();
 522   TFile *file= (TFile*)tree->GetCurrentFile();
 523
 524   if (!file) {
 525     AliFatal("Current file not found");
 526     return; // coverity
 527   }
 528
 529   //find pass from file name (UGLY, but not stored in ESD... )
 530   TString fileName(file->GetName());
 531   if (fileName.Contains("pass1") ) {
 532     fRecoPass=1;
 533   } else if (fileName.Contains("pass2") ) {
 534     fRecoPass=2;
 535   } else if (fileName.Contains("pass3") ) {
 536     fRecoPass=3;
 537   } else if (fileName.Contains("pass4") ) {
 538     fRecoPass=4;
 539   } else if (fileName.Contains("pass5") ) {
 540     fRecoPass=5;
 541   } else if (fileName.Contains("pass6") ) {
 542     fRecoPass=6;
 543   }
 544   if (fRecoPass == 0) {
 545     AliInfo(Form("From file name %s reco pass cannot be detected",fileName.Data()));
 546     AliInfo("Change file name or use SetUserRecoPass method");
 547     AliFatal("------------- TOF tender cannot run with reco pass unspecified, issuing FATAL error ---------- ");
 548   }
 549 }
 550
 551
 552 //______________________________________________________________________________
 553 void AliTOFTenderSupply::InitGeom()
 554 {
 555
 556   if (fGeomSet == kTRUE) return;
 557
 558   //  Printf("\n \n ----- calling InitGeom to fix TRD Bug ----- \n \n");
 559   AliCDBManager * man = AliCDBManager::Instance();
 560   man->SetDefaultStorage("raw://");
 561   fCDBkey = man->SetLock(kFALSE, fCDBkey);
 562   Int_t run = fTender->GetRun();
 563   //  Printf(" ---------> run is %d",run);
 564   man->SetRun(run);
 565   AliCDBEntry *entry = (AliCDBEntry*)man->Get("GRP/Geometry/Data");
 566   if (entry) {
 567     AliGeomManager::LoadGeometry();
 568     AliGeomManager::ApplyAlignObjsFromCDB("ITS TPC TRD TOF");
 569     //    fCDBkey = man->SetLock(kTRUE, fCDBkey);
 570     //    Printf("\n \n ----- Geometry loaded ------ \n \n");
 571   }
 572   fGeomSet=kTRUE;
 573
 574 }
 575
 576
 577 //______________________________________________________________________________
 578 void AliTOFTenderSupply::FixTRDBug(AliESDEvent* event)
 579 //
 580 // recompute texp fixing wrong dE/dx from TRD (LHC10b,c pass3)
 581 //
 582 {
 583
 584   if (fGeomSet == kFALSE) InitGeom();
 585
 586   //  Printf("Running FixTRD bug ");
 587   /* loop over tracks */
 588   AliESDtrack *track = NULL;
 589   for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) {
 590     track = event->GetTrack(itrk);
 591     FixTRDBug(track);
 592   }
 593 }
 594
 595
 596 //_____________________________________________________
 597 void AliTOFTenderSupply::FixTRDBug(AliESDtrack *track)
 598 {
 599   //
 600   //
 601   //
 602
 603
 604     ULong_t status=track->GetStatus();
 605     if (!( ( (status & AliVTrack::kITSrefit)==AliVTrack::kITSrefit ) &&
 606            ( (status & AliVTrack::kTPCrefit)==AliVTrack::kTPCrefit ) &&
 607            ( (status & AliVTrack::kTPCout)==AliVTrack::kTPCout ) &&
 608            ( (status & AliVTrack::kTOFout)==AliVTrack::kTOFout ) &&
 609            ( (status & AliVTrack::kTIME)==AliVTrack::kTIME ) ) ) return;
 610
 611     fIsEnteringInTRD=kFALSE;
 612     fInTRD=kFALSE;
 613     fIsComingOutTRD=kFALSE;
 614     fOutTRD=kFALSE;
 615
 616     //    Printf("Track reached TOF %f",track->P());
 617     Double_t correctionTimes[AliPID::kSPECIES] = {0.,0.,0.,0.,0.}; // to be added to the expected times
 618     FindTRDFix(track, correctionTimes);
 619     Double_t expectedTimes[AliPID::kSPECIESC] = {0.,0.,0.,0.,0.,0.,0.,0.,0.};
 620     track->GetIntegratedTimes(expectedTimes,AliPID::kSPECIESC);
 621     //    Printf("Exp. times: %f %f %f %f %f",
 622     //     expectedTimes[0],expectedTimes[1],expectedTimes[2],expectedTimes[3],expectedTimes[4]);
 623     //    Printf("Corr. times: %f %f %f %f %f",
 624     //     correctionTimes[0],correctionTimes[1],correctionTimes[2],correctionTimes[3],correctionTimes[4]);
 625
 626     for (Int_t jj=0; jj<AliPID::kSPECIES; jj++) expectedTimes[jj]+=correctionTimes[jj];
 627     // we only correct up to protons, for this horridous LHC10d pass2 light nuclei were not supported
 628     // so expected times will remain zero
 629     track->SetIntegratedTimes(expectedTimes);
 630
 631 }
 632
 633
 634 //________________________________________________________________________
 635 void AliTOFTenderSupply::FindTRDFix(AliESDtrack *track,Double_t *corrections)
 636 {
 637
 638   Double_t pT = track->Pt();
 639   ULong_t status=track->GetStatus();
 640   Bool_t isTRDout = (status & AliVTrack::kTRDout)==AliVTrack::kTRDout;
 641
 642   Double_t length = 0.;
 643
 644   Double_t xyzIN[3]={0.,0.,0.};
 645   fIsEnteringInTRD = track->GetXYZAt(fRhoTRDin,fMagField,xyzIN);
 646
 647   Double_t xyzOUT[3]={0.,0.,0.};
 648   fIsComingOutTRD = track->GetXYZAt(fRhoTRDout,fMagField,xyzOUT);
 649
 650   if (fIsEnteringInTRD && fIsComingOutTRD) {
 651
 652
 653     Double_t phiIN = TMath::Pi()+TMath::ATan2(-xyzIN[1],-xyzIN[0]);
 654     phiIN *= TMath::RadToDeg();
 655     fInTRD = ( (phiIN>=  0. && phiIN<= 40.) ||
 656                (phiIN>=140. && phiIN<=220.) ||
 657                (phiIN>=340. && phiIN<=360.) ); // TRD SMs installed @ 2010
 658
 659     Double_t phiOUT = TMath::Pi()+TMath::ATan2(-xyzOUT[1],-xyzOUT[0]);
 660     phiOUT *= TMath::RadToDeg();
 661     fOutTRD = ( (phiOUT>=  0. && phiOUT<= 40.) ||
 662                 (phiOUT>=140. && phiOUT<=220.) ||
 663                 (phiOUT>=340. && phiOUT<=360.) ); // TRD SMs installed @ 2010
 664
 665     length = 0.;
 666
 667     if (fInTRD || fOutTRD) {
 668
 669       if ( ( fInTRD && fOutTRD ) || ( fInTRD && !fOutTRD ) ) {
 670         length = EstimateLengthInTRD1(track);
 671       } else if ( !fInTRD && fOutTRD ) {
 672         length = EstimateLengthInTRD2(track);
 673       }
 674
 675     } else { // ( !fInTRD && !fOutTRD )
 676
 677       length = EstimateLengthOutTRD(track);
 678
 679     }
 680
 681   }
 682   //  Printf("estimated length in TRD %f [isTRDout %d]",length,isTRDout);
 683   CorrectDeltaTimes(pT,length,isTRDout,corrections);
 684
 685 }
 686
 687 //________________________________________________________________________
 688 void AliTOFTenderSupply::CorrectDeltaTimes(Double_t pT,
 689                                                 Double_t length,
 690                                                 Bool_t flagTRDout,
 691                                                 Double_t *corrections)
 692 {
 693
 694   corrections[2] = CorrectExpectedPionTime(pT,length,flagTRDout);
 695   corrections[0] = corrections[2]; // x electrons used pion corrections
 696   corrections[1] = corrections[2]; // x muons used pion corrections
 697   corrections[3] = CorrectExpectedKaonTime(pT,length,flagTRDout);
 698   corrections[4] = CorrectExpectedProtonTime(pT,length,flagTRDout);
 699
 700 }
 701
 702 //________________________________________________________________________
 703 Double_t AliTOFTenderSupply::CorrectExpectedPionTime(Double_t pT,
 704                                                           Double_t length,
 705                                                           Bool_t isTRDout)
 706 {
 707   // correction for expected time for pions
 708
 709   Double_t delta=0.;
 710
 711   //  Printf("Flags Ent In ComingOut Out %d %d %d %d",fIsEnteringInTRD,fInTRD,fIsComingOutTRD,fOutTRD);
 712   if (!fIsEnteringInTRD || !fIsComingOutTRD) { // zone 5
 713
 714     Float_t p[2]={0.,0.};
 715
 716     if (isTRDout) {
 717
 718       if (pT<0.30) {
 719         p[0] = 180.; p[1] = 0.;
 720       } else if (pT>=0.30 && pT<0.35) {
 721         p[0] = 740.; p[1] = -1800.;
 722       } else if (pT>=0.35 && pT<0.40) {
 723         p[0] = 488.; p[1] =-1080.;
 724       } else if (pT>=0.40 && pT<0.45) {
 725         p[0] = 179.; p[1] = -307.;
 726       } else if (pT>=0.45 && pT<0.50) {
 727         p[0] =  97.; p[1] = -123.;
 728       } else { //if (pT>=0.50)
 729         p[0] = 120.; p[1] = -172.;
 730       }
 731
 732     } else {
 733
 734       if (pT<0.30) {
 735         p[0] =  70.; p[1] =    0.;
 736       } else if (pT>=0.30 && pT<0.35) {
 737         p[0] = 339.; p[1] = -927.;
 738       } else if (pT>=0.35 && pT<0.40) {
 739         p[0] =  59.; p[1] = -121.;
 740       } else if (pT>=0.40 && pT<0.50) {
 741         p[0] =  21.; p[1] =  -24.;
 742       } else { //if (pT>=0.50)
 743         p[0] =  42.; p[1] =  -67.;
 744       }
 745
 746     }
 747
 748     delta = p[0]+p[1]*pT;
 749     //    Printf("Pion time: %f %f %f %f",p[0],p[1],length,delta);
 750
 751   } else {
 752
 753     Float_t p[2] = {0.,0.};
 754
 755     if ( fInTRD &&  fOutTRD) { // zone 1
 756
 757       if (isTRDout) {
 758
 759         if (length<130.) {
 760           p[0] = 0.; p[1] = 0.;
 761         } else if (length>=130. && length<170.) {
 762           p[0] = -20.5; p[1] = 0.25;
 763         } else {//if (length>=170.)
 764           p[0] = 22.; p[1] = 0.;
 765         }
 766
 767       } else { // !isTRDout
 768
 769         p[0] = 20.; p[1] = 0.;
 770
 771       }
 772
 773     } else if (!fInTRD && !fOutTRD) { // zone 2
 774
 775       p[0] = 0.; p[1] = 0.;
 776
 777     } else if ( fInTRD &&  !fOutTRD) { // zone 3
 778
 779       if (isTRDout) {
 780
 781         if (length< 75.) {
 782           p[0] = 17.; p[1] =  0.;
 783         } else if (length>= 75. && length< 95.) {
 784           p[0] = 81.; p[1] = -0.85;
 785         } else if (length>= 95. && length<155.) {
 786           p[0] =  0.; p[1] =  0.;
 787         } else {//if (length>=155.)
 788           p[0] = 10.; p[1] =  0.;
 789         }
 790
 791       } else { // !isTRDout
 792
 793         p[0] = 0.; p[1] = 0.;
 794
 795       }
 796
 797     } else if (!fInTRD &&  fOutTRD) { // zone 4
 798
 799       if (isTRDout) {
 800
 801         if (length<80.) {
 802           p[0] =  0.; p[1] =  0.;
 803         } else {//if (length>=80.)
 804           p[0] = 10.; p[1] =  0.;
 805         }
 806
 807       } else { // !isTRDout
 808
 809         if (length<30.) {
 810           p[0] =  0.; p[1] =  0.;
 811         } else {//if (length>=30.)
 812           p[0] =  6.; p[1] =  0.;
 813         }
 814
 815       }
 816
 817     }
 818
 819     delta = p[0]+p[1]*length;
 820     //    Printf("Pion time: %f %f %f %f",p[0],p[1],length,delta);
 821
 822   }
 823
 824   return delta;
 825
 826 }
 827
 828 //________________________________________________________________________
 829 Double_t AliTOFTenderSupply::CorrectExpectedKaonTime(Double_t pT,
 830                                                           Double_t length,
 831                                                           Bool_t isTRDout)
 832 {
 833   // correction for expected time for kaons
 834
 835   Double_t delta=0.;
 836   //  Printf("Flags Ent In ComingOut Out %d %d %d %d",fIsEnteringInTRD,fInTRD,fIsComingOutTRD,fOutTRD);
 837
 838   if (!fIsEnteringInTRD || !fIsComingOutTRD) { // zone 5
 839
 840     Float_t p[2]={0.,0.};
 841
 842     if (isTRDout) {
 843
 844       if (pT<0.4) {
 845         p[0] =  900.; p[1] = 0.;
 846       } else if (pT>=0.40 && pT<0.45) {
 847         p[0] = 3100.; p[1] = -6000.;
 848       } else if (pT>=0.45 && pT<0.50) {
 849         p[0] = 1660.; p[1] = -2800.;
 850       } else if (pT>=0.50 && pT<0.55) {
 851         p[0] =  860.; p[1] = -1200.;
 852       } else { //if (pT>=0.55)
 853         p[0] =  200.; p[1] = 0.;
 854       }
 855
 856     } else {
 857
 858       if (pT<0.30) {
 859         p[0] =   0.; p[1] =    0.;
 860       } else if (pT>=0.30 && pT<0.32) {
 861         p[0] = 570.; p[1] =    0.;
 862       } else if (pT>=0.32 && pT<0.35) {
 863         p[0] = 3171.; p[1] = -8133.;
 864       } else if (pT>=0.35 && pT<0.40) {
 865         p[0] = 1815.; p[1] = -4260.;
 866       } else if (pT>=0.40 && pT<0.45) {
 867         p[0] =  715.; p[1] =  -1471.;
 868       } else if (pT>=0.45 && pT<0.50) {
 869         p[0] =  233.; p[1] =  -407.;
 870       } else if (pT>=0.50 && pT<0.55) {
 871         p[0] =  408.; p[1] =  -752.;
 872       } else { //if (pT>=0.55)
 873         p[0] =  408.-752.*0.55; p[1] = 0.;
 874       }
 875
 876     }
 877
 878     delta = p[0]+p[1]*pT;
 879     //    Printf("Kaon time: %f %f %f %f",p[0],p[1],length,delta);
 880
 881   } else {
 882
 883     Float_t p[2] = {0.,0.};
 884
 885     if ( fInTRD &&  fOutTRD) { // zone 1
 886
 887       if (isTRDout) {
 888
 889         if (length<95.) {
 890           p[0] = 20.; p[1] = 0.;
 891         } else if (length>=95. && length<195.) {
 892           p[0] = -24.0; p[1] = 0.10+0.0041*length;
 893         } else {//if (length>=195.)
 894           p[0] =  150.; p[1] = 0.;
 895         }
 896
 897       } else { // !isTRDout
 898
 899         p[0] = 40.; p[1] = 0.;
 900
 901       }
 902
 903     } else if (!fInTRD && !fOutTRD) { // zone 2
 904
 905       p[0] = 0.; p[1] = 0.;
 906
 907     } else if ( fInTRD &&  !fOutTRD) { // zone 3
 908
 909       if (isTRDout) {
 910
 911         if (length< 15.) {
 912           p[0] = 180.; p[1] =  0.;
 913         } else if (length>= 15. && length< 55.) {
 914           p[0] = 215.; p[1] = -2.5;
 915         } else {//if (length>=55.)
 916           p[0] = 78.; p[1] =  0.;
 917         }
 918
 919       } else { // !isTRDout
 920
 921         p[0] = 0.; p[1] = 0.;
 922
 923       }
 924
 925     } else if (!fInTRD &&  fOutTRD) { // zone 4
 926
 927       if (isTRDout) {
 928
 929         if (length< 55.) {
 930           p[0] =  0.; p[1] =  0.;
 931         } else if (length>= 55. && length<115.) {
 932           p[0] = -85.; p[1] = 1.9;
 933         } else {//if (length>=115.)
 934           p[0] = 100.; p[1] =  0.;
 935         }
 936
 937       } else { // !isTRDout
 938
 939         p[0] =  0.; p[1] =  0.;
 940
 941       }
 942
 943     }
 944
 945     delta = p[0]+p[1]*length;
 946     //    Printf("Kaon time: %f %f %f %f",p[0],p[1],length,delta);
 947
 948   }
 949
 950   return delta;
 951
 952 }
 953
 954 //________________________________________________________________________
 955 Double_t AliTOFTenderSupply::CorrectExpectedProtonTime(Double_t pT,
 956                                                             Double_t length,
 957                                                             Bool_t isTRDout)
 958 {
 959   // correction for expected time for protons
 960
 961   Double_t delta=0.;
 962   //  Printf("Flags Ent In ComingOut Out %d %d %d %d",fIsEnteringInTRD,fInTRD,fIsComingOutTRD,fOutTRD);
 963
 964   if (!fIsEnteringInTRD || !fIsComingOutTRD) { // zone 5
 965     Float_t p[2]={0.,0.};
 966
 967
 968     if (isTRDout) {
 969
 970       if (pT<0.375) {
 971         p[0] = 1000.; p[1] = 0.;
 972       } else if (pT>=0.375 && pT<0.45) {
 973         p[0] = 1500.; p[1] = 0.;
 974       } else if (pT>=0.45 && pT<0.50) {
 975         p[0] = 4650.; p[1] = -7000.;
 976       } else if (pT>=0.50 && pT<0.55) {
 977         p[0] = 3150.; p[1] = -4000.;
 978       } else { //if (pT>=0.55)
 979         p[0] = 3150. -4000.*0.55; p[1] = 0.;
 980       }
 981
 982     } else {
 983
 984       if (pT<0.32) {
 985         p[0] = 2963.-5670.*0.032; p[1] = 0.;
 986       } else if (pT>=0.32 && pT<0.35) {
 987         p[0] = 2963.; p[1] =  -5670.;
 988       } else if (pT>=0.35 && pT<0.40) {
 989         p[0] = 4270.; p[1] =  -9400.;
 990       } else if (pT>=0.40 && pT<0.45) {
 991         p[0] = 1550.; p[1] =  -2600.;
 992       } else if (pT>=0.45 && pT<0.50) {
 993         p[0] = 1946.; p[1] =  -3480.;
 994       } else if (pT>=0.50 && pT<0.55) {
 995         p[0] = 1193.; p[1] =  -1974.;
 996       } else { //if (pT>=0.55)
 997         p[0] = 1193.-1974.*0.55; p[1] = 0.;
 998       }
 999
1000     }
1001
1002     delta = p[0]+p[1]*pT;
1003     //    Printf("Proton time: %f %f %f %f",p[0],p[1],length,delta);
1004
1005   } else {
1006
1007     Float_t p[2] = {0.,0.};
1008
1009     if ( fInTRD &&  fOutTRD) { // zone 1
1010
1011       if (isTRDout) {
1012
1013         if (length<90.) {
1014           p[0] = 0.; p[1] = 0.;
1015         } else if (length>=90. && length<200.) {
1016           p[0] = 1063.; p[1] = -32.+0.30*length-0.00072*length*length;
1017         } else {//if (length>=200.)
1018           p[0] =  900.; p[1] = 0.;
1019         }
1020
1021       } else { // !isTRDout
1022
1023         p[0] = 80.; p[1] = 0.;
1024
1025       }
1026
1027     } else if (!fInTRD && !fOutTRD) { // zone 2
1028
1029       if (isTRDout) {
1030         p[0] = 0.; p[1] = 0.;
1031       } else {
1032         if (length<125.) {
1033           p[0] = 0.; p[1] = 0.;
1034         } else if (length>=125. && length<180.) {
1035           p[0] = -132.; p[1] = 1.3;
1036         } else {
1037           p[0] = 100.; p[1] = 0.;
1038         }
1039
1040       }
1041
1042     } else if ( fInTRD &&  !fOutTRD) { // zone 3
1043
1044       if (isTRDout) {
1045
1046         if (length< 30.) {
1047           p[0] = 670.; p[1] =  0.;
1048         } else if (length>= 30. && length<155.) {
1049           p[0] = 944.; p[1] = -11.+0.064*length;
1050         } else {//if (length>=155.)
1051           p[0] = 780.; p[1] =  0.;
1052         }
1053
1054       } else { // !isTRDout
1055
1056         if (length< 30.) {
1057           p[0] = 140.; p[1] = -4.5;
1058         } else {
1059           p[0] = 0.; p[1] = 0.;
1060         }
1061
1062       }
1063
1064     } else if (!fInTRD &&  fOutTRD) { // zone 4
1065
1066       if (isTRDout) {
1067
1068         if (length< 45.) {
1069           p[0] = 130.; p[1] =  0.;
1070         } else if (length>= 45. && length<120.) {
1071           p[0] = -190.; p[1] = 6.5;
1072         } else {//if (length>=120.)
1073           p[0] = 750.; p[1] =  0.;
1074         }
1075
1076       } else { // !isTRDout
1077
1078         if (length<75.5) {
1079           p[0] =    0.; p[1] =  0.;
1080         } else if (length>= 75.5 && length<90.) {
1081           p[0] = -830.; p[1] = 11.;
1082         } else {
1083           p[0] =  160.; p[1] =  0.;
1084         }
1085
1086       }
1087
1088     }
1089
1090     delta = p[0]+p[1]*length;
1091     //    Printf("Proton time: %f %f %f %f",p[0],p[1],length,delta);
1092
1093   }
1094
1095   return delta;
1096
1097 }
1098
1099 //________________________________________________________________________
1100 Double_t AliTOFTenderSupply::EstimateLengthInTRD1(AliESDtrack *track)
1101 {
1102
1103   Double_t xyz0[3]={0.,0.,0.};
1104   Bool_t stayInTRD = track->GetXYZAt(fRhoTRDin,fMagField,xyz0);
1105
1106   Double_t phi0 = TMath::Pi()+TMath::ATan2(-xyz0[1],-xyz0[0]);
1107   phi0 *= TMath::RadToDeg();
1108   stayInTRD = stayInTRD && ( (phi0>=  0. && phi0<= 40.) ||
1109                              (phi0>=140. && phi0<=220.) ||
1110                              (phi0>=340. && phi0<=360.) );
1111
1112   Double_t trackLengthInTRD = 0.;
1113   Int_t iStep=0;
1114
1115   Double_t b[3];track->GetBxByBz(b);
1116
1117   Double_t xyz1[3]={0.,0.,0.};
1118   Double_t rho = fRhoTRDin;
1119   while (stayInTRD && rho<=fRhoTRDout) {
1120     iStep++;
1121     rho += fStep;
1122
1123     for (Int_t ii=0; ii<3; ii++) xyz1[ii]=0.;
1124     stayInTRD = track->GetXYZAt(rho,fMagField,xyz1);
1125     Double_t phi1 = TMath::Pi()+TMath::ATan2(-xyz1[1],-xyz1[0]);
1126     phi1 *= TMath::RadToDeg();
1127     stayInTRD = stayInTRD && ( (phi1>=  0. && phi1<= 40.) ||
1128                                (phi1>=140. && phi1<=220.) ||
1129                                (phi1>=340. && phi1<=360.) );
1130
1131     Double_t l2  = TMath::Sqrt((xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0]) +
1132                                (xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1]) +
1133                                (xyz1[2]-xyz0[2])*(xyz1[2]-xyz0[2]));
1134     trackLengthInTRD += l2;
1135
1136     for (Int_t ii=0; ii<3; ii++) xyz0[ii]=xyz1[ii];
1137   }
1138
1139   return trackLengthInTRD;
1140
1141 }
1142
1143 //________________________________________________________________________
1144 Double_t AliTOFTenderSupply::EstimateLengthInTRD2(AliESDtrack *track)
1145 {
1146
1147   Double_t xyz0[3]={0.,0.,0.};
1148   Bool_t stayInTRD = track->GetXYZAt(fRhoTRDout,fMagField,xyz0);
1149
1150   Double_t phi0 = TMath::Pi()+TMath::ATan2(-xyz0[1],-xyz0[0]);
1151   phi0 *= TMath::RadToDeg();
1152   stayInTRD = stayInTRD && ( (phi0>=  0. && phi0<= 40.) ||
1153                              (phi0>=140. && phi0<=220.) ||
1154                              (phi0>=340. && phi0<=360.) );
1155
1156   Double_t trackLengthInTRD = 0.;
1157   Int_t iStep=0;
1158
1159   Double_t b[3];track->GetBxByBz(b);
1160
1161   Double_t xyz1[3]={0.,0.,0.};
1162   Double_t rho = fRhoTRDout;
1163   while (stayInTRD && rho>=fRhoTRDin) {
1164     iStep++;
1165     rho -= fStep;
1166
1167     for (Int_t ii=0; ii<3; ii++) xyz1[ii]=0.;
1168     stayInTRD = track->GetXYZAt(rho,fMagField,xyz1);
1169     Double_t phi1 = TMath::Pi()+TMath::ATan2(-xyz1[1],-xyz1[0]);
1170     phi1 *= TMath::RadToDeg();
1171     stayInTRD = stayInTRD && ( (phi1>=  0. && phi1<= 40.) ||
1172                                (phi1>=140. && phi1<=220.) ||
1173                                (phi1>=340. && phi1<=360.) );
1174
1175     Double_t l2  = TMath::Sqrt((xyz0[0]-xyz1[0])*(xyz0[0]-xyz1[0]) +
1176                                (xyz0[1]-xyz1[1])*(xyz0[1]-xyz1[1]) +
1177                                (xyz0[2]-xyz1[2])*(xyz0[2]-xyz1[2]));
1178     trackLengthInTRD += l2;
1179
1180     for (Int_t ii=0; ii<3; ii++) xyz0[ii]=xyz1[ii];
1181   }
1182
1183   return trackLengthInTRD;
1184
1185 }
1186
1187 //________________________________________________________________________
1188 Double_t AliTOFTenderSupply::EstimateLengthOutTRD(AliESDtrack *track)
1189 {
1190
1191   Double_t xyz0[3]={0.,0.,0.};
1192   Bool_t stayInTRD = track->GetXYZAt(fRhoTRDin,fMagField,xyz0);
1193
1194   Double_t phi0 = TMath::Pi()+TMath::ATan2(-xyz0[1],-xyz0[0]);
1195   phi0 *= TMath::RadToDeg();
1196   stayInTRD = stayInTRD && !( (phi0>=  0. && phi0<= 40.) ||
1197                               (phi0>=140. && phi0<=220.) ||
1198                               (phi0>=340. && phi0<=360.) );
1199
1200   Double_t trackLengthInTRD = 0.;
1201   Int_t iStep=0;
1202
1203   Double_t b[3];track->GetBxByBz(b);
1204
1205   Double_t xyz1[3]={0.,0.,0.};
1206   Double_t rho = fRhoTRDin;
1207   while (stayInTRD && rho<=fRhoTRDout) {
1208     iStep++;
1209     rho += fStep;
1210
1211     for (Int_t ii=0; ii<3; ii++) xyz1[ii]=0.;
1212     stayInTRD = track->GetXYZAt(rho,fMagField,xyz1);
1213     Double_t phi1 = TMath::Pi()+TMath::ATan2(-xyz1[1],-xyz1[0]);
1214     phi1 *= TMath::RadToDeg();
1215     stayInTRD = stayInTRD && !( (phi1>=  0. && phi1<= 40.) ||
1216                                 (phi1>=140. && phi1<=220.) ||
1217                                 (phi1>=340. && phi1<=360.) );
1218
1219     Double_t l2  = TMath::Sqrt((xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0]) +
1220                                (xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1]) +
1221                                (xyz1[2]-xyz0[2])*(xyz1[2]-xyz0[2]));
1222     trackLengthInTRD += l2;
1223
1224     for (Int_t ii=0; ii<3; ii++) xyz0[ii]=xyz1[ii];
1225   }
1226
1227   return trackLengthInTRD;
1228
1229 }
1230
1231 //________________________________________________________________________
1232 Int_t AliTOFTenderSupply::GetOCDBVersion(Int_t runNo)
1233 {
1234   Int_t verNo = -1;
1235   if ( (runNo>=118503 && runNo <=121040) )  { // LHC10C
1236     if (fRecoPass == 2) {                     // on pass2
1237       if (runNo >= 119159 && runNo <= 119163) verNo=3;
1238       else if (runNo >= 119837 && runNo <= 119934) verNo=4;
1239       else if (runNo >= 120067 && runNo <= 120244) verNo=4;
1240       else if (runNo >= 120503 && runNo <= 120505) verNo=4;
1241       else if (runNo >= 120616 && runNo <= 120671) verNo=4;
1242       else if (runNo >= 120741 && runNo <= 120829) verNo=4;
1243     }
1244   }
1245   return verNo;
1246 }
1247
1248
1249 //__________________________________________________________________________
1250 void AliTOFTenderSupply::LoadTOFPIDParams(Int_t runNumber)
1251 {
1252   //
1253   // Load the TOF pid params from the OADB
1254   //
1255
1256   if (fTOFPIDParams) delete fTOFPIDParams;
1257   fTOFPIDParams=0x0;
1258
1259   TFile *oadbf = new TFile("$ALICE_ROOT/OADB/COMMON/PID/data/TOFPIDParams.root");
1260   if (oadbf && oadbf->IsOpen()) {
1261     AliInfo("Loading TOF Params from $ALICE_ROOT/OADB/COMMON/PID/data/TOFPIDParams.root");
1262     AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
1263     if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(runNumber,"TOFparams"));
1264     oadbf->Close();
1265     delete oadbc;
1266   }
1267   delete oadbf;
1268
1269   if (!fTOFPIDParams) {
1270     AliError("TOFPIDParams.root not found in $ALICE_ROOT/OADB/COMMON/PID/data !!");
1271     fTOFPIDParams = new AliTOFPIDParams;
1272     fTOFPIDParams->SetTOFresolution(90.);
1273     fTOFPIDParams->SetStartTimeMethod(AliESDpid::kTOF_T0);
1274   }
1275 }
1276
1277
1278 //__________________________________________________________________________
1279 Double_t AliTOFTenderSupply::SampleT0Signal(Int_t side, Double_t zvertex, Double_t tracklets) const
1280 {
1281   Double_t p = 0.;
1282   Double_t signal = 99999.;
1283   if (TMath::Abs(zvertex) > 10.) return signal;
1284   Double_t resolution[2] = {75.,65.};
1285   if (side == 0) {
1286     if (zvertex >= 5. && zvertex <= 10.) {
1287 //      p = 0.84 - exp(-1.03 - 0.31*tracklets);
1288         p = 0.88 - exp(-1.1 - 0.32*tracklets);
1289     }
1290     else if (zvertex >=-10. && zvertex <5.) {
1291 //      p = 0.82 - exp(-0.81 - 0.25*tracklets);
1292         p = 0.77 - exp(-0.88 - 0.25*tracklets);
1293    }
1294   } else if (side == 1) {
1295     if (zvertex >= -10. && zvertex < -5.) {
1296       p = 0.99 - exp(-0.74 - 0.34*tracklets);
1297     }
1298     else if (zvertex >=-5. && zvertex <10.) {
1299       p = 0.96 - exp(-0.51 - 0.27*tracklets);
1300     }
1301   } else {
1302     return signal;
1303   }
1304   Double_t pu = gRandom->Rndm();
1305   if (fDebugLevel > 1) {
1306     printf(" TofTender: T0 simul [side %d zvt %f track %f] %f [pu: %f]",side,zvertex,tracklets,p,pu);
1307     if (pu<p) printf(" --> signal will be generated: ");
1308     else printf(" --> signal wil not be generated: ");
1309   }
1310   if (pu < p) signal = gRandom->Gaus(0.,resolution[side]);
1311   if (fDebugLevel >1) Printf(" %f ",signal);
1312   return signal;
1313 }
1314
1315 void AliTOFTenderSupply::GetTrackletsForT0(AliESDEvent* event, Double_t *trkA, Double_t *trkC) const
1316 {
1317   Double_t minetaA = 0.7;
1318   Double_t maxetaA = 1.4;
1319   Double_t minetaC = -1.4;
1320   Double_t maxetaC = -0.7;
1321   AliMultiplicity *alimult = (AliMultiplicity *)event->GetMultiplicity();
1322   Int_t nTr=alimult->GetNumberOfTracklets();
1323   if (fDebugLevel > 1) Printf(" TofTender: T0 simul number of tracklets %d",nTr);
1324   Int_t nTrackletsA=0, nTrackletsC=0;
1325   for(Int_t iTr=0; iTr<nTr; iTr++){
1326     Double_t eta=alimult->GetEta(iTr);
1327     if(eta>minetaA && eta<maxetaA) nTrackletsA++;
1328     if(eta>minetaC && eta<maxetaC) nTrackletsC++;
1329     if (fDebugLevel > 1) Printf(" TofTender: T0 simul [tracklet # %d] ETA: %f %d %d",nTr,eta,nTrackletsA,nTrackletsC);
1330   }
1331   *trkA=(Double_t)nTrackletsA;
1332   *trkC=(Double_t)nTrackletsC;
1333 }
1334
1335