]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSEmcRecPoint.cxx
more debug output and bugfix for sub event number
[u/mrichter/AliRoot.git] / PHOS / AliPHOSEmcRecPoint.cxx
index 9f2c2735edb51be9be339e1f9aa0e6546575ce21..87467a8184b2ad5417aa704b0b2247032ada881a 100644 (file)
 
 
 // --- ROOT system ---
-#include "TPad.h"
 #include "TH2.h"
 #include "TMath.h" 
 #include "TCanvas.h" 
+#include "TGraph.h"
 
 // --- Standard library ---
 
 // --- AliRoot header files ---
-
- #include "AliGenerator.h"
+#include "AliPHOSLoader.h"
+#include "AliGenerator.h"
 #include "AliPHOSGeometry.h"
+#include "AliPHOSDigit.h"
 #include "AliPHOSEmcRecPoint.h"
-#include "AliRun.h"
-#include "AliPHOSGetter.h"
-
 ClassImp(AliPHOSEmcRecPoint)
 
 //____________________________________________________________________________
@@ -52,6 +51,7 @@ AliPHOSEmcRecPoint::AliPHOSEmcRecPoint() : AliPHOSRecPoint()
   fNExMax     = 0 ;   //Not unfolded yet
   fTime = -1. ;
   fLocPos.SetX(1000000.)  ;      //Local position should be evaluated
+  fDebug=0;
    
 }
 
@@ -67,9 +67,26 @@ AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(const char * opt) : AliPHOSRecPoint(opt)
   fEnergyList = 0 ;
   fTime = -1. ;
   fLocPos.SetX(1000000.)  ;      //Local position should be evaluated
+  fDebug=0;
   
 }
 
+//____________________________________________________________________________
+AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(const AliPHOSEmcRecPoint & rp) : AliPHOSRecPoint(rp)
+{
+  // cpy ctor
+
+  fMulDigit   = rp.fMulDigit ;  
+  fAmp        = rp.fAmp ;   
+  fCoreEnergy = rp.fCoreEnergy ; 
+  fEnergyList = new Float_t[rp.fMulDigit] ;
+  Int_t index ; 
+  for(index = 0 ; index < fMulDigit ; index++) 
+    fEnergyList[index] = rp.fEnergyList[index] ; 
+  fNExMax     = rp.fNExMax ;  
+  fTime       = rp.fTime ;   
+}
+
 //____________________________________________________________________________
 AliPHOSEmcRecPoint::~AliPHOSEmcRecPoint()
 {
@@ -129,8 +146,7 @@ Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * d
   
   Bool_t aren = kFALSE ;
   
-  AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
-  AliPHOSGeometry * phosgeom =  (AliPHOSGeometry*)gime->PHOSGeometry();
+  AliPHOSGeometry * phosgeom =  AliPHOSLoader::GetPHOSGeometry();
 
   Int_t relid1[4] ; 
   phosgeom->AbsToRelNumbering(digit1->GetId(), relid1) ; 
@@ -190,7 +206,7 @@ Int_t AliPHOSEmcRecPoint::Compare(const TObject * obj) const
   return rv ; 
 }
 //______________________________________________________________________________
-void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py) const
+void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t, Int_t) const
 {
   
   // Execute action corresponding to one event
@@ -199,10 +215,8 @@ void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py) const
   //  If Left button is clicked on AliPHOSRecPoint, the digits are switched on    
   //  and switched off when the mouse button is released.
   
-    
-  AliPHOSGetter * gime =  AliPHOSGetter::GetInstance() ; 
-  if(!gime) return ;
-  AliPHOSGeometry * phosgeom =  (AliPHOSGeometry*)gime->PHOSGeometry();
+  
+  AliPHOSGeometry * phosgeom =  AliPHOSLoader::GetPHOSGeometry();
   
   static TGraph *  digitgraph = 0 ;
   
@@ -211,6 +225,22 @@ void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py) const
   TH2F * histo = 0 ;
   TCanvas * histocanvas ; 
 
+  
+  //try to get run loader from default event folder
+  AliRunLoader* rn = AliRunLoader::GetRunLoader(AliConfig::fgkDefaultEventFolderName);
+  if (rn == 0x0) 
+    {
+      Error("ExecuteEvent","Can not find Run Loader in Default Event Folder");
+      return;
+    }
+  AliPHOSLoader* gime = dynamic_cast<AliPHOSLoader*>(rn->GetLoader("PHOSLoader"));
+  if (gime == 0x0) 
+    {
+      Error("ExecuteEvent","Can not find PHOS Loader from Run Loader");
+      return;
+    }
+  
+  
   const TClonesArray * digits = gime->Digits() ;
   
   switch (event) {
@@ -312,9 +342,7 @@ void  AliPHOSEmcRecPoint::EvalDispersion(Float_t logWeight,TClonesArray * digits
 
   AliPHOSDigit * digit ;
  
-  AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
-  AliPHOSGeometry * phosgeom =  (AliPHOSGeometry*)gime->PHOSGeometry();
-  
+  AliPHOSGeometry * phosgeom =  AliPHOSLoader::GetPHOSGeometry();
 
  // Calculates the center of gravity in the local PHOS-module coordinates 
   
@@ -373,9 +401,8 @@ void AliPHOSEmcRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits
 
   AliPHOSDigit * digit ;
 
-  AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
-  AliPHOSGeometry * phosgeom =  (AliPHOSGeometry*)gime->PHOSGeometry();
-    
+  AliPHOSGeometry * phosgeom =  AliPHOSLoader::GetPHOSGeometry();
+
   Int_t iDigit;
 
 // Calculates the center of gravity in the local PHOS-module coordinates 
@@ -424,8 +451,7 @@ void  AliPHOSEmcRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
 
   AliPHOSDigit * digit ;
 
-  AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
-  AliPHOSGeometry * phosgeom =  (AliPHOSGeometry*)gime->PHOSGeometry();
+  AliPHOSGeometry * phosgeom =  AliPHOSLoader::GetPHOSGeometry();
 
   Int_t iDigit;
 
@@ -457,7 +483,7 @@ void  AliPHOSEmcRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
 //   //Apply correction due to non-perpendicular incidence
 //   Double_t CosX ;
 //   Double_t CosZ ;
-//   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
+//   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
 //   AliPHOSGeometry * phosgeom =  (AliPHOSGeometry*)gime->PHOSGeometry();
   //   Double_t DistanceToIP= (Double_t ) phosgeom->GetIPtoCrystalSurface() ;
   
@@ -480,6 +506,132 @@ void  AliPHOSEmcRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
     fLambda[1]= 0. ;
 }
 
+//____________________________________________________________________________
+void  AliPHOSEmcRecPoint::EvalMoments(Float_t logWeight,TClonesArray * digits)
+{
+  // Calculate the shower moments in the eigen reference system
+  // M2x, M2z, M3x, M4z
+  // Calculate the angle between the shower position vector and the eigen vector
+
+  Double_t wtot = 0. ;
+  Double_t x    = 0.;
+  Double_t z    = 0.;
+  Double_t dxx  = 0.;
+  Double_t dzz  = 0.;
+  Double_t dxz  = 0.;
+  Double_t lambda0=0, lambda1=0;
+
+  AliPHOSDigit * digit ;
+
+  AliPHOSGeometry * phosgeom =  AliPHOSLoader::GetPHOSGeometry();
+
+  Int_t iDigit;
+
+  // 1) Find covariance matrix elements:
+  //    || dxx dxz ||
+  //    || dxz dzz ||
+
+  Float_t xi ;
+  Float_t zi ;
+  Int_t relid[4] ;
+  Double_t w;
+  for(iDigit=0; iDigit<fMulDigit; iDigit++) {
+    digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit])  ;
+    phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+    phosgeom->RelPosInModule(relid, xi, zi);
+    w     = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
+    x    += w * xi ;
+    z    += w * zi ; 
+    dxx  += w * xi * xi ;
+    dzz  += w * zi * zi ;
+    dxz  += w * xi * zi ; 
+    wtot += w ;
+  }
+  x   /= wtot ;
+  z   /= wtot ;
+  dxx /= wtot ;
+  dzz /= wtot ;
+  dxz /= wtot ;
+  dxx -= x * x ;
+  dzz -= z * z ;
+  dxz -= x * z ;
+
+  // 2) Find covariance matrix eigen values lambda0 and lambda1
+
+  lambda0 =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
+  lambda1 =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
+
+  // 3) Find covariance matrix eigen vectors e0 and e1
+
+  TVector2 e0,e1;
+  if (dxz != 0)
+    e0.Set(1.,(lambda0-dxx)/dxz);
+  else
+    e0.Set(0.,1.);
+
+  e0 = e0.Unit();
+  e1.Set(-e0.Y(),e0.X());
+
+  // 4) Rotate cluster tensor from (x,z) to (e0,e1) system
+  //    and calculate moments M3x and M4z
+
+  Float_t cosPhi = e0.X();
+  Float_t sinPhi = e0.Y();
+
+  Float_t xiPHOS ;
+  Float_t ziPHOS ;
+  Double_t dx3, dz3, dz4;
+  wtot = 0.;
+  x    = 0.;
+  z    = 0.;
+  dxx  = 0.;
+  dzz  = 0.;
+  dxz  = 0.;
+  dx3  = 0.;
+  dz3  = 0.;
+  dz4  = 0.;
+  for(iDigit=0; iDigit<fMulDigit; iDigit++) {
+    digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit])  ;
+    phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+    phosgeom->RelPosInModule(relid, xiPHOS, ziPHOS);
+    xi    = xiPHOS*cosPhi + ziPHOS*sinPhi;
+    zi    = ziPHOS*cosPhi - xiPHOS*sinPhi;
+    w     = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
+    x    += w * xi ;
+    z    += w * zi ; 
+    dxx  += w * xi * xi ;
+    dzz  += w * zi * zi ;
+    dxz  += w * xi * zi ; 
+    dx3  += w * xi * xi * xi;
+    dz3  += w * zi * zi * zi ;
+    dz4  += w * zi * zi * zi * zi ;
+    wtot += w ;
+  }
+  x   /= wtot ;
+  z   /= wtot ;
+  dxx /= wtot ;
+  dzz /= wtot ;
+  dxz /= wtot ;
+  dx3 /= wtot ;
+  dz3 /= wtot ;
+  dz4 /= wtot ;
+  dx3 += -3*dxx*x + 2*x*x*x;
+  dz4 += -4*dz3*z + 6*dzz*z*z -3*z*z*z*z;
+  dxx -= x * x ;
+  dzz -= z * z ;
+  dxz -= x * z ;
+
+  // 5) Find an angle between cluster center vector and eigen vector e0
+
+  Float_t phi = TMath::ACos ((x*e0.X() + z*e0.Y()) / sqrt(x*x + z*z));
+
+  fM2x   = lambda0;
+  fM2z   = lambda1;
+  fM3x   = dx3;
+  fM4z   = dz4;
+  fPhixe = phi;
+
+}
 //____________________________________________________________________________
 void AliPHOSEmcRecPoint::EvalAll(Float_t logWeight, TClonesArray * digits )
 {
@@ -487,10 +639,11 @@ void AliPHOSEmcRecPoint::EvalAll(Float_t logWeight, TClonesArray * digits )
 
   EvalLocalPosition(logWeight, digits) ;
   EvalElipsAxis(logWeight, digits) ;
+  EvalMoments(logWeight, digits) ;
   EvalDispersion(logWeight, digits) ;
   EvalCoreEnergy(logWeight, digits);
   EvalTime(digits) ;
-  AliPHOSRecPoint::EvalAll(logWeight,digits) ;
+  AliPHOSRecPoint::EvalAll(digits) ;
 }
 //____________________________________________________________________________
 void AliPHOSEmcRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digits)
@@ -505,8 +658,7 @@ void AliPHOSEmcRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * dig
   
   AliPHOSDigit * digit ;
 
-  AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
-  AliPHOSGeometry * phosgeom =  (AliPHOSGeometry*)gime->PHOSGeometry();
+  AliPHOSGeometry * phosgeom =  AliPHOSLoader::GetPHOSGeometry();
 
   Int_t iDigit;
 
@@ -679,9 +831,11 @@ void AliPHOSEmcRecPoint::Purify(Float_t threshold){
   fDigitsList = new (Int_t[fMulDigit]) ;
   fEnergyList = new ( Float_t[fMulDigit]) ;
 
+  fAmp = 0.; //Recalculate total energy
   for(Int_t iDigit=0;iDigit< fMulDigit ;iDigit++){
-      fDigitsList[iDigit] = tempo[iDigit];
-      fEnergyList[iDigit] = tempoE[iDigit] ;
+     fDigitsList[iDigit] = tempo[iDigit];
+     fEnergyList[iDigit] = tempoE[iDigit] ;
+     fAmp+=tempoE[iDigit];
   }
       
   delete [] tempo ;
@@ -689,7 +843,7 @@ void AliPHOSEmcRecPoint::Purify(Float_t threshold){
 
 }
 //____________________________________________________________________________
-void AliPHOSEmcRecPoint::Print(Option_t * option) const
+void AliPHOSEmcRecPoint::Print(Option_t *) const
 {
   // Print the list of digits belonging to the cluster