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