]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MFT/AliMFT.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Support and analysis code modified
[u/mrichter/AliRoot.git] / MFT / AliMFT.cxx
CommitLineData
820b4d9e 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
48ClassImp(AliMFT)
49
50//====================================================================================================================================================
51
52AliMFT::AliMFT():
53 AliDetector(),
54 fVersion(1),
55 fNPlanes(0),
58b93f36 56 fNSlices(1),
820b4d9e 57 fSDigitsPerPlane(0),
58 fDigitsPerPlane(0),
59 fRecPointsPerPlane(0),
60 fSideDigits(0),
61 fSegmentation(0),
62 fNameGeomFile(0),
58b93f36 63 fChargeDispersion(25.e-4),
820b4d9e 64 fSingleStepForChargeDispersion(0),
58b93f36 65 fNStepForChargeDispersion(4),
66 fDensitySupportOverSi(0.15)
820b4d9e 67{
68
69 // default constructor
70
71}
72
73//====================================================================================================================================================
74
75AliMFT::AliMFT(const Char_t *name, const Char_t *title):
76 AliDetector(name, title),
77 fVersion(1),
78 fNPlanes(0),
58b93f36 79 fNSlices(1),
820b4d9e 80 fSDigitsPerPlane(0),
81 fDigitsPerPlane(0),
82 fRecPointsPerPlane(0),
83 fSideDigits(0),
84 fSegmentation(0),
85 fNameGeomFile(0),
58b93f36 86 fChargeDispersion(25.e-4),
820b4d9e 87 fSingleStepForChargeDispersion(0),
58b93f36 88 fNStepForChargeDispersion(4),
89 fDensitySupportOverSi(0.15)
820b4d9e 90{
91
92 fNameGeomFile = "AliMFTGeometry.root";
93
94 SetGeometry();
95
96 Init();
97
98}
99
100//====================================================================================================================================================
101
102AliMFT::AliMFT(const Char_t *name, const Char_t *title, Char_t *nameGeomFile):
103 AliDetector(name, title),
104 fVersion(1),
105 fNPlanes(0),
58b93f36 106 fNSlices(1),
820b4d9e 107 fSDigitsPerPlane(0),
108 fDigitsPerPlane(0),
109 fRecPointsPerPlane(0),
110 fSideDigits(0),
111 fSegmentation(0),
112 fNameGeomFile(0),
58b93f36 113 fChargeDispersion(25.e-4),
820b4d9e 114 fSingleStepForChargeDispersion(0),
58b93f36 115 fNStepForChargeDispersion(4),
116 fDensitySupportOverSi(0.15)
820b4d9e 117{
118
119 fNameGeomFile = nameGeomFile;
120
121 SetGeometry();
122
123 Init();
124
125}
126
127//====================================================================================================================================================
128
129AliMFT::~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
139void AliMFT::CreateMaterials() {
140
141 // Definition of MFT materials - to be updated to the most recent values
142
143 AliInfo("Start MFT materials");
144
026547c6 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
820b4d9e 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
026547c6 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]
820b4d9e 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
026547c6 189 AliMixture(++matId,"Air", aAir, zAir, dAir, nAir, wAir);
820b4d9e 190 AliMedium(kAir, "Air", matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
191
026547c6 192 AliMaterial(++matId, "Si", aSi, zSi, dSi, radSi, absSi);
820b4d9e 193 AliMedium(kSi, "Si", matId, sens, isxfld, sxmgmx, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
194
026547c6 195 AliMaterial(++matId, "Readout", aSi, zSi, dSi, radSi, absSi);
820b4d9e 196 AliMedium(kReadout, "Readout", matId, unsens, isxfld, sxmgmx, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
197
026547c6 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);
4cc511f4 200
026547c6 201 AliMaterial(++matId, "Carbon", aCarb, zCarb, dCarb, radCarb, absCarb );
202 AliMedium(kCarbon, "Carbon", matId, unsens, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
4cc511f4 203
026547c6 204 AliMaterial(++matId, "Alu", aAlu, zAlu, dAlu, radAlu, absAlu);
205 AliMedium(kAlu, "Alu", matId, unsens, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
4cc511f4 206
026547c6 207 AliMixture(++matId, "Water", aWater, zWater, dWater, nWater, wWater);
208 AliMedium(kWater, "Water", matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
4cc511f4 209
026547c6 210 AliMixture(++matId, "SiO2", aSiO2, zSiO2, dSiO2, nSiO2, wSiO2);
4cc511f4 211 AliMedium(kSiO2, "SiO2", matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
026547c6 212
213 AliMixture(++matId, "Inox", aInox, zInox, dInox, nInox, wInox);
214 AliMedium(kInox, "Inox", matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
4cc511f4 215
820b4d9e 216 AliInfo("End MFT materials");
217
218}
219
220//====================================================================================================================================================
221
222void AliMFT::CreateGeometry() {
223
224 // Creates detailed geometry simulation (currently GEANT volumes tree)
225
226 AliInfo("Start MFT preliminary version building");
227 if(!gMC->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
a2b7dc2a 239void 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
820b4d9e 256void AliMFT::StepManager() {
257
258 // Full Step Manager
259
d4643a10 260 AliDebug(2, Form("Entering StepManager: gMC->CurrentVolName() = %s", gMC->CurrentVolName()));
261
820b4d9e 262 if (!fSegmentation) AliFatal("No segmentation available"); // DO WE HAVE A SEGMENTATION???
263
264 if (!(this->IsActive())) return;
265 if (!(gMC->TrackCharge())) return;
266
267 TString planeNumber = gMC->CurrentVolName();
268 TString detElemNumber = gMC->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);
026547c6 275 Int_t detElemID = fSegmentation->GetDetElemGlobalID(planeNumber.Atoi(), detElemNumber.Atoi());
820b4d9e 276
277 if (gMC->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 (gMC->IsTrackInside()) status += 1;
288 if (gMC->IsTrackEntering()) status += 2;
289 if (gMC->IsTrackExiting()) status += 4;
290 if (gMC->IsTrackOut()) status += 8;
291 if (gMC->IsTrackDisappeared()) status += 16;
292 if (gMC->IsTrackStop()) status += 32;
293 if (gMC->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 gMC->TrackPosition(position);
302 gMC->TrackMomentum(momentum);
303
d4643a10 304 AliDebug(1, Form("AliMFT::StepManager()->%s Hit #%06d (x=%f, y=%f, z=%f) belongs to track %02d\n",
305 gMC->CurrentVolName(), fNhits, position.X(), position.Y(), position.Z(), gAlice->GetMCApp()->GetCurrentTrackNumber()));
820b4d9e 306
307 hit.SetPosition(position);
308 hit.SetTOF(gMC->TrackTime());
309 hit.SetMomentum(momentum);
310 hit.SetStatus(status);
311 hit.SetEloss(gMC->Edep());
026547c6 312 // hit.SetShunt(GetIshunt());
313// if (gMC->IsTrackEntering()) {
314// hit.SetStartPosition(position);
315// hit.SetStartTime(gMC->TrackTime());
316// hit.SetStartStatus(status);
317// return; // don't save entering hit.
318// }
820b4d9e 319
320 // Fill hit structure with this new hit.
321 new ((*fHits)[fNhits++]) AliMFTHit(hit);
322
026547c6 323 // Save old position... for next hit.
324// hit.SetStartPosition(position);
325// hit.SetStartTime(gMC->TrackTime());
326// hit.SetStartStatus(status);
327
820b4d9e 328 return;
329
330}
331
332//====================================================================================================================================================
333
334TGeoVolumeAssembly* 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");
026547c6 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");
4cc511f4 349
026547c6 350 // ---- Cage & Services Description --------------------------------------------
4cc511f4 351 // R. Tieulent - 17/01/2014 - Basic description for ITS/TPC matching studies
026547c6 352 // -----------------------------------------------------------------------------
4cc511f4
AU		 353
354 TGeoVolumeAssembly *cageNservices = new TGeoVolumeAssembly("MFT_cageNservices");
355
356 // cage definition
026547c6 357 Float_t cageDz = 150./2.;
358 Float_t cageRMax = 49.5;
359 Float_t cageRMin = cageRMax - 0.120; // 0.64% of X0
4cc511f4 360
026547c6 361 TGeoVolume *cage = gGeoManager->MakeTube("MFT_cage", carbon, cageRMin, cageRMax, cageDz);
362 cage->SetLineColor(kBlue);
363
4cc511f4
AU		 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
026547c6 370 Float_t busBarDz = 150.;
371 Float_t busBarThick = 0.1 ;
372 Float_t busBarWidth = 1.;
4cc511f4 373
026547c6 374 TGeoVolume *aluBusBar = gGeoManager->MakeBox("MFT_busbar", alu, busBarWidth/2., busBarThick/2., busBarDz/2.);
375 aluBusBar->SetLineColor(kYellow);
4cc511f4
AU		 376
377 Int_t nBusBar = 30;
026547c6 378 Float_t dPhiBusBar = 2.*TMath::Pi() / nBusBar;
379 Float_t dShift = cageRMin - busBarThick;
4cc511f4 380
026547c6 381 TGeoRotation *rot;
4cc511f4
AU		 382
383 for (Int_t iBusBar=0; iBusBar<nBusBar; iBusBar++) {
026547c6 384 Float_t phi = dPhiBusBar*iBusBar;
385 Float_t xp = dShift*TMath::Cos(phi);
386 Float_t yp = dShift*TMath::Sin(phi);
4cc511f4 387 rot = new TGeoRotation();
026547c6 388 rot->RotateZ(phi*TMath::RadToDeg()+90.);
389 services->AddNode(aluBusBar, iBusBar+1, new TGeoCombiTrans(xp,yp,0,rot));
4cc511f4
AU		 390 }
391
026547c6 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
4cc511f4 398
026547c6 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 ));
4cc511f4
AU		 407
408 Int_t nCooling = 18;
026547c6 409 dShift = cageRMin - coolingR ;
4cc511f4
AU		 410 Float_t phi0 = 0.02;
411
412 for (Int_t iCooling=0; iCooling<nCooling; iCooling++) {
026547c6 413 Float_t phi ;
414 if (iCooling<nCooling/2) phi = dPhiBusBar*(iCooling+3) + phi0;
415 else phi = dPhiBusBar*(iCooling+9) + phi0;
4cc511f4
AU		 416 Float_t xp = dShift*TMath::Cos(phi);
417 Float_t yp = dShift*TMath::Sin(phi);
026547c6 418 services->AddNode(cooling, iCooling+1, new TGeoTranslation(xp,yp,0 ));
4cc511f4
AU		 419 }
420
421 // Optical Fibers
026547c6 422 Float_t fiberDz = 150.;
423 Float_t fiberRadius = 0.0125 /2. ; // 0.125mm in diameter
4cc511f4 424
026547c6 425 TGeoVolume *fiber = gGeoManager->MakeTube("MFT_fiber", si02, 0., fiberRadius, fiberDz/2.);
426 fiber->SetLineColor(kCyan);
427
4cc511f4 428 Int_t nFiber = 340;
026547c6 429 dShift = cageRMin - 2*fiberRadius;
4cc511f4
AU		 430 phi0 = 0.03;
431
432 for (Int_t iFiber=0; iFiber<nFiber; iFiber++) {
026547c6 433 Float_t phi = dPhiBusBar*(Int_t)(iFiber/11) - phi0-(iFiber%11)*2.*TMath::ATan(fiberRadius/dShift);
4cc511f4
AU		 434 Float_t xp = dShift*TMath::Cos(phi);
435 Float_t yp = dShift*TMath::Sin(phi);
026547c6 436 services->AddNode(fiber, iFiber+1, new TGeoTranslation(xp,yp,0 ));
4cc511f4
AU		 437 }
438
026547c6 439 cageNservices->AddNode(services, 1, new TGeoTranslation(0., 0., 0. ));
440
4cc511f4
AU		 441 vol->AddNode(cageNservices,1,new TGeoTranslation(0., 0., 0. ));
442
026547c6 443 // ------------------- Creating volumes for MFT planes --------------------------------
820b4d9e 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()) );
d4643a10 463 AliDebug(2, Form("Created vol %s", supportElem->GetName()));
464 supportElem->SetLineColor(kCyan-9);
820b4d9e 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);
d4643a10 484 AliDebug(2, Form("Created vol %s", activeElem->GetName()));
485 activeElem->SetLineColor(kGreen);
820b4d9e 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);
d4643a10 506 AliDebug(2, Form("Created vol %s", readoutElem->GetName()));
507 readoutElem->SetLineColor(kRed);
820b4d9e 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
522void 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
564void AliMFT::Hits2SDigitsLocal(TClonesArray *hits, const TObjArray *pSDig, Int_t track) {
565
566 // Add sdigits of these hits to the list
567
58b93f36 568 AliDebug(1, "Entering Hits2SDigitsLocal");
820b4d9e 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
58b93f36 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());
820b4d9e 591
592 Int_t xPixel = -1;
593 Int_t yPixel = -1;
58b93f36 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);
820b4d9e 603 AliDebug(1, Form("Created new sdigit (%f, %f, %f) from hit (%f, %f, %f)",
58b93f36 604 mainSDigit->GetPixelCenterX(), mainSDigit->GetPixelCenterY(), mainSDigit->GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
820b4d9e 605 }
606
58b93f36 607 // creating "side digits" to simulate the effect of charge dispersion
820b4d9e 608
820b4d9e 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);
58b93f36 615 if (fSegmentation->Hit2PixelID(hit->X()+shiftX, hit->Y()+shiftY, hit->GetDetElemID(), xPixel, yPixel)) {
820b4d9e 616 Bool_t digitExists = kFALSE;
58b93f36 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;
820b4d9e 622 }
623 }
624 if (!digitExists) {
58b93f36 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);
820b4d9e 638 }
639 }
640 }
641 }
58b93f36 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));
bcaf50eb 648 if (mySDig->GetPixelX()%2 == mySDig->GetPixelY()%2) { // both pair or both odd
58b93f36 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);
bcaf50eb 664 }
665 }
666 else { // pair-odd
58b93f36 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);
bcaf50eb 682 }
683 }
684 }
685 }
58b93f36 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
bcaf50eb 730 }
731
58b93f36 732 AliDebug(1,"Exiting Hits2SDigitsLocal");
820b4d9e 733
734}
735
736//====================================================================================================================================================
737
738void AliMFT::MakeBranch(Option_t *option) {
739
820b4d9e 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
776void 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
819void AliMFT::SetGeometry() {
820
ffc53def 821 AliInfo("AliMFT::SetGeometry\n");
820b4d9e 822
823 fSegmentation = new AliMFTSegmentation(fNameGeomFile.Data());
824
825 fNPlanes = fSegmentation->GetNPlanes();
826
827}
828
829//====================================================================================================================================================
830
831void 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
844void 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
860void 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
874void 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//====================================================================================================================================================