monte carlo flag
[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)
25///////////////////////////////////////////////////////////////////////////////
26
27// --- standard c ---
28
29// standard C++ includes
30//#include <Riostream.h>
31
32// ROOT includes
094786cc 33#include <TGeoManager.h>
34#include <TGeoMatrix.h>
35#include <TGeoBBox.h>
d9b3567c 36
37// STEER includes
38#include "AliEMCALRecoUtils.h"
094786cc 39#include "AliEMCALGeometry.h"
d9b3567c 40#include "AliVCluster.h"
41#include "AliVCaloCells.h"
bd8c7aef 42#include "AliVEvent.h"
d9b3567c 43#include "AliLog.h"
83bfd77a 44#include "AliEMCALPIDUtils.h"
45#include "AliPID.h"
bd8c7aef 46#include "AliESDEvent.h"
47#include "AliESDtrack.h"
48#include "AliEMCALTrack.h"
d9b3567c 49
50ClassImp(AliEMCALRecoUtils)
51
52//______________________________________________
53AliEMCALRecoUtils::AliEMCALRecoUtils():
fd6df01c 54 fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
55 fPosAlgo(kUnchanged),fW0(4.),
56 fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
bd8c7aef 57 fRemoveBadChannels(kFALSE), fEMCALBadChannelMap(),
58 fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
1e4723b3 59 fMatchedClusterIndex(0x0), fResidualZ(0x0), fResidualR(0x0), fCutR(20), fCutZ(20),
bd8c7aef 60 fCutMinNClusterTPC(0), fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
61 fCutRequireTPCRefit(0), fCutRequireITSRefit(0), fCutAcceptKinkDaughters(0),
62 fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0),fCutDCAToVertex2D(0),
63 fPIDUtils()
d9b3567c 64{
65//
66 // Constructor.
67 // Initialize all constant values which have to be used
68 // during Reco algorithm execution
69 //
70
fd6df01c 71 for(Int_t i = 0; i < 15 ; i++) {
72 fMisalTransShift[i] = 0.;
73 fMisalRotShift[i] = 0.;
74 }
d9b3567c 75 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = 0.;
fd6df01c 76 //For kPi0GammaGamma case, but default is no correction
d9b3567c 77 fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
78 fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
79 fNonLinearityParams[2] = 1.046;
bd8c7aef 80
81 fMatchedClusterIndex = new TArrayI();
82 fResidualZ = new TArrayF();
83 fResidualR = new TArrayF();
d9b3567c 84
83bfd77a 85 fPIDUtils = new AliEMCALPIDUtils();
86
bd8c7aef 87 InitTrackCuts();
88
d9b3567c 89}
90
91//______________________________________________________________________
92AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
094786cc 93: TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction),
fd6df01c 94 fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0),
95 fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
96 fRemoveBadChannels(reco.fRemoveBadChannels),fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
83bfd77a 97 fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
bd8c7aef 98 fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
99 fResidualZ(reco.fResidualZ?new TArrayF(*reco.fResidualZ):0x0),
100 fResidualR(reco.fResidualR?new TArrayF(*reco.fResidualR):0x0),
101 fCutR(reco.fCutR),fCutZ(reco.fCutZ),
102 fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS),
103 fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
104 fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
105 fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters),
106 fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY), fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ),fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
83bfd77a 107 fPIDUtils(reco.fPIDUtils)
fd6df01c 108
d9b3567c 109{
110 //Copy ctor
111
fd6df01c 112 for(Int_t i = 0; i < 15 ; i++) {
113 fMisalRotShift[i] = reco.fMisalRotShift[i];
114 fMisalTransShift[i] = reco.fMisalTransShift[i];
115 }
d9b3567c 116 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
bd8c7aef 117
d9b3567c 118}
119
120
121//______________________________________________________________________
122AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
123{
124 //Assignment operator
125
126 if(this == &reco)return *this;
127 ((TNamed *)this)->operator=(reco);
128
fd6df01c 129 fNonLinearityFunction = reco.fNonLinearityFunction;
130 fParticleType = reco.fParticleType;
131 fPosAlgo = reco.fPosAlgo;
132 fW0 = reco.fW0;
133 fRecalibration = reco.fRecalibration;
094786cc 134 fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
fd6df01c 135 fRemoveBadChannels = reco.fRemoveBadChannels;
136 fEMCALBadChannelMap = reco.fEMCALBadChannelMap;
137 fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
138 fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
bd8c7aef 139
83bfd77a 140
2a71e873 141 for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
d9b3567c 142 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
143
bd8c7aef 144 fCutR = reco.fCutR;
145 fCutZ = reco.fCutZ;
146
147 fCutMinNClusterTPC = reco.fCutMinNClusterTPC;
148 fCutMinNClusterITS = reco.fCutMinNClusterITS;
149 fCutMaxChi2PerClusterTPC = reco.fCutMaxChi2PerClusterTPC;
150 fCutMaxChi2PerClusterITS = reco.fCutMaxChi2PerClusterITS;
151 fCutRequireTPCRefit = reco.fCutRequireTPCRefit;
152 fCutRequireITSRefit = reco.fCutRequireITSRefit;
153 fCutAcceptKinkDaughters = reco.fCutAcceptKinkDaughters;
154 fCutMaxDCAToVertexXY = reco.fCutMaxDCAToVertexXY;
155 fCutMaxDCAToVertexZ = reco.fCutMaxDCAToVertexZ;
156 fCutDCAToVertex2D = reco.fCutDCAToVertex2D;
157
158 fPIDUtils = reco.fPIDUtils;
159
160
161 if(reco.fResidualR){
162 // assign or copy construct
163 if(fResidualR){
164 *fResidualR = *reco.fResidualR;
165 }
166 else fResidualR = new TArrayF(*reco.fResidualR);
167 }
168 else{
169 if(fResidualR)delete fResidualR;
170 fResidualR = 0;
171 }
172
173 if(reco.fResidualZ){
174 // assign or copy construct
175 if(fResidualZ){
176 *fResidualZ = *reco.fResidualZ;
177 }
178 else fResidualZ = new TArrayF(*reco.fResidualZ);
179 }
180 else{
181 if(fResidualZ)delete fResidualZ;
182 fResidualZ = 0;
183 }
184
185
186 if(reco.fMatchedClusterIndex){
187 // assign or copy construct
188 if(fMatchedClusterIndex){
189 *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
190 }
191 else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
192 }
193 else{
194 if(fMatchedClusterIndex)delete fMatchedClusterIndex;
195 fMatchedClusterIndex = 0;
196 }
197
198
d9b3567c 199 return *this;
200}
201
202
203//__________________________________________________
094786cc 204AliEMCALRecoUtils::~AliEMCALRecoUtils()
205{
206 //Destructor.
207
208 if(fEMCALRecalibrationFactors) {
209 fEMCALRecalibrationFactors->Clear();
210 delete fEMCALRecalibrationFactors;
211 }
fd6df01c 212
213 if(fEMCALBadChannelMap) {
214 fEMCALBadChannelMap->Clear();
215 delete fEMCALBadChannelMap;
216 }
bd8c7aef 217
218 if(fMatchedClusterIndex) {delete fMatchedClusterIndex; fMatchedClusterIndex=0;}
219 if(fResidualR) {delete fResidualR; fResidualR=0;}
220 if(fResidualZ) {delete fResidualZ; fResidualZ=0;}
221
094786cc 222}
223
fd6df01c 224//_______________________________________________________________
225Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells)
226{
227 // Given the list of AbsId of the cluster, get the maximum cell and
228 // check if there are fNCellsFromBorder from the calorimeter border
229
230 //If the distance to the border is 0 or negative just exit accept all clusters
231 if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
232
cb231979 233 Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
234 Bool_t shared = kFALSE;
235 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
fd6df01c 236
83bfd77a 237 AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
238 absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
fd6df01c 239
240 if(absIdMax==-1) return kFALSE;
241
242 //Check if the cell is close to the borders:
243 Bool_t okrow = kFALSE;
244 Bool_t okcol = kFALSE;
245
246 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
247 AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
248 iSM,ieta,iphi));
249 }
250
251 //Check rows/phi
252 if(iSM < 10){
253 if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
254 }
255 else{
256 if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
257 }
258
259 //Check columns/eta
260 if(!fNoEMCALBorderAtEta0){
261 if(ieta > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
262 }
263 else{
264 if(iSM%2==0){
265 if(ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
266 }
267 else {
268 if(ieta < 48-fNCellsFromEMCALBorder) okcol = kTRUE;
269 }
270 }//eta 0 not checked
271
83bfd77a 272 AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
fd6df01c 273 fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
274
83bfd77a 275 if (okcol && okrow) {
276 //printf("Accept\n");
277 return kTRUE;
278 }
279 else {
280 //printf("Reject\n");
281 AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
282 return kFALSE;
283 }
fd6df01c 284
285}
286
287
288//_________________________________________________________________________________________________________
289Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells){
290 // Check that in the cluster cells, there is no bad channel of those stored
291 // in fEMCALBadChannelMap or fPHOSBadChannelMap
292
293 if(!fRemoveBadChannels) return kFALSE;
294 if(!fEMCALBadChannelMap) return kFALSE;
295
296 Int_t icol = -1;
297 Int_t irow = -1;
298 Int_t imod = -1;
299 for(Int_t iCell = 0; iCell<nCells; iCell++){
300
301 //Get the column and row
302 Int_t iTower = -1, iIphi = -1, iIeta = -1;
303 geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
304 if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
305 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
83bfd77a 306 if(GetEMCALChannelStatus(imod, icol, irow)){
307 AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
308 return kTRUE;
309 }
fd6df01c 310
311 }// cell cluster loop
312
313 return kFALSE;
314
315}
094786cc 316
317//__________________________________________________
d9b3567c 318Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
319// Correct cluster energy from non linearity functions
320 Float_t energy = cluster->E();
321
322 switch (fNonLinearityFunction) {
323
324 case kPi0MC:
325 //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
326 //Double_t par0 = 1.001;
327 //Double_t par1 = -0.01264;
328 //Double_t par2 = -0.03632;
329 //Double_t par3 = 0.1798;
330 //Double_t par4 = -0.522;
331 energy /= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
332 ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
333 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
334 break;
335
336 case kPi0GammaGamma:
337
338 //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
339 //Double_t par0 = 0.1457;
340 //Double_t par1 = -0.02024;
341 //Double_t par2 = 1.046;
342 energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
343 break;
344
345 case kPi0GammaConversion:
346
347 //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
348 //Double_t C = 0.139393/0.1349766;
349 //Double_t a = 0.0566186;
350 //Double_t b = 0.982133;
351 energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
352
353 break;
354
355 case kNoCorrection:
356 AliDebug(2,"No correction on the energy\n");
357 break;
358
359 }
360
361 return energy;
362
363}
364
365//__________________________________________________
094786cc 366Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
367{
368 //Calculate shower depth for a given cluster energy and particle type
369
370 // parameters
cb231979 371 Float_t x0 = 1.31;
094786cc 372 Float_t ecr = 8;
373 Float_t depth = 0;
374
375 switch ( iParticle )
376 {
377 case kPhoton:
fd6df01c 378 depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 379 break;
380
381 case kElectron:
fd6df01c 382 depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 383 break;
384
385 case kHadron:
386 // hadron
387 // boxes anc. here
388 if(gGeoManager){
389 gGeoManager->cd("ALIC_1/XEN1_1");
390 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
391 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
fd6df01c 392 if(geoSM){
393 TGeoVolume *geoSMVol = geoSM->GetVolume();
394 TGeoShape *geoSMShape = geoSMVol->GetShape();
395 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
396 if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
397 else AliFatal("Null GEANT box");
398 }else AliFatal("NULL GEANT node matrix");
094786cc 399 }
400 else{//electron
fd6df01c 401 depth = x0 * (TMath::Log(energy*1000 / ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 402 }
403
404 break;
405
406 default://photon
fd6df01c 407 depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
094786cc 408 }
409
410 return depth;
411
412}
413
414//__________________________________________________
cb231979 415void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
416 Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
d9b3567c 417{
418 //For a given CaloCluster gets the absId of the cell
419 //with maximum energy deposit.
420
421 Double_t eMax = -1.;
422 Double_t eCell = -1.;
094786cc 423 Float_t fraction = 1.;
424 Float_t recalFactor = 1.;
d9b3567c 425 Int_t cellAbsId = -1 ;
094786cc 426
d9b3567c 427 Int_t iTower = -1;
428 Int_t iIphi = -1;
429 Int_t iIeta = -1;
cb231979 430 Int_t iSupMod0= -1;
83bfd77a 431 //printf("---Max?\n");
d9b3567c 432 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 433 cellAbsId = clu->GetCellAbsId(iDig);
434 fraction = clu->GetCellAmplitudeFraction(iDig);
83bfd77a 435 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
094786cc 436 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
cb231979 437 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
438 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
439 if(iDig==0) iSupMod0=iSupMod;
440 else if(iSupMod0!=iSupMod) {
441 shared = kTRUE;
442 //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
443 }
094786cc 444 if(IsRecalibrationOn()) {
094786cc 445 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
446 }
447 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
83bfd77a 448 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
094786cc 449 if(eCell > eMax) {
d9b3567c 450 eMax = eCell;
451 absId = cellAbsId;
452 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
453 }
454 }// cell loop
455
456 //Get from the absid the supermodule, tower and eta/phi numbers
457 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
458 //Gives SuperModule and Tower numbers
459 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
83bfd77a 460 iIphi, iIeta,iphi,ieta);
461 //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
462 //printf("Max end---\n");
d9b3567c 463
464}
465
094786cc 466//________________________________________________________________
467void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
468 //Init EMCAL recalibration factors
469 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
470 //In order to avoid rewriting the same histograms
471 Bool_t oldStatus = TH1::AddDirectoryStatus();
472 TH1::AddDirectory(kFALSE);
473
cb231979 474 fEMCALRecalibrationFactors = new TObjArray(10);
094786cc 475 for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
476 //Init the histograms with 1
477 for (Int_t sm = 0; sm < 12; sm++) {
478 for (Int_t i = 0; i < 48; i++) {
479 for (Int_t j = 0; j < 24; j++) {
480 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
481 }
482 }
483 }
484 fEMCALRecalibrationFactors->SetOwner(kTRUE);
485 fEMCALRecalibrationFactors->Compress();
486
487 //In order to avoid rewriting the same histograms
488 TH1::AddDirectory(oldStatus);
489}
490
491
492//________________________________________________________________
fd6df01c 493void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
494 //Init EMCAL bad channels map
495 AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
496 //In order to avoid rewriting the same histograms
497 Bool_t oldStatus = TH1::AddDirectoryStatus();
498 TH1::AddDirectory(kFALSE);
499
cb231979 500 fEMCALBadChannelMap = new TObjArray(10);
fd6df01c 501 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
502 for (int i = 0; i < 12; i++) {
503 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
504 }
505
506 //delete hTemp;
507
508 fEMCALBadChannelMap->SetOwner(kTRUE);
509 fEMCALBadChannelMap->Compress();
510
511 //In order to avoid rewriting the same histograms
512 TH1::AddDirectory(oldStatus);
513}
514
515//________________________________________________________________
094786cc 516void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
517 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
518
519 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
520 UShort_t * index = cluster->GetCellsAbsId() ;
521 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
522 Int_t ncells = cluster->GetNCells();
523
524 //Initialize some used variables
525 Float_t energy = 0;
526 Int_t absId = -1;
527 Int_t icol = -1, irow = -1, imod=1;
528 Float_t factor = 1, frac = 0;
529
530 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
531 for(Int_t icell = 0; icell < ncells; icell++){
532 absId = index[icell];
533 frac = fraction[icell];
534 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
535 Int_t iTower = -1, iIphi = -1, iIeta = -1;
536 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
537 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
538 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
539 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
540 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
541 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
542
543 energy += cells->GetCellAmplitude(absId)*factor*frac;
544 }
545
546
547 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
548
549 cluster->SetE(energy);
550
551}
552
553
d9b3567c 554//__________________________________________________
094786cc 555void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
d9b3567c 556{
557 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
558
094786cc 559 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
560 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
fd6df01c 561 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
094786cc 562
563}
564
565//__________________________________________________
566void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
567{
568 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
569 // The algorithm is a copy of what is done in AliEMCALRecPoint
570
571 Double_t eCell = 0.;
572 Float_t fraction = 1.;
573 Float_t recalFactor = 1.;
574
575 Int_t absId = -1;
576 Int_t iTower = -1, iIphi = -1, iIeta = -1;
577 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
578 Float_t weight = 0., totalWeight=0.;
579 Float_t newPos[3] = {0,0,0};
580 Double_t pLocal[3], pGlobal[3];
cb231979 581 Bool_t shared = kFALSE;
582
094786cc 583 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
cb231979 584 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 585 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
586
83bfd77a 587 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
588
094786cc 589 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
590 absId = clu->GetCellAbsId(iDig);
591 fraction = clu->GetCellAmplitudeFraction(iDig);
592 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
593 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
594 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
595
596 if(IsRecalibrationOn()) {
597 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
598 }
599 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
600
601 weight = GetCellWeight(eCell,clEnergy);
83bfd77a 602 //printf("cell energy %f, weight %f\n",eCell,weight);
094786cc 603 totalWeight += weight;
604 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
83bfd77a 605 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
094786cc 606 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
83bfd77a 607 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
608
094786cc 609 for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
610
611 }// cell loop
612
613 if(totalWeight>0){
614 for(int i=0; i<3; i++ ) newPos[i] /= totalWeight;
615 }
616
094786cc 617 //Float_t pos[]={0,0,0};
618 //clu->GetPosition(pos);
619 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
83bfd77a 620 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
094786cc 621
622 if(iSupModMax > 1) {//sector 1
623 newPos[0] +=fMisalTransShift[3];//-=3.093;
624 newPos[1] +=fMisalTransShift[4];//+=6.82;
625 newPos[2] +=fMisalTransShift[5];//+=1.635;
83bfd77a 626 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
627
094786cc 628 }
629 else {//sector 0
630 newPos[0] +=fMisalTransShift[0];//+=1.134;
631 newPos[1] +=fMisalTransShift[1];//+=8.2;
632 newPos[2] +=fMisalTransShift[2];//+=1.197;
83bfd77a 633 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
634
094786cc 635 }
83bfd77a 636 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
637
094786cc 638 clu->SetPosition(newPos);
639
094786cc 640}
641
642//__________________________________________________
643void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
644{
645 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
646 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
647
648 Double_t eCell = 1.;
649 Float_t fraction = 1.;
650 Float_t recalFactor = 1.;
651
652 Int_t absId = -1;
d9b3567c 653 Int_t iTower = -1;
094786cc 654 Int_t iIphi = -1, iIeta = -1;
655 Int_t iSupMod = -1, iSupModMax = -1;
d9b3567c 656 Int_t iphi = -1, ieta =-1;
cb231979 657 Bool_t shared = kFALSE;
658
d9b3567c 659 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
cb231979 660 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
094786cc 661 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
662
d9b3567c 663 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
094786cc 664 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
665 Int_t startingSM = -1;
d9b3567c 666
667 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
094786cc 668 absId = clu->GetCellAbsId(iDig);
669 fraction = clu->GetCellAmplitudeFraction(iDig);
670 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
671 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
d9b3567c 672 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
094786cc 673
d9b3567c 674 if (iDig==0) startingSM = iSupMod;
675 else if(iSupMod != startingSM) areInSameSM = kFALSE;
094786cc 676
677 eCell = cells->GetCellAmplitude(absId);
d9b3567c 678
094786cc 679 if(IsRecalibrationOn()) {
680 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
681 }
682 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
d9b3567c 683
094786cc 684 weight = GetCellWeight(eCell,clEnergy);
d9b3567c 685 if(weight < 0) weight = 0;
686 totalWeight += weight;
687 weightedCol += ieta*weight;
688 weightedRow += iphi*weight;
689
690 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
691
094786cc 692 }// cell loop
693
d9b3567c 694 Float_t xyzNew[]={0.,0.,0.};
695 if(areInSameSM == kTRUE) {
696 //printf("In Same SM\n");
697 weightedCol = weightedCol/totalWeight;
698 weightedRow = weightedRow/totalWeight;
094786cc 699 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 700 }
701 else {
702 //printf("In Different SM\n");
094786cc 703 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
d9b3567c 704 }
d9b3567c 705
094786cc 706 clu->SetPosition(xyzNew);
d9b3567c 707
708}
709
83bfd77a 710//____________________________________________________________________________
cb231979 711void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
712
713 //re-evaluate distance to bad channel with updated bad map
714
715 if(!fRemoveBadChannels) return;
716
717 //Get channels map of the supermodule where the cluster is.
718
719 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
720 Bool_t shared = kFALSE;
721 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
722 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
723
724 Int_t dRrow, dRcol;
725 Float_t minDist = 10000.;
726 Float_t dist = 0.;
727
728 //Loop on tower status map
729 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
730 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
731 //Check if tower is bad.
732 if(hMap->GetBinContent(icol,irow)==0) continue;
733 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
734 // iSupMod,icol, irow, icolM,irowM);
735
736 dRrow=TMath::Abs(irowM-irow);
737 dRcol=TMath::Abs(icolM-icol);
738 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
739 if(dist < minDist){
740 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
741 minDist = dist;
742 }
743
744 }
745 }
746
747 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
748 if (shared) {
749 TH2D* hMap2 = 0;
750 Int_t iSupMod2 = -1;
751
752 //The only possible combinations are (0,1), (2,3) ... (8,9)
753 if(iSupMod%2) iSupMod2 = iSupMod-1;
754 else iSupMod2 = iSupMod+1;
755 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
756
757 //Loop on tower status map of second super module
758 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
759 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
760 //Check if tower is bad.
761 if(hMap2->GetBinContent(icol,irow)==0) continue;
762 //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",
763 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
764
765 dRrow=TMath::Abs(irow-irowM);
766
767 if(iSupMod%2) {
768 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
769 }
770 else {
771 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
772 }
773
774 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
775 if(dist < minDist) minDist = dist;
776
777 }
778 }
779
780 }// shared cluster in 2 SuperModules
781 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - Distance to Bad Channel %2.2f\n",minDist);
782 cluster->SetDistanceToBadChannel(minDist);
783
784}
785
786//____________________________________________________________________________
83bfd77a 787void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
788
789 //re-evaluate identification parameters with bayesian
790
791 if ( cluster->GetM02() != 0)
792 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
793
794 Float_t pidlist[AliPID::kSPECIESN+1];
795 for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
796
797 cluster->SetPID(pidlist);
798
799}
800
801//____________________________________________________________________________
802void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
803{
804 // Calculates new center of gravity in the local EMCAL-module coordinates
805 // and tranfers into global ALICE coordinates
806 // Calculates Dispersion and main axis
807
808 Int_t nstat = 0;
809 Float_t wtot = 0. ;
810 Double_t eCell = 0.;
811 Float_t fraction = 1.;
812 Float_t recalFactor = 1.;
813
814 Int_t iSupMod = -1;
815 Int_t iTower = -1;
816 Int_t iIphi = -1;
817 Int_t iIeta = -1;
818 Int_t iphi = -1;
819 Int_t ieta = -1;
820 Double_t etai = -1.;
821 Double_t phii = -1.;
822
823 Double_t w = 0.;
824 Double_t d = 0.;
825 Double_t dxx = 0.;
826 Double_t dzz = 0.;
827 Double_t dxz = 0.;
828 Double_t xmean = 0.;
829 Double_t zmean = 0.;
830
831 //Loop on cells
832 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
833
834 //Get from the absid the supermodule, tower and eta/phi numbers
835 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
836 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
837
838 //Get the cell energy, if recalibration is on, apply factors
839 fraction = cluster->GetCellAmplitudeFraction(iDigit);
840 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
841 if(IsRecalibrationOn()) {
842 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
843 }
844 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
845
846 if(cluster->E() > 0 && eCell > 0){
847
848 w = GetCellWeight(eCell,cluster->E());
849
850 etai=(Double_t)ieta;
851 phii=(Double_t)iphi;
852 if(w > 0.0) {
853 wtot += w ;
854 nstat++;
855 //Shower shape
856 dxx += w * etai * etai ;
857 xmean+= w * etai ;
858 dzz += w * phii * phii ;
859 zmean+= w * phii ;
860 dxz += w * etai * phii ;
861 }
862 }
863 else
864 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
865 }//cell loop
866
867 //Normalize to the weight
868 if (wtot > 0) {
869 xmean /= wtot ;
870 zmean /= wtot ;
871 }
872 else
873 AliError(Form("Wrong weight %f\n", wtot));
874
875 //Calculate dispersion
876 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
877
878 //Get from the absid the supermodule, tower and eta/phi numbers
879 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
880 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
881
882 //Get the cell energy, if recalibration is on, apply factors
883 fraction = cluster->GetCellAmplitudeFraction(iDigit);
884 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
885 if(IsRecalibrationOn()) {
886 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
887 }
888 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
889
890 if(cluster->E() > 0 && eCell > 0){
891
892 w = GetCellWeight(eCell,cluster->E());
893
894 etai=(Double_t)ieta;
895 phii=(Double_t)iphi;
896 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
897 }
898 else
899 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
900 }// cell loop
901
902 //Normalize to the weigth and set shower shape parameters
903 if (wtot > 0 && nstat > 1) {
904 d /= wtot ;
905 dxx /= wtot ;
906 dzz /= wtot ;
907 dxz /= wtot ;
908 dxx -= xmean * xmean ;
909 dzz -= zmean * zmean ;
910 dxz -= xmean * zmean ;
911 cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
912 cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
913 }
914 else{
915 d=0. ;
916 cluster->SetM20(0.) ;
917 cluster->SetM02(0.) ;
918 }
919
920 if (d>=0)
921 cluster->SetDispersion(TMath::Sqrt(d)) ;
922 else
923 cluster->SetDispersion(0) ;
924
925}
926
bd8c7aef 927//__________________________________________________
928void AliEMCALRecoUtils::FindMatches(AliVEvent *event)
929{
930 //This function should be called before the cluster loop
931 //Before call this function, please recalculate the cluster positions
932 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
933 //Store matched cluster indexes and residuals
934 //It only works with ESDs, not AODs
935
936 fMatchedClusterIndex->Reset();
937 fResidualZ->Reset();
938 fResidualR->Reset();
939
940 fMatchedClusterIndex->Set(100);
941 fResidualZ->Set(100);
942 fResidualR->Set(100);
943
944 Int_t matched=0;
945 Float_t clsPos[3];
946 Double_t trkPos[3];
947 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
948 {
949 AliESDtrack *track = ((AliESDEvent*)event)->GetTrack(itr);
950 if(!track || !IsAccepted(track)) continue;
951
952 Float_t dRMax = fCutR, dZMax = fCutZ;
953 Int_t index = -1;
954 AliEMCALTrack *emctrack = new AliEMCALTrack(*track);
955 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
956 {
957 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
958 if(!cluster->IsEMCAL()) continue;
959 cluster->GetPosition(clsPos); //Has been recalculated
960 if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) continue;
961 emctrack->GetXYZ(trkPos);
962 Float_t tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
963 Float_t tmpZ = TMath::Abs(clsPos[2]-trkPos[2]);
964
965 if(tmpR<dRMax)
966 {
967 dRMax=tmpR;
968 dZMax=tmpZ;
969 index=icl;
970 }
971
972 }//cluser loop
973
974 if(index>-1)
975 {
976 fMatchedClusterIndex->AddAt(index,matched);
977 fResidualZ->AddAt(dZMax,matched);
978 fResidualR->AddAt(dRMax,matched);
979 matched++;
980 }
981 delete emctrack;
982 }//track loop
983 fMatchedClusterIndex->Set(matched);
984 fResidualZ->Set(matched);
985 fResidualR->Set(matched);
986
987 //printf("Number of matched pairs: %d\n",matched);
988}
989
990//__________________________________________________
991void AliEMCALRecoUtils::GetMatchedResiduals(Int_t index, Float_t &dR, Float_t &dZ)
992{
993 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
994 //Get the residuals dR and dZ for this cluster
995 //It only works with ESDs, not AODs
996
81efb149 997 if( FindMatchedPos(index) >= 999 )
bd8c7aef 998 {
999 AliDebug(2,"No matched tracks found!\n");
1000 dR=999.;
1001 dZ=999.;
1002 return;
1003 }
1004 dR = fResidualR->At(FindMatchedPos(index));
1005 dZ = fResidualZ->At(FindMatchedPos(index));
1006}
1007
1008//__________________________________________________
1009Bool_t AliEMCALRecoUtils::IsMatched(Int_t index)
1010{
1011 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1012 //Returns if cluster has a match
81efb149 1013 if(FindMatchedPos(index) < 999)
82d09e74 1014 return kTRUE;
bd8c7aef 1015 else
1016 return kFALSE;
1017}
1018//__________________________________________________
81efb149 1019UInt_t AliEMCALRecoUtils::FindMatchedPos(Int_t index) const
bd8c7aef 1020{
1021 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
1022 //Returns the position of the match in the fMatchedClusterIndex array
1023 Float_t tmpR = fCutR;
81efb149 1024 UInt_t pos = 999;
bd8c7aef 1025
1026 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
1027 {
1028 if(fMatchedClusterIndex->At(i)==index && fResidualR->At(i)<tmpR)
1029 {
1030 pos=i;
1031 tmpR=fResidualR->At(i);
1032 }
82d09e74 1033 AliDebug(3,Form("Matched cluster pos: %d, index: %d, dR: %2.4f, dZ: %2.4f.\n",i,fMatchedClusterIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
bd8c7aef 1034 }
1035 return pos;
1036}
1037
1038Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1039{
1040 // Given a esd track, return whether the track survive all the cuts
1041
1042 // The different quality parameter are first
1043 // retrieved from the track. then it is found out what cuts the
1044 // track did not survive and finally the cuts are imposed.
1045
1046 UInt_t status = esdTrack->GetStatus();
1047
1048 Int_t nClustersITS = esdTrack->GetITSclusters(0);
1049 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1050
1051 Float_t chi2PerClusterITS = -1;
1052 Float_t chi2PerClusterTPC = -1;
1053 if (nClustersITS!=0)
1054 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1055 if (nClustersTPC!=0)
1056 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
82d09e74 1057
1058
1059 //DCA cuts
1060 Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1061 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1062 SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1063
1064
bd8c7aef 1065 Float_t b[2];
1066 Float_t bCov[3];
1067 esdTrack->GetImpactParameters(b,bCov);
1068 if (bCov[0]<=0 || bCov[2]<=0) {
1069 AliDebug(1, "Estimated b resolution lower or equal zero!");
1070 bCov[0]=0; bCov[2]=0;
1071 }
1072
1073 Float_t dcaToVertexXY = b[0];
1074 Float_t dcaToVertexZ = b[1];
1075 Float_t dcaToVertex = -1;
1076
1077 if (fCutDCAToVertex2D)
1078 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1079 else
1080 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1081
1082 // cut the track?
1083
1084 Bool_t cuts[kNCuts];
1085 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1086
1087 // track quality cuts
1088 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1089 cuts[0]=kTRUE;
1090 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1091 cuts[1]=kTRUE;
1092 if (nClustersTPC<fCutMinNClusterTPC)
1093 cuts[2]=kTRUE;
1094 if (nClustersITS<fCutMinNClusterITS)
1095 cuts[3]=kTRUE;
1096 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1097 cuts[4]=kTRUE;
1098 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1099 cuts[5]=kTRUE;
1100 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1101 cuts[6]=kTRUE;
1102 if (fCutDCAToVertex2D && dcaToVertex > 1)
1103 cuts[7] = kTRUE;
1104 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1105 cuts[8] = kTRUE;
1106 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1107 cuts[9] = kTRUE;
1108
82d09e74 1109 //Require at least one SPD point + anything else in ITS
1110 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1111 cuts[10] = kTRUE;
1112
bd8c7aef 1113 Bool_t cut=kFALSE;
1114 for (Int_t i=0; i<kNCuts; i++)
1115 if (cuts[i]) {cut = kTRUE;}
1116
1117 // cut the track
1118 if (cut)
1119 return kFALSE;
1120 else
1121 return kTRUE;
1122}
1123//__________________________________________________
1124void AliEMCALRecoUtils::InitTrackCuts()
1125{
1126 //Intilize the track cut criteria
82d09e74 1127 //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
bd8c7aef 1128 //Also you can customize the cuts using the setters
82d09e74 1129
1130 //TPC
1131 SetMinNClustersTPC(70);
bd8c7aef 1132 SetMaxChi2PerClusterTPC(4);
1133 SetAcceptKinkDaughters(kFALSE);
82d09e74 1134 SetRequireTPCRefit(kTRUE);
1135
1136 //ITS
1137 SetRequireITSRefit(kTRUE);
1138 SetMaxDCAToVertexZ(2);
1139 SetDCAToVertex2D(kFALSE);
bd8c7aef 1140}
83bfd77a 1141
d9b3567c 1142//__________________________________________________
1143void AliEMCALRecoUtils::Print(const Option_t *) const
1144{
1145 // Print Parameters
1146
1147 printf("AliEMCALRecoUtils Settings: \n");
1148 printf("Misalignment shifts\n");
2a71e873 1149 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,
1150 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1151 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
d9b3567c 1152 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1153 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
094786cc 1154
1155 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
bd8c7aef 1156
1157 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1158
1159 printf("Track cuts: \n");
1160 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1161 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1162 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1163 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1164 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1165
d9b3567c 1166
1167}