added run dependent corrections on T for calibration, removed old method
[u/mrichter/AliRoot.git] / ANALYSIS / TenderSupplies / AliEMCALTenderSupply.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 //                                                                           //
18 //  EMCAL tender, apply corrections to EMCAL clusters and do track matching. //
19 //                                                                           //
20 //  Author: Deepa Thomas (Utrecht University)                                // 
21 //  Later mods/rewrite: Jiri Kral (University of Jyvaskyla)                  //
22 //  S. Aiola, C. Loizides : Make it work for AODs                            //
23 //                                                                           //
24 ///////////////////////////////////////////////////////////////////////////////
25
26 #include <TROOT.h>
27 #include <TFile.h>
28 #include <TTree.h>
29 #include <TInterpreter.h>
30 #include <TObjArray.h>
31 #include <TClonesArray.h>
32 #include <TGeoGlobalMagField.h>
33 #include "AliEMCALAfterBurnerUF.h"
34 #include "AliEMCALClusterizer.h"
35 #include "AliEMCALClusterizerNxN.h"
36 #include "AliEMCALClusterizerv1.h"
37 #include "AliEMCALClusterizerv2.h"
38 #include "AliEMCALDigit.h"
39 #include "AliEMCALGeometry.h"
40 #include "AliEMCALRecParam.h"
41 #include "AliEMCALRecParam.h"
42 #include "AliEMCALRecPoint.h"
43 #include "AliEMCALRecoUtils.h"
44 #include "AliEMCALTenderSupply.h"
45 #include "AliESDCaloCluster.h"
46 #include "AliMagF.h"
47 #include "AliOADBContainer.h"
48 #include "AliAODEvent.h"
49 #include "AliAnalysisManager.h"
50 #include "AliESDEvent.h"
51 #include "AliLog.h"
52 #include "AliTender.h"
53
54 ClassImp(AliEMCALTenderSupply)
55
56 AliEMCALTenderSupply::AliEMCALTenderSupply() :
57 AliTenderSupply()
58 ,fTask(0)
59 ,fRun(0)
60 ,fEMCALGeo(0x0)
61 ,fEMCALGeoName("")
62 ,fEMCALRecoUtils(0)
63 ,fConfigName("")
64 ,fDebugLevel(0)
65 ,fNonLinearFunc(-1) 
66 ,fNonLinearThreshold(-1)
67 ,fReCalibCluster(kFALSE)
68 ,fUpdateCell(kFALSE)  
69 ,fCalibrateEnergy(kFALSE)
70 ,fCalibrateTime(kFALSE)
71 ,fDoNonLinearity(kFALSE)
72 ,fBadCellRemove(kFALSE)
73 ,fRejectExoticCells(kFALSE)
74 ,fRejectExoticClusters(kFALSE)
75 ,fClusterBadChannelCheck(kFALSE)
76 ,fRecalClusPos(kFALSE)
77 ,fFiducial(kFALSE) 
78 ,fNCellsFromEMCALBorder(-1)
79 ,fRecalDistToBadChannels(kFALSE)
80 ,fRecalShowerShape(kFALSE)
81 ,fInputTree(0)
82 ,fInputFile(0)
83 ,fFilepass(0) 
84 ,fMass(-1)
85 ,fStep(-1)
86 ,fCutEtaPhiSum(kTRUE)
87 ,fCutEtaPhiSeparate(kFALSE)
88 ,fRcut(-1)
89 ,fEtacut(-1)
90 ,fPhicut(-1)
91 ,fBasePath("")
92 ,fReClusterize(kFALSE)
93 ,fClusterizer(0)
94 ,fGeomMatrixSet(kFALSE)
95 ,fLoadGeomMatrices(kFALSE)
96 ,fRecParam(0x0)
97 ,fDoTrackMatch(kFALSE)
98 ,fDoUpdateOnly(kFALSE)
99 ,fUnfolder(0)
100 ,fDigitsArr(0)
101 ,fClusterArr(0)
102 ,fMisalignSurvey(kdefault)  
103 ,fExoticCellFraction(-1)
104 ,fExoticCellDiffTime(-1)
105 ,fExoticCellMinAmplitude(-1)
106 {
107   // Default constructor.
108
109   for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
110 }
111
112 //_____________________________________________________
113 AliEMCALTenderSupply::AliEMCALTenderSupply(const char *name, const AliTender *tender) :
114 AliTenderSupply(name,tender)
115 ,fTask(0)
116 ,fRun(0)
117 ,fEMCALGeo(0x0)
118 ,fEMCALGeoName("")
119 ,fEMCALRecoUtils(0)
120 ,fConfigName("")
121 ,fDebugLevel(0)
122 ,fNonLinearFunc(-1) 
123 ,fNonLinearThreshold(-1)        
124 ,fReCalibCluster(kFALSE)  
125 ,fUpdateCell(kFALSE)  
126 ,fCalibrateEnergy(kFALSE)
127 ,fCalibrateTime(kFALSE)
128 ,fDoNonLinearity(kFALSE)
129 ,fBadCellRemove(kFALSE)
130 ,fRejectExoticCells(kFALSE)
131 ,fRejectExoticClusters(kFALSE)
132 ,fClusterBadChannelCheck(kFALSE)
133 ,fRecalClusPos(kFALSE)
134 ,fFiducial(kFALSE) 
135 ,fNCellsFromEMCALBorder(-1)  
136 ,fRecalDistToBadChannels(kFALSE)  
137 ,fRecalShowerShape(kFALSE)
138 ,fInputTree(0)  
139 ,fInputFile(0)
140 ,fFilepass(0) 
141 ,fMass(-1)
142 ,fStep(-1)
143 ,fCutEtaPhiSum(kTRUE)
144 ,fCutEtaPhiSeparate(kFALSE)
145 ,fRcut(-1)  
146 ,fEtacut(-1)  
147 ,fPhicut(-1)  
148 ,fBasePath("")
149 ,fReClusterize(kFALSE)
150 ,fClusterizer(0)
151 ,fGeomMatrixSet(kFALSE)
152 ,fLoadGeomMatrices(kFALSE)
153 ,fRecParam(0x0)
154 ,fDoTrackMatch(kFALSE)
155 ,fDoUpdateOnly(kFALSE)
156 ,fUnfolder(0)
157 ,fDigitsArr(0)
158 ,fClusterArr(0)
159 ,fMisalignSurvey(kdefault)  
160 ,fExoticCellFraction(-1)
161 ,fExoticCellDiffTime(-1)
162 ,fExoticCellMinAmplitude(-1)
163 {
164   // Named constructor
165   
166   for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
167 }
168
169 //_____________________________________________________
170 AliEMCALTenderSupply::AliEMCALTenderSupply(const char *name, AliAnalysisTaskSE *task) :
171 AliTenderSupply(name)
172 ,fTask(task)
173 ,fRun(0)
174 ,fEMCALGeo(0x0)
175 ,fEMCALGeoName("")
176 ,fEMCALRecoUtils(0)
177 ,fConfigName("")
178 ,fDebugLevel(0)
179 ,fNonLinearFunc(-1) 
180 ,fNonLinearThreshold(-1)        
181 ,fReCalibCluster(kFALSE)  
182 ,fUpdateCell(kFALSE)  
183 ,fCalibrateEnergy(kFALSE)
184 ,fCalibrateTime(kFALSE)
185 ,fDoNonLinearity(kFALSE)
186 ,fBadCellRemove(kFALSE)
187 ,fRejectExoticCells(kFALSE)
188 ,fRejectExoticClusters(kFALSE)
189 ,fClusterBadChannelCheck(kFALSE)
190 ,fRecalClusPos(kFALSE)
191 ,fFiducial(kFALSE) 
192 ,fNCellsFromEMCALBorder(-1)  
193 ,fRecalDistToBadChannels(kFALSE)  
194 ,fRecalShowerShape(kFALSE)
195 ,fInputTree(0)  
196 ,fInputFile(0)
197 ,fFilepass(0) 
198 ,fMass(-1)
199 ,fStep(-1)
200 ,fCutEtaPhiSum(kTRUE)
201 ,fCutEtaPhiSeparate(kFALSE)
202 ,fRcut(-1)  
203 ,fEtacut(-1)  
204 ,fPhicut(-1)  
205 ,fBasePath("")
206 ,fReClusterize(kFALSE)
207 ,fClusterizer(0)
208 ,fGeomMatrixSet(kFALSE)
209 ,fLoadGeomMatrices(kFALSE)
210 ,fRecParam(0x0)
211 ,fDoTrackMatch(kFALSE)
212 ,fDoUpdateOnly(kFALSE)
213 ,fUnfolder(0)
214 ,fDigitsArr(0)
215 ,fClusterArr(0)
216 ,fMisalignSurvey(kdefault)  
217 ,fExoticCellFraction(-1)
218 ,fExoticCellDiffTime(-1)
219 ,fExoticCellMinAmplitude(-1)
220 {
221   // Named constructor.
222   
223   for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
224 }
225
226 //_____________________________________________________
227 AliEMCALTenderSupply::~AliEMCALTenderSupply()
228 {
229   //Destructor
230   
231   delete fEMCALRecoUtils;
232   delete fRecParam;
233   delete fUnfolder;
234   if (!fClusterizer) {
235     fDigitsArr->Clear("C");
236     delete fDigitsArr; 
237   } else {
238     delete fClusterizer;
239     fDigitsArr = 0;
240   }
241 }
242
243 //_____________________________________________________
244 void AliEMCALTenderSupply::SetDefaults()
245 {
246   // Set default settings.
247
248   SwitchOnBadCellRemove();
249   SwitchOnExoticCellRemove();
250   SwitchOnCalibrateEnergy();
251   SwitchOnCalibrateTime();
252   SwitchOnUpdateCell();
253   SwitchOnReclustering();
254   SwitchOnClusterBadChannelCheck();
255   SwitchOnClusterExoticChannelCheck();
256   SwitchOnTrackMatch();
257 }
258
259 //_____________________________________________________
260 Bool_t AliEMCALTenderSupply::RunChanged() const
261 {
262   // Get run number.
263
264   return (fTender && fTender->RunChanged()) || (fTask && fRun != fTask->InputEvent()->GetRunNumber()); 
265 }
266
267 //_____________________________________________________
268 void AliEMCALTenderSupply::Init()
269 {
270   // Initialise EMCAL tender.
271
272   if (fDebugLevel>0) 
273     AliWarning("Init EMCAL Tender supply"); 
274   
275   if (fConfigName.Length()>0 && gROOT->LoadMacro(fConfigName) >=0) {
276     AliDebug(1, Form("Loading settings from macro %s", fConfigName.Data()));
277     AliEMCALTenderSupply *tender = (AliEMCALTenderSupply*)gInterpreter->ProcessLine("ConfigEMCALTenderSupply()");
278     fDebugLevel             = tender->fDebugLevel;
279     fEMCALGeoName           = tender->fEMCALGeoName; 
280     fEMCALRecoUtils         = tender->fEMCALRecoUtils; 
281     fConfigName             = tender->fConfigName;
282     fNonLinearFunc          = tender->fNonLinearFunc;
283     fNonLinearThreshold     = tender->fNonLinearThreshold;
284     fReCalibCluster         = tender->fReCalibCluster;
285     fUpdateCell             = tender->fUpdateCell;
286     fRecalClusPos           = tender->fRecalClusPos;
287     fCalibrateEnergy        = tender->fCalibrateEnergy;
288     fCalibrateTime          = tender->fCalibrateTime;
289     fFiducial               = tender->fFiducial;
290     fNCellsFromEMCALBorder  = tender->fNCellsFromEMCALBorder;
291     fRecalDistToBadChannels = tender->fRecalDistToBadChannels;    
292     fRecalShowerShape       = tender->fRecalShowerShape;
293     fClusterBadChannelCheck = tender->fClusterBadChannelCheck;
294     fBadCellRemove          = tender->fBadCellRemove;
295     fRejectExoticCells      = tender->fRejectExoticCells;
296     fRejectExoticClusters   = tender->fRejectExoticClusters;
297     fMass                   = tender->fMass;
298     fStep                   = tender->fStep;
299     fCutEtaPhiSum           = tender->fCutEtaPhiSum;
300     fCutEtaPhiSeparate      = tender->fCutEtaPhiSeparate;
301     fRcut                   = tender->fRcut;
302     fEtacut                 = tender->fEtacut;
303     fPhicut                 = tender->fPhicut;
304     fReClusterize           = tender->fReClusterize;
305     fLoadGeomMatrices       = tender->fLoadGeomMatrices;
306     fRecParam               = tender->fRecParam;
307     fDoNonLinearity         = tender->fDoNonLinearity;
308     fDoTrackMatch           = tender->fDoTrackMatch;
309     fDoUpdateOnly           = tender->fDoUpdateOnly;
310     fMisalignSurvey         = tender->fMisalignSurvey;
311     fExoticCellFraction     = tender->fExoticCellFraction;
312     fExoticCellDiffTime     = tender->fExoticCellDiffTime;
313     fExoticCellMinAmplitude = tender->fExoticCellMinAmplitude;
314
315     for(Int_t i = 0; i < 12; i++) 
316       fEMCALMatrix[i] = tender->fEMCALMatrix[i] ;
317   }
318   
319   if (fDebugLevel>0){
320     AliInfo( "Emcal Tender settings: ======================================" ); 
321     AliInfo( "------------ Switches --------------------------" ); 
322     AliInfo( Form( "BadCellRemove : %d", fBadCellRemove )); 
323     AliInfo( Form( "ExoticCellRemove : %d", fRejectExoticCells )); 
324     AliInfo( Form( "CalibrateEnergy : %d", fCalibrateEnergy )); 
325     AliInfo( Form( "CalibrateTime : %d", fCalibrateTime )); 
326     AliInfo( Form( "UpdateCell : %d", fUpdateCell )); 
327     AliInfo( Form( "DoUpdateOnly : %d", fDoUpdateOnly )); 
328     AliInfo( Form( "Reclustering : %d", fReClusterize )); 
329     AliInfo( Form( "ClusterBadChannelCheck : %d", fClusterBadChannelCheck )); 
330     AliInfo( Form( "ClusterExoticChannelCheck : %d", fRejectExoticClusters )); 
331     AliInfo( Form( "CellFiducialRegion : %d", fFiducial )); 
332     AliInfo( Form( "ReCalibrateCluster : %d", fReCalibCluster )); 
333     AliInfo( Form( "RecalculateClusPos : %d", fRecalClusPos )); 
334     AliInfo( Form( "RecalShowerShape : %d", fRecalShowerShape )); 
335     AliInfo( Form( "NonLinearityCorrection : %d", fDoNonLinearity )); 
336     AliInfo( Form( "RecalDistBadChannel : %d", fRecalDistToBadChannels )); 
337     AliInfo( Form( "TrackMatch : %d", fDoTrackMatch )); 
338     AliInfo( "------------ Variables -------------------------" ); 
339     AliInfo( Form( "DebugLevel : %d", fDebugLevel )); 
340     AliInfo( Form( "BasePath : %s", fBasePath.Data() )); 
341     AliInfo( Form( "ConfigFileName : %s", fConfigName.Data() )); 
342     AliInfo( Form( "EMCALGeometryName : %s", fEMCALGeoName.Data() )); 
343     AliInfo( Form( "NonLinearityFunction : %d", fNonLinearFunc )); 
344     AliInfo( Form( "NonLinearityThreshold : %d", fNonLinearThreshold )); 
345     AliInfo( Form( "MisalignmentMatrixSurvey : %d", fMisalignSurvey )); 
346     AliInfo( Form( "NumberOfCellsFromEMCALBorder : %d", fNCellsFromEMCALBorder )); 
347     AliInfo( Form( "RCut : %f", fRcut )); 
348     AliInfo( Form( "Mass : %f", fMass )); 
349     AliInfo( Form( "Step : %f", fStep )); 
350     AliInfo( Form( "EtaCut : %f", fEtacut )); 
351     AliInfo( Form( "PhiCut : %f", fPhicut )); 
352     AliInfo( Form( "ExoticCellFraction : %f", fExoticCellFraction )); 
353     AliInfo( Form( "ExoticCellDiffTime : %f", fExoticCellDiffTime )); 
354     AliInfo( Form( "ExoticCellMinAmplitude : %f", fExoticCellMinAmplitude )); 
355     AliInfo( "=============================================================" ); 
356   }
357
358   // init reco utils
359   fEMCALRecoUtils  = new AliEMCALRecoUtils;
360
361   // init geometry if requested
362   if (fEMCALGeoName.Length()>0) 
363     fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
364
365   // digits array
366   fDigitsArr       = new TClonesArray("AliEMCALDigit",1000);
367
368   // initialising non-linearity parameters
369   if( fNonLinearThreshold != -1 )
370     fEMCALRecoUtils->SetNonLinearityThreshold(fNonLinearThreshold);
371   if( fNonLinearFunc != -1 )
372     fEMCALRecoUtils->SetNonLinearityFunction(fNonLinearFunc);
373
374   // missalignment function
375   fEMCALRecoUtils->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerGlobal);
376
377   // fiducial cut
378   // do not do the eta0 fiducial cut
379   if( fNCellsFromEMCALBorder != -1 )
380     fEMCALRecoUtils->SetNumberOfCellsFromEMCALBorder(fNCellsFromEMCALBorder);
381   fEMCALRecoUtils->SwitchOnNoFiducialBorderInEMCALEta0();
382     
383   // exotic cell rejection
384   if( fExoticCellFraction != -1 )
385     fEMCALRecoUtils->SetExoticCellFractionCut( fExoticCellFraction );
386   if( fExoticCellDiffTime != -1 )
387     fEMCALRecoUtils->SetExoticCellDiffTimeCut( fExoticCellDiffTime );
388   if( fExoticCellMinAmplitude != -1 )
389     fEMCALRecoUtils->SetExoticCellMinAmplitudeCut( fExoticCellMinAmplitude );
390
391   // setting track matching parameters ... mass, step size and residual cut
392   if( fMass != -1 )
393     fEMCALRecoUtils->SetMass(fMass);
394   if( fStep != -1 )
395     fEMCALRecoUtils->SetStep(fStep);
396   
397   // spatial cut based on separate eta/phi or common processing
398   if(fCutEtaPhiSum){ 
399     fEMCALRecoUtils->SwitchOnCutEtaPhiSum(); 
400     if( fRcut != -1 )
401       fEMCALRecoUtils->SetCutR(fRcut);
402   } else if (fCutEtaPhiSeparate) {
403     fEMCALRecoUtils->SwitchOnCutEtaPhiSeparate();
404     if( fEtacut != -1 )
405       fEMCALRecoUtils->SetCutEta(fEtacut);
406     if( fPhicut != -1 )
407       fEMCALRecoUtils->SetCutPhi(fPhicut);
408   }
409 }
410
411 //_____________________________________________________
412 AliVEvent* AliEMCALTenderSupply::GetEvent()
413 {
414   // Return the event pointer.
415   
416   if (fTender) {
417     return fTender->GetEvent();
418   }
419   else if (fTask) {
420     return fTask->InputEvent();
421   }
422   
423   return 0;
424 }
425
426 //_____________________________________________________
427 void AliEMCALTenderSupply::ProcessEvent()
428 {
429   // Event loop.
430   
431   AliVEvent *event = GetEvent();
432   
433   if (!event) {
434     AliError("Event ptr = 0, returning");
435     return;
436   }
437   
438  // Initialising parameters once per run number
439   if (RunChanged()) { 
440
441     AliWarning( "Run changed, initializing parameters" );
442     fRun = event->GetRunNumber();
443
444     // init geometry if not already done
445     if (fEMCALGeoName.Length()==0) {
446       fEMCALGeoName = "EMCAL_FIRSTYEARV1";
447       if (fRun>139517) {
448         fEMCALGeoName = "EMCAL_COMPLETEV1";
449       } 
450       if (fRun>170593) {
451         fEMCALGeoName = "EMCAL_COMPLETE12SMV1";
452       }
453       fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
454       if (!fEMCALGeo) {
455         AliFatal(Form("Can not create geometry: %s", fEMCALGeoName.Data()));
456         return;
457       }
458     } 
459
460     // get pass
461     GetPass();
462
463     // define what recalib parameters are needed for various switches
464     // this is based on implementation in AliEMCALRecoUtils
465     Bool_t needRecoParam   = fReClusterize;
466     Bool_t needBadChannels = fBadCellRemove   | fClusterBadChannelCheck | fRecalDistToBadChannels | fReClusterize;
467     Bool_t needRecalib     = fCalibrateEnergy | fReClusterize;
468     Bool_t needTimecalib   = fCalibrateTime   | fReClusterize;
469     Bool_t needMisalign    = fRecalClusPos    | fReClusterize;
470     Bool_t needClusterizer = fReClusterize;
471
472     // initiate reco params with some defaults 
473     // will not overwrite, if those have been provided by user
474     if( needRecoParam ){
475       Int_t initRC = InitRecParam();
476       
477       if( initRC == 0 )
478         AliInfo("Defaults reco params loaded.");
479       if( initRC > 1 )
480         AliWarning("User defined reco params.");
481     }
482
483     // init bad channels
484     if( needBadChannels ){
485       Int_t fInitBC = InitBadChannels();
486       if (fInitBC==0)
487         AliError("InitBadChannels returned false, returning");
488       if (fInitBC==1)
489         AliWarning("InitBadChannels OK");
490       if (fInitBC>1)
491         AliWarning(Form("No external hot channel set: %d - %s", event->GetRunNumber(), fFilepass.Data()));
492     }
493
494     // init recalibration factors
495     if( needRecalib ) 
496     { 
497       Int_t fInitRecalib = InitRecalib();
498       if (fInitRecalib==0)
499         AliError("InitRecalib returned false, returning");
500       if (fInitRecalib==1)
501         AliWarning("InitRecalib OK");
502       if (fInitRecalib>1)
503         AliWarning(Form("No recalibration available: %d - %s", event->GetRunNumber(), fFilepass.Data()));
504       
505       Int_t fInitRunDepRecalib = InitRunDepRecalib();
506       if (fInitRunDepRecalib==0)
507         AliError("InitrunDepRecalib returned false, returning");
508       if (fInitRunDepRecalib==1)
509         AliWarning("InitRecalib OK");
510       if (fInitRunDepRecalib>1)
511         AliWarning(Form("No Temperature recalibration available: %d - %s", event->GetRunNumber(), fFilepass.Data()));
512       
513       fReCalibCluster = kFALSE;
514       
515     }
516     
517     
518     
519     // init time calibration
520     if( needTimecalib ){
521       Int_t initTC = InitTimeCalibration();
522       if ( !initTC ) 
523         AliError("InitTimeCalibration returned false, returning");
524       if (initTC==1)
525         AliWarning("InitTimeCalib OK");
526       if( initTC > 1 )
527         AliWarning(Form("No external time calibration set: %d - %s", event->GetRunNumber(), fFilepass.Data()));
528     }
529
530     // init misalignment matrix
531     if( needMisalign ) { 
532       if (!InitMisalignMatrix())
533         AliError("InitMisalignmentMatrix returned false, returning");
534       else
535         AliWarning("InitMisalignMatrix OK");
536     }
537     
538     // init clusterizer
539     if( needClusterizer ) {
540       if (!InitClusterization()) 
541         AliError("InitClusterization returned false, returning");
542       else
543         AliWarning("InitClusterization OK");
544     }
545     
546     if(fDebugLevel>1) 
547       fEMCALRecoUtils->Print("");
548   }
549   
550   // disable implied switches -------------------------------------------------
551   // AliEMCALRecoUtils or clusterizer functions alredy include some
552   // recalibration so based on those implied callibration te switches
553   // here are disabled to avoid duplication
554     
555   // clusterizer does cluster energy recalibration, position recomputation
556   // and shower shape
557   if( fReClusterize ){
558     fReCalibCluster   = kFALSE;
559     fRecalClusPos     = kFALSE;
560     fRecalShowerShape = kFALSE;
561   }
562   
563   // CONFIGURE THE RECO UTILS -------------------------------------------------
564   // configure the reco utils
565   // this option does energy recalibration
566   if( fCalibrateEnergy )
567     fEMCALRecoUtils->SwitchOnRecalibration();
568   else
569     fEMCALRecoUtils->SwitchOffRecalibration();
570   
571   // allows time calibration
572   if( fCalibrateTime )
573     fEMCALRecoUtils->SwitchOnTimeRecalibration();
574   else
575     fEMCALRecoUtils->SwitchOffTimeRecalibration();
576
577   // allows to zero bad cells
578   if( fBadCellRemove )
579     fEMCALRecoUtils->SwitchOnBadChannelsRemoval();
580   else
581     fEMCALRecoUtils->SwitchOffBadChannelsRemoval();
582   
583   // distance to bad channel recalibration
584   if( fRecalDistToBadChannels )
585     fEMCALRecoUtils->SwitchOnDistToBadChannelRecalculation();
586   else
587     fEMCALRecoUtils->SwitchOffDistToBadChannelRecalculation();
588
589   // exclude exotic cells
590   if( fRejectExoticCells )
591     fEMCALRecoUtils->SwitchOnRejectExoticCell();
592   else
593     fEMCALRecoUtils->SwitchOffRejectExoticCell();
594   
595   // exclude clusters with exotic cells
596   if( fRejectExoticClusters )
597     fEMCALRecoUtils->SwitchOnRejectExoticCluster();
598   else
599     fEMCALRecoUtils->SwitchOffRejectExoticCluster();
600
601   // TODO: not implemented switches
602   // SwitchOnClusterEnergySmearing
603   // SwitchOnRunDepCorrection
604
605   // START PROCESSING ---------------------------------------------------------
606   // Test if cells present
607   AliVCaloCells *cells= event->GetEMCALCells();
608   if (cells->GetNumberOfCells()<=0) 
609   {
610     if(fDebugLevel>1) 
611       AliWarning(Form("Number of EMCAL cells = %d, returning", cells->GetNumberOfCells()));
612     return;
613   }
614   
615   if (fDebugLevel>2)
616     AliInfo(Form("Re-calibrate cluster %d\n",fReCalibCluster));
617
618   // mark the cells not recalibrated in case of selected
619   // time, energy recalibration or bad channel removal
620   if( fCalibrateEnergy || fCalibrateTime || fBadCellRemove )
621     fEMCALRecoUtils->ResetCellsCalibrated();
622   
623  // CELL RECALIBRATION -------------------------------------------------------
624   // cell objects will be updated
625   // the cell calibrations are also processed locally any time those are needed
626   // in case that the cell objects are not to be updated here for later use
627   if( fUpdateCell )
628   {
629     // do the update
630     // includes exotic cell check in the UpdateCells function - is not provided
631     // by the reco utils
632     UpdateCells();
633
634     // switch off recalibrations so those are not done multiple times
635     // this is just for safety, the recalibrated flag of cell object
636     // should not allow for farther processing anyways
637     fEMCALRecoUtils->SwitchOffRecalibration();
638     fEMCALRecoUtils->SwitchOffTimeRecalibration();  
639   
640     if (fDoUpdateOnly)
641       return;
642   }
643
644   // RECLUSTERIZATION ---------------------------------------------------------
645   if (fReClusterize)
646   {
647     FillDigitsArray();
648     Clusterize();
649     UpdateClusters();
650   }
651
652   // Store good clusters
653   TClonesArray *clusArr = dynamic_cast<TClonesArray*>(event->FindListObject("caloClusters"));
654   if (!clusArr) 
655     clusArr = dynamic_cast<TClonesArray*>(event->FindListObject("CaloClusters"));
656   if (!clusArr) {
657     AliWarning(Form("No cluster array, number of cells in event = %d, returning", cells->GetNumberOfCells()));
658     return;
659   }
660
661   // PROCESS SINGLE CLUSTER RECALIBRATION -------------------------------------
662   // now go through clusters one by one and process separate correction
663   // as those were defined or not
664   Int_t nclusters = clusArr->GetEntriesFast();
665   for (Int_t icluster=0; icluster < nclusters; ++icluster) 
666   { 
667     AliVCluster *clust = static_cast<AliVCluster*>(clusArr->At(icluster));
668     if (!clust) 
669       continue;
670     if  (!clust->IsEMCAL()) 
671       continue;
672
673     // REMOVE CLUSTERS WITH BAD CELLS -----------------------------
674     if( fClusterBadChannelCheck )
675     {
676       // careful, the the ClusterContainsBadChannel is dependent on
677       // SwitchOnBadChannelsRemoval, switching it ON automatically
678       // and returning to original value after processing
679       Bool_t badRemoval = fEMCALRecoUtils->IsBadChannelsRemovalSwitchedOn();
680       fEMCALRecoUtils->SwitchOnBadChannelsRemoval();
681       
682       Bool_t badResult = fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo, clust->GetCellsAbsId(), clust->GetNCells());
683
684       // switch the bad channels removal back
685       if( ! badRemoval )
686         fEMCALRecoUtils->SwitchOffBadChannelsRemoval();
687       
688       if( badResult )
689       {
690         delete clusArr->RemoveAt(icluster);
691         continue; //TODO is it really needed to remove it? Or should we flag it?
692       }
693     }
694     
695     // REMOVE EXOTIC CLUSTERS -------------------------------------
696     // does process local cell recalibration energy and time without replacing
697     // the global cell values, in case of no cell recalib done yet
698     if( fRejectExoticClusters )
699     {
700       // careful, the IsExoticCluster is dependent on
701       // SwitchOnRejectExoticCell, switching it ON automatically
702       // and returning to original value after processing
703       Bool_t exRemoval = fEMCALRecoUtils->IsRejectExoticCell();
704       fEMCALRecoUtils->SwitchOnRejectExoticCell();
705
706       // get bunch crossing
707       Int_t bunchCrossNo = event->GetBunchCrossNumber();
708
709       Bool_t exResult = fEMCALRecoUtils->IsExoticCluster(clust, cells, bunchCrossNo );
710
711       // switch the exotic channels removal back
712       if( ! exRemoval )
713         fEMCALRecoUtils->SwitchOffRejectExoticCell();
714       
715       if( exResult )
716       {
717         delete clusArr->RemoveAt(icluster);
718         continue; //TODO is it really needed to remove it? Or should we flag it?
719       }
720     }
721     
722     // FIDUCIAL CUT -----------------------------------------------
723     if (fFiducial)
724     {
725       // depends on SetNumberOfCellsFromEMCALBorder
726       // SwitchOnNoFiducialBorderInEMCALEta0
727       if (!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, clust, cells)){
728         delete clusArr->RemoveAt(icluster);
729         continue; //TODO it would be nice to store the distance
730       }
731     }
732     
733     // CLUSTER ENERGY ---------------------------------------------
734     // does process local cell recalibration energy and time without replacing
735     // the global cell values, in case of no cell recalib done yet
736     if( fReCalibCluster ) 
737       fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, clust, cells);
738
739     // CLUSTER POSITION -------------------------------------------
740     // does process local cell energy recalibration, if enabled and cells
741     // not calibrated yet
742     if( fRecalClusPos ) 
743       fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, cells, clust);
744     
745     // SHOWER SHAPE -----------------------------------------------
746     if( fRecalShowerShape )
747       fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, cells, clust);  
748
749     // NONLINEARITY -----------------------------------------------
750     if( fDoNonLinearity )
751     {
752       Float_t correctedEnergy = fEMCALRecoUtils->CorrectClusterEnergyLinearity(clust);
753       clust->SetE(correctedEnergy);
754     }
755
756     // DISTANCE TO BAD CHANNELS -----------------------------------
757     if( fRecalDistToBadChannels )
758       fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, cells, clust);  
759   }
760
761   clusArr->Compress();
762
763   if (!fDoTrackMatch)
764     return;
765
766   // TRACK MATCHING -----------------------------------------------------------
767   if (!TGeoGlobalMagField::Instance()->GetField()) {
768     AliESDEvent *esd = dynamic_cast<AliESDEvent*>(event);
769     if (esd)
770       esd->InitMagneticField();
771     else {
772       AliAODEvent *aod = dynamic_cast<AliAODEvent*>(event);
773       Double_t curSol = 30000*aod->GetMagneticField()/5.00668;
774       Double_t curDip = 6000 *aod->GetMuonMagFieldScale();
775       AliMagF *field  = AliMagF::CreateFieldMap(curSol,curDip);
776       TGeoGlobalMagField::Instance()->SetField(field);
777     }
778   }
779
780   fEMCALRecoUtils->FindMatches(event,0x0,fEMCALGeo);
781   fEMCALRecoUtils->SetClusterMatchedToTrack(event);
782   fEMCALRecoUtils->SetTracksMatchedToCluster(event);
783 }
784
785 //_____________________________________________________
786 Bool_t AliEMCALTenderSupply::InitMisalignMatrix()
787 {
788   // Initialising misalignment matrices
789
790   AliVEvent *event = GetEvent();
791   
792   if (!event) 
793     return kFALSE;
794   
795   if (fGeomMatrixSet) 
796   {
797     AliInfo("Misalignment matrix already set");  
798     return kTRUE;
799   }
800   
801   if (fDebugLevel>0) 
802     AliInfo("Initialising misalignment matrix");  
803   
804   if (fLoadGeomMatrices) {
805     for(Int_t mod=0; mod < fEMCALGeo->GetNumberOfSuperModules(); ++mod)
806     {
807       if (fEMCALMatrix[mod]){
808         if(fDebugLevel > 2) 
809           fEMCALMatrix[mod]->Print();
810         fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod);  
811       }
812     }
813     fGeomMatrixSet = kTRUE;
814     return kTRUE;
815   }
816   
817   Int_t runGM = event->GetRunNumber();
818   TObjArray *mobj = 0;
819
820   if(fMisalignSurvey == kdefault)
821   { //take default alignment corresponding to run no
822     AliOADBContainer emcalgeoCont(Form("emcal"));
823     emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
824     mobj=(TObjArray*)emcalgeoCont.GetObject(runGM,"EmcalMatrices");
825   }
826   
827   if(fMisalignSurvey == kSurveybyS)
828   { //take alignment at sector level
829     if (runGM <= 140000) { //2010 data
830       AliOADBContainer emcalgeoCont(Form("emcal2010"));
831       emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
832       mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey10");
833     } 
834     else if (runGM>140000)
835     { // 2011 LHC11a pass1 data
836       AliOADBContainer emcalgeoCont(Form("emcal2011"));
837       emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
838       mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey11byS");      
839     }
840   }
841
842   if(fMisalignSurvey == kSurveybyM)
843   { //take alignment at module level
844     if (runGM <= 140000) { //2010 data
845       AliOADBContainer emcalgeoCont(Form("emcal2010"));
846       emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
847       mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey10");
848     } 
849     else if (runGM>140000) 
850     { // 2011 LHC11a pass1 data
851       AliOADBContainer emcalgeoCont(Form("emcal2011"));
852       emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
853       mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey11byM");      
854     }
855   }
856
857   if(!mobj)
858   {
859     AliFatal("Geometry matrix array not found");
860     return kFALSE;
861   }
862   
863   for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
864   {
865     fEMCALMatrix[mod] = (TGeoHMatrix*) mobj->At(mod);
866     fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod); 
867     fEMCALMatrix[mod]->Print();
868   }
869   
870   return kTRUE;
871 }
872
873 //_____________________________________________________
874 Int_t AliEMCALTenderSupply::InitBadChannels()
875 {
876   // Initialising bad channel maps
877
878   AliVEvent *event = GetEvent();
879
880   if (!event) 
881     return 0;
882   
883   if (fDebugLevel>0) 
884     AliInfo("Initialising Bad channel map");
885   
886   // init default maps first
887   if( !fEMCALRecoUtils->GetEMCALBadChannelStatusMapArray() )
888     fEMCALRecoUtils->InitEMCALBadChannelStatusMap() ;
889   
890   Int_t runBC = event->GetRunNumber();
891   
892   AliOADBContainer *contBC = new AliOADBContainer("");
893   if (fBasePath!="")
894   { //if fBasePath specified in the ->SetBasePath()
895     if (fDebugLevel>0) AliInfo(Form("Loading Bad Channels OADB from given path %s",fBasePath.Data()));
896     
897     TFile *fbad=new TFile(Form("%s/EMCALBadChannels.root",fBasePath.Data()),"read");
898     if (!fbad || fbad->IsZombie())
899     {
900       AliFatal(Form("EMCALBadChannels.root was not found in the path provided: %s",fBasePath.Data()));
901       return 0;
902     }  
903     
904     if (fbad) delete fbad;
905     
906     contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fBasePath.Data()),"AliEMCALBadChannels");    
907   } 
908   else 
909   { // Else choose the one in the $ALICE_ROOT directory
910     if (fDebugLevel>0) AliInfo("Loading Bad Channels OADB from $ALICE_ROOT/OADB/EMCAL");
911     
912     TFile *fbad=new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root","read");
913     if (!fbad || fbad->IsZombie())
914     {
915       AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root was not found");
916       return 0;
917     }  
918       
919     if (fbad) delete fbad;
920     
921     contBC->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root","AliEMCALBadChannels"); 
922   }
923   
924   TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runBC);
925   if (!arrayBC)
926   {
927     AliError(Form("No external hot channel set for run number: %d", runBC));
928     return 2; 
929   }
930
931   Int_t sms = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules();
932   for (Int_t i=0; i<sms; ++i) 
933   {
934     TH2I *h = fEMCALRecoUtils->GetEMCALChannelStatusMap(i);
935     if (h)
936       delete h;
937     h=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
938
939     if (!h) 
940     {
941       AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
942       continue;
943     }
944     h->SetDirectory(0);
945     fEMCALRecoUtils->SetEMCALChannelStatusMap(i,h);
946   }
947   return 1;  
948 }
949
950 //_____________________________________________________
951 Int_t AliEMCALTenderSupply::InitRecalib()
952 {
953   // Initialising recalibration factors.
954   
955   AliVEvent *event = GetEvent();
956
957   if (!event) 
958     return 0;
959   
960   if (fDebugLevel>0) 
961     AliInfo("Initialising recalibration factors");
962   
963   // init default maps first
964   if( !fEMCALRecoUtils->GetEMCALRecalibrationFactorsArray() )
965     fEMCALRecoUtils->InitEMCALRecalibrationFactors() ;
966
967   Int_t runRC = event->GetRunNumber();
968       
969   AliOADBContainer *contRF=new AliOADBContainer("");
970   if (fBasePath!="") 
971   { //if fBasePath specified in the ->SetBasePath()
972     if (fDebugLevel>0)  AliInfo(Form("Loading Recalib OADB from given path %s",fBasePath.Data()));
973     
974     TFile *fRecalib= new TFile(Form("%s/EMCALRecalib.root",fBasePath.Data()),"read");
975     if (!fRecalib || fRecalib->IsZombie()) 
976     {
977       AliFatal(Form("EMCALRecalib.root not found in %s",fBasePath.Data()));
978       return 0;
979     }
980     
981     if (fRecalib) delete fRecalib;
982     
983     contRF->InitFromFile(Form("%s/EMCALRecalib.root",fBasePath.Data()),"AliEMCALRecalib");
984   }
985   else
986   { // Else choose the one in the $ALICE_ROOT directory
987     if (fDebugLevel>0)  AliInfo("Loading Recalib OADB from $ALICE_ROOT/OADB/EMCAL");
988     
989     TFile *fRecalib= new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root","read");
990     if (!fRecalib || fRecalib->IsZombie()) 
991     {
992       AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root was not found");
993       return 0;
994     }
995     
996     if (fRecalib) delete fRecalib;
997       
998     contRF->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root","AliEMCALRecalib");     
999   }
1000
1001   TObjArray *recal=(TObjArray*)contRF->GetObject(runRC);
1002   if (!recal)
1003   {
1004     AliError(Form("No Objects for run: %d",runRC));
1005     return 2;
1006   } 
1007
1008   TObjArray *recalpass=(TObjArray*)recal->FindObject(fFilepass);
1009   if (!recalpass)
1010   {
1011     AliError(Form("No Objects for run: %d - %s",runRC,fFilepass.Data()));
1012     return 2;
1013   }
1014
1015   TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
1016   if (!recalib)
1017   {
1018     AliError(Form("No Recalib histos found for  %d - %s",runRC,fFilepass.Data())); 
1019     return 2;
1020   }
1021
1022   if (fDebugLevel>0) recalib->Print();
1023
1024   Int_t sms = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules();
1025   for (Int_t i=0; i<sms; ++i) 
1026   {
1027     TH2F *h = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactors(i);
1028     if (h)
1029       delete h;
1030     h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
1031     if (!h) 
1032     {
1033       AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
1034       continue;
1035     }
1036     h->SetDirectory(0);
1037     fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
1038   }
1039   return 1;
1040 }
1041
1042 //_____________________________________________________
1043 Int_t AliEMCALTenderSupply::InitRunDepRecalib()
1044 {
1045   // Initialising recalibration factors.
1046   
1047   AliVEvent *event = GetEvent();
1048   
1049   if (!event) 
1050     return 0;
1051   
1052   if (fDebugLevel>0) 
1053     AliInfo("Initialising recalibration factors");
1054   
1055   // init default maps first
1056   if( !fEMCALRecoUtils->GetEMCALRecalibrationFactorsArray() )
1057     fEMCALRecoUtils->InitEMCALRecalibrationFactors() ;
1058   
1059   Int_t runRC = event->GetRunNumber();
1060   
1061   AliOADBContainer *contRF=new AliOADBContainer("");
1062   if (fBasePath!="") 
1063   { //if fBasePath specified in the ->SetBasePath()
1064     if (fDebugLevel>0)  AliInfo(Form("Loading Recalib OADB from given path %s",fBasePath.Data()));
1065     
1066     TFile *fRunDepRecalib= new TFile(Form("%s/EMCALTemperatureCorrCalib.root",fBasePath.Data()),"read");
1067     if (!fRunDepRecalib || fRunDepRecalib->IsZombie()) 
1068     {
1069       AliFatal(Form("EMCALTemperatureCorrCalib.root not found in %s",fBasePath.Data()));
1070       return 0;
1071     }
1072     
1073     if (fRunDepRecalib) delete fRunDepRecalib;
1074     
1075     contRF->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fBasePath.Data()),"AliEMCALRunDepTempCalibCorrections");
1076   }
1077   else
1078   { // Else choose the one in the $ALICE_ROOT directory
1079     if (fDebugLevel>0)  AliInfo("Loading Recalib OADB from $ALICE_ROOT/OADB/EMCAL");
1080     
1081     TFile *fRunDepRecalib= new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALTemperatureCorrCalib.root","read");
1082     if (!fRunDepRecalib || fRunDepRecalib->IsZombie()) 
1083     {
1084       AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALTemperatureCorrCalib.root was not found");
1085       return 0;
1086     }
1087     
1088     if (fRunDepRecalib) delete fRunDepRecalib;
1089     
1090     contRF->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALTemperatureCorrCalib.root","AliEMCALRunDepTempCalibCorrections");     
1091   }
1092   
1093   TH1S *rundeprecal=(TH1S*)contRF->GetObject(runRC);
1094   if (!rundeprecal)
1095   {
1096     AliError(Form("No Objects for run: %d",runRC));
1097     return 2;
1098   } 
1099   
1100   if (fDebugLevel>0) rundeprecal->Print();
1101   
1102   Int_t nSM = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules();
1103   
1104   for (Int_t ism=0; ism<nSM; ++ism) 
1105   {        
1106     for (Int_t icol=0; icol<48; ++icol) 
1107     {        
1108       for (Int_t irow=0; irow<24; ++irow) 
1109       {
1110         Float_t factor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
1111         
1112         Int_t absID = fEMCALGeo->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
1113         factor *= rundeprecal->GetBinContent(absID) / 10000. ; // correction dependent on T
1114         //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID, 
1115         //      GetEMCALChannelRecalibrationFactor(ism,icol,irow) , rundeprecal->GetBinContent(absID) / 10000., factor);
1116         fEMCALRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
1117       } // columns
1118     } // rows 
1119   } // SM loop
1120   
1121   return 1;
1122 }
1123
1124
1125 //_____________________________________________________
1126 Int_t AliEMCALTenderSupply::InitTimeCalibration()
1127 {
1128   // Initialising bad channel maps
1129   AliVEvent *event = GetEvent();
1130
1131   if (!event) 
1132     return 0;
1133   
1134   if (fDebugLevel>0) 
1135     AliInfo("Initialising time calibration map");
1136   
1137   // init default maps first
1138   if ( !fEMCALRecoUtils->GetEMCALTimeRecalibrationFactorsArray() )
1139     fEMCALRecoUtils->InitEMCALTimeRecalibrationFactors() ;
1140
1141   Int_t runBC = event->GetRunNumber();
1142   
1143   AliOADBContainer *contBC = new AliOADBContainer("");
1144   if (fBasePath!="")
1145   { //if fBasePath specified in the ->SetBasePath()
1146     if (fDebugLevel>0) AliInfo(Form("Loading time calibration OADB from given path %s",fBasePath.Data()));
1147     
1148     TFile *fbad=new TFile(Form("%s/EMCALTimeCalib.root",fBasePath.Data()),"read");
1149     if (!fbad || fbad->IsZombie())
1150     {
1151       AliFatal(Form("EMCALTimeCalib.root was not found in the path provided: %s",fBasePath.Data()));
1152       return 0;
1153     }  
1154     
1155     if (fbad) delete fbad;
1156     
1157     contBC->InitFromFile(Form("%s/EMCALTimeCalib.root",fBasePath.Data()),"AliEMCALTimeCalib");    
1158   } 
1159   else 
1160   { // Else choose the one in the $ALICE_ROOT directory
1161     if (fDebugLevel>0) AliInfo("Loading time calibration OADB from $ALICE_ROOT/OADB/EMCAL");
1162     
1163     TFile *fbad=new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALTimeCalib.root","read");
1164     if (!fbad || fbad->IsZombie())
1165     {
1166       AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALTimeCalib.root was not found");
1167       return 0;
1168     }  
1169       
1170     if (fbad) delete fbad;
1171     
1172     contBC->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALTimeCalib.root","AliEMCALTimeCalib"); 
1173   }
1174   
1175   TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runBC);
1176   if (!arrayBC)
1177   {
1178     AliError(Form("No external time calibration set for run number: %d", runBC));
1179     return 2; 
1180   }
1181   
1182   // Here, it looks for a specific pass
1183   TObjArray *arrayBCpass=(TObjArray*)arrayBC->FindObject(fFilepass); 
1184   if (!arrayBCpass)
1185   {
1186     AliError(Form("No external time calibration set for: %d -%s", runBC,fFilepass.Data()));
1187     return 2; 
1188   }
1189
1190   if (fDebugLevel>0) arrayBCpass->Print();
1191
1192   for( Int_t i = 0; i < 4; i++ )
1193   {
1194     TH1F *h = fEMCALRecoUtils->GetEMCALChannelTimeRecalibrationFactors( i );
1195     if( h )
1196       delete h;
1197     
1198     h = (TH1F*)arrayBCpass->FindObject(Form("hAllTimeAvBC%d",i));
1199     
1200     if (!h)
1201     {
1202       AliError(Form("Can not get hAllTimeAvBC%d",i));
1203       continue;
1204     }
1205     h->SetDirectory(0);
1206     fEMCALRecoUtils->SetEMCALChannelTimeRecalibrationFactors(i,h);
1207   }
1208   return 1;  
1209 }
1210
1211 //_____________________________________________________
1212 void AliEMCALTenderSupply::UpdateCells()
1213 {
1214   //Remove bad cells from the cell list
1215   //Recalibrate energy and time cells 
1216   //This is required for later reclusterization
1217
1218   AliVEvent *event = GetEvent();
1219
1220   if(!event) return ;
1221   
1222   AliVCaloCells *cells = event->GetEMCALCells();
1223   Int_t bunchCrossNo = event->GetBunchCrossNumber();
1224
1225   fEMCALRecoUtils->RecalibrateCells(cells, bunchCrossNo); 
1226   
1227   // remove exotic cells - loop through cells and zero the exotic ones
1228   // just like with bad cell rejection in reco utils (inside RecalibrateCells)
1229   if( fRejectExoticCells )
1230   {
1231     Short_t  absId  =-1;
1232     Double_t ecell = 0;
1233     Double_t tcell = 0;
1234     Double_t efrac = 0;
1235     Short_t  mclabel = -1;
1236     Bool_t   isExot = kFALSE;
1237   
1238     // loop through cells
1239     Int_t nEMcell  = cells->GetNumberOfCells() ;  
1240     for (Int_t iCell = 0; iCell < nEMcell; iCell++) 
1241     { 
1242       cells->GetCell( iCell, absId, ecell, tcell, mclabel, efrac );
1243     
1244       isExot = fEMCALRecoUtils->IsExoticCell( absId, cells, bunchCrossNo ); 
1245       // zero if exotic
1246       if( isExot )
1247         cells->SetCell( iCell, absId, 0.0, -1.0, mclabel, efrac );
1248     } // cell loop
1249   } // reject exotic cells
1250
1251   cells->Sort();
1252 }
1253
1254 //_____________________________________________________
1255 TString AliEMCALTenderSupply::GetBeamType()
1256 {
1257   // Get beam type : pp-AA-pA
1258   // ESDs have it directly, AODs get it from hardcoded run number ranges
1259   
1260   AliVEvent *event = GetEvent();
1261
1262   if (!event) { 
1263     AliError("Couldn't retrieve event!");
1264     return "";
1265   }
1266
1267   TString beamType;
1268
1269   AliESDEvent *esd = dynamic_cast<AliESDEvent*>(event);
1270   if (esd) {
1271     const AliESDRun *run = esd->GetESDRun();
1272     beamType = run->GetBeamType();
1273   }
1274   else
1275   {
1276     Int_t runNumber = event->GetRunNumber();
1277     if ((runNumber >= 136851 && runNumber <= 139517)  // LHC10h
1278   || (runNumber >= 166529 && runNumber <= 170593))  // LHC11h
1279     {
1280       beamType = "A-A";
1281     }
1282     else 
1283     {
1284       beamType = "p-p";
1285     }
1286   }
1287
1288   return beamType;    
1289 }
1290
1291 //_____________________________________________________
1292 Int_t AliEMCALTenderSupply::InitRecParam()
1293 {
1294   // Init reco params if not yet exist (probably shipped by the user already)
1295
1296   if( fRecParam != 0 )
1297     return 2;
1298
1299   TString beamType = GetBeamType();
1300
1301   // set some default reco params
1302   fRecParam = new AliEMCALRecParam();
1303   fRecParam->SetClusteringThreshold(0.100);
1304   fRecParam->SetMinECut(0.050);
1305   fRecParam->SetTimeCut(250);
1306   fRecParam->SetTimeMin(425);
1307   fRecParam->SetTimeMax(825);
1308   if ( beamType == "A-A") {
1309     fRecParam->SetClusterizerFlag(AliEMCALRecParam::kClusterizerv2);
1310   } 
1311   else {
1312     fRecParam->SetClusterizerFlag(AliEMCALRecParam::kClusterizerv1);
1313   }
1314
1315   return 0;
1316 }
1317
1318 //_____________________________________________________
1319 Bool_t AliEMCALTenderSupply::InitClusterization()
1320 {
1321   // Initialing clusterization/unfolding algorithm and set all the needed parameters.
1322   
1323   AliVEvent *event = GetEvent();
1324
1325   if (!event) 
1326     return kFALSE;
1327   
1328   if (fDebugLevel>0) 
1329     AliInfo(Form("Initialising reclustering parameters: Clusterizer type: %d",fRecParam->GetClusterizerFlag()));
1330   
1331   //---setup clusterizer
1332   delete fClusterizer;
1333   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1334     fClusterizer = new AliEMCALClusterizerv1 (fEMCALGeo);
1335   else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2) 
1336     fClusterizer = new AliEMCALClusterizerv2(fEMCALGeo);
1337   else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) 
1338   {
1339     AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fEMCALGeo);
1340     clusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
1341     clusterizer->SetNColDiff(fRecParam->GetNColDiff());
1342     fClusterizer = clusterizer;
1343   } 
1344   else 
1345   {
1346     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1347     return kFALSE;
1348   }
1349   
1350   // Set the clustering parameters
1351   fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
1352   fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
1353   fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );    
1354   fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
1355   fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
1356   fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
1357   fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
1358   fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
1359   fClusterizer->SetInputCalibrated       ( kTRUE                               );
1360   fClusterizer->SetJustClusters          ( kTRUE                               );  
1361   
1362   // In case of unfolding after clusterization is requested, set the corresponding parameters
1363   if (fRecParam->GetUnfold()) 
1364   {
1365     for (Int_t i = 0; i < 8; ++i) 
1366     {
1367       fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
1368     }
1369     for (Int_t i = 0; i < 3; ++i)
1370     {
1371       fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
1372       fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
1373     }
1374     fClusterizer->InitClusterUnfolding();
1375   }
1376   
1377   fClusterizer->SetDigitsArr(fDigitsArr);
1378   fClusterizer->SetOutput(0);
1379   fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1380   return kTRUE;
1381 }
1382
1383 //_____________________________________________________
1384 void AliEMCALTenderSupply::FillDigitsArray()
1385 {
1386   // Fill digits from cells to a TClonesArray.
1387   
1388   AliVEvent *event = GetEvent();
1389
1390  if (!event)
1391     return;
1392   
1393   fDigitsArr->Clear("C");
1394   AliVCaloCells *cells = event->GetEMCALCells();
1395   Int_t ncells = cells->GetNumberOfCells();
1396   for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) 
1397   {
1398     Double_t cellAmplitude=0, cellTime=0, efrac = 0;
1399     Short_t  cellNumber=0, mcLabel=-1;
1400
1401     if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime, mcLabel, efrac) != kTRUE)
1402       break;
1403
1404     // Do not add if already too low (some cells set to 0 if bad channels)
1405     if (cellAmplitude < fRecParam->GetMinECut() ) 
1406       continue;
1407
1408     // If requested, do not include exotic cells
1409    if (fEMCALRecoUtils->IsExoticCell(cellNumber,cells,event->GetBunchCrossNumber())) 
1410       continue;
1411         
1412     AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
1413     digit->SetId(cellNumber);
1414     digit->SetTime((Float_t)cellTime);
1415     digit->SetTimeR((Float_t)cellTime);
1416     digit->SetIndexInList(idigit);
1417     digit->SetType(AliEMCALDigit::kHG);
1418     digit->SetAmplitude((Float_t)cellAmplitude);
1419     idigit++;
1420   }
1421 }
1422
1423 //_____________________________________________________
1424 void AliEMCALTenderSupply::Clusterize()
1425 {
1426   // Clusterize.
1427   
1428   fClusterizer->Digits2Clusters("");
1429 }
1430
1431 //_____________________________________________________
1432 void AliEMCALTenderSupply::UpdateClusters()
1433 {
1434   // Update ESD cluster list.
1435   
1436   AliVEvent *event = GetEvent();
1437
1438   if (!event)
1439     return;
1440   
1441   TClonesArray *clus = dynamic_cast<TClonesArray*>(event->FindListObject("caloClusters"));
1442   if (!clus) 
1443     clus = dynamic_cast<TClonesArray*>(event->FindListObject("CaloClusters"));
1444   if (!clus) 
1445   {
1446     AliError(" Null pointer to calo clusters array, returning");
1447     return;
1448   }
1449   
1450   Int_t nents = clus->GetEntriesFast();
1451   for (Int_t i=0; i < nents; ++i) 
1452   {
1453     AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(i));
1454     if (!c)
1455       continue;
1456     if (c->IsEMCAL())
1457       delete clus->RemoveAt(i);
1458   }
1459   
1460   clus->Compress();
1461   
1462   RecPoints2Clusters(clus);
1463 }
1464
1465 //_____________________________________________________
1466 void AliEMCALTenderSupply::RecPoints2Clusters(TClonesArray *clus)
1467 {
1468   // Convert AliEMCALRecoPoints to AliESDCaloClusters/AliAODCaloClusters.
1469   // Cluster energy, global position, cells and their amplitude fractions are restored.
1470   
1471   AliVEvent *event = GetEvent();
1472
1473   if (!event)
1474     return;
1475
1476   Int_t ncls = fClusterArr->GetEntriesFast();
1477   for(Int_t i=0, nout=clus->GetEntriesFast(); i < ncls; ++i) 
1478   {
1479     AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
1480     
1481     Int_t ncellsTrue = 0;
1482     const Int_t ncells = recpoint->GetMultiplicity();
1483     UShort_t   absIds[ncells];  
1484     Double32_t ratios[ncells];
1485     Int_t *dlist = recpoint->GetDigitsList();
1486     Float_t *elist = recpoint->GetEnergiesList();
1487     for (Int_t c = 0; c < ncells; ++c) 
1488     {
1489       AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
1490       absIds[ncellsTrue] = digit->GetId();
1491       ratios[ncellsTrue] = elist[c]/digit->GetAmplitude();
1492       if (ratios[ncellsTrue] < 0.001) 
1493         continue;
1494       ++ncellsTrue;
1495     }
1496     
1497     if (ncellsTrue < 1) 
1498     {
1499       AliWarning("Skipping cluster with no cells");
1500       continue;
1501     }
1502     
1503     // calculate new cluster position
1504     TVector3 gpos;
1505     recpoint->GetGlobalPosition(gpos);
1506     Float_t g[3];
1507     gpos.GetXYZ(g);
1508     
1509     AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
1510     c->SetID(nout-1); 
1511     c->SetType(AliVCluster::kEMCALClusterv1);
1512     c->SetE(recpoint->GetEnergy());
1513     c->SetPosition(g);
1514     c->SetNCells(ncellsTrue);
1515     c->SetDispersion(recpoint->GetDispersion());
1516     c->SetEmcCpvDistance(-1);            //not yet implemented
1517     c->SetChi2(-1);                      //not yet implemented
1518     c->SetTOF(recpoint->GetTime()) ;     //time-of-flight
1519     c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
1520     Float_t elipAxis[2];
1521     recpoint->GetElipsAxis(elipAxis);
1522     c->SetM02(elipAxis[0]*elipAxis[0]) ;
1523     c->SetM20(elipAxis[1]*elipAxis[1]) ;
1524     c->SetCellsAbsId(absIds);
1525     c->SetCellsAmplitudeFraction(ratios);
1526   }
1527 }
1528
1529 //_____________________________________________________
1530 void AliEMCALTenderSupply::GetPass()
1531 {
1532   // Get passx from filename.
1533   
1534   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1535   fInputTree = mgr->GetTree();
1536   
1537   if (!fInputTree) 
1538   {
1539     AliError("Pointer to tree = 0, returning");
1540     return;
1541   }
1542   
1543   fInputFile = fInputTree->GetCurrentFile();
1544   if (!fInputFile) {
1545     AliError("Null pointer input file, returning");
1546     return;
1547   }
1548   
1549   TString fname(fInputFile->GetName());
1550   if      (fname.Contains("pass1")) fFilepass = TString("pass1");
1551   else if (fname.Contains("pass2")) fFilepass = TString("pass2");
1552   else if (fname.Contains("pass3")) fFilepass = TString("pass3");
1553   else 
1554   {
1555     AliError(Form("Pass number string not found: %s", fname.Data()));
1556     return;            
1557   }
1558 }