]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSv1.cxx
Inheritance from TObject. Automatic streamers.
[u/mrichter/AliRoot.git] / PHOS / AliPHOSv1.cxx
index a86a32406fc74bb383490d19ac769596a0323d4c..5975ba27d06aa300facbfb8740cc8ef6b0b73a16 100644 (file)
@@ -36,7 +36,7 @@
 #include "TNode.h"
 #include "TRandom.h"
 #include "TTree.h"
-
+#include "TParticle.h"
 
 // --- Standard library ---
 
@@ -54,6 +54,9 @@
 #include "AliConst.h"
 #include "AliMC.h"
 #include "AliPHOSGeometry.h"
+#include "AliPHOSQAIntCheckable.h"
+#include "AliPHOSQAFloatCheckable.h"
+#include "AliPHOSQAMeanChecker.h"
 
 ClassImp(AliPHOSv1)
 
@@ -61,18 +64,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) 
+ 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
@@ -85,14 +97,65 @@ 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->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) ; 
+  fQATotEnerB  = new TClonesArray("AliPHOSQAFloatCheckable", nb); 
+  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) ;
+  }
+
+  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) ; 
+  }
+
 }
 
 //____________________________________________________________________________
@@ -105,215 +168,281 @@ AliPHOSv1::~AliPHOSv1()
     delete fHits ;
     fHits = 0 ; 
   }
+  if (fTreeQA) 
+    delete fTreeQA ; 
 }
 
 //____________________________________________________________________________
 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) ;
 
   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) ;
+    // 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++ ;
   }
 
   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
+//   if ( fQAHitsMul ) 
+//     fQAHitsMul->Update( fHits->GetEntriesFast() ) ; 
+
+}
+
+//____________________________________________________________________________
+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() ) {
+      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() ;
+      }
+    } 
+  }
+  
+  // check the total multiplicity 
+  
+  if ( fQAHitsMul ) {
+    if ( fQAHitsMul->HasChanged() ) { 
+      fQAHitsMul->CheckMe() ; 
+      fQAHitsMul->Reset() ; 
+    }
+  } 
+}
 //____________________________________________________________________________
 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
+  TString name      =  GetGeometry()->GetName() ; 
 
-               fGeom->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 ;
+  
+  if( gMC->CurrentVolID(copy) == gMC->VolId("PCPQ") &&
+      (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,moduleNumber,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
+      primary = -1;                                             // No need in primary for CPV
+      AddHit(fIshunt, primary, tracknumber, 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) ;
-         }
-       }
-      }
-      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);
-
-       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 ) {  // 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();      
+  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.002 ||
+         xyzd[1] < -GetGeometry()->GetCrystalSize(1)/2+0.002) {
+       TParticle * part = 0 ; 
+       Int_t parent = gAlice->CurrentTrack() ; 
+       while ( parent != -1 ) {
+         part = gAlice->Particle(parent) ; 
+         part->SetBit(kKeepBit);
+         parent = part->GetFirstMother() ; 
+       }
+      }
+    }
+    if ( lostenergy != 0 ) {  // Track is inside the crystal and deposits some energy 
+      xyzte[3] = gMC->TrackTime() ;     
+      
+      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  ;
       
-      // get the absolute Id number
-      fGeom->RelToAbsNumbering(relid, absid) ; 
-
       // add current hit to the hit list
-       AddHit(fIshunt, primary,tracknumber, absid, xyze);
-
-
+      //cout << "AliPHOSv1::StepManager " << primary << " " << tracknumber << endl ; 
+      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
+  
 }
 
 //____________________________________________________________________________
@@ -329,7 +458,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
@@ -352,13 +481,13 @@ void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, Int_t moduleNumbe
 
 //    cout << "CPVDigitize: YVK : "<<hitX<<" "<<hitZ<<" | "<<pX<<" "<<pZ<<" "<<pNorm<<endl;
 
-  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 +547,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 +597,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));