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