1 // **************************************************************************
2 // * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 // * Author: The ALICE Off-line Project. *
5 // * Contributors are mentioned in the code where appropriate. *
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 // **************************************************************************
16 //====================================================================================================================================================
18 // Main class of the ALICE Muon Forward Tracker
20 // Contact author: antonio.uras@cern.ch
22 //====================================================================================================================================================
26 #include "TGeoManager.h"
27 #include "TGeoVolume.h"
28 #include "TGeoMatrix.h"
29 #include "TVirtualMC.h"
30 #include "TClonesArray.h"
31 #include "TGeoGlobalMagField.h"
33 #include "AliLoader.h"
34 #include "AliDetector.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"
46 #include "TObjArray.h"
50 //====================================================================================================================================================
59 fRecPointsPerPlane(0),
63 fChargeDispersion(25.e-4),
64 fSingleStepForChargeDispersion(0),
65 fNStepForChargeDispersion(4),
66 fDensitySupportOverSi(0.15)
69 // default constructor
73 //====================================================================================================================================================
75 AliMFT::AliMFT(const Char_t *name, const Char_t *title):
76 AliDetector(name, title),
82 fRecPointsPerPlane(0),
86 fChargeDispersion(25.e-4),
87 fSingleStepForChargeDispersion(0),
88 fNStepForChargeDispersion(4),
89 fDensitySupportOverSi(0.15)
92 fNameGeomFile = "AliMFTGeometry.root";
100 //====================================================================================================================================================
102 AliMFT::AliMFT(const Char_t *name, const Char_t *title, Char_t *nameGeomFile):
103 AliDetector(name, title),
109 fRecPointsPerPlane(0),
113 fChargeDispersion(25.e-4),
114 fSingleStepForChargeDispersion(0),
115 fNStepForChargeDispersion(4),
116 fDensitySupportOverSi(0.15)
119 fNameGeomFile = nameGeomFile;
127 //====================================================================================================================================================
131 if (fSDigitsPerPlane) { fSDigitsPerPlane->Delete(); delete fSDigitsPerPlane; }
132 if (fDigitsPerPlane) { fDigitsPerPlane->Delete(); delete fDigitsPerPlane; }
133 if (fRecPointsPerPlane) { fRecPointsPerPlane->Delete(); delete fRecPointsPerPlane; }
137 //====================================================================================================================================================
139 void AliMFT::CreateMaterials() {
141 // Definition of MFT materials - to be updated to the most recent values
143 AliInfo("Start MFT materials");
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
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
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)
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)
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
169 AliMixture(++matId,"Air", aAir, zAir, dAir, nAir, wAir);
170 AliMedium(kAir, "Air", matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
172 AliMaterial(++matId, "Si", aSi, zSi, dSi, radSi, absSi );
173 AliMedium(kSi, "Si", matId, sens, isxfld, sxmgmx, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
175 AliMaterial(++matId, "Readout", aSi, zSi, dSi, radSi, absSi );
176 AliMedium(kReadout, "Readout", matId, unsens, isxfld, sxmgmx, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
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);
181 AliInfo("End MFT materials");
185 //====================================================================================================================================================
187 void AliMFT::CreateGeometry() {
189 // Creates detailed geometry simulation (currently GEANT volumes tree)
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");
198 if (fNStepForChargeDispersion) fSingleStepForChargeDispersion = fChargeDispersion/Double_t(fNStepForChargeDispersion);
202 //====================================================================================================================================================
204 void AliMFT::StepManager() {
208 AliDebug(2, Form("Entering StepManager: gMC->CurrentVolName() = %s", gMC->CurrentVolName()));
210 if (!fSegmentation) AliFatal("No segmentation available"); // DO WE HAVE A SEGMENTATION???
212 if (!(this->IsActive())) return;
213 if (!(gMC->TrackCharge())) return;
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());
225 if (gMC->IsTrackExiting()) {
226 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kMFT);
229 static TLorentzVector position, momentum;
230 static AliMFTHit hit;
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;
243 // ---------- Fill hit structure
245 hit.SetDetElemID(detElemID);
246 hit.SetPlane(planeNumber.Atoi());
247 hit.SetTrack(gAlice->GetMCApp()->GetCurrentTrackNumber());
249 gMC->TrackPosition(position);
250 gMC->TrackMomentum(momentum);
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()));
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.
268 // Fill hit structure with this new hit.
269 new ((*fHits)[fNhits++]) AliMFTHit(hit);
271 // Save old position... for next hit.
272 // hit.SetStartPosition(position);
273 // hit.SetStartTime(gMC->TrackTime());
274 // hit.SetStartStatus(status);
280 //====================================================================================================================================================
282 TGeoVolumeAssembly* AliMFT::CreateVol() {
284 // method to create the MFT Geometry (silicon circular planes)
286 if (!fSegmentation) CreateGeometry();
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)));
296 Double_t origin[3] = {0};
298 for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
300 AliDebug(1, Form("Creating volumes for MFT plane %02d",iPlane));
302 AliMFTPlane *plane = fSegmentation->GetPlane(iPlane);
304 // --------- support element(s)
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]));
318 AliDebug(1, "support elements created!");
320 // --------- active elements
322 for (Int_t iActive=0; iActive<plane->GetNActiveElements(); iActive++) {
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);
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());
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]));
342 AliDebug(1, "active elements created!");
344 // --------- readout elements
346 for (Int_t iReadout=0; iReadout<plane->GetNReadoutElements(); iReadout++) {
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());
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());
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]));
363 AliDebug(1, "readout elements created!");
371 //====================================================================================================================================================
373 void AliMFT::Hits2SDigits(){
375 // Interface method invoked from AliSimulation to create a list of sdigits corresponding to list of hits. Every hit generates one sdigit.
377 AliDebug(1,"Start Hits2SDigits.");
379 if (!fSegmentation) CreateGeometry();
381 if (!fLoader->TreeH()) fLoader->LoadHits();
383 if (!fLoader->TreeS()) {
385 for (Int_t iEvt=0;iEvt<fLoader->GetRunLoader()->GetNumberOfEvents(); iEvt++) {
387 fLoader->GetRunLoader()->GetEvent(iEvt);
388 fLoader->MakeTree("S");
392 AliDebug(1, Form("Event %03d: fLoader->TreeH()->GetEntries() = %2d", iEvt, Int_t(fLoader->TreeH()->GetEntries())));
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
399 fLoader->TreeS()->Fill();
400 fLoader->WriteSDigits("OVERWRITE");
406 fLoader->UnloadHits();
407 fLoader->UnloadSDigits();
409 AliDebug(1,"Stop Hits2SDigits.");
413 //====================================================================================================================================================
415 void AliMFT::Hits2SDigitsLocal(TClonesArray *hits, const TObjArray *pSDig, Int_t track) {
417 // Add sdigits of these hits to the list
419 AliDebug(1, "Entering Hits2SDigitsLocal");
421 if (!fSegmentation) CreateGeometry();
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");
431 for (Int_t iHit=0; iHit<hits->GetEntries(); iHit++) {
433 AliMFTHit *hit = (AliMFTHit*) hits->At(iHit);
435 // Creating "main digit"
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());
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()));
458 // creating "side digits" to simulate the effect of charge dispersion
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()) {
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);
494 // ------------ In case we should simulate a rectangular pattern of pixel...
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);
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);
538 // -------- checking which pixels switched on have their diode actually within the charge dispersion radius
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);
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);
579 fSideDigits->Delete();
583 AliDebug(1,"Exiting Hits2SDigitsLocal");
587 //====================================================================================================================================================
589 void AliMFT::MakeBranch(Option_t *option) {
591 printf("AliMFT::MakeBranch(...)\n");
593 // Create Tree branches
594 AliDebug(1, Form("Start with option= %s.",option));
596 const Int_t kBufSize = 4000;
598 const Char_t *cH = strstr(option,"H");
599 const Char_t *cD = strstr(option,"D");
600 const Char_t *cS = strstr(option,"S");
602 if (cH && fLoader->TreeH()) {
604 MakeBranchInTree(fLoader->TreeH(), "MFT", &fHits, kBufSize, 0);
607 if (cS && fLoader->TreeS()) {
609 for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) MakeBranchInTree(fLoader->TreeS(),
610 Form("Plane_%02d",iPlane),
611 &((*fSDigitsPerPlane)[iPlane]),
615 if (cD && fLoader->TreeD()) {
617 for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) MakeBranchInTree(fLoader->TreeD(),
618 Form("Plane_%02d",iPlane),
619 &((*fDigitsPerPlane)[iPlane]),
627 //====================================================================================================================================================
629 void AliMFT::SetTreeAddress() {
631 AliDebug(1, "AliMFT::SetTreeAddress()");
633 //Set branch address for the Hits and Digits Tree.
634 AliDebug(1, "Start.");
636 AliDebug(1, Form("AliMFT::SetTreeAddress Hits fLoader->TreeH() = %p\n", fLoader->TreeH()));
637 if (fLoader->TreeH() && fLoader->TreeH()->GetBranch("MFT")) {
639 fLoader->TreeH()->SetBranchAddress("MFT", &fHits);
642 AliDebug(1, Form("AliMFT::SetTreeAddress SDigits fLoader->TreeS() = %p\n", fLoader->TreeS()));
643 if (fLoader->TreeS() && fLoader->TreeS()->GetBranch("Plane_00")) {
645 for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
646 fLoader->TreeS()->SetBranchAddress(Form("Plane_%02d",iPlane), &((*fSDigitsPerPlane)[iPlane]));
650 AliDebug(1, Form("AliMFT::SetTreeAddress Digits fLoader->TreeD() = %p\n", fLoader->TreeD()));
651 if (fLoader->TreeD() && fLoader->TreeD()->GetBranch("Plane_00")) {
653 for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
654 fLoader->TreeD()->SetBranchAddress(Form("Plane_%02d",iPlane), &((*fDigitsPerPlane)[iPlane]));
658 AliDebug(1, Form("AliMFT::SetTreeAddress RecPoints fLoader->TreeR() = %p\n", fLoader->TreeR()));
659 if (fLoader->TreeR() && fLoader->TreeR()->GetBranch("Plane_00")) {
661 for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
662 fLoader->TreeR()->SetBranchAddress(Form("Plane_%02d",iPlane), &((*fRecPointsPerPlane)[iPlane]));
670 //====================================================================================================================================================
672 void AliMFT::SetGeometry() {
674 printf("AliMFT::SetGeometry\n");
676 fSegmentation = new AliMFTSegmentation(fNameGeomFile.Data());
678 fNPlanes = fSegmentation->GetNPlanes();
682 //====================================================================================================================================================
684 void AliMFT::CreateHits() {
686 // create array of hits
688 AliDebug(1, "AliMFT::CreateHits()");
691 fHits = new TClonesArray("AliMFTHit");
695 //====================================================================================================================================================
697 void AliMFT::CreateSDigits() {
699 // create sdigits list
701 AliDebug(1, "AliMFT::CreateSDigits()");
703 if (fSDigitsPerPlane) return;
704 fSDigitsPerPlane = new TObjArray(fNPlanes);
705 for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) fSDigitsPerPlane->AddAt(new TClonesArray("AliMFTDigit"), iPlane);
707 fSideDigits = new TClonesArray("AliMFTDigit");
711 //====================================================================================================================================================
713 void AliMFT::CreateDigits() {
715 // create digits list
717 AliDebug(1, "AliMFT::CreateDigits()");
719 if (fDigitsPerPlane) return;
720 fDigitsPerPlane = new TObjArray(fNPlanes);
721 for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) fDigitsPerPlane->AddAt(new TClonesArray("AliMFTDigit"), iPlane);
725 //====================================================================================================================================================
727 void AliMFT::CreateRecPoints() {
729 // create recPoints list
731 AliDebug(1, "AliMFT::CreateRecPoints()");
733 if (fRecPointsPerPlane) return;
734 fRecPointsPerPlane = new TObjArray(fNPlanes);
735 for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) fRecPointsPerPlane->AddAt(new TClonesArray("AliMFTCluster"), iPlane);
739 //====================================================================================================================================================