Add the number of local boards
[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 "AliLog.h"
29
30 #include "AliTRDcalibDB.h"
31 #include "AliTRDCommonParam.h"
32 #include "AliTRDpadPlane.h"
33 #include "AliTRDgeometry.h"
34 #include "AliTRDmcmTracklet.h"
35
36 ClassImp(AliTRDmcmTracklet)
37
38 //_____________________________________________________________________________
39 AliTRDmcmTracklet::AliTRDmcmTracklet() 
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)
54 {
55   //
56   // AliTRDmcmTracklet default constructor
57   //
58
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
70 }
71
72 //_____________________________________________________________________________
73 AliTRDmcmTracklet::AliTRDmcmTracklet(Int_t det, Int_t row, Int_t n) 
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)
88 {
89   //
90   // AliTRDmcmTracklet default constructor
91   //
92
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
104   fGPos = new TGraph(0);
105   fGAmp = new TGraph(0);
106
107 }
108
109 //_____________________________________________________________________________
110 AliTRDmcmTracklet::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   }
140
141 }
142
143 //_____________________________________________________________________________
144 AliTRDmcmTracklet::~AliTRDmcmTracklet() 
145 {
146   //
147   // AliTRDmcmTracklet destructor
148   //
149
150   if (fGPos != 0) delete fGPos;
151   if (fGAmp != 0) delete fGAmp;
152
153 }
154
155 //_____________________________________________________________________________
156 AliTRDmcmTracklet &AliTRDmcmTracklet::operator=(const AliTRDmcmTracklet &t)
157 {
158   //
159   // Assignment operator
160   //
161
162   if (this != &t) ((AliTRDmcmTracklet &) t).Copy(*this); 
163   return *this;
164
165 }
166
167 //_____________________________________________________________________________
168 void AliTRDmcmTracklet::Copy(TObject &t) const
169 {
170   //
171   // Copy function
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 //_____________________________________________________________________________
202 void AliTRDmcmTracklet::Reset() 
203 {
204   //
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 //_____________________________________________________________________________
239 void AliTRDmcmTracklet::AddCluster(Int_t icol, Int_t itb, Float_t *adc, Int_t *track) 
240 {
241   //
242   // Add a cluster to the tracklet
243   //
244  
245   if (fNclusters >= kNtimeBins) return;
246
247   for (Int_t icl = 0; icl < kNclsPads; icl++) {
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++;
258
259 }
260
261 //_____________________________________________________________________________
262 void AliTRDmcmTracklet::MakeTrackletGraph(AliTRDgeometry *geo, Float_t field) 
263 {
264   //
265   // Tracklet graph of positions (global coordinates, rotated [cm])
266   //
267   
268   if (!geo) {
269     AliError("No geometry.");
270     return;
271   }
272   
273   AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
274   if (!commonParam) {
275     AliError("No common parameters.");
276     return;
277   }
278
279   AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
280   if (!calibration) {
281     AliError("No instance of AliTRDcalibDB.");
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
292   Float_t samplFreq = commonParam->GetSamplingFrequency();
293
294   Int_t time, col;
295   Float_t amp[3];
296   Float_t xpos, ypos, xzero, colSize, timeBinSize;
297   Float_t vDrift, omegaTau, lorentzAngle, thetaSlope;
298   Float_t tiltingAngle;
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
313     ypos = GetClusY(amp,iplan);
314
315     colSize = padPlane->GetColSize(col);
316     vDrift = calibration->GetVdrift(fDetector,col,fRow);
317     timeBinSize = vDrift/samplFreq;
318     
319     // From v4-03-Release to HEAD28Mar06 the sign has changed from "-" to "+" 
320     // due to a change in the digitizer
321     omegaTau     = TMath::Sign(1.0,(Double_t)field)*GetOmegaTau(vDrift,TMath::Abs(field));
322     lorentzAngle = TMath::ATan(omegaTau)*180.0/TMath::Pi();
323     
324     xpos = (time+0.5) * timeBinSize;
325     xpos = geo->GetTime0(iplan) - xpos;
326
327     ypos = padPlane->GetColPos(col) - (ypos + 0.5) * colSize;
328
329     // ExB correction
330     xzero = geo->GetTime0(iplan);
331     ypos = ypos + (xpos-xzero) * omegaTau;
332
333     // tilted pads correction
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);
337
338     fGPos->SetPoint(npg,(Double_t)xpos,(Double_t)ypos);
339     npg++;
340
341   }
342
343   fGPos->Set(npg);
344   
345   fTime0 = geo->GetTime0(iplan) - AliTRDgeometry::CdrHght() - 0.5*AliTRDgeometry::CamHght();
346   fRowz = padPlane->GetRowPos(fRow) - padPlane->GetRowSize(fRow)/2.0;
347
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;
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   
364   Float_t fx = fTime0;
365   Float_t fy = fOffset;
366   
367   Float_t infSlope = TMath::ATan(fy/fx)/TMath::Pi()*180.0;    
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()));
370
371   fPt = 0.3 * field * 0.01 * r;
372   
373   return;
374
375 }
376
377 //_____________________________________________________________________________
378 void AliTRDmcmTracklet::MakeClusAmpGraph() 
379 {
380   //
381   // Tracklet graph of cluster charges
382   //
383
384   Int_t   time;
385   Float_t amp[3];
386   Int_t   npg = 0;
387
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 //_____________________________________________________________________________
413 Float_t AliTRDmcmTracklet::GetClusY(Float_t *adc, Int_t pla) const
414 {
415   //
416   // Cluster position in the phi direction in pad units (relative to the pad border)
417   //
418
419   Float_t ypos = 0.0;
420
421   Float_t a0 = adc[0];
422   Float_t a1 = adc[1];
423   Float_t a2 = adc[2];
424
425   Float_t t1, t2, w1 ,w2;
426
427   Float_t w = 1.0;  // pad units
428
429   Float_t sigma = 0.0;
430
431   switch(pla) {
432   case 0:
433     sigma = 0.515; break;
434   case 1:
435     sigma = 0.501; break;
436   case 2:
437     sigma = 0.491; break;
438   case 3:
439     sigma = 0.481; break;
440   case 4:
441     sigma = 0.471; break;
442   case 5:
443     sigma = 0.463; break;
444   default:
445     AliError("Wrong plane number.");
446     return 0.0;
447   }
448
449   sigma *= w;
450
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);
457   }
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);
464   }
465
466   ypos = w*(t1+t2)/(w1+w2);  // range: -0.5*w ... +0.5*w
467
468   return ypos;
469
470 }
471
472 //_____________________________________________________________________________
473 void 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) {
510           AliWarning("Too many tracks for this tracklet.");
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 //_____________________________________________________________________________
530 Float_t AliTRDmcmTracklet::GetOmegaTau(Float_t vdrift, Float_t field) const
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] 
553                  + p1[ib] * vdrift
554                  + p2[ib] * vdrift*vdrift
555                  + p3[ib] * vdrift*vdrift*vdrift;
556
557   return TMath::Tan(alphaL);
558
559 }