added run dependent corrections on T for calibration, removed old method
[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 #include <TFile.h>
32
33 // --- ANALYSIS system ---
34 #include "AliCalorimeterUtils.h"
35 #include "AliESDEvent.h"
36 #include "AliMCEvent.h"
37 #include "AliStack.h"
38 #include "AliAODPWG4Particle.h"
39 #include "AliVCluster.h"
40 #include "AliVCaloCells.h"
41 #include "AliMixedEvent.h"
42 #include "AliAODCaloCluster.h"
43 #include "AliOADBContainer.h"
44 #include "AliAnalysisManager.h"
45
46 // --- Detector ---
47 #include "AliEMCALGeometry.h"
48 #include "AliPHOSGeoUtils.h"
49
50 ClassImp(AliCalorimeterUtils)
51   
52   
53 //____________________________________________
54   AliCalorimeterUtils::AliCalorimeterUtils() : 
55     TObject(), fDebug(0), 
56     fEMCALGeoName(""),
57     fPHOSGeoName (""), 
58     fEMCALGeo(0x0),                   fPHOSGeo(0x0), 
59     fEMCALGeoMatrixSet(kFALSE),       fPHOSGeoMatrixSet(kFALSE), 
60     fLoadEMCALMatrices(kFALSE),       fLoadPHOSMatrices(kFALSE),
61     fRemoveBadChannels(kFALSE),       fPHOSBadChannelMap(0x0), 
62     fNCellsFromPHOSBorder(0),
63     fNMaskCellColumns(0),             fMaskCellColumns(0x0),
64     fRecalibration(kFALSE),           fRunDependentCorrection(kFALSE), fPHOSRecalibrationFactors(),
65     fEMCALRecoUtils(new AliEMCALRecoUtils),
66     fRecalculatePosition(kFALSE),     fCorrectELinearity(kFALSE),
67     fRecalculateMatching(kFALSE),
68     fCutR(20),                        fCutZ(20),
69     fCutEta(20),                      fCutPhi(20),
70     fLocMaxCutE(0),                   fLocMaxCutEDiff(0),
71     fPlotCluster(0),                  fOADBSet(kFALSE),
72     fOADBForEMCAL(kFALSE),            fOADBForPHOS(kFALSE),
73     fOADBFilePathEMCAL(""),           fOADBFilePathPHOS("")
74 {
75   //Ctor
76   
77   //Initialize parameters
78   InitParameters();
79   for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
80   for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix [i] = 0 ;
81   
82 }
83
84 //_________________________________________
85 AliCalorimeterUtils::~AliCalorimeterUtils()
86 {
87   //Dtor
88   
89   //if(fPHOSGeo)  delete fPHOSGeo  ;
90   if(fEMCALGeo) delete fEMCALGeo ;
91   
92   if(fPHOSBadChannelMap) { 
93     fPHOSBadChannelMap->Clear();
94     delete  fPHOSBadChannelMap;
95   }
96         
97   if(fPHOSRecalibrationFactors) { 
98     fPHOSRecalibrationFactors->Clear();
99     delete  fPHOSRecalibrationFactors;
100   }
101         
102   if(fEMCALRecoUtils)   delete fEMCALRecoUtils ;
103   if(fNMaskCellColumns) delete [] fMaskCellColumns;
104   
105 }
106
107 //____________________________________________________
108 void AliCalorimeterUtils::AccessOADB(AliVEvent* event)
109 {
110   // Set the AODB calibration, bad channels etc. parameters at least once
111   // alignment matrices from OADB done in SetGeometryMatrices
112   
113   //Set it only once
114   if(fOADBSet) return ; 
115   
116   Int_t   runnumber = event->GetRunNumber() ;
117   TString pass      = GetPass();
118   
119   // EMCAL
120   if(fOADBForEMCAL)
121   {
122     printf("AliCalorimeterUtils::SetOADBParameters() - Get AODB parameters from EMCAL in %s for run %d, and <%s> \n",fOADBFilePathEMCAL.Data(),runnumber,pass.Data());
123     
124     Int_t nSM = fEMCALGeo->GetNumberOfSuperModules();
125     
126     // Bad map
127     if(fRemoveBadChannels)
128     {
129       AliOADBContainer *contBC=new AliOADBContainer("");
130       contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePathEMCAL.Data()),"AliEMCALBadChannels"); 
131       
132       TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber);
133       
134       if(arrayBC)
135       {
136         SwitchOnDistToBadChannelRecalculation();
137         printf("AliCalorimeterUtils::SetOADBParameters() - Remove EMCAL bad cells \n");
138         
139         for (Int_t i=0; i<nSM; ++i) 
140         {
141           TH2I *hbm = GetEMCALChannelStatusMap(i);
142           
143           if (hbm)
144             delete hbm;
145           
146           hbm=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
147           
148           if (!hbm) 
149           {
150             AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
151             continue;
152           }
153           
154           hbm->SetDirectory(0);
155           SetEMCALChannelStatusMap(i,hbm);
156           
157         } // loop
158       } else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT remove EMCAL bad channels\n"); // run array
159     }  // Remove bad
160     
161     // Energy Recalibration
162     if(fRecalibration)
163     {
164       AliOADBContainer *contRF=new AliOADBContainer("");
165       
166       contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALRecalib");
167       
168       TObjArray *recal=(TObjArray*)contRF->GetObject(runnumber); 
169       
170       if(recal)
171       {
172         TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
173         
174         if(recalpass)
175         {
176           TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
177           
178           if(recalib)
179           {
180             printf("AliCalorimeterUtils::SetOADBParameters() - Recalibrate EMCAL \n");
181             for (Int_t i=0; i<nSM; ++i) 
182             {
183               TH2F *h = GetEMCALChannelRecalibrationFactors(i);
184               
185               if (h)
186                 delete h;
187               
188               h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
189               
190               if (!h) 
191               {
192                 AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
193                 continue;
194               }
195               
196               h->SetDirectory(0);
197               
198               SetEMCALChannelRecalibrationFactors(i,h);
199             } // SM loop
200           }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL, no params object array\n"); // array ok
201         }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for pass\n"); // array pass ok
202       }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for run\n");  // run number array ok
203       
204       // once set, apply run dependent corrections if requested
205       //fEMCALRecoUtils->SetRunDependentCorrections(runnumber);
206             
207     } // Recalibration on
208     
209     // Energy Recalibration, apply on top of previous calibration factors
210     if(fRunDependentCorrection)
211     {
212       AliOADBContainer *contRFTD=new AliOADBContainer("");
213       
214       contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALRunDepTempCalibCorrections");
215       
216       TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber); 
217       
218       if(htd)
219       {
220         printf("AliCalorimeterUtils::SetOADBParameters() - Recalibrate (Temperature) EMCAL \n");
221         
222         for (Int_t ism=0; ism<nSM; ++ism) 
223         {        
224           for (Int_t icol=0; icol<48; ++icol) 
225           {        
226             for (Int_t irow=0; irow<24; ++irow) 
227             {
228               Float_t factor = GetEMCALChannelRecalibrationFactor(ism,icol,irow);
229               
230               Int_t absID = fEMCALGeo->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
231               factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
232               //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID, 
233               //      GetEMCALChannelRecalibrationFactor(ism,icol,irow) , htd->GetBinContent(absID) / 10000., factor);
234               SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
235             } // columns
236           } // rows 
237         } // SM loop
238       }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL with T variations, no params TH1 \n"); 
239     } // Run by Run T calibration    
240     
241     // Time Recalibration
242     if(fEMCALRecoUtils->IsTimeRecalibrationOn())
243     {
244       AliOADBContainer *contTRF=new AliOADBContainer("");
245       
246       contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALTimeCalib");
247       
248       TObjArray *trecal=(TObjArray*)contTRF->GetObject(runnumber); 
249       
250       if(trecal)
251       {
252         TObjArray *trecalpass=(TObjArray*)trecal->FindObject(pass);
253         
254         if(trecalpass)
255         {
256           printf("AliCalorimeterUtils::SetOADBParameters() - Time Recalibrate EMCAL \n");
257           for (Int_t ibc = 0; ibc < 4; ++ibc) 
258           {
259             TH1F *h = GetEMCALChannelTimeRecalibrationFactors(ibc);
260             
261             if (h)
262               delete h;
263             
264             h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
265             
266             if (!h) 
267             {
268               AliError(Form("Could not load hAllTimeAvBC%d",ibc));
269               continue;
270             }
271             
272             h->SetDirectory(0);
273             
274             SetEMCALChannelTimeRecalibrationFactors(ibc,h);
275           } // bunch crossing loop
276         }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for pass\n"); // array pass ok
277       }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for run\n");  // run number array ok
278       
279     } // Recalibration on    
280     
281   }// EMCAL
282   
283   // PHOS
284   if(fOADBForPHOS)
285   {
286     printf("AliCalorimeterUtils::SetOADBParameters() - Get AODB parameters from PHOS in %s for run %d, and <%s> \n",fOADBFilePathPHOS.Data(),runnumber,pass.Data());
287     
288     // Bad map
289     if(fRemoveBadChannels)
290     {
291       AliOADBContainer badmapContainer(Form("phosBadMap"));
292       TString fileName="$ALICE_ROOT/OADB/PHOS/PHOSBadMaps.root";
293       badmapContainer.InitFromFile(Form("%s/PHOSBadMaps.root",fOADBFilePathPHOS.Data()),"phosBadMap");
294       
295       //Use a fixed run number from year 2010, this year not available yet.
296       TObjArray *maps = (TObjArray*)badmapContainer.GetObject(139000,"phosBadMap");
297       if(!maps)
298       {
299         printf("AliCalorimeterUtils::SetOADBParameters() - Can not read PHOS bad map for run %d.\n",runnumber) ;    
300       }
301       else
302       {
303         printf("AliCalorimeterUtils::SetOADBParameters() - Setting PHOS bad map with name %s \n",maps->GetName()) ;
304         for(Int_t mod=1; mod<5;mod++)
305         {
306           TH2I *hbmPH = GetPHOSChannelStatusMap(mod);
307           
308           if(hbmPH) 
309             delete hbmPH ;  
310           
311           hbmPH = (TH2I*)maps->At(mod);
312           
313           if(hbmPH) hbmPH->SetDirectory(0);
314           
315           SetPHOSChannelStatusMap(mod-1,hbmPH);
316           
317         } // modules loop  
318       } // maps exist
319     } // Remove bad channels
320   } // PHOS
321   
322   // Parameters already set once, so do not it again
323   fOADBSet = kTRUE;
324   
325 }  
326
327 //_____________________________________________________________
328 void AliCalorimeterUtils::AccessGeometry(AliVEvent* inputEvent) 
329 {
330   //Set the calorimeters transformation matrices and init geometry
331   
332   // First init the geometry, a priory not done before
333   Int_t runnumber = inputEvent->GetRunNumber() ;
334   InitPHOSGeometry (runnumber);
335   InitEMCALGeometry(runnumber);
336   
337   //Get the EMCAL transformation geometry matrices from ESD 
338   if(!fEMCALGeoMatrixSet && fEMCALGeo)
339   {
340     if(fLoadEMCALMatrices)
341     {
342       printf("AliCalorimeterUtils::AccessGeometry() - Load user defined EMCAL geometry matrices\n");
343       
344       // OADB if available
345       AliOADBContainer emcGeoMat("AliEMCALgeo");
346       emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePathEMCAL.Data()),"AliEMCALgeo");
347       TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
348       
349       for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
350       {
351         if (!fEMCALMatrix[mod]) // Get it from OADB
352         {
353           if(fDebug > 1 ) 
354             printf("AliCalorimeterUtils::AccessGeometry() - EMCAL matrices SM %d, %p\n",
355                    mod,((TGeoHMatrix*) matEMCAL->At(mod)));
356           //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
357           
358           fEMCALMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
359         }
360         
361         if(fEMCALMatrix[mod])
362         {
363           if(fDebug > 1) 
364             fEMCALMatrix[mod]->Print();
365           
366           fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;  
367         }
368       }//SM loop
369       
370       fEMCALGeoMatrixSet = kTRUE;//At least one, so good
371       
372     }//Load matrices
373     else if (!gGeoManager) 
374     { 
375       if(fDebug > 1) 
376         printf(" AliCalorimeterUtils::AccessGeometry() - Load EMCAL misalignment matrices. \n");
377       if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  
378       {
379         for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
380         { 
381           //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
382           if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod)) 
383           {
384             fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
385           }
386         }// loop over super modules     
387         
388         fEMCALGeoMatrixSet = kTRUE;//At least one, so good
389         
390       }//ESD as input
391       else 
392       {
393         if(fDebug > 1)
394           printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
395       }//AOD as input
396     }//Get matrix from data
397     else if(gGeoManager)
398     {
399       fEMCALGeoMatrixSet = kTRUE;
400     }
401   }//EMCAL geo && no geoManager
402   
403         //Get the PHOS transformation geometry matrices from ESD 
404   if(!fPHOSGeoMatrixSet && fPHOSGeo)
405   {
406     if(fLoadPHOSMatrices)
407     {
408       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined PHOS geometry matrices\n");
409       
410       // OADB if available
411       AliOADBContainer geomContainer("phosGeo");
412       geomContainer.InitFromFile(Form("%s/PHOSGeometry.root",fOADBFilePathPHOS.Data()),"PHOSRotationMatrixes");
413       TObjArray *matPHOS = (TObjArray*)geomContainer.GetObject(139000,"PHOSRotationMatrixes");    
414       
415       for(Int_t mod = 0 ; mod < 5 ; mod++)
416       {
417         if (!fPHOSMatrix[mod]) // Get it from OADB
418         {
419           if(fDebug > 1 ) 
420             printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - PHOS matrices module %d, %p\n",
421                    mod,((TGeoHMatrix*) matPHOS->At(mod)));
422           //((TGeoHMatrix*) matPHOS->At(mod))->Print();
423           
424           fPHOSMatrix[mod] = (TGeoHMatrix*) matPHOS->At(mod) ;
425         }
426         
427         // Set it, if it exists
428         if(fPHOSMatrix[mod])
429         {
430           if(fDebug > 1 ) 
431             fPHOSMatrix[mod]->Print();
432           
433           fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;  
434         }      
435       }// SM loop
436       
437       fPHOSGeoMatrixSet = kTRUE;//At least one, so good
438       
439     }//Load matrices
440     else if (!gGeoManager) 
441     { 
442       if(fDebug > 1) 
443         printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
444                         if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  
445       {
446                                 for(Int_t mod = 0; mod < 5; mod++)
447         { 
448                                         if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod)) 
449           {
450                                                 //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
451                                                 fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
452                                         }
453                                 }// loop over modules
454         fPHOSGeoMatrixSet  = kTRUE; //At least one so good
455                         }//ESD as input
456                         else 
457       {
458                                 if(fDebug > 1) 
459                                         printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
460       }//AOD as input
461     }// get matrix from data
462     else if(gGeoManager)
463     {
464       fPHOSGeoMatrixSet = kTRUE;
465     }
466         }//PHOS geo     and  geoManager was not set
467   
468 }
469
470 //______________________________________________________________________________________
471 Bool_t AliCalorimeterUtils::AreNeighbours(const TString calo, 
472                                           const Int_t absId1, const Int_t absId2 ) const
473 {
474   // Tells if (true) or not (false) two cells are neighbours
475   // A neighbour is defined as being two cells which share a side or corner
476         
477   Bool_t areNeighbours = kFALSE ;
478   
479   Int_t iRCU1 = -1, irow1 = -1, icol1 = -1;
480   Int_t iRCU2 = -1, irow2 = -1, icol2 = -1;
481   
482   Int_t rowdiff =  0, coldiff =  0;
483   
484   Int_t nSupMod1 = GetModuleNumberCellIndexes(absId1, calo, icol1, irow1, iRCU1); 
485   Int_t nSupMod2 = GetModuleNumberCellIndexes(absId2, calo, icol2, irow2, iRCU2); 
486   
487   if(calo=="EMCAL" && nSupMod1!=nSupMod2)
488   {
489     // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
490     // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
491     if(nSupMod1%2) icol1+=AliEMCALGeoParams::fgkEMCALCols;
492     else           icol2+=AliEMCALGeoParams::fgkEMCALCols;    
493         }
494   
495   rowdiff = TMath::Abs( irow1 - irow2 ) ;  
496   coldiff = TMath::Abs( icol1 - icol2 ) ;  
497   
498   if (( coldiff <= 1 )  && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0)) 
499     areNeighbours = kTRUE ;
500   
501   return areNeighbours;
502   
503 }
504
505
506 //_____________________________________________________________________________________
507 Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster, 
508                                                     AliVCaloCells* cells, 
509                                                     AliVEvent * event, Int_t iev) const 
510 {
511   
512         // Given the list of AbsId of the cluster, get the maximum cell and 
513         // check if there are fNCellsFromBorder from the calorimeter border
514         
515   //If the distance to the border is 0 or negative just exit accept all clusters
516         if(cells->GetType()==AliVCaloCells::kEMCALCell && fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder() <= 0 ) return kTRUE;
517         if(cells->GetType()==AliVCaloCells::kPHOSCell  && fNCellsFromPHOSBorder  <= 0 ) return kTRUE;
518   
519   Int_t absIdMax        = -1;
520         Float_t ampMax  = -1;
521         
522   AliMixedEvent * mixEvent = dynamic_cast<AliMixedEvent*> (event);
523   Int_t nMixedEvents = 0 ; 
524   Int_t * cellsCumul = NULL ;
525   Int_t numberOfCells = 0 ;  
526   if (mixEvent){
527     nMixedEvents = mixEvent->GetNumberOfEvents() ; 
528     if (cells->GetType()==AliVCaloCells::kEMCALCell) {
529       cellsCumul =  mixEvent->GetEMCALCellsCumul() ; 
530       numberOfCells = mixEvent->GetNumberOfEMCALCells() ;
531     } 
532     
533     else if (cells->GetType()==AliVCaloCells::kPHOSCell) {
534       cellsCumul =  mixEvent->GetPHOSCellsCumul() ; 
535       numberOfCells = mixEvent->GetNumberOfPHOSCells() ;
536     } 
537     
538     if(cellsCumul){
539       
540       Int_t startCell = cellsCumul[iev] ; 
541       Int_t endCell   = (iev+1 < nMixedEvents)?cellsCumul[iev+1]:numberOfCells;
542       //Find cells with maximum amplitude
543       for(Int_t i = 0; i < cluster->GetNCells() ; i++){
544         Int_t absId = cluster->GetCellAbsId(i) ;
545         for (Int_t j = startCell; j < endCell ;  j++) {
546           Short_t cellNumber, mclabel; 
547           Double_t amp, time, efrac; 
548           cells->GetCell(j, cellNumber, amp, time,mclabel,efrac) ; 
549           if (absId == cellNumber) {
550             if(amp > ampMax){
551               ampMax   = amp;
552               absIdMax = absId;
553             }        
554           }
555         }
556       }//loop on cluster cells
557     }// cells cumul available
558     else {
559       printf("AliCalorimeterUtils::CheckCellFiducialRegion() - CellsCumul is NULL!!!\n");
560       abort();
561     }
562   } else {//Normal SE Events
563     for(Int_t i = 0; i < cluster->GetNCells() ; i++){
564       Int_t absId = cluster->GetCellAbsId(i) ;
565       Float_t amp       = cells->GetCellAmplitude(absId);
566       if(amp > ampMax){
567         ampMax   = amp;
568         absIdMax = absId;
569       }
570     }
571   }
572         
573         if(fDebug > 1)
574                 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n", 
575            absIdMax, ampMax, cluster->E());
576         
577         if(absIdMax==-1) return kFALSE;
578         
579         //Check if the cell is close to the borders:
580         Bool_t okrow = kFALSE;
581         Bool_t okcol = kFALSE;
582   
583         if(cells->GetType()==AliVCaloCells::kEMCALCell){
584     
585                 Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1; 
586                 fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta); 
587                 fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
588                 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
589       Fatal("CheckCellFidutialRegion","Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",iSM,ieta,iphi);
590     }
591     
592                 //Check rows/phi
593     Int_t nborder = fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder();
594                 if(iSM < 10)
595     {
596                         if(iphi >= nborder && iphi < 24-nborder) okrow =kTRUE; 
597     }
598                 else
599     {
600       if(fEMCALGeoName.Contains("12SM")) // 1/3 SM
601       {
602         if(iphi >= nborder && iphi < 8-nborder) okrow =kTRUE; 
603       }
604       else // 1/2 SM
605       {
606         if(iphi >= nborder && iphi <12-nborder) okrow =kTRUE; 
607       }
608                 }
609                 
610                 //Check columns/eta
611                 if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0())
612     {
613                         if(ieta  > nborder && ieta < 48-nborder) okcol =kTRUE; 
614                 }
615                 else
616     {
617                         if(iSM%2==0)
618       {
619                                 if(ieta >= nborder)     okcol = kTRUE;  
620                         }
621                         else 
622       {
623                                 if(ieta <  48-nborder)  okcol = kTRUE;  
624                         }
625                 }//eta 0 not checked
626     
627                 if(fDebug > 1)
628                 {
629                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
630              nborder, ieta, iphi, iSM);
631                         if (okcol && okrow ) printf(" YES \n");
632                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
633                 }
634         }//EMCAL
635         else if(cells->GetType()==AliVCaloCells::kPHOSCell){
636                 Int_t relId[4];
637                 Int_t irow = -1, icol = -1;
638                 fPHOSGeo->AbsToRelNumbering(absIdMax,relId);
639                 irow = relId[2];
640                 icol = relId[3];
641                 //imod = relId[0]-1;
642                 if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE; 
643                 if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE; 
644                 if(fDebug > 1)
645                 {
646                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - PHOS Cluster in %d cells fiducial volume: icol %d, irow %d, Module %d?",
647              fNCellsFromPHOSBorder, icol, irow, relId[0]-1);
648                         if (okcol && okrow ) printf(" YES \n");
649                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
650                 }
651         }//PHOS
652         
653         if (okcol && okrow) return kTRUE; 
654         else                return kFALSE;
655         
656 }       
657
658 //_________________________________________________________________________________________________________
659 Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter,UShort_t* cellList, Int_t nCells)
660 {
661         // Check that in the cluster cells, there is no bad channel of those stored 
662         // in fEMCALBadChannelMap or fPHOSBadChannelMap
663         
664         if (!fRemoveBadChannels) return kFALSE;
665         //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
666         if(calorimeter == "EMCAL" && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
667         if(calorimeter == "PHOS"  && !fPHOSBadChannelMap)  return kFALSE;
668   
669         Int_t icol = -1;
670         Int_t irow = -1;
671         Int_t imod = -1;
672         for(Int_t iCell = 0; iCell<nCells; iCell++){
673     
674                 //Get the column and row
675                 if(calorimeter == "EMCAL"){
676       return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
677                 }
678                 else if(calorimeter=="PHOS"){
679                         Int_t    relId[4];
680                         fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
681                         irow = relId[2];
682                         icol = relId[3];
683                         imod = relId[0]-1;
684                         if(fPHOSBadChannelMap->GetEntries() <= imod)continue;
685       //printf("PHOS bad channels imod %d, icol %d, irow %d\n",imod, irow, icol);
686                         if(GetPHOSChannelStatus(imod, irow, icol)) return kTRUE;
687                 }
688                 else return kFALSE;
689                 
690         }// cell cluster loop
691   
692         return kFALSE;
693   
694 }
695
696 //_______________________________________________________________
697 void AliCalorimeterUtils::CorrectClusterEnergy(AliVCluster *clus)
698 {
699   // Correct cluster energy non linearity
700   
701   clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
702
703 }
704
705 //________________________________________________________________________________________
706 Int_t  AliCalorimeterUtils::GetMaxEnergyCell(AliVCaloCells* cells, const AliVCluster* clu, 
707                                              Float_t & clusterFraction) const 
708 {
709   
710   //For a given CaloCluster gets the absId of the cell 
711   //with maximum energy deposit.
712   
713   if( !clu || !cells ){
714     AliInfo("Cluster or cells pointer is null!");
715     return -1;
716   }
717   
718   Double_t eMax        =-1.;
719   Double_t eTot        = 0.;
720   Double_t eCell       =-1.;
721   Float_t  fraction    = 1.;
722   Float_t  recalFactor = 1.;
723   Int_t    cellAbsId   =-1 , absId =-1 ;
724   Int_t    iSupMod     =-1 , ieta  =-1 , iphi = -1, iRCU = -1;
725   
726   TString           calo = "EMCAL";
727   if(clu->IsPHOS()) calo = "PHOS";
728   
729   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
730     
731     cellAbsId = clu->GetCellAbsId(iDig);
732     
733     fraction  = clu->GetCellAmplitudeFraction(iDig);
734     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
735     
736     iSupMod = GetModuleNumberCellIndexes(cellAbsId, calo, ieta, iphi, iRCU);
737     
738     if(IsRecalibrationOn()) {
739       if(calo=="EMCAL") recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
740       else              recalFactor = GetPHOSChannelRecalibrationFactor (iSupMod,iphi,ieta);
741     }
742     
743     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
744     
745     if(eCell > eMax)  { 
746       eMax  = eCell; 
747       absId = cellAbsId;
748     }
749     
750     eTot+=eCell;
751     
752   }// cell loop
753   
754   clusterFraction = (eTot-eMax)/eTot; //Do not use cluster energy in case it was corrected for non linearity.
755   
756   return absId;
757   
758 }
759
760 //__________________________________________________________________________
761 AliVTrack * AliCalorimeterUtils::GetMatchedTrack(const AliVCluster* cluster, 
762                                                  const AliVEvent* event, 
763                                                  const Int_t index) const
764 {
765   // Get the matched track given its index, usually just the first match
766   // Since it is different for ESDs and AODs here it is a wrap method to do it
767   
768   AliVTrack *track = 0;
769   
770   // EMCAL case only when matching is recalculated
771   if(cluster->IsEMCAL() && IsRecalculationOfClusterTrackMatchingOn())
772   {
773     Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(cluster->GetID());
774     //printf("track index %d, cluster ID %d \n ",trackIndex,cluster->GetID());
775     
776     if(trackIndex < 0 )
777     { 
778       printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d, from recalculation\n", trackIndex);
779     }
780     else 
781     {
782       track = dynamic_cast<AliVTrack*> (event->GetTrack(trackIndex));
783     }
784
785     return track ;
786     
787   }   
788   
789   // Normal case, get info from ESD or AOD
790   // ESDs
791   if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
792   {
793     Int_t iESDtrack = cluster->GetTrackMatchedIndex();
794     
795     if(iESDtrack < 0 )
796     { 
797       printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d\n", index);
798       return 0x0;
799     }
800     
801     track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
802     
803   }
804   else // AODs
805   {
806     if(cluster->GetNTracksMatched() > 0 )
807       track = dynamic_cast<AliVTrack*>(cluster->GetTrackMatched(index));
808   }
809   
810   return track ;
811   
812 }
813
814 //_____________________________________________________________________________________________________
815 Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
816 {
817         //Get the EMCAL/PHOS module number that corresponds to this particle
818         
819         Int_t absId = -1;
820         if(particle->GetDetector()=="EMCAL"){
821                 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
822                 if(fDebug > 2) 
823                   printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
824              particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
825                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
826         }//EMCAL
827         else if(particle->GetDetector()=="PHOS")
828   {
829     // In case we use the MC reader, the input are TParticles, 
830     // in this case use the corresponing method in PHOS Geometry to get the particle.
831     if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
832     {
833       Int_t mod =-1;
834       Double_t z = 0., x=0.;
835       TParticle* primary = 0x0;
836       AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
837       if(stack) {
838         primary = stack->Particle(particle->GetCaloLabel(0));
839       }
840       else {
841         Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
842       }
843       
844       if(primary){
845         fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
846       }
847       else{
848         Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
849       }
850       return mod;
851     }
852     // Input are ESDs or AODs, get the PHOS module number like this.
853     else{
854       //FIXME
855       //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
856       //return GetModuleNumber(cluster);
857       //MEFIX
858       return -1;
859     }
860         }//PHOS
861         
862         return -1;
863 }
864
865 //_____________________________________________________________________
866 Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
867 {
868         //Get the EMCAL/PHOS module number that corresponds to this cluster
869         TLorentzVector lv;
870         Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
871   if(!cluster)
872   {
873     if(fDebug > 1) printf("AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
874     return -1;
875   }
876   
877         cluster->GetMomentum(lv,v);
878         Float_t phi = lv.Phi();
879         if(phi < 0) phi+=TMath::TwoPi();        
880         Int_t absId = -1;
881         if(cluster->IsEMCAL()){
882                 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
883                 if(fDebug > 2) 
884                   printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
885              lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
886                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
887         }//EMCAL
888         else if(cluster->IsPHOS()) 
889   {
890                 Int_t    relId[4];
891                 if ( cluster->GetNCells() > 0) 
892     {
893                         absId = cluster->GetCellAbsId(0);
894                         if(fDebug > 2) 
895                                 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
896                lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
897                 }
898                 else return -1;
899                 
900                 if ( absId >= 0) 
901     {
902                         fPHOSGeo->AbsToRelNumbering(absId,relId);
903                         if(fDebug > 2) 
904                           printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
905                         return relId[0]-1;
906                 }
907                 else return -1;
908         }//PHOS
909         
910         return -1;
911 }
912
913 //___________________________________________________________________________________________________
914 Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(const Int_t absId, const TString calo, 
915                                                       Int_t & icol, Int_t & irow, Int_t & iRCU) const
916 {
917         //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
918         Int_t imod = -1;
919         if ( absId >= 0) 
920   {
921                 if(calo=="EMCAL")
922     {
923                         Int_t iTower = -1, iIphi = -1, iIeta = -1; 
924                         fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
925                         fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
926       if(imod < 0 || irow < 0 || icol < 0 ) 
927       {
928         Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
929       }
930       
931                         //RCU0
932       if(imod < 10 )
933       {
934         if      (0<=irow&&irow<8)                       iRCU=0; // first cable row
935         else if (8<=irow&&irow<16 &&  0<=icol&&icol<24) iRCU=0; // first half; 
936         //second cable row
937         //RCU1
938         else if (8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half; 
939         //second cable row
940         else if (16<=irow&&irow<24)                     iRCU=1; // third cable row
941         
942         if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
943       }
944       else 
945       {
946         // Last 2 SM have one single SRU, just assign RCU 0
947         iRCU = 0 ;
948       }
949
950                         if (iRCU<0) 
951       {
952                                 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
953                         }                       
954                         
955                         return imod ;
956       
957                 }//EMCAL
958                 else //PHOS
959     {
960                         Int_t    relId[4];
961                         fPHOSGeo->AbsToRelNumbering(absId,relId);
962                         irow = relId[2];
963                         icol = relId[3];
964                         imod = relId[0]-1;
965                         iRCU= (Int_t)(relId[2]-1)/16 ;
966                         //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
967                         if (iRCU >= 4) 
968       {
969                                 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
970                         }                       
971                         return imod;
972                 }//PHOS 
973         }
974         
975         return -1;
976 }
977
978 //___________________________________________________________________________________________
979 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells) 
980 {
981   // Find local maxima in cluster
982   
983   const Int_t   nc = cluster->GetNCells();
984   Int_t   absIdList[nc]; 
985   Float_t maxEList[nc]; 
986   
987   Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList);
988   
989   return nMax;
990   
991 }
992
993 //___________________________________________________________________________________________
994 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
995                                                   Int_t *absIdList,     Float_t *maxEList) 
996 {
997   // Find local maxima in cluster
998   
999   Int_t iDigitN = 0 ;
1000   Int_t iDigit  = 0 ;
1001   Int_t absId1 = -1 ;
1002   Int_t absId2 = -1 ;
1003   const Int_t nCells = cluster->GetNCells();
1004   
1005   TString calorimeter = "EMCAL";
1006   if(!cluster->IsEMCAL()) calorimeter = "PHOS";
1007   
1008   //printf("cluster : ncells %d \n",nCells);
1009   
1010   Float_t emax  = 0;
1011   Int_t   idmax =-1;
1012   for(iDigit = 0; iDigit < nCells ; iDigit++)
1013   {
1014     absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit]  ; 
1015     Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1016     RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);  
1017     if( en > emax )
1018     {
1019       emax  = en ;
1020       idmax = absIdList[iDigit] ;
1021     }
1022     //Int_t icol = -1, irow = -1, iRCU = -1;
1023     //Int_t sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1024     //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
1025   }
1026   
1027   for(iDigit = 0 ; iDigit < nCells; iDigit++) 
1028   {   
1029     if(absIdList[iDigit]>=0) 
1030     {
1031       absId1 = cluster->GetCellsAbsId()[iDigit];
1032       
1033       Float_t en1 = cells->GetCellAmplitude(absId1);
1034       RecalibrateCellAmplitude(en1,calorimeter,absId1);  
1035       
1036       //printf("%d : absIDi %d, E %f\n",iDigit, absId1,en1);
1037       
1038       for(iDigitN = 0; iDigitN < nCells; iDigitN++) 
1039       { 
1040         absId2 = cluster->GetCellsAbsId()[iDigitN] ;
1041         
1042         if(absId2==-1 || absId2==absId1) continue;
1043         
1044         //printf("\t %d : absIDj %d\n",iDigitN, absId2);
1045         
1046         Float_t en2 = cells->GetCellAmplitude(absId2);
1047         RecalibrateCellAmplitude(en2,calorimeter,absId2);
1048         
1049         //printf("\t %d : absIDj %d, E %f\n",iDigitN, absId2,en2);
1050         
1051         if ( AreNeighbours(calorimeter, absId1, absId2) ) 
1052         {
1053           // printf("\t \t Neighbours \n");
1054           if (en1 > en2 ) 
1055           {    
1056             absIdList[iDigitN] = -1 ;
1057             //printf("\t \t indexN %d not local max\n",iDigitN);
1058             // but may be digit too is not local max ?
1059             if(en1 < en2 + fLocMaxCutEDiff) {
1060               //printf("\t \t index %d not local max cause locMaxCutEDiff\n",iDigit);
1061               absIdList[iDigit] = -1 ;
1062             }
1063           }
1064           else 
1065           {
1066             absIdList[iDigit] = -1 ;
1067             //printf("\t \t index %d not local max\n",iDigitN);
1068             // but may be digitN too is not local max ?
1069             if(en1 > en2 - fLocMaxCutEDiff) 
1070             {
1071               absIdList[iDigitN] = -1 ; 
1072               //printf("\t \t indexN %d not local max cause locMaxCutEDiff\n",iDigit);
1073             }
1074           } 
1075         } // if Are neighbours
1076         //else printf("\t \t NOT Neighbours \n");
1077       } // while digitN
1078     } // slot not empty
1079   } // while digit
1080   
1081   iDigitN = 0 ;
1082   for(iDigit = 0; iDigit < nCells; iDigit++) 
1083   { 
1084     if(absIdList[iDigit]>=0 )
1085     {
1086       absIdList[iDigitN] = absIdList[iDigit] ;
1087       Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1088       RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);  
1089       if(en < fLocMaxCutE) continue; // Maxima only with seed energy at least
1090       maxEList[iDigitN] = en ;
1091       //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
1092       iDigitN++ ; 
1093     }
1094   }
1095   
1096   if(iDigitN == 0)
1097   {
1098     if(fDebug > 0) 
1099       printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - No local maxima found, assign highest energy cell as maxima, id %d, en cell %2.2f, en cluster %2.2f\n",
1100              idmax,emax,cluster->E());
1101     iDigitN      = 1     ;
1102     maxEList[0]  = emax  ;
1103     absIdList[0] = idmax ; 
1104   }
1105   
1106   if(fDebug > 0) 
1107   {    
1108     printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - In cluster E %2.2f, M02 %2.2f, M20 %2.2f, N maxima %d \n", 
1109            cluster->E(),cluster->GetM02(),cluster->GetM20(), iDigitN);
1110   
1111     if(fDebug > 1) for(Int_t imax = 0; imax < iDigitN; imax++) 
1112     {
1113       printf(" \t i %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
1114     }
1115   }
1116   
1117   return iDigitN ;
1118   
1119 }
1120
1121 //____________________________________
1122 TString AliCalorimeterUtils::GetPass()
1123 {
1124   // Get passx from filename.
1125     
1126   if (!AliAnalysisManager::GetAnalysisManager()->GetTree()) 
1127   {
1128     AliError("AliCalorimeterUtils::GetPass() - Pointer to tree = 0, returning null\n");
1129     return TString("");
1130   }
1131   
1132   if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()) 
1133   {
1134     AliError("AliCalorimeterUtils::GetPass() - Null pointer input file, returning null\n");
1135     return TString("");
1136   }
1137   
1138   TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1139   if      (pass.Contains("ass1")) return TString("pass1");
1140   else if (pass.Contains("ass2")) return TString("pass2");
1141   else if (pass.Contains("ass3")) return TString("pass3");
1142
1143   // No condition fullfilled, give a default value
1144   printf("AliCalorimeterUtils::GetPass() - Pass number string not found \n");
1145   return TString("");            
1146   
1147 }
1148
1149 //________________________________________
1150 void AliCalorimeterUtils::InitParameters()
1151 {
1152   //Initialize the parameters of the analysis.
1153   
1154   fEMCALGeoName         = "";
1155   fPHOSGeoName          = "";
1156   
1157   fEMCALGeoMatrixSet    = kFALSE;
1158   fPHOSGeoMatrixSet     = kFALSE;
1159   
1160   fRemoveBadChannels    = kFALSE;
1161   
1162   fNCellsFromPHOSBorder = 0;
1163   
1164   fLocMaxCutE           = 0.1 ;
1165   fLocMaxCutEDiff       = 0.0 ;
1166
1167   //  fMaskCellColumns = new Int_t[fNMaskCellColumns];
1168   //  fMaskCellColumns[0] = 6 ;  fMaskCellColumns[1] = 7 ;  fMaskCellColumns[2] = 8 ; 
1169   //  fMaskCellColumns[3] = 35;  fMaskCellColumns[4] = 36;  fMaskCellColumns[5] = 37; 
1170   //  fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
1171   //  fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols; 
1172   //  fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols; 
1173   
1174   fOADBSet      = kFALSE;
1175   fOADBForEMCAL = kTRUE ;            
1176   fOADBForPHOS  = kFALSE;
1177   
1178   fOADBFilePathEMCAL = "$ALICE_ROOT/OADB/EMCAL" ;          
1179   fOADBFilePathPHOS  = "$ALICE_ROOT/OADB/PHOS"  ;     
1180   
1181 }
1182
1183
1184 //_____________________________________________________
1185 void AliCalorimeterUtils::InitPHOSBadChannelStatusMap()
1186 {
1187   //Init PHOS bad channels map
1188   if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
1189   //In order to avoid rewriting the same histograms
1190   Bool_t oldStatus = TH1::AddDirectoryStatus();
1191   TH1::AddDirectory(kFALSE);
1192   
1193   fPHOSBadChannelMap = new TObjArray(5);        
1194   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));
1195   
1196   fPHOSBadChannelMap->SetOwner(kTRUE);
1197   fPHOSBadChannelMap->Compress();
1198   
1199   //In order to avoid rewriting the same histograms
1200   TH1::AddDirectory(oldStatus);         
1201 }
1202
1203 //______________________________________________________
1204 void AliCalorimeterUtils::InitPHOSRecalibrationFactors()
1205 {
1206         //Init EMCAL recalibration factors
1207         if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
1208         //In order to avoid rewriting the same histograms
1209         Bool_t oldStatus = TH1::AddDirectoryStatus();
1210         TH1::AddDirectory(kFALSE);
1211   
1212         fPHOSRecalibrationFactors = new TObjArray(5);
1213         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));
1214         //Init the histograms with 1
1215         for (Int_t m = 0; m < 5; m++) {
1216                 for (Int_t i = 0; i < 56; i++) {
1217                         for (Int_t j = 0; j < 64; j++) {
1218                                 SetPHOSChannelRecalibrationFactor(m,j,i,1.);
1219                         }
1220                 }
1221         }
1222         fPHOSRecalibrationFactors->SetOwner(kTRUE);
1223         fPHOSRecalibrationFactors->Compress();
1224         
1225         //In order to avoid rewriting the same histograms
1226         TH1::AddDirectory(oldStatus);           
1227 }
1228
1229
1230 //__________________________________________________________
1231 void AliCalorimeterUtils::InitEMCALGeometry(Int_t runnumber)
1232 {
1233         //Initialize EMCAL geometry if it did not exist previously
1234   
1235         if (!fEMCALGeo)
1236   {
1237     if(fEMCALGeoName=="")
1238     {
1239       if     (runnumber <  140000 && 
1240               runnumber >= 100000)   fEMCALGeoName = "EMCAL_FIRSTYEARV1";
1241       else if(runnumber >= 140000 &&
1242               runnumber <  171000)   fEMCALGeoName = "EMCAL_COMPLETEV1";
1243       else                           fEMCALGeoName = "EMCAL_COMPLETE12SMV1";  
1244       printf("AliCalorimeterUtils::InitEMCALGeometry() - Set EMCAL geometry name to <%s> for run %d\n",fEMCALGeoName.Data(),runnumber);
1245     }
1246     
1247                 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
1248     
1249                 if(fDebug > 0)
1250     {
1251                         printf("AliCalorimeterUtils::InitEMCALGeometry(run=%d)",runnumber);
1252                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1253                         printf("\n");
1254                 }
1255         }
1256 }
1257
1258 //_________________________________________________________
1259 void AliCalorimeterUtils::InitPHOSGeometry(Int_t runnumber)
1260 {
1261         //Initialize PHOS geometry if it did not exist previously
1262   
1263         if (!fPHOSGeo)
1264   {
1265     if(fPHOSGeoName=="") fPHOSGeoName = "PHOSgeo";
1266       
1267                 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName); 
1268     
1269                 if(fDebug > 0)
1270     {
1271                         printf("AliCalorimeterUtils::InitPHOSGeometry(run=%d)",runnumber);
1272                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1273                         printf("\n");
1274                 }       
1275         }       
1276 }
1277
1278 //__________________________________________________________________
1279 Bool_t AliCalorimeterUtils::MaskFrameCluster(const Int_t iSM,  
1280                                              const Int_t ieta) const 
1281 {
1282   //Check if cell is in one of the regions where we have significant amount 
1283   //of material in front. Only EMCAL
1284   
1285   Int_t icol = ieta;
1286   if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
1287   
1288   if (fNMaskCellColumns && fMaskCellColumns) 
1289   {
1290     for (Int_t imask = 0; imask < fNMaskCellColumns; imask++) 
1291     {
1292       if(icol==fMaskCellColumns[imask]) return kTRUE;
1293     }
1294   }
1295   
1296   return kFALSE;
1297   
1298 }
1299
1300 //_________________________________________________________
1301 void AliCalorimeterUtils::Print(const Option_t * opt) const
1302 {
1303   
1304   //Print some relevant parameters set for the analysis
1305   if(! opt)
1306     return;
1307   
1308   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1309   printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
1310   printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
1311          fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
1312   if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
1313   printf("Recalibrate Clusters? %d, run by run  %d\n",fRecalibration,fRunDependentCorrection);
1314   printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
1315   printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
1316   printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1317   
1318   printf("Loc. Max. E > %2.2f\n",       fLocMaxCutE);
1319   printf("Loc. Max. E Diff > %2.2f\n",  fLocMaxCutEDiff);
1320   
1321   printf("    \n") ;
1322
1323
1324 //__________________________________________________________________________________________
1325 void AliCalorimeterUtils::RecalibrateCellAmplitude(Float_t & amp,
1326                                                    const TString calo, const Int_t id) const
1327 {
1328   //Recaculate cell energy if recalibration factor
1329   
1330   Int_t icol     = -1; Int_t irow     = -1; Int_t iRCU     = -1;
1331   Int_t nModule  = GetModuleNumberCellIndexes(id,calo, icol, irow, iRCU);
1332   
1333   if (IsRecalibrationOn()) 
1334   {
1335     if(calo == "PHOS") 
1336     {
1337       amp *= GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
1338     }
1339     else                                   
1340     {
1341       amp *= GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
1342     }
1343   }
1344 }
1345
1346 //_________________________________________________________________________________
1347 void AliCalorimeterUtils::RecalibrateCellTime(Double_t & time, 
1348                                               const TString calo, 
1349                                               const Int_t id, const Int_t bc) const
1350 {
1351   // Recalculate time if time recalibration available for EMCAL
1352   // not ready for PHOS
1353   
1354   if(calo == "EMCAL" && GetEMCALRecoUtils()->IsTimeRecalibrationOn()) 
1355   {
1356     GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
1357   }
1358   
1359 }
1360
1361 //__________________________________________________________________________
1362 Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster, 
1363                                                       AliVCaloCells * cells)
1364 {
1365         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
1366   
1367   //Initialize some used variables
1368         Float_t frac  = 0., energy = 0.;  
1369   
1370         if(cells) 
1371   {
1372     //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1373     
1374     UShort_t * index    = cluster->GetCellsAbsId() ;
1375     Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1376     
1377     Int_t ncells     = cluster->GetNCells();    
1378     
1379     TString calo     = "EMCAL";
1380     if(cluster->IsPHOS()) 
1381       calo = "PHOS";
1382     
1383     //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1384     for(Int_t icell = 0; icell < ncells; icell++){
1385       
1386       Int_t absId = index[icell];
1387       
1388       frac =  fraction[icell];
1389       if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
1390       
1391       Float_t amp = cells->GetCellAmplitude(absId);
1392       RecalibrateCellAmplitude(amp,calo, absId);
1393       
1394       if(fDebug>2)
1395         printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, cell fraction %f, cell energy %f\n", 
1396                calo.Data(),frac,cells->GetCellAmplitude(absId));
1397       
1398       energy += amp*frac;
1399     }
1400     
1401     if(fDebug>1)
1402       printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
1403     
1404         }// cells available
1405   else
1406   {
1407     Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
1408   }
1409   
1410         return energy;
1411 }
1412
1413 //__________________________________________________________________________________________
1414 void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
1415 {
1416   
1417   //Recalculate EMCAL cluster position
1418   
1419   fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
1420   
1421 }
1422
1423 //________________________________________________________________________________
1424 void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event, 
1425                                                           TObjArray* clusterArray) 
1426
1427   //Recalculate track matching
1428   
1429   if (fRecalculateMatching) 
1430   {
1431     fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo)   ; 
1432     //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1433     //if(esdevent){
1434     //  fEMCALRecoUtils->SetClusterMatchedToTrack(esdevent)  ;
1435     //  fEMCALRecoUtils->SetTracksMatchedToCluster(esdevent) ; 
1436     //}
1437   }
1438 }
1439
1440 //___________________________________________________________________________
1441 void AliCalorimeterUtils::SplitEnergy(const Int_t absId1, const Int_t absId2,
1442                                       AliVCluster* cluster, 
1443                                       AliVCaloCells* cells,
1444                                       //Float_t & e1, Float_t & e2,
1445                                       AliAODCaloCluster* cluster1,
1446                                       AliAODCaloCluster* cluster2,
1447                                       const Int_t nMax,
1448                                       const Int_t eventNumber)
1449 {
1450   
1451   // Split energy of cluster between the 2 local maxima, sum energy on 3x3, and if the 2 
1452   // maxima are too close and have common cells, split the energy between the 2
1453   
1454   TH2F* hClusterMap    = 0 ;
1455   TH2F* hClusterLocMax = 0 ;
1456   TH2F* hCluster1      = 0 ;
1457   TH2F* hCluster2      = 0 ;
1458   
1459   if(fPlotCluster)
1460   {
1461     hClusterMap    = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
1462     hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
1463     hCluster1      = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
1464     hCluster2      = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
1465     hClusterMap    ->SetXTitle("column");
1466     hClusterMap    ->SetYTitle("row");
1467     hClusterLocMax ->SetXTitle("column");
1468     hClusterLocMax ->SetYTitle("row");
1469     hCluster1      ->SetXTitle("column");
1470     hCluster1      ->SetYTitle("row");
1471     hCluster2      ->SetXTitle("column");
1472     hCluster2      ->SetYTitle("row");
1473   }
1474   
1475   TString calorimeter = "EMCAL";
1476   if(cluster->IsPHOS())
1477   {
1478     calorimeter="PHOS";
1479     printf("AliCalorimeterUtils::SplitEnerg() Not supported for PHOS yet \n");
1480     return;
1481   }
1482   
1483   const Int_t ncells  = cluster->GetNCells();  
1484   Int_t absIdList[ncells]; 
1485   
1486   Float_t e1 = 0,  e2   = 0 ;
1487   Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;  
1488   Float_t eCluster = 0;
1489   Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1; 
1490   for(Int_t iDigit  = 0; iDigit < ncells; iDigit++ ) 
1491   {
1492     absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
1493     
1494     Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
1495     RecalibrateCellAmplitude(ec,calorimeter, absIdList[iDigit]);
1496     eCluster+=ec;
1497     
1498     if(fPlotCluster) 
1499     {
1500       //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
1501       sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1502       if(icol > maxCol) maxCol = icol;
1503       if(icol < minCol) minCol = icol;
1504       if(irow > maxRow) maxRow = irow;
1505       if(irow < minRow) minRow = irow;
1506       hClusterMap->Fill(icol,irow,ec);
1507     }
1508     
1509   }
1510   
1511   // Init counters and variables
1512   Int_t ncells1 = 1 ;
1513   UShort_t absIdList1[9] ;  
1514   Double_t fracList1 [9] ;  
1515   absIdList1[0] = absId1 ;
1516   fracList1 [0] = 1. ;
1517   
1518   Float_t ecell1 = cells->GetCellAmplitude(absId1);
1519   RecalibrateCellAmplitude(ecell1, calorimeter, absId1);
1520   e1 =  ecell1;  
1521   
1522   Int_t ncells2 = 1 ;
1523   UShort_t absIdList2[9] ;  
1524   Double_t fracList2 [9] ; 
1525   absIdList2[0] = absId2 ;
1526   fracList2 [0] = 1. ;
1527   
1528   Float_t ecell2 = cells->GetCellAmplitude(absId2);
1529   RecalibrateCellAmplitude(ecell2, calorimeter, absId2);
1530   e2 =  ecell2;  
1531   
1532   if(fPlotCluster)
1533   {
1534     Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
1535     sm = GetModuleNumberCellIndexes(absId1, calorimeter, icol1, irow1, iRCU) ;
1536     hClusterLocMax->Fill(icol1,irow1,ecell1);
1537     sm = GetModuleNumberCellIndexes(absId2, calorimeter, icol2, irow2, iRCU) ;
1538     hClusterLocMax->Fill(icol2,irow2,ecell2);
1539   }
1540   
1541   // Very rough way to share the cluster energy
1542   Float_t eRemain = (eCluster-ecell1-ecell2)/2;
1543   Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
1544   Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
1545   
1546   for(Int_t iDigit = 0; iDigit < ncells; iDigit++)
1547   {
1548     Int_t absId = absIdList[iDigit];
1549     
1550     if(absId==absId1 || absId==absId2 || absId < 0) continue;
1551     
1552     Float_t ecell = cells->GetCellAmplitude(absId);
1553     RecalibrateCellAmplitude(ecell, calorimeter, absId);
1554     
1555     if(AreNeighbours(calorimeter, absId1,absId ))
1556     { 
1557       absIdList1[ncells1]= absId;
1558       
1559       if(AreNeighbours(calorimeter, absId2,absId ))
1560       { 
1561         fracList1[ncells1] = shareFraction1; 
1562         e1 += ecell*shareFraction1;
1563       }
1564       else 
1565       {
1566         fracList1[ncells1] = 1.; 
1567         e1 += ecell;
1568       }
1569       
1570       ncells1++;
1571       
1572     } // neigbour to cell1
1573     
1574     if(AreNeighbours(calorimeter, absId2,absId ))
1575     { 
1576       absIdList2[ncells2]= absId;
1577       
1578       if(AreNeighbours(calorimeter, absId1,absId ))
1579       { 
1580         fracList2[ncells2] = shareFraction2; 
1581         e2 += ecell*shareFraction2;
1582       }
1583       else
1584       { 
1585         fracList2[ncells2] = 1.; 
1586         e2 += ecell;
1587       }
1588       
1589       ncells2++;
1590       
1591     } // neigbour to cell2
1592     
1593   }
1594   
1595   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",
1596                             nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2);
1597   
1598   cluster1->SetE(e1);
1599   cluster2->SetE(e2);  
1600   
1601   cluster1->SetNCells(ncells1);
1602   cluster2->SetNCells(ncells2);  
1603   
1604   cluster1->SetCellsAbsId(absIdList1);
1605   cluster2->SetCellsAbsId(absIdList2);
1606   
1607   cluster1->SetCellsAmplitudeFraction(fracList1);
1608   cluster2->SetCellsAmplitudeFraction(fracList2);
1609   
1610   if(calorimeter=="EMCAL")
1611   {
1612     GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
1613     GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
1614   }
1615   
1616   if(fPlotCluster)
1617   {
1618     //printf("Cells of cluster1: ");
1619     for(Int_t iDigit  = 0; iDigit < ncells1; iDigit++ ) 
1620     {
1621       //printf(" %d ",absIdList1[iDigit]);
1622       
1623       sm = GetModuleNumberCellIndexes(absIdList1[iDigit], calorimeter, icol, irow, iRCU) ;
1624       
1625       if( AreNeighbours(calorimeter, absId2,absIdList1[iDigit]) )
1626         hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit])*shareFraction1);
1627       else 
1628         hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit]));
1629     }
1630     
1631     //printf(" \n ");
1632     //printf("Cells of cluster2: ");
1633     
1634     for(Int_t iDigit  = 0; iDigit < ncells2; iDigit++ ) 
1635     {
1636       //printf(" %d ",absIdList2[iDigit]);
1637       
1638       sm = GetModuleNumberCellIndexes(absIdList2[iDigit], calorimeter, icol, irow, iRCU) ;
1639       if( AreNeighbours(calorimeter, absId1,absIdList2[iDigit]) )
1640         hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit])*shareFraction2);
1641       else
1642         hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit]));
1643       
1644     }
1645     //printf(" \n ");
1646     
1647     gStyle->SetPadRightMargin(0.1);
1648     gStyle->SetPadLeftMargin(0.1);
1649     gStyle->SetOptStat(0);
1650     gStyle->SetOptFit(000000);
1651     
1652     if(maxCol-minCol > maxRow-minRow)
1653     {
1654       maxRow+= maxCol-minCol;
1655     }
1656     else 
1657     {
1658       maxCol+= maxRow-minRow;
1659     }
1660     
1661     TCanvas  * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
1662     c->Divide(2,2);  
1663     c->cd(1);
1664     gPad->SetGridy();
1665     gPad->SetGridx();
1666     hClusterMap    ->SetAxisRange(minCol, maxCol,"X");
1667     hClusterMap    ->SetAxisRange(minRow, maxRow,"Y");
1668     hClusterMap    ->Draw("colz");
1669     c->cd(2);
1670     gPad->SetGridy();
1671     gPad->SetGridx();
1672     hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
1673     hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
1674     hClusterLocMax ->Draw("colz");
1675     c->cd(3);
1676     gPad->SetGridy();
1677     gPad->SetGridx();
1678     hCluster1      ->SetAxisRange(minCol, maxCol,"X");
1679     hCluster1      ->SetAxisRange(minRow, maxRow,"Y");
1680     hCluster1      ->Draw("colz");
1681     c->cd(4);
1682     gPad->SetGridy();
1683     gPad->SetGridx();
1684     hCluster2      ->SetAxisRange(minCol, maxCol,"X");
1685     hCluster2      ->SetAxisRange(minRow, maxRow,"Y");
1686     hCluster2      ->Draw("colz");
1687     
1688     if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",
1689                                    eventNumber,cluster->E(),nMax,ncells1,ncells2));
1690     
1691     delete c;
1692     delete hClusterMap;
1693     delete hClusterLocMax;
1694     delete hCluster1;
1695     delete hCluster2;
1696   }
1697 }
1698