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