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