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