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