Increased error checking and possibility to extract some diagnostics
[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 #if 0
434 //____________________________________________________________________
435 void  
436 AliFMD::SetTrackingParameters(Int_t imed, 
437                               Float_t gamma,                 
438                               Float_t electron, 
439                               Float_t neutral_hadron, 
440                               Float_t charged_hadron, 
441                               Float_t muon,
442                               Float_t electron_bremstrahlung, 
443                               Float_t muon__bremstrahlung, 
444                               Float_t electron_delta,
445                               Float_t muon_delta,
446                               Float_t muon_pair,
447                               Int_t   annihilation, 
448                               Int_t   bremstrahlung, 
449                               Int_t   compton_scattering, 
450                               Int_t   decay,
451                               Int_t   delta_ray, 
452                               Int_t   hadronic, 
453                               Int_t   energy_loss, 
454                               Int_t   multiple_scattering, 
455                               Int_t   pair_production, 
456                               Int_t   photon_production, 
457                               Int_t   rayleigh_scattering)
458 {
459   // Disabled by request of FCA, kept for reference only
460   if (!gMC) return;
461   TArrayI& idtmed = *(GetIdtmed());
462   Int_t    iimed  = idtmed[imed];
463   // gMC->Gstpar(iimed, "CUTGAM",       gamma);
464   // gMC->Gstpar(iimed, "CUTELE",       electron);
465   // gMC->Gstpar(iimed, "CUTNEU",       neutral_hadron);
466   // gMC->Gstpar(iimed, "CUTHAD",       charged_hadron);
467   // gMC->Gstpar(iimed, "CUTMUO",       muon);
468   // gMC->Gstpar(iimed, "BCUTE",        electron_bremstrahlung);
469   // gMC->Gstpar(iimed, "BCUTM",        muon__bremstrahlung);
470   // gMC->Gstpar(iimed, "DCUTE",        electron_delta);
471   // gMC->Gstpar(iimed, "DCUTM",        muon_delta);
472   // gMC->Gstpar(iimed, "PPCUTM",       muon_pair);
473   // gMC->Gstpar(iimed, "ANNI", Float_t(annihilation));
474   // gMC->Gstpar(iimed, "BREM", Float_t(bremstrahlung));
475   // gMC->Gstpar(iimed, "COMP", Float_t(compton_scattering));
476   // gMC->Gstpar(iimed, "DCAY", Float_t(decay));
477   // gMC->Gstpar(iimed, "DRAY", Float_t(delta_ray));
478   // gMC->Gstpar(iimed, "HADR", Float_t(hadronic));
479   // gMC->Gstpar(iimed, "LOSS", Float_t(energy_loss));
480   // gMC->Gstpar(iimed, "MULS", Float_t(multiple_scattering));
481   // gMC->Gstpar(iimed, "PAIR", Float_t(pair_production));
482   // gMC->Gstpar(iimed, "PHOT", Float_t(photon_production));
483   // gMC->Gstpar(iimed, "RAYL", Float_t(rayleigh_scattering));
484 }
485 #endif
486
487 //____________________________________________________________________
488 void  
489 AliFMD::Init()
490 {
491   // Initialize the detector 
492   // 
493   AliFMDDebug(1, ("Initialising FMD detector object"));
494   TVirtualMC*      mc     = TVirtualMC::GetMC();
495   AliFMDGeometry*  fmd    = AliFMDGeometry::Instance();
496   TArrayI          actGeo = fmd->ActiveIds();
497   bool             valid  = true;
498   if (actGeo.fN <= 0) valid = false;
499   else { 
500     for (int i = 0; i < actGeo.fN; i++) {
501       if (actGeo[i] < 0) { 
502         valid = false;
503         break;
504       }
505     }
506   }
507   if (!valid) { 
508     AliFMDDebug(1, ("Extracting geometry info from loaded geometry"));
509     fmd->ExtractGeomInfo();
510     actGeo = fmd->ActiveIds();
511   }
512   TArrayI          actVmc(actGeo.fN);
513   for (Int_t i = 0; i < actGeo.fN; i++) {
514     if (actGeo[i] < 0) { 
515       AliError(Form("Invalid id: %d", actGeo[i]));
516       continue;
517     }
518     TGeoVolume *sens = gGeoManager->GetVolume(actGeo[i]);
519     if (!sens) {
520       AliError(Form("No TGeo volume for sensitive volume ID=%d",actGeo[i]));
521       continue;
522     }   
523     actVmc[i] = mc->VolId(sens->GetName());
524     AliFMDDebug(1, ("Active vol id # %d: %d changed to %d", 
525                     i, actGeo[i], actVmc[i]));
526   }
527   fmd->SetActive(actVmc.fArray, actVmc.fN);
528   // fmd->InitTransformations();
529 }
530
531 //____________________________________________________________________
532 void
533 AliFMD::FinishEvent()
534 {
535   // Called at the end of the an event in simulations.  If the debug
536   // level is high enough, then the `bad' hits are printed.
537   // 
538   if (AliLog::GetDebugLevel("FMD", "AliFMD") < 10) return;
539   if (fBad && fBad->GetEntries() > 0) {
540     AliWarning(Form("got %d 'bad' hits", fBad->GetEntries()));
541     TIter next(fBad);
542     AliFMDHit* hit;
543     while ((hit = static_cast<AliFMDHit*>(next()))) hit->Print("D");
544     fBad->Clear();
545   }
546 }
547
548
549
550 //====================================================================
551 //
552 // Hit and Digit managment 
553 //
554 //____________________________________________________________________
555 void 
556 AliFMD::MakeBranch(Option_t * option)
557 {
558   // Create Tree branches for the FMD.
559   //
560   // Options:
561   //
562   //    H          Make a branch of TClonesArray of AliFMDHit's
563   //    D          Make a branch of TClonesArray of AliFMDDigit's
564   //    S          Make a branch of TClonesArray of AliFMDSDigit's
565   // 
566   const Int_t kBufferSize = 16000;
567   TString branchname(GetName());
568   TString opt(option);
569   
570   if (opt.Contains("H", TString::kIgnoreCase)) {
571     HitsArray();
572     AliDetector::MakeBranch(option); 
573   }
574   if (opt.Contains("D", TString::kIgnoreCase)) { 
575     DigitsArray();
576     MakeBranchInTree(fLoader->TreeD(), branchname.Data(),
577                      &fDigits, kBufferSize, 0);
578   }
579   if (opt.Contains("S", TString::kIgnoreCase)) { 
580     SDigitsArray();
581     MakeBranchInTree(fLoader->TreeS(), branchname.Data(),
582                      &fSDigits, kBufferSize, 0);
583   }
584 }
585
586 //____________________________________________________________________
587 void 
588 AliFMD::SetTreeAddress()
589 {
590   // Set branch address for the Hits, Digits, and SDigits Tree.
591   if (fLoader->TreeH()) HitsArray();
592   AliDetector::SetTreeAddress();
593
594   TTree *treeD = fLoader->TreeD();
595   if (treeD) {
596     DigitsArray();
597     TBranch* branch = treeD->GetBranch ("FMD");
598     if (branch) branch->SetAddress(&fDigits);
599   }
600
601   TTree *treeS = fLoader->TreeS();
602   if (treeS) {
603     SDigitsArray();
604     TBranch* branch = treeS->GetBranch ("FMD");
605     if (branch) branch->SetAddress(&fSDigits);
606   }
607 }
608
609 //____________________________________________________________________
610 void 
611 AliFMD::SetHitsAddressBranch(TBranch *b)
612 {
613   // Set the TClonesArray to read hits into. 
614   b->SetAddress(&fHits);
615 }
616 //____________________________________________________________________
617 void 
618 AliFMD::SetSDigitsAddressBranch(TBranch *b)
619 {
620   // Set the TClonesArray to read hits into. 
621   b->SetAddress(&fSDigits);
622 }
623
624 //____________________________________________________________________
625 void 
626 AliFMD::AddHit(Int_t track, Int_t *vol, Float_t *hits) 
627 {
628   // Add a hit to the hits tree 
629   // 
630   // The information of the two arrays are decoded as 
631   // 
632   // Parameters
633   //    track                Track #
634   //    ivol[0]  [UShort_t ] Detector # 
635   //    ivol[1]  [Char_t   ] Ring ID 
636   //    ivol[2]  [UShort_t ] Sector #
637   //    ivol[3]  [UShort_t ] Strip # 
638   //    hits[0]  [Float_t  ] Track's X-coordinate at hit 
639   //    hits[1]  [Float_t  ] Track's Y-coordinate at hit
640   //    hits[3]  [Float_t  ] Track's Z-coordinate at hit
641   //    hits[4]  [Float_t  ] X-component of track's momentum             
642   //    hits[5]  [Float_t  ] Y-component of track's momentum             
643   //    hits[6]  [Float_t  ] Z-component of track's momentum            
644   //    hits[7]  [Float_t  ] Energy deposited by track                  
645   //    hits[8]  [Int_t    ] Track's particle Id # 
646   //    hits[9]  [Float_t  ] Time when the track hit
647   // 
648   // 
649   AddHitByFields(track, 
650                  UShort_t(vol[0]),  // Detector # 
651                  Char_t(vol[1]),    // Ring ID
652                  UShort_t(vol[2]),  // Sector # 
653                  UShort_t(vol[3]),  // Strip # 
654                  hits[0],           // X
655                  hits[1],           // Y
656                  hits[2],           // Z
657                  hits[3],           // Px
658                  hits[4],           // Py
659                  hits[5],           // Pz
660                  hits[6],           // Energy loss 
661                  Int_t(hits[7]),    // PDG 
662                  hits[8]);          // Time
663 }
664
665 //____________________________________________________________________
666 AliFMDHit*
667 AliFMD::AddHitByFields(Int_t    track, 
668                        UShort_t detector, 
669                        Char_t   ring, 
670                        UShort_t sector, 
671                        UShort_t strip, 
672                        Float_t  x, 
673                        Float_t  y, 
674                        Float_t  z,
675                        Float_t  px, 
676                        Float_t  py, 
677                        Float_t  pz,
678                        Float_t  edep,
679                        Int_t    pdg,
680                        Float_t  t, 
681                        Float_t  l, 
682                        Bool_t   stop)
683 {
684   // Add a hit to the list
685   //
686   // Parameters:
687   // 
688   //    track     Track #
689   //    detector  Detector # (1, 2, or 3)                      
690   //    ring      Ring ID ('I' or 'O')
691   //    sector    Sector # (For inner/outer rings: 0-19/0-39)
692   //    strip     Strip # (For inner/outer rings: 0-511/0-255)
693   //    x         Track's X-coordinate at hit
694   //    y         Track's Y-coordinate at hit
695   //    z         Track's Z-coordinate at hit
696   //    px        X-component of track's momentum 
697   //    py        Y-component of track's momentum
698   //    pz        Z-component of track's momentum
699   //    edep      Energy deposited by track
700   //    pdg       Track's particle Id #
701   //    t         Time when the track hit 
702   //    l         Track length through the material. 
703   //    stop      Whether track was stopped or disappeared
704   // 
705   TClonesArray& a = *(HitsArray());
706   // Search through the list of already registered hits, and see if we
707   // find a hit with the same parameters.  If we do, then don't create
708   // a new hit, but rather update the energy deposited in the hit.
709   // This is done, so that a FLUKA based simulation will get the
710   // number of hits right, not just the enerrgy deposition. 
711   AliFMDHit* hit = 0;
712   for (Int_t i = 0; i < fNhits; i++) {
713     if (!a.At(i)) continue;
714     hit = static_cast<AliFMDHit*>(a.At(i));
715     if (hit->Detector() == detector 
716         && hit->Ring() == ring
717         && hit->Sector() == sector 
718         && hit->Strip() == strip
719         && hit->Track() == track) {
720       AliFMDDebug(1, ("already had a hit in FMD%d%c[%2d,%3d] for track # %d,"
721                        " adding energy (%f) to that hit (%f) -> %f", 
722                        detector, ring, sector, strip, track, edep, hit->Edep(),
723                        hit->Edep() + edep));
724       hit->SetEdep(hit->Edep() + edep);
725       return hit;
726     }
727   }
728   // If hit wasn't already registered, do so know. 
729   hit = new (a[fNhits]) AliFMDHit(fIshunt, track, detector, ring, sector, 
730                                   strip, x, y, z, px, py, pz, edep, pdg, t, 
731                                   l, stop);
732   // gMC->AddTrackReference(track, 12);
733   fNhits++;
734   
735   //Reference track
736
737   AliMC *mcApplication = (AliMC*)gAlice->GetMCApp();
738   
739   AliTrackReference* trackRef = 
740     AddTrackReference(mcApplication->GetCurrentTrackNumber(), 
741                       AliTrackReference::kFMD); 
742   UInt_t stripId = AliFMDStripIndex::Pack(detector,ring,sector,strip);
743   trackRef->SetUserId(stripId);
744   
745   
746   
747   return hit;
748 }
749
750 //____________________________________________________________________
751 void 
752 AliFMD::AddDigit(Int_t* digits, Int_t*)
753 {
754   // Add a digit to the Digit tree 
755   // 
756   // Paramters 
757   //
758   //    digits[0]  [UShort_t] Detector #
759   //    digits[1]  [Char_t]   Ring ID
760   //    digits[2]  [UShort_t] Sector #
761   //    digits[3]  [UShort_t] Strip #
762   //    digits[4]  [UShort_t] ADC Count 
763   //    digits[5]  [Short_t]  ADC Count, -1 if not used
764   //    digits[6]  [Short_t]  ADC Count, -1 if not used 
765   // 
766   AddDigitByFields(UShort_t(digits[0]),  // Detector #
767                    Char_t(digits[1]),    // Ring ID
768                    UShort_t(digits[2]),  // Sector #
769                    UShort_t(digits[3]),  // Strip #
770                    UShort_t(digits[4]),  // ADC Count1 
771                    Short_t(digits[5]),   // ADC Count2 
772                    Short_t(digits[6]),   // ADC Count3 
773                    Short_t(digits[7])); 
774 }
775
776 //____________________________________________________________________
777 void 
778 AliFMD::AddDigitByFields(UShort_t       detector, 
779                          Char_t         ring, 
780                          UShort_t       sector, 
781                          UShort_t       strip, 
782                          UShort_t       count1, 
783                          Short_t        count2,
784                          Short_t        count3, 
785                          Short_t        count4,
786                          UShort_t       nrefs,
787                          Int_t*         refs)
788 {
789   // add a real digit - as coming from data
790   // 
791   // Parameters 
792   //
793   //    detector  Detector # (1, 2, or 3)                      
794   //    ring      Ring ID ('I' or 'O')
795   //    sector    Sector # (For inner/outer rings: 0-19/0-39)
796   //    strip     Strip # (For inner/outer rings: 0-511/0-255)
797   //    count1    ADC count (a 10-bit word)
798   //    count2    ADC count (a 10-bit word), or -1 if not used
799   //    count3    ADC count (a 10-bit word), or -1 if not used
800   TClonesArray& a = *(DigitsArray());
801   
802   AliFMDDebug(15, ("Adding digit # %5d/%5d for FMD%d%c[%2d,%3d]"
803                    "=(%d,%d,%d,%d) with %d tracks",
804                    fNdigits-1, a.GetEntriesFast(),
805                    detector, ring, sector, strip, 
806                    count1, count2, count3, count4, nrefs));
807   new (a[fNdigits++]) 
808     AliFMDDigit(detector, ring, sector, strip, 
809                 count1, count2, count3, count4, nrefs, refs);
810   
811 }
812
813 //____________________________________________________________________
814 void 
815 AliFMD::AddSDigit(Int_t* digits)
816 {
817   // Add a digit to the SDigit tree 
818   // 
819   // Paramters 
820   //
821   //    digits[0]  [UShort_t] Detector #
822   //    digits[1]  [Char_t]   Ring ID
823   //    digits[2]  [UShort_t] Sector #
824   //    digits[3]  [UShort_t] Strip #
825   //    digits[4]  [Float_t]  Total energy deposited 
826   //    digits[5]  [UShort_t] ADC Count 
827   //    digits[6]  [Short_t]  ADC Count, -1 if not used
828   //    digits[7]  [Short_t]  ADC Count, -1 if not used 
829   // 
830   AddSDigitByFields(UShort_t(digits[0]),   // Detector #
831                     Char_t(digits[1]),     // Ring ID
832                     UShort_t(digits[2]),   // Sector #
833                     UShort_t(digits[3]),   // Strip #
834                     Float_t(digits[4]),    // Edep
835                     UShort_t(digits[5]),   // ADC Count1 
836                     Short_t(digits[6]),    // ADC Count2 
837                     Short_t(digits[7]),    // ADC Count3 
838                     Short_t(digits[8]),    // ADC Count4
839                     UShort_t(digits[9]),   // N particles
840                     UShort_t(digits[10])); // N primaries
841 }
842
843 //____________________________________________________________________
844 void 
845 AliFMD::AddSDigitByFields(UShort_t       detector, 
846                           Char_t         ring, 
847                           UShort_t       sector, 
848                           UShort_t       strip, 
849                           Float_t        edep,
850                           UShort_t       count1, 
851                           Short_t        count2,
852                           Short_t        count3, 
853                           Short_t        count4, 
854                           UShort_t       ntot, 
855                           UShort_t       nprim,
856                           Int_t*         refs)
857 {
858   // add a summable digit
859   // 
860   // Parameters 
861   //
862   //    detector  Detector # (1, 2, or 3)                      
863   //    ring      Ring ID ('I' or 'O')
864   //    sector    Sector # (For inner/outer rings: 0-19/0-39)
865   //    strip     Strip # (For inner/outer rings: 0-511/0-255)
866   //    edep      Total energy deposited
867   //    count1    ADC count (a 10-bit word)
868   //    count2    ADC count (a 10-bit word), or -1 if not used
869   //    count3    ADC count (a 10-bit word), or -1 if not used
870   //
871   TClonesArray& a = *(SDigitsArray());
872   // AliFMDDebug(0, ("Adding sdigit # %d", fNsdigits));
873   
874   AliFMDDebug(15, ("Adding sdigit # %5d/%5d for FMD%d%c[%2d,%3d]"
875                    "=(%d,%d,%d,%d) with %d tracks %d primaries (%p)",
876                    fNsdigits-1, a.GetEntriesFast(),
877                    detector, ring, sector, strip, 
878                    count1, count2, count3, count4, ntot, nprim, refs));
879   new (a[fNsdigits++]) 
880     AliFMDSDigit(detector, ring, sector, strip, edep, 
881                  count1, count2, count3, count4, ntot, nprim, refs);
882 }
883
884 //____________________________________________________________________
885 void 
886 AliFMD::ResetSDigits()
887 {
888   // Reset number of digits and the digits array for this detector. 
889   //
890   fNsdigits   = 0;
891   if (fSDigits) fSDigits->Clear();
892 }
893
894
895 //____________________________________________________________________
896 TClonesArray*
897 AliFMD::HitsArray() 
898 {
899   // Initialize hit array if not already, and return pointer to it. 
900   if (!fHits) { 
901     fHits = new TClonesArray("AliFMDHit", 1000);
902     fNhits = 0;
903     if (gAlice && gAlice->GetMCApp() && gAlice->GetMCApp()->GetHitLists()) 
904       gAlice->GetMCApp()->AddHitList(fHits);
905   }
906   return fHits;
907 }
908
909 //____________________________________________________________________
910 TClonesArray*
911 AliFMD::DigitsArray() 
912 {
913   // Initialize digit array if not already, and return pointer to it. 
914   if (!fDigits) { 
915     fDigits = new TClonesArray("AliFMDDigit", 1000);
916     fNdigits = 0;
917   }
918   return fDigits;
919 }
920
921 //____________________________________________________________________
922 TClonesArray*
923 AliFMD::SDigitsArray() 
924 {
925   // Initialize digit array if not already, and return pointer to it. 
926   if (!fSDigits) { 
927     fSDigits = new TClonesArray("AliFMDSDigit", 1000);
928     fNsdigits = 0;
929   }
930   return fSDigits;
931 }
932
933 //====================================================================
934 //
935 // Digitization 
936 //
937 //____________________________________________________________________
938 void 
939 AliFMD::Hits2Digits() 
940 {
941   // Create AliFMDDigit's from AliFMDHit's.  This is done by making a
942   // AliFMDDigitizer, and executing that code.
943   // 
944   AliFMDHitDigitizer digitizer(this, AliFMDHitDigitizer::kDigits);
945   digitizer.Init();
946   digitizer.Exec("");
947 }
948
949 //____________________________________________________________________
950 void 
951 AliFMD::Hits2SDigits() 
952 {
953   // Create AliFMDSDigit's from AliFMDHit's.  This is done by creating
954   // an AliFMDSDigitizer object, and executing it. 
955   // 
956   AliFMDHitDigitizer digitizer(this, AliFMDHitDigitizer::kSDigits);
957   digitizer.Init();
958   digitizer.Exec("");
959 }
960
961   
962 //____________________________________________________________________
963 AliDigitizer* 
964 AliFMD::CreateDigitizer(AliRunDigitizer* manager) const
965 {
966   // Create a digitizer object 
967   
968   /* This is what we probably _should_ do */
969   AliFMDBaseDigitizer* digitizer = 0;
970   
971 #ifdef USE_SSDIGITIZER
972   digitizer = new AliFMDSSDigitizer(manager);
973 #else 
974   /* This is what we actually do, and will work */
975 #if 0
976   AliInfo("SDigit->Digit conversion not really supported, "
977           "doing Hit->Digit conversion instead");
978 #endif
979   digitizer = new AliFMDDigitizer(manager);
980 #endif
981   return digitizer;
982 }
983
984 //====================================================================
985 //
986 // Raw data simulation 
987 //
988 //__________________________________________________________________
989 void 
990 AliFMD::Digits2Raw() 
991 {
992   // Turn digits into raw data. 
993   // 
994   // This uses the class AliFMDRawWriter to do the job.   Please refer
995   // to that class for more information. 
996   AliFMDRawWriter writer(this);
997   writer.Exec();
998 }
999
1000 //====================================================================
1001 //
1002 // Raw data reading 
1003 //
1004 //__________________________________________________________________
1005 Bool_t
1006 AliFMD::Raw2SDigits(AliRawReader* reader) 
1007 {
1008   // Turn digits into raw data. 
1009   // 
1010   // This uses the class AliFMDRawWriter to do the job.   Please refer
1011   // to that class for more information. 
1012   AliFMDParameters::Instance()->Init();
1013   MakeTree("S");
1014   MakeBranch("S");
1015   
1016   TClonesArray*       sdigits = SDigits();
1017   AliFMDReconstructor rec;
1018   
1019   // The two boolean arguments
1020   //   Make sdigits instead of digits 
1021   //   Subtract the pedestal off the signal
1022   rec.Digitize(reader, sdigits);
1023   // 
1024   // Bool_t ret = fmdReader.ReadAdcs(sdigits, kTRUE, kTRUE);
1025   // sdigits->ls();
1026   UShort_t ns = sdigits->GetEntriesFast();
1027   for (UShort_t i = 0; i < ns; i++) 
1028     sdigits->At(i)->Print("pl");
1029   
1030   AliFMDDebug(1, ("Got a total of %d SDigits", ns));
1031
1032   fLoader->TreeS()->Fill();
1033   ResetSDigits();
1034   fLoader->WriteSDigits("OVERWRITE");
1035
1036   return kTRUE;
1037 }
1038
1039
1040 //====================================================================
1041 //
1042 // Utility 
1043 //
1044 //__________________________________________________________________
1045 void 
1046 AliFMD::Browse(TBrowser* b) 
1047 {
1048   // Browse this object. 
1049   //
1050   AliFMDDebug(30, ("\tBrowsing the FMD"));
1051   AliDetector::Browse(b);
1052   b->Add(AliFMDGeometry::Instance());
1053 }
1054
1055 //____________________________________________________________________  
1056 void
1057 AliFMD::AddAlignableVolumes() const
1058 {
1059   //
1060   // Create entries for alignable volumes associating the symbolic volume
1061   // name with the corresponding volume path. Needs to be syncronized with
1062   // eventual changes in the geometry.
1063   // 
1064   // This code was made by Raffaele Grosso <rgrosso@mail.cern.ch>.  I
1065   // (cholm) will probably want to change it.   For one, I think it
1066   // should be the job of the geometry manager to deal with this. 
1067   AliInfo("Add FMD alignable volumes");
1068   AliFMDGeometry::Instance()->SetAlignableVolumes();
1069 #if 0  
1070   for(size_t f = 1; f <= 3; f++){ // Detector 1,2,3
1071     for(size_t tb =  0; tb <2 ; tb++){ // Top/Bottom 
1072       char     stb = tb == 0 ? 'T' : 'B';
1073       unsigned min = tb == 0 ? 0   : 5;
1074
1075       TString halfVol(Form("/ALIC_1/F%dM%c_%d", f, stb, f));
1076       TString halfSym(halfVol);
1077       if(!gGeoManager->SetAlignableEntry(halfSym.Data(),halfVol.Data()))
1078         AliFatal(Form("Alignable entry %s not created. "
1079                       "Volume path %s not valid", 
1080                       halfSym.Data(),halfVol.Data()));
1081       for(size_t io = 0; io < 2; io++){ // inner, outer 
1082         if (f==1 && io==1) continue; // Only one ring in FMD1 
1083         if(tb == 1 && io==1) min=10;
1084         char     sio = (io == 0 ? 'I' : 'O');
1085         unsigned nio = (io == 0 ? 3   : 9);
1086         unsigned max = (io == 0 ? 5   : 10) + min;
1087         
1088         for(size_t i = min; i < max; i++) { // Modules
1089           TString modVol(Form("%s/F%c%cV_7%d/F%cSE_%d", halfVol.Data(), 
1090                               sio, stb, nio, sio, i));
1091           TString modSym(modVol);
1092           if(!gGeoManager->SetAlignableEntry(modSym.Data(),modVol.Data()))
1093             AliFatal(Form("Alignable entry %s not created. "
1094                           "Volume path %s not valid", 
1095                           modSym.Data(), modVol.Data()));
1096         }
1097       }
1098     }
1099   }
1100 #endif
1101 }
1102 //___________________________________________________________________
1103 //
1104 // EOF
1105 //