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