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