coverity fix
[u/mrichter/AliRoot.git] / ANALYSIS / TenderSupplies / AliTOFTenderSupply.cxx
CommitLineData
ee981ab3 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// //
dc369c13 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 //
ee981ab3 25///////////////////////////////////////////////////////////////////////////////
fd611537 26#include <TFile.h>
27#include <TChain.h>
28
dc369c13 29#include <TMath.h>
30#include <TRandom.h>
ee981ab3 31#include <AliLog.h>
32#include <AliESDEvent.h>
88ccd4cf 33#include <AliESDtrack.h>
ee981ab3 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
fd611537 42#include <AliGeomManager.h>
ee981ab3 43#include <AliCDBManager.h>
44#include <AliCDBEntry.h>
45#include <AliT0CalibSeasonTimeShift.h>
46
47#include "AliTOFTenderSupply.h"
48
e26aa0bb 49ClassImp(AliTOFTenderSupply)
50
dc369c13 51Float_t AliTOFTenderSupply::fgT0Aresolution = 75.;
52Float_t AliTOFTenderSupply::fgT0Cresolution = 65.;
53
ee981ab3 54AliTOFTenderSupply::AliTOFTenderSupply() :
55 AliTenderSupply(),
56 fESDpid(0x0),
57 fIsMC(kFALSE),
dc369c13 58 fTimeZeroType(AliESDpid::kBest_T0),
ee981ab3 59 fCorrectExpTimes(kTRUE),
fd611537 60 fCorrectTRDBug(kFALSE),
88ccd4cf 61 fLHC10dPatch(kFALSE),
3eb060be 62 fT0DetectorAdjust(kFALSE),
dc369c13 63 fDebugLevel(0),
63f693ef 64 fAutomaticSettings(kTRUE),
fd611537 65 fRecoPass(0),
66 fUserRecoPass(0),
ee981ab3 67 fTOFCalib(0x0),
68 fTOFT0maker(0x0),
63f693ef 69 fTOFres(100.),
fd611537 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
ee981ab3 82
dc369c13 83
ee981ab3 84{
85 //
86 // default ctor
87 //
ee981ab3 88 fT0shift[0] = 0;
89 fT0shift[1] = 0;
90 fT0shift[2] = 0;
91 fT0shift[3] = 0;
92}
93
94//_____________________________________________________
95AliTOFTenderSupply::AliTOFTenderSupply(const char *name, const AliTender *tender) :
96 AliTenderSupply(name,tender),
97 fESDpid(0x0),
98 fIsMC(kFALSE),
dc369c13 99 fTimeZeroType(AliESDpid::kBest_T0),
ee981ab3 100 fCorrectExpTimes(kTRUE),
fd611537 101 fCorrectTRDBug(kFALSE),
88ccd4cf 102 fLHC10dPatch(kFALSE),
3eb060be 103 fT0DetectorAdjust(kFALSE),
dc369c13 104 fDebugLevel(0),
63f693ef 105 fAutomaticSettings(kTRUE),
fd611537 106 fRecoPass(0),
107 fUserRecoPass(0),
ee981ab3 108 fTOFCalib(0x0),
109 fTOFT0maker(0x0),
63f693ef 110 fTOFres(100.),
fd611537 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)
ee981ab3 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//_____________________________________________________
135void AliTOFTenderSupply::Init()
136{
ee981ab3 137
3eb060be 138 Bool_t tenderUnsupported = kFALSE;
dc369c13 139 // Initialise TOF tender (this is called at each detected run change)
b3e3d195 140 AliLog::SetClassDebugLevel("AliTOFTenderSupply",10);
141
dc369c13 142 // Setup PID object, check for MC, set AliTOFcalib and TOFT0 maker conf
63f693ef 143 Int_t run = fTender->GetRun();
144 if (run == 0) return; // to skip first init, when we don't have yet a run number
ee981ab3 145
fd611537 146 fGeomSet=kFALSE;
147
63f693ef 148 if (fAutomaticSettings) {
fd611537 149 if (fUserRecoPass == 0) DetectRecoPass();
150 else fRecoPass = fUserRecoPass;
3eb060be 151 if (run<114737) {
152 tenderUnsupported = kTRUE;
153 }
154 else if (run>=114737&&run<=117223) { //period="LHC10B";
fd611537 155 if (fRecoPass == 2) {fCorrectExpTimes=kTRUE; fCorrectTRDBug=kFALSE;}
156 else if (fRecoPass == 3) {fCorrectExpTimes=kFALSE; fCorrectTRDBug=kTRUE;}
63f693ef 157 fLHC10dPatch=kFALSE;
158 fTOFres=100.;
159 fTimeZeroType=AliESDpid::kTOF_T0;
160 fT0IntercalibrationShift = 0;
3eb060be 161 fT0DetectorAdjust=kTRUE;
63f693ef 162 }
163 else if (run>=118503&&run<=121040) { //period="LHC10C";
fd611537 164 if (fRecoPass == 2) {fCorrectExpTimes=kTRUE; fCorrectTRDBug=kFALSE;}
165 else if (fRecoPass == 3) {fCorrectExpTimes=kFALSE; fCorrectTRDBug=kTRUE;}
63f693ef 166 fLHC10dPatch=kFALSE;
167 fTOFres=100.;
168 fTimeZeroType=AliESDpid::kTOF_T0;
169 fT0IntercalibrationShift = 0;
3eb060be 170 fT0DetectorAdjust=kFALSE;
63f693ef 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;
3eb060be 178 fT0DetectorAdjust=kTRUE;
63f693ef 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.;
3eb060be 186 fT0DetectorAdjust=kTRUE;
63f693ef 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.;
3eb060be 194 fT0DetectorAdjust=kTRUE;
63f693ef 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.;
3eb060be 203 fT0DetectorAdjust=kTRUE;
63f693ef 204 AliWarning("TOF tender not supported for LHC10G period!! Settings are just a guess!!");
205 }
0715e776 206 else if (run>=136851&&run<=139517) { //period="LHC10H" - pass2;
63f693ef 207 fCorrectExpTimes=kFALSE;
0715e776 208 fLHC10dPatch=kFALSE;
63f693ef 209 fTOFres=90.;
210 fTimeZeroType=AliESDpid::kTOF_T0;
211 fT0IntercalibrationShift = 0.;
3eb060be 212 fT0DetectorAdjust=kTRUE;
63f693ef 213 }
214 else if (run>=139699) { //period="LHC11A";
3eb060be 215 /*
63f693ef 216 fCorrectExpTimes=kFALSE;
217 fLHC10dPatch=kFALSE;
218 fTOFres=100.;
219 fTimeZeroType=AliESDpid::kBest_T0;
220 fT0IntercalibrationShift = 0.;
3eb060be 221 fT0DetectorAdjust=kFALSE;
63f693ef 222 AliWarning("TOF tender not supported for LHC11A period!! Settings are just a guess!!");
3eb060be 223 */
224 AliError("TOF tender not supported for 2011 data!!!!!");
225 tenderUnsupported = kTRUE;
63f693ef 226 }
227 }
ee981ab3 228
3eb060be 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
fd611537 241 // Check if another tender wagon already created the esd pid object
ee981ab3 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
dc369c13 249 // Even if the user didn't set fIsMC, we force it on if we find the MC handler
ee981ab3 250 AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
dc369c13 251 if (mgr->GetMCtruthEventHandler() && !(fIsMC) ) {
252 AliWarning("This ESD is MC, fIsMC found OFF: fIsMC turned ON");
253 fIsMC=kTRUE;
254 }
ee981ab3 255
dc369c13 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
fd611537 261 // (this is done for LHC10b/LHC10c pass2)
dc369c13 262
263 // Configure TOFT0 maker class
74830383 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
dc369c13 266 fTOFT0maker->SetTimeResolution(fTOFres); // set TOF resolution for the PID
fd611537 267 fTOFT0maker->SetTOFT0algorithm(2);
ee981ab3 268
dc369c13 269 AliInfo("|******************************************************|");
270 AliInfo(Form("| Alice TOF Tender Initialisation (Run %d) |",fTender->GetRun()));
271 AliInfo("| Settings: |");
272 AliInfo(Form("| Correct Exp Times : %d |",fCorrectExpTimes));
fd611537 273 AliInfo(Form("| Correct TRD Bug : %d |",fCorrectTRDBug));
dc369c13 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));
3eb060be 278 AliInfo(Form("| T0 detector offsets applied : %d |",fT0DetectorAdjust));
63f693ef 279 AliInfo(Form("| TOF/T0 intecalibration shift : %5.2f (ps) |",fT0IntercalibrationShift));
dc369c13 280 AliInfo("|******************************************************|");
e4e3a90d 281
282
ee981ab3 283}
284
285//_____________________________________________________
286void AliTOFTenderSupply::ProcessEvent()
287{
288 //
dc369c13 289 // Use updated calibrations for TOF and T0, reapply PID information
290 // For MC: timeZero sampling and additional smearing for T0
ee981ab3 291
dc369c13 292 if (fDebugLevel > 1) AliInfo("process event");
ee981ab3 293
294 AliESDEvent *event=fTender->GetEvent();
295 if (!event) return;
dc369c13 296 if (fDebugLevel > 1) AliInfo("event read");
297
e4e3a90d 298
ee981ab3 299
dc369c13 300 if (fTender->RunChanged()){
301
63f693ef 302 Init();
303
ee981ab3 304 fTOFCalib->Init(fTender->GetRun());
305
306 if(event->GetT0TOF()){ // read T0 detector correction from OCDB
307 // OCDB instance
3eb060be 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 }
dc369c13 324 } else {
63f693ef 325 for (Int_t i=0;i<4;i++) fT0shift[i]=0;
dc369c13 326 }
ee981ab3 327 }
328 }
88ccd4cf 329
88ccd4cf 330
dc369c13 331 fTOFCalib->CalibrateESD(event); //recalculate TOF signal (no harm for MC, see settings inside init)
fd611537 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
dc369c13 337
338 Double_t startTime = 0.;
fd611537 339 if(fIsMC) startTime = fTOFCalib->TuneForMC(event,fTOFres); // this is for old MC when we didn't jitter startTime in MC
dc369c13 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));
fd611537 343
dc369c13 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.);
ee981ab3 350
dc369c13 351 if(!fIsMC){ // data: apply shifts to align around Zero
220433fd 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]);
dc369c13 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));
220433fd 379 }
dc369c13 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.);
ee981ab3 384 }
dc369c13 385 if (fDebugLevel > 1) Printf(" TofTender: T0 time (FINAL) %f %f %f",event->GetT0TOF(0),event->GetT0TOF(1),event->GetT0TOF(2));
ee981ab3 386
dc369c13 387 //compute timeZero of the event via TOF-TO
ee981ab3 388 fTOFT0maker->ComputeT0TOF(event);
389 fTOFT0maker->WriteInESD(event);
390
fd611537 391 //set preferred startTime
ee981ab3 392 fESDpid->SetTOFResponse(event, (AliESDpid::EStartTimeType_t)fTimeZeroType);
dc369c13 393
fd611537 394 // recalculate PID probabilities
dc369c13 395 // this is for safety, especially if the user doesn't attach a PID tender after TOF tender
ee981ab3 396 Int_t ntracks=event->GetNumberOfTracks();
fd611537 397 // AliESDtrack *track = NULL;
398 // Float_t tzeroTrack = 0;
ee981ab3 399 for(Int_t itrack = 0; itrack < ntracks; itrack++){
fd611537 400 // track = event->GetTrack(itrack);
401 // tzeroTrack = fESDpid->GetTOFResponse().GetStartTime(track->P());
402 // Printf("================> Track # %d mom: %f tzeroTrack %f",itrack,track->P(),tzeroTrack);
dc369c13 403 fESDpid->MakeTOFPID(event->GetTrack(itrack),0);
ee981ab3 404 }
405
406
407}
408
409
88ccd4cf 410//_____________________________________________________
411void 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//_____________________________________________________
429void 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;
ee981ab3 440
88ccd4cf 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}
dc369c13 459
460
fd611537 461//______________________________________________________________________________
462void 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");
ac93546e 481 return; // coverity
fd611537 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//______________________________________________________________________________
502void 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//______________________________________________________________________________
527void 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//_____________________________________________________
546void 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//________________________________________________________________________
581void 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//________________________________________________________________________
634void 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//________________________________________________________________________
649Double_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//________________________________________________________________________
775Double_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//________________________________________________________________________
901Double_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//________________________________________________________________________
1046Double_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//________________________________________________________________________
1090Double_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//________________________________________________________________________
1134Double_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