MFT/AliMFT.cxx

   1 // **************************************************************************
   2 // * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 // *                                                                        *
   4 // * Author: The ALICE Off-line Project.                                    *
   5 // * Contributors are mentioned in the code where appropriate.              *
   6 // *                                                                        *
   7 // * Permission to use, copy, modify and distribute this software and its   *
   8 // * documentation strictly for non-commercial purposes is hereby granted   *
   9 // * without fee, provided that the above copyright notice appears in all   *
  10 // * copies and that both the copyright notice and this permission notice   *
  11 // * appear in the supporting documentation. The authors make no claims     *
  12 // * about the suitability of this software for any purpose. It is          *
  13 // * provided "as is" without express or implied warranty.                  *
  14 // **************************************************************************
  15
  16 //====================================================================================================================================================
  17 //
  18 //      Main class of the ALICE Muon Forward Tracker
  19 //
  20 //      Contact author: antonio.uras@cern.ch
  21 //
  22 //====================================================================================================================================================
  23
  24 #include "AliLog.h"
  25 #include "TFile.h"
  26 #include "TGeoManager.h"
  27 #include "TGeoVolume.h"
  28 #include "TGeoMatrix.h"
  29 #include "TVirtualMC.h"
  30 #include "TClonesArray.h"
  31 #include "TGeoGlobalMagField.h"
  32 #include "AliRun.h"
  33 #include "AliLoader.h"
  34 #include "AliDetector.h"
  35 #include "AliMC.h"
  36 #include "AliMagF.h"
  37 #include "AliMFT.h"
  38 #include "AliMFTHit.h"
  39 #include "AliMFTDigit.h"
  40 #include "AliMFTCluster.h"
  41 #include "AliTrackReference.h"
  42 #include "AliMFTSegmentation.h"
  43 #include "AliMFTDigitizer.h"
  44 #include "AliMFTPlane.h"
  45 #include "TString.h"
  46 #include "TObjArray.h"
  47
  48 ClassImp(AliMFT)
  49
  50 //====================================================================================================================================================
  51
  52 AliMFT::AliMFT():
  53   AliDetector(),
  54   fVersion(1),
  55   fNPlanes(0),
  56   fNSlices(1),
  57   fSDigitsPerPlane(0),
  58   fDigitsPerPlane(0),
  59   fRecPointsPerPlane(0),
  60   fSideDigits(0),
  61   fSegmentation(0),
  62   fNameGeomFile(0),
  63   fChargeDispersion(25.e-4),
  64   fSingleStepForChargeDispersion(0),
  65   fNStepForChargeDispersion(4),
  66   fDensitySupportOverSi(0.15)
  67 {
  68
  69   // default constructor
  70
  71 }
  72
  73 //====================================================================================================================================================
  74
  75 AliMFT::AliMFT(const Char_t *name, const Char_t *title):
  76   AliDetector(name, title),
  77   fVersion(1),
  78   fNPlanes(0),
  79   fNSlices(1),
  80   fSDigitsPerPlane(0),
  81   fDigitsPerPlane(0),
  82   fRecPointsPerPlane(0),
  83   fSideDigits(0),
  84   fSegmentation(0),
  85   fNameGeomFile(0),
  86   fChargeDispersion(25.e-4),
  87   fSingleStepForChargeDispersion(0),
  88   fNStepForChargeDispersion(4),
  89   fDensitySupportOverSi(0.15)
  90 {
  91
  92   fNameGeomFile = "AliMFTGeometry.root";
  93
  94   SetGeometry();
  95
  96   Init();
  97
  98 }
  99
 100 //====================================================================================================================================================
 101
 102 AliMFT::AliMFT(const Char_t *name, const Char_t *title, Char_t *nameGeomFile):
 103   AliDetector(name, title),
 104   fVersion(1),
 105   fNPlanes(0),
 106   fNSlices(1),
 107   fSDigitsPerPlane(0),
 108   fDigitsPerPlane(0),
 109   fRecPointsPerPlane(0),
 110   fSideDigits(0),
 111   fSegmentation(0),
 112   fNameGeomFile(0),
 113   fChargeDispersion(25.e-4),
 114   fSingleStepForChargeDispersion(0),
 115   fNStepForChargeDispersion(4),
 116   fDensitySupportOverSi(0.15)
 117 {
 118
 119   fNameGeomFile = nameGeomFile;
 120
 121   SetGeometry();
 122
 123   Init();
 124
 125 }
 126
 127 //====================================================================================================================================================
 128
 129 AliMFT::~AliMFT() {
 130
 131   if (fSDigitsPerPlane)   { fSDigitsPerPlane->Delete();    delete fSDigitsPerPlane;   }
 132   if (fDigitsPerPlane)    { fDigitsPerPlane->Delete();     delete fDigitsPerPlane;    }
 133   if (fRecPointsPerPlane) { fRecPointsPerPlane->Delete();  delete fRecPointsPerPlane; }
 134
 135 }
 136
 137 //====================================================================================================================================================
 138
 139 void AliMFT::CreateMaterials() {
 140
 141   // Definition of MFT materials  - to be updated to the most recent values
 142
 143   AliInfo("Start MFT materials");
 144
 145   // data from PDG booklet 2002     density [gr/cm^3] rad len [cm] abs len [cm]
 146   Float_t   aAir[4]={12,14,16,36} ,   zAir[4]={6,7,8,18} ,   wAir[4]={0.000124,0.755267,0.231781,0.012827} , dAir=0.00120479; Int_t nAir=4;   // Air mixture
 147   Float_t   aSi = 28.085 ,            zSi   = 14 ,           dSi   =  2.329    ,   radSi   =  21.82/dSi , absSi   = 108.4/dSi  ;              // Silicon
 148
 149   Int_t   matId  = 0;                        // tmp material id number
 150   Int_t   unsens = 0, sens=1;                // sensitive or unsensitive medium
 151   Int_t   itgfld = 3;                        // type of field intergration 0 no field -1 user in guswim 1 Runge Kutta 2 helix 3 const field along z
 152   Float_t maxfld = 5.;                       // max field value
 153
 154   Float_t tmaxfd = -10.0;                    // max deflection angle due to magnetic field in one step
 155   Float_t stemax =  0.001;                   // max step allowed [cm]
 156   Float_t deemax = -0.2;                     // maximum fractional energy loss in one step 0<deemax<=1
 157   Float_t epsil  =  0.001;                   // tracking precision [cm]
 158   Float_t stmin  = -0.001;                   // minimum step due to continuous processes [cm] (negative value: choose it automatically)
 159
 160   Float_t tmaxfdSi =  0.1;                    // max deflection angle due to magnetic field in one step
 161   Float_t stemaxSi =  5.0e-4;                 // maximum step allowed [cm]
 162   Float_t deemaxSi =  0.1;                    // maximum fractional energy loss in one step 0<deemax<=1
 163   Float_t epsilSi  =  0.5e-4;                 // tracking precision [cm]
 164   Float_t stminSi  = -0.001;                 // minimum step due to continuous processes [cm] (negative value: choose it automatically)
 165
 166   Int_t   isxfld = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Integ();   // from CreateMaterials in STRUCT/AliPIPEv3.cxx
 167   Float_t sxmgmx = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Max();     // from CreateMaterials in STRUCT/AliPIPEv3.cxx
 168
 169   AliMixture(++matId,"Air", aAir,  zAir,   dAir,   nAir,   wAir);
 170   AliMedium(kAir,    "Air", matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
 171
 172   AliMaterial(++matId, "Si", aSi,   zSi,  dSi,    radSi,  absSi  );
 173   AliMedium(kSi,       "Si", matId, sens, isxfld, sxmgmx, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
 174
 175   AliMaterial(++matId, "Readout", aSi,   zSi,    dSi,    radSi,  absSi  );
 176   AliMedium(kReadout,  "Readout", matId, unsens, isxfld, sxmgmx, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
 177
 178   AliMaterial(++matId, "Support", aSi,   zSi,    dSi*fDensitySupportOverSi, radSi/fDensitySupportOverSi, absSi/fDensitySupportOverSi);
 179   AliMedium(kSupport,  "Support", matId, unsens, isxfld,  sxmgmx,    tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
 180
 181   AliInfo("End MFT materials");
 182
 183 }
 184
 185 //====================================================================================================================================================
 186
 187 void AliMFT::CreateGeometry() {
 188
 189   // Creates detailed geometry simulation (currently GEANT volumes tree)
 190
 191   AliInfo("Start MFT preliminary version building");
 192   if(!gMC->IsRootGeometrySupported()) return;
 193   TGeoVolumeAssembly *vol = CreateVol();
 194   AliInfo("TGeoVolumeAssembly created!");
 195   gGeoManager->GetVolume("ALIC")->AddNode(vol,0);
 196   AliInfo("Stop MFT preliminary version building");
 197
 198   if (fNStepForChargeDispersion) fSingleStepForChargeDispersion = fChargeDispersion/Double_t(fNStepForChargeDispersion);
 199
 200 }
 201
 202 //====================================================================================================================================================
 203
 204 void AliMFT::AddAlignableVolumes() {
 205
 206   // Create entries for alignable volumes associating the symbolic volume
 207   // name with the corresponding volume path. Needs to be syncronized with
 208   // eventual changes in the geometry.
 209
 210   TString sysName = "MFT";
 211   TString volPath = "/ALIC_1/MFT_0";
 212
 213   if (!gGeoManager->SetAlignableEntry(sysName.Data(),volPath.Data())) {
 214     AliFatal(Form("Alignable entry %s not created. Volume path %s not valid", sysName.Data(), volPath.Data()));
 215   }
 216
 217 }
 218
 219 //====================================================================================================================================================
 220
 221 void AliMFT::StepManager() {
 222
 223   // Full Step Manager
 224
 225   AliDebug(2, Form("Entering StepManager: gMC->CurrentVolName() = %s", gMC->CurrentVolName()));
 226
 227   if (!fSegmentation) AliFatal("No segmentation available");    // DO WE HAVE A SEGMENTATION???
 228
 229   if (!(this->IsActive())) return;
 230   if (!(gMC->TrackCharge())) return;
 231
 232   TString planeNumber   = gMC->CurrentVolName();
 233   TString detElemNumber = gMC->CurrentVolName();
 234   if (planeNumber.Contains("support")) return;
 235   if (planeNumber.Contains("readout")) return;
 236   planeNumber.Remove(0,9);
 237   detElemNumber.Remove(0,19);
 238   planeNumber.Remove(2);
 239   detElemNumber.Remove(3);
 240   Int_t detElemID = fSegmentation->GetDetElemID(planeNumber.Atoi(), detElemNumber.Atoi());
 241
 242   if (gMC->IsTrackExiting()) {
 243     AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kMFT);
 244   }
 245
 246   static TLorentzVector position, momentum;
 247   static AliMFTHit hit;
 248
 249   Int_t  status = 0;
 250
 251   // Track status
 252   if (gMC->IsTrackInside())      status +=  1;
 253   if (gMC->IsTrackEntering())    status +=  2;
 254   if (gMC->IsTrackExiting())     status +=  4;
 255   if (gMC->IsTrackOut())         status +=  8;
 256   if (gMC->IsTrackDisappeared()) status += 16;
 257   if (gMC->IsTrackStop())        status += 32;
 258   if (gMC->IsTrackAlive())       status += 64;
 259
 260   // ---------- Fill hit structure
 261
 262   hit.SetDetElemID(detElemID);
 263   hit.SetPlane(planeNumber.Atoi());
 264   hit.SetTrack(gAlice->GetMCApp()->GetCurrentTrackNumber());
 265
 266   gMC->TrackPosition(position);
 267   gMC->TrackMomentum(momentum);
 268
 269   AliDebug(1, Form("AliMFT::StepManager()->%s Hit #%06d (x=%f, y=%f, z=%f) belongs to track %02d\n",
 270                    gMC->CurrentVolName(), fNhits, position.X(), position.Y(), position.Z(), gAlice->GetMCApp()->GetCurrentTrackNumber()));
 271
 272   hit.SetPosition(position);
 273   hit.SetTOF(gMC->TrackTime());
 274   hit.SetMomentum(momentum);
 275   hit.SetStatus(status);
 276   hit.SetEloss(gMC->Edep());
 277   //  hit.SetShunt(GetIshunt());
 278 //   if (gMC->IsTrackEntering()) {
 279 //     hit.SetStartPosition(position);
 280 //     hit.SetStartTime(gMC->TrackTime());
 281 //     hit.SetStartStatus(status);
 282 //     return; // don't save entering hit.
 283 //   }
 284
 285   // Fill hit structure with this new hit.
 286   new ((*fHits)[fNhits++]) AliMFTHit(hit);
 287
 288   // Save old position... for next hit.
 289 //   hit.SetStartPosition(position);
 290 //   hit.SetStartTime(gMC->TrackTime());
 291 //   hit.SetStartStatus(status);
 292
 293   return;
 294
 295 }
 296
 297 //====================================================================================================================================================
 298
 299 TGeoVolumeAssembly* AliMFT::CreateVol() {
 300
 301   // method to create the MFT Geometry (silicon circular planes)
 302
 303   if (!fSegmentation) CreateGeometry();
 304
 305   TGeoVolumeAssembly *vol = new TGeoVolumeAssembly("MFT");
 306   TGeoMedium *silicon = gGeoManager->GetMedium("MFT_Si");
 307   TGeoMedium *readout = gGeoManager->GetMedium("MFT_Readout");
 308   TGeoMedium *support = gGeoManager->GetMedium("MFT_Support");
 309   for (Int_t iPar=0; iPar<8; iPar++) AliInfo(Form("silicon->GetParam(%d) = %f", iPar, silicon->GetParam(iPar)));
 310   for (Int_t iPar=0; iPar<8; iPar++) AliInfo(Form("readout->GetParam(%d) = %f", iPar, readout->GetParam(iPar)));
 311   for (Int_t iPar=0; iPar<8; iPar++) AliInfo(Form("support->GetParam(%d) = %f", iPar, support->GetParam(iPar)));
 312
 313   Double_t origin[3] = {0};
 314
 315   for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
 316
 317     AliDebug(1, Form("Creating volumes for MFT plane %02d",iPlane));
 318
 319     AliMFTPlane *plane = fSegmentation->GetPlane(iPlane);
 320
 321     // --------- support element(s)
 322
 323     origin[0] =  0.5*(plane->GetSupportElement(0)->GetAxis(0)->GetXmax() + plane->GetSupportElement(0)->GetAxis(0)->GetXmin());
 324     origin[1] =  0.5*(plane->GetSupportElement(0)->GetAxis(1)->GetXmax() + plane->GetSupportElement(0)->GetAxis(1)->GetXmin());
 325     origin[2] = -0.5*(plane->GetSupportElement(0)->GetAxis(2)->GetXmax() + plane->GetSupportElement(0)->GetAxis(2)->GetXmin());
 326     TGeoVolume *supportElem = gGeoManager->MakeTube(Form("MFT_plane%02d_support", iPlane), support,
 327                                                     plane->GetRMinSupport(),
 328                                                     plane->GetRMaxSupport(),
 329                                                     0.5*(plane->GetSupportElement(0)->GetAxis(2)->GetXmax() -
 330                                                          plane->GetSupportElement(0)->GetAxis(2)->GetXmin()) );
 331     AliDebug(2, Form("Created vol %s", supportElem->GetName()));
 332     supportElem->SetLineColor(kCyan-9);
 333     vol -> AddNode(supportElem, 0, new TGeoTranslation(origin[0], origin[1], origin[2]));
 334
 335     AliDebug(1, "support elements created!");
 336
 337     // --------- active elements
 338
 339     for (Int_t iActive=0; iActive<plane->GetNActiveElements(); iActive++) {
 340
 341       Double_t dx = 0.5*TMath::Abs(plane->GetActiveElement(iActive)->GetAxis(0)->GetXmax() - plane->GetActiveElement(iActive)->GetAxis(0)->GetXmin());
 342       Double_t dy = 0.5*TMath::Abs(plane->GetActiveElement(iActive)->GetAxis(1)->GetXmax() - plane->GetActiveElement(iActive)->GetAxis(1)->GetXmin());
 343       Double_t dz = 0.5*TMath::Abs(plane->GetActiveElement(iActive)->GetAxis(2)->GetXmax() - plane->GetActiveElement(iActive)->GetAxis(2)->GetXmin());
 344       dz /= Double_t(fNSlices);
 345
 346       origin[0] =  0.5*(plane->GetActiveElement(iActive)->GetAxis(0)->GetXmax() + plane->GetActiveElement(iActive)->GetAxis(0)->GetXmin());
 347       origin[1] =  0.5*(plane->GetActiveElement(iActive)->GetAxis(1)->GetXmax() + plane->GetActiveElement(iActive)->GetAxis(1)->GetXmin());
 348
 349       for (Int_t iSlice=0; iSlice<fNSlices; iSlice++) {
 350         origin[2] = -0.5*(plane->GetActiveElement(iActive)->GetAxis(2)->GetXmin() + 2*dz*(iSlice+1) + plane->GetActiveElement(iActive)->GetAxis(2)->GetXmin() + 2*dz*(iSlice) );
 351         TGeoVolume *activeElem = gGeoManager->MakeBox(Form("MFT_plane%02d_active%03d_slice%02d", iPlane, iActive, iSlice), silicon, dx, dy, dz);
 352         AliDebug(2, Form("Created vol %s", activeElem->GetName()));
 353         activeElem->SetLineColor(kGreen);
 354         vol -> AddNode(activeElem, 0, new TGeoTranslation(origin[0], origin[1], origin[2]));
 355       }
 356
 357     }
 358
 359     AliDebug(1, "active elements created!");
 360
 361     // --------- readout elements
 362
 363     for (Int_t iReadout=0; iReadout<plane->GetNReadoutElements(); iReadout++) {
 364
 365       Double_t dx = 0.5*TMath::Abs(plane->GetReadoutElement(iReadout)->GetAxis(0)->GetXmax() - plane->GetReadoutElement(iReadout)->GetAxis(0)->GetXmin());
 366       Double_t dy = 0.5*TMath::Abs(plane->GetReadoutElement(iReadout)->GetAxis(1)->GetXmax() - plane->GetReadoutElement(iReadout)->GetAxis(1)->GetXmin());
 367       Double_t dz = 0.5*TMath::Abs(plane->GetReadoutElement(iReadout)->GetAxis(2)->GetXmax() - plane->GetReadoutElement(iReadout)->GetAxis(2)->GetXmin());
 368
 369       origin[0] =  0.5*(plane->GetReadoutElement(iReadout)->GetAxis(0)->GetXmax() + plane->GetReadoutElement(iReadout)->GetAxis(0)->GetXmin());
 370       origin[1] =  0.5*(plane->GetReadoutElement(iReadout)->GetAxis(1)->GetXmax() + plane->GetReadoutElement(iReadout)->GetAxis(1)->GetXmin());
 371       origin[2] = -0.5*(plane->GetReadoutElement(iReadout)->GetAxis(2)->GetXmax() + plane->GetReadoutElement(iReadout)->GetAxis(2)->GetXmin());
 372
 373       TGeoVolume *readoutElem = gGeoManager->MakeBox(Form("MFT_plane%02d_readout%03d", iPlane, iReadout), readout, dx, dy, dz);
 374       AliDebug(2, Form("Created vol %s", readoutElem->GetName()));
 375       readoutElem->SetLineColor(kRed);
 376       vol -> AddNode(readoutElem, 0, new TGeoTranslation(origin[0], origin[1], origin[2]));
 377
 378     }
 379
 380     AliDebug(1, "readout elements created!");
 381
 382   }
 383
 384   return vol;
 385
 386 }
 387
 388 //====================================================================================================================================================
 389
 390 void AliMFT::Hits2SDigits(){
 391
 392   // Interface method invoked from AliSimulation to create a list of sdigits corresponding to list of hits. Every hit generates one sdigit.
 393
 394   AliDebug(1,"Start Hits2SDigits.");
 395
 396   if (!fSegmentation) CreateGeometry();
 397
 398   if (!fLoader->TreeH()) fLoader->LoadHits();
 399
 400   if (!fLoader->TreeS()) {
 401
 402     for (Int_t iEvt=0;iEvt<fLoader->GetRunLoader()->GetNumberOfEvents(); iEvt++) {
 403
 404       fLoader->GetRunLoader()->GetEvent(iEvt);
 405       fLoader->MakeTree("S");
 406       MakeBranch("S");
 407       SetTreeAddress();
 408
 409       AliDebug(1, Form("Event %03d: fLoader->TreeH()->GetEntries() = %2d", iEvt, Int_t(fLoader->TreeH()->GetEntries())));
 410
 411       for (Int_t iTrack=0; iTrack<fLoader->TreeH()->GetEntries(); iTrack++) {
 412         fLoader->TreeH()->GetEntry(iTrack);
 413         Hits2SDigitsLocal(Hits(), GetSDigitsList(), iTrack);    // convert these hits to a list of sdigits
 414       }
 415
 416       fLoader->TreeS()->Fill();
 417       fLoader->WriteSDigits("OVERWRITE");
 418       ResetSDigits();
 419
 420     }
 421   }
 422
 423   fLoader->UnloadHits();
 424   fLoader->UnloadSDigits();
 425
 426   AliDebug(1,"Stop Hits2SDigits.");
 427
 428 }
 429
 430 //====================================================================================================================================================
 431
 432 void AliMFT::Hits2SDigitsLocal(TClonesArray *hits, const TObjArray *pSDig, Int_t track) {
 433
 434   //  Add sdigits of these hits to the list
 435
 436   AliDebug(1, "Entering Hits2SDigitsLocal");
 437
 438   if (!fSegmentation) CreateGeometry();
 439
 440   TClonesArray *pSDigList[fNMaxPlanes];
 441   for (Int_t iPlane=0; iPlane<fNMaxPlanes; iPlane++) pSDigList[iPlane] = NULL;
 442   for (Int_t iPlane=0; iPlane<fNPlanes;    iPlane++) {
 443     pSDigList[iPlane] = (TClonesArray*) (*pSDig)[iPlane];
 444     AliDebug(1,Form("Entries of pSDigList %3d; plane: %02d,",pSDigList[iPlane]->GetEntries(),iPlane));
 445     if (!track && pSDigList[iPlane]->GetEntries()!=0) AliErrorClass("Some of sdigits lists is not empty");
 446   }
 447
 448   for (Int_t iHit=0; iHit<hits->GetEntries(); iHit++) {
 449
 450     AliMFTHit *hit = (AliMFTHit*) hits->At(iHit);
 451
 452     // Creating "main digit"
 453
 454     AliMFTDigit *mainSDigit = new AliMFTDigit();
 455     mainSDigit->SetEloss(hit->GetEloss());
 456     mainSDigit->SetDetElemID(hit->GetDetElemID());
 457     mainSDigit->SetPlane(hit->GetPlane());
 458     mainSDigit->AddMCLabel(hit->GetTrack());
 459
 460     Int_t xPixel = -1;
 461     Int_t yPixel = -1;
 462     if (fSegmentation->Hit2PixelID(hit->X(), hit->Y(), mainSDigit->GetDetElemID(), xPixel, yPixel)) {
 463       mainSDigit->SetPixID(xPixel, yPixel, 0);
 464       mainSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(mainSDigit->GetDetElemID()),
 465                           fSegmentation->GetPixelSizeY(mainSDigit->GetDetElemID()),
 466                           fSegmentation->GetPixelSizeZ(mainSDigit->GetDetElemID()));
 467       mainSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(mainSDigit->GetDetElemID(), xPixel),
 468                            fSegmentation->GetPixelCenterY(mainSDigit->GetDetElemID(), yPixel),
 469                            fSegmentation->GetPixelCenterZ(mainSDigit->GetDetElemID(), 0));
 470       new ((*fSideDigits)[fSideDigits->GetEntries()]) AliMFTDigit(*mainSDigit);
 471       AliDebug(1, Form("Created new sdigit (%f, %f, %f) from hit (%f, %f, %f)",
 472                        mainSDigit->GetPixelCenterX(), mainSDigit->GetPixelCenterY(), mainSDigit->GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
 473     }
 474
 475     // creating "side digits" to simulate the effect of charge dispersion
 476
 477     Double_t pi4 = TMath::Pi()/4.;
 478     for (Int_t iStep=0; iStep<fNStepForChargeDispersion; iStep++) {
 479       Double_t shift = (iStep+1) * fSingleStepForChargeDispersion;
 480       for (Int_t iAngle=0; iAngle<8; iAngle++) {
 481         Double_t shiftX = shift*TMath::Cos(iAngle*pi4);
 482         Double_t shiftY = shift*TMath::Sin(iAngle*pi4);
 483         if (fSegmentation->Hit2PixelID(hit->X()+shiftX, hit->Y()+shiftY, hit->GetDetElemID(), xPixel, yPixel)) {
 484           Bool_t digitExists = kFALSE;
 485           for (Int_t iSideDigit=0; iSideDigit<fSideDigits->GetEntries(); iSideDigit++) {
 486             if (xPixel==((AliMFTDigit*) fSideDigits->At(iSideDigit))->GetPixelX() &&
 487                 yPixel==((AliMFTDigit*) fSideDigits->At(iSideDigit))->GetPixelY()) {
 488               digitExists = kTRUE;
 489               break;
 490             }
 491           }
 492           if (!digitExists) {
 493             AliMFTDigit *sideSDigit = new AliMFTDigit();
 494             sideSDigit->SetEloss(0.);
 495             sideSDigit->SetDetElemID(hit->GetDetElemID());
 496             sideSDigit->SetPlane(hit->GetPlane());
 497             sideSDigit->AddMCLabel(hit->GetTrack());
 498             sideSDigit->SetPixID(xPixel, yPixel, 0);
 499             sideSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(sideSDigit->GetDetElemID()),
 500                                     fSegmentation->GetPixelSizeY(sideSDigit->GetDetElemID()),
 501                                     fSegmentation->GetPixelSizeZ(sideSDigit->GetDetElemID()));
 502             sideSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(sideSDigit->GetDetElemID(), xPixel),
 503                                      fSegmentation->GetPixelCenterY(sideSDigit->GetDetElemID(), yPixel),
 504                                      fSegmentation->GetPixelCenterZ(sideSDigit->GetDetElemID(), 0));
 505             new ((*fSideDigits)[fSideDigits->GetEntries()]) AliMFTDigit(*sideSDigit);
 506           }
 507         }
 508       }
 509     }
 510
 511     // ------------ In case we should simulate a rectangular pattern of pixel...
 512
 513     if (fSegmentation->GetPlane(mainSDigit->GetPlane())->HasPixelRectangularPatternAlongY()) {
 514       for (Int_t iSDigit=0; iSDigit<fSideDigits->GetEntries(); iSDigit++) {
 515         AliMFTDigit *mySDig = (AliMFTDigit*) (fSideDigits->At(iSDigit));
 516         if (mySDig->GetPixelX()%2 == mySDig->GetPixelY()%2) {   // both pair or both odd
 517           xPixel = mySDig->GetPixelX();
 518           yPixel = mySDig->GetPixelY()+1;
 519           if (fSegmentation->DoesPixelExist(mySDig->GetDetElemID(), xPixel, yPixel)) {
 520             AliMFTDigit *newSDigit = new AliMFTDigit();
 521             newSDigit->SetEloss(0.);
 522             newSDigit->SetDetElemID(mySDig->GetDetElemID());
 523             newSDigit->SetPlane(mySDig->GetDetElemID());
 524             newSDigit->SetPixID(xPixel, yPixel, 0);
 525             newSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(newSDigit->GetDetElemID()),
 526                                    fSegmentation->GetPixelSizeY(newSDigit->GetDetElemID()),
 527                                    fSegmentation->GetPixelSizeZ(newSDigit->GetDetElemID()));
 528             newSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(newSDigit->GetDetElemID(), xPixel),
 529                                     fSegmentation->GetPixelCenterY(newSDigit->GetDetElemID(), yPixel),
 530                                     fSegmentation->GetPixelCenterZ(newSDigit->GetDetElemID(), 0));
 531             new ((*fSideDigits)[fSideDigits->GetEntries()]) AliMFTDigit(*newSDigit);
 532           }
 533         }
 534         else {   // pair-odd
 535           xPixel = mySDig->GetPixelX();
 536           yPixel = mySDig->GetPixelY()-1;
 537           if (fSegmentation->DoesPixelExist(mySDig->GetDetElemID(), xPixel, yPixel)) {
 538             AliMFTDigit *newSDigit = new AliMFTDigit();
 539             newSDigit->SetEloss(0.);
 540             newSDigit->SetDetElemID(mySDig->GetDetElemID());
 541             newSDigit->SetPlane(mySDig->GetPlane());
 542             newSDigit->SetPixID(xPixel, yPixel, 0);
 543             newSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(newSDigit->GetDetElemID()),
 544                                    fSegmentation->GetPixelSizeY(newSDigit->GetDetElemID()),
 545                                    fSegmentation->GetPixelSizeZ(newSDigit->GetDetElemID()));
 546             newSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(newSDigit->GetDetElemID(), xPixel),
 547                                     fSegmentation->GetPixelCenterY(newSDigit->GetDetElemID(), yPixel),
 548                                     fSegmentation->GetPixelCenterZ(newSDigit->GetDetElemID(), 0));
 549             new ((*fSideDigits)[fSideDigits->GetEntries()]) AliMFTDigit(*newSDigit);
 550           }
 551         }
 552       }
 553     }
 554
 555     // -------- checking which pixels switched on have their diode actually within the charge dispersion radius
 556
 557     for (Int_t iSDigit=0; iSDigit<fSideDigits->GetEntries(); iSDigit++) {
 558       AliMFTDigit *mySDig = (AliMFTDigit*) (fSideDigits->At(iSDigit));
 559       Double_t distance = TMath::Sqrt(TMath::Power(mySDig->GetPixelCenterX()-hit->X(),2) + TMath::Power(mySDig->GetPixelCenterY()-hit->Y(),2));
 560       if (fSegmentation->GetPlane(mySDig->GetPlane())->HasPixelRectangularPatternAlongY()) {
 561         if (mySDig->GetPixelX()%2 == mySDig->GetPixelY()%2) {  // both pair or both odd
 562           if (distance<fChargeDispersion) {
 563             AliDebug(1, Form("Created new sdigit (%f, %f, %f) from hit (%f, %f, %f)",
 564                              mySDig->GetPixelCenterX(), mySDig->GetPixelCenterY(), mySDig->GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
 565             new ((*pSDigList[mySDig->GetPlane()])[pSDigList[mySDig->GetPlane()]->GetEntries()]) AliMFTDigit(*mySDig);
 566             xPixel = mySDig->GetPixelX();
 567             yPixel = mySDig->GetPixelY()+1;
 568             if (fSegmentation->DoesPixelExist(mySDig->GetDetElemID(), xPixel, yPixel)) {
 569               AliMFTDigit *newSDigit = new AliMFTDigit();
 570               newSDigit->SetEloss(0.);
 571               newSDigit->SetDetElemID(mySDig->GetDetElemID());
 572               newSDigit->SetPlane(mySDig->GetPlane());
 573               newSDigit->SetPixID(xPixel, yPixel, 0);
 574               newSDigit->SetPixWidth(fSegmentation->GetPixelSizeX(newSDigit->GetDetElemID()),
 575                                      fSegmentation->GetPixelSizeY(newSDigit->GetDetElemID()),
 576                                      fSegmentation->GetPixelSizeZ(newSDigit->GetDetElemID()));
 577               newSDigit->SetPixCenter(fSegmentation->GetPixelCenterX(newSDigit->GetDetElemID(), xPixel),
 578                                       fSegmentation->GetPixelCenterY(newSDigit->GetDetElemID(), yPixel),
 579                                       fSegmentation->GetPixelCenterZ(newSDigit->GetDetElemID(), 0));
 580               AliDebug(1, Form("Created new sdigit (%f, %f, %f) from hit (%f, %f, %f)",
 581                                newSDigit->GetPixelCenterX(), newSDigit->GetPixelCenterY(), newSDigit->GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
 582               new ((*pSDigList[newSDigit->GetPlane()])[pSDigList[newSDigit->GetPlane()]->GetEntries()]) AliMFTDigit(*newSDigit);
 583             }
 584           }
 585         }
 586       }
 587       else {
 588         if (distance<fChargeDispersion) {
 589           AliDebug(1, Form("Created new sdigit (%f, %f, %f) from hit (%f, %f, %f)",
 590                            mySDig->GetPixelCenterX(), mySDig->GetPixelCenterY(), mySDig->GetPixelCenterZ(), hit->X(), hit->Y(), hit->Z()));
 591           new ((*pSDigList[mySDig->GetPlane()])[pSDigList[mySDig->GetPlane()]->GetEntries()]) AliMFTDigit(*mySDig);
 592         }
 593       }
 594     }
 595
 596     fSideDigits->Delete();
 597
 598   }
 599
 600   AliDebug(1,"Exiting Hits2SDigitsLocal");
 601
 602 }
 603
 604 //====================================================================================================================================================
 605
 606 void AliMFT::MakeBranch(Option_t *option) {
 607
 608   // Create Tree branches
 609   AliDebug(1, Form("Start with option= %s.",option));
 610
 611   const Int_t kBufSize = 4000;
 612
 613   const Char_t *cH = strstr(option,"H");
 614   const Char_t *cD = strstr(option,"D");
 615   const Char_t *cS = strstr(option,"S");
 616
 617   if (cH && fLoader->TreeH()) {
 618     CreateHits();
 619     MakeBranchInTree(fLoader->TreeH(), "MFT", &fHits, kBufSize, 0);
 620   }
 621
 622   if (cS && fLoader->TreeS()) {
 623     CreateSDigits();
 624     for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) MakeBranchInTree(fLoader->TreeS(),
 625                                                                     Form("Plane_%02d",iPlane),
 626                                                                     &((*fSDigitsPerPlane)[iPlane]),
 627                                                                     kBufSize, 0);
 628   }
 629
 630   if (cD && fLoader->TreeD()) {
 631     CreateDigits();
 632     for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) MakeBranchInTree(fLoader->TreeD(),
 633                                                                     Form("Plane_%02d",iPlane),
 634                                                                     &((*fDigitsPerPlane)[iPlane]),
 635                                                                     kBufSize, 0);
 636   }
 637
 638   AliDebug(1,"Stop.");
 639
 640 }
 641
 642 //====================================================================================================================================================
 643
 644 void AliMFT::SetTreeAddress() {
 645
 646   AliDebug(1, "AliMFT::SetTreeAddress()");
 647
 648   //Set branch address for the Hits and Digits Tree.
 649   AliDebug(1, "Start.");
 650
 651   AliDebug(1, Form("AliMFT::SetTreeAddress Hits  fLoader->TreeH() = %p\n", fLoader->TreeH()));
 652   if (fLoader->TreeH() && fLoader->TreeH()->GetBranch("MFT")) {
 653     CreateHits();
 654     fLoader->TreeH()->SetBranchAddress("MFT", &fHits);
 655   }
 656
 657   AliDebug(1, Form("AliMFT::SetTreeAddress SDigits  fLoader->TreeS() = %p\n", fLoader->TreeS()));
 658   if (fLoader->TreeS() && fLoader->TreeS()->GetBranch("Plane_00")) {
 659     CreateSDigits();
 660     for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
 661       fLoader->TreeS()->SetBranchAddress(Form("Plane_%02d",iPlane), &((*fSDigitsPerPlane)[iPlane]));
 662     }
 663   }
 664
 665   AliDebug(1, Form("AliMFT::SetTreeAddress Digits  fLoader->TreeD() = %p\n", fLoader->TreeD()));
 666   if (fLoader->TreeD() && fLoader->TreeD()->GetBranch("Plane_00")) {
 667     CreateDigits();
 668     for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
 669       fLoader->TreeD()->SetBranchAddress(Form("Plane_%02d",iPlane), &((*fDigitsPerPlane)[iPlane]));
 670     }
 671   }
 672
 673   AliDebug(1, Form("AliMFT::SetTreeAddress RecPoints  fLoader->TreeR() = %p\n", fLoader->TreeR()));
 674   if (fLoader->TreeR() && fLoader->TreeR()->GetBranch("Plane_00")) {
 675     CreateRecPoints();
 676     for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
 677       fLoader->TreeR()->SetBranchAddress(Form("Plane_%02d",iPlane), &((*fRecPointsPerPlane)[iPlane]));
 678     }
 679   }
 680
 681   AliDebug(1,"Stop.");
 682
 683 }
 684
 685 //====================================================================================================================================================
 686
 687 void AliMFT::SetGeometry() {
 688
 689   AliInfo("AliMFT::SetGeometry\n");
 690
 691   fSegmentation = new AliMFTSegmentation(fNameGeomFile.Data());
 692
 693   fNPlanes = fSegmentation->GetNPlanes();
 694
 695 }
 696
 697 //====================================================================================================================================================
 698
 699 void AliMFT::CreateHits() {
 700
 701   // create array of hits
 702
 703   AliDebug(1, "AliMFT::CreateHits()");
 704
 705   if (fHits) return;
 706   fHits = new TClonesArray("AliMFTHit");
 707
 708 }
 709
 710 //====================================================================================================================================================
 711
 712 void AliMFT::CreateSDigits() {
 713
 714   // create sdigits list
 715
 716   AliDebug(1, "AliMFT::CreateSDigits()");
 717
 718   if (fSDigitsPerPlane) return;
 719   fSDigitsPerPlane = new TObjArray(fNPlanes);
 720   for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) fSDigitsPerPlane->AddAt(new TClonesArray("AliMFTDigit"), iPlane);
 721
 722   fSideDigits = new TClonesArray("AliMFTDigit");
 723
 724 }
 725
 726 //====================================================================================================================================================
 727
 728 void AliMFT::CreateDigits() {
 729
 730   // create digits list
 731
 732   AliDebug(1, "AliMFT::CreateDigits()");
 733
 734   if (fDigitsPerPlane) return;
 735   fDigitsPerPlane = new TObjArray(fNPlanes);
 736   for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) fDigitsPerPlane->AddAt(new TClonesArray("AliMFTDigit"), iPlane);
 737
 738 }
 739
 740 //====================================================================================================================================================
 741
 742 void AliMFT::CreateRecPoints() {
 743
 744   // create recPoints list
 745
 746   AliDebug(1, "AliMFT::CreateRecPoints()");
 747
 748   if (fRecPointsPerPlane) return;
 749   fRecPointsPerPlane = new TObjArray(fNPlanes);
 750   for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) fRecPointsPerPlane->AddAt(new TClonesArray("AliMFTCluster"), iPlane);
 751
 752 }
 753
 754 //====================================================================================================================================================