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 }