[u/mrichter/AliRoot.git] / PHOS / AliPHOSAodCluster.cxx

/**************************************************************************
 * 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.                  *
 **************************************************************************/

//_________________________________________________________________________
//  AliAODCaloCluster extension for PHOS to recalculate cluster 
//  parameters in case of recalibration.
//*--
//*-- Author: Dmitri Peressounko (RRC KI)


// --- ROOT system ---
#include "TVector3.h"
#include "TMath.h"

// --- Standard library ---

// --- AliRoot header files ---
#include "AliLog.h" 
#include "AliPHOSGeometry.h" 
#include "AliPHOSPIDv1.h" 
#include "AliPHOSReconstructor.h" 
#include "AliPHOSAodCluster.h" 
#include "AliPHOSCalibData.h"
#include "AliAODCaloCells.h"

ClassImp(AliPHOSAodCluster)

//____________________________________________________________________________
AliPHOSAodCluster::AliPHOSAodCluster() : 
  AliAODCaloCluster(),fRecalibrated(0)
{
  // ctor
}
//____________________________________________________________________________
AliPHOSAodCluster::AliPHOSAodCluster(const AliAODCaloCluster & clu) : 
  AliAODCaloCluster(clu),fRecalibrated(0)
{
  // cpy ctor
}

//____________________________________________________________________________
AliPHOSAodCluster::~AliPHOSAodCluster()
{
  // dtor
}
//____________________________________________________________________________
void AliPHOSAodCluster::Recalibrate(AliPHOSCalibData * calibData,AliAODCaloCells *phsCells){
  //If not done yet, apply recalibration coefficients to energies list
  //NOTE that after recalibration fCellsAmpFraction contains not FRACTION but FULL energy 
  
  if(fRecalibrated)
   return ;
  
  if(!calibData)
    return ;

  AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
  if(!phosgeom)
    AliFatal("AliPHOSGeometry was not contructed\n") ;
	
  Double32_t * cellsAmpFraction = GetCellsAmplitudeFraction(); 
	
  for(Int_t i=0; i < GetNCells(); i++){
    Int_t relId[4];
    phosgeom->AbsToRelNumbering(GetCellAbsId(i),relId) ;
    Int_t   module = relId[0];
    Int_t   column = relId[3];
    Int_t   row    = relId[2];
    Double_t energy = phsCells->GetCellAmplitude(GetCellAbsId(i)) ;
    cellsAmpFraction[i]*=energy*calibData->GetADCchannelEmc(module,column,row);
  }
	
  SetCellsAmplitudeFraction(cellsAmpFraction);
  fRecalibrated=kTRUE; 
}
//____________________________________________________________________________
void  AliPHOSAodCluster::EvalAll(Float_t logWeight, TVector3 &vtx){
    //If recalibrated - recalculate all cluster parameters
  if(!fRecalibrated)
    return ;

  EvalEnergy() ; //Energy should be evaluated first
  EvalCoord(logWeight, vtx) ;
  
}
//____________________________________________________________________________
void AliPHOSAodCluster::EvalEnergy(){
  if(!fRecalibrated) // no need to recalibrate
    return ;
    
  Float_t energy=0. ;
  for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
    energy+=GetCellAmplitudeFraction(iDigit) ;
  }
  
  SetE(energy);
  //Correct for nonlinearity later
   
}
//____________________________________________________________________________
void AliPHOSAodCluster::EnergyCorrection(){
  //apply nonlinearity correction same as in AliPHOSPIDv1.
  SetE(AliPHOSReconstructor::CorrectNonlinearity(E())) ;
}
//____________________________________________________________________________
void AliPHOSAodCluster::EvalPID(AliPHOSPIDv1 * /*pid*/){           
  //re-evaluate identification parameters
//  pid->CalculatePID(E(),GetDispersion(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;  
//  pid->CalculatePID(E(),GetDispersion(),GetM20(),GetM02(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
}
//____________________________________________________________________________
void AliPHOSAodCluster::EvalCoord(Float_t logWeight, TVector3 &vtx)
{
  // Calculates new center of gravity in the local PHOS-module coordinates 
  // and tranfers into global ALICE coordinates
  // Calculates Dispersion and main axis
  if(!fRecalibrated) // no need to recalibrate
    return ;
 
  Float_t wtot = 0. ;
  Int_t relid[4] ;
  Int_t phosMod=0 ;
  Float_t xMean = 0. ;
  Float_t zMean = 0. ;

  AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
  if(!phosgeom)
    AliFatal("AliPHOSGeometry was not contructed\n") ;

  for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
    Float_t xi ;
    Float_t zi ;
    phosgeom->AbsToRelNumbering(GetCellAbsId(iDigit), relid) ;
    phosgeom->RelPosInModule(relid, xi, zi);
    phosMod=relid[0] ;
    Double_t ei=GetCellAmplitudeFraction(iDigit) ;
    if (E() > 0 && ei > 0) {
      Float_t w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
      xMean+= xi * w ;
      zMean+= zi * w ;
      wtot += w ;
    }
    else
      AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
  }
  if (wtot>0) {
    xMean /= wtot ;
    zMean /= wtot ;
  }
  else
    AliError(Form("Wrong weight %f\n", wtot));


// Calculates the dispersion and second momenta
  Double_t d=0. ;
  Double_t dxx  = 0.;
  Double_t dzz  = 0.;
  Double_t dxz  = 0.;
  for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
    Float_t xi ;
    Float_t zi ;
    phosgeom->AbsToRelNumbering(GetCellAbsId(iDigit), relid) ;
    phosgeom->RelPosInModule(relid, xi, zi);
    Double_t ei=GetCellAmplitudeFraction(iDigit) ;
    if (E() > 0 && ei > 0) {
      Float_t w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
      d += w*((xi-xMean)*(xi-xMean) + (zi-zMean)*(zi-zMean) ) ; 
      dxx  += w * xi * xi ;
      dzz  += w * zi * zi ;
      dxz  += w * xi * zi ; 
   }
    else
      AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
  }
  
  if (wtot>0) {
    d /= wtot ;
    dxx /= wtot ;
    dzz /= wtot ;
    dxz /= wtot ;
    dxx -= xMean * xMean ;
    dzz -= zMean * zMean ;
    dxz -= xMean * zMean ;
    SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )) ;
    SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
  }
  else{
    AliError(Form("Wrong weight %f\n", wtot));
    d=0. ;
    SetM20(0.) ;
    SetM02(0.) ;
  }

  if (d>=0)
    SetDispersion(TMath::Sqrt(d)) ;
  else    
    SetDispersion(0) ;


  // Correction for the depth of the shower starting point (TDR p 127)  
  Float_t para = 0.925 ; 
  Float_t parb = 6.52 ; 

  TVector3 vInc ;
  phosgeom->GetIncidentVector(vtx,phosMod,xMean,zMean,vInc) ;

  Float_t depthx = 0.; 
  Float_t depthz = 0.;
  if (E()>0&&vInc.Y()!=0.) {
    depthx = ( para * TMath::Log(E()) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ;
    depthz = ( para * TMath::Log(E()) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ;
  }
  else 
    AliError(Form("Wrong amplitude %f\n", E()));

  xMean-= depthx  ;
  zMean-= depthz  ;

  //Go to the global system
  TVector3 gps ;
  phosgeom->Local2Global(phosMod, xMean, zMean, gps) ;
  SetPositionAt(gps[0],0) ;
  SetPositionAt(gps[1],1) ;  
  SetPositionAt(gps[2],2) ;
}
Commit	Line	Data
914d1bb7	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//_________________________________________________________________________
	17	// AliAODCaloCluster extension for PHOS to recalculate cluster
	18	// parameters in case of recalibration.
	19	//*--
	20	//*-- Author: Dmitri Peressounko (RRC KI)
	21
	22
	23	// --- ROOT system ---
	24	#include "TVector3.h"
	25	#include "TMath.h"
	26
	27	// --- Standard library ---
	28
	29	// --- AliRoot header files ---
	30	#include "AliLog.h"
	31	#include "AliPHOSGeometry.h"
	32	#include "AliPHOSPIDv1.h"
25312a8e	33	#include "AliPHOSReconstructor.h"
914d1bb7	34	#include "AliPHOSAodCluster.h"
	35	#include "AliPHOSCalibData.h"
	36	#include "AliAODCaloCells.h"
	37
	38	ClassImp(AliPHOSAodCluster)
	39
	40	//____________________________________________________________________________
	41	AliPHOSAodCluster::AliPHOSAodCluster() :
	42	AliAODCaloCluster(),fRecalibrated(0)
	43	{
	44	// ctor
	45	}
	46	//____________________________________________________________________________
	47	AliPHOSAodCluster::AliPHOSAodCluster(const AliAODCaloCluster & clu) :
	48	AliAODCaloCluster(clu),fRecalibrated(0)
	49	{
	50	// cpy ctor
	51	}
	52
	53	//____________________________________________________________________________
	54	AliPHOSAodCluster::~AliPHOSAodCluster()
	55	{
	56	// dtor
	57	}
	58	//____________________________________________________________________________
	59	void AliPHOSAodCluster::Recalibrate(AliPHOSCalibData * calibData,AliAODCaloCells *phsCells){
	60	//If not done yet, apply recalibration coefficients to energies list
	61	//NOTE that after recalibration fCellsAmpFraction contains not FRACTION but FULL energy
	62
	63	if(fRecalibrated)
	64	return ;
	65
	66	if(!calibData)
	67	return ;
	68
	69	AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ;
	70	if(!phosgeom)
	71	AliFatal("AliPHOSGeometry was not contructed\n") ;
	72
	73	Double32_t * cellsAmpFraction = GetCellsAmplitudeFraction();
	74
	75	for(Int_t i=0; i < GetNCells(); i++){
	76	Int_t relId[4];
	77	phosgeom->AbsToRelNumbering(GetCellAbsId(i),relId) ;
	78	Int_t module = relId[0];
	79	Int_t column = relId[3];
	80	Int_t row = relId[2];
	81	Double_t energy = phsCells->GetCellAmplitude(GetCellAbsId(i)) ;
	82	cellsAmpFraction[i]=energycalibData->GetADCchannelEmc(module,column,row);
	83	}
	84
	85	SetCellsAmplitudeFraction(cellsAmpFraction);
	86	fRecalibrated=kTRUE;
	87	}
	88	//____________________________________________________________________________
	89	void AliPHOSAodCluster::EvalAll(Float_t logWeight, TVector3 &vtx){
	90	//If recalibrated - recalculate all cluster parameters
	91	if(!fRecalibrated)
	92	return ;
	93
	94	EvalEnergy() ; //Energy should be evaluated first
	95	EvalCoord(logWeight, vtx) ;
	96
	97	}
98	//____________________________________________________________________________
99	void AliPHOSAodCluster::EvalEnergy(){
100	if(!fRecalibrated) // no need to recalibrate
101	return ;
102
103	Float_t energy=0. ;
104	for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
105	energy+=GetCellAmplitudeFraction(iDigit) ;
106	}
107
108	SetE(energy);
109	//Correct for nonlinearity later
110
111	}
112	//____________________________________________________________________________
25312a8e	113	void AliPHOSAodCluster::EnergyCorrection(){
914d1bb7	114	//apply nonlinearity correction same as in AliPHOSPIDv1.
25312a8e	115	SetE(AliPHOSReconstructor::CorrectNonlinearity(E())) ;
914d1bb7	116	}
	117	//____________________________________________________________________________
	118	void AliPHOSAodCluster::EvalPID(AliPHOSPIDv1 * /pid/){
	119	//re-evaluate identification parameters
	120	// pid->CalculatePID(E(),GetDispersion(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
	121	// pid->CalculatePID(E(),GetDispersion(),GetM20(),GetM02(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
	122	}
	123	//____________________________________________________________________________
	124	void AliPHOSAodCluster::EvalCoord(Float_t logWeight, TVector3 &vtx)
	125	{
	126	// Calculates new center of gravity in the local PHOS-module coordinates
	127	// and tranfers into global ALICE coordinates
	128	// Calculates Dispersion and main axis
	129	if(!fRecalibrated) // no need to recalibrate
	130	return ;
	131
	132	Float_t wtot = 0. ;
	133	Int_t relid[4] ;
	134	Int_t phosMod=0 ;
	135	Float_t xMean = 0. ;
	136	Float_t zMean = 0. ;
	137
	138	AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ;
	139	if(!phosgeom)
	140	AliFatal("AliPHOSGeometry was not contructed\n") ;
	141
	142	for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
	143	Float_t xi ;
	144	Float_t zi ;
	145	phosgeom->AbsToRelNumbering(GetCellAbsId(iDigit), relid) ;
	146	phosgeom->RelPosInModule(relid, xi, zi);
	147	phosMod=relid[0] ;
	148	Double_t ei=GetCellAmplitudeFraction(iDigit) ;
	149	if (E() > 0 && ei > 0) {
	150	Float_t w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
	151	xMean+= xi * w ;
	152	zMean+= zi * w ;
	153	wtot += w ;
	154	}
	155	else
	156	AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
	157	}
	158	if (wtot>0) {
	159	xMean /= wtot ;
	160	zMean /= wtot ;
	161	}
	162	else
	163	AliError(Form("Wrong weight %f\n", wtot));
	164
	165
	166	// Calculates the dispersion and second momenta
	167	Double_t d=0. ;
	168	Double_t dxx = 0.;
	169	Double_t dzz = 0.;
	170	Double_t dxz = 0.;
	171	for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
	172	Float_t xi ;
	173	Float_t zi ;
	174	phosgeom->AbsToRelNumbering(GetCellAbsId(iDigit), relid) ;
	175	phosgeom->RelPosInModule(relid, xi, zi);
	176	Double_t ei=GetCellAmplitudeFraction(iDigit) ;
	177	if (E() > 0 && ei > 0) {
	178	Float_t w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
	179	d += w((xi-xMean)(xi-xMean) + (zi-zMean)*(zi-zMean) ) ;
180	dxx += w * xi * xi ;
181	dzz += w * zi * zi ;
182	dxz += w * xi * zi ;
183	}
184	else
185	AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
186	}
187
188	if (wtot>0) {
189	d /= wtot ;
190	dxx /= wtot ;
191	dzz /= wtot ;
192	dxz /= wtot ;
193	dxx -= xMean * xMean ;
194	dzz -= zMean * zMean ;
195	dxz -= xMean * zMean ;
196	SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )) ;
197	SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
198	}
199	else{
200	AliError(Form("Wrong weight %f\n", wtot));
201	d=0. ;
202	SetM20(0.) ;
203	SetM02(0.) ;
204	}
205
206	if (d>=0)
207	SetDispersion(TMath::Sqrt(d)) ;
208	else
209	SetDispersion(0) ;
210
211
212	// Correction for the depth of the shower starting point (TDR p 127)
213	Float_t para = 0.925 ;
214	Float_t parb = 6.52 ;
215
216	TVector3 vInc ;
217	phosgeom->GetIncidentVector(vtx,phosMod,xMean,zMean,vInc) ;
218
219	Float_t depthx = 0.;
220	Float_t depthz = 0.;
221	if (E()>0&&vInc.Y()!=0.) {
222	depthx = ( para * TMath::Log(E()) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ;
223	depthz = ( para * TMath::Log(E()) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ;
224	}
225	else
226	AliError(Form("Wrong amplitude %f\n", E()));
227
228	xMean-= depthx ;
229	zMean-= depthz ;
230
231	//Go to the global system
232	TVector3 gps ;
233	phosgeom->Local2Global(phosMod, xMean, zMean, gps) ;
c8fe2783	234	SetPositionAt(gps[0],0) ;
	235	SetPositionAt(gps[1],1) ;
	236	SetPositionAt(gps[2],2) ;
914d1bb7	237	}