coverity fix: histograms initialized in the member initialization list
[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>
55d66f31 31#include <TFile.h>
765d44e7 32
c5693f62 33// --- ANALYSIS system ---
765d44e7 34#include "AliCalorimeterUtils.h"
4b892846 35#include "AliESDEvent.h"
46a3cde6 36#include "AliMCEvent.h"
37#include "AliStack.h"
765d44e7 38#include "AliAODPWG4Particle.h"
c8fe2783 39#include "AliVCluster.h"
40#include "AliVCaloCells.h"
41#include "AliMixedEvent.h"
3c1d9afb 42#include "AliAODCaloCluster.h"
55d66f31 43#include "AliOADBContainer.h"
44#include "AliAnalysisManager.h"
765d44e7 45
c5693f62 46// --- Detector ---
47#include "AliEMCALGeometry.h"
48#include "AliPHOSGeoUtils.h"
49
765d44e7 50ClassImp(AliCalorimeterUtils)
51
52
cb5780f4 53//____________________________________________
765d44e7 54 AliCalorimeterUtils::AliCalorimeterUtils() :
55 TObject(), fDebug(0),
12e91266 56 fEMCALGeoName(""),
57 fPHOSGeoName (""),
9e8998b1 58 fEMCALGeo(0x0), fPHOSGeo(0x0),
59 fEMCALGeoMatrixSet(kFALSE), fPHOSGeoMatrixSet(kFALSE),
60 fLoadEMCALMatrices(kFALSE), fLoadPHOSMatrices(kFALSE),
61 fRemoveBadChannels(kFALSE), fPHOSBadChannelMap(0x0),
a5fb4114 62 fNCellsFromPHOSBorder(0),
9e8998b1 63 fNMaskCellColumns(0), fMaskCellColumns(0x0),
64 fRecalibration(kFALSE), fPHOSRecalibrationFactors(),
247abff4 65 fEMCALRecoUtils(new AliEMCALRecoUtils),
9e8998b1 66 fRecalculatePosition(kFALSE), fCorrectELinearity(kFALSE),
67 fRecalculateMatching(kFALSE),
68 fCutR(20), fCutZ(20),
7db7dcb6 69 fCutEta(20), fCutPhi(20),
3c1d9afb 70 fLocMaxCutE(0), fLocMaxCutEDiff(0),
55d66f31 71 fPlotCluster(0), fOADBSet(kFALSE),
72 fOADBForEMCAL(kFALSE), fOADBForPHOS(kFALSE),
73 fOADBFilePathEMCAL(""), fOADBFilePathPHOS("")
765d44e7 74{
75 //Ctor
76
77 //Initialize parameters
78 InitParameters();
90e32961 79 for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
80 for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix [i] = 0 ;
7cc7d3f8 81
765d44e7 82}
765d44e7 83
cb5780f4 84//_________________________________________
85AliCalorimeterUtils::~AliCalorimeterUtils()
86{
765d44e7 87 //Dtor
88
4df35693 89 //if(fPHOSGeo) delete fPHOSGeo ;
765d44e7 90 if(fEMCALGeo) delete fEMCALGeo ;
7cc7d3f8 91
765d44e7 92 if(fPHOSBadChannelMap) {
93 fPHOSBadChannelMap->Clear();
94 delete fPHOSBadChannelMap;
95 }
96
09e819c9 97 if(fPHOSRecalibrationFactors) {
78219bac 98 fPHOSRecalibrationFactors->Clear();
99 delete fPHOSRecalibrationFactors;
09e819c9 100 }
101
a5fb4114 102 if(fEMCALRecoUtils) delete fEMCALRecoUtils ;
103 if(fNMaskCellColumns) delete [] fMaskCellColumns;
9584c261 104
765d44e7 105}
106
55d66f31 107//____________________________________________________
108void AliCalorimeterUtils::AccessOADB(AliVEvent* event)
109{
110 // Set the AODB calibration, bad channels etc. parameters at least once
111 // alignment matrices from OADB done in SetGeometryMatrices
112
113 //Set it only once
114 if(fOADBSet) return ;
115
116 Int_t runnumber = event->GetRunNumber() ;
117 TString pass = GetPass();
118
119 // EMCAL
120 if(fOADBForEMCAL)
121 {
122 printf("AliCalorimeterUtils::SetOADBParameters() - Get AODB parameters from EMCAL in %s for run %d, and <%s> \n",fOADBFilePathEMCAL.Data(),runnumber,pass.Data());
123
124 Int_t nSM = fEMCALGeo->GetNumberOfSuperModules();
125
126 // Bad map
127 if(fRemoveBadChannels)
128 {
129 AliOADBContainer *contBC=new AliOADBContainer("");
130 contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePathEMCAL.Data()),"AliEMCALBadChannels");
131
132 TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber);
133
134 if(arrayBC)
135 {
136
137 TObjArray * arrayBCpass = 0x0;
138 if(pass!="")arrayBCpass = (TObjArray*)arrayBC->FindObject(pass);
139 // There are no passes for simulation, in order to get the bad map put pass1
140 else arrayBCpass = (TObjArray*)arrayBC->FindObject("pass1");
141
142 if(arrayBCpass)
143 {
144 SwitchOnDistToBadChannelRecalculation();
145 printf("AliCalorimeterUtils::SetOADBParameters() - Remove EMCAL bad cells \n");
146
147 for (Int_t i=0; i<nSM; ++i)
148 {
149 TH2I *hbm = GetEMCALChannelStatusMap(i);
150
151 if (hbm)
152 delete hbm;
153
154 hbm=(TH2I*)arrayBCpass->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
155
156 if (!hbm)
157 {
158 AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
159 continue;
160 }
161
162 hbm->SetDirectory(0);
163 SetEMCALChannelStatusMap(i,hbm);
164
165 } // loop
166 } else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT remove EMCAL bad channels 1\n"); // pass array
167 } else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT remove EMCAL bad channels 2\n"); // run array
168 } // Remove bad
169
170 // Energy Recalibration
171 if(fRecalibration)
172 {
173 AliOADBContainer *contRF=new AliOADBContainer("");
174
175 contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALRecalib");
176
177 TObjArray *recal=(TObjArray*)contRF->GetObject(runnumber);
178
179 if(recal)
180 {
181 TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
182
183 if(recalpass)
184 {
185 TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
186
187 if(recalib)
188 {
189 printf("AliCalorimeterUtils::SetOADBParameters() - Recalibrate EMCAL \n");
190 for (Int_t i=0; i<nSM; ++i)
191 {
192 TH2F *h = GetEMCALChannelRecalibrationFactors(i);
193
194 if (h)
195 delete h;
196
197 h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
198
199 if (!h)
200 {
201 AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
202 continue;
203 }
204
205 h->SetDirectory(0);
206
207 SetEMCALChannelRecalibrationFactors(i,h);
208 } // SM loop
209 }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL 1\n"); // array ok
210 }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL 2\n"); // array pass ok
211 }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL 3\n"); // run number array ok
212
213 // once set, apply run dependent corrections if requested
214 fEMCALRecoUtils->SetRunDependentCorrections(runnumber);
215
216 } // Recalibration on
217
218
219 // Time Recalibration
220 if(fEMCALRecoUtils->IsTimeRecalibrationOn())
221 {
222 AliOADBContainer *contTRF=new AliOADBContainer("");
223
f46af216 224 contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALTimeCalib");
55d66f31 225
226 TObjArray *trecal=(TObjArray*)contTRF->GetObject(runnumber);
227
228 if(trecal)
229 {
230 TObjArray *trecalpass=(TObjArray*)trecal->FindObject(pass);
231
232 if(trecalpass)
233 {
f46af216 234 printf("AliCalorimeterUtils::SetOADBParameters() - Time Recalibrate EMCAL \n");
235 for (Int_t ibc = 0; ibc < 4; ++ibc)
55d66f31 236 {
f46af216 237 TH1F *h = GetEMCALChannelTimeRecalibrationFactors(ibc);
238
239 if (h)
240 delete h;
241
242 h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
243
244 if (!h)
55d66f31 245 {
f46af216 246 AliError(Form("Could not load hAllTimeAvBC%d",ibc));
247 continue;
248 }
249
250 h->SetDirectory(0);
251
252 SetEMCALChannelTimeRecalibrationFactors(ibc,h);
253 } // bunch crossing loop
254 }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate time EMCAL 1\n"); // array pass ok
255 }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate time EMCAL 2\n"); // run number array ok
55d66f31 256
257 } // Recalibration on
258
259 }// EMCAL
260
261 // PHOS
262 if(fOADBForPHOS)
263 {
264 printf("AliCalorimeterUtils::SetOADBParameters() - Get AODB parameters from PHOS in %s for run %d, and <%s> \n",fOADBFilePathPHOS.Data(),runnumber,pass.Data());
265
266 // Bad map
267 if(fRemoveBadChannels)
268 {
269 AliOADBContainer badmapContainer(Form("phosBadMap"));
270 TString fileName="$ALICE_ROOT/OADB/PHOS/PHOSBadMaps.root";
271 badmapContainer.InitFromFile(Form("%s/PHOSBadMaps.root",fOADBFilePathPHOS.Data()),"phosBadMap");
272
273 //Use a fixed run number from year 2010, this year not available yet.
274 TObjArray *maps = (TObjArray*)badmapContainer.GetObject(139000,"phosBadMap");
275 if(!maps)
276 {
277 printf("AliCalorimeterUtils::SetOADBParameters() - Can not read PHOS bad map for run %d.\n",runnumber) ;
278 }
279 else
280 {
281 printf("AliCalorimeterUtils::SetOADBParameters() - Setting PHOS bad map with name %s \n",maps->GetName()) ;
282 for(Int_t mod=1; mod<5;mod++)
283 {
284 TH2I *hbmPH = GetPHOSChannelStatusMap(mod);
285
286 if(hbmPH)
287 delete hbmPH ;
288
289 hbmPH = (TH2I*)maps->At(mod);
290
291 if(hbmPH) hbmPH->SetDirectory(0);
292
293 SetPHOSChannelStatusMap(mod-1,hbmPH);
294
295 } // modules loop
296 } // maps exist
297 } // Remove bad channels
298 } // PHOS
299
300 // Parameters already set once, so do not it again
301 fOADBSet = kTRUE;
302
303}
304
305//_____________________________________________________________
306void AliCalorimeterUtils::AccessGeometry(AliVEvent* inputEvent)
307{
308 //Set the calorimeters transformation matrices and init geometry
309
310 // First init the geometry, a priory not done before
311 Int_t runnumber = inputEvent->GetRunNumber() ;
312 InitPHOSGeometry (runnumber);
313 InitEMCALGeometry(runnumber);
314
315 //Get the EMCAL transformation geometry matrices from ESD
316 if(!fEMCALGeoMatrixSet && fEMCALGeo)
317 {
318 if(fLoadEMCALMatrices)
319 {
320 printf("AliCalorimeterUtils::AccessGeometry() - Load user defined EMCAL geometry matrices\n");
321
322 // OADB if available
323 AliOADBContainer emcGeoMat("AliEMCALgeo");
324 emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePathEMCAL.Data()),"AliEMCALgeo");
325 TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
326
327 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
328 {
329 if (!fEMCALMatrix[mod]) // Get it from OADB
330 {
331 if(fDebug > 1 )
332 printf("AliCalorimeterUtils::AccessGeometry() - EMCAL matrices SM %d, %p\n",
333 mod,((TGeoHMatrix*) matEMCAL->At(mod)));
334 //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
335
336 fEMCALMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
337 }
338
339 if(fEMCALMatrix[mod])
340 {
341 if(fDebug > 1)
342 fEMCALMatrix[mod]->Print();
343
344 fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;
345 }
346 }//SM loop
347
348 fEMCALGeoMatrixSet = kTRUE;//At least one, so good
349
350 }//Load matrices
351 else if (!gGeoManager)
352 {
353 if(fDebug > 1)
354 printf(" AliCalorimeterUtils::AccessGeometry() - Load EMCAL misalignment matrices. \n");
355 if(!strcmp(inputEvent->GetName(),"AliESDEvent"))
356 {
357 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
358 {
359 //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
360 if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod))
361 {
362 fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
363 }
364 }// loop over super modules
365
366 fEMCALGeoMatrixSet = kTRUE;//At least one, so good
367
368 }//ESD as input
369 else
370 {
371 if(fDebug > 1)
372 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
373 }//AOD as input
374 }//Get matrix from data
375 else if(gGeoManager)
376 {
377 fEMCALGeoMatrixSet = kTRUE;
378 }
379 }//EMCAL geo && no geoManager
380
381 //Get the PHOS transformation geometry matrices from ESD
382 if(!fPHOSGeoMatrixSet && fPHOSGeo)
383 {
384 if(fLoadPHOSMatrices)
385 {
386 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined PHOS geometry matrices\n");
387
388 // OADB if available
389 AliOADBContainer geomContainer("phosGeo");
390 geomContainer.InitFromFile(Form("%s/PHOSGeometry.root",fOADBFilePathPHOS.Data()),"PHOSRotationMatrixes");
391 TObjArray *matPHOS = (TObjArray*)geomContainer.GetObject(139000,"PHOSRotationMatrixes");
392
393 for(Int_t mod = 0 ; mod < 5 ; mod++)
394 {
395 if (!fPHOSMatrix[mod]) // Get it from OADB
396 {
397 if(fDebug > 1 )
398 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - PHOS matrices module %d, %p\n",
399 mod,((TGeoHMatrix*) matPHOS->At(mod)));
400 //((TGeoHMatrix*) matPHOS->At(mod))->Print();
401
402 fPHOSMatrix[mod] = (TGeoHMatrix*) matPHOS->At(mod) ;
403 }
404
405 // Set it, if it exists
406 if(fPHOSMatrix[mod])
407 {
408 if(fDebug > 1 )
409 fPHOSMatrix[mod]->Print();
410
411 fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;
412 }
413 }// SM loop
414
415 fPHOSGeoMatrixSet = kTRUE;//At least one, so good
416
417 }//Load matrices
418 else if (!gGeoManager)
419 {
420 if(fDebug > 1)
421 printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
422 if(!strcmp(inputEvent->GetName(),"AliESDEvent"))
423 {
424 for(Int_t mod = 0; mod < 5; mod++)
425 {
426 if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod))
427 {
428 //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
429 fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
430 }
431 }// loop over modules
432 fPHOSGeoMatrixSet = kTRUE; //At least one so good
433 }//ESD as input
434 else
435 {
436 if(fDebug > 1)
437 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
438 }//AOD as input
439 }// get matrix from data
440 else if(gGeoManager)
441 {
442 fPHOSGeoMatrixSet = kTRUE;
443 }
444 }//PHOS geo and geoManager was not set
445
446}
447
7db7dcb6 448//______________________________________________________________________________________
449Bool_t AliCalorimeterUtils::AreNeighbours(const TString calo,
450 const Int_t absId1, const Int_t absId2 ) const
451{
452 // Tells if (true) or not (false) two cells are neighbours
453 // A neighbour is defined as being two cells which share a side or corner
454
455 Bool_t areNeighbours = kFALSE ;
456
457 Int_t iRCU1 = -1, irow1 = -1, icol1 = -1;
458 Int_t iRCU2 = -1, irow2 = -1, icol2 = -1;
459
460 Int_t rowdiff = 0, coldiff = 0;
461
9369a2b1 462 Int_t nSupMod1 = GetModuleNumberCellIndexes(absId1, calo, icol1, irow1, iRCU1);
463 Int_t nSupMod2 = GetModuleNumberCellIndexes(absId2, calo, icol2, irow2, iRCU2);
7db7dcb6 464
465 if(calo=="EMCAL" && nSupMod1!=nSupMod2)
466 {
467 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
468 // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
469 if(nSupMod1%2) icol1+=AliEMCALGeoParams::fgkEMCALCols;
9369a2b1 470 else icol2+=AliEMCALGeoParams::fgkEMCALCols;
7db7dcb6 471 }
472
473 rowdiff = TMath::Abs( irow1 - irow2 ) ;
474 coldiff = TMath::Abs( icol1 - icol2 ) ;
475
476 if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
477 areNeighbours = kTRUE ;
478
479 return areNeighbours;
480
481}
482
483
cb5780f4 484//_____________________________________________________________________________________
485Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster,
486 AliVCaloCells* cells,
487 AliVEvent * event, Int_t iev) const
488{
7cc7d3f8 489
765d44e7 490 // Given the list of AbsId of the cluster, get the maximum cell and
491 // check if there are fNCellsFromBorder from the calorimeter border
492
7cc7d3f8 493 //If the distance to the border is 0 or negative just exit accept all clusters
247abff4 494 if(cells->GetType()==AliVCaloCells::kEMCALCell && fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder() <= 0 ) return kTRUE;
c8fe2783 495 if(cells->GetType()==AliVCaloCells::kPHOSCell && fNCellsFromPHOSBorder <= 0 ) return kTRUE;
7cc7d3f8 496
c8fe2783 497 Int_t absIdMax = -1;
765d44e7 498 Float_t ampMax = -1;
c8fe2783 499
500 AliMixedEvent * mixEvent = dynamic_cast<AliMixedEvent*> (event);
501 Int_t nMixedEvents = 0 ;
502 Int_t * cellsCumul = NULL ;
503 Int_t numberOfCells = 0 ;
504 if (mixEvent){
505 nMixedEvents = mixEvent->GetNumberOfEvents() ;
506 if (cells->GetType()==AliVCaloCells::kEMCALCell) {
507 cellsCumul = mixEvent->GetEMCALCellsCumul() ;
508 numberOfCells = mixEvent->GetNumberOfEMCALCells() ;
509 }
7cc7d3f8 510
c8fe2783 511 else if (cells->GetType()==AliVCaloCells::kPHOSCell) {
512 cellsCumul = mixEvent->GetPHOSCellsCumul() ;
513 numberOfCells = mixEvent->GetNumberOfPHOSCells() ;
514 }
898c9d44 515
516 if(cellsCumul){
517
518 Int_t startCell = cellsCumul[iev] ;
519 Int_t endCell = (iev+1 < nMixedEvents)?cellsCumul[iev+1]:numberOfCells;
c8fe2783 520 //Find cells with maximum amplitude
898c9d44 521 for(Int_t i = 0; i < cluster->GetNCells() ; i++){
522 Int_t absId = cluster->GetCellAbsId(i) ;
523 for (Int_t j = startCell; j < endCell ; j++) {
524 Short_t cellNumber;
525 Double_t amp ;
526 Double_t time;
527 cells->GetCell(j, cellNumber, amp, time) ;
528 if (absId == cellNumber) {
529 if(amp > ampMax){
530 ampMax = amp;
531 absIdMax = absId;
532 }
533 }
c8fe2783 534 }
898c9d44 535 }//loop on cluster cells
536 }// cells cumul available
537 else {
538 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - CellsCumul is NULL!!!\n");
539 abort();
c8fe2783 540 }
898c9d44 541 } else {//Normal SE Events
c8fe2783 542 for(Int_t i = 0; i < cluster->GetNCells() ; i++){
543 Int_t absId = cluster->GetCellAbsId(i) ;
544 Float_t amp = cells->GetCellAmplitude(absId);
545 if(amp > ampMax){
546 ampMax = amp;
547 absIdMax = absId;
548 }
549 }
550 }
765d44e7 551
552 if(fDebug > 1)
280e6766 553 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n",
7cc7d3f8 554 absIdMax, ampMax, cluster->E());
765d44e7 555
556 if(absIdMax==-1) return kFALSE;
557
558 //Check if the cell is close to the borders:
559 Bool_t okrow = kFALSE;
560 Bool_t okcol = kFALSE;
7cc7d3f8 561
c8fe2783 562 if(cells->GetType()==AliVCaloCells::kEMCALCell){
7cc7d3f8 563
765d44e7 564 Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1;
565 fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta);
566 fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
280e6766 567 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
568 Fatal("CheckCellFidutialRegion","Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",iSM,ieta,iphi);
569 }
7cc7d3f8 570
765d44e7 571 //Check rows/phi
247abff4 572 Int_t nborder = fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder();
5408e59e 573 if(iSM < 10)
574 {
247abff4 575 if(iphi >= nborder && iphi < 24-nborder) okrow =kTRUE;
7cc7d3f8 576 }
5408e59e 577 else
578 {
579 if(fEMCALGeoName.Contains("12SM")) // 1/3 SM
580 if(iphi >= nborder && iphi < 8-nborder) okrow =kTRUE;
581 else // 1/2 SM
582 if(iphi >= nborder && iphi <12-nborder) okrow =kTRUE;
765d44e7 583 }
584
247abff4 585 //Check columns/eta
5408e59e 586 if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0())
587 {
247abff4 588 if(ieta > nborder && ieta < 48-nborder) okcol =kTRUE;
765d44e7 589 }
5408e59e 590 else
591 {
592 if(iSM%2==0)
593 {
247abff4 594 if(ieta >= nborder) okcol = kTRUE;
765d44e7 595 }
5408e59e 596 else
597 {
247abff4 598 if(ieta < 48-nborder) okcol = kTRUE;
765d44e7 599 }
600 }//eta 0 not checked
5408e59e 601
765d44e7 602 if(fDebug > 1)
603 {
280e6766 604 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
7cc7d3f8 605 nborder, ieta, iphi, iSM);
765d44e7 606 if (okcol && okrow ) printf(" YES \n");
607 else printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
608 }
609 }//EMCAL
c8fe2783 610 else if(cells->GetType()==AliVCaloCells::kPHOSCell){
765d44e7 611 Int_t relId[4];
612 Int_t irow = -1, icol = -1;
613 fPHOSGeo->AbsToRelNumbering(absIdMax,relId);
614 irow = relId[2];
615 icol = relId[3];
616 //imod = relId[0]-1;
617 if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE;
618 if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE;
619 if(fDebug > 1)
620 {
280e6766 621 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - PHOS Cluster in %d cells fiducial volume: icol %d, irow %d, Module %d?",
7cc7d3f8 622 fNCellsFromPHOSBorder, icol, irow, relId[0]-1);
765d44e7 623 if (okcol && okrow ) printf(" YES \n");
624 else printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
625 }
626 }//PHOS
627
628 if (okcol && okrow) return kTRUE;
629 else return kFALSE;
630
631}
632
765d44e7 633//_________________________________________________________________________________________________________
cb5780f4 634Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter,UShort_t* cellList, Int_t nCells)
635{
765d44e7 636 // Check that in the cluster cells, there is no bad channel of those stored
637 // in fEMCALBadChannelMap or fPHOSBadChannelMap
638
639 if (!fRemoveBadChannels) return kFALSE;
36037088 640 //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
247abff4 641 if(calorimeter == "EMCAL" && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
78219bac 642 if(calorimeter == "PHOS" && !fPHOSBadChannelMap) return kFALSE;
7cc7d3f8 643
765d44e7 644 Int_t icol = -1;
645 Int_t irow = -1;
646 Int_t imod = -1;
647 for(Int_t iCell = 0; iCell<nCells; iCell++){
7cc7d3f8 648
765d44e7 649 //Get the column and row
650 if(calorimeter == "EMCAL"){
247abff4 651 return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
765d44e7 652 }
653 else if(calorimeter=="PHOS"){
654 Int_t relId[4];
655 fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
656 irow = relId[2];
657 icol = relId[3];
658 imod = relId[0]-1;
659 if(fPHOSBadChannelMap->GetEntries() <= imod)continue;
c5693f62 660 //printf("PHOS bad channels imod %d, icol %d, irow %d\n",imod, irow, icol);
661 if(GetPHOSChannelStatus(imod, irow, icol)) return kTRUE;
765d44e7 662 }
663 else return kFALSE;
664
665 }// cell cluster loop
7cc7d3f8 666
765d44e7 667 return kFALSE;
7cc7d3f8 668
765d44e7 669}
670
cb5780f4 671//_______________________________________________________________
672void AliCalorimeterUtils::CorrectClusterEnergy(AliVCluster *clus)
673{
9584c261 674 // Correct cluster energy non linearity
13cd2872 675
9584c261 676 clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
13cd2872 677
678}
679
c5693f62 680//________________________________________________________________________________________
681Int_t AliCalorimeterUtils::GetMaxEnergyCell(AliVCaloCells* cells, const AliVCluster* clu,
682 Float_t & clusterFraction) const
cb5780f4 683{
13cd2872 684
685 //For a given CaloCluster gets the absId of the cell
686 //with maximum energy deposit.
687
688 if( !clu || !cells ){
689 AliInfo("Cluster or cells pointer is null!");
690 return -1;
691 }
692
693 Double_t eMax =-1.;
694 Double_t eTot = 0.;
695 Double_t eCell =-1.;
696 Float_t fraction = 1.;
697 Float_t recalFactor = 1.;
698 Int_t cellAbsId =-1 , absId =-1 ;
699 Int_t iSupMod =-1 , ieta =-1 , iphi = -1, iRCU = -1;
700
701 TString calo = "EMCAL";
702 if(clu->IsPHOS()) calo = "PHOS";
703
704 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
705
706 cellAbsId = clu->GetCellAbsId(iDig);
707
708 fraction = clu->GetCellAmplitudeFraction(iDig);
709 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
710
b08dd6d5 711 iSupMod = GetModuleNumberCellIndexes(cellAbsId, calo, ieta, iphi, iRCU);
13cd2872 712
713 if(IsRecalibrationOn()) {
714 if(calo=="EMCAL") recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
c5693f62 715 else recalFactor = GetPHOSChannelRecalibrationFactor (iSupMod,iphi,ieta);
13cd2872 716 }
717
718 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
719
720 if(eCell > eMax) {
721 eMax = eCell;
722 absId = cellAbsId;
723 }
724
725 eTot+=eCell;
726
727 }// cell loop
728
729 clusterFraction = (eTot-eMax)/eTot; //Do not use cluster energy in case it was corrected for non linearity.
730
731 return absId;
732
9584c261 733}
734
d832b695 735//__________________________________________________________________________
736AliVTrack * AliCalorimeterUtils::GetMatchedTrack(const AliVCluster* cluster,
737 const AliVEvent* event,
738 const Int_t index) const
739{
740 // Get the matched track given its index, usually just the first match
741 // Since it is different for ESDs and AODs here it is a wrap method to do it
742
743 AliVTrack *track = 0;
744
745 // EMCAL case only when matching is recalculated
746 if(cluster->IsEMCAL() && IsRecalculationOfClusterTrackMatchingOn())
747 {
748 Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(cluster->GetID());
749 //printf("track index %d, cluster ID %d \n ",trackIndex,cluster->GetID());
750
751 if(trackIndex < 0 )
752 {
753 printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d, from recalculation\n", trackIndex);
754 }
755 else
756 {
757 track = dynamic_cast<AliVTrack*> (event->GetTrack(trackIndex));
758 }
759
760 return track ;
761
762 }
763
764 // Normal case, get info from ESD or AOD
765 // ESDs
766 if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
767 {
768 Int_t iESDtrack = cluster->GetTrackMatchedIndex();
769
770 if(iESDtrack < 0 )
771 {
772 printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d\n", index);
773 return 0x0;
774 }
775
776 track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
777
778 }
779 else // AODs
780 {
781 if(cluster->GetNTracksMatched() > 0 )
782 track = dynamic_cast<AliVTrack*>(cluster->GetTrackMatched(index));
783 }
784
785 return track ;
786
787}
788
cb5780f4 789//_____________________________________________________________________________________________________
765d44e7 790Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
791{
792 //Get the EMCAL/PHOS module number that corresponds to this particle
793
794 Int_t absId = -1;
795 if(particle->GetDetector()=="EMCAL"){
796 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
797 if(fDebug > 2)
798 printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
7cc7d3f8 799 particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
765d44e7 800 return fEMCALGeo->GetSuperModuleNumber(absId) ;
801 }//EMCAL
55d66f31 802 else if(particle->GetDetector()=="PHOS")
803 {
46a3cde6 804 // In case we use the MC reader, the input are TParticles,
805 // in this case use the corresponing method in PHOS Geometry to get the particle.
806 if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
807 {
808 Int_t mod =-1;
809 Double_t z = 0., x=0.;
810 TParticle* primary = 0x0;
811 AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
812 if(stack) {
813 primary = stack->Particle(particle->GetCaloLabel(0));
814 }
815 else {
280e6766 816 Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
46a3cde6 817 }
898c9d44 818
819 if(primary){
820 fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
821 }
822 else{
823 Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
824 }
46a3cde6 825 return mod;
826 }
827 // Input are ESDs or AODs, get the PHOS module number like this.
828 else{
2206e918 829 //FIXME
830 //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
831 //return GetModuleNumber(cluster);
832 //MEFIX
833 return -1;
46a3cde6 834 }
765d44e7 835 }//PHOS
836
837 return -1;
838}
839
cb5780f4 840//_____________________________________________________________________
c8fe2783 841Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
765d44e7 842{
280e6766 843 //Get the EMCAL/PHOS module number that corresponds to this cluster
765d44e7 844 TLorentzVector lv;
845 Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
55d66f31 846 if(!cluster)
847 {
9cbbc28b 848 if(fDebug > 1) printf("AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
849 return -1;
850 }
55d66f31 851
765d44e7 852 cluster->GetMomentum(lv,v);
853 Float_t phi = lv.Phi();
854 if(phi < 0) phi+=TMath::TwoPi();
855 Int_t absId = -1;
c8fe2783 856 if(cluster->IsEMCAL()){
765d44e7 857 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
858 if(fDebug > 2)
280e6766 859 printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
7cc7d3f8 860 lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
765d44e7 861 return fEMCALGeo->GetSuperModuleNumber(absId) ;
862 }//EMCAL
55d66f31 863 else if(cluster->IsPHOS())
864 {
765d44e7 865 Int_t relId[4];
55d66f31 866 if ( cluster->GetNCells() > 0)
867 {
765d44e7 868 absId = cluster->GetCellAbsId(0);
869 if(fDebug > 2)
280e6766 870 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
7cc7d3f8 871 lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
765d44e7 872 }
873 else return -1;
874
55d66f31 875 if ( absId >= 0)
876 {
765d44e7 877 fPHOSGeo->AbsToRelNumbering(absId,relId);
878 if(fDebug > 2)
280e6766 879 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
765d44e7 880 return relId[0]-1;
881 }
882 else return -1;
883 }//PHOS
884
885 return -1;
886}
887
cb5780f4 888//___________________________________________________________________________________________________
889Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(const Int_t absId, const TString calo,
890 Int_t & icol, Int_t & irow, Int_t & iRCU) const
765d44e7 891{
892 //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
893 Int_t imod = -1;
55d66f31 894 if ( absId >= 0)
895 {
896 if(calo=="EMCAL")
897 {
765d44e7 898 Int_t iTower = -1, iIphi = -1, iIeta = -1;
899 fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
900 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
55d66f31 901 if(imod < 0 || irow < 0 || icol < 0 )
902 {
280e6766 903 Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
904 }
905
765d44e7 906 //RCU0
55d66f31 907 if(imod < 10 )
908 {
909 if (0<=irow&&irow<8) iRCU=0; // first cable row
910 else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
911 //second cable row
912 //RCU1
913 else if (8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
914 //second cable row
915 else if (16<=irow&&irow<24) iRCU=1; // third cable row
916
917 if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
918 }
919 else
920 {
921 // Last 2 SM have one single SRU, just assign RCU 0
922 iRCU = 0 ;
923 }
924
925 if (iRCU<0)
926 {
280e6766 927 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
765d44e7 928 }
929
930 return imod ;
55d66f31 931
765d44e7 932 }//EMCAL
55d66f31 933 else //PHOS
934 {
765d44e7 935 Int_t relId[4];
936 fPHOSGeo->AbsToRelNumbering(absId,relId);
937 irow = relId[2];
938 icol = relId[3];
939 imod = relId[0]-1;
940 iRCU= (Int_t)(relId[2]-1)/16 ;
941 //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
55d66f31 942 if (iRCU >= 4)
943 {
280e6766 944 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
765d44e7 945 }
946 return imod;
947 }//PHOS
948 }
949
950 return -1;
951}
952
7db7dcb6 953//___________________________________________________________________________________________
71e3889f 954Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells)
955{
956 // Find local maxima in cluster
957
958 const Int_t nc = cluster->GetNCells();
2bf17171 959 Int_t absIdList[nc];
960 Float_t maxEList[nc];
71e3889f 961
962 Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList);
963
71e3889f 964 return nMax;
965
966}
967
968//___________________________________________________________________________________________
7db7dcb6 969Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
970 Int_t *absIdList, Float_t *maxEList)
971{
972 // Find local maxima in cluster
f241ecc3 973
7db7dcb6 974 Int_t iDigitN = 0 ;
975 Int_t iDigit = 0 ;
976 Int_t absId1 = -1 ;
977 Int_t absId2 = -1 ;
978 const Int_t nCells = cluster->GetNCells();
979
980 TString calorimeter = "EMCAL";
981 if(!cluster->IsEMCAL()) calorimeter = "PHOS";
982
983 //printf("cluster : ncells %d \n",nCells);
f241ecc3 984
985 Float_t emax = 0;
986 Int_t idmax =-1;
9369a2b1 987 for(iDigit = 0; iDigit < nCells ; iDigit++)
988 {
7db7dcb6 989 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit] ;
f241ecc3 990 Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
991 RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
992 if( en > emax )
993 {
994 emax = en ;
995 idmax = absIdList[iDigit] ;
996 }
997 //Int_t icol = -1, irow = -1, iRCU = -1;
998 //Int_t sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
999 //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
7db7dcb6 1000 }
1001
9369a2b1 1002 for(iDigit = 0 ; iDigit < nCells; iDigit++)
1003 {
1004 if(absIdList[iDigit]>=0)
1005 {
1006 absId1 = cluster->GetCellsAbsId()[iDigit];
7db7dcb6 1007
1008 Float_t en1 = cells->GetCellAmplitude(absId1);
1009 RecalibrateCellAmplitude(en1,calorimeter,absId1);
f241ecc3 1010
9369a2b1 1011 //printf("%d : absIDi %d, E %f\n",iDigit, absId1,en1);
7db7dcb6 1012
9369a2b1 1013 for(iDigitN = 0; iDigitN < nCells; iDigitN++)
1014 {
1015 absId2 = cluster->GetCellsAbsId()[iDigitN] ;
7db7dcb6 1016
9369a2b1 1017 if(absId2==-1 || absId2==absId1) continue;
7db7dcb6 1018
9369a2b1 1019 //printf("\t %d : absIDj %d\n",iDigitN, absId2);
7db7dcb6 1020
1021 Float_t en2 = cells->GetCellAmplitude(absId2);
1022 RecalibrateCellAmplitude(en2,calorimeter,absId2);
f241ecc3 1023
9369a2b1 1024 //printf("\t %d : absIDj %d, E %f\n",iDigitN, absId2,en2);
7db7dcb6 1025
9369a2b1 1026 if ( AreNeighbours(calorimeter, absId1, absId2) )
1027 {
1028 // printf("\t \t Neighbours \n");
1029 if (en1 > en2 )
1030 {
7db7dcb6 1031 absIdList[iDigitN] = -1 ;
1032 //printf("\t \t indexN %d not local max\n",iDigitN);
1033 // but may be digit too is not local max ?
1034 if(en1 < en2 + fLocMaxCutEDiff) {
9369a2b1 1035 //printf("\t \t index %d not local max cause locMaxCutEDiff\n",iDigit);
7db7dcb6 1036 absIdList[iDigit] = -1 ;
1037 }
1038 }
9369a2b1 1039 else
1040 {
7db7dcb6 1041 absIdList[iDigit] = -1 ;
1042 //printf("\t \t index %d not local max\n",iDigitN);
1043 // but may be digitN too is not local max ?
1044 if(en1 > en2 - fLocMaxCutEDiff)
1045 {
1046 absIdList[iDigitN] = -1 ;
9369a2b1 1047 //printf("\t \t indexN %d not local max cause locMaxCutEDiff\n",iDigit);
7db7dcb6 1048 }
1049 }
9369a2b1 1050 } // if Are neighbours
1051 //else printf("\t \t NOT Neighbours \n");
7db7dcb6 1052 } // while digitN
1053 } // slot not empty
1054 } // while digit
1055
1056 iDigitN = 0 ;
9369a2b1 1057 for(iDigit = 0; iDigit < nCells; iDigit++)
1058 {
1059 if(absIdList[iDigit]>=0 )
1060 {
7db7dcb6 1061 absIdList[iDigitN] = absIdList[iDigit] ;
1062 Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1063 RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1064 if(en < fLocMaxCutE) continue; // Maxima only with seed energy at least
1065 maxEList[iDigitN] = en ;
1066 //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
1067 iDigitN++ ;
1068 }
1069 }
1070
f241ecc3 1071 if(iDigitN == 0)
1072 {
1073 if(fDebug > 0)
1074 printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - No local maxima found, assign highest energy cell as maxima, id %d, en cell %2.2f, en cluster %2.2f\n",
1075 idmax,emax,cluster->E());
1076 iDigitN = 1 ;
1077 maxEList[0] = emax ;
1078 absIdList[0] = idmax ;
1079 }
1080
1081 if(fDebug > 0)
1082 {
1083 printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - In cluster E %2.2f, M02 %2.2f, M20 %2.2f, N maxima %d \n",
1084 cluster->E(),cluster->GetM02(),cluster->GetM20(), iDigitN);
1085
1086 if(fDebug > 1) for(Int_t imax = 0; imax < iDigitN; imax++)
1087 {
1088 printf(" \t i %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
1089 }
1090 }
7db7dcb6 1091
1092 return iDigitN ;
1093
1094}
1095
55d66f31 1096//____________________________________
1097TString AliCalorimeterUtils::GetPass()
1098{
1099 // Get passx from filename.
1100
1101 if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1102 {
1103 AliError("AliCalorimeterUtils::GetPass() - Pointer to tree = 0, returning null\n");
1104 return TString("");
1105 }
1106
1107 if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1108 {
1109 AliError("AliCalorimeterUtils::GetPass() - Null pointer input file, returning null\n");
1110 return TString("");
1111 }
1112
1113 TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1114 if (pass.Contains("ass1")) return TString("pass1");
1115 else if (pass.Contains("ass2")) return TString("pass2");
1116 else if (pass.Contains("ass3")) return TString("pass3");
1117
1118 // No condition fullfilled, give a default value
1119 printf("AliCalorimeterUtils::GetPass() - Pass number string not found \n");
1120 return TString("");
1121
1122}
1123
cb5780f4 1124//________________________________________
765d44e7 1125void AliCalorimeterUtils::InitParameters()
1126{
1127 //Initialize the parameters of the analysis.
7db7dcb6 1128
12e91266 1129 fEMCALGeoName = "";
1130 fPHOSGeoName = "";
7db7dcb6 1131
7cc7d3f8 1132 fEMCALGeoMatrixSet = kFALSE;
1133 fPHOSGeoMatrixSet = kFALSE;
7db7dcb6 1134
7cc7d3f8 1135 fRemoveBadChannels = kFALSE;
7db7dcb6 1136
7cc7d3f8 1137 fNCellsFromPHOSBorder = 0;
a5fb4114 1138
7db7dcb6 1139 fLocMaxCutE = 0.1 ;
1140 fLocMaxCutEDiff = 0.0 ;
1141
7cc7d3f8 1142 // fMaskCellColumns = new Int_t[fNMaskCellColumns];
1143 // fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
1144 // fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
1145 // fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
1146 // fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
1147 // fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
a5fb4114 1148
55d66f31 1149 fOADBSet = kFALSE;
1150 fOADBForEMCAL = kTRUE ;
1151 fOADBForPHOS = kFALSE;
1152
1153 fOADBFilePathEMCAL = "$ALICE_ROOT/OADB/EMCAL" ;
1154 fOADBFilePathPHOS = "$ALICE_ROOT/OADB/PHOS" ;
1155
765d44e7 1156}
1157
765d44e7 1158
cb5780f4 1159//_____________________________________________________
1160void AliCalorimeterUtils::InitPHOSBadChannelStatusMap()
1161{
765d44e7 1162 //Init PHOS bad channels map
1163 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
1164 //In order to avoid rewriting the same histograms
1165 Bool_t oldStatus = TH1::AddDirectoryStatus();
1166 TH1::AddDirectory(kFALSE);
7cc7d3f8 1167
78219bac 1168 fPHOSBadChannelMap = new TObjArray(5);
c5693f62 1169 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 1170
765d44e7 1171 fPHOSBadChannelMap->SetOwner(kTRUE);
1172 fPHOSBadChannelMap->Compress();
1173
1174 //In order to avoid rewriting the same histograms
1175 TH1::AddDirectory(oldStatus);
1176}
1177
cb5780f4 1178//______________________________________________________
1179void AliCalorimeterUtils::InitPHOSRecalibrationFactors()
1180{
09e819c9 1181 //Init EMCAL recalibration factors
1182 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
1183 //In order to avoid rewriting the same histograms
1184 Bool_t oldStatus = TH1::AddDirectoryStatus();
1185 TH1::AddDirectory(kFALSE);
7cc7d3f8 1186
78219bac 1187 fPHOSRecalibrationFactors = new TObjArray(5);
c5693f62 1188 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 1189 //Init the histograms with 1
1190 for (Int_t m = 0; m < 5; m++) {
1191 for (Int_t i = 0; i < 56; i++) {
1192 for (Int_t j = 0; j < 64; j++) {
c5693f62 1193 SetPHOSChannelRecalibrationFactor(m,j,i,1.);
09e819c9 1194 }
1195 }
1196 }
1197 fPHOSRecalibrationFactors->SetOwner(kTRUE);
1198 fPHOSRecalibrationFactors->Compress();
1199
1200 //In order to avoid rewriting the same histograms
1201 TH1::AddDirectory(oldStatus);
1202}
1203
1204
55d66f31 1205//__________________________________________________________
1206void AliCalorimeterUtils::InitEMCALGeometry(Int_t runnumber)
765d44e7 1207{
1208 //Initialize EMCAL geometry if it did not exist previously
55d66f31 1209
1210 if (!fEMCALGeo)
1211 {
1212 if(fEMCALGeoName=="")
1213 {
5408e59e 1214 if (runnumber < 140000 &&
1215 runnumber >= 100000) fEMCALGeoName = "EMCAL_FIRSTYEARV1";
1216 else if(runnumber >= 140000 &&
1217 runnumber < 171000) fEMCALGeoName = "EMCAL_COMPLETEV1";
1218 else fEMCALGeoName = "EMCAL_COMPLETE12SMV1";
55d66f31 1219 printf("AliCalorimeterUtils::InitEMCALGeometry() - Set EMCAL geometry name to <%s> for run %d\n",fEMCALGeoName.Data(),runnumber);
1220 }
1221
a38a48f2 1222 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
55d66f31 1223
1224 if(fDebug > 0)
1225 {
1226 printf("AliCalorimeterUtils::InitEMCALGeometry(run=%d)",runnumber);
765d44e7 1227 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1228 printf("\n");
1229 }
1230 }
1231}
1232
55d66f31 1233//_________________________________________________________
1234void AliCalorimeterUtils::InitPHOSGeometry(Int_t runnumber)
765d44e7 1235{
1236 //Initialize PHOS geometry if it did not exist previously
55d66f31 1237
1238 if (!fPHOSGeo)
1239 {
1240 if(fPHOSGeoName=="") fPHOSGeoName = "PHOSgeo";
1241
765d44e7 1242 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
55d66f31 1243
1244 if(fDebug > 0)
1245 {
1246 printf("AliCalorimeterUtils::InitPHOSGeometry(run=%d)",runnumber);
765d44e7 1247 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1248 printf("\n");
1249 }
1250 }
1251}
1252
55d66f31 1253//__________________________________________________________________
1254Bool_t AliCalorimeterUtils::MaskFrameCluster(const Int_t iSM,
1255 const Int_t ieta) const
1256{
1257 //Check if cell is in one of the regions where we have significant amount
1258 //of material in front. Only EMCAL
1259
1260 Int_t icol = ieta;
1261 if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
1262
1263 if (fNMaskCellColumns && fMaskCellColumns)
1264 {
1265 for (Int_t imask = 0; imask < fNMaskCellColumns; imask++)
1266 {
1267 if(icol==fMaskCellColumns[imask]) return kTRUE;
1268 }
1269 }
1270
1271 return kFALSE;
1272
1273}
1274
cb5780f4 1275//_________________________________________________________
765d44e7 1276void AliCalorimeterUtils::Print(const Option_t * opt) const
1277{
7cc7d3f8 1278
765d44e7 1279 //Print some relevant parameters set for the analysis
1280 if(! opt)
1281 return;
7cc7d3f8 1282
765d44e7 1283 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1284 printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
1285 printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
7cc7d3f8 1286 fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
247abff4 1287 if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
09e819c9 1288 printf("Recalibrate Clusters? %d\n",fRecalibration);
9584c261 1289 printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
1290 printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
f2ccb5b8 1291 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
7cc7d3f8 1292
7db7dcb6 1293 printf("Loc. Max. E > %2.2f\n", fLocMaxCutE);
1294 printf("Loc. Max. E Diff > %2.2f\n", fLocMaxCutEDiff);
1295
765d44e7 1296 printf(" \n") ;
1297}
1298
9369a2b1 1299//__________________________________________________________________________________________
7db7dcb6 1300void AliCalorimeterUtils::RecalibrateCellAmplitude(Float_t & amp,
9369a2b1 1301 const TString calo, const Int_t id) const
7db7dcb6 1302{
1303 //Recaculate cell energy if recalibration factor
1304
1305 Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
1306 Int_t nModule = GetModuleNumberCellIndexes(id,calo, icol, irow, iRCU);
1307
1308 if (IsRecalibrationOn())
1309 {
1310 if(calo == "PHOS")
1311 {
1312 amp *= GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
1313 }
1314 else
1315 {
1316 amp *= GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
1317 }
1318 }
1319}
1320
9369a2b1 1321//_________________________________________________________________________________
7db7dcb6 1322void AliCalorimeterUtils::RecalibrateCellTime(Double_t & time,
1323 const TString calo,
9369a2b1 1324 const Int_t id, const Int_t bc) const
7db7dcb6 1325{
1326 // Recalculate time if time recalibration available for EMCAL
1327 // not ready for PHOS
1328
1329 if(calo == "EMCAL" && GetEMCALRecoUtils()->IsTimeRecalibrationOn())
1330 {
1331 GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
1332 }
1333
1334}
1335
cb5780f4 1336//__________________________________________________________________________
1337Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster,
1338 AliVCaloCells * cells)
1339{
09e819c9 1340 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
7cc7d3f8 1341
898c9d44 1342 //Initialize some used variables
7db7dcb6 1343 Float_t frac = 0., energy = 0.;
898c9d44 1344
7db7dcb6 1345 if(cells)
1346 {
7cc7d3f8 1347 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
7db7dcb6 1348
7cc7d3f8 1349 UShort_t * index = cluster->GetCellsAbsId() ;
1350 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
7db7dcb6 1351
7cc7d3f8 1352 Int_t ncells = cluster->GetNCells();
7db7dcb6 1353
7cc7d3f8 1354 TString calo = "EMCAL";
7db7dcb6 1355 if(cluster->IsPHOS())
1356 calo = "PHOS";
7cc7d3f8 1357
1358 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1359 for(Int_t icell = 0; icell < ncells; icell++){
7db7dcb6 1360
1361 Int_t absId = index[icell];
1362
7cc7d3f8 1363 frac = fraction[icell];
1364 if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
7db7dcb6 1365
1366 Float_t amp = cells->GetCellAmplitude(absId);
1367 RecalibrateCellAmplitude(amp,calo, absId);
1368
7cc7d3f8 1369 if(fDebug>2)
7db7dcb6 1370 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, cell fraction %f, cell energy %f\n",
1371 calo.Data(),frac,cells->GetCellAmplitude(absId));
7cc7d3f8 1372
7db7dcb6 1373 energy += amp*frac;
7cc7d3f8 1374 }
1375
1376 if(fDebug>1)
1377 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
898c9d44 1378
1379 }// cells available
7db7dcb6 1380 else
1381 {
898c9d44 1382 Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
1383 }
1384
09e819c9 1385 return energy;
09e819c9 1386}
1387
cb5780f4 1388//__________________________________________________________________________________________
1389void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
1390{
9584c261 1391
1392 //Recalculate EMCAL cluster position
1393
19db8f8c 1394 fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
7cc7d3f8 1395
9584c261 1396}
cb5780f4 1397
55d66f31 1398//________________________________________________________________________________
1399void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event,
1400 TObjArray* clusterArray)
1401{
1402 //Recalculate track matching
1403
1404 if (fRecalculateMatching)
1405 {
1406 fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo) ;
1407 //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1408 //if(esdevent){
1409 // fEMCALRecoUtils->SetClusterMatchedToTrack(esdevent) ;
1410 // fEMCALRecoUtils->SetTracksMatchedToCluster(esdevent) ;
1411 //}
1412 }
1413}
1414
3c1d9afb 1415//___________________________________________________________________________
1416void AliCalorimeterUtils::SplitEnergy(const Int_t absId1, const Int_t absId2,
1417 AliVCluster* cluster,
1418 AliVCaloCells* cells,
1419 //Float_t & e1, Float_t & e2,
1420 AliAODCaloCluster* cluster1,
1421 AliAODCaloCluster* cluster2,
1422 const Int_t nMax,
1423 const Int_t eventNumber)
1424{
1425
1426 // Split energy of cluster between the 2 local maxima, sum energy on 3x3, and if the 2
1427 // maxima are too close and have common cells, split the energy between the 2
1428
1429 TH2F* hClusterMap = 0 ;
1430 TH2F* hClusterLocMax = 0 ;
1431 TH2F* hCluster1 = 0 ;
1432 TH2F* hCluster2 = 0 ;
1433
1434 if(fPlotCluster)
1435 {
1436 hClusterMap = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
1437 hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
1438 hCluster1 = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
1439 hCluster2 = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
1440 hClusterMap ->SetXTitle("column");
1441 hClusterMap ->SetYTitle("row");
1442 hClusterLocMax ->SetXTitle("column");
1443 hClusterLocMax ->SetYTitle("row");
1444 hCluster1 ->SetXTitle("column");
1445 hCluster1 ->SetYTitle("row");
1446 hCluster2 ->SetXTitle("column");
1447 hCluster2 ->SetYTitle("row");
1448 }
1449
1450 TString calorimeter = "EMCAL";
1451 if(cluster->IsPHOS())
1452 {
1453 calorimeter="PHOS";
1454 printf("AliCalorimeterUtils::SplitEnerg() Not supported for PHOS yet \n");
1455 return;
1456 }
1457
1458 const Int_t ncells = cluster->GetNCells();
1459 Int_t absIdList[ncells];
1460
1461 Float_t e1 = 0, e2 = 0 ;
1462 Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;
1463 Float_t eCluster = 0;
1464 Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1;
1465 for(Int_t iDigit = 0; iDigit < ncells; iDigit++ )
1466 {
1467 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
1468
1469 Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
1470 RecalibrateCellAmplitude(ec,calorimeter, absIdList[iDigit]);
1471 eCluster+=ec;
1472
1473 if(fPlotCluster)
1474 {
1475 //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
1476 sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1477 if(icol > maxCol) maxCol = icol;
1478 if(icol < minCol) minCol = icol;
1479 if(irow > maxRow) maxRow = irow;
1480 if(irow < minRow) minRow = irow;
1481 hClusterMap->Fill(icol,irow,ec);
1482 }
1483
1484 }
1485
1486 // Init counters and variables
1487 Int_t ncells1 = 1 ;
1488 UShort_t absIdList1[9] ;
1489 Double_t fracList1 [9] ;
1490 absIdList1[0] = absId1 ;
1491 fracList1 [0] = 1. ;
1492
1493 Float_t ecell1 = cells->GetCellAmplitude(absId1);
1494 RecalibrateCellAmplitude(ecell1, calorimeter, absId1);
1495 e1 = ecell1;
1496
1497 Int_t ncells2 = 1 ;
1498 UShort_t absIdList2[9] ;
1499 Double_t fracList2 [9] ;
1500 absIdList2[0] = absId2 ;
1501 fracList2 [0] = 1. ;
1502
1503 Float_t ecell2 = cells->GetCellAmplitude(absId2);
1504 RecalibrateCellAmplitude(ecell2, calorimeter, absId2);
1505 e2 = ecell2;
1506
1507 if(fPlotCluster)
1508 {
1509 Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
1510 sm = GetModuleNumberCellIndexes(absId1, calorimeter, icol1, irow1, iRCU) ;
1511 hClusterLocMax->Fill(icol1,irow1,ecell1);
1512 sm = GetModuleNumberCellIndexes(absId2, calorimeter, icol2, irow2, iRCU) ;
1513 hClusterLocMax->Fill(icol2,irow2,ecell2);
1514 }
1515
1516 // Very rough way to share the cluster energy
1517 Float_t eRemain = (eCluster-ecell1-ecell2)/2;
1518 Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
1519 Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
1520
1521 for(Int_t iDigit = 0; iDigit < ncells; iDigit++)
1522 {
1523 Int_t absId = absIdList[iDigit];
1524
1525 if(absId==absId1 || absId==absId2 || absId < 0) continue;
1526
1527 Float_t ecell = cells->GetCellAmplitude(absId);
1528 RecalibrateCellAmplitude(ecell, calorimeter, absId);
1529
1530 if(AreNeighbours(calorimeter, absId1,absId ))
1531 {
1532 absIdList1[ncells1]= absId;
1533
1534 if(AreNeighbours(calorimeter, absId2,absId ))
1535 {
1536 fracList1[ncells1] = shareFraction1;
1537 e1 += ecell*shareFraction1;
1538 }
1539 else
1540 {
1541 fracList1[ncells1] = 1.;
1542 e1 += ecell;
1543 }
1544
1545 ncells1++;
1546
1547 } // neigbour to cell1
1548
1549 if(AreNeighbours(calorimeter, absId2,absId ))
1550 {
1551 absIdList2[ncells2]= absId;
1552
1553 if(AreNeighbours(calorimeter, absId1,absId ))
1554 {
1555 fracList2[ncells2] = shareFraction2;
1556 e2 += ecell*shareFraction2;
1557 }
1558 else
1559 {
1560 fracList2[ncells2] = 1.;
1561 e2 += ecell;
1562 }
1563
1564 ncells2++;
1565
1566 } // neigbour to cell2
1567
1568 }
1569
1570 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",
1571 nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2);
1572
1573 cluster1->SetE(e1);
1574 cluster2->SetE(e2);
1575
1576 cluster1->SetNCells(ncells1);
1577 cluster2->SetNCells(ncells2);
1578
1579 cluster1->SetCellsAbsId(absIdList1);
1580 cluster2->SetCellsAbsId(absIdList2);
1581
1582 cluster1->SetCellsAmplitudeFraction(fracList1);
1583 cluster2->SetCellsAmplitudeFraction(fracList2);
1584
1585 if(calorimeter=="EMCAL")
1586 {
1587 GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
1588 GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
1589 }
1590
1591 if(fPlotCluster)
1592 {
1593 //printf("Cells of cluster1: ");
1594 for(Int_t iDigit = 0; iDigit < ncells1; iDigit++ )
1595 {
1596 //printf(" %d ",absIdList1[iDigit]);
1597
1598 sm = GetModuleNumberCellIndexes(absIdList1[iDigit], calorimeter, icol, irow, iRCU) ;
1599
1600 if( AreNeighbours(calorimeter, absId2,absIdList1[iDigit]) )
1601 hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit])*shareFraction1);
1602 else
1603 hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit]));
1604 }
1605
1606 //printf(" \n ");
1607 //printf("Cells of cluster2: ");
1608
1609 for(Int_t iDigit = 0; iDigit < ncells2; iDigit++ )
1610 {
1611 //printf(" %d ",absIdList2[iDigit]);
1612
1613 sm = GetModuleNumberCellIndexes(absIdList2[iDigit], calorimeter, icol, irow, iRCU) ;
1614 if( AreNeighbours(calorimeter, absId1,absIdList2[iDigit]) )
1615 hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit])*shareFraction2);
1616 else
1617 hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit]));
1618
1619 }
1620 //printf(" \n ");
1621
1622 gStyle->SetPadRightMargin(0.1);
1623 gStyle->SetPadLeftMargin(0.1);
1624 gStyle->SetOptStat(0);
1625 gStyle->SetOptFit(000000);
1626
1627 if(maxCol-minCol > maxRow-minRow)
1628 {
1629 maxRow+= maxCol-minCol;
1630 }
1631 else
1632 {
1633 maxCol+= maxRow-minRow;
1634 }
1635
1636 TCanvas * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
1637 c->Divide(2,2);
1638 c->cd(1);
1639 gPad->SetGridy();
1640 gPad->SetGridx();
1641 hClusterMap ->SetAxisRange(minCol, maxCol,"X");
1642 hClusterMap ->SetAxisRange(minRow, maxRow,"Y");
1643 hClusterMap ->Draw("colz");
1644 c->cd(2);
1645 gPad->SetGridy();
1646 gPad->SetGridx();
1647 hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
1648 hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
1649 hClusterLocMax ->Draw("colz");
1650 c->cd(3);
1651 gPad->SetGridy();
1652 gPad->SetGridx();
1653 hCluster1 ->SetAxisRange(minCol, maxCol,"X");
1654 hCluster1 ->SetAxisRange(minRow, maxRow,"Y");
1655 hCluster1 ->Draw("colz");
1656 c->cd(4);
1657 gPad->SetGridy();
1658 gPad->SetGridx();
1659 hCluster2 ->SetAxisRange(minCol, maxCol,"X");
1660 hCluster2 ->SetAxisRange(minRow, maxRow,"Y");
1661 hCluster2 ->Draw("colz");
1662
1663 if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",
1664 eventNumber,cluster->E(),nMax,ncells1,ncells2));
1665
1666 delete c;
1667 delete hClusterMap;
1668 delete hClusterLocMax;
1669 delete hCluster1;
1670 delete hCluster2;
1671 }
1672}
1673