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