coverity
[u/mrichter/AliRoot.git] / PWG / CaloTrackCorrBase / 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
16 //_________________________________________________________________________
17 // Class utility for Calorimeter specific selection methods                ///
18 //
19 //
20 //
21 //-- Author: Gustavo Conesa (LPSC-Grenoble) 
22 //////////////////////////////////////////////////////////////////////////////
23
24
25 // --- ROOT system ---
26 #include <TGeoManager.h>
27 #include <TH2F.h>
28 #include <TCanvas.h>
29 #include <TStyle.h>
30 #include <TPad.h>
31
32 // --- ANALYSIS system ---
33 #include "AliCalorimeterUtils.h"
34 #include "AliESDEvent.h"
35 #include "AliMCEvent.h"
36 #include "AliStack.h"
37 #include "AliAODPWG4Particle.h"
38 #include "AliVCluster.h"
39 #include "AliVCaloCells.h"
40 #include "AliMixedEvent.h"
41 #include "AliAODCaloCluster.h"
42
43 // --- Detector ---
44 #include "AliEMCALGeometry.h"
45 #include "AliPHOSGeoUtils.h"
46
47 ClassImp(AliCalorimeterUtils)
48   
49   
50 //____________________________________________
51   AliCalorimeterUtils::AliCalorimeterUtils() : 
52     TObject(), fDebug(0), 
53     fEMCALGeoName("EMCAL_COMPLETE12SMV1"),
54     fPHOSGeoName ("PHOSgeo"), 
55     fEMCALGeo(0x0),                   fPHOSGeo(0x0), 
56     fEMCALGeoMatrixSet(kFALSE),       fPHOSGeoMatrixSet(kFALSE), 
57     fLoadEMCALMatrices(kFALSE),       fLoadPHOSMatrices(kFALSE),
58     fRemoveBadChannels(kFALSE),       fPHOSBadChannelMap(0x0), 
59     fNCellsFromPHOSBorder(0),
60     fNMaskCellColumns(0),             fMaskCellColumns(0x0),
61     fRecalibration(kFALSE),           fPHOSRecalibrationFactors(),
62     fEMCALRecoUtils(new AliEMCALRecoUtils),
63     fRecalculatePosition(kFALSE),     fCorrectELinearity(kFALSE),
64     fRecalculateMatching(kFALSE),
65     fCutR(20),                        fCutZ(20),
66     fCutEta(20),                      fCutPhi(20),
67     fLocMaxCutE(0),                   fLocMaxCutEDiff(0),
68     fPlotCluster(0)
69 {
70   //Ctor
71   
72   //Initialize parameters
73   InitParameters();
74   for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
75   for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix [i] = 0 ;
76   
77 }
78
79 //_________________________________________
80 AliCalorimeterUtils::~AliCalorimeterUtils()
81 {
82   //Dtor
83   
84   //if(fPHOSGeo)  delete fPHOSGeo  ;
85   if(fEMCALGeo) delete fEMCALGeo ;
86   
87   if(fPHOSBadChannelMap) { 
88     fPHOSBadChannelMap->Clear();
89     delete  fPHOSBadChannelMap;
90   }
91         
92   if(fPHOSRecalibrationFactors) { 
93     fPHOSRecalibrationFactors->Clear();
94     delete  fPHOSRecalibrationFactors;
95   }
96         
97   if(fEMCALRecoUtils)   delete fEMCALRecoUtils ;
98   if(fNMaskCellColumns) delete [] fMaskCellColumns;
99   
100 }
101
102 //______________________________________________________________________________________
103 Bool_t AliCalorimeterUtils::AreNeighbours(const TString calo, 
104                                           const Int_t absId1, const Int_t absId2 ) const
105 {
106   // Tells if (true) or not (false) two cells are neighbours
107   // A neighbour is defined as being two cells which share a side or corner
108         
109   Bool_t areNeighbours = kFALSE ;
110   
111   Int_t iRCU1 = -1, irow1 = -1, icol1 = -1;
112   Int_t iRCU2 = -1, irow2 = -1, icol2 = -1;
113   
114   Int_t rowdiff =  0, coldiff =  0;
115   
116   Int_t nSupMod1 = GetModuleNumberCellIndexes(absId1, calo, icol1, irow1, iRCU1); 
117   Int_t nSupMod2 = GetModuleNumberCellIndexes(absId2, calo, icol2, irow2, iRCU2); 
118   
119   if(calo=="EMCAL" && nSupMod1!=nSupMod2)
120   {
121     // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
122     // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
123     if(nSupMod1%2) icol1+=AliEMCALGeoParams::fgkEMCALCols;
124     else           icol2+=AliEMCALGeoParams::fgkEMCALCols;    
125         }
126   
127   rowdiff = TMath::Abs( irow1 - irow2 ) ;  
128   coldiff = TMath::Abs( icol1 - icol2 ) ;  
129   
130   if (( coldiff <= 1 )  && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0)) 
131     areNeighbours = kTRUE ;
132   
133   return areNeighbours;
134   
135 }
136
137
138 //_____________________________________________________________________________________
139 Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster, 
140                                                     AliVCaloCells* cells, 
141                                                     AliVEvent * event, Int_t iev) const 
142 {
143   
144         // Given the list of AbsId of the cluster, get the maximum cell and 
145         // check if there are fNCellsFromBorder from the calorimeter border
146         
147   //If the distance to the border is 0 or negative just exit accept all clusters
148         if(cells->GetType()==AliVCaloCells::kEMCALCell && fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder() <= 0 ) return kTRUE;
149         if(cells->GetType()==AliVCaloCells::kPHOSCell  && fNCellsFromPHOSBorder  <= 0 ) return kTRUE;
150   
151   Int_t absIdMax        = -1;
152         Float_t ampMax  = -1;
153         
154   AliMixedEvent * mixEvent = dynamic_cast<AliMixedEvent*> (event);
155   Int_t nMixedEvents = 0 ; 
156   Int_t * cellsCumul = NULL ;
157   Int_t numberOfCells = 0 ;  
158   if (mixEvent){
159     nMixedEvents = mixEvent->GetNumberOfEvents() ; 
160     if (cells->GetType()==AliVCaloCells::kEMCALCell) {
161       cellsCumul =  mixEvent->GetEMCALCellsCumul() ; 
162       numberOfCells = mixEvent->GetNumberOfEMCALCells() ;
163     } 
164     
165     else if (cells->GetType()==AliVCaloCells::kPHOSCell) {
166       cellsCumul =  mixEvent->GetPHOSCellsCumul() ; 
167       numberOfCells = mixEvent->GetNumberOfPHOSCells() ;
168     } 
169     
170     if(cellsCumul){
171       
172       Int_t startCell = cellsCumul[iev] ; 
173       Int_t endCell   = (iev+1 < nMixedEvents)?cellsCumul[iev+1]:numberOfCells;
174       //Find cells with maximum amplitude
175       for(Int_t i = 0; i < cluster->GetNCells() ; i++){
176         Int_t absId = cluster->GetCellAbsId(i) ;
177         for (Int_t j = startCell; j < endCell ;  j++) {
178           Short_t cellNumber; 
179           Double_t amp ; 
180           Double_t time; 
181           cells->GetCell(j, cellNumber, amp, time) ; 
182           if (absId == cellNumber) {
183             if(amp > ampMax){
184               ampMax   = amp;
185               absIdMax = absId;
186             }        
187           }
188         }
189       }//loop on cluster cells
190     }// cells cumul available
191     else {
192       printf("AliCalorimeterUtils::CheckCellFiducialRegion() - CellsCumul is NULL!!!\n");
193       abort();
194     }
195   } else {//Normal SE Events
196     for(Int_t i = 0; i < cluster->GetNCells() ; i++){
197       Int_t absId = cluster->GetCellAbsId(i) ;
198       Float_t amp       = cells->GetCellAmplitude(absId);
199       if(amp > ampMax){
200         ampMax   = amp;
201         absIdMax = absId;
202       }
203     }
204   }
205         
206         if(fDebug > 1)
207                 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n", 
208            absIdMax, ampMax, cluster->E());
209         
210         if(absIdMax==-1) return kFALSE;
211         
212         //Check if the cell is close to the borders:
213         Bool_t okrow = kFALSE;
214         Bool_t okcol = kFALSE;
215   
216         if(cells->GetType()==AliVCaloCells::kEMCALCell){
217     
218                 Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1; 
219                 fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta); 
220                 fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
221                 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
222       Fatal("CheckCellFidutialRegion","Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",iSM,ieta,iphi);
223     }
224     
225                 //Check rows/phi
226     Int_t nborder = fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder();
227                 if(iSM < 10){
228                         if(iphi >= nborder && iphi < 24-nborder) okrow =kTRUE; 
229     }
230                 else{
231                         if(iphi >= nborder && iphi < 12-nborder) okrow =kTRUE; 
232                 }
233                 
234                 //Check columns/eta
235                 if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0()){
236                         if(ieta  > nborder && ieta < 48-nborder) okcol =kTRUE; 
237                 }
238                 else{
239                         if(iSM%2==0){
240                                 if(ieta >= nborder)     okcol = kTRUE;  
241                         }
242                         else {
243                                 if(ieta <  48-nborder)  okcol = kTRUE;  
244                         }
245                 }//eta 0 not checked
246                 if(fDebug > 1)
247                 {
248                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
249              nborder, ieta, iphi, iSM);
250                         if (okcol && okrow ) printf(" YES \n");
251                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
252                 }
253         }//EMCAL
254         else if(cells->GetType()==AliVCaloCells::kPHOSCell){
255                 Int_t relId[4];
256                 Int_t irow = -1, icol = -1;
257                 fPHOSGeo->AbsToRelNumbering(absIdMax,relId);
258                 irow = relId[2];
259                 icol = relId[3];
260                 //imod = relId[0]-1;
261                 if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE; 
262                 if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE; 
263                 if(fDebug > 1)
264                 {
265                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - PHOS Cluster in %d cells fiducial volume: icol %d, irow %d, Module %d?",
266              fNCellsFromPHOSBorder, icol, irow, relId[0]-1);
267                         if (okcol && okrow ) printf(" YES \n");
268                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
269                 }
270         }//PHOS
271         
272         if (okcol && okrow) return kTRUE; 
273         else                return kFALSE;
274         
275 }       
276
277 //_________________________________________________________________________________________________________
278 Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter,UShort_t* cellList, Int_t nCells)
279 {
280         // Check that in the cluster cells, there is no bad channel of those stored 
281         // in fEMCALBadChannelMap or fPHOSBadChannelMap
282         
283         if (!fRemoveBadChannels) return kFALSE;
284         //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
285         if(calorimeter == "EMCAL" && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
286         if(calorimeter == "PHOS"  && !fPHOSBadChannelMap)  return kFALSE;
287   
288         Int_t icol = -1;
289         Int_t irow = -1;
290         Int_t imod = -1;
291         for(Int_t iCell = 0; iCell<nCells; iCell++){
292     
293                 //Get the column and row
294                 if(calorimeter == "EMCAL"){
295       return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
296                 }
297                 else if(calorimeter=="PHOS"){
298                         Int_t    relId[4];
299                         fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
300                         irow = relId[2];
301                         icol = relId[3];
302                         imod = relId[0]-1;
303                         if(fPHOSBadChannelMap->GetEntries() <= imod)continue;
304       //printf("PHOS bad channels imod %d, icol %d, irow %d\n",imod, irow, icol);
305                         if(GetPHOSChannelStatus(imod, irow, icol)) return kTRUE;
306                 }
307                 else return kFALSE;
308                 
309         }// cell cluster loop
310   
311         return kFALSE;
312   
313 }
314
315 //_______________________________________________________________
316 void AliCalorimeterUtils::CorrectClusterEnergy(AliVCluster *clus)
317 {
318   // Correct cluster energy non linearity
319   
320   clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
321
322 }
323
324 //________________________________________________________________________________________
325 Int_t  AliCalorimeterUtils::GetMaxEnergyCell(AliVCaloCells* cells, const AliVCluster* clu, 
326                                              Float_t & clusterFraction) const 
327 {
328   
329   //For a given CaloCluster gets the absId of the cell 
330   //with maximum energy deposit.
331   
332   if( !clu || !cells ){
333     AliInfo("Cluster or cells pointer is null!");
334     return -1;
335   }
336   
337   Double_t eMax        =-1.;
338   Double_t eTot        = 0.;
339   Double_t eCell       =-1.;
340   Float_t  fraction    = 1.;
341   Float_t  recalFactor = 1.;
342   Int_t    cellAbsId   =-1 , absId =-1 ;
343   Int_t    iSupMod     =-1 , ieta  =-1 , iphi = -1, iRCU = -1;
344   
345   TString           calo = "EMCAL";
346   if(clu->IsPHOS()) calo = "PHOS";
347   
348   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
349     
350     cellAbsId = clu->GetCellAbsId(iDig);
351     
352     fraction  = clu->GetCellAmplitudeFraction(iDig);
353     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
354     
355     iSupMod = GetModuleNumberCellIndexes(cellAbsId, calo, ieta, iphi, iRCU);
356     
357     if(IsRecalibrationOn()) {
358       if(calo=="EMCAL") recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
359       else              recalFactor = GetPHOSChannelRecalibrationFactor (iSupMod,iphi,ieta);
360     }
361     
362     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
363     
364     if(eCell > eMax)  { 
365       eMax  = eCell; 
366       absId = cellAbsId;
367     }
368     
369     eTot+=eCell;
370     
371   }// cell loop
372   
373   clusterFraction = (eTot-eMax)/eTot; //Do not use cluster energy in case it was corrected for non linearity.
374   
375   return absId;
376   
377 }
378
379 //__________________________________________________________________________
380 AliVTrack * AliCalorimeterUtils::GetMatchedTrack(const AliVCluster* cluster, 
381                                                  const AliVEvent* event, 
382                                                  const Int_t index) const
383 {
384   // Get the matched track given its index, usually just the first match
385   // Since it is different for ESDs and AODs here it is a wrap method to do it
386   
387   AliVTrack *track = 0;
388   
389   // EMCAL case only when matching is recalculated
390   if(cluster->IsEMCAL() && IsRecalculationOfClusterTrackMatchingOn())
391   {
392     Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(cluster->GetID());
393     //printf("track index %d, cluster ID %d \n ",trackIndex,cluster->GetID());
394     
395     if(trackIndex < 0 )
396     { 
397       printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d, from recalculation\n", trackIndex);
398     }
399     else 
400     {
401       track = dynamic_cast<AliVTrack*> (event->GetTrack(trackIndex));
402     }
403
404     return track ;
405     
406   }   
407   
408   // Normal case, get info from ESD or AOD
409   // ESDs
410   if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
411   {
412     Int_t iESDtrack = cluster->GetTrackMatchedIndex();
413     
414     if(iESDtrack < 0 )
415     { 
416       printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d\n", index);
417       return 0x0;
418     }
419     
420     track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
421     
422   }
423   else // AODs
424   {
425     if(cluster->GetNTracksMatched() > 0 )
426       track = dynamic_cast<AliVTrack*>(cluster->GetTrackMatched(index));
427   }
428   
429   return track ;
430   
431 }
432
433 //_____________________________________________________________________________________________________
434 Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
435 {
436         //Get the EMCAL/PHOS module number that corresponds to this particle
437         
438         Int_t absId = -1;
439         if(particle->GetDetector()=="EMCAL"){
440                 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
441                 if(fDebug > 2) 
442                   printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
443              particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
444                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
445         }//EMCAL
446         else if(particle->GetDetector()=="PHOS"){
447     // In case we use the MC reader, the input are TParticles, 
448     // in this case use the corresponing method in PHOS Geometry to get the particle.
449     if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
450     {
451       Int_t mod =-1;
452       Double_t z = 0., x=0.;
453       TParticle* primary = 0x0;
454       AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
455       if(stack) {
456         primary = stack->Particle(particle->GetCaloLabel(0));
457       }
458       else {
459         Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
460       }
461       
462       if(primary){
463         fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
464       }
465       else{
466         Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
467       }
468       return mod;
469     }
470     // Input are ESDs or AODs, get the PHOS module number like this.
471     else{
472       //FIXME
473       //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
474       //return GetModuleNumber(cluster);
475       //MEFIX
476       return -1;
477     }
478         }//PHOS
479         
480         return -1;
481 }
482
483 //_____________________________________________________________________
484 Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
485 {
486         //Get the EMCAL/PHOS module number that corresponds to this cluster
487         TLorentzVector lv;
488         Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
489   if(!cluster){
490     if(fDebug > 1) printf("AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
491     return -1;
492   }
493         cluster->GetMomentum(lv,v);
494         Float_t phi = lv.Phi();
495         if(phi < 0) phi+=TMath::TwoPi();        
496         Int_t absId = -1;
497         if(cluster->IsEMCAL()){
498                 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
499                 if(fDebug > 2) 
500                   printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
501              lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
502                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
503         }//EMCAL
504         else if(cluster->IsPHOS()) {
505                 Int_t    relId[4];
506                 if ( cluster->GetNCells() > 0) {
507                         absId = cluster->GetCellAbsId(0);
508                         if(fDebug > 2) 
509                                 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
510                lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
511                 }
512                 else return -1;
513                 
514                 if ( absId >= 0) {
515                         fPHOSGeo->AbsToRelNumbering(absId,relId);
516                         if(fDebug > 2) 
517                           printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
518                         return relId[0]-1;
519                 }
520                 else return -1;
521         }//PHOS
522         
523         return -1;
524 }
525
526 //___________________________________________________________________________________________________
527 Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(const Int_t absId, const TString calo, 
528                                                       Int_t & icol, Int_t & irow, Int_t & iRCU) const
529 {
530         //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
531         Int_t imod = -1;
532         if ( absId >= 0) {
533                 if(calo=="EMCAL"){
534                         Int_t iTower = -1, iIphi = -1, iIeta = -1; 
535                         fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
536                         fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
537       if(imod < 0 || irow < 0 || icol < 0 ) {
538         Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
539       }
540       
541                         //RCU0
542                         if (0<=irow&&irow<8) iRCU=0; // first cable row
543                         else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half; 
544                         //second cable row
545                         //RCU1
546                         else if(8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half; 
547                         //second cable row
548                         else if(16<=irow&&irow<24) iRCU=1; // third cable row
549                         
550                         if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
551                         if (iRCU<0) {
552                                 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
553                         }                       
554                         
555                         return imod ;
556                 }//EMCAL
557                 else{//PHOS
558                         Int_t    relId[4];
559                         fPHOSGeo->AbsToRelNumbering(absId,relId);
560                         irow = relId[2];
561                         icol = relId[3];
562                         imod = relId[0]-1;
563                         iRCU= (Int_t)(relId[2]-1)/16 ;
564                         //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
565                         if (iRCU >= 4) {
566                                 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
567                         }                       
568                         return imod;
569                 }//PHOS 
570         }
571         
572         return -1;
573 }
574
575 //___________________________________________________________________________________________
576 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells) 
577 {
578   // Find local maxima in cluster
579   
580   const Int_t   nc = cluster->GetNCells();
581   Int_t   absIdList[nc]; 
582   Float_t maxEList[nc]; 
583   
584   Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList);
585   
586   return nMax;
587   
588 }
589
590 //___________________________________________________________________________________________
591 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
592                                                   Int_t *absIdList,     Float_t *maxEList) 
593 {
594   // Find local maxima in cluster
595   
596   Int_t iDigitN = 0 ;
597   Int_t iDigit  = 0 ;
598   Int_t absId1 = -1 ;
599   Int_t absId2 = -1 ;
600   const Int_t nCells = cluster->GetNCells();
601   
602   TString calorimeter = "EMCAL";
603   if(!cluster->IsEMCAL()) calorimeter = "PHOS";
604   
605   //printf("cluster : ncells %d \n",nCells);
606   
607   Float_t emax  = 0;
608   Int_t   idmax =-1;
609   for(iDigit = 0; iDigit < nCells ; iDigit++)
610   {
611     absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit]  ; 
612     Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
613     RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);  
614     if( en > emax )
615     {
616       emax  = en ;
617       idmax = absIdList[iDigit] ;
618     }
619     //Int_t icol = -1, irow = -1, iRCU = -1;
620     //Int_t sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
621     //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
622   }
623   
624   for(iDigit = 0 ; iDigit < nCells; iDigit++) 
625   {   
626     if(absIdList[iDigit]>=0) 
627     {
628       absId1 = cluster->GetCellsAbsId()[iDigit];
629       
630       Float_t en1 = cells->GetCellAmplitude(absId1);
631       RecalibrateCellAmplitude(en1,calorimeter,absId1);  
632       
633       //printf("%d : absIDi %d, E %f\n",iDigit, absId1,en1);
634       
635       for(iDigitN = 0; iDigitN < nCells; iDigitN++) 
636       { 
637         absId2 = cluster->GetCellsAbsId()[iDigitN] ;
638         
639         if(absId2==-1 || absId2==absId1) continue;
640         
641         //printf("\t %d : absIDj %d\n",iDigitN, absId2);
642         
643         Float_t en2 = cells->GetCellAmplitude(absId2);
644         RecalibrateCellAmplitude(en2,calorimeter,absId2);
645         
646         //printf("\t %d : absIDj %d, E %f\n",iDigitN, absId2,en2);
647         
648         if ( AreNeighbours(calorimeter, absId1, absId2) ) 
649         {
650           // printf("\t \t Neighbours \n");
651           if (en1 > en2 ) 
652           {    
653             absIdList[iDigitN] = -1 ;
654             //printf("\t \t indexN %d not local max\n",iDigitN);
655             // but may be digit too is not local max ?
656             if(en1 < en2 + fLocMaxCutEDiff) {
657               //printf("\t \t index %d not local max cause locMaxCutEDiff\n",iDigit);
658               absIdList[iDigit] = -1 ;
659             }
660           }
661           else 
662           {
663             absIdList[iDigit] = -1 ;
664             //printf("\t \t index %d not local max\n",iDigitN);
665             // but may be digitN too is not local max ?
666             if(en1 > en2 - fLocMaxCutEDiff) 
667             {
668               absIdList[iDigitN] = -1 ; 
669               //printf("\t \t indexN %d not local max cause locMaxCutEDiff\n",iDigit);
670             }
671           } 
672         } // if Are neighbours
673         //else printf("\t \t NOT Neighbours \n");
674       } // while digitN
675     } // slot not empty
676   } // while digit
677   
678   iDigitN = 0 ;
679   for(iDigit = 0; iDigit < nCells; iDigit++) 
680   { 
681     if(absIdList[iDigit]>=0 )
682     {
683       absIdList[iDigitN] = absIdList[iDigit] ;
684       Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
685       RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);  
686       if(en < fLocMaxCutE) continue; // Maxima only with seed energy at least
687       maxEList[iDigitN] = en ;
688       //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
689       iDigitN++ ; 
690     }
691   }
692   
693   if(iDigitN == 0)
694   {
695     if(fDebug > 0) 
696       printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - No local maxima found, assign highest energy cell as maxima, id %d, en cell %2.2f, en cluster %2.2f\n",
697              idmax,emax,cluster->E());
698     iDigitN      = 1     ;
699     maxEList[0]  = emax  ;
700     absIdList[0] = idmax ; 
701   }
702   
703   if(fDebug > 0) 
704   {    
705     printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - In cluster E %2.2f, M02 %2.2f, M20 %2.2f, N maxima %d \n", 
706            cluster->E(),cluster->GetM02(),cluster->GetM20(), iDigitN);
707   
708     if(fDebug > 1) for(Int_t imax = 0; imax < iDigitN; imax++) 
709     {
710       printf(" \t i %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
711     }
712   }
713   
714   return iDigitN ;
715   
716 }
717
718 //________________________________________
719 void AliCalorimeterUtils::InitParameters()
720 {
721   //Initialize the parameters of the analysis.
722   
723   fEMCALGeoName         = "EMCAL_COMPLETEV1";
724   fPHOSGeoName          = "PHOSgeo";
725   
726   fEMCALGeoMatrixSet    = kFALSE;
727   fPHOSGeoMatrixSet     = kFALSE;
728   
729   fRemoveBadChannels    = kFALSE;
730   
731   fNCellsFromPHOSBorder = 0;
732   
733   fLocMaxCutE           = 0.1 ;
734   fLocMaxCutEDiff       = 0.0 ;
735
736   //  fMaskCellColumns = new Int_t[fNMaskCellColumns];
737   //  fMaskCellColumns[0] = 6 ;  fMaskCellColumns[1] = 7 ;  fMaskCellColumns[2] = 8 ; 
738   //  fMaskCellColumns[3] = 35;  fMaskCellColumns[4] = 36;  fMaskCellColumns[5] = 37; 
739   //  fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
740   //  fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols; 
741   //  fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols; 
742   
743 }
744
745
746 //_____________________________________________________
747 void AliCalorimeterUtils::InitPHOSBadChannelStatusMap()
748 {
749   //Init PHOS bad channels map
750   if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
751   //In order to avoid rewriting the same histograms
752   Bool_t oldStatus = TH1::AddDirectoryStatus();
753   TH1::AddDirectory(kFALSE);
754   
755   fPHOSBadChannelMap = new TObjArray(5);        
756   for (int i = 0; i < 5; i++)fPHOSBadChannelMap->Add(new TH2I(Form("PHOS_BadMap_mod%d",i),Form("PHOS_BadMap_mod%d",i), 64, 0, 64, 56, 0, 56));
757   
758   fPHOSBadChannelMap->SetOwner(kTRUE);
759   fPHOSBadChannelMap->Compress();
760   
761   //In order to avoid rewriting the same histograms
762   TH1::AddDirectory(oldStatus);         
763 }
764
765 //______________________________________________________
766 void AliCalorimeterUtils::InitPHOSRecalibrationFactors()
767 {
768         //Init EMCAL recalibration factors
769         if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
770         //In order to avoid rewriting the same histograms
771         Bool_t oldStatus = TH1::AddDirectoryStatus();
772         TH1::AddDirectory(kFALSE);
773   
774         fPHOSRecalibrationFactors = new TObjArray(5);
775         for (int i = 0; i < 5; i++)fPHOSRecalibrationFactors->Add(new TH2F(Form("PHOSRecalFactors_Mod%d",i),Form("PHOSRecalFactors_Mod%d",i), 64, 0, 64, 56, 0, 56));
776         //Init the histograms with 1
777         for (Int_t m = 0; m < 5; m++) {
778                 for (Int_t i = 0; i < 56; i++) {
779                         for (Int_t j = 0; j < 64; j++) {
780                                 SetPHOSChannelRecalibrationFactor(m,j,i,1.);
781                         }
782                 }
783         }
784         fPHOSRecalibrationFactors->SetOwner(kTRUE);
785         fPHOSRecalibrationFactors->Compress();
786         
787         //In order to avoid rewriting the same histograms
788         TH1::AddDirectory(oldStatus);           
789 }
790
791
792 //___________________________________________
793 void AliCalorimeterUtils::InitEMCALGeometry()
794 {
795         //Initialize EMCAL geometry if it did not exist previously
796         if (!fEMCALGeo){
797                 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
798                 if(fDebug > 0){
799                         printf("AliCalorimeterUtils::InitEMCALGeometry()");
800                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
801                         printf("\n");
802                 }
803         }
804 }
805
806 //__________________________________________
807 void AliCalorimeterUtils::InitPHOSGeometry()
808 {
809         //Initialize PHOS geometry if it did not exist previously
810         if (!fPHOSGeo){
811                 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName); 
812                 if(fDebug > 0){
813                         printf("AliCalorimeterUtils::InitPHOSGeometry()");
814                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
815                         printf("\n");
816                 }       
817         }       
818 }
819
820 //_________________________________________________________
821 void AliCalorimeterUtils::Print(const Option_t * opt) const
822 {
823   
824   //Print some relevant parameters set for the analysis
825   if(! opt)
826     return;
827   
828   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
829   printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
830   printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
831          fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
832   if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
833   printf("Recalibrate Clusters? %d\n",fRecalibration);
834   printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
835   printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
836   printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
837   
838   printf("Loc. Max. E > %2.2f\n",       fLocMaxCutE);
839   printf("Loc. Max. E Diff > %2.2f\n",  fLocMaxCutEDiff);
840   
841   printf("    \n") ;
842
843
844 //__________________________________________________________________
845 Bool_t AliCalorimeterUtils::MaskFrameCluster(const Int_t iSM,  
846                                              const Int_t ieta) const 
847 {
848   //Check if cell is in one of the regions where we have significant amount 
849   //of material in front. Only EMCAL
850   
851   Int_t icol = ieta;
852   if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
853   
854   if (fNMaskCellColumns && fMaskCellColumns) {
855     for (Int_t imask = 0; imask < fNMaskCellColumns; imask++) {
856       if(icol==fMaskCellColumns[imask]) return kTRUE;
857     }
858   }
859   
860   return kFALSE;
861   
862 }
863
864 //__________________________________________________________________________________________
865 void AliCalorimeterUtils::RecalibrateCellAmplitude(Float_t & amp,
866                                                    const TString calo, const Int_t id) const
867 {
868   //Recaculate cell energy if recalibration factor
869   
870   Int_t icol     = -1; Int_t irow     = -1; Int_t iRCU     = -1;
871   Int_t nModule  = GetModuleNumberCellIndexes(id,calo, icol, irow, iRCU);
872   
873   if (IsRecalibrationOn()) 
874   {
875     if(calo == "PHOS") 
876     {
877       amp *= GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
878     }
879     else                                   
880     {
881       amp *= GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
882     }
883   }
884 }
885
886 //_________________________________________________________________________________
887 void AliCalorimeterUtils::RecalibrateCellTime(Double_t & time, 
888                                               const TString calo, 
889                                               const Int_t id, const Int_t bc) const
890 {
891   // Recalculate time if time recalibration available for EMCAL
892   // not ready for PHOS
893   
894   if(calo == "EMCAL" && GetEMCALRecoUtils()->IsTimeRecalibrationOn()) 
895   {
896     GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
897   }
898   
899 }
900
901
902 //__________________________________________________________________________
903 Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster, 
904                                                       AliVCaloCells * cells)
905 {
906         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
907   
908   //Initialize some used variables
909         Float_t frac  = 0., energy = 0.;  
910   
911         if(cells) 
912   {
913     //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
914     
915     UShort_t * index    = cluster->GetCellsAbsId() ;
916     Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
917     
918     Int_t ncells     = cluster->GetNCells();    
919     
920     TString calo     = "EMCAL";
921     if(cluster->IsPHOS()) 
922       calo = "PHOS";
923     
924     //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
925     for(Int_t icell = 0; icell < ncells; icell++){
926       
927       Int_t absId = index[icell];
928       
929       frac =  fraction[icell];
930       if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
931       
932       Float_t amp = cells->GetCellAmplitude(absId);
933       RecalibrateCellAmplitude(amp,calo, absId);
934       
935       if(fDebug>2)
936         printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, cell fraction %f, cell energy %f\n", 
937                calo.Data(),frac,cells->GetCellAmplitude(absId));
938       
939       energy += amp*frac;
940     }
941     
942     if(fDebug>1)
943       printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
944     
945         }// cells available
946   else
947   {
948     Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
949   }
950   
951         return energy;
952 }
953
954 //________________________________________________________________________________
955 void AliCalorimeterUtils::SetGeometryTransformationMatrices(AliVEvent* inputEvent) 
956 {
957   //Set the calorimeters transformation matrices
958   
959   //Get the EMCAL transformation geometry matrices from ESD 
960   if(!fEMCALGeoMatrixSet && fEMCALGeo){
961     if(fLoadEMCALMatrices){
962       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined EMCAL geometry matrices\n");
963       for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
964         if(fEMCALMatrix[mod]){
965           if(fDebug > 1) 
966             fEMCALMatrix[mod]->Print();
967           fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;  
968         }
969       }//SM loop
970       fEMCALGeoMatrixSet = kTRUE;//At least one, so good
971       
972     }//Load matrices
973     else if (!gGeoManager) { 
974       
975       if(fDebug > 1) 
976         printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load EMCAL misalignment matrices. \n");
977       if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  {
978         for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ 
979           //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
980           if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod)) {
981             fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
982           }
983         }// loop over super modules     
984         
985         fEMCALGeoMatrixSet = kTRUE;//At least one, so good
986         
987       }//ESD as input
988       else {
989         if(fDebug > 1)
990           printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
991       }//AOD as input
992     }//Get matrix from data
993     else if(gGeoManager){
994       fEMCALGeoMatrixSet = kTRUE;
995     }
996   }//EMCAL geo && no geoManager
997   
998         //Get the PHOS transformation geometry matrices from ESD 
999   if(!fPHOSGeoMatrixSet && fPHOSGeo){
1000
1001     if(fLoadPHOSMatrices){
1002       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined PHOS geometry matrices\n");
1003       for(Int_t mod = 0 ; mod < 5 ; mod++){
1004         if(fPHOSMatrix[mod]){
1005           if(fDebug > 1 ) 
1006             fPHOSMatrix[mod]->Print();
1007           fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;  
1008         }
1009       }//SM loop
1010       fPHOSGeoMatrixSet = kTRUE;//At least one, so good
1011
1012     }//Load matrices
1013     else if (!gGeoManager) { 
1014       if(fDebug > 1) 
1015         printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
1016                         if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  {
1017                                 for(Int_t mod = 0; mod < 5; mod++){ 
1018                                         if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod)) {
1019                                                 //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
1020                                                 fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
1021                                         }
1022                                 }// loop over modules
1023         fPHOSGeoMatrixSet  = kTRUE; //At least one so good
1024                         }//ESD as input
1025                         else {
1026                                 if(fDebug > 1) 
1027                                         printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
1028       }//AOD as input
1029     }// get matrix from data
1030     else if(gGeoManager){
1031       fPHOSGeoMatrixSet = kTRUE;
1032     }
1033         }//PHOS geo     and  geoManager was not set
1034   
1035 }
1036
1037 //__________________________________________________________________________________________
1038 void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
1039 {
1040   
1041   //Recalculate EMCAL cluster position
1042   
1043   fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
1044   
1045 }
1046
1047 //___________________________________________________________________________
1048 void AliCalorimeterUtils::SplitEnergy(const Int_t absId1, const Int_t absId2,
1049                                       AliVCluster* cluster, 
1050                                       AliVCaloCells* cells,
1051                                       //Float_t & e1, Float_t & e2,
1052                                       AliAODCaloCluster* cluster1,
1053                                       AliAODCaloCluster* cluster2,
1054                                       const Int_t nMax,
1055                                       const Int_t eventNumber)
1056 {
1057   
1058   // Split energy of cluster between the 2 local maxima, sum energy on 3x3, and if the 2 
1059   // maxima are too close and have common cells, split the energy between the 2
1060   
1061   TH2F* hClusterMap    = 0 ;
1062   TH2F* hClusterLocMax = 0 ;
1063   TH2F* hCluster1      = 0 ;
1064   TH2F* hCluster2      = 0 ;
1065   
1066   if(fPlotCluster)
1067   {
1068     hClusterMap    = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
1069     hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
1070     hCluster1      = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
1071     hCluster2      = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
1072     hClusterMap    ->SetXTitle("column");
1073     hClusterMap    ->SetYTitle("row");
1074     hClusterLocMax ->SetXTitle("column");
1075     hClusterLocMax ->SetYTitle("row");
1076     hCluster1      ->SetXTitle("column");
1077     hCluster1      ->SetYTitle("row");
1078     hCluster2      ->SetXTitle("column");
1079     hCluster2      ->SetYTitle("row");
1080   }
1081   
1082   TString calorimeter = "EMCAL";
1083   if(cluster->IsPHOS())
1084   {
1085     calorimeter="PHOS";
1086     printf("AliCalorimeterUtils::SplitEnerg() Not supported for PHOS yet \n");
1087     return;
1088   }
1089   
1090   const Int_t ncells  = cluster->GetNCells();  
1091   Int_t absIdList[ncells]; 
1092   
1093   Float_t e1 = 0,  e2   = 0 ;
1094   Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;  
1095   Float_t eCluster = 0;
1096   Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1; 
1097   for(Int_t iDigit  = 0; iDigit < ncells; iDigit++ ) 
1098   {
1099     absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
1100     
1101     Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
1102     RecalibrateCellAmplitude(ec,calorimeter, absIdList[iDigit]);
1103     eCluster+=ec;
1104     
1105     if(fPlotCluster) 
1106     {
1107       //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
1108       sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1109       if(icol > maxCol) maxCol = icol;
1110       if(icol < minCol) minCol = icol;
1111       if(irow > maxRow) maxRow = irow;
1112       if(irow < minRow) minRow = irow;
1113       hClusterMap->Fill(icol,irow,ec);
1114     }
1115     
1116   }
1117   
1118   // Init counters and variables
1119   Int_t ncells1 = 1 ;
1120   UShort_t absIdList1[9] ;  
1121   Double_t fracList1 [9] ;  
1122   absIdList1[0] = absId1 ;
1123   fracList1 [0] = 1. ;
1124   
1125   Float_t ecell1 = cells->GetCellAmplitude(absId1);
1126   RecalibrateCellAmplitude(ecell1, calorimeter, absId1);
1127   e1 =  ecell1;  
1128   
1129   Int_t ncells2 = 1 ;
1130   UShort_t absIdList2[9] ;  
1131   Double_t fracList2 [9] ; 
1132   absIdList2[0] = absId2 ;
1133   fracList2 [0] = 1. ;
1134   
1135   Float_t ecell2 = cells->GetCellAmplitude(absId2);
1136   RecalibrateCellAmplitude(ecell2, calorimeter, absId2);
1137   e2 =  ecell2;  
1138   
1139   if(fPlotCluster)
1140   {
1141     Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
1142     sm = GetModuleNumberCellIndexes(absId1, calorimeter, icol1, irow1, iRCU) ;
1143     hClusterLocMax->Fill(icol1,irow1,ecell1);
1144     sm = GetModuleNumberCellIndexes(absId2, calorimeter, icol2, irow2, iRCU) ;
1145     hClusterLocMax->Fill(icol2,irow2,ecell2);
1146   }
1147   
1148   // Very rough way to share the cluster energy
1149   Float_t eRemain = (eCluster-ecell1-ecell2)/2;
1150   Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
1151   Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
1152   
1153   for(Int_t iDigit = 0; iDigit < ncells; iDigit++)
1154   {
1155     Int_t absId = absIdList[iDigit];
1156     
1157     if(absId==absId1 || absId==absId2 || absId < 0) continue;
1158     
1159     Float_t ecell = cells->GetCellAmplitude(absId);
1160     RecalibrateCellAmplitude(ecell, calorimeter, absId);
1161     
1162     if(AreNeighbours(calorimeter, absId1,absId ))
1163     { 
1164       absIdList1[ncells1]= absId;
1165       
1166       if(AreNeighbours(calorimeter, absId2,absId ))
1167       { 
1168         fracList1[ncells1] = shareFraction1; 
1169         e1 += ecell*shareFraction1;
1170       }
1171       else 
1172       {
1173         fracList1[ncells1] = 1.; 
1174         e1 += ecell;
1175       }
1176       
1177       ncells1++;
1178       
1179     } // neigbour to cell1
1180     
1181     if(AreNeighbours(calorimeter, absId2,absId ))
1182     { 
1183       absIdList2[ncells2]= absId;
1184       
1185       if(AreNeighbours(calorimeter, absId1,absId ))
1186       { 
1187         fracList2[ncells2] = shareFraction2; 
1188         e2 += ecell*shareFraction2;
1189       }
1190       else
1191       { 
1192         fracList2[ncells2] = 1.; 
1193         e2 += ecell;
1194       }
1195       
1196       ncells2++;
1197       
1198     } // neigbour to cell2
1199     
1200   }
1201   
1202   if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::SplitEnergy() - n Local Max %d, Cluster energy  = %f, Ecell1 = %f, Ecell2 = %f, Enew1 = %f, Enew2 = %f, Remain %f, \n ncells %d, ncells1 %d, ncells2 %d, f1 %f, f2  %f, sum f12 = %f \n",
1203                             nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2);
1204   
1205   cluster1->SetE(e1);
1206   cluster2->SetE(e2);  
1207   
1208   cluster1->SetNCells(ncells1);
1209   cluster2->SetNCells(ncells2);  
1210   
1211   cluster1->SetCellsAbsId(absIdList1);
1212   cluster2->SetCellsAbsId(absIdList2);
1213   
1214   cluster1->SetCellsAmplitudeFraction(fracList1);
1215   cluster2->SetCellsAmplitudeFraction(fracList2);
1216   
1217   if(calorimeter=="EMCAL")
1218   {
1219     GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
1220     GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
1221   }
1222   
1223   if(fPlotCluster)
1224   {
1225     //printf("Cells of cluster1: ");
1226     for(Int_t iDigit  = 0; iDigit < ncells1; iDigit++ ) 
1227     {
1228       //printf(" %d ",absIdList1[iDigit]);
1229       
1230       sm = GetModuleNumberCellIndexes(absIdList1[iDigit], calorimeter, icol, irow, iRCU) ;
1231       
1232       if( AreNeighbours(calorimeter, absId2,absIdList1[iDigit]) )
1233         hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit])*shareFraction1);
1234       else 
1235         hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit]));
1236     }
1237     
1238     //printf(" \n ");
1239     //printf("Cells of cluster2: ");
1240     
1241     for(Int_t iDigit  = 0; iDigit < ncells2; iDigit++ ) 
1242     {
1243       //printf(" %d ",absIdList2[iDigit]);
1244       
1245       sm = GetModuleNumberCellIndexes(absIdList2[iDigit], calorimeter, icol, irow, iRCU) ;
1246       if( AreNeighbours(calorimeter, absId1,absIdList2[iDigit]) )
1247         hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit])*shareFraction2);
1248       else
1249         hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit]));
1250       
1251     }
1252     //printf(" \n ");
1253     
1254     gStyle->SetPadRightMargin(0.1);
1255     gStyle->SetPadLeftMargin(0.1);
1256     gStyle->SetOptStat(0);
1257     gStyle->SetOptFit(000000);
1258     
1259     if(maxCol-minCol > maxRow-minRow)
1260     {
1261       maxRow+= maxCol-minCol;
1262     }
1263     else 
1264     {
1265       maxCol+= maxRow-minRow;
1266     }
1267     
1268     TCanvas  * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
1269     c->Divide(2,2);  
1270     c->cd(1);
1271     gPad->SetGridy();
1272     gPad->SetGridx();
1273     hClusterMap    ->SetAxisRange(minCol, maxCol,"X");
1274     hClusterMap    ->SetAxisRange(minRow, maxRow,"Y");
1275     hClusterMap    ->Draw("colz");
1276     c->cd(2);
1277     gPad->SetGridy();
1278     gPad->SetGridx();
1279     hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
1280     hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
1281     hClusterLocMax ->Draw("colz");
1282     c->cd(3);
1283     gPad->SetGridy();
1284     gPad->SetGridx();
1285     hCluster1      ->SetAxisRange(minCol, maxCol,"X");
1286     hCluster1      ->SetAxisRange(minRow, maxRow,"Y");
1287     hCluster1      ->Draw("colz");
1288     c->cd(4);
1289     gPad->SetGridy();
1290     gPad->SetGridx();
1291     hCluster2      ->SetAxisRange(minCol, maxCol,"X");
1292     hCluster2      ->SetAxisRange(minRow, maxRow,"Y");
1293     hCluster2      ->Draw("colz");
1294     
1295     if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",
1296                                    eventNumber,cluster->E(),nMax,ncells1,ncells2));
1297     
1298     delete c;
1299     delete hClusterMap;
1300     delete hClusterLocMax;
1301     delete hCluster1;
1302     delete hCluster2;
1303   }
1304 }
1305
1306 //________________________________________________________________________________
1307 void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event, 
1308                                                           TObjArray* clusterArray) 
1309
1310   //Recalculate track matching
1311   
1312   if (fRecalculateMatching) 
1313   {
1314     fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo)   ; 
1315     //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1316     //if(esdevent){
1317     //  fEMCALRecoUtils->SetClusterMatchedToTrack(esdevent)  ;
1318     //  fEMCALRecoUtils->SetTracksMatchedToCluster(esdevent) ; 
1319     //}
1320   }
1321 }