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