Ignoring coverity to compile
[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   aAir[4]={12,14,16,36} ,   zAir[4]={6,7,8,18} ,   wAir[4]={0.000124,0.755267,0.231781,0.012827} , dAir=0.00120479; Int_t nAir=4;   // Air mixture
147   Float_t   aSi = 28.085 ,            zSi   = 14 ,           dSi   =  2.329    ,   radSi   =  21.82/dSi , absSi   = 108.4/dSi  ;              // Silicon
148
149   Int_t   matId  = 0;                        // tmp material id number
150   Int_t   unsens = 0, sens=1;                // sensitive or unsensitive medium
151   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
152   Float_t maxfld = 5.;                       // max field value
153
154   Float_t tmaxfd = -10.0;                    // max deflection angle due to magnetic field in one step
155   Float_t stemax =  0.001;                   // max step allowed [cm]
156   Float_t deemax = -0.2;                     // maximum fractional energy loss in one step 0<deemax<=1  
157   Float_t epsil  =  0.001;                   // tracking precision [cm]   
158   Float_t stmin  = -0.001;                   // minimum step due to continuous processes [cm] (negative value: choose it automatically)
159
160   Float_t tmaxfdSi =  0.1;                    // max deflection angle due to magnetic field in one step
161   Float_t stemaxSi =  5.0e-4;                 // maximum step allowed [cm]
162   Float_t deemaxSi =  0.1;                    // maximum fractional energy loss in one step 0<deemax<=1
163   Float_t epsilSi  =  0.5e-4;                 // tracking precision [cm]
164   Float_t stminSi  = -0.001;                 // minimum step due to continuous processes [cm] (negative value: choose it automatically)
165
166   Int_t   isxfld = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Integ();   // from CreateMaterials in STRUCT/AliPIPEv3.cxx
167   Float_t sxmgmx = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Max();     // from CreateMaterials in STRUCT/AliPIPEv3.cxx
168       
169   AliMixture(++matId,"Air", aAir,  zAir,   dAir,   nAir,   wAir); 
170   AliMedium(kAir,    "Air", matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
171   
172   AliMaterial(++matId, "Si", aSi,   zSi,  dSi,    radSi,  absSi  );  
173   AliMedium(kSi,       "Si", matId, sens, isxfld, sxmgmx, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
174
175   AliMaterial(++matId, "Readout", aSi,   zSi,    dSi,    radSi,  absSi  );  
176   AliMedium(kReadout,  "Readout", matId, unsens, isxfld, sxmgmx, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
177
178   AliMaterial(++matId, "Support", aSi,   zSi,    dSi*fDensitySupportOverSi, radSi/fDensitySupportOverSi, absSi/fDensitySupportOverSi);  
179   AliMedium(kSupport,  "Support", matId, unsens, isxfld,  sxmgmx,    tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
180     
181   AliInfo("End MFT materials");
182           
183 }
184
185 //====================================================================================================================================================
186
187 void AliMFT::CreateGeometry() {
188
189   // Creates detailed geometry simulation (currently GEANT volumes tree)        
190
191   AliInfo("Start MFT preliminary version building");
192   if(!gMC->IsRootGeometrySupported()) return;                
193   TGeoVolumeAssembly *vol = CreateVol();
194   AliInfo("TGeoVolumeAssembly created!");
195   gGeoManager->GetVolume("ALIC")->AddNode(vol,0);
196   AliInfo("Stop MFT preliminary version building");
197
198   if (fNStepForChargeDispersion) fSingleStepForChargeDispersion = fChargeDispersion/Double_t(fNStepForChargeDispersion);
199
200
201
202 //====================================================================================================================================================
203
204 void AliMFT::AddAlignableVolumes() {
205
206   // Create entries for alignable volumes associating the symbolic volume
207   // name with the corresponding volume path. Needs to be syncronized with
208   // eventual changes in the geometry.
209
210   TString sysName = "MFT";
211   TString volPath = "/ALIC_1/MFT_0";
212   
213   if (!gGeoManager->SetAlignableEntry(sysName.Data(),volPath.Data())) {
214     AliFatal(Form("Alignable entry %s not created. Volume path %s not valid", sysName.Data(), volPath.Data()));
215   }  
216
217 }
218
219 //====================================================================================================================================================
220
221 void AliMFT::StepManager() {
222
223   // Full Step Manager
224
225   AliDebug(2, Form("Entering StepManager: gMC->CurrentVolName() = %s", gMC->CurrentVolName()));
226
227   if (!fSegmentation) AliFatal("No segmentation available");    // DO WE HAVE A SEGMENTATION???
228
229   if (!(this->IsActive())) return;
230   if (!(gMC->TrackCharge())) return;
231
232   TString planeNumber   = gMC->CurrentVolName();
233   TString detElemNumber = gMC->CurrentVolName();
234   if (planeNumber.Contains("support")) return;
235   if (planeNumber.Contains("readout")) return;
236   planeNumber.Remove(0,9);
237   detElemNumber.Remove(0,19);
238   planeNumber.Remove(2);
239   detElemNumber.Remove(3);
240   Int_t detElemID = fSegmentation->GetDetElemID(planeNumber.Atoi(), detElemNumber.Atoi());
241
242   if (gMC->IsTrackExiting()) {
243     AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kMFT);
244   }
245
246   static TLorentzVector position, momentum;
247   static AliMFTHit hit;
248   
249   Int_t  status = 0;
250   
251   // Track status
252   if (gMC->IsTrackInside())      status +=  1;
253   if (gMC->IsTrackEntering())    status +=  2;
254   if (gMC->IsTrackExiting())     status +=  4;
255   if (gMC->IsTrackOut())         status +=  8;
256   if (gMC->IsTrackDisappeared()) status += 16;
257   if (gMC->IsTrackStop())        status += 32;
258   if (gMC->IsTrackAlive())       status += 64;
259
260   // ---------- Fill hit structure
261
262   hit.SetDetElemID(detElemID);
263   hit.SetPlane(planeNumber.Atoi());
264   hit.SetTrack(gAlice->GetMCApp()->GetCurrentTrackNumber());
265     
266   gMC->TrackPosition(position);
267   gMC->TrackMomentum(momentum);
268
269   AliDebug(1, Form("AliMFT::StepManager()->%s Hit #%06d (x=%f, y=%f, z=%f) belongs to track %02d\n", 
270                    gMC->CurrentVolName(), fNhits, position.X(), position.Y(), position.Z(), gAlice->GetMCApp()->GetCurrentTrackNumber())); 
271
272   hit.SetPosition(position);
273   hit.SetTOF(gMC->TrackTime());
274   hit.SetMomentum(momentum);
275   hit.SetStatus(status);
276   hit.SetEloss(gMC->Edep());
277   //  hit.SetShunt(GetIshunt());
278 //   if (gMC->IsTrackEntering()) {
279 //     hit.SetStartPosition(position);
280 //     hit.SetStartTime(gMC->TrackTime());
281 //     hit.SetStartStatus(status);
282 //     return; // don't save entering hit.
283 //   } 
284
285   // Fill hit structure with this new hit.
286   new ((*fHits)[fNhits++]) AliMFTHit(hit);
287
288   // Save old position... for next hit.
289 //   hit.SetStartPosition(position);
290 //   hit.SetStartTime(gMC->TrackTime());
291 //   hit.SetStartStatus(status);
292
293   return;
294
295 }
296
297 //====================================================================================================================================================
298
299 TGeoVolumeAssembly* AliMFT::CreateVol() {
300
301   // method to create the MFT Geometry (silicon circular planes)
302
303   if (!fSegmentation) CreateGeometry();
304   
305   TGeoVolumeAssembly *vol = new TGeoVolumeAssembly("MFT");
306   TGeoMedium *silicon = gGeoManager->GetMedium("MFT_Si");
307   TGeoMedium *readout = gGeoManager->GetMedium("MFT_Readout");
308   TGeoMedium *support = gGeoManager->GetMedium("MFT_Support");
309   for (Int_t iPar=0; iPar<8; iPar++) AliInfo(Form("silicon->GetParam(%d) = %f", iPar, silicon->GetParam(iPar)));
310   for (Int_t iPar=0; iPar<8; iPar++) AliInfo(Form("readout->GetParam(%d) = %f", iPar, readout->GetParam(iPar)));
311   for (Int_t iPar=0; iPar<8; iPar++) AliInfo(Form("support->GetParam(%d) = %f", iPar, support->GetParam(iPar)));
312
313   Double_t origin[3] = {0};
314
315   for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
316
317     AliDebug(1, Form("Creating volumes for MFT plane %02d",iPlane));
318
319     AliMFTPlane *plane = fSegmentation->GetPlane(iPlane);
320
321     // --------- support element(s)
322     
323     origin[0] =  0.5*(plane->GetSupportElement(0)->GetAxis(0)->GetXmax() + plane->GetSupportElement(0)->GetAxis(0)->GetXmin());
324     origin[1] =  0.5*(plane->GetSupportElement(0)->GetAxis(1)->GetXmax() + plane->GetSupportElement(0)->GetAxis(1)->GetXmin());
325     origin[2] = -0.5*(plane->GetSupportElement(0)->GetAxis(2)->GetXmax() + plane->GetSupportElement(0)->GetAxis(2)->GetXmin());
326     TGeoVolume *supportElem = gGeoManager->MakeTube(Form("MFT_plane%02d_support", iPlane), support, 
327                                                     plane->GetRMinSupport(), 
328                                                     plane->GetRMaxSupport(),
329                                                     0.5*(plane->GetSupportElement(0)->GetAxis(2)->GetXmax() - 
330                                                          plane->GetSupportElement(0)->GetAxis(2)->GetXmin()) );
331     AliDebug(2, Form("Created vol %s", supportElem->GetName()));
332     supportElem->SetLineColor(kCyan-9);
333     vol -> AddNode(supportElem, 0, new TGeoTranslation(origin[0], origin[1], origin[2]));
334
335     AliDebug(1, "support elements created!");
336     
337     // --------- active elements
338
339     for (Int_t iActive=0; iActive<plane->GetNActiveElements(); iActive++) {
340       
341       Double_t dx = 0.5*TMath::Abs(plane->GetActiveElement(iActive)->GetAxis(0)->GetXmax() - plane->GetActiveElement(iActive)->GetAxis(0)->GetXmin());
342       Double_t dy = 0.5*TMath::Abs(plane->GetActiveElement(iActive)->GetAxis(1)->GetXmax() - plane->GetActiveElement(iActive)->GetAxis(1)->GetXmin());
343       Double_t dz = 0.5*TMath::Abs(plane->GetActiveElement(iActive)->GetAxis(2)->GetXmax() - plane->GetActiveElement(iActive)->GetAxis(2)->GetXmin());
344       dz /= Double_t(fNSlices);
345
346       origin[0] =  0.5*(plane->GetActiveElement(iActive)->GetAxis(0)->GetXmax() + plane->GetActiveElement(iActive)->GetAxis(0)->GetXmin());
347       origin[1] =  0.5*(plane->GetActiveElement(iActive)->GetAxis(1)->GetXmax() + plane->GetActiveElement(iActive)->GetAxis(1)->GetXmin());
348
349       for (Int_t iSlice=0; iSlice<fNSlices; iSlice++) {
350         origin[2] = -0.5*(plane->GetActiveElement(iActive)->GetAxis(2)->GetXmin() + 2*dz*(iSlice+1) + plane->GetActiveElement(iActive)->GetAxis(2)->GetXmin() + 2*dz*(iSlice) );
351         TGeoVolume *activeElem = gGeoManager->MakeBox(Form("MFT_plane%02d_active%03d_slice%02d", iPlane, iActive, iSlice), silicon, dx, dy, dz);
352         AliDebug(2, Form("Created vol %s", activeElem->GetName()));
353         activeElem->SetLineColor(kGreen);
354         vol -> AddNode(activeElem, 0, new TGeoTranslation(origin[0], origin[1], origin[2]));
355       }
356
357     }
358
359     AliDebug(1, "active elements created!");
360
361     // --------- readout elements
362
363     for (Int_t iReadout=0; iReadout<plane->GetNReadoutElements(); iReadout++) {
364       
365       Double_t dx = 0.5*TMath::Abs(plane->GetReadoutElement(iReadout)->GetAxis(0)->GetXmax() - plane->GetReadoutElement(iReadout)->GetAxis(0)->GetXmin());
366       Double_t dy = 0.5*TMath::Abs(plane->GetReadoutElement(iReadout)->GetAxis(1)->GetXmax() - plane->GetReadoutElement(iReadout)->GetAxis(1)->GetXmin());
367       Double_t dz = 0.5*TMath::Abs(plane->GetReadoutElement(iReadout)->GetAxis(2)->GetXmax() - plane->GetReadoutElement(iReadout)->GetAxis(2)->GetXmin());
368
369       origin[0] =  0.5*(plane->GetReadoutElement(iReadout)->GetAxis(0)->GetXmax() + plane->GetReadoutElement(iReadout)->GetAxis(0)->GetXmin());
370       origin[1] =  0.5*(plane->GetReadoutElement(iReadout)->GetAxis(1)->GetXmax() + plane->GetReadoutElement(iReadout)->GetAxis(1)->GetXmin());
371       origin[2] = -0.5*(plane->GetReadoutElement(iReadout)->GetAxis(2)->GetXmax() + plane->GetReadoutElement(iReadout)->GetAxis(2)->GetXmin());
372
373       TGeoVolume *readoutElem = gGeoManager->MakeBox(Form("MFT_plane%02d_readout%03d", iPlane, iReadout), readout, dx, dy, dz);
374       AliDebug(2, Form("Created vol %s", readoutElem->GetName()));
375       readoutElem->SetLineColor(kRed);
376       vol -> AddNode(readoutElem, 0, new TGeoTranslation(origin[0], origin[1], origin[2]));
377       
378     }
379
380     AliDebug(1, "readout elements created!");
381
382   }
383
384   return vol;
385
386 }
387
388 //====================================================================================================================================================
389
390 void AliMFT::Hits2SDigits(){
391   
392   // Interface method invoked from AliSimulation to create a list of sdigits corresponding to list of hits. Every hit generates one sdigit.
393
394   AliDebug(1,"Start Hits2SDigits.");
395   
396   if (!fSegmentation) CreateGeometry();
397  
398   if (!fLoader->TreeH()) fLoader->LoadHits();
399
400   if (!fLoader->TreeS()) {
401
402     for (Int_t iEvt=0;iEvt<fLoader->GetRunLoader()->GetNumberOfEvents(); iEvt++) {
403
404       fLoader->GetRunLoader()->GetEvent(iEvt);
405       fLoader->MakeTree("S");
406       MakeBranch("S");
407       SetTreeAddress();
408
409       AliDebug(1, Form("Event %03d: fLoader->TreeH()->GetEntries() = %2d", iEvt, Int_t(fLoader->TreeH()->GetEntries())));
410
411       for (Int_t iTrack=0; iTrack<fLoader->TreeH()->GetEntries(); iTrack++) {
412         fLoader->TreeH()->GetEntry(iTrack);     
413         Hits2SDigitsLocal(Hits(), GetSDigitsList(), iTrack);    // convert these hits to a list of sdigits  
414       }
415       
416       fLoader->TreeS()->Fill();
417       fLoader->WriteSDigits("OVERWRITE");
418       ResetSDigits();
419
420     }
421   }
422
423   fLoader->UnloadHits();
424   fLoader->UnloadSDigits();
425
426   AliDebug(1,"Stop Hits2SDigits.");
427   
428 }
429
430 //====================================================================================================================================================
431
432 void AliMFT::Hits2SDigitsLocal(TClonesArray *hits, const TObjArray *pSDig, Int_t track) {
433
434   //  Add sdigits of these hits to the list
435   
436   AliDebug(1, "Entering Hits2SDigitsLocal");
437   
438   if (!fSegmentation) CreateGeometry();
439   
440   TClonesArray *pSDigList[fNMaxPlanes];
441   for (Int_t iPlane=0; iPlane<fNMaxPlanes; iPlane++) pSDigList[iPlane] = NULL; 
442   for (Int_t iPlane=0; iPlane<fNPlanes;    iPlane++) { 
443     pSDigList[iPlane] = (TClonesArray*) (*pSDig)[iPlane];
444     AliDebug(1,Form("Entries of pSDigList %3d; plane: %02d,",pSDigList[iPlane]->GetEntries(),iPlane));
445     if (!track && pSDigList[iPlane]->GetEntries()!=0) AliErrorClass("Some of sdigits lists is not empty");
446   }
447   
448   for (Int_t iHit=0; iHit<hits->GetEntries(); iHit++) {
449
450     AliMFTHit *hit = (AliMFTHit*) hits->At(iHit);
451
452     // Creating "main digit"
453
454     AliMFTDigit *mainSDigit = new AliMFTDigit();
455     mainSDigit->SetEloss(hit->GetEloss());
456     mainSDigit->SetDetElemID(hit->GetDetElemID());
457     mainSDigit->SetPlane(hit->GetPlane());
458     mainSDigit->AddMCLabel(hit->GetTrack()); 
459
460     Int_t xPixel = -1;
461     Int_t yPixel = -1;
462     if (fSegmentation->Hit2PixelID(hit->X(), hit->Y(), mainSDigit->GetDetElemID(), xPixel, yPixel)) {
463       mainSDigit->SetPixID(xPixel, yPixel, 0);
464       mainSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(mainSDigit->GetDetElemID()), 
465                           fSegmentation->GetPixelSizeY(mainSDigit->GetDetElemID()),
466                           fSegmentation->GetPixelSizeZ(mainSDigit->GetDetElemID()));  
467       mainSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(mainSDigit->GetDetElemID(), xPixel), 
468                            fSegmentation->GetPixelCenterY(mainSDigit->GetDetElemID(), yPixel),
469                            fSegmentation->GetPixelCenterZ(mainSDigit->GetDetElemID(), 0));
470       new ((*fSideDigits)[fSideDigits->GetEntries()]) AliMFTDigit(*mainSDigit);
471       AliDebug(1, Form("Created new sdigit (%f, %f, %f) from hit (%f, %f, %f)",
472                        mainSDigit->GetPixelCenterX(), mainSDigit->GetPixelCenterY(), mainSDigit->GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
473     }
474
475     // creating "side digits" to simulate the effect of charge dispersion
476
477     Double_t pi4 = TMath::Pi()/4.;
478     for (Int_t iStep=0; iStep<fNStepForChargeDispersion; iStep++) {
479       Double_t shift = (iStep+1) * fSingleStepForChargeDispersion;
480       for (Int_t iAngle=0; iAngle<8; iAngle++) {
481         Double_t shiftX = shift*TMath::Cos(iAngle*pi4);
482         Double_t shiftY = shift*TMath::Sin(iAngle*pi4);
483         if (fSegmentation->Hit2PixelID(hit->X()+shiftX, hit->Y()+shiftY, hit->GetDetElemID(), xPixel, yPixel)) {
484           Bool_t digitExists = kFALSE;
485           for (Int_t iSideDigit=0; iSideDigit<fSideDigits->GetEntries(); iSideDigit++) {
486             if (xPixel==((AliMFTDigit*) fSideDigits->At(iSideDigit))->GetPixelX() && 
487                 yPixel==((AliMFTDigit*) fSideDigits->At(iSideDigit))->GetPixelY()) {
488               digitExists = kTRUE;
489               break;
490             }
491           }
492           if (!digitExists) {
493             AliMFTDigit *sideSDigit = new AliMFTDigit();
494             sideSDigit->SetEloss(0.);
495             sideSDigit->SetDetElemID(hit->GetDetElemID());
496             sideSDigit->SetPlane(hit->GetPlane());
497             sideSDigit->AddMCLabel(hit->GetTrack());
498             sideSDigit->SetPixID(xPixel, yPixel, 0);
499             sideSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(sideSDigit->GetDetElemID()), 
500                                     fSegmentation->GetPixelSizeY(sideSDigit->GetDetElemID()),
501                                     fSegmentation->GetPixelSizeZ(sideSDigit->GetDetElemID()));  
502             sideSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(sideSDigit->GetDetElemID(), xPixel), 
503                                      fSegmentation->GetPixelCenterY(sideSDigit->GetDetElemID(), yPixel),
504                                      fSegmentation->GetPixelCenterZ(sideSDigit->GetDetElemID(), 0)); 
505             new ((*fSideDigits)[fSideDigits->GetEntries()]) AliMFTDigit(*sideSDigit);
506           }
507         }
508       }
509     }
510     
511     // ------------ In case we should simulate a rectangular pattern of pixel...
512     
513     if (fSegmentation->GetPlane(mainSDigit->GetPlane())->HasPixelRectangularPatternAlongY()) {
514       for (Int_t iSDigit=0; iSDigit<fSideDigits->GetEntries(); iSDigit++) {
515         AliMFTDigit *mySDig = (AliMFTDigit*) (fSideDigits->At(iSDigit));
516         if (mySDig->GetPixelX()%2 == mySDig->GetPixelY()%2) {   // both pair or both odd
517           xPixel = mySDig->GetPixelX();
518           yPixel = mySDig->GetPixelY()+1;
519           if (fSegmentation->DoesPixelExist(mySDig->GetDetElemID(), xPixel, yPixel)) {
520             AliMFTDigit *newSDigit = new AliMFTDigit();
521             newSDigit->SetEloss(0.);
522             newSDigit->SetDetElemID(mySDig->GetDetElemID());
523             newSDigit->SetPlane(mySDig->GetDetElemID());
524             newSDigit->SetPixID(xPixel, yPixel, 0);
525             newSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(newSDigit->GetDetElemID()), 
526                                    fSegmentation->GetPixelSizeY(newSDigit->GetDetElemID()),
527                                    fSegmentation->GetPixelSizeZ(newSDigit->GetDetElemID()));  
528             newSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(newSDigit->GetDetElemID(), xPixel), 
529                                     fSegmentation->GetPixelCenterY(newSDigit->GetDetElemID(), yPixel),
530                                     fSegmentation->GetPixelCenterZ(newSDigit->GetDetElemID(), 0)); 
531             new ((*fSideDigits)[fSideDigits->GetEntries()]) AliMFTDigit(*newSDigit);
532           }
533         }
534         else {   // pair-odd
535           xPixel = mySDig->GetPixelX();
536           yPixel = mySDig->GetPixelY()-1;
537           if (fSegmentation->DoesPixelExist(mySDig->GetDetElemID(), xPixel, yPixel)) {
538             AliMFTDigit *newSDigit = new AliMFTDigit();
539             newSDigit->SetEloss(0.);
540             newSDigit->SetDetElemID(mySDig->GetDetElemID());
541             newSDigit->SetPlane(mySDig->GetPlane());
542             newSDigit->SetPixID(xPixel, yPixel, 0);
543             newSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(newSDigit->GetDetElemID()), 
544                                    fSegmentation->GetPixelSizeY(newSDigit->GetDetElemID()),
545                                    fSegmentation->GetPixelSizeZ(newSDigit->GetDetElemID()));  
546             newSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(newSDigit->GetDetElemID(), xPixel), 
547                                     fSegmentation->GetPixelCenterY(newSDigit->GetDetElemID(), yPixel),
548                                     fSegmentation->GetPixelCenterZ(newSDigit->GetDetElemID(), 0)); 
549             new ((*fSideDigits)[fSideDigits->GetEntries()]) AliMFTDigit(*newSDigit);
550           }
551         }
552       }
553     }
554
555     // -------- checking which pixels switched on have their diode actually within the charge dispersion radius
556
557     for (Int_t iSDigit=0; iSDigit<fSideDigits->GetEntries(); iSDigit++) {
558       AliMFTDigit *mySDig = (AliMFTDigit*) (fSideDigits->At(iSDigit));
559       Double_t distance = TMath::Sqrt(TMath::Power(mySDig->GetPixelCenterX()-hit->X(),2) + TMath::Power(mySDig->GetPixelCenterY()-hit->Y(),2));
560       if (fSegmentation->GetPlane(mySDig->GetPlane())->HasPixelRectangularPatternAlongY()) {
561         if (mySDig->GetPixelX()%2 == mySDig->GetPixelY()%2) {  // both pair or both odd
562           if (distance<fChargeDispersion) {
563             AliDebug(1, Form("Created new sdigit (%f, %f, %f) from hit (%f, %f, %f)",
564                              mySDig->GetPixelCenterX(), mySDig->GetPixelCenterY(), mySDig->GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
565             new ((*pSDigList[mySDig->GetPlane()])[pSDigList[mySDig->GetPlane()]->GetEntries()]) AliMFTDigit(*mySDig);
566             xPixel = mySDig->GetPixelX();
567             yPixel = mySDig->GetPixelY()+1;
568             if (fSegmentation->DoesPixelExist(mySDig->GetDetElemID(), xPixel, yPixel)) {
569               AliMFTDigit *newSDigit = new AliMFTDigit();
570               newSDigit->SetEloss(0.);
571               newSDigit->SetDetElemID(mySDig->GetDetElemID());
572               newSDigit->SetPlane(mySDig->GetPlane());
573               newSDigit->SetPixID(xPixel, yPixel, 0);
574               newSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(newSDigit->GetDetElemID()), 
575                                      fSegmentation->GetPixelSizeY(newSDigit->GetDetElemID()),
576                                      fSegmentation->GetPixelSizeZ(newSDigit->GetDetElemID()));  
577               newSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(newSDigit->GetDetElemID(), xPixel), 
578                                       fSegmentation->GetPixelCenterY(newSDigit->GetDetElemID(), yPixel),
579                                       fSegmentation->GetPixelCenterZ(newSDigit->GetDetElemID(), 0));
580               AliDebug(1, Form("Created new sdigit (%f, %f, %f) from hit (%f, %f, %f)",
581                                newSDigit->GetPixelCenterX(), newSDigit->GetPixelCenterY(), newSDigit->GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
582               new ((*pSDigList[newSDigit->GetPlane()])[pSDigList[newSDigit->GetPlane()]->GetEntries()]) AliMFTDigit(*newSDigit);
583             }
584           }
585         }
586       }
587       else {
588         if (distance<fChargeDispersion) {
589           AliDebug(1, Form("Created new sdigit (%f, %f, %f) from hit (%f, %f, %f)",
590                            mySDig->GetPixelCenterX(), mySDig->GetPixelCenterY(), mySDig->GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
591           new ((*pSDigList[mySDig->GetPlane()])[pSDigList[mySDig->GetPlane()]->GetEntries()]) AliMFTDigit(*mySDig);
592         }
593       }
594     }
595
596     fSideDigits->Delete(); 
597
598   }
599
600   AliDebug(1,"Exiting Hits2SDigitsLocal");
601
602 }
603
604 //====================================================================================================================================================
605
606 void AliMFT::MakeBranch(Option_t *option) {
607
608   // Create Tree branches 
609   AliDebug(1, Form("Start with option= %s.",option));
610   
611   const Int_t kBufSize = 4000;
612   
613   const Char_t *cH = strstr(option,"H");
614   const Char_t *cD = strstr(option,"D");
615   const Char_t *cS = strstr(option,"S");
616
617   if (cH && fLoader->TreeH()) {
618     CreateHits();
619     MakeBranchInTree(fLoader->TreeH(), "MFT", &fHits, kBufSize, 0);   
620   }
621
622   if (cS && fLoader->TreeS()) {
623     CreateSDigits();
624     for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) MakeBranchInTree(fLoader->TreeS(), 
625                                                                     Form("Plane_%02d",iPlane), 
626                                                                     &((*fSDigitsPerPlane)[iPlane]), 
627                                                                     kBufSize, 0);
628   }
629   
630   if (cD && fLoader->TreeD()) {
631     CreateDigits();
632     for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) MakeBranchInTree(fLoader->TreeD(), 
633                                                                     Form("Plane_%02d",iPlane), 
634                                                                     &((*fDigitsPerPlane)[iPlane]),
635                                                                     kBufSize, 0);
636   }
637
638   AliDebug(1,"Stop.");
639
640 }
641
642 //====================================================================================================================================================
643
644 void AliMFT::SetTreeAddress() {
645
646   AliDebug(1, "AliMFT::SetTreeAddress()");
647
648   //Set branch address for the Hits and Digits Tree.
649   AliDebug(1, "Start.");
650
651   AliDebug(1, Form("AliMFT::SetTreeAddress Hits  fLoader->TreeH() = %p\n", fLoader->TreeH()));
652   if (fLoader->TreeH() && fLoader->TreeH()->GetBranch("MFT")) {
653     CreateHits();
654     fLoader->TreeH()->SetBranchAddress("MFT", &fHits);
655   }
656     
657   AliDebug(1, Form("AliMFT::SetTreeAddress SDigits  fLoader->TreeS() = %p\n", fLoader->TreeS()));
658   if (fLoader->TreeS() && fLoader->TreeS()->GetBranch("Plane_00")) {
659     CreateSDigits();
660     for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
661       fLoader->TreeS()->SetBranchAddress(Form("Plane_%02d",iPlane), &((*fSDigitsPerPlane)[iPlane]));
662     }
663   }
664     
665   AliDebug(1, Form("AliMFT::SetTreeAddress Digits  fLoader->TreeD() = %p\n", fLoader->TreeD()));
666   if (fLoader->TreeD() && fLoader->TreeD()->GetBranch("Plane_00")) {
667     CreateDigits(); 
668     for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
669       fLoader->TreeD()->SetBranchAddress(Form("Plane_%02d",iPlane), &((*fDigitsPerPlane)[iPlane]));
670     }
671   }
672
673   AliDebug(1, Form("AliMFT::SetTreeAddress RecPoints  fLoader->TreeR() = %p\n", fLoader->TreeR()));
674   if (fLoader->TreeR() && fLoader->TreeR()->GetBranch("Plane_00")) {
675     CreateRecPoints(); 
676     for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
677       fLoader->TreeR()->SetBranchAddress(Form("Plane_%02d",iPlane), &((*fRecPointsPerPlane)[iPlane]));
678     }
679   }
680
681   AliDebug(1,"Stop.");
682
683 }
684
685 //====================================================================================================================================================
686
687 void AliMFT::SetGeometry() {
688
689   AliInfo("AliMFT::SetGeometry\n");
690
691   fSegmentation = new AliMFTSegmentation(fNameGeomFile.Data());
692
693   fNPlanes = fSegmentation->GetNPlanes();
694
695 }
696
697 //====================================================================================================================================================
698
699 void AliMFT::CreateHits() { 
700
701   // create array of hits
702
703   AliDebug(1, "AliMFT::CreateHits()");
704
705   if (fHits) return;    
706   fHits = new TClonesArray("AliMFTHit");  
707
708 }
709
710 //====================================================================================================================================================
711
712 void AliMFT::CreateSDigits() { 
713  
714   // create sdigits list
715
716   AliDebug(1, "AliMFT::CreateSDigits()");
717  
718   if (fSDigitsPerPlane) return; 
719   fSDigitsPerPlane = new TObjArray(fNPlanes); 
720   for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) fSDigitsPerPlane->AddAt(new TClonesArray("AliMFTDigit"), iPlane);
721
722   fSideDigits = new TClonesArray("AliMFTDigit");
723
724 }
725
726 //====================================================================================================================================================
727
728 void AliMFT::CreateDigits() {
729
730   // create digits list
731
732   AliDebug(1, "AliMFT::CreateDigits()");
733
734   if (fDigitsPerPlane) return; 
735   fDigitsPerPlane = new TObjArray(fNPlanes);
736   for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) fDigitsPerPlane->AddAt(new TClonesArray("AliMFTDigit"), iPlane);
737
738 }
739
740 //====================================================================================================================================================
741
742 void AliMFT::CreateRecPoints() {
743
744   // create recPoints list
745
746   AliDebug(1, "AliMFT::CreateRecPoints()");
747
748   if (fRecPointsPerPlane) return; 
749   fRecPointsPerPlane = new TObjArray(fNPlanes);
750   for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) fRecPointsPerPlane->AddAt(new TClonesArray("AliMFTCluster"), iPlane);
751
752 }
753
754 //====================================================================================================================================================