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