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