Added the AliTrackReference::kFMD to enable the TrackRefs to be stored
[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 <TRotMatrix.h>         // ROOT_TRotMatrix
90 #include <TTree.h>              // ROOT_TTree
91 #include <TBrowser.h>           // ROOT_TBrowser
92 #include <TVirtualMC.h>         // ROOT_TVirtualMC
93 #include <TVector2.h>           // ROOT_TVector2 
94 #include <TGeoManager.h>        // ROOT_TGeoManager
95
96 #include <AliRunDigitizer.h>    // ALIRUNDIGITIZER_H
97 #include <AliLoader.h>          // ALILOADER_H
98 #include <AliRun.h>             // ALIRUN_H
99 #include <AliMC.h>              // ALIMC_H
100 #include <AliMagF.h>            // ALIMAGF_H
101 // #include <AliLog.h>          // ALILOG_H
102 #include "AliFMDDebug.h" // Better debug macros
103 #include "AliFMD.h"             // ALIFMD_H
104 #include "AliFMDDigit.h"        // ALIFMDDIGIT_H
105 #include "AliFMDSDigit.h"       // ALIFMDSDIGIT_H
106 #include "AliFMDHit.h"          // ALIFMDHIT_H
107 #include "AliFMDGeometry.h"     // ALIFMDGEOMETRY_H
108 #include "AliFMDDetector.h"     // ALIFMDDETECTOR_H
109 #include "AliFMDRing.h"         // ALIFMDRING_H
110 #include "AliFMDDigitizer.h"    // ALIFMDDIGITIZER_H
111 #include "AliFMDHitDigitizer.h" // ALIFMDSDIGITIZER_H
112 // #define USE_SSDIGITIZER 
113 //#ifdef USE_SSDIGITIZER
114 //# include "AliFMDSSDigitizer.h"       // ALIFMDSDIGITIZER_H
115 //#endif
116 // #include "AliFMDGeometryBuilder.h"
117 #include "AliFMDRawWriter.h"    // ALIFMDRAWWRITER_H
118 #include "AliTrackReference.h" 
119 #include "AliFMDStripIndex.h"
120
121 //____________________________________________________________________
122 ClassImp(AliFMD)
123 #if 0
124   ; // This is to keep Emacs from indenting the next line 
125 #endif 
126
127 //____________________________________________________________________
128 AliFMD::AliFMD()
129   : AliDetector(),
130     fSDigits(0), 
131     fNsdigits(0),
132     fDetailed(kTRUE),
133     fUseOld(kFALSE),
134     fUseAssembly(kTRUE),
135     fBad(0) 
136 {
137   //
138   // Default constructor for class AliFMD
139   //
140   AliFMDDebug(10, ("\tDefault CTOR"));
141   fHits        = 0;
142   fDigits      = 0;
143   fIshunt      = 0;
144   fBad         = new TClonesArray("AliFMDHit");
145 }
146
147 //____________________________________________________________________
148 AliFMD::AliFMD(const char *name, const char *title)
149   : AliDetector (name, title),
150     fSDigits(0),
151     fNsdigits(0),
152     fDetailed(kTRUE),
153     fUseOld(kFALSE),
154     fUseAssembly(kFALSE),
155     fBad(0)
156 {
157   //
158   // Standard constructor for Forward Multiplicity Detector
159   //
160   AliFMDDebug(10, ("\tStandard CTOR"));
161   fBad         = new TClonesArray("AliFMDHit");
162   
163   // Initialise Hit array
164   HitsArray();
165   gAlice->GetMCApp()->AddHitList(fHits);
166
167   // (S)Digits for the detectors disk
168   DigitsArray();
169   SDigitsArray();
170   
171   // CHC: What is this?
172   fIshunt = 0;
173   //PH  SetMarkerColor(kRed);
174   //PH  SetLineColor(kYellow);
175 }
176
177 //____________________________________________________________________
178 AliFMD::~AliFMD ()
179 {
180   // Destructor for base class AliFMD
181   if (fHits) {
182     fHits->Delete();
183     delete fHits;
184     fHits = 0;
185   }
186   if (fDigits) {
187     fDigits->Delete();
188     delete fDigits;
189     fDigits = 0;
190   }
191   if (fSDigits) {
192     fSDigits->Delete();
193     delete fSDigits;
194     fSDigits = 0;
195   }
196   if (fBad) {
197     fBad->Delete();
198     delete fBad;
199     fBad = 0;
200   }
201 }
202
203
204 //====================================================================
205 //
206 // GEometry ANd Traking
207 //
208 //____________________________________________________________________
209 void 
210 AliFMD::CreateGeometry()
211 {
212   //
213   // Create the geometry of Forward Multiplicity Detector.  The actual
214   // construction of the geometry is delegated to the class
215   // AliFMDGeometryBuilder, invoked by the singleton manager
216   // AliFMDGeometry. 
217   //
218   AliFMDGeometry*  fmd = AliFMDGeometry::Instance();
219   fmd->SetDetailed(fDetailed);
220   fmd->UseAssembly(fUseAssembly);
221   fmd->Build();
222 }    
223
224 //____________________________________________________________________
225 void AliFMD::CreateMaterials() 
226 {
227   // Define the materials and tracking mediums needed by the FMD
228   // simulation.   These mediums are made by sending the messages
229   // AliMaterial, AliMixture, and AliMedium to the passed AliModule
230   // object module.   The defined mediums are 
231   // 
232   //    FMD Si$         Silicon (active medium in sensors)
233   //    FMD C$          Carbon fibre (support cone for FMD3 and vacuum pipe)
234   //    FMD Al$         Aluminium (honeycomb support plates)
235   //    FMD PCB$        Printed Circuit Board (FEE board with VA1_3)
236   //    FMD Chip$       Electronics chips (currently not used)
237   //    FMD Air$        Air (Air in the FMD)
238   //    FMD Plastic$    Plastic (Support legs for the hybrid cards)
239   //
240   // The geometry builder should really be the one that creates the
241   // materials, but the architecture of AliROOT makes that design
242   // akward.  What should happen, was that the AliFMDGeometryBuilder
243   // made the mediums, and that this class retrives pointers from the
244   // TGeoManager, and registers the mediums here.  Alas, it's not
245   // really that easy. 
246   //
247   AliFMDDebug(10, ("\tCreating materials"));
248   // Get pointer to geometry singleton object. 
249   AliFMDGeometry* geometry = AliFMDGeometry::Instance();
250   geometry->Init();
251 #if 0
252   if (gGeoManager && gGeoManager->GetMedium("FMD Si$")) {
253     // We need to figure out the some stuff about the geometry
254     fmd->ExtractGeomInfo();
255     return;
256   }
257 #endif  
258   Int_t    id;
259   Double_t a                = 0;
260   Double_t z                = 0;
261   Double_t density          = 0;
262   Double_t radiationLength  = 0;
263   Double_t absorbtionLength = 999;
264   Int_t    fieldType        = gAlice->Field()->Integ();     // Field type 
265   Double_t maxField         = gAlice->Field()->Max();     // Field max.
266   Double_t maxBending       = 0;     // Max Angle
267   Double_t maxStepSize      = 0.001; // Max step size 
268   Double_t maxEnergyLoss    = 1;     // Max Delta E
269   Double_t precision        = 0.001; // Precision
270   Double_t minStepSize      = 0.001; // Minimum step size 
271  
272   // Silicon 
273   a                = 28.0855;
274   z                = 14.;
275   density          = geometry->GetSiDensity();
276   radiationLength  = 9.36;
277   maxBending       = 1;
278   maxStepSize      = .001;
279   precision        = .001;
280   minStepSize      = .001;
281   id               = kSiId;
282   AliMaterial(id, "Si$", a, z, density, radiationLength, absorbtionLength);
283   AliMedium(kSiId, "Si$", id,1,fieldType,maxField,maxBending,
284             maxStepSize,maxEnergyLoss,precision,minStepSize);
285   
286
287   // Carbon 
288   a                = 12.011;
289   z                = 6.;
290   density          = 2.265;
291   radiationLength  = 18.8;
292   maxBending       = 10;
293   maxStepSize      = .01;
294   precision        = .003;
295   minStepSize      = .003;
296   id               = kCarbonId;
297   AliMaterial(id, "Carbon$", a, z, density, radiationLength, absorbtionLength);
298   AliMedium(kCarbonId, "Carbon$", id,0,fieldType,maxField,maxBending,
299                     maxStepSize,maxEnergyLoss,precision,minStepSize);
300
301   // Aluminum
302   a                = 26.981539;
303   z                = 13.;
304   density          = 2.7;
305   radiationLength  = 8.9;
306   id               = kAlId;
307   AliMaterial(id, "Aluminum$",a,z, density, radiationLength, absorbtionLength);
308   AliMedium(kAlId, "Aluminum$", id, 0, fieldType, maxField, maxBending,
309             maxStepSize, maxEnergyLoss, precision, minStepSize);
310   
311   
312   // Copper 
313   a                = 63.546;
314   z                = 29;
315   density          =  8.96;
316   radiationLength  =  1.43;
317   id               = kCopperId;
318   AliMaterial(id, "Copper$", 
319                       a, z, density, radiationLength, absorbtionLength);
320   AliMedium(kCopperId, "Copper$", id, 0, fieldType, maxField, maxBending,
321             maxStepSize, maxEnergyLoss, precision, minStepSize);
322   
323
324   // Silicon chip 
325   {
326     Float_t as[] = { 12.0107,      14.0067,      15.9994,
327                       1.00794,     28.0855,     107.8682 };
328     Float_t zs[] = {  6.,           7.,           8.,
329                       1.,          14.,          47. };
330     Float_t ws[] = {  0.039730642,  0.001396798,  0.01169634,
331                       0.004367771,  0.844665,     0.09814344903 };
332     density          = 2.36436;
333     maxBending       = 10;
334     maxStepSize      = .01;
335     precision        = .003;
336     minStepSize      = .003;
337     id               = kSiChipId;
338     AliMixture(id, "Si Chip$", as, zs, density, 6, ws);
339     AliMedium(kSiChipId, "Si Chip$",  id, 0, fieldType, maxField, maxBending, 
340               maxStepSize, maxEnergyLoss, precision, minStepSize);
341   }
342   
343   // Kaption
344   {
345     Float_t as[] = { 1.00794,  12.0107,  14.010,   15.9994};
346     Float_t zs[] = { 1.,        6.,       7.,       8.};
347     Float_t ws[] = { 0.026362,  0.69113,  0.07327,  0.209235};
348     density          = 1.42;
349     maxBending       = 1;
350     maxStepSize      = .001;
351     precision        = .001;
352     minStepSize      = .001;
353     id               = kKaptonId;
354     AliMixture(id, "Kaption$", as, zs, density, 4, ws);
355     AliMedium(kKaptonId, "Kaption$", id,0,fieldType,maxField,maxBending,
356               maxStepSize,maxEnergyLoss,precision,minStepSize);
357   }
358
359   // Air
360   {
361     Float_t as[] = { 12.0107, 14.0067,   15.9994,  39.948 };
362     Float_t zs[] = {  6.,      7.,       8.,       18. };
363     Float_t ws[] = { 0.000124, 0.755267, 0.231781, 0.012827 }; 
364     density      = .00120479;
365     maxBending   = 1;
366     maxStepSize  = .001;
367     precision    = .001;
368     minStepSize  = .001;
369     id           = kAirId;
370     AliMixture(id, "Air$", as, zs, density, 4, ws);
371     AliMedium(kAirId, "Air$", id,0,fieldType,maxField,maxBending,
372               maxStepSize,maxEnergyLoss,precision,minStepSize);
373   }
374   
375   // PCB
376   {
377     Float_t zs[] = { 14.,         20.,         13.,         12.,
378                       5.,         22.,         11.,         19.,
379                      26.,          9.,          8.,          6.,
380                       7.,          1.};
381     Float_t as[] = { 28.0855,     40.078,      26.981538,   24.305, 
382                      10.811,      47.867,      22.98977,    39.0983,
383                      55.845,      18.9984,     15.9994,     12.0107,
384                      14.0067,      1.00794};
385     Float_t ws[] = {  0.15144894,  0.08147477,  0.04128158,  0.00904554, 
386                       0.01397570,  0.00287685,  0.00445114,  0.00498089,
387                       0.00209828,  0.00420000,  0.36043788,  0.27529426,
388                       0.01415852,  0.03427566};
389     density      = 1.8;
390     maxBending   = 1;
391     maxStepSize  = .001;
392     precision    = .001;
393     minStepSize  = .001;
394     id           = kPcbId;
395     AliMixture(id, "PCB$", as, zs, density, 14, ws);
396     AliMedium(kPcbId, "PCB$", id,0,fieldType,maxField,maxBending,
397               maxStepSize,maxEnergyLoss,precision,minStepSize);
398   }
399   
400   // Stainless steel
401   {
402     Float_t as[] = { 55.847, 51.9961, 58.6934, 28.0855 };
403     Float_t zs[] = { 26.,    24.,     28.,     14.     };
404     Float_t ws[] = { .715,   .18,     .1,      .005    };
405     density      = 7.88;
406     id           = kSteelId;
407     AliMixture(id, "Steel$", as, zs, density, 4, ws);
408     AliMedium(kSteelId, "Steel$", id, 0, fieldType, maxField, maxBending, 
409               maxStepSize, maxEnergyLoss, precision, minStepSize);
410   }
411   // Plastic 
412   {
413     Float_t as[] = { 1.01, 12.01 };
414     Float_t zs[] = { 1.,   6.    };
415     Float_t ws[] = { 1.,   1.    };
416     density      = 1.03;
417     maxBending   = 10;
418     maxStepSize  = .01;
419     precision    = .003;
420     minStepSize  = .003;
421     id           = kPlasticId;
422     AliMixture(id, "Plastic$", as, zs, density, -2, ws);
423     AliMedium(kPlasticId, "Plastic$", id,0,fieldType,maxField,maxBending,
424               maxStepSize,maxEnergyLoss,precision,minStepSize);
425   }
426
427 }
428
429 //____________________________________________________________________
430 void  
431 AliFMD::Init()
432 {
433   // Initialize the detector 
434   // 
435   AliFMDDebug(1, ("Initialising FMD detector object"));
436   TVirtualMC*      mc     = TVirtualMC::GetMC();
437   AliFMDGeometry*  fmd    = AliFMDGeometry::Instance();
438   const TArrayI&   actGeo = fmd->ActiveIds();
439   TArrayI          actVmc(actGeo.fN);
440   for (Int_t i = 0; i < actGeo.fN; i++) {
441     TGeoVolume *sens = gGeoManager->GetVolume(actGeo[i]);
442     if (!sens) {
443       AliError(Form("No TGeo volume for sensitive volume ID=%d",actGeo[i]));
444       continue;
445     }   
446     actVmc[i] = mc->VolId(sens->GetName());
447     AliFMDDebug(1, ("Active vol id # %d: %d changed to %d", 
448                     i, actGeo[i], actVmc[i]));
449   }
450   fmd->SetActive(actVmc.fArray, actVmc.fN);
451   // fmd->InitTransformations();
452 }
453
454 //____________________________________________________________________
455 void
456 AliFMD::FinishEvent()
457 {
458   // Called at the end of the an event in simulations.  If the debug
459   // level is high enough, then the `bad' hits are printed.
460   // 
461   if (AliLog::GetDebugLevel("FMD", "AliFMD") < 10) return;
462   if (fBad && fBad->GetEntries() > 0) {
463     AliWarning((Form("EndEvent", "got %d 'bad' hits", fBad->GetEntries())));
464     TIter next(fBad);
465     AliFMDHit* hit;
466     while ((hit = static_cast<AliFMDHit*>(next()))) hit->Print("D");
467     fBad->Clear();
468   }
469 }
470
471
472
473 //====================================================================
474 //
475 // Hit and Digit managment 
476 //
477 //____________________________________________________________________
478 void 
479 AliFMD::MakeBranch(Option_t * option)
480 {
481   // Create Tree branches for the FMD.
482   //
483   // Options:
484   //
485   //    H          Make a branch of TClonesArray of AliFMDHit's
486   //    D          Make a branch of TClonesArray of AliFMDDigit's
487   //    S          Make a branch of TClonesArray of AliFMDSDigit's
488   // 
489   const Int_t kBufferSize = 16000;
490   TString branchname(GetName());
491   TString opt(option);
492   
493   if (opt.Contains("H", TString::kIgnoreCase)) {
494     HitsArray();
495     AliDetector::MakeBranch(option); 
496   }
497   if (opt.Contains("D", TString::kIgnoreCase)) { 
498     DigitsArray();
499     MakeBranchInTree(fLoader->TreeD(), branchname.Data(),
500                      &fDigits, kBufferSize, 0);
501   }
502   if (opt.Contains("S", TString::kIgnoreCase)) { 
503     SDigitsArray();
504     MakeBranchInTree(fLoader->TreeS(), branchname.Data(),
505                      &fSDigits, kBufferSize, 0);
506   }
507 }
508
509 //____________________________________________________________________
510 void 
511 AliFMD::SetTreeAddress()
512 {
513   // Set branch address for the Hits, Digits, and SDigits Tree.
514   if (fLoader->TreeH()) HitsArray();
515   AliDetector::SetTreeAddress();
516
517   TTree *treeD = fLoader->TreeD();
518   if (treeD) {
519     DigitsArray();
520     TBranch* branch = treeD->GetBranch ("FMD");
521     if (branch) branch->SetAddress(&fDigits);
522   }
523
524   TTree *treeS = fLoader->TreeS();
525   if (treeS) {
526     SDigitsArray();
527     TBranch* branch = treeS->GetBranch ("FMD");
528     if (branch) branch->SetAddress(&fSDigits);
529   }
530 }
531
532 //____________________________________________________________________
533 void 
534 AliFMD::SetHitsAddressBranch(TBranch *b)
535 {
536   // Set the TClonesArray to read hits into. 
537   b->SetAddress(&fHits);
538 }
539 //____________________________________________________________________
540 void 
541 AliFMD::SetSDigitsAddressBranch(TBranch *b)
542 {
543   // Set the TClonesArray to read hits into. 
544   b->SetAddress(&fSDigits);
545 }
546
547 //____________________________________________________________________
548 void 
549 AliFMD::AddHit(Int_t track, Int_t *vol, Float_t *hits) 
550 {
551   // Add a hit to the hits tree 
552   // 
553   // The information of the two arrays are decoded as 
554   // 
555   // Parameters
556   //    track                Track #
557   //    ivol[0]  [UShort_t ] Detector # 
558   //    ivol[1]  [Char_t   ] Ring ID 
559   //    ivol[2]  [UShort_t ] Sector #
560   //    ivol[3]  [UShort_t ] Strip # 
561   //    hits[0]  [Float_t  ] Track's X-coordinate at hit 
562   //    hits[1]  [Float_t  ] Track's Y-coordinate at hit
563   //    hits[3]  [Float_t  ] Track's Z-coordinate at hit
564   //    hits[4]  [Float_t  ] X-component of track's momentum             
565   //    hits[5]  [Float_t  ] Y-component of track's momentum             
566   //    hits[6]  [Float_t  ] Z-component of track's momentum            
567   //    hits[7]  [Float_t  ] Energy deposited by track                  
568   //    hits[8]  [Int_t    ] Track's particle Id # 
569   //    hits[9]  [Float_t  ] Time when the track hit
570   // 
571   // 
572   AddHitByFields(track, 
573                  UShort_t(vol[0]),  // Detector # 
574                  Char_t(vol[1]),    // Ring ID
575                  UShort_t(vol[2]),  // Sector # 
576                  UShort_t(vol[3]),  // Strip # 
577                  hits[0],           // X
578                  hits[1],           // Y
579                  hits[2],           // Z
580                  hits[3],           // Px
581                  hits[4],           // Py
582                  hits[5],           // Pz
583                  hits[6],           // Energy loss 
584                  Int_t(hits[7]),    // PDG 
585                  hits[8]);          // Time
586 }
587
588 //____________________________________________________________________
589 AliFMDHit*
590 AliFMD::AddHitByFields(Int_t    track, 
591                        UShort_t detector, 
592                        Char_t   ring, 
593                        UShort_t sector, 
594                        UShort_t strip, 
595                        Float_t  x, 
596                        Float_t  y, 
597                        Float_t  z,
598                        Float_t  px, 
599                        Float_t  py, 
600                        Float_t  pz,
601                        Float_t  edep,
602                        Int_t    pdg,
603                        Float_t  t, 
604                        Float_t  l, 
605                        Bool_t   stop)
606 {
607   // Add a hit to the list
608   //
609   // Parameters:
610   // 
611   //    track     Track #
612   //    detector  Detector # (1, 2, or 3)                      
613   //    ring      Ring ID ('I' or 'O')
614   //    sector    Sector # (For inner/outer rings: 0-19/0-39)
615   //    strip     Strip # (For inner/outer rings: 0-511/0-255)
616   //    x         Track's X-coordinate at hit
617   //    y         Track's Y-coordinate at hit
618   //    z         Track's Z-coordinate at hit
619   //    px        X-component of track's momentum 
620   //    py        Y-component of track's momentum
621   //    pz        Z-component of track's momentum
622   //    edep      Energy deposited by track
623   //    pdg       Track's particle Id #
624   //    t         Time when the track hit 
625   //    l         Track length through the material. 
626   //    stop      Whether track was stopped or disappeared
627   // 
628   TClonesArray& a = *(HitsArray());
629   // Search through the list of already registered hits, and see if we
630   // find a hit with the same parameters.  If we do, then don't create
631   // a new hit, but rather update the energy deposited in the hit.
632   // This is done, so that a FLUKA based simulation will get the
633   // number of hits right, not just the enerrgy deposition. 
634   AliFMDHit* hit = 0;
635   for (Int_t i = 0; i < fNhits; i++) {
636     if (!a.At(i)) continue;
637     hit = static_cast<AliFMDHit*>(a.At(i));
638     if (hit->Detector() == detector 
639         && hit->Ring() == ring
640         && hit->Sector() == sector 
641         && hit->Strip() == strip
642         && hit->Track() == track) {
643       AliFMDDebug(1, ("already had a hit in FMD%d%c[%2d,%3d] for track # %d,"
644                        " adding energy (%f) to that hit (%f) -> %f", 
645                        detector, ring, sector, strip, track, edep, hit->Edep(),
646                        hit->Edep() + edep));
647       hit->SetEdep(hit->Edep() + edep);
648       return hit;
649     }
650   }
651   // If hit wasn't already registered, do so know. 
652   hit = new (a[fNhits]) AliFMDHit(fIshunt, track, detector, ring, sector, 
653                                   strip, x, y, z, px, py, pz, edep, pdg, t, 
654                                   l, stop);
655   fNhits++;
656   
657   //Reference track
658
659   AliMC *mcApplication = (AliMC*)gAlice->GetMCApp();
660   
661   AliTrackReference* trackRef = AddTrackReference(mcApplication->GetCurrentTrackNumber(), AliTrackReference::kFMD); 
662   UInt_t stripId = AliFMDStripIndex::Pack(detector,ring,sector,strip);
663   trackRef->SetUserId(stripId);
664   
665   
666   
667   return hit;
668 }
669
670 //____________________________________________________________________
671 void 
672 AliFMD::AddDigit(Int_t* digits, Int_t*)
673 {
674   // Add a digit to the Digit tree 
675   // 
676   // Paramters 
677   //
678   //    digits[0]  [UShort_t] Detector #
679   //    digits[1]  [Char_t]   Ring ID
680   //    digits[2]  [UShort_t] Sector #
681   //    digits[3]  [UShort_t] Strip #
682   //    digits[4]  [UShort_t] ADC Count 
683   //    digits[5]  [Short_t]  ADC Count, -1 if not used
684   //    digits[6]  [Short_t]  ADC Count, -1 if not used 
685   // 
686   AddDigitByFields(UShort_t(digits[0]),  // Detector #
687                    Char_t(digits[1]),    // Ring ID
688                    UShort_t(digits[2]),  // Sector #
689                    UShort_t(digits[3]),  // Strip #
690                    UShort_t(digits[4]),  // ADC Count1 
691                    Short_t(digits[5]),   // ADC Count2 
692                    Short_t(digits[6]),   // ADC Count3 
693                    Short_t(digits[7])); 
694 }
695
696 //____________________________________________________________________
697 void 
698 AliFMD::AddDigitByFields(UShort_t detector, 
699                          Char_t   ring, 
700                          UShort_t sector, 
701                          UShort_t strip, 
702                          UShort_t count1, 
703                          Short_t  count2,
704                          Short_t  count3, 
705                          Short_t  count4)
706 {
707   // add a real digit - as coming from data
708   // 
709   // Parameters 
710   //
711   //    detector  Detector # (1, 2, or 3)                      
712   //    ring      Ring ID ('I' or 'O')
713   //    sector    Sector # (For inner/outer rings: 0-19/0-39)
714   //    strip     Strip # (For inner/outer rings: 0-511/0-255)
715   //    count1    ADC count (a 10-bit word)
716   //    count2    ADC count (a 10-bit word), or -1 if not used
717   //    count3    ADC count (a 10-bit word), or -1 if not used
718   TClonesArray& a = *(DigitsArray());
719   
720   new (a[fNdigits++]) 
721     AliFMDDigit(detector, ring, sector, strip, count1, count2, count3, count4);
722   AliFMDDebug(15, ("Adding digit # %5d/%5d for FMD%d%c[%2d,%3d]=(%d,%d,%d,%d)",
723                    fNdigits-1, a.GetEntriesFast(),
724                    detector, ring, sector, strip, 
725                    count1, count2, count3, count4));
726   
727 }
728
729 //____________________________________________________________________
730 void 
731 AliFMD::AddSDigit(Int_t* digits)
732 {
733   // Add a digit to the SDigit tree 
734   // 
735   // Paramters 
736   //
737   //    digits[0]  [UShort_t] Detector #
738   //    digits[1]  [Char_t]   Ring ID
739   //    digits[2]  [UShort_t] Sector #
740   //    digits[3]  [UShort_t] Strip #
741   //    digits[4]  [Float_t]  Total energy deposited 
742   //    digits[5]  [UShort_t] ADC Count 
743   //    digits[6]  [Short_t]  ADC Count, -1 if not used
744   //    digits[7]  [Short_t]  ADC Count, -1 if not used 
745   // 
746   AddSDigitByFields(UShort_t(digits[0]),   // Detector #
747                     Char_t(digits[1]),     // Ring ID
748                     UShort_t(digits[2]),   // Sector #
749                     UShort_t(digits[3]),   // Strip #
750                     Float_t(digits[4]),    // Edep
751                     UShort_t(digits[5]),   // ADC Count1 
752                     Short_t(digits[6]),    // ADC Count2 
753                     Short_t(digits[7]),    // ADC Count3 
754                     Short_t(digits[8]),    // ADC Count4
755                     UShort_t(digits[9]),   // N particles
756                     UShort_t(digits[10])); // N primaries
757    
758 }
759
760 //____________________________________________________________________
761 void 
762 AliFMD::AddSDigitByFields(UShort_t detector, 
763                           Char_t   ring, 
764                           UShort_t sector, 
765                           UShort_t strip, 
766                           Float_t  edep,
767                           UShort_t count1, 
768                           Short_t  count2,
769                           Short_t  count3, 
770                           Short_t  count4, 
771                           UShort_t ntot, 
772                           UShort_t nprim)
773 {
774   // add a summable digit
775   // 
776   // Parameters 
777   //
778   //    detector  Detector # (1, 2, or 3)                      
779   //    ring      Ring ID ('I' or 'O')
780   //    sector    Sector # (For inner/outer rings: 0-19/0-39)
781   //    strip     Strip # (For inner/outer rings: 0-511/0-255)
782   //    edep      Total energy deposited
783   //    count1    ADC count (a 10-bit word)
784   //    count2    ADC count (a 10-bit word), or -1 if not used
785   //    count3    ADC count (a 10-bit word), or -1 if not used
786   //
787   TClonesArray& a = *(SDigitsArray());
788   // AliFMDDebug(0, ("Adding sdigit # %d", fNsdigits));
789   
790   new (a[fNsdigits++]) 
791     AliFMDSDigit(detector, ring, sector, strip, edep, 
792                  count1, count2, count3, count4, ntot, nprim);
793 }
794
795 //____________________________________________________________________
796 void 
797 AliFMD::ResetSDigits()
798 {
799   // Reset number of digits and the digits array for this detector. 
800   //
801   fNsdigits   = 0;
802   if (fSDigits) fSDigits->Clear();
803 }
804
805
806 //____________________________________________________________________
807 TClonesArray*
808 AliFMD::HitsArray() 
809 {
810   // Initialize hit array if not already, and return pointer to it. 
811   if (!fHits) { 
812     fHits = new TClonesArray("AliFMDHit", 1000);
813     fNhits = 0;
814   }
815   return fHits;
816 }
817
818 //____________________________________________________________________
819 TClonesArray*
820 AliFMD::DigitsArray() 
821 {
822   // Initialize digit array if not already, and return pointer to it. 
823   if (!fDigits) { 
824     fDigits = new TClonesArray("AliFMDDigit", 1000);
825     fNdigits = 0;
826   }
827   return fDigits;
828 }
829
830 //____________________________________________________________________
831 TClonesArray*
832 AliFMD::SDigitsArray() 
833 {
834   // Initialize digit array if not already, and return pointer to it. 
835   if (!fSDigits) { 
836     fSDigits = new TClonesArray("AliFMDSDigit", 1000);
837     fNsdigits = 0;
838   }
839   return fSDigits;
840 }
841
842 //====================================================================
843 //
844 // Digitization 
845 //
846 //____________________________________________________________________
847 void 
848 AliFMD::Hits2Digits() 
849 {
850   // Create AliFMDDigit's from AliFMDHit's.  This is done by making a
851   // AliFMDDigitizer, and executing that code.
852   // 
853   AliFMDHitDigitizer digitizer(this, AliFMDHitDigitizer::kDigits);
854   digitizer.Init();
855   digitizer.Exec("");
856 }
857
858 //____________________________________________________________________
859 void 
860 AliFMD::Hits2SDigits() 
861 {
862   // Create AliFMDSDigit's from AliFMDHit's.  This is done by creating
863   // an AliFMDSDigitizer object, and executing it. 
864   // 
865   AliFMDHitDigitizer digitizer(this, AliFMDHitDigitizer::kSDigits);
866   digitizer.Init();
867   digitizer.Exec("");
868 }
869
870   
871 //____________________________________________________________________
872 AliDigitizer* 
873 AliFMD::CreateDigitizer(AliRunDigitizer* manager) const
874 {
875   // Create a digitizer object 
876   
877   /* This is what we probably _should_ do */
878   AliFMDBaseDigitizer* digitizer = 0;
879   
880 #ifdef USE_SSDIGITIZER
881   digitizer = new AliFMDSSDigitizer(manager);
882 #else 
883   /* This is what we actually do, and will work */
884 #if 0
885   AliInfo("SDigit->Digit conversion not really supported, "
886           "doing Hit->Digit conversion instead");
887 #endif
888   digitizer = new AliFMDDigitizer(manager);
889 #endif
890   return digitizer;
891 }
892
893 //====================================================================
894 //
895 // Raw data simulation 
896 //
897 //__________________________________________________________________
898 void 
899 AliFMD::Digits2Raw() 
900 {
901   // Turn digits into raw data. 
902   // 
903   // This uses the class AliFMDRawWriter to do the job.   Please refer
904   // to that class for more information. 
905   AliFMDRawWriter writer(this);
906   writer.Exec();
907 }
908
909
910 //====================================================================
911 //
912 // Utility 
913 //
914 //__________________________________________________________________
915 void 
916 AliFMD::Browse(TBrowser* b) 
917 {
918   // Browse this object. 
919   //
920   AliFMDDebug(30, ("\tBrowsing the FMD"));
921   AliDetector::Browse(b);
922   b->Add(AliFMDGeometry::Instance());
923 }
924
925 //____________________________________________________________________  
926 void
927 AliFMD::AddAlignableVolumes() const
928 {
929   //
930   // Create entries for alignable volumes associating the symbolic volume
931   // name with the corresponding volume path. Needs to be syncronized with
932   // eventual changes in the geometry.
933   // 
934   // This code was made by Raffaele Grosso <rgrosso@mail.cern.ch>.  I
935   // (cholm) will probably want to change it.   For one, I think it
936   // should be the job of the geometry manager to deal with this. 
937   AliInfo("Add FMD alignable volumes");
938   AliFMDGeometry::Instance()->SetAlignableVolumes();
939 #if 0  
940   for(size_t f = 1; f <= 3; f++){ // Detector 1,2,3
941     for(size_t tb =  0; tb <2 ; tb++){ // Top/Bottom 
942       char     stb = tb == 0 ? 'T' : 'B';
943       unsigned min = tb == 0 ? 0   : 5;
944
945       TString halfVol(Form("/ALIC_1/F%dM%c_%d", f, stb, f));
946       TString halfSym(halfVol);
947       if(!gGeoManager->SetAlignableEntry(halfSym.Data(),halfVol.Data()))
948         AliFatal(Form("Alignable entry %s not created. "
949                       "Volume path %s not valid", 
950                       halfSym.Data(),halfVol.Data()));
951       for(size_t io = 0; io < 2; io++){ // inner, outer 
952         if (f==1 && io==1) continue; // Only one ring in FMD1 
953         if(tb == 1 && io==1) min=10;
954         char     sio = (io == 0 ? 'I' : 'O');
955         unsigned nio = (io == 0 ? 3   : 9);
956         unsigned max = (io == 0 ? 5   : 10) + min;
957         
958         for(size_t i = min; i < max; i++) { // Modules
959           TString modVol(Form("%s/F%c%cV_7%d/F%cSE_%d", halfVol.Data(), 
960                               sio, stb, nio, sio, i));
961           TString modSym(modVol);
962           if(!gGeoManager->SetAlignableEntry(modSym.Data(),modVol.Data()))
963             AliFatal(Form("Alignable entry %s not created. "
964                           "Volume path %s not valid", 
965                           modSym.Data(), modVol.Data()));
966         }
967       }
968     }
969   }
970 #endif
971 }
972 //___________________________________________________________________
973 //
974 // EOF
975 //