MFT/AliMFT.cxx

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