Added the class AliFMDGeometryBuilder (and derived
[u/mrichter/AliRoot.git] / 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 //