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