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