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