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