TRD/AliTRDv0.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.11  1999/11/01 20:41:51  fca
  19 Added protections against using the wrong version of FRAME
  20
  21 Revision 1.10  1999/09/29 09:24:35  fca
  22 Introduction of the Copyright and cvs Log
  23
  24 */
  25
  26 ///////////////////////////////////////////////////////////////////////////////
  27 //                                                                           //
  28 //  Transition Radiation Detector version 0 -- fast simulator                //
  29 //                                                                           //
  30 //Begin_Html
  31 /*
  32 <img src="picts/AliTRDfullClass.gif">
  33 */
  34 //End_Html
  35 //                                                                           //
  36 //                                                                           //
  37 ///////////////////////////////////////////////////////////////////////////////
  38
  39 #include <TMath.h>
  40 #include <TRandom.h>
  41 #include <TVector.h>
  42
  43 #include "AliTRDv0.h"
  44 #include "AliRun.h"
  45 #include "AliMC.h"
  46 #include "AliConst.h"
  47
  48 ClassImp(AliTRDv0)
  49
  50 //_____________________________________________________________________________
  51 AliTRDv0::AliTRDv0(const char *name, const char *title)
  52          :AliTRD(name, title)
  53 {
  54   //
  55   // Standard constructor for Transition Radiation Detector version 0
  56   //
  57
  58   fIdSens     = 0;
  59   fHitsOn     = 0;
  60
  61   fIdChamber1 = 0;
  62   fIdChamber2 = 0;
  63   fIdChamber3 = 0;
  64
  65   fRphiSigma  = 0;
  66   fRphiDist   = 0;
  67
  68 }
  69
  70 //_____________________________________________________________________________
  71 void AliTRDv0::CreateGeometry()
  72 {
  73   //
  74   // Create the GEANT geometry for the Transition Radiation Detector - Version 0
  75   // This version covers the full azimuth.
  76   //
  77
  78   // Check that FRAME is there otherwise we have no place where to put the TRD
  79   AliModule* FRAME = gAlice->GetModule("FRAME");
  80   if (!FRAME) return;
  81
  82   // Define the chambers
  83   AliTRD::CreateGeometry();
  84
  85 }
  86
  87 //_____________________________________________________________________________
  88 void AliTRDv0::CreateMaterials()
  89 {
  90   //
  91   // Create materials for the Transition Radiation Detector
  92   //
  93
  94   AliTRD::CreateMaterials();
  95
  96 }
  97
  98 //_____________________________________________________________________________
  99 void AliTRDv0::Hits2Clusters()
 100 {
 101   // A simple cluster generator. It takes the hits from the
 102   // fast simulator (one hit per plane) and transforms them
 103   // into cluster, by applying position smearing and merging
 104   // of nearby cluster. The searing is done uniformly in z-direction
 105   // over the length of a readout pad. In rphi-direction a Gaussian
 106   // smearing is applied with a sigma given by fRphiSigma.
 107   // Clusters are considered as overlapping when they are closer in
 108   // rphi-direction than the value defined in fRphiDist.
 109   // Use the macro fastClusterCreate.C to create the cluster.
 110
 111   printf("AliTRDv0::Hits2Clusters -- Start creating cluster\n");
 112
 113   Int_t nBytes = 0;
 114
 115   AliTRDhit *TRDhit;
 116
 117   // Get the pointer to the hit tree
 118   TTree *HitTree     = gAlice->TreeH();
 119   // Get the pointer to the reconstruction tree
 120   TTree *ClusterTree = gAlice->TreeD();
 121
 122   TObjArray *Chamber = new TObjArray();
 123
 124   // Get the number of entries in the hit tree
 125   // (Number of primary particles creating a hit somewhere)
 126   Int_t nTrack = (Int_t) HitTree->GetEntries();
 127
 128   // Loop through all the chambers
 129   for (Int_t icham = 0; icham < kNcham; icham++) {
 130     for (Int_t iplan = 0; iplan < kNplan; iplan++) {
 131       for (Int_t isect = 0; isect < kNsect; isect++) {
 132
 133         // Loop through all entries in the tree
 134         for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
 135
 136           gAlice->ResetHits();
 137           nBytes += HitTree->GetEvent(iTrack);
 138
 139           // Get the number of hits in the TRD created by this particle
 140           Int_t nHit = fHits->GetEntriesFast();
 141
 142           // Loop through the TRD hits
 143           for (Int_t iHit = 0; iHit < nHit; iHit++) {
 144
 145             if (!(TRDhit = (AliTRDhit *) fHits->UncheckedAt(iHit)))
 146               continue;
 147
 148             Float_t x       = TRDhit->fX;
 149             Float_t y       = TRDhit->fY;
 150             Float_t z       = TRDhit->fZ;
 151             Int_t   track   = TRDhit->fTrack;
 152             Int_t   plane   = TRDhit->fPlane;
 153             Int_t   sector  = TRDhit->fSector;
 154             Int_t   chamber = TRDhit->fChamber;
 155
 156             if ((sector  != isect+1) ||
 157                 (plane   != iplan+1) ||
 158                 (chamber != icham+1))
 159               continue;
 160
 161             // Rotate the sectors on top of each other
 162             Float_t phi  = 2.0 * kPI /  (Float_t) kNsect
 163                                * ((Float_t) sector - 0.5);
 164             Float_t xRot = -x * TMath::Cos(phi) + y * TMath::Sin(phi);
 165             Float_t yRot =  x * TMath::Sin(phi) + y * TMath::Cos(phi);
 166             Float_t zRot =  z;
 167
 168             // Add this cluster to the temporary cluster-array for this chamber
 169             Int_t   tracks[3];
 170             tracks[0] = track;
 171             Int_t   clusters[5];
 172             clusters[0] = sector;
 173             clusters[1] = chamber;
 174             clusters[2] = plane;
 175             clusters[3] = 0;
 176             clusters[4] = 0;
 177             Float_t position[3];
 178             position[0] = zRot;
 179             position[1] = yRot;
 180             position[2] = xRot;
 181             AliTRDcluster *Cluster = new AliTRDcluster(tracks,clusters,position);
 182             Chamber->Add(Cluster);
 183
 184           }
 185
 186         }
 187
 188         // Loop through the temporary cluster-array
 189         for (Int_t iClus1 = 0; iClus1 < Chamber->GetEntries(); iClus1++) {
 190
 191           AliTRDcluster *Cluster1 = (AliTRDcluster *) Chamber->UncheckedAt(iClus1);
 192           Float_t x1 = Cluster1->fX;
 193           Float_t y1 = Cluster1->fY;
 194           Float_t z1 = Cluster1->fZ;
 195
 196           if (!(z1)) continue;             // Skip marked cluster
 197
 198           const Int_t nSave = 2;
 199           Int_t idxSave[nSave];
 200           Int_t iSave = 0;
 201
 202           Int_t tracks[3];
 203           tracks[0] = Cluster1->fTracks[0];
 204
 205           // Check the other cluster to see, whether there are close ones
 206           for (Int_t iClus2 = iClus1 + 1; iClus2 < Chamber->GetEntries(); iClus2++) {
 207             AliTRDcluster *Cluster2 = (AliTRDcluster *) Chamber->UncheckedAt(iClus2);
 208             Float_t x2 = Cluster2->fX;
 209             Float_t y2 = Cluster2->fY;
 210             if ((TMath::Abs(x1 - x2) < fRowPadSize) ||
 211                 (TMath::Abs(y1 - y2) <   fRphiDist)) {
 212               if (iSave == nSave) {
 213                 printf("AliTRDv0::Hits2Clusters -- Boundary error: iSave = %d, nSave = %d\n"
 214                       ,iSave,nSave);
 215               }
 216               else {
 217                 idxSave[iSave]  = iClus2;
 218                 tracks[iSave+1] = Cluster2->fTracks[0];
 219               }
 220               iSave++;
 221             }
 222           }
 223
 224           // Merge close cluster
 225           Float_t yMerge = y1;
 226           Float_t xMerge = x1;
 227           if (iSave) {
 228             for (Int_t iMerge = 0; iMerge < iSave; iMerge++) {
 229               AliTRDcluster *Cluster2 = (AliTRDcluster *) Chamber->UncheckedAt(idxSave[iMerge]);
 230               xMerge += Cluster2->fX;
 231               yMerge += Cluster2->fY;
 232               Cluster2->fZ = 0;            // Mark merged cluster
 233             }
 234             xMerge /= (iSave + 1);
 235             yMerge /= (iSave + 1);
 236           }
 237
 238           // The position smearing in z-direction (uniform over pad width)
 239           Int_t row = (Int_t) ((xMerge - fRow0[iplan][icham][isect]) / fRowPadSize);
 240           Float_t xSmear = (row + gRandom->Rndm()) * fRowPadSize
 241                          + fRow0[iplan][icham][isect];
 242
 243           // The position smearing in rphi-direction (Gaussian)
 244           Float_t ySmear = 0;
 245           do
 246             ySmear = gRandom->Gaus(yMerge,fRphiSigma);
 247           while ((ySmear < fCol0[iplan])                               ||
 248                  (ySmear > fCol0[iplan] + fColMax[iplan] * fColPadSize));
 249
 250           // Time direction stays unchanged
 251           Float_t zSmear = z1;
 252
 253           Int_t   clusters[5];
 254           clusters[0] = Cluster1->fSector;
 255           clusters[1] = Cluster1->fChamber;
 256           clusters[2] = Cluster1->fPlane;
 257           clusters[3] = 0;
 258           clusters[4] = 0;
 259           Float_t position[3];
 260           // Rotate the sectors back into their real position
 261           Float_t phi = 2*kPI / kNsect * ((Float_t) Cluster1->fSector - 0.5);
 262           position[0] = -zSmear * TMath::Cos(phi) + ySmear * TMath::Sin(phi);
 263           position[1] =  zSmear * TMath::Sin(phi) + ySmear * TMath::Cos(phi);
 264           position[2] =  xSmear;
 265
 266           // Add the smeared cluster to the output array
 267           AddCluster(tracks,clusters,position);
 268
 269         }
 270
 271         // Clear the temporary cluster-array and delete the cluster
 272         Chamber->Delete();
 273
 274       }
 275     }
 276   }
 277
 278   printf("AliTRDv0::Hits2Clusters -- Found %d cluster\n",fClusters->GetEntries());
 279   printf("AliTRDv0::Hits2Clusters -- Fill the cluster tree\n");
 280   ClusterTree->Fill();
 281
 282 }
 283
 284 //_____________________________________________________________________________
 285 void AliTRDv0::Init()
 286 {
 287   //
 288   // Initialise Transition Radiation Detector after geometry is built
 289   //
 290
 291   AliTRD::Init();
 292
 293   // Identifier of the sensitive volume (amplification region)
 294   fIdSens     = gMC->VolId("UL06");
 295
 296   // Identifier of the TRD-driftchambers
 297   fIdChamber1 = gMC->VolId("UCIO");
 298   fIdChamber2 = gMC->VolId("UCIM");
 299   fIdChamber3 = gMC->VolId("UCII");
 300
 301   // Parameter for Hits2Cluster
 302
 303   // Position resolution in rphi-direction
 304   fRphiSigma  = 0.02;
 305   // Minimum distance of non-overlapping cluster
 306   fRphiDist   = 1.0;
 307
 308   printf("          Fast simulator\n");
 309   for (Int_t i = 0; i < 80; i++) printf("*");
 310   printf("\n");
 311
 312 }
 313
 314 //_____________________________________________________________________________
 315 void AliTRDv0::StepManager()
 316 {
 317   //
 318   // Procedure called at every step in the TRD
 319   // Fast simulator. If switched on, a hit is produced when a track
 320   // crosses the border between amplification region and pad plane.
 321   //
 322
 323   Int_t   vol[3];
 324   Int_t   iIdSens, icSens;
 325   Int_t   iIdChamber, icChamber;
 326
 327   Float_t hits[4];
 328
 329   TLorentzVector p;
 330   TClonesArray  &lhits = *fHits;
 331
 332   // Writing out hits enabled?
 333   if (!(fHitsOn)) return;
 334
 335   // Use only charged tracks and count them only once per volume
 336   if (gMC->TrackCharge()    &&
 337       gMC->IsTrackExiting()) {
 338
 339     // Check on sensitive volume
 340     iIdSens = gMC->CurrentVolID(icSens);
 341     if (iIdSens == fIdSens) {
 342
 343       gMC->TrackPosition(p);
 344       for (Int_t i = 0; i < 3; i++) hits[i] = p[i];
 345       // No charge created
 346       hits[3] = 0;
 347
 348       // The sector number
 349       Float_t phi = hits[1] != 0 ? TMath::Atan2(hits[0],hits[1]) : (hits[0] > 0 ? 180. : 0.);
 350       vol[0] = ((Int_t) (phi / 20)) + 1;
 351
 352       // The chamber number
 353       //   1: outer left
 354       //   2: middle left
 355       //   3: inner
 356       //   4: middle right
 357       //   5: outer right
 358       iIdChamber = gMC->CurrentVolOffID(1,icChamber);
 359       if      (iIdChamber == fIdChamber1)
 360         vol[1] = (hits[2] < 0 ? 1 : 5);
 361       else if (iIdChamber == fIdChamber2)
 362         vol[1] = (hits[2] < 0 ? 2 : 4);
 363       else if (iIdChamber == fIdChamber3)
 364         vol[1] = 3;
 365
 366       // The plane number
 367       vol[2] = icChamber - TMath::Nint((Float_t) (icChamber / 7)) * 6;
 368
 369       new(lhits[fNhits++]) AliTRDhit(fIshunt,gAlice->CurrentTrack(),vol,hits);
 370
 371     }
 372
 373   }
 374
 375 }