]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSv1.cxx
update from Marco
[u/mrichter/AliRoot.git] / PHOS / AliPHOSv1.cxx
index 47c0ebebc26de2f99e52946bc56378a91d340af4..a7b1d29b5f80909103a1d045ace54e9f93d8d920 100644 (file)
 
 /* $Id$ */
 
+/* History of cvs commits:
+ *
+ * $Log$
+ * Revision 1.111  2007/07/24 09:41:19  morsch
+ * AliStack included for kKeepBit.
+ *
+ * Revision 1.110  2007/03/10 08:58:52  kharlov
+ * Protection for noCPV geometry
+ *
+ * Revision 1.109  2007/03/01 11:37:37  kharlov
+ * Strip units changed from 8x1 to 8x2 (T.Pocheptsov)
+ *
+ * Revision 1.108  2007/02/02 09:40:50  alibrary
+ * Includes required by ROOT head
+ *
+ * Revision 1.107  2007/02/01 10:34:47  hristov
+ * Removing warnings on Solaris x86
+ *
+ * Revision 1.106  2006/11/14 17:11:15  hristov
+ * Removing inheritances from TAttLine, TAttMarker and AliRndm in AliModule. The copy constructor and assignment operators are moved to the private part of the class and not implemented. The corresponding changes are propagated to the detectors
+ *
+ * Revision 1.105  2006/09/13 07:31:01  kharlov
+ * Effective C++ corrections (T.Pocheptsov)
+ *
+ * Revision 1.104  2005/05/28 14:19:05  schutz
+ * Compilation warnings fixed by T.P.
+ *
+ */
+
 //_________________________________________________________________________
 // Implementation version v1 of PHOS Manager class 
 //---
-// Layout EMC + PPSD has name GPS2:
-// Produces cumulated hits
 //---
 // Layout EMC + CPV  has name IHEP:
 // Produces hits for CPV, cumulated hits
 //---
-// Layout EMC + CPV + PPSD has name GPS:
-// Produces hits for CPV, cumulated hits
 //---
 //*-- Author: Yves Schutz (SUBATECH)
 
 
 // --- ROOT system ---
-
-#include "TBRIK.h"
-#include "TNode.h"
-#include "TRandom.h"
-#include "TTree.h"
-#include "TParticle.h"
+#include <TClonesArray.h>
+#include <TParticle.h>
+#include <TVirtualMC.h>
 
 // --- Standard library ---
 
-#include <string.h>
-#include <stdlib.h>
-#include <strstream.h>
 
 // --- AliRoot header files ---
-
-#include "AliPHOSv1.h"
-#include "AliPHOSHit.h"
 #include "AliPHOSCPVDigit.h"
+#include "AliPHOSGeometry.h"
+#include "AliPHOSHit.h"
+#include "AliPHOSv1.h"
 #include "AliRun.h"
-#include "AliConst.h"
 #include "AliMC.h"
-#include "AliPHOSGeometry.h"
-#include "AliPHOSQAIntCheckable.h"
-#include "AliPHOSQAFloatCheckable.h"
-#include "AliPHOSQAMeanChecker.h"
+#include "AliStack.h"
+#include "AliPHOSSimParam.h"
 
 ClassImp(AliPHOSv1)
 
 //____________________________________________________________________________
-AliPHOSv1::AliPHOSv1():
-AliPHOSv0()
+AliPHOSv1::AliPHOSv1() : fCPVDigits("AliPHOSCPVDigit",20)
 {
-  // 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. ; 
-
+  //Def ctor.
 }
 
 //____________________________________________________________________________
 AliPHOSv1::AliPHOSv1(const char *name, const char *title):
- AliPHOSv0(name,title) 
+  AliPHOSv0(name,title), fCPVDigits("AliPHOSCPVDigit",20)
 {
   //
   // We store hits :
@@ -99,91 +100,27 @@ 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->AddHitList(fHits) ; 
+//   fCPVDigits("AliPHOSCPVDigit",20);
+  gAlice->GetMCApp()->AddHitList(fHits) ; 
 
   fNhits = 0 ;
 
   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->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) ;
-  }
-
-  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) ; 
-  }
-
 }
 
 //____________________________________________________________________________
 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::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits)
+void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t Id, Float_t * hits)
 {
   // Add a hit to the hit list.
   // A PHOS hit is the sum of all hits in a single crystal from one primary and within some time gate
@@ -194,10 +131,10 @@ void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id,
   Bool_t deja = kFALSE ;
   AliPHOSGeometry * geom = GetGeometry() ; 
 
-  newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits) ;
+  newHit = new AliPHOSHit(shunt, primary, Id, hits) ;
 
   for ( hitCounter = fNhits-1 ; hitCounter >= 0 && !deja ; hitCounter-- ) {
-    curHit = dynamic_cast<AliPHOSHit*>((*fHits)[hitCounter]) ;
+    curHit = static_cast<AliPHOSHit*>((*fHits)[hitCounter]) ;
     if(curHit->GetPrimary() != primary) break ; 
            // We add hits with the same primary, while GEANT treats primaries succesively 
     if( *curHit == *newHit ) {
@@ -211,14 +148,10 @@ void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id,
     // 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;
 }
 
@@ -228,8 +161,6 @@ 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() ) ; 
 
 }
 
@@ -239,38 +170,8 @@ 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() ; 
-    }
-  } 
+  AliDetector::FinishEvent(); 
 }
 //____________________________________________________________________________
 void AliPHOSv1::StepManager(void)
@@ -283,23 +184,24 @@ void AliPHOSv1::StepManager(void)
   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      =  GetGeometry()->GetName() ; 
-
   Int_t moduleNumber ;
   
-  if( gMC->CurrentVolID(copy) == gMC->VolId("PCPQ") &&
-      (gMC->IsTrackEntering() ) &&
-      gMC->TrackCharge() != 0) {      
-    
-    gMC -> TrackPosition(pos);
+  static Int_t idPCPQ = -1;
+  if (strstr(fTitle.Data(),"noCPV") == 0) 
+    idPCPQ = TVirtualMC::GetMC()->VolId("PCPQ");
+
+  if( TVirtualMC::GetMC()->CurrentVolID(copy) == idPCPQ &&
+     (TVirtualMC::GetMC()->IsTrackEntering() ) &&
+      TVirtualMC::GetMC()->TrackCharge() != 0) {      
     
+    TVirtualMC::GetMC() -> 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
-    
+    TVirtualMC::GetMC() -> 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];
@@ -307,12 +209,12 @@ void AliPHOSv1::StepManager(void)
     
     // Current momentum of the hit's track in the local ref. system
     TLorentzVector pmom     ;        //momentum of the particle initiated hit
-    gMC -> TrackMomentum(pmom);
+    TVirtualMC::GetMC() -> 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
+    TVirtualMC::GetMC() -> Gmtod (pm, pd, 2);        // transform 3-momentum from master to daughter system
     pmom[0] = pd[0];
     pmom[1] =-pd[1];
     pmom[2] =-pd[2];
@@ -320,11 +222,11 @@ void AliPHOSv1::StepManager(void)
     // Digitize the current CPV hit:
     
     // 1. find pad response and    
-    gMC->CurrentVolOffID(3,moduleNumber);
+    TVirtualMC::GetMC()->CurrentVolOffID(3,moduleNumber);
     moduleNumber--;
     
-    TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0);   // array of digits for current hit
-    CPVDigitize(pmom,xyd,moduleNumber,cpvDigits);
+//     TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0);   // array of digits for current hit
+    CPVDigitize(pmom,xyd,&fCPVDigits);
       
     Float_t xmean = 0;
     Float_t zmean = 0;
@@ -333,29 +235,29 @@ void AliPHOSv1::StepManager(void)
     
     // 2. go through the current digit list and sum digits in pads
     
-    ndigits = cpvDigits->GetEntriesFast();
+    ndigits = fCPVDigits.GetEntriesFast();
     for (idigit=0; idigit<ndigits-1; idigit++) {
-      AliPHOSCPVDigit  *cpvDigit1 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
+      AliPHOSCPVDigit  *cpvDigit1 = static_cast<AliPHOSCPVDigit*>(fCPVDigits.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));
+       AliPHOSCPVDigit  *cpvDigit2 = static_cast<AliPHOSCPVDigit*>(fCPVDigits.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) ;
+         Float_t qsumpad = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
+         cpvDigit2->SetQpad(qsumpad) ;
+         fCPVDigits.RemoveAt(idigit) ;
        }
       }
     }
-    cpvDigits->Compress() ;
+    fCPVDigits.Compress() ;
     
     // 3. add digits to temporary hit list fTmpHits
     
-    ndigits = cpvDigits->GetEntriesFast();
+    ndigits = fCPVDigits.GetEntriesFast();
     for (idigit=0; idigit<ndigits; idigit++) {
-      AliPHOSCPVDigit  *cpvDigit = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
+      AliPHOSCPVDigit  *cpvDigit = static_cast<AliPHOSCPVDigit*>(fCPVDigits.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
@@ -366,10 +268,11 @@ void AliPHOSv1::StepManager(void)
       
       // add current digit to the temporary hit list
 
-      xyzte[3] = gMC->TrackTime() ;
+      xyzte[3] = TVirtualMC::GetMC()->TrackTime() ;
       xyzte[4] = cpvDigit->GetQpad() ;                          // amplitude in a pad
-      primary = -1;                                             // No need in primary for CPV
-      AddHit(fIshunt, primary, tracknumber, absid, xyzte);
+
+      Int_t primary  =  gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() ); 
+      AddHit(fIshunt, primary, absid, xyzte);  
       
       if (cpvDigit->GetQpad() > 0.02) {
        xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
@@ -377,86 +280,103 @@ void AliPHOSv1::StepManager(void)
        qsum  += cpvDigit->GetQpad();
       }
     }
-    if (cpvDigits) {
-      cpvDigits->Delete();
-      delete cpvDigits;
-      cpvDigits=0;
-    }
+    fCPVDigits.Clear();
   }
 
  
-  
-  if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) { //  We are inside a PBWO crystal
+  static Int_t idPXTL = TVirtualMC::GetMC()->VolId("PXTL");  
+  if(TVirtualMC::GetMC()->CurrentVolID(copy) == idPXTL ) { //  We are inside a PBWO crystal
 
-    gMC->TrackPosition(pos) ;
+    TVirtualMC::GetMC()->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(); 
+    Float_t lostenergy = TVirtualMC::GetMC()->Edep(); 
     
     //Put in the TreeK particle entering PHOS and all its parents
-    if ( gMC->IsTrackEntering() ){
+    if ( TVirtualMC::GetMC()->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) ; 
+      TVirtualMC::GetMC() -> 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() ;
+       TVirtualMC::GetMC() -> Gmtod (vert, vertd, 1);    // transform coordinate from master to daughter system
+       if(vertd[1]<-GetGeometry()->GetCrystalSize(1)/2.-0.1){ //Particle is created in foront of PHOS 
+                                                              //0.1 to get rid of numerical errors 
          part->SetBit(kKeepBit);
-         parent = part->GetFirstMother() ; 
+         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() ;     
+      xyzte[3] = TVirtualMC::GetMC()->TrackTime() ;     
       
-      gMC->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ;
+      TVirtualMC::GetMC()->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ;
       
       Int_t strip ;
-      gMC->CurrentVolOffID(3, strip);
+      TVirtualMC::GetMC()->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 ;
-      
+      TVirtualMC::GetMC()->CurrentVolOffID(2, cell);
+
+      //Old formula for row is wrong. For example, I have strip 56 (28 for 2 x 8), row must be 1.
+      //But row == 1 + 56 - 56 % 56 == 57 (row == 1 + 28 - 28 % 28 == 29)
+      //Int_t row = 1 + GetGeometry()->GetEMCAGeometry()->GetNStripZ() - strip % (GetGeometry()->GetEMCAGeometry()->GetNStripZ()) ;
+      Int_t row = GetGeometry()->GetEMCAGeometry()->GetNStripZ() - (strip - 1) % (GetGeometry()->GetEMCAGeometry()->GetNStripZ()) ;
+      Int_t col = (Int_t) TMath::Ceil((Double_t) strip/(GetGeometry()->GetEMCAGeometry()->GetNStripZ())) -1 ;
+
+      // Absid for 8x2-strips. Looks nice :) 
       absid = (moduleNumber-1)*GetGeometry()->GetNCristalsInModule() + 
-       row + (col*GetGeometry()->GetEMCAGeometry()->GetNCellsInStrip() + cell-1)*GetGeometry()->GetNZ() ;
-      
-      gMC->Gmtod(global, local, 1) ;
+                   row * 2 + (col*GetGeometry()->GetEMCAGeometry()->GetNCellsXInStrip() + (cell - 1) / 2)*GetGeometry()->GetNZ() - (cell & 1 ? 1 : 0);
+
       
       //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 ))
-                                           ) ;
+      //There is no dependence of reponce on distance from energy deposition to APD
+      Float_t lightYield = gRandom->Poisson(AliPHOSSimParam::GetInstance()->GetLightFactor() * lostenergy) ;
 
       //Calculates de energy deposited in the crystal  
-      xyzte[4] = fAPDFactor * lightYield  ;
+      xyzte[4] = AliPHOSSimParam::GetInstance()->GetAPDFactor() * lightYield ;
       
-      // add current hit to the hit list
-      //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  
+      Int_t primary ;
+      if(fIshunt == 2){
+       primary = gAlice->GetMCApp()->GetCurrentTrackNumber() ;
+       TParticle * part = gAlice->GetMCApp()->Particle(primary) ;
+       while ( !part->TestBit(kKeepBit) ) {
+         primary = part->GetFirstMother() ;
+         if(primary == -1){        
+           primary  =  gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() ); 
+           break ; //there is a possibility that particle passed e.g. thermal isulator and hits a side 
+         //surface of the crystal. In this case it may have no primary at all. 
+         //We can not easily separate this case from the case when this is part of the shower, 
+         //developed in the neighboring crystal.
+         }
+         part = gAlice->GetMCApp()->Particle(primary) ;
+       }
+      }
+      else{
+       primary  =  gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() ); 
+      }
       
+      // add current hit to the hit list
+      // Info("StepManager","%d %d", primary, tracknumber) ; 
+      AddHit(fIshunt, primary, absid, xyzte);
+        
     } // there is deposited energy
   } // we are inside a PHOS Xtal
   
 }
 
 //____________________________________________________________________________
-void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, Int_t moduleNumber, TClonesArray *cpvDigits)
+void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, TClonesArray *cpvDigits)
 {
   // ------------------------------------------------------------------------
   // Digitize one CPV hit:
@@ -489,8 +409,6 @@ 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;
-
   Float_t dZY   = pZ/pNorm * GetGeometry()->GetCPVGasThickness();
   Float_t dXY   = pX/pNorm * GetGeometry()->GetCPVGasThickness();
   gRandom->Rannor(rnor1,rnor2);