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