]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWG4/PartCorrBase/AliCalorimeterUtils.cxx
Comparison between different detectors in linear scale and in different canvases...
[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 "AliESDEvent.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
39 ClassImp(AliCalorimeterUtils)
40   
41   
42 //____________________________________________________________________________
43   AliCalorimeterUtils::AliCalorimeterUtils() : 
44     TObject(), fDebug(0), 
45     fEMCALGeoName("EMCAL_FIRSTYEARV1"),fPHOSGeoName("PHOSgeo"), 
46     fEMCALGeo(0x0), fPHOSGeo(0x0), 
47     fEMCALGeoMatrixSet(kFALSE), fPHOSGeoMatrixSet(kFALSE), 
48     fLoadEMCALMatrices(kFALSE), fLoadPHOSMatrices(kFALSE),
49     fRemoveBadChannels(kFALSE),fPHOSBadChannelMap(0x0), 
50     fNCellsFromPHOSBorder(0), fRecalibration(kFALSE),
51     fPHOSRecalibrationFactors(),
52     fEMCALRecoUtils(new AliEMCALRecoUtils),
53     fRecalculatePosition(kFALSE),fCorrectELinearity(kFALSE),
54     fRecalculateMatching(kFALSE),fCutR(20), fCutZ(20)
55 {
56   //Ctor
57   
58   //Initialize parameters
59   InitParameters();
60   for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
61   for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix[i]  = 0 ;
62
63 }
64
65 //_________________________________
66 AliCalorimeterUtils::~AliCalorimeterUtils() {
67   //Dtor
68   
69   //if(fPHOSGeo)  delete fPHOSGeo  ;
70   if(fEMCALGeo) delete fEMCALGeo ;
71
72   if(fPHOSBadChannelMap) { 
73     fPHOSBadChannelMap->Clear();
74     delete  fPHOSBadChannelMap;
75   }
76         
77   if(fPHOSRecalibrationFactors) { 
78     fPHOSRecalibrationFactors->Clear();
79     delete  fPHOSRecalibrationFactors;
80   }
81         
82   if(fEMCALRecoUtils) delete fEMCALRecoUtils ;
83   
84 }
85
86 //_______________________________________________________________
87 Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster, AliVCaloCells* cells, AliVEvent * event, Int_t iev) const {
88
89         // Given the list of AbsId of the cluster, get the maximum cell and 
90         // check if there are fNCellsFromBorder from the calorimeter border
91         
92     //If the distance to the border is 0 or negative just exit accept all clusters
93         if(cells->GetType()==AliVCaloCells::kEMCALCell && fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder() <= 0 ) return kTRUE;
94         if(cells->GetType()==AliVCaloCells::kPHOSCell  && fNCellsFromPHOSBorder  <= 0 ) return kTRUE;
95
96   Int_t absIdMax        = -1;
97         Float_t ampMax  = -1;
98         
99   AliMixedEvent * mixEvent = dynamic_cast<AliMixedEvent*> (event);
100   Int_t nMixedEvents = 0 ; 
101   Int_t * cellsCumul = NULL ;
102   Int_t numberOfCells = 0 ;  
103   if (mixEvent){
104     nMixedEvents = mixEvent->GetNumberOfEvents() ; 
105     if (cells->GetType()==AliVCaloCells::kEMCALCell) {
106       cellsCumul =  mixEvent->GetEMCALCellsCumul() ; 
107       numberOfCells = mixEvent->GetNumberOfEMCALCells() ;
108     } 
109       
110     else if (cells->GetType()==AliVCaloCells::kPHOSCell) {
111       cellsCumul =  mixEvent->GetPHOSCellsCumul() ; 
112       numberOfCells = mixEvent->GetNumberOfPHOSCells() ;
113     } 
114     
115     if(cellsCumul){
116       
117       Int_t startCell = cellsCumul[iev] ; 
118       Int_t endCell   = (iev+1 < nMixedEvents)?cellsCumul[iev+1]:numberOfCells;
119       //Find cells with maximum amplitude
120       for(Int_t i = 0; i < cluster->GetNCells() ; i++){
121         Int_t absId = cluster->GetCellAbsId(i) ;
122         for (Int_t j = startCell; j < endCell ;  j++) {
123           Short_t cellNumber; 
124           Double_t amp ; 
125           Double_t time; 
126           cells->GetCell(j, cellNumber, amp, time) ; 
127           if (absId == cellNumber) {
128             if(amp > ampMax){
129               ampMax   = amp;
130               absIdMax = absId;
131             }        
132           }
133         }
134       }//loop on cluster cells
135     }// cells cumul available
136     else {
137       printf("AliCalorimeterUtils::CheckCellFiducialRegion() - CellsCumul is NULL!!!\n");
138       abort();
139     }
140   } else {//Normal SE Events
141     for(Int_t i = 0; i < cluster->GetNCells() ; i++){
142       Int_t absId = cluster->GetCellAbsId(i) ;
143       Float_t amp       = cells->GetCellAmplitude(absId);
144       if(amp > ampMax){
145         ampMax   = amp;
146         absIdMax = absId;
147       }
148     }
149   }
150         
151         if(fDebug > 1)
152                 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n", 
153                    absIdMax, ampMax, cluster->E());
154         
155         if(absIdMax==-1) return kFALSE;
156         
157         //Check if the cell is close to the borders:
158         Bool_t okrow = kFALSE;
159         Bool_t okcol = kFALSE;
160
161         if(cells->GetType()==AliVCaloCells::kEMCALCell){
162         
163                 Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1; 
164                 fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta); 
165                 fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
166                 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
167       Fatal("CheckCellFidutialRegion","Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",iSM,ieta,iphi);
168     }
169       
170                 //Check rows/phi
171     Int_t nborder = fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder();
172                 if(iSM < 10){
173                         if(iphi >= nborder && iphi < 24-nborder) okrow =kTRUE; 
174             }
175                 else{
176                         if(iphi >= nborder && iphi < 12-nborder) okrow =kTRUE; 
177                 }
178                 
179                 //Check columns/eta
180                 if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0()){
181                         if(ieta  > nborder && ieta < 48-nborder) okcol =kTRUE; 
182                 }
183                 else{
184                         if(iSM%2==0){
185                                 if(ieta >= nborder)     okcol = kTRUE;  
186                         }
187                         else {
188                                 if(ieta <  48-nborder)  okcol = kTRUE;  
189                         }
190                 }//eta 0 not checked
191                 if(fDebug > 1)
192                 {
193                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
194                                    nborder, ieta, iphi, iSM);
195                         if (okcol && okrow ) printf(" YES \n");
196                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
197                 }
198         }//EMCAL
199         else if(cells->GetType()==AliVCaloCells::kPHOSCell){
200                 Int_t relId[4];
201                 Int_t irow = -1, icol = -1;
202                 fPHOSGeo->AbsToRelNumbering(absIdMax,relId);
203                 irow = relId[2];
204                 icol = relId[3];
205                 //imod = relId[0]-1;
206                 if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE; 
207                 if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE; 
208                 if(fDebug > 1)
209                 {
210                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - PHOS Cluster in %d cells fiducial volume: icol %d, irow %d, Module %d?",
211                                    fNCellsFromPHOSBorder, icol, irow, relId[0]-1);
212                         if (okcol && okrow ) printf(" YES \n");
213                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
214                 }
215         }//PHOS
216         
217         if (okcol && okrow) return kTRUE; 
218         else                return kFALSE;
219         
220 }       
221
222 //_________________________________________________________________________________________________________
223 Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter,UShort_t* cellList, Int_t nCells){
224         // Check that in the cluster cells, there is no bad channel of those stored 
225         // in fEMCALBadChannelMap or fPHOSBadChannelMap
226         
227         if (!fRemoveBadChannels) return kFALSE;
228         //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
229         if(calorimeter == "EMCAL" && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
230         if(calorimeter == "PHOS"  && !fPHOSBadChannelMap)  return kFALSE;
231
232         Int_t icol = -1;
233         Int_t irow = -1;
234         Int_t imod = -1;
235         for(Int_t iCell = 0; iCell<nCells; iCell++){
236         
237                 //Get the column and row
238                 if(calorimeter == "EMCAL"){
239       return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
240                 }
241                 else if(calorimeter=="PHOS"){
242                         Int_t    relId[4];
243                         fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
244                         irow = relId[2];
245                         icol = relId[3];
246                         imod = relId[0]-1;
247                         if(fPHOSBadChannelMap->GetEntries() <= imod)continue;
248                         if(GetPHOSChannelStatus(imod, icol, irow)) return kTRUE;
249                 }
250                 else return kFALSE;
251                 
252         }// cell cluster loop
253
254         return kFALSE;
255
256 }
257
258 //____________________________________________________________________________________________________________________________________________________
259 void AliCalorimeterUtils::CorrectClusterEnergy(AliVCluster *clus){
260   // Correct cluster energy non linearity
261   clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
262 }
263
264 //____________________________________________________________________________________________________________________________________________________
265 Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
266 {
267         //Get the EMCAL/PHOS module number that corresponds to this particle
268         
269         Int_t absId = -1;
270         if(particle->GetDetector()=="EMCAL"){
271                 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
272                 if(fDebug > 2) 
273                   printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
274                          particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
275                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
276         }//EMCAL
277         else if(particle->GetDetector()=="PHOS"){
278     // In case we use the MC reader, the input are TParticles, 
279     // in this case use the corresponing method in PHOS Geometry to get the particle.
280     if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
281     {
282       Int_t mod =-1;
283       Double_t z = 0., x=0.;
284       TParticle* primary = 0x0;
285       AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
286       if(stack) {
287         primary = stack->Particle(particle->GetCaloLabel(0));
288       }
289       else {
290         Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
291       }
292       
293       if(primary){
294         fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
295       }
296       else{
297         Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
298       }
299       return mod;
300     }
301     // Input are ESDs or AODs, get the PHOS module number like this.
302     else{
303       AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
304       return GetModuleNumber(cluster);
305     }
306         }//PHOS
307         
308         return -1;
309 }
310
311 //____________________________________________________________________________________________________________________________________________________
312 Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
313 {
314         //Get the EMCAL/PHOS module number that corresponds to this cluster
315         TLorentzVector lv;
316         Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
317         cluster->GetMomentum(lv,v);
318         Float_t phi = lv.Phi();
319         if(phi < 0) phi+=TMath::TwoPi();        
320         Int_t absId = -1;
321         if(cluster->IsEMCAL()){
322                 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
323                 if(fDebug > 2) 
324                   printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
325                          lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
326                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
327         }//EMCAL
328         else if(cluster->IsPHOS()) {
329                 Int_t    relId[4];
330                 if ( cluster->GetNCells() > 0) {
331                         absId = cluster->GetCellAbsId(0);
332                         if(fDebug > 2) 
333                                 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
334                                            lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
335                 }
336                 else return -1;
337                 
338                 if ( absId >= 0) {
339                         fPHOSGeo->AbsToRelNumbering(absId,relId);
340                         if(fDebug > 2) 
341                           printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
342                         return relId[0]-1;
343                 }
344                 else return -1;
345         }//PHOS
346         
347         return -1;
348 }
349
350 //_____________________________________________________________________________________________________________
351 Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(const Int_t absId, const TString calo, Int_t & icol, Int_t & irow, Int_t & iRCU) const
352 {
353         //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
354         Int_t imod = -1;
355         if ( absId >= 0) {
356                 if(calo=="EMCAL"){
357                         Int_t iTower = -1, iIphi = -1, iIeta = -1; 
358                         fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
359                         fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
360       if(imod < 0 || irow < 0 || icol < 0 ) {
361         Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
362       }
363       
364                         //RCU0
365                         if (0<=irow&&irow<8) iRCU=0; // first cable row
366                         else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half; 
367                         //second cable row
368                         //RCU1
369                         else if(8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half; 
370                         //second cable row
371                         else if(16<=irow&&irow<24) iRCU=1; // third cable row
372                         
373                         if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
374                         if (iRCU<0) {
375                                 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
376                         }                       
377                         
378                         return imod ;
379                 }//EMCAL
380                 else{//PHOS
381                         Int_t    relId[4];
382                         fPHOSGeo->AbsToRelNumbering(absId,relId);
383                         irow = relId[2];
384                         icol = relId[3];
385                         imod = relId[0]-1;
386                         iRCU= (Int_t)(relId[2]-1)/16 ;
387                         //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
388                         if (iRCU >= 4) {
389                                 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
390                         }                       
391                         return imod;
392                 }//PHOS 
393         }
394         
395         return -1;
396 }
397
398 //_______________________________________________________________
399 //void AliCalorimeterUtils::Init()
400 //{
401 //      //Init reader. Method to be called in AliAnaPartCorrMaker
402 //      
403 //      fEMCALBadChannelMap->SetName(Form("EMCALBadMap_%s",fTaskName.Data()));
404 //      fPHOSBadChannelMap->SetName(Form("PHOSBadMap_%s",fTaskName.Data()));
405 //}
406
407 //_______________________________________________________________
408 void AliCalorimeterUtils::InitParameters()
409 {
410   //Initialize the parameters of the analysis.
411   fEMCALGeoName = "EMCAL_FIRSTYEARV1";
412   fPHOSGeoName  = "PHOSgeo";
413         
414   if(gGeoManager) {// geoManager was set
415         if(fDebug > 2)printf("AliCalorimeterUtils::InitParameters() - Geometry manager available\n");
416         fEMCALGeoMatrixSet = kTRUE;      
417         fPHOSGeoMatrixSet  = kTRUE;      
418   }
419   else{
420         fEMCALGeoMatrixSet = kFALSE;
421         fPHOSGeoMatrixSet  = kFALSE;
422   }
423                 
424   fRemoveBadChannels   = kFALSE;
425         
426   fNCellsFromPHOSBorder    = 0;
427 }
428
429
430 //________________________________________________________________
431 void AliCalorimeterUtils::InitPHOSBadChannelStatusMap(){
432   //Init PHOS bad channels map
433   if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
434   //In order to avoid rewriting the same histograms
435   Bool_t oldStatus = TH1::AddDirectoryStatus();
436   TH1::AddDirectory(kFALSE);
437
438   fPHOSBadChannelMap = new TObjArray(5);        
439   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));
440
441   fPHOSBadChannelMap->SetOwner(kTRUE);
442   fPHOSBadChannelMap->Compress();
443   
444   //In order to avoid rewriting the same histograms
445   TH1::AddDirectory(oldStatus);         
446 }
447
448 //________________________________________________________________
449 void AliCalorimeterUtils::InitPHOSRecalibrationFactors(){
450         //Init EMCAL recalibration factors
451         if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
452         //In order to avoid rewriting the same histograms
453         Bool_t oldStatus = TH1::AddDirectoryStatus();
454         TH1::AddDirectory(kFALSE);
455
456         fPHOSRecalibrationFactors = new TObjArray(5);
457         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));
458         //Init the histograms with 1
459         for (Int_t m = 0; m < 5; m++) {
460                 for (Int_t i = 0; i < 56; i++) {
461                         for (Int_t j = 0; j < 64; j++) {
462                                 SetPHOSChannelRecalibrationFactor(m,i,j,1.);
463                         }
464                 }
465         }
466         fPHOSRecalibrationFactors->SetOwner(kTRUE);
467         fPHOSRecalibrationFactors->Compress();
468         
469         //In order to avoid rewriting the same histograms
470         TH1::AddDirectory(oldStatus);           
471 }
472
473
474 //________________________________________________________________
475 void AliCalorimeterUtils::InitEMCALGeometry()
476 {
477         //Initialize EMCAL geometry if it did not exist previously
478         if (!fEMCALGeo){
479                 fEMCALGeo = new AliEMCALGeoUtils(fEMCALGeoName);
480                 if(fDebug > 0){
481                         printf("AliCalorimeterUtils::InitEMCALGeometry()");
482                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
483                         printf("\n");
484                 }
485         }
486 }
487
488 //________________________________________________________________
489 void AliCalorimeterUtils::InitPHOSGeometry()
490 {
491         //Initialize PHOS geometry if it did not exist previously
492         if (!fPHOSGeo){
493                 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName); 
494                 if(fDebug > 0){
495                         printf("AliCalorimeterUtils::InitPHOSGeometry()");
496                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
497                         printf("\n");
498                 }       
499         }       
500 }
501
502 //________________________________________________________________
503 void AliCalorimeterUtils::Print(const Option_t * opt) const
504 {
505
506   //Print some relevant parameters set for the analysis
507   if(! opt)
508     return;
509
510   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
511   printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
512   printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
513          fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
514   if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
515   printf("Recalibrate Clusters? %d\n",fRecalibration);
516   printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
517   printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
518   printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
519
520   printf("    \n") ;
521
522
523 //________________________________________________________________
524 Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster, AliVCaloCells * cells){
525         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
526
527   //Initialize some used variables
528         Float_t energy   = 0;
529         Int_t absId      = -1;
530         Int_t icol = -1, irow = -1, iRCU = -1, module=1;
531         Float_t factor = 1, frac = 0;  
532   
533         if(cells) {
534         
535         //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
536         UShort_t * index    = cluster->GetCellsAbsId() ;
537         Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
538         Int_t ncells     = cluster->GetNCells();        
539         TString calo     = "EMCAL";
540         if(cluster->IsPHOS()) calo = "PHOS";
541         
542         //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
543         for(Int_t icell = 0; icell < ncells; icell++){
544                 absId = index[icell];
545                 frac =  fraction[icell];
546                 if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
547         module = GetModuleNumberCellIndexes(absId,calo,icol,irow,iRCU);
548                 if(cluster->IsPHOS()) factor = GetPHOSChannelRecalibrationFactor (module,icol,irow);
549                 else                  factor = GetEMCALChannelRecalibrationFactor(module,icol,irow);
550                 if(fDebug>2)
551                         printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, module %d, col %d, row %d, cell fraction %f, recalibration factor %f, cell energy %f\n", 
552                                 calo.Data(),module,icol,irow,frac,factor,cells->GetCellAmplitude(absId));
553                 
554                 energy += cells->GetCellAmplitude(absId)*factor*frac;
555         }
556         
557         if(fDebug>1)
558                 printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
559     
560         }// cells available
561   else{
562     Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
563   }
564   
565         return energy;
566 }
567
568 //________________________________________________________________
569 void AliCalorimeterUtils::SetGeometryTransformationMatrices(AliVEvent* inputEvent) 
570 {
571   //Set the calorimeters transformation matrices
572
573   //Get the EMCAL transformation geometry matrices from ESD 
574   if(!fEMCALGeoMatrixSet && fEMCALGeo){
575     if(fLoadEMCALMatrices){
576       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined geometry matrices\n");
577       for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
578         if(fEMCALMatrix[mod]){
579           if(fDebug > 1) 
580             fEMCALMatrix[mod]->Print();
581           fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;  
582         }
583       }//SM loop
584       fEMCALGeoMatrixSet = kTRUE;//At least one, so good
585       
586     }//Load matrices
587     else if (!gGeoManager) { 
588
589       if(fDebug > 1) 
590         printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load EMCAL misalignment matrices. \n");
591       if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  {
592         for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ 
593           //printf("Load matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
594           if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod)) {
595             fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
596           }
597         }// loop over super modules     
598         fEMCALGeoMatrixSet = kTRUE;//At least one, so good
599         
600       }//ESD as input
601       else {
602         if(fDebug > 1)
603           printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
604       }//AOD as input
605     }//Get matrix from data
606   }//EMCAL geo && no geoManager
607         
608         //Get the PHOS transformation geometry matrices from ESD 
609   if(!fPHOSGeoMatrixSet && fPHOSGeo){
610     if(fLoadPHOSMatrices){
611       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined geometry matrices\n");
612       for(Int_t mod=0; mod < 5; mod++){
613         if(fPHOSMatrix[mod]){
614           if(fDebug > 1) 
615             fPHOSMatrix[mod]->Print();
616           fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;  
617         }
618       }//SM loop
619       fPHOSGeoMatrixSet = kTRUE;//At least one, so good
620     }//Load matrices
621     else if (!gGeoManager) { 
622       if(fDebug > 1) 
623         printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
624                         if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  {
625                                 for(Int_t mod=0; mod < 5; mod++){ 
626                                         if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod)) {
627                                                 //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
628                                                 fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
629                                         }
630                                 }// loop over modules
631         fPHOSGeoMatrixSet  = kTRUE; //At least one so good
632                         }//ESD as input
633                         else {
634                                 if(fDebug > 1) 
635                                         printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
636       }//AOD as input
637     }// get matrix from data
638         }//PHOS geo     and  geoManager was not set
639   
640 }
641
642 //________________________________________________________________
643 void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu){
644   
645   //Recalculate EMCAL cluster position
646   
647   fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
648
649 }