]> 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
Updating LHC data file test files: replacing LHC_Active_Injection_Scheme with LHC_Dat...
[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///////////////////////////////////////////////////////////////////////////////
28
29// --- standard c ---
30
31// standard C++ includes
32//#include <Riostream.h>
33
34// ROOT includes
094786cc 35#include <TGeoManager.h>
36#include <TGeoMatrix.h>
37#include <TGeoBBox.h>
d9b3567c 38
39// STEER includes
d9b3567c 40#include "AliVCluster.h"
41#include "AliVCaloCells.h"
bd8c7aef 42#include "AliVEvent.h"
d9b3567c 43#include "AliLog.h"
83bfd77a 44#include "AliPID.h"
bd8c7aef 45#include "AliESDEvent.h"
bb6f5f0b 46#include "AliAODEvent.h"
bd8c7aef 47#include "AliESDtrack.h"
bb6f5f0b 48#include "AliAODTrack.h"
49#include "AliExternalTrackParam.h"
50#include "AliESDfriendTrack.h"
51#include "AliTrackerBase.h"
b540d03f 52
53// EMCAL includes
54#include "AliEMCALRecoUtils.h"
55#include "AliEMCALGeometry.h"
bd8c7aef 56#include "AliEMCALTrack.h"
96957075 57#include "AliEMCALCalibTimeDepCorrection.h"
b540d03f 58#include "AliEMCALPIDUtils.h"
d9b3567c 59
60ClassImp(AliEMCALRecoUtils)
61
62//______________________________________________
63AliEMCALRecoUtils::AliEMCALRecoUtils():
fd6df01c 64 fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
65 fPosAlgo(kUnchanged),fW0(4.),
66 fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
78467229 67 fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
bd8c7aef 68 fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
bb6f5f0b 69 fAODFilterMask(32),
b540d03f 70 fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
bb6f5f0b 71 fResidualZ(0x0), fResidualR(0x0), fCutR(10), fCutZ(10), fMass(0.139), fStep(1),
b5078f5d 72 fRejectExoticCluster(kFALSE),
bd8c7aef 73 fCutMinNClusterTPC(0), fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
74 fCutRequireTPCRefit(0), fCutRequireITSRefit(0), fCutAcceptKinkDaughters(0),
96957075 75 fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0),fCutDCAToVertex2D(0),fPIDUtils(),
76 fUseTimeCorrectionFactors(kFALSE), fTimeCorrectionFactorsSet(kFALSE)
d9b3567c 77{
78//
79 // Constructor.
80 // Initialize all constant values which have to be used
81 // during Reco algorithm execution
82 //
83
b540d03f 84 //Misalignment matrices
fd6df01c 85 for(Int_t i = 0; i < 15 ; i++) {
86 fMisalTransShift[i] = 0.;
b540d03f 87 fMisalRotShift[i] = 0.;
fd6df01c 88 }
b540d03f 89
90 //Non linearity
dff9e2e3 91 for(Int_t i = 0; i < 7 ; i++) fNonLinearityParams[i] = 0.;
92
93 //For kBeamTestCorrected case, but default is no correction
94 fNonLinearityParams[0] = 0.99078;
95 fNonLinearityParams[1] = 0.161499;
96 fNonLinearityParams[2] = 0.655166;
97 fNonLinearityParams[3] = 0.134101;
98 fNonLinearityParams[4] = 163.282;
99 fNonLinearityParams[5] = 23.6904;
100 fNonLinearityParams[6] = 0.978;
101
102 //For kPi0GammaGamma case
103 //fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
104 //fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
105 //fNonLinearityParams[2] = 1.046;
bd8c7aef 106
b540d03f 107 //Track matching
108 fMatchedTrackIndex = new TArrayI();
bd8c7aef 109 fMatchedClusterIndex = new TArrayI();
b540d03f 110 fResidualZ = new TArrayF();
111 fResidualR = new TArrayF();
d9b3567c 112
bd8c7aef 113 InitTrackCuts();
b540d03f 114
115 fPIDUtils = new AliEMCALPIDUtils();
116
bd8c7aef 117
d9b3567c 118}
119
120//______________________________________________________________________
121AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
094786cc 122: TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction),
fd6df01c 123 fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0),
124 fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
78467229 125 fRemoveBadChannels(reco.fRemoveBadChannels),fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
126 fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
83bfd77a 127 fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
bb6f5f0b 128 fAODFilterMask(reco.fAODFilterMask),
b540d03f 129 fMatchedTrackIndex(reco.fMatchedTrackIndex?new TArrayI(*reco.fMatchedTrackIndex):0x0),
bd8c7aef 130 fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
131 fResidualZ(reco.fResidualZ?new TArrayF(*reco.fResidualZ):0x0),
132 fResidualR(reco.fResidualR?new TArrayF(*reco.fResidualR):0x0),
bb6f5f0b 133 fCutR(reco.fCutR),fCutZ(reco.fCutZ),fMass(reco.fMass), fStep(reco.fStep),
b5078f5d 134 fRejectExoticCluster(reco.fRejectExoticCluster),
bd8c7aef 135 fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS),
136 fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
137 fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
138 fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters),
139 fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY), fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ),fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
96957075 140 fPIDUtils(reco.fPIDUtils),
141 fUseTimeCorrectionFactors(reco.fUseTimeCorrectionFactors), fTimeCorrectionFactorsSet(reco.fTimeCorrectionFactorsSet)
d9b3567c 142{
143 //Copy ctor
144
fd6df01c 145 for(Int_t i = 0; i < 15 ; i++) {
146 fMisalRotShift[i] = reco.fMisalRotShift[i];
147 fMisalTransShift[i] = reco.fMisalTransShift[i];
148 }
dff9e2e3 149 for(Int_t i = 0; i < 7 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
bd8c7aef 150
d9b3567c 151}
152
153
154//______________________________________________________________________
155AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
156{
157 //Assignment operator
158
159 if(this == &reco)return *this;
160 ((TNamed *)this)->operator=(reco);
161
96957075 162 fNonLinearityFunction = reco.fNonLinearityFunction;
163 fParticleType = reco.fParticleType;
164 fPosAlgo = reco.fPosAlgo;
165 fW0 = reco.fW0;
166 fRecalibration = reco.fRecalibration;
094786cc 167 fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
96957075 168 fRemoveBadChannels = reco.fRemoveBadChannels;
169 fRecalDistToBadChannels = reco.fRecalDistToBadChannels;
170 fEMCALBadChannelMap = reco.fEMCALBadChannelMap;
171 fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
172 fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
bd8c7aef 173
83bfd77a 174
2a71e873 175 for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
dff9e2e3 176 for(Int_t i = 0; i < 7 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
bb6f5f0b 177
178 fAODFilterMask = reco.fAODFilterMask;
d9b3567c 179
96957075 180 fCutR = reco.fCutR;
181 fCutZ = reco.fCutZ;
bb6f5f0b 182 fMass = reco.fMass;
183 fStep = reco.fStep;
b5078f5d 184 fRejectExoticCluster = reco.fRejectExoticCluster;
bd8c7aef 185
96957075 186 fCutMinNClusterTPC = reco.fCutMinNClusterTPC;
187 fCutMinNClusterITS = reco.fCutMinNClusterITS;
188 fCutMaxChi2PerClusterTPC = reco.fCutMaxChi2PerClusterTPC;
189 fCutMaxChi2PerClusterITS = reco.fCutMaxChi2PerClusterITS;
190 fCutRequireTPCRefit = reco.fCutRequireTPCRefit;
191 fCutRequireITSRefit = reco.fCutRequireITSRefit;
192 fCutAcceptKinkDaughters = reco.fCutAcceptKinkDaughters;
193 fCutMaxDCAToVertexXY = reco.fCutMaxDCAToVertexXY;
194 fCutMaxDCAToVertexZ = reco.fCutMaxDCAToVertexZ;
195 fCutDCAToVertex2D = reco.fCutDCAToVertex2D;
bd8c7aef 196
96957075 197 fPIDUtils = reco.fPIDUtils;
bd8c7aef 198
96957075 199 fUseTimeCorrectionFactors = reco.fUseTimeCorrectionFactors;
200 fTimeCorrectionFactorsSet = reco.fTimeCorrectionFactorsSet;
201
bd8c7aef 202
203 if(reco.fResidualR){
204 // assign or copy construct
205 if(fResidualR){
206 *fResidualR = *reco.fResidualR;
207 }
208 else fResidualR = new TArrayF(*reco.fResidualR);
209 }
210 else{
211 if(fResidualR)delete fResidualR;
212 fResidualR = 0;
213 }
214
215 if(reco.fResidualZ){
216 // assign or copy construct
217 if(fResidualZ){
218 *fResidualZ = *reco.fResidualZ;
219 }
220 else fResidualZ = new TArrayF(*reco.fResidualZ);
221 }
222 else{
223 if(fResidualZ)delete fResidualZ;
224 fResidualZ = 0;
225 }
226
b540d03f 227 if(reco.fMatchedTrackIndex){
228 // assign or copy construct
229 if(fMatchedTrackIndex){
230 *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
231 }
232 else fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
233 }
234 else{
235 if(fMatchedTrackIndex)delete fMatchedTrackIndex;
236 fMatchedTrackIndex = 0;
237 }
bd8c7aef 238
239 if(reco.fMatchedClusterIndex){
240 // assign or copy construct
241 if(fMatchedClusterIndex){
242 *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
243 }
244 else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
245 }
246 else{
247 if(fMatchedClusterIndex)delete fMatchedClusterIndex;
248 fMatchedClusterIndex = 0;
249 }
250
251
d9b3567c 252 return *this;
253}
254
255
094786cc 256//__________________________________________________
257AliEMCALRecoUtils::~AliEMCALRecoUtils()
258{
259 //Destructor.
260
261 if(fEMCALRecalibrationFactors) {
262 fEMCALRecalibrationFactors->Clear();
263 delete fEMCALRecalibrationFactors;
264 }
fd6df01c 265
266 if(fEMCALBadChannelMap) {
267 fEMCALBadChannelMap->Clear();
268 delete fEMCALBadChannelMap;
269 }
bd8c7aef 270
b540d03f 271 if(fMatchedTrackIndex) {delete fMatchedTrackIndex; fMatchedTrackIndex=0;}
bd8c7aef 272 if(fMatchedClusterIndex) {delete fMatchedClusterIndex; fMatchedClusterIndex=0;}
273 if(fResidualR) {delete fResidualR; fResidualR=0;}
274 if(fResidualZ) {delete fResidualZ; fResidualZ=0;}
275
094786cc 276}
277
fd6df01c 278//_______________________________________________________________
279Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells)
280{
281 // Given the list of AbsId of the cluster, get the maximum cell and
282 // check if there are fNCellsFromBorder from the calorimeter border
283
284 //If the distance to the border is 0 or negative just exit accept all clusters
285 if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
286
cb231979 287 Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
288 Bool_t shared = kFALSE;
289 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
fd6df01c 290
83bfd77a 291 AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
292 absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
fd6df01c 293
294 if(absIdMax==-1) return kFALSE;
295
296 //Check if the cell is close to the borders:
297 Bool_t okrow = kFALSE;
298 Bool_t okcol = kFALSE;
299
300 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
301 AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
302 iSM,ieta,iphi));
303 }
304
305 //Check rows/phi
306 if(iSM < 10){
307 if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
308 }
309 else{
310 if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
311 }
312
313 //Check columns/eta
314 if(!fNoEMCALBorderAtEta0){
315 if(ieta > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
316 }
317 else{
318 if(iSM%2==0){
319 if(ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
320 }
321 else {
322 if(ieta < 48-fNCellsFromEMCALBorder) okcol = kTRUE;
323 }
324 }//eta 0 not checked
325
83bfd77a 326 AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
fd6df01c 327 fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
328
83bfd77a 329 if (okcol && okrow) {
330 //printf("Accept\n");
331 return kTRUE;
332 }
333 else {
334 //printf("Reject\n");
335 AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
336 return kFALSE;
337 }
fd6df01c 338
339}
340
341
342//_________________________________________________________________________________________________________
343Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells){
344 // Check that in the cluster cells, there is no bad channel of those stored
345 // in fEMCALBadChannelMap or fPHOSBadChannelMap
346
347 if(!fRemoveBadChannels) return kFALSE;
348 if(!fEMCALBadChannelMap) return kFALSE;
349
350 Int_t icol = -1;
351 Int_t irow = -1;
352 Int_t imod = -1;
353 for(Int_t iCell = 0; iCell<nCells; iCell++){
354
355 //Get the column and row
356 Int_t iTower = -1, iIphi = -1, iIeta = -1;
357 geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
358 if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
359 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
83bfd77a 360 if(GetEMCALChannelStatus(imod, icol, irow)){
361 AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
362 return kTRUE;
363 }
fd6df01c 364
365 }// cell cluster loop
366
367 return kFALSE;
368
369}
094786cc 370
b5078f5d 371//_________________________________________________
372Bool_t AliEMCALRecoUtils::IsExoticCluster(AliVCluster *cluster){
373 // Check if the cluster has high energy but small number of cells
374 // The criteria comes from Gustavo's study
375 //
376
377 if(cluster->GetNCells()<(1+cluster->E()/3.))
378 return kTRUE;
379 else
380 return kFALSE;
381}
382
d9b3567c 383//__________________________________________________
384Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
385// Correct cluster energy from non linearity functions
386 Float_t energy = cluster->E();
387
388 switch (fNonLinearityFunction) {
389
390 case kPi0MC:
871aee7a 391 {
d9b3567c 392 //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
dff9e2e3 393 //Double_t fNonLinearityParams[0] = 1.014;
394 //Double_t fNonLinearityParams[1] = -0.03329;
395 //Double_t fNonLinearityParams[2] = -0.3853;
396 //Double_t fNonLinearityParams[3] = 0.5423;
397 //Double_t fNonLinearityParams[4] = -0.4335;
8cdd1f1f 398 energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
d9b3567c 399 ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
400 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
401 break;
871aee7a 402 }
dff9e2e3 403
d9b3567c 404 case kPi0GammaGamma:
871aee7a 405 {
d9b3567c 406 //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
96957075 407 //Double_t fNonLinearityParams[0] = 1.04;
408 //Double_t fNonLinearityParams[1] = -0.1445;
871aee7a 409 //Double_t fNonLinearityParams[2] = 1.046;
d9b3567c 410 energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
411 break;
871aee7a 412 }
d9b3567c 413
414 case kPi0GammaConversion:
871aee7a 415 {
d9b3567c 416 //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
871aee7a 417 //fNonLinearityParams[0] = 0.139393/0.1349766;
418 //fNonLinearityParams[1] = 0.0566186;
419 //fNonLinearityParams[2] = 0.982133;
d9b3567c 420 energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
421
422 break;
871aee7a 423 }
424
425 case kBeamTest:
426 {
427 //From beam test, Alexei's results, for different ZS thresholds
428 // th=30 MeV; th = 45 MeV; th = 75 MeV
96957075 429 //fNonLinearityParams[0] = 1.007; 1.003; 1.002
871aee7a 430 //fNonLinearityParams[1] = 0.894; 0.719; 0.797
431 //fNonLinearityParams[2] = 0.246; 0.334; 0.358
96957075 432 //Rescale the param[0] with 1.03
871aee7a 433 energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
434
435 break;
436 }
dff9e2e3 437
4b58ac4f 438 case kBeamTestCorrected:
439 {
440 //From beam test, corrected for material between beam and EMCAL
dff9e2e3 441 //fNonLinearityParams[0] = 0.99078
442 //fNonLinearityParams[1] = 0.161499;
443 //fNonLinearityParams[2] = 0.655166;
444 //fNonLinearityParams[3] = 0.134101;
445 //fNonLinearityParams[4] = 163.282;
446 //fNonLinearityParams[5] = 23.6904;
447 //fNonLinearityParams[6] = 0.978;
448 energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
449
4b58ac4f 450 break;
451 }
d9b3567c 452
453 case kNoCorrection:
454 AliDebug(2,"No correction on the energy\n");
455 break;
456
457 }
458
459 return energy;
460
461}
d9b3567c 462//__________________________________________________
094786cc 463Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
464{
465 //Calculate shower depth for a given cluster energy and particle type
466
467 // parameters
cb231979 468 Float_t x0 = 1.31;
094786cc 469 Float_t ecr = 8;
470 Float_t depth = 0;
471
472 switch ( iParticle )
473 {
474 case kPhoton:
fd6df01c 475 depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 476 break;
477
478 case kElectron:
fd6df01c 479 depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 480 break;
481
482 case kHadron:
483 // hadron
484 // boxes anc. here
485 if(gGeoManager){
486 gGeoManager->cd("ALIC_1/XEN1_1");
487 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
488 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
fd6df01c 489 if(geoSM){
490 TGeoVolume *geoSMVol = geoSM->GetVolume();
491 TGeoShape *geoSMShape = geoSMVol->GetShape();
492 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
493 if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
494 else AliFatal("Null GEANT box");
495 }else AliFatal("NULL GEANT node matrix");
094786cc 496 }
497 else{//electron
fd6df01c 498 depth = x0 * (TMath::Log(energy*1000 / ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 499 }
500
501 break;
502
503 default://photon
fd6df01c 504 depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 505 }
506
507 return depth;
508
509}
510
511//__________________________________________________
cb231979 512void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
513 Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
d9b3567c 514{
515 //For a given CaloCluster gets the absId of the cell
516 //with maximum energy deposit.
517
518 Double_t eMax = -1.;
519 Double_t eCell = -1.;
094786cc 520 Float_t fraction = 1.;
521 Float_t recalFactor = 1.;
d9b3567c 522 Int_t cellAbsId = -1 ;
094786cc 523
d9b3567c 524 Int_t iTower = -1;
525 Int_t iIphi = -1;
526 Int_t iIeta = -1;
cb231979 527 Int_t iSupMod0= -1;
83bfd77a 528 //printf("---Max?\n");
d9b3567c 529 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 530 cellAbsId = clu->GetCellAbsId(iDig);
531 fraction = clu->GetCellAmplitudeFraction(iDig);
83bfd77a 532 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
094786cc 533 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
cb231979 534 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
535 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
536 if(iDig==0) iSupMod0=iSupMod;
537 else if(iSupMod0!=iSupMod) {
538 shared = kTRUE;
539 //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
540 }
094786cc 541 if(IsRecalibrationOn()) {
094786cc 542 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
543 }
544 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
83bfd77a 545 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
094786cc 546 if(eCell > eMax) {
d9b3567c 547 eMax = eCell;
548 absId = cellAbsId;
549 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
550 }
551 }// cell loop
552
553 //Get from the absid the supermodule, tower and eta/phi numbers
554 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
555 //Gives SuperModule and Tower numbers
556 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
83bfd77a 557 iIphi, iIeta,iphi,ieta);
558 //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
559 //printf("Max end---\n");
d9b3567c 560
561}
562
094786cc 563//________________________________________________________________
564void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
565 //Init EMCAL recalibration factors
566 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
567 //In order to avoid rewriting the same histograms
568 Bool_t oldStatus = TH1::AddDirectoryStatus();
569 TH1::AddDirectory(kFALSE);
570
cb231979 571 fEMCALRecalibrationFactors = new TObjArray(10);
d7194a80 572 for (int i = 0; i < 10; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
094786cc 573 //Init the histograms with 1
d7194a80 574 for (Int_t sm = 0; sm < 10; sm++) {
094786cc 575 for (Int_t i = 0; i < 48; i++) {
576 for (Int_t j = 0; j < 24; j++) {
577 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
578 }
579 }
580 }
581 fEMCALRecalibrationFactors->SetOwner(kTRUE);
582 fEMCALRecalibrationFactors->Compress();
583
584 //In order to avoid rewriting the same histograms
585 TH1::AddDirectory(oldStatus);
586}
587
588
fd6df01c 589//________________________________________________________________
590void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
591 //Init EMCAL bad channels map
592 AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
593 //In order to avoid rewriting the same histograms
594 Bool_t oldStatus = TH1::AddDirectoryStatus();
595 TH1::AddDirectory(kFALSE);
596
cb231979 597 fEMCALBadChannelMap = new TObjArray(10);
fd6df01c 598 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
6fe0e6d0 599 for (int i = 0; i < 10; i++) {
fd6df01c 600 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
601 }
602
603 //delete hTemp;
604
605 fEMCALBadChannelMap->SetOwner(kTRUE);
606 fEMCALBadChannelMap->Compress();
607
608 //In order to avoid rewriting the same histograms
609 TH1::AddDirectory(oldStatus);
610}
611
094786cc 612//________________________________________________________________
613void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
614 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
615
616 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
617 UShort_t * index = cluster->GetCellsAbsId() ;
618 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
619 Int_t ncells = cluster->GetNCells();
620
621 //Initialize some used variables
622 Float_t energy = 0;
623 Int_t absId = -1;
624 Int_t icol = -1, irow = -1, imod=1;
625 Float_t factor = 1, frac = 0;
626
627 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
628 for(Int_t icell = 0; icell < ncells; icell++){
629 absId = index[icell];
630 frac = fraction[icell];
631 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
632 Int_t iTower = -1, iIphi = -1, iIeta = -1;
633 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
634 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
635 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
636 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
637 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
638 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
639
640 energy += cells->GetCellAmplitude(absId)*factor*frac;
641 }
642
643
644 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
645
646 cluster->SetE(energy);
647
648}
649
650
d9b3567c 651//__________________________________________________
094786cc 652void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
d9b3567c 653{
654 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
655
094786cc 656 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
657 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
fd6df01c 658 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
094786cc 659
660}
661
662//__________________________________________________
663void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
664{
665 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
666 // The algorithm is a copy of what is done in AliEMCALRecPoint
667
668 Double_t eCell = 0.;
669 Float_t fraction = 1.;
670 Float_t recalFactor = 1.;
671
672 Int_t absId = -1;
673 Int_t iTower = -1, iIphi = -1, iIeta = -1;
674 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
675 Float_t weight = 0., totalWeight=0.;
676 Float_t newPos[3] = {0,0,0};
677 Double_t pLocal[3], pGlobal[3];
cb231979 678 Bool_t shared = kFALSE;
679
094786cc 680 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
cb231979 681 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 682 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
683
83bfd77a 684 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
685
094786cc 686 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
687 absId = clu->GetCellAbsId(iDig);
688 fraction = clu->GetCellAmplitudeFraction(iDig);
689 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
690 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
691 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
692
693 if(IsRecalibrationOn()) {
694 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
695 }
696 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
697
698 weight = GetCellWeight(eCell,clEnergy);
83bfd77a 699 //printf("cell energy %f, weight %f\n",eCell,weight);
094786cc 700 totalWeight += weight;
701 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
83bfd77a 702 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
094786cc 703 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
83bfd77a 704 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
705
094786cc 706 for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
707
708 }// cell loop
709
710 if(totalWeight>0){
711 for(int i=0; i<3; i++ ) newPos[i] /= totalWeight;
712 }
713
094786cc 714 //Float_t pos[]={0,0,0};
715 //clu->GetPosition(pos);
716 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
83bfd77a 717 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
094786cc 718
719 if(iSupModMax > 1) {//sector 1
720 newPos[0] +=fMisalTransShift[3];//-=3.093;
721 newPos[1] +=fMisalTransShift[4];//+=6.82;
722 newPos[2] +=fMisalTransShift[5];//+=1.635;
83bfd77a 723 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
724
094786cc 725 }
726 else {//sector 0
727 newPos[0] +=fMisalTransShift[0];//+=1.134;
728 newPos[1] +=fMisalTransShift[1];//+=8.2;
729 newPos[2] +=fMisalTransShift[2];//+=1.197;
83bfd77a 730 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
731
094786cc 732 }
83bfd77a 733 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
734
094786cc 735 clu->SetPosition(newPos);
736
094786cc 737}
738
739//__________________________________________________
740void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
741{
742 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
743 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
744
745 Double_t eCell = 1.;
746 Float_t fraction = 1.;
747 Float_t recalFactor = 1.;
748
749 Int_t absId = -1;
d9b3567c 750 Int_t iTower = -1;
094786cc 751 Int_t iIphi = -1, iIeta = -1;
752 Int_t iSupMod = -1, iSupModMax = -1;
d9b3567c 753 Int_t iphi = -1, ieta =-1;
cb231979 754 Bool_t shared = kFALSE;
755
d9b3567c 756 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
cb231979 757 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 758 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
759
d9b3567c 760 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
094786cc 761 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
762 Int_t startingSM = -1;
d9b3567c 763
764 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 765 absId = clu->GetCellAbsId(iDig);
766 fraction = clu->GetCellAmplitudeFraction(iDig);
767 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
768 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
d9b3567c 769 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
094786cc 770
d9b3567c 771 if (iDig==0) startingSM = iSupMod;
772 else if(iSupMod != startingSM) areInSameSM = kFALSE;
094786cc 773
774 eCell = cells->GetCellAmplitude(absId);
d9b3567c 775
094786cc 776 if(IsRecalibrationOn()) {
777 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
778 }
779 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
d9b3567c 780
094786cc 781 weight = GetCellWeight(eCell,clEnergy);
d9b3567c 782 if(weight < 0) weight = 0;
783 totalWeight += weight;
784 weightedCol += ieta*weight;
785 weightedRow += iphi*weight;
786
787 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
788
094786cc 789 }// cell loop
790
d9b3567c 791 Float_t xyzNew[]={0.,0.,0.};
792 if(areInSameSM == kTRUE) {
793 //printf("In Same SM\n");
794 weightedCol = weightedCol/totalWeight;
795 weightedRow = weightedRow/totalWeight;
094786cc 796 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 797 }
798 else {
799 //printf("In Different SM\n");
094786cc 800 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 801 }
d9b3567c 802
094786cc 803 clu->SetPosition(xyzNew);
d9b3567c 804
805}
806
cb231979 807//____________________________________________________________________________
808void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
809
810 //re-evaluate distance to bad channel with updated bad map
811
78467229 812 if(!fRecalDistToBadChannels) return;
cb231979 813
814 //Get channels map of the supermodule where the cluster is.
cb231979 815 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
816 Bool_t shared = kFALSE;
817 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
818 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
819
820 Int_t dRrow, dRcol;
821 Float_t minDist = 10000.;
822 Float_t dist = 0.;
823
824 //Loop on tower status map
825 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
826 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
827 //Check if tower is bad.
828 if(hMap->GetBinContent(icol,irow)==0) continue;
829 //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 830 // iSupMod,icol, irow, icolM,irowM);
cb231979 831
832 dRrow=TMath::Abs(irowM-irow);
833 dRcol=TMath::Abs(icolM-icol);
834 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
835 if(dist < minDist){
836 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
837 minDist = dist;
838 }
839
840 }
841 }
842
843 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
844 if (shared) {
845 TH2D* hMap2 = 0;
846 Int_t iSupMod2 = -1;
847
848 //The only possible combinations are (0,1), (2,3) ... (8,9)
849 if(iSupMod%2) iSupMod2 = iSupMod-1;
850 else iSupMod2 = iSupMod+1;
851 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
852
853 //Loop on tower status map of second super module
854 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
855 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
856 //Check if tower is bad.
857 if(hMap2->GetBinContent(icol,irow)==0) continue;
858 //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",
859 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
860
861 dRrow=TMath::Abs(irow-irowM);
862
863 if(iSupMod%2) {
864 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
865 }
866 else {
867 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
868 }
869
870 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
871 if(dist < minDist) minDist = dist;
872
873 }
874 }
875
876 }// shared cluster in 2 SuperModules
78467229 877
6fe0e6d0 878 AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
879 cluster->SetDistanceToBadChannel(minDist);
cb231979 880
881}
882
83bfd77a 883//____________________________________________________________________________
884void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
885
886 //re-evaluate identification parameters with bayesian
887
888 if ( cluster->GetM02() != 0)
889 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
890
891 Float_t pidlist[AliPID::kSPECIESN+1];
892 for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
893
894 cluster->SetPID(pidlist);
895
896}
897
898//____________________________________________________________________________
899void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
900{
901 // Calculates new center of gravity in the local EMCAL-module coordinates
902 // and tranfers into global ALICE coordinates
903 // Calculates Dispersion and main axis
904
905 Int_t nstat = 0;
906 Float_t wtot = 0. ;
907 Double_t eCell = 0.;
908 Float_t fraction = 1.;
909 Float_t recalFactor = 1.;
910
911 Int_t iSupMod = -1;
912 Int_t iTower = -1;
913 Int_t iIphi = -1;
914 Int_t iIeta = -1;
915 Int_t iphi = -1;
916 Int_t ieta = -1;
917 Double_t etai = -1.;
918 Double_t phii = -1.;
919
920 Double_t w = 0.;
921 Double_t d = 0.;
922 Double_t dxx = 0.;
923 Double_t dzz = 0.;
924 Double_t dxz = 0.;
925 Double_t xmean = 0.;
926 Double_t zmean = 0.;
927
928 //Loop on cells
929 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
930
931 //Get from the absid the supermodule, tower and eta/phi numbers
932 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
933 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
934
935 //Get the cell energy, if recalibration is on, apply factors
936 fraction = cluster->GetCellAmplitudeFraction(iDigit);
937 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
938 if(IsRecalibrationOn()) {
939 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
940 }
941 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
942
943 if(cluster->E() > 0 && eCell > 0){
944
945 w = GetCellWeight(eCell,cluster->E());
946
947 etai=(Double_t)ieta;
948 phii=(Double_t)iphi;
949 if(w > 0.0) {
950 wtot += w ;
951 nstat++;
952 //Shower shape
953 dxx += w * etai * etai ;
954 xmean+= w * etai ;
955 dzz += w * phii * phii ;
956 zmean+= w * phii ;
957 dxz += w * etai * phii ;
958 }
959 }
960 else
961 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
962 }//cell loop
963
964 //Normalize to the weight
965 if (wtot > 0) {
966 xmean /= wtot ;
967 zmean /= wtot ;
968 }
969 else
970 AliError(Form("Wrong weight %f\n", wtot));
971
972 //Calculate dispersion
973 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
974
975 //Get from the absid the supermodule, tower and eta/phi numbers
976 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
977 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
978
979 //Get the cell energy, if recalibration is on, apply factors
980 fraction = cluster->GetCellAmplitudeFraction(iDigit);
981 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
982 if(IsRecalibrationOn()) {
983 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
984 }
985 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
986
987 if(cluster->E() > 0 && eCell > 0){
988
989 w = GetCellWeight(eCell,cluster->E());
990
991 etai=(Double_t)ieta;
992 phii=(Double_t)iphi;
993 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
994 }
995 else
996 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
997 }// cell loop
998
999 //Normalize to the weigth and set shower shape parameters
1000 if (wtot > 0 && nstat > 1) {
1001 d /= wtot ;
1002 dxx /= wtot ;
1003 dzz /= wtot ;
1004 dxz /= wtot ;
1005 dxx -= xmean * xmean ;
1006 dzz -= zmean * zmean ;
1007 dxz -= xmean * zmean ;
1008 cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1009 cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1010 }
1011 else{
1012 d=0. ;
1013 cluster->SetM20(0.) ;
1014 cluster->SetM02(0.) ;
1015 }
1016
1017 if (d>=0)
1018 cluster->SetDispersion(TMath::Sqrt(d)) ;
1019 else
1020 cluster->SetDispersion(0) ;
83bfd77a 1021}
1022
b540d03f 1023//____________________________________________________________________________
b5078f5d 1024void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr, AliEMCALGeometry *geom)
bd8c7aef 1025{
1026 //This function should be called before the cluster loop
1027 //Before call this function, please recalculate the cluster positions
1028 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
1029 //Store matched cluster indexes and residuals
bb6f5f0b 1030
b540d03f 1031 fMatchedTrackIndex ->Reset();
bd8c7aef 1032 fMatchedClusterIndex->Reset();
b540d03f 1033 fResidualZ ->Reset();
1034 fResidualR ->Reset();
bd8c7aef 1035
67241b7e 1036 fMatchedTrackIndex ->Set(500);
1037 fMatchedClusterIndex->Set(500);
1038 fResidualZ ->Set(500);
1039 fResidualR ->Set(500);
bd8c7aef 1040
1c7a2bf4 1041 AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1042 AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
1043
bd8c7aef 1044 Int_t matched=0;
bb6f5f0b 1045 Double_t cv[21];
1046 for (Int_t i=0; i<21;i++) cv[i]=0;
bd8c7aef 1047 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1048 {
bb6f5f0b 1049 AliExternalTrackParam *trackParam=0;
1050
1051 //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
1c7a2bf4 1052 if(esdevent)
bb6f5f0b 1053 {
1c7a2bf4 1054 AliESDtrack *esdTrack = esdevent->GetTrack(itr);
bb6f5f0b 1055 if(!esdTrack || !IsAccepted(esdTrack)) continue;
1056 const AliESDfriendTrack* friendTrack = esdTrack->GetFriendTrack();
1057 if(friendTrack && friendTrack->GetTPCOut())
1058 {
1059 //Use TPC Out as starting point if it is available
1060 trackParam= new AliExternalTrackParam(*friendTrack->GetTPCOut());
1061 }
1062 else
1063 {
1064 //Otherwise use TPC inner
1065 trackParam = new AliExternalTrackParam(*esdTrack->GetInnerParam());
1066 }
1067 }
1068
1069 //If the input event is AOD, the starting point for extrapolation is at vertex
1070 //AOD tracks are selected according to its bit.
1c7a2bf4 1071 else if(aodevent)
bb6f5f0b 1072 {
1c7a2bf4 1073 AliAODTrack *aodTrack = aodevent->GetTrack(itr);
bb6f5f0b 1074 if(!aodTrack) continue;
1075 if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
1076 Double_t pos[3],mom[3];
1077 aodTrack->GetXYZ(pos);
1078 aodTrack->GetPxPyPz(mom);
1079 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()));
1080 trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1081 }
bd8c7aef 1082
bb6f5f0b 1083 //Return if the input data is not "AOD" or "ESD"
1084 else
1085 {
b5078f5d 1086 printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
bb6f5f0b 1087 return;
1088 }
1089
1090 if(!trackParam) continue;
1091
1092 Float_t dRMax = fCutR, dZMax=fCutZ;
bd8c7aef 1093 Int_t index = -1;
b540d03f 1094 if(!clusterArr){// get clusters from event
1095 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1096 {
bb6f5f0b 1097 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
b540d03f 1098 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
b5078f5d 1099 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
bb6f5f0b 1100 Float_t tmpR=-1, tmpZ=-1;
1101 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
b540d03f 1102 if(tmpR<dRMax)
1103 {
1104 dRMax=tmpR;
1105 dZMax=tmpZ;
1106 index=icl;
1107 }
bb6f5f0b 1108 delete trkPamTmp;
b540d03f 1109 }//cluster loop
1110 } else { // external cluster array, not from ESD event
1111 for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1112 {
bb6f5f0b 1113 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
b540d03f 1114 AliVCluster *cluster = (AliVCluster*) clusterArr->At(icl);
1115 if(!cluster->IsEMCAL()) continue;
bb6f5f0b 1116 Float_t tmpR=-1, tmpZ=-1;
1117 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
b540d03f 1118 if(tmpR<dRMax)
1119 {
1120 dRMax=tmpR;
1121 dZMax=tmpZ;
1122 index=icl;
1123 }
bb6f5f0b 1124 delete trkPamTmp;
b540d03f 1125 }//cluster loop
1126 }// external list of clusters
1127
bd8c7aef 1128 if(index>-1)
1129 {
b540d03f 1130 fMatchedTrackIndex ->AddAt(itr,matched);
bd8c7aef 1131 fMatchedClusterIndex->AddAt(index,matched);
b540d03f 1132 fResidualZ ->AddAt(dZMax,matched);
1133 fResidualR ->AddAt(dRMax,matched);
bd8c7aef 1134 matched++;
1135 }
bb6f5f0b 1136 delete trackParam;
bd8c7aef 1137 }//track loop
b540d03f 1138
1139 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1140
1141 fMatchedTrackIndex ->Set(matched);
bd8c7aef 1142 fMatchedClusterIndex->Set(matched);
b540d03f 1143 fResidualZ ->Set(matched);
1144 fResidualR ->Set(matched);
bd8c7aef 1145}
1146
b540d03f 1147//________________________________________________________________________________
b5078f5d 1148Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom)
bb6f5f0b 1149{
1150 //
1151 // This function returns the index of matched cluster to input track
1152 // Cut on match is dR<10cm by default. Returns -1 if no match is found
1153
1154
1155 Float_t dRMax = fCutR;
1156 Int_t index = -1;
1157
1158 AliExternalTrackParam *trackParam=0;
1159 const AliESDfriendTrack* friendTrack = track->GetFriendTrack();
1160 if(friendTrack && friendTrack->GetTPCOut())
1161 trackParam= const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
1162 else
1163 trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
1164
1165 if(!trackParam) return index;
1166 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1167 {
1168 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
1169 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
b5078f5d 1170 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
bb6f5f0b 1171 Float_t tmpR=-1, tmpZ=-1;
1172 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
1173 if(tmpR>-1 && tmpR<dRMax)
1174 {
1175 dRMax=tmpR;
1176 index=icl;
1177 }
1178 delete trkPamTmp;
1179 }//cluster loop
1180 return index;
1181}
1182
1183//________________________________________________________________________________
1184Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpR, Float_t &tmpZ)
1185{
1186 //
1187 //Return the residual by extrapolating a track to a cluster
1188 //
1189 if(!trkParam || !cluster) return kFALSE;
1190 Float_t clsPos[3];
1191 Double_t trkPos[3];
1192 cluster->GetPosition(clsPos); //Has been recalculated
1193 TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
1194 Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
1195 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
1196 trkParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
1197 if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kFALSE)) return kFALSE;
1198 trkParam->GetXYZ(trkPos); //Get the extrapolated global position
1199 tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
1200 tmpZ = clsPos[2]-trkPos[2];
1201 return kTRUE;
1202}
1203
1204//________________________________________________________________________________
1205void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dR, Float_t &dZ)
bd8c7aef 1206{
bb6f5f0b 1207 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1208 //Get the residuals dR and dZ for this cluster to the closest track
1209 //Works with ESDs and AODs
bd8c7aef 1210
bb6f5f0b 1211 if( FindMatchedPosForCluster(clsIndex) >= 999 )
bd8c7aef 1212 {
1213 AliDebug(2,"No matched tracks found!\n");
1214 dR=999.;
1215 dZ=999.;
1216 return;
1217 }
bb6f5f0b 1218 dR = fResidualR->At(FindMatchedPosForCluster(clsIndex));
1219 dZ = fResidualZ->At(FindMatchedPosForCluster(clsIndex));
1220}
1221//________________________________________________________________________________
1222void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dR, Float_t &dZ)
1223{
1224 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1225 //Get the residuals dR and dZ for this track to the closest cluster
1226 //Works with ESDs and AODs
1227
1228 if( FindMatchedPosForTrack(trkIndex) >= 999 )
1229 {
1230 AliDebug(2,"No matched cluster found!\n");
1231 dR=999.;
1232 dZ=999.;
1233 return;
1234 }
1235 dR = fResidualR->At(FindMatchedPosForTrack(trkIndex));
1236 dZ = fResidualZ->At(FindMatchedPosForTrack(trkIndex));
1237}
1238
1239//__________________________________________________________
1240Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
1241{
1242 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1243 //Get the index of matched track to this cluster
1244 //Works with ESDs and AODs
1245
1246 if(IsClusterMatched(clsIndex))
1247 return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
1248 else
1249 return -1;
bd8c7aef 1250}
1251
b540d03f 1252//__________________________________________________________
bb6f5f0b 1253Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
b540d03f 1254{
bb6f5f0b 1255 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1256 //Get the index of matched cluster to this track
1257 //Works with ESDs and AODs
b540d03f 1258
bb6f5f0b 1259 if(IsTrackMatched(trkIndex))
1260 return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
b540d03f 1261 else
1262 return -1;
1263}
1264
bb6f5f0b 1265//__________________________________________________
1266Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex)
1267{
1268 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1269 //Returns if the cluster has a match
1270 if(FindMatchedPosForCluster(clsIndex) < 999)
1271 return kTRUE;
1272 else
1273 return kFALSE;
1274}
b540d03f 1275
bd8c7aef 1276//__________________________________________________
bb6f5f0b 1277Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex)
bd8c7aef 1278{
bb6f5f0b 1279 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1280 //Returns if the track has a match
1281 if(FindMatchedPosForTrack(trkIndex) < 999)
82d09e74 1282 return kTRUE;
bd8c7aef 1283 else
1284 return kFALSE;
1285}
bb6f5f0b 1286
b540d03f 1287//__________________________________________________________
bb6f5f0b 1288UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
bd8c7aef 1289{
bb6f5f0b 1290 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
bd8c7aef 1291 //Returns the position of the match in the fMatchedClusterIndex array
1292 Float_t tmpR = fCutR;
81efb149 1293 UInt_t pos = 999;
b540d03f 1294
bd8c7aef 1295 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
b540d03f 1296 {
bb6f5f0b 1297 if(fMatchedClusterIndex->At(i)==clsIndex && fResidualR->At(i)<tmpR)
1298 {
1299 pos=i;
1300 tmpR=fResidualR->At(i);
1301 AliDebug(3,Form("Matched cluster index: index: %d, dR: %2.4f, dZ: %2.4f.\n",fMatchedClusterIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
1302 }
1303 }
1304 return pos;
1305}
1306
1307//__________________________________________________________
1308UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
1309{
1310 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1311 //Returns the position of the match in the fMatchedTrackIndex array
1312 Float_t tmpR = fCutR;
1313 UInt_t pos = 999;
1314
1315 for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
1316 {
1317 if(fMatchedTrackIndex->At(i)==trkIndex && fResidualR->At(i)<tmpR)
bd8c7aef 1318 {
b540d03f 1319 pos=i;
1320 tmpR=fResidualR->At(i);
bb6f5f0b 1321 AliDebug(3,Form("Matched track index: index: %d, dR: %2.4f, dZ: %2.4f.\n",fMatchedTrackIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
bd8c7aef 1322 }
b540d03f 1323 }
bd8c7aef 1324 return pos;
1325}
1326
b5078f5d 1327//__________________________________________________________
1328Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells)
1329{
1330 // check if the cluster survives some quality cut
1331 //
1332 //
1333 Bool_t isGood=kTRUE;
08ff636b 1334 if(!cluster || !cluster->IsEMCAL()) return kFALSE;
b5078f5d 1335 if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) isGood=kFALSE;
1336 if(!CheckCellFiducialRegion(geom,cluster,cells)) isGood=kFALSE;
1337 if(fRejectExoticCluster && IsExoticCluster(cluster)) isGood=kFALSE;
1338
1339 return isGood;
1340}
1341
b540d03f 1342//__________________________________________________________
bd8c7aef 1343Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1344{
1345 // Given a esd track, return whether the track survive all the cuts
1346
1347 // The different quality parameter are first
1348 // retrieved from the track. then it is found out what cuts the
1349 // track did not survive and finally the cuts are imposed.
1350
1351 UInt_t status = esdTrack->GetStatus();
1352
1353 Int_t nClustersITS = esdTrack->GetITSclusters(0);
1354 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1355
1356 Float_t chi2PerClusterITS = -1;
1357 Float_t chi2PerClusterTPC = -1;
1358 if (nClustersITS!=0)
1359 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1360 if (nClustersTPC!=0)
1361 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
82d09e74 1362
1363
1364 //DCA cuts
1365 Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1366 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1367 SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1368
1369
bd8c7aef 1370 Float_t b[2];
1371 Float_t bCov[3];
1372 esdTrack->GetImpactParameters(b,bCov);
1373 if (bCov[0]<=0 || bCov[2]<=0) {
1374 AliDebug(1, "Estimated b resolution lower or equal zero!");
1375 bCov[0]=0; bCov[2]=0;
1376 }
1377
1378 Float_t dcaToVertexXY = b[0];
1379 Float_t dcaToVertexZ = b[1];
1380 Float_t dcaToVertex = -1;
1381
1382 if (fCutDCAToVertex2D)
1383 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1384 else
1385 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1386
1387 // cut the track?
1388
1389 Bool_t cuts[kNCuts];
1390 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1391
1392 // track quality cuts
1393 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1394 cuts[0]=kTRUE;
1395 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1396 cuts[1]=kTRUE;
1397 if (nClustersTPC<fCutMinNClusterTPC)
1398 cuts[2]=kTRUE;
1399 if (nClustersITS<fCutMinNClusterITS)
1400 cuts[3]=kTRUE;
1401 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1402 cuts[4]=kTRUE;
1403 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1404 cuts[5]=kTRUE;
1405 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1406 cuts[6]=kTRUE;
1407 if (fCutDCAToVertex2D && dcaToVertex > 1)
1408 cuts[7] = kTRUE;
1409 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1410 cuts[8] = kTRUE;
1411 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1412 cuts[9] = kTRUE;
1413
82d09e74 1414 //Require at least one SPD point + anything else in ITS
1415 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1416 cuts[10] = kTRUE;
1417
bd8c7aef 1418 Bool_t cut=kFALSE;
1419 for (Int_t i=0; i<kNCuts; i++)
1420 if (cuts[i]) {cut = kTRUE;}
1421
1422 // cut the track
1423 if (cut)
1424 return kFALSE;
1425 else
1426 return kTRUE;
1427}
1428//__________________________________________________
1429void AliEMCALRecoUtils::InitTrackCuts()
1430{
1431 //Intilize the track cut criteria
82d09e74 1432 //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
bd8c7aef 1433 //Also you can customize the cuts using the setters
82d09e74 1434
1435 //TPC
1436 SetMinNClustersTPC(70);
bd8c7aef 1437 SetMaxChi2PerClusterTPC(4);
1438 SetAcceptKinkDaughters(kFALSE);
82d09e74 1439 SetRequireTPCRefit(kTRUE);
1440
1441 //ITS
1442 SetRequireITSRefit(kTRUE);
1443 SetMaxDCAToVertexZ(2);
1444 SetDCAToVertex2D(kFALSE);
171be15c 1445 SetMaxChi2PerClusterITS(); //which by default sets the value to 1e10.
1446 SetMinNClustersITS();
bd8c7aef 1447}
83bfd77a 1448
b540d03f 1449//___________________________________________________
d9b3567c 1450void AliEMCALRecoUtils::Print(const Option_t *) const
1451{
1452 // Print Parameters
1453
1454 printf("AliEMCALRecoUtils Settings: \n");
1455 printf("Misalignment shifts\n");
2a71e873 1456 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,
1457 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1458 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
d9b3567c 1459 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1460 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
094786cc 1461
1462 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
bd8c7aef 1463
1464 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
bb6f5f0b 1465 printf("Mass hypothesis = %2.3f[GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
bd8c7aef 1466
1467 printf("Track cuts: \n");
bb6f5f0b 1468 printf("AOD track selection mask: %d\n",fAODFilterMask);
bd8c7aef 1469 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1470 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1471 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1472 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1473 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1474
d9b3567c 1475}
96957075 1476
b540d03f 1477//_____________________________________________________________________
96957075 1478void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1479 //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1480 //Do it only once and only if it is requested
1481
1482 if(!fUseTimeCorrectionFactors) return;
1483 if(fTimeCorrectionFactorsSet) return;
1484
1485 printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1486
1487 AliEMCALCalibTimeDepCorrection *corr = new AliEMCALCalibTimeDepCorrection();
1488 corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1489
1490 SwitchOnRecalibration();
1491 for(Int_t ism = 0; ism < 4; ism++){
1492 for(Int_t icol = 0; icol < 48; icol++){
1493 for(Int_t irow = 0; irow < 24; irow++){
1494 Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1495 Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1496 GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1497 //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
1498 // orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1499 // (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1500 }
1501 }
1502 }
1503 fTimeCorrectionFactorsSet = kTRUE;
1504}
1505