]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/CaloTrackCorrBase/AliCalorimeterUtils.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Additional modifications for Ntrue-efficiency acquisition at final selection level...
[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 {
4898bb87 276 TString passTmp = pass;
e07209ec 277 if(pass!="pass1" && pass!="pass2") passTmp = "pass2"; // TEMPORARY FIX FOR LHC11a analysis
278
4898bb87 279 TObjArray *trecalpass=(TObjArray*)trecal->FindObject(passTmp);
55d66f31 280
281 if(trecalpass)
282 {
f46af216 283 printf("AliCalorimeterUtils::SetOADBParameters() - Time Recalibrate EMCAL \n");
284 for (Int_t ibc = 0; ibc < 4; ++ibc)
55d66f31 285 {
f46af216 286 TH1F *h = GetEMCALChannelTimeRecalibrationFactors(ibc);
287
288 if (h)
289 delete h;
290
291 h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
292
293 if (!h)
55d66f31 294 {
f46af216 295 AliError(Form("Could not load hAllTimeAvBC%d",ibc));
296 continue;
297 }
298
299 h->SetDirectory(0);
300
301 SetEMCALChannelTimeRecalibrationFactors(ibc,h);
302 } // bunch crossing loop
6f557f7f 303 }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for pass\n"); // array pass ok
304 }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for run\n"); // run number array ok
55d66f31 305
306 } // Recalibration on
307
308 }// EMCAL
309
310 // PHOS
311 if(fOADBForPHOS)
312 {
313 printf("AliCalorimeterUtils::SetOADBParameters() - Get AODB parameters from PHOS in %s for run %d, and <%s> \n",fOADBFilePathPHOS.Data(),runnumber,pass.Data());
314
315 // Bad map
316 if(fRemoveBadChannels)
317 {
318 AliOADBContainer badmapContainer(Form("phosBadMap"));
319 TString fileName="$ALICE_ROOT/OADB/PHOS/PHOSBadMaps.root";
320 badmapContainer.InitFromFile(Form("%s/PHOSBadMaps.root",fOADBFilePathPHOS.Data()),"phosBadMap");
321
322 //Use a fixed run number from year 2010, this year not available yet.
323 TObjArray *maps = (TObjArray*)badmapContainer.GetObject(139000,"phosBadMap");
324 if(!maps)
325 {
326 printf("AliCalorimeterUtils::SetOADBParameters() - Can not read PHOS bad map for run %d.\n",runnumber) ;
327 }
328 else
329 {
330 printf("AliCalorimeterUtils::SetOADBParameters() - Setting PHOS bad map with name %s \n",maps->GetName()) ;
331 for(Int_t mod=1; mod<5;mod++)
332 {
333 TH2I *hbmPH = GetPHOSChannelStatusMap(mod);
334
335 if(hbmPH)
336 delete hbmPH ;
337
338 hbmPH = (TH2I*)maps->At(mod);
339
340 if(hbmPH) hbmPH->SetDirectory(0);
341
342 SetPHOSChannelStatusMap(mod-1,hbmPH);
343
344 } // modules loop
345 } // maps exist
346 } // Remove bad channels
347 } // PHOS
348
349 // Parameters already set once, so do not it again
350 fOADBSet = kTRUE;
351
352}
353
354//_____________________________________________________________
355void AliCalorimeterUtils::AccessGeometry(AliVEvent* inputEvent)
356{
357 //Set the calorimeters transformation matrices and init geometry
358
359 // First init the geometry, a priory not done before
360 Int_t runnumber = inputEvent->GetRunNumber() ;
361 InitPHOSGeometry (runnumber);
362 InitEMCALGeometry(runnumber);
363
364 //Get the EMCAL transformation geometry matrices from ESD
365 if(!fEMCALGeoMatrixSet && fEMCALGeo)
366 {
367 if(fLoadEMCALMatrices)
368 {
369 printf("AliCalorimeterUtils::AccessGeometry() - Load user defined EMCAL geometry matrices\n");
370
371 // OADB if available
372 AliOADBContainer emcGeoMat("AliEMCALgeo");
373 emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePathEMCAL.Data()),"AliEMCALgeo");
374 TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
375
376 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
377 {
378 if (!fEMCALMatrix[mod]) // Get it from OADB
379 {
380 if(fDebug > 1 )
381 printf("AliCalorimeterUtils::AccessGeometry() - EMCAL matrices SM %d, %p\n",
382 mod,((TGeoHMatrix*) matEMCAL->At(mod)));
383 //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
384
385 fEMCALMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
386 }
387
388 if(fEMCALMatrix[mod])
389 {
390 if(fDebug > 1)
391 fEMCALMatrix[mod]->Print();
392
393 fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;
394 }
395 }//SM loop
396
397 fEMCALGeoMatrixSet = kTRUE;//At least one, so good
398
399 }//Load matrices
400 else if (!gGeoManager)
401 {
402 if(fDebug > 1)
403 printf(" AliCalorimeterUtils::AccessGeometry() - Load EMCAL misalignment matrices. \n");
404 if(!strcmp(inputEvent->GetName(),"AliESDEvent"))
405 {
406 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
407 {
408 //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
409 if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod))
410 {
411 fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
412 }
413 }// loop over super modules
414
415 fEMCALGeoMatrixSet = kTRUE;//At least one, so good
416
417 }//ESD as input
418 else
419 {
420 if(fDebug > 1)
421 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
422 }//AOD as input
423 }//Get matrix from data
424 else if(gGeoManager)
425 {
426 fEMCALGeoMatrixSet = kTRUE;
427 }
428 }//EMCAL geo && no geoManager
429
430 //Get the PHOS transformation geometry matrices from ESD
431 if(!fPHOSGeoMatrixSet && fPHOSGeo)
432 {
433 if(fLoadPHOSMatrices)
434 {
435 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined PHOS geometry matrices\n");
436
437 // OADB if available
438 AliOADBContainer geomContainer("phosGeo");
439 geomContainer.InitFromFile(Form("%s/PHOSGeometry.root",fOADBFilePathPHOS.Data()),"PHOSRotationMatrixes");
440 TObjArray *matPHOS = (TObjArray*)geomContainer.GetObject(139000,"PHOSRotationMatrixes");
441
442 for(Int_t mod = 0 ; mod < 5 ; mod++)
443 {
444 if (!fPHOSMatrix[mod]) // Get it from OADB
445 {
446 if(fDebug > 1 )
447 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - PHOS matrices module %d, %p\n",
448 mod,((TGeoHMatrix*) matPHOS->At(mod)));
449 //((TGeoHMatrix*) matPHOS->At(mod))->Print();
450
451 fPHOSMatrix[mod] = (TGeoHMatrix*) matPHOS->At(mod) ;
452 }
453
454 // Set it, if it exists
455 if(fPHOSMatrix[mod])
456 {
457 if(fDebug > 1 )
458 fPHOSMatrix[mod]->Print();
459
460 fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;
461 }
462 }// SM loop
463
464 fPHOSGeoMatrixSet = kTRUE;//At least one, so good
465
466 }//Load matrices
467 else if (!gGeoManager)
468 {
469 if(fDebug > 1)
470 printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
471 if(!strcmp(inputEvent->GetName(),"AliESDEvent"))
472 {
473 for(Int_t mod = 0; mod < 5; mod++)
474 {
475 if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod))
476 {
477 //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
478 fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
479 }
480 }// loop over modules
481 fPHOSGeoMatrixSet = kTRUE; //At least one so good
482 }//ESD as input
483 else
484 {
485 if(fDebug > 1)
486 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
487 }//AOD as input
488 }// get matrix from data
489 else if(gGeoManager)
490 {
491 fPHOSGeoMatrixSet = kTRUE;
492 }
493 }//PHOS geo and geoManager was not set
494
495}
496
7db7dcb6 497//______________________________________________________________________________________
498Bool_t AliCalorimeterUtils::AreNeighbours(const TString calo,
499 const Int_t absId1, const Int_t absId2 ) const
500{
501 // Tells if (true) or not (false) two cells are neighbours
502 // A neighbour is defined as being two cells which share a side or corner
503
504 Bool_t areNeighbours = kFALSE ;
505
506 Int_t iRCU1 = -1, irow1 = -1, icol1 = -1;
507 Int_t iRCU2 = -1, irow2 = -1, icol2 = -1;
508
509 Int_t rowdiff = 0, coldiff = 0;
510
9369a2b1 511 Int_t nSupMod1 = GetModuleNumberCellIndexes(absId1, calo, icol1, irow1, iRCU1);
512 Int_t nSupMod2 = GetModuleNumberCellIndexes(absId2, calo, icol2, irow2, iRCU2);
7db7dcb6 513
514 if(calo=="EMCAL" && nSupMod1!=nSupMod2)
515 {
516 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
517 // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
518 if(nSupMod1%2) icol1+=AliEMCALGeoParams::fgkEMCALCols;
9369a2b1 519 else icol2+=AliEMCALGeoParams::fgkEMCALCols;
7db7dcb6 520 }
521
522 rowdiff = TMath::Abs( irow1 - irow2 ) ;
523 coldiff = TMath::Abs( icol1 - icol2 ) ;
524
525 if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
526 areNeighbours = kTRUE ;
527
528 return areNeighbours;
529
530}
531
532
cb5780f4 533//_____________________________________________________________________________________
534Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster,
535 AliVCaloCells* cells,
536 AliVEvent * event, Int_t iev) const
537{
7cc7d3f8 538
765d44e7 539 // Given the list of AbsId of the cluster, get the maximum cell and
540 // check if there are fNCellsFromBorder from the calorimeter border
541
7cc7d3f8 542 //If the distance to the border is 0 or negative just exit accept all clusters
247abff4 543 if(cells->GetType()==AliVCaloCells::kEMCALCell && fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder() <= 0 ) return kTRUE;
c8fe2783 544 if(cells->GetType()==AliVCaloCells::kPHOSCell && fNCellsFromPHOSBorder <= 0 ) return kTRUE;
7cc7d3f8 545
c8fe2783 546 Int_t absIdMax = -1;
765d44e7 547 Float_t ampMax = -1;
c8fe2783 548
549 AliMixedEvent * mixEvent = dynamic_cast<AliMixedEvent*> (event);
550 Int_t nMixedEvents = 0 ;
551 Int_t * cellsCumul = NULL ;
552 Int_t numberOfCells = 0 ;
553 if (mixEvent){
554 nMixedEvents = mixEvent->GetNumberOfEvents() ;
555 if (cells->GetType()==AliVCaloCells::kEMCALCell) {
556 cellsCumul = mixEvent->GetEMCALCellsCumul() ;
557 numberOfCells = mixEvent->GetNumberOfEMCALCells() ;
558 }
7cc7d3f8 559
c8fe2783 560 else if (cells->GetType()==AliVCaloCells::kPHOSCell) {
561 cellsCumul = mixEvent->GetPHOSCellsCumul() ;
562 numberOfCells = mixEvent->GetNumberOfPHOSCells() ;
563 }
898c9d44 564
565 if(cellsCumul){
566
567 Int_t startCell = cellsCumul[iev] ;
568 Int_t endCell = (iev+1 < nMixedEvents)?cellsCumul[iev+1]:numberOfCells;
c8fe2783 569 //Find cells with maximum amplitude
898c9d44 570 for(Int_t i = 0; i < cluster->GetNCells() ; i++){
571 Int_t absId = cluster->GetCellAbsId(i) ;
572 for (Int_t j = startCell; j < endCell ; j++) {
77e93dc2 573 Short_t cellNumber, mclabel;
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
781 clusterFraction = (eTot-eMax)/eTot; //Do not use cluster energy in case it was corrected for non linearity.
782
783 return absId;
784
9584c261 785}
786
d832b695 787//__________________________________________________________________________
788AliVTrack * AliCalorimeterUtils::GetMatchedTrack(const AliVCluster* cluster,
789 const AliVEvent* event,
790 const Int_t index) const
791{
792 // Get the matched track given its index, usually just the first match
793 // Since it is different for ESDs and AODs here it is a wrap method to do it
794
795 AliVTrack *track = 0;
796
797 // EMCAL case only when matching is recalculated
798 if(cluster->IsEMCAL() && IsRecalculationOfClusterTrackMatchingOn())
799 {
800 Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(cluster->GetID());
801 //printf("track index %d, cluster ID %d \n ",trackIndex,cluster->GetID());
802
803 if(trackIndex < 0 )
804 {
805 printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d, from recalculation\n", trackIndex);
806 }
807 else
808 {
809 track = dynamic_cast<AliVTrack*> (event->GetTrack(trackIndex));
810 }
811
812 return track ;
813
814 }
815
816 // Normal case, get info from ESD or AOD
817 // ESDs
818 if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
819 {
820 Int_t iESDtrack = cluster->GetTrackMatchedIndex();
821
822 if(iESDtrack < 0 )
823 {
824 printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d\n", index);
825 return 0x0;
826 }
827
828 track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
829
830 }
831 else // AODs
832 {
833 if(cluster->GetNTracksMatched() > 0 )
834 track = dynamic_cast<AliVTrack*>(cluster->GetTrackMatched(index));
835 }
836
837 return track ;
838
839}
840
cb5780f4 841//_____________________________________________________________________________________________________
765d44e7 842Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
843{
844 //Get the EMCAL/PHOS module number that corresponds to this particle
845
846 Int_t absId = -1;
847 if(particle->GetDetector()=="EMCAL"){
848 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
849 if(fDebug > 2)
850 printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
7cc7d3f8 851 particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
765d44e7 852 return fEMCALGeo->GetSuperModuleNumber(absId) ;
853 }//EMCAL
55d66f31 854 else if(particle->GetDetector()=="PHOS")
855 {
46a3cde6 856 // In case we use the MC reader, the input are TParticles,
857 // in this case use the corresponing method in PHOS Geometry to get the particle.
858 if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
859 {
860 Int_t mod =-1;
861 Double_t z = 0., x=0.;
862 TParticle* primary = 0x0;
863 AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
864 if(stack) {
865 primary = stack->Particle(particle->GetCaloLabel(0));
866 }
867 else {
280e6766 868 Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
46a3cde6 869 }
898c9d44 870
871 if(primary){
872 fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
873 }
874 else{
875 Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
876 }
46a3cde6 877 return mod;
878 }
879 // Input are ESDs or AODs, get the PHOS module number like this.
880 else{
2206e918 881 //FIXME
882 //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
883 //return GetModuleNumber(cluster);
884 //MEFIX
885 return -1;
46a3cde6 886 }
765d44e7 887 }//PHOS
888
889 return -1;
890}
891
cb5780f4 892//_____________________________________________________________________
c8fe2783 893Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
765d44e7 894{
280e6766 895 //Get the EMCAL/PHOS module number that corresponds to this cluster
765d44e7 896 TLorentzVector lv;
897 Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
55d66f31 898 if(!cluster)
899 {
9cbbc28b 900 if(fDebug > 1) printf("AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
901 return -1;
902 }
55d66f31 903
765d44e7 904 cluster->GetMomentum(lv,v);
905 Float_t phi = lv.Phi();
906 if(phi < 0) phi+=TMath::TwoPi();
907 Int_t absId = -1;
c8fe2783 908 if(cluster->IsEMCAL()){
765d44e7 909 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
910 if(fDebug > 2)
280e6766 911 printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
7cc7d3f8 912 lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
765d44e7 913 return fEMCALGeo->GetSuperModuleNumber(absId) ;
914 }//EMCAL
55d66f31 915 else if(cluster->IsPHOS())
916 {
765d44e7 917 Int_t relId[4];
55d66f31 918 if ( cluster->GetNCells() > 0)
919 {
765d44e7 920 absId = cluster->GetCellAbsId(0);
921 if(fDebug > 2)
280e6766 922 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
7cc7d3f8 923 lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
765d44e7 924 }
925 else return -1;
926
55d66f31 927 if ( absId >= 0)
928 {
765d44e7 929 fPHOSGeo->AbsToRelNumbering(absId,relId);
930 if(fDebug > 2)
280e6766 931 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
765d44e7 932 return relId[0]-1;
933 }
934 else return -1;
935 }//PHOS
936
937 return -1;
938}
939
cb5780f4 940//___________________________________________________________________________________________________
941Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(const Int_t absId, const TString calo,
942 Int_t & icol, Int_t & irow, Int_t & iRCU) const
765d44e7 943{
944 //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
945 Int_t imod = -1;
55d66f31 946 if ( absId >= 0)
947 {
948 if(calo=="EMCAL")
949 {
765d44e7 950 Int_t iTower = -1, iIphi = -1, iIeta = -1;
951 fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
952 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
55d66f31 953 if(imod < 0 || irow < 0 || icol < 0 )
954 {
280e6766 955 Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
956 }
957
765d44e7 958 //RCU0
55d66f31 959 if(imod < 10 )
960 {
961 if (0<=irow&&irow<8) iRCU=0; // first cable row
962 else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
963 //second cable row
964 //RCU1
965 else if (8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
966 //second cable row
967 else if (16<=irow&&irow<24) iRCU=1; // third cable row
968
969 if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
970 }
971 else
972 {
973 // Last 2 SM have one single SRU, just assign RCU 0
974 iRCU = 0 ;
975 }
976
977 if (iRCU<0)
978 {
280e6766 979 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
765d44e7 980 }
981
982 return imod ;
55d66f31 983
765d44e7 984 }//EMCAL
55d66f31 985 else //PHOS
986 {
765d44e7 987 Int_t relId[4];
988 fPHOSGeo->AbsToRelNumbering(absId,relId);
989 irow = relId[2];
990 icol = relId[3];
991 imod = relId[0]-1;
992 iRCU= (Int_t)(relId[2]-1)/16 ;
993 //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
55d66f31 994 if (iRCU >= 4)
995 {
280e6766 996 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
765d44e7 997 }
998 return imod;
999 }//PHOS
1000 }
1001
1002 return -1;
1003}
1004
71e3889f 1005//___________________________________________________________________________________________
1006Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells)
1007{
1008 // Find local maxima in cluster
1009
1010 const Int_t nc = cluster->GetNCells();
2bf17171 1011 Int_t absIdList[nc];
1012 Float_t maxEList[nc];
71e3889f 1013
1014 Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList);
1015
71e3889f 1016 return nMax;
1017
1018}
1019
7db7dcb6 1020//___________________________________________________________________________________________
1021Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
1022 Int_t *absIdList, Float_t *maxEList)
1023{
1024 // Find local maxima in cluster
f241ecc3 1025
7db7dcb6 1026 Int_t iDigitN = 0 ;
1027 Int_t iDigit = 0 ;
1028 Int_t absId1 = -1 ;
1029 Int_t absId2 = -1 ;
1030 const Int_t nCells = cluster->GetNCells();
1031
1032 TString calorimeter = "EMCAL";
1033 if(!cluster->IsEMCAL()) calorimeter = "PHOS";
1034
1035 //printf("cluster : ncells %d \n",nCells);
f241ecc3 1036
1037 Float_t emax = 0;
1038 Int_t idmax =-1;
9369a2b1 1039 for(iDigit = 0; iDigit < nCells ; iDigit++)
1040 {
7db7dcb6 1041 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit] ;
f241ecc3 1042 Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1043 RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1044 if( en > emax )
1045 {
1046 emax = en ;
1047 idmax = absIdList[iDigit] ;
1048 }
1049 //Int_t icol = -1, irow = -1, iRCU = -1;
1050 //Int_t sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1051 //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
7db7dcb6 1052 }
1053
9369a2b1 1054 for(iDigit = 0 ; iDigit < nCells; iDigit++)
1055 {
1056 if(absIdList[iDigit]>=0)
1057 {
1058 absId1 = cluster->GetCellsAbsId()[iDigit];
7db7dcb6 1059
1060 Float_t en1 = cells->GetCellAmplitude(absId1);
1061 RecalibrateCellAmplitude(en1,calorimeter,absId1);
f241ecc3 1062
9369a2b1 1063 //printf("%d : absIDi %d, E %f\n",iDigit, absId1,en1);
7db7dcb6 1064
9369a2b1 1065 for(iDigitN = 0; iDigitN < nCells; iDigitN++)
1066 {
1067 absId2 = cluster->GetCellsAbsId()[iDigitN] ;
7db7dcb6 1068
9369a2b1 1069 if(absId2==-1 || absId2==absId1) continue;
7db7dcb6 1070
9369a2b1 1071 //printf("\t %d : absIDj %d\n",iDigitN, absId2);
7db7dcb6 1072
1073 Float_t en2 = cells->GetCellAmplitude(absId2);
1074 RecalibrateCellAmplitude(en2,calorimeter,absId2);
f241ecc3 1075
9369a2b1 1076 //printf("\t %d : absIDj %d, E %f\n",iDigitN, absId2,en2);
7db7dcb6 1077
9369a2b1 1078 if ( AreNeighbours(calorimeter, absId1, absId2) )
1079 {
1080 // printf("\t \t Neighbours \n");
1081 if (en1 > en2 )
1082 {
7db7dcb6 1083 absIdList[iDigitN] = -1 ;
1084 //printf("\t \t indexN %d not local max\n",iDigitN);
1085 // but may be digit too is not local max ?
1086 if(en1 < en2 + fLocMaxCutEDiff) {
9369a2b1 1087 //printf("\t \t index %d not local max cause locMaxCutEDiff\n",iDigit);
7db7dcb6 1088 absIdList[iDigit] = -1 ;
1089 }
1090 }
9369a2b1 1091 else
1092 {
7db7dcb6 1093 absIdList[iDigit] = -1 ;
1094 //printf("\t \t index %d not local max\n",iDigitN);
1095 // but may be digitN too is not local max ?
1096 if(en1 > en2 - fLocMaxCutEDiff)
1097 {
1098 absIdList[iDigitN] = -1 ;
9369a2b1 1099 //printf("\t \t indexN %d not local max cause locMaxCutEDiff\n",iDigit);
7db7dcb6 1100 }
1101 }
9369a2b1 1102 } // if Are neighbours
1103 //else printf("\t \t NOT Neighbours \n");
7db7dcb6 1104 } // while digitN
1105 } // slot not empty
1106 } // while digit
1107
1108 iDigitN = 0 ;
9369a2b1 1109 for(iDigit = 0; iDigit < nCells; iDigit++)
1110 {
1111 if(absIdList[iDigit]>=0 )
1112 {
7db7dcb6 1113 absIdList[iDigitN] = absIdList[iDigit] ;
1114 Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1115 RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1116 if(en < fLocMaxCutE) continue; // Maxima only with seed energy at least
1117 maxEList[iDigitN] = en ;
1118 //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
1119 iDigitN++ ;
1120 }
1121 }
1122
f241ecc3 1123 if(iDigitN == 0)
1124 {
1125 if(fDebug > 0)
1126 printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - No local maxima found, assign highest energy cell as maxima, id %d, en cell %2.2f, en cluster %2.2f\n",
1127 idmax,emax,cluster->E());
1128 iDigitN = 1 ;
1129 maxEList[0] = emax ;
1130 absIdList[0] = idmax ;
1131 }
1132
1133 if(fDebug > 0)
1134 {
1135 printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - In cluster E %2.2f, M02 %2.2f, M20 %2.2f, N maxima %d \n",
1136 cluster->E(),cluster->GetM02(),cluster->GetM20(), iDigitN);
1137
1138 if(fDebug > 1) for(Int_t imax = 0; imax < iDigitN; imax++)
1139 {
1140 printf(" \t i %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
1141 }
1142 }
7db7dcb6 1143
1144 return iDigitN ;
1145
1146}
1147
55d66f31 1148//____________________________________
1149TString AliCalorimeterUtils::GetPass()
1150{
1151 // Get passx from filename.
1152
1153 if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1154 {
1155 AliError("AliCalorimeterUtils::GetPass() - Pointer to tree = 0, returning null\n");
1156 return TString("");
1157 }
1158
1159 if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1160 {
1161 AliError("AliCalorimeterUtils::GetPass() - Null pointer input file, returning null\n");
1162 return TString("");
1163 }
1164
1165 TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1166 if (pass.Contains("ass1")) return TString("pass1");
1167 else if (pass.Contains("ass2")) return TString("pass2");
1168 else if (pass.Contains("ass3")) return TString("pass3");
4898bb87 1169 else if (pass.Contains("ass4")) return TString("pass4");
1170 else if (pass.Contains("ass5")) return TString("pass5");
55d66f31 1171
1172 // No condition fullfilled, give a default value
1173 printf("AliCalorimeterUtils::GetPass() - Pass number string not found \n");
1174 return TString("");
1175
1176}
1177
cb5780f4 1178//________________________________________
765d44e7 1179void AliCalorimeterUtils::InitParameters()
1180{
1181 //Initialize the parameters of the analysis.
7db7dcb6 1182
12e91266 1183 fEMCALGeoName = "";
1184 fPHOSGeoName = "";
7db7dcb6 1185
7cc7d3f8 1186 fEMCALGeoMatrixSet = kFALSE;
1187 fPHOSGeoMatrixSet = kFALSE;
7db7dcb6 1188
7cc7d3f8 1189 fRemoveBadChannels = kFALSE;
7db7dcb6 1190
7cc7d3f8 1191 fNCellsFromPHOSBorder = 0;
a5fb4114 1192
7db7dcb6 1193 fLocMaxCutE = 0.1 ;
1194 fLocMaxCutEDiff = 0.0 ;
1195
7cc7d3f8 1196 // fMaskCellColumns = new Int_t[fNMaskCellColumns];
1197 // fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
1198 // fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
1199 // fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
1200 // fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
1201 // fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
a5fb4114 1202
55d66f31 1203 fOADBSet = kFALSE;
1204 fOADBForEMCAL = kTRUE ;
1205 fOADBForPHOS = kFALSE;
1206
1207 fOADBFilePathEMCAL = "$ALICE_ROOT/OADB/EMCAL" ;
1208 fOADBFilePathPHOS = "$ALICE_ROOT/OADB/PHOS" ;
1209
765d44e7 1210}
1211
765d44e7 1212
cb5780f4 1213//_____________________________________________________
1214void AliCalorimeterUtils::InitPHOSBadChannelStatusMap()
1215{
765d44e7 1216 //Init PHOS bad channels map
1217 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
1218 //In order to avoid rewriting the same histograms
1219 Bool_t oldStatus = TH1::AddDirectoryStatus();
1220 TH1::AddDirectory(kFALSE);
7cc7d3f8 1221
78219bac 1222 fPHOSBadChannelMap = new TObjArray(5);
c5693f62 1223 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 1224
765d44e7 1225 fPHOSBadChannelMap->SetOwner(kTRUE);
1226 fPHOSBadChannelMap->Compress();
1227
1228 //In order to avoid rewriting the same histograms
1229 TH1::AddDirectory(oldStatus);
1230}
1231
cb5780f4 1232//______________________________________________________
1233void AliCalorimeterUtils::InitPHOSRecalibrationFactors()
1234{
09e819c9 1235 //Init EMCAL recalibration factors
1236 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
1237 //In order to avoid rewriting the same histograms
1238 Bool_t oldStatus = TH1::AddDirectoryStatus();
1239 TH1::AddDirectory(kFALSE);
7cc7d3f8 1240
78219bac 1241 fPHOSRecalibrationFactors = new TObjArray(5);
c5693f62 1242 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 1243 //Init the histograms with 1
1244 for (Int_t m = 0; m < 5; m++) {
1245 for (Int_t i = 0; i < 56; i++) {
1246 for (Int_t j = 0; j < 64; j++) {
c5693f62 1247 SetPHOSChannelRecalibrationFactor(m,j,i,1.);
09e819c9 1248 }
1249 }
1250 }
1251 fPHOSRecalibrationFactors->SetOwner(kTRUE);
1252 fPHOSRecalibrationFactors->Compress();
1253
1254 //In order to avoid rewriting the same histograms
1255 TH1::AddDirectory(oldStatus);
1256}
1257
1258
55d66f31 1259//__________________________________________________________
1260void AliCalorimeterUtils::InitEMCALGeometry(Int_t runnumber)
765d44e7 1261{
1262 //Initialize EMCAL geometry if it did not exist previously
55d66f31 1263
1264 if (!fEMCALGeo)
1265 {
1266 if(fEMCALGeoName=="")
1267 {
5408e59e 1268 if (runnumber < 140000 &&
1269 runnumber >= 100000) fEMCALGeoName = "EMCAL_FIRSTYEARV1";
1270 else if(runnumber >= 140000 &&
1271 runnumber < 171000) fEMCALGeoName = "EMCAL_COMPLETEV1";
1272 else fEMCALGeoName = "EMCAL_COMPLETE12SMV1";
55d66f31 1273 printf("AliCalorimeterUtils::InitEMCALGeometry() - Set EMCAL geometry name to <%s> for run %d\n",fEMCALGeoName.Data(),runnumber);
1274 }
1275
a38a48f2 1276 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
55d66f31 1277
1278 if(fDebug > 0)
1279 {
1280 printf("AliCalorimeterUtils::InitEMCALGeometry(run=%d)",runnumber);
765d44e7 1281 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1282 printf("\n");
1283 }
1284 }
1285}
1286
55d66f31 1287//_________________________________________________________
1288void AliCalorimeterUtils::InitPHOSGeometry(Int_t runnumber)
765d44e7 1289{
1290 //Initialize PHOS geometry if it did not exist previously
55d66f31 1291
1292 if (!fPHOSGeo)
1293 {
1294 if(fPHOSGeoName=="") fPHOSGeoName = "PHOSgeo";
1295
765d44e7 1296 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
55d66f31 1297
1298 if(fDebug > 0)
1299 {
1300 printf("AliCalorimeterUtils::InitPHOSGeometry(run=%d)",runnumber);
765d44e7 1301 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1302 printf("\n");
1303 }
1304 }
1305}
1306
55d66f31 1307//__________________________________________________________________
1308Bool_t AliCalorimeterUtils::MaskFrameCluster(const Int_t iSM,
1309 const Int_t ieta) const
1310{
1311 //Check if cell is in one of the regions where we have significant amount
1312 //of material in front. Only EMCAL
1313
1314 Int_t icol = ieta;
1315 if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
1316
1317 if (fNMaskCellColumns && fMaskCellColumns)
1318 {
1319 for (Int_t imask = 0; imask < fNMaskCellColumns; imask++)
1320 {
1321 if(icol==fMaskCellColumns[imask]) return kTRUE;
1322 }
1323 }
1324
1325 return kFALSE;
1326
1327}
1328
cb5780f4 1329//_________________________________________________________
765d44e7 1330void AliCalorimeterUtils::Print(const Option_t * opt) const
1331{
7cc7d3f8 1332
765d44e7 1333 //Print some relevant parameters set for the analysis
1334 if(! opt)
1335 return;
7cc7d3f8 1336
765d44e7 1337 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1338 printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
1339 printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
7cc7d3f8 1340 fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
247abff4 1341 if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
7bf608c9 1342 printf("Recalibrate Clusters? %d, run by run %d\n",fRecalibration,fRunDependentCorrection);
9584c261 1343 printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
1344 printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
f2ccb5b8 1345 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
7cc7d3f8 1346
7db7dcb6 1347 printf("Loc. Max. E > %2.2f\n", fLocMaxCutE);
1348 printf("Loc. Max. E Diff > %2.2f\n", fLocMaxCutEDiff);
1349
765d44e7 1350 printf(" \n") ;
1351}
1352
9369a2b1 1353//__________________________________________________________________________________________
7db7dcb6 1354void AliCalorimeterUtils::RecalibrateCellAmplitude(Float_t & amp,
9369a2b1 1355 const TString calo, const Int_t id) const
7db7dcb6 1356{
1357 //Recaculate cell energy if recalibration factor
1358
1359 Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
1360 Int_t nModule = GetModuleNumberCellIndexes(id,calo, icol, irow, iRCU);
1361
1362 if (IsRecalibrationOn())
1363 {
1364 if(calo == "PHOS")
1365 {
1366 amp *= GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
1367 }
1368 else
1369 {
1370 amp *= GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
1371 }
1372 }
1373}
1374
9369a2b1 1375//_________________________________________________________________________________
7db7dcb6 1376void AliCalorimeterUtils::RecalibrateCellTime(Double_t & time,
1377 const TString calo,
9369a2b1 1378 const Int_t id, const Int_t bc) const
7db7dcb6 1379{
1380 // Recalculate time if time recalibration available for EMCAL
1381 // not ready for PHOS
1382
1383 if(calo == "EMCAL" && GetEMCALRecoUtils()->IsTimeRecalibrationOn())
1384 {
1385 GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
1386 }
1387
1388}
1389
cb5780f4 1390//__________________________________________________________________________
1391Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster,
1392 AliVCaloCells * cells)
1393{
09e819c9 1394 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
7cc7d3f8 1395
898c9d44 1396 //Initialize some used variables
7db7dcb6 1397 Float_t frac = 0., energy = 0.;
898c9d44 1398
7db7dcb6 1399 if(cells)
1400 {
7cc7d3f8 1401 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
7db7dcb6 1402
7cc7d3f8 1403 UShort_t * index = cluster->GetCellsAbsId() ;
1404 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
7db7dcb6 1405
7cc7d3f8 1406 Int_t ncells = cluster->GetNCells();
7db7dcb6 1407
7cc7d3f8 1408 TString calo = "EMCAL";
7db7dcb6 1409 if(cluster->IsPHOS())
1410 calo = "PHOS";
7cc7d3f8 1411
1412 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1413 for(Int_t icell = 0; icell < ncells; icell++){
7db7dcb6 1414
1415 Int_t absId = index[icell];
1416
7cc7d3f8 1417 frac = fraction[icell];
1418 if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
7db7dcb6 1419
1420 Float_t amp = cells->GetCellAmplitude(absId);
1421 RecalibrateCellAmplitude(amp,calo, absId);
1422
7cc7d3f8 1423 if(fDebug>2)
7db7dcb6 1424 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, cell fraction %f, cell energy %f\n",
1425 calo.Data(),frac,cells->GetCellAmplitude(absId));
7cc7d3f8 1426
7db7dcb6 1427 energy += amp*frac;
7cc7d3f8 1428 }
1429
1430 if(fDebug>1)
1431 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
898c9d44 1432
1433 }// cells available
7db7dcb6 1434 else
1435 {
898c9d44 1436 Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
1437 }
1438
09e819c9 1439 return energy;
09e819c9 1440}
1441
cb5780f4 1442//__________________________________________________________________________________________
1443void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
1444{
9584c261 1445
1446 //Recalculate EMCAL cluster position
1447
19db8f8c 1448 fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
7cc7d3f8 1449
9584c261 1450}
cb5780f4 1451
55d66f31 1452//________________________________________________________________________________
1453void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event,
1454 TObjArray* clusterArray)
1455{
1456 //Recalculate track matching
1457
1458 if (fRecalculateMatching)
1459 {
1460 fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo) ;
1461 //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1462 //if(esdevent){
1463 // fEMCALRecoUtils->SetClusterMatchedToTrack(esdevent) ;
1464 // fEMCALRecoUtils->SetTracksMatchedToCluster(esdevent) ;
1465 //}
1466 }
1467}
1468
3c1d9afb 1469//___________________________________________________________________________
1470void AliCalorimeterUtils::SplitEnergy(const Int_t absId1, const Int_t absId2,
1471 AliVCluster* cluster,
1472 AliVCaloCells* cells,
1473 //Float_t & e1, Float_t & e2,
1474 AliAODCaloCluster* cluster1,
1475 AliAODCaloCluster* cluster2,
1476 const Int_t nMax,
1477 const Int_t eventNumber)
1478{
1479
1480 // Split energy of cluster between the 2 local maxima, sum energy on 3x3, and if the 2
1481 // maxima are too close and have common cells, split the energy between the 2
1482
1483 TH2F* hClusterMap = 0 ;
1484 TH2F* hClusterLocMax = 0 ;
1485 TH2F* hCluster1 = 0 ;
1486 TH2F* hCluster2 = 0 ;
1487
1488 if(fPlotCluster)
1489 {
1490 hClusterMap = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
1491 hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
1492 hCluster1 = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
1493 hCluster2 = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
1494 hClusterMap ->SetXTitle("column");
1495 hClusterMap ->SetYTitle("row");
1496 hClusterLocMax ->SetXTitle("column");
1497 hClusterLocMax ->SetYTitle("row");
1498 hCluster1 ->SetXTitle("column");
1499 hCluster1 ->SetYTitle("row");
1500 hCluster2 ->SetXTitle("column");
1501 hCluster2 ->SetYTitle("row");
1502 }
1503
1504 TString calorimeter = "EMCAL";
1505 if(cluster->IsPHOS())
1506 {
1507 calorimeter="PHOS";
1508 printf("AliCalorimeterUtils::SplitEnerg() Not supported for PHOS yet \n");
1509 return;
1510 }
1511
1512 const Int_t ncells = cluster->GetNCells();
1513 Int_t absIdList[ncells];
1514
1515 Float_t e1 = 0, e2 = 0 ;
1516 Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;
1517 Float_t eCluster = 0;
1518 Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1;
1519 for(Int_t iDigit = 0; iDigit < ncells; iDigit++ )
1520 {
1521 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
1522
1523 Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
1524 RecalibrateCellAmplitude(ec,calorimeter, absIdList[iDigit]);
1525 eCluster+=ec;
1526
1527 if(fPlotCluster)
1528 {
1529 //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
1530 sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1531 if(icol > maxCol) maxCol = icol;
1532 if(icol < minCol) minCol = icol;
1533 if(irow > maxRow) maxRow = irow;
1534 if(irow < minRow) minRow = irow;
1535 hClusterMap->Fill(icol,irow,ec);
1536 }
1537
1538 }
1539
1540 // Init counters and variables
1541 Int_t ncells1 = 1 ;
1542 UShort_t absIdList1[9] ;
1543 Double_t fracList1 [9] ;
1544 absIdList1[0] = absId1 ;
1545 fracList1 [0] = 1. ;
1546
1547 Float_t ecell1 = cells->GetCellAmplitude(absId1);
1548 RecalibrateCellAmplitude(ecell1, calorimeter, absId1);
1549 e1 = ecell1;
1550
1551 Int_t ncells2 = 1 ;
1552 UShort_t absIdList2[9] ;
1553 Double_t fracList2 [9] ;
1554 absIdList2[0] = absId2 ;
1555 fracList2 [0] = 1. ;
1556
1557 Float_t ecell2 = cells->GetCellAmplitude(absId2);
1558 RecalibrateCellAmplitude(ecell2, calorimeter, absId2);
1559 e2 = ecell2;
1560
1561 if(fPlotCluster)
1562 {
1563 Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
1564 sm = GetModuleNumberCellIndexes(absId1, calorimeter, icol1, irow1, iRCU) ;
1565 hClusterLocMax->Fill(icol1,irow1,ecell1);
1566 sm = GetModuleNumberCellIndexes(absId2, calorimeter, icol2, irow2, iRCU) ;
1567 hClusterLocMax->Fill(icol2,irow2,ecell2);
1568 }
1569
1570 // Very rough way to share the cluster energy
1571 Float_t eRemain = (eCluster-ecell1-ecell2)/2;
1572 Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
1573 Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
1574
1575 for(Int_t iDigit = 0; iDigit < ncells; iDigit++)
1576 {
1577 Int_t absId = absIdList[iDigit];
1578
1579 if(absId==absId1 || absId==absId2 || absId < 0) continue;
1580
1581 Float_t ecell = cells->GetCellAmplitude(absId);
1582 RecalibrateCellAmplitude(ecell, calorimeter, absId);
1583
1584 if(AreNeighbours(calorimeter, absId1,absId ))
1585 {
1586 absIdList1[ncells1]= absId;
1587
1588 if(AreNeighbours(calorimeter, absId2,absId ))
1589 {
1590 fracList1[ncells1] = shareFraction1;
1591 e1 += ecell*shareFraction1;
1592 }
1593 else
1594 {
1595 fracList1[ncells1] = 1.;
1596 e1 += ecell;
1597 }
1598
1599 ncells1++;
1600
1601 } // neigbour to cell1
1602
1603 if(AreNeighbours(calorimeter, absId2,absId ))
1604 {
1605 absIdList2[ncells2]= absId;
1606
1607 if(AreNeighbours(calorimeter, absId1,absId ))
1608 {
1609 fracList2[ncells2] = shareFraction2;
1610 e2 += ecell*shareFraction2;
1611 }
1612 else
1613 {
1614 fracList2[ncells2] = 1.;
1615 e2 += ecell;
1616 }
1617
1618 ncells2++;
1619
1620 } // neigbour to cell2
1621
1622 }
1623
715d1872 1624 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 1625 nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2);
1626
1627 cluster1->SetE(e1);
1628 cluster2->SetE(e2);
1629
1630 cluster1->SetNCells(ncells1);
1631 cluster2->SetNCells(ncells2);
1632
1633 cluster1->SetCellsAbsId(absIdList1);
1634 cluster2->SetCellsAbsId(absIdList2);
1635
1636 cluster1->SetCellsAmplitudeFraction(fracList1);
1637 cluster2->SetCellsAmplitudeFraction(fracList2);
1638
715d1872 1639 //Correct linearity
1640 CorrectClusterEnergy(cluster1) ;
1641 CorrectClusterEnergy(cluster2) ;
1642
3c1d9afb 1643 if(calorimeter=="EMCAL")
1644 {
1645 GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
1646 GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
1647 }
1648
1649 if(fPlotCluster)
1650 {
1651 //printf("Cells of cluster1: ");
1652 for(Int_t iDigit = 0; iDigit < ncells1; iDigit++ )
1653 {
1654 //printf(" %d ",absIdList1[iDigit]);
1655
1656 sm = GetModuleNumberCellIndexes(absIdList1[iDigit], calorimeter, icol, irow, iRCU) ;
1657
1658 if( AreNeighbours(calorimeter, absId2,absIdList1[iDigit]) )
1659 hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit])*shareFraction1);
1660 else
1661 hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit]));
1662 }
1663
1664 //printf(" \n ");
1665 //printf("Cells of cluster2: ");
1666
1667 for(Int_t iDigit = 0; iDigit < ncells2; iDigit++ )
1668 {
1669 //printf(" %d ",absIdList2[iDigit]);
1670
1671 sm = GetModuleNumberCellIndexes(absIdList2[iDigit], calorimeter, icol, irow, iRCU) ;
1672 if( AreNeighbours(calorimeter, absId1,absIdList2[iDigit]) )
1673 hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit])*shareFraction2);
1674 else
1675 hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit]));
1676
1677 }
1678 //printf(" \n ");
1679
1680 gStyle->SetPadRightMargin(0.1);
1681 gStyle->SetPadLeftMargin(0.1);
1682 gStyle->SetOptStat(0);
1683 gStyle->SetOptFit(000000);
1684
1685 if(maxCol-minCol > maxRow-minRow)
1686 {
1687 maxRow+= maxCol-minCol;
1688 }
1689 else
1690 {
1691 maxCol+= maxRow-minRow;
1692 }
1693
1694 TCanvas * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
1695 c->Divide(2,2);
1696 c->cd(1);
1697 gPad->SetGridy();
1698 gPad->SetGridx();
1699 hClusterMap ->SetAxisRange(minCol, maxCol,"X");
1700 hClusterMap ->SetAxisRange(minRow, maxRow,"Y");
1701 hClusterMap ->Draw("colz");
1702 c->cd(2);
1703 gPad->SetGridy();
1704 gPad->SetGridx();
1705 hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
1706 hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
1707 hClusterLocMax ->Draw("colz");
1708 c->cd(3);
1709 gPad->SetGridy();
1710 gPad->SetGridx();
1711 hCluster1 ->SetAxisRange(minCol, maxCol,"X");
1712 hCluster1 ->SetAxisRange(minRow, maxRow,"Y");
1713 hCluster1 ->Draw("colz");
1714 c->cd(4);
1715 gPad->SetGridy();
1716 gPad->SetGridx();
1717 hCluster2 ->SetAxisRange(minCol, maxCol,"X");
1718 hCluster2 ->SetAxisRange(minRow, maxRow,"Y");
1719 hCluster2 ->Draw("colz");
1720
1721 if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",
1722 eventNumber,cluster->E(),nMax,ncells1,ncells2));
1723
1724 delete c;
1725 delete hClusterMap;
1726 delete hClusterLocMax;
1727 delete hCluster1;
1728 delete hCluster2;
1729 }
1730}
1731