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