X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALRecoUtils.cxx;h=705746eb818c95c5ce9ef6cd5cecac91c198ac87;hb=c1d7cf8d81b67cf2d2617d0f75bbd43fd43b0450;hp=58734f723d4d1eba6f3352809b6802dbabfbb102;hpb=841dbf607485c6a5f6d7810cdecd2cf1b477fb1c;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALRecoUtils.cxx b/EMCAL/AliEMCALRecoUtils.cxx index 58734f723d4..705746eb818 100644 --- a/EMCAL/AliEMCALRecoUtils.cxx +++ b/EMCAL/AliEMCALRecoUtils.cxx @@ -13,7 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* $Id: AliEMCALRecoUtils.cxx 33808 2009-07-15 09:48:08Z gconesab $ */ +/* $Id: AliEMCALRecoUtils.cxx | Sun Dec 8 06:56:48 2013 +0100 | Constantin Loizides $ */ /////////////////////////////////////////////////////////////////////////////// // @@ -56,10 +56,8 @@ #include "AliEMCALRecoUtils.h" #include "AliEMCALGeometry.h" #include "AliTrackerBase.h" -#include "AliEMCALCalibTimeDepCorrection.h" // Run dependent #include "AliEMCALPIDUtils.h" - ClassImp(AliEMCALRecoUtils) //_____________________________________ @@ -68,22 +66,24 @@ AliEMCALRecoUtils::AliEMCALRecoUtils(): fNonLinearityFunction(0), fNonLinearThreshold(0), fSmearClusterEnergy(kFALSE), fRandom(), fCellsRecalibrated(kFALSE), fRecalibration(kFALSE), fEMCALRecalibrationFactors(), - fTimeRecalibration(kFALSE), fEMCALTimeRecalibrationFactors(), - fUseRunCorrectionFactors(kFALSE), fRunCorrectionFactorsSet(kFALSE), + fTimeRecalibration(kFALSE), fEMCALTimeRecalibrationFactors(), fUseRunCorrectionFactors(kFALSE), fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(), fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE), fRejectExoticCluster(kFALSE), fRejectExoticCells(kFALSE), fExoticCellFraction(0), fExoticCellDiffTime(0), fExoticCellMinAmplitude(0), fPIDUtils(), fAODFilterMask(0), + fAODHybridTracks(0), fAODTPCOnlyTracks(0), fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0), fResidualEta(0x0), fResidualPhi(0x0), fCutEtaPhiSum(kFALSE), fCutEtaPhiSeparate(kFALSE), fCutR(0), fCutEta(0), fCutPhi(0), fClusterWindow(0), fMass(0), fStepSurface(0), fStepCluster(0), + fITSTrackSA(kFALSE), fEMCalSurfaceDistance(440.), fTrackCutsType(0), fCutMinTrackPt(0), fCutMinNClusterTPC(0), fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0), fCutRequireTPCRefit(kFALSE), fCutRequireITSRefit(kFALSE), fCutAcceptKinkDaughters(kFALSE), - fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0), fCutDCAToVertex2D(kFALSE) + fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0), fCutDCAToVertex2D(kFALSE), + fCutRequireITSStandAlone(kFALSE), fCutRequireITSpureSA(kFALSE) { // // Constructor. @@ -101,7 +101,6 @@ AliEMCALRecoUtils::AliEMCALRecoUtils(): fResidualEta = new TArrayF(); fPIDUtils = new AliEMCALPIDUtils(); - InitTrackCuts(); } //______________________________________________________________________ @@ -113,7 +112,7 @@ AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco) fCellsRecalibrated(reco.fCellsRecalibrated), fRecalibration(reco.fRecalibration), fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors), fTimeRecalibration(reco.fTimeRecalibration), fEMCALTimeRecalibrationFactors(reco.fEMCALTimeRecalibrationFactors), - fUseRunCorrectionFactors(reco.fUseRunCorrectionFactors), fRunCorrectionFactorsSet(reco.fRunCorrectionFactorsSet), + fUseRunCorrectionFactors(reco.fUseRunCorrectionFactors), fRemoveBadChannels(reco.fRemoveBadChannels), fRecalDistToBadChannels(reco.fRecalDistToBadChannels), fEMCALBadChannelMap(reco.fEMCALBadChannelMap), fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder), fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0), @@ -121,6 +120,7 @@ AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco) fExoticCellFraction(reco.fExoticCellFraction), fExoticCellDiffTime(reco.fExoticCellDiffTime), fExoticCellMinAmplitude(reco.fExoticCellMinAmplitude), fPIDUtils(reco.fPIDUtils), fAODFilterMask(reco.fAODFilterMask), + fAODHybridTracks(reco.fAODHybridTracks), fAODTPCOnlyTracks(reco.fAODTPCOnlyTracks), fMatchedTrackIndex( reco.fMatchedTrackIndex? new TArrayI(*reco.fMatchedTrackIndex):0x0), fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0), fResidualEta( reco.fResidualEta? new TArrayF(*reco.fResidualEta):0x0), @@ -129,19 +129,21 @@ AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco) fCutR(reco.fCutR), fCutEta(reco.fCutEta), fCutPhi(reco.fCutPhi), fClusterWindow(reco.fClusterWindow), fMass(reco.fMass), fStepSurface(reco.fStepSurface), fStepCluster(reco.fStepCluster), + fITSTrackSA(reco.fITSTrackSA), fEMCalSurfaceDistance(440.), fTrackCutsType(reco.fTrackCutsType), fCutMinTrackPt(reco.fCutMinTrackPt), fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS), fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS), fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit), fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters), fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY), - fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ), fCutDCAToVertex2D(reco.fCutDCAToVertex2D) + fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ), fCutDCAToVertex2D(reco.fCutDCAToVertex2D), + fCutRequireITSStandAlone(reco.fCutRequireITSStandAlone), fCutRequireITSpureSA(reco.fCutRequireITSpureSA) { //Copy ctor - for(Int_t i = 0; i < 15 ; i++) { fMisalRotShift[i] = reco.fMisalRotShift[i] ; - fMisalTransShift[i] = reco.fMisalTransShift[i] ; } - for(Int_t i = 0; i < 7 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; } - for(Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; } + for (Int_t i = 0; i < 15 ; i++) { fMisalRotShift[i] = reco.fMisalRotShift[i] ; + fMisalTransShift[i] = reco.fMisalTransShift[i] ; } + for (Int_t i = 0; i < 7 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; } + for (Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; } } @@ -151,13 +153,13 @@ AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & rec { //Assignment operator - if(this == &reco)return *this; + if (this == &reco)return *this; ((TNamed *)this)->operator=(reco); - for(Int_t i = 0; i < 15 ; i++) { fMisalTransShift[i] = reco.fMisalTransShift[i] ; + for (Int_t i = 0; i < 15 ; i++) { fMisalTransShift[i] = reco.fMisalTransShift[i] ; fMisalRotShift[i] = reco.fMisalRotShift[i] ; } - for(Int_t i = 0; i < 7 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; } - for(Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; } + for (Int_t i = 0; i < 7 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; } + for (Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; } fParticleType = reco.fParticleType; fPosAlgo = reco.fPosAlgo; @@ -175,7 +177,6 @@ AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & rec fEMCALTimeRecalibrationFactors = reco.fEMCALTimeRecalibrationFactors; fUseRunCorrectionFactors = reco.fUseRunCorrectionFactors; - fRunCorrectionFactorsSet = reco.fRunCorrectionFactorsSet; fRemoveBadChannels = reco.fRemoveBadChannels; fRecalDistToBadChannels = reco.fRecalDistToBadChannels; @@ -193,6 +194,8 @@ AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & rec fPIDUtils = reco.fPIDUtils; fAODFilterMask = reco.fAODFilterMask; + fAODHybridTracks = reco.fAODHybridTracks; + fAODTPCOnlyTracks = reco.fAODTPCOnlyTracks; fCutEtaPhiSum = reco.fCutEtaPhiSum; fCutEtaPhiSeparate = reco.fCutEtaPhiSeparate; @@ -203,7 +206,9 @@ AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & rec fMass = reco.fMass; fStepSurface = reco.fStepSurface; fStepCluster = reco.fStepCluster; - + fITSTrackSA = reco.fITSTrackSA; + fEMCalSurfaceDistance = reco.fEMCalSurfaceDistance; + fTrackCutsType = reco.fTrackCutsType; fCutMinTrackPt = reco.fCutMinTrackPt; fCutMinNClusterTPC = reco.fCutMinNClusterTPC; @@ -216,75 +221,75 @@ AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & rec fCutMaxDCAToVertexXY = reco.fCutMaxDCAToVertexXY; fCutMaxDCAToVertexZ = reco.fCutMaxDCAToVertexZ; fCutDCAToVertex2D = reco.fCutDCAToVertex2D; - - if(reco.fResidualEta){ + fCutRequireITSStandAlone = reco.fCutRequireITSStandAlone; + fCutRequireITSpureSA = reco.fCutRequireITSpureSA; + if (reco.fResidualEta) { // assign or copy construct - if(fResidualEta){ + if (fResidualEta) { *fResidualEta = *reco.fResidualEta; + } else { + fResidualEta = new TArrayF(*reco.fResidualEta); } - else fResidualEta = new TArrayF(*reco.fResidualEta); - } - else{ - if(fResidualEta)delete fResidualEta; + } else { + if (fResidualEta) delete fResidualEta; fResidualEta = 0; } - if(reco.fResidualPhi){ + if (reco.fResidualPhi) { // assign or copy construct - if(fResidualPhi){ + if (fResidualPhi) { *fResidualPhi = *reco.fResidualPhi; + } else { + fResidualPhi = new TArrayF(*reco.fResidualPhi); } - else fResidualPhi = new TArrayF(*reco.fResidualPhi); - } - else{ - if(fResidualPhi)delete fResidualPhi; + } else { + if (fResidualPhi) delete fResidualPhi; fResidualPhi = 0; } - if(reco.fMatchedTrackIndex){ + if (reco.fMatchedTrackIndex) { // assign or copy construct - if(fMatchedTrackIndex){ + if (fMatchedTrackIndex) { *fMatchedTrackIndex = *reco.fMatchedTrackIndex; + } else { + fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex); } - else fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex); - } - else{ - if(fMatchedTrackIndex)delete fMatchedTrackIndex; + } else { + if (fMatchedTrackIndex) delete fMatchedTrackIndex; fMatchedTrackIndex = 0; } - if(reco.fMatchedClusterIndex){ + if (reco.fMatchedClusterIndex) { // assign or copy construct - if(fMatchedClusterIndex){ + if (fMatchedClusterIndex) { *fMatchedClusterIndex = *reco.fMatchedClusterIndex; + } else { + fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex); } - else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex); - } - else{ - if(fMatchedClusterIndex)delete fMatchedClusterIndex; + } else { + if (fMatchedClusterIndex) delete fMatchedClusterIndex; fMatchedClusterIndex = 0; } return *this; } - //_____________________________________ AliEMCALRecoUtils::~AliEMCALRecoUtils() { //Destructor. - - if(fEMCALRecalibrationFactors) { + + if (fEMCALRecalibrationFactors) { fEMCALRecalibrationFactors->Clear(); delete fEMCALRecalibrationFactors; - } + } - if(fEMCALTimeRecalibrationFactors) { + if (fEMCALTimeRecalibrationFactors) { fEMCALTimeRecalibrationFactors->Clear(); delete fEMCALTimeRecalibrationFactors; - } + } - if(fEMCALBadChannelMap) { + if (fEMCALBadChannelMap) { fEMCALBadChannelMap->Clear(); delete fEMCALBadChannelMap; } @@ -299,7 +304,7 @@ AliEMCALRecoUtils::~AliEMCALRecoUtils() } //_______________________________________________________________________________ -Bool_t AliEMCALRecoUtils::AcceptCalibrateCell(const Int_t absID, const Int_t bc, +Bool_t AliEMCALRecoUtils::AcceptCalibrateCell(Int_t absID, Int_t bc, Float_t & amp, Double_t & time, AliVCaloCells* cells) { @@ -307,23 +312,29 @@ Bool_t AliEMCALRecoUtils::AcceptCalibrateCell(const Int_t absID, const Int_t bc, AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance(); - if(absID < 0 || absID >= 24*48*geom->GetNumberOfSuperModules()) return kFALSE; + if (absID < 0 || absID >= 24*48*geom->GetNumberOfSuperModules()) + return kFALSE; Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; - geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta); - geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta); + if (!geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta)) { + // cell absID does not exist + amp=0; time = 1.e9; + return kFALSE; + } + + geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta); + // Do not include bad channels found in analysis, - if( IsBadChannelsRemovalSwitchedOn() && GetEMCALChannelStatus(imod, ieta, iphi)) { + if (IsBadChannelsRemovalSwitchedOn() && GetEMCALChannelStatus(imod, ieta, iphi)) { return kFALSE; } //Recalibrate energy amp = cells->GetCellAmplitude(absID); - if(!fCellsRecalibrated && IsRecalibrationOn()) + if (!fCellsRecalibrated && IsRecalibrationOn()) amp *= GetEMCALChannelRecalibrationFactor(imod,ieta,iphi); - // Recalibrate time time = cells->GetCellTime(absID); @@ -332,162 +343,180 @@ Bool_t AliEMCALRecoUtils::AcceptCalibrateCell(const Int_t absID, const Int_t bc, return kTRUE; } -//_______________________________________________________________ -Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells) +//_____________________________________________________________________________ +Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(const AliEMCALGeometry* geom, + const AliVCluster* cluster, + AliVCaloCells* cells) { - // Given the list of AbsId of the cluster, get the maximum cell and - // check if there are fNCellsFromBorder from the calorimeter border - - if(!cluster){ + // Given the list of AbsId of the cluster, get the maximum cell and + // check if there are fNCellsFromBorder from the calorimeter border + + if (!cluster) + { AliInfo("Cluster pointer null!"); return kFALSE; } //If the distance to the border is 0 or negative just exit accept all clusters - if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE; + if (cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) + return kTRUE; - Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1; + Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1; Bool_t shared = kFALSE; GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared); - + AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n", - absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0)); - - if(absIdMax==-1) return kFALSE; - - //Check if the cell is close to the borders: - Bool_t okrow = kFALSE; - Bool_t okcol = kFALSE; - - if(iSM < 0 || iphi < 0 || ieta < 0 ) { + absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0)); + + if (absIdMax==-1) return kFALSE; + + //Check if the cell is close to the borders: + Bool_t okrow = kFALSE; + Bool_t okcol = kFALSE; + + if (iSM < 0 || iphi < 0 || ieta < 0 ) { AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n", iSM,ieta,iphi)); } //Check rows/phi - if(iSM < 10){ - if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE; - } - else{ - if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE; - } + Int_t iPhiLast = 24; + if( geom->GetSMType(iSM) == AliEMCALGeometry::kEMCAL_Half ) iPhiLast /= 2; + else if ( geom->GetSMType(iSM) == AliEMCALGeometry::kEMCAL_3rd ) iPhiLast /= 3;// 1/3 sm case + if(iphi >= fNCellsFromEMCALBorder && iphi < iPhiLast - fNCellsFromEMCALBorder) okrow = kTRUE; + //Check columns/eta - if(!fNoEMCALBorderAtEta0){ - if(ieta > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE; - } - else{ - if(iSM%2==0){ - if(ieta >= fNCellsFromEMCALBorder) okcol = kTRUE; - } - else { - if(ieta < 48-fNCellsFromEMCALBorder) okcol = kTRUE; + Int_t iEtaLast = 48; + if(!fNoEMCALBorderAtEta0 || geom->IsDCALSM(iSM)) {// conside inner border + if( geom->GetSMType(iSM) == AliEMCALGeometry::kDCAL_Standard ) iEtaLast = iEtaLast*2/3; + if(ieta > fNCellsFromEMCALBorder && ieta < iEtaLast-fNCellsFromEMCALBorder) okcol = kTRUE; + } else { + if (iSM%2==0) { + if (ieta >= fNCellsFromEMCALBorder) okcol = kTRUE; + } else { + if(ieta < iEtaLast-fNCellsFromEMCALBorder) okcol = kTRUE; } }//eta 0 not checked - + AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq", - fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow)); - - if (okcol && okrow) { + fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow)); + + if (okcol && okrow) { //printf("Accept\n"); return kTRUE; - } - else { + } else { //printf("Reject\n"); AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM)); return kFALSE; } - -} - +} -//_________________________________________________________________________________________________________ -Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(const AliEMCALGeometry* geom, const UShort_t* cellList, const Int_t nCells) +//_______________________________________________________________________________ +Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(const AliEMCALGeometry* geom, + const UShort_t* cellList, + Int_t nCells) { // Check that in the cluster cells, there is no bad channel of those stored // in fEMCALBadChannelMap or fPHOSBadChannelMap - - if(!fRemoveBadChannels) return kFALSE; - if(!fEMCALBadChannelMap) return kFALSE; - + + if (!fRemoveBadChannels) return kFALSE; + if (!fEMCALBadChannelMap) return kFALSE; + Int_t icol = -1; Int_t irow = -1; Int_t imod = -1; - for(Int_t iCell = 0; iCellGetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta); - if(fEMCALBadChannelMap->GetEntries() <= imod) continue; - geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol); - if(GetEMCALChannelStatus(imod, icol, irow)){ + if (fEMCALBadChannelMap->GetEntries() <= imod) continue; + geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol); + if (GetEMCALChannelStatus(imod, icol, irow)) { AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow)); return kTRUE; } - }// cell cluster loop - + return kFALSE; } -//_______________________________________________________________________ -Bool_t AliEMCALRecoUtils::IsExoticCell(const Int_t absID, AliVCaloCells* cells, const Int_t bc) -{ - // Look to cell neighbourhood and reject if it seems exotic - // Do before recalibrating the cells - - if(!fRejectExoticCells) return kFALSE; +//___________________________________________________________________________ +Float_t AliEMCALRecoUtils::GetECross(Int_t absID, Double_t tcell, + AliVCaloCells* cells, Int_t bc) +{ + //Calculate the energy in the cross around the energy given cell + AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance(); Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta); - geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta); + geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta); //Get close cells index, energy and time, not in corners - Int_t absID1 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta); - Int_t absID2 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta); - Int_t absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1); - Int_t absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1); - Float_t ecell = 0, ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0; - Double_t tcell = 0, tcell1 = 0, tcell2 = 0, tcell3 = 0, tcell4 = 0; - Bool_t accept = 0, accept1 = 0, accept2 = 0, accept3 = 0, accept4 = 0; + Int_t absID1 = -1; + Int_t absID2 = -1; - accept = AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells); - - if(!accept) return kTRUE; // reject this cell - - if(ecell < fExoticCellMinAmplitude) return kFALSE; // do not reject low energy cells - - accept1 = AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells); - accept2 = AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells); - accept3 = AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells); - accept4 = AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells); - - /* - printf("Cell absID %d \n",absID); - printf("\t accept1 %d, accept2 %d, accept3 %d, accept4 %d\n", - accept1,accept2,accept3,accept4); - printf("\t id %d: id1 %d, id2 %d, id3 %d, id4 %d\n", - absID,absID1,absID2,absID3,absID4); - printf("\t e %f: e1 %f, e2 %f, e3 %f, e4 %f\n", - ecell,ecell1,ecell2,ecell3,ecell4); - printf("\t t %f: t1 %f, t2 %f, t3 %f, t4 %f;\n dt1 %f, dt2 %f, dt3 %f, dt4 %f\n", - tcell*1.e9,tcell1*1.e9,tcell2*1.e9,tcell3*1.e9,tcell4*1.e9, - TMath::Abs(tcell-tcell1)*1.e9, TMath::Abs(tcell-tcell2)*1.e9, TMath::Abs(tcell-tcell3)*1.e9, TMath::Abs(tcell-tcell4)*1.e9); - */ - - if(TMath::Abs(tcell-tcell1)*1.e9 > fExoticCellDiffTime) ecell1 = 0 ; - if(TMath::Abs(tcell-tcell2)*1.e9 > fExoticCellDiffTime) ecell2 = 0 ; - if(TMath::Abs(tcell-tcell3)*1.e9 > fExoticCellDiffTime) ecell3 = 0 ; - if(TMath::Abs(tcell-tcell4)*1.e9 > fExoticCellDiffTime) ecell4 = 0 ; - - Float_t eCross = ecell1+ecell2+ecell3+ecell4; - - //printf("\t eCell %f, eCross %f, 1-eCross/eCell %f\n",ecell,eCross,1-eCross/ecell); - - if(1-eCross/ecell > fExoticCellFraction) { + if ( iphi < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta); + if ( iphi > 0 ) absID2 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta); + + // In case of cell in eta = 0 border, depending on SM shift the cross cell index + + Int_t absID3 = -1; + Int_t absID4 = -1; + + if ( ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2) ) { + absID3 = geom-> GetAbsCellIdFromCellIndexes(imod+1, iphi, 0); + absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1); + } else if ( ieta == 0 && imod%2 ) { + absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1); + absID4 = geom-> GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1); + } else { + if ( ieta < AliEMCALGeoParams::fgkEMCALCols-1 ) + absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1); + if ( ieta > 0 ) + absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1); + } + + //printf("IMOD %d, AbsId %d, a %d, b %d, c %d e %d \n",imod,absID,absID1,absID2,absID3,absID4); + + Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0; + Double_t tcell1 = 0, tcell2 = 0, tcell3 = 0, tcell4 = 0; + + AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells); + AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells); + AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells); + AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells); + + if (TMath::Abs(tcell-tcell1)*1.e9 > fExoticCellDiffTime) ecell1 = 0 ; + if (TMath::Abs(tcell-tcell2)*1.e9 > fExoticCellDiffTime) ecell2 = 0 ; + if (TMath::Abs(tcell-tcell3)*1.e9 > fExoticCellDiffTime) ecell3 = 0 ; + if (TMath::Abs(tcell-tcell4)*1.e9 > fExoticCellDiffTime) ecell4 = 0 ; + + return ecell1+ecell2+ecell3+ecell4; +} + +//_____________________________________________________________________________________________ +Bool_t AliEMCALRecoUtils::IsExoticCell(Int_t absID, AliVCaloCells* cells, Int_t bc) +{ + // Look to cell neighbourhood and reject if it seems exotic + // Do before recalibrating the cells + + if (!fRejectExoticCells) return kFALSE; + + Float_t ecell = 0; + Double_t tcell = 0; + Bool_t accept = AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells); + + if (!accept) return kTRUE; // reject this cell + + if (ecell < fExoticCellMinAmplitude) return kFALSE; // do not reject low energy cells + + Float_t eCross = GetECross(absID,tcell,cells,bc); + + if (1-eCross/ecell > fExoticCellFraction) { AliDebug(2,Form("AliEMCALRecoUtils::IsExoticCell() - EXOTIC CELL id %d, eCell %f, eCross %f, 1-eCross/eCell %f\n", absID,ecell,eCross,1-eCross/ecell)); return kTRUE; @@ -496,40 +525,42 @@ Bool_t AliEMCALRecoUtils::IsExoticCell(const Int_t absID, AliVCaloCells* cells, return kFALSE; } -//_________________________________________________ -Bool_t AliEMCALRecoUtils::IsExoticCluster(AliVCluster *cluster, AliVCaloCells *cells, const Int_t bc) +//___________________________________________________________________ +Bool_t AliEMCALRecoUtils::IsExoticCluster(const AliVCluster *cluster, + AliVCaloCells *cells, + Int_t bc) { // Check if the cluster highest energy tower is exotic - if(!cluster){ + if (!cluster) { AliInfo("Cluster pointer null!"); return kFALSE; } - if(!fRejectExoticCluster) return kFALSE; + if (!fRejectExoticCluster) return kFALSE; // Get highest energy tower AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance(); Int_t iSupMod = -1, absId = -1, ieta = -1, iphi = -1; Bool_t shared = kFALSE; GetMaxEnergyCell(geom, cells, cluster, absId, iSupMod, ieta, iphi, shared); - + return IsExoticCell(absId,cells,bc); } -//__________________________________________________ +//_______________________________________________________________________ Float_t AliEMCALRecoUtils::SmearClusterEnergy(const AliVCluster* cluster) { //In case of MC analysis, smear energy to match resolution/calibration in real data - if(!cluster){ + if (!cluster) { AliInfo("Cluster pointer null!"); return 0; } Float_t energy = cluster->E() ; Float_t rdmEnergy = energy ; - if(fSmearClusterEnergy){ + if (fSmearClusterEnergy) { rdmEnergy = fRandom.Gaus(energy,fSmearClusterParam[0] * TMath::Sqrt(energy) + fSmearClusterParam[1] * energy + fSmearClusterParam[2] ); @@ -539,40 +570,73 @@ Float_t AliEMCALRecoUtils::SmearClusterEnergy(const AliVCluster* cluster) return rdmEnergy; } -//__________________________________________________ +//____________________________________________________________________________ Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster) { // Correct cluster energy from non linearity functions - if(!cluster){ + if (!cluster) { AliInfo("Cluster pointer null!"); return 0; } Float_t energy = cluster->E(); - switch (fNonLinearityFunction) { - + if (energy < 0.05) { + // Clusters with less than 50 MeV or negative are not possible + AliInfo(Form("Too Low Cluster energy!, E = %f < 0.05 GeV",energy)); + return 0; + } + + switch (fNonLinearityFunction) + { case kPi0MC: { //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2))))) - //Double_t fNonLinearityParams[0] = 1.014; - //Double_t fNonLinearityParams[1] = -0.03329; - //Double_t fNonLinearityParams[2] = -0.3853; - //Double_t fNonLinearityParams[3] = 0.5423; - //Double_t fNonLinearityParams[4] = -0.4335; + //fNonLinearityParams[0] = 1.014; + //fNonLinearityParams[1] =-0.03329; + //fNonLinearityParams[2] =-0.3853; + //fNonLinearityParams[3] = 0.5423; + //fNonLinearityParams[4] =-0.4335; energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+ ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())* exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3])))); break; } + case kPi0MCv2: + { + //Non-Linearity correction (from MC with function [0]/((x+[1])^[2]))+1; + //fNonLinearityParams[0] = 3.11111e-02; + //fNonLinearityParams[1] =-5.71666e-02; + //fNonLinearityParams[2] = 5.67995e-01; + + energy *= fNonLinearityParams[0]/TMath::Power(energy+fNonLinearityParams[1],fNonLinearityParams[2])+1; + break; + } + + case kPi0MCv3: + { + //Same as beam test corrected, change parameters + //fNonLinearityParams[0] = 9.81039e-01 + //fNonLinearityParams[1] = 1.13508e-01; + //fNonLinearityParams[2] = 1.00173e+00; + //fNonLinearityParams[3] = 9.67998e-02; + //fNonLinearityParams[4] = 2.19381e+02; + //fNonLinearityParams[5] = 6.31604e+01; + //fNonLinearityParams[6] = 1; + energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5])))); + + break; + } + + case kPi0GammaGamma: { //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E)) - //Double_t fNonLinearityParams[0] = 1.04; - //Double_t fNonLinearityParams[1] = -0.1445; - //Double_t fNonLinearityParams[2] = 1.046; + //fNonLinearityParams[0] = 1.04; + //fNonLinearityParams[1] = -0.1445; + //fNonLinearityParams[2] = 1.046; energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function break; } @@ -615,6 +679,21 @@ Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster) break; } + + case kBeamTestCorrectedv2: + { + //From beam test, corrected for material between beam and EMCAL + //fNonLinearityParams[0] = 0.983504; + //fNonLinearityParams[1] = 0.210106; + //fNonLinearityParams[2] = 0.897274; + //fNonLinearityParams[3] = 0.0829064; + //fNonLinearityParams[4] = 152.299; + //fNonLinearityParams[5] = 31.5028; + //fNonLinearityParams[6] = 0.968; + energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5])))); + + break; + } case kNoCorrection: AliDebug(2,"No correction on the energy\n"); @@ -629,8 +708,8 @@ Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster) void AliEMCALRecoUtils::InitNonLinearityParam() { //Initialising Non Linearity Parameters - - if(fNonLinearityFunction == kPi0MC) { + + if (fNonLinearityFunction == kPi0MC) { fNonLinearityParams[0] = 1.014; fNonLinearityParams[1] = -0.03329; fNonLinearityParams[2] = -0.3853; @@ -638,37 +717,53 @@ void AliEMCALRecoUtils::InitNonLinearityParam() fNonLinearityParams[4] = -0.4335; } - if(fNonLinearityFunction == kPi0GammaGamma) { + if (fNonLinearityFunction == kPi0MCv2) { + fNonLinearityParams[0] = 3.11111e-02; + fNonLinearityParams[1] =-5.71666e-02; + fNonLinearityParams[2] = 5.67995e-01; + } + + if (fNonLinearityFunction == kPi0MCv3) { + fNonLinearityParams[0] = 9.81039e-01; + fNonLinearityParams[1] = 1.13508e-01; + fNonLinearityParams[2] = 1.00173e+00; + fNonLinearityParams[3] = 9.67998e-02; + fNonLinearityParams[4] = 2.19381e+02; + fNonLinearityParams[5] = 6.31604e+01; + fNonLinearityParams[6] = 1; + } + + if (fNonLinearityFunction == kPi0GammaGamma) { fNonLinearityParams[0] = 1.04; fNonLinearityParams[1] = -0.1445; fNonLinearityParams[2] = 1.046; - } + } - if(fNonLinearityFunction == kPi0GammaConversion) { + if (fNonLinearityFunction == kPi0GammaConversion) { fNonLinearityParams[0] = 0.139393; fNonLinearityParams[1] = 0.0566186; fNonLinearityParams[2] = 0.982133; - } + } - if(fNonLinearityFunction == kBeamTest) { - if(fNonLinearThreshold == 30) { + if (fNonLinearityFunction == kBeamTest) { + if (fNonLinearThreshold == 30) { fNonLinearityParams[0] = 1.007; fNonLinearityParams[1] = 0.894; fNonLinearityParams[2] = 0.246; } - if(fNonLinearThreshold == 45) { + if (fNonLinearThreshold == 45) { fNonLinearityParams[0] = 1.003; fNonLinearityParams[1] = 0.719; fNonLinearityParams[2] = 0.334; } - if(fNonLinearThreshold == 75) { + if (fNonLinearThreshold == 75) { fNonLinearityParams[0] = 1.002; fNonLinearityParams[1] = 0.797; fNonLinearityParams[2] = 0.358; } } - if(fNonLinearityFunction == kBeamTestCorrected) { + if (fNonLinearityFunction == kBeamTestCorrected) { fNonLinearityParams[0] = 0.99078; fNonLinearityParams[1] = 0.161499; fNonLinearityParams[2] = 0.655166; @@ -677,10 +772,22 @@ void AliEMCALRecoUtils::InitNonLinearityParam() fNonLinearityParams[5] = 23.6904; fNonLinearityParams[6] = 0.978; } + + if (fNonLinearityFunction == kBeamTestCorrectedv2) { + fNonLinearityParams[0] = 0.983504; + fNonLinearityParams[1] = 0.210106; + fNonLinearityParams[2] = 0.897274; + fNonLinearityParams[3] = 0.0829064; + fNonLinearityParams[4] = 152.299; + fNonLinearityParams[5] = 31.5028; + fNonLinearityParams[6] = 0.968; + } } -//__________________________________________________ -Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const +//_________________________________________________________ +Float_t AliEMCALRecoUtils::GetDepth(Float_t energy, + Int_t iParticle, + Int_t iSM) const { //Calculate shower depth for a given cluster energy and particle type @@ -688,40 +795,55 @@ Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle Float_t x0 = 1.31; Float_t ecr = 8; Float_t depth = 0; + Float_t arg = energy*1000/ ecr; //Multiply energy by 1000 to transform to MeV switch ( iParticle ) { case kPhoton: - depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV + if (arg < 1) + depth = 0; + else + depth = x0 * (TMath::Log(arg) + 0.5); break; case kElectron: - depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV + if (arg < 1) + depth = 0; + else + depth = x0 * (TMath::Log(arg) - 0.5); break; case kHadron: // hadron // boxes anc. here - if(gGeoManager){ + if (gGeoManager) { gGeoManager->cd("ALIC_1/XEN1_1"); TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode(); TGeoNodeMatrix *geoSM = dynamic_cast(geoXEn1->GetDaughter(iSM)); - if(geoSM){ + if (geoSM) { TGeoVolume *geoSMVol = geoSM->GetVolume(); TGeoShape *geoSMShape = geoSMVol->GetShape(); TGeoBBox *geoBox = dynamic_cast(geoSMShape); - if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ; + if (geoBox) depth = 0.5 * geoBox->GetDX()*2 ; else AliFatal("Null GEANT box"); - }else AliFatal("NULL GEANT node matrix"); + } + else AliFatal("NULL GEANT node matrix"); } - else{//electron - depth = x0 * (TMath::Log(energy*1000 / ecr) - 0.5); //Multiply energy by 1000 to transform to MeV + else + {//electron + if (arg < 1) + depth = 0; + else + depth = x0 * (TMath::Log(arg) - 0.5); } break; default://photon - depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV + if (arg < 1) + depth = 0; + else + depth = x0 * (TMath::Log(arg) + 0.5); } return depth; @@ -731,11 +853,11 @@ Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle void AliEMCALRecoUtils::GetMaxEnergyCell(const AliEMCALGeometry *geom, AliVCaloCells* cells, const AliVCluster* clu, - Int_t & absId, - Int_t& iSupMod, - Int_t& ieta, - Int_t& iphi, - Bool_t &shared) + Int_t & absId, + Int_t & iSupMod, + Int_t & ieta, + Int_t & iphi, + Bool_t & shared) { //For a given CaloCluster gets the absId of the cell //with maximum energy deposit. @@ -751,7 +873,7 @@ void AliEMCALRecoUtils::GetMaxEnergyCell(const AliEMCALGeometry *geom, Int_t iIeta = -1; Int_t iSupMod0= -1; - if(!clu){ + if (!clu) { AliInfo("Cluster pointer null!"); absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1; return; @@ -761,20 +883,21 @@ void AliEMCALRecoUtils::GetMaxEnergyCell(const AliEMCALGeometry *geom, cellAbsId = clu->GetCellAbsId(iDig); fraction = clu->GetCellAmplitudeFraction(iDig); //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction); - if(fraction < 1e-4) fraction = 1.; // in case unfolding is off + if (fraction < 1e-4) fraction = 1.; // in case unfolding is off geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta); geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta); - if(iDig==0) iSupMod0=iSupMod; - else if(iSupMod0!=iSupMod) { + if (iDig==0) { + iSupMod0=iSupMod; + } else if (iSupMod0!=iSupMod) { shared = kTRUE; //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n"); } - if(!fCellsRecalibrated && IsRecalibrationOn()) { + if (!fCellsRecalibrated && IsRecalibrationOn()) { recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi); } eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor; //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction); - if(eCell > eMax) { + if (eCell > eMax) { eMax = eCell; absId = cellAbsId; //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell); @@ -806,7 +929,9 @@ void AliEMCALRecoUtils::InitParameters() fExoticCellDiffTime = 1e6; fExoticCellMinAmplitude = 0.5; - fAODFilterMask = 32; + fAODFilterMask = 128; + fAODHybridTracks = kFALSE; + fAODTPCOnlyTracks = kTRUE; fCutEtaPhiSum = kTRUE; fCutEtaPhiSeparate = kFALSE; @@ -837,29 +962,27 @@ void AliEMCALRecoUtils::InitParameters() fCutMaxDCAToVertexZ = 1e10; fCutDCAToVertex2D = kFALSE; + fCutRequireITSStandAlone = kFALSE; //MARCEL + fCutRequireITSpureSA = kFALSE; //Marcel //Misalignment matrices - for(Int_t i = 0; i < 15 ; i++) { + for (Int_t i = 0; i < 15 ; i++) + { fMisalTransShift[i] = 0.; fMisalRotShift[i] = 0.; } //Non linearity - for(Int_t i = 0; i < 7 ; i++) fNonLinearityParams[i] = 0.; - - //For kBeamTestCorrected case, but default is no correction - fNonLinearityParams[0] = 0.99078; - fNonLinearityParams[1] = 0.161499; - fNonLinearityParams[2] = 0.655166; - fNonLinearityParams[3] = 0.134101; - fNonLinearityParams[4] = 163.282; - fNonLinearityParams[5] = 23.6904; - fNonLinearityParams[6] = 0.978; - - //For kPi0GammaGamma case - //fNonLinearityParams[0] = 0.1457/0.1349766/1.038; - //fNonLinearityParams[1] = -0.02024/0.1349766/1.038; - //fNonLinearityParams[2] = 1.046; + for (Int_t i = 0; i < 7 ; i++) fNonLinearityParams[i] = 0.; + + //For kBeamTestCorrectedv2 case, but default is no correction + fNonLinearityParams[0] = 0.983504; + fNonLinearityParams[1] = 0.210106; + fNonLinearityParams[2] = 0.897274; + fNonLinearityParams[3] = 0.0829064; + fNonLinearityParams[4] = 152.299; + fNonLinearityParams[5] = 31.5028; + fNonLinearityParams[6] = 0.968; //Cluster energy smearing fSmearClusterEnergy = kFALSE; @@ -873,30 +996,34 @@ void AliEMCALRecoUtils::InitEMCALRecalibrationFactors() { //Init EMCAL recalibration factors AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()"); - //In order to avoid rewriting the same histograms + //In order to avoid rewriting the same histograms Bool_t oldStatus = TH1::AddDirectoryStatus(); TH1::AddDirectory(kFALSE); - fEMCALRecalibrationFactors = new TObjArray(10); - for (int i = 0; i < 10; i++) + fEMCALRecalibrationFactors = new TObjArray(12); + 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 < 10; sm++) { - for (Int_t i = 0; i < 48; i++) { - for (Int_t j = 0; j < 24; j++) { + 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); + TH1::AddDirectory(oldStatus); } -//________________________________________________________________ +//_________________________________________________________ void AliEMCALRecoUtils::InitEMCALTimeRecalibrationFactors() { //Init EMCAL recalibration factors @@ -909,10 +1036,11 @@ void AliEMCALRecoUtils::InitEMCALTimeRecalibrationFactors() for (int i = 0; i < 4; i++) fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvBC%d",i), Form("hAllTimeAvBC%d",i), - 48*24*10,0.,48*24*10) ); + 48*24*12,0.,48*24*12) ); //Init the histograms with 1 - for (Int_t bc = 0; bc < 4; bc++) { - for (Int_t i = 0; i < 48*24*10; i++) + for (Int_t bc = 0; bc < 4; bc++) + { + for (Int_t i = 0; i < 48*24*12; i++) SetEMCALChannelTimeRecalibrationFactor(bc,i,0.); } @@ -920,10 +1048,10 @@ void AliEMCALRecoUtils::InitEMCALTimeRecalibrationFactors() fEMCALTimeRecalibrationFactors->Compress(); //In order to avoid rewriting the same histograms - TH1::AddDirectory(oldStatus); + TH1::AddDirectory(oldStatus); } -//________________________________________________________________ +//____________________________________________________ void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap() { //Init EMCAL bad channels map @@ -932,9 +1060,10 @@ void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap() Bool_t oldStatus = TH1::AddDirectoryStatus(); TH1::AddDirectory(kFALSE); - fEMCALBadChannelMap = new TObjArray(10); + fEMCALBadChannelMap = new TObjArray(12); //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24); - for (int i = 0; i < 10; i++) { + for (int i = 0; i < 12; i++) + { fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24)); } @@ -942,20 +1071,20 @@ void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap() fEMCALBadChannelMap->Compress(); //In order to avoid rewriting the same histograms - TH1::AddDirectory(oldStatus); + TH1::AddDirectory(oldStatus); } //____________________________________________________________________________ void AliEMCALRecoUtils::RecalibrateClusterEnergy(const AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells, - const Int_t bc) + Int_t bc) { // Recalibrate the cluster energy and Time, considering the recalibration map // and the energy of the cells and time that compose the cluster. // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber(); - - if(!cluster){ + + if (!cluster) { AliInfo("Cluster pointer null!"); return; } @@ -964,7 +1093,7 @@ void AliEMCALRecoUtils::RecalibrateClusterEnergy(const AliEMCALGeometry* geom, UShort_t * index = cluster->GetCellsAbsId() ; Double_t * fraction = cluster->GetCellsAmplitudeFraction() ; Int_t ncells = cluster->GetNCells(); - + //Initialize some used variables Float_t energy = 0; Int_t absId =-1; @@ -974,17 +1103,19 @@ void AliEMCALRecoUtils::RecalibrateClusterEnergy(const AliEMCALGeometry* geom, Float_t emax = 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++){ + for (Int_t icell = 0; icell < ncells; icell++) + { absId = index[icell]; frac = fraction[icell]; - if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off + if (frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off - if(!fCellsRecalibrated && IsRecalibrationOn()) { + if (!fCellsRecalibrated && IsRecalibrationOn()) { // Energy Int_t iTower = -1, iIphi = -1, iIeta = -1; geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta); - if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue; - geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol); + if (fEMCALRecalibrationFactors->GetEntries() <= imod) + continue; + geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol); factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow); AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n", @@ -994,124 +1125,103 @@ void AliEMCALRecoUtils::RecalibrateClusterEnergy(const AliEMCALGeometry* geom, energy += cells->GetCellAmplitude(absId)*factor*frac; - if(emax < cells->GetCellAmplitude(absId)*factor*frac){ + if (emax < cells->GetCellAmplitude(absId)*factor*frac) { emax = cells->GetCellAmplitude(absId)*factor*frac; absIdMax = absId; } } - - cluster->SetE(energy); + + AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f \n",cluster->E(),energy)); - AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy)); + cluster->SetE(energy); - // Recalculate time of cluster only for ESDs - if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName()))) { + // Recalculate time of cluster + Double_t timeorg = cluster->GetTOF(); - // Time - Double_t weightedTime = 0; - Double_t weight = 0; - Double_t weightTot = 0; - Double_t maxcellTime = 0; - for(Int_t icell = 0; icell < ncells; icell++){ - absId = index[icell]; - frac = fraction[icell]; - if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off + Double_t time = cells->GetCellTime(absIdMax); + if (!fCellsRecalibrated && IsTimeRecalibrationOn()) + RecalibrateCellTime(absIdMax,bc,time); - Double_t celltime = cells->GetCellTime(absId); - RecalibrateCellTime(absId, bc, celltime); - if(absId == absIdMax) maxcellTime = celltime; + cluster->SetTOF(time); - if(!fCellsRecalibrated){ - - Int_t iTower = -1, iIphi = -1, iIeta = -1; - geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta); - if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue; - geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol); - factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow); - - AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell:" - " module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n", - imod,icol,irow,frac,factor,cells->GetCellTime(absId))); - - } - - weight = GetCellWeight(cells->GetCellAmplitude(absId)*factor*frac , energy ); - weightTot += weight; - weightedTime += celltime * weight; - } - - if(weightTot > 0) - cluster->SetTOF(weightedTime/weightTot); - else - cluster->SetTOF(maxcellTime); - } + AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Time before %f, after %f \n",timeorg,cluster->GetTOF())); } -//________________________________________________________________ -void AliEMCALRecoUtils::RecalibrateCells(AliVCaloCells * cells, Int_t bc) +//_____________________________________________________________ +void AliEMCALRecoUtils::RecalibrateCells(AliVCaloCells * cells, + Int_t bc) { // Recalibrate the cells time and energy, considering the recalibration map and the energy // of the cells that compose the cluster. // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber(); - if(!IsRecalibrationOn() && !IsTimeRecalibrationOn()) return; + if (!IsRecalibrationOn() && !IsTimeRecalibrationOn() && !IsBadChannelsRemovalSwitchedOn()) + return; - if(!cells){ + if (!cells) { AliInfo("Cells pointer null!"); return; } - Int_t absId =-1; + Short_t absId =-1; Bool_t accept = kFALSE; Float_t ecell = 0; Double_t tcell = 0; + Double_t ecellin = 0; + Double_t tcellin = 0; + Int_t mclabel = -1; + Double_t efrac = 0; Int_t nEMcell = cells->GetNumberOfCells() ; - for (Int_t iCell = 0; iCell < nEMcell; iCell++) { - - absId = cells->GetCellNumber(iCell); + for (Int_t iCell = 0; iCell < nEMcell; iCell++) + { + cells->GetCell( iCell, absId, ecellin, tcellin, mclabel, efrac ); accept = AcceptCalibrateCell(absId, bc, ecell ,tcell ,cells); - if(!accept) { + if (!accept) { ecell = 0; - tcell = 0; + tcell = -1; } //Set new values - cells->SetCell(iCell,absId,ecell, tcell); + cells->SetCell(iCell,absId,ecell, tcell, mclabel, efrac); } fCellsRecalibrated = kTRUE; } //_______________________________________________________________________________________________________ -void AliEMCALRecoUtils::RecalibrateCellTime(const Int_t absId, const Int_t bc, Double_t & celltime) const +void AliEMCALRecoUtils::RecalibrateCellTime(Int_t absId, Int_t bc, Double_t & celltime) const { // Recalibrate time of cell with absID considering the recalibration map // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber(); - if(!fCellsRecalibrated && IsTimeRecalibrationOn() && bc >= 0){ + if (!fCellsRecalibrated && IsTimeRecalibrationOn() && bc >= 0) { celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId)*1.e-9; ; } } -//__________________________________________________ -void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu) +//______________________________________________________________________________ +void AliEMCALRecoUtils::RecalculateClusterPosition(const AliEMCALGeometry *geom, + AliVCaloCells* cells, + AliVCluster* clu) { //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster. - if(!clu){ + if (!clu) { AliInfo("Cluster pointer null!"); return; } - if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu); - else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu); - else AliDebug(2,"Algorithm to recalculate position not selected, do nothing."); + if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu); + else if (fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu); + else AliDebug(2,"Algorithm to recalculate position not selected, do nothing."); } -//__________________________________________________ -void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu) +//_____________________________________________________________________________________________ +void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(const AliEMCALGeometry *geom, + AliVCaloCells* cells, + AliVCluster* clu) { // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster. // The algorithm is a copy of what is done in AliEMCALRecPoint @@ -1129,22 +1239,23 @@ void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeomet Bool_t shared = kFALSE; Float_t clEnergy = clu->E(); //Energy already recalibrated previously + if (clEnergy <= 0) + return; GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared); Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ; //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth); - for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) { - + for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) + { absId = clu->GetCellAbsId(iDig); fraction = clu->GetCellAmplitudeFraction(iDig); - if(fraction < 1e-4) fraction = 1.; // in case unfolding is off + if (fraction < 1e-4) fraction = 1.; // in case unfolding is off - if (!fCellsRecalibrated){ + if (!fCellsRecalibrated) { geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta); - geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta); - - if(IsRecalibrationOn()) { + geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta); + if (IsRecalibrationOn()) { recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi); } } @@ -1159,11 +1270,11 @@ void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeomet geom->GetGlobal(pLocal,pGlobal,iSupModMax); //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]); - for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]); + for (int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]); }// cell loop - if(totalWeight>0){ - for(int i=0; i<3; i++ ) newPos[i] /= totalWeight; + if (totalWeight>0) { + for (int i=0; i<3; i++ ) newPos[i] /= totalWeight; } //Float_t pos[]={0,0,0}; @@ -1171,12 +1282,12 @@ void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeomet //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]); //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]); - if(iSupModMax > 1) {//sector 1 + if (iSupModMax > 1) { //sector 1 newPos[0] +=fMisalTransShift[3];//-=3.093; newPos[1] +=fMisalTransShift[4];//+=6.82; newPos[2] +=fMisalTransShift[5];//+=1.635; //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]); - } else {//sector 0 + } else { //sector 0 newPos[0] +=fMisalTransShift[0];//+=1.134; newPos[1] +=fMisalTransShift[1];//+=8.2; newPos[2] +=fMisalTransShift[2];//+=1.197; @@ -1187,8 +1298,10 @@ void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeomet clu->SetPosition(newPos); } -//__________________________________________________ -void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu) +//____________________________________________________________________________________________ +void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(const AliEMCALGeometry *geom, + AliVCaloCells* cells, + AliVCluster* clu) { // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster. // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position @@ -1205,6 +1318,9 @@ void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometr Bool_t shared = kFALSE; Float_t clEnergy = clu->E(); //Energy already recalibrated previously. + + if (clEnergy <= 0) + return; GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared); Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ; @@ -1212,21 +1328,23 @@ void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometr Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now Int_t startingSM = -1; - for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) { + for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) + { absId = clu->GetCellAbsId(iDig); fraction = clu->GetCellAmplitudeFraction(iDig); - if(fraction < 1e-4) fraction = 1.; // in case unfolding is off + if (fraction < 1e-4) fraction = 1.; // in case unfolding is off - if (iDig==0) startingSM = iSupMod; - else if(iSupMod != startingSM) areInSameSM = kFALSE; + if (iDig==0) startingSM = iSupMod; + else if (iSupMod != startingSM) areInSameSM = kFALSE; eCell = cells->GetCellAmplitude(absId); geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta); - geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta); + geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta); - if (!fCellsRecalibrated){ - if(IsRecalibrationOn()) { + if (!fCellsRecalibrated) + { + if (IsRecalibrationOn()) { recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi); } } @@ -1234,7 +1352,7 @@ void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometr eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor; weight = GetCellWeight(eCell,clEnergy); - if(weight < 0) weight = 0; + if (weight < 0) weight = 0; totalWeight += weight; weightedCol += ieta*weight; weightedRow += iphi*weight; @@ -1243,12 +1361,14 @@ void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometr }// cell loop Float_t xyzNew[]={0.,0.,0.}; - if(areInSameSM == kTRUE) { + if (areInSameSM == kTRUE) { //printf("In Same SM\n"); weightedCol = weightedCol/totalWeight; weightedRow = weightedRow/totalWeight; geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew); - } else { + } + else + { //printf("In Different SM\n"); geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew); } @@ -1256,40 +1376,46 @@ void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometr clu->SetPosition(xyzNew); } -//____________________________________________________________________________ -void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster) +//___________________________________________________________________________________________ +void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(const AliEMCALGeometry * geom, + AliVCaloCells* cells, + AliVCluster * cluster) { //re-evaluate distance to bad channel with updated bad map - if(!fRecalDistToBadChannels) return; + if (!fRecalDistToBadChannels) return; - if(!cluster){ + if (!cluster) + { AliInfo("Cluster pointer null!"); return; } //Get channels map of the supermodule where the cluster is. - Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1; + Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1; Bool_t shared = kFALSE; GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared); TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod); - Int_t dRrow, dRcol; + Int_t dRrow, dRcol; Float_t minDist = 10000.; Float_t dist = 0.; //Loop on tower status map - for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){ - for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){ + for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++) + { + for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++) + { //Check if tower is bad. - if(hMap->GetBinContent(icol,irow)==0) continue; + if (hMap->GetBinContent(icol,irow)==0) continue; //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n", // iSupMod,icol, irow, icolM,irowM); dRrow=TMath::Abs(irowM-irow); dRcol=TMath::Abs(icolM-icol); dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol); - if(dist < minDist){ + if (dist < minDist) + { //printf("MIN DISTANCE TO BAD %2.2f\n",dist); minDist = dist; } @@ -1297,32 +1423,36 @@ void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry } //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module - if (shared) { + if (shared) + { TH2D* hMap2 = 0; Int_t iSupMod2 = -1; //The only possible combinations are (0,1), (2,3) ... (8,9) - if(iSupMod%2) iSupMod2 = iSupMod-1; - else iSupMod2 = iSupMod+1; + if (iSupMod%2) iSupMod2 = iSupMod-1; + else iSupMod2 = iSupMod+1; hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2); //Loop on tower status map of second super module - for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){ - for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){ + for (Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++) + { + for (Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++) + { //Check if tower is bad. - if(hMap2->GetBinContent(icol,irow)==0) continue; + if (hMap2->GetBinContent(icol,irow)==0) + continue; //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels(shared) - \n \t Bad channel in SM %d, col %d, row %d \n \t Cluster max in SM %d, col %d, row %d\n", // iSupMod2,icol, irow,iSupMod,icolM,irowM); dRrow=TMath::Abs(irow-irowM); - if(iSupMod%2) { + if (iSupMod%2) { dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM)); } else { dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM); } dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol); - if(dist < minDist) minDist = dist; + if (dist < minDist) minDist = dist; } } }// shared cluster in 2 SuperModules @@ -1331,160 +1461,223 @@ void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry cluster->SetDistanceToBadChannel(minDist); } -//____________________________________________________________________________ +//__________________________________________________________________ void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster) { //re-evaluate identification parameters with bayesian - if(!cluster){ + if (!cluster) { AliInfo("Cluster pointer null!"); return; } - if ( cluster->GetM02() != 0) + if (cluster->GetM02() != 0) fPIDUtils->ComputePID(cluster->E(),cluster->GetM02()); - Float_t pidlist[AliPID::kSPECIESN+1]; - for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i); + Float_t pidlist[AliPID::kSPECIESCN+1]; + for (Int_t i = 0; i < AliPID::kSPECIESCN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i); cluster->SetPID(pidlist); } -//____________________________________________________________________________ -void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster) +//___________________________________________________________________________________________________________________ +void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom, + AliVCaloCells* cells, + AliVCluster * cluster, + Float_t & l0, Float_t & l1, + Float_t & disp, Float_t & dEta, Float_t & dPhi, + Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi) { // Calculates new center of gravity in the local EMCAL-module coordinates // and tranfers into global ALICE coordinates // Calculates Dispersion and main axis - if(!cluster){ + if (!cluster) { AliInfo("Cluster pointer null!"); return; } - Int_t nstat = 0; - Float_t wtot = 0. ; Double_t eCell = 0.; Float_t fraction = 1.; Float_t recalFactor = 1.; - Int_t iSupMod = -1; - Int_t iTower = -1; - Int_t iIphi = -1; - Int_t iIeta = -1; - Int_t iphi = -1; - Int_t ieta = -1; - Double_t etai = -1.; - Double_t phii = -1.; - - Double_t w = 0.; - Double_t d = 0.; - Double_t dxx = 0.; - Double_t dzz = 0.; - Double_t dxz = 0.; - Double_t xmean = 0.; - Double_t zmean = 0.; + Int_t iSupMod = -1; + Int_t iTower = -1; + Int_t iIphi = -1; + Int_t iIeta = -1; + Int_t iphi = -1; + Int_t ieta = -1; + Double_t etai = -1.; + Double_t phii = -1.; + + Int_t nstat = 0 ; + Float_t wtot = 0.; + Double_t w = 0.; + Double_t etaMean = 0.; + Double_t phiMean = 0.; + + //Loop on cells, calculate the cluster energy, in case a cut on cell energy is added + // and to check if the cluster is between 2 SM in eta + Int_t iSM0 = -1; + Bool_t shared = kFALSE; + Float_t energy = 0; + + for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) + { + //Get from the absid the supermodule, tower and eta/phi numbers + geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); + geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); - //Loop on cells - for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) { + //Check if there are cells of different SM + if (iDigit == 0 ) iSM0 = iSupMod; + else if (iSupMod!= iSM0) shared = kTRUE; + + //Get the cell energy, if recalibration is on, apply factors + fraction = cluster->GetCellAmplitudeFraction(iDigit); + if (fraction < 1e-4) fraction = 1.; // in case unfolding is off + + if (IsRecalibrationOn()) { + recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi); + } + + eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor; + + energy += eCell; + }//cell loop + + //Loop on cells + for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) + { //Get from the absid the supermodule, tower and eta/phi numbers geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); //Get the cell energy, if recalibration is on, apply factors fraction = cluster->GetCellAmplitudeFraction(iDigit); - if(fraction < 1e-4) fraction = 1.; // in case unfolding is off + if (fraction < 1e-4) fraction = 1.; // in case unfolding is off - if (!fCellsRecalibrated){ - if(IsRecalibrationOn()) { + if (!fCellsRecalibrated) { + if (IsRecalibrationOn()) { recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi); } } eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor; - if(cluster->E() > 0 && eCell > 0){ - + // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2 + // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0 + if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols; + + if (cluster->E() > 0 && eCell > 0) { w = GetCellWeight(eCell,cluster->E()); etai=(Double_t)ieta; - phii=(Double_t)iphi; - if(w > 0.0) { + phii=(Double_t)iphi; + + if (w > 0.0) { wtot += w ; - nstat++; + nstat++; //Shower shape - dxx += w * etai * etai ; - xmean+= w * etai ; - dzz += w * phii * phii ; - zmean+= w * phii ; - dxz += w * etai * phii ; + sEta += w * etai * etai ; + etaMean += w * etai ; + sPhi += w * phii * phii ; + phiMean += w * phii ; + sEtaPhi += w * etai * phii ; } } else AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E())); }//cell loop - //Normalize to the weight + //Normalize to the weight if (wtot > 0) { - xmean /= wtot ; - zmean /= wtot ; - } - else + etaMean /= wtot ; + phiMean /= wtot ; + } else AliError(Form("Wrong weight %f\n", wtot)); - //Calculate dispersion - for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) { - + //Calculate dispersion + for (Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) + { //Get from the absid the supermodule, tower and eta/phi numbers geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta); geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta); //Get the cell energy, if recalibration is on, apply factors fraction = cluster->GetCellAmplitudeFraction(iDigit); - if(fraction < 1e-4) fraction = 1.; // in case unfolding is off + if (fraction < 1e-4) fraction = 1.; // in case unfolding is off if (IsRecalibrationOn()) { recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi); } eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor; - if(cluster->E() > 0 && eCell > 0){ - + // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2 + // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0 + if (shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols; + + if (cluster->E() > 0 && eCell > 0) { w = GetCellWeight(eCell,cluster->E()); etai=(Double_t)ieta; - phii=(Double_t)iphi; - if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean)); - } - else + phii=(Double_t)iphi; + if (w > 0.0) { + disp += w *((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean)); + dEta += w * (etai-etaMean)*(etai-etaMean) ; + dPhi += w * (phii-phiMean)*(phii-phiMean) ; + } + } else AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E())); }// cell loop //Normalize to the weigth and set shower shape parameters if (wtot > 0 && nstat > 1) { - d /= wtot ; - dxx /= wtot ; - dzz /= wtot ; - dxz /= wtot ; - dxx -= xmean * xmean ; - dzz -= zmean * zmean ; - dxz -= xmean * zmean ; - cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )); - cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )); - } - else{ - d=0. ; - cluster->SetM20(0.) ; - cluster->SetM02(0.) ; - } - - if (d>=0) - cluster->SetDispersion(TMath::Sqrt(d)) ; - else - cluster->SetDispersion(0) ; + disp /= wtot ; + dEta /= wtot ; + dPhi /= wtot ; + sEta /= wtot ; + sPhi /= wtot ; + sEtaPhi /= wtot ; + + sEta -= etaMean * etaMean ; + sPhi -= phiMean * phiMean ; + sEtaPhi -= etaMean * phiMean ; + + l0 = (0.5 * (sEta + sPhi) + TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi )); + l1 = (0.5 * (sEta + sPhi) - TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi )); + } else { + l0 = 0. ; + l1 = 0. ; + dEta = 0. ; dPhi = 0. ; disp = 0. ; + sEta = 0. ; sPhi = 0. ; sEtaPhi = 0. ; + } } +//____________________________________________________________________________________________ +void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom, + AliVCaloCells* cells, + AliVCluster * cluster) +{ + // Calculates new center of gravity in the local EMCAL-module coordinates + // and tranfers into global ALICE coordinates + // Calculates Dispersion and main axis and puts them into the cluster + + Float_t l0 = 0., l1 = 0.; + Float_t disp = 0., dEta = 0., dPhi = 0.; + Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.; + + AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(geom,cells,cluster,l0,l1,disp, + dEta, dPhi, sEta, sPhi, sEtaPhi); + + cluster->SetM02(l0); + cluster->SetM20(l1); + if (disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ; + +} + //____________________________________________________________________________ -void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr, AliEMCALGeometry *geom) +void AliEMCALRecoUtils::FindMatches(AliVEvent *event, + TObjArray * clusterArr, + const AliEMCALGeometry *geom) { //This function should be called before the cluster loop //Before call this function, please recalculate the cluster positions @@ -1505,139 +1698,142 @@ void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr, Al AliAODEvent* aodevent = dynamic_cast (event); // init the magnetic field if not already on - if(!TGeoGlobalMagField::Instance()->GetField()){ - AliInfo("Init the magnetic field\n"); - if (esdevent) - { - esdevent->InitMagneticField(); - } - else if(aodevent) - { - Double_t curSol = 30000*aodevent->GetMagneticField()/5.00668; - Double_t curDip = 6000 *aodevent->GetMuonMagFieldScale(); - AliMagF *field = AliMagF::CreateFieldMap(curSol,curDip); - TGeoGlobalMagField::Instance()->SetField(field); - } - else + if (!TGeoGlobalMagField::Instance()->GetField()) { + if (!event->InitMagneticField()) { AliInfo("Mag Field not initialized, null esd/aod evetn pointers"); } - } // Init mag field + if (esdevent) { + UInt_t mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ(); + UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco(); + Bool_t desc1 = (mask1 >> 3) & 0x1; + Bool_t desc2 = (mask2 >> 3) & 0x1; + if (desc1==0 || desc2==0) { +// AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)", +// mask1, esdevent->GetESDRun()->GetDetectorsInReco(), +// mask2, esdevent->GetESDRun()->GetDetectorsInDAQ())); + fITSTrackSA=kTRUE; + } + } + TObjArray *clusterArray = 0x0; - if(!clusterArr) + if (!clusterArr) { + clusterArray = new TObjArray(event->GetNumberOfCaloClusters()); + for (Int_t icl=0; iclGetNumberOfCaloClusters(); icl++) { - clusterArray = new TObjArray(event->GetNumberOfCaloClusters()); - for(Int_t icl=0; iclGetNumberOfCaloClusters(); icl++) - { - AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl); - if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue; - clusterArray->AddAt(cluster,icl); - } + AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl); + if (geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) + continue; + clusterArray->AddAt(cluster,icl); } + } Int_t matched=0; Double_t cv[21]; for (Int_t i=0; i<21;i++) cv[i]=0; - for(Int_t itr=0; itrGetNumberOfTracks(); itr++) + for (Int_t itr=0; itrGetNumberOfTracks(); itr++) { AliExternalTrackParam *trackParam = 0; //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner AliESDtrack *esdTrack = 0; AliAODTrack *aodTrack = 0; - if(esdevent) - { + if (esdevent) { esdTrack = esdevent->GetTrack(itr); - if(!esdTrack) continue; - if(!IsAccepted(esdTrack)) continue; - if(esdTrack->Pt()Pt()Phi()*TMath::RadToDeg(); - if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue; - trackParam = const_cast(esdTrack->GetInnerParam()); + if (TMath::Abs(esdTrack->Eta())>0.9 || phi <= 10 || phi >= 250 ) continue; + if (!fITSTrackSA) + trackParam = const_cast(esdTrack->GetInnerParam()); // if TPC Available + else + trackParam = new AliExternalTrackParam(*esdTrack); // If ITS Track Standing alone } //If the input event is AOD, the starting point for extrapolation is at vertex //AOD tracks are selected according to its filterbit. - else if(aodevent) - { + else if (aodevent) { aodTrack = aodevent->GetTrack(itr); - if(!aodTrack) continue; - if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010() - if(aodTrack->Pt() \n", aodTrack->IsTPCOnly(), aodTrack->TestFilterMask(128)); + if (!aodTrack->IsTPCOnly()) continue ; + } else if (fAODHybridTracks) { // Match with hybrid tracks + //printf("Match with Hybrid tracks, accept? %d \n", aodTrack->IsHybridGlobalConstrainedGlobal()); + if (!aodTrack->IsHybridGlobalConstrainedGlobal()) continue ; + } else { // Match with tracks on a mask + //printf("Match with tracks having filter bit mask %d, accept? %d \n",fAODFilterMask,aodTrack->TestFilterMask(fAODFilterMask)); + if (!aodTrack->TestFilterMask(fAODFilterMask) ) continue; //Select AOD tracks + } + + if (aodTrack->Pt()Phi()*TMath::RadToDeg(); - if(TMath::Abs(aodTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue; + if (TMath::Abs(aodTrack->Eta())>0.9 || phi <= 10 || phi >= 250 ) + continue; Double_t pos[3],mom[3]; aodTrack->GetXYZ(pos); aodTrack->GetPxPyPz(mom); AliDebug(5,Form("aod track: i=%d | pos=(%5.4f,%5.4f,%5.4f) | mom=(%5.4f,%5.4f,%5.4f) | charge=%d\n",itr,pos[0],pos[1],pos[2],mom[0],mom[1],mom[2],aodTrack->Charge())); + trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge()); } //Return if the input data is not "AOD" or "ESD" - else - { + else { printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n"); - if(clusterArray) - { - clusterArray->Clear(); - delete clusterArray; - } + if (clusterArray) { + clusterArray->Clear(); + delete clusterArray; + } return; } - if(!trackParam) continue; + if (!trackParam) continue; //Extrapolate the track to EMCal surface AliExternalTrackParam emcalParam(*trackParam); - Float_t eta, phi; - if(!ExtrapolateTrackToEMCalSurface(&emcalParam, 430., fMass, fStepSurface, eta, phi)) - { - if(aodevent && trackParam) delete trackParam; - continue; - } - - if(esdevent) - { - esdTrack->SetOuterParam(&emcalParam,AliExternalTrackParam::kMultSec); - } - - if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) - { - if(aodevent && trackParam) delete trackParam; - continue; - } + Float_t eta, phi, pt; + if (!ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt)) { + if (aodevent && trackParam) delete trackParam; + if (fITSTrackSA && trackParam) delete trackParam; + continue; + } + if (TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) { + if (aodevent && trackParam) delete trackParam; + if (fITSTrackSA && trackParam) delete trackParam; + continue; + } //Find matched clusters Int_t index = -1; Float_t dEta = -999, dPhi = -999; - if(!clusterArr) - { - index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi); - } - else - { - index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi); - } + if (!clusterArr) { + index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi); + } else { + index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi); + } - if(index>-1) - { + if (index>-1) { fMatchedTrackIndex ->AddAt(itr,matched); fMatchedClusterIndex ->AddAt(index,matched); fResidualEta ->AddAt(dEta,matched); fResidualPhi ->AddAt(dPhi,matched); matched++; } - if(aodevent && trackParam) delete trackParam; + if (aodevent && trackParam) delete trackParam; + if (fITSTrackSA && trackParam) delete trackParam; }//track loop - if(clusterArray) - { - clusterArray->Clear(); - delete clusterArray; - } + if (clusterArray) { + clusterArray->Clear(); + delete clusterArray; + } AliDebug(2,Form("Number of matched pairs = %d !\n",matched)); @@ -1648,87 +1844,101 @@ void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr, Al } //________________________________________________________________________________ -Int_t AliEMCALRecoUtils::FindMatchedClusterInEvent(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom, Float_t &dEta, Float_t &dPhi) +Int_t AliEMCALRecoUtils::FindMatchedClusterInEvent(const AliESDtrack *track, + const AliVEvent *event, + const AliEMCALGeometry *geom, + Float_t &dEta, Float_t &dPhi) { // // This function returns the index of matched cluster to input track // Returns -1 if no match is found Int_t index = -1; Double_t phiV = track->Phi()*TMath::RadToDeg(); - if(TMath::Abs(track->Eta())>0.8 || phiV <= 20 || phiV >= 240 ) return index; - AliExternalTrackParam *trackParam = const_cast(track->GetInnerParam()); - if(!trackParam) return index; + if (TMath::Abs(track->Eta())>0.9 || phiV <= 10 || phiV >= 250 ) return index; + AliExternalTrackParam *trackParam = 0; + if (!fITSTrackSA) + trackParam = const_cast(track->GetInnerParam()); // If TPC + else + trackParam = new AliExternalTrackParam(*track); + + if (!trackParam) return index; AliExternalTrackParam emcalParam(*trackParam); - Float_t eta, phi; - if(!ExtrapolateTrackToEMCalSurface(&emcalParam, 430., fMass, fStepSurface, eta, phi)) return index; - if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) return index; + Float_t eta, phi, pt; + if (!ExtrapolateTrackToEMCalSurface(&emcalParam, fEMCalSurfaceDistance, fMass, fStepSurface, eta, phi, pt)) { + if (fITSTrackSA) delete trackParam; + return index; + } + if (TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) { + if (fITSTrackSA) delete trackParam; + return index; + } + TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters()); - for(Int_t icl=0; iclGetNumberOfCaloClusters(); icl++) + for (Int_t icl=0; iclGetNumberOfCaloClusters(); icl++) { AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl); - if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue; + if (geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue; clusterArr->AddAt(cluster,icl); } index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi); clusterArr->Clear(); delete clusterArr; - + if (fITSTrackSA) delete trackParam; + return index; } -//________________________________________________________________________________ -Int_t AliEMCALRecoUtils::FindMatchedClusterInClusterArr(AliExternalTrackParam *emcalParam, AliExternalTrackParam *trkParam, TObjArray * clusterArr, Float_t &dEta, Float_t &dPhi) +//_______________________________________________________________________________________________ +Int_t AliEMCALRecoUtils::FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam, + AliExternalTrackParam *trkParam, + const TObjArray * clusterArr, + Float_t &dEta, Float_t &dPhi) { + // Find matched cluster in array + dEta=-999, dPhi=-999; Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi; Int_t index = -1; Float_t tmpEta=-999, tmpPhi=-999; Double_t exPos[3] = {0.,0.,0.}; - if(!emcalParam->GetXYZ(exPos)) return index; + if (!emcalParam->GetXYZ(exPos)) return index; Float_t clsPos[3] = {0.,0.,0.}; - for(Int_t icl=0; iclGetEntriesFast(); icl++) - { - AliVCluster *cluster = dynamic_cast (clusterArr->At(icl)) ; - if(!cluster || !cluster->IsEMCAL()) continue; - cluster->GetPosition(clsPos); - Double_t dR = TMath::Sqrt(TMath::Power(exPos[0]-clsPos[0],2)+TMath::Power(exPos[1]-clsPos[1],2)+TMath::Power(exPos[2]-clsPos[2],2)); - if(dR > fClusterWindow) continue; - - AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation - if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue; - if(fCutEtaPhiSum) - { - Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi); - if(tmpRGetEntriesFast(); icl++) + { + AliVCluster *cluster = dynamic_cast (clusterArr->At(icl)) ; + if (!cluster || !cluster->IsEMCAL()) continue; + cluster->GetPosition(clsPos); + Double_t dR = TMath::Sqrt(TMath::Power(exPos[0]-clsPos[0],2)+TMath::Power(exPos[1]-clsPos[1],2)+TMath::Power(exPos[2]-clsPos[2],2)); + if (dR > fClusterWindow) continue; + + AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation + if (!ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue; + if (fCutEtaPhiSum) { + Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi); + if (tmpRSetTrackPhiEtaPtOnEMCal(-999, -999, -999); + + if (track->Pt()<0.350) //todo: discuss with Marta, everywhere else we use pT=0 + return kFALSE; + + Double_t phi = track->Phi()*TMath::RadToDeg(); + if (TMath::Abs(track->Eta())>0.9 || phi <= 10 || phi >= 250) + return kFALSE; + + AliESDtrack *esdt = dynamic_cast(track); + AliAODTrack *aodt = 0; + if (!esdt) { + aodt = dynamic_cast(track); + if (!aodt) + return kFALSE; + } + + if (mass<0) { + Bool_t onlyTPC = kFALSE; + if (mass==-99) + onlyTPC=kTRUE; + if (esdt) + mass = esdt->GetMass(onlyTPC); + else + mass = aodt->M(); + } + + AliExternalTrackParam *trackParam = 0; + if (esdt) { + const AliExternalTrackParam *in = esdt->GetInnerParam(); + if (!in) + return kFALSE; + trackParam = new AliExternalTrackParam(*in); + } else { + Double_t xyz[3] = {0}, pxpypz[3] = {0}, cv[21] = {0}; + aodt->PxPyPz(pxpypz); + aodt->XvYvZv(xyz); + aodt->GetCovarianceXYZPxPyPz(cv); + trackParam = new AliExternalTrackParam(xyz,pxpypz,cv,aodt->Charge()); + } + if (!trackParam) + return kFALSE; + + Float_t etaout=-999, phiout=-999, ptout=-999; + Bool_t ret = ExtrapolateTrackToEMCalSurface(trackParam, + emcalR, + mass, + step, + etaout, + phiout, + ptout); + delete trackParam; + if (!ret) + return kFALSE; + if (TMath::Abs(etaout)>0.75 || (phiout<70*TMath::DegToRad()) || (phiout>190*TMath::DegToRad())) + return kFALSE; + track->SetTrackPhiEtaPtOnEMCal(phiout, etaout, ptout); + return kTRUE; +} + + //------------------------------------------------------------------------------------ Bool_t AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(AliExternalTrackParam *trkParam, Double_t emcalR, Double_t mass, Double_t step, Float_t &eta, - Float_t &phi) + Float_t &phi, + Float_t &pt) { //Extrapolate track to EMCAL surface - eta = -999, phi = -999; - if(!trkParam) return kFALSE; - if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, emcalR, mass, step, kTRUE, 0.8, -1)) return kFALSE; + eta = -999, phi = -999, pt = -999; + if (!trkParam) return kFALSE; + if (!AliTrackerBase::PropagateTrackToBxByBz(trkParam, emcalR, mass, step, kTRUE, 0.8, -1)) return kFALSE; Double_t trkPos[3] = {0.,0.,0.}; - if(!trkParam->GetXYZ(trkPos)) return kFALSE; + if (!trkParam->GetXYZ(trkPos)) return kFALSE; TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]); eta = trkPosVec.Eta(); phi = trkPosVec.Phi(); - if(phi<0) + pt = trkParam->Pt(); + if (phi<0) phi += 2*TMath::Pi(); return kTRUE; @@ -1773,13 +2052,13 @@ Bool_t AliEMCALRecoUtils::ExtrapolateTrackToPosition(AliExternalTrackParam *trkP // tmpEta = -999; tmpPhi = -999; - if(!trkParam) return kFALSE; + if (!trkParam) return kFALSE; Double_t trkPos[3] = {0.,0.,0.}; TVector3 vec(clsPos[0],clsPos[1],clsPos[2]); Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad(); vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system - if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), mass, step,kTRUE, 0.8, -1)) return kFALSE; - if(!trkParam->GetXYZ(trkPos)) return kFALSE; //Get the extrapolated global position + if (!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), mass, step,kTRUE, 0.8, -1)) return kFALSE; + if (!trkParam->GetXYZ(trkPos)) return kFALSE; //Get the extrapolated global position TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]); TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]); @@ -1793,7 +2072,7 @@ Bool_t AliEMCALRecoUtils::ExtrapolateTrackToPosition(AliExternalTrackParam *trkP //---------------------------------------------------------------------------------- Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, - AliVCluster *cluster, + const AliVCluster *cluster, Double_t mass, Double_t step, Float_t &tmpEta, @@ -1804,7 +2083,8 @@ Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkPa // tmpEta = -999; tmpPhi = -999; - if(!cluster || !trkParam) return kFALSE; + if (!cluster || !trkParam) + return kFALSE; Float_t clsPos[3] = {0.,0.,0.}; cluster->GetPosition(clsPos); @@ -1814,7 +2094,7 @@ Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkPa //--------------------------------------------------------------------------------- Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, - AliVCluster *cluster, + const AliVCluster *cluster, Float_t &tmpEta, Float_t &tmpPhi) { @@ -1825,15 +2105,15 @@ Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkPa return ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, tmpEta, tmpPhi); } -//_______________________________________________________________________________________ -void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi) +//_______________________________________________________________________ +void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, + Float_t &dEta, Float_t &dPhi) { //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex) //Get the residuals dEta and dPhi for this cluster to the closest track //Works with ESDs and AODs - if( FindMatchedPosForCluster(clsIndex) >= 999 ) - { + if (FindMatchedPosForCluster(clsIndex) >= 999) { AliDebug(2,"No matched tracks found!\n"); dEta=999.; dPhi=999.; @@ -1850,8 +2130,7 @@ void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta //Get the residuals dEta and dPhi for this track to the closest cluster //Works with ESDs and AODs - if( FindMatchedPosForTrack(trkIndex) >= 999 ) - { + if (FindMatchedPosForTrack(trkIndex) >= 999) { AliDebug(2,"No matched cluster found!\n"); dEta=999.; dPhi=999.; @@ -1868,7 +2147,7 @@ Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex) //Get the index of matched track to this cluster //Works with ESDs and AODs - if(IsClusterMatched(clsIndex)) + if (IsClusterMatched(clsIndex)) return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex)); else return -1; @@ -1881,35 +2160,35 @@ Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex) //Get the index of matched cluster to this track //Works with ESDs and AODs - if(IsTrackMatched(trkIndex)) + if (IsTrackMatched(trkIndex)) return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex)); else return -1; } -//__________________________________________________ +//______________________________________________________________ Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const { //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex) //Returns if the cluster has a match - if(FindMatchedPosForCluster(clsIndex) < 999) + if (FindMatchedPosForCluster(clsIndex) < 999) return kTRUE; else return kFALSE; } -//__________________________________________________ +//____________________________________________________________ Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const { //Given a track index as in AliESDEvent::GetTrack(trkIndex) //Returns if the track has a match - if(FindMatchedPosForTrack(trkIndex) < 999) + if (FindMatchedPosForTrack(trkIndex) < 999) return kTRUE; else return kFALSE; } -//__________________________________________________________ +//______________________________________________________________________ UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const { //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex) @@ -1917,20 +2196,22 @@ UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const Float_t tmpR = fCutR; UInt_t pos = 999; - for(Int_t i=0; iGetSize(); i++) { - if(fMatchedClusterIndex->At(i)==clsIndex) { + for (Int_t i=0; iGetSize(); i++) + { + if (fMatchedClusterIndex->At(i)==clsIndex) { Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i)); - if(rAt(i),fResidualEta->At(i),fResidualPhi->At(i))); + AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n", + fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i))); } } } return pos; } -//__________________________________________________________ +//____________________________________________________________________ UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const { //Given a track index as in AliESDEvent::GetTrack(trkIndex) @@ -1938,35 +2219,35 @@ UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const Float_t tmpR = fCutR; UInt_t pos = 999; - for(Int_t i=0; iGetSize(); i++) { - if(fMatchedTrackIndex->At(i)==trkIndex) { + for (Int_t i=0; iGetSize(); i++) + { + if (fMatchedTrackIndex->At(i)==trkIndex) { Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i)); - if(rAt(i),fResidualEta->At(i),fResidualPhi->At(i))); + AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n", + fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i))); } } } return pos; } -//___________________________________________________________________________________ -Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, - AliVCaloCells* cells,const Int_t bc) +//__________________________________________________________________________ +Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, + const AliEMCALGeometry *geom, + AliVCaloCells* cells, Int_t bc) { // check if the cluster survives some quality cut // // Bool_t isGood=kTRUE; - - if(!cluster || !cluster->IsEMCAL()) return kFALSE; - - if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE; - - if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE; - - if(IsExoticCluster(cluster, cells,bc)) return kFALSE; + + if (!cluster || !cluster->IsEMCAL()) return kFALSE; + if (ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE; + if (!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE; + if (IsExoticCluster(cluster, cells,bc)) return kFALSE; return isGood; } @@ -1994,13 +2275,11 @@ Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack) //DCA cuts - if(fTrackCutsType==kGlobalCut) - { - Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010() - //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep)); - SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane - } - + if (fTrackCutsType==kGlobalCut) { + Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010() + //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep)); + SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane + } Float_t b[2]; Float_t bCov[3]; @@ -2046,16 +2325,30 @@ Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack) if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ) cuts[9] = kTRUE; - if(fTrackCutsType==kGlobalCut) - { - //Require at least one SPD point + anything else in ITS - if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE) - cuts[10] = kTRUE; - } + if (fTrackCutsType==kGlobalCut) { + //Require at least one SPD point + anything else in ITS + if ( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE) + cuts[10] = kTRUE; + } + // ITS + if (fCutRequireITSStandAlone || fCutRequireITSpureSA) { + if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin)) { + // TPC tracks + cuts[11] = kTRUE; + } else { + // ITS standalone tracks + if (fCutRequireITSStandAlone && !fCutRequireITSpureSA) { + if (status & AliESDtrack::kITSpureSA) cuts[11] = kTRUE; + } else if (fCutRequireITSpureSA) { + if (!(status & AliESDtrack::kITSpureSA)) cuts[11] = kTRUE; + } + } + } + Bool_t cut=kFALSE; for (Int_t i=0; iGetNumberOfTracks(); - for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack) { - AliESDtrack* track = event->GetTrack(iTrack); - if (!track) { + for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack) + { + AliVTrack* track = dynamic_cast(event->GetTrack(iTrack)); + if (!track) + { AliWarning(Form("Could not receive track %d", iTrack)); continue; } - Int_t matchClusIndex = GetMatchedClusterIndex(iTrack); + + Int_t matchClusIndex = GetMatchedClusterIndex(iTrack); track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual - if(matchClusIndex != -1) - track->SetStatus(AliESDtrack::kEMCALmatch); - else - track->ResetStatus(AliESDtrack::kEMCALmatch); + /*the following can be done better if AliVTrack::SetStatus will be there. Patch pending with Andreas/Peter*/ + AliESDtrack* esdtrack = dynamic_cast(track); + if (esdtrack) { + if (matchClusIndex != -1) + esdtrack->SetStatus(AliESDtrack::kEMCALmatch); + else + esdtrack->ResetStatus(AliESDtrack::kEMCALmatch); + } else { + AliAODTrack* aodtrack = dynamic_cast(track); + if (matchClusIndex != -1) + aodtrack->SetStatus(AliESDtrack::kEMCALmatch); + else + aodtrack->ResetStatus(AliESDtrack::kEMCALmatch); + } } - AliDebug(2,"Track matched to closest cluster"); + AliDebug(2,"Track matched to closest cluster"); } //_________________________________________________________________________ -void AliEMCALRecoUtils::SetTracksMatchedToCluster(const AliESDEvent *event) +void AliEMCALRecoUtils::SetTracksMatchedToCluster(const AliVEvent *event) { // Checks if EMC clusters are matched to ESD track. // Adds track indexes of all the tracks matched to a cluster withing residuals in ESDCalocluster. - for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus) { - AliESDCaloCluster *cluster = event->GetCaloCluster(iClus); + for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus) + { + AliVCluster *cluster = event->GetCaloCluster(iClus); if (!cluster->IsEMCAL()) continue; @@ -2160,16 +2477,18 @@ void AliEMCALRecoUtils::SetTracksMatchedToCluster(const AliESDEvent *event) // Get the closest track matched to the cluster Int_t nMatched = 0; Int_t matchTrackIndex = GetMatchedTrackIndex(iClus); - if (matchTrackIndex != -1) { + if (matchTrackIndex != -1) + { arrayTrackMatched[nMatched] = matchTrackIndex; nMatched++; } // Get all other tracks matched to the cluster - for(Int_t iTrk=0; iTrkGetTrack(iTrk); - if(iTrk == matchTrackIndex) continue; - if(track->GetEMCALcluster() == iClus){ + for (Int_t iTrk=0; iTrk(event->GetTrack(iTrk)); + if (iTrk == matchTrackIndex) continue; + if (track->GetEMCALcluster() == iClus) { arrayTrackMatched[nMatched] = iTrk; ++nMatched; } @@ -2178,7 +2497,14 @@ void AliEMCALRecoUtils::SetTracksMatchedToCluster(const AliESDEvent *event) //printf("Tender::SetTracksMatchedToCluster - cluster E %f, N matches %d, first match %d\n",cluster->E(),nMatched,arrayTrackMatched[0]); arrayTrackMatched.Set(nMatched); - cluster->AddTracksMatched(arrayTrackMatched); + AliESDCaloCluster *esdcluster = dynamic_cast(cluster); + if (esdcluster) + esdcluster->AddTracksMatched(arrayTrackMatched); + else if (nMatched>0) { + AliAODCaloCluster *aodcluster = dynamic_cast(cluster); + if (aodcluster) + aodcluster->AddTrackMatched(event->GetTrack(arrayTrackMatched.At(0))); + } Float_t eta= -999, phi = -999; if (matchTrackIndex != -1) @@ -2186,10 +2512,9 @@ void AliEMCALRecoUtils::SetTracksMatchedToCluster(const AliESDEvent *event) cluster->SetTrackDistance(phi, eta); } - AliDebug(2,"Cluster matched to tracks"); + AliDebug(2,"Cluster matched to tracks"); } - //___________________________________________________ void AliEMCALRecoUtils::Print(const Option_t *) const { @@ -2197,70 +2522,35 @@ void AliEMCALRecoUtils::Print(const Option_t *) const printf("AliEMCALRecoUtils Settings: \n"); printf("Misalignment shifts\n"); - for(Int_t i=0; i<5; i++) printf("\t sector %d, traslation (x,y,z)=(%f,%f,%f), rotation (x,y,z)=(%f,%f,%f)\n",i, + for (Int_t i=0; i<5; i++) printf("\t sector %d, traslation (x,y,z)=(%f,%f,%f), rotation (x,y,z)=(%f,%f,%f)\n",i, fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2], fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] ); printf("Non linearity function %d, parameters:\n", fNonLinearityFunction); - for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]); + for (Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]); printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration); printf("Matching criteria: "); - if(fCutEtaPhiSum) - { - printf("sqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR); - } - else if(fCutEtaPhiSeparate) - { - printf("dEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi); - } - else - { - printf("Error\n"); - printf("please specify your cut criteria\n"); - printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n"); - printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n"); - } + if (fCutEtaPhiSum) { + printf("sqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR); + } else if (fCutEtaPhiSeparate) { + printf("dEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi); + } else { + printf("Error\n"); + printf("please specify your cut criteria\n"); + printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n"); + printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n"); + } printf("Mass hypothesis = %2.3f [GeV/c^2], extrapolation step to surface = %2.2f[cm], step to cluster = %2.2f[cm]\n",fMass,fStepSurface, fStepCluster); printf("Cluster selection window: dR < %2.0f\n",fClusterWindow); printf("Track cuts: \n"); printf("Minimum track pT: %1.2f\n",fCutMinTrackPt); - printf("AOD track selection mask: %d\n",fAODFilterMask); + printf("AOD track selection: tpc only %d, or hybrid %d, or mask: %d\n",fAODTPCOnlyTracks,fAODHybridTracks, fAODFilterMask); printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit); printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters); printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS); printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS); printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ); } - -//_____________________________________________________________________ -void AliEMCALRecoUtils::SetRunDependentCorrections(Int_t runnumber) -{ - //Get EMCAL time dependent corrections from file and put them in the recalibration histograms - //Do it only once and only if it is requested - - if(!fUseRunCorrectionFactors) return; - if(fRunCorrectionFactorsSet) return; - - AliInfo(Form("AliEMCALRecoUtils::GetRunDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber)); - - AliEMCALCalibTimeDepCorrection *corr = new AliEMCALCalibTimeDepCorrection(); - corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber)); - - SwitchOnRecalibration(); - for(Int_t ism = 0; ism < 4; ism++){ - for(Int_t icol = 0; icol < 48; icol++){ - for(Int_t irow = 0; irow < 24; irow++){ - Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow); - Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0); - GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor); - //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow, - // orgRecalFactor, corr->GetCorrection(ism, icol,irow,0), - // (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor); - } - } - } - fRunCorrectionFactorsSet = kTRUE; -}