Added a lot of Doxygen documentation
[u/mrichter/AliRoot.git] / FMD / AliFMD.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 // is the driver for especially simulation. 
22 //
23 // The Forward Multiplicity Detector consists of 3 sub-detectors FMD1,
24 // FMD2, and FMD3, each of which has 1 or 2 rings of silicon sensors. 
25 //                                                       
26 // This is the base class for all FMD manager classes. 
27 //                    
28 // The actual code is done by various separate classes.   Below is
29 // diagram showing the relationship between the various FMD classes
30 // that handles the simulation
31 //
32 //
33 //       +----------+   +----------+   
34 //       | AliFMDv1 |   | AliFMDv0 |   
35 //       +----------+   +----------+   
36 //            |              |                    +-----------------+
37 //       +----+--------------+                 +--| AliFMDDigitizer |
38 //       |                                     |  +-----------------+
39 //       |           +---------------------+   |
40 //       |        +--| AliFMDBaseDigitizer |<--+
41 //       V     1  |  +---------------------+   |
42 //  +--------+<>--+                            |  +------------------+
43 //  | AliFMD |                                 +--| AliFMDSDigitizer |    
44 //  +--------+<>--+                               +------------------+       
45 //             1  |  +---------------------+
46 //                +--| AliFMDReconstructor |
47 //                   +---------------------+
48 //
49 // *  AliFMD 
50 //    This defines the interface for the various parts of AliROOT that
51 //    uses the FMD, like AliFMDSimulator, AliFMDDigitizer, 
52 //    AliFMDReconstructor, and so on. 
53 //
54 // *  AliFMDv0
55 //    This is a concrete implementation of the AliFMD interface. 
56 //    It is the responsibility of this class to create the FMD
57 //    geometry.
58 //
59 // *  AliFMDv1 
60 //    This is a concrete implementation of the AliFMD interface. 
61 //    It is the responsibility of this class to create the FMD
62 //    geometry, process hits in the FMD, and serve hits and digits to
63 //    the various clients. 
64 //  
65 // *  AliFMDSimulator
66 //    This is the base class for the FMD simulation tasks.   The
67 //    simulator tasks are responsible to implment the geoemtry, and
68 //    process hits. 
69 //                                                                          
70 // *  AliFMDReconstructor
71 //    This is a concrete implementation of the AliReconstructor that
72 //    reconstructs pseudo-inclusive-multiplicities from digits (raw or
73 //    from simulation)
74 //
75 // Calibration and geometry parameters are managed by separate
76 // singleton managers.  These are AliFMDGeometry and
77 // AliFMDParameters.  Please refer to these classes for more
78 // information on these.
79 //
80
81 // These files are not in the same directory, so there's no reason to
82 // ask the preprocessor to search in the current directory for these
83 // files by including them with `#include "..."' 
84 #include <math.h>               // __CMATH__
85 #include <TClonesArray.h>       // ROOT_TClonesArray
86 #include <TGeometry.h>          // ROOT_TGeomtry
87 #include <TNode.h>              // ROOT_TNode
88 #include <TXTRU.h>              // ROOT_TXTRU
89 #include <TRotMatrix.h>         // ROOT_TRotMatrix
90 #include <TTUBE.h>              // ROOT_TTUBE
91 #include <TTree.h>              // ROOT_TTree
92 #include <TBrowser.h>           // ROOT_TBrowser
93 #include <TMath.h>              // ROOT_TMath
94 #include <TVirtualMC.h>         // ROOT_TVirtualMC
95 #include <TVector2.h>
96 #include <TVector3.h>
97 #include <TMarker3DBox.h>
98
99 #include <AliRunDigitizer.h>    // ALIRUNDIGITIZER_H
100 #include <AliLoader.h>          // ALILOADER_H
101 #include <AliRun.h>             // ALIRUN_H
102 #include <AliMC.h>              // ALIMC_H
103 #include "AliMagF.h"            // ALIMAGF_H
104 #include <AliLog.h>             // ALILOG_H
105 #include "AliFMD.h"             // ALIFMD_H
106 #include "AliFMDDigit.h"        // ALIFMDDIGIG_H
107 #include "AliFMDHit.h"          // ALIFMDHIT_H
108 #include "AliFMDGeometry.h"     // ALIFMDGEOMETRY_H
109 #include "AliFMDDetector.h"     // ALIFMDDETECTOR_H
110 #include "AliFMDRing.h"         // ALIFMDRING_H
111 #include "AliFMDDigitizer.h"    // ALIFMDDIGITIZER_H
112 #include "AliPoints.h"          // ALIPOINTS_H
113 #include "AliFMDGeometryBuilder.h"
114 #include "AliFMDRawWriter.h"    // ALIFMDRAWWRITER_H
115
116
117 class AliFMDPoints : public AliPoints 
118 {
119 public:
120   /** Constructor 
121       @param hit   Hit to draw
122       @param color Color of hit */
123   AliFMDPoints(AliFMDHit* hit, UInt_t color) 
124     : AliPoints(1), fMarker(0)
125   {
126     if (!hit) return;
127     Float_t  size  = TMath::Min(TMath::Max(hit->Edep() * .1, .1), 1.);
128     TVector3 p(hit->Px(), hit->Py(), hit->Pz());
129     fMarker = new TMarker3DBox(hit->X(), hit->Y(), hit->Z(), size, size, size,
130                                p.Theta(), p.Phi());
131     fMarker->SetLineColor(color);
132     fMarker->SetRefObject(this);
133     fP[0] = hit->X();
134     fP[1] = hit->Y();
135     fP[2] = hit->Z();
136   }
137   /** Destructor  */
138   ~AliFMDPoints() 
139   {
140     // if (fMarker) delete  fMarker;
141   }
142   void SetXYZ(Double_t x, Double_t y, Double_t z)
143   {
144     if (fMarker) fMarker->SetPosition(x, y, z);
145   }
146   Int_t DistancetoPrimitive(Int_t px, Int_t py) 
147   {
148     return fMarker->DistancetoPrimitive(px, py);
149   }
150   void Draw(Option_t* option) 
151   {
152     if (fMarker) fMarker->Draw(option);
153   }
154   void Paint(Option_t* option)
155   {
156     if (fMarker) fMarker->Paint(option);
157   }
158   void SetMarkerColor(Color_t colour) 
159   {
160     if (fMarker) fMarker->SetLineColor(colour);
161   }
162 private:
163   TMarker3DBox* fMarker;
164 };
165
166
167 //____________________________________________________________________
168 ClassImp(AliFMD)
169 #if 0
170   ; // This is to keep Emacs from indenting the next line 
171 #endif 
172
173 //____________________________________________________________________
174 AliFMD::AliFMD()
175   : AliDetector(),
176     fSDigits(0), 
177     fNsdigits(0),
178     fDetailed(kTRUE),
179     fBad(0)
180 {
181   //
182   // Default constructor for class AliFMD
183   //
184   AliDebug(10, "\tDefault CTOR");
185   fHits        = 0;
186   fDigits      = 0;
187   fIshunt      = 0;
188   fUseOld      = kFALSE;
189   fUseAssembly = kTRUE;
190   fBad         = new TClonesArray("AliFMDHit");
191 }
192
193 //____________________________________________________________________
194 AliFMD::AliFMD(const AliFMD& other)
195   : AliDetector(other),
196     fSDigits(other.fSDigits), 
197     fNsdigits(other.fNsdigits),
198     fDetailed(other.fDetailed),
199     fBad(other.fBad)
200 {
201   // Copy constructor 
202   fUseOld      = other.fUseOld;
203   fUseAssembly = other.fUseAssembly;
204 }
205
206 //____________________________________________________________________
207 AliFMD::AliFMD(const char *name, const char *title)
208   : AliDetector (name, title),
209     fSDigits(0),
210     fNsdigits(0),
211     fDetailed(kTRUE),
212     fBad(0)
213 {
214   //
215   // Standard constructor for Forward Multiplicity Detector
216   //
217   AliDebug(10, "\tStandard CTOR");
218   fUseOld      = kFALSE;
219   fUseAssembly = kFALSE;
220   fBad         = new TClonesArray("AliFMDHit");
221   
222   // Initialise Hit array
223   HitsArray();
224   gAlice->GetMCApp()->AddHitList(fHits);
225
226   // (S)Digits for the detectors disk
227   DigitsArray();
228   SDigitsArray();
229   
230   // CHC: What is this?
231   fIshunt = 0;
232   SetMarkerColor(kRed);
233   SetLineColor(kYellow);
234 }
235
236 //____________________________________________________________________
237 AliFMD::~AliFMD ()
238 {
239   // Destructor for base class AliFMD
240   if (fHits) {
241     fHits->Delete();
242     delete fHits;
243     fHits = 0;
244   }
245   if (fDigits) {
246     fDigits->Delete();
247     delete fDigits;
248     fDigits = 0;
249   }
250   if (fSDigits) {
251     fSDigits->Delete();
252     delete fSDigits;
253     fSDigits = 0;
254   }
255   if (fBad) {
256     fBad->Delete();
257     delete fBad;
258     fBad = 0;
259   }
260 }
261
262 //____________________________________________________________________
263 AliFMD&
264 AliFMD::operator=(const AliFMD& other)
265 {
266   // Assignment operator
267   AliDetector::operator=(other);
268   fSDigits              = other.fSDigits; 
269   fNsdigits             = other.fNsdigits;
270   fDetailed             = other.fDetailed;
271   fBad                  = other.fBad;
272   return *this;
273 }
274
275 //====================================================================
276 //
277 // GEometry ANd Traking
278 //
279 //____________________________________________________________________
280 void 
281 AliFMD::CreateGeometry()
282 {
283   //
284   // Create the geometry of Forward Multiplicity Detector.  The actual
285   // construction of the geometry is delegated to the class
286   // AliFMDGeometryBuilder, invoked by the singleton manager
287   // AliFMDGeometry. 
288   //
289   AliFMDGeometry*  fmd = AliFMDGeometry::Instance();
290   fmd->SetDetailed(fDetailed);
291   fmd->UseAssembly(fUseAssembly);
292   fmd->Build();
293 }    
294
295 //____________________________________________________________________
296 void AliFMD::CreateMaterials() 
297 {
298   // Define the materials and tracking mediums needed by the FMD
299   // simulation.   These mediums are made by sending the messages
300   // AliMaterial, AliMixture, and AliMedium to the passed AliModule
301   // object module.   The defined mediums are 
302   // 
303   //    FMD Si$         Silicon (active medium in sensors)
304   //    FMD C$          Carbon fibre (support cone for FMD3 and vacuum pipe)
305   //    FMD Al$         Aluminium (honeycomb support plates)
306   //    FMD PCB$        Printed Circuit Board (FEE board with VA1_3)
307   //    FMD Chip$       Electronics chips (currently not used)
308   //    FMD Air$        Air (Air in the FMD)
309   //    FMD Plastic$    Plastic (Support legs for the hybrid cards)
310   //
311   // The geometry builder should really be the one that creates the
312   // materials, but the architecture of AliROOT makes that design
313   // akward.  What should happen, was that the AliFMDGeometryBuilder
314   // made the mediums, and that this class retrives pointers from the
315   // TGeoManager, and registers the mediums here.  Alas, it's not
316   // really that easy. 
317   //
318   AliDebug(10, "\tCreating materials");
319   // Get pointer to geometry singleton object. 
320   AliFMDGeometry* geometry = AliFMDGeometry::Instance();
321   geometry->Init();
322 #if 0
323   if (gGeoManager && gGeoManager->GetMedium("FMD Si$")) {
324     // We need to figure out the some stuff about the geometry
325     fmd->ExtractGeomInfo();
326     return;
327   }
328 #endif  
329   Int_t    id;
330   Double_t a                = 0;
331   Double_t z                = 0;
332   Double_t density          = 0;
333   Double_t radiationLength  = 0;
334   Double_t absorbtionLength = 999;
335   Int_t    fieldType        = gAlice->Field()->Integ();     // Field type 
336   Double_t maxField         = gAlice->Field()->Max();     // Field max.
337   Double_t maxBending       = 0;     // Max Angle
338   Double_t maxStepSize      = 0.001; // Max step size 
339   Double_t maxEnergyLoss    = 1;     // Max Delta E
340   Double_t precision        = 0.001; // Precision
341   Double_t minStepSize      = 0.001; // Minimum step size 
342  
343   // Silicon 
344   a                = 28.0855;
345   z                = 14.;
346   density          = geometry->GetSiDensity();
347   radiationLength  = 9.36;
348   maxBending       = 1;
349   maxStepSize      = .001;
350   precision        = .001;
351   minStepSize      = .001;
352   id               = kSiId;
353   AliMaterial(id, "Si$", a, z, density, radiationLength, absorbtionLength);
354   AliMedium(kSiId, "Si$", id,1,fieldType,maxField,maxBending,
355             maxStepSize,maxEnergyLoss,precision,minStepSize);
356   
357
358   // Carbon 
359   a                = 12.011;
360   z                = 6.;
361   density          = 2.265;
362   radiationLength  = 18.8;
363   maxBending       = 10;
364   maxStepSize      = .01;
365   precision        = .003;
366   minStepSize      = .003;
367   id               = kCarbonId;
368   AliMaterial(id, "Carbon$", a, z, density, radiationLength, absorbtionLength);
369   AliMedium(kCarbonId, "Carbon$", id,0,fieldType,maxField,maxBending,
370                     maxStepSize,maxEnergyLoss,precision,minStepSize);
371
372   // Aluminum
373   a                = 26.981539;
374   z                = 13.;
375   density          = 2.7;
376   radiationLength  = 8.9;
377   id               = kAlId;
378   AliMaterial(id, "Aluminum$",a,z, density, radiationLength, absorbtionLength);
379   AliMedium(kAlId, "Aluminum$", id, 0, fieldType, maxField, maxBending,
380             maxStepSize, maxEnergyLoss, precision, minStepSize);
381   
382   
383   // Copper 
384   a                = 63.546;
385   z                = 29;
386   density          =  8.96;
387   radiationLength  =  1.43;
388   id               = kCopperId;
389   AliMaterial(id, "Copper$", 
390                       a, z, density, radiationLength, absorbtionLength);
391   AliMedium(kCopperId, "Copper$", id, 0, fieldType, maxField, maxBending,
392             maxStepSize, maxEnergyLoss, precision, minStepSize);
393   
394
395   // Silicon chip 
396   {
397     Float_t as[] = { 12.0107,      14.0067,      15.9994,
398                       1.00794,     28.0855,     107.8682 };
399     Float_t zs[] = {  6.,           7.,           8.,
400                       1.,          14.,          47. };
401     Float_t ws[] = {  0.039730642,  0.001396798,  0.01169634,
402                       0.004367771,  0.844665,     0.09814344903 };
403     density          = 2.36436;
404     maxBending       = 10;
405     maxStepSize      = .01;
406     precision        = .003;
407     minStepSize      = .003;
408     id               = kSiChipId;
409     AliMixture(id, "Si Chip$", as, zs, density, 6, ws);
410     AliMedium(kSiChipId, "Si Chip$",  id, 0, fieldType, maxField, maxBending, 
411               maxStepSize, maxEnergyLoss, precision, minStepSize);
412   }
413   
414   // Kaption
415   {
416     Float_t as[] = { 1.00794,  12.0107,  14.010,   15.9994};
417     Float_t zs[] = { 1.,        6.,       7.,       8.};
418     Float_t ws[] = { 0.026362,  0.69113,  0.07327,  0.209235};
419     density          = 1.42;
420     maxBending       = 1;
421     maxStepSize      = .001;
422     precision        = .001;
423     minStepSize      = .001;
424     id               = kKaptonId;
425     AliMixture(id, "Kaption$", as, zs, density, 4, ws);
426     AliMedium(kKaptonId, "Kaption$", id,0,fieldType,maxField,maxBending,
427               maxStepSize,maxEnergyLoss,precision,minStepSize);
428   }
429
430   // Air
431   {
432     Float_t as[] = { 12.0107, 14.0067,   15.9994,  39.948 };
433     Float_t zs[] = {  6.,      7.,       8.,       18. };
434     Float_t ws[] = { 0.000124, 0.755267, 0.231781, 0.012827 }; 
435     density      = .00120479;
436     maxBending   = 1;
437     maxStepSize  = .001;
438     precision    = .001;
439     minStepSize  = .001;
440     id           = kAirId;
441     AliMixture(id, "Air$", as, zs, density, 4, ws);
442     AliMedium(kAirId, "Air$", id,0,fieldType,maxField,maxBending,
443               maxStepSize,maxEnergyLoss,precision,minStepSize);
444   }
445   
446   // PCB
447   {
448     Float_t zs[] = { 14.,         20.,         13.,         12.,
449                       5.,         22.,         11.,         19.,
450                      26.,          9.,          8.,          6.,
451                       7.,          1.};
452     Float_t as[] = { 28.0855,     40.078,      26.981538,   24.305, 
453                      10.811,      47.867,      22.98977,    39.0983,
454                      55.845,      18.9984,     15.9994,     12.0107,
455                      14.0067,      1.00794};
456     Float_t ws[] = {  0.15144894,  0.08147477,  0.04128158,  0.00904554, 
457                       0.01397570,  0.00287685,  0.00445114,  0.00498089,
458                       0.00209828,  0.00420000,  0.36043788,  0.27529426,
459                       0.01415852,  0.03427566};
460     density      = 1.8;
461     maxBending   = 1;
462     maxStepSize  = .001;
463     precision    = .001;
464     minStepSize  = .001;
465     id           = kPcbId;
466     AliMixture(id, "PCB$", as, zs, density, 14, ws);
467     AliMedium(kPcbId, "PCB$", id,0,fieldType,maxField,maxBending,
468               maxStepSize,maxEnergyLoss,precision,minStepSize);
469   }
470   
471   // Plastic 
472   {
473     Float_t as[] = { 1.01, 12.01 };
474     Float_t zs[] = { 1.,   6.    };
475     Float_t ws[] = { 1.,   1.    };
476     density      = 1.03;
477     maxBending   = 10;
478     maxStepSize  = .01;
479     precision    = .003;
480     minStepSize  = .003;
481     id           = kPlasticId;
482     AliMixture(id, "Plastic$", as, zs, density, -2, ws);
483     AliMedium(kPlasticId, "Plastic$", id,0,fieldType,maxField,maxBending,
484               maxStepSize,maxEnergyLoss,precision,minStepSize);
485   }
486 }
487
488 //____________________________________________________________________
489 void  
490 AliFMD::Init()
491 {
492   // Initialize the detector 
493   // 
494   AliDebug(1, "Initialising FMD detector object");
495   // AliFMDGeometry*  fmd = AliFMDGeometry::Instance();
496   // fmd->InitTransformations();
497 }
498
499 //____________________________________________________________________
500 void
501 AliFMD::FinishEvent()
502 {
503   // Called at the end of the an event in simulations.  If the debug
504   // level is high enough, then the `bad' hits are printed.
505   // 
506   if (AliLog::GetDebugLevel("FMD", "AliFMD") < 10) return;
507   if (fBad && fBad->GetEntries() > 0) {
508     AliWarning((Form("EndEvent", "got %d 'bad' hits", fBad->GetEntries())));
509     TIter next(fBad);
510     AliFMDHit* hit;
511     while ((hit = static_cast<AliFMDHit*>(next()))) hit->Print("D");
512     fBad->Clear();
513   }
514 }
515
516
517 //====================================================================
518 //
519 // Graphics and event display
520 //
521 //____________________________________________________________________
522 void 
523 AliFMD::BuildGeometry()
524 {
525   //
526   // Build simple ROOT TNode geometry for event display. With the new
527   // geometry modeller, TGeoManager, this seems rather redundant. 
528   AliDebug(10, "\tCreating a simplified geometry");
529
530   AliFMDGeometry* fmd = AliFMDGeometry::Instance();
531   
532   static TXTRU*     innerShape = 0;
533   static TXTRU*     outerShape = 0;
534   static TObjArray* innerRot   = 0;
535   static TObjArray* outerRot   = 0;
536
537   if (!innerShape || !outerShape) {
538     // Make the shapes for the modules 
539     for (Int_t i = 0; i < 2; i++) {
540       AliFMDRing* r = 0;
541       switch (i) {
542       case 0: r = fmd->GetRing('I'); break;
543       case 1: r = fmd->GetRing('O'); break;
544       }
545       if (!r) {
546         AliError(Form("no ring found for i=%d", i));
547         return;
548       }
549       Double_t    siThick  = r->GetSiThickness();
550       const Int_t nv       = r->GetNVerticies();
551       Double_t    theta    = r->GetTheta();
552       Int_t       nmod     = r->GetNModules();
553       
554       TXTRU* shape = new TXTRU(r->GetName(), r->GetTitle(), "void", nv, 2);
555       for (Int_t j = 0; j < nv; j++) {
556         TVector2* vv = r->GetVertex(nv - 1 - j);
557         shape->DefineVertex(j, vv->X(), vv->Y());
558       }
559       shape->DefineSection(0, -siThick / 2, 1, 0, 0);
560       shape->DefineSection(1, +siThick / 2, 1, 0, 0);
561       shape->SetLineColor(GetLineColor());
562       
563       TObjArray* rots = new TObjArray(nmod);
564       for (Int_t j = 0; j < nmod; j++) {
565         Double_t th = (j + .5) * theta * 2;
566         TString name(Form("FMD_ring_%c_rot_%02d", r->GetId(), j));
567         TString title(Form("FMD Ring %c Rotation # %d", r->GetId(), j));
568         TRotMatrix* rot = new TRotMatrix(name.Data(), title.Data(),
569                                          90, th, 90, fmod(90+th,360), 0, 0);
570         rots->AddAt(rot, j);
571       }
572       
573       switch (r->GetId()) {
574       case 'i':
575       case 'I': innerShape = shape; innerRot = rots; break;
576       case 'o':
577       case 'O': outerShape = shape; outerRot = rots; break;
578       }
579     }
580   }
581   
582   TNode* top = gAlice->GetGeometry()->GetNode("alice");
583   
584   for (Int_t i = 1; i <= 3; i++) {
585     AliFMDDetector* det = fmd->GetDetector(i);
586     if (!det) {
587       Warning("BuildGeometry", "FMD%d seems to be disabled", i);
588       continue;
589     }
590     Double_t w  = 0;
591     Double_t rh = det->GetRing('I')->GetHighR();
592     Char_t   id = 'I';
593     if (det->GetRing('O')) {
594       w  = TMath::Abs(det->GetRingZ('O') - det->GetRingZ('I'));
595       id = (TMath::Abs(det->GetRingZ('O')) 
596             > TMath::Abs(det->GetRingZ('I')) ? 'O' : 'I');
597       rh = det->GetRing('O')->GetHighR();
598     }
599     w += (det->GetRing(id)->GetModuleSpacing() +
600           det->GetRing(id)->GetSiThickness());
601     TShape* shape = new TTUBE(det->GetName(), det->GetTitle(), "void",
602                               det->GetRing('I')->GetLowR(), rh, w / 2);
603     Double_t z = (det->GetRingZ('I') - w / 2);
604     if (z > 0) z += det->GetRing(id)->GetModuleSpacing();
605     top->cd();
606     TNode* node = new TNode(det->GetName(), det->GetTitle(), shape, 
607                             0, 0, z, 0);
608     fNodes->Add(node);
609     
610     for (Int_t j = 0; j < 2; j++) {
611       AliFMDRing* r      = 0;
612       TShape*     rshape = 0;
613       TObjArray*  rots   = 0;
614       switch (j) {
615       case 0: 
616         r = det->GetRing('I'); rshape = innerShape; rots = innerRot; break;
617       case 1: 
618         r = det->GetRing('O'); rshape = outerShape; rots = outerRot; break;
619       }
620       if (!r) continue;
621       
622       Double_t    siThick  = r->GetSiThickness();
623       Int_t       nmod     = r->GetNModules();
624       Double_t    modspace = r->GetModuleSpacing();
625       Double_t    rz       = - (z - det->GetRingZ(r->GetId()));
626       
627       for (Int_t k = 0; k < nmod; k++) {
628         node->cd();
629         Double_t    offz    = (k % 2 == 1 ? modspace : 0);
630         TRotMatrix* rot     = static_cast<TRotMatrix*>(rots->At(k));
631         TString name(Form("%s%c_module_%02d", det->GetName(), r->GetId(),k));
632         TString title(Form("%s%c Module %d", det->GetName(), r->GetId(),k));
633         TNode* mnod = new TNode(name.Data(), title.Data(), rshape, 
634                                 0, 0, rz - siThick / 2 
635                                 + TMath::Sign(offz,z), rot);
636         mnod->SetLineColor(GetLineColor());
637         fNodes->Add(mnod);
638       } // for (Int_t k = 0 ; ...)
639     } // for (Int_t j = 0 ; ...)
640   } // for (Int_t i = 1 ; ...)
641 }
642
643 //____________________________________________________________________
644 void 
645 AliFMD::LoadPoints(Int_t /* track */) 
646 {
647   // Store x, y, z of all hits in memory for display. 
648   // 
649   // Normally, the hits are drawn using TPolyMarker3D - however, that
650   // is not very useful for the FMD.  Therefor, this member function
651   // is overloaded to make TMarker3D, via the class AliFMDPoints.
652   // AliFMDPoints is a local class. 
653   //
654   if (!fHits) {
655     AliError(Form("fHits == 0. Name is %s",GetName()));
656     return;
657   }
658   Int_t nHits = fHits->GetEntriesFast();
659   if (nHits == 0) {
660     return;
661   }
662   Int_t tracks = gAlice->GetMCApp()->GetNtrack();
663   if (fPoints == 0) fPoints = new TObjArray(2 * tracks);
664
665   // Get geometry 
666   AliFMDGeometry* geom = AliFMDGeometry::Instance();
667   geom->Init();
668   geom->InitTransformations();
669
670   // Now make markers for each hit  
671   // AliInfo(Form("Drawing %d hits (have %d points) for track %d", 
672   //              nHits, fPoints->GetEntriesFast(), track));
673   for (Int_t ihit = 0; ihit < nHits; ihit++) {
674     AliFMDHit* hit = static_cast<AliFMDHit*>(fHits->At(ihit));
675     if (!hit) continue;
676     Double_t edep    = hit->Edep();
677     Double_t m       = hit->M();
678     Double_t poverm  = (m == 0 ? 0 : hit->P());
679     Double_t absQ    = TMath::Abs(hit->Q());
680     Bool_t   bad     = kFALSE;
681     // This `if' is to debug abnormal energy depositions.  We trigger on
682     // p/m approx larger than or equal to a MIP, and a large edep - more 
683     // than 1 keV - a MIP is 100 eV. 
684     if (edep > absQ * absQ && poverm > 1) bad = kTRUE;
685
686     AliFMDPoints* p1 = new AliFMDPoints(hit, GetMarkerColor());
687     // AliPoints* p1 = new AliPoints();
688     // p1->SetMarkerColor(GetMarkerColor());
689     // p1->SetMarkerSize(GetMarkerSize());
690     // p1->SetPoint(0, hit->X(), hit->Y(), hit->Z());
691     p1->SetDetector(this);
692     p1->SetParticle(hit->GetTrack());
693     fPoints->AddAt(p1, hit->GetTrack());
694     if (bad) {
695       p1->SetMarkerColor(4);
696       // p1->SetMarkerSize(2 * GetMarkerSize());
697     }
698     
699     Double_t x, y, z;
700     geom->Detector2XYZ(hit->Detector(), hit->Ring(), hit->Sector(), 
701                        hit->Strip(), x, y, z);
702     AliFMDPoints* p = new AliFMDPoints(hit, 3);
703     // AliPoints* p = new AliPoints();
704     // p->SetMarkerColor(3);
705     // p->SetMarkerSize(GetMarkerSize());
706     // p->SetPoint(0, x, y, z);
707     p->SetDetector(this);
708     p->SetParticle(hit->GetTrack());
709     p->SetXYZ(x, y, z);
710     p->SetMarkerColor(3);
711     fPoints->AddAt(p, tracks+hit->GetTrack());
712     if (bad) {
713       p->SetMarkerColor(5);
714       // p->SetMarkerSize(2 * GetMarkerSize());
715     }
716     // AliInfo(Form("Adding point at %d", tracks+hit->GetTrack()));
717   }
718 }
719
720 //____________________________________________________________________
721 void 
722 AliFMD::DrawDetector()
723 {
724   // Draw a shaded view of the Forward multiplicity detector.  This
725   // isn't really useful anymore. 
726   AliDebug(10, "\tDraw detector");
727 }
728
729 //____________________________________________________________________
730 Int_t 
731 AliFMD::DistanceToPrimitive(Int_t, Int_t)
732 {
733   // Calculate the distance from the mouse to the FMD on the screen
734   // Dummy routine.
735   //
736   return 9999;
737 }
738
739 //====================================================================
740 //
741 // Hit and Digit managment 
742 //
743 //____________________________________________________________________
744 void 
745 AliFMD::MakeBranch(Option_t * option)
746 {
747   // Create Tree branches for the FMD.
748   //
749   // Options:
750   //
751   //    H          Make a branch of TClonesArray of AliFMDHit's
752   //    D          Make a branch of TClonesArray of AliFMDDigit's
753   //    S          Make a branch of TClonesArray of AliFMDSDigit's
754   // 
755   const Int_t kBufferSize = 16000;
756   TString branchname(GetName());
757   TString opt(option);
758   
759   if (opt.Contains("H", TString::kIgnoreCase)) {
760     HitsArray();
761     AliDetector::MakeBranch(option); 
762   }
763   if (opt.Contains("D", TString::kIgnoreCase)) { 
764     DigitsArray();
765     MakeBranchInTree(fLoader->TreeD(), branchname.Data(),
766                      &fDigits, kBufferSize, 0);
767   }
768   if (opt.Contains("S", TString::kIgnoreCase)) { 
769     SDigitsArray();
770     MakeBranchInTree(fLoader->TreeS(), branchname.Data(),
771                      &fSDigits, kBufferSize, 0);
772   }
773 }
774
775 //____________________________________________________________________
776 void 
777 AliFMD::SetTreeAddress()
778 {
779   // Set branch address for the Hits, Digits, and SDigits Tree.
780   if (fLoader->TreeH()) HitsArray();
781   AliDetector::SetTreeAddress();
782
783   TTree *treeD = fLoader->TreeD();
784   if (treeD) {
785     DigitsArray();
786     TBranch* branch = treeD->GetBranch ("FMD");
787     if (branch) branch->SetAddress(&fDigits);
788   }
789
790   TTree *treeS = fLoader->TreeS();
791   if (treeS) {
792     SDigitsArray();
793     TBranch* branch = treeS->GetBranch ("FMD");
794     if (branch) branch->SetAddress(&fSDigits);
795   }
796 }
797
798 //____________________________________________________________________
799 void 
800 AliFMD::SetHitsAddressBranch(TBranch *b)
801 {
802   // Set the TClonesArray to read hits into. 
803   b->SetAddress(&fHits);
804 }
805
806 //____________________________________________________________________
807 void 
808 AliFMD::AddHit(Int_t track, Int_t *vol, Float_t *hits) 
809 {
810   // Add a hit to the hits tree 
811   // 
812   // The information of the two arrays are decoded as 
813   // 
814   // Parameters
815   //    track                Track #
816   //    ivol[0]  [UShort_t ] Detector # 
817   //    ivol[1]  [Char_t   ] Ring ID 
818   //    ivol[2]  [UShort_t ] Sector #
819   //    ivol[3]  [UShort_t ] Strip # 
820   //    hits[0]  [Float_t  ] Track's X-coordinate at hit 
821   //    hits[1]  [Float_t  ] Track's Y-coordinate at hit
822   //    hits[3]  [Float_t  ] Track's Z-coordinate at hit
823   //    hits[4]  [Float_t  ] X-component of track's momentum             
824   //    hits[5]  [Float_t  ] Y-component of track's momentum             
825   //    hits[6]  [Float_t  ] Z-component of track's momentum            
826   //    hits[7]  [Float_t  ] Energy deposited by track                  
827   //    hits[8]  [Int_t    ] Track's particle Id # 
828   //    hits[9]  [Float_t  ] Time when the track hit
829   // 
830   // 
831   AddHitByFields(track, 
832                  UShort_t(vol[0]),  // Detector # 
833                  Char_t(vol[1]),    // Ring ID
834                  UShort_t(vol[2]),  // Sector # 
835                  UShort_t(vol[3]),  // Strip # 
836                  hits[0],           // X
837                  hits[1],           // Y
838                  hits[2],           // Z
839                  hits[3],           // Px
840                  hits[4],           // Py
841                  hits[5],           // Pz
842                  hits[6],           // Energy loss 
843                  Int_t(hits[7]),    // PDG 
844                  hits[8]);          // Time
845 }
846
847 //____________________________________________________________________
848 AliFMDHit*
849 AliFMD::AddHitByFields(Int_t    track, 
850                        UShort_t detector, 
851                        Char_t   ring, 
852                        UShort_t sector, 
853                        UShort_t strip, 
854                        Float_t  x, 
855                        Float_t  y, 
856                        Float_t  z,
857                        Float_t  px, 
858                        Float_t  py, 
859                        Float_t  pz,
860                        Float_t  edep,
861                        Int_t    pdg,
862                        Float_t  t, 
863                        Float_t  l, 
864                        Bool_t   stop)
865 {
866   // Add a hit to the list
867   //
868   // Parameters:
869   // 
870   //    track     Track #
871   //    detector  Detector # (1, 2, or 3)                      
872   //    ring      Ring ID ('I' or 'O')
873   //    sector    Sector # (For inner/outer rings: 0-19/0-39)
874   //    strip     Strip # (For inner/outer rings: 0-511/0-255)
875   //    x         Track's X-coordinate at hit
876   //    y         Track's Y-coordinate at hit
877   //    z         Track's Z-coordinate at hit
878   //    px        X-component of track's momentum 
879   //    py        Y-component of track's momentum
880   //    pz        Z-component of track's momentum
881   //    edep      Energy deposited by track
882   //    pdg       Track's particle Id #
883   //    t         Time when the track hit 
884   //    l         Track length through the material. 
885   //    stop      Whether track was stopped or disappeared
886   // 
887   TClonesArray& a = *(HitsArray());
888   // Search through the list of already registered hits, and see if we
889   // find a hit with the same parameters.  If we do, then don't create
890   // a new hit, but rather update the energy deposited in the hit.
891   // This is done, so that a FLUKA based simulation will get the
892   // number of hits right, not just the enerrgy deposition. 
893   AliFMDHit* hit = 0;
894   for (Int_t i = 0; i < fNhits; i++) {
895     if (!a.At(i)) continue;
896     hit = static_cast<AliFMDHit*>(a.At(i));
897     if (hit->Detector() == detector 
898         && hit->Ring() == ring
899         && hit->Sector() == sector 
900         && hit->Strip() == strip
901         && hit->Track() == track) {
902       AliDebug(1, Form("already had a hit in FMD%d%c[%2d,%3d] for track # %d,"
903                        " adding energy (%f) to that hit (%f) -> %f", 
904                        detector, ring, sector, strip, track, edep, hit->Edep(),
905                        hit->Edep() + edep));
906       hit->SetEdep(hit->Edep() + edep);
907       return hit;
908     }
909   }
910   // If hit wasn't already registered, do so know. 
911   hit = new (a[fNhits]) AliFMDHit(fIshunt, track, detector, ring, sector, 
912                                   strip, x, y, z, px, py, pz, edep, pdg, t, 
913                                   l, stop);
914   fNhits++;
915   return hit;
916 }
917
918 //____________________________________________________________________
919 void 
920 AliFMD::AddDigit(Int_t* digits, Int_t*)
921 {
922   // Add a digit to the Digit tree 
923   // 
924   // Paramters 
925   //
926   //    digits[0]  [UShort_t] Detector #
927   //    digits[1]  [Char_t]   Ring ID
928   //    digits[2]  [UShort_t] Sector #
929   //    digits[3]  [UShort_t] Strip #
930   //    digits[4]  [UShort_t] ADC Count 
931   //    digits[5]  [Short_t]  ADC Count, -1 if not used
932   //    digits[6]  [Short_t]  ADC Count, -1 if not used 
933   // 
934   AddDigitByFields(UShort_t(digits[0]),  // Detector #
935                    Char_t(digits[1]),    // Ring ID
936                    UShort_t(digits[2]),  // Sector #
937                    UShort_t(digits[3]),  // Strip #
938                    UShort_t(digits[4]),  // ADC Count1 
939                    Short_t(digits[5]),   // ADC Count2 
940                    Short_t(digits[6]));  // ADC Count3 
941 }
942
943 //____________________________________________________________________
944 void 
945 AliFMD::AddDigitByFields(UShort_t detector, 
946                          Char_t   ring, 
947                          UShort_t sector, 
948                          UShort_t strip, 
949                          UShort_t count1, 
950                          Short_t  count2,
951                          Short_t  count3)
952 {
953   // add a real digit - as coming from data
954   // 
955   // Parameters 
956   //
957   //    detector  Detector # (1, 2, or 3)                      
958   //    ring      Ring ID ('I' or 'O')
959   //    sector    Sector # (For inner/outer rings: 0-19/0-39)
960   //    strip     Strip # (For inner/outer rings: 0-511/0-255)
961   //    count1    ADC count (a 10-bit word)
962   //    count2    ADC count (a 10-bit word), or -1 if not used
963   //    count3    ADC count (a 10-bit word), or -1 if not used
964   TClonesArray& a = *(DigitsArray());
965   
966   new (a[fNdigits++]) 
967     AliFMDDigit(detector, ring, sector, strip, count1, count2, count3);
968 }
969
970 //____________________________________________________________________
971 void 
972 AliFMD::AddSDigit(Int_t* digits)
973 {
974   // Add a digit to the SDigit tree 
975   // 
976   // Paramters 
977   //
978   //    digits[0]  [UShort_t] Detector #
979   //    digits[1]  [Char_t]   Ring ID
980   //    digits[2]  [UShort_t] Sector #
981   //    digits[3]  [UShort_t] Strip #
982   //    digits[4]  [Float_t]  Total energy deposited 
983   //    digits[5]  [UShort_t] ADC Count 
984   //    digits[6]  [Short_t]  ADC Count, -1 if not used
985   //    digits[7]  [Short_t]  ADC Count, -1 if not used 
986   // 
987   AddSDigitByFields(UShort_t(digits[0]),  // Detector #
988                     Char_t(digits[1]),    // Ring ID
989                     UShort_t(digits[2]),  // Sector #
990                     UShort_t(digits[3]),  // Strip #
991                     Float_t(digits[4]),   // Edep
992                     UShort_t(digits[5]),  // ADC Count1 
993                     Short_t(digits[6]),   // ADC Count2 
994                     Short_t(digits[7]));  // ADC Count3 
995 }
996
997 //____________________________________________________________________
998 void 
999 AliFMD::AddSDigitByFields(UShort_t detector, 
1000                           Char_t   ring, 
1001                           UShort_t sector, 
1002                           UShort_t strip, 
1003                           Float_t  edep,
1004                           UShort_t count1, 
1005                           Short_t  count2,
1006                           Short_t  count3)
1007 {
1008   // add a summable digit
1009   // 
1010   // Parameters 
1011   //
1012   //    detector  Detector # (1, 2, or 3)                      
1013   //    ring      Ring ID ('I' or 'O')
1014   //    sector    Sector # (For inner/outer rings: 0-19/0-39)
1015   //    strip     Strip # (For inner/outer rings: 0-511/0-255)
1016   //    edep      Total energy deposited
1017   //    count1    ADC count (a 10-bit word)
1018   //    count2    ADC count (a 10-bit word), or -1 if not used
1019   //    count3    ADC count (a 10-bit word), or -1 if not used
1020   //
1021   TClonesArray& a = *(SDigitsArray());
1022   
1023   new (a[fNsdigits++]) 
1024     AliFMDSDigit(detector, ring, sector, strip, edep, count1, count2, count3);
1025 }
1026
1027 //____________________________________________________________________
1028 void 
1029 AliFMD::ResetSDigits()
1030 {
1031   // Reset number of digits and the digits array for this detector. 
1032   //
1033   fNsdigits   = 0;
1034   if (fSDigits) fSDigits->Clear();
1035 }
1036
1037
1038 //____________________________________________________________________
1039 TClonesArray*
1040 AliFMD::HitsArray() 
1041 {
1042   // Initialize hit array if not already, and return pointer to it. 
1043   if (!fHits) { 
1044     fHits = new TClonesArray("AliFMDHit", 1000);
1045     fNhits = 0;
1046   }
1047   return fHits;
1048 }
1049
1050 //____________________________________________________________________
1051 TClonesArray*
1052 AliFMD::DigitsArray() 
1053 {
1054   // Initialize digit array if not already, and return pointer to it. 
1055   if (!fDigits) { 
1056     fDigits = new TClonesArray("AliFMDDigit", 1000);
1057     fNdigits = 0;
1058   }
1059   return fDigits;
1060 }
1061
1062 //____________________________________________________________________
1063 TClonesArray*
1064 AliFMD::SDigitsArray() 
1065 {
1066   // Initialize digit array if not already, and return pointer to it. 
1067   if (!fSDigits) { 
1068     fSDigits = new TClonesArray("AliFMDSDigit", 1000);
1069     fNsdigits = 0;
1070   }
1071   return fSDigits;
1072 }
1073
1074 //====================================================================
1075 //
1076 // Digitization 
1077 //
1078 //____________________________________________________________________
1079 void 
1080 AliFMD::Hits2Digits() 
1081 {
1082   // Create AliFMDDigit's from AliFMDHit's.  This is done by making a
1083   // AliFMDDigitizer, and executing that code.
1084   // 
1085   Warning("Hits2Digits", "Try not to use this method.\n"
1086           "Instead, use AliSimulator");
1087   AliRunDigitizer* manager = new AliRunDigitizer(1, 1);
1088   manager->SetInputStream(0, "galice.root");
1089   manager->SetOutputFile("H2Dfile");
1090   
1091   /* AliDigitizer* dig =*/ CreateDigitizer(manager);
1092   manager->Exec("");
1093   delete manager;
1094 }
1095
1096 //____________________________________________________________________
1097 void 
1098 AliFMD::Hits2SDigits() 
1099 {
1100   // Create AliFMDSDigit's from AliFMDHit's.  This is done by creating
1101   // an AliFMDSDigitizer object, and executing it. 
1102   // 
1103   AliFMDSDigitizer* digitizer = new AliFMDSDigitizer("galice.root");
1104   digitizer->Exec("");
1105   delete digitizer;
1106 }
1107
1108   
1109 //____________________________________________________________________
1110 AliDigitizer* 
1111 AliFMD::CreateDigitizer(AliRunDigitizer* manager) const
1112 {
1113   // Create a digitizer object 
1114   AliFMDDigitizer* digitizer = new AliFMDDigitizer(manager);
1115   return digitizer;
1116 }
1117
1118 //====================================================================
1119 //
1120 // Raw data simulation 
1121 //
1122 //__________________________________________________________________
1123 void 
1124 AliFMD::Digits2Raw() 
1125 {
1126   // Turn digits into raw data. 
1127   // 
1128   // This uses the class AliFMDRawWriter to do the job.   Please refer
1129   // to that class for more information. 
1130   AliFMDRawWriter writer(this);
1131   writer.Exec();
1132 }
1133
1134
1135 //====================================================================
1136 //
1137 // Utility 
1138 //
1139 //__________________________________________________________________
1140 void 
1141 AliFMD::Browse(TBrowser* b) 
1142 {
1143   // Browse this object. 
1144   //
1145   AliDebug(30, "\tBrowsing the FMD");
1146   AliDetector::Browse(b);
1147   b->Add(AliFMDGeometry::Instance());
1148 }
1149
1150 //___________________________________________________________________
1151 //
1152 // EOF
1153 //