Test for Coverity
[u/mrichter/AliRoot.git] / TPC / AliTPCTrackHitsV2.cxx
CommitLineData
f641f6bd 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
88cb7938 16/* $Id$ */
19364939 17
f641f6bd 18///////////////////////////////////////////////////////////////////////////////
19// //
20// Time Projection Chamber track hits object //
21//
22// Origin: Marian Ivanov , GSI Darmstadt
23//
24// AliTPCTrackHitsV2
25// Container for Track Hits - based on standard TClonesArray -
26// fArray of AliTPCTrackHitsParamV2
27// In AliTPCTrackHitsParamV2 - parameterization of the track segment is stored
28// for each of the track segment - relative position ( distance between hits) and
29// charge of the hits is stored - comparing to classical TClonesArray of AliTPChit -
30// comperssion factor of 5-7 (depending on the required precision) -
31// In future release AliTPCTrackHitsV2 - will replace old AliTPCTrackHits - which were not
32// based on standard ROOT containers
33// Basic function:
34// // during building Container
35// AddHitKartez(Int_t volumeID, Int_t trackID, Double_t x, Double_t y, Double_t z,Int_t q)
36// void SetHitPrecision(Double_t prec) {fPrecision=prec;}
37// void SetStepPrecision(Double_t prec) {fStep=prec;}
38// Bool_t FlushHitStack(Bool_t force=kTRUE);
39// //at the end necessary to have Container in consistent state
40//
41// // looping over Container
42// Bool_t First(), Bool_t Next() - iterators - return status of the operation
43// AliTPChit * GetHit(); - return current hit
44
45
46//Begin_Html
47/*
48<img src="gif/AliTPCTrackHitsV2.gif">
49*/
50//End_Html
51// //
52// //
53///////////////////////////////////////////////////////////////////////////////
018a927a 54//
f641f6bd 55
3010c308 56#include <TClonesArray.h>
57#include <TMath.h>
58
f641f6bd 59#include "AliTPCTrackHitsV2.h"
f641f6bd 60#include "AliTPC.h"
61
f641f6bd 62
63
64ClassImp(AliTPCTrackHitsV2)
65ClassImp(AliTrackHitsParamV2)
66
67 //
68Int_t AliTrackHitsParamV2::fgCounter1 =0;
69Int_t AliTrackHitsParamV2::fgCounter2 =0;
70//
71Int_t AliTPCTrackHitsV2::fgCounter1 =0;
72Int_t AliTPCTrackHitsV2::fgCounter2 =0;
73//
74const Double_t AliTPCTrackHitsV2::fgkPrecision=1e-6; //precision
75const Double_t AliTPCTrackHitsV2::fgkPrecision2=1e-20; //precision
e61fd20d 76const Double_t AliTPCTrackHitsV2::fgkTimePrecision=20.e-9; //hit time precision
f641f6bd 77
78
79
80
b9671574 81class AliTPCTempHitInfoV2 {
82public:
f641f6bd 83 AliTPCTempHitInfoV2();
e61fd20d 84 void SetHit(Double_t r, Double_t z, Double_t fi, Int_t q, Float_t time);
b9671574 85 UInt_t GetStackIndex() const {return fStackIndex;}
86 void SetStackIndex(UInt_t i) {fStackIndex=i;}
87 UInt_t GetParamIndex() const {return fParamIndex;}
88 void SetParamIndex(UInt_t i) {fParamIndex=i;}
89 Float_t GetTimeStack(Int_t i) const {return fTimeStack[i];}
90 const Float_t* GetTimeStackP(Int_t i) const {return &fTimeStack[i];}
91 UInt_t GetQStack(Int_t i) const {return fQStack[i];}
92 const UInt_t* GetQStackP(Int_t i) const {return &fQStack[i];}
f641f6bd 93 Double_t * GetPosition(Int_t index){return &fPositionStack[index*3];}
b9671574 94 Double_t GetOldR() const {return fOldR;}
95 void SetOldR(Double_t r) {fOldR=r;}
96
97
98 AliTrackHitsParamV2 * GetParam() const {return fParam;}
99 void SetParam(AliTrackHitsParamV2 * p) {fParam=p;}
100 void UpdateParam(Double_t maxdelta); //recal
101 void NewParam(Double_t r, Double_t z, Double_t fi, Int_t q, Float_t time);
102 enum {kStackSize = 10000};
103
104protected:
105 AliTPCTempHitInfoV2(const AliTPCTempHitInfoV2 &hit);
106 AliTPCTempHitInfoV2& operator = (const AliTPCTempHitInfoV2 &hit);
f641f6bd 107 void Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2,
108 Double_t fSumX, Double_t fSumX2, Double_t fSumX3,
109 Double_t fSumX4, Int_t n,
110 Double_t &a, Double_t &b, Double_t &c);
111 void Fit(AliTrackHitsParamV2 * param);
53110b9d 112 Double_t fSumDr; // Sum of Dr
113 Double_t fSumDr2; // Square of sum of Dr
114 Double_t fSumDr3; // Cube of sum of Dr
115 Double_t fSumDr4; // Fourth power of sum of Dr
116 Double_t fSumDFi; // Sum of DFi
117 Double_t fSumDFiDr; // Sum of DFiDr
118 Double_t fSumDFiDr2;// Sum of square of DFiDr
119 Double_t fSumDZ; // Sum of DZ
120 Double_t fSumDZDr; // Sum of DZDr
121 Double_t fSumDZDr2; // Sum of square of DZDr
f641f6bd 122 Double_t fOldR; //previos r
53110b9d 123 Double_t fPositionStack[3*kStackSize]; //position stack
124 UInt_t fQStack[kStackSize]; //Q stack
e61fd20d 125 Float_t fTimeStack[kStackSize]; //time stack
f641f6bd 126 UInt_t fStackIndex; //current stack index
127 // UInt_t fInfoIndex; //current track info index
128 UInt_t fParamIndex; //current track parameters index
129 // AliTrackHitsInfo * fInfo; //current track info
130 AliTrackHitsParamV2 * fParam; //current track param
131};
132
133
134AliTPCTempHitInfoV2::AliTPCTempHitInfoV2()
179c6296 135 :fSumDr(0.),
136 fSumDr2(0.),
137 fSumDr3(0.),
138 fSumDr4(0.),
139 fSumDFi(0.),
140 fSumDFiDr(0.),
141 fSumDFiDr2(0.),
142 fSumDZ(0.),
143 fSumDZDr(0.),
144 fSumDZDr2(0.),
145 fOldR(0.),
146 fStackIndex(0),
147 fParamIndex(0),
148 fParam(0)
f641f6bd 149{
150 //
53110b9d 151 // Standard constructor
152 // set to default value
153 //
f641f6bd 154 fSumDr=fSumDr2=fSumDr3=fSumDr4=
155 fSumDFi=fSumDFiDr=fSumDFiDr2=
156 fSumDZ=fSumDZDr=fSumDZDr2=0;
157 fStackIndex = 0;
158 // fInfoIndex = 0;
159 fParamIndex = 0;
712976a6 160 for(Int_t i=0;i<3*kStackSize;i++) fPositionStack[i]=0.;
161 for(Int_t i=0;i<kStackSize;i++){
162 fQStack[i]=0;
163 fTimeStack[i]=0.;
164 }
f641f6bd 165}
166
167
e61fd20d 168void AliTPCTempHitInfoV2::NewParam(Double_t r, Double_t z, Double_t fi, Int_t q, Float_t time)
f641f6bd 169{
170 //
171 //reset stack and sum parameters
172 //store line initial point
53110b9d 173 //
f641f6bd 174 fSumDr=fSumDr2=fSumDr3=fSumDr4=
175 fSumDFi=fSumDFiDr=fSumDFiDr2=
176 fSumDZ=fSumDZDr=fSumDZDr2=0;
177 fStackIndex=0;
b9671574 178 fParam->SetR(r);
f641f6bd 179 fOldR = r;
b9671574 180 fParam->SetZ(z);
181 fParam->SetFi(fi);
182 fParam->SetAn(0.);
183 fParam->SetAd(0.);
184 fParam->SetTheta(0.);
185 fParam->SetThetaD(0.);
e61fd20d 186 SetHit(r,z,fi,q,time);
f641f6bd 187}
188
e61fd20d 189void AliTPCTempHitInfoV2::SetHit(Double_t r, Double_t z, Double_t fi, Int_t q, Float_t time)
f641f6bd 190{
191 //
192 //add hit to the stack
193 //recalculate new estimete of line parameters
194 Double_t *f = GetPosition(fStackIndex);
195 f[0] = r;
196 f[1] = z;
197 f[2] = fi;
198 fQStack[fStackIndex]=q;
e61fd20d 199 fTimeStack[fStackIndex]=time;
f641f6bd 200 if (fStackIndex==0) return;
b9671574 201 Double_t dr = (r-fParam->GetR());
202 if (TMath::Abs(dr)<AliTPCTrackHitsV2::GetKPrecision()) dr =AliTPCTrackHitsV2::GetKPrecision();
203 Double_t dfi = fi-fParam->GetFi();
204 Double_t dz = z -fParam->GetZ();
f641f6bd 205 Double_t dr2 =dr*dr;
206 Double_t dr3 =dr2*dr;
207 Double_t dr4 =dr3*dr;
208 fSumDr +=dr;
209 fSumDr2+=dr2;
210 fSumDr3+=dr3;
211 fSumDr4+=dr4;
212 fSumDFi +=dfi;
213 fSumDFiDr+=dfi*dr;
214 fSumDFiDr2+=dfi*dr2;
215 fSumDZ +=dz;
216 fSumDZDr+=dz*dr;
217 fSumDZDr2+=dz*dr2;
218
219 //update fit parameters
220 //
221 Double_t det = fSumDr2*fSumDr4-fSumDr3*fSumDr3;
b9671574 222 if (TMath::Abs(det)<AliTPCTrackHitsV2::GetKPrecision2()) return;
f641f6bd 223 if ( ( fStackIndex>1 ) ){
b9671574 224 fParam->SetAn((fSumDr4*fSumDFiDr-fSumDr3*fSumDFiDr2)/det);
225 fParam->SetAd((fSumDr2*fSumDFiDr2-fSumDr3*fSumDFiDr)/det);
f641f6bd 226 }
227 else
b9671574 228 fParam->SetAn(fSumDFiDr/fSumDr2);
f641f6bd 229 if ( ( fStackIndex>1 ) ){
b9671574 230 fParam->SetTheta((fSumDr4*fSumDZDr-fSumDr3*fSumDZDr2)/det);
231 fParam->SetThetaD((fSumDr2*fSumDZDr2-fSumDr3*fSumDZDr)/det);
f641f6bd 232 }
233 else
b9671574 234 fParam->SetTheta(fSumDZDr/fSumDr2);
f641f6bd 235}
236
237
238void AliTPCTempHitInfoV2::UpdateParam(Double_t maxdelta)
239{
53110b9d 240 //
241 // recalc parameters not fixing origin point
242 //
f641f6bd 243 if (fStackIndex>5){
244 Double_t a,b,c;
245 a=b=c=0;
246 Fit2(fSumDFi, fSumDFiDr, fSumDFiDr2, fSumDr,fSumDr2,fSumDr3,fSumDr4,
247 fStackIndex, a,b,c);
248 if (TMath::Abs(a)<maxdelta){
b9671574 249 fParam->SetFi(fParam->GetFi()+a/fParam->GetR());
250 fParam->SetAn(b);
251 fParam->SetAd(c);
f641f6bd 252 }
253 Fit2(fSumDZ, fSumDZDr, fSumDZDr2, fSumDr,fSumDr2,fSumDr3,fSumDr4,
254 fStackIndex, a,b,c) ;
255 if (TMath::Abs(a)<maxdelta){
b9671574 256 fParam->SetZ(fParam->GetZ()+a);
257 fParam->SetTheta(b);
258 fParam->SetThetaD(c);
f641f6bd 259 }
260 }
261
262}
53110b9d 263
f641f6bd 264void AliTPCTempHitInfoV2::Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2,
265 Double_t fSumX, Double_t fSumX2, Double_t fSumX3,
266 Double_t fSumX4, Int_t n,
267 Double_t &a, Double_t &b, Double_t &c)
268{
53110b9d 269 //
270 // fit of second order
271 //
f641f6bd 272 Double_t det =
273 n* (fSumX2*fSumX4-fSumX3*fSumX3) -
274 fSumX* (fSumX*fSumX4-fSumX3*fSumX2)+
275 fSumX2* (fSumX*fSumX3-fSumX2*fSumX2);
276
b9671574 277 if (TMath::Abs(det)> AliTPCTrackHitsV2::GetKPrecision()) {
f641f6bd 278 a =
279 (fSumY * (fSumX2*fSumX4-fSumX3*fSumX3)-
280 fSumX *(fSumYX*fSumX4-fSumYX2*fSumX3)+
281 fSumX2*(fSumYX*fSumX3-fSumYX2*fSumX2))/det;
282 b=
283 (n*(fSumYX*fSumX4-fSumX3*fSumYX2)-
284 fSumY*(fSumX*fSumX4-fSumX3*fSumX2)+
285 fSumX2*(fSumX*fSumYX2-fSumYX*fSumX2))/det;
286 c=
287 (n*(fSumX2*fSumYX2-fSumYX*fSumX3)-
288 fSumX*(fSumX*fSumYX2-fSumYX*fSumX2)+
289 fSumY*(fSumX*fSumX3-fSumX2*fSumX2))/det;
290 }
291}
292
293void AliTPCTempHitInfoV2::Fit(AliTrackHitsParamV2 * param)
294{
53110b9d 295 //
f641f6bd 296 // fit fixing first and the last point
53110b9d 297 // result stored in new param
298 //
b9671574 299 Double_t dx2 = (GetPosition(fStackIndex))[0]-fParam->GetR();
f641f6bd 300 Double_t det = fSumDr4+dx2*fSumDr2-2*dx2*fSumDr3;
b9671574 301 if ( (TMath::Abs(det)> AliTPCTrackHitsV2::GetKPrecision()) &&
302 ((TMath::Abs(dx2)> AliTPCTrackHitsV2::GetKPrecision()))){
303 Double_t dfi2 = (GetPosition(fStackIndex))[1]-fParam->GetFi();
304 param->SetAd((fSumDFiDr2+dfi2*fSumDr-dx2*fSumDFiDr-dfi2*fSumDr3/dx2)/det);
305 param->SetAn((dfi2-param->GetAd()*dx2*dx2)/dx2);
f641f6bd 306
b9671574 307 Double_t dz2 = (GetPosition(fStackIndex))[1]-fParam->GetZ();
308 param->SetTheta((fSumDZDr2+dz2*fSumDr-dx2*fSumDZDr-dz2*fSumDr3/dx2)/det);
309 param->SetTheta((dz2-param->GetAd()*dx2*dx2)/dx2);
f641f6bd 310 }
311
312}
313
179c6296 314AliTrackHitsParamV2::AliTrackHitsParamV2():TObject(),
315 fTrackID(0),
316 fVolumeID(0),
317 fR(0.),
318 fZ(0.),
319 fFi(0.),
320 fAn(0.),
321 fAd(0.),
322 fTheta(0.),
323 fThetaD(0.),
324 fNHits(0),
325 fHitDistance(0),
326 fCharge(0),
327 fTime(0)
f641f6bd 328{
53110b9d 329 //
330 // default constructor
331 //
f641f6bd 332 fgCounter1++;
333 fgCounter2++;
f641f6bd 334}
53110b9d 335
f641f6bd 336AliTrackHitsParamV2::~AliTrackHitsParamV2()
337{
53110b9d 338 //
339 // Standard destructor
340 //
f641f6bd 341 fgCounter1--;
342 if (fHitDistance) {
343 delete[]fHitDistance;
344 fHitDistance=0;
345 }
346 if (fCharge){
347 delete[]fCharge;
348 fCharge =0;
349 }
e61fd20d 350 if (fTime){
351 delete[]fTime;
352 fTime =0;
353 }
f641f6bd 354}
355
af05e6c8 356Float_t AliTrackHitsParamV2::Eta() const
357{
358 Float_t ctg = fZ / fR;
359 Float_t eta = -TMath::Log(TMath::Hypot(1,ctg)-TMath::Abs(ctg));
360 if(ctg < 0) eta = -eta;
361 return eta;
362}
363
f641f6bd 364
179c6296 365AliTPCTrackHitsV2::AliTPCTrackHitsV2():TObject(),
366 fArray(0),
367 fSize(0),
368 fPrecision(0.),
369 fStep(0.),
370 fMaxDistance(0),
371 fNVolumes(0),
372 fVolumes(0),
373 fTempInfo(0),
374 fCurrentHit(0),
375 fHit(0)
f641f6bd 376{
377 //
378 //default constructor
379 //
380 const Float_t kHitPrecision=0.002; //default precision for hit position in cm
381 const Float_t kStep =0.003; //30 mum step
382 const UShort_t kMaxDistance =100; //maximum distance 100
383
384 fPrecision=kHitPrecision; //precision in cm
385 fStep = kStep; //step size
386 fMaxDistance = kMaxDistance; //maximum distance
179c6296 387
f641f6bd 388 //fTrackHitsInfo = new AliObjectArray("AliTrackHitsInfo");
389 //fTrackHitsParam = new AliObjectArray("AliTrackHitsParamV2");
390 //fHitsPosAndQ = new TArrayOfArrayVStack("AliHitInfo");
391 fArray = new TClonesArray("AliTrackHitsParamV2");
392 fCurrentHit = new AliTPCCurrentHitV2;
f641f6bd 393 fgCounter1++;
394 fgCounter2++;
395
396}
397
398AliTPCTrackHitsV2::~AliTPCTrackHitsV2()
399{
400 //
401 //default destructor
402 //
403 // if (fTrackHitsInfo) delete fTrackHitsInfo;
404 if (fArray) {
405 delete fArray;
406 fArray =0;
407 }
408 //if (fHitsPosAndQ) delete fHitsPosAndQ;
409 if (fCurrentHit) delete fCurrentHit;
410 if (fTempInfo) delete fTempInfo;
411 if (fVolumes) {
412 delete [] fVolumes;
413 fVolumes =0;
414 fNVolumes=0;
415 }
fce0283f 416 if (fHit){
417 delete fHit;
418 fHit=0;
419 }
f641f6bd 420 fgCounter1--;
421}
422
5f6e636c 423void AliTPCTrackHitsV2::Clear(Option_t * /*option*/)
f641f6bd 424{
425 //
53110b9d 426 // clear object
427 //
6e493dba 428 fSize = 0;
f641f6bd 429 if (fArray){
430 for (Int_t i=0;i<fArray->GetEntriesFast();i++){
431 AliTrackHitsParamV2 * par = (AliTrackHitsParamV2 *)fArray->UncheckedAt(i);
432 par->~AliTrackHitsParamV2(); // delete object
433 }
434 fArray->Clear();
435 }
436 if (fTempInfo){
437 delete fTempInfo;
fce0283f 438 delete fHit;
439 fHit =0;
f641f6bd 440 fTempInfo =0;
441 }
442 if (fVolumes){
443 delete [] fVolumes;
444 fVolumes=0;
445 fNVolumes=0;
446 }
447}
448
449
450void AliTPCTrackHitsV2::AddHitKartez(Int_t volumeID, Int_t trackID, Double_t x,
e61fd20d 451 Double_t y, Double_t z,Int_t q, Float_t time)
f641f6bd 452{
453 //
53110b9d 454 // add hit to the container - it add hit at the end - input in global coordinata
455 //
f641f6bd 456 Double_t r = TMath::Sqrt(x*x+y*y);
457 Double_t fi = TMath::ACos(x/r);
458 if (y<0) fi*=-1.;
e61fd20d 459 AddHit(volumeID,trackID,r,z,fi,q,time);
f641f6bd 460}
461
462
463void AliTPCTrackHitsV2::AddHit(Int_t volumeID, Int_t trackID,
e61fd20d 464 Double_t r, Double_t z, Double_t fi, Int_t q, Float_t time)
f641f6bd 465{
466 //
53110b9d 467 // Adding one hit
468 //
f641f6bd 469 fSize++;
470 Bool_t diff=kFALSE;
471 if (!fTempInfo) { //initialisation of track - initialisation of parameters
472 fTempInfo = new AliTPCTempHitInfoV2;
b9671574 473 fTempInfo->SetParam(new((*fArray)[0]) AliTrackHitsParamV2);
474 fTempInfo->GetParam()->SetVolumeID(volumeID);
475 fTempInfo->GetParam()->SetTrackID(trackID);
f641f6bd 476 AddVolume(volumeID);
477 //
b9671574 478 fTempInfo->SetParamIndex(0);
e61fd20d 479 fTempInfo->NewParam(r,z,fi,q,time);
f641f6bd 480 return;
481 }
482
483 // if new volume or new trackID
b9671574 484 if ( (volumeID!=fTempInfo->GetParam()->GetVolumeID()) ||
485 (trackID!=fTempInfo->GetParam()->GetTrackID())){
486 if (volumeID!=fTempInfo->GetParam()->GetVolumeID()) AddVolume(volumeID);
f641f6bd 487 diff=kTRUE;
488 FlushHitStack(kTRUE);
489
b9671574 490 fTempInfo->SetParamIndex(fTempInfo->GetParamIndex()+1);
491 fTempInfo->SetParam(new((*fArray)[fTempInfo->GetParamIndex()]) AliTrackHitsParamV2);
492 fTempInfo->GetParam()->SetVolumeID(volumeID);
493 fTempInfo->GetParam()->SetTrackID(trackID);
e61fd20d 494 fTempInfo->NewParam(r,z,fi,q,time);
f641f6bd 495 return;
496 }
497
498 //calculate current fit precission to next point
b9671574 499 AliTrackHitsParamV2 &param = *(fTempInfo->GetParam());
f641f6bd 500 Double_t dd=0;
501 Double_t dl=0;
502 Double_t ratio=0;
503 Double_t dr,dz,dfi,ddz,ddfi;
504 Double_t drhit,ddl;
505 dr=dz=dfi=ddz=ddfi=0;
b9671574 506 drhit = r-fTempInfo->GetOldR();
f641f6bd 507 {
508 //Double_t dfi2 = param.fAn+2*param.fAd*(r-param.fR);
b9671574 509 Double_t dfi2 = param.GetAn();
510 dfi2*=dfi2*fTempInfo->GetOldR()*fTempInfo->GetOldR();
f641f6bd 511 //Double_t ddz2 = param.fTheta+2*param.fThetaD*(r-param.fR);
b9671574 512 Double_t ddz2 = param.GetTheta();
f641f6bd 513 ddz2*=ddz2;
514 ratio = TMath::Sqrt(1.+ dfi2+ ddz2);
515 }
c4a237ae 516 //
517 // dl = fStep * Short_t(TMath::Nint(drhit*ratio/fStep)); // MI change - range check
518 dl = drhit*ratio/fStep;
519 if (TMath::Abs(dl)>32765) dl =0;
520 dl = fStep * Short_t(TMath::Nint(dl));
521 //
f641f6bd 522 ddl = dl - drhit*ratio;
b9671574 523 fTempInfo->SetOldR(fTempInfo->GetOldR()+dl/ratio);
524
525 if (fTempInfo->GetStackIndex()>2){
526 dr = r-param.GetR();
527 dz = z-param.GetZ();
528 dfi = fi-param.GetFi();
529 ddz = dr*param.GetTheta()+dr*dr*param.GetThetaD()-dz;
530 ddfi= dr*param.GetAn()+dr*dr*param.GetAd()-dfi;
f641f6bd 531 dd = TMath::Sqrt(ddz*ddz+r*r*ddfi*ddfi+ddl*ddl);
532 //
533 }
534 //safety factor 1.25
535 if ( ( (dd*1.25>fPrecision) ) ||
b9671574 536 (fTempInfo->GetStackIndex()+4>fTempInfo->kStackSize) ||
f641f6bd 537 (TMath::Abs(dl/fStep)>fMaxDistance) )
538 diff=kTRUE;
539 else{ // if precision OK
b9671574 540 fTempInfo->SetStackIndex(fTempInfo->GetStackIndex()+1);
e61fd20d 541 fTempInfo->SetHit(r,z,fi,q,time);
f641f6bd 542 return;
543 }
544
545
546 //if parameter changed
547 if (FlushHitStack(kFALSE)){ //if full buffer flushed
b9671574 548 fTempInfo->SetParamIndex(fTempInfo->GetParamIndex()+1);
549 fTempInfo->SetParam(new((*fArray)[fTempInfo->GetParamIndex()]) AliTrackHitsParamV2);
550 fTempInfo->GetParam()->SetVolumeID(volumeID);
551 fTempInfo->GetParam()->SetTrackID(trackID);
e61fd20d 552 fTempInfo->NewParam(r,z,fi,q,time);
f641f6bd 553 }
554 else{
b9671574 555 fTempInfo->SetStackIndex(fTempInfo->GetStackIndex()+1);
e61fd20d 556 fTempInfo->SetHit(r,z,fi,q,time);
f641f6bd 557 }
558}
559
560Bool_t AliTPCTrackHitsV2::FlushHitStack(Bool_t force)
561{
562 //
53110b9d 563 // write fHitsPosAndQ information from the stack to te arrays
564 //
f641f6bd 565 if (!fTempInfo) return kFALSE;
566
b9671574 567 AliTrackHitsParamV2 & param = *(fTempInfo->GetParam());
f641f6bd 568 //recalculate track parameter not fixing first point
569 fTempInfo->UpdateParam(fStep/4.);
570 //fTempInfo->Fit(fTempInfo->fParam); //- fixing the first and the last point
571
b9671574 572 Double_t oldr = param.GetR();
f641f6bd 573 UInt_t i;
574 Double_t dd;
b9671574 575 param.SetNHits(fTempInfo->GetStackIndex()+1);
576 // if (param.fHitDistance) delete []param.fHitDistance;
577 // if (param.fCharge) delete []param.fCharge;
578 // if (param.fTime) delete []param.fTime;
579 param.SetHitDistance(param.GetNHits());
580 param.SetCharge(param.GetNHits());
581 param.SetTime(param.GetNHits());
f641f6bd 582
583
b9671574 584 for (i=0; i <= fTempInfo->GetStackIndex(); i++){
f641f6bd 585 Double_t * position = fTempInfo->GetPosition(i);
586 Double_t dr = position[0]-oldr;
587 Double_t ratio;
588 {
589 //Double_t dfi2 = param.fAn+2*param.fAd*(position[0]-param.fR);
b9671574 590 Double_t dfi2 = param.GetAn();
f641f6bd 591 dfi2*=dfi2*oldr*oldr;
592 //Double_t ddz2 = param.fTheta+2*param.fThetaD*(position[0]-param.fR);
b9671574 593 Double_t ddz2 = param.GetTheta();
f641f6bd 594 ddz2*=ddz2;
595 ratio = TMath::Sqrt(1.+ dfi2+ ddz2);
596 }
597
595c3700 598 // Double_t dl = fStep*(Short_t)TMath::Nint(dr*ratio/fStep); //MI change
599 Double_t dl = dr*ratio/fStep;
600 if (TMath::Abs(dl)>32765) dl =0;
601 dl = fStep * Short_t(TMath::Nint(dl));
602
f641f6bd 603 dr = dl/ratio;
604 oldr+=dr;
605 //calculate precission
77f88633 606 AliTrackHitsParamV2 &paraml = *(fTempInfo->GetParam());
f641f6bd 607 //real deltas
77f88633 608 Double_t dr1= position[0]-paraml.GetR();
609 Double_t dz = position[1]-paraml.GetZ();
610 Double_t dfi = position[2]-paraml.GetFi();
f641f6bd 611 //extrapolated deltas
77f88633 612 Double_t dr2 = oldr-paraml.GetR();
f641f6bd 613 Double_t ddr = dr2-dr1;
77f88633 614 Double_t ddz = dr2*paraml.GetTheta()+dr2*dr2*paraml.GetThetaD()-dz;
615 Double_t ddfi= dr2*paraml.GetAn()+dr2*dr2*paraml.GetAd()-dfi;
f641f6bd 616 dd = TMath::Sqrt(ddz*ddz+oldr*oldr*ddfi*ddfi+ddr*ddr);
617
618
619 if ( (dd>fPrecision) ){
620 //if ( (dd<0) ){
621 if (i==0){
77f88633 622 paraml.SetAn(0);
623 paraml.SetAd(0);
624 paraml.SetTheta(0);
625 paraml.SetThetaD(0);
626 Double_t ddz1 = dr2*paraml.GetTheta()+dr2*dr2*paraml.GetThetaD()-dz;
627 Double_t ddfi1= dr2*paraml.GetAn()+dr2*dr2*paraml.GetAd()-dfi;
f641f6bd 628 dl = 0;
77f88633 629 dd = TMath::Sqrt(ddz1*ddz1+oldr*oldr*ddfi1*ddfi1+ddr*ddr);
f641f6bd 630 }
631 else
632 break;
633 }
634
77f88633 635 paraml.HitDistance(i)= Short_t(TMath::Nint(dl/fStep));
636 paraml.Charge(i)= Short_t(fTempInfo->GetQStack(i));
637 paraml.Time(i)= Short_t(fTempInfo->GetTimeStack(i)/AliTPCTrackHitsV2::fgkTimePrecision);
f641f6bd 638 }
639
b9671574 640 if (i<=fTempInfo->GetStackIndex()){ //if previous iteration not succesfull
641 // Short_t * charge = new Short_t[i];
642 // Short_t * time = new Short_t[i];
643 // Short_t * hitDistance= new Short_t[i];
644 // memcpy(charge, param.fCharge,sizeof(Short_t)*i);
645 // memcpy(time, param.fTime,sizeof(Short_t)*i);
646 // memcpy(hitDistance, param.fHitDistance,sizeof(Short_t)*i);
647 // delete [] param.fCharge;
648 // delete [] param.fTime;
649 // delete [] param.fHitDistance;
650 param.SetNHits(i);
651 param.ResizeCharge(i);
652 param.ResizeTime(i);
653 param.ResizeHitDistance(i);
f641f6bd 654 //
b9671574 655 Int_t volumeID = fTempInfo->GetParam()->GetVolumeID();
656 Int_t trackID =fTempInfo->GetParam()->GetTrackID();
657 fTempInfo->SetParamIndex(fTempInfo->GetParamIndex()+1);
658 fTempInfo->SetParam(new((*fArray)[fTempInfo->GetParamIndex()]) AliTrackHitsParamV2);
f641f6bd 659 Double_t * p = fTempInfo->GetPosition(i);
b9671574 660 UInt_t index2 = fTempInfo->GetStackIndex();
661 fTempInfo->NewParam(p[0],p[1],p[2],fTempInfo->GetQStack(i),fTempInfo->GetTimeStack(i));
662 fTempInfo->GetParam()->SetVolumeID(volumeID);
663 fTempInfo->GetParam()->SetTrackID(trackID);
f641f6bd 664 if (i+1<=index2) FlushHitStack2(i+1,index2);
665
666 if (force) return FlushHitStack(kTRUE);
667 return kFALSE;
668 }
669 return kTRUE;
670}
671
672
673void AliTPCTrackHitsV2::FlushHitStack2(Int_t index1, Int_t index2)
674{
675 //
676 // second iteration flush stack
53110b9d 677 // call only for hits where first iteration were not succesfully interpolated
678 //
f641f6bd 679 Double_t * positionstack = new Double_t[3*(index2-index1+1)];
680 UInt_t * qstack = new UInt_t[index2-index1+1];
e61fd20d 681 Float_t * timestack = new Float_t[index2-index1+1];
b9671574 682 memcpy(positionstack, fTempInfo->GetPosition(index1),
f641f6bd 683 (3*(index2-index1+1))*sizeof(Double_t));
b9671574 684 memcpy(qstack, fTempInfo->GetQStackP(index1),(index2-index1+1)*sizeof(UInt_t));
685 memcpy(timestack, fTempInfo->GetTimeStackP(index1),(index2-index1+1)*sizeof(Float_t));
f641f6bd 686 Double_t *p = positionstack;
687 for (Int_t j=0; j<=index2-index1;j++){
b9671574 688 fTempInfo->SetStackIndex(fTempInfo->GetStackIndex()+1);
e61fd20d 689 fTempInfo->SetHit(p[3*j+0],p[3*j+1],p[3*j+2],qstack[j],timestack[j]);
f641f6bd 690 }
691 delete []positionstack;
692 delete []qstack;
e61fd20d 693 delete []timestack;
f641f6bd 694}
695
696
697void AliTPCTrackHitsV2::AddVolume(Int_t volume)
698{
699 //
700 //add volumes to tthe list of volumes
53110b9d 701 //
f641f6bd 702 Int_t * volumes = new Int_t[fNVolumes+1];
f55b9162 703 if (fVolumes) memcpy(volumes,fVolumes,(fNVolumes)*sizeof(Int_t));
f641f6bd 704 volumes[fNVolumes]=volume;
705 fNVolumes++;
706 if (fVolumes) delete []fVolumes;
707 fVolumes = volumes;
708}
709
710
f641f6bd 711Bool_t AliTPCTrackHitsV2::First()
712{
713 //
714 //set Current hit for the first hit
715 //
6e493dba 716
717 if (fArray->GetSize()<=0) {
b9671574 718 fCurrentHit->SetStatus(kFALSE);
6e493dba 719 return kFALSE;
720 }
721
f641f6bd 722 AliTrackHitsParamV2 *param = (AliTrackHitsParamV2 *)fArray->At(0);
fce0283f 723 if (!fHit) fHit = new AliTPChit;
f641f6bd 724 if (!(param) ) {
b9671574 725 fCurrentHit->SetStatus(kFALSE);
f641f6bd 726 return kFALSE;
727 }
728 //
b9671574 729 fCurrentHit->SetParamIndex(0);
730 fCurrentHit->SetStackIndex(0);
f641f6bd 731 //
fce0283f 732 //
b9671574 733 ((AliTPChit*)fHit)->fSector = param->GetVolumeID();
734 ((AliTPChit*)fHit)->SetTrack(param->GetTrackID());
735 ((AliTPChit*)fHit)->SetX(param->GetR()*TMath::Cos(param->GetFi()));
736 ((AliTPChit*)fHit)->SetY(param->GetR()*TMath::Sin(param->GetFi()));
737 ((AliTPChit*)fHit)->SetZ(param->GetZ());
738 ((AliTPChit*)fHit)->fQ = param->Charge(0);
739 ((AliTPChit*)fHit)->fTime = (Float_t)(param->Time(0)*AliTPCTrackHitsV2::fgkTimePrecision);
fce0283f 740 /*
741 fCurrentHit->fHit.fSector = param->fVolumeID;
742 fCurrentHit->fHit.SetTrack(param->fTrackID);
743 fCurrentHit->fHit.SetX(param->fR*TMath::Cos(param->fFi));
744 fCurrentHit->fHit.SetY(param->fR*TMath::Sin(param->fFi));
745 fCurrentHit->fHit.SetZ(param->fZ);
746 fCurrentHit->fHit.fQ = param->fCharge[0];
e61fd20d 747 fCurrentHit->fHit.fTime = (Float_t)(param->fTime[0]*AliTPCTrackHitsV2::fgkTimePrecision);
fce0283f 748 */
b9671574 749 fCurrentHit->SetR(param->GetR());
f641f6bd 750
b9671574 751 fCurrentHit->SetStatus(kTRUE);
752 return fCurrentHit->GetStatus();
f641f6bd 753}
754
755Bool_t AliTPCTrackHitsV2::Next()
756{
757 //
53110b9d 758 // Hit iterator
759 //
b9671574 760 if (!(fCurrentHit->GetStatus()))
f641f6bd 761 return kFALSE;
762
b9671574 763 fCurrentHit->SetStackIndex(fCurrentHit->GetStackIndex()+1);
f641f6bd 764
b9671574 765 AliTrackHitsParamV2 *param = (AliTrackHitsParamV2 *)fArray->At(fCurrentHit->GetParamIndex());
766 if (fCurrentHit->GetStackIndex()>=param->GetNHits()){
767 fCurrentHit->SetParamIndex(fCurrentHit->GetParamIndex()+1);
768 if (fCurrentHit->GetParamIndex()>=fArray->GetEntriesFast()){
769 fCurrentHit->SetStatus(kFALSE);
f641f6bd 770 return kFALSE;
771 }
b9671574 772 param = (AliTrackHitsParamV2 *)fArray->At(fCurrentHit->GetParamIndex());
773 fCurrentHit->SetStackIndex(0);
774 fCurrentHit->SetR(param->GetR());
f641f6bd 775 }
776
777
778
779 Double_t ratio;
780 {
781 // Double_t dfi2 = param->fAn+2*param->fAd*(fCurrentHit->fR-param->fR);
b9671574 782 Double_t dfi2 = param->GetAn();
783 dfi2*=dfi2*fCurrentHit->GetR()*fCurrentHit->GetR();
f641f6bd 784 // Double_t ddz2 = param->fTheta+2*param->fThetaD*(fCurrentHit->fR-param->fR);
b9671574 785 Double_t ddz2 = param->GetTheta();
f641f6bd 786 ddz2*=ddz2;
787 ratio = TMath::Sqrt(1.+ dfi2+ ddz2);
788 }
789
b9671574 790 fCurrentHit->SetR(fCurrentHit->GetR()+fStep*param->HitDistance(fCurrentHit->GetStackIndex())/ratio);
f641f6bd 791
b9671574 792 Double_t dR = fCurrentHit->GetR() - param->GetR();
793 Double_t fi = param->GetFi() + (param->GetAn()*dR+param->GetAd()*dR*dR);
794 Double_t z = param->GetZ() + (param->GetTheta()*dR+param->GetThetaD()*dR*dR);
fce0283f 795 /*
f641f6bd 796 fCurrentHit->fHit.fQ = param->fCharge[fCurrentHit->fStackIndex];
e61fd20d 797 fCurrentHit->fHit.fTime = (Float_t)(param->fTime[fCurrentHit->fStackIndex]*AliTPCTrackHitsV2::fgkTimePrecision);
f641f6bd 798 fCurrentHit->fHit.SetX(fCurrentHit->fR*TMath::Cos(fi));
799 fCurrentHit->fHit.SetY(fCurrentHit->fR*TMath::Sin(fi));
800 fCurrentHit->fHit.SetZ(z);
801 fCurrentHit->fHit.fSector = param->fVolumeID;
802 fCurrentHit->fHit.SetTrack(param->fTrackID);
fce0283f 803 */
b9671574 804 ((AliTPChit*)fHit)->fQ = param->Charge(fCurrentHit->GetStackIndex());
805 ((AliTPChit*)fHit)->fTime = (Float_t)(param->Time(fCurrentHit->GetStackIndex())*AliTPCTrackHitsV2::fgkTimePrecision);
806 ((AliTPChit*)fHit)->SetX(fCurrentHit->GetR()*TMath::Cos(fi));
807 ((AliTPChit*)fHit)->SetY(fCurrentHit->GetR()*TMath::Sin(fi));
fce0283f 808 ((AliTPChit*)fHit)->SetZ(z);
b9671574 809 ((AliTPChit*)fHit)->fSector = param->GetVolumeID();
810 ((AliTPChit*)fHit)->SetTrack(param->GetTrackID());
fce0283f 811
f641f6bd 812 return kTRUE;
813}
814
53110b9d 815AliHit * AliTPCTrackHitsV2::GetHit() const
f641f6bd 816{
817 //
53110b9d 818 // Return one hit
819 //
b9671574 820 return (fCurrentHit->GetStatus())? fHit:0;
f641f6bd 821 //return &fCurrentHit->fHit;
822
823}
824
825AliTrackHitsParamV2 * AliTPCTrackHitsV2::GetParam()
826{
53110b9d 827 //
828 // Return current parameters
829 //
b9671574 830 return (fCurrentHit->GetStatus())?
831 (AliTrackHitsParamV2 *)fArray->At(fCurrentHit->GetParamIndex()):0;
f641f6bd 832}
53110b9d 833