[u/mrichter/AliRoot.git] / PMD / AliPMDDigitizer.cxx

//-----------------------------------------------------//
//                                                     //
//  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"

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