]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/PartCorrBase/AliCalorimeterUtils.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Allow position recaculation for EMCAL rotational and translation shifts, add asymmetr...
[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / 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 **************************************************************************/
15/* $Id: $ */
16
17//_________________________________________________________________________
18// Class utility for Calorimeter specific selection methods ///
19//
20//
21//
22//-- Author: Gustavo Conesa (LPSC-Grenoble)
23//////////////////////////////////////////////////////////////////////////////
24
25
26// --- ROOT system ---
27#include "TGeoManager.h"
28
29//---- ANALYSIS system ----
30#include "AliCalorimeterUtils.h"
c8fe2783 31#include "AliVEvent.h"
46a3cde6 32#include "AliMCEvent.h"
33#include "AliStack.h"
765d44e7 34#include "AliAODPWG4Particle.h"
c8fe2783 35#include "AliVCluster.h"
36#include "AliVCaloCells.h"
37#include "AliMixedEvent.h"
9584c261 38#include "AliEMCALRecoUtils.h"
765d44e7 39
40
41ClassImp(AliCalorimeterUtils)
42
43
44//____________________________________________________________________________
45 AliCalorimeterUtils::AliCalorimeterUtils() :
46 TObject(), fDebug(0),
c8fe2783 47 fEMCALGeoName("EMCAL_FIRSTYEAR"),fPHOSGeoName("PHOSgeo"),
765d44e7 48 fEMCALGeo(0x0), fPHOSGeo(0x0),
49 fEMCALGeoMatrixSet(kFALSE), fPHOSGeoMatrixSet(kFALSE),
50 fRemoveBadChannels(kFALSE),
78219bac 51 fEMCALBadChannelMap(0x0),fPHOSBadChannelMap(0x0),
765d44e7 52 fNCellsFromEMCALBorder(0), fNCellsFromPHOSBorder(0),
09e819c9 53 fNoEMCALBorderAtEta0(kFALSE),fRecalibration(kFALSE),
9584c261 54 fEMCALRecalibrationFactors(), fPHOSRecalibrationFactors(),
55 fEMCALRecoUtils(new AliEMCALRecoUtils),fRecalculatePosition(kFALSE),fCorrectELinearity(kFALSE)
56
765d44e7 57{
58 //Ctor
59
60 //Initialize parameters
61 InitParameters();
62}
765d44e7 63
64//_________________________________
65AliCalorimeterUtils::~AliCalorimeterUtils() {
66 //Dtor
67
4df35693 68 //if(fPHOSGeo) delete fPHOSGeo ;
765d44e7 69 if(fEMCALGeo) delete fEMCALGeo ;
70
71 if(fEMCALBadChannelMap) {
72 fEMCALBadChannelMap->Clear();
73 delete fEMCALBadChannelMap;
74 }
75 if(fPHOSBadChannelMap) {
76 fPHOSBadChannelMap->Clear();
77 delete fPHOSBadChannelMap;
78 }
79
09e819c9 80 if(fEMCALRecalibrationFactors) {
78219bac 81 fEMCALRecalibrationFactors->Clear();
82 delete fEMCALRecalibrationFactors;
09e819c9 83 }
84 if(fPHOSRecalibrationFactors) {
78219bac 85 fPHOSRecalibrationFactors->Clear();
86 delete fPHOSRecalibrationFactors;
09e819c9 87 }
88
9584c261 89 if(fEMCALRecoUtils) delete fEMCALRecoUtils ;
90
765d44e7 91}
92
93//_______________________________________________________________
c8fe2783 94Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster, AliVCaloCells* cells, AliVEvent * event, Int_t iev) const {
95
765d44e7 96 // Given the list of AbsId of the cluster, get the maximum cell and
97 // check if there are fNCellsFromBorder from the calorimeter border
98
c8fe2783 99 //If the distance to the border is 0 or negative just exit accept all clusters
100 if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
101 if(cells->GetType()==AliVCaloCells::kPHOSCell && fNCellsFromPHOSBorder <= 0 ) return kTRUE;
765d44e7 102
c8fe2783 103 Int_t absIdMax = -1;
765d44e7 104 Float_t ampMax = -1;
c8fe2783 105
106 AliMixedEvent * mixEvent = dynamic_cast<AliMixedEvent*> (event);
107 Int_t nMixedEvents = 0 ;
108 Int_t * cellsCumul = NULL ;
109 Int_t numberOfCells = 0 ;
110 if (mixEvent){
111 nMixedEvents = mixEvent->GetNumberOfEvents() ;
112 if (cells->GetType()==AliVCaloCells::kEMCALCell) {
113 cellsCumul = mixEvent->GetEMCALCellsCumul() ;
114 numberOfCells = mixEvent->GetNumberOfEMCALCells() ;
115 }
116
117 else if (cells->GetType()==AliVCaloCells::kPHOSCell) {
118 cellsCumul = mixEvent->GetPHOSCellsCumul() ;
119 numberOfCells = mixEvent->GetNumberOfPHOSCells() ;
120 }
898c9d44 121
122 if(cellsCumul){
123
124 Int_t startCell = cellsCumul[iev] ;
125 Int_t endCell = (iev+1 < nMixedEvents)?cellsCumul[iev+1]:numberOfCells;
c8fe2783 126 //Find cells with maximum amplitude
898c9d44 127 for(Int_t i = 0; i < cluster->GetNCells() ; i++){
128 Int_t absId = cluster->GetCellAbsId(i) ;
129 for (Int_t j = startCell; j < endCell ; j++) {
130 Short_t cellNumber;
131 Double_t amp ;
132 Double_t time;
133 cells->GetCell(j, cellNumber, amp, time) ;
134 if (absId == cellNumber) {
135 if(amp > ampMax){
136 ampMax = amp;
137 absIdMax = absId;
138 }
139 }
c8fe2783 140 }
898c9d44 141 }//loop on cluster cells
142 }// cells cumul available
143 else {
144 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - CellsCumul is NULL!!!\n");
145 abort();
c8fe2783 146 }
898c9d44 147 } else {//Normal SE Events
c8fe2783 148 for(Int_t i = 0; i < cluster->GetNCells() ; i++){
149 Int_t absId = cluster->GetCellAbsId(i) ;
150 Float_t amp = cells->GetCellAmplitude(absId);
151 if(amp > ampMax){
152 ampMax = amp;
153 absIdMax = absId;
154 }
155 }
156 }
765d44e7 157
158 if(fDebug > 1)
280e6766 159 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n",
765d44e7 160 absIdMax, ampMax, cluster->E());
161
162 if(absIdMax==-1) return kFALSE;
163
164 //Check if the cell is close to the borders:
165 Bool_t okrow = kFALSE;
166 Bool_t okcol = kFALSE;
167
c8fe2783 168 if(cells->GetType()==AliVCaloCells::kEMCALCell){
765d44e7 169
170 Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1;
171 fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta);
172 fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
280e6766 173 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
174 Fatal("CheckCellFidutialRegion","Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",iSM,ieta,iphi);
175 }
176
765d44e7 177 //Check rows/phi
178 if(iSM < 10){
179 if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
180 }
181 else{
182 if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
183 }
184
185 //Check collumns/eta
186 if(!fNoEMCALBorderAtEta0){
187 if(ieta > fNCellsFromEMCALBorder && ieta < 48-fNCellsFromEMCALBorder) okcol =kTRUE;
188 }
189 else{
190 if(iSM%2==0){
191 if(ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
192 }
193 else {
194 if(ieta < 48-fNCellsFromEMCALBorder) okcol = kTRUE;
195 }
196 }//eta 0 not checked
197 if(fDebug > 1)
198 {
280e6766 199 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
765d44e7 200 fNCellsFromEMCALBorder, ieta, iphi, iSM);
201 if (okcol && okrow ) printf(" YES \n");
202 else printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
203 }
204 }//EMCAL
c8fe2783 205 else if(cells->GetType()==AliVCaloCells::kPHOSCell){
765d44e7 206 Int_t relId[4];
207 Int_t irow = -1, icol = -1;
208 fPHOSGeo->AbsToRelNumbering(absIdMax,relId);
209 irow = relId[2];
210 icol = relId[3];
211 //imod = relId[0]-1;
212 if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE;
213 if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE;
214 if(fDebug > 1)
215 {
280e6766 216 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - PHOS Cluster in %d cells fiducial volume: icol %d, irow %d, Module %d?",
765d44e7 217 fNCellsFromPHOSBorder, icol, irow, relId[0]-1);
218 if (okcol && okrow ) printf(" YES \n");
219 else printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
220 }
221 }//PHOS
222
223 if (okcol && okrow) return kTRUE;
224 else return kFALSE;
225
226}
227
765d44e7 228//_________________________________________________________________________________________________________
229Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter,UShort_t* cellList, Int_t nCells){
230 // Check that in the cluster cells, there is no bad channel of those stored
231 // in fEMCALBadChannelMap or fPHOSBadChannelMap
232
233 if (!fRemoveBadChannels) return kFALSE;
36037088 234 //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
78219bac 235 if(calorimeter == "EMCAL" && !fEMCALBadChannelMap) return kFALSE;
236 if(calorimeter == "PHOS" && !fPHOSBadChannelMap) return kFALSE;
765d44e7 237
238 Int_t icol = -1;
239 Int_t irow = -1;
240 Int_t imod = -1;
241 for(Int_t iCell = 0; iCell<nCells; iCell++){
242
243 //Get the column and row
244 if(calorimeter == "EMCAL"){
245 Int_t iTower = -1, iIphi = -1, iIeta = -1;
246 fEMCALGeo->GetCellIndex(cellList[iCell],imod,iTower,iIphi,iIeta);
247 if(fEMCALBadChannelMap->GetEntries() <= imod) continue;
280e6766 248 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
249
250 if(imod < 0 || irow < 0 || icol < 0 ) {
251 Fatal("ClusterContainsBadChannel","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
252 }
253
765d44e7 254 if(GetEMCALChannelStatus(imod, icol, irow))return kTRUE;
255 }
256 else if(calorimeter=="PHOS"){
257 Int_t relId[4];
258 fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
259 irow = relId[2];
260 icol = relId[3];
261 imod = relId[0]-1;
262 if(fPHOSBadChannelMap->GetEntries() <= imod)continue;
263 if(GetPHOSChannelStatus(imod, icol, irow)) return kTRUE;
264 }
265 else return kFALSE;
266
267 }// cell cluster loop
268
269 return kFALSE;
270
271}
272
9584c261 273//____________________________________________________________________________________________________________________________________________________
274void AliCalorimeterUtils::CorrectClusterEnergy(AliVCluster *clus){
275 // Correct cluster energy non linearity
276 clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
277}
278
765d44e7 279//____________________________________________________________________________________________________________________________________________________
280Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
281{
282 //Get the EMCAL/PHOS module number that corresponds to this particle
283
284 Int_t absId = -1;
285 if(particle->GetDetector()=="EMCAL"){
286 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
287 if(fDebug > 2)
288 printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
289 particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
290 return fEMCALGeo->GetSuperModuleNumber(absId) ;
291 }//EMCAL
292 else if(particle->GetDetector()=="PHOS"){
46a3cde6 293 // In case we use the MC reader, the input are TParticles,
294 // in this case use the corresponing method in PHOS Geometry to get the particle.
295 if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
296 {
297 Int_t mod =-1;
298 Double_t z = 0., x=0.;
299 TParticle* primary = 0x0;
300 AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
301 if(stack) {
302 primary = stack->Particle(particle->GetCaloLabel(0));
303 }
304 else {
280e6766 305 Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
46a3cde6 306 }
898c9d44 307
308 if(primary){
309 fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
310 }
311 else{
312 Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
313 }
46a3cde6 314 return mod;
315 }
316 // Input are ESDs or AODs, get the PHOS module number like this.
317 else{
318 AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
319 return GetModuleNumber(cluster);
320 }
765d44e7 321 }//PHOS
322
323 return -1;
324}
325
326//____________________________________________________________________________________________________________________________________________________
c8fe2783 327Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
765d44e7 328{
280e6766 329 //Get the EMCAL/PHOS module number that corresponds to this cluster
765d44e7 330 TLorentzVector lv;
331 Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
332 cluster->GetMomentum(lv,v);
333 Float_t phi = lv.Phi();
334 if(phi < 0) phi+=TMath::TwoPi();
335 Int_t absId = -1;
c8fe2783 336 if(cluster->IsEMCAL()){
765d44e7 337 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
338 if(fDebug > 2)
280e6766 339 printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
765d44e7 340 lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
341 return fEMCALGeo->GetSuperModuleNumber(absId) ;
342 }//EMCAL
c8fe2783 343 else if(cluster->IsPHOS()) {
765d44e7 344 Int_t relId[4];
345 if ( cluster->GetNCells() > 0) {
346 absId = cluster->GetCellAbsId(0);
347 if(fDebug > 2)
280e6766 348 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
765d44e7 349 lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
350 }
351 else return -1;
352
353 if ( absId >= 0) {
354 fPHOSGeo->AbsToRelNumbering(absId,relId);
355 if(fDebug > 2)
280e6766 356 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
765d44e7 357 return relId[0]-1;
358 }
359 else return -1;
360 }//PHOS
361
362 return -1;
363}
364
765d44e7 365//_____________________________________________________________________________________________________________
366Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(const Int_t absId, const TString calo, Int_t & icol, Int_t & irow, Int_t & iRCU) const
367{
368 //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
369 Int_t imod = -1;
370 if ( absId >= 0) {
371 if(calo=="EMCAL"){
372 Int_t iTower = -1, iIphi = -1, iIeta = -1;
373 fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
374 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
280e6766 375 if(imod < 0 || irow < 0 || icol < 0 ) {
376 Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
377 }
378
765d44e7 379 //RCU0
380 if (0<=irow&&irow<8) iRCU=0; // first cable row
381 else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
382 //second cable row
383 //RCU1
384 else if(8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
385 //second cable row
386 else if(16<=irow&&irow<24) iRCU=1; // third cable row
387
388 if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
389 if (iRCU<0) {
280e6766 390 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
765d44e7 391 }
392
393 return imod ;
394 }//EMCAL
395 else{//PHOS
396 Int_t relId[4];
397 fPHOSGeo->AbsToRelNumbering(absId,relId);
398 irow = relId[2];
399 icol = relId[3];
400 imod = relId[0]-1;
401 iRCU= (Int_t)(relId[2]-1)/16 ;
402 //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
403 if (iRCU >= 4) {
280e6766 404 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
765d44e7 405 }
406 return imod;
407 }//PHOS
408 }
409
410 return -1;
411}
412
413//_______________________________________________________________
414//void AliCalorimeterUtils::Init()
415//{
416// //Init reader. Method to be called in AliAnaPartCorrMaker
417//
418// fEMCALBadChannelMap->SetName(Form("EMCALBadMap_%s",fTaskName.Data()));
419// fPHOSBadChannelMap->SetName(Form("PHOSBadMap_%s",fTaskName.Data()));
420//}
421
422//_______________________________________________________________
423void AliCalorimeterUtils::InitParameters()
424{
425 //Initialize the parameters of the analysis.
c8fe2783 426 fEMCALGeoName = "EMCAL_FIRSTYEAR";
765d44e7 427 fPHOSGeoName = "PHOSgeo";
428
429 if(gGeoManager) {// geoManager was set
430 if(fDebug > 2)printf("AliCalorimeterUtils::InitParameters() - Geometry manager available\n");
431 fEMCALGeoMatrixSet = kTRUE;
432 fPHOSGeoMatrixSet = kTRUE;
433 }
434 else{
435 fEMCALGeoMatrixSet = kFALSE;
436 fPHOSGeoMatrixSet = kFALSE;
437 }
438
439 fRemoveBadChannels = kFALSE;
440
441 fNCellsFromEMCALBorder = 0;
442 fNCellsFromPHOSBorder = 0;
443 fNoEMCALBorderAtEta0 = kFALSE;
444}
445
446//________________________________________________________________
447void AliCalorimeterUtils::InitEMCALBadChannelStatusMap(){
448 //Init EMCAL bad channels map
449 if(fDebug > 0 )printf("AliCalorimeterUtils::InitEMCALBadChannelStatusMap()\n");
450 //In order to avoid rewriting the same histograms
451 Bool_t oldStatus = TH1::AddDirectoryStatus();
452 TH1::AddDirectory(kFALSE);
78219bac 453
454 fEMCALBadChannelMap = new TObjArray(12);
455 //TH2F * hTemp = new TH2I("EMCALBadChannelMap","EMCAL SuperModule bad channel map", 48, 0, 48, 24, 0, 24);
456 for (int i = 0; i < 12; i++) {
457 fEMCALBadChannelMap->Add(new TH2I(Form("EMCALBadChannelMap_Mod%d",i),Form("EMCALBadChannelMap_Mod%d",i), 48, 0, 48, 24, 0, 24));
458 //fEMCALBadChannelMap->Add((TH2I*)hTemp->Clone(Form("EMCALBadChannelMap_Mod%d",i)));
459 }
460
461 //delete hTemp;
765d44e7 462
463 fEMCALBadChannelMap->SetOwner(kTRUE);
464 fEMCALBadChannelMap->Compress();
78219bac 465
765d44e7 466 //In order to avoid rewriting the same histograms
467 TH1::AddDirectory(oldStatus);
468}
469
470//________________________________________________________________
471void AliCalorimeterUtils::InitPHOSBadChannelStatusMap(){
472 //Init PHOS bad channels map
473 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
474 //In order to avoid rewriting the same histograms
475 Bool_t oldStatus = TH1::AddDirectoryStatus();
476 TH1::AddDirectory(kFALSE);
78219bac 477
478 fPHOSBadChannelMap = new TObjArray(5);
765d44e7 479 for (int i = 0; i < 5; i++)fPHOSBadChannelMap->Add(new TH2I(Form("PHOSBadChannelMap_Mod%d",i),Form("PHOSBadChannelMap_Mod%d",i), 56, 0, 56, 64, 0, 64));
480
481 fPHOSBadChannelMap->SetOwner(kTRUE);
482 fPHOSBadChannelMap->Compress();
483
484 //In order to avoid rewriting the same histograms
485 TH1::AddDirectory(oldStatus);
486}
487
09e819c9 488//________________________________________________________________
489void AliCalorimeterUtils::InitEMCALRecalibrationFactors(){
490 //Init EMCAL recalibration factors
491 if(fDebug > 0 )printf("AliCalorimeterUtils::InitEMCALRecalibrationFactors()\n");
492 //In order to avoid rewriting the same histograms
493 Bool_t oldStatus = TH1::AddDirectoryStatus();
494 TH1::AddDirectory(kFALSE);
78219bac 495
496 fEMCALRecalibrationFactors = new TObjArray(12);
09e819c9 497 for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
498 //Init the histograms with 1
499 for (Int_t sm = 0; sm < 12; sm++) {
500 for (Int_t i = 0; i < 48; i++) {
501 for (Int_t j = 0; j < 24; j++) {
502 SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
503 }
504 }
505 }
506 fEMCALRecalibrationFactors->SetOwner(kTRUE);
507 fEMCALRecalibrationFactors->Compress();
508
509 //In order to avoid rewriting the same histograms
510 TH1::AddDirectory(oldStatus);
511}
512
513//________________________________________________________________
514void AliCalorimeterUtils::InitPHOSRecalibrationFactors(){
515 //Init EMCAL recalibration factors
516 if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
517 //In order to avoid rewriting the same histograms
518 Bool_t oldStatus = TH1::AddDirectoryStatus();
519 TH1::AddDirectory(kFALSE);
78219bac 520
521 fPHOSRecalibrationFactors = new TObjArray(5);
09e819c9 522 for (int i = 0; i < 5; i++)fPHOSRecalibrationFactors->Add(new TH2F(Form("PHOSRecalFactors_Mod%d",i),Form("PHOSRecalFactors_Mod%d",i), 56, 0, 56, 64, 0, 64));
523 //Init the histograms with 1
524 for (Int_t m = 0; m < 5; m++) {
525 for (Int_t i = 0; i < 56; i++) {
526 for (Int_t j = 0; j < 64; j++) {
527 SetPHOSChannelRecalibrationFactor(m,i,j,1.);
528 }
529 }
530 }
531 fPHOSRecalibrationFactors->SetOwner(kTRUE);
532 fPHOSRecalibrationFactors->Compress();
533
534 //In order to avoid rewriting the same histograms
535 TH1::AddDirectory(oldStatus);
536}
537
538
765d44e7 539//________________________________________________________________
540void AliCalorimeterUtils::InitEMCALGeometry()
541{
542 //Initialize EMCAL geometry if it did not exist previously
543 if (!fEMCALGeo){
544 fEMCALGeo = new AliEMCALGeoUtils(fEMCALGeoName);
545 if(fDebug > 0){
546 printf("AliCalorimeterUtils::InitEMCALGeometry()");
547 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
548 printf("\n");
549 }
550 }
551}
552
553//________________________________________________________________
554void AliCalorimeterUtils::InitPHOSGeometry()
555{
556 //Initialize PHOS geometry if it did not exist previously
557 if (!fPHOSGeo){
558 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
559 if(fDebug > 0){
560 printf("AliCalorimeterUtils::InitPHOSGeometry()");
561 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
562 printf("\n");
563 }
564 }
565}
566
567//________________________________________________________________
568void AliCalorimeterUtils::Print(const Option_t * opt) const
569{
570
571 //Print some relevant parameters set for the analysis
572 if(! opt)
573 return;
574
575 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
576 printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
577 printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
578 fNCellsFromEMCALBorder, fNCellsFromPHOSBorder);
579 if(fNoEMCALBorderAtEta0) printf("Do not remove EMCAL clusters at Eta = 0\n");
09e819c9 580 printf("Recalibrate Clusters? %d\n",fRecalibration);
9584c261 581 printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
582 printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
765d44e7 583
584 printf(" \n") ;
585}
586
09e819c9 587//________________________________________________________________
c8fe2783 588Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster, AliVCaloCells * cells){
09e819c9 589 // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
09e819c9 590
898c9d44 591 //Initialize some used variables
592 Float_t energy = 0;
593 Int_t absId = -1;
594 Int_t icol = -1, irow = -1, iRCU = -1, module=1;
595 Float_t factor = 1, frac = 0;
596
597 if(cells) {
598
09e819c9 599 //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
600 UShort_t * index = cluster->GetCellsAbsId() ;
601 Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
602 Int_t ncells = cluster->GetNCells();
603 TString calo = "EMCAL";
604 if(cluster->IsPHOS()) calo = "PHOS";
09e819c9 605
606 //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
607 for(Int_t icell = 0; icell < ncells; icell++){
608 absId = index[icell];
609 frac = fraction[icell];
610 if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
611 module = GetModuleNumberCellIndexes(absId,calo,icol,irow,iRCU);
612 if(cluster->IsPHOS()) factor = GetPHOSChannelRecalibrationFactor (module,icol,irow);
613 else factor = GetEMCALChannelRecalibrationFactor(module,icol,irow);
614 if(fDebug>2)
280e6766 615 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, module %d, col %d, row %d, cell fraction %f, recalibration factor %f, cell energy %f\n",
09e819c9 616 calo.Data(),module,icol,irow,frac,factor,cells->GetCellAmplitude(absId));
617
618 energy += cells->GetCellAmplitude(absId)*factor*frac;
619 }
620
621 if(fDebug>1)
280e6766 622 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
898c9d44 623
624 }// cells available
625 else{
626 Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
627 }
628
09e819c9 629 return energy;
09e819c9 630}
631
765d44e7 632//________________________________________________________________
633void AliCalorimeterUtils::SetGeometryTransformationMatrices(AliVEvent* inputEvent)
634{
635 //Set the calorimeters transformation matrices
636
637 //Get the EMCAL transformation geometry matrices from ESD
36037088 638 if (!gGeoManager && fEMCALGeo && !fEMCALGeoMatrixSet) {
765d44e7 639 if(fDebug > 1)
640 printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load EMCAL misalignment matrices. \n");
641 if(!strcmp(inputEvent->GetName(),"AliESDEvent")) {
642 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
c8fe2783 643 if(inputEvent->GetEMCALMatrix(mod)) {
765d44e7 644 //printf("EMCAL: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
c8fe2783 645 fEMCALGeo->SetMisalMatrix(inputEvent->GetEMCALMatrix(mod),mod) ;
765d44e7 646 fEMCALGeoMatrixSet = kTRUE;//At least one, so good
647 }
648 }// loop over super modules
649 }//ESD as input
650 else {
651 if(fDebug > 1)
652 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
653 }//AOD as input
654 }//EMCAL geo && no geoManager
655
656 //Get the PHOS transformation geometry matrices from ESD
657 if (!gGeoManager && fPHOSGeo && !fPHOSGeoMatrixSet) {
658 if(fDebug > 1)
659 printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
660 if(!strcmp(inputEvent->GetName(),"AliESDEvent")) {
661 for(Int_t mod=0; mod < 5; mod++){
c8fe2783 662 if(inputEvent->GetPHOSMatrix(mod)) {
765d44e7 663 //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
c8fe2783 664 fPHOSGeo->SetMisalMatrix(inputEvent->GetPHOSMatrix(mod),mod) ;
765d44e7 665 fPHOSGeoMatrixSet = kTRUE; //At least one so good
666 }
667 }// loop over modules
668 }//ESD as input
669 else {
670 if(fDebug > 1)
671 printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
672 }//AOD as input
673 }//PHOS geo and geoManager was not set
674
675}
676
9584c261 677//________________________________________________________________
678void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu){
679
680 //Recalculate EMCAL cluster position
681
682 Double_t eMax = -1.;
683 Double_t eCell = -1.;
684 Float_t fraction = 1.;
685 Int_t cellAbsId = -1;
686 Float_t recalFactor = 1.;
687
688 Int_t maxId = -1;
689 Int_t imod = -1, iphi = -1, ieta =-1;
690 Int_t iTower = -1, iIphi = -1, iIeta =-1;
691
692 Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
693 Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
694 Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
695 Int_t startingSM = -1;
696
697 for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
698 cellAbsId = clu->GetCellAbsId(iDig);
699 fraction = clu->GetCellAmplitudeFraction(iDig);
700 if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
701 fEMCALGeo->GetCellIndex(cellAbsId,imod,iTower,iIphi,iIeta);
702 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
703 if (iDig==0) startingSM = imod;
704 else if(imod != startingSM) areInSameSM = kFALSE;
705
706 if(IsRecalibrationOn()) {
707 recalFactor = GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
708 }
709 eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
710
711 weight = TMath::Log(eCell/clEnergy) + 4;
712 if(weight < 0) weight = 0;
713 totalWeight += weight;
714 weightedCol += ieta*weight;
715 weightedRow += iphi*weight;
716
717 //printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
718
719 if(eCell > eMax) {
720 eMax = eCell;
721 maxId = cellAbsId;
722 //printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
723 }
724 }// cell loop
725
726 //Get from the absid the supermodule, tower and eta/phi numbers
727 fEMCALGeo->GetCellIndex(maxId,imod,iTower,iIphi,iIeta);
728 //Gives SuperModule and Tower numbers
729 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,
730 iIphi, iIeta,iphi,ieta);
731
732 Float_t xyzNew[3];
733 if(areInSameSM == kTRUE) {
734 //printf("In Same SM\n");
735 weightedCol = weightedCol/totalWeight;
736 weightedRow = weightedRow/totalWeight;
737
738 //Float_t *xyzNew = RecalculatePosition(weightedRow, weightedCol, clEnergy, 0, iSupMod); //1 = electrons, 0 photons
739 fEMCALGeo->RecalculateTowerPosition(weightedRow, weightedCol, imod, clEnergy, 0, //1 = electrons, 0 photons
9fdaff9a 740 fEMCALRecoUtils->GetMisalTransShiftArray(), fEMCALRecoUtils->GetMisalRotShiftArray(), xyzNew);
9584c261 741 }
742 else {
743 //printf("In Different SM\n");
744 //Float_t *xyzNew = RecalculatePosition(iphi, ieta, clEnergy, 0, iSupMod); //1 = electrons, 0 photons
745 fEMCALGeo->RecalculateTowerPosition(iphi, ieta, imod, clEnergy, 0, //1 = electrons, 0 photons
9fdaff9a 746 fEMCALRecoUtils->GetMisalTransShiftArray(), fEMCALRecoUtils->GetMisalRotShiftArray(), xyzNew);
9584c261 747
748 }
749
750 clu->SetPosition(xyzNew);
751
752 //printf("\t Max : cell %d, iSupMod %d, ieta %d, iphi %d \n",maxId,iSupMod, ieta,iphi);
753
754}