]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPC.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / TPC / AliTPC.cxx
CommitLineData
fe4da5cc 1///////////////////////////////////////////////////////////////////////////////
2// //
3// Time Projection Chamber //
4// This class contains the basic functions for the Time Projection Chamber //
5// detector. Functions specific to one particular geometry are //
6// contained in the derived classes //
7// //
8//Begin_Html
9/*
10<img src="gif/AliTPCClass.gif">
11*/
12//End_Html
13// //
14// //
15///////////////////////////////////////////////////////////////////////////////
16
17#include <TMath.h>
18#include <TRandom.h>
19#include <TVector.h>
20#include <TGeometry.h>
21#include <TNode.h>
22#include <TTUBS.h>
23#include <TObjectTable.h>
24#include "GParticle.h"
25#include "AliTPC.h"
26#include "AliRun.h"
27#include <iostream.h>
28#include <fstream.h>
29#include "AliMC.h"
30
31ClassImp(AliTPC)
32
33//_____________________________________________________________________________
34AliTPC::AliTPC()
35{
36 //
37 // Default constructor
38 //
39 fIshunt = 0;
40 fClusters = 0;
41 fHits = 0;
42 fDigits = 0;
43 fTracks = 0;
44 fNsectors = 0;
45 fNtracks = 0;
46 fNclusters= 0;
47}
48
49//_____________________________________________________________________________
50AliTPC::AliTPC(const char *name, const char *title)
51 : AliDetector(name,title)
52{
53 //
54 // Standard constructor
55 //
56
57 //
58 // Initialise arrays of hits and digits
59 fHits = new TClonesArray("AliTPChit", 176);
60 fDigits = new TClonesArray("AliTPCdigit",10000);
61 //
62 // Initialise counters
63 fClusters = 0;
64 fTracks = 0;
65 fNsectors = 72;
66 fNtracks = 0;
67 fNclusters= 0;
68 fDigitsIndex = new Int_t[fNsectors+1];
69 fClustersIndex = new Int_t[fNsectors+1];
70 //
71 fIshunt = 0;
72 //
73 // Initialise color attributes
74 SetMarkerColor(kYellow);
75}
76
77//_____________________________________________________________________________
78AliTPC::~AliTPC()
79{
80 //
81 // TPC destructor
82 //
83 fIshunt = 0;
84 delete fHits;
85 delete fDigits;
86 delete fClusters;
87 delete fTracks;
88 if (fDigitsIndex) delete [] fDigitsIndex;
89 if (fClustersIndex) delete [] fClustersIndex;
90}
91
92//_____________________________________________________________________________
93void AliTPC::AddCluster(Float_t *hits, Int_t *tracks)
94{
95 //
96 // Add a simulated cluster to the list
97 //
98 if(!fClusters) fClusters=new TClonesArray("AliTPCcluster",10000);
99 TClonesArray &lclusters = *fClusters;
100 new(lclusters[fNclusters++]) AliTPCcluster(hits,tracks);
101}
102
103//_____________________________________________________________________________
104void AliTPC::AddCluster(const AliTPCcluster &c)
105{
106 //
107 // Add a simulated cluster copy to the list
108 //
109 if(!fClusters) fClusters=new TClonesArray("AliTPCcluster",10000);
110 TClonesArray &lclusters = *fClusters;
111 new(lclusters[fNclusters++]) AliTPCcluster(c);
112}
113
114//_____________________________________________________________________________
115void AliTPC::AddDigit(Int_t *tracks, Int_t *digits)
116{
117 //
118 // Add a TPC digit to the list
119 //
120 TClonesArray &ldigits = *fDigits;
121 new(ldigits[fNdigits++]) AliTPCdigit(tracks,digits);
122}
123
124//_____________________________________________________________________________
125void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
126{
127 //
128 // Add a hit to the list
129 //
130 TClonesArray &lhits = *fHits;
131 new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
132}
133
134//_____________________________________________________________________________
135void AliTPC::AddTrack(Float_t *hits)
136{
137 //
138 // Add a track to the list of tracks
139 //
140 TClonesArray &ltracks = *fTracks;
141 new(ltracks[fNtracks++]) AliTPCtrack(hits);
142}
143
144//_____________________________________________________________________________
145void AliTPC::AddTrack(const AliTPCtrack& t)
146{
147 //
148 // Add a track copy to the list of tracks
149 //
150 if(!fTracks) fTracks=new TClonesArray("AliTPCtrack",10000);
151 TClonesArray &ltracks = *fTracks;
152 new(ltracks[fNtracks++]) AliTPCtrack(t);
153}
154
155//_____________________________________________________________________________
156void AliTPC::BuildGeometry()
157{
158 //
159 // Build TPC ROOT TNode geometry for the event display
160 //
161 TNode *Node, *Top;
162 TTUBS *tubs;
163 Int_t i;
164 const int kColorTPC=19;
165 char name[5], title[18];
166 const Double_t kDegrad=TMath::Pi()/180;
167 const Double_t loAng=30;
168 const Double_t hiAng=15;
169 const Int_t nLo = Int_t (360/loAng+0.5);
170 const Int_t nHi = Int_t (360/hiAng+0.5);
171 const Double_t loCorr = 1/TMath::Cos(0.5*loAng*kDegrad);
172 const Double_t hiCorr = 1/TMath::Cos(0.5*hiAng*kDegrad);
173 //
174 // Get ALICE top node
175 Top=gAlice->GetGeometry()->GetNode("alice");
176 //
177 // Inner sectors
178 for(i=0;i<nLo;i++) {
179 sprintf(name,"LS%2.2d",i);
180 sprintf(title,"TPC low sector %d",i);
181 tubs = new TTUBS(name,title,"void",88*loCorr,136*loCorr,250,loAng*(i-0.5),loAng*(i+0.5));
182 tubs->SetNumberOfDivisions(1);
183 Top->cd();
184 Node = new TNode(name,title,name,0,0,0,"");
185 Node->SetLineColor(kColorTPC);
186 fNodes->Add(Node);
187 }
188 // Outer sectors
189 for(i=0;i<nHi;i++) {
190 sprintf(name,"US%2.2d",i);
191 sprintf(title,"TPC upper sector %d",i);
192 tubs = new TTUBS(name,title,"void",142*hiCorr,250*hiCorr,250,hiAng*(i-0.5),hiAng*(i+0.5));
193 tubs->SetNumberOfDivisions(1);
194 Top->cd();
195 Node = new TNode(name,title,name,0,0,0,"");
196 Node->SetLineColor(kColorTPC);
197 fNodes->Add(Node);
198 }
199}
200
201
202//_____________________________________________________________________________
203void AliTPC::CreateList(Int_t *tracks,Float_t signal[][MAXTPCTBK+1],
204 Int_t ntr,Int_t time)
205{
206 //
207 // Creates list of tracks contributing to a given digit
208 // Only the 3 most significant tracks are taken into account
209 //
210
211 Int_t i,j;
212
213 for(i=0;i<3;i++) tracks[i]=-1;
214
215 //
216 // Loop over signals, only 3 times
217 //
218
219 Float_t qmax;
220 Int_t jmax;
221 Int_t jout[3] = {-1,-1,-1};
222
223 for(i=0;i<3;i++){
224 qmax=0.;
225 jmax=0;
226
227 for(j=0;j<ntr;j++){
228
229 if((i == 1 && j == jout[i-1])
230 ||(i == 2 && (j == jout[i-1] || j == jout[i-2]))) continue;
231
232 if(signal[j][time] > qmax) {
233 qmax = signal[j][time];
234 jmax=j;
235 }
236 }
237
238 if(qmax > 0.) {
239 tracks[i]=jmax;
240 jout[i]=jmax;
241 }
242
243 }
244
245 for(i=0;i<3;i++){
246 if(tracks[i] < 0){
247 tracks[i]=0;
248 }
249 else {
250 tracks[i]=(Int_t) signal[tracks[i]][0]; // first element is a track number
251 }
252 }
253}
254
255//_____________________________________________________________________________
256void AliTPC::DigSignal(Int_t isec,Int_t irow,TObjArray *pointer)
257{
258 //
259 // Digitalise TPC signal
260 //
261 Int_t pad_c,n_of_pads;
262 Int_t pad_number;
263
264 n_of_pads = (isec < 25) ? npads_low[irow] : npads_up[irow];
265 pad_c=(n_of_pads+1)/2; // this is the "central" pad for a row
266
267 Int_t track,idx;
268 Int_t entries;
269 TVector *pp;
270 TVector *ppp;
271 Float_t y,yy,z;
272 Float_t pad_signal = 0;
273 Float_t signal[MAXTPCTBK]; // Integrated signal over all tracks
274 Float_t signal_tr[100][MAXTPCTBK+1]; // contribution from less than 50 tracks
275 Int_t flag; // flag indicating a track contributing to a pad signal
276
277 //
278
279 Int_t ntracks = pointer->GetEntriesFast();
280
281 if(ntracks == 0) return; // no signal on this row!
282
283 //--------------------------------------------------------------
284 // For each electron calculate the pad number and the avalanche
285 // This is only once per pad row
286 //--------------------------------------------------------------
287
288 TVector **el = new TVector* [ntracks]; // each track is a new vector
289
290 TObjArray *arr = new TObjArray; // array of tracks for this row
291
292 for(track=0;track<ntracks;track++) {
293 pp = (TVector*) pointer->At(track);
294 entries = pp->GetNrows();
295 el[track] = new TVector(entries-1);
296 TVector &v1 = *el[track];
297 TVector &v2 = *pp;
298
299 for(idx=0;idx<entries-2;idx+=2)
300 {
301 y=v2(idx+1);
302 yy=TMath::Abs(y);
303
304 Float_t range=((n_of_pads-1)/2 + 0.5)*pad_pitch_w;
305 //
306 // Pad number and pad range
307 //
308 if(yy < 0.5*pad_pitch_w){
309 pad_number=pad_c;
310 }
311 else if (yy < range){
312 pad_number=(Int_t) ((yy-0.5*pad_pitch_w)/pad_pitch_w +1.);
313 pad_number=(Int_t) (pad_number*TMath::Sign(1.,(double) y)+pad_c);
314 }
315 else{
316 pad_number=0;
317 }
318
319 v1(idx) = (Float_t) pad_number;
320
321 // Avalanche, taking the fluctuations into account
322
323 Int_t gain_fluct = (Int_t) (-gas_gain*TMath::Log(gRandom->Rndm()));
324 v1(idx+1)= (Float_t) gain_fluct;
325
326 } // end of loop over electrons
327
328 arr->Add(el[track]); // add the vector to the array
329
330 }// end of loop over tracks
331
332 delete [] el; // delete an array of pointers
333
334 //-------------------------------------------------------------
335 // Calculation of signal for every pad
336 //-------------------------------------------------------------
337
338 //-------------------------------------------------------------
339 // Loop over pads
340 //-------------------------------------------------------------
341
342
343 for(Int_t np=1;np<n_of_pads+1;np++)
344 {
345 for(Int_t l =0;l<MAXTPCTBK;l++) signal[l]=0.; // set signals for this pad to 0
346
347 for(Int_t k1=0;k1<100;k1++){
348 for(Int_t k2=0;k2<MAXTPCTBK+1;k2++) signal_tr[k1][k2]=0.;
349 }
350 Int_t track_counter=0;
351 //
352
353 //---------------------------------------------------------
354 // Loop over tracks
355 // --------------------------------------------------------
356
357 for(track=0;track<ntracks;track++)
358 {
359 flag = 0;
360 pp = (TVector*) pointer->At(track);
361 ppp = (TVector*) arr->At(track);
362
363 TVector &v1 = *pp;
364 TVector &v2 = *ppp;
365
366 entries = pp->GetNrows();
367
368
369 //----------------------------------------------------
370 // Loop over electrons
371 //----------------------------------------------------
372
373 for(idx=0;idx<entries-2;idx+=2)
374 {
375
376 pad_number = (Int_t) v2(idx);
377
378 if(pad_number == 0) continue; // skip electrons outside range
379
380 Int_t pad_dist = pad_number-np;
381 Int_t abs_dist = TMath::Abs(pad_dist);
382
383 if(abs_dist > 3) continue; // beyond signal range
384
385 y= v1(idx+1);
386 z = v1(idx+2);
387
388 Float_t dist = y-(pad_number-pad_c)*pad_pitch_w;
389
390 //----------------------------------------------
391 // Calculate the signal induced on a pad "np"
392 //----------------------------------------------
393
394 if(pad_dist < 0) dist = -dist;
395
396 switch((Int_t) abs_dist){
397 case 0 : pad_signal = P4(dist); // electron is on pad "np"
398 break;
399 case 1 : pad_signal = P3(dist); // electron is 1 pad away
400 break;
401 case 2 : pad_signal = P2(dist); // electron is 2 pads away
402 break;
403 case 3 : pad_signal = P1(dist); // electron is 3 pads away
404 }
405
406 //---------------------------------
407 // Multiply by a gas gain
408 //---------------------------------
409
410 pad_signal=pad_signal*v2(idx+1);
411
412 flag = 1;
413
414
415 //-----------------------------------------------
416 // Sample the pad signal in time
417 //-----------------------------------------------
418
419 Float_t t_drift = (z_end-TMath::Abs(z))/v_drift; // drift time
420
421 Float_t t_offset = t_drift-t_sample*(Int_t)(t_drift/t_sample);
422 Int_t first_bucket = (Int_t) (t_drift/t_sample+1.);
423
424 for(Int_t bucket = 1;bucket<6;bucket++){
425 Float_t time = (bucket-1)*t_sample+t_offset;
426 Int_t time_idx = first_bucket+bucket-1;
427 Float_t ampl = pad_signal*TimeRes(time);
428 if (time_idx > MAXTPCTBK) break; //Y.Belikov
429 if (track_counter >=100) break; //Y.Belikov
430
431 signal_tr[track_counter][time_idx] += ampl; // single track only
432 signal[time_idx-1] += ampl; // fill a signal array for this pad
433
434 } // end of time sampling
435
436 } // end of loop over electrons for a current track
437
438 //-----------------------------------------------
439 // add the track number and update the counter
440 // if it gave a nonzero contribution to the pad
441 //-----------------------------------------------
442 if(flag != 0 && track_counter < 100){
443 signal_tr[track_counter][0] = v1(0);
444 track_counter++; // counter is looking at the NEXT track!
445 }
446
447 } // end of loop over tracks
448
449 //----------------------------------------------
450 // Fill the Digits for this pad
451 //----------------------------------------------
452
453 Int_t tracks[3];
454 Int_t digits[5];
455
456 digits[0]=isec; // sector number
457 digits[1]=irow+1; // row number
458 digits[2]=np; // pad number
459
460 Float_t q;
461
462 for(Int_t time = 0;time<MAXTPCTBK;time++){
463 digits[3] = time+1; // time bucket
464
465 q = signal[time];
466 q = gRandom->Gaus(q,sigma_noise); // apply noise
467
468 q *= (q_el*1.e15); // convert to fC
469 q *= chip_gain; // convert to mV
470 q *= (adc_sat/dyn_range); // convert to ADC counts
471
472 if(q < zero_sup) continue; // do not fill "zeros"
473 if(q > adc_sat) q = adc_sat; // saturation
474 digits[4] = (Int_t) q; // ADC counts
475
476 //--------------------------------------------------
477 // "Real signal" or electronics noise
478 //--------------------------------------------------
479
480 if(signal[time] > 0.){
481
482 //-----------------------------------------------------
483 // Create a list of tracks contributing to this digit
484 // If the digit results from a noise, track number is 0
485 //-----------------------------------------------------
486
487 CreateList(tracks,signal_tr,track_counter,time);
488 }
489 else {
490 for(Int_t ii=0;ii<3;ii++) tracks[ii]=0;
491 }
492
493 AddDigit(tracks,digits);
494
495
496 } // end of digits for this pad
497
498 } // end of loop over pads
499
500 arr->Delete(); // delete objects in this array
501
502 delete arr; // deletes the TObjArray itselves
503
504}
505
506//_____________________________________________________________________________
507Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t )
508{
509 //
510 // Calculate distance from TPC to mouse on the display
511 // Dummy procedure
512 //
513 return 9999;
514}
515
516//_____________________________________________________________________________
517const int MAX_CLUSTER=nrow_low+nrow_up;
518const int S_MAXSEC=24;
519const int L_MAXSEC=48;
520const int ROWS_TO_SKIP=21;
521const Float_t MAX_CHI2=15.;
522const Float_t THRESHOLD=8*zero_sup;
523
524//_____________________________________________________________________________
525static Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
526{
527 //
528 // Calculate spread in Y
529 //
530 pt=TMath::Abs(pt)*1000.;
531 Double_t x=r/pt;
532 tgl=TMath::Abs(tgl);
533 Double_t s=a_rphi - b_rphi*r*tgl + c_rphi*x*x + d_rphi*x;
534 if (s<0.4e-3) s=0.4e-3;
535 return s;
536}
537
538//_____________________________________________________________________________
539static Double_t SigmaZ2(Double_t r, Double_t tgl)
540{
541 //
542 // Calculate spread in Z
543 //
544 tgl=TMath::Abs(tgl);
545 Double_t s=a_z - b_z*r*tgl + c_z*tgl*tgl;
546 if (s<0.4e-3) s=0.4e-3;
547 return s;
548}
549
550//_____________________________________________________________________________
551inline Double_t f1(Double_t x1,Double_t y1, //C
552 Double_t x2,Double_t y2,
553 Double_t x3,Double_t y3)
554{
555 //
556 // Function f1
557 //
558 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
559 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
560 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
561 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
562 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
563 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
564
565 return -xr*yr/sqrt(xr*xr+yr*yr);
566}
567
568
569//_____________________________________________________________________________
570inline Double_t f2(Double_t x1,Double_t y1, //eta=C*x0
571 Double_t x2,Double_t y2,
572 Double_t x3,Double_t y3)
573{
574 //
575 // Function f2
576 //
577 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
578 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
579 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
580 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
581 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
582 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
583
584 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
585}
586
587//_____________________________________________________________________________
588inline Double_t f3(Double_t x1,Double_t y1, //tgl
589 Double_t x2,Double_t y2,
590 Double_t z1,Double_t z2)
591{
592 //
593 // Function f3
594 //
595 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
596}
597
598//_____________________________________________________________________________
599static int FindProlongation(AliTPCtrack& t, const AliTPCSector *sec,
600 int s, int ri, int rf=0)
601{
602 //
603 // Propagate track
604 //
605 int try_again=ROWS_TO_SKIP;
606 Double_t alpha=sec->GetAlpha();
607 int ns=int(2*TMath::Pi()/alpha)+1;
608 for (int nr=ri; nr>=rf; nr--) {
609 Double_t x=sec[s].GetX(nr), ymax=sec[s].GetMaxY(nr);
610 if (!t.PropagateTo(x)) break;
611
612 const AliTPCcluster *cl=0;
613 Double_t max_chi2=MAX_CHI2;
614 const AliTPCRow& row=sec[s][nr];
615 Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
616 Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
617 Double_t road=3.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
618
619 if (road>30) {
620 if (t>3) cerr<<t<<" AliTPCtrack warning: Too broad road !\n";
621 break;
622 }
623
624 if (row)
625 for (int i=row.Find(y-road); i<row; i++) {
626 AliTPCcluster* c=(AliTPCcluster*)(row[i]);
627 if (c->fY > y+road) break;
628 if (c->IsUsed()) continue;
629 if ((c->fZ - z)*(c->fZ - z) > 9.*(t.GetSigmaZ2() + sz2)) continue;
630 Double_t chi2=t.GetPredictedChi2(c);
631 if (chi2 > max_chi2) continue;
632 max_chi2=chi2;
633 cl=c;
634 }
635 if (cl) {
636 t.Update(cl,max_chi2);
637 try_again=ROWS_TO_SKIP;
638 } else {
639 if (try_again--) {
640 if (y > ymax) {
641 s = (s+1) % ns;
642 if (!t.Rotate(alpha)) break;
643 } else
644 if (y <-ymax) {
645 s = (s-1+ns) % ns;
646 if (!t.Rotate(-alpha)) break;
647 };
648 continue;
649 }
650 break;
651 }
652 }
653 return s;
654}
655
656//_____________________________________________________________________________
657static void MakeSeeds(TObjArray& seeds,const AliTPCSector* sec,int i1,int i2)
658{
659 //
660 // Find seed for tracking
661 //
662 TMatrix C(5,5); TVector x(5);
663 int max_sec=L_MAXSEC/2;
664 for (int ns=0; ns<max_sec; ns++) {
665 int nl=sec[(ns-1+max_sec)%max_sec][i2];
666 int nm=sec[ns][i2];
667 int nu=sec[(ns+1)%max_sec][i2];
668 Double_t alpha=sec[ns].GetAlpha();
669 const AliTPCRow& r1=sec[ns][i1];
670 for (int is=0; is < r1; is++) {
671 Double_t x1=sec[ns].GetX(i1), y1=r1[is]->fY, z1=r1[is]->fZ;
672 for (int js=0; js < nl+nm+nu; js++) {
673 const AliTPCcluster *cl;
674 Double_t cs,sn;
675 //int ks;
676
677 if (js<nl) {
678 //ks=(ns-1+max_sec)%max_sec;
679 const AliTPCRow& r2=sec[(ns-1+max_sec)%max_sec][i2];
680 cl=r2[js];
681 cs=cos(alpha); sn=sin(alpha);
682 } else
683 if (js<nl+nm) {
684 //ks=ns;
685 const AliTPCRow& r2=sec[ns][i2];
686 cl=r2[js-nl];
687 cs=1; sn=0.;
688 } else {
689 //ks=(ns+1)%max_sec;
690 const AliTPCRow& r2=sec[(ns+1)%max_sec][i2];
691 cl=r2[js-nl-nm];
692 cs=cos(alpha); sn=-sin(alpha);
693 }
694 Double_t x2=sec[ns].GetX(i2), y2=cl->fY, z2=cl->fZ;
695 //if (z1*z2 < 0) continue;
696 //if (TMath::Abs(z1) < TMath::Abs(z2)) continue;
697
698 Double_t tmp= x2*cs+y2*sn;
699 y2 =-x2*sn+y2*cs;
700 x2=tmp;
701
702 x(0)=y1;
703 x(1)=z1;
704 x(2)=f1(x1,y1,x2,y2,0.,0.);
705 x(3)=f2(x1,y1,x2,y2,0.,0.);
706 x(4)=f3(x1,y1,x2,y2,z1,z2);
707
708 if (TMath::Abs(x(2)*x1-x(3)) >= 0.999) continue;
709
710 if (TMath::Abs(x(4)) > 1.2) continue;
711
712 Double_t a=asin(x(3));
713 /*
714 Double_t tgl1=z1*x(2)/(a+asin(x(2)*x1-x(3)));
715 Double_t tgl2=z2*x(2)/(a+asin(x(2)*x2-x(3)));
716 Double_t ratio=2*(tgl1-tgl2)/(tgl1+tgl2);
717 if (TMath::Abs(ratio)>0.0170) continue; //or > 0.005
718 */
719 Double_t zv=z1 - x(4)/x(2)*(a+asin(x(2)*x1-x(3)));
720 if (TMath::Abs(zv)>33.) continue;
721
722
723 TMatrix X(6,6); X=0.;
724 X(0,0)=r1[is]->fSigmaY2; X(1,1)=r1[is]->fSigmaZ2;
725 X(2,2)=cl->fSigmaY2; X(3,3)=cl->fSigmaZ2;
726 X(4,4)=3./12.; X(5,5)=3./12.;
727 TMatrix F(5,6); F.UnitMatrix();
728 Double_t sy=sqrt(X(0,0)), sz=sqrt(X(1,1));
729 F(2,0)=(f1(x1,y1+sy,x2,y2,0.,0.)-x(2))/sy;
730 F(2,2)=(f1(x1,y1,x2,y2+sy,0.,0.)-x(2))/sy;
731 F(2,4)=(f1(x1,y1,x2,y2,0.,0.+sy)-x(2))/sy;
732 F(3,0)=(f2(x1,y1+sy,x2,y2,0.,0.)-x(3))/sy;
733 F(3,2)=(f2(x1,y1,x2,y2+sy,0.,0.)-x(3))/sy;
734 F(3,4)=(f2(x1,y1,x2,y2,0.,0.+sy)-x(3))/sy;
735 F(4,0)=(f3(x1,y1+sy,x2,y2,z1,z2)-x(4))/sy;
736 F(4,1)=(f3(x1,y1,x2,y2,z1+sz,z2)-x(4))/sz;
737 F(4,2)=(f3(x1,y1,x2,y2+sy,z1,z2)-x(4))/sy;
738 F(4,3)=(f3(x1,y1,x2,y2,z1,z2+sz)-x(4))/sz;
739 F(4,4)=0;
740 F(3,3)=0;
741
742 TMatrix t(F,TMatrix::kMult,X);
743 C.Mult(t,TMatrix(TMatrix::kTransposed,F));
744
745 TrackSeed *track=new TrackSeed(*(r1[is]),x,C);
746 FindProlongation(*track,sec,ns,i1-1,i2);
747 int ii=(i1-i2)/2;
748 if (*track >= ii) {seeds.AddLast(track); continue;}
749 else delete track;
750 }
751 }
752 }
753}
754
755//_____________________________________________________________________________
756void AliTPC::Clusters2Tracks()
757{
758 //
759 // TPC Track finder from clusters.
760 //
761 if (!fClusters) return;
762 AliTPCSSector ssec[S_MAXSEC/2];
763 AliTPCLSector lsec[L_MAXSEC/2];
764 int ncl=fClusters->GetEntriesFast();
765 while (ncl--) {
766 AliTPCcluster *c=(AliTPCcluster*)fClusters->UncheckedAt(ncl);
767 int sec=int(c->fSector), row=int(c->fPadRow);
768
769 if (sec<24) {
770 if (row<0 || row>nrow_low) {cerr<<"low !!!"<<row<<endl; continue;}
771 ssec[sec%12][row].InsertCluster(c);
772 } else {
773 if (row<0 || row>nrow_up ) {cerr<<"up !!!"<<row<<endl; continue;}
774 sec -= 24;
775 lsec[sec%24][row].InsertCluster(c);
776 }
777 }
778
779
780 TObjArray seeds(20000);
781 MakeSeeds(seeds,lsec,nrow_up-1,nrow_up-1-8);
782 MakeSeeds(seeds,lsec,nrow_up-1-4,nrow_up-1-4-8);
783
784 seeds.Sort();
785
786 int found=0;
787 int nseed=seeds.GetEntriesFast();
788
789 for (int s=0; s<nseed; s++) {
790 AliTPCtrack& t=*((AliTPCtrack*)seeds.UncheckedAt(s));
791 Double_t alpha=t.GetAlpha();
792 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
793 if (alpha < 0. ) alpha += 2.*TMath::Pi();
794 int ns=int(alpha/lsec->GetAlpha() + 0.5);
795
796 Double_t x=t.GetX();
797 int nr;
798 if (x<pad_row_up[nrow_up-1-4-7]) nr=nrow_up-1-4-8;
799 else if (x<pad_row_up[nrow_up-1-7]) nr=nrow_up-1-8;
800 else {cerr<<x<<" =x !!!\n"; continue;}
801
802 int ls=FindProlongation(t,lsec,ns,nr-1);
803
804 // if (t < 25) continue;
805
806 x=t.GetX(); alpha=lsec[ls].GetAlpha(); //
807 Double_t phi=ls*alpha + atan(t.GetY()/x); // Find S-sector
808 int ss=int(0.5*(phi/alpha+1)); //
809 alpha *= 2*(ss-0.5*ls); // and rotation angle
810 if (!t.Rotate(alpha)) continue; //
811 ss %= (S_MAXSEC/2); //
812
813 ss=FindProlongation(t,ssec,ss,nrow_low-1);
814 if (t < 30) continue;
815
816 AddTrack(t);
817 t.UseClusters();
818 cerr<<found++<<'\r';
819 }
820}
821
822//_____________________________________________________________________________
823void AliTPC::CreateMaterials()
824{
825 //
826 // Create Materials for for TPC
827 // Origin M.Kowalski
828 //
829
830 AliMC* pMC = AliMC::GetMC();
831
832 Int_t ISXFLD=gAlice->Field()->Integ();
833 Float_t SXMGMX=gAlice->Field()->Max();
834
835 Float_t absl, radl, a, d, z;
836 Float_t dg;
837 Float_t x0ne;
838 Float_t buf[1];
839 Int_t nbuf;
840
841 // --- Methane (CH4) ---
842 Float_t am[2] = { 12.,1. };
843 Float_t zm[2] = { 6.,1. };
844 Float_t wm[2] = { 1.,4. };
845 Float_t dm = 7.17e-4;
846 // --- The Neon CO2 90/10 mixture ---
847 Float_t ag[2] = { 20.18 };
848 Float_t zg[2] = { 10. };
849 Float_t wg[2] = { .8,.2 };
850 Float_t dne = 9e-4; // --- Neon density in g/cm3 ---
851 // --- Mylar (C5H4O2) ---
852 Float_t amy[3] = { 12.,1.,16. };
853 Float_t zmy[3] = { 6.,1.,8. };
854 Float_t wmy[3] = { 5.,4.,2. };
855 Float_t dmy = 1.39;
856 // --- CO2 ---
857 Float_t ac[2] = { 12.,16. };
858 Float_t zc[2] = { 6.,8. };
859 Float_t wc[2] = { 1.,2. };
860 Float_t dc = .001977;
861 // --- Carbon density and radiation length ---
862 Float_t densc = 2.265;
863 Float_t radlc = 18.8;
864 // --- Silicon ---
865 Float_t asi = 28.09;
866 Float_t zsi = 14.;
867 Float_t desi = 2.33;
868 Float_t radsi = 9.36;
869
870 // --- Define the various materials for GEANT ---
871 AliMaterial(0, "Al $", 26.98, 13., 2.7, 8.9, 37.2);
872 x0ne = 28.94 / dne;
873 AliMaterial(1, "Ne $", 20.18, 10., dne, x0ne, 999.);
874
875 // -- Methane, defined by the proportions of atoms
876
877 AliMixture(2, "Methane$", am, zm, dm, -2, wm);
878
879 // --- CO2, defined by the proportion of atoms
880
881 AliMixture(7, "CO2$", ac, zc, dc, -2, wc);
882
883 // -- Get A,Z etc. for CO2
884
885 char namate[21];
886 pMC->Gfmate((*fIdmate)[7], namate, a, z, d, radl, absl, buf, nbuf);
887
888 ag[1] = a;
889 zg[1] = z;
890 dg = dne * .9 + dc * .1;
891
892 //-------------------------------------
893
894 // -- Create Ne/CO2 90/10 mixture
895
896 AliMixture(3, "Gas-mixt $", ag, zg, dg, 2, wg);
897 AliMixture(4, "Gas-mixt $", ag, zg, dg, 2, wg);
898
899 AliMaterial(5, "G10$", 20., 10., 1.7, 19.4, 999.);
900 AliMixture(6, "Mylar$", amy, zmy, dmy, -3, wmy);
901
902 a = ac[0];
903 z = zc[0];
904 AliMaterial(8, "Carbon", a, z, densc, radlc, 999.);
905
906 AliMaterial(9, "Silicon", asi, zsi, desi, radsi, 999.);
907 AliMaterial(99, "Air$", 14.61, 7.3, .001205, 30420., 67500.);
908
909 AliMedium(400, "Al wall$", 0, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
910 AliMedium(402, "Gas mix1$", 3, 0, ISXFLD, SXMGMX, 10., .01,.1, .001, .01);
911 AliMedium(403, "Gas mix2$", 3, 0, ISXFLD, SXMGMX, 10., .01,.1, .001, .01);
912 AliMedium(404, "Gas mix3$", 4, 1, ISXFLD, SXMGMX, 10., .01,.1, .001, .01);
913 AliMedium(405, "G10 pln$", 5, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1 );
914 AliMedium(406, "Mylar $", 6, 0, ISXFLD, SXMGMX, 10., .01,.1, .001, .01);
915 AliMedium(407, "CO2 $", 7, 0, ISXFLD, SXMGMX, 10., .01,.1, .01, .01);
916 AliMedium(408, "Carbon $", 8, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1 );
917 AliMedium(409, "Silicon$", 9, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1 );
918 AliMedium(499, "Air gap$", 99, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1 );
919}
920
921//_____________________________________________________________________________
922struct Bin {
923 const AliTPCdigit *dig;
924 int idx;
925 Bin() {dig=0; idx=-1;}
926};
927
928struct PreCluster : public AliTPCcluster {
929 const AliTPCdigit* summit;
930 int idx;
931 int cut;
932 int npeaks;
933 PreCluster() : AliTPCcluster() {cut=npeaks=0;}
934};
935
936//_____________________________________________________________________________
937static void FindCluster(int i, int j, Bin bins[][MAXTPCTBK+2], PreCluster &c)
938{
939 //
940 // Find clusters
941 //
942 Bin& b=bins[i][j];
943 int q=b.dig->fSignal;
944
945 if (q<0) { // digit is at the edge of the pad row
946 q=-q;
947 c.cut=1;
948 }
949 if (b.idx >= 0 && b.idx != c.idx) {
950 c.idx=b.idx;
951 c.npeaks++;
952 }
953
954 if (q > TMath::Abs(c.summit->fSignal)) c.summit=b.dig;
955
956 c.fY += i*q;
957 c.fZ += j*q;
958 c.fSigmaY2 += i*i*q;
959 c.fSigmaZ2 += j*j*q;
960 c.fQ += q;
961
962 b.dig = 0; b.idx = c.idx;
963
964 if (bins[i-1][j].dig) FindCluster(i-1,j,bins,c);
965 if (bins[i][j-1].dig) FindCluster(i,j-1,bins,c);
966 if (bins[i+1][j].dig) FindCluster(i+1,j,bins,c);
967 if (bins[i][j+1].dig) FindCluster(i,j+1,bins,c);
968
969}
970
971//_____________________________________________________________________________
972void AliTPC::Digits2Clusters()
973{
974 //
975 // simple TPC cluster finder from digits.
976 //
977 const Int_t MAX_PAD=200+2, MAX_BUCKET=MAXTPCTBK+2;
978 const Int_t Q_min=200;//75;
979
980 TTree *t=gAlice->TreeD();
981 t->GetBranch("TPC")->SetAddress(&fDigits);
982 Int_t sectors_by_rows=(Int_t)t->GetEntries();
983
984 int ncls=0;
985
986 for (Int_t n=0; n<sectors_by_rows; n++) {
987 if (!t->GetEvent(n)) continue;
988 Bin bins[MAX_PAD][MAX_BUCKET];
989 AliTPCdigit *dig=(AliTPCdigit*)fDigits->UncheckedAt(0);
990 Int_t nsec=dig->fSector, nrow=dig->fPadRow;
991 Int_t ndigits=fDigits->GetEntriesFast();
992
993 int npads; int sign_z;
994 if (nsec<25) {
995 sign_z=(nsec<13) ? 1 : -1;
996 npads=npads_low[nrow-1];
997 } else {
998 sign_z=(nsec<49) ? 1 : -1;
999 npads=npads_up[nrow-1];
1000 }
1001
1002 int ndig;
1003 for (ndig=0; ndig<ndigits; ndig++) {
1004 dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig);
1005 int i=dig->fPad, j=dig->fTime;
1006 if (dig->fSignal >= THRESHOLD) bins[i][j].dig=dig;
1007 if (i==1 || i==npads || j==1 || j==MAXTPCTBK) dig->fSignal*=-1;
1008 }
1009
1010 int ncl=0;
1011 int i,j;
1012 for (i=1; i<MAX_PAD-1; i++) {
1013 for (j=1; j<MAX_BUCKET-1; j++) {
1014 if (bins[i][j].dig == 0) continue;
1015 PreCluster c; c.summit=bins[i][j].dig; c.idx=ncls;
1016 FindCluster(i, j, bins, c);
1017 //if (c.fQ <= Q_min) continue; //noise cluster
1018 c.fY /= c.fQ;
1019 c.fZ /= c.fQ;
1020 c.fSigmaY2 = c.fSigmaY2/c.fQ - c.fY*c.fY + 1./12.;
1021 c.fSigmaZ2 = c.fSigmaZ2/c.fQ - c.fZ*c.fZ + 1./12.;
1022 c.fSigmaY2 *= pad_pitch_w*pad_pitch_w;
1023 c.fSigmaZ2 *= z_end/MAXTPCTBK*z_end/MAXTPCTBK;
1024 c.fSigmaY2 *= 0.022*8;
1025 c.fSigmaZ2 *= 0.068*4;
1026 c.fY = (c.fY - 0.5 - 0.5*npads)*pad_pitch_w;
1027 c.fZ = z_end/MAXTPCTBK*c.fZ;
1028 c.fZ -= 3.*fwhm/2.35482*v_drift; // PASA delay
1029 c.fZ = sign_z*(z_end - c.fZ);
1030 c.fSector=nsec-1;
1031 c.fPadRow=nrow-1;
1032 c.fTracks[0]=c.summit->fTracks[0];
1033 c.fTracks[1]=c.summit->fTracks[1];
1034 c.fTracks[2]=c.summit->fTracks[2];
1035
1036 if (c.cut) {
1037 c.fSigmaY2 *= 25.;
1038 c.fSigmaZ2 *= 4.;
1039 }
1040
1041 AddCluster(c); ncls++; ncl++;
1042 }
1043 }
1044
1045 for (ndig=0; ndig<ndigits; ndig++) {
1046 dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig);
1047 if (TMath::Abs(dig->fSignal) >= THRESHOLD/3)
1048 bins[dig->fPad][dig->fTime].dig=dig;
1049 }
1050
1051 for (i=1; i<MAX_PAD-1; i++) {
1052 for (j=1; j<MAX_BUCKET-1; j++) {
1053 if (bins[i][j].dig == 0) continue;
1054 PreCluster c; c.summit=bins[i][j].dig; c.idx=ncls;
1055 FindCluster(i, j, bins, c);
1056 if (c.fQ <= Q_min) continue; //noise cluster
1057 if (c.npeaks>1) continue; //overlapped cluster
1058 c.fY /= c.fQ;
1059 c.fZ /= c.fQ;
1060 c.fSigmaY2 = c.fSigmaY2/c.fQ - c.fY*c.fY + 1./12.;
1061 c.fSigmaZ2 = c.fSigmaZ2/c.fQ - c.fZ*c.fZ + 1./12.;
1062 c.fSigmaY2 *= pad_pitch_w*pad_pitch_w;
1063 c.fSigmaZ2 *= z_end/MAXTPCTBK*z_end/MAXTPCTBK;
1064 c.fSigmaY2 *= 0.022*4;
1065 c.fSigmaZ2 *= 0.068*4;
1066 c.fY = (c.fY - 0.5 - 0.5*npads)*pad_pitch_w;
1067 c.fZ = z_end/MAXTPCTBK*c.fZ;
1068 c.fZ -= 3.*fwhm/2.35482*v_drift; // PASA delay
1069 c.fZ = sign_z*(z_end - c.fZ);
1070 c.fSector=nsec-1;
1071 c.fPadRow=nrow-1;
1072 c.fTracks[0]=c.summit->fTracks[0];
1073 c.fTracks[1]=c.summit->fTracks[1];
1074 c.fTracks[2]=c.summit->fTracks[2];
1075
1076 if (c.cut) {
1077 c.fSigmaY2 *= 25.;
1078 c.fSigmaZ2 *= 4.;
1079 }
1080
1081 if (c.npeaks==0) {AddCluster(c); ncls++; ncl++;}
1082 else {
1083 new ((*fClusters)[c.idx]) AliTPCcluster(c);
1084 }
1085 }
1086 }
1087
1088 cerr<<"sector, row, digits, clusters: "
1089 <<nsec<<' '<<nrow<<' '<<ndigits<<' '<<ncl<<" \r";
1090
1091 fDigits->Clear();
1092
1093 }
1094}
1095
1096//_____________________________________________________________________________
1097void AliTPC::ElDiff(Float_t *xyz)
1098{
1099 //
1100 // calculates the diffusion of a single electron
1101 //
1102
1103 Float_t driftl;
1104 //
1105 Float_t z0=xyz[2];
1106 driftl=z_end-TMath::Abs(xyz[2]);
1107 if(driftl<0.01) driftl=0.01;
1108 driftl=TMath::Sqrt(driftl);
1109 Float_t sig_t = driftl*diff_t;
1110 Float_t sig_l = driftl*diff_l;
1111 //
1112 xyz[0]=gRandom->Gaus(xyz[0],sig_t);
1113 xyz[1]=gRandom->Gaus(xyz[1],sig_t);
1114 xyz[2]=gRandom->Gaus(xyz[2],sig_l);
1115 //
1116 if (TMath::Abs(xyz[2])>z_end){
1117 xyz[2]=z_end*TMath::Sign(1.,(double) z0);
1118 }
1119 if(xyz[2]*z0 < 0.){
1120 xyz[2]=0.0001*TMath::Sign(1.,(double) z0);
1121 }
1122}
1123
1124//_____________________________________________________________________________
1125void AliTPC::Hits2Clusters()
1126{
1127 //
1128 // TPC simple cluster generator from hits
1129 // obtained from the TPC Fast Simulator
1130 //
1131 Float_t sigma_rphi,sigma_z,cl_rphi,cl_z;
1132 //
1133 GParticle *particle; // pointer to a given particle
1134 AliTPChit *tpcHit; // pointer to a sigle TPC hit
1135 TClonesArray *Particles; //pointer to the particle list
1136 Int_t sector,nhits;
1137 Int_t ipart;
1138 Float_t xyz[5];
1139 Float_t pl,pt,tanth,rpad,ratio;
1140 Float_t rot_angle;
1141 Float_t cph,sph;
1142
1143 //---------------------------------------------------------------
1144 // Get the access to the tracks
1145 //---------------------------------------------------------------
1146
1147 TTree *TH = gAlice->TreeH();
1148 Stat_t ntracks = TH->GetEntries();
1149
1150 //------------------------------------------------------------
1151 // Loop over all sectors (72 sectors)
1152 // Sectors 1-24 are lower sectors, 1-12 z>0, 13-24 z<0
1153 // Sectors 25-72 are upper sectors, 25-48 z>0, 49-72 z<0
1154 //
1155 // First cluster for sector 1 starts at "0"
1156 //------------------------------------------------------------
1157
1158
1159 fClustersIndex[0] = 0;
1160
1161 //
1162 for(Int_t isec=1;isec<fNsectors+1;isec++){
1163 //
1164 if(isec < 25){
1165 rot_angle = (isec < 13) ? (isec-1)*alpha_low : (isec-13)*alpha_low;
1166 }
1167 else {
1168 rot_angle = (isec < 49) ? (isec-25)*alpha_up : (isec-49)*alpha_up;
1169 }
1170
1171 cph=TMath::Cos(rot_angle);
1172 sph=TMath::Sin(rot_angle);
1173
1174 //------------------------------------------------------------
1175 // Loop over tracks
1176 //------------------------------------------------------------
1177
1178 for(Int_t track=0;track<ntracks;track++){
1179 ResetHits();
1180 TH->GetEvent(track);
1181 //
1182 // Get number of the TPC hits and a pointer
1183 // to the particles
1184 //
1185 nhits=fHits->GetEntriesFast();
1186 Particles=gAlice->Particles();
1187 //
1188 // Loop over hits
1189 //
1190 for(Int_t hit=0;hit<nhits;hit++){
1191 tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
1192 sector=tpcHit->fSector; // sector number
1193 if(sector != isec) continue; //terminate iteration
1194 ipart=tpcHit->fTrack;
1195 particle=(GParticle*)Particles->UncheckedAt(ipart);
1196 pl=particle->GetPz();
1197 pt=particle->GetPT();
1198 if(pt < 1.e-9) pt=1.e-9;
1199 tanth=pl/pt;
1200 tanth = TMath::Abs(tanth);
1201 rpad=TMath::Sqrt(tpcHit->fX*tpcHit->fX + tpcHit->fY*tpcHit->fY);
1202 ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
1203
1204 // space-point resolutions
1205
1206 sigma_rphi=SigmaY2(rpad,tanth,pt);
1207 sigma_z =SigmaZ2(rpad,tanth );
1208
1209 // cluster widths
1210
1211 cl_rphi=ac_rphi-bc_rphi*rpad*tanth+cc_rphi*ratio*ratio;
1212 cl_z=ac_z-bc_z*rpad*tanth+cc_z*tanth*tanth;
1213
1214 // temporary protection
1215
1216 if(sigma_rphi < 0.) sigma_rphi=0.4e-3;
1217 if(sigma_z < 0.) sigma_z=0.4e-3;
1218 if(cl_rphi < 0.) cl_rphi=2.5e-3;
1219 if(cl_z < 0.) cl_z=2.5e-5;
1220
1221 //
1222
1223 //
1224 // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
1225 // then the inaccuracy in a X-Y plane is only along Y (pad row)!
1226 //
1227 Float_t xprim= tpcHit->fX*cph + tpcHit->fY*sph;
1228 Float_t yprim=-tpcHit->fX*sph + tpcHit->fY*cph;
1229 xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigma_rphi)); // y
1230 xyz[1]=gRandom->Gaus(tpcHit->fZ,TMath::Sqrt(sigma_z)); // z
1231 xyz[2]=tpcHit->fQ; // q
1232 xyz[3]=sigma_rphi; // fSigmaY2
1233 xyz[4]=sigma_z; // fSigmaZ2
1234
1235 //find row number
1236 int row;
1237 if (xprim > 0.5*(pad_row_up[0]+pad_row_low[nrow_low-1])) {
1238 for (row=0; row<nrow_up; row++) if (xprim < pad_row_up[row]) break;
1239 } else {
1240 for (row=0; row<nrow_low; row++) if (xprim < pad_row_low[row]) break;
1241 }
1242
1243 // and finally add the cluster
1244 Int_t tracks[5]={tpcHit->fTrack+1, 0, 0, sector-1, row-1};
1245 AddCluster(xyz,tracks);
1246
1247 } // end of loop over hits
1248 } // end of loop over tracks
1249
1250 fClustersIndex[isec] = fNclusters; // update clusters index
1251
1252 } // end of loop over sectors
1253
1254 fClustersIndex[fNsectors]--; // set end of the clusters buffer
1255
1256}
1257
1258//_____________________________________________________________________________
1259void AliTPC::Hits2Digits()
1260{
1261 //
1262 // TPC conversion from hits to digits.
1263 //
1264
1265 Int_t i;
1266 //
1267 AliTPChit *tpcHit; // pointer to a sigle TPC hit
1268 //
1269 Float_t xyz[3];
1270 Float_t rot_angle;
1271 //-------------------------------------------------------
1272 // Get the access to the tracks
1273 //-------------------------------------------------------
1274 TTree *TH = gAlice->TreeH();
1275 Stat_t ntracks = TH->GetEntries();
1276
1277 //----------------------------------------------------
1278 // Loop over all sectors (72 sectors)
1279 // Sectors 1-24 are lower sectors, 1-12 z>0, 13-24 z<0
1280 // Sectors 25-72 are upper sectors, 25-48 z>0, 49-72 z<0
1281 //----------------------------------------------------
1282
1283 for(Int_t isec=1;isec<fNsectors+1;isec++){
1284
1285 //
1286 printf("*** Processing sector number %d ***\n",isec);
1287
1288 if(isec < 25){
1289 rot_angle = (isec < 13) ? (isec-1)*alpha_low : (isec-13)*alpha_low;
1290 }
1291 else {
1292 rot_angle = (isec < 49) ? (isec-25)*alpha_up : (isec-49)*alpha_up;
1293 }
1294
1295 Int_t nrows = (isec<25) ? 23 : 52;
1296
1297
1298 Float_t cph=TMath::Cos(rot_angle);
1299 Float_t sph=TMath::Sin(rot_angle);
1300
1301
1302
1303 //----------------------------------------------
1304 // Create TObjArray-s, one for each row
1305 //----------------------------------------------
1306
1307 TObjArray **row = new TObjArray* [nrows];
1308 for(i=0; i<nrows; i++){
1309 row[i] = new TObjArray;
1310 }
1311
1312 //----------------------------------------------
1313 // Loop over tracks
1314 //----------------------------------------------
1315 for(Int_t track=0;track<ntracks;track++){
1316 ResetHits();
1317 TH->GetEvent(track);
1318
1319 //------------------------------------------------
1320 // Get number of the TPC hits and a pointer
1321 // to the particles
1322 //------------------------------------------------
1323 Int_t nhits=fHits->GetEntriesFast();
1324 if(nhits == 0) continue;
1325 //-----------------------------------------------
1326 // Create vectors for storing the track information,
1327 // one vector per track per row,
1328 // first element is a track number and then
1329 // there are (y,z) pairs * number of electrons
1330 //----------------------------------------------
1331
1332 TVector **tr = new TVector* [nrows];
1333 Int_t *n_of_electrons= new int [nrows]; // electron counter
1334 for(i=0;i<nrows;i++){
1335 tr[i] = new TVector(241); // 120 electrons initialy
1336 n_of_electrons[i]=0;
1337 }
1338 //-----------------------------------------------------
1339 // Loop over hits
1340 //------------------------------------------------------
1341 for(Int_t hit=0;hit<nhits;hit++){
1342 tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
1343 Int_t sector=tpcHit->fSector; // sector number
1344 if(sector != isec) continue; //terminate iteration
1345
1346 xyz[0]=tpcHit->fX;
1347 xyz[1]=tpcHit->fY;
1348 xyz[2]=tpcHit->fZ;
1349 Int_t QI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
1350
1351 //-----------------------------------------------
1352 // Rotate the electron cluster to sector 1,13,25,49
1353 //-----------------------------------------------
1354 Float_t xprim=xyz[0]*cph+xyz[1]*sph;
1355 Float_t yprim=-xyz[0]*sph+xyz[1]*cph;
1356 Float_t zprim=xyz[2];
1357
1358 //-------------------------------------
1359 // Loop over electrons
1360 //-------------------------------------
1361 for(Int_t nel=0;nel<QI;nel++){
1362 xyz[0]=xprim; //
1363 xyz[1]=yprim; // Keep the initial cluster position!
1364 xyz[2]=zprim; //
1365
1366 ElDiff(xyz); // Appply the diffusion
1367
1368 Float_t row_first;
1369 Int_t row_number;
1370 row_first = (isec<25) ? pad_row_low[0] : pad_row_up[0];
1371
1372 row_number=(Int_t) ((xyz[0]-row_first+0.5*pad_pitch_l)/pad_pitch_l);
1373
1374 // Check if out of range
1375
1376 if((xyz[0]-row_first+0.5*pad_pitch_l) < 0
1377 || row_number > (nrows-1)) continue;
1378
1379 n_of_electrons[row_number]++;
1380
1381 //
1382 // Expand vector if necessary
1383 //
1384 if(n_of_electrons[row_number]>120){
1385 Int_t range = tr[row_number]->GetNrows();
1386 if(n_of_electrons[row_number] > (range-1)/2){
1387 tr[row_number]->ResizeTo(range+30); // Add 15 electrons
1388 }
1389 }
1390
1391 //---------------------------------
1392 // E x B effect at the wires
1393 //---------------------------------
1394 Int_t nw;
1395 nw=(nwires+1)/2;
1396 Float_t xx,dx;
1397 for (Int_t nwire=1;nwire<=nwires;nwire++){
1398 xx=(nwire-nw)*ww_pitch+
1399 ((isec<13) ? pad_row_low[row_number]:pad_row_up[row_number]);
1400 dx=xx-xyz[0];
1401 if(TMath::Abs(dx) < 0.5*ww_pitch) {
1402 xyz[1]=dx*omega_tau+xyz[1];
1403 break;
1404 }
1405 } // end of loop over the wires
1406
1407 TVector &v = *tr[row_number];
1408 Int_t idx = 2*n_of_electrons[row_number]-1;
1409 v(idx)=xyz[1];
1410 v(idx+1)=xyz[2];
1411 } // end of loop over electrons
1412 } // end of loop over hits
1413
1414 //
1415 // The track is finished
1416 //
1417 int trk=((AliTPChit*)fHits->UncheckedAt(0))->fTrack; //Y.Belikov
1418 for(i=0;i<nrows;i++){
1419 TVector &v = *tr[i];
1420 if(n_of_electrons[i] >0) {
1421 // v(0)=(Float_t)(track+1); // track number
1422 v(0)=(Float_t)(trk+1); // Y.Belikov
1423 tr[i]->ResizeTo(2*n_of_electrons[i]+1); // shrink if necessary
1424 row[i]->Add(tr[i]); // add to the row-array
1425 }
1426 else{
1427 delete tr[i]; // delete TVector if empty
1428 }
1429 }
1430 delete [] tr; // delete track pointers
1431 delete [] n_of_electrons; // delete n_of_electrons array
1432 } //end of loop over tracks
1433 //---------------------------------------------------
1434 // Digitize the sector data, row by row
1435 //---------------------------------------------------
1436
1437 printf("*** Digitizing sector number %d ***\n",isec);
1438
1439 for(i=0;i<nrows;i++){
1440
1441 DigSignal(isec,i,row[i]);
1442 gAlice->TreeD()->Fill();
1443 int ndig=fDigits->GetEntriesFast();
1444 printf("*** Sector, row, digits %d %d %d ***\n",isec,i+1,ndig);
1445
1446 ResetDigits();
1447
1448 row[i]->Delete(); // delete objects in this array
1449 delete row[i]; // delete the TObjArray itselves
1450
1451 } // end of digitization
1452
1453 delete [] row; // delete vectors of pointers to sectors
1454
1455 } //end of loop over sectors
1456
1457}
1458
1459//_____________________________________________________________________________
1460void AliTPC::Init()
1461{
1462 //
1463 // Initialise TPC detector after definition of geometry
1464 //
1465 Int_t i;
1466 //
1467 printf("\n");
1468 for(i=0;i<35;i++) printf("*");
1469 printf(" TPC_INIT ");
1470 for(i=0;i<35;i++) printf("*");
1471 printf("\n");
1472 //
1473 for(i=0;i<80;i++) printf("*");
1474 printf("\n");
1475}
1476
1477//_____________________________________________________________________________
1478void AliTPC::MakeBranch(Option_t* option)
1479{
1480 //
1481 // Create Tree branches for the TPC.
1482 //
1483 Int_t buffersize = 4000;
1484 char branchname[10];
1485 sprintf(branchname,"%s",GetName());
1486
1487 AliDetector::MakeBranch(option);
1488
1489 char *D = strstr(option,"D");
1490
1491 if (fDigits && gAlice->TreeD() && D) {
1492 gAlice->TreeD()->Branch(branchname,&fDigits, buffersize);
1493 printf("Making Branch %s for digits\n",branchname);
1494 }
1495
1496 char *R = strstr(option,"R");
1497
1498 if (fClusters && gAlice->TreeR() && R) {
1499 gAlice->TreeR()->Branch(branchname,&fClusters, buffersize);
1500 printf("Making Branch %s for Clusters\n",branchname);
1501 }
1502}
1503
1504//_____________________________________________________________________________
1505void AliTPC::ResetDigits()
1506{
1507 //
1508 // Reset number of digits and the digits array for this detector
1509 // reset clusters
1510 //
1511 fNdigits = 0;
1512 if (fDigits) fDigits->Clear();
1513 fNclusters = 0;
1514 if (fClusters) fClusters->Clear();
1515}
1516
1517//_____________________________________________________________________________
1518void AliTPC::SetSecAL(Int_t sec)
1519{
1520 //
1521 // Activate/deactivate selection for lower sectors
1522 //
1523 fSecAL = sec;
1524}
1525
1526//_____________________________________________________________________________
1527void AliTPC::SetSecAU(Int_t sec)
1528{
1529 //
1530 // Activate/deactivate selection for upper sectors
1531 //
1532 fSecAU = sec;
1533}
1534
1535//_____________________________________________________________________________
1536void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6)
1537{
1538 //
1539 // Select active lower sectors
1540 //
1541 fSecLows[0] = s1;
1542 fSecLows[1] = s2;
1543 fSecLows[2] = s3;
1544 fSecLows[3] = s4;
1545 fSecLows[4] = s5;
1546 fSecLows[5] = s6;
1547}
1548
1549//_____________________________________________________________________________
1550void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6,
1551 Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10,
1552 Int_t s11 , Int_t s12)
1553{
1554 //
1555 // Select active upper sectors
1556 //
1557 fSecUps[0] = s1;
1558 fSecUps[1] = s2;
1559 fSecUps[2] = s3;
1560 fSecUps[3] = s4;
1561 fSecUps[4] = s5;
1562 fSecUps[5] = s6;
1563 fSecUps[6] = s7;
1564 fSecUps[7] = s8;
1565 fSecUps[8] = s9;
1566 fSecUps[9] = s10;
1567 fSecUps[10] = s11;
1568 fSecUps[11] = s12;
1569}
1570
1571//_____________________________________________________________________________
1572void AliTPC::SetSens(Int_t sens)
1573{
1574 fSens = sens;
1575}
1576
1577//_____________________________________________________________________________
1578void AliTPC::Streamer(TBuffer &R__b)
1579{
1580 //
1581 // Stream an object of class AliTPC.
1582 //
1583 if (R__b.IsReading()) {
1584 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1585 AliDetector::Streamer(R__b);
1586 if (R__v < 2) return;
1587 R__b >> fNsectors;
1588 R__b >> fNclusters;
1589 R__b >> fNtracks;
1590 fClustersIndex = new Int_t[fNsectors+1];
1591 fDigitsIndex = new Int_t[fNsectors+1];
1592 } else {
1593 R__b.WriteVersion(AliTPC::IsA());
1594 AliDetector::Streamer(R__b);
1595 R__b << fNsectors;
1596 R__b << fNclusters;
1597 R__b << fNtracks;
1598 }
1599}
1600
1601ClassImp(AliTPCcluster)
1602
1603//_____________________________________________________________________________
1604AliTPCcluster::AliTPCcluster(Float_t *hits, Int_t *lab)
1605{
1606 //
1607 // Creates a simulated cluster for the TPC
1608 //
1609 fTracks[0] = lab[0];
1610 fTracks[1] = lab[1];
1611 fTracks[2] = lab[2];
1612 fSector = lab[3];
1613 fPadRow = lab[4];
1614 fY = hits[0];
1615 fZ = hits[1];
1616 fQ = hits[2];
1617 fSigmaY2 = hits[3];
1618 fSigmaZ2 = hits[4];
1619}
1620
1621//_____________________________________________________________________________
1622void AliTPCcluster::GetXYZ(Double_t& x, Double_t& y, Double_t& z) const
1623{
1624 //
1625 // Returns centroid for of a simulated TPC cluster
1626 //
1627 Double_t alpha,xrow;
1628 if (fSector<24) {
1629 alpha=(fSector<12) ? fSector*alpha_low : (fSector-12)*alpha_low;
1630 xrow=pad_row_low[fPadRow];
1631 } else {
1632 alpha=(fSector<48) ? (fSector-24)*alpha_up : (fSector-48)*alpha_up;
1633 xrow=pad_row_up[fPadRow];
1634 }
1635 x=xrow*cos(alpha) - fY*sin(alpha);
1636 y=xrow*sin(alpha) + fY*cos(alpha);
1637 z=fZ;
1638}
1639
1640ClassImp(AliTPCdigit)
1641
1642//_____________________________________________________________________________
1643AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits):
1644 AliDigit(tracks)
1645{
1646 //
1647 // Creates a TPC digit object
1648 //
1649 fSector = digits[0];
1650 fPadRow = digits[1];
1651 fPad = digits[2];
1652 fTime = digits[3];
1653 fSignal = digits[4];
1654}
1655
1656
1657ClassImp(AliTPChit)
1658
1659//_____________________________________________________________________________
1660AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
1661AliHit(shunt,track)
1662{
1663 //
1664 // Creates a TPC hit object
1665 //
1666 fSector = vol[0];
1667 fPadRow = vol[1];
1668 fX = hits[0];
1669 fY = hits[1];
1670 fZ = hits[2];
1671 fQ = hits[3];
1672}
1673
1674
1675ClassImp(AliTPCtrack)
1676
1677//_____________________________________________________________________________
1678AliTPCtrack::AliTPCtrack(Float_t *hits)
1679{
1680 //
1681 // Default creator for a TPC reconstructed track object
1682 //
1683 ref=hits[0]; // This is dummy code !
1684}
1685
1686AliTPCtrack::AliTPCtrack(const AliTPCcluster& c,const TVector& xx,
1687 const TMatrix& CC):
1688 x(xx),C(CC),clusters(MAX_CLUSTER)
1689{
1690 //
1691 // Standard creator for a TPC reconstructed track object
1692 //
1693 chi2=0.;
1694 int sec=c.fSector, row=c.fPadRow;
1695 if (sec<24) {
1696 fAlpha=(sec%12)*alpha_low; ref=pad_row_low[0] + row*pad_pitch_l;
1697 } else {
1698 fAlpha=(sec%24)*alpha_up; ref=pad_row_up[0] + row*pad_pitch_l;
1699 }
1700 clusters.AddLast((AliTPCcluster*)(&c));
1701}
1702
1703//_____________________________________________________________________________
1704AliTPCtrack::AliTPCtrack(const AliTPCtrack& t) : x(t.x), C(t.C),
1705 clusters(t.clusters.GetEntriesFast())
1706{
1707 //
1708 // Copy creator for a TPC reconstructed track
1709 //
1710 ref=t.ref;
1711 chi2=t.chi2;
1712 fAlpha=t.fAlpha;
1713 int n=t.clusters.GetEntriesFast();
1714 for (int i=0; i<n; i++) clusters.AddLast(t.clusters.UncheckedAt(i));
1715}
1716
1717//_____________________________________________________________________________
1718Double_t AliTPCtrack::GetY(Double_t xk) const
1719{
1720 //
1721 //
1722 //
1723 Double_t c2=x(2)*xk - x(3);
1724 if (TMath::Abs(c2) >= 0.999) {
1725 if (*this>10) cerr<<*this<<" AliTPCtrack warning: No y for this x !\n";
1726 return 0.;
1727 }
1728 Double_t c1=x(2)*ref - x(3);
1729 Double_t r1=sqrt(1.-c1*c1), r2=sqrt(1.-c2*c2);
1730 Double_t dx=xk-ref;
1731 return x(0) + dx*(c1+c2)/(r1+r2);
1732}
1733
1734//_____________________________________________________________________________
1735int AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm)
1736{
1737 //
1738 // Propagate a TPC reconstructed track
1739 //
1740 if (TMath::Abs(x(2)*xk - x(3)) >= 0.999) {
1741 if (*this>3) cerr<<*this<<" AliTPCtrack warning: Propagation failed !\n";
1742 return 0;
1743 }
1744
1745 Double_t x1=ref, x2=x1+0.5*(xk-x1), dx=x2-x1, y1=x(0), z1=x(1);
1746 Double_t c1=x(2)*x1 - x(3), r1=sqrt(1.- c1*c1);
1747 Double_t c2=x(2)*x2 - x(3), r2=sqrt(1.- c2*c2);
1748
1749 x(0) += dx*(c1+c2)/(r1+r2);
1750 x(1) += dx*(c1+c2)/(c1*r2 + c2*r1)*x(4);
1751
1752 TMatrix F(5,5); F.UnitMatrix();
1753 Double_t rr=r1+r2, cc=c1+c2, xx=x1+x2;
1754 F(0,2)= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr);
1755 F(0,3)=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr);
1756 Double_t cr=c1*r2+c2*r1;
1757 F(1,2)= dx*x(4)*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr);
1758 F(1,3)=-dx*x(4)*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr);
1759 F(1,4)= dx*cc/cr;
1760 TMatrix tmp(F,TMatrix::kMult,C);
1761 C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
1762
1763 ref=x2;
1764
1765 //Multiple scattering******************
1766 Double_t ey=x(2)*ref - x(3);
1767 Double_t ex=sqrt(1-ey*ey);
1768 Double_t ez=x(4);
1769 TMatrix Q(5,5); Q=0.;
1770 Q(2,2)=ez*ez+ey*ey; Q(2,3)=-ex*ey; Q(2,4)=-ex*ez;
1771 Q(3,2)=Q(2,3); Q(3,3)= ez*ez+ex*ex; Q(3,4)=-ey*ez;
1772 Q(4,2)=Q(2,4); Q(4,3)= Q(3,4); Q(4,4)=1.;
1773
1774 F=0;
1775 F(2,2)=-x(2)*ex; F(2,3)=-x(2)*ey;
1776 F(3,2)=-ex*(x(2)*ref-ey); F(3,3)=-(1.+ x(2)*ref*ey - ey*ey);
1777 F(4,2)=-ez*ex; F(4,3)=-ez*ey; F(4,4)=1.;
1778
1779 tmp.Mult(F,Q);
1780 Q.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
1781
1782 Double_t p2=GetPt()*GetPt()*(1.+x(4)*x(4));
1783 Double_t beta2=p2/(p2 + pm*pm);
1784 Double_t d=sqrt((x1-ref)*(x1-ref)+(y1-x(0))*(y1-x(0))+(z1-x(1))*(z1-x(1)));
1785 d*=2.;
1786 Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho;
1787 Q*=theta2;
1788 C+=Q;
1789
1790 //Energy losses************************
1791 Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d*rho;
1792 if (x1 < x2) dE=-dE;
1793 x(2)*=(1.- sqrt(p2+pm*pm)/p2*dE);
1794
1795 x1=ref; x2=xk; y1=x(0); z1=x(1);
1796 c1=x(2)*x1 - x(3); r1=sqrt(1.- c1*c1);
1797 c2=x(2)*x2 - x(3); r2=sqrt(1.- c2*c2);
1798
1799 x(0) += dx*(c1+c2)/(r1+r2);
1800 x(1) += dx*(c1+c2)/(c1*r2 + c2*r1)*x(4);
1801
1802 F.UnitMatrix();
1803 rr=r1+r2; cc=c1+c2; xx=x1+x2;
1804 F(0,2)= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr);
1805 F(0,3)=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr);
1806 cr=c1*r2+c2*r1;
1807 F(1,2)= dx*x(4)*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr);
1808 F(1,3)=-dx*x(4)*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr);
1809 F(1,4)= dx*cc/cr;
1810 tmp.Mult(F,C);
1811 C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
1812
1813 ref=x2;
1814
1815 return 1;
1816}
1817
1818//_____________________________________________________________________________
1819void AliTPCtrack::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm)
1820{
1821 //
1822 // Propagate a reconstructed track from the vertex
1823 //
1824 Double_t c=x(2)*ref - x(3);
1825 Double_t tgf=-x(3)/(x(2)*x(0) + sqrt(1-c*c));
1826 Double_t snf=tgf/sqrt(1.+ tgf*tgf);
1827 Double_t xv=(x(3)+snf)/x(2);
1828 PropagateTo(xv,x0,rho,pm);
1829}
1830
1831//_____________________________________________________________________________
1832void AliTPCtrack::Update(const AliTPCcluster *c, Double_t chisq)
1833{
1834 //
1835 // Update statistics for a reconstructed TPC track
1836 //
1837 TMatrix H(2,5); H.UnitMatrix();
1838 TMatrix Ht(TMatrix::kTransposed,H);
1839 TVector m(2); m(0)=c->fY; m(1)=c->fZ;
1840 TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2;
1841
1842 TMatrix tmp(H,TMatrix::kMult,C);
1843 TMatrix R(tmp,TMatrix::kMult,Ht); R+=V;
1844
1845 Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0);
1846 R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1);
1847 R(1,0)*=-1; R(0,1)=R(1,0);
1848 R*=1./det;
1849
1850 //R.Invert();
1851
1852 TMatrix K(C,TMatrix::kMult,Ht); K*=R;
1853
1854 TVector savex=x;
1855 x*=H; x-=m; x*=-1; x*=K; x+=savex;
1856 if (TMath::Abs(x(2)*ref-x(3)) >= 0.999) {
1857 if (*this>3) cerr<<*this<<" AliTPCtrack warning: Filtering failed !\n";
1858 x=savex;
1859 return;
1860 }
1861
1862 TMatrix saveC=C;
1863 C.Mult(K,tmp); C-=saveC; C*=-1;
1864
1865 clusters.AddLast((AliTPCcluster*)c);
1866 chi2 += chisq;
1867}
1868
1869//_____________________________________________________________________________
1870int AliTPCtrack::Rotate(Double_t alpha)
1871{
1872 //
1873 // Rotate a reconstructed TPC track
1874 //
1875 fAlpha += alpha;
1876
1877 Double_t x1=ref, y1=x(0);
1878 Double_t ca=cos(alpha), sa=sin(alpha);
1879 Double_t r1=x(2)*ref - x(3);
1880
1881 ref = x1*ca + y1*sa;
1882 x(0)=-x1*sa + y1*ca;
1883 x(3)=x(3)*ca + (x(2)*y1 + sqrt(1.- r1*r1))*sa;
1884
1885 Double_t r2=x(2)*ref - x(3);
1886 if (TMath::Abs(r2) >= 0.999) {
1887 if (*this>3) cerr<<*this<<" AliTPCtrack warning: Rotation failed !\n";
1888 return 0;
1889 }
1890
1891 Double_t y0=x(0) + sqrt(1.- r2*r2)/x(2);
1892 if ((x(0)-y0)*x(2) >= 0.) {
1893 if (*this>3) cerr<<*this<<" AliTPCtrack warning: Rotation failed !!!\n";
1894 return 0;
1895 }
1896
1897 TMatrix F(5,5); F.UnitMatrix();
1898 F(0,0)=ca;
1899 F(3,0)=x(2)*sa;
1900 F(3,2)=(y1 - r1*x1/sqrt(1.- r1*r1))*sa;
1901 F(3,3)= ca + sa*r1/sqrt(1.- r1*r1);
1902 TMatrix tmp(F,TMatrix::kMult,C);
1903 // Double_t dy2=C(0,0);
1904 C.Mult(tmp,TMatrix(TMatrix::kTransposed,F));
1905 // C(0,0)+=dy2*sa*sa*r1*r1/(1.- r1*r1);
1906 // C(1,1)+=dy2*sa*sa*x(4)*x(4)/(1.- r1*r1);
1907
1908 return 1;
1909}
1910
1911//_____________________________________________________________________________
1912void AliTPCtrack::UseClusters() const
1913{
1914 //
1915 //
1916 //
1917 int num_of_clusters=clusters.GetEntriesFast();
1918 for (int i=0; i<num_of_clusters; i++) {
1919 //if (i<=14) continue;
1920 AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(i);
1921 c->Use();
1922 }
1923}
1924
1925//_____________________________________________________________________________
1926Double_t AliTPCtrack::GetPredictedChi2(const AliTPCcluster *c) const
1927{
1928 //
1929 // Calculate chi2 for a reconstructed TPC track
1930 //
1931 TMatrix H(2,5); H.UnitMatrix();
1932 TVector m(2); m(0)=c->fY; m(1)=c->fZ;
1933 TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2;
1934 TVector res=x; res*=H; res-=m; //res*=-1;
1935 TMatrix tmp(H,TMatrix::kMult,C);
1936 TMatrix R(tmp,TMatrix::kMult,TMatrix(TMatrix::kTransposed,H)); R+=V;
1937
1938 Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0);
1939 if (TMath::Abs(det) < 1.e-10) {
1940 if (*this>3) cerr<<*this<<" AliTPCtrack warning: Singular matrix !\n";
1941 return 1e10;
1942 }
1943 R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1);
1944 R(1,0)*=-1; R(0,1)=R(1,0);
1945 R*=1./det;
1946
1947 //R.Invert();
1948
1949 TVector r=res;
1950 res*=R;
1951 return r*res;
1952}
1953
1954//_____________________________________________________________________________
1955int AliTPCtrack::GetLab() const
1956{
1957 //
1958 //
1959 //
1960 int lab = 0;
1961 struct {
1962 int lab;
1963 int max;
1964 } s[MAX_CLUSTER]={{0,0}};
1965
1966 int i;
1967 int num_of_clusters=clusters.GetEntriesFast();
1968 for (i=0; i<num_of_clusters; i++) {
1969 AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(i);
1970 lab=c->fTracks[0]; if (lab<0) lab=-lab;
1971 int j;
1972 for (j=0; j<MAX_CLUSTER; j++)
1973 if (s[j].lab==lab || s[j].max==0) break;
1974 s[j].lab=lab;
1975 s[j].max++;
1976 }
1977
1978 int max=0;
1979 for (i=0; i<num_of_clusters; i++)
1980 if (s[i].max>max) {max=s[i].max; lab=s[i].lab;}
1981 if (lab>0) lab--;
1982
1983 if (1.-float(max)/num_of_clusters > 0.10) return -lab;
1984
1985 if (num_of_clusters < 6) return lab;
1986
1987 max=0;
1988 for (i=1; i<=6; i++) {
1989 AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(num_of_clusters-i);
1990 if (lab != TMath::Abs(c->fTracks[0])-1
1991 && lab != TMath::Abs(c->fTracks[1])-1
1992 && lab != TMath::Abs(c->fTracks[2])-1
1993 ) max++;
1994 }
1995 if (max>3) return -lab;
1996
1997 return lab;
1998}
1999
2000//_____________________________________________________________________________
2001void AliTPCtrack::GetPxPyPz(Double_t& px, Double_t& py, Double_t& pz) const
2002{
2003 //
2004 // Get reconstructed TPC track momentum
2005 //
2006 Double_t pt=0.3*FIELD/TMath::Abs(x(2))/100; // GeV/c
2007 Double_t r=x(2)*ref-x(3);
2008 Double_t y0=x(0) + sqrt(1.- r*r)/x(2);
2009 px=-pt*(x(0)-y0)*x(2); //cos(phi);
2010 py=-pt*(x(3)-ref*x(2)); //sin(phi);
2011 pz=pt*x(4);
2012 Double_t tmp=px*TMath::Cos(fAlpha) - py*TMath::Sin(fAlpha);
2013 py=px*TMath::Sin(fAlpha) + py*TMath::Cos(fAlpha);
2014 px=tmp;
2015}
2016
2017//_____________________________________________________________________________
2018//
2019// Classes for internal tracking use
2020//
2021
2022//_____________________________________________________________________________
2023void AliTPCRow::InsertCluster(const AliTPCcluster* c)
2024{
2025 //
2026 // Insert a cluster in the list
2027 //
2028 if (num_of_clusters==MAX_CLUSTER_PER_ROW) {
2029 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
2030 }
2031 if (num_of_clusters==0) {clusters[num_of_clusters++]=c; return;}
2032 int i=Find(c->fY);
2033 memmove(clusters+i+1 ,clusters+i,(num_of_clusters-i)*sizeof(AliTPCcluster*));
2034 clusters[i]=c; num_of_clusters++;
2035}
2036
2037//_____________________________________________________________________________
2038int AliTPCRow::Find(Double_t y) const
2039{
2040 //
2041 //
2042 //
2043 if (y <= clusters[0]->fY) return 0;
2044 if (y > clusters[num_of_clusters-1]->fY) return num_of_clusters;
2045 int b=0, e=num_of_clusters-1, m=(b+e)/2;
2046 for (; b<e; m=(b+e)/2) {
2047 if (y > clusters[m]->fY) b=m+1;
2048 else e=m;
2049 }
2050 return m;
2051}