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