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