Enabled vertexing QA
[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 ---
26#include "TGeoManager.h"
27
c5693f62 28// --- ANALYSIS system ---
765d44e7 29#include "AliCalorimeterUtils.h"
4b892846 30#include "AliESDEvent.h"
46a3cde6 31#include "AliMCEvent.h"
32#include "AliStack.h"
765d44e7 33#include "AliAODPWG4Particle.h"
c8fe2783 34#include "AliVCluster.h"
35#include "AliVCaloCells.h"
36#include "AliMixedEvent.h"
765d44e7 37
c5693f62 38// --- Detector ---
39#include "AliEMCALGeometry.h"
40#include "AliPHOSGeoUtils.h"
41
765d44e7 42ClassImp(AliCalorimeterUtils)
43
44
cb5780f4 45//____________________________________________
765d44e7 46 AliCalorimeterUtils::AliCalorimeterUtils() :
47 TObject(), fDebug(0),
90e32961 48 fEMCALGeoName("EMCAL_COMPLETE12SMV1"),
49 fPHOSGeoName ("PHOSgeo"),
9e8998b1 50 fEMCALGeo(0x0), fPHOSGeo(0x0),
51 fEMCALGeoMatrixSet(kFALSE), fPHOSGeoMatrixSet(kFALSE),
52 fLoadEMCALMatrices(kFALSE), fLoadPHOSMatrices(kFALSE),
53 fRemoveBadChannels(kFALSE), fPHOSBadChannelMap(0x0),
a5fb4114 54 fNCellsFromPHOSBorder(0),
9e8998b1 55 fNMaskCellColumns(0), fMaskCellColumns(0x0),
56 fRecalibration(kFALSE), fPHOSRecalibrationFactors(),
247abff4 57 fEMCALRecoUtils(new AliEMCALRecoUtils),
9e8998b1 58 fRecalculatePosition(kFALSE), fCorrectELinearity(kFALSE),
59 fRecalculateMatching(kFALSE),
60 fCutR(20), fCutZ(20),
61 fCutEta(20), fCutPhi(20)
765d44e7 62{
63 //Ctor
64
65 //Initialize parameters
66 InitParameters();
90e32961 67 for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
68 for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix [i] = 0 ;
7cc7d3f8 69
765d44e7 70}
765d44e7 71
cb5780f4 72//_________________________________________
73AliCalorimeterUtils::~AliCalorimeterUtils()
74{
765d44e7 75 //Dtor
76
4df35693 77 //if(fPHOSGeo) delete fPHOSGeo ;
765d44e7 78 if(fEMCALGeo) delete fEMCALGeo ;
7cc7d3f8 79
765d44e7 80 if(fPHOSBadChannelMap) {
81 fPHOSBadChannelMap->Clear();
82 delete fPHOSBadChannelMap;
83 }
84
09e819c9 85 if(fPHOSRecalibrationFactors) {
78219bac 86 fPHOSRecalibrationFactors->Clear();
87 delete fPHOSRecalibrationFactors;
09e819c9 88 }
89
a5fb4114 90 if(fEMCALRecoUtils) delete fEMCALRecoUtils ;
91 if(fNMaskCellColumns) delete [] fMaskCellColumns;
9584c261 92
765d44e7 93}
94
cb5780f4 95//_____________________________________________________________________________________
96Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster,
97 AliVCaloCells* cells,
98 AliVEvent * event, Int_t iev) const
99{
7cc7d3f8 100
765d44e7 101 // Given the list of AbsId of the cluster, get the maximum cell and
102 // check if there are fNCellsFromBorder from the calorimeter border
103
7cc7d3f8 104 //If the distance to the border is 0 or negative just exit accept all clusters
247abff4 105 if(cells->GetType()==AliVCaloCells::kEMCALCell && fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder() <= 0 ) return kTRUE;
c8fe2783 106 if(cells->GetType()==AliVCaloCells::kPHOSCell && fNCellsFromPHOSBorder <= 0 ) return kTRUE;
7cc7d3f8 107
c8fe2783 108 Int_t absIdMax = -1;
765d44e7 109 Float_t ampMax = -1;
c8fe2783 110
111 AliMixedEvent * mixEvent = dynamic_cast<AliMixedEvent*> (event);
112 Int_t nMixedEvents = 0 ;
113 Int_t * cellsCumul = NULL ;
114 Int_t numberOfCells = 0 ;
115 if (mixEvent){
116 nMixedEvents = mixEvent->GetNumberOfEvents() ;
117 if (cells->GetType()==AliVCaloCells::kEMCALCell) {
118 cellsCumul = mixEvent->GetEMCALCellsCumul() ;
119 numberOfCells = mixEvent->GetNumberOfEMCALCells() ;
120 }
7cc7d3f8 121
c8fe2783 122 else if (cells->GetType()==AliVCaloCells::kPHOSCell) {
123 cellsCumul = mixEvent->GetPHOSCellsCumul() ;
124 numberOfCells = mixEvent->GetNumberOfPHOSCells() ;
125 }
898c9d44 126
127 if(cellsCumul){
128
129 Int_t startCell = cellsCumul[iev] ;
130 Int_t endCell = (iev+1 < nMixedEvents)?cellsCumul[iev+1]:numberOfCells;
c8fe2783 131 //Find cells with maximum amplitude
898c9d44 132 for(Int_t i = 0; i < cluster->GetNCells() ; i++){
133 Int_t absId = cluster->GetCellAbsId(i) ;
134 for (Int_t j = startCell; j < endCell ; j++) {
135 Short_t cellNumber;
136 Double_t amp ;
137 Double_t time;
138 cells->GetCell(j, cellNumber, amp, time) ;
139 if (absId == cellNumber) {
140 if(amp > ampMax){
141 ampMax = amp;
142 absIdMax = absId;
143 }
144 }
c8fe2783 145 }
898c9d44 146 }//loop on cluster cells
147 }// cells cumul available
148 else {
149 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - CellsCumul is NULL!!!\n");
150 abort();
c8fe2783 151 }
898c9d44 152 } else {//Normal SE Events
c8fe2783 153 for(Int_t i = 0; i < cluster->GetNCells() ; i++){
154 Int_t absId = cluster->GetCellAbsId(i) ;
155 Float_t amp = cells->GetCellAmplitude(absId);
156 if(amp > ampMax){
157 ampMax = amp;
158 absIdMax = absId;
159 }
160 }
161 }
765d44e7 162
163 if(fDebug > 1)
280e6766 164 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n",
7cc7d3f8 165 absIdMax, ampMax, cluster->E());
765d44e7 166
167 if(absIdMax==-1) return kFALSE;
168
169 //Check if the cell is close to the borders:
170 Bool_t okrow = kFALSE;
171 Bool_t okcol = kFALSE;
7cc7d3f8 172
c8fe2783 173 if(cells->GetType()==AliVCaloCells::kEMCALCell){
7cc7d3f8 174
765d44e7 175 Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1;
176 fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta);
177 fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
280e6766 178 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
179 Fatal("CheckCellFidutialRegion","Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",iSM,ieta,iphi);
180 }
7cc7d3f8 181
765d44e7 182 //Check rows/phi
247abff4 183 Int_t nborder = fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder();
765d44e7 184 if(iSM < 10){
247abff4 185 if(iphi >= nborder && iphi < 24-nborder) okrow =kTRUE;
7cc7d3f8 186 }
765d44e7 187 else{
247abff4 188 if(iphi >= nborder && iphi < 12-nborder) okrow =kTRUE;
765d44e7 189 }
190
247abff4 191 //Check columns/eta
192 if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0()){
193 if(ieta > nborder && ieta < 48-nborder) okcol =kTRUE;
765d44e7 194 }
195 else{
196 if(iSM%2==0){
247abff4 197 if(ieta >= nborder) okcol = kTRUE;
765d44e7 198 }
199 else {
247abff4 200 if(ieta < 48-nborder) okcol = kTRUE;
765d44e7 201 }
202 }//eta 0 not checked
203 if(fDebug > 1)
204 {
280e6766 205 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
7cc7d3f8 206 nborder, ieta, iphi, iSM);
765d44e7 207 if (okcol && okrow ) printf(" YES \n");
208 else printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
209 }
210 }//EMCAL
c8fe2783 211 else if(cells->GetType()==AliVCaloCells::kPHOSCell){
765d44e7 212 Int_t relId[4];
213 Int_t irow = -1, icol = -1;
214 fPHOSGeo->AbsToRelNumbering(absIdMax,relId);
215 irow = relId[2];
216 icol = relId[3];
217 //imod = relId[0]-1;
218 if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE;
219 if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE;
220 if(fDebug > 1)
221 {
280e6766 222 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - PHOS Cluster in %d cells fiducial volume: icol %d, irow %d, Module %d?",
7cc7d3f8 223 fNCellsFromPHOSBorder, icol, irow, relId[0]-1);
765d44e7 224 if (okcol && okrow ) printf(" YES \n");
225 else printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
226 }
227 }//PHOS
228
229 if (okcol && okrow) return kTRUE;
230 else return kFALSE;
231
232}
233
765d44e7 234//_________________________________________________________________________________________________________
cb5780f4 235Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter,UShort_t* cellList, Int_t nCells)
236{
765d44e7 237 // Check that in the cluster cells, there is no bad channel of those stored
238 // in fEMCALBadChannelMap or fPHOSBadChannelMap
239
240 if (!fRemoveBadChannels) return kFALSE;
36037088 241 //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
247abff4 242 if(calorimeter == "EMCAL" && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
78219bac 243 if(calorimeter == "PHOS" && !fPHOSBadChannelMap) return kFALSE;
7cc7d3f8 244
765d44e7 245 Int_t icol = -1;
246 Int_t irow = -1;
247 Int_t imod = -1;
248 for(Int_t iCell = 0; iCell<nCells; iCell++){
7cc7d3f8 249
765d44e7 250 //Get the column and row
251 if(calorimeter == "EMCAL"){
247abff4 252 return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
765d44e7 253 }
254 else if(calorimeter=="PHOS"){
255 Int_t relId[4];
256 fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
257 irow = relId[2];
258 icol = relId[3];
259 imod = relId[0]-1;
260 if(fPHOSBadChannelMap->GetEntries() <= imod)continue;
c5693f62 261 //printf("PHOS bad channels imod %d, icol %d, irow %d\n",imod, irow, icol);
262 if(GetPHOSChannelStatus(imod, irow, icol)) return kTRUE;
765d44e7 263 }
264 else return kFALSE;
265
266 }// cell cluster loop
7cc7d3f8 267
765d44e7 268 return kFALSE;
7cc7d3f8 269
765d44e7 270}
271
cb5780f4 272//_______________________________________________________________
273void AliCalorimeterUtils::CorrectClusterEnergy(AliVCluster *clus)
274{
9584c261 275 // Correct cluster energy non linearity
13cd2872 276
9584c261 277 clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
13cd2872 278
279}
280
c5693f62 281//________________________________________________________________________________________
282Int_t AliCalorimeterUtils::GetMaxEnergyCell(AliVCaloCells* cells, const AliVCluster* clu,
283 Float_t & clusterFraction) const
cb5780f4 284{
13cd2872 285
286 //For a given CaloCluster gets the absId of the cell
287 //with maximum energy deposit.
288
289 if( !clu || !cells ){
290 AliInfo("Cluster or cells pointer is null!");
291 return -1;
292 }
293
294 Double_t eMax =-1.;
295 Double_t eTot = 0.;
296 Double_t eCell =-1.;
297 Float_t fraction = 1.;
298 Float_t recalFactor = 1.;
299 Int_t cellAbsId =-1 , absId =-1 ;
300 Int_t iSupMod =-1 , ieta =-1 , iphi = -1, iRCU = -1;
301
302 TString calo = "EMCAL";
303 if(clu->IsPHOS()) calo = "PHOS";
304
305 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
306
307 cellAbsId = clu->GetCellAbsId(iDig);
308
309 fraction = clu->GetCellAmplitudeFraction(iDig);
310 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
311
b08dd6d5 312 iSupMod = GetModuleNumberCellIndexes(cellAbsId, calo, ieta, iphi, iRCU);
13cd2872 313
314 if(IsRecalibrationOn()) {
315 if(calo=="EMCAL") recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
c5693f62 316 else recalFactor = GetPHOSChannelRecalibrationFactor (iSupMod,iphi,ieta);
13cd2872 317 }
318
319 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
320
321 if(eCell > eMax) {
322 eMax = eCell;
323 absId = cellAbsId;
324 }
325
326 eTot+=eCell;
327
328 }// cell loop
329
330 clusterFraction = (eTot-eMax)/eTot; //Do not use cluster energy in case it was corrected for non linearity.
331
332 return absId;
333
9584c261 334}
335
d832b695 336//__________________________________________________________________________
337AliVTrack * AliCalorimeterUtils::GetMatchedTrack(const AliVCluster* cluster,
338 const AliVEvent* event,
339 const Int_t index) const
340{
341 // Get the matched track given its index, usually just the first match
342 // Since it is different for ESDs and AODs here it is a wrap method to do it
343
344 AliVTrack *track = 0;
345
346 // EMCAL case only when matching is recalculated
347 if(cluster->IsEMCAL() && IsRecalculationOfClusterTrackMatchingOn())
348 {
349 Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(cluster->GetID());
350 //printf("track index %d, cluster ID %d \n ",trackIndex,cluster->GetID());
351
352 if(trackIndex < 0 )
353 {
354 printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d, from recalculation\n", trackIndex);
355 }
356 else
357 {
358 track = dynamic_cast<AliVTrack*> (event->GetTrack(trackIndex));
359 }
360
361 return track ;
362
363 }
364
365 // Normal case, get info from ESD or AOD
366 // ESDs
367 if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
368 {
369 Int_t iESDtrack = cluster->GetTrackMatchedIndex();
370
371 if(iESDtrack < 0 )
372 {
373 printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d\n", index);
374 return 0x0;
375 }
376
377 track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
378
379 }
380 else // AODs
381 {
382 if(cluster->GetNTracksMatched() > 0 )
383 track = dynamic_cast<AliVTrack*>(cluster->GetTrackMatched(index));
384 }
385
386 return track ;
387
388}
389
cb5780f4 390//_____________________________________________________________________________________________________
765d44e7 391Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
392{
393 //Get the EMCAL/PHOS module number that corresponds to this particle
394
395 Int_t absId = -1;
396 if(particle->GetDetector()=="EMCAL"){
397 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
398 if(fDebug > 2)
399 printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
7cc7d3f8 400 particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
765d44e7 401 return fEMCALGeo->GetSuperModuleNumber(absId) ;
402 }//EMCAL
403 else if(particle->GetDetector()=="PHOS"){
46a3cde6 404 // In case we use the MC reader, the input are TParticles,
405 // in this case use the corresponing method in PHOS Geometry to get the particle.
406 if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
407 {
408 Int_t mod =-1;
409 Double_t z = 0., x=0.;
410 TParticle* primary = 0x0;
411 AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
412 if(stack) {
413 primary = stack->Particle(particle->GetCaloLabel(0));
414 }
415 else {
280e6766 416 Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
46a3cde6 417 }
898c9d44 418
419 if(primary){
420 fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
421 }
422 else{
423 Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
424 }
46a3cde6 425 return mod;
426 }
427 // Input are ESDs or AODs, get the PHOS module number like this.
428 else{
2206e918 429 //FIXME
430 //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
431 //return GetModuleNumber(cluster);
432 //MEFIX
433 return -1;
46a3cde6 434 }
765d44e7 435 }//PHOS
436
437 return -1;
438}
439
cb5780f4 440//_____________________________________________________________________
c8fe2783 441Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
765d44e7 442{
280e6766 443 //Get the EMCAL/PHOS module number that corresponds to this cluster
765d44e7 444 TLorentzVector lv;
445 Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
9cbbc28b 446 if(!cluster){
447 if(fDebug > 1) printf("AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
448 return -1;
449 }
765d44e7 450 cluster->GetMomentum(lv,v);
451 Float_t phi = lv.Phi();
452 if(phi < 0) phi+=TMath::TwoPi();
453 Int_t absId = -1;
c8fe2783 454 if(cluster->IsEMCAL()){
765d44e7 455 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
456 if(fDebug > 2)
280e6766 457 printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
7cc7d3f8 458 lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
765d44e7 459 return fEMCALGeo->GetSuperModuleNumber(absId) ;
460 }//EMCAL
c8fe2783 461 else if(cluster->IsPHOS()) {
765d44e7 462 Int_t relId[4];
463 if ( cluster->GetNCells() > 0) {
464 absId = cluster->GetCellAbsId(0);
465 if(fDebug > 2)
280e6766 466 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
7cc7d3f8 467 lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
765d44e7 468 }
469 else return -1;
470
471 if ( absId >= 0) {
472 fPHOSGeo->AbsToRelNumbering(absId,relId);
473 if(fDebug > 2)
280e6766 474 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
765d44e7 475 return relId[0]-1;
476 }
477 else return -1;
478 }//PHOS
479
480 return -1;
481}
482
cb5780f4 483//___________________________________________________________________________________________________
484Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(const Int_t absId, const TString calo,
485 Int_t & icol, Int_t & irow, Int_t & iRCU) const
765d44e7 486{
487 //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
488 Int_t imod = -1;
489 if ( absId >= 0) {
490 if(calo=="EMCAL"){
491 Int_t iTower = -1, iIphi = -1, iIeta = -1;
492 fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
493 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
280e6766 494 if(imod < 0 || irow < 0 || icol < 0 ) {
495 Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
496 }
497
765d44e7 498 //RCU0
499 if (0<=irow&&irow<8) iRCU=0; // first cable row
500 else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
501 //second cable row
502 //RCU1
503 else if(8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
504 //second cable row
505 else if(16<=irow&&irow<24) iRCU=1; // third cable row
506
507 if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
508 if (iRCU<0) {
280e6766 509 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
765d44e7 510 }
511
512 return imod ;
513 }//EMCAL
514 else{//PHOS
515 Int_t relId[4];
516 fPHOSGeo->AbsToRelNumbering(absId,relId);
517 irow = relId[2];
518 icol = relId[3];
519 imod = relId[0]-1;
520 iRCU= (Int_t)(relId[2]-1)/16 ;
521 //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
522 if (iRCU >= 4) {
280e6766 523 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
765d44e7 524 }
525 return imod;
526 }//PHOS
527 }
528
529 return -1;
530}
531
cb5780f4 532//________________________________________
765d44e7 533void AliCalorimeterUtils::InitParameters()
534{
535 //Initialize the parameters of the analysis.
7cc7d3f8 536 fEMCALGeoName = "EMCAL_COMPLETEV1";
537 fPHOSGeoName = "PHOSgeo";
538 fEMCALGeoMatrixSet = kFALSE;
539 fPHOSGeoMatrixSet = kFALSE;
540 fRemoveBadChannels = kFALSE;
541 fNCellsFromPHOSBorder = 0;
a5fb4114 542
7cc7d3f8 543 // fMaskCellColumns = new Int_t[fNMaskCellColumns];
544 // fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
545 // fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
546 // fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
547 // fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
548 // fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
a5fb4114 549
765d44e7 550}
551
765d44e7 552
cb5780f4 553//_____________________________________________________
554void AliCalorimeterUtils::InitPHOSBadChannelStatusMap()
555{
765d44e7 556 //Init PHOS bad channels map
557 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
558 //In order to avoid rewriting the same histograms
559 Bool_t oldStatus = TH1::AddDirectoryStatus();
560 TH1::AddDirectory(kFALSE);
7cc7d3f8 561
78219bac 562 fPHOSBadChannelMap = new TObjArray(5);
c5693f62 563 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 564
765d44e7 565 fPHOSBadChannelMap->SetOwner(kTRUE);
566 fPHOSBadChannelMap->Compress();
567
568 //In order to avoid rewriting the same histograms
569 TH1::AddDirectory(oldStatus);
570}
571
cb5780f4 572//______________________________________________________
573void AliCalorimeterUtils::InitPHOSRecalibrationFactors()
574{
09e819c9 575 //Init EMCAL recalibration factors
576 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
577 //In order to avoid rewriting the same histograms
578 Bool_t oldStatus = TH1::AddDirectoryStatus();
579 TH1::AddDirectory(kFALSE);
7cc7d3f8 580
78219bac 581 fPHOSRecalibrationFactors = new TObjArray(5);
c5693f62 582 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 583 //Init the histograms with 1
584 for (Int_t m = 0; m < 5; m++) {
585 for (Int_t i = 0; i < 56; i++) {
586 for (Int_t j = 0; j < 64; j++) {
c5693f62 587 SetPHOSChannelRecalibrationFactor(m,j,i,1.);
09e819c9 588 }
589 }
590 }
591 fPHOSRecalibrationFactors->SetOwner(kTRUE);
592 fPHOSRecalibrationFactors->Compress();
593
594 //In order to avoid rewriting the same histograms
595 TH1::AddDirectory(oldStatus);
596}
597
598
cb5780f4 599//___________________________________________
765d44e7 600void AliCalorimeterUtils::InitEMCALGeometry()
601{
602 //Initialize EMCAL geometry if it did not exist previously
603 if (!fEMCALGeo){
a38a48f2 604 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
765d44e7 605 if(fDebug > 0){
606 printf("AliCalorimeterUtils::InitEMCALGeometry()");
607 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
608 printf("\n");
609 }
610 }
611}
612
cb5780f4 613//__________________________________________
765d44e7 614void AliCalorimeterUtils::InitPHOSGeometry()
615{
616 //Initialize PHOS geometry if it did not exist previously
617 if (!fPHOSGeo){
618 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
619 if(fDebug > 0){
620 printf("AliCalorimeterUtils::InitPHOSGeometry()");
621 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
622 printf("\n");
623 }
624 }
625}
626
cb5780f4 627//_________________________________________________________
765d44e7 628void AliCalorimeterUtils::Print(const Option_t * opt) const
629{
7cc7d3f8 630
765d44e7 631 //Print some relevant parameters set for the analysis
632 if(! opt)
633 return;
7cc7d3f8 634
765d44e7 635 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
636 printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
637 printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
7cc7d3f8 638 fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
247abff4 639 if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
09e819c9 640 printf("Recalibrate Clusters? %d\n",fRecalibration);
9584c261 641 printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
642 printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
f2ccb5b8 643 printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
7cc7d3f8 644
765d44e7 645 printf(" \n") ;
646}
647
cb5780f4 648//__________________________________________________________________
649Bool_t AliCalorimeterUtils::MaskFrameCluster(const Int_t iSM,
650 const Int_t ieta) const
651{
a5fb4114 652 //Check if cell is in one of the regions where we have significant amount
653 //of material in front. Only EMCAL
654
655 Int_t icol = ieta;
656 if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
657
658 if (fNMaskCellColumns && fMaskCellColumns) {
659 for (Int_t imask = 0; imask < fNMaskCellColumns; imask++) {
660 if(icol==fMaskCellColumns[imask]) return kTRUE;
661 }
662 }
663
664 return kFALSE;
665
666}
667
668
cb5780f4 669//__________________________________________________________________________
670Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster,
671 AliVCaloCells * cells)
672{
09e819c9 673 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
7cc7d3f8 674
898c9d44 675 //Initialize some used variables
676 Float_t energy = 0;
677 Int_t absId = -1;
678 Int_t icol = -1, irow = -1, iRCU = -1, module=1;
679 Float_t factor = 1, frac = 0;
680
681 if(cells) {
7cc7d3f8 682
683 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
684 UShort_t * index = cluster->GetCellsAbsId() ;
685 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
686 Int_t ncells = cluster->GetNCells();
687 TString calo = "EMCAL";
688 if(cluster->IsPHOS()) calo = "PHOS";
689
690 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
691 for(Int_t icell = 0; icell < ncells; icell++){
692 absId = index[icell];
693 frac = fraction[icell];
694 if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
695 module = GetModuleNumberCellIndexes(absId,calo,icol,irow,iRCU);
c5693f62 696 if(cluster->IsPHOS()) factor = GetPHOSChannelRecalibrationFactor (module,irow,icol);
7cc7d3f8 697 else factor = GetEMCALChannelRecalibrationFactor(module,icol,irow);
698 if(fDebug>2)
699 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, module %d, col %d, row %d, cell fraction %f, recalibration factor %f, cell energy %f\n",
700 calo.Data(),module,icol,irow,frac,factor,cells->GetCellAmplitude(absId));
701
702 energy += cells->GetCellAmplitude(absId)*factor*frac;
703 }
704
705 if(fDebug>1)
706 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
898c9d44 707
708 }// cells available
709 else{
710 Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
711 }
712
09e819c9 713 return energy;
09e819c9 714}
715
cb5780f4 716//________________________________________________________________________________
765d44e7 717void AliCalorimeterUtils::SetGeometryTransformationMatrices(AliVEvent* inputEvent)
718{
719 //Set the calorimeters transformation matrices
7cc7d3f8 720
765d44e7 721 //Get the EMCAL transformation geometry matrices from ESD
3b13c34c 722 if(!fEMCALGeoMatrixSet && fEMCALGeo){
723 if(fLoadEMCALMatrices){
c5693f62 724 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined EMCAL geometry matrices\n");
3b13c34c 725 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
726 if(fEMCALMatrix[mod]){
727 if(fDebug > 1)
728 fEMCALMatrix[mod]->Print();
729 fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;
730 }
731 }//SM loop
732 fEMCALGeoMatrixSet = kTRUE;//At least one, so good
733
734 }//Load matrices
735 else if (!gGeoManager) {
7cc7d3f8 736
3b13c34c 737 if(fDebug > 1)
738 printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load EMCAL misalignment matrices. \n");
739 if(!strcmp(inputEvent->GetName(),"AliESDEvent")) {
740 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
7cc7d3f8 741 //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
4b892846 742 if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod)) {
743 fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
3b13c34c 744 }
745 }// loop over super modules
c5693f62 746
3b13c34c 747 fEMCALGeoMatrixSet = kTRUE;//At least one, so good
748
749 }//ESD as input
750 else {
751 if(fDebug > 1)
752 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
753 }//AOD as input
754 }//Get matrix from data
04131edb 755 else if(gGeoManager){
756 fEMCALGeoMatrixSet = kTRUE;
757 }
765d44e7 758 }//EMCAL geo && no geoManager
c5693f62 759
765d44e7 760 //Get the PHOS transformation geometry matrices from ESD
3b13c34c 761 if(!fPHOSGeoMatrixSet && fPHOSGeo){
c5693f62 762
3b13c34c 763 if(fLoadPHOSMatrices){
c5693f62 764 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined PHOS geometry matrices\n");
765 for(Int_t mod = 0 ; mod < 5 ; mod++){
3b13c34c 766 if(fPHOSMatrix[mod]){
c5693f62 767 if(fDebug > 1 )
3b13c34c 768 fPHOSMatrix[mod]->Print();
dbf54f1e 769 fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;
3b13c34c 770 }
771 }//SM loop
772 fPHOSGeoMatrixSet = kTRUE;//At least one, so good
c5693f62 773
3b13c34c 774 }//Load matrices
775 else if (!gGeoManager) {
776 if(fDebug > 1)
777 printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
765d44e7 778 if(!strcmp(inputEvent->GetName(),"AliESDEvent")) {
dbf54f1e 779 for(Int_t mod = 0; mod < 5; mod++){
4b892846 780 if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod)) {
765d44e7 781 //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
4b892846 782 fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
765d44e7 783 }
3b13c34c 784 }// loop over modules
785 fPHOSGeoMatrixSet = kTRUE; //At least one so good
765d44e7 786 }//ESD as input
787 else {
788 if(fDebug > 1)
789 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
3b13c34c 790 }//AOD as input
791 }// get matrix from data
04131edb 792 else if(gGeoManager){
793 fPHOSGeoMatrixSet = kTRUE;
794 }
765d44e7 795 }//PHOS geo and geoManager was not set
c5693f62 796
765d44e7 797}
798
cb5780f4 799//__________________________________________________________________________________________
800void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
801{
9584c261 802
803 //Recalculate EMCAL cluster position
804
19db8f8c 805 fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
7cc7d3f8 806
9584c261 807}
cb5780f4 808
809//________________________________________________________________________________
810void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event,
811 TObjArray* clusterArray)
812{
813 //Recalculate track matching
814
815 if (fRecalculateMatching) {
816 fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo) ;
817 //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
818 //if(esdevent){
819 // fEMCALRecoUtils->SetClusterMatchedToTrack(esdevent) ;
820 // fEMCALRecoUtils->SetTracksMatchedToCluster(esdevent) ;
821 //}
822 }
823}