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 ///////////////////////////////////////////////////////////////////////////////
31 // standard C++ includes
32 //#include <Riostream.h>
35 #include <TGeoManager.h>
36 #include <TGeoMatrix.h>
40 #include "AliVCluster.h"
41 #include "AliVCaloCells.h"
42 #include "AliVEvent.h"
45 #include "AliESDEvent.h"
46 #include "AliAODEvent.h"
47 #include "AliESDtrack.h"
48 #include "AliAODTrack.h"
49 #include "AliExternalTrackParam.h"
50 #include "AliESDfriendTrack.h"
51 #include "AliTrackerBase.h"
54 #include "AliEMCALRecoUtils.h"
55 #include "AliEMCALGeometry.h"
56 #include "AliEMCALTrack.h"
57 #include "AliEMCALCalibTimeDepCorrection.h"
58 #include "AliEMCALPIDUtils.h"
60 ClassImp(AliEMCALRecoUtils)
62 //______________________________________________
63 AliEMCALRecoUtils::AliEMCALRecoUtils():
64 fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
65 fPosAlgo(kUnchanged),fW0(4.), fNonLinearThreshold(30),
66 fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
67 fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
68 fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
70 fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
71 fResidualZ(0x0), fResidualR(0x0), fCutR(10), fCutZ(10), fMass(0.139), fStep(1),
72 fRejectExoticCluster(kFALSE),
73 fCutMinNClusterTPC(0), fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
74 fCutRequireTPCRefit(0), fCutRequireITSRefit(0), fCutAcceptKinkDaughters(0),
75 fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0),fCutDCAToVertex2D(0),fPIDUtils(),
76 fUseTimeCorrectionFactors(kFALSE), fTimeCorrectionFactorsSet(kFALSE)
80 // Initialize all constant values which have to be used
81 // during Reco algorithm execution
84 //Misalignment matrices
85 for(Int_t i = 0; i < 15 ; i++) {
86 fMisalTransShift[i] = 0.;
87 fMisalRotShift[i] = 0.;
91 for(Int_t i = 0; i < 7 ; i++) fNonLinearityParams[i] = 0.;
93 //For kBeamTestCorrected case, but default is no correction
94 fNonLinearityParams[0] = 0.99078;
95 fNonLinearityParams[1] = 0.161499;
96 fNonLinearityParams[2] = 0.655166;
97 fNonLinearityParams[3] = 0.134101;
98 fNonLinearityParams[4] = 163.282;
99 fNonLinearityParams[5] = 23.6904;
100 fNonLinearityParams[6] = 0.978;
102 //For kPi0GammaGamma case
103 //fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
104 //fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
105 //fNonLinearityParams[2] = 1.046;
108 fMatchedTrackIndex = new TArrayI();
109 fMatchedClusterIndex = new TArrayI();
110 fResidualZ = new TArrayF();
111 fResidualR = new TArrayF();
115 fPIDUtils = new AliEMCALPIDUtils();
120 //______________________________________________________________________
121 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
122 : TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction),
123 fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0), fNonLinearThreshold(reco.fNonLinearThreshold),
124 fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
125 fRemoveBadChannels(reco.fRemoveBadChannels),fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
126 fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
127 fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
128 fAODFilterMask(reco.fAODFilterMask),
129 fMatchedTrackIndex(reco.fMatchedTrackIndex?new TArrayI(*reco.fMatchedTrackIndex):0x0),
130 fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
131 fResidualZ(reco.fResidualZ?new TArrayF(*reco.fResidualZ):0x0),
132 fResidualR(reco.fResidualR?new TArrayF(*reco.fResidualR):0x0),
133 fCutR(reco.fCutR),fCutZ(reco.fCutZ),fMass(reco.fMass), fStep(reco.fStep),
134 fRejectExoticCluster(reco.fRejectExoticCluster),
135 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;
185 fRejectExoticCluster = reco.fRejectExoticCluster;
187 fCutMinNClusterTPC = reco.fCutMinNClusterTPC;
188 fCutMinNClusterITS = reco.fCutMinNClusterITS;
189 fCutMaxChi2PerClusterTPC = reco.fCutMaxChi2PerClusterTPC;
190 fCutMaxChi2PerClusterITS = reco.fCutMaxChi2PerClusterITS;
191 fCutRequireTPCRefit = reco.fCutRequireTPCRefit;
192 fCutRequireITSRefit = reco.fCutRequireITSRefit;
193 fCutAcceptKinkDaughters = reco.fCutAcceptKinkDaughters;
194 fCutMaxDCAToVertexXY = reco.fCutMaxDCAToVertexXY;
195 fCutMaxDCAToVertexZ = reco.fCutMaxDCAToVertexZ;
196 fCutDCAToVertex2D = reco.fCutDCAToVertex2D;
198 fPIDUtils = reco.fPIDUtils;
200 fUseTimeCorrectionFactors = reco.fUseTimeCorrectionFactors;
201 fTimeCorrectionFactorsSet = reco.fTimeCorrectionFactorsSet;
205 // assign or copy construct
207 *fResidualR = *reco.fResidualR;
209 else fResidualR = new TArrayF(*reco.fResidualR);
212 if(fResidualR)delete fResidualR;
217 // assign or copy construct
219 *fResidualZ = *reco.fResidualZ;
221 else fResidualZ = new TArrayF(*reco.fResidualZ);
224 if(fResidualZ)delete fResidualZ;
228 if(reco.fMatchedTrackIndex){
229 // assign or copy construct
230 if(fMatchedTrackIndex){
231 *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
233 else fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
236 if(fMatchedTrackIndex)delete fMatchedTrackIndex;
237 fMatchedTrackIndex = 0;
240 if(reco.fMatchedClusterIndex){
241 // assign or copy construct
242 if(fMatchedClusterIndex){
243 *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
245 else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
248 if(fMatchedClusterIndex)delete fMatchedClusterIndex;
249 fMatchedClusterIndex = 0;
257 //__________________________________________________
258 AliEMCALRecoUtils::~AliEMCALRecoUtils()
262 if(fEMCALRecalibrationFactors) {
263 fEMCALRecalibrationFactors->Clear();
264 delete fEMCALRecalibrationFactors;
267 if(fEMCALBadChannelMap) {
268 fEMCALBadChannelMap->Clear();
269 delete fEMCALBadChannelMap;
272 if(fMatchedTrackIndex) {delete fMatchedTrackIndex; fMatchedTrackIndex=0;}
273 if(fMatchedClusterIndex) {delete fMatchedClusterIndex; fMatchedClusterIndex=0;}
274 if(fResidualR) {delete fResidualR; fResidualR=0;}
275 if(fResidualZ) {delete fResidualZ; fResidualZ=0;}
279 //_______________________________________________________________
280 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells)
282 // Given the list of AbsId of the cluster, get the maximum cell and
283 // check if there are fNCellsFromBorder from the calorimeter border
285 //If the distance to the border is 0 or negative just exit accept all clusters
286 if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
288 Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
289 Bool_t shared = kFALSE;
290 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
292 AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
293 absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
295 if(absIdMax==-1) return kFALSE;
297 //Check if the cell is close to the borders:
298 Bool_t okrow = kFALSE;
299 Bool_t okcol = kFALSE;
301 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
302 AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
308 if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
311 if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
315 if(!fNoEMCALBorderAtEta0){
316 if(ieta > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
320 if(ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
323 if(ieta < 48-fNCellsFromEMCALBorder) okcol = kTRUE;
327 AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
328 fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
330 if (okcol && okrow) {
331 //printf("Accept\n");
335 //printf("Reject\n");
336 AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
343 //_________________________________________________________________________________________________________
344 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells){
345 // Check that in the cluster cells, there is no bad channel of those stored
346 // in fEMCALBadChannelMap or fPHOSBadChannelMap
348 if(!fRemoveBadChannels) return kFALSE;
349 if(!fEMCALBadChannelMap) return kFALSE;
354 for(Int_t iCell = 0; iCell<nCells; iCell++){
356 //Get the column and row
357 Int_t iTower = -1, iIphi = -1, iIeta = -1;
358 geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
359 if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
360 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
361 if(GetEMCALChannelStatus(imod, icol, irow)){
362 AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
366 }// cell cluster loop
372 //_________________________________________________
373 Bool_t AliEMCALRecoUtils::IsExoticCluster(AliVCluster *cluster){
374 // Check if the cluster has high energy but small number of cells
375 // The criteria comes from Gustavo's study
378 if(cluster->GetNCells()<(1+cluster->E()/3.))
384 //__________________________________________________
385 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
386 // Correct cluster energy from non linearity functions
387 Float_t energy = cluster->E();
389 switch (fNonLinearityFunction) {
393 //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
394 //Double_t fNonLinearityParams[0] = 1.014;
395 //Double_t fNonLinearityParams[1] = -0.03329;
396 //Double_t fNonLinearityParams[2] = -0.3853;
397 //Double_t fNonLinearityParams[3] = 0.5423;
398 //Double_t fNonLinearityParams[4] = -0.4335;
399 energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
400 ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
401 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
407 //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
408 //Double_t fNonLinearityParams[0] = 1.04;
409 //Double_t fNonLinearityParams[1] = -0.1445;
410 //Double_t fNonLinearityParams[2] = 1.046;
411 energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
415 case kPi0GammaConversion:
417 //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
418 //fNonLinearityParams[0] = 0.139393/0.1349766;
419 //fNonLinearityParams[1] = 0.0566186;
420 //fNonLinearityParams[2] = 0.982133;
421 energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
428 //From beam test, Alexei's results, for different ZS thresholds
429 // th=30 MeV; th = 45 MeV; th = 75 MeV
430 //fNonLinearityParams[0] = 1.007; 1.003; 1.002
431 //fNonLinearityParams[1] = 0.894; 0.719; 0.797
432 //fNonLinearityParams[2] = 0.246; 0.334; 0.358
433 //Rescale the param[0] with 1.03
434 energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
439 case kBeamTestCorrected:
441 //From beam test, corrected for material between beam and EMCAL
442 //fNonLinearityParams[0] = 0.99078
443 //fNonLinearityParams[1] = 0.161499;
444 //fNonLinearityParams[2] = 0.655166;
445 //fNonLinearityParams[3] = 0.134101;
446 //fNonLinearityParams[4] = 163.282;
447 //fNonLinearityParams[5] = 23.6904;
448 //fNonLinearityParams[6] = 0.978;
449 energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
455 AliDebug(2,"No correction on the energy\n");
463 //__________________________________________________
464 void AliEMCALRecoUtils::InitNonLinearityParam()
466 //Initialising Non Linearity Parameters
468 if(fNonLinearityFunction == kPi0MC)
470 fNonLinearityParams[0] = 1.014;
471 fNonLinearityParams[1] = -0.03329;
472 fNonLinearityParams[2] = -0.3853;
473 fNonLinearityParams[3] = 0.5423;
474 fNonLinearityParams[4] = -0.4335;
477 if(fNonLinearityFunction == kPi0GammaGamma)
479 fNonLinearityParams[0] = 1.04;
480 fNonLinearityParams[1] = -0.1445;
481 fNonLinearityParams[2] = 1.046;
484 if(fNonLinearityFunction == kPi0GammaConversion)
486 fNonLinearityParams[0] = 0.139393;
487 fNonLinearityParams[1] = 0.0566186;
488 fNonLinearityParams[2] = 0.982133;
491 if(fNonLinearityFunction == kBeamTest)
493 if(fNonLinearThreshold == 30)
495 fNonLinearityParams[0] = 1.007;
496 fNonLinearityParams[1] = 0.894;
497 fNonLinearityParams[2] = 0.246;
499 if(fNonLinearThreshold == 45)
501 fNonLinearityParams[0] = 1.003;
502 fNonLinearityParams[1] = 0.719;
503 fNonLinearityParams[2] = 0.334;
505 if(fNonLinearThreshold == 75)
507 fNonLinearityParams[0] = 1.002;
508 fNonLinearityParams[1] = 0.797;
509 fNonLinearityParams[2] = 0.358;
513 if(fNonLinearityFunction == kBeamTestCorrected)
515 fNonLinearityParams[0] = 0.99078;
516 fNonLinearityParams[1] = 0.161499;
517 fNonLinearityParams[2] = 0.655166;
518 fNonLinearityParams[3] = 0.134101;
519 fNonLinearityParams[4] = 163.282;
520 fNonLinearityParams[5] = 23.6904;
521 fNonLinearityParams[6] = 0.978;
525 //__________________________________________________
526 Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
528 //Calculate shower depth for a given cluster energy and particle type
538 depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
542 depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
549 gGeoManager->cd("ALIC_1/XEN1_1");
550 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
551 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
553 TGeoVolume *geoSMVol = geoSM->GetVolume();
554 TGeoShape *geoSMShape = geoSMVol->GetShape();
555 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
556 if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
557 else AliFatal("Null GEANT box");
558 }else AliFatal("NULL GEANT node matrix");
561 depth = x0 * (TMath::Log(energy*1000 / ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
567 depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
574 //__________________________________________________
575 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
576 Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
578 //For a given CaloCluster gets the absId of the cell
579 //with maximum energy deposit.
582 Double_t eCell = -1.;
583 Float_t fraction = 1.;
584 Float_t recalFactor = 1.;
585 Int_t cellAbsId = -1 ;
591 //printf("---Max?\n");
592 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
593 cellAbsId = clu->GetCellAbsId(iDig);
594 fraction = clu->GetCellAmplitudeFraction(iDig);
595 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
596 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
597 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
598 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
599 if(iDig==0) iSupMod0=iSupMod;
600 else if(iSupMod0!=iSupMod) {
602 //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
604 if(IsRecalibrationOn()) {
605 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
607 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
608 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
612 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
616 //Get from the absid the supermodule, tower and eta/phi numbers
617 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
618 //Gives SuperModule and Tower numbers
619 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
620 iIphi, iIeta,iphi,ieta);
621 //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
622 //printf("Max end---\n");
626 //________________________________________________________________
627 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
628 //Init EMCAL recalibration factors
629 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
630 //In order to avoid rewriting the same histograms
631 Bool_t oldStatus = TH1::AddDirectoryStatus();
632 TH1::AddDirectory(kFALSE);
634 fEMCALRecalibrationFactors = new TObjArray(10);
635 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));
636 //Init the histograms with 1
637 for (Int_t sm = 0; sm < 10; sm++) {
638 for (Int_t i = 0; i < 48; i++) {
639 for (Int_t j = 0; j < 24; j++) {
640 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
644 fEMCALRecalibrationFactors->SetOwner(kTRUE);
645 fEMCALRecalibrationFactors->Compress();
647 //In order to avoid rewriting the same histograms
648 TH1::AddDirectory(oldStatus);
652 //________________________________________________________________
653 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
654 //Init EMCAL bad channels map
655 AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
656 //In order to avoid rewriting the same histograms
657 Bool_t oldStatus = TH1::AddDirectoryStatus();
658 TH1::AddDirectory(kFALSE);
660 fEMCALBadChannelMap = new TObjArray(10);
661 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
662 for (int i = 0; i < 10; i++) {
663 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
668 fEMCALBadChannelMap->SetOwner(kTRUE);
669 fEMCALBadChannelMap->Compress();
671 //In order to avoid rewriting the same histograms
672 TH1::AddDirectory(oldStatus);
675 //________________________________________________________________
676 void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
677 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
679 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
680 UShort_t * index = cluster->GetCellsAbsId() ;
681 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
682 Int_t ncells = cluster->GetNCells();
684 //Initialize some used variables
687 Int_t icol = -1, irow = -1, imod=1;
688 Float_t factor = 1, frac = 0;
690 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
691 for(Int_t icell = 0; icell < ncells; icell++){
692 absId = index[icell];
693 frac = fraction[icell];
694 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
695 Int_t iTower = -1, iIphi = -1, iIeta = -1;
696 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
697 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
698 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
699 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
700 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
701 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
703 energy += cells->GetCellAmplitude(absId)*factor*frac;
707 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
709 cluster->SetE(energy);
714 //__________________________________________________
715 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
717 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
719 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
720 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
721 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
725 //__________________________________________________
726 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
728 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
729 // The algorithm is a copy of what is done in AliEMCALRecPoint
732 Float_t fraction = 1.;
733 Float_t recalFactor = 1.;
736 Int_t iTower = -1, iIphi = -1, iIeta = -1;
737 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
738 Float_t weight = 0., totalWeight=0.;
739 Float_t newPos[3] = {0,0,0};
740 Double_t pLocal[3], pGlobal[3];
741 Bool_t shared = kFALSE;
743 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
744 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
745 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
747 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
749 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
750 absId = clu->GetCellAbsId(iDig);
751 fraction = clu->GetCellAmplitudeFraction(iDig);
752 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
753 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
754 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
756 if(IsRecalibrationOn()) {
757 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
759 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
761 weight = GetCellWeight(eCell,clEnergy);
762 //printf("cell energy %f, weight %f\n",eCell,weight);
763 totalWeight += weight;
764 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
765 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
766 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
767 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
769 for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
774 for(int i=0; i<3; i++ ) newPos[i] /= totalWeight;
777 //Float_t pos[]={0,0,0};
778 //clu->GetPosition(pos);
779 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
780 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
782 if(iSupModMax > 1) {//sector 1
783 newPos[0] +=fMisalTransShift[3];//-=3.093;
784 newPos[1] +=fMisalTransShift[4];//+=6.82;
785 newPos[2] +=fMisalTransShift[5];//+=1.635;
786 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
790 newPos[0] +=fMisalTransShift[0];//+=1.134;
791 newPos[1] +=fMisalTransShift[1];//+=8.2;
792 newPos[2] +=fMisalTransShift[2];//+=1.197;
793 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
796 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
798 clu->SetPosition(newPos);
802 //__________________________________________________
803 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
805 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
806 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
809 Float_t fraction = 1.;
810 Float_t recalFactor = 1.;
814 Int_t iIphi = -1, iIeta = -1;
815 Int_t iSupMod = -1, iSupModMax = -1;
816 Int_t iphi = -1, ieta =-1;
817 Bool_t shared = kFALSE;
819 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
820 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
821 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
823 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
824 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
825 Int_t startingSM = -1;
827 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
828 absId = clu->GetCellAbsId(iDig);
829 fraction = clu->GetCellAmplitudeFraction(iDig);
830 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
831 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
832 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
834 if (iDig==0) startingSM = iSupMod;
835 else if(iSupMod != startingSM) areInSameSM = kFALSE;
837 eCell = cells->GetCellAmplitude(absId);
839 if(IsRecalibrationOn()) {
840 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
842 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
844 weight = GetCellWeight(eCell,clEnergy);
845 if(weight < 0) weight = 0;
846 totalWeight += weight;
847 weightedCol += ieta*weight;
848 weightedRow += iphi*weight;
850 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
854 Float_t xyzNew[]={0.,0.,0.};
855 if(areInSameSM == kTRUE) {
856 //printf("In Same SM\n");
857 weightedCol = weightedCol/totalWeight;
858 weightedRow = weightedRow/totalWeight;
859 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
862 //printf("In Different SM\n");
863 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
866 clu->SetPosition(xyzNew);
870 //____________________________________________________________________________
871 void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
873 //re-evaluate distance to bad channel with updated bad map
875 if(!fRecalDistToBadChannels) return;
877 //Get channels map of the supermodule where the cluster is.
878 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
879 Bool_t shared = kFALSE;
880 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
881 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
884 Float_t minDist = 10000.;
887 //Loop on tower status map
888 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
889 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
890 //Check if tower is bad.
891 if(hMap->GetBinContent(icol,irow)==0) continue;
892 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
893 // iSupMod,icol, irow, icolM,irowM);
895 dRrow=TMath::Abs(irowM-irow);
896 dRcol=TMath::Abs(icolM-icol);
897 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
899 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
906 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
911 //The only possible combinations are (0,1), (2,3) ... (8,9)
912 if(iSupMod%2) iSupMod2 = iSupMod-1;
913 else iSupMod2 = iSupMod+1;
914 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
916 //Loop on tower status map of second super module
917 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
918 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
919 //Check if tower is bad.
920 if(hMap2->GetBinContent(icol,irow)==0) continue;
921 //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",
922 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
924 dRrow=TMath::Abs(irow-irowM);
927 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
930 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
933 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
934 if(dist < minDist) minDist = dist;
939 }// shared cluster in 2 SuperModules
941 AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
942 cluster->SetDistanceToBadChannel(minDist);
946 //____________________________________________________________________________
947 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
949 //re-evaluate identification parameters with bayesian
951 if ( cluster->GetM02() != 0)
952 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
954 Float_t pidlist[AliPID::kSPECIESN+1];
955 for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
957 cluster->SetPID(pidlist);
961 //____________________________________________________________________________
962 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
964 // Calculates new center of gravity in the local EMCAL-module coordinates
965 // and tranfers into global ALICE coordinates
966 // Calculates Dispersion and main axis
971 Float_t fraction = 1.;
972 Float_t recalFactor = 1.;
992 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
994 //Get from the absid the supermodule, tower and eta/phi numbers
995 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
996 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
998 //Get the cell energy, if recalibration is on, apply factors
999 fraction = cluster->GetCellAmplitudeFraction(iDigit);
1000 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1001 if(IsRecalibrationOn()) {
1002 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1004 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1006 if(cluster->E() > 0 && eCell > 0){
1008 w = GetCellWeight(eCell,cluster->E());
1010 etai=(Double_t)ieta;
1011 phii=(Double_t)iphi;
1016 dxx += w * etai * etai ;
1018 dzz += w * phii * phii ;
1020 dxz += w * etai * phii ;
1024 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1027 //Normalize to the weight
1033 AliError(Form("Wrong weight %f\n", wtot));
1035 //Calculate dispersion
1036 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
1038 //Get from the absid the supermodule, tower and eta/phi numbers
1039 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1040 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1042 //Get the cell energy, if recalibration is on, apply factors
1043 fraction = cluster->GetCellAmplitudeFraction(iDigit);
1044 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1045 if(IsRecalibrationOn()) {
1046 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1048 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1050 if(cluster->E() > 0 && eCell > 0){
1052 w = GetCellWeight(eCell,cluster->E());
1054 etai=(Double_t)ieta;
1055 phii=(Double_t)iphi;
1056 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
1059 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1062 //Normalize to the weigth and set shower shape parameters
1063 if (wtot > 0 && nstat > 1) {
1068 dxx -= xmean * xmean ;
1069 dzz -= zmean * zmean ;
1070 dxz -= xmean * zmean ;
1071 cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1072 cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1076 cluster->SetM20(0.) ;
1077 cluster->SetM02(0.) ;
1081 cluster->SetDispersion(TMath::Sqrt(d)) ;
1083 cluster->SetDispersion(0) ;
1086 //____________________________________________________________________________
1087 void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr, AliEMCALGeometry *geom)
1089 //This function should be called before the cluster loop
1090 //Before call this function, please recalculate the cluster positions
1091 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
1092 //Store matched cluster indexes and residuals
1094 fMatchedTrackIndex ->Reset();
1095 fMatchedClusterIndex->Reset();
1096 fResidualZ ->Reset();
1097 fResidualR ->Reset();
1099 fMatchedTrackIndex ->Set(500);
1100 fMatchedClusterIndex->Set(500);
1101 fResidualZ ->Set(500);
1102 fResidualR ->Set(500);
1104 AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1105 AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
1109 for (Int_t i=0; i<21;i++) cv[i]=0;
1110 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1112 AliExternalTrackParam *trackParam=0;
1114 //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
1117 AliESDtrack *esdTrack = esdevent->GetTrack(itr);
1118 if(!esdTrack || !IsAccepted(esdTrack)) continue;
1119 const AliESDfriendTrack* friendTrack = esdTrack->GetFriendTrack();
1120 if(friendTrack && friendTrack->GetTPCOut())
1122 //Use TPC Out as starting point if it is available
1123 trackParam= new AliExternalTrackParam(*friendTrack->GetTPCOut());
1127 //Otherwise use TPC inner
1128 trackParam = new AliExternalTrackParam(*esdTrack->GetInnerParam());
1132 //If the input event is AOD, the starting point for extrapolation is at vertex
1133 //AOD tracks are selected according to its bit.
1136 AliAODTrack *aodTrack = aodevent->GetTrack(itr);
1137 if(!aodTrack) continue;
1138 if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
1139 Double_t pos[3],mom[3];
1140 aodTrack->GetXYZ(pos);
1141 aodTrack->GetPxPyPz(mom);
1142 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()));
1143 trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1146 //Return if the input data is not "AOD" or "ESD"
1149 printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
1153 if(!trackParam) continue;
1155 Float_t dRMax = fCutR, dZMax=fCutZ;
1157 if(!clusterArr){// get clusters from event
1158 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1160 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
1161 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1162 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1163 Float_t tmpR=-1, tmpZ=-1;
1164 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
1173 } else { // external cluster array, not from ESD event
1174 for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1176 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
1177 AliVCluster *cluster = (AliVCluster*) clusterArr->At(icl);
1178 if(!cluster->IsEMCAL()) continue;
1179 Float_t tmpR=-1, tmpZ=-1;
1180 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
1189 }// external list of clusters
1193 fMatchedTrackIndex ->AddAt(itr,matched);
1194 fMatchedClusterIndex->AddAt(index,matched);
1195 fResidualZ ->AddAt(dZMax,matched);
1196 fResidualR ->AddAt(dRMax,matched);
1202 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1204 fMatchedTrackIndex ->Set(matched);
1205 fMatchedClusterIndex->Set(matched);
1206 fResidualZ ->Set(matched);
1207 fResidualR ->Set(matched);
1210 //________________________________________________________________________________
1211 Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom)
1214 // This function returns the index of matched cluster to input track
1215 // Cut on match is dR<10cm by default. Returns -1 if no match is found
1218 Float_t dRMax = fCutR;
1221 AliExternalTrackParam *trackParam=0;
1222 const AliESDfriendTrack* friendTrack = track->GetFriendTrack();
1223 if(friendTrack && friendTrack->GetTPCOut())
1224 trackParam= const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
1226 trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
1228 if(!trackParam) return index;
1229 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1231 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
1232 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1233 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1234 Float_t tmpR=-1, tmpZ=-1;
1235 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
1236 if(tmpR>-1 && tmpR<dRMax)
1246 //________________________________________________________________________________
1247 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpR, Float_t &tmpZ)
1250 //Return the residual by extrapolating a track to a cluster
1252 if(!trkParam || !cluster) return kFALSE;
1255 cluster->GetPosition(clsPos); //Has been recalculated
1256 TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
1257 Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
1258 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
1259 trkParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
1260 if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kFALSE)) return kFALSE;
1261 trkParam->GetXYZ(trkPos); //Get the extrapolated global position
1262 tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
1263 tmpZ = clsPos[2]-trkPos[2];
1267 //________________________________________________________________________________
1268 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dR, Float_t &dZ)
1270 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1271 //Get the residuals dR and dZ for this cluster to the closest track
1272 //Works with ESDs and AODs
1274 if( FindMatchedPosForCluster(clsIndex) >= 999 )
1276 AliDebug(2,"No matched tracks found!\n");
1281 dR = fResidualR->At(FindMatchedPosForCluster(clsIndex));
1282 dZ = fResidualZ->At(FindMatchedPosForCluster(clsIndex));
1284 //________________________________________________________________________________
1285 void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dR, Float_t &dZ)
1287 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1288 //Get the residuals dR and dZ for this track to the closest cluster
1289 //Works with ESDs and AODs
1291 if( FindMatchedPosForTrack(trkIndex) >= 999 )
1293 AliDebug(2,"No matched cluster found!\n");
1298 dR = fResidualR->At(FindMatchedPosForTrack(trkIndex));
1299 dZ = fResidualZ->At(FindMatchedPosForTrack(trkIndex));
1302 //__________________________________________________________
1303 Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
1305 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1306 //Get the index of matched track to this cluster
1307 //Works with ESDs and AODs
1309 if(IsClusterMatched(clsIndex))
1310 return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
1315 //__________________________________________________________
1316 Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
1318 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1319 //Get the index of matched cluster to this track
1320 //Works with ESDs and AODs
1322 if(IsTrackMatched(trkIndex))
1323 return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
1328 //__________________________________________________
1329 Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex)
1331 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1332 //Returns if the cluster has a match
1333 if(FindMatchedPosForCluster(clsIndex) < 999)
1339 //__________________________________________________
1340 Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex)
1342 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1343 //Returns if the track has a match
1344 if(FindMatchedPosForTrack(trkIndex) < 999)
1350 //__________________________________________________________
1351 UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
1353 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1354 //Returns the position of the match in the fMatchedClusterIndex array
1355 Float_t tmpR = fCutR;
1358 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
1360 if(fMatchedClusterIndex->At(i)==clsIndex && fResidualR->At(i)<tmpR)
1363 tmpR=fResidualR->At(i);
1364 AliDebug(3,Form("Matched cluster index: index: %d, dR: %2.4f, dZ: %2.4f.\n",fMatchedClusterIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
1370 //__________________________________________________________
1371 UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
1373 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1374 //Returns the position of the match in the fMatchedTrackIndex array
1375 Float_t tmpR = fCutR;
1378 for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
1380 if(fMatchedTrackIndex->At(i)==trkIndex && fResidualR->At(i)<tmpR)
1383 tmpR=fResidualR->At(i);
1384 AliDebug(3,Form("Matched track index: index: %d, dR: %2.4f, dZ: %2.4f.\n",fMatchedTrackIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
1390 //__________________________________________________________
1391 Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells)
1393 // check if the cluster survives some quality cut
1396 Bool_t isGood=kTRUE;
1397 if(!cluster || !cluster->IsEMCAL()) return kFALSE;
1398 if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) isGood=kFALSE;
1399 if(!CheckCellFiducialRegion(geom,cluster,cells)) isGood=kFALSE;
1400 if(fRejectExoticCluster && IsExoticCluster(cluster)) isGood=kFALSE;
1405 //__________________________________________________________
1406 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1408 // Given a esd track, return whether the track survive all the cuts
1410 // The different quality parameter are first
1411 // retrieved from the track. then it is found out what cuts the
1412 // track did not survive and finally the cuts are imposed.
1414 UInt_t status = esdTrack->GetStatus();
1416 Int_t nClustersITS = esdTrack->GetITSclusters(0);
1417 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1419 Float_t chi2PerClusterITS = -1;
1420 Float_t chi2PerClusterTPC = -1;
1421 if (nClustersITS!=0)
1422 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1423 if (nClustersTPC!=0)
1424 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1428 Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1429 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1430 SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1435 esdTrack->GetImpactParameters(b,bCov);
1436 if (bCov[0]<=0 || bCov[2]<=0) {
1437 AliDebug(1, "Estimated b resolution lower or equal zero!");
1438 bCov[0]=0; bCov[2]=0;
1441 Float_t dcaToVertexXY = b[0];
1442 Float_t dcaToVertexZ = b[1];
1443 Float_t dcaToVertex = -1;
1445 if (fCutDCAToVertex2D)
1446 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1448 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1452 Bool_t cuts[kNCuts];
1453 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1455 // track quality cuts
1456 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1458 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1460 if (nClustersTPC<fCutMinNClusterTPC)
1462 if (nClustersITS<fCutMinNClusterITS)
1464 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1466 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1468 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1470 if (fCutDCAToVertex2D && dcaToVertex > 1)
1472 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1474 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1477 //Require at least one SPD point + anything else in ITS
1478 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1482 for (Int_t i=0; i<kNCuts; i++)
1483 if (cuts[i]) {cut = kTRUE;}
1491 //__________________________________________________
1492 void AliEMCALRecoUtils::InitTrackCuts()
1494 //Intilize the track cut criteria
1495 //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1496 //Also you can customize the cuts using the setters
1499 SetMinNClustersTPC(70);
1500 SetMaxChi2PerClusterTPC(4);
1501 SetAcceptKinkDaughters(kFALSE);
1502 SetRequireTPCRefit(kTRUE);
1505 SetRequireITSRefit(kTRUE);
1506 SetMaxDCAToVertexZ(2);
1507 SetDCAToVertex2D(kFALSE);
1508 SetMaxChi2PerClusterITS(); //which by default sets the value to 1e10.
1509 SetMinNClustersITS();
1512 //___________________________________________________
1513 void AliEMCALRecoUtils::Print(const Option_t *) const
1517 printf("AliEMCALRecoUtils Settings: \n");
1518 printf("Misalignment shifts\n");
1519 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,
1520 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1521 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
1522 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1523 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
1525 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
1527 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1528 printf("Mass hypothesis = %2.3f[GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
1530 printf("Track cuts: \n");
1531 printf("AOD track selection mask: %d\n",fAODFilterMask);
1532 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1533 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1534 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1535 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1536 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1540 //_____________________________________________________________________
1541 void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1542 //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1543 //Do it only once and only if it is requested
1545 if(!fUseTimeCorrectionFactors) return;
1546 if(fTimeCorrectionFactorsSet) return;
1548 printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1550 AliEMCALCalibTimeDepCorrection *corr = new AliEMCALCalibTimeDepCorrection();
1551 corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1553 SwitchOnRecalibration();
1554 for(Int_t ism = 0; ism < 4; ism++){
1555 for(Int_t icol = 0; icol < 48; icol++){
1556 for(Int_t irow = 0; irow < 24; irow++){
1557 Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1558 Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1559 GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1560 //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
1561 // orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1562 // (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1566 fTimeCorrectionFactorsSet = kTRUE;