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