revert previous commit on this class
[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
297ee714 580 {
5408e59e 581 if(iphi >= nborder && iphi < 8-nborder) okrow =kTRUE;
297ee714 582 }
5408e59e 583 else // 1/2 SM
297ee714 584 {
5408e59e 585 if(iphi >= nborder && iphi <12-nborder) okrow =kTRUE;
297ee714 586 }
765d44e7 587 }
588
247abff4 589 //Check columns/eta
5408e59e 590 if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0())
591 {
247abff4 592 if(ieta > nborder && ieta < 48-nborder) okcol =kTRUE;
765d44e7 593 }
5408e59e 594 else
595 {
596 if(iSM%2==0)
597 {
247abff4 598 if(ieta >= nborder) okcol = kTRUE;
765d44e7 599 }
5408e59e 600 else
601 {
247abff4 602 if(ieta < 48-nborder) okcol = kTRUE;
765d44e7 603 }
604 }//eta 0 not checked
5408e59e 605
765d44e7 606 if(fDebug > 1)
607 {
280e6766 608 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
7cc7d3f8 609 nborder, ieta, iphi, iSM);
765d44e7 610 if (okcol && okrow ) printf(" YES \n");
611 else printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
612 }
613 }//EMCAL
c8fe2783 614 else if(cells->GetType()==AliVCaloCells::kPHOSCell){
765d44e7 615 Int_t relId[4];
616 Int_t irow = -1, icol = -1;
617 fPHOSGeo->AbsToRelNumbering(absIdMax,relId);
618 irow = relId[2];
619 icol = relId[3];
620 //imod = relId[0]-1;
621 if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE;
622 if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE;
623 if(fDebug > 1)
624 {
280e6766 625 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - PHOS Cluster in %d cells fiducial volume: icol %d, irow %d, Module %d?",
7cc7d3f8 626 fNCellsFromPHOSBorder, icol, irow, relId[0]-1);
765d44e7 627 if (okcol && okrow ) printf(" YES \n");
628 else printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
629 }
630 }//PHOS
631
632 if (okcol && okrow) return kTRUE;
633 else return kFALSE;
634
635}
636
765d44e7 637//_________________________________________________________________________________________________________
cb5780f4 638Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter,UShort_t* cellList, Int_t nCells)
639{
765d44e7 640 // Check that in the cluster cells, there is no bad channel of those stored
641 // in fEMCALBadChannelMap or fPHOSBadChannelMap
642
643 if (!fRemoveBadChannels) return kFALSE;
36037088 644 //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
247abff4 645 if(calorimeter == "EMCAL" && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
78219bac 646 if(calorimeter == "PHOS" && !fPHOSBadChannelMap) return kFALSE;
7cc7d3f8 647
765d44e7 648 Int_t icol = -1;
649 Int_t irow = -1;
650 Int_t imod = -1;
651 for(Int_t iCell = 0; iCell<nCells; iCell++){
7cc7d3f8 652
765d44e7 653 //Get the column and row
654 if(calorimeter == "EMCAL"){
247abff4 655 return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
765d44e7 656 }
657 else if(calorimeter=="PHOS"){
658 Int_t relId[4];
659 fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
660 irow = relId[2];
661 icol = relId[3];
662 imod = relId[0]-1;
663 if(fPHOSBadChannelMap->GetEntries() <= imod)continue;
c5693f62 664 //printf("PHOS bad channels imod %d, icol %d, irow %d\n",imod, irow, icol);
665 if(GetPHOSChannelStatus(imod, irow, icol)) return kTRUE;
765d44e7 666 }
667 else return kFALSE;
668
669 }// cell cluster loop
7cc7d3f8 670
765d44e7 671 return kFALSE;
7cc7d3f8 672
765d44e7 673}
674
cb5780f4 675//_______________________________________________________________
676void AliCalorimeterUtils::CorrectClusterEnergy(AliVCluster *clus)
677{
9584c261 678 // Correct cluster energy non linearity
13cd2872 679
9584c261 680 clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
13cd2872 681
682}
683
c5693f62 684//________________________________________________________________________________________
685Int_t AliCalorimeterUtils::GetMaxEnergyCell(AliVCaloCells* cells, const AliVCluster* clu,
686 Float_t & clusterFraction) const
cb5780f4 687{
13cd2872 688
689 //For a given CaloCluster gets the absId of the cell
690 //with maximum energy deposit.
691
692 if( !clu || !cells ){
693 AliInfo("Cluster or cells pointer is null!");
694 return -1;
695 }
696
697 Double_t eMax =-1.;
698 Double_t eTot = 0.;
699 Double_t eCell =-1.;
700 Float_t fraction = 1.;
701 Float_t recalFactor = 1.;
702 Int_t cellAbsId =-1 , absId =-1 ;
703 Int_t iSupMod =-1 , ieta =-1 , iphi = -1, iRCU = -1;
704
705 TString calo = "EMCAL";
706 if(clu->IsPHOS()) calo = "PHOS";
707
708 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
709
710 cellAbsId = clu->GetCellAbsId(iDig);
711
712 fraction = clu->GetCellAmplitudeFraction(iDig);
713 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
714
b08dd6d5 715 iSupMod = GetModuleNumberCellIndexes(cellAbsId, calo, ieta, iphi, iRCU);
13cd2872 716
717 if(IsRecalibrationOn()) {
718 if(calo=="EMCAL") recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
c5693f62 719 else recalFactor = GetPHOSChannelRecalibrationFactor (iSupMod,iphi,ieta);
13cd2872 720 }
721
722 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
723
724 if(eCell > eMax) {
725 eMax = eCell;
726 absId = cellAbsId;
727 }
728
729 eTot+=eCell;
730
731 }// cell loop
732
733 clusterFraction = (eTot-eMax)/eTot; //Do not use cluster energy in case it was corrected for non linearity.
734
735 return absId;
736
9584c261 737}
738
d832b695 739//__________________________________________________________________________
740AliVTrack * AliCalorimeterUtils::GetMatchedTrack(const AliVCluster* cluster,
741 const AliVEvent* event,
742 const Int_t index) const
743{
744 // Get the matched track given its index, usually just the first match
745 // Since it is different for ESDs and AODs here it is a wrap method to do it
746
747 AliVTrack *track = 0;
748
749 // EMCAL case only when matching is recalculated
750 if(cluster->IsEMCAL() && IsRecalculationOfClusterTrackMatchingOn())
751 {
752 Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(cluster->GetID());
753 //printf("track index %d, cluster ID %d \n ",trackIndex,cluster->GetID());
754
755 if(trackIndex < 0 )
756 {
757 printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d, from recalculation\n", trackIndex);
758 }
759 else
760 {
761 track = dynamic_cast<AliVTrack*> (event->GetTrack(trackIndex));
762 }
763
764 return track ;
765
766 }
767
768 // Normal case, get info from ESD or AOD
769 // ESDs
770 if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
771 {
772 Int_t iESDtrack = cluster->GetTrackMatchedIndex();
773
774 if(iESDtrack < 0 )
775 {
776 printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d\n", index);
777 return 0x0;
778 }
779
780 track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
781
782 }
783 else // AODs
784 {
785 if(cluster->GetNTracksMatched() > 0 )
786 track = dynamic_cast<AliVTrack*>(cluster->GetTrackMatched(index));
787 }
788
789 return track ;
790
791}
792
cb5780f4 793//_____________________________________________________________________________________________________
765d44e7 794Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
795{
796 //Get the EMCAL/PHOS module number that corresponds to this particle
797
798 Int_t absId = -1;
799 if(particle->GetDetector()=="EMCAL"){
800 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
801 if(fDebug > 2)
802 printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
7cc7d3f8 803 particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
765d44e7 804 return fEMCALGeo->GetSuperModuleNumber(absId) ;
805 }//EMCAL
55d66f31 806 else if(particle->GetDetector()=="PHOS")
807 {
46a3cde6 808 // In case we use the MC reader, the input are TParticles,
809 // in this case use the corresponing method in PHOS Geometry to get the particle.
810 if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
811 {
812 Int_t mod =-1;
813 Double_t z = 0., x=0.;
814 TParticle* primary = 0x0;
815 AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
816 if(stack) {
817 primary = stack->Particle(particle->GetCaloLabel(0));
818 }
819 else {
280e6766 820 Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
46a3cde6 821 }
898c9d44 822
823 if(primary){
824 fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
825 }
826 else{
827 Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
828 }
46a3cde6 829 return mod;
830 }
831 // Input are ESDs or AODs, get the PHOS module number like this.
832 else{
2206e918 833 //FIXME
834 //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
835 //return GetModuleNumber(cluster);
836 //MEFIX
837 return -1;
46a3cde6 838 }
765d44e7 839 }//PHOS
840
841 return -1;
842}
843
cb5780f4 844//_____________________________________________________________________
c8fe2783 845Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
765d44e7 846{
280e6766 847 //Get the EMCAL/PHOS module number that corresponds to this cluster
765d44e7 848 TLorentzVector lv;
849 Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
55d66f31 850 if(!cluster)
851 {
9cbbc28b 852 if(fDebug > 1) printf("AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
853 return -1;
854 }
55d66f31 855
765d44e7 856 cluster->GetMomentum(lv,v);
857 Float_t phi = lv.Phi();
858 if(phi < 0) phi+=TMath::TwoPi();
859 Int_t absId = -1;
c8fe2783 860 if(cluster->IsEMCAL()){
765d44e7 861 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
862 if(fDebug > 2)
280e6766 863 printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
7cc7d3f8 864 lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
765d44e7 865 return fEMCALGeo->GetSuperModuleNumber(absId) ;
866 }//EMCAL
55d66f31 867 else if(cluster->IsPHOS())
868 {
765d44e7 869 Int_t relId[4];
55d66f31 870 if ( cluster->GetNCells() > 0)
871 {
765d44e7 872 absId = cluster->GetCellAbsId(0);
873 if(fDebug > 2)
280e6766 874 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
7cc7d3f8 875 lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
765d44e7 876 }
877 else return -1;
878
55d66f31 879 if ( absId >= 0)
880 {
765d44e7 881 fPHOSGeo->AbsToRelNumbering(absId,relId);
882 if(fDebug > 2)
280e6766 883 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
765d44e7 884 return relId[0]-1;
885 }
886 else return -1;
887 }//PHOS
888
889 return -1;
890}
891
cb5780f4 892//___________________________________________________________________________________________________
893Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(const Int_t absId, const TString calo,
894 Int_t & icol, Int_t & irow, Int_t & iRCU) const
765d44e7 895{
896 //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
897 Int_t imod = -1;
55d66f31 898 if ( absId >= 0)
899 {
900 if(calo=="EMCAL")
901 {
765d44e7 902 Int_t iTower = -1, iIphi = -1, iIeta = -1;
903 fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
904 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
55d66f31 905 if(imod < 0 || irow < 0 || icol < 0 )
906 {
280e6766 907 Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
908 }
909
765d44e7 910 //RCU0
55d66f31 911 if(imod < 10 )
912 {
913 if (0<=irow&&irow<8) iRCU=0; // first cable row
914 else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
915 //second cable row
916 //RCU1
917 else if (8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
918 //second cable row
919 else if (16<=irow&&irow<24) iRCU=1; // third cable row
920
921 if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
922 }
923 else
924 {
925 // Last 2 SM have one single SRU, just assign RCU 0
926 iRCU = 0 ;
927 }
928
929 if (iRCU<0)
930 {
280e6766 931 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
765d44e7 932 }
933
934 return imod ;
55d66f31 935
765d44e7 936 }//EMCAL
55d66f31 937 else //PHOS
938 {
765d44e7 939 Int_t relId[4];
940 fPHOSGeo->AbsToRelNumbering(absId,relId);
941 irow = relId[2];
942 icol = relId[3];
943 imod = relId[0]-1;
944 iRCU= (Int_t)(relId[2]-1)/16 ;
945 //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
55d66f31 946 if (iRCU >= 4)
947 {
280e6766 948 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
765d44e7 949 }
950 return imod;
951 }//PHOS
952 }
953
954 return -1;
955}
956
7db7dcb6 957//___________________________________________________________________________________________
71e3889f 958Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells)
959{
960 // Find local maxima in cluster
961
962 const Int_t nc = cluster->GetNCells();
2bf17171 963 Int_t absIdList[nc];
964 Float_t maxEList[nc];
71e3889f 965
966 Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList);
967
71e3889f 968 return nMax;
969
970}
971
972//___________________________________________________________________________________________
7db7dcb6 973Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
974 Int_t *absIdList, Float_t *maxEList)
975{
976 // Find local maxima in cluster
f241ecc3 977
7db7dcb6 978 Int_t iDigitN = 0 ;
979 Int_t iDigit = 0 ;
980 Int_t absId1 = -1 ;
981 Int_t absId2 = -1 ;
982 const Int_t nCells = cluster->GetNCells();
983
984 TString calorimeter = "EMCAL";
985 if(!cluster->IsEMCAL()) calorimeter = "PHOS";
986
987 //printf("cluster : ncells %d \n",nCells);
f241ecc3 988
989 Float_t emax = 0;
990 Int_t idmax =-1;
9369a2b1 991 for(iDigit = 0; iDigit < nCells ; iDigit++)
992 {
7db7dcb6 993 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit] ;
f241ecc3 994 Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
995 RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
996 if( en > emax )
997 {
998 emax = en ;
999 idmax = absIdList[iDigit] ;
1000 }
1001 //Int_t icol = -1, irow = -1, iRCU = -1;
1002 //Int_t sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1003 //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
7db7dcb6 1004 }
1005
9369a2b1 1006 for(iDigit = 0 ; iDigit < nCells; iDigit++)
1007 {
1008 if(absIdList[iDigit]>=0)
1009 {
1010 absId1 = cluster->GetCellsAbsId()[iDigit];
7db7dcb6 1011
1012 Float_t en1 = cells->GetCellAmplitude(absId1);
1013 RecalibrateCellAmplitude(en1,calorimeter,absId1);
f241ecc3 1014
9369a2b1 1015 //printf("%d : absIDi %d, E %f\n",iDigit, absId1,en1);
7db7dcb6 1016
9369a2b1 1017 for(iDigitN = 0; iDigitN < nCells; iDigitN++)
1018 {
1019 absId2 = cluster->GetCellsAbsId()[iDigitN] ;
7db7dcb6 1020
9369a2b1 1021 if(absId2==-1 || absId2==absId1) continue;
7db7dcb6 1022
9369a2b1 1023 //printf("\t %d : absIDj %d\n",iDigitN, absId2);
7db7dcb6 1024
1025 Float_t en2 = cells->GetCellAmplitude(absId2);
1026 RecalibrateCellAmplitude(en2,calorimeter,absId2);
f241ecc3 1027
9369a2b1 1028 //printf("\t %d : absIDj %d, E %f\n",iDigitN, absId2,en2);
7db7dcb6 1029
9369a2b1 1030 if ( AreNeighbours(calorimeter, absId1, absId2) )
1031 {
1032 // printf("\t \t Neighbours \n");
1033 if (en1 > en2 )
1034 {
7db7dcb6 1035 absIdList[iDigitN] = -1 ;
1036 //printf("\t \t indexN %d not local max\n",iDigitN);
1037 // but may be digit too is not local max ?
1038 if(en1 < en2 + fLocMaxCutEDiff) {
9369a2b1 1039 //printf("\t \t index %d not local max cause locMaxCutEDiff\n",iDigit);
7db7dcb6 1040 absIdList[iDigit] = -1 ;
1041 }
1042 }
9369a2b1 1043 else
1044 {
7db7dcb6 1045 absIdList[iDigit] = -1 ;
1046 //printf("\t \t index %d not local max\n",iDigitN);
1047 // but may be digitN too is not local max ?
1048 if(en1 > en2 - fLocMaxCutEDiff)
1049 {
1050 absIdList[iDigitN] = -1 ;
9369a2b1 1051 //printf("\t \t indexN %d not local max cause locMaxCutEDiff\n",iDigit);
7db7dcb6 1052 }
1053 }
9369a2b1 1054 } // if Are neighbours
1055 //else printf("\t \t NOT Neighbours \n");
7db7dcb6 1056 } // while digitN
1057 } // slot not empty
1058 } // while digit
1059
1060 iDigitN = 0 ;
9369a2b1 1061 for(iDigit = 0; iDigit < nCells; iDigit++)
1062 {
1063 if(absIdList[iDigit]>=0 )
1064 {
7db7dcb6 1065 absIdList[iDigitN] = absIdList[iDigit] ;
1066 Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1067 RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1068 if(en < fLocMaxCutE) continue; // Maxima only with seed energy at least
1069 maxEList[iDigitN] = en ;
1070 //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
1071 iDigitN++ ;
1072 }
1073 }
1074
f241ecc3 1075 if(iDigitN == 0)
1076 {
1077 if(fDebug > 0)
1078 printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - No local maxima found, assign highest energy cell as maxima, id %d, en cell %2.2f, en cluster %2.2f\n",
1079 idmax,emax,cluster->E());
1080 iDigitN = 1 ;
1081 maxEList[0] = emax ;
1082 absIdList[0] = idmax ;
1083 }
1084
1085 if(fDebug > 0)
1086 {
1087 printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - In cluster E %2.2f, M02 %2.2f, M20 %2.2f, N maxima %d \n",
1088 cluster->E(),cluster->GetM02(),cluster->GetM20(), iDigitN);
1089
1090 if(fDebug > 1) for(Int_t imax = 0; imax < iDigitN; imax++)
1091 {
1092 printf(" \t i %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
1093 }
1094 }
7db7dcb6 1095
1096 return iDigitN ;
1097
1098}
1099
55d66f31 1100//____________________________________
1101TString AliCalorimeterUtils::GetPass()
1102{
1103 // Get passx from filename.
1104
1105 if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1106 {
1107 AliError("AliCalorimeterUtils::GetPass() - Pointer to tree = 0, returning null\n");
1108 return TString("");
1109 }
1110
1111 if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1112 {
1113 AliError("AliCalorimeterUtils::GetPass() - Null pointer input file, returning null\n");
1114 return TString("");
1115 }
1116
1117 TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1118 if (pass.Contains("ass1")) return TString("pass1");
1119 else if (pass.Contains("ass2")) return TString("pass2");
1120 else if (pass.Contains("ass3")) return TString("pass3");
1121
1122 // No condition fullfilled, give a default value
1123 printf("AliCalorimeterUtils::GetPass() - Pass number string not found \n");
1124 return TString("");
1125
1126}
1127
cb5780f4 1128//________________________________________
765d44e7 1129void AliCalorimeterUtils::InitParameters()
1130{
1131 //Initialize the parameters of the analysis.
7db7dcb6 1132
12e91266 1133 fEMCALGeoName = "";
1134 fPHOSGeoName = "";
7db7dcb6 1135
7cc7d3f8 1136 fEMCALGeoMatrixSet = kFALSE;
1137 fPHOSGeoMatrixSet = kFALSE;
7db7dcb6 1138
7cc7d3f8 1139 fRemoveBadChannels = kFALSE;
7db7dcb6 1140
7cc7d3f8 1141 fNCellsFromPHOSBorder = 0;
a5fb4114 1142
7db7dcb6 1143 fLocMaxCutE = 0.1 ;
1144 fLocMaxCutEDiff = 0.0 ;
1145
7cc7d3f8 1146 // fMaskCellColumns = new Int_t[fNMaskCellColumns];
1147 // fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
1148 // fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
1149 // fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
1150 // fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
1151 // fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
a5fb4114 1152
55d66f31 1153 fOADBSet = kFALSE;
1154 fOADBForEMCAL = kTRUE ;
1155 fOADBForPHOS = kFALSE;
1156
1157 fOADBFilePathEMCAL = "$ALICE_ROOT/OADB/EMCAL" ;
1158 fOADBFilePathPHOS = "$ALICE_ROOT/OADB/PHOS" ;
1159
765d44e7 1160}
1161
765d44e7 1162
cb5780f4 1163//_____________________________________________________
1164void AliCalorimeterUtils::InitPHOSBadChannelStatusMap()
1165{
765d44e7 1166 //Init PHOS bad channels map
1167 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
1168 //In order to avoid rewriting the same histograms
1169 Bool_t oldStatus = TH1::AddDirectoryStatus();
1170 TH1::AddDirectory(kFALSE);
7cc7d3f8 1171
78219bac 1172 fPHOSBadChannelMap = new TObjArray(5);
c5693f62 1173 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 1174
765d44e7 1175 fPHOSBadChannelMap->SetOwner(kTRUE);
1176 fPHOSBadChannelMap->Compress();
1177
1178 //In order to avoid rewriting the same histograms
1179 TH1::AddDirectory(oldStatus);
1180}
1181
cb5780f4 1182//______________________________________________________
1183void AliCalorimeterUtils::InitPHOSRecalibrationFactors()
1184{
09e819c9 1185 //Init EMCAL recalibration factors
1186 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
1187 //In order to avoid rewriting the same histograms
1188 Bool_t oldStatus = TH1::AddDirectoryStatus();
1189 TH1::AddDirectory(kFALSE);
7cc7d3f8 1190
78219bac 1191 fPHOSRecalibrationFactors = new TObjArray(5);
c5693f62 1192 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 1193 //Init the histograms with 1
1194 for (Int_t m = 0; m < 5; m++) {
1195 for (Int_t i = 0; i < 56; i++) {
1196 for (Int_t j = 0; j < 64; j++) {
c5693f62 1197 SetPHOSChannelRecalibrationFactor(m,j,i,1.);
09e819c9 1198 }
1199 }
1200 }
1201 fPHOSRecalibrationFactors->SetOwner(kTRUE);
1202 fPHOSRecalibrationFactors->Compress();
1203
1204 //In order to avoid rewriting the same histograms
1205 TH1::AddDirectory(oldStatus);
1206}
1207
1208
55d66f31 1209//__________________________________________________________
1210void AliCalorimeterUtils::InitEMCALGeometry(Int_t runnumber)
765d44e7 1211{
1212 //Initialize EMCAL geometry if it did not exist previously
55d66f31 1213
1214 if (!fEMCALGeo)
1215 {
1216 if(fEMCALGeoName=="")
1217 {
5408e59e 1218 if (runnumber < 140000 &&
1219 runnumber >= 100000) fEMCALGeoName = "EMCAL_FIRSTYEARV1";
1220 else if(runnumber >= 140000 &&
1221 runnumber < 171000) fEMCALGeoName = "EMCAL_COMPLETEV1";
1222 else fEMCALGeoName = "EMCAL_COMPLETE12SMV1";
55d66f31 1223 printf("AliCalorimeterUtils::InitEMCALGeometry() - Set EMCAL geometry name to <%s> for run %d\n",fEMCALGeoName.Data(),runnumber);
1224 }
1225
a38a48f2 1226 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
55d66f31 1227
1228 if(fDebug > 0)
1229 {
1230 printf("AliCalorimeterUtils::InitEMCALGeometry(run=%d)",runnumber);
765d44e7 1231 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1232 printf("\n");
1233 }
1234 }
1235}
1236
55d66f31 1237//_________________________________________________________
1238void AliCalorimeterUtils::InitPHOSGeometry(Int_t runnumber)
765d44e7 1239{
1240 //Initialize PHOS geometry if it did not exist previously
55d66f31 1241
1242 if (!fPHOSGeo)
1243 {
1244 if(fPHOSGeoName=="") fPHOSGeoName = "PHOSgeo";
1245
765d44e7 1246 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
55d66f31 1247
1248 if(fDebug > 0)
1249 {
1250 printf("AliCalorimeterUtils::InitPHOSGeometry(run=%d)",runnumber);
765d44e7 1251 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1252 printf("\n");
1253 }
1254 }
1255}
1256
55d66f31 1257//__________________________________________________________________
1258Bool_t AliCalorimeterUtils::MaskFrameCluster(const Int_t iSM,
1259 const Int_t ieta) const
1260{
1261 //Check if cell is in one of the regions where we have significant amount
1262 //of material in front. Only EMCAL
1263
1264 Int_t icol = ieta;
1265 if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
1266
1267 if (fNMaskCellColumns && fMaskCellColumns)
1268 {
1269 for (Int_t imask = 0; imask < fNMaskCellColumns; imask++)
1270 {
1271 if(icol==fMaskCellColumns[imask]) return kTRUE;
1272 }
1273 }
1274
1275 return kFALSE;
1276
1277}
1278
cb5780f4 1279//_________________________________________________________
765d44e7 1280void AliCalorimeterUtils::Print(const Option_t * opt) const
1281{
7cc7d3f8 1282
765d44e7 1283 //Print some relevant parameters set for the analysis
1284 if(! opt)
1285 return;
7cc7d3f8 1286
765d44e7 1287 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1288 printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
1289 printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
7cc7d3f8 1290 fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
247abff4 1291 if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
09e819c9 1292 printf("Recalibrate Clusters? %d\n",fRecalibration);
9584c261 1293 printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
1294 printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
f2ccb5b8 1295 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
7cc7d3f8 1296
7db7dcb6 1297 printf("Loc. Max. E > %2.2f\n", fLocMaxCutE);
1298 printf("Loc. Max. E Diff > %2.2f\n", fLocMaxCutEDiff);
1299
765d44e7 1300 printf(" \n") ;
1301}
1302
9369a2b1 1303//__________________________________________________________________________________________
7db7dcb6 1304void AliCalorimeterUtils::RecalibrateCellAmplitude(Float_t & amp,
9369a2b1 1305 const TString calo, const Int_t id) const
7db7dcb6 1306{
1307 //Recaculate cell energy if recalibration factor
1308
1309 Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
1310 Int_t nModule = GetModuleNumberCellIndexes(id,calo, icol, irow, iRCU);
1311
1312 if (IsRecalibrationOn())
1313 {
1314 if(calo == "PHOS")
1315 {
1316 amp *= GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
1317 }
1318 else
1319 {
1320 amp *= GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
1321 }
1322 }
1323}
1324
9369a2b1 1325//_________________________________________________________________________________
7db7dcb6 1326void AliCalorimeterUtils::RecalibrateCellTime(Double_t & time,
1327 const TString calo,
9369a2b1 1328 const Int_t id, const Int_t bc) const
7db7dcb6 1329{
1330 // Recalculate time if time recalibration available for EMCAL
1331 // not ready for PHOS
1332
1333 if(calo == "EMCAL" && GetEMCALRecoUtils()->IsTimeRecalibrationOn())
1334 {
1335 GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
1336 }
1337
1338}
1339
cb5780f4 1340//__________________________________________________________________________
1341Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster,
1342 AliVCaloCells * cells)
1343{
09e819c9 1344 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
7cc7d3f8 1345
898c9d44 1346 //Initialize some used variables
7db7dcb6 1347 Float_t frac = 0., energy = 0.;
898c9d44 1348
7db7dcb6 1349 if(cells)
1350 {
7cc7d3f8 1351 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
7db7dcb6 1352
7cc7d3f8 1353 UShort_t * index = cluster->GetCellsAbsId() ;
1354 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
7db7dcb6 1355
7cc7d3f8 1356 Int_t ncells = cluster->GetNCells();
7db7dcb6 1357
7cc7d3f8 1358 TString calo = "EMCAL";
7db7dcb6 1359 if(cluster->IsPHOS())
1360 calo = "PHOS";
7cc7d3f8 1361
1362 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1363 for(Int_t icell = 0; icell < ncells; icell++){
7db7dcb6 1364
1365 Int_t absId = index[icell];
1366
7cc7d3f8 1367 frac = fraction[icell];
1368 if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
7db7dcb6 1369
1370 Float_t amp = cells->GetCellAmplitude(absId);
1371 RecalibrateCellAmplitude(amp,calo, absId);
1372
7cc7d3f8 1373 if(fDebug>2)
7db7dcb6 1374 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, cell fraction %f, cell energy %f\n",
1375 calo.Data(),frac,cells->GetCellAmplitude(absId));
7cc7d3f8 1376
7db7dcb6 1377 energy += amp*frac;
7cc7d3f8 1378 }
1379
1380 if(fDebug>1)
1381 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
898c9d44 1382
1383 }// cells available
7db7dcb6 1384 else
1385 {
898c9d44 1386 Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
1387 }
1388
09e819c9 1389 return energy;
09e819c9 1390}
1391
cb5780f4 1392//__________________________________________________________________________________________
1393void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
1394{
9584c261 1395
1396 //Recalculate EMCAL cluster position
1397
19db8f8c 1398 fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
7cc7d3f8 1399
9584c261 1400}
cb5780f4 1401
55d66f31 1402//________________________________________________________________________________
1403void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event,
1404 TObjArray* clusterArray)
1405{
1406 //Recalculate track matching
1407
1408 if (fRecalculateMatching)
1409 {
1410 fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo) ;
1411 //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1412 //if(esdevent){
1413 // fEMCALRecoUtils->SetClusterMatchedToTrack(esdevent) ;
1414 // fEMCALRecoUtils->SetTracksMatchedToCluster(esdevent) ;
1415 //}
1416 }
1417}
1418
3c1d9afb 1419//___________________________________________________________________________
1420void AliCalorimeterUtils::SplitEnergy(const Int_t absId1, const Int_t absId2,
1421 AliVCluster* cluster,
1422 AliVCaloCells* cells,
1423 //Float_t & e1, Float_t & e2,
1424 AliAODCaloCluster* cluster1,
1425 AliAODCaloCluster* cluster2,
1426 const Int_t nMax,
1427 const Int_t eventNumber)
1428{
1429
1430 // Split energy of cluster between the 2 local maxima, sum energy on 3x3, and if the 2
1431 // maxima are too close and have common cells, split the energy between the 2
1432
1433 TH2F* hClusterMap = 0 ;
1434 TH2F* hClusterLocMax = 0 ;
1435 TH2F* hCluster1 = 0 ;
1436 TH2F* hCluster2 = 0 ;
1437
1438 if(fPlotCluster)
1439 {
1440 hClusterMap = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
1441 hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
1442 hCluster1 = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
1443 hCluster2 = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
1444 hClusterMap ->SetXTitle("column");
1445 hClusterMap ->SetYTitle("row");
1446 hClusterLocMax ->SetXTitle("column");
1447 hClusterLocMax ->SetYTitle("row");
1448 hCluster1 ->SetXTitle("column");
1449 hCluster1 ->SetYTitle("row");
1450 hCluster2 ->SetXTitle("column");
1451 hCluster2 ->SetYTitle("row");
1452 }
1453
1454 TString calorimeter = "EMCAL";
1455 if(cluster->IsPHOS())
1456 {
1457 calorimeter="PHOS";
1458 printf("AliCalorimeterUtils::SplitEnerg() Not supported for PHOS yet \n");
1459 return;
1460 }
1461
1462 const Int_t ncells = cluster->GetNCells();
1463 Int_t absIdList[ncells];
1464
1465 Float_t e1 = 0, e2 = 0 ;
1466 Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;
1467 Float_t eCluster = 0;
1468 Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1;
1469 for(Int_t iDigit = 0; iDigit < ncells; iDigit++ )
1470 {
1471 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
1472
1473 Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
1474 RecalibrateCellAmplitude(ec,calorimeter, absIdList[iDigit]);
1475 eCluster+=ec;
1476
1477 if(fPlotCluster)
1478 {
1479 //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
1480 sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1481 if(icol > maxCol) maxCol = icol;
1482 if(icol < minCol) minCol = icol;
1483 if(irow > maxRow) maxRow = irow;
1484 if(irow < minRow) minRow = irow;
1485 hClusterMap->Fill(icol,irow,ec);
1486 }
1487
1488 }
1489
1490 // Init counters and variables
1491 Int_t ncells1 = 1 ;
1492 UShort_t absIdList1[9] ;
1493 Double_t fracList1 [9] ;
1494 absIdList1[0] = absId1 ;
1495 fracList1 [0] = 1. ;
1496
1497 Float_t ecell1 = cells->GetCellAmplitude(absId1);
1498 RecalibrateCellAmplitude(ecell1, calorimeter, absId1);
1499 e1 = ecell1;
1500
1501 Int_t ncells2 = 1 ;
1502 UShort_t absIdList2[9] ;
1503 Double_t fracList2 [9] ;
1504 absIdList2[0] = absId2 ;
1505 fracList2 [0] = 1. ;
1506
1507 Float_t ecell2 = cells->GetCellAmplitude(absId2);
1508 RecalibrateCellAmplitude(ecell2, calorimeter, absId2);
1509 e2 = ecell2;
1510
1511 if(fPlotCluster)
1512 {
1513 Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
1514 sm = GetModuleNumberCellIndexes(absId1, calorimeter, icol1, irow1, iRCU) ;
1515 hClusterLocMax->Fill(icol1,irow1,ecell1);
1516 sm = GetModuleNumberCellIndexes(absId2, calorimeter, icol2, irow2, iRCU) ;
1517 hClusterLocMax->Fill(icol2,irow2,ecell2);
1518 }
1519
1520 // Very rough way to share the cluster energy
1521 Float_t eRemain = (eCluster-ecell1-ecell2)/2;
1522 Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
1523 Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
1524
1525 for(Int_t iDigit = 0; iDigit < ncells; iDigit++)
1526 {
1527 Int_t absId = absIdList[iDigit];
1528
1529 if(absId==absId1 || absId==absId2 || absId < 0) continue;
1530
1531 Float_t ecell = cells->GetCellAmplitude(absId);
1532 RecalibrateCellAmplitude(ecell, calorimeter, absId);
1533
1534 if(AreNeighbours(calorimeter, absId1,absId ))
1535 {
1536 absIdList1[ncells1]= absId;
1537
1538 if(AreNeighbours(calorimeter, absId2,absId ))
1539 {
1540 fracList1[ncells1] = shareFraction1;
1541 e1 += ecell*shareFraction1;
1542 }
1543 else
1544 {
1545 fracList1[ncells1] = 1.;
1546 e1 += ecell;
1547 }
1548
1549 ncells1++;
1550
1551 } // neigbour to cell1
1552
1553 if(AreNeighbours(calorimeter, absId2,absId ))
1554 {
1555 absIdList2[ncells2]= absId;
1556
1557 if(AreNeighbours(calorimeter, absId1,absId ))
1558 {
1559 fracList2[ncells2] = shareFraction2;
1560 e2 += ecell*shareFraction2;
1561 }
1562 else
1563 {
1564 fracList2[ncells2] = 1.;
1565 e2 += ecell;
1566 }
1567
1568 ncells2++;
1569
1570 } // neigbour to cell2
1571
1572 }
1573
1574 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",
1575 nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2);
1576
1577 cluster1->SetE(e1);
1578 cluster2->SetE(e2);
1579
1580 cluster1->SetNCells(ncells1);
1581 cluster2->SetNCells(ncells2);
1582
1583 cluster1->SetCellsAbsId(absIdList1);
1584 cluster2->SetCellsAbsId(absIdList2);
1585
1586 cluster1->SetCellsAmplitudeFraction(fracList1);
1587 cluster2->SetCellsAmplitudeFraction(fracList2);
1588
1589 if(calorimeter=="EMCAL")
1590 {
1591 GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
1592 GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
1593 }
1594
1595 if(fPlotCluster)
1596 {
1597 //printf("Cells of cluster1: ");
1598 for(Int_t iDigit = 0; iDigit < ncells1; iDigit++ )
1599 {
1600 //printf(" %d ",absIdList1[iDigit]);
1601
1602 sm = GetModuleNumberCellIndexes(absIdList1[iDigit], calorimeter, icol, irow, iRCU) ;
1603
1604 if( AreNeighbours(calorimeter, absId2,absIdList1[iDigit]) )
1605 hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit])*shareFraction1);
1606 else
1607 hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit]));
1608 }
1609
1610 //printf(" \n ");
1611 //printf("Cells of cluster2: ");
1612
1613 for(Int_t iDigit = 0; iDigit < ncells2; iDigit++ )
1614 {
1615 //printf(" %d ",absIdList2[iDigit]);
1616
1617 sm = GetModuleNumberCellIndexes(absIdList2[iDigit], calorimeter, icol, irow, iRCU) ;
1618 if( AreNeighbours(calorimeter, absId1,absIdList2[iDigit]) )
1619 hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit])*shareFraction2);
1620 else
1621 hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit]));
1622
1623 }
1624 //printf(" \n ");
1625
1626 gStyle->SetPadRightMargin(0.1);
1627 gStyle->SetPadLeftMargin(0.1);
1628 gStyle->SetOptStat(0);
1629 gStyle->SetOptFit(000000);
1630
1631 if(maxCol-minCol > maxRow-minRow)
1632 {
1633 maxRow+= maxCol-minCol;
1634 }
1635 else
1636 {
1637 maxCol+= maxRow-minRow;
1638 }
1639
1640 TCanvas * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
1641 c->Divide(2,2);
1642 c->cd(1);
1643 gPad->SetGridy();
1644 gPad->SetGridx();
1645 hClusterMap ->SetAxisRange(minCol, maxCol,"X");
1646 hClusterMap ->SetAxisRange(minRow, maxRow,"Y");
1647 hClusterMap ->Draw("colz");
1648 c->cd(2);
1649 gPad->SetGridy();
1650 gPad->SetGridx();
1651 hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
1652 hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
1653 hClusterLocMax ->Draw("colz");
1654 c->cd(3);
1655 gPad->SetGridy();
1656 gPad->SetGridx();
1657 hCluster1 ->SetAxisRange(minCol, maxCol,"X");
1658 hCluster1 ->SetAxisRange(minRow, maxRow,"Y");
1659 hCluster1 ->Draw("colz");
1660 c->cd(4);
1661 gPad->SetGridy();
1662 gPad->SetGridx();
1663 hCluster2 ->SetAxisRange(minCol, maxCol,"X");
1664 hCluster2 ->SetAxisRange(minRow, maxRow,"Y");
1665 hCluster2 ->Draw("colz");
1666
1667 if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",
1668 eventNumber,cluster->E(),nMax,ncells1,ncells2));
1669
1670 delete c;
1671 delete hClusterMap;
1672 delete hClusterLocMax;
1673 delete hCluster1;
1674 delete hCluster2;
1675 }
1676}
1677