d70ac4591c6a9c064ce99b9e7db664eb8e3ad23a
[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 // *  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, "FMD Si$", 
200                       a, z, density, radiationLength, absorbtionLength);
201   fFMD->AliMedium(kSiId, "FMD 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, "FMD Carbon$", 
217                       a, z, density, radiationLength, absorbtionLength);
218   fFMD->AliMedium(kCarbonId, "FMD 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, "FMD Aluminum$", 
229                       a, z, density, radiationLength, absorbtionLength);
230   fFMD->AliMedium(kAlId, "FMD 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, "FMD Si Chip$", as, zs, density, 6, ws);
250     fFMD->AliMedium(kSiChipId, "FMD 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, "FMD Kaption$", as, zs, density, 4, ws);
268     fFMD->AliMedium(kAlId, "FMD 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, "FMD Air$", as, zs, density, 4, ws);
286     fFMD->AliMedium(kAirId, "FMD 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, "FMD PCB$", as, zs, density, 14, ws);
312     fFMD->AliMedium(kPcbId, "FMD 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, "FMD Plastic$", as, zs, density, -2, ws);
329     fFMD->AliMedium(kPlasticId, "FMD 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   
371   if (!mc->IsTrackAlive()) return;
372   if (TMath::Abs(mc->TrackCharge()) <= 0) return;
373   
374   Int_t copy;
375   Int_t vol = mc->CurrentVolID(copy);
376   if (vol != fInnerId && vol != fOuterId) {
377     AliDebug(15, Form("Not an FMD volume %d '%s' (%d or %d)", 
378                       vol, mc->CurrentVolName(), fInnerId, fOuterId));
379     return;
380   }
381
382   // Check that the track is actually within the active area 
383   Bool_t entering = mc->IsTrackEntering();
384   Bool_t inside   = mc->IsTrackInside();
385   Bool_t out      = (mc->IsTrackExiting()|| mc->IsTrackDisappeared()||
386                      mc->IsTrackStop());
387
388   // Reset the energy deposition for this track, and update some of
389   // our parameters.
390   if (entering) {
391     AliDebug(15, "Entering active FMD volume");
392     fCurrentDeltaE = 0;
393
394     // Get production vertex and momentum of the track 
395     mc->TrackMomentum(fCurrentP);
396     mc->TrackPosition(fCurrentV);
397     fCurrentPdg = mc->IdFromPDG(mc->TrackPid());
398   }
399   
400   // If the track is inside, then update the energy deposition
401   if (inside && fCurrentDeltaE >= 0) 
402     AliDebug(15, "Inside active FMD volume");
403     fCurrentDeltaE += 1000 * mc->Edep();
404
405   // The track exits the volume, or it disappeared in the volume, or
406   // the track is stopped because it no longer fulfills the cuts
407   // defined, then we create a hit. 
408   if (out && fCurrentDeltaE >= 0) {
409     AliDebug(15, Form("Leaving active FMD volume %s", mc->CurrentVolPath()));
410
411     Int_t strip = copy - 1;
412     Int_t sectordiv;
413     mc->CurrentVolOffID(fSectorOff, sectordiv);
414     Int_t module;
415     mc->CurrentVolOffID(fModuleOff, module);
416     Int_t sector = 2 * module + sectordiv;
417     Int_t iring;
418     mc->CurrentVolOffID(fRingOff, iring);
419     Char_t ring = Char_t(iring);
420     Int_t detector;
421     mc->CurrentVolOffID(fDetectorOff, detector);
422
423
424     AliFMDGeometry* fmd = AliFMDGeometry::Instance();
425     Double_t rz = fmd->GetDetector(detector)->GetRingZ(ring);
426     Int_t    n  = fmd->GetDetector(detector)->GetRing(ring)->GetNSectors();
427     if (rz < 0) {
428       Int_t s = ((n - sector + n / 2) % n) + 1;
429       AliDebug(40, Form("Recalculating sector to %d (=%d-%d+%d/2%%%d+1 z=%f)", 
430                         s, n, sector, n, n, rz));
431       sector = s;
432     }
433     if (sector < 1 || sector > n) {
434       Warning("Step", "sector # %d out of range (0-%d)", sector-1, n-1);
435       return;
436     }
437     sector--;
438     fCurrentDeltaE += 1000 * mc->Edep();
439
440     AliDebug(20, Form("Processing hit in FMD%d%c[%2d,%3d]: %f", 
441                       detector, ring, sector, strip, fCurrentDeltaE));
442     
443     fFMD->AddHitByFields(gAlice->GetMCApp()->GetCurrentTrackNumber(),
444                          UShort_t(detector), ring, UShort_t(sector), 
445                          UShort_t(strip),
446                          fCurrentV.X(), fCurrentV.Y(), fCurrentV.Z(),
447                          fCurrentP.X(), fCurrentP.Y(), fCurrentP.Z(), 
448                          fCurrentDeltaE, fCurrentPdg, fCurrentV.T());
449     fCurrentDeltaE = -1;
450   }
451 }
452
453
454
455 //____________________________________________________________________
456 //
457 // EOF
458 //