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