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