/************************************************************************** * 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. * **************************************************************************/ /////////////////////////////////////////////////////////////////////////////// // // // 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 // // // /////////////////////////////////////////////////////////////////////////////// #include #include #include #include #include #include #include #include "AliAODEvent.h" #include "AliAODMCParticle.h" #include "AliAnalysisManager.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 "AliESDCaloCluster.h" #include "AliESDEvent.h" #include "AliLog.h" #include "AliMagF.h" #include "AliOADBContainer.h" #include "AliTender.h" #include "AliEMCALTenderSupply.h" ClassImp(AliEMCALTenderSupply) AliEMCALTenderSupply::AliEMCALTenderSupply() : AliTenderSupply() ,fTask(0) ,fRun(0) ,fEMCALGeo(0x0) ,fEMCALGeoName("") ,fEMCALRecoUtils(0) ,fConfigName("") ,fDebugLevel(0) ,fNonLinearFunc(-1) ,fNonLinearThreshold(-1) ,fReCalibCluster(kFALSE) ,fUpdateCell(kFALSE) ,fCalibrateEnergy(kFALSE) ,fCalibrateTime(kFALSE) ,fCalibrateTimeParamAvailable(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) ,fGetPassFromFileName(kTRUE) ,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) ,fSetCellMCLabelFromCluster(0) ,fTempClusterArr(0) ,fRemapMCLabelForAODs(0) ,fUseAutomaticRecalib(1) ,fUseAutomaticRunDepRecalib(1) ,fUseAutomaticTimeCalib(1) ,fUseAutomaticRecParam(1) { // Default constructor. for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ; for(Int_t j = 0; j < 12672; j++) fOrgClusterCellId[j] = -1; } //_____________________________________________________ AliEMCALTenderSupply::AliEMCALTenderSupply(const char *name, const AliTender *tender) : AliTenderSupply(name,tender) ,fTask(0) ,fRun(0) ,fEMCALGeo(0x0) ,fEMCALGeoName("") ,fEMCALRecoUtils(0) ,fConfigName("") ,fDebugLevel(0) ,fNonLinearFunc(-1) ,fNonLinearThreshold(-1) ,fReCalibCluster(kFALSE) ,fUpdateCell(kFALSE) ,fCalibrateEnergy(kFALSE) ,fCalibrateTime(kFALSE) ,fCalibrateTimeParamAvailable(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) ,fGetPassFromFileName(kTRUE) ,fFilepass("") ,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) ,fSetCellMCLabelFromCluster(0) ,fTempClusterArr(0) ,fRemapMCLabelForAODs(0) ,fUseAutomaticRecalib(1) ,fUseAutomaticRunDepRecalib(1) ,fUseAutomaticTimeCalib(1) ,fUseAutomaticRecParam(1) { // Named constructor for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ; for(Int_t j = 0; j < 12672; j++) fOrgClusterCellId[j] = -1; } //_____________________________________________________ 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) ,fCalibrateTimeParamAvailable(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) ,fGetPassFromFileName(kTRUE) ,fFilepass("") ,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) ,fSetCellMCLabelFromCluster(0) ,fTempClusterArr(0) ,fRemapMCLabelForAODs(0) ,fUseAutomaticRecalib(1) ,fUseAutomaticRunDepRecalib(1) ,fUseAutomaticTimeCalib(1) ,fUseAutomaticRecParam(1) { // Named constructor. for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ; for(Int_t j = 0; j < 12672; j++) fOrgClusterCellId[j] = -1; } //_____________________________________________________ AliEMCALTenderSupply::~AliEMCALTenderSupply() { //Destructor if (!AliAnalysisManager::GetAnalysisManager()) return; if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { delete fEMCALRecoUtils; delete fRecParam; delete fUnfolder; if (!fClusterizer) { if (fDigitsArr) { fDigitsArr->Clear("C"); delete fDigitsArr; } } else { delete fClusterizer; fDigitsArr = 0; } } } //_____________________________________________________ void AliEMCALTenderSupply::SetDefaults() { // Set default settings. SwitchOnReclustering(); 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) AliWarning("Init EMCAL Tender supply"); 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; fEMCALGeoName = tender->fEMCALGeoName; fEMCALRecoUtils = tender->fEMCALRecoUtils; fConfigName = tender->fConfigName; fNonLinearFunc = tender->fNonLinearFunc; fNonLinearThreshold = tender->fNonLinearThreshold; fReCalibCluster = tender->fReCalibCluster; fUpdateCell = tender->fUpdateCell; fRecalClusPos = tender->fRecalClusPos; fCalibrateEnergy = tender->fCalibrateEnergy; fCalibrateTime = tender->fCalibrateTime; fCalibrateTimeParamAvailable = tender->fCalibrateTimeParamAvailable; 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; 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] ; } 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 if (!fEMCALRecoUtils) 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); // spatial cut based on separate eta/phi or common processing if (fCutEtaPhiSum) { fEMCALRecoUtils->SwitchOnCutEtaPhiSum(); if (fRcut != -1) fEMCALRecoUtils->SetCutR(fRcut); } else if (fCutEtaPhiSeparate) { fEMCALRecoUtils->SwitchOnCutEtaPhiSeparate(); 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. AliVEvent *event = GetEvent(); if (!event) { AliError("Event ptr = 0, returning"); return; } // Initialising parameters once per run number if (RunChanged()) { fRun = event->GetRunNumber(); AliWarning(Form("Run changed, initializing parameters for %d", fRun)); if (dynamic_cast(event)) { AliWarning("============================================================="); AliWarning("=== Running on AOD is not equivalent to running on ESD! ==="); AliWarning("============================================================="); } // 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 if (fGetPassFromFileName) GetPass(); // 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; // 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())); } // init recalibration factors if (needRecalib) { if(fUseAutomaticRecalib) { Int_t fInitRecalib = InitRecalib(); if (fInitRecalib==0) AliError("InitRecalib returned false, returning"); if (fInitRecalib==1) AliWarning("InitRecalib OK"); if (fInitRecalib>1) AliWarning(Form("No recalibration available: %d - %s", event->GetRunNumber(), fFilepass.Data())); } if(fUseAutomaticRunDepRecalib) { 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())); } } // init time calibration if (needTimecalib && fUseAutomaticTimeCalib) { Int_t initTC = InitTimeCalibration(); if (!initTC) AliError("InitTimeCalibration returned false, returning"); if (initTC==1) { fCalibrateTimeParamAvailable = kTRUE; AliWarning("InitTimeCalib OK"); } if (initTC > 1) AliWarning(Form("No external time calibration available: %d - %s", event->GetRunNumber(), fFilepass.Data())); } // init misalignment matrix if (needMisalign) { if (!InitMisalignMatrix()) AliError("InitMisalignmentMatrix returned false, returning"); else AliWarning("InitMisalignMatrix OK"); } // initiate reco params with some defaults // will not overwrite, if those have been provided by user if (needRecoParam && fUseAutomaticRecParam) { Int_t initRC = InitRecParam(); if (initRC == 0) AliInfo("Defaults reco params loaded."); if (initRC > 1) AliWarning("User defined reco params."); } // init clusterizer if (needClusterizer) { if (!InitClusterization()) AliError("InitClusterization returned false, returning"); 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 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) 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(); // 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(); // 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(); if (fDoUpdateOnly) return; } // RECLUSTERIZATION --------------------------------------------------------- if (fReClusterize) { FillDigitsArray(); Clusterize(); UpdateClusters(); } // Store good clusters TClonesArray *clusArr = dynamic_cast(event->FindListObject("caloClusters")); if (!clusArr) clusArr = dynamic_cast(event->FindListObject("CaloClusters")); if (!clusArr) { 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) { AliVCluster *clust = static_cast(clusArr->At(icluster)); if (!clust) continue; if (!clust->IsEMCAL()) continue; // 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? } } // 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 } } // 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 (clust->E() < 1e-9) { delete clusArr->RemoveAt(icluster); continue; } } // CLUSTER POSITION ------------------------------------------- // does process local cell energy recalibration, if enabled and cells // not calibrated yet if (fRecalClusPos) fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, cells, clust); // 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(); if (!fDoTrackMatch) return; // TRACK MATCHING ----------------------------------------------------------- if (!TGeoGlobalMagField::Instance()->GetField()) { event->InitMagneticField(); } fEMCALRecoUtils->FindMatches(event,0x0,fEMCALGeo); fEMCALRecoUtils->SetClusterMatchedToTrack(event); fEMCALRecoUtils->SetTracksMatchedToCluster(event); } //_____________________________________________________ Bool_t AliEMCALTenderSupply::InitMisalignMatrix() { // Initialising misalignment matrices AliVEvent *event = GetEvent(); if (!event) return kFALSE; if (fGeomMatrixSet) { AliInfo("Misalignment matrix already set"); return kTRUE; } if (fDebugLevel>0) AliInfo("Initialising misalignment matrix"); if (fLoadGeomMatrices) { for(Int_t mod=0; mod < fEMCALGeo->GetNumberOfSuperModules(); ++mod) { if (fEMCALMatrix[mod]){ if (fDebugLevel > 2) fEMCALMatrix[mod]->Print(); fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod); } } fGeomMatrixSet = kTRUE; return kTRUE; } Int_t runGM = event->GetRunNumber(); TObjArray *mobj = 0; 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 == 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 (!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; } //_____________________________________________________ Int_t AliEMCALTenderSupply::InitBadChannels() { // Initialising bad channel maps AliVEvent *event = GetEvent(); if (!event) return 0; if (fDebugLevel>0) AliInfo("Initialising Bad channel map"); // 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() 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 || fbad->IsZombie()) { AliFatal(Form("EMCALBadChannels.root was not found in the path provided: %s",fBasePath.Data())); return 0; } 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"); } TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runBC); if (!arrayBC) { AliError(Form("No external hot channel set for run number: %d", runBC)); return 2; } Int_t sms = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules(); for (Int_t i=0; iGetEMCALChannelStatusMap(i); if (h) delete h; h=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i)); 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. 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 *fRecalib= new TFile(Form("%s/EMCALRecalib.root",fBasePath.Data()),"read"); if (!fRecalib || fRecalib->IsZombie()) { AliFatal(Form("EMCALRecalib.root not found in %s",fBasePath.Data())); return 0; } 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 || 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); if (!recal) { AliError(Form("No Objects for run: %d",runRC)); return 2; } TObjArray *recalpass=(TObjArray*)recal->FindObject(fFilepass); if (!recalpass) { AliError(Form("No Objects for run: %d - %s",runRC,fFilepass.Data())); return 2; } TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib"); if (!recalib) { AliError(Form("No Recalib histos found for %d - %s",runRC,fFilepass.Data())); return 2; } 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) { AliError(Form("Could not load EMCALRecalFactors_SM%d",i)); continue; } h->SetDirectory(0); fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(i,h); } 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) { AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runRC)); // let's get the closest runnumber instead then.. Int_t lower = 0; Int_t ic = 0; Int_t maxEntry = contRF->GetNumberOfEntries(); while ((ic < maxEntry) && (contRF->UpperLimit(ic) < runRC)) { lower = ic; ic++; } Int_t closest = lower; if ((icLowerLimit(ic)-runRC) < (runRC - contRF->UpperLimit(lower))) { closest = ic; } AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRF->LowerLimit(closest))); rundeprecal = (TH1S*) contRF->GetObjectByIndex(closest); } 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 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 TString pass = fFilepass; if (fFilepass=="calo_spc") pass ="pass1"; TObjArray *arrayBCpass=(TObjArray*)arrayBC->FindObject(pass); if (!arrayBCpass) { AliError(Form("No external time calibration set for: %d -%s", runBC,pass.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() { //Remove bad cells from the cell list //Recalibrate energy and time cells //This is required for later reclusterization AliVEvent *event = GetEvent(); if (!event) return ; AliVCaloCells *cells = event->GetEMCALCells(); Int_t bunchCrossNo = event->GetBunchCrossNumber(); 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; Int_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 cells->Sort(); } //_____________________________________________________ TString AliEMCALTenderSupply::GetBeamType() { // 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); if (!fCalibrateTimeParamAvailable) { fRecParam->SetTimeCut(250*1.e-9); fRecParam->SetTimeMin(425*1.e-9); fRecParam->SetTimeMax(825*1.e-9); } else { fRecParam->SetTimeCut(100*1.e-9); fRecParam->SetTimeMin(-50*1.e-9); fRecParam->SetTimeMax(50*1.e-9); } 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. AliVEvent *event = GetEvent(); if (!event) return kFALSE; if (fDebugLevel>0) AliInfo(Form("Initialising reclustering parameters: Clusterizer type: %d",fRecParam->GetClusterizerFlag())); //---setup clusterizer if (fClusterizer) { // avoid deleting fDigitsArr fClusterizer->SetDigitsArr(0); delete fClusterizer; fClusterizer = 0; } if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1) fClusterizer = new AliEMCALClusterizerv1 (fEMCALGeo); else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2) fClusterizer = new AliEMCALClusterizerv2(fEMCALGeo); else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) { AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fEMCALGeo); clusterizer->SetNRowDiff(fRecParam->GetNRowDiff()); clusterizer->SetNColDiff(fRecParam->GetNColDiff()); fClusterizer = clusterizer; } else { AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag())); return kFALSE; } // Set the clustering parameters fClusterizer->SetECAClusteringThreshold(fRecParam->GetClusteringThreshold()); fClusterizer->SetECALogWeight (fRecParam->GetW0() ); fClusterizer->SetMinECut (fRecParam->GetMinECut() ); fClusterizer->SetUnfolding (fRecParam->GetUnfold() ); fClusterizer->SetECALocalMaxCut (fRecParam->GetLocMaxCut() ); fClusterizer->SetTimeCut (fRecParam->GetTimeCut() ); fClusterizer->SetTimeMin (fRecParam->GetTimeMin() ); fClusterizer->SetTimeMax (fRecParam->GetTimeMax() ); fClusterizer->SetInputCalibrated (kTRUE ); 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) { fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i)); } for (Int_t i = 0; i < 3; ++i) { fClusterizer->SetPar5 (i, fRecParam->GetPar5(i)); fClusterizer->SetPar6 (i, fRecParam->GetPar6(i)); } fClusterizer->InitClusterUnfolding(); } fClusterizer->SetDigitsArr(fDigitsArr); fClusterizer->SetOutput(0); fClusterArr = const_cast(fClusterizer->GetRecPoints()); return kTRUE; } //_____________________________________________________ void AliEMCALTenderSupply::FillDigitsArray() { // Fill digits from cells to a TClonesArray. AliVEvent *event = GetEvent(); if (!event) return; // In case of MC productions done before aliroot tag v5-02-Rev09 // assing the cluster label to all the cells belonging to this cluster // very rough Int_t cellLabels[12672]; if (fSetCellMCLabelFromCluster) { for (Int_t i = 0; i < 12672; i++) { cellLabels [i] = 0 ; fOrgClusterCellId[i] =-1 ; } Int_t nClusters = event->GetNumberOfCaloClusters(); for (Int_t i = 0; i < nClusters; i++) { AliVCluster *clus = event->GetCaloCluster(i); if (!clus) continue; if (!clus->IsEMCAL()) continue ; Int_t label = clus->GetLabel(); UShort_t * index = clus->GetCellsAbsId() ; for(Int_t icell=0; icell < clus->GetNCells(); icell++) { cellLabels [index[icell]] = label; fOrgClusterCellId[index[icell]] = i ; // index of the original cluster } }// cluster loop } fDigitsArr->Clear("C"); 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, efrac = 0; Short_t cellNumber=0; Int_t mcLabel=-1; if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime, mcLabel, efrac) != kTRUE) break; // Do not add if energy already too low (some cells set to 0 if bad channels) if (cellAmplitude < fRecParam->GetMinECut()) continue; // If requested, do not include exotic cells if (fEMCALRecoUtils->IsExoticCell(cellNumber,cells,event->GetBunchCrossNumber())) continue; if (fSetCellMCLabelFromCluster) mcLabel = cellLabels[cellNumber]; else if (fRemapMCLabelForAODs ) RemapMCLabelForAODs(mcLabel); if (mcLabel > 0 && efrac < 1e-6) efrac = 1; new((*fDigitsArr)[idigit]) AliEMCALDigit(mcLabel, mcLabel, cellNumber, (Float_t)cellAmplitude, (Float_t)cellTime, AliEMCALDigit::kHG,idigit, 0, 0, efrac*cellAmplitude); idigit++; } } //_____________________________________________________ void AliEMCALTenderSupply::Clusterize() { // Clusterize. fClusterizer->Digits2Clusters(""); } //_____________________________________________________ void AliEMCALTenderSupply::UpdateClusters() { // Update ESD cluster list. AliVEvent *event = GetEvent(); if (!event) return; TClonesArray *clus = dynamic_cast(event->FindListObject("caloClusters")); if (!clus) clus = dynamic_cast(event->FindListObject("CaloClusters")); if (!clus) { AliError(" Null pointer to calo clusters array, returning"); return; } // Before destroying the orignal list, assign to the rec points the MC labels // of the original clusters, if requested if (fSetCellMCLabelFromCluster == 2) SetClustersMCLabelFromOriginalClusters() ; Int_t nents = clus->GetEntriesFast(); 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/AliAODCaloClusters. // Cluster energy, global position, cells and their amplitude fractions are restored. AliVEvent *event = GetEvent(); if (!event) return; Int_t ncls = fClusterArr->GetEntriesFast(); for(Int_t i=0, nout=clus->GetEntriesFast(); i < ncls; ++i) { AliEMCALRecPoint *recpoint = static_cast(fClusterArr->At(i)); 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) { AliEMCALDigit *digit = static_cast(fDigitsArr->At(dlist[c])); absIds[ncellsTrue] = digit->GetId(); ratios[ncellsTrue] = elist[c]/digit->GetAmplitude(); if (ratios[ncellsTrue] < 0.001) continue; ++ncellsTrue; } if (ncellsTrue < 1) { AliWarning("Skipping cluster with no cells"); continue; } // calculate new cluster position TVector3 gpos; recpoint->GetGlobalPosition(gpos); Float_t g[3]; gpos.GetXYZ(g); AliVCluster *c = static_cast(clus->New(nout++)); c->SetID(nout-1); c->SetType(AliVCluster::kEMCALClusterv1); c->SetE(recpoint->GetEnergy()); c->SetPosition(g); c->SetNCells(ncellsTrue); c->SetDispersion(recpoint->GetDispersion()); c->SetEmcCpvDistance(-1); //not yet implemented c->SetChi2(-1); //not yet implemented c->SetTOF(recpoint->GetTime()) ; //time-of-flight c->SetNExMax(recpoint->GetNExMax()); //number of local maxima Float_t elipAxis[2]; recpoint->GetElipsAxis(elipAxis); c->SetM02(elipAxis[0]*elipAxis[0]) ; c->SetM20(elipAxis[1]*elipAxis[1]) ; c->SetCellsAbsId(absIds); c->SetCellsAmplitudeFraction(ratios); //MC labels Int_t parentMult = 0; Int_t *parentList = recpoint->GetParents(parentMult); if (parentMult > 0) c->SetLabel(parentList, parentMult); } } //___________________________________________________________ void AliEMCALTenderSupply::RemapMCLabelForAODs(Int_t & label) { // MC label for Cells not remapped after ESD filtering, do it here. if (label < 0) return; AliAODEvent * evt = dynamic_cast (GetEvent()) ; if (!evt) return ; TClonesArray * arr = dynamic_cast(evt->FindListObject("mcparticles")) ; if (!arr) return ; if (label < arr->GetEntriesFast()) { AliAODMCParticle * particle = dynamic_cast(arr->At(label)); if (!particle) return ; if (label == particle->Label()) return ; // label already OK } // loop on the particles list and check if there is one with the same label for (Int_t ind = 0; ind < arr->GetEntriesFast(); ind++) { AliAODMCParticle * particle = dynamic_cast(arr->At(ind)); if (!particle) continue ; if (label == particle->Label()) { label = ind; return; } } label = -1; } //_____________________________________________________________________________________________ void AliEMCALTenderSupply::SetClustersMCLabelFromOriginalClusters() { // Get the original clusters that contribute to the new rec point cluster, // assign the labels of such clusters to the new rec point cluster. // Only approximatedly valid when output are V1 clusters, or input/output clusterizer // are the same handle with care // Copy from same method in AliAnalysisTaskEMCALClusterize, but here modify the recpoint and // not the output calocluster Int_t ncls = fClusterArr->GetEntriesFast(); for(Int_t irp=0; irp < ncls; ++irp) { TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ... clArray.Reset(); Int_t nClu = 0; Int_t nLabTotOrg = 0; Float_t emax = -1; Int_t idMax = -1; AliEMCALRecPoint *clus = static_cast(fClusterArr->At(irp)); //Find the clusters that originally had the cells const Int_t ncells = clus->GetMultiplicity(); Int_t *digList = clus->GetDigitsList(); for (Int_t iLoopCell = 0 ; iLoopCell < ncells ; iLoopCell++) { AliEMCALDigit *digit = static_cast(fDigitsArr->At(digList[iLoopCell])); Int_t idCell = digit->GetId(); if (idCell>=0) { Int_t idCluster = fOrgClusterCellId[idCell]; Bool_t set = kTRUE; for (Int_t icl =0; icl < nClu; icl++) { if (((Int_t)clArray.GetAt(icl))==-1) continue; if (idCluster == ((Int_t)clArray.GetAt(icl))) set = kFALSE; } if (set && idCluster >= 0) { clArray.SetAt(idCluster,nClu++); nLabTotOrg+=(GetEvent()->GetCaloCluster(idCluster))->GetNLabels(); //Search highest E cluster AliVCluster * clOrg = GetEvent()->GetCaloCluster(idCluster); if (emax < clOrg->E()) { emax = clOrg->E(); idMax = idCluster; } } } }// cell loop // Put the first in the list the cluster with highest energy if (idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position { Int_t maxIndex = -1; Int_t firstCluster = ((Int_t)clArray.GetAt(0)); for (Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++) { if (idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster; } if (firstCluster >=0 && idMax >=0) { clArray.SetAt(idMax,0); clArray.SetAt(firstCluster,maxIndex); } } // Get the labels list in the original clusters, assign all to the new cluster TArrayI clMCArray(nLabTotOrg) ; clMCArray.Reset(); Int_t nLabTot = 0; for (Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++) { Int_t idCluster = (Int_t) clArray.GetAt(iLoopCluster); AliVCluster * clOrg = GetEvent()->GetCaloCluster(idCluster); Int_t nLab = clOrg->GetNLabels(); for (Int_t iLab = 0 ; iLab < nLab ; iLab++) { Int_t lab = clOrg->GetLabelAt(iLab) ; if (lab>=0) { Bool_t set = kTRUE; for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++) { if (lab == ((Int_t)clMCArray.GetAt(iLabTot))) set = kFALSE; } if (set) clMCArray.SetAt(lab,nLabTot++); } } }// cluster loop // Set the final list of labels to rec point Int_t *labels = new Int_t[nLabTot]; for(Int_t il = 0; il < nLabTot; il++) labels[il] = clMCArray.GetArray()[il]; clus->SetParents(nLabTot,labels); } // rec point array } //_____________________________________________________ void AliEMCALTenderSupply::GetPass() { // Get passx from filename AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); fInputTree = mgr->GetTree(); if (!fInputTree) { AliError("Pointer to tree = 0, returning"); return; } fInputFile = fInputTree->GetCurrentFile(); if (!fInputFile) { AliError("Null pointer input file, returning"); return; } 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("pass4")) fFilepass = TString("pass4"); else if (fname.Contains("pass5")) fFilepass = TString("pass5"); else if (fname.Contains("LHC11c") && fname.Contains("spc_calo")) fFilepass = TString("spc_calo"); else if (fname.Contains("calo") || fname.Contains("high_lumi")) { //printf("AliEMCALTenderSupply::GetPass() - Path contains or , set as \n"); fFilepass = TString("pass1"); } else if (fname.Contains("LHC14a1a")) fFilepass = TString("LHC14a1a"); else { AliError(Form("Pass number string not found: %s", fname.Data())); return; } }