]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
PMD digitization program
authorbnandi <bnandi@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 13 Sep 2003 20:43:01 +0000 (20:43 +0000)
committerbnandi <bnandi@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 13 Sep 2003 20:43:01 +0000 (20:43 +0000)
PMD/AliPMDDigitizer.cxx [new file with mode: 0644]
PMD/AliPMDDigitizer.h [new file with mode: 0644]

diff --git a/PMD/AliPMDDigitizer.cxx b/PMD/AliPMDDigitizer.cxx
new file mode 100644 (file)
index 0000000..153f1e3
--- /dev/null
@@ -0,0 +1,934 @@
+//-----------------------------------------------------//
+//                                                     //
+//  Source File : PMDDigitization.cxx, Version 00      //
+//                                                     //
+//  Date   : September 20 2002                         //
+//                                                     //
+//-----------------------------------------------------//
+
+#include <TBRIK.h>
+#include <TNode.h>
+#include <TTree.h>
+#include <TGeometry.h>
+#include <TObjArray.h>
+#include <TClonesArray.h>
+#include <TFile.h>
+#include <TNtuple.h>
+#include <TParticle.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 "AliHeader.h"
+
+#include "AliPMDcell.h"
+#include "AliPMDsdigit.h"
+#include "AliPMDdigit.h"
+#include "AliPMDDigitizer.h"
+#include "AliPMDClustering.h"
+#include "AliPMDContainer.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++)
+    {
+      for (Int_t j = 0; j < fNCell; j++)
+       {
+         for (Int_t k = 0; k < fNCell; k++)
+           {
+             fCPV[i][j][k] = 0.; 
+             fPMD[i][j][k] = 0.; 
+           }
+       }
+    }
+
+  if (!fCell) fCell = new TObjArray();
+
+  fZPos = 361.5; // in units of cm, This is the default position of PMD
+
+}
+AliPMDDigitizer::~AliPMDDigitizer()
+{
+  delete fSDigits;
+  delete fDigits;
+  delete fCell;
+}
+//
+// Member functions
+//
+void AliPMDDigitizer::OpengAliceFile(Char_t *file, Option_t *option)
+{
+
+  fRunLoader = AliRunLoader::Open(file,AliConfig::fgkDefaultEventFolderName,
+                                 "UPDATE");
+  
+  if (!fRunLoader)
+   {
+     Error("Open","Can not open session for file %s.",file);
+   }
+  
+  fRunLoader->LoadgAlice();
+  fRunLoader->LoadHeader();
+  fRunLoader->LoadKinematics();
+
+  gAlice = fRunLoader->GetAliRun();
+  
+  if (gAlice)
+    {
+      printf("<AliPMDdigitizer::Open> ");
+      printf("AliRun object found on file.\n");
+    }
+  else
+    {
+      printf("<AliPMDdigitizer::Open> ");
+      printf("Could not find AliRun object.\n");
+    }
+
+  PMD  = (AliPMD*)gAlice->GetDetector("PMD");
+  pmdloader = fRunLoader->GetLoader("PMDLoader");
+  if (pmdloader == 0x0)
+    {
+      cerr<<"Hits2Digits : Can not find PMD or PMDLoader\n";
+    }
+
+  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");
+    }
+  else if (cHD)
+    {
+      pmdloader->LoadHits("READ");
+      pmdloader->LoadDigits("recreate");
+    }
+  else if (cSD)
+    {
+      pmdloader->LoadSDigits("READ");
+      pmdloader->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;
+  Float_t vx = -999.0, vy = -999.0, vz = -999.0;
+
+  
+  ResetSDigit();
+
+  printf("Event Number =  %d \n",ievt); 
+  Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
+  printf("Number of Particles = %d \n", nparticles);
+  fRunLoader->GetEvent(ievt);
+  Particles = gAlice->Particles();
+  // ------------------------------------------------------- //
+  // Pointer to specific detector hits.
+  // Get pointers to Alice detectors and Hits containers
+
+  treeH = pmdloader->TreeH();
+  
+  Int_t ntracks    = (Int_t) treeH->GetEntries();
+  printf("Number of Tracks in the TreeH = %d \n", ntracks);
+
+  treeS = pmdloader->TreeS();
+  if (treeS == 0x0)
+    {
+      pmdloader->MakeTree("S");
+      treeS = pmdloader->TreeS();
+    }
+  Int_t bufsize = 16000;
+  treeS->Branch("PMDSDigit", &fSDigits, bufsize); 
+  
+  if (PMD) PMDhits   = PMD->Hits();
+
+  // Start loop on tracks in the hits containers
+
+  
+  for (Int_t track=0; track<ntracks;track++) 
+    {
+      gAlice->ResetHits();
+      treeH->GetEvent(track);
+      
+      if (PMD) 
+       {
+         npmd = PMDhits->GetEntriesFast();
+         for (int ipmd = 0; ipmd < npmd; ipmd++) 
+           {
+             pmdHit = (AliPMDhit*) PMDhits->UncheckedAt(ipmd);
+             trackno = pmdHit->GetTrack();
+
+             //  get kinematics of the particles
+             
+             particle = gAlice->Particle(trackno);
+             trackpid  = particle->GetPdgCode();
+
+             Int_t igatr = -999;
+             Int_t ichtr = -999;
+
+             Int_t imo;
+             Int_t igen = 0;
+             Int_t id_mo = -999;
+             TParticle*  mparticle = particle;
+             
+             while((imo = mparticle->GetFirstMother()) >= 0)
+               {
+                 igen++;
+                 mparticle =  gAlice->Particle(imo);
+                 id_mo = mparticle->GetPdgCode();
+                 
+                 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.)
+                   {
+                     igatr = imo;
+                   }
+                 ichtr = imo;
+               }
+
+             if (id_mo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
+               {
+                 mtrackno = igatr;
+                 mtrackpid = kGamma;
+               }
+             else
+               {
+                 mtrackno  = ichtr;
+                 mtrackpid = id_mo;
+               }
+
+             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 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;
+               }
+             smnumber = Convert2RealSMNumber(smnumber1);
+             ypad = (cellnumber - 1)/fNCell + 1;
+             xpad = cellnumber - (ypad-1)*fNCell;
+
+             //cout << "-zpos = " << -zPos << endl;
+
+             if (-zPos < fZPos)
+               {
+                 // CPV
+                 fDetNo = 1;
+               }
+             else if (-zPos > fZPos)
+               {
+                 // PMD
+                 fDetNo = 0;
+               }
+
+             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;
+
+                 pmdcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
+
+                 fCell->Add(pmdcell);
+               }
+             else if(fDetNo == 1)
+               {
+                 fCPV[smnumber-1][xpad-1][ypad-1] += edep;
+                 fCPVTrackNo[smnumber-1][xpad-1][ypad-1] = mtrackno;
+               }
+           }
+       }
+    } // Track Loop ended
+  
+  
+  TrackAssignment2Cell();
+  
+  ResetCell();
+
+  Float_t deltaE      = 0.;
+  Int_t   detno       = 0;
+  Int_t   count_digit = 0;
+  Int_t   trno        = -1;
+  Int_t   cellno      = 0;
+
+  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;
+                 if (idet == 0)
+                   {
+                     deltaE = fPMD[ism][jrow][kcol];
+                     trno   = fPMDTrackNo[ism][jrow][kcol];
+                     detno = 0;
+                   }
+                 else if (idet == 1)
+                   {
+                     deltaE = fCPV[ism][jrow][kcol];
+                     trno   = fCPVTrackNo[ism][jrow][kcol];
+                     detno = 1;
+                   }
+                 if (deltaE > 0.)
+                   {
+                     count_digit++;
+                     AddSDigit(trno,detno,ism,cellno,deltaE);
+                   }
+               }
+           }
+         treeS->Fill();
+         ResetSDigit();
+       }
+    }
+  pmdloader->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;
+  Float_t vx = -999.0, vy = -999.0, vz = -999.0;
+
+  
+  ResetDigit();
+
+  printf("Event Number =  %d \n",ievt); 
+
+  Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
+  printf("Number of Particles = %d \n", 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");
+
+  if (pmdloader == 0x0)
+    {
+      cerr<<"Hits2Digits method : Can not find PMD or PMDLoader\n";
+    }
+  treeH = pmdloader->TreeH();
+  Int_t ntracks    = (Int_t) treeH->GetEntries();
+  printf("Number of Tracks in the TreeH = %d \n", ntracks);
+  pmdloader->LoadDigits("recreate");
+  treeD = pmdloader->TreeD();
+  if (treeD == 0x0)
+    {
+      pmdloader->MakeTree("D");
+      treeD = pmdloader->TreeD();
+    }
+  Int_t bufsize = 16000;
+  treeD->Branch("PMDDigit", &fDigits, bufsize); 
+  
+  if (PMD) PMDhits   = PMD->Hits();
+
+  // Start loop on tracks in the hits containers
+
+  for (Int_t track=0; track<ntracks;track++) 
+    {
+      gAlice->ResetHits();
+      treeH->GetEvent(track);
+      
+      if (PMD) 
+       {
+         npmd = PMDhits->GetEntriesFast();
+         for (int ipmd = 0; ipmd < npmd; ipmd++) 
+           {
+             pmdHit = (AliPMDhit*) PMDhits->UncheckedAt(ipmd);
+             trackno = pmdHit->GetTrack();
+             
+             //  get kinematics of the particles
+             
+             particle = gAlice->Particle(trackno);
+             trackpid  = particle->GetPdgCode();
+
+             Int_t igatr = -999;
+             Int_t ichtr = -999;
+
+             Int_t imo;
+             Int_t igen = 0;
+             Int_t id_mo = -999;
+             TParticle*  mparticle = particle;
+             
+             while((imo = mparticle->GetFirstMother()) >= 0)
+               {
+                 igen++;
+                 mparticle =  gAlice->Particle(imo);
+                 id_mo = mparticle->GetPdgCode();
+                 
+                 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.)
+                   {
+                     igatr = imo;
+                   }
+                 ichtr = imo;
+               }
+
+             if (id_mo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
+               {
+                 mtrackno = igatr;
+                 mtrackpid = kGamma;
+               }
+             else
+               {
+                 mtrackno  = ichtr;
+                 mtrackpid = id_mo;
+               }
+             
+             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 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;
+               }
+             
+             smnumber = Convert2RealSMNumber(smnumber1);
+             ypad = (cellnumber - 1)/fNCell + 1;
+             xpad = cellnumber - (ypad-1)*fNCell;
+
+             //cout << "-zpos = " << -zPos << endl;
+             if (-zPos < fZPos)
+               {
+                 // CPV
+                 fDetNo = 1;
+               }
+             else if (-zPos > fZPos)
+               {
+                 // PMD
+                 fDetNo = 0;
+               }
+             
+             if(fDetNo == 1)
+               {
+                 fCPV[smnumber-1][xpad-1][ypad-1] += edep;
+               }
+             else if (fDetNo == 0)
+               {
+                 fPMD[smnumber-1][xpad-1][ypad-1] += edep;
+                 fPMDCounter[smnumber-1][xpad-1][ypad-1]++;
+               }
+           }
+       }
+    } // Track Loop ended
+
+  TrackAssignment2Cell();
+  ResetCell();
+
+  Float_t deltaE = 0.;
+  Int_t detno = 0;
+  Int_t count_digit = 0;
+  Int_t trno = 1;
+  Int_t cellno;
+
+  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;
+                 if (idet == 0)
+                   {
+                     deltaE = fPMD[ism][jrow][kcol];
+                     detno = 0;
+                   }
+                 else if (idet == 1)
+                   {
+                     deltaE = fCPV[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();
+  
+  //  cout << " -------- End of Hits2Digit ----------- " << endl;
+}
+
+
+void AliPMDDigitizer::SDigits2Digits(Int_t ievt)
+{
+  //  cout << " -------- Beginning of SDigits2Digit ----------- " << endl;
+  fRunLoader->GetEvent(ievt);
+
+  treeS = pmdloader->TreeS();
+  AliPMDsdigit  *pmdsdigit;
+  TBranch *branch = treeS->GetBranch("PMDSDigit");
+  branch->SetAddress(&fSDigits);
+
+  treeD = pmdloader->TreeD();
+  if (treeD == 0x0)
+    {
+      pmdloader->MakeTree("D");
+      treeD = pmdloader->TreeD();
+    }
+  Int_t bufsize = 16000;
+  treeD->Branch("PMDDigit", &fDigits, bufsize); 
+
+  Int_t   trno, det, smn;
+  Int_t   cellno;
+  Float_t edep, adc;
+
+  Int_t nmodules = (Int_t) treeS->GetEntries();
+
+  for (Int_t imodule = 0; imodule < nmodules; imodule++)
+    {
+      treeS->GetEntry(imodule); 
+      Int_t nentries = fSDigits->GetLast();
+      //cout << " nentries = " << nentries << endl;
+      for (Int_t ient = 0; ient < nentries+1; ient++)
+       {
+         pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
+         trno   = pmdsdigit->GetTrackNumber();
+         det    = pmdsdigit->GetDetector();
+         smn    = pmdsdigit->GetSMNumber();
+         cellno = pmdsdigit->GetCellNumber();
+         edep   = pmdsdigit->GetCellEdep();
+
+         MeV2ADC(edep,adc);
+
+         AddDigit(trno,det,smn,cellno,adc);      
+       }
+      treeD->Fill();
+      ResetDigit();
+    }
+  pmdloader->WriteDigits("OVERWRITE");
+  //  cout << " -------- End of SDigits2Digit ----------- " << endl;
+}
+
+void AliPMDDigitizer::TrackAssignment2Cell()
+{
+  // To be checked again
+  // 
+  // This blocks assign the cell id when there is 
+  // multiple tracks in a cell according to the
+  // energy deposition
+  //
+
+  Int_t i, j, k;
+
+  Float_t *frac_edp;
+  Float_t *tr_edp;
+  Int_t *status;
+  Int_t   ****PMDTrack;
+  Float_t ****PMDEdep;
+
+  PMDTrack = new Int_t ***[27];
+  PMDEdep  = new Float_t ***[27];
+  for (i=0; i<fTotSM; i++)
+    {
+      PMDTrack[i] = new Int_t **[72];
+      PMDEdep[i]  = new Float_t **[72];
+    }
+
+  for (i = 0; i < fTotSM; i++)
+    {
+      for (j = 0; j < fNCell; j++)
+       {
+         PMDTrack[i][j] = new Int_t *[72];
+         PMDEdep[i][j]  = new Float_t *[72];
+       }
+    }
+  
+  for (i = 0; i < fTotSM; i++)
+    {
+      for (j = 0; j < fNCell; j++)
+       {
+         for (k = 0; k < fNCell; k++)
+           {
+             Int_t nn = fPMDCounter[i][j][k];
+             if(nn > 0)
+               {
+                 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;
+           }
+       }
+    }
+
+
+  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;
+      
+      PMDTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
+      PMDEdep[ism][ixp][iyp][nn] = edep;
+      fPMDCounter[ism][ixp][iyp]++;
+    }
+  
+  Int_t iz, il;
+  Int_t im, ix, iy;
+  Int_t nn;
+  
+  for (im=0; im<27; im++)
+    {
+      for (ix=0; ix<72; ix++)
+       {
+         for (iy=0; iy<72; iy++)
+           {
+             nn = fPMDCounter[im][ix][iy];
+             if (nn > 1)
+               {
+                 // This block handles if a cell is fired
+                 // many times by many tracks
+                 
+                 status = new Int_t[nn];
+                 for (iz = 0; iz < nn; iz++)
+                   {
+                     status[iz] = PMDTrack[im][ix][iy][iz];
+                   }
+                 sort(status,status+nn);
+                 Int_t track_old = -99999;
+                 Int_t track, tr_count = 0;
+                 for (iz = 0; iz < nn; iz++)
+                   {
+                     track = status[iz];
+                     if (track_old != track)
+                       {
+                         tr_count++;
+                         vjunkTRN.push_back(track);
+                       }                             
+                     track_old = 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++)
+                   {
+                     tr_edp[il] = 0.;
+                     track = vjunkTRN[il];
+                     for (iz = 0; iz < nn; iz++)
+                       {
+                         if (track == PMDTrack[im][ix][iy][iz])
+                           {
+                             tr_edp[il] += PMDEdep[im][ix][iy][iz];
+                           }
+                       }
+                     tot_edp += tr_edp[il];
+                   }
+                 vjunkTRN.clear();
+                 Int_t il_old = 0;
+                 Float_t frac_old = 0.;
+                 
+                 for (il = 0; il < tr_count; il++)
+                   {
+                     frac_edp[il] = tr_edp[il]/tot_edp;
+                     if (frac_old < frac_edp[il])
+                       {
+                         frac_old = frac_edp[il];
+                         il_old = il;
+                       }
+                   }
+                 delete frac_edp;
+                 delete tr_edp;
+                 
+                 fPMDTrackNo[im][ix][iy] = vjunkTRN[il_old];
+               }
+             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];
+                 
+               }
+             else if (nn ==0)
+               {
+                 // This is if no cell is fired
+                 fPMDTrackNo[im][ix][iy] = -999;
+               }
+           } // end of iy
+       } // end of ix
+    } // end of im
+  
+  // Delete all the pointers
+  
+  for (i = 0; i < fTotSM; i++)
+    {
+      for (j = 0; j < fNCell; j++)
+       {
+         for (k = 0; k < fNCell; k++)
+           {
+             delete [] PMDTrack[i][j][k];
+             delete [] PMDEdep[i][j][k];
+           }
+       }
+    }
+  
+  for (i = 0; i < fTotSM; i++)
+    {
+      for (j = 0; j < fNCell; j++)
+       {
+         delete [] PMDTrack[i][j];
+         delete [] PMDEdep[i][j];
+       }
+    }
+  
+  for (i = 0; i < fTotSM; i++)
+    {
+      delete [] PMDTrack[i];
+      delete [] PMDEdep[i];
+    }
+  delete PMDTrack;
+  delete PMDEdep;
+  // 
+  // End of the cell id assignment
+  //
+}
+
+
+void AliPMDDigitizer::MeV2ADC(Float_t mev, Float_t & adc)
+{
+  // To be done
+
+  adc = mev*1.;
+}
+void AliPMDDigitizer::AddSDigit(Int_t trnumber, Int_t det, Int_t smnumber, 
+  Int_t cellnumber, Float_t adc)
+{
+  TClonesArray &lsdigits = *fSDigits;
+  AliPMDsdigit *newcell;
+  newcell = new AliPMDsdigit(trnumber,det,smnumber,cellnumber,adc);
+  new(lsdigits[fNsdigit++]) AliPMDsdigit(newcell);
+  delete newcell;
+}
+
+void AliPMDDigitizer::AddDigit(Int_t trnumber, Int_t det, Int_t smnumber, 
+  Int_t cellnumber, Float_t adc)
+{
+  TClonesArray &ldigits = *fDigits;
+  AliPMDdigit *newcell;
+  newcell = new AliPMDdigit(trnumber,det,smnumber,cellnumber,adc);
+  new(ldigits[fNdigit++]) AliPMDdigit(newcell);
+  delete newcell;
+}
+
+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++)
+    {
+      for (Int_t j = 0; j < fNCell; j++)
+       {
+         for (Int_t k = 0; k < fNCell; k++)
+           {
+             fPMDCounter[i][j][k] = 0; 
+           }
+       }
+    }
+}
+void AliPMDDigitizer::ResetSDigit()
+{
+  fNsdigit = 0;
+  if (fSDigits) fSDigits->Clear();
+}
+void AliPMDDigitizer::ResetDigit()
+{
+  fNdigit = 0;
+  if (fDigits) fDigits->Clear();
+}
+
+void AliPMDDigitizer::ResetCellADC()
+{
+  for (Int_t i = 0; i < fTotSM; i++)
+    {
+      for (Int_t j = 0; j < fNCell; j++)
+       {
+         for (Int_t k = 0; k < fNCell; k++)
+           {
+             fCPV[i][j][k] = 0.; 
+             fPMD[i][j][k] = 0.; 
+           }
+       }
+    }
+}
+
+void AliPMDDigitizer::UnLoad(Option_t *option)
+{
+  const char *cS = strstr(option,"S");
+  const char *cD = strstr(option,"D");
+
+  fRunLoader->UnloadgAlice();
+  fRunLoader->UnloadHeader();
+  fRunLoader->UnloadKinematics();
+
+  if (cS)
+    {
+      pmdloader->UnloadHits();
+    }
+  if (cD)
+    {
+      pmdloader->UnloadHits();
+      pmdloader->UnloadSDigits();
+    }
+}
diff --git a/PMD/AliPMDDigitizer.h b/PMD/AliPMDDigitizer.h
new file mode 100644 (file)
index 0000000..d62fe78
--- /dev/null
@@ -0,0 +1,111 @@
+#ifndef PMDDigitizer_H
+#define PMDDigitizer_H
+//-----------------------------------------------------//
+//                                                     //
+//  Header File : PMDDigitization.h, Version 00        //
+//                                                     //
+//  Date   : September 20 2002                         //
+//                                                     //
+//-----------------------------------------------------//
+
+#include <Riostream.h>
+#include <stdlib.h>
+#include <math.h>
+#include <TMath.h>
+#include <vector>
+#include <algorithm>
+
+class TClonesArray;
+class TFile;
+
+class TObjArray;
+class TParticle;
+class TTree;
+class TNtuple;
+
+class AliLoader;
+class AliRunLoader;
+class AliRun;
+class AliDetector;
+class AliPMDhit;
+class AliHit;
+class AliHeader;
+
+class AliPMDcell;
+class AliPMDsdigit;
+class AliPMDdigit;
+class AliPMDClustering;
+class AliPMDContainer;
+class AliPMDrecpoint;
+
+using namespace std;
+
+class AliPMDDigitizer
+{
+ protected:
+  AliRunLoader *fRunLoader;
+  AliRun       *gAlice;
+  AliPMDhit    *pmdHit;   /* Pointer to specific detector hits. */
+  AliDetector  *PMD;      /* Get pointers to Alice detectors 
+                            and Hits containers */
+  AliLoader    *pmdloader;
+
+  TClonesArray *PMDhits;
+  TObjArray    *Particles;
+  TParticle    *particle;
+
+  TTree        *treeH;
+  TTree        *treeS;
+  TTree        *treeD;
+  TTree        *treeR;
+
+  TClonesArray *fSDigits;
+  TClonesArray *fDigits;
+  TClonesArray *fRecpoints;
+
+  TObjArray    *fCell;
+  AliPMDcell   *pmdcell;
+  vector<int> vjunkTRN;
+
+  Int_t fNsdigit;
+  Int_t fNdigit;
+  Int_t fDetNo;
+  Float_t fZPos;
+
+  static const Int_t fTotSM = 27;
+  static const Int_t fNCell = 72;
+  Float_t fCPV[fTotSM][fNCell][fNCell];
+  Float_t fPMD[fTotSM][fNCell][fNCell];
+  Int_t   fPMDCounter[fTotSM][fNCell][fNCell];
+  Int_t   fPMDTrackNo[fTotSM][fNCell][fNCell];
+  Int_t   fCPVTrackNo[fTotSM][fNCell][fNCell];
+
+ public:
+
+  AliPMDDigitizer();
+  virtual ~AliPMDDigitizer();
+
+  void OpengAliceFile(char * /* galice.root */, Option_t * /* option */);
+
+  void Hits2SDigits(Int_t /* ievt */);
+  void Hits2Digits(Int_t /* ievt */);
+  void SDigits2Digits(Int_t /* ievt */);
+  void TrackAssignment2Cell();
+  void MeV2ADC(Float_t /* mev */, Float_t & /* adc */);
+  void AddSDigit(Int_t /* trnumber */, Int_t /* det */, Int_t /* smnumber */, 
+                Int_t /* cellnumber */, Float_t /* adc */);
+  void AddDigit(Int_t /* trnumber */, Int_t /* det */, Int_t /* smnumber */, 
+               Int_t /* cellnumber */, Float_t /* adc */);
+  Int_t Convert2RealSMNumber(Int_t /* smnumber1 */ );
+  void SetZPosition(Float_t /* zpos */);
+  Float_t GetZPosition() const;
+  void ResetCell();
+  void ResetSDigit();
+  void ResetDigit();
+  void ResetCellADC();
+  void UnLoad(Option_t * /* option */);
+
+  ClassDef(AliPMDDigitizer,1)
+};
+#endif
+