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