Consistent pT binning for Xi's
[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 ---
3c1d9afb 26#include <TGeoManager.h>
27#include <TH2F.h>
28#include <TCanvas.h>
29#include <TStyle.h>
30#include <TPad.h>
765d44e7 31
c5693f62 32// --- ANALYSIS system ---
765d44e7 33#include "AliCalorimeterUtils.h"
4b892846 34#include "AliESDEvent.h"
46a3cde6 35#include "AliMCEvent.h"
36#include "AliStack.h"
765d44e7 37#include "AliAODPWG4Particle.h"
c8fe2783 38#include "AliVCluster.h"
39#include "AliVCaloCells.h"
40#include "AliMixedEvent.h"
3c1d9afb 41#include "AliAODCaloCluster.h"
765d44e7 42
c5693f62 43// --- Detector ---
44#include "AliEMCALGeometry.h"
45#include "AliPHOSGeoUtils.h"
46
765d44e7 47ClassImp(AliCalorimeterUtils)
48
49
cb5780f4 50//____________________________________________
765d44e7 51 AliCalorimeterUtils::AliCalorimeterUtils() :
52 TObject(), fDebug(0),
90e32961 53 fEMCALGeoName("EMCAL_COMPLETE12SMV1"),
54 fPHOSGeoName ("PHOSgeo"),
9e8998b1 55 fEMCALGeo(0x0), fPHOSGeo(0x0),
56 fEMCALGeoMatrixSet(kFALSE), fPHOSGeoMatrixSet(kFALSE),
57 fLoadEMCALMatrices(kFALSE), fLoadPHOSMatrices(kFALSE),
58 fRemoveBadChannels(kFALSE), fPHOSBadChannelMap(0x0),
a5fb4114 59 fNCellsFromPHOSBorder(0),
9e8998b1 60 fNMaskCellColumns(0), fMaskCellColumns(0x0),
61 fRecalibration(kFALSE), fPHOSRecalibrationFactors(),
247abff4 62 fEMCALRecoUtils(new AliEMCALRecoUtils),
9e8998b1 63 fRecalculatePosition(kFALSE), fCorrectELinearity(kFALSE),
64 fRecalculateMatching(kFALSE),
65 fCutR(20), fCutZ(20),
7db7dcb6 66 fCutEta(20), fCutPhi(20),
3c1d9afb 67 fLocMaxCutE(0), fLocMaxCutEDiff(0),
68 fPlotCluster(0)
765d44e7 69{
70 //Ctor
71
72 //Initialize parameters
73 InitParameters();
90e32961 74 for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
75 for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix [i] = 0 ;
7cc7d3f8 76
765d44e7 77}
765d44e7 78
cb5780f4 79//_________________________________________
80AliCalorimeterUtils::~AliCalorimeterUtils()
81{
765d44e7 82 //Dtor
83
4df35693 84 //if(fPHOSGeo) delete fPHOSGeo ;
765d44e7 85 if(fEMCALGeo) delete fEMCALGeo ;
7cc7d3f8 86
765d44e7 87 if(fPHOSBadChannelMap) {
88 fPHOSBadChannelMap->Clear();
89 delete fPHOSBadChannelMap;
90 }
91
09e819c9 92 if(fPHOSRecalibrationFactors) {
78219bac 93 fPHOSRecalibrationFactors->Clear();
94 delete fPHOSRecalibrationFactors;
09e819c9 95 }
96
a5fb4114 97 if(fEMCALRecoUtils) delete fEMCALRecoUtils ;
98 if(fNMaskCellColumns) delete [] fMaskCellColumns;
9584c261 99
765d44e7 100}
101
7db7dcb6 102//______________________________________________________________________________________
103Bool_t AliCalorimeterUtils::AreNeighbours(const TString calo,
104 const Int_t absId1, const Int_t absId2 ) const
105{
106 // Tells if (true) or not (false) two cells are neighbours
107 // A neighbour is defined as being two cells which share a side or corner
108
109 Bool_t areNeighbours = kFALSE ;
110
111 Int_t iRCU1 = -1, irow1 = -1, icol1 = -1;
112 Int_t iRCU2 = -1, irow2 = -1, icol2 = -1;
113
114 Int_t rowdiff = 0, coldiff = 0;
115
9369a2b1 116 Int_t nSupMod1 = GetModuleNumberCellIndexes(absId1, calo, icol1, irow1, iRCU1);
117 Int_t nSupMod2 = GetModuleNumberCellIndexes(absId2, calo, icol2, irow2, iRCU2);
7db7dcb6 118
119 if(calo=="EMCAL" && nSupMod1!=nSupMod2)
120 {
121 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
122 // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
123 if(nSupMod1%2) icol1+=AliEMCALGeoParams::fgkEMCALCols;
9369a2b1 124 else icol2+=AliEMCALGeoParams::fgkEMCALCols;
7db7dcb6 125 }
126
127 rowdiff = TMath::Abs( irow1 - irow2 ) ;
128 coldiff = TMath::Abs( icol1 - icol2 ) ;
129
130 if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
131 areNeighbours = kTRUE ;
132
133 return areNeighbours;
134
135}
136
137
cb5780f4 138//_____________________________________________________________________________________
139Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster,
140 AliVCaloCells* cells,
141 AliVEvent * event, Int_t iev) const
142{
7cc7d3f8 143
765d44e7 144 // Given the list of AbsId of the cluster, get the maximum cell and
145 // check if there are fNCellsFromBorder from the calorimeter border
146
7cc7d3f8 147 //If the distance to the border is 0 or negative just exit accept all clusters
247abff4 148 if(cells->GetType()==AliVCaloCells::kEMCALCell && fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder() <= 0 ) return kTRUE;
c8fe2783 149 if(cells->GetType()==AliVCaloCells::kPHOSCell && fNCellsFromPHOSBorder <= 0 ) return kTRUE;
7cc7d3f8 150
c8fe2783 151 Int_t absIdMax = -1;
765d44e7 152 Float_t ampMax = -1;
c8fe2783 153
154 AliMixedEvent * mixEvent = dynamic_cast<AliMixedEvent*> (event);
155 Int_t nMixedEvents = 0 ;
156 Int_t * cellsCumul = NULL ;
157 Int_t numberOfCells = 0 ;
158 if (mixEvent){
159 nMixedEvents = mixEvent->GetNumberOfEvents() ;
160 if (cells->GetType()==AliVCaloCells::kEMCALCell) {
161 cellsCumul = mixEvent->GetEMCALCellsCumul() ;
162 numberOfCells = mixEvent->GetNumberOfEMCALCells() ;
163 }
7cc7d3f8 164
c8fe2783 165 else if (cells->GetType()==AliVCaloCells::kPHOSCell) {
166 cellsCumul = mixEvent->GetPHOSCellsCumul() ;
167 numberOfCells = mixEvent->GetNumberOfPHOSCells() ;
168 }
898c9d44 169
170 if(cellsCumul){
171
172 Int_t startCell = cellsCumul[iev] ;
173 Int_t endCell = (iev+1 < nMixedEvents)?cellsCumul[iev+1]:numberOfCells;
c8fe2783 174 //Find cells with maximum amplitude
898c9d44 175 for(Int_t i = 0; i < cluster->GetNCells() ; i++){
176 Int_t absId = cluster->GetCellAbsId(i) ;
177 for (Int_t j = startCell; j < endCell ; j++) {
178 Short_t cellNumber;
179 Double_t amp ;
180 Double_t time;
181 cells->GetCell(j, cellNumber, amp, time) ;
182 if (absId == cellNumber) {
183 if(amp > ampMax){
184 ampMax = amp;
185 absIdMax = absId;
186 }
187 }
c8fe2783 188 }
898c9d44 189 }//loop on cluster cells
190 }// cells cumul available
191 else {
192 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - CellsCumul is NULL!!!\n");
193 abort();
c8fe2783 194 }
898c9d44 195 } else {//Normal SE Events
c8fe2783 196 for(Int_t i = 0; i < cluster->GetNCells() ; i++){
197 Int_t absId = cluster->GetCellAbsId(i) ;
198 Float_t amp = cells->GetCellAmplitude(absId);
199 if(amp > ampMax){
200 ampMax = amp;
201 absIdMax = absId;
202 }
203 }
204 }
765d44e7 205
206 if(fDebug > 1)
280e6766 207 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n",
7cc7d3f8 208 absIdMax, ampMax, cluster->E());
765d44e7 209
210 if(absIdMax==-1) return kFALSE;
211
212 //Check if the cell is close to the borders:
213 Bool_t okrow = kFALSE;
214 Bool_t okcol = kFALSE;
7cc7d3f8 215
c8fe2783 216 if(cells->GetType()==AliVCaloCells::kEMCALCell){
7cc7d3f8 217
765d44e7 218 Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1;
219 fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta);
220 fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
280e6766 221 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
222 Fatal("CheckCellFidutialRegion","Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",iSM,ieta,iphi);
223 }
7cc7d3f8 224
765d44e7 225 //Check rows/phi
247abff4 226 Int_t nborder = fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder();
765d44e7 227 if(iSM < 10){
247abff4 228 if(iphi >= nborder && iphi < 24-nborder) okrow =kTRUE;
7cc7d3f8 229 }
765d44e7 230 else{
247abff4 231 if(iphi >= nborder && iphi < 12-nborder) okrow =kTRUE;
765d44e7 232 }
233
247abff4 234 //Check columns/eta
235 if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0()){
236 if(ieta > nborder && ieta < 48-nborder) okcol =kTRUE;
765d44e7 237 }
238 else{
239 if(iSM%2==0){
247abff4 240 if(ieta >= nborder) okcol = kTRUE;
765d44e7 241 }
242 else {
247abff4 243 if(ieta < 48-nborder) okcol = kTRUE;
765d44e7 244 }
245 }//eta 0 not checked
246 if(fDebug > 1)
247 {
280e6766 248 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
7cc7d3f8 249 nborder, ieta, iphi, iSM);
765d44e7 250 if (okcol && okrow ) printf(" YES \n");
251 else printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
252 }
253 }//EMCAL
c8fe2783 254 else if(cells->GetType()==AliVCaloCells::kPHOSCell){
765d44e7 255 Int_t relId[4];
256 Int_t irow = -1, icol = -1;
257 fPHOSGeo->AbsToRelNumbering(absIdMax,relId);
258 irow = relId[2];
259 icol = relId[3];
260 //imod = relId[0]-1;
261 if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE;
262 if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE;
263 if(fDebug > 1)
264 {
280e6766 265 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - PHOS Cluster in %d cells fiducial volume: icol %d, irow %d, Module %d?",
7cc7d3f8 266 fNCellsFromPHOSBorder, icol, irow, relId[0]-1);
765d44e7 267 if (okcol && okrow ) printf(" YES \n");
268 else printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
269 }
270 }//PHOS
271
272 if (okcol && okrow) return kTRUE;
273 else return kFALSE;
274
275}
276
765d44e7 277//_________________________________________________________________________________________________________
cb5780f4 278Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter,UShort_t* cellList, Int_t nCells)
279{
765d44e7 280 // Check that in the cluster cells, there is no bad channel of those stored
281 // in fEMCALBadChannelMap or fPHOSBadChannelMap
282
283 if (!fRemoveBadChannels) return kFALSE;
36037088 284 //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
247abff4 285 if(calorimeter == "EMCAL" && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
78219bac 286 if(calorimeter == "PHOS" && !fPHOSBadChannelMap) return kFALSE;
7cc7d3f8 287
765d44e7 288 Int_t icol = -1;
289 Int_t irow = -1;
290 Int_t imod = -1;
291 for(Int_t iCell = 0; iCell<nCells; iCell++){
7cc7d3f8 292
765d44e7 293 //Get the column and row
294 if(calorimeter == "EMCAL"){
247abff4 295 return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
765d44e7 296 }
297 else if(calorimeter=="PHOS"){
298 Int_t relId[4];
299 fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
300 irow = relId[2];
301 icol = relId[3];
302 imod = relId[0]-1;
303 if(fPHOSBadChannelMap->GetEntries() <= imod)continue;
c5693f62 304 //printf("PHOS bad channels imod %d, icol %d, irow %d\n",imod, irow, icol);
305 if(GetPHOSChannelStatus(imod, irow, icol)) return kTRUE;
765d44e7 306 }
307 else return kFALSE;
308
309 }// cell cluster loop
7cc7d3f8 310
765d44e7 311 return kFALSE;
7cc7d3f8 312
765d44e7 313}
314
cb5780f4 315//_______________________________________________________________
316void AliCalorimeterUtils::CorrectClusterEnergy(AliVCluster *clus)
317{
9584c261 318 // Correct cluster energy non linearity
13cd2872 319
9584c261 320 clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
13cd2872 321
322}
323
c5693f62 324//________________________________________________________________________________________
325Int_t AliCalorimeterUtils::GetMaxEnergyCell(AliVCaloCells* cells, const AliVCluster* clu,
326 Float_t & clusterFraction) const
cb5780f4 327{
13cd2872 328
329 //For a given CaloCluster gets the absId of the cell
330 //with maximum energy deposit.
331
332 if( !clu || !cells ){
333 AliInfo("Cluster or cells pointer is null!");
334 return -1;
335 }
336
337 Double_t eMax =-1.;
338 Double_t eTot = 0.;
339 Double_t eCell =-1.;
340 Float_t fraction = 1.;
341 Float_t recalFactor = 1.;
342 Int_t cellAbsId =-1 , absId =-1 ;
343 Int_t iSupMod =-1 , ieta =-1 , iphi = -1, iRCU = -1;
344
345 TString calo = "EMCAL";
346 if(clu->IsPHOS()) calo = "PHOS";
347
348 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
349
350 cellAbsId = clu->GetCellAbsId(iDig);
351
352 fraction = clu->GetCellAmplitudeFraction(iDig);
353 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
354
b08dd6d5 355 iSupMod = GetModuleNumberCellIndexes(cellAbsId, calo, ieta, iphi, iRCU);
13cd2872 356
357 if(IsRecalibrationOn()) {
358 if(calo=="EMCAL") recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
c5693f62 359 else recalFactor = GetPHOSChannelRecalibrationFactor (iSupMod,iphi,ieta);
13cd2872 360 }
361
362 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
363
364 if(eCell > eMax) {
365 eMax = eCell;
366 absId = cellAbsId;
367 }
368
369 eTot+=eCell;
370
371 }// cell loop
372
373 clusterFraction = (eTot-eMax)/eTot; //Do not use cluster energy in case it was corrected for non linearity.
374
375 return absId;
376
9584c261 377}
378
d832b695 379//__________________________________________________________________________
380AliVTrack * AliCalorimeterUtils::GetMatchedTrack(const AliVCluster* cluster,
381 const AliVEvent* event,
382 const Int_t index) const
383{
384 // Get the matched track given its index, usually just the first match
385 // Since it is different for ESDs and AODs here it is a wrap method to do it
386
387 AliVTrack *track = 0;
388
389 // EMCAL case only when matching is recalculated
390 if(cluster->IsEMCAL() && IsRecalculationOfClusterTrackMatchingOn())
391 {
392 Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(cluster->GetID());
393 //printf("track index %d, cluster ID %d \n ",trackIndex,cluster->GetID());
394
395 if(trackIndex < 0 )
396 {
397 printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d, from recalculation\n", trackIndex);
398 }
399 else
400 {
401 track = dynamic_cast<AliVTrack*> (event->GetTrack(trackIndex));
402 }
403
404 return track ;
405
406 }
407
408 // Normal case, get info from ESD or AOD
409 // ESDs
410 if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
411 {
412 Int_t iESDtrack = cluster->GetTrackMatchedIndex();
413
414 if(iESDtrack < 0 )
415 {
416 printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d\n", index);
417 return 0x0;
418 }
419
420 track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
421
422 }
423 else // AODs
424 {
425 if(cluster->GetNTracksMatched() > 0 )
426 track = dynamic_cast<AliVTrack*>(cluster->GetTrackMatched(index));
427 }
428
429 return track ;
430
431}
432
cb5780f4 433//_____________________________________________________________________________________________________
765d44e7 434Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
435{
436 //Get the EMCAL/PHOS module number that corresponds to this particle
437
438 Int_t absId = -1;
439 if(particle->GetDetector()=="EMCAL"){
440 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
441 if(fDebug > 2)
442 printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
7cc7d3f8 443 particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
765d44e7 444 return fEMCALGeo->GetSuperModuleNumber(absId) ;
445 }//EMCAL
446 else if(particle->GetDetector()=="PHOS"){
46a3cde6 447 // In case we use the MC reader, the input are TParticles,
448 // in this case use the corresponing method in PHOS Geometry to get the particle.
449 if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
450 {
451 Int_t mod =-1;
452 Double_t z = 0., x=0.;
453 TParticle* primary = 0x0;
454 AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
455 if(stack) {
456 primary = stack->Particle(particle->GetCaloLabel(0));
457 }
458 else {
280e6766 459 Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
46a3cde6 460 }
898c9d44 461
462 if(primary){
463 fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
464 }
465 else{
466 Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
467 }
46a3cde6 468 return mod;
469 }
470 // Input are ESDs or AODs, get the PHOS module number like this.
471 else{
2206e918 472 //FIXME
473 //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
474 //return GetModuleNumber(cluster);
475 //MEFIX
476 return -1;
46a3cde6 477 }
765d44e7 478 }//PHOS
479
480 return -1;
481}
482
cb5780f4 483//_____________________________________________________________________
c8fe2783 484Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
765d44e7 485{
280e6766 486 //Get the EMCAL/PHOS module number that corresponds to this cluster
765d44e7 487 TLorentzVector lv;
488 Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
9cbbc28b 489 if(!cluster){
490 if(fDebug > 1) printf("AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
491 return -1;
492 }
765d44e7 493 cluster->GetMomentum(lv,v);
494 Float_t phi = lv.Phi();
495 if(phi < 0) phi+=TMath::TwoPi();
496 Int_t absId = -1;
c8fe2783 497 if(cluster->IsEMCAL()){
765d44e7 498 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
499 if(fDebug > 2)
280e6766 500 printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
7cc7d3f8 501 lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
765d44e7 502 return fEMCALGeo->GetSuperModuleNumber(absId) ;
503 }//EMCAL
c8fe2783 504 else if(cluster->IsPHOS()) {
765d44e7 505 Int_t relId[4];
506 if ( cluster->GetNCells() > 0) {
507 absId = cluster->GetCellAbsId(0);
508 if(fDebug > 2)
280e6766 509 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
7cc7d3f8 510 lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
765d44e7 511 }
512 else return -1;
513
514 if ( absId >= 0) {
515 fPHOSGeo->AbsToRelNumbering(absId,relId);
516 if(fDebug > 2)
280e6766 517 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
765d44e7 518 return relId[0]-1;
519 }
520 else return -1;
521 }//PHOS
522
523 return -1;
524}
525
cb5780f4 526//___________________________________________________________________________________________________
527Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(const Int_t absId, const TString calo,
528 Int_t & icol, Int_t & irow, Int_t & iRCU) const
765d44e7 529{
530 //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
531 Int_t imod = -1;
532 if ( absId >= 0) {
533 if(calo=="EMCAL"){
534 Int_t iTower = -1, iIphi = -1, iIeta = -1;
535 fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
536 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
280e6766 537 if(imod < 0 || irow < 0 || icol < 0 ) {
538 Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
539 }
540
765d44e7 541 //RCU0
542 if (0<=irow&&irow<8) iRCU=0; // first cable row
543 else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
544 //second cable row
545 //RCU1
546 else if(8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
547 //second cable row
548 else if(16<=irow&&irow<24) iRCU=1; // third cable row
549
550 if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
551 if (iRCU<0) {
280e6766 552 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
765d44e7 553 }
554
555 return imod ;
556 }//EMCAL
557 else{//PHOS
558 Int_t relId[4];
559 fPHOSGeo->AbsToRelNumbering(absId,relId);
560 irow = relId[2];
561 icol = relId[3];
562 imod = relId[0]-1;
563 iRCU= (Int_t)(relId[2]-1)/16 ;
564 //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
565 if (iRCU >= 4) {
280e6766 566 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
765d44e7 567 }
568 return imod;
569 }//PHOS
570 }
571
572 return -1;
573}
574
7db7dcb6 575//___________________________________________________________________________________________
71e3889f 576Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells)
577{
578 // Find local maxima in cluster
579
580 const Int_t nc = cluster->GetNCells();
581 Int_t *absIdList = new Int_t [nc];
582 Float_t *maxEList = new Float_t[nc];
583
584 Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList);
585
586 delete [] absIdList;
587 delete [] maxEList;
588
589 return nMax;
590
591}
592
593//___________________________________________________________________________________________
7db7dcb6 594Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
595 Int_t *absIdList, Float_t *maxEList)
596{
597 // Find local maxima in cluster
f241ecc3 598
7db7dcb6 599 Int_t iDigitN = 0 ;
600 Int_t iDigit = 0 ;
601 Int_t absId1 = -1 ;
602 Int_t absId2 = -1 ;
603 const Int_t nCells = cluster->GetNCells();
604
605 TString calorimeter = "EMCAL";
606 if(!cluster->IsEMCAL()) calorimeter = "PHOS";
607
608 //printf("cluster : ncells %d \n",nCells);
f241ecc3 609
610 Float_t emax = 0;
611 Int_t idmax =-1;
9369a2b1 612 for(iDigit = 0; iDigit < nCells ; iDigit++)
613 {
7db7dcb6 614 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit] ;
f241ecc3 615 Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
616 RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
617 if( en > emax )
618 {
619 emax = en ;
620 idmax = absIdList[iDigit] ;
621 }
622 //Int_t icol = -1, irow = -1, iRCU = -1;
623 //Int_t sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
624 //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
7db7dcb6 625 }
626
9369a2b1 627 for(iDigit = 0 ; iDigit < nCells; iDigit++)
628 {
629 if(absIdList[iDigit]>=0)
630 {
631 absId1 = cluster->GetCellsAbsId()[iDigit];
7db7dcb6 632
633 Float_t en1 = cells->GetCellAmplitude(absId1);
634 RecalibrateCellAmplitude(en1,calorimeter,absId1);
f241ecc3 635
9369a2b1 636 //printf("%d : absIDi %d, E %f\n",iDigit, absId1,en1);
7db7dcb6 637
9369a2b1 638 for(iDigitN = 0; iDigitN < nCells; iDigitN++)
639 {
640 absId2 = cluster->GetCellsAbsId()[iDigitN] ;
7db7dcb6 641
9369a2b1 642 if(absId2==-1 || absId2==absId1) continue;
7db7dcb6 643
9369a2b1 644 //printf("\t %d : absIDj %d\n",iDigitN, absId2);
7db7dcb6 645
646 Float_t en2 = cells->GetCellAmplitude(absId2);
647 RecalibrateCellAmplitude(en2,calorimeter,absId2);
f241ecc3 648
9369a2b1 649 //printf("\t %d : absIDj %d, E %f\n",iDigitN, absId2,en2);
7db7dcb6 650
9369a2b1 651 if ( AreNeighbours(calorimeter, absId1, absId2) )
652 {
653 // printf("\t \t Neighbours \n");
654 if (en1 > en2 )
655 {
7db7dcb6 656 absIdList[iDigitN] = -1 ;
657 //printf("\t \t indexN %d not local max\n",iDigitN);
658 // but may be digit too is not local max ?
659 if(en1 < en2 + fLocMaxCutEDiff) {
9369a2b1 660 //printf("\t \t index %d not local max cause locMaxCutEDiff\n",iDigit);
7db7dcb6 661 absIdList[iDigit] = -1 ;
662 }
663 }
9369a2b1 664 else
665 {
7db7dcb6 666 absIdList[iDigit] = -1 ;
667 //printf("\t \t index %d not local max\n",iDigitN);
668 // but may be digitN too is not local max ?
669 if(en1 > en2 - fLocMaxCutEDiff)
670 {
671 absIdList[iDigitN] = -1 ;
9369a2b1 672 //printf("\t \t indexN %d not local max cause locMaxCutEDiff\n",iDigit);
7db7dcb6 673 }
674 }
9369a2b1 675 } // if Are neighbours
676 //else printf("\t \t NOT Neighbours \n");
7db7dcb6 677 } // while digitN
678 } // slot not empty
679 } // while digit
680
681 iDigitN = 0 ;
9369a2b1 682 for(iDigit = 0; iDigit < nCells; iDigit++)
683 {
684 if(absIdList[iDigit]>=0 )
685 {
7db7dcb6 686 absIdList[iDigitN] = absIdList[iDigit] ;
687 Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
688 RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
689 if(en < fLocMaxCutE) continue; // Maxima only with seed energy at least
690 maxEList[iDigitN] = en ;
691 //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
692 iDigitN++ ;
693 }
694 }
695
f241ecc3 696 if(iDigitN == 0)
697 {
698 if(fDebug > 0)
699 printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - No local maxima found, assign highest energy cell as maxima, id %d, en cell %2.2f, en cluster %2.2f\n",
700 idmax,emax,cluster->E());
701 iDigitN = 1 ;
702 maxEList[0] = emax ;
703 absIdList[0] = idmax ;
704 }
705
706 if(fDebug > 0)
707 {
708 printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - In cluster E %2.2f, M02 %2.2f, M20 %2.2f, N maxima %d \n",
709 cluster->E(),cluster->GetM02(),cluster->GetM20(), iDigitN);
710
711 if(fDebug > 1) for(Int_t imax = 0; imax < iDigitN; imax++)
712 {
713 printf(" \t i %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
714 }
715 }
7db7dcb6 716
717 return iDigitN ;
718
719}
720
cb5780f4 721//________________________________________
765d44e7 722void AliCalorimeterUtils::InitParameters()
723{
724 //Initialize the parameters of the analysis.
7db7dcb6 725
7cc7d3f8 726 fEMCALGeoName = "EMCAL_COMPLETEV1";
727 fPHOSGeoName = "PHOSgeo";
7db7dcb6 728
7cc7d3f8 729 fEMCALGeoMatrixSet = kFALSE;
730 fPHOSGeoMatrixSet = kFALSE;
7db7dcb6 731
7cc7d3f8 732 fRemoveBadChannels = kFALSE;
7db7dcb6 733
7cc7d3f8 734 fNCellsFromPHOSBorder = 0;
a5fb4114 735
7db7dcb6 736 fLocMaxCutE = 0.1 ;
737 fLocMaxCutEDiff = 0.0 ;
738
7cc7d3f8 739 // fMaskCellColumns = new Int_t[fNMaskCellColumns];
740 // fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
741 // fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
742 // fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
743 // fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
744 // fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
a5fb4114 745
765d44e7 746}
747
765d44e7 748
cb5780f4 749//_____________________________________________________
750void AliCalorimeterUtils::InitPHOSBadChannelStatusMap()
751{
765d44e7 752 //Init PHOS bad channels map
753 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
754 //In order to avoid rewriting the same histograms
755 Bool_t oldStatus = TH1::AddDirectoryStatus();
756 TH1::AddDirectory(kFALSE);
7cc7d3f8 757
78219bac 758 fPHOSBadChannelMap = new TObjArray(5);
c5693f62 759 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 760
765d44e7 761 fPHOSBadChannelMap->SetOwner(kTRUE);
762 fPHOSBadChannelMap->Compress();
763
764 //In order to avoid rewriting the same histograms
765 TH1::AddDirectory(oldStatus);
766}
767
cb5780f4 768//______________________________________________________
769void AliCalorimeterUtils::InitPHOSRecalibrationFactors()
770{
09e819c9 771 //Init EMCAL recalibration factors
772 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
773 //In order to avoid rewriting the same histograms
774 Bool_t oldStatus = TH1::AddDirectoryStatus();
775 TH1::AddDirectory(kFALSE);
7cc7d3f8 776
78219bac 777 fPHOSRecalibrationFactors = new TObjArray(5);
c5693f62 778 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 779 //Init the histograms with 1
780 for (Int_t m = 0; m < 5; m++) {
781 for (Int_t i = 0; i < 56; i++) {
782 for (Int_t j = 0; j < 64; j++) {
c5693f62 783 SetPHOSChannelRecalibrationFactor(m,j,i,1.);
09e819c9 784 }
785 }
786 }
787 fPHOSRecalibrationFactors->SetOwner(kTRUE);
788 fPHOSRecalibrationFactors->Compress();
789
790 //In order to avoid rewriting the same histograms
791 TH1::AddDirectory(oldStatus);
792}
793
794
cb5780f4 795//___________________________________________
765d44e7 796void AliCalorimeterUtils::InitEMCALGeometry()
797{
798 //Initialize EMCAL geometry if it did not exist previously
799 if (!fEMCALGeo){
a38a48f2 800 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
765d44e7 801 if(fDebug > 0){
802 printf("AliCalorimeterUtils::InitEMCALGeometry()");
803 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
804 printf("\n");
805 }
806 }
807}
808
cb5780f4 809//__________________________________________
765d44e7 810void AliCalorimeterUtils::InitPHOSGeometry()
811{
812 //Initialize PHOS geometry if it did not exist previously
813 if (!fPHOSGeo){
814 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
815 if(fDebug > 0){
816 printf("AliCalorimeterUtils::InitPHOSGeometry()");
817 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
818 printf("\n");
819 }
820 }
821}
822
cb5780f4 823//_________________________________________________________
765d44e7 824void AliCalorimeterUtils::Print(const Option_t * opt) const
825{
7cc7d3f8 826
765d44e7 827 //Print some relevant parameters set for the analysis
828 if(! opt)
829 return;
7cc7d3f8 830
765d44e7 831 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
832 printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
833 printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
7cc7d3f8 834 fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
247abff4 835 if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
09e819c9 836 printf("Recalibrate Clusters? %d\n",fRecalibration);
9584c261 837 printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
838 printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
f2ccb5b8 839 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
7cc7d3f8 840
7db7dcb6 841 printf("Loc. Max. E > %2.2f\n", fLocMaxCutE);
842 printf("Loc. Max. E Diff > %2.2f\n", fLocMaxCutEDiff);
843
765d44e7 844 printf(" \n") ;
845}
846
cb5780f4 847//__________________________________________________________________
848Bool_t AliCalorimeterUtils::MaskFrameCluster(const Int_t iSM,
849 const Int_t ieta) const
850{
a5fb4114 851 //Check if cell is in one of the regions where we have significant amount
852 //of material in front. Only EMCAL
853
854 Int_t icol = ieta;
855 if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
856
857 if (fNMaskCellColumns && fMaskCellColumns) {
858 for (Int_t imask = 0; imask < fNMaskCellColumns; imask++) {
859 if(icol==fMaskCellColumns[imask]) return kTRUE;
860 }
861 }
862
863 return kFALSE;
864
865}
866
9369a2b1 867//__________________________________________________________________________________________
7db7dcb6 868void AliCalorimeterUtils::RecalibrateCellAmplitude(Float_t & amp,
9369a2b1 869 const TString calo, const Int_t id) const
7db7dcb6 870{
871 //Recaculate cell energy if recalibration factor
872
873 Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
874 Int_t nModule = GetModuleNumberCellIndexes(id,calo, icol, irow, iRCU);
875
876 if (IsRecalibrationOn())
877 {
878 if(calo == "PHOS")
879 {
880 amp *= GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
881 }
882 else
883 {
884 amp *= GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
885 }
886 }
887}
888
9369a2b1 889//_________________________________________________________________________________
7db7dcb6 890void AliCalorimeterUtils::RecalibrateCellTime(Double_t & time,
891 const TString calo,
9369a2b1 892 const Int_t id, const Int_t bc) const
7db7dcb6 893{
894 // Recalculate time if time recalibration available for EMCAL
895 // not ready for PHOS
896
897 if(calo == "EMCAL" && GetEMCALRecoUtils()->IsTimeRecalibrationOn())
898 {
899 GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
900 }
901
902}
903
a5fb4114 904
cb5780f4 905//__________________________________________________________________________
906Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster,
907 AliVCaloCells * cells)
908{
09e819c9 909 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
7cc7d3f8 910
898c9d44 911 //Initialize some used variables
7db7dcb6 912 Float_t frac = 0., energy = 0.;
898c9d44 913
7db7dcb6 914 if(cells)
915 {
7cc7d3f8 916 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
7db7dcb6 917
7cc7d3f8 918 UShort_t * index = cluster->GetCellsAbsId() ;
919 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
7db7dcb6 920
7cc7d3f8 921 Int_t ncells = cluster->GetNCells();
7db7dcb6 922
7cc7d3f8 923 TString calo = "EMCAL";
7db7dcb6 924 if(cluster->IsPHOS())
925 calo = "PHOS";
7cc7d3f8 926
927 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
928 for(Int_t icell = 0; icell < ncells; icell++){
7db7dcb6 929
930 Int_t absId = index[icell];
931
7cc7d3f8 932 frac = fraction[icell];
933 if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
7db7dcb6 934
935 Float_t amp = cells->GetCellAmplitude(absId);
936 RecalibrateCellAmplitude(amp,calo, absId);
937
7cc7d3f8 938 if(fDebug>2)
7db7dcb6 939 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, cell fraction %f, cell energy %f\n",
940 calo.Data(),frac,cells->GetCellAmplitude(absId));
7cc7d3f8 941
7db7dcb6 942 energy += amp*frac;
7cc7d3f8 943 }
944
945 if(fDebug>1)
946 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
898c9d44 947
948 }// cells available
7db7dcb6 949 else
950 {
898c9d44 951 Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
952 }
953
09e819c9 954 return energy;
09e819c9 955}
956
cb5780f4 957//________________________________________________________________________________
765d44e7 958void AliCalorimeterUtils::SetGeometryTransformationMatrices(AliVEvent* inputEvent)
959{
960 //Set the calorimeters transformation matrices
7cc7d3f8 961
765d44e7 962 //Get the EMCAL transformation geometry matrices from ESD
3b13c34c 963 if(!fEMCALGeoMatrixSet && fEMCALGeo){
964 if(fLoadEMCALMatrices){
c5693f62 965 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined EMCAL geometry matrices\n");
3b13c34c 966 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
967 if(fEMCALMatrix[mod]){
968 if(fDebug > 1)
969 fEMCALMatrix[mod]->Print();
970 fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;
971 }
972 }//SM loop
973 fEMCALGeoMatrixSet = kTRUE;//At least one, so good
974
975 }//Load matrices
976 else if (!gGeoManager) {
7cc7d3f8 977
3b13c34c 978 if(fDebug > 1)
979 printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load EMCAL misalignment matrices. \n");
980 if(!strcmp(inputEvent->GetName(),"AliESDEvent")) {
981 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
7cc7d3f8 982 //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
4b892846 983 if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod)) {
984 fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
3b13c34c 985 }
986 }// loop over super modules
c5693f62 987
3b13c34c 988 fEMCALGeoMatrixSet = kTRUE;//At least one, so good
989
990 }//ESD as input
991 else {
992 if(fDebug > 1)
993 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
994 }//AOD as input
995 }//Get matrix from data
04131edb 996 else if(gGeoManager){
997 fEMCALGeoMatrixSet = kTRUE;
998 }
765d44e7 999 }//EMCAL geo && no geoManager
c5693f62 1000
765d44e7 1001 //Get the PHOS transformation geometry matrices from ESD
3b13c34c 1002 if(!fPHOSGeoMatrixSet && fPHOSGeo){
c5693f62 1003
3b13c34c 1004 if(fLoadPHOSMatrices){
c5693f62 1005 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined PHOS geometry matrices\n");
1006 for(Int_t mod = 0 ; mod < 5 ; mod++){
3b13c34c 1007 if(fPHOSMatrix[mod]){
c5693f62 1008 if(fDebug > 1 )
3b13c34c 1009 fPHOSMatrix[mod]->Print();
dbf54f1e 1010 fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;
3b13c34c 1011 }
1012 }//SM loop
1013 fPHOSGeoMatrixSet = kTRUE;//At least one, so good
c5693f62 1014
3b13c34c 1015 }//Load matrices
1016 else if (!gGeoManager) {
1017 if(fDebug > 1)
1018 printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
765d44e7 1019 if(!strcmp(inputEvent->GetName(),"AliESDEvent")) {
dbf54f1e 1020 for(Int_t mod = 0; mod < 5; mod++){
4b892846 1021 if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod)) {
765d44e7 1022 //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
4b892846 1023 fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
765d44e7 1024 }
3b13c34c 1025 }// loop over modules
1026 fPHOSGeoMatrixSet = kTRUE; //At least one so good
765d44e7 1027 }//ESD as input
1028 else {
1029 if(fDebug > 1)
1030 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
3b13c34c 1031 }//AOD as input
1032 }// get matrix from data
04131edb 1033 else if(gGeoManager){
1034 fPHOSGeoMatrixSet = kTRUE;
1035 }
765d44e7 1036 }//PHOS geo and geoManager was not set
c5693f62 1037
765d44e7 1038}
1039
cb5780f4 1040//__________________________________________________________________________________________
1041void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
1042{
9584c261 1043
1044 //Recalculate EMCAL cluster position
1045
19db8f8c 1046 fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
7cc7d3f8 1047
9584c261 1048}
cb5780f4 1049
3c1d9afb 1050//___________________________________________________________________________
1051void AliCalorimeterUtils::SplitEnergy(const Int_t absId1, const Int_t absId2,
1052 AliVCluster* cluster,
1053 AliVCaloCells* cells,
1054 //Float_t & e1, Float_t & e2,
1055 AliAODCaloCluster* cluster1,
1056 AliAODCaloCluster* cluster2,
1057 const Int_t nMax,
1058 const Int_t eventNumber)
1059{
1060
1061 // Split energy of cluster between the 2 local maxima, sum energy on 3x3, and if the 2
1062 // maxima are too close and have common cells, split the energy between the 2
1063
1064 TH2F* hClusterMap = 0 ;
1065 TH2F* hClusterLocMax = 0 ;
1066 TH2F* hCluster1 = 0 ;
1067 TH2F* hCluster2 = 0 ;
1068
1069 if(fPlotCluster)
1070 {
1071 hClusterMap = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
1072 hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
1073 hCluster1 = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
1074 hCluster2 = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
1075 hClusterMap ->SetXTitle("column");
1076 hClusterMap ->SetYTitle("row");
1077 hClusterLocMax ->SetXTitle("column");
1078 hClusterLocMax ->SetYTitle("row");
1079 hCluster1 ->SetXTitle("column");
1080 hCluster1 ->SetYTitle("row");
1081 hCluster2 ->SetXTitle("column");
1082 hCluster2 ->SetYTitle("row");
1083 }
1084
1085 TString calorimeter = "EMCAL";
1086 if(cluster->IsPHOS())
1087 {
1088 calorimeter="PHOS";
1089 printf("AliCalorimeterUtils::SplitEnerg() Not supported for PHOS yet \n");
1090 return;
1091 }
1092
1093 const Int_t ncells = cluster->GetNCells();
1094 Int_t absIdList[ncells];
1095
1096 Float_t e1 = 0, e2 = 0 ;
1097 Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;
1098 Float_t eCluster = 0;
1099 Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1;
1100 for(Int_t iDigit = 0; iDigit < ncells; iDigit++ )
1101 {
1102 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
1103
1104 Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
1105 RecalibrateCellAmplitude(ec,calorimeter, absIdList[iDigit]);
1106 eCluster+=ec;
1107
1108 if(fPlotCluster)
1109 {
1110 //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
1111 sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1112 if(icol > maxCol) maxCol = icol;
1113 if(icol < minCol) minCol = icol;
1114 if(irow > maxRow) maxRow = irow;
1115 if(irow < minRow) minRow = irow;
1116 hClusterMap->Fill(icol,irow,ec);
1117 }
1118
1119 }
1120
1121 // Init counters and variables
1122 Int_t ncells1 = 1 ;
1123 UShort_t absIdList1[9] ;
1124 Double_t fracList1 [9] ;
1125 absIdList1[0] = absId1 ;
1126 fracList1 [0] = 1. ;
1127
1128 Float_t ecell1 = cells->GetCellAmplitude(absId1);
1129 RecalibrateCellAmplitude(ecell1, calorimeter, absId1);
1130 e1 = ecell1;
1131
1132 Int_t ncells2 = 1 ;
1133 UShort_t absIdList2[9] ;
1134 Double_t fracList2 [9] ;
1135 absIdList2[0] = absId2 ;
1136 fracList2 [0] = 1. ;
1137
1138 Float_t ecell2 = cells->GetCellAmplitude(absId2);
1139 RecalibrateCellAmplitude(ecell2, calorimeter, absId2);
1140 e2 = ecell2;
1141
1142 if(fPlotCluster)
1143 {
1144 Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
1145 sm = GetModuleNumberCellIndexes(absId1, calorimeter, icol1, irow1, iRCU) ;
1146 hClusterLocMax->Fill(icol1,irow1,ecell1);
1147 sm = GetModuleNumberCellIndexes(absId2, calorimeter, icol2, irow2, iRCU) ;
1148 hClusterLocMax->Fill(icol2,irow2,ecell2);
1149 }
1150
1151 // Very rough way to share the cluster energy
1152 Float_t eRemain = (eCluster-ecell1-ecell2)/2;
1153 Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
1154 Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
1155
1156 for(Int_t iDigit = 0; iDigit < ncells; iDigit++)
1157 {
1158 Int_t absId = absIdList[iDigit];
1159
1160 if(absId==absId1 || absId==absId2 || absId < 0) continue;
1161
1162 Float_t ecell = cells->GetCellAmplitude(absId);
1163 RecalibrateCellAmplitude(ecell, calorimeter, absId);
1164
1165 if(AreNeighbours(calorimeter, absId1,absId ))
1166 {
1167 absIdList1[ncells1]= absId;
1168
1169 if(AreNeighbours(calorimeter, absId2,absId ))
1170 {
1171 fracList1[ncells1] = shareFraction1;
1172 e1 += ecell*shareFraction1;
1173 }
1174 else
1175 {
1176 fracList1[ncells1] = 1.;
1177 e1 += ecell;
1178 }
1179
1180 ncells1++;
1181
1182 } // neigbour to cell1
1183
1184 if(AreNeighbours(calorimeter, absId2,absId ))
1185 {
1186 absIdList2[ncells2]= absId;
1187
1188 if(AreNeighbours(calorimeter, absId1,absId ))
1189 {
1190 fracList2[ncells2] = shareFraction2;
1191 e2 += ecell*shareFraction2;
1192 }
1193 else
1194 {
1195 fracList2[ncells2] = 1.;
1196 e2 += ecell;
1197 }
1198
1199 ncells2++;
1200
1201 } // neigbour to cell2
1202
1203 }
1204
1205 if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::SplitEnergy() - n Local Max %d, Cluster energy = %f, Ecell1 = %f, Ecell2 = %f, Enew1 = %f, Enew2 = %f, Remain %f, \n ncells %d, ncells1 %d, ncells2 %d, f1 %f, f2 %f, sum f12 = %f \n",
1206 nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2);
1207
1208 cluster1->SetE(e1);
1209 cluster2->SetE(e2);
1210
1211 cluster1->SetNCells(ncells1);
1212 cluster2->SetNCells(ncells2);
1213
1214 cluster1->SetCellsAbsId(absIdList1);
1215 cluster2->SetCellsAbsId(absIdList2);
1216
1217 cluster1->SetCellsAmplitudeFraction(fracList1);
1218 cluster2->SetCellsAmplitudeFraction(fracList2);
1219
1220 if(calorimeter=="EMCAL")
1221 {
1222 GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
1223 GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
1224 }
1225
1226 if(fPlotCluster)
1227 {
1228 //printf("Cells of cluster1: ");
1229 for(Int_t iDigit = 0; iDigit < ncells1; iDigit++ )
1230 {
1231 //printf(" %d ",absIdList1[iDigit]);
1232
1233 sm = GetModuleNumberCellIndexes(absIdList1[iDigit], calorimeter, icol, irow, iRCU) ;
1234
1235 if( AreNeighbours(calorimeter, absId2,absIdList1[iDigit]) )
1236 hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit])*shareFraction1);
1237 else
1238 hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit]));
1239 }
1240
1241 //printf(" \n ");
1242 //printf("Cells of cluster2: ");
1243
1244 for(Int_t iDigit = 0; iDigit < ncells2; iDigit++ )
1245 {
1246 //printf(" %d ",absIdList2[iDigit]);
1247
1248 sm = GetModuleNumberCellIndexes(absIdList2[iDigit], calorimeter, icol, irow, iRCU) ;
1249 if( AreNeighbours(calorimeter, absId1,absIdList2[iDigit]) )
1250 hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit])*shareFraction2);
1251 else
1252 hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit]));
1253
1254 }
1255 //printf(" \n ");
1256
1257 gStyle->SetPadRightMargin(0.1);
1258 gStyle->SetPadLeftMargin(0.1);
1259 gStyle->SetOptStat(0);
1260 gStyle->SetOptFit(000000);
1261
1262 if(maxCol-minCol > maxRow-minRow)
1263 {
1264 maxRow+= maxCol-minCol;
1265 }
1266 else
1267 {
1268 maxCol+= maxRow-minRow;
1269 }
1270
1271 TCanvas * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
1272 c->Divide(2,2);
1273 c->cd(1);
1274 gPad->SetGridy();
1275 gPad->SetGridx();
1276 hClusterMap ->SetAxisRange(minCol, maxCol,"X");
1277 hClusterMap ->SetAxisRange(minRow, maxRow,"Y");
1278 hClusterMap ->Draw("colz");
1279 c->cd(2);
1280 gPad->SetGridy();
1281 gPad->SetGridx();
1282 hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
1283 hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
1284 hClusterLocMax ->Draw("colz");
1285 c->cd(3);
1286 gPad->SetGridy();
1287 gPad->SetGridx();
1288 hCluster1 ->SetAxisRange(minCol, maxCol,"X");
1289 hCluster1 ->SetAxisRange(minRow, maxRow,"Y");
1290 hCluster1 ->Draw("colz");
1291 c->cd(4);
1292 gPad->SetGridy();
1293 gPad->SetGridx();
1294 hCluster2 ->SetAxisRange(minCol, maxCol,"X");
1295 hCluster2 ->SetAxisRange(minRow, maxRow,"Y");
1296 hCluster2 ->Draw("colz");
1297
1298 if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",
1299 eventNumber,cluster->E(),nMax,ncells1,ncells2));
1300
1301 delete c;
1302 delete hClusterMap;
1303 delete hClusterLocMax;
1304 delete hCluster1;
1305 delete hCluster2;
1306 }
1307}
1308
cb5780f4 1309//________________________________________________________________________________
1310void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event,
1311 TObjArray* clusterArray)
1312{
1313 //Recalculate track matching
1314
3c1d9afb 1315 if (fRecalculateMatching)
1316 {
cb5780f4 1317 fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo) ;
1318 //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1319 //if(esdevent){
1320 // fEMCALRecoUtils->SetClusterMatchedToTrack(esdevent) ;
1321 // fEMCALRecoUtils->SetTracksMatchedToCluster(esdevent) ;
1322 //}
1323 }
1324}