1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* $Id: AliEMCALRecoUtils.cxx 33808 2009-07-15 09:48:08Z gconesab $ */
18 ///////////////////////////////////////////////////////////////////////////////
20 // Class AliEMCALRecoUtils
21 // Some utilities to recalculate the cluster position or energy linearity
24 // Author: Gustavo Conesa (LPSC- Grenoble)
25 // Track matching part: Rongrong Ma (Yale)
27 ///////////////////////////////////////////////////////////////////////////////
30 // standard C++ includes
31 //#include <Riostream.h>
34 #include <TGeoManager.h>
35 #include <TGeoMatrix.h>
39 #include "AliVCluster.h"
40 #include "AliVCaloCells.h"
41 #include "AliVEvent.h"
44 #include "AliESDEvent.h"
45 #include "AliAODEvent.h"
46 #include "AliESDtrack.h"
47 #include "AliAODTrack.h"
48 #include "AliExternalTrackParam.h"
49 #include "AliESDfriendTrack.h"
50 #include "AliTrackerBase.h"
53 #include "AliEMCALRecoUtils.h"
54 #include "AliEMCALGeometry.h"
55 #include "AliEMCALTrack.h"
56 #include "AliEMCALCalibTimeDepCorrection.h"
57 #include "AliEMCALPIDUtils.h"
59 ClassImp(AliEMCALRecoUtils)
61 //______________________________________________
62 AliEMCALRecoUtils::AliEMCALRecoUtils():
63 fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
64 fPosAlgo(kUnchanged),fW0(4.), fNonLinearThreshold(30),
65 fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
66 fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
67 fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
69 fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
70 fResidualEta(0x0), fResidualPhi(0x0), fCutEtaPhiSum(kTRUE), fCutEtaPhiSeparate(kFALSE), fCutR(0.1), fCutEta(0.02), fCutPhi(0.04), fMass(0.139), fStep(1),
71 fRejectExoticCluster(kFALSE),
72 fCutMinTrackPt(0), fCutMinNClusterTPC(0), fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
73 fCutRequireTPCRefit(0), fCutRequireITSRefit(0), fCutAcceptKinkDaughters(0),
74 fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0),fCutDCAToVertex2D(0),fPIDUtils(),
75 fUseTimeCorrectionFactors(kFALSE), fTimeCorrectionFactorsSet(kFALSE)
79 // Initialize all constant values which have to be used
80 // during Reco algorithm execution
83 //Misalignment matrices
84 for(Int_t i = 0; i < 15 ; i++) {
85 fMisalTransShift[i] = 0.;
86 fMisalRotShift[i] = 0.;
90 for(Int_t i = 0; i < 7 ; i++) fNonLinearityParams[i] = 0.;
92 //For kBeamTestCorrected case, but default is no correction
93 fNonLinearityParams[0] = 0.99078;
94 fNonLinearityParams[1] = 0.161499;
95 fNonLinearityParams[2] = 0.655166;
96 fNonLinearityParams[3] = 0.134101;
97 fNonLinearityParams[4] = 163.282;
98 fNonLinearityParams[5] = 23.6904;
99 fNonLinearityParams[6] = 0.978;
101 //For kPi0GammaGamma case
102 //fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
103 //fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
104 //fNonLinearityParams[2] = 1.046;
107 fMatchedTrackIndex = new TArrayI();
108 fMatchedClusterIndex = new TArrayI();
109 fResidualPhi = new TArrayF();
110 fResidualEta = new TArrayF();
114 fPIDUtils = new AliEMCALPIDUtils();
119 //______________________________________________________________________
120 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
121 : TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction),
122 fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0), fNonLinearThreshold(reco.fNonLinearThreshold),
123 fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
124 fRemoveBadChannels(reco.fRemoveBadChannels),fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
125 fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
126 fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
127 fAODFilterMask(reco.fAODFilterMask),
128 fMatchedTrackIndex(reco.fMatchedTrackIndex?new TArrayI(*reco.fMatchedTrackIndex):0x0),
129 fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
130 fResidualEta(reco.fResidualEta?new TArrayF(*reco.fResidualEta):0x0),
131 fResidualPhi(reco.fResidualPhi?new TArrayF(*reco.fResidualPhi):0x0),
132 fCutEtaPhiSum(reco.fCutEtaPhiSum), fCutEtaPhiSeparate(reco.fCutEtaPhiSeparate), fCutR(reco.fCutR), fCutEta(reco.fCutEta), fCutPhi(reco.fCutPhi),
133 fMass(reco.fMass), fStep(reco.fStep),
134 fRejectExoticCluster(reco.fRejectExoticCluster),
135 fCutMinTrackPt(reco.fCutMinTrackPt), fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS),
136 fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
137 fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
138 fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters),
139 fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY), fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ),fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
140 fPIDUtils(reco.fPIDUtils),
141 fUseTimeCorrectionFactors(reco.fUseTimeCorrectionFactors), fTimeCorrectionFactorsSet(reco.fTimeCorrectionFactorsSet)
145 for(Int_t i = 0; i < 15 ; i++) {
146 fMisalRotShift[i] = reco.fMisalRotShift[i];
147 fMisalTransShift[i] = reco.fMisalTransShift[i];
149 for(Int_t i = 0; i < 7 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
154 //______________________________________________________________________
155 AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
157 //Assignment operator
159 if(this == &reco)return *this;
160 ((TNamed *)this)->operator=(reco);
162 fNonLinearityFunction = reco.fNonLinearityFunction;
163 fParticleType = reco.fParticleType;
164 fPosAlgo = reco.fPosAlgo;
166 fNonLinearThreshold = reco.fNonLinearThreshold;
167 fRecalibration = reco.fRecalibration;
168 fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
169 fRemoveBadChannels = reco.fRemoveBadChannels;
170 fRecalDistToBadChannels = reco.fRecalDistToBadChannels;
171 fEMCALBadChannelMap = reco.fEMCALBadChannelMap;
172 fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
173 fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
176 for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
177 for(Int_t i = 0; i < 7 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
179 fAODFilterMask = reco.fAODFilterMask;
181 fCutEtaPhiSum = reco.fCutEtaPhiSum;
182 fCutEtaPhiSeparate = reco.fCutEtaPhiSeparate;
184 fCutEta = reco.fCutEta;
185 fCutPhi = reco.fCutPhi;
188 fRejectExoticCluster = reco.fRejectExoticCluster;
190 fCutMinTrackPt = reco.fCutMinTrackPt;
191 fCutMinNClusterTPC = reco.fCutMinNClusterTPC;
192 fCutMinNClusterITS = reco.fCutMinNClusterITS;
193 fCutMaxChi2PerClusterTPC = reco.fCutMaxChi2PerClusterTPC;
194 fCutMaxChi2PerClusterITS = reco.fCutMaxChi2PerClusterITS;
195 fCutRequireTPCRefit = reco.fCutRequireTPCRefit;
196 fCutRequireITSRefit = reco.fCutRequireITSRefit;
197 fCutAcceptKinkDaughters = reco.fCutAcceptKinkDaughters;
198 fCutMaxDCAToVertexXY = reco.fCutMaxDCAToVertexXY;
199 fCutMaxDCAToVertexZ = reco.fCutMaxDCAToVertexZ;
200 fCutDCAToVertex2D = reco.fCutDCAToVertex2D;
202 fPIDUtils = reco.fPIDUtils;
204 fUseTimeCorrectionFactors = reco.fUseTimeCorrectionFactors;
205 fTimeCorrectionFactorsSet = reco.fTimeCorrectionFactorsSet;
208 if(reco.fResidualEta){
209 // assign or copy construct
211 *fResidualEta = *reco.fResidualEta;
213 else fResidualEta = new TArrayF(*reco.fResidualEta);
216 if(fResidualEta)delete fResidualEta;
220 if(reco.fResidualPhi){
221 // assign or copy construct
223 *fResidualPhi = *reco.fResidualPhi;
225 else fResidualPhi = new TArrayF(*reco.fResidualPhi);
228 if(fResidualPhi)delete fResidualPhi;
232 if(reco.fMatchedTrackIndex){
233 // assign or copy construct
234 if(fMatchedTrackIndex){
235 *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
237 else fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
240 if(fMatchedTrackIndex)delete fMatchedTrackIndex;
241 fMatchedTrackIndex = 0;
244 if(reco.fMatchedClusterIndex){
245 // assign or copy construct
246 if(fMatchedClusterIndex){
247 *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
249 else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
252 if(fMatchedClusterIndex)delete fMatchedClusterIndex;
253 fMatchedClusterIndex = 0;
261 //__________________________________________________
262 AliEMCALRecoUtils::~AliEMCALRecoUtils()
266 if(fEMCALRecalibrationFactors) {
267 fEMCALRecalibrationFactors->Clear();
268 delete fEMCALRecalibrationFactors;
271 if(fEMCALBadChannelMap) {
272 fEMCALBadChannelMap->Clear();
273 delete fEMCALBadChannelMap;
276 if(fMatchedTrackIndex) {delete fMatchedTrackIndex; fMatchedTrackIndex=0;}
277 if(fMatchedClusterIndex) {delete fMatchedClusterIndex; fMatchedClusterIndex=0;}
278 if(fResidualEta) {delete fResidualEta; fResidualEta=0;}
279 if(fResidualPhi) {delete fResidualPhi; fResidualPhi=0;}
283 //_______________________________________________________________
284 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells)
286 // Given the list of AbsId of the cluster, get the maximum cell and
287 // check if there are fNCellsFromBorder from the calorimeter border
289 //If the distance to the border is 0 or negative just exit accept all clusters
290 if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
292 Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
293 Bool_t shared = kFALSE;
294 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
296 AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
297 absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
299 if(absIdMax==-1) return kFALSE;
301 //Check if the cell is close to the borders:
302 Bool_t okrow = kFALSE;
303 Bool_t okcol = kFALSE;
305 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
306 AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
312 if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
315 if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
319 if(!fNoEMCALBorderAtEta0){
320 if(ieta > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
324 if(ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
327 if(ieta < 48-fNCellsFromEMCALBorder) okcol = kTRUE;
331 AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
332 fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
334 if (okcol && okrow) {
335 //printf("Accept\n");
339 //printf("Reject\n");
340 AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
347 //_________________________________________________________________________________________________________
348 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells){
349 // Check that in the cluster cells, there is no bad channel of those stored
350 // in fEMCALBadChannelMap or fPHOSBadChannelMap
352 if(!fRemoveBadChannels) return kFALSE;
353 if(!fEMCALBadChannelMap) return kFALSE;
358 for(Int_t iCell = 0; iCell<nCells; iCell++){
360 //Get the column and row
361 Int_t iTower = -1, iIphi = -1, iIeta = -1;
362 geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
363 if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
364 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
365 if(GetEMCALChannelStatus(imod, icol, irow)){
366 AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
370 }// cell cluster loop
376 //_________________________________________________
377 Bool_t AliEMCALRecoUtils::IsExoticCluster(AliVCluster *cluster){
378 // Check if the cluster has high energy but small number of cells
379 // The criteria comes from Gustavo's study
382 if(cluster->GetNCells()<(1+cluster->E()/3.))
388 //__________________________________________________
389 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
390 // Correct cluster energy from non linearity functions
391 Float_t energy = cluster->E();
393 switch (fNonLinearityFunction) {
397 //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
398 //Double_t fNonLinearityParams[0] = 1.014;
399 //Double_t fNonLinearityParams[1] = -0.03329;
400 //Double_t fNonLinearityParams[2] = -0.3853;
401 //Double_t fNonLinearityParams[3] = 0.5423;
402 //Double_t fNonLinearityParams[4] = -0.4335;
403 energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
404 ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
405 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
411 //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
412 //Double_t fNonLinearityParams[0] = 1.04;
413 //Double_t fNonLinearityParams[1] = -0.1445;
414 //Double_t fNonLinearityParams[2] = 1.046;
415 energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
419 case kPi0GammaConversion:
421 //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
422 //fNonLinearityParams[0] = 0.139393/0.1349766;
423 //fNonLinearityParams[1] = 0.0566186;
424 //fNonLinearityParams[2] = 0.982133;
425 energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
432 //From beam test, Alexei's results, for different ZS thresholds
433 // th=30 MeV; th = 45 MeV; th = 75 MeV
434 //fNonLinearityParams[0] = 1.007; 1.003; 1.002
435 //fNonLinearityParams[1] = 0.894; 0.719; 0.797
436 //fNonLinearityParams[2] = 0.246; 0.334; 0.358
437 //Rescale the param[0] with 1.03
438 energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
443 case kBeamTestCorrected:
445 //From beam test, corrected for material between beam and EMCAL
446 //fNonLinearityParams[0] = 0.99078
447 //fNonLinearityParams[1] = 0.161499;
448 //fNonLinearityParams[2] = 0.655166;
449 //fNonLinearityParams[3] = 0.134101;
450 //fNonLinearityParams[4] = 163.282;
451 //fNonLinearityParams[5] = 23.6904;
452 //fNonLinearityParams[6] = 0.978;
453 energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
459 AliDebug(2,"No correction on the energy\n");
467 //__________________________________________________
468 void AliEMCALRecoUtils::InitNonLinearityParam()
470 //Initialising Non Linearity Parameters
472 if(fNonLinearityFunction == kPi0MC)
474 fNonLinearityParams[0] = 1.014;
475 fNonLinearityParams[1] = -0.03329;
476 fNonLinearityParams[2] = -0.3853;
477 fNonLinearityParams[3] = 0.5423;
478 fNonLinearityParams[4] = -0.4335;
481 if(fNonLinearityFunction == kPi0GammaGamma)
483 fNonLinearityParams[0] = 1.04;
484 fNonLinearityParams[1] = -0.1445;
485 fNonLinearityParams[2] = 1.046;
488 if(fNonLinearityFunction == kPi0GammaConversion)
490 fNonLinearityParams[0] = 0.139393;
491 fNonLinearityParams[1] = 0.0566186;
492 fNonLinearityParams[2] = 0.982133;
495 if(fNonLinearityFunction == kBeamTest)
497 if(fNonLinearThreshold == 30)
499 fNonLinearityParams[0] = 1.007;
500 fNonLinearityParams[1] = 0.894;
501 fNonLinearityParams[2] = 0.246;
503 if(fNonLinearThreshold == 45)
505 fNonLinearityParams[0] = 1.003;
506 fNonLinearityParams[1] = 0.719;
507 fNonLinearityParams[2] = 0.334;
509 if(fNonLinearThreshold == 75)
511 fNonLinearityParams[0] = 1.002;
512 fNonLinearityParams[1] = 0.797;
513 fNonLinearityParams[2] = 0.358;
517 if(fNonLinearityFunction == kBeamTestCorrected)
519 fNonLinearityParams[0] = 0.99078;
520 fNonLinearityParams[1] = 0.161499;
521 fNonLinearityParams[2] = 0.655166;
522 fNonLinearityParams[3] = 0.134101;
523 fNonLinearityParams[4] = 163.282;
524 fNonLinearityParams[5] = 23.6904;
525 fNonLinearityParams[6] = 0.978;
529 //__________________________________________________
530 Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
532 //Calculate shower depth for a given cluster energy and particle type
542 depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
546 depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
553 gGeoManager->cd("ALIC_1/XEN1_1");
554 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
555 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
557 TGeoVolume *geoSMVol = geoSM->GetVolume();
558 TGeoShape *geoSMShape = geoSMVol->GetShape();
559 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
560 if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
561 else AliFatal("Null GEANT box");
562 }else AliFatal("NULL GEANT node matrix");
565 depth = x0 * (TMath::Log(energy*1000 / ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
571 depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
578 //__________________________________________________
579 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
580 Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
582 //For a given CaloCluster gets the absId of the cell
583 //with maximum energy deposit.
586 Double_t eCell = -1.;
587 Float_t fraction = 1.;
588 Float_t recalFactor = 1.;
589 Int_t cellAbsId = -1 ;
595 //printf("---Max?\n");
596 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
597 cellAbsId = clu->GetCellAbsId(iDig);
598 fraction = clu->GetCellAmplitudeFraction(iDig);
599 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
600 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
601 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
602 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
603 if(iDig==0) iSupMod0=iSupMod;
604 else if(iSupMod0!=iSupMod) {
606 //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
608 if(IsRecalibrationOn()) {
609 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
611 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
612 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
616 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
620 //Get from the absid the supermodule, tower and eta/phi numbers
621 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
622 //Gives SuperModule and Tower numbers
623 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
624 iIphi, iIeta,iphi,ieta);
625 //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
626 //printf("Max end---\n");
630 //________________________________________________________________
631 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
632 //Init EMCAL recalibration factors
633 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
634 //In order to avoid rewriting the same histograms
635 Bool_t oldStatus = TH1::AddDirectoryStatus();
636 TH1::AddDirectory(kFALSE);
638 fEMCALRecalibrationFactors = new TObjArray(10);
639 for (int i = 0; i < 10; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
640 //Init the histograms with 1
641 for (Int_t sm = 0; sm < 10; sm++) {
642 for (Int_t i = 0; i < 48; i++) {
643 for (Int_t j = 0; j < 24; j++) {
644 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
648 fEMCALRecalibrationFactors->SetOwner(kTRUE);
649 fEMCALRecalibrationFactors->Compress();
651 //In order to avoid rewriting the same histograms
652 TH1::AddDirectory(oldStatus);
656 //________________________________________________________________
657 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
658 //Init EMCAL bad channels map
659 AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
660 //In order to avoid rewriting the same histograms
661 Bool_t oldStatus = TH1::AddDirectoryStatus();
662 TH1::AddDirectory(kFALSE);
664 fEMCALBadChannelMap = new TObjArray(10);
665 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
666 for (int i = 0; i < 10; i++) {
667 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
672 fEMCALBadChannelMap->SetOwner(kTRUE);
673 fEMCALBadChannelMap->Compress();
675 //In order to avoid rewriting the same histograms
676 TH1::AddDirectory(oldStatus);
679 //________________________________________________________________
680 void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
681 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
683 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
684 UShort_t * index = cluster->GetCellsAbsId() ;
685 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
686 Int_t ncells = cluster->GetNCells();
688 //Initialize some used variables
691 Int_t icol = -1, irow = -1, imod=1;
692 Float_t factor = 1, frac = 0;
694 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
695 for(Int_t icell = 0; icell < ncells; icell++){
696 absId = index[icell];
697 frac = fraction[icell];
698 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
699 Int_t iTower = -1, iIphi = -1, iIeta = -1;
700 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
701 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
702 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
703 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
704 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
705 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
707 energy += cells->GetCellAmplitude(absId)*factor*frac;
711 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
713 cluster->SetE(energy);
718 //__________________________________________________
719 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
721 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
723 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
724 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
725 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
729 //__________________________________________________
730 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
732 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
733 // The algorithm is a copy of what is done in AliEMCALRecPoint
736 Float_t fraction = 1.;
737 Float_t recalFactor = 1.;
740 Int_t iTower = -1, iIphi = -1, iIeta = -1;
741 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
742 Float_t weight = 0., totalWeight=0.;
743 Float_t newPos[3] = {0,0,0};
744 Double_t pLocal[3], pGlobal[3];
745 Bool_t shared = kFALSE;
747 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
748 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
749 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
751 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
753 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
754 absId = clu->GetCellAbsId(iDig);
755 fraction = clu->GetCellAmplitudeFraction(iDig);
756 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
757 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
758 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
760 if(IsRecalibrationOn()) {
761 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
763 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
765 weight = GetCellWeight(eCell,clEnergy);
766 //printf("cell energy %f, weight %f\n",eCell,weight);
767 totalWeight += weight;
768 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
769 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
770 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
771 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
773 for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
778 for(int i=0; i<3; i++ ) newPos[i] /= totalWeight;
781 //Float_t pos[]={0,0,0};
782 //clu->GetPosition(pos);
783 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
784 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
786 if(iSupModMax > 1) {//sector 1
787 newPos[0] +=fMisalTransShift[3];//-=3.093;
788 newPos[1] +=fMisalTransShift[4];//+=6.82;
789 newPos[2] +=fMisalTransShift[5];//+=1.635;
790 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
794 newPos[0] +=fMisalTransShift[0];//+=1.134;
795 newPos[1] +=fMisalTransShift[1];//+=8.2;
796 newPos[2] +=fMisalTransShift[2];//+=1.197;
797 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
800 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
802 clu->SetPosition(newPos);
806 //__________________________________________________
807 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
809 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
810 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
813 Float_t fraction = 1.;
814 Float_t recalFactor = 1.;
818 Int_t iIphi = -1, iIeta = -1;
819 Int_t iSupMod = -1, iSupModMax = -1;
820 Int_t iphi = -1, ieta =-1;
821 Bool_t shared = kFALSE;
823 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
824 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
825 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
827 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
828 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
829 Int_t startingSM = -1;
831 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
832 absId = clu->GetCellAbsId(iDig);
833 fraction = clu->GetCellAmplitudeFraction(iDig);
834 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
835 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
836 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
838 if (iDig==0) startingSM = iSupMod;
839 else if(iSupMod != startingSM) areInSameSM = kFALSE;
841 eCell = cells->GetCellAmplitude(absId);
843 if(IsRecalibrationOn()) {
844 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
846 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
848 weight = GetCellWeight(eCell,clEnergy);
849 if(weight < 0) weight = 0;
850 totalWeight += weight;
851 weightedCol += ieta*weight;
852 weightedRow += iphi*weight;
854 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
858 Float_t xyzNew[]={0.,0.,0.};
859 if(areInSameSM == kTRUE) {
860 //printf("In Same SM\n");
861 weightedCol = weightedCol/totalWeight;
862 weightedRow = weightedRow/totalWeight;
863 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
866 //printf("In Different SM\n");
867 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
870 clu->SetPosition(xyzNew);
874 //____________________________________________________________________________
875 void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
877 //re-evaluate distance to bad channel with updated bad map
879 if(!fRecalDistToBadChannels) return;
881 //Get channels map of the supermodule where the cluster is.
882 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
883 Bool_t shared = kFALSE;
884 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
885 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
888 Float_t minDist = 10000.;
891 //Loop on tower status map
892 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
893 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
894 //Check if tower is bad.
895 if(hMap->GetBinContent(icol,irow)==0) continue;
896 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
897 // iSupMod,icol, irow, icolM,irowM);
899 dRrow=TMath::Abs(irowM-irow);
900 dRcol=TMath::Abs(icolM-icol);
901 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
903 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
910 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
915 //The only possible combinations are (0,1), (2,3) ... (8,9)
916 if(iSupMod%2) iSupMod2 = iSupMod-1;
917 else iSupMod2 = iSupMod+1;
918 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
920 //Loop on tower status map of second super module
921 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
922 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
923 //Check if tower is bad.
924 if(hMap2->GetBinContent(icol,irow)==0) continue;
925 //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",
926 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
928 dRrow=TMath::Abs(irow-irowM);
931 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
934 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
937 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
938 if(dist < minDist) minDist = dist;
943 }// shared cluster in 2 SuperModules
945 AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
946 cluster->SetDistanceToBadChannel(minDist);
950 //____________________________________________________________________________
951 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
953 //re-evaluate identification parameters with bayesian
955 if ( cluster->GetM02() != 0)
956 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
958 Float_t pidlist[AliPID::kSPECIESN+1];
959 for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
961 cluster->SetPID(pidlist);
965 //____________________________________________________________________________
966 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
968 // Calculates new center of gravity in the local EMCAL-module coordinates
969 // and tranfers into global ALICE coordinates
970 // Calculates Dispersion and main axis
975 Float_t fraction = 1.;
976 Float_t recalFactor = 1.;
996 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
998 //Get from the absid the supermodule, tower and eta/phi numbers
999 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1000 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1002 //Get the cell energy, if recalibration is on, apply factors
1003 fraction = cluster->GetCellAmplitudeFraction(iDigit);
1004 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1005 if(IsRecalibrationOn()) {
1006 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1008 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1010 if(cluster->E() > 0 && eCell > 0){
1012 w = GetCellWeight(eCell,cluster->E());
1014 etai=(Double_t)ieta;
1015 phii=(Double_t)iphi;
1020 dxx += w * etai * etai ;
1022 dzz += w * phii * phii ;
1024 dxz += w * etai * phii ;
1028 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1031 //Normalize to the weight
1037 AliError(Form("Wrong weight %f\n", wtot));
1039 //Calculate dispersion
1040 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
1042 //Get from the absid the supermodule, tower and eta/phi numbers
1043 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1044 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1046 //Get the cell energy, if recalibration is on, apply factors
1047 fraction = cluster->GetCellAmplitudeFraction(iDigit);
1048 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1049 if(IsRecalibrationOn()) {
1050 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1052 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1054 if(cluster->E() > 0 && eCell > 0){
1056 w = GetCellWeight(eCell,cluster->E());
1058 etai=(Double_t)ieta;
1059 phii=(Double_t)iphi;
1060 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
1063 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1066 //Normalize to the weigth and set shower shape parameters
1067 if (wtot > 0 && nstat > 1) {
1072 dxx -= xmean * xmean ;
1073 dzz -= zmean * zmean ;
1074 dxz -= xmean * zmean ;
1075 cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1076 cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1080 cluster->SetM20(0.) ;
1081 cluster->SetM02(0.) ;
1085 cluster->SetDispersion(TMath::Sqrt(d)) ;
1087 cluster->SetDispersion(0) ;
1090 //____________________________________________________________________________
1091 void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr, AliEMCALGeometry *geom)
1093 //This function should be called before the cluster loop
1094 //Before call this function, please recalculate the cluster positions
1095 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
1096 //Store matched cluster indexes and residuals
1098 fMatchedTrackIndex->Reset();
1099 fMatchedClusterIndex->Reset();
1100 fResidualPhi->Reset();
1101 fResidualEta->Reset();
1103 fMatchedTrackIndex->Set(500);
1104 fMatchedClusterIndex->Set(500);
1105 fResidualPhi->Set(500);
1106 fResidualEta->Set(500);
1108 AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1109 AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
1113 for (Int_t i=0; i<21;i++) cv[i]=0;
1114 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1116 AliExternalTrackParam *trackParam=0;
1118 //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
1121 AliESDtrack *esdTrack = esdevent->GetTrack(itr);
1122 if(!esdTrack || !IsAccepted(esdTrack)) continue;
1123 if(esdTrack->Pt()<fCutMinTrackPt) continue;
1124 const AliESDfriendTrack* friendTrack = esdTrack->GetFriendTrack();
1125 if(friendTrack && friendTrack->GetTPCOut())
1127 //Use TPC Out as starting point if it is available
1128 trackParam= new AliExternalTrackParam(*friendTrack->GetTPCOut());
1132 //Otherwise use TPC inner
1133 trackParam = new AliExternalTrackParam(*esdTrack->GetInnerParam());
1137 //If the input event is AOD, the starting point for extrapolation is at vertex
1138 //AOD tracks are selected according to its bit.
1141 AliAODTrack *aodTrack = aodevent->GetTrack(itr);
1142 if(!aodTrack) continue;
1143 if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
1144 if(aodTrack->Pt()<fCutMinTrackPt) continue;
1145 Double_t pos[3],mom[3];
1146 aodTrack->GetXYZ(pos);
1147 aodTrack->GetPxPyPz(mom);
1148 AliDebug(5,Form("aod track: i=%d | pos=(%5.4f,%5.4f,%5.4f) | mom=(%5.4f,%5.4f,%5.4f) | charge=%d\n",itr,pos[0],pos[1],pos[2],mom[0],mom[1],mom[2],aodTrack->Charge()));
1149 trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1152 //Return if the input data is not "AOD" or "ESD"
1155 printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
1159 if(!trackParam) continue;
1161 Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
1163 if(!clusterArr){// get clusters from event
1164 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1166 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
1167 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1168 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1169 Float_t tmpEta=-999, tmpPhi=-999;
1170 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1173 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1182 else if(fCutEtaPhiSeparate)
1184 if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1193 printf("Error: please specify your cut criteria\n");
1194 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1195 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1200 } else { // external cluster array, not from ESD event
1201 for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1203 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
1204 AliVCluster *cluster = (AliVCluster*) clusterArr->At(icl);
1205 if(!cluster->IsEMCAL()) continue;
1206 Float_t tmpEta=-999, tmpPhi=-999;
1207 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1210 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1219 else if(fCutEtaPhiSeparate)
1221 if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1230 printf("Error: please specify your cut criteria\n");
1231 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1232 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1237 }// external list of clusters
1241 fMatchedTrackIndex ->AddAt(itr,matched);
1242 fMatchedClusterIndex->AddAt(index,matched);
1243 fResidualEta ->AddAt(dEtaMax,matched);
1244 fResidualPhi ->AddAt(dPhiMax,matched);
1250 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1252 fMatchedTrackIndex ->Set(matched);
1253 fMatchedClusterIndex->Set(matched);
1254 fResidualPhi ->Set(matched);
1255 fResidualEta ->Set(matched);
1258 //________________________________________________________________________________
1259 Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom)
1262 // This function returns the index of matched cluster to input track
1263 // Returns -1 if no match is found
1266 Float_t dRMax = fCutR, dEtaMax = fCutEta, dPhiMax = fCutPhi;
1269 AliExternalTrackParam *trackParam=0;
1270 const AliESDfriendTrack* friendTrack = track->GetFriendTrack();
1271 if(friendTrack && friendTrack->GetTPCOut())
1272 trackParam= const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
1274 trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
1276 if(!trackParam) return index;
1277 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1279 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
1280 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1281 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1282 Float_t tmpEta=-999, tmpPhi=-999;
1283 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1286 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1295 else if(fCutEtaPhiSeparate)
1297 if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1306 printf("Error: please specify your cut criteria\n");
1307 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1308 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1316 //________________________________________________________________________________
1317 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpEta, Float_t &tmpPhi)
1320 //Return the residual by extrapolating a track to a cluster
1322 if(!trkParam || !cluster) return kFALSE;
1325 cluster->GetPosition(clsPos); //Has been recalculated
1326 TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
1327 Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
1328 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
1329 trkParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
1330 if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kFALSE, 0.8, -1)) return kFALSE;
1331 trkParam->GetXYZ(trkPos); //Get the extrapolated global position
1333 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1334 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
1336 Float_t clsPhi = (Float_t)clsPosVec.Phi();
1337 if(clsPhi<0) clsPhi+=2*TMath::Pi();
1338 Float_t trkPhi = (Float_t)trkPosVec.Phi();
1339 if(trkPhi<0) trkPhi+=2*TMath::Pi();
1340 tmpPhi = clsPhi-trkPhi; // track cluster matching
1341 tmpEta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
1346 //________________________________________________________________________________
1347 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
1349 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1350 //Get the residuals dEta and dPhi for this cluster to the closest track
1351 //Works with ESDs and AODs
1353 if( FindMatchedPosForCluster(clsIndex) >= 999 )
1355 AliDebug(2,"No matched tracks found!\n");
1360 dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
1361 dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
1363 //________________________________________________________________________________
1364 void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
1366 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1367 //Get the residuals dEta and dPhi for this track to the closest cluster
1368 //Works with ESDs and AODs
1370 if( FindMatchedPosForTrack(trkIndex) >= 999 )
1372 AliDebug(2,"No matched cluster found!\n");
1377 dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
1378 dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
1381 //__________________________________________________________
1382 Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
1384 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1385 //Get the index of matched track to this cluster
1386 //Works with ESDs and AODs
1388 if(IsClusterMatched(clsIndex))
1389 return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
1394 //__________________________________________________________
1395 Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
1397 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1398 //Get the index of matched cluster to this track
1399 //Works with ESDs and AODs
1401 if(IsTrackMatched(trkIndex))
1402 return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
1407 //__________________________________________________
1408 Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex)
1410 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1411 //Returns if the cluster has a match
1412 if(FindMatchedPosForCluster(clsIndex) < 999)
1418 //__________________________________________________
1419 Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex)
1421 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1422 //Returns if the track has a match
1423 if(FindMatchedPosForTrack(trkIndex) < 999)
1429 //__________________________________________________________
1430 UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
1432 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1433 //Returns the position of the match in the fMatchedClusterIndex array
1434 Float_t tmpR = fCutR;
1437 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
1439 if(fMatchedClusterIndex->At(i)==clsIndex)
1441 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1446 AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1453 //__________________________________________________________
1454 UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
1456 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1457 //Returns the position of the match in the fMatchedTrackIndex array
1458 Float_t tmpR = fCutR;
1461 for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
1463 if(fMatchedTrackIndex->At(i)==trkIndex)
1465 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1470 AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1477 //__________________________________________________________
1478 Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells)
1480 // check if the cluster survives some quality cut
1483 Bool_t isGood=kTRUE;
1484 if(!cluster || !cluster->IsEMCAL()) return kFALSE;
1485 if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
1486 if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
1487 if(fRejectExoticCluster && IsExoticCluster(cluster)) return kFALSE;
1492 //__________________________________________________________
1493 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1495 // Given a esd track, return whether the track survive all the cuts
1497 // The different quality parameter are first
1498 // retrieved from the track. then it is found out what cuts the
1499 // track did not survive and finally the cuts are imposed.
1501 UInt_t status = esdTrack->GetStatus();
1503 Int_t nClustersITS = esdTrack->GetITSclusters(0);
1504 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1506 Float_t chi2PerClusterITS = -1;
1507 Float_t chi2PerClusterTPC = -1;
1508 if (nClustersITS!=0)
1509 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1510 if (nClustersTPC!=0)
1511 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1515 Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1516 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1517 SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1522 esdTrack->GetImpactParameters(b,bCov);
1523 if (bCov[0]<=0 || bCov[2]<=0) {
1524 AliDebug(1, "Estimated b resolution lower or equal zero!");
1525 bCov[0]=0; bCov[2]=0;
1528 Float_t dcaToVertexXY = b[0];
1529 Float_t dcaToVertexZ = b[1];
1530 Float_t dcaToVertex = -1;
1532 if (fCutDCAToVertex2D)
1533 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1535 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1539 Bool_t cuts[kNCuts];
1540 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1542 // track quality cuts
1543 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1545 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1547 if (nClustersTPC<fCutMinNClusterTPC)
1549 if (nClustersITS<fCutMinNClusterITS)
1551 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1553 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1555 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1557 if (fCutDCAToVertex2D && dcaToVertex > 1)
1559 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1561 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1564 //Require at least one SPD point + anything else in ITS
1565 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1569 for (Int_t i=0; i<kNCuts; i++)
1570 if (cuts[i]) {cut = kTRUE;}
1578 //__________________________________________________
1579 void AliEMCALRecoUtils::InitTrackCuts()
1581 //Intilize the track cut criteria
1582 //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1583 //Also you can customize the cuts using the setters
1586 SetMinNClustersTPC(70);
1587 SetMaxChi2PerClusterTPC(4);
1588 SetAcceptKinkDaughters(kFALSE);
1589 SetRequireTPCRefit(kTRUE);
1592 SetRequireITSRefit(kTRUE);
1593 SetMaxDCAToVertexZ(2);
1594 SetDCAToVertex2D(kFALSE);
1595 SetMaxChi2PerClusterITS(); //which by default sets the value to 1e10.
1596 SetMinNClustersITS();
1599 //___________________________________________________
1600 void AliEMCALRecoUtils::Print(const Option_t *) const
1604 printf("AliEMCALRecoUtils Settings: \n");
1605 printf("Misalignment shifts\n");
1606 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,
1607 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1608 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
1609 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1610 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
1612 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
1614 printf("Matching criteria: ");
1617 printf("sqrt(dEta^2+dPhi^2)<%2.2f\n",fCutR);
1619 else if(fCutEtaPhiSeparate)
1621 printf("dEta<%2.2f, dPhi<%2.2f\n",fCutEta,fCutPhi);
1626 printf("please specify your cut criteria\n");
1627 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1628 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1631 printf("Mass hypothesis = %2.3f[GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
1633 printf("Track cuts: \n");
1634 printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
1635 printf("AOD track selection mask: %d\n",fAODFilterMask);
1636 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1637 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1638 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1639 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1640 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1644 //_____________________________________________________________________
1645 void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1646 //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1647 //Do it only once and only if it is requested
1649 if(!fUseTimeCorrectionFactors) return;
1650 if(fTimeCorrectionFactorsSet) return;
1652 printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1654 AliEMCALCalibTimeDepCorrection *corr = new AliEMCALCalibTimeDepCorrection();
1655 corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1657 SwitchOnRecalibration();
1658 for(Int_t ism = 0; ism < 4; ism++){
1659 for(Int_t icol = 0; icol < 48; icol++){
1660 for(Int_t irow = 0; irow < 24; irow++){
1661 Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1662 Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1663 GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1664 //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
1665 // orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1666 // (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1670 fTimeCorrectionFactorsSet = kTRUE;