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 }