]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSv1.cxx
Bug correction
[u/mrichter/AliRoot.git] / PHOS / AliPHOSv1.cxx
index 582d834f6f62cf0a7216f89031244904b08a336c..d8cc5f0d9dce57d17248f06ca5f0f7242fca3997 100644 (file)
 //_________________________________________________________________________
 // 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 "AliRun.h"
-#include "AliConst.h"
-#include "AliMC.h"
 #include "AliPHOSGeometry.h"
-#include "AliPHOSQAIntCheckable.h"
+#include "AliPHOSHit.h"
 #include "AliPHOSQAFloatCheckable.h"
+#include "AliPHOSQAIntCheckable.h"
 #include "AliPHOSQAMeanChecker.h"
+#include "AliPHOSv1.h"
+#include "AliRun.h"
+#include "AliMC.h"
 
 ClassImp(AliPHOSv1)
 
@@ -64,18 +50,27 @@ ClassImp(AliPHOSv1)
 AliPHOSv1::AliPHOSv1():
 AliPHOSv0()
 {
-  // ctor
+  // default ctor: initialze data memebers
+  fQAHitsMul  = 0 ;
+  fQAHitsMulB = 0 ; 
+  fQATotEner  = 0 ; 
+  fQATotEnerB = 0 ; 
+
+  fLightYieldMean         = 0. ;         
+  fIntrinsicPINEfficiency = 0. ; 
+  fLightYieldAttenuation  = 0. ;  
+  fRecalibrationFactor    = 0. ;    
+  fElectronsPerGeV        = 0. ;
+  fAPDGain                = 0. ;  
+  fLightFactor            = 0. ; 
+  fAPDFactor              = 0. ; 
+
 }
 
 //____________________________________________________________________________
 AliPHOSv1::AliPHOSv1(const char *name, const char *title):
  AliPHOSv0(name,title) 
 {
-  // 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
@@ -91,18 +86,42 @@ AliPHOSv1::AliPHOSv1(const char *name, const char *title):
   // 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 ; 
+
 
   Int_t nb   = GetGeometry()->GetNModules() ; 
   
   // create checkables 
   fQAHitsMul   = new AliPHOSQAIntCheckable("HitsM") ; 
   fQATotEner   = new AliPHOSQAFloatCheckable("TotEn") ; 
-  fQAHitsMulB  = new TClonesArray("AliPHOSQAIntCheckable",nb) ; 
+  fQAHitsMulB  = new TClonesArray("AliPHOSQAIntCheckable",nb) ;
+  fQAHitsMulB->SetOwner() ; 
   fQATotEnerB  = new TClonesArray("AliPHOSQAFloatCheckable", nb); 
+  fQATotEnerB->SetOwner() ; 
   char tempo[20]  ; 
   Int_t i ; 
   for ( i = 0 ; i < nb ; i++ ) {
@@ -121,29 +140,54 @@ AliPHOSv1::AliPHOSv1(const char *name, const char *title):
   fQAHitsMul->AddChecker(hmc) ; 
   fQATotEner->AddChecker(emc) ; 
   for ( i = 0 ; i < nb ; i++ ) {
-    ((AliPHOSQAIntCheckable*)(*fQAHitsMulB)[i])->AddChecker(bhmc) ;
-    ((AliPHOSQAFloatCheckable*)(*fQATotEnerB)[i])->AddChecker(bemc) ; 
+    (static_cast<AliPHOSQAIntCheckable*>((*fQAHitsMulB)[i]))->AddChecker(bhmc) ;
+    (static_cast<AliPHOSQAFloatCheckable*>((*fQATotEnerB)[i]))->AddChecker(bemc) ; 
   }
+
 }
 
 //____________________________________________________________________________
 AliPHOSv1::~AliPHOSv1()
 {
   // dtor
-
-  if ( fHits) {
+ if ( fHits) {
     fHits->Delete() ; 
     delete fHits ;
     fHits = 0 ; 
   }
+  
+  if ( fQAHitsMulB ) {
+    fQAHitsMulB->Delete() ;
+    delete fQAHitsMulB ; 
+  }
+
+  if ( fQATotEnerB ) {
+    fQATotEnerB->Delete() ;
+    delete fQATotEnerB ; 
+  }
+}
+
+//____________________________________________________________________________
+void AliPHOSv1::Copy(AliPHOSv1 & phos)
+{
+  TObject::Copy(phos) ; 
+  AliPHOSv0::Copy(phos) ; 
+  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)
 {
   // 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 ;
@@ -154,10 +198,11 @@ void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id,
   newHit = new AliPHOSHit(shunt, primary, tracknumber, 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 ;
     }
   }
@@ -165,14 +210,13 @@ void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id,
   if ( !deja ) {
     new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;
     // get the block Id number
-    Int_t * relid = new Int_t[geom->GetNModules()] ; 
+    Int_t relid[4] ;
     geom->AbsToRelNumbering(Id, relid) ;
     // and fill the relevant QA checkable (only if in PbW04)
     if ( relid[1] == 0 ) {
       fQAHitsMul->Update(1) ; 
-      ((AliPHOSQAIntCheckable*)(*fQAHitsMulB)[relid[0]-1])->Update(1) ;
+      (static_cast<AliPHOSQAIntCheckable*>((*fQAHitsMulB)[relid[0]-1]))->Update(1) ;
     } 
-    delete relid ; 
     fNhits++ ;
   }
 
@@ -196,6 +240,7 @@ 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
+  
 
   if ( fQATotEner ) { 
     if ( fQATotEner->HasChanged() ) {
@@ -207,8 +252,8 @@ void AliPHOSv1::FinishEvent()
   Int_t i ; 
   if ( fQAHitsMulB && fQATotEnerB ) {
     for (i = 0 ; i < GetGeometry()->GetNModules() ; i++) {
-      AliPHOSQAIntCheckable * ci = (AliPHOSQAIntCheckable*)(*fQAHitsMulB)[i] ;  
-      AliPHOSQAFloatCheckable* cf = (AliPHOSQAFloatCheckable*)(*fQATotEnerB)[i] ; 
+      AliPHOSQAIntCheckable * ci = static_cast<AliPHOSQAIntCheckable*>((*fQAHitsMulB)[i]) ;  
+      AliPHOSQAFloatCheckable* cf = static_cast<AliPHOSQAFloatCheckable*>((*fQATotEnerB)[i]) ; 
       if ( ci->HasChanged() ) { 
        ci->CheckMe() ;  
        ci->Reset() ;
@@ -228,192 +273,216 @@ void AliPHOSv1::FinishEvent()
       fQAHitsMul->Reset() ; 
     }
   } 
-}
 
+  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() ); 
+  Int_t tracknumber =  gAlice->GetMCApp()->GetCurrentTrackNumber() ; 
   TString name      =  GetGeometry()->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] += GetGeometry()->GetNModules() - GetGeometry()->GetNPPSDModules();
-       }
-               gMC->CurrentVolOffID(3, relid[1]) ;  // get the Micromegas Module number 
-      // 1-> GetGeometry()->GetNumberOfModulesPhi() * GetGeometry()->GetNumberOfModulesZ() upper
-      //   > GetGeometry()->GetNumberOfModulesPhi() * GetGeometry()->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
-
-               GetGeometry()->RelToAbsNumbering(relid, absid) ; 
-
-       // add current hit to the hit list      
-         AddHit(fIshunt, primary, tracknumber, absid, xyze);
-
-
-      } // 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
 
-      TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0);   // array of digits for current hit
-      CPVDigitize(pmom,xyd,moduleNumber,cpvDigits);
+      xyzte[3] = gMC->TrackTime() ;
+      xyzte[4] = cpvDigit->GetQpad() ;                          // amplitude in a pad
+      AddHit(fIshunt, -1, tracknumber, absid, xyzte);      // -1: No need in primary for CPV
       
-      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
-       GetGeometry()->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 
+         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 ( lostenergy != 0 ) {  // Track is inside the crystal and deposits some energy 
+      xyzte[3] = gMC->TrackTime() ;     
       
-      // fill the relevant QA Checkables
-      fQATotEner->Update( xyze[3] ) ;                                             // total energy in PHOS
-      ((AliPHOSQAFloatCheckable*)(*fQATotEnerB)[relid[0]-1])->Update( xyze[3] ) ; // energy in this block  
-
-      if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(10),"PHO1") == 0 )
-       relid[0] += GetGeometry()->GetNModules() - GetGeometry()->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
+      gMC->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ;
       
-      // get the absolute Id number
-      GetGeometry()->RelToAbsNumbering(relid, absid) ; 
+      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 =-1 ;
+      if(fIshunt == 1)
+        primary  =  gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() ); 
+      else if(fIshunt == 2){
+         primary = gAlice->GetMCApp()->GetCurrentTrackNumber() ;
+         TParticle * part = gAlice->GetMCApp()->Particle(primary) ;
+         while ( !part->TestBit(kKeepBit) ) {
+           primary = part->GetFirstMother() ;
+           if(primary == -1) 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) ;
+         }
+      }
 
+        
       // add current hit to the hit list
-       AddHit(fIshunt, primary,tracknumber, absid, xyze);
-
-
+      // Info("StepManager","%d %d", primary, tracknumber) ; 
+      AddHit(fIshunt, primary,tracknumber, absid, xyzte);
+      
+      // fill the relevant QA Checkables
+      fQATotEner->Update( xyzte[4] ) ;                                             // total energy in PHOS
+      (static_cast<AliPHOSQAFloatCheckable*>((*fQATotEnerB)[moduleNumber-1]))->Update( xyzte[4] ) ; // energy in this block  
+      
     } // 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:
@@ -446,7 +515,7 @@ 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 * GetGeometry()->GetCPVGasThickness();
   Float_t dXY   = pX/pNorm * GetGeometry()->GetCPVGasThickness();
@@ -519,7 +588,7 @@ void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, Int_t moduleNumbe
   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++) {