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