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