f87948e85b56d20f0ec9c3a9d6a3be591529a515
[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       Int_t fInitRecalib = InitRecalib();
497       if (fInitRecalib==0)
498         AliError("InitRecalib returned false, returning");
499       if (fInitRecalib==1)
500         AliWarning("InitRecalib OK");
501       if (fInitRecalib>1)
502         AliWarning(Form("No recalibration available: %d - %s", event->GetRunNumber(), fFilepass.Data()));
503         fReCalibCluster = kFALSE;
504     }
505     
506     // init time calibration
507     if( needTimecalib ){
508       Int_t initTC = InitTimeCalibration();
509       if ( !initTC ) 
510         AliError("InitTimeCalibration returned false, returning");
511       if (initTC==1)
512         AliWarning("InitTimeCalib OK");
513       if( initTC > 1 )
514         AliWarning(Form("No external time calibration set: %d - %s", event->GetRunNumber(), fFilepass.Data()));
515     }
516
517     // init misalignment matrix
518     if( needMisalign ) { 
519       if (!InitMisalignMatrix())
520         AliError("InitMisalignmentMatrix returned false, returning");
521       else
522         AliWarning("InitMisalignMatrix OK");
523     }
524     
525     // init clusterizer
526     if( needClusterizer ) {
527       if (!InitClusterization()) 
528         AliError("InitClusterization returned false, returning");
529       else
530         AliWarning("InitClusterization OK");
531     }
532     
533     if(fDebugLevel>1) 
534       fEMCALRecoUtils->Print("");
535   }
536   
537   // disable implied switches -------------------------------------------------
538   // AliEMCALRecoUtils or clusterizer functions alredy include some
539   // recalibration so based on those implied callibration te switches
540   // here are disabled to avoid duplication
541     
542   // clusterizer does cluster energy recalibration, position recomputation
543   // and shower shape
544   if( fReClusterize ){
545     fReCalibCluster   = kFALSE;
546     fRecalClusPos     = kFALSE;
547     fRecalShowerShape = kFALSE;
548   }
549   
550   // CONFIGURE THE RECO UTILS -------------------------------------------------
551   // configure the reco utils
552   // this option does energy recalibration
553   if( fCalibrateEnergy )
554     fEMCALRecoUtils->SwitchOnRecalibration();
555   else
556     fEMCALRecoUtils->SwitchOffRecalibration();
557   
558   // allows time calibration
559   if( fCalibrateTime )
560     fEMCALRecoUtils->SwitchOnTimeRecalibration();
561   else
562     fEMCALRecoUtils->SwitchOffTimeRecalibration();
563
564   // allows to zero bad cells
565   if( fBadCellRemove )
566     fEMCALRecoUtils->SwitchOnBadChannelsRemoval();
567   else
568     fEMCALRecoUtils->SwitchOffBadChannelsRemoval();
569   
570   // distance to bad channel recalibration
571   if( fRecalDistToBadChannels )
572     fEMCALRecoUtils->SwitchOnDistToBadChannelRecalculation();
573   else
574     fEMCALRecoUtils->SwitchOffDistToBadChannelRecalculation();
575
576   // exclude exotic cells
577   if( fRejectExoticCells )
578     fEMCALRecoUtils->SwitchOnRejectExoticCell();
579   else
580     fEMCALRecoUtils->SwitchOffRejectExoticCell();
581   
582   // exclude clusters with exotic cells
583   if( fRejectExoticClusters )
584     fEMCALRecoUtils->SwitchOnRejectExoticCluster();
585   else
586     fEMCALRecoUtils->SwitchOffRejectExoticCluster();
587
588   // TODO: not implemented switches
589   // SwitchOnClusterEnergySmearing
590   // SwitchOnRunDepCorrection
591
592   // START PROCESSING ---------------------------------------------------------
593   // Test if cells present
594   AliVCaloCells *cells= event->GetEMCALCells();
595   if (cells->GetNumberOfCells()<=0) 
596   {
597     if(fDebugLevel>1) 
598       AliWarning(Form("Number of EMCAL cells = %d, returning", cells->GetNumberOfCells()));
599     return;
600   }
601   
602   if (fDebugLevel>2)
603     AliInfo(Form("Re-calibrate cluster %d\n",fReCalibCluster));
604
605   // mark the cells not recalibrated in case of selected
606   // time, energy recalibration or bad channel removal
607   if( fCalibrateEnergy || fCalibrateTime || fBadCellRemove )
608     fEMCALRecoUtils->ResetCellsCalibrated();
609   
610  // CELL RECALIBRATION -------------------------------------------------------
611   // cell objects will be updated
612   // the cell calibrations are also processed locally any time those are needed
613   // in case that the cell objects are not to be updated here for later use
614   if( fUpdateCell )
615   {
616     // do the update
617     // includes exotic cell check in the UpdateCells function - is not provided
618     // by the reco utils
619     UpdateCells();
620
621     // switch off recalibrations so those are not done multiple times
622     // this is just for safety, the recalibrated flag of cell object
623     // should not allow for farther processing anyways
624     fEMCALRecoUtils->SwitchOffRecalibration();
625     fEMCALRecoUtils->SwitchOffTimeRecalibration();  
626   
627     if (fDoUpdateOnly)
628       return;
629   }
630
631   // RECLUSTERIZATION ---------------------------------------------------------
632   if (fReClusterize)
633   {
634     FillDigitsArray();
635     Clusterize();
636     UpdateClusters();
637   }
638
639   // Store good clusters
640   TClonesArray *clusArr = dynamic_cast<TClonesArray*>(event->FindListObject("caloClusters"));
641   if (!clusArr) 
642     clusArr = dynamic_cast<TClonesArray*>(event->FindListObject("CaloClusters"));
643   if (!clusArr) {
644     AliWarning(Form("No cluster array, number of cells in event = %d, returning", cells->GetNumberOfCells()));
645     return;
646   }
647
648   // PROCESS SINGLE CLUSTER RECALIBRATION -------------------------------------
649   // now go through clusters one by one and process separate correction
650   // as those were defined or not
651   Int_t nclusters = clusArr->GetEntriesFast();
652   for (Int_t icluster=0; icluster < nclusters; ++icluster) 
653   { 
654     AliVCluster *clust = static_cast<AliVCluster*>(clusArr->At(icluster));
655     if (!clust) 
656       continue;
657     if  (!clust->IsEMCAL()) 
658       continue;
659
660     // REMOVE CLUSTERS WITH BAD CELLS -----------------------------
661     if( fClusterBadChannelCheck )
662     {
663       // careful, the the ClusterContainsBadChannel is dependent on
664       // SwitchOnBadChannelsRemoval, switching it ON automatically
665       // and returning to original value after processing
666       Bool_t badRemoval = fEMCALRecoUtils->IsBadChannelsRemovalSwitchedOn();
667       fEMCALRecoUtils->SwitchOnBadChannelsRemoval();
668       
669       Bool_t badResult = fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo, clust->GetCellsAbsId(), clust->GetNCells());
670
671       // switch the bad channels removal back
672       if( ! badRemoval )
673         fEMCALRecoUtils->SwitchOffBadChannelsRemoval();
674       
675       if( badResult )
676       {
677         delete clusArr->RemoveAt(icluster);
678         continue; //TODO is it really needed to remove it? Or should we flag it?
679       }
680     }
681     
682     // REMOVE EXOTIC CLUSTERS -------------------------------------
683     // does process local cell recalibration energy and time without replacing
684     // the global cell values, in case of no cell recalib done yet
685     if( fRejectExoticClusters )
686     {
687       // careful, the IsExoticCluster is dependent on
688       // SwitchOnRejectExoticCell, switching it ON automatically
689       // and returning to original value after processing
690       Bool_t exRemoval = fEMCALRecoUtils->IsRejectExoticCell();
691       fEMCALRecoUtils->SwitchOnRejectExoticCell();
692
693       // get bunch crossing
694       Int_t bunchCrossNo = event->GetBunchCrossNumber();
695
696       Bool_t exResult = fEMCALRecoUtils->IsExoticCluster(clust, cells, bunchCrossNo );
697
698       // switch the exotic channels removal back
699       if( ! exRemoval )
700         fEMCALRecoUtils->SwitchOffRejectExoticCell();
701       
702       if( exResult )
703       {
704         delete clusArr->RemoveAt(icluster);
705         continue; //TODO is it really needed to remove it? Or should we flag it?
706       }
707     }
708     
709     // FIDUCIAL CUT -----------------------------------------------
710     if (fFiducial)
711     {
712       // depends on SetNumberOfCellsFromEMCALBorder
713       // SwitchOnNoFiducialBorderInEMCALEta0
714       if (!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, clust, cells)){
715         delete clusArr->RemoveAt(icluster);
716         continue; //TODO it would be nice to store the distance
717       }
718     }
719     
720     // CLUSTER ENERGY ---------------------------------------------
721     // does process local cell recalibration energy and time without replacing
722     // the global cell values, in case of no cell recalib done yet
723     if( fReCalibCluster ) 
724       fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, clust, cells);
725
726     // CLUSTER POSITION -------------------------------------------
727     // does process local cell energy recalibration, if enabled and cells
728     // not calibrated yet
729     if( fRecalClusPos ) 
730       fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, cells, clust);
731     
732     // SHOWER SHAPE -----------------------------------------------
733     if( fRecalShowerShape )
734       fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, cells, clust);  
735
736     // NONLINEARITY -----------------------------------------------
737     if( fDoNonLinearity )
738     {
739       Float_t correctedEnergy = fEMCALRecoUtils->CorrectClusterEnergyLinearity(clust);
740       clust->SetE(correctedEnergy);
741     }
742
743     // DISTANCE TO BAD CHANNELS -----------------------------------
744     if( fRecalDistToBadChannels )
745       fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, cells, clust);  
746   }
747
748   clusArr->Compress();
749
750   if (!fDoTrackMatch)
751     return;
752
753   // TRACK MATCHING -----------------------------------------------------------
754   if (!TGeoGlobalMagField::Instance()->GetField()) {
755     AliESDEvent *esd = dynamic_cast<AliESDEvent*>(event);
756     if (esd)
757       esd->InitMagneticField();
758     else {
759       AliAODEvent *aod = dynamic_cast<AliAODEvent*>(event);
760       Double_t curSol = 30000*aod->GetMagneticField()/5.00668;
761       Double_t curDip = 6000 *aod->GetMuonMagFieldScale();
762       AliMagF *field  = AliMagF::CreateFieldMap(curSol,curDip);
763       TGeoGlobalMagField::Instance()->SetField(field);
764     }
765   }
766
767   fEMCALRecoUtils->FindMatches(event,0x0,fEMCALGeo);
768   fEMCALRecoUtils->SetClusterMatchedToTrack(event);
769   fEMCALRecoUtils->SetTracksMatchedToCluster(event);
770 }
771
772 //_____________________________________________________
773 Bool_t AliEMCALTenderSupply::InitMisalignMatrix()
774 {
775   // Initialising misalignment matrices
776
777   AliVEvent *event = GetEvent();
778   
779   if (!event) 
780     return kFALSE;
781   
782   if (fGeomMatrixSet) 
783   {
784     AliInfo("Misalignment matrix already set");  
785     return kTRUE;
786   }
787   
788   if (fDebugLevel>0) 
789     AliInfo("Initialising misalignment matrix");  
790   
791   if (fLoadGeomMatrices) {
792     for(Int_t mod=0; mod < fEMCALGeo->GetNumberOfSuperModules(); ++mod)
793     {
794       if (fEMCALMatrix[mod]){
795         if(fDebugLevel > 2) 
796           fEMCALMatrix[mod]->Print();
797         fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod);  
798       }
799     }
800     fGeomMatrixSet = kTRUE;
801     return kTRUE;
802   }
803   
804   Int_t runGM = event->GetRunNumber();
805   TObjArray *mobj = 0;
806
807   if(fMisalignSurvey == kdefault)
808   { //take default alignment corresponding to run no
809     AliOADBContainer emcalgeoCont(Form("emcal"));
810     emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
811     mobj=(TObjArray*)emcalgeoCont.GetObject(runGM,"EmcalMatrices");
812   }
813   
814   if(fMisalignSurvey == kSurveybyS)
815   { //take alignment at sector level
816     if (runGM <= 140000) { //2010 data
817       AliOADBContainer emcalgeoCont(Form("emcal2010"));
818       emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
819       mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey10");
820     } 
821     else if (runGM>140000)
822     { // 2011 LHC11a pass1 data
823       AliOADBContainer emcalgeoCont(Form("emcal2011"));
824       emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
825       mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey11byS");      
826     }
827   }
828
829   if(fMisalignSurvey == kSurveybyM)
830   { //take alignment at module level
831     if (runGM <= 140000) { //2010 data
832       AliOADBContainer emcalgeoCont(Form("emcal2010"));
833       emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
834       mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey10");
835     } 
836     else if (runGM>140000) 
837     { // 2011 LHC11a pass1 data
838       AliOADBContainer emcalgeoCont(Form("emcal2011"));
839       emcalgeoCont.InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root",Form("AliEMCALgeo"));
840       mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey11byM");      
841     }
842   }
843
844   if(!mobj)
845   {
846     AliFatal("Geometry matrix array not found");
847     return kFALSE;
848   }
849   
850   for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
851   {
852     fEMCALMatrix[mod] = (TGeoHMatrix*) mobj->At(mod);
853     fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod); 
854     fEMCALMatrix[mod]->Print();
855   }
856   
857   return kTRUE;
858 }
859
860 //_____________________________________________________
861 Int_t AliEMCALTenderSupply::InitBadChannels()
862 {
863   // Initialising bad channel maps
864
865   AliVEvent *event = GetEvent();
866
867   if (!event) 
868     return 0;
869   
870   if (fDebugLevel>0) 
871     AliInfo("Initialising Bad channel map");
872   
873   // init default maps first
874   if( !fEMCALRecoUtils->GetEMCALBadChannelStatusMapArray() )
875     fEMCALRecoUtils->InitEMCALBadChannelStatusMap() ;
876   
877   Int_t runBC = event->GetRunNumber();
878   
879   AliOADBContainer *contBC = new AliOADBContainer("");
880   if (fBasePath!="")
881   { //if fBasePath specified in the ->SetBasePath()
882     if (fDebugLevel>0) AliInfo(Form("Loading Bad Channels OADB from given path %s",fBasePath.Data()));
883     
884     TFile *fbad=new TFile(Form("%s/EMCALBadChannels.root",fBasePath.Data()),"read");
885     if (!fbad || fbad->IsZombie())
886     {
887       AliFatal(Form("EMCALBadChannels.root was not found in the path provided: %s",fBasePath.Data()));
888       return 0;
889     }  
890     
891     if (fbad) delete fbad;
892     
893     contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fBasePath.Data()),"AliEMCALBadChannels");    
894   } 
895   else 
896   { // Else choose the one in the $ALICE_ROOT directory
897     if (fDebugLevel>0) AliInfo("Loading Bad Channels OADB from $ALICE_ROOT/OADB/EMCAL");
898     
899     TFile *fbad=new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root","read");
900     if (!fbad || fbad->IsZombie())
901     {
902       AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root was not found");
903       return 0;
904     }  
905       
906     if (fbad) delete fbad;
907     
908     contBC->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALBadChannels.root","AliEMCALBadChannels"); 
909   }
910   
911   TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runBC);
912   if (!arrayBC)
913   {
914     AliError(Form("No external hot channel set for run number: %d", runBC));
915     return 2; 
916   }
917
918   Int_t sms = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules();
919   for (Int_t i=0; i<sms; ++i) 
920   {
921     TH2I *h = fEMCALRecoUtils->GetEMCALChannelStatusMap(i);
922     if (h)
923       delete h;
924     h=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
925
926     if (!h) 
927     {
928       AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
929       continue;
930     }
931     h->SetDirectory(0);
932     fEMCALRecoUtils->SetEMCALChannelStatusMap(i,h);
933   }
934   return 1;  
935 }
936
937 //_____________________________________________________
938 Int_t AliEMCALTenderSupply::InitRecalib()
939 {
940   // Initialising recalibration factors.
941   
942   AliVEvent *event = GetEvent();
943
944   if (!event) 
945     return 0;
946   
947   if (fDebugLevel>0) 
948     AliInfo("Initialising recalibration factors");
949   
950   // init default maps first
951   if( !fEMCALRecoUtils->GetEMCALRecalibrationFactorsArray() )
952     fEMCALRecoUtils->InitEMCALRecalibrationFactors() ;
953
954   Int_t runRC = event->GetRunNumber();
955       
956   AliOADBContainer *contRF=new AliOADBContainer("");
957   if (fBasePath!="") 
958   { //if fBasePath specified in the ->SetBasePath()
959     if (fDebugLevel>0)  AliInfo(Form("Loading Recalib OADB from given path %s",fBasePath.Data()));
960     
961     TFile *fRecalib= new TFile(Form("%s/EMCALRecalib.root",fBasePath.Data()),"read");
962     if (!fRecalib || fRecalib->IsZombie()) 
963     {
964       AliFatal(Form("EMCALRecalib.root not found in %s",fBasePath.Data()));
965       return 0;
966     }
967     
968     if (fRecalib) delete fRecalib;
969     
970     contRF->InitFromFile(Form("%s/EMCALRecalib.root",fBasePath.Data()),"AliEMCALRecalib");
971   }
972   else
973   { // Else choose the one in the $ALICE_ROOT directory
974     if (fDebugLevel>0)  AliInfo("Loading Recalib OADB from $ALICE_ROOT/OADB/EMCAL");
975     
976     TFile *fRecalib= new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root","read");
977     if (!fRecalib || fRecalib->IsZombie()) 
978     {
979       AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root was not found");
980       return 0;
981     }
982     
983     if (fRecalib) delete fRecalib;
984       
985     contRF->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALRecalib.root","AliEMCALRecalib");     
986   }
987
988   TObjArray *recal=(TObjArray*)contRF->GetObject(runRC);
989   if (!recal)
990   {
991     AliError(Form("No Objects for run: %d",runRC));
992     return 2;
993   } 
994
995   TObjArray *recalpass=(TObjArray*)recal->FindObject(fFilepass);
996   if (!recalpass)
997   {
998     AliError(Form("No Objects for run: %d - %s",runRC,fFilepass.Data()));
999     return 2;
1000   }
1001
1002   TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
1003   if (!recalib)
1004   {
1005     AliError(Form("No Recalib histos found for  %d - %s",runRC,fFilepass.Data())); 
1006     return 2;
1007   }
1008
1009   if (fDebugLevel>0) recalib->Print();
1010
1011   Int_t sms = fEMCALGeo->GetEMCGeometry()->GetNumberOfSuperModules();
1012   for (Int_t i=0; i<sms; ++i) 
1013   {
1014     TH2F *h = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactors(i);
1015     if (h)
1016       delete h;
1017     h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
1018     if (!h) 
1019     {
1020       AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
1021       continue;
1022     }
1023     h->SetDirectory(0);
1024     fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
1025   }
1026   return 1;
1027 }
1028
1029 //_____________________________________________________
1030 Int_t AliEMCALTenderSupply::InitTimeCalibration()
1031 {
1032   // Initialising bad channel maps
1033   AliVEvent *event = GetEvent();
1034
1035   if (!event) 
1036     return 0;
1037   
1038   if (fDebugLevel>0) 
1039     AliInfo("Initialising time calibration map");
1040   
1041   // init default maps first
1042   if ( !fEMCALRecoUtils->GetEMCALTimeRecalibrationFactorsArray() )
1043     fEMCALRecoUtils->InitEMCALTimeRecalibrationFactors() ;
1044
1045   Int_t runBC = event->GetRunNumber();
1046   
1047   AliOADBContainer *contBC = new AliOADBContainer("");
1048   if (fBasePath!="")
1049   { //if fBasePath specified in the ->SetBasePath()
1050     if (fDebugLevel>0) AliInfo(Form("Loading time calibration OADB from given path %s",fBasePath.Data()));
1051     
1052     TFile *fbad=new TFile(Form("%s/EMCALTimeCalib.root",fBasePath.Data()),"read");
1053     if (!fbad || fbad->IsZombie())
1054     {
1055       AliFatal(Form("EMCALTimeCalib.root was not found in the path provided: %s",fBasePath.Data()));
1056       return 0;
1057     }  
1058     
1059     if (fbad) delete fbad;
1060     
1061     contBC->InitFromFile(Form("%s/EMCALTimeCalib.root",fBasePath.Data()),"AliEMCALTimeCalib");    
1062   } 
1063   else 
1064   { // Else choose the one in the $ALICE_ROOT directory
1065     if (fDebugLevel>0) AliInfo("Loading time calibration OADB from $ALICE_ROOT/OADB/EMCAL");
1066     
1067     TFile *fbad=new TFile("$ALICE_ROOT/OADB/EMCAL/EMCALTimeCalib.root","read");
1068     if (!fbad || fbad->IsZombie())
1069     {
1070       AliFatal("$ALICE_ROOT/OADB/EMCAL/EMCALTimeCalib.root was not found");
1071       return 0;
1072     }  
1073       
1074     if (fbad) delete fbad;
1075     
1076     contBC->InitFromFile("$ALICE_ROOT/OADB/EMCAL/EMCALTimeCalib.root","AliEMCALTimeCalib"); 
1077   }
1078   
1079   TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runBC);
1080   if (!arrayBC)
1081   {
1082     AliError(Form("No external time calibration set for run number: %d", runBC));
1083     return 2; 
1084   }
1085   
1086   // Here, it looks for a specific pass
1087   TObjArray *arrayBCpass=(TObjArray*)arrayBC->FindObject(fFilepass); 
1088   if (!arrayBCpass)
1089   {
1090     AliError(Form("No external time calibration set for: %d -%s", runBC,fFilepass.Data()));
1091     return 2; 
1092   }
1093
1094   if (fDebugLevel>0) arrayBCpass->Print();
1095
1096   for( Int_t i = 0; i < 4; i++ )
1097   {
1098     TH1F *h = fEMCALRecoUtils->GetEMCALChannelTimeRecalibrationFactors( i );
1099     if( h )
1100       delete h;
1101     
1102     h = (TH1F*)arrayBCpass->FindObject(Form("hAllTimeAvBC%d",i));
1103     
1104     if (!h)
1105     {
1106       AliError(Form("Can not get hAllTimeAvBC%d",i));
1107       continue;
1108     }
1109     h->SetDirectory(0);
1110     fEMCALRecoUtils->SetEMCALChannelTimeRecalibrationFactors(i,h);
1111   }
1112   return 1;  
1113 }
1114
1115 //_____________________________________________________
1116 void AliEMCALTenderSupply::UpdateCells()
1117 {
1118   //Remove bad cells from the cell list
1119   //Recalibrate energy and time cells 
1120   //This is required for later reclusterization
1121
1122   AliVEvent *event = GetEvent();
1123
1124   if(!event) return ;
1125   
1126   AliVCaloCells *cells = event->GetEMCALCells();
1127   Int_t bunchCrossNo = event->GetBunchCrossNumber();
1128
1129   fEMCALRecoUtils->RecalibrateCells(cells, bunchCrossNo); 
1130   
1131   // remove exotic cells - loop through cells and zero the exotic ones
1132   // just like with bad cell rejection in reco utils (inside RecalibrateCells)
1133   if( fRejectExoticCells )
1134   {
1135     Short_t  absId  =-1;
1136     Double_t ecell = 0;
1137     Double_t tcell = 0;
1138     Double_t efrac = 0;
1139     Short_t  mclabel = -1;
1140     Bool_t   isExot = kFALSE;
1141   
1142     // loop through cells
1143     Int_t nEMcell  = cells->GetNumberOfCells() ;  
1144     for (Int_t iCell = 0; iCell < nEMcell; iCell++) 
1145     { 
1146       cells->GetCell( iCell, absId, ecell, tcell, mclabel, efrac );
1147     
1148       isExot = fEMCALRecoUtils->IsExoticCell( absId, cells, bunchCrossNo ); 
1149       // zero if exotic
1150       if( isExot )
1151         cells->SetCell( iCell, absId, 0.0, -1.0, mclabel, efrac );
1152     } // cell loop
1153   } // reject exotic cells
1154
1155   cells->Sort();
1156 }
1157
1158 //_____________________________________________________
1159 TString AliEMCALTenderSupply::GetBeamType()
1160 {
1161   // Get beam type : pp-AA-pA
1162   // ESDs have it directly, AODs get it from hardcoded run number ranges
1163   
1164   AliVEvent *event = GetEvent();
1165
1166   if (!event) { 
1167     AliError("Couldn't retrieve event!");
1168     return "";
1169   }
1170
1171   TString beamType;
1172
1173   AliESDEvent *esd = dynamic_cast<AliESDEvent*>(event);
1174   if (esd) {
1175     const AliESDRun *run = esd->GetESDRun();
1176     beamType = run->GetBeamType();
1177   }
1178   else
1179   {
1180     Int_t runNumber = event->GetRunNumber();
1181     if ((runNumber >= 136851 && runNumber <= 139517)  // LHC10h
1182   || (runNumber >= 166529 && runNumber <= 170593))  // LHC11h
1183     {
1184       beamType = "A-A";
1185     }
1186     else 
1187     {
1188       beamType = "p-p";
1189     }
1190   }
1191
1192   return beamType;    
1193 }
1194
1195 //_____________________________________________________
1196 Int_t AliEMCALTenderSupply::InitRecParam()
1197 {
1198   // Init reco params if not yet exist (probably shipped by the user already)
1199
1200   if( fRecParam != 0 )
1201     return 2;
1202
1203   TString beamType = GetBeamType();
1204
1205   // set some default reco params
1206   fRecParam = new AliEMCALRecParam();
1207   fRecParam->SetClusteringThreshold(0.100);
1208   fRecParam->SetMinECut(0.050);
1209   fRecParam->SetTimeCut(250);
1210   fRecParam->SetTimeMin(425);
1211   fRecParam->SetTimeMax(825);
1212   if ( beamType == "A-A") {
1213     fRecParam->SetClusterizerFlag(AliEMCALRecParam::kClusterizerv2);
1214   } 
1215   else {
1216     fRecParam->SetClusterizerFlag(AliEMCALRecParam::kClusterizerv1);
1217   }
1218
1219   return 0;
1220 }
1221
1222 //_____________________________________________________
1223 Bool_t AliEMCALTenderSupply::InitClusterization()
1224 {
1225   // Initialing clusterization/unfolding algorithm and set all the needed parameters.
1226   
1227   AliVEvent *event = GetEvent();
1228
1229   if (!event) 
1230     return kFALSE;
1231   
1232   if (fDebugLevel>0) 
1233     AliInfo(Form("Initialising reclustering parameters: Clusterizer type: %d",fRecParam->GetClusterizerFlag()));
1234   
1235   //---setup clusterizer
1236   delete fClusterizer;
1237   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1238     fClusterizer = new AliEMCALClusterizerv1 (fEMCALGeo);
1239   else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2) 
1240     fClusterizer = new AliEMCALClusterizerv2(fEMCALGeo);
1241   else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) 
1242   {
1243     AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fEMCALGeo);
1244     clusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
1245     clusterizer->SetNColDiff(fRecParam->GetNColDiff());
1246     fClusterizer = clusterizer;
1247   } 
1248   else 
1249   {
1250     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1251     return kFALSE;
1252   }
1253   
1254   // Set the clustering parameters
1255   fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
1256   fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
1257   fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );    
1258   fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
1259   fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
1260   fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
1261   fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
1262   fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
1263   fClusterizer->SetInputCalibrated       ( kTRUE                               );
1264   fClusterizer->SetJustClusters          ( kTRUE                               );  
1265   
1266   // In case of unfolding after clusterization is requested, set the corresponding parameters
1267   if (fRecParam->GetUnfold()) 
1268   {
1269     for (Int_t i = 0; i < 8; ++i) 
1270     {
1271       fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
1272     }
1273     for (Int_t i = 0; i < 3; ++i)
1274     {
1275       fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
1276       fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
1277     }
1278     fClusterizer->InitClusterUnfolding();
1279   }
1280   
1281   fClusterizer->SetDigitsArr(fDigitsArr);
1282   fClusterizer->SetOutput(0);
1283   fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1284   return kTRUE;
1285 }
1286
1287 //_____________________________________________________
1288 void AliEMCALTenderSupply::FillDigitsArray()
1289 {
1290   // Fill digits from cells to a TClonesArray.
1291   
1292   AliVEvent *event = GetEvent();
1293
1294  if (!event)
1295     return;
1296   
1297   fDigitsArr->Clear("C");
1298   AliVCaloCells *cells = event->GetEMCALCells();
1299   Int_t ncells = cells->GetNumberOfCells();
1300   for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) 
1301   {
1302     Double_t cellAmplitude=0, cellTime=0, efrac = 0;
1303     Short_t  cellNumber=0, mcLabel=-1;
1304
1305     if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime, mcLabel, efrac) != kTRUE)
1306       break;
1307
1308     // Do not add if already too low (some cells set to 0 if bad channels)
1309     if (cellAmplitude < fRecParam->GetMinECut() ) 
1310       continue;
1311
1312     // If requested, do not include exotic cells
1313    if (fEMCALRecoUtils->IsExoticCell(cellNumber,cells,event->GetBunchCrossNumber())) 
1314       continue;
1315         
1316     AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
1317     digit->SetId(cellNumber);
1318     digit->SetTime((Float_t)cellTime);
1319     digit->SetTimeR((Float_t)cellTime);
1320     digit->SetIndexInList(idigit);
1321     digit->SetType(AliEMCALDigit::kHG);
1322     digit->SetAmplitude((Float_t)cellAmplitude);
1323     idigit++;
1324   }
1325 }
1326
1327 //_____________________________________________________
1328 void AliEMCALTenderSupply::Clusterize()
1329 {
1330   // Clusterize.
1331   
1332   fClusterizer->Digits2Clusters("");
1333 }
1334
1335 //_____________________________________________________
1336 void AliEMCALTenderSupply::UpdateClusters()
1337 {
1338   // Update ESD cluster list.
1339   
1340   AliVEvent *event = GetEvent();
1341
1342   if (!event)
1343     return;
1344   
1345   TClonesArray *clus = dynamic_cast<TClonesArray*>(event->FindListObject("caloClusters"));
1346   if (!clus) 
1347     clus = dynamic_cast<TClonesArray*>(event->FindListObject("CaloClusters"));
1348   if (!clus) 
1349   {
1350     AliError(" Null pointer to calo clusters array, returning");
1351     return;
1352   }
1353   
1354   Int_t nents = clus->GetEntriesFast();
1355   for (Int_t i=0; i < nents; ++i) 
1356   {
1357     AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(i));
1358     if (!c)
1359       continue;
1360     if (c->IsEMCAL())
1361       delete clus->RemoveAt(i);
1362   }
1363   
1364   clus->Compress();
1365   
1366   RecPoints2Clusters(clus);
1367 }
1368
1369 //_____________________________________________________
1370 void AliEMCALTenderSupply::RecPoints2Clusters(TClonesArray *clus)
1371 {
1372   // Convert AliEMCALRecoPoints to AliESDCaloClusters/AliAODCaloClusters.
1373   // Cluster energy, global position, cells and their amplitude fractions are restored.
1374   
1375   AliVEvent *event = GetEvent();
1376
1377   if (!event)
1378     return;
1379
1380   Int_t ncls = fClusterArr->GetEntriesFast();
1381   for(Int_t i=0, nout=clus->GetEntriesFast(); i < ncls; ++i) 
1382   {
1383     AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
1384     
1385     Int_t ncellsTrue = 0;
1386     const Int_t ncells = recpoint->GetMultiplicity();
1387     UShort_t   absIds[ncells];  
1388     Double32_t ratios[ncells];
1389     Int_t *dlist = recpoint->GetDigitsList();
1390     Float_t *elist = recpoint->GetEnergiesList();
1391     for (Int_t c = 0; c < ncells; ++c) 
1392     {
1393       AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
1394       absIds[ncellsTrue] = digit->GetId();
1395       ratios[ncellsTrue] = elist[c]/digit->GetAmplitude();
1396       if (ratios[ncellsTrue] < 0.001) 
1397         continue;
1398       ++ncellsTrue;
1399     }
1400     
1401     if (ncellsTrue < 1) 
1402     {
1403       AliWarning("Skipping cluster with no cells");
1404       continue;
1405     }
1406     
1407     // calculate new cluster position
1408     TVector3 gpos;
1409     recpoint->GetGlobalPosition(gpos);
1410     Float_t g[3];
1411     gpos.GetXYZ(g);
1412     
1413     AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
1414     c->SetID(nout-1); 
1415     c->SetType(AliVCluster::kEMCALClusterv1);
1416     c->SetE(recpoint->GetEnergy());
1417     c->SetPosition(g);
1418     c->SetNCells(ncellsTrue);
1419     c->SetDispersion(recpoint->GetDispersion());
1420     c->SetEmcCpvDistance(-1);            //not yet implemented
1421     c->SetChi2(-1);                      //not yet implemented
1422     c->SetTOF(recpoint->GetTime()) ;     //time-of-flight
1423     c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
1424     Float_t elipAxis[2];
1425     recpoint->GetElipsAxis(elipAxis);
1426     c->SetM02(elipAxis[0]*elipAxis[0]) ;
1427     c->SetM20(elipAxis[1]*elipAxis[1]) ;
1428     c->SetCellsAbsId(absIds);
1429     c->SetCellsAmplitudeFraction(ratios);
1430   }
1431 }
1432
1433 //_____________________________________________________
1434 void AliEMCALTenderSupply::GetPass()
1435 {
1436   // Get passx from filename.
1437   
1438   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1439   fInputTree = mgr->GetTree();
1440   
1441   if (!fInputTree) 
1442   {
1443     AliError("Pointer to tree = 0, returning");
1444     return;
1445   }
1446   
1447   fInputFile = fInputTree->GetCurrentFile();
1448   if (!fInputFile) {
1449     AliError("Null pointer input file, returning");
1450     return;
1451   }
1452   
1453   TString fname(fInputFile->GetName());
1454   if      (fname.Contains("pass1")) fFilepass = TString("pass1");
1455   else if (fname.Contains("pass2")) fFilepass = TString("pass2");
1456   else if (fname.Contains("pass3")) fFilepass = TString("pass3");
1457   else 
1458   {
1459     AliError(Form("Pass number string not found: %s", fname.Data()));
1460     return;            
1461   }
1462 }