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