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