Transition PWG4 --> PWGGA
[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 /* $Id: $ */
16
17 //_________________________________________________________________________
18 // This analysis provides a new list of clusters to be used in other analysis
19 //
20 // Author: Gustavo Conesa Balbastre,
21 //         Adapted from analysis class from Deepa Thomas
22 //
23 //
24 //_________________________________________________________________________
25
26 // --- Root ---
27 #include "TString.h"
28 #include "TRefArray.h"
29 #include "TClonesArray.h"
30 #include "TTree.h"
31 #include "TGeoManager.h"
32 #include "TROOT.h"
33 #include "TInterpreter.h"
34 #include "TFile.h"
35
36 // --- AliRoot Analysis Steering
37 #include "AliAnalysisTask.h"
38 #include "AliAnalysisManager.h"
39 #include "AliESDEvent.h"
40 #include "AliGeomManager.h"
41 #include "AliVCaloCells.h"
42 #include "AliAODCaloCluster.h"
43 #include "AliCDBManager.h"
44 #include "AliCDBStorage.h"
45 #include "AliCDBEntry.h"
46 #include "AliLog.h"
47 #include "AliVEventHandler.h"
48 #include "AliAODInputHandler.h"
49
50 // --- EMCAL
51 #include "AliEMCALRecParam.h"
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 #include "AliEMCALRecoUtils.h"
62
63 #include "AliAnalysisTaskEMCALClusterize.h"
64
65 ClassImp(AliAnalysisTaskEMCALClusterize)
66
67 //______________________________________________________________________________
68 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name) 
69 : AliAnalysisTaskSE(name)
70 , fEvent(0)
71 , fGeom(0),               fGeomName("EMCAL_COMPLETEV1") 
72 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
73 , fCalibData(0),          fPedestalData(0)
74 , fOCDBpath("raw://"),    fAccessOCDB(kFALSE)
75 , fDigitsArr(0),          fClusterArr(0),             fCaloClusterArr(0)
76 , fRecParam(0),           fClusterizer(0)
77 , fUnfolder(0),           fJustUnfold(kFALSE) 
78 , fOutputAODBranch(0),    fOutputAODBranchName("newEMCALClusters")
79 , fFillAODFile(kTRUE),    fFillAODHeader(0)
80 , fFillAODCaloCells(0),   fRun(-1)
81 , fRecoUtils(0),          fConfigName("")
82 , fCellLabels(),          fCellSecondLabels(),        fCellTime()
83 , fCellMatchdEta(),       fCellMatchdPhi()
84 , fMaxEvent(1000000000),  fDoTrackMatching(kFALSE)
85 , fSelectCell(kFALSE),    fSelectCellMinE(0.005),     fSelectCellMinFrac(0.001)
86 , fRemoveLEDEvents(kFALSE),        fRemoveExoticEvents(kFALSE)
87 , fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("")
88 {
89   //ctor
90   for(Int_t i = 0; i < 10;    i++)  fGeomMatrix[i] =  0;
91   for(Int_t j = 0; j < 24*48*11; j++)  {
92     fCellLabels[j]       = -1;
93     fCellSecondLabels[j] = -1;
94     fCellTime[j]         =  0.;    
95     fCellMatchdEta[j]    = -999;
96     fCellMatchdPhi[j]    = -999;
97   }  
98   
99   fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
100   fClusterArr      = new TObjArray(100);
101   fCaloClusterArr  = new TObjArray(1000);
102   fRecParam        = new AliEMCALRecParam;
103   fBranchNames     = "ESD:AliESDHeader.,EMCALCells.";
104   fRecoUtils       = new AliEMCALRecoUtils();
105   
106 }
107
108 //______________________________________________________________
109 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize() 
110 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
111 , fEvent(0)
112 , fGeom(0),                 fGeomName("EMCAL_COMPLETEV1") 
113 , fGeomMatrixSet(kFALSE),   fLoadGeomMatrices(kFALSE)
114 , fCalibData(0),            fPedestalData(0)
115 , fOCDBpath("raw://"),      fAccessOCDB(kFALSE)
116 , fDigitsArr(0),            fClusterArr(0),             fCaloClusterArr(0)
117 , fRecParam(0),             fClusterizer(0)
118 , fUnfolder(0),             fJustUnfold(kFALSE) 
119 , fOutputAODBranch(0),      fOutputAODBranchName("newEMCALClusters")
120 , fFillAODFile(kTRUE),      fFillAODHeader(0)
121 , fFillAODCaloCells(0),     fRun(-1)
122 , fRecoUtils(0),            fConfigName("")
123 , fCellLabels(),            fCellSecondLabels(),        fCellTime()
124 , fCellMatchdEta(),         fCellMatchdPhi()
125 , fMaxEvent(1000000000),    fDoTrackMatching(kFALSE)
126 , fSelectCell(kFALSE),      fSelectCellMinE(0.005),     fSelectCellMinFrac(0.001)
127 , fRemoveLEDEvents(kFALSE), fRemoveExoticEvents(kFALSE)
128 , fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("")
129 {
130   // Constructor
131   for(Int_t i = 0; i < 10;    i++)  fGeomMatrix[i] =  0;
132   for(Int_t j = 0; j < 24*48*11; j++)  {
133     fCellLabels[j]       = -1;
134     fCellSecondLabels[j] = -1;
135     fCellTime[j]         =  0.; 
136     fCellMatchdEta[j]    = -999;
137     fCellMatchdPhi[j]    = -999;
138   }
139   fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
140   fClusterArr      = new TObjArray(100);
141   fCaloClusterArr  = new TObjArray(100);
142   fRecParam        = new AliEMCALRecParam;
143   fBranchNames     = "ESD:AliESDHeader.,EMCALCells.";
144   fRecoUtils       = new AliEMCALRecoUtils();
145 }
146
147
148 //_______________________________________________________________
149 AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
150 {
151   //dtor 
152   
153   if (fDigitsArr){
154     fDigitsArr->Clear("C");
155     delete fDigitsArr; 
156   }
157   
158   if (fClusterArr){
159     fClusterArr->Delete();
160     delete fClusterArr;
161   }
162   
163   if (fCaloClusterArr){
164     fCaloClusterArr->Delete();
165     delete fCaloClusterArr; 
166   }
167   
168   if(fClusterizer) delete fClusterizer;
169   if(fUnfolder)    delete fUnfolder;   
170   if(fRecoUtils)   delete fRecoUtils;
171   
172 }
173
174 //_________________________________________________
175 Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
176 {
177   //Access to OCDB stuff
178   
179   fEvent = InputEvent();
180   if (!fEvent)
181   {
182     Warning("AccessOCDB","Event not available!!!");
183     return kFALSE;
184   }
185   
186   if (fEvent->GetRunNumber()==fRun)
187     return kTRUE;
188   fRun = fEvent->GetRunNumber();
189   
190   if(DebugLevel() > 1 )
191     printf("AliAnalysisTaksEMCALClusterize::AccessODCD() - Begin");
192   
193   //fGeom = AliEMCALGeometry::GetInstance(fGeomName);
194   
195   AliCDBManager *cdb = AliCDBManager::Instance();
196   
197   
198   if (fOCDBpath.Length()){
199     cdb->SetDefaultStorage(fOCDBpath.Data());
200     printf("AliAnalysisTaksEMCALClusterize::AccessOCDB() - Default storage %s",fOCDBpath.Data());
201   }
202   
203   cdb->SetRun(fEvent->GetRunNumber());
204   
205   //
206   // EMCAL from RAW OCDB
207   if (fOCDBpath.Contains("alien:"))
208   {
209     cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
210     cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
211   }
212   
213   TString path = cdb->GetDefaultStorage()->GetBaseFolder();
214   
215   // init parameters:
216   
217   //Get calibration parameters  
218   if(!fCalibData)
219   {
220     AliCDBEntry *entry = (AliCDBEntry*) 
221     AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
222     
223     if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
224   }
225   
226   if(!fCalibData)
227     AliFatal("Calibration parameters not found in CDB!");
228   
229   //Get calibration parameters  
230   if(!fPedestalData)
231   {
232     AliCDBEntry *entry = (AliCDBEntry*) 
233     AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
234     
235     if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
236   }
237   
238   if(!fPedestalData)
239     AliFatal("Dead map not found in CDB!");
240   
241   return kTRUE;
242 }
243
244 //_____________________________________________________
245 void AliAnalysisTaskEMCALClusterize::CheckAndGetEvent()
246 {
247   // Get the input event, it can depend in embedded events what you want to get
248   // Also check if the quality of the event is good if not reject it
249   
250   fEvent = 0x0;
251   
252   AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
253   Int_t eventN = Entry();
254   if(aodIH) eventN = aodIH->GetReadEntry(); 
255   
256   if (eventN > fMaxEvent) 
257     return ;
258   
259   //printf("Clusterizer --- Event %d-- \n",eventN);
260
261   //Check if input event are embedded events
262   //If so, take output event
263   if (aodIH && aodIH->GetMergeEvents()) {
264     fEvent  = AODEvent();
265     
266     if(!aodIH->GetMergeEMCALCells()) 
267       AliFatal("Events merged but not EMCAL cells, check analysis settings!");
268     
269     if(DebugLevel() > 1){
270      printf("AliAnalysisTaksEMCALClusterize::UserExec() - Use embedded events\n");
271         printf("\t InputEvent  N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
272                InputEvent()->GetEMCALCells()->GetNumberOfCells());
273         printf("\t MergedEvent  N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(),
274                aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells());
275         for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++) {
276             AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
277             if(sigCluster->IsEMCAL()) printf("\t \t Signal cluster: i %d, E  %f\n",icl,sigCluster->E());
278         }
279         printf("\t OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
280                AODEvent()->GetEMCALCells()->GetNumberOfCells());
281     }
282   }
283   else {
284     fEvent =  InputEvent();
285     if(fFillAODCaloCells) FillAODCaloCells();   
286     if(fFillAODHeader)    FillAODHeader();
287   }
288   
289   if (!fEvent) {
290     Error("UserExec","Event not available");
291     return ;
292   }
293   
294   //-------------------------------------------------------------------------------------
295   // Reject events if LED was firing, use only for LHC11a data 
296   // Reject event if triggered by exotic cell and remove exotic cells if not triggered
297   //-------------------------------------------------------------------------------------
298   
299   if(IsLEDEvent   ()) { fEvent = 0x0 ; return ; }
300   
301   if(IsExoticEvent()) { fEvent = 0x0 ; return ; }
302   
303   //Magic line to write events to AOD filem put after event rejection
304   AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
305   
306 }
307
308 //____________________________________________________
309 void AliAnalysisTaskEMCALClusterize::ClusterizeCells()
310
311   // Recluster calocells, transform them into digits, 
312   // feed the clusterizer with them and get new list of clusters
313   
314   //In case of MC, first loop on the clusters and fill MC label to array  
315   Int_t nClusters     = fEvent->GetNumberOfCaloClusters();
316   Int_t nClustersOrg  = 0;
317
318   AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
319   if(aodIH && aodIH->GetEventToMerge())  //Embedding
320     nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
321   
322   for (Int_t i = 0; i < nClusters; i++)
323   {
324     AliVCluster *clus = 0;
325     if(aodIH && aodIH->GetEventToMerge()) //Embedding
326       clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
327     else      
328       clus = fEvent->GetCaloCluster(i);
329     
330     if(!clus) return;
331     
332     if(clus->IsEMCAL()){   
333       Int_t label = clus->GetLabel();
334       Int_t label2 = -1 ;
335       //printf("Org cluster E %f, Time  %e, Id = ", clus->E(), clus->GetTOF() );
336       if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
337       UShort_t * index    = clus->GetCellsAbsId() ;
338       for(Int_t icell=0; icell < clus->GetNCells(); icell++ ){
339         fCellLabels[index[icell]]       = label;
340         fCellSecondLabels[index[icell]] = label2;
341         fCellTime[index[icell]]         = clus->GetTOF();    
342         fCellMatchdEta[index[icell]]    = clus->GetTrackDz();
343         fCellMatchdPhi[index[icell]]    = clus->GetTrackDx();
344         //printf(" %d,", index[icell] );
345       }
346       nClustersOrg++;
347     }
348      // printf("\n");
349   } 
350   
351   // Transform CaloCells into Digits
352
353   Int_t    idigit =  0;
354   Int_t    id     = -1;
355   Float_t  amp    = -1; 
356   Double_t time   = -1; 
357   
358   AliVCaloCells *cells   = fEvent->GetEMCALCells();
359     
360   TTree *digitsTree = new TTree("digitstree","digitstree");
361   digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
362   Int_t bc = InputEvent()->GetBunchCrossNumber();
363   for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
364   {
365     
366     // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
367     id = cells->GetCellNumber(icell);
368     Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,cells);
369     
370     // Do not include cells with too low energy, nor exotic cell
371     if(amp < fRecParam->GetMinECut() ) accept = kFALSE;
372     
373     // In case of AOD analysis cell time is 0, approximate replacing by time of the cluster the digit belongs.
374     if (time*1e9 < 1.) { 
375       time = fCellTime[id];
376       //printf("cell %d time cluster %e\n",id, time);
377       fRecoUtils->RecalibrateCellTime(id,bc,time);
378     }
379     
380     if(  accept && fRecoUtils->IsExoticCell(id,cells,bc)){
381       accept = kFALSE;
382     }
383
384     if( !accept ){
385       fCellLabels[id]      =-1; //reset the entry in the array for next event
386       fCellSecondLabels[id]=-1; //reset the entry in the array for next event
387       fCellTime[id]        = 0.; 
388       fCellMatchdEta[id]   =-999;
389       fCellMatchdPhi[id]   =-999;
390       if( DebugLevel() > 2 ) printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - Remove channel absId %d, index %d of %d, amp %f, time %f\n",
391                                     id,icell, cells->GetNumberOfCells(), amp, time*1.e9);
392       continue;
393     }
394             
395     //Create the digit, put a fake primary deposited energy to trick the clusterizer when checking the most likely primary
396     new((*fDigitsArr)[idigit]) AliEMCALDigit( fCellLabels[id], fCellLabels[id],id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, 1); 
397  
398     fCellLabels[id]      =-1; //reset the entry in the array for next event
399     
400     idigit++;
401   }
402   
403   //Fill the tree with digits
404   digitsTree->Fill();
405   
406   //-------------------------------------------------------------------------------------
407   //Do the clusterization
408   //-------------------------------------------------------------------------------------        
409   TTree *clustersTree = new TTree("clustertree","clustertree");
410   
411   fClusterizer->SetInput(digitsTree);
412   fClusterizer->SetOutput(clustersTree);
413   fClusterizer->Digits2Clusters("");
414   
415   //-------------------------------------------------------------------------------------
416   //Transform the recpoints into AliVClusters
417   //-------------------------------------------------------------------------------------
418   
419   clustersTree->SetBranchStatus("*",0); //disable all branches
420   clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
421   
422   TBranch *branch = clustersTree->GetBranch("EMCALECARP");
423   branch->SetAddress(&fClusterArr);
424   branch->GetEntry(0);
425   
426   RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
427   
428   if(!fCaloClusterArr){ 
429     printf("AliAnalysisTaksEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast());
430     return;    
431   }
432   
433   if( DebugLevel() > 0 ){
434     
435     printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - N clusters: before recluster %d, after recluster %d\n",nClustersOrg, fCaloClusterArr->GetEntriesFast());
436
437     if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast()){
438       printf("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
439              fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
440              fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
441     }
442   }
443   
444   //Reset the array with second labels for this event
445   memset(fCellSecondLabels, -1, sizeof(fCellSecondLabels));
446   
447   //---CLEAN UP-----
448   fClusterizer->Clear();
449   fDigitsArr  ->Clear("C");
450   fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
451   
452   clustersTree->Delete("all");
453   digitsTree  ->Delete("all");
454   
455 }
456
457 //_____________________________________________________
458 void AliAnalysisTaskEMCALClusterize::ClusterUnfolding()
459 {
460   // Take the event clusters and unfold them
461   
462   AliVCaloCells *cells   = fEvent->GetEMCALCells();
463   Double_t cellAmplitude = 0;
464   Double_t cellTime      = 0;
465   Short_t  cellNumber    = 0;
466   Int_t    nClustersOrg  = 0;
467   
468   // Fill the array with the EMCAL clusters, copy them
469   for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++)
470   {
471     AliVCluster *clus = fEvent->GetCaloCluster(i);
472     if(clus->IsEMCAL()){        
473       
474       //recalibrate/remove bad channels/etc if requested
475       if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells())){
476         continue;
477       } 
478       
479       if(fRecoUtils->IsRecalibrationOn()){
480         
481         //Calibrate cluster
482         fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
483         
484         //CalibrateCells
485         for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
486         {
487           if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
488             break;
489           
490           Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
491           fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta); 
492           fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);       
493           
494           //Do not include bad channels found in analysis?
495           if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() && 
496              fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
497             fCellLabels[cellNumber]      =-1; //reset the entry in the array for next event
498             fCellSecondLabels[cellNumber]=-1; //reset the entry in the array for next event
499             fCellTime[cellNumber]        = 0.;   
500             fCellMatchdEta[cellNumber]   =-999;
501             fCellMatchdPhi[cellNumber]   =-999;
502             continue;
503           }
504           
505           cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
506           
507         }// cells loop            
508       }// recalibrate
509       
510       //Cast to ESD or AOD, needed to create the cluster array
511       AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
512       AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
513       if     (esdCluster){
514         fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );   
515       }//ESD
516       else if(aodCluster){
517         fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );   
518       }//AOD
519       else 
520         Warning("UserExec()"," - Wrong CaloCluster type?");
521       nClustersOrg++;
522     }
523   }
524   
525   //Do the unfolding
526   fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
527   
528   //CLEAN-UP
529   fUnfolder->Clear();
530   
531 }
532
533 //_____________________________________________________
534 void AliAnalysisTaskEMCALClusterize::FillAODCaloCells() 
535 {
536   // Put calo cells in standard branch  
537   AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
538   Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
539   
540   AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
541   aodEMcells.CreateContainer(nEMcell);
542   aodEMcells.SetType(AliVCaloCells::kEMCALCell);
543   Double_t calibFactor = 1.;   
544   for (Int_t iCell = 0; iCell < nEMcell; iCell++) { 
545     Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
546     fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta); 
547     fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);     
548     
549     if(fRecoUtils->IsRecalibrationOn()){ 
550       calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
551     }
552     
553     if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ //Channel is not declared as bad
554       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor);
555     }
556     else {
557       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0);
558     }
559   }
560   aodEMcells.Sort();
561   
562 }
563
564 //__________________________________________________
565 void AliAnalysisTaskEMCALClusterize::FillAODHeader() 
566 {
567   //Put event header information in standard AOD branch
568   
569   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent);
570   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent);
571   
572   Double_t pos[3]   ;
573   Double_t covVtx[6];
574   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
575   
576   AliAODHeader* header = AODEvent()->GetHeader();
577   header->SetRunNumber(fEvent->GetRunNumber());
578   
579   if(esdevent){
580     TTree* tree = fInputHandler->GetTree();
581     if (tree) {
582       TFile* file = tree->GetCurrentFile();
583       if (file) header->SetESDFileName(file->GetName());
584     }
585   }
586   else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName());
587   
588   header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
589   header->SetOrbitNumber(fEvent->GetOrbitNumber());
590   header->SetPeriodNumber(fEvent->GetPeriodNumber());
591   header->SetEventType(fEvent->GetEventType());
592   
593   //Centrality
594   if(fEvent->GetCentrality()){
595     header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
596   }
597   else{
598     header->SetCentrality(0);
599   }
600   
601   //Trigger  
602   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
603   if      (esdevent) header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses());
604   else if (aodevent) header->SetFiredTriggerClasses(aodevent->GetFiredTriggerClasses());
605   header->SetTriggerMask(fEvent->GetTriggerMask()); 
606   header->SetTriggerCluster(fEvent->GetTriggerCluster());
607   if(esdevent){
608     header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());    
609     header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());    
610     header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());    
611   }
612   else if (aodevent){
613     header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());    
614     header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());    
615     header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());    
616   }
617   
618   header->SetMagneticField(fEvent->GetMagneticField());
619   //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.); 
620   
621   header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
622   header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
623   header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
624   header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
625   header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
626   
627   Float_t diamxy[2]={fEvent->GetDiamondX(),fEvent->GetDiamondY()};
628   Float_t diamcov[3];
629   fEvent->GetDiamondCovXY(diamcov);
630   header->SetDiamond(diamxy,diamcov);
631   if      (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
632   else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
633   
634   //
635   //
636   Int_t nVertices = 1 ;/* = prim. vtx*/;
637   Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
638   
639   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
640   
641   // Access to the AOD container of vertices
642   TClonesArray &vertices = *(AODEvent()->GetVertices());
643   Int_t jVertices=0;
644   
645   // Add primary vertex. The primary tracks will be defined
646   // after the loops on the composite objects (V0, cascades, kinks)
647   fEvent->GetPrimaryVertex()->GetXYZ(pos);
648   Float_t chi = 0;
649   if      (esdevent){
650     esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
651     chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
652   }
653   else if (aodevent){
654     aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
655     chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
656   }
657   
658   AliAODVertex * primary = new(vertices[jVertices++])
659   AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
660   primary->SetName(fEvent->GetPrimaryVertex()->GetName());
661   primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
662   
663 }
664
665 //_________________________________________
666 void AliAnalysisTaskEMCALClusterize::Init()
667 {
668   //Init analysis with configuration macro if available
669   
670   if(gROOT->LoadMacro(fConfigName) >=0){
671     printf("AliAnalysisTaksEMCALClusterize::Init() - Configure analysis with %s\n",fConfigName.Data());
672     AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
673     fGeomName         = clus->fGeomName; 
674     fLoadGeomMatrices = clus->fLoadGeomMatrices;
675     fOCDBpath         = clus->fOCDBpath;   
676     fAccessOCDB       = clus->fAccessOCDB;
677     fRecParam         = clus->fRecParam;
678     fJustUnfold       = clus->fJustUnfold;
679     fFillAODFile      = clus->fFillAODFile;
680     fRecoUtils        = clus->fRecoUtils; 
681     fConfigName       = clus->fConfigName;
682     fMaxEvent         = clus->fMaxEvent;
683     fDoTrackMatching  = clus->fDoTrackMatching;
684     fOutputAODBranchName = clus->fOutputAODBranchName;
685     for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
686     
687   }
688   
689   // Init geometry, I do not like much to do it like this ...
690   if(fImportGeometryFromFile && !gGeoManager) {
691     printf("AliAnalysisTaskEMCALClusterize::Init() - Import geometry.root file\n");
692     TGeoManager::Import(Form("%s/geometry.root", fImportGeometryFilePath.Data())) ; // default need file "geometry.root" in local dir!!!!
693   }
694   
695 }  
696
697 //_______________________________________________________
698 void AliAnalysisTaskEMCALClusterize::InitClusterization()
699 {
700   //Select clusterization/unfolding algorithm and set all the needed parameters
701   
702   if (fJustUnfold){
703     // init the unfolding afterburner 
704     delete fUnfolder;
705     fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
706     return;
707   }
708   
709   //First init the clusterizer
710   delete fClusterizer;
711   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
712     fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
713   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2) 
714     fClusterizer = new AliEMCALClusterizerv2(fGeom, fCalibData, fPedestalData);
715   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN){ 
716     fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
717     fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
718     fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
719   } else {
720     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
721   }
722   
723   //Now set the parameters
724   fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
725   fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
726   fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );    
727   fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
728   fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
729   fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
730   fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
731   fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
732   fClusterizer->SetInputCalibrated       ( kTRUE                               );
733   fClusterizer->SetJustClusters          ( kTRUE                               );  
734   //In case of unfolding after clusterization is requested, set the corresponding parameters
735   if(fRecParam->GetUnfold()){
736     Int_t i=0;
737     for (i = 0; i < 8; i++) {
738       fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
739     }//end of loop over parameters
740     for (i = 0; i < 3; i++) {
741       fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
742       fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
743     }//end of loop over parameters
744     
745     fClusterizer->InitClusterUnfolding();
746     
747   }// to unfold
748 }
749
750 //_________________________________________________
751 void AliAnalysisTaskEMCALClusterize::InitGeometry()
752 {
753   // Set the geometry matrix, for the first event, skip the rest
754   // Also set once the run dependent calibrations
755   
756   if(!fGeomMatrixSet){
757     if(fLoadGeomMatrices){
758       for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
759         if(fGeomMatrix[mod]){
760           if(DebugLevel() > 1) 
761             fGeomMatrix[mod]->Print();
762           fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;  
763         }
764         fGeomMatrixSet=kTRUE;
765       }//SM loop
766     }//Load matrices
767     else if(!gGeoManager){
768       printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - Get geo matrices from data");
769       //Still not implemented in AOD, just a workaround to be able to work at least with ESDs   
770       if(!strcmp(fEvent->GetName(),"AliAODEvent")) {
771         if(DebugLevel() > 1) 
772           Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
773       }//AOD
774       else {    
775         if(DebugLevel() > 1) 
776           printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
777         AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ;
778         if(!esd) {
779           Error("InitGeometry"," - This event does not contain ESDs?");
780           AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
781           return;
782         }
783         for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
784           if(DebugLevel() > 1) 
785             esd->GetEMCALMatrix(mod)->Print();
786           if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
787         } 
788         fGeomMatrixSet=kTRUE;
789       }//ESD
790     }//Load matrices from Data 
791     
792     //Recover time dependent corrections, put then in recalibration histograms. Do it once
793     fRecoUtils->SetRunDependentCorrections(InputEvent()->GetRunNumber());
794     
795   }//first event  
796   
797 }
798
799 //____________________________________________________
800 Bool_t AliAnalysisTaskEMCALClusterize::IsExoticEvent()
801 {
802   
803   // Check if event is exotic, get an exotic cell and compare with triggered patch
804   // If there is a match remove event ... to be completed, filled with something provisional
805   
806   if(!fRemoveExoticEvents) return kFALSE;
807   
808   // Loop on cells
809   AliVCaloCells * cells = fEvent->GetEMCALCells();
810   Float_t totCellE = 0;
811   Int_t bc = InputEvent()->GetBunchCrossNumber();
812   for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++){
813     
814     Float_t  ecell = 0 ;
815     Double_t tcell = 0 ;
816     
817     Int_t absID   = cells->GetCellNumber(icell);
818     Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,cells);
819     if(accept && !fRecoUtils->IsExoticCell(absID,cells,bc)) totCellE += ecell;
820   }
821   
822   //  TString triggerclasses = "";
823   //  if(esdevent) triggerclasses = esdevent             ->GetFiredTriggerClasses();
824   //  else         triggerclasses = ((AliAODEvent*)event)->GetFiredTriggerClasses();
825   //    //  
826   //    printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster  - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
827   //    AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
828   //    return;
829   //  
830   
831   //printf("TotE cell %f\n",totCellE);
832   if(totCellE < 1) return kTRUE;
833   
834   return kFALSE;
835   
836
837
838 //_________________________________________________
839 Bool_t AliAnalysisTaskEMCALClusterize::IsLEDEvent()
840 {
841   //Check if event is LED
842   
843   if(!fRemoveLEDEvents) return kFALSE;
844     
845   // Count number of cells with energy larger than 0.1 in SM3, cut on this number
846   Int_t ncellsSM3 = 0;
847   Int_t ncellsSM4 = 0;
848   AliVCaloCells * cells = fEvent->GetEMCALCells();
849   for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++){
850     if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
851     if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==4) ncellsSM4++;      
852   }
853   
854   TString triggerclasses = "";
855   
856   AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(fEvent);
857   if(esdevent) triggerclasses = esdevent              ->GetFiredTriggerClasses();
858   else         triggerclasses = ((AliAODEvent*)fEvent)->GetFiredTriggerClasses();
859   
860   Int_t ncellcut = 21;
861   if(triggerclasses.Contains("EMC")) ncellcut = 35;
862   
863   if( ncellsSM3 >= ncellcut || ncellsSM4 >= 100 ){
864     printf("AliAnalysisTaksEMCALClusterize::IsLEDEvent() - reject event %d with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
865     AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
866     return kTRUE;
867   }
868   
869   return kFALSE;
870   
871
872
873 //______________________________________________________________________________
874 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, 
875                                                         TObjArray *recPoints, 
876                                                         TObjArray *clusArray)
877 {
878   // Restore clusters from recPoints
879   // Cluster energy, global position, cells and their amplitude fractions are restored
880   Int_t j = 0;
881   for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++)
882   {
883     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i);
884     
885     const Int_t ncells = recPoint->GetMultiplicity();
886     Int_t ncellsTrue = 0;
887     
888     if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
889     
890     // cells and their amplitude fractions
891     UShort_t   absIds[ncells];  
892     Double32_t ratios[ncells];
893     
894     //For later check embedding
895     AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
896     
897     Float_t clusterE = 0; 
898     for (Int_t c = 0; c < ncells; c++) {
899       AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]);
900       
901       absIds[ncellsTrue] = digit->GetId();
902       ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
903       
904       // In case of unfolding, remove digits with unfolded energy too low      
905       if(fSelectCell){
906         if     (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)  {
907           
908           if(DebugLevel() > 1)  {
909             printf("AliAnalysisTaksEMCALClusterize::RecPoints2Clusters() - Too small energy in cell of cluster: cluster cell %f, digit %f\n",
910                    recPoint->GetEnergiesList()[c],digit->GetAmplitude());
911           }
912           
913           continue;
914           
915         } // if cuts
916       }// Select cells
917       
918       //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
919       ncellsTrue++;
920       clusterE  +=recPoint->GetEnergiesList()[c];
921       
922       // In case of embedding, fill ratio with amount of signal, 
923       if (aodIH && aodIH->GetMergeEvents()) {
924         
925         //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
926         AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
927         AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
928         
929         Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
930         //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
931         Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
932         //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
933         
934         if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
935         //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
936         
937       }//Embedding
938       
939     }// cluster cell loop
940     
941     if (ncellsTrue < 1) {
942       if (DebugLevel() > 1) 
943         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Skipping cluster with no cells avobe threshold E = %f, ncells %d\n",
944                recPoint->GetEnergy(), ncells);
945       continue;
946     }
947     
948     //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
949     
950     if(clusterE <  fRecParam->GetClusteringThreshold()) {
951       if (DebugLevel()>1)
952         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Remove cluster with energy below seed threshold %f\n",clusterE);
953       continue;
954     }
955     
956     TVector3 gpos;
957     Float_t g[3];
958     
959     // calculate new cluster position
960     recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr);
961     recPoint->GetGlobalPosition(gpos);
962     gpos.GetXYZ(g);
963     
964     // create a new cluster
965     (*clusArray)[j] = new AliAODCaloCluster() ;
966     AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( clusArray->At(j) ) ;
967     j++;
968     clus->SetType(AliVCluster::kEMCALClusterv1);
969     clus->SetE(clusterE);
970     clus->SetPosition(g);
971     clus->SetNCells(ncellsTrue);
972     clus->SetCellsAbsId(absIds);
973     clus->SetCellsAmplitudeFraction(ratios);
974     clus->SetChi2(-1); //not yet implemented
975     clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
976     clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
977     clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower()); 
978
979     if(ncells == ncellsTrue){
980       Float_t elipAxis[2];
981       recPoint->GetElipsAxis(elipAxis);
982       clus->SetM02(elipAxis[0]*elipAxis[0]) ;
983       clus->SetM20(elipAxis[1]*elipAxis[1]) ;
984       clus->SetDispersion(recPoint->GetDispersion());
985     }
986     else if(fSelectCell){
987       // In case some cells rejected, in unfolding case, recalculate
988       // shower shape parameters and position
989       AliVCaloCells* cells = 0x0; 
990       if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent()  ->GetEMCALCells();
991       else                                  cells = InputEvent()->GetEMCALCells();
992       fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus);
993       fRecoUtils->RecalculateClusterPID(clus);
994       fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus); 
995       
996     }
997     
998     //MC
999     Int_t  parentMult = 0;
1000     Int_t *parentList = recPoint->GetParents(parentMult);
1001     clus->SetLabel(parentList, parentMult); 
1002     
1003     //Write the second major contributor to each MC cluster.
1004     Int_t iNewLabel ;
1005     for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ ){
1006       
1007       Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1008       if(idCell>=0){
1009         iNewLabel = 1 ; //iNewLabel makes sure we  don't write twice the same label.
1010         for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
1011         {
1012           if ( fCellSecondLabels[idCell] == -1 )  iNewLabel = 0;  // -1 is never a good second label.
1013           if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) )  iNewLabel = 0;
1014         }
1015         if (iNewLabel == 1) 
1016         {
1017           Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
1018           for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ ){
1019             newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
1020           }
1021           
1022           newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
1023           clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
1024           delete [] newLabelArray;
1025         }
1026       }//positive cell number
1027     }
1028     
1029   } // recPoints loop
1030   
1031 }
1032
1033 //____________________________________________________________
1034 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
1035 {
1036   // Init geometry, create list of output clusters
1037   
1038   fGeom =  AliEMCALGeometry::GetInstance(fGeomName) ;   
1039   if(fOutputAODBranchName.Length()!=0){
1040     fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1041     fOutputAODBranch->SetName(fOutputAODBranchName);
1042     //fOutputAODBranch->SetOwner(kFALSE);
1043     AddAODBranch("TClonesArray", &fOutputAODBranch);
1044   }
1045   else {
1046     AliFatal("fOutputAODBranchName not set\n");
1047   }
1048   
1049   //PostData(0,fOutputAODBranch);
1050   
1051 }
1052
1053 //_______________________________________________________
1054 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *) 
1055 {
1056   // Main loop
1057   // Called for each event
1058     
1059   //Remove the contents of output list set in the previous event 
1060   fOutputAODBranch->Clear("C");
1061   
1062   LoadBranches();
1063   
1064   //Get the event, do some checks and settings
1065   CheckAndGetEvent() ;
1066   
1067   if (!fEvent) {
1068     if(DebugLevel() > 0 ) printf("AliAnalysisTaksEMCALClusterize::UserExec() - Skip Event %d", (Int_t) Entry());
1069     return ;
1070   }
1071   
1072   //Init pointers, clusterizer, ocdb
1073   if(fAccessOCDB) AccessOCDB();
1074   
1075   InitClusterization();
1076   
1077   InitGeometry(); // only once
1078   
1079   // Make clusters
1080   if (fJustUnfold) ClusterUnfolding();
1081   else             ClusterizeCells() ;
1082   
1083   //Recalculate track-matching for the new clusters
1084   if(fDoTrackMatching) fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom);
1085   
1086   //Put the new clusters in the AOD list
1087   
1088   Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntriesFast();
1089   for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
1090     AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
1091     
1092     //Add matched track
1093     if(fDoTrackMatching){
1094       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
1095       if(trackIndex >= 0){
1096         newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex));
1097         if(DebugLevel() > 1) 
1098           printf("AliAnalysisTaksEMCALClusterize::UserExec() - Matched Track index %d to new cluster %d \n",trackIndex,i);
1099       }
1100       
1101       Float_t dR = 999., dZ = 999.;
1102       fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dR,dZ);
1103       newCluster->SetTrackDistance(dR,dZ);
1104       
1105     }
1106     else {// Assign previously assigned matched track in reco, very very rough
1107       Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster
1108       newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]);
1109     }
1110
1111     //printf("New cluster E %f, Time  %e, Id = ", newCluster->E(), newCluster->GetTOF() );
1112     //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] );
1113     //printf("\n");
1114
1115     //In case of new bad channels, recalculate distance to bad channels
1116     if(fRecoUtils->IsBadChannelsRemovalSwitchedOn())
1117       fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fEvent->GetEMCALCells(), newCluster);
1118
1119     newCluster->SetID(i);
1120     new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
1121         
1122     if(DebugLevel() > 1 )    
1123       printf("AliAnalysisTaksEMCALClusterize::UserExec() - New cluster %d of %d, energy %f\n",newCluster->GetID(), kNumberOfCaloClusters, newCluster->E());
1124     
1125   } // cluster loop
1126   
1127   fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
1128   
1129   // Clean up
1130   fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
1131   
1132   //PostData(0,fOutputAODBranch);
1133   
1134 }      
1135