change default value of centrality percentile
[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
d832b695 897
cb5780f4 898//_____________________________________________________________________________________________________
765d44e7 899Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
900{
901 //Get the EMCAL/PHOS module number that corresponds to this particle
902
903 Int_t absId = -1;
904 if(particle->GetDetector()=="EMCAL"){
905 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
906 if(fDebug > 2)
907 printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
7cc7d3f8 908 particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
765d44e7 909 return fEMCALGeo->GetSuperModuleNumber(absId) ;
910 }//EMCAL
55d66f31 911 else if(particle->GetDetector()=="PHOS")
912 {
46a3cde6 913 // In case we use the MC reader, the input are TParticles,
914 // in this case use the corresponing method in PHOS Geometry to get the particle.
915 if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
916 {
917 Int_t mod =-1;
918 Double_t z = 0., x=0.;
919 TParticle* primary = 0x0;
920 AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
921 if(stack) {
922 primary = stack->Particle(particle->GetCaloLabel(0));
923 }
924 else {
280e6766 925 Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
46a3cde6 926 }
898c9d44 927
928 if(primary){
929 fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
930 }
931 else{
932 Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
933 }
46a3cde6 934 return mod;
935 }
936 // Input are ESDs or AODs, get the PHOS module number like this.
937 else{
2206e918 938 //FIXME
939 //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
940 //return GetModuleNumber(cluster);
941 //MEFIX
942 return -1;
46a3cde6 943 }
765d44e7 944 }//PHOS
945
946 return -1;
947}
948
cb5780f4 949//_____________________________________________________________________
c8fe2783 950Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
765d44e7 951{
280e6766 952 //Get the EMCAL/PHOS module number that corresponds to this cluster
765d44e7 953 TLorentzVector lv;
954 Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
55d66f31 955 if(!cluster)
956 {
9cbbc28b 957 if(fDebug > 1) printf("AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
958 return -1;
959 }
55d66f31 960
765d44e7 961 cluster->GetMomentum(lv,v);
962 Float_t phi = lv.Phi();
963 if(phi < 0) phi+=TMath::TwoPi();
964 Int_t absId = -1;
c8fe2783 965 if(cluster->IsEMCAL()){
765d44e7 966 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
967 if(fDebug > 2)
280e6766 968 printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
7cc7d3f8 969 lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
765d44e7 970 return fEMCALGeo->GetSuperModuleNumber(absId) ;
971 }//EMCAL
55d66f31 972 else if(cluster->IsPHOS())
973 {
765d44e7 974 Int_t relId[4];
55d66f31 975 if ( cluster->GetNCells() > 0)
976 {
765d44e7 977 absId = cluster->GetCellAbsId(0);
978 if(fDebug > 2)
280e6766 979 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
7cc7d3f8 980 lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
765d44e7 981 }
982 else return -1;
983
55d66f31 984 if ( absId >= 0)
985 {
765d44e7 986 fPHOSGeo->AbsToRelNumbering(absId,relId);
987 if(fDebug > 2)
280e6766 988 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
765d44e7 989 return relId[0]-1;
990 }
991 else return -1;
992 }//PHOS
993
994 return -1;
995}
996
cb5780f4 997//___________________________________________________________________________________________________
8a2dbbff 998Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(Int_t absId, TString calo,
cb5780f4 999 Int_t & icol, Int_t & irow, Int_t & iRCU) const
765d44e7 1000{
1001 //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
1002 Int_t imod = -1;
55d66f31 1003 if ( absId >= 0)
1004 {
1005 if(calo=="EMCAL")
1006 {
765d44e7 1007 Int_t iTower = -1, iIphi = -1, iIeta = -1;
1008 fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
1009 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
55d66f31 1010 if(imod < 0 || irow < 0 || icol < 0 )
1011 {
280e6766 1012 Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
1013 }
1014
765d44e7 1015 //RCU0
55d66f31 1016 if(imod < 10 )
1017 {
1018 if (0<=irow&&irow<8) iRCU=0; // first cable row
1019 else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
1020 //second cable row
1021 //RCU1
1022 else if (8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
1023 //second cable row
1024 else if (16<=irow&&irow<24) iRCU=1; // third cable row
1025
1026 if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
1027 }
1028 else
1029 {
1030 // Last 2 SM have one single SRU, just assign RCU 0
1031 iRCU = 0 ;
1032 }
1033
1034 if (iRCU<0)
1035 {
280e6766 1036 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
765d44e7 1037 }
1038
1039 return imod ;
55d66f31 1040
765d44e7 1041 }//EMCAL
55d66f31 1042 else //PHOS
1043 {
765d44e7 1044 Int_t relId[4];
1045 fPHOSGeo->AbsToRelNumbering(absId,relId);
1046 irow = relId[2];
1047 icol = relId[3];
1048 imod = relId[0]-1;
1049 iRCU= (Int_t)(relId[2]-1)/16 ;
1050 //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
55d66f31 1051 if (iRCU >= 4)
1052 {
280e6766 1053 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
765d44e7 1054 }
1055 return imod;
1056 }//PHOS
1057 }
1058
1059 return -1;
1060}
1061
7db7dcb6 1062//___________________________________________________________________________________________
71e3889f 1063Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells)
1064{
1065 // Find local maxima in cluster
1066
1067 const Int_t nc = cluster->GetNCells();
2bf17171 1068 Int_t absIdList[nc];
1069 Float_t maxEList[nc];
71e3889f 1070
1071 Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList);
1072
71e3889f 1073 return nMax;
1074
1075}
1076
1077//___________________________________________________________________________________________
7db7dcb6 1078Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
1079 Int_t *absIdList, Float_t *maxEList)
1080{
1081 // Find local maxima in cluster
f241ecc3 1082
7db7dcb6 1083 Int_t iDigitN = 0 ;
1084 Int_t iDigit = 0 ;
1085 Int_t absId1 = -1 ;
1086 Int_t absId2 = -1 ;
1087 const Int_t nCells = cluster->GetNCells();
1088
e1d24b6a 1089 Float_t eCluster = RecalibrateClusterEnergy(cluster, cells);// recalculate cluster energy, avoid non lin correction.
1090
1091 Float_t simuTotWeight = 0;
1092 if(fMCECellClusFracCorrOn)
1093 {
1094 simuTotWeight = RecalibrateClusterEnergyWeightCell(cluster, cells,eCluster);// same but apply a weight
1095 simuTotWeight/= eCluster;
1096 }
1097
7db7dcb6 1098 TString calorimeter = "EMCAL";
1099 if(!cluster->IsEMCAL()) calorimeter = "PHOS";
1100
1101 //printf("cluster : ncells %d \n",nCells);
f241ecc3 1102
1103 Float_t emax = 0;
1104 Int_t idmax =-1;
9369a2b1 1105 for(iDigit = 0; iDigit < nCells ; iDigit++)
1106 {
7db7dcb6 1107 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit] ;
f241ecc3 1108 Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
e1d24b6a 1109 RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1110
1111 if(fMCECellClusFracCorrOn)
1112 en*=GetMCECellClusFracCorrection(en,eCluster)/simuTotWeight;
1113
f241ecc3 1114 if( en > emax )
1115 {
1116 emax = en ;
1117 idmax = absIdList[iDigit] ;
1118 }
1119 //Int_t icol = -1, irow = -1, iRCU = -1;
1120 //Int_t sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1121 //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
7db7dcb6 1122 }
1123
9369a2b1 1124 for(iDigit = 0 ; iDigit < nCells; iDigit++)
1125 {
1126 if(absIdList[iDigit]>=0)
1127 {
1128 absId1 = cluster->GetCellsAbsId()[iDigit];
7db7dcb6 1129
1130 Float_t en1 = cells->GetCellAmplitude(absId1);
1131 RecalibrateCellAmplitude(en1,calorimeter,absId1);
f241ecc3 1132
e1d24b6a 1133 if(fMCECellClusFracCorrOn)
1134 en1*=GetMCECellClusFracCorrection(en1,eCluster)/simuTotWeight;
1135
9369a2b1 1136 //printf("%d : absIDi %d, E %f\n",iDigit, absId1,en1);
7db7dcb6 1137
9369a2b1 1138 for(iDigitN = 0; iDigitN < nCells; iDigitN++)
1139 {
1140 absId2 = cluster->GetCellsAbsId()[iDigitN] ;
7db7dcb6 1141
9369a2b1 1142 if(absId2==-1 || absId2==absId1) continue;
7db7dcb6 1143
9369a2b1 1144 //printf("\t %d : absIDj %d\n",iDigitN, absId2);
7db7dcb6 1145
1146 Float_t en2 = cells->GetCellAmplitude(absId2);
1147 RecalibrateCellAmplitude(en2,calorimeter,absId2);
f241ecc3 1148
e1d24b6a 1149 if(fMCECellClusFracCorrOn)
1150 en2*=GetMCECellClusFracCorrection(en2,eCluster)/simuTotWeight;
1151
9369a2b1 1152 //printf("\t %d : absIDj %d, E %f\n",iDigitN, absId2,en2);
7db7dcb6 1153
9369a2b1 1154 if ( AreNeighbours(calorimeter, absId1, absId2) )
1155 {
1156 // printf("\t \t Neighbours \n");
1157 if (en1 > en2 )
1158 {
7db7dcb6 1159 absIdList[iDigitN] = -1 ;
1160 //printf("\t \t indexN %d not local max\n",iDigitN);
1161 // but may be digit too is not local max ?
1162 if(en1 < en2 + fLocMaxCutEDiff) {
9369a2b1 1163 //printf("\t \t index %d not local max cause locMaxCutEDiff\n",iDigit);
7db7dcb6 1164 absIdList[iDigit] = -1 ;
1165 }
1166 }
9369a2b1 1167 else
1168 {
7db7dcb6 1169 absIdList[iDigit] = -1 ;
1170 //printf("\t \t index %d not local max\n",iDigitN);
1171 // but may be digitN too is not local max ?
1172 if(en1 > en2 - fLocMaxCutEDiff)
1173 {
1174 absIdList[iDigitN] = -1 ;
9369a2b1 1175 //printf("\t \t indexN %d not local max cause locMaxCutEDiff\n",iDigit);
7db7dcb6 1176 }
1177 }
9369a2b1 1178 } // if Are neighbours
1179 //else printf("\t \t NOT Neighbours \n");
7db7dcb6 1180 } // while digitN
1181 } // slot not empty
1182 } // while digit
1183
1184 iDigitN = 0 ;
9369a2b1 1185 for(iDigit = 0; iDigit < nCells; iDigit++)
1186 {
1187 if(absIdList[iDigit]>=0 )
1188 {
7db7dcb6 1189 absIdList[iDigitN] = absIdList[iDigit] ;
e1d24b6a 1190
7db7dcb6 1191 Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
e1d24b6a 1192 RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1193
1194 if(fMCECellClusFracCorrOn)
1195 en*=GetMCECellClusFracCorrection(en,eCluster)/simuTotWeight;
1196
7db7dcb6 1197 if(en < fLocMaxCutE) continue; // Maxima only with seed energy at least
e1d24b6a 1198
7db7dcb6 1199 maxEList[iDigitN] = en ;
e1d24b6a 1200
7db7dcb6 1201 //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
1202 iDigitN++ ;
1203 }
1204 }
1205
f241ecc3 1206 if(iDigitN == 0)
1207 {
1208 if(fDebug > 0)
1209 printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - No local maxima found, assign highest energy cell as maxima, id %d, en cell %2.2f, en cluster %2.2f\n",
1210 idmax,emax,cluster->E());
1211 iDigitN = 1 ;
1212 maxEList[0] = emax ;
1213 absIdList[0] = idmax ;
1214 }
1215
1216 if(fDebug > 0)
1217 {
e1d24b6a 1218 printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - In cluster E %2.2f (wth non lin. %2.2f), M02 %2.2f, M20 %2.2f, N maxima %d \n",
1219 cluster->E(),eCluster, cluster->GetM02(),cluster->GetM20(), iDigitN);
f241ecc3 1220
1221 if(fDebug > 1) for(Int_t imax = 0; imax < iDigitN; imax++)
1222 {
1223 printf(" \t i %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
1224 }
1225 }
7db7dcb6 1226
1227 return iDigitN ;
1228
1229}
1230
55d66f31 1231//____________________________________
1232TString AliCalorimeterUtils::GetPass()
1233{
1234 // Get passx from filename.
1235
1236 if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1237 {
1238 AliError("AliCalorimeterUtils::GetPass() - Pointer to tree = 0, returning null\n");
1239 return TString("");
1240 }
1241
1242 if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1243 {
1244 AliError("AliCalorimeterUtils::GetPass() - Null pointer input file, returning null\n");
1245 return TString("");
1246 }
1247
1248 TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1249 if (pass.Contains("ass1")) return TString("pass1");
1250 else if (pass.Contains("ass2")) return TString("pass2");
1251 else if (pass.Contains("ass3")) return TString("pass3");
4898bb87 1252 else if (pass.Contains("ass4")) return TString("pass4");
1253 else if (pass.Contains("ass5")) return TString("pass5");
9783326d 1254 else if (pass.Contains("LHC11c") && pass.Contains("spc_calo") ) return TString("spc_calo");
fca053d8 1255 else if (pass.Contains("calo") || pass.Contains("high_lumi"))
1256 {
1257 printf("AliCalorimeterUtils::GetPass() - Path contains <calo> or <high-lumi>, set as <pass1>\n");
1258 return TString("pass1");
1259 }
55d66f31 1260
1261 // No condition fullfilled, give a default value
1262 printf("AliCalorimeterUtils::GetPass() - Pass number string not found \n");
1263 return TString("");
1264
1265}
1266
cb5780f4 1267//________________________________________
765d44e7 1268void AliCalorimeterUtils::InitParameters()
1269{
1270 //Initialize the parameters of the analysis.
7db7dcb6 1271
12e91266 1272 fEMCALGeoName = "";
1273 fPHOSGeoName = "";
7db7dcb6 1274
7cc7d3f8 1275 fEMCALGeoMatrixSet = kFALSE;
1276 fPHOSGeoMatrixSet = kFALSE;
7db7dcb6 1277
7cc7d3f8 1278 fRemoveBadChannels = kFALSE;
7db7dcb6 1279
7cc7d3f8 1280 fNCellsFromPHOSBorder = 0;
a5fb4114 1281
7db7dcb6 1282 fLocMaxCutE = 0.1 ;
1283 fLocMaxCutEDiff = 0.0 ;
1284
7cc7d3f8 1285 // fMaskCellColumns = new Int_t[fNMaskCellColumns];
1286 // fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
1287 // fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
1288 // fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
1289 // fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
1290 // fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
a5fb4114 1291
55d66f31 1292 fOADBSet = kFALSE;
1293 fOADBForEMCAL = kTRUE ;
1294 fOADBForPHOS = kFALSE;
1295
1296 fOADBFilePathEMCAL = "$ALICE_ROOT/OADB/EMCAL" ;
1297 fOADBFilePathPHOS = "$ALICE_ROOT/OADB/PHOS" ;
1298
9e536695 1299 fImportGeometryFromFile = kTRUE;
1300 fImportGeometryFilePath = "";
190f6f2d 1301
1302 fNSuperModulesUsed = 22;
9e536695 1303
e1d24b6a 1304 fMCECellClusFracCorrParam[0] = 0.78;
1305 fMCECellClusFracCorrParam[1] =-1.8;
1306 fMCECellClusFracCorrParam[2] =-6.3;
1307 fMCECellClusFracCorrParam[3] = 0.014;
1308
765d44e7 1309}
1310
765d44e7 1311
cb5780f4 1312//_____________________________________________________
1313void AliCalorimeterUtils::InitPHOSBadChannelStatusMap()
1314{
765d44e7 1315 //Init PHOS bad channels map
1316 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
1317 //In order to avoid rewriting the same histograms
1318 Bool_t oldStatus = TH1::AddDirectoryStatus();
1319 TH1::AddDirectory(kFALSE);
7cc7d3f8 1320
78219bac 1321 fPHOSBadChannelMap = new TObjArray(5);
c5693f62 1322 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 1323
765d44e7 1324 fPHOSBadChannelMap->SetOwner(kTRUE);
1325 fPHOSBadChannelMap->Compress();
1326
1327 //In order to avoid rewriting the same histograms
1328 TH1::AddDirectory(oldStatus);
1329}
1330
cb5780f4 1331//______________________________________________________
1332void AliCalorimeterUtils::InitPHOSRecalibrationFactors()
1333{
09e819c9 1334 //Init EMCAL recalibration factors
1335 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
1336 //In order to avoid rewriting the same histograms
1337 Bool_t oldStatus = TH1::AddDirectoryStatus();
1338 TH1::AddDirectory(kFALSE);
7cc7d3f8 1339
78219bac 1340 fPHOSRecalibrationFactors = new TObjArray(5);
c5693f62 1341 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 1342 //Init the histograms with 1
1343 for (Int_t m = 0; m < 5; m++) {
1344 for (Int_t i = 0; i < 56; i++) {
1345 for (Int_t j = 0; j < 64; j++) {
c5693f62 1346 SetPHOSChannelRecalibrationFactor(m,j,i,1.);
09e819c9 1347 }
1348 }
1349 }
1350 fPHOSRecalibrationFactors->SetOwner(kTRUE);
1351 fPHOSRecalibrationFactors->Compress();
1352
1353 //In order to avoid rewriting the same histograms
1354 TH1::AddDirectory(oldStatus);
1355}
1356
1357
55d66f31 1358//__________________________________________________________
1359void AliCalorimeterUtils::InitEMCALGeometry(Int_t runnumber)
765d44e7 1360{
1361 //Initialize EMCAL geometry if it did not exist previously
55d66f31 1362
1363 if (!fEMCALGeo)
1364 {
1365 if(fEMCALGeoName=="")
1366 {
5408e59e 1367 if (runnumber < 140000 &&
1368 runnumber >= 100000) fEMCALGeoName = "EMCAL_FIRSTYEARV1";
1369 else if(runnumber >= 140000 &&
1370 runnumber < 171000) fEMCALGeoName = "EMCAL_COMPLETEV1";
1371 else fEMCALGeoName = "EMCAL_COMPLETE12SMV1";
55d66f31 1372 printf("AliCalorimeterUtils::InitEMCALGeometry() - Set EMCAL geometry name to <%s> for run %d\n",fEMCALGeoName.Data(),runnumber);
1373 }
1374
a38a48f2 1375 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
55d66f31 1376
9e536695 1377 // Init geometry, I do not like much to do it like this ...
1378 if(fImportGeometryFromFile && !gGeoManager)
1379 {
1380 if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
1381 {
1382 // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1383 if (runnumber < 140000 &&
1384 runnumber >= 100000) fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2010.root";
1385 if (runnumber >= 140000 &&
1386 runnumber < 171000) fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2011.root";
1387 else fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2012.root"; // 2012-2013
1388
1389 }
1390 printf("AliCalorimeterUtils::InitEMCALGeometry() - Import %s\n",fImportGeometryFilePath.Data());
1391 TGeoManager::Import(fImportGeometryFilePath) ; // default need file "geometry.root" in local dir!!!!
1392 }
1393
1394
55d66f31 1395 if(fDebug > 0)
1396 {
1397 printf("AliCalorimeterUtils::InitEMCALGeometry(run=%d)",runnumber);
765d44e7 1398 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1399 printf("\n");
1400 }
1401 }
1402}
1403
55d66f31 1404//_________________________________________________________
1405void AliCalorimeterUtils::InitPHOSGeometry(Int_t runnumber)
765d44e7 1406{
1407 //Initialize PHOS geometry if it did not exist previously
55d66f31 1408
1409 if (!fPHOSGeo)
1410 {
1411 if(fPHOSGeoName=="") fPHOSGeoName = "PHOSgeo";
1412
765d44e7 1413 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
55d66f31 1414
1415 if(fDebug > 0)
1416 {
1417 printf("AliCalorimeterUtils::InitPHOSGeometry(run=%d)",runnumber);
765d44e7 1418 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1419 printf("\n");
1420 }
1421 }
1422}
1423
4eb65103 1424//_______________________________________________________________________________________________
1425Bool_t AliCalorimeterUtils::IsMCParticleInCalorimeterAcceptance(TString calo, TParticle* particle)
1426{
1427 // Check that a MC ESD is in the calorimeter acceptance
1428
1429 if( (!IsPHOSGeoMatrixSet() && calo == "PHOS" ) ||
1430 (!IsPHOSGeoMatrixSet() && calo == "EMCAL") )
1431 {
1432 AliFatal(Form("Careful Geo Matrix for %s is not set, use AliFidutialCut instead \n",calo.Data()));
1433 return kFALSE ;
1434 }
1435
1436 if(calo == "PHOS" )
1437 {
1438 Int_t mod = 0 ;
1439 Double_t x = 0, z = 0 ;
1440
1441 return GetPHOSGeometry()->ImpactOnEmc( particle, mod, z, x);
1442 }
1443 else if(calo == "EMCAL")
1444 {
1445 Int_t absID = 0 ;
1446
1447 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi( particle->Eta(), particle->Phi(), absID);
1448
1449 if( absID >= 0) return kTRUE ;
1450 else return kFALSE ;
1451 }
1452
1453 return kFALSE ;
1454}
1455
1456//______________________________________________________________________________________________________
1457Bool_t AliCalorimeterUtils::IsMCParticleInCalorimeterAcceptance(TString calo, AliAODMCParticle* particle)
1458{
1459 // Check that a MC AOD is in the calorimeter acceptance
1460
1461 if( (!IsPHOSGeoMatrixSet() && calo == "PHOS" ) ||
1462 (!IsPHOSGeoMatrixSet() && calo == "EMCAL") )
1463 {
1464 AliFatal(Form("Careful Geo Matrix for %s is not set, use AliFidutialCut instead \n",calo.Data()));
1465 return kFALSE ;
1466 }
1467
1468 if(calo == "PHOS" )
1469 {
1470 Int_t mod = 0 ;
1471 Double_t x = 0, z = 0 ;
1472
1473 Double_t vtx[]={ particle->Xv(), particle->Yv(), particle->Zv() } ;
1474
1475 return GetPHOSGeometry()->ImpactOnEmc(vtx, particle->Theta(), particle->Phi(), mod, z, x) ;
1476 }
1477 else if(calo == "EMCAL")
1478 {
1479 Int_t absID = 0 ;
1480
1481 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(),absID);
1482
1483 if( absID >= 0) return kTRUE ;
1484 else return kFALSE ;
1485 }
1486
1487 return kFALSE ;
1488}
1489
8a2dbbff 1490//_______________________________________________________________________
1491Bool_t AliCalorimeterUtils::MaskFrameCluster(Int_t iSM, Int_t ieta) const
55d66f31 1492{
1493 //Check if cell is in one of the regions where we have significant amount
1494 //of material in front. Only EMCAL
1495
1496 Int_t icol = ieta;
1497 if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
1498
1499 if (fNMaskCellColumns && fMaskCellColumns)
1500 {
1501 for (Int_t imask = 0; imask < fNMaskCellColumns; imask++)
1502 {
1503 if(icol==fMaskCellColumns[imask]) return kTRUE;
1504 }
1505 }
1506
1507 return kFALSE;
1508
1509}
1510
cb5780f4 1511//_________________________________________________________
765d44e7 1512void AliCalorimeterUtils::Print(const Option_t * opt) const
1513{
7cc7d3f8 1514
765d44e7 1515 //Print some relevant parameters set for the analysis
1516 if(! opt)
1517 return;
7cc7d3f8 1518
765d44e7 1519 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1520 printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
1521 printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
7cc7d3f8 1522 fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
247abff4 1523 if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
7bf608c9 1524 printf("Recalibrate Clusters? %d, run by run %d\n",fRecalibration,fRunDependentCorrection);
9584c261 1525 printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
1526 printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
f2ccb5b8 1527 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
7cc7d3f8 1528
7db7dcb6 1529 printf("Loc. Max. E > %2.2f\n", fLocMaxCutE);
1530 printf("Loc. Max. E Diff > %2.2f\n", fLocMaxCutEDiff);
1531
765d44e7 1532 printf(" \n") ;
1533}
1534
8a2dbbff 1535//_____________________________________________________________________________________________
1536void AliCalorimeterUtils::RecalibrateCellAmplitude(Float_t & amp, TString calo, Int_t id) const
7db7dcb6 1537{
1538 //Recaculate cell energy if recalibration factor
1539
1540 Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
1541 Int_t nModule = GetModuleNumberCellIndexes(id,calo, icol, irow, iRCU);
1542
1543 if (IsRecalibrationOn())
1544 {
1545 if(calo == "PHOS")
1546 {
1547 amp *= GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
1548 }
1549 else
1550 {
1551 amp *= GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
1552 }
1553 }
1554}
1555
8a2dbbff 1556//____________________________________________________________________________________________________
1557void AliCalorimeterUtils::RecalibrateCellTime(Double_t & time, TString calo, Int_t id, Int_t bc) const
7db7dcb6 1558{
1559 // Recalculate time if time recalibration available for EMCAL
1560 // not ready for PHOS
1561
1562 if(calo == "EMCAL" && GetEMCALRecoUtils()->IsTimeRecalibrationOn())
1563 {
1564 GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
1565 }
1566
1567}
1568
cb5780f4 1569//__________________________________________________________________________
1570Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster,
1571 AliVCaloCells * cells)
1572{
09e819c9 1573 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
7cc7d3f8 1574
898c9d44 1575 //Initialize some used variables
7db7dcb6 1576 Float_t frac = 0., energy = 0.;
898c9d44 1577
7db7dcb6 1578 if(cells)
1579 {
7cc7d3f8 1580 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
7db7dcb6 1581
7cc7d3f8 1582 UShort_t * index = cluster->GetCellsAbsId() ;
1583 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
7db7dcb6 1584
7cc7d3f8 1585 Int_t ncells = cluster->GetNCells();
7db7dcb6 1586
7cc7d3f8 1587 TString calo = "EMCAL";
7db7dcb6 1588 if(cluster->IsPHOS())
1589 calo = "PHOS";
7cc7d3f8 1590
1591 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1592 for(Int_t icell = 0; icell < ncells; icell++){
7db7dcb6 1593
1594 Int_t absId = index[icell];
1595
7cc7d3f8 1596 frac = fraction[icell];
1597 if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
7db7dcb6 1598
1599 Float_t amp = cells->GetCellAmplitude(absId);
1600 RecalibrateCellAmplitude(amp,calo, absId);
1601
7cc7d3f8 1602 if(fDebug>2)
7db7dcb6 1603 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, cell fraction %f, cell energy %f\n",
1604 calo.Data(),frac,cells->GetCellAmplitude(absId));
7cc7d3f8 1605
7db7dcb6 1606 energy += amp*frac;
7cc7d3f8 1607 }
1608
1609 if(fDebug>1)
1610 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
898c9d44 1611
1612 }// cells available
7db7dcb6 1613 else
1614 {
898c9d44 1615 Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
1616 }
1617
09e819c9 1618 return energy;
09e819c9 1619}
1620
e1d24b6a 1621//__________________________________________________________________________
1622Float_t AliCalorimeterUtils::RecalibrateClusterEnergyWeightCell(AliVCluster * cluster,
1623 AliVCaloCells * cells, Float_t energyOrg)
1624{
1625 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
1626 // Also consider reweighting of cells energy
1627
1628 //Initialize some used variables
1629 Float_t frac = 0., energy = 0.;
1630
1631 if(cells)
1632 {
1633 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1634
1635 UShort_t * index = cluster->GetCellsAbsId() ;
1636 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1637
1638 Int_t ncells = cluster->GetNCells();
1639
1640 TString calo = "EMCAL";
1641 if(cluster->IsPHOS())
1642 calo = "PHOS";
1643
1644 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1645 for(Int_t icell = 0; icell < ncells; icell++){
1646
1647 Int_t absId = index[icell];
1648
1649 frac = fraction[icell];
1650 if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
1651
1652 Float_t amp = cells->GetCellAmplitude(absId);
1653 RecalibrateCellAmplitude(amp,calo, absId);
1654
1655 amp*=GetMCECellClusFracCorrection(amp,energyOrg);
1656
1657 if(fDebug>2)
1658 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, cell fraction %f, cell energy %f\n",
1659 calo.Data(),frac,cells->GetCellAmplitude(absId));
1660
1661 energy += amp*frac;
1662 }
1663
1664 if(fDebug>1)
1665 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
1666
1667 }// cells available
1668 else
1669 {
1670 Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
1671 }
1672
1673 return energy;
1674}
1675
1676
cb5780f4 1677//__________________________________________________________________________________________
1678void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
1679{
9584c261 1680
1681 //Recalculate EMCAL cluster position
1682
19db8f8c 1683 fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
7cc7d3f8 1684
9584c261 1685}
cb5780f4 1686
55d66f31 1687//________________________________________________________________________________
1688void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event,
1689 TObjArray* clusterArray)
1690{
1691 //Recalculate track matching
1692
1693 if (fRecalculateMatching)
1694 {
1695 fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo) ;
1696 //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1697 //if(esdevent){
1698 // fEMCALRecoUtils->SetClusterMatchedToTrack(esdevent) ;
1699 // fEMCALRecoUtils->SetTracksMatchedToCluster(esdevent) ;
1700 //}
1701 }
1702}
1703
3c1d9afb 1704//___________________________________________________________________________
8a2dbbff 1705void AliCalorimeterUtils::SplitEnergy(Int_t absId1, Int_t absId2,
1706 AliVCluster* cluster,
3c1d9afb 1707 AliVCaloCells* cells,
1708 //Float_t & e1, Float_t & e2,
1709 AliAODCaloCluster* cluster1,
1710 AliAODCaloCluster* cluster2,
8a2dbbff 1711 Int_t nMax, Int_t eventNumber)
3c1d9afb 1712{
1713
1714 // Split energy of cluster between the 2 local maxima, sum energy on 3x3, and if the 2
1715 // maxima are too close and have common cells, split the energy between the 2
1716
1717 TH2F* hClusterMap = 0 ;
1718 TH2F* hClusterLocMax = 0 ;
1719 TH2F* hCluster1 = 0 ;
1720 TH2F* hCluster2 = 0 ;
1721
1722 if(fPlotCluster)
1723 {
1724 hClusterMap = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
1725 hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
1726 hCluster1 = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
1727 hCluster2 = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
1728 hClusterMap ->SetXTitle("column");
1729 hClusterMap ->SetYTitle("row");
1730 hClusterLocMax ->SetXTitle("column");
1731 hClusterLocMax ->SetYTitle("row");
1732 hCluster1 ->SetXTitle("column");
1733 hCluster1 ->SetYTitle("row");
1734 hCluster2 ->SetXTitle("column");
1735 hCluster2 ->SetYTitle("row");
1736 }
1737
1738 TString calorimeter = "EMCAL";
1739 if(cluster->IsPHOS())
1740 {
1741 calorimeter="PHOS";
1742 printf("AliCalorimeterUtils::SplitEnerg() Not supported for PHOS yet \n");
1743 return;
1744 }
1745
1746 const Int_t ncells = cluster->GetNCells();
1747 Int_t absIdList[ncells];
1748
1749 Float_t e1 = 0, e2 = 0 ;
1750 Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;
1751 Float_t eCluster = 0;
1752 Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1;
1753 for(Int_t iDigit = 0; iDigit < ncells; iDigit++ )
1754 {
1755 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
1756
1757 Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
1758 RecalibrateCellAmplitude(ec,calorimeter, absIdList[iDigit]);
1759 eCluster+=ec;
1760
1761 if(fPlotCluster)
1762 {
1763 //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
1764 sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
8a871c4f 1765 if(sm > -1 && sm < 12) // just to avoid compilation warning
1766 {
1767 if(icol > maxCol) maxCol = icol;
1768 if(icol < minCol) minCol = icol;
1769 if(irow > maxRow) maxRow = irow;
1770 if(irow < minRow) minRow = irow;
1771 hClusterMap->Fill(icol,irow,ec);
1772 }
3c1d9afb 1773 }
1774
1775 }
1776
1777 // Init counters and variables
1778 Int_t ncells1 = 1 ;
1779 UShort_t absIdList1[9] ;
1780 Double_t fracList1 [9] ;
1781 absIdList1[0] = absId1 ;
1782 fracList1 [0] = 1. ;
1783
1784 Float_t ecell1 = cells->GetCellAmplitude(absId1);
1785 RecalibrateCellAmplitude(ecell1, calorimeter, absId1);
1786 e1 = ecell1;
1787
1788 Int_t ncells2 = 1 ;
1789 UShort_t absIdList2[9] ;
1790 Double_t fracList2 [9] ;
1791 absIdList2[0] = absId2 ;
1792 fracList2 [0] = 1. ;
1793
1794 Float_t ecell2 = cells->GetCellAmplitude(absId2);
1795 RecalibrateCellAmplitude(ecell2, calorimeter, absId2);
1796 e2 = ecell2;
1797
1798 if(fPlotCluster)
1799 {
1800 Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
1801 sm = GetModuleNumberCellIndexes(absId1, calorimeter, icol1, irow1, iRCU) ;
1802 hClusterLocMax->Fill(icol1,irow1,ecell1);
1803 sm = GetModuleNumberCellIndexes(absId2, calorimeter, icol2, irow2, iRCU) ;
1804 hClusterLocMax->Fill(icol2,irow2,ecell2);
1805 }
1806
1807 // Very rough way to share the cluster energy
1808 Float_t eRemain = (eCluster-ecell1-ecell2)/2;
1809 Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
1810 Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
1811
1812 for(Int_t iDigit = 0; iDigit < ncells; iDigit++)
1813 {
1814 Int_t absId = absIdList[iDigit];
1815
1816 if(absId==absId1 || absId==absId2 || absId < 0) continue;
1817
1818 Float_t ecell = cells->GetCellAmplitude(absId);
1819 RecalibrateCellAmplitude(ecell, calorimeter, absId);
1820
1821 if(AreNeighbours(calorimeter, absId1,absId ))
1822 {
1823 absIdList1[ncells1]= absId;
1824
1825 if(AreNeighbours(calorimeter, absId2,absId ))
1826 {
1827 fracList1[ncells1] = shareFraction1;
1828 e1 += ecell*shareFraction1;
1829 }
1830 else
1831 {
1832 fracList1[ncells1] = 1.;
1833 e1 += ecell;
1834 }
1835
1836 ncells1++;
1837
1838 } // neigbour to cell1
1839
1840 if(AreNeighbours(calorimeter, absId2,absId ))
1841 {
1842 absIdList2[ncells2]= absId;
1843
1844 if(AreNeighbours(calorimeter, absId1,absId ))
1845 {
1846 fracList2[ncells2] = shareFraction2;
1847 e2 += ecell*shareFraction2;
1848 }
1849 else
1850 {
1851 fracList2[ncells2] = 1.;
1852 e2 += ecell;
1853 }
1854
1855 ncells2++;
1856
1857 } // neigbour to cell2
1858
1859 }
1860
715d1872 1861 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 1862 nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2);
1863
1864 cluster1->SetE(e1);
1865 cluster2->SetE(e2);
1866
1867 cluster1->SetNCells(ncells1);
1868 cluster2->SetNCells(ncells2);
1869
1870 cluster1->SetCellsAbsId(absIdList1);
1871 cluster2->SetCellsAbsId(absIdList2);
1872
1873 cluster1->SetCellsAmplitudeFraction(fracList1);
1874 cluster2->SetCellsAmplitudeFraction(fracList2);
1875
715d1872 1876 //Correct linearity
1877 CorrectClusterEnergy(cluster1) ;
1878 CorrectClusterEnergy(cluster2) ;
1879
3c1d9afb 1880 if(calorimeter=="EMCAL")
1881 {
1882 GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
1883 GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
1884 }
1885
1886 if(fPlotCluster)
1887 {
1888 //printf("Cells of cluster1: ");
1889 for(Int_t iDigit = 0; iDigit < ncells1; iDigit++ )
1890 {
1891 //printf(" %d ",absIdList1[iDigit]);
1892
1893 sm = GetModuleNumberCellIndexes(absIdList1[iDigit], calorimeter, icol, irow, iRCU) ;
1894
8cf5a81b 1895 Float_t ecell = cells->GetCellAmplitude(absIdList1[iDigit]);
1896 RecalibrateCellAmplitude(ecell, calorimeter, absIdList1[iDigit]);
1897
1898 if( AreNeighbours(calorimeter, absId2,absIdList1[iDigit]) && absId1!=absIdList1[iDigit])
1899 hCluster1->Fill(icol,irow,ecell*shareFraction1);
3c1d9afb 1900 else
8cf5a81b 1901 hCluster1->Fill(icol,irow,ecell);
3c1d9afb 1902 }
1903
1904 //printf(" \n ");
1905 //printf("Cells of cluster2: ");
1906
1907 for(Int_t iDigit = 0; iDigit < ncells2; iDigit++ )
1908 {
1909 //printf(" %d ",absIdList2[iDigit]);
1910
1911 sm = GetModuleNumberCellIndexes(absIdList2[iDigit], calorimeter, icol, irow, iRCU) ;
8cf5a81b 1912
1913 Float_t ecell = cells->GetCellAmplitude(absIdList2[iDigit]);
1914 RecalibrateCellAmplitude(ecell, calorimeter, absIdList2[iDigit]);
1915
1916 if( AreNeighbours(calorimeter, absId1,absIdList2[iDigit]) && absId2!=absIdList2[iDigit])
1917 hCluster2->Fill(icol,irow,ecell*shareFraction2);
3c1d9afb 1918 else
8cf5a81b 1919 hCluster2->Fill(icol,irow,ecell);
3c1d9afb 1920
1921 }
1922 //printf(" \n ");
1923
1924 gStyle->SetPadRightMargin(0.1);
1925 gStyle->SetPadLeftMargin(0.1);
1926 gStyle->SetOptStat(0);
1927 gStyle->SetOptFit(000000);
1928
1929 if(maxCol-minCol > maxRow-minRow)
1930 {
1931 maxRow+= maxCol-minCol;
1932 }
1933 else
1934 {
1935 maxCol+= maxRow-minRow;
1936 }
1937
1938 TCanvas * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
1939 c->Divide(2,2);
1940 c->cd(1);
1941 gPad->SetGridy();
1942 gPad->SetGridx();
8cf5a81b 1943 gPad->SetLogz();
3c1d9afb 1944 hClusterMap ->SetAxisRange(minCol, maxCol,"X");
1945 hClusterMap ->SetAxisRange(minRow, maxRow,"Y");
8cf5a81b 1946 hClusterMap ->Draw("colz TEXT");
3c1d9afb 1947 c->cd(2);
1948 gPad->SetGridy();
1949 gPad->SetGridx();
8cf5a81b 1950 gPad->SetLogz();
3c1d9afb 1951 hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
1952 hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
8cf5a81b 1953 hClusterLocMax ->Draw("colz TEXT");
3c1d9afb 1954 c->cd(3);
1955 gPad->SetGridy();
1956 gPad->SetGridx();
8cf5a81b 1957 gPad->SetLogz();
3c1d9afb 1958 hCluster1 ->SetAxisRange(minCol, maxCol,"X");
1959 hCluster1 ->SetAxisRange(minRow, maxRow,"Y");
8cf5a81b 1960 hCluster1 ->Draw("colz TEXT");
3c1d9afb 1961 c->cd(4);
1962 gPad->SetGridy();
1963 gPad->SetGridx();
8cf5a81b 1964 gPad->SetLogz();
3c1d9afb 1965 hCluster2 ->SetAxisRange(minCol, maxCol,"X");
1966 hCluster2 ->SetAxisRange(minRow, maxRow,"Y");
8cf5a81b 1967 hCluster2 ->Draw("colz TEXT");
3c1d9afb 1968
1969 if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",
1970 eventNumber,cluster->E(),nMax,ncells1,ncells2));
1971
1972 delete c;
1973 delete hClusterMap;
1974 delete hClusterLocMax;
1975 delete hCluster1;
1976 delete hCluster2;
1977 }
1978}
1979