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