Line data Source code
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 : /* $Id$ */
17 :
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 : ///////////////////////////////////////////////////////////////////////////////
54 : //
55 :
56 : #include <TClonesArray.h>
57 : #include <TMath.h>
58 :
59 : #include "AliTPCTrackHitsV2.h"
60 : #include "AliTPC.h"
61 :
62 :
63 :
64 12 : ClassImp(AliTPCTrackHitsV2)
65 12 : ClassImp(AliTrackHitsParamV2)
66 :
67 : //
68 : Int_t AliTrackHitsParamV2::fgCounter1 =0;
69 : Int_t AliTrackHitsParamV2::fgCounter2 =0;
70 : //
71 : Int_t AliTPCTrackHitsV2::fgCounter1 =0;
72 : Int_t AliTPCTrackHitsV2::fgCounter2 =0;
73 : //
74 : const Double_t AliTPCTrackHitsV2::fgkPrecision=1e-6; //precision
75 : const Double_t AliTPCTrackHitsV2::fgkPrecision2=1e-20; //precision
76 : const Double_t AliTPCTrackHitsV2::fgkTimePrecision=20.e-9; //hit time precision
77 :
78 :
79 :
80 :
81 : class AliTPCTempHitInfoV2 {
82 : public:
83 : AliTPCTempHitInfoV2();
84 : void SetHit(Double_t r, Double_t z, Double_t fi, Int_t q, Float_t time);
85 9170704 : UInt_t GetStackIndex() const {return fStackIndex;}
86 4307844 : void SetStackIndex(UInt_t i) {fStackIndex=i;}
87 281640 : UInt_t GetParamIndex() const {return fParamIndex;}
88 140942 : void SetParamIndex(UInt_t i) {fParamIndex=i;}
89 1588828 : Float_t GetTimeStack(Int_t i) const {return fTimeStack[i];}
90 107110 : const Float_t* GetTimeStackP(Int_t i) const {return &fTimeStack[i];}
91 1588828 : UInt_t GetQStack(Int_t i) const {return fQStack[i];}
92 107110 : const UInt_t* GetQStackP(Int_t i) const {return &fQStack[i];}
93 6253364 : Double_t * GetPosition(Int_t index){return &fPositionStack[index*3];}
94 5915944 : Double_t GetOldR() const {return fOldR;}
95 1478986 : void SetOldR(Double_t r) {fOldR=r;}
96 :
97 :
98 6668536 : AliTrackHitsParamV2 * GetParam() const {return fParam;}
99 140942 : 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 :
104 : protected:
105 : AliTPCTempHitInfoV2(const AliTPCTempHitInfoV2 &hit);
106 : AliTPCTempHitInfoV2& operator = (const AliTPCTempHitInfoV2 &hit);
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);
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
122 : Double_t fOldR; //previos r
123 : Double_t fPositionStack[3*kStackSize]; //position stack
124 : UInt_t fQStack[kStackSize]; //Q stack
125 : Float_t fTimeStack[kStackSize]; //time stack
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 :
134 : AliTPCTempHitInfoV2::AliTPCTempHitInfoV2()
135 122 : :fSumDr(0.),
136 61 : fSumDr2(0.),
137 61 : fSumDr3(0.),
138 61 : fSumDr4(0.),
139 61 : fSumDFi(0.),
140 61 : fSumDFiDr(0.),
141 61 : fSumDFiDr2(0.),
142 61 : fSumDZ(0.),
143 61 : fSumDZDr(0.),
144 61 : fSumDZDr2(0.),
145 61 : fOldR(0.),
146 61 : fStackIndex(0),
147 61 : fParamIndex(0),
148 61 : fParam(0)
149 122 : {
150 : //
151 : // Standard constructor
152 : // set to default value
153 : //
154 61 : fSumDr=fSumDr2=fSumDr3=fSumDr4=
155 61 : fSumDFi=fSumDFiDr=fSumDFiDr2=
156 61 : fSumDZ=fSumDZDr=fSumDZDr2=0;
157 61 : fStackIndex = 0;
158 : // fInfoIndex = 0;
159 61 : fParamIndex = 0;
160 3660122 : for(Int_t i=0;i<3*kStackSize;i++) fPositionStack[i]=0.;
161 1220122 : for(Int_t i=0;i<kStackSize;i++){
162 610000 : fQStack[i]=0;
163 610000 : fTimeStack[i]=0.;
164 : }
165 122 : }
166 :
167 :
168 : void AliTPCTempHitInfoV2::NewParam(Double_t r, Double_t z, Double_t fi, Int_t q, Float_t time)
169 : {
170 : //
171 : //reset stack and sum parameters
172 : //store line initial point
173 : //
174 70471 : fSumDr=fSumDr2=fSumDr3=fSumDr4=
175 70471 : fSumDFi=fSumDFiDr=fSumDFiDr2=
176 140942 : fSumDZ=fSumDZDr=fSumDZDr2=0;
177 70471 : fStackIndex=0;
178 70471 : fParam->SetR(r);
179 70471 : fOldR = r;
180 70471 : fParam->SetZ(z);
181 70471 : fParam->SetFi(fi);
182 70471 : fParam->SetAn(0.);
183 70471 : fParam->SetAd(0.);
184 70471 : fParam->SetTheta(0.);
185 70471 : fParam->SetThetaD(0.);
186 70471 : SetHit(r,z,fi,q,time);
187 70471 : }
188 :
189 : void AliTPCTempHitInfoV2::SetHit(Double_t r, Double_t z, Double_t fi, Int_t q, Float_t time)
190 : {
191 : //
192 : //add hit to the stack
193 : //recalculate new estimete of line parameters
194 4448786 : Double_t *f = GetPosition(fStackIndex);
195 2224393 : f[0] = r;
196 2224393 : f[1] = z;
197 2224393 : f[2] = fi;
198 2224393 : fQStack[fStackIndex]=q;
199 2224393 : fTimeStack[fStackIndex]=time;
200 2294864 : if (fStackIndex==0) return;
201 2153922 : Double_t dr = (r-fParam->GetR());
202 2154341 : if (TMath::Abs(dr)<AliTPCTrackHitsV2::GetKPrecision()) dr =AliTPCTrackHitsV2::GetKPrecision();
203 2153922 : Double_t dfi = fi-fParam->GetFi();
204 2153922 : Double_t dz = z -fParam->GetZ();
205 2153922 : Double_t dr2 =dr*dr;
206 2153922 : Double_t dr3 =dr2*dr;
207 2153922 : Double_t dr4 =dr3*dr;
208 2153922 : fSumDr +=dr;
209 2153922 : fSumDr2+=dr2;
210 2153922 : fSumDr3+=dr3;
211 2153922 : fSumDr4+=dr4;
212 2153922 : fSumDFi +=dfi;
213 2153922 : fSumDFiDr+=dfi*dr;
214 2153922 : fSumDFiDr2+=dfi*dr2;
215 2153922 : fSumDZ +=dz;
216 2153922 : fSumDZDr+=dz*dr;
217 2153922 : fSumDZDr2+=dz*dr2;
218 :
219 : //update fit parameters
220 : //
221 2153922 : Double_t det = fSumDr2*fSumDr4-fSumDr3*fSumDr3;
222 2303429 : if (TMath::Abs(det)<AliTPCTrackHitsV2::GetKPrecision2()) return;
223 4008830 : if ( ( fStackIndex>1 ) ){
224 4008798 : fParam->SetAn((fSumDr4*fSumDFiDr-fSumDr3*fSumDFiDr2)/det);
225 2004383 : fParam->SetAd((fSumDr2*fSumDFiDr2-fSumDr3*fSumDFiDr)/det);
226 2004383 : }
227 : else
228 32 : fParam->SetAn(fSumDFiDr/fSumDr2);
229 4008830 : if ( ( fStackIndex>1 ) ){
230 4008798 : fParam->SetTheta((fSumDr4*fSumDZDr-fSumDr3*fSumDZDr2)/det);
231 2004383 : fParam->SetThetaD((fSumDr2*fSumDZDr2-fSumDr3*fSumDZDr)/det);
232 2004383 : }
233 : else
234 32 : fParam->SetTheta(fSumDZDr/fSumDr2);
235 4228808 : }
236 :
237 :
238 : void AliTPCTempHitInfoV2::UpdateParam(Double_t maxdelta)
239 : {
240 : //
241 : // recalc parameters not fixing origin point
242 : //
243 140942 : if (fStackIndex>5){
244 60900 : Double_t a,b,c;
245 60900 : a=b=c=0;
246 121800 : Fit2(fSumDFi, fSumDFiDr, fSumDFiDr2, fSumDr,fSumDr2,fSumDr3,fSumDr4,
247 60900 : fStackIndex, a,b,c);
248 60900 : if (TMath::Abs(a)<maxdelta){
249 60884 : fParam->SetFi(fParam->GetFi()+a/fParam->GetR());
250 60884 : fParam->SetAn(b);
251 60884 : fParam->SetAd(c);
252 60884 : }
253 121800 : Fit2(fSumDZ, fSumDZDr, fSumDZDr2, fSumDr,fSumDr2,fSumDr3,fSumDr4,
254 60900 : fStackIndex, a,b,c) ;
255 60900 : if (TMath::Abs(a)<maxdelta){
256 54084 : fParam->SetZ(fParam->GetZ()+a);
257 54084 : fParam->SetTheta(b);
258 54084 : fParam->SetThetaD(c);
259 54084 : }
260 60900 : }
261 :
262 70471 : }
263 :
264 : void 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 : {
269 : //
270 : // fit of second order
271 : //
272 : Double_t det =
273 365400 : n* (fSumX2*fSumX4-fSumX3*fSumX3) -
274 243600 : fSumX* (fSumX*fSumX4-fSumX3*fSumX2)+
275 121800 : fSumX2* (fSumX*fSumX3-fSumX2*fSumX2);
276 :
277 121800 : if (TMath::Abs(det)> AliTPCTrackHitsV2::GetKPrecision()) {
278 25286 : a =
279 50572 : (fSumY * (fSumX2*fSumX4-fSumX3*fSumX3)-
280 50572 : fSumX *(fSumYX*fSumX4-fSumYX2*fSumX3)+
281 50572 : fSumX2*(fSumYX*fSumX3-fSumYX2*fSumX2))/det;
282 25286 : b=
283 50572 : (n*(fSumYX*fSumX4-fSumX3*fSumYX2)-
284 50572 : fSumY*(fSumX*fSumX4-fSumX3*fSumX2)+
285 50572 : fSumX2*(fSumX*fSumYX2-fSumYX*fSumX2))/det;
286 25286 : c=
287 50572 : (n*(fSumX2*fSumYX2-fSumYX*fSumX3)-
288 50572 : fSumX*(fSumX*fSumYX2-fSumYX*fSumX2)+
289 50572 : fSumY*(fSumX*fSumX3-fSumX2*fSumX2))/det;
290 25286 : }
291 121800 : }
292 :
293 : void AliTPCTempHitInfoV2::Fit(AliTrackHitsParamV2 * param)
294 : {
295 : //
296 : // fit fixing first and the last point
297 : // result stored in new param
298 : //
299 0 : Double_t dx2 = (GetPosition(fStackIndex))[0]-fParam->GetR();
300 0 : Double_t det = fSumDr4+dx2*fSumDr2-2*dx2*fSumDr3;
301 0 : if ( (TMath::Abs(det)> AliTPCTrackHitsV2::GetKPrecision()) &&
302 0 : ((TMath::Abs(dx2)> AliTPCTrackHitsV2::GetKPrecision()))){
303 0 : Double_t dfi2 = (GetPosition(fStackIndex))[1]-fParam->GetFi();
304 0 : param->SetAd((fSumDFiDr2+dfi2*fSumDr-dx2*fSumDFiDr-dfi2*fSumDr3/dx2)/det);
305 0 : param->SetAn((dfi2-param->GetAd()*dx2*dx2)/dx2);
306 :
307 0 : Double_t dz2 = (GetPosition(fStackIndex))[1]-fParam->GetZ();
308 0 : param->SetTheta((fSumDZDr2+dz2*fSumDr-dx2*fSumDZDr-dz2*fSumDr3/dx2)/det);
309 0 : param->SetTheta((dz2-param->GetAd()*dx2*dx2)/dx2);
310 0 : }
311 :
312 0 : }
313 :
314 758310 : AliTrackHitsParamV2::AliTrackHitsParamV2():TObject(),
315 758310 : fTrackID(0),
316 758310 : fVolumeID(0),
317 758310 : fR(0.),
318 758310 : fZ(0.),
319 758310 : fFi(0.),
320 758310 : fAn(0.),
321 758310 : fAd(0.),
322 758310 : fTheta(0.),
323 758310 : fThetaD(0.),
324 758310 : fNHits(0),
325 758310 : fHitDistance(0),
326 758310 : fCharge(0),
327 758310 : fTime(0)
328 3791550 : {
329 : //
330 : // default constructor
331 : //
332 758310 : fgCounter1++;
333 758310 : fgCounter2++;
334 1516620 : }
335 :
336 : AliTrackHitsParamV2::~AliTrackHitsParamV2()
337 3076508 : {
338 : //
339 : // Standard destructor
340 : //
341 758310 : fgCounter1--;
342 758310 : if (fHitDistance) {
343 1516618 : delete[]fHitDistance;
344 758309 : fHitDistance=0;
345 758309 : }
346 758310 : if (fCharge){
347 1516618 : delete[]fCharge;
348 758309 : fCharge =0;
349 758309 : }
350 758310 : if (fTime){
351 1516618 : delete[]fTime;
352 758309 : fTime =0;
353 758309 : }
354 1538254 : }
355 :
356 : Float_t AliTrackHitsParamV2::Eta() const
357 : {
358 0 : Float_t ctg = fZ / fR;
359 0 : Float_t eta = -TMath::Log(TMath::Hypot(1,ctg)-TMath::Abs(ctg));
360 0 : if(ctg < 0) eta = -eta;
361 0 : return eta;
362 : }
363 :
364 :
365 3 : AliTPCTrackHitsV2::AliTPCTrackHitsV2():TObject(),
366 3 : fArray(0),
367 3 : fSize(0),
368 3 : fPrecision(0.),
369 3 : fStep(0.),
370 3 : fMaxDistance(0),
371 3 : fNVolumes(0),
372 3 : fVolumes(0),
373 3 : fTempInfo(0),
374 3 : fCurrentHit(0),
375 3 : fHit(0)
376 15 : {
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 3 : fPrecision=kHitPrecision; //precision in cm
385 3 : fStep = kStep; //step size
386 3 : fMaxDistance = kMaxDistance; //maximum distance
387 :
388 : //fTrackHitsInfo = new AliObjectArray("AliTrackHitsInfo");
389 : //fTrackHitsParam = new AliObjectArray("AliTrackHitsParamV2");
390 : //fHitsPosAndQ = new TArrayOfArrayVStack("AliHitInfo");
391 9 : fArray = new TClonesArray("AliTrackHitsParamV2");
392 6 : fCurrentHit = new AliTPCCurrentHitV2;
393 3 : fgCounter1++;
394 3 : fgCounter2++;
395 :
396 6 : }
397 :
398 : AliTPCTrackHitsV2::~AliTPCTrackHitsV2()
399 18 : {
400 : //
401 : //default destructor
402 : //
403 : // if (fTrackHitsInfo) delete fTrackHitsInfo;
404 3 : if (fArray) {
405 6 : delete fArray;
406 3 : fArray =0;
407 3 : }
408 : //if (fHitsPosAndQ) delete fHitsPosAndQ;
409 9 : if (fCurrentHit) delete fCurrentHit;
410 5 : if (fTempInfo) delete fTempInfo;
411 3 : if (fVolumes) {
412 6 : delete [] fVolumes;
413 3 : fVolumes =0;
414 3 : fNVolumes=0;
415 3 : }
416 3 : if (fHit){
417 4 : delete fHit;
418 2 : fHit=0;
419 2 : }
420 3 : fgCounter1--;
421 9 : }
422 :
423 : void AliTPCTrackHitsV2::Clear(Option_t * /*option*/)
424 : {
425 : //
426 : // clear object
427 : //
428 16360 : fSize = 0;
429 8180 : if (fArray){
430 1489712 : for (Int_t i=0;i<fArray->GetEntriesFast();i++){
431 736676 : AliTrackHitsParamV2 * par = (AliTrackHitsParamV2 *)fArray->UncheckedAt(i);
432 736676 : par->~AliTrackHitsParamV2(); // delete object
433 : }
434 8180 : fArray->Clear();
435 8180 : }
436 8180 : if (fTempInfo){
437 120 : delete fTempInfo;
438 60 : delete fHit;
439 60 : fHit =0;
440 60 : fTempInfo =0;
441 60 : }
442 8180 : if (fVolumes){
443 8902 : delete [] fVolumes;
444 4451 : fVolumes=0;
445 4451 : fNVolumes=0;
446 4451 : }
447 8180 : }
448 :
449 :
450 : void AliTPCTrackHitsV2::AddHitKartez(Int_t volumeID, Int_t trackID, Double_t x,
451 : Double_t y, Double_t z,Int_t q, Float_t time)
452 : {
453 : //
454 : // add hit to the container - it add hit at the end - input in global coordinata
455 : //
456 1480188 : Double_t r = TMath::Sqrt(x*x+y*y);
457 740094 : Double_t fi = TMath::ACos(x/r);
458 964609 : if (y<0) fi*=-1.;
459 740094 : AddHit(volumeID,trackID,r,z,fi,q,time);
460 740094 : }
461 :
462 :
463 : void AliTPCTrackHitsV2::AddHit(Int_t volumeID, Int_t trackID,
464 : Double_t r, Double_t z, Double_t fi, Int_t q, Float_t time)
465 : {
466 : //
467 : // Adding one hit
468 : //
469 1480188 : fSize++;
470 : Bool_t diff=kFALSE;
471 740094 : if (!fTempInfo) { //initialisation of track - initialisation of parameters
472 61 : fTempInfo = new AliTPCTempHitInfoV2;
473 122 : fTempInfo->SetParam(new((*fArray)[0]) AliTrackHitsParamV2);
474 61 : fTempInfo->GetParam()->SetVolumeID(volumeID);
475 61 : fTempInfo->GetParam()->SetTrackID(trackID);
476 61 : AddVolume(volumeID);
477 : //
478 61 : fTempInfo->SetParamIndex(0);
479 61 : fTempInfo->NewParam(r,z,fi,q,time);
480 61 : return;
481 : }
482 :
483 : // if new volume or new trackID
484 1479768 : if ( (volumeID!=fTempInfo->GetParam()->GetVolumeID()) ||
485 739735 : (trackID!=fTempInfo->GetParam()->GetTrackID())){
486 838 : if (volumeID!=fTempInfo->GetParam()->GetVolumeID()) AddVolume(volumeID);
487 : diff=kTRUE;
488 540 : FlushHitStack(kTRUE);
489 :
490 540 : fTempInfo->SetParamIndex(fTempInfo->GetParamIndex()+1);
491 1080 : fTempInfo->SetParam(new((*fArray)[fTempInfo->GetParamIndex()]) AliTrackHitsParamV2);
492 540 : fTempInfo->GetParam()->SetVolumeID(volumeID);
493 540 : fTempInfo->GetParam()->SetTrackID(trackID);
494 540 : fTempInfo->NewParam(r,z,fi,q,time);
495 540 : return;
496 : }
497 :
498 : //calculate current fit precission to next point
499 739493 : AliTrackHitsParamV2 ¶m = *(fTempInfo->GetParam());
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;
506 739493 : drhit = r-fTempInfo->GetOldR();
507 : {
508 : //Double_t dfi2 = param.fAn+2*param.fAd*(r-param.fR);
509 739493 : Double_t dfi2 = param.GetAn();
510 739493 : dfi2*=dfi2*fTempInfo->GetOldR()*fTempInfo->GetOldR();
511 : //Double_t ddz2 = param.fTheta+2*param.fThetaD*(r-param.fR);
512 739493 : Double_t ddz2 = param.GetTheta();
513 739493 : ddz2*=ddz2;
514 739493 : ratio = TMath::Sqrt(1.+ dfi2+ ddz2);
515 : }
516 : //
517 : // dl = fStep * Short_t(TMath::Nint(drhit*ratio/fStep)); // MI change - range check
518 739493 : dl = drhit*ratio/fStep;
519 739555 : if (TMath::Abs(dl)>32765) dl =0;
520 739493 : dl = fStep * Short_t(TMath::Nint(dl));
521 : //
522 739493 : ddl = dl - drhit*ratio;
523 739493 : fTempInfo->SetOldR(fTempInfo->GetOldR()+dl/ratio);
524 :
525 739493 : if (fTempInfo->GetStackIndex()>2){
526 688609 : dr = r-param.GetR();
527 688609 : dz = z-param.GetZ();
528 688609 : dfi = fi-param.GetFi();
529 688609 : ddz = dr*param.GetTheta()+dr*dr*param.GetThetaD()-dz;
530 688609 : ddfi= dr*param.GetAn()+dr*dr*param.GetAd()-dfi;
531 688609 : dd = TMath::Sqrt(ddz*ddz+r*r*ddfi*ddfi+ddl*ddl);
532 : //
533 688609 : }
534 : //safety factor 1.25
535 1425603 : if ( ( (dd*1.25>fPrecision) ) ||
536 686110 : (fTempInfo->GetStackIndex()+4>fTempInfo->kStackSize) ||
537 686110 : (TMath::Abs(dl/fStep)>fMaxDistance) )
538 : diff=kTRUE;
539 : else{ // if precision OK
540 670300 : fTempInfo->SetStackIndex(fTempInfo->GetStackIndex()+1);
541 670300 : fTempInfo->SetHit(r,z,fi,q,time);
542 670300 : return;
543 : }
544 :
545 :
546 : //if parameter changed
547 138386 : if (FlushHitStack(kFALSE)){ //if full buffer flushed
548 84743 : fTempInfo->SetParamIndex(fTempInfo->GetParamIndex()+1);
549 31100 : fTempInfo->SetParam(new((*fArray)[fTempInfo->GetParamIndex()]) AliTrackHitsParamV2);
550 15550 : fTempInfo->GetParam()->SetVolumeID(volumeID);
551 15550 : fTempInfo->GetParam()->SetTrackID(trackID);
552 15550 : fTempInfo->NewParam(r,z,fi,q,time);
553 15550 : }
554 : else{
555 53643 : fTempInfo->SetStackIndex(fTempInfo->GetStackIndex()+1);
556 53643 : fTempInfo->SetHit(r,z,fi,q,time);
557 : }
558 809287 : }
559 :
560 : Bool_t AliTPCTrackHitsV2::FlushHitStack(Bool_t force)
561 : {
562 : //
563 : // write fHitsPosAndQ information from the stack to te arrays
564 : //
565 70573 : if (!fTempInfo) return kFALSE;
566 :
567 70471 : AliTrackHitsParamV2 & param = *(fTempInfo->GetParam());
568 : //recalculate track parameter not fixing first point
569 70471 : fTempInfo->UpdateParam(fStep/4.);
570 : //fTempInfo->Fit(fTempInfo->fParam); //- fixing the first and the last point
571 :
572 70471 : Double_t oldr = param.GetR();
573 : UInt_t i;
574 : Double_t dd;
575 70471 : 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 70471 : param.SetHitDistance(param.GetNHits());
580 70471 : param.SetCharge(param.GetNHits());
581 70471 : param.SetTime(param.GetNHits());
582 :
583 :
584 1621130 : for (i=0; i <= fTempInfo->GetStackIndex(); i++){
585 794414 : Double_t * position = fTempInfo->GetPosition(i);
586 794414 : Double_t dr = position[0]-oldr;
587 : Double_t ratio;
588 : {
589 : //Double_t dfi2 = param.fAn+2*param.fAd*(position[0]-param.fR);
590 794414 : Double_t dfi2 = param.GetAn();
591 794414 : dfi2*=dfi2*oldr*oldr;
592 : //Double_t ddz2 = param.fTheta+2*param.fThetaD*(position[0]-param.fR);
593 794414 : Double_t ddz2 = param.GetTheta();
594 794414 : ddz2*=ddz2;
595 794414 : ratio = TMath::Sqrt(1.+ dfi2+ ddz2);
596 : }
597 :
598 : // Double_t dl = fStep*(Short_t)TMath::Nint(dr*ratio/fStep); //MI change
599 794414 : Double_t dl = dr*ratio/fStep;
600 794417 : if (TMath::Abs(dl)>32765) dl =0;
601 794414 : dl = fStep * Short_t(TMath::Nint(dl));
602 :
603 794414 : dr = dl/ratio;
604 794414 : oldr+=dr;
605 : //calculate precission
606 794414 : AliTrackHitsParamV2 ¶ml = *(fTempInfo->GetParam());
607 : //real deltas
608 794414 : Double_t dr1= position[0]-paraml.GetR();
609 794414 : Double_t dz = position[1]-paraml.GetZ();
610 794414 : Double_t dfi = position[2]-paraml.GetFi();
611 : //extrapolated deltas
612 794414 : Double_t dr2 = oldr-paraml.GetR();
613 794414 : Double_t ddr = dr2-dr1;
614 794414 : Double_t ddz = dr2*paraml.GetTheta()+dr2*dr2*paraml.GetThetaD()-dz;
615 794414 : Double_t ddfi= dr2*paraml.GetAn()+dr2*dr2*paraml.GetAd()-dfi;
616 794414 : dd = TMath::Sqrt(ddz*ddz+oldr*oldr*ddfi*ddfi+ddr*ddr);
617 :
618 :
619 794414 : if ( (dd>fPrecision) ){
620 : //if ( (dd<0) ){
621 54321 : if (i==0){
622 1 : paraml.SetAn(0);
623 1 : paraml.SetAd(0);
624 1 : paraml.SetTheta(0);
625 1 : paraml.SetThetaD(0);
626 1 : Double_t ddz1 = dr2*paraml.GetTheta()+dr2*dr2*paraml.GetThetaD()-dz;
627 1 : Double_t ddfi1= dr2*paraml.GetAn()+dr2*dr2*paraml.GetAd()-dfi;
628 : dl = 0;
629 1 : dd = TMath::Sqrt(ddz1*ddz1+oldr*oldr*ddfi1*ddfi1+ddr*ddr);
630 : }
631 : else
632 54320 : break;
633 1 : }
634 :
635 740094 : paraml.HitDistance(i)= Short_t(TMath::Nint(dl/fStep));
636 740094 : paraml.Charge(i)= Short_t(fTempInfo->GetQStack(i));
637 740094 : paraml.Time(i)= Short_t(fTempInfo->GetTimeStack(i)/AliTPCTrackHitsV2::fgkTimePrecision);
638 740094 : }
639 :
640 70471 : 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 54320 : param.SetNHits(i);
651 54320 : param.ResizeCharge(i);
652 54320 : param.ResizeTime(i);
653 54320 : param.ResizeHitDistance(i);
654 : //
655 54320 : Int_t volumeID = fTempInfo->GetParam()->GetVolumeID();
656 54320 : Int_t trackID =fTempInfo->GetParam()->GetTrackID();
657 54320 : fTempInfo->SetParamIndex(fTempInfo->GetParamIndex()+1);
658 108640 : fTempInfo->SetParam(new((*fArray)[fTempInfo->GetParamIndex()]) AliTrackHitsParamV2);
659 54320 : Double_t * p = fTempInfo->GetPosition(i);
660 54320 : UInt_t index2 = fTempInfo->GetStackIndex();
661 54320 : fTempInfo->NewParam(p[0],p[1],p[2],fTempInfo->GetQStack(i),fTempInfo->GetTimeStack(i));
662 54320 : fTempInfo->GetParam()->SetVolumeID(volumeID);
663 54320 : fTempInfo->GetParam()->SetTrackID(trackID);
664 107875 : if (i+1<=index2) FlushHitStack2(i+1,index2);
665 :
666 54997 : if (force) return FlushHitStack(kTRUE);
667 53643 : return kFALSE;
668 : }
669 16151 : return kTRUE;
670 70522 : }
671 :
672 :
673 : void AliTPCTrackHitsV2::FlushHitStack2(Int_t index1, Int_t index2)
674 : {
675 : //
676 : // second iteration flush stack
677 : // call only for hits where first iteration were not succesfully interpolated
678 : //
679 107110 : Double_t * positionstack = new Double_t[3*(index2-index1+1)];
680 53555 : UInt_t * qstack = new UInt_t[index2-index1+1];
681 53555 : Float_t * timestack = new Float_t[index2-index1+1];
682 107110 : memcpy(positionstack, fTempInfo->GetPosition(index1),
683 53555 : (3*(index2-index1+1))*sizeof(Double_t));
684 53555 : memcpy(qstack, fTempInfo->GetQStackP(index1),(index2-index1+1)*sizeof(UInt_t));
685 53555 : memcpy(timestack, fTempInfo->GetTimeStackP(index1),(index2-index1+1)*sizeof(Float_t));
686 : Double_t *p = positionstack;
687 2967068 : for (Int_t j=0; j<=index2-index1;j++){
688 1429979 : fTempInfo->SetStackIndex(fTempInfo->GetStackIndex()+1);
689 1429979 : fTempInfo->SetHit(p[3*j+0],p[3*j+1],p[3*j+2],qstack[j],timestack[j]);
690 : }
691 107110 : delete []positionstack;
692 107110 : delete []qstack;
693 107110 : delete []timestack;
694 53555 : }
695 :
696 :
697 : void AliTPCTrackHitsV2::AddVolume(Int_t volume)
698 : {
699 : //
700 : //add volumes to tthe list of volumes
701 : //
702 718 : Int_t * volumes = new Int_t[fNVolumes+1];
703 657 : if (fVolumes) memcpy(volumes,fVolumes,(fNVolumes)*sizeof(Int_t));
704 359 : volumes[fNVolumes]=volume;
705 359 : fNVolumes++;
706 955 : if (fVolumes) delete []fVolumes;
707 359 : fVolumes = volumes;
708 359 : }
709 :
710 :
711 : Bool_t AliTPCTrackHitsV2::First()
712 : {
713 : //
714 : //set Current hit for the first hit
715 : //
716 :
717 642 : if (fArray->GetSize()<=0) {
718 0 : fCurrentHit->SetStatus(kFALSE);
719 0 : return kFALSE;
720 : }
721 :
722 321 : AliTrackHitsParamV2 *param = (AliTrackHitsParamV2 *)fArray->At(0);
723 325 : if (!fHit) fHit = new AliTPChit;
724 642 : if (!(param) ) {
725 51 : fCurrentHit->SetStatus(kFALSE);
726 51 : return kFALSE;
727 : }
728 : //
729 591 : fCurrentHit->SetParamIndex(0);
730 270 : fCurrentHit->SetStackIndex(0);
731 : //
732 : //
733 270 : ((AliTPChit*)fHit)->fSector = param->GetVolumeID();
734 270 : ((AliTPChit*)fHit)->SetTrack(param->GetTrackID());
735 270 : ((AliTPChit*)fHit)->SetX(param->GetR()*TMath::Cos(param->GetFi()));
736 270 : ((AliTPChit*)fHit)->SetY(param->GetR()*TMath::Sin(param->GetFi()));
737 270 : ((AliTPChit*)fHit)->SetZ(param->GetZ());
738 270 : ((AliTPChit*)fHit)->fQ = param->Charge(0);
739 270 : ((AliTPChit*)fHit)->fTime = (Float_t)(param->Time(0)*AliTPCTrackHitsV2::fgkTimePrecision);
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];
747 : fCurrentHit->fHit.fTime = (Float_t)(param->fTime[0]*AliTPCTrackHitsV2::fgkTimePrecision);
748 : */
749 270 : fCurrentHit->SetR(param->GetR());
750 :
751 270 : fCurrentHit->SetStatus(kTRUE);
752 270 : return fCurrentHit->GetStatus();
753 321 : }
754 :
755 : Bool_t AliTPCTrackHitsV2::Next()
756 : {
757 : //
758 : // Hit iterator
759 : //
760 14999712 : if (!(fCurrentHit->GetStatus()))
761 0 : return kFALSE;
762 :
763 7499856 : fCurrentHit->SetStackIndex(fCurrentHit->GetStackIndex()+1);
764 :
765 7499856 : AliTrackHitsParamV2 *param = (AliTrackHitsParamV2 *)fArray->At(fCurrentHit->GetParamIndex());
766 7499856 : if (fCurrentHit->GetStackIndex()>=param->GetNHits()){
767 736732 : fCurrentHit->SetParamIndex(fCurrentHit->GetParamIndex()+1);
768 736732 : if (fCurrentHit->GetParamIndex()>=fArray->GetEntriesFast()){
769 270 : fCurrentHit->SetStatus(kFALSE);
770 270 : return kFALSE;
771 : }
772 736462 : param = (AliTrackHitsParamV2 *)fArray->At(fCurrentHit->GetParamIndex());
773 736462 : fCurrentHit->SetStackIndex(0);
774 736462 : fCurrentHit->SetR(param->GetR());
775 736462 : }
776 :
777 :
778 :
779 : Double_t ratio;
780 : {
781 : // Double_t dfi2 = param->fAn+2*param->fAd*(fCurrentHit->fR-param->fR);
782 7499586 : Double_t dfi2 = param->GetAn();
783 7499586 : dfi2*=dfi2*fCurrentHit->GetR()*fCurrentHit->GetR();
784 : // Double_t ddz2 = param->fTheta+2*param->fThetaD*(fCurrentHit->fR-param->fR);
785 7499586 : Double_t ddz2 = param->GetTheta();
786 7499586 : ddz2*=ddz2;
787 7499586 : ratio = TMath::Sqrt(1.+ dfi2+ ddz2);
788 : }
789 :
790 7499586 : fCurrentHit->SetR(fCurrentHit->GetR()+fStep*param->HitDistance(fCurrentHit->GetStackIndex())/ratio);
791 :
792 7499586 : Double_t dR = fCurrentHit->GetR() - param->GetR();
793 7499586 : Double_t fi = param->GetFi() + (param->GetAn()*dR+param->GetAd()*dR*dR);
794 7499586 : Double_t z = param->GetZ() + (param->GetTheta()*dR+param->GetThetaD()*dR*dR);
795 : /*
796 : fCurrentHit->fHit.fQ = param->fCharge[fCurrentHit->fStackIndex];
797 : fCurrentHit->fHit.fTime = (Float_t)(param->fTime[fCurrentHit->fStackIndex]*AliTPCTrackHitsV2::fgkTimePrecision);
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);
803 : */
804 7499586 : ((AliTPChit*)fHit)->fQ = param->Charge(fCurrentHit->GetStackIndex());
805 7499586 : ((AliTPChit*)fHit)->fTime = (Float_t)(param->Time(fCurrentHit->GetStackIndex())*AliTPCTrackHitsV2::fgkTimePrecision);
806 7499586 : ((AliTPChit*)fHit)->SetX(fCurrentHit->GetR()*TMath::Cos(fi));
807 7499586 : ((AliTPChit*)fHit)->SetY(fCurrentHit->GetR()*TMath::Sin(fi));
808 7499586 : ((AliTPChit*)fHit)->SetZ(z);
809 7499586 : ((AliTPChit*)fHit)->fSector = param->GetVolumeID();
810 7499586 : ((AliTPChit*)fHit)->SetTrack(param->GetTrackID());
811 :
812 : return kTRUE;
813 7499856 : }
814 :
815 : AliHit * AliTPCTrackHitsV2::GetHit() const
816 : {
817 : //
818 : // Return one hit
819 : //
820 30000387 : return (fCurrentHit->GetStatus())? fHit:0;
821 : //return &fCurrentHit->fHit;
822 :
823 : }
824 :
825 : AliTrackHitsParamV2 * AliTPCTrackHitsV2::GetParam()
826 : {
827 : //
828 : // Return current parameters
829 : //
830 0 : return (fCurrentHit->GetStatus())?
831 0 : (AliTrackHitsParamV2 *)fArray->At(fCurrentHit->GetParamIndex()):0;
832 : }
833 :
|
