coverity fix
[u/mrichter/AliRoot.git] / ANALYSIS / 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);
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