Reduced memory consumption (Bogdan)
[u/mrichter/AliRoot.git] / TRD / AliTRDmcmTracklet.cxx
CommitLineData
0ee00e25 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
e3b2b5e5 16///////////////////////////////////////////////////////
17// //
18// //
19// Tracklet object created in the local tracking //
20// //
21// //
22///////////////////////////////////////////////////////
23
0ee00e25 24#include <TGraph.h>
25#include <TMath.h>
26#include <TF1.h>
27
28#include "AliTRDcalibDB.h"
29#include "AliTRDCommonParam.h"
30#include "AliTRDpadPlane.h"
31#include "AliTRDgeometry.h"
32
33#include "AliTRDmcmTracklet.h"
34
35ClassImp(AliTRDmcmTracklet)
36
37//_____________________________________________________________________________
38AliTRDmcmTracklet::AliTRDmcmTracklet()
39{
0ee00e25 40 //
41 // AliTRDmcmTracklet default constructor
42 //
43
44 fDetector = -1;
45 fRow = -1;
46
47 for (Int_t time = 0; time < kNtimeBins; time++) {
48 for (Int_t icl = 0; icl < kNclsPads; icl++) {
49 fADC[time][icl] = 0;
50 }
51 for (Int_t it = 0; it < kNdict; it++) {
52 fTrack[time][it] = -1;
53 }
54 fTime[time] = 0;
55 fCol[time] = 0;
56 }
57
58 fNclusters = 0;
59 fN = 0;
60 fTrackLabel = -1;
61
62 fGPos = 0;
63 fGAmp = 0;
64
65 fSlope = 0.0;
66 fOffset = 0.0;
67 fTime0 = 0.0;
68 fRowz = 0.0;
69 fPt = 0.0;
70 fdQdl = 0.0;
71
72}
73
74//_____________________________________________________________________________
75AliTRDmcmTracklet::AliTRDmcmTracklet(Int_t det, Int_t row, Int_t n)
76{
0ee00e25 77 //
78 // AliTRDmcmTracklet default constructor
79 //
80
81 fDetector = det;
82 fRow = row;
83
84 for (Int_t time = 0; time < kNtimeBins; time++) {
85 for (Int_t icl = 0; icl < kNclsPads; icl++) {
86 fADC[time][icl] = 0;
87 }
88 for (Int_t it = 0; it < kNdict; it++) {
89 fTrack[time][it] = -1;
90 }
91 fTime[time] = 0;
92 fCol[time] = 0;
93 }
94
95 fNclusters = 0;
96
97 fN = n;
98
99 fTrackLabel = -1;
100
101 fGPos = new TGraph(0);
102 fGAmp = new TGraph(0);
103
104 fSlope = 0.0;
105 fOffset = 0.0;
106 fTime0 = 0.0;
107 fRowz = 0.0;
108 fPt = 0.0;
109 fdQdl = 0.0;
110
111}
112
113//_____________________________________________________________________________
114AliTRDmcmTracklet::~AliTRDmcmTracklet()
115{
0ee00e25 116 //
117 // AliTRDmcmTracklet destructor
118 //
3fcb9908 119
120 if (fGPos != 0) delete fGPos;
121 if (fGAmp != 0) delete fGAmp;
122
e3b2b5e5 123}
124
125//_____________________________________________________________________________
126AliTRDmcmTracklet &AliTRDmcmTracklet::operator=(const AliTRDmcmTracklet &t)
127{
128 //
129 // assignment operator
130 //
131
132 if (this != &t) ((AliTRDmcmTracklet &) t).Copy(*this);
133 return *this;
134
135}
136
137//_____________________________________________________________________________
138void AliTRDmcmTracklet::Copy(TObject &t) const
139{
140 //
141 // copy function
142 //
143
144 ((AliTRDmcmTracklet &) t).fDetector = fDetector;
145 ((AliTRDmcmTracklet &) t).fRow = fRow;
146 ((AliTRDmcmTracklet &) t).fTrackLabel = fTrackLabel;
147 ((AliTRDmcmTracklet &) t).fNclusters = fNclusters;
148 ((AliTRDmcmTracklet &) t).fN = fN;
149 ((AliTRDmcmTracklet &) t).fGPos = NULL;
150 ((AliTRDmcmTracklet &) t).fGAmp = NULL;
151 ((AliTRDmcmTracklet &) t).fTime0 = fTime0;
152 ((AliTRDmcmTracklet &) t).fRowz = fRowz;
153 ((AliTRDmcmTracklet &) t).fSlope = fSlope;
154 ((AliTRDmcmTracklet &) t).fOffset = fOffset;
155 ((AliTRDmcmTracklet &) t).fPt = fPt;
156 ((AliTRDmcmTracklet &) t).fdQdl = fdQdl;
157
158 for (Int_t time = 0; time < kNtimeBins; time++) {
159 for (Int_t icl = 0; icl < kNclsPads; icl++) {
160 ((AliTRDmcmTracklet &) t).fADC[time][icl] = 0;
161 }
162 for (Int_t it = 0; it < kNdict; it++) {
163 ((AliTRDmcmTracklet &) t).fTrack[time][it] = -1;
164 }
165 ((AliTRDmcmTracklet &) t).fTime[time] = 0;
166 ((AliTRDmcmTracklet &) t).fCol[time] = 0;
167 }
168
169}
170
171//_____________________________________________________________________________
3fcb9908 172void AliTRDmcmTracklet::Reset()
e3b2b5e5 173{
174 //
3fcb9908 175 // Reset the tracklet information
176 //
177
178 fDetector = -1;
179 fRow = -1;
180
181 for (Int_t time = 0; time < kNtimeBins; time++) {
182 for (Int_t icl = 0; icl < kNclsPads; icl++) {
183 fADC[time][icl] = 0;
184 }
185 for (Int_t it = 0; it < kNdict; it++) {
186 fTrack[time][it] = -1;
187 }
188 fTime[time] = 0;
189 fCol[time] = 0;
190 }
191
192 fNclusters = 0;
193 fN = 0;
194 fTrackLabel = -1;
195
196 fGPos->Set(0);
197 fGAmp->Set(0);
198
199 fSlope = 0.0;
200 fOffset = 0.0;
201 fTime0 = 0.0;
202 fRowz = 0.0;
203 fPt = 0.0;
204 fdQdl = 0.0;
205
206}
207
208//_____________________________________________________________________________
209void AliTRDmcmTracklet::AddCluster(Int_t icol, Int_t itb, Float_t *adc, Int_t *track)
210{
e3b2b5e5 211 //
3fcb9908 212 // Add a cluster to the tracklet
e3b2b5e5 213 //
214
215 if (fNclusters >= kNtimeBins) return;
216
217 for (Int_t icl = 0; icl < kNclsPads; icl++) {
218 //fADC[fNclusters][icl] = (Int_t)adc[icl];
219 fADC[fNclusters][icl] = adc[icl];
220 }
221
222 fTrack[fNclusters][0] = track[0];
223 fTrack[fNclusters][1] = track[1];
224 fTrack[fNclusters][2] = track[2];
225 fTime[fNclusters] = itb;
226 fCol[fNclusters] = icol;
227
228 fNclusters++;
0ee00e25 229
230}
231
232//_____________________________________________________________________________
233void AliTRDmcmTracklet::MakeTrackletGraph(AliTRDgeometry *geo, Float_t field)
234{
0ee00e25 235 //
236 // Tracklet graph of positions (global coordinates, rotated [cm])
237 //
238
239 if (!geo) {
240 Error("MakeTrackletGraph","No geometry.");
241 return;
242 }
243
244 AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
245 if (!commonParam)
246 {
247 Error("MakeTrackletGraph","No common params.");
248 return;
249 }
250
251 AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
252 if (!calibration)
253 {
254 Error("MakeTrackletGraph","No instance of AliTRDcalibDB.");
255 return;
256 }
257
258 Int_t iplan, icham;
259
260 iplan = geo->GetPlane(fDetector);
261 icham = geo->GetChamber(fDetector);
262
263 AliTRDpadPlane *padPlane = commonParam->GetPadPlane(iplan,icham);
264
e3b2b5e5 265 Float_t samplFreq = calibration->GetSamplingFrequency();
0ee00e25 266
267 Int_t time, col;
268 Float_t amp[3];
e3b2b5e5 269 Float_t xpos, ypos, xzero, colSize, timeBinSize;
270 Float_t vDrift, omegaTau, lorentzAngle, thetaSlope;
271 Float_t tiltingAngle;
0ee00e25 272 Int_t npg = 0;
273
274 for (Int_t icl = 0; icl < fNclusters; icl++) {
275
276 time = GetClusterTime(icl);
277
278 amp[0] = GetClusterADC(icl)[0];
279 amp[1] = GetClusterADC(icl)[1];
280 amp[2] = GetClusterADC(icl)[2];
281
282 col = GetClusterCol(icl);
283
284 if (amp[0] < 0.0 || amp[1] < 0.0 || amp[2] < 0.0) continue;
285
e3b2b5e5 286 ypos = GetClusY(amp,iplan);
0ee00e25 287
e3b2b5e5 288 colSize = padPlane->GetColSize(col);
0ee00e25 289 vDrift = calibration->GetVdrift(fDetector,col,fRow);
e3b2b5e5 290 timeBinSize = vDrift/samplFreq;
0ee00e25 291
292 // From v4-03-Release to HEAD28Mar06 the sign has changed from "-" to "+"
293 // due to a change in the digitizer
e3b2b5e5 294 omegaTau = +TMath::Sign(1.0,(Double_t)field)*GetOmegaTau(vDrift,TMath::Abs(field));
295 lorentzAngle = TMath::ATan(omegaTau)*180.0/TMath::Pi();
0ee00e25 296
e3b2b5e5 297 xpos = (time+0.5) * timeBinSize;
298 xpos = geo->GetTime0(iplan) - xpos;
0ee00e25 299
e3b2b5e5 300 ypos = padPlane->GetColPos(col) - (ypos + 0.5) * colSize;
0ee00e25 301
302 // ExB correction
e3b2b5e5 303 xzero = geo->GetTime0(iplan);
304 ypos = ypos + (xpos-xzero) * omegaTau;
0ee00e25 305
306 // tilted pads correction
e3b2b5e5 307 thetaSlope = - padPlane->GetRowPos(fRow)/geo->GetTime0(iplan);
308 tiltingAngle = padPlane->GetTiltingAngle()/180.0*TMath::Pi();
309 ypos = ypos - (xpos-xzero) * thetaSlope * TMath::Sin(tiltingAngle);
0ee00e25 310
e3b2b5e5 311 fGPos->SetPoint(npg,(Double_t)xpos,(Double_t)ypos);
0ee00e25 312 npg++;
313
314 }
315
316 fGPos->Set(npg);
317
318 fTime0 = geo->GetTime0(iplan) - AliTRDgeometry::CdrHght() - 0.5*AliTRDgeometry::CamHght();
3fcb9908 319 fRowz = 0.5*(padPlane->GetRowPos(fRow) + padPlane->GetRowPos(fRow+1));
0ee00e25 320
321 Double_t xMin = 0, xMax = 0, x, y;
322 fGPos->GetPoint(0 ,x,y); xMax = x + 0.1;
323 fGPos->GetPoint(npg-1,x,y); xMin = x - 0.1;
324
325 TF1 *line = new TF1("line","[0]+x*[1]",xMin,xMax);
326 fGPos->Fit(line,"WRQ0");
327
328 fOffset = line->Eval(fTime0);
329 fSlope = TMath::ATan(line->GetParameter(1))*180.0/TMath::Pi();
330
331 line->Delete();
332
e3b2b5e5 333 Float_t fx = fTime0;
334 Float_t fy = fOffset;
0ee00e25 335
e3b2b5e5 336 Float_t infSlope = TMath::ATan(fy/fx)/TMath::Pi()*180.0;
0ee00e25 337 Float_t alpha = fSlope - infSlope;
e3b2b5e5 338 Float_t r = TMath::Sqrt(fx*fx + fy*fy)/(2.0*TMath::Sin(alpha/180.0*TMath::Pi()));
0ee00e25 339
e3b2b5e5 340 fPt = 0.3 * field * 0.01 * r;
0ee00e25 341
342 return;
343
344}
345
346//_____________________________________________________________________________
347void AliTRDmcmTracklet::MakeClusAmpGraph()
348{
349 //
350 // Tracklet graph of cluster charges
351 //
352
0ee00e25 353 Int_t time;
354 Float_t amp[3];
355 Int_t npg = 0;
356 fdQdl = 0.0;
357 for (Int_t icl = 0; icl < fNclusters; icl++) {
358
359 time = GetClusterTime(icl);
360
361 amp[0] = GetClusterADC(icl)[0];
362 amp[1] = GetClusterADC(icl)[1];
363 amp[2] = GetClusterADC(icl)[2];
364
365 fGAmp->SetPoint(npg,(Double_t)(time+0.5),(Double_t)(amp[0]+amp[1]+amp[2]));
366 npg++;
367
368 fdQdl += amp[0]+amp[1]+amp[2];
369
370 }
371
372 fGAmp->Set(npg);
373
374 fdQdl /= (Float_t)npg;
375
376 return;
377
378}
379
380//_____________________________________________________________________________
e3b2b5e5 381Float_t AliTRDmcmTracklet::GetClusY(Float_t *adc, Int_t pla) const
0ee00e25 382{
383 //
384 // Cluster position in the phi direction in pad units (relative to the pad border)
385 //
386
e3b2b5e5 387 Float_t ypos = 0.0;
0ee00e25 388
e3b2b5e5 389 Float_t a0 = adc[0];
390 Float_t a1 = adc[1];
391 Float_t a2 = adc[2];
0ee00e25 392
e3b2b5e5 393 Float_t t1, t2, w1 ,w2;
0ee00e25 394
e3b2b5e5 395 Float_t w = 1.0; // pad units
0ee00e25 396
e3b2b5e5 397 Float_t sigma = 0.0;
0ee00e25 398
399 switch(pla) {
400 case 0:
e3b2b5e5 401 sigma = 0.515; break;
0ee00e25 402 case 1:
e3b2b5e5 403 sigma = 0.501; break;
0ee00e25 404 case 2:
e3b2b5e5 405 sigma = 0.491; break;
0ee00e25 406 case 3:
e3b2b5e5 407 sigma = 0.481; break;
0ee00e25 408 case 4:
e3b2b5e5 409 sigma = 0.471; break;
0ee00e25 410 case 5:
e3b2b5e5 411 sigma = 0.463; break;
0ee00e25 412 default:
413 Error("GetClusY","Wrong plane number.");
414 return 0.0;
415 }
416
e3b2b5e5 417 sigma *= w;
0ee00e25 418
e3b2b5e5 419 t1 = 0.0;
420 w1 = 0.0;
421 if( a0 > 0 ) {
422 w1 = a0*a0;
423 //w1 = a0;
424 t1 = w1*((sigma*sigma)/w*TMath::Log(a1/a0)-0.5*w);
0ee00e25 425 }
e3b2b5e5 426 t2 = 0.0;
427 w2 = 0.0;
428 if( a2 > 0 ) {
429 w2 = a2*a2;
430 //w2 = a2;
431 t2 = w2*((sigma*sigma)/w*TMath::Log(a2/a1)+0.5*w);
0ee00e25 432 }
433
e3b2b5e5 434 ypos = w*(t1+t2)/(w1+w2); // range: -0.5*w ... +0.5*w
0ee00e25 435
e3b2b5e5 436 return ypos;
0ee00e25 437
438}
439
440//_____________________________________________________________________________
441void AliTRDmcmTracklet::CookLabel(Float_t frac)
442{
443 //
444 // Cook the track label from cluster labels
445 //
446
447 const Int_t kMaxTracks = 10;
448 Int_t trackLabel[kMaxTracks];
449 Int_t trackCount[kMaxTracks];
450 for (Int_t it = 0; it < kMaxTracks; it++) {
451 trackLabel[it] = -1;
452 trackCount[it] = 0;
453 }
454
455 Bool_t counted;
456 Int_t label, nTracks = 0;
457 for (Int_t icl = 0; icl < fNclusters; icl++) {
458
459 for (Int_t id = 0; id < kNdict; id++) {
460
461 if (fTrack[icl][id] == -1) continue;
462
463 label = fTrack[icl][id];
464
465 counted = kFALSE;
466 for (Int_t it = 0; it < nTracks; it++) {
467 if (label == trackLabel[it]) {
468 trackCount[it]++;
469 counted = kTRUE;
470 break;
471 }
472 }
473 if (!counted) {
474 trackLabel[nTracks] = label;
475 trackCount[nTracks]++;
476 nTracks++;
477 if (nTracks == kMaxTracks) {
478 Warning("CookLabel","Too many tracks for this tracklet.");
479 nTracks--;
480 break;
481 }
482 }
483
484 }
485
486 }
487
488 for (Int_t it = 0; it < kMaxTracks; it++) {
489 if (trackCount[it] >= (Int_t)(frac*fNclusters)) {
490 fTrackLabel = trackLabel[it];
491 break;
492 }
493 }
494
495}
496
497//_____________________________________________________________________________
e3b2b5e5 498Float_t AliTRDmcmTracklet::GetOmegaTau(Float_t vdrift, Float_t field) const
0ee00e25 499{
500 //
501 // Returns omega*tau (tan(Lorentz-angle)) for a given drift velocity <vd>
502 // and a B-field <b> for Xe/CO2 (15%).
503 // The values are according to a GARFIELD simulation.
504 //
505
506 //
507 // Copy of the "AliTRDcalibDB" function, taking as argument the magnetic field too
508 //
509
510 const Int_t kNb = 5;
511 Float_t p0[kNb] = { 0.004810, 0.007412, 0.010252, 0.013409, 0.016888 };
512 Float_t p1[kNb] = { 0.054875, 0.081534, 0.107333, 0.131983, 0.155455 };
513 Float_t p2[kNb] = { -0.008682, -0.012896, -0.016987, -0.020880, -0.024623 };
514 Float_t p3[kNb] = { 0.000155, 0.000238, 0.000330, 0.000428, 0.000541 };
515
516 Int_t ib = ((Int_t) (10 * (field - 0.15)));
517 ib = TMath::Max( 0,ib);
518 ib = TMath::Min(kNb,ib);
519
520 Float_t alphaL = p0[ib]
521 + p1[ib] * vdrift
522 + p2[ib] * vdrift*vdrift
523 + p3[ib] * vdrift*vdrift*vdrift;
524
525 return TMath::Tan(alphaL);
526
527}