]>
Commit | Line | Data |
---|---|---|
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 | /* | |
1439f98e | 10 | <img src="picts/AliTPCClass.gif"> |
fe4da5cc | 11 | */ |
12 | //End_Html | |
13 | // // | |
8c555625 | 14 | // // |
fe4da5cc | 15 | /////////////////////////////////////////////////////////////////////////////// |
16 | ||
17 | #include <TMath.h> | |
18 | #include <TRandom.h> | |
19 | #include <TVector.h> | |
8c555625 | 20 | #include <TMatrix.h> |
fe4da5cc | 21 | #include <TGeometry.h> |
22 | #include <TNode.h> | |
23 | #include <TTUBS.h> | |
24 | #include <TObjectTable.h> | |
1578254f | 25 | #include "TParticle.h" |
fe4da5cc | 26 | #include "AliTPC.h" |
27 | #include "AliRun.h" | |
28 | #include <iostream.h> | |
29 | #include <fstream.h> | |
30 | #include "AliMC.h" | |
31 | ||
8c555625 | 32 | //MI change |
33 | #include "AliTPCParam.h" | |
34 | #include "AliTPCD.h" | |
35 | #include "AliTPCPRF2D.h" | |
36 | #include "AliTPCRF1D.h" | |
37 | ||
38 | ||
1283eee5 | 39 | |
fe4da5cc | 40 | ClassImp(AliTPC) |
41 | ||
42 | //_____________________________________________________________________________ | |
43 | AliTPC::AliTPC() | |
44 | { | |
45 | // | |
46 | // Default constructor | |
47 | // | |
48 | fIshunt = 0; | |
49 | fClusters = 0; | |
50 | fHits = 0; | |
51 | fDigits = 0; | |
52 | fTracks = 0; | |
53 | fNsectors = 0; | |
54 | fNtracks = 0; | |
55 | fNclusters= 0; | |
8c555625 | 56 | //MI changes |
57 | fDigParam= new AliTPCD(); | |
58 | fDigits = fDigParam->GetArray(); | |
fe4da5cc | 59 | } |
60 | ||
61 | //_____________________________________________________________________________ | |
62 | AliTPC::AliTPC(const char *name, const char *title) | |
63 | : AliDetector(name,title) | |
64 | { | |
65 | // | |
66 | // Standard constructor | |
67 | // | |
68 | ||
69 | // | |
70 | // Initialise arrays of hits and digits | |
71 | fHits = new TClonesArray("AliTPChit", 176); | |
8c555625 | 72 | // fDigits = new TClonesArray("AliTPCdigit",10000); |
73 | //MI change | |
74 | fDigParam= new AliTPCD; | |
75 | fDigits = fDigParam->GetArray(); | |
fe4da5cc | 76 | // |
77 | // Initialise counters | |
78 | fClusters = 0; | |
79 | fTracks = 0; | |
80 | fNsectors = 72; | |
81 | fNtracks = 0; | |
82 | fNclusters= 0; | |
83 | fDigitsIndex = new Int_t[fNsectors+1]; | |
84 | fClustersIndex = new Int_t[fNsectors+1]; | |
85 | // | |
86 | fIshunt = 0; | |
87 | // | |
88 | // Initialise color attributes | |
89 | SetMarkerColor(kYellow); | |
90 | } | |
91 | ||
92 | //_____________________________________________________________________________ | |
93 | AliTPC::~AliTPC() | |
94 | { | |
95 | // | |
96 | // TPC destructor | |
97 | // | |
98 | fIshunt = 0; | |
99 | delete fHits; | |
100 | delete fDigits; | |
101 | delete fClusters; | |
102 | delete fTracks; | |
8c555625 | 103 | delete fDigParam; |
fe4da5cc | 104 | if (fDigitsIndex) delete [] fDigitsIndex; |
105 | if (fClustersIndex) delete [] fClustersIndex; | |
106 | } | |
107 | ||
108 | //_____________________________________________________________________________ | |
109 | void AliTPC::AddCluster(Float_t *hits, Int_t *tracks) | |
110 | { | |
111 | // | |
112 | // Add a simulated cluster to the list | |
113 | // | |
114 | if(!fClusters) fClusters=new TClonesArray("AliTPCcluster",10000); | |
115 | TClonesArray &lclusters = *fClusters; | |
116 | new(lclusters[fNclusters++]) AliTPCcluster(hits,tracks); | |
117 | } | |
118 | ||
119 | //_____________________________________________________________________________ | |
120 | void AliTPC::AddCluster(const AliTPCcluster &c) | |
121 | { | |
122 | // | |
123 | // Add a simulated cluster copy to the list | |
124 | // | |
125 | if(!fClusters) fClusters=new TClonesArray("AliTPCcluster",10000); | |
126 | TClonesArray &lclusters = *fClusters; | |
127 | new(lclusters[fNclusters++]) AliTPCcluster(c); | |
128 | } | |
129 | ||
130 | //_____________________________________________________________________________ | |
131 | void AliTPC::AddDigit(Int_t *tracks, Int_t *digits) | |
132 | { | |
133 | // | |
134 | // Add a TPC digit to the list | |
135 | // | |
8c555625 | 136 | // TClonesArray &ldigits = *fDigits; |
137 | //MI change | |
138 | TClonesArray &ldigits = *fDigParam->GetArray(); | |
fe4da5cc | 139 | new(ldigits[fNdigits++]) AliTPCdigit(tracks,digits); |
140 | } | |
141 | ||
142 | //_____________________________________________________________________________ | |
143 | void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits) | |
144 | { | |
145 | // | |
146 | // Add a hit to the list | |
147 | // | |
148 | TClonesArray &lhits = *fHits; | |
149 | new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits); | |
150 | } | |
151 | ||
152 | //_____________________________________________________________________________ | |
153 | void AliTPC::AddTrack(Float_t *hits) | |
154 | { | |
155 | // | |
156 | // Add a track to the list of tracks | |
157 | // | |
158 | TClonesArray <racks = *fTracks; | |
159 | new(ltracks[fNtracks++]) AliTPCtrack(hits); | |
160 | } | |
161 | ||
162 | //_____________________________________________________________________________ | |
163 | void AliTPC::AddTrack(const AliTPCtrack& t) | |
164 | { | |
165 | // | |
166 | // Add a track copy to the list of tracks | |
167 | // | |
168 | if(!fTracks) fTracks=new TClonesArray("AliTPCtrack",10000); | |
169 | TClonesArray <racks = *fTracks; | |
170 | new(ltracks[fNtracks++]) AliTPCtrack(t); | |
171 | } | |
172 | ||
173 | //_____________________________________________________________________________ | |
174 | void AliTPC::BuildGeometry() | |
175 | { | |
176 | // | |
177 | // Build TPC ROOT TNode geometry for the event display | |
178 | // | |
179 | TNode *Node, *Top; | |
180 | TTUBS *tubs; | |
181 | Int_t i; | |
182 | const int kColorTPC=19; | |
1283eee5 | 183 | char name[5], title[25]; |
fe4da5cc | 184 | const Double_t kDegrad=TMath::Pi()/180; |
1283eee5 | 185 | const Double_t kRaddeg=180./TMath::Pi(); |
186 | ||
187 | AliTPCParam * fTPCParam = &(fDigParam->GetParam()); | |
188 | ||
189 | Float_t InnerOpenAngle = fTPCParam->GetInnerAngle(); | |
190 | Float_t OuterOpenAngle = fTPCParam->GetOuterAngle(); | |
191 | ||
192 | Float_t InnerAngleShift = fTPCParam->GetInnerAngleShift(); | |
193 | Float_t OuterAngleShift = fTPCParam->GetOuterAngleShift(); | |
194 | ||
195 | Int_t nLo = fTPCParam->GetNInnerSector()/2; | |
196 | Int_t nHi = fTPCParam->GetNOuterSector()/2; | |
197 | ||
198 | const Double_t loAng = (Double_t)TMath::Nint(InnerOpenAngle*kRaddeg); | |
199 | const Double_t hiAng = (Double_t)TMath::Nint(OuterOpenAngle*kRaddeg); | |
200 | const Double_t loAngSh = (Double_t)TMath::Nint(InnerAngleShift*kRaddeg); | |
201 | const Double_t hiAngSh = (Double_t)TMath::Nint(OuterAngleShift*kRaddeg); | |
202 | ||
203 | ||
fe4da5cc | 204 | const Double_t loCorr = 1/TMath::Cos(0.5*loAng*kDegrad); |
205 | const Double_t hiCorr = 1/TMath::Cos(0.5*hiAng*kDegrad); | |
1283eee5 | 206 | |
207 | Double_t rl,ru; | |
208 | ||
209 | ||
fe4da5cc | 210 | // |
211 | // Get ALICE top node | |
fe4da5cc | 212 | // |
1283eee5 | 213 | |
214 | Top=gAlice->GetGeometry()->GetNode("alice"); | |
215 | ||
216 | // inner sectors | |
217 | ||
218 | rl = fTPCParam->GetInSecLowEdge(); | |
219 | ru = fTPCParam->GetInSecUpEdge(); | |
220 | ||
221 | ||
fe4da5cc | 222 | for(i=0;i<nLo;i++) { |
223 | sprintf(name,"LS%2.2d",i); | |
1283eee5 | 224 | name[4]='\0'; |
225 | sprintf(title,"TPC low sector %3d",i); | |
226 | title[24]='\0'; | |
227 | ||
228 | tubs = new TTUBS(name,title,"void",rl*loCorr,ru*loCorr,250., | |
229 | loAng*(i-0.5)+loAngSh,loAng*(i+0.5)+loAngSh); | |
fe4da5cc | 230 | tubs->SetNumberOfDivisions(1); |
231 | Top->cd(); | |
232 | Node = new TNode(name,title,name,0,0,0,""); | |
233 | Node->SetLineColor(kColorTPC); | |
234 | fNodes->Add(Node); | |
235 | } | |
1283eee5 | 236 | |
fe4da5cc | 237 | // Outer sectors |
1283eee5 | 238 | |
239 | rl = fTPCParam->GetOuSecLowEdge(); | |
240 | ru = fTPCParam->GetOuSecUpEdge(); | |
241 | ||
fe4da5cc | 242 | for(i=0;i<nHi;i++) { |
243 | sprintf(name,"US%2.2d",i); | |
1283eee5 | 244 | name[4]='\0'; |
fe4da5cc | 245 | sprintf(title,"TPC upper sector %d",i); |
1283eee5 | 246 | title[24]='\0'; |
247 | tubs = new TTUBS(name,title,"void",rl*hiCorr,ru*hiCorr,250, | |
248 | hiAng*(i-0.5)+hiAngSh,hiAng*(i+0.5)+hiAngSh); | |
fe4da5cc | 249 | tubs->SetNumberOfDivisions(1); |
250 | Top->cd(); | |
251 | Node = new TNode(name,title,name,0,0,0,""); | |
252 | Node->SetLineColor(kColorTPC); | |
253 | fNodes->Add(Node); | |
254 | } | |
1283eee5 | 255 | } |
256 | ||
257 | ||
258 | ||
fe4da5cc | 259 | //_____________________________________________________________________________ |
260 | Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t ) | |
261 | { | |
262 | // | |
263 | // Calculate distance from TPC to mouse on the display | |
264 | // Dummy procedure | |
265 | // | |
266 | return 9999; | |
267 | } | |
268 | ||
269 | //_____________________________________________________________________________ | |
8c555625 | 270 | //const int MAX_CLUSTER=nrow_low+nrow_up; |
271 | const int MAX_CLUSTER=200; | |
fe4da5cc | 272 | const int S_MAXSEC=24; |
273 | const int L_MAXSEC=48; | |
274 | const int ROWS_TO_SKIP=21; | |
8c555625 | 275 | const Float_t MAX_CHI2=12.; |
276 | ||
fe4da5cc | 277 | |
278 | //_____________________________________________________________________________ | |
279 | static Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt) | |
280 | { | |
281 | // | |
282 | // Calculate spread in Y | |
283 | // | |
284 | pt=TMath::Abs(pt)*1000.; | |
285 | Double_t x=r/pt; | |
286 | tgl=TMath::Abs(tgl); | |
287 | Double_t s=a_rphi - b_rphi*r*tgl + c_rphi*x*x + d_rphi*x; | |
288 | if (s<0.4e-3) s=0.4e-3; | |
289 | return s; | |
290 | } | |
291 | ||
292 | //_____________________________________________________________________________ | |
293 | static Double_t SigmaZ2(Double_t r, Double_t tgl) | |
294 | { | |
295 | // | |
296 | // Calculate spread in Z | |
297 | // | |
298 | tgl=TMath::Abs(tgl); | |
299 | Double_t s=a_z - b_z*r*tgl + c_z*tgl*tgl; | |
300 | if (s<0.4e-3) s=0.4e-3; | |
301 | return s; | |
302 | } | |
303 | ||
304 | //_____________________________________________________________________________ | |
305 | inline Double_t f1(Double_t x1,Double_t y1, //C | |
306 | Double_t x2,Double_t y2, | |
307 | Double_t x3,Double_t y3) | |
308 | { | |
309 | // | |
310 | // Function f1 | |
311 | // | |
312 | Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1); | |
313 | Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)- | |
314 | (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2)); | |
315 | Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)- | |
316 | (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1)); | |
8c555625 | 317 | |
fe4da5cc | 318 | Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b); |
8c555625 | 319 | |
fe4da5cc | 320 | return -xr*yr/sqrt(xr*xr+yr*yr); |
321 | } | |
322 | ||
323 | ||
324 | //_____________________________________________________________________________ | |
325 | inline Double_t f2(Double_t x1,Double_t y1, //eta=C*x0 | |
326 | Double_t x2,Double_t y2, | |
327 | Double_t x3,Double_t y3) | |
328 | { | |
329 | // | |
330 | // Function f2 | |
331 | // | |
332 | Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1); | |
333 | Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)- | |
334 | (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2)); | |
335 | Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)- | |
336 | (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1)); | |
8c555625 | 337 | |
fe4da5cc | 338 | Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b); |
339 | ||
340 | return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr); | |
341 | } | |
342 | ||
343 | //_____________________________________________________________________________ | |
344 | inline Double_t f3(Double_t x1,Double_t y1, //tgl | |
345 | Double_t x2,Double_t y2, | |
346 | Double_t z1,Double_t z2) | |
347 | { | |
348 | // | |
349 | // Function f3 | |
350 | // | |
351 | return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); | |
352 | } | |
353 | ||
354 | //_____________________________________________________________________________ | |
355 | static int FindProlongation(AliTPCtrack& t, const AliTPCSector *sec, | |
356 | int s, int ri, int rf=0) | |
357 | { | |
358 | // | |
359 | // Propagate track | |
360 | // | |
361 | int try_again=ROWS_TO_SKIP; | |
362 | Double_t alpha=sec->GetAlpha(); | |
363 | int ns=int(2*TMath::Pi()/alpha)+1; | |
8c555625 | 364 | |
fe4da5cc | 365 | for (int nr=ri; nr>=rf; nr--) { |
366 | Double_t x=sec[s].GetX(nr), ymax=sec[s].GetMaxY(nr); | |
8c555625 | 367 | if (!t.PropagateTo(x)) return -1; |
368 | ||
fe4da5cc | 369 | const AliTPCcluster *cl=0; |
370 | Double_t max_chi2=MAX_CHI2; | |
371 | const AliTPCRow& row=sec[s][nr]; | |
372 | Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),t.GetPt()); | |
373 | Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl()); | |
8c555625 | 374 | Double_t road=3.*sqrt(t.GetSigmaY2() + 4*sy2), y=t.GetY(), z=t.GetZ(); |
375 | ||
fe4da5cc | 376 | if (road>30) { |
377 | if (t>3) cerr<<t<<" AliTPCtrack warning: Too broad road !\n"; | |
8c555625 | 378 | return -1; |
fe4da5cc | 379 | } |
8c555625 | 380 | |
381 | if (row) { | |
fe4da5cc | 382 | for (int i=row.Find(y-road); i<row; i++) { |
383 | AliTPCcluster* c=(AliTPCcluster*)(row[i]); | |
384 | if (c->fY > y+road) break; | |
385 | if (c->IsUsed()) continue; | |
8c555625 | 386 | if ((c->fZ - z)*(c->fZ - z) > 9.*(t.GetSigmaZ2() + 4*sz2)) continue; |
fe4da5cc | 387 | Double_t chi2=t.GetPredictedChi2(c); |
388 | if (chi2 > max_chi2) continue; | |
389 | max_chi2=chi2; | |
390 | cl=c; | |
391 | } | |
8c555625 | 392 | } |
fe4da5cc | 393 | if (cl) { |
394 | t.Update(cl,max_chi2); | |
395 | try_again=ROWS_TO_SKIP; | |
396 | } else { | |
8c555625 | 397 | if (try_again==0) break; |
398 | if (y > ymax) { | |
399 | s = (s+1) % ns; | |
400 | if (!t.Rotate(alpha)) return -1; | |
401 | } else if (y <-ymax) { | |
402 | s = (s-1+ns) % ns; | |
403 | if (!t.Rotate(-alpha)) return -1; | |
fe4da5cc | 404 | } |
8c555625 | 405 | try_again--; |
fe4da5cc | 406 | } |
407 | } | |
8c555625 | 408 | |
fe4da5cc | 409 | return s; |
410 | } | |
411 | ||
8c555625 | 412 | |
fe4da5cc | 413 | //_____________________________________________________________________________ |
8c555625 | 414 | static void MakeSeeds(TObjArray& seeds,const AliTPCSector* sec,int i1,int i2, |
415 | const AliTPCParam *p) | |
fe4da5cc | 416 | { |
417 | // | |
418 | // Find seed for tracking | |
419 | // | |
420 | TMatrix C(5,5); TVector x(5); | |
421 | int max_sec=L_MAXSEC/2; | |
422 | for (int ns=0; ns<max_sec; ns++) { | |
423 | int nl=sec[(ns-1+max_sec)%max_sec][i2]; | |
424 | int nm=sec[ns][i2]; | |
425 | int nu=sec[(ns+1)%max_sec][i2]; | |
426 | Double_t alpha=sec[ns].GetAlpha(); | |
427 | const AliTPCRow& r1=sec[ns][i1]; | |
428 | for (int is=0; is < r1; is++) { | |
429 | Double_t x1=sec[ns].GetX(i1), y1=r1[is]->fY, z1=r1[is]->fZ; | |
430 | for (int js=0; js < nl+nm+nu; js++) { | |
431 | const AliTPCcluster *cl; | |
432 | Double_t cs,sn; | |
8c555625 | 433 | int ks; |
fe4da5cc | 434 | |
435 | if (js<nl) { | |
8c555625 | 436 | ks=(ns-1+max_sec)%max_sec; |
fe4da5cc | 437 | const AliTPCRow& r2=sec[(ns-1+max_sec)%max_sec][i2]; |
438 | cl=r2[js]; | |
439 | cs=cos(alpha); sn=sin(alpha); | |
440 | } else | |
441 | if (js<nl+nm) { | |
8c555625 | 442 | ks=ns; |
fe4da5cc | 443 | const AliTPCRow& r2=sec[ns][i2]; |
444 | cl=r2[js-nl]; | |
445 | cs=1; sn=0.; | |
446 | } else { | |
8c555625 | 447 | ks=(ns+1)%max_sec; |
fe4da5cc | 448 | const AliTPCRow& r2=sec[(ns+1)%max_sec][i2]; |
449 | cl=r2[js-nl-nm]; | |
450 | cs=cos(alpha); sn=-sin(alpha); | |
451 | } | |
452 | Double_t x2=sec[ns].GetX(i2), y2=cl->fY, z2=cl->fZ; | |
453 | //if (z1*z2 < 0) continue; | |
454 | //if (TMath::Abs(z1) < TMath::Abs(z2)) continue; | |
455 | ||
456 | Double_t tmp= x2*cs+y2*sn; | |
457 | y2 =-x2*sn+y2*cs; | |
458 | x2=tmp; | |
459 | ||
460 | x(0)=y1; | |
461 | x(1)=z1; | |
462 | x(2)=f1(x1,y1,x2,y2,0.,0.); | |
463 | x(3)=f2(x1,y1,x2,y2,0.,0.); | |
464 | x(4)=f3(x1,y1,x2,y2,z1,z2); | |
465 | ||
466 | if (TMath::Abs(x(2)*x1-x(3)) >= 0.999) continue; | |
467 | ||
468 | if (TMath::Abs(x(4)) > 1.2) continue; | |
8c555625 | 469 | |
fe4da5cc | 470 | Double_t a=asin(x(3)); |
fe4da5cc | 471 | Double_t zv=z1 - x(4)/x(2)*(a+asin(x(2)*x1-x(3))); |
472 | if (TMath::Abs(zv)>33.) continue; | |
8c555625 | 473 | |
fe4da5cc | 474 | TMatrix X(6,6); X=0.; |
475 | X(0,0)=r1[is]->fSigmaY2; X(1,1)=r1[is]->fSigmaZ2; | |
476 | X(2,2)=cl->fSigmaY2; X(3,3)=cl->fSigmaZ2; | |
477 | X(4,4)=3./12.; X(5,5)=3./12.; | |
478 | TMatrix F(5,6); F.UnitMatrix(); | |
479 | Double_t sy=sqrt(X(0,0)), sz=sqrt(X(1,1)); | |
480 | F(2,0)=(f1(x1,y1+sy,x2,y2,0.,0.)-x(2))/sy; | |
481 | F(2,2)=(f1(x1,y1,x2,y2+sy,0.,0.)-x(2))/sy; | |
482 | F(2,4)=(f1(x1,y1,x2,y2,0.,0.+sy)-x(2))/sy; | |
483 | F(3,0)=(f2(x1,y1+sy,x2,y2,0.,0.)-x(3))/sy; | |
484 | F(3,2)=(f2(x1,y1,x2,y2+sy,0.,0.)-x(3))/sy; | |
485 | F(3,4)=(f2(x1,y1,x2,y2,0.,0.+sy)-x(3))/sy; | |
486 | F(4,0)=(f3(x1,y1+sy,x2,y2,z1,z2)-x(4))/sy; | |
487 | F(4,1)=(f3(x1,y1,x2,y2,z1+sz,z2)-x(4))/sz; | |
488 | F(4,2)=(f3(x1,y1,x2,y2+sy,z1,z2)-x(4))/sy; | |
489 | F(4,3)=(f3(x1,y1,x2,y2,z1,z2+sz)-x(4))/sz; | |
490 | F(4,4)=0; | |
491 | F(3,3)=0; | |
492 | ||
493 | TMatrix t(F,TMatrix::kMult,X); | |
494 | C.Mult(t,TMatrix(TMatrix::kTransposed,F)); | |
8c555625 | 495 | |
496 | TrackSeed *track=new TrackSeed(*(r1[is]),x,C,p); | |
497 | int rc=FindProlongation(*track,sec,ns,i1-1,i2); | |
498 | if (rc<0 || *track<(i1-i2)/2) delete track; | |
499 | else seeds.AddLast(track); | |
fe4da5cc | 500 | } |
501 | } | |
502 | } | |
503 | } | |
504 | ||
505 | //_____________________________________________________________________________ | |
506 | void AliTPC::Clusters2Tracks() | |
507 | { | |
508 | // | |
509 | // TPC Track finder from clusters. | |
510 | // | |
511 | if (!fClusters) return; | |
8c555625 | 512 | |
513 | AliTPCParam *p=&fDigParam->GetParam(); | |
514 | Int_t nrow_low=p->GetNRowLow(); | |
515 | Int_t nrow_up=p->GetNRowUp(); | |
516 | ||
fe4da5cc | 517 | AliTPCSSector ssec[S_MAXSEC/2]; |
8c555625 | 518 | for (int i=0; i<S_MAXSEC/2; i++) ssec[i].SetUp(p); |
519 | ||
fe4da5cc | 520 | AliTPCLSector lsec[L_MAXSEC/2]; |
8c555625 | 521 | for (int j=0; j<L_MAXSEC/2; j++) lsec[j].SetUp(p); |
522 | ||
fe4da5cc | 523 | int ncl=fClusters->GetEntriesFast(); |
524 | while (ncl--) { | |
525 | AliTPCcluster *c=(AliTPCcluster*)fClusters->UncheckedAt(ncl); | |
8c555625 | 526 | |
527 | int sec=int(c->fSector)-1, row=int(c->fPadRow)-1; | |
fe4da5cc | 528 | |
529 | if (sec<24) { | |
530 | if (row<0 || row>nrow_low) {cerr<<"low !!!"<<row<<endl; continue;} | |
531 | ssec[sec%12][row].InsertCluster(c); | |
532 | } else { | |
533 | if (row<0 || row>nrow_up ) {cerr<<"up !!!"<<row<<endl; continue;} | |
534 | sec -= 24; | |
535 | lsec[sec%24][row].InsertCluster(c); | |
536 | } | |
537 | } | |
538 | ||
539 | ||
540 | TObjArray seeds(20000); | |
8c555625 | 541 | MakeSeeds(seeds,lsec,nrow_up-1,nrow_up-1-8,p); |
542 | MakeSeeds(seeds,lsec,nrow_up-1-4,nrow_up-1-4-8,p); | |
543 | ||
fe4da5cc | 544 | seeds.Sort(); |
545 | ||
546 | int found=0; | |
547 | int nseed=seeds.GetEntriesFast(); | |
548 | ||
549 | for (int s=0; s<nseed; s++) { | |
550 | AliTPCtrack& t=*((AliTPCtrack*)seeds.UncheckedAt(s)); | |
551 | Double_t alpha=t.GetAlpha(); | |
552 | if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); | |
553 | if (alpha < 0. ) alpha += 2.*TMath::Pi(); | |
554 | int ns=int(alpha/lsec->GetAlpha() + 0.5); | |
555 | ||
556 | Double_t x=t.GetX(); | |
557 | int nr; | |
8c555625 | 558 | if (x<p->GetPadRowRadiiUp(nrow_up-1-4-7)) nr=nrow_up-1-4-8; |
559 | else if (x<p->GetPadRowRadiiUp(nrow_up-1-7)) nr=nrow_up-1-8; | |
fe4da5cc | 560 | else {cerr<<x<<" =x !!!\n"; continue;} |
8c555625 | 561 | |
fe4da5cc | 562 | int ls=FindProlongation(t,lsec,ns,nr-1); |
8c555625 | 563 | if (ls<0) continue; |
fe4da5cc | 564 | x=t.GetX(); alpha=lsec[ls].GetAlpha(); // |
565 | Double_t phi=ls*alpha + atan(t.GetY()/x); // Find S-sector | |
566 | int ss=int(0.5*(phi/alpha+1)); // | |
567 | alpha *= 2*(ss-0.5*ls); // and rotation angle | |
568 | if (!t.Rotate(alpha)) continue; // | |
569 | ss %= (S_MAXSEC/2); // | |
570 | ||
8c555625 | 571 | if (FindProlongation(t,ssec,ss,nrow_low-1)<0) continue; |
fe4da5cc | 572 | if (t < 30) continue; |
573 | ||
574 | AddTrack(t); | |
575 | t.UseClusters(); | |
576 | cerr<<found++<<'\r'; | |
577 | } | |
578 | } | |
579 | ||
580 | //_____________________________________________________________________________ | |
581 | void AliTPC::CreateMaterials() | |
582 | { | |
8c555625 | 583 | //----------------------------------------------- |
fe4da5cc | 584 | // Create Materials for for TPC |
8c555625 | 585 | //----------------------------------------------- |
586 | ||
587 | //----------------------------------------------------------------- | |
588 | // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl | |
589 | //----------------------------------------------------------------- | |
1283eee5 | 590 | |
fe4da5cc | 591 | Int_t ISXFLD=gAlice->Field()->Integ(); |
592 | Float_t SXMGMX=gAlice->Field()->Max(); | |
1283eee5 | 593 | |
594 | Float_t amat[5]; // atomic numbers | |
595 | Float_t zmat[5]; // z | |
596 | Float_t wmat[5]; // proportions | |
597 | ||
598 | Float_t density; | |
599 | ||
600 | // ********************* Gases ******************* | |
601 | ||
602 | //-------------------------------------------------------------- | |
603 | // pure gases | |
604 | //-------------------------------------------------------------- | |
605 | ||
606 | // Ne | |
607 | ||
608 | ||
609 | Float_t a_ne = 20.18; | |
610 | Float_t z_ne = 10.; | |
fe4da5cc | 611 | |
1283eee5 | 612 | density = 0.0009; |
613 | ||
614 | AliMaterial(20,"Ne",a_ne,z_ne,density,999.,999.); | |
615 | ||
616 | // Ar | |
617 | ||
618 | Float_t a_ar = 39.948; | |
619 | Float_t z_ar = 18.; | |
620 | ||
621 | density = 0.001782; | |
622 | ||
623 | AliMaterial(21,"Ar",a_ar,z_ar,density,999.,999.); | |
624 | ||
625 | Float_t a_pure[2]; | |
fe4da5cc | 626 | |
1283eee5 | 627 | a_pure[0] = a_ne; |
628 | a_pure[1] = a_ar; | |
fe4da5cc | 629 | |
1283eee5 | 630 | |
631 | //-------------------------------------------------------------- | |
632 | // gases - compounds | |
633 | //-------------------------------------------------------------- | |
634 | ||
635 | Float_t amol[3]; | |
636 | ||
637 | // CO2 | |
638 | ||
639 | amat[0]=12.011; | |
640 | amat[1]=15.9994; | |
641 | ||
642 | zmat[0]=6.; | |
643 | zmat[1]=8.; | |
644 | ||
645 | wmat[0]=1.; | |
646 | wmat[1]=2.; | |
647 | ||
648 | density=0.001977; | |
649 | ||
650 | amol[0] = amat[0]*wmat[0]+amat[1]*wmat[1]; | |
651 | ||
652 | AliMixture(10,"CO2",amat,zmat,density,-2,wmat); | |
653 | ||
654 | // CF4 | |
655 | ||
656 | amat[0]=12.011; | |
657 | amat[1]=18.998; | |
658 | ||
659 | zmat[0]=6.; | |
660 | zmat[1]=9.; | |
661 | ||
662 | wmat[0]=1.; | |
663 | wmat[1]=4.; | |
664 | ||
665 | density=0.003034; | |
666 | ||
667 | amol[1] = amat[0]*wmat[0]+amat[1]*wmat[1]; | |
668 | ||
669 | AliMixture(11,"CF4",amat,zmat,density,-2,wmat); | |
670 | ||
671 | // CH4 | |
672 | ||
673 | amat[0]=12.011; | |
674 | amat[1]=1.; | |
675 | ||
676 | zmat[0]=6.; | |
677 | zmat[1]=1.; | |
678 | ||
679 | wmat[0]=1.; | |
680 | wmat[1]=4.; | |
681 | ||
682 | density=0.000717; | |
683 | ||
684 | amol[2] = amat[0]*wmat[0]+amat[1]*wmat[1]; | |
685 | ||
686 | AliMixture(12,"CH4",amat,zmat,density,-2,wmat); | |
687 | ||
688 | //---------------------------------------------------------------- | |
689 | // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds | |
690 | //---------------------------------------------------------------- | |
691 | ||
fe4da5cc | 692 | char namate[21]; |
1283eee5 | 693 | |
694 | density = 0.; | |
695 | Float_t am=0; | |
696 | Int_t nc; | |
697 | ||
698 | Float_t a,z,rho,absl,X0,buf[1]; | |
699 | Int_t nbuf; | |
700 | ||
701 | for(nc = 0;nc<fNoComp;nc++) | |
702 | { | |
703 | ||
704 | // retrive material constants | |
705 | ||
706 | gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,X0,absl,buf,nbuf); | |
707 | ||
708 | amat[nc] = a; | |
709 | zmat[nc] = z; | |
710 | ||
711 | Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10; | |
712 | ||
713 | am += fMixtProp[nc]*((fMixtComp[nc]>=20) ? a_pure[nnc] : amol[nnc]); | |
714 | density += fMixtProp[nc]*rho; // density of the mixture | |
715 | ||
716 | } | |
717 | ||
718 | // mixture proportions by weight! | |
719 | ||
720 | for(nc = 0;nc<fNoComp;nc++) | |
721 | { | |
722 | ||
723 | Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10; | |
724 | ||
725 | wmat[nc] = fMixtProp[nc]*((fMixtComp[nc]>=20) ? a_pure[nnc] : amol[nnc])/am; | |
726 | ||
727 | } | |
fe4da5cc | 728 | |
1283eee5 | 729 | AliMixture(31,"Drift gas 1",amat,zmat,density,fNoComp,wmat); |
730 | AliMixture(32,"Drift gas 2",amat,zmat,density,fNoComp,wmat); | |
731 | AliMixture(33,"Drift gas 3",amat,zmat,density,fNoComp,wmat); | |
732 | ||
733 | AliMedium(2, "Drift gas 1", 31, 0, ISXFLD, SXMGMX, 10., 999.,.1, .001, .001); | |
734 | AliMedium(3, "Drift gas 2", 32, 0, ISXFLD, SXMGMX, 10., 999.,.1, .001, .001); | |
735 | AliMedium(4, "Drift gas 3", 33, 1, ISXFLD, SXMGMX, 10., 999.,.1, .001, .001); | |
736 | ||
737 | // Air | |
738 | ||
739 | AliMaterial(24, "Air", 14.61, 7.3, .001205, 30420., 67500.); | |
740 | ||
741 | AliMedium(24, "Air", 24, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1); | |
742 | ||
743 | //---------------------------------------------------------------------- | |
744 | // solid materials | |
745 | //---------------------------------------------------------------------- | |
746 | ||
747 | // Al | |
748 | ||
749 | AliMaterial(30, "Al", 26.98, 13., 2.7, 8.9, 37.2); | |
750 | ||
751 | AliMedium(0, "Al",30, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1); | |
752 | ||
753 | // Si | |
754 | ||
755 | AliMaterial(31, "Si", 28.086, 14.,2.33, 9.36, 999.); | |
756 | ||
757 | AliMedium(7, "Al",31, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1); | |
fe4da5cc | 758 | |
1283eee5 | 759 | |
760 | // Mylar C5H4O2 | |
761 | ||
762 | amat[0]=12.011; | |
763 | amat[1]=1.; | |
764 | amat[2]=15.9994; | |
765 | ||
766 | zmat[0]=6.; | |
767 | zmat[1]=1.; | |
768 | zmat[2]=8.; | |
769 | ||
770 | wmat[0]=5.; | |
771 | wmat[1]=4.; | |
772 | wmat[2]=2.; | |
773 | ||
774 | density = 1.39; | |
fe4da5cc | 775 | |
1283eee5 | 776 | AliMixture(32, "Mylar",amat,zmat,density,-3,wmat); |
777 | ||
778 | AliMedium(5, "Mylar",32, 0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01); | |
779 | ||
780 | ||
781 | ||
782 | ||
783 | // Carbon (normal) | |
784 | ||
785 | AliMaterial(33,"C normal",12.011,6.,2.265,18.8,999.); | |
786 | ||
787 | AliMedium(6,"C normal",33,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01); | |
788 | ||
789 | // G10 for inner and outr field cage | |
790 | // G10 is 60% SiO2 + 40% epoxy, right now I use A and Z for SiO2 | |
791 | ||
792 | Float_t rhoFactor; | |
793 | ||
794 | amat[0]=28.086; | |
795 | amat[1]=15.9994; | |
796 | ||
797 | zmat[0]=14.; | |
798 | zmat[1]=8.; | |
799 | ||
800 | wmat[0]=1.; | |
801 | wmat[1]=2.; | |
802 | ||
803 | density = 1.7; | |
fe4da5cc | 804 | |
1283eee5 | 805 | |
806 | AliMixture(34,"G10 aux.",amat,zmat,density,-2,wmat); | |
807 | ||
808 | ||
809 | gMC->Gfmate((*fIdmate)[34],namate,a,z,rho,X0,absl,buf,nbuf); | |
810 | ||
811 | Float_t thickX0 = 0.0052; // field cage in X0 units | |
fe4da5cc | 812 | |
1283eee5 | 813 | Float_t thick = 2.; // in cm |
814 | ||
815 | X0=19.4; // G10 | |
816 | ||
817 | rhoFactor = X0*thickX0/thick; | |
818 | density = rho*rhoFactor; | |
819 | ||
820 | AliMaterial(35,"G10-fc",a,z,density,999.,999.); | |
821 | ||
822 | AliMedium(8,"G10-fc",35,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01); | |
823 | ||
824 | thickX0 = 0.0027; // inner vessel (eta <0.9) | |
825 | thick=0.5; | |
826 | rhoFactor = X0*thickX0/thick; | |
827 | density = rho*rhoFactor; | |
828 | ||
829 | AliMaterial(36,"G10-iv",a,z,density,999.,999.); | |
830 | ||
831 | AliMedium(9,"G10-iv",36,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01); | |
832 | ||
833 | // Carbon fibre | |
fe4da5cc | 834 | |
1283eee5 | 835 | gMC->Gfmate((*fIdmate)[33],namate,a,z,rho,X0,absl,buf,nbuf); |
836 | ||
837 | thickX0 = 0.0133; // outer vessel | |
838 | thick=3.0; | |
839 | rhoFactor = X0*thickX0/thick; | |
840 | density = rho*rhoFactor; | |
841 | ||
842 | ||
843 | AliMaterial(37,"C-ov",a,z,density,999.,999.); | |
844 | ||
845 | AliMedium(10,"C-ov",37,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01); | |
846 | ||
847 | thickX0=0.015; // inner vessel (cone, eta > 0.9) | |
848 | thick=1.5; | |
849 | rhoFactor = X0*thickX0/thick; | |
850 | density = rho*rhoFactor; | |
851 | ||
852 | AliMaterial(38,"C-ivc",a,z,density,999.,999.); | |
853 | ||
854 | AliMedium(11,"C-ivc",38,0, ISXFLD, SXMGMX, 10., .1, .1, .001, .01); | |
855 | ||
856 | // | |
857 | ||
858 | AliMedium(12,"CO2",10,0, ISXFLD, SXMGMX, 10., 999.,.1, .001, .001); | |
859 | ||
860 | ||
861 | ||
fe4da5cc | 862 | } |
863 | ||
864 | //_____________________________________________________________________________ | |
865 | struct Bin { | |
866 | const AliTPCdigit *dig; | |
867 | int idx; | |
868 | Bin() {dig=0; idx=-1;} | |
869 | }; | |
870 | ||
871 | struct PreCluster : public AliTPCcluster { | |
872 | const AliTPCdigit* summit; | |
873 | int idx; | |
874 | int cut; | |
875 | int npeaks; | |
8c555625 | 876 | PreCluster(); |
fe4da5cc | 877 | }; |
8c555625 | 878 | PreCluster::PreCluster() : AliTPCcluster() {cut=npeaks=0;} |
879 | ||
fe4da5cc | 880 | |
881 | //_____________________________________________________________________________ | |
882 | static void FindCluster(int i, int j, Bin bins[][MAXTPCTBK+2], PreCluster &c) | |
883 | { | |
884 | // | |
885 | // Find clusters | |
886 | // | |
887 | Bin& b=bins[i][j]; | |
8c555625 | 888 | double q=double(b.dig->fSignal); |
889 | ||
fe4da5cc | 890 | if (q<0) { // digit is at the edge of the pad row |
891 | q=-q; | |
892 | c.cut=1; | |
893 | } | |
894 | if (b.idx >= 0 && b.idx != c.idx) { | |
895 | c.idx=b.idx; | |
896 | c.npeaks++; | |
897 | } | |
898 | ||
899 | if (q > TMath::Abs(c.summit->fSignal)) c.summit=b.dig; | |
900 | ||
901 | c.fY += i*q; | |
902 | c.fZ += j*q; | |
903 | c.fSigmaY2 += i*i*q; | |
904 | c.fSigmaZ2 += j*j*q; | |
905 | c.fQ += q; | |
8c555625 | 906 | |
fe4da5cc | 907 | b.dig = 0; b.idx = c.idx; |
908 | ||
909 | if (bins[i-1][j].dig) FindCluster(i-1,j,bins,c); | |
910 | if (bins[i][j-1].dig) FindCluster(i,j-1,bins,c); | |
911 | if (bins[i+1][j].dig) FindCluster(i+1,j,bins,c); | |
912 | if (bins[i][j+1].dig) FindCluster(i,j+1,bins,c); | |
913 | ||
914 | } | |
915 | ||
916 | //_____________________________________________________________________________ | |
917 | void AliTPC::Digits2Clusters() | |
918 | { | |
919 | // | |
920 | // simple TPC cluster finder from digits. | |
921 | // | |
8c555625 | 922 | // |
923 | AliTPCParam * fTPCParam = &(fDigParam->GetParam()); | |
924 | ||
fe4da5cc | 925 | const Int_t MAX_PAD=200+2, MAX_BUCKET=MAXTPCTBK+2; |
8c555625 | 926 | const Int_t Q_min=60; |
927 | const Int_t THRESHOLD=20; | |
fe4da5cc | 928 | |
8c555625 | 929 | TTree *t=(TTree*)gDirectory->Get("TreeD0_Param1"); |
930 | t->GetBranch("Digits")->SetAddress(&fDigits); | |
fe4da5cc | 931 | Int_t sectors_by_rows=(Int_t)t->GetEntries(); |
932 | ||
933 | int ncls=0; | |
934 | ||
935 | for (Int_t n=0; n<sectors_by_rows; n++) { | |
936 | if (!t->GetEvent(n)) continue; | |
937 | Bin bins[MAX_PAD][MAX_BUCKET]; | |
938 | AliTPCdigit *dig=(AliTPCdigit*)fDigits->UncheckedAt(0); | |
939 | Int_t nsec=dig->fSector, nrow=dig->fPadRow; | |
940 | Int_t ndigits=fDigits->GetEntriesFast(); | |
941 | ||
942 | int npads; int sign_z; | |
943 | if (nsec<25) { | |
944 | sign_z=(nsec<13) ? 1 : -1; | |
8c555625 | 945 | npads=fTPCParam->GetNPadsLow(nrow-1); |
fe4da5cc | 946 | } else { |
947 | sign_z=(nsec<49) ? 1 : -1; | |
8c555625 | 948 | npads=fTPCParam->GetNPadsUp(nrow-1); |
fe4da5cc | 949 | } |
950 | ||
951 | int ndig; | |
952 | for (ndig=0; ndig<ndigits; ndig++) { | |
953 | dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig); | |
954 | int i=dig->fPad, j=dig->fTime; | |
955 | if (dig->fSignal >= THRESHOLD) bins[i][j].dig=dig; | |
956 | if (i==1 || i==npads || j==1 || j==MAXTPCTBK) dig->fSignal*=-1; | |
957 | } | |
958 | ||
959 | int ncl=0; | |
960 | int i,j; | |
8c555625 | 961 | |
fe4da5cc | 962 | for (i=1; i<MAX_PAD-1; i++) { |
963 | for (j=1; j<MAX_BUCKET-1; j++) { | |
964 | if (bins[i][j].dig == 0) continue; | |
965 | PreCluster c; c.summit=bins[i][j].dig; c.idx=ncls; | |
966 | FindCluster(i, j, bins, c); | |
fe4da5cc | 967 | c.fY /= c.fQ; |
968 | c.fZ /= c.fQ; | |
8c555625 | 969 | |
970 | double s2 = c.fSigmaY2/c.fQ - c.fY*c.fY; | |
971 | c.fSigmaY2 = s2 + 1./12.; | |
972 | c.fSigmaY2 *= fTPCParam->GetPadPitchWidth()* | |
973 | fTPCParam->GetPadPitchWidth(); | |
974 | if (s2 != 0.) c.fSigmaY2 *= 0.022*8*4; | |
975 | ||
976 | s2 = c.fSigmaZ2/c.fQ - c.fZ*c.fZ; | |
977 | c.fSigmaZ2 = s2 + 1./12.; | |
978 | c.fSigmaZ2 *= fTPCParam->GetZWidth()*fTPCParam->GetZWidth(); | |
979 | if (s2 != 0.) c.fSigmaZ2 *= 0.068*4*4; | |
980 | ||
981 | c.fY = (c.fY - 0.5 - 0.5*npads)*fTPCParam->GetPadPitchWidth(); | |
982 | c.fZ = fTPCParam->GetZWidth()*(c.fZ+1); | |
983 | c.fZ -= 3.*fTPCParam->GetZSigma(); // PASA delay | |
fe4da5cc | 984 | c.fZ = sign_z*(z_end - c.fZ); |
8c555625 | 985 | //c.fZ += 0.023; |
986 | c.fSector=nsec; | |
987 | c.fPadRow=nrow; | |
fe4da5cc | 988 | c.fTracks[0]=c.summit->fTracks[0]; |
989 | c.fTracks[1]=c.summit->fTracks[1]; | |
990 | c.fTracks[2]=c.summit->fTracks[2]; | |
8c555625 | 991 | |
fe4da5cc | 992 | if (c.cut) { |
993 | c.fSigmaY2 *= 25.; | |
994 | c.fSigmaZ2 *= 4.; | |
995 | } | |
996 | ||
997 | AddCluster(c); ncls++; ncl++; | |
998 | } | |
999 | } | |
1000 | ||
1001 | for (ndig=0; ndig<ndigits; ndig++) { | |
1002 | dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig); | |
8c555625 | 1003 | if (TMath::Abs(dig->fSignal) >= 0) |
fe4da5cc | 1004 | bins[dig->fPad][dig->fTime].dig=dig; |
1005 | } | |
1006 | ||
1007 | for (i=1; i<MAX_PAD-1; i++) { | |
1008 | for (j=1; j<MAX_BUCKET-1; j++) { | |
1009 | if (bins[i][j].dig == 0) continue; | |
1010 | PreCluster c; c.summit=bins[i][j].dig; c.idx=ncls; | |
1011 | FindCluster(i, j, bins, c); | |
1012 | if (c.fQ <= Q_min) continue; //noise cluster | |
8c555625 | 1013 | if (c.npeaks>1) continue; //overlapped cluster |
fe4da5cc | 1014 | c.fY /= c.fQ; |
1015 | c.fZ /= c.fQ; | |
8c555625 | 1016 | |
1017 | double s2 = c.fSigmaY2/c.fQ - c.fY*c.fY; | |
1018 | c.fSigmaY2 = s2 + 1./12.; | |
1019 | c.fSigmaY2 *= fTPCParam->GetPadPitchWidth()* | |
1020 | fTPCParam->GetPadPitchWidth(); | |
1021 | if (s2 != 0.) c.fSigmaY2 *= 0.022*4*0.6*4; | |
1022 | ||
1023 | s2 = c.fSigmaZ2/c.fQ - c.fZ*c.fZ; | |
1024 | c.fSigmaZ2 = s2 + 1./12.; | |
1025 | c.fSigmaZ2 *= fTPCParam->GetZWidth()*fTPCParam->GetZWidth(); | |
1026 | if (s2 != 0.) c.fSigmaZ2 *= 0.068*4*0.4; | |
1027 | ||
1028 | c.fY = (c.fY - 0.5 - 0.5*npads)*fTPCParam->GetPadPitchWidth(); | |
1029 | c.fZ = fTPCParam->GetZWidth()*(c.fZ+1); | |
1030 | c.fZ -= 3.*fTPCParam->GetZSigma(); // PASA delay | |
fe4da5cc | 1031 | c.fZ = sign_z*(z_end - c.fZ); |
8c555625 | 1032 | //c.fZ += 0.023; |
1033 | c.fSector=nsec; | |
1034 | c.fPadRow=nrow; | |
fe4da5cc | 1035 | c.fTracks[0]=c.summit->fTracks[0]; |
1036 | c.fTracks[1]=c.summit->fTracks[1]; | |
1037 | c.fTracks[2]=c.summit->fTracks[2]; | |
1038 | ||
1039 | if (c.cut) { | |
1040 | c.fSigmaY2 *= 25.; | |
1041 | c.fSigmaZ2 *= 4.; | |
1042 | } | |
1043 | ||
1044 | if (c.npeaks==0) {AddCluster(c); ncls++; ncl++;} | |
1045 | else { | |
1046 | new ((*fClusters)[c.idx]) AliTPCcluster(c); | |
1047 | } | |
1048 | } | |
1049 | } | |
1050 | ||
1051 | cerr<<"sector, row, digits, clusters: " | |
1052 | <<nsec<<' '<<nrow<<' '<<ndigits<<' '<<ncl<<" \r"; | |
1053 | ||
1054 | fDigits->Clear(); | |
1055 | ||
1056 | } | |
1057 | } | |
1058 | ||
1059 | //_____________________________________________________________________________ | |
1060 | void AliTPC::ElDiff(Float_t *xyz) | |
1061 | { | |
8c555625 | 1062 | //-------------------------------------------------- |
fe4da5cc | 1063 | // calculates the diffusion of a single electron |
8c555625 | 1064 | //-------------------------------------------------- |
1065 | ||
1066 | //----------------------------------------------------------------- | |
1067 | // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl | |
1068 | //----------------------------------------------------------------- | |
1069 | AliTPCParam * fTPCParam = &(fDigParam->GetParam()); | |
fe4da5cc | 1070 | Float_t driftl; |
8c555625 | 1071 | |
fe4da5cc | 1072 | Float_t z0=xyz[2]; |
8c555625 | 1073 | |
fe4da5cc | 1074 | driftl=z_end-TMath::Abs(xyz[2]); |
8c555625 | 1075 | |
fe4da5cc | 1076 | if(driftl<0.01) driftl=0.01; |
8c555625 | 1077 | |
1078 | // check the attachment | |
1079 | ||
fe4da5cc | 1080 | driftl=TMath::Sqrt(driftl); |
8c555625 | 1081 | |
1082 | // Float_t sig_t = driftl*diff_t; | |
1083 | //Float_t sig_l = driftl*diff_l; | |
1084 | Float_t sig_t = driftl*fTPCParam->GetDiffT(); | |
1085 | Float_t sig_l = driftl*fTPCParam->GetDiffL(); | |
fe4da5cc | 1086 | xyz[0]=gRandom->Gaus(xyz[0],sig_t); |
1087 | xyz[1]=gRandom->Gaus(xyz[1],sig_t); | |
1088 | xyz[2]=gRandom->Gaus(xyz[2],sig_l); | |
8c555625 | 1089 | |
fe4da5cc | 1090 | if (TMath::Abs(xyz[2])>z_end){ |
8c555625 | 1091 | xyz[2]=TMath::Sign(z_end,z0); |
fe4da5cc | 1092 | } |
1093 | if(xyz[2]*z0 < 0.){ | |
8c555625 | 1094 | Float_t eps = 0.0001; |
1095 | xyz[2]=TMath::Sign(eps,z0); | |
fe4da5cc | 1096 | } |
1097 | } | |
1098 | ||
1099 | //_____________________________________________________________________________ | |
1100 | void AliTPC::Hits2Clusters() | |
1101 | { | |
8c555625 | 1102 | //-------------------------------------------------------- |
fe4da5cc | 1103 | // TPC simple cluster generator from hits |
1104 | // obtained from the TPC Fast Simulator | |
8c555625 | 1105 | // The point errors are taken from the parametrization |
1106 | //-------------------------------------------------------- | |
1107 | ||
1108 | //----------------------------------------------------------------- | |
1109 | // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl | |
1110 | //----------------------------------------------------------------- | |
1111 | AliTPCParam * fTPCParam = &(fDigParam->GetParam()); | |
fe4da5cc | 1112 | Float_t sigma_rphi,sigma_z,cl_rphi,cl_z; |
1113 | // | |
1578254f | 1114 | TParticle *particle; // pointer to a given particle |
fe4da5cc | 1115 | AliTPChit *tpcHit; // pointer to a sigle TPC hit |
1116 | TClonesArray *Particles; //pointer to the particle list | |
1117 | Int_t sector,nhits; | |
1118 | Int_t ipart; | |
1119 | Float_t xyz[5]; | |
1120 | Float_t pl,pt,tanth,rpad,ratio; | |
fe4da5cc | 1121 | Float_t cph,sph; |
1122 | ||
1123 | //--------------------------------------------------------------- | |
1124 | // Get the access to the tracks | |
1125 | //--------------------------------------------------------------- | |
1126 | ||
1127 | TTree *TH = gAlice->TreeH(); | |
1128 | Stat_t ntracks = TH->GetEntries(); | |
1129 | ||
1130 | //------------------------------------------------------------ | |
1131 | // Loop over all sectors (72 sectors) | |
1132 | // Sectors 1-24 are lower sectors, 1-12 z>0, 13-24 z<0 | |
1133 | // Sectors 25-72 are upper sectors, 25-48 z>0, 49-72 z<0 | |
1134 | // | |
1135 | // First cluster for sector 1 starts at "0" | |
1136 | //------------------------------------------------------------ | |
1137 | ||
1138 | ||
1139 | fClustersIndex[0] = 0; | |
1140 | ||
1141 | // | |
1142 | for(Int_t isec=1;isec<fNsectors+1;isec++){ | |
8c555625 | 1143 | //MI change |
1144 | fTPCParam->AdjustAngles(isec,cph,sph); | |
fe4da5cc | 1145 | |
1146 | //------------------------------------------------------------ | |
1147 | // Loop over tracks | |
1148 | //------------------------------------------------------------ | |
1149 | ||
1150 | for(Int_t track=0;track<ntracks;track++){ | |
1151 | ResetHits(); | |
1152 | TH->GetEvent(track); | |
1153 | // | |
1154 | // Get number of the TPC hits and a pointer | |
1155 | // to the particles | |
1156 | // | |
1157 | nhits=fHits->GetEntriesFast(); | |
1158 | Particles=gAlice->Particles(); | |
1159 | // | |
1160 | // Loop over hits | |
1161 | // | |
1162 | for(Int_t hit=0;hit<nhits;hit++){ | |
1163 | tpcHit=(AliTPChit*)fHits->UncheckedAt(hit); | |
1164 | sector=tpcHit->fSector; // sector number | |
1165 | if(sector != isec) continue; //terminate iteration | |
1166 | ipart=tpcHit->fTrack; | |
1578254f | 1167 | particle=(TParticle*)Particles->UncheckedAt(ipart); |
1168 | pl=particle->Pz(); | |
1169 | pt=particle->Pt(); | |
fe4da5cc | 1170 | if(pt < 1.e-9) pt=1.e-9; |
1171 | tanth=pl/pt; | |
1172 | tanth = TMath::Abs(tanth); | |
1173 | rpad=TMath::Sqrt(tpcHit->fX*tpcHit->fX + tpcHit->fY*tpcHit->fY); | |
1174 | ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason | |
1175 | ||
1176 | // space-point resolutions | |
1177 | ||
1178 | sigma_rphi=SigmaY2(rpad,tanth,pt); | |
1179 | sigma_z =SigmaZ2(rpad,tanth ); | |
1180 | ||
1181 | // cluster widths | |
1182 | ||
1183 | cl_rphi=ac_rphi-bc_rphi*rpad*tanth+cc_rphi*ratio*ratio; | |
1184 | cl_z=ac_z-bc_z*rpad*tanth+cc_z*tanth*tanth; | |
1185 | ||
1186 | // temporary protection | |
1187 | ||
1188 | if(sigma_rphi < 0.) sigma_rphi=0.4e-3; | |
1189 | if(sigma_z < 0.) sigma_z=0.4e-3; | |
1190 | if(cl_rphi < 0.) cl_rphi=2.5e-3; | |
1191 | if(cl_z < 0.) cl_z=2.5e-5; | |
1192 | ||
1193 | // | |
1194 | ||
1195 | // | |
1196 | // smearing --> rotate to the 1 (13) or to the 25 (49) sector, | |
1197 | // then the inaccuracy in a X-Y plane is only along Y (pad row)! | |
1198 | // | |
1199 | Float_t xprim= tpcHit->fX*cph + tpcHit->fY*sph; | |
1200 | Float_t yprim=-tpcHit->fX*sph + tpcHit->fY*cph; | |
1201 | xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigma_rphi)); // y | |
8c555625 | 1202 | Double_t alpha=(sector<25) ? alpha_low : alpha_up; |
1203 | if (TMath::Abs(xyz[0]/xprim) > TMath::Tan(0.5*alpha)) xyz[0]=yprim; | |
fe4da5cc | 1204 | xyz[1]=gRandom->Gaus(tpcHit->fZ,TMath::Sqrt(sigma_z)); // z |
8c555625 | 1205 | if (TMath::Abs(xyz[1]) > 250) xyz[1]=tpcHit->fZ; |
fe4da5cc | 1206 | xyz[2]=tpcHit->fQ; // q |
1207 | xyz[3]=sigma_rphi; // fSigmaY2 | |
1208 | xyz[4]=sigma_z; // fSigmaZ2 | |
1283eee5 | 1209 | |
fe4da5cc | 1210 | // and finally add the cluster |
1283eee5 | 1211 | Int_t tracks[5]={tpcHit->fTrack, -1, -1, sector, tpcHit->fPadRow}; |
fe4da5cc | 1212 | AddCluster(xyz,tracks); |
1213 | ||
1214 | } // end of loop over hits | |
1215 | } // end of loop over tracks | |
1216 | ||
1217 | fClustersIndex[isec] = fNclusters; // update clusters index | |
1218 | ||
1219 | } // end of loop over sectors | |
1220 | ||
1221 | fClustersIndex[fNsectors]--; // set end of the clusters buffer | |
1222 | ||
1223 | } | |
1224 | ||
8c555625 | 1225 | |
1226 | void AliTPC::Hits2Digits() | |
1227 | { | |
1228 | ||
1229 | //---------------------------------------------------- | |
1230 | // Loop over all sectors (72 sectors) | |
1231 | // Sectors 1-24 are lower sectors, 1-12 z>0, 13-24 z<0 | |
1232 | // Sectors 25-72 are upper sectors, 25-48 z>0, 49-72 z<0 | |
1233 | //---- | |
1283eee5 | 1234 | for(Int_t isec=0;isec<fNsectors;isec++) Hits2DigitsSector(isec); |
8c555625 | 1235 | } |
1236 | ||
1237 | ||
fe4da5cc | 1238 | //_____________________________________________________________________________ |
8c555625 | 1239 | void AliTPC::Hits2DigitsSector(Int_t isec) |
fe4da5cc | 1240 | { |
8c555625 | 1241 | //------------------------------------------------------------------- |
fe4da5cc | 1242 | // TPC conversion from hits to digits. |
8c555625 | 1243 | //------------------------------------------------------------------- |
1244 | ||
1245 | //----------------------------------------------------------------- | |
1246 | // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl | |
1247 | //----------------------------------------------------------------- | |
1248 | ||
fe4da5cc | 1249 | //------------------------------------------------------- |
8c555625 | 1250 | // Get the access to the track hits |
fe4da5cc | 1251 | //------------------------------------------------------- |
8c555625 | 1252 | |
1253 | AliTPCParam * fTPCParam = &(fDigParam->GetParam()); | |
1254 | TTree *TH = gAlice->TreeH(); // pointer to the hits tree | |
fe4da5cc | 1255 | Stat_t ntracks = TH->GetEntries(); |
8c555625 | 1256 | |
1257 | if( ntracks > 0){ | |
1258 | ||
1259 | //------------------------------------------- | |
1260 | // Only if there are any tracks... | |
1261 | //------------------------------------------- | |
1262 | ||
fe4da5cc | 1263 | |
8c555625 | 1264 | // TObjArrays for three neighouring pad-rows |
1265 | ||
1266 | TObjArray **rowTriplet = new TObjArray* [3]; | |
fe4da5cc | 1267 | |
8c555625 | 1268 | // TObjArray-s for each pad-row |
1269 | ||
1270 | TObjArray **row; | |
1271 | ||
1272 | for (Int_t trip=0;trip<3;trip++){ | |
1273 | rowTriplet[trip]=new TObjArray; | |
fe4da5cc | 1274 | } |
8c555625 | 1275 | |
1276 | ||
fe4da5cc | 1277 | |
8c555625 | 1278 | printf("*** Processing sector number %d ***\n",isec); |
1279 | ||
1280 | Int_t nrows =fTPCParam->GetNRow(isec); | |
1281 | ||
1282 | row= new TObjArray* [nrows]; | |
fe4da5cc | 1283 | |
8c555625 | 1284 | MakeSector(isec,nrows,TH,ntracks,row); |
1285 | ||
1286 | //-------------------------------------------------------- | |
1287 | // Digitize this sector, row by row | |
1288 | // row[i] is the pointer to the TObjArray of TVectors, | |
1289 | // each one containing electrons accepted on this | |
1290 | // row, assigned into tracks | |
1291 | //-------------------------------------------------------- | |
1292 | ||
1293 | Int_t i; | |
1294 | ||
1295 | for (i=0;i<nrows;i++){ | |
1296 | ||
1297 | // Triplets for i = 0 and i=1 are identical! | |
1298 | // The same for i = nrows-1 and nrows! | |
1299 | ||
1300 | if(i != 1 && i != nrows-1){ | |
1301 | MakeTriplet(i,rowTriplet,row); | |
1302 | } | |
1303 | ||
1304 | DigitizeRow(i,isec,rowTriplet); | |
1305 | ||
1306 | fDigParam->Fill(); | |
1307 | ||
1308 | Int_t ndig=fDigParam->GetArray()->GetEntriesFast(); | |
1309 | ||
1310 | printf("*** Sector, row, digits %d %d %d ***\n",isec,i,ndig); | |
1311 | ||
1312 | ResetDigits(); // reset digits for this row after storing them | |
1313 | ||
1314 | } // end of the sector digitization | |
1315 | ||
1316 | // delete the last triplet | |
1317 | ||
1318 | for (i=0;i<3;i++) rowTriplet[i]->Delete(); | |
1319 | ||
1320 | delete [] row; // delete the array of pointers to TObjArray-s | |
1321 | ||
1322 | } // ntracks >0 | |
1323 | } // end of Hits2Digits | |
1324 | //_____________________________________________________________________________ | |
1325 | void AliTPC::MakeTriplet(Int_t row, | |
1326 | TObjArray **rowTriplet, TObjArray **prow) | |
1327 | { | |
1328 | //------------------------------------------------------------------ | |
1329 | // Makes the "triplet" of the neighbouring pad-row for the | |
1330 | // digitization including the cross-talk between the pad-rows | |
1331 | //------------------------------------------------------------------ | |
1332 | ||
1333 | //----------------------------------------------------------------- | |
1334 | // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl | |
1335 | //----------------------------------------------------------------- | |
1336 | ||
1337 | AliTPCParam * fTPCParam = &(fDigParam->GetParam()); | |
1338 | Float_t gasgain = fTPCParam->GetGasGain(); | |
1339 | Int_t nTracks[3]; | |
1340 | ||
1341 | Int_t nElements,nElectrons; | |
1342 | ||
1343 | TVector *pv; | |
1344 | TVector *tv; | |
1345 | ||
1346 | //------------------------------------------------------------------- | |
1347 | // pv is an "old" track, i.e. label + triplets of (x,y,z) | |
1348 | // for each electron | |
1349 | // | |
1350 | //------------------------------------------------------------------- | |
1351 | ||
1352 | ||
1353 | Int_t i1,i2; | |
1354 | Int_t nel,nt; | |
1355 | ||
1356 | if(row == 0 || row == 1){ | |
1357 | ||
1358 | // create entire triplet for the first AND the second row | |
1359 | ||
1360 | nTracks[0] = prow[0]->GetEntries(); | |
1361 | nTracks[1] = prow[1]->GetEntries(); | |
1362 | nTracks[2] = prow[2]->GetEntries(); | |
1363 | ||
1364 | for(i1=0;i1<3;i1++){ | |
1365 | nt = nTracks[i1]; // number of tracks for this row | |
1366 | ||
1367 | for(i2=0;i2<nt;i2++){ | |
1368 | pv = (TVector*)prow[i1]->At(i2); | |
1369 | TVector &v1 = *pv; | |
1370 | nElements = pv->GetNrows(); | |
1371 | nElectrons = (nElements-1)/3; | |
1372 | ||
1373 | tv = new TVector(4*nElectrons+1); // create TVector for a modified track | |
1374 | TVector &v2 = *tv; | |
1375 | v2(0)=v1(0); //track label | |
1376 | ||
1377 | for(nel=0;nel<nElectrons;nel++){ | |
1378 | Int_t idx1 = nel*3; | |
1379 | Int_t idx2 = nel*4; | |
1380 | // Avalanche, including fluctuations | |
1381 | Int_t aval = (Int_t) (-gasgain*TMath::Log(gRandom->Rndm())); | |
1382 | v2(idx2+1)= v1(idx1+1); | |
1383 | v2(idx2+2)= v1(idx1+2); | |
1384 | v2(idx2+3)= v1(idx1+3); | |
1385 | v2(idx2+4)= (Float_t)aval; // in number of electrons! | |
1386 | } // end of loop over electrons | |
1387 | // | |
1388 | // Add this track to a row | |
1389 | // | |
1390 | ||
1391 | rowTriplet[i1]->Add(tv); | |
1392 | ||
1393 | ||
1394 | } // end of loop over tracks for this row | |
1395 | ||
1396 | prow[i1]->Delete(); // remove "old" tracks | |
1397 | delete prow[i1]; // delete TObjArray for this row | |
1398 | prow[i1]=0; // set pointer to NULL | |
1399 | ||
1400 | } // end of loop over row triplets | |
1401 | ||
1402 | ||
1403 | } | |
1404 | else{ | |
1405 | ||
1406 | rowTriplet[0]->Delete(); // remove old lower row | |
1407 | ||
1408 | nTracks[0]=rowTriplet[1]->GetEntries(); // previous middle row | |
1409 | nTracks[1]=rowTriplet[2]->GetEntries(); // previous upper row | |
1410 | nTracks[2]=prow[row+1]->GetEntries(); // next row | |
1411 | ||
1412 | ||
1413 | //------------------------------------------- | |
1414 | // shift new tracks downwards | |
1415 | //------------------------------------------- | |
1416 | ||
1417 | for(i1=0;i1<nTracks[0];i1++){ | |
1418 | pv=(TVector*)rowTriplet[1]->At(i1); | |
1419 | rowTriplet[0]->Add(pv); | |
1420 | } | |
1421 | ||
1422 | rowTriplet[1]->Clear(); // leave tracks on the heap!!! | |
1423 | ||
1424 | for(i1=0;i1<nTracks[1];i1++){ | |
1425 | pv=(TVector*)rowTriplet[2]->At(i1); | |
1426 | rowTriplet[1]->Add(pv); | |
1427 | } | |
1428 | ||
1429 | rowTriplet[2]->Clear(); // leave tracks on the heap!!! | |
1430 | ||
1431 | //--------------------------------------------- | |
1432 | // Create new upper row | |
1433 | //--------------------------------------------- | |
1434 | ||
1435 | ||
1436 | ||
1437 | for(i1=0;i1<nTracks[2];i1++){ | |
1438 | pv = (TVector*)prow[row+1]->At(i1); | |
1439 | TVector &v1 = *pv; | |
1440 | nElements = pv->GetNrows(); | |
1441 | nElectrons = (nElements-1)/3; | |
1442 | ||
1443 | tv = new TVector(4*nElectrons+1); // create TVector for a modified track | |
1444 | TVector &v2 = *tv; | |
1445 | v2(0)=v1(0); //track label | |
1446 | ||
1447 | for(nel=0;nel<nElectrons;nel++){ | |
1448 | ||
1449 | Int_t idx1 = nel*3; | |
1450 | Int_t idx2 = nel*4; | |
1451 | // Avalanche, including fluctuations | |
1452 | Int_t aval = (Int_t) | |
1453 | (-gasgain*TMath::Log(gRandom->Rndm())); | |
1454 | ||
1455 | v2(idx2+1)= v1(idx1+1); | |
1456 | v2(idx2+2)= v1(idx1+2); | |
1457 | v2(idx2+3)= v1(idx1+3); | |
1458 | v2(idx2+4)= (Float_t)aval; // in number of electrons! | |
1459 | } // end of loop over electrons | |
1460 | ||
1461 | rowTriplet[2]->Add(tv); | |
1462 | ||
1463 | } // end of loop over tracks | |
1464 | ||
1465 | prow[row+1]->Delete(); // delete tracks for this row | |
1466 | delete prow[row+1]; // delete TObjArray for this row | |
1467 | prow[row+1]=0; // set a pointer to NULL | |
1468 | ||
1469 | } | |
1470 | ||
1471 | } // end of MakeTriplet | |
1472 | //_____________________________________________________________________________ | |
1473 | void AliTPC::ExB(Float_t *xyz) | |
1474 | { | |
1475 | //------------------------------------------------ | |
1476 | // ExB at the wires and wire number calulation | |
1477 | //------------------------------------------------ | |
1478 | ||
1479 | //----------------------------------------------------------------- | |
1480 | // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl | |
1481 | //----------------------------------------------------------------- | |
1482 | AliTPCParam * fTPCParam = &(fDigParam->GetParam()); | |
1483 | ||
1484 | Float_t x1=xyz[0]; | |
1485 | fTPCParam->GetWire(x1); //calculate nearest wire position | |
1486 | Float_t dx=xyz[0]-x1; | |
1487 | xyz[1]+=dx*fTPCParam->GetOmegaTau(); | |
1488 | ||
1489 | } // end of ExB | |
1490 | //_____________________________________________________________________________ | |
1491 | void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rowTriplet) | |
1492 | { | |
1493 | //----------------------------------------------------------- | |
1494 | // Single row digitization, coupling from the neighbouring | |
1495 | // rows taken into account | |
1496 | //----------------------------------------------------------- | |
1497 | ||
1498 | //----------------------------------------------------------------- | |
1499 | // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl | |
1500 | //----------------------------------------------------------------- | |
1501 | ||
1502 | ||
1503 | AliTPCParam * fTPCParam = &(fDigParam->GetParam()); | |
1504 | Float_t chipgain= fTPCParam->GetChipGain(); | |
1505 | Float_t zerosup = fTPCParam->GetZeroSup(); | |
1506 | Int_t nrows =fTPCParam->GetNRow(isec); | |
1507 | ||
1508 | Int_t nTracks[3]; | |
1509 | Int_t n_of_pads[3]; | |
1510 | Int_t IndexRange[4]; | |
1511 | Int_t i1; | |
1512 | Int_t iFlag; | |
1513 | ||
1514 | // | |
1515 | // iFlag = 0 -> inner row, iFlag = 1 -> the middle one, iFlag = 2 -> the outer one | |
1516 | // | |
1517 | ||
1518 | nTracks[0]=rowTriplet[0]->GetEntries(); // lower row | |
1519 | nTracks[1]=rowTriplet[1]->GetEntries(); // middle row | |
1520 | nTracks[2]=rowTriplet[2]->GetEntries(); // upper row | |
1521 | ||
1522 | ||
1523 | if(irow == 0){ | |
1524 | iFlag=0; | |
1525 | n_of_pads[0]=fTPCParam->GetNPads(isec,0); | |
1526 | n_of_pads[1]=fTPCParam->GetNPads(isec,1); | |
1527 | } | |
1528 | else if(irow == nrows-1){ | |
1529 | iFlag=2; | |
1530 | n_of_pads[1]=fTPCParam->GetNPads(isec,irow-1); | |
1531 | n_of_pads[2]=fTPCParam->GetNPads(isec,irow); | |
1532 | } | |
1533 | else { | |
1534 | iFlag=1; | |
1535 | for(i1=0;i1<3;i1++){ | |
1536 | n_of_pads[i1]=fTPCParam->GetNPads(isec,irow-1+i1); | |
1537 | } | |
1538 | } | |
1539 | ||
1540 | // | |
1541 | // Integrated signal for this row | |
1542 | // and a single track signal | |
1543 | // | |
1544 | ||
1283eee5 | 1545 | TMatrix *m1 = new TMatrix(0,n_of_pads[iFlag]-1,0,MAXTPCTBK-1); // integrated |
1546 | TMatrix *m2 = new TMatrix(0,n_of_pads[iFlag]-1,0,MAXTPCTBK-1); // single | |
8c555625 | 1547 | |
1548 | // | |
1549 | ||
1550 | TMatrix &Total = *m1; | |
1551 | ||
1552 | // Array of pointers to the label-signal list | |
1553 | ||
1554 | Int_t NofDigits = n_of_pads[iFlag]*MAXTPCTBK; // number of digits for this row | |
1555 | ||
1556 | Float_t **pList = new Float_t* [NofDigits]; | |
1557 | ||
1558 | Int_t lp; | |
1559 | ||
1560 | for(lp=0;lp<NofDigits;lp++)pList[lp]=0; // set all pointers to NULL | |
1561 | ||
1562 | // | |
1563 | // Straight signal and cross-talk, cross-talk is integrated over all | |
1564 | // tracks and added to the signal at the very end | |
1565 | // | |
1566 | ||
1567 | ||
1568 | for (i1=0;i1<nTracks[iFlag];i1++){ | |
1569 | ||
1570 | m2->Zero(); // clear single track signal matrix | |
1571 | ||
1572 | Float_t TrackLabel = | |
1573 | GetSignal(rowTriplet[iFlag],i1,n_of_pads[iFlag],m2,m1,IndexRange); | |
1574 | ||
1575 | GetList(TrackLabel,n_of_pads[iFlag],m2,IndexRange,pList); | |
1576 | ||
1577 | } | |
1578 | ||
1579 | // | |
1580 | // Cross talk from the neighbouring pad-rows | |
1581 | // | |
1582 | ||
1283eee5 | 1583 | TMatrix *m3 = new TMatrix(0,n_of_pads[iFlag]-1,0,MAXTPCTBK-1); // cross-talk |
8c555625 | 1584 | |
1585 | TMatrix &Cross = *m3; | |
1586 | ||
1587 | if(iFlag == 0){ | |
1588 | ||
1589 | // cross-talk from the outer row only (first pad row) | |
1590 | ||
1591 | GetCrossTalk(0,rowTriplet[1],nTracks[1],n_of_pads,m3); | |
1592 | ||
1593 | } | |
1594 | else if(iFlag == 2){ | |
1595 | ||
1596 | // cross-talk from the inner row only (last pad row) | |
1597 | ||
1598 | GetCrossTalk(2,rowTriplet[1],nTracks[1],n_of_pads,m3); | |
1599 | ||
1600 | } | |
1601 | else{ | |
1602 | ||
1603 | // cross-talk from both inner and outer rows | |
1604 | ||
1605 | GetCrossTalk(3,rowTriplet[0],nTracks[0],n_of_pads,m3); // inner | |
1606 | GetCrossTalk(4,rowTriplet[2],nTracks[2],n_of_pads,m3); //outer | |
1607 | } | |
1608 | ||
1609 | Total += Cross; // add the cross-talk | |
1610 | ||
1611 | // | |
1612 | // Convert analog signal to ADC counts | |
1613 | // | |
1614 | ||
1615 | Int_t tracks[3]; | |
1616 | Int_t digits[5]; | |
1617 | ||
1618 | ||
1283eee5 | 1619 | for(Int_t ip=0;ip<n_of_pads[iFlag];ip++){ |
1620 | for(Int_t it=0;it<MAXTPCTBK;it++){ | |
8c555625 | 1621 | |
1622 | Float_t q = Total(ip,it); | |
1623 | ||
1283eee5 | 1624 | Int_t gi =it*n_of_pads[iFlag]+ip; // global index |
8c555625 | 1625 | |
1626 | q = gRandom->Gaus(q,fTPCParam->GetNoise()); // apply noise | |
1627 | q *= (q_el*1.e15); // convert to fC | |
1628 | q *= chipgain; // convert to mV | |
1629 | q *= (adc_sat/dyn_range); // convert to ADC counts | |
1630 | ||
1631 | if(q <zerosup) continue; // do not fill zeros | |
1632 | if(q > adc_sat) q = adc_sat; // saturation | |
1633 | ||
1634 | // | |
1635 | // "real" signal or electronic noise (list = -1)? | |
1636 | // | |
1637 | ||
1638 | for(Int_t j1=0;j1<3;j1++){ | |
1639 | tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -1; | |
1640 | } | |
1641 | ||
1642 | digits[0]=isec; | |
1283eee5 | 1643 | digits[1]=irow; |
8c555625 | 1644 | digits[2]=ip; |
1645 | digits[3]=it; | |
1646 | digits[4]= (Int_t)q; | |
1647 | ||
1648 | // Add this digit | |
1649 | ||
1650 | AddDigit(tracks,digits); | |
1651 | ||
1652 | } // end of loop over time buckets | |
1653 | } // end of lop over pads | |
1654 | ||
1655 | // | |
1656 | // This row has been digitized, delete nonused stuff | |
1657 | // | |
1658 | ||
1659 | for(lp=0;lp<NofDigits;lp++){ | |
1660 | if(pList[lp]) delete [] pList[lp]; | |
1661 | } | |
1662 | ||
1663 | delete [] pList; | |
1664 | ||
1665 | delete m1; | |
1666 | delete m2; | |
1667 | delete m3; | |
1668 | ||
1669 | } // end of DigitizeRow | |
1670 | //_____________________________________________________________________________ | |
1671 | Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, Int_t np, TMatrix *m1, TMatrix *m2, | |
1672 | Int_t *IndexRange) | |
1673 | { | |
1674 | ||
1675 | //--------------------------------------------------------------- | |
1676 | // Calculates 2-D signal (pad,time) for a single track, | |
1677 | // returns a pointer to the signal matrix and the track label | |
1678 | // No digitization is performed at this level!!! | |
1679 | //--------------------------------------------------------------- | |
1680 | ||
1681 | //----------------------------------------------------------------- | |
1682 | // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl | |
1683 | //----------------------------------------------------------------- | |
1684 | ||
1685 | TVector *tv; | |
1686 | AliTPCParam * fTPCParam = &(fDigParam->GetParam()); | |
1687 | AliTPCPRF2D * fPRF2D = &(fDigParam->GetPRF2D()); | |
1688 | AliTPCRF1D * fRF = &(fDigParam->GetRF()); | |
1689 | ||
1690 | //to make the code faster we put parameters to the stack | |
1691 | ||
1692 | Float_t zwidth = fTPCParam->GetZWidth(); | |
7ee081c2 | 1693 | Float_t zwidthm1 =1./zwidth; |
8c555625 | 1694 | |
1695 | tv = (TVector*)p1->At(ntr); // pointer to a track | |
1696 | TVector &v = *tv; | |
1697 | ||
1698 | Float_t label = v(0); | |
1699 | ||
1283eee5 | 1700 | Int_t CentralPad = (np-1)/2; |
8c555625 | 1701 | Int_t PadNumber; |
1702 | Int_t nElectrons = (tv->GetNrows()-1)/4; | |
1703 | Float_t range=((np-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth(); // pad range | |
8c555625 | 1704 | |
1283eee5 | 1705 | range -= 0.5; // dead zone, 5mm from the edge, according to H.G. Fischer |
1706 | ||
8c555625 | 1707 | Float_t IneffFactor = 0.5; // inefficiency in the gain close to the edge, as above |
1708 | ||
1709 | ||
1710 | Float_t PadSignal[7]; // signal from a single electron | |
1711 | ||
1712 | TMatrix &signal = *m1; | |
1713 | TMatrix &total = *m2; | |
1714 | ||
1715 | ||
1716 | IndexRange[0]=9999; // min pad | |
1717 | IndexRange[1]=-1; // max pad | |
1718 | IndexRange[2]=9999; //min time | |
1719 | IndexRange[3]=-1; // max time | |
1720 | ||
1721 | // | |
1722 | // Loop over all electrons | |
1723 | // | |
1724 | ||
1725 | for(Int_t nel=0; nel<nElectrons; nel++){ | |
1726 | Int_t idx=nel*4; | |
1727 | Float_t xwire = v(idx+1); | |
1728 | Float_t y = v(idx+2); | |
1729 | Float_t z = v(idx+3); | |
1730 | ||
1731 | ||
1732 | Float_t absy=TMath::Abs(y); | |
1733 | ||
1734 | if(absy < 0.5*fTPCParam->GetPadPitchWidth()){ | |
1735 | PadNumber=CentralPad; | |
1736 | } | |
1737 | else if (absy < range){ | |
8022212f | 1738 | PadNumber=(Int_t) ((absy-0.5*fTPCParam->GetPadPitchWidth())/fTPCParam->GetPadPitchWidth()+1.); |
8c555625 | 1739 | PadNumber=(Int_t) (TMath::Sign((Float_t)PadNumber, y)+CentralPad); |
1740 | } | |
1741 | else continue; // electron out of pad-range , lost at the sector edge | |
1742 | ||
1743 | Float_t aval = (absy<range-0.5) ? v(idx+4):v(idx+4)*IneffFactor; | |
1744 | ||
1745 | ||
1746 | Float_t dist = y - (Float_t)(PadNumber-CentralPad)*fTPCParam->GetPadPitchWidth(); | |
a03540ae | 1747 | for (Int_t i=0;i<7;i++){ |
b6738b36 | 1748 | PadSignal[i]=fPRF2D->GetPRF(-dist+(i-3)*fTPCParam->GetPadPitchWidth(),xwire)*aval; |
a03540ae | 1749 | PadSignal[i] *= fTPCParam->GetPadCoupling(); |
1750 | } | |
8c555625 | 1751 | |
1283eee5 | 1752 | Int_t LeftPad = TMath::Max(0,PadNumber-3); |
1753 | Int_t RightPad = TMath::Min(np-1,PadNumber+3); | |
8c555625 | 1754 | |
1755 | Int_t pmin=LeftPad-PadNumber+3; // lower index of the pad_signal vector | |
1756 | Int_t pmax=RightPad-PadNumber+3; // upper index | |
1757 | ||
1283eee5 | 1758 | Float_t z_drift = z*zwidthm1; |
8c555625 | 1759 | Float_t z_offset = z_drift-(Int_t)z_drift; |
1760 | //distance to the centre of nearest time bin (in time bin units) | |
1283eee5 | 1761 | Int_t FirstBucket = (Int_t)z_drift; |
8c555625 | 1762 | |
1763 | ||
1764 | // loop over time bins (4 bins is enough - 3 sigma truncated Gaussian) | |
1765 | for (Int_t i2=0;i2<4;i2++){ | |
1766 | Int_t TrueTime = FirstBucket+i2; // current time bucket | |
1767 | Float_t dz = (Float_t(i2)+z_offset)*zwidth; | |
1768 | Float_t ampl = fRF->GetRF(dz); | |
1283eee5 | 1769 | if( (TrueTime>MAXTPCTBK-1) ) break; // beyond the time range |
8c555625 | 1770 | |
1771 | IndexRange[2]=TMath::Min(IndexRange[2],TrueTime); // min time | |
1772 | IndexRange[3]=TMath::Max(IndexRange[3],TrueTime); // max time | |
1773 | ||
1774 | // loop over pads, from pmin to pmax | |
1775 | for(Int_t i3=pmin;i3<=pmax;i3++){ | |
1776 | Int_t TruePad = LeftPad+i3-pmin; | |
1777 | IndexRange[0]=TMath::Min(IndexRange[0],TruePad); // min pad | |
1778 | IndexRange[1]=TMath::Max(IndexRange[1],TruePad); // max pad | |
1779 | signal(TruePad,TrueTime)+=(PadSignal[i3]*ampl); // not converted to charge!!! | |
1780 | total(TruePad,TrueTime)+=(PadSignal[i3]*ampl); // not converted to charge!!! | |
1781 | } // end of pads loop | |
1782 | } // end of time loop | |
1783 | } // end of loop over electrons | |
1784 | ||
1785 | return label; // returns track label when finished | |
1786 | } | |
1787 | ||
1788 | //_____________________________________________________________________________ | |
1789 | void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *IndexRange, | |
1790 | Float_t **pList) | |
1791 | { | |
1792 | //---------------------------------------------------------------------- | |
1793 | // Updates the list of tracks contributing to digits for a given row | |
1794 | //---------------------------------------------------------------------- | |
1795 | ||
1796 | //----------------------------------------------------------------- | |
1797 | // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl | |
1798 | //----------------------------------------------------------------- | |
1799 | ||
1800 | TMatrix &signal = *m; | |
1801 | ||
1802 | // lop over nonzero digits | |
1803 | ||
1804 | for(Int_t it=IndexRange[2];it<IndexRange[3]+1;it++){ | |
1805 | for(Int_t ip=IndexRange[0];ip<IndexRange[1]+1;ip++){ | |
1806 | ||
1807 | ||
1283eee5 | 1808 | Int_t GlobalIndex = it*np+ip; // GlobalIndex starts from 0! |
8c555625 | 1809 | |
1810 | if(!pList[GlobalIndex]){ | |
1811 | ||
1812 | // | |
1813 | // Create new list (6 elements - 3 signals and 3 labels), | |
1814 | // but only if the signal is > 0. | |
1815 | // | |
1816 | ||
1817 | if(signal(ip,it)>0.){ | |
1818 | ||
1819 | pList[GlobalIndex] = new Float_t [6]; | |
1820 | ||
1821 | // set list to -1 | |
1822 | ||
1823 | *pList[GlobalIndex] = -1.; | |
1824 | *(pList[GlobalIndex]+1) = -1.; | |
1825 | *(pList[GlobalIndex]+2) = -1.; | |
1826 | *(pList[GlobalIndex]+3) = -1.; | |
1827 | *(pList[GlobalIndex]+4) = -1.; | |
1828 | *(pList[GlobalIndex]+5) = -1.; | |
1829 | ||
1830 | ||
1831 | *pList[GlobalIndex] = label; | |
1832 | *(pList[GlobalIndex]+3) = signal(ip,it);} | |
1833 | } | |
1834 | else{ | |
1835 | ||
1836 | // check the signal magnitude | |
1837 | ||
1838 | Float_t highest = *(pList[GlobalIndex]+3); | |
1839 | Float_t middle = *(pList[GlobalIndex]+4); | |
1840 | Float_t lowest = *(pList[GlobalIndex]+5); | |
1841 | ||
1842 | // | |
1843 | // compare the new signal with already existing list | |
1844 | // | |
1845 | ||
1846 | if(signal(ip,it)<lowest) continue; // neglect this track | |
1847 | ||
1848 | // | |
1849 | ||
1850 | if (signal(ip,it)>highest){ | |
1851 | *(pList[GlobalIndex]+5) = middle; | |
1852 | *(pList[GlobalIndex]+4) = highest; | |
1853 | *(pList[GlobalIndex]+3) = signal(ip,it); | |
1854 | ||
1855 | *(pList[GlobalIndex]+2) = *(pList[GlobalIndex]+1); | |
1856 | *(pList[GlobalIndex]+1) = *pList[GlobalIndex]; | |
1857 | *pList[GlobalIndex] = label; | |
1858 | } | |
1859 | else if (signal(ip,it)>middle){ | |
1860 | *(pList[GlobalIndex]+5) = middle; | |
1861 | *(pList[GlobalIndex]+4) = signal(ip,it); | |
1862 | ||
1863 | *(pList[GlobalIndex]+2) = *(pList[GlobalIndex]+1); | |
1864 | *(pList[GlobalIndex]+1) = label; | |
1865 | } | |
1866 | else{ | |
1867 | *(pList[GlobalIndex]+5) = signal(ip,it); | |
1868 | *(pList[GlobalIndex]+2) = label; | |
1869 | } | |
1870 | } | |
1871 | ||
1872 | } // end of loop over pads | |
1873 | } // end of loop over time bins | |
1874 | ||
1875 | ||
1876 | ||
1877 | ||
1878 | }//end of GetList | |
1879 | //___________________________________________________________________ | |
1880 | void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH, | |
1881 | Stat_t ntracks,TObjArray **row) | |
1882 | { | |
1883 | ||
1884 | //----------------------------------------------------------------- | |
1885 | // Prepares the sector digitization, creates the vectors of | |
1886 | // tracks for each row of this sector. The track vector | |
1887 | // contains the track label and the position of electrons. | |
1888 | //----------------------------------------------------------------- | |
1889 | ||
1890 | //----------------------------------------------------------------- | |
1891 | // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl | |
1892 | //----------------------------------------------------------------- | |
1893 | ||
1894 | AliTPCParam * fTPCParam = &(fDigParam->GetParam()); | |
1895 | Int_t i; | |
1896 | Float_t xyz[3]; | |
1897 | ||
1898 | AliTPChit *tpcHit; // pointer to a sigle TPC hit | |
1899 | ||
1900 | //---------------------------------------------- | |
1901 | // Create TObjArray-s, one for each row, | |
1902 | // each TObjArray will store the TVectors | |
1903 | // of electrons, one TVector per each track. | |
1904 | //---------------------------------------------- | |
1905 | ||
1906 | for(i=0; i<nrows; i++){ | |
1907 | row[i] = new TObjArray; | |
1908 | } | |
1909 | Int_t *n_of_electrons = new Int_t [nrows]; // electron counter for each row | |
1910 | TVector **tr = new TVector* [nrows]; //pointers to the track vectors | |
1911 | ||
1912 | //-------------------------------------------------------------------- | |
1913 | // Loop over tracks, the "track" contains the full history | |
1914 | //-------------------------------------------------------------------- | |
1915 | ||
1916 | Int_t previousTrack,currentTrack; | |
1917 | previousTrack = -1; // nothing to store so far! | |
1918 | ||
1919 | for(Int_t track=0;track<ntracks;track++){ | |
1920 | ||
1921 | ResetHits(); | |
1922 | ||
1923 | TH->GetEvent(track); // get next track | |
1924 | Int_t nhits = fHits->GetEntriesFast(); // get number of hits for this track | |
1925 | ||
1926 | if(nhits == 0) continue; // no hits in the TPC for this track | |
1927 | ||
1928 | //-------------------------------------------------------------- | |
1929 | // Loop over hits | |
1930 | //-------------------------------------------------------------- | |
1931 | ||
1932 | for(Int_t hit=0;hit<nhits;hit++){ | |
1933 | ||
1934 | tpcHit = (AliTPChit*)fHits->UncheckedAt(hit); // get a pointer to a hit | |
1935 | ||
1936 | Int_t sector=tpcHit->fSector; // sector number | |
1937 | if(sector != isec) continue; //terminate iteration | |
1938 | ||
1939 | currentTrack = tpcHit->fTrack; // track number | |
1940 | if(currentTrack != previousTrack){ | |
1941 | ||
1942 | // store already filled fTrack | |
1943 | ||
1944 | for(i=0;i<nrows;i++){ | |
1945 | if(previousTrack != -1){ | |
1946 | if(n_of_electrons[i]>0){ | |
1947 | TVector &v = *tr[i]; | |
1948 | v(0) = previousTrack; | |
1949 | tr[i]->ResizeTo(3*n_of_electrons[i]+1); // shrink if necessary | |
1950 | row[i]->Add(tr[i]); | |
1951 | } | |
1952 | else{ | |
1953 | delete tr[i]; // delete empty TVector | |
1954 | tr[i]=0; | |
1955 | } | |
1956 | } | |
1957 | ||
1958 | n_of_electrons[i]=0; | |
1959 | tr[i] = new TVector(361); // TVectors for the next fTrack | |
1960 | ||
1961 | } // end of loop over rows | |
1962 | ||
1963 | previousTrack=currentTrack; // update track label | |
1964 | } | |
1965 | ||
1966 | Int_t QI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons) | |
1967 | ||
1968 | //--------------------------------------------------- | |
1969 | // Calculate the electron attachment probability | |
1970 | //--------------------------------------------------- | |
1971 | ||
1972 | Float_t time = 1.e6*(z_end-TMath::Abs(tpcHit->fZ))/fTPCParam->GetDriftV(); | |
1973 | // in microseconds! | |
1974 | Float_t AttProb = fTPCParam->GetAttCoef()* | |
1975 | fTPCParam->GetOxyCont()*time; // fraction! | |
1976 | ||
1977 | //----------------------------------------------- | |
1978 | // Loop over electrons | |
1979 | //----------------------------------------------- | |
1980 | ||
1981 | for(Int_t nel=0;nel<QI;nel++){ | |
1982 | // skip if electron lost due to the attachment | |
1983 | if((gRandom->Rndm(0)) < AttProb) continue; // electron lost! | |
1984 | xyz[0]=tpcHit->fX; | |
1985 | xyz[1]=tpcHit->fY; | |
1986 | xyz[2]=tpcHit->fZ; | |
1987 | ElDiff(xyz); // Appply the diffusion | |
1988 | Int_t row_number; | |
1989 | fTPCParam->XYZtoCRXYZ(xyz,isec,row_number,3); | |
1990 | ||
1991 | //transform position to local coordinates | |
1992 | //option 3 means that we calculate x position relative to | |
1993 | //nearest pad row | |
1994 | ||
1995 | if ((row_number<0)||row_number>=fTPCParam->GetNRow(isec)) continue; | |
1996 | ExB(xyz); // ExB effect at the sense wires | |
1997 | n_of_electrons[row_number]++; | |
1998 | //---------------------------------- | |
1999 | // Expand vector if necessary | |
2000 | //---------------------------------- | |
fe4da5cc | 2001 | if(n_of_electrons[row_number]>120){ |
2002 | Int_t range = tr[row_number]->GetNrows(); | |
8c555625 | 2003 | if(n_of_electrons[row_number] > (range-1)/3){ |
2004 | tr[row_number]->ResizeTo(range+150); // Add 50 electrons | |
fe4da5cc | 2005 | } |
2006 | } | |
2007 | ||
fe4da5cc | 2008 | TVector &v = *tr[row_number]; |
8c555625 | 2009 | Int_t idx = 3*n_of_electrons[row_number]-2; |
2010 | ||
2011 | v(idx)= xyz[0]; // X | |
2012 | v(idx+1)=xyz[1]; // Y (along the pad-row) | |
2013 | v(idx+2)=xyz[2]; // Z | |
2014 | ||
2015 | } // end of loop over electrons | |
2016 | ||
2017 | } // end of loop over hits | |
2018 | } // end of loop over tracks | |
2019 | ||
2020 | // | |
2021 | // store remaining track (the last one) if not empty | |
2022 | // | |
2023 | ||
2024 | for(i=0;i<nrows;i++){ | |
2025 | if(n_of_electrons[i]>0){ | |
2026 | TVector &v = *tr[i]; | |
2027 | v(0) = previousTrack; | |
2028 | tr[i]->ResizeTo(3*n_of_electrons[i]+1); // shrink if necessary | |
2029 | row[i]->Add(tr[i]); | |
fe4da5cc | 2030 | } |
2031 | else{ | |
8c555625 | 2032 | delete tr[i]; |
2033 | tr[i]=0; | |
2034 | } | |
2035 | } | |
2036 | ||
2037 | delete [] tr; | |
2038 | delete [] n_of_electrons; | |
2039 | ||
2040 | } // end of MakeSector | |
2041 | //_____________________________________________________________________________ | |
2042 | void AliTPC::GetCrossTalk (Int_t iFlag,TObjArray *p,Int_t ntracks,Int_t *npads, | |
2043 | TMatrix *m) | |
2044 | { | |
2045 | ||
2046 | //------------------------------------------------------------- | |
2047 | // Calculates the cross-talk from one row (inner or outer) | |
2048 | //------------------------------------------------------------- | |
2049 | ||
2050 | //----------------------------------------------------------------- | |
2051 | // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl | |
2052 | //----------------------------------------------------------------- | |
2053 | ||
2054 | // | |
2055 | // iFlag=2 & 3 --> cross-talk from the inner row | |
2056 | // iFlag=0 & 4 --> cross-talk from the outer row | |
2057 | // | |
2058 | AliTPCParam * fTPCParam = &(fDigParam->GetParam()); | |
2059 | AliTPCPRF2D * fPRF2D = &(fDigParam->GetPRF2D()); | |
2060 | AliTPCRF1D * fRF = &(fDigParam->GetRF()); | |
2061 | ||
2062 | //to make code faster | |
2063 | ||
2064 | Float_t zwidth = fTPCParam->GetZWidth(); | |
2065 | Float_t zwidthm1 =1/fTPCParam->GetZWidth(); | |
2066 | ||
2067 | Int_t PadNumber; | |
2068 | Float_t xwire; | |
2069 | ||
2070 | Int_t nPadsSignal; // for this pads the signal is calculated | |
2071 | Float_t range; // sense wire range | |
2072 | Int_t nPadsDiff; | |
2073 | ||
2074 | Float_t IneffFactor=0.5; // gain inefficiency close to the edges | |
2075 | ||
2076 | if(iFlag == 0){ | |
2077 | // 1-->0 | |
2078 | nPadsSignal = npads[1]; | |
2079 | range = ((npads[1]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth(); | |
2080 | nPadsDiff = (npads[1]-npads[0])/2; | |
2081 | } | |
2082 | else if (iFlag == 2){ | |
2083 | // 1-->2 | |
2084 | nPadsSignal = npads[2]; | |
2085 | range = ((npads[1]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth(); | |
2086 | nPadsDiff = 0; | |
2087 | } | |
2088 | else if (iFlag == 3){ | |
2089 | // 0-->1 | |
2090 | nPadsSignal = npads[1]; | |
2091 | range = ((npads[0]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth(); | |
2092 | nPadsDiff = 0; | |
2093 | } | |
2094 | else{ | |
2095 | // 2-->1 | |
2096 | nPadsSignal = npads[2]; | |
2097 | range = ((npads[2]-1)/2 + 0.5)*fTPCParam->GetPadPitchWidth(); | |
2098 | nPadsDiff = (npads[2]-npads[1])/2; | |
2099 | } | |
2100 | ||
2101 | range-=0.5; // dead zone close to the edges | |
2102 | ||
2103 | TVector *tv; | |
2104 | TMatrix &signal = *m; | |
2105 | ||
1283eee5 | 2106 | Int_t CentralPad = (nPadsSignal-1)/2; |
8c555625 | 2107 | Float_t PadSignal[7]; // signal from a single electron |
2108 | // Loop over tracks | |
2109 | for(Int_t nt=0;nt<ntracks;nt++){ | |
2110 | tv=(TVector*)p->At(nt); // pointer to a track | |
2111 | TVector &v = *tv; | |
2112 | Int_t nElectrons = (tv->GetNrows()-1)/4; | |
2113 | // Loop over electrons | |
2114 | for(Int_t nel=0; nel<nElectrons; nel++){ | |
2115 | Int_t idx=nel*4; | |
2116 | xwire=v(idx+1); | |
2117 | ||
2118 | if (iFlag==0) xwire+=fTPCParam->GetPadPitchLength(); | |
2119 | if (iFlag==2) xwire-=fTPCParam->GetPadPitchLength(); | |
2120 | if (iFlag==3) xwire-=fTPCParam->GetPadPitchLength(); | |
2121 | if (iFlag==4) xwire+=fTPCParam->GetPadPitchLength(); | |
2122 | ||
2123 | // electron acceptance for the cross-talk (only the closest wire) | |
2124 | ||
2125 | Float_t dxMax = fTPCParam->GetPadPitchLength()*0.5+fTPCParam->GetWWPitch(); | |
2126 | if(TMath::Abs(xwire)>dxMax) continue; | |
2127 | ||
2128 | Float_t y = v(idx+2); | |
2129 | Float_t z = v(idx+3); | |
2130 | Float_t absy=TMath::Abs(y); | |
2131 | ||
2132 | if(absy < 0.5*fTPCParam->GetPadPitchWidth()){ | |
2133 | PadNumber=CentralPad; | |
2134 | } | |
2135 | else if (absy < range){ | |
8022212f | 2136 | PadNumber=(Int_t) ((absy-0.5*fTPCParam->GetPadPitchWidth())/fTPCParam->GetPadPitchWidth() +1.); |
8c555625 | 2137 | PadNumber=(Int_t) (TMath::Sign((Float_t)PadNumber, y)+CentralPad); |
2138 | } | |
2139 | else continue; // electron out of sense wire range, lost at the sector edge | |
2140 | ||
2141 | Float_t aval = (absy<range-0.5) ? v(idx+4):v(idx+4)*IneffFactor; | |
2142 | ||
2143 | Float_t dist = y - (Float_t)(PadNumber-CentralPad)*fTPCParam->GetPadPitchWidth(); | |
2144 | ||
a03540ae | 2145 | for (Int_t i=0;i<7;i++){ |
b6738b36 | 2146 | PadSignal[i]=fPRF2D->GetPRF(-dist+(i-3)*fTPCParam->GetPadPitchWidth(),xwire)*aval; |
8c555625 | 2147 | |
a03540ae | 2148 | PadSignal[i] *= fTPCParam->GetPadCoupling(); |
2149 | } | |
8c555625 | 2150 | // real pad range |
2151 | ||
1283eee5 | 2152 | Int_t LeftPad = TMath::Max(0,PadNumber-3); |
2153 | Int_t RightPad = TMath::Min(nPadsSignal-1,PadNumber+3); | |
8c555625 | 2154 | |
2155 | Int_t pmin=LeftPad-PadNumber+3; // lower index of the pad_signal vector | |
2156 | Int_t pmax=RightPad-PadNumber+3; // upper index | |
2157 | ||
2158 | ||
1283eee5 | 2159 | Float_t z_drift = z*zwidthm1; |
8c555625 | 2160 | Float_t z_offset = z_drift-(Int_t)z_drift; |
2161 | //distance to the centre of nearest time bin (in time bin units) | |
1283eee5 | 2162 | Int_t FirstBucket = (Int_t)z_drift; |
8c555625 | 2163 | // MI check it --time offset |
2164 | for (Int_t i2=0;i2<4;i2++){ | |
2165 | Int_t TrueTime = FirstBucket+i2; // current time bucket | |
2166 | Float_t dz = (Float_t(i2)+z_offset)*zwidth; | |
2167 | Float_t ampl = fRF->GetRF(dz); | |
1283eee5 | 2168 | if((TrueTime>MAXTPCTBK-1)) break; // beyond the time range |
8c555625 | 2169 | |
2170 | ||
2171 | // loop over pads, from pmin to pmax | |
2172 | ||
2173 | for(Int_t i3=pmin;i3<pmax+1;i3++){ | |
2174 | Int_t TruePad = LeftPad+i3-pmin; | |
2175 | ||
1283eee5 | 2176 | if(TruePad<nPadsDiff || TruePad > nPadsSignal-nPadsDiff-1) continue; |
8c555625 | 2177 | |
2178 | TruePad -= nPadsDiff; | |
2179 | signal(TruePad,TrueTime)+=(PadSignal[i3]*ampl); // not converted to charge! | |
2180 | ||
2181 | } // end of loop over pads | |
2182 | } // end of loop over time bins | |
2183 | ||
2184 | } // end of loop over electrons | |
2185 | ||
2186 | } // end of loop over tracks | |
2187 | ||
2188 | } // end of GetCrossTalk | |
2189 | ||
2190 | ||
fe4da5cc | 2191 | |
2192 | //_____________________________________________________________________________ | |
2193 | void AliTPC::Init() | |
2194 | { | |
2195 | // | |
2196 | // Initialise TPC detector after definition of geometry | |
2197 | // | |
2198 | Int_t i; | |
2199 | // | |
2200 | printf("\n"); | |
2201 | for(i=0;i<35;i++) printf("*"); | |
2202 | printf(" TPC_INIT "); | |
2203 | for(i=0;i<35;i++) printf("*"); | |
2204 | printf("\n"); | |
2205 | // | |
2206 | for(i=0;i<80;i++) printf("*"); | |
2207 | printf("\n"); | |
2208 | } | |
2209 | ||
2210 | //_____________________________________________________________________________ | |
2211 | void AliTPC::MakeBranch(Option_t* option) | |
2212 | { | |
2213 | // | |
2214 | // Create Tree branches for the TPC. | |
2215 | // | |
2216 | Int_t buffersize = 4000; | |
2217 | char branchname[10]; | |
2218 | sprintf(branchname,"%s",GetName()); | |
2219 | ||
2220 | AliDetector::MakeBranch(option); | |
2221 | ||
2222 | char *D = strstr(option,"D"); | |
2223 | ||
2224 | if (fDigits && gAlice->TreeD() && D) { | |
2225 | gAlice->TreeD()->Branch(branchname,&fDigits, buffersize); | |
2226 | printf("Making Branch %s for digits\n",branchname); | |
2227 | } | |
2228 | ||
2229 | char *R = strstr(option,"R"); | |
2230 | ||
2231 | if (fClusters && gAlice->TreeR() && R) { | |
2232 | gAlice->TreeR()->Branch(branchname,&fClusters, buffersize); | |
2233 | printf("Making Branch %s for Clusters\n",branchname); | |
2234 | } | |
2235 | } | |
2236 | ||
2237 | //_____________________________________________________________________________ | |
2238 | void AliTPC::ResetDigits() | |
2239 | { | |
2240 | // | |
2241 | // Reset number of digits and the digits array for this detector | |
2242 | // reset clusters | |
2243 | // | |
2244 | fNdigits = 0; | |
8c555625 | 2245 | // if (fDigits) fDigits->Clear(); |
2246 | if (fDigParam->GetArray()!=0) fDigParam->GetArray()->Clear(); | |
fe4da5cc | 2247 | fNclusters = 0; |
2248 | if (fClusters) fClusters->Clear(); | |
2249 | } | |
2250 | ||
2251 | //_____________________________________________________________________________ | |
2252 | void AliTPC::SetSecAL(Int_t sec) | |
2253 | { | |
8c555625 | 2254 | //--------------------------------------------------- |
fe4da5cc | 2255 | // Activate/deactivate selection for lower sectors |
8c555625 | 2256 | //--------------------------------------------------- |
2257 | ||
2258 | //----------------------------------------------------------------- | |
2259 | // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl | |
2260 | //----------------------------------------------------------------- | |
2261 | ||
fe4da5cc | 2262 | fSecAL = sec; |
2263 | } | |
2264 | ||
2265 | //_____________________________________________________________________________ | |
2266 | void AliTPC::SetSecAU(Int_t sec) | |
2267 | { | |
8c555625 | 2268 | //---------------------------------------------------- |
fe4da5cc | 2269 | // Activate/deactivate selection for upper sectors |
8c555625 | 2270 | //--------------------------------------------------- |
2271 | ||
2272 | //----------------------------------------------------------------- | |
2273 | // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl | |
2274 | //----------------------------------------------------------------- | |
2275 | ||
fe4da5cc | 2276 | fSecAU = sec; |
2277 | } | |
2278 | ||
2279 | //_____________________________________________________________________________ | |
2280 | void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6) | |
2281 | { | |
8c555625 | 2282 | //---------------------------------------- |
fe4da5cc | 2283 | // Select active lower sectors |
8c555625 | 2284 | //---------------------------------------- |
2285 | ||
2286 | //----------------------------------------------------------------- | |
2287 | // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl | |
2288 | //----------------------------------------------------------------- | |
2289 | ||
fe4da5cc | 2290 | fSecLows[0] = s1; |
2291 | fSecLows[1] = s2; | |
2292 | fSecLows[2] = s3; | |
2293 | fSecLows[3] = s4; | |
2294 | fSecLows[4] = s5; | |
2295 | fSecLows[5] = s6; | |
2296 | } | |
2297 | ||
2298 | //_____________________________________________________________________________ | |
2299 | void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6, | |
2300 | Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10, | |
2301 | Int_t s11 , Int_t s12) | |
2302 | { | |
8c555625 | 2303 | //-------------------------------- |
fe4da5cc | 2304 | // Select active upper sectors |
8c555625 | 2305 | //-------------------------------- |
2306 | ||
2307 | //----------------------------------------------------------------- | |
2308 | // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl | |
2309 | //----------------------------------------------------------------- | |
2310 | ||
fe4da5cc | 2311 | fSecUps[0] = s1; |
2312 | fSecUps[1] = s2; | |
2313 | fSecUps[2] = s3; | |
2314 | fSecUps[3] = s4; | |
2315 | fSecUps[4] = s5; | |
2316 | fSecUps[5] = s6; | |
2317 | fSecUps[6] = s7; | |
2318 | fSecUps[7] = s8; | |
2319 | fSecUps[8] = s9; | |
2320 | fSecUps[9] = s10; | |
2321 | fSecUps[10] = s11; | |
2322 | fSecUps[11] = s12; | |
2323 | } | |
2324 | ||
2325 | //_____________________________________________________________________________ | |
2326 | void AliTPC::SetSens(Int_t sens) | |
2327 | { | |
8c555625 | 2328 | |
2329 | //------------------------------------------------------------- | |
2330 | // Activates/deactivates the sensitive strips at the center of | |
2331 | // the pad row -- this is for the space-point resolution calculations | |
2332 | //------------------------------------------------------------- | |
2333 | ||
2334 | //----------------------------------------------------------------- | |
2335 | // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl | |
2336 | //----------------------------------------------------------------- | |
2337 | ||
fe4da5cc | 2338 | fSens = sens; |
2339 | } | |
4b0fdcad | 2340 | |
2341 | void AliTPC::SetSide(Float_t side) | |
2342 | { | |
2343 | fSide = side; | |
2344 | ||
2345 | } | |
1283eee5 | 2346 | //____________________________________________________________________________ |
2347 | void AliTPC::SetGasMixt(Int_t nc,Int_t c1,Int_t c2,Int_t c3,Float_t p1, | |
2348 | Float_t p2,Float_t p3) | |
2349 | { | |
fe4da5cc | 2350 | |
1283eee5 | 2351 | fNoComp = nc; |
2352 | ||
2353 | fMixtComp[0]=c1; | |
2354 | fMixtComp[1]=c2; | |
2355 | fMixtComp[2]=c3; | |
2356 | ||
2357 | fMixtProp[0]=p1; | |
2358 | fMixtProp[1]=p2; | |
2359 | fMixtProp[2]=p3; | |
2360 | ||
2361 | ||
2362 | } | |
fe4da5cc | 2363 | //_____________________________________________________________________________ |
2364 | void AliTPC::Streamer(TBuffer &R__b) | |
2365 | { | |
2366 | // | |
2367 | // Stream an object of class AliTPC. | |
2368 | // | |
2369 | if (R__b.IsReading()) { | |
2370 | Version_t R__v = R__b.ReadVersion(); if (R__v) { } | |
2371 | AliDetector::Streamer(R__b); | |
2372 | if (R__v < 2) return; | |
2373 | R__b >> fNsectors; | |
2374 | R__b >> fNclusters; | |
2375 | R__b >> fNtracks; | |
2376 | fClustersIndex = new Int_t[fNsectors+1]; | |
2377 | fDigitsIndex = new Int_t[fNsectors+1]; | |
2378 | } else { | |
2379 | R__b.WriteVersion(AliTPC::IsA()); | |
2380 | AliDetector::Streamer(R__b); | |
2381 | R__b << fNsectors; | |
2382 | R__b << fNclusters; | |
2383 | R__b << fNtracks; | |
2384 | } | |
2385 | } | |
2386 | ||
2387 | ClassImp(AliTPCcluster) | |
2388 | ||
2389 | //_____________________________________________________________________________ | |
2390 | AliTPCcluster::AliTPCcluster(Float_t *hits, Int_t *lab) | |
2391 | { | |
2392 | // | |
2393 | // Creates a simulated cluster for the TPC | |
2394 | // | |
2395 | fTracks[0] = lab[0]; | |
2396 | fTracks[1] = lab[1]; | |
2397 | fTracks[2] = lab[2]; | |
2398 | fSector = lab[3]; | |
2399 | fPadRow = lab[4]; | |
2400 | fY = hits[0]; | |
2401 | fZ = hits[1]; | |
2402 | fQ = hits[2]; | |
2403 | fSigmaY2 = hits[3]; | |
2404 | fSigmaZ2 = hits[4]; | |
2405 | } | |
2406 | ||
2407 | //_____________________________________________________________________________ | |
8c555625 | 2408 | void AliTPCcluster::GetXYZ(Float_t *x, const AliTPCParam *par) const |
fe4da5cc | 2409 | { |
2410 | // | |
8c555625 | 2411 | // Transformation from local to global coordinate system |
fe4da5cc | 2412 | // |
8c555625 | 2413 | x[0]=par->GetPadRowRadii(fSector,fPadRow-1); |
2414 | x[1]=fY; | |
2415 | x[2]=fZ; | |
2416 | par->CRXYZtoXYZ(x,fSector,fPadRow-1,1); | |
fe4da5cc | 2417 | } |
2418 | ||
8c555625 | 2419 | //_____________________________________________________________________________ |
2420 | Int_t AliTPCcluster::Compare(TObject * o) | |
2421 | { | |
2422 | // | |
2423 | // compare two clusters according y coordinata | |
2424 | AliTPCcluster *cl= (AliTPCcluster *)o; | |
2425 | if (fY<cl->fY) return -1; | |
2426 | if (fY==cl->fY) return 0; | |
2427 | return 1; | |
2428 | } | |
2429 | ||
2430 | Bool_t AliTPCcluster::IsSortable() const | |
2431 | { | |
2432 | // | |
2433 | //make AliTPCcluster sortabale | |
2434 | return kTRUE; | |
2435 | } | |
2436 | ||
2437 | ||
2438 | ||
fe4da5cc | 2439 | ClassImp(AliTPCdigit) |
2440 | ||
2441 | //_____________________________________________________________________________ | |
2442 | AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits): | |
2443 | AliDigit(tracks) | |
2444 | { | |
2445 | // | |
2446 | // Creates a TPC digit object | |
2447 | // | |
2448 | fSector = digits[0]; | |
2449 | fPadRow = digits[1]; | |
2450 | fPad = digits[2]; | |
2451 | fTime = digits[3]; | |
2452 | fSignal = digits[4]; | |
2453 | } | |
2454 | ||
2455 | ||
2456 | ClassImp(AliTPChit) | |
2457 | ||
2458 | //_____________________________________________________________________________ | |
2459 | AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits): | |
2460 | AliHit(shunt,track) | |
2461 | { | |
2462 | // | |
2463 | // Creates a TPC hit object | |
2464 | // | |
2465 | fSector = vol[0]; | |
2466 | fPadRow = vol[1]; | |
2467 | fX = hits[0]; | |
2468 | fY = hits[1]; | |
2469 | fZ = hits[2]; | |
2470 | fQ = hits[3]; | |
2471 | } | |
2472 | ||
2473 | ||
2474 | ClassImp(AliTPCtrack) | |
2475 | ||
2476 | //_____________________________________________________________________________ | |
2477 | AliTPCtrack::AliTPCtrack(Float_t *hits) | |
2478 | { | |
2479 | // | |
2480 | // Default creator for a TPC reconstructed track object | |
2481 | // | |
2482 | ref=hits[0]; // This is dummy code ! | |
2483 | } | |
2484 | ||
2485 | AliTPCtrack::AliTPCtrack(const AliTPCcluster& c,const TVector& xx, | |
8c555625 | 2486 | const TMatrix& CC, const AliTPCParam *p): |
fe4da5cc | 2487 | x(xx),C(CC),clusters(MAX_CLUSTER) |
2488 | { | |
2489 | // | |
2490 | // Standard creator for a TPC reconstructed track object | |
2491 | // | |
2492 | chi2=0.; | |
8c555625 | 2493 | int sec=c.fSector-1, row=c.fPadRow-1; |
2494 | ref=p->GetPadRowRadii(sec+1,row); | |
2495 | ||
fe4da5cc | 2496 | if (sec<24) { |
8c555625 | 2497 | fAlpha=(sec%12)*alpha_low; |
fe4da5cc | 2498 | } else { |
8c555625 | 2499 | fAlpha=(sec%24)*alpha_up; |
fe4da5cc | 2500 | } |
2501 | clusters.AddLast((AliTPCcluster*)(&c)); | |
2502 | } | |
2503 | ||
2504 | //_____________________________________________________________________________ | |
2505 | AliTPCtrack::AliTPCtrack(const AliTPCtrack& t) : x(t.x), C(t.C), | |
2506 | clusters(t.clusters.GetEntriesFast()) | |
2507 | { | |
2508 | // | |
2509 | // Copy creator for a TPC reconstructed track | |
2510 | // | |
2511 | ref=t.ref; | |
2512 | chi2=t.chi2; | |
2513 | fAlpha=t.fAlpha; | |
2514 | int n=t.clusters.GetEntriesFast(); | |
2515 | for (int i=0; i<n; i++) clusters.AddLast(t.clusters.UncheckedAt(i)); | |
2516 | } | |
2517 | ||
2518 | //_____________________________________________________________________________ | |
2519 | Double_t AliTPCtrack::GetY(Double_t xk) const | |
2520 | { | |
2521 | // | |
2522 | // | |
2523 | // | |
2524 | Double_t c2=x(2)*xk - x(3); | |
2525 | if (TMath::Abs(c2) >= 0.999) { | |
2526 | if (*this>10) cerr<<*this<<" AliTPCtrack warning: No y for this x !\n"; | |
2527 | return 0.; | |
2528 | } | |
2529 | Double_t c1=x(2)*ref - x(3); | |
2530 | Double_t r1=sqrt(1.-c1*c1), r2=sqrt(1.-c2*c2); | |
2531 | Double_t dx=xk-ref; | |
2532 | return x(0) + dx*(c1+c2)/(r1+r2); | |
2533 | } | |
2534 | ||
2535 | //_____________________________________________________________________________ | |
2536 | int AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm) | |
2537 | { | |
2538 | // | |
2539 | // Propagate a TPC reconstructed track | |
2540 | // | |
2541 | if (TMath::Abs(x(2)*xk - x(3)) >= 0.999) { | |
2542 | if (*this>3) cerr<<*this<<" AliTPCtrack warning: Propagation failed !\n"; | |
2543 | return 0; | |
2544 | } | |
8c555625 | 2545 | |
fe4da5cc | 2546 | Double_t x1=ref, x2=x1+0.5*(xk-x1), dx=x2-x1, y1=x(0), z1=x(1); |
2547 | Double_t c1=x(2)*x1 - x(3), r1=sqrt(1.- c1*c1); | |
2548 | Double_t c2=x(2)*x2 - x(3), r2=sqrt(1.- c2*c2); | |
2549 | ||
2550 | x(0) += dx*(c1+c2)/(r1+r2); | |
2551 | x(1) += dx*(c1+c2)/(c1*r2 + c2*r1)*x(4); | |
8c555625 | 2552 | |
fe4da5cc | 2553 | TMatrix F(5,5); F.UnitMatrix(); |
2554 | Double_t rr=r1+r2, cc=c1+c2, xx=x1+x2; | |
2555 | F(0,2)= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr); | |
2556 | F(0,3)=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr); | |
2557 | Double_t cr=c1*r2+c2*r1; | |
2558 | F(1,2)= dx*x(4)*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr); | |
2559 | F(1,3)=-dx*x(4)*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr); | |
2560 | F(1,4)= dx*cc/cr; | |
2561 | TMatrix tmp(F,TMatrix::kMult,C); | |
2562 | C.Mult(tmp,TMatrix(TMatrix::kTransposed,F)); | |
2563 | ||
2564 | ref=x2; | |
2565 | ||
2566 | //Multiple scattering****************** | |
2567 | Double_t ey=x(2)*ref - x(3); | |
2568 | Double_t ex=sqrt(1-ey*ey); | |
2569 | Double_t ez=x(4); | |
2570 | TMatrix Q(5,5); Q=0.; | |
2571 | Q(2,2)=ez*ez+ey*ey; Q(2,3)=-ex*ey; Q(2,4)=-ex*ez; | |
2572 | Q(3,2)=Q(2,3); Q(3,3)= ez*ez+ex*ex; Q(3,4)=-ey*ez; | |
2573 | Q(4,2)=Q(2,4); Q(4,3)= Q(3,4); Q(4,4)=1.; | |
2574 | ||
2575 | F=0; | |
2576 | F(2,2)=-x(2)*ex; F(2,3)=-x(2)*ey; | |
2577 | F(3,2)=-ex*(x(2)*ref-ey); F(3,3)=-(1.+ x(2)*ref*ey - ey*ey); | |
2578 | F(4,2)=-ez*ex; F(4,3)=-ez*ey; F(4,4)=1.; | |
2579 | ||
2580 | tmp.Mult(F,Q); | |
2581 | Q.Mult(tmp,TMatrix(TMatrix::kTransposed,F)); | |
2582 | ||
2583 | Double_t p2=GetPt()*GetPt()*(1.+x(4)*x(4)); | |
2584 | Double_t beta2=p2/(p2 + pm*pm); | |
2585 | Double_t d=sqrt((x1-ref)*(x1-ref)+(y1-x(0))*(y1-x(0))+(z1-x(1))*(z1-x(1))); | |
2586 | d*=2.; | |
2587 | Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho; | |
2588 | Q*=theta2; | |
2589 | C+=Q; | |
2590 | ||
2591 | //Energy losses************************ | |
2592 | Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d*rho; | |
2593 | if (x1 < x2) dE=-dE; | |
2594 | x(2)*=(1.- sqrt(p2+pm*pm)/p2*dE); | |
8c555625 | 2595 | //x(3)*=(1.- sqrt(p2+pm*pm)/p2*dE); |
fe4da5cc | 2596 | |
2597 | x1=ref; x2=xk; y1=x(0); z1=x(1); | |
2598 | c1=x(2)*x1 - x(3); r1=sqrt(1.- c1*c1); | |
2599 | c2=x(2)*x2 - x(3); r2=sqrt(1.- c2*c2); | |
2600 | ||
2601 | x(0) += dx*(c1+c2)/(r1+r2); | |
2602 | x(1) += dx*(c1+c2)/(c1*r2 + c2*r1)*x(4); | |
2603 | ||
2604 | F.UnitMatrix(); | |
2605 | rr=r1+r2; cc=c1+c2; xx=x1+x2; | |
2606 | F(0,2)= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr); | |
2607 | F(0,3)=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr); | |
2608 | cr=c1*r2+c2*r1; | |
2609 | F(1,2)= dx*x(4)*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr); | |
2610 | F(1,3)=-dx*x(4)*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr); | |
2611 | F(1,4)= dx*cc/cr; | |
2612 | tmp.Mult(F,C); | |
2613 | C.Mult(tmp,TMatrix(TMatrix::kTransposed,F)); | |
2614 | ||
2615 | ref=x2; | |
2616 | ||
2617 | return 1; | |
2618 | } | |
2619 | ||
2620 | //_____________________________________________________________________________ | |
2621 | void AliTPCtrack::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm) | |
2622 | { | |
2623 | // | |
2624 | // Propagate a reconstructed track from the vertex | |
2625 | // | |
2626 | Double_t c=x(2)*ref - x(3); | |
2627 | Double_t tgf=-x(3)/(x(2)*x(0) + sqrt(1-c*c)); | |
2628 | Double_t snf=tgf/sqrt(1.+ tgf*tgf); | |
2629 | Double_t xv=(x(3)+snf)/x(2); | |
2630 | PropagateTo(xv,x0,rho,pm); | |
2631 | } | |
2632 | ||
2633 | //_____________________________________________________________________________ | |
2634 | void AliTPCtrack::Update(const AliTPCcluster *c, Double_t chisq) | |
2635 | { | |
2636 | // | |
2637 | // Update statistics for a reconstructed TPC track | |
2638 | // | |
2639 | TMatrix H(2,5); H.UnitMatrix(); | |
2640 | TMatrix Ht(TMatrix::kTransposed,H); | |
2641 | TVector m(2); m(0)=c->fY; m(1)=c->fZ; | |
2642 | TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2; | |
2643 | ||
2644 | TMatrix tmp(H,TMatrix::kMult,C); | |
2645 | TMatrix R(tmp,TMatrix::kMult,Ht); R+=V; | |
2646 | ||
2647 | Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0); | |
2648 | R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1); | |
2649 | R(1,0)*=-1; R(0,1)=R(1,0); | |
2650 | R*=1./det; | |
2651 | ||
2652 | //R.Invert(); | |
2653 | ||
2654 | TMatrix K(C,TMatrix::kMult,Ht); K*=R; | |
2655 | ||
2656 | TVector savex=x; | |
2657 | x*=H; x-=m; x*=-1; x*=K; x+=savex; | |
2658 | if (TMath::Abs(x(2)*ref-x(3)) >= 0.999) { | |
2659 | if (*this>3) cerr<<*this<<" AliTPCtrack warning: Filtering failed !\n"; | |
2660 | x=savex; | |
2661 | return; | |
2662 | } | |
2663 | ||
2664 | TMatrix saveC=C; | |
2665 | C.Mult(K,tmp); C-=saveC; C*=-1; | |
2666 | ||
2667 | clusters.AddLast((AliTPCcluster*)c); | |
2668 | chi2 += chisq; | |
2669 | } | |
2670 | ||
2671 | //_____________________________________________________________________________ | |
2672 | int AliTPCtrack::Rotate(Double_t alpha) | |
2673 | { | |
2674 | // | |
2675 | // Rotate a reconstructed TPC track | |
2676 | // | |
2677 | fAlpha += alpha; | |
2678 | ||
2679 | Double_t x1=ref, y1=x(0); | |
2680 | Double_t ca=cos(alpha), sa=sin(alpha); | |
2681 | Double_t r1=x(2)*ref - x(3); | |
2682 | ||
2683 | ref = x1*ca + y1*sa; | |
2684 | x(0)=-x1*sa + y1*ca; | |
2685 | x(3)=x(3)*ca + (x(2)*y1 + sqrt(1.- r1*r1))*sa; | |
2686 | ||
2687 | Double_t r2=x(2)*ref - x(3); | |
2688 | if (TMath::Abs(r2) >= 0.999) { | |
2689 | if (*this>3) cerr<<*this<<" AliTPCtrack warning: Rotation failed !\n"; | |
2690 | return 0; | |
2691 | } | |
2692 | ||
2693 | Double_t y0=x(0) + sqrt(1.- r2*r2)/x(2); | |
2694 | if ((x(0)-y0)*x(2) >= 0.) { | |
2695 | if (*this>3) cerr<<*this<<" AliTPCtrack warning: Rotation failed !!!\n"; | |
2696 | return 0; | |
2697 | } | |
2698 | ||
2699 | TMatrix F(5,5); F.UnitMatrix(); | |
2700 | F(0,0)=ca; | |
2701 | F(3,0)=x(2)*sa; | |
2702 | F(3,2)=(y1 - r1*x1/sqrt(1.- r1*r1))*sa; | |
2703 | F(3,3)= ca + sa*r1/sqrt(1.- r1*r1); | |
2704 | TMatrix tmp(F,TMatrix::kMult,C); | |
2705 | // Double_t dy2=C(0,0); | |
2706 | C.Mult(tmp,TMatrix(TMatrix::kTransposed,F)); | |
2707 | // C(0,0)+=dy2*sa*sa*r1*r1/(1.- r1*r1); | |
2708 | // C(1,1)+=dy2*sa*sa*x(4)*x(4)/(1.- r1*r1); | |
2709 | ||
2710 | return 1; | |
2711 | } | |
2712 | ||
2713 | //_____________________________________________________________________________ | |
2714 | void AliTPCtrack::UseClusters() const | |
2715 | { | |
2716 | // | |
2717 | // | |
2718 | // | |
2719 | int num_of_clusters=clusters.GetEntriesFast(); | |
2720 | for (int i=0; i<num_of_clusters; i++) { | |
2721 | //if (i<=14) continue; | |
2722 | AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(i); | |
2723 | c->Use(); | |
2724 | } | |
2725 | } | |
2726 | ||
2727 | //_____________________________________________________________________________ | |
2728 | Double_t AliTPCtrack::GetPredictedChi2(const AliTPCcluster *c) const | |
2729 | { | |
2730 | // | |
2731 | // Calculate chi2 for a reconstructed TPC track | |
2732 | // | |
2733 | TMatrix H(2,5); H.UnitMatrix(); | |
2734 | TVector m(2); m(0)=c->fY; m(1)=c->fZ; | |
2735 | TMatrix V(2,2); V(0,0)=c->fSigmaY2; V(0,1)=0.; V(1,0)=0.; V(1,1)=c->fSigmaZ2; | |
2736 | TVector res=x; res*=H; res-=m; //res*=-1; | |
2737 | TMatrix tmp(H,TMatrix::kMult,C); | |
2738 | TMatrix R(tmp,TMatrix::kMult,TMatrix(TMatrix::kTransposed,H)); R+=V; | |
2739 | ||
2740 | Double_t det=(Double_t)R(0,0)*R(1,1) - (Double_t)R(0,1)*R(1,0); | |
2741 | if (TMath::Abs(det) < 1.e-10) { | |
2742 | if (*this>3) cerr<<*this<<" AliTPCtrack warning: Singular matrix !\n"; | |
2743 | return 1e10; | |
2744 | } | |
2745 | R(0,1)=R(0,0); R(0,0)=R(1,1); R(1,1)=R(0,1); | |
2746 | R(1,0)*=-1; R(0,1)=R(1,0); | |
2747 | R*=1./det; | |
2748 | ||
2749 | //R.Invert(); | |
2750 | ||
2751 | TVector r=res; | |
2752 | res*=R; | |
2753 | return r*res; | |
2754 | } | |
2755 | ||
2756 | //_____________________________________________________________________________ | |
2757 | int AliTPCtrack::GetLab() const | |
2758 | { | |
2759 | // | |
2760 | // | |
2761 | // | |
7ee081c2 | 2762 | int lab=123456789; |
fe4da5cc | 2763 | struct { |
2764 | int lab; | |
2765 | int max; | |
2766 | } s[MAX_CLUSTER]={{0,0}}; | |
2767 | ||
2768 | int i; | |
2769 | int num_of_clusters=clusters.GetEntriesFast(); | |
2770 | for (i=0; i<num_of_clusters; i++) { | |
2771 | AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(i); | |
8c555625 | 2772 | lab=TMath::Abs(c->fTracks[0]); |
fe4da5cc | 2773 | int j; |
2774 | for (j=0; j<MAX_CLUSTER; j++) | |
2775 | if (s[j].lab==lab || s[j].max==0) break; | |
2776 | s[j].lab=lab; | |
2777 | s[j].max++; | |
2778 | } | |
2779 | ||
2780 | int max=0; | |
2781 | for (i=0; i<num_of_clusters; i++) | |
2782 | if (s[i].max>max) {max=s[i].max; lab=s[i].lab;} | |
8c555625 | 2783 | |
2784 | for (i=0; i<num_of_clusters; i++) { | |
2785 | AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(i); | |
2786 | if (TMath::Abs(c->fTracks[1]) == lab || | |
2787 | TMath::Abs(c->fTracks[2]) == lab ) max++; | |
2788 | } | |
fe4da5cc | 2789 | |
2790 | if (1.-float(max)/num_of_clusters > 0.10) return -lab; | |
2791 | ||
2792 | if (num_of_clusters < 6) return lab; | |
2793 | ||
2794 | max=0; | |
2795 | for (i=1; i<=6; i++) { | |
2796 | AliTPCcluster *c=(AliTPCcluster*)clusters.UncheckedAt(num_of_clusters-i); | |
8c555625 | 2797 | if (lab == TMath::Abs(c->fTracks[0]) || |
2798 | lab == TMath::Abs(c->fTracks[1]) || | |
2799 | lab == TMath::Abs(c->fTracks[2])) max++; | |
fe4da5cc | 2800 | } |
8c555625 | 2801 | if (max<3) return -lab; |
fe4da5cc | 2802 | |
2803 | return lab; | |
2804 | } | |
2805 | ||
2806 | //_____________________________________________________________________________ | |
2807 | void AliTPCtrack::GetPxPyPz(Double_t& px, Double_t& py, Double_t& pz) const | |
2808 | { | |
2809 | // | |
2810 | // Get reconstructed TPC track momentum | |
2811 | // | |
2812 | Double_t pt=0.3*FIELD/TMath::Abs(x(2))/100; // GeV/c | |
2813 | Double_t r=x(2)*ref-x(3); | |
2814 | Double_t y0=x(0) + sqrt(1.- r*r)/x(2); | |
2815 | px=-pt*(x(0)-y0)*x(2); //cos(phi); | |
2816 | py=-pt*(x(3)-ref*x(2)); //sin(phi); | |
2817 | pz=pt*x(4); | |
2818 | Double_t tmp=px*TMath::Cos(fAlpha) - py*TMath::Sin(fAlpha); | |
2819 | py=px*TMath::Sin(fAlpha) + py*TMath::Cos(fAlpha); | |
2820 | px=tmp; | |
2821 | } | |
2822 | ||
2823 | //_____________________________________________________________________________ | |
2824 | // | |
2825 | // Classes for internal tracking use | |
2826 | // | |
2827 | ||
2828 | //_____________________________________________________________________________ | |
2829 | void AliTPCRow::InsertCluster(const AliTPCcluster* c) | |
2830 | { | |
2831 | // | |
2832 | // Insert a cluster in the list | |
2833 | // | |
2834 | if (num_of_clusters==MAX_CLUSTER_PER_ROW) { | |
2835 | cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return; | |
2836 | } | |
2837 | if (num_of_clusters==0) {clusters[num_of_clusters++]=c; return;} | |
2838 | int i=Find(c->fY); | |
2839 | memmove(clusters+i+1 ,clusters+i,(num_of_clusters-i)*sizeof(AliTPCcluster*)); | |
2840 | clusters[i]=c; num_of_clusters++; | |
2841 | } | |
2842 | ||
2843 | //_____________________________________________________________________________ | |
2844 | int AliTPCRow::Find(Double_t y) const | |
2845 | { | |
2846 | // | |
2847 | // | |
2848 | // | |
2849 | if (y <= clusters[0]->fY) return 0; | |
2850 | if (y > clusters[num_of_clusters-1]->fY) return num_of_clusters; | |
2851 | int b=0, e=num_of_clusters-1, m=(b+e)/2; | |
2852 | for (; b<e; m=(b+e)/2) { | |
2853 | if (y > clusters[m]->fY) b=m+1; | |
2854 | else e=m; | |
2855 | } | |
2856 | return m; | |
2857 | } | |
8c555625 | 2858 |