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 < 7 ){return fNonLinearityParams[i]; }
 110     else { AliInfo(Form("Index %d larger than 7, do nothing\n",i)); return 0.;}
 111   }
 112   void     SetNonLinearityParam(const Int_t i, const Float_t param) {
 113     if(i < 7 ){fNonLinearityParams[i] = param; }
 114     else { AliInfo(Form("Index %d larger than 7, do nothing\n",i));}
 115   }
 116   void  InitNonLinearityParam();
 117
 118   Int_t GetNonLinearityFunction() const    { return fNonLinearityFunction ; }
 119   void  SetNonLinearityFunction(Int_t fun) { fNonLinearityFunction = fun  ; InitNonLinearityParam() ; }
 120
 121   void  SetNonLinearityThreshold(Int_t threshold) {fNonLinearThreshold = threshold ;} //only for Alexie's non linearity correction
 122   Int_t GetNonLinearityThreshold() const   {return fNonLinearThreshold ;}
 123
 124
 125   //-----------------------------------------------------
 126   //Recalibration
 127   //-----------------------------------------------------
 128
 129   void RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells * cells);
 130
 131   Bool_t IsRecalibrationOn()           const { return fRecalibration ; }
 132   void   SwitchOnRecalibration()             { fRecalibration = kTRUE ; if(!fEMCALRecalibrationFactors)InitEMCALRecalibrationFactors();}
 133   void   SwitchOffRecalibration()            { fRecalibration = kFALSE ; }
 134   void   InitEMCALRecalibrationFactors() ;
 135
 136   //Recalibrate channels with time dependent corrections
 137   void SwitchOnTimeDepCorrection()          { fUseTimeCorrectionFactors = kTRUE ; SwitchOnRecalibration();}
 138   void SwitchOffTimeDepCorrection()         { fUseTimeCorrectionFactors = kFALSE;}
 139   void SetTimeDependentCorrections(Int_t runnumber);
 140
 141   Float_t GetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow) const {
 142     if(fEMCALRecalibrationFactors) return (Float_t) ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->GetBinContent(iCol,iRow);
 143     else return 1;}
 144
 145   void SetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow, Double_t c = 1) {
 146     if(!fEMCALRecalibrationFactors) InitEMCALRecalibrationFactors();
 147     ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->SetBinContent(iCol,iRow,c);}
 148
 149   TH2F * GetEMCALChannelRecalibrationFactors(Int_t iSM)   const { return (TH2F*)fEMCALRecalibrationFactors->At(iSM) ;}
 150   void SetEMCALChannelRecalibrationFactors(TObjArray *map)      { fEMCALRecalibrationFactors = map                  ;}
 151   void SetEMCALChannelRecalibrationFactors(Int_t iSM , TH2F* h) { fEMCALRecalibrationFactors->AddAt(h,iSM)          ;}
 152
 153   //-----------------------------------------------------
 154   //Modules fiducial region, remove clusters in borders
 155   //-----------------------------------------------------
 156
 157   Bool_t CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells) ;
 158   void   SetNumberOfCellsFromEMCALBorder(Int_t n) { fNCellsFromEMCALBorder = n    ;}
 159   Int_t  GetNumberOfCellsFromEMCALBorder() const  { return fNCellsFromEMCALBorder ;}
 160
 161   void   SwitchOnNoFiducialBorderInEMCALEta0()    { fNoEMCALBorderAtEta0 = kTRUE  ;}
 162   void   SwitchOffNoFiducialBorderInEMCALEta0()   { fNoEMCALBorderAtEta0 = kFALSE ;}
 163   Bool_t IsEMCALNoBorderAtEta0()                  { return fNoEMCALBorderAtEta0   ;}
 164
 165   //-----------------------------------------------------
 166   // Bad channels
 167   //-----------------------------------------------------
 168
 169   Bool_t IsBadChannelsRemovalSwitchedOn()     const { return fRemoveBadChannels       ;}
 170   void SwitchOnBadChannelsRemoval ()                { fRemoveBadChannels = kTRUE ; if(!fEMCALBadChannelMap)InitEMCALBadChannelStatusMap();}
 171   void SwitchOffBadChannelsRemoval()                { fRemoveBadChannels = kFALSE     ;}
 172
 173   Bool_t IsDistanceToBadChannelRecalculated() const { return fRecalDistToBadChannels  ;}
 174   void SwitchOnDistToBadChannelRecalculation()      { fRecalDistToBadChannels = kTRUE  ; if(!fEMCALBadChannelMap)InitEMCALBadChannelStatusMap();}
 175   void SwitchOffDistToBadChannelRecalculation()     { fRecalDistToBadChannels = kFALSE ;}
 176
 177   void InitEMCALBadChannelStatusMap() ;
 178
 179   Int_t GetEMCALChannelStatus(Int_t iSM , Int_t iCol, Int_t iRow) const {
 180     if(fEMCALBadChannelMap) return (Int_t) ((TH2I*)fEMCALBadChannelMap->At(iSM))->GetBinContent(iCol,iRow);
 181     else return 0;}//Channel is ok by default
 182
 183   void SetEMCALChannelStatus(Int_t iSM , Int_t iCol, Int_t iRow, Double_t c = 1) {
 184     if(!fEMCALBadChannelMap)InitEMCALBadChannelStatusMap() ;
 185     ((TH2I*)fEMCALBadChannelMap->At(iSM))->SetBinContent(iCol,iRow,c);}
 186
 187   TH2I * GetEMCALChannelStatusMap(Int_t iSM) const {return (TH2I*)fEMCALBadChannelMap->At(iSM);}
 188   void   SetEMCALChannelStatusMap(TObjArray *map)  {fEMCALBadChannelMap = map;}
 189   void   SetEMCALChannelStatusMap(Int_t iSM , TH2I* h) {fEMCALBadChannelMap->AddAt(h,iSM);}
 190
 191   Bool_t ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells);
 192
 193   //-----------------------------------------------------
 194   // Recalculate other cluster parameters
 195   //-----------------------------------------------------
 196
 197   void RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster);
 198   void RecalculateClusterPID(AliVCluster * cluster);
 199
 200   AliEMCALPIDUtils * GetPIDUtils() { return fPIDUtils;}
 201
 202   void RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster);
 203
 204   //----------------------------------------------------
 205   // Track matching
 206   //----------------------------------------------------
 207
 208   void    FindMatches(AliVEvent *event, TObjArray * clusterArr=0x0, AliEMCALGeometry *geom=0x0);
 209   Int_t   FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom);
 210   Bool_t  ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpEta, Float_t &tmpPhi);
 211   void    GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi);
 212   void    GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi);
 213   Int_t   GetMatchedTrackIndex(Int_t clsIndex);
 214   Int_t   GetMatchedClusterIndex(Int_t trkIndex);
 215   Bool_t  IsClusterMatched(Int_t clsIndex);
 216   Bool_t  IsTrackMatched(Int_t trkIndex);
 217   UInt_t  FindMatchedPosForCluster(Int_t clsIndex) const;
 218   UInt_t  FindMatchedPosForTrack(Int_t trkIndex) const;
 219
 220   void     SwitchOnCutEtaPhiSum()      {fCutEtaPhiSum=kTRUE;fCutEtaPhiSeparate=kFALSE;}
 221   void     SwitchOnCutEtaPhiSeparate() {fCutEtaPhiSeparate=kTRUE;fCutEtaPhiSum=kFALSE;}
 222
 223   Float_t  GetCutR()           const { return fCutR     ;}
 224   Float_t  GetCutEta()         const { return fCutEta   ;}
 225   Float_t  GetCutPhi()         const { return fCutPhi   ;}
 226   void     SetCutR(Float_t cutR)     { fCutR=cutR       ;}
 227   void     SetCutEta(Float_t cutEta) { fCutEta=cutEta   ;}
 228   void     SetCutPhi(Float_t cutPhi) { fCutPhi=cutPhi   ;}
 229   void     SetCutZ(Float_t cutZ)     {printf("Obsolete fucntion of cutZ=%1.1f\n",cutZ);} //Obsolete
 230
 231   Double_t GetMass()       const { return fMass ;}
 232   Double_t GetStep()       const { return fStep ;}
 233   void     SetMass(Double_t mass){ fMass=mass   ;}
 234   void     SetStep(Double_t step){ fStep=step   ;}
 235
 236   //Cluster cut
 237   Bool_t IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells);
 238   Bool_t IsExoticCluster(AliVCluster *cluster);
 239
 240   void SwitchOnRejectExoticCluster()  { fRejectExoticCluster=kTRUE  ;}
 241   void SwitchOffRejectExoticCluster() { fRejectExoticCluster=kFALSE ;}
 242   Bool_t IsRejectExoticCluster()      { return fRejectExoticCluster ;}
 243
 244
 245   //Track Cuts
 246   Bool_t  IsAccepted(AliESDtrack *track);
 247   void    InitTrackCuts();
 248
 249   // track quality cut setters
 250   void    SetMinTrackPt(Double_t pt=0)              { fCutMinTrackPt           = pt   ;}
 251   void    SetMinNClustersTPC(Int_t min=-1)          { fCutMinNClusterTPC       = min  ;}
 252   void    SetMinNClustersITS(Int_t min=-1)          { fCutMinNClusterITS       = min  ;}
 253   void    SetMaxChi2PerClusterTPC(Float_t max=1e10) { fCutMaxChi2PerClusterTPC = max  ;}
 254   void    SetMaxChi2PerClusterITS(Float_t max=1e10) { fCutMaxChi2PerClusterITS = max  ;}
 255   void    SetRequireTPCRefit(Bool_t b=kFALSE)       { fCutRequireTPCRefit      = b    ;}
 256   void    SetRequireITSRefit(Bool_t b=kFALSE)       { fCutRequireITSRefit      = b    ;}
 257   void    SetAcceptKinkDaughters(Bool_t b=kTRUE)    { fCutAcceptKinkDaughters  = b    ;}
 258   void    SetMaxDCAToVertexXY(Float_t dist=1e10)    { fCutMaxDCAToVertexXY     = dist ;}
 259   void    SetMaxDCAToVertexZ(Float_t dist=1e10)     { fCutMaxDCAToVertexZ      = dist ;}
 260   void    SetDCAToVertex2D(Bool_t b=kFALSE)         { fCutDCAToVertex2D        = b    ;}
 261
 262   // getters
 263   Double_t GetMinTrackPt()                  const   { return fCutMinTrackPt           ;}
 264   Int_t   GetMinNClusterTPC()               const   { return fCutMinNClusterTPC       ;}
 265   Int_t   GetMinNClustersITS()              const   { return fCutMinNClusterITS       ;}
 266   Float_t GetMaxChi2PerClusterTPC()         const   { return fCutMaxChi2PerClusterTPC ;}
 267   Float_t GetMaxChi2PerClusterITS()         const   { return fCutMaxChi2PerClusterITS ;}
 268   Bool_t  GetRequireTPCRefit()              const   { return fCutRequireTPCRefit      ;}
 269   Bool_t  GetRequireITSRefit()              const   { return fCutRequireITSRefit      ;}
 270   Bool_t  GetAcceptKinkDaughters()          const   { return fCutAcceptKinkDaughters  ;}
 271   Float_t GetMaxDCAToVertexXY()             const   { return fCutMaxDCAToVertexXY     ;}
 272   Float_t GetMaxDCAToVertexZ()              const   { return fCutMaxDCAToVertexZ      ;}
 273   Bool_t  GetDCAToVertex2D()                const   { return fCutDCAToVertex2D        ;}
 274
 275
 276 private:
 277
 278   //Position recalculation
 279   Float_t    fMisalTransShift[15];       // Shift parameters
 280   Float_t    fMisalRotShift[15];         // Shift parameters
 281   Int_t      fNonLinearityFunction;      // Non linearity function choice
 282   Float_t    fNonLinearityParams[7];     // Parameters for the non linearity function
 283   Int_t      fParticleType;              // Particle type for depth calculation
 284   Int_t      fPosAlgo;                   // Position recalculation algorithm
 285   Float_t    fW0;                        // Weight0
 286   Int_t      fNonLinearThreshold;        // Non linearity threshold value for kBeamTesh non linearity function
 287
 288   // Recalibration
 289   Bool_t     fRecalibration;             // Switch on or off the recalibration
 290   TObjArray* fEMCALRecalibrationFactors; // Array of histograms with map of recalibration factors, EMCAL
 291
 292   // Bad Channels
 293   Bool_t     fRemoveBadChannels;         // Check the channel status provided and remove clusters with bad channels
 294   Bool_t     fRecalDistToBadChannels;    // Calculate distance from highest energy tower of cluster to closes bad channel
 295   TObjArray* fEMCALBadChannelMap;        // Array of histograms with map of bad channels, EMCAL
 296
 297   // Border cells
 298   Int_t      fNCellsFromEMCALBorder;     // Number of cells from EMCAL border the cell with maximum amplitude has to be.
 299   Bool_t     fNoEMCALBorderAtEta0;       // Do fiducial cut in EMCAL region eta = 0?
 300
 301   //Track matching
 302   UInt_t     fAODFilterMask;             // Filter mask to select AOD tracks. Refer to $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
 303   TArrayI  * fMatchedTrackIndex;         // Array that stores indexes of matched tracks
 304   TArrayI  * fMatchedClusterIndex;       // Array that stores indexes of matched clusters
 305   TArrayF  * fResidualEta;               // Array that stores the residual eta
 306   TArrayF  * fResidualPhi;               // Array that stores the residual phi
 307   Bool_t     fCutEtaPhiSum;              // Place cut on sqrt(dEta^2+dPhi^2)
 308   Bool_t     fCutEtaPhiSeparate;         // Cut on dEta and dPhi separately
 309   Float_t    fCutR;                      // sqrt(dEta^2+dPhi^2) cut on matching
 310   Float_t    fCutEta;                    // dEta cut on matching
 311   Float_t    fCutPhi;                    // dPhi cut on matching
 312   Double_t   fMass;                      // Mass hypothesis of the track
 313   Double_t   fStep;                      // Length of each step used in extrapolation in the unit of cm.
 314
 315   // Cluster cuts
 316   Bool_t     fRejectExoticCluster;      // Switch on or off exotic cluster rejection
 317
 318   // Track cuts
 319   Double_t   fCutMinTrackPt;             // Cut on track pT
 320   Int_t      fCutMinNClusterTPC;         // Min number of tpc clusters
 321   Int_t      fCutMinNClusterITS;         // Min number of its clusters
 322   Float_t    fCutMaxChi2PerClusterTPC;   // Max tpc fit chi2 per tpc cluster
 323   Float_t    fCutMaxChi2PerClusterITS;   // Max its fit chi2 per its cluster
 324   Bool_t     fCutRequireTPCRefit;        // Require TPC refit
 325   Bool_t     fCutRequireITSRefit;        // Require ITS refit
 326   Bool_t     fCutAcceptKinkDaughters;    // Accepting kink daughters?
 327   Float_t    fCutMaxDCAToVertexXY;       // Track-to-vertex cut in max absolute distance in xy-plane
 328   Float_t    fCutMaxDCAToVertexZ;        // Track-to-vertex cut in max absolute distance in z-plane
 329   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
 330
 331   //PID
 332   AliEMCALPIDUtils * fPIDUtils;          // Recalculate PID parameters
 333
 334   //Time Correction
 335   Bool_t     fUseTimeCorrectionFactors;  // Use Time Dependent Correction
 336   Bool_t     fTimeCorrectionFactorsSet;  // Time Correction set at leat once
 337
 338   ClassDef(AliEMCALRecoUtils, 12)
 339
 340 };
 341
 342 #endif // ALIEMCALRECOUTILS_H
 343
 344