Changes to compile with Root6 on macosx64
[u/mrichter/AliRoot.git] / MFT / AliMFT.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
16 //====================================================================================================================================================
17 //
18 //      Main class of the ALICE Muon Forward Tracker
19 //
20 //      Contact author: antonio.uras@cern.ch
21 //
22 //====================================================================================================================================================
23
24 #include "AliLog.h"
25 #include "TFile.h"  
26 #include "TGeoManager.h"    
27 #include "TGeoVolume.h"
28 #include "TGeoMatrix.h"
29 #include "TVirtualMC.h"
30 #include "TClonesArray.h"
31 #include "TGeoGlobalMagField.h"
32 #include "AliRun.h"
33 #include "AliLoader.h"
34 #include "AliDetector.h"
35 #include "AliMC.h"
36 #include "AliMagF.h"
37 #include "AliMFT.h"
38 #include "AliMFTHit.h"
39 #include "AliMFTDigit.h"
40 #include "AliMFTCluster.h"
41 #include "AliTrackReference.h"
42 #include "AliMFTSegmentation.h"
43 #include "AliMFTDigitizer.h"
44 #include "AliMFTPlane.h"
45 #include "TString.h"
46 #include "TObjArray.h"
47
48 ClassImp(AliMFT) 
49
50 //====================================================================================================================================================
51
52 AliMFT::AliMFT():
53   AliDetector(),
54   fVersion(1),
55   fNPlanes(0),
56   fNSlices(1),
57   fSDigitsPerPlane(0),
58   fDigitsPerPlane(0),
59   fRecPointsPerPlane(0),
60   fSideDigits(0),
61   fSegmentation(0),
62   fNameGeomFile(0),
63   fChargeDispersion(25.e-4),
64   fSingleStepForChargeDispersion(0),
65   fNStepForChargeDispersion(4),
66   fDensitySupportOverSi(0.15)
67 {
68
69   // default constructor
70
71 }
72
73 //====================================================================================================================================================
74
75 AliMFT::AliMFT(const Char_t *name, const Char_t *title):
76   AliDetector(name, title),
77   fVersion(1),
78   fNPlanes(0),
79   fNSlices(1),
80   fSDigitsPerPlane(0),
81   fDigitsPerPlane(0),
82   fRecPointsPerPlane(0),
83   fSideDigits(0),
84   fSegmentation(0),
85   fNameGeomFile(0),
86   fChargeDispersion(25.e-4),
87   fSingleStepForChargeDispersion(0),
88   fNStepForChargeDispersion(4),
89   fDensitySupportOverSi(0.15)
90 {
91
92   fNameGeomFile = "AliMFTGeometry.root";
93
94   SetGeometry();
95
96   Init();
97
98 }
99
100 //====================================================================================================================================================
101
102 AliMFT::AliMFT(const Char_t *name, const Char_t *title, Char_t *nameGeomFile):
103   AliDetector(name, title),
104   fVersion(1),
105   fNPlanes(0),
106   fNSlices(1),
107   fSDigitsPerPlane(0),
108   fDigitsPerPlane(0),
109   fRecPointsPerPlane(0),
110   fSideDigits(0),
111   fSegmentation(0),
112   fNameGeomFile(0),
113   fChargeDispersion(25.e-4),
114   fSingleStepForChargeDispersion(0),
115   fNStepForChargeDispersion(4),
116   fDensitySupportOverSi(0.15)
117 {
118
119   fNameGeomFile = nameGeomFile;
120
121   SetGeometry();
122
123   Init();
124
125 }
126
127 //====================================================================================================================================================
128
129 AliMFT::~AliMFT() {
130
131   if (fSDigitsPerPlane)   { fSDigitsPerPlane->Delete();    delete fSDigitsPerPlane;   }
132   if (fDigitsPerPlane)    { fDigitsPerPlane->Delete();     delete fDigitsPerPlane;    }
133   if (fRecPointsPerPlane) { fRecPointsPerPlane->Delete();  delete fRecPointsPerPlane; }
134
135 }
136
137 //====================================================================================================================================================
138
139 void AliMFT::CreateMaterials() {
140
141   // Definition of MFT materials  - to be updated to the most recent values
142
143   AliInfo("Start MFT materials");
144
145   // data from PDG booklet 2002                 density [gr/cm^3]     rad len [cm]           abs len [cm]    
146   Float_t   aSi = 28.085 ,    zSi   = 14. ,     dSi      =  2.329 ,   radSi   =  21.82/dSi , absSi   = 108.4/dSi  ;    // Silicon
147   Float_t   aCarb = 12.01 ,   zCarb =  6. ,     dCarb    =  2.265 ,   radCarb =  18.8 ,      absCarb = 49.9       ;    // Carbon
148   Float_t   aAlu = 26.98 ,    zAlu  = 13. ,     dAlu     =  2.70  ,   radAlu  =  8.897 ,     absAlu  = 39.70      ;    // Aluminum
149
150   // Air mixture
151   const Int_t nAir = 4;   
152   Float_t   aAir[nAir] = {12, 14, 16, 36} ,  zAir[nAir] = {6, 7, 8, 18} ,   wAir[nAir]={0.000124, 0.755267, 0.231781, 0.012827} , dAir=0.00120479;   
153
154   // Water mixture
155   const Int_t nWater = 2;   
156   Float_t   aWater[nWater] = {1.00794, 15.9994} ,  zWater[nWater] = {1, 8} ,   wWater[nWater] = {0.111894, 0.888106} , dWater=1.;   
157
158   // SiO2 mixture
159   const Int_t nSiO2 = 2; 
160   Float_t   aSiO2[nSiO2] = {15.9994, 28.0855} ,   zSiO2[nSiO2] = {8., 14.} ,   wSiO2[nSiO2] = {0.532565, 0.467435} , dSiO2 = 2.20;  
161
162   // Inox mixture
163   const Int_t nInox = 9;
164   Float_t   aInox[nInox] = {12.0107, 54.9380, 28.0855, 30.9738, 32.0660, 58.6928, 51.9961, 95.9400, 55.8450} ;   
165   Float_t   zInox[nInox] = { 6,      25,      14,      15,      16,      28,      24,      42,      26     } ;   
166   Float_t   wInox[nInox] = {0.0003,  0.02,    0.01,    0.00045, 0.0003,  0.12,    0.17,    0.025,   0.65395} ;
167   Float_t   dInox = 8.03; 
168   
169   Int_t   matId  = 0;                        // tmp material id number
170   Int_t   unsens = 0, sens=1;                // sensitive or unsensitive medium
171   Int_t   itgfld = 3;                        // type of field intergration 0 no field -1 user in guswim 1 Runge Kutta 2 helix 3 const field along z
172   Float_t maxfld = 5.;                       // max field value
173
174   Float_t tmaxfd = -10.0;                    // max deflection angle due to magnetic field in one step
175   Float_t stemax =  0.001;                   // max step allowed [cm]
176   Float_t deemax = -0.2;                     // maximum fractional energy loss in one step 0<deemax<=1  
177   Float_t epsil  =  0.001;                   // tracking precision [cm]   
178   Float_t stmin  = -0.001;                   // minimum step due to continuous processes [cm] (negative value: choose it automatically)
179
180   Float_t tmaxfdSi =  0.1;                   // max deflection angle due to magnetic field in one step
181   Float_t stemaxSi =  5.0e-4;                // maximum step allowed [cm]
182   Float_t deemaxSi =  0.1;                   // maximum fractional energy loss in one step 0<deemax<=1
183   Float_t epsilSi  =  0.5e-4;                // tracking precision [cm]
184   Float_t stminSi  = -0.001;                 // minimum step due to continuous processes [cm] (negative value: choose it automatically)
185
186   Int_t   isxfld = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Integ();   // from CreateMaterials in STRUCT/AliPIPEv3.cxx
187   Float_t sxmgmx = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Max();     // from CreateMaterials in STRUCT/AliPIPEv3.cxx
188       
189   AliMixture(++matId,"Air", aAir, zAir, dAir, nAir, wAir); 
190   AliMedium(kAir,    "Air", matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
191   
192   AliMaterial(++matId, "Si", aSi, zSi, dSi, radSi, absSi);  
193   AliMedium(kSi,       "Si", matId, sens, isxfld, sxmgmx, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
194
195   AliMaterial(++matId, "Readout", aSi, zSi, dSi, radSi, absSi);  
196   AliMedium(kReadout,  "Readout", matId, unsens, isxfld, sxmgmx, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
197
198   AliMaterial(++matId, "Support", aSi, zSi, dSi*fDensitySupportOverSi, radSi/fDensitySupportOverSi, absSi/fDensitySupportOverSi);  
199   AliMedium(kSupport,  "Support", matId, unsens, isxfld, sxmgmx, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
200   
201   AliMaterial(++matId, "Carbon", aCarb, zCarb, dCarb, radCarb, absCarb );
202   AliMedium(kCarbon,   "Carbon", matId, unsens, isxfld,  sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
203   
204   AliMaterial(++matId, "Alu", aAlu, zAlu, dAlu, radAlu, absAlu);
205   AliMedium(kAlu,      "Alu", matId, unsens, isxfld,  sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
206   
207   AliMixture(++matId, "Water", aWater, zWater, dWater, nWater, wWater);
208   AliMedium(kWater,   "Water", matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
209   
210   AliMixture(++matId, "SiO2", aSiO2, zSiO2, dSiO2, nSiO2, wSiO2);
211   AliMedium(kSiO2,    "SiO2", matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
212   
213   AliMixture(++matId, "Inox", aInox, zInox, dInox, nInox, wInox);
214   AliMedium(kInox,    "Inox", matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
215
216   AliInfo("End MFT materials");
217           
218 }
219
220 //====================================================================================================================================================
221
222 void AliMFT::CreateGeometry() {
223
224   // Creates detailed geometry simulation (currently GEANT volumes tree)        
225
226   AliInfo("Start MFT preliminary version building");
227   if(!TVirtualMC::GetMC()->IsRootGeometrySupported()) return;                
228   TGeoVolumeAssembly *vol = CreateVol();
229   AliInfo("TGeoVolumeAssembly created!");
230   gGeoManager->GetVolume("ALIC")->AddNode(vol,0);
231   AliInfo("Stop MFT preliminary version building");
232
233   if (fNStepForChargeDispersion) fSingleStepForChargeDispersion = fChargeDispersion/Double_t(fNStepForChargeDispersion);
234
235
236
237 //====================================================================================================================================================
238
239 void AliMFT::AddAlignableVolumes() {
240
241   // Create entries for alignable volumes associating the symbolic volume
242   // name with the corresponding volume path. Needs to be syncronized with
243   // eventual changes in the geometry.
244
245   TString sysName = "MFT";
246   TString volPath = "/ALIC_1/MFT_0";
247   
248   if (!gGeoManager->SetAlignableEntry(sysName.Data(),volPath.Data())) {
249     AliFatal(Form("Alignable entry %s not created. Volume path %s not valid", sysName.Data(), volPath.Data()));
250   }  
251
252 }
253
254 //====================================================================================================================================================
255
256 void AliMFT::StepManager() {
257
258   // Full Step Manager
259
260   AliDebug(2, Form("Entering StepManager: TVirtualMC::GetMC()->CurrentVolName() = %s", TVirtualMC::GetMC()->CurrentVolName()));
261
262   if (!fSegmentation) AliFatal("No segmentation available");    // DO WE HAVE A SEGMENTATION???
263
264   if (!(this->IsActive())) return;
265   if (!(TVirtualMC::GetMC()->TrackCharge())) return;
266
267   TString planeNumber   = TVirtualMC::GetMC()->CurrentVolName();
268   TString detElemNumber = TVirtualMC::GetMC()->CurrentVolName();
269   if (planeNumber.Contains("support")) return;
270   if (planeNumber.Contains("readout")) return;
271   planeNumber.Remove(0,9);
272   detElemNumber.Remove(0,19);
273   planeNumber.Remove(2);
274   detElemNumber.Remove(3);
275   Int_t detElemID = fSegmentation->GetDetElemGlobalID(planeNumber.Atoi(), detElemNumber.Atoi());
276
277   if (TVirtualMC::GetMC()->IsTrackExiting()) {
278     AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kMFT);
279   }
280
281   static TLorentzVector position, momentum;
282   static AliMFTHit hit;
283   
284   Int_t  status = 0;
285   
286   // Track status
287   if (TVirtualMC::GetMC()->IsTrackInside())      status +=  1;
288   if (TVirtualMC::GetMC()->IsTrackEntering())    status +=  2;
289   if (TVirtualMC::GetMC()->IsTrackExiting())     status +=  4;
290   if (TVirtualMC::GetMC()->IsTrackOut())         status +=  8;
291   if (TVirtualMC::GetMC()->IsTrackDisappeared()) status += 16;
292   if (TVirtualMC::GetMC()->IsTrackStop())        status += 32;
293   if (TVirtualMC::GetMC()->IsTrackAlive())       status += 64;
294
295   // ---------- Fill hit structure
296
297   hit.SetDetElemID(detElemID);
298   hit.SetPlane(planeNumber.Atoi());
299   hit.SetTrack(gAlice->GetMCApp()->GetCurrentTrackNumber());
300     
301   TVirtualMC::GetMC()->TrackPosition(position);
302   TVirtualMC::GetMC()->TrackMomentum(momentum);
303
304   AliDebug(1, Form("AliMFT::StepManager()->%s Hit #%06d (x=%f, y=%f, z=%f) belongs to track %02d\n", 
305                    TVirtualMC::GetMC()->CurrentVolName(), fNhits, position.X(), position.Y(), position.Z(), gAlice->GetMCApp()->GetCurrentTrackNumber())); 
306
307   hit.SetPosition(position);
308   hit.SetTOF(TVirtualMC::GetMC()->TrackTime());
309   hit.SetMomentum(momentum);
310   hit.SetStatus(status);
311   hit.SetEloss(TVirtualMC::GetMC()->Edep());
312   //  hit.SetShunt(GetIshunt());
313 //   if (TVirtualMC::GetMC()->IsTrackEntering()) {
314 //     hit.SetStartPosition(position);
315 //     hit.SetStartTime(TVirtualMC::GetMC()->TrackTime());
316 //     hit.SetStartStatus(status);
317 //     return; // don't save entering hit.
318 //   } 
319
320   // Fill hit structure with this new hit.
321   new ((*fHits)[fNhits++]) AliMFTHit(hit);
322
323   // Save old position... for next hit.
324 //   hit.SetStartPosition(position);
325 //   hit.SetStartTime(TVirtualMC::GetMC()->TrackTime());
326 //   hit.SetStartStatus(status);
327
328   return;
329
330 }
331
332 //====================================================================================================================================================
333
334 TGeoVolumeAssembly* AliMFT::CreateVol() {
335
336   // method to create the MFT Geometry (silicon circular planes)
337
338   if (!fSegmentation) CreateGeometry();
339   
340   TGeoVolumeAssembly *vol = new TGeoVolumeAssembly("MFT");
341   TGeoMedium *silicon = gGeoManager->GetMedium("MFT_Si");
342   TGeoMedium *readout = gGeoManager->GetMedium("MFT_Readout");
343   TGeoMedium *support = gGeoManager->GetMedium("MFT_Support");
344   TGeoMedium *carbon = gGeoManager->GetMedium("MFT_Carbon");
345   TGeoMedium *alu = gGeoManager->GetMedium("MFT_Alu");
346   TGeoMedium *water = gGeoManager->GetMedium("MFT_Water");
347   TGeoMedium *si02 = gGeoManager->GetMedium("MFT_SiO2");
348   TGeoMedium *inox = gGeoManager->GetMedium("MFT_Inox");
349
350   // ---- Cage & Services Description --------------------------------------------
351   // R. Tieulent - 17/01/2014 - Basic description for ITS/TPC matching studies
352   // -----------------------------------------------------------------------------
353
354   TGeoVolumeAssembly *cageNservices = new TGeoVolumeAssembly("MFT_cageNservices");
355   
356   // cage definition
357   Float_t cageDz   = 150./2.;
358   Float_t cageRMax = 49.5;
359   Float_t cageRMin = cageRMax - 0.120; // 0.64% of X0
360   
361   TGeoVolume *cage = gGeoManager->MakeTube("MFT_cage", carbon, cageRMin, cageRMax, cageDz);
362   cage->SetLineColor(kBlue);
363
364   cageNservices->AddNode(cage,1,new TGeoTranslation(0., 0., 0. ));
365
366   // Services definition
367   TGeoVolumeAssembly *services = new TGeoVolumeAssembly("MFT_services");
368
369   // Aluminum bus-Bar
370   Float_t busBarDz = 150.;
371   Float_t busBarThick = 0.1 ;
372   Float_t busBarWidth = 1.;
373   
374   TGeoVolume *aluBusBar = gGeoManager->MakeBox("MFT_busbar", alu, busBarWidth/2., busBarThick/2., busBarDz/2.);
375   aluBusBar->SetLineColor(kYellow);
376
377   Int_t nBusBar = 30;
378   Float_t dPhiBusBar = 2.*TMath::Pi() / nBusBar;
379   Float_t dShift = cageRMin - busBarThick;
380   
381   TGeoRotation *rot;
382
383   for (Int_t iBusBar=0; iBusBar<nBusBar; iBusBar++) {
384     Float_t phi =  dPhiBusBar*iBusBar;
385     Float_t xp = dShift*TMath::Cos(phi);
386     Float_t yp = dShift*TMath::Sin(phi);
387     rot = new TGeoRotation();
388     rot->RotateZ(phi*TMath::RadToDeg()+90.);
389     services->AddNode(aluBusBar, iBusBar+1, new TGeoCombiTrans(xp,yp,0,rot));
390   }
391   
392   // Cooling Services  definition
393   TGeoVolumeAssembly *cooling = new TGeoVolumeAssembly("MFT_cooling");
394   // Cooling Water
395   Float_t coolingDz = 150.;
396   Float_t coolingR  = 0.3 /2. ; // 3mm in diameter
397   Float_t coolingDR = 0.02 ;    // Thickness of the pipe 0.2mm
398   
399   TGeoVolume *coolingWater = gGeoManager->MakeTube("MFT_coolingWater", water, 0., coolingR, coolingDz/2.);
400   coolingWater->SetLineColor(kCyan);
401   cooling->AddNode(coolingWater, 1, new TGeoTranslation(0,0,0 ));
402
403   // Cooling Pipes
404   TGeoVolume *coolingPipes = gGeoManager->MakeTube("MFT_coolingPipes", inox, coolingR, coolingR+coolingDR, coolingDz/2.);
405   coolingPipes->SetLineColor(kGray);
406   cooling->AddNode(coolingPipes,1,new TGeoTranslation(0,0,0 ));
407   
408   Int_t nCooling = 18;
409   dShift = cageRMin - coolingR ;
410   Float_t phi0 = 0.02;
411   
412   for (Int_t iCooling=0; iCooling<nCooling; iCooling++) {
413     Float_t phi ;
414     if (iCooling<nCooling/2) phi = dPhiBusBar*(iCooling+3) + phi0;
415     else phi = dPhiBusBar*(iCooling+9) + phi0;
416     Float_t xp = dShift*TMath::Cos(phi);
417     Float_t yp = dShift*TMath::Sin(phi);
418     services->AddNode(cooling, iCooling+1, new TGeoTranslation(xp,yp,0 ));
419   }
420   
421   // Optical Fibers
422   Float_t fiberDz = 150.;
423   Float_t fiberRadius = 0.0125 /2. ; // 0.125mm in diameter
424   
425   TGeoVolume *fiber = gGeoManager->MakeTube("MFT_fiber", si02, 0., fiberRadius, fiberDz/2.);
426   fiber->SetLineColor(kCyan);
427
428   Int_t nFiber = 340;
429   dShift = cageRMin - 2*fiberRadius;
430   phi0 = 0.03;
431   
432   for (Int_t iFiber=0; iFiber<nFiber; iFiber++) {
433     Float_t phi = dPhiBusBar*(Int_t)(iFiber/11) - phi0-(iFiber%11)*2.*TMath::ATan(fiberRadius/dShift);
434     Float_t xp  = dShift*TMath::Cos(phi);
435     Float_t yp  = dShift*TMath::Sin(phi);
436     services->AddNode(fiber, iFiber+1, new TGeoTranslation(xp,yp,0 ));
437   }
438
439   cageNservices->AddNode(services, 1, new TGeoTranslation(0., 0., 0. ));
440   
441   vol->AddNode(cageNservices,1,new TGeoTranslation(0., 0., 0. ));
442   
443   // ------------------- Creating volumes for MFT planes --------------------------------
444
445   Double_t origin[3] = {0};
446
447   for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
448
449     AliDebug(1, Form("Creating volumes for MFT plane %02d",iPlane));
450
451     AliMFTPlane *plane = fSegmentation->GetPlane(iPlane);
452
453     // --------- support element(s)
454     
455     origin[0] =  0.5*(plane->GetSupportElement(0)->GetAxis(0)->GetXmax() + plane->GetSupportElement(0)->GetAxis(0)->GetXmin());
456     origin[1] =  0.5*(plane->GetSupportElement(0)->GetAxis(1)->GetXmax() + plane->GetSupportElement(0)->GetAxis(1)->GetXmin());
457     origin[2] = -0.5*(plane->GetSupportElement(0)->GetAxis(2)->GetXmax() + plane->GetSupportElement(0)->GetAxis(2)->GetXmin());
458     TGeoVolume *supportElem = gGeoManager->MakeTube(Form("MFT_plane%02d_support", iPlane), support, 
459                                                     plane->GetRMinSupport(), 
460                                                     plane->GetRMaxSupport(),
461                                                     0.5*(plane->GetSupportElement(0)->GetAxis(2)->GetXmax() - 
462                                                          plane->GetSupportElement(0)->GetAxis(2)->GetXmin()) );
463     AliDebug(2, Form("Created vol %s", supportElem->GetName()));
464     supportElem->SetLineColor(kCyan-9);
465     vol -> AddNode(supportElem, 0, new TGeoTranslation(origin[0], origin[1], origin[2]));
466
467     AliDebug(1, "support elements created!");
468     
469     // --------- active elements
470
471     for (Int_t iActive=0; iActive<plane->GetNActiveElements(); iActive++) {
472       
473       Double_t dx = 0.5*TMath::Abs(plane->GetActiveElement(iActive)->GetAxis(0)->GetXmax() - plane->GetActiveElement(iActive)->GetAxis(0)->GetXmin());
474       Double_t dy = 0.5*TMath::Abs(plane->GetActiveElement(iActive)->GetAxis(1)->GetXmax() - plane->GetActiveElement(iActive)->GetAxis(1)->GetXmin());
475       Double_t dz = 0.5*TMath::Abs(plane->GetActiveElement(iActive)->GetAxis(2)->GetXmax() - plane->GetActiveElement(iActive)->GetAxis(2)->GetXmin());
476       dz /= Double_t(fNSlices);
477
478       origin[0] =  0.5*(plane->GetActiveElement(iActive)->GetAxis(0)->GetXmax() + plane->GetActiveElement(iActive)->GetAxis(0)->GetXmin());
479       origin[1] =  0.5*(plane->GetActiveElement(iActive)->GetAxis(1)->GetXmax() + plane->GetActiveElement(iActive)->GetAxis(1)->GetXmin());
480
481       for (Int_t iSlice=0; iSlice<fNSlices; iSlice++) {
482         origin[2] = -0.5*(plane->GetActiveElement(iActive)->GetAxis(2)->GetXmin() + 2*dz*(iSlice+1) + plane->GetActiveElement(iActive)->GetAxis(2)->GetXmin() + 2*dz*(iSlice) );
483         TGeoVolume *activeElem = gGeoManager->MakeBox(Form("MFT_plane%02d_active%03d_slice%02d", iPlane, iActive, iSlice), silicon, dx, dy, dz);
484         AliDebug(2, Form("Created vol %s", activeElem->GetName()));
485         activeElem->SetLineColor(kGreen);
486         vol -> AddNode(activeElem, 0, new TGeoTranslation(origin[0], origin[1], origin[2]));
487       }
488
489     }
490
491     AliDebug(1, "active elements created!");
492
493     // --------- readout elements
494
495     for (Int_t iReadout=0; iReadout<plane->GetNReadoutElements(); iReadout++) {
496       
497       Double_t dx = 0.5*TMath::Abs(plane->GetReadoutElement(iReadout)->GetAxis(0)->GetXmax() - plane->GetReadoutElement(iReadout)->GetAxis(0)->GetXmin());
498       Double_t dy = 0.5*TMath::Abs(plane->GetReadoutElement(iReadout)->GetAxis(1)->GetXmax() - plane->GetReadoutElement(iReadout)->GetAxis(1)->GetXmin());
499       Double_t dz = 0.5*TMath::Abs(plane->GetReadoutElement(iReadout)->GetAxis(2)->GetXmax() - plane->GetReadoutElement(iReadout)->GetAxis(2)->GetXmin());
500
501       origin[0] =  0.5*(plane->GetReadoutElement(iReadout)->GetAxis(0)->GetXmax() + plane->GetReadoutElement(iReadout)->GetAxis(0)->GetXmin());
502       origin[1] =  0.5*(plane->GetReadoutElement(iReadout)->GetAxis(1)->GetXmax() + plane->GetReadoutElement(iReadout)->GetAxis(1)->GetXmin());
503       origin[2] = -0.5*(plane->GetReadoutElement(iReadout)->GetAxis(2)->GetXmax() + plane->GetReadoutElement(iReadout)->GetAxis(2)->GetXmin());
504
505       TGeoVolume *readoutElem = gGeoManager->MakeBox(Form("MFT_plane%02d_readout%03d", iPlane, iReadout), readout, dx, dy, dz);
506       AliDebug(2, Form("Created vol %s", readoutElem->GetName()));
507       readoutElem->SetLineColor(kRed);
508       vol -> AddNode(readoutElem, 0, new TGeoTranslation(origin[0], origin[1], origin[2]));
509       
510     }
511
512     AliDebug(1, "readout elements created!");
513
514   }
515
516   return vol;
517
518 }
519
520 //====================================================================================================================================================
521
522 void AliMFT::Hits2SDigits(){
523   
524   // Interface method invoked from AliSimulation to create a list of sdigits corresponding to list of hits. Every hit generates one sdigit.
525
526   AliDebug(1,"Start Hits2SDigits.");
527   
528   if (!fSegmentation) CreateGeometry();
529  
530   if (!fLoader->TreeH()) fLoader->LoadHits();
531
532   if (!fLoader->TreeS()) {
533
534     for (Int_t iEvt=0;iEvt<fLoader->GetRunLoader()->GetNumberOfEvents(); iEvt++) {
535
536       fLoader->GetRunLoader()->GetEvent(iEvt);
537       fLoader->MakeTree("S");
538       MakeBranch("S");
539       SetTreeAddress();
540
541       AliDebug(1, Form("Event %03d: fLoader->TreeH()->GetEntries() = %2d", iEvt, Int_t(fLoader->TreeH()->GetEntries())));
542
543       for (Int_t iTrack=0; iTrack<fLoader->TreeH()->GetEntries(); iTrack++) {
544         fLoader->TreeH()->GetEntry(iTrack);     
545         Hits2SDigitsLocal(Hits(), GetSDigitsList(), iTrack);    // convert these hits to a list of sdigits  
546       }
547       
548       fLoader->TreeS()->Fill();
549       fLoader->WriteSDigits("OVERWRITE");
550       ResetSDigits();
551
552     }
553   }
554
555   fLoader->UnloadHits();
556   fLoader->UnloadSDigits();
557
558   AliDebug(1,"Stop Hits2SDigits.");
559   
560 }
561
562 //====================================================================================================================================================
563
564 void AliMFT::Hits2SDigitsLocal(TClonesArray *hits, const TObjArray *pSDig, Int_t track) {
565
566   //  Add sdigits of these hits to the list
567   
568   AliDebug(1, "Entering Hits2SDigitsLocal");
569   
570   if (!fSegmentation) CreateGeometry();
571   
572   TClonesArray *pSDigList[fNMaxPlanes];
573   for (Int_t iPlane=0; iPlane<fNMaxPlanes; iPlane++) pSDigList[iPlane] = NULL; 
574   for (Int_t iPlane=0; iPlane<fNPlanes;    iPlane++) { 
575     pSDigList[iPlane] = (TClonesArray*) (*pSDig)[iPlane];
576     AliDebug(1,Form("Entries of pSDigList %3d; plane: %02d,",pSDigList[iPlane]->GetEntries(),iPlane));
577     if (!track && pSDigList[iPlane]->GetEntries()!=0) AliErrorClass("Some of sdigits lists is not empty");
578   }
579   
580   for (Int_t iHit=0; iHit<hits->GetEntries(); iHit++) {
581
582     AliMFTHit *hit = (AliMFTHit*) hits->At(iHit);
583
584     // Creating "main digit"
585
586     AliMFTDigit *mainSDigit = new AliMFTDigit();
587     mainSDigit->SetEloss(hit->GetEloss());
588     mainSDigit->SetDetElemID(hit->GetDetElemID());
589     mainSDigit->SetPlane(hit->GetPlane());
590     mainSDigit->AddMCLabel(hit->GetTrack()); 
591
592     Int_t xPixel = -1;
593     Int_t yPixel = -1;
594     if (fSegmentation->Hit2PixelID(hit->X(), hit->Y(), mainSDigit->GetDetElemID(), xPixel, yPixel)) {
595       mainSDigit->SetPixID(xPixel, yPixel, 0);
596       mainSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(mainSDigit->GetDetElemID()), 
597                           fSegmentation->GetPixelSizeY(mainSDigit->GetDetElemID()),
598                           fSegmentation->GetPixelSizeZ(mainSDigit->GetDetElemID()));  
599       mainSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(mainSDigit->GetDetElemID(), xPixel), 
600                            fSegmentation->GetPixelCenterY(mainSDigit->GetDetElemID(), yPixel),
601                            fSegmentation->GetPixelCenterZ(mainSDigit->GetDetElemID(), 0));
602       new ((*fSideDigits)[fSideDigits->GetEntries()]) AliMFTDigit(*mainSDigit);
603       AliDebug(1, Form("Created new sdigit (%f, %f, %f) from hit (%f, %f, %f)",
604                        mainSDigit->GetPixelCenterX(), mainSDigit->GetPixelCenterY(), mainSDigit->GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
605     }
606
607     // creating "side digits" to simulate the effect of charge dispersion
608
609     Double_t pi4 = TMath::Pi()/4.;
610     for (Int_t iStep=0; iStep<fNStepForChargeDispersion; iStep++) {
611       Double_t shift = (iStep+1) * fSingleStepForChargeDispersion;
612       for (Int_t iAngle=0; iAngle<8; iAngle++) {
613         Double_t shiftX = shift*TMath::Cos(iAngle*pi4);
614         Double_t shiftY = shift*TMath::Sin(iAngle*pi4);
615         if (fSegmentation->Hit2PixelID(hit->X()+shiftX, hit->Y()+shiftY, hit->GetDetElemID(), xPixel, yPixel)) {
616           Bool_t digitExists = kFALSE;
617           for (Int_t iSideDigit=0; iSideDigit<fSideDigits->GetEntries(); iSideDigit++) {
618             if (xPixel==((AliMFTDigit*) fSideDigits->At(iSideDigit))->GetPixelX() && 
619                 yPixel==((AliMFTDigit*) fSideDigits->At(iSideDigit))->GetPixelY()) {
620               digitExists = kTRUE;
621               break;
622             }
623           }
624           if (!digitExists) {
625             AliMFTDigit *sideSDigit = new AliMFTDigit();
626             sideSDigit->SetEloss(0.);
627             sideSDigit->SetDetElemID(hit->GetDetElemID());
628             sideSDigit->SetPlane(hit->GetPlane());
629             sideSDigit->AddMCLabel(hit->GetTrack());
630             sideSDigit->SetPixID(xPixel, yPixel, 0);
631             sideSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(sideSDigit->GetDetElemID()), 
632                                     fSegmentation->GetPixelSizeY(sideSDigit->GetDetElemID()),
633                                     fSegmentation->GetPixelSizeZ(sideSDigit->GetDetElemID()));  
634             sideSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(sideSDigit->GetDetElemID(), xPixel), 
635                                      fSegmentation->GetPixelCenterY(sideSDigit->GetDetElemID(), yPixel),
636                                      fSegmentation->GetPixelCenterZ(sideSDigit->GetDetElemID(), 0)); 
637             new ((*fSideDigits)[fSideDigits->GetEntries()]) AliMFTDigit(*sideSDigit);
638           }
639         }
640       }
641     }
642     
643     // ------------ In case we should simulate a rectangular pattern of pixel...
644     
645     if (fSegmentation->GetPlane(mainSDigit->GetPlane())->HasPixelRectangularPatternAlongY()) {
646       for (Int_t iSDigit=0; iSDigit<fSideDigits->GetEntries(); iSDigit++) {
647         AliMFTDigit *mySDig = (AliMFTDigit*) (fSideDigits->At(iSDigit));
648         if (mySDig->GetPixelX()%2 == mySDig->GetPixelY()%2) {   // both pair or both odd
649           xPixel = mySDig->GetPixelX();
650           yPixel = mySDig->GetPixelY()+1;
651           if (fSegmentation->DoesPixelExist(mySDig->GetDetElemID(), xPixel, yPixel)) {
652             AliMFTDigit *newSDigit = new AliMFTDigit();
653             newSDigit->SetEloss(0.);
654             newSDigit->SetDetElemID(mySDig->GetDetElemID());
655             newSDigit->SetPlane(mySDig->GetDetElemID());
656             newSDigit->SetPixID(xPixel, yPixel, 0);
657             newSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(newSDigit->GetDetElemID()), 
658                                    fSegmentation->GetPixelSizeY(newSDigit->GetDetElemID()),
659                                    fSegmentation->GetPixelSizeZ(newSDigit->GetDetElemID()));  
660             newSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(newSDigit->GetDetElemID(), xPixel), 
661                                     fSegmentation->GetPixelCenterY(newSDigit->GetDetElemID(), yPixel),
662                                     fSegmentation->GetPixelCenterZ(newSDigit->GetDetElemID(), 0)); 
663             new ((*fSideDigits)[fSideDigits->GetEntries()]) AliMFTDigit(*newSDigit);
664           }
665         }
666         else {   // pair-odd
667           xPixel = mySDig->GetPixelX();
668           yPixel = mySDig->GetPixelY()-1;
669           if (fSegmentation->DoesPixelExist(mySDig->GetDetElemID(), xPixel, yPixel)) {
670             AliMFTDigit *newSDigit = new AliMFTDigit();
671             newSDigit->SetEloss(0.);
672             newSDigit->SetDetElemID(mySDig->GetDetElemID());
673             newSDigit->SetPlane(mySDig->GetPlane());
674             newSDigit->SetPixID(xPixel, yPixel, 0);
675             newSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(newSDigit->GetDetElemID()), 
676                                    fSegmentation->GetPixelSizeY(newSDigit->GetDetElemID()),
677                                    fSegmentation->GetPixelSizeZ(newSDigit->GetDetElemID()));  
678             newSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(newSDigit->GetDetElemID(), xPixel), 
679                                     fSegmentation->GetPixelCenterY(newSDigit->GetDetElemID(), yPixel),
680                                     fSegmentation->GetPixelCenterZ(newSDigit->GetDetElemID(), 0)); 
681             new ((*fSideDigits)[fSideDigits->GetEntries()]) AliMFTDigit(*newSDigit);
682           }
683         }
684       }
685     }
686
687     // -------- checking which pixels switched on have their diode actually within the charge dispersion radius
688
689     for (Int_t iSDigit=0; iSDigit<fSideDigits->GetEntries(); iSDigit++) {
690       AliMFTDigit *mySDig = (AliMFTDigit*) (fSideDigits->At(iSDigit));
691       Double_t distance = TMath::Sqrt(TMath::Power(mySDig->GetPixelCenterX()-hit->X(),2) + TMath::Power(mySDig->GetPixelCenterY()-hit->Y(),2));
692       if (fSegmentation->GetPlane(mySDig->GetPlane())->HasPixelRectangularPatternAlongY()) {
693         if (mySDig->GetPixelX()%2 == mySDig->GetPixelY()%2) {  // both pair or both odd
694           if (distance<fChargeDispersion) {
695             AliDebug(1, Form("Created new sdigit (%f, %f, %f) from hit (%f, %f, %f)",
696                              mySDig->GetPixelCenterX(), mySDig->GetPixelCenterY(), mySDig->GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
697             new ((*pSDigList[mySDig->GetPlane()])[pSDigList[mySDig->GetPlane()]->GetEntries()]) AliMFTDigit(*mySDig);
698             xPixel = mySDig->GetPixelX();
699             yPixel = mySDig->GetPixelY()+1;
700             if (fSegmentation->DoesPixelExist(mySDig->GetDetElemID(), xPixel, yPixel)) {
701               AliMFTDigit *newSDigit = new AliMFTDigit();
702               newSDigit->SetEloss(0.);
703               newSDigit->SetDetElemID(mySDig->GetDetElemID());
704               newSDigit->SetPlane(mySDig->GetPlane());
705               newSDigit->SetPixID(xPixel, yPixel, 0);
706               newSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(newSDigit->GetDetElemID()), 
707                                      fSegmentation->GetPixelSizeY(newSDigit->GetDetElemID()),
708                                      fSegmentation->GetPixelSizeZ(newSDigit->GetDetElemID()));  
709               newSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(newSDigit->GetDetElemID(), xPixel), 
710                                       fSegmentation->GetPixelCenterY(newSDigit->GetDetElemID(), yPixel),
711                                       fSegmentation->GetPixelCenterZ(newSDigit->GetDetElemID(), 0));
712               AliDebug(1, Form("Created new sdigit (%f, %f, %f) from hit (%f, %f, %f)",
713                                newSDigit->GetPixelCenterX(), newSDigit->GetPixelCenterY(), newSDigit->GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
714               new ((*pSDigList[newSDigit->GetPlane()])[pSDigList[newSDigit->GetPlane()]->GetEntries()]) AliMFTDigit(*newSDigit);
715             }
716           }
717         }
718       }
719       else {
720         if (distance<fChargeDispersion) {
721           AliDebug(1, Form("Created new sdigit (%f, %f, %f) from hit (%f, %f, %f)",
722                            mySDig->GetPixelCenterX(), mySDig->GetPixelCenterY(), mySDig->GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
723           new ((*pSDigList[mySDig->GetPlane()])[pSDigList[mySDig->GetPlane()]->GetEntries()]) AliMFTDigit(*mySDig);
724         }
725       }
726     }
727
728     fSideDigits->Delete(); 
729
730   }
731
732   AliDebug(1,"Exiting Hits2SDigitsLocal");
733
734 }
735
736 //====================================================================================================================================================
737
738 void AliMFT::MakeBranch(Option_t *option) {
739
740   // Create Tree branches 
741   AliDebug(1, Form("Start with option= %s.",option));
742   
743   const Int_t kBufSize = 4000;
744   
745   const Char_t *cH = strstr(option,"H");
746   const Char_t *cD = strstr(option,"D");
747   const Char_t *cS = strstr(option,"S");
748
749   if (cH && fLoader->TreeH()) {
750     CreateHits();
751     MakeBranchInTree(fLoader->TreeH(), "MFT", &fHits, kBufSize, 0);   
752   }
753
754   if (cS && fLoader->TreeS()) {
755     CreateSDigits();
756     for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) MakeBranchInTree(fLoader->TreeS(), 
757                                                                     Form("Plane_%02d",iPlane), 
758                                                                     &((*fSDigitsPerPlane)[iPlane]), 
759                                                                     kBufSize, 0);
760   }
761   
762   if (cD && fLoader->TreeD()) {
763     CreateDigits();
764     for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) MakeBranchInTree(fLoader->TreeD(), 
765                                                                     Form("Plane_%02d",iPlane), 
766                                                                     &((*fDigitsPerPlane)[iPlane]),
767                                                                     kBufSize, 0);
768   }
769
770   AliDebug(1,"Stop.");
771
772 }
773
774 //====================================================================================================================================================
775
776 void AliMFT::SetTreeAddress() {
777
778   AliDebug(1, "AliMFT::SetTreeAddress()");
779
780   //Set branch address for the Hits and Digits Tree.
781   AliDebug(1, "Start.");
782
783   AliDebug(1, Form("AliMFT::SetTreeAddress Hits  fLoader->TreeH() = %p\n", fLoader->TreeH()));
784   if (fLoader->TreeH() && fLoader->TreeH()->GetBranch("MFT")) {
785     CreateHits();
786     fLoader->TreeH()->SetBranchAddress("MFT", &fHits);
787   }
788     
789   AliDebug(1, Form("AliMFT::SetTreeAddress SDigits  fLoader->TreeS() = %p\n", fLoader->TreeS()));
790   if (fLoader->TreeS() && fLoader->TreeS()->GetBranch("Plane_00")) {
791     CreateSDigits();
792     for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
793       fLoader->TreeS()->SetBranchAddress(Form("Plane_%02d",iPlane), &((*fSDigitsPerPlane)[iPlane]));
794     }
795   }
796     
797   AliDebug(1, Form("AliMFT::SetTreeAddress Digits  fLoader->TreeD() = %p\n", fLoader->TreeD()));
798   if (fLoader->TreeD() && fLoader->TreeD()->GetBranch("Plane_00")) {
799     CreateDigits(); 
800     for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
801       fLoader->TreeD()->SetBranchAddress(Form("Plane_%02d",iPlane), &((*fDigitsPerPlane)[iPlane]));
802     }
803   }
804
805   AliDebug(1, Form("AliMFT::SetTreeAddress RecPoints  fLoader->TreeR() = %p\n", fLoader->TreeR()));
806   if (fLoader->TreeR() && fLoader->TreeR()->GetBranch("Plane_00")) {
807     CreateRecPoints(); 
808     for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
809       fLoader->TreeR()->SetBranchAddress(Form("Plane_%02d",iPlane), &((*fRecPointsPerPlane)[iPlane]));
810     }
811   }
812
813   AliDebug(1,"Stop.");
814
815 }
816
817 //====================================================================================================================================================
818
819 void AliMFT::SetGeometry() {
820
821   AliInfo("AliMFT::SetGeometry\n");
822
823   fSegmentation = new AliMFTSegmentation(fNameGeomFile.Data());
824
825   fNPlanes = fSegmentation->GetNPlanes();
826
827 }
828
829 //====================================================================================================================================================
830
831 void AliMFT::CreateHits() { 
832
833   // create array of hits
834
835   AliDebug(1, "AliMFT::CreateHits()");
836
837   if (fHits) return;    
838   fHits = new TClonesArray("AliMFTHit");  
839
840 }
841
842 //====================================================================================================================================================
843
844 void AliMFT::CreateSDigits() { 
845  
846   // create sdigits list
847
848   AliDebug(1, "AliMFT::CreateSDigits()");
849  
850   if (fSDigitsPerPlane) return; 
851   fSDigitsPerPlane = new TObjArray(fNPlanes); 
852   for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) fSDigitsPerPlane->AddAt(new TClonesArray("AliMFTDigit"), iPlane);
853
854   fSideDigits = new TClonesArray("AliMFTDigit");
855
856 }
857
858 //====================================================================================================================================================
859
860 void AliMFT::CreateDigits() {
861
862   // create digits list
863
864   AliDebug(1, "AliMFT::CreateDigits()");
865
866   if (fDigitsPerPlane) return; 
867   fDigitsPerPlane = new TObjArray(fNPlanes);
868   for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) fDigitsPerPlane->AddAt(new TClonesArray("AliMFTDigit"), iPlane);
869
870 }
871
872 //====================================================================================================================================================
873
874 void AliMFT::CreateRecPoints() {
875
876   // create recPoints list
877
878   AliDebug(1, "AliMFT::CreateRecPoints()");
879
880   if (fRecPointsPerPlane) return; 
881   fRecPointsPerPlane = new TObjArray(fNPlanes);
882   for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) fRecPointsPerPlane->AddAt(new TClonesArray("AliMFTCluster"), iPlane);
883
884 }
885
886 //====================================================================================================================================================