Moving to the new VMC naming convention
[u/mrichter/AliRoot.git] / PHOS / AliPHOSv3.cxx
index 87a718c..d8050c9 100644 (file)
@@ -24,7 +24,7 @@
 //  PIN due to MIPS particle and electronic noise.
 // This is done in the StepManager 
 //                  
-//*-- Author:  Odd Harald Oddland & Gines Martinez (SUBATECH)
+//*-- Author:  Odd Harald Odland & Gines Martinez (SUBATECH)
 
 
 // --- ROOT system ---
@@ -32,7 +32,6 @@
 
 // --- Standard library ---
 
-#include <stdio.h>
 #include <string.h>
 #include <stdlib.h>
 #include <strstream.h>
 
 #include "AliPHOSv3.h"
 #include "AliPHOSHit.h"
-#include "AliPHOSDigit.h"
+#include "AliPHOSCPVDigit.h"
 #include "AliRun.h"
 #include "AliConst.h"
 #include "AliMC.h"
+#include "AliPHOSGeometry.h"
 
 ClassImp(AliPHOSv3)
 
 //____________________________________________________________________________
-  AliPHOSv3::AliPHOSv3(const char *name, const char *title):
-AliPHOSv1(name,title)
+  AliPHOSv3::AliPHOSv3(void) : AliPHOSv1() 
 {
-  // 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 from Odd Harald Odland work
-  // k_0 is 0.0045 from Valery Antonenko 
-
-
-  fLightYieldMean = 700000. ;
-  fIntrinsicPINEfficiency = 0.1875 ;
-  fLightYieldAttenuation = 0.0045 ;
-  fRecalibrationFactor = 6.2 / fLightYieldMean ;
-  fElectronsPerGeV = 2.77e+8 ; 
+  // default ctor: initialize daa members
+
+  fLightYieldMean         = 0. ;         
+  fIntrinsicPINEfficiency = 0. ; 
+  fLightYieldAttenuation  = 0. ;  
+  fRecalibrationFactor    = 0. ;    
+  fElectronsPerGeV        = 0. ;
+  fAPDGain                = 0. ;    
+  
 }
 
 //____________________________________________________________________________
-AliPHOSv3::AliPHOSv3(AliPHOSReconstructioner * Reconstructioner, const char *name, const char *title):
-  AliPHOSv1(Reconstructioner,name,title)
+  AliPHOSv3::AliPHOSv3(const char *name, const char *title):
+AliPHOSv1(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 from Odd Harald Odland work
+  // 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 
-
-  fLightYieldMean = 700000.;
-  fIntrinsicPINEfficiency = 0.1875 ;
-  fLightYieldAttenuation = 0.0045 ;
-  fRecalibrationFactor = 6.2 / fLightYieldMean ;
+  // 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. ;
 }
+
+// //____________________________________________________________________________
+// 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)
@@ -99,54 +116,29 @@ void AliPHOSv3::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;
+//      Info("StepManager", "Track enters the volume %d", gMC->CurrentVolName()) ;
 //    if (gMC->IsTrackExiting())
-//      cout << "Track leaves the volume " << gMC->CurrentVolName() << endl;
+//      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
+  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[2]   ;      //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[2];
-      
-      // 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
+  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("GCEL") ) // We are inside a gas cell 
+    if( gMC->CurrentVolID(copy) == gMC->VolId("PPCE") ) // We are inside a gas cell 
     {
       gMC->TrackPosition(pos) ;
       xyze[0] = pos[0] ;
@@ -154,23 +146,23 @@ void AliPHOSv3::StepManager(void)
       xyze[2] = pos[2] ;
       xyze[3] = gMC->Edep() ; 
 
-      if ( (xyze[3] != 0) || entered ) { // there is deposited energy or new particle entering  PPSD
+      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();
+         relid[0] += GetGeometry()->GetNModules() - GetGeometry()->GetNPPSDModules();
        }
                gMC->CurrentVolOffID(3, relid[1]) ;  // get the Micromegas Module number 
-      // 1-> fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() upper
-      //   > fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() lower
+      // 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
 
-               fGeom->RelToAbsNumbering(relid, absid) ; 
+               GetGeometry()->RelToAbsNumbering(relid, absid) ; 
 
        // add current hit to the hit list      
-         AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid, pmom, xyd);
+         AddHit(fIshunt, primary, tracknumber, absid, xyze);
 
 
       } // there is deposited energy 
@@ -181,9 +173,31 @@ void AliPHOSv3::StepManager(void)
 
     // Yuri Kharlov, 28 September 2000
 
-    if( gMC->CurrentVolID(copy) == gMC->VolId("CPVQ") &&
-       entered &&
+    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:
 
@@ -206,11 +220,11 @@ void AliPHOSv3::StepManager(void)
 
       ndigits = cpvDigits->GetEntriesFast();
       for (idigit=0; idigit<ndigits-1; idigit++) {
-       AliPHOSCPVDigit  *cpvDigit1 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(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 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(jdigit);
+         AliPHOSCPVDigit  *cpvDigit2 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(jdigit));
          Float_t x2 = cpvDigit2->GetXpad() ;
          Float_t z2 = cpvDigit2->GetYpad() ;
          if (x1==x2 && z1==z2) {
@@ -226,14 +240,14 @@ void AliPHOSv3::StepManager(void)
 
       ndigits = cpvDigits->GetEntriesFast();
       for (idigit=0; idigit<ndigits; idigit++) {
-       AliPHOSCPVDigit  *cpvDigit = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(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
-       fGeom->RelToAbsNumbering(relid, absid) ; 
+       GetGeometry()->RelToAbsNumbering(relid, absid) ; 
 
        // add current digit to the temporary hit list
        xyze[0] = 0. ;
@@ -241,7 +255,7 @@ void AliPHOSv3::StepManager(void)
        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);
+       AddHit(fIshunt, primary, tracknumber, absid, xyze);
 
        if (cpvDigit->GetQpad() > 0.02) {
          xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
@@ -254,63 +268,111 @@ void AliPHOSv3::StepManager(void)
   } // end of IHEP configuration
   
 
-  if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) { //  We are inside a PBWO crystal
+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) || entered ) {  // Track is inside the crystal and deposits some energy or just entered 
+  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();      
+       relid[0] += GetGeometry()->GetNModules() - GetGeometry()->GetNPPSDModules();      
 
       relid[1] = 0   ;                    // means PBW04
-      gMC->CurrentVolOffID(4, relid[2]) ; // get the row number inside the module
-      gMC->CurrentVolOffID(3, relid[3]) ; // get the cell number inside the module
+   gMC->CurrentVolOffID(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
-      fGeom->RelToAbsNumbering(relid, absid) ; 
+
+   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
-       AddHit(fIshunt, primary,tracknumber, absid, xyze, trackpid,pmom, xyd);
 
+   if (xyze[3] != 0.) AddHit(fIshunt, primary,tracknumber, absid, xyze);
+  
 
     } // there is deposited energy
   } // we are inside a PHOS Xtal
 
-  if(gMC->CurrentVolID(copy) == gMC->VolId("PPIN") ) // We are inside de PIN diode 
-    {
-      gMC->TrackPosition(pos) ;
-      xyze[0] = pos[0] ;
-      xyze[1] = pos[1] ;
-      xyze[2] = pos[2] ;
-      Float_t lostenergy = gMC->Edep() ;
-      xyze[3] = gMC->Edep() ;
+//  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
+//       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
        
-       Int_t absid ; 
-       fGeom->RelToAbsNumbering(relid,absid) ;
+//     // get the absolute Id number
        
-       // calculating number of electrons in the PIN diode asociated to this hit
-         Float_t nElectrons = lostenergy * fElectronsPerGeV ;
-         xyze[3] = nElectrons * fRecalibrationFactor ;
+//     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, trackpid,pmom,xyd);
-         //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
+//       // 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
 }