reverting coverity fix
[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 //          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