X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ANALYSIS%2FTenderSupplies%2FAliEMCALTenderSupply.cxx;h=5bd83bf8f57d20cf38c31db1d1af308ae6ee4c16;hb=7bf608c977247d84476857687aa73b7079bfb1bc;hp=9f704e6e42d77d30e68d23d58bc3ba84c6db7f7b;hpb=94f7c83694bce64e30bc557d5f03ba7b165fbcf7;p=u%2Fmrichter%2FAliRoot.git diff --git a/ANALYSIS/TenderSupplies/AliEMCALTenderSupply.cxx b/ANALYSIS/TenderSupplies/AliEMCALTenderSupply.cxx index 9f704e6e42d..5bd83bf8f57 100644 --- a/ANALYSIS/TenderSupplies/AliEMCALTenderSupply.cxx +++ b/ANALYSIS/TenderSupplies/AliEMCALTenderSupply.cxx @@ -13,12 +13,13 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ - /////////////////////////////////////////////////////////////////////////////// // // -// EMCAL tender, apply corrections to EMCAl clusters // -// and do track matching. // +// EMCAL tender, apply corrections to EMCAL clusters and do track matching. // +// // // Author: Deepa Thomas (Utrecht University) // +// Later mods/rewrite: Jiri Kral (University of Jyvaskyla) // +// S. Aiola, C. Loizides : Make it work for AODs // // // /////////////////////////////////////////////////////////////////////////////// @@ -29,118 +30,197 @@ #include #include #include -#include -#include -#include -#include -#include "AliOADBContainer.h" -#include "AliCDBManager.h" -#include "AliCDBStorage.h" -#include "AliCDBEntry.h" -#include "AliMagF.h" -#include "AliESDCaloCluster.h" -#include "AliEMCALTenderSupply.h" -#include "AliEMCALGeometry.h" -#include "AliEMCALRecoUtils.h" -#include "AliEMCALClusterizer.h" -#include "AliEMCALRecParam.h" -#include "AliEMCALRecPoint.h" #include "AliEMCALAfterBurnerUF.h" +#include "AliEMCALClusterizer.h" #include "AliEMCALClusterizerNxN.h" #include "AliEMCALClusterizerv1.h" #include "AliEMCALClusterizerv2.h" #include "AliEMCALDigit.h" +#include "AliEMCALGeometry.h" +#include "AliEMCALRecParam.h" +#include "AliEMCALRecParam.h" +#include "AliEMCALRecPoint.h" +#include "AliEMCALRecoUtils.h" +#include "AliEMCALTenderSupply.h" +#include "AliESDCaloCluster.h" +#include "AliMagF.h" +#include "AliOADBContainer.h" +#include "AliAODEvent.h" +#include "AliAnalysisManager.h" +#include "AliESDEvent.h" +#include "AliLog.h" +#include "AliTender.h" ClassImp(AliEMCALTenderSupply) AliEMCALTenderSupply::AliEMCALTenderSupply() : AliTenderSupply() +,fTask(0) +,fRun(0) ,fEMCALGeo(0x0) -,fEMCALGeoName("EMCAL_COMPLETEV1") +,fEMCALGeoName("") ,fEMCALRecoUtils(0) ,fConfigName("") ,fDebugLevel(0) -,fNonLinearFunc(AliEMCALRecoUtils::kNoCorrection) -,fNonLinearThreshold(30) -,fReCalibCluster(kFALSE) +,fNonLinearFunc(-1) +,fNonLinearThreshold(-1) +,fReCalibCluster(kFALSE) ,fUpdateCell(kFALSE) +,fCalibrateEnergy(kFALSE) +,fCalibrateTime(kFALSE) +,fDoNonLinearity(kFALSE) +,fBadCellRemove(kFALSE) +,fRejectExoticCells(kFALSE) +,fRejectExoticClusters(kFALSE) +,fClusterBadChannelCheck(kFALSE) ,fRecalClusPos(kFALSE) ,fFiducial(kFALSE) -,fNCellsFromEMCALBorder(1) -,fRecalDistToBadChannels(kFALSE) -,fInputTree(0) +,fNCellsFromEMCALBorder(-1) +,fRecalDistToBadChannels(kFALSE) +,fRecalShowerShape(kFALSE) +,fInputTree(0) ,fInputFile(0) ,fFilepass(0) -,fMass(0.139) -,fStep(50) +,fMass(-1) +,fStep(-1) ,fCutEtaPhiSum(kTRUE) ,fCutEtaPhiSeparate(kFALSE) -,fRcut(0.05) -,fEtacut(0.025) -,fPhicut(0.05) +,fRcut(-1) +,fEtacut(-1) +,fPhicut(-1) ,fBasePath("") ,fReClusterize(kFALSE) ,fClusterizer(0) ,fGeomMatrixSet(kFALSE) ,fLoadGeomMatrices(kFALSE) ,fRecParam(0x0) -,fRecParamSet(kFALSE) -,fOCDBpath(" ") +,fDoTrackMatch(kFALSE) +,fDoUpdateOnly(kFALSE) ,fUnfolder(0) ,fDigitsArr(0) ,fClusterArr(0) -,fMisalignSurvey(kdefault) +,fMisalignSurvey(kdefault) +,fExoticCellFraction(-1) +,fExoticCellDiffTime(-1) +,fExoticCellMinAmplitude(-1) { // Default constructor. - for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ; + + for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ; } //_____________________________________________________ AliEMCALTenderSupply::AliEMCALTenderSupply(const char *name, const AliTender *tender) : AliTenderSupply(name,tender) +,fTask(0) +,fRun(0) ,fEMCALGeo(0x0) -,fEMCALGeoName("EMCAL_COMPLETEV1") +,fEMCALGeoName("") ,fEMCALRecoUtils(0) -,fConfigName("") +,fConfigName("") ,fDebugLevel(0) -,fNonLinearFunc(AliEMCALRecoUtils::kNoCorrection) -,fNonLinearThreshold(30) -,fReCalibCluster(kFALSE) +,fNonLinearFunc(-1) +,fNonLinearThreshold(-1) +,fReCalibCluster(kFALSE) ,fUpdateCell(kFALSE) +,fCalibrateEnergy(kFALSE) +,fCalibrateTime(kFALSE) +,fDoNonLinearity(kFALSE) +,fBadCellRemove(kFALSE) +,fRejectExoticCells(kFALSE) +,fRejectExoticClusters(kFALSE) +,fClusterBadChannelCheck(kFALSE) ,fRecalClusPos(kFALSE) ,fFiducial(kFALSE) -,fNCellsFromEMCALBorder(1) -,fRecalDistToBadChannels(kFALSE) -,fInputTree(0) +,fNCellsFromEMCALBorder(-1) +,fRecalDistToBadChannels(kFALSE) +,fRecalShowerShape(kFALSE) +,fInputTree(0) ,fInputFile(0) -,fFilepass(0) -,fMass(0.139) -,fStep(50) +,fFilepass(0) +,fMass(-1) +,fStep(-1) ,fCutEtaPhiSum(kTRUE) ,fCutEtaPhiSeparate(kFALSE) -,fRcut(0.05) -,fEtacut(0.025) -,fPhicut(0.05) +,fRcut(-1) +,fEtacut(-1) +,fPhicut(-1) ,fBasePath("") ,fReClusterize(kFALSE) ,fClusterizer(0) ,fGeomMatrixSet(kFALSE) ,fLoadGeomMatrices(kFALSE) ,fRecParam(0x0) -,fRecParamSet(kFALSE) -,fOCDBpath(" ") +,fDoTrackMatch(kFALSE) +,fDoUpdateOnly(kFALSE) ,fUnfolder(0) ,fDigitsArr(0) ,fClusterArr(0) -,fMisalignSurvey(kdefault) +,fMisalignSurvey(kdefault) +,fExoticCellFraction(-1) +,fExoticCellDiffTime(-1) +,fExoticCellMinAmplitude(-1) { // Named constructor - for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ; + for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ; +} + +//_____________________________________________________ +AliEMCALTenderSupply::AliEMCALTenderSupply(const char *name, AliAnalysisTaskSE *task) : +AliTenderSupply(name) +,fTask(task) +,fRun(0) +,fEMCALGeo(0x0) +,fEMCALGeoName("") +,fEMCALRecoUtils(0) +,fConfigName("") +,fDebugLevel(0) +,fNonLinearFunc(-1) +,fNonLinearThreshold(-1) +,fReCalibCluster(kFALSE) +,fUpdateCell(kFALSE) +,fCalibrateEnergy(kFALSE) +,fCalibrateTime(kFALSE) +,fDoNonLinearity(kFALSE) +,fBadCellRemove(kFALSE) +,fRejectExoticCells(kFALSE) +,fRejectExoticClusters(kFALSE) +,fClusterBadChannelCheck(kFALSE) +,fRecalClusPos(kFALSE) +,fFiducial(kFALSE) +,fNCellsFromEMCALBorder(-1) +,fRecalDistToBadChannels(kFALSE) +,fRecalShowerShape(kFALSE) +,fInputTree(0) +,fInputFile(0) +,fFilepass(0) +,fMass(-1) +,fStep(-1) +,fCutEtaPhiSum(kTRUE) +,fCutEtaPhiSeparate(kFALSE) +,fRcut(-1) +,fEtacut(-1) +,fPhicut(-1) +,fBasePath("") +,fReClusterize(kFALSE) +,fClusterizer(0) +,fGeomMatrixSet(kFALSE) +,fLoadGeomMatrices(kFALSE) +,fRecParam(0x0) +,fDoTrackMatch(kFALSE) +,fDoUpdateOnly(kFALSE) +,fUnfolder(0) +,fDigitsArr(0) +,fClusterArr(0) +,fMisalignSurvey(kdefault) +,fExoticCellFraction(-1) +,fExoticCellDiffTime(-1) +,fExoticCellMinAmplitude(-1) +{ + // Named constructor. - fEMCALRecoUtils = new AliEMCALRecoUtils(); - fRecParam = new AliEMCALRecParam; - fDigitsArr = new TClonesArray("AliEMCALDigit",1000); + for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ; } //_____________________________________________________ @@ -150,23 +230,49 @@ AliEMCALTenderSupply::~AliEMCALTenderSupply() delete fEMCALRecoUtils; delete fRecParam; - delete fClusterizer; delete fUnfolder; - if (fDigitsArr){ + if (!fClusterizer) { fDigitsArr->Clear("C"); delete fDigitsArr; + } else { + delete fClusterizer; + fDigitsArr = 0; } } +//_____________________________________________________ +void AliEMCALTenderSupply::SetDefaults() +{ + // Set default settings. + + SwitchOnBadCellRemove(); + SwitchOnExoticCellRemove(); + SwitchOnCalibrateEnergy(); + SwitchOnCalibrateTime(); + SwitchOnUpdateCell(); + SwitchOnReclustering(); + SwitchOnClusterBadChannelCheck(); + SwitchOnClusterExoticChannelCheck(); + SwitchOnTrackMatch(); +} + +//_____________________________________________________ +Bool_t AliEMCALTenderSupply::RunChanged() const +{ + // Get run number. + + return (fTender && fTender->RunChanged()) || (fTask && fRun != fTask->InputEvent()->GetRunNumber()); +} + //_____________________________________________________ void AliEMCALTenderSupply::Init() { // Initialise EMCAL tender. if (fDebugLevel>0) - AliInfo("Init EMCAL Tender supply"); + AliWarning("Init EMCAL Tender supply"); - if(fConfigName.Length()>0 && gROOT->LoadMacro(fConfigName) >=0) { + if (fConfigName.Length()>0 && gROOT->LoadMacro(fConfigName) >=0) { AliDebug(1, Form("Loading settings from macro %s", fConfigName.Data())); AliEMCALTenderSupply *tender = (AliEMCALTenderSupply*)gInterpreter->ProcessLine("ConfigEMCALTenderSupply()"); fDebugLevel = tender->fDebugLevel; @@ -176,132 +282,373 @@ void AliEMCALTenderSupply::Init() fNonLinearFunc = tender->fNonLinearFunc; fNonLinearThreshold = tender->fNonLinearThreshold; fReCalibCluster = tender->fReCalibCluster; + fUpdateCell = tender->fUpdateCell; fRecalClusPos = tender->fRecalClusPos; - fFiducial = tender->fFiducial; + fCalibrateEnergy = tender->fCalibrateEnergy; + fCalibrateTime = tender->fCalibrateTime; + fFiducial = tender->fFiducial; fNCellsFromEMCALBorder = tender->fNCellsFromEMCALBorder; fRecalDistToBadChannels = tender->fRecalDistToBadChannels; + fRecalShowerShape = tender->fRecalShowerShape; + fClusterBadChannelCheck = tender->fClusterBadChannelCheck; + fBadCellRemove = tender->fBadCellRemove; + fRejectExoticCells = tender->fRejectExoticCells; + fRejectExoticClusters = tender->fRejectExoticClusters; fMass = tender->fMass; fStep = tender->fStep; + fCutEtaPhiSum = tender->fCutEtaPhiSum; + fCutEtaPhiSeparate = tender->fCutEtaPhiSeparate; fRcut = tender->fRcut; fEtacut = tender->fEtacut; fPhicut = tender->fPhicut; fReClusterize = tender->fReClusterize; fLoadGeomMatrices = tender->fLoadGeomMatrices; fRecParam = tender->fRecParam; - fOCDBpath = tender->fOCDBpath; - for(Int_t i = 0; i < 10; i++) + fDoNonLinearity = tender->fDoNonLinearity; + fDoTrackMatch = tender->fDoTrackMatch; + fDoUpdateOnly = tender->fDoUpdateOnly; + fMisalignSurvey = tender->fMisalignSurvey; + fExoticCellFraction = tender->fExoticCellFraction; + fExoticCellDiffTime = tender->fExoticCellDiffTime; + fExoticCellMinAmplitude = tender->fExoticCellMinAmplitude; + + for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = tender->fEMCALMatrix[i] ; } - // Init goemetry - fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ; - - // Initialising Non linearity parameters - fEMCALRecoUtils->SetNonLinearityThreshold(fNonLinearThreshold); - fEMCALRecoUtils->SetNonLinearityFunction(fNonLinearFunc); + if (fDebugLevel>0){ + AliInfo( "Emcal Tender settings: ======================================" ); + AliInfo( "------------ Switches --------------------------" ); + AliInfo( Form( "BadCellRemove : %d", fBadCellRemove )); + AliInfo( Form( "ExoticCellRemove : %d", fRejectExoticCells )); + AliInfo( Form( "CalibrateEnergy : %d", fCalibrateEnergy )); + AliInfo( Form( "CalibrateTime : %d", fCalibrateTime )); + AliInfo( Form( "UpdateCell : %d", fUpdateCell )); + AliInfo( Form( "DoUpdateOnly : %d", fDoUpdateOnly )); + AliInfo( Form( "Reclustering : %d", fReClusterize )); + AliInfo( Form( "ClusterBadChannelCheck : %d", fClusterBadChannelCheck )); + AliInfo( Form( "ClusterExoticChannelCheck : %d", fRejectExoticClusters )); + AliInfo( Form( "CellFiducialRegion : %d", fFiducial )); + AliInfo( Form( "ReCalibrateCluster : %d", fReCalibCluster )); + AliInfo( Form( "RecalculateClusPos : %d", fRecalClusPos )); + AliInfo( Form( "RecalShowerShape : %d", fRecalShowerShape )); + AliInfo( Form( "NonLinearityCorrection : %d", fDoNonLinearity )); + AliInfo( Form( "RecalDistBadChannel : %d", fRecalDistToBadChannels )); + AliInfo( Form( "TrackMatch : %d", fDoTrackMatch )); + AliInfo( "------------ Variables -------------------------" ); + AliInfo( Form( "DebugLevel : %d", fDebugLevel )); + AliInfo( Form( "BasePath : %s", fBasePath.Data() )); + AliInfo( Form( "ConfigFileName : %s", fConfigName.Data() )); + AliInfo( Form( "EMCALGeometryName : %s", fEMCALGeoName.Data() )); + AliInfo( Form( "NonLinearityFunction : %d", fNonLinearFunc )); + AliInfo( Form( "NonLinearityThreshold : %d", fNonLinearThreshold )); + AliInfo( Form( "MisalignmentMatrixSurvey : %d", fMisalignSurvey )); + AliInfo( Form( "NumberOfCellsFromEMCALBorder : %d", fNCellsFromEMCALBorder )); + AliInfo( Form( "RCut : %f", fRcut )); + AliInfo( Form( "Mass : %f", fMass )); + AliInfo( Form( "Step : %f", fStep )); + AliInfo( Form( "EtaCut : %f", fEtacut )); + AliInfo( Form( "PhiCut : %f", fPhicut )); + AliInfo( Form( "ExoticCellFraction : %f", fExoticCellFraction )); + AliInfo( Form( "ExoticCellDiffTime : %f", fExoticCellDiffTime )); + AliInfo( Form( "ExoticCellMinAmplitude : %f", fExoticCellMinAmplitude )); + AliInfo( "=============================================================" ); + } + + // init reco utils + fEMCALRecoUtils = new AliEMCALRecoUtils; + + // init geometry if requested + if (fEMCALGeoName.Length()>0) + fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ; + + // digits array + fDigitsArr = new TClonesArray("AliEMCALDigit",1000); + + // initialising non-linearity parameters + if( fNonLinearThreshold != -1 ) + fEMCALRecoUtils->SetNonLinearityThreshold(fNonLinearThreshold); + if( fNonLinearFunc != -1 ) + fEMCALRecoUtils->SetNonLinearityFunction(fNonLinearFunc); + + // missalignment function + fEMCALRecoUtils->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerGlobal); + + // fiducial cut + // do not do the eta0 fiducial cut + if( fNCellsFromEMCALBorder != -1 ) + fEMCALRecoUtils->SetNumberOfCellsFromEMCALBorder(fNCellsFromEMCALBorder); + fEMCALRecoUtils->SwitchOnNoFiducialBorderInEMCALEta0(); + + // exotic cell rejection + if( fExoticCellFraction != -1 ) + fEMCALRecoUtils->SetExoticCellFractionCut( fExoticCellFraction ); + if( fExoticCellDiffTime != -1 ) + fEMCALRecoUtils->SetExoticCellDiffTimeCut( fExoticCellDiffTime ); + if( fExoticCellMinAmplitude != -1 ) + fEMCALRecoUtils->SetExoticCellMinAmplitudeCut( fExoticCellMinAmplitude ); + + // setting track matching parameters ... mass, step size and residual cut + if( fMass != -1 ) + fEMCALRecoUtils->SetMass(fMass); + if( fStep != -1 ) + fEMCALRecoUtils->SetStep(fStep); - // Setting mass, step size and residual cut - fEMCALRecoUtils->SetMass(fMass); - fEMCALRecoUtils->SetStep(fStep); + // spatial cut based on separate eta/phi or common processing if(fCutEtaPhiSum){ fEMCALRecoUtils->SwitchOnCutEtaPhiSum(); - fEMCALRecoUtils->SetCutR(fRcut); + if( fRcut != -1 ) + fEMCALRecoUtils->SetCutR(fRcut); } else if (fCutEtaPhiSeparate) { fEMCALRecoUtils->SwitchOnCutEtaPhiSeparate(); - fEMCALRecoUtils->SetCutEta(fEtacut); - fEMCALRecoUtils->SetCutPhi(fPhicut); + if( fEtacut != -1 ) + fEMCALRecoUtils->SetCutEta(fEtacut); + if( fPhicut != -1 ) + fEMCALRecoUtils->SetCutPhi(fPhicut); } } +//_____________________________________________________ +AliVEvent* AliEMCALTenderSupply::GetEvent() +{ + // Return the event pointer. + + if (fTender) { + return fTender->GetEvent(); + } + else if (fTask) { + return fTask->InputEvent(); + } + + return 0; +} + //_____________________________________________________ void AliEMCALTenderSupply::ProcessEvent() { // Event loop. - AliESDEvent *event = fTender->GetEvent(); + AliVEvent *event = GetEvent(); + if (!event) { - AliError("ESD event ptr = 0, returning"); + AliError("Event ptr = 0, returning"); return; } - // Initialising parameters once per run number - if(fTender->RunChanged()){ + // Initialising parameters once per run number + if (RunChanged()) { + + AliWarning( "Run changed, initializing parameters" ); + fRun = event->GetRunNumber(); + + // init geometry if not already done + if (fEMCALGeoName.Length()==0) { + fEMCALGeoName = "EMCAL_FIRSTYEARV1"; + if (fRun>139517) { + fEMCALGeoName = "EMCAL_COMPLETEV1"; + } + if (fRun>170593) { + fEMCALGeoName = "EMCAL_COMPLETE12SMV1"; + } + fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName); + if (!fEMCALGeo) { + AliFatal(Form("Can not create geometry: %s", fEMCALGeoName.Data())); + return; + } + } + + // get pass GetPass(); - // Init bad channels - Int_t fInitBC=InitBadChannels(); - - if (fInitBC==0) - { - AliError("InitBadChannels returned false, returning"); - return; + + // define what recalib parameters are needed for various switches + // this is based on implementation in AliEMCALRecoUtils + Bool_t needRecoParam = fReClusterize; + Bool_t needBadChannels = fBadCellRemove | fClusterBadChannelCheck | fRecalDistToBadChannels | fReClusterize; + Bool_t needRecalib = fCalibrateEnergy | fReClusterize; + Bool_t needTimecalib = fCalibrateTime | fReClusterize; + Bool_t needMisalign = fRecalClusPos | fReClusterize; + Bool_t needClusterizer = fReClusterize; + + // initiate reco params with some defaults + // will not overwrite, if those have been provided by user + if( needRecoParam ){ + Int_t initRC = InitRecParam(); + + if( initRC == 0 ) + AliInfo("Defaults reco params loaded."); + if( initRC > 1 ) + AliWarning("User defined reco params."); } - if(fInitBC>1) - { - AliInfo(Form("No external hot channel set: %d - %s", event->GetRunNumber(), fFilepass.Data())); + + // init bad channels + if( needBadChannels ){ + Int_t fInitBC = InitBadChannels(); + if (fInitBC==0) + AliError("InitBadChannels returned false, returning"); + if (fInitBC==1) + AliWarning("InitBadChannels OK"); + if (fInitBC>1) + AliWarning(Form("No external hot channel set: %d - %s", event->GetRunNumber(), fFilepass.Data())); } - - if (fReCalibCluster || fUpdateCell) { - Int_t fInitRecalib=InitRecalib(); + + // init recalibration factors + if( needRecalib ) + { + Int_t fInitRecalib = InitRecalib(); if (fInitRecalib==0) - { AliError("InitRecalib returned false, returning"); - return; - } - if(fInitRecalib >1) - { - AliInfo(Form("No recalibration available: %d - %s", event->GetRunNumber(), fFilepass.Data())); - fUpdateCell = kFALSE; - fReCalibCluster = kFALSE; - } + if (fInitRecalib==1) + AliWarning("InitRecalib OK"); + if (fInitRecalib>1) + AliWarning(Form("No recalibration available: %d - %s", event->GetRunNumber(), fFilepass.Data())); + + Int_t fInitRunDepRecalib = InitRunDepRecalib(); + if (fInitRunDepRecalib==0) + AliError("InitrunDepRecalib returned false, returning"); + if (fInitRunDepRecalib==1) + AliWarning("InitRecalib OK"); + if (fInitRunDepRecalib>1) + AliWarning(Form("No Temperature recalibration available: %d - %s", event->GetRunNumber(), fFilepass.Data())); + + fReCalibCluster = kFALSE; + } - if (fRecalClusPos || fReClusterize || fUpdateCell) { - if (!InitMisalignMatrix()) { + + + // init time calibration + if( needTimecalib ){ + Int_t initTC = InitTimeCalibration(); + if ( !initTC ) + AliError("InitTimeCalibration returned false, returning"); + if (initTC==1) + AliWarning("InitTimeCalib OK"); + if( initTC > 1 ) + AliWarning(Form("No external time calibration set: %d - %s", event->GetRunNumber(), fFilepass.Data())); + } + + // init misalignment matrix + if( needMisalign ) { + if (!InitMisalignMatrix()) AliError("InitMisalignmentMatrix returned false, returning"); - return; - } + else + AliWarning("InitMisalignMatrix OK"); } - if (fReClusterize || fUpdateCell) { - if (!InitClusterization()) { + // init clusterizer + if( needClusterizer ) { + if (!InitClusterization()) AliError("InitClusterization returned false, returning"); - return; - } + else + AliWarning("InitClusterization OK"); } if(fDebugLevel>1) fEMCALRecoUtils->Print(""); + } + + // disable implied switches ------------------------------------------------- + // AliEMCALRecoUtils or clusterizer functions alredy include some + // recalibration so based on those implied callibration te switches + // here are disabled to avoid duplication + // clusterizer does cluster energy recalibration, position recomputation + // and shower shape + if( fReClusterize ){ + fReCalibCluster = kFALSE; + fRecalClusPos = kFALSE; + fRecalShowerShape = kFALSE; } + // CONFIGURE THE RECO UTILS ------------------------------------------------- + // configure the reco utils + // this option does energy recalibration + if( fCalibrateEnergy ) + fEMCALRecoUtils->SwitchOnRecalibration(); + else + fEMCALRecoUtils->SwitchOffRecalibration(); + + // allows time calibration + if( fCalibrateTime ) + fEMCALRecoUtils->SwitchOnTimeRecalibration(); + else + fEMCALRecoUtils->SwitchOffTimeRecalibration(); + + // allows to zero bad cells + if( fBadCellRemove ) + fEMCALRecoUtils->SwitchOnBadChannelsRemoval(); + else + fEMCALRecoUtils->SwitchOffBadChannelsRemoval(); + + // distance to bad channel recalibration + if( fRecalDistToBadChannels ) + fEMCALRecoUtils->SwitchOnDistToBadChannelRecalculation(); + else + fEMCALRecoUtils->SwitchOffDistToBadChannelRecalculation(); + + // exclude exotic cells + if( fRejectExoticCells ) + fEMCALRecoUtils->SwitchOnRejectExoticCell(); + else + fEMCALRecoUtils->SwitchOffRejectExoticCell(); + + // exclude clusters with exotic cells + if( fRejectExoticClusters ) + fEMCALRecoUtils->SwitchOnRejectExoticCluster(); + else + fEMCALRecoUtils->SwitchOffRejectExoticCluster(); + + // TODO: not implemented switches + // SwitchOnClusterEnergySmearing + // SwitchOnRunDepCorrection + + // START PROCESSING --------------------------------------------------------- // Test if cells present - AliESDCaloCells *cells= event->GetEMCALCells(); - if (cells->GetNumberOfCells()<=0) { + AliVCaloCells *cells= event->GetEMCALCells(); + if (cells->GetNumberOfCells()<=0) + { if(fDebugLevel>1) AliWarning(Form("Number of EMCAL cells = %d, returning", cells->GetNumberOfCells())); return; } - if(fDebugLevel>2) + if (fDebugLevel>2) AliInfo(Form("Re-calibrate cluster %d\n",fReCalibCluster)); + + // mark the cells not recalibrated in case of selected + // time, energy recalibration or bad channel removal + if( fCalibrateEnergy || fCalibrateTime || fBadCellRemove ) + fEMCALRecoUtils->ResetCellsCalibrated(); - // Recalibrate cells - if(fUpdateCell) + // CELL RECALIBRATION ------------------------------------------------------- + // cell objects will be updated + // the cell calibrations are also processed locally any time those are needed + // in case that the cell objects are not to be updated here for later use + if( fUpdateCell ) { + // do the update + // includes exotic cell check in the UpdateCells function - is not provided + // by the reco utils UpdateCells(); - fReClusterize = kTRUE; - } + + // switch off recalibrations so those are not done multiple times + // this is just for safety, the recalibrated flag of cell object + // should not allow for farther processing anyways + fEMCALRecoUtils->SwitchOffRecalibration(); + fEMCALRecoUtils->SwitchOffTimeRecalibration(); - // Reclusterize - if(fReClusterize) + if (fDoUpdateOnly) + return; + } + + // RECLUSTERIZATION --------------------------------------------------------- + if (fReClusterize) { - fRecalClusPos = kFALSE; // Done in reclusterization, do nothing - fReCalibCluster= kFALSE; // Done in reclusterization, do nothing FillDigitsArray(); Clusterize(); UpdateClusters(); } - + // Store good clusters TClonesArray *clusArr = dynamic_cast(event->FindListObject("caloClusters")); if (!clusArr) @@ -310,82 +657,153 @@ void AliEMCALTenderSupply::ProcessEvent() AliWarning(Form("No cluster array, number of cells in event = %d, returning", cells->GetNumberOfCells())); return; } - + + // PROCESS SINGLE CLUSTER RECALIBRATION ------------------------------------- + // now go through clusters one by one and process separate correction + // as those were defined or not Int_t nclusters = clusArr->GetEntriesFast(); - - for (Int_t icluster=0; icluster < nclusters; ++icluster) { + for (Int_t icluster=0; icluster < nclusters; ++icluster) + { AliVCluster *clust = static_cast(clusArr->At(icluster)); if (!clust) continue; if (!clust->IsEMCAL()) continue; - if (fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo, clust->GetCellsAbsId(), clust->GetNCells())) { - delete clusArr->RemoveAt(icluster); - continue; //todo is it really needed to remove it? Or should we flag it? + + // REMOVE CLUSTERS WITH BAD CELLS ----------------------------- + if( fClusterBadChannelCheck ) + { + // careful, the the ClusterContainsBadChannel is dependent on + // SwitchOnBadChannelsRemoval, switching it ON automatically + // and returning to original value after processing + Bool_t badRemoval = fEMCALRecoUtils->IsBadChannelsRemovalSwitchedOn(); + fEMCALRecoUtils->SwitchOnBadChannelsRemoval(); + + Bool_t badResult = fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo, clust->GetCellsAbsId(), clust->GetNCells()); + + // switch the bad channels removal back + if( ! badRemoval ) + fEMCALRecoUtils->SwitchOffBadChannelsRemoval(); + + if( badResult ) + { + delete clusArr->RemoveAt(icluster); + continue; //TODO is it really needed to remove it? Or should we flag it? + } } - if (fFiducial){ + // REMOVE EXOTIC CLUSTERS ------------------------------------- + // does process local cell recalibration energy and time without replacing + // the global cell values, in case of no cell recalib done yet + if( fRejectExoticClusters ) + { + // careful, the IsExoticCluster is dependent on + // SwitchOnRejectExoticCell, switching it ON automatically + // and returning to original value after processing + Bool_t exRemoval = fEMCALRecoUtils->IsRejectExoticCell(); + fEMCALRecoUtils->SwitchOnRejectExoticCell(); + + // get bunch crossing + Int_t bunchCrossNo = event->GetBunchCrossNumber(); + + Bool_t exResult = fEMCALRecoUtils->IsExoticCluster(clust, cells, bunchCrossNo ); + + // switch the exotic channels removal back + if( ! exRemoval ) + fEMCALRecoUtils->SwitchOffRejectExoticCell(); + + if( exResult ) + { + delete clusArr->RemoveAt(icluster); + continue; //TODO is it really needed to remove it? Or should we flag it? + } + } + + // FIDUCIAL CUT ----------------------------------------------- + if (fFiducial) + { + // depends on SetNumberOfCellsFromEMCALBorder + // SwitchOnNoFiducialBorderInEMCALEta0 if (!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, clust, cells)){ delete clusArr->RemoveAt(icluster); - continue; // todo it would be nice to store the distance + continue; //TODO it would be nice to store the distance } } - if(fRecalDistToBadChannels) - fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, cells, clust); - if(fReCalibCluster) + // CLUSTER ENERGY --------------------------------------------- + // does process local cell recalibration energy and time without replacing + // the global cell values, in case of no cell recalib done yet + if( fReCalibCluster ) fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, clust, cells); - if(fRecalClusPos) + + // CLUSTER POSITION ------------------------------------------- + // does process local cell energy recalibration, if enabled and cells + // not calibrated yet + if( fRecalClusPos ) fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, cells, clust); - Float_t correctedEnergy = fEMCALRecoUtils->CorrectClusterEnergyLinearity(clust); - clust->SetE(correctedEnergy); - - } + // SHOWER SHAPE ----------------------------------------------- + if( fRecalShowerShape ) + fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, cells, clust); + // NONLINEARITY ----------------------------------------------- + if( fDoNonLinearity ) + { + Float_t correctedEnergy = fEMCALRecoUtils->CorrectClusterEnergyLinearity(clust); + clust->SetE(correctedEnergy); + } + + // DISTANCE TO BAD CHANNELS ----------------------------------- + if( fRecalDistToBadChannels ) + fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, cells, clust); + } clusArr->Compress(); - - // Track matching + + if (!fDoTrackMatch) + return; + + // TRACK MATCHING ----------------------------------------------------------- if (!TGeoGlobalMagField::Instance()->GetField()) { - const AliESDRun *erun = event->GetESDRun(); - AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(), - erun->GetCurrentDip(), - AliMagF::kConvLHC, - kFALSE, - erun->GetBeamEnergy(), - erun->GetBeamType()); - TGeoGlobalMagField::Instance()->SetField(field); + AliESDEvent *esd = dynamic_cast(event); + if (esd) + esd->InitMagneticField(); + else { + AliAODEvent *aod = dynamic_cast(event); + Double_t curSol = 30000*aod->GetMagneticField()/5.00668; + Double_t curDip = 6000 *aod->GetMuonMagFieldScale(); + AliMagF *field = AliMagF::CreateFieldMap(curSol,curDip); + TGeoGlobalMagField::Instance()->SetField(field); + } } - + fEMCALRecoUtils->FindMatches(event,0x0,fEMCALGeo); - fEMCALRecoUtils->SetClusterMatchedToTrack(event); fEMCALRecoUtils->SetTracksMatchedToCluster(event); - } //_____________________________________________________ Bool_t AliEMCALTenderSupply::InitMisalignMatrix() { // Initialising misalignment matrices + + AliVEvent *event = GetEvent(); - AliESDEvent *event = fTender->GetEvent(); if (!event) return kFALSE; - if (fGeomMatrixSet) { - AliInfo("Misalignment matrix already set"); + if (fGeomMatrixSet) + { + AliInfo("Misalignment matrix already set"); return kTRUE; } if (fDebugLevel>0) - AliInfo("Initialising misalignment matrix"); - - fEMCALRecoUtils->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerGlobal); + AliInfo("Initialising misalignment matrix"); if (fLoadGeomMatrices) { - for(Int_t mod=0; mod < fEMCALGeo->GetNumberOfSuperModules(); ++mod) { + for(Int_t mod=0; mod < fEMCALGeo->GetNumberOfSuperModules(); ++mod) + { if (fEMCALMatrix[mod]){ if(fDebugLevel > 2) fEMCALMatrix[mod]->Print(); @@ -399,44 +817,56 @@ Bool_t AliEMCALTenderSupply::InitMisalignMatrix() Int_t runGM = event->GetRunNumber(); TObjArray *mobj = 0; - if(fMisalignSurvey == kdefault){ //take default alignment corresponding to run no + if(fMisalignSurvey == kdefault) + { //take default alignment corresponding to run no AliOADBContainer emcalgeoCont(Form("emcal")); emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo")); mobj=(TObjArray*)emcalgeoCont.GetObject(runGM,"EmcalMatrices"); - } - - if(fMisalignSurvey == kSurveybyS){ //take alignment at sector level - if (runGM <= 140000) { //2010 data - AliOADBContainer emcalgeoCont(Form("emcal2010")); - emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo")); - mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey10"); - - } else if (runGM>140000) { // 2011 LHC11a pass1 data - AliOADBContainer emcalgeoCont(Form("emcal2011")); - emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo")); - mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey11byS"); } - } + + if(fMisalignSurvey == kSurveybyS) + { //take alignment at sector level + if (runGM <= 140000) { //2010 data + AliOADBContainer emcalgeoCont(Form("emcal2010")); + emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo")); + mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey10"); + } + else if (runGM>140000) + { // 2011 LHC11a pass1 data + AliOADBContainer emcalgeoCont(Form("emcal2011")); + emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo")); + mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey11byS"); + } + } - if(fMisalignSurvey == kSurveybyM){ //take alignment at module level - if (runGM <= 140000) { //2010 data - AliOADBContainer emcalgeoCont(Form("emcal2010")); - emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo")); - mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey10"); - - } else if (runGM>140000) { // 2011 LHC11a pass1 data - - AliOADBContainer emcalgeoCont(Form("emcal2011")); - emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo")); - mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey11byM"); + if(fMisalignSurvey == kSurveybyM) + { //take alignment at module level + if (runGM <= 140000) { //2010 data + AliOADBContainer emcalgeoCont(Form("emcal2010")); + emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo")); + mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey10"); + } + else if (runGM>140000) + { // 2011 LHC11a pass1 data + AliOADBContainer emcalgeoCont(Form("emcal2011")); + emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo")); + mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey11byM"); + } } - } - for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ - fEMCALMatrix[mod] = (TGeoHMatrix*) mobj->At(mod); - fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod); - fEMCALMatrix[mod]->Print(); - } + if(!mobj) + { + AliFatal("Geometry matrix array not found"); + return kFALSE; + } + + for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++) + { + fEMCALMatrix[mod] = (TGeoHMatrix*) mobj->At(mod); + fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod); + fEMCALMatrix[mod]->Print(); + } + return kTRUE; } @@ -444,143 +874,162 @@ Bool_t AliEMCALTenderSupply::InitMisalignMatrix() Int_t AliEMCALTenderSupply::InitBadChannels() { // Initialising bad channel maps - AliESDEvent *event = fTender->GetEvent(); + + AliVEvent *event = GetEvent(); + if (!event) return 0; if (fDebugLevel>0) AliInfo("Initialising Bad channel map"); - if (fFiducial){ - fEMCALRecoUtils->SetNumberOfCellsFromEMCALBorder(fNCellsFromEMCALBorder); - fEMCALRecoUtils->SwitchOnNoFiducialBorderInEMCALEta0(); - } - - fEMCALRecoUtils->SwitchOnBadChannelsRemoval(); - if (fRecalDistToBadChannels) - fEMCALRecoUtils->SwitchOnDistToBadChannelRecalculation(); + // init default maps first + if( !fEMCALRecoUtils->GetEMCALBadChannelStatusMapArray() ) + fEMCALRecoUtils->InitEMCALBadChannelStatusMap() ; Int_t runBC = event->GetRunNumber(); - AliOADBContainer *contBC=new AliOADBContainer(""); - if(fBasePath!=""){ //if fBasePath specified in the ->SetBasePath() + AliOADBContainer *contBC = new AliOADBContainer(""); + if (fBasePath!="") + { //if fBasePath specified in the ->SetBasePath() if (fDebugLevel>0) AliInfo(Form("Loading Bad Channels OADB from given path %s",fBasePath.Data())); + TFile *fbad=new TFile(Form("%s/EMCALBadChannels.root",fBasePath.Data()),"read"); - if(fbad->IsZombie()){ - AliFatal(Form("EMCALBadChannels.root was not found in the path provided: %s, aborting",fBasePath.Data())); + if (!fbad || fbad->IsZombie()) + { + AliFatal(Form("EMCALBadChannels.root was not found in the path provided: %s",fBasePath.Data())); return 0; } - if(fbad) delete fbad; + + if (fbad) delete fbad; + contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fBasePath.Data()),"AliEMCALBadChannels"); + } + else + { // Else choose the one in the $ALICE_ROOT directory + if (fDebugLevel>0) AliInfo("Loading Bad Channels OADB from $ALICE_ROOT/OADB/EMCAL"); + + TFile *fbad=new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root","read"); + if (!fbad || fbad->IsZombie()) + { + AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root was not found"); + return 0; + } + + if (fbad) delete fbad; + + contBC->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root","AliEMCALBadChannels"); } - else { // Else choose the one in the $ALICE_ROOT directory - if (fDebugLevel>0) AliInfo("Loading Bad Channels OADB from $ALICE_ROOT/OADB/EMCAL"); - TFile *fbad=new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root","read"); - if(fbad->IsZombie()){ - AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root was not found, aborting"); - return 0; - } - if(fbad) delete fbad; - contBC->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root","AliEMCALBadChannels"); - } TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runBC); - if(!arrayBC){ + if (!arrayBC) + { AliError(Form("No external hot channel set for run number: %d", runBC)); return 2; } - - TObjArray *arrayBCpass=(TObjArray*)arrayBC->FindObject(fFilepass); // Here, it looks for a specific pass - if(!arrayBCpass){ - AliError(Form("No external hot channel set for: %d -%s", runBC,fFilepass.Data())); - return 2; - } - - if(fDebugLevel>0) arrayBCpass->Print(); Int_t sms = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules(); - for (Int_t i=0; iGetEMCALChannelStatusMap(i); if (h) delete h; - h=(TH2I*)arrayBCpass->FindObject(Form("EMCALBadChannelMap_Mod%d",i)); + h=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i)); - if (!h) { + if (!h) + { AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i)); continue; } h->SetDirectory(0); fEMCALRecoUtils->SetEMCALChannelStatusMap(i,h); } - return 1; } - //_____________________________________________________ Int_t AliEMCALTenderSupply::InitRecalib() { // Initialising recalibration factors. - AliESDEvent *event = fTender->GetEvent(); + AliVEvent *event = GetEvent(); + if (!event) return 0; if (fDebugLevel>0) AliInfo("Initialising recalibration factors"); - fEMCALRecoUtils->SwitchOnRecalibration(); + // init default maps first + if( !fEMCALRecoUtils->GetEMCALRecalibrationFactorsArray() ) + fEMCALRecoUtils->InitEMCALRecalibrationFactors() ; + Int_t runRC = event->GetRunNumber(); AliOADBContainer *contRF=new AliOADBContainer(""); - if(fBasePath!="") { //if fBasePath specified in the ->SetBasePath() - if (fDebugLevel>0) AliInfo(Form("Loading Recalib OADB from given path %s",fBasePath.Data())); + if (fBasePath!="") + { //if fBasePath specified in the ->SetBasePath() + if (fDebugLevel>0) AliInfo(Form("Loading Recalib OADB from given path %s",fBasePath.Data())); + TFile *fRecalib= new TFile(Form("%s/EMCALRecalib.root",fBasePath.Data()),"read"); - if (fRecalib->IsZombie()) { - AliFatal(Form("EMCALRecalib.root not found in %s, aborting",fBasePath.Data())); + if (!fRecalib || fRecalib->IsZombie()) + { + AliFatal(Form("EMCALRecalib.root not found in %s",fBasePath.Data())); return 0; } - if(fRecalib) delete fRecalib; + + if (fRecalib) delete fRecalib; + contRF->InitFromFile(Form("%s/EMCALRecalib.root",fBasePath.Data()),"AliEMCALRecalib"); } - else{ // Else choose the one in the $ALICE_ROOT directory - if (fDebugLevel>0) AliInfo("Loading Recalib OADB from $ALICE_ROOT/OADB/EMCAL"); - TFile *fRecalib= new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root","read"); - if (fRecalib->IsZombie()) { - AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root was not found, aborting"); - return 0; - } - if(fRecalib) delete fRecalib; - contRF->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root","AliEMCALRecalib"); + else + { // Else choose the one in the $ALICE_ROOT directory + if (fDebugLevel>0) AliInfo("Loading Recalib OADB from $ALICE_ROOT/OADB/EMCAL"); + + TFile *fRecalib= new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root","read"); + if (!fRecalib || fRecalib->IsZombie()) + { + AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root was not found"); + return 0; } + + if (fRecalib) delete fRecalib; + + contRF->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root","AliEMCALRecalib"); + } - TObjArray *recal=(TObjArray*)contRF->GetObject(runRC); //GetObject(int runnumber) - if(!recal){ + TObjArray *recal=(TObjArray*)contRF->GetObject(runRC); + if (!recal) + { AliError(Form("No Objects for run: %d",runRC)); return 2; } TObjArray *recalpass=(TObjArray*)recal->FindObject(fFilepass); - if(!recalpass){ + if (!recalpass) + { AliError(Form("No Objects for run: %d - %s",runRC,fFilepass.Data())); return 2; } TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib"); - if(!recalib){ + if (!recalib) + { AliError(Form("No Recalib histos found for %d - %s",runRC,fFilepass.Data())); return 2; } - if(fDebugLevel>0) recalib->Print(); + if (fDebugLevel>0) recalib->Print(); Int_t sms = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules(); - for (Int_t i=0; iGetEMCALChannelRecalibrationFactors(i); if (h) delete h; h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i)); - if (!h) { + if (!h) + { AliError(Form("Could not load EMCALRecalFactors_SM%d",i)); continue; } @@ -590,6 +1039,175 @@ Int_t AliEMCALTenderSupply::InitRecalib() return 1; } +//_____________________________________________________ +Int_t AliEMCALTenderSupply::InitRunDepRecalib() +{ + // Initialising recalibration factors. + + AliVEvent *event = GetEvent(); + + if (!event) + return 0; + + if (fDebugLevel>0) + AliInfo("Initialising recalibration factors"); + + // init default maps first + if( !fEMCALRecoUtils->GetEMCALRecalibrationFactorsArray() ) + fEMCALRecoUtils->InitEMCALRecalibrationFactors() ; + + Int_t runRC = event->GetRunNumber(); + + AliOADBContainer *contRF=new AliOADBContainer(""); + if (fBasePath!="") + { //if fBasePath specified in the ->SetBasePath() + if (fDebugLevel>0) AliInfo(Form("Loading Recalib OADB from given path %s",fBasePath.Data())); + + TFile *fRunDepRecalib= new TFile(Form("%s/EMCALTemperatureCorrCalib.root",fBasePath.Data()),"read"); + if (!fRunDepRecalib || fRunDepRecalib->IsZombie()) + { + AliFatal(Form("EMCALTemperatureCorrCalib.root not found in %s",fBasePath.Data())); + return 0; + } + + if (fRunDepRecalib) delete fRunDepRecalib; + + contRF->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fBasePath.Data()),"AliEMCALRunDepTempCalibCorrections"); + } + else + { // Else choose the one in the $ALICE_ROOT directory + if (fDebugLevel>0) AliInfo("Loading Recalib OADB from $ALICE_ROOT/OADB/EMCAL"); + + TFile *fRunDepRecalib= new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALTemperatureCorrCalib.root","read"); + if (!fRunDepRecalib || fRunDepRecalib->IsZombie()) + { + AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALTemperatureCorrCalib.root was not found"); + return 0; + } + + if (fRunDepRecalib) delete fRunDepRecalib; + + contRF->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALTemperatureCorrCalib.root","AliEMCALRunDepTempCalibCorrections"); + } + + TH1S *rundeprecal=(TH1S*)contRF->GetObject(runRC); + if (!rundeprecal) + { + AliError(Form("No Objects for run: %d",runRC)); + return 2; + } + + if (fDebugLevel>0) rundeprecal->Print(); + + Int_t nSM = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules(); + + for (Int_t ism=0; ismGetEMCALChannelRecalibrationFactor(ism,icol,irow); + + Int_t absID = fEMCALGeo->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor + factor *= rundeprecal->GetBinContent(absID) / 10000. ; // correction dependent on T + //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID, + // GetEMCALChannelRecalibrationFactor(ism,icol,irow) , rundeprecal->GetBinContent(absID) / 10000., factor); + fEMCALRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor); + } // columns + } // rows + } // SM loop + + return 1; +} + + +//_____________________________________________________ +Int_t AliEMCALTenderSupply::InitTimeCalibration() +{ + // Initialising bad channel maps + AliVEvent *event = GetEvent(); + + if (!event) + return 0; + + if (fDebugLevel>0) + AliInfo("Initialising time calibration map"); + + // init default maps first + if ( !fEMCALRecoUtils->GetEMCALTimeRecalibrationFactorsArray() ) + fEMCALRecoUtils->InitEMCALTimeRecalibrationFactors() ; + + Int_t runBC = event->GetRunNumber(); + + AliOADBContainer *contBC = new AliOADBContainer(""); + if (fBasePath!="") + { //if fBasePath specified in the ->SetBasePath() + if (fDebugLevel>0) AliInfo(Form("Loading time calibration OADB from given path %s",fBasePath.Data())); + + TFile *fbad=new TFile(Form("%s/EMCALTimeCalib.root",fBasePath.Data()),"read"); + if (!fbad || fbad->IsZombie()) + { + AliFatal(Form("EMCALTimeCalib.root was not found in the path provided: %s",fBasePath.Data())); + return 0; + } + + if (fbad) delete fbad; + + contBC->InitFromFile(Form("%s/EMCALTimeCalib.root",fBasePath.Data()),"AliEMCALTimeCalib"); + } + else + { // Else choose the one in the $ALICE_ROOT directory + if (fDebugLevel>0) AliInfo("Loading time calibration OADB from $ALICE_ROOT/OADB/EMCAL"); + + TFile *fbad=new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALTimeCalib.root","read"); + if (!fbad || fbad->IsZombie()) + { + AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALTimeCalib.root was not found"); + return 0; + } + + if (fbad) delete fbad; + + contBC->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALTimeCalib.root","AliEMCALTimeCalib"); + } + + TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runBC); + if (!arrayBC) + { + AliError(Form("No external time calibration set for run number: %d", runBC)); + return 2; + } + + // Here, it looks for a specific pass + TObjArray *arrayBCpass=(TObjArray*)arrayBC->FindObject(fFilepass); + if (!arrayBCpass) + { + AliError(Form("No external time calibration set for: %d -%s", runBC,fFilepass.Data())); + return 2; + } + + if (fDebugLevel>0) arrayBCpass->Print(); + + for( Int_t i = 0; i < 4; i++ ) + { + TH1F *h = fEMCALRecoUtils->GetEMCALChannelTimeRecalibrationFactors( i ); + if( h ) + delete h; + + h = (TH1F*)arrayBCpass->FindObject(Form("hAllTimeAvBC%d",i)); + + if (!h) + { + AliError(Form("Can not get hAllTimeAvBC%d",i)); + continue; + } + h->SetDirectory(0); + fEMCALRecoUtils->SetEMCALChannelTimeRecalibrationFactors(i,h); + } + return 1; +} + //_____________________________________________________ void AliEMCALTenderSupply::UpdateCells() { @@ -597,49 +1215,134 @@ void AliEMCALTenderSupply::UpdateCells() //Recalibrate energy and time cells //This is required for later reclusterization - AliESDCaloCells *cells = fTender->GetEvent()->GetEMCALCells(); - Int_t bunchCrossNo = fTender->GetEvent()->GetBunchCrossNumber(); + AliVEvent *event = GetEvent(); + + if(!event) return ; + + AliVCaloCells *cells = event->GetEMCALCells(); + Int_t bunchCrossNo = event->GetBunchCrossNumber(); - fEMCALRecoUtils->SwitchOnRecalibration(); - fEMCALRecoUtils->SwitchOnTimeRecalibration(); + fEMCALRecoUtils->RecalibrateCells(cells, bunchCrossNo); + + // remove exotic cells - loop through cells and zero the exotic ones + // just like with bad cell rejection in reco utils (inside RecalibrateCells) + if( fRejectExoticCells ) + { + Short_t absId =-1; + Double_t ecell = 0; + Double_t tcell = 0; + Double_t efrac = 0; + Short_t mclabel = -1; + Bool_t isExot = kFALSE; + + // loop through cells + Int_t nEMcell = cells->GetNumberOfCells() ; + for (Int_t iCell = 0; iCell < nEMcell; iCell++) + { + cells->GetCell( iCell, absId, ecell, tcell, mclabel, efrac ); + + isExot = fEMCALRecoUtils->IsExoticCell( absId, cells, bunchCrossNo ); + // zero if exotic + if( isExot ) + cells->SetCell( iCell, absId, 0.0, -1.0, mclabel, efrac ); + } // cell loop + } // reject exotic cells - fEMCALRecoUtils->RecalibrateCells(cells, bunchCrossNo); - + cells->Sort(); } //_____________________________________________________ -void AliEMCALTenderSupply::InitRecParam() +TString AliEMCALTenderSupply::GetBeamType() { - if (fDebugLevel>0) - AliInfo("Initialize the recParam"); - fRecParam = new AliEMCALRecParam; + // Get beam type : pp-AA-pA + // ESDs have it directly, AODs get it from hardcoded run number ranges + + AliVEvent *event = GetEvent(); + + if (!event) { + AliError("Couldn't retrieve event!"); + return ""; + } + + TString beamType; + + AliESDEvent *esd = dynamic_cast(event); + if (esd) { + const AliESDRun *run = esd->GetESDRun(); + beamType = run->GetBeamType(); + } + else + { + Int_t runNumber = event->GetRunNumber(); + if ((runNumber >= 136851 && runNumber <= 139517) // LHC10h + || (runNumber >= 166529 && runNumber <= 170593)) // LHC11h + { + beamType = "A-A"; + } + else + { + beamType = "p-p"; + } + } + + return beamType; } +//_____________________________________________________ +Int_t AliEMCALTenderSupply::InitRecParam() +{ + // Init reco params if not yet exist (probably shipped by the user already) + + if( fRecParam != 0 ) + return 2; + + TString beamType = GetBeamType(); + + // set some default reco params + fRecParam = new AliEMCALRecParam(); + fRecParam->SetClusteringThreshold(0.100); + fRecParam->SetMinECut(0.050); + fRecParam->SetTimeCut(250); + fRecParam->SetTimeMin(425); + fRecParam->SetTimeMax(825); + if ( beamType == "A-A") { + fRecParam->SetClusterizerFlag(AliEMCALRecParam::kClusterizerv2); + } + else { + fRecParam->SetClusterizerFlag(AliEMCALRecParam::kClusterizerv1); + } + + return 0; +} //_____________________________________________________ Bool_t AliEMCALTenderSupply::InitClusterization() { // Initialing clusterization/unfolding algorithm and set all the needed parameters. - AliESDEvent *event=fTender->GetEvent(); + AliVEvent *event = GetEvent(); + if (!event) return kFALSE; if (fDebugLevel>0) - AliInfo(Form("Initialising reclustering parameters: Clusterizer-%d",fRecParam->GetClusterizerFlag())); + AliInfo(Form("Initialising reclustering parameters: Clusterizer type: %d",fRecParam->GetClusterizerFlag())); //---setup clusterizer delete fClusterizer; if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1) fClusterizer = new AliEMCALClusterizerv1 (fEMCALGeo); - else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2) + else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2) fClusterizer = new AliEMCALClusterizerv2(fEMCALGeo); - else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) { + else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) + { AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fEMCALGeo); clusterizer->SetNRowDiff(fRecParam->GetNRowDiff()); clusterizer->SetNColDiff(fRecParam->GetNColDiff()); fClusterizer = clusterizer; - } else { + } + else + { AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag())); return kFALSE; } @@ -657,11 +1360,14 @@ Bool_t AliEMCALTenderSupply::InitClusterization() fClusterizer->SetJustClusters ( kTRUE ); // In case of unfolding after clusterization is requested, set the corresponding parameters - if (fRecParam->GetUnfold()) { - for (Int_t i = 0; i < 8; ++i) { + if (fRecParam->GetUnfold()) + { + for (Int_t i = 0; i < 8; ++i) + { fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i)); } - for (Int_t i = 0; i < 3; ++i) { + for (Int_t i = 0; i < 3; ++i) + { fClusterizer->SetPar5 (i, fRecParam->GetPar5(i)); fClusterizer->SetPar6 (i, fRecParam->GetPar6(i)); } @@ -679,34 +1385,37 @@ void AliEMCALTenderSupply::FillDigitsArray() { // Fill digits from cells to a TClonesArray. - AliESDEvent *event = fTender->GetEvent(); - if (!event) + AliVEvent *event = GetEvent(); + + if (!event) return; fDigitsArr->Clear("C"); - AliESDCaloCells *cells = event->GetEMCALCells(); + AliVCaloCells *cells = event->GetEMCALCells(); Int_t ncells = cells->GetNumberOfCells(); - for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) { - Double_t cellAmplitude=0, cellTime=0; - Short_t cellNumber=0; - if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE) + for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) + { + Double_t cellAmplitude=0, cellTime=0, efrac = 0; + Short_t cellNumber=0, mcLabel=-1; + + if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime, mcLabel, efrac) != kTRUE) break; // Do not add if already too low (some cells set to 0 if bad channels) - if(cellAmplitude < fRecParam->GetMinECut() ) + if (cellAmplitude < fRecParam->GetMinECut() ) continue; - - // If requested, do not include exotic clusters - if(fEMCALRecoUtils->IsExoticCell(cellNumber,cells,event->GetBunchCrossNumber())) + + // If requested, do not include exotic cells + if (fEMCALRecoUtils->IsExoticCell(cellNumber,cells,event->GetBunchCrossNumber())) continue; AliEMCALDigit *digit = static_cast(fDigitsArr->New(idigit)); digit->SetId(cellNumber); - digit->SetTime(cellTime); - digit->SetTimeR(cellTime); + digit->SetTime((Float_t)cellTime); + digit->SetTimeR((Float_t)cellTime); digit->SetIndexInList(idigit); digit->SetType(AliEMCALDigit::kHG); - digit->SetAmplitude(cellAmplitude); + digit->SetAmplitude((Float_t)cellAmplitude); idigit++; } } @@ -724,60 +1433,69 @@ void AliEMCALTenderSupply::UpdateClusters() { // Update ESD cluster list. - AliESDEvent *event = fTender->GetEvent(); + AliVEvent *event = GetEvent(); + if (!event) return; TClonesArray *clus = dynamic_cast(event->FindListObject("caloClusters")); if (!clus) clus = dynamic_cast(event->FindListObject("CaloClusters")); - if (!clus) { + if (!clus) + { AliError(" Null pointer to calo clusters array, returning"); return; } Int_t nents = clus->GetEntriesFast(); - for (Int_t i=0; i < nents; ++i) { - AliESDCaloCluster *c = static_cast(clus->At(i)); + for (Int_t i=0; i < nents; ++i) + { + AliVCluster *c = dynamic_cast(clus->At(i)); if (!c) continue; if (c->IsEMCAL()) delete clus->RemoveAt(i); } + clus->Compress(); + RecPoints2Clusters(clus); } //_____________________________________________________ void AliEMCALTenderSupply::RecPoints2Clusters(TClonesArray *clus) { - // Convert AliEMCALRecoPoints to AliESDCaloClusters. + // Convert AliEMCALRecoPoints to AliESDCaloClusters/AliAODCaloClusters. // Cluster energy, global position, cells and their amplitude fractions are restored. - AliESDEvent *event = fTender->GetEvent(); + AliVEvent *event = GetEvent(); + if (!event) return; Int_t ncls = fClusterArr->GetEntriesFast(); - for(Int_t i=0, nout=clus->GetEntriesFast(); i < ncls; ++i) { + for(Int_t i=0, nout=clus->GetEntriesFast(); i < ncls; ++i) + { AliEMCALRecPoint *recpoint = static_cast(fClusterArr->At(i)); - Int_t ncells_true = 0; + Int_t ncellsTrue = 0; const Int_t ncells = recpoint->GetMultiplicity(); UShort_t absIds[ncells]; Double32_t ratios[ncells]; Int_t *dlist = recpoint->GetDigitsList(); Float_t *elist = recpoint->GetEnergiesList(); - for (Int_t c = 0; c < ncells; ++c) { + for (Int_t c = 0; c < ncells; ++c) + { AliEMCALDigit *digit = static_cast(fDigitsArr->At(dlist[c])); - absIds[ncells_true] = digit->GetId(); - ratios[ncells_true] = elist[c]/digit->GetAmplitude(); - if (ratios[ncells_true] < 0.001) + absIds[ncellsTrue] = digit->GetId(); + ratios[ncellsTrue] = elist[c]/digit->GetAmplitude(); + if (ratios[ncellsTrue] < 0.001) continue; - ++ncells_true; + ++ncellsTrue; } - if (ncells_true < 1) { + if (ncellsTrue < 1) + { AliWarning("Skipping cluster with no cells"); continue; } @@ -788,12 +1506,12 @@ void AliEMCALTenderSupply::RecPoints2Clusters(TClonesArray *clus) Float_t g[3]; gpos.GetXYZ(g); - AliESDCaloCluster *c = static_cast(clus->New(nout++)); + AliVCluster *c = static_cast(clus->New(nout++)); c->SetID(nout-1); c->SetType(AliVCluster::kEMCALClusterv1); c->SetE(recpoint->GetEnergy()); c->SetPosition(g); - c->SetNCells(ncells_true); + c->SetNCells(ncellsTrue); c->SetDispersion(recpoint->GetDispersion()); c->SetEmcCpvDistance(-1); //not yet implemented c->SetChi2(-1); //not yet implemented @@ -803,10 +1521,8 @@ void AliEMCALTenderSupply::RecPoints2Clusters(TClonesArray *clus) recpoint->GetElipsAxis(elipAxis); c->SetM02(elipAxis[0]*elipAxis[0]) ; c->SetM20(elipAxis[1]*elipAxis[1]) ; - AliESDCaloCluster *cesd = static_cast(c); - cesd->SetCellsAbsId(absIds); - cesd->SetCellsAmplitudeFraction(ratios); - + c->SetCellsAbsId(absIds); + c->SetCellsAmplitudeFraction(ratios); } } @@ -818,7 +1534,8 @@ void AliEMCALTenderSupply::GetPass() AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); fInputTree = mgr->GetTree(); - if (!fInputTree) { + if (!fInputTree) + { AliError("Pointer to tree = 0, returning"); return; } @@ -830,12 +1547,12 @@ void AliEMCALTenderSupply::GetPass() } TString fname(fInputFile->GetName()); - if (fname.Contains("pass1")) fFilepass = TString("pass1"); - else if(fname.Contains("pass2")) fFilepass = TString("pass2"); - else if(fname.Contains("pass3")) fFilepass = TString("pass3"); - else { + if (fname.Contains("pass1")) fFilepass = TString("pass1"); + else if (fname.Contains("pass2")) fFilepass = TString("pass2"); + else if (fname.Contains("pass3")) fFilepass = TString("pass3"); + else + { AliError(Form("Pass number string not found: %s", fname.Data())); return; } } -