in case of LHC11c spc_calo production, data is calibrated, avoid using pass1 as the...
[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");
9783326d 1172 else if (pass.Contains("LHC11c") && pass.Contains("spc_calo") ) return TString("spc_calo");
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
9e536695 1217 fImportGeometryFromFile = kTRUE;
1218 fImportGeometryFilePath = "";
1219
765d44e7 1220}
1221
765d44e7 1222
cb5780f4 1223//_____________________________________________________
1224void AliCalorimeterUtils::InitPHOSBadChannelStatusMap()
1225{
765d44e7 1226 //Init PHOS bad channels map
1227 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
1228 //In order to avoid rewriting the same histograms
1229 Bool_t oldStatus = TH1::AddDirectoryStatus();
1230 TH1::AddDirectory(kFALSE);
7cc7d3f8 1231
78219bac 1232 fPHOSBadChannelMap = new TObjArray(5);
c5693f62 1233 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 1234
765d44e7 1235 fPHOSBadChannelMap->SetOwner(kTRUE);
1236 fPHOSBadChannelMap->Compress();
1237
1238 //In order to avoid rewriting the same histograms
1239 TH1::AddDirectory(oldStatus);
1240}
1241
cb5780f4 1242//______________________________________________________
1243void AliCalorimeterUtils::InitPHOSRecalibrationFactors()
1244{
09e819c9 1245 //Init EMCAL recalibration factors
1246 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
1247 //In order to avoid rewriting the same histograms
1248 Bool_t oldStatus = TH1::AddDirectoryStatus();
1249 TH1::AddDirectory(kFALSE);
7cc7d3f8 1250
78219bac 1251 fPHOSRecalibrationFactors = new TObjArray(5);
c5693f62 1252 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 1253 //Init the histograms with 1
1254 for (Int_t m = 0; m < 5; m++) {
1255 for (Int_t i = 0; i < 56; i++) {
1256 for (Int_t j = 0; j < 64; j++) {
c5693f62 1257 SetPHOSChannelRecalibrationFactor(m,j,i,1.);
09e819c9 1258 }
1259 }
1260 }
1261 fPHOSRecalibrationFactors->SetOwner(kTRUE);
1262 fPHOSRecalibrationFactors->Compress();
1263
1264 //In order to avoid rewriting the same histograms
1265 TH1::AddDirectory(oldStatus);
1266}
1267
1268
55d66f31 1269//__________________________________________________________
1270void AliCalorimeterUtils::InitEMCALGeometry(Int_t runnumber)
765d44e7 1271{
1272 //Initialize EMCAL geometry if it did not exist previously
55d66f31 1273
1274 if (!fEMCALGeo)
1275 {
1276 if(fEMCALGeoName=="")
1277 {
5408e59e 1278 if (runnumber < 140000 &&
1279 runnumber >= 100000) fEMCALGeoName = "EMCAL_FIRSTYEARV1";
1280 else if(runnumber >= 140000 &&
1281 runnumber < 171000) fEMCALGeoName = "EMCAL_COMPLETEV1";
1282 else fEMCALGeoName = "EMCAL_COMPLETE12SMV1";
55d66f31 1283 printf("AliCalorimeterUtils::InitEMCALGeometry() - Set EMCAL geometry name to <%s> for run %d\n",fEMCALGeoName.Data(),runnumber);
1284 }
1285
a38a48f2 1286 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
55d66f31 1287
9e536695 1288 // Init geometry, I do not like much to do it like this ...
1289 if(fImportGeometryFromFile && !gGeoManager)
1290 {
1291 if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
1292 {
1293 // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1294 if (runnumber < 140000 &&
1295 runnumber >= 100000) fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2010.root";
1296 if (runnumber >= 140000 &&
1297 runnumber < 171000) fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2011.root";
1298 else fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2012.root"; // 2012-2013
1299
1300 }
1301 printf("AliCalorimeterUtils::InitEMCALGeometry() - Import %s\n",fImportGeometryFilePath.Data());
1302 TGeoManager::Import(fImportGeometryFilePath) ; // default need file "geometry.root" in local dir!!!!
1303 }
1304
1305
55d66f31 1306 if(fDebug > 0)
1307 {
1308 printf("AliCalorimeterUtils::InitEMCALGeometry(run=%d)",runnumber);
765d44e7 1309 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1310 printf("\n");
1311 }
1312 }
1313}
1314
55d66f31 1315//_________________________________________________________
1316void AliCalorimeterUtils::InitPHOSGeometry(Int_t runnumber)
765d44e7 1317{
1318 //Initialize PHOS geometry if it did not exist previously
55d66f31 1319
1320 if (!fPHOSGeo)
1321 {
1322 if(fPHOSGeoName=="") fPHOSGeoName = "PHOSgeo";
1323
765d44e7 1324 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
55d66f31 1325
1326 if(fDebug > 0)
1327 {
1328 printf("AliCalorimeterUtils::InitPHOSGeometry(run=%d)",runnumber);
765d44e7 1329 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1330 printf("\n");
1331 }
1332 }
1333}
1334
8a2dbbff 1335//_______________________________________________________________________
1336Bool_t AliCalorimeterUtils::MaskFrameCluster(Int_t iSM, Int_t ieta) const
55d66f31 1337{
1338 //Check if cell is in one of the regions where we have significant amount
1339 //of material in front. Only EMCAL
1340
1341 Int_t icol = ieta;
1342 if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
1343
1344 if (fNMaskCellColumns && fMaskCellColumns)
1345 {
1346 for (Int_t imask = 0; imask < fNMaskCellColumns; imask++)
1347 {
1348 if(icol==fMaskCellColumns[imask]) return kTRUE;
1349 }
1350 }
1351
1352 return kFALSE;
1353
1354}
1355
cb5780f4 1356//_________________________________________________________
765d44e7 1357void AliCalorimeterUtils::Print(const Option_t * opt) const
1358{
7cc7d3f8 1359
765d44e7 1360 //Print some relevant parameters set for the analysis
1361 if(! opt)
1362 return;
7cc7d3f8 1363
765d44e7 1364 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1365 printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
1366 printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
7cc7d3f8 1367 fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
247abff4 1368 if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
7bf608c9 1369 printf("Recalibrate Clusters? %d, run by run %d\n",fRecalibration,fRunDependentCorrection);
9584c261 1370 printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
1371 printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
f2ccb5b8 1372 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
7cc7d3f8 1373
7db7dcb6 1374 printf("Loc. Max. E > %2.2f\n", fLocMaxCutE);
1375 printf("Loc. Max. E Diff > %2.2f\n", fLocMaxCutEDiff);
1376
765d44e7 1377 printf(" \n") ;
1378}
1379
8a2dbbff 1380//_____________________________________________________________________________________________
1381void AliCalorimeterUtils::RecalibrateCellAmplitude(Float_t & amp, TString calo, Int_t id) const
7db7dcb6 1382{
1383 //Recaculate cell energy if recalibration factor
1384
1385 Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
1386 Int_t nModule = GetModuleNumberCellIndexes(id,calo, icol, irow, iRCU);
1387
1388 if (IsRecalibrationOn())
1389 {
1390 if(calo == "PHOS")
1391 {
1392 amp *= GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
1393 }
1394 else
1395 {
1396 amp *= GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
1397 }
1398 }
1399}
1400
8a2dbbff 1401//____________________________________________________________________________________________________
1402void AliCalorimeterUtils::RecalibrateCellTime(Double_t & time, TString calo, Int_t id, Int_t bc) const
7db7dcb6 1403{
1404 // Recalculate time if time recalibration available for EMCAL
1405 // not ready for PHOS
1406
1407 if(calo == "EMCAL" && GetEMCALRecoUtils()->IsTimeRecalibrationOn())
1408 {
1409 GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
1410 }
1411
1412}
1413
cb5780f4 1414//__________________________________________________________________________
1415Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster,
1416 AliVCaloCells * cells)
1417{
09e819c9 1418 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
7cc7d3f8 1419
898c9d44 1420 //Initialize some used variables
7db7dcb6 1421 Float_t frac = 0., energy = 0.;
898c9d44 1422
7db7dcb6 1423 if(cells)
1424 {
7cc7d3f8 1425 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
7db7dcb6 1426
7cc7d3f8 1427 UShort_t * index = cluster->GetCellsAbsId() ;
1428 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
7db7dcb6 1429
7cc7d3f8 1430 Int_t ncells = cluster->GetNCells();
7db7dcb6 1431
7cc7d3f8 1432 TString calo = "EMCAL";
7db7dcb6 1433 if(cluster->IsPHOS())
1434 calo = "PHOS";
7cc7d3f8 1435
1436 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1437 for(Int_t icell = 0; icell < ncells; icell++){
7db7dcb6 1438
1439 Int_t absId = index[icell];
1440
7cc7d3f8 1441 frac = fraction[icell];
1442 if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
7db7dcb6 1443
1444 Float_t amp = cells->GetCellAmplitude(absId);
1445 RecalibrateCellAmplitude(amp,calo, absId);
1446
7cc7d3f8 1447 if(fDebug>2)
7db7dcb6 1448 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, cell fraction %f, cell energy %f\n",
1449 calo.Data(),frac,cells->GetCellAmplitude(absId));
7cc7d3f8 1450
7db7dcb6 1451 energy += amp*frac;
7cc7d3f8 1452 }
1453
1454 if(fDebug>1)
1455 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
898c9d44 1456
1457 }// cells available
7db7dcb6 1458 else
1459 {
898c9d44 1460 Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
1461 }
1462
09e819c9 1463 return energy;
09e819c9 1464}
1465
cb5780f4 1466//__________________________________________________________________________________________
1467void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
1468{
9584c261 1469
1470 //Recalculate EMCAL cluster position
1471
19db8f8c 1472 fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
7cc7d3f8 1473
9584c261 1474}
cb5780f4 1475
55d66f31 1476//________________________________________________________________________________
1477void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event,
1478 TObjArray* clusterArray)
1479{
1480 //Recalculate track matching
1481
1482 if (fRecalculateMatching)
1483 {
1484 fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo) ;
1485 //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1486 //if(esdevent){
1487 // fEMCALRecoUtils->SetClusterMatchedToTrack(esdevent) ;
1488 // fEMCALRecoUtils->SetTracksMatchedToCluster(esdevent) ;
1489 //}
1490 }
1491}
1492
3c1d9afb 1493//___________________________________________________________________________
8a2dbbff 1494void AliCalorimeterUtils::SplitEnergy(Int_t absId1, Int_t absId2,
1495 AliVCluster* cluster,
3c1d9afb 1496 AliVCaloCells* cells,
1497 //Float_t & e1, Float_t & e2,
1498 AliAODCaloCluster* cluster1,
1499 AliAODCaloCluster* cluster2,
8a2dbbff 1500 Int_t nMax, Int_t eventNumber)
3c1d9afb 1501{
1502
1503 // Split energy of cluster between the 2 local maxima, sum energy on 3x3, and if the 2
1504 // maxima are too close and have common cells, split the energy between the 2
1505
1506 TH2F* hClusterMap = 0 ;
1507 TH2F* hClusterLocMax = 0 ;
1508 TH2F* hCluster1 = 0 ;
1509 TH2F* hCluster2 = 0 ;
1510
1511 if(fPlotCluster)
1512 {
1513 hClusterMap = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
1514 hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
1515 hCluster1 = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
1516 hCluster2 = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
1517 hClusterMap ->SetXTitle("column");
1518 hClusterMap ->SetYTitle("row");
1519 hClusterLocMax ->SetXTitle("column");
1520 hClusterLocMax ->SetYTitle("row");
1521 hCluster1 ->SetXTitle("column");
1522 hCluster1 ->SetYTitle("row");
1523 hCluster2 ->SetXTitle("column");
1524 hCluster2 ->SetYTitle("row");
1525 }
1526
1527 TString calorimeter = "EMCAL";
1528 if(cluster->IsPHOS())
1529 {
1530 calorimeter="PHOS";
1531 printf("AliCalorimeterUtils::SplitEnerg() Not supported for PHOS yet \n");
1532 return;
1533 }
1534
1535 const Int_t ncells = cluster->GetNCells();
1536 Int_t absIdList[ncells];
1537
1538 Float_t e1 = 0, e2 = 0 ;
1539 Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;
1540 Float_t eCluster = 0;
1541 Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1;
1542 for(Int_t iDigit = 0; iDigit < ncells; iDigit++ )
1543 {
1544 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
1545
1546 Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
1547 RecalibrateCellAmplitude(ec,calorimeter, absIdList[iDigit]);
1548 eCluster+=ec;
1549
1550 if(fPlotCluster)
1551 {
1552 //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
1553 sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
8a871c4f 1554 if(sm > -1 && sm < 12) // just to avoid compilation warning
1555 {
1556 if(icol > maxCol) maxCol = icol;
1557 if(icol < minCol) minCol = icol;
1558 if(irow > maxRow) maxRow = irow;
1559 if(irow < minRow) minRow = irow;
1560 hClusterMap->Fill(icol,irow,ec);
1561 }
3c1d9afb 1562 }
1563
1564 }
1565
1566 // Init counters and variables
1567 Int_t ncells1 = 1 ;
1568 UShort_t absIdList1[9] ;
1569 Double_t fracList1 [9] ;
1570 absIdList1[0] = absId1 ;
1571 fracList1 [0] = 1. ;
1572
1573 Float_t ecell1 = cells->GetCellAmplitude(absId1);
1574 RecalibrateCellAmplitude(ecell1, calorimeter, absId1);
1575 e1 = ecell1;
1576
1577 Int_t ncells2 = 1 ;
1578 UShort_t absIdList2[9] ;
1579 Double_t fracList2 [9] ;
1580 absIdList2[0] = absId2 ;
1581 fracList2 [0] = 1. ;
1582
1583 Float_t ecell2 = cells->GetCellAmplitude(absId2);
1584 RecalibrateCellAmplitude(ecell2, calorimeter, absId2);
1585 e2 = ecell2;
1586
1587 if(fPlotCluster)
1588 {
1589 Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
1590 sm = GetModuleNumberCellIndexes(absId1, calorimeter, icol1, irow1, iRCU) ;
1591 hClusterLocMax->Fill(icol1,irow1,ecell1);
1592 sm = GetModuleNumberCellIndexes(absId2, calorimeter, icol2, irow2, iRCU) ;
1593 hClusterLocMax->Fill(icol2,irow2,ecell2);
1594 }
1595
1596 // Very rough way to share the cluster energy
1597 Float_t eRemain = (eCluster-ecell1-ecell2)/2;
1598 Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
1599 Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
1600
1601 for(Int_t iDigit = 0; iDigit < ncells; iDigit++)
1602 {
1603 Int_t absId = absIdList[iDigit];
1604
1605 if(absId==absId1 || absId==absId2 || absId < 0) continue;
1606
1607 Float_t ecell = cells->GetCellAmplitude(absId);
1608 RecalibrateCellAmplitude(ecell, calorimeter, absId);
1609
1610 if(AreNeighbours(calorimeter, absId1,absId ))
1611 {
1612 absIdList1[ncells1]= absId;
1613
1614 if(AreNeighbours(calorimeter, absId2,absId ))
1615 {
1616 fracList1[ncells1] = shareFraction1;
1617 e1 += ecell*shareFraction1;
1618 }
1619 else
1620 {
1621 fracList1[ncells1] = 1.;
1622 e1 += ecell;
1623 }
1624
1625 ncells1++;
1626
1627 } // neigbour to cell1
1628
1629 if(AreNeighbours(calorimeter, absId2,absId ))
1630 {
1631 absIdList2[ncells2]= absId;
1632
1633 if(AreNeighbours(calorimeter, absId1,absId ))
1634 {
1635 fracList2[ncells2] = shareFraction2;
1636 e2 += ecell*shareFraction2;
1637 }
1638 else
1639 {
1640 fracList2[ncells2] = 1.;
1641 e2 += ecell;
1642 }
1643
1644 ncells2++;
1645
1646 } // neigbour to cell2
1647
1648 }
1649
715d1872 1650 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 1651 nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2);
1652
1653 cluster1->SetE(e1);
1654 cluster2->SetE(e2);
1655
1656 cluster1->SetNCells(ncells1);
1657 cluster2->SetNCells(ncells2);
1658
1659 cluster1->SetCellsAbsId(absIdList1);
1660 cluster2->SetCellsAbsId(absIdList2);
1661
1662 cluster1->SetCellsAmplitudeFraction(fracList1);
1663 cluster2->SetCellsAmplitudeFraction(fracList2);
1664
715d1872 1665 //Correct linearity
1666 CorrectClusterEnergy(cluster1) ;
1667 CorrectClusterEnergy(cluster2) ;
1668
3c1d9afb 1669 if(calorimeter=="EMCAL")
1670 {
1671 GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
1672 GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
1673 }
1674
1675 if(fPlotCluster)
1676 {
1677 //printf("Cells of cluster1: ");
1678 for(Int_t iDigit = 0; iDigit < ncells1; iDigit++ )
1679 {
1680 //printf(" %d ",absIdList1[iDigit]);
1681
1682 sm = GetModuleNumberCellIndexes(absIdList1[iDigit], calorimeter, icol, irow, iRCU) ;
1683
8cf5a81b 1684 Float_t ecell = cells->GetCellAmplitude(absIdList1[iDigit]);
1685 RecalibrateCellAmplitude(ecell, calorimeter, absIdList1[iDigit]);
1686
1687 if( AreNeighbours(calorimeter, absId2,absIdList1[iDigit]) && absId1!=absIdList1[iDigit])
1688 hCluster1->Fill(icol,irow,ecell*shareFraction1);
3c1d9afb 1689 else
8cf5a81b 1690 hCluster1->Fill(icol,irow,ecell);
3c1d9afb 1691 }
1692
1693 //printf(" \n ");
1694 //printf("Cells of cluster2: ");
1695
1696 for(Int_t iDigit = 0; iDigit < ncells2; iDigit++ )
1697 {
1698 //printf(" %d ",absIdList2[iDigit]);
1699
1700 sm = GetModuleNumberCellIndexes(absIdList2[iDigit], calorimeter, icol, irow, iRCU) ;
8cf5a81b 1701
1702 Float_t ecell = cells->GetCellAmplitude(absIdList2[iDigit]);
1703 RecalibrateCellAmplitude(ecell, calorimeter, absIdList2[iDigit]);
1704
1705 if( AreNeighbours(calorimeter, absId1,absIdList2[iDigit]) && absId2!=absIdList2[iDigit])
1706 hCluster2->Fill(icol,irow,ecell*shareFraction2);
3c1d9afb 1707 else
8cf5a81b 1708 hCluster2->Fill(icol,irow,ecell);
3c1d9afb 1709
1710 }
1711 //printf(" \n ");
1712
1713 gStyle->SetPadRightMargin(0.1);
1714 gStyle->SetPadLeftMargin(0.1);
1715 gStyle->SetOptStat(0);
1716 gStyle->SetOptFit(000000);
1717
1718 if(maxCol-minCol > maxRow-minRow)
1719 {
1720 maxRow+= maxCol-minCol;
1721 }
1722 else
1723 {
1724 maxCol+= maxRow-minRow;
1725 }
1726
1727 TCanvas * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
1728 c->Divide(2,2);
1729 c->cd(1);
1730 gPad->SetGridy();
1731 gPad->SetGridx();
8cf5a81b 1732 gPad->SetLogz();
3c1d9afb 1733 hClusterMap ->SetAxisRange(minCol, maxCol,"X");
1734 hClusterMap ->SetAxisRange(minRow, maxRow,"Y");
8cf5a81b 1735 hClusterMap ->Draw("colz TEXT");
3c1d9afb 1736 c->cd(2);
1737 gPad->SetGridy();
1738 gPad->SetGridx();
8cf5a81b 1739 gPad->SetLogz();
3c1d9afb 1740 hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
1741 hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
8cf5a81b 1742 hClusterLocMax ->Draw("colz TEXT");
3c1d9afb 1743 c->cd(3);
1744 gPad->SetGridy();
1745 gPad->SetGridx();
8cf5a81b 1746 gPad->SetLogz();
3c1d9afb 1747 hCluster1 ->SetAxisRange(minCol, maxCol,"X");
1748 hCluster1 ->SetAxisRange(minRow, maxRow,"Y");
8cf5a81b 1749 hCluster1 ->Draw("colz TEXT");
3c1d9afb 1750 c->cd(4);
1751 gPad->SetGridy();
1752 gPad->SetGridx();
8cf5a81b 1753 gPad->SetLogz();
3c1d9afb 1754 hCluster2 ->SetAxisRange(minCol, maxCol,"X");
1755 hCluster2 ->SetAxisRange(minRow, maxRow,"Y");
8cf5a81b 1756 hCluster2 ->Draw("colz TEXT");
3c1d9afb 1757
1758 if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",
1759 eventNumber,cluster->E(),nMax,ncells1,ncells2));
1760
1761 delete c;
1762 delete hClusterMap;
1763 delete hClusterLocMax;
1764 delete hCluster1;
1765 delete hCluster2;
1766 }
1767}
1768