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