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