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