[u/mrichter/AliRoot.git] / EMCAL / AliEMCALAodCluster.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 EMCAL to recalculate cluster 
//  parameters in case of recalibration.
//  Copy-paste from methods in AliEMCALRecPoint.
//*--
//*-- Author: Dmitri Peressounko (RRC KI) for PHOS
//*-- Adapted for EMCAL: Gustavo Conesa (INFN-LNF)

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

// --- Standard library ---

// --- AliRoot header files ---
#include "AliLog.h" 
#include "AliEMCALGeometry.h" 
#include "AliEMCALPID.h" 
#include "AliEMCALAodCluster.h" 
#include "AliEMCALCalibData.h"
#include "AliAODCaloCells.h"

ClassImp(AliEMCALAodCluster)

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

//____________________________________________________________________________
AliEMCALAodCluster::~AliEMCALAodCluster()
{
  // dtor
}
//____________________________________________________________________________
void AliEMCALAodCluster::Recalibrate(AliEMCALCalibData * calibData, AliAODCaloCells *emcCells, TString emcalGeoName){
  //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 ;

  AliEMCALGeometry * emcalgeo =  AliEMCALGeometry::GetInstance(emcalGeoName) ;
  if(emcalgeo){
	
  Double32_t * cellsAmpFraction = GetCellsAmplitudeFraction(); 
  Int_t iSupMod = -1;
  Int_t iTower  = -1;
  Int_t iIphi   = -1;
  Int_t iIeta   = -1;
  Int_t iphi    = -1;
  Int_t ieta    = -1;
  	
  for(Int_t i=0; i < GetNCells(); i++){
    
    //Get from the absid the supermodule, tower and eta/phi numbers
    emcalgeo->GetCellIndex(GetCellAbsId(i),iSupMod,iTower,iIphi,iIeta); 
    //Gives SuperModule and Tower numbers
    emcalgeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
					  iIphi, iIeta,iphi,ieta);   	
    
    Double_t energy = emcCells->GetCellAmplitude(GetCellAbsId(i)) ;
    //AliDebug(2,Form("Recalibrate: cell %f, calib %f, fraction %f\n",energy,calibData->GetADCchannel(iSupMod,ieta,iphi),cellsAmpFraction[i]));
	if(cellsAmpFraction[i]< 1e-4) cellsAmpFraction[i] = 1; //Unfolding off, took all the cell energy 
	cellsAmpFraction[i]*=energy*calibData->GetADCchannel(iSupMod,ieta,iphi);
  }
	
  SetCellsAmplitudeFraction(cellsAmpFraction);
  fRecalibrated=kTRUE; 
    
  }
  else  AliFatal("AliEMCALGeometry was not constructed\n") ;

}
//____________________________________________________________________________
void  AliEMCALAodCluster::EvalAll(Float_t logWeight, TString geoname){
    //If recalibrated - recalculate all cluster parameters
  if(!fRecalibrated)
    return ;
  //printf("EvalAll e org %f\n",E());
  EvalEnergy() ; //Energy should be evaluated first
  //printf("EvalAll e2 %f\n",E());
  EvalPositionAndShowerShape(logWeight, geoname) ;
  //printf("EvalAll e3 %f\n",E());
  EvalPID() ; //Should be evaluated after energy and shower shape recalculation
  //printf("EvalAll e4 %f\n",E());
}
//____________________________________________________________________________
void AliEMCALAodCluster::EvalEnergy(){
  //Evaluate energy
  if(!fRecalibrated) // no need to recalibrate
    return ;
    
  Float_t energy=0. ;
  for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
    energy+=GetCellAmplitudeFraction(iDigit) ;
  }
  //printf("EvalEnergy: e %f\n", energy);
  SetE(energy);
  
   
}
////____________________________________________________________________________
//void AliEMCALAodCluster::EnergyCorrection(AliEMCALPID * pid){
//  //apply nonlinearity correction same as in AliEMCALPID.
//  SetE(pid->GetCalibratedEnergy(E())) ;
//}

//____________________________________________________________________________
void AliEMCALAodCluster::EvalPID(){           
	
  //re-evaluate identification parameters
//  pid->CalculatePID(E(),GetDispersion(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;  
//  pid->CalculatePID(E(),GetDispersion(),GetM20(),GetM02(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;

  //With bayesian
  AliEMCALPID *pid = new AliEMCALPID(kFALSE);
  pid->SetLowFluxParam(); // Need to be fixed
  Float_t pidlist[AliPID::kSPECIESN+1];
	for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = pid->GetPIDFinal(i);	
  SetPIDFromESD(pidlist);
	
}

//____________________________________________________________________________
void AliEMCALAodCluster::EvalPositionAndShowerShape(Float_t logWeight, TString emcalGeoName)
{
  // Calculates new center of gravity in the local EMCAL-module coordinates 
  // and tranfers into global ALICE coordinates
  // Calculates Dispersion and main axis
  if(!fRecalibrated) // no need to recalibrate
    return ;
  
  Int_t nstat  = 0;
  Float_t wtot = 0. ;

  Int_t iSupMod = -1;
  Int_t iTower  = -1;
  Int_t iIphi   = -1;
  Int_t iIeta   = -1;
  Int_t iphi    = -1;
  Int_t ieta    = -1;
  Double_t etai = -1.;
  Double_t phii = -1.;
	
  Double_t clXYZ[3] ={0.,0.,0.}; 
  Double_t xyzi[3]  ={0.,0.,0.};
 
  Double_t d     = 0.;
  Double_t dxx   = 0.;
  Double_t dzz   = 0.;
  Double_t dxz   = 0.;  
  Double_t xmean = 0.;
  Double_t zmean = 0.;

  AliEMCALGeometry * emcalgeo =  AliEMCALGeometry::GetInstance(emcalGeoName) ;
  if(emcalgeo){
	
  Double_t dist  = TmaxInCm(E(),0);
  for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
	
	//Get from the absid the supermodule, tower and eta/phi numbers
	emcalgeo->GetCellIndex(GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); 
	emcalgeo->RelPosCellInSModule(GetCellAbsId(iDigit), dist, xyzi[0], xyzi[1], xyzi[2]);
	emcalgeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);

    Double_t ei = GetCellAmplitudeFraction(iDigit) ;
    if (E() > 0 && ei > 0) {
		Float_t w = ei;  
		if(logWeight > 0) w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
		etai=(Double_t)ieta;
		phii=(Double_t)iphi;		
		if(w > 0.0) {
			wtot += w ;
			nstat++;		
			for(Int_t i = 0; i < 3; i++ ) clXYZ[i]    += (w*xyzi[i]);
	
			//Shower shape
			dxx  += w * etai * etai ;
			xmean+= w * etai ;
			dzz  += w * phii * phii ;
			zmean+= w * phii ; 
			dxz  += w * etai * phii ; 
		}
    }
    else
      AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
  }

  //Normalize to the weight	
  if (wtot > 0) {
		for(Int_t i=0; i<3; i++ ) clXYZ[i] /= wtot;
		xmean /= wtot ;
		zmean /= wtot ;
  }
  else
    AliError(Form("Wrong weight %f\n", wtot));

  //Put cluster position in the global system
  TVector3 gpos ;
  emcalgeo->GetGlobal(clXYZ, gpos, iSupMod);

  SetPositionAt(gpos[0],0) ;
  SetPositionAt(gpos[1],1) ;  
  SetPositionAt(gpos[2],2) ;

  //Calculate dispersion	
  for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
		//Get from the absid the supermodule, tower and eta/phi numbers
		emcalgeo->GetCellIndex(GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); 
		emcalgeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
		
		Double_t ei=GetCellAmplitudeFraction(iDigit) ;
		if (E() > 0 && ei > 0) {
			Float_t w = ei;  
			if(logWeight > 0) w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
			etai=(Double_t)ieta;
			phii=(Double_t)iphi;		
			if(w > 0.0)  d +=  w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean)); 
		}
		else
			AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
  }
	
  //Normalize to the weigth and set shower shape parameters
  if (wtot > 0 && nstat > 1) {
    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{
    d=0. ;
    SetM20(0.) ;
    SetM02(0.) ;
  }	
	
  if (d>=0)
	  SetDispersion(TMath::Sqrt(d)) ;
  else    
	  SetDispersion(0) ;
  }
  else {
    AliFatal("AliEMCALGeometry was not constructed\n") ;
  }


}

//_____________________________________________________________________
Double_t AliEMCALAodCluster::TmaxInCm(const Double_t e , const Int_t key) const
{ 
	// e energy in GeV)
	// key  =  0(gamma, default)
	//     !=  0(electron)
	static Double_t ca = 4.82;  // shower max parameter - first guess; ca=TMath::Log(1000./8.07)
	static Double_t x0 = 1.23;  // radiation lenght (cm)
	static Double_t tmax = 0.;   // position of electromagnetic shower max in cm
	
	tmax = 0.0;
	if(e>0.1) {
		tmax = TMath::Log(e) + ca;
		if      (key==0) tmax += 0.5; 
		else             tmax -= 0.5;
		tmax *= x0; // convert to cm
	}
	return tmax;
}
Commit	Line	Data
0c5b726e	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 EMCAL to recalculate cluster
	18	// parameters in case of recalibration.
	19	// Copy-paste from methods in AliEMCALRecPoint.
	20	//*--
	21	//*-- Author: Dmitri Peressounko (RRC KI) for PHOS
	22	//*-- Adapted for EMCAL: Gustavo Conesa (INFN-LNF)
	23
	24	// --- ROOT system ---
	25	#include "TVector3.h"
	26	#include "TMath.h"
	27
	28	// --- Standard library ---
	29
	30	// --- AliRoot header files ---
	31	#include "AliLog.h"
	32	#include "AliEMCALGeometry.h"
	33	#include "AliEMCALPID.h"
	34	#include "AliEMCALAodCluster.h"
	35	#include "AliEMCALCalibData.h"
	36	#include "AliAODCaloCells.h"
	37
	38	ClassImp(AliEMCALAodCluster)
	39
	40	//____________________________________________________________________________
	41	AliEMCALAodCluster::AliEMCALAodCluster() :
	42	AliAODCaloCluster(),fRecalibrated(0)
	43	{
	44	// ctor
	45	}
	46	//____________________________________________________________________________
	47	AliEMCALAodCluster::AliEMCALAodCluster(const AliAODCaloCluster & clu) :
	48	AliAODCaloCluster(clu),fRecalibrated(0)
	49	{
	50	// cpy ctor
	51	}
	52
	53	//____________________________________________________________________________
	54	AliEMCALAodCluster::~AliEMCALAodCluster()
	55	{
	56	// dtor
	57	}
	58	//____________________________________________________________________________
	59	void AliEMCALAodCluster::Recalibrate(AliEMCALCalibData * calibData, AliAODCaloCells *emcCells, TString emcalGeoName){
	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	AliEMCALGeometry * emcalgeo = AliEMCALGeometry::GetInstance(emcalGeoName) ;
7e1d9a9b	70	if(emcalgeo){
0c5b726e	71
	72	Double32_t * cellsAmpFraction = GetCellsAmplitudeFraction();
	73	Int_t iSupMod = -1;
	74	Int_t iTower = -1;
	75	Int_t iIphi = -1;
	76	Int_t iIeta = -1;
	77	Int_t iphi = -1;
	78	Int_t ieta = -1;
	79
	80	for(Int_t i=0; i < GetNCells(); i++){
	81
	82	//Get from the absid the supermodule, tower and eta/phi numbers
	83	emcalgeo->GetCellIndex(GetCellAbsId(i),iSupMod,iTower,iIphi,iIeta);
	84	//Gives SuperModule and Tower numbers
	85	emcalgeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
	86	iIphi, iIeta,iphi,ieta);
	87
	88	Double_t energy = emcCells->GetCellAmplitude(GetCellAbsId(i)) ;
60ea57cf	89	//AliDebug(2,Form("Recalibrate: cell %f, calib %f, fraction %f\n",energy,calibData->GetADCchannel(iSupMod,ieta,iphi),cellsAmpFraction[i]));
60ea57cf	90	if(cellsAmpFraction[i]< 1e-4) cellsAmpFraction[i] = 1; //Unfolding off, took all the cell energy
0c5b726e	91	cellsAmpFraction[i]=energycalibData->GetADCchannel(iSupMod,ieta,iphi);
	92	}
	93
	94	SetCellsAmplitudeFraction(cellsAmpFraction);
	95	fRecalibrated=kTRUE;
7e1d9a9b	96
	97	}
	98	else AliFatal("AliEMCALGeometry was not constructed\n") ;
	99
0c5b726e	100	}
	101	//____________________________________________________________________________
	102	void AliEMCALAodCluster::EvalAll(Float_t logWeight, TString geoname){
	103	//If recalibrated - recalculate all cluster parameters
	104	if(!fRecalibrated)
	105	return ;
	106	//printf("EvalAll e org %f\n",E());
	107	EvalEnergy() ; //Energy should be evaluated first
	108	//printf("EvalAll e2 %f\n",E());
	109	EvalPositionAndShowerShape(logWeight, geoname) ;
	110	//printf("EvalAll e3 %f\n",E());
	111	EvalPID() ; //Should be evaluated after energy and shower shape recalculation
	112	//printf("EvalAll e4 %f\n",E());
	113	}
	114	//____________________________________________________________________________
	115	void AliEMCALAodCluster::EvalEnergy(){
	116	//Evaluate energy
	117	if(!fRecalibrated) // no need to recalibrate
	118	return ;
	119
	120	Float_t energy=0. ;
	121	for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
	122	energy+=GetCellAmplitudeFraction(iDigit) ;
	123	}
	124	//printf("EvalEnergy: e %f\n", energy);
	125	SetE(energy);
	126
	127
	128	}
	129	////____________________________________________________________________________
	130	//void AliEMCALAodCluster::EnergyCorrection(AliEMCALPID * pid){
	131	// //apply nonlinearity correction same as in AliEMCALPID.
	132	// SetE(pid->GetCalibratedEnergy(E())) ;
	133	//}
	134
	135	//____________________________________________________________________________
	136	void AliEMCALAodCluster::EvalPID(){
	137
	138	//re-evaluate identification parameters
	139	// pid->CalculatePID(E(),GetDispersion(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
	140	// pid->CalculatePID(E(),GetDispersion(),GetM20(),GetM02(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
	141
	142	//With bayesian
	143	AliEMCALPID *pid = new AliEMCALPID(kFALSE);
	144	pid->SetLowFluxParam(); // Need to be fixed
	145	Float_t pidlist[AliPID::kSPECIESN+1];
	146	for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = pid->GetPIDFinal(i);
	147	SetPIDFromESD(pidlist);
	148
	149	}
	150
	151	//____________________________________________________________________________
	152	void AliEMCALAodCluster::EvalPositionAndShowerShape(Float_t logWeight, TString emcalGeoName)
	153	{
	154	// Calculates new center of gravity in the local EMCAL-module coordinates
	155	// and tranfers into global ALICE coordinates
	156	// Calculates Dispersion and main axis
	157	if(!fRecalibrated) // no need to recalibrate
	158	return ;
	159
	160	Int_t nstat = 0;
	161	Float_t wtot = 0. ;
	162
0c5b726e	163	Int_t iSupMod = -1;
	164	Int_t iTower = -1;
	165	Int_t iIphi = -1;
	166	Int_t iIeta = -1;
	167	Int_t iphi = -1;
	168	Int_t ieta = -1;
	169	Double_t etai = -1.;
	170	Double_t phii = -1.;
	171
	172	Double_t clXYZ[3] ={0.,0.,0.};
	173	Double_t xyzi[3] ={0.,0.,0.};
	174
	175	Double_t d = 0.;
	176	Double_t dxx = 0.;
	177	Double_t dzz = 0.;
	178	Double_t dxz = 0.;
	179	Double_t xmean = 0.;
	180	Double_t zmean = 0.;
	181
	182	AliEMCALGeometry * emcalgeo = AliEMCALGeometry::GetInstance(emcalGeoName) ;
7e1d9a9b	183	if(emcalgeo){
0c5b726e	184
60ea57cf	185	Double_t dist = TmaxInCm(E(),0);
0c5b726e	186	for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
7cfcebd3	187
0c5b726e	188	//Get from the absid the supermodule, tower and eta/phi numbers
0c5b726e	189	emcalgeo->GetCellIndex(GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
7cfcebd3	190	emcalgeo->RelPosCellInSModule(GetCellAbsId(iDigit), dist, xyzi[0], xyzi[1], xyzi[2]);
0c5b726e	191	emcalgeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
	192
	193	Double_t ei = GetCellAmplitudeFraction(iDigit) ;
	194	if (E() > 0 && ei > 0) {
	195	Float_t w = ei;
	196	if(logWeight > 0) w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
	197	etai=(Double_t)ieta;
	198	phii=(Double_t)iphi;
	199	if(w > 0.0) {
	200	wtot += w ;
	201	nstat++;
	202	for(Int_t i = 0; i < 3; i++ ) clXYZ[i] += (w*xyzi[i]);
	203
	204	//Shower shape
	205	dxx += w * etai * etai ;
	206	xmean+= w * etai ;
	207	dzz += w * phii * phii ;
	208	zmean+= w * phii ;
	209	dxz += w * etai * phii ;
	210	}
	211	}
	212	else
	213	AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
	214	}
	215
	216	//Normalize to the weight
	217	if (wtot > 0) {
	218	for(Int_t i=0; i<3; i++ ) clXYZ[i] /= wtot;
	219	xmean /= wtot ;
	220	zmean /= wtot ;
	221	}
	222	else
	223	AliError(Form("Wrong weight %f\n", wtot));
	224
	225	//Put cluster position in the global system
	226	TVector3 gpos ;
	227	emcalgeo->GetGlobal(clXYZ, gpos, iSupMod);
	228
c8fe2783	229	SetPositionAt(gpos[0],0) ;
	230	SetPositionAt(gpos[1],1) ;
	231	SetPositionAt(gpos[2],2) ;
0c5b726e	232
	233	//Calculate dispersion
	234	for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
	235	//Get from the absid the supermodule, tower and eta/phi numbers
	236	emcalgeo->GetCellIndex(GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
	237	emcalgeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
	238
	239	Double_t ei=GetCellAmplitudeFraction(iDigit) ;
	240	if (E() > 0 && ei > 0) {
	241	Float_t w = ei;
	242	if(logWeight > 0) w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
	243	etai=(Double_t)ieta;
	244	phii=(Double_t)iphi;
	245	if(w > 0.0) d += w((etai-xmean)(etai-xmean)+(phii-zmean)*(phii-zmean));
	246	}
	247	else
	248	AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
	249	}
	250
	251	//Normalize to the weigth and set shower shape parameters
	252	if (wtot > 0 && nstat > 1) {
	253	d /= wtot ;
	254	dxx /= wtot ;
	255	dzz /= wtot ;
	256	dxz /= wtot ;
	257	dxx -= xmean * xmean ;
	258	dzz -= zmean * zmean ;
	259	dxz -= xmean * zmean ;
	260	SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )) ;
	261	SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
	262	}
	263	else{
	264	d=0. ;
	265	SetM20(0.) ;
	266	SetM02(0.) ;
	267	}
	268
	269	if (d>=0)
	270	SetDispersion(TMath::Sqrt(d)) ;
	271	else
	272	SetDispersion(0) ;
7e1d9a9b	273	}
	274	else {
	275	AliFatal("AliEMCALGeometry was not constructed\n") ;
	276	}
	277
0c5b726e	278
	279	}
	280
	281	//_____________________________________________________________________
	282	Double_t AliEMCALAodCluster::TmaxInCm(const Double_t e , const Int_t key) const
	283	{
	284	// e energy in GeV)
	285	// key = 0(gamma, default)
	286	// != 0(electron)
	287	static Double_t ca = 4.82; // shower max parameter - first guess; ca=TMath::Log(1000./8.07)
	288	static Double_t x0 = 1.23; // radiation lenght (cm)
	289	static Double_t tmax = 0.; // position of electromagnetic shower max in cm
	290
	291	tmax = 0.0;
	292	if(e>0.1) {
	293	tmax = TMath::Log(e) + ca;
	294	if (key==0) tmax += 0.5;
	295	else tmax -= 0.5;
	296	tmax *= x0; // convert to cm
	297	}
	298	return tmax;
	299	}
	300