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 ///////////////////////////////////////////////////////////////////////////////
29 // standard C++ includes
30 //#include <Riostream.h>
33 #include <TGeoManager.h>
34 #include <TGeoMatrix.h>
38 #include "AliEMCALRecoUtils.h"
39 #include "AliEMCALGeometry.h"
40 #include "AliVCluster.h"
41 #include "AliVCaloCells.h"
42 #include "AliVEvent.h"
44 #include "AliEMCALPIDUtils.h"
46 #include "AliESDEvent.h"
47 #include "AliESDtrack.h"
48 #include "AliEMCALTrack.h"
50 ClassImp(AliEMCALRecoUtils)
52 //______________________________________________
53 AliEMCALRecoUtils::AliEMCALRecoUtils():
54 fNonLinearityFunction (kNoCorrection), fParticleType(kPhoton),
55 fPosAlgo(kUnchanged),fW0(4.),
56 fRecalibration(kFALSE), fEMCALRecalibrationFactors(),
57 fRemoveBadChannels(kFALSE), fEMCALBadChannelMap(),
58 fNCellsFromEMCALBorder(0), fNoEMCALBorderAtEta0(kTRUE),
59 fMatchedClusterIndex(0x0), fResidualZ(0x0), fResidualR(0x0), fCutR(20), fCutZ(20),
60 fCutMinNClusterTPC(0), fCutMinNClusterITS(0), fCutMaxChi2PerClusterTPC(0), fCutMaxChi2PerClusterITS(0),
61 fCutRequireTPCRefit(0), fCutRequireITSRefit(0), fCutAcceptKinkDaughters(0),
62 fCutMaxDCAToVertexXY(0), fCutMaxDCAToVertexZ(0),fCutDCAToVertex2D(0),
67 // Initialize all constant values which have to be used
68 // during Reco algorithm execution
71 for(Int_t i = 0; i < 15 ; i++) {
72 fMisalTransShift[i] = 0.;
73 fMisalRotShift[i] = 0.;
75 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = 0.;
76 //For kPi0GammaGamma case, but default is no correction
77 fNonLinearityParams[0] = 0.1457/0.1349766/1.038;
78 fNonLinearityParams[1] = -0.02024/0.1349766/1.038;
79 fNonLinearityParams[2] = 1.046;
81 fMatchedClusterIndex = new TArrayI();
82 fResidualZ = new TArrayF();
83 fResidualR = new TArrayF();
85 fPIDUtils = new AliEMCALPIDUtils();
91 //______________________________________________________________________
92 AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
93 : TNamed(reco), fNonLinearityFunction(reco.fNonLinearityFunction),
94 fParticleType(reco.fParticleType), fPosAlgo(reco.fPosAlgo), fW0(reco.fW0),
95 fRecalibration(reco.fRecalibration),fEMCALRecalibrationFactors(reco.fEMCALRecalibrationFactors),
96 fRemoveBadChannels(reco.fRemoveBadChannels),fEMCALBadChannelMap(reco.fEMCALBadChannelMap),
97 fNCellsFromEMCALBorder(reco.fNCellsFromEMCALBorder),fNoEMCALBorderAtEta0(reco.fNoEMCALBorderAtEta0),
98 fMatchedClusterIndex(reco.fMatchedClusterIndex?new TArrayI(*reco.fMatchedClusterIndex):0x0),
99 fResidualZ(reco.fResidualZ?new TArrayF(*reco.fResidualZ):0x0),
100 fResidualR(reco.fResidualR?new TArrayF(*reco.fResidualR):0x0),
101 fCutR(reco.fCutR),fCutZ(reco.fCutZ),
102 fCutMinNClusterTPC(reco.fCutMinNClusterTPC), fCutMinNClusterITS(reco.fCutMinNClusterITS),
103 fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC), fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
104 fCutRequireTPCRefit(reco.fCutRequireTPCRefit), fCutRequireITSRefit(reco.fCutRequireITSRefit),
105 fCutAcceptKinkDaughters(reco.fCutAcceptKinkDaughters),
106 fCutMaxDCAToVertexXY(reco.fCutMaxDCAToVertexXY), fCutMaxDCAToVertexZ(reco.fCutMaxDCAToVertexZ),fCutDCAToVertex2D(reco.fCutDCAToVertex2D),
107 fPIDUtils(reco.fPIDUtils)
112 for(Int_t i = 0; i < 15 ; i++) {
113 fMisalRotShift[i] = reco.fMisalRotShift[i];
114 fMisalTransShift[i] = reco.fMisalTransShift[i];
116 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
121 //______________________________________________________________________
122 AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & reco)
124 //Assignment operator
126 if(this == &reco)return *this;
127 ((TNamed *)this)->operator=(reco);
129 fNonLinearityFunction = reco.fNonLinearityFunction;
130 fParticleType = reco.fParticleType;
131 fPosAlgo = reco.fPosAlgo;
133 fRecalibration = reco.fRecalibration;
134 fEMCALRecalibrationFactors = reco.fEMCALRecalibrationFactors;
135 fRemoveBadChannels = reco.fRemoveBadChannels;
136 fEMCALBadChannelMap = reco.fEMCALBadChannelMap;
137 fNCellsFromEMCALBorder = reco.fNCellsFromEMCALBorder;
138 fNoEMCALBorderAtEta0 = reco.fNoEMCALBorderAtEta0;
141 for(Int_t i = 0; i < 15 ; i++) {fMisalTransShift[i] = reco.fMisalTransShift[i]; fMisalRotShift[i] = reco.fMisalRotShift[i];}
142 for(Int_t i = 0; i < 6 ; i++) fNonLinearityParams[i] = reco.fNonLinearityParams[i];
147 fCutMinNClusterTPC = reco.fCutMinNClusterTPC;
148 fCutMinNClusterITS = reco.fCutMinNClusterITS;
149 fCutMaxChi2PerClusterTPC = reco.fCutMaxChi2PerClusterTPC;
150 fCutMaxChi2PerClusterITS = reco.fCutMaxChi2PerClusterITS;
151 fCutRequireTPCRefit = reco.fCutRequireTPCRefit;
152 fCutRequireITSRefit = reco.fCutRequireITSRefit;
153 fCutAcceptKinkDaughters = reco.fCutAcceptKinkDaughters;
154 fCutMaxDCAToVertexXY = reco.fCutMaxDCAToVertexXY;
155 fCutMaxDCAToVertexZ = reco.fCutMaxDCAToVertexZ;
156 fCutDCAToVertex2D = reco.fCutDCAToVertex2D;
158 fPIDUtils = reco.fPIDUtils;
162 // assign or copy construct
164 *fResidualR = *reco.fResidualR;
166 else fResidualR = new TArrayF(*reco.fResidualR);
169 if(fResidualR)delete fResidualR;
174 // assign or copy construct
176 *fResidualZ = *reco.fResidualZ;
178 else fResidualZ = new TArrayF(*reco.fResidualZ);
181 if(fResidualZ)delete fResidualZ;
186 if(reco.fMatchedClusterIndex){
187 // assign or copy construct
188 if(fMatchedClusterIndex){
189 *fMatchedClusterIndex = *reco.fMatchedClusterIndex;
191 else fMatchedClusterIndex = new TArrayI(*reco.fMatchedClusterIndex);
194 if(fMatchedClusterIndex)delete fMatchedClusterIndex;
195 fMatchedClusterIndex = 0;
203 //__________________________________________________
204 AliEMCALRecoUtils::~AliEMCALRecoUtils()
208 if(fEMCALRecalibrationFactors) {
209 fEMCALRecalibrationFactors->Clear();
210 delete fEMCALRecalibrationFactors;
213 if(fEMCALBadChannelMap) {
214 fEMCALBadChannelMap->Clear();
215 delete fEMCALBadChannelMap;
218 if(fMatchedClusterIndex) {delete fMatchedClusterIndex; fMatchedClusterIndex=0;}
219 if(fResidualR) {delete fResidualR; fResidualR=0;}
220 if(fResidualZ) {delete fResidualZ; fResidualZ=0;}
224 //_______________________________________________________________
225 Bool_t AliEMCALRecoUtils::CheckCellFiducialRegion(AliEMCALGeometry* geom, AliVCluster* cluster, AliVCaloCells* cells)
227 // Given the list of AbsId of the cluster, get the maximum cell and
228 // check if there are fNCellsFromBorder from the calorimeter border
230 //If the distance to the border is 0 or negative just exit accept all clusters
231 if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
233 Int_t absIdMax = -1, iSM =-1, ieta = -1, iphi = -1;
234 GetMaxEnergyCell(geom, cells, cluster, absIdMax, iSM, ieta, iphi);
236 AliDebug(2,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f, Ncells from border %d, EMCAL eta=0 %d\n",
237 absIdMax, cells->GetCellAmplitude(absIdMax), cluster->E(), fNCellsFromEMCALBorder, fNoEMCALBorderAtEta0));
239 if(absIdMax==-1) return kFALSE;
241 //Check if the cell is close to the borders:
242 Bool_t okrow = kFALSE;
243 Bool_t okcol = kFALSE;
245 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
246 AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
252 if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
255 if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
259 if(!fNoEMCALBorderAtEta0){
260 if(ieta > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
264 if(ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
267 if(ieta < 48-fNCellsFromEMCALBorder) okcol = kTRUE;
271 AliDebug(2,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d: column? %d, row? %d\nq",
272 fNCellsFromEMCALBorder, ieta, iphi, iSM, okcol, okrow));
274 if (okcol && okrow) {
275 //printf("Accept\n");
279 //printf("Reject\n");
280 AliDebug(2,Form("Reject cluster in border, max cell : ieta %d, iphi %d, SM %d\n",ieta, iphi, iSM));
287 //_________________________________________________________________________________________________________
288 Bool_t AliEMCALRecoUtils::ClusterContainsBadChannel(AliEMCALGeometry* geom, UShort_t* cellList, Int_t nCells){
289 // Check that in the cluster cells, there is no bad channel of those stored
290 // in fEMCALBadChannelMap or fPHOSBadChannelMap
292 if(!fRemoveBadChannels) return kFALSE;
293 if(!fEMCALBadChannelMap) return kFALSE;
298 for(Int_t iCell = 0; iCell<nCells; iCell++){
300 //Get the column and row
301 Int_t iTower = -1, iIphi = -1, iIeta = -1;
302 geom->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
303 if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
304 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
305 if(GetEMCALChannelStatus(imod, icol, irow)){
306 AliDebug(2,Form("Cluster with bad channel: SM %d, col %d, row %d\n",imod, icol, irow));
310 }// cell cluster loop
316 //__________________________________________________
317 Float_t AliEMCALRecoUtils::CorrectClusterEnergyLinearity(AliVCluster* cluster){
318 // Correct cluster energy from non linearity functions
319 Float_t energy = cluster->E();
321 switch (fNonLinearityFunction) {
324 //Non-Linearity correction (from MC with function ([0]*exp(-[1]/E))+(([2]/([3]*2.*TMath::Pi())*exp(-(E-[4])^2/(2.*[3]^2)))))
325 //Double_t par0 = 1.001;
326 //Double_t par1 = -0.01264;
327 //Double_t par2 = -0.03632;
328 //Double_t par3 = 0.1798;
329 //Double_t par4 = -0.522;
330 energy /= (fNonLinearityParams[0]*exp(-fNonLinearityParams[1]/energy))+
331 ((fNonLinearityParams[2]/(fNonLinearityParams[3]*2.*TMath::Pi())*
332 exp(-(energy-fNonLinearityParams[4])*(energy-fNonLinearityParams[4])/(2.*fNonLinearityParams[3]*fNonLinearityParams[3]))));
337 //Non-Linearity correction (from Olga Data with function p0+p1*exp(-p2*E))
338 //Double_t par0 = 0.1457;
339 //Double_t par1 = -0.02024;
340 //Double_t par2 = 1.046;
341 energy /= (fNonLinearityParams[0]+fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy)); //Olga function
344 case kPi0GammaConversion:
346 //Non-Linearity correction (Nicolas from Dimitri Data with function C*[1-a*exp(-b*E)])
347 //Double_t C = 0.139393/0.1349766;
348 //Double_t a = 0.0566186;
349 //Double_t b = 0.982133;
350 energy /= fNonLinearityParams[0]*(1-fNonLinearityParams[1]*exp(-fNonLinearityParams[2]*energy));
355 AliDebug(2,"No correction on the energy\n");
364 //__________________________________________________
365 Float_t AliEMCALRecoUtils::GetDepth(const Float_t energy, const Int_t iParticle, const Int_t iSM) const
367 //Calculate shower depth for a given cluster energy and particle type
377 depth = x0 * (TMath::Log(energy*1000/ ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
381 depth = x0 * (TMath::Log(energy*1000/ ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
388 gGeoManager->cd("ALIC_1/XEN1_1");
389 TGeoNode *geoXEn1 = gGeoManager->GetCurrentNode();
390 TGeoNodeMatrix *geoSM = dynamic_cast<TGeoNodeMatrix *>(geoXEn1->GetDaughter(iSM));
392 TGeoVolume *geoSMVol = geoSM->GetVolume();
393 TGeoShape *geoSMShape = geoSMVol->GetShape();
394 TGeoBBox *geoBox = dynamic_cast<TGeoBBox *>(geoSMShape);
395 if(geoBox) depth = 0.5 * geoBox->GetDX()*2 ;
396 else AliFatal("Null GEANT box");
397 }else AliFatal("NULL GEANT node matrix");
400 depth = x0 * (TMath::Log(energy*1000 / ecr) - 0.5); //Multiply energy by 1000 to transform to MeV
406 depth = x0 * (TMath::Log(energy*1000 / ecr) + 0.5); //Multiply energy by 1000 to transform to MeV
413 //__________________________________________________
414 void AliEMCALRecoUtils::GetMaxEnergyCell(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu, Int_t & absId, Int_t& iSupMod, Int_t& ieta, Int_t& iphi)
416 //For a given CaloCluster gets the absId of the cell
417 //with maximum energy deposit.
420 Double_t eCell = -1.;
421 Float_t fraction = 1.;
422 Float_t recalFactor = 1.;
423 Int_t cellAbsId = -1 ;
428 //printf("---Max?\n");
429 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
430 cellAbsId = clu->GetCellAbsId(iDig);
431 fraction = clu->GetCellAmplitudeFraction(iDig);
432 //printf("a Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,cells->GetCellAmplitude(cellAbsId),fraction);
433 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
434 if(IsRecalibrationOn()) {
435 geom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
436 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
437 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
439 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
440 //printf("b Cell %d, id, %d, amp %f, fraction %f\n",iDig,cellAbsId,eCell,fraction);
444 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
448 //Get from the absid the supermodule, tower and eta/phi numbers
449 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
450 //Gives SuperModule and Tower numbers
451 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
452 iIphi, iIeta,iphi,ieta);
453 //printf("Max id %d, iSM %d, col %d, row %d\n",absId,iSupMod,ieta,iphi);
454 //printf("Max end---\n");
458 //________________________________________________________________
459 void AliEMCALRecoUtils::InitEMCALRecalibrationFactors(){
460 //Init EMCAL recalibration factors
461 AliDebug(2,"AliCalorimeterUtils::InitEMCALRecalibrationFactors()");
462 //In order to avoid rewriting the same histograms
463 Bool_t oldStatus = TH1::AddDirectoryStatus();
464 TH1::AddDirectory(kFALSE);
466 fEMCALRecalibrationFactors = new TObjArray(12);
467 for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
468 //Init the histograms with 1
469 for (Int_t sm = 0; sm < 12; sm++) {
470 for (Int_t i = 0; i < 48; i++) {
471 for (Int_t j = 0; j < 24; j++) {
472 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
476 fEMCALRecalibrationFactors->SetOwner(kTRUE);
477 fEMCALRecalibrationFactors->Compress();
479 //In order to avoid rewriting the same histograms
480 TH1::AddDirectory(oldStatus);
484 //________________________________________________________________
485 void AliEMCALRecoUtils::InitEMCALBadChannelStatusMap(){
486 //Init EMCAL bad channels map
487 AliDebug(2,"AliEMCALRecoUtils::InitEMCALBadChannelStatusMap()");
488 //In order to avoid rewriting the same histograms
489 Bool_t oldStatus = TH1::AddDirectoryStatus();
490 TH1::AddDirectory(kFALSE);
492 fEMCALBadChannelMap = new TObjArray(12);
493 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
494 for (int i = 0; i < 12; i++) {
495 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
500 fEMCALBadChannelMap->SetOwner(kTRUE);
501 fEMCALBadChannelMap->Compress();
503 //In order to avoid rewriting the same histograms
504 TH1::AddDirectory(oldStatus);
507 //________________________________________________________________
508 void AliEMCALRecoUtils::RecalibrateClusterEnergy(AliEMCALGeometry* geom, AliVCluster * cluster, AliVCaloCells * cells){
509 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
511 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
512 UShort_t * index = cluster->GetCellsAbsId() ;
513 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
514 Int_t ncells = cluster->GetNCells();
516 //Initialize some used variables
519 Int_t icol = -1, irow = -1, imod=1;
520 Float_t factor = 1, frac = 0;
522 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
523 for(Int_t icell = 0; icell < ncells; icell++){
524 absId = index[icell];
525 frac = fraction[icell];
526 if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
527 Int_t iTower = -1, iIphi = -1, iIeta = -1;
528 geom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
529 if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
530 geom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
531 factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
532 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
533 imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId)));
535 energy += cells->GetCellAmplitude(absId)*factor*frac;
539 AliDebug(2,Form("AliEMCALRecoUtils::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy));
541 cluster->SetE(energy);
546 //__________________________________________________
547 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
549 //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
551 if (fPosAlgo==kPosTowerGlobal) RecalculateClusterPositionFromTowerGlobal( geom, cells, clu);
552 else if(fPosAlgo==kPosTowerIndex) RecalculateClusterPositionFromTowerIndex ( geom, cells, clu);
553 else AliDebug(2,"Algorithm to recalculate position not selected, do nothing.");
557 //__________________________________________________
558 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerGlobal(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
560 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
561 // The algorithm is a copy of what is done in AliEMCALRecPoint
564 Float_t fraction = 1.;
565 Float_t recalFactor = 1.;
568 Int_t iTower = -1, iIphi = -1, iIeta = -1;
569 Int_t iSupModMax = -1, iSM=-1, iphi = -1, ieta = -1;
570 Float_t weight = 0., totalWeight=0.;
571 Float_t newPos[3] = {0,0,0};
572 Double_t pLocal[3], pGlobal[3];
574 Float_t clEnergy = clu->E(); //Energy already recalibrated previously
575 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi);
576 Double_t depth = GetDepth(clEnergy,fParticleType,iSupModMax) ;
578 //printf("** Cluster energy %f, ncells %d, depth %f\n",clEnergy,clu->GetNCells(),depth);
580 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
581 absId = clu->GetCellAbsId(iDig);
582 fraction = clu->GetCellAmplitudeFraction(iDig);
583 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
584 geom->GetCellIndex(absId,iSM,iTower,iIphi,iIeta);
585 geom->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
587 if(IsRecalibrationOn()) {
588 recalFactor = GetEMCALChannelRecalibrationFactor(iSM,ieta,iphi);
590 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
592 weight = GetCellWeight(eCell,clEnergy);
593 //printf("cell energy %f, weight %f\n",eCell,weight);
594 totalWeight += weight;
595 geom->RelPosCellInSModule(absId,depth,pLocal[0],pLocal[1],pLocal[2]);
596 //printf("pLocal (%f,%f,%f), SM %d, absId %d\n",pLocal[0],pLocal[1],pLocal[2],iSupModMax,absId);
597 geom->GetGlobal(pLocal,pGlobal,iSupModMax);
598 //printf("pLocal (%f,%f,%f)\n",pGlobal[0],pGlobal[1],pGlobal[2]);
600 for(int i=0; i<3; i++ ) newPos[i] += (weight*pGlobal[i]);
605 for(int i=0; i<3; i++ ) newPos[i] /= totalWeight;
608 //Float_t pos[]={0,0,0};
609 //clu->GetPosition(pos);
610 //printf("OldPos : %2.3f,%2.3f,%2.3f\n",pos[0],pos[1],pos[2]);
611 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
613 if(iSupModMax > 1) {//sector 1
614 newPos[0] +=fMisalTransShift[3];//-=3.093;
615 newPos[1] +=fMisalTransShift[4];//+=6.82;
616 newPos[2] +=fMisalTransShift[5];//+=1.635;
617 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[3],fMisalTransShift[4],fMisalTransShift[5]);
621 newPos[0] +=fMisalTransShift[0];//+=1.134;
622 newPos[1] +=fMisalTransShift[1];//+=8.2;
623 newPos[2] +=fMisalTransShift[2];//+=1.197;
624 //printf(" + : %2.3f,%2.3f,%2.3f\n",fMisalTransShift[0],fMisalTransShift[1],fMisalTransShift[2]);
627 //printf("NewPos : %2.3f,%2.3f,%2.3f\n",newPos[0],newPos[1],newPos[2]);
629 clu->SetPosition(newPos);
633 //__________________________________________________
634 void AliEMCALRecoUtils::RecalculateClusterPositionFromTowerIndex(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
636 // For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
637 // The algorithm works with the tower indeces, averages the indeces and from them it calculates the global position
640 Float_t fraction = 1.;
641 Float_t recalFactor = 1.;
645 Int_t iIphi = -1, iIeta = -1;
646 Int_t iSupMod = -1, iSupModMax = -1;
647 Int_t iphi = -1, ieta =-1;
649 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
650 GetMaxEnergyCell(geom, cells, clu, absId, iSupModMax, ieta, iphi);
651 Float_t depth = GetDepth(clEnergy,fParticleType,iSupMod) ;
653 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
654 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
655 Int_t startingSM = -1;
657 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
658 absId = clu->GetCellAbsId(iDig);
659 fraction = clu->GetCellAmplitudeFraction(iDig);
660 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
661 geom->GetCellIndex(absId,iSupMod,iTower,iIphi,iIeta);
662 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
664 if (iDig==0) startingSM = iSupMod;
665 else if(iSupMod != startingSM) areInSameSM = kFALSE;
667 eCell = cells->GetCellAmplitude(absId);
669 if(IsRecalibrationOn()) {
670 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
672 eCell = cells->GetCellAmplitude(absId)*fraction*recalFactor;
674 weight = GetCellWeight(eCell,clEnergy);
675 if(weight < 0) weight = 0;
676 totalWeight += weight;
677 weightedCol += ieta*weight;
678 weightedRow += iphi*weight;
680 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
684 Float_t xyzNew[]={0.,0.,0.};
685 if(areInSameSM == kTRUE) {
686 //printf("In Same SM\n");
687 weightedCol = weightedCol/totalWeight;
688 weightedRow = weightedRow/totalWeight;
689 geom->RecalculateTowerPosition(weightedRow, weightedCol, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
692 //printf("In Different SM\n");
693 geom->RecalculateTowerPosition(iphi, ieta, iSupModMax, depth, fMisalTransShift, fMisalRotShift, xyzNew);
696 clu->SetPosition(xyzNew);
700 //____________________________________________________________________________
701 void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster){
703 //re-evaluate identification parameters with bayesian
705 if ( cluster->GetM02() != 0)
706 fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
708 Float_t pidlist[AliPID::kSPECIESN+1];
709 for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
711 cluster->SetPID(pidlist);
715 //____________________________________________________________________________
716 void AliEMCALRecoUtils::RecalculateClusterShowerShapeParameters(AliEMCALGeometry * geom, AliVCaloCells* cells, AliVCluster * cluster)
718 // Calculates new center of gravity in the local EMCAL-module coordinates
719 // and tranfers into global ALICE coordinates
720 // Calculates Dispersion and main axis
725 Float_t fraction = 1.;
726 Float_t recalFactor = 1.;
746 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
748 //Get from the absid the supermodule, tower and eta/phi numbers
749 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
750 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
752 //Get the cell energy, if recalibration is on, apply factors
753 fraction = cluster->GetCellAmplitudeFraction(iDigit);
754 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
755 if(IsRecalibrationOn()) {
756 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
758 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
760 if(cluster->E() > 0 && eCell > 0){
762 w = GetCellWeight(eCell,cluster->E());
770 dxx += w * etai * etai ;
772 dzz += w * phii * phii ;
774 dxz += w * etai * phii ;
778 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
781 //Normalize to the weight
787 AliError(Form("Wrong weight %f\n", wtot));
789 //Calculate dispersion
790 for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++) {
792 //Get from the absid the supermodule, tower and eta/phi numbers
793 geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
794 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
796 //Get the cell energy, if recalibration is on, apply factors
797 fraction = cluster->GetCellAmplitudeFraction(iDigit);
798 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
799 if(IsRecalibrationOn()) {
800 recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
802 eCell = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
804 if(cluster->E() > 0 && eCell > 0){
806 w = GetCellWeight(eCell,cluster->E());
810 if(w > 0.0) d += w*((etai-xmean)*(etai-xmean)+(phii-zmean)*(phii-zmean));
813 AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, cluster->E()));
816 //Normalize to the weigth and set shower shape parameters
817 if (wtot > 0 && nstat > 1) {
822 dxx -= xmean * xmean ;
823 dzz -= zmean * zmean ;
824 dxz -= xmean * zmean ;
825 cluster->SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
826 cluster->SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
830 cluster->SetM20(0.) ;
831 cluster->SetM02(0.) ;
835 cluster->SetDispersion(TMath::Sqrt(d)) ;
837 cluster->SetDispersion(0) ;
841 //__________________________________________________
842 void AliEMCALRecoUtils::FindMatches(AliVEvent *event)
844 //This function should be called before the cluster loop
845 //Before call this function, please recalculate the cluster positions
846 //Given the input event, loop over all the tracks, select the closest cluster as matched with fCutR
847 //Store matched cluster indexes and residuals
848 //It only works with ESDs, not AODs
850 fMatchedClusterIndex->Reset();
854 fMatchedClusterIndex->Set(100);
855 fResidualZ->Set(100);
856 fResidualR->Set(100);
861 for(Int_t itr=0; itr<event->GetNumberOfTracks(); itr++)
863 AliESDtrack *track = ((AliESDEvent*)event)->GetTrack(itr);
864 if(!track || !IsAccepted(track)) continue;
866 Float_t dRMax = fCutR, dZMax = fCutZ;
868 AliEMCALTrack *emctrack = new AliEMCALTrack(*track);
869 for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
871 AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
872 if(!cluster->IsEMCAL()) continue;
873 cluster->GetPosition(clsPos); //Has been recalculated
874 if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) continue;
875 emctrack->GetXYZ(trkPos);
876 Float_t tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
877 Float_t tmpZ = TMath::Abs(clsPos[2]-trkPos[2]);
890 fMatchedClusterIndex->AddAt(index,matched);
891 fResidualZ->AddAt(dZMax,matched);
892 fResidualR->AddAt(dRMax,matched);
897 fMatchedClusterIndex->Set(matched);
898 fResidualZ->Set(matched);
899 fResidualR->Set(matched);
901 //printf("Number of matched pairs: %d\n",matched);
904 //__________________________________________________
905 void AliEMCALRecoUtils::GetMatchedResiduals(Int_t index, Float_t &dR, Float_t &dZ)
907 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
908 //Get the residuals dR and dZ for this cluster
909 //It only works with ESDs, not AODs
911 if( FindMatchedPos(index)==-1 )
913 AliDebug(2,"No matched tracks found!\n");
918 dR = fResidualR->At(FindMatchedPos(index));
919 dZ = fResidualZ->At(FindMatchedPos(index));
922 //__________________________________________________
923 Bool_t AliEMCALRecoUtils::IsMatched(Int_t index)
925 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
926 //Returns if cluster has a match
927 if(FindMatchedPos(index)>-1)
932 //__________________________________________________
933 Int_t AliEMCALRecoUtils::FindMatchedPos(Int_t index) const
935 //Given a cluster index as in AliESDEvent::GetCaloCluster(index)
936 //Returns the position of the match in the fMatchedClusterIndex array
937 Float_t tmpR = fCutR;
940 for(Int_t i=0; i<fMatchedClusterIndex->GetSize(); i++)
942 if(fMatchedClusterIndex->At(i)==index && fResidualR->At(i)<tmpR)
945 tmpR=fResidualR->At(i);
947 AliDebug(3,Form("Matched cluster pos: %d, index: %d, dR: %2.4f, dZ: %2.4f.\n",i,fMatchedClusterIndex->At(i),fResidualR->At(i),fResidualZ->At(i)));
952 Bool_t AliEMCALRecoUtils::IsAccepted(AliESDtrack *esdTrack)
954 // Given a esd track, return whether the track survive all the cuts
956 // The different quality parameter are first
957 // retrieved from the track. then it is found out what cuts the
958 // track did not survive and finally the cuts are imposed.
960 UInt_t status = esdTrack->GetStatus();
962 Int_t nClustersITS = esdTrack->GetITSclusters(0);
963 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
965 Float_t chi2PerClusterITS = -1;
966 Float_t chi2PerClusterTPC = -1;
968 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
970 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
974 Float_t MaxDCAToVertexXYPtDep = 0.0182 + 0.0350/TMath::Power(esdTrack->Pt(),1.01); //This expression comes from AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
975 //AliDebug(3,Form("Track pT = %f, DCAtoVertexXY = %f",esdTrack->Pt(),MaxDCAToVertexXYPtDep));
976 SetMaxDCAToVertexXY(MaxDCAToVertexXYPtDep); //Set pT dependent DCA cut to vertex in x-y plane
981 esdTrack->GetImpactParameters(b,bCov);
982 if (bCov[0]<=0 || bCov[2]<=0) {
983 AliDebug(1, "Estimated b resolution lower or equal zero!");
984 bCov[0]=0; bCov[2]=0;
987 Float_t dcaToVertexXY = b[0];
988 Float_t dcaToVertexZ = b[1];
989 Float_t dcaToVertex = -1;
991 if (fCutDCAToVertex2D)
992 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
994 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
999 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
1001 // track quality cuts
1002 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
1004 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
1006 if (nClustersTPC<fCutMinNClusterTPC)
1008 if (nClustersITS<fCutMinNClusterITS)
1010 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
1012 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
1014 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
1016 if (fCutDCAToVertex2D && dcaToVertex > 1)
1018 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
1020 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
1023 //Require at least one SPD point + anything else in ITS
1024 if( (esdTrack->HasPointOnITSLayer(0) || esdTrack->HasPointOnITSLayer(1)) == kFALSE)
1028 for (Int_t i=0; i<kNCuts; i++)
1029 if (cuts[i]) {cut = kTRUE;}
1037 //__________________________________________________
1038 void AliEMCALRecoUtils::InitTrackCuts()
1040 //Intilize the track cut criteria
1041 //By default these cuts are set according to AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()
1042 //Also you can customize the cuts using the setters
1045 SetMinNClustersTPC(70);
1046 SetMaxChi2PerClusterTPC(4);
1047 SetAcceptKinkDaughters(kFALSE);
1048 SetRequireTPCRefit(kTRUE);
1051 SetRequireITSRefit(kTRUE);
1052 SetMaxDCAToVertexZ(2);
1053 SetDCAToVertex2D(kFALSE);
1056 //__________________________________________________
1057 void AliEMCALRecoUtils::Print(const Option_t *) const
1061 printf("AliEMCALRecoUtils Settings: \n");
1062 printf("Misalignment shifts\n");
1063 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,
1064 fMisalTransShift[i*3],fMisalTransShift[i*3+1],fMisalTransShift[i*3+2],
1065 fMisalRotShift[i*3], fMisalRotShift[i*3+1], fMisalRotShift[i*3+2] );
1066 printf("Non linearity function %d, parameters:\n", fNonLinearityFunction);
1067 for(Int_t i=0; i<6; i++) printf("param[%d]=%f\n",i, fNonLinearityParams[i]);
1069 printf("Position Recalculation option %d, Particle Type %d, fW0 %2.2f, Recalibrate Data %d \n",fPosAlgo,fParticleType,fW0, fRecalibration);
1071 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1073 printf("Track cuts: \n");
1074 printf("TPCRefit = %d, ITSRefit = %d\n",fCutRequireTPCRefit,fCutRequireITSRefit);
1075 printf("AcceptKinks = %d\n",fCutAcceptKinkDaughters);
1076 printf("MinNCulsterTPC = %d, MinNClusterITS = %d\n",fCutMinNClusterTPC,fCutMinNClusterITS);
1077 printf("MaxChi2TPC = %2.2f, MaxChi2ITS = %2.2f\n",fCutMaxChi2PerClusterTPC,fCutMaxChi2PerClusterITS);
1078 printf("DCSToVertex2D = %d, MaxDCAToVertexXY = %2.2f, MaxDCAToVertexZ = %2.2f\n",fCutDCAToVertex2D,fCutMaxDCAToVertexXY,fCutMaxDCAToVertexZ);