TRD/AliTRDclusterizerV0.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.10.6.1  2002/07/24 10:09:30  alibrary
  19 Updating VirtualMC
  20
  21 Revision 1.11  2002/06/12 09:54:35  cblume
  22 Update of tracking code provided by Sergei
  23
  24 Revision 1.10  2001/12/05 15:04:34  hristov
  25 Changes related to the corrections of AliRecPoint
  26
  27 Revision 1.9  2001/05/28 17:07:58  hristov
  28 Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh)
  29
  30 Revision 1.8  2001/05/07 08:06:44  cblume
  31 Speedup of the code. Create only AliTRDcluster
  32
  33 Revision 1.7  2001/02/14 18:22:26  cblume
  34 Change in the geometry of the padplane
  35
  36 Revision 1.6  2000/11/01 14:53:20  cblume
  37 Merge with TRD-develop
  38
  39 Revision 1.1.4.6  2000/10/16 01:16:53  cblume
  40 Changed timebin 0 to be the one closest to the readout
  41
  42 Revision 1.1.4.5  2000/10/15 23:40:01  cblume
  43 Remove AliTRDconst
  44
  45 Revision 1.1.4.4  2000/10/06 16:49:46  cblume
  46 Made Getters const
  47
  48 Revision 1.1.4.3  2000/10/04 16:34:58  cblume
  49 Replace include files by forward declarations
  50
  51 Revision 1.1.4.2  2000/09/22 14:49:49  cblume
  52 Adapted to tracking code
  53
  54 Revision 1.5  2000/10/02 21:28:19  fca
  55 Removal of useless dependecies via forward declarations
  56
  57 Revision 1.4  2000/06/08 18:32:58  cblume
  58 Make code compliant to coding conventions
  59
  60 Revision 1.3  2000/06/07 16:27:01  cblume
  61 Try to remove compiler warnings on Sun and HP
  62
  63 Revision 1.2  2000/05/08 16:17:27  cblume
  64 Merge TRD-develop
  65
  66 Revision 1.1.4.1  2000/05/08 15:08:41  cblume
  67 Replace AliTRDcluster by AliTRDrecPoint
  68
  69 Revision 1.4  2000/06/08 18:32:58  cblume
  70 Make code compliant to coding conventions
  71
  72 Revision 1.3  2000/06/07 16:27:01  cblume
  73 Try to remove compiler warnings on Sun and HP
  74
  75 Revision 1.2  2000/05/08 16:17:27  cblume
  76 Merge TRD-develop
  77
  78 Revision 1.1.4.1  2000/05/08 15:08:41  cblume
  79 Replace AliTRDcluster by AliTRDrecPoint
  80
  81 Revision 1.1  2000/02/28 18:58:33  cblume
  82 Add new TRD classes
  83
  84 */
  85
  86 ///////////////////////////////////////////////////////////////////////////////
  87 //                                                                           //
  88 // TRD cluster finder for the fast simulator. It takes the hits from the     //
  89 // fast simulator (one hit per plane) and transforms them                    //
  90 // into cluster, by applying position smearing and merging                   //
  91 // of nearby cluster. The searing is done uniformly in z-direction           //
  92 // over the length of a readout pad. In rphi-direction a Gaussian            //
  93 // smearing is applied with a sigma given by fRphiSigma.                     //
  94 // Clusters are considered as overlapping when they are closer in            //
  95 // rphi-direction than the value defined in fRphiDist.                       //
  96 // Use the macro fastClusterCreate.C to create the cluster.                  //
  97 //                                                                           //
  98 ///////////////////////////////////////////////////////////////////////////////
  99
 100 #include <TRandom.h>
 101 #include <TTree.h>
 102
 103 #include "AliRun.h"
 104
 105 #include "AliTRD.h"
 106 #include "AliTRDclusterizerV0.h"
 107 #include "AliTRDhit.h"
 108 #include "AliTRDgeometry.h"
 109 #include "AliTRDrecPoint.h"
 110 #include "AliTRDparameter.h"
 111
 112 ClassImp(AliTRDclusterizerV0)
 113
 114 //_____________________________________________________________________________
 115 AliTRDclusterizerV0::AliTRDclusterizerV0():AliTRDclusterizer()
 116 {
 117   //
 118   // AliTRDclusterizerV0 default constructor
 119   //
 120
 121 }
 122
 123 //_____________________________________________________________________________
 124 AliTRDclusterizerV0::AliTRDclusterizerV0(const Text_t* name, const Text_t* title)
 125                     :AliTRDclusterizer(name,title)
 126 {
 127   //
 128   // AliTRDclusterizerV0 default constructor
 129   //
 130
 131   Init();
 132
 133 }
 134
 135 //_____________________________________________________________________________
 136 AliTRDclusterizerV0::~AliTRDclusterizerV0()
 137 {
 138   //
 139   // AliTRDclusterizerV0 destructor
 140   //
 141
 142 }
 143
 144 //_____________________________________________________________________________
 145 void AliTRDclusterizerV0::Init()
 146 {
 147   //
 148   // Initializes the cluster finder
 149   //
 150
 151   // Position resolution in rphi-direction
 152   fRphiSigma  = 0.02;
 153   // Minimum distance of non-overlapping cluster
 154   fRphiDist   = 1.0;
 155
 156 }
 157
 158 //_____________________________________________________________________________
 159 Bool_t AliTRDclusterizerV0::MakeClusters()
 160 {
 161   //
 162   // Generates the cluster
 163   //
 164
 165   if (fTRD->IsVersion() != 0) {
 166     printf("AliTRDclusterizerV0::MakeCluster -- ");
 167     printf("TRD must be version 0 (fast simulator).\n");
 168     return kFALSE;
 169   }
 170
 171   // Get the geometry
 172   AliTRDgeometry *geo = fTRD->GetGeometry();
 173
 174   // Create a default parameter class if none is defined
 175   if (!fPar) {
 176     fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter");
 177     if (fVerbose > 0) {
 178       printf("<AliTRDclusterizerV0::MakeCluster> ");
 179       printf("Create the default parameter object.\n");
 180     }
 181   }
 182
 183   printf("AliTRDclusterizerV0::MakeCluster -- ");
 184   printf("Start creating cluster.\n");
 185
 186   Int_t nBytes = 0;
 187
 188   AliTRDhit *hit;
 189
 190   // Get the pointer to the hit tree
 191   TTree     *hitTree      = gAlice->TreeH();
 192   // Get the pointer to the reconstruction tree
 193   TTree     *clusterTree  = gAlice->TreeR();
 194
 195   TObjArray *chamberArray = new TObjArray();
 196
 197   // Get the number of entries in the hit tree
 198   // (Number of primary particles creating a hit somewhere)
 199   Int_t nTrack = (Int_t) hitTree->GetEntries();
 200
 201   // Loop through all the chambers
 202   for (Int_t icham = 0; icham < AliTRDgeometry::Ncham(); icham++) {
 203     for (Int_t iplan = 0; iplan < AliTRDgeometry::Nplan(); iplan++) {
 204       for (Int_t isect = 0; isect < AliTRDgeometry::Nsect(); isect++) {
 205
 206         Int_t   nColMax     = fPar->GetColMax(iplan);
 207         Float_t row0        = fPar->GetRow0(iplan,icham,isect);
 208         Float_t col0        = fPar->GetCol0(iplan);
 209         Float_t time0       = fPar->GetTime0(iplan);
 210
 211         Float_t rowPadSize  = fPar->GetRowPadSize(iplan,icham,isect);
 212         Float_t colPadSize  = fPar->GetColPadSize(iplan);
 213         Float_t timeBinSize = fPar->GetTimeBinSize();
 214
 215         // Loop through all entries in the tree
 216         for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
 217
 218           gAlice->ResetHits();
 219           nBytes += hitTree->GetEvent(iTrack);
 220
 221           // Get the number of hits in the TRD created by this particle
 222           Int_t nHit = fTRD->Hits()->GetEntriesFast();
 223
 224           // Loop through the TRD hits
 225           for (Int_t iHit = 0; iHit < nHit; iHit++) {
 226
 227             if (!(hit = (AliTRDhit *) fTRD->Hits()->UncheckedAt(iHit)))
 228               continue;
 229
 230             Float_t pos[3];
 231                     pos[0]   = hit->X();
 232                     pos[1]   = hit->Y();
 233                     pos[2]   = hit->Z();
 234             Int_t   track    = hit->Track();
 235             Int_t   detector = hit->GetDetector();
 236             Int_t   plane    = geo->GetPlane(detector);
 237             Int_t   sector   = geo->GetSector(detector);
 238             Int_t   chamber  = geo->GetChamber(detector);
 239
 240             if ((sector  != isect) ||
 241                 (plane   != iplan) ||
 242                 (chamber != icham))
 243               continue;
 244
 245             // Rotate the sectors on top of each other
 246             Float_t rot[3];
 247             geo->Rotate(detector,pos,rot);
 248
 249             // Add this recPoint to the temporary array for this chamber
 250             AliTRDrecPoint *recPoint = new AliTRDrecPoint("");
 251             recPoint->SetLocalRow(rot[2]);
 252             recPoint->SetLocalCol(rot[1]);
 253             recPoint->SetLocalTime(rot[0]);
 254             recPoint->SetEnergy(0);
 255             recPoint->SetDetector(detector);
 256             recPoint->AddDigit(track);
 257             chamberArray->Add(recPoint);
 258
 259           }
 260
 261         }
 262
 263         // Loop through the temporary cluster-array
 264         for (Int_t iClus1 = 0; iClus1 < chamberArray->GetEntries(); iClus1++) {
 265
 266           AliTRDrecPoint *recPoint1 = (AliTRDrecPoint *)
 267                                       chamberArray->UncheckedAt(iClus1);
 268           Float_t row1  = recPoint1->GetLocalRow();
 269           Float_t col1  = recPoint1->GetLocalCol();
 270           Float_t time1 = recPoint1->GetLocalTime();
 271
 272           if (recPoint1->GetEnergy() < 0) continue;        // Skip marked cluster
 273
 274           const Int_t kNsave  = 5;
 275           Int_t idxSave[kNsave];
 276           Int_t iSave = 0;
 277
 278           const Int_t kNsaveTrack = 3;
 279           Int_t tracks[kNsaveTrack];
 280           tracks[0] = recPoint1->GetDigit(0);
 281
 282           // Check the other cluster to see, whether there are close ones
 283           for (Int_t iClus2 = iClus1 + 1; iClus2 < chamberArray->GetEntries(); iClus2++) {
 284
 285             AliTRDrecPoint *recPoint2 = (AliTRDrecPoint *)
 286                                         chamberArray->UncheckedAt(iClus2);
 287             Float_t row2 = recPoint2->GetLocalRow();
 288             Float_t col2 = recPoint2->GetLocalCol();
 289
 290             if ((TMath::Abs(row1 - row2) < rowPadSize) ||
 291                 (TMath::Abs(col1 - col2) <  fRphiDist)) {
 292               if (iSave == kNsave) {
 293                 printf("AliTRDclusterizerV0::MakeCluster -- ");
 294                 printf("Boundary error: iSave = %d, kNsave = %d.\n"
 295                       ,iSave,kNsave);
 296               }
 297               else {
 298                 idxSave[iSave]  = iClus2;
 299                 iSave++;
 300                 if (iSave < kNsaveTrack) tracks[iSave] = recPoint2->GetDigit(0);
 301               }
 302             }
 303           }
 304
 305           // Merge close cluster
 306           Float_t rowMerge = row1;
 307           Float_t colMerge = col1;
 308           if (iSave) {
 309             for (Int_t iMerge = 0; iMerge < iSave; iMerge++) {
 310               AliTRDrecPoint *recPoint2 =
 311                 (AliTRDrecPoint *) chamberArray->UncheckedAt(idxSave[iMerge]);
 312               rowMerge += recPoint2->GetLocalRow();
 313               colMerge += recPoint2->GetLocalCol();
 314               recPoint2->SetEnergy(-1);     // Mark merged cluster
 315             }
 316             rowMerge /= (iSave + 1);
 317             colMerge /= (iSave + 1);
 318           }
 319
 320           Float_t smear[3];
 321
 322           // The position smearing in row-direction (uniform over pad width)
 323           Int_t row = (Int_t) ((rowMerge - row0) / rowPadSize);
 324           smear[0]  = (row + gRandom->Rndm()) * rowPadSize + row0;
 325
 326           // The position smearing in rphi-direction (Gaussian)
 327           smear[1] = 0;
 328           do
 329             smear[1] = gRandom->Gaus(colMerge,fRphiSigma);
 330           while ((smear[1] < col0                        ) ||
 331                  (smear[1] > col0 + nColMax * colPadSize));
 332
 333           // Time direction stays unchanged
 334           smear[2] = time1;
 335
 336           // Transform into local coordinates
 337           smear[0] = (Int_t) ((smear[0] -  row0) /  rowPadSize);
 338           smear[1] = (Int_t) ((smear[1] -  col0) /  colPadSize);
 339           smear[2] = (Int_t) ((time0 - smear[2]) / timeBinSize);
 340
 341           // Add the smeared cluster to the output array
 342           Int_t   detector  = recPoint1->GetDetector();
 343           Int_t   tr[9]     = { -1   };
 344           Float_t pos[3];
 345           Float_t sigma[2]  = {  0.0 };
 346           pos[0] = smear[1];
 347           pos[1] = smear[0];
 348           pos[2] = (time0 - smear[2]) / timeBinSize;
 349           fTRD->AddCluster(pos,detector,0.0,tr,sigma,0);
 350
 351         }
 352
 353         // Clear the temporary cluster-array and delete the cluster
 354         chamberArray->Delete();
 355
 356       }
 357     }
 358   }
 359
 360   printf("AliTRDclusterizerV0::MakeCluster -- ");
 361   printf("Found %d points.\n",fTRD->RecPoints()->GetEntries());
 362   printf("AliTRDclusterizerV0::MakeCluster -- ");
 363   printf("Fill the cluster tree.\n");
 364   clusterTree->Fill();
 365
 366   return kTRUE;
 367
 368 }