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