]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSv1.cxx
ReaderESDtree, MUON analysis, reading MUON data froESD in ReaderESD (Christian FINCK)
[u/mrichter/AliRoot.git] / PHOS / AliPHOSv1.cxx
index bb4a802a6a7a282ee554ef8bc9a3306566ed76ad..1a7fec91aabe2c887ee8fd3747f7b59e1abf0437 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 "AliPHOSCPVDigit.h"
+#include "AliPHOSGeometry.h"
 #include "AliPHOSHit.h"
-#include "AliPHOSDigit.h"
-#include "AliPHOSReconstructioner.h"
+#include "AliPHOSQAFloatCheckable.h"
+#include "AliPHOSQAIntCheckable.h"
+#include "AliPHOSQAMeanChecker.h"
+#include "AliPHOSv1.h"
 #include "AliRun.h"
-#include "AliConst.h"
 #include "AliMC.h"
 
 ClassImp(AliPHOSv1)
 
 //____________________________________________________________________________
-AliPHOSv1::AliPHOSv1()
+AliPHOSv1::AliPHOSv1():
+AliPHOSv0()
 {
-  // ctor
-
-  // Create an empty array of AliPHOSCPVModule to satisfy
-  // AliPHOSv1::Streamer when reading root file
-  fReconstructioner  = 0;
-  fTrackSegmentMaker = 0;
-
-  fHits = new TClonesArray("AliPHOSHit",1000) ;
-  //  if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
-  //    Error("AliPHOSv1","Can not create array of EMC modules");
-  //    exit(1);
-  //  }
-
-  //  if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
-  //    Error("AliPHOSv1","Can not create array of CPV modules");
-  //    exit(1);
-  //  }
+  // 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) 
+ 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
@@ -97,94 +78,71 @@ AliPHOSv0(name,title)
   //     (this array is reset after each primary has been tracked).
   //
 
-  fPinElectronicNoise = 0.010 ;
-  fDigitThreshold      = 0.01 ;   // 1 GeV 
-  fDigitizeA= 0. ;             
-  fDigitizeB = 10000000. ;    
 
 
   // 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 ;
 
-  fReconstructioner  = 0;
-  fTrackSegmentMaker = 0;
-
-  fIshunt     =  1 ; // All hits are associated with primary particles
-  
-  // Create array of EMC modules of the size of PHOS modules number
+  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() ; 
   
-//  if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNModules())) ) {
-//    Error("AliPHOSv1","Can not create array of EMC modules");
-//    exit(1);
-//  }
-//  TClonesArray &lemcmodule = *fEMCModules;
-//  for (Int_t i=0; i<fGeom->GetNModules(); i++) new(lemcmodule[i]) AliPHOSCPVModule();
-
-  // Create array of CPV modules for the IHEP's version of CPV
-
-//  if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
-//    // Create array of CPV modules of the size of PHOS modules number
-
-//    if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNCPVModules())) ) {
-//      Error("AliPHOSv1","Can not create array of CPV modules");
-//      exit(1);
-//    }
-//    TClonesArray &lcpvmodule = *fCPVModules;
-//    for (Int_t i=0; i<fGeom->GetNCPVModules(); i++) new(lcpvmodule[i]) AliPHOSCPVModule();
-//  }
-//  else {
-//    // Create an empty array of AliPHOSCPVModule to satisfy
-//    // AliPHOSv1::Streamer when writing root file
-
-//    fCPVModules=new TClonesArray("AliPHOSCPVModule",0);
-
-//  }
-}
-
-//____________________________________________________________________________
-AliPHOSv1::AliPHOSv1(AliPHOSReconstructioner * Reconstructioner, const char *name, const char *title):
-  AliPHOSv0(name,title)
-{
-  // ctor : title is used to identify the layout
-  //        GPS2 = 5 modules (EMC + PPSD)   
-  // We use 2 arrays of hits :
-  //
-  //   - fHits (the "normal" one), which retains the hits associated with
-  //     the current primary particle being tracked
-  //     (this array is reset after each primary has been tracked).
-  //
-  //   - fTmpHits, which retains all the hits of the current event. It 
-  //     is used for the digitization part.
-
-  fPinElectronicNoise = 0.010 ;
-
-  // We do not want to save in TreeH the raw hits
-  //fHits   = new TClonesArray("AliPHOSHit",100) ;
-
-  fDigits = 0 ;
-  fHits= new TClonesArray("AliPHOSHit",1000) ;
-
-  fNhits = 0 ;
-
-  fIshunt     =  1 ; // All hits are associated with primary particles
-  // gets an instance of the geometry parameters class  
-  fGeom =  AliPHOSGeometry::GetInstance(title, "") ; 
-
-  if (fGeom->IsInitialized() ) 
-    cout << "AliPHOS" << Version() << " : PHOS geometry intialized for " << fGeom->GetName() << endl ;
-  else
-    cout << "AliPHOS" << Version() << " : PHOS geometry initialization failed !" << endl ;   
+  // create checkables 
+  fQAHitsMul   = new AliPHOSQAIntCheckable("HitsM") ; 
+  fQATotEner   = new AliPHOSQAFloatCheckable("TotEn") ; 
+  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++ ) {
+    sprintf(tempo, "HitsMB%d", i+1) ; 
+    new( (*fQAHitsMulB)[i]) AliPHOSQAIntCheckable(tempo) ; 
+    sprintf(tempo, "TotEnB%d", i+1) ; 
+    new( (*fQATotEnerB)[i] ) AliPHOSQAFloatCheckable(tempo) ;
+  }
 
-  // Defining the PHOS Reconstructioner
- fReconstructioner = Reconstructioner ;
+  AliPHOSQAMeanChecker * hmc  = new AliPHOSQAMeanChecker("HitsMul", 100. ,25.) ; 
+  AliPHOSQAMeanChecker * emc  = new AliPHOSQAMeanChecker("TotEner", 10. ,5.) ; 
+  AliPHOSQAMeanChecker * bhmc = new AliPHOSQAMeanChecker("HitsMulB", 100. ,5.) ; 
+  AliPHOSQAMeanChecker * bemc = new AliPHOSQAMeanChecker("TotEnerB", 2. ,.5) ; 
+
+  // associate checkables and checkers 
+  fQAHitsMul->AddChecker(hmc) ; 
+  fQATotEner->AddChecker(emc) ; 
+  for ( i = 0 ; i < nb ; i++ ) {
+    (static_cast<AliPHOSQAIntCheckable*>((*fQAHitsMulB)[i]))->AddChecker(bhmc) ;
+    (static_cast<AliPHOSQAFloatCheckable*>((*fQATotEnerB)[i]))->AddChecker(bemc) ; 
+  }
 
 }
 
@@ -192,69 +150,73 @@ AliPHOSv1::AliPHOSv1(AliPHOSReconstructioner * Reconstructioner, const char *nam
 AliPHOSv1::~AliPHOSv1()
 {
   // dtor
-
-  if ( fHits) {
+ if ( fHits) {
     fHits->Delete() ; 
     delete fHits ;
     fHits = 0 ; 
   }
-
-  if ( fSDigits) {
-    fSDigits->Delete() ; 
-    delete fSDigits ;
-    fSDigits = 0 ; 
+  
+  if ( fQAHitsMulB ) {
+    fQAHitsMulB->Delete() ;
+    delete fQAHitsMulB ; 
   }
 
-  if ( fDigits) {
-    fDigits->Delete() ; 
-    delete fDigits ;
-    fDigits = 0 ; 
+  if ( fQATotEnerB ) {
+    fQATotEnerB->Delete() ;
+    delete fQATotEnerB ; 
   }
+}
 
-  if ( fEmcRecPoints ) {
-    fEmcRecPoints->Delete() ; 
-    delete fEmcRecPoints ; 
-    fEmcRecPoints = 0 ; 
-  }
-  
-  if ( fPpsdRecPoints ) { 
-    fPpsdRecPoints->Delete() ;
-    delete fPpsdRecPoints ;
-    fPpsdRecPoints = 0 ; 
-  }
-    
-  if ( fTrackSegments ) {
-    fTrackSegments->Delete() ; 
-    delete fTrackSegments ;
-    fTrackSegments = 0 ; 
-  }
-  
+//____________________________________________________________________________
+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, Int_t trackpid, TLorentzVector p, Float_t * lpos)
+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 ;
   AliPHOSHit *curHit ;
   Bool_t deja = kFALSE ;
+  AliPHOSGeometry * geom = GetGeometry() ; 
 
-  newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits, trackpid, p, lpos) ;
+  newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits) ;
 
   for ( hitCounter = fNhits-1 ; hitCounter >= 0 && !deja ; hitCounter-- ) {
-    curHit = (AliPHOSHit*) (*fHits)[hitCounter] ;
+    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) ;
+    // and fill the relevant QA checkable (only if in PbW04)
+    if ( relid[1] == 0 ) {
+      fQAHitsMul->Update(1) ; 
+      (static_cast<AliPHOSQAIntCheckable*>((*fQAHitsMulB)[relid[0]-1]))->Update(1) ;
+    } 
     fNhits++ ;
   }
 
@@ -262,552 +224,257 @@ void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id,
 }
 
 //____________________________________________________________________________
-void AliPHOSv1::Hits2SDigits(){
-  //Collects all hits in the same active volume into digit
-
-  Int_t i ;
-  Int_t j ; 
-  AliPHOSHit   * hit ;
-  AliPHOSDigit * newdigit ;
-  AliPHOSDigit * curdigit ;
-  Bool_t deja = kFALSE ; 
-  
+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
+//   if ( fQAHitsMul ) 
+//     fQAHitsMul->Update( fHits->GetEntriesFast() ) ; 
 
-  Int_t itrack ;
-  for (itrack=0; itrack<gAlice->GetNtrack(); itrack++){
-        
-    //=========== Get the Hits Tree for the Primary track itrack
-    gAlice->ResetHits();    
-    gAlice->TreeH()->GetEvent(itrack);
-      
+}
 
-    for ( i = 0 ; i < fHits->GetEntries() ; i++ ) {
-      hit = (AliPHOSHit*)fHits->At(i) ;
-    
-      // Assign primary number only if contribution is significant
-      if( hit->GetEnergy() > fDigitThreshold)
-       newdigit = new AliPHOSDigit( hit->GetPrimary(), hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
-      else
-       newdigit = new AliPHOSDigit( -1               , hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
-      deja =kFALSE ;
-
-      for ( j = 0 ; j < fnSdigits ;  j++) { 
-       curdigit = (AliPHOSDigit*) fSDigits->At(j) ;
-       if ( *curdigit == *newdigit) {
-         *curdigit = *curdigit + *newdigit ; 
-         deja = kTRUE ; 
-       }
-      }
+//____________________________________________________________________________
+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 ( !deja ) {
-       new((*fSDigits)[fnSdigits]) AliPHOSDigit(* newdigit) ;
-       fnSdigits++ ;  
+  if ( fQATotEner ) { 
+    if ( fQATotEner->HasChanged() ) {
+      fQATotEner->CheckMe() ; 
+      fQATotEner->Reset() ; 
+    }
+  }
+  
+  Int_t i ; 
+  if ( fQAHitsMulB && fQATotEnerB ) {
+    for (i = 0 ; i < GetGeometry()->GetNModules() ; i++) {
+      AliPHOSQAIntCheckable * ci = static_cast<AliPHOSQAIntCheckable*>((*fQAHitsMulB)[i]) ;  
+      AliPHOSQAFloatCheckable* cf = static_cast<AliPHOSQAFloatCheckable*>((*fQATotEnerB)[i]) ; 
+      if ( ci->HasChanged() ) { 
+       ci->CheckMe() ;  
+       ci->Reset() ;
+      } 
+      if ( cf->HasChanged() ) { 
+       cf->CheckMe() ; 
+       cf->Reset() ;
       }
-      
-      delete newdigit ;    
     } 
-
-  } // loop over tracks
-
-  fSDigits->Sort() ;
-
-  fnSdigits = fSDigits->GetEntries() ;
-  fSDigits->Expand(fnSdigits) ;
-
-  for (i = 0 ; i < fnSdigits ; i++) { 
-    AliPHOSDigit * digit = (AliPHOSDigit *) fSDigits->At(i) ; 
-    digit->SetIndexInList(i) ;     
   }
-
-  gAlice->TreeS()->Fill() ;
   
-  TTree *ts = gAlice->TreeS();
-  TDirectory *dir = ts->GetDirectory() ; 
-  cout << "dir name is " << dir->GetName() << endl ;  
-  //  ts->Write(0,TObject::kOverwrite) ;
-
-gAlice->TreeS()->GetBranch("PHOS")->Fill();
-gAlice->TreeS()->GetBranch("PHOS")->Write();
-
+  // check the total multiplicity 
+  
+  if ( fQAHitsMul ) {
+    if ( fQAHitsMul->HasChanged() ) { 
+      fQAHitsMul->CheckMe() ; 
+      fQAHitsMul->Reset() ; 
+    }
+  } 
 
+  AliDetector::FinishEvent(); 
 }
 //____________________________________________________________________________
-void AliPHOSv1::SDigits2Digits(){
-  //Adds noise to the summable digits and removes everething below thresholds
-  //Note, that sDigits should be SORTED in accordance with abs ID.
+void AliPHOSv1::StepManager(void)
+{
+   // 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        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     ;
 
-  gAlice->TreeS()->GetEvent(0) ;
+  Int_t tracknumber =  gAlice->GetMCApp()->GetCurrentTrackNumber() ; 
+  TString name      =  GetGeometry()->GetName() ; 
 
-  cout << "fSdigits " <<  fSDigits << " " << fSDigits->GetEntries() << endl ;
+  Int_t moduleNumber ;
   
-
-  // Noise induced by the PIN diode of the PbWO crystals
-  Int_t iCurSDigit = 0 ;
-  //we assume, that there is al least one EMC digit...
-  Int_t idCurSDigit = ((AliPHOSDigit *)fSDigits->At(0))->GetId() ;
-
-  cout << "fDigits " << fDigits << " " <<idCurSDigit<<  endl ;
-
-  Int_t absID ;
-  for(absID = 1; absID < fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ(); absID++){
+  if( gMC->CurrentVolID(copy) == gMC->VolId("PCPQ") &&
+      (gMC->IsTrackEntering() ) &&
+      gMC->TrackCharge() != 0) {      
     
-    cout << "absID " << absID << " " << idCurSDigit << endl ; 
-
-    Float_t noise = gRandom->Gaus(0., fPinElectronicNoise) ; 
-    if(absID < idCurSDigit ){ 
-      cout << "In < idC  " << noise << " " << fDigitThreshold << endl ;
-      if(noise >fDigitThreshold ){
-       cout << "noise " << absID << " " << noise << endl;
-       new((*fDigits)[fNdigits]) AliPHOSDigit( -1,absID,Digitize(noise) ) ;
-       cout << "FN " << fNdigits << endl ;
-       fNdigits++ ;  
-      }
-    }
-    else{ //add noise and may be remove the true hit
-      cout << "correcting digit " << iCurSDigit << endl ;
-      Float_t signal = noise + Calibrate(((AliPHOSDigit *)fSDigits->At(iCurSDigit))->GetAmp()) ;
-      cout << "signal " << signal << endl ;
-      if( signal >fDigitThreshold ){
-       cout << "signal " << signal << endl ;
-       AliPHOSDigit * digit = (AliPHOSDigit*) fSDigits->At(iCurSDigit) ;
-       new((*fDigits)[fNdigits]) AliPHOSDigit( *digit ) ;
-       ((AliPHOSDigit *)fDigits->At(fNdigits))->SetAmp(Digitize(signal));
-       cout << "fNdigits " << fNdigits << endl ;
-       fNdigits++ ;  
-      }
-
-      if(iCurSDigit < fSDigits->GetEntries()-1){
-       iCurSDigit++ ;
-       idCurSDigit = ((AliPHOSDigit*)fSDigits->At(iCurSDigit))->GetId() ;
-      }
-      else
-       idCurSDigit = 10000000; //no real hits left    
-    }
+    gMC -> TrackPosition(pos);
     
-  }  
-
-  //remove PPSD/CPV digits below thresholds
-  Int_t idigit ;
-  for(idigit = iCurSDigit; idigit < fSDigits->GetEntries() ; idigit++){  //loop over CPV/PPSD digits
+    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
     
-    AliPHOSDigit * digit = (AliPHOSDigit *) fSDigits->At(idigit) ; 
-    Float_t ene = Calibrate(digit->GetAmp()) ;
+    Float_t        xyd[3]={0,0,0}   ;   //local position of the entering
+    xyd[0]  = xyzd[0];
+    xyd[1]  =-xyzd[2];
+    xyd[2]  =-xyzd[1];
     
-    Int_t relid[4] ; 
-    fGeom->AbsToRelNumbering(digit->GetId(), relid) ; 
-    if ( relid[0] > fGeom->GetNCPVModules() ){ //ppsd
-      if ( ( (relid[1] > 0) && (ene > fPpsdEnergyThreshold)) ||    //PPSD digit
-          ( (relid[1] < 0) && (ene > fCpvEnergyThreshold ) ) )    //CPV digit 
-       new((*fDigits)[fNdigits]) AliPHOSDigit( *digit ) ;
-       fNdigits++ ;
-    }
-  }    
+    // 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];
     
-  fDigits->Compress() ;  
-  
-  fNdigits = fDigits->GetEntries() ;
-  fDigits->Expand(fNdigits) ;
-
-  Int_t i ;
-  for (i = 0 ; i < fNdigits ; i++) { 
-    AliPHOSDigit * digit = (AliPHOSDigit *) fDigits->At(i) ; 
-    digit->SetIndexInList(i) ;     
-  }
-
-
-  gAlice->TreeD()->Fill() ;
-
-  gAlice->TreeD()->Write(0,TObject::kOverwrite) ;
-}
-
-//___________________________________________________________________________
-void AliPHOSv1::MakeBranch(Option_t* opt, char *file)
-{ 
-
-
-  char *cH ; 
-  // Create new branche in the current Root Tree in the digit Tree
-  AliDetector::MakeBranch(opt) ;
+    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];
 
-
-  cH = strstr(opt,"S");
-  //Create a branch for SDigits
-  if( cH ){
-    char branchname[20];
-    sprintf(branchname,"%s",GetName());  
-    if(fSDigits)
-      fSDigits->Clear();
-    else
-      fSDigits = new TClonesArray("AliPHOSDigit",1000);
-    fnSdigits = 0 ;
-    cout << " AliPHOSv1::MakeBranch : " << file << endl ; 
-    gAlice->MakeBranchInTree(gAlice->TreeS(),branchname,&fSDigits,fBufferSize,file);  
-  }
-  
-  cH = strstr(opt,"D");
-  //Create a branch for Digits
-  if( cH ){
-    char branchname[20];
-    sprintf(branchname,"%s",GetName());  
-    if(fSDigits)
-      fDigits->Clear();
-    else
-      fDigits = new TClonesArray("AliPHOSDigit",1000);
-    fNdigits = 0 ;
-    
-    gAlice->MakeBranchInTree(gAlice->TreeD(),branchname,&fSDigits,fBufferSize,file);  
-  }
-
-  cH = strstr(opt,"R");
-  //Create a branch for Reconstruction
-  if( cH ){
-    char branchname[20];
-
-    Int_t splitlevel = 0 ; 
+    // Digitize the current CPV hit:
     
-    fEmcRecPoints->Delete() ; 
-
-    if ( fEmcRecPoints && gAlice->TreeR() ) {
-      sprintf(branchname,"%sEmcRP",GetName()) ;
-      gAlice->TreeR()->Branch(branchname, "TObjArray", &fEmcRecPoints, fBufferSize, splitlevel) ; 
-    }
-
-    fPpsdRecPoints->Delete() ; 
-
-    if ( fPpsdRecPoints && gAlice->TreeR() ) {
-      sprintf(branchname,"%sPpsdRP",GetName()) ;
-      gAlice->TreeR()->Branch(branchname, "TObjArray", &fPpsdRecPoints, fBufferSize, splitlevel) ;
-    }
-
-    fTrackSegments->Delete() ; 
+    // 1. find pad response and    
+    gMC->CurrentVolOffID(3,moduleNumber);
+    moduleNumber--;
     
-    if ( fTrackSegments && gAlice->TreeR() ) { 
-      sprintf(branchname,"%sTS",GetName()) ;
-      gAlice->TreeR()->Branch(branchname, &fTrackSegments, fBufferSize) ;
-    }
+    TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0);   // array of digits for current hit
+    CPVDigitize(pmom,xyd,cpvDigits);
+      
+    Float_t xmean = 0;
+    Float_t zmean = 0;
+    Float_t qsum  = 0;
+    Int_t   idigit,ndigits;
     
-    fRecParticles->Delete() ; 
+    // 2. go through the current digit list and sum digits in pads
     
-    if ( fRecParticles && gAlice->TreeR() ) { 
-      sprintf(branchname,"%sRP",GetName()) ;
-      gAlice->TreeR()->Branch(branchname, &fRecParticles, fBufferSize) ;
+    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() ;
     
-  }
-
-
-
-}
-
-//_____________________________________________________________________________
-void AliPHOSv1::Reconstruction(AliPHOSReconstructioner * Reconstructioner)
-{ 
-  // 1. Reinitializes the existing RecPoint, TrackSegment, and RecParticles Lists and 
-  // 2. Creates TreeR with a branch for each list
-  // 3. Steers the reconstruction processes
-  // 4. Saves the 3 lists in TreeR
-  // 5. Write the Tree to File
-  
-  fReconstructioner = Reconstructioner ;
+    // 3. add digits to temporary hit list fTmpHits
     
-  // 1.
-
-  //  gAlice->MakeTree("R") ; 
-  
-  MakeBranch("R") ;
-  
-  // 3.
-
-  fReconstructioner->Make(fDigits, fEmcRecPoints, fPpsdRecPoints, fTrackSegments, fRecParticles);
-
-  printf("Reconstruction: %d %d %d %d\n",
-        fEmcRecPoints->GetEntries(),fPpsdRecPoints->GetEntries(),
-        fTrackSegments->GetEntries(),fRecParticles->GetEntries());
-
-  // 4. Expand or Shrink the arrays to the proper size
-  
-  Int_t size ;
-  
-  size = fEmcRecPoints->GetEntries() ;
-  fEmcRecPoints->Expand(size) ;
-
-  size = fPpsdRecPoints->GetEntries() ;
-  fPpsdRecPoints->Expand(size) ;
-
-  size = fTrackSegments->GetEntries() ;
-  fTrackSegments->Expand(size) ;
-
-  size = fRecParticles->GetEntries() ;
-  fRecParticles->Expand(size) ;
-
-  gAlice->TreeR()->Fill() ;
-  // 5.
-
-  gAlice->TreeR()->Write(0,TObject::kOverwrite) ;
-  // Deleting reconstructed objects
-  ResetReconstruction();
-  
-}
-
-//____________________________________________________________________________
-void AliPHOSv1::ResetHits() 
-{              
-  // Reset hit tree for EMC and CPV
-  // Yuri Kharlov, 28 September 2000
-
-  AliDetector::ResetHits();
-  //  for (Int_t i=0; i<fGeom->GetNModules(); i++) ((AliPHOSCPVModule*)(*fEMCModules)[i]) -> Clear();
-  //  if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
-  //    for (Int_t i=0; i<fGeom->GetNCPVModules(); i++) ((AliPHOSCPVModule*)(*fCPVModules)[i]) -> Clear();
-  //  }
-  
-}  
-//____________________________________________________________________________
-void AliPHOSv1::ResetReconstruction() 
-{ 
-  // Deleting reconstructed objects
-
-  if ( fEmcRecPoints  )  fEmcRecPoints ->Delete();
-  if ( fPpsdRecPoints )  fPpsdRecPoints->Delete();
-  if ( fTrackSegments )  fTrackSegments->Delete();
-  if ( fRecParticles  )  fRecParticles ->Delete();
-  
-}
-
-//____________________________________________________________________________
-//void AliPHOSv1::SDigits2Digits() {
-//  // Adds the noise to the summable digits and keeps digits above a threshold 
-//  // To make a digit
-//}
-
-//____________________________________________________________________________
-void AliPHOSv1::SetTreeAddress()
-{ 
-  //  TBranch *branch;
-  AliPHOS::SetTreeAddress();
-
-//  //Branch address for TreeR: RecPpsdRecPoint
-//   TTree *treeR = gAlice->TreeR();
-//   if ( treeR && fPpsdRecPoints ) {
-//     branch = treeR->GetBranch("PHOSPpsdRP");
-//     if (branch) branch->SetAddress(&fPpsdRecPoints) ;
-//  }
-
-  // Set branch address for the Hits Tree for hits in EMC modules
-  // Yuri Kharlov, 23 November 2000.
-
-  //  for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetEMCModule(i).SetTreeAddress("EMC",i+1);
-
-  // Set branch address for the Hits Tree for hits in CPV modules for IHEP geometry
-  // Yuri Kharlov, 28 September 2000.
-
-  //  if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
-  //    for( Int_t i=0; i<fGeom->GetNCPVModules(); i++ ) GetCPVModule(i).SetTreeAddress("CPV",i+1);
-  //  }
-
-}
-
-//____________________________________________________________________________
-
-void AliPHOSv1::StepManager(void)
-{
-  // Accumulates hits as long as the track stays in a single crystal or PPSD gas Cell
-
-//    if (gMC->IsTrackEntering())
-//      cout << "Track enters the volume " << gMC->CurrentVolName() << endl;
-//    if (gMC->IsTrackExiting())
-//      cout << "Track leaves the volume " << gMC->CurrentVolName() << endl;
-
-  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
-  TLorentzVector pos      ;      // Lorentz vector of the track current position
-  TLorentzVector pmom     ;      //momentum of the particle initiated hit
-  Float_t        xyd[3]   ;      //local posiiton of the entering
-  Bool_t         entered = kFALSE    ;  
-  Int_t          copy     ;
-
-  Int_t tracknumber =  gAlice->CurrentTrack() ; 
-  Int_t primary     =  gAlice->GetPrimary( gAlice->CurrentTrack() ); 
-  TString name      =  fGeom->GetName() ; 
-  Int_t trackpid    =  0  ; 
-
-  if( gMC->IsTrackEntering() ){ // create hit with position and momentum of new particle, 
-                                // but may be without energy deposition
-
-    // Current position of the hit in the local ref. system
-      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
-      xyd[0]  = xyzd[0];
-      xyd[1]  =-xyzd[1];
-      xyd[2]  =-xyzd[2];
-
+    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
-      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];
-
-      trackpid = gMC->TrackPid();
-      entered = kTRUE ;      // Mark to create hit even withou energy deposition
-
-  }
-
-
-  if ( name == "GPS2" || name == "MIXT" ) {            // ======> CPV is a GPS' PPSD
-
-    if( gMC->CurrentVolID(copy) == gMC->VolId("GCEL") ) // 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) || entered ) { // there is deposited energy or new particle entering  PPSD
-               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, trackpid, pmom, xyd);
-
-
-      } // 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("CPVQ") &&
-       entered &&
-       gMC->TrackCharge() != 0) {      
-      
-      // 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) ;
-         }
-       }
-      }
-      cpvDigits->Compress() ;
-
-      // 3. add digits to temporary hit list fTmpHits
-
-      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, trackpid, pmom, xyd);
-
-       if (cpvDigit->GetQpad() > 0.02) {
-         xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
-         zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
-         qsum  += cpvDigit->GetQpad();
-       }
+      if (cpvDigit->GetQpad() > 0.02) {
+       xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
+       zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
+       qsum  += cpvDigit->GetQpad();
       }
+    }
+    if (cpvDigits) {
+      cpvDigits->Delete();
       delete cpvDigits;
+      cpvDigits=0;
     }
-  } // 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) || entered ) {  // Track is inside the crystal and deposits some energy or just entered 
-
-      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();      
+  if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) { //  We are inside a PBWO crystal
 
-      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->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  
+       TParticle * part = 0 ; 
+       Int_t parent = gAlice->GetMCApp()->GetCurrentTrackNumber() ; 
+       while ( parent != -1 ) {
+         part = gAlice->GetMCApp()->Particle(parent) ; 
+         part->SetBit(kKeepBit);
+         parent = part->GetFirstMother() ; 
+       }
+      }
+    }
+    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 =-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, trackpid,pmom, xyd);
-
-
+      // 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:
@@ -819,7 +486,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
@@ -840,15 +507,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;
 
@@ -908,20 +575,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;
@@ -958,10 +625,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));