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