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 = 0.5*(padPlane->GetRowPos(fRow) + padPlane->GetRowPos(fRow+1));
 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 }