additional DPs for new supermodules 2,3,6
[u/mrichter/AliRoot.git] / PWG / CaloTrackCorrBase / AliCalorimeterUtils.cxx
CommitLineData
765d44e7 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
765d44e7 15
16//_________________________________________________________________________
17// Class utility for Calorimeter specific selection methods ///
18//
19//
20//
21//-- Author: Gustavo Conesa (LPSC-Grenoble)
22//////////////////////////////////////////////////////////////////////////////
23
24
25// --- ROOT system ---
26#include "TGeoManager.h"
27
c5693f62 28// --- ANALYSIS system ---
765d44e7 29#include "AliCalorimeterUtils.h"
4b892846 30#include "AliESDEvent.h"
46a3cde6 31#include "AliMCEvent.h"
32#include "AliStack.h"
765d44e7 33#include "AliAODPWG4Particle.h"
c8fe2783 34#include "AliVCluster.h"
35#include "AliVCaloCells.h"
36#include "AliMixedEvent.h"
765d44e7 37
c5693f62 38// --- Detector ---
39#include "AliEMCALGeometry.h"
40#include "AliPHOSGeoUtils.h"
41
765d44e7 42ClassImp(AliCalorimeterUtils)
43
44
cb5780f4 45//____________________________________________
765d44e7 46 AliCalorimeterUtils::AliCalorimeterUtils() :
47 TObject(), fDebug(0),
90e32961 48 fEMCALGeoName("EMCAL_COMPLETE12SMV1"),
49 fPHOSGeoName ("PHOSgeo"),
9e8998b1 50 fEMCALGeo(0x0), fPHOSGeo(0x0),
51 fEMCALGeoMatrixSet(kFALSE), fPHOSGeoMatrixSet(kFALSE),
52 fLoadEMCALMatrices(kFALSE), fLoadPHOSMatrices(kFALSE),
53 fRemoveBadChannels(kFALSE), fPHOSBadChannelMap(0x0),
a5fb4114 54 fNCellsFromPHOSBorder(0),
9e8998b1 55 fNMaskCellColumns(0), fMaskCellColumns(0x0),
56 fRecalibration(kFALSE), fPHOSRecalibrationFactors(),
247abff4 57 fEMCALRecoUtils(new AliEMCALRecoUtils),
9e8998b1 58 fRecalculatePosition(kFALSE), fCorrectELinearity(kFALSE),
59 fRecalculateMatching(kFALSE),
60 fCutR(20), fCutZ(20),
7db7dcb6 61 fCutEta(20), fCutPhi(20),
62 fLocMaxCutE(0), fLocMaxCutEDiff(0)
765d44e7 63{
64 //Ctor
65
66 //Initialize parameters
67 InitParameters();
90e32961 68 for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
69 for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix [i] = 0 ;
7cc7d3f8 70
765d44e7 71}
765d44e7 72
cb5780f4 73//_________________________________________
74AliCalorimeterUtils::~AliCalorimeterUtils()
75{
765d44e7 76 //Dtor
77
4df35693 78 //if(fPHOSGeo) delete fPHOSGeo ;
765d44e7 79 if(fEMCALGeo) delete fEMCALGeo ;
7cc7d3f8 80
765d44e7 81 if(fPHOSBadChannelMap) {
82 fPHOSBadChannelMap->Clear();
83 delete fPHOSBadChannelMap;
84 }
85
09e819c9 86 if(fPHOSRecalibrationFactors) {
78219bac 87 fPHOSRecalibrationFactors->Clear();
88 delete fPHOSRecalibrationFactors;
09e819c9 89 }
90
a5fb4114 91 if(fEMCALRecoUtils) delete fEMCALRecoUtils ;
92 if(fNMaskCellColumns) delete [] fMaskCellColumns;
9584c261 93
765d44e7 94}
95
7db7dcb6 96//______________________________________________________________________________________
97Bool_t AliCalorimeterUtils::AreNeighbours(const TString calo,
98 const Int_t absId1, const Int_t absId2 ) const
99{
100 // Tells if (true) or not (false) two cells are neighbours
101 // A neighbour is defined as being two cells which share a side or corner
102
103 Bool_t areNeighbours = kFALSE ;
104
105 Int_t iRCU1 = -1, irow1 = -1, icol1 = -1;
106 Int_t iRCU2 = -1, irow2 = -1, icol2 = -1;
107
108 Int_t rowdiff = 0, coldiff = 0;
109
110 Int_t nSupMod1 = GetModuleNumberCellIndexes(absId1, calo, irow1, icol1, iRCU1);
111 Int_t nSupMod2 = GetModuleNumberCellIndexes(absId2, calo, irow2, icol2, iRCU2);
112
113 if(calo=="EMCAL" && nSupMod1!=nSupMod2)
114 {
115 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
116 // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
117 if(nSupMod1%2) icol1+=AliEMCALGeoParams::fgkEMCALCols;
118 else icol2+=AliEMCALGeoParams::fgkEMCALCols;
119 }
120
121 rowdiff = TMath::Abs( irow1 - irow2 ) ;
122 coldiff = TMath::Abs( icol1 - icol2 ) ;
123
124 if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
125 areNeighbours = kTRUE ;
126
127 return areNeighbours;
128
129}
130
131
cb5780f4 132//_____________________________________________________________________________________
133Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster,
134 AliVCaloCells* cells,
135 AliVEvent * event, Int_t iev) const
136{
7cc7d3f8 137
765d44e7 138 // Given the list of AbsId of the cluster, get the maximum cell and
139 // check if there are fNCellsFromBorder from the calorimeter border
140
7cc7d3f8 141 //If the distance to the border is 0 or negative just exit accept all clusters
247abff4 142 if(cells->GetType()==AliVCaloCells::kEMCALCell && fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder() <= 0 ) return kTRUE;
c8fe2783 143 if(cells->GetType()==AliVCaloCells::kPHOSCell && fNCellsFromPHOSBorder <= 0 ) return kTRUE;
7cc7d3f8 144
c8fe2783 145 Int_t absIdMax = -1;
765d44e7 146 Float_t ampMax = -1;
c8fe2783 147
148 AliMixedEvent * mixEvent = dynamic_cast<AliMixedEvent*> (event);
149 Int_t nMixedEvents = 0 ;
150 Int_t * cellsCumul = NULL ;
151 Int_t numberOfCells = 0 ;
152 if (mixEvent){
153 nMixedEvents = mixEvent->GetNumberOfEvents() ;
154 if (cells->GetType()==AliVCaloCells::kEMCALCell) {
155 cellsCumul = mixEvent->GetEMCALCellsCumul() ;
156 numberOfCells = mixEvent->GetNumberOfEMCALCells() ;
157 }
7cc7d3f8 158
c8fe2783 159 else if (cells->GetType()==AliVCaloCells::kPHOSCell) {
160 cellsCumul = mixEvent->GetPHOSCellsCumul() ;
161 numberOfCells = mixEvent->GetNumberOfPHOSCells() ;
162 }
898c9d44 163
164 if(cellsCumul){
165
166 Int_t startCell = cellsCumul[iev] ;
167 Int_t endCell = (iev+1 < nMixedEvents)?cellsCumul[iev+1]:numberOfCells;
c8fe2783 168 //Find cells with maximum amplitude
898c9d44 169 for(Int_t i = 0; i < cluster->GetNCells() ; i++){
170 Int_t absId = cluster->GetCellAbsId(i) ;
171 for (Int_t j = startCell; j < endCell ; j++) {
172 Short_t cellNumber;
173 Double_t amp ;
174 Double_t time;
175 cells->GetCell(j, cellNumber, amp, time) ;
176 if (absId == cellNumber) {
177 if(amp > ampMax){
178 ampMax = amp;
179 absIdMax = absId;
180 }
181 }
c8fe2783 182 }
898c9d44 183 }//loop on cluster cells
184 }// cells cumul available
185 else {
186 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - CellsCumul is NULL!!!\n");
187 abort();
c8fe2783 188 }
898c9d44 189 } else {//Normal SE Events
c8fe2783 190 for(Int_t i = 0; i < cluster->GetNCells() ; i++){
191 Int_t absId = cluster->GetCellAbsId(i) ;
192 Float_t amp = cells->GetCellAmplitude(absId);
193 if(amp > ampMax){
194 ampMax = amp;
195 absIdMax = absId;
196 }
197 }
198 }
765d44e7 199
200 if(fDebug > 1)
280e6766 201 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n",
7cc7d3f8 202 absIdMax, ampMax, cluster->E());
765d44e7 203
204 if(absIdMax==-1) return kFALSE;
205
206 //Check if the cell is close to the borders:
207 Bool_t okrow = kFALSE;
208 Bool_t okcol = kFALSE;
7cc7d3f8 209
c8fe2783 210 if(cells->GetType()==AliVCaloCells::kEMCALCell){
7cc7d3f8 211
765d44e7 212 Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1;
213 fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta);
214 fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
280e6766 215 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
216 Fatal("CheckCellFidutialRegion","Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",iSM,ieta,iphi);
217 }
7cc7d3f8 218
765d44e7 219 //Check rows/phi
247abff4 220 Int_t nborder = fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder();
765d44e7 221 if(iSM < 10){
247abff4 222 if(iphi >= nborder && iphi < 24-nborder) okrow =kTRUE;
7cc7d3f8 223 }
765d44e7 224 else{
247abff4 225 if(iphi >= nborder && iphi < 12-nborder) okrow =kTRUE;
765d44e7 226 }
227
247abff4 228 //Check columns/eta
229 if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0()){
230 if(ieta > nborder && ieta < 48-nborder) okcol =kTRUE;
765d44e7 231 }
232 else{
233 if(iSM%2==0){
247abff4 234 if(ieta >= nborder) okcol = kTRUE;
765d44e7 235 }
236 else {
247abff4 237 if(ieta < 48-nborder) okcol = kTRUE;
765d44e7 238 }
239 }//eta 0 not checked
240 if(fDebug > 1)
241 {
280e6766 242 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
7cc7d3f8 243 nborder, ieta, iphi, iSM);
765d44e7 244 if (okcol && okrow ) printf(" YES \n");
245 else printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
246 }
247 }//EMCAL
c8fe2783 248 else if(cells->GetType()==AliVCaloCells::kPHOSCell){
765d44e7 249 Int_t relId[4];
250 Int_t irow = -1, icol = -1;
251 fPHOSGeo->AbsToRelNumbering(absIdMax,relId);
252 irow = relId[2];
253 icol = relId[3];
254 //imod = relId[0]-1;
255 if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE;
256 if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE;
257 if(fDebug > 1)
258 {
280e6766 259 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - PHOS Cluster in %d cells fiducial volume: icol %d, irow %d, Module %d?",
7cc7d3f8 260 fNCellsFromPHOSBorder, icol, irow, relId[0]-1);
765d44e7 261 if (okcol && okrow ) printf(" YES \n");
262 else printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
263 }
264 }//PHOS
265
266 if (okcol && okrow) return kTRUE;
267 else return kFALSE;
268
269}
270
765d44e7 271//_________________________________________________________________________________________________________
cb5780f4 272Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter,UShort_t* cellList, Int_t nCells)
273{
765d44e7 274 // Check that in the cluster cells, there is no bad channel of those stored
275 // in fEMCALBadChannelMap or fPHOSBadChannelMap
276
277 if (!fRemoveBadChannels) return kFALSE;
36037088 278 //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
247abff4 279 if(calorimeter == "EMCAL" && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
78219bac 280 if(calorimeter == "PHOS" && !fPHOSBadChannelMap) return kFALSE;
7cc7d3f8 281
765d44e7 282 Int_t icol = -1;
283 Int_t irow = -1;
284 Int_t imod = -1;
285 for(Int_t iCell = 0; iCell<nCells; iCell++){
7cc7d3f8 286
765d44e7 287 //Get the column and row
288 if(calorimeter == "EMCAL"){
247abff4 289 return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
765d44e7 290 }
291 else if(calorimeter=="PHOS"){
292 Int_t relId[4];
293 fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
294 irow = relId[2];
295 icol = relId[3];
296 imod = relId[0]-1;
297 if(fPHOSBadChannelMap->GetEntries() <= imod)continue;
c5693f62 298 //printf("PHOS bad channels imod %d, icol %d, irow %d\n",imod, irow, icol);
299 if(GetPHOSChannelStatus(imod, irow, icol)) return kTRUE;
765d44e7 300 }
301 else return kFALSE;
302
303 }// cell cluster loop
7cc7d3f8 304
765d44e7 305 return kFALSE;
7cc7d3f8 306
765d44e7 307}
308
cb5780f4 309//_______________________________________________________________
310void AliCalorimeterUtils::CorrectClusterEnergy(AliVCluster *clus)
311{
9584c261 312 // Correct cluster energy non linearity
13cd2872 313
9584c261 314 clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
13cd2872 315
316}
317
c5693f62 318//________________________________________________________________________________________
319Int_t AliCalorimeterUtils::GetMaxEnergyCell(AliVCaloCells* cells, const AliVCluster* clu,
320 Float_t & clusterFraction) const
cb5780f4 321{
13cd2872 322
323 //For a given CaloCluster gets the absId of the cell
324 //with maximum energy deposit.
325
326 if( !clu || !cells ){
327 AliInfo("Cluster or cells pointer is null!");
328 return -1;
329 }
330
331 Double_t eMax =-1.;
332 Double_t eTot = 0.;
333 Double_t eCell =-1.;
334 Float_t fraction = 1.;
335 Float_t recalFactor = 1.;
336 Int_t cellAbsId =-1 , absId =-1 ;
337 Int_t iSupMod =-1 , ieta =-1 , iphi = -1, iRCU = -1;
338
339 TString calo = "EMCAL";
340 if(clu->IsPHOS()) calo = "PHOS";
341
342 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
343
344 cellAbsId = clu->GetCellAbsId(iDig);
345
346 fraction = clu->GetCellAmplitudeFraction(iDig);
347 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
348
b08dd6d5 349 iSupMod = GetModuleNumberCellIndexes(cellAbsId, calo, ieta, iphi, iRCU);
13cd2872 350
351 if(IsRecalibrationOn()) {
352 if(calo=="EMCAL") recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
c5693f62 353 else recalFactor = GetPHOSChannelRecalibrationFactor (iSupMod,iphi,ieta);
13cd2872 354 }
355
356 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
357
358 if(eCell > eMax) {
359 eMax = eCell;
360 absId = cellAbsId;
361 }
362
363 eTot+=eCell;
364
365 }// cell loop
366
367 clusterFraction = (eTot-eMax)/eTot; //Do not use cluster energy in case it was corrected for non linearity.
368
369 return absId;
370
9584c261 371}
372
d832b695 373//__________________________________________________________________________
374AliVTrack * AliCalorimeterUtils::GetMatchedTrack(const AliVCluster* cluster,
375 const AliVEvent* event,
376 const Int_t index) const
377{
378 // Get the matched track given its index, usually just the first match
379 // Since it is different for ESDs and AODs here it is a wrap method to do it
380
381 AliVTrack *track = 0;
382
383 // EMCAL case only when matching is recalculated
384 if(cluster->IsEMCAL() && IsRecalculationOfClusterTrackMatchingOn())
385 {
386 Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(cluster->GetID());
387 //printf("track index %d, cluster ID %d \n ",trackIndex,cluster->GetID());
388
389 if(trackIndex < 0 )
390 {
391 printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d, from recalculation\n", trackIndex);
392 }
393 else
394 {
395 track = dynamic_cast<AliVTrack*> (event->GetTrack(trackIndex));
396 }
397
398 return track ;
399
400 }
401
402 // Normal case, get info from ESD or AOD
403 // ESDs
404 if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
405 {
406 Int_t iESDtrack = cluster->GetTrackMatchedIndex();
407
408 if(iESDtrack < 0 )
409 {
410 printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d\n", index);
411 return 0x0;
412 }
413
414 track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
415
416 }
417 else // AODs
418 {
419 if(cluster->GetNTracksMatched() > 0 )
420 track = dynamic_cast<AliVTrack*>(cluster->GetTrackMatched(index));
421 }
422
423 return track ;
424
425}
426
cb5780f4 427//_____________________________________________________________________________________________________
765d44e7 428Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
429{
430 //Get the EMCAL/PHOS module number that corresponds to this particle
431
432 Int_t absId = -1;
433 if(particle->GetDetector()=="EMCAL"){
434 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
435 if(fDebug > 2)
436 printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
7cc7d3f8 437 particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
765d44e7 438 return fEMCALGeo->GetSuperModuleNumber(absId) ;
439 }//EMCAL
440 else if(particle->GetDetector()=="PHOS"){
46a3cde6 441 // In case we use the MC reader, the input are TParticles,
442 // in this case use the corresponing method in PHOS Geometry to get the particle.
443 if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
444 {
445 Int_t mod =-1;
446 Double_t z = 0., x=0.;
447 TParticle* primary = 0x0;
448 AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
449 if(stack) {
450 primary = stack->Particle(particle->GetCaloLabel(0));
451 }
452 else {
280e6766 453 Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
46a3cde6 454 }
898c9d44 455
456 if(primary){
457 fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
458 }
459 else{
460 Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
461 }
46a3cde6 462 return mod;
463 }
464 // Input are ESDs or AODs, get the PHOS module number like this.
465 else{
2206e918 466 //FIXME
467 //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
468 //return GetModuleNumber(cluster);
469 //MEFIX
470 return -1;
46a3cde6 471 }
765d44e7 472 }//PHOS
473
474 return -1;
475}
476
cb5780f4 477//_____________________________________________________________________
c8fe2783 478Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
765d44e7 479{
280e6766 480 //Get the EMCAL/PHOS module number that corresponds to this cluster
765d44e7 481 TLorentzVector lv;
482 Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
9cbbc28b 483 if(!cluster){
484 if(fDebug > 1) printf("AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
485 return -1;
486 }
765d44e7 487 cluster->GetMomentum(lv,v);
488 Float_t phi = lv.Phi();
489 if(phi < 0) phi+=TMath::TwoPi();
490 Int_t absId = -1;
c8fe2783 491 if(cluster->IsEMCAL()){
765d44e7 492 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
493 if(fDebug > 2)
280e6766 494 printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
7cc7d3f8 495 lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
765d44e7 496 return fEMCALGeo->GetSuperModuleNumber(absId) ;
497 }//EMCAL
c8fe2783 498 else if(cluster->IsPHOS()) {
765d44e7 499 Int_t relId[4];
500 if ( cluster->GetNCells() > 0) {
501 absId = cluster->GetCellAbsId(0);
502 if(fDebug > 2)
280e6766 503 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
7cc7d3f8 504 lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
765d44e7 505 }
506 else return -1;
507
508 if ( absId >= 0) {
509 fPHOSGeo->AbsToRelNumbering(absId,relId);
510 if(fDebug > 2)
280e6766 511 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
765d44e7 512 return relId[0]-1;
513 }
514 else return -1;
515 }//PHOS
516
517 return -1;
518}
519
cb5780f4 520//___________________________________________________________________________________________________
521Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(const Int_t absId, const TString calo,
522 Int_t & icol, Int_t & irow, Int_t & iRCU) const
765d44e7 523{
524 //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
525 Int_t imod = -1;
526 if ( absId >= 0) {
527 if(calo=="EMCAL"){
528 Int_t iTower = -1, iIphi = -1, iIeta = -1;
529 fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
530 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
280e6766 531 if(imod < 0 || irow < 0 || icol < 0 ) {
532 Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
533 }
534
765d44e7 535 //RCU0
536 if (0<=irow&&irow<8) iRCU=0; // first cable row
537 else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
538 //second cable row
539 //RCU1
540 else if(8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
541 //second cable row
542 else if(16<=irow&&irow<24) iRCU=1; // third cable row
543
544 if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
545 if (iRCU<0) {
280e6766 546 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
765d44e7 547 }
548
549 return imod ;
550 }//EMCAL
551 else{//PHOS
552 Int_t relId[4];
553 fPHOSGeo->AbsToRelNumbering(absId,relId);
554 irow = relId[2];
555 icol = relId[3];
556 imod = relId[0]-1;
557 iRCU= (Int_t)(relId[2]-1)/16 ;
558 //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
559 if (iRCU >= 4) {
280e6766 560 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
765d44e7 561 }
562 return imod;
563 }//PHOS
564 }
565
566 return -1;
567}
568
7db7dcb6 569//___________________________________________________________________________________________
570Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
571 Int_t *absIdList, Float_t *maxEList)
572{
573 // Find local maxima in cluster
574
575 Int_t iDigitN = 0 ;
576 Int_t iDigit = 0 ;
577 Int_t absId1 = -1 ;
578 Int_t absId2 = -1 ;
579 const Int_t nCells = cluster->GetNCells();
580
581 TString calorimeter = "EMCAL";
582 if(!cluster->IsEMCAL()) calorimeter = "PHOS";
583
584 //printf("cluster : ncells %d \n",nCells);
585 for(iDigit = 0; iDigit < nCells ; iDigit++){
586 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit] ;
587 /*
588 Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
589 RecalibrateCellAmplitude(en,calo,absIdList[iDigit]);
590 Int_t icol = -1, irow = -1, iRCU = -1;
591 Int_t sm = GetCaloUtils()->GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
592
593 printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
594 */
595 }
596
597
598 for(iDigit = 0 ; iDigit < nCells; iDigit++) {
599 if(absIdList[iDigit]>=0) {
600
601 absId1 = absIdList[iDigit] ;
602 //printf("%d : absID111 %d, %s\n",iDigit, absId1,calorimeter.Data());
603
604 Float_t en1 = cells->GetCellAmplitude(absId1);
605 RecalibrateCellAmplitude(en1,calorimeter,absId1);
606
607 for(iDigitN = 0; iDigitN < nCells; iDigitN++) {
608
609 absId2 = absIdList[iDigitN] ;
610
611 if(absId2==-1) continue;
612
613 //printf("\t %d : absID222 %d, %s\n",iDigitN, absId2,calorimeter.Data());
614
615 Float_t en2 = cells->GetCellAmplitude(absId2);
616 RecalibrateCellAmplitude(en2,calorimeter,absId2);
617
618 if ( AreNeighbours(calorimeter, absId1, absId2) ) {
619
620 if (en1 > en2 ) {
621 absIdList[iDigitN] = -1 ;
622 //printf("\t \t indexN %d not local max\n",iDigitN);
623 // but may be digit too is not local max ?
624 if(en1 < en2 + fLocMaxCutEDiff) {
625 //printf("\t \t index %d not local max cause locMaxCut\n",iDigit);
626 absIdList[iDigit] = -1 ;
627 }
628 }
629 else {
630 absIdList[iDigit] = -1 ;
631 //printf("\t \t index %d not local max\n",iDigitN);
632 // but may be digitN too is not local max ?
633 if(en1 > en2 - fLocMaxCutEDiff)
634 {
635 absIdList[iDigitN] = -1 ;
636 //printf("\t \t indexN %d not local max cause locMaxCut\n",iDigit);
637 }
638 }
639 } // if Areneighbours
640 } // while digitN
641 } // slot not empty
642 } // while digit
643
644 iDigitN = 0 ;
645 for(iDigit = 0; iDigit < nCells; iDigit++) {
646 if(absIdList[iDigit]>=0 ){
647 absIdList[iDigitN] = absIdList[iDigit] ;
648 Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
649 RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
650 if(en < fLocMaxCutE) continue; // Maxima only with seed energy at least
651 maxEList[iDigitN] = en ;
652 //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
653 iDigitN++ ;
654 }
655 }
656
657 //printf("N maxima %d \n",iDigitN);
658 //for(Int_t imax = 0; imax < iDigitN; imax++) printf("imax %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
659
660 return iDigitN ;
661
662}
663
cb5780f4 664//________________________________________
765d44e7 665void AliCalorimeterUtils::InitParameters()
666{
667 //Initialize the parameters of the analysis.
7db7dcb6 668
7cc7d3f8 669 fEMCALGeoName = "EMCAL_COMPLETEV1";
670 fPHOSGeoName = "PHOSgeo";
7db7dcb6 671
7cc7d3f8 672 fEMCALGeoMatrixSet = kFALSE;
673 fPHOSGeoMatrixSet = kFALSE;
7db7dcb6 674
7cc7d3f8 675 fRemoveBadChannels = kFALSE;
7db7dcb6 676
7cc7d3f8 677 fNCellsFromPHOSBorder = 0;
a5fb4114 678
7db7dcb6 679 fLocMaxCutE = 0.1 ;
680 fLocMaxCutEDiff = 0.0 ;
681
7cc7d3f8 682 // fMaskCellColumns = new Int_t[fNMaskCellColumns];
683 // fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
684 // fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
685 // fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
686 // fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
687 // fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
a5fb4114 688
765d44e7 689}
690
765d44e7 691
cb5780f4 692//_____________________________________________________
693void AliCalorimeterUtils::InitPHOSBadChannelStatusMap()
694{
765d44e7 695 //Init PHOS bad channels map
696 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
697 //In order to avoid rewriting the same histograms
698 Bool_t oldStatus = TH1::AddDirectoryStatus();
699 TH1::AddDirectory(kFALSE);
7cc7d3f8 700
78219bac 701 fPHOSBadChannelMap = new TObjArray(5);
c5693f62 702 for (int i = 0; i < 5; i++)fPHOSBadChannelMap->Add(new TH2I(Form("PHOS_BadMap_mod%d",i),Form("PHOS_BadMap_mod%d",i), 64, 0, 64, 56, 0, 56));
7cc7d3f8 703
765d44e7 704 fPHOSBadChannelMap->SetOwner(kTRUE);
705 fPHOSBadChannelMap->Compress();
706
707 //In order to avoid rewriting the same histograms
708 TH1::AddDirectory(oldStatus);
709}
710
cb5780f4 711//______________________________________________________
712void AliCalorimeterUtils::InitPHOSRecalibrationFactors()
713{
09e819c9 714 //Init EMCAL recalibration factors
715 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
716 //In order to avoid rewriting the same histograms
717 Bool_t oldStatus = TH1::AddDirectoryStatus();
718 TH1::AddDirectory(kFALSE);
7cc7d3f8 719
78219bac 720 fPHOSRecalibrationFactors = new TObjArray(5);
c5693f62 721 for (int i = 0; i < 5; i++)fPHOSRecalibrationFactors->Add(new TH2F(Form("PHOSRecalFactors_Mod%d",i),Form("PHOSRecalFactors_Mod%d",i), 64, 0, 64, 56, 0, 56));
09e819c9 722 //Init the histograms with 1
723 for (Int_t m = 0; m < 5; m++) {
724 for (Int_t i = 0; i < 56; i++) {
725 for (Int_t j = 0; j < 64; j++) {
c5693f62 726 SetPHOSChannelRecalibrationFactor(m,j,i,1.);
09e819c9 727 }
728 }
729 }
730 fPHOSRecalibrationFactors->SetOwner(kTRUE);
731 fPHOSRecalibrationFactors->Compress();
732
733 //In order to avoid rewriting the same histograms
734 TH1::AddDirectory(oldStatus);
735}
736
737
cb5780f4 738//___________________________________________
765d44e7 739void AliCalorimeterUtils::InitEMCALGeometry()
740{
741 //Initialize EMCAL geometry if it did not exist previously
742 if (!fEMCALGeo){
a38a48f2 743 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
765d44e7 744 if(fDebug > 0){
745 printf("AliCalorimeterUtils::InitEMCALGeometry()");
746 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
747 printf("\n");
748 }
749 }
750}
751
cb5780f4 752//__________________________________________
765d44e7 753void AliCalorimeterUtils::InitPHOSGeometry()
754{
755 //Initialize PHOS geometry if it did not exist previously
756 if (!fPHOSGeo){
757 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
758 if(fDebug > 0){
759 printf("AliCalorimeterUtils::InitPHOSGeometry()");
760 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
761 printf("\n");
762 }
763 }
764}
765
cb5780f4 766//_________________________________________________________
765d44e7 767void AliCalorimeterUtils::Print(const Option_t * opt) const
768{
7cc7d3f8 769
765d44e7 770 //Print some relevant parameters set for the analysis
771 if(! opt)
772 return;
7cc7d3f8 773
765d44e7 774 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
775 printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
776 printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
7cc7d3f8 777 fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
247abff4 778 if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
09e819c9 779 printf("Recalibrate Clusters? %d\n",fRecalibration);
9584c261 780 printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
781 printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
f2ccb5b8 782 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
7cc7d3f8 783
7db7dcb6 784 printf("Loc. Max. E > %2.2f\n", fLocMaxCutE);
785 printf("Loc. Max. E Diff > %2.2f\n", fLocMaxCutEDiff);
786
765d44e7 787 printf(" \n") ;
788}
789
cb5780f4 790//__________________________________________________________________
791Bool_t AliCalorimeterUtils::MaskFrameCluster(const Int_t iSM,
792 const Int_t ieta) const
793{
a5fb4114 794 //Check if cell is in one of the regions where we have significant amount
795 //of material in front. Only EMCAL
796
797 Int_t icol = ieta;
798 if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
799
800 if (fNMaskCellColumns && fMaskCellColumns) {
801 for (Int_t imask = 0; imask < fNMaskCellColumns; imask++) {
802 if(icol==fMaskCellColumns[imask]) return kTRUE;
803 }
804 }
805
806 return kFALSE;
807
808}
809
7db7dcb6 810//____________________________________________________________________________________
811void AliCalorimeterUtils::RecalibrateCellAmplitude(Float_t & amp,
812 const TString calo, const Int_t id)
813{
814 //Recaculate cell energy if recalibration factor
815
816 Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
817 Int_t nModule = GetModuleNumberCellIndexes(id,calo, icol, irow, iRCU);
818
819 if (IsRecalibrationOn())
820 {
821 if(calo == "PHOS")
822 {
823 amp *= GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
824 }
825 else
826 {
827 amp *= GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
828 }
829 }
830}
831
832//___________________________________________________________________________
833void AliCalorimeterUtils::RecalibrateCellTime(Double_t & time,
834 const TString calo,
835 const Int_t id, const Int_t bc)
836{
837 // Recalculate time if time recalibration available for EMCAL
838 // not ready for PHOS
839
840 if(calo == "EMCAL" && GetEMCALRecoUtils()->IsTimeRecalibrationOn())
841 {
842 GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
843 }
844
845}
846
a5fb4114 847
cb5780f4 848//__________________________________________________________________________
849Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster,
850 AliVCaloCells * cells)
851{
09e819c9 852 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
7cc7d3f8 853
898c9d44 854 //Initialize some used variables
7db7dcb6 855 Float_t frac = 0., energy = 0.;
898c9d44 856
7db7dcb6 857 if(cells)
858 {
7cc7d3f8 859 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
7db7dcb6 860
7cc7d3f8 861 UShort_t * index = cluster->GetCellsAbsId() ;
862 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
7db7dcb6 863
7cc7d3f8 864 Int_t ncells = cluster->GetNCells();
7db7dcb6 865
7cc7d3f8 866 TString calo = "EMCAL";
7db7dcb6 867 if(cluster->IsPHOS())
868 calo = "PHOS";
7cc7d3f8 869
870 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
871 for(Int_t icell = 0; icell < ncells; icell++){
7db7dcb6 872
873 Int_t absId = index[icell];
874
7cc7d3f8 875 frac = fraction[icell];
876 if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
7db7dcb6 877
878 Float_t amp = cells->GetCellAmplitude(absId);
879 RecalibrateCellAmplitude(amp,calo, absId);
880
7cc7d3f8 881 if(fDebug>2)
7db7dcb6 882 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, cell fraction %f, cell energy %f\n",
883 calo.Data(),frac,cells->GetCellAmplitude(absId));
7cc7d3f8 884
7db7dcb6 885 energy += amp*frac;
7cc7d3f8 886 }
887
888 if(fDebug>1)
889 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
898c9d44 890
891 }// cells available
7db7dcb6 892 else
893 {
898c9d44 894 Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
895 }
896
09e819c9 897 return energy;
09e819c9 898}
899
cb5780f4 900//________________________________________________________________________________
765d44e7 901void AliCalorimeterUtils::SetGeometryTransformationMatrices(AliVEvent* inputEvent)
902{
903 //Set the calorimeters transformation matrices
7cc7d3f8 904
765d44e7 905 //Get the EMCAL transformation geometry matrices from ESD
3b13c34c 906 if(!fEMCALGeoMatrixSet && fEMCALGeo){
907 if(fLoadEMCALMatrices){
c5693f62 908 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined EMCAL geometry matrices\n");
3b13c34c 909 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
910 if(fEMCALMatrix[mod]){
911 if(fDebug > 1)
912 fEMCALMatrix[mod]->Print();
913 fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;
914 }
915 }//SM loop
916 fEMCALGeoMatrixSet = kTRUE;//At least one, so good
917
918 }//Load matrices
919 else if (!gGeoManager) {
7cc7d3f8 920
3b13c34c 921 if(fDebug > 1)
922 printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load EMCAL misalignment matrices. \n");
923 if(!strcmp(inputEvent->GetName(),"AliESDEvent")) {
924 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
7cc7d3f8 925 //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
4b892846 926 if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod)) {
927 fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
3b13c34c 928 }
929 }// loop over super modules
c5693f62 930
3b13c34c 931 fEMCALGeoMatrixSet = kTRUE;//At least one, so good
932
933 }//ESD as input
934 else {
935 if(fDebug > 1)
936 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
937 }//AOD as input
938 }//Get matrix from data
04131edb 939 else if(gGeoManager){
940 fEMCALGeoMatrixSet = kTRUE;
941 }
765d44e7 942 }//EMCAL geo && no geoManager
c5693f62 943
765d44e7 944 //Get the PHOS transformation geometry matrices from ESD
3b13c34c 945 if(!fPHOSGeoMatrixSet && fPHOSGeo){
c5693f62 946
3b13c34c 947 if(fLoadPHOSMatrices){
c5693f62 948 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined PHOS geometry matrices\n");
949 for(Int_t mod = 0 ; mod < 5 ; mod++){
3b13c34c 950 if(fPHOSMatrix[mod]){
c5693f62 951 if(fDebug > 1 )
3b13c34c 952 fPHOSMatrix[mod]->Print();
dbf54f1e 953 fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;
3b13c34c 954 }
955 }//SM loop
956 fPHOSGeoMatrixSet = kTRUE;//At least one, so good
c5693f62 957
3b13c34c 958 }//Load matrices
959 else if (!gGeoManager) {
960 if(fDebug > 1)
961 printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
765d44e7 962 if(!strcmp(inputEvent->GetName(),"AliESDEvent")) {
dbf54f1e 963 for(Int_t mod = 0; mod < 5; mod++){
4b892846 964 if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod)) {
765d44e7 965 //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
4b892846 966 fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
765d44e7 967 }
3b13c34c 968 }// loop over modules
969 fPHOSGeoMatrixSet = kTRUE; //At least one so good
765d44e7 970 }//ESD as input
971 else {
972 if(fDebug > 1)
973 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
3b13c34c 974 }//AOD as input
975 }// get matrix from data
04131edb 976 else if(gGeoManager){
977 fPHOSGeoMatrixSet = kTRUE;
978 }
765d44e7 979 }//PHOS geo and geoManager was not set
c5693f62 980
765d44e7 981}
982
cb5780f4 983//__________________________________________________________________________________________
984void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
985{
9584c261 986
987 //Recalculate EMCAL cluster position
988
19db8f8c 989 fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
7cc7d3f8 990
9584c261 991}
cb5780f4 992
993//________________________________________________________________________________
994void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event,
995 TObjArray* clusterArray)
996{
997 //Recalculate track matching
998
999 if (fRecalculateMatching) {
1000 fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo) ;
1001 //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1002 //if(esdevent){
1003 // fEMCALRecoUtils->SetClusterMatchedToTrack(esdevent) ;
1004 // fEMCALRecoUtils->SetTracksMatchedToCluster(esdevent) ;
1005 //}
1006 }
1007}