improved params for T0A efficiency
[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::kSPECIES];
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   /* set integrated times */
497   track->SetIntegratedTimes(texp);
498
499 }
500
501
502 //______________________________________________________________________________
503 void AliTOFTenderSupply::DetectRecoPass()
504 {
505   //
506   // Detect reconstruction information
507   //
508   
509   //reset information
510   fRecoPass=0;
511   
512   //Get the current file to check the reconstruction pass (UGLY, but not stored in ESD... )
513   AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
514   AliVEventHandler *inputHandler=mgr->GetInputEventHandler();
515   if (!inputHandler) return;
516   
517   TTree *tree= (TTree*)inputHandler->GetTree();
518   TFile *file= (TFile*)tree->GetCurrentFile();
519   
520   if (!file) {
521     AliFatal("Current file not found");
522     return; // coverity
523   }
524   
525   //find pass from file name (UGLY, but not stored in ESD... )
526   TString fileName(file->GetName());
527   if (fileName.Contains("pass1") ) {
528     fRecoPass=1;
529   } else if (fileName.Contains("pass2") ) {
530     fRecoPass=2;
531   } else if (fileName.Contains("pass3") ) {
532     fRecoPass=3;
533   } else if (fileName.Contains("pass4") ) {
534     fRecoPass=4;
535   } else if (fileName.Contains("pass5") ) {
536     fRecoPass=5;
537   } else if (fileName.Contains("pass6") ) {
538     fRecoPass=6;
539   }
540   if (fRecoPass == 0) {
541     AliInfo(Form("From file name %s reco pass cannot be detected",fileName.Data()));
542     AliInfo("Change file name or use SetUserRecoPass method");
543     AliFatal("------------- TOF tender cannot run with reco pass unspecified, issuing FATAL error ---------- ");
544   }
545 }
546
547
548 //______________________________________________________________________________
549 void AliTOFTenderSupply::InitGeom()
550 {
551
552   if (fGeomSet == kTRUE) return;
553
554   //  Printf("\n \n ----- calling InitGeom to fix TRD Bug ----- \n \n");
555   AliCDBManager * man = AliCDBManager::Instance();
556   man->SetDefaultStorage("raw://");
557   fCDBkey = man->SetLock(kFALSE, fCDBkey);
558   Int_t run = fTender->GetRun();
559   //  Printf(" ---------> run is %d",run);
560   man->SetRun(run);
561   AliCDBEntry *entry = (AliCDBEntry*)man->Get("GRP/Geometry/Data");
562   if (entry) {
563     AliGeomManager::LoadGeometry();
564     AliGeomManager::ApplyAlignObjsFromCDB("ITS TPC TRD TOF");
565     //    fCDBkey = man->SetLock(kTRUE, fCDBkey);
566     //    Printf("\n \n ----- Geometry loaded ------ \n \n");
567   }
568   fGeomSet=kTRUE;
569
570 }
571
572
573 //______________________________________________________________________________
574 void AliTOFTenderSupply::FixTRDBug(AliESDEvent* event)
575 //
576 // recompute texp fixing wrong dE/dx from TRD (LHC10b,c pass3)
577 //
578 {
579
580   if (fGeomSet == kFALSE) InitGeom();
581
582   //  Printf("Running FixTRD bug ");
583   /* loop over tracks */
584   AliESDtrack *track = NULL;
585   for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) {
586     track = event->GetTrack(itrk);
587     FixTRDBug(track);
588   }
589 }
590
591
592 //_____________________________________________________
593 void AliTOFTenderSupply::FixTRDBug(AliESDtrack *track)
594 {
595   // 
596   //
597   //
598
599
600     ULong_t status=track->GetStatus();
601     if (!( ( (status & AliVTrack::kITSrefit)==AliVTrack::kITSrefit ) &&
602            ( (status & AliVTrack::kTPCrefit)==AliVTrack::kTPCrefit ) &&
603            ( (status & AliVTrack::kTPCout)==AliVTrack::kTPCout ) &&
604            ( (status & AliVTrack::kTOFout)==AliVTrack::kTOFout ) &&
605            ( (status & AliVTrack::kTIME)==AliVTrack::kTIME ) ) ) return;
606
607     fIsEnteringInTRD=kFALSE;
608     fInTRD=kFALSE;
609     fIsComingOutTRD=kFALSE;
610     fOutTRD=kFALSE;
611
612     //    Printf("Track reached TOF %f",track->P());
613     Double_t correctionTimes[5] = {0.,0.,0.,0.,0.}; // to be added to the expected times
614     FindTRDFix(track, correctionTimes);
615     Double_t expectedTimes[5] = {0.,0.,0.,0.,0.}; track->GetIntegratedTimes(expectedTimes);
616     //    Printf("Exp. times: %f %f %f %f %f",
617     //     expectedTimes[0],expectedTimes[1],expectedTimes[2],expectedTimes[3],expectedTimes[4]);
618     //    Printf("Corr. times: %f %f %f %f %f",
619     //     correctionTimes[0],correctionTimes[1],correctionTimes[2],correctionTimes[3],correctionTimes[4]);
620
621     for (Int_t jj=0; jj<5; jj++) expectedTimes[jj]+=correctionTimes[jj];
622     track->SetIntegratedTimes(expectedTimes);
623   
624 }
625
626
627 //________________________________________________________________________
628 void AliTOFTenderSupply::FindTRDFix(AliESDtrack *track,Double_t *corrections)
629 {
630
631   Double_t pT = track->Pt();
632   ULong_t status=track->GetStatus();
633   Bool_t isTRDout = (status & AliVTrack::kTRDout)==AliVTrack::kTRDout;
634
635   Double_t length = 0.;
636
637   Double_t xyzIN[3]={0.,0.,0.};
638   fIsEnteringInTRD = track->GetXYZAt(fRhoTRDin,fMagField,xyzIN);
639
640   Double_t xyzOUT[3]={0.,0.,0.};
641   fIsComingOutTRD = track->GetXYZAt(fRhoTRDout,fMagField,xyzOUT);
642
643   if (fIsEnteringInTRD && fIsComingOutTRD) {
644
645
646     Double_t phiIN = TMath::Pi()+TMath::ATan2(-xyzIN[1],-xyzIN[0]);
647     phiIN *= TMath::RadToDeg();
648     fInTRD = ( (phiIN>=  0. && phiIN<= 40.) ||
649                (phiIN>=140. && phiIN<=220.) ||
650                (phiIN>=340. && phiIN<=360.) ); // TRD SMs installed @ 2010
651
652     Double_t phiOUT = TMath::Pi()+TMath::ATan2(-xyzOUT[1],-xyzOUT[0]);
653     phiOUT *= TMath::RadToDeg();
654     fOutTRD = ( (phiOUT>=  0. && phiOUT<= 40.) ||
655                 (phiOUT>=140. && phiOUT<=220.) ||
656                 (phiOUT>=340. && phiOUT<=360.) ); // TRD SMs installed @ 2010
657
658     length = 0.;
659
660     if (fInTRD || fOutTRD) {
661
662       if ( ( fInTRD && fOutTRD ) || ( fInTRD && !fOutTRD ) ) {
663         length = EstimateLengthInTRD1(track);
664       } else if ( !fInTRD && fOutTRD ) {
665         length = EstimateLengthInTRD2(track);
666       }
667
668     } else { // ( !fInTRD && !fOutTRD )
669
670       length = EstimateLengthOutTRD(track);
671
672     }
673
674   }
675   //  Printf("estimated length in TRD %f [isTRDout %d]",length,isTRDout);
676   CorrectDeltaTimes(pT,length,isTRDout,corrections);
677
678 }
679
680 //________________________________________________________________________
681 void AliTOFTenderSupply::CorrectDeltaTimes(Double_t pT,
682                                                 Double_t length,
683                                                 Bool_t flagTRDout,
684                                                 Double_t *corrections)
685 {
686
687   corrections[2] = CorrectExpectedPionTime(pT,length,flagTRDout);
688   corrections[0] = corrections[2]; // x electrons used pion corrections
689   corrections[1] = corrections[2]; // x muons used pion corrections
690   corrections[3] = CorrectExpectedKaonTime(pT,length,flagTRDout);
691   corrections[4] = CorrectExpectedProtonTime(pT,length,flagTRDout);
692
693 }
694
695 //________________________________________________________________________
696 Double_t AliTOFTenderSupply::CorrectExpectedPionTime(Double_t pT,
697                                                           Double_t length,
698                                                           Bool_t isTRDout)
699 {
700   // correction for expected time for pions
701
702   Double_t delta=0.;
703
704   //  Printf("Flags Ent In ComingOut Out %d %d %d %d",fIsEnteringInTRD,fInTRD,fIsComingOutTRD,fOutTRD);
705   if (!fIsEnteringInTRD || !fIsComingOutTRD) { // zone 5
706
707     Float_t p[2]={0.,0.};
708
709     if (isTRDout) {
710
711       if (pT<0.30) {
712         p[0] = 180.; p[1] = 0.;
713       } else if (pT>=0.30 && pT<0.35) {
714         p[0] = 740.; p[1] = -1800.;
715       } else if (pT>=0.35 && pT<0.40) {
716         p[0] = 488.; p[1] =-1080.;
717       } else if (pT>=0.40 && pT<0.45) {
718         p[0] = 179.; p[1] = -307.;
719       } else if (pT>=0.45 && pT<0.50) {
720         p[0] =  97.; p[1] = -123.;
721       } else { //if (pT>=0.50)
722         p[0] = 120.; p[1] = -172.;
723       }
724
725     } else {
726
727       if (pT<0.30) {
728         p[0] =  70.; p[1] =    0.;
729       } else if (pT>=0.30 && pT<0.35) {
730         p[0] = 339.; p[1] = -927.;
731       } else if (pT>=0.35 && pT<0.40) {
732         p[0] =  59.; p[1] = -121.;
733       } else if (pT>=0.40 && pT<0.50) {
734         p[0] =  21.; p[1] =  -24.;
735       } else { //if (pT>=0.50)
736         p[0] =  42.; p[1] =  -67.;
737       }
738
739     }
740
741     delta = p[0]+p[1]*pT;
742     //    Printf("Pion time: %f %f %f %f",p[0],p[1],length,delta);
743
744   } else {
745
746     Float_t p[2] = {0.,0.};
747
748     if ( fInTRD &&  fOutTRD) { // zone 1
749
750       if (isTRDout) {
751
752         if (length<130.) {
753           p[0] = 0.; p[1] = 0.;
754         } else if (length>=130. && length<170.) {
755           p[0] = -20.5; p[1] = 0.25;
756         } else {//if (length>=170.)
757           p[0] = 22.; p[1] = 0.;
758         }
759
760       } else { // !isTRDout
761
762         p[0] = 20.; p[1] = 0.;
763
764       }
765
766     } else if (!fInTRD && !fOutTRD) { // zone 2
767
768       p[0] = 0.; p[1] = 0.;
769
770     } else if ( fInTRD &&  !fOutTRD) { // zone 3
771
772       if (isTRDout) {
773
774         if (length< 75.) {
775           p[0] = 17.; p[1] =  0.;
776         } else if (length>= 75. && length< 95.) {
777           p[0] = 81.; p[1] = -0.85;
778         } else if (length>= 95. && length<155.) {
779           p[0] =  0.; p[1] =  0.;
780         } else {//if (length>=155.)
781           p[0] = 10.; p[1] =  0.;
782         }
783
784       } else { // !isTRDout
785
786         p[0] = 0.; p[1] = 0.;
787
788       }
789
790     } else if (!fInTRD &&  fOutTRD) { // zone 4
791
792       if (isTRDout) {
793
794         if (length<80.) {
795           p[0] =  0.; p[1] =  0.;
796         } else {//if (length>=80.)
797           p[0] = 10.; p[1] =  0.;
798         }
799
800       } else { // !isTRDout
801
802         if (length<30.) {
803           p[0] =  0.; p[1] =  0.;
804         } else {//if (length>=30.)
805           p[0] =  6.; p[1] =  0.;
806         }
807
808       }
809
810     }
811
812     delta = p[0]+p[1]*length;
813     //    Printf("Pion time: %f %f %f %f",p[0],p[1],length,delta);
814
815   }
816
817   return delta;
818
819 }
820
821 //________________________________________________________________________
822 Double_t AliTOFTenderSupply::CorrectExpectedKaonTime(Double_t pT,
823                                                           Double_t length,
824                                                           Bool_t isTRDout)
825 {
826   // correction for expected time for kaons
827
828   Double_t delta=0.;
829   //  Printf("Flags Ent In ComingOut Out %d %d %d %d",fIsEnteringInTRD,fInTRD,fIsComingOutTRD,fOutTRD);
830
831   if (!fIsEnteringInTRD || !fIsComingOutTRD) { // zone 5
832
833     Float_t p[2]={0.,0.};
834
835     if (isTRDout) {
836
837       if (pT<0.4) {
838         p[0] =  900.; p[1] = 0.;
839       } else if (pT>=0.40 && pT<0.45) {
840         p[0] = 3100.; p[1] = -6000.;
841       } else if (pT>=0.45 && pT<0.50) {
842         p[0] = 1660.; p[1] = -2800.;
843       } else if (pT>=0.50 && pT<0.55) {
844         p[0] =  860.; p[1] = -1200.;
845       } else { //if (pT>=0.55)
846         p[0] =  200.; p[1] = 0.;
847       }
848
849     } else {
850
851       if (pT<0.30) {
852         p[0] =   0.; p[1] =    0.;
853       } else if (pT>=0.30 && pT<0.32) {
854         p[0] = 570.; p[1] =    0.;
855       } else if (pT>=0.32 && pT<0.35) {
856         p[0] = 3171.; p[1] = -8133.;
857       } else if (pT>=0.35 && pT<0.40) {
858         p[0] = 1815.; p[1] = -4260.;
859       } else if (pT>=0.40 && pT<0.45) {
860         p[0] =  715.; p[1] =  -1471.;
861       } else if (pT>=0.45 && pT<0.50) {
862         p[0] =  233.; p[1] =  -407.;
863       } else if (pT>=0.50 && pT<0.55) {
864         p[0] =  408.; p[1] =  -752.;
865       } else { //if (pT>=0.55)
866         p[0] =  408.-752.*0.55; p[1] = 0.;
867       }
868
869     }
870
871     delta = p[0]+p[1]*pT;
872     //    Printf("Kaon time: %f %f %f %f",p[0],p[1],length,delta);
873
874   } else {
875
876     Float_t p[2] = {0.,0.};
877
878     if ( fInTRD &&  fOutTRD) { // zone 1
879
880       if (isTRDout) {
881
882         if (length<95.) {
883           p[0] = 20.; p[1] = 0.;
884         } else if (length>=95. && length<195.) {
885           p[0] = -24.0; p[1] = 0.10+0.0041*length;
886         } else {//if (length>=195.)
887           p[0] =  150.; p[1] = 0.;
888         }
889
890       } else { // !isTRDout
891
892         p[0] = 40.; p[1] = 0.;
893
894       }
895
896     } else if (!fInTRD && !fOutTRD) { // zone 2
897
898       p[0] = 0.; p[1] = 0.;
899
900     } else if ( fInTRD &&  !fOutTRD) { // zone 3
901
902       if (isTRDout) {
903
904         if (length< 15.) {
905           p[0] = 180.; p[1] =  0.;
906         } else if (length>= 15. && length< 55.) {
907           p[0] = 215.; p[1] = -2.5;
908         } else {//if (length>=55.)
909           p[0] = 78.; p[1] =  0.;
910         }
911
912       } else { // !isTRDout
913
914         p[0] = 0.; p[1] = 0.;
915
916       }
917
918     } else if (!fInTRD &&  fOutTRD) { // zone 4
919
920       if (isTRDout) {
921
922         if (length< 55.) {
923           p[0] =  0.; p[1] =  0.;
924         } else if (length>= 55. && length<115.) {
925           p[0] = -85.; p[1] = 1.9;
926         } else {//if (length>=115.)
927           p[0] = 100.; p[1] =  0.;
928         }
929
930       } else { // !isTRDout
931
932         p[0] =  0.; p[1] =  0.;
933
934       }
935
936     }
937
938     delta = p[0]+p[1]*length;
939     //    Printf("Kaon time: %f %f %f %f",p[0],p[1],length,delta);
940
941   }
942
943   return delta;
944
945 }
946
947 //________________________________________________________________________
948 Double_t AliTOFTenderSupply::CorrectExpectedProtonTime(Double_t pT,
949                                                             Double_t length,
950                                                             Bool_t isTRDout)
951 {
952   // correction for expected time for protons
953
954   Double_t delta=0.;
955   //  Printf("Flags Ent In ComingOut Out %d %d %d %d",fIsEnteringInTRD,fInTRD,fIsComingOutTRD,fOutTRD);
956
957   if (!fIsEnteringInTRD || !fIsComingOutTRD) { // zone 5
958     Float_t p[2]={0.,0.};
959
960
961     if (isTRDout) {
962
963       if (pT<0.375) {
964         p[0] = 1000.; p[1] = 0.;
965       } else if (pT>=0.375 && pT<0.45) {
966         p[0] = 1500.; p[1] = 0.;
967       } else if (pT>=0.45 && pT<0.50) {
968         p[0] = 4650.; p[1] = -7000.;
969       } else if (pT>=0.50 && pT<0.55) {
970         p[0] = 3150.; p[1] = -4000.;
971       } else { //if (pT>=0.55)
972         p[0] = 3150. -4000.*0.55; p[1] = 0.;
973       }
974
975     } else {
976
977       if (pT<0.32) {
978         p[0] = 2963.-5670.*0.032; p[1] = 0.;
979       } else if (pT>=0.32 && pT<0.35) {
980         p[0] = 2963.; p[1] =  -5670.;
981       } else if (pT>=0.35 && pT<0.40) {
982         p[0] = 4270.; p[1] =  -9400.;
983       } else if (pT>=0.40 && pT<0.45) {
984         p[0] = 1550.; p[1] =  -2600.;
985       } else if (pT>=0.45 && pT<0.50) {
986         p[0] = 1946.; p[1] =  -3480.;
987       } else if (pT>=0.50 && pT<0.55) {
988         p[0] = 1193.; p[1] =  -1974.;
989       } else { //if (pT>=0.55)
990         p[0] = 1193.-1974.*0.55; p[1] = 0.;
991       }
992
993     }
994
995     delta = p[0]+p[1]*pT;
996     //    Printf("Proton time: %f %f %f %f",p[0],p[1],length,delta);
997
998   } else {
999
1000     Float_t p[2] = {0.,0.};
1001
1002     if ( fInTRD &&  fOutTRD) { // zone 1
1003
1004       if (isTRDout) {
1005
1006         if (length<90.) {
1007           p[0] = 0.; p[1] = 0.;
1008         } else if (length>=90. && length<200.) {
1009           p[0] = 1063.; p[1] = -32.+0.30*length-0.00072*length*length;
1010         } else {//if (length>=200.)
1011           p[0] =  900.; p[1] = 0.;
1012         }
1013
1014       } else { // !isTRDout
1015
1016         p[0] = 80.; p[1] = 0.;
1017
1018       }
1019
1020     } else if (!fInTRD && !fOutTRD) { // zone 2
1021
1022       if (isTRDout) {
1023         p[0] = 0.; p[1] = 0.;
1024       } else {
1025         if (length<125.) {
1026           p[0] = 0.; p[1] = 0.;
1027         } else if (length>=125. && length<180.) {
1028           p[0] = -132.; p[1] = 1.3;
1029         } else {
1030           p[0] = 100.; p[1] = 0.;
1031         }
1032
1033       }
1034
1035     } else if ( fInTRD &&  !fOutTRD) { // zone 3
1036
1037       if (isTRDout) {
1038
1039         if (length< 30.) {
1040           p[0] = 670.; p[1] =  0.;
1041         } else if (length>= 30. && length<155.) {
1042           p[0] = 944.; p[1] = -11.+0.064*length;
1043         } else {//if (length>=155.)
1044           p[0] = 780.; p[1] =  0.;
1045         }
1046
1047       } else { // !isTRDout
1048
1049         if (length< 30.) {
1050           p[0] = 140.; p[1] = -4.5;
1051         } else {
1052           p[0] = 0.; p[1] = 0.;
1053         }
1054       
1055       }
1056
1057     } else if (!fInTRD &&  fOutTRD) { // zone 4
1058
1059       if (isTRDout) {
1060
1061         if (length< 45.) {
1062           p[0] = 130.; p[1] =  0.;
1063         } else if (length>= 45. && length<120.) {
1064           p[0] = -190.; p[1] = 6.5;
1065         } else {//if (length>=120.)
1066           p[0] = 750.; p[1] =  0.;
1067         }
1068
1069       } else { // !isTRDout
1070
1071         if (length<75.5) {
1072           p[0] =    0.; p[1] =  0.;
1073         } else if (length>= 75.5 && length<90.) {
1074           p[0] = -830.; p[1] = 11.;
1075         } else {
1076           p[0] =  160.; p[1] =  0.;
1077         }
1078
1079       }
1080
1081     }
1082
1083     delta = p[0]+p[1]*length;
1084     //    Printf("Proton time: %f %f %f %f",p[0],p[1],length,delta);
1085
1086   }
1087
1088   return delta;
1089
1090 }
1091
1092 //________________________________________________________________________
1093 Double_t AliTOFTenderSupply::EstimateLengthInTRD1(AliESDtrack *track)
1094 {
1095
1096   Double_t xyz0[3]={0.,0.,0.};
1097   Bool_t stayInTRD = track->GetXYZAt(fRhoTRDin,fMagField,xyz0);
1098
1099   Double_t phi0 = TMath::Pi()+TMath::ATan2(-xyz0[1],-xyz0[0]);
1100   phi0 *= TMath::RadToDeg();
1101   stayInTRD = stayInTRD && ( (phi0>=  0. && phi0<= 40.) ||
1102                              (phi0>=140. && phi0<=220.) ||
1103                              (phi0>=340. && phi0<=360.) );
1104
1105   Double_t trackLengthInTRD = 0.;
1106   Int_t iStep=0;
1107
1108   Double_t b[3];track->GetBxByBz(b);
1109
1110   Double_t xyz1[3]={0.,0.,0.};
1111   Double_t rho = fRhoTRDin;
1112   while (stayInTRD && rho<=fRhoTRDout) {
1113     iStep++;
1114     rho += fStep;
1115
1116     for (Int_t ii=0; ii<3; ii++) xyz1[ii]=0.;
1117     stayInTRD = track->GetXYZAt(rho,fMagField,xyz1);
1118     Double_t phi1 = TMath::Pi()+TMath::ATan2(-xyz1[1],-xyz1[0]);
1119     phi1 *= TMath::RadToDeg();
1120     stayInTRD = stayInTRD && ( (phi1>=  0. && phi1<= 40.) ||
1121                                (phi1>=140. && phi1<=220.) ||
1122                                (phi1>=340. && phi1<=360.) );
1123
1124     Double_t l2  = TMath::Sqrt((xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0]) +
1125                                (xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1]) +
1126                                (xyz1[2]-xyz0[2])*(xyz1[2]-xyz0[2]));
1127     trackLengthInTRD += l2;
1128
1129     for (Int_t ii=0; ii<3; ii++) xyz0[ii]=xyz1[ii];
1130   }
1131
1132   return trackLengthInTRD;
1133
1134 }
1135
1136 //________________________________________________________________________
1137 Double_t AliTOFTenderSupply::EstimateLengthInTRD2(AliESDtrack *track)
1138 {
1139
1140   Double_t xyz0[3]={0.,0.,0.};
1141   Bool_t stayInTRD = track->GetXYZAt(fRhoTRDout,fMagField,xyz0);
1142
1143   Double_t phi0 = TMath::Pi()+TMath::ATan2(-xyz0[1],-xyz0[0]);
1144   phi0 *= TMath::RadToDeg();
1145   stayInTRD = stayInTRD && ( (phi0>=  0. && phi0<= 40.) ||
1146                              (phi0>=140. && phi0<=220.) ||
1147                              (phi0>=340. && phi0<=360.) );
1148
1149   Double_t trackLengthInTRD = 0.;
1150   Int_t iStep=0;
1151
1152   Double_t b[3];track->GetBxByBz(b);
1153
1154   Double_t xyz1[3]={0.,0.,0.};
1155   Double_t rho = fRhoTRDout;
1156   while (stayInTRD && rho>=fRhoTRDin) {
1157     iStep++;
1158     rho -= fStep;
1159
1160     for (Int_t ii=0; ii<3; ii++) xyz1[ii]=0.;
1161     stayInTRD = track->GetXYZAt(rho,fMagField,xyz1);
1162     Double_t phi1 = TMath::Pi()+TMath::ATan2(-xyz1[1],-xyz1[0]);
1163     phi1 *= TMath::RadToDeg();
1164     stayInTRD = stayInTRD && ( (phi1>=  0. && phi1<= 40.) ||
1165                                (phi1>=140. && phi1<=220.) ||
1166                                (phi1>=340. && phi1<=360.) );
1167
1168     Double_t l2  = TMath::Sqrt((xyz0[0]-xyz1[0])*(xyz0[0]-xyz1[0]) +
1169                                (xyz0[1]-xyz1[1])*(xyz0[1]-xyz1[1]) +
1170                                (xyz0[2]-xyz1[2])*(xyz0[2]-xyz1[2]));
1171     trackLengthInTRD += l2;
1172
1173     for (Int_t ii=0; ii<3; ii++) xyz0[ii]=xyz1[ii];
1174   }
1175
1176   return trackLengthInTRD;
1177
1178 }
1179
1180 //________________________________________________________________________
1181 Double_t AliTOFTenderSupply::EstimateLengthOutTRD(AliESDtrack *track)
1182 {
1183
1184   Double_t xyz0[3]={0.,0.,0.};
1185   Bool_t stayInTRD = track->GetXYZAt(fRhoTRDin,fMagField,xyz0);
1186
1187   Double_t phi0 = TMath::Pi()+TMath::ATan2(-xyz0[1],-xyz0[0]);
1188   phi0 *= TMath::RadToDeg();
1189   stayInTRD = stayInTRD && !( (phi0>=  0. && phi0<= 40.) ||
1190                               (phi0>=140. && phi0<=220.) ||
1191                               (phi0>=340. && phi0<=360.) );
1192
1193   Double_t trackLengthInTRD = 0.;
1194   Int_t iStep=0;
1195
1196   Double_t b[3];track->GetBxByBz(b);
1197
1198   Double_t xyz1[3]={0.,0.,0.};
1199   Double_t rho = fRhoTRDin;
1200   while (stayInTRD && rho<=fRhoTRDout) {
1201     iStep++;
1202     rho += fStep;
1203
1204     for (Int_t ii=0; ii<3; ii++) xyz1[ii]=0.;
1205     stayInTRD = track->GetXYZAt(rho,fMagField,xyz1);
1206     Double_t phi1 = TMath::Pi()+TMath::ATan2(-xyz1[1],-xyz1[0]);
1207     phi1 *= TMath::RadToDeg();
1208     stayInTRD = stayInTRD && !( (phi1>=  0. && phi1<= 40.) ||
1209                                 (phi1>=140. && phi1<=220.) ||
1210                                 (phi1>=340. && phi1<=360.) );
1211
1212     Double_t l2  = TMath::Sqrt((xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0]) +
1213                                (xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1]) +
1214                                (xyz1[2]-xyz0[2])*(xyz1[2]-xyz0[2]));
1215     trackLengthInTRD += l2;
1216
1217     for (Int_t ii=0; ii<3; ii++) xyz0[ii]=xyz1[ii];
1218   }
1219
1220   return trackLengthInTRD;
1221
1222 }
1223
1224 //________________________________________________________________________
1225 Int_t AliTOFTenderSupply::GetOCDBVersion(Int_t runNo)
1226 {
1227   Int_t verNo = -1;
1228   if ( (runNo>=118503 && runNo <=121040) )  { // LHC10C
1229     if (fRecoPass == 2) {                     // on pass2
1230       if (runNo >= 119159 && runNo <= 119163) verNo=3;
1231       else if (runNo >= 119837 && runNo <= 119934) verNo=4;
1232       else if (runNo >= 120067 && runNo <= 120244) verNo=4;
1233       else if (runNo >= 120503 && runNo <= 120505) verNo=4;
1234       else if (runNo >= 120616 && runNo <= 120671) verNo=4;
1235       else if (runNo >= 120741 && runNo <= 120829) verNo=4;
1236     }
1237   }
1238   return verNo;
1239 }
1240
1241
1242 //__________________________________________________________________________
1243 void AliTOFTenderSupply::LoadTOFPIDParams(Int_t runNumber)
1244 {
1245   //
1246   // Load the TOF pid params from the OADB
1247   //
1248
1249   if (fTOFPIDParams) delete fTOFPIDParams;
1250   fTOFPIDParams=0x0;
1251   
1252   TFile *oadbf = new TFile("$ALICE_ROOT/OADB/COMMON/PID/data/TOFPIDParams.root");
1253   if (oadbf && oadbf->IsOpen()) {
1254     AliInfo("Loading TOF Params from $ALICE_ROOT/OADB/COMMON/PID/data/TOFPIDParams.root");
1255     AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
1256     if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(runNumber,"TOFparams"));
1257     oadbf->Close();
1258     delete oadbc;
1259   }
1260   delete oadbf;
1261
1262   if (!fTOFPIDParams) {
1263     AliError("TOFPIDParams.root not found in $ALICE_ROOT/OADB/COMMON/PID/data !!");
1264     fTOFPIDParams = new AliTOFPIDParams;
1265     fTOFPIDParams->SetTOFresolution(90.);
1266     fTOFPIDParams->SetStartTimeMethod(AliESDpid::kTOF_T0);
1267   }  
1268 }
1269
1270
1271 //__________________________________________________________________________
1272 Double_t AliTOFTenderSupply::SampleT0Signal(Int_t side, Double_t zvertex, Double_t tracklets) const
1273 {
1274   Double_t p = 0.;
1275   Double_t signal = 99999.;
1276   if (TMath::Abs(zvertex) > 10.) return signal;
1277   Double_t resolution[2] = {75.,65.};
1278   if (side == 0) {
1279     if (zvertex >= 5. && zvertex <= 10.) {
1280 //      p = 0.84 - exp(-1.03 - 0.31*tracklets);
1281         p = 0.88 - exp(-1.1 - 0.32*tracklets);
1282     }
1283     else if (zvertex >=-10. && zvertex <5.) {
1284 //      p = 0.82 - exp(-0.81 - 0.25*tracklets);
1285         p = 0.77 - exp(-0.88 - 0.25*tracklets);    
1286    }
1287   } else if (side == 1) {
1288     if (zvertex >= -10. && zvertex < -5.) {
1289       p = 0.99 - exp(-0.74 - 0.34*tracklets);
1290     }
1291     else if (zvertex >=-5. && zvertex <10.) {
1292       p = 0.96 - exp(-0.51 - 0.27*tracklets);
1293     }
1294   } else {
1295     return signal;
1296   }
1297   Double_t pu = gRandom->Rndm();
1298   if (fDebugLevel > 1) {
1299     printf(" TofTender: T0 simul [side %d zvt %f track %f] %f [pu: %f]",side,zvertex,tracklets,p,pu);
1300     if (pu<p) printf(" --> signal will be generated: ");
1301     else printf(" --> signal wil not be generated: ");
1302   }
1303   if (pu < p) signal = gRandom->Gaus(0.,resolution[side]);
1304   if (fDebugLevel >1) Printf(" %f ",signal);
1305   return signal;
1306 }
1307
1308 void AliTOFTenderSupply::GetTrackletsForT0(AliESDEvent* event, Double_t *trkA, Double_t *trkC) const
1309 {
1310   Double_t minetaA = 0.7;
1311   Double_t maxetaA = 1.4;
1312   Double_t minetaC = -1.4;
1313   Double_t maxetaC = -0.7; 
1314   AliMultiplicity *alimult = (AliMultiplicity *)event->GetMultiplicity(); 
1315   Int_t nTr=alimult->GetNumberOfTracklets();
1316   if (fDebugLevel > 1) Printf(" TofTender: T0 simul number of tracklets %d",nTr);
1317   Int_t nTrackletsA=0, nTrackletsC=0;
1318   for(Int_t iTr=0; iTr<nTr; iTr++){
1319     Double_t eta=alimult->GetEta(iTr);
1320     if(eta>minetaA && eta<maxetaA) nTrackletsA++;
1321     if(eta>minetaC && eta<maxetaC) nTrackletsC++;
1322     if (fDebugLevel > 1) Printf(" TofTender: T0 simul [tracklet # %d] ETA: %f %d %d",nTr,eta,nTrackletsA,nTrackletsC);
1323   }
1324   *trkA=(Double_t)nTrackletsA;
1325   *trkC=(Double_t)nTrackletsC;
1326 }
1327
1328