FMD/AliFMDSimulator.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 // Forward Multiplicity Detector based on Silicon wafers. This class
  21 // contains the base procedures for the Forward Multiplicity detector
  22 // Detector consists of 3 sub-detectors FMD1, FMD2, and FMD3, each of
  23 // which has 1 or 2 rings of silicon sensors.
  24 //
  25 // This is the base class for all FMD manager classes.
  26 //
  27 // The actual code is done by various separate classes.   Below is
  28 // diagram showing the relationship between the various FMD classes
  29 // that handles the simulation
  30 //
  31 //      +--------+ 1     +-----------------+
  32 //      | AliFMD |<>-----| AliFMDSimulator |
  33 //      +--------+       +-----------------+
  34 //                               ^
  35 //                               |
  36 //                 +-------------+-------------+
  37 //                 |                           |
  38 //        +--------------------+   +-------------------+
  39 //        | AliFMDGeoSimulator |   | AliFMDG3Simulator |
  40 //        +--------------------+   +---------+---------+
  41 //
  42 //
  43 // *  AliFMD
  44 //    This defines the interface for the various parts of AliROOT that
  45 //    uses the FMD, like AliFMDSimulator, AliFMDDigitizer,
  46 //    AliFMDReconstructor, and so on.
  47 //
  48 // *  AliFMDSimulator
  49 //    This is the base class for the FMD simulation tasks.   The
  50 //    simulator tasks are responsible to implment the geoemtry, and
  51 //    process hits.
  52 //
  53 // *  AliFMDGeoSimulator
  54 //    This is a concrete implementation of the AliFMDSimulator that
  55 //    uses the TGeo classes directly only.  This defines the active
  56 //    volume as an ONLY XTRU shape with a divided MANY TUBS shape
  57 //    inside to implement the particular shape of the silicon
  58 //    sensors.
  59 //
  60 // *  AliFMDG3Simulator
  61 //    This is a concrete implementation of the AliFMDSimulator that
  62 //    uses the TVirtualMC interface with GEANT 3.21-like messages.
  63 //    This implements the active volume as a divided TUBS shape.  Hits
  64 //    in the corners should be cut away at run time (but currently
  65 //    isn't).
  66 //
  67 #include "AliFMDSimulator.h"    // ALIFMDSIMULATOR_H
  68 #include "AliFMDGeometry.h"     // ALIFMDGEOMETRY_H
  69 #include "AliFMDDetector.h"     // ALIFMDDETECTOR_H
  70 #include "AliFMDRing.h"         // ALIFMDRING_H
  71 #include "AliFMD1.h"            // ALIFMD1_H
  72 #include "AliFMD2.h"            // ALIFMD2_H
  73 #include "AliFMD3.h"            // ALIFMD3_H
  74 #include "AliFMD.h"             // ALIFMD_H
  75 #include <AliRun.h>             // ALIRUN_H
  76 #include <AliMC.h>              // ALIMC_H
  77 #include <AliMagF.h>            // ALIMAGF_H
  78 #include <AliLog.h>             // ALILOG_H
  79 #include <TGeoVolume.h>         // ROOT_TGeoVolume
  80 #include <TGeoTube.h>           // ROOT_TGeoTube
  81 #include <TGeoPcon.h>           // ROOT_TGeoPcon
  82 #include <TGeoMaterial.h>       // ROOT_TGeoMaterial
  83 #include <TGeoMedium.h>         // ROOT_TGeoMedium
  84 #include <TGeoXtru.h>           // ROOT_TGeoXtru
  85 #include <TGeoPolygon.h>        // ROOT_TGeoPolygon
  86 #include <TGeoTube.h>           // ROOT_TGeoTube
  87 #include <TGeoManager.h>        // ROOT_TGeoManager
  88 #include <TTree.h>              // ROOT_TTree
  89 #include <TParticle.h>          // ROOT_TParticle
  90 #include <TLorentzVector.h>     // ROOT_TLorentzVector
  91 #include <TVector2.h>           // ROOT_TVector2
  92 #include <TVector3.h>           // ROOT_TVector3
  93 #include <TVirtualMC.h>         // ROOT_TVirtualMC
  94 #include <TArrayD.h>            // ROOT_TArrayD
  95
  96 //====================================================================
  97 ClassImp(AliFMDSimulator)
  98 #if 0
  99   ; // This is here to keep Emacs for indenting the next line
 100 #endif
 101
 102 //____________________________________________________________________
 103 const Char_t* AliFMDSimulator::fgkActiveName    = "F%cAC";
 104 const Char_t* AliFMDSimulator::fgkSectorName    = "F%cSE";
 105 const Char_t* AliFMDSimulator::fgkStripName     = "F%cST";
 106 const Char_t* AliFMDSimulator::fgkModuleName    = "F%cMO";
 107 const Char_t* AliFMDSimulator::fgkPCBName       = "F%cP%c";
 108 const Char_t* AliFMDSimulator::fgkLongLegName   = "F%cLL";
 109 const Char_t* AliFMDSimulator::fgkShortLegName  = "F%cSL";
 110 const Char_t* AliFMDSimulator::fgkFrontVName    = "F%cFV";
 111 const Char_t* AliFMDSimulator::fgkBackVName     = "F%cBV";
 112 const Char_t* AliFMDSimulator::fgkRingName      = "FMD%c";
 113 const Char_t* AliFMDSimulator::fgkTopHCName     = "F%d%cI";
 114 const Char_t* AliFMDSimulator::fgkBotHCName     = "F%d%cJ";
 115 const Char_t* AliFMDSimulator::fgkTopIHCName    = "F%d%cK";
 116 const Char_t* AliFMDSimulator::fgkBotIHCName    = "F%d%cL";
 117 const Char_t* AliFMDSimulator::fgkNoseName      = "F3SN";
 118 const Char_t* AliFMDSimulator::fgkBackName      = "F3SB";
 119 const Char_t* AliFMDSimulator::fgkBeamName      = "F3SL";
 120 const Char_t* AliFMDSimulator::fgkFlangeName    = "F3SF";
 121
 122 //____________________________________________________________________
 123 AliFMDSimulator::AliFMDSimulator()
 124   : fFMD(0),
 125     fDetailed(kFALSE),
 126     fInnerId(-1),
 127     fOuterId(-1)
 128 {
 129   // Default constructor
 130 }
 131
 132 //____________________________________________________________________
 133 AliFMDSimulator::AliFMDSimulator(AliFMD* fmd, Bool_t detailed)
 134   : TTask("FMDsimulator", "Forward Multiplicity Detector Simulator"),
 135     fFMD(fmd),
 136     fDetailed(detailed),
 137     fInnerId(-1),
 138     fOuterId(-1)
 139 {
 140   // Normal constructor
 141   //
 142   // Parameters:
 143   //
 144   //      fmd           Pointer to AliFMD object
 145   //      detailed      Whether to make a detailed simulation or not
 146   //
 147 }
 148
 149
 150 //____________________________________________________________________
 151 void
 152 AliFMDSimulator::DefineMaterials()
 153 {
 154   // Define the materials and tracking mediums needed by the FMD
 155   // simulation.   These mediums are made by sending the messages
 156   // AliMaterial, AliMixture, and AliMedium to the passed AliModule
 157   // object module.   The defined mediums are
 158   //
 159   //    FMD Si$         Silicon (active medium in sensors)
 160   //    FMD C$          Carbon fibre (support cone for FMD3 and vacuum pipe)
 161   //    FMD Al$         Aluminium (honeycomb support plates)
 162   //    FMD PCB$        Printed Circuit Board (FEE board with VA1_ALICE)
 163   //    FMD Chip$       Electronics chips (currently not used)
 164   //    FMD Air$        Air (Air in the FMD)
 165   //    FMD Plastic$    Plastic (Support legs for the hybrid cards)
 166   //
 167   // Pointers to TGeoMedium objects are retrived from the TGeoManager
 168   // singleton.  These pointers are later used when setting up the
 169   // geometry
 170   AliDebug(10, "\tCreating materials");
 171   // Get pointer to geometry singleton object.
 172   AliFMDGeometry* geometry = AliFMDGeometry::Instance();
 173   geometry->Init();
 174
 175   Int_t    id;
 176   Double_t a                = 0;
 177   Double_t z                = 0;
 178   Double_t density          = 0;
 179   Double_t radiationLength  = 0;
 180   Double_t absorbtionLength = 999;
 181   Int_t    fieldType        = gAlice->Field()->Integ();     // Field type
 182   Double_t maxField         = gAlice->Field()->Max();     // Field max.
 183   Double_t maxBending       = 0;     // Max Angle
 184   Double_t maxStepSize      = 0.001; // Max step size
 185   Double_t maxEnergyLoss    = 1;     // Max Delta E
 186   Double_t precision        = 0.001; // Precision
 187   Double_t minStepSize      = 0.001; // Minimum step size
 188
 189   // Silicon
 190   a                = 28.0855;
 191   z                = 14.;
 192   density          = geometry->GetSiDensity();
 193   radiationLength  = 9.36;
 194   maxBending       = 1;
 195   maxStepSize      = .001;
 196   precision        = .001;
 197   minStepSize      = .001;
 198   id               = kSiId;
 199   fFMD->AliMaterial(id, "Si$",
 200                       a, z, density, radiationLength, absorbtionLength);
 201   fFMD->AliMedium(kSiId, "Si$",
 202                     id,1,fieldType,maxField,maxBending,
 203                     maxStepSize,maxEnergyLoss,precision,minStepSize);
 204
 205
 206   // Carbon
 207   a                = 12.011;
 208   z                = 6.;
 209   density          = 2.265;
 210   radiationLength  = 18.8;
 211   maxBending       = 10;
 212   maxStepSize      = .01;
 213   precision        = .003;
 214   minStepSize      = .003;
 215   id               = kCarbonId;
 216   fFMD->AliMaterial(id, "Carbon$",
 217                       a, z, density, radiationLength, absorbtionLength);
 218   fFMD->AliMedium(kCarbonId, "Carbon$",
 219                     id,0,fieldType,maxField,maxBending,
 220                     maxStepSize,maxEnergyLoss,precision,minStepSize);
 221
 222   // Aluminum
 223   a                = 26.981539;
 224   z                = 13.;
 225   density          = 2.7;
 226   radiationLength  = 8.9;
 227   id               = kAlId;
 228   fFMD->AliMaterial(id, "Aluminum$",
 229                       a, z, density, radiationLength, absorbtionLength);
 230   fFMD->AliMedium(kAlId, "Aluminum$",
 231                     id, 0, fieldType, maxField, maxBending,
 232                     maxStepSize, maxEnergyLoss, precision, minStepSize);
 233
 234
 235   // Silicon chip
 236   {
 237     Float_t as[] = { 12.0107,      14.0067,      15.9994,
 238                      1.00794,      28.0855,     107.8682 };
 239     Float_t zs[] = {  6.,           7.,           8.,
 240                       1.,          14.,          47. };
 241     Float_t ws[] = {  0.039730642,  0.001396798,  0.01169634,
 242                       0.004367771,  0.844665,     0.09814344903 };
 243     density = 2.36436;
 244     maxBending       = 10;
 245     maxStepSize      = .01;
 246     precision        = .003;
 247     minStepSize      = .003;
 248     id = kSiChipId;
 249     fFMD->AliMixture(id, "Si Chip$", as, zs, density, 6, ws);
 250     fFMD->AliMedium(kSiChipId, "Si Chip$",
 251                       id, 0, fieldType, maxField, maxBending,
 252                       maxStepSize, maxEnergyLoss, precision, minStepSize);
 253   }
 254
 255 #if 0
 256   // Kaption
 257   {
 258     Float_t as[] = { 1.00794,  12.0107,  14.010,   15.9994};
 259     Float_t zs[] = { 1.,        6.,       7.,       8.};
 260     Float_t ws[] = { 0.026362,  0.69113,  0.07327,  0.209235};
 261     density          = 1.42;
 262     maxBending       = 1;
 263     maxStepSize      = .001;
 264     precision        = .001;
 265     minStepSize      = .001;
 266     id               = KaptionId;
 267     fFMD->AliMixture(id, "Kaption$", as, zs, density, 4, ws);
 268     fFMD->AliMedium(kAlId, "Kaption$",
 269                       id,0,fieldType,maxField,maxBending,
 270                       maxStepSize,maxEnergyLoss,precision,minStepSize);
 271   }
 272 #endif
 273
 274   // Air
 275   {
 276     Float_t as[] = { 12.0107, 14.0067,   15.9994,  39.948 };
 277     Float_t zs[] = {  6.,      7.,       8.,       18. };
 278     Float_t ws[] = { 0.000124, 0.755267, 0.231781, 0.012827 };
 279     density      = .00120479;
 280     maxBending   = 1;
 281     maxStepSize  = .001;
 282     precision    = .001;
 283     minStepSize  = .001;
 284     id           = kAirId;
 285     fFMD->AliMixture(id, "Air$", as, zs, density, 4, ws);
 286     fFMD->AliMedium(kAirId, "Air$",
 287                       id,0,fieldType,maxField,maxBending,
 288                       maxStepSize,maxEnergyLoss,precision,minStepSize);
 289   }
 290
 291   // PCB
 292   {
 293     Float_t zs[] = { 14.,         20.,         13.,         12.,
 294                       5.,         22.,         11.,         19.,
 295                      26.,          9.,          8.,          6.,
 296                       7.,          1.};
 297     Float_t as[] = { 28.0855,     40.078,      26.981538,   24.305,
 298                      10.811,      47.867,      22.98977,    39.0983,
 299                      55.845,      18.9984,     15.9994,     12.0107,
 300                      14.0067,      1.00794};
 301     Float_t ws[] = {  0.15144894,  0.08147477,  0.04128158,  0.00904554,
 302                       0.01397570,  0.00287685,  0.00445114,  0.00498089,
 303                       0.00209828,  0.00420000,  0.36043788,  0.27529426,
 304                       0.01415852,  0.03427566};
 305     density      = 1.8;
 306     maxBending   = 1;
 307     maxStepSize  = .001;
 308     precision    = .001;
 309     minStepSize  = .001;
 310     id           = kPcbId;
 311     fFMD->AliMixture(id, "PCB$", as, zs, density, 14, ws);
 312     fFMD->AliMedium(kPcbId, "PCB$",
 313                       id,0,fieldType,maxField,maxBending,
 314                       maxStepSize,maxEnergyLoss,precision,minStepSize);
 315   }
 316
 317   // Plastic
 318   {
 319     Float_t as[] = { 1.01, 12.01 };
 320     Float_t zs[] = { 1.,   6.    };
 321     Float_t ws[] = { 1.,   1.    };
 322     density      = 1.03;
 323     maxBending   = 10;
 324     maxStepSize  = .01;
 325     precision    = .003;
 326     minStepSize  = .003;
 327     id           = kPlasticId;
 328     fFMD->AliMixture(id, "Plastic$", as, zs, density, -2, ws);
 329     fFMD->AliMedium(kPlasticId, "Plastic$",
 330                       id,0,fieldType,maxField,maxBending,
 331                       maxStepSize,maxEnergyLoss,precision,minStepSize);
 332   }
 333 }
 334
 335 //____________________________________________________________________
 336 void
 337 AliFMDSimulator::Exec(Option_t* /* option */)
 338 {
 339   // Member function that is executed each time a hit is made in the
 340   // FMD.  None-charged particles are ignored.   Dead tracks  are
 341   // ignored.
 342   //
 343   // The procedure is as follows:
 344   //
 345   //   - IF NOT track is alive THEN RETURN ENDIF
 346   //   - IF NOT particle is charged THEN RETURN ENDIF
 347   //   - IF NOT volume name is "STRI" or "STRO" THEN RETURN ENDIF
 348   //   - Get strip number (volume copy # minus 1)
 349   //   - Get phi division number (mother volume copy #)
 350   //   - Get module number (grand-mother volume copy #)
 351   //   - section # = 2 * module # + phi division # - 1
 352   //   - Get ring Id from volume name
 353   //   - Get detector # from grand-grand-grand-mother volume name
 354   //   - Get pointer to sub-detector object.
 355   //   - Get track position
 356   //   - IF track is entering volume AND track is inside real shape THEN
 357   //   -   Reset energy deposited
 358   //   -   Get track momentum
 359   //   -   Get particle ID #
 360   ///  - ENDIF
 361   //   - IF track is inside volume AND inside real shape THEN
 362   ///  -   Update energy deposited
 363   //   - ENDIF
 364   //   - IF track is inside real shape AND (track is leaving volume,
 365   //         or it died, or it is stopped  THEN
 366   //   -   Create a hit
 367   //   - ENDIF
 368   //
 369   TVirtualMC* mc = TVirtualMC::GetMC();
 370   static int nCall = 0;
 371
 372   if (!mc->IsTrackAlive()) return;
 373   if (TMath::Abs(mc->TrackCharge()) <= 0) return;
 374
 375   Int_t copy;
 376   Int_t vol = mc->CurrentVolID(copy);
 377   if (vol != fInnerId && vol != fOuterId) {
 378     AliDebug(25, Form("Not an FMD volume %d '%s' (%d or %d)",
 379                       vol, mc->CurrentVolName(), fInnerId, fOuterId));
 380     return;
 381   }
 382
 383   nCall++;
 384   // Check that the track is actually within the active area
 385   Bool_t entering = mc->IsTrackEntering();
 386   Bool_t inside   = mc->IsTrackInside();
 387   Bool_t out      = (mc->IsTrackExiting()|| mc->IsTrackDisappeared()||
 388                      mc->IsTrackStop());
 389
 390   // Reset the energy deposition for this track, and update some of
 391   // our parameters.
 392   if (entering) {
 393     AliDebug(15, Form("Track # %8d entering active FMD volume %s: "
 394                       "Edep=%f (call # %d)",
 395                       gAlice->GetMCApp()->GetCurrentTrackNumber(),
 396                       mc->CurrentVolPath(), 1000 * mc->Edep(), nCall));
 397     fCurrentDeltaE = 0;
 398
 399     // Get production vertex and momentum of the track
 400     mc->TrackMomentum(fCurrentP);
 401     mc->TrackPosition(fCurrentV);
 402     fCurrentPdg = mc->IdFromPDG(mc->TrackPid());
 403   }
 404
 405   // If the track is inside, then update the energy deposition
 406   if (inside && fCurrentDeltaE >= 0) {
 407     fCurrentDeltaE += 1000 * mc->Edep();
 408     AliDebug(15, Form("Track # %8d inside active FMD volume %s: Edep=%f, "
 409                       "Accumulated Edep=%f (call # %d)",
 410                       gAlice->GetMCApp()->GetCurrentTrackNumber(),
 411                       mc->CurrentVolPath(), 1000 * mc->Edep(),
 412                       fCurrentDeltaE, nCall));
 413   }
 414
 415   // The track exits the volume, or it disappeared in the volume, or
 416   // the track is stopped because it no longer fulfills the cuts
 417   // defined, then we create a hit.
 418   if (out && fCurrentDeltaE >= 0) {
 419
 420     Int_t strip = copy - 1;
 421     Int_t sectordiv;
 422     mc->CurrentVolOffID(fSectorOff, sectordiv);
 423     Int_t module;
 424     mc->CurrentVolOffID(fModuleOff, module);
 425     Int_t sector = 2 * module + sectordiv;
 426     Int_t iring;
 427     mc->CurrentVolOffID(fRingOff, iring);
 428     Char_t ring = Char_t(iring);
 429     Int_t detector;
 430     mc->CurrentVolOffID(fDetectorOff, detector);
 431
 432
 433     AliFMDGeometry* fmd = AliFMDGeometry::Instance();
 434     Double_t rz = fmd->GetDetector(detector)->GetRingZ(ring);
 435     Int_t    n  = fmd->GetDetector(detector)->GetRing(ring)->GetNSectors();
 436     if (rz < 0) {
 437       Int_t s = ((n - sector + n / 2) % n) + 1;
 438       AliDebug(40, Form("Recalculating sector to %d (=%d-%d+%d/2%%%d+1 z=%f)",
 439                         s, n, sector, n, n, rz));
 440       sector = s;
 441     }
 442     if (sector < 1 || sector > n) {
 443       Warning("Step", "sector # %d out of range (0-%d)", sector-1, n-1);
 444       return;
 445     }
 446     sector--;
 447     fCurrentDeltaE += 1000 * mc->Edep();
 448     AliDebug(15, Form("Track # %8d leaving active FMD volume %s: Edep=%f, "
 449                       "Accumulated Edep=%f (call # %d nCall)",
 450                       gAlice->GetMCApp()->GetCurrentTrackNumber(),
 451                       mc->CurrentVolPath(),
 452                       1000 * mc->Edep(), fCurrentDeltaE, nCall));
 453
 454     AliDebug(20, Form("Processing hit in FMD%d%c[%2d,%3d]: %f",
 455                       detector, ring, sector, strip, fCurrentDeltaE));
 456
 457     fFMD->AddHitByFields(gAlice->GetMCApp()->GetCurrentTrackNumber(),
 458                          UShort_t(detector), ring, UShort_t(sector),
 459                          UShort_t(strip),
 460                          fCurrentV.X(), fCurrentV.Y(), fCurrentV.Z(),
 461                          fCurrentP.X(), fCurrentP.Y(), fCurrentP.Z(),
 462                          fCurrentDeltaE, fCurrentPdg, fCurrentV.T());
 463     fCurrentDeltaE = -1;
 464   }
 465 }
 466
 467
 468
 469 //____________________________________________________________________
 470 //
 471 // EOF
 472 //