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