EMCAL/AliEMCALRecoUtils.h

   1 #ifndef ALIEMCALRECOUTILS_H
   2 #define ALIEMCALRECOUTILS_H
   3
   4 /* $Id: AliEMCALRecoUtils.h 33808 2009-07-15 09:48:08Z gconesab $ */
   5
   6 ///////////////////////////////////////////////////////////////////////////////
   7 //
   8 // Class AliEMCALRecoUtils
   9 // Some utilities to recalculate the cluster position or energy linearity
  10 //
  11 //
  12 // Author:  Gustavo Conesa (LPSC- Grenoble)
  13 ///////////////////////////////////////////////////////////////////////////////
  14
  15 //Root includes
  16 #include "TNamed.h"
  17 #include "TMath.h"
  18 #include "TObjArray.h"
  19 #include "TArrayI.h"
  20 #include "TArrayF.h"
  21 #include "TH2F.h"
  22
  23 //AliRoot includes
  24 class AliVCluster;
  25 class AliVCaloCells;
  26 class AliVEvent;
  27 #include "AliLog.h"
  28 class AliEMCALGeometry;
  29 class AliEMCALPIDUtils;
  30 class AliESDtrack;
  31
  32 class AliEMCALRecoUtils : public TNamed {
  33
  34 public:
  35
  36   AliEMCALRecoUtils();
  37   AliEMCALRecoUtils(const AliEMCALRecoUtils&);
  38   AliEMCALRecoUtils& operator=(const AliEMCALRecoUtils&);
  39   virtual ~AliEMCALRecoUtils() ;
  40
  41   enum NonlinearityFunctions{kPi0MC=0,kPi0GammaGamma=1,kPi0GammaConversion=2,kNoCorrection=3,kBeamTest=4};
  42   enum PositionAlgorithms{kUnchanged=-1,kPosTowerIndex=0, kPosTowerGlobal=1};
  43   enum ParticleType{kPhoton=0, kElectron=1,kHadron =2, kUnknown=-1};
  44
  45   //Position recalculation
  46   void     RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu);
  47   void     RecalculateClusterPositionFromTowerIndex (AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu);
  48   void     RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu);
  49
  50   Float_t  GetCellWeight(const Float_t eCell, const Float_t eCluster) const { return TMath::Max( 0., fW0 + TMath::Log( eCell / eCluster ));}
  51
  52   Float_t  GetDepth(const Float_t eCluster, const Int_t iParticle, const Int_t iSM) const ;
  53
  54   void     GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
  55                             Int_t & absId,  Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared);
  56
  57   Float_t  GetMisalTransShift(const Int_t i) const {
  58     if(i < 15 ){return fMisalTransShift[i]; }
  59     else { AliInfo(Form("Index %d larger than 15, do nothing\n",i)); return 0.;}
  60   }
  61   Float_t* GetMisalTransShiftArray() {return fMisalTransShift; }
  62
  63   void     SetMisalTransShift(const Int_t i, const Float_t shift) {
  64     if(i < 15 ){fMisalTransShift[i] = shift; }
  65     else { AliInfo(Form("Index %d larger than 15, do nothing\n",i));}
  66   }
  67   void     SetMisalTransShiftArray(Float_t * misal)
  68   { for(Int_t i = 0; i < 15; i++)fMisalTransShift[i] = misal[i]; }
  69
  70   Float_t  GetMisalRotShift(const Int_t i) const {
  71     if(i < 15 ){return fMisalRotShift[i]; }
  72     else { AliInfo(Form("Index %d larger than 15, do nothing\n",i)); return 0.;}
  73   }
  74   Float_t* GetMisalRotShiftArray() {return fMisalRotShift; }
  75
  76   void     SetMisalRotShift(const Int_t i, const Float_t shift) {
  77     if(i < 15 ){fMisalRotShift[i] = shift; }
  78     else { AliInfo(Form("Index %d larger than 15, do nothing\n",i));}
  79   }
  80   void     SetMisalRotShiftArray(Float_t * misal)
  81   { for(Int_t i = 0; i < 15; i++)fMisalRotShift[i] = misal[i]; }
  82
  83   Int_t    GetParticleType()        const  { return  fParticleType    ;}
  84   void     SetParticleType(Int_t particle) { fParticleType = particle ;}
  85
  86   Int_t    GetPositionAlgorithm()   const  { return fPosAlgo          ;}
  87   void     SetPositionAlgorithm(Int_t alg) { fPosAlgo = alg           ;}
  88
  89   Float_t  GetW0()                 const   { return fW0               ;}
  90   void     SetW0(Float_t w0)               { fW0  = w0                ;}
  91
  92   //Non Linearity
  93
  94   Float_t CorrectClusterEnergyLinearity(AliVCluster* clu);
  95
  96   Float_t  GetNonLinearityParam(const Int_t i) const {
  97     if(i < 6 ){return fNonLinearityParams[i]; }
  98     else { AliInfo(Form("Index %d larger than 6, do nothing\n",i)); return 0.;}
  99   }
 100   void     SetNonLinearityParam(const Int_t i, const Float_t param) {
 101     if(i < 6 ){fNonLinearityParams[i] = param; }
 102     else { AliInfo(Form("Index %d larger than 6, do nothing\n",i));}
 103   }
 104
 105   Int_t GetNonLinearityFunction() const    { return fNonLinearityFunction ;}
 106   void  SetNonLinearityFunction(Int_t fun) { fNonLinearityFunction = fun  ;}
 107
 108   void Print(const Option_t*) const;
 109
 110   //Recalibration
 111   void RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells * cells);
 112
 113   Bool_t IsRecalibrationOn()  const        { return fRecalibration ; }
 114   void SwitchOnRecalibration()             { fRecalibration = kTRUE ; if(!fEMCALRecalibrationFactors)InitEMCALRecalibrationFactors();}
 115   void SwitchOffRecalibration()            { fRecalibration = kFALSE ; }
 116   void InitEMCALRecalibrationFactors() ;
 117
 118   //Recalibrate channels with time dependent corrections
 119   void SwitchOnTimeDepCorrection()        { fUseTimeCorrectionFactors = kTRUE ; SwitchOnRecalibration();}
 120   void SwitchOffTimeDepCorrection()       { fUseTimeCorrectionFactors = kFALSE;}
 121   void SetTimeDependentCorrections(Int_t runnumber);
 122
 123   Float_t GetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow) const {
 124     if(fEMCALRecalibrationFactors) return (Float_t) ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->GetBinContent(iCol,iRow);
 125     else return 1;}
 126
 127   void SetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow, Double_t c = 1) {
 128     if(!fEMCALRecalibrationFactors) InitEMCALRecalibrationFactors();
 129     ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->SetBinContent(iCol,iRow,c);}
 130
 131   TH2F * GetEMCALChannelRecalibrationFactors(Int_t iSM)   const { return (TH2F*)fEMCALRecalibrationFactors->At(iSM) ;}
 132   void SetEMCALChannelRecalibrationFactors(TObjArray *map)      { fEMCALRecalibrationFactors = map                  ;}
 133   void SetEMCALChannelRecalibrationFactors(Int_t iSM , TH2F* h) { fEMCALRecalibrationFactors->AddAt(h,iSM)          ;}
 134
 135   //Modules fiducial region, remove clusters in borders
 136   Bool_t CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells) ;
 137   void   SetNumberOfCellsFromEMCALBorder(Int_t n) { fNCellsFromEMCALBorder = n    ;}
 138   Int_t  GetNumberOfCellsFromEMCALBorder() const  { return fNCellsFromEMCALBorder ;}
 139
 140   void   SwitchOnNoFiducialBorderInEMCALEta0()    { fNoEMCALBorderAtEta0 = kTRUE  ;}
 141   void   SwitchOffNoFiducialBorderInEMCALEta0()   { fNoEMCALBorderAtEta0 = kFALSE ;}
 142   Bool_t IsEMCALNoBorderAtEta0()                  { return fNoEMCALBorderAtEta0   ;}
 143
 144   // Bad channels
 145   Bool_t IsBadChannelsRemovalSwitchedOn()  const { return fRemoveBadChannels      ;}
 146   void SwitchOnBadChannelsRemoval ()             { fRemoveBadChannels = kTRUE ; if(!fEMCALBadChannelMap)InitEMCALBadChannelStatusMap();}
 147   void SwitchOffBadChannelsRemoval()             { fRemoveBadChannels = kFALSE    ;}
 148
 149   Bool_t IsDistanceToBadChannelRecalculated() const { return fRecalDistToBadChannels;}
 150   void SwitchOnDistToBadChannelRecalculation()   { fRecalDistToBadChannels = kTRUE  ; if(!fEMCALBadChannelMap)InitEMCALBadChannelStatusMap();}
 151   void SwitchOffDistToBadChannelRecalculation()  { fRecalDistToBadChannels = kFALSE ;}
 152
 153   void InitEMCALBadChannelStatusMap() ;
 154
 155   Int_t GetEMCALChannelStatus(Int_t iSM , Int_t iCol, Int_t iRow) const {
 156     if(fEMCALBadChannelMap) return (Int_t) ((TH2I*)fEMCALBadChannelMap->At(iSM))->GetBinContent(iCol,iRow);
 157     else return 0;}//Channel is ok by default
 158
 159   void SetEMCALChannelStatus(Int_t iSM , Int_t iCol, Int_t iRow, Double_t c = 1) {
 160     if(!fEMCALBadChannelMap)InitEMCALBadChannelStatusMap() ;
 161     ((TH2I*)fEMCALBadChannelMap->At(iSM))->SetBinContent(iCol,iRow,c);}
 162
 163   TH2I * GetEMCALChannelStatusMap(Int_t iSM) const {return (TH2I*)fEMCALBadChannelMap->At(iSM);}
 164   void   SetEMCALChannelStatusMap(TObjArray *map)  {fEMCALBadChannelMap = map;}
 165   void   SetEMCALChannelStatusMap(Int_t iSM , TH2I* h) {fEMCALBadChannelMap->AddAt(h,iSM);}
 166
 167   Bool_t ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells);
 168
 169   //Recalculate other cluster parameters
 170   void RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster);
 171   void RecalculateClusterPID(AliVCluster * cluster);
 172
 173   AliEMCALPIDUtils * GetPIDUtils() { return fPIDUtils;}
 174
 175   void RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster);
 176
 177   //Track matching
 178   void FindMatches(AliVEvent *event);
 179   void GetMatchedResiduals(Int_t index, Float_t &dR, Float_t &dZ);
 180   Bool_t IsMatched(Int_t index);
 181   UInt_t FindMatchedPos(Int_t index) const;
 182
 183   Float_t GetCutR()    const { return fCutR ;}
 184   Float_t GetCutZ()    const { return fCutZ ;}
 185   void SetCutR(Float_t cutR) { fCutR=cutR   ;}
 186   void SetCutZ(Float_t cutZ) { fCutZ=cutZ   ;}
 187
 188   //Track Cuts
 189   Bool_t IsAccepted(AliESDtrack *track);
 190   void InitTrackCuts();
 191
 192   // track quality cut setters
 193   void SetMinNClustersTPC(Int_t min=-1)          {fCutMinNClusterTPC       = min  ;}
 194   void SetMinNClustersITS(Int_t min=-1)          {fCutMinNClusterITS       = min  ;}
 195   void SetMaxChi2PerClusterTPC(Float_t max=1e10) {fCutMaxChi2PerClusterTPC = max  ;}
 196   void SetMaxChi2PerClusterITS(Float_t max=1e10) {fCutMaxChi2PerClusterITS = max  ;}
 197   void SetRequireTPCRefit(Bool_t b=kFALSE)       {fCutRequireTPCRefit      = b    ;}
 198   void SetRequireITSRefit(Bool_t b=kFALSE)       {fCutRequireITSRefit      = b    ;}
 199   void SetAcceptKinkDaughters(Bool_t b=kTRUE)    {fCutAcceptKinkDaughters  = b    ;}
 200   void SetMaxDCAToVertexXY(Float_t dist=1e10)    {fCutMaxDCAToVertexXY     = dist ;}
 201   void SetMaxDCAToVertexZ(Float_t dist=1e10)     {fCutMaxDCAToVertexZ      = dist ;}
 202   void SetDCAToVertex2D(Bool_t b=kFALSE)         {fCutDCAToVertex2D        = b    ;}
 203
 204   // getters
 205
 206   Int_t   GetMinNClusterTPC()        const   { return fCutMinNClusterTPC;}
 207   Int_t   GetMinNClustersITS()       const   { return fCutMinNClusterITS;}
 208   Float_t GetMaxChi2PerClusterTPC()  const   { return fCutMaxChi2PerClusterTPC;}
 209   Float_t GetMaxChi2PerClusterITS()  const   { return fCutMaxChi2PerClusterITS;}
 210   Bool_t  GetRequireTPCRefit()       const   { return fCutRequireTPCRefit;}
 211   Bool_t  GetRequireITSRefit()       const   { return fCutRequireITSRefit;}
 212   Bool_t  GetAcceptKinkDaughters()   const   { return fCutAcceptKinkDaughters;}
 213   Float_t GetMaxDCAToVertexXY()      const   { return fCutMaxDCAToVertexXY;}
 214   Float_t GetMaxDCAToVertexZ()       const   { return fCutMaxDCAToVertexZ;}
 215   Bool_t  GetDCAToVertex2D()         const   { return fCutDCAToVertex2D;}
 216
 217
 218 private:
 219
 220   Float_t    fMisalTransShift[15];       // Shift parameters
 221   Float_t    fMisalRotShift[15];         // Shift parameters
 222   Int_t      fNonLinearityFunction;      // Non linearity function choice
 223   Float_t    fNonLinearityParams[6];     // Parameters for the non linearity function
 224   Int_t      fParticleType;              // Particle type for depth calculation
 225   Int_t      fPosAlgo;                   // Position recalculation algorithm
 226   Float_t    fW0;                        // Weight0
 227
 228   Bool_t     fRecalibration;             // Switch on or off the recalibration
 229   TObjArray* fEMCALRecalibrationFactors; // Array of histograms with map of recalibration factors, EMCAL
 230   Bool_t     fRemoveBadChannels;         // Check the channel status provided and remove clusters with bad channels
 231   Bool_t     fRecalDistToBadChannels;    // Calculate distance from highest energy tower of cluster to closes bad channel
 232   TObjArray* fEMCALBadChannelMap;        // Array of histograms with map of bad channels, EMCAL
 233   Int_t      fNCellsFromEMCALBorder;     // Number of cells from EMCAL border the cell with maximum amplitude has to be.
 234   Bool_t     fNoEMCALBorderAtEta0;       // Do fiducial cut in EMCAL region eta = 0?
 235
 236   TArrayI  * fMatchedClusterIndex;       // Array that stores indexes of matched clusters
 237   TArrayF  * fResidualZ;                 // Array that stores the residual z
 238   TArrayF  * fResidualR;                 // Array that stores the residual r
 239   Float_t    fCutR;                      // dR cut on matching
 240   Float_t    fCutZ;                      // dZ cut on matching
 241
 242   enum { kNCuts = 11 };
 243   Int_t      fCutMinNClusterTPC;         // Min number of tpc clusters
 244   Int_t      fCutMinNClusterITS;         // Min number of its clusters
 245   Float_t    fCutMaxChi2PerClusterTPC;   // Max tpc fit chi2 per tpc cluster
 246   Float_t    fCutMaxChi2PerClusterITS;   // Max its fit chi2 per its cluster
 247   Bool_t     fCutRequireTPCRefit;        // Require TPC refit
 248   Bool_t     fCutRequireITSRefit;        // Require ITS refit
 249   Bool_t     fCutAcceptKinkDaughters;    // Accepting kink daughters?
 250   Float_t    fCutMaxDCAToVertexXY;       // Track-to-vertex cut in max absolute distance in xy-plane
 251   Float_t    fCutMaxDCAToVertexZ;        // Track-to-vertex cut in max absolute distance in z-plane
 252   Bool_t     fCutDCAToVertex2D;          // If true a 2D DCA cut is made. Tracks are accepted if sqrt((DCAXY / fCutMaxDCAToVertexXY)^2 + (DCAZ / fCutMaxDCAToVertexZ)^2) < 1 AND sqrt((DCAXY / fCutMinDCAToVertexXY)^2 + (DCAZ / fCutMinDCAToVertexZ)^2) > 1
 253
 254   AliEMCALPIDUtils * fPIDUtils;          // Recalculate PID parameters
 255
 256   //Time Correction
 257   Bool_t     fUseTimeCorrectionFactors;  // Use Time Dependent Correction
 258   Bool_t     fTimeCorrectionFactorsSet;  // Time Correction set at leat once
 259
 260   ClassDef(AliEMCALRecoUtils, 6)
 261
 262 };
 263
 264 #endif // ALIEMCALRECOUTILS_H
 265
 266