]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PMD/AliPMDDigitizer.cxx
Removing the tasks from the digitization (Ruben)
[u/mrichter/AliRoot.git] / PMD / AliPMDDigitizer.cxx
index a47041a7837a9d95b4a40db88a8e9cf53c1e18e6..d77c940d4097a40d4e2ba42434259466fcbee49c 100644 (file)
@@ -1,11 +1,26 @@
+/***************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
 //-----------------------------------------------------//
 //                                                     //
-//  Source File : PMDDigitization.cxx, Version 00      //
+//  Source File : PMDDigitizer.cxx, Version 00         //
 //                                                     //
 //  Date   : September 20 2002                         //
 //                                                     //
 //-----------------------------------------------------//
 
+#include <Riostream.h>
 #include <TBRIK.h>
 #include <TNode.h>
 #include <TTree.h>
 #include <TFile.h>
 #include <TNtuple.h>
 #include <TParticle.h>
+#include <TRandom.h>
 
+#include "AliLog.h"
 #include "AliRun.h"
-#include "AliPMD.h"
 #include "AliHit.h"
 #include "AliDetector.h"
 #include "AliRunLoader.h"
 #include "AliLoader.h"
 #include "AliConfig.h"
 #include "AliMagF.h"
-#include "AliRunDigitizer.h"
+#include "AliDigitizationInput.h"
+#include "AliDigitizer.h"
 #include "AliHeader.h"
+#include "AliCDBManager.h"
+#include "AliCDBStorage.h"
+#include "AliCDBEntry.h"
+#include "AliMC.h"
 
+#include "AliPMD.h"
+#include "AliPMDhit.h"
 #include "AliPMDcell.h"
 #include "AliPMDsdigit.h"
 #include "AliPMDdigit.h"
+#include "AliPMDCalibData.h"
+#include "AliPMDPedestal.h"
 #include "AliPMDDigitizer.h"
-#include "AliPMDClustering.h"
-#include "AliPMDrecpoint.h"
 
 
 ClassImp(AliPMDDigitizer)
-//
-// Constructor
-//
-AliPMDDigitizer::AliPMDDigitizer()
-{
-  if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);  
-  fNsdigit = 0;
-  if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);  
-  fNdigit = 0;
 
-  for (Int_t i = 0; i < fTotSM; i++)
+AliPMDDigitizer::AliPMDDigitizer() :
+  fRunLoader(0),
+  fPMDHit(0),
+  fPMD(0),
+  fPMDLoader(0),
+  fCalibGain(GetCalibGain()),
+  fCalibPed(GetCalibPed()),
+  fSDigits(0),
+  fDigits(0),
+  fCPVCell(0),
+  fCell(0),
+  fNsdigit(0),
+  fNdigit(0),
+  fDetNo(0),
+  fZPos(361.5)   // in units of cm, default position of PMD
+{
+  // Default Constructor
+  //
+  for (Int_t i = 0; i < fgkTotUM; i++)
     {
-      for (Int_t j = 0; j < fNCell; j++)
+      for (Int_t j = 0; j < fgkRow; j++)
        {
-         for (Int_t k = 0; k < fNCell; k++)
+         for (Int_t k = 0; k < fgkCol; k++)
            {
-             fCPV[i][j][k] = 0.; 
-             fPMD[i][j][k] = 0.; 
+             fCPV[i][j][k]         = 0.;
+             fPRE[i][j][k]         = 0.;
+             fCPVCounter[i][j][k]  =  0; 
+             fPRECounter[i][j][k]  =  0;
+             fCPVTrackNo[i][j][k]  = -1;
+             fPRETrackNo[i][j][k]  = -1;
+             fCPVTrackPid[i][j][k] = -1;
+             fPRETrackPid[i][j][k] = -1;
            }
        }
     }
 
-  if (!fCell) fCell = new TObjArray();
+}
+//____________________________________________________________________________
+AliPMDDigitizer::AliPMDDigitizer(const AliPMDDigitizer& digitizer):
+  AliDigitizer(digitizer),
+  fRunLoader(0),
+  fPMDHit(0),
+  fPMD(0),
+  fPMDLoader(0),
+  fCalibGain(GetCalibGain()),
+  fCalibPed(GetCalibPed()),
+  fSDigits(0),
+  fDigits(0),
+  fCPVCell(0),
+  fCell(0),
+  fNsdigit(0),
+  fNdigit(0),
+  fDetNo(0),
+  fZPos(361.5)   // in units of cm, default position of PMD
+{
+  // copy constructor
+  AliError("Copy constructor not allowed ");
+  
+}
+//____________________________________________________________________________
+AliPMDDigitizer & AliPMDDigitizer::operator=(const AliPMDDigitizer& /*digitizer*/)
+{
+  // Assignment operator
+  AliError("Assignement operator not allowed ");
+
+  return *this;
+}
+//____________________________________________________________________________
+AliPMDDigitizer::AliPMDDigitizer(AliDigitizationInput* digInput):
+  AliDigitizer(digInput),
+  fRunLoader(0),
+  fPMDHit(0),
+  fPMD(0),
+  fPMDLoader(0),
+  fCalibGain(GetCalibGain()),
+  fCalibPed(GetCalibPed()),
+  fSDigits(new TClonesArray("AliPMDsdigit", 1000)),
+  fDigits(new TClonesArray("AliPMDdigit", 1000)),
+  fCPVCell(0),
+  fCell(0),
+  fNsdigit(0),
+  fNdigit(0),
+  fDetNo(0),
+  fZPos(361.5)// in units of cm, This is the default position of PMD
+{
+  // ctor which should be used
 
-  fZPos = 361.5; // in units of cm, This is the default position of PMD
 
+  for (Int_t i = 0; i < fgkTotUM; i++)
+    {
+      for (Int_t j = 0; j < fgkRow; j++)
+       {
+         for (Int_t k = 0; k < fgkCol; k++)
+           {
+             fCPV[i][j][k]         = 0.;
+             fPRE[i][j][k]         = 0.;
+             fCPVCounter[i][j][k]  =  0; 
+             fPRECounter[i][j][k]  =  0;
+             fCPVTrackNo[i][j][k]  = -1;
+             fPRETrackNo[i][j][k]  = -1;
+             fCPVTrackPid[i][j][k] = -1;
+             fPRETrackPid[i][j][k] = -1;
+           }
+       }
+    }
 }
+
+//____________________________________________________________________________
 AliPMDDigitizer::~AliPMDDigitizer()
 {
-  delete fSDigits;
-  delete fDigits;
-  delete fCell;
+  // Default Destructor
+  //
+  if (fSDigits) {
+    fSDigits->Delete();
+    delete fSDigits;
+    fSDigits=0;
+  }
+  if (fDigits) {
+    fDigits->Delete();
+    delete fDigits;
+    fDigits=0;
+  }
+  fCPVCell.Delete();
+  fCell.Delete();
 }
 //
 // Member functions
 //
-void AliPMDDigitizer::OpengAliceFile(Char_t *file, Option_t *option)
+//____________________________________________________________________________
+void AliPMDDigitizer::OpengAliceFile(const char *file, Option_t *option)
 {
+  // Loads galice.root file and corresponding header, kinematics
+  // hits and sdigits or digits depending on the option
+  //
 
-  fRunLoader = AliRunLoader::Open(file,AliConfig::fgkDefaultEventFolderName,
-                                 "UPDATE");
+  TString evfoldname = AliConfig::GetDefaultEventFolderName();
+  fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
+  if (!fRunLoader)
+      fRunLoader = AliRunLoader::Open(file,AliConfig::GetDefaultEventFolderName(), "UPDATE");
   
   if (!fRunLoader)
    {
-     Error("Open","Can not open session for file %s.",file);
+     AliError(Form("Can not open session for file %s.",file));
    }
-  
-  fRunLoader->LoadgAlice();
-  fRunLoader->LoadHeader();
-  fRunLoader->LoadKinematics();
 
-  gAlice = fRunLoader->GetAliRun();
+  const char *cHS = strstr(option,"HS");
+  const char *cHD = strstr(option,"HD");
+  const char *cSD = strstr(option,"SD");
   
-  if (gAlice)
-    {
-      printf("<AliPMDdigitizer::Open> ");
-      printf("AliRun object found on file.\n");
-    }
-  else
+  if(cHS || cHD)
     {
-      printf("<AliPMDdigitizer::Open> ");
-      printf("Could not find AliRun object.\n");
+      if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
+      if (!fRunLoader->TreeE()) fRunLoader->LoadHeader();
+      if (!fRunLoader->TreeK()) fRunLoader->LoadKinematics();
+  
+      gAlice = fRunLoader->GetAliRun();
+  
+      if (gAlice)
+       {
+         AliDebug(1,"Alirun object found");
+       }
+      else
+       {
+         AliError("Could not found Alirun object");
+       }
+  
+      fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
     }
 
-  PMD  = (AliPMD*)gAlice->GetDetector("PMD");
-  pmdloader = fRunLoader->GetLoader("PMDLoader");
-  if (pmdloader == 0x0)
+  fPMDLoader = fRunLoader->GetLoader("PMDLoader");
+  if (fPMDLoader == 0x0)
     {
-      cerr<<"Hits2Digits : Can not find PMD or PMDLoader\n";
+      AliError("Can not find PMDLoader");
     }
 
-  const char *cHS = strstr(option,"HS");
-  const char *cHD = strstr(option,"HD");
-  const char *cSD = strstr(option,"SD");
 
   if (cHS)
     {
-      pmdloader->LoadHits("READ");
-      pmdloader->LoadSDigits("recreate");
+      fPMDLoader->LoadHits("READ");
+      fPMDLoader->LoadSDigits("recreate");
     }
   else if (cHD)
     {
-      pmdloader->LoadHits("READ");
-      pmdloader->LoadDigits("recreate");
+      fPMDLoader->LoadHits("READ");
+      fPMDLoader->LoadDigits("recreate");
     }
   else if (cSD)
     {
-      pmdloader->LoadSDigits("READ");
-      pmdloader->LoadDigits("recreate");
+      fPMDLoader->LoadSDigits("READ");
+      fPMDLoader->LoadDigits("recreate");
     }
-
 }
+//____________________________________________________________________________
 void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
 {
-  cout << " -------- Beginning of Hits2SDigits ----------- " << endl;
-
-  Int_t kPi0 = 111;
-  Int_t kGamma = 22;
-  Int_t npmd;
-  Int_t trackno;
-  Int_t cellnumber;
-  Int_t smnumber1;
-  Int_t smnumber;
-  Int_t trackpid;
-  Int_t mtrackno;
-  Int_t mtrackpid;
-
-  Float_t xPos, yPos, zPos;
-  Int_t xpad, ypad;
-  Float_t edep;
+  // This reads the PMD Hits tree and assigns the right track number
+  // to a cell and stores in the summable digits tree
+  //
+
+  const Int_t kPi0 = 111;
+  const Int_t kGamma = 22;
+  Int_t npmd = 0;
+  Int_t trackno = 0;
+  Int_t smnumber = 0;
+  Int_t trackpid = 0;
+  Int_t mtrackno = 0;
+  Int_t mtrackpid = 0;
+
+  Float_t xPos = 0., yPos = 0., zPos = 0.;
+  Int_t xpad = -1, ypad = -1;
+  Float_t edep = 0.;
   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
 
-  
+
+  if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
   ResetSDigit();
 
-  printf("Event Number =  %d \n",ievt); 
+  AliDebug(1,Form("Event Number = %d",ievt));
   Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
-  printf("Number of Particles = %d \n", nparticles);
+  AliDebug(1,Form("Number of Particles = %d",nparticles));
+
+  //
+
   fRunLoader->GetEvent(ievt);
-  Particles = gAlice->Particles();
   // ------------------------------------------------------- //
   // Pointer to specific detector hits.
   // Get pointers to Alice detectors and Hits containers
 
-  treeH = pmdloader->TreeH();
+  TTree* treeH = fPMDLoader->TreeH();
   
   Int_t ntracks    = (Int_t) treeH->GetEntries();
-  printf("Number of Tracks in the TreeH = %d \n", ntracks);
-
-  treeS = pmdloader->TreeS();
+  AliDebug(1,Form("Number of Tracks in the TreeH = %d", ntracks));
+  TTree* treeS = fPMDLoader->TreeS();
   if (treeS == 0x0)
     {
-      pmdloader->MakeTree("S");
-      treeS = pmdloader->TreeS();
+      fPMDLoader->MakeTree("S");
+      treeS = fPMDLoader->TreeS();
     }
   Int_t bufsize = 16000;
-  treeS->Branch("PMDSDigit", &fSDigits, bufsize); 
+  treeS->Branch("PMDSDigit", &fSDigits, bufsize);
   
-  if (PMD) PMDhits   = PMD->Hits();
+  TClonesArray* hits = 0;
+  if (fPMD) hits = fPMD->Hits();
 
   // Start loop on tracks in the hits containers
 
-  
-  for (Int_t track=0; track<ntracks;track++) 
+  for (Int_t track=0; track<ntracks;track++)
     {
-      gAlice->ResetHits();
+      gAlice->GetMCApp()->ResetHits();
       treeH->GetEvent(track);
-      
-      if (PMD) 
+      if (fPMD)
        {
-         npmd = PMDhits->GetEntriesFast();
-         for (int ipmd = 0; ipmd < npmd; ipmd++) 
+         npmd = hits->GetEntriesFast();
+         for (Int_t ipmd = 0; ipmd < npmd; ipmd++)
            {
-             pmdHit = (AliPMDhit*) PMDhits->UncheckedAt(ipmd);
-             trackno = pmdHit->GetTrack();
-
+             fPMDHit = (AliPMDhit*) hits->UncheckedAt(ipmd);
+             trackno = fPMDHit->GetTrack();
              //  get kinematics of the particles
-             
-             particle = gAlice->Particle(trackno);
-             trackpid  = particle->GetPdgCode();
-
-             Int_t igatr = -999;
-             Int_t ichtr = -999;
 
+             TParticle* mparticle = gAlice->GetMCApp()->Particle(trackno);
+             trackpid  = mparticle->GetPdgCode();
+             Int_t  ks = mparticle->GetStatusCode();
              Int_t imo;
-             Int_t igen = 0;
-             Int_t id_mo = -999;
-             TParticle*  mparticle = particle;
+             Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
              
-             while((imo = mparticle->GetFirstMother()) >= 0)
+             if (mparticle->GetFirstMother() == -1)
                {
-                 igen++;
-                 mparticle =  gAlice->Particle(imo);
-                 id_mo = mparticle->GetPdgCode();
-                 
+                 tracknoOld  = trackno;
+                 trackpidOld = trackpid;
+                 statusOld   = -1;
+               }
+             Int_t igstatus = 0;
+
+             Int_t trnotemp = trackno;    // Modified on 25th Nov 2009
+             if(ks==1||(imo = mparticle->GetFirstMother())<0 ){
+               vx = mparticle->Vx();
+               vy = mparticle->Vy();
+               vz = mparticle->Vz();
+               
+               if(trackpid==kGamma||trackpid==11||trackpid==-11||
+                  trackpid==kPi0)igstatus=1;
+             }
+             
+             
+             while(((imo = mparticle->GetFirstMother()) >= 0)&& 
+                   (ks = mparticle->GetStatusCode() <1) )
+               {
+                 mparticle =  gAlice->GetMCApp()->Particle(imo);
+                 trackpid = mparticle->GetPdgCode();
+                 ks = mparticle->GetStatusCode();
                  vx = mparticle->Vx();
                  vy = mparticle->Vy();
                  vz = mparticle->Vz();
-               
-                 //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, id_mo, vx, vy, vz);
-                 //fprintf(ftest1,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, id_mo, vx, vy, vz);
                  
-                 if (id_mo == kGamma && vx == 0. && vy == 0. && vz == 0.)
+                 // Modified on 25th Nov 2009
+
+                 trnotemp = trackno;
+                  if(trackpid == 111)
+                   {
+                     trackno = trnotemp;
+                   }
+                  if(trackpid != 111)
                    {
-                     igatr = imo;
+                     trackno=imo;
                    }
-                 ichtr = imo;
+                 // end of modification on 25th Nov 2009
                }
+        
+             if(trackpid==kGamma||trackpid==11||trackpid==-11||
+                trackpid==kPi0)igstatus=1;
+             mtrackpid=trackpid;
+             mtrackno=trackno;
+             trackpid=trackpidOld;
+             trackno=tracknoOld;
+             
+             //-----------------end of modification----------------
+             Float_t ptime = fPMDHit->GetTime()*1e6;   // time in microsec
+             if (ptime < 0. || ptime > 1.2) continue;  
+
+             xPos = fPMDHit->X();
+             yPos = fPMDHit->Y();
+             zPos = fPMDHit->Z();
+
+             edep       = fPMDHit->GetEnergy();
+             Int_t vol1 = fPMDHit->GetVolume(1); // Column
+             Int_t vol2 = fPMDHit->GetVolume(2); // Row
+             Int_t vol7 = fPMDHit->GetVolume(4); // Serial Module No
 
-             if (id_mo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
+
+             // -----------------------------------------//
+             // In new geometry after adding electronics //
+             // For Super Module 1 & 2                   //
+             //  nrow = 48, ncol = 96                    //
+             // For Super Module 3 & 4                   //
+             //  nrow = 96, ncol = 48                    //
+             // -----------------------------------------//
+
+             if (vol7 < 24)
                {
-                 mtrackno = igatr;
-                 mtrackpid = kGamma;
+                 smnumber = vol7;
                }
              else
                {
-                 mtrackno  = ichtr;
-                 mtrackpid = id_mo;
+                 smnumber = vol7 - 24;
                }
+             Int_t vol8 = smnumber/6 + 1;   // fake supermodule
 
-             xPos = pmdHit->X();
-             yPos = pmdHit->Y();
-             zPos = pmdHit->Z();
-             cellnumber = pmdHit->fVolume[1];
-             smnumber1  = pmdHit->fVolume[4];
-             edep       = pmdHit->fEnergy;
-             
-             if (smnumber1 > 3 && smnumber1 <= 6)
+             if (vol8 == 1 || vol8 == 2)
                {
-                 Int_t ny = (cellnumber-1)/48 + 1;
-                 Int_t nx = cellnumber - (ny-1)*48;
-                 Int_t nx1 = nx + 24;
-                 Int_t ny1 = ny + 48;
-                 Int_t ncell1 = (ny1 - 1)* 72 + nx1;
-                 cellnumber = ncell1;
+                 xpad = vol2;
+                 ypad = vol1;
+               }
+             else if (vol8 == 3 || vol8 == 4)
+               {
+                 xpad = vol1;
+                 ypad = vol2;
                }
-             smnumber = Convert2RealSMNumber(smnumber1);
-             ypad = (cellnumber - 1)/fNCell + 1;
-             xpad = cellnumber - (ypad-1)*fNCell;
 
-             //cout << "-zpos = " << -zPos << endl;
+             AliDebug(2,Form("Zposition = %f Edeposition = %f",zPos,edep));
 
-             if (-zPos < fZPos)
+             if (vol7 < 24)
                {
-                 // CPV
-                 fDetNo = 1;
+                 // PRE
+                 fDetNo = 0;
                }
-             else if (-zPos > fZPos)
+             else
                {
-                 // PMD
-                 fDetNo = 0;
+                 // CPV
+                 fDetNo = 1;
                }
 
+             Int_t smn = smnumber;
+             Int_t ixx = xpad     - 1;
+             Int_t iyy = ypad     - 1;
              if (fDetNo == 0)
                {
-                 fPMD[smnumber-1][xpad-1][ypad-1] += edep;
-                 fPMDCounter[smnumber-1][xpad-1][ypad-1]++;
-                 Int_t smn = smnumber - 1;
-                 Int_t ixx = xpad     - 1;
-                 Int_t iyy = ypad     - 1;
+                 fPRE[smn][ixx][iyy] += edep;
+                 fPRECounter[smn][ixx][iyy]++;
 
-                 pmdcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
-
-                 fCell->Add(pmdcell);
+                 AliPMDcell* cell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
+                 fCell.Add(cell);
                }
              else if(fDetNo == 1)
                {
-                 fCPV[smnumber-1][xpad-1][ypad-1] += edep;
-                 fCPVTrackNo[smnumber-1][xpad-1][ypad-1] = mtrackno;
+                 fCPV[smn][ixx][iyy] += edep;
+                 fCPVCounter[smn][ixx][iyy]++;
+                 AliPMDcell* cpvcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep); 
+                 fCPVCell.Add(cpvcell);
                }
            }
        }
     } // Track Loop ended
-  
-  
+  TrackAssignment2CPVCell();
   TrackAssignment2Cell();
-  
   ResetCell();
 
   Float_t deltaE      = 0.;
   Int_t   detno       = 0;
-  Int_t   count_digit = 0;
   Int_t   trno        = -1;
-  Int_t   cellno      = 0;
+  Int_t   trpid       = -99;
 
   for (Int_t idet = 0; idet < 2; idet++)
     {
-      for (Int_t ism = 0; ism < fTotSM; ism++)
+      for (Int_t ism = 0; ism < fgkTotUM; ism++)
        {
-         for (Int_t jrow = 0; jrow < fNCell; jrow++)
+         for (Int_t jrow = 0; jrow < fgkRow; jrow++)
            {
-             for (Int_t kcol = 0; kcol < fNCell; kcol++)
+             for (Int_t kcol = 0; kcol < fgkCol; kcol++)
                {
-                 cellno = jrow + kcol*fNCell;
                  if (idet == 0)
                    {
-                     deltaE = fPMD[ism][jrow][kcol];
-                     trno   = fPMDTrackNo[ism][jrow][kcol];
+                     deltaE = fPRE[ism][jrow][kcol];
+                     trno   = fPRETrackNo[ism][jrow][kcol];
                      detno = 0;
                    }
                  else if (idet == 1)
@@ -324,8 +483,10 @@ void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
                    }
                  if (deltaE > 0.)
                    {
-                     count_digit++;
-                     AddSDigit(trno,detno,ism,cellno,deltaE);
+                     // Natasha
+                     TParticle *mparticle = gAlice->GetMCApp()->Particle(trno);
+                     trpid = mparticle->GetPdgCode();
+                     AddSDigit(trno,trpid,detno,ism,jrow,kcol,deltaE);
                    }
                }
            }
@@ -333,432 +494,903 @@ void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
          ResetSDigit();
        }
     }
-  pmdloader->WriteSDigits("OVERWRITE");
-
+  fPMDLoader->WriteSDigits("OVERWRITE");
   ResetCellADC();
-
-  //  cout << " -------- End of Hits2SDigit ----------- " << endl;
 }
+//____________________________________________________________________________
 
 void AliPMDDigitizer::Hits2Digits(Int_t ievt)
 {
-  Int_t kPi0 = 111;
-  Int_t kGamma = 22;
-  Int_t npmd;
-  Int_t trackno;
-  Int_t cellnumber;
-  Int_t smnumber1;
-  Int_t smnumber;
-  Int_t trackpid;
-  Int_t mtrackno;
-  Int_t mtrackpid;
-
-  Float_t xPos, yPos, zPos;
-  Int_t xpad, ypad;
-  Float_t edep;
+  // This reads the PMD Hits tree and assigns the right track number
+  // to a cell and stores in the digits tree
+  //
+  const Int_t kPi0 = 111;
+  const Int_t kGamma = 22;
+  Int_t npmd      = 0;
+  Int_t trackno   = 0;
+  Int_t smnumber  = 0;
+  Int_t trackpid  = 0;
+  Int_t mtrackno  = 0;
+  Int_t mtrackpid = 0;
+
+  Float_t xPos = 0., yPos = 0., zPos = 0.;
+  Int_t xpad = -1, ypad = -1;
+  Float_t edep = 0.;
   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
 
-  
+  if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
   ResetDigit();
 
-  printf("Event Number =  %d \n",ievt); 
-
+  AliDebug(1,Form("Event Number =  %d",ievt));
   Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
-  printf("Number of Particles = %d \n", nparticles);
+  AliDebug(1,Form("Number of Particles = %d", nparticles));
+
   fRunLoader->GetEvent(ievt);
-  Particles = gAlice->Particles();
   // ------------------------------------------------------- //
   // Pointer to specific detector hits.
   // Get pointers to Alice detectors and Hits containers
 
-  PMD  = (AliPMD*)gAlice->GetDetector("PMD");
-  pmdloader = fRunLoader->GetLoader("PMDLoader");
+  fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
+  fPMDLoader = fRunLoader->GetLoader("PMDLoader");
 
-  if (pmdloader == 0x0)
+  if (fPMDLoader == 0x0)
     {
-      cerr<<"Hits2Digits method : Can not find PMD or PMDLoader\n";
+      AliError("Can not find PMD or PMDLoader");
     }
-  treeH = pmdloader->TreeH();
+  TTree* treeH = fPMDLoader->TreeH();
   Int_t ntracks    = (Int_t) treeH->GetEntries();
-  printf("Number of Tracks in the TreeH = %d \n", ntracks);
-  pmdloader->LoadDigits("recreate");
-  treeD = pmdloader->TreeD();
+  AliDebug(1,Form("Number of Tracks in the TreeH = %d", ntracks));
+  fPMDLoader->LoadDigits("recreate");
+  TTree* treeD = fPMDLoader->TreeD();
   if (treeD == 0x0)
     {
-      pmdloader->MakeTree("D");
-      treeD = pmdloader->TreeD();
+      fPMDLoader->MakeTree("D");
+      treeD = fPMDLoader->TreeD();
     }
   Int_t bufsize = 16000;
-  treeD->Branch("PMDDigit", &fDigits, bufsize); 
+  treeD->Branch("PMDDigit", &fDigits, bufsize);
   
-  if (PMD) PMDhits   = PMD->Hits();
+  TClonesArray* hits = 0;
+  if (fPMD) hits = fPMD->Hits();
 
   // Start loop on tracks in the hits containers
 
-  for (Int_t track=0; track<ntracks;track++) 
+  for (Int_t track=0; track<ntracks;track++)
     {
-      gAlice->ResetHits();
+      gAlice->GetMCApp()->ResetHits();
       treeH->GetEvent(track);
       
-      if (PMD) 
+      if (fPMD)
        {
-         npmd = PMDhits->GetEntriesFast();
-         for (int ipmd = 0; ipmd < npmd; ipmd++) 
+         npmd = hits->GetEntriesFast();
+         for (Int_t ipmd = 0; ipmd < npmd; ipmd++)
            {
-             pmdHit = (AliPMDhit*) PMDhits->UncheckedAt(ipmd);
-             trackno = pmdHit->GetTrack();
+             fPMDHit = (AliPMDhit*) hits->UncheckedAt(ipmd);
+             trackno = fPMDHit->GetTrack();
              
              //  get kinematics of the particles
              
-             particle = gAlice->Particle(trackno);
-             trackpid  = particle->GetPdgCode();
+             TParticle* mparticle = gAlice->GetMCApp()->Particle(trackno);
+             trackpid  = mparticle->GetPdgCode();
+             Int_t  ks = mparticle->GetStatusCode();
+             Int_t imo;
+             Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
+             if (mparticle->GetFirstMother() == -1)
+               {
+                 tracknoOld  = trackno;
+                 trackpidOld = trackpid;
+                 statusOld   = -1;
+               }
 
-             Int_t igatr = -999;
-             Int_t ichtr = -999;
+             Int_t igstatus = 0;
 
-             Int_t imo;
-             Int_t igen = 0;
-             Int_t id_mo = -999;
-             TParticle*  mparticle = particle;
+             Int_t trnotemp = trackno;   // modified on 25th Nov 2009
+             if(ks==1||(imo = mparticle->GetFirstMother())<0 ){
+               vx = mparticle->Vx();
+               vy = mparticle->Vy();
+               vz = mparticle->Vz();
+               
+               if(trackpid==kGamma||trackpid==11||trackpid==-11||trackpid==kPi0)
+                 igstatus=1;
+             }
+             
              
-             while((imo = mparticle->GetFirstMother()) >= 0)
+             while(((imo = mparticle->GetFirstMother()) >= 0)&& 
+                   (ks = mparticle->GetStatusCode() <1) )
                {
-                 igen++;
-                 mparticle =  gAlice->Particle(imo);
-                 id_mo = mparticle->GetPdgCode();
-                 
+                 mparticle =  gAlice->GetMCApp()->Particle(imo);
+                 trackpid = mparticle->GetPdgCode();
+                 ks = mparticle->GetStatusCode();
                  vx = mparticle->Vx();
                  vy = mparticle->Vy();
                  vz = mparticle->Vz();
-               
-                 //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, id_mo, vx, vy, vz);
-                 //fprintf(ftest1,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, id_mo, vx, vy, vz);
                  
-                 if (id_mo == kGamma && vx == 0. && vy == 0. && vz == 0.)
+                 // Modified on 25th Nov 2009
+                 
+                 trnotemp = trackno;
+                  if(trackpid == 111)
                    {
-                     igatr = imo;
+                     trackno = trnotemp;
+                   }
+                  if(trackpid != 111)
+                   {
+                     trackno=imo;
                    }
-                 ichtr = imo;
                }
-
-             if (id_mo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
+        
+             if(trackpid==kGamma||trackpid==11||trackpid==-11||trackpid==kPi0)
+               igstatus=1;
+             mtrackpid=trackpid;
+             mtrackno=trackno;
+             trackpid=trackpidOld;
+             trackno=tracknoOld;
+
+             Float_t ptime = fPMDHit->GetTime()*1e6;
+             if (ptime < 0. || ptime > 1.2) continue;
+             
+             xPos = fPMDHit->X();
+             yPos = fPMDHit->Y();
+             zPos = fPMDHit->Z();
+             edep       = fPMDHit->GetEnergy();
+             Int_t vol1 = fPMDHit->GetVolume(1); // Column
+             Int_t vol2 = fPMDHit->GetVolume(2); // Row
+             Int_t vol7 = fPMDHit->GetVolume(4); // Serial Module No
+
+             // -----------------------------------------//
+             // In new geometry after adding electronics //
+             // For Super Module 1 & 2                   //
+             //  nrow = 48, ncol = 96                    //
+             // For Super Module 3 & 4                   //
+             //  nrow = 96, ncol = 48                    //
+             // -----------------------------------------//
+             
+             if (vol7 < 24)
                {
-                 mtrackno = igatr;
-                 mtrackpid = kGamma;
+                 smnumber = vol7;
                }
              else
                {
-                 mtrackno  = ichtr;
-                 mtrackpid = id_mo;
+                 smnumber = vol7 - 24;
                }
-             
-             xPos = pmdHit->X();
-             yPos = pmdHit->Y();
-             zPos = pmdHit->Z();
-             cellnumber = pmdHit->fVolume[1];
-             smnumber1  = pmdHit->fVolume[4];
-             edep       = pmdHit->fEnergy;
-             
-             if (smnumber1 > 3 && smnumber1 <= 6)
+             Int_t vol8 = smnumber/6 + 1;    // fake supermodule
+
+             if (vol8 == 1 || vol8 == 2)
                {
-                 Int_t ny = (cellnumber-1)/48 + 1;
-                 Int_t nx = cellnumber - (ny-1)*48;
-                 Int_t nx1 = nx + 24;
-                 Int_t ny1 = ny + 48;
-                 Int_t ncell1 = (ny1 - 1)* 72 + nx1;
-                 cellnumber = ncell1;
+                 xpad = vol2;
+                 ypad = vol1;
                }
-             
-             smnumber = Convert2RealSMNumber(smnumber1);
-             ypad = (cellnumber - 1)/fNCell + 1;
-             xpad = cellnumber - (ypad-1)*fNCell;
-
-             //cout << "-zpos = " << -zPos << endl;
-             if (-zPos < fZPos)
+             else if (vol8 == 3 || vol8 == 4)
                {
-                 // CPV
-                 fDetNo = 1;
+                 xpad = vol1;
+                 ypad = vol2;
                }
-             else if (-zPos > fZPos)
+
+             AliDebug(2,Form("ZPosition = %f Edeposition = %f",zPos,edep));
+
+             if (vol7 < 24)
                {
-                 // PMD
+                 // PRE
                  fDetNo = 0;
                }
-             
-             if(fDetNo == 1)
+             else
                {
-                 fCPV[smnumber-1][xpad-1][ypad-1] += edep;
+                 fDetNo = 1;
                }
-             else if (fDetNo == 0)
+
+             Int_t smn = smnumber;
+             Int_t ixx = xpad     - 1;
+             Int_t iyy = ypad     - 1;
+             if (fDetNo == 0)
                {
-                 fPMD[smnumber-1][xpad-1][ypad-1] += edep;
-                 fPMDCounter[smnumber-1][xpad-1][ypad-1]++;
+                 fPRE[smn][ixx][iyy] += edep;
+                 fPRECounter[smn][ixx][iyy]++;
+
+                 AliPMDcell* cell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
+
+                 fCell.Add(cell);
+               }
+             else if(fDetNo == 1)
+               {
+                 fCPV[smn][ixx][iyy] += edep;
+                 fCPVCounter[smn][ixx][iyy]++;
+                 AliPMDcell* cpvcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep); 
+                 fCPVCell.Add(cpvcell);
                }
            }
        }
     } // Track Loop ended
-
+  TrackAssignment2CPVCell();
   TrackAssignment2Cell();
   ResetCell();
 
+  Float_t gain1  = 1.;
+  Float_t adc    = 0. ;
   Float_t deltaE = 0.;
-  Int_t detno = 0;
-  Int_t count_digit = 0;
-  Int_t trno = 1;
-  Int_t cellno;
+  Int_t detno    = 0;
+  Int_t trno     = 1;
+  Int_t trpid    = -99;
 
   for (Int_t idet = 0; idet < 2; idet++)
-    {
-      for (Int_t ism = 0; ism < fTotSM; ism++)
-       {
-         for (Int_t jrow = 0; jrow < fNCell; jrow++)
-           {
-             for (Int_t kcol = 0; kcol < fNCell; kcol++)
-               {
-                 cellno = jrow + kcol*fNCell;
+  {
+      for (Int_t ism = 0; ism < fgkTotUM; ism++)
+      {
+         for (Int_t jrow = 0; jrow < fgkRow; jrow++)
+         {
+             for (Int_t kcol = 0; kcol < fgkCol; kcol++)
+             {
                  if (idet == 0)
-                   {
-                     deltaE = fPMD[ism][jrow][kcol];
-                     detno = 0;
-                   }
+                 {
+                     deltaE = fPRE[ism][jrow][kcol];
+                     trno   = fPRETrackNo[ism][jrow][kcol];
+                     detno  = 0;
+                 }
                  else if (idet == 1)
-                   {
+                 {
                      deltaE = fCPV[ism][jrow][kcol];
-                     detno = 1;
-                   }
+                     trno   = fCPVTrackNo[ism][jrow][kcol];
+                     detno  = 1;
+                 }
                  if (deltaE > 0.)
-                   {
-                     count_digit++;
-                     AddDigit(trno,detno,ism,cellno,deltaE);
-                   }
-               } // column loop
-           } // row    loop
-       } // supermodule loop
-      treeD->Fill();
-      ResetDigit();
-    } // detector loop
-
-  pmdloader->WriteDigits("OVERWRITE");
-
-  ResetCellADC();
+                 {
+                     MeV2ADC(deltaE,adc);
+
+                     // To decalibrate the adc values
+                     //
+                     gain1 = Gain(idet,ism,jrow,kcol);
+                     if (gain1 != 0.)
+                     {
+                         Int_t adcDecalib = (Int_t)(adc/gain1);
+                         adc = (Float_t) adcDecalib;
+                     }
+                     else if(gain1 == 0.)
+                     {
+                         adc = 0.;
+                     }
+
+                     // Pedestal Decalibration
+                     Int_t   pedmeanrms = 
+                         fCalibPed->GetPedMeanRms(idet,ism,jrow,kcol);
+                     Int_t   pedrms1    = (Int_t) pedmeanrms%100;
+                     Float_t pedrms     = (Float_t)pedrms1/10.;
+                     Float_t pedmean    = 
+                         (Float_t) (pedmeanrms - pedrms1)/1000.0;
+                     if (adc > 0.)
+                     {
+                         adc += (pedmean + 3.0*pedrms);
+                         TParticle *mparticle
+                           = gAlice->GetMCApp()->Particle(trno);
+                         trpid = mparticle->GetPdgCode();
+                         
+                         AddDigit(trno,trpid,detno,ism,jrow,kcol,adc);
+                     }
+                 }
+             } // column loop
+         } // row    loop
+         treeD->Fill();
+         ResetDigit();
+      } // supermodule loop
+  } // detector loop
   
-  //  cout << " -------- End of Hits2Digit ----------- " << endl;
+  fPMDLoader->WriteDigits("OVERWRITE");
+  ResetCellADC();
+
 }
+//____________________________________________________________________________
 
 
 void AliPMDDigitizer::SDigits2Digits(Int_t ievt)
 {
-  //  cout << " -------- Beginning of SDigits2Digit ----------- " << endl;
+  // This reads the PMD sdigits tree and converts energy deposition
+  // in a cell to ADC and stores in the digits tree
+  //
+
   fRunLoader->GetEvent(ievt);
 
-  treeS = pmdloader->TreeS();
+  TTree* treeS = fPMDLoader->TreeS();
   AliPMDsdigit  *pmdsdigit;
   TBranch *branch = treeS->GetBranch("PMDSDigit");
+  if(!branch)
+    {
+      AliError("PMD Sdigit branch does not exist");
+      return;
+    }
+  if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
   branch->SetAddress(&fSDigits);
 
-  treeD = pmdloader->TreeD();
+  TTree* treeD = fPMDLoader->TreeD();
   if (treeD == 0x0)
     {
-      pmdloader->MakeTree("D");
-      treeD = pmdloader->TreeD();
+      fPMDLoader->MakeTree("D");
+      treeD = fPMDLoader->TreeD();
     }
   Int_t bufsize = 16000;
-  treeD->Branch("PMDDigit", &fDigits, bufsize); 
+  if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
+  treeD->Branch("PMDDigit", &fDigits, bufsize);
 
-  Int_t   trno, det, smn;
-  Int_t   cellno;
-  Float_t edep, adc;
+  Int_t   trno = 1, trpid = 0, det = 0, smn = 0;
+  Int_t   irow = 0, icol = 0;
+  Float_t edep = 0., adc = 0.;
 
   Int_t nmodules = (Int_t) treeS->GetEntries();
+  AliDebug(1,Form("Number of modules = %d",nmodules));
 
   for (Int_t imodule = 0; imodule < nmodules; imodule++)
     {
-      treeS->GetEntry(imodule); 
+      treeS->GetEntry(imodule);
       Int_t nentries = fSDigits->GetLast();
-      //cout << " nentries = " << nentries << endl;
+      AliDebug(2,Form("Number of entries per module = %d",nentries+1));
       for (Int_t ient = 0; ient < nentries+1; ient++)
        {
          pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
          trno   = pmdsdigit->GetTrackNumber();
+         trpid  = pmdsdigit->GetTrackPid();
          det    = pmdsdigit->GetDetector();
          smn    = pmdsdigit->GetSMNumber();
-         cellno = pmdsdigit->GetCellNumber();
+         irow   = pmdsdigit->GetRow();
+         icol   = pmdsdigit->GetColumn();
          edep   = pmdsdigit->GetCellEdep();
 
          MeV2ADC(edep,adc);
 
-         AddDigit(trno,det,smn,cellno,adc);      
+         // To decalibrte the adc values
+         //
+         Float_t gain1 = Gain(det,smn,irow,icol);
+         if (gain1 != 0.)
+         {
+           Int_t adcDecalib = (Int_t)(adc/gain1);
+           adc = (Float_t) adcDecalib;
+         }
+         else if(gain1 == 0.)
+         {
+             adc = 0.;
+         }
+         // Pedestal Decalibration
+         Int_t   pedmeanrms = fCalibPed->GetPedMeanRms(det,smn,irow,icol);
+         Int_t   pedrms1    = (Int_t) pedmeanrms%100;
+         Float_t pedrms     = (Float_t)pedrms1/10.;
+         Float_t pedmean    = (Float_t) (pedmeanrms - pedrms1)/1000.0;
+         if(adc > 0.)
+         {
+             adc += (pedmean + 3.0*pedrms);
+             AddDigit(trno,trpid,det,smn,irow,icol,adc);
+         }
+
        }
       treeD->Fill();
       ResetDigit();
     }
-  pmdloader->WriteDigits("OVERWRITE");
-  //  cout << " -------- End of SDigits2Digit ----------- " << endl;
+  fPMDLoader->WriteDigits("OVERWRITE");
+
+}
+//____________________________________________________________________________
+void AliPMDDigitizer::Digitize(Option_t *option)
+{
+  // Does the event merging and digitization
+  const char *cdeb = strstr(option,"deb");
+  if(cdeb)
+    {
+      AliDebug(100," *** PMD Exec is called ***");
+    }
+
+  Int_t ninputs = fDigInput->GetNinputs();
+  AliDebug(1,Form("Number of files to be processed = %d",ninputs));
+  ResetCellADC();
+
+  for (Int_t i = 0; i < ninputs; i++)
+    {
+      Int_t troffset = fDigInput->GetMask(i);
+      MergeSDigits(i, troffset);
+    }
+
+  fRunLoader = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
+  fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
+  fPMDLoader = fRunLoader->GetLoader("PMDLoader");
+  if (fPMDLoader == 0x0)
+    {
+      AliError("Can not find PMD or PMDLoader");
+    }
+  fPMDLoader->LoadDigits("update");
+  TTree* treeD = fPMDLoader->TreeD();
+  if (treeD == 0x0)
+    {
+      fPMDLoader->MakeTree("D");
+      treeD = fPMDLoader->TreeD();
+    }
+  Int_t bufsize = 16000;
+  if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
+  treeD->Branch("PMDDigit", &fDigits, bufsize);
+
+  Float_t adc    = 0.;
+  Float_t deltaE = 0.;
+  Int_t detno    = 0;
+  Int_t trno     = 1;
+  Int_t trpid    = -99;
+
+  for (Int_t idet = 0; idet < 2; idet++)
+    {
+      for (Int_t ism = 0; ism < fgkTotUM; ism++)
+       {
+         for (Int_t jrow = 0; jrow < fgkRow; jrow++)
+           {
+             for (Int_t kcol = 0; kcol < fgkCol; kcol++)
+               {
+                 if (idet == 0)
+                   {
+                     deltaE = fPRE[ism][jrow][kcol];
+                     trno   = fPRETrackNo[ism][jrow][kcol];
+                     trpid  = fPRETrackPid[ism][jrow][kcol];
+                     detno  = 0;
+                   }
+                 else if (idet == 1)
+                   {
+                     deltaE = fCPV[ism][jrow][kcol];
+                     trno   = fCPVTrackNo[ism][jrow][kcol];
+                     trpid  = fCPVTrackPid[ism][jrow][kcol];
+                     detno  = 1;
+                   }
+                 if (deltaE > 0.)
+                   {
+                     MeV2ADC(deltaE,adc);
+
+                      //
+                     // Gain decalibration
+                     //
+                     Float_t gain1 = Gain(idet,ism,jrow,kcol);
+
+                     if (gain1 != 0.)
+                     {
+                         Int_t adcDecalib = (Int_t)(adc/gain1);
+                         adc = (Float_t) adcDecalib;
+                     }
+                     else if(gain1 == 0.)
+                     {
+                         adc = 0.;
+                     }
+                     // Pedestal Decalibration
+                     Int_t   pedmeanrms = 
+                         fCalibPed->GetPedMeanRms(idet,ism,jrow,kcol);
+                     Int_t   pedrms1    = (Int_t) pedmeanrms%100;
+                     Float_t pedrms     = (Float_t)pedrms1/10.;
+                     Float_t pedmean    = 
+                         (Float_t) (pedmeanrms - pedrms1)/1000.0;
+                     if (adc > 0.)
+                     {
+                         adc += (pedmean + 3.0*pedrms);
+                         AddDigit(trno,trpid,detno,ism,jrow,kcol,adc);
+                     }
+
+                   }
+               } // column loop
+           } // row    loop
+         treeD->Fill();
+         ResetDigit();
+       } // supermodule loop
+    } // detector loop
+  fPMDLoader->WriteDigits("OVERWRITE");
+  fPMDLoader->UnloadDigits();
+  ResetCellADC();
+}
+//____________________________________________________________________________
+void AliPMDDigitizer::TrackAssignment2CPVCell()
+{
+  // This block assigns the cell id when there are
+  // multiple tracks in a cell according to the
+  // energy deposition
+  // This method added by Ajay
+  Bool_t jsort = false;
+
+  Int_t i = 0, j = 0, k = 0;
+
+  Int_t   *status1;
+  Int_t   *status2;
+  Int_t   *trnarray;  
+  Float_t *fracEdp;
+  Float_t *trEdp;
+  
+  Int_t   ****cpvTrack;
+  Float_t ****cpvEdep;
+
+  cpvTrack = new Int_t ***[fgkTotUM];
+  cpvEdep  = new Float_t ***[fgkTotUM];
+  for (i=0; i<fgkTotUM; i++)
+    {
+      cpvTrack[i] = new Int_t **[fgkRow];
+      cpvEdep[i]  = new Float_t **[fgkRow];
+    }
+
+  for (i = 0; i < fgkTotUM; i++)
+    {
+      for (j = 0; j < fgkRow; j++)
+       {
+         cpvTrack[i][j] = new Int_t *[fgkCol];
+         cpvEdep[i][j]  = new Float_t *[fgkCol];
+       }
+    }
+  for (i = 0; i < fgkTotUM; i++)
+    {
+      for (j = 0; j < fgkRow; j++)
+       {
+         for (k = 0; k < fgkCol; k++)
+           {
+             Int_t nn = fCPVCounter[i][j][k];
+             if(nn > 0)
+               {
+                 cpvTrack[i][j][k] = new Int_t[nn];
+                 cpvEdep[i][j][k] = new Float_t[nn];
+               }
+             else
+               {
+                 nn = 1;
+                 cpvTrack[i][j][k] = new Int_t[nn];
+                 cpvEdep[i][j][k] = new Float_t[nn];
+               }                     
+             fCPVCounter[i][j][k] = 0;
+           }
+       }
+    }
+
+
+  Int_t nentries = fCPVCell.GetEntries();
+  Int_t   mtrackno = 0, ism = 0, ixp = 0, iyp = 0;
+  Float_t edep = 0.;
+  for (i = 0; i < nentries; i++)
+    {
+      AliPMDcell* cpvcell = (AliPMDcell*)fCPVCell.UncheckedAt(i);
+      
+      mtrackno = cpvcell->GetTrackNumber();
+      ism      = cpvcell->GetSMNumber();
+      ixp      = cpvcell->GetX();
+      iyp      = cpvcell->GetY();
+      edep     = cpvcell->GetEdep();
+      Int_t nn = fCPVCounter[ism][ixp][iyp];
+      cpvTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
+      cpvEdep[ism][ixp][iyp][nn] = edep;
+      fCPVCounter[ism][ixp][iyp]++;
+    }
+  
+  Int_t iz = 0, il = 0;
+  Int_t im = 0, ix = 0, iy = 0;
+  Int_t nn = 0;
+  for (im=0; im<fgkTotUM; im++)
+    {
+      for (ix=0; ix<fgkRow; ix++)
+       {
+         for (iy=0; iy<fgkCol; iy++)
+           {
+             nn = fCPVCounter[im][ix][iy];
+             if (nn > 1)
+               {
+                 // This block handles if a cell is fired
+                 // many times by many tracks
+                 status1  = new Int_t[nn];
+                 status2  = new Int_t[2*nn];
+                 trnarray = new Int_t[nn];
+                 for (iz = 0; iz < nn; iz++)
+                   {
+                     status1[iz] = cpvTrack[im][ix][iy][iz];
+                   }
+                 TMath::Sort(nn,status1,status2,jsort);
+                 Int_t trackOld = -99999;
+                 Int_t track, trCount = 0;
+                 for (iz = 0; iz < nn; iz++)
+                   {
+                     track = status1[status2[iz]];
+                     if (trackOld != track)
+                       {
+                         trnarray[trCount] = track;
+                         trCount++;
+                       }                             
+                     trackOld = track;
+                   }
+                 delete [] status1;
+                 delete [] status2;
+                 Float_t totEdp = 0.;
+                 trEdp = new Float_t[trCount];
+                 fracEdp = new Float_t[trCount];
+                 for (il = 0; il < trCount; il++)
+                   {
+                     trEdp[il] = 0.;
+                     track = trnarray[il];
+                     for (iz = 0; iz < nn; iz++)
+                       {
+                         if (track == cpvTrack[im][ix][iy][iz])
+                           {
+                             trEdp[il] += cpvEdep[im][ix][iy][iz];
+                           }
+                       }
+                     totEdp += trEdp[il];
+                   }
+                 Int_t ilOld = 0;
+                 Float_t fracOld = 0.;
+                 
+                 for (il = 0; il < trCount; il++)
+                   {
+                     fracEdp[il] = trEdp[il]/totEdp;
+                     if (fracOld < fracEdp[il])
+                       {
+                         fracOld = fracEdp[il];
+                         ilOld = il;
+                       }
+                   }
+                 fCPVTrackNo[im][ix][iy] = trnarray[ilOld];
+                 delete [] fracEdp;
+                 delete [] trEdp;
+                 delete [] trnarray;
+               }
+             else if (nn == 1)
+               {
+                 // This only handles if a cell is fired
+                 // by only one track
+                 
+                 fCPVTrackNo[im][ix][iy] = cpvTrack[im][ix][iy][0];
+                 
+               }
+             else if (nn ==0)
+               {
+                 // This is if no cell is fired
+                 fCPVTrackNo[im][ix][iy] = -999;
+               }
+           } // end of iy
+       } // end of ix
+    } // end of im
+  
+  // Delete all the pointers
+  
+ for (i = 0; i < fgkTotUM; i++)
+    {
+      for (j = 0; j < fgkRow; j++)
+       {
+         for (k = 0; k < fgkCol; k++)
+           {
+             delete []cpvTrack[i][j][k];
+             delete []cpvEdep[i][j][k];
+           }
+       }
+    }
+  for (i = 0; i < fgkTotUM; i++)
+    {
+      for (j = 0; j < fgkRow; j++)
+       {
+         delete [] cpvTrack[i][j];
+         delete [] cpvEdep[i][j];
+       }
+    }
+  
+  for (i = 0; i < fgkTotUM; i++)
+    {
+      delete [] cpvTrack[i];
+      delete [] cpvEdep[i];
+    }
+  delete [] cpvTrack;
+  delete [] cpvEdep;
+  
+  // 
+  // End of the cell id assignment
+  //
 }
+//____________________________________________________________________________
+
+void AliPMDDigitizer::MergeSDigits(Int_t filenumber, Int_t troffset)
+{
+  // merging sdigits
+  fRunLoader = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(filenumber));
+  fPMDLoader = fRunLoader->GetLoader("PMDLoader");
+  fPMDLoader->LoadSDigits("read");
+  TTree* treeS = fPMDLoader->TreeS();
+  AliPMDsdigit  *pmdsdigit;
+  TBranch *branch = treeS->GetBranch("PMDSDigit");
+  if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
+  branch->SetAddress(&fSDigits);
 
+  Int_t   itrackno = 1, itrackpid = 0, idet = 0, ism = 0;
+  Int_t   ixp = 0, iyp = 0;
+  Float_t edep = 0.;
+  Int_t nmodules = (Int_t) treeS->GetEntries();
+  AliDebug(1,Form("Number of Modules in the treeS = %d",nmodules));
+  AliDebug(1,Form("Track Offset = %d",troffset));
+  for (Int_t imodule = 0; imodule < nmodules; imodule++)
+    {
+      treeS->GetEntry(imodule);
+      Int_t nentries = fSDigits->GetLast();
+      AliDebug(2,Form("Number of Entries per Module = %d",nentries));
+      for (Int_t ient = 0; ient < nentries+1; ient++)
+       {
+         pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
+         itrackno  = pmdsdigit->GetTrackNumber();
+         itrackpid = pmdsdigit->GetTrackPid();
+         idet      = pmdsdigit->GetDetector();
+         ism       = pmdsdigit->GetSMNumber();
+         ixp       = pmdsdigit->GetRow();
+         iyp       = pmdsdigit->GetColumn();
+         edep      = pmdsdigit->GetCellEdep();
+         if (idet == 0)
+           {
+             if (fPRE[ism][ixp][iyp] < edep)
+               {
+                 fPRETrackNo[ism][ixp][iyp] = troffset + itrackno;
+                 fPRETrackPid[ism][ixp][iyp] = itrackpid;
+               }
+             fPRE[ism][ixp][iyp] += edep;
+           }
+         else if (idet == 1)
+           {
+             if (fCPV[ism][ixp][iyp] < edep)
+               {
+                 fCPVTrackNo[ism][ixp][iyp] = troffset + itrackno;
+                 fCPVTrackPid[ism][ixp][iyp] = itrackpid;
+               }
+             fCPV[ism][ixp][iyp] += edep;
+           }
+       }
+    }
+
+}
+// ----------------------------------------------------------------------
 void AliPMDDigitizer::TrackAssignment2Cell()
 {
-  // To be checked again
   // 
-  // This blocks assign the cell id when there is 
+  // This block assigns the cell id when there are
   // multiple tracks in a cell according to the
   // energy deposition
   //
+  Bool_t jsort = false;
 
-  Int_t i, j, k;
+  Int_t i = 0, j = 0, k = 0;
 
-  Float_t *frac_edp;
-  Float_t *tr_edp;
-  Int_t *status;
-  Int_t   ****PMDTrack;
-  Float_t ****PMDEdep;
+  Int_t   *status1;
+  Int_t   *status2;
+  Int_t   *trnarray;
+  Float_t *fracEdp;
+  Float_t *trEdp;
+  
+  Int_t   ****pmdTrack;
+  Float_t ****pmdEdep;
 
-  PMDTrack = new Int_t ***[27];
-  PMDEdep  = new Float_t ***[27];
-  for (i=0; i<fTotSM; i++)
+  pmdTrack = new Int_t ***[fgkTotUM];
+  pmdEdep  = new Float_t ***[fgkTotUM];
+  for (i=0; i<fgkTotUM; i++)
     {
-      PMDTrack[i] = new Int_t **[72];
-      PMDEdep[i]  = new Float_t **[72];
+      pmdTrack[i] = new Int_t **[fgkRow];
+      pmdEdep[i]  = new Float_t **[fgkRow];
     }
 
-  for (i = 0; i < fTotSM; i++)
+  for (i = 0; i < fgkTotUM; i++)
     {
-      for (j = 0; j < fNCell; j++)
+      for (j = 0; j < fgkRow; j++)
        {
-         PMDTrack[i][j] = new Int_t *[72];
-         PMDEdep[i][j]  = new Float_t *[72];
+         pmdTrack[i][j] = new Int_t *[fgkCol];
+         pmdEdep[i][j]  = new Float_t *[fgkCol];
        }
     }
   
-  for (i = 0; i < fTotSM; i++)
+  for (i = 0; i < fgkTotUM; i++)
     {
-      for (j = 0; j < fNCell; j++)
+      for (j = 0; j < fgkRow; j++)
        {
-         for (k = 0; k < fNCell; k++)
+         for (k = 0; k < fgkCol; k++)
            {
-             Int_t nn = fPMDCounter[i][j][k];
+             Int_t nn = fPRECounter[i][j][k];
              if(nn > 0)
                {
-                 PMDTrack[i][j][k] = new Int_t[nn];
-                 PMDEdep[i][j][k] = new Float_t[nn];
+                 pmdTrack[i][j][k] = new Int_t[nn];
+                 pmdEdep[i][j][k] = new Float_t[nn];
                }
              else
                {
                  nn = 1;
-                 PMDTrack[i][j][k] = new Int_t[nn];
-                 PMDEdep[i][j][k] = new Float_t[nn];
-               }                     
-             fPMDCounter[i][j][k] = 0;
+                 pmdTrack[i][j][k] = new Int_t[nn];
+                 pmdEdep[i][j][k] = new Float_t[nn];
+               }
+             fPRECounter[i][j][k] = 0;
            }
        }
     }
 
 
-  Int_t nentries = fCell->GetEntries();
+  Int_t nentries = fCell.GetEntries();
 
   Int_t   mtrackno, ism, ixp, iyp;
   Float_t edep;
 
   for (i = 0; i < nentries; i++)
     {
-      pmdcell = (AliPMDcell*)fCell->UncheckedAt(i);
-      
-      mtrackno = pmdcell->GetTrackNumber();
-      ism = pmdcell->GetSMNumber();
-      ixp = pmdcell->GetX();
-      iyp = pmdcell->GetY();
-      edep = pmdcell->GetEdep();
-
-      Int_t nn = fPMDCounter[ism][ixp][iyp];
-
-      //      cout << " nn = " << nn << endl;
+      AliPMDcell* cell = (AliPMDcell*)fCell.UncheckedAt(i);
       
-      PMDTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
-      PMDEdep[ism][ixp][iyp][nn] = edep;
-      fPMDCounter[ism][ixp][iyp]++;
+      mtrackno = cell->GetTrackNumber();
+      ism      = cell->GetSMNumber();
+      ixp      = cell->GetX();
+      iyp      = cell->GetY();
+      edep     = cell->GetEdep();
+      Int_t nn = fPRECounter[ism][ixp][iyp];
+      pmdTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
+      pmdEdep[ism][ixp][iyp][nn] = edep;
+      fPRECounter[ism][ixp][iyp]++;
     }
   
-  Int_t iz, il;
-  Int_t im, ix, iy;
-  Int_t nn;
+  Int_t iz = 0, il = 0;
+  Int_t im = 0, ix = 0, iy = 0;
+  Int_t nn = 0;
   
-  for (im=0; im<27; im++)
+  for (im=0; im<fgkTotUM; im++)
     {
-      for (ix=0; ix<72; ix++)
+      for (ix=0; ix<fgkRow; ix++)
        {
-         for (iy=0; iy<72; iy++)
+         for (iy=0; iy<fgkCol; iy++)
            {
-             nn = fPMDCounter[im][ix][iy];
+             nn = fPRECounter[im][ix][iy];
              if (nn > 1)
                {
                  // This block handles if a cell is fired
                  // many times by many tracks
-                 
-                 status = new Int_t[nn];
+                 status1  = new Int_t[nn];
+                 status2  = new Int_t[2*nn];
+                 trnarray = new Int_t[nn];
                  for (iz = 0; iz < nn; iz++)
                    {
-                     status[iz] = PMDTrack[im][ix][iy][iz];
+                     status1[iz] = pmdTrack[im][ix][iy][iz];
                    }
-                 sort(status,status+nn);
-                 Int_t track_old = -99999;
-                 Int_t track, tr_count = 0;
+                 TMath::Sort(nn,status1,status2,jsort);
+                 Int_t trackOld = -99999;
+                 Int_t track, trCount = 0;
                  for (iz = 0; iz < nn; iz++)
                    {
-                     track = status[iz];
-                     if (track_old != track)
+                     track = status1[status2[iz]];
+                     if (trackOld != track)
                        {
-                         tr_count++;
-                         vjunkTRN.push_back(track);
-                       }                             
-                     track_old = track;
+                         trnarray[trCount] = track;
+                         trCount++;
+                       }
+                     trackOld = track;
                    }
-                 delete status;
-                 Float_t tot_edp = 0.;
-                 tr_edp = new Float_t[tr_count];
-                 frac_edp = new Float_t[tr_count];
-                 for (il = 0; il < tr_count; il++)
+                 delete [] status1;
+                 delete [] status2;
+                 Float_t totEdp = 0.;
+                 trEdp = new Float_t[trCount];
+                 fracEdp = new Float_t[trCount];
+                 for (il = 0; il < trCount; il++)
                    {
-                     tr_edp[il] = 0.;
-                     track = vjunkTRN[il];
+                     trEdp[il] = 0.;
+                     track = trnarray[il];
                      for (iz = 0; iz < nn; iz++)
                        {
-                         if (track == PMDTrack[im][ix][iy][iz])
+                         if (track == pmdTrack[im][ix][iy][iz])
                            {
-                             tr_edp[il] += PMDEdep[im][ix][iy][iz];
+                             trEdp[il] += pmdEdep[im][ix][iy][iz];
                            }
                        }
-                     tot_edp += tr_edp[il];
+                     totEdp += trEdp[il];
                    }
-                 vjunkTRN.clear();
-                 Int_t il_old = 0;
-                 Float_t frac_old = 0.;
+                 Int_t ilOld = 0;
+                 Float_t fracOld = 0.;
                  
-                 for (il = 0; il < tr_count; il++)
+                 for (il = 0; il < trCount; il++)
                    {
-                     frac_edp[il] = tr_edp[il]/tot_edp;
-                     if (frac_old < frac_edp[il])
+                     fracEdp[il] = trEdp[il]/totEdp;
+                     if (fracOld < fracEdp[il])
                        {
-                         frac_old = frac_edp[il];
-                         il_old = il;
+                         fracOld = fracEdp[il];
+                         ilOld = il;
                        }
                    }
-                 delete frac_edp;
-                 delete tr_edp;
-                 
-                 fPMDTrackNo[im][ix][iy] = vjunkTRN[il_old];
+                 fPRETrackNo[im][ix][iy] = trnarray[ilOld];
+                 delete [] fracEdp;
+                 delete [] trEdp;
+                 delete [] trnarray;
                }
              else if (nn == 1)
                {
                  // This only handles if a cell is fired
                  // by only one track
                  
-                 fPMDTrackNo[im][ix][iy] = PMDTrack[im][ix][iy][0];
+                 fPRETrackNo[im][ix][iy] = pmdTrack[im][ix][iy][0];
                  
                }
              else if (nn ==0)
                {
                  // This is if no cell is fired
-                 fPMDTrackNo[im][ix][iy] = -999;
+                 fPRETrackNo[im][ix][iy] = -999;
                }
            } // end of iy
        } // end of ix
@@ -766,154 +1398,170 @@ void AliPMDDigitizer::TrackAssignment2Cell()
   
   // Delete all the pointers
   
-  for (i = 0; i < fTotSM; i++)
+  for (i = 0; i < fgkTotUM; i++)
     {
-      for (j = 0; j < fNCell; j++)
+      for (j = 0; j < fgkRow; j++)
        {
-         for (k = 0; k < fNCell; k++)
+         for (k = 0; k < fgkCol; k++)
            {
-             delete [] PMDTrack[i][j][k];
-             delete [] PMDEdep[i][j][k];
+             delete [] pmdTrack[i][j][k];
+             delete [] pmdEdep[i][j][k];
            }
        }
     }
   
-  for (i = 0; i < fTotSM; i++)
+  for (i = 0; i < fgkTotUM; i++)
     {
-      for (j = 0; j < fNCell; j++)
+      for (j = 0; j < fgkRow; j++)
        {
-         delete [] PMDTrack[i][j];
-         delete [] PMDEdep[i][j];
+         delete [] pmdTrack[i][j];
+         delete [] pmdEdep[i][j];
        }
     }
   
-  for (i = 0; i < fTotSM; i++)
+  for (i = 0; i < fgkTotUM; i++)
     {
-      delete [] PMDTrack[i];
-      delete [] PMDEdep[i];
+      delete [] pmdTrack[i];
+      delete [] pmdEdep[i];
     }
-  delete PMDTrack;
-  delete PMDEdep;
+  delete [] pmdTrack;
+  delete [] pmdEdep;
   // 
   // End of the cell id assignment
   //
 }
+//____________________________________________________________________________
+void AliPMDDigitizer::MeV2ADC(Float_t mev, Float_t & adc) const
+{
+  // This converts the simulated edep to ADC according to the
+  // Test Beam Data
+  // PS Test in June 2010, Voltage @ 1300 V
+  // KeV - ADC conversion for 12bit ADC
+  // MPV data used for the fit and taken here
 
+  // constants are from Test Beam 2010
+  
+  const Float_t kConstant   = 0.612796;
+  const Float_t kSlope      = 130.158;
+  
+  Float_t adc12bit = kSlope*mev*0.001 + kConstant;
+  if (adc12bit < 0.) adc12bit = 0.;
 
-void AliPMDDigitizer::MeV2ADC(Float_t mev, Float_t & adc)
-{
-  // To be done
+  //Introducing Readout Resolution for ALICE-PMD
+
+  Float_t sigrr     = 0.605016 - 0.000273*adc12bit + 6.54e-8*adc12bit*adc12bit;
+  Float_t adcwithrr = gRandom->Gaus(adc12bit,sigrr);
+
+  if(adcwithrr < 0.)
+    {
+      adc = 0.;
+    }
+  else if(adcwithrr >= 0. && adcwithrr < 1600.0)
+    {
+      adc = adcwithrr;
+    }
+  else if (adcwithrr >= 1600.0)
+    {
+      adc = 1600.0;
+    }
 
-  adc = mev*1.;
 }
-void AliPMDDigitizer::AddSDigit(Int_t trnumber, Int_t det, Int_t smnumber, 
-  Int_t cellnumber, Float_t adc)
+//____________________________________________________________________________
+void AliPMDDigitizer::AddSDigit(Int_t trnumber, Int_t trpid, Int_t det,
+                               Int_t smnumber, Int_t irow, Int_t icol,
+                               Float_t adc)
 {
+  // Add SDigit
+  //
+  if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
   TClonesArray &lsdigits = *fSDigits;
-  AliPMDsdigit *newcell;
-  newcell = new AliPMDsdigit(trnumber,det,smnumber,cellnumber,adc);
-  new(lsdigits[fNsdigit++]) AliPMDsdigit(newcell);
-  delete newcell;
+  new(lsdigits[fNsdigit++])  AliPMDsdigit(trnumber,trpid,det,smnumber,irow,icol,adc);
 }
+//____________________________________________________________________________
 
-void AliPMDDigitizer::AddDigit(Int_t trnumber, Int_t det, Int_t smnumber, 
-  Int_t cellnumber, Float_t adc)
+void AliPMDDigitizer::AddDigit(Int_t trnumber, Int_t trpid, Int_t det,
+                              Int_t smnumber, Int_t irow, Int_t icol,
+                              Float_t adc)
 {
+  // Add Digit
+  //
+  if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
   TClonesArray &ldigits = *fDigits;
-  AliPMDdigit *newcell;
-  newcell = new AliPMDdigit(trnumber,det,smnumber,cellnumber,adc);
-  new(ldigits[fNdigit++]) AliPMDdigit(newcell);
-  delete newcell;
+  new(ldigits[fNdigit++]) AliPMDdigit(trnumber,trpid, det,smnumber,irow,icol,adc);
 }
+//____________________________________________________________________________
 
-Int_t AliPMDDigitizer::Convert2RealSMNumber(Int_t smnumber1)
-{
-  Int_t smnumber = -999;
-
-  if (smnumber1==1)  smnumber =  1;
-  if (smnumber1==2)  smnumber = 10;
-  if (smnumber1==3)  smnumber = 19;
-  if (smnumber1==4)  smnumber =  1;
-  if (smnumber1==5)  smnumber = 10;
-  if (smnumber1==6)  smnumber = 19;
-  if (smnumber1==7)  smnumber =  2;
-  if (smnumber1==8)  smnumber =  3;
-  if (smnumber1==9)  smnumber =  4;
-  if (smnumber1==10) smnumber =  5;
-  if (smnumber1==11) smnumber =  6;
-  if (smnumber1==12) smnumber =  7;
-  if (smnumber1==13) smnumber =  8;
-  if (smnumber1==14) smnumber =  9;
-  if (smnumber1==15) smnumber = 11;
-  if (smnumber1==16) smnumber = 12;
-  if (smnumber1==17) smnumber = 13;
-  if (smnumber1==18) smnumber = 14;
-  if (smnumber1==19) smnumber = 15;
-  if (smnumber1==20) smnumber = 16;
-  if (smnumber1==21) smnumber = 17;
-  if (smnumber1==22) smnumber = 18;
-  if (smnumber1==23) smnumber = 20;
-  if (smnumber1==24) smnumber = 21;
-  if (smnumber1==25) smnumber = 22;
-  if (smnumber1==26) smnumber = 23;
-  if (smnumber1==27) smnumber = 24;
-  if (smnumber1==28) smnumber = 25;
-  if (smnumber1==29) smnumber = 26;
-  if (smnumber1==30) smnumber = 27;
-
-  return smnumber;
-}
 void AliPMDDigitizer::SetZPosition(Float_t zpos)
 {
   fZPos = zpos;
 }
+//____________________________________________________________________________
 Float_t AliPMDDigitizer::GetZPosition() const
 {
   return fZPos;
 }
+//____________________________________________________________________________
 
 void AliPMDDigitizer::ResetCell()
 {
-  fCell->Clear();
-  for (Int_t i = 0; i < fTotSM; i++)
+  // clears the cell array and also the counter
+  //  for each cell
+  //
+  fCPVCell.Delete();
+  fCell.Delete();
+  for (Int_t i = 0; i < fgkTotUM; i++)
     {
-      for (Int_t j = 0; j < fNCell; j++)
+      for (Int_t j = 0; j < fgkRow; j++)
        {
-         for (Int_t k = 0; k < fNCell; k++)
+         for (Int_t k = 0; k < fgkCol; k++)
            {
-             fPMDCounter[i][j][k] = 0; 
+             fCPVCounter[i][j][k] = 0; 
+             fPRECounter[i][j][k] = 0;
            }
        }
     }
 }
+//____________________________________________________________________________
 void AliPMDDigitizer::ResetSDigit()
 {
+  // Clears SDigits
   fNsdigit = 0;
-  if (fSDigits) fSDigits->Clear();
+  if (fSDigits) fSDigits->Delete();
 }
+//____________________________________________________________________________
 void AliPMDDigitizer::ResetDigit()
 {
+  // Clears Digits
   fNdigit = 0;
-  if (fDigits) fDigits->Clear();
+  if (fDigits) fDigits->Delete();
 }
+//____________________________________________________________________________
 
 void AliPMDDigitizer::ResetCellADC()
 {
-  for (Int_t i = 0; i < fTotSM; i++)
+  // Clears individual cells edep and track number
+  for (Int_t i = 0; i < fgkTotUM; i++)
     {
-      for (Int_t j = 0; j < fNCell; j++)
+      for (Int_t j = 0; j < fgkRow; j++)
        {
-         for (Int_t k = 0; k < fNCell; k++)
+         for (Int_t k = 0; k < fgkCol; k++)
            {
-             fCPV[i][j][k] = 0.; 
-             fPMD[i][j][k] = 0.; 
+             fCPV[i][j][k]         = 0.;
+             fPRE[i][j][k]         = 0.;
+             fCPVTrackNo[i][j][k]  = 0;
+             fPRETrackNo[i][j][k]  = 0;
+             fCPVTrackPid[i][j][k] = -1;
+             fPRETrackPid[i][j][k] = -1;
            }
        }
     }
 }
+//____________________________________________________________________________
 
 void AliPMDDigitizer::UnLoad(Option_t *option)
 {
+  // Unloads all the root files
+  //
   const char *cS = strstr(option,"S");
   const char *cD = strstr(option,"D");
 
@@ -923,11 +1571,65 @@ void AliPMDDigitizer::UnLoad(Option_t *option)
 
   if (cS)
     {
-      pmdloader->UnloadHits();
+      fPMDLoader->UnloadHits();
     }
   if (cD)
     {
-      pmdloader->UnloadHits();
-      pmdloader->UnloadSDigits();
+      fPMDLoader->UnloadHits();
+      fPMDLoader->UnloadSDigits();
     }
 }
+
+//----------------------------------------------------------------------
+Float_t AliPMDDigitizer::Gain(Int_t det, Int_t smn, Int_t row, Int_t col) const
+{
+  // returns of the gain of the cell
+  // Added this method by ZA
+
+  //cout<<" I am here in gain "<<fCalibData<< "smn,row, col "<<smn
+  //<<" "<<row<<" "<<col<<endl;
+
+  if(!fCalibGain) {
+    AliError("No calibration data loaded from CDB!!!");
+    return 1;
+  }
+
+  Float_t GainFact;
+  GainFact = fCalibGain->GetGainFact(det,smn,row,col);
+  return GainFact;
+}
+//----------------------------------------------------------------------
+AliPMDCalibData* AliPMDDigitizer::GetCalibGain() const
+{
+  // The run number will be centralized in AliCDBManager,
+  // you don't need to set it here!
+  // Added this method by ZA
+  // Cleaned up by Alberto
+  AliCDBEntry  *entry = AliCDBManager::Instance()->Get("PMD/Calib/Gain");
+  
+  if(!entry) AliFatal("Calibration object retrieval failed!");
+  
+  AliPMDCalibData *calibdata=0;
+  if (entry) calibdata = (AliPMDCalibData*) entry->GetObject();
+  
+  if (!calibdata)  AliFatal("No calibration data from calibration database !");
+  
+  return calibdata;
+}
+//----------------------------------------------------------------------
+AliPMDPedestal* AliPMDDigitizer::GetCalibPed() const
+{
+  // The run number will be centralized in AliCDBManager,
+  // you don't need to set it here!
+
+  AliCDBEntry  *entry = AliCDBManager::Instance()->Get("PMD/Calib/Ped");
+  
+  if(!entry) AliFatal("Pedestal object retrieval failed!");
+  
+  AliPMDPedestal *pedestal=0;
+  if (entry) pedestal = (AliPMDPedestal*) entry->GetObject();
+  
+  if (!pedestal)  AliFatal("No pedestal data from calibration database !");
+  
+  return pedestal;
+}