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