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