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