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.6  2000/06/09 11:10:07  cblume
  19 Compiler warnings and coding conventions, next round
  20
  21 Revision 1.5  2000/06/08 18:32:58  cblume
  22 Make code compliant to coding conventions
  23
  24 Revision 1.4  2000/06/07 16:27:01  cblume
  25 Try to remove compiler warnings on Sun and HP
  26
  27 Revision 1.3  2000/05/08 16:17:27  cblume
  28 Merge TRD-develop
  29
  30 Revision 1.1.4.1  2000/05/08 15:09:01  cblume
  31 Introduce AliTRDdigitsManager
  32
  33 Revision 1.1  2000/02/28 18:58:54  cblume
  34 Add new TRD classes
  35
  36 */
  37
  38 ///////////////////////////////////////////////////////////////////////////////
  39 //                                                                           //
  40 // TRD cluster finder for the slow simulator.
  41 //                                                                           //
  42 ///////////////////////////////////////////////////////////////////////////////
  43
  44 #include <TF1.h>
  45
  46 #include "AliTRDclusterizerV1.h"
  47 #include "AliTRDmatrix.h"
  48 #include "AliTRDgeometry.h"
  49 #include "AliTRDdigitizer.h"
  50 #include "AliTRDrecPoint.h"
  51 #include "AliTRDdataArrayF.h"
  52
  53 ClassImp(AliTRDclusterizerV1)
  54
  55 //_____________________________________________________________________________
  56 AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
  57 {
  58   //
  59   // AliTRDclusterizerV1 default constructor
  60   //
  61
  62   fDigitsManager = NULL;
  63
  64 }
  65
  66 //_____________________________________________________________________________
  67 AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title)
  68                     :AliTRDclusterizer(name,title)
  69 {
  70   //
  71   // AliTRDclusterizerV1 default constructor
  72   //
  73
  74   fDigitsManager = new AliTRDdigitsManager();
  75
  76   Init();
  77
  78 }
  79
  80 //_____________________________________________________________________________
  81 AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
  82 {
  83   //
  84   // AliTRDclusterizerV1 copy constructor
  85   //
  86
  87   ((AliTRDclusterizerV1 &) c).Copy(*this);
  88
  89 }
  90
  91 //_____________________________________________________________________________
  92 AliTRDclusterizerV1::~AliTRDclusterizerV1()
  93 {
  94   //
  95   // AliTRDclusterizerV1 destructor
  96   //
  97
  98   if (fDigitsManager) {
  99     delete fDigitsManager;
 100   }
 101
 102 }
 103
 104 //_____________________________________________________________________________
 105 AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c)
 106 {
 107   //
 108   // Assignment operator
 109   //
 110
 111   if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this);
 112   return *this;
 113
 114 }
 115
 116 //_____________________________________________________________________________
 117 void AliTRDclusterizerV1::Copy(TObject &c)
 118 {
 119   //
 120   // Copy function
 121   //
 122
 123   ((AliTRDclusterizerV1 &) c).fClusMaxThresh = fClusMaxThresh;
 124   ((AliTRDclusterizerV1 &) c).fClusSigThresh = fClusSigThresh;
 125   ((AliTRDclusterizerV1 &) c).fClusMethod    = fClusMethod;
 126   ((AliTRDclusterizerV1 &) c).fDigitsManager = NULL;
 127
 128   AliTRDclusterizer::Copy(c);
 129
 130 }
 131
 132 //_____________________________________________________________________________
 133 void AliTRDclusterizerV1::Init()
 134 {
 135   //
 136   // Initializes the cluster finder
 137   //
 138
 139   // The default parameter for the clustering
 140   fClusMaxThresh = 5.0;
 141   fClusSigThresh = 2.0;
 142   fClusMethod    = 1;
 143
 144 }
 145
 146 //_____________________________________________________________________________
 147 Bool_t AliTRDclusterizerV1::ReadDigits()
 148 {
 149   //
 150   // Reads the digits arrays from the input aliroot file
 151   //
 152
 153   if (!fInputFile) {
 154     printf("AliTRDclusterizerV1::ReadDigits -- ");
 155     printf("No input file open\n");
 156     return kFALSE;
 157   }
 158
 159   // Read in the digit arrays
 160   return (fDigitsManager->ReadDigits());
 161
 162 }
 163
 164 //_____________________________________________________________________________
 165 Bool_t AliTRDclusterizerV1::MakeCluster()
 166 {
 167   //
 168   // Generates the cluster.
 169   //
 170
 171   Int_t row, col, time;
 172
 173   // Get the pointer to the detector class and check for version 1
 174   AliTRD *trd = (AliTRD*) gAlice->GetDetector("TRD");
 175   if (trd->IsVersion() != 1) {
 176     printf("AliTRDclusterizerV1::MakeCluster -- ");
 177     printf("TRD must be version 1 (slow simulator).\n");
 178     return kFALSE;
 179   }
 180
 181   // Get the geometry
 182   AliTRDgeometry *geo = trd->GetGeometry();
 183
 184   printf("AliTRDclusterizerV1::MakeCluster -- ");
 185   printf("Start creating clusters.\n");
 186
 187   AliTRDdataArrayI *digits;
 188
 189   // Parameters
 190   Float_t maxThresh        = fClusMaxThresh;   // threshold value for maximum
 191   Float_t signalThresh     = fClusSigThresh;   // threshold value for digit signal
 192   Int_t   clusteringMethod = fClusMethod;      // clustering method option (for testing)
 193
 194   // Iteration limit for unfolding procedure
 195   const Float_t kEpsilon = 0.01;
 196
 197   const Int_t   kNclus   = 3;
 198   const Int_t   kNsig    = 5;
 199
 200   Int_t chamBeg = 0;
 201   Int_t chamEnd = kNcham;
 202   if (trd->GetSensChamber()  >= 0) {
 203     chamBeg = trd->GetSensChamber();
 204     chamEnd = chamBeg + 1;
 205   }
 206   Int_t planBeg = 0;
 207   Int_t planEnd = kNplan;
 208   if (trd->GetSensPlane()    >= 0) {
 209     planBeg = trd->GetSensPlane();
 210     planEnd = planBeg + 1;
 211   }
 212   Int_t sectBeg = 0;
 213   Int_t sectEnd = kNsect;
 214
 215   // *** Start clustering *** in every chamber
 216   for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
 217     for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
 218       for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
 219
 220         if (trd->GetSensSector() >= 0) {
 221           Int_t sens1 = trd->GetSensSector();
 222           Int_t sens2 = sens1 + trd->GetSensSectorRange();
 223           sens2 -= ((Int_t) (sens2 / kNsect)) * kNsect;
 224           if (sens1 < sens2) {
 225             if ((isect < sens1) || (isect >= sens2)) continue;
 226           }
 227           else {
 228             if ((isect < sens1) && (isect >= sens2)) continue;
 229           }
 230         }
 231
 232         Int_t idet = geo->GetDetector(iplan,icham,isect);
 233
 234         Int_t nClusters = 0;
 235         printf("AliTRDclusterizerV1::MakeCluster -- ");
 236         printf("Analyzing chamber %d, plane %d, sector %d.\n"
 237                ,icham,iplan,isect);
 238
 239         Int_t   nRowMax  = geo->GetRowMax(iplan,icham,isect);
 240         Int_t   nColMax  = geo->GetColMax(iplan);
 241         Int_t   nTimeMax = geo->GetTimeMax();
 242
 243         // Create a detector matrix to keep maxima
 244         AliTRDmatrix *digitMatrix  = new AliTRDmatrix(nRowMax,nColMax,nTimeMax
 245                                                      ,isect,icham,iplan);
 246         // Create a matrix to contain maximum flags
 247         AliTRDmatrix *maximaMatrix = new AliTRDmatrix(nRowMax,nColMax,nTimeMax
 248                                                      ,isect,icham,iplan);
 249
 250         // Read in the digits
 251         digits = fDigitsManager->GetDigits(idet);
 252
 253         // Loop through the detector pixel
 254         for (time = 0; time < nTimeMax; time++) {
 255           for ( col = 0;  col <  nColMax;  col++) {
 256             for ( row = 0;  row <  nRowMax;  row++) {
 257
 258               Int_t signal = digits->GetData(row,col,time);
 259               Int_t index  = digits->GetIndex(row,col,time);
 260
 261               // Fill the detector matrix
 262               if (signal > signalThresh) {
 263                 // Store the signal amplitude
 264                 digitMatrix->SetSignal(row,col,time,signal);
 265                 // Store the digits number
 266                 digitMatrix->AddTrack(row,col,time,index);
 267               }
 268
 269             }
 270           }
 271         }
 272
 273         // Loop chamber and find maxima in digitMatrix
 274         for ( row = 0;  row <  nRowMax;  row++) {
 275           for ( col = 1;  col <  nColMax;  col++) {
 276             for (time = 0; time < nTimeMax; time++) {
 277
 278               if (digitMatrix->GetSignal(row,col,time)
 279                   < digitMatrix->GetSignal(row,col - 1,time)) {
 280                 // really maximum?
 281                 if (col > 1) {
 282                   if (digitMatrix->GetSignal(row,col - 2,time)
 283                       < digitMatrix->GetSignal(row,col - 1,time)) {
 284                     // yes, so set maximum flag
 285                     maximaMatrix->SetSignal(row,col - 1,time,1);
 286                   }
 287                   else maximaMatrix->SetSignal(row,col - 1,time,0);
 288                 }
 289               }
 290
 291             }   // time
 292           }     // col
 293         }       // row
 294
 295         // now check maxima and calculate cluster position
 296         for ( row = 0;  row <  nRowMax;  row++) {
 297           for ( col = 1;  col <  nColMax;  col++) {
 298             for (time = 0; time < nTimeMax; time++) {
 299
 300               if ((maximaMatrix->GetSignal(row,col,time) > 0)
 301                   && (digitMatrix->GetSignal(row,col,time) > maxThresh)) {
 302
 303                 // Ratio resulting from unfolding
 304                 Float_t ratio                 =  0;
 305                 // Signals on max and neighbouring pads
 306                 Float_t padSignal[kNsig]      = {0};
 307                 // Signals from cluster
 308                 Float_t clusterSignal[kNclus] = {0};
 309                 // Cluster pad info
 310                 Float_t clusterPads[kNclus]   = {0};
 311                 // Cluster digit info
 312                 Int_t   clusterDigit[kNclus]  = {0};
 313
 314                 Int_t iPad;
 315                 for (iPad = 0; iPad < kNclus; iPad++) {
 316                   clusterSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time);
 317                   clusterDigit[iPad]  = digitMatrix->GetTrack(row,col-1+iPad,time,0);
 318                 }
 319
 320                 // neighbouring maximum on right side?
 321                 if (col < nColMax - 2) {
 322                   if (maximaMatrix->GetSignal(row,col + 2,time) > 0) {
 323
 324                     for (iPad = 0; iPad < 5; iPad++) {
 325                       padSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time);
 326                     }
 327
 328                     // unfold:
 329                     ratio = Unfold(kEpsilon, padSignal);
 330
 331                     // set signal on overlapping pad to ratio
 332                     clusterSignal[2] *= ratio;
 333
 334                   }
 335                 }
 336
 337                 // Calculate the position of the cluster
 338                 switch (clusteringMethod) {
 339                 case 1:
 340                   // method 1: simply center of mass
 341                   clusterPads[0] = row + 0.5;
 342                   clusterPads[1] = col - 0.5 + (clusterSignal[2] - clusterSignal[0]) /
 343                                    (clusterSignal[0] + clusterSignal[1] + clusterSignal[2]);
 344                   clusterPads[2] = time + 0.5;
 345
 346                   nClusters++;
 347                   break;
 348                 case 2:
 349                   // method 2: integral gauss fit on 3 pads
 350                   TH1F *hPadCharges = new TH1F("hPadCharges", "Charges on center 3 pads"
 351                                                             , 5, -1.5, 3.5);
 352                   for (Int_t iCol = -1; iCol <= 3; iCol++) {
 353                     if (clusterSignal[iCol] < 1) clusterSignal[iCol] = 1;
 354                     hPadCharges->Fill(iCol, clusterSignal[iCol]);
 355                   }
 356                   hPadCharges->Fit("gaus", "IQ", "SAME", -0.5, 2.5);
 357                   TF1     *fPadChargeFit = hPadCharges->GetFunction("gaus");
 358                   Double_t  colMean = fPadChargeFit->GetParameter(1);
 359
 360                   clusterPads[0] = row + 0.5;
 361                   clusterPads[1] = col - 1.5 + colMean;
 362                   clusterPads[2] = time + 0.5;
 363
 364                   delete hPadCharges;
 365
 366                   nClusters++;
 367                   break;
 368                 }
 369
 370                 Float_t clusterCharge =   clusterSignal[0]
 371                                         + clusterSignal[1]
 372                                         + clusterSignal[2];
 373
 374                 // Add the cluster to the output array
 375                 trd->AddRecPoint(clusterPads,clusterDigit,idet,clusterCharge);
 376
 377               }
 378             }  // time
 379           }    // col
 380         }      // row
 381
 382         printf("AliTRDclusterizerV1::MakeCluster -- ");
 383         printf("Number of clusters found: %d\n",nClusters);
 384
 385         delete digitMatrix;
 386         delete maximaMatrix;
 387
 388       }          // isect
 389     }            // iplan
 390   }              // icham
 391
 392   printf("AliTRDclusterizerV1::MakeCluster -- ");
 393   printf("Total number of points found: %d\n"
 394         ,trd->RecPoints()->GetEntries());
 395
 396   // Get the pointer to the cluster branch
 397   TTree *clusterTree = gAlice->TreeR();
 398
 399   // Fill the cluster-branch
 400   printf("AliTRDclusterizerV1::MakeCluster -- ");
 401   printf("Fill the cluster tree.\n");
 402   clusterTree->Fill();
 403   printf("AliTRDclusterizerV1::MakeCluster -- ");
 404   printf("Done.\n");
 405
 406   return kTRUE;
 407
 408 }
 409
 410 //_____________________________________________________________________________
 411 Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal)
 412 {
 413   //
 414   // Method to unfold neighbouring maxima.
 415   // The charge ratio on the overlapping pad is calculated
 416   // until there is no more change within the range given by eps.
 417   // The resulting ratio is then returned to the calling method.
 418   //
 419
 420   Int_t   itStep            = 0;      // count iteration steps
 421
 422   Float_t ratio             = 0.5;    // start value for ratio
 423   Float_t prevRatio         = 0;      // store previous ratio
 424
 425   Float_t newLeftSignal[3]  = {0};    // array to store left cluster signal
 426   Float_t newRightSignal[3] = {0};    // array to store right cluster signal
 427
 428   // start iteration:
 429   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
 430
 431     itStep++;
 432     prevRatio = ratio;
 433
 434     // cluster position according to charge ratio
 435     Float_t maxLeft  = (ratio*padSignal[2] - padSignal[0]) /
 436                        (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
 437     Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) /
 438                        ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
 439
 440     // set cluster charge ratio
 441     Float_t ampLeft  = padSignal[1];
 442     Float_t ampRight = padSignal[3];
 443
 444     // apply pad response to parameters
 445     newLeftSignal[0] = ampLeft*PadResponse(-1 - maxLeft);
 446     newLeftSignal[1] = ampLeft*PadResponse( 0 - maxLeft);
 447     newLeftSignal[2] = ampLeft*PadResponse( 1 - maxLeft);
 448
 449     newRightSignal[0] = ampRight*PadResponse(-1 - maxRight);
 450     newRightSignal[1] = ampRight*PadResponse( 0 - maxRight);
 451     newRightSignal[2] = ampRight*PadResponse( 1 - maxRight);
 452
 453     // calculate new overlapping ratio
 454     ratio = newLeftSignal[2]/(newLeftSignal[2] + newRightSignal[0]);
 455
 456   }
 457
 458   return ratio;
 459
 460 }
 461
 462 //_____________________________________________________________________________
 463 Float_t AliTRDclusterizerV1::PadResponse(Float_t x)
 464 {
 465   //
 466   // The pad response for the chevron pads.
 467   // We use a simple Gaussian approximation which should be good
 468   // enough for our purpose.
 469   //
 470
 471   // The parameters for the response function
 472   const Float_t kA  =  0.8872;
 473   const Float_t kB  = -0.00573;
 474   const Float_t kC  =  0.454;
 475   const Float_t kC2 =  kC*kC;
 476
 477   Float_t pr = kA * (kB + TMath::Exp(-x*x / (2. * kC2)));
 478
 479   return (pr);
 480
 481 }