TRD/AliTRDclusterizerV1.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 $Log$
  18 Revision 1.9  2000/11/01 14:53:20  cblume
  19 Merge with TRD-develop
  20
  21 Revision 1.1.4.5  2000/10/15 23:40:01  cblume
  22 Remove AliTRDconst
  23
  24 Revision 1.1.4.4  2000/10/06 16:49:46  cblume
  25 Made Getters const
  26
  27 Revision 1.1.4.3  2000/10/04 16:34:58  cblume
  28 Replace include files by forward declarations
  29
  30 Revision 1.1.4.2  2000/09/22 14:49:49  cblume
  31 Adapted to tracking code
  32
  33 Revision 1.8  2000/10/02 21:28:19  fca
  34 Removal of useless dependecies via forward declarations
  35
  36 Revision 1.7  2000/06/27 13:08:50  cblume
  37 Changed to Copy(TObject &A) to appease the HP-compiler
  38
  39 Revision 1.6  2000/06/09 11:10:07  cblume
  40 Compiler warnings and coding conventions, next round
  41
  42 Revision 1.5  2000/06/08 18:32:58  cblume
  43 Make code compliant to coding conventions
  44
  45 Revision 1.4  2000/06/07 16:27:01  cblume
  46 Try to remove compiler warnings on Sun and HP
  47
  48 Revision 1.3  2000/05/08 16:17:27  cblume
  49 Merge TRD-develop
  50
  51 Revision 1.1.4.1  2000/05/08 15:09:01  cblume
  52 Introduce AliTRDdigitsManager
  53
  54 Revision 1.1  2000/02/28 18:58:54  cblume
  55 Add new TRD classes
  56
  57 */
  58
  59 ///////////////////////////////////////////////////////////////////////////////
  60 //                                                                           //
  61 // TRD cluster finder for the slow simulator.
  62 //                                                                           //
  63 ///////////////////////////////////////////////////////////////////////////////
  64
  65 #include <TF1.h>
  66 #include <TTree.h>
  67 #include <TH1.h>
  68
  69 #include "AliRun.h"
  70
  71 #include "AliTRD.h"
  72 #include "AliTRDclusterizerV1.h"
  73 #include "AliTRDmatrix.h"
  74 #include "AliTRDgeometry.h"
  75 #include "AliTRDdigitizer.h"
  76 #include "AliTRDdataArrayF.h"
  77 #include "AliTRDdataArrayI.h"
  78 #include "AliTRDdigitsManager.h"
  79
  80 ClassImp(AliTRDclusterizerV1)
  81
  82 //_____________________________________________________________________________
  83 AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
  84 {
  85   //
  86   // AliTRDclusterizerV1 default constructor
  87   //
  88
  89   fDigitsManager = NULL;
  90
  91   fClusMaxThresh = 0;
  92   fClusSigThresh = 0;
  93
  94 }
  95
  96 //_____________________________________________________________________________
  97 AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title)
  98                     :AliTRDclusterizer(name,title)
  99 {
 100   //
 101   // AliTRDclusterizerV1 default constructor
 102   //
 103
 104   fDigitsManager = new AliTRDdigitsManager();
 105
 106   Init();
 107
 108 }
 109
 110 //_____________________________________________________________________________
 111 AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
 112 {
 113   //
 114   // AliTRDclusterizerV1 copy constructor
 115   //
 116
 117   ((AliTRDclusterizerV1 &) c).Copy(*this);
 118
 119 }
 120
 121 //_____________________________________________________________________________
 122 AliTRDclusterizerV1::~AliTRDclusterizerV1()
 123 {
 124   //
 125   // AliTRDclusterizerV1 destructor
 126   //
 127
 128   if (fDigitsManager) {
 129     delete fDigitsManager;
 130   }
 131
 132 }
 133
 134 //_____________________________________________________________________________
 135 AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c)
 136 {
 137   //
 138   // Assignment operator
 139   //
 140
 141   if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this);
 142   return *this;
 143
 144 }
 145
 146 //_____________________________________________________________________________
 147 void AliTRDclusterizerV1::Copy(TObject &c)
 148 {
 149   //
 150   // Copy function
 151   //
 152
 153   ((AliTRDclusterizerV1 &) c).fClusMaxThresh = fClusMaxThresh;
 154   ((AliTRDclusterizerV1 &) c).fClusSigThresh = fClusSigThresh;
 155   ((AliTRDclusterizerV1 &) c).fDigitsManager = NULL;
 156
 157   AliTRDclusterizer::Copy(c);
 158
 159 }
 160
 161 //_____________________________________________________________________________
 162 void AliTRDclusterizerV1::Init()
 163 {
 164   //
 165   // Initializes the cluster finder
 166   //
 167
 168   // The default parameter for the clustering
 169   fClusMaxThresh = 3;
 170   fClusSigThresh = 1;
 171
 172 }
 173
 174 //_____________________________________________________________________________
 175 Bool_t AliTRDclusterizerV1::ReadDigits()
 176 {
 177   //
 178   // Reads the digits arrays from the input aliroot file
 179   //
 180
 181   if (!fInputFile) {
 182     printf("AliTRDclusterizerV1::ReadDigits -- ");
 183     printf("No input file open\n");
 184     return kFALSE;
 185   }
 186
 187   // Read in the digit arrays
 188   return (fDigitsManager->ReadDigits());
 189
 190 }
 191
 192 //_____________________________________________________________________________
 193 Bool_t AliTRDclusterizerV1::MakeClusters()
 194 {
 195   //
 196   // Generates the cluster.
 197   //
 198
 199   Int_t row, col, time;
 200
 201   if (fTRD->IsVersion() != 1) {
 202     printf("AliTRDclusterizerV1::MakeCluster -- ");
 203     printf("TRD must be version 1 (slow simulator).\n");
 204     return kFALSE;
 205   }
 206
 207   // Get the geometry
 208   AliTRDgeometry *geo = fTRD->GetGeometry();
 209
 210   printf("AliTRDclusterizerV1::MakeCluster -- ");
 211   printf("Start creating clusters.\n");
 212
 213   AliTRDdataArrayI *digits;
 214   AliTRDdataArrayI *track0;
 215   AliTRDdataArrayI *track1;
 216   AliTRDdataArrayI *track2;
 217
 218   // Threshold value for the maximum
 219   Int_t maxThresh = fClusMaxThresh;
 220   // Threshold value for the digit signal
 221   Int_t sigThresh = fClusSigThresh;
 222
 223   // Iteration limit for unfolding procedure
 224   const Float_t kEpsilon = 0.01;
 225
 226   const Int_t   kNclus   = 3;
 227   const Int_t   kNsig    = 5;
 228   const Int_t   kNtrack  = 3 * kNclus;
 229
 230   Float_t padSignal[kNsig];
 231   Float_t clusterSignal[kNclus];
 232   Float_t clusterPads[kNclus];
 233   Int_t   clusterDigit[kNclus];
 234   Int_t   clusterTracks[kNtrack];
 235
 236   Int_t chamBeg = 0;
 237   Int_t chamEnd = AliTRDgeometry::Ncham();
 238   if (fTRD->GetSensChamber()  >= 0) {
 239     chamBeg = fTRD->GetSensChamber();
 240     chamEnd = chamBeg + 1;
 241   }
 242   Int_t planBeg = 0;
 243   Int_t planEnd = AliTRDgeometry::Nplan();
 244   if (fTRD->GetSensPlane()    >= 0) {
 245     planBeg = fTRD->GetSensPlane();
 246     planEnd = planBeg + 1;
 247   }
 248   Int_t sectBeg = 0;
 249   Int_t sectEnd = AliTRDgeometry::Nsect();
 250
 251   // Start clustering in every chamber
 252   for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
 253     for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
 254       for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
 255
 256         if (fTRD->GetSensSector() >= 0) {
 257           Int_t sens1 = fTRD->GetSensSector();
 258           Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
 259           sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
 260                  * AliTRDgeometry::Nsect();
 261           if (sens1 < sens2) {
 262             if ((isect < sens1) || (isect >= sens2)) continue;
 263           }
 264           else {
 265             if ((isect < sens1) && (isect >= sens2)) continue;
 266           }
 267         }
 268
 269         Int_t idet = geo->GetDetector(iplan,icham,isect);
 270
 271         Int_t nClusters       = 0;
 272         Int_t nClustersUnfold = 0;
 273
 274         printf("AliTRDclusterizerV1::MakeCluster -- ");
 275         printf("Analyzing chamber %d, plane %d, sector %d.\n"
 276               ,icham,iplan,isect);
 277
 278         Int_t nRowMax     = geo->GetRowMax(iplan,icham,isect);
 279         Int_t nColMax     = geo->GetColMax(iplan);
 280         Int_t nTimeBefore = geo->GetTimeBefore();
 281         Int_t nTimeTotal  = geo->GetTimeTotal();
 282
 283         // Get the digits
 284         digits = fDigitsManager->GetDigits(idet);
 285         digits->Expand();
 286         track0 = fDigitsManager->GetDictionary(idet,0);
 287         track0->Expand();
 288         track1 = fDigitsManager->GetDictionary(idet,1);
 289         track1->Expand();
 290         track2 = fDigitsManager->GetDictionary(idet,2);
 291         track2->Expand();
 292
 293         // Loop through the chamber and find the maxima
 294         for ( row = 0;  row <  nRowMax;    row++) {
 295           for ( col = 2;  col <  nColMax;    col++) {
 296             for (time = 0; time < nTimeTotal; time++) {
 297
 298               Int_t signalL = digits->GetDataUnchecked(row,col  ,time);
 299               Int_t signalM = digits->GetDataUnchecked(row,col-1,time);
 300               Int_t signalR = digits->GetDataUnchecked(row,col-2,time);
 301
 302               // Look for the maximum
 303               if ((signalM >= maxThresh) &&
 304                   (signalL >= sigThresh) &&
 305                   (signalR >= sigThresh)) {
 306                 if ((signalL < signalM) && (signalR < signalM)) {
 307                   // Maximum found, mark the position by a negative signal
 308                   digits->SetDataUnchecked(row,col-1,time,-signalM);
 309                 }
 310               }
 311
 312             }
 313           }
 314         }
 315
 316         // Now check the maxima and calculate the cluster position
 317         for ( row = 0;  row <  nRowMax  ;  row++) {
 318           for ( col = 1;  col <  nColMax-1;  col++) {
 319             for (time = 0; time < nTimeTotal; time++) {
 320
 321               // Maximum found ?
 322               if (digits->GetDataUnchecked(row,col,time) < 0) {
 323
 324                 Int_t iPad;
 325                 for (iPad = 0; iPad < kNclus; iPad++) {
 326                   Int_t iPadCol = col - 1 + iPad;
 327                   clusterSignal[iPad]     = TMath::Abs(digits->GetDataUnchecked(row
 328                                                                                ,iPadCol
 329                                                                                ,time));
 330                   clusterDigit[iPad]      = digits->GetIndexUnchecked(row,iPadCol,time);
 331                   clusterTracks[3*iPad  ] = track0->GetDataUnchecked(row,iPadCol,time) - 1;
 332                   clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1;
 333                   clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1;
 334                 }
 335
 336                 // Neighbouring maximum on right side?
 337                 Int_t iType = 0;
 338                 if (col < nColMax-3) {
 339                   if (digits->GetDataUnchecked(row,col+2,time) < 0) {
 340                     for (iPad = 0; iPad < kNsig; iPad++) {
 341                       padSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
 342                                                                            ,col-1+iPad
 343                                                                            ,time));
 344                     }
 345                     // Unfold the two maxima and set the signal on
 346                     // the overlapping pad to the ratio
 347                     clusterSignal[2] *= Unfold(kEpsilon,padSignal);
 348                     iType = 1;
 349                   }
 350                 }
 351
 352                 Float_t clusterCharge = clusterSignal[0]
 353                                       + clusterSignal[1]
 354                                       + clusterSignal[2];
 355
 356                 // Calculate the position of the cluster by using the
 357                 // center of gravity method
 358                 clusterPads[0] = row + 0.5;
 359                 clusterPads[1] = col + 0.5
 360                                + (clusterSignal[2] - clusterSignal[0])
 361                                / clusterCharge;
 362                 // Take the shift of the additional time bins into account
 363                 clusterPads[2] = time - nTimeBefore + 0.5;
 364
 365                 if (fVerbose) {
 366                   printf("-----------------------------------------------------------\n");
 367                   printf("Create cluster no. %d\n",nClusters);
 368                   printf("Position: row = %f, col = %f, time = %f\n",clusterPads[0]
 369                                                                     ,clusterPads[1]
 370                                                                     ,clusterPads[2]);
 371                   printf("Indices: %d, %d, %d\n",clusterDigit[0]
 372                                                 ,clusterDigit[1]
 373                                                 ,clusterDigit[2]);
 374                   printf("Total charge = %f\n",clusterCharge);
 375                   printf("Tracks: pad0 %d, %d, %d\n",clusterTracks[0]
 376                                                     ,clusterTracks[1]
 377                                                     ,clusterTracks[2]);
 378                   printf("        pad1 %d, %d, %d\n",clusterTracks[3]
 379                                                     ,clusterTracks[4]
 380                                                     ,clusterTracks[5]);
 381                   printf("        pad2 %d, %d, %d\n",clusterTracks[6]
 382                                                     ,clusterTracks[7]
 383                                                     ,clusterTracks[8]);
 384                 }
 385
 386                 nClusters++;
 387                 if (iType == 1) nClustersUnfold++;
 388
 389                 // Add the cluster to the output array
 390                 fTRD->AddCluster(clusterPads
 391                                 ,clusterDigit
 392                                 ,idet
 393                                 ,clusterCharge
 394                                 ,clusterTracks
 395                                 ,iType);
 396
 397               }
 398             }
 399           }
 400         }
 401
 402         // Compress the arrays
 403         digits->Compress(1,0);
 404         track0->Compress(1,0);
 405         track1->Compress(1,0);
 406         track2->Compress(1,0);
 407
 408         // Write the cluster and reset the array
 409         WriteClusters(idet);
 410         fTRD->ResetRecPoints();
 411
 412         printf("AliTRDclusterizerV1::MakeCluster -- ");
 413         printf("Found %d (%d unfolded) clusters.\n"
 414               ,nClusters,nClustersUnfold);
 415
 416       }
 417     }
 418   }
 419
 420   printf("AliTRDclusterizerV1::MakeCluster -- ");
 421   printf("Done.\n");
 422
 423   return kTRUE;
 424
 425 }
 426
 427 //_____________________________________________________________________________
 428 Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal)
 429 {
 430   //
 431   // Method to unfold neighbouring maxima.
 432   // The charge ratio on the overlapping pad is calculated
 433   // until there is no more change within the range given by eps.
 434   // The resulting ratio is then returned to the calling method.
 435   //
 436
 437   Int_t   itStep            = 0;      // Count iteration steps
 438
 439   Float_t ratio             = 0.5;    // Start value for ratio
 440   Float_t prevRatio         = 0;      // Store previous ratio
 441
 442   Float_t newLeftSignal[3]  = {0};    // Array to store left cluster signal
 443   Float_t newRightSignal[3] = {0};    // Array to store right cluster signal
 444
 445   // Start the iteration
 446   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
 447
 448     itStep++;
 449     prevRatio = ratio;
 450
 451     // Cluster position according to charge ratio
 452     Float_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
 453                      / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
 454     Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
 455                      / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
 456
 457     // Set cluster charge ratio
 458     Float_t ampLeft  = padSignal[1];
 459     Float_t ampRight = padSignal[3];
 460
 461     // Apply pad response to parameters
 462     newLeftSignal[0]  = ampLeft  * PadResponse(-1 - maxLeft);
 463     newLeftSignal[1]  = ampLeft  * PadResponse( 0 - maxLeft);
 464     newLeftSignal[2]  = ampLeft  * PadResponse( 1 - maxLeft);
 465
 466     newRightSignal[0] = ampRight * PadResponse(-1 - maxRight);
 467     newRightSignal[1] = ampRight * PadResponse( 0 - maxRight);
 468     newRightSignal[2] = ampRight * PadResponse( 1 - maxRight);
 469
 470     // Calculate new overlapping ratio
 471     ratio = newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]);
 472
 473   }
 474
 475   return ratio;
 476
 477 }
 478
 479 //_____________________________________________________________________________
 480 Float_t AliTRDclusterizerV1::PadResponse(Float_t x)
 481 {
 482   //
 483   // The pad response for the chevron pads.
 484   // We use a simple Gaussian approximation which should be good
 485   // enough for our purpose.
 486   // Updated for new PRF 1/5/01.
 487   //
 488
 489   // The parameters for the response function
 490   const Float_t kA  =  0.8303;
 491   const Float_t kB  = -0.00392;
 492   const Float_t kC  =  0.472 * 0.472;
 493   const Float_t kD  =  2.19;
 494
 495   Float_t pr = kA * (kB + TMath::Exp(-TMath::Power(x*x,kD) / (2.*kC)));
 496
 497   return (pr);
 498
 499 }