]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGGA/EMCALTasks/AliAnalysisTaskEMCALClusterize.cxx
submitted AddTasks for GG pPb
[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliAnalysisTaskEMCALClusterize.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 // This analysis provides a new list of clusters to be used in other analysis
18 //
19 // Author: Gustavo Conesa Balbastre,
20 //         Adapted from analysis class from Deepa Thomas
21 //
22 // $Id$
23 //_________________________________________________________________________
24
25 // --- Root ---
26 #include <TString.h>
27 #include <TRefArray.h>
28 #include <TClonesArray.h>
29 #include <TTree.h>
30 #include <TGeoManager.h>
31 #include <TROOT.h>
32 #include <TInterpreter.h>
33 #include <TFile.h>
34
35 // --- AliRoot Analysis Steering
36 #include "AliAnalysisTask.h"
37 #include "AliAnalysisManager.h"
38 #include "AliESDEvent.h"
39 #include "AliGeomManager.h"
40 #include "AliVCaloCells.h"
41 #include "AliAODCaloCluster.h"
42 #include "AliCDBManager.h"
43 #include "AliCDBStorage.h"
44 #include "AliCDBEntry.h"
45 #include "AliLog.h"
46 #include "AliVEventHandler.h"
47 #include "AliAODInputHandler.h"
48 #include "AliOADBContainer.h"
49 #include "AliAODMCParticle.h"
50
51 // --- EMCAL
52 #include "AliEMCALAfterBurnerUF.h"
53 #include "AliEMCALGeometry.h"
54 #include "AliEMCALClusterizerNxN.h"
55 #include "AliEMCALClusterizerv1.h"
56 #include "AliEMCALClusterizerv2.h"
57 #include "AliEMCALRecPoint.h"
58 #include "AliEMCALDigit.h"
59 #include "AliCaloCalibPedestal.h"
60 #include "AliEMCALCalibData.h"
61
62 #include "AliAnalysisTaskEMCALClusterize.h"
63
64 ClassImp(AliAnalysisTaskEMCALClusterize)
65
66 //______________________________________________________________________________
67 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name) 
68 : AliAnalysisTaskSE(name)
69 , fEvent(0)
70 , fGeom(0),               fGeomName("") 
71 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
72 , fCalibData(0),          fPedestalData(0)
73 , fOCDBpath(""),          fAccessOCDB(kFALSE)
74 , fDigitsArr(0),          fClusterArr(0),             fCaloClusterArr(0)
75 , fRecParam(0),           fClusterizer(0)
76 , fUnfolder(0),           fJustUnfold(kFALSE) 
77 , fOutputAODBranch(0),    fOutputAODBranchName(""),   fOutputAODBranchSet(0)
78 , fFillAODFile(kFALSE),   fFillAODHeader(0)
79 , fFillAODCaloCells(0),   fRun(-1)
80 , fRecoUtils(0),          fConfigName("")
81 , fOrgClusterCellId()
82 , fCellLabels(),          fCellSecondLabels(),        fCellTime()
83 , fCellMatchdEta(),       fCellMatchdPhi()
84 , fRecalibrateWithClusterTime(0)
85 , fMaxEvent(0),           fDoTrackMatching(kFALSE)
86 , fSelectCell(kFALSE),    fSelectCellMinE(0),         fSelectCellMinFrac(0)
87 , fRemoveLEDEvents(kTRUE),fRemoveExoticEvents(kFALSE)
88 , fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("") 
89 , fOADBSet(kFALSE),       fAccessOADB(kTRUE),         fOADBFilePath("")
90 , fCentralityClass(""),   fSelectEMCALEvent(0)
91 , fEMCALEnergyCut(0.),    fEMCALNcellsCut (0)
92 , fSetCellMCLabelFromCluster(0)
93 , fRemapMCLabelForAODs(0)
94 , fInputFromFilter(0)
95 {
96   // Constructor
97   
98   for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
99   
100   ResetArrays();
101   
102   fCentralityBin[0] = fCentralityBin[1]=-1;
103   
104 }
105
106 //______________________________________________________________
107 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize() 
108 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
109 , fEvent(0)
110 , fGeom(0),                 fGeomName("") 
111 , fGeomMatrixSet(kFALSE),   fLoadGeomMatrices(kFALSE)
112 , fCalibData(0),            fPedestalData(0)
113 , fOCDBpath(""),            fAccessOCDB(kFALSE)
114 , fDigitsArr(0),            fClusterArr(0),             fCaloClusterArr(0)
115 , fRecParam(0),             fClusterizer(0)
116 , fUnfolder(0),             fJustUnfold(kFALSE) 
117 , fOutputAODBranch(0),      fOutputAODBranchName(""),   fOutputAODBranchSet(0)
118 , fFillAODFile(kFALSE),     fFillAODHeader(0)
119 , fFillAODCaloCells(0),     fRun(-1)
120 , fRecoUtils(0),            fConfigName("")
121 , fOrgClusterCellId()
122 , fCellLabels(),            fCellSecondLabels(),        fCellTime()
123 , fCellMatchdEta(),         fCellMatchdPhi()
124 , fRecalibrateWithClusterTime(0)
125 , fMaxEvent(0),             fDoTrackMatching(kFALSE)
126 , fSelectCell(kFALSE),      fSelectCellMinE(0),         fSelectCellMinFrac(0)
127 , fRemoveLEDEvents(kTRUE),  fRemoveExoticEvents(kFALSE)
128 , fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("")
129 , fOADBSet(kFALSE),         fAccessOADB(kTRUE),        fOADBFilePath("")
130 , fCentralityClass(""),     fSelectEMCALEvent(0)
131 , fEMCALEnergyCut(0.),      fEMCALNcellsCut (0)
132 , fSetCellMCLabelFromCluster(0)
133 , fRemapMCLabelForAODs(0)
134 , fInputFromFilter(0)
135 {
136   // Constructor
137   
138   for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
139   
140   ResetArrays();
141   
142   fCentralityBin[0] = fCentralityBin[1]=-1;
143
144 }
145
146
147 //_______________________________________________________________
148 AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
149 {
150   //dtor 
151   
152   if (fDigitsArr)
153   {
154     fDigitsArr->Clear("C");
155     delete fDigitsArr; 
156   }
157   
158   if (fClusterArr)
159   {
160     fClusterArr->Delete();
161     delete fClusterArr;
162   }
163   
164   if (fCaloClusterArr)
165   {
166     fCaloClusterArr->Delete();
167     delete fCaloClusterArr; 
168   }
169   
170   if(fClusterizer) delete fClusterizer;
171   if(fUnfolder)    delete fUnfolder;   
172   if(fRecoUtils)   delete fRecoUtils;
173   
174 }
175
176 //_______________________________________________________
177 Bool_t AliAnalysisTaskEMCALClusterize::AcceptEventEMCAL()
178 {
179   // Accept event given there is a EMCAL cluster with enough energy, and not noisy, exotic
180   
181   if(!fSelectEMCALEvent)   return kTRUE; // accept
182   
183   if(fEMCALEnergyCut <= 0) return kTRUE; // accept
184   
185   Int_t           nCluster = fEvent -> GetNumberOfCaloClusters();
186   AliVCaloCells * caloCell = fEvent -> GetEMCALCells();
187   Int_t           bc       = fEvent -> GetBunchCrossNumber();
188   
189   for(Int_t icalo = 0; icalo < nCluster; icalo++)
190   {
191     AliVCluster *clus = (AliVCluster*) (fEvent->GetCaloCluster(icalo));
192     
193     if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fEMCALNcellsCut ) && ( clus->E() > fEMCALEnergyCut ) &&
194        fRecoUtils->IsGoodCluster(clus,fGeom,caloCell,bc))
195     {
196       
197       if (fDebug > 0)
198         printf("AliAnalysisTaskEMCALClusterize::AcceptEventEMCAL() - Accept :  E %2.2f > %2.2f, nCells %d > %d \n",
199                clus->E(), fEMCALEnergyCut, clus->GetNCells(), fEMCALNcellsCut);
200       
201       return kTRUE;
202     }
203     
204   }// loop
205   
206   if (fDebug > 0)
207     printf("AliAnalysisTaskEMCALClusterize::AcceptEventEMCAL() - Reject \n");
208   
209   return kFALSE;
210   
211 }
212
213 //_______________________________________________
214 void AliAnalysisTaskEMCALClusterize::AccessOADB()
215 {
216   // Set the AODB calibration, bad channels etc. parameters at least once
217   // alignment matrices from OADB done in SetGeometryMatrices
218   
219   //Set it only once
220   if(fOADBSet) return ; 
221   
222   Int_t   runnumber = InputEvent()->GetRunNumber() ;
223   TString pass      = GetPass();
224   
225   printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Get AODB parameters from EMCAL in %s for run %d, and <%s> \n",fOADBFilePath.Data(),runnumber,pass.Data());
226   
227   Int_t nSM = fGeom->GetNumberOfSuperModules();
228   
229   // Bad map
230   if(fRecoUtils->IsBadChannelsRemovalSwitchedOn())
231   {
232     AliOADBContainer *contBC=new AliOADBContainer("");
233     contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePath.Data()),"AliEMCALBadChannels"); 
234     
235     TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber);
236     
237     if(arrayBC)
238     {
239         fRecoUtils->SwitchOnDistToBadChannelRecalculation();
240         printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Remove EMCAL bad cells \n");
241         
242         for (Int_t i=0; i<nSM; ++i) 
243         {
244           TH2I *hbm = fRecoUtils->GetEMCALChannelStatusMap(i);
245           
246           if (hbm)
247             delete hbm;
248           
249           hbm=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
250           
251           if (!hbm) 
252           {
253             AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
254             continue;
255           }
256           
257           hbm->SetDirectory(0);
258           fRecoUtils->SetEMCALChannelStatusMap(i,hbm);
259           
260         } // loop
261     } else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT remove EMCAL bad channels\n"); // run array
262   }  // Remove bad
263   
264   // Energy Recalibration
265   if(fRecoUtils->IsRecalibrationOn())
266   {
267     AliOADBContainer *contRF=new AliOADBContainer("");
268     
269     contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePath.Data()),"AliEMCALRecalib");
270     
271     TObjArray *recal=(TObjArray*)contRF->GetObject(runnumber); 
272     
273     if(recal)
274     {
275       TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
276       
277       if(recalpass)
278       {
279         TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
280         
281         if(recalib)
282         {
283           printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Recalibrate EMCAL \n");
284           for (Int_t i=0; i<nSM; ++i) 
285           {
286             TH2F *h = fRecoUtils->GetEMCALChannelRecalibrationFactors(i);
287             
288             if (h)
289               delete h;
290             
291             h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
292             
293             if (!h) 
294             {
295               AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
296               continue;
297             }
298             
299             h->SetDirectory(0);
300             
301             fRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
302           } // SM loop
303         }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params object array \n"); // array ok
304       }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for pass\n"); // array pass ok
305     }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for run\n");  // run number array ok
306         
307   } // Recalibration on
308   
309   // Energy Recalibration, apply on top of previous calibration factors
310   if(fRecoUtils->IsRunDepRecalibrationOn())
311   {
312     AliOADBContainer *contRFTD=new AliOADBContainer("");
313     
314     contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePath.Data()),"AliEMCALRunDepTempCalibCorrections");
315     
316     TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber); 
317     
318     //If it did not exist for this run, get closes one
319     if (!htd)
320     {
321       AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runnumber));
322       // let's get the closest runnumber instead then..
323       Int_t lower = 0;
324       Int_t ic = 0;
325       Int_t maxEntry = contRFTD->GetNumberOfEntries();
326       
327       while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < runnumber) ) {
328         lower = ic;
329         ic++;
330       }
331       
332       Int_t closest = lower;
333       if ( (ic<maxEntry) &&
334           (contRFTD->LowerLimit(ic)-runnumber) < (runnumber - contRFTD->UpperLimit(lower)) ) {
335         closest = ic;
336       }
337       
338       AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRFTD->LowerLimit(closest)));
339       htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
340     } 
341
342     if(htd)
343     {
344       printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Recalibrate (Temperature) EMCAL \n");
345       
346       for (Int_t ism=0; ism<nSM; ++ism) 
347       {        
348         for (Int_t icol=0; icol<48; ++icol) 
349         {        
350           for (Int_t irow=0; irow<24; ++irow) 
351           {
352             Float_t factor = fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
353             
354             Int_t absID = fGeom->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
355             factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
356             //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID, 
357             //      GetEMCALChannelRecalibrationFactor(ism,icol,irow) , htd->GetBinContent(absID) / 10000., factor);
358             fRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
359           } // columns
360         } // rows 
361       } // SM loop
362     }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL with T variations, no params TH1 \n"); 
363   } // Run by Run T calibration
364   
365   // Time Recalibration
366   if(fRecoUtils->IsTimeRecalibrationOn())
367   {
368     AliOADBContainer *contTRF=new AliOADBContainer("");
369     
370     contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePath.Data()),"AliEMCALTimeCalib");
371     
372     TObjArray *trecal=(TObjArray*)contTRF->GetObject(runnumber); 
373     
374     if(trecal)
375     {
376       TObjArray *trecalpass=(TObjArray*)trecal->FindObject(pass);
377
378       if(trecalpass)
379       {
380         printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Time Recalibrate EMCAL \n");
381         for (Int_t ibc = 0; ibc < 4; ++ibc) 
382         {
383           TH1F *h = fRecoUtils->GetEMCALChannelTimeRecalibrationFactors(ibc);
384           
385           if (h)
386             delete h;
387           
388           h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
389           
390           if (!h) 
391           {
392             AliError(Form("Could not load hAllTimeAvBC%d",ibc));
393             continue;
394           }
395           
396           h->SetDirectory(0);
397           
398           fRecoUtils->SetEMCALChannelTimeRecalibrationFactors(ibc,h);
399         } // bunch crossing loop
400       }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for pass\n"); // array pass ok
401     }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for run\n");  // run number array ok
402     
403   } // Time recalibration on    
404   
405   // Parameters already set once, so do not it again
406   fOADBSet = kTRUE;
407   
408 }  
409
410 //_________________________________________________
411 Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
412 {
413   //Access to OCDB stuff
414   
415   fEvent = InputEvent();
416   if (!fEvent)
417   {
418     Warning("AccessOCDB","Event not available!!!");
419     return kFALSE;
420   }
421   
422   if (fEvent->GetRunNumber()==fRun)
423     return kTRUE;
424   fRun = fEvent->GetRunNumber();
425   
426   if(DebugLevel() > 1 )
427     printf("AliAnalysisTaksEMCALClusterize::AccessOCDB() - Begin");
428   
429   AliCDBManager *cdb = AliCDBManager::Instance();
430   
431   
432   if (fOCDBpath.Length())
433   {
434     cdb->SetDefaultStorage(fOCDBpath.Data());
435     printf("AliAnalysisTaksEMCALClusterize::AccessOCDB() - Default storage %s",fOCDBpath.Data());
436   }
437   
438   cdb->SetRun(fEvent->GetRunNumber());
439   
440   //
441   // EMCAL from RAW OCDB
442   if (fOCDBpath.Contains("alien:"))
443   {
444     cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
445     cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
446   }
447   
448   TString path = cdb->GetDefaultStorage()->GetBaseFolder();
449   
450   // init parameters:
451   
452   //Get calibration parameters  
453   if(!fCalibData)
454   {
455     AliCDBEntry *entry = (AliCDBEntry*) 
456     AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
457     
458     if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
459   }
460   
461   if(!fCalibData)
462     AliFatal("Calibration parameters not found in CDB!");
463   
464   //Get calibration parameters  
465   if(!fPedestalData)
466   {
467     AliCDBEntry *entry = (AliCDBEntry*) 
468     AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
469     
470     if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
471   }
472   
473   if(!fPedestalData)
474     AliFatal("Dead map not found in CDB!");
475   
476   return kTRUE;
477 }
478
479 //_____________________________________________________
480 void AliAnalysisTaskEMCALClusterize::CheckAndGetEvent()
481 {
482   // Get the input event, it can depend in embedded events what you want to get
483   // Also check if the quality of the event is good if not reject it
484   
485   fEvent = 0x0;
486       
487   AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
488   Int_t eventN = Entry();
489   if(aodIH) eventN = aodIH->GetReadEntry(); 
490   
491   if (eventN > fMaxEvent) 
492     return ;
493   
494   //printf("Clusterizer --- Event %d-- \n",eventN);
495   
496   //Check if input event are embedded events
497   //If so, take output event
498   if (aodIH && aodIH->GetMergeEvents()) 
499   {
500     fEvent  = AODEvent();
501     
502     if(!aodIH->GetMergeEMCALCells()) 
503       AliFatal("Events merged but not EMCAL cells, check analysis settings!");
504     
505     if(DebugLevel() > 1)
506     {
507       printf("AliAnalysisTaksEMCALClusterize::UserExec() - Use embedded events\n");
508       
509       printf("\t InputEvent  N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
510              InputEvent()->GetEMCALCells()->GetNumberOfCells());
511       
512       printf("\t MergedEvent  N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(),
513              aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells());
514       
515       for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++) 
516       {
517         AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
518         if(sigCluster->IsEMCAL()) printf("\t \t Signal cluster: i %d, E  %f\n",icl,sigCluster->E());
519       }
520       
521       printf("\t OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
522              AODEvent()->GetEMCALCells()->GetNumberOfCells());
523     }
524   }
525   else if(fInputFromFilter)
526   {
527     //printf("Get Input From Filtered AOD\n");
528     fEvent =  AODEvent();
529   }
530   else 
531   {
532     fEvent =  InputEvent();
533     if(fFillAODCaloCells) FillAODCaloCells();   
534     if(fFillAODHeader)    FillAODHeader();
535   }
536   
537   if (!fEvent) 
538   {
539     Error("UserExec","Event not available");
540     return ;
541   }
542   
543   //Process events if there is a high energy cluster
544   if(!AcceptEventEMCAL())  { fEvent = 0x0 ; return ; }
545   
546   //-------------------------------------------------------------------------------------
547   // Reject events if LED was firing, use only for LHC11a data 
548   // Reject event if triggered by exotic cell and remove exotic cells if not triggered
549   //-------------------------------------------------------------------------------------
550   
551   if( IsLEDEvent( InputEvent()->GetRunNumber() ) ) { fEvent = 0x0 ; return ; }
552   
553   if( IsExoticEvent() )                            { fEvent = 0x0 ; return ; }
554   
555   //-------------------------------------------------------------------------------------
556   // Set the cluster array in the event (output or input)
557   //-------------------------------------------------------------------------------------
558   
559   if     ( fFillAODFile ) 
560   {
561     //Magic line to write events to AOD filem put after event rejection
562     AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
563   }
564   else if( !fOutputAODBranchSet )
565   {
566     // Create array and put it in the input event, if output AOD not selected, only once
567     InputEvent()->AddObject(fOutputAODBranch);
568     fOutputAODBranchSet = kTRUE;
569     printf("AliAnalysisTaskEMCALClusterize::UserExec() - Add AOD branch <%s> to input event\n",fOutputAODBranchName.Data());
570   }
571   
572 }
573
574 //____________________________________________________
575 void AliAnalysisTaskEMCALClusterize::ClusterizeCells()
576
577   // Recluster calocells, transform them into digits, 
578   // feed the clusterizer with them and get new list of clusters
579   
580   //In case of MC, first loop on the clusters and fill MC label to array  
581   Int_t nClusters     = fEvent->GetNumberOfCaloClusters();
582   Int_t nClustersOrg  = 0;
583   
584   AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
585   if(aodIH && aodIH->GetEventToMerge())  //Embedding
586     nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
587
588   ResetArrays();
589    
590   for (Int_t i = 0; i < nClusters; i++)
591   {
592     AliVCluster *clus = 0;
593     if(aodIH && aodIH->GetEventToMerge()) //Embedding
594       clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
595     else      
596       clus = fEvent->GetCaloCluster(i);
597     
598     if(!clus) return;
599     
600     if(clus->IsEMCAL())
601     {
602       Int_t label = clus->GetLabel();
603       Int_t label2 = -1 ;
604       //printf("Org cluster E %f, Time  %e, Index = %d, ID %d, MC label %d\n", clus->E(), clus->GetTOF(),i, clus->GetID(),label );
605       //printf("Original list of labels from old cluster : \n");
606       //for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc));
607       
608       if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
609       UShort_t * index    = clus->GetCellsAbsId() ;
610       for(Int_t icell=0; icell < clus->GetNCells(); icell++ )
611       {
612         //printf("\t cell %d, MC label %d\n",index[icell],fEvent->GetEMCALCells()->GetCellMCLabel(index[icell]));
613         fOrgClusterCellId[index[icell]] = i;
614         fCellLabels[index[icell]]       = label;
615         fCellSecondLabels[index[icell]] = label2;
616         fCellTime[index[icell]]         = clus->GetTOF();
617         fCellMatchdEta[index[icell]]    = clus->GetTrackDz();
618         fCellMatchdPhi[index[icell]]    = clus->GetTrackDx();
619       }
620       nClustersOrg++;
621     }
622     // printf("\n");
623   } 
624   
625   // Transform CaloCells into Digits
626   
627   Int_t    idigit =  0;
628   Int_t    id     = -1;
629   Float_t  amp    = -1; 
630   Double_t time   = -1; 
631   
632   AliVCaloCells *cells = fEvent->GetEMCALCells();
633   
634   Int_t bc = InputEvent()->GetBunchCrossNumber();
635
636   for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
637   {
638     // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
639     id = cells->GetCellNumber(icell);
640     Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,cells);
641     
642     // Do not include cells with too low energy, nor exotic cell
643     if( amp  < fRecParam->GetMinECut() ||
644         time > fRecParam->GetTimeMax() ||
645         time < fRecParam->GetTimeMin()    ) accept = kFALSE;
646     
647     // In case of old AOD analysis cell time is -1 s, approximate replacing by time of the cluster the digit belongs.
648     if (fRecalibrateWithClusterTime)
649     { 
650       time = fCellTime[id];
651       //printf("cell %d time org %f - ",id, time*1.e9);
652       fRecoUtils->RecalibrateCellTime(id,bc,time);
653       //printf("recal %f\n",time*1.e9);
654     }
655     
656     //Exotic?
657     if (accept && fRecoUtils->IsExoticCell(id,cells,bc))
658         accept = kFALSE;
659     
660     if( !accept )
661     {
662       if( DebugLevel() > 2 )
663         printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - Remove channel absId %d, index %d of %d, amp %f, time %f\n",
664                id,icell, cells->GetNumberOfCells(), amp, time*1.e9);
665       continue;
666     }
667     
668     Int_t mcLabel = cells->GetMCLabel(icell);
669     //printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - cell %d, mc label %d\n",id,mcLabel);
670
671     //if(fCellLabels[id]!=mcLabel)printf("mcLabel %d - %d\n",mcLabel,fCellLabels[id]);
672     if     ( fSetCellMCLabelFromCluster == 1 ) mcLabel = fCellLabels[id]; // Older aliroot MC productions
673     else if( fSetCellMCLabelFromCluster == 0 && fRemapMCLabelForAODs) RemapMCLabelForAODs(mcLabel);
674     else mcLabel = -1; // found later
675     
676     //printf("\t new label %d\n",mcLabel);
677     
678     // Create the digit, put a fake primary deposited energy to trick the clusterizer
679     // when checking the most likely primary
680     
681     Float_t efrac = cells->GetEFraction(icell);
682     
683     //When checking the MC of digits, give weight to cells with embedded signal
684     if (mcLabel > 0 && efrac < 1.e-6) efrac = 1;
685     
686     //printf("******* Cell %d, id %d, e %f,  fraction %f, MC label %d, used MC label %d\n",icell,id,amp,cells->GetEFraction(icell),cells->GetMCLabel(icell),mcLabel);
687     
688     new((*fDigitsArr)[idigit]) AliEMCALDigit( mcLabel, mcLabel, id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, amp*efrac);
689     // Last parameter should be MC deposited energy, since it is not available, add just the cell amplitude so that
690     // we give more weight to the MC label of the cell with highest energy in the cluster
691         
692     idigit++;
693   }
694   
695   fDigitsArr->Sort();
696   
697   //-------------------------------------------------------------------------------------
698   //Do the clusterization
699   //-------------------------------------------------------------------------------------        
700   
701   fClusterizer->Digits2Clusters("");
702   
703   //-------------------------------------------------------------------------------------
704   //Transform the recpoints into AliVClusters
705   //-------------------------------------------------------------------------------------
706   
707   RecPoints2Clusters();
708   
709   if(!fCaloClusterArr)
710   { 
711     printf("AliAnalysisTaksEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast());
712     return;    
713   }
714   
715   if( DebugLevel() > 0 )
716   {
717     printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - N clusters: before recluster %d, after recluster %d\n",nClustersOrg, fCaloClusterArr->GetEntriesFast());
718     
719     if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast())
720     {
721       printf("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
722              fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
723              fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
724     }
725   }
726 }
727
728 //_____________________________________________________
729 void AliAnalysisTaskEMCALClusterize::ClusterUnfolding()
730 {
731   // Take the event clusters and unfold them
732   
733   AliVCaloCells *cells   = fEvent->GetEMCALCells();
734   Double_t cellAmplitude = 0;
735   Double_t cellTime      = 0;
736   Short_t  cellNumber    = 0;
737   Int_t    cellMCLabel   = 0;
738   Double_t cellEFrac     = 0;
739   Int_t    nClustersOrg  = 0;
740   
741   // Fill the array with the EMCAL clusters, copy them
742   for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++)
743   {
744     AliVCluster *clus = fEvent->GetCaloCluster(i);
745     if(clus->IsEMCAL())
746     {        
747       //recalibrate/remove bad channels/etc if requested
748       if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells()))
749       {
750         continue;
751       } 
752       
753       if(fRecoUtils->IsRecalibrationOn())
754       {
755         //Calibrate cluster
756         fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
757         
758         //CalibrateCells
759         for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
760         {
761           if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime, cellMCLabel, cellEFrac) != kTRUE)
762             break;
763           
764           Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
765           fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta); 
766           fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);       
767           
768           //Do not include bad channels found in analysis?
769           if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() && 
770               fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
771             continue;
772           
773           cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
774           
775         }// cells loop            
776       }// recalibrate
777       
778       //Cast to ESD or AOD, needed to create the cluster array
779       AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
780       AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
781       
782       if     (esdCluster)
783       {
784         fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );   
785       }//ESD
786       else if(aodCluster)
787       {
788         fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );   
789       }//AOD
790       else 
791         Warning("UserExec()"," - Wrong CaloCluster type?");
792       
793       nClustersOrg++;
794     }
795   }
796   
797   //Do the unfolding
798   fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
799   
800   //CLEAN-UP
801   fUnfolder->Clear();
802   
803 }
804
805 //_____________________________________________________
806 void AliAnalysisTaskEMCALClusterize::FillAODCaloCells() 
807 {
808   // Put calo cells in standard branch  
809   AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
810   Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
811   
812   AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
813   aodEMcells.CreateContainer(nEMcell);
814   aodEMcells.SetType(AliVCaloCells::kEMCALCell);
815   Double_t calibFactor = 1.;   
816   for (Int_t iCell = 0; iCell < nEMcell; iCell++) 
817   { 
818     Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
819     fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta); 
820     fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);     
821     
822     if(fRecoUtils->IsRecalibrationOn())
823     { 
824       calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
825     }
826     
827     if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
828     { //Channel is not declared as bad
829       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor,
830                          eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell));
831     }
832     else 
833     {
834       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
835     }
836   }
837   aodEMcells.Sort();
838   
839 }
840
841 //__________________________________________________
842 void AliAnalysisTaskEMCALClusterize::FillAODHeader() 
843 {
844   //Put event header information in standard AOD branch
845   
846   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent);
847   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent);
848   
849   Double_t pos[3]   ;
850   Double_t covVtx[6];
851   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
852   
853   AliAODHeader* header = AODEvent()->GetHeader();
854   header->SetRunNumber(fEvent->GetRunNumber());
855   
856   if(esdevent)
857   {
858     TTree* tree = fInputHandler->GetTree();
859     if (tree) 
860     {
861       TFile* file = tree->GetCurrentFile();
862       if (file) header->SetESDFileName(file->GetName());
863     }
864   }
865   else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName());
866   
867   header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
868   header->SetOrbitNumber(fEvent->GetOrbitNumber());
869   header->SetPeriodNumber(fEvent->GetPeriodNumber());
870   header->SetEventType(fEvent->GetEventType());
871   
872   //Centrality
873   if(fEvent->GetCentrality())
874     header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
875   else
876     header->SetCentrality(0);
877   
878   //Trigger  
879   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
880
881   header->SetFiredTriggerClasses(fEvent->GetFiredTriggerClasses());
882   
883   header->SetTriggerMask(fEvent->GetTriggerMask()); 
884   header->SetTriggerCluster(fEvent->GetTriggerCluster());
885   
886   if      (esdevent)
887   {
888     header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());    
889     header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());    
890     header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());    
891   }
892   else if (aodevent)
893   {
894     header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());    
895     header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());    
896     header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());    
897   }
898   
899   header->SetMagneticField(fEvent->GetMagneticField());
900   //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.); 
901   
902   header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
903   header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
904   header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
905   header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
906   header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
907   
908   Float_t diamxy[2]={fEvent->GetDiamondX(),fEvent->GetDiamondY()};
909   Float_t diamcov[3];
910   fEvent->GetDiamondCovXY(diamcov);
911   header->SetDiamond(diamxy,diamcov);
912   if      (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
913   else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
914   
915   //
916   Int_t nVertices = 1 ;/* = prim. vtx*/;
917   Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
918   
919   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
920   
921   // Access to the AOD container of vertices
922   TClonesArray &vertices = *(AODEvent()->GetVertices());
923   Int_t jVertices=0;
924   
925   // Add primary vertex. The primary tracks will be defined
926   // after the loops on the composite objects (V0, cascades, kinks)
927   fEvent->GetPrimaryVertex()->GetXYZ(pos);
928   Float_t chi = 0;
929   if      (esdevent)
930   {
931     esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
932     chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
933   }
934   else if (aodevent)
935   {
936     aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
937     chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
938   }
939   
940   AliAODVertex * primary = new(vertices[jVertices++])
941   AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
942   primary->SetName(fEvent->GetPrimaryVertex()->GetName());
943   primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
944   
945 }
946
947 //___________________________________________________________
948 void AliAnalysisTaskEMCALClusterize::FillCaloClusterInEvent()
949 {
950   // Get the CaloClusters array, do some final calculations 
951   // and put the clusters in the output or input event 
952   // as a separate branch 
953   
954   //First recalculate track-matching for the new clusters
955   if(fDoTrackMatching) 
956   {
957     fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom);
958   }
959   //Put the new clusters in the AOD list
960   
961   Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntriesFast();
962
963   for(Int_t i = 0; i < kNumberOfCaloClusters; i++)
964   {
965     AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
966     
967     newCluster->SetID(i);
968     
969     // Correct cluster energy non linearity    
970     newCluster->SetE(fRecoUtils->CorrectClusterEnergyLinearity(newCluster));
971     
972     //Add matched track
973     if(fDoTrackMatching)
974     {
975       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
976       if(trackIndex >= 0)
977       {
978         newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex));
979         if(DebugLevel() > 1) 
980           printf("AliAnalysisTaksEMCALClusterize::UserExec() - Matched Track index %d to new cluster %d \n",trackIndex,i);
981       }
982       
983       Float_t dR = 999., dZ = 999.;
984       fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dR,dZ);
985       newCluster->SetTrackDistance(dR,dZ);
986       
987     }
988     else 
989     {// Assign previously assigned matched track in reco, very very rough
990       Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster
991       newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]);
992     }
993     
994     //printf("New cluster E %f, Time  %e, Id = ", newCluster->E(), newCluster->GetTOF() );
995     //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] );
996     //printf("\n");
997     
998     // Calculate distance to bad channel for new cluster. Make sure you give the list of bad channels.
999     fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fEvent->GetEMCALCells(), newCluster);
1000     
1001     new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
1002     
1003     if(DebugLevel() > 1 )
1004       printf("AliAnalysisTaksEMCALClusterize::UserExec() - New cluster %d of %d, energy %f, mc label %d \n",newCluster->GetID(), kNumberOfCaloClusters, newCluster->E(), newCluster->GetLabel());
1005     
1006   } // cluster loop
1007   
1008   fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
1009   
1010   
1011 }
1012
1013 //_______________________________________________
1014 TString AliAnalysisTaskEMCALClusterize::GetPass()
1015 {
1016   // Get passx from filename.
1017   
1018   if (!AliAnalysisManager::GetAnalysisManager()->GetTree()) 
1019   {
1020     AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Pointer to tree = 0, returning null\n");
1021     return TString("");
1022   }
1023   
1024   if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()) 
1025   {
1026     AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Null pointer input file, returning null\n");
1027     return TString("");
1028   }
1029   
1030   TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1031   if      (pass.Contains("ass1")) return TString("pass1");
1032   else if (pass.Contains("ass2")) return TString("pass2");
1033   else if (pass.Contains("ass3")) return TString("pass3");
1034   else if (pass.Contains("ass4")) return TString("pass4");
1035   else if (pass.Contains("ass5")) return TString("pass5");
1036   else if (pass.Contains("calo") || pass.Contains("high_lumi"))
1037   {
1038     printf("AliAnalysisTaskEMCALClusterize::GetPass() - Path contains <calo> or <high-lumi>, set as <pass1>\n");
1039     return TString("pass1");
1040   }
1041   
1042   // No condition fullfilled, give a default value
1043   printf("AliAnalysisTaskEMCALClusterize::GetPass() - Pass number string not found \n");
1044   return TString("");            
1045   
1046 }
1047
1048 //_________________________________________
1049 void AliAnalysisTaskEMCALClusterize::Init()
1050 {
1051   //Init analysis with configuration macro if available
1052   
1053   fOADBSet           = kFALSE;
1054   if(fOADBFilePath == "") fOADBFilePath = "$ALICE_ROOT/OADB/EMCAL" ;          
1055   
1056   fBranchNames       = "ESD:AliESDHeader.,EMCALCells.";
1057   
1058   if(!fRecParam)     fRecParam  = new AliEMCALRecParam;
1059   if(!fRecoUtils)    fRecoUtils = new AliEMCALRecoUtils();  
1060   
1061   if(fMaxEvent          <= 0) fMaxEvent          = 1000000000;
1062   if(fSelectCellMinE    <= 0) fSelectCellMinE    = 0.005;     
1063   if(fSelectCellMinFrac <= 0) fSelectCellMinFrac = 0.001;
1064   
1065   //Centrality
1066   if(fCentralityClass  == "") fCentralityClass  = "V0M";
1067   
1068   if (fOCDBpath            == "") fOCDBpath            = "raw://" ;
1069   if (fOutputAODBranchName == "") fOutputAODBranchName = "newEMCALClusters" ;
1070   
1071   if(gROOT->LoadMacro(fConfigName) >=0)
1072   {
1073     printf("AliAnalysisTaksEMCALClusterize::Init() - Configure analysis with %s\n",fConfigName.Data());
1074     AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
1075     fGeomName         = clus->fGeomName; 
1076     fLoadGeomMatrices = clus->fLoadGeomMatrices;
1077     fOCDBpath         = clus->fOCDBpath;   
1078     fAccessOCDB       = clus->fAccessOCDB;
1079     fRecParam         = clus->fRecParam;
1080     fJustUnfold       = clus->fJustUnfold;
1081     fFillAODFile      = clus->fFillAODFile;
1082     fRecoUtils        = clus->fRecoUtils; 
1083     fConfigName       = clus->fConfigName;
1084     fMaxEvent         = clus->fMaxEvent;
1085     fDoTrackMatching  = clus->fDoTrackMatching;
1086     fOutputAODBranchName = clus->fOutputAODBranchName;
1087     for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
1088     fCentralityClass  = clus->fCentralityClass;
1089     fCentralityBin[0] = clus->fCentralityBin[0];
1090     fCentralityBin[1] = clus->fCentralityBin[1];
1091   }
1092   
1093   // Init geometry, I do not like much to do it like this ...
1094   if(fImportGeometryFromFile && !gGeoManager) 
1095   {
1096     if (fImportGeometryFilePath == "") fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2011.root" ; // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1097     printf("AliAnalysisTaskEMCALClusterize::Init() - Import %s\n",fImportGeometryFilePath.Data());
1098     TGeoManager::Import(fImportGeometryFilePath) ; 
1099   }
1100   
1101 }  
1102
1103 //_______________________________________________________
1104 void AliAnalysisTaskEMCALClusterize::InitClusterization()
1105 {
1106   //Select clusterization/unfolding algorithm and set all the needed parameters
1107   
1108   if (fJustUnfold)
1109   {
1110     // init the unfolding afterburner 
1111     delete fUnfolder;
1112     fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
1113     return;
1114   }
1115   
1116   //First init the clusterizer
1117   delete fClusterizer;
1118   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1119     fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
1120   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2) 
1121     fClusterizer = new AliEMCALClusterizerv2(fGeom, fCalibData, fPedestalData);
1122   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
1123   { 
1124     fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
1125     fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
1126     fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
1127   } 
1128   else 
1129   {
1130     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1131   }
1132   
1133   //Now set the parameters
1134   fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
1135   fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
1136   fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );    
1137   fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
1138   fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
1139   fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
1140   fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
1141   fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
1142   fClusterizer->SetInputCalibrated       ( kTRUE                               );
1143   fClusterizer->SetJustClusters          ( kTRUE                               );  
1144   
1145   // Initialize the cluster rec points and digits arrays and get them.
1146   fClusterizer->SetOutput(0);
1147   fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1148   fDigitsArr  = fClusterizer->GetDigits();
1149   
1150   //In case of unfolding after clusterization is requested, set the corresponding parameters
1151   if(fRecParam->GetUnfold())
1152   {
1153     Int_t i=0;
1154     for (i = 0; i < 8; i++) 
1155     {
1156       fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
1157     }//end of loop over parameters
1158     
1159     for (i = 0; i < 3; i++) 
1160     {
1161       fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
1162       fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
1163     }//end of loop over parameters
1164     
1165     fClusterizer->InitClusterUnfolding();
1166     
1167   }// to unfold
1168 }
1169
1170 //_________________________________________________
1171 void AliAnalysisTaskEMCALClusterize::InitGeometry()
1172 {
1173   // Init geometry and set the geometry matrix, for the first event, skip the rest
1174   // Also set once the run dependent calibrations
1175   
1176   if(fGeomMatrixSet) return;
1177   
1178   Int_t runnumber = InputEvent()->GetRunNumber() ;
1179   if (!fGeom)
1180   {
1181     if(fGeomName=="")
1182     {
1183       if     (runnumber < 140000) fGeomName = "EMCAL_FIRSTYEARV1";
1184       else if(runnumber < 171000) fGeomName = "EMCAL_COMPLETEV1";
1185       else                        fGeomName = "EMCAL_COMPLETE12SMV1";  
1186       printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Set EMCAL geometry name to <%s> for run %d\n",fGeomName.Data(),runnumber);
1187     }
1188     
1189                 fGeom = AliEMCALGeometry::GetInstance(fGeomName);
1190     
1191                 if(fDebug > 0)
1192     {
1193                         printf("AliAnalysisTaskEMCALClusterize::InitGeometry(run=%d)",runnumber);
1194                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1195                         printf("\n");
1196                 }
1197         } // geometry pointer did not exist before
1198   
1199   if(fLoadGeomMatrices)
1200   {
1201     printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Load user defined EMCAL geometry matrices\n");
1202     
1203     // OADB if available
1204     AliOADBContainer emcGeoMat("AliEMCALgeo");
1205     emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
1206     TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
1207     
1208     
1209     for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1210     {
1211       
1212       if (!fGeomMatrix[mod]) // Get it from OADB
1213       {
1214         if(fDebug > 1 ) 
1215           printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - EMCAL matrices SM %d, %p\n",
1216                  mod,((TGeoHMatrix*) matEMCAL->At(mod)));
1217         //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
1218         
1219         fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
1220       }        
1221       
1222       if(fGeomMatrix[mod])
1223       {
1224         if(DebugLevel() > 1) 
1225           fGeomMatrix[mod]->Print();
1226         
1227         fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;  
1228       }
1229       
1230       fGeomMatrixSet=kTRUE;
1231       
1232     }//SM loop
1233   }//Load matrices
1234   else if(!gGeoManager)
1235   {
1236     printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - Get geo matrices from data");
1237     //Still not implemented in AOD, just a workaround to be able to work at least with ESDs     
1238     if(!strcmp(fEvent->GetName(),"AliAODEvent")) 
1239     {
1240       if(DebugLevel() > 1) 
1241         Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
1242     }//AOD
1243     else 
1244     {   
1245       if(DebugLevel() > 1) 
1246         printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
1247       
1248       AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ;
1249       
1250       if(!esd)
1251       {
1252         Error("InitGeometry"," - This event does not contain ESDs?");
1253         if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
1254         return;
1255       }
1256       
1257       for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1258       {
1259         if(DebugLevel() > 1) 
1260           esd->GetEMCALMatrix(mod)->Print();
1261         
1262         if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
1263         
1264       } 
1265       
1266       fGeomMatrixSet=kTRUE;
1267       
1268     }//ESD
1269   }//Load matrices from Data 
1270   
1271 }
1272
1273 //____________________________________________________
1274 Bool_t AliAnalysisTaskEMCALClusterize::IsExoticEvent()
1275 {
1276   
1277   // Check if event is exotic, get an exotic cell and compare with triggered patch
1278   // If there is a match remove event ... to be completed, filled with something provisional
1279   
1280   if(!fRemoveExoticEvents) return kFALSE;
1281   
1282   // Loop on cells
1283   AliVCaloCells * cells = fEvent->GetEMCALCells();
1284   Float_t totCellE = 0;
1285   Int_t bc = InputEvent()->GetBunchCrossNumber();
1286   for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1287   {
1288     
1289     Float_t  ecell = 0 ;
1290     Double_t tcell = 0 ;
1291     
1292     Int_t absID   = cells->GetCellNumber(icell);
1293     Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,cells);
1294     if(accept && !fRecoUtils->IsExoticCell(absID,cells,bc)) totCellE += ecell;
1295   }
1296   
1297   //  TString triggerclasses = event->GetFiredTriggerClasses();
1298   //    printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster  - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
1299   //    if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1300   //    return;
1301   //  
1302   
1303   //printf("TotE cell %f\n",totCellE);
1304   if(totCellE < 1) return kTRUE;
1305   
1306   return kFALSE;
1307   
1308
1309
1310 //________________________________________________________________
1311 Bool_t AliAnalysisTaskEMCALClusterize::IsLEDEvent(const Int_t run)
1312 {
1313   //Check if event is LED
1314   
1315   if(!fRemoveLEDEvents) return kFALSE;
1316   
1317   //check events of LHC11a period
1318   if(run < 146858 || run > 146860) return kFALSE ;
1319   
1320   // Count number of cells with energy larger than 0.1 in SM3, cut on this number
1321   Int_t ncellsSM3 = 0;
1322   AliVCaloCells * cells = fEvent->GetEMCALCells();
1323   for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1324   {
1325     if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
1326   }
1327   
1328   TString triggerclasses = fEvent->GetFiredTriggerClasses();
1329   
1330   Int_t ncellcut = 21;
1331   if(triggerclasses.Contains("EMC")) ncellcut = 35;
1332   
1333   if( ncellsSM3 >= ncellcut)
1334   {
1335     printf("AliAnalysisTaksEMCALClusterize::IsLEDEvent() - reject event %d with ncells in SM3 %d\n",(Int_t)Entry(),ncellsSM3);
1336     if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1337     return kTRUE;
1338   }
1339   
1340   return kFALSE;
1341   
1342
1343
1344 //_______________________________________________________
1345 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters()
1346 {
1347   // Restore clusters from recPoints
1348   // Cluster energy, global position, cells and their amplitude 
1349   // fractions are restored
1350   
1351   Int_t j = 0;
1352   for(Int_t i = 0; i < fClusterArr->GetEntriesFast(); i++)
1353   {
1354     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fClusterArr->At(i);
1355     
1356     const Int_t ncells = recPoint->GetMultiplicity();
1357     Int_t ncellsTrue = 0;
1358     
1359     if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
1360     
1361     // cells and their amplitude fractions
1362     UShort_t   absIds[ncells];  
1363     Double32_t ratios[ncells];
1364     
1365     //For later check embedding
1366     AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1367     
1368     Float_t clusterE = 0; 
1369     for (Int_t c = 0; c < ncells; c++) 
1370     {
1371       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
1372       
1373       absIds[ncellsTrue] = digit->GetId();
1374       ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
1375             
1376       // In case of unfolding, remove digits with unfolded energy too low      
1377       if(fSelectCell)
1378       {
1379         if     (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)  
1380         {
1381           if(DebugLevel() > 1)
1382           {
1383             printf("AliAnalysisTaksEMCALClusterize::RecPoints2Clusters() - Too small energy in cell of cluster: cluster cell %f, digit %f\n",
1384                    recPoint->GetEnergiesList()[c],digit->GetAmplitude());
1385           }
1386           
1387           continue;
1388           
1389         } // if cuts
1390       }// Select cells
1391       
1392       //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
1393       clusterE  +=recPoint->GetEnergiesList()[c];
1394       
1395       // In case of embedding, fill ratio with amount of signal, 
1396       if (aodIH && aodIH->GetMergeEvents()) 
1397       {
1398         //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
1399         AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
1400         AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
1401         
1402         Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1403         //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1404         Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1405         //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
1406         
1407         if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
1408         //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
1409         
1410       }//Embedding
1411       
1412       ncellsTrue++;
1413       
1414     }// cluster cell loop
1415     
1416     if (ncellsTrue < 1) 
1417     {
1418       if (DebugLevel() > 1) 
1419         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Skipping cluster with no cells avobe threshold E = %f, ncells %d\n",
1420                recPoint->GetEnergy(), ncells);
1421       continue;
1422     }
1423     
1424     //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
1425     
1426     if(clusterE <  fRecParam->GetClusteringThreshold()) 
1427     {
1428       if (DebugLevel()>1)
1429         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Remove cluster with energy below seed threshold %f\n",clusterE);
1430       continue;
1431     }
1432     
1433     TVector3 gpos;
1434     Float_t g[3];
1435     
1436     // calculate new cluster position
1437     
1438     recPoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
1439     recPoint->GetGlobalPosition(gpos);
1440     gpos.GetXYZ(g);
1441     
1442     // create a new cluster
1443     
1444     (*fCaloClusterArr)[j] = new AliAODCaloCluster() ;
1445     AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( fCaloClusterArr->At(j) ) ;
1446     j++;
1447     clus->SetType(AliVCluster::kEMCALClusterv1);
1448     clus->SetE(clusterE);
1449     clus->SetPosition(g);
1450     clus->SetNCells(ncellsTrue);
1451     clus->SetCellsAbsId(absIds);
1452     clus->SetCellsAmplitudeFraction(ratios);
1453     clus->SetChi2(-1); //not yet implemented
1454     clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
1455     clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
1456     clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower()); 
1457     
1458     if(ncells == ncellsTrue)
1459     {
1460       Float_t elipAxis[2];
1461       recPoint->GetElipsAxis(elipAxis);
1462       clus->SetM02(elipAxis[0]*elipAxis[0]) ;
1463       clus->SetM20(elipAxis[1]*elipAxis[1]) ;
1464       clus->SetDispersion(recPoint->GetDispersion());
1465     }
1466     else if(fSelectCell)
1467     {
1468       // In case some cells rejected, in unfolding case, recalculate
1469       // shower shape parameters and position
1470       if(DebugLevel() > 1) 
1471         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape\n",ncells,ncellsTrue);
1472       
1473       AliVCaloCells* cells = 0x0; 
1474       if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent()  ->GetEMCALCells();
1475       else                                  cells = InputEvent()->GetEMCALCells();
1476       
1477       fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus);
1478       fRecoUtils->RecalculateClusterPID(clus);
1479       fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus); 
1480       
1481     }
1482     
1483     // MC
1484
1485     if     ( fSetCellMCLabelFromCluster == 1 ) SetClustersMCLabelFrom2SelectedLabels(recPoint,clus) ;
1486     else if( fSetCellMCLabelFromCluster == 2 ) SetClustersMCLabelFromOriginalClusters(clus) ;
1487     else
1488     {
1489       // Normal case, trust what the clusterizer has found
1490       Int_t  parentMult = 0;
1491       Int_t *parentList = recPoint->GetParents(parentMult);
1492       clus->SetLabel(parentList, parentMult);
1493 //      printf("Label list : ");
1494 //      for(Int_t ilabel = 0; ilabel < parentMult; ilabel++ ) printf(" %d ",parentList[ilabel]);
1495 //      printf("\n");
1496     }
1497     
1498   } // recPoints loop
1499   
1500 }
1501
1502 //___________________________________________________________________________
1503 void AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs(Int_t & label)
1504 {
1505   // MC label for Cells not remapped after ESD filtering, do it here.
1506
1507   if(label < 0) return ;
1508   
1509   AliAODEvent  * evt = dynamic_cast<AliAODEvent*> (fEvent) ;
1510   if(!evt) return ;
1511   
1512   TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
1513   if(!arr) return ;
1514     
1515   if(label < arr->GetEntriesFast())
1516   {
1517     AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
1518     if(!particle) return ;
1519         
1520     if(label == particle->Label()) return ; // label already OK
1521     //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label  %d - AOD stack %d \n",label, particle->Label());
1522   }
1523   //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label  %d > AOD labels %d \n",label, arr->GetEntriesFast());
1524   
1525   // loop on the particles list and check if there is one with the same label
1526   for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
1527   {
1528     AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
1529     if(!particle) continue ;
1530
1531     if(label == particle->Label())
1532     {
1533       label = ind;
1534       //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index  %d \n",label);
1535       return;
1536     }
1537   }
1538   
1539   label = -1;
1540   
1541   //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label not found set to -1 \n");
1542  
1543 }
1544
1545 //________________________________________________
1546 void AliAnalysisTaskEMCALClusterize::ResetArrays()
1547 {
1548   // Reset arrays containing information for all possible cells
1549   for(Int_t j = 0; j < 12672; j++)
1550   {
1551     fOrgClusterCellId[j] = -1;
1552     fCellLabels[j]       = -1;
1553     fCellSecondLabels[j] = -1;
1554     fCellTime[j]         =  0.;
1555     fCellMatchdEta[j]    = -999;
1556     fCellMatchdPhi[j]    = -999;
1557   }
1558 }
1559
1560 //_____________________________________________________________________________________________________
1561 void AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFrom2SelectedLabels(AliEMCALRecPoint  * recPoint,
1562                                                                            AliAODCaloCluster * clus)
1563 {
1564   // Set the cluster MC label, the digizer was filled with most likely MC label for all cells in original cluster
1565   // Now check the second most likely MC label and add it to the new cluster
1566   
1567   Int_t  parentMult = 0;
1568   Int_t *parentList = recPoint->GetParents(parentMult);
1569   clus->SetLabel(parentList, parentMult);
1570   
1571   //Write the second major contributor to each MC cluster.
1572   Int_t iNewLabel ;
1573   for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1574   {
1575     
1576     Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1577     if(idCell>=0)
1578     {
1579       iNewLabel = 1 ; //iNewLabel makes sure we  don't write twice the same label.
1580       for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
1581       {
1582         if ( fCellSecondLabels[idCell] == -1 )  iNewLabel = 0;  // -1 is never a good second label.
1583           if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) )  iNewLabel = 0;
1584             }
1585       if (iNewLabel == 1)
1586       {
1587         Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
1588         for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ )
1589         {
1590           newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
1591         }
1592         
1593         newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
1594         clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
1595         delete [] newLabelArray;
1596       }
1597     }//positive cell number
1598   }
1599 }
1600
1601 //___________________________________________________________________________________________________
1602 void AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFromOriginalClusters(AliAODCaloCluster * clus)
1603 {
1604   // Get the original clusters that contribute to the new cluster, assign the labels of such clusters
1605   // to the new cluster.
1606   // Only approximatedly valid  when output are V1 clusters, handle with care
1607     
1608   TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ...
1609   clArray.Reset();
1610   Int_t nClu = 0;
1611   Int_t nLabTotOrg = 0;
1612   Float_t emax = -1;
1613   Int_t idMax = -1;
1614   
1615   AliVEvent * event = fEvent;
1616   
1617   //In case of embedding the MC signal is in the added event
1618   AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1619   if(aodIH && aodIH->GetEventToMerge())  //Embedding
1620     event = aodIH->GetEventToMerge(); //Get clusters directly from embedded signal
1621
1622   
1623   //Find the clusters that originally had the cells
1624   for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1625   {
1626     Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1627     
1628     if(idCell>=0)
1629     {
1630       Int_t idCluster = fOrgClusterCellId[idCell];
1631       
1632       Bool_t set = kTRUE;
1633       for(Int_t icl =0; icl < nClu; icl++)
1634       {
1635         if(((Int_t)clArray.GetAt(icl))==-1 ) continue;
1636         if( idCluster == ((Int_t)clArray.GetAt(icl)) ) set = kFALSE;
1637         //   printf("\t \t icell %d  Cluster in array %d, IdCluster %d, in array %d, set %d\n",
1638         //          iLoopCell,                 icl,  idCluster,((Int_t)clArray.GetAt(icl)),set);
1639       }
1640       if( set && idCluster >= 0)
1641       {
1642         clArray.SetAt(idCluster,nClu++);
1643         //printf("******** idCluster %d \n",idCluster);
1644         nLabTotOrg+=(event->GetCaloCluster(idCluster))->GetNLabels();
1645
1646         //printf("Cluster in array %d, IdCluster %d\n",nClu-1,  idCluster);
1647
1648         //Search highest E cluster
1649         AliVCluster * clOrg = event->GetCaloCluster(idCluster);
1650         //printf("\t E %f\n",clOrg->E());
1651         if(emax < clOrg->E())
1652         {
1653           emax  = clOrg->E();
1654           idMax = idCluster;
1655         }        
1656       }
1657     }
1658   }// cell loop
1659
1660   if(nClu==0 || nLabTotOrg == 0)
1661   {
1662     //if(clus->E() > 0.25) printf("AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFromOriginalClusters() - Check: N org clusters %d, n tot labels %d, cluster E %f,  n cells %d\n",nClu,nLabTotOrg,clus->E(), clus->GetNCells());
1663     //for(Int_t icell = 0; icell < clus->GetNCells(); icell++) printf("\t cell %d",clus->GetCellsAbsId()[icell]);
1664     //printf("\n");
1665   }
1666   
1667   // Put the first in the list the cluster with highest energy
1668   if(idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position
1669   {
1670     Int_t maxIndex = -1;
1671     Int_t firstCluster = ((Int_t)clArray.GetAt(0));
1672     for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
1673     {
1674       if(idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster;
1675     }
1676     
1677     if(firstCluster >=0 && idMax >=0)
1678     {
1679       clArray.SetAt(idMax,0);
1680       clArray.SetAt(firstCluster,maxIndex);
1681     }
1682   }
1683   
1684   // Get the labels list in the original clusters, assign all to the new cluster
1685   TArrayI clMCArray(nLabTotOrg) ;
1686   clMCArray.Reset();
1687   
1688   Int_t nLabTot = 0;
1689   for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
1690   {
1691     Int_t idCluster = (Int_t) clArray.GetAt(iLoopCluster);
1692     //printf("New Cluster in Array %d,  idCluster %d \n",iLoopCluster,idCluster);
1693     AliVCluster * clOrg = event->GetCaloCluster(idCluster);
1694     Int_t nLab = clOrg->GetNLabels();
1695     
1696     for ( Int_t iLab = 0 ; iLab < nLab ; iLab++ )
1697     {
1698       Int_t lab = clOrg->GetLabelAt(iLab) ;
1699       if(lab>=0)
1700       {
1701         Bool_t set = kTRUE;
1702         //printf("\t \t Set Label %d \n", lab);
1703         for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++)
1704         {
1705           if( lab == ((Int_t)clMCArray.GetAt(iLabTot)) ) set = kFALSE;
1706           //printf("iLoopCluster %d, Label ID in Org Cluster %d,label %d  Label ID in array %d, label in array %d, set %d\n",
1707           //       iLoopCluster,                           iLab,     lab,           iLabTot,  ((Int_t)clMCArray.GetAt(iLabTot)),set);
1708         }
1709         if( set ) clMCArray.SetAt(lab,nLabTot++);
1710       }
1711     }
1712   }// cluster loop
1713   
1714   // Set the final list of labels
1715   
1716   clus->SetLabel(clMCArray.GetArray(), nLabTot);
1717   
1718 //  printf("Final list of labels for new cluster : \n");
1719 //  for(Int_t ice = 0; ice < clus->GetNCells() ; ice++)
1720 //  {
1721 //    printf("\t Cell %d ",clus->GetCellsAbsId()[ice]);
1722 //    Int_t label = InputEvent()->GetEMCALCells()->GetCellMCLabel(clus->GetCellsAbsId()[ice]);
1723 //    printf(" org %d ",label);
1724 //    RemapMCLabelForAODs(label);
1725 //    printf(" new %d \n",label);
1726 //  }
1727 //  for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc));
1728 }
1729
1730
1731 //____________________________________________________________
1732 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
1733 {
1734   // Init geometry, create list of output clusters
1735   
1736
1737   fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1738
1739   if(fOutputAODBranchName.Length()==0)
1740   {
1741     fOutputAODBranchName = "newEMCALClustersArray";
1742     printf("Cluster branch name not set, set it to newEMCALClustersArray \n");
1743   }
1744   
1745   fOutputAODBranch->SetName(fOutputAODBranchName);
1746   
1747   if( fFillAODFile )
1748   {
1749     //fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1750
1751     
1752     //fOutputAODBranch->SetOwner(kFALSE);
1753     
1754     AddAODBranch("TClonesArray", &fOutputAODBranch);
1755   }
1756   
1757 }
1758
1759 //_______________________________________________________
1760 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *) 
1761 {
1762   // Do clusterization event by event, execute different steps
1763   // 1) Do some checks on the kind of events (ESD, AOD) or if some filtering is needed, initializations
1764   //    load things and clear arrays
1765   // 2) Clusterize a) just unfolding existing clusters (fJustUnfold)
1766   //               b) recluster cells
1767   //                   +  convert cells into digits (calibrating them)
1768   //                   +  recluster digits into recPoints with the chosen clusterizer (V1, V1+Unfold,V2, NxN) 
1769   //                      with methods in AliEMCALClusterizer
1770   //                   +  transform recPoints into CaloClusters
1771   // 3) Do some final calculations in the found clusters (track-matching) and put them in an AOD branch
1772   
1773   //-------
1774   // Step 1
1775   
1776
1777   //Remove the contents of AOD branch output list set in the previous event 
1778   fOutputAODBranch->Clear("C");
1779
1780   LoadBranches();
1781   
1782   //Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
1783   //If we need a centrality bin, we select only those events in the corresponding bin.
1784   if( GetCentrality() && fCentralityBin[0] >= 0 && fCentralityBin[1] >= 0 )
1785   {
1786     Float_t cen = GetEventCentrality();
1787     if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return ; //reject events out of bin.
1788   }  
1789     
1790   // intermediate array with new clusters : init the array only once or clear from previous event
1791   if(!fCaloClusterArr) fCaloClusterArr    = new TObjArray(10000);
1792   else                 fCaloClusterArr->Delete();//Clear("C"); it leaks?
1793
1794   InitGeometry(); // only once, must be done before OADB, geo OADB accessed here
1795   
1796   //Get the event, do some checks and settings
1797   CheckAndGetEvent() ;
1798   
1799   if (!fEvent) 
1800   {
1801     if(DebugLevel() > 0 ) printf("AliAnalysisTaksEMCALClusterize::UserExec() - Skip Event %d", (Int_t) Entry());
1802     return ;
1803   }
1804
1805   //Init pointers, geometry, clusterizer, ocdb, aodb
1806   
1807   
1808   if(fAccessOCDB) AccessOCDB();
1809   if(fAccessOADB) AccessOADB(); // only once
1810   
1811   InitClusterization();
1812   
1813   //-------
1814   // Step 2
1815   
1816   // Make clusters
1817   if (fJustUnfold) ClusterUnfolding();
1818   else             ClusterizeCells() ;
1819   
1820   //-------
1821   // Step 3
1822   
1823   FillCaloClusterInEvent();
1824   
1825 }      
1826