]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSv1.cxx
Effective C++ corrections (T.Pocheptsov)
[u/mrichter/AliRoot.git] / PHOS / AliPHOSv1.cxx
index a86a32406fc74bb383490d19ac769596a0323d4c..6d221a3920bcbf51a3c14f427731542d1003a632 100644 (file)
 
 /* $Id$ */
 
+/* History of cvs commits:
+ *
+ * $Log$
+ * Revision 1.104  2005/05/28 14:19:05  schutz
+ * Compilation warnings fixed by T.P.
+ *
+ */
+
 //_________________________________________________________________________
 // Implementation version v1 of PHOS Manager class 
 //---
-// Layout EMC + PPSD has name GPS2:
-// Produces cumulated hits
 //---
 // Layout EMC + CPV  has name IHEP:
 // Produces hits for CPV, cumulated hits
 //---
-// Layout EMC + CPV + PPSD has name GPS:
-// Produces hits for CPV, cumulated hits
 //---
 //*-- Author: Yves Schutz (SUBATECH)
 
 
 // --- ROOT system ---
-
-#include "TBRIK.h"
-#include "TNode.h"
-#include "TRandom.h"
-#include "TTree.h"
-
+#include <TParticle.h>
+#include <TVirtualMC.h>
 
 // --- Standard library ---
 
-#include <stdio.h>
-#include <string.h>
-#include <stdlib.h>
-#include <strstream.h>
 
 // --- AliRoot header files ---
-
-#include "AliPHOSv1.h"
-#include "AliPHOSHit.h"
 #include "AliPHOSCPVDigit.h"
+#include "AliPHOSGeometry.h"
+#include "AliPHOSHit.h"
+#include "AliPHOSv1.h"
 #include "AliRun.h"
-#include "AliConst.h"
 #include "AliMC.h"
-#include "AliPHOSGeometry.h"
 
 ClassImp(AliPHOSv1)
 
 //____________________________________________________________________________
 AliPHOSv1::AliPHOSv1():
-AliPHOSv0()
+  fLightYieldMean(0.),
+  fIntrinsicPINEfficiency(0.),
+  fLightYieldAttenuation(0.),
+  fRecalibrationFactor(0.),
+  fElectronsPerGeV(0.),
+  fAPDGain(0.),
+  fLightFactor(0.),
+  fAPDFactor(0.)
 {
-  // ctor
+  //Def ctor.
 }
 
 //____________________________________________________________________________
 AliPHOSv1::AliPHOSv1(const char *name, const char *title):
-AliPHOSv0(name,title) 
+  AliPHOSv0(name,title),
+  fLightYieldMean(0.),
+  fIntrinsicPINEfficiency(0.),
+  fLightYieldAttenuation(0.),
+  fRecalibrationFactor(0.),
+  fElectronsPerGeV(0.),
+  fAPDGain(0.),
+  fLightFactor(0.),
+  fAPDFactor(0.)
 {
-  // ctor : title is used to identify the layout
-  //        GPS2 = 5 modules (EMC + PPSD)
-  //        IHEP = 5 modules (EMC + CPV )
-  //        MIXT = 4 modules (EMC + CPV ) and 1 module (EMC + PPSD)
   //
   // We store hits :
   //   - fHits (the "normal" one), which retains the hits associated with
@@ -85,239 +89,343 @@ AliPHOSv0(name,title)
   // We do not want to save in TreeH the raw hits
   // But save the cumulated hits instead (need to create the branch myself)
   // It is put in the Digit Tree because the TreeH is filled after each primary
-  // and the TreeD at the end of the event (branch is set in FinishEvent() ).
+  // and the TreeD at the end of the event (branch is set in FinishEvent() ). 
   
   fHits= new TClonesArray("AliPHOSHit",1000) ;
+  gAlice->GetMCApp()->AddHitList(fHits) ; 
 
   fNhits = 0 ;
 
-  fIshunt     =  1 ; // All hits are associated with primary particles
-  
+  fIshunt     =  2 ; // All hits are associated with primary particles
+
+  //Photoelectron statistics:
+  // The light yield is a poissonian distribution of the number of
+  // photons created in the PbWo4 crystal, calculated using following formula
+  // NumberOfPhotons = EnergyLost * LightYieldMean* APDEfficiency *
+  //              exp (-LightYieldAttenuation * DistanceToPINdiodeFromTheHit);
+  // LightYieldMean is parameter calculated to be over 47000 photons per GeV
+  // APDEfficiency is 0.02655
+  // k_0 is 0.0045 from Valery Antonenko
+  // The number of electrons created in the APD is
+  // NumberOfElectrons = APDGain * LightYield
+  // The APD Gain is 300
+  fLightYieldMean = 47000;
+  fIntrinsicPINEfficiency = 0.02655 ; //APD= 0.1875/0.1271 * 0.018 (PIN)
+  fLightYieldAttenuation  = 0.0045 ; 
+  fRecalibrationFactor    = 13.418/ fLightYieldMean ;
+  fElectronsPerGeV        = 2.77e+8 ;
+  fAPDGain                = 300. ;
+  fLightFactor            = fLightYieldMean * fIntrinsicPINEfficiency ; 
+  fAPDFactor              = (fRecalibrationFactor/100.) * fAPDGain ;   
+}
+
+AliPHOSv1::AliPHOSv1(AliPHOSv1 & phos) :
+  AliPHOSv0(phos),
+  fLightYieldMean(0.),
+  fIntrinsicPINEfficiency(0.),
+  fLightYieldAttenuation(0.),
+  fRecalibrationFactor(0.),
+  fElectronsPerGeV(0.),
+  fAPDGain(0.),
+  fLightFactor(0.),
+  fAPDFactor(0.)
+{
+  //Copy ctor. Can be wrong.
+  phos.Copy(*this) ; 
 }
 
+
 //____________________________________________________________________________
 AliPHOSv1::~AliPHOSv1()
 {
   // dtor
-
-  if ( fHits) {
+ if ( fHits) {
     fHits->Delete() ; 
     delete fHits ;
     fHits = 0 ; 
-  }
+ }
+}
+
+//____________________________________________________________________________
+void AliPHOSv1::Copy(TObject & base)const
+{
+  TObject::Copy(base) ; 
+  AliPHOSv0::Copy(base) ;
+  AliPHOSv1 &phos = static_cast<AliPHOSv1 &>(base); 
+  phos.fLightYieldMean         = fLightYieldMean ; 
+  phos.fIntrinsicPINEfficiency = fIntrinsicPINEfficiency ; 
+  phos.fLightYieldAttenuation  = fLightYieldAttenuation ; 
+  phos.fRecalibrationFactor    = fRecalibrationFactor ; 
+  phos.fElectronsPerGeV        = fElectronsPerGeV ; 
+  phos.fAPDGain                = fAPDGain ; 
+  phos.fLightFactor            = fLightFactor ; 
+  phos.fAPDFactor              = fAPDFactor ; 
 }
 
 //____________________________________________________________________________
-void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits)
+void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t Id, Float_t * hits)
 {
   // Add a hit to the hit list.
-  // A PHOS hit is the sum of all hits in a single crystal
-  //   or in a single PPSD gas cell
+  // A PHOS hit is the sum of all hits in a single crystal from one primary and within some time gate
 
   Int_t hitCounter ;
   AliPHOSHit *newHit ;
   AliPHOSHit *curHit ;
   Bool_t deja = kFALSE ;
+  AliPHOSGeometry * geom = GetGeometry() ; 
 
-  newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits) ;
+  newHit = new AliPHOSHit(shunt, primary, Id, hits) ;
 
   for ( hitCounter = fNhits-1 ; hitCounter >= 0 && !deja ; hitCounter-- ) {
-    curHit = (AliPHOSHit*) (*fHits)[hitCounter] ;
-    if(curHit->GetPrimary() != primary) break ; // We add hits with the same primary, while GEANT treats primaries succesively 
+    curHit = dynamic_cast<AliPHOSHit*>((*fHits)[hitCounter]) ;
+    if(curHit->GetPrimary() != primary) break ; 
+           // We add hits with the same primary, while GEANT treats primaries succesively 
     if( *curHit == *newHit ) {
-      *curHit = *curHit + *newHit ;
+      *curHit + *newHit ;
       deja = kTRUE ;
     }
   }
          
   if ( !deja ) {
     new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;
+    // get the block Id number
+    Int_t relid[4] ;
+    geom->AbsToRelNumbering(Id, relid) ;
+
     fNhits++ ;
   }
-
+  
   delete newHit;
 }
 
+//____________________________________________________________________________
+void AliPHOSv1::FinishPrimary() 
+{
+  // called at the end of each track (primary) by AliRun
+  // hits are reset for each new track
+  // accumulate the total hit-multiplicity
+
+}
+
+//____________________________________________________________________________
+void AliPHOSv1::FinishEvent() 
+{
+  // called at the end of each event by AliRun
+  // accumulate the hit-multiplicity and total energy per block 
+  // if the values have been updated check it
+  
+  AliDetector::FinishEvent(); 
+}
 //____________________________________________________________________________
 void AliPHOSv1::StepManager(void)
 {
-  // Accumulates hits as long as the track stays in a single crystal or PPSD gas Cell
+   // Accumulates hits as long as the track stays in a single crystal or CPV gas Cell
 
   Int_t          relid[4] ;           // (box, layer, row, column) indices
   Int_t          absid    ;           // absolute cell ID number
-  Float_t        xyze[4]={0,0,0,0}  ; // position wrt MRS and energy deposited
+  Float_t        xyzte[5]={-1000.,-1000.,-1000.,0.,0.}  ; // position wrt MRS, time and energy deposited
   TLorentzVector pos      ;           // Lorentz vector of the track current position
   Int_t          copy     ;
 
-  Int_t tracknumber =  gAlice->CurrentTrack() ; 
-  Int_t primary     =  gAlice->GetPrimary( gAlice->CurrentTrack() ); 
-  TString name      =  fGeom->GetName() ; 
-
-
-  if ( name == "GPS2" || name == "MIXT" ) {            // ======> CPV is a GPS' PPSD
-
-    if( gMC->CurrentVolID(copy) == gMC->VolId("PPCE") ) // We are inside a gas cell 
-    {
-      gMC->TrackPosition(pos) ;
-      xyze[0] = pos[0] ;
-      xyze[1] = pos[1] ;
-      xyze[2] = pos[2] ;
-      xyze[3] = gMC->Edep() ; 
-
-      if ( xyze[3] != 0 ) { // there is deposited energy 
-               gMC->CurrentVolOffID(5, relid[0]) ;  // get the PHOS Module number
-       if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(5),"PHO1") == 0 ){
-         relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules();
-       }
-               gMC->CurrentVolOffID(3, relid[1]) ;  // get the Micromegas Module number 
-      // 1-> fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() upper
-      //   > fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() lower
-               gMC->CurrentVolOffID(1, relid[2]) ;  // get the row number of the cell
-        gMC->CurrentVolID(relid[3]) ;        // get the column number 
-
-       // get the absolute Id number
-
-               fGeom->RelToAbsNumbering(relid, absid) ; 
-
-       // add current hit to the hit list      
-         AddHit(fIshunt, primary, tracknumber, absid, xyze);
-
+  TString name      =  GetGeometry()->GetName() ; 
 
-      } // there is deposited energy 
-    } // We are inside the gas of the CPV  
-  } // GPS2 configuration
-
-  if ( name == "IHEP" || name == "MIXT" ) {       // ======> CPV is a IHEP's one
-
-    // Yuri Kharlov, 28 September 2000
-
-    if( gMC->CurrentVolID(copy) == gMC->VolId("PCPQ") &&
-       (gMC->IsTrackEntering() ) &&
-       gMC->TrackCharge() != 0) {      
+  Int_t moduleNumber ;
+  
+  static Int_t idPCPQ = gMC->VolId("PCPQ");
+  if( gMC->CurrentVolID(copy) == idPCPQ &&
+      (gMC->IsTrackEntering() ) &&
+      gMC->TrackCharge() != 0) {      
+    
+    gMC -> TrackPosition(pos);
+    
+    Float_t xyzm[3], xyzd[3] ;
+    Int_t i;
+    for (i=0; i<3; i++) xyzm[i] = pos[i];
+    gMC -> Gmtod (xyzm, xyzd, 1);    // transform coordinate from master to daughter system
+    
+    Float_t        xyd[3]={0,0,0}   ;   //local position of the entering
+    xyd[0]  = xyzd[0];
+    xyd[1]  =-xyzd[2];
+    xyd[2]  =-xyzd[1];
+    
+    // Current momentum of the hit's track in the local ref. system
+    TLorentzVector pmom     ;        //momentum of the particle initiated hit
+    gMC -> TrackMomentum(pmom);
+    Float_t pm[3], pd[3];
+    for (i=0; i<3; i++)  
+      pm[i]   = pmom[i];
+    
+    gMC -> Gmtod (pm, pd, 2);        // transform 3-momentum from master to daughter system
+    pmom[0] = pd[0];
+    pmom[1] =-pd[1];
+    pmom[2] =-pd[2];
+
+    // Digitize the current CPV hit:
+    
+    // 1. find pad response and    
+    gMC->CurrentVolOffID(3,moduleNumber);
+    moduleNumber--;
+    
+    TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0);   // array of digits for current hit
+    CPVDigitize(pmom,xyd,cpvDigits);
       
-      gMC -> TrackPosition(pos);
-      Float_t xyzm[3], xyzd[3] ;
-      Int_t i;
-      for (i=0; i<3; i++) xyzm[i] = pos[i];
-      gMC -> Gmtod (xyzm, xyzd, 1);    // transform coordinate from master to daughter system
-
-      Float_t        xyd[3]={0,0,0}   ;   //local posiiton of the entering
-      xyd[0]  = xyzd[0];
-      xyd[1]  =-xyzd[1];
-      xyd[2]  =-xyzd[2];
-
+    Float_t xmean = 0;
+    Float_t zmean = 0;
+    Float_t qsum  = 0;
+    Int_t   idigit,ndigits;
+    
+    // 2. go through the current digit list and sum digits in pads
+    
+    ndigits = cpvDigits->GetEntriesFast();
+    for (idigit=0; idigit<ndigits-1; idigit++) {
+      AliPHOSCPVDigit  *cpvDigit1 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
+      Float_t x1 = cpvDigit1->GetXpad() ;
+      Float_t z1 = cpvDigit1->GetYpad() ;
+      for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
+       AliPHOSCPVDigit  *cpvDigit2 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(jdigit));
+       Float_t x2 = cpvDigit2->GetXpad() ;
+       Float_t z2 = cpvDigit2->GetYpad() ;
+       if (x1==x2 && z1==z2) {
+         Float_t qsum = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
+         cpvDigit2->SetQpad(qsum) ;
+         cpvDigits->RemoveAt(idigit) ;
+       }
+      }
+    }
+    cpvDigits->Compress() ;
+    
+    // 3. add digits to temporary hit list fTmpHits
+    
+    ndigits = cpvDigits->GetEntriesFast();
+    for (idigit=0; idigit<ndigits; idigit++) {
+      AliPHOSCPVDigit  *cpvDigit = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
+      relid[0] = moduleNumber + 1 ;                             // CPV (or PHOS) module number
+      relid[1] =-1 ;                                            // means CPV
+      relid[2] = cpvDigit->GetXpad() ;                          // column number of a pad
+      relid[3] = cpvDigit->GetYpad() ;                          // row    number of a pad
       
-      // Current momentum of the hit's track in the local ref. system
-        TLorentzVector pmom     ;        //momentum of the particle initiated hit
-      gMC -> TrackMomentum(pmom);
-      Float_t pm[3], pd[3];
-      for (i=0; i<3; i++) pm[i]   = pmom[i];
-      gMC -> Gmtod (pm, pd, 2);        // transform 3-momentum from master to daughter system
-      pmom[0] = pd[0];
-      pmom[1] =-pd[1];
-      pmom[2] =-pd[2];
-
-      // Digitize the current CPV hit:
-
-      // 1. find pad response and
+      // get the absolute Id number
+      GetGeometry()->RelToAbsNumbering(relid, absid) ; 
       
-      Int_t moduleNumber;
-      gMC->CurrentVolOffID(3,moduleNumber);
-      moduleNumber--;
+      // add current digit to the temporary hit list
 
+      xyzte[3] = gMC->TrackTime() ;
+      xyzte[4] = cpvDigit->GetQpad() ;                          // amplitude in a pad
 
-      TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0);   // array of digits for current hit
-      CPVDigitize(pmom,xyd,moduleNumber,cpvDigits);
+      Int_t primary  =  gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() ); 
+      AddHit(fIshunt, primary, absid, xyzte);  
       
-      Float_t xmean = 0;
-      Float_t zmean = 0;
-      Float_t qsum  = 0;
-      Int_t   idigit,ndigits;
-
-      // 2. go through the current digit list and sum digits in pads
-
-      ndigits = cpvDigits->GetEntriesFast();
-      for (idigit=0; idigit<ndigits-1; idigit++) {
-       AliPHOSCPVDigit  *cpvDigit1 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
-       Float_t x1 = cpvDigit1->GetXpad() ;
-       Float_t z1 = cpvDigit1->GetYpad() ;
-       for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
-         AliPHOSCPVDigit  *cpvDigit2 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(jdigit);
-         Float_t x2 = cpvDigit2->GetXpad() ;
-         Float_t z2 = cpvDigit2->GetYpad() ;
-         if (x1==x2 && z1==z2) {
-           Float_t qsum = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
-           cpvDigit2->SetQpad(qsum) ;
-           cpvDigits->RemoveAt(idigit) ;
-         }
-       }
+      if (cpvDigit->GetQpad() > 0.02) {
+       xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
+       zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
+       qsum  += cpvDigit->GetQpad();
       }
-      cpvDigits->Compress() ;
+    }
+    if (cpvDigits) {
+      cpvDigits->Delete();
+      delete cpvDigits;
+      cpvDigits=0;
+    }
+  }
 
-      // 3. add digits to temporary hit list fTmpHits
+  static Int_t idPXTL = gMC->VolId("PXTL");  
+  if(gMC->CurrentVolID(copy) == idPXTL ) { //  We are inside a PBWO crystal
 
-      ndigits = cpvDigits->GetEntriesFast();
-      for (idigit=0; idigit<ndigits; idigit++) {
-       AliPHOSCPVDigit  *cpvDigit = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
-       relid[0] = moduleNumber + 1 ;                             // CPV (or PHOS) module number
-       relid[1] =-1 ;                                            // means CPV
-       relid[2] = cpvDigit->GetXpad() ;                          // column number of a pad
-       relid[3] = cpvDigit->GetYpad() ;                          // row    number of a pad
-       
-       // get the absolute Id number
-       fGeom->RelToAbsNumbering(relid, absid) ; 
-
-       // add current digit to the temporary hit list
-       xyze[0] = 0. ;
-       xyze[1] = 0. ;
-       xyze[2] = 0. ;
-       xyze[3] = cpvDigit->GetQpad() ;                           // amplitude in a pad
-       primary = -1;                                             // No need in primary for CPV
-       AddHit(fIshunt, primary, tracknumber, absid, xyze);
-
-       if (cpvDigit->GetQpad() > 0.02) {
-         xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
-         zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
-         qsum  += cpvDigit->GetQpad();
+    gMC->TrackPosition(pos) ;
+    xyzte[0] = pos[0] ;
+    xyzte[1] = pos[1] ;
+    xyzte[2] = pos[2] ;
+
+    Float_t global[3], local[3] ;
+    global[0] = pos[0] ;
+    global[1] = pos[1] ;
+    global[2] = pos[2] ;
+    Float_t lostenergy = gMC->Edep(); 
+    
+    //Put in the TreeK particle entering PHOS and all its parents
+    if ( gMC->IsTrackEntering() ){
+      Float_t xyzd[3] ;
+      gMC -> Gmtod (xyzte, xyzd, 1);    // transform coordinate from master to daughter system    
+      if (xyzd[1] < -GetGeometry()->GetCrystalSize(1)/2.+0.1){   //Entered close to forward surface  
+       Int_t parent = gAlice->GetMCApp()->GetCurrentTrackNumber() ; 
+       TParticle * part = gAlice->GetMCApp()->Particle(parent) ; 
+       Float_t vert[3],vertd[3] ;
+       vert[0]=part->Vx() ;
+       vert[1]=part->Vy() ;
+       vert[2]=part->Vz() ;
+       gMC -> Gmtod (vert, vertd, 1);    // transform coordinate from master to daughter system
+       if(vertd[1]<-GetGeometry()->GetCrystalSize(1)/2.-0.1){ //Particle is created in foront of PHOS 
+                                                              //0.1 to get rid of numerical errors 
+         part->SetBit(kKeepBit);
+         while ( parent != -1 ) {
+           part = gAlice->GetMCApp()->Particle(parent) ; 
+           part->SetBit(kKeepBit);
+           parent = part->GetFirstMother() ; 
+         }
        }
       }
-      delete cpvDigits;
     }
-  } // end of IHEP configuration
-  
-
-  if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) { //  We are inside a PBWO crystal
-    gMC->TrackPosition(pos) ;
-    xyze[0] = pos[0] ;
-    xyze[1] = pos[1] ;
-    xyze[2] = pos[2] ;
-    xyze[3] = gMC->Edep() ;
-
-  
-    if ( xyze[3] != 0 ) {  // Track is inside the crystal and deposits some energy 
-
-      gMC->CurrentVolOffID(10, relid[0]) ; // get the PHOS module number ;
-
-      if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(10),"PHO1") == 0 )
-       relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules();      
-
-      relid[1] = 0   ;                    // means PBW04
-      gMC->CurrentVolOffID(4, relid[2]) ; // get the row number inside the module
-      gMC->CurrentVolOffID(3, relid[3]) ; // get the cell number inside the module
+    if ( lostenergy != 0 ) {  // Track is inside the crystal and deposits some energy 
+      xyzte[3] = gMC->TrackTime() ;     
       
-      // get the absolute Id number
-      fGeom->RelToAbsNumbering(relid, absid) ; 
+      gMC->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ;
+      
+      Int_t strip ;
+      gMC->CurrentVolOffID(3, strip);
+      Int_t cell ;
+      gMC->CurrentVolOffID(2, cell);
+      
+      Int_t row = 1 + GetGeometry()->GetNZ() - strip % GetGeometry()->GetNZ() ;
+      Int_t col = (Int_t) TMath::Ceil((Double_t) strip/GetGeometry()->GetNZ()) -1 ;
+      
+      absid = (moduleNumber-1)*GetGeometry()->GetNCristalsInModule() + 
+       row + (col*GetGeometry()->GetEMCAGeometry()->GetNCellsInStrip() + cell-1)*GetGeometry()->GetNZ() ;
+      
+      gMC->Gmtod(global, local, 1) ;
+      
+      //Calculates the light yield, the number of photons produced in the
+      //crystal 
+      Float_t lightYield = gRandom->Poisson(fLightFactor * lostenergy *
+                                           exp(-fLightYieldAttenuation *
+                                               (local[1]+GetGeometry()->GetCrystalSize(1)/2.0 ))
+                                           ) ;
+
+      //Calculates de energy deposited in the crystal  
+      xyzte[4] = fAPDFactor * lightYield  ;
+      
+      Int_t primary ;
+      if(fIshunt == 2){
+       primary = gAlice->GetMCApp()->GetCurrentTrackNumber() ;
+       TParticle * part = gAlice->GetMCApp()->Particle(primary) ;
+       while ( !part->TestBit(kKeepBit) ) {
+         primary = part->GetFirstMother() ;
+         if(primary == -1){        
+           primary  =  gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() ); 
+           break ; //there is a possibility that particle passed e.g. thermal isulator and hits a side 
+         //surface of the crystal. In this case it may have no primary at all. 
+         //We can not easily separate this case from the case when this is part of the shower, 
+         //developed in the neighboring crystal.
+         }
+         part = gAlice->GetMCApp()->Particle(primary) ;
+       }
+      }
+      else
+       primary  =  gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() ); 
 
+      
+      
       // add current hit to the hit list
-       AddHit(fIshunt, primary,tracknumber, absid, xyze);
-
-
+      // Info("StepManager","%d %d", primary, tracknumber) ; 
+      AddHit(fIshunt, primary, absid, xyzte);
+        
     } // there is deposited energy
   } // we are inside a PHOS Xtal
+  
 }
 
 //____________________________________________________________________________
-void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, Int_t moduleNumber, TClonesArray *cpvDigits)
+void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, TClonesArray *cpvDigits)
 {
   // ------------------------------------------------------------------------
   // Digitize one CPV hit:
@@ -329,7 +437,7 @@ void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, Int_t moduleNumbe
   // 2 October 2000
   // ------------------------------------------------------------------------
 
-  const Float_t kCelWr  = fGeom->GetPadSizePhi()/2;  // Distance between wires (2 wires above 1 pad)
+  const Float_t kCelWr  = GetGeometry()->GetPadSizePhi()/2;  // Distance between wires (2 wires above 1 pad)
   const Float_t kDetR   = 0.1;     // Relative energy fluctuation in track for 100 e-
   const Float_t kdEdx   = 4.0;     // Average energy loss in CPV;
   const Int_t   kNgamz  = 5;       // Ionization size in Z
@@ -350,15 +458,15 @@ void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, Int_t moduleNumbe
   Float_t pNorm = p.Py();
   Float_t eloss = kdEdx;
 
-//    cout << "CPVDigitize: YVK : "<<hitX<<" "<<hitZ<<" | "<<pX<<" "<<pZ<<" "<<pNorm<<endl;
+//Info("CPVDigitize", "YVK : %f %f | %f %f %d", hitX, hitZ, pX, pZ, pNorm) ;
 
-  Float_t dZY   = pZ/pNorm * fGeom->GetCPVGasThickness();
-  Float_t dXY   = pX/pNorm * fGeom->GetCPVGasThickness();
+  Float_t dZY   = pZ/pNorm * GetGeometry()->GetCPVGasThickness();
+  Float_t dXY   = pX/pNorm * GetGeometry()->GetCPVGasThickness();
   gRandom->Rannor(rnor1,rnor2);
   eloss *= (1 + kDetR*rnor1) *
-           TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(fGeom->GetCPVGasThickness(),2)));
-  Float_t zhit1 = hitZ + fGeom->GetCPVActiveSize(1)/2 - dZY/2;
-  Float_t xhit1 = hitX + fGeom->GetCPVActiveSize(0)/2 - dXY/2;
+           TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(GetGeometry()->GetCPVGasThickness(),2)));
+  Float_t zhit1 = hitZ + GetGeometry()->GetCPVActiveSize(1)/2 - dZY/2;
+  Float_t xhit1 = hitX + GetGeometry()->GetCPVActiveSize(0)/2 - dXY/2;
   Float_t zhit2 = zhit1 + dZY;
   Float_t xhit2 = xhit1 + dXY;
 
@@ -418,20 +526,20 @@ void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, Int_t moduleNumbe
 
   // Finite size of ionization region
 
-  Int_t nCellZ  = fGeom->GetNumberOfCPVPadsZ();
-  Int_t nCellX  = fGeom->GetNumberOfCPVPadsPhi();
+  Int_t nCellZ  = GetGeometry()->GetNumberOfCPVPadsZ();
+  Int_t nCellX  = GetGeometry()->GetNumberOfCPVPadsPhi();
   Int_t nz3     = (kNgamz+1)/2;
   Int_t nx3     = (kNgamx+1)/2;
   cpvDigits->Expand(nIter*kNgamx*kNgamz);
-  TClonesArray &ldigits = *(TClonesArray *)cpvDigits;
+  TClonesArray &ldigits = *(static_cast<TClonesArray *>(cpvDigits));
 
   for (Int_t iter=0; iter<nIter; iter++) {
 
     Float_t zhit = zxe[0][iter];
     Float_t xhit = zxe[1][iter];
     Float_t qhit = zxe[2][iter];
-    Float_t zcell = zhit / fGeom->GetPadSizeZ();
-    Float_t xcell = xhit / fGeom->GetPadSizePhi();
+    Float_t zcell = zhit / GetGeometry()->GetPadSizeZ();
+    Float_t xcell = xhit / GetGeometry()->GetPadSizePhi();
     if ( zcell<=0      || xcell<=0 ||
         zcell>=nCellZ || xcell>=nCellX) return;
     Int_t izcell = (Int_t) zcell;
@@ -468,10 +576,10 @@ Float_t AliPHOSv1::CPVPadResponseFunction(Float_t qhit, Float_t zhit, Float_t xh
   // 3 October 2000
   // ------------------------------------------------------------------------
 
-  Double_t dz = fGeom->GetPadSizeZ()   / 2;
-  Double_t dx = fGeom->GetPadSizePhi() / 2;
-  Double_t z  = zhit * fGeom->GetPadSizeZ();
-  Double_t x  = xhit * fGeom->GetPadSizePhi();
+  Double_t dz = GetGeometry()->GetPadSizeZ()   / 2;
+  Double_t dx = GetGeometry()->GetPadSizePhi() / 2;
+  Double_t z  = zhit * GetGeometry()->GetPadSizeZ();
+  Double_t x  = xhit * GetGeometry()->GetPadSizePhi();
   Double_t amplitude = qhit *
     (CPVCumulPadResponse(z+dz,x+dx) - CPVCumulPadResponse(z+dz,x-dx) -
      CPVCumulPadResponse(z-dz,x+dx) + CPVCumulPadResponse(z-dz,x-dx));