add set and getter for neutral energy fraction
[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 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