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