Adding TRD trigger (B.Vulpescu)
[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 #include <TGraph.h>
17 #include <TMath.h>
18 #include <TF1.h>
19
20 #include "AliTRDcalibDB.h"
21 #include "AliTRDCommonParam.h"
22 #include "AliTRDpadPlane.h"
23 #include "AliTRDgeometry.h"
24
25 #include "AliTRDmcmTracklet.h"
26
27 ClassImp(AliTRDmcmTracklet)
28
29 //_____________________________________________________________________________
30 AliTRDmcmTracklet::AliTRDmcmTracklet() 
31 {
32
33   //
34   // AliTRDmcmTracklet default constructor
35   //
36
37   fDetector = -1;
38   fRow      = -1;
39
40   for (Int_t time = 0; time < kNtimeBins; time++) {
41     for (Int_t icl = 0; icl < kNclsPads; icl++) {
42       fADC[time][icl] = 0;
43     }
44     for (Int_t it = 0; it < kNdict; it++) {
45       fTrack[time][it] = -1;
46     }
47     fTime[time]   = 0;
48     fCol[time]    = 0;
49   }
50
51   fNclusters  =  0;
52   fN          =  0;
53   fTrackLabel = -1;
54
55   fGPos = 0;
56   fGAmp = 0;
57
58   fSlope  = 0.0;
59   fOffset = 0.0;
60   fTime0  = 0.0;
61   fRowz   = 0.0;
62   fPt     = 0.0;
63   fdQdl   = 0.0;
64
65 }
66
67 //_____________________________________________________________________________
68 AliTRDmcmTracklet::AliTRDmcmTracklet(Int_t det, Int_t row, Int_t n) 
69 {
70
71   //
72   // AliTRDmcmTracklet default constructor
73   //
74
75   fDetector = det;
76   fRow      = row;
77
78   for (Int_t time = 0; time < kNtimeBins; time++) {
79     for (Int_t icl = 0; icl < kNclsPads; icl++) {
80       fADC[time][icl] = 0;
81     }
82     for (Int_t it = 0; it < kNdict; it++) {
83       fTrack[time][it] = -1;
84     }
85     fTime[time]   = 0;
86     fCol[time]    = 0;
87   }
88
89   fNclusters = 0;
90
91   fN = n;
92
93   fTrackLabel = -1;
94
95   fGPos = new TGraph(0);
96   fGAmp = new TGraph(0);
97
98   fSlope  = 0.0;
99   fOffset = 0.0;
100   fTime0  = 0.0;
101   fRowz   = 0.0;
102   fPt     = 0.0;
103   fdQdl   = 0.0;
104
105 }
106
107 //_____________________________________________________________________________
108 AliTRDmcmTracklet::~AliTRDmcmTracklet() 
109 {
110
111   //
112   // AliTRDmcmTracklet destructor
113   //
114
115 }
116
117 //_____________________________________________________________________________
118 void AliTRDmcmTracklet::MakeTrackletGraph(AliTRDgeometry *geo, Float_t field) 
119 {
120
121   //
122   // Tracklet graph of positions (global coordinates, rotated [cm])
123   //
124   
125   if (!geo) {
126     Error("MakeTrackletGraph","No geometry.");
127     return;
128   }
129   
130   AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
131   if (!commonParam)
132   {
133     Error("MakeTrackletGraph","No common params.");
134     return;
135   }
136
137   AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
138   if (!calibration)
139   {
140     Error("MakeTrackletGraph","No instance of AliTRDcalibDB.");
141     return;
142   }
143
144   Int_t iplan, icham;
145
146   iplan = geo->GetPlane(fDetector);
147   icham = geo->GetChamber(fDetector);
148
149   AliTRDpadPlane *padPlane = commonParam->GetPadPlane(iplan,icham);
150
151   Float_t SamplFreq = calibration->GetSamplingFrequency();
152
153   Int_t time, col;
154   Float_t amp[3];
155   Float_t Xpos, Ypos, Xzero, ColSize, TimeBinSize;
156   Float_t vDrift, OmegaTau, LorentzAngle, ThetaSlope;
157   Float_t TiltingAngle;
158   Int_t npg = 0;
159
160   for (Int_t icl = 0; icl < fNclusters; icl++) {
161
162     time   = GetClusterTime(icl);
163
164     amp[0] = GetClusterADC(icl)[0];
165     amp[1] = GetClusterADC(icl)[1];
166     amp[2] = GetClusterADC(icl)[2];
167
168     col    = GetClusterCol(icl);
169
170     if (amp[0] < 0.0 || amp[1] < 0.0 || amp[2] < 0.0) continue;
171
172     Ypos = GetClusY(amp,iplan);
173
174     ColSize = padPlane->GetColSize(col);
175     vDrift = calibration->GetVdrift(fDetector,col,fRow);
176     TimeBinSize = vDrift/SamplFreq;
177     
178     // From v4-03-Release to HEAD28Mar06 the sign has changed from "-" to "+" 
179     // due to a change in the digitizer
180     OmegaTau = +TMath::Sign(1.0,(Double_t)field)*GetOmegaTau(vDrift,TMath::Abs(field));
181     LorentzAngle = TMath::ATan(OmegaTau)*180.0/TMath::Pi();
182     
183     Xpos = (time+0.5) * TimeBinSize;
184     Xpos = geo->GetTime0(iplan) - Xpos;
185
186     Ypos = padPlane->GetColPos(col) - (Ypos + 0.5) * ColSize;
187
188     // ExB correction
189     Xzero = geo->GetTime0(iplan);
190     Ypos = Ypos + (Xpos-Xzero) * OmegaTau;
191
192     // tilted pads correction
193     ThetaSlope = - padPlane->GetRowPos(fRow)/geo->GetTime0(iplan);
194     TiltingAngle = padPlane->GetTiltingAngle()/180.0*TMath::Pi();
195     Ypos = Ypos - (Xpos-Xzero) * ThetaSlope * TMath::Sin(TiltingAngle);
196
197     fGPos->SetPoint(npg,(Double_t)Xpos,(Double_t)Ypos);
198     npg++;
199
200   }
201
202   fGPos->Set(npg);
203   
204   fTime0 = geo->GetTime0(iplan) - AliTRDgeometry::CdrHght() - 0.5*AliTRDgeometry::CamHght();
205   fRowz = 0.5*(padPlane->GetRowPos(fRow) + padPlane->GetRowPos(fRow+1));
206
207   Double_t xMin = 0, xMax = 0, x, y;
208   fGPos->GetPoint(0    ,x,y); xMax = x + 0.1;
209   fGPos->GetPoint(npg-1,x,y); xMin = x - 0.1;
210   
211   TF1 *line = new TF1("line","[0]+x*[1]",xMin,xMax);
212   fGPos->Fit(line,"WRQ0");
213
214   fOffset = line->Eval(fTime0);
215   fSlope  = TMath::ATan(line->GetParameter(1))*180.0/TMath::Pi();
216
217   line->Delete();
218   
219   Float_t fX = fTime0;
220   Float_t fY = fOffset;
221   
222   Float_t infSlope = TMath::ATan(fY/fX)/TMath::Pi()*180.0;    
223   Float_t alpha = fSlope - infSlope;
224   Float_t R = TMath::Sqrt(fX*fX + fY*fY)/(2.0*TMath::Sin(alpha/180.0*TMath::Pi()));
225
226   fPt = 0.3 * field * 0.01 * R;
227   
228   return;
229
230 }
231
232 //_____________________________________________________________________________
233 void AliTRDmcmTracklet::MakeClusAmpGraph() 
234 {
235   //
236   // Tracklet graph of cluster charges
237   //
238
239   AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
240   if (!calibration)
241   {
242     Error("MakeClusAmpGraph","No instance of AliTRDcalibDB.");
243     return;
244   }
245
246   Int_t time;
247   Float_t amp[3];
248   Int_t npg = 0;
249   fdQdl = 0.0;
250   for (Int_t icl = 0; icl < fNclusters; icl++) {
251
252     time   = GetClusterTime(icl);
253
254     amp[0] = GetClusterADC(icl)[0];
255     amp[1] = GetClusterADC(icl)[1];
256     amp[2] = GetClusterADC(icl)[2];
257
258     fGAmp->SetPoint(npg,(Double_t)(time+0.5),(Double_t)(amp[0]+amp[1]+amp[2]));
259     npg++;
260
261     fdQdl += amp[0]+amp[1]+amp[2];
262
263   }
264
265   fGAmp->Set(npg);
266
267   fdQdl /= (Float_t)npg;
268
269   return;
270
271 }
272
273 //_____________________________________________________________________________
274 Float_t AliTRDmcmTracklet::GetClusY(Float_t *adc, Int_t pla) 
275 {
276   //
277   // Cluster position in the phi direction in pad units (relative to the pad border)
278   //
279
280   Float_t Ypos = 0.0;
281
282   Float_t A0 = adc[0];
283   Float_t A1 = adc[1];
284   Float_t A2 = adc[2];
285
286   Float_t T1, T2, W1 ,W2;
287
288   Float_t W = 1.0;  // pad units
289
290   Float_t Sigma = 0.0;
291
292   switch(pla) {
293   case 0:
294     Sigma = 0.515; break;
295   case 1:
296     Sigma = 0.501; break;
297   case 2:
298     Sigma = 0.491; break;
299   case 3:
300     Sigma = 0.481; break;
301   case 4:
302     Sigma = 0.471; break;
303   case 5:
304     Sigma = 0.463; break;
305   default:
306     Error("GetClusY","Wrong plane number.");
307     return 0.0;
308   }
309
310   Sigma *= W;
311
312   T1 = 0.0;
313   W1 = 0.0;
314   if( A0 > 0 ) {
315     W1 = A0*A0;
316     //W1 = A0;
317     T1 = W1*((Sigma*Sigma)/W*TMath::Log(A1/A0)-0.5*W);
318   }
319   T2 = 0.0;
320   W2 = 0.0;
321   if( A2 > 0 ) {
322     W2 = A2*A2;
323     //W2 = A2;
324     T2 = W2*((Sigma*Sigma)/W*TMath::Log(A2/A1)+0.5*W);
325   }
326
327   Ypos = W*(T1+T2)/(W1+W2);  // range: -0.5*W ... +0.5*W
328
329   return Ypos;
330
331 }
332
333 //_____________________________________________________________________________
334 void AliTRDmcmTracklet::CookLabel(Float_t frac) 
335 {
336   //
337   // Cook the track label from cluster labels
338   //
339
340   const Int_t kMaxTracks = 10;
341   Int_t trackLabel[kMaxTracks];
342   Int_t trackCount[kMaxTracks];
343   for (Int_t it = 0; it < kMaxTracks; it++) {
344     trackLabel[it] = -1;
345     trackCount[it] = 0;
346   }
347
348   Bool_t counted;
349   Int_t label, nTracks = 0;
350   for (Int_t icl = 0; icl < fNclusters; icl++) {
351
352     for (Int_t id = 0; id < kNdict; id++) {
353
354       if (fTrack[icl][id] == -1) continue;
355
356       label = fTrack[icl][id];
357
358       counted = kFALSE;
359       for (Int_t it = 0; it < nTracks; it++) {
360         if (label == trackLabel[it]) {
361           trackCount[it]++;
362           counted = kTRUE;
363           break;
364         }
365       }
366       if (!counted) {
367         trackLabel[nTracks] = label;
368         trackCount[nTracks]++;
369         nTracks++;
370         if (nTracks == kMaxTracks) {
371           Warning("CookLabel","Too many tracks for this tracklet.");
372           nTracks--;
373           break;
374         }
375       }
376
377     }
378
379   }
380
381   for (Int_t it = 0; it < kMaxTracks; it++) {
382     if (trackCount[it] >= (Int_t)(frac*fNclusters)) {
383       fTrackLabel = trackLabel[it];
384       break;
385     }
386   }
387
388 }
389
390 //_____________________________________________________________________________
391 Float_t AliTRDmcmTracklet::GetOmegaTau(Float_t vdrift, Float_t field)
392 {
393   //
394   // Returns omega*tau (tan(Lorentz-angle)) for a given drift velocity <vd> 
395   // and a B-field <b> for Xe/CO2 (15%).
396   // The values are according to a GARFIELD simulation.
397   //
398
399   //
400   // Copy of the "AliTRDcalibDB" function, taking as argument the magnetic field too
401   //
402   
403   const Int_t kNb = 5;
404   Float_t p0[kNb] = {  0.004810,  0.007412,  0.010252,  0.013409,  0.016888 };
405   Float_t p1[kNb] = {  0.054875,  0.081534,  0.107333,  0.131983,  0.155455 };
406   Float_t p2[kNb] = { -0.008682, -0.012896, -0.016987, -0.020880, -0.024623 };
407   Float_t p3[kNb] = {  0.000155,  0.000238,  0.000330,  0.000428,  0.000541 };
408
409   Int_t ib = ((Int_t) (10 * (field - 0.15)));
410   ib       = TMath::Max(  0,ib);
411   ib       = TMath::Min(kNb,ib);
412
413   Float_t alphaL = p0[ib] 
414       + p1[ib] * vdrift
415       + p2[ib] * vdrift*vdrift
416       + p3[ib] * vdrift*vdrift*vdrift;
417
418   return TMath::Tan(alphaL);
419
420 }