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