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