Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[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),
5bddb161 64 fRecalibration(kFALSE), fRunDependentCorrection(kFALSE),
65 fPHOSRecalibrationFactors(), 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),
5bddb161 73 fOADBFilePathEMCAL(""), fOADBFilePathPHOS(""),
17a9cbf5 74 fImportGeometryFromFile(0), fImportGeometryFilePath("")
765d44e7 75{
76 //Ctor
77
78 //Initialize parameters
79 InitParameters();
90e32961 80 for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
81 for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix [i] = 0 ;
7cc7d3f8 82
765d44e7 83}
765d44e7 84
cb5780f4 85//_________________________________________
86AliCalorimeterUtils::~AliCalorimeterUtils()
87{
765d44e7 88 //Dtor
89
4df35693 90 //if(fPHOSGeo) delete fPHOSGeo ;
765d44e7 91 if(fEMCALGeo) delete fEMCALGeo ;
7cc7d3f8 92
765d44e7 93 if(fPHOSBadChannelMap) {
94 fPHOSBadChannelMap->Clear();
95 delete fPHOSBadChannelMap;
96 }
97
09e819c9 98 if(fPHOSRecalibrationFactors) {
78219bac 99 fPHOSRecalibrationFactors->Clear();
100 delete fPHOSRecalibrationFactors;
09e819c9 101 }
102
a5fb4114 103 if(fEMCALRecoUtils) delete fEMCALRecoUtils ;
104 if(fNMaskCellColumns) delete [] fMaskCellColumns;
9584c261 105
765d44e7 106}
107
55d66f31 108//____________________________________________________
109void AliCalorimeterUtils::AccessOADB(AliVEvent* event)
110{
111 // Set the AODB calibration, bad channels etc. parameters at least once
112 // alignment matrices from OADB done in SetGeometryMatrices
113
114 //Set it only once
115 if(fOADBSet) return ;
116
117 Int_t runnumber = event->GetRunNumber() ;
118 TString pass = GetPass();
119
120 // EMCAL
121 if(fOADBForEMCAL)
122 {
123 printf("AliCalorimeterUtils::SetOADBParameters() - Get AODB parameters from EMCAL in %s for run %d, and <%s> \n",fOADBFilePathEMCAL.Data(),runnumber,pass.Data());
124
125 Int_t nSM = fEMCALGeo->GetNumberOfSuperModules();
126
127 // Bad map
128 if(fRemoveBadChannels)
129 {
130 AliOADBContainer *contBC=new AliOADBContainer("");
131 contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePathEMCAL.Data()),"AliEMCALBadChannels");
132
133 TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber);
134
135 if(arrayBC)
136 {
6f557f7f 137 SwitchOnDistToBadChannelRecalculation();
138 printf("AliCalorimeterUtils::SetOADBParameters() - Remove EMCAL bad cells \n");
55d66f31 139
6f557f7f 140 for (Int_t i=0; i<nSM; ++i)
55d66f31 141 {
6f557f7f 142 TH2I *hbm = GetEMCALChannelStatusMap(i);
143
144 if (hbm)
145 delete hbm;
55d66f31 146
6f557f7f 147 hbm=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
148
149 if (!hbm)
55d66f31 150 {
6f557f7f 151 AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
152 continue;
153 }
154
155 hbm->SetDirectory(0);
156 SetEMCALChannelStatusMap(i,hbm);
157
158 } // loop
159 } else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT remove EMCAL bad channels\n"); // run array
55d66f31 160 } // Remove bad
161
162 // Energy Recalibration
163 if(fRecalibration)
164 {
165 AliOADBContainer *contRF=new AliOADBContainer("");
166
167 contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALRecalib");
168
169 TObjArray *recal=(TObjArray*)contRF->GetObject(runnumber);
170
171 if(recal)
172 {
173 TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
174
175 if(recalpass)
176 {
177 TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
178
179 if(recalib)
180 {
181 printf("AliCalorimeterUtils::SetOADBParameters() - Recalibrate EMCAL \n");
182 for (Int_t i=0; i<nSM; ++i)
183 {
184 TH2F *h = GetEMCALChannelRecalibrationFactors(i);
185
186 if (h)
187 delete h;
188
189 h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
190
191 if (!h)
192 {
193 AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
194 continue;
195 }
196
197 h->SetDirectory(0);
198
199 SetEMCALChannelRecalibrationFactors(i,h);
200 } // SM loop
6f557f7f 201 }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL, no params object array\n"); // array ok
202 }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for pass\n"); // array pass ok
203 }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for run\n"); // run number array ok
55d66f31 204
205 // once set, apply run dependent corrections if requested
7bf608c9 206 //fEMCALRecoUtils->SetRunDependentCorrections(runnumber);
207
55d66f31 208 } // Recalibration on
209
7bf608c9 210 // Energy Recalibration, apply on top of previous calibration factors
211 if(fRunDependentCorrection)
212 {
213 AliOADBContainer *contRFTD=new AliOADBContainer("");
214
215 contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALRunDepTempCalibCorrections");
216
217 TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber);
218
74d2310f 219 //If it did not exist for this run, get closes one
220 if (!htd)
221 {
222 AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runnumber));
223 // let's get the closest runnumber instead then..
224 Int_t lower = 0;
225 Int_t ic = 0;
226 Int_t maxEntry = contRFTD->GetNumberOfEntries();
227
228 while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < runnumber) ) {
229 lower = ic;
230 ic++;
231 }
232
233 Int_t closest = lower;
234 if ( (ic<maxEntry) &&
235 (contRFTD->LowerLimit(ic)-runnumber) < (runnumber - contRFTD->UpperLimit(lower)) ) {
236 closest = ic;
237 }
238
239 AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRFTD->LowerLimit(closest)));
240 htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
241 }
242
7bf608c9 243 if(htd)
244 {
245 printf("AliCalorimeterUtils::SetOADBParameters() - Recalibrate (Temperature) EMCAL \n");
246
247 for (Int_t ism=0; ism<nSM; ++ism)
248 {
249 for (Int_t icol=0; icol<48; ++icol)
250 {
251 for (Int_t irow=0; irow<24; ++irow)
252 {
253 Float_t factor = GetEMCALChannelRecalibrationFactor(ism,icol,irow);
254
255 Int_t absID = fEMCALGeo->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
256 factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
257 //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID,
258 // GetEMCALChannelRecalibrationFactor(ism,icol,irow) , htd->GetBinContent(absID) / 10000., factor);
259 SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
260 } // columns
261 } // rows
262 } // SM loop
263 }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL with T variations, no params TH1 \n");
264 } // Run by Run T calibration
55d66f31 265
266 // Time Recalibration
267 if(fEMCALRecoUtils->IsTimeRecalibrationOn())
268 {
269 AliOADBContainer *contTRF=new AliOADBContainer("");
270
f46af216 271 contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALTimeCalib");
55d66f31 272
273 TObjArray *trecal=(TObjArray*)contTRF->GetObject(runnumber);
274
275 if(trecal)
276 {
ff0dfe2e 277 TObjArray *trecalpass=(TObjArray*)trecal->FindObject(pass);
55d66f31 278
279 if(trecalpass)
280 {
f46af216 281 printf("AliCalorimeterUtils::SetOADBParameters() - Time Recalibrate EMCAL \n");
282 for (Int_t ibc = 0; ibc < 4; ++ibc)
55d66f31 283 {
f46af216 284 TH1F *h = GetEMCALChannelTimeRecalibrationFactors(ibc);
285
286 if (h)
287 delete h;
288
289 h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
290
291 if (!h)
55d66f31 292 {
f46af216 293 AliError(Form("Could not load hAllTimeAvBC%d",ibc));
294 continue;
295 }
296
297 h->SetDirectory(0);
298
299 SetEMCALChannelTimeRecalibrationFactors(ibc,h);
300 } // bunch crossing loop
6f557f7f 301 }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for pass\n"); // array pass ok
302 }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for run\n"); // run number array ok
55d66f31 303
304 } // Recalibration on
305
306 }// EMCAL
307
308 // PHOS
309 if(fOADBForPHOS)
310 {
311 printf("AliCalorimeterUtils::SetOADBParameters() - Get AODB parameters from PHOS in %s for run %d, and <%s> \n",fOADBFilePathPHOS.Data(),runnumber,pass.Data());
312
313 // Bad map
314 if(fRemoveBadChannels)
315 {
316 AliOADBContainer badmapContainer(Form("phosBadMap"));
317 TString fileName="$ALICE_ROOT/OADB/PHOS/PHOSBadMaps.root";
318 badmapContainer.InitFromFile(Form("%s/PHOSBadMaps.root",fOADBFilePathPHOS.Data()),"phosBadMap");
319
320 //Use a fixed run number from year 2010, this year not available yet.
321 TObjArray *maps = (TObjArray*)badmapContainer.GetObject(139000,"phosBadMap");
322 if(!maps)
323 {
324 printf("AliCalorimeterUtils::SetOADBParameters() - Can not read PHOS bad map for run %d.\n",runnumber) ;
325 }
326 else
327 {
328 printf("AliCalorimeterUtils::SetOADBParameters() - Setting PHOS bad map with name %s \n",maps->GetName()) ;
329 for(Int_t mod=1; mod<5;mod++)
330 {
331 TH2I *hbmPH = GetPHOSChannelStatusMap(mod);
332
333 if(hbmPH)
334 delete hbmPH ;
335
336 hbmPH = (TH2I*)maps->At(mod);
337
338 if(hbmPH) hbmPH->SetDirectory(0);
339
340 SetPHOSChannelStatusMap(mod-1,hbmPH);
341
342 } // modules loop
343 } // maps exist
344 } // Remove bad channels
345 } // PHOS
346
347 // Parameters already set once, so do not it again
348 fOADBSet = kTRUE;
349
350}
351
352//_____________________________________________________________
353void AliCalorimeterUtils::AccessGeometry(AliVEvent* inputEvent)
354{
355 //Set the calorimeters transformation matrices and init geometry
356
357 // First init the geometry, a priory not done before
358 Int_t runnumber = inputEvent->GetRunNumber() ;
359 InitPHOSGeometry (runnumber);
360 InitEMCALGeometry(runnumber);
361
362 //Get the EMCAL transformation geometry matrices from ESD
363 if(!fEMCALGeoMatrixSet && fEMCALGeo)
364 {
365 if(fLoadEMCALMatrices)
366 {
367 printf("AliCalorimeterUtils::AccessGeometry() - Load user defined EMCAL geometry matrices\n");
368
369 // OADB if available
370 AliOADBContainer emcGeoMat("AliEMCALgeo");
371 emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePathEMCAL.Data()),"AliEMCALgeo");
372 TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
373
374 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
375 {
376 if (!fEMCALMatrix[mod]) // Get it from OADB
377 {
378 if(fDebug > 1 )
379 printf("AliCalorimeterUtils::AccessGeometry() - EMCAL matrices SM %d, %p\n",
380 mod,((TGeoHMatrix*) matEMCAL->At(mod)));
381 //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
382
383 fEMCALMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
384 }
385
386 if(fEMCALMatrix[mod])
387 {
388 if(fDebug > 1)
389 fEMCALMatrix[mod]->Print();
390
391 fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;
392 }
393 }//SM loop
394
395 fEMCALGeoMatrixSet = kTRUE;//At least one, so good
396
397 }//Load matrices
398 else if (!gGeoManager)
399 {
400 if(fDebug > 1)
401 printf(" AliCalorimeterUtils::AccessGeometry() - Load EMCAL misalignment matrices. \n");
402 if(!strcmp(inputEvent->GetName(),"AliESDEvent"))
403 {
404 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
405 {
406 //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
407 if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod))
408 {
409 fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
410 }
411 }// loop over super modules
412
413 fEMCALGeoMatrixSet = kTRUE;//At least one, so good
414
415 }//ESD as input
416 else
417 {
418 if(fDebug > 1)
419 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
420 }//AOD as input
421 }//Get matrix from data
422 else if(gGeoManager)
423 {
424 fEMCALGeoMatrixSet = kTRUE;
425 }
426 }//EMCAL geo && no geoManager
427
428 //Get the PHOS transformation geometry matrices from ESD
429 if(!fPHOSGeoMatrixSet && fPHOSGeo)
430 {
431 if(fLoadPHOSMatrices)
432 {
433 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined PHOS geometry matrices\n");
434
435 // OADB if available
436 AliOADBContainer geomContainer("phosGeo");
437 geomContainer.InitFromFile(Form("%s/PHOSGeometry.root",fOADBFilePathPHOS.Data()),"PHOSRotationMatrixes");
438 TObjArray *matPHOS = (TObjArray*)geomContainer.GetObject(139000,"PHOSRotationMatrixes");
439
440 for(Int_t mod = 0 ; mod < 5 ; mod++)
441 {
442 if (!fPHOSMatrix[mod]) // Get it from OADB
443 {
444 if(fDebug > 1 )
445 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - PHOS matrices module %d, %p\n",
446 mod,((TGeoHMatrix*) matPHOS->At(mod)));
447 //((TGeoHMatrix*) matPHOS->At(mod))->Print();
448
449 fPHOSMatrix[mod] = (TGeoHMatrix*) matPHOS->At(mod) ;
450 }
451
452 // Set it, if it exists
453 if(fPHOSMatrix[mod])
454 {
455 if(fDebug > 1 )
456 fPHOSMatrix[mod]->Print();
457
458 fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;
459 }
460 }// SM loop
461
462 fPHOSGeoMatrixSet = kTRUE;//At least one, so good
463
464 }//Load matrices
465 else if (!gGeoManager)
466 {
467 if(fDebug > 1)
468 printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
469 if(!strcmp(inputEvent->GetName(),"AliESDEvent"))
470 {
471 for(Int_t mod = 0; mod < 5; mod++)
472 {
473 if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod))
474 {
475 //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
476 fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
477 }
478 }// loop over modules
479 fPHOSGeoMatrixSet = kTRUE; //At least one so good
480 }//ESD as input
481 else
482 {
483 if(fDebug > 1)
484 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
485 }//AOD as input
486 }// get matrix from data
487 else if(gGeoManager)
488 {
489 fPHOSGeoMatrixSet = kTRUE;
490 }
491 }//PHOS geo and geoManager was not set
492
493}
494
8a2dbbff 495//________________________________________________________________________________________
496Bool_t AliCalorimeterUtils::AreNeighbours(TString calo, Int_t absId1, Int_t absId2 ) const
7db7dcb6 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
8a2dbbff 684//__________________________________________________________________________________________________________
685Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter, UShort_t* cellList, Int_t nCells)
cb5780f4 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
8a2dbbff 789//___________________________________________________________________________________
790AliVTrack * AliCalorimeterUtils::GetMatchedTrack(AliVCluster* cluster,
791 AliVEvent* event, Int_t index) const
d832b695 792{
793 // Get the matched track given its index, usually just the first match
794 // Since it is different for ESDs and AODs here it is a wrap method to do it
795
796 AliVTrack *track = 0;
797
798 // EMCAL case only when matching is recalculated
799 if(cluster->IsEMCAL() && IsRecalculationOfClusterTrackMatchingOn())
800 {
801 Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(cluster->GetID());
802 //printf("track index %d, cluster ID %d \n ",trackIndex,cluster->GetID());
803
804 if(trackIndex < 0 )
805 {
806 printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d, from recalculation\n", trackIndex);
807 }
808 else
809 {
810 track = dynamic_cast<AliVTrack*> (event->GetTrack(trackIndex));
811 }
812
813 return track ;
814
815 }
816
817 // Normal case, get info from ESD or AOD
818 // ESDs
819 if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
820 {
821 Int_t iESDtrack = cluster->GetTrackMatchedIndex();
822
823 if(iESDtrack < 0 )
824 {
825 printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d\n", index);
826 return 0x0;
827 }
828
829 track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
830
831 }
832 else // AODs
833 {
834 if(cluster->GetNTracksMatched() > 0 )
835 track = dynamic_cast<AliVTrack*>(cluster->GetTrackMatched(index));
836 }
837
838 return track ;
839
840}
841
cb5780f4 842//_____________________________________________________________________________________________________
765d44e7 843Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
844{
845 //Get the EMCAL/PHOS module number that corresponds to this particle
846
847 Int_t absId = -1;
848 if(particle->GetDetector()=="EMCAL"){
849 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
850 if(fDebug > 2)
851 printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
7cc7d3f8 852 particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
765d44e7 853 return fEMCALGeo->GetSuperModuleNumber(absId) ;
854 }//EMCAL
55d66f31 855 else if(particle->GetDetector()=="PHOS")
856 {
46a3cde6 857 // In case we use the MC reader, the input are TParticles,
858 // in this case use the corresponing method in PHOS Geometry to get the particle.
859 if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
860 {
861 Int_t mod =-1;
862 Double_t z = 0., x=0.;
863 TParticle* primary = 0x0;
864 AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
865 if(stack) {
866 primary = stack->Particle(particle->GetCaloLabel(0));
867 }
868 else {
280e6766 869 Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
46a3cde6 870 }
898c9d44 871
872 if(primary){
873 fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
874 }
875 else{
876 Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
877 }
46a3cde6 878 return mod;
879 }
880 // Input are ESDs or AODs, get the PHOS module number like this.
881 else{
2206e918 882 //FIXME
883 //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
884 //return GetModuleNumber(cluster);
885 //MEFIX
886 return -1;
46a3cde6 887 }
765d44e7 888 }//PHOS
889
890 return -1;
891}
892
cb5780f4 893//_____________________________________________________________________
c8fe2783 894Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
765d44e7 895{
280e6766 896 //Get the EMCAL/PHOS module number that corresponds to this cluster
765d44e7 897 TLorentzVector lv;
898 Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
55d66f31 899 if(!cluster)
900 {
9cbbc28b 901 if(fDebug > 1) printf("AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
902 return -1;
903 }
55d66f31 904
765d44e7 905 cluster->GetMomentum(lv,v);
906 Float_t phi = lv.Phi();
907 if(phi < 0) phi+=TMath::TwoPi();
908 Int_t absId = -1;
c8fe2783 909 if(cluster->IsEMCAL()){
765d44e7 910 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
911 if(fDebug > 2)
280e6766 912 printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
7cc7d3f8 913 lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
765d44e7 914 return fEMCALGeo->GetSuperModuleNumber(absId) ;
915 }//EMCAL
55d66f31 916 else if(cluster->IsPHOS())
917 {
765d44e7 918 Int_t relId[4];
55d66f31 919 if ( cluster->GetNCells() > 0)
920 {
765d44e7 921 absId = cluster->GetCellAbsId(0);
922 if(fDebug > 2)
280e6766 923 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
7cc7d3f8 924 lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
765d44e7 925 }
926 else return -1;
927
55d66f31 928 if ( absId >= 0)
929 {
765d44e7 930 fPHOSGeo->AbsToRelNumbering(absId,relId);
931 if(fDebug > 2)
280e6766 932 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
765d44e7 933 return relId[0]-1;
934 }
935 else return -1;
936 }//PHOS
937
938 return -1;
939}
940
cb5780f4 941//___________________________________________________________________________________________________
8a2dbbff 942Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(Int_t absId, TString calo,
cb5780f4 943 Int_t & icol, Int_t & irow, Int_t & iRCU) const
765d44e7 944{
945 //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
946 Int_t imod = -1;
55d66f31 947 if ( absId >= 0)
948 {
949 if(calo=="EMCAL")
950 {
765d44e7 951 Int_t iTower = -1, iIphi = -1, iIeta = -1;
952 fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
953 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
55d66f31 954 if(imod < 0 || irow < 0 || icol < 0 )
955 {
280e6766 956 Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
957 }
958
765d44e7 959 //RCU0
55d66f31 960 if(imod < 10 )
961 {
962 if (0<=irow&&irow<8) iRCU=0; // first cable row
963 else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
964 //second cable row
965 //RCU1
966 else if (8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
967 //second cable row
968 else if (16<=irow&&irow<24) iRCU=1; // third cable row
969
970 if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
971 }
972 else
973 {
974 // Last 2 SM have one single SRU, just assign RCU 0
975 iRCU = 0 ;
976 }
977
978 if (iRCU<0)
979 {
280e6766 980 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
765d44e7 981 }
982
983 return imod ;
55d66f31 984
765d44e7 985 }//EMCAL
55d66f31 986 else //PHOS
987 {
765d44e7 988 Int_t relId[4];
989 fPHOSGeo->AbsToRelNumbering(absId,relId);
990 irow = relId[2];
991 icol = relId[3];
992 imod = relId[0]-1;
993 iRCU= (Int_t)(relId[2]-1)/16 ;
994 //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
55d66f31 995 if (iRCU >= 4)
996 {
280e6766 997 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
765d44e7 998 }
999 return imod;
1000 }//PHOS
1001 }
1002
1003 return -1;
1004}
1005
7db7dcb6 1006//___________________________________________________________________________________________
71e3889f 1007Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells)
1008{
1009 // Find local maxima in cluster
1010
1011 const Int_t nc = cluster->GetNCells();
2bf17171 1012 Int_t absIdList[nc];
1013 Float_t maxEList[nc];
71e3889f 1014
1015 Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList);
1016
71e3889f 1017 return nMax;
1018
1019}
1020
1021//___________________________________________________________________________________________
7db7dcb6 1022Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
1023 Int_t *absIdList, Float_t *maxEList)
1024{
1025 // Find local maxima in cluster
f241ecc3 1026
7db7dcb6 1027 Int_t iDigitN = 0 ;
1028 Int_t iDigit = 0 ;
1029 Int_t absId1 = -1 ;
1030 Int_t absId2 = -1 ;
1031 const Int_t nCells = cluster->GetNCells();
1032
1033 TString calorimeter = "EMCAL";
1034 if(!cluster->IsEMCAL()) calorimeter = "PHOS";
1035
1036 //printf("cluster : ncells %d \n",nCells);
f241ecc3 1037
1038 Float_t emax = 0;
1039 Int_t idmax =-1;
9369a2b1 1040 for(iDigit = 0; iDigit < nCells ; iDigit++)
1041 {
7db7dcb6 1042 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit] ;
f241ecc3 1043 Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1044 RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1045 if( en > emax )
1046 {
1047 emax = en ;
1048 idmax = absIdList[iDigit] ;
1049 }
1050 //Int_t icol = -1, irow = -1, iRCU = -1;
1051 //Int_t sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1052 //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
7db7dcb6 1053 }
1054
9369a2b1 1055 for(iDigit = 0 ; iDigit < nCells; iDigit++)
1056 {
1057 if(absIdList[iDigit]>=0)
1058 {
1059 absId1 = cluster->GetCellsAbsId()[iDigit];
7db7dcb6 1060
1061 Float_t en1 = cells->GetCellAmplitude(absId1);
1062 RecalibrateCellAmplitude(en1,calorimeter,absId1);
f241ecc3 1063
9369a2b1 1064 //printf("%d : absIDi %d, E %f\n",iDigit, absId1,en1);
7db7dcb6 1065
9369a2b1 1066 for(iDigitN = 0; iDigitN < nCells; iDigitN++)
1067 {
1068 absId2 = cluster->GetCellsAbsId()[iDigitN] ;
7db7dcb6 1069
9369a2b1 1070 if(absId2==-1 || absId2==absId1) continue;
7db7dcb6 1071
9369a2b1 1072 //printf("\t %d : absIDj %d\n",iDigitN, absId2);
7db7dcb6 1073
1074 Float_t en2 = cells->GetCellAmplitude(absId2);
1075 RecalibrateCellAmplitude(en2,calorimeter,absId2);
f241ecc3 1076
9369a2b1 1077 //printf("\t %d : absIDj %d, E %f\n",iDigitN, absId2,en2);
7db7dcb6 1078
9369a2b1 1079 if ( AreNeighbours(calorimeter, absId1, absId2) )
1080 {
1081 // printf("\t \t Neighbours \n");
1082 if (en1 > en2 )
1083 {
7db7dcb6 1084 absIdList[iDigitN] = -1 ;
1085 //printf("\t \t indexN %d not local max\n",iDigitN);
1086 // but may be digit too is not local max ?
1087 if(en1 < en2 + fLocMaxCutEDiff) {
9369a2b1 1088 //printf("\t \t index %d not local max cause locMaxCutEDiff\n",iDigit);
7db7dcb6 1089 absIdList[iDigit] = -1 ;
1090 }
1091 }
9369a2b1 1092 else
1093 {
7db7dcb6 1094 absIdList[iDigit] = -1 ;
1095 //printf("\t \t index %d not local max\n",iDigitN);
1096 // but may be digitN too is not local max ?
1097 if(en1 > en2 - fLocMaxCutEDiff)
1098 {
1099 absIdList[iDigitN] = -1 ;
9369a2b1 1100 //printf("\t \t indexN %d not local max cause locMaxCutEDiff\n",iDigit);
7db7dcb6 1101 }
1102 }
9369a2b1 1103 } // if Are neighbours
1104 //else printf("\t \t NOT Neighbours \n");
7db7dcb6 1105 } // while digitN
1106 } // slot not empty
1107 } // while digit
1108
1109 iDigitN = 0 ;
9369a2b1 1110 for(iDigit = 0; iDigit < nCells; iDigit++)
1111 {
1112 if(absIdList[iDigit]>=0 )
1113 {
7db7dcb6 1114 absIdList[iDigitN] = absIdList[iDigit] ;
1115 Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1116 RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1117 if(en < fLocMaxCutE) continue; // Maxima only with seed energy at least
1118 maxEList[iDigitN] = en ;
1119 //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
1120 iDigitN++ ;
1121 }
1122 }
1123
f241ecc3 1124 if(iDigitN == 0)
1125 {
1126 if(fDebug > 0)
1127 printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - No local maxima found, assign highest energy cell as maxima, id %d, en cell %2.2f, en cluster %2.2f\n",
1128 idmax,emax,cluster->E());
1129 iDigitN = 1 ;
1130 maxEList[0] = emax ;
1131 absIdList[0] = idmax ;
1132 }
1133
1134 if(fDebug > 0)
1135 {
1136 printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - In cluster E %2.2f, M02 %2.2f, M20 %2.2f, N maxima %d \n",
1137 cluster->E(),cluster->GetM02(),cluster->GetM20(), iDigitN);
1138
1139 if(fDebug > 1) for(Int_t imax = 0; imax < iDigitN; imax++)
1140 {
1141 printf(" \t i %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
1142 }
1143 }
7db7dcb6 1144
1145 return iDigitN ;
1146
1147}
1148
55d66f31 1149//____________________________________
1150TString AliCalorimeterUtils::GetPass()
1151{
1152 // Get passx from filename.
1153
1154 if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1155 {
1156 AliError("AliCalorimeterUtils::GetPass() - Pointer to tree = 0, returning null\n");
1157 return TString("");
1158 }
1159
1160 if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1161 {
1162 AliError("AliCalorimeterUtils::GetPass() - Null pointer input file, returning null\n");
1163 return TString("");
1164 }
1165
1166 TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1167 if (pass.Contains("ass1")) return TString("pass1");
1168 else if (pass.Contains("ass2")) return TString("pass2");
1169 else if (pass.Contains("ass3")) return TString("pass3");
4898bb87 1170 else if (pass.Contains("ass4")) return TString("pass4");
1171 else if (pass.Contains("ass5")) return TString("pass5");
fca053d8 1172 else if (pass.Contains("calo") || pass.Contains("high_lumi"))
1173 {
1174 printf("AliCalorimeterUtils::GetPass() - Path contains <calo> or <high-lumi>, set as <pass1>\n");
1175 return TString("pass1");
1176 }
55d66f31 1177
1178 // No condition fullfilled, give a default value
1179 printf("AliCalorimeterUtils::GetPass() - Pass number string not found \n");
1180 return TString("");
1181
1182}
1183
cb5780f4 1184//________________________________________
765d44e7 1185void AliCalorimeterUtils::InitParameters()
1186{
1187 //Initialize the parameters of the analysis.
7db7dcb6 1188
12e91266 1189 fEMCALGeoName = "";
1190 fPHOSGeoName = "";
7db7dcb6 1191
7cc7d3f8 1192 fEMCALGeoMatrixSet = kFALSE;
1193 fPHOSGeoMatrixSet = kFALSE;
7db7dcb6 1194
7cc7d3f8 1195 fRemoveBadChannels = kFALSE;
7db7dcb6 1196
7cc7d3f8 1197 fNCellsFromPHOSBorder = 0;
a5fb4114 1198
7db7dcb6 1199 fLocMaxCutE = 0.1 ;
1200 fLocMaxCutEDiff = 0.0 ;
1201
7cc7d3f8 1202 // fMaskCellColumns = new Int_t[fNMaskCellColumns];
1203 // fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
1204 // fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
1205 // fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
1206 // fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
1207 // fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
a5fb4114 1208
55d66f31 1209 fOADBSet = kFALSE;
1210 fOADBForEMCAL = kTRUE ;
1211 fOADBForPHOS = kFALSE;
1212
1213 fOADBFilePathEMCAL = "$ALICE_ROOT/OADB/EMCAL" ;
1214 fOADBFilePathPHOS = "$ALICE_ROOT/OADB/PHOS" ;
1215
9e536695 1216 fImportGeometryFromFile = kTRUE;
1217 fImportGeometryFilePath = "";
1218
765d44e7 1219}
1220
765d44e7 1221
cb5780f4 1222//_____________________________________________________
1223void AliCalorimeterUtils::InitPHOSBadChannelStatusMap()
1224{
765d44e7 1225 //Init PHOS bad channels map
1226 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
1227 //In order to avoid rewriting the same histograms
1228 Bool_t oldStatus = TH1::AddDirectoryStatus();
1229 TH1::AddDirectory(kFALSE);
7cc7d3f8 1230
78219bac 1231 fPHOSBadChannelMap = new TObjArray(5);
c5693f62 1232 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 1233
765d44e7 1234 fPHOSBadChannelMap->SetOwner(kTRUE);
1235 fPHOSBadChannelMap->Compress();
1236
1237 //In order to avoid rewriting the same histograms
1238 TH1::AddDirectory(oldStatus);
1239}
1240
cb5780f4 1241//______________________________________________________
1242void AliCalorimeterUtils::InitPHOSRecalibrationFactors()
1243{
09e819c9 1244 //Init EMCAL recalibration factors
1245 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
1246 //In order to avoid rewriting the same histograms
1247 Bool_t oldStatus = TH1::AddDirectoryStatus();
1248 TH1::AddDirectory(kFALSE);
7cc7d3f8 1249
78219bac 1250 fPHOSRecalibrationFactors = new TObjArray(5);
c5693f62 1251 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 1252 //Init the histograms with 1
1253 for (Int_t m = 0; m < 5; m++) {
1254 for (Int_t i = 0; i < 56; i++) {
1255 for (Int_t j = 0; j < 64; j++) {
c5693f62 1256 SetPHOSChannelRecalibrationFactor(m,j,i,1.);
09e819c9 1257 }
1258 }
1259 }
1260 fPHOSRecalibrationFactors->SetOwner(kTRUE);
1261 fPHOSRecalibrationFactors->Compress();
1262
1263 //In order to avoid rewriting the same histograms
1264 TH1::AddDirectory(oldStatus);
1265}
1266
1267
55d66f31 1268//__________________________________________________________
1269void AliCalorimeterUtils::InitEMCALGeometry(Int_t runnumber)
765d44e7 1270{
1271 //Initialize EMCAL geometry if it did not exist previously
55d66f31 1272
1273 if (!fEMCALGeo)
1274 {
1275 if(fEMCALGeoName=="")
1276 {
5408e59e 1277 if (runnumber < 140000 &&
1278 runnumber >= 100000) fEMCALGeoName = "EMCAL_FIRSTYEARV1";
1279 else if(runnumber >= 140000 &&
1280 runnumber < 171000) fEMCALGeoName = "EMCAL_COMPLETEV1";
1281 else fEMCALGeoName = "EMCAL_COMPLETE12SMV1";
55d66f31 1282 printf("AliCalorimeterUtils::InitEMCALGeometry() - Set EMCAL geometry name to <%s> for run %d\n",fEMCALGeoName.Data(),runnumber);
1283 }
1284
a38a48f2 1285 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
55d66f31 1286
9e536695 1287 // Init geometry, I do not like much to do it like this ...
1288 if(fImportGeometryFromFile && !gGeoManager)
1289 {
1290 if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
1291 {
1292 // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1293 if (runnumber < 140000 &&
1294 runnumber >= 100000) fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2010.root";
1295 if (runnumber >= 140000 &&
1296 runnumber < 171000) fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2011.root";
1297 else fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2012.root"; // 2012-2013
1298
1299 }
1300 printf("AliCalorimeterUtils::InitEMCALGeometry() - Import %s\n",fImportGeometryFilePath.Data());
1301 TGeoManager::Import(fImportGeometryFilePath) ; // default need file "geometry.root" in local dir!!!!
1302 }
1303
1304
55d66f31 1305 if(fDebug > 0)
1306 {
1307 printf("AliCalorimeterUtils::InitEMCALGeometry(run=%d)",runnumber);
765d44e7 1308 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1309 printf("\n");
1310 }
1311 }
1312}
1313
55d66f31 1314//_________________________________________________________
1315void AliCalorimeterUtils::InitPHOSGeometry(Int_t runnumber)
765d44e7 1316{
1317 //Initialize PHOS geometry if it did not exist previously
55d66f31 1318
1319 if (!fPHOSGeo)
1320 {
1321 if(fPHOSGeoName=="") fPHOSGeoName = "PHOSgeo";
1322
765d44e7 1323 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
55d66f31 1324
1325 if(fDebug > 0)
1326 {
1327 printf("AliCalorimeterUtils::InitPHOSGeometry(run=%d)",runnumber);
765d44e7 1328 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1329 printf("\n");
1330 }
1331 }
1332}
1333
8a2dbbff 1334//_______________________________________________________________________
1335Bool_t AliCalorimeterUtils::MaskFrameCluster(Int_t iSM, Int_t ieta) const
55d66f31 1336{
1337 //Check if cell is in one of the regions where we have significant amount
1338 //of material in front. Only EMCAL
1339
1340 Int_t icol = ieta;
1341 if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
1342
1343 if (fNMaskCellColumns && fMaskCellColumns)
1344 {
1345 for (Int_t imask = 0; imask < fNMaskCellColumns; imask++)
1346 {
1347 if(icol==fMaskCellColumns[imask]) return kTRUE;
1348 }
1349 }
1350
1351 return kFALSE;
1352
1353}
1354
cb5780f4 1355//_________________________________________________________
765d44e7 1356void AliCalorimeterUtils::Print(const Option_t * opt) const
1357{
7cc7d3f8 1358
765d44e7 1359 //Print some relevant parameters set for the analysis
1360 if(! opt)
1361 return;
7cc7d3f8 1362
765d44e7 1363 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1364 printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
1365 printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
7cc7d3f8 1366 fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
247abff4 1367 if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
7bf608c9 1368 printf("Recalibrate Clusters? %d, run by run %d\n",fRecalibration,fRunDependentCorrection);
9584c261 1369 printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
1370 printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
f2ccb5b8 1371 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
7cc7d3f8 1372
7db7dcb6 1373 printf("Loc. Max. E > %2.2f\n", fLocMaxCutE);
1374 printf("Loc. Max. E Diff > %2.2f\n", fLocMaxCutEDiff);
1375
765d44e7 1376 printf(" \n") ;
1377}
1378
8a2dbbff 1379//_____________________________________________________________________________________________
1380void AliCalorimeterUtils::RecalibrateCellAmplitude(Float_t & amp, TString calo, Int_t id) const
7db7dcb6 1381{
1382 //Recaculate cell energy if recalibration factor
1383
1384 Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
1385 Int_t nModule = GetModuleNumberCellIndexes(id,calo, icol, irow, iRCU);
1386
1387 if (IsRecalibrationOn())
1388 {
1389 if(calo == "PHOS")
1390 {
1391 amp *= GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
1392 }
1393 else
1394 {
1395 amp *= GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
1396 }
1397 }
1398}
1399
8a2dbbff 1400//____________________________________________________________________________________________________
1401void AliCalorimeterUtils::RecalibrateCellTime(Double_t & time, TString calo, Int_t id, Int_t bc) const
7db7dcb6 1402{
1403 // Recalculate time if time recalibration available for EMCAL
1404 // not ready for PHOS
1405
1406 if(calo == "EMCAL" && GetEMCALRecoUtils()->IsTimeRecalibrationOn())
1407 {
1408 GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
1409 }
1410
1411}
1412
cb5780f4 1413//__________________________________________________________________________
1414Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster,
1415 AliVCaloCells * cells)
1416{
09e819c9 1417 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
7cc7d3f8 1418
898c9d44 1419 //Initialize some used variables
7db7dcb6 1420 Float_t frac = 0., energy = 0.;
898c9d44 1421
7db7dcb6 1422 if(cells)
1423 {
7cc7d3f8 1424 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
7db7dcb6 1425
7cc7d3f8 1426 UShort_t * index = cluster->GetCellsAbsId() ;
1427 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
7db7dcb6 1428
7cc7d3f8 1429 Int_t ncells = cluster->GetNCells();
7db7dcb6 1430
7cc7d3f8 1431 TString calo = "EMCAL";
7db7dcb6 1432 if(cluster->IsPHOS())
1433 calo = "PHOS";
7cc7d3f8 1434
1435 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1436 for(Int_t icell = 0; icell < ncells; icell++){
7db7dcb6 1437
1438 Int_t absId = index[icell];
1439
7cc7d3f8 1440 frac = fraction[icell];
1441 if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
7db7dcb6 1442
1443 Float_t amp = cells->GetCellAmplitude(absId);
1444 RecalibrateCellAmplitude(amp,calo, absId);
1445
7cc7d3f8 1446 if(fDebug>2)
7db7dcb6 1447 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, cell fraction %f, cell energy %f\n",
1448 calo.Data(),frac,cells->GetCellAmplitude(absId));
7cc7d3f8 1449
7db7dcb6 1450 energy += amp*frac;
7cc7d3f8 1451 }
1452
1453 if(fDebug>1)
1454 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
898c9d44 1455
1456 }// cells available
7db7dcb6 1457 else
1458 {
898c9d44 1459 Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
1460 }
1461
09e819c9 1462 return energy;
09e819c9 1463}
1464
cb5780f4 1465//__________________________________________________________________________________________
1466void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
1467{
9584c261 1468
1469 //Recalculate EMCAL cluster position
1470
19db8f8c 1471 fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
7cc7d3f8 1472
9584c261 1473}
cb5780f4 1474
55d66f31 1475//________________________________________________________________________________
1476void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event,
1477 TObjArray* clusterArray)
1478{
1479 //Recalculate track matching
1480
1481 if (fRecalculateMatching)
1482 {
1483 fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo) ;
1484 //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1485 //if(esdevent){
1486 // fEMCALRecoUtils->SetClusterMatchedToTrack(esdevent) ;
1487 // fEMCALRecoUtils->SetTracksMatchedToCluster(esdevent) ;
1488 //}
1489 }
1490}
1491
3c1d9afb 1492//___________________________________________________________________________
8a2dbbff 1493void AliCalorimeterUtils::SplitEnergy(Int_t absId1, Int_t absId2,
1494 AliVCluster* cluster,
3c1d9afb 1495 AliVCaloCells* cells,
1496 //Float_t & e1, Float_t & e2,
1497 AliAODCaloCluster* cluster1,
1498 AliAODCaloCluster* cluster2,
8a2dbbff 1499 Int_t nMax, Int_t eventNumber)
3c1d9afb 1500{
1501
1502 // Split energy of cluster between the 2 local maxima, sum energy on 3x3, and if the 2
1503 // maxima are too close and have common cells, split the energy between the 2
1504
1505 TH2F* hClusterMap = 0 ;
1506 TH2F* hClusterLocMax = 0 ;
1507 TH2F* hCluster1 = 0 ;
1508 TH2F* hCluster2 = 0 ;
1509
1510 if(fPlotCluster)
1511 {
1512 hClusterMap = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
1513 hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
1514 hCluster1 = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
1515 hCluster2 = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
1516 hClusterMap ->SetXTitle("column");
1517 hClusterMap ->SetYTitle("row");
1518 hClusterLocMax ->SetXTitle("column");
1519 hClusterLocMax ->SetYTitle("row");
1520 hCluster1 ->SetXTitle("column");
1521 hCluster1 ->SetYTitle("row");
1522 hCluster2 ->SetXTitle("column");
1523 hCluster2 ->SetYTitle("row");
1524 }
1525
1526 TString calorimeter = "EMCAL";
1527 if(cluster->IsPHOS())
1528 {
1529 calorimeter="PHOS";
1530 printf("AliCalorimeterUtils::SplitEnerg() Not supported for PHOS yet \n");
1531 return;
1532 }
1533
1534 const Int_t ncells = cluster->GetNCells();
1535 Int_t absIdList[ncells];
1536
1537 Float_t e1 = 0, e2 = 0 ;
1538 Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;
1539 Float_t eCluster = 0;
1540 Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1;
1541 for(Int_t iDigit = 0; iDigit < ncells; iDigit++ )
1542 {
1543 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
1544
1545 Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
1546 RecalibrateCellAmplitude(ec,calorimeter, absIdList[iDigit]);
1547 eCluster+=ec;
1548
1549 if(fPlotCluster)
1550 {
1551 //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
1552 sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1553 if(icol > maxCol) maxCol = icol;
1554 if(icol < minCol) minCol = icol;
1555 if(irow > maxRow) maxRow = irow;
1556 if(irow < minRow) minRow = irow;
1557 hClusterMap->Fill(icol,irow,ec);
1558 }
1559
1560 }
1561
1562 // Init counters and variables
1563 Int_t ncells1 = 1 ;
1564 UShort_t absIdList1[9] ;
1565 Double_t fracList1 [9] ;
1566 absIdList1[0] = absId1 ;
1567 fracList1 [0] = 1. ;
1568
1569 Float_t ecell1 = cells->GetCellAmplitude(absId1);
1570 RecalibrateCellAmplitude(ecell1, calorimeter, absId1);
1571 e1 = ecell1;
1572
1573 Int_t ncells2 = 1 ;
1574 UShort_t absIdList2[9] ;
1575 Double_t fracList2 [9] ;
1576 absIdList2[0] = absId2 ;
1577 fracList2 [0] = 1. ;
1578
1579 Float_t ecell2 = cells->GetCellAmplitude(absId2);
1580 RecalibrateCellAmplitude(ecell2, calorimeter, absId2);
1581 e2 = ecell2;
1582
1583 if(fPlotCluster)
1584 {
1585 Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
1586 sm = GetModuleNumberCellIndexes(absId1, calorimeter, icol1, irow1, iRCU) ;
1587 hClusterLocMax->Fill(icol1,irow1,ecell1);
1588 sm = GetModuleNumberCellIndexes(absId2, calorimeter, icol2, irow2, iRCU) ;
1589 hClusterLocMax->Fill(icol2,irow2,ecell2);
1590 }
1591
1592 // Very rough way to share the cluster energy
1593 Float_t eRemain = (eCluster-ecell1-ecell2)/2;
1594 Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
1595 Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
1596
1597 for(Int_t iDigit = 0; iDigit < ncells; iDigit++)
1598 {
1599 Int_t absId = absIdList[iDigit];
1600
1601 if(absId==absId1 || absId==absId2 || absId < 0) continue;
1602
1603 Float_t ecell = cells->GetCellAmplitude(absId);
1604 RecalibrateCellAmplitude(ecell, calorimeter, absId);
1605
1606 if(AreNeighbours(calorimeter, absId1,absId ))
1607 {
1608 absIdList1[ncells1]= absId;
1609
1610 if(AreNeighbours(calorimeter, absId2,absId ))
1611 {
1612 fracList1[ncells1] = shareFraction1;
1613 e1 += ecell*shareFraction1;
1614 }
1615 else
1616 {
1617 fracList1[ncells1] = 1.;
1618 e1 += ecell;
1619 }
1620
1621 ncells1++;
1622
1623 } // neigbour to cell1
1624
1625 if(AreNeighbours(calorimeter, absId2,absId ))
1626 {
1627 absIdList2[ncells2]= absId;
1628
1629 if(AreNeighbours(calorimeter, absId1,absId ))
1630 {
1631 fracList2[ncells2] = shareFraction2;
1632 e2 += ecell*shareFraction2;
1633 }
1634 else
1635 {
1636 fracList2[ncells2] = 1.;
1637 e2 += ecell;
1638 }
1639
1640 ncells2++;
1641
1642 } // neigbour to cell2
1643
1644 }
1645
715d1872 1646 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 1647 nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2);
1648
1649 cluster1->SetE(e1);
1650 cluster2->SetE(e2);
1651
1652 cluster1->SetNCells(ncells1);
1653 cluster2->SetNCells(ncells2);
1654
1655 cluster1->SetCellsAbsId(absIdList1);
1656 cluster2->SetCellsAbsId(absIdList2);
1657
1658 cluster1->SetCellsAmplitudeFraction(fracList1);
1659 cluster2->SetCellsAmplitudeFraction(fracList2);
1660
715d1872 1661 //Correct linearity
1662 CorrectClusterEnergy(cluster1) ;
1663 CorrectClusterEnergy(cluster2) ;
1664
3c1d9afb 1665 if(calorimeter=="EMCAL")
1666 {
1667 GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
1668 GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
1669 }
1670
1671 if(fPlotCluster)
1672 {
1673 //printf("Cells of cluster1: ");
1674 for(Int_t iDigit = 0; iDigit < ncells1; iDigit++ )
1675 {
1676 //printf(" %d ",absIdList1[iDigit]);
1677
1678 sm = GetModuleNumberCellIndexes(absIdList1[iDigit], calorimeter, icol, irow, iRCU) ;
1679
8cf5a81b 1680 Float_t ecell = cells->GetCellAmplitude(absIdList1[iDigit]);
1681 RecalibrateCellAmplitude(ecell, calorimeter, absIdList1[iDigit]);
1682
1683 if( AreNeighbours(calorimeter, absId2,absIdList1[iDigit]) && absId1!=absIdList1[iDigit])
1684 hCluster1->Fill(icol,irow,ecell*shareFraction1);
3c1d9afb 1685 else
8cf5a81b 1686 hCluster1->Fill(icol,irow,ecell);
3c1d9afb 1687 }
1688
1689 //printf(" \n ");
1690 //printf("Cells of cluster2: ");
1691
1692 for(Int_t iDigit = 0; iDigit < ncells2; iDigit++ )
1693 {
1694 //printf(" %d ",absIdList2[iDigit]);
1695
1696 sm = GetModuleNumberCellIndexes(absIdList2[iDigit], calorimeter, icol, irow, iRCU) ;
8cf5a81b 1697
1698 Float_t ecell = cells->GetCellAmplitude(absIdList2[iDigit]);
1699 RecalibrateCellAmplitude(ecell, calorimeter, absIdList2[iDigit]);
1700
1701 if( AreNeighbours(calorimeter, absId1,absIdList2[iDigit]) && absId2!=absIdList2[iDigit])
1702 hCluster2->Fill(icol,irow,ecell*shareFraction2);
3c1d9afb 1703 else
8cf5a81b 1704 hCluster2->Fill(icol,irow,ecell);
3c1d9afb 1705
1706 }
1707 //printf(" \n ");
1708
1709 gStyle->SetPadRightMargin(0.1);
1710 gStyle->SetPadLeftMargin(0.1);
1711 gStyle->SetOptStat(0);
1712 gStyle->SetOptFit(000000);
1713
1714 if(maxCol-minCol > maxRow-minRow)
1715 {
1716 maxRow+= maxCol-minCol;
1717 }
1718 else
1719 {
1720 maxCol+= maxRow-minRow;
1721 }
1722
1723 TCanvas * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
1724 c->Divide(2,2);
1725 c->cd(1);
1726 gPad->SetGridy();
1727 gPad->SetGridx();
8cf5a81b 1728 gPad->SetLogz();
3c1d9afb 1729 hClusterMap ->SetAxisRange(minCol, maxCol,"X");
1730 hClusterMap ->SetAxisRange(minRow, maxRow,"Y");
8cf5a81b 1731 hClusterMap ->Draw("colz TEXT");
3c1d9afb 1732 c->cd(2);
1733 gPad->SetGridy();
1734 gPad->SetGridx();
8cf5a81b 1735 gPad->SetLogz();
3c1d9afb 1736 hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
1737 hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
8cf5a81b 1738 hClusterLocMax ->Draw("colz TEXT");
3c1d9afb 1739 c->cd(3);
1740 gPad->SetGridy();
1741 gPad->SetGridx();
8cf5a81b 1742 gPad->SetLogz();
3c1d9afb 1743 hCluster1 ->SetAxisRange(minCol, maxCol,"X");
1744 hCluster1 ->SetAxisRange(minRow, maxRow,"Y");
8cf5a81b 1745 hCluster1 ->Draw("colz TEXT");
3c1d9afb 1746 c->cd(4);
1747 gPad->SetGridy();
1748 gPad->SetGridx();
8cf5a81b 1749 gPad->SetLogz();
3c1d9afb 1750 hCluster2 ->SetAxisRange(minCol, maxCol,"X");
1751 hCluster2 ->SetAxisRange(minRow, maxRow,"Y");
8cf5a81b 1752 hCluster2 ->Draw("colz TEXT");
3c1d9afb 1753
1754 if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",
1755 eventNumber,cluster->E(),nMax,ncells1,ncells2));
1756
1757 delete c;
1758 delete hClusterMap;
1759 delete hClusterLocMax;
1760 delete hCluster1;
1761 delete hCluster2;
1762 }
1763}
1764