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
 117   void InitEMCALRecalibrationFactors() ;
 118
 119   Float_t GetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow) const {
 120     if(fEMCALRecalibrationFactors) return (Float_t) ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->GetBinContent(iCol,iRow);
 121     else return 1;}
 122
 123   void SetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow, Double_t c = 1) {
 124     if(!fEMCALRecalibrationFactors) InitEMCALRecalibrationFactors();
 125     ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->SetBinContent(iCol,iRow,c);}
 126
 127   TH2F * GetEMCALChannelRecalibrationFactors(Int_t iSM) const {return (TH2F*)fEMCALRecalibrationFactors->At(iSM);}
 128   void SetEMCALChannelRecalibrationFactors(TObjArray *map)      {fEMCALRecalibrationFactors = map;}
 129   void SetEMCALChannelRecalibrationFactors(Int_t iSM , TH2F* h) {fEMCALRecalibrationFactors->AddAt(h,iSM);}
 130
 131   //Modules fiducial region, remove clusters in borders
 132   Bool_t CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells) ;
 133   void   SetNumberOfCellsFromEMCALBorder(Int_t n) {fNCellsFromEMCALBorder = n; }
 134   Int_t  GetNumberOfCellsFromEMCALBorder() const  {return fNCellsFromEMCALBorder; }
 135
 136   void   SwitchOnNoFiducialBorderInEMCALEta0()  {fNoEMCALBorderAtEta0 = kTRUE; }
 137   void   SwitchOffNoFiducialBorderInEMCALEta0() {fNoEMCALBorderAtEta0 = kFALSE; }
 138   Bool_t IsEMCALNoBorderAtEta0()                {return fNoEMCALBorderAtEta0;}
 139
 140   // Bad channels
 141   Bool_t IsBadChannelsRemovalSwitchedOn()  const { return fRemoveBadChannels ; }
 142   void SwitchOnBadChannelsRemoval ()  {fRemoveBadChannels = kTRUE  ; if(!fEMCALBadChannelMap)InitEMCALBadChannelStatusMap();}
 143   void SwitchOffBadChannelsRemoval()  {fRemoveBadChannels = kFALSE ; }
 144
 145   Bool_t IsDistanceToBadChannelRecalculated()  const { return fRecalDistToBadChannels ; }
 146   void SwitchOnDistToBadChannelRecalculation()   {fRecalDistToBadChannels = kTRUE  ; if(!fEMCALBadChannelMap)InitEMCALBadChannelStatusMap();}
 147   void SwitchOffDistToBadChannelRecalculation()  {fRecalDistToBadChannels = kFALSE ; }
 148
 149   void InitEMCALBadChannelStatusMap() ;
 150
 151   Int_t GetEMCALChannelStatus(Int_t iSM , Int_t iCol, Int_t iRow) const {
 152     if(fEMCALBadChannelMap) return (Int_t) ((TH2I*)fEMCALBadChannelMap->At(iSM))->GetBinContent(iCol,iRow);
 153     else return 0;}//Channel is ok by default
 154
 155   void SetEMCALChannelStatus(Int_t iSM , Int_t iCol, Int_t iRow, Double_t c = 1) {
 156     if(!fEMCALBadChannelMap)InitEMCALBadChannelStatusMap() ;
 157     ((TH2I*)fEMCALBadChannelMap->At(iSM))->SetBinContent(iCol,iRow,c);}
 158
 159   TH2I * GetEMCALChannelStatusMap(Int_t iSM) const {return (TH2I*)fEMCALBadChannelMap->At(iSM);}
 160   void   SetEMCALChannelStatusMap(TObjArray *map)  {fEMCALBadChannelMap = map;}
 161   void   SetEMCALChannelStatusMap(Int_t iSM , TH2I* h) {fEMCALBadChannelMap->AddAt(h,iSM);}
 162
 163   Bool_t ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells);
 164
 165   //Recalculate other cluster parameters
 166   void RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster);
 167   void RecalculateClusterPID(AliVCluster * cluster);
 168
 169   AliEMCALPIDUtils * GetPIDUtils() { return fPIDUtils;}
 170
 171   void RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster);
 172
 173   //Track matching
 174   void FindMatches(AliVEvent *event);
 175   void GetMatchedResiduals(Int_t index, Float_t &dR, Float_t &dZ);
 176   Bool_t IsMatched(Int_t index);
 177   UInt_t FindMatchedPos(Int_t index) const;
 178
 179   Float_t GetCutR() const { return fCutR; }
 180   Float_t GetCutZ() const { return fCutZ; }
 181
 182   void SetCutR(Float_t cutR) { fCutR=cutR; }
 183   void SetCutZ(Float_t cutZ) { fCutZ=cutZ; }
 184
 185   //Track Cuts
 186   Bool_t IsAccepted(AliESDtrack *track);
 187   void InitTrackCuts();
 188
 189   // track quality cut setters
 190   void SetMinNClustersTPC(Int_t min=-1)          {fCutMinNClusterTPC=min;}
 191   void SetMinNClustersITS(Int_t min=-1)          {fCutMinNClusterITS=min;}
 192   void SetMaxChi2PerClusterTPC(Float_t max=1e10) {fCutMaxChi2PerClusterTPC=max;}
 193   void SetMaxChi2PerClusterITS(Float_t max=1e10) {fCutMaxChi2PerClusterITS=max;}
 194   void SetRequireTPCRefit(Bool_t b=kFALSE)       {fCutRequireTPCRefit=b;}
 195   void SetRequireITSRefit(Bool_t b=kFALSE)       {fCutRequireITSRefit=b;}
 196   void SetAcceptKinkDaughters(Bool_t b=kTRUE)    {fCutAcceptKinkDaughters=b;}
 197   void SetMaxDCAToVertexXY(Float_t dist=1e10)         {fCutMaxDCAToVertexXY = dist;}
 198   void SetMaxDCAToVertexZ(Float_t dist=1e10)          {fCutMaxDCAToVertexZ = dist;}
 199   void SetDCAToVertex2D(Bool_t b=kFALSE)              {fCutDCAToVertex2D = b;}
 200
 201   // getters
 202
 203   Int_t   GetMinNClusterTPC()        const   { return fCutMinNClusterTPC;}
 204   Int_t   GetMinNClustersITS()       const   { return fCutMinNClusterITS;}
 205   Float_t GetMaxChi2PerClusterTPC()  const   { return fCutMaxChi2PerClusterTPC;}
 206   Float_t GetMaxChi2PerClusterITS()  const   { return fCutMaxChi2PerClusterITS;}
 207   Bool_t  GetRequireTPCRefit()       const   { return fCutRequireTPCRefit;}
 208   Bool_t  GetRequireITSRefit()       const   { return fCutRequireITSRefit;}
 209   Bool_t  GetAcceptKinkDaughters()   const   { return fCutAcceptKinkDaughters;}
 210   Float_t GetMaxDCAToVertexXY()      const   { return fCutMaxDCAToVertexXY;}
 211   Float_t GetMaxDCAToVertexZ()       const   { return fCutMaxDCAToVertexZ;}
 212   Bool_t  GetDCAToVertex2D()         const   { return fCutDCAToVertex2D;}
 213
 214
 215 private:
 216
 217   Float_t fMisalTransShift[15];   // Shift parameters
 218   Float_t fMisalRotShift[15];     // Shift parameters
 219   Int_t   fNonLinearityFunction;  // Non linearity function choice
 220   Float_t fNonLinearityParams[6]; // Parameters for the non linearity function
 221   Int_t   fParticleType;          // Particle type for depth calculation
 222   Int_t   fPosAlgo;               // Position recalculation algorithm
 223   Float_t fW0;                    // Weight0
 224
 225   Bool_t     fRecalibration;             // Switch on or off the recalibration
 226   TObjArray* fEMCALRecalibrationFactors; // Array of histograms with map of recalibration factors, EMCAL
 227   Bool_t     fRemoveBadChannels;         // Check the channel status provided and remove clusters with bad channels
 228   Bool_t     fRecalDistToBadChannels;    // Calculate distance from highest energy tower of cluster to closes bad channel
 229   TObjArray* fEMCALBadChannelMap;        // Array of histograms with map of bad channels, EMCAL
 230   Int_t      fNCellsFromEMCALBorder;     // Number of cells from EMCAL border the cell with maximum amplitude has to be.
 231   Bool_t     fNoEMCALBorderAtEta0;       // Do fiducial cut in EMCAL region eta = 0?
 232
 233   TArrayI *fMatchedClusterIndex;  //Array that stores indexes of matched clusters
 234   TArrayF *fResidualZ;            //Array that stores the residual z
 235   TArrayF *fResidualR;            //Array that stores the residual r
 236   Float_t fCutR; //dR cut on matching
 237   Float_t fCutZ; //dZ cut on matching
 238
 239   enum { kNCuts = 11 };
 240   Int_t   fCutMinNClusterTPC;         // min number of tpc clusters
 241   Int_t   fCutMinNClusterITS;         // min number of its clusters
 242   Float_t fCutMaxChi2PerClusterTPC;   // max tpc fit chi2 per tpc cluster
 243   Float_t fCutMaxChi2PerClusterITS;   // max its fit chi2 per its cluster
 244   Bool_t  fCutRequireTPCRefit;        // require TPC refit
 245   Bool_t  fCutRequireITSRefit;        // require ITS refit
 246   Bool_t  fCutAcceptKinkDaughters;    // accepting kink daughters?
 247   Float_t fCutMaxDCAToVertexXY;       // track-to-vertex cut in max absolute distance in xy-plane
 248   Float_t fCutMaxDCAToVertexZ;        // track-to-vertex cut in max absolute distance in z-plane
 249   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
 250
 251   AliEMCALPIDUtils * fPIDUtils;               // Recalculate PID parameters
 252
 253   ClassDef(AliEMCALRecoUtils, 6)
 254
 255 };
 256
 257 #endif // ALIEMCALRECOUTILS_H
 258
 259