add switch for track matching rejection for merged pi0 cluters analysis
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRecoUtils.cxx
CommitLineData
d9b3567c 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id: AliEMCALRecoUtils.cxx 33808 2009-07-15 09:48:08Z gconesab $ */
17
18///////////////////////////////////////////////////////////////////////////////
19//
20// Class AliEMCALRecoUtils
21// Some utilities to recalculate the cluster position or energy linearity
22//
23//
24// Author: Gustavo Conesa (LPSC- Grenoble)
b540d03f 25// Track matching part: Rongrong Ma (Yale)
26
d9b3567c 27///////////////////////////////////////////////////////////////////////////////
d9b3567c 28// --- standard c ---
29
30// standard C++ includes
31//#include <Riostream.h>
32
33// ROOT includes
094786cc 34#include <TGeoManager.h>
35#include <TGeoMatrix.h>
36#include <TGeoBBox.h>
7cdec71f 37#include <TH2F.h>
38#include <TArrayI.h>
39#include <TArrayF.h>
01d44f1f 40#include <TObjArray.h>
d9b3567c 41
42// STEER includes
d9b3567c 43#include "AliVCluster.h"
44#include "AliVCaloCells.h"
45#include "AliLog.h"
83bfd77a 46#include "AliPID.h"
a520bcd0 47#include "AliESDEvent.h"
bb6f5f0b 48#include "AliAODEvent.h"
bd8c7aef 49#include "AliESDtrack.h"
bb6f5f0b 50#include "AliAODTrack.h"
51#include "AliExternalTrackParam.h"
52#include "AliESDfriendTrack.h"
53#include "AliTrackerBase.h"
b540d03f 54
55// EMCAL includes
56#include "AliEMCALRecoUtils.h"
57#include "AliEMCALGeometry.h"
ee602376 58#include "AliTrackerBase.h"
b540d03f 59#include "AliEMCALPIDUtils.h"
ee602376 60
d9b3567c 61ClassImp(AliEMCALRecoUtils)
62
88b96ad8 63//_____________________________________
d9b3567c 64AliEMCALRecoUtils::AliEMCALRecoUtils():
88b96ad8 65 fParticleType(0), fPosAlgo(0), fW0(0),
66 fNonLinearityFunction(0), fNonLinearThreshold(0),
01d44f1f 67 fSmearClusterEnergy(kFALSE), fRandom(),
3bfc4732 68 fCellsRecalibrated(kFALSE), fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
7bf608c9 69 fTimeRecalibration(kFALSE), fEMCALTimeRecalibrationFactors(), fUseRunCorrectionFactors(kFALSE),
01d44f1f 70 fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
71 fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
a7e5a381 72 fRejectExoticCluster(kFALSE), fRejectExoticCells(kFALSE),
88b96ad8 73 fExoticCellFraction(0), fExoticCellDiffTime(0), fExoticCellMinAmplitude(0),
1af378e6 74 fPIDUtils(), fAODFilterMask(0), fAODHybridTracks(0),
01d44f1f 75 fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
88b96ad8 76 fResidualEta(0x0), fResidualPhi(0x0), fCutEtaPhiSum(kFALSE), fCutEtaPhiSeparate(kFALSE),
77 fCutR(0), fCutEta(0), fCutPhi(0),
78 fClusterWindow(0), fMass(0),
79 fStepSurface(0), fStepCluster(0),
42ceff04 80 fITSTrackSA(kFALSE),
88b96ad8 81 fTrackCutsType(0), fCutMinTrackPt(0), fCutMinNClusterTPC(0),
82 fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
01d44f1f 83 fCutRequireTPCRefit(kFALSE), fCutRequireITSRefit(kFALSE), fCutAcceptKinkDaughters(kFALSE),
42ceff04 84 fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0), fCutDCAToVertex2D(kFALSE),
85 fCutRequireITSStandAlone(kFALSE), fCutRequireITSpureSA(kFALSE)
d9b3567c 86{
87//
88 // Constructor.
89 // Initialize all constant values which have to be used
90 // during Reco algorithm execution
91 //
92
88b96ad8 93 // Init parameters
94 InitParameters();
01d44f1f 95
b540d03f 96 //Track matching
7cdec71f 97 fMatchedTrackIndex = new TArrayI();
98 fMatchedClusterIndex = new TArrayI();
bd36717e 99 fResidualPhi = new TArrayF();
100 fResidualEta = new TArrayF();
7cdec71f 101 fPIDUtils = new AliEMCALPIDUtils();
01d44f1f 102
d9b3567c 103}
104
105//______________________________________________________________________
106AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
01d44f1f 107: TNamed(reco),
108 fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0),
109 fNonLinearityFunction(reco.fNonLinearityFunction), fNonLinearThreshold(reco.fNonLinearThreshold),
110 fSmearClusterEnergy(reco.fSmearClusterEnergy), fRandom(),
3bfc4732 111 fCellsRecalibrated(reco.fCellsRecalibrated),
01d44f1f 112 fRecalibration(reco.fRecalibration), fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
3bfc4732 113 fTimeRecalibration(reco.fTimeRecalibration), fEMCALTimeRecalibrationFactors(reco.fEMCALTimeRecalibrationFactors),
7bf608c9 114 fUseRunCorrectionFactors(reco.fUseRunCorrectionFactors),
01d44f1f 115 fRemoveBadChannels(reco.fRemoveBadChannels), fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
78467229 116 fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
01d44f1f 117 fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder), fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
a7e5a381 118 fRejectExoticCluster(reco.fRejectExoticCluster), fRejectExoticCells(reco.fRejectExoticCells),
119 fExoticCellFraction(reco.fExoticCellFraction), fExoticCellDiffTime(reco.fExoticCellDiffTime),
120 fExoticCellMinAmplitude(reco.fExoticCellMinAmplitude),
121 fPIDUtils(reco.fPIDUtils), fAODFilterMask(reco.fAODFilterMask),
1af378e6 122 fAODHybridTracks(reco.fAODHybridTracks),
01d44f1f 123 fMatchedTrackIndex( reco.fMatchedTrackIndex? new TArrayI(*reco.fMatchedTrackIndex):0x0),
bd8c7aef 124 fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
01d44f1f 125 fResidualEta( reco.fResidualEta? new TArrayF(*reco.fResidualEta):0x0),
126 fResidualPhi( reco.fResidualPhi? new TArrayF(*reco.fResidualPhi):0x0),
127 fCutEtaPhiSum(reco.fCutEtaPhiSum), fCutEtaPhiSeparate(reco.fCutEtaPhiSeparate),
128 fCutR(reco.fCutR), fCutEta(reco.fCutEta), fCutPhi(reco.fCutPhi),
8fc351e3 129 fClusterWindow(reco.fClusterWindow),
130 fMass(reco.fMass), fStepSurface(reco.fStepSurface), fStepCluster(reco.fStepCluster),
42ceff04 131 fITSTrackSA(reco.fITSTrackSA),
01d44f1f 132 fTrackCutsType(reco.fTrackCutsType), fCutMinTrackPt(reco.fCutMinTrackPt),
133 fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS),
134 fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
135 fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
136 fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters), fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY),
42ceff04 137 fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ), fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
138 fCutRequireITSStandAlone(reco.fCutRequireITSStandAlone), fCutRequireITSpureSA(reco.fCutRequireITSpureSA)
d9b3567c 139{
140 //Copy ctor
141
01d44f1f 142 for(Int_t i = 0; i < 15 ; i++) { fMisalRotShift[i] = reco.fMisalRotShift[i] ;
143 fMisalTransShift[i] = reco.fMisalTransShift[i] ; }
144 for(Int_t i = 0; i < 7 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
145 for(Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
bd8c7aef 146
d9b3567c 147}
148
149
150//______________________________________________________________________
151AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
152{
153 //Assignment operator
154
155 if(this == &reco)return *this;
156 ((TNamed *)this)->operator=(reco);
157
01d44f1f 158 for(Int_t i = 0; i < 15 ; i++) { fMisalTransShift[i] = reco.fMisalTransShift[i] ;
159 fMisalRotShift[i] = reco.fMisalRotShift[i] ; }
160 for(Int_t i = 0; i < 7 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
161 for(Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
162
96957075 163 fParticleType = reco.fParticleType;
164 fPosAlgo = reco.fPosAlgo;
165 fW0 = reco.fW0;
01d44f1f 166
167 fNonLinearityFunction = reco.fNonLinearityFunction;
7e0ecb89 168 fNonLinearThreshold = reco.fNonLinearThreshold;
01d44f1f 169 fSmearClusterEnergy = reco.fSmearClusterEnergy;
170
3bfc4732 171 fCellsRecalibrated = reco.fCellsRecalibrated;
96957075 172 fRecalibration = reco.fRecalibration;
094786cc 173 fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
3bfc4732 174
175 fTimeRecalibration = reco.fTimeRecalibration;
176 fEMCALTimeRecalibrationFactors = reco.fEMCALTimeRecalibrationFactors;
177
178 fUseRunCorrectionFactors = reco.fUseRunCorrectionFactors;
01d44f1f 179
96957075 180 fRemoveBadChannels = reco.fRemoveBadChannels;
181 fRecalDistToBadChannels = reco.fRecalDistToBadChannels;
182 fEMCALBadChannelMap = reco.fEMCALBadChannelMap;
01d44f1f 183
96957075 184 fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
185 fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
a7e5a381 186
01d44f1f 187 fRejectExoticCluster = reco.fRejectExoticCluster;
a7e5a381 188 fRejectExoticCells = reco.fRejectExoticCells;
189 fExoticCellFraction = reco.fExoticCellFraction;
190 fExoticCellDiffTime = reco.fExoticCellDiffTime;
191 fExoticCellMinAmplitude = reco.fExoticCellMinAmplitude;
192
01d44f1f 193 fPIDUtils = reco.fPIDUtils;
83bfd77a 194
01d44f1f 195 fAODFilterMask = reco.fAODFilterMask;
1af378e6 196 fAODHybridTracks = reco.fAODHybridTracks;
d9b3567c 197
fa4287a2 198 fCutEtaPhiSum = reco.fCutEtaPhiSum;
199 fCutEtaPhiSeparate = reco.fCutEtaPhiSeparate;
96957075 200 fCutR = reco.fCutR;
fa4287a2 201 fCutEta = reco.fCutEta;
202 fCutPhi = reco.fCutPhi;
8fc351e3 203 fClusterWindow = reco.fClusterWindow;
bb6f5f0b 204 fMass = reco.fMass;
8fc351e3 205 fStepSurface = reco.fStepSurface;
206 fStepCluster = reco.fStepCluster;
42ceff04 207 fITSTrackSA = reco.fITSTrackSA;
208
5f7714ad 209 fTrackCutsType = reco.fTrackCutsType;
fa4287a2 210 fCutMinTrackPt = reco.fCutMinTrackPt;
96957075 211 fCutMinNClusterTPC = reco.fCutMinNClusterTPC;
212 fCutMinNClusterITS = reco.fCutMinNClusterITS;
213 fCutMaxChi2PerClusterTPC = reco.fCutMaxChi2PerClusterTPC;
214 fCutMaxChi2PerClusterITS = reco.fCutMaxChi2PerClusterITS;
215 fCutRequireTPCRefit = reco.fCutRequireTPCRefit;
216 fCutRequireITSRefit = reco.fCutRequireITSRefit;
217 fCutAcceptKinkDaughters = reco.fCutAcceptKinkDaughters;
218 fCutMaxDCAToVertexXY = reco.fCutMaxDCAToVertexXY;
219 fCutMaxDCAToVertexZ = reco.fCutMaxDCAToVertexZ;
220 fCutDCAToVertex2D = reco.fCutDCAToVertex2D;
42ceff04 221 fCutRequireITSStandAlone = reco.fCutRequireITSStandAlone;
222 fCutRequireITSpureSA = reco.fCutRequireITSpureSA;
7f5392da 223 if(reco.fResidualEta)
224 {
bd8c7aef 225 // assign or copy construct
7f5392da 226 if(fResidualEta)
227 {
fa4287a2 228 *fResidualEta = *reco.fResidualEta;
bd8c7aef 229 }
7f5392da 230 else
231 {
232 fResidualEta = new TArrayF(*reco.fResidualEta);
233 }
bd8c7aef 234 }
7f5392da 235 else
236 {
fa4287a2 237 if(fResidualEta)delete fResidualEta;
238 fResidualEta = 0;
bd8c7aef 239 }
240
7f5392da 241 if(reco.fResidualPhi)
242 {
bd8c7aef 243 // assign or copy construct
7f5392da 244 if(fResidualPhi)
245 {
fa4287a2 246 *fResidualPhi = *reco.fResidualPhi;
bd8c7aef 247 }
7f5392da 248 else
249 {
250 fResidualPhi = new TArrayF(*reco.fResidualPhi);
251 }
bd8c7aef 252 }
7f5392da 253 else
254 {
fa4287a2 255 if(fResidualPhi)delete fResidualPhi;
256 fResidualPhi = 0;
bd8c7aef 257 }
258
7f5392da 259 if(reco.fMatchedTrackIndex)
260 {
b540d03f 261 // assign or copy construct
7f5392da 262 if(fMatchedTrackIndex)
263 {
b540d03f 264 *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
265 }
7f5392da 266 else
267 {
268 fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
269 }
b540d03f 270 }
7f5392da 271 else
272 {
b540d03f 273 if(fMatchedTrackIndex)delete fMatchedTrackIndex;
274 fMatchedTrackIndex = 0;
275 }
bd8c7aef 276
7f5392da 277 if(reco.fMatchedClusterIndex)
278 {
bd8c7aef 279 // assign or copy construct
7f5392da 280 if(fMatchedClusterIndex)
281 {
bd8c7aef 282 *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
283 }
7f5392da 284 else
285 {
286 fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
287 }
bd8c7aef 288 }
7f5392da 289 else
290 {
bd8c7aef 291 if(fMatchedClusterIndex)delete fMatchedClusterIndex;
292 fMatchedClusterIndex = 0;
293 }
5f7714ad 294
d9b3567c 295 return *this;
296}
297
298
a7e5a381 299//_____________________________________
094786cc 300AliEMCALRecoUtils::~AliEMCALRecoUtils()
301{
302 //Destructor.
fbddd006 303
7f5392da 304 if(fEMCALRecalibrationFactors)
305 {
b6557fd1 306 fEMCALRecalibrationFactors->Clear();
841dbf60 307 delete fEMCALRecalibrationFactors;
fbddd006 308 }
fd6df01c 309
7f5392da 310 if(fEMCALTimeRecalibrationFactors)
311 {
841dbf60 312 fEMCALTimeRecalibrationFactors->Clear();
313 delete fEMCALTimeRecalibrationFactors;
fbddd006 314 }
3bfc4732 315
7f5392da 316 if(fEMCALBadChannelMap)
317 {
b6557fd1 318 fEMCALBadChannelMap->Clear();
841dbf60 319 delete fEMCALBadChannelMap;
b6557fd1 320 }
bd8c7aef 321
7cdec71f 322 delete fMatchedTrackIndex ;
323 delete fMatchedClusterIndex ;
324 delete fResidualEta ;
325 delete fResidualPhi ;
b6557fd1 326 delete fPIDUtils ;
bd8c7aef 327
b6557fd1 328 InitTrackCuts();
094786cc 329}
330
a7e5a381 331//_______________________________________________________________________________
332Bool_t AliEMCALRecoUtils::AcceptCalibrateCell(const Int_t absID, const Int_t bc,
333 Float_t & amp, Double_t & time,
334 AliVCaloCells* cells)
335{
336 // Reject cell if criteria not passed and calibrate it
337
338 AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
339
340 if(absID < 0 || absID >= 24*48*geom->GetNumberOfSuperModules()) return kFALSE;
341
342 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
c3d9f926 343
344 if(!geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta))
345 {
346 // cell absID does not exist
347 amp=0; time = 1.e9;
348 return kFALSE;
349 }
350
fbddd006 351 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
84cb7b3f 352
a7e5a381 353 // Do not include bad channels found in analysis,
c3d9f926 354 if( IsBadChannelsRemovalSwitchedOn() && GetEMCALChannelStatus(imod, ieta, iphi))
355 {
a7e5a381 356 return kFALSE;
357 }
358
359 //Recalibrate energy
360 amp = cells->GetCellAmplitude(absID);
841dbf60 361 if(!fCellsRecalibrated && IsRecalibrationOn())
a7e5a381 362 amp *= GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
363
364
365 // Recalibrate time
366 time = cells->GetCellTime(absID);
367
368 RecalibrateCellTime(absID,bc,time);
369
370 return kTRUE;
371}
372
a520bcd0 373//_____________________________________________________________________________
374Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(const AliEMCALGeometry* geom,
375 const AliVCluster* cluster,
376 AliVCaloCells* cells)
fd6df01c 377{
fbddd006 378 // Given the list of AbsId of the cluster, get the maximum cell and
379 // check if there are fNCellsFromBorder from the calorimeter border
380
7f5392da 381 if(!cluster)
382 {
2aeb4226 383 AliInfo("Cluster pointer null!");
384 return kFALSE;
385 }
386
fd6df01c 387 //If the distance to the border is 0 or negative just exit accept all clusters
fbddd006 388 if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
fd6df01c 389
fbddd006 390 Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
cb231979 391 Bool_t shared = kFALSE;
392 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
7f5392da 393
83bfd77a 394 AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
7f5392da 395 absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
fbddd006 396
397 if(absIdMax==-1) return kFALSE;
398
399 //Check if the cell is close to the borders:
400 Bool_t okrow = kFALSE;
401 Bool_t okcol = kFALSE;
fd6df01c 402
7f5392da 403 if(iSM < 0 || iphi < 0 || ieta < 0 )
404 {
fd6df01c 405 AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
406 iSM,ieta,iphi));
407 }
408
409 //Check rows/phi
7f5392da 410 if(iSM < 10)
411 {
fd6df01c 412 if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
413 }
7f5392da 414 else if (iSM >=10 && ( ( geom->GetEMCGeometry()->GetGeoName()).Contains("12SMV1")))
415 {
416 if(iphi >= fNCellsFromEMCALBorder && iphi < 8-fNCellsFromEMCALBorder) okrow =kTRUE; //1/3 sm case
417 }
418 else
419 {
420 if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE; // half SM case
421 }
fd6df01c 422
423 //Check columns/eta
7f5392da 424 if(!fNoEMCALBorderAtEta0)
425 {
fd6df01c 426 if(ieta > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
427 }
7f5392da 428 else
429 {
430 if(iSM%2==0)
431 {
fbddd006 432 if(ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
fd6df01c 433 }
7f5392da 434 else
435 {
fbddd006 436 if(ieta < 48-fNCellsFromEMCALBorder) okcol = kTRUE;
fd6df01c 437 }
438 }//eta 0 not checked
7f5392da 439
83bfd77a 440 AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
7f5392da 441 fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
fbddd006 442
443 if (okcol && okrow)
7f5392da 444 {
83bfd77a 445 //printf("Accept\n");
446 return kTRUE;
447 }
fbddd006 448 else
7f5392da 449 {
83bfd77a 450 //printf("Reject\n");
451 AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
452 return kFALSE;
453 }
fbddd006 454
455}
fd6df01c 456
457
a520bcd0 458//_______________________________________________________________________________
459Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(const AliEMCALGeometry* geom,
460 const UShort_t* cellList,
461 const Int_t nCells)
841dbf60 462{
463 // Check that in the cluster cells, there is no bad channel of those stored
464 // in fEMCALBadChannelMap or fPHOSBadChannelMap
fbddd006 465
841dbf60 466 if(!fRemoveBadChannels) return kFALSE;
467 if(!fEMCALBadChannelMap) return kFALSE;
fbddd006 468
841dbf60 469 Int_t icol = -1;
470 Int_t irow = -1;
471 Int_t imod = -1;
7f5392da 472 for(Int_t iCell = 0; iCell<nCells; iCell++)
473 {
841dbf60 474 //Get the column and row
fd6df01c 475 Int_t iTower = -1, iIphi = -1, iIeta = -1;
476 geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
477 if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
fbddd006 478 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
7f5392da 479 if(GetEMCALChannelStatus(imod, icol, irow))
480 {
83bfd77a 481 AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
482 return kTRUE;
483 }
fbddd006 484
841dbf60 485 }// cell cluster loop
fbddd006 486
841dbf60 487 return kFALSE;
fd6df01c 488}
094786cc 489
841dbf60 490
ba19aaf1 491//___________________________________________________________________________
492Float_t AliEMCALRecoUtils::GetECross(const Int_t absID, const Double_t tcell,
493 AliVCaloCells* cells, const Int_t bc)
494{
495 //Calculate the energy in the cross around the energy given cell
7f5392da 496
a7e5a381 497 AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
498
499 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
500 geom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
fbddd006 501 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
a7e5a381 502
503 //Get close cells index, energy and time, not in corners
ba19aaf1 504
84cb7b3f 505 Int_t absID1 = -1;
506 Int_t absID2 = -1;
7f5392da 507
84cb7b3f 508 if( iphi < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
509 if( iphi > 0 ) absID2 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
a7e5a381 510
7f5392da 511 // In case of cell in eta = 0 border, depending on SM shift the cross cell index
ba19aaf1 512
7f5392da 513 Int_t absID3 = -1;
514 Int_t absID4 = -1;
515
84cb7b3f 516 if ( ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2) )
7f5392da 517 {
84cb7b3f 518 absID3 = geom-> GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
519 absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
7f5392da 520 }
521 else if( ieta == 0 && imod%2 )
522 {
84cb7b3f 523 absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
524 absID4 = geom-> GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
7f5392da 525 }
526 else
527 {
84cb7b3f 528 if( ieta < AliEMCALGeoParams::fgkEMCALCols-1 )
529 absID3 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
530 if( ieta > 0 )
531 absID4 = geom-> GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
7f5392da 532 }
a7e5a381 533
ba19aaf1 534 //printf("IMOD %d, AbsId %d, a %d, b %d, c %d e %d \n",imod,absID,absID1,absID2,absID3,absID4);
a7e5a381 535
ba19aaf1 536 Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
537 Double_t tcell1 = 0, tcell2 = 0, tcell3 = 0, tcell4 = 0;
538 Bool_t accept1 = 0, accept2 = 0, accept3 = 0, accept4 = 0;
a7e5a381 539
540 accept1 = AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells);
541 accept2 = AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells);
542 accept3 = AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells);
543 accept4 = AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells);
544
545 /*
ba19aaf1 546 printf("Cell absID %d \n",absID);
547 printf("\t accept1 %d, accept2 %d, accept3 %d, accept4 %d\n",
548 accept1,accept2,accept3,accept4);
549 printf("\t id %d: id1 %d, id2 %d, id3 %d, id4 %d\n",
550 absID,absID1,absID2,absID3,absID4);
551 printf("\t e %f: e1 %f, e2 %f, e3 %f, e4 %f\n",
552 ecell,ecell1,ecell2,ecell3,ecell4);
553 printf("\t t %f: t1 %f, t2 %f, t3 %f, t4 %f;\n dt1 %f, dt2 %f, dt3 %f, dt4 %f\n",
554 tcell*1.e9,tcell1*1.e9,tcell2*1.e9,tcell3*1.e9,tcell4*1.e9,
555 TMath::Abs(tcell-tcell1)*1.e9, TMath::Abs(tcell-tcell2)*1.e9, TMath::Abs(tcell-tcell3)*1.e9, TMath::Abs(tcell-tcell4)*1.e9);
556 */
a7e5a381 557
558 if(TMath::Abs(tcell-tcell1)*1.e9 > fExoticCellDiffTime) ecell1 = 0 ;
559 if(TMath::Abs(tcell-tcell2)*1.e9 > fExoticCellDiffTime) ecell2 = 0 ;
560 if(TMath::Abs(tcell-tcell3)*1.e9 > fExoticCellDiffTime) ecell3 = 0 ;
561 if(TMath::Abs(tcell-tcell4)*1.e9 > fExoticCellDiffTime) ecell4 = 0 ;
ba19aaf1 562
563 return ecell1+ecell2+ecell3+ecell4;
564
565}
a7e5a381 566
ba19aaf1 567//_____________________________________________________________________________________________
568Bool_t AliEMCALRecoUtils::IsExoticCell(const Int_t absID, AliVCaloCells* cells, const Int_t bc)
569{
570 // Look to cell neighbourhood and reject if it seems exotic
571 // Do before recalibrating the cells
572
573 if(!fRejectExoticCells) return kFALSE;
574
575 Float_t ecell = 0;
576 Double_t tcell = 0;
577 Bool_t accept = AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells);
578
579 if(!accept) return kTRUE; // reject this cell
580
581 if(ecell < fExoticCellMinAmplitude) return kFALSE; // do not reject low energy cells
a7e5a381 582
ba19aaf1 583 Float_t eCross = GetECross(absID,tcell,cells,bc);
a7e5a381 584
7f5392da 585 if(1-eCross/ecell > fExoticCellFraction)
586 {
841dbf60 587 AliDebug(2,Form("AliEMCALRecoUtils::IsExoticCell() - EXOTIC CELL id %d, eCell %f, eCross %f, 1-eCross/eCell %f\n",
588 absID,ecell,eCross,1-eCross/ecell));
a7e5a381 589 return kTRUE;
590 }
841dbf60 591
a7e5a381 592 return kFALSE;
a7e5a381 593}
594
a520bcd0 595//___________________________________________________________________
596Bool_t AliEMCALRecoUtils::IsExoticCluster(const AliVCluster *cluster,
597 AliVCaloCells *cells,
598 const Int_t bc)
841dbf60 599{
a7e5a381 600 // Check if the cluster highest energy tower is exotic
2aeb4226 601
7f5392da 602 if(!cluster)
603 {
2aeb4226 604 AliInfo("Cluster pointer null!");
605 return kFALSE;
606 }
45516c1f 607
a7e5a381 608 if(!fRejectExoticCluster) return kFALSE;
45516c1f 609
a7e5a381 610 // Get highest energy tower
611 AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
612 Int_t iSupMod = -1, absId = -1, ieta = -1, iphi = -1;
613 Bool_t shared = kFALSE;
614 GetMaxEnergyCell(geom, cells, cluster, absId, iSupMod, ieta, iphi, shared);
7f5392da 615
a7e5a381 616 return IsExoticCell(absId,cells,bc);
7f5392da 617
b5078f5d 618}
619
a520bcd0 620//_______________________________________________________________________
841dbf60 621Float_t AliEMCALRecoUtils::SmearClusterEnergy(const AliVCluster* cluster)
622{
01d44f1f 623 //In case of MC analysis, smear energy to match resolution/calibration in real data
624
7f5392da 625 if(!cluster)
626 {
01d44f1f 627 AliInfo("Cluster pointer null!");
628 return 0;
629 }
630
631 Float_t energy = cluster->E() ;
632 Float_t rdmEnergy = energy ;
7f5392da 633 if(fSmearClusterEnergy)
634 {
01d44f1f 635 rdmEnergy = fRandom.Gaus(energy,fSmearClusterParam[0] * TMath::Sqrt(energy) +
636 fSmearClusterParam[1] * energy +
637 fSmearClusterParam[2] );
638 AliDebug(2, Form("Energy: original %f, smeared %f\n", energy, rdmEnergy));
639 }
640
841dbf60 641 return rdmEnergy;
01d44f1f 642}
643
a520bcd0 644//____________________________________________________________________________
841dbf60 645Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster)
646{
647 // Correct cluster energy from non linearity functions
2aeb4226 648
7f5392da 649 if(!cluster)
650 {
2aeb4226 651 AliInfo("Cluster pointer null!");
652 return 0;
653 }
654
d9b3567c 655 Float_t energy = cluster->E();
57131575 656
24f680c7 657 if(energy < 0.05)
5890e234 658 {
24f680c7 659 // Clusters with less than 50 MeV or negative are not possible
660 AliInfo(Form("Too Low Cluster energy!, E = %f < 0.05 GeV",energy));
5890e234 661 return 0;
662 }
663
7f5392da 664 switch (fNonLinearityFunction)
665 {
d9b3567c 666
667 case kPi0MC:
871aee7a 668 {
d9b3567c 669 //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
4439cfaf 670 //fNonLinearityParams[0] = 1.014;
671 //fNonLinearityParams[1] =-0.03329;
672 //fNonLinearityParams[2] =-0.3853;
673 //fNonLinearityParams[3] = 0.5423;
674 //fNonLinearityParams[4] =-0.4335;
8cdd1f1f 675 energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
d9b3567c 676 ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
677 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
678 break;
871aee7a 679 }
dff9e2e3 680
4439cfaf 681 case kPi0MCv2:
682 {
683 //Non-Linearity correction (from MC with function [0]/((x+[1])^[2]))+1;
684 //fNonLinearityParams[0] = 3.11111e-02;
685 //fNonLinearityParams[1] =-5.71666e-02;
686 //fNonLinearityParams[2] = 5.67995e-01;
687
688 energy *= fNonLinearityParams[0]/TMath::Power(energy+fNonLinearityParams[1],fNonLinearityParams[2])+1;
689 break;
690 }
e33e0c4a 691
692 case kPi0MCv3:
693 {
694 //Same as beam test corrected, change parameters
695 //fNonLinearityParams[0] = 9.81039e-01
696 //fNonLinearityParams[1] = 1.13508e-01;
697 //fNonLinearityParams[2] = 1.00173e+00;
698 //fNonLinearityParams[3] = 9.67998e-02;
699 //fNonLinearityParams[4] = 2.19381e+02;
700 //fNonLinearityParams[5] = 6.31604e+01;
701 //fNonLinearityParams[6] = 1;
702 energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
703
704 break;
705 }
706
4439cfaf 707
d9b3567c 708 case kPi0GammaGamma:
871aee7a 709 {
d9b3567c 710 //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
4439cfaf 711 //fNonLinearityParams[0] = 1.04;
712 //fNonLinearityParams[1] = -0.1445;
713 //fNonLinearityParams[2] = 1.046;
d9b3567c 714 energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
715 break;
871aee7a 716 }
d9b3567c 717
718 case kPi0GammaConversion:
871aee7a 719 {
d9b3567c 720 //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
871aee7a 721 //fNonLinearityParams[0] = 0.139393/0.1349766;
722 //fNonLinearityParams[1] = 0.0566186;
723 //fNonLinearityParams[2] = 0.982133;
d9b3567c 724 energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
725
726 break;
871aee7a 727 }
728
729 case kBeamTest:
730 {
731 //From beam test, Alexei's results, for different ZS thresholds
732 // th=30 MeV; th = 45 MeV; th = 75 MeV
96957075 733 //fNonLinearityParams[0] = 1.007; 1.003; 1.002
871aee7a 734 //fNonLinearityParams[1] = 0.894; 0.719; 0.797
735 //fNonLinearityParams[2] = 0.246; 0.334; 0.358
96957075 736 //Rescale the param[0] with 1.03
871aee7a 737 energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
738
739 break;
740 }
dff9e2e3 741
4b58ac4f 742 case kBeamTestCorrected:
743 {
744 //From beam test, corrected for material between beam and EMCAL
dff9e2e3 745 //fNonLinearityParams[0] = 0.99078
746 //fNonLinearityParams[1] = 0.161499;
747 //fNonLinearityParams[2] = 0.655166;
748 //fNonLinearityParams[3] = 0.134101;
749 //fNonLinearityParams[4] = 163.282;
750 //fNonLinearityParams[5] = 23.6904;
751 //fNonLinearityParams[6] = 0.978;
752 energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
753
4b58ac4f 754 break;
755 }
d9b3567c 756
757 case kNoCorrection:
758 AliDebug(2,"No correction on the energy\n");
759 break;
760
761 }
57131575 762
d9b3567c 763 return energy;
d9b3567c 764}
841dbf60 765
d9b3567c 766//__________________________________________________
7e0ecb89 767void AliEMCALRecoUtils::InitNonLinearityParam()
768{
841dbf60 769 //Initialising Non Linearity Parameters
fbddd006 770
7f5392da 771 if(fNonLinearityFunction == kPi0MC)
772 {
841dbf60 773 fNonLinearityParams[0] = 1.014;
774 fNonLinearityParams[1] = -0.03329;
775 fNonLinearityParams[2] = -0.3853;
776 fNonLinearityParams[3] = 0.5423;
777 fNonLinearityParams[4] = -0.4335;
778 }
779
e33e0c4a 780 if(fNonLinearityFunction == kPi0MCv2)
781 {
782 fNonLinearityParams[0] = 3.11111e-02;
783 fNonLinearityParams[1] =-5.71666e-02;
784 fNonLinearityParams[2] = 5.67995e-01;
785 }
786
787 if(fNonLinearityFunction == kPi0MCv3)
788 {
789 fNonLinearityParams[0] = 9.81039e-01;
790 fNonLinearityParams[1] = 1.13508e-01;
791 fNonLinearityParams[2] = 1.00173e+00;
792 fNonLinearityParams[3] = 9.67998e-02;
793 fNonLinearityParams[4] = 2.19381e+02;
794 fNonLinearityParams[5] = 6.31604e+01;
795 fNonLinearityParams[6] = 1;
796 }
797
7f5392da 798 if(fNonLinearityFunction == kPi0GammaGamma)
799 {
841dbf60 800 fNonLinearityParams[0] = 1.04;
801 fNonLinearityParams[1] = -0.1445;
802 fNonLinearityParams[2] = 1.046;
fbddd006 803 }
841dbf60 804
7f5392da 805 if(fNonLinearityFunction == kPi0GammaConversion)
806 {
841dbf60 807 fNonLinearityParams[0] = 0.139393;
808 fNonLinearityParams[1] = 0.0566186;
809 fNonLinearityParams[2] = 0.982133;
fbddd006 810 }
841dbf60 811
7f5392da 812 if(fNonLinearityFunction == kBeamTest)
813 {
814 if(fNonLinearThreshold == 30)
815 {
841dbf60 816 fNonLinearityParams[0] = 1.007;
817 fNonLinearityParams[1] = 0.894;
818 fNonLinearityParams[2] = 0.246;
819 }
7f5392da 820 if(fNonLinearThreshold == 45)
821 {
841dbf60 822 fNonLinearityParams[0] = 1.003;
823 fNonLinearityParams[1] = 0.719;
824 fNonLinearityParams[2] = 0.334;
825 }
7f5392da 826 if(fNonLinearThreshold == 75)
827 {
841dbf60 828 fNonLinearityParams[0] = 1.002;
829 fNonLinearityParams[1] = 0.797;
830 fNonLinearityParams[2] = 0.358;
831 }
832 }
833
7f5392da 834 if(fNonLinearityFunction == kBeamTestCorrected)
835 {
841dbf60 836 fNonLinearityParams[0] = 0.99078;
837 fNonLinearityParams[1] = 0.161499;
838 fNonLinearityParams[2] = 0.655166;
839 fNonLinearityParams[3] = 0.134101;
840 fNonLinearityParams[4] = 163.282;
841 fNonLinearityParams[5] = 23.6904;
842 fNonLinearityParams[6] = 0.978;
843 }
7e0ecb89 844}
845
a520bcd0 846//_________________________________________________________
847Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy,
848 const Int_t iParticle,
849 const Int_t iSM) const
094786cc 850{
851 //Calculate shower depth for a given cluster energy and particle type
852
853 // parameters
cb231979 854 Float_t x0 = 1.31;
094786cc 855 Float_t ecr = 8;
856 Float_t depth = 0;
a840d589 857 Float_t arg = energy*1000/ ecr; //Multiply energy by 1000 to transform to MeV
094786cc 858
859 switch ( iParticle )
860 {
861 case kPhoton:
a840d589 862 if (arg < 1)
863 depth = 0;
864 else
865 depth = x0 * (TMath::Log(arg) + 0.5);
094786cc 866 break;
867
868 case kElectron:
a840d589 869 if (arg < 1)
870 depth = 0;
871 else
872 depth = x0 * (TMath::Log(arg) - 0.5);
094786cc 873 break;
874
875 case kHadron:
876 // hadron
877 // boxes anc. here
7f5392da 878 if(gGeoManager)
879 {
094786cc 880 gGeoManager->cd("ALIC_1/XEN1_1");
881 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
882 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
7f5392da 883 if(geoSM)
884 {
fd6df01c 885 TGeoVolume *geoSMVol = geoSM->GetVolume();
886 TGeoShape *geoSMShape = geoSMVol->GetShape();
887 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
888 if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
889 else AliFatal("Null GEANT box");
7f5392da 890 }
891 else AliFatal("NULL GEANT node matrix");
094786cc 892 }
7f5392da 893 else
894 {//electron
a840d589 895 if (arg < 1)
896 depth = 0;
897 else
898 depth = x0 * (TMath::Log(arg) - 0.5);
094786cc 899 }
900
901 break;
902
903 default://photon
a840d589 904 if (arg < 1)
905 depth = 0;
906 else
907 depth = x0 * (TMath::Log(arg) + 0.5);
094786cc 908 }
909
910 return depth;
094786cc 911}
912
88b96ad8 913//____________________________________________________________________
914void AliEMCALRecoUtils::GetMaxEnergyCell(const AliEMCALGeometry *geom,
915 AliVCaloCells* cells,
916 const AliVCluster* clu,
a520bcd0 917 Int_t & absId,
918 Int_t & iSupMod,
919 Int_t & ieta,
920 Int_t & iphi,
921 Bool_t & shared)
d9b3567c 922{
923 //For a given CaloCluster gets the absId of the cell
924 //with maximum energy deposit.
925
926 Double_t eMax = -1.;
927 Double_t eCell = -1.;
094786cc 928 Float_t fraction = 1.;
929 Float_t recalFactor = 1.;
d9b3567c 930 Int_t cellAbsId = -1 ;
094786cc 931
d9b3567c 932 Int_t iTower = -1;
933 Int_t iIphi = -1;
934 Int_t iIeta = -1;
cb231979 935 Int_t iSupMod0= -1;
2aeb4226 936
7f5392da 937 if(!clu)
938 {
2aeb4226 939 AliInfo("Cluster pointer null!");
940 absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1;
941 return;
942 }
943
7f5392da 944 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
945 {
094786cc 946 cellAbsId = clu->GetCellAbsId(iDig);
947 fraction = clu->GetCellAmplitudeFraction(iDig);
83bfd77a 948 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
094786cc 949 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
cb231979 950 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
951 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
7f5392da 952 if (iDig==0)
953 {
954 iSupMod0=iSupMod;
955 }
956 else if(iSupMod0!=iSupMod)
957 {
cb231979 958 shared = kTRUE;
959 //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
960 }
7f5392da 961 if(!fCellsRecalibrated && IsRecalibrationOn())
962 {
094786cc 963 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
964 }
965 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
83bfd77a 966 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
7f5392da 967 if(eCell > eMax)
968 {
d9b3567c 969 eMax = eCell;
970 absId = cellAbsId;
971 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
972 }
973 }// cell loop
974
975 //Get from the absid the supermodule, tower and eta/phi numbers
976 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
977 //Gives SuperModule and Tower numbers
978 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
83bfd77a 979 iIphi, iIeta,iphi,ieta);
980 //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
981 //printf("Max end---\n");
d9b3567c 982}
983
88b96ad8 984//______________________________________
985void AliEMCALRecoUtils::InitParameters()
986{
987 // Initialize data members with default values
988
989 fParticleType = kPhoton;
990 fPosAlgo = kUnchanged;
991 fW0 = 4.5;
992
993 fNonLinearityFunction = kNoCorrection;
994 fNonLinearThreshold = 30;
995
996 fExoticCellFraction = 0.97;
997 fExoticCellDiffTime = 1e6;
998 fExoticCellMinAmplitude = 0.5;
999
1af378e6 1000 fAODFilterMask = 32;
1001 fAODHybridTracks = kFALSE;
88b96ad8 1002
1003 fCutEtaPhiSum = kTRUE;
1004 fCutEtaPhiSeparate = kFALSE;
1005
1006 fCutR = 0.05;
1007 fCutEta = 0.025;
1008 fCutPhi = 0.05;
1009
1010 fClusterWindow = 100;
1011 fMass = 0.139;
1012
1013 fStepSurface = 20.;
1014 fStepCluster = 5.;
1015 fTrackCutsType = kLooseCut;
1016
1017 fCutMinTrackPt = 0;
1018 fCutMinNClusterTPC = -1;
1019 fCutMinNClusterITS = -1;
1020
1021 fCutMaxChi2PerClusterTPC = 1e10;
1022 fCutMaxChi2PerClusterITS = 1e10;
1023
1024 fCutRequireTPCRefit = kFALSE;
1025 fCutRequireITSRefit = kFALSE;
1026 fCutAcceptKinkDaughters = kFALSE;
1027
1028 fCutMaxDCAToVertexXY = 1e10;
1029 fCutMaxDCAToVertexZ = 1e10;
1030 fCutDCAToVertex2D = kFALSE;
1031
42ceff04 1032 fCutRequireITSStandAlone = kFALSE; //MARCEL
1033 fCutRequireITSpureSA = kFALSE; //Marcel
88b96ad8 1034
1035 //Misalignment matrices
7f5392da 1036 for(Int_t i = 0; i < 15 ; i++)
1037 {
88b96ad8 1038 fMisalTransShift[i] = 0.;
1039 fMisalRotShift[i] = 0.;
1040 }
1041
1042 //Non linearity
1043 for(Int_t i = 0; i < 7 ; i++) fNonLinearityParams[i] = 0.;
1044
1045 //For kBeamTestCorrected case, but default is no correction
1046 fNonLinearityParams[0] = 0.99078;
1047 fNonLinearityParams[1] = 0.161499;
1048 fNonLinearityParams[2] = 0.655166;
1049 fNonLinearityParams[3] = 0.134101;
1050 fNonLinearityParams[4] = 163.282;
1051 fNonLinearityParams[5] = 23.6904;
1052 fNonLinearityParams[6] = 0.978;
1053
1054 //For kPi0GammaGamma case
1055 //fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
1056 //fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
1057 //fNonLinearityParams[2] = 1.046;
1058
1059 //Cluster energy smearing
1060 fSmearClusterEnergy = kFALSE;
1061 fSmearClusterParam[0] = 0.07; // * sqrt E term
1062 fSmearClusterParam[1] = 0.00; // * E term
1063 fSmearClusterParam[2] = 0.00; // constant
88b96ad8 1064}
1065
1066//_____________________________________________________
1067void AliEMCALRecoUtils::InitEMCALRecalibrationFactors()
1068{
841dbf60 1069 //Init EMCAL recalibration factors
1070 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
fbddd006 1071 //In order to avoid rewriting the same histograms
841dbf60 1072 Bool_t oldStatus = TH1::AddDirectoryStatus();
1073 TH1::AddDirectory(kFALSE);
1074
c3d9f926 1075 fEMCALRecalibrationFactors = new TObjArray(12);
1076 for (int i = 0; i < 12; i++)
841dbf60 1077 fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),
1078 Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
1079 //Init the histograms with 1
7f5392da 1080 for (Int_t sm = 0; sm < 12; sm++)
1081 {
1082 for (Int_t i = 0; i < 48; i++)
1083 {
1084 for (Int_t j = 0; j < 24; j++)
1085 {
841dbf60 1086 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
1087 }
1088 }
1089 }
7f5392da 1090
841dbf60 1091 fEMCALRecalibrationFactors->SetOwner(kTRUE);
1092 fEMCALRecalibrationFactors->Compress();
fbddd006 1093
841dbf60 1094 //In order to avoid rewriting the same histograms
fbddd006 1095 TH1::AddDirectory(oldStatus);
094786cc 1096}
1097
a520bcd0 1098//_________________________________________________________
841dbf60 1099void AliEMCALRecoUtils::InitEMCALTimeRecalibrationFactors()
1100{
1101 //Init EMCAL recalibration factors
1102 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
1103 //In order to avoid rewriting the same histograms
1104 Bool_t oldStatus = TH1::AddDirectoryStatus();
1105 TH1::AddDirectory(kFALSE);
1106
1107 fEMCALTimeRecalibrationFactors = new TObjArray(4);
1108 for (int i = 0; i < 4; i++)
3bfc4732 1109 fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvBC%d",i),
1110 Form("hAllTimeAvBC%d",i),
c3d9f926 1111 48*24*12,0.,48*24*12) );
841dbf60 1112 //Init the histograms with 1
7f5392da 1113 for (Int_t bc = 0; bc < 4; bc++)
1114 {
c3d9f926 1115 for (Int_t i = 0; i < 48*24*12; i++)
841dbf60 1116 SetEMCALChannelTimeRecalibrationFactor(bc,i,0.);
3bfc4732 1117 }
841dbf60 1118
1119 fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
1120 fEMCALTimeRecalibrationFactors->Compress();
1121
1122 //In order to avoid rewriting the same histograms
fbddd006 1123 TH1::AddDirectory(oldStatus);
3bfc4732 1124}
094786cc 1125
a520bcd0 1126//____________________________________________________
841dbf60 1127void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()
1128{
1129 //Init EMCAL bad channels map
1130 AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
1131 //In order to avoid rewriting the same histograms
1132 Bool_t oldStatus = TH1::AddDirectoryStatus();
608c80a3 1133 TH1::AddDirectory(kFALSE);
1134
c3d9f926 1135 fEMCALBadChannelMap = new TObjArray(12);
841dbf60 1136 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
7f5392da 1137 for (int i = 0; i < 12; i++)
1138 {
841dbf60 1139 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
1140 }
1141
1142 fEMCALBadChannelMap->SetOwner(kTRUE);
1143 fEMCALBadChannelMap->Compress();
1144
1145 //In order to avoid rewriting the same histograms
fbddd006 1146 TH1::AddDirectory(oldStatus);
fd6df01c 1147}
1148
88b96ad8 1149//____________________________________________________________________________
1150void AliEMCALRecoUtils::RecalibrateClusterEnergy(const AliEMCALGeometry* geom,
1151 AliVCluster * cluster,
1152 AliVCaloCells * cells,
1153 const Int_t bc)
1154{
841dbf60 1155 // Recalibrate the cluster energy and Time, considering the recalibration map
3bfc4732 1156 // and the energy of the cells and time that compose the cluster.
1157 // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
fbddd006 1158
7f5392da 1159 if(!cluster)
1160 {
2aeb4226 1161 AliInfo("Cluster pointer null!");
1162 return;
1163 }
1164
841dbf60 1165 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1166 UShort_t * index = cluster->GetCellsAbsId() ;
1167 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1168 Int_t ncells = cluster->GetNCells();
fbddd006 1169
841dbf60 1170 //Initialize some used variables
1171 Float_t energy = 0;
1172 Int_t absId =-1;
3bfc4732 1173 Int_t icol =-1, irow =-1, imod=1;
841dbf60 1174 Float_t factor = 1, frac = 0;
3bfc4732 1175 Int_t absIdMax = -1;
1176 Float_t emax = 0;
1177
841dbf60 1178 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
7f5392da 1179 for(Int_t icell = 0; icell < ncells; icell++)
1180 {
841dbf60 1181 absId = index[icell];
1182 frac = fraction[icell];
1183 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
3bfc4732 1184
7f5392da 1185 if(!fCellsRecalibrated && IsRecalibrationOn())
1186 {
3bfc4732 1187 // Energy
1188 Int_t iTower = -1, iIphi = -1, iIeta = -1;
1189 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
1190 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
fbddd006 1191 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
3bfc4732 1192 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
1193
1194 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
1195 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
1196
1197 }
1198
841dbf60 1199 energy += cells->GetCellAmplitude(absId)*factor*frac;
3bfc4732 1200
7f5392da 1201 if(emax < cells->GetCellAmplitude(absId)*factor*frac)
1202 {
3bfc4732 1203 emax = cells->GetCellAmplitude(absId)*factor*frac;
1204 absIdMax = absId;
1205 }
841dbf60 1206 }
fbddd006 1207
3a5708cd 1208 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f \n",cluster->E(),energy));
3bfc4732 1209
3a5708cd 1210 cluster->SetE(energy);
3bfc4732 1211
bc8330e9 1212 // Recalculate time of cluster
3a5708cd 1213 Double_t timeorg = cluster->GetTOF();
bc8330e9 1214
1215 Double_t time = cells->GetCellTime(absIdMax);
3a5708cd 1216 if(!fCellsRecalibrated && IsTimeRecalibrationOn())
3a5708cd 1217 RecalibrateCellTime(absIdMax,bc,time);
3bfc4732 1218
bc8330e9 1219 cluster->SetTOF(time);
3bfc4732 1220
bc8330e9 1221 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Time before %f, after %f \n",timeorg,cluster->GetTOF()));
3bfc4732 1222}
1223
a520bcd0 1224//_____________________________________________________________
1225void AliEMCALRecoUtils::RecalibrateCells(AliVCaloCells * cells,
1226 const Int_t bc)
841dbf60 1227{
1228 // Recalibrate the cells time and energy, considering the recalibration map and the energy
3bfc4732 1229 // of the cells that compose the cluster.
1230 // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
1231
87284c04 1232 if(!IsRecalibrationOn() && !IsTimeRecalibrationOn() && !IsBadChannelsRemovalSwitchedOn()) return;
3bfc4732 1233
7f5392da 1234 if(!cells)
1235 {
3bfc4732 1236 AliInfo("Cells pointer null!");
1237 return;
1238 }
1239
fbddd006 1240 Short_t absId =-1;
a7e5a381 1241 Bool_t accept = kFALSE;
1242 Float_t ecell = 0;
1243 Double_t tcell = 0;
fbddd006 1244 Double_t ecellin = 0;
1245 Double_t tcellin = 0;
60d77596 1246 Int_t mclabel = -1;
fbddd006 1247 Double_t efrac = 0;
3bfc4732 1248
841dbf60 1249 Int_t nEMcell = cells->GetNumberOfCells() ;
7f5392da 1250 for (Int_t iCell = 0; iCell < nEMcell; iCell++)
1251 {
fbddd006 1252 cells->GetCell( iCell, absId, ecellin, tcellin, mclabel, efrac );
3bfc4732 1253
a7e5a381 1254 accept = AcceptCalibrateCell(absId, bc, ecell ,tcell ,cells);
7f5392da 1255 if(!accept)
1256 {
a7e5a381 1257 ecell = 0;
fbddd006 1258 tcell = -1;
a7e5a381 1259 }
3bfc4732 1260
1261 //Set new values
fbddd006 1262 cells->SetCell(iCell,absId,ecell, tcell, mclabel, efrac);
3bfc4732 1263 }
841dbf60 1264
1265 fCellsRecalibrated = kTRUE;
094786cc 1266}
1267
88b96ad8 1268//_______________________________________________________________________________________________________
1269void AliEMCALRecoUtils::RecalibrateCellTime(const Int_t absId, const Int_t bc, Double_t & celltime) const
7d692da6 1270{
841dbf60 1271 // Recalibrate time of cell with absID considering the recalibration map
3bfc4732 1272 // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
7d692da6 1273
7f5392da 1274 if(!fCellsRecalibrated && IsTimeRecalibrationOn() && bc >= 0)
1275 {
7d692da6 1276 celltime -= GetEMCALChannelTimeRecalibrationFactor(bc%4,absId)*1.e-9; ;
3bfc4732 1277 }
3bfc4732 1278}
1279
a520bcd0 1280//______________________________________________________________________________
1281void AliEMCALRecoUtils::RecalculateClusterPosition(const AliEMCALGeometry *geom,
1282 AliVCaloCells* cells,
1283 AliVCluster* clu)
d9b3567c 1284{
1285 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1286
7f5392da 1287 if(!clu)
1288 {
2aeb4226 1289 AliInfo("Cluster pointer null!");
1290 return;
1291 }
1292
094786cc 1293 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
1294 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
fd6df01c 1295 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
094786cc 1296}
1297
a520bcd0 1298//_____________________________________________________________________________________________
1299void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(const AliEMCALGeometry *geom,
1300 AliVCaloCells* cells,
1301 AliVCluster* clu)
094786cc 1302{
1303 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1304 // The algorithm is a copy of what is done in AliEMCALRecPoint
1305
1306 Double_t eCell = 0.;
1307 Float_t fraction = 1.;
1308 Float_t recalFactor = 1.;
1309
1310 Int_t absId = -1;
1311 Int_t iTower = -1, iIphi = -1, iIeta = -1;
1312 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
1313 Float_t weight = 0., totalWeight=0.;
1314 Float_t newPos[3] = {0,0,0};
1315 Double_t pLocal[3], pGlobal[3];
cb231979 1316 Bool_t shared = kFALSE;
1317
094786cc 1318 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
a840d589 1319 if (clEnergy <= 0)
1320 return;
cb231979 1321 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 1322 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
1323
83bfd77a 1324 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
1325
7f5392da 1326 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1327 {
841dbf60 1328 absId = clu->GetCellAbsId(iDig);
1329 fraction = clu->GetCellAmplitudeFraction(iDig);
1330 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
3bfc4732 1331
7f5392da 1332 if (!fCellsRecalibrated)
1333 {
3bfc4732 1334 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
fbddd006 1335 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
3bfc4732 1336
7f5392da 1337 if(IsRecalibrationOn())
1338 {
3bfc4732 1339 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
1340 }
094786cc 1341 }
3bfc4732 1342
094786cc 1343 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1344
1345 weight = GetCellWeight(eCell,clEnergy);
1346 totalWeight += weight;
3bfc4732 1347
094786cc 1348 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
83bfd77a 1349 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
094786cc 1350 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
83bfd77a 1351 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
1352
094786cc 1353 for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
094786cc 1354 }// cell loop
1355
7f5392da 1356 if(totalWeight>0)
1357 {
094786cc 1358 for(int i=0; i<3; i++ ) newPos[i] /= totalWeight;
1359 }
1360
094786cc 1361 //Float_t pos[]={0,0,0};
1362 //clu->GetPosition(pos);
1363 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
83bfd77a 1364 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
094786cc 1365
7f5392da 1366 if(iSupModMax > 1) //sector 1
1367 {
841dbf60 1368 newPos[0] +=fMisalTransShift[3];//-=3.093;
1369 newPos[1] +=fMisalTransShift[4];//+=6.82;
1370 newPos[2] +=fMisalTransShift[5];//+=1.635;
83bfd77a 1371 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
7f5392da 1372 } else //sector 0
1373 {
841dbf60 1374 newPos[0] +=fMisalTransShift[0];//+=1.134;
1375 newPos[1] +=fMisalTransShift[1];//+=8.2;
1376 newPos[2] +=fMisalTransShift[2];//+=1.197;
83bfd77a 1377 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
841dbf60 1378 }
83bfd77a 1379 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1380
094786cc 1381 clu->SetPosition(newPos);
094786cc 1382}
1383
a520bcd0 1384//____________________________________________________________________________________________
1385void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(const AliEMCALGeometry *geom,
1386 AliVCaloCells* cells,
1387 AliVCluster* clu)
094786cc 1388{
1389 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1390 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
1391
1392 Double_t eCell = 1.;
1393 Float_t fraction = 1.;
1394 Float_t recalFactor = 1.;
1395
1396 Int_t absId = -1;
d9b3567c 1397 Int_t iTower = -1;
094786cc 1398 Int_t iIphi = -1, iIeta = -1;
841dbf60 1399 Int_t iSupMod = -1, iSupModMax = -1;
d9b3567c 1400 Int_t iphi = -1, ieta =-1;
cb231979 1401 Bool_t shared = kFALSE;
1402
d9b3567c 1403 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
a1c0f640 1404
a840d589 1405 if (clEnergy <= 0)
1406 return;
cb231979 1407 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 1408 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
1409
d9b3567c 1410 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
094786cc 1411 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
1412 Int_t startingSM = -1;
d9b3567c 1413
7f5392da 1414 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1415 {
094786cc 1416 absId = clu->GetCellAbsId(iDig);
1417 fraction = clu->GetCellAmplitudeFraction(iDig);
1418 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
3bfc4732 1419
d9b3567c 1420 if (iDig==0) startingSM = iSupMod;
1421 else if(iSupMod != startingSM) areInSameSM = kFALSE;
094786cc 1422
1423 eCell = cells->GetCellAmplitude(absId);
d9b3567c 1424
3bfc4732 1425 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
fbddd006 1426 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
3bfc4732 1427
7f5392da 1428 if (!fCellsRecalibrated)
1429 {
1430 if(IsRecalibrationOn())
1431 {
3bfc4732 1432 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
3bfc4732 1433 }
094786cc 1434 }
3bfc4732 1435
094786cc 1436 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
d9b3567c 1437
094786cc 1438 weight = GetCellWeight(eCell,clEnergy);
d9b3567c 1439 if(weight < 0) weight = 0;
1440 totalWeight += weight;
1441 weightedCol += ieta*weight;
1442 weightedRow += iphi*weight;
1443
1444 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
841dbf60 1445 }// cell loop
094786cc 1446
d9b3567c 1447 Float_t xyzNew[]={0.,0.,0.};
7f5392da 1448 if(areInSameSM == kTRUE)
1449 {
d9b3567c 1450 //printf("In Same SM\n");
1451 weightedCol = weightedCol/totalWeight;
1452 weightedRow = weightedRow/totalWeight;
094786cc 1453 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
7f5392da 1454 }
1455 else
1456 {
d9b3567c 1457 //printf("In Different SM\n");
094786cc 1458 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 1459 }
d9b3567c 1460
094786cc 1461 clu->SetPosition(xyzNew);
d9b3567c 1462}
1463
a520bcd0 1464//___________________________________________________________________________________________
1465void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(const AliEMCALGeometry * geom,
1466 AliVCaloCells* cells,
1467 AliVCluster * cluster)
841dbf60 1468{
cb231979 1469 //re-evaluate distance to bad channel with updated bad map
1470
78467229 1471 if(!fRecalDistToBadChannels) return;
cb231979 1472
7f5392da 1473 if(!cluster)
1474 {
2aeb4226 1475 AliInfo("Cluster pointer null!");
1476 return;
1477 }
1478
841dbf60 1479 //Get channels map of the supermodule where the cluster is.
fbddd006 1480 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
cb231979 1481 Bool_t shared = kFALSE;
1482 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
1483 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
1484
fbddd006 1485 Int_t dRrow, dRcol;
841dbf60 1486 Float_t minDist = 10000.;
1487 Float_t dist = 0.;
cb231979 1488
1489 //Loop on tower status map
7f5392da 1490 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
1491 {
1492 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
1493 {
841dbf60 1494 //Check if tower is bad.
1495 if(hMap->GetBinContent(icol,irow)==0) continue;
cb231979 1496 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
6fe0e6d0 1497 // iSupMod,icol, irow, icolM,irowM);
cb231979 1498
1499 dRrow=TMath::Abs(irowM-irow);
1500 dRcol=TMath::Abs(icolM-icol);
1501 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
7f5392da 1502 if(dist < minDist)
1503 {
cb231979 1504 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
1505 minDist = dist;
1506 }
841dbf60 1507 }
1508 }
cb231979 1509
841dbf60 1510 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
7f5392da 1511 if (shared)
1512 {
841dbf60 1513 TH2D* hMap2 = 0;
1514 Int_t iSupMod2 = -1;
cb231979 1515
841dbf60 1516 //The only possible combinations are (0,1), (2,3) ... (8,9)
1517 if(iSupMod%2) iSupMod2 = iSupMod-1;
1518 else iSupMod2 = iSupMod+1;
1519 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
cb231979 1520
841dbf60 1521 //Loop on tower status map of second super module
7f5392da 1522 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
1523 {
1524 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
1525 {
841dbf60 1526 //Check if tower is bad.
1527 if(hMap2->GetBinContent(icol,irow)==0) continue;
1528 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels(shared) - \n \t Bad channel in SM %d, col %d, row %d \n \t Cluster max in SM %d, col %d, row %d\n",
1529 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
cb231979 1530 dRrow=TMath::Abs(irow-irowM);
1531
7f5392da 1532 if(iSupMod%2)
1533 {
841dbf60 1534 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
7f5392da 1535 } else
1536 {
cb231979 1537 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
841dbf60 1538 }
cb231979 1539
841dbf60 1540 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
cb231979 1541 if(dist < minDist) minDist = dist;
841dbf60 1542 }
1543 }
1544 }// shared cluster in 2 SuperModules
78467229 1545
6fe0e6d0 1546 AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
1547 cluster->SetDistanceToBadChannel(minDist);
cb231979 1548}
1549
a520bcd0 1550//__________________________________________________________________
841dbf60 1551void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster)
1552{
83bfd77a 1553 //re-evaluate identification parameters with bayesian
2aeb4226 1554
7f5392da 1555 if(!cluster)
1556 {
2aeb4226 1557 AliInfo("Cluster pointer null!");
1558 return;
1559 }
1560
841dbf60 1561 if ( cluster->GetM02() != 0)
83bfd77a 1562 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
1563
00a38d07 1564 Float_t pidlist[AliPID::kSPECIESCN+1];
1565 for(Int_t i = 0; i < AliPID::kSPECIESCN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
841dbf60 1566
83bfd77a 1567 cluster->SetPID(pidlist);
83bfd77a 1568}
1569
f0e9e976 1570//___________________________________________________________________________________________________________________
a520bcd0 1571void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom,
1572 AliVCaloCells* cells,
f0e9e976 1573 AliVCluster * cluster,
1574 Float_t & l0, Float_t & l1,
1575 Float_t & disp, Float_t & dEta, Float_t & dPhi,
1576 Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
83bfd77a 1577{
1578 // Calculates new center of gravity in the local EMCAL-module coordinates
1579 // and tranfers into global ALICE coordinates
1580 // Calculates Dispersion and main axis
1581
7f5392da 1582 if(!cluster)
1583 {
2aeb4226 1584 AliInfo("Cluster pointer null!");
1585 return;
1586 }
1587
83bfd77a 1588 Double_t eCell = 0.;
1589 Float_t fraction = 1.;
1590 Float_t recalFactor = 1.;
1591
f0e9e976 1592 Int_t iSupMod = -1;
1593 Int_t iTower = -1;
1594 Int_t iIphi = -1;
1595 Int_t iIeta = -1;
1596 Int_t iphi = -1;
1597 Int_t ieta = -1;
1598 Double_t etai = -1.;
1599 Double_t phii = -1.;
1600
1601 Int_t nstat = 0 ;
1602 Float_t wtot = 0.;
1603 Double_t w = 0.;
1604 Double_t etaMean = 0.;
1605 Double_t phiMean = 0.;
a1c0f640 1606
1607 //Loop on cells, calculate the cluster energy, in case a cut on cell energy is added
1608 // and to check if the cluster is between 2 SM in eta
1609 Int_t iSM0 = -1;
1610 Bool_t shared = kFALSE;
1611 Float_t energy = 0;
1612
1613 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
1614 {
1615 //Get from the absid the supermodule, tower and eta/phi numbers
1616 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1617 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1618
1619 //Check if there are cells of different SM
1620 if (iDigit == 0 ) iSM0 = iSupMod;
1621 else if(iSupMod!= iSM0) shared = kTRUE;
1622
1623 //Get the cell energy, if recalibration is on, apply factors
1624 fraction = cluster->GetCellAmplitudeFraction(iDigit);
1625 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
83bfd77a 1626
a1c0f640 1627 if(IsRecalibrationOn())
1628 {
1629 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1630 }
1631
1632 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1633
1634 energy += eCell;
1635
1636 }//cell loop
1637
83bfd77a 1638 //Loop on cells
7f5392da 1639 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
1640 {
83bfd77a 1641 //Get from the absid the supermodule, tower and eta/phi numbers
1642 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1643 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1644
1645 //Get the cell energy, if recalibration is on, apply factors
1646 fraction = cluster->GetCellAmplitudeFraction(iDigit);
1647 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
3bfc4732 1648
7f5392da 1649 if (!fCellsRecalibrated)
1650 {
1651 if(IsRecalibrationOn())
1652 {
3bfc4732 1653 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1654 }
83bfd77a 1655 }
3bfc4732 1656
83bfd77a 1657 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1658
a1c0f640 1659 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
1660 // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
1661 if(shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
1662
7f5392da 1663 if(cluster->E() > 0 && eCell > 0)
1664 {
83bfd77a 1665 w = GetCellWeight(eCell,cluster->E());
1666
1667 etai=(Double_t)ieta;
fbddd006 1668 phii=(Double_t)iphi;
f0e9e976 1669
7f5392da 1670 if(w > 0.0)
1671 {
83bfd77a 1672 wtot += w ;
fbddd006 1673 nstat++;
83bfd77a 1674 //Shower shape
f0e9e976 1675 sEta += w * etai * etai ;
1676 etaMean += w * etai ;
1677 sPhi += w * phii * phii ;
1678 phiMean += w * phii ;
1679 sEtaPhi += w * etai * phii ;
83bfd77a 1680 }
7f5392da 1681 }
1682 else
83bfd77a 1683 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1684 }//cell loop
1685
fbddd006 1686 //Normalize to the weight
7f5392da 1687 if (wtot > 0)
1688 {
f0e9e976 1689 etaMean /= wtot ;
1690 phiMean /= wtot ;
83bfd77a 1691 }
1692 else
1693 AliError(Form("Wrong weight %f\n", wtot));
1694
fbddd006 1695 //Calculate dispersion
7f5392da 1696 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
1697 {
83bfd77a 1698 //Get from the absid the supermodule, tower and eta/phi numbers
1699 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1700 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1701
1702 //Get the cell energy, if recalibration is on, apply factors
1703 fraction = cluster->GetCellAmplitudeFraction(iDigit);
1704 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
7f5392da 1705 if (IsRecalibrationOn())
1706 {
83bfd77a 1707 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1708 }
1709 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1710
a1c0f640 1711 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
1712 // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
1713 if(shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
1714
7f5392da 1715 if(cluster->E() > 0 && eCell > 0)
1716 {
83bfd77a 1717 w = GetCellWeight(eCell,cluster->E());
1718
1719 etai=(Double_t)ieta;
fbddd006 1720 phii=(Double_t)iphi;
f0e9e976 1721 if(w > 0.0)
1722 {
1723 disp += w *((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
1724 dEta += w * (etai-etaMean)*(etai-etaMean) ;
1725 dPhi += w * (phii-phiMean)*(phii-phiMean) ;
1726 }
83bfd77a 1727 }
1728 else
1729 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1730 }// cell loop
1731
1732 //Normalize to the weigth and set shower shape parameters
7f5392da 1733 if (wtot > 0 && nstat > 1)
1734 {
f0e9e976 1735 disp /= wtot ;
1736 dEta /= wtot ;
1737 dPhi /= wtot ;
1738 sEta /= wtot ;
1739 sPhi /= wtot ;
1740 sEtaPhi /= wtot ;
1741
1742 sEta -= etaMean * etaMean ;
1743 sPhi -= phiMean * phiMean ;
1744 sEtaPhi -= etaMean * phiMean ;
1745
1746 l0 = (0.5 * (sEta + sPhi) + TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
1747 l1 = (0.5 * (sEta + sPhi) - TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
83bfd77a 1748 }
7f5392da 1749 else
1750 {
f0e9e976 1751 l0 = 0. ;
1752 l1 = 0. ;
1753 dEta = 0. ; dPhi = 0. ; disp = 0. ;
1754 sEta = 0. ; sPhi = 0. ; sEtaPhi = 0. ;
fbddd006 1755 }
83bfd77a 1756
83bfd77a 1757}
1758
f0e9e976 1759//____________________________________________________________________________________________
1760void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom,
1761 AliVCaloCells* cells,
1762 AliVCluster * cluster)
1763{
1764 // Calculates new center of gravity in the local EMCAL-module coordinates
1765 // and tranfers into global ALICE coordinates
1766 // Calculates Dispersion and main axis and puts them into the cluster
1767
1768 Float_t l0 = 0., l1 = 0.;
1769 Float_t disp = 0., dEta = 0., dPhi = 0.;
1770 Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
1771
1772 AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(geom,cells,cluster,l0,l1,disp,
1773 dEta, dPhi, sEta, sPhi, sEtaPhi);
1774
1775 cluster->SetM02(l0);
1776 cluster->SetM20(l1);
1777 if(disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
1778
1779}
1780
b540d03f 1781//____________________________________________________________________________
a520bcd0 1782void AliEMCALRecoUtils::FindMatches(AliVEvent *event,
1783 TObjArray * clusterArr,
1784 const AliEMCALGeometry *geom)
bd8c7aef 1785{
1786 //This function should be called before the cluster loop
1787 //Before call this function, please recalculate the cluster positions
1788 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
1789 //Store matched cluster indexes and residuals
61160f1f 1790
88b96ad8 1791 fMatchedTrackIndex ->Reset();
bd8c7aef 1792 fMatchedClusterIndex->Reset();
fa4287a2 1793 fResidualPhi->Reset();
1794 fResidualEta->Reset();
bd8c7aef 1795
841dbf60 1796 fMatchedTrackIndex ->Set(1000);
1797 fMatchedClusterIndex->Set(1000);
1798 fResidualPhi->Set(1000);
1799 fResidualEta->Set(1000);
bd8c7aef 1800
1c7a2bf4 1801 AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1802 AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
57131575 1803
88b96ad8 1804 // init the magnetic field if not already on
7f5392da 1805 if(!TGeoGlobalMagField::Instance()->GetField())
1806 {
88b96ad8 1807 AliInfo("Init the magnetic field\n");
1808 if (esdevent)
1809 {
1810 esdevent->InitMagneticField();
1811 }
1812 else if(aodevent)
1813 {
1814 Double_t curSol = 30000*aodevent->GetMagneticField()/5.00668;
1815 Double_t curDip = 6000 *aodevent->GetMuonMagFieldScale();
1816 AliMagF *field = AliMagF::CreateFieldMap(curSol,curDip);
1817 TGeoGlobalMagField::Instance()->SetField(field);
1818 }
1819 else
1820 {
1821 AliInfo("Mag Field not initialized, null esd/aod evetn pointers");
1822 }
1823
1824 } // Init mag field
1825
42ceff04 1826 if (esdevent) {
1827 UInt_t mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ();
1828 UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco();
1829 Bool_t desc1 = (mask1 >> 3) & 0x1;
1830 Bool_t desc2 = (mask2 >> 3) & 0x1;
1831 if (desc1==0 || desc2==0) {
30e29b2a 1832// AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
1833// mask1, esdevent->GetESDRun()->GetDetectorsInReco(),
1834// mask2, esdevent->GetESDRun()->GetDetectorsInDAQ()));
42ceff04 1835 fITSTrackSA=kTRUE;
1836 }
1837 }
1838
8fc351e3 1839 TObjArray *clusterArray = 0x0;
1840 if(!clusterArr)
1841 {
1842 clusterArray = new TObjArray(event->GetNumberOfCaloClusters());
1843 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
fbddd006 1844 {
1845 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1846 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1847 clusterArray->AddAt(cluster,icl);
1848 }
8fc351e3 1849 }
61160f1f 1850
bd8c7aef 1851 Int_t matched=0;
bb6f5f0b 1852 Double_t cv[21];
1853 for (Int_t i=0; i<21;i++) cv[i]=0;
bd8c7aef 1854 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1855 {
456126ad 1856 AliExternalTrackParam *trackParam = 0;
61160f1f 1857
bb6f5f0b 1858 //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
8fc351e3 1859 AliESDtrack *esdTrack = 0;
1860 AliAODTrack *aodTrack = 0;
1c7a2bf4 1861 if(esdevent)
61160f1f 1862 {
8fc351e3 1863 esdTrack = esdevent->GetTrack(itr);
1864 if(!esdTrack) continue;
1865 if(!IsAccepted(esdTrack)) continue;
61160f1f 1866 if(esdTrack->Pt()<fCutMinTrackPt) continue;
8fc351e3 1867 Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
1868 if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
42ceff04 1869 if(!fITSTrackSA)
1870 trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam()); // if TPC Available
1871 else
1872 trackParam = new AliExternalTrackParam(*esdTrack); // If ITS Track Standing alone
61160f1f 1873 }
bb6f5f0b 1874
1875 //If the input event is AOD, the starting point for extrapolation is at vertex
8fc351e3 1876 //AOD tracks are selected according to its filterbit.
1c7a2bf4 1877 else if(aodevent)
61160f1f 1878 {
8fc351e3 1879 aodTrack = aodevent->GetTrack(itr);
61160f1f 1880 if(!aodTrack) continue;
1af378e6 1881
1882 //Check mask if not hybrid
1883 if(!fAODHybridTracks && !aodTrack->TestFilterMask(fAODFilterMask) ) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
1884
1885 //Check hybrid
1886 if( fAODHybridTracks && !aodTrack->IsHybridGlobalConstrainedGlobal()) continue ;
1887
61160f1f 1888 if(aodTrack->Pt()<fCutMinTrackPt) continue;
1af378e6 1889
8fc351e3 1890 Double_t phi = aodTrack->Phi()*TMath::RadToDeg();
1891 if(TMath::Abs(aodTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
61160f1f 1892 Double_t pos[3],mom[3];
1893 aodTrack->GetXYZ(pos);
1894 aodTrack->GetPxPyPz(mom);
1895 AliDebug(5,Form("aod track: i=%d | pos=(%5.4f,%5.4f,%5.4f) | mom=(%5.4f,%5.4f,%5.4f) | charge=%d\n",itr,pos[0],pos[1],pos[2],mom[0],mom[1],mom[2],aodTrack->Charge()));
1af378e6 1896
61160f1f 1897 trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1898 }
bd8c7aef 1899
bb6f5f0b 1900 //Return if the input data is not "AOD" or "ESD"
1901 else
61160f1f 1902 {
1903 printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
8fc351e3 1904 if(clusterArray)
fbddd006 1905 {
1906 clusterArray->Clear();
1907 delete clusterArray;
1908 }
61160f1f 1909 return;
1910 }
1911
bb6f5f0b 1912 if(!trackParam) continue;
8fc351e3 1913
1914 //Extrapolate the track to EMCal surface
1915 AliExternalTrackParam emcalParam(*trackParam);
ee602376 1916 Float_t eta, phi;
1917 if(!ExtrapolateTrackToEMCalSurface(&emcalParam, 430., fMass, fStepSurface, eta, phi))
8fc351e3 1918 {
fbddd006 1919 if(aodevent && trackParam) delete trackParam;
30e29b2a 1920 if(fITSTrackSA && trackParam) delete trackParam;
fbddd006 1921 continue;
8fc351e3 1922 }
1923
150f4870 1924// if(esdevent)
1925// {
fbddd006 1926// esdTrack->SetOuterParam(&emcalParam,AliExternalTrackParam::kMultSec);
150f4870 1927// }
8fc351e3 1928
1929 if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad())
1930 {
fbddd006 1931 if(aodevent && trackParam) delete trackParam;
30e29b2a 1932 if(fITSTrackSA && trackParam) delete trackParam;
fbddd006 1933 continue;
8fc351e3 1934 }
1935
1936
1937 //Find matched clusters
bd8c7aef 1938 Int_t index = -1;
8fc351e3 1939 Float_t dEta = -999, dPhi = -999;
1940 if(!clusterArr)
61160f1f 1941 {
fbddd006 1942 index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi);
8fc351e3 1943 }
1944 else
61160f1f 1945 {
fbddd006 1946 index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
8fc351e3 1947 }
61160f1f 1948
bd8c7aef 1949 if(index>-1)
1950 {
8fc351e3 1951 fMatchedTrackIndex ->AddAt(itr,matched);
1952 fMatchedClusterIndex ->AddAt(index,matched);
1953 fResidualEta ->AddAt(dEta,matched);
1954 fResidualPhi ->AddAt(dPhi,matched);
bd8c7aef 1955 matched++;
1956 }
456126ad 1957 if(aodevent && trackParam) delete trackParam;
30e29b2a 1958 if(fITSTrackSA && trackParam) delete trackParam;
bd8c7aef 1959 }//track loop
8fc351e3 1960
1961 if(clusterArray)
1962 {
1963 clusterArray->Clear();
1964 delete clusterArray;
1965 }
b540d03f 1966
1967 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1968
8fc351e3 1969 fMatchedTrackIndex ->Set(matched);
1970 fMatchedClusterIndex ->Set(matched);
1971 fResidualPhi ->Set(matched);
1972 fResidualEta ->Set(matched);
bd8c7aef 1973}
1974
b540d03f 1975//________________________________________________________________________________
a520bcd0 1976Int_t AliEMCALRecoUtils::FindMatchedClusterInEvent(const AliESDtrack *track,
1977 const AliVEvent *event,
1978 const AliEMCALGeometry *geom,
1979 Float_t &dEta, Float_t &dPhi)
bb6f5f0b 1980{
1981 //
1982 // This function returns the index of matched cluster to input track
fa4287a2 1983 // Returns -1 if no match is found
bb6f5f0b 1984 Int_t index = -1;
8fc351e3 1985 Double_t phiV = track->Phi()*TMath::RadToDeg();
1986 if(TMath::Abs(track->Eta())>0.8 || phiV <= 20 || phiV >= 240 ) return index;
42ceff04 1987 AliExternalTrackParam *trackParam = 0;
1988 if(!fITSTrackSA)
1989 trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam()); // If TPC
1990 else
1991 trackParam = new AliExternalTrackParam(*track);
1992
8fc351e3 1993 if(!trackParam) return index;
1994 AliExternalTrackParam emcalParam(*trackParam);
ee602376 1995 Float_t eta, phi;
8fc351e3 1996
30e29b2a 1997 if(!ExtrapolateTrackToEMCalSurface(&emcalParam, 430., fMass, fStepSurface, eta, phi)) {
1998 if(fITSTrackSA) delete trackParam;
1999 return index;
2000 }
2001 if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()){
2002 if(fITSTrackSA) delete trackParam;
2003 return index;
2004 }
2005
8fc351e3 2006 TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters());
2007
bb6f5f0b 2008 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
61160f1f 2009 {
2010 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
8fc351e3 2011 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
2012 clusterArr->AddAt(cluster,icl);
2013 }
2014
2015 index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
2016 clusterArr->Clear();
2017 delete clusterArr;
30e29b2a 2018 if(fITSTrackSA) delete trackParam;
2019
8fc351e3 2020 return index;
2021}
2022
7f5392da 2023//_______________________________________________________________________________________________
2024Int_t AliEMCALRecoUtils::FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam,
a520bcd0 2025 AliExternalTrackParam *trkParam,
7f5392da 2026 const TObjArray * clusterArr,
a520bcd0 2027 Float_t &dEta, Float_t &dPhi)
8fc351e3 2028{
7f5392da 2029 // Find matched cluster in array
2030
8fc351e3 2031 dEta=-999, dPhi=-999;
2032 Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
2033 Int_t index = -1;
ee602376 2034 Float_t tmpEta=-999, tmpPhi=-999;
8fc351e3 2035
2036 Double_t exPos[3] = {0.,0.,0.};
2037 if(!emcalParam->GetXYZ(exPos)) return index;
2038
2039 Float_t clsPos[3] = {0.,0.,0.};
2040 for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
bb6f5f0b 2041 {
8fc351e3 2042 AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
2043 if(!cluster || !cluster->IsEMCAL()) continue;
2044 cluster->GetPosition(clsPos);
2045 Double_t dR = TMath::Sqrt(TMath::Power(exPos[0]-clsPos[0],2)+TMath::Power(exPos[1]-clsPos[1],2)+TMath::Power(exPos[2]-clsPos[2],2));
2046 if(dR > fClusterWindow) continue;
2047
2048 AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation
ee602376 2049 if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue;
8fc351e3 2050 if(fCutEtaPhiSum)
2051 {
2052 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
2053 if(tmpR<dRMax)
fbddd006 2054 {
2055 dRMax=tmpR;
2056 dEtaMax=tmpEta;
2057 dPhiMax=tmpPhi;
2058 index=icl;
2059 }
8fc351e3 2060 }
2061 else if(fCutEtaPhiSeparate)
2062 {
2063 if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
fbddd006 2064 {
2065 dEtaMax = tmpEta;
2066 dPhiMax = tmpPhi;
2067 index=icl;
2068 }
8fc351e3 2069 }
2070 else
2071 {
2072 printf("Error: please specify your cut criteria\n");
2073 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
2074 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
2075 return index;
2076 }
61160f1f 2077 }
8fc351e3 2078
2079 dEta=dEtaMax;
2080 dPhi=dPhiMax;
2081
bb6f5f0b 2082 return index;
2083}
2084
88b96ad8 2085//------------------------------------------------------------------------------------
2086Bool_t AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(AliExternalTrackParam *trkParam,
a520bcd0 2087 const Double_t emcalR,
2088 const Double_t mass,
2089 const Double_t step,
88b96ad8 2090 Float_t &eta,
2091 Float_t &phi)
ee602376 2092{
88b96ad8 2093 //Extrapolate track to EMCAL surface
2094
ee602376 2095 eta = -999, phi = -999;
2096 if(!trkParam) return kFALSE;
2097 if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, emcalR, mass, step, kTRUE, 0.8, -1)) return kFALSE;
2098 Double_t trkPos[3] = {0.,0.,0.};
2099 if(!trkParam->GetXYZ(trkPos)) return kFALSE;
2100 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
2101 eta = trkPosVec.Eta();
2102 phi = trkPosVec.Phi();
2103 if(phi<0)
2104 phi += 2*TMath::Pi();
2105
2106 return kTRUE;
2107}
2108
88b96ad8 2109//-----------------------------------------------------------------------------------
2110Bool_t AliEMCALRecoUtils::ExtrapolateTrackToPosition(AliExternalTrackParam *trkParam,
2111 const Float_t *clsPos,
2112 Double_t mass,
2113 Double_t step,
2114 Float_t &tmpEta,
2115 Float_t &tmpPhi)
ee602376 2116{
2117 //
2118 //Return the residual by extrapolating a track param to a global position
2119 //
2120 tmpEta = -999;
2121 tmpPhi = -999;
2122 if(!trkParam) return kFALSE;
2123 Double_t trkPos[3] = {0.,0.,0.};
2124 TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
2125 Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
2126 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
2127 if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), mass, step,kTRUE, 0.8, -1)) return kFALSE;
2128 if(!trkParam->GetXYZ(trkPos)) return kFALSE; //Get the extrapolated global position
2129
2130 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
2131 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
2132
2133 // track cluster matching
2134 tmpPhi = clsPosVec.DeltaPhi(trkPosVec); // tmpPhi is between -pi and pi
2135 tmpEta = clsPosVec.Eta()-trkPosVec.Eta();
2136
2137 return kTRUE;
2138}
2139
88b96ad8 2140//----------------------------------------------------------------------------------
2141Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
7f5392da 2142 const AliVCluster *cluster,
a520bcd0 2143 const Double_t mass,
2144 const Double_t step,
88b96ad8 2145 Float_t &tmpEta,
2146 Float_t &tmpPhi)
ee602376 2147{
2148 //
2149 //Return the residual by extrapolating a track param to a cluster
2150 //
2151 tmpEta = -999;
2152 tmpPhi = -999;
2153 if(!cluster || !trkParam) return kFALSE;
2154
2155 Float_t clsPos[3] = {0.,0.,0.};
2156 cluster->GetPosition(clsPos);
2157
2158 return ExtrapolateTrackToPosition(trkParam, clsPos, mass, step, tmpEta, tmpPhi);
2159}
2160
88b96ad8 2161//---------------------------------------------------------------------------------
2162Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
7f5392da 2163 const AliVCluster *cluster,
88b96ad8 2164 Float_t &tmpEta,
2165 Float_t &tmpPhi)
bb6f5f0b 2166{
2167 //
ee602376 2168 //Return the residual by extrapolating a track param to a clusterfStepCluster
bb6f5f0b 2169 //
8fc351e3 2170
ee602376 2171 return ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, tmpEta, tmpPhi);
bb6f5f0b 2172}
2173
a520bcd0 2174//_______________________________________________________________________
2175void AliEMCALRecoUtils::GetMatchedResiduals(const Int_t clsIndex,
2176 Float_t &dEta, Float_t &dPhi)
bd8c7aef 2177{
bb6f5f0b 2178 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
fa4287a2 2179 //Get the residuals dEta and dPhi for this cluster to the closest track
bb6f5f0b 2180 //Works with ESDs and AODs
bd8c7aef 2181
bb6f5f0b 2182 if( FindMatchedPosForCluster(clsIndex) >= 999 )
bd8c7aef 2183 {
2184 AliDebug(2,"No matched tracks found!\n");
fa4287a2 2185 dEta=999.;
2186 dPhi=999.;
bd8c7aef 2187 return;
2188 }
fa4287a2 2189 dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
2190 dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
bb6f5f0b 2191}
841dbf60 2192
88b96ad8 2193//______________________________________________________________________________________________
fa4287a2 2194void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
bb6f5f0b 2195{
2196 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
fa4287a2 2197 //Get the residuals dEta and dPhi for this track to the closest cluster
bb6f5f0b 2198 //Works with ESDs and AODs
2199
2200 if( FindMatchedPosForTrack(trkIndex) >= 999 )
2201 {
2202 AliDebug(2,"No matched cluster found!\n");
fa4287a2 2203 dEta=999.;
2204 dPhi=999.;
bb6f5f0b 2205 return;
2206 }
fa4287a2 2207 dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
2208 dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
bb6f5f0b 2209}
2210
2211//__________________________________________________________
2212Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
2213{
2214 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2215 //Get the index of matched track to this cluster
2216 //Works with ESDs and AODs
2217
2218 if(IsClusterMatched(clsIndex))
2219 return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
2220 else
2221 return -1;
bd8c7aef 2222}
2223
b540d03f 2224//__________________________________________________________
bb6f5f0b 2225Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
b540d03f 2226{
bb6f5f0b 2227 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2228 //Get the index of matched cluster to this track
2229 //Works with ESDs and AODs
b540d03f 2230
bb6f5f0b 2231 if(IsTrackMatched(trkIndex))
2232 return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
b540d03f 2233 else
2234 return -1;
2235}
2236
7f5392da 2237//______________________________________________________________
7cdec71f 2238Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
bb6f5f0b 2239{
2240 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2241 //Returns if the cluster has a match
2242 if(FindMatchedPosForCluster(clsIndex) < 999)
2243 return kTRUE;
2244 else
2245 return kFALSE;
2246}
b540d03f 2247
7f5392da 2248//____________________________________________________________
7cdec71f 2249Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const
bd8c7aef 2250{
bb6f5f0b 2251 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2252 //Returns if the track has a match
2253 if(FindMatchedPosForTrack(trkIndex) < 999)
82d09e74 2254 return kTRUE;
bd8c7aef 2255 else
2256 return kFALSE;
2257}
bb6f5f0b 2258
7f5392da 2259//______________________________________________________________________
bb6f5f0b 2260UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
bd8c7aef 2261{
bb6f5f0b 2262 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
bd8c7aef 2263 //Returns the position of the match in the fMatchedClusterIndex array
2264 Float_t tmpR = fCutR;
81efb149 2265 UInt_t pos = 999;
b540d03f 2266
7f5392da 2267 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
2268 {
2269 if(fMatchedClusterIndex->At(i)==clsIndex)
2270 {
841dbf60 2271 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
7f5392da 2272 if(r<tmpR)
2273 {
841dbf60 2274 pos=i;
2275 tmpR=r;
7f5392da 2276 AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
2277 fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
fa4287a2 2278 }
841dbf60 2279 }
bb6f5f0b 2280 }
2281 return pos;
2282}
2283
7f5392da 2284//____________________________________________________________________
bb6f5f0b 2285UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
2286{
2287 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2288 //Returns the position of the match in the fMatchedTrackIndex array
2289 Float_t tmpR = fCutR;
2290 UInt_t pos = 999;
2291
7f5392da 2292 for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
2293 {
2294 if(fMatchedTrackIndex->At(i)==trkIndex)
2295 {
841dbf60 2296 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
7f5392da 2297 if(r<tmpR)
2298 {
841dbf60 2299 pos=i;
2300 tmpR=r;
7f5392da 2301 AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
2302 fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
fa4287a2 2303 }
841dbf60 2304 }
b540d03f 2305 }
bd8c7aef 2306 return pos;
2307}
2308
a520bcd0 2309//__________________________________________________________________________
2310Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster,
2311 const AliEMCALGeometry *geom,
a7e5a381 2312 AliVCaloCells* cells,const Int_t bc)
b5078f5d 2313{
2314 // check if the cluster survives some quality cut
2315 //
2316 //
2317 Bool_t isGood=kTRUE;
7f5392da 2318
a7e5a381 2319 if(!cluster || !cluster->IsEMCAL()) return kFALSE;
2320
fa4287a2 2321 if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
a7e5a381 2322
fa4287a2 2323 if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
a7e5a381 2324
2325 if(IsExoticCluster(cluster, cells,bc)) return kFALSE;
b5078f5d 2326
2327 return isGood;
2328}
2329
2330//__________________________________________________________
bd8c7aef 2331Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
2332{
2333 // Given a esd track, return whether the track survive all the cuts
2334
2335 // The different quality parameter are first
2336 // retrieved from the track. then it is found out what cuts the
2337 // track did not survive and finally the cuts are imposed.
2338
2339 UInt_t status = esdTrack->GetStatus();
2340
2341 Int_t nClustersITS = esdTrack->GetITSclusters(0);
2342 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
2343
2344 Float_t chi2PerClusterITS = -1;
2345 Float_t chi2PerClusterTPC = -1;
2346 if (nClustersITS!=0)
2347 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
2348 if (nClustersTPC!=0)
2349 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
82d09e74 2350
2351
2352 //DCA cuts
827f9f23 2353 if(fTrackCutsType==kGlobalCut)
2354 {
2355 Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
2356 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
2357 SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
2358 }
82d09e74 2359
2360
bd8c7aef 2361 Float_t b[2];
2362 Float_t bCov[3];
2363 esdTrack->GetImpactParameters(b,bCov);
7f5392da 2364 if (bCov[0]<=0 || bCov[2]<=0)
2365 {
bd8c7aef 2366 AliDebug(1, "Estimated b resolution lower or equal zero!");
2367 bCov[0]=0; bCov[2]=0;
2368 }
2369
2370 Float_t dcaToVertexXY = b[0];
2371 Float_t dcaToVertexZ = b[1];
2372 Float_t dcaToVertex = -1;
2373
2374 if (fCutDCAToVertex2D)
2375 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
2376 else
2377 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
2378
2379 // cut the track?
2380
2381 Bool_t cuts[kNCuts];
2382 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
2383
2384 // track quality cuts
2385 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
2386 cuts[0]=kTRUE;
2387 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
2388 cuts[1]=kTRUE;
2389 if (nClustersTPC<fCutMinNClusterTPC)
2390 cuts[2]=kTRUE;
2391 if (nClustersITS<fCutMinNClusterITS)
2392 cuts[3]=kTRUE;
2393 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
2394 cuts[4]=kTRUE;
2395 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
2396 cuts[5]=kTRUE;
2397 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
2398 cuts[6]=kTRUE;
2399 if (fCutDCAToVertex2D && dcaToVertex > 1)
2400 cuts[7] = kTRUE;
2401 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
2402 cuts[8] = kTRUE;
2403 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
2404 cuts[9] = kTRUE;
2405
827f9f23 2406 if(fTrackCutsType==kGlobalCut)
2407 {
2408 //Require at least one SPD point + anything else in ITS
2409 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
fbddd006 2410 cuts[10] = kTRUE;
827f9f23 2411 }
82d09e74 2412
42ceff04 2413 // ITS
2414 if(fCutRequireITSStandAlone || fCutRequireITSpureSA){
2415 if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin)){
2416 // TPC tracks
2417 cuts[11] = kTRUE;
2418 }else{
2419 // ITS standalone tracks
2420 if(fCutRequireITSStandAlone && !fCutRequireITSpureSA){
2421 if(status & AliESDtrack::kITSpureSA) cuts[11] = kTRUE;
2422 }else if(fCutRequireITSpureSA){
2423 if(!(status & AliESDtrack::kITSpureSA)) cuts[11] = kTRUE;
2424 }
2425 }
2426 }
2427
bd8c7aef 2428 Bool_t cut=kFALSE;
827f9f23 2429 for (Int_t i=0; i<kNCuts; i++)
7f5392da 2430 if (cuts[i]) { cut = kTRUE ; }
bd8c7aef 2431
2432 // cut the track
2433 if (cut)
2434 return kFALSE;
2435 else
2436 return kTRUE;
2437}
827f9f23 2438
88b96ad8 2439//_____________________________________
bd8c7aef 2440void AliEMCALRecoUtils::InitTrackCuts()
2441{
2442 //Intilize the track cut criteria
5f7714ad 2443 //By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
bd8c7aef 2444 //Also you can customize the cuts using the setters
82d09e74 2445
5f7714ad 2446 switch (fTrackCutsType)
88b96ad8 2447 {
5f7714ad 2448 case kTPCOnlyCut:
88b96ad8 2449 {
2450 AliInfo(Form("Track cuts for matching: GetStandardTPCOnlyTrackCuts()"));
2451 //TPC
2452 SetMinNClustersTPC(70);
2453 SetMaxChi2PerClusterTPC(4);
2454 SetAcceptKinkDaughters(kFALSE);
2455 SetRequireTPCRefit(kFALSE);
2456
2457 //ITS
2458 SetRequireITSRefit(kFALSE);
2459 SetMaxDCAToVertexZ(3.2);
2460 SetMaxDCAToVertexXY(2.4);
2461 SetDCAToVertex2D(kTRUE);
2462
2463 break;
2464 }
2465
5f7714ad 2466 case kGlobalCut:
88b96ad8 2467 {
2468 AliInfo(Form("Track cuts for matching: GetStandardITSTPCTrackCuts2010(kTURE)"));
2469 //TPC
2470 SetMinNClustersTPC(70);
2471 SetMaxChi2PerClusterTPC(4);
2472 SetAcceptKinkDaughters(kFALSE);
2473 SetRequireTPCRefit(kTRUE);
2474
2475 //ITS
2476 SetRequireITSRefit(kTRUE);
2477 SetMaxDCAToVertexZ(2);
2478 SetMaxDCAToVertexXY();
2479 SetDCAToVertex2D(kFALSE);
2480
2481 break;
2482 }
2483
0e7de35b 2484 case kLooseCut:
88b96ad8 2485 {
2486 AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
2487 SetMinNClustersTPC(50);
2488 SetAcceptKinkDaughters(kTRUE);
2489
2490 break;
5f7714ad 2491 }
42ceff04 2492
2493 case kITSStandAlone:
2494 {
2495 AliInfo(Form("Track cuts for matching: ITS Stand Alone tracks cut w/o DCA cut"));
2496 SetRequireITSRefit(kTRUE);
2497 SetRequireITSStandAlone(kTRUE);
2498 SetITSTrackSA(kTRUE);
2499 break;
2500 }
2501
88b96ad8 2502 }
bd8c7aef 2503}
83bfd77a 2504
57131575 2505
88b96ad8 2506//________________________________________________________________________
dda65b42 2507void AliEMCALRecoUtils::SetClusterMatchedToTrack(const AliVEvent *event)
57131575 2508{
2509 // Checks if tracks are matched to EMC clusters and set the matched EMCAL cluster index to ESD track.
dda65b42 2510
57131575 2511 Int_t nTracks = event->GetNumberOfTracks();
7f5392da 2512 for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack)
2513 {
dda65b42 2514 AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrack));
7f5392da 2515 if (!track)
2516 {
57131575 2517 AliWarning(Form("Could not receive track %d", iTrack));
2518 continue;
2519 }
7f5392da 2520
fbddd006 2521 Int_t matchClusIndex = GetMatchedClusterIndex(iTrack);
57131575 2522 track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
dda65b42 2523 /*the following can be done better if AliVTrack::SetStatus will be there. Patch pending with Andreas/Peter*/
2524 AliESDtrack* esdtrack = dynamic_cast<AliESDtrack*>(track);
2525 if (esdtrack) {
2526 if(matchClusIndex != -1)
2527 esdtrack->SetStatus(AliESDtrack::kEMCALmatch);
2528 else
2529 esdtrack->ResetStatus(AliESDtrack::kEMCALmatch);
2530 } else {
2531 AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
2532 if(matchClusIndex != -1)
2533 aodtrack->SetStatus(AliESDtrack::kEMCALmatch);
2534 else
2535 aodtrack->ResetStatus(AliESDtrack::kEMCALmatch);
2536 }
2537
57131575 2538 }
fbddd006 2539 AliDebug(2,"Track matched to closest cluster");
57131575 2540}
2541
88b96ad8 2542//_________________________________________________________________________
dda65b42 2543void AliEMCALRecoUtils::SetTracksMatchedToCluster(const AliVEvent *event)
57131575 2544{
2545 // Checks if EMC clusters are matched to ESD track.
2546 // Adds track indexes of all the tracks matched to a cluster withing residuals in ESDCalocluster.
2547
7f5392da 2548 for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus)
2549 {
dda65b42 2550 AliVCluster *cluster = event->GetCaloCluster(iClus);
57131575 2551 if (!cluster->IsEMCAL())
2552 continue;
2553
2554 Int_t nTracks = event->GetNumberOfTracks();
2555 TArrayI arrayTrackMatched(nTracks);
2556
2557 // Get the closest track matched to the cluster
2558 Int_t nMatched = 0;
2559 Int_t matchTrackIndex = GetMatchedTrackIndex(iClus);
7f5392da 2560 if (matchTrackIndex != -1)
2561 {
57131575 2562 arrayTrackMatched[nMatched] = matchTrackIndex;
2563 nMatched++;
2564 }
2565
2566 // Get all other tracks matched to the cluster
7f5392da 2567 for(Int_t iTrk=0; iTrk<nTracks; ++iTrk)
2568 {
dda65b42 2569 AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrk));
57131575 2570 if(iTrk == matchTrackIndex) continue;
7f5392da 2571 if(track->GetEMCALcluster() == iClus)
2572 {
57131575 2573 arrayTrackMatched[nMatched] = iTrk;
2574 ++nMatched;
2575 }
2576 }
2577
2578 //printf("Tender::SetTracksMatchedToCluster - cluster E %f, N matches %d, first match %d\n",cluster->E(),nMatched,arrayTrackMatched[0]);
2579
2580 arrayTrackMatched.Set(nMatched);
dda65b42 2581 AliESDCaloCluster *esdcluster = dynamic_cast<AliESDCaloCluster*>(cluster);
2582 if (esdcluster)
2583 esdcluster->AddTracksMatched(arrayTrackMatched);
2584 else if (nMatched>0) {
2585 AliAODCaloCluster *aodcluster = dynamic_cast<AliAODCaloCluster*>(cluster);
2586 if (aodcluster)
2587 aodcluster->AddTrackMatched(event->GetTrack(arrayTrackMatched.At(0)));
2588 }
57131575 2589
2590 Float_t eta= -999, phi = -999;
2591 if (matchTrackIndex != -1)
2592 GetMatchedResiduals(iClus, eta, phi);
2593 cluster->SetTrackDistance(phi, eta);
2594 }
2595
fbddd006 2596 AliDebug(2,"Cluster matched to tracks");
57131575 2597}
2598
b540d03f 2599//___________________________________________________
d9b3567c 2600void AliEMCALRecoUtils::Print(const Option_t *) const
2601{
2602 // Print Parameters
2603
2604 printf("AliEMCALRecoUtils Settings: \n");
2605 printf("Misalignment shifts\n");
2a71e873 2606 for(Int_t i=0; i<5; i++) printf("\t sector %d, traslation (x,y,z)=(%f,%f,%f), rotation (x,y,z)=(%f,%f,%f)\n",i,
2607 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
2608 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
d9b3567c 2609 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
2610 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
094786cc 2611
2612 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
bd8c7aef 2613
fa4287a2 2614 printf("Matching criteria: ");
2615 if(fCutEtaPhiSum)
2616 {
8fc351e3 2617 printf("sqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR);
fa4287a2 2618 }
2619 else if(fCutEtaPhiSeparate)
2620 {
8fc351e3 2621 printf("dEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi);
fa4287a2 2622 }
2623 else
2624 {
2625 printf("Error\n");
2626 printf("please specify your cut criteria\n");
2627 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
2628 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
2629 }
2630
8fc351e3 2631 printf("Mass hypothesis = %2.3f [GeV/c^2], extrapolation step to surface = %2.2f[cm], step to cluster = %2.2f[cm]\n",fMass,fStepSurface, fStepCluster);
2632 printf("Cluster selection window: dR < %2.0f\n",fClusterWindow);
bd8c7aef 2633
2634 printf("Track cuts: \n");
fa4287a2 2635 printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
bb6f5f0b 2636 printf("AOD track selection mask: %d\n",fAODFilterMask);
bd8c7aef 2637 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
2638 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
2639 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
2640 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
2641 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
d9b3567c 2642}
96957075 2643