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