From 09e819c98763a50c84f16c283a3550cdca98f0e1 Mon Sep 17 00:00:00 2001 From: gconesab Date: Sat, 12 Jun 2010 13:48:26 +0000 Subject: [PATCH] Added possibility to recalibrate CaloClusters during the analysis --- PWG4/PartCorrBase/AliCaloTrackAODReader.cxx | 14 +- PWG4/PartCorrBase/AliCaloTrackESDReader.cxx | 7 + PWG4/PartCorrBase/AliCaloTrackReader.cxx | 8 +- PWG4/PartCorrBase/AliCalorimeterUtils.cxx | 155 +++++++++++++++++++- PWG4/PartCorrBase/AliCalorimeterUtils.h | 67 +++++++-- 5 files changed, 230 insertions(+), 21 deletions(-) diff --git a/PWG4/PartCorrBase/AliCaloTrackAODReader.cxx b/PWG4/PartCorrBase/AliCaloTrackAODReader.cxx index 4ada8ac2ede..9347e4a71f4 100755 --- a/PWG4/PartCorrBase/AliCaloTrackAODReader.cxx +++ b/PWG4/PartCorrBase/AliCaloTrackAODReader.cxx @@ -167,6 +167,12 @@ void AliCaloTrackAODReader::FillInputEMCAL() { if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackAODReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n", momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta()); + //Recalibrate the cluster energy + if(GetCaloUtils()->IsRecalibrationOn()) { + Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetEMCALCells()); + clus->SetE(energy); + } + if(fWriteOutputStdAOD){ AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*clus); fAODEMCAL->Add(newclus); @@ -239,7 +245,13 @@ void AliCaloTrackAODReader::FillInputPHOS() { if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackAODReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n", momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta()); - + + //Recalibrate the cluster energy + if(GetCaloUtils()->IsRecalibrationOn()) { + Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells()); + clus->SetE(energy); + } + if(fWriteOutputStdAOD){ AliAODCaloCluster * newclus = new((*(fOutputEvent->GetCaloClusters()))[naod++])AliAODCaloCluster(*clus); fAODPHOS->Add(newclus); diff --git a/PWG4/PartCorrBase/AliCaloTrackESDReader.cxx b/PWG4/PartCorrBase/AliCaloTrackESDReader.cxx index f1b73a6f1a1..fb975eaee39 100755 --- a/PWG4/PartCorrBase/AliCaloTrackESDReader.cxx +++ b/PWG4/PartCorrBase/AliCaloTrackESDReader.cxx @@ -222,6 +222,10 @@ void AliCaloTrackESDReader::FillInputEMCAL() { Float_t energy = clus->E(); Char_t ttype= AliAODCluster::kEMCALClusterv1; + + //Recalibrate the cluster energy + if(GetCaloUtils()->IsRecalibrationOn()) energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliESDCaloCells*)GetEMCALCells()); + //Put new aod object in file in AOD calo clusters array AliAODCaloCluster *caloCluster = new((*(fOutputEvent->GetCaloClusters()))[naod++]) @@ -320,6 +324,9 @@ void AliCaloTrackESDReader::FillInputPHOS() { pid[AliPID::kProton]+pid[AliPID::kKaon]+pid[AliPID::kPion]; if( wCharged > wNeutral) ttype= AliAODCluster::kPHOSCharged; + //Recalibrate the cluster energy + if(GetCaloUtils()->IsRecalibrationOn()) energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliESDCaloCells*)GetPHOSCells()); + //Put new aod object in file in AOD calo clusters array AliAODCaloCluster *caloCluster = new((*(fOutputEvent->GetCaloClusters()))[naod++]) AliAODCaloCluster(id,nLabel,labels,energy, pos, NULL, ttype, 0); diff --git a/PWG4/PartCorrBase/AliCaloTrackReader.cxx b/PWG4/PartCorrBase/AliCaloTrackReader.cxx index 859d3874f54..4ea9a64f077 100755 --- a/PWG4/PartCorrBase/AliCaloTrackReader.cxx +++ b/PWG4/PartCorrBase/AliCaloTrackReader.cxx @@ -475,12 +475,14 @@ Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * curre } } - + + if(fFillEMCALCells) FillInputEMCALCells(); + if(fFillPHOSCells) FillInputPHOSCells(); + if(fFillCTS) FillInputCTS(); if(fFillEMCAL) FillInputEMCAL(); if(fFillPHOS) FillInputPHOS(); - if(fFillEMCALCells) FillInputEMCALCells(); - if(fFillPHOSCells) FillInputPHOSCells(); + return kTRUE ; } diff --git a/PWG4/PartCorrBase/AliCalorimeterUtils.cxx b/PWG4/PartCorrBase/AliCalorimeterUtils.cxx index b205280c750..8524c32e507 100755 --- a/PWG4/PartCorrBase/AliCalorimeterUtils.cxx +++ b/PWG4/PartCorrBase/AliCalorimeterUtils.cxx @@ -49,7 +49,8 @@ ClassImp(AliCalorimeterUtils) fRemoveBadChannels(kFALSE), fEMCALBadChannelMap(new TObjArray()),fPHOSBadChannelMap(new TObjArray()), fNCellsFromEMCALBorder(0), fNCellsFromPHOSBorder(0), - fNoEMCALBorderAtEta0(kFALSE) + fNoEMCALBorderAtEta0(kFALSE),fRecalibration(kFALSE), + fEMCALRecalibrationFactors(new TObjArray()), fPHOSRecalibrationFactors(new TObjArray()) { //Ctor @@ -70,7 +71,10 @@ AliCalorimeterUtils::AliCalorimeterUtils(const AliCalorimeterUtils & calo) : fPHOSBadChannelMap(new TObjArray(*calo.fPHOSBadChannelMap)), fNCellsFromEMCALBorder(calo.fNCellsFromEMCALBorder), fNCellsFromPHOSBorder(calo.fNCellsFromPHOSBorder), - fNoEMCALBorderAtEta0(calo.fNoEMCALBorderAtEta0) + fNoEMCALBorderAtEta0(calo.fNoEMCALBorderAtEta0), + fRecalibration(calo.fRecalibration), + fEMCALRecalibrationFactors(new TObjArray(*calo.fEMCALRecalibrationFactors)), + fPHOSRecalibrationFactors(new TObjArray(*calo.fEMCALRecalibrationFactors)) { // cpy ctor } @@ -92,6 +96,15 @@ AliCalorimeterUtils::~AliCalorimeterUtils() { delete fPHOSBadChannelMap; } + if(fEMCALRecalibrationFactors) { + fEMCALRecalibrationFactors->Clear(); + delete fEMCALRecalibrationFactors; + } + if(fPHOSRecalibrationFactors) { + fPHOSRecalibrationFactors->Clear(); + delete fPHOSRecalibrationFactors; + } + } //_______________________________________________________________ @@ -531,6 +544,54 @@ void AliCalorimeterUtils::InitPHOSBadChannelStatusMap(){ TH1::AddDirectory(oldStatus); } +//________________________________________________________________ +void AliCalorimeterUtils::InitEMCALRecalibrationFactors(){ + //Init EMCAL recalibration factors + if(fDebug > 0 )printf("AliCalorimeterUtils::InitEMCALRecalibrationFactors()\n"); + //In order to avoid rewriting the same histograms + Bool_t oldStatus = TH1::AddDirectoryStatus(); + TH1::AddDirectory(kFALSE); + + for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24)); + //Init the histograms with 1 + for (Int_t sm = 0; sm < 12; sm++) { + for (Int_t i = 0; i < 48; i++) { + for (Int_t j = 0; j < 24; j++) { + SetEMCALChannelRecalibrationFactor(sm,i,j,1.); + } + } + } + fEMCALRecalibrationFactors->SetOwner(kTRUE); + fEMCALRecalibrationFactors->Compress(); + + //In order to avoid rewriting the same histograms + TH1::AddDirectory(oldStatus); +} + +//________________________________________________________________ +void AliCalorimeterUtils::InitPHOSRecalibrationFactors(){ + //Init EMCAL recalibration factors + if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n"); + //In order to avoid rewriting the same histograms + Bool_t oldStatus = TH1::AddDirectoryStatus(); + TH1::AddDirectory(kFALSE); + for (int i = 0; i < 5; i++)fPHOSRecalibrationFactors->Add(new TH2F(Form("PHOSRecalFactors_Mod%d",i),Form("PHOSRecalFactors_Mod%d",i), 56, 0, 56, 64, 0, 64)); + //Init the histograms with 1 + for (Int_t m = 0; m < 5; m++) { + for (Int_t i = 0; i < 56; i++) { + for (Int_t j = 0; j < 64; j++) { + SetPHOSChannelRecalibrationFactor(m,i,j,1.); + } + } + } + fPHOSRecalibrationFactors->SetOwner(kTRUE); + fPHOSRecalibrationFactors->Compress(); + + //In order to avoid rewriting the same histograms + TH1::AddDirectory(oldStatus); +} + + //________________________________________________________________ void AliCalorimeterUtils::InitEMCALGeometry() { @@ -572,10 +633,100 @@ void AliCalorimeterUtils::Print(const Option_t * opt) const printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n", fNCellsFromEMCALBorder, fNCellsFromPHOSBorder); if(fNoEMCALBorderAtEta0) printf("Do not remove EMCAL clusters at Eta = 0\n"); + printf("Recalibrate Clusters? %d\n",fRecalibration); printf(" \n") ; } +//________________________________________________________________ +Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliESDCaloCluster * cluster, AliESDCaloCells * cells){ + // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster. + // ESD case + + if(!cells) { + printf("AliCalorimeterUtils::RecalibrateClusterEnergy(ESD) - Cells pointer does not exist, stop!"); + abort(); + } + //Get the cluster number of cells and list of absId, check what kind of cluster do we have. + UShort_t * index = cluster->GetCellsAbsId() ; + Double_t * fraction = cluster->GetCellsAmplitudeFraction() ; + Int_t ncells = cluster->GetNCells(); + TString calo = "EMCAL"; + if(cluster->IsPHOS()) calo = "PHOS"; + //Initialize some used variables + Float_t energy = 0; + Int_t absId = -1; + Int_t icol = -1, irow = -1, iRCU = -1, module=1; + Float_t factor = 1, frac = 0; + + //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy + for(Int_t icell = 0; icell < ncells; icell++){ + absId = index[icell]; + frac = fraction[icell]; + if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used. + module = GetModuleNumberCellIndexes(absId,calo,icol,irow,iRCU); + if(cluster->IsPHOS()) factor = GetPHOSChannelRecalibrationFactor (module,icol,irow); + else factor = GetEMCALChannelRecalibrationFactor(module,icol,irow); + if(fDebug>2) + printf("AliCalorimeterUtils::RecalibrateClusterEnergy(ESD) - recalibrate cell: %s, module %d, col %d, row %d, cell fraction %f, recalibration factor %f, cell energy %f\n", + calo.Data(),module,icol,irow,frac,factor,cells->GetCellAmplitude(absId)); + + energy += cells->GetCellAmplitude(absId)*factor*frac; + } + + if(fDebug>1) + printf("AliCalorimeterUtils::RecalibrateClusterEnergy(ESD) - Energy before %f, after %f\n",cluster->E(),energy); + + return energy; + +} + +//________________________________________________________________ +Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliAODCaloCluster * cluster, AliAODCaloCells * cells){ + // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster. + // AOD case + + if(!cells) { + printf("AliCalorimeterUtils::RecalibrateClusterEnergy(AOD) - Cells pointer does not exist, stop!"); + abort(); + } + + //Get the cluster number of cells and list of absId, check what kind of cluster do we have. + UShort_t * index = cluster->GetCellsAbsId() ; + Double_t * fraction = cluster->GetCellsAmplitudeFraction() ; + Int_t ncells = cluster->GetNCells(); + TString calo = "EMCAL"; + if(cluster->IsPHOSCluster()) calo = "PHOS"; + + //Initialize some used variables + Float_t energy = 0; + Int_t absId = -1; + Int_t icol = -1, irow = -1, iRCU = -1, module=1; + Float_t factor = 1, frac = 0; + + //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy + for(Int_t icell = 0; icell < ncells; icell++){ + absId = index[icell]; + frac = fraction[icell]; + if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used. + module = GetModuleNumberCellIndexes(absId,calo,icol,irow,iRCU); + if(cluster->IsPHOSCluster()) factor = GetPHOSChannelRecalibrationFactor (module,icol,irow); + else factor = GetEMCALChannelRecalibrationFactor(module,icol,irow); + if(fDebug>2) + printf("AliCalorimeterUtils::RecalibrateClusterEnergy(ESD) - recalibrate cell: %s, module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n", + calo.Data(),module,icol,irow,frac,factor,cells->GetCellAmplitude(absId)); + + energy += cells->GetCellAmplitude(absId)*factor*frac; + } + + if(fDebug>1) + printf("AliCalorimeterUtils::RecalibrateClusterEnergy(ESD) - Energy before %f, after %f\n",cluster->E(),energy); + + return energy; + +} + + //________________________________________________________________ void AliCalorimeterUtils::SetGeometryTransformationMatrices(AliVEvent* inputEvent) { diff --git a/PWG4/PartCorrBase/AliCalorimeterUtils.h b/PWG4/PartCorrBase/AliCalorimeterUtils.h index 8f57336ecde..815f98d0882 100755 --- a/PWG4/PartCorrBase/AliCalorimeterUtils.h +++ b/PWG4/PartCorrBase/AliCalorimeterUtils.h @@ -115,24 +115,61 @@ public: void SwitchOnNoFiducialBorderInEMCALEta0() {fNoEMCALBorderAtEta0 = kTRUE; } void SwitchOffNoFiducialBorderInEMCALEta0() {fNoEMCALBorderAtEta0 = kFALSE; } + + // Recalibration + Bool_t IsRecalibrationOn() const { return fRecalibration ; } + void SwitchOnRecalibration() {fRecalibration = kTRUE ; InitEMCALRecalibrationFactors(); InitPHOSRecalibrationFactors();} + void SwitchOffRecalibration() {fRecalibration = kFALSE ; } + + void InitEMCALRecalibrationFactors() ; + void InitPHOSRecalibrationFactors () ; + + Float_t GetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow) const { + if(fEMCALRecalibrationFactors->GetEntries()>0) return (Float_t) ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->GetBinContent(iCol,iRow); + else return 1;} + + Float_t GetPHOSChannelRecalibrationFactor (Int_t imod, Int_t iCol, Int_t iRow) const { + if(fPHOSRecalibrationFactors->GetEntries()>0)return (Float_t) ((TH2F*)fPHOSRecalibrationFactors->At(imod))->GetBinContent(iCol,iRow); + else return 1;} + + void SetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow, Double_t c = 1) { + ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->SetBinContent(iCol,iRow,c);} + + void SetPHOSChannelRecalibrationFactor (Int_t imod, Int_t iCol, Int_t iRow, Double_t c = 1) { + ((TH2F*)fPHOSRecalibrationFactors->At(imod))->SetBinContent(iCol,iRow,c);} + + void SetEMCALChannelRecalibrationFactors(Int_t iSM , TH2F* h) {fEMCALRecalibrationFactors->AddAt(h,iSM);} + void SetPHOSChannelRecalibrationFactors(Int_t imod , TH2F* h) {fPHOSRecalibrationFactors ->AddAt(h,imod);} + + TH2F * GetEMCALChannelRecalibrationFactors(Int_t iSM) const {return (TH2F*)fEMCALRecalibrationFactors->At(iSM);} + TH2F * GetPHOSChannelRecalibrationFactors(Int_t imod) const {return (TH2F*)fPHOSRecalibrationFactors->At(imod);} + + void SetEMCALChannelRecalibrationFactors(TObjArray *map) {fEMCALRecalibrationFactors = map;} + void SetPHOSChannelRecalibrationFactors (TObjArray *map) {fPHOSRecalibrationFactors = map;} + + Float_t RecalibrateClusterEnergy(AliESDCaloCluster* cluster, AliESDCaloCells * cells); + Float_t RecalibrateClusterEnergy(AliAODCaloCluster* cluster, AliAODCaloCells * cells); private: - Int_t fDebug; // Debugging level - TString fEMCALGeoName; // Name of geometry to use for EMCAL. - TString fPHOSGeoName; // Name of geometry to use for PHOS. - AliEMCALGeoUtils * fEMCALGeo ; //! EMCAL geometry pointer - AliPHOSGeoUtils * fPHOSGeo ; //! PHOS geometry pointer - Bool_t fEMCALGeoMatrixSet; // Check if the transformation matrix is set for EMCAL - Bool_t fPHOSGeoMatrixSet ; // Check if the transformation matrix is set for PHOS - Bool_t fRemoveBadChannels; // Check the channel status provided and remove clusters with bad channels - TObjArray * fEMCALBadChannelMap; //! Array of histograms with map of bad channels, EMCAL - TObjArray * fPHOSBadChannelMap; //! Array of histograms with map of bad channels, PHOS - Int_t fNCellsFromEMCALBorder; // Number of cells from EMCAL border the cell with maximum amplitude has to be. - Int_t fNCellsFromPHOSBorder; // Number of cells from PHOS border the cell with maximum amplitude has to be. - Bool_t fNoEMCALBorderAtEta0;// Do fidutial cut in EMCAL region eta = 0? - - ClassDef(AliCalorimeterUtils,1) + Int_t fDebug; // Debugging level + TString fEMCALGeoName; // Name of geometry to use for EMCAL. + TString fPHOSGeoName; // Name of geometry to use for PHOS. + AliEMCALGeoUtils * fEMCALGeo ; //! EMCAL geometry pointer + AliPHOSGeoUtils * fPHOSGeo ; //! PHOS geometry pointer + Bool_t fEMCALGeoMatrixSet; // Check if the transformation matrix is set for EMCAL + Bool_t fPHOSGeoMatrixSet ; // Check if the transformation matrix is set for PHOS + Bool_t fRemoveBadChannels; // Check the channel status provided and remove clusters with bad channels + TObjArray * fEMCALBadChannelMap; //! Array of histograms with map of bad channels, EMCAL + TObjArray * fPHOSBadChannelMap; //! Array of histograms with map of bad channels, PHOS + Int_t fNCellsFromEMCALBorder; // Number of cells from EMCAL border the cell with maximum amplitude has to be. + Int_t fNCellsFromPHOSBorder; // Number of cells from PHOS border the cell with maximum amplitude has to be. + Bool_t fNoEMCALBorderAtEta0; // Do fiducial cut in EMCAL region eta = 0? + Bool_t fRecalibration; // Switch on or off the recalibration + TObjArray * fEMCALRecalibrationFactors; //! Array of histograms with map of recalibration factors, EMCAL + TObjArray * fPHOSRecalibrationFactors; //! Array of histograms with map of recalibration factors, PHOS + + ClassDef(AliCalorimeterUtils,2) } ; -- 2.39.3