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