Moving to the new VMC naming convention
[u/mrichter/AliRoot.git] / PHOS / AliPHOSv3.cxx
index 6c4ff74..d8050c9 100644 (file)
-/////////////////////////////////////////////////////////
-//  Manager and hits classes for set:PHOS version 3    //
-/////////////////////////////////////////////////////////
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//_________________________________________________________________________
+// Implementation version v0 of PHOS Manager class 
+// Layout EMC + PPSD has name GPS2  
+// The main goal of this version of AliPHOS is to calculte the 
+//  induced charged in the PIN diode, taking into account light
+//  tracking in the PbWO4 crystal, induced signal in the 
+//  PIN due to MIPS particle and electronic noise.
+// This is done in the StepManager 
+//                  
+//*-- Author:  Odd Harald Odland & Gines Martinez (SUBATECH)
+
+
 // --- ROOT system ---
-#include "TH1.h"
 #include "TRandom.h"
-#include "TFile.h"
-#include "TTree.h"
-#include "TBRIK.h"
-#include "TNode.h"
 
-// --- galice header files ---
+// --- Standard library ---
+
+#include <string.h>
+#include <stdlib.h>
+#include <strstream.h>
+
+// --- AliRoot header files ---
+
 #include "AliPHOSv3.h"
+#include "AliPHOSHit.h"
+#include "AliPHOSCPVDigit.h"
 #include "AliRun.h"
-#include "TGeant3.h" 
+#include "AliConst.h"
+#include "AliMC.h"
+#include "AliPHOSGeometry.h"
 
 ClassImp(AliPHOSv3)
 
-//______________________________________________________________________________
-
-
-AliPHOSv3::AliPHOSv3()
+//____________________________________________________________________________
+  AliPHOSv3::AliPHOSv3(void) : AliPHOSv1() 
 {
+  // default ctor: initialize daa members
+
+  fLightYieldMean         = 0. ;         
+  fIntrinsicPINEfficiency = 0. ; 
+  fLightYieldAttenuation  = 0. ;  
+  fRecalibrationFactor    = 0. ;    
+  fElectronsPerGeV        = 0. ;
+  fAPDGain                = 0. ;    
+  
 }
-//______________________________________________________________________________
 
-AliPHOSv3::AliPHOSv3(const char *name, const char *title)
-          : AliPHOS(name, title)
+//____________________________________________________________________________
+  AliPHOSv3::AliPHOSv3(const char *name, const char *title):
+AliPHOSv1(name,title)
 {
+  // ctor 
+
+  // 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. ;
 }
-//___________________________________________
-void AliPHOSv3::CreateGeometry()
-{
 
-  cout << "AliPHOSv3::CreateGeometry() PHOS creation\n";
-
-  AliMC* pMC = AliMC::GetMC();
-
-  AliPHOS *PHOS_tmp = (AliPHOS*)gAlice->GetModule("PHOS");
-  if( NULL==PHOS_tmp )
-  {
-    printf("There isn't PHOS detector!\n");
-    return;
-  }
-//  AliPHOS &PHOS = *PHOS_tmp;
-
-  //////////////////////////////////////////////////////////////////////////////
-
-  Int_t                 rotation_matrix_number=0;
-  Float_t               par[11],
-                        x,y,z;
-
-  const float           cell_length             = GetCrystalLength()+GetAirThickness()+GetWrapThickness()+GetPIN_Length(),
-                        cell_side_size          = GetCrystalSideSize()+2*GetAirThickness()+2*GetWrapThickness(),
-//                        cell_angle              = 180/kPI * 2 * atan(cell_side_size/2 / GetRadius()),        // radians
-                        cradle_thikness         = cell_length + GetCPV_Thickness() + GetCPV_PHOS_Distance(),
-                        distance_to_CPV         = GetRadius() - GetCPV_Thickness() - GetCPV_PHOS_Distance();
-
-  //////////////////////////////////////////////////////////////////////////////
-  // CELL volume and subvolumes creation
-  //////////////////////////////////////////////////////////////////////////////
-
-  par[0] = GetCrystalSideSize()/2 + GetWrapThickness();
-  par[1] = GetCrystalSideSize()/2 + GetWrapThickness();
-  par[2] = GetCrystalLength()  /2 + GetWrapThickness()/2;
-  pMC->Gsvolu("WRAP","BOX ",GetPHOS_IDTMED_Tyvek(),par,3);
-
-  par[0] = GetCrystalSideSize()/2;
-  par[1] = GetCrystalSideSize()/2;
-  par[2] = GetCrystalLength()/2;
-  pMC->Gsvolu("CRST","BOX ",GetPHOS_IDTMED_PbWO4(),par,3);
-
-  // PIN
-  par[0] = GetPIN_SideSize()/2;
-  par[1] = GetPIN_SideSize()/2;
-  par[2] = GetPIN_Length()/2;
-  pMC->Gsvolu("PIN ","BOX ",GetPHOS_IDTMED_PIN(),par,3);
-
-  //////////////////////////////////////////////////////////////////////////////
-  // CRADLE,CPV creation.
-  //////////////////////////////////////////////////////////////////////////////
-
-  par[0] = cell_side_size/2 * GetNz();
-  par[1] = cell_side_size/2 * GetNphi();
-  par[2] = cradle_thikness/2;
-  pMC->Gsvolu("PHOS","BOX ",GetPHOS_IDTMED_AIR(),par,3);
-
-//par[0] : the same as above
-//par[1] : the same as above
-  par[2] = GetCPV_Thickness()/2;
-  pMC->Gsvolu("CPV ","BOX ",GetPHOS_IDTMED_CPV(),par,3);
-
-  x = 0;
-  y = 0;
-  z = (cell_length+GetCPV_PHOS_Distance())/2;
-  pMC->Gspos("CPV ",1,"PHOS",x,y,z,0,"ONLY");
-
-  par[0] = cell_side_size/2 * GetNz();
-  par[1] = cell_side_size/2 * GetNphi();
-  par[2] = cell_length/2;
-  pMC->Gsvolu("CRS0","BOX ",GetPHOS_IDTMED_AIR(),par,3);
-
-  x = 0;
-  y = 0;
-  z = -(cradle_thikness-cell_length)/2;
-  pMC->Gspos("CRS0",1,"PHOS",x,y,z,0,"ONLY");
-
-  pMC->Gsdvn("CRS1","CRS0",GetNphi(),2);
-  pMC->Gsdvn("CELL","CRS1",GetNz()  ,1);
-
-  //////////////////////////////////////////////////////////////////////////////
-  // CELL creation
-  //////////////////////////////////////////////////////////////////////////////
-
-  x = 0;
-  y = 0;
-  z = -GetWrapThickness()/2;
-  pMC->Gspos("CRST",1,"WRAP",x,y,z,0,"ONLY");
-
-  x = 0;
-  y = 0;
-  z = GetPIN_Length()/2;
-  pMC->Gspos("WRAP",1,"CELL",x,y,z,0,"ONLY");
-
-  x = 0;
-  y = 0;
-  z = -GetCrystalLength()/2-GetWrapThickness()/2;
-  pMC->Gspos("PIN ",1,"CELL",x,y,z,0,"ONLY");
-
-  //////////////////////////////////////////////////////////////////////////////
-  // CELL has been created.
-  //////////////////////////////////////////////////////////////////////////////
-
-//   int n=0;
-//   z = -(GetCPV_Thickness()+GetCPV_PHOS_Distance())/2;
-//
-//   for( int iy=0; iy<GetNphi(); iy++ )
-//   {
-//     y = (iy-(GetNphi()-1)/2.)*cell_side_size;
-//     for( int ix=0; ix<GetNz(); ix++ )
-//     {
-//       x = (ix-(GetNz()-1)/2.)*cell_side_size;
-//       pMC->Gspos("CELL",++n,"PHOS",x,y,z,0,"ONLY");
-//     }
-//   }
-
-  //////////////////////////////////////////////////////////////////////////////
-  // End of CRADLE creation.
-  //////////////////////////////////////////////////////////////////////////////
-
-
-  //////////////////////////////////////////////////////////////////////////////
-  // PHOS creation
-  //////////////////////////////////////////////////////////////////////////////
-
-  for( int i=0; i<GetCradlesAmount(); i++ )
-  {
-    float c                = distance_to_CPV,           // Distance to CPV
-          l                = cell_side_size*GetNphi()/2,      // Cradle half size around beam (for rect. geom.)
-          cradle_angle     = 360/kPI*atan(l/c),
-          cradle_angle_pos = -90+(i-(GetCradlesAmount()-1)/2.) * (cradle_angle+GetAngleBetweenCradles());
-    // Cradles are numerated in clock reversed order. (general way of angle increment)
-
-    float   r       = GetRadius() + cradle_thikness/2;
-    x = r*cos(cradle_angle_pos*kPI/180);
-    y = r*sin(cradle_angle_pos*kPI/180);
-    z = 0;
-    AliMatrix(rotation_matrix_number, 0,0 , 90,90+cradle_angle_pos , 90,180+cradle_angle_pos);
-    pMC->Gspos("PHOS",i+1,"ALIC",x,y,z,rotation_matrix_number,"ONLY");
-
-    GetCradleAngle(i) = cradle_angle_pos;
-//
-//    int n = PHOS.fCradles->GetEntries();
-//    PHOS.fCradles->Add(new AliPHOSCradle( 1,            // geometry.
-//                                          GetCrystalSideSize    (),
-//                                          GetCrystalLength      (),
-//                                          GetWrapThickness      (),
-//                                          GetAirThickness       (),
-//                                          GetPIN_SideSize       (),
-//                                          GetPIN_Length         (),
-//                                          GetRadius             (),
-//                                          GetCPV_Thickness      (),
-//                                          GetCPV_PHOS_Distance  (),
-//                                          GetNz                 (),
-//                                          GetNphi               (),
-//                                          cradle_angle_pos      ));
-//
-//    if( n+1 != PHOS.fCradles->GetEntries() ||
-//        NULL == PHOS.fCradles->At(n) )
-//    {
-//      cout << "  Can not create or add AliPHOSCradle.\n";
-//      exit(1);
-//    }
-  }
-  AddPHOSCradles();
-
-  //////////////////////////////////////////////////////////////////////////////
-  // All is done.
-  // Print some information.
-  //////////////////////////////////////////////////////////////////////////////
-}
-
-void AliPHOSv3::StepManager()
+// //____________________________________________________________________________
+// AliPHOSv3::AliPHOSv3(AliPHOSReconstructioner * Reconstructioner, const char *name, const char *title):
+//   AliPHOSv1(Reconstructioner,name,title)
+// {
+//   // ctor 
+
+//   // Number of electrons created in the PIN due to light collected in the PbWo4 crystal is calculated using 
+//   // following formula
+//   // NumberOfElectrons = EnergyLost * LightYield * PINEfficiency * 
+//   //                     exp (-LightYieldAttenuation * DistanceToPINdiodeFromTheHit) *
+//   //                     RecalibrationFactor ;
+//   // LightYield is obtained as a Poissonian distribution with a mean at 700000 photons per GeV fromValery Antonenko
+//   // PINEfficiency is 0.1875 
+//   // k_0 is 0.0045 from Valery Antonenko 
+
+//   fLightYieldMean = 700000.;
+//   fIntrinsicPINEfficiency = 0.1875 ;
+//   fLightYieldAttenuation = 0.0045 ;
+//   fRecalibrationFactor = 6.2 / fLightYieldMean ;
+//   fRecalibrationFactor = 5.67 /fLightYieldMean ;//25.04.2001 OHO
+//   fElectronsPerGeV = 2.77e+8 ;
+// }
+//____________________________________________________________________________
+
+void AliPHOSv3::StepManager(void)
 {
-  static Bool_t inwold=0;   // Status of previous ctrak->inwvol
-  AliMC *MC = AliMC::GetMC();
-  Int_t copy;
-
-//   if( MC->TrackEntering() ) {
-//     Int_t Volume_ID = MC->CurrentVol(Volume_name, copy);
-//     cout << "AliPHOSv3::StepManager() entered to PHOS to the volume " << Volume_name << "!\n";
-//   }
-
-  int cradle_number, cell_Z, cell_Phi;  // Variables that describe cell position.
-
-  if( MC->GetMedium()==GetPHOS_IDTMED_PIN() && MC->TrackEntering() && MC->TrackCharge()!=0 )
-  {
-    // GEANT particle just have entered into PIN diode.
-
-    AliPHOS &PHOS = *(AliPHOS*)gAlice->GetModule("PHOS");
-
-    MC->CurrentVolOff(4,0,copy);
-    cradle_number  = copy-1;
-    MC->CurrentVolOff(1,0,copy);
-    cell_Z         = copy-1;
-    MC->CurrentVolOff(2,0,copy);
-    cell_Phi       = copy-1;
-
-    TH2S &h = PHOS.GetCradle(cradle_number).fChargedTracksInPIN;
-    h.AddBinContent(h.GetBin(cell_Z,cell_Phi));
-
-//     cout << "AliPHOSv3::StepManager() entered to PHOS pin diode\n";
-//     cout << "   cradle_nimber = " << cradle_number << endl;
-//     cout << "   cell_z        = " << cell_Z << endl;
-//     cout << "   cell_Phi      = " << cell_Phi << endl;
-  }
-
-  //////////////////////////////////////////////////////////////////////////////
-
-  if( MC->GetMedium() == GetPHOS_IDTMED_PbWO4() )
-  {
-    // GEANT particle into crystal.
-
-    AliPHOS &PHOS = *(AliPHOS*)gAlice->GetModule("PHOS");
-
-    MC->CurrentVolOff(5,0,copy);
-    cradle_number  = copy-1;
-    MC->CurrentVolOff(2,0,copy);
-    cell_Z         = copy-1;
-    MC->CurrentVolOff(3,0,copy);
-    cell_Phi       = copy-1;
-
-    TH2F &h = PHOS.GetCradle(cradle_number).fCellEnergy;
-    h.AddBinContent(h.GetBin(cell_Z,cell_Phi),MC->Edep());
-  }
-
-  //////////////////////////////////////////////////////////////////////////////
-
-  if( MC->GetMedium()==GetPHOS_IDTMED_CPV() && MC->TrackEntering() )
-  {
-    // GEANT particle just have entered into CPV detector.
-
-    AliPHOS &PHOS = *(AliPHOS*)gAlice->GetModule("PHOS");
-
-    MC->CurrentVolOff(1,0,cradle_number);
-    cradle_number--;
-
-    // Save CPV x,y hits position of charged particles.
-
-    AliPHOSCradle  &cradle = PHOS.GetCradle(cradle_number);
-
-    Float_t   xyz[3];
-    MC->TrackPosition(xyz);
-    TVector3          p(xyz[0],xyz[1],xyz[2]),v;
-
-    float x,y,l;
-    float R = cradle.GetRadius() - cradle.GetCPV_PHOS_Distance() - cradle.GetCPV_Thikness();
-    cradle.GetXY(p,v,R,x,y,l);
-    if( PHOS.fDebugLevel>0 )
-      if( l<0 )
-        printf("PHOS_STEP:  warning: negative distance to CPV!! %f\n", l);
-
-    // Store current particle in the list of Cradle particles.
-    Float_t  pmom[4];
-    MC->TrackMomentum(pmom);
-    Float_t Px    = pmom[0] * pmom[3],
-            Py    = pmom[1] * pmom[3],
-            Pz    = pmom[2] * pmom[3];
-    Float_t Getot = MC->Etot();
-    Int_t   Ipart = MC->TrackPid();
-
-    cradle.GetParticles().Add(new AliPHOSgamma(x,y,Getot,Px,Py,Pz,Ipart));
+  // Accumulates hits as long as the track stays in a single crystal or PPSD gas Cell
+
+//    if (gMC->IsTrackEntering())
+//      Info("StepManager", "Track enters the volume %d", gMC->CurrentVolName()) ;
+//    if (gMC->IsTrackExiting())
+//      Info("StepManager", "Track leaves the volume %d", gMC->CurrentVolName()) ;
+
+  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
+  Int_t          copy     ;
+  Float_t        fLightYield ;   // Light Yield per GeV
+
+  Int_t tracknumber =  gAlice->GetCurrentTrackNumber() ; 
+  Int_t primary     =  gAlice->GetPrimary( gAlice->GetCurrentTrackNumber() ); 
+  TString name      =  GetGeometry()->GetName() ; 
+  Float_t        lostenergy ;
+  Float_t        global[3] ;
+  Float_t        local[3] ;
+
+
+  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) {      
+
+      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];
+
+      
+      // 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
+      
+      Int_t moduleNumber;
+      gMC->CurrentVolOffID(3,moduleNumber);
+      moduleNumber--;
+
+
+      TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0);   // array of digits for current hit
+      CPVDigitize(pmom,xyd,moduleNumber,cpvDigits);
+      
+      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
+       
+       // 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();
+       }
+      }
+      delete cpvDigits;
+    }
+  } // end of IHEP configuration
+  
+
+if(gMC->CurrentVolID(copy)==gMC->VolId("PXTL")){// We are inside a PBWO4 crystal
+    gMC->TrackPosition(pos) ;
+    xyze[0] = pos[0] ;
+    xyze[1] = pos[1] ;
+    xyze[2] = pos[2] ;
+    global[0] = pos[0] ;
+    global[1] = pos[1] ;
+    global[2] = pos[2] ;
+    lostenergy = gMC->Edep(); 
+    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] += 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
+      
+      // get the absolute Id number
+
+   GetGeometry()->RelToAbsNumbering(relid, absid) ; 
+
+   gMC->Gmtod(global, local, 1) ;
+
+   //Calculates the light yield, the number of photns produced in the
+   //crystal 
+   fLightYield = gRandom->Poisson(
+                                 fLightYieldMean * lostenergy *
+                                 fIntrinsicPINEfficiency * 
+                                 exp(-fLightYieldAttenuation *
+                                     (local[1]+GetGeometry()->GetCrystalSize(1)/2.0 ))
+                                 ) ;
+   //Calculates de energy deposited in the crystal  
+   xyze[3] = (fRecalibrationFactor/100.) * fAPDGain * fLightYield  ;
+  
+    
+   // Info("StepManager", "xyze[3]: %f", xyze[3]) ;
+   // Info("StepManager", "lostenergy: %f", lostenergy) ; 
+
+
+        
+
+      // add current hit to the hit list
 
-    //    printf ("Cradle %i, x,y = %8.3f, %8.3f cm, E,Px,Py,Pz = %8.3f, %8.3f, %8.3f GeV, %8.3f, Ipart = %i\n",
-    //       cradle_number,x,y,Getot,Px,Py,Pz,Ipart);
+   if (xyze[3] != 0.) AddHit(fIshunt, primary,tracknumber, absid, xyze);
+  
 
-  }
+    } // there is deposited energy
+  } // we are inside a PHOS Xtal
 
-  inwold=MC->TrackEntering();         // Save current status of GEANT variable.
+//  if(gMC->CurrentVolID(copy) == gMC->VolId("PPIN"))//We are inside the PIN diode 
+//     {
+//       Info("StepManager", "Inside PIN";
+//       gMC->TrackPosition(pos) ;
+//       global[0] = pos[0] ;
+//       global[1] = pos[1] ;
+//       global[2] = pos[2] ;
+//       xyze[0] = pos[0] ;
+//       xyze[1] = pos[1] ;
+//       xyze[2] = pos[2] ;
+//       lostenergy = gMC->Edep() ;
+//       xyze[3] = gMC->Edep() ;
+      
+//       if ( xyze[3] != 0 ) {
+//     gMC->CurrentVolOffID(11, relid[0]) ; // get the PHOS module number ;
+//     relid[1] = 0   ;                    // means PW04 and PIN
+//     gMC->CurrentVolOffID(5, relid[2]) ; // get the row number inside the module
+//     gMC->CurrentVolOffID(4, relid[3]) ; // get the cell number inside the module
+       
+//     // get the absolute Id number
+       
+//     Int_t absid ; 
+//     GetGeometry()->RelToAbsNumbering(relid,absid) ;
+//     gMC->Gmtod(global, local, 1) ;
+
+// // calculating number of electrons in the PIN diode asociated to this hit
+
+//       //nElectrons = lostenergy * fElectronsPerGeV ;
+//        //  xyze[3] = nElectrons * fRecalibrationFactor ;
+//           apdgain = gRandom->Poisson(300.) ;
+//           // apdgain = 300.;
+//    //if(local[1]<-0.0045) xyze[3] = apdgain * nElectrons * fRecalibrationFactor/10000.;
+
+//     if(local[1]<-0.0045) xyze[3] = apdgain * lostenergy * fElectronsPerGeV* (fRecalibrationFactor/100.);
+
+//    //if((local[1]>-0.0045)&&(gMC->TrackPid()==-11)) xyze[3] = apdgain * nElectrons * fRecalibrationFactor/10000.;
+
+//      if((local[1]>-0.0045)&&(gMC->TrackPid()==-11)) xyze[3] = apdgain * lostenergy * fElectronsPerGeV * (fRecalibrationFactor/100.);
+
+//    //if(local[1]>-0.0045) xyze[3] = nElectrons * fRecalibrationFactor/10000.;
+
+//      if(local[1]>-0.0045) xyze[3] = lostenergy * fElectronsPerGeV * (fRecalibrationFactor/100.);
+         
+//       // add current hit to the hit list
+
+//       AddHit(fIshunt, primary, tracknumber, absid, xyze);
+
+//       //printf("PIN volume is  %d, %d, %d, %d \n",relid[0],relid[1],relid[2],relid[3]);
+//       printf("Lost energy in the PIN is %f \n",lostenergy) ;
+//       } // there is deposited energy
+//     } // we are inside a PHOS XtalPHOS PIN diode
 }