9ad1e0d4cb71f8d3687691b53d5ac28098d36af3
[u/mrichter/AliRoot.git] / TRD / AliTRDmcmTracklet.cxx
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 //                                                   //
18 //                                                   //
19 //  Tracklet object created in the local tracking    //
20 //                                                   //
21 //                                                   //
22 ///////////////////////////////////////////////////////
23
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
35 ClassImp(AliTRDmcmTracklet)
36
37 //_____________________________________________________________________________
38 AliTRDmcmTracklet::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 //_____________________________________________________________________________
76 AliTRDmcmTracklet::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 //_____________________________________________________________________________
116 AliTRDmcmTracklet::~AliTRDmcmTracklet() 
117 {
118
119   //
120   // AliTRDmcmTracklet destructor
121   //
122   delete fGPos;
123   delete fGAmp;
124 }
125
126 //_____________________________________________________________________________
127 AliTRDmcmTracklet &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 //_____________________________________________________________________________
139 void 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 //_____________________________________________________________________________
173 void 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++;
193
194 }
195
196 //_____________________________________________________________________________
197 void 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
230   Float_t samplFreq = calibration->GetSamplingFrequency();
231
232   Int_t time, col;
233   Float_t amp[3];
234   Float_t xpos, ypos, xzero, colSize, timeBinSize;
235   Float_t vDrift, omegaTau, lorentzAngle, thetaSlope;
236   Float_t tiltingAngle;
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
251     ypos = GetClusY(amp,iplan);
252
253     colSize = padPlane->GetColSize(col);
254     vDrift = calibration->GetVdrift(fDetector,col,fRow);
255     timeBinSize = vDrift/samplFreq;
256     
257     // From v4-03-Release to HEAD28Mar06 the sign has changed from "-" to "+" 
258     // due to a change in the digitizer
259     omegaTau = +TMath::Sign(1.0,(Double_t)field)*GetOmegaTau(vDrift,TMath::Abs(field));
260     lorentzAngle = TMath::ATan(omegaTau)*180.0/TMath::Pi();
261     
262     xpos = (time+0.5) * timeBinSize;
263     xpos = geo->GetTime0(iplan) - xpos;
264
265     ypos = padPlane->GetColPos(col) - (ypos + 0.5) * colSize;
266
267     // ExB correction
268     xzero = geo->GetTime0(iplan);
269     ypos = ypos + (xpos-xzero) * omegaTau;
270
271     // tilted pads correction
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);
275
276     fGPos->SetPoint(npg,(Double_t)xpos,(Double_t)ypos);
277     npg++;
278
279   }
280
281   fGPos->Set(npg);
282   
283   fTime0 = geo->GetTime0(iplan) - AliTRDgeometry::CdrHght() - 0.5*AliTRDgeometry::CamHght();
284   fRowz = padPlane->GetRowPos(fRow) - padPlane->GetRowSize(fRow)/2.0;
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   
298   Float_t fx = fTime0;
299   Float_t fy = fOffset;
300   
301   Float_t infSlope = TMath::ATan(fy/fx)/TMath::Pi()*180.0;    
302   Float_t alpha = fSlope - infSlope;
303   Float_t r = TMath::Sqrt(fx*fx + fy*fy)/(2.0*TMath::Sin(alpha/180.0*TMath::Pi()));
304
305   fPt = 0.3 * field * 0.01 * r;
306   
307   return;
308
309 }
310
311 //_____________________________________________________________________________
312 void 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 //_____________________________________________________________________________
353 Float_t AliTRDmcmTracklet::GetClusY(Float_t *adc, Int_t pla) const
354 {
355   //
356   // Cluster position in the phi direction in pad units (relative to the pad border)
357   //
358
359   Float_t ypos = 0.0;
360
361   Float_t a0 = adc[0];
362   Float_t a1 = adc[1];
363   Float_t a2 = adc[2];
364
365   Float_t t1, t2, w1 ,w2;
366
367   Float_t w = 1.0;  // pad units
368
369   Float_t sigma = 0.0;
370
371   switch(pla) {
372   case 0:
373     sigma = 0.515; break;
374   case 1:
375     sigma = 0.501; break;
376   case 2:
377     sigma = 0.491; break;
378   case 3:
379     sigma = 0.481; break;
380   case 4:
381     sigma = 0.471; break;
382   case 5:
383     sigma = 0.463; break;
384   default:
385     Error("GetClusY","Wrong plane number.");
386     return 0.0;
387   }
388
389   sigma *= w;
390
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);
397   }
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);
404   }
405
406   ypos = w*(t1+t2)/(w1+w2);  // range: -0.5*w ... +0.5*w
407
408   return ypos;
409
410 }
411
412 //_____________________________________________________________________________
413 void 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 //_____________________________________________________________________________
470 Float_t AliTRDmcmTracklet::GetOmegaTau(Float_t vdrift, Float_t field) const
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 }