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.),
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),
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 fRecalibration = reco.fRecalibration;
167 fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
168 fRemoveBadChannels = reco.fRemoveBadChannels;
169 fRecalDistToBadChannels = reco.fRecalDistToBadChannels;
170 fEMCALBadChannelMap = reco.fEMCALBadChannelMap;
171 fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
172 fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
175 for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
176 for(Int_t i = 0; i < 7 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
178 fAODFilterMask = reco.fAODFilterMask;
184 fRejectExoticCluster = reco.fRejectExoticCluster;
186 fCutMinNClusterTPC = reco.fCutMinNClusterTPC;
187 fCutMinNClusterITS = reco.fCutMinNClusterITS;
188 fCutMaxChi2PerClusterTPC = reco.fCutMaxChi2PerClusterTPC;
189 fCutMaxChi2PerClusterITS = reco.fCutMaxChi2PerClusterITS;
190 fCutRequireTPCRefit = reco.fCutRequireTPCRefit;
191 fCutRequireITSRefit = reco.fCutRequireITSRefit;
192 fCutAcceptKinkDaughters = reco.fCutAcceptKinkDaughters;
193 fCutMaxDCAToVertexXY = reco.fCutMaxDCAToVertexXY;
194 fCutMaxDCAToVertexZ = reco.fCutMaxDCAToVertexZ;
195 fCutDCAToVertex2D = reco.fCutDCAToVertex2D;
197 fPIDUtils = reco.fPIDUtils;
199 fUseTimeCorrectionFactors = reco.fUseTimeCorrectionFactors;
200 fTimeCorrectionFactorsSet = reco.fTimeCorrectionFactorsSet;
204 // assign or copy construct
206 *fResidualR = *reco.fResidualR;
208 else fResidualR = new TArrayF(*reco.fResidualR);
211 if(fResidualR)delete fResidualR;
216 // assign or copy construct
218 *fResidualZ = *reco.fResidualZ;
220 else fResidualZ = new TArrayF(*reco.fResidualZ);
223 if(fResidualZ)delete fResidualZ;
227 if(reco.fMatchedTrackIndex){
228 // assign or copy construct
229 if(fMatchedTrackIndex){
230 *fMatchedTrackIndex = *reco.fMatchedTrackIndex;
232 else fMatchedTrackIndex = new TArrayI(*reco.fMatchedTrackIndex);
235 if(fMatchedTrackIndex)delete fMatchedTrackIndex;
236 fMatchedTrackIndex = 0;
239 if(reco.fMatchedClusterIndex){
240 // assign or copy construct
241 if(fMatchedClusterIndex){
242 *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
244 else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
247 if(fMatchedClusterIndex)delete fMatchedClusterIndex;
248 fMatchedClusterIndex = 0;
256 //__________________________________________________
257 AliEMCALRecoUtils::~AliEMCALRecoUtils()
261 if(fEMCALRecalibrationFactors) {
262 fEMCALRecalibrationFactors->Clear();
263 delete fEMCALRecalibrationFactors;
266 if(fEMCALBadChannelMap) {
267 fEMCALBadChannelMap->Clear();
268 delete fEMCALBadChannelMap;
271 if(fMatchedTrackIndex) {delete fMatchedTrackIndex; fMatchedTrackIndex=0;}
272 if(fMatchedClusterIndex) {delete fMatchedClusterIndex; fMatchedClusterIndex=0;}
273 if(fResidualR) {delete fResidualR; fResidualR=0;}
274 if(fResidualZ) {delete fResidualZ; fResidualZ=0;}
278 //_______________________________________________________________
279 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells)
281 // Given the list of AbsId of the cluster, get the maximum cell and
282 // check if there are fNCellsFromBorder from the calorimeter border
284 //If the distance to the border is 0 or negative just exit accept all clusters
285 if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
287 Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
288 Bool_t shared = kFALSE;
289 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi, shared);
291 AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
292 absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
294 if(absIdMax==-1) return kFALSE;
296 //Check if the cell is close to the borders:
297 Bool_t okrow = kFALSE;
298 Bool_t okcol = kFALSE;
300 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
301 AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
307 if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
310 if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
314 if(!fNoEMCALBorderAtEta0){
315 if(ieta > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
319 if(ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
322 if(ieta < 48-fNCellsFromEMCALBorder) okcol = kTRUE;
326 AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
327 fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
329 if (okcol && okrow) {
330 //printf("Accept\n");
334 //printf("Reject\n");
335 AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
342 //_________________________________________________________________________________________________________
343 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells){
344 // Check that in the cluster cells, there is no bad channel of those stored
345 // in fEMCALBadChannelMap or fPHOSBadChannelMap
347 if(!fRemoveBadChannels) return kFALSE;
348 if(!fEMCALBadChannelMap) return kFALSE;
353 for(Int_t iCell = 0; iCell<nCells; iCell++){
355 //Get the column and row
356 Int_t iTower = -1, iIphi = -1, iIeta = -1;
357 geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
358 if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
359 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
360 if(GetEMCALChannelStatus(imod, icol, irow)){
361 AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
365 }// cell cluster loop
371 //_________________________________________________
372 Bool_t AliEMCALRecoUtils::IsExoticCluster(AliVCluster *cluster){
373 // Check if the cluster has high energy but small number of cells
374 // The criteria comes from Gustavo's study
377 if(cluster->GetNCells()<(1+cluster->E()/3.))
383 //__________________________________________________
384 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
385 // Correct cluster energy from non linearity functions
386 Float_t energy = cluster->E();
388 switch (fNonLinearityFunction) {
392 //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
393 //Double_t fNonLinearityParams[0] = 1.014;
394 //Double_t fNonLinearityParams[1] = -0.03329;
395 //Double_t fNonLinearityParams[2] = -0.3853;
396 //Double_t fNonLinearityParams[3] = 0.5423;
397 //Double_t fNonLinearityParams[4] = -0.4335;
398 energy *= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
399 ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
400 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
406 //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
407 //Double_t fNonLinearityParams[0] = 1.04;
408 //Double_t fNonLinearityParams[1] = -0.1445;
409 //Double_t fNonLinearityParams[2] = 1.046;
410 energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
414 case kPi0GammaConversion:
416 //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
417 //fNonLinearityParams[0] = 0.139393/0.1349766;
418 //fNonLinearityParams[1] = 0.0566186;
419 //fNonLinearityParams[2] = 0.982133;
420 energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
427 //From beam test, Alexei's results, for different ZS thresholds
428 // th=30 MeV; th = 45 MeV; th = 75 MeV
429 //fNonLinearityParams[0] = 1.007; 1.003; 1.002
430 //fNonLinearityParams[1] = 0.894; 0.719; 0.797
431 //fNonLinearityParams[2] = 0.246; 0.334; 0.358
432 //Rescale the param[0] with 1.03
433 energy /= fNonLinearityParams[0]/(1+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]));
438 case kBeamTestCorrected:
440 //From beam test, corrected for material between beam and EMCAL
441 //fNonLinearityParams[0] = 0.99078
442 //fNonLinearityParams[1] = 0.161499;
443 //fNonLinearityParams[2] = 0.655166;
444 //fNonLinearityParams[3] = 0.134101;
445 //fNonLinearityParams[4] = 163.282;
446 //fNonLinearityParams[5] = 23.6904;
447 //fNonLinearityParams[6] = 0.978;
448 energy *= fNonLinearityParams[6]/(fNonLinearityParams[0]*(1./(1.+fNonLinearityParams[1]*exp(-energy/fNonLinearityParams[2]))*1./(1.+fNonLinearityParams[3]*exp((energy-fNonLinearityParams[4])/fNonLinearityParams[5]))));
454 AliDebug(2,"No correction on the energy\n");
462 //__________________________________________________
463 Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
465 //Calculate shower depth for a given cluster energy and particle type
475 depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
479 depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
486 gGeoManager->cd("ALIC_1/XEN1_1");
487 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
488 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
490 TGeoVolume *geoSMVol = geoSM->GetVolume();
491 TGeoShape *geoSMShape = geoSMVol->GetShape();
492 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
493 if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
494 else AliFatal("Null GEANT box");
495 }else AliFatal("NULL GEANT node matrix");
498 depth = x0 * (TMath::Log(energy*1000 / ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
504 depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
511 //__________________________________________________
512 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu,
513 Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi, Bool_t &shared)
515 //For a given CaloCluster gets the absId of the cell
516 //with maximum energy deposit.
519 Double_t eCell = -1.;
520 Float_t fraction = 1.;
521 Float_t recalFactor = 1.;
522 Int_t cellAbsId = -1 ;
528 //printf("---Max?\n");
529 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
530 cellAbsId = clu->GetCellAbsId(iDig);
531 fraction = clu->GetCellAmplitudeFraction(iDig);
532 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
533 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
534 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
535 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
536 if(iDig==0) iSupMod0=iSupMod;
537 else if(iSupMod0!=iSupMod) {
539 //printf("AliEMCALRecoUtils::GetMaxEnergyCell() - SHARED CLUSTER\n");
541 if(IsRecalibrationOn()) {
542 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
544 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
545 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
549 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
553 //Get from the absid the supermodule, tower and eta/phi numbers
554 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
555 //Gives SuperModule and Tower numbers
556 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
557 iIphi, iIeta,iphi,ieta);
558 //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
559 //printf("Max end---\n");
563 //________________________________________________________________
564 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
565 //Init EMCAL recalibration factors
566 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
567 //In order to avoid rewriting the same histograms
568 Bool_t oldStatus = TH1::AddDirectoryStatus();
569 TH1::AddDirectory(kFALSE);
571 fEMCALRecalibrationFactors = new TObjArray(10);
572 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));
573 //Init the histograms with 1
574 for (Int_t sm = 0; sm < 10; sm++) {
575 for (Int_t i = 0; i < 48; i++) {
576 for (Int_t j = 0; j < 24; j++) {
577 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
581 fEMCALRecalibrationFactors->SetOwner(kTRUE);
582 fEMCALRecalibrationFactors->Compress();
584 //In order to avoid rewriting the same histograms
585 TH1::AddDirectory(oldStatus);
589 //________________________________________________________________
590 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
591 //Init EMCAL bad channels map
592 AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
593 //In order to avoid rewriting the same histograms
594 Bool_t oldStatus = TH1::AddDirectoryStatus();
595 TH1::AddDirectory(kFALSE);
597 fEMCALBadChannelMap = new TObjArray(10);
598 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
599 for (int i = 0; i < 10; i++) {
600 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
605 fEMCALBadChannelMap->SetOwner(kTRUE);
606 fEMCALBadChannelMap->Compress();
608 //In order to avoid rewriting the same histograms
609 TH1::AddDirectory(oldStatus);
612 //________________________________________________________________
613 void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
614 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
616 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
617 UShort_t * index = cluster->GetCellsAbsId() ;
618 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
619 Int_t ncells = cluster->GetNCells();
621 //Initialize some used variables
624 Int_t icol = -1, irow = -1, imod=1;
625 Float_t factor = 1, frac = 0;
627 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
628 for(Int_t icell = 0; icell < ncells; icell++){
629 absId = index[icell];
630 frac = fraction[icell];
631 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
632 Int_t iTower = -1, iIphi = -1, iIeta = -1;
633 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
634 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
635 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
636 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
637 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
638 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
640 energy += cells->GetCellAmplitude(absId)*factor*frac;
644 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
646 cluster->SetE(energy);
651 //__________________________________________________
652 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
654 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
656 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
657 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
658 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
662 //__________________________________________________
663 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
665 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
666 // The algorithm is a copy of what is done in AliEMCALRecPoint
669 Float_t fraction = 1.;
670 Float_t recalFactor = 1.;
673 Int_t iTower = -1, iIphi = -1, iIeta = -1;
674 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
675 Float_t weight = 0., totalWeight=0.;
676 Float_t newPos[3] = {0,0,0};
677 Double_t pLocal[3], pGlobal[3];
678 Bool_t shared = kFALSE;
680 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
681 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
682 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
684 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
686 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
687 absId = clu->GetCellAbsId(iDig);
688 fraction = clu->GetCellAmplitudeFraction(iDig);
689 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
690 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
691 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
693 if(IsRecalibrationOn()) {
694 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
696 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
698 weight = GetCellWeight(eCell,clEnergy);
699 //printf("cell energy %f, weight %f\n",eCell,weight);
700 totalWeight += weight;
701 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
702 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
703 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
704 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
706 for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
711 for(int i=0; i<3; i++ ) newPos[i] /= totalWeight;
714 //Float_t pos[]={0,0,0};
715 //clu->GetPosition(pos);
716 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
717 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
719 if(iSupModMax > 1) {//sector 1
720 newPos[0] +=fMisalTransShift[3];//-=3.093;
721 newPos[1] +=fMisalTransShift[4];//+=6.82;
722 newPos[2] +=fMisalTransShift[5];//+=1.635;
723 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
727 newPos[0] +=fMisalTransShift[0];//+=1.134;
728 newPos[1] +=fMisalTransShift[1];//+=8.2;
729 newPos[2] +=fMisalTransShift[2];//+=1.197;
730 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
733 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
735 clu->SetPosition(newPos);
739 //__________________________________________________
740 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
742 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
743 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
746 Float_t fraction = 1.;
747 Float_t recalFactor = 1.;
751 Int_t iIphi = -1, iIeta = -1;
752 Int_t iSupMod = -1, iSupModMax = -1;
753 Int_t iphi = -1, ieta =-1;
754 Bool_t shared = kFALSE;
756 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
757 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi,shared);
758 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
760 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
761 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
762 Int_t startingSM = -1;
764 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
765 absId = clu->GetCellAbsId(iDig);
766 fraction = clu->GetCellAmplitudeFraction(iDig);
767 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
768 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
769 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
771 if (iDig==0) startingSM = iSupMod;
772 else if(iSupMod != startingSM) areInSameSM = kFALSE;
774 eCell = cells->GetCellAmplitude(absId);
776 if(IsRecalibrationOn()) {
777 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
779 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
781 weight = GetCellWeight(eCell,clEnergy);
782 if(weight < 0) weight = 0;
783 totalWeight += weight;
784 weightedCol += ieta*weight;
785 weightedRow += iphi*weight;
787 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
791 Float_t xyzNew[]={0.,0.,0.};
792 if(areInSameSM == kTRUE) {
793 //printf("In Same SM\n");
794 weightedCol = weightedCol/totalWeight;
795 weightedRow = weightedRow/totalWeight;
796 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
799 //printf("In Different SM\n");
800 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
803 clu->SetPosition(xyzNew);
807 //____________________________________________________________________________
808 void AliEMCALRecoUtils::RecalculateClusterDistanceToBadChannel(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster){
810 //re-evaluate distance to bad channel with updated bad map
812 if(!fRecalDistToBadChannels) return;
814 //Get channels map of the supermodule where the cluster is.
815 Int_t absIdMax = -1, iSupMod =-1, icolM = -1, irowM = -1;
816 Bool_t shared = kFALSE;
817 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSupMod, icolM, irowM, shared);
818 TH2D* hMap = (TH2D*)fEMCALBadChannelMap->At(iSupMod);
821 Float_t minDist = 10000.;
824 //Loop on tower status map
825 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
826 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
827 //Check if tower is bad.
828 if(hMap->GetBinContent(icol,irow)==0) continue;
829 //printf("AliEMCALRecoUtils::RecalculateDistanceToBadChannels() - \n \t Bad channel in SM %d, col %d, row %d, \n \t Cluster max in col %d, row %d\n",
830 // iSupMod,icol, irow, icolM,irowM);
832 dRrow=TMath::Abs(irowM-irow);
833 dRcol=TMath::Abs(icolM-icol);
834 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
836 //printf("MIN DISTANCE TO BAD %2.2f\n",dist);
843 //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module
848 //The only possible combinations are (0,1), (2,3) ... (8,9)
849 if(iSupMod%2) iSupMod2 = iSupMod-1;
850 else iSupMod2 = iSupMod+1;
851 hMap2 = (TH2D*)fEMCALBadChannelMap->At(iSupMod2);
853 //Loop on tower status map of second super module
854 for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){
855 for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){
856 //Check if tower is bad.
857 if(hMap2->GetBinContent(icol,irow)==0) continue;
858 //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",
859 // iSupMod2,icol, irow,iSupMod,icolM,irowM);
861 dRrow=TMath::Abs(irow-irowM);
864 dRcol=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+icolM));
867 dRcol=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-icolM);
870 dist=TMath::Sqrt(dRrow*dRrow+dRcol*dRcol);
871 if(dist < minDist) minDist = dist;
876 }// shared cluster in 2 SuperModules
878 AliDebug(2,Form("Max cluster cell (SM,col,row)=(%d %d %d) - Distance to Bad Channel %2.2f",iSupMod, icolM, irowM, minDist));
879 cluster->SetDistanceToBadChannel(minDist);
883 //____________________________________________________________________________
884 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
886 //re-evaluate identification parameters with bayesian
888 if ( cluster->GetM02() != 0)
889 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
891 Float_t pidlist[AliPID::kSPECIESN+1];
892 for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
894 cluster->SetPID(pidlist);
898 //____________________________________________________________________________
899 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
901 // Calculates new center of gravity in the local EMCAL-module coordinates
902 // and tranfers into global ALICE coordinates
903 // Calculates Dispersion and main axis
908 Float_t fraction = 1.;
909 Float_t recalFactor = 1.;
929 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
931 //Get from the absid the supermodule, tower and eta/phi numbers
932 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
933 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
935 //Get the cell energy, if recalibration is on, apply factors
936 fraction = cluster->GetCellAmplitudeFraction(iDigit);
937 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
938 if(IsRecalibrationOn()) {
939 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
941 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
943 if(cluster->E() > 0 && eCell > 0){
945 w = GetCellWeight(eCell,cluster->E());
953 dxx += w * etai * etai ;
955 dzz += w * phii * phii ;
957 dxz += w * etai * phii ;
961 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
964 //Normalize to the weight
970 AliError(Form("Wrong weight %f\n", wtot));
972 //Calculate dispersion
973 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
975 //Get from the absid the supermodule, tower and eta/phi numbers
976 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
977 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
979 //Get the cell energy, if recalibration is on, apply factors
980 fraction = cluster->GetCellAmplitudeFraction(iDigit);
981 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
982 if(IsRecalibrationOn()) {
983 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
985 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
987 if(cluster->E() > 0 && eCell > 0){
989 w = GetCellWeight(eCell,cluster->E());
993 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
996 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
999 //Normalize to the weigth and set shower shape parameters
1000 if (wtot > 0 && nstat > 1) {
1005 dxx -= xmean * xmean ;
1006 dzz -= zmean * zmean ;
1007 dxz -= xmean * zmean ;
1008 cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1009 cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
1013 cluster->SetM20(0.) ;
1014 cluster->SetM02(0.) ;
1018 cluster->SetDispersion(TMath::Sqrt(d)) ;
1020 cluster->SetDispersion(0) ;
1023 //____________________________________________________________________________
1024 void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr, AliEMCALGeometry *geom)
1026 //This function should be called before the cluster loop
1027 //Before call this function, please recalculate the cluster positions
1028 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
1029 //Store matched cluster indexes and residuals
1031 fMatchedTrackIndex ->Reset();
1032 fMatchedClusterIndex->Reset();
1033 fResidualZ ->Reset();
1034 fResidualR ->Reset();
1036 fMatchedTrackIndex ->Set(500);
1037 fMatchedClusterIndex->Set(500);
1038 fResidualZ ->Set(500);
1039 fResidualR ->Set(500);
1043 for (Int_t i=0; i<21;i++) cv[i]=0;
1044 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
1046 AliExternalTrackParam *trackParam=0;
1048 //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner
1049 if(dynamic_cast<AliESDEvent*> (event))
1051 AliESDtrack *esdTrack = dynamic_cast<AliESDEvent*> (event)->GetTrack(itr);
1052 if(!esdTrack || !IsAccepted(esdTrack)) continue;
1053 const AliESDfriendTrack* friendTrack = esdTrack->GetFriendTrack();
1054 if(friendTrack && friendTrack->GetTPCOut())
1056 //Use TPC Out as starting point if it is available
1057 trackParam= new AliExternalTrackParam(*friendTrack->GetTPCOut());
1061 //Otherwise use TPC inner
1062 trackParam = new AliExternalTrackParam(*esdTrack->GetInnerParam());
1066 //If the input event is AOD, the starting point for extrapolation is at vertex
1067 //AOD tracks are selected according to its bit.
1068 else if(dynamic_cast<AliAODEvent*> (event))
1070 AliAODTrack *aodTrack = dynamic_cast<AliAODEvent*> (event)->GetTrack(itr);
1071 if(!aodTrack) continue;
1072 if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
1073 Double_t pos[3],mom[3];
1074 aodTrack->GetXYZ(pos);
1075 aodTrack->GetPxPyPz(mom);
1076 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()));
1077 trackParam= new AliExternalTrackParam(pos,mom,cv,aodTrack->Charge());
1080 //Return if the input data is not "AOD" or "ESD"
1083 printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
1087 if(!trackParam) continue;
1089 Float_t dRMax = fCutR, dZMax=fCutZ;
1091 if(!clusterArr){// get clusters from event
1092 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1094 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
1095 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1096 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1097 Float_t tmpR=-1, tmpZ=-1;
1098 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
1107 } else { // external cluster array, not from ESD event
1108 for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
1110 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
1111 AliVCluster *cluster = (AliVCluster*) clusterArr->At(icl);
1112 if(!cluster->IsEMCAL()) continue;
1113 Float_t tmpR=-1, tmpZ=-1;
1114 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
1123 }// external list of clusters
1127 fMatchedTrackIndex ->AddAt(itr,matched);
1128 fMatchedClusterIndex->AddAt(index,matched);
1129 fResidualZ ->AddAt(dZMax,matched);
1130 fResidualR ->AddAt(dRMax,matched);
1136 AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
1138 fMatchedTrackIndex ->Set(matched);
1139 fMatchedClusterIndex->Set(matched);
1140 fResidualZ ->Set(matched);
1141 fResidualR ->Set(matched);
1144 //________________________________________________________________________________
1145 Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom)
1148 // This function returns the index of matched cluster to input track
1149 // Cut on match is dR<10cm by default. Returns -1 if no match is found
1152 Float_t dRMax = fCutR;
1155 AliExternalTrackParam *trackParam=0;
1156 const AliESDfriendTrack* friendTrack = track->GetFriendTrack();
1157 if(friendTrack && friendTrack->GetTPCOut())
1158 trackParam= const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
1160 trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
1162 if(!trackParam) return index;
1163 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
1165 AliExternalTrackParam *trkPamTmp = new AliExternalTrackParam(*trackParam);//Retrieve the starting point every time before the extrapolation
1166 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
1167 if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
1168 Float_t tmpR=-1, tmpZ=-1;
1169 if(!ExtrapolateTrackToCluster(trkPamTmp, cluster, tmpR, tmpZ)) continue;
1170 if(tmpR>-1 && tmpR<dRMax)
1180 //________________________________________________________________________________
1181 Bool_t AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpR, Float_t &tmpZ)
1184 //Return the residual by extrapolating a track to a cluster
1186 if(!trkParam || !cluster) return kFALSE;
1189 cluster->GetPosition(clsPos); //Has been recalculated
1190 TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
1191 Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
1192 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
1193 trkParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
1194 if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kFALSE)) return kFALSE;
1195 trkParam->GetXYZ(trkPos); //Get the extrapolated global position
1196 tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
1197 tmpZ = clsPos[2]-trkPos[2];
1201 //________________________________________________________________________________
1202 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t clsIndex, Float_t &dR, Float_t &dZ)
1204 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1205 //Get the residuals dR and dZ for this cluster to the closest track
1206 //Works with ESDs and AODs
1208 if( FindMatchedPosForCluster(clsIndex) >= 999 )
1210 AliDebug(2,"No matched tracks found!\n");
1215 dR = fResidualR->At(FindMatchedPosForCluster(clsIndex));
1216 dZ = fResidualZ->At(FindMatchedPosForCluster(clsIndex));
1218 //________________________________________________________________________________
1219 void AliEMCALRecoUtils::GetMatchedClusterResiduals(Int_t trkIndex, Float_t &dR, Float_t &dZ)
1221 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1222 //Get the residuals dR and dZ for this track to the closest cluster
1223 //Works with ESDs and AODs
1225 if( FindMatchedPosForTrack(trkIndex) >= 999 )
1227 AliDebug(2,"No matched cluster found!\n");
1232 dR = fResidualR->At(FindMatchedPosForTrack(trkIndex));
1233 dZ = fResidualZ->At(FindMatchedPosForTrack(trkIndex));
1236 //__________________________________________________________
1237 Int_t AliEMCALRecoUtils::GetMatchedTrackIndex(Int_t clsIndex)
1239 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1240 //Get the index of matched track to this cluster
1241 //Works with ESDs and AODs
1243 if(IsClusterMatched(clsIndex))
1244 return fMatchedTrackIndex->At(FindMatchedPosForCluster(clsIndex));
1249 //__________________________________________________________
1250 Int_t AliEMCALRecoUtils::GetMatchedClusterIndex(Int_t trkIndex)
1252 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1253 //Get the index of matched cluster to this track
1254 //Works with ESDs and AODs
1256 if(IsTrackMatched(trkIndex))
1257 return fMatchedClusterIndex->At(FindMatchedPosForTrack(trkIndex));
1262 //__________________________________________________
1263 Bool_t AliEMCALRecoUtils::IsClusterMatched(Int_t clsIndex)
1265 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1266 //Returns if the cluster has a match
1267 if(FindMatchedPosForCluster(clsIndex) < 999)
1273 //__________________________________________________
1274 Bool_t AliEMCALRecoUtils::IsTrackMatched(Int_t trkIndex)
1276 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1277 //Returns if the track has a match
1278 if(FindMatchedPosForTrack(trkIndex) < 999)
1284 //__________________________________________________________
1285 UInt_t AliEMCALRecoUtils::FindMatchedPosForCluster(Int_t clsIndex) const
1287 //Given a cluster index as in AliESDEvent::GetCaloCluster(clsIndex)
1288 //Returns the position of the match in the fMatchedClusterIndex array
1289 Float_t tmpR = fCutR;
1292 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
1294 if(fMatchedClusterIndex->At(i)==clsIndex && fResidualR->At(i)<tmpR)
1297 tmpR=fResidualR->At(i);
1298 AliDebug(3,Form("Matched cluster index: index: %d, dR: %2.4f, dZ: %2.4f.\n",fMatchedClusterIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
1304 //__________________________________________________________
1305 UInt_t AliEMCALRecoUtils::FindMatchedPosForTrack(Int_t trkIndex) const
1307 //Given a track index as in AliESDEvent::GetTrack(trkIndex)
1308 //Returns the position of the match in the fMatchedTrackIndex array
1309 Float_t tmpR = fCutR;
1312 for(Int_t i=0; i<fMatchedTrackIndex->GetSize(); i++)
1314 if(fMatchedTrackIndex->At(i)==trkIndex && fResidualR->At(i)<tmpR)
1317 tmpR=fResidualR->At(i);
1318 AliDebug(3,Form("Matched track index: index: %d, dR: %2.4f, dZ: %2.4f.\n",fMatchedTrackIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
1324 //__________________________________________________________
1325 Bool_t AliEMCALRecoUtils::IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells)
1327 // check if the cluster survives some quality cut
1330 Bool_t isGood=kTRUE;
1331 if(!cluster || !cluster->IsEMCAL()) return kFALSE;
1332 if(ClusterContainsBadChannel(geom,cluster->GetCellsAbsId(),cluster->GetNCells())) isGood=kFALSE;
1333 if(!CheckCellFiducialRegion(geom,cluster,cells)) isGood=kFALSE;
1334 if(fRejectExoticCluster && IsExoticCluster(cluster)) isGood=kFALSE;
1339 //__________________________________________________________
1340 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
1342 // Given a esd track, return whether the track survive all the cuts
1344 // The different quality parameter are first
1345 // retrieved from the track. then it is found out what cuts the
1346 // track did not survive and finally the cuts are imposed.
1348 UInt_t status = esdTrack->GetStatus();
1350 Int_t nClustersITS = esdTrack->GetITSclusters(0);
1351 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
1353 Float_t chi2PerClusterITS = -1;
1354 Float_t chi2PerClusterTPC = -1;
1355 if (nClustersITS!=0)
1356 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
1357 if (nClustersTPC!=0)
1358 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
1362 Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1363 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
1364 SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
1369 esdTrack->GetImpactParameters(b,bCov);
1370 if (bCov[0]<=0 || bCov[2]<=0) {
1371 AliDebug(1, "Estimated b resolution lower or equal zero!");
1372 bCov[0]=0; bCov[2]=0;
1375 Float_t dcaToVertexXY = b[0];
1376 Float_t dcaToVertexZ = b[1];
1377 Float_t dcaToVertex = -1;
1379 if (fCutDCAToVertex2D)
1380 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1382 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
1386 Bool_t cuts[kNCuts];
1387 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1389 // track quality cuts
1390 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1392 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1394 if (nClustersTPC<fCutMinNClusterTPC)
1396 if (nClustersITS<fCutMinNClusterITS)
1398 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1400 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1402 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1404 if (fCutDCAToVertex2D && dcaToVertex > 1)
1406 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1408 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1411 //Require at least one SPD point + anything else in ITS
1412 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1416 for (Int_t i=0; i<kNCuts; i++)
1417 if (cuts[i]) {cut = kTRUE;}
1425 //__________________________________________________
1426 void AliEMCALRecoUtils::InitTrackCuts()
1428 //Intilize the track cut criteria
1429 //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1430 //Also you can customize the cuts using the setters
1433 SetMinNClustersTPC(70);
1434 SetMaxChi2PerClusterTPC(4);
1435 SetAcceptKinkDaughters(kFALSE);
1436 SetRequireTPCRefit(kTRUE);
1439 SetRequireITSRefit(kTRUE);
1440 SetMaxDCAToVertexZ(2);
1441 SetDCAToVertex2D(kFALSE);
1442 SetMaxChi2PerClusterITS(); //which by default sets the value to 1e10.
1443 SetMinNClustersITS();
1446 //___________________________________________________
1447 void AliEMCALRecoUtils::Print(const Option_t *) const
1451 printf("AliEMCALRecoUtils Settings: \n");
1452 printf("Misalignment shifts\n");
1453 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,
1454 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1455 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
1456 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1457 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
1459 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
1461 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1462 printf("Mass hypothesis = %2.3f[GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
1464 printf("Track cuts: \n");
1465 printf("AOD track selection mask: %d\n",fAODFilterMask);
1466 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1467 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1468 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1469 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1470 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);
1474 //_____________________________________________________________________
1475 void AliEMCALRecoUtils::SetTimeDependentCorrections(Int_t runnumber){
1476 //Get EMCAL time dependent corrections from file and put them in the recalibration histograms
1477 //Do it only once and only if it is requested
1479 if(!fUseTimeCorrectionFactors) return;
1480 if(fTimeCorrectionFactorsSet) return;
1482 printf("AliEMCALRecoUtils::GetTimeDependentCorrections() - Get Correction Factors for Run number %d\n",runnumber);
1484 AliEMCALCalibTimeDepCorrection *corr = new AliEMCALCalibTimeDepCorrection();
1485 corr->ReadRootInfo(Form("CorrectionFiles/Run%d_Correction.root",runnumber));
1487 SwitchOnRecalibration();
1488 for(Int_t ism = 0; ism < 4; ism++){
1489 for(Int_t icol = 0; icol < 48; icol++){
1490 for(Int_t irow = 0; irow < 24; irow++){
1491 Float_t orgRecalFactor = GetEMCALChannelRecalibrationFactors(ism)->GetBinContent(icol,irow);
1492 Float_t newRecalFactor = orgRecalFactor*corr->GetCorrection(ism, icol,irow,0);
1493 GetEMCALChannelRecalibrationFactors(ism)->SetBinContent(icol,irow,newRecalFactor);
1494 //printf("ism %d, icol %d, irow %d, corrections : org %f, time dep %f, final %f (org*time %f)\n",ism, icol, irow,
1495 // orgRecalFactor, corr->GetCorrection(ism, icol,irow,0),
1496 // (GetEMCALChannelRecalibrationFactors(ism))->GetBinContent(icol,irow),newRecalFactor);
1500 fTimeCorrectionFactorsSet = kTRUE;