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