[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 "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(AliPHOSPIDv1 * pid){
  //apply nonlinearity correction same as in AliPHOSPIDv1.
  SetE(pid->GetCalibratedEnergy(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) ;
  SetPosition(0, gps[0]) ;
  SetPosition(1, gps[1]) ;  
  SetPosition(2, gps[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"
	33	#include "AliPHOSAodCluster.h"
	34	#include "AliPHOSCalibData.h"
	35	#include "AliAODCaloCells.h"
	36
	37	ClassImp(AliPHOSAodCluster)
	38
	39	//____________________________________________________________________________
	40	AliPHOSAodCluster::AliPHOSAodCluster() :
	41	AliAODCaloCluster(),fRecalibrated(0)
	42	{
	43	// ctor
	44	}
	45	//____________________________________________________________________________
	46	AliPHOSAodCluster::AliPHOSAodCluster(const AliAODCaloCluster & clu) :
	47	AliAODCaloCluster(clu),fRecalibrated(0)
	48	{
	49	// cpy ctor
	50	}
	51
	52	//____________________________________________________________________________
	53	AliPHOSAodCluster::~AliPHOSAodCluster()
	54	{
	55	// dtor
	56	}
	57	//____________________________________________________________________________
	58	void AliPHOSAodCluster::Recalibrate(AliPHOSCalibData * calibData,AliAODCaloCells *phsCells){
	59	//If not done yet, apply recalibration coefficients to energies list
	60	//NOTE that after recalibration fCellsAmpFraction contains not FRACTION but FULL energy
	61
	62	if(fRecalibrated)
	63	return ;
	64
65	if(!calibData)
66	return ;
67
68	AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ;
69	if(!phosgeom)
70	AliFatal("AliPHOSGeometry was not contructed\n") ;
71
72	Double32_t * cellsAmpFraction = GetCellsAmplitudeFraction();
73
74	for(Int_t i=0; i < GetNCells(); i++){
75	Int_t relId[4];
76	phosgeom->AbsToRelNumbering(GetCellAbsId(i),relId) ;
77	Int_t module = relId[0];
78	Int_t column = relId[3];
79	Int_t row = relId[2];
80	Double_t energy = phsCells->GetCellAmplitude(GetCellAbsId(i)) ;
81	cellsAmpFraction[i]=energycalibData->GetADCchannelEmc(module,column,row);
82	}
83
84	SetCellsAmplitudeFraction(cellsAmpFraction);
85	fRecalibrated=kTRUE;
86	}
87	//____________________________________________________________________________
88	void AliPHOSAodCluster::EvalAll(Float_t logWeight, TVector3 &vtx){
89	//If recalibrated - recalculate all cluster parameters
90	if(!fRecalibrated)
91	return ;
92
93	EvalEnergy() ; //Energy should be evaluated first
94	EvalCoord(logWeight, vtx) ;
95
96	}
97	//____________________________________________________________________________
98	void AliPHOSAodCluster::EvalEnergy(){
99	if(!fRecalibrated) // no need to recalibrate
100	return ;
101
102	Float_t energy=0. ;
103	for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
104	energy+=GetCellAmplitudeFraction(iDigit) ;
105	}
106
107	SetE(energy);
108	//Correct for nonlinearity later
109
110	}
111	//____________________________________________________________________________
112	void AliPHOSAodCluster::EnergyCorrection(AliPHOSPIDv1 * pid){
113	//apply nonlinearity correction same as in AliPHOSPIDv1.
114	SetE(pid->GetCalibratedEnergy(E())) ;
115	}
116	//____________________________________________________________________________
117	void AliPHOSAodCluster::EvalPID(AliPHOSPIDv1 * /pid/){
118	//re-evaluate identification parameters
119	// pid->CalculatePID(E(),GetDispersion(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
120	// pid->CalculatePID(E(),GetDispersion(),GetM20(),GetM02(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
121	}
122	//____________________________________________________________________________
123	void AliPHOSAodCluster::EvalCoord(Float_t logWeight, TVector3 &vtx)
124	{
125	// Calculates new center of gravity in the local PHOS-module coordinates
126	// and tranfers into global ALICE coordinates
127	// Calculates Dispersion and main axis
128	if(!fRecalibrated) // no need to recalibrate
129	return ;
130
131	Float_t wtot = 0. ;
132	Int_t relid[4] ;
133	Int_t phosMod=0 ;
134	Float_t xMean = 0. ;
135	Float_t zMean = 0. ;
136
137	AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ;
138	if(!phosgeom)
139	AliFatal("AliPHOSGeometry was not contructed\n") ;
140
141	for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
142	Float_t xi ;
143	Float_t zi ;
144	phosgeom->AbsToRelNumbering(GetCellAbsId(iDigit), relid) ;
145	phosgeom->RelPosInModule(relid, xi, zi);
146	phosMod=relid[0] ;
147	Double_t ei=GetCellAmplitudeFraction(iDigit) ;
148	if (E() > 0 && ei > 0) {
149	Float_t w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
150	xMean+= xi * w ;
151	zMean+= zi * w ;
152	wtot += w ;
153	}
154	else
155	AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
156	}
157	if (wtot>0) {
158	xMean /= wtot ;
159	zMean /= wtot ;
160	}
161	else
162	AliError(Form("Wrong weight %f\n", wtot));
163
164
165	// Calculates the dispersion and second momenta
166	Double_t d=0. ;
167	Double_t dxx = 0.;
168	Double_t dzz = 0.;
169	Double_t dxz = 0.;
170	for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
171	Float_t xi ;
172	Float_t zi ;
173	phosgeom->AbsToRelNumbering(GetCellAbsId(iDigit), relid) ;
174	phosgeom->RelPosInModule(relid, xi, zi);
175	Double_t ei=GetCellAmplitudeFraction(iDigit) ;
176	if (E() > 0 && ei > 0) {
177	Float_t w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
178	d += w((xi-xMean)(xi-xMean) + (zi-zMean)*(zi-zMean) ) ;
179	dxx += w * xi * xi ;
180	dzz += w * zi * zi ;
181	dxz += w * xi * zi ;
182	}
183	else
184	AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
185	}
186
187	if (wtot>0) {
188	d /= wtot ;
189	dxx /= wtot ;
190	dzz /= wtot ;
191	dxz /= wtot ;
192	dxx -= xMean * xMean ;
193	dzz -= zMean * zMean ;
194	dxz -= xMean * zMean ;
195	SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )) ;
196	SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
197	}
198	else{
199	AliError(Form("Wrong weight %f\n", wtot));
200	d=0. ;
201	SetM20(0.) ;
202	SetM02(0.) ;
203	}
204
205	if (d>=0)
206	SetDispersion(TMath::Sqrt(d)) ;
207	else
208	SetDispersion(0) ;
209
210
211	// Correction for the depth of the shower starting point (TDR p 127)
212	Float_t para = 0.925 ;
213	Float_t parb = 6.52 ;
214
215	TVector3 vInc ;
216	phosgeom->GetIncidentVector(vtx,phosMod,xMean,zMean,vInc) ;
217
218	Float_t depthx = 0.;
219	Float_t depthz = 0.;
220	if (E()>0&&vInc.Y()!=0.) {
221	depthx = ( para * TMath::Log(E()) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ;
222	depthz = ( para * TMath::Log(E()) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ;
223	}
224	else
225	AliError(Form("Wrong amplitude %f\n", E()));
226
227	xMean-= depthx ;
228	zMean-= depthz ;
229
230	//Go to the global system
231	TVector3 gps ;
232	phosgeom->Local2Global(phosMod, xMean, zMean, gps) ;
233	SetPosition(0, gps[0]) ;
234	SetPosition(1, gps[1]) ;
235	SetPosition(2, gps[2]) ;
236	}