PMD/AliPMDDigitizer.cxx

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