PMD/AliPMDDigitizer.cxx

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