]> 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
change in AddTask macro for PID fluctuations (Zubayer Ahammed <Zubayer.Ahammed@cern...
[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;
1243 Short_t mclabel = -1;
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.
a840d589 1401 if (clEnergy <= 0)
1402 return;
cb231979 1403 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 1404 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
1405
d9b3567c 1406 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
094786cc 1407 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
1408 Int_t startingSM = -1;
d9b3567c 1409
7f5392da 1410 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1411 {
094786cc 1412 absId = clu->GetCellAbsId(iDig);
1413 fraction = clu->GetCellAmplitudeFraction(iDig);
1414 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
3bfc4732 1415
d9b3567c 1416 if (iDig==0) startingSM = iSupMod;
1417 else if(iSupMod != startingSM) areInSameSM = kFALSE;
094786cc 1418
1419 eCell = cells->GetCellAmplitude(absId);
d9b3567c 1420
3bfc4732 1421 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
fbddd006 1422 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
3bfc4732 1423
7f5392da 1424 if (!fCellsRecalibrated)
1425 {
1426 if(IsRecalibrationOn())
1427 {
3bfc4732 1428 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
3bfc4732 1429 }
094786cc 1430 }
3bfc4732 1431
094786cc 1432 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
d9b3567c 1433
094786cc 1434 weight = GetCellWeight(eCell,clEnergy);
d9b3567c 1435 if(weight < 0) weight = 0;
1436 totalWeight += weight;
1437 weightedCol += ieta*weight;
1438 weightedRow += iphi*weight;
1439
1440 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
841dbf60 1441 }// cell loop
094786cc 1442
d9b3567c 1443 Float_t xyzNew[]={0.,0.,0.};
7f5392da 1444 if(areInSameSM == kTRUE)
1445 {
d9b3567c 1446 //printf("In Same SM\n");
1447 weightedCol = weightedCol/totalWeight;
1448 weightedRow = weightedRow/totalWeight;
094786cc 1449 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
7f5392da 1450 }
1451 else
1452 {
d9b3567c 1453 //printf("In Different SM\n");
094786cc 1454 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 1455 }
d9b3567c 1456
094786cc 1457 clu->SetPosition(xyzNew);
d9b3567c 1458}
1459
a520bcd0 1460//___________________________________________________________________________________________
1461void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(const AliEMCALGeometry * geom,
1462 AliVCaloCells* cells,
1463 AliVCluster * cluster)
841dbf60 1464{
cb231979 1465 //re-evaluate distance to bad channel with updated bad map
1466
78467229 1467 if(!fRecalDistToBadChannels) return;
cb231979 1468
7f5392da 1469 if(!cluster)
1470 {
2aeb4226 1471 AliInfo("Cluster pointer null!");
1472 return;
1473 }
1474
841dbf60 1475 //Get channels map of the supermodule where the cluster is.
fbddd006 1476 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
cb231979 1477 Bool_t shared = kFALSE;
1478 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
1479 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
1480
fbddd006 1481 Int_t dRrow, dRcol;
841dbf60 1482 Float_t minDist = 10000.;
1483 Float_t dist = 0.;
cb231979 1484
1485 //Loop on tower status map
7f5392da 1486 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
1487 {
1488 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
1489 {
841dbf60 1490 //Check if tower is bad.
1491 if(hMap->GetBinContent(icol,irow)==0) continue;
cb231979 1492 //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 1493 // iSupMod,icol, irow, icolM,irowM);
cb231979 1494
1495 dRrow=TMath::Abs(irowM-irow);
1496 dRcol=TMath::Abs(icolM-icol);
1497 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
7f5392da 1498 if(dist < minDist)
1499 {
cb231979 1500 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
1501 minDist = dist;
1502 }
841dbf60 1503 }
1504 }
cb231979 1505
841dbf60 1506 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
7f5392da 1507 if (shared)
1508 {
841dbf60 1509 TH2D* hMap2 = 0;
1510 Int_t iSupMod2 = -1;
cb231979 1511
841dbf60 1512 //The only possible combinations are (0,1), (2,3) ... (8,9)
1513 if(iSupMod%2) iSupMod2 = iSupMod-1;
1514 else iSupMod2 = iSupMod+1;
1515 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
cb231979 1516
841dbf60 1517 //Loop on tower status map of second super module
7f5392da 1518 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++)
1519 {
1520 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++)
1521 {
841dbf60 1522 //Check if tower is bad.
1523 if(hMap2->GetBinContent(icol,irow)==0) continue;
1524 //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",
1525 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
cb231979 1526 dRrow=TMath::Abs(irow-irowM);
1527
7f5392da 1528 if(iSupMod%2)
1529 {
841dbf60 1530 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
7f5392da 1531 } else
1532 {
cb231979 1533 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
841dbf60 1534 }
cb231979 1535
841dbf60 1536 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
cb231979 1537 if(dist < minDist) minDist = dist;
841dbf60 1538 }
1539 }
1540 }// shared cluster in 2 SuperModules
78467229 1541
6fe0e6d0 1542 AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
1543 cluster->SetDistanceToBadChannel(minDist);
cb231979 1544}
1545
a520bcd0 1546//__________________________________________________________________
841dbf60 1547void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster)
1548{
83bfd77a 1549 //re-evaluate identification parameters with bayesian
2aeb4226 1550
7f5392da 1551 if(!cluster)
1552 {
2aeb4226 1553 AliInfo("Cluster pointer null!");
1554 return;
1555 }
1556
841dbf60 1557 if ( cluster->GetM02() != 0)
83bfd77a 1558 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
1559
00a38d07 1560 Float_t pidlist[AliPID::kSPECIESCN+1];
1561 for(Int_t i = 0; i < AliPID::kSPECIESCN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
841dbf60 1562
83bfd77a 1563 cluster->SetPID(pidlist);
83bfd77a 1564}
1565
f0e9e976 1566//___________________________________________________________________________________________________________________
a520bcd0 1567void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom,
1568 AliVCaloCells* cells,
f0e9e976 1569 AliVCluster * cluster,
1570 Float_t & l0, Float_t & l1,
1571 Float_t & disp, Float_t & dEta, Float_t & dPhi,
1572 Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi)
83bfd77a 1573{
1574 // Calculates new center of gravity in the local EMCAL-module coordinates
1575 // and tranfers into global ALICE coordinates
1576 // Calculates Dispersion and main axis
1577
7f5392da 1578 if(!cluster)
1579 {
2aeb4226 1580 AliInfo("Cluster pointer null!");
1581 return;
1582 }
1583
83bfd77a 1584 Double_t eCell = 0.;
1585 Float_t fraction = 1.;
1586 Float_t recalFactor = 1.;
1587
f0e9e976 1588 Int_t iSupMod = -1;
1589 Int_t iTower = -1;
1590 Int_t iIphi = -1;
1591 Int_t iIeta = -1;
1592 Int_t iphi = -1;
1593 Int_t ieta = -1;
1594 Double_t etai = -1.;
1595 Double_t phii = -1.;
1596
1597 Int_t nstat = 0 ;
1598 Float_t wtot = 0.;
1599 Double_t w = 0.;
1600 Double_t etaMean = 0.;
1601 Double_t phiMean = 0.;
83bfd77a 1602
1603 //Loop on cells
7f5392da 1604 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
1605 {
83bfd77a 1606 //Get from the absid the supermodule, tower and eta/phi numbers
1607 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1608 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1609
1610 //Get the cell energy, if recalibration is on, apply factors
1611 fraction = cluster->GetCellAmplitudeFraction(iDigit);
1612 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
3bfc4732 1613
7f5392da 1614 if (!fCellsRecalibrated)
1615 {
1616 if(IsRecalibrationOn())
1617 {
3bfc4732 1618 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1619 }
83bfd77a 1620 }
3bfc4732 1621
83bfd77a 1622 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1623
7f5392da 1624 if(cluster->E() > 0 && eCell > 0)
1625 {
83bfd77a 1626 w = GetCellWeight(eCell,cluster->E());
1627
1628 etai=(Double_t)ieta;
fbddd006 1629 phii=(Double_t)iphi;
f0e9e976 1630
7f5392da 1631 if(w > 0.0)
1632 {
83bfd77a 1633 wtot += w ;
fbddd006 1634 nstat++;
83bfd77a 1635 //Shower shape
f0e9e976 1636 sEta += w * etai * etai ;
1637 etaMean += w * etai ;
1638 sPhi += w * phii * phii ;
1639 phiMean += w * phii ;
1640 sEtaPhi += w * etai * phii ;
83bfd77a 1641 }
7f5392da 1642 }
1643 else
83bfd77a 1644 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1645 }//cell loop
1646
fbddd006 1647 //Normalize to the weight
7f5392da 1648 if (wtot > 0)
1649 {
f0e9e976 1650 etaMean /= wtot ;
1651 phiMean /= wtot ;
83bfd77a 1652 }
1653 else
1654 AliError(Form("Wrong weight %f\n", wtot));
1655
fbddd006 1656 //Calculate dispersion
7f5392da 1657 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
1658 {
83bfd77a 1659 //Get from the absid the supermodule, tower and eta/phi numbers
1660 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1661 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1662
1663 //Get the cell energy, if recalibration is on, apply factors
1664 fraction = cluster->GetCellAmplitudeFraction(iDigit);
1665 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
7f5392da 1666 if (IsRecalibrationOn())
1667 {
83bfd77a 1668 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1669 }
1670 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1671
7f5392da 1672 if(cluster->E() > 0 && eCell > 0)
1673 {
83bfd77a 1674 w = GetCellWeight(eCell,cluster->E());
1675
1676 etai=(Double_t)ieta;
fbddd006 1677 phii=(Double_t)iphi;
f0e9e976 1678 if(w > 0.0)
1679 {
1680 disp += w *((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
1681 dEta += w * (etai-etaMean)*(etai-etaMean) ;
1682 dPhi += w * (phii-phiMean)*(phii-phiMean) ;
1683 }
83bfd77a 1684 }
1685 else
1686 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1687 }// cell loop
1688
1689 //Normalize to the weigth and set shower shape parameters
7f5392da 1690 if (wtot > 0 && nstat > 1)
1691 {
f0e9e976 1692 disp /= wtot ;
1693 dEta /= wtot ;
1694 dPhi /= wtot ;
1695 sEta /= wtot ;
1696 sPhi /= wtot ;
1697 sEtaPhi /= wtot ;
1698
1699 sEta -= etaMean * etaMean ;
1700 sPhi -= phiMean * phiMean ;
1701 sEtaPhi -= etaMean * phiMean ;
1702
1703 l0 = (0.5 * (sEta + sPhi) + TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
1704 l1 = (0.5 * (sEta + sPhi) - TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
83bfd77a 1705 }
7f5392da 1706 else
1707 {
f0e9e976 1708 l0 = 0. ;
1709 l1 = 0. ;
1710 dEta = 0. ; dPhi = 0. ; disp = 0. ;
1711 sEta = 0. ; sPhi = 0. ; sEtaPhi = 0. ;
fbddd006 1712 }
83bfd77a 1713
83bfd77a 1714}
1715
f0e9e976 1716//____________________________________________________________________________________________
1717void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(const AliEMCALGeometry * geom,
1718 AliVCaloCells* cells,
1719 AliVCluster * cluster)
1720{
1721 // Calculates new center of gravity in the local EMCAL-module coordinates
1722 // and tranfers into global ALICE coordinates
1723 // Calculates Dispersion and main axis and puts them into the cluster
1724
1725 Float_t l0 = 0., l1 = 0.;
1726 Float_t disp = 0., dEta = 0., dPhi = 0.;
1727 Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
1728
1729 AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(geom,cells,cluster,l0,l1,disp,
1730 dEta, dPhi, sEta, sPhi, sEtaPhi);
1731
1732 cluster->SetM02(l0);
1733 cluster->SetM20(l1);
1734 if(disp > 0. ) cluster->SetDispersion(TMath::Sqrt(disp)) ;
1735
1736}
1737
b540d03f 1738//____________________________________________________________________________
a520bcd0 1739void AliEMCALRecoUtils::FindMatches(AliVEvent *event,
1740 TObjArray * clusterArr,
1741 const AliEMCALGeometry *geom)
bd8c7aef 1742{
1743 //This function should be called before the cluster loop
1744 //Before call this function, please recalculate the cluster positions
1745 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
1746 //Store matched cluster indexes and residuals
61160f1f 1747
88b96ad8 1748 fMatchedTrackIndex ->Reset();
bd8c7aef 1749 fMatchedClusterIndex->Reset();
fa4287a2 1750 fResidualPhi->Reset();
1751 fResidualEta->Reset();
bd8c7aef 1752
841dbf60 1753 fMatchedTrackIndex ->Set(1000);
1754 fMatchedClusterIndex->Set(1000);
1755 fResidualPhi->Set(1000);
1756 fResidualEta->Set(1000);
bd8c7aef 1757
1c7a2bf4 1758 AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1759 AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
57131575 1760
88b96ad8 1761 // init the magnetic field if not already on
7f5392da 1762 if(!TGeoGlobalMagField::Instance()->GetField())
1763 {
88b96ad8 1764 AliInfo("Init the magnetic field\n");
1765 if (esdevent)
1766 {
1767 esdevent->InitMagneticField();
1768 }
1769 else if(aodevent)
1770 {
1771 Double_t curSol = 30000*aodevent->GetMagneticField()/5.00668;
1772 Double_t curDip = 6000 *aodevent->GetMuonMagFieldScale();
1773 AliMagF *field = AliMagF::CreateFieldMap(curSol,curDip);
1774 TGeoGlobalMagField::Instance()->SetField(field);
1775 }
1776 else
1777 {
1778 AliInfo("Mag Field not initialized, null esd/aod evetn pointers");
1779 }
1780
1781 } // Init mag field
1782
42ceff04 1783 if (esdevent) {
1784 UInt_t mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ();
1785 UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco();
1786 Bool_t desc1 = (mask1 >> 3) & 0x1;
1787 Bool_t desc2 = (mask2 >> 3) & 0x1;
1788 if (desc1==0 || desc2==0) {
1789 AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
1790 mask1, esdevent->GetESDRun()->GetDetectorsInReco(),
1791 mask2, esdevent->GetESDRun()->GetDetectorsInDAQ()));
1792 fITSTrackSA=kTRUE;
1793 }
1794 }
1795
8fc351e3 1796 TObjArray *clusterArray = 0x0;
1797 if(!clusterArr)
1798 {
1799 clusterArray = new TObjArray(event->GetNumberOfCaloClusters());
1800 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
fbddd006 1801 {
1802 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1803 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1804 clusterArray->AddAt(cluster,icl);
1805 }
8fc351e3 1806 }
61160f1f 1807
bd8c7aef 1808 Int_t matched=0;
bb6f5f0b 1809 Double_t cv[21];
1810 for (Int_t i=0; i<21;i++) cv[i]=0;
bd8c7aef 1811 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1812 {
456126ad 1813 AliExternalTrackParam *trackParam = 0;
61160f1f 1814
bb6f5f0b 1815 //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
8fc351e3 1816 AliESDtrack *esdTrack = 0;
1817 AliAODTrack *aodTrack = 0;
1c7a2bf4 1818 if(esdevent)
61160f1f 1819 {
8fc351e3 1820 esdTrack = esdevent->GetTrack(itr);
1821 if(!esdTrack) continue;
1822 if(!IsAccepted(esdTrack)) continue;
61160f1f 1823 if(esdTrack->Pt()<fCutMinTrackPt) continue;
8fc351e3 1824 Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
1825 if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
42ceff04 1826 if(!fITSTrackSA)
1827 trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam()); // if TPC Available
1828 else
1829 trackParam = new AliExternalTrackParam(*esdTrack); // If ITS Track Standing alone
61160f1f 1830 }
bb6f5f0b 1831
1832 //If the input event is AOD, the starting point for extrapolation is at vertex
8fc351e3 1833 //AOD tracks are selected according to its filterbit.
1c7a2bf4 1834 else if(aodevent)
61160f1f 1835 {
8fc351e3 1836 aodTrack = aodevent->GetTrack(itr);
61160f1f 1837 if(!aodTrack) continue;
1838 if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
1839 if(aodTrack->Pt()<fCutMinTrackPt) continue;
8fc351e3 1840 Double_t phi = aodTrack->Phi()*TMath::RadToDeg();
1841 if(TMath::Abs(aodTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
61160f1f 1842 Double_t pos[3],mom[3];
1843 aodTrack->GetXYZ(pos);
1844 aodTrack->GetPxPyPz(mom);
1845 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()));
1846 trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1847 }
bd8c7aef 1848
bb6f5f0b 1849 //Return if the input data is not "AOD" or "ESD"
1850 else
61160f1f 1851 {
1852 printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
8fc351e3 1853 if(clusterArray)
fbddd006 1854 {
1855 clusterArray->Clear();
1856 delete clusterArray;
1857 }
61160f1f 1858 return;
1859 }
1860
bb6f5f0b 1861 if(!trackParam) continue;
8fc351e3 1862
1863 //Extrapolate the track to EMCal surface
1864 AliExternalTrackParam emcalParam(*trackParam);
ee602376 1865 Float_t eta, phi;
1866 if(!ExtrapolateTrackToEMCalSurface(&emcalParam, 430., fMass, fStepSurface, eta, phi))
8fc351e3 1867 {
fbddd006 1868 if(aodevent && trackParam) delete trackParam;
1869 continue;
8fc351e3 1870 }
1871
150f4870 1872// if(esdevent)
1873// {
fbddd006 1874// esdTrack->SetOuterParam(&emcalParam,AliExternalTrackParam::kMultSec);
150f4870 1875// }
8fc351e3 1876
1877 if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad())
1878 {
fbddd006 1879 if(aodevent && trackParam) delete trackParam;
1880 continue;
8fc351e3 1881 }
1882
1883
1884 //Find matched clusters
bd8c7aef 1885 Int_t index = -1;
8fc351e3 1886 Float_t dEta = -999, dPhi = -999;
1887 if(!clusterArr)
61160f1f 1888 {
fbddd006 1889 index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi);
8fc351e3 1890 }
1891 else
61160f1f 1892 {
fbddd006 1893 index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
8fc351e3 1894 }
61160f1f 1895
bd8c7aef 1896 if(index>-1)
1897 {
8fc351e3 1898 fMatchedTrackIndex ->AddAt(itr,matched);
1899 fMatchedClusterIndex ->AddAt(index,matched);
1900 fResidualEta ->AddAt(dEta,matched);
1901 fResidualPhi ->AddAt(dPhi,matched);
bd8c7aef 1902 matched++;
1903 }
456126ad 1904 if(aodevent && trackParam) delete trackParam;
bd8c7aef 1905 }//track loop
8fc351e3 1906
1907 if(clusterArray)
1908 {
1909 clusterArray->Clear();
1910 delete clusterArray;
1911 }
b540d03f 1912
1913 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1914
8fc351e3 1915 fMatchedTrackIndex ->Set(matched);
1916 fMatchedClusterIndex ->Set(matched);
1917 fResidualPhi ->Set(matched);
1918 fResidualEta ->Set(matched);
bd8c7aef 1919}
1920
b540d03f 1921//________________________________________________________________________________
a520bcd0 1922Int_t AliEMCALRecoUtils::FindMatchedClusterInEvent(const AliESDtrack *track,
1923 const AliVEvent *event,
1924 const AliEMCALGeometry *geom,
1925 Float_t &dEta, Float_t &dPhi)
bb6f5f0b 1926{
1927 //
1928 // This function returns the index of matched cluster to input track
fa4287a2 1929 // Returns -1 if no match is found
bb6f5f0b 1930 Int_t index = -1;
8fc351e3 1931 Double_t phiV = track->Phi()*TMath::RadToDeg();
1932 if(TMath::Abs(track->Eta())>0.8 || phiV <= 20 || phiV >= 240 ) return index;
42ceff04 1933 AliExternalTrackParam *trackParam = 0;
1934 if(!fITSTrackSA)
1935 trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam()); // If TPC
1936 else
1937 trackParam = new AliExternalTrackParam(*track);
1938
8fc351e3 1939 if(!trackParam) return index;
1940 AliExternalTrackParam emcalParam(*trackParam);
ee602376 1941 Float_t eta, phi;
1942 if(!ExtrapolateTrackToEMCalSurface(&emcalParam, 430., fMass, fStepSurface, eta, phi)) return index;
8fc351e3 1943 if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) return index;
1944
1945 TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters());
1946
bb6f5f0b 1947 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
61160f1f 1948 {
1949 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
8fc351e3 1950 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1951 clusterArr->AddAt(cluster,icl);
1952 }
1953
1954 index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);
1955 clusterArr->Clear();
1956 delete clusterArr;
1957
1958 return index;
1959}
1960
7f5392da 1961//_______________________________________________________________________________________________
1962Int_t AliEMCALRecoUtils::FindMatchedClusterInClusterArr(const AliExternalTrackParam *emcalParam,
a520bcd0 1963 AliExternalTrackParam *trkParam,
7f5392da 1964 const TObjArray * clusterArr,
a520bcd0 1965 Float_t &dEta, Float_t &dPhi)
8fc351e3 1966{
7f5392da 1967 // Find matched cluster in array
1968
8fc351e3 1969 dEta=-999, dPhi=-999;
1970 Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
1971 Int_t index = -1;
ee602376 1972 Float_t tmpEta=-999, tmpPhi=-999;
8fc351e3 1973
1974 Double_t exPos[3] = {0.,0.,0.};
1975 if(!emcalParam->GetXYZ(exPos)) return index;
1976
1977 Float_t clsPos[3] = {0.,0.,0.};
1978 for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
bb6f5f0b 1979 {
8fc351e3 1980 AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
1981 if(!cluster || !cluster->IsEMCAL()) continue;
1982 cluster->GetPosition(clsPos);
1983 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));
1984 if(dR > fClusterWindow) continue;
1985
1986 AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation
ee602376 1987 if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue;
8fc351e3 1988 if(fCutEtaPhiSum)
1989 {
1990 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1991 if(tmpR<dRMax)
fbddd006 1992 {
1993 dRMax=tmpR;
1994 dEtaMax=tmpEta;
1995 dPhiMax=tmpPhi;
1996 index=icl;
1997 }
8fc351e3 1998 }
1999 else if(fCutEtaPhiSeparate)
2000 {
2001 if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
fbddd006 2002 {
2003 dEtaMax = tmpEta;
2004 dPhiMax = tmpPhi;
2005 index=icl;
2006 }
8fc351e3 2007 }
2008 else
2009 {
2010 printf("Error: please specify your cut criteria\n");
2011 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
2012 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
2013 return index;
2014 }
61160f1f 2015 }
8fc351e3 2016
2017 dEta=dEtaMax;
2018 dPhi=dPhiMax;
2019
bb6f5f0b 2020 return index;
2021}
2022
88b96ad8 2023//------------------------------------------------------------------------------------
2024Bool_t AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(AliExternalTrackParam *trkParam,
a520bcd0 2025 const Double_t emcalR,
2026 const Double_t mass,
2027 const Double_t step,
88b96ad8 2028 Float_t &eta,
2029 Float_t &phi)
ee602376 2030{
88b96ad8 2031 //Extrapolate track to EMCAL surface
2032
ee602376 2033 eta = -999, phi = -999;
2034 if(!trkParam) return kFALSE;
2035 if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, emcalR, mass, step, kTRUE, 0.8, -1)) return kFALSE;
2036 Double_t trkPos[3] = {0.,0.,0.};
2037 if(!trkParam->GetXYZ(trkPos)) return kFALSE;
2038 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
2039 eta = trkPosVec.Eta();
2040 phi = trkPosVec.Phi();
2041 if(phi<0)
2042 phi += 2*TMath::Pi();
2043
2044 return kTRUE;
2045}
2046
88b96ad8 2047//-----------------------------------------------------------------------------------
2048Bool_t AliEMCALRecoUtils::ExtrapolateTrackToPosition(AliExternalTrackParam *trkParam,
2049 const Float_t *clsPos,
2050 Double_t mass,
2051 Double_t step,
2052 Float_t &tmpEta,
2053 Float_t &tmpPhi)
ee602376 2054{
2055 //
2056 //Return the residual by extrapolating a track param to a global position
2057 //
2058 tmpEta = -999;
2059 tmpPhi = -999;
2060 if(!trkParam) return kFALSE;
2061 Double_t trkPos[3] = {0.,0.,0.};
2062 TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
2063 Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
2064 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
2065 if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), mass, step,kTRUE, 0.8, -1)) return kFALSE;
2066 if(!trkParam->GetXYZ(trkPos)) return kFALSE; //Get the extrapolated global position
2067
2068 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
2069 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
2070
2071 // track cluster matching
2072 tmpPhi = clsPosVec.DeltaPhi(trkPosVec); // tmpPhi is between -pi and pi
2073 tmpEta = clsPosVec.Eta()-trkPosVec.Eta();
2074
2075 return kTRUE;
2076}
2077
88b96ad8 2078//----------------------------------------------------------------------------------
2079Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
7f5392da 2080 const AliVCluster *cluster,
a520bcd0 2081 const Double_t mass,
2082 const Double_t step,
88b96ad8 2083 Float_t &tmpEta,
2084 Float_t &tmpPhi)
ee602376 2085{
2086 //
2087 //Return the residual by extrapolating a track param to a cluster
2088 //
2089 tmpEta = -999;
2090 tmpPhi = -999;
2091 if(!cluster || !trkParam) return kFALSE;
2092
2093 Float_t clsPos[3] = {0.,0.,0.};
2094 cluster->GetPosition(clsPos);
2095
2096 return ExtrapolateTrackToPosition(trkParam, clsPos, mass, step, tmpEta, tmpPhi);
2097}
2098
88b96ad8 2099//---------------------------------------------------------------------------------
2100Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam,
7f5392da 2101 const AliVCluster *cluster,
88b96ad8 2102 Float_t &tmpEta,
2103 Float_t &tmpPhi)
bb6f5f0b 2104{
2105 //
ee602376 2106 //Return the residual by extrapolating a track param to a clusterfStepCluster
bb6f5f0b 2107 //
8fc351e3 2108
ee602376 2109 return ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, tmpEta, tmpPhi);
bb6f5f0b 2110}
2111
a520bcd0 2112//_______________________________________________________________________
2113void AliEMCALRecoUtils::GetMatchedResiduals(const Int_t clsIndex,
2114 Float_t &dEta, Float_t &dPhi)
bd8c7aef 2115{
bb6f5f0b 2116 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
fa4287a2 2117 //Get the residuals dEta and dPhi for this cluster to the closest track
bb6f5f0b 2118 //Works with ESDs and AODs
bd8c7aef 2119
bb6f5f0b 2120 if( FindMatchedPosForCluster(clsIndex) >= 999 )
bd8c7aef 2121 {
2122 AliDebug(2,"No matched tracks found!\n");
fa4287a2 2123 dEta=999.;
2124 dPhi=999.;
bd8c7aef 2125 return;
2126 }
fa4287a2 2127 dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
2128 dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
bb6f5f0b 2129}
841dbf60 2130
88b96ad8 2131//______________________________________________________________________________________________
fa4287a2 2132void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
bb6f5f0b 2133{
2134 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
fa4287a2 2135 //Get the residuals dEta and dPhi for this track to the closest cluster
bb6f5f0b 2136 //Works with ESDs and AODs
2137
2138 if( FindMatchedPosForTrack(trkIndex) >= 999 )
2139 {
2140 AliDebug(2,"No matched cluster found!\n");
fa4287a2 2141 dEta=999.;
2142 dPhi=999.;
bb6f5f0b 2143 return;
2144 }
fa4287a2 2145 dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
2146 dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
bb6f5f0b 2147}
2148
2149//__________________________________________________________
2150Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
2151{
2152 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2153 //Get the index of matched track to this cluster
2154 //Works with ESDs and AODs
2155
2156 if(IsClusterMatched(clsIndex))
2157 return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
2158 else
2159 return -1;
bd8c7aef 2160}
2161
b540d03f 2162//__________________________________________________________
bb6f5f0b 2163Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
b540d03f 2164{
bb6f5f0b 2165 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2166 //Get the index of matched cluster to this track
2167 //Works with ESDs and AODs
b540d03f 2168
bb6f5f0b 2169 if(IsTrackMatched(trkIndex))
2170 return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
b540d03f 2171 else
2172 return -1;
2173}
2174
7f5392da 2175//______________________________________________________________
7cdec71f 2176Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
bb6f5f0b 2177{
2178 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
2179 //Returns if the cluster has a match
2180 if(FindMatchedPosForCluster(clsIndex) < 999)
2181 return kTRUE;
2182 else
2183 return kFALSE;
2184}
b540d03f 2185
7f5392da 2186//____________________________________________________________
7cdec71f 2187Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const
bd8c7aef 2188{
bb6f5f0b 2189 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2190 //Returns if the track has a match
2191 if(FindMatchedPosForTrack(trkIndex) < 999)
82d09e74 2192 return kTRUE;
bd8c7aef 2193 else
2194 return kFALSE;
2195}
bb6f5f0b 2196
7f5392da 2197//______________________________________________________________________
bb6f5f0b 2198UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
bd8c7aef 2199{
bb6f5f0b 2200 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
bd8c7aef 2201 //Returns the position of the match in the fMatchedClusterIndex array
2202 Float_t tmpR = fCutR;
81efb149 2203 UInt_t pos = 999;
b540d03f 2204
7f5392da 2205 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
2206 {
2207 if(fMatchedClusterIndex->At(i)==clsIndex)
2208 {
841dbf60 2209 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
7f5392da 2210 if(r<tmpR)
2211 {
841dbf60 2212 pos=i;
2213 tmpR=r;
7f5392da 2214 AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
2215 fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
fa4287a2 2216 }
841dbf60 2217 }
bb6f5f0b 2218 }
2219 return pos;
2220}
2221
7f5392da 2222//____________________________________________________________________
bb6f5f0b 2223UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
2224{
2225 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
2226 //Returns the position of the match in the fMatchedTrackIndex array
2227 Float_t tmpR = fCutR;
2228 UInt_t pos = 999;
2229
7f5392da 2230 for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
2231 {
2232 if(fMatchedTrackIndex->At(i)==trkIndex)
2233 {
841dbf60 2234 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
7f5392da 2235 if(r<tmpR)
2236 {
841dbf60 2237 pos=i;
2238 tmpR=r;
7f5392da 2239 AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",
2240 fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
fa4287a2 2241 }
841dbf60 2242 }
b540d03f 2243 }
bd8c7aef 2244 return pos;
2245}
2246
a520bcd0 2247//__________________________________________________________________________
2248Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster,
2249 const AliEMCALGeometry *geom,
a7e5a381 2250 AliVCaloCells* cells,const Int_t bc)
b5078f5d 2251{
2252 // check if the cluster survives some quality cut
2253 //
2254 //
2255 Bool_t isGood=kTRUE;
7f5392da 2256
a7e5a381 2257 if(!cluster || !cluster->IsEMCAL()) return kFALSE;
2258
fa4287a2 2259 if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
a7e5a381 2260
fa4287a2 2261 if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
a7e5a381 2262
2263 if(IsExoticCluster(cluster, cells,bc)) return kFALSE;
b5078f5d 2264
2265 return isGood;
2266}
2267
b540d03f 2268//__________________________________________________________
bd8c7aef 2269Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
2270{
2271 // Given a esd track, return whether the track survive all the cuts
2272
2273 // The different quality parameter are first
2274 // retrieved from the track. then it is found out what cuts the
2275 // track did not survive and finally the cuts are imposed.
2276
2277 UInt_t status = esdTrack->GetStatus();
2278
2279 Int_t nClustersITS = esdTrack->GetITSclusters(0);
2280 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
2281
2282 Float_t chi2PerClusterITS = -1;
2283 Float_t chi2PerClusterTPC = -1;
2284 if (nClustersITS!=0)
2285 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
2286 if (nClustersTPC!=0)
2287 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
82d09e74 2288
2289
2290 //DCA cuts
827f9f23 2291 if(fTrackCutsType==kGlobalCut)
2292 {
2293 Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
2294 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
2295 SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
2296 }
82d09e74 2297
2298
bd8c7aef 2299 Float_t b[2];
2300 Float_t bCov[3];
2301 esdTrack->GetImpactParameters(b,bCov);
7f5392da 2302 if (bCov[0]<=0 || bCov[2]<=0)
2303 {
bd8c7aef 2304 AliDebug(1, "Estimated b resolution lower or equal zero!");
2305 bCov[0]=0; bCov[2]=0;
2306 }
2307
2308 Float_t dcaToVertexXY = b[0];
2309 Float_t dcaToVertexZ = b[1];
2310 Float_t dcaToVertex = -1;
2311
2312 if (fCutDCAToVertex2D)
2313 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
2314 else
2315 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
2316
2317 // cut the track?
2318
2319 Bool_t cuts[kNCuts];
2320 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
2321
2322 // track quality cuts
2323 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
2324 cuts[0]=kTRUE;
2325 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
2326 cuts[1]=kTRUE;
2327 if (nClustersTPC<fCutMinNClusterTPC)
2328 cuts[2]=kTRUE;
2329 if (nClustersITS<fCutMinNClusterITS)
2330 cuts[3]=kTRUE;
2331 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
2332 cuts[4]=kTRUE;
2333 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
2334 cuts[5]=kTRUE;
2335 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
2336 cuts[6]=kTRUE;
2337 if (fCutDCAToVertex2D && dcaToVertex > 1)
2338 cuts[7] = kTRUE;
2339 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
2340 cuts[8] = kTRUE;
2341 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
2342 cuts[9] = kTRUE;
2343
827f9f23 2344 if(fTrackCutsType==kGlobalCut)
2345 {
2346 //Require at least one SPD point + anything else in ITS
2347 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
fbddd006 2348 cuts[10] = kTRUE;
827f9f23 2349 }
82d09e74 2350
42ceff04 2351 // ITS
2352 if(fCutRequireITSStandAlone || fCutRequireITSpureSA){
2353 if ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin)){
2354 // TPC tracks
2355 cuts[11] = kTRUE;
2356 }else{
2357 // ITS standalone tracks
2358 if(fCutRequireITSStandAlone && !fCutRequireITSpureSA){
2359 if(status & AliESDtrack::kITSpureSA) cuts[11] = kTRUE;
2360 }else if(fCutRequireITSpureSA){
2361 if(!(status & AliESDtrack::kITSpureSA)) cuts[11] = kTRUE;
2362 }
2363 }
2364 }
2365
bd8c7aef 2366 Bool_t cut=kFALSE;
827f9f23 2367 for (Int_t i=0; i<kNCuts; i++)
7f5392da 2368 if (cuts[i]) { cut = kTRUE ; }
bd8c7aef 2369
2370 // cut the track
2371 if (cut)
2372 return kFALSE;
2373 else
2374 return kTRUE;
2375}
827f9f23 2376
88b96ad8 2377//_____________________________________
bd8c7aef 2378void AliEMCALRecoUtils::InitTrackCuts()
2379{
2380 //Intilize the track cut criteria
5f7714ad 2381 //By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
bd8c7aef 2382 //Also you can customize the cuts using the setters
82d09e74 2383
5f7714ad 2384 switch (fTrackCutsType)
88b96ad8 2385 {
5f7714ad 2386 case kTPCOnlyCut:
88b96ad8 2387 {
2388 AliInfo(Form("Track cuts for matching: GetStandardTPCOnlyTrackCuts()"));
2389 //TPC
2390 SetMinNClustersTPC(70);
2391 SetMaxChi2PerClusterTPC(4);
2392 SetAcceptKinkDaughters(kFALSE);
2393 SetRequireTPCRefit(kFALSE);
2394
2395 //ITS
2396 SetRequireITSRefit(kFALSE);
2397 SetMaxDCAToVertexZ(3.2);
2398 SetMaxDCAToVertexXY(2.4);
2399 SetDCAToVertex2D(kTRUE);
2400
2401 break;
2402 }
2403
5f7714ad 2404 case kGlobalCut:
88b96ad8 2405 {
2406 AliInfo(Form("Track cuts for matching: GetStandardITSTPCTrackCuts2010(kTURE)"));
2407 //TPC
2408 SetMinNClustersTPC(70);
2409 SetMaxChi2PerClusterTPC(4);
2410 SetAcceptKinkDaughters(kFALSE);
2411 SetRequireTPCRefit(kTRUE);
2412
2413 //ITS
2414 SetRequireITSRefit(kTRUE);
2415 SetMaxDCAToVertexZ(2);
2416 SetMaxDCAToVertexXY();
2417 SetDCAToVertex2D(kFALSE);
2418
2419 break;
2420 }
2421
0e7de35b 2422 case kLooseCut:
88b96ad8 2423 {
2424 AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
2425 SetMinNClustersTPC(50);
2426 SetAcceptKinkDaughters(kTRUE);
2427
2428 break;
5f7714ad 2429 }
42ceff04 2430
2431 case kITSStandAlone:
2432 {
2433 AliInfo(Form("Track cuts for matching: ITS Stand Alone tracks cut w/o DCA cut"));
2434 SetRequireITSRefit(kTRUE);
2435 SetRequireITSStandAlone(kTRUE);
2436 SetITSTrackSA(kTRUE);
2437 break;
2438 }
2439
88b96ad8 2440 }
bd8c7aef 2441}
83bfd77a 2442
57131575 2443
88b96ad8 2444//________________________________________________________________________
dda65b42 2445void AliEMCALRecoUtils::SetClusterMatchedToTrack(const AliVEvent *event)
57131575 2446{
2447 // Checks if tracks are matched to EMC clusters and set the matched EMCAL cluster index to ESD track.
dda65b42 2448
57131575 2449 Int_t nTracks = event->GetNumberOfTracks();
7f5392da 2450 for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack)
2451 {
dda65b42 2452 AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrack));
7f5392da 2453 if (!track)
2454 {
57131575 2455 AliWarning(Form("Could not receive track %d", iTrack));
2456 continue;
2457 }
7f5392da 2458
fbddd006 2459 Int_t matchClusIndex = GetMatchedClusterIndex(iTrack);
57131575 2460 track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
dda65b42 2461 /*the following can be done better if AliVTrack::SetStatus will be there. Patch pending with Andreas/Peter*/
2462 AliESDtrack* esdtrack = dynamic_cast<AliESDtrack*>(track);
2463 if (esdtrack) {
2464 if(matchClusIndex != -1)
2465 esdtrack->SetStatus(AliESDtrack::kEMCALmatch);
2466 else
2467 esdtrack->ResetStatus(AliESDtrack::kEMCALmatch);
2468 } else {
2469 AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
2470 if(matchClusIndex != -1)
2471 aodtrack->SetStatus(AliESDtrack::kEMCALmatch);
2472 else
2473 aodtrack->ResetStatus(AliESDtrack::kEMCALmatch);
2474 }
2475
57131575 2476 }
fbddd006 2477 AliDebug(2,"Track matched to closest cluster");
57131575 2478}
2479
88b96ad8 2480//_________________________________________________________________________
dda65b42 2481void AliEMCALRecoUtils::SetTracksMatchedToCluster(const AliVEvent *event)
57131575 2482{
2483 // Checks if EMC clusters are matched to ESD track.
2484 // Adds track indexes of all the tracks matched to a cluster withing residuals in ESDCalocluster.
2485
7f5392da 2486 for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); ++iClus)
2487 {
dda65b42 2488 AliVCluster *cluster = event->GetCaloCluster(iClus);
57131575 2489 if (!cluster->IsEMCAL())
2490 continue;
2491
2492 Int_t nTracks = event->GetNumberOfTracks();
2493 TArrayI arrayTrackMatched(nTracks);
2494
2495 // Get the closest track matched to the cluster
2496 Int_t nMatched = 0;
2497 Int_t matchTrackIndex = GetMatchedTrackIndex(iClus);
7f5392da 2498 if (matchTrackIndex != -1)
2499 {
57131575 2500 arrayTrackMatched[nMatched] = matchTrackIndex;
2501 nMatched++;
2502 }
2503
2504 // Get all other tracks matched to the cluster
7f5392da 2505 for(Int_t iTrk=0; iTrk<nTracks; ++iTrk)
2506 {
dda65b42 2507 AliVTrack* track = dynamic_cast<AliVTrack*>(event->GetTrack(iTrk));
57131575 2508 if(iTrk == matchTrackIndex) continue;
7f5392da 2509 if(track->GetEMCALcluster() == iClus)
2510 {
57131575 2511 arrayTrackMatched[nMatched] = iTrk;
2512 ++nMatched;
2513 }
2514 }
2515
2516 //printf("Tender::SetTracksMatchedToCluster - cluster E %f, N matches %d, first match %d\n",cluster->E(),nMatched,arrayTrackMatched[0]);
2517
2518 arrayTrackMatched.Set(nMatched);
dda65b42 2519 AliESDCaloCluster *esdcluster = dynamic_cast<AliESDCaloCluster*>(cluster);
2520 if (esdcluster)
2521 esdcluster->AddTracksMatched(arrayTrackMatched);
2522 else if (nMatched>0) {
2523 AliAODCaloCluster *aodcluster = dynamic_cast<AliAODCaloCluster*>(cluster);
2524 if (aodcluster)
2525 aodcluster->AddTrackMatched(event->GetTrack(arrayTrackMatched.At(0)));
2526 }
57131575 2527
2528 Float_t eta= -999, phi = -999;
2529 if (matchTrackIndex != -1)
2530 GetMatchedResiduals(iClus, eta, phi);
2531 cluster->SetTrackDistance(phi, eta);
2532 }
2533
fbddd006 2534 AliDebug(2,"Cluster matched to tracks");
57131575 2535}
2536
b540d03f 2537//___________________________________________________
d9b3567c 2538void AliEMCALRecoUtils::Print(const Option_t *) const
2539{
2540 // Print Parameters
2541
2542 printf("AliEMCALRecoUtils Settings: \n");
2543 printf("Misalignment shifts\n");
2a71e873 2544 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,
2545 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
2546 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
d9b3567c 2547 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
2548 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
094786cc 2549
2550 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
bd8c7aef 2551
fa4287a2 2552 printf("Matching criteria: ");
2553 if(fCutEtaPhiSum)
2554 {
8fc351e3 2555 printf("sqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR);
fa4287a2 2556 }
2557 else if(fCutEtaPhiSeparate)
2558 {
8fc351e3 2559 printf("dEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi);
fa4287a2 2560 }
2561 else
2562 {
2563 printf("Error\n");
2564 printf("please specify your cut criteria\n");
2565 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
2566 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
2567 }
2568
8fc351e3 2569 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);
2570 printf("Cluster selection window: dR < %2.0f\n",fClusterWindow);
bd8c7aef 2571
2572 printf("Track cuts: \n");
fa4287a2 2573 printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
bb6f5f0b 2574 printf("AOD track selection mask: %d\n",fAODFilterMask);
bd8c7aef 2575 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
2576 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
2577 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
2578 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
2579 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
d9b3567c 2580}
96957075 2581