]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ANALYSIS/TenderSupplies/AliTOFTenderSupply.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
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");
481 }
482
483 //find pass from file name (UGLY, but not stored in ESD... )
484 TString fileName(file->GetName());
485 if (fileName.Contains("pass1") ) {
486 fRecoPass=1;
487 } else if (fileName.Contains("pass2") ) {
488 fRecoPass=2;
489 } else if (fileName.Contains("pass3") ) {
490 fRecoPass=3;
491 }
492 if (fRecoPass == 0) {
493 AliInfo(Form("From file name %s reco pass cannot be detected",fileName.Data()));
494 AliInfo("Change file name or use SetUserRecoPass method");
495 AliFatal("------------- TOF tender cannot run with reco pass unspecified, issuing FATAL error ---------- ");
496 }
497}
498
499
500//______________________________________________________________________________
501void AliTOFTenderSupply::InitGeom()
502{
503
504 if (fGeomSet == kTRUE) return;
505
506 // Printf("\n \n ----- calling InitGeom to fix TRD Bug ----- \n \n");
507 AliCDBManager * man = AliCDBManager::Instance();
508 man->SetDefaultStorage("raw://");
509 fCDBkey = man->SetLock(kFALSE, fCDBkey);
510 Int_t run = fTender->GetRun();
511 // Printf(" ---------> run is %d",run);
512 man->SetRun(run);
513 AliCDBEntry *entry = (AliCDBEntry*)man->Get("GRP/Geometry/Data");
514 if (entry) {
515 AliGeomManager::LoadGeometry();
516 AliGeomManager::ApplyAlignObjsFromCDB("ITS TPC TRD TOF");
517 // fCDBkey = man->SetLock(kTRUE, fCDBkey);
518 // Printf("\n \n ----- Geometry loaded ------ \n \n");
519 }
520 fGeomSet=kTRUE;
521
522}
523
524
525//______________________________________________________________________________
526void AliTOFTenderSupply::FixTRDBug(AliESDEvent* event)
527//
528// recompute texp fixing wrong dE/dx from TRD (LHC10b,c pass3)
529//
530{
531
532 if (fGeomSet == kFALSE) InitGeom();
533
534 // Printf("Running FixTRD bug ");
535 /* loop over tracks */
536 AliESDtrack *track = NULL;
537 for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) {
538 track = event->GetTrack(itrk);
539 FixTRDBug(track);
540 }
541}
542
543
544//_____________________________________________________
545void AliTOFTenderSupply::FixTRDBug(AliESDtrack *track)
546{
547 //
548 //
549 //
550
551
552 ULong_t status=track->GetStatus();
553 if (!( ( (status & AliVTrack::kITSrefit)==AliVTrack::kITSrefit ) &&
554 ( (status & AliVTrack::kTPCrefit)==AliVTrack::kTPCrefit ) &&
555 ( (status & AliVTrack::kTPCout)==AliVTrack::kTPCout ) &&
556 ( (status & AliVTrack::kTOFout)==AliVTrack::kTOFout ) &&
557 ( (status & AliVTrack::kTIME)==AliVTrack::kTIME ) ) ) return;
558
559 fIsEnteringInTRD=kFALSE;
560 fInTRD=kFALSE;
561 fIsComingOutTRD=kFALSE;
562 fOutTRD=kFALSE;
563
564 // Printf("Track reached TOF %f",track->P());
565 Double_t correctionTimes[5] = {0.,0.,0.,0.,0.}; // to be added to the expected times
566 FindTRDFix(track, correctionTimes);
567 Double_t expectedTimes[5] = {0.,0.,0.,0.,0.}; track->GetIntegratedTimes(expectedTimes);
568 // Printf("Exp. times: %f %f %f %f %f",
569 // expectedTimes[0],expectedTimes[1],expectedTimes[2],expectedTimes[3],expectedTimes[4]);
570 // Printf("Corr. times: %f %f %f %f %f",
571 // correctionTimes[0],correctionTimes[1],correctionTimes[2],correctionTimes[3],correctionTimes[4]);
572
573 for (Int_t jj=0; jj<5; jj++) expectedTimes[jj]+=correctionTimes[jj];
574 track->SetIntegratedTimes(expectedTimes);
575
576}
577
578
579//________________________________________________________________________
580void AliTOFTenderSupply::FindTRDFix(AliESDtrack *track,Double_t *corrections)
581{
582
583 Double_t pT = track->Pt();
584 ULong_t status=track->GetStatus();
585 Bool_t isTRDout = (status & AliVTrack::kTRDout)==AliVTrack::kTRDout;
586
587 Double_t length = 0.;
588
589 Double_t xyzIN[3]={0.,0.,0.};
590 fIsEnteringInTRD = track->GetXYZAt(fRhoTRDin,fMagField,xyzIN);
591
592 Double_t xyzOUT[3]={0.,0.,0.};
593 fIsComingOutTRD = track->GetXYZAt(fRhoTRDout,fMagField,xyzOUT);
594
595 if (fIsEnteringInTRD && fIsComingOutTRD) {
596
597
598 Double_t phiIN = TMath::Pi()+TMath::ATan2(-xyzIN[1],-xyzIN[0]);
599 phiIN *= TMath::RadToDeg();
600 fInTRD = ( (phiIN>= 0. && phiIN<= 40.) ||
601 (phiIN>=140. && phiIN<=220.) ||
602 (phiIN>=340. && phiIN<=360.) ); // TRD SMs installed @ 2010
603
604 Double_t phiOUT = TMath::Pi()+TMath::ATan2(-xyzOUT[1],-xyzOUT[0]);
605 phiOUT *= TMath::RadToDeg();
606 fOutTRD = ( (phiOUT>= 0. && phiOUT<= 40.) ||
607 (phiOUT>=140. && phiOUT<=220.) ||
608 (phiOUT>=340. && phiOUT<=360.) ); // TRD SMs installed @ 2010
609
610 length = 0.;
611
612 if (fInTRD || fOutTRD) {
613
614 if ( ( fInTRD && fOutTRD ) || ( fInTRD && !fOutTRD ) ) {
615 length = EstimateLengthInTRD1(track);
616 } else if ( !fInTRD && fOutTRD ) {
617 length = EstimateLengthInTRD2(track);
618 }
619
620 } else { // ( !fInTRD && !fOutTRD )
621
622 length = EstimateLengthOutTRD(track);
623
624 }
625
626 }
627 // Printf("estimated length in TRD %f [isTRDout %d]",length,isTRDout);
628 CorrectDeltaTimes(pT,length,isTRDout,corrections);
629
630}
631
632//________________________________________________________________________
633void AliTOFTenderSupply::CorrectDeltaTimes(Double_t pT,
634 Double_t length,
635 Bool_t flagTRDout,
636 Double_t *corrections)
637{
638
639 corrections[2] = CorrectExpectedPionTime(pT,length,flagTRDout);
640 corrections[0] = corrections[2]; // x electrons used pion corrections
641 corrections[1] = corrections[2]; // x muons used pion corrections
642 corrections[3] = CorrectExpectedKaonTime(pT,length,flagTRDout);
643 corrections[4] = CorrectExpectedProtonTime(pT,length,flagTRDout);
644
645}
646
647//________________________________________________________________________
648Double_t AliTOFTenderSupply::CorrectExpectedPionTime(Double_t pT,
649 Double_t length,
650 Bool_t isTRDout)
651{
652 // correction for expected time for pions
653
654 Double_t delta=0.;
655
656 // Printf("Flags Ent In ComingOut Out %d %d %d %d",fIsEnteringInTRD,fInTRD,fIsComingOutTRD,fOutTRD);
657 if (!fIsEnteringInTRD || !fIsComingOutTRD) { // zone 5
658
659 Float_t p[2]={0.,0.};
660
661 if (isTRDout) {
662
663 if (pT<0.30) {
664 p[0] = 180.; p[1] = 0.;
665 } else if (pT>=0.30 && pT<0.35) {
666 p[0] = 740.; p[1] = -1800.;
667 } else if (pT>=0.35 && pT<0.40) {
668 p[0] = 488.; p[1] =-1080.;
669 } else if (pT>=0.40 && pT<0.45) {
670 p[0] = 179.; p[1] = -307.;
671 } else if (pT>=0.45 && pT<0.50) {
672 p[0] = 97.; p[1] = -123.;
673 } else { //if (pT>=0.50)
674 p[0] = 120.; p[1] = -172.;
675 }
676
677 } else {
678
679 if (pT<0.30) {
680 p[0] = 70.; p[1] = 0.;
681 } else if (pT>=0.30 && pT<0.35) {
682 p[0] = 339.; p[1] = -927.;
683 } else if (pT>=0.35 && pT<0.40) {
684 p[0] = 59.; p[1] = -121.;
685 } else if (pT>=0.40 && pT<0.50) {
686 p[0] = 21.; p[1] = -24.;
687 } else { //if (pT>=0.50)
688 p[0] = 42.; p[1] = -67.;
689 }
690
691 }
692
693 delta = p[0]+p[1]*pT;
694 // Printf("Pion time: %f %f %f %f",p[0],p[1],length,delta);
695
696 } else {
697
698 Float_t p[2] = {0.,0.};
699
700 if ( fInTRD && fOutTRD) { // zone 1
701
702 if (isTRDout) {
703
704 if (length<130.) {
705 p[0] = 0.; p[1] = 0.;
706 } else if (length>=130. && length<170.) {
707 p[0] = -20.5; p[1] = 0.25;
708 } else {//if (length>=170.)
709 p[0] = 22.; p[1] = 0.;
710 }
711
712 } else { // !isTRDout
713
714 p[0] = 20.; p[1] = 0.;
715
716 }
717
718 } else if (!fInTRD && !fOutTRD) { // zone 2
719
720 p[0] = 0.; p[1] = 0.;
721
722 } else if ( fInTRD && !fOutTRD) { // zone 3
723
724 if (isTRDout) {
725
726 if (length< 75.) {
727 p[0] = 17.; p[1] = 0.;
728 } else if (length>= 75. && length< 95.) {
729 p[0] = 81.; p[1] = -0.85;
730 } else if (length>= 95. && length<155.) {
731 p[0] = 0.; p[1] = 0.;
732 } else {//if (length>=155.)
733 p[0] = 10.; p[1] = 0.;
734 }
735
736 } else { // !isTRDout
737
738 p[0] = 0.; p[1] = 0.;
739
740 }
741
742 } else if (!fInTRD && fOutTRD) { // zone 4
743
744 if (isTRDout) {
745
746 if (length<80.) {
747 p[0] = 0.; p[1] = 0.;
748 } else {//if (length>=80.)
749 p[0] = 10.; p[1] = 0.;
750 }
751
752 } else { // !isTRDout
753
754 if (length<30.) {
755 p[0] = 0.; p[1] = 0.;
756 } else {//if (length>=30.)
757 p[0] = 6.; p[1] = 0.;
758 }
759
760 }
761
762 }
763
764 delta = p[0]+p[1]*length;
765 // Printf("Pion time: %f %f %f %f",p[0],p[1],length,delta);
766
767 }
768
769 return delta;
770
771}
772
773//________________________________________________________________________
774Double_t AliTOFTenderSupply::CorrectExpectedKaonTime(Double_t pT,
775 Double_t length,
776 Bool_t isTRDout)
777{
778 // correction for expected time for kaons
779
780 Double_t delta=0.;
781 // Printf("Flags Ent In ComingOut Out %d %d %d %d",fIsEnteringInTRD,fInTRD,fIsComingOutTRD,fOutTRD);
782
783 if (!fIsEnteringInTRD || !fIsComingOutTRD) { // zone 5
784
785 Float_t p[2]={0.,0.};
786
787 if (isTRDout) {
788
789 if (pT<0.4) {
790 p[0] = 900.; p[1] = 0.;
791 } else if (pT>=0.40 && pT<0.45) {
792 p[0] = 3100.; p[1] = -6000.;
793 } else if (pT>=0.45 && pT<0.50) {
794 p[0] = 1660.; p[1] = -2800.;
795 } else if (pT>=0.50 && pT<0.55) {
796 p[0] = 860.; p[1] = -1200.;
797 } else { //if (pT>=0.55)
798 p[0] = 200.; p[1] = 0.;
799 }
800
801 } else {
802
803 if (pT<0.30) {
804 p[0] = 0.; p[1] = 0.;
805 } else if (pT>=0.30 && pT<0.32) {
806 p[0] = 570.; p[1] = 0.;
807 } else if (pT>=0.32 && pT<0.35) {
808 p[0] = 3171.; p[1] = -8133.;
809 } else if (pT>=0.35 && pT<0.40) {
810 p[0] = 1815.; p[1] = -4260.;
811 } else if (pT>=0.40 && pT<0.45) {
812 p[0] = 715.; p[1] = -1471.;
813 } else if (pT>=0.45 && pT<0.50) {
814 p[0] = 233.; p[1] = -407.;
815 } else if (pT>=0.50 && pT<0.55) {
816 p[0] = 408.; p[1] = -752.;
817 } else { //if (pT>=0.55)
818 p[0] = 408.-752.*0.55; p[1] = 0.;
819 }
820
821 }
822
823 delta = p[0]+p[1]*pT;
824 // Printf("Kaon time: %f %f %f %f",p[0],p[1],length,delta);
825
826 } else {
827
828 Float_t p[2] = {0.,0.};
829
830 if ( fInTRD && fOutTRD) { // zone 1
831
832 if (isTRDout) {
833
834 if (length<95.) {
835 p[0] = 20.; p[1] = 0.;
836 } else if (length>=95. && length<195.) {
837 p[0] = -24.0; p[1] = 0.10+0.0041*length;
838 } else {//if (length>=195.)
839 p[0] = 150.; p[1] = 0.;
840 }
841
842 } else { // !isTRDout
843
844 p[0] = 40.; p[1] = 0.;
845
846 }
847
848 } else if (!fInTRD && !fOutTRD) { // zone 2
849
850 p[0] = 0.; p[1] = 0.;
851
852 } else if ( fInTRD && !fOutTRD) { // zone 3
853
854 if (isTRDout) {
855
856 if (length< 15.) {
857 p[0] = 180.; p[1] = 0.;
858 } else if (length>= 15. && length< 55.) {
859 p[0] = 215.; p[1] = -2.5;
860 } else {//if (length>=55.)
861 p[0] = 78.; p[1] = 0.;
862 }
863
864 } else { // !isTRDout
865
866 p[0] = 0.; p[1] = 0.;
867
868 }
869
870 } else if (!fInTRD && fOutTRD) { // zone 4
871
872 if (isTRDout) {
873
874 if (length< 55.) {
875 p[0] = 0.; p[1] = 0.;
876 } else if (length>= 55. && length<115.) {
877 p[0] = -85.; p[1] = 1.9;
878 } else {//if (length>=115.)
879 p[0] = 100.; p[1] = 0.;
880 }
881
882 } else { // !isTRDout
883
884 p[0] = 0.; p[1] = 0.;
885
886 }
887
888 }
889
890 delta = p[0]+p[1]*length;
891 // Printf("Kaon time: %f %f %f %f",p[0],p[1],length,delta);
892
893 }
894
895 return delta;
896
897}
898
899//________________________________________________________________________
900Double_t AliTOFTenderSupply::CorrectExpectedProtonTime(Double_t pT,
901 Double_t length,
902 Bool_t isTRDout)
903{
904 // correction for expected time for protons
905
906 Double_t delta=0.;
907 // Printf("Flags Ent In ComingOut Out %d %d %d %d",fIsEnteringInTRD,fInTRD,fIsComingOutTRD,fOutTRD);
908
909 if (!fIsEnteringInTRD || !fIsComingOutTRD) { // zone 5
910 Float_t p[2]={0.,0.};
911
912
913 if (isTRDout) {
914
915 if (pT<0.375) {
916 p[0] = 1000.; p[1] = 0.;
917 } else if (pT>=0.375 && pT<0.45) {
918 p[0] = 1500.; p[1] = 0.;
919 } else if (pT>=0.45 && pT<0.50) {
920 p[0] = 4650.; p[1] = -7000.;
921 } else if (pT>=0.50 && pT<0.55) {
922 p[0] = 3150.; p[1] = -4000.;
923 } else { //if (pT>=0.55)
924 p[0] = 3150. -4000.*0.55; p[1] = 0.;
925 }
926
927 } else {
928
929 if (pT<0.32) {
930 p[0] = 2963.-5670.*0.032; p[1] = 0.;
931 } else if (pT>=0.32 && pT<0.35) {
932 p[0] = 2963.; p[1] = -5670.;
933 } else if (pT>=0.35 && pT<0.40) {
934 p[0] = 4270.; p[1] = -9400.;
935 } else if (pT>=0.40 && pT<0.45) {
936 p[0] = 1550.; p[1] = -2600.;
937 } else if (pT>=0.45 && pT<0.50) {
938 p[0] = 1946.; p[1] = -3480.;
939 } else if (pT>=0.50 && pT<0.55) {
940 p[0] = 1193.; p[1] = -1974.;
941 } else { //if (pT>=0.55)
942 p[0] = 1193.-1974.*0.55; p[1] = 0.;
943 }
944
945 }
946
947 delta = p[0]+p[1]*pT;
948 // Printf("Proton time: %f %f %f %f",p[0],p[1],length,delta);
949
950 } else {
951
952 Float_t p[2] = {0.,0.};
953
954 if ( fInTRD && fOutTRD) { // zone 1
955
956 if (isTRDout) {
957
958 if (length<90.) {
959 p[0] = 0.; p[1] = 0.;
960 } else if (length>=90. && length<200.) {
961 p[0] = 1063.; p[1] = -32.+0.30*length-0.00072*length*length;
962 } else {//if (length>=200.)
963 p[0] = 900.; p[1] = 0.;
964 }
965
966 } else { // !isTRDout
967
968 p[0] = 80.; p[1] = 0.;
969
970 }
971
972 } else if (!fInTRD && !fOutTRD) { // zone 2
973
974 if (isTRDout) {
975 p[0] = 0.; p[1] = 0.;
976 } else {
977 if (length<125.) {
978 p[0] = 0.; p[1] = 0.;
979 } else if (length>=125. && length<180.) {
980 p[0] = -132.; p[1] = 1.3;
981 } else {
982 p[0] = 100.; p[1] = 0.;
983 }
984
985 }
986
987 } else if ( fInTRD && !fOutTRD) { // zone 3
988
989 if (isTRDout) {
990
991 if (length< 30.) {
992 p[0] = 670.; p[1] = 0.;
993 } else if (length>= 30. && length<155.) {
994 p[0] = 944.; p[1] = -11.+0.064*length;
995 } else {//if (length>=155.)
996 p[0] = 780.; p[1] = 0.;
997 }
998
999 } else { // !isTRDout
1000
1001 if (length< 30.) {
1002 p[0] = 140.; p[1] = -4.5;
1003 } else {
1004 p[0] = 0.; p[1] = 0.;
1005 }
1006
1007 }
1008
1009 } else if (!fInTRD && fOutTRD) { // zone 4
1010
1011 if (isTRDout) {
1012
1013 if (length< 45.) {
1014 p[0] = 130.; p[1] = 0.;
1015 } else if (length>= 45. && length<120.) {
1016 p[0] = -190.; p[1] = 6.5;
1017 } else {//if (length>=120.)
1018 p[0] = 750.; p[1] = 0.;
1019 }
1020
1021 } else { // !isTRDout
1022
1023 if (length<75.5) {
1024 p[0] = 0.; p[1] = 0.;
1025 } else if (length>= 75.5 && length<90.) {
1026 p[0] = -830.; p[1] = 11.;
1027 } else {
1028 p[0] = 160.; p[1] = 0.;
1029 }
1030
1031 }
1032
1033 }
1034
1035 delta = p[0]+p[1]*length;
1036 // Printf("Proton time: %f %f %f %f",p[0],p[1],length,delta);
1037
1038 }
1039
1040 return delta;
1041
1042}
1043
1044//________________________________________________________________________
1045Double_t AliTOFTenderSupply::EstimateLengthInTRD1(AliESDtrack *track)
1046{
1047
1048 Double_t xyz0[3]={0.,0.,0.};
1049 Bool_t stayInTRD = track->GetXYZAt(fRhoTRDin,fMagField,xyz0);
1050
1051 Double_t phi0 = TMath::Pi()+TMath::ATan2(-xyz0[1],-xyz0[0]);
1052 phi0 *= TMath::RadToDeg();
1053 stayInTRD = stayInTRD && ( (phi0>= 0. && phi0<= 40.) ||
1054 (phi0>=140. && phi0<=220.) ||
1055 (phi0>=340. && phi0<=360.) );
1056
1057 Double_t trackLengthInTRD = 0.;
1058 Int_t iStep=0;
1059
1060 Double_t b[3];track->GetBxByBz(b);
1061
1062 Double_t xyz1[3]={0.,0.,0.};
1063 Double_t rho = fRhoTRDin;
1064 while (stayInTRD && rho<=fRhoTRDout) {
1065 iStep++;
1066 rho += fStep;
1067
1068 for (Int_t ii=0; ii<3; ii++) xyz1[ii]=0.;
1069 stayInTRD = track->GetXYZAt(rho,fMagField,xyz1);
1070 Double_t phi1 = TMath::Pi()+TMath::ATan2(-xyz1[1],-xyz1[0]);
1071 phi1 *= TMath::RadToDeg();
1072 stayInTRD = stayInTRD && ( (phi1>= 0. && phi1<= 40.) ||
1073 (phi1>=140. && phi1<=220.) ||
1074 (phi1>=340. && phi1<=360.) );
1075
1076 Double_t l2 = TMath::Sqrt((xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0]) +
1077 (xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1]) +
1078 (xyz1[2]-xyz0[2])*(xyz1[2]-xyz0[2]));
1079 trackLengthInTRD += l2;
1080
1081 for (Int_t ii=0; ii<3; ii++) xyz0[ii]=xyz1[ii];
1082 }
1083
1084 return trackLengthInTRD;
1085
1086}
1087
1088//________________________________________________________________________
1089Double_t AliTOFTenderSupply::EstimateLengthInTRD2(AliESDtrack *track)
1090{
1091
1092 Double_t xyz0[3]={0.,0.,0.};
1093 Bool_t stayInTRD = track->GetXYZAt(fRhoTRDout,fMagField,xyz0);
1094
1095 Double_t phi0 = TMath::Pi()+TMath::ATan2(-xyz0[1],-xyz0[0]);
1096 phi0 *= TMath::RadToDeg();
1097 stayInTRD = stayInTRD && ( (phi0>= 0. && phi0<= 40.) ||
1098 (phi0>=140. && phi0<=220.) ||
1099 (phi0>=340. && phi0<=360.) );
1100
1101 Double_t trackLengthInTRD = 0.;
1102 Int_t iStep=0;
1103
1104 Double_t b[3];track->GetBxByBz(b);
1105
1106 Double_t xyz1[3]={0.,0.,0.};
1107 Double_t rho = fRhoTRDout;
1108 while (stayInTRD && rho>=fRhoTRDin) {
1109 iStep++;
1110 rho -= fStep;
1111
1112 for (Int_t ii=0; ii<3; ii++) xyz1[ii]=0.;
1113 stayInTRD = track->GetXYZAt(rho,fMagField,xyz1);
1114 Double_t phi1 = TMath::Pi()+TMath::ATan2(-xyz1[1],-xyz1[0]);
1115 phi1 *= TMath::RadToDeg();
1116 stayInTRD = stayInTRD && ( (phi1>= 0. && phi1<= 40.) ||
1117 (phi1>=140. && phi1<=220.) ||
1118 (phi1>=340. && phi1<=360.) );
1119
1120 Double_t l2 = TMath::Sqrt((xyz0[0]-xyz1[0])*(xyz0[0]-xyz1[0]) +
1121 (xyz0[1]-xyz1[1])*(xyz0[1]-xyz1[1]) +
1122 (xyz0[2]-xyz1[2])*(xyz0[2]-xyz1[2]));
1123 trackLengthInTRD += l2;
1124
1125 for (Int_t ii=0; ii<3; ii++) xyz0[ii]=xyz1[ii];
1126 }
1127
1128 return trackLengthInTRD;
1129
1130}
1131
1132//________________________________________________________________________
1133Double_t AliTOFTenderSupply::EstimateLengthOutTRD(AliESDtrack *track)
1134{
1135
1136 Double_t xyz0[3]={0.,0.,0.};
1137 Bool_t stayInTRD = track->GetXYZAt(fRhoTRDin,fMagField,xyz0);
1138
1139 Double_t phi0 = TMath::Pi()+TMath::ATan2(-xyz0[1],-xyz0[0]);
1140 phi0 *= TMath::RadToDeg();
1141 stayInTRD = stayInTRD && !( (phi0>= 0. && phi0<= 40.) ||
1142 (phi0>=140. && phi0<=220.) ||
1143 (phi0>=340. && phi0<=360.) );
1144
1145 Double_t trackLengthInTRD = 0.;
1146 Int_t iStep=0;
1147
1148 Double_t b[3];track->GetBxByBz(b);
1149
1150 Double_t xyz1[3]={0.,0.,0.};
1151 Double_t rho = fRhoTRDin;
1152 while (stayInTRD && rho<=fRhoTRDout) {
1153 iStep++;
1154 rho += fStep;
1155
1156 for (Int_t ii=0; ii<3; ii++) xyz1[ii]=0.;
1157 stayInTRD = track->GetXYZAt(rho,fMagField,xyz1);
1158 Double_t phi1 = TMath::Pi()+TMath::ATan2(-xyz1[1],-xyz1[0]);
1159 phi1 *= TMath::RadToDeg();
1160 stayInTRD = stayInTRD && !( (phi1>= 0. && phi1<= 40.) ||
1161 (phi1>=140. && phi1<=220.) ||
1162 (phi1>=340. && phi1<=360.) );
1163
1164 Double_t l2 = TMath::Sqrt((xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0]) +
1165 (xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1]) +
1166 (xyz1[2]-xyz0[2])*(xyz1[2]-xyz0[2]));
1167 trackLengthInTRD += l2;
1168
1169 for (Int_t ii=0; ii<3; ii++) xyz0[ii]=xyz1[ii];
1170 }
1171
1172 return trackLengthInTRD;
1173
1174}
1175