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 //                                 +----------------------+
  44 //                                 | AliFMDG3OldSimulator |
  45 //                                 +----------------------+
  46 //
  47 // *  AliFMD
  48 //    This defines the interface for the various parts of AliROOT that
  49 //    uses the FMD, like AliFMDSimulator, AliFMDDigitizer,
  50 //    AliFMDReconstructor, and so on.
  51 //
  52 // *  AliFMDSimulator
  53 //    This is the base class for the FMD simulation tasks.   The
  54 //    simulator tasks are responsible to implment the geoemtry, and
  55 //    process hits.
  56 //
  57 // *  AliFMDGeoSimulator
  58 //    This is a concrete implementation of the AliFMDSimulator that
  59 //    uses the TGeo classes directly only.  This defines the active
  60 //    volume as an ONLY XTRU shape with a divided MANY TUBS shape
  61 //    inside to implement the particular shape of the silicon
  62 //    sensors.
  63 //
  64 // *  AliFMDG3Simulator
  65 //    This is a concrete implementation of the AliFMDSimulator that
  66 //    uses the TVirtualMC interface with GEANT 3.21-like messages.
  67 //    This implements the active volume as a divided TUBS shape.  Hits
  68 //    in the corners should be cut away at run time (but currently
  69 //    isn't).
  70 //
  71 // *  AliFMDG3OldSimulator
  72 //    This is a concrete implementation of AliFMDSimulator.   It
  73 //    approximates the of the rings as segmented disks.
  74 //
  75 #include "AliFMDSimulator.h"    // ALIFMDSIMULATOR_H
  76 #include "AliFMDGeometry.h"     // ALIFMDGEOMETRY_H
  77 #include "AliFMDDetector.h"     // ALIFMDDETECTOR_H
  78 #include "AliFMDRing.h"         // ALIFMDRING_H
  79 #include "AliFMD1.h"            // ALIFMD1_H
  80 #include "AliFMD2.h"            // ALIFMD2_H
  81 #include "AliFMD3.h"            // ALIFMD3_H
  82 #include "AliFMD.h"             // ALIFMD_H
  83 #include "AliFMDHit.h"          // ALIFMDHIT_H
  84 #include <AliRun.h>             // ALIRUN_H
  85 #include <AliMC.h>              // ALIMC_H
  86 #include <AliMagF.h>            // ALIMAGF_H
  87 #include <AliLog.h>             // ALILOG_H
  88 #include <TGeoVolume.h>         // ROOT_TGeoVolume
  89 #include <TGeoTube.h>           // ROOT_TGeoTube
  90 #include <TGeoPcon.h>           // ROOT_TGeoPcon
  91 #include <TGeoMaterial.h>       // ROOT_TGeoMaterial
  92 #include <TGeoMedium.h>         // ROOT_TGeoMedium
  93 #include <TGeoXtru.h>           // ROOT_TGeoXtru
  94 #include <TGeoPolygon.h>        // ROOT_TGeoPolygon
  95 #include <TGeoTube.h>           // ROOT_TGeoTube
  96 #include <TGeoManager.h>        // ROOT_TGeoManager
  97 #include <TTree.h>              // ROOT_TTree
  98 #include <TParticle.h>          // ROOT_TParticle
  99 #include <TLorentzVector.h>     // ROOT_TLorentzVector
 100 #include <TVector2.h>           // ROOT_TVector2
 101 #include <TVector3.h>           // ROOT_TVector3
 102 #include <TVirtualMC.h>         // ROOT_TVirtualMC
 103 #include <TArrayD.h>            // ROOT_TArrayD
 104
 105
 106 //====================================================================
 107 ClassImp(AliFMDSimulator)
 108 #if 0
 109   ; // This is here to keep Emacs for indenting the next line
 110 #endif
 111
 112 //____________________________________________________________________
 113 const Char_t* AliFMDSimulator::fgkActiveName    = "F%cAC";
 114 const Char_t* AliFMDSimulator::fgkSectorName    = "F%cSE";
 115 const Char_t* AliFMDSimulator::fgkStripName     = "F%cST";
 116 const Char_t* AliFMDSimulator::fgkModuleName    = "F%cMO";
 117 const Char_t* AliFMDSimulator::fgkPCBName       = "F%cP%c";
 118 const Char_t* AliFMDSimulator::fgkLongLegName   = "F%cLL";
 119 const Char_t* AliFMDSimulator::fgkShortLegName  = "F%cSL";
 120 const Char_t* AliFMDSimulator::fgkFrontVName    = "F%cFV";
 121 const Char_t* AliFMDSimulator::fgkBackVName     = "F%cBV";
 122 const Char_t* AliFMDSimulator::fgkRingName      = "FMD%c";
 123 const Char_t* AliFMDSimulator::fgkTopHCName     = "F%d%cI";
 124 const Char_t* AliFMDSimulator::fgkBotHCName     = "F%d%cJ";
 125 const Char_t* AliFMDSimulator::fgkTopIHCName    = "F%d%cK";
 126 const Char_t* AliFMDSimulator::fgkBotIHCName    = "F%d%cL";
 127 const Char_t* AliFMDSimulator::fgkNoseName      = "F3SN";
 128 const Char_t* AliFMDSimulator::fgkBackName      = "F3SB";
 129 const Char_t* AliFMDSimulator::fgkBeamName      = "F3SL";
 130 const Char_t* AliFMDSimulator::fgkFlangeName    = "F3SF";
 131
 132 //____________________________________________________________________
 133 AliFMDSimulator::AliFMDSimulator()
 134   : fFMD(0),
 135     fActiveId(4),
 136     fDetailed(kFALSE),
 137     fUseDivided(kFALSE),
 138     fUseAssembly(kTRUE),
 139     fBad(0)
 140 {
 141   // Default constructor
 142 }
 143
 144 //____________________________________________________________________
 145 AliFMDSimulator::AliFMDSimulator(AliFMD* fmd, Bool_t detailed)
 146   : TTask("FMDSimulator", "Forward Multiplicity Detector Simulator"),
 147     fFMD(fmd),
 148     fActiveId(4),
 149     fDetailed(detailed),
 150     fUseDivided(kFALSE),
 151     fUseAssembly(kTRUE),
 152     fBad(0)
 153 {
 154   // Normal constructor
 155   //
 156   // Parameters:
 157   //
 158   //      fmd           Pointer to AliFMD object
 159   //      detailed      Whether to make a detailed simulation or not
 160   //
 161   fBad = new TClonesArray("AliFMDHit");
 162 }
 163
 164
 165 //____________________________________________________________________
 166 void
 167 AliFMDSimulator::DefineMaterials()
 168 {
 169   // Define the materials and tracking mediums needed by the FMD
 170   // simulation.   These mediums are made by sending the messages
 171   // AliMaterial, AliMixture, and AliMedium to the passed AliModule
 172   // object module.   The defined mediums are
 173   //
 174   //    FMD Si$         Silicon (active medium in sensors)
 175   //    FMD C$          Carbon fibre (support cone for FMD3 and vacuum pipe)
 176   //    FMD Al$         Aluminium (honeycomb support plates)
 177   //    FMD PCB$        Printed Circuit Board (FEE board with VA1_3)
 178   //    FMD Chip$       Electronics chips (currently not used)
 179   //    FMD Air$        Air (Air in the FMD)
 180   //    FMD Plastic$    Plastic (Support legs for the hybrid cards)
 181   //
 182   // Pointers to TGeoMedium objects are retrived from the TGeoManager
 183   // singleton.  These pointers are later used when setting up the
 184   // geometry
 185   AliDebug(10, "\tCreating materials");
 186   AliDebug(1,  Form("\tGeometry options: %s, %s, %s",
 187                     (fDetailed    ? "detailed" : "coarse"),
 188                     (fUseDivided  ? "divided into strips" : "one volume"),
 189                     (fUseAssembly ? "within assemblies" : "in real volumes")));
 190   // Get pointer to geometry singleton object.
 191   AliFMDGeometry* geometry = AliFMDGeometry::Instance();
 192   geometry->Init();
 193
 194   Int_t    id;
 195   Double_t a                = 0;
 196   Double_t z                = 0;
 197   Double_t density          = 0;
 198   Double_t radiationLength  = 0;
 199   Double_t absorbtionLength = 999;
 200   Int_t    fieldType        = gAlice->Field()->Integ();     // Field type
 201   Double_t maxField         = gAlice->Field()->Max();     // Field max.
 202   Double_t maxBending       = 0;     // Max Angle
 203   Double_t maxStepSize      = 0.001; // Max step size
 204   Double_t maxEnergyLoss    = 1;     // Max Delta E
 205   Double_t precision        = 0.001; // Precision
 206   Double_t minStepSize      = 0.001; // Minimum step size
 207
 208   // Silicon
 209   a                = 28.0855;
 210   z                = 14.;
 211   density          = geometry->GetSiDensity();
 212   radiationLength  = 9.36;
 213   maxBending       = 1;
 214   maxStepSize      = .001;
 215   precision        = .001;
 216   minStepSize      = .001;
 217   id               = kSiId;
 218   fFMD->AliMaterial(id, "Si$",
 219                       a, z, density, radiationLength, absorbtionLength);
 220   fFMD->AliMedium(kSiId, "Si$",
 221                     id,1,fieldType,maxField,maxBending,
 222                     maxStepSize,maxEnergyLoss,precision,minStepSize);
 223
 224
 225   // Carbon
 226   a                = 12.011;
 227   z                = 6.;
 228   density          = 2.265;
 229   radiationLength  = 18.8;
 230   maxBending       = 10;
 231   maxStepSize      = .01;
 232   precision        = .003;
 233   minStepSize      = .003;
 234   id               = kCarbonId;
 235   fFMD->AliMaterial(id, "Carbon$",
 236                       a, z, density, radiationLength, absorbtionLength);
 237   fFMD->AliMedium(kCarbonId, "Carbon$",
 238                     id,0,fieldType,maxField,maxBending,
 239                     maxStepSize,maxEnergyLoss,precision,minStepSize);
 240
 241   // Aluminum
 242   a                = 26.981539;
 243   z                = 13.;
 244   density          = 2.7;
 245   radiationLength  = 8.9;
 246   id               = kAlId;
 247   fFMD->AliMaterial(id, "Aluminum$",
 248                       a, z, density, radiationLength, absorbtionLength);
 249   fFMD->AliMedium(kAlId, "Aluminum$",
 250                     id, 0, fieldType, maxField, maxBending,
 251                     maxStepSize, maxEnergyLoss, precision, minStepSize);
 252
 253
 254   // Copper
 255   a                = 63.546;
 256   z                = 29;
 257   density          =  8.96;
 258   radiationLength  =  1.43;
 259   id               = kCopperId;
 260   fFMD->AliMaterial(id, "Copper$",
 261                       a, z, density, radiationLength, absorbtionLength);
 262   fFMD->AliMedium(kCopperId, "Copper$",
 263                     id, 0, fieldType, maxField, maxBending,
 264                     maxStepSize, maxEnergyLoss, precision, minStepSize);
 265
 266
 267   // Silicon chip
 268   {
 269     Float_t as[] = { 12.0107,      14.0067,      15.9994,
 270                       1.00794,     28.0855,     107.8682 };
 271     Float_t zs[] = {  6.,           7.,           8.,
 272                       1.,          14.,          47. };
 273     Float_t ws[] = {  0.039730642,  0.001396798,  0.01169634,
 274                       0.004367771,  0.844665,     0.09814344903 };
 275     density          = 2.36436;
 276     maxBending       = 10;
 277     maxStepSize      = .01;
 278     precision        = .003;
 279     minStepSize      = .003;
 280     id               = kSiChipId;
 281     fFMD->AliMixture(id, "Si Chip$", as, zs, density, 6, ws);
 282     fFMD->AliMedium(kSiChipId, "Si Chip$",
 283                       id, 0, fieldType, maxField, maxBending,
 284                       maxStepSize, maxEnergyLoss, precision, minStepSize);
 285   }
 286
 287   // Kaption
 288   {
 289     Float_t as[] = { 1.00794,  12.0107,  14.010,   15.9994};
 290     Float_t zs[] = { 1.,        6.,       7.,       8.};
 291     Float_t ws[] = { 0.026362,  0.69113,  0.07327,  0.209235};
 292     density          = 1.42;
 293     maxBending       = 1;
 294     maxStepSize      = .001;
 295     precision        = .001;
 296     minStepSize      = .001;
 297     id               = kKaptonId;
 298     fFMD->AliMixture(id, "Kaption$", as, zs, density, 4, ws);
 299     fFMD->AliMedium(kKaptonId, "Kaption$",
 300                       id,0,fieldType,maxField,maxBending,
 301                       maxStepSize,maxEnergyLoss,precision,minStepSize);
 302   }
 303
 304   // Air
 305   {
 306     Float_t as[] = { 12.0107, 14.0067,   15.9994,  39.948 };
 307     Float_t zs[] = {  6.,      7.,       8.,       18. };
 308     Float_t ws[] = { 0.000124, 0.755267, 0.231781, 0.012827 };
 309     density      = .00120479;
 310     maxBending   = 1;
 311     maxStepSize  = .001;
 312     precision    = .001;
 313     minStepSize  = .001;
 314     id           = kAirId;
 315     fFMD->AliMixture(id, "Air$", as, zs, density, 4, ws);
 316     fFMD->AliMedium(kAirId, "Air$",
 317                       id,0,fieldType,maxField,maxBending,
 318                       maxStepSize,maxEnergyLoss,precision,minStepSize);
 319   }
 320
 321   // PCB
 322   {
 323     Float_t zs[] = { 14.,         20.,         13.,         12.,
 324                       5.,         22.,         11.,         19.,
 325                      26.,          9.,          8.,          6.,
 326                       7.,          1.};
 327     Float_t as[] = { 28.0855,     40.078,      26.981538,   24.305,
 328                      10.811,      47.867,      22.98977,    39.0983,
 329                      55.845,      18.9984,     15.9994,     12.0107,
 330                      14.0067,      1.00794};
 331     Float_t ws[] = {  0.15144894,  0.08147477,  0.04128158,  0.00904554,
 332                       0.01397570,  0.00287685,  0.00445114,  0.00498089,
 333                       0.00209828,  0.00420000,  0.36043788,  0.27529426,
 334                       0.01415852,  0.03427566};
 335     density      = 1.8;
 336     maxBending   = 1;
 337     maxStepSize  = .001;
 338     precision    = .001;
 339     minStepSize  = .001;
 340     id           = kPcbId;
 341     fFMD->AliMixture(id, "PCB$", as, zs, density, 14, ws);
 342     fFMD->AliMedium(kPcbId, "PCB$",
 343                       id,0,fieldType,maxField,maxBending,
 344                       maxStepSize,maxEnergyLoss,precision,minStepSize);
 345   }
 346
 347   // Plastic
 348   {
 349     Float_t as[] = { 1.01, 12.01 };
 350     Float_t zs[] = { 1.,   6.    };
 351     Float_t ws[] = { 1.,   1.    };
 352     density      = 1.03;
 353     maxBending   = 10;
 354     maxStepSize  = .01;
 355     precision    = .003;
 356     minStepSize  = .003;
 357     id           = kPlasticId;
 358     fFMD->AliMixture(id, "Plastic$", as, zs, density, -2, ws);
 359     fFMD->AliMedium(kPlasticId, "Plastic$",
 360                       id,0,fieldType,maxField,maxBending,
 361                       maxStepSize,maxEnergyLoss,precision,minStepSize);
 362   }
 363 }
 364
 365 //____________________________________________________________________
 366 Bool_t
 367 AliFMDSimulator::IsActive(Int_t volId) const
 368 {
 369   for (Int_t i = 0; i < fActiveId.fN; i++)
 370     if  (volId == fActiveId[i]) return kTRUE;
 371   return kFALSE;
 372 }
 373
 374 //____________________________________________________________________
 375 Bool_t
 376 AliFMDSimulator::VMC2FMD(TLorentzVector& v, UShort_t& detector,
 377                          Char_t& ring, UShort_t& sector, UShort_t& strip)
 378 {
 379   TVirtualMC* mc = TVirtualMC::GetMC();
 380
 381   // Get track position
 382   mc->TrackPosition(v);
 383   Int_t moduleno; mc->CurrentVolOffID(fModuleOff, moduleno);
 384   Int_t iring;    mc->CurrentVolOffID(fRingOff, iring);   ring = Char_t(iring);
 385   Int_t det;      mc->CurrentVolOffID(fDetectorOff, det); detector = det;
 386
 387
 388   // Get the ring geometry
 389   AliFMDGeometry*  fmd = AliFMDGeometry::Instance();
 390   //Int_t     nsec = fmd->GetDetector(detector)->GetRing(ring)->GetNSectors();
 391   Int_t     nstr = fmd->GetDetector(detector)->GetRing(ring)->GetNStrips();
 392   Double_t  lowr = fmd->GetDetector(detector)->GetRing(ring)->GetLowR();
 393   Double_t  highr= fmd->GetDetector(detector)->GetRing(ring)->GetHighR();
 394   Double_t  theta= fmd->GetDetector(detector)->GetRing(ring)->GetTheta();
 395
 396   // Figure out the strip number
 397   Double_t r     = TMath::Sqrt(v.X() * v.X() + v.Y() * v.Y());
 398   Double_t pitch = (highr - lowr) / nstr;
 399   Int_t    str   = Int_t((r - lowr) / pitch);
 400   if (str < 0 || str >= nstr) return kFALSE;
 401   strip          = str;
 402
 403   // Figure out the sector number
 404   Double_t phi    = TMath::ATan2(v.Y(), v.X()) * 180. / TMath::Pi();
 405   if (phi < 0) phi = 360. + phi;
 406   Double_t t      = phi - 2 * moduleno * theta;
 407   sector          = 2 * moduleno;
 408   if (t < 0 || t > 2 * theta) return kFALSE;
 409   else if (t > theta)         sector += 1;
 410
 411   AliDebug(40, Form("<1> Inside an active FMD volume FMD%d%c[%2d,%3d] %s",
 412                     detector, ring, sector, strip, mc->CurrentVolPath()));
 413   return kTRUE;
 414 }
 415
 416 //____________________________________________________________________
 417 Bool_t
 418 AliFMDSimulator::VMC2FMD(Int_t copy, TLorentzVector& v,
 419                          UShort_t& detector, Char_t& ring,
 420                          UShort_t& sector, UShort_t& strip)
 421 {
 422   TVirtualMC* mc = TVirtualMC::GetMC();
 423
 424   strip = copy - 1;
 425   Int_t sectordiv; mc->CurrentVolOffID(fSectorOff, sectordiv);
 426   if (fModuleOff >= 0) {
 427     Int_t module;    mc->CurrentVolOffID(fModuleOff, module);
 428     sector = 2 * module + sectordiv;
 429   }
 430   else
 431     sector = sectordiv;
 432   Int_t iring;     mc->CurrentVolOffID(fRingOff, iring); ring = Char_t(iring);
 433   Int_t det;       mc->CurrentVolOffID(fDetectorOff, det); detector = det;
 434
 435   AliFMDGeometry* fmd = AliFMDGeometry::Instance();
 436   //Double_t  rz  = fmd->GetDetector(detector)->GetRingZ(ring);
 437   Int_t     n   = fmd->GetDetector(detector)->GetRing(ring)->GetNSectors();
 438 #if 0
 439   if (rz < 0) {
 440     Int_t s = ((n - sector + n / 2) % n) + 1;
 441     AliDebug(1, Form("Recalculating sector to %d (=%d-%d+%d/2%%%d+1 z=%f)",
 442                      s, n, sector, n, n, rz));
 443     sector = s;
 444   }
 445 #endif
 446   if (sector < 1 || sector > n) {
 447     Warning("Step", "sector # %d out of range (0-%d)", sector-1, n-1);
 448     return kFALSE;
 449   }
 450   sector--;
 451   // Get track position
 452   mc->TrackPosition(v);
 453   AliDebug(15, Form("<2> Inside an active FMD volume FMD%d%c[%2d,%3d] %s",
 454                     detector, ring, sector, strip, mc->CurrentVolPath()));
 455
 456   return kTRUE;
 457 }
 458
 459 //____________________________________________________________________
 460 void
 461 AliFMDSimulator::Exec(Option_t* /* option */)
 462 {
 463   // Member function that is executed each time a hit is made in the
 464   // FMD.  None-charged particles are ignored.   Dead tracks  are
 465   // ignored.
 466   //
 467   // The procedure is as follows:
 468   //
 469   //   - IF NOT track is alive THEN RETURN ENDIF
 470   //   - IF NOT particle is charged THEN RETURN ENDIF
 471   //   - IF NOT volume name is "STRI" or "STRO" THEN RETURN ENDIF
 472   //   - Get strip number (volume copy # minus 1)
 473   //   - Get phi division number (mother volume copy #)
 474   //   - Get module number (grand-mother volume copy #)
 475   //   - section # = 2 * module # + phi division # - 1
 476   //   - Get ring Id from volume name
 477   //   - Get detector # from grand-grand-grand-mother volume name
 478   //   - Get pointer to sub-detector object.
 479   //   - Get track position
 480   //   - IF track is entering volume AND track is inside real shape THEN
 481   //   -   Reset energy deposited
 482   //   -   Get track momentum
 483   //   -   Get particle ID #
 484   ///  - ENDIF
 485   //   - IF track is inside volume AND inside real shape THEN
 486   ///  -   Update energy deposited
 487   //   - ENDIF
 488   //   - IF track is inside real shape AND (track is leaving volume,
 489   //         or it died, or it is stopped  THEN
 490   //   -   Create a hit
 491   //   - ENDIF
 492   //
 493   TVirtualMC* mc = TVirtualMC::GetMC();
 494   if (!mc->IsTrackAlive()) return;
 495   Double_t absQ = TMath::Abs(mc->TrackCharge());
 496   if (absQ <= 0) return;
 497
 498   Int_t copy;
 499   Int_t vol = mc->CurrentVolID(copy);
 500   if (!IsActive(vol)) {
 501     AliDebug(50, Form("Not an FMD volume %d '%s'",vol,mc->CurrentVolName()));
 502     return;
 503   }
 504   TLorentzVector v;
 505   UShort_t       detector;
 506   Char_t         ring;
 507   UShort_t       sector;
 508   UShort_t       strip;
 509
 510   if (fUseDivided) {
 511     if (!VMC2FMD(copy, v, detector, ring, sector, strip)) return;
 512   } else {
 513     if (!VMC2FMD(v, detector, ring, sector, strip)) return;
 514   }
 515   TLorentzVector p;
 516   mc->TrackMomentum(p);
 517   Int_t    trackno = gAlice->GetMCApp()->GetCurrentTrackNumber();
 518   Int_t    pdg     = mc->TrackPid();
 519   Double_t mass    = mc->TrackMass();
 520   Double_t edep    = mc->Edep() * 1000; // keV
 521   Double_t poverm  = (mass == 0 ? 0 : p.P() / mass);
 522   Bool_t   isBad   = kFALSE;
 523
 524   // This `if' is to debug abnormal energy depositions.  We trigger on
 525   // p/m approx larger than or equal to a MIP, and a large edep - more
 526   // than 1 keV - a MIP is 100 eV.
 527   if (edep > absQ * absQ && poverm > 1) {
 528     isBad = kTRUE;
 529     TArrayI procs;
 530     mc->StepProcesses(procs);
 531     TString processes;
 532     for (Int_t ip = 0; ip < procs.fN; ip++) {
 533       if (ip != 0) processes.Append(",");
 534       processes.Append(TMCProcessName[procs.fArray[ip]]);
 535     }
 536     TDatabasePDG* pdgDB        = TDatabasePDG::Instance();
 537     TParticlePDG* particleType = pdgDB->GetParticle(pdg);
 538     TString pname(particleType ? particleType->GetName() : "???");
 539     TString what;
 540     if (mc->IsTrackEntering())    what.Append("entering ");
 541     if (mc->IsTrackExiting())     what.Append("exiting ");
 542     if (mc->IsTrackInside())      what.Append("inside ");
 543     if (mc->IsTrackDisappeared()) what.Append("disappeared ");
 544     if (mc->IsTrackStop())        what.Append("stopped ");
 545     if (mc->IsNewTrack())         what.Append("new ");
 546     if (mc->IsTrackAlive())       what.Append("alive ");
 547     if (mc->IsTrackOut())         what.Append("out ");
 548
 549     Int_t mother = gAlice->GetMCApp()->GetPrimary(trackno);
 550     Warning("Step", "Track # %5d deposits a lot of energy\n"
 551             "  Volume:    %s\n"
 552             "  Momentum:  (%7.4f,%7.4f,%7.4f)\n"
 553             "  PDG:       %d (%s)\n"
 554             "  Edep:      %-14.7f keV (mother %d)\n"
 555             "  p/m:       %-7.4f/%-7.4f = %-14.7f\n"
 556             "  Processes: %s\n"
 557             "  What:      %s\n",
 558             trackno, mc->CurrentVolPath(), p.X(), p.Y(), p.Z(),
 559             pdg, pname.Data(), edep, mother, p.P(), mass,
 560             poverm, processes.Data(), what.Data());
 561   }
 562
 563   // Check that the track is actually within the active area
 564   Bool_t entering = mc->IsTrackEntering();
 565   Bool_t inside   = mc->IsTrackInside();
 566   Bool_t out      = (mc->IsTrackExiting()|| mc->IsTrackDisappeared()||
 567                      mc->IsTrackStop());
 568   // Reset the energy deposition for this track, and update some of
 569   // our parameters.
 570   if (entering) {
 571     AliDebug(15, Form("Track # %8d entering active FMD volume %s: "
 572                       "Edep=%f (%f,%f,%f)", trackno, mc->CurrentVolPath(),
 573                       edep, v.X(), v.Y(), v.Z()));
 574     fCurrentP      = p;
 575     fCurrentV      = v;
 576     fCurrentDeltaE = edep;
 577     fCurrentPdg    = pdg; // mc->IdFromPDG(pdg);
 578   }
 579   // If the track is inside, then update the energy deposition
 580   if (inside && fCurrentDeltaE >= 0) {
 581     fCurrentDeltaE += edep;
 582     AliDebug(15, Form("Track # %8d inside active FMD volume %s: Edep=%f, "
 583                       "Accumulated Edep=%f  (%f,%f,%f)", trackno,
 584                       mc->CurrentVolPath(), edep, fCurrentDeltaE,
 585                       v.X(), v.Y(), v.Z()));
 586   }
 587   // The track exits the volume, or it disappeared in the volume, or
 588   // the track is stopped because it no longer fulfills the cuts
 589   // defined, then we create a hit.
 590   if (out) {
 591     if (fCurrentDeltaE >= 0) {
 592       fCurrentDeltaE += edep;
 593       AliDebug(15, Form("Track # %8d exiting active FMD volume %s: Edep=%g, "
 594                         "Accumulated Edep=%g (%f,%f,%f)", trackno,
 595                         mc->CurrentVolPath(), edep, fCurrentDeltaE,
 596                         v.X(), v.Y(), v.Z()));
 597       AliFMDHit* h =
 598         fFMD->AddHitByFields(trackno, detector, ring, sector, strip,
 599                              fCurrentV.X(),  fCurrentV.Y(), fCurrentV.Z(),
 600                              fCurrentP.X(),  fCurrentP.Y(), fCurrentP.Z(),
 601                              fCurrentDeltaE, fCurrentPdg,   fCurrentV.T());
 602       // Add a copy
 603       if (isBad && fBad) {
 604         new ((*fBad)[fBad->GetEntries()]) AliFMDHit(*h);
 605       }
 606     }
 607     fCurrentDeltaE = -1;
 608   }
 609 }
 610
 611 //____________________________________________________________________
 612 void
 613 AliFMDSimulator::EndEvent()
 614 {
 615   if (fBad && fBad->GetEntries() > 0) {
 616     Warning("EndEvent", "got %d 'bad' hits", fBad->GetEntries());
 617     TIter next(fBad);
 618     AliFMDHit* hit;
 619     while ((hit = static_cast<AliFMDHit*>(next())))
 620       hit->Print("D");
 621     fBad->Clear();
 622   }
 623 }
 624
 625
 626 //____________________________________________________________________
 627 //
 628 // EOF
 629 //