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