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>
40 #include <TObjArray.h>
43 #include "AliVCluster.h"
44 #include "AliVCaloCells.h"
45 #include "AliVEvent.h"
48 #include "AliESDEvent.h"
49 #include "AliAODEvent.h"
50 #include "AliESDtrack.h"
51 #include "AliAODTrack.h"
52 #include "AliExternalTrackParam.h"
53 #include "AliESDfriendTrack.h"
54 #include "AliTrackerBase.h"
57 #include "AliEMCALRecoUtils.h"
58 #include "AliEMCALGeometry.h"
59 #include "AliEMCALTrack.h"
60 #include "AliEMCALCalibTimeDepCorrection.h"
61 #include "AliEMCALPIDUtils.h"
63 ClassImp(AliEMCALRecoUtils)
65 //______________________________________________
66 AliEMCALRecoUtils::AliEMCALRecoUtils():
67 fParticleType(kPhoton), fPosAlgo(kUnchanged), fW0(4.),
68 fNonLinearityFunction(kNoCorrection), fNonLinearThreshold(30),
69 fSmearClusterEnergy(kFALSE), fRandom(),
70 fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
71 fUseTimeCorrectionFactors(kFALSE), fTimeCorrectionFactorsSet(kFALSE),
72 fRemoveBadChannels(kFALSE), fRecalDistToBadChannels(kFALSE), fEMCALBadChannelMap(),
73 fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
74 fRejectExoticCluster(kFALSE), fPIDUtils(), fAODFilterMask(32),
75 fMatchedTrackIndex(0x0), fMatchedClusterIndex(0x0),
76 fResidualEta(0x0), fResidualPhi(0x0), fCutEtaPhiSum(kTRUE), fCutEtaPhiSeparate(kFALSE),
77 fCutR(0.1), fCutEta(0.02), fCutPhi(0.04),
78 fMass(0.139), fStep(50),
79 fTrackCutsType(kTPCOnlyCut), fCutMinTrackPt(0), fCutMinNClusterTPC(-1),
80 fCutMinNClusterITS(-1), fCutMaxChi2PerClusterTPC(1e10), fCutMaxChi2PerClusterITS(1e10),
81 fCutRequireTPCRefit(kFALSE), fCutRequireITSRefit(kFALSE), fCutAcceptKinkDaughters(kFALSE),
82 fCutMaxDCAToVertexXY(1e10), fCutMaxDCAToVertexZ(1e10), fCutDCAToVertex2D(kFALSE)
86 // Initialize all constant values which have to be used
87 // during Reco algorithm execution
90 //Misalignment matrices
91 for(Int_t i = 0; i < 15 ; i++) {
92 fMisalTransShift[i] = 0.;
93 fMisalRotShift[i] = 0.;
97 for(Int_t i = 0; i < 7 ; i++) fNonLinearityParams[i] = 0.;
99 //For kBeamTestCorrected case, but default is no correction
100 fNonLinearityParams[0] = 0.99078;
101 fNonLinearityParams[1] = 0.161499;
102 fNonLinearityParams[2] = 0.655166;
103 fNonLinearityParams[3] = 0.134101;
104 fNonLinearityParams[4] = 163.282;
105 fNonLinearityParams[5] = 23.6904;
106 fNonLinearityParams[6] = 0.978;
108 //For kPi0GammaGamma case
109 //fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
110 //fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
111 //fNonLinearityParams[2] = 1.046;
113 //Cluster energy smearing
114 fSmearClusterEnergy = kFALSE;
115 fSmearClusterParam[0] = 0.07; // * sqrt E term
116 fSmearClusterParam[1] = 0.00; // * E term
117 fSmearClusterParam[2] = 0.00; // constant
120 fMatchedTrackIndex = new TArrayI();
121 fMatchedClusterIndex = new TArrayI();
122 fResidualPhi = new TArrayF();
123 fResidualEta = new TArrayF();
124 fPIDUtils = new AliEMCALPIDUtils();
129 //______________________________________________________________________
130 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
132 fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0),
133 fNonLinearityFunction(reco.fNonLinearityFunction), fNonLinearThreshold(reco.fNonLinearThreshold),
134 fSmearClusterEnergy(reco.fSmearClusterEnergy), fRandom(),
135 fRecalibration(reco.fRecalibration), fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
136 fUseTimeCorrectionFactors(reco.fUseTimeCorrectionFactors), fTimeCorrectionFactorsSet(reco.fTimeCorrectionFactorsSet),
137 fRemoveBadChannels(reco.fRemoveBadChannels), fRecalDistToBadChannels(reco.fRecalDistToBadChannels),
138 fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
139 fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder), fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
140 fRejectExoticCluster(reco.fRejectExoticCluster), fPIDUtils(reco.fPIDUtils),
141 fAODFilterMask(reco.fAODFilterMask),
142 fMatchedTrackIndex( reco.fMatchedTrackIndex? new TArrayI(*reco.fMatchedTrackIndex):0x0),
143 fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
144 fResidualEta( reco.fResidualEta? new TArrayF(*reco.fResidualEta):0x0),
145 fResidualPhi( reco.fResidualPhi? new TArrayF(*reco.fResidualPhi):0x0),
146 fCutEtaPhiSum(reco.fCutEtaPhiSum), fCutEtaPhiSeparate(reco.fCutEtaPhiSeparate),
147 fCutR(reco.fCutR), fCutEta(reco.fCutEta), fCutPhi(reco.fCutPhi),
148 fMass(reco.fMass), fStep(reco.fStep),
149 fTrackCutsType(reco.fTrackCutsType), fCutMinTrackPt(reco.fCutMinTrackPt),
150 fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS),
151 fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
152 fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
153 fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters), fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY),
154 fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ), fCutDCAToVertex2D(reco.fCutDCAToVertex2D)
158 for(Int_t i = 0; i < 15 ; i++) { fMisalRotShift[i] = reco.fMisalRotShift[i] ;
159 fMisalTransShift[i] = reco.fMisalTransShift[i] ; }
160 for(Int_t i = 0; i < 7 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
161 for(Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
166 //______________________________________________________________________
167 AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
169 //Assignment operator
171 if(this == &reco)return *this;
172 ((TNamed *)this)->operator=(reco);
174 for(Int_t i = 0; i < 15 ; i++) { fMisalTransShift[i] = reco.fMisalTransShift[i] ;
175 fMisalRotShift[i] = reco.fMisalRotShift[i] ; }
176 for(Int_t i = 0; i < 7 ; i++) { fNonLinearityParams[i] = reco.fNonLinearityParams[i] ; }
177 for(Int_t i = 0; i < 3 ; i++) { fSmearClusterParam[i] = reco.fSmearClusterParam[i] ; }
179 fParticleType = reco.fParticleType;
180 fPosAlgo = reco.fPosAlgo;
183 fNonLinearityFunction = reco.fNonLinearityFunction;
184 fNonLinearThreshold = reco.fNonLinearThreshold;
185 fSmearClusterEnergy = reco.fSmearClusterEnergy;
187 fRecalibration = reco.fRecalibration;
188 fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
189 fUseTimeCorrectionFactors = reco.fUseTimeCorrectionFactors;
190 fTimeCorrectionFactorsSet = reco.fTimeCorrectionFactorsSet;
192 fRemoveBadChannels = reco.fRemoveBadChannels;
193 fRecalDistToBadChannels = reco.fRecalDistToBadChannels;
194 fEMCALBadChannelMap = reco.fEMCALBadChannelMap;
196 fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
197 fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
198 fRejectExoticCluster = reco.fRejectExoticCluster;
200 fPIDUtils = reco.fPIDUtils;
202 fAODFilterMask = reco.fAODFilterMask;
204 fCutEtaPhiSum = reco.fCutEtaPhiSum;
205 fCutEtaPhiSeparate = reco.fCutEtaPhiSeparate;
207 fCutEta = reco.fCutEta;
208 fCutPhi = reco.fCutPhi;
211 fRejectExoticCluster = reco.fRejectExoticCluster;
213 fTrackCutsType = reco.fTrackCutsType;
214 fCutMinTrackPt = reco.fCutMinTrackPt;
215 fCutMinNClusterTPC = reco.fCutMinNClusterTPC;
216 fCutMinNClusterITS = reco.fCutMinNClusterITS;
217 fCutMaxChi2PerClusterTPC = reco.fCutMaxChi2PerClusterTPC;
218 fCutMaxChi2PerClusterITS = reco.fCutMaxChi2PerClusterITS;
219 fCutRequireTPCRefit = reco.fCutRequireTPCRefit;
220 fCutRequireITSRefit = reco.fCutRequireITSRefit;
221 fCutAcceptKinkDaughters = reco.fCutAcceptKinkDaughters;
222 fCutMaxDCAToVertexXY = reco.fCutMaxDCAToVertexXY;
223 fCutMaxDCAToVertexZ = reco.fCutMaxDCAToVertexZ;
224 fCutDCAToVertex2D = reco.fCutDCAToVertex2D;
226 if(reco.fResidualEta){
227 // assign or copy construct
229 *fResidualEta = *reco.fResidualEta;
231 else fResidualEta = new TArrayF(*reco.fResidualEta);
234 if(fResidualEta)delete fResidualEta;
238 if(reco.fResidualPhi){
239 // assign or copy construct
241 *fResidualPhi = *reco.fResidualPhi;
243 else fResidualPhi = new TArrayF(*reco.fResidualPhi);
246 if(fResidualPhi)delete fResidualPhi;
250 if(reco.fMatchedTrackIndex){
251 // assign or copy construct
252 if(fMatchedTrackIndex){
253 *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
255 else fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
258 if(fMatchedTrackIndex)delete fMatchedTrackIndex;
259 fMatchedTrackIndex = 0;
262 if(reco.fMatchedClusterIndex){
263 // assign or copy construct
264 if(fMatchedClusterIndex){
265 *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
267 else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
270 if(fMatchedClusterIndex)delete fMatchedClusterIndex;
271 fMatchedClusterIndex = 0;
278 //__________________________________________________
279 AliEMCALRecoUtils::~AliEMCALRecoUtils()
283 if(fEMCALRecalibrationFactors) {
284 fEMCALRecalibrationFactors->Clear();
285 delete fEMCALRecalibrationFactors;
288 if(fEMCALBadChannelMap) {
289 fEMCALBadChannelMap->Clear();
290 delete fEMCALBadChannelMap;
293 delete fMatchedTrackIndex ;
294 delete fMatchedClusterIndex ;
295 delete fResidualEta ;
296 delete fResidualPhi ;
300 //_______________________________________________________________
301 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells)
303 // Given the list of AbsId of the cluster, get the maximum cell and
304 // check if there are fNCellsFromBorder from the calorimeter border
307 AliInfo("Cluster pointer null!");
311 //If the distance to the border is 0 or negative just exit accept all clusters
312 if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
314 Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
315 Bool_t shared = kFALSE;
316 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
318 AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
319 absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
321 if(absIdMax==-1) return kFALSE;
323 //Check if the cell is close to the borders:
324 Bool_t okrow = kFALSE;
325 Bool_t okcol = kFALSE;
327 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
328 AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
334 if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
337 if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
341 if(!fNoEMCALBorderAtEta0){
342 if(ieta > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
346 if(ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
349 if(ieta < 48-fNCellsFromEMCALBorder) okcol = kTRUE;
353 AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
354 fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
356 if (okcol && okrow) {
357 //printf("Accept\n");
361 //printf("Reject\n");
362 AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
369 //_________________________________________________________________________________________________________
370 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, const Int_t nCells){
371 // Check that in the cluster cells, there is no bad channel of those stored
372 // in fEMCALBadChannelMap or fPHOSBadChannelMap
374 if(!fRemoveBadChannels) return kFALSE;
375 if(!fEMCALBadChannelMap) return kFALSE;
380 for(Int_t iCell = 0; iCell<nCells; iCell++){
382 //Get the column and row
383 Int_t iTower = -1, iIphi = -1, iIeta = -1;
384 geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
385 if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
386 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
387 if(GetEMCALChannelStatus(imod, icol, irow)){
388 AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
392 }// cell cluster loop
398 //_________________________________________________
399 Bool_t AliEMCALRecoUtils::IsExoticCluster(AliVCluster *cluster) const {
400 // Check if the cluster has high energy but small number of cells
401 // The criteria comes from Gustavo's study
405 AliInfo("Cluster pointer null!");
409 Int_t nc = cluster->GetNCells() ;
411 if ( nc > 8 ) return kFALSE ; // Good cluster, needed for 3x3 clusterizer
412 else if ( nc < 1 + cluster->E()/3. ) return kTRUE ; // Bad cluster
413 else return kFALSE ; // Good cluster
417 //__________________________________________________
418 Float_t AliEMCALRecoUtils::SmearClusterEnergy(AliVCluster* cluster) {
420 //In case of MC analysis, smear energy to match resolution/calibration in real data
423 AliInfo("Cluster pointer null!");
427 Float_t energy = cluster->E() ;
428 Float_t rdmEnergy = energy ;
429 if(fSmearClusterEnergy){
430 rdmEnergy = fRandom.Gaus(energy,fSmearClusterParam[0] * TMath::Sqrt(energy) +
431 fSmearClusterParam[1] * energy +
432 fSmearClusterParam[2] );
433 AliDebug(2, Form("Energy: original %f, smeared %f\n", energy, rdmEnergy));
440 //__________________________________________________
441 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
442 // Correct cluster energy from non linearity functions
445 AliInfo("Cluster pointer null!");
449 Float_t energy = cluster->E();
451 switch (fNonLinearityFunction) {
455 //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
456 //Double_t fNonLinearityParams[0] = 1.014;
457 //Double_t fNonLinearityParams[1] = -0.03329;
458 //Double_t fNonLinearityParams[2] = -0.3853;
459 //Double_t fNonLinearityParams[3] = 0.5423;
460 //Double_t fNonLinearityParams[4] = -0.4335;
461 energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
462 ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
463 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
469 //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
470 //Double_t fNonLinearityParams[0] = 1.04;
471 //Double_t fNonLinearityParams[1] = -0.1445;
472 //Double_t fNonLinearityParams[2] = 1.046;
473 energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
477 case kPi0GammaConversion:
479 //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
480 //fNonLinearityParams[0] = 0.139393/0.1349766;
481 //fNonLinearityParams[1] = 0.0566186;
482 //fNonLinearityParams[2] = 0.982133;
483 energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
490 //From beam test, Alexei's results, for different ZS thresholds
491 // th=30 MeV; th = 45 MeV; th = 75 MeV
492 //fNonLinearityParams[0] = 1.007; 1.003; 1.002
493 //fNonLinearityParams[1] = 0.894; 0.719; 0.797
494 //fNonLinearityParams[2] = 0.246; 0.334; 0.358
495 //Rescale the param[0] with 1.03
496 energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
501 case kBeamTestCorrected:
503 //From beam test, corrected for material between beam and EMCAL
504 //fNonLinearityParams[0] = 0.99078
505 //fNonLinearityParams[1] = 0.161499;
506 //fNonLinearityParams[2] = 0.655166;
507 //fNonLinearityParams[3] = 0.134101;
508 //fNonLinearityParams[4] = 163.282;
509 //fNonLinearityParams[5] = 23.6904;
510 //fNonLinearityParams[6] = 0.978;
511 energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
517 AliDebug(2,"No correction on the energy\n");
525 //__________________________________________________
526 void AliEMCALRecoUtils::InitNonLinearityParam()
528 //Initialising Non Linearity Parameters
530 if(fNonLinearityFunction == kPi0MC)
532 fNonLinearityParams[0] = 1.014;
533 fNonLinearityParams[1] = -0.03329;
534 fNonLinearityParams[2] = -0.3853;
535 fNonLinearityParams[3] = 0.5423;
536 fNonLinearityParams[4] = -0.4335;
539 if(fNonLinearityFunction == kPi0GammaGamma)
541 fNonLinearityParams[0] = 1.04;
542 fNonLinearityParams[1] = -0.1445;
543 fNonLinearityParams[2] = 1.046;
546 if(fNonLinearityFunction == kPi0GammaConversion)
548 fNonLinearityParams[0] = 0.139393;
549 fNonLinearityParams[1] = 0.0566186;
550 fNonLinearityParams[2] = 0.982133;
553 if(fNonLinearityFunction == kBeamTest)
555 if(fNonLinearThreshold == 30)
557 fNonLinearityParams[0] = 1.007;
558 fNonLinearityParams[1] = 0.894;
559 fNonLinearityParams[2] = 0.246;
561 if(fNonLinearThreshold == 45)
563 fNonLinearityParams[0] = 1.003;
564 fNonLinearityParams[1] = 0.719;
565 fNonLinearityParams[2] = 0.334;
567 if(fNonLinearThreshold == 75)
569 fNonLinearityParams[0] = 1.002;
570 fNonLinearityParams[1] = 0.797;
571 fNonLinearityParams[2] = 0.358;
575 if(fNonLinearityFunction == kBeamTestCorrected)
577 fNonLinearityParams[0] = 0.99078;
578 fNonLinearityParams[1] = 0.161499;
579 fNonLinearityParams[2] = 0.655166;
580 fNonLinearityParams[3] = 0.134101;
581 fNonLinearityParams[4] = 163.282;
582 fNonLinearityParams[5] = 23.6904;
583 fNonLinearityParams[6] = 0.978;
587 //__________________________________________________
588 Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
590 //Calculate shower depth for a given cluster energy and particle type
600 depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
604 depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
611 gGeoManager->cd("ALIC_1/XEN1_1");
612 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
613 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
615 TGeoVolume *geoSMVol = geoSM->GetVolume();
616 TGeoShape *geoSMShape = geoSMVol->GetShape();
617 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
618 if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
619 else AliFatal("Null GEANT box");
620 }else AliFatal("NULL GEANT node matrix");
623 depth = x0 * (TMath::Log(energy*1000 / ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
629 depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
636 //__________________________________________________
637 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
638 Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
640 //For a given CaloCluster gets the absId of the cell
641 //with maximum energy deposit.
644 Double_t eCell = -1.;
645 Float_t fraction = 1.;
646 Float_t recalFactor = 1.;
647 Int_t cellAbsId = -1 ;
655 AliInfo("Cluster pointer null!");
656 absId=-1; iSupMod0=-1, ieta = -1; iphi = -1; shared = -1;
660 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
661 cellAbsId = clu->GetCellAbsId(iDig);
662 fraction = clu->GetCellAmplitudeFraction(iDig);
663 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
664 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
665 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
666 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
667 if(iDig==0) iSupMod0=iSupMod;
668 else if(iSupMod0!=iSupMod) {
670 //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
672 if(IsRecalibrationOn()) {
673 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
675 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
676 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
680 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
684 //Get from the absid the supermodule, tower and eta/phi numbers
685 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
686 //Gives SuperModule and Tower numbers
687 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
688 iIphi, iIeta,iphi,ieta);
689 //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
690 //printf("Max end---\n");
694 //________________________________________________________________
695 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
696 //Init EMCAL recalibration factors
697 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
698 //In order to avoid rewriting the same histograms
699 Bool_t oldStatus = TH1::AddDirectoryStatus();
700 TH1::AddDirectory(kFALSE);
702 fEMCALRecalibrationFactors = new TObjArray(10);
703 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));
704 //Init the histograms with 1
705 for (Int_t sm = 0; sm < 10; sm++) {
706 for (Int_t i = 0; i < 48; i++) {
707 for (Int_t j = 0; j < 24; j++) {
708 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
712 fEMCALRecalibrationFactors->SetOwner(kTRUE);
713 fEMCALRecalibrationFactors->Compress();
715 //In order to avoid rewriting the same histograms
716 TH1::AddDirectory(oldStatus);
720 //________________________________________________________________
721 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
722 //Init EMCAL bad channels map
723 AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
724 //In order to avoid rewriting the same histograms
725 Bool_t oldStatus = TH1::AddDirectoryStatus();
726 TH1::AddDirectory(kFALSE);
728 fEMCALBadChannelMap = new TObjArray(10);
729 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
730 for (int i = 0; i < 10; i++) {
731 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
736 fEMCALBadChannelMap->SetOwner(kTRUE);
737 fEMCALBadChannelMap->Compress();
739 //In order to avoid rewriting the same histograms
740 TH1::AddDirectory(oldStatus);
743 //________________________________________________________________
744 void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
745 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
748 AliInfo("Cluster pointer null!");
752 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
753 UShort_t * index = cluster->GetCellsAbsId() ;
754 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
755 Int_t ncells = cluster->GetNCells();
757 //Initialize some used variables
760 Int_t icol = -1, irow = -1, imod=1;
761 Float_t factor = 1, frac = 0;
763 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
764 for(Int_t icell = 0; icell < ncells; icell++){
765 absId = index[icell];
766 frac = fraction[icell];
767 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
768 Int_t iTower = -1, iIphi = -1, iIeta = -1;
769 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
770 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
771 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
772 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
773 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
774 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
776 energy += cells->GetCellAmplitude(absId)*factor*frac;
780 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
782 cluster->SetE(energy);
787 //__________________________________________________
788 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
790 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
793 AliInfo("Cluster pointer null!");
797 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
798 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
799 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
803 //__________________________________________________
804 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
806 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
807 // The algorithm is a copy of what is done in AliEMCALRecPoint
810 Float_t fraction = 1.;
811 Float_t recalFactor = 1.;
814 Int_t iTower = -1, iIphi = -1, iIeta = -1;
815 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
816 Float_t weight = 0., totalWeight=0.;
817 Float_t newPos[3] = {0,0,0};
818 Double_t pLocal[3], pGlobal[3];
819 Bool_t shared = kFALSE;
821 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
822 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
823 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
825 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
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,iSM,iTower,iIphi,iIeta);
832 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
834 if(IsRecalibrationOn()) {
835 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
837 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
839 weight = GetCellWeight(eCell,clEnergy);
840 //printf("cell energy %f, weight %f\n",eCell,weight);
841 totalWeight += weight;
842 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
843 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
844 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
845 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
847 for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
852 for(int i=0; i<3; i++ ) newPos[i] /= totalWeight;
855 //Float_t pos[]={0,0,0};
856 //clu->GetPosition(pos);
857 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
858 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
860 if(iSupModMax > 1) {//sector 1
861 newPos[0] +=fMisalTransShift[3];//-=3.093;
862 newPos[1] +=fMisalTransShift[4];//+=6.82;
863 newPos[2] +=fMisalTransShift[5];//+=1.635;
864 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
868 newPos[0] +=fMisalTransShift[0];//+=1.134;
869 newPos[1] +=fMisalTransShift[1];//+=8.2;
870 newPos[2] +=fMisalTransShift[2];//+=1.197;
871 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
874 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
876 clu->SetPosition(newPos);
880 //__________________________________________________
881 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
883 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
884 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
887 Float_t fraction = 1.;
888 Float_t recalFactor = 1.;
892 Int_t iIphi = -1, iIeta = -1;
893 Int_t iSupMod = -1, iSupModMax = -1;
894 Int_t iphi = -1, ieta =-1;
895 Bool_t shared = kFALSE;
897 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
898 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
899 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
901 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
902 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
903 Int_t startingSM = -1;
905 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
906 absId = clu->GetCellAbsId(iDig);
907 fraction = clu->GetCellAmplitudeFraction(iDig);
908 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
909 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
910 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
912 if (iDig==0) startingSM = iSupMod;
913 else if(iSupMod != startingSM) areInSameSM = kFALSE;
915 eCell = cells->GetCellAmplitude(absId);
917 if(IsRecalibrationOn()) {
918 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
920 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
922 weight = GetCellWeight(eCell,clEnergy);
923 if(weight < 0) weight = 0;
924 totalWeight += weight;
925 weightedCol += ieta*weight;
926 weightedRow += iphi*weight;
928 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
932 Float_t xyzNew[]={0.,0.,0.};
933 if(areInSameSM == kTRUE) {
934 //printf("In Same SM\n");
935 weightedCol = weightedCol/totalWeight;
936 weightedRow = weightedRow/totalWeight;
937 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
940 //printf("In Different SM\n");
941 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
944 clu->SetPosition(xyzNew);
948 //____________________________________________________________________________
949 void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
951 //re-evaluate distance to bad channel with updated bad map
953 if(!fRecalDistToBadChannels) return;
956 AliInfo("Cluster pointer null!");
960 //Get channels map of the supermodule where the cluster is.
961 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
962 Bool_t shared = kFALSE;
963 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
964 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
967 Float_t minDist = 10000.;
970 //Loop on tower status map
971 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
972 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
973 //Check if tower is bad.
974 if(hMap->GetBinContent(icol,irow)==0) continue;
975 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
976 // iSupMod,icol, irow, icolM,irowM);
978 dRrow=TMath::Abs(irowM-irow);
979 dRcol=TMath::Abs(icolM-icol);
980 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
982 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
989 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
994 //The only possible combinations are (0,1), (2,3) ... (8,9)
995 if(iSupMod%2) iSupMod2 = iSupMod-1;
996 else iSupMod2 = iSupMod+1;
997 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
999 //Loop on tower status map of second super module
1000 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
1001 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
1002 //Check if tower is bad.
1003 if(hMap2->GetBinContent(icol,irow)==0) continue;
1004 //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",
1005 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
1007 dRrow=TMath::Abs(irow-irowM);
1010 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
1013 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
1016 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
1017 if(dist < minDist) minDist = dist;
1022 }// shared cluster in 2 SuperModules
1024 AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
1025 cluster->SetDistanceToBadChannel(minDist);
1029 //____________________________________________________________________________
1030 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
1032 //re-evaluate identification parameters with bayesian
1035 AliInfo("Cluster pointer null!");
1039 if ( cluster->GetM02() != 0)
1040 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
1042 Float_t pidlist[AliPID::kSPECIESN+1];
1043 for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
1045 cluster->SetPID(pidlist);
1049 //____________________________________________________________________________
1050 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
1052 // Calculates new center of gravity in the local EMCAL-module coordinates
1053 // and tranfers into global ALICE coordinates
1054 // Calculates Dispersion and main axis
1057 AliInfo("Cluster pointer null!");
1063 Double_t eCell = 0.;
1064 Float_t fraction = 1.;
1065 Float_t recalFactor = 1.;
1073 Double_t etai = -1.;
1074 Double_t phii = -1.;
1081 Double_t xmean = 0.;
1082 Double_t zmean = 0.;
1085 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
1087 //Get from the absid the supermodule, tower and eta/phi numbers
1088 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1089 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1091 //Get the cell energy, if recalibration is on, apply factors
1092 fraction = cluster->GetCellAmplitudeFraction(iDigit);
1093 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1094 if(IsRecalibrationOn()) {
1095 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1097 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1099 if(cluster->E() > 0 && eCell > 0){
1101 w = GetCellWeight(eCell,cluster->E());
1103 etai=(Double_t)ieta;
1104 phii=(Double_t)iphi;
1109 dxx += w * etai * etai ;
1111 dzz += w * phii * phii ;
1113 dxz += w * etai * phii ;
1117 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1120 //Normalize to the weight
1126 AliError(Form("Wrong weight %f\n", wtot));
1128 //Calculate dispersion
1129 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
1131 //Get from the absid the supermodule, tower and eta/phi numbers
1132 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
1133 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
1135 //Get the cell energy, if recalibration is on, apply factors
1136 fraction = cluster->GetCellAmplitudeFraction(iDigit);
1137 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1138 if(IsRecalibrationOn()) {
1139 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1141 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
1143 if(cluster->E() > 0 && eCell > 0){
1145 w = GetCellWeight(eCell,cluster->E());
1147 etai=(Double_t)ieta;
1148 phii=(Double_t)iphi;
1149 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
1152 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
1155 //Normalize to the weigth and set shower shape parameters
1156 if (wtot > 0 && nstat > 1) {
1161 dxx -= xmean * xmean ;
1162 dzz -= zmean * zmean ;
1163 dxz -= xmean * zmean ;
1164 cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1165 cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1169 cluster->SetM20(0.) ;
1170 cluster->SetM02(0.) ;
1174 cluster->SetDispersion(TMath::Sqrt(d)) ;
1176 cluster->SetDispersion(0) ;
1179 //____________________________________________________________________________
1180 void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr, AliEMCALGeometry *geom)
1182 //This function should be called before the cluster loop
1183 //Before call this function, please recalculate the cluster positions
1184 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
1185 //Store matched cluster indexes and residuals
1187 fMatchedTrackIndex->Reset();
1188 fMatchedClusterIndex->Reset();
1189 fResidualPhi->Reset();
1190 fResidualEta->Reset();
1192 fMatchedTrackIndex->Set(500);
1193 fMatchedClusterIndex->Set(500);
1194 fResidualPhi->Set(500);
1195 fResidualEta->Set(500);
1197 AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1198 AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
1202 for (Int_t i=0; i<21;i++) cv[i]=0;
1203 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1205 AliExternalTrackParam *trackParam = 0;
1207 //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
1210 AliESDtrack *esdTrack = esdevent->GetTrack(itr);
1211 if(!esdTrack || !IsAccepted(esdTrack)) continue;
1212 if(esdTrack->Pt()<fCutMinTrackPt) continue;
1213 trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
1216 //If the input event is AOD, the starting point for extrapolation is at vertex
1217 //AOD tracks are selected according to its bit.
1220 AliAODTrack *aodTrack = aodevent->GetTrack(itr);
1221 if(!aodTrack) continue;
1222 if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
1223 if(aodTrack->Pt()<fCutMinTrackPt) continue;
1224 Double_t pos[3],mom[3];
1225 aodTrack->GetXYZ(pos);
1226 aodTrack->GetPxPyPz(mom);
1227 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()));
1228 trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1231 //Return if the input data is not "AOD" or "ESD"
1234 printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
1238 if(!trackParam) continue;
1240 Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
1242 if(!clusterArr){// get clusters from event
1243 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1245 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1246 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1247 AliExternalTrackParam trkPamTmp(*trackParam);//Retrieve the starting point every time before the extrapolation
1248 Float_t tmpEta=-999, tmpPhi=-999;
1249 if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1252 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1261 else if(fCutEtaPhiSeparate)
1263 if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1272 printf("Error: please specify your cut criteria\n");
1273 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1274 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1275 if(aodevent && trackParam) delete trackParam;
1280 else { // external cluster array, not from ESD event
1281 for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1283 AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
1285 AliInfo("Cluster not found!!!");
1288 if(!cluster->IsEMCAL()) continue;
1289 AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
1290 Float_t tmpEta=-999, tmpPhi=-999;
1291 if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1294 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1303 else if(fCutEtaPhiSeparate)
1305 if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1314 printf("Error: please specify your cut criteria\n");
1315 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1316 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1317 if(aodevent && trackParam) delete trackParam;
1321 }// external list of clusters
1325 fMatchedTrackIndex ->AddAt(itr,matched);
1326 fMatchedClusterIndex->AddAt(index,matched);
1327 fResidualEta ->AddAt(dEtaMax,matched);
1328 fResidualPhi ->AddAt(dPhiMax,matched);
1331 if(aodevent && trackParam) delete trackParam;
1334 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1336 fMatchedTrackIndex ->Set(matched);
1337 fMatchedClusterIndex->Set(matched);
1338 fResidualPhi ->Set(matched);
1339 fResidualEta ->Set(matched);
1342 //________________________________________________________________________________
1343 Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom)
1346 // This function returns the index of matched cluster to input track
1347 // Returns -1 if no match is found
1350 Float_t dRMax = fCutR, dEtaMax = fCutEta, dPhiMax = fCutPhi;
1353 AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
1355 if(!trackParam) return index;
1356 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1358 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1359 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1360 AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
1361 Float_t tmpEta=-999, tmpPhi=-999;
1362 if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
1365 Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
1374 else if(fCutEtaPhiSeparate)
1376 if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
1385 printf("Error: please specify your cut criteria\n");
1386 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1387 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1394 //________________________________________________________________________________
1395 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpEta, Float_t &tmpPhi)
1398 //Return the residual by extrapolating a track to a cluster
1400 if(!trkParam || !cluster) return kFALSE;
1403 cluster->GetPosition(clsPos); //Has been recalculated
1404 TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
1405 Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
1406 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
1407 if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kTRUE, 0.99, -1)) return kFALSE;
1408 trkParam->GetXYZ(trkPos); //Get the extrapolated global position
1410 TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
1411 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
1413 // track cluster matching
1414 tmpPhi = clsPosVec.DeltaPhi(trkPosVec); // tmpPhi is between -pi and pi
1415 tmpEta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
1420 //________________________________________________________________________________
1421 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
1423 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1424 //Get the residuals dEta and dPhi for this cluster to the closest track
1425 //Works with ESDs and AODs
1427 if( FindMatchedPosForCluster(clsIndex) >= 999 )
1429 AliDebug(2,"No matched tracks found!\n");
1434 dEta = fResidualEta->At(FindMatchedPosForCluster(clsIndex));
1435 dPhi = fResidualPhi->At(FindMatchedPosForCluster(clsIndex));
1437 //________________________________________________________________________________
1438 void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dEta, Float_t &dPhi)
1440 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1441 //Get the residuals dEta and dPhi for this track to the closest cluster
1442 //Works with ESDs and AODs
1444 if( FindMatchedPosForTrack(trkIndex) >= 999 )
1446 AliDebug(2,"No matched cluster found!\n");
1451 dEta = fResidualEta->At(FindMatchedPosForTrack(trkIndex));
1452 dPhi = fResidualPhi->At(FindMatchedPosForTrack(trkIndex));
1455 //__________________________________________________________
1456 Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
1458 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1459 //Get the index of matched track to this cluster
1460 //Works with ESDs and AODs
1462 if(IsClusterMatched(clsIndex))
1463 return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
1468 //__________________________________________________________
1469 Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
1471 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1472 //Get the index of matched cluster to this track
1473 //Works with ESDs and AODs
1475 if(IsTrackMatched(trkIndex))
1476 return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
1481 //__________________________________________________
1482 Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex) const
1484 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1485 //Returns if the cluster has a match
1486 if(FindMatchedPosForCluster(clsIndex) < 999)
1492 //__________________________________________________
1493 Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex) const
1495 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1496 //Returns if the track has a match
1497 if(FindMatchedPosForTrack(trkIndex) < 999)
1503 //__________________________________________________________
1504 UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
1506 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1507 //Returns the position of the match in the fMatchedClusterIndex array
1508 Float_t tmpR = fCutR;
1511 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
1513 if(fMatchedClusterIndex->At(i)==clsIndex)
1515 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1520 AliDebug(3,Form("Matched cluster index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedClusterIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1527 //__________________________________________________________
1528 UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
1530 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1531 //Returns the position of the match in the fMatchedTrackIndex array
1532 Float_t tmpR = fCutR;
1535 for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
1537 if(fMatchedTrackIndex->At(i)==trkIndex)
1539 Float_t r = TMath::Sqrt(fResidualEta->At(i)*fResidualEta->At(i) + fResidualPhi->At(i)*fResidualPhi->At(i));
1544 AliDebug(3,Form("Matched track index: index: %d, dEta: %2.4f, dPhi: %2.4f.\n",fMatchedTrackIndex->At(i),fResidualEta->At(i),fResidualPhi->At(i)));
1551 //__________________________________________________________
1552 Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells)
1554 // check if the cluster survives some quality cut
1557 Bool_t isGood=kTRUE;
1558 if(!cluster || !cluster->IsEMCAL()) return kFALSE;
1559 if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
1560 if(!CheckCellFiducialRegion(geom,cluster,cells)) return kFALSE;
1561 if(fRejectExoticCluster && IsExoticCluster(cluster)) return kFALSE;
1566 //__________________________________________________________
1567 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1569 // Given a esd track, return whether the track survive all the cuts
1571 // The different quality parameter are first
1572 // retrieved from the track. then it is found out what cuts the
1573 // track did not survive and finally the cuts are imposed.
1575 UInt_t status = esdTrack->GetStatus();
1577 Int_t nClustersITS = esdTrack->GetITSclusters(0);
1578 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1580 Float_t chi2PerClusterITS = -1;
1581 Float_t chi2PerClusterTPC = -1;
1582 if (nClustersITS!=0)
1583 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1584 if (nClustersTPC!=0)
1585 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1589 Float_t maxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1590 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1591 SetMaxDCAToVertexXY(maxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1596 esdTrack->GetImpactParameters(b,bCov);
1597 if (bCov[0]<=0 || bCov[2]<=0) {
1598 AliDebug(1, "Estimated b resolution lower or equal zero!");
1599 bCov[0]=0; bCov[2]=0;
1602 Float_t dcaToVertexXY = b[0];
1603 Float_t dcaToVertexZ = b[1];
1604 Float_t dcaToVertex = -1;
1606 if (fCutDCAToVertex2D)
1607 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1609 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1613 Bool_t cuts[kNCuts];
1614 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1616 // track quality cuts
1617 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1619 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1621 if (nClustersTPC<fCutMinNClusterTPC)
1623 if (nClustersITS<fCutMinNClusterITS)
1625 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1627 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1629 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1631 if (fCutDCAToVertex2D && dcaToVertex > 1)
1633 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1635 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1638 //Require at least one SPD point + anything else in ITS
1639 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1643 for (Int_t i=0; i<kNCuts; i++)
1644 if (cuts[i]) {cut = kTRUE;}
1652 //__________________________________________________
1653 void AliEMCALRecoUtils::InitTrackCuts()
1655 //Intilize the track cut criteria
1656 //By default these cuts are set according to AliESDtrackCuts::GetStandardTPCOnlyTrackCuts()
1657 //Also you can customize the cuts using the setters
1659 switch (fTrackCutsType)
1663 AliInfo(Form("Track cuts for matching: GetStandardTPCOnlyTrackCuts()"));
1665 SetMinNClustersTPC(70);
1666 SetMaxChi2PerClusterTPC(4);
1667 SetAcceptKinkDaughters(kFALSE);
1668 SetRequireTPCRefit(kFALSE);
1671 SetRequireITSRefit(kFALSE);
1672 SetMaxDCAToVertexZ(3.2);
1673 SetMaxDCAToVertexXY(2.4);
1674 SetDCAToVertex2D(kTRUE);
1681 AliInfo(Form("Track cuts for matching: GetStandardITSTPCTrackCuts2010(kTURE)"));
1683 SetMinNClustersTPC(70);
1684 SetMaxChi2PerClusterTPC(4);
1685 SetAcceptKinkDaughters(kFALSE);
1686 SetRequireTPCRefit(kTRUE);
1689 SetRequireITSRefit(kTRUE);
1690 SetMaxDCAToVertexZ(2);
1691 SetMaxDCAToVertexXY();
1692 SetDCAToVertex2D(kFALSE);
1699 AliInfo(Form("Track cuts for matching: Loose cut w/o DCA cut"));
1700 SetMinNClustersTPC(50);
1701 SetAcceptKinkDaughters(kFALSE);
1708 //___________________________________________________
1709 void AliEMCALRecoUtils::Print(const Option_t *) const
1713 printf("AliEMCALRecoUtils Settings: \n");
1714 printf("Misalignment shifts\n");
1715 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,
1716 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1717 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
1718 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1719 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
1721 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
1723 printf("Matching criteria: ");
1726 printf("sqrt(dEta^2+dPhi^2)<%2.2f\n",fCutR);
1728 else if(fCutEtaPhiSeparate)
1730 printf("dEta<%2.2f, dPhi<%2.2f\n",fCutEta,fCutPhi);
1735 printf("please specify your cut criteria\n");
1736 printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
1737 printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
1740 printf("Mass hypothesis = %2.3f [GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
1742 printf("Track cuts: \n");
1743 printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
1744 printf("AOD track selection mask: %d\n",fAODFilterMask);
1745 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1746 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1747 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1748 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1749 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1753 //_____________________________________________________________________
1754 void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1755 //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1756 //Do it only once and only if it is requested
1758 if(!fUseTimeCorrectionFactors) return;
1759 if(fTimeCorrectionFactorsSet) return;
1761 printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1763 AliEMCALCalibTimeDepCorrection *corr = new AliEMCALCalibTimeDepCorrection();
1764 corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1766 SwitchOnRecalibration();
1767 for(Int_t ism = 0; ism < 4; ism++){
1768 for(Int_t icol = 0; icol < 48; icol++){
1769 for(Int_t irow = 0; irow < 24; irow++){
1770 Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1771 Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1772 GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1773 //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
1774 // orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1775 // (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1779 fTimeCorrectionFactorsSet = kTRUE;