PMD/AliPMDDigitizer.cxx

   1 //-----------------------------------------------------//
   2 //                                                     //
   3 //  Source File : PMDDigitization.cxx, Version 00      //
   4 //                                                     //
   5 //  Date   : September 20 2002                         //
   6 //                                                     //
   7 //-----------------------------------------------------//
   8
   9 #include <TBRIK.h>
  10 #include <TNode.h>
  11 #include <TTree.h>
  12 #include <TGeometry.h>
  13 #include <TObjArray.h>
  14 #include <TClonesArray.h>
  15 #include <TFile.h>
  16 #include <TNtuple.h>
  17 #include <TParticle.h>
  18
  19 #include "AliRun.h"
  20 #include "AliPMD.h"
  21 #include "AliHit.h"
  22 #include "AliDetector.h"
  23 #include "AliRunLoader.h"
  24 #include "AliLoader.h"
  25 #include "AliConfig.h"
  26 #include "AliMagF.h"
  27 #include "AliRunDigitizer.h"
  28 #include "AliHeader.h"
  29
  30 #include "AliPMDcell.h"
  31 #include "AliPMDsdigit.h"
  32 #include "AliPMDdigit.h"
  33 #include "AliPMDDigitizer.h"
  34 #include "AliPMDClustering.h"
  35 #include "AliPMDContainer.h"
  36 #include "AliPMDrecpoint.h"
  37
  38
  39 ClassImp(AliPMDDigitizer)
  40 //
  41 // Constructor
  42 //
  43 AliPMDDigitizer::AliPMDDigitizer()
  44 {
  45   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
  46   fNsdigit = 0;
  47   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
  48   fNdigit = 0;
  49
  50   for (Int_t i = 0; i < fTotSM; i++)
  51     {
  52       for (Int_t j = 0; j < fNCell; j++)
  53         {
  54           for (Int_t k = 0; k < fNCell; k++)
  55             {
  56               fCPV[i][j][k] = 0.;
  57               fPMD[i][j][k] = 0.;
  58             }
  59         }
  60     }
  61
  62   if (!fCell) fCell = new TObjArray();
  63
  64   fZPos = 361.5; // in units of cm, This is the default position of PMD
  65
  66 }
  67 AliPMDDigitizer::~AliPMDDigitizer()
  68 {
  69   delete fSDigits;
  70   delete fDigits;
  71   delete fCell;
  72 }
  73 //
  74 // Member functions
  75 //
  76 void AliPMDDigitizer::OpengAliceFile(Char_t *file, Option_t *option)
  77 {
  78
  79   fRunLoader = AliRunLoader::Open(file,AliConfig::fgkDefaultEventFolderName,
  80                                   "UPDATE");
  81
  82   if (!fRunLoader)
  83    {
  84      Error("Open","Can not open session for file %s.",file);
  85    }
  86
  87   fRunLoader->LoadgAlice();
  88   fRunLoader->LoadHeader();
  89   fRunLoader->LoadKinematics();
  90
  91   gAlice = fRunLoader->GetAliRun();
  92
  93   if (gAlice)
  94     {
  95       printf("<AliPMDdigitizer::Open> ");
  96       printf("AliRun object found on file.\n");
  97     }
  98   else
  99     {
 100       printf("<AliPMDdigitizer::Open> ");
 101       printf("Could not find AliRun object.\n");
 102     }
 103
 104   PMD  = (AliPMD*)gAlice->GetDetector("PMD");
 105   pmdloader = fRunLoader->GetLoader("PMDLoader");
 106   if (pmdloader == 0x0)
 107     {
 108       cerr<<"Hits2Digits : Can not find PMD or PMDLoader\n";
 109     }
 110
 111   const char *cHS = strstr(option,"HS");
 112   const char *cHD = strstr(option,"HD");
 113   const char *cSD = strstr(option,"SD");
 114
 115   if (cHS)
 116     {
 117       pmdloader->LoadHits("READ");
 118       pmdloader->LoadSDigits("recreate");
 119     }
 120   else if (cHD)
 121     {
 122       pmdloader->LoadHits("READ");
 123       pmdloader->LoadDigits("recreate");
 124     }
 125   else if (cSD)
 126     {
 127       pmdloader->LoadSDigits("READ");
 128       pmdloader->LoadDigits("recreate");
 129     }
 130
 131 }
 132 void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
 133 {
 134   cout << " -------- Beginning of Hits2SDigits ----------- " << endl;
 135
 136   Int_t kPi0 = 111;
 137   Int_t kGamma = 22;
 138   Int_t npmd;
 139   Int_t trackno;
 140   Int_t cellnumber;
 141   Int_t smnumber1;
 142   Int_t smnumber;
 143   Int_t trackpid;
 144   Int_t mtrackno;
 145   Int_t mtrackpid;
 146
 147   Float_t xPos, yPos, zPos;
 148   Int_t xpad, ypad;
 149   Float_t edep;
 150   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
 151
 152
 153   ResetSDigit();
 154
 155   printf("Event Number =  %d \n",ievt);
 156   Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
 157   printf("Number of Particles = %d \n", nparticles);
 158   fRunLoader->GetEvent(ievt);
 159   Particles = gAlice->Particles();
 160   // ------------------------------------------------------- //
 161   // Pointer to specific detector hits.
 162   // Get pointers to Alice detectors and Hits containers
 163
 164   treeH = pmdloader->TreeH();
 165
 166   Int_t ntracks    = (Int_t) treeH->GetEntries();
 167   printf("Number of Tracks in the TreeH = %d \n", ntracks);
 168
 169   treeS = pmdloader->TreeS();
 170   if (treeS == 0x0)
 171     {
 172       pmdloader->MakeTree("S");
 173       treeS = pmdloader->TreeS();
 174     }
 175   Int_t bufsize = 16000;
 176   treeS->Branch("PMDSDigit", &fSDigits, bufsize);
 177
 178   if (PMD) PMDhits   = PMD->Hits();
 179
 180   // Start loop on tracks in the hits containers
 181
 182
 183   for (Int_t track=0; track<ntracks;track++)
 184     {
 185       gAlice->ResetHits();
 186       treeH->GetEvent(track);
 187
 188       if (PMD)
 189         {
 190           npmd = PMDhits->GetEntriesFast();
 191           for (int ipmd = 0; ipmd < npmd; ipmd++)
 192             {
 193               pmdHit = (AliPMDhit*) PMDhits->UncheckedAt(ipmd);
 194               trackno = pmdHit->GetTrack();
 195
 196               //  get kinematics of the particles
 197
 198               particle = gAlice->Particle(trackno);
 199               trackpid  = particle->GetPdgCode();
 200
 201               Int_t igatr = -999;
 202               Int_t ichtr = -999;
 203
 204               Int_t imo;
 205               Int_t igen = 0;
 206               Int_t id_mo = -999;
 207               TParticle*  mparticle = particle;
 208
 209               while((imo = mparticle->GetFirstMother()) >= 0)
 210                 {
 211                   igen++;
 212                   mparticle =  gAlice->Particle(imo);
 213                   id_mo = mparticle->GetPdgCode();
 214
 215                   vx = mparticle->Vx();
 216                   vy = mparticle->Vy();
 217                   vz = mparticle->Vz();
 218
 219                   //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, id_mo, vx, vy, vz);
 220                   //fprintf(ftest1,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, id_mo, vx, vy, vz);
 221
 222                   if (id_mo == kGamma && vx == 0. && vy == 0. && vz == 0.)
 223                     {
 224                       igatr = imo;
 225                     }
 226                   ichtr = imo;
 227                 }
 228
 229               if (id_mo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
 230                 {
 231                   mtrackno = igatr;
 232                   mtrackpid = kGamma;
 233                 }
 234               else
 235                 {
 236                   mtrackno  = ichtr;
 237                   mtrackpid = id_mo;
 238                 }
 239
 240               xPos = pmdHit->X();
 241               yPos = pmdHit->Y();
 242               zPos = pmdHit->Z();
 243               cellnumber = pmdHit->fVolume[1];
 244               smnumber1  = pmdHit->fVolume[4];
 245               edep       = pmdHit->fEnergy;
 246
 247               if (smnumber1 > 3 && smnumber1 <= 6)
 248                 {
 249                   Int_t ny = (cellnumber-1)/48 + 1;
 250                   Int_t nx = cellnumber - (ny-1)*48;
 251                   Int_t nx1 = nx + 24;
 252                   Int_t ny1 = ny + 48;
 253                   Int_t ncell1 = (ny1 - 1)* 72 + nx1;
 254                   cellnumber = ncell1;
 255                 }
 256               smnumber = Convert2RealSMNumber(smnumber1);
 257               ypad = (cellnumber - 1)/fNCell + 1;
 258               xpad = cellnumber - (ypad-1)*fNCell;
 259
 260               //cout << "-zpos = " << -zPos << endl;
 261
 262               if (-zPos < fZPos)
 263                 {
 264                   // CPV
 265                   fDetNo = 1;
 266                 }
 267               else if (-zPos > fZPos)
 268                 {
 269                   // PMD
 270                   fDetNo = 0;
 271                 }
 272
 273               if (fDetNo == 0)
 274                 {
 275                   fPMD[smnumber-1][xpad-1][ypad-1] += edep;
 276                   fPMDCounter[smnumber-1][xpad-1][ypad-1]++;
 277                   Int_t smn = smnumber - 1;
 278                   Int_t ixx = xpad     - 1;
 279                   Int_t iyy = ypad     - 1;
 280
 281                   pmdcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
 282
 283                   fCell->Add(pmdcell);
 284                 }
 285               else if(fDetNo == 1)
 286                 {
 287                   fCPV[smnumber-1][xpad-1][ypad-1] += edep;
 288                   fCPVTrackNo[smnumber-1][xpad-1][ypad-1] = mtrackno;
 289                 }
 290             }
 291         }
 292     } // Track Loop ended
 293
 294
 295   TrackAssignment2Cell();
 296
 297   ResetCell();
 298
 299   Float_t deltaE      = 0.;
 300   Int_t   detno       = 0;
 301   Int_t   count_digit = 0;
 302   Int_t   trno        = -1;
 303   Int_t   cellno      = 0;
 304
 305   for (Int_t idet = 0; idet < 2; idet++)
 306     {
 307       for (Int_t ism = 0; ism < fTotSM; ism++)
 308         {
 309           for (Int_t jrow = 0; jrow < fNCell; jrow++)
 310             {
 311               for (Int_t kcol = 0; kcol < fNCell; kcol++)
 312                 {
 313                   cellno = jrow + kcol*fNCell;
 314                   if (idet == 0)
 315                     {
 316                       deltaE = fPMD[ism][jrow][kcol];
 317                       trno   = fPMDTrackNo[ism][jrow][kcol];
 318                       detno = 0;
 319                     }
 320                   else if (idet == 1)
 321                     {
 322                       deltaE = fCPV[ism][jrow][kcol];
 323                       trno   = fCPVTrackNo[ism][jrow][kcol];
 324                       detno = 1;
 325                     }
 326                   if (deltaE > 0.)
 327                     {
 328                       count_digit++;
 329                       AddSDigit(trno,detno,ism,cellno,deltaE);
 330                     }
 331                 }
 332             }
 333           treeS->Fill();
 334           ResetSDigit();
 335         }
 336     }
 337   pmdloader->WriteSDigits("OVERWRITE");
 338
 339   ResetCellADC();
 340
 341   //  cout << " -------- End of Hits2SDigit ----------- " << endl;
 342 }
 343
 344 void AliPMDDigitizer::Hits2Digits(Int_t ievt)
 345 {
 346   Int_t kPi0 = 111;
 347   Int_t kGamma = 22;
 348   Int_t npmd;
 349   Int_t trackno;
 350   Int_t cellnumber;
 351   Int_t smnumber1;
 352   Int_t smnumber;
 353   Int_t trackpid;
 354   Int_t mtrackno;
 355   Int_t mtrackpid;
 356
 357   Float_t xPos, yPos, zPos;
 358   Int_t xpad, ypad;
 359   Float_t edep;
 360   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
 361
 362
 363   ResetDigit();
 364
 365   printf("Event Number =  %d \n",ievt);
 366
 367   Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
 368   printf("Number of Particles = %d \n", nparticles);
 369   fRunLoader->GetEvent(ievt);
 370   Particles = gAlice->Particles();
 371   // ------------------------------------------------------- //
 372   // Pointer to specific detector hits.
 373   // Get pointers to Alice detectors and Hits containers
 374
 375   PMD  = (AliPMD*)gAlice->GetDetector("PMD");
 376   pmdloader = fRunLoader->GetLoader("PMDLoader");
 377
 378   if (pmdloader == 0x0)
 379     {
 380       cerr<<"Hits2Digits method : Can not find PMD or PMDLoader\n";
 381     }
 382   treeH = pmdloader->TreeH();
 383   Int_t ntracks    = (Int_t) treeH->GetEntries();
 384   printf("Number of Tracks in the TreeH = %d \n", ntracks);
 385   pmdloader->LoadDigits("recreate");
 386   treeD = pmdloader->TreeD();
 387   if (treeD == 0x0)
 388     {
 389       pmdloader->MakeTree("D");
 390       treeD = pmdloader->TreeD();
 391     }
 392   Int_t bufsize = 16000;
 393   treeD->Branch("PMDDigit", &fDigits, bufsize);
 394
 395   if (PMD) PMDhits   = PMD->Hits();
 396
 397   // Start loop on tracks in the hits containers
 398
 399   for (Int_t track=0; track<ntracks;track++)
 400     {
 401       gAlice->ResetHits();
 402       treeH->GetEvent(track);
 403
 404       if (PMD)
 405         {
 406           npmd = PMDhits->GetEntriesFast();
 407           for (int ipmd = 0; ipmd < npmd; ipmd++)
 408             {
 409               pmdHit = (AliPMDhit*) PMDhits->UncheckedAt(ipmd);
 410               trackno = pmdHit->GetTrack();
 411
 412               //  get kinematics of the particles
 413
 414               particle = gAlice->Particle(trackno);
 415               trackpid  = particle->GetPdgCode();
 416
 417               Int_t igatr = -999;
 418               Int_t ichtr = -999;
 419
 420               Int_t imo;
 421               Int_t igen = 0;
 422               Int_t id_mo = -999;
 423               TParticle*  mparticle = particle;
 424
 425               while((imo = mparticle->GetFirstMother()) >= 0)
 426                 {
 427                   igen++;
 428                   mparticle =  gAlice->Particle(imo);
 429                   id_mo = mparticle->GetPdgCode();
 430
 431                   vx = mparticle->Vx();
 432                   vy = mparticle->Vy();
 433                   vz = mparticle->Vz();
 434
 435                   //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, id_mo, vx, vy, vz);
 436                   //fprintf(ftest1,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, id_mo, vx, vy, vz);
 437
 438                   if (id_mo == kGamma && vx == 0. && vy == 0. && vz == 0.)
 439                     {
 440                       igatr = imo;
 441                     }
 442                   ichtr = imo;
 443                 }
 444
 445               if (id_mo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
 446                 {
 447                   mtrackno = igatr;
 448                   mtrackpid = kGamma;
 449                 }
 450               else
 451                 {
 452                   mtrackno  = ichtr;
 453                   mtrackpid = id_mo;
 454                 }
 455
 456               xPos = pmdHit->X();
 457               yPos = pmdHit->Y();
 458               zPos = pmdHit->Z();
 459               cellnumber = pmdHit->fVolume[1];
 460               smnumber1  = pmdHit->fVolume[4];
 461               edep       = pmdHit->fEnergy;
 462
 463               if (smnumber1 > 3 && smnumber1 <= 6)
 464                 {
 465                   Int_t ny = (cellnumber-1)/48 + 1;
 466                   Int_t nx = cellnumber - (ny-1)*48;
 467                   Int_t nx1 = nx + 24;
 468                   Int_t ny1 = ny + 48;
 469                   Int_t ncell1 = (ny1 - 1)* 72 + nx1;
 470                   cellnumber = ncell1;
 471                 }
 472
 473               smnumber = Convert2RealSMNumber(smnumber1);
 474               ypad = (cellnumber - 1)/fNCell + 1;
 475               xpad = cellnumber - (ypad-1)*fNCell;
 476
 477               //cout << "-zpos = " << -zPos << endl;
 478               if (-zPos < fZPos)
 479                 {
 480                   // CPV
 481                   fDetNo = 1;
 482                 }
 483               else if (-zPos > fZPos)
 484                 {
 485                   // PMD
 486                   fDetNo = 0;
 487                 }
 488
 489               if(fDetNo == 1)
 490                 {
 491                   fCPV[smnumber-1][xpad-1][ypad-1] += edep;
 492                 }
 493               else if (fDetNo == 0)
 494                 {
 495                   fPMD[smnumber-1][xpad-1][ypad-1] += edep;
 496                   fPMDCounter[smnumber-1][xpad-1][ypad-1]++;
 497                 }
 498             }
 499         }
 500     } // Track Loop ended
 501
 502   TrackAssignment2Cell();
 503   ResetCell();
 504
 505   Float_t deltaE = 0.;
 506   Int_t detno = 0;
 507   Int_t count_digit = 0;
 508   Int_t trno = 1;
 509   Int_t cellno;
 510
 511   for (Int_t idet = 0; idet < 2; idet++)
 512     {
 513       for (Int_t ism = 0; ism < fTotSM; ism++)
 514         {
 515           for (Int_t jrow = 0; jrow < fNCell; jrow++)
 516             {
 517               for (Int_t kcol = 0; kcol < fNCell; kcol++)
 518                 {
 519                   cellno = jrow + kcol*fNCell;
 520                   if (idet == 0)
 521                     {
 522                       deltaE = fPMD[ism][jrow][kcol];
 523                       detno = 0;
 524                     }
 525                   else if (idet == 1)
 526                     {
 527                       deltaE = fCPV[ism][jrow][kcol];
 528                       detno = 1;
 529                     }
 530                   if (deltaE > 0.)
 531                     {
 532                       count_digit++;
 533                       AddDigit(trno,detno,ism,cellno,deltaE);
 534                     }
 535                 } // column loop
 536             } // row    loop
 537         } // supermodule loop
 538       treeD->Fill();
 539       ResetDigit();
 540     } // detector loop
 541
 542   pmdloader->WriteDigits("OVERWRITE");
 543
 544   ResetCellADC();
 545
 546   //  cout << " -------- End of Hits2Digit ----------- " << endl;
 547 }
 548
 549
 550 void AliPMDDigitizer::SDigits2Digits(Int_t ievt)
 551 {
 552   //  cout << " -------- Beginning of SDigits2Digit ----------- " << endl;
 553   fRunLoader->GetEvent(ievt);
 554
 555   treeS = pmdloader->TreeS();
 556   AliPMDsdigit  *pmdsdigit;
 557   TBranch *branch = treeS->GetBranch("PMDSDigit");
 558   branch->SetAddress(&fSDigits);
 559
 560   treeD = pmdloader->TreeD();
 561   if (treeD == 0x0)
 562     {
 563       pmdloader->MakeTree("D");
 564       treeD = pmdloader->TreeD();
 565     }
 566   Int_t bufsize = 16000;
 567   treeD->Branch("PMDDigit", &fDigits, bufsize);
 568
 569   Int_t   trno, det, smn;
 570   Int_t   cellno;
 571   Float_t edep, adc;
 572
 573   Int_t nmodules = (Int_t) treeS->GetEntries();
 574
 575   for (Int_t imodule = 0; imodule < nmodules; imodule++)
 576     {
 577       treeS->GetEntry(imodule);
 578       Int_t nentries = fSDigits->GetLast();
 579       //cout << " nentries = " << nentries << endl;
 580       for (Int_t ient = 0; ient < nentries+1; ient++)
 581         {
 582           pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
 583           trno   = pmdsdigit->GetTrackNumber();
 584           det    = pmdsdigit->GetDetector();
 585           smn    = pmdsdigit->GetSMNumber();
 586           cellno = pmdsdigit->GetCellNumber();
 587           edep   = pmdsdigit->GetCellEdep();
 588
 589           MeV2ADC(edep,adc);
 590
 591           AddDigit(trno,det,smn,cellno,adc);
 592         }
 593       treeD->Fill();
 594       ResetDigit();
 595     }
 596   pmdloader->WriteDigits("OVERWRITE");
 597   //  cout << " -------- End of SDigits2Digit ----------- " << endl;
 598 }
 599
 600 void AliPMDDigitizer::TrackAssignment2Cell()
 601 {
 602   // To be checked again
 603   //
 604   // This blocks assign the cell id when there is
 605   // multiple tracks in a cell according to the
 606   // energy deposition
 607   //
 608
 609   Int_t i, j, k;
 610
 611   Float_t *frac_edp;
 612   Float_t *tr_edp;
 613   Int_t *status;
 614   Int_t   ****PMDTrack;
 615   Float_t ****PMDEdep;
 616
 617   PMDTrack = new Int_t ***[27];
 618   PMDEdep  = new Float_t ***[27];
 619   for (i=0; i<fTotSM; i++)
 620     {
 621       PMDTrack[i] = new Int_t **[72];
 622       PMDEdep[i]  = new Float_t **[72];
 623     }
 624
 625   for (i = 0; i < fTotSM; i++)
 626     {
 627       for (j = 0; j < fNCell; j++)
 628         {
 629           PMDTrack[i][j] = new Int_t *[72];
 630           PMDEdep[i][j]  = new Float_t *[72];
 631         }
 632     }
 633
 634   for (i = 0; i < fTotSM; i++)
 635     {
 636       for (j = 0; j < fNCell; j++)
 637         {
 638           for (k = 0; k < fNCell; k++)
 639             {
 640               Int_t nn = fPMDCounter[i][j][k];
 641               if(nn > 0)
 642                 {
 643                   PMDTrack[i][j][k] = new Int_t[nn];
 644                   PMDEdep[i][j][k] = new Float_t[nn];
 645                 }
 646               else
 647                 {
 648                   nn = 1;
 649                   PMDTrack[i][j][k] = new Int_t[nn];
 650                   PMDEdep[i][j][k] = new Float_t[nn];
 651                 }
 652               fPMDCounter[i][j][k] = 0;
 653             }
 654         }
 655     }
 656
 657
 658   Int_t nentries = fCell->GetEntries();
 659
 660   Int_t   mtrackno, ism, ixp, iyp;
 661   Float_t edep;
 662
 663   for (i = 0; i < nentries; i++)
 664     {
 665       pmdcell = (AliPMDcell*)fCell->UncheckedAt(i);
 666
 667       mtrackno = pmdcell->GetTrackNumber();
 668       ism = pmdcell->GetSMNumber();
 669       ixp = pmdcell->GetX();
 670       iyp = pmdcell->GetY();
 671       edep = pmdcell->GetEdep();
 672
 673       Int_t nn = fPMDCounter[ism][ixp][iyp];
 674
 675       //      cout << " nn = " << nn << endl;
 676
 677       PMDTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
 678       PMDEdep[ism][ixp][iyp][nn] = edep;
 679       fPMDCounter[ism][ixp][iyp]++;
 680     }
 681
 682   Int_t iz, il;
 683   Int_t im, ix, iy;
 684   Int_t nn;
 685
 686   for (im=0; im<27; im++)
 687     {
 688       for (ix=0; ix<72; ix++)
 689         {
 690           for (iy=0; iy<72; iy++)
 691             {
 692               nn = fPMDCounter[im][ix][iy];
 693               if (nn > 1)
 694                 {
 695                   // This block handles if a cell is fired
 696                   // many times by many tracks
 697
 698                   status = new Int_t[nn];
 699                   for (iz = 0; iz < nn; iz++)
 700                     {
 701                       status[iz] = PMDTrack[im][ix][iy][iz];
 702                     }
 703                   sort(status,status+nn);
 704                   Int_t track_old = -99999;
 705                   Int_t track, tr_count = 0;
 706                   for (iz = 0; iz < nn; iz++)
 707                     {
 708                       track = status[iz];
 709                       if (track_old != track)
 710                         {
 711                           tr_count++;
 712                           vjunkTRN.push_back(track);
 713                         }
 714                       track_old = track;
 715                     }
 716                   delete status;
 717                   Float_t tot_edp = 0.;
 718                   tr_edp = new Float_t[tr_count];
 719                   frac_edp = new Float_t[tr_count];
 720                   for (il = 0; il < tr_count; il++)
 721                     {
 722                       tr_edp[il] = 0.;
 723                       track = vjunkTRN[il];
 724                       for (iz = 0; iz < nn; iz++)
 725                         {
 726                           if (track == PMDTrack[im][ix][iy][iz])
 727                             {
 728                               tr_edp[il] += PMDEdep[im][ix][iy][iz];
 729                             }
 730                         }
 731                       tot_edp += tr_edp[il];
 732                     }
 733                   vjunkTRN.clear();
 734                   Int_t il_old = 0;
 735                   Float_t frac_old = 0.;
 736
 737                   for (il = 0; il < tr_count; il++)
 738                     {
 739                       frac_edp[il] = tr_edp[il]/tot_edp;
 740                       if (frac_old < frac_edp[il])
 741                         {
 742                           frac_old = frac_edp[il];
 743                           il_old = il;
 744                         }
 745                     }
 746                   delete frac_edp;
 747                   delete tr_edp;
 748
 749                   fPMDTrackNo[im][ix][iy] = vjunkTRN[il_old];
 750                 }
 751               else if (nn == 1)
 752                 {
 753                   // This only handles if a cell is fired
 754                   // by only one track
 755
 756                   fPMDTrackNo[im][ix][iy] = PMDTrack[im][ix][iy][0];
 757
 758                 }
 759               else if (nn ==0)
 760                 {
 761                   // This is if no cell is fired
 762                   fPMDTrackNo[im][ix][iy] = -999;
 763                 }
 764             } // end of iy
 765         } // end of ix
 766     } // end of im
 767
 768   // Delete all the pointers
 769
 770   for (i = 0; i < fTotSM; i++)
 771     {
 772       for (j = 0; j < fNCell; j++)
 773         {
 774           for (k = 0; k < fNCell; k++)
 775             {
 776               delete [] PMDTrack[i][j][k];
 777               delete [] PMDEdep[i][j][k];
 778             }
 779         }
 780     }
 781
 782   for (i = 0; i < fTotSM; i++)
 783     {
 784       for (j = 0; j < fNCell; j++)
 785         {
 786           delete [] PMDTrack[i][j];
 787           delete [] PMDEdep[i][j];
 788         }
 789     }
 790
 791   for (i = 0; i < fTotSM; i++)
 792     {
 793       delete [] PMDTrack[i];
 794       delete [] PMDEdep[i];
 795     }
 796   delete PMDTrack;
 797   delete PMDEdep;
 798   //
 799   // End of the cell id assignment
 800   //
 801 }
 802
 803
 804 void AliPMDDigitizer::MeV2ADC(Float_t mev, Float_t & adc)
 805 {
 806   // To be done
 807
 808   adc = mev*1.;
 809 }
 810 void AliPMDDigitizer::AddSDigit(Int_t trnumber, Int_t det, Int_t smnumber,
 811   Int_t cellnumber, Float_t adc)
 812 {
 813   TClonesArray &lsdigits = *fSDigits;
 814   AliPMDsdigit *newcell;
 815   newcell = new AliPMDsdigit(trnumber,det,smnumber,cellnumber,adc);
 816   new(lsdigits[fNsdigit++]) AliPMDsdigit(newcell);
 817   delete newcell;
 818 }
 819
 820 void AliPMDDigitizer::AddDigit(Int_t trnumber, Int_t det, Int_t smnumber,
 821   Int_t cellnumber, Float_t adc)
 822 {
 823   TClonesArray &ldigits = *fDigits;
 824   AliPMDdigit *newcell;
 825   newcell = new AliPMDdigit(trnumber,det,smnumber,cellnumber,adc);
 826   new(ldigits[fNdigit++]) AliPMDdigit(newcell);
 827   delete newcell;
 828 }
 829
 830 Int_t AliPMDDigitizer::Convert2RealSMNumber(Int_t smnumber1)
 831 {
 832   Int_t smnumber = -999;
 833
 834   if (smnumber1==1)  smnumber =  1;
 835   if (smnumber1==2)  smnumber = 10;
 836   if (smnumber1==3)  smnumber = 19;
 837   if (smnumber1==4)  smnumber =  1;
 838   if (smnumber1==5)  smnumber = 10;
 839   if (smnumber1==6)  smnumber = 19;
 840   if (smnumber1==7)  smnumber =  2;
 841   if (smnumber1==8)  smnumber =  3;
 842   if (smnumber1==9)  smnumber =  4;
 843   if (smnumber1==10) smnumber =  5;
 844   if (smnumber1==11) smnumber =  6;
 845   if (smnumber1==12) smnumber =  7;
 846   if (smnumber1==13) smnumber =  8;
 847   if (smnumber1==14) smnumber =  9;
 848   if (smnumber1==15) smnumber = 11;
 849   if (smnumber1==16) smnumber = 12;
 850   if (smnumber1==17) smnumber = 13;
 851   if (smnumber1==18) smnumber = 14;
 852   if (smnumber1==19) smnumber = 15;
 853   if (smnumber1==20) smnumber = 16;
 854   if (smnumber1==21) smnumber = 17;
 855   if (smnumber1==22) smnumber = 18;
 856   if (smnumber1==23) smnumber = 20;
 857   if (smnumber1==24) smnumber = 21;
 858   if (smnumber1==25) smnumber = 22;
 859   if (smnumber1==26) smnumber = 23;
 860   if (smnumber1==27) smnumber = 24;
 861   if (smnumber1==28) smnumber = 25;
 862   if (smnumber1==29) smnumber = 26;
 863   if (smnumber1==30) smnumber = 27;
 864
 865   return smnumber;
 866 }
 867 void AliPMDDigitizer::SetZPosition(Float_t zpos)
 868 {
 869   fZPos = zpos;
 870 }
 871 Float_t AliPMDDigitizer::GetZPosition() const
 872 {
 873   return fZPos;
 874 }
 875
 876 void AliPMDDigitizer::ResetCell()
 877 {
 878   fCell->Clear();
 879   for (Int_t i = 0; i < fTotSM; i++)
 880     {
 881       for (Int_t j = 0; j < fNCell; j++)
 882         {
 883           for (Int_t k = 0; k < fNCell; k++)
 884             {
 885               fPMDCounter[i][j][k] = 0;
 886             }
 887         }
 888     }
 889 }
 890 void AliPMDDigitizer::ResetSDigit()
 891 {
 892   fNsdigit = 0;
 893   if (fSDigits) fSDigits->Clear();
 894 }
 895 void AliPMDDigitizer::ResetDigit()
 896 {
 897   fNdigit = 0;
 898   if (fDigits) fDigits->Clear();
 899 }
 900
 901 void AliPMDDigitizer::ResetCellADC()
 902 {
 903   for (Int_t i = 0; i < fTotSM; i++)
 904     {
 905       for (Int_t j = 0; j < fNCell; j++)
 906         {
 907           for (Int_t k = 0; k < fNCell; k++)
 908             {
 909               fCPV[i][j][k] = 0.;
 910               fPMD[i][j][k] = 0.;
 911             }
 912         }
 913     }
 914 }
 915
 916 void AliPMDDigitizer::UnLoad(Option_t *option)
 917 {
 918   const char *cS = strstr(option,"S");
 919   const char *cD = strstr(option,"D");
 920
 921   fRunLoader->UnloadgAlice();
 922   fRunLoader->UnloadHeader();
 923   fRunLoader->UnloadKinematics();
 924
 925   if (cS)
 926     {
 927       pmdloader->UnloadHits();
 928     }
 929   if (cD)
 930     {
 931       pmdloader->UnloadHits();
 932       pmdloader->UnloadSDigits();
 933     }
 934 }