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