588d08d8ec94fb5c2d2d9b69e9850d06850bc510
[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliCalorimeterUtils.cxx
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"
31 #include "AliVEvent.h"
32 #include "AliMCEvent.h"
33 #include "AliStack.h"
34 #include "AliAODPWG4Particle.h"
35 #include "AliVCluster.h"
36 #include "AliVCaloCells.h"
37 #include "AliMixedEvent.h"
38 #include "AliEMCALRecoUtils.h"
39
40
41 ClassImp(AliCalorimeterUtils)
42   
43   
44 //____________________________________________________________________________
45   AliCalorimeterUtils::AliCalorimeterUtils() : 
46     TObject(), fDebug(0), 
47     fEMCALGeoName("EMCAL_FIRSTYEAR"),fPHOSGeoName("PHOSgeo"), 
48     fEMCALGeo(0x0), fPHOSGeo(0x0), 
49     fEMCALGeoMatrixSet(kFALSE), fPHOSGeoMatrixSet(kFALSE), 
50     fRemoveBadChannels(kFALSE),
51     fEMCALBadChannelMap(0x0),fPHOSBadChannelMap(0x0), 
52     fNCellsFromEMCALBorder(0), fNCellsFromPHOSBorder(0), 
53     fNoEMCALBorderAtEta0(kFALSE),fRecalibration(kFALSE),
54     fEMCALRecalibrationFactors(), fPHOSRecalibrationFactors(),
55     fEMCALRecoUtils(new AliEMCALRecoUtils),fRecalculatePosition(kFALSE),fCorrectELinearity(kFALSE)
56
57 {
58   //Ctor
59   
60   //Initialize parameters
61   InitParameters();
62 }
63
64 //_________________________________
65 AliCalorimeterUtils::~AliCalorimeterUtils() {
66   //Dtor
67   
68   //if(fPHOSGeo)  delete fPHOSGeo  ;
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         
80   if(fEMCALRecalibrationFactors) { 
81     fEMCALRecalibrationFactors->Clear();
82     delete  fEMCALRecalibrationFactors;
83   }
84   if(fPHOSRecalibrationFactors) { 
85     fPHOSRecalibrationFactors->Clear();
86     delete  fPHOSRecalibrationFactors;
87   }
88         
89   if(fEMCALRecoUtils) delete fEMCALRecoUtils ;
90   
91 }
92
93 //_______________________________________________________________
94 Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster, AliVCaloCells* cells, AliVEvent * event, Int_t iev) const {
95
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         
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;
102
103   Int_t absIdMax        = -1;
104         Float_t ampMax  = -1;
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     } 
121     
122     if(cellsCumul){
123       
124       Int_t startCell = cellsCumul[iev] ; 
125       Int_t endCell   = (iev+1 < nMixedEvents)?cellsCumul[iev+1]:numberOfCells;
126       //Find cells with maximum amplitude
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           }
140         }
141       }//loop on cluster cells
142     }// cells cumul available
143     else {
144       printf("AliCalorimeterUtils::CheckCellFiducialRegion() - CellsCumul is NULL!!!\n");
145       abort();
146     }
147   } else {//Normal SE Events
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   }
157         
158         if(fDebug > 1)
159                 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n", 
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
168         if(cells->GetType()==AliVCaloCells::kEMCALCell){
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);
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       
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                 {
199                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
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
205         else if(cells->GetType()==AliVCaloCells::kPHOSCell){
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                 {
216                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - PHOS Cluster in %d cells fiducial volume: icol %d, irow %d, Module %d?",
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
228 //_________________________________________________________________________________________________________
229 Bool_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;
234         //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
235         if(calorimeter == "EMCAL" && !fEMCALBadChannelMap) return kFALSE;
236         if(calorimeter == "PHOS"  && !fPHOSBadChannelMap)  return kFALSE;
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;
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       
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
273 //____________________________________________________________________________________________________________________________________________________
274 void AliCalorimeterUtils::CorrectClusterEnergy(AliVCluster *clus){
275   // Correct cluster energy non linearity
276   clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
277 }
278
279 //____________________________________________________________________________________________________________________________________________________
280 Int_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"){
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 {
305         Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
306       }
307       
308       if(primary){
309         fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
310       }
311       else{
312         Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
313       }
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     }
321         }//PHOS
322         
323         return -1;
324 }
325
326 //____________________________________________________________________________________________________________________________________________________
327 Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
328 {
329         //Get the EMCAL/PHOS module number that corresponds to this cluster
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;
336         if(cluster->IsEMCAL()){
337                 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
338                 if(fDebug > 2) 
339                   printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
340                          lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
341                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
342         }//EMCAL
343         else if(cluster->IsPHOS()) {
344                 Int_t    relId[4];
345                 if ( cluster->GetNCells() > 0) {
346                         absId = cluster->GetCellAbsId(0);
347                         if(fDebug > 2) 
348                                 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
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) 
356                           printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
357                         return relId[0]-1;
358                 }
359                 else return -1;
360         }//PHOS
361         
362         return -1;
363 }
364
365 //_____________________________________________________________________________________________________________
366 Int_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);
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       
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) {
390                                 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
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) {
404                                 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
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 //_______________________________________________________________
423 void AliCalorimeterUtils::InitParameters()
424 {
425   //Initialize the parameters of the analysis.
426   fEMCALGeoName = "EMCAL_FIRSTYEAR";
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 //________________________________________________________________
447 void 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);
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;
462
463   fEMCALBadChannelMap->SetOwner(kTRUE);
464   fEMCALBadChannelMap->Compress();
465
466   //In order to avoid rewriting the same histograms
467   TH1::AddDirectory(oldStatus);         
468 }
469
470 //________________________________________________________________
471 void 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);
477
478   fPHOSBadChannelMap = new TObjArray(5);        
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
488 //________________________________________________________________
489 void 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);
495
496         fEMCALRecalibrationFactors = new TObjArray(12);
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 //________________________________________________________________
514 void 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);
520
521         fPHOSRecalibrationFactors = new TObjArray(5);
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
539 //________________________________________________________________
540 void 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 //________________________________________________________________
554 void 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 //________________________________________________________________
568 void 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");
580   printf("Recalibrate Clusters? %d\n",fRecalibration);
581   printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
582   printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
583
584   printf("    \n") ;
585
586
587 //________________________________________________________________
588 Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster, AliVCaloCells * cells){
589         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
590
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         
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";
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)
615                         printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, module %d, col %d, row %d, cell fraction %f, recalibration factor %f, cell energy %f\n", 
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)
622                 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
623     
624         }// cells available
625   else{
626     Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
627   }
628   
629         return energy;
630 }
631
632 //________________________________________________________________
633 void AliCalorimeterUtils::SetGeometryTransformationMatrices(AliVEvent* inputEvent) 
634 {
635   //Set the calorimeters transformation matrices
636         
637   //Get the EMCAL transformation geometry matrices from ESD 
638   if (!gGeoManager && fEMCALGeo && !fEMCALGeoMatrixSet) { 
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++){ 
643                                 if(inputEvent->GetEMCALMatrix(mod)) {
644                                         //printf("EMCAL: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
645                                         fEMCALGeo->SetMisalMatrix(inputEvent->GetEMCALMatrix(mod),mod) ;
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++){ 
662                                         if(inputEvent->GetPHOSMatrix(mod)) {
663                                                 //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
664                                                 fPHOSGeo->SetMisalMatrix(inputEvent->GetPHOSMatrix(mod),mod) ;
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
677 //________________________________________________________________
678 void 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
740                                         fEMCALRecoUtils->GetMisalShiftArray(), xyzNew);
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
746                                         fEMCALRecoUtils->GetMisalShiftArray(), xyzNew);
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 }