#include "AliEMCALRecParam.h"
#include "AliEMCALReconstructor.h"
#include "AliCDBManager.h"
+#include "AliCaloCalibPedestal.h"
class AliCDBStorage;
#include "AliCDBEntry.h"
fGeom(0),
fDefaultInit(kFALSE),
fToUnfold(kFALSE),
- fNumberOfECAClusters(0),fCalibData(0),
+ fNumberOfECAClusters(0),fCalibData(0),fCaloPed(0),
fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.),
fECAW0(0.),fTimeCut(0.),fMinECut(0.)
{
fGeom(geometry),
fDefaultInit(kFALSE),
fToUnfold(kFALSE),
- fNumberOfECAClusters(0),fCalibData(0),
+ fNumberOfECAClusters(0),fCalibData(0), fCaloPed(0),
fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.),
fECAW0(0.),fTimeCut(0.),fMinECut(0.)
{
}
+//____________________________________________________________________________
+AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
+: AliEMCALClusterizer(),
+fGeom(geometry),
+fDefaultInit(kFALSE),
+fToUnfold(kFALSE),
+fNumberOfECAClusters(0),fCalibData(calib), fCaloPed(caloped),
+fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.),
+fECAW0(0.),fTimeCut(0.),fMinECut(0.)
+{
+ // ctor, geometry and calibration are initialized elsewhere.
+
+ if (!fGeom)
+ AliFatal("Geometry not initialized.");
+
+ if(!gMinuit)
+ gMinuit = new TMinuit(100) ;
+
+}
+
+
//____________________________________________________________________________
AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
{
}
fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
-
+
+ // Check if channel is bad (dead, hot ...), in this case return 0.
+ // Gustavo: 15-12-09 In case of RAW data this selection is already done, but not in simulation.
+ // for the moment keep it here but remember to do the selection at the sdigitizer level
+ // and remove it from here
+ if(fCaloPed->IsBadChannel(iSupMod,ieta,iphi)) {
+ AliDebug(2,Form("Tower from SM %d, ieta %d, iphi %d is BAD!!!",iSupMod,ieta,iphi));
+ return 0;
+ }
+
fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
//Get calibration parameters from file or digitizer default values.
GetCalibrationParameters() ;
-
+ //Get dead channel map from file or digitizer default values.
+ GetCaloCalibPedestal() ;
+
fNumberOfECAClusters = 0;
MakeClusters() ; //only the real clusters
//Evaluate position, dispersion and other RecPoint properties for EC section
for(index = 0; index < fRecPoints->GetEntries(); index++) {
dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->EvalAll(fECAW0,fDigitsArr) ;
+ //For each rec.point set the distance to the nearest bad crystal
+ dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->EvalDistanceToBadChannels(fCaloPed);
}
fRecPoints->Sort() ;
AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
+ fRecPoints->Delete();
+
if(strstr(option,"tim")){
gBenchmark->Stop("EMCALClusterizer");
printf("Exec took %f seconds for Clusterizing",
}
//____________________________________________________________________________
-Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * RecPoint, AliEMCALDigit ** maxAt,
- Float_t* maxAtEnergy,
+Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * recPoint, AliEMCALDigit ** maxAt,
+ const Float_t* maxAtEnergy,
Int_t nPar, Float_t * fitparameters) const
{
// Calls TMinuit to fit the energy distribution of a cluster with several maxima
gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ;
// To set the address of the minimization function
TList * toMinuit = new TList();
- toMinuit->AddAt(RecPoint,0) ;
+ toMinuit->AddAt(recPoint,0) ;
toMinuit->AddAt(fDigitsArr,1) ;
toMinuit->AddAt(fGeom,2) ;
//Check if calibration is stored in data base
- if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
+ if(!fCalibData)
{
AliCDBEntry *entry = (AliCDBEntry*)
AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
}
+//____________________________________________________________________________
+void AliEMCALClusterizerv1::GetCaloCalibPedestal()
+{
+ // Set calibration parameters:
+ // if calibration database exists, they are read from database,
+ // otherwise, they are taken from digitizer.
+ //
+ // It is a user responsilibity to open CDB before reconstruction,
+ // for example:
+ // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
+
+ //Check if calibration is stored in data base
+
+ if(!fCaloPed)
+ {
+ AliCDBEntry *entry = (AliCDBEntry*)
+ AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
+ if (entry) fCaloPed = (AliCaloCalibPedestal*) entry->GetObject();
+ }
+
+ if(!fCaloPed)
+ AliFatal("Pedestal info not found in CDB!");
+
+}
+
+
//____________________________________________________________________________
void AliEMCALClusterizerv1::Init()
{
// Make all memory allocations which can not be done in default constructor.
// Attach the Clusterizer task to the list of EMCAL tasks
- AliRunLoader *rl = AliRunLoader::GetRunLoader();
+ AliRunLoader *rl = AliRunLoader::Instance();
if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL"))
fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
else
fTimeCut = 300e-9 ; // 300 ns time cut (to be tuned)
fCalibData = 0 ;
-
+ fCaloPed = 0 ;
+
const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
if(!recParam) {
AliFatal("Reconstruction parameters for EMCAL not set!");
Int_t index ;
for(index = 0 ; index < numberofNotUnfolded ; index++){
- AliEMCALRecPoint * RecPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;
+ AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;
TVector3 gpos;
Int_t absId;
- RecPoint->GetGlobalPosition(gpos);
+ recPoint->GetGlobalPosition(gpos);
fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);
if(absId > nModulesToUnfold)
break ;
- Int_t nMultipl = RecPoint->GetMultiplicity() ;
+ Int_t nMultipl = recPoint->GetMultiplicity() ;
AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;
Float_t * maxAtEnergy = new Float_t[nMultipl] ;
- Int_t nMax = RecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;
+ Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;
if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
- UnfoldCluster(RecPoint, nMax, maxAt, maxAtEnergy) ;
- fRecPoints->Remove(RecPoint);
+ UnfoldCluster(recPoint, nMax, maxAt, maxAtEnergy) ;
+ fRecPoints->Remove(recPoint);
fRecPoints->Compress() ;
index-- ;
fNumberOfECAClusters-- ;
numberofNotUnfolded-- ;
}
else{
- RecPoint->SetNExMax(1) ; //Only one local maximum
+ recPoint->SetNExMax(1) ; //Only one local maximum
}
delete[] maxAt ;
Float_t xpar=0.,zpar=0.,epar=0. ;
AliEMCALDigit * digit = 0 ;
- Int_t * Digits = iniTower->GetDigitsList() ;
+ Int_t * digitsList = iniTower->GetDigitsList() ;
Int_t iparam ;
Int_t iDigit ;
for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
- digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(Digits[iDigit] ) ) ;
+ digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;
fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit);
efit[iDigit] = 0;
// so that energy deposited in each cell is distributed between new clusters proportionally
// to its contribution to efit
- Float_t * Energies = iniTower->GetEnergiesList() ;
+ Float_t * energiesList = iniTower->GetEnergiesList() ;
Float_t ratio ;
iparam = 0 ;
epar = fitparameters[iparam+2] ;
iparam += 3 ;
- AliEMCALRecPoint * RecPoint = 0 ;
+ AliEMCALRecPoint * recPoint = 0 ;
if(fNumberOfECAClusters >= fRecPoints->GetSize())
fRecPoints->Expand(2*fNumberOfECAClusters) ;
(*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;
- RecPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
+ recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
fNumberOfECAClusters++ ;
- RecPoint->SetNExMax((Int_t)nPar/3) ;
+ recPoint->SetNExMax((Int_t)nPar/3) ;
Float_t eDigit ;
for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
- digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( Digits[iDigit] ) ) ;
+ digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit);
ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
- eDigit = Energies[iDigit] * ratio ;
- RecPoint->AddDigit( *digit, eDigit ) ;
+ eDigit = energiesList[iDigit] * ratio ;
+ recPoint->AddDigit( *digit, eDigit ) ;
}
}
TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
- AliEMCALRecPoint * RecPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) ) ;
+ AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) ) ;
TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
// A bit buggy way to get an access to the geometry
// To be revised!
AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2));
- Int_t * Digits = RecPoint->GetDigitsList() ;
+ Int_t * digitsList = recPoint->GetDigitsList() ;
- Int_t nOdigits = RecPoint->GetDigitsMultiplicity() ;
+ Int_t nOdigits = recPoint->GetDigitsMultiplicity() ;
- Float_t * Energies = RecPoint->GetEnergiesList() ;
+ Float_t * energiesList = recPoint->GetEnergiesList() ;
fret = 0. ;
Int_t iparam ;
for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
- digit = dynamic_cast<AliEMCALDigit*>( digits->At( Digits[iDigit] ) );
+ digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) );
Double_t xDigit=0 ;
Double_t zDigit=0 ;
efit += x[iParam] * ShowerShape(dx,dz) ;
iParam++ ;
}
- Double_t sum = 2. * (efit - Energies[iDigit]) / Energies[iDigit] ; // Here we assume, that sigma = sqrt(E)
+ Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E)
iParam = 0 ;
while(iParam < nPar ){
Double_t xpar = x[iParam] ;
efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
}
- fret += (efit-Energies[iDigit])*(efit-Energies[iDigit])/Energies[iDigit] ;
+ fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ;
// Here we assume, that sigma = sqrt(E)
}
}
TString taskName(Version()) ;
printf("--------------- ");
- printf(taskName.Data()) ;
+ printf("%s",taskName.Data()) ;
printf(" ");
printf("Clusterizing digits: ");
printf("\n ECA Local Maximum cut = %f", fECALocMaxCut);