]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSv1.cxx
Smaller changes
[u/mrichter/AliRoot.git] / PHOS / AliPHOSv1.cxx
index d0db7854ebdbb6cf23402ab23aba8811898d66a8..6155bacd7490d72f4353fc9c58aba01a452d9883 100644 (file)
 #include "TNode.h"
 #include "TRandom.h"
 #include "TTree.h"
-
+#include "TParticle.h"
 
 // --- Standard library ---
 
-#include <stdio.h>
 #include <string.h>
 #include <stdlib.h>
 #include <strstream.h>
@@ -75,7 +74,9 @@ AliPHOSv0()
   fLightYieldAttenuation  = 0. ;  
   fRecalibrationFactor    = 0. ;    
   fElectronsPerGeV        = 0. ;
-  fAPDGain                = 0. ;    
+  fAPDGain                = 0. ;  
+  fLightFactor            = 0. ; 
+  fAPDFactor              = 0. ; 
 
 }
 
@@ -98,10 +99,11 @@ 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) ; 
 
   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
@@ -116,18 +118,23 @@ AliPHOSv1::AliPHOSv1(const char *name, const char *title):
   // 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. ;
+  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  = 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++ ) {
@@ -162,15 +169,26 @@ AliPHOSv1::~AliPHOSv1()
     delete fHits ;
     fHits = 0 ; 
   }
-  if (fTreeQA) 
-    delete fTreeQA ; 
+  
+  delete fTreeQA ; 
+
+  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)
 {
   // 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 soem taime gate
+  // 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 ;
@@ -263,7 +281,7 @@ void AliPHOSv1::StepManager(void)
 
   Int_t          relid[4] ;           // (box, layer, row, column) indices
   Int_t          absid    ;           // absolute cell ID number
-  Float_t        xyze[5]={-1000,-1000,-1000,0,0}  ; // position wrt MRS, time 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     ;
 
@@ -276,31 +294,31 @@ void AliPHOSv1::StepManager(void)
   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
+    Float_t        xyd[3]={0,0,0}   ;   //local position of the entering
     xyd[0]  = xyzd[0];
-    xyd[1]  =-xyzd[1];
-    xyd[2]  =-xyzd[2];
-  
+    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    
@@ -350,10 +368,10 @@ void AliPHOSv1::StepManager(void)
       
       // add current digit to the temporary hit list
 
-      xyze[3] = gMC->TrackTime() ;
-      xyze[4] = cpvDigit->GetQpad() ;                          // amplitude in a pad
+      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, xyze);
+      AddHit(fIshunt, primary, tracknumber, absid, xyzte);
       
       if (cpvDigit->GetQpad() > 0.02) {
        xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
@@ -361,7 +379,11 @@ void AliPHOSv1::StepManager(void)
        qsum  += cpvDigit->GetQpad();
       }
     }
-    delete cpvDigits;
+    if (cpvDigits) {
+      cpvDigits->Delete();
+      delete cpvDigits;
+      cpvDigits=0;
+    }
   }
 
  
@@ -369,55 +391,70 @@ void AliPHOSv1::StepManager(void)
   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] ;
+    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() ;     
       
-      xyze[3] = gMC->TrackTime() ;     
-
       gMC->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ;
       
-      // fill the relevant QA Checkables
-      fQATotEner->Update( xyze[4] ) ;                                             // total energy in PHOS
-      (static_cast<AliPHOSQAFloatCheckable*>((*fQATotEnerB)[moduleNumber-1]))->Update( xyze[4] ) ; // energy in this block  
-
       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 photns produced in the
+      //Calculates the light yield, the number of photons produced in the
       //crystal 
-      Float_t lightYield = gRandom->Poisson(fLightYieldMean * lostenergy *
-                                           fIntrinsicPINEfficiency * 
+      Float_t lightYield = gRandom->Poisson(fLightFactor * lostenergy *
                                            exp(-fLightYieldAttenuation *
                                                (local[1]+GetGeometry()->GetCrystalSize(1)/2.0 ))
                                            ) ;
+
       //Calculates de energy deposited in the crystal  
-      xyze[4] = (fRecalibrationFactor/100.) * fAPDGain * lightYield  ;
+      xyzte[4] = fAPDFactor * lightYield  ;
       
       // 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
-
+  
 }
 
 //____________________________________________________________________________