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