]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - EMCAL/AliEMCALRecoUtils.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
bug fix - update to track cuts by Rongrong
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRecoUtils.cxx
... / ...
CommitLineData
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)
25// Track matching part: Rongrong Ma (Yale)
26
27///////////////////////////////////////////////////////////////////////////////
28// --- standard c ---
29
30// standard C++ includes
31//#include <Riostream.h>
32
33// ROOT includes
34#include <TGeoManager.h>
35#include <TGeoMatrix.h>
36#include <TGeoBBox.h>
37#include <TH2F.h>
38#include <TArrayI.h>
39#include <TArrayF.h>
40#include <TObjArray.h>
41
42// STEER includes
43#include "AliVCluster.h"
44#include "AliVCaloCells.h"
45#include "AliVEvent.h"
46#include "AliLog.h"
47#include "AliPID.h"
48#include "AliESDEvent.h"
49#include "AliAODEvent.h"
50#include "AliESDtrack.h"
51#include "AliAODTrack.h"
52#include "AliExternalTrackParam.h"
53#include "AliESDfriendTrack.h"
54#include "AliTrackerBase.h"
55
56// EMCAL includes
57#include "AliEMCALRecoUtils.h"
58#include "AliEMCALGeometry.h"
59#include "AliEMCALTrack.h"
60#include "AliEMCALCalibTimeDepCorrection.h" // Run dependent
61#include "AliEMCALPIDUtils.h"
62
63ClassImp(AliEMCALRecoUtils)
64
65//______________________________________________
66AliEMCALRecoUtils::AliEMCALRecoUtils():
67 fParticleType(kPhoton), fPosAlgo(kUnchanged), fW0(4.),
68 fNonLinearityFunction(kNoCorrection), fNonLinearThreshold(30),
69 fSmearClusterEnergy(kFALSE), fRandom(),
70 fCellsRecalibrated(kFALSE), fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
71 fTimeRecalibration(kFALSE), fEMCALTimeRecalibrationFactors(),
72 fUseRunCorrectionFactors(kFALSE), fRunCorrectionFactorsSet(kFALSE),
73 fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
74 fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
75 fRejectExoticCluster(kFALSE), fPIDUtils(), fAODFilterMask(32),
76 fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
77 fResidualEta(0x0), fResidualPhi(0x0), fCutEtaPhiSum(kTRUE), fCutEtaPhiSeparate(kFALSE),
78 fCutR(0.1), fCutEta(0.025), fCutPhi(0.05),
79 fMass(0.139), fStep(10),
80 fTrackCutsType(kLooseCut), fCutMinTrackPt(0), fCutMinNClusterTPC(-1),
81 fCutMinNClusterITS(-1), fCutMaxChi2PerClusterTPC(1e10), fCutMaxChi2PerClusterITS(1e10),
82 fCutRequireTPCRefit(kFALSE), fCutRequireITSRefit(kFALSE), fCutAcceptKinkDaughters(kFALSE),
83 fCutMaxDCAToVertexXY(1e10), fCutMaxDCAToVertexZ(1e10), fCutDCAToVertex2D(kFALSE)
84{
85//
86 // Constructor.
87 // Initialize all constant values which have to be used
88 // during Reco algorithm execution
89 //
90
91 //Misalignment matrices
92 for(Int_t i = 0; i < 15 ; i++) {
93 fMisalTransShift[i] = 0.;
94 fMisalRotShift[i] = 0.;
95 }
96
97 //Non linearity
98 for(Int_t i = 0; i < 7 ; i++) fNonLinearityParams[i] = 0.;
99
100 //For kBeamTestCorrected case, but default is no correction
101 fNonLinearityParams[0] = 0.99078;
102 fNonLinearityParams[1] = 0.161499;
103 fNonLinearityParams[2] = 0.655166;
104 fNonLinearityParams[3] = 0.134101;
105 fNonLinearityParams[4] = 163.282;
106 fNonLinearityParams[5] = 23.6904;
107 fNonLinearityParams[6] = 0.978;
108
109 //For kPi0GammaGamma case
110 //fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
111 //fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
112 //fNonLinearityParams[2] = 1.046;
113
114 //Cluster energy smearing
115 fSmearClusterEnergy = kFALSE;
116 fSmearClusterParam[0] = 0.07; // * sqrt E term
117 fSmearClusterParam[1] = 0.00; // * E term
118 fSmearClusterParam[2] = 0.00; // constant
119
120 //Track matching
121 fMatchedTrackIndex = new TArrayI();
122 fMatchedClusterIndex = new TArrayI();
123 fResidualPhi = new TArrayF();
124 fResidualEta = new TArrayF();
125 fPIDUtils = new AliEMCALPIDUtils();
126
127 InitTrackCuts();
128}
129
130//______________________________________________________________________
131AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
132: TNamed(reco),
133 fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0),
134 fNonLinearityFunction(reco.fNonLinearityFunction), fNonLinearThreshold(reco.fNonLinearThreshold),
135 fSmearClusterEnergy(reco.fSmearClusterEnergy), fRandom(),
136 fCellsRecalibrated(reco.fCellsRecalibrated),
137 fRecalibration(reco.fRecalibration), fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
138 fTimeRecalibration(reco.fTimeRecalibration), fEMCALTimeRecalibrationFactors(reco.fEMCALTimeRecalibrationFactors),
139 fUseRunCorrectionFactors(reco.fUseRunCorrectionFactors), fRunCorrectionFactorsSet(reco.fRunCorrectionFactorsSet),
140 fRemoveBadChannels(reco.fRemoveBadChannels), fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
141 fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
142 fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder), fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
143 fRejectExoticCluster(reco.fRejectExoticCluster), fPIDUtils(reco.fPIDUtils),
144 fAODFilterMask(reco.fAODFilterMask),
145 fMatchedTrackIndex( reco.fMatchedTrackIndex? new TArrayI(*reco.fMatchedTrackIndex):0x0),
146 fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
147 fResidualEta( reco.fResidualEta? new TArrayF(*reco.fResidualEta):0x0),
148 fResidualPhi( reco.fResidualPhi? new TArrayF(*reco.fResidualPhi):0x0),
149 fCutEtaPhiSum(reco.fCutEtaPhiSum), fCutEtaPhiSeparate(reco.fCutEtaPhiSeparate),
150 fCutR(reco.fCutR), fCutEta(reco.fCutEta), fCutPhi(reco.fCutPhi),
151 fMass(reco.fMass), fStep(reco.fStep),
152 fTrackCutsType(reco.fTrackCutsType), fCutMinTrackPt(reco.fCutMinTrackPt),
153 fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS),
154 fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
155 fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
156 fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters), fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY),
157 fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ), fCutDCAToVertex2D(reco.fCutDCAToVertex2D)
158{
159 //Copy ctor
160
161 for(Int_t i = 0; i < 15 ; i++) { fMisalRotShift[i] = reco.fMisalRotShift[i] ;
162 fMisalTransShift[i] = reco.fMisalTransShift[i] ; }
163 for(Int_t i = 0; i < 7 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
164 for(Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
165
166}
167
168
169//______________________________________________________________________
170AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
171{
172 //Assignment operator
173
174 if(this == &reco)return *this;
175 ((TNamed *)this)->operator=(reco);
176
177 for(Int_t i = 0; i < 15 ; i++) { fMisalTransShift[i] = reco.fMisalTransShift[i] ;
178 fMisalRotShift[i] = reco.fMisalRotShift[i] ; }
179 for(Int_t i = 0; i < 7 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
180 for(Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
181
182 fParticleType = reco.fParticleType;
183 fPosAlgo = reco.fPosAlgo;
184 fW0 = reco.fW0;
185
186 fNonLinearityFunction = reco.fNonLinearityFunction;
187 fNonLinearThreshold = reco.fNonLinearThreshold;
188 fSmearClusterEnergy = reco.fSmearClusterEnergy;
189
190 fCellsRecalibrated = reco.fCellsRecalibrated;
191 fRecalibration = reco.fRecalibration;
192 fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
193
194 fTimeRecalibration = reco.fTimeRecalibration;
195 fEMCALTimeRecalibrationFactors = reco.fEMCALTimeRecalibrationFactors;
196
197 fUseRunCorrectionFactors = reco.fUseRunCorrectionFactors;
198 fRunCorrectionFactorsSet = reco.fRunCorrectionFactorsSet;
199
200 fRemoveBadChannels = reco.fRemoveBadChannels;
201 fRecalDistToBadChannels = reco.fRecalDistToBadChannels;
202 fEMCALBadChannelMap = reco.fEMCALBadChannelMap;
203
204 fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
205 fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
206 fRejectExoticCluster = reco.fRejectExoticCluster;
207
208 fPIDUtils = reco.fPIDUtils;
209
210 fAODFilterMask = reco.fAODFilterMask;
211
212 fCutEtaPhiSum = reco.fCutEtaPhiSum;
213 fCutEtaPhiSeparate = reco.fCutEtaPhiSeparate;
214 fCutR = reco.fCutR;
215 fCutEta = reco.fCutEta;
216 fCutPhi = reco.fCutPhi;
217 fMass = reco.fMass;
218 fStep = reco.fStep;
219 fRejectExoticCluster = reco.fRejectExoticCluster;
220
221 fTrackCutsType = reco.fTrackCutsType;
222 fCutMinTrackPt = reco.fCutMinTrackPt;
223 fCutMinNClusterTPC = reco.fCutMinNClusterTPC;
224 fCutMinNClusterITS = reco.fCutMinNClusterITS;
225 fCutMaxChi2PerClusterTPC = reco.fCutMaxChi2PerClusterTPC;
226 fCutMaxChi2PerClusterITS = reco.fCutMaxChi2PerClusterITS;
227 fCutRequireTPCRefit = reco.fCutRequireTPCRefit;
228 fCutRequireITSRefit = reco.fCutRequireITSRefit;
229 fCutAcceptKinkDaughters = reco.fCutAcceptKinkDaughters;
230 fCutMaxDCAToVertexXY = reco.fCutMaxDCAToVertexXY;
231 fCutMaxDCAToVertexZ = reco.fCutMaxDCAToVertexZ;
232 fCutDCAToVertex2D = reco.fCutDCAToVertex2D;
233
234 if(reco.fResidualEta){
235 // assign or copy construct
236 if(fResidualEta){
237 *fResidualEta = *reco.fResidualEta;
238 }
239 else fResidualEta = new TArrayF(*reco.fResidualEta);
240 }
241 else{
242 if(fResidualEta)delete fResidualEta;
243 fResidualEta = 0;
244 }
245
246 if(reco.fResidualPhi){
247 // assign or copy construct
248 if(fResidualPhi){
249 *fResidualPhi = *reco.fResidualPhi;
250 }
251 else fResidualPhi = new TArrayF(*reco.fResidualPhi);
252 }
253 else{
254 if(fResidualPhi)delete fResidualPhi;
255 fResidualPhi = 0;
256 }
257
258 if(reco.fMatchedTrackIndex){
259 // assign or copy construct
260 if(fMatchedTrackIndex){
261 *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
262 }
263 else fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
264 }
265 else{
266 if(fMatchedTrackIndex)delete fMatchedTrackIndex;
267 fMatchedTrackIndex = 0;
268 }
269
270 if(reco.fMatchedClusterIndex){
271 // assign or copy construct
272 if(fMatchedClusterIndex){
273 *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
274 }
275 else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
276 }
277 else{
278 if(fMatchedClusterIndex)delete fMatchedClusterIndex;
279 fMatchedClusterIndex = 0;
280 }
281
282 return *this;
283}
284
285
286//__________________________________________________
287AliEMCALRecoUtils::~AliEMCALRecoUtils()
288{
289 //Destructor.
290
291 if(fEMCALRecalibrationFactors) {
292 fEMCALRecalibrationFactors->Clear();
293 delete fEMCALRecalibrationFactors;
294 }
295
296 if(fEMCALTimeRecalibrationFactors) {
297 fEMCALTimeRecalibrationFactors->Clear();
298 delete fEMCALTimeRecalibrationFactors;
299 }
300
301 if(fEMCALBadChannelMap) {
302 fEMCALBadChannelMap->Clear();
303 delete fEMCALBadChannelMap;
304 }
305
306 delete fMatchedTrackIndex ;
307 delete fMatchedClusterIndex ;
308 delete fResidualEta ;
309 delete fResidualPhi ;
310 delete fPIDUtils ;
311
312 InitTrackCuts();
313}
314
315//_______________________________________________________________
316Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells)
317{
318 // Given the list of AbsId of the cluster, get the maximum cell and
319 // check if there are fNCellsFromBorder from the calorimeter border
320
321 if(!cluster){
322 AliInfo("Cluster pointer null!");
323 return kFALSE;
324 }
325
326 //If the distance to the border is 0 or negative just exit accept all clusters
327 if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
328
329 Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
330 Bool_t shared = kFALSE;
331 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
332
333 AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
334 absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
335
336 if(absIdMax==-1) return kFALSE;
337
338 //Check if the cell is close to the borders:
339 Bool_t okrow = kFALSE;
340 Bool_t okcol = kFALSE;
341
342 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
343 AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
344 iSM,ieta,iphi));
345 }
346
347 //Check rows/phi
348 if(iSM < 10){
349 if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
350 }
351 else{
352 if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
353 }
354
355 //Check columns/eta
356 if(!fNoEMCALBorderAtEta0){
357 if(ieta > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
358 }
359 else{
360 if(iSM%2==0){
361 if(ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
362 }
363 else {
364 if(ieta < 48-fNCellsFromEMCALBorder) okcol = kTRUE;
365 }
366 }//eta 0 not checked
367
368 AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
369 fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
370
371 if (okcol && okrow) {
372 //printf("Accept\n");
373 return kTRUE;
374 }
375 else {
376 //printf("Reject\n");
377 AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
378 return kFALSE;
379 }
380
381}
382
383
384//_________________________________________________________________________________________________________
385Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, const Int_t nCells){
386 // Check that in the cluster cells, there is no bad channel of those stored
387 // in fEMCALBadChannelMap or fPHOSBadChannelMap
388
389 if(!fRemoveBadChannels) return kFALSE;
390 if(!fEMCALBadChannelMap) return kFALSE;
391
392 Int_t icol = -1;
393 Int_t irow = -1;
394 Int_t imod = -1;
395 for(Int_t iCell = 0; iCell<nCells; iCell++){
396
397 //Get the column and row
398 Int_t iTower = -1, iIphi = -1, iIeta = -1;
399 geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
400 if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
401 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
402 if(GetEMCALChannelStatus(imod, icol, irow)){
403 AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
404 return kTRUE;
405 }
406
407 }// cell cluster loop
408
409 return kFALSE;
410
411}
412
413//_________________________________________________
414Bool_t AliEMCALRecoUtils::IsExoticCluster(AliVCluster *cluster) const {
415 // Check if the cluster has high energy but small number of cells
416 // The criteria comes from Gustavo's study
417 //
418
419 if(!cluster){
420 AliInfo("Cluster pointer null!");
421 return kFALSE;
422 }
423
424 Int_t nc = cluster->GetNCells() ;
425
426 if ( nc > 8 ) return kFALSE ; // Good cluster, needed for 3x3 clusterizer
427 else if ( nc < 1 + cluster->E()/3. ) return kTRUE ; // Bad cluster
428 else return kFALSE ; // Good cluster
429
430}
431
432//__________________________________________________
433Float_t AliEMCALRecoUtils::SmearClusterEnergy(AliVCluster* cluster) {
434
435 //In case of MC analysis, smear energy to match resolution/calibration in real data
436
437 if(!cluster){
438 AliInfo("Cluster pointer null!");
439 return 0;
440 }
441
442 Float_t energy = cluster->E() ;
443 Float_t rdmEnergy = energy ;
444 if(fSmearClusterEnergy){
445 rdmEnergy = fRandom.Gaus(energy,fSmearClusterParam[0] * TMath::Sqrt(energy) +
446 fSmearClusterParam[1] * energy +
447 fSmearClusterParam[2] );
448 AliDebug(2, Form("Energy: original %f, smeared %f\n", energy, rdmEnergy));
449 }
450
451 return rdmEnergy ;
452
453}
454
455//__________________________________________________
456Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
457// Correct cluster energy from non linearity functions
458
459 if(!cluster){
460 AliInfo("Cluster pointer null!");
461 return 0;
462 }
463
464 Float_t energy = cluster->E();
465
466 switch (fNonLinearityFunction) {
467
468 case kPi0MC:
469 {
470 //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
471 //Double_t fNonLinearityParams[0] = 1.014;
472 //Double_t fNonLinearityParams[1] = -0.03329;
473 //Double_t fNonLinearityParams[2] = -0.3853;
474 //Double_t fNonLinearityParams[3] = 0.5423;
475 //Double_t fNonLinearityParams[4] = -0.4335;
476 energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
477 ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
478 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
479 break;
480 }
481
482 case kPi0GammaGamma:
483 {
484 //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
485 //Double_t fNonLinearityParams[0] = 1.04;
486 //Double_t fNonLinearityParams[1] = -0.1445;
487 //Double_t fNonLinearityParams[2] = 1.046;
488 energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
489 break;
490 }
491
492 case kPi0GammaConversion:
493 {
494 //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
495 //fNonLinearityParams[0] = 0.139393/0.1349766;
496 //fNonLinearityParams[1] = 0.0566186;
497 //fNonLinearityParams[2] = 0.982133;
498 energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
499
500 break;
501 }
502
503 case kBeamTest:
504 {
505 //From beam test, Alexei's results, for different ZS thresholds
506 // th=30 MeV; th = 45 MeV; th = 75 MeV
507 //fNonLinearityParams[0] = 1.007; 1.003; 1.002
508 //fNonLinearityParams[1] = 0.894; 0.719; 0.797
509 //fNonLinearityParams[2] = 0.246; 0.334; 0.358
510 //Rescale the param[0] with 1.03
511 energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
512
513 break;
514 }
515
516 case kBeamTestCorrected:
517 {
518 //From beam test, corrected for material between beam and EMCAL
519 //fNonLinearityParams[0] = 0.99078
520 //fNonLinearityParams[1] = 0.161499;
521 //fNonLinearityParams[2] = 0.655166;
522 //fNonLinearityParams[3] = 0.134101;
523 //fNonLinearityParams[4] = 163.282;
524 //fNonLinearityParams[5] = 23.6904;
525 //fNonLinearityParams[6] = 0.978;
526 energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
527
528 break;
529 }
530
531 case kNoCorrection:
532 AliDebug(2,"No correction on the energy\n");
533 break;
534
535 }
536
537 return energy;
538
539}
540//__________________________________________________
541void AliEMCALRecoUtils::InitNonLinearityParam()
542{
543 //Initialising Non Linearity Parameters
544
545 if(fNonLinearityFunction == kPi0MC)
546 {
547 fNonLinearityParams[0] = 1.014;
548 fNonLinearityParams[1] = -0.03329;
549 fNonLinearityParams[2] = -0.3853;
550 fNonLinearityParams[3] = 0.5423;
551 fNonLinearityParams[4] = -0.4335;
552 }
553
554 if(fNonLinearityFunction == kPi0GammaGamma)
555 {
556 fNonLinearityParams[0] = 1.04;
557 fNonLinearityParams[1] = -0.1445;
558 fNonLinearityParams[2] = 1.046;
559 }
560
561 if(fNonLinearityFunction == kPi0GammaConversion)
562 {
563 fNonLinearityParams[0] = 0.139393;
564 fNonLinearityParams[1] = 0.0566186;
565 fNonLinearityParams[2] = 0.982133;
566 }
567
568 if(fNonLinearityFunction == kBeamTest)
569 {
570 if(fNonLinearThreshold == 30)
571 {
572 fNonLinearityParams[0] = 1.007;
573 fNonLinearityParams[1] = 0.894;
574 fNonLinearityParams[2] = 0.246;
575 }
576 if(fNonLinearThreshold == 45)
577 {
578 fNonLinearityParams[0] = 1.003;
579 fNonLinearityParams[1] = 0.719;
580 fNonLinearityParams[2] = 0.334;
581 }
582 if(fNonLinearThreshold == 75)
583 {
584 fNonLinearityParams[0] = 1.002;
585 fNonLinearityParams[1] = 0.797;
586 fNonLinearityParams[2] = 0.358;
587 }
588 }
589
590 if(fNonLinearityFunction == kBeamTestCorrected)
591 {
592 fNonLinearityParams[0] = 0.99078;
593 fNonLinearityParams[1] = 0.161499;
594 fNonLinearityParams[2] = 0.655166;
595 fNonLinearityParams[3] = 0.134101;
596 fNonLinearityParams[4] = 163.282;
597 fNonLinearityParams[5] = 23.6904;
598 fNonLinearityParams[6] = 0.978;
599 }
600}
601
602//__________________________________________________
603Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
604{
605 //Calculate shower depth for a given cluster energy and particle type
606
607 // parameters
608 Float_t x0 = 1.31;
609 Float_t ecr = 8;
610 Float_t depth = 0;
611
612 switch ( iParticle )
613 {
614 case kPhoton:
615 depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
616 break;
617
618 case kElectron:
619 depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
620 break;
621
622 case kHadron:
623 // hadron
624 // boxes anc. here
625 if(gGeoManager){
626 gGeoManager->cd("ALIC_1/XEN1_1");
627 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
628 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
629 if(geoSM){
630 TGeoVolume *geoSMVol = geoSM->GetVolume();
631 TGeoShape *geoSMShape = geoSMVol->GetShape();
632 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
633 if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
634 else AliFatal("Null GEANT box");
635 }else AliFatal("NULL GEANT node matrix");
636 }
637 else{//electron
638 depth = x0 * (TMath::Log(energy*1000 / ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
639 }
640
641 break;
642
643 default://photon
644 depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
645 }
646
647 return depth;
648
649}
650
651//__________________________________________________
652void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
653 Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
654{
655 //For a given CaloCluster gets the absId of the cell
656 //with maximum energy deposit.
657
658 Double_t eMax = -1.;
659 Double_t eCell = -1.;
660 Float_t fraction = 1.;
661 Float_t recalFactor = 1.;
662 Int_t cellAbsId = -1 ;
663
664 Int_t iTower = -1;
665 Int_t iIphi = -1;
666 Int_t iIeta = -1;
667 Int_t iSupMod0= -1;
668
669 if(!clu){
670 AliInfo("Cluster pointer null!");
671 absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1;
672 return;
673 }
674
675 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
676 cellAbsId = clu->GetCellAbsId(iDig);
677 fraction = clu->GetCellAmplitudeFraction(iDig);
678 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
679 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
680 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
681 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
682 if(iDig==0) iSupMod0=iSupMod;
683 else if(iSupMod0!=iSupMod) {
684 shared = kTRUE;
685 //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
686 }
687 if(IsRecalibrationOn()) {
688 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
689 }
690 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
691 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
692 if(eCell > eMax) {
693 eMax = eCell;
694 absId = cellAbsId;
695 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
696 }
697 }// cell loop
698
699 //Get from the absid the supermodule, tower and eta/phi numbers
700 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
701 //Gives SuperModule and Tower numbers
702 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
703 iIphi, iIeta,iphi,ieta);
704 //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
705 //printf("Max end---\n");
706
707}
708
709//________________________________________________________________
710void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
711 //Init EMCAL recalibration factors
712 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
713 //In order to avoid rewriting the same histograms
714 Bool_t oldStatus = TH1::AddDirectoryStatus();
715 TH1::AddDirectory(kFALSE);
716
717 fEMCALRecalibrationFactors = new TObjArray(10);
718 for (int i = 0; i < 10; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
719 //Init the histograms with 1
720 for (Int_t sm = 0; sm < 10; sm++) {
721 for (Int_t i = 0; i < 48; i++) {
722 for (Int_t j = 0; j < 24; j++) {
723 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
724 }
725 }
726 }
727 fEMCALRecalibrationFactors->SetOwner(kTRUE);
728 fEMCALRecalibrationFactors->Compress();
729
730 //In order to avoid rewriting the same histograms
731 TH1::AddDirectory(oldStatus);
732}
733
734//________________________________________________________________
735void AliEMCALRecoUtils::InitEMCALTimeRecalibrationFactors(){
736 //Init EMCAL recalibration factors
737 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
738 //In order to avoid rewriting the same histograms
739 Bool_t oldStatus = TH1::AddDirectoryStatus();
740 TH1::AddDirectory(kFALSE);
741
742 fEMCALTimeRecalibrationFactors = new TObjArray(4);
743 for (int i = 0; i < 4; i++)
744 fEMCALTimeRecalibrationFactors->Add(new TH1F(Form("hAllTimeAvBC%d",i),
745 Form("hAllTimeAvBC%d",i),
746 48*24*10,0.,48*24*10) );
747 //Init the histograms with 1
748 for (Int_t bc = 0; bc < 4; bc++) {
749 for (Int_t i = 0; i < 48*24*10; i++)
750 SetEMCALChannelTimeRecalibrationFactor(bc,i,0.);
751 }
752
753 fEMCALTimeRecalibrationFactors->SetOwner(kTRUE);
754 fEMCALTimeRecalibrationFactors->Compress();
755
756 //In order to avoid rewriting the same histograms
757 TH1::AddDirectory(oldStatus);
758}
759
760//________________________________________________________________
761void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
762 //Init EMCAL bad channels map
763 AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
764 //In order to avoid rewriting the same histograms
765 Bool_t oldStatus = TH1::AddDirectoryStatus();
766 TH1::AddDirectory(kFALSE);
767
768 fEMCALBadChannelMap = new TObjArray(10);
769 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
770 for (int i = 0; i < 10; i++) {
771 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
772 }
773
774 fEMCALBadChannelMap->SetOwner(kTRUE);
775 fEMCALBadChannelMap->Compress();
776
777 //In order to avoid rewriting the same histograms
778 TH1::AddDirectory(oldStatus);
779}
780
781//________________________________________________________________
782void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells, const Int_t bc){
783 // Recalibrate the cluster energy and Time, considering the recalibration map
784 // and the energy of the cells and time that compose the cluster.
785 // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
786
787 if(!cluster){
788 AliInfo("Cluster pointer null!");
789 return;
790 }
791
792 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
793 UShort_t * index = cluster->GetCellsAbsId() ;
794 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
795 Int_t ncells = cluster->GetNCells();
796
797 //Initialize some used variables
798 Float_t energy = 0;
799 Int_t absId =-1;
800 Int_t icol =-1, irow =-1, imod=1;
801 Float_t factor = 1, frac = 0;
802 Int_t absIdMax = -1;
803 Float_t emax = 0;
804
805 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
806 for(Int_t icell = 0; icell < ncells; icell++){
807 absId = index[icell];
808 frac = fraction[icell];
809 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
810
811 if(!fCellsRecalibrated && IsRecalibrationOn()){
812
813 // Energy
814 Int_t iTower = -1, iIphi = -1, iIeta = -1;
815 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
816 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
817 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
818 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
819
820 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
821 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
822
823 }
824
825 energy += cells->GetCellAmplitude(absId)*factor*frac;
826
827 if(emax < cells->GetCellAmplitude(absId)*factor*frac){
828 emax = cells->GetCellAmplitude(absId)*factor*frac;
829 absIdMax = absId;
830 }
831
832 }
833
834 cluster->SetE(energy);
835
836 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
837
838 // Recalculate time of cluster only for ESDs
839 if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName()))){
840
841 // Time
842 Double_t weightedTime = 0;
843 Double_t weight = 0;
844 Double_t weightTot = 0;
845 Double_t maxcellTime = 0;
846 for(Int_t icell = 0; icell < ncells; icell++){
847 absId = index[icell];
848 frac = fraction[icell];
849 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
850
851 Double_t celltime = cells->GetCellTime(absId);
852 RecalibrateCellTime(absId, bc, celltime);
853 if(absId == absIdMax) maxcellTime = celltime;
854
855 if(!fCellsRecalibrated){
856
857 Int_t iTower = -1, iIphi = -1, iIeta = -1;
858 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
859 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
860 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
861 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
862
863 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
864 imod,icol,irow,frac,factor,cells->GetCellTime(absId)));
865
866 }
867
868 weight = GetCellWeight(cells->GetCellAmplitude(absId)*factor*frac , energy );
869 weightTot += weight;
870 weightedTime += celltime * weight;
871
872 }
873
874 if(weightTot > 0)
875 cluster->SetTOF(weightedTime/weightTot);
876 else
877 cluster->SetTOF(maxcellTime);
878
879 }
880}
881
882//________________________________________________________________
883void AliEMCALRecoUtils::RecalibrateCells(AliEMCALGeometry* geom, AliVCaloCells * cells, Int_t bc){
884 // Recalibrate the cells time and energy, considering the recalibration map and the energy
885 // of the cells that compose the cluster.
886 // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
887
888 if(!IsRecalibrationOn()) return;
889
890 if(!cells){
891 AliInfo("Cells pointer null!");
892 return;
893 }
894
895 fCellsRecalibrated = kTRUE;
896
897 Int_t absId =-1;
898 Int_t icol =-1, irow =-1, imod = 1;
899 Int_t iTower =-1, iIeta =-1, iIphi =-1;
900
901 Int_t nEMcell = cells->GetNumberOfCells() ;
902
903 for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
904
905 absId = cells->GetCellNumber(iCell);
906
907 // Energy
908 Float_t factor = 1;
909 if(IsRecalibrationOn()){
910 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
911 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
912 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
913 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
914 }
915
916 Float_t cellE = cells->GetAmplitude(iCell) * factor ;
917
918 //Time
919 Double_t celltime = cells->GetCellTime(absId);
920 RecalibrateCellTime(absId, bc, celltime);
921
922 //Set new values
923 cells->SetCell(iCell,cells->GetCellNumber(iCell),cellE, celltime);
924
925 }
926
927}
928
929//________________________________________________________________
930void AliEMCALRecoUtils::RecalibrateCellTime(const Int_t absId, const Int_t bc, Double_t & celltime){
931 // Recalibrate time of cell with absID considering the recalibration map
932 // bc= bunch crossing number returned by esdevent->GetBunchCrossNumber();
933
934 if(!fCellsRecalibrated && IsTimeRecalibrationOn()){
935// printf("cell time org %g, ",celltime);
936
937 Double_t timeBCoffset = 0.;
938 if( bc%4 ==0 || bc%4==1) timeBCoffset = 100.*1.e-9; //in ns
939
940 Double_t celloffset = GetEMCALChannelTimeRecalibrationFactor(bc%4,absId)*1.e-9;
941
942// printf("absId %d, time %f bc %d-%d: bc0 %f, bc1 %f, bc2 %f, bc3 %f \n", absId, celltime*1.e9,bc, bc%4,
943// GetEMCALChannelTimeRecalibrationFactor(0,absId),GetEMCALChannelTimeRecalibrationFactor(1,absId),
944// GetEMCALChannelTimeRecalibrationFactor(2,absId),GetEMCALChannelTimeRecalibrationFactor(3,absId));
945
946 celltime -= timeBCoffset ;
947 celltime -= celloffset ;
948// printf("new %g\n",celltime);
949 }
950
951}
952
953//__________________________________________________
954void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
955{
956 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
957
958 if(!clu){
959 AliInfo("Cluster pointer null!");
960 return;
961 }
962
963 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
964 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
965 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
966
967}
968
969//__________________________________________________
970void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
971{
972 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
973 // The algorithm is a copy of what is done in AliEMCALRecPoint
974
975 Double_t eCell = 0.;
976 Float_t fraction = 1.;
977 Float_t recalFactor = 1.;
978
979 Int_t absId = -1;
980 Int_t iTower = -1, iIphi = -1, iIeta = -1;
981 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
982 Float_t weight = 0., totalWeight=0.;
983 Float_t newPos[3] = {0,0,0};
984 Double_t pLocal[3], pGlobal[3];
985 Bool_t shared = kFALSE;
986
987 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
988 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
989 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
990
991 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
992
993 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
994
995 absId = clu->GetCellAbsId(iDig);
996 fraction = clu->GetCellAmplitudeFraction(iDig);
997 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
998
999 if(!fCellsRecalibrated){
1000
1001 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
1002 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
1003
1004 if(IsRecalibrationOn()) {
1005 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
1006 }
1007 }
1008
1009 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1010
1011 weight = GetCellWeight(eCell,clEnergy);
1012 totalWeight += weight;
1013
1014 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
1015 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
1016 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
1017 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
1018
1019 for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
1020
1021 }// cell loop
1022
1023 if(totalWeight>0){
1024 for(int i=0; i<3; i++ ) newPos[i] /= totalWeight;
1025 }
1026
1027 //Float_t pos[]={0,0,0};
1028 //clu->GetPosition(pos);
1029 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
1030 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1031
1032 if(iSupModMax > 1) {//sector 1
1033 newPos[0] +=fMisalTransShift[3];//-=3.093;
1034 newPos[1] +=fMisalTransShift[4];//+=6.82;
1035 newPos[2] +=fMisalTransShift[5];//+=1.635;
1036 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
1037
1038 }
1039 else {//sector 0
1040 newPos[0] +=fMisalTransShift[0];//+=1.134;
1041 newPos[1] +=fMisalTransShift[1];//+=8.2;
1042 newPos[2] +=fMisalTransShift[2];//+=1.197;
1043 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
1044
1045 }
1046 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
1047
1048 clu->SetPosition(newPos);
1049
1050}
1051
1052//__________________________________________________
1053void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
1054{
1055 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
1056 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
1057
1058 Double_t eCell = 1.;
1059 Float_t fraction = 1.;
1060 Float_t recalFactor = 1.;
1061
1062 Int_t absId = -1;
1063 Int_t iTower = -1;
1064 Int_t iIphi = -1, iIeta = -1;
1065 Int_t iSupMod = -1, iSupModMax = -1;
1066 Int_t iphi = -1, ieta =-1;
1067 Bool_t shared = kFALSE;
1068
1069 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
1070 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
1071 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
1072
1073 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
1074 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
1075 Int_t startingSM = -1;
1076
1077 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
1078 absId = clu->GetCellAbsId(iDig);
1079 fraction = clu->GetCellAmplitudeFraction(iDig);
1080 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1081
1082 if (iDig==0) startingSM = iSupMod;
1083 else if(iSupMod != startingSM) areInSameSM = kFALSE;
1084
1085 eCell = cells->GetCellAmplitude(absId);
1086
1087 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
1088 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1089
1090 if(!fCellsRecalibrated){
1091
1092 if(IsRecalibrationOn()) {
1093
1094 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1095
1096 }
1097 }
1098
1099 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
1100
1101 weight = GetCellWeight(eCell,clEnergy);
1102 if(weight < 0) weight = 0;
1103 totalWeight += weight;
1104 weightedCol += ieta*weight;
1105 weightedRow += iphi*weight;
1106
1107 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
1108
1109 }// cell loop
1110
1111 Float_t xyzNew[]={0.,0.,0.};
1112 if(areInSameSM == kTRUE) {
1113 //printf("In Same SM\n");
1114 weightedCol = weightedCol/totalWeight;
1115 weightedRow = weightedRow/totalWeight;
1116 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
1117 }
1118 else {
1119 //printf("In Different SM\n");
1120 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
1121 }
1122
1123 clu->SetPosition(xyzNew);
1124
1125}
1126
1127//____________________________________________________________________________
1128void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
1129
1130 //re-evaluate distance to bad channel with updated bad map
1131
1132 if(!fRecalDistToBadChannels) return;
1133
1134 if(!cluster){
1135 AliInfo("Cluster pointer null!");
1136 return;
1137 }
1138
1139 //Get channels map of the supermodule where the cluster is.
1140 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
1141 Bool_t shared = kFALSE;
1142 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
1143 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
1144
1145 Int_t dRrow, dRcol;
1146 Float_t minDist = 10000.;
1147 Float_t dist = 0.;
1148
1149 //Loop on tower status map
1150 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
1151 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
1152 //Check if tower is bad.
1153 if(hMap->GetBinContent(icol,irow)==0) continue;
1154 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
1155 // iSupMod,icol, irow, icolM,irowM);
1156
1157 dRrow=TMath::Abs(irowM-irow);
1158 dRcol=TMath::Abs(icolM-icol);
1159 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
1160 if(dist < minDist){
1161 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
1162 minDist = dist;
1163 }
1164
1165 }
1166 }
1167
1168 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
1169 if (shared) {
1170 TH2D* hMap2 = 0;
1171 Int_t iSupMod2 = -1;
1172
1173 //The only possible combinations are (0,1), (2,3) ... (8,9)
1174 if(iSupMod%2) iSupMod2 = iSupMod-1;
1175 else iSupMod2 = iSupMod+1;
1176 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
1177
1178 //Loop on tower status map of second super module
1179 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
1180 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
1181 //Check if tower is bad.
1182 if(hMap2->GetBinContent(icol,irow)==0) continue;
1183 //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",
1184 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
1185
1186 dRrow=TMath::Abs(irow-irowM);
1187
1188 if(iSupMod%2) {
1189 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
1190 }
1191 else {
1192 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
1193 }
1194
1195 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
1196 if(dist < minDist) minDist = dist;
1197
1198 }
1199 }
1200
1201 }// shared cluster in 2 SuperModules
1202
1203 AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
1204 cluster->SetDistanceToBadChannel(minDist);
1205
1206}
1207
1208//____________________________________________________________________________
1209void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
1210
1211 //re-evaluate identification parameters with bayesian
1212
1213 if(!cluster){
1214 AliInfo("Cluster pointer null!");
1215 return;
1216 }
1217
1218 if ( cluster->GetM02() != 0)
1219 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
1220
1221 Float_t pidlist[AliPID::kSPECIESN+1];
1222 for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
1223
1224 cluster->SetPID(pidlist);
1225
1226}
1227
1228//____________________________________________________________________________
1229void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
1230{
1231 // Calculates new center of gravity in the local EMCAL-module coordinates
1232 // and tranfers into global ALICE coordinates
1233 // Calculates Dispersion and main axis
1234
1235 if(!cluster){
1236 AliInfo("Cluster pointer null!");
1237 return;
1238 }
1239
1240 Int_t nstat = 0;
1241 Float_t wtot = 0. ;
1242 Double_t eCell = 0.;
1243 Float_t fraction = 1.;
1244 Float_t recalFactor = 1.;
1245
1246 Int_t iSupMod = -1;
1247 Int_t iTower = -1;
1248 Int_t iIphi = -1;
1249 Int_t iIeta = -1;
1250 Int_t iphi = -1;
1251 Int_t ieta = -1;
1252 Double_t etai = -1.;
1253 Double_t phii = -1.;
1254
1255 Double_t w = 0.;
1256 Double_t d = 0.;
1257 Double_t dxx = 0.;
1258 Double_t dzz = 0.;
1259 Double_t dxz = 0.;
1260 Double_t xmean = 0.;
1261 Double_t zmean = 0.;
1262
1263 //Loop on cells
1264 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
1265
1266 //Get from the absid the supermodule, tower and eta/phi numbers
1267 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1268 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1269
1270 //Get the cell energy, if recalibration is on, apply factors
1271 fraction = cluster->GetCellAmplitudeFraction(iDigit);
1272 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1273
1274 if(!fCellsRecalibrated){
1275
1276 if(IsRecalibrationOn()) {
1277 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1278 }
1279
1280 }
1281
1282 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1283
1284 if(cluster->E() > 0 && eCell > 0){
1285
1286 w = GetCellWeight(eCell,cluster->E());
1287
1288 etai=(Double_t)ieta;
1289 phii=(Double_t)iphi;
1290 if(w > 0.0) {
1291 wtot += w ;
1292 nstat++;
1293 //Shower shape
1294 dxx += w * etai * etai ;
1295 xmean+= w * etai ;
1296 dzz += w * phii * phii ;
1297 zmean+= w * phii ;
1298 dxz += w * etai * phii ;
1299 }
1300 }
1301 else
1302 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1303 }//cell loop
1304
1305 //Normalize to the weight
1306 if (wtot > 0) {
1307 xmean /= wtot ;
1308 zmean /= wtot ;
1309 }
1310 else
1311 AliError(Form("Wrong weight %f\n", wtot));
1312
1313 //Calculate dispersion
1314 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
1315
1316 //Get from the absid the supermodule, tower and eta/phi numbers
1317 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1318 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1319
1320 //Get the cell energy, if recalibration is on, apply factors
1321 fraction = cluster->GetCellAmplitudeFraction(iDigit);
1322 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1323 if(IsRecalibrationOn()) {
1324 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1325 }
1326 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1327
1328 if(cluster->E() > 0 && eCell > 0){
1329
1330 w = GetCellWeight(eCell,cluster->E());
1331
1332 etai=(Double_t)ieta;
1333 phii=(Double_t)iphi;
1334 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
1335 }
1336 else
1337 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1338 }// cell loop
1339
1340 //Normalize to the weigth and set shower shape parameters
1341 if (wtot > 0 && nstat > 1) {
1342 d /= wtot ;
1343 dxx /= wtot ;
1344 dzz /= wtot ;
1345 dxz /= wtot ;
1346 dxx -= xmean * xmean ;
1347 dzz -= zmean * zmean ;
1348 dxz -= xmean * zmean ;
1349 cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1350 cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1351 }
1352 else{
1353 d=0. ;
1354 cluster->SetM20(0.) ;
1355 cluster->SetM02(0.) ;
1356 }
1357
1358 if (d>=0)
1359 cluster->SetDispersion(TMath::Sqrt(d)) ;
1360 else
1361 cluster->SetDispersion(0) ;
1362}
1363
1364//____________________________________________________________________________
1365void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr, AliEMCALGeometry *geom)
1366{
1367 //This function should be called before the cluster loop
1368 //Before call this function, please recalculate the cluster positions
1369 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
1370 //Store matched cluster indexes and residuals
1371
1372 fMatchedTrackIndex->Reset();
1373 fMatchedClusterIndex->Reset();
1374 fResidualPhi->Reset();
1375 fResidualEta->Reset();
1376
1377 fMatchedTrackIndex->Set(500);
1378 fMatchedClusterIndex->Set(500);
1379 fResidualPhi->Set(500);
1380 fResidualEta->Set(500);
1381
1382 AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1383 AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
1384
1385 Int_t matched=0;
1386 Double_t cv[21];
1387 for (Int_t i=0; i<21;i++) cv[i]=0;
1388 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1389 {
1390 AliExternalTrackParam *trackParam = 0;
1391
1392 //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
1393 if(esdevent)
1394 {
1395 AliESDtrack *esdTrack = esdevent->GetTrack(itr);
1396 if(!esdTrack || !IsAccepted(esdTrack)) continue;
1397 if(esdTrack->Pt()<fCutMinTrackPt) continue;
1398 trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
1399 }
1400
1401 //If the input event is AOD, the starting point for extrapolation is at vertex
1402 //AOD tracks are selected according to its bit.
1403 else if(aodevent)
1404 {
1405 AliAODTrack *aodTrack = aodevent->GetTrack(itr);
1406 if(!aodTrack) continue;
1407 if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
1408 if(aodTrack->Pt()<fCutMinTrackPt) continue;
1409 Double_t pos[3],mom[3];
1410 aodTrack->GetXYZ(pos);
1411 aodTrack->GetPxPyPz(mom);
1412 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()));
1413 trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1414 }
1415
1416 //Return if the input data is not "AOD" or "ESD"
1417 else
1418 {
1419 printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
1420 return;
1421 }
1422
1423 if(!trackParam) continue;
1424
1425 Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
1426 Int_t index = -1;
1427 if(!clusterArr){// get clusters from event
1428 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1429 {
1430 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1431 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1432 AliExternalTrackParam trkPamTmp(*trackParam);//Retrieve the starting point every time before the extrapolation
1433 Float_t tmpEta=-999, tmpPhi=-999;
1434 if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1435 if(fCutEtaPhiSum)
1436 {
1437 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1438 if(tmpR<dRMax)
1439 {
1440 dRMax=tmpR;
1441 dEtaMax=tmpEta;
1442 dPhiMax=tmpPhi;
1443 index=icl;
1444 }
1445 }
1446 else if(fCutEtaPhiSeparate)
1447 {
1448 if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1449 {
1450 dEtaMax = tmpEta;
1451 dPhiMax = tmpPhi;
1452 index=icl;
1453 }
1454 }
1455 else
1456 {
1457 printf("Error: please specify your cut criteria\n");
1458 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1459 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1460 if(aodevent && trackParam) delete trackParam;
1461 return;
1462 }
1463 }//cluster loop
1464 }
1465 else { // external cluster array, not from ESD event
1466 for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1467 {
1468 AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
1469 if(!cluster){
1470 AliInfo("Cluster not found!!!");
1471 continue;
1472 }
1473 if(!cluster->IsEMCAL()) continue;
1474 AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
1475 Float_t tmpEta=-999, tmpPhi=-999;
1476 if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1477 if(fCutEtaPhiSum)
1478 {
1479 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1480 if(tmpR<dRMax)
1481 {
1482 dRMax=tmpR;
1483 dEtaMax=tmpEta;
1484 dPhiMax=tmpPhi;
1485 index=icl;
1486 }
1487 }
1488 else if(fCutEtaPhiSeparate)
1489 {
1490 if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1491 {
1492 dEtaMax = tmpEta;
1493 dPhiMax = tmpPhi;
1494 index=icl;
1495 }
1496 }
1497 else
1498 {
1499 printf("Error: please specify your cut criteria\n");
1500 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1501 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1502 if(aodevent && trackParam) delete trackParam;
1503 return;
1504 }
1505 }//cluster loop
1506 }// external list of clusters
1507
1508 if(index>-1)
1509 {
1510 fMatchedTrackIndex ->AddAt(itr,matched);
1511 fMatchedClusterIndex->AddAt(index,matched);
1512 fResidualEta ->AddAt(dEtaMax,matched);
1513 fResidualPhi ->AddAt(dPhiMax,matched);
1514 matched++;
1515 }
1516 if(aodevent && trackParam) delete trackParam;
1517 }//track loop
1518
1519 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1520
1521 fMatchedTrackIndex ->Set(matched);
1522 fMatchedClusterIndex->Set(matched);
1523 fResidualPhi ->Set(matched);
1524 fResidualEta ->Set(matched);
1525}
1526
1527//________________________________________________________________________________
1528Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom)
1529{
1530 //
1531 // This function returns the index of matched cluster to input track
1532 // Returns -1 if no match is found
1533
1534
1535 Float_t dRMax = fCutR, dEtaMax = fCutEta, dPhiMax = fCutPhi;
1536 Int_t index = -1;
1537
1538 AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
1539
1540 if(!trackParam) return index;
1541 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1542 {
1543 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1544 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1545 AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
1546 Float_t tmpEta=-999, tmpPhi=-999;
1547 if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1548 if(fCutEtaPhiSum)
1549 {
1550 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1551 if(tmpR<dRMax)
1552 {
1553 dRMax=tmpR;
1554 dEtaMax=tmpEta;
1555 dPhiMax=tmpPhi;
1556 index=icl;
1557 }
1558 }
1559 else if(fCutEtaPhiSeparate)
1560 {
1561 if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1562 {
1563 dEtaMax = tmpEta;
1564 dPhiMax = tmpPhi;
1565 index=icl;
1566 }
1567 }
1568 else
1569 {
1570 printf("Error: please specify your cut criteria\n");
1571 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1572 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1573 return -1;
1574 }
1575 }//cluster loop
1576 return index;
1577}
1578
1579//________________________________________________________________________________
1580Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpEta, Float_t &tmpPhi)
1581{
1582 //
1583 //Return the residual by extrapolating a track to a cluster
1584 //
1585 if(!trkParam || !cluster) return kFALSE;
1586 Float_t clsPos[3];
1587 Double_t trkPos[3];
1588 cluster->GetPosition(clsPos); //Has been recalculated
1589 TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
1590 Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
1591 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
1592 if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kTRUE, 0.8, -1)) return kFALSE;
1593 trkParam->GetXYZ(trkPos); //Get the extrapolated global position
1594
1595 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1596 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
1597
1598 // track cluster matching
1599 tmpPhi = clsPosVec.DeltaPhi(trkPosVec); // tmpPhi is between -pi and pi
1600 tmpEta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
1601
1602 return kTRUE;
1603}
1604
1605//________________________________________________________________________________
1606void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
1607{
1608 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1609 //Get the residuals dEta and dPhi for this cluster to the closest track
1610 //Works with ESDs and AODs
1611
1612 if( FindMatchedPosForCluster(clsIndex) >= 999 )
1613 {
1614 AliDebug(2,"No matched tracks found!\n");
1615 dEta=999.;
1616 dPhi=999.;
1617 return;
1618 }
1619 dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
1620 dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
1621}
1622//________________________________________________________________________________
1623void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
1624{
1625 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1626 //Get the residuals dEta and dPhi for this track to the closest cluster
1627 //Works with ESDs and AODs
1628
1629 if( FindMatchedPosForTrack(trkIndex) >= 999 )
1630 {
1631 AliDebug(2,"No matched cluster found!\n");
1632 dEta=999.;
1633 dPhi=999.;
1634 return;
1635 }
1636 dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
1637 dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
1638}
1639
1640//__________________________________________________________
1641Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
1642{
1643 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1644 //Get the index of matched track to this cluster
1645 //Works with ESDs and AODs
1646
1647 if(IsClusterMatched(clsIndex))
1648 return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
1649 else
1650 return -1;
1651}
1652
1653//__________________________________________________________
1654Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
1655{
1656 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1657 //Get the index of matched cluster to this track
1658 //Works with ESDs and AODs
1659
1660 if(IsTrackMatched(trkIndex))
1661 return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
1662 else
1663 return -1;
1664}
1665
1666//__________________________________________________
1667Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
1668{
1669 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1670 //Returns if the cluster has a match
1671 if(FindMatchedPosForCluster(clsIndex) < 999)
1672 return kTRUE;
1673 else
1674 return kFALSE;
1675}
1676
1677//__________________________________________________
1678Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const
1679{
1680 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1681 //Returns if the track has a match
1682 if(FindMatchedPosForTrack(trkIndex) < 999)
1683 return kTRUE;
1684 else
1685 return kFALSE;
1686}
1687
1688//__________________________________________________________
1689UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
1690{
1691 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1692 //Returns the position of the match in the fMatchedClusterIndex array
1693 Float_t tmpR = fCutR;
1694 UInt_t pos = 999;
1695
1696 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
1697 {
1698 if(fMatchedClusterIndex->At(i)==clsIndex)
1699 {
1700 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1701 if(r<tmpR)
1702 {
1703 pos=i;
1704 tmpR=r;
1705 AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1706 }
1707 }
1708 }
1709 return pos;
1710}
1711
1712//__________________________________________________________
1713UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
1714{
1715 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1716 //Returns the position of the match in the fMatchedTrackIndex array
1717 Float_t tmpR = fCutR;
1718 UInt_t pos = 999;
1719
1720 for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
1721 {
1722 if(fMatchedTrackIndex->At(i)==trkIndex)
1723 {
1724 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1725 if(r<tmpR)
1726 {
1727 pos=i;
1728 tmpR=r;
1729 AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1730 }
1731 }
1732 }
1733 return pos;
1734}
1735
1736//__________________________________________________________
1737Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells)
1738{
1739 // check if the cluster survives some quality cut
1740 //
1741 //
1742 Bool_t isGood=kTRUE;
1743 if(!cluster || !cluster->IsEMCAL()) return kFALSE;
1744 if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
1745 if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
1746 if(fRejectExoticCluster && IsExoticCluster(cluster)) return kFALSE;
1747
1748 return isGood;
1749}
1750
1751//__________________________________________________________
1752Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1753{
1754 // Given a esd track, return whether the track survive all the cuts
1755
1756 // The different quality parameter are first
1757 // retrieved from the track. then it is found out what cuts the
1758 // track did not survive and finally the cuts are imposed.
1759
1760 UInt_t status = esdTrack->GetStatus();
1761
1762 Int_t nClustersITS = esdTrack->GetITSclusters(0);
1763 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1764
1765 Float_t chi2PerClusterITS = -1;
1766 Float_t chi2PerClusterTPC = -1;
1767 if (nClustersITS!=0)
1768 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1769 if (nClustersTPC!=0)
1770 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1771
1772
1773 //DCA cuts
1774 if(fTrackCutsType==kGlobalCut)
1775 {
1776 Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1777 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1778 SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1779 }
1780
1781
1782 Float_t b[2];
1783 Float_t bCov[3];
1784 esdTrack->GetImpactParameters(b,bCov);
1785 if (bCov[0]<=0 || bCov[2]<=0) {
1786 AliDebug(1, "Estimated b resolution lower or equal zero!");
1787 bCov[0]=0; bCov[2]=0;
1788 }
1789
1790 Float_t dcaToVertexXY = b[0];
1791 Float_t dcaToVertexZ = b[1];
1792 Float_t dcaToVertex = -1;
1793
1794 if (fCutDCAToVertex2D)
1795 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1796 else
1797 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1798
1799 // cut the track?
1800
1801 Bool_t cuts[kNCuts];
1802 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1803
1804 // track quality cuts
1805 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1806 cuts[0]=kTRUE;
1807 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1808 cuts[1]=kTRUE;
1809 if (nClustersTPC<fCutMinNClusterTPC)
1810 cuts[2]=kTRUE;
1811 if (nClustersITS<fCutMinNClusterITS)
1812 cuts[3]=kTRUE;
1813 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1814 cuts[4]=kTRUE;
1815 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1816 cuts[5]=kTRUE;
1817 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1818 cuts[6]=kTRUE;
1819 if (fCutDCAToVertex2D && dcaToVertex > 1)
1820 cuts[7] = kTRUE;
1821 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1822 cuts[8] = kTRUE;
1823 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1824 cuts[9] = kTRUE;
1825
1826 if(fTrackCutsType==kGlobalCut)
1827 {
1828 //Require at least one SPD point + anything else in ITS
1829 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1830 cuts[10] = kTRUE;
1831 }
1832
1833 Bool_t cut=kFALSE;
1834 for (Int_t i=0; i<kNCuts; i++)
1835 if (cuts[i]) {cut = kTRUE;}
1836
1837 // cut the track
1838 if (cut)
1839 return kFALSE;
1840 else
1841 return kTRUE;
1842}
1843
1844
1845//__________________________________________________
1846void AliEMCALRecoUtils::InitTrackCuts()
1847{
1848 //Intilize the track cut criteria
1849 //By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
1850 //Also you can customize the cuts using the setters
1851
1852 switch (fTrackCutsType)
1853 {
1854 case kTPCOnlyCut:
1855 {
1856 AliInfo(Form("Track cuts for matching: GetStandardTPCOnlyTrackCuts()"));
1857 //TPC
1858 SetMinNClustersTPC(70);
1859 SetMaxChi2PerClusterTPC(4);
1860 SetAcceptKinkDaughters(kFALSE);
1861 SetRequireTPCRefit(kFALSE);
1862
1863 //ITS
1864 SetRequireITSRefit(kFALSE);
1865 SetMaxDCAToVertexZ(3.2);
1866 SetMaxDCAToVertexXY(2.4);
1867 SetDCAToVertex2D(kTRUE);
1868
1869 break;
1870 }
1871
1872 case kGlobalCut:
1873 {
1874 AliInfo(Form("Track cuts for matching: GetStandardITSTPCTrackCuts2010(kTURE)"));
1875 //TPC
1876 SetMinNClustersTPC(70);
1877 SetMaxChi2PerClusterTPC(4);
1878 SetAcceptKinkDaughters(kFALSE);
1879 SetRequireTPCRefit(kTRUE);
1880
1881 //ITS
1882 SetRequireITSRefit(kTRUE);
1883 SetMaxDCAToVertexZ(2);
1884 SetMaxDCAToVertexXY();
1885 SetDCAToVertex2D(kFALSE);
1886
1887 break;
1888 }
1889
1890 case kLooseCut:
1891 {
1892 AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
1893 SetMinNClustersTPC(50);
1894 SetAcceptKinkDaughters(kTRUE);
1895
1896 break;
1897 }
1898 }
1899}
1900
1901//___________________________________________________
1902void AliEMCALRecoUtils::Print(const Option_t *) const
1903{
1904 // Print Parameters
1905
1906 printf("AliEMCALRecoUtils Settings: \n");
1907 printf("Misalignment shifts\n");
1908 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,
1909 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1910 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
1911 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1912 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
1913
1914 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
1915
1916 printf("Matching criteria: ");
1917 if(fCutEtaPhiSum)
1918 {
1919 printf("sqrt(dEta^2+dPhi^2)<%2.2f\n",fCutR);
1920 }
1921 else if(fCutEtaPhiSeparate)
1922 {
1923 printf("dEta<%2.2f, dPhi<%2.2f\n",fCutEta,fCutPhi);
1924 }
1925 else
1926 {
1927 printf("Error\n");
1928 printf("please specify your cut criteria\n");
1929 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1930 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1931 }
1932
1933 printf("Mass hypothesis = %2.3f [GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
1934
1935 printf("Track cuts: \n");
1936 printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
1937 printf("AOD track selection mask: %d\n",fAODFilterMask);
1938 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1939 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1940 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1941 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1942 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1943
1944}
1945
1946//_____________________________________________________________________
1947void AliEMCALRecoUtils::SetRunDependentCorrections(Int_t runnumber){
1948 //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1949 //Do it only once and only if it is requested
1950
1951 if(!fUseRunCorrectionFactors) return;
1952 if(fRunCorrectionFactorsSet) return;
1953
1954 AliInfo(Form("AliEMCALRecoUtils::GetRunDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber));
1955
1956 AliEMCALCalibTimeDepCorrection *corr = new AliEMCALCalibTimeDepCorrection();
1957 corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1958
1959 SwitchOnRecalibration();
1960 for(Int_t ism = 0; ism < 4; ism++){
1961 for(Int_t icol = 0; icol < 48; icol++){
1962 for(Int_t irow = 0; irow < 24; irow++){
1963 Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1964 Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1965 GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1966 //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
1967 // orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1968 // (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1969 }
1970 }
1971 }
1972 fRunCorrectionFactorsSet = kTRUE;
1973}
1974