| 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 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame