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