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