TRD/AliTRDclusterizerMI.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 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 // TRD cluster finder for the slow simulator.
  21 //                                                                           //
  22 ///////////////////////////////////////////////////////////////////////////////
  23
  24 #include <TF1.h>
  25 #include <TTree.h>
  26 #include <TH1.h>
  27 #include <TFile.h>
  28
  29 #include "AliRun.h"
  30 #include "AliRunLoader.h"
  31 #include "AliLoader.h"
  32
  33 #include "AliTRDclusterizerMI.h"
  34 #include "AliTRDmatrix.h"
  35 #include "AliTRDgeometry.h"
  36 #include "AliTRDdataArrayF.h"
  37 #include "AliTRDdataArrayI.h"
  38 #include "AliTRDdigitsManager.h"
  39 #include "AliTRDparameter.h"
  40 #include "AliTRDclusterMI.h"
  41
  42 ClassImp(AliTRDclusterizerMI)
  43
  44 //_____________________________________________________________________________
  45 AliTRDclusterizerMI::AliTRDclusterizerMI():AliTRDclusterizerV1()
  46 {
  47   //
  48   // AliTRDclusterizerMI default constructor
  49   //
  50 }
  51
  52 //_____________________________________________________________________________
  53 AliTRDclusterizerMI::AliTRDclusterizerMI(const Text_t* name, const Text_t* title)
  54                     :AliTRDclusterizerV1(name,title)
  55 {
  56   //
  57   // AliTRDclusterizerMI default constructor
  58   //
  59 }
  60
  61 //_____________________________________________________________________________
  62 AliTRDclusterizerMI::~AliTRDclusterizerMI()
  63 {
  64   //
  65   // AliTRDclusterizerMI destructor
  66   //
  67 }
  68
  69
  70 AliTRDclusterMI *  AliTRDclusterizerMI::AddCluster()
  71 {
  72   AliTRDclusterMI *c = new AliTRDclusterMI();
  73   fClusterContainer->Add(c);
  74   return c;
  75 }
  76
  77 void AliTRDclusterizerMI::SetCluster(AliTRDclusterMI * cl,Float_t *pos, Int_t det, Float_t amp
  78                   ,Int_t *tracks, Float_t *sig, Int_t iType, Float_t sigmay, Float_t relpad)
  79 {
  80   //
  81   //
  82   //
  83   cl->SetDetector(det);
  84   cl->AddTrackIndex(tracks);
  85   cl->SetQ(amp);
  86   cl->SetY(pos[0]);
  87   cl->SetZ(pos[1]);
  88   cl->SetSigmaY2(sig[0]);
  89   cl->SetSigmaZ2(sig[1]);
  90   cl->SetLocalTimeBin(((Int_t) pos[2]));
  91   cl->SetNPads(iType);
  92   cl->SetRelPos(relpad);
  93   cl->fRmsY = sigmay;
  94 }
  95
  96 void AliTRDclusterizerMI::MakeCluster(Float_t * padSignal, Float_t * pos, Float_t &sigma, Float_t & relpad)
  97 {
  98   //
  99   //
 100   Float_t sum   = 0;
 101   Float_t sumx  = 0;
 102   Float_t sumx2 = 0;
 103   Float_t signal[3]={padSignal[0],padSignal[1],padSignal[2]};
 104   if ( signal[0]<2){
 105     signal[0] = 0.015*(signal[1]*signal[1])/(signal[2]+0.5);
 106     if (signal[0]>2) signal[0]=2;
 107   }
 108   if ( signal[2]<2){
 109     signal[2] = 0.015*(signal[1]*signal[1])/(signal[0]+0.5);
 110     if (signal[2]>2) signal[2]=2;
 111   }
 112
 113   for (Int_t i=-1;i<=1;i++){
 114     sum   +=signal[i+1];
 115     sumx  +=signal[i+1]*float(i);
 116     sumx2 +=signal[i+1]*float(i)*float(i);
 117   }
 118
 119   pos[0] = sumx/sum;
 120   sigma  = sumx2/sum-pos[0]*pos[0];
 121   relpad = pos[0];
 122 }
 123
 124 //_____________________________________________________________________________
 125 Bool_t AliTRDclusterizerMI::MakeClusters()
 126 {
 127   //
 128   // Generates the cluster.
 129   //
 130
 131   //////////////////////
 132   //STUPIDITY to be fixed later
 133   fClusterContainer = RecPoints();
 134
 135   Int_t row, col, time;
 136
 137   /*
 138   if (fTRD->IsVersion() != 1) {
 139     printf("<AliTRDclusterizerMI::MakeCluster> ");
 140     printf("TRD must be version 1 (slow simulator).\n");
 141     return kFALSE;
 142   }
 143   */
 144
 145   // Get the geometry
 146   AliTRDgeometry *geo = AliTRDgeometry::GetGeometry(fRunLoader);
 147
 148   // Create a default parameter class if none is defined
 149   if (!fPar) {
 150     fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter");
 151     printf("<AliTRDclusterizerMI::MakeCluster> ");
 152     printf("Create the default parameter object.\n");
 153   }
 154
 155   Float_t timeBinSize = fPar->GetTimeBinSize();
 156   // Half of ampl.region
 157   const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.;
 158
 159   Float_t omegaTau = fPar->GetOmegaTau();
 160   if (fVerbose > 0) {
 161     printf("<AliTRDclusterizerMI::MakeCluster> ");
 162     printf("OmegaTau = %f \n",omegaTau);
 163     printf("<AliTRDclusterizerMI::MakeCluster> ");
 164     printf("Start creating clusters.\n");
 165   }
 166
 167   AliTRDdataArrayI *digits;
 168   AliTRDdataArrayI *track0;
 169   AliTRDdataArrayI *track1;
 170   AliTRDdataArrayI *track2;
 171
 172   // Threshold value for the maximum
 173   Int_t maxThresh = fPar->GetClusMaxThresh();
 174   // Threshold value for the digit signal
 175   Int_t sigThresh = fPar->GetClusSigThresh();
 176
 177   // Iteration limit for unfolding procedure
 178   const Float_t kEpsilon = 0.01;
 179
 180   const Int_t   kNclus   = 3;
 181   const Int_t   kNsig    = 5;
 182   const Int_t   kNtrack  = 3 * kNclus;
 183
 184   Int_t   iType          = 0;
 185   Int_t   iUnfold        = 0;
 186
 187   Float_t ratioLeft      = 1.0;
 188   Float_t ratioRight     = 1.0;
 189
 190   Float_t padSignal[kNsig];
 191   Float_t clusterSignal[kNclus];
 192   Float_t clusterPads[kNclus];
 193   Int_t   clusterDigit[kNclus];
 194   Int_t   clusterTracks[kNtrack];
 195
 196   Int_t chamBeg = 0;
 197   Int_t chamEnd = AliTRDgeometry::Ncham();
 198   Int_t planBeg = 0;
 199   Int_t planEnd = AliTRDgeometry::Nplan();
 200   Int_t sectBeg = 0;
 201   Int_t sectEnd = AliTRDgeometry::Nsect();
 202
 203   // Start clustering in every chamber
 204   for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
 205     for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
 206       for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
 207
 208         Int_t idet = geo->GetDetector(iplan,icham,isect);
 209
 210         Int_t nClusters      = 0;
 211         Int_t nClusters2pad  = 0;
 212         Int_t nClusters3pad  = 0;
 213         Int_t nClusters4pad  = 0;
 214         Int_t nClusters5pad  = 0;
 215         Int_t nClustersLarge = 0;
 216
 217         if (fVerbose > 0) {
 218           printf("<AliTRDclusterizerMI::MakeCluster> ");
 219           printf("Analyzing chamber %d, plane %d, sector %d.\n"
 220                 ,icham,iplan,isect);
 221         }
 222
 223         Int_t   nRowMax     = fPar->GetRowMax(iplan,icham,isect);
 224         Int_t   nColMax     = fPar->GetColMax(iplan);
 225         Int_t   nTimeBefore = fPar->GetTimeBefore();
 226         Int_t   nTimeTotal  = fPar->GetTimeTotal();
 227
 228         Float_t row0        = fPar->GetRow0(iplan,icham,isect);
 229         Float_t col0        = fPar->GetCol0(iplan);
 230         Float_t rowSize     = fPar->GetRowPadSize(iplan,icham,isect);
 231         Float_t colSize     = fPar->GetColPadSize(iplan);
 232
 233         // Get the digits
 234         digits = fDigitsManager->GetDigits(idet);
 235         digits->Expand();
 236         track0 = fDigitsManager->GetDictionary(idet,0);
 237         track0->Expand();
 238         track1 = fDigitsManager->GetDictionary(idet,1);
 239         track1->Expand();
 240         track2 = fDigitsManager->GetDictionary(idet,2);
 241         track2->Expand();
 242
 243         // Loop through the chamber and find the maxima
 244         for ( row = 0;  row <  nRowMax;    row++) {
 245           for ( col = 2;  col <  nColMax;    col++) {
 246             for (time = 0; time < nTimeTotal; time++) {
 247
 248               Int_t signalL = TMath::Abs(digits->GetDataUnchecked(row,col  ,time));
 249               Int_t signalM = TMath::Abs(digits->GetDataUnchecked(row,col-1,time));
 250               Int_t signalR = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
 251
 252               // Look for the maximum
 253               if (signalM >= maxThresh) {
 254                 if (((signalL >= sigThresh) &&
 255                      (signalL <  signalM))  ||
 256                     ((signalR >= sigThresh) &&
 257                      (signalR <  signalM))) {
 258                   // Maximum found, mark the position by a negative signal
 259                   digits->SetDataUnchecked(row,col-1,time,-signalM);
 260                 }
 261               }
 262
 263             }
 264           }
 265         }
 266
 267         // Now check the maxima and calculate the cluster position
 268         for ( row = 0;  row <  nRowMax  ;  row++) {
 269           for (time = 0; time < nTimeTotal; time++) {
 270             for ( col = 1;  col <  nColMax-1;  col++) {
 271
 272               // Maximum found ?
 273               if (digits->GetDataUnchecked(row,col,time) < 0) {
 274
 275                 Int_t iPad;
 276                 for (iPad = 0; iPad < kNclus; iPad++) {
 277                   Int_t iPadCol = col - 1 + iPad;
 278                   clusterSignal[iPad]     = TMath::Abs(digits->GetDataUnchecked(row
 279                                                                                ,iPadCol
 280                                                                                ,time));
 281                   clusterDigit[iPad]      = digits->GetIndexUnchecked(row,iPadCol,time);
 282                   clusterTracks[3*iPad  ] = track0->GetDataUnchecked(row,iPadCol,time) - 1;
 283                   clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1;
 284                   clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1;
 285                 }
 286
 287                 // Count the number of pads in the cluster
 288                 Int_t nPadCount = 0;
 289                 Int_t ii        = 0;
 290                 while (TMath::Abs(digits->GetDataUnchecked(row,col-ii  ,time))
 291                                                                   >= sigThresh) {
 292                   nPadCount++;
 293                   ii++;
 294                   if (col-ii   <        0) break;
 295                 }
 296                 ii = 0;
 297                 while (TMath::Abs(digits->GetDataUnchecked(row,col+ii+1,time))
 298                                                                   >= sigThresh) {
 299                   nPadCount++;
 300                   ii++;
 301                   if (col+ii+1 >= nColMax) break;
 302                 }
 303
 304                 nClusters++;
 305                 switch (nPadCount) {
 306                 case 2:
 307                   iType = 0;
 308                   nClusters2pad++;
 309                   break;
 310                 case 3:
 311                   iType = 1;
 312                   nClusters3pad++;
 313                   break;
 314                 case 4:
 315                   iType = 2;
 316                   nClusters4pad++;
 317                   break;
 318                 case 5:
 319                   iType = 3;
 320                   nClusters5pad++;
 321                   break;
 322                 default:
 323                   iType = 4;
 324                   nClustersLarge++;
 325                   break;
 326                 };
 327
 328                 // Don't analyze large clusters
 329                 //if (iType == 4) continue;
 330
 331                 // Look for 5 pad cluster with minimum in the middle
 332                 Bool_t fivePadCluster = kFALSE;
 333                 if (col < nColMax-3) {
 334                   if (digits->GetDataUnchecked(row,col+2,time) < 0) {
 335                     fivePadCluster = kTRUE;
 336                   }
 337                   if ((fivePadCluster) && (col < nColMax-5)) {
 338                     if (digits->GetDataUnchecked(row,col+4,time) >= sigThresh) {
 339                       fivePadCluster = kFALSE;
 340                     }
 341                   }
 342                   if ((fivePadCluster) && (col >         1)) {
 343                     if (digits->GetDataUnchecked(row,col-2,time) >= sigThresh) {
 344                       fivePadCluster = kFALSE;
 345                     }
 346                   }
 347                 }
 348
 349                 // 5 pad cluster
 350                 // Modify the signal of the overlapping pad for the left part
 351                 // of the cluster which remains from a previous unfolding
 352                 if (iUnfold) {
 353                   clusterSignal[0] *= ratioLeft;
 354                   iType   = 3;
 355                   iUnfold = 0;
 356                 }
 357
 358                 // Unfold the 5 pad cluster
 359                 if (fivePadCluster) {
 360                   for (iPad = 0; iPad < kNsig; iPad++) {
 361                     padSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
 362                                                                          ,col-1+iPad
 363                                                                          ,time));
 364                   }
 365                   // Unfold the two maxima and set the signal on
 366                   // the overlapping pad to the ratio
 367                   ratioRight        = Unfold(kEpsilon,iplan,padSignal);
 368                   ratioLeft         = 1.0 - ratioRight;
 369                   clusterSignal[2] *= ratioRight;
 370                   iType   = 3;
 371                   iUnfold = 1;
 372                 }
 373
 374                 Float_t clusterCharge = clusterSignal[0]
 375                                       + clusterSignal[1]
 376                                       + clusterSignal[2];
 377
 378                 // The position of the cluster
 379                 clusterPads[0] = row + 0.5;
 380                 // Take the shift of the additional time bins into account
 381                 clusterPads[2] = time - nTimeBefore + 0.5;
 382
 383                 if (fPar->LUTOn()) {
 384
 385                   // Calculate the position of the cluster by using the
 386                   // lookup table method
 387                   clusterPads[1] = col + 0.5
 388                                  + fPar->LUTposition(iplan,clusterSignal[0]
 389                                                           ,clusterSignal[1]
 390                                                           ,clusterSignal[2]);
 391
 392                 }
 393                 else {
 394
 395                   // Calculate the position of the cluster by using the
 396                   // center of gravity method
 397                   clusterPads[1] = col + 0.5
 398                                  + (clusterSignal[2] - clusterSignal[0])
 399                                  / clusterCharge;
 400
 401                 }
 402
 403                 Float_t q0 = clusterSignal[0];
 404                 Float_t q1 = clusterSignal[1];
 405                 Float_t q2 = clusterSignal[2];
 406                 Float_t clusterSigmaY2 = (q1*(q0+q2)+4*q0*q2) /
 407                                          (clusterCharge*clusterCharge);
 408
 409                 // Correct for ExB displacement
 410                 if (fPar->ExBOn()) {
 411                   Int_t   local_time_bin = (Int_t) clusterPads[2];
 412                   Float_t driftLength    = local_time_bin * timeBinSize + kAmWidth;
 413                   Float_t colSize        = fPar->GetColPadSize(iplan);
 414                   Float_t deltaY         = omegaTau*driftLength/colSize;
 415                   clusterPads[1]         = clusterPads[1] - deltaY;
 416                 }
 417
 418
 419                 // Calculate the position and the error
 420                 Float_t clusterPos[3];
 421                 clusterPos[0] = clusterPads[1] * colSize + col0;
 422                 clusterPos[1] = clusterPads[0] * rowSize + row0;
 423                 clusterPos[2] = clusterPads[2];
 424                 Float_t clusterSig[2];
 425                 clusterSig[0] = (clusterSigmaY2 + 1./12.) * colSize*colSize;
 426                 clusterSig[1] = rowSize * rowSize / 12.;
 427                 //
 428                 //
 429                 AliTRDclusterMI * cluster = AddCluster();
 430                 Float_t sigma, relpos;
 431                 MakeCluster(clusterSignal, clusterPos, sigma,relpos);
 432
 433                 clusterPos[0] = clusterPads[1] * colSize + col0;
 434                 clusterPos[1] = clusterPads[0] * rowSize + row0;
 435                 clusterPos[2] = clusterPads[2];
 436                 SetCluster(cluster, clusterPos,idet,clusterCharge,clusterTracks,clusterSig,iType,sigma,relpos);
 437                 // Add the cluster to the output array
 438                 //                fTRD->AddCluster(clusterPos
 439                 //                ,idet
 440                 //                ,clusterCharge
 441                 //                ,clusterTracks
 442                 //              ,clusterSig
 443                 //                ,iType);
 444
 445               }
 446             }
 447           }
 448         }
 449
 450         // Compress the arrays
 451         digits->Compress(1,0);
 452         track0->Compress(1,0);
 453         track1->Compress(1,0);
 454         track2->Compress(1,0);
 455
 456         // Write the cluster and reset the array
 457         WriteClusters(idet);
 458         ResetRecPoints();
 459
 460         if (fVerbose > 0) {
 461           printf("<AliTRDclusterizerMI::MakeCluster> ");
 462           printf("Found %d clusters in total.\n"
 463                 ,nClusters);
 464           printf("                                    2pad:  %d\n",nClusters2pad);
 465           printf("                                    3pad:  %d\n",nClusters3pad);
 466           printf("                                    4pad:  %d\n",nClusters4pad);
 467           printf("                                    5pad:  %d\n",nClusters5pad);
 468           printf("                                    Large: %d\n",nClustersLarge);
 469         }
 470
 471       }
 472     }
 473   }
 474
 475   if (fVerbose > 0) {
 476     printf("<AliTRDclusterizerMI::MakeCluster> ");
 477     printf("Done.\n");
 478   }
 479
 480   return kTRUE;
 481
 482 }
 483
 484 //_____________________________________________________________________________
 485 Float_t AliTRDclusterizerMI::Unfold(Float_t eps, Int_t plane, Float_t* padSignal)
 486 {
 487   //
 488   // Method to unfold neighbouring maxima.
 489   // The charge ratio on the overlapping pad is calculated
 490   // until there is no more change within the range given by eps.
 491   // The resulting ratio is then returned to the calling method.
 492   //
 493
 494   Int_t   irc               = 0;
 495   Int_t   itStep            = 0;      // Count iteration steps
 496
 497   Float_t ratio             = 0.5;    // Start value for ratio
 498   Float_t prevRatio         = 0;      // Store previous ratio
 499
 500   Float_t newLeftSignal[3]  = {0};    // Array to store left cluster signal
 501   Float_t newRightSignal[3] = {0};    // Array to store right cluster signal
 502   Float_t newSignal[3]      = {0};
 503
 504   // Start the iteration
 505   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
 506
 507     itStep++;
 508     prevRatio = ratio;
 509
 510     // Cluster position according to charge ratio
 511     Float_t maxLeft  = (ratio*padSignal[2] - padSignal[0])
 512                      / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
 513     Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
 514                      / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
 515
 516     // Set cluster charge ratio
 517     irc = fPar->PadResponse(1.0,maxLeft ,plane,newSignal);
 518     Float_t ampLeft  = padSignal[1] / newSignal[1];
 519     irc = fPar->PadResponse(1.0,maxRight,plane,newSignal);
 520     Float_t ampRight = padSignal[3] / newSignal[1];
 521
 522     // Apply pad response to parameters
 523     irc = fPar->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
 524     irc = fPar->PadResponse(ampRight,maxRight,plane,newRightSignal);
 525
 526     // Calculate new overlapping ratio
 527     ratio = TMath::Min((Float_t)1.0,newLeftSignal[2] /
 528                           (newLeftSignal[2] + newRightSignal[0]));
 529
 530   }
 531
 532   return ratio;
 533
 534 }
 535
 536
 537