From 5bd0f028e015264d585b60dfea0417e6cba801c0 Mon Sep 17 00:00:00 2001
From: prsnko <Dmitri.Peressounko@cern.ch>
Date: Sat, 12 Apr 2014 14:32:14 +0400
Subject: [PATCH] CoreEnergy calculation added

---
 .../TenderSupplies/AliPHOSTenderSupply.cxx    | 59 ++++++++++++++++++-
 ANALYSIS/TenderSupplies/AliPHOSTenderSupply.h |  1 +
 2 files changed, 58 insertions(+), 2 deletions(-)

diff --git a/ANALYSIS/TenderSupplies/AliPHOSTenderSupply.cxx b/ANALYSIS/TenderSupplies/AliPHOSTenderSupply.cxx
index 5ae72bfe1ba..72bfbaf30cc 100644
--- a/ANALYSIS/TenderSupplies/AliPHOSTenderSupply.cxx
+++ b/ANALYSIS/TenderSupplies/AliPHOSTenderSupply.cxx
@@ -309,6 +309,9 @@ void AliPHOSTenderSupply::ProcessEvent()
       cluPHOS.EvalAll(logWeight,vertex);         // recalculate the cluster parameters
       cluPHOS.SetE(CorrectNonlinearity(cluPHOS.E()));// Users's nonlinearity
 
+      Double_t ecore=CoreEnergy(&cluPHOS) ; 
+      ecore=CorrectNonlinearity(ecore) ;
+      
       //Correct Misalignment
 //      CorrectPHOSMisalignment(global,mod) ;
 //      position[0]=global.X() ;
@@ -324,7 +327,9 @@ void AliPHOSTenderSupply::ProcessEvent()
       Float_t  xyz[3];
       cluPHOS.GetPosition(xyz);
       clu->SetPosition(xyz);                       //rec.point position in MARS
-      clu->SetE(cluPHOS.E());                           //total or core particle energy
+      clu->SetE(cluPHOS.E());                      //total particle energy
+      clu->SetMCEnergyFraction(ecore);                            //core particle energy
+      
       //    clu->SetDispersion(cluPHOS.GetDispersion());  //cluster dispersion
       //    ec->SetPID(rp->GetPID()) ;            //array of particle identification
       clu->SetM02(cluPHOS.GetM02()) ;               //second moment M2x
@@ -378,6 +383,10 @@ void AliPHOSTenderSupply::ProcessEvent()
       cluPHOS.EvalAll(logWeight,vertex);         // recalculate the cluster parameters
       cluPHOS.SetE(CorrectNonlinearity(cluPHOS.E()));// Users's nonlinearity
 
+      Double_t ecore=CoreEnergy(&cluPHOS) ; 
+      ecore=CorrectNonlinearity(ecore) ;
+
+      
       //Correct Misalignment
 //      cluPHOS.GetPosition(position);
 //      global.SetXYZ(position[0],position[1],position[2]);
@@ -387,7 +396,8 @@ void AliPHOSTenderSupply::ProcessEvent()
 //      position[2]=global.Z() ;
 
       clu->SetPosition(position);                       //rec.point position in MARS
-      clu->SetE(cluPHOS.E());                           //total or core particle energy
+      clu->SetE(cluPHOS.E());                           //total particle energy
+      clu->SetMCEnergyFraction(ecore);                  //core particle energy
       clu->SetDispersion(cluPHOS.GetDispersion());  //cluster dispersion
       //    ec->SetPID(rp->GetPID()) ;            //array of particle identification
       clu->SetM02(cluPHOS.GetM02()) ;               //second moment M2x
@@ -748,6 +758,51 @@ void AliPHOSTenderSupply::EvalLambdas(AliVCluster * clu, Double_t &m02, Double_t
   }
 
 }
+//____________________________________________________________________________
+Double_t  AliPHOSTenderSupply::CoreEnergy(AliVCluster * clu){  
+  //calculate energy of the cluster in the circle with radius distanceCut around the maximum
+  
+  //Can not use already calculated coordinates?
+  //They have incidence correction...
+  const Double_t distanceCut =3.5 ;
+  const Double_t logWeight=4.5 ;
+  
+  Double32_t * elist = clu->GetCellsAmplitudeFraction() ;  
+// Calculates the center of gravity in the local PHOS-module coordinates
+  Float_t wtot = 0;
+  Double_t xc[100]={0} ;
+  Double_t zc[100]={0} ;
+  Double_t x = 0 ;
+  Double_t z = 0 ;
+  Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ;
+  for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
+    Int_t relid[4] ;
+    Float_t xi ;
+    Float_t zi ;
+    fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
+    fPHOSGeo->RelPosInModule(relid, xi, zi);
+    xc[iDigit]=xi ;
+    zc[iDigit]=zi ;
+    if (clu->E()>0 && elist[iDigit]>0) {
+      Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
+      x    += xc[iDigit] * w ;
+      z    += zc[iDigit] * w ;
+      wtot += w ;
+    }
+  }
+  if (wtot>0) {
+    x /= wtot ;
+    z /= wtot ;
+  }
+  Double_t coreE=0. ;
+  for(Int_t iDigit=0; iDigit < mulDigit; iDigit++) {
+    Double_t distance = TMath::Sqrt((xc[iDigit]-x)*(xc[iDigit]-x)+(zc[iDigit]-z)*(zc[iDigit]-z)) ;
+    if(distance < distanceCut)
+      coreE += elist[iDigit] ;
+  }
+  //Apply non-linearity correction
+  return coreE ;
+}
 //________________________________________________________________________
 Double_t AliPHOSTenderSupply::EvalTOF(AliVCluster * clu,AliVCaloCells * cells){ 
   //Evaluate TOF of the cluster after re-calibration
diff --git a/ANALYSIS/TenderSupplies/AliPHOSTenderSupply.h b/ANALYSIS/TenderSupplies/AliPHOSTenderSupply.h
index 7c5460a9ceb..fe9a1c80622 100644
--- a/ANALYSIS/TenderSupplies/AliPHOSTenderSupply.h
+++ b/ANALYSIS/TenderSupplies/AliPHOSTenderSupply.h
@@ -60,6 +60,7 @@ protected:
   Bool_t IsGoodChannel(Int_t mod, Int_t ix, Int_t iz) ;
   void   CorrectPHOSMisalignment(TVector3 & globalPos, Int_t module);
   void   EvalLambdas(AliVCluster * clu, Double_t &m02, Double_t &m20) ;
+  Double_t CoreEnergy(AliVCluster * clu) ;  
   Double_t EvalTOF(AliVCluster * clu,AliVCaloCells * cells); 
   Double_t CalibrateTOF(Double_t tof, Int_t absId, Bool_t isHG); 
   void DistanceToBadChannel(Int_t mod, TVector3 * locPos, Double_t &minDist) ;
-- 
2.43.0