remove print
[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 //
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
50 // --- EMCAL
51 #include "AliEMCALAfterBurnerUF.h"
52 #include "AliEMCALGeometry.h"
53 #include "AliEMCALClusterizerNxN.h"
54 #include "AliEMCALClusterizerv1.h"
55 #include "AliEMCALClusterizerv2.h"
56 #include "AliEMCALRecPoint.h"
57 #include "AliEMCALDigit.h"
58 #include "AliCaloCalibPedestal.h"
59 #include "AliEMCALCalibData.h"
60
61 #include "AliAnalysisTaskEMCALClusterize.h"
62
63 ClassImp(AliAnalysisTaskEMCALClusterize)
64
65 //______________________________________________________________________________
66 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name) 
67 : AliAnalysisTaskSE(name)
68 , fEvent(0)
69 , fGeom(0),               fGeomName("") 
70 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
71 , fCalibData(0),          fPedestalData(0)
72 , fOCDBpath(""),          fAccessOCDB(kFALSE)
73 , fDigitsArr(0),          fClusterArr(0),             fCaloClusterArr(0)
74 , fRecParam(0),           fClusterizer(0)
75 , fUnfolder(0),           fJustUnfold(kFALSE) 
76 , fOutputAODBranch(0),    fOutputAODBranchName("")
77 , fFillAODFile(kTRUE),    fFillAODHeader(0)
78 , fFillAODCaloCells(0),   fRun(-1)
79 , fRecoUtils(0),          fConfigName("")
80 , fCellLabels(),          fCellSecondLabels(),        fCellTime()
81 , fCellMatchdEta(),       fCellMatchdPhi()
82 , fMaxEvent(0),           fDoTrackMatching(kFALSE)
83 , fSelectCell(kFALSE),    fSelectCellMinE(0),         fSelectCellMinFrac(0)
84 , fRemoveLEDEvents(kTRUE), fRemoveExoticEvents(kFALSE)
85 , fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("") 
86 , fOADBSet(kFALSE),       fAccessOADB(kTRUE),         fOADBFilePath("")           
87 {
88   //ctor
89   for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
90   for(Int_t j = 0; j < 24*48*11; j++)  
91   {
92     fCellLabels[j]       = -1;
93     fCellSecondLabels[j] = -1;
94     fCellTime[j]         =  0.;    
95     fCellMatchdEta[j]    = -999;
96     fCellMatchdPhi[j]    = -999;
97   }  
98   
99   
100 }
101
102 //______________________________________________________________
103 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize() 
104 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
105 , fEvent(0)
106 , fGeom(0),                 fGeomName("") 
107 , fGeomMatrixSet(kFALSE),   fLoadGeomMatrices(kFALSE)
108 , fCalibData(0),            fPedestalData(0)
109 , fOCDBpath(""),            fAccessOCDB(kFALSE)
110 , fDigitsArr(0),            fClusterArr(0),             fCaloClusterArr(0)
111 , fRecParam(0),             fClusterizer(0)
112 , fUnfolder(0),             fJustUnfold(kFALSE) 
113 , fOutputAODBranch(0),      fOutputAODBranchName("")
114 , fFillAODFile(kTRUE),      fFillAODHeader(0)
115 , fFillAODCaloCells(0),     fRun(-1)
116 , fRecoUtils(0),            fConfigName("")
117 , fCellLabels(),            fCellSecondLabels(),        fCellTime()
118 , fCellMatchdEta(),         fCellMatchdPhi()
119 , fMaxEvent(0),             fDoTrackMatching(kFALSE)
120 , fSelectCell(kFALSE),      fSelectCellMinE(0),         fSelectCellMinFrac(0)
121 , fRemoveLEDEvents(kTRUE), fRemoveExoticEvents(kFALSE)
122 , fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("")
123 , fOADBSet(kFALSE),         fAccessOADB(kTRUE),        fOADBFilePath("")           
124
125 {
126   // Constructor
127   for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
128   for(Int_t j = 0; j < 24*48*11; j++)  
129   {
130     fCellLabels[j]       = -1;
131     fCellSecondLabels[j] = -1;
132     fCellTime[j]         =  0.; 
133     fCellMatchdEta[j]    = -999;
134     fCellMatchdPhi[j]    = -999;
135   }
136
137 }
138
139
140 //_______________________________________________________________
141 AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
142 {
143   //dtor 
144   
145   if (fDigitsArr){
146     fDigitsArr->Clear("C");
147     delete fDigitsArr; 
148   }
149   
150   if (fClusterArr){
151     fClusterArr->Delete();
152     delete fClusterArr;
153   }
154   
155   if (fCaloClusterArr){
156     fCaloClusterArr->Delete();
157     delete fCaloClusterArr; 
158   }
159   
160   if(fClusterizer) delete fClusterizer;
161   if(fUnfolder)    delete fUnfolder;   
162   if(fRecoUtils)   delete fRecoUtils;
163   
164 }
165
166 //_______________________________________________
167 void AliAnalysisTaskEMCALClusterize::AccessOADB()
168 {
169   // Set the AODB calibration, bad channels etc. parameters at least once
170   // alignment matrices from OADB done in SetGeometryMatrices
171   
172   //Set it only once
173   if(fOADBSet) return ; 
174   
175   Int_t   runnumber = InputEvent()->GetRunNumber() ;
176   TString pass      = GetPass();
177   
178   printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Get AODB parameters from EMCAL in %s for run %d, and <%s> \n",fOADBFilePath.Data(),runnumber,pass.Data());
179   
180   Int_t nSM = fGeom->GetNumberOfSuperModules();
181   
182   // Bad map
183   if(fRecoUtils->IsBadChannelsRemovalSwitchedOn())
184   {
185     AliOADBContainer *contBC=new AliOADBContainer("");
186     contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePath.Data()),"AliEMCALBadChannels"); 
187     
188     TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber);
189     
190     if(arrayBC)
191     {
192       TObjArray * arrayBCpass = 0x0; 
193       if(pass!="")arrayBCpass = (TObjArray*)arrayBC->FindObject(pass);
194       // There are no passes for simulation, in order to get the bad map put pass1
195       else        arrayBCpass = (TObjArray*)arrayBC->FindObject("pass1");
196       
197       if(arrayBCpass)
198       {
199         fRecoUtils->SwitchOnDistToBadChannelRecalculation();
200         printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Remove EMCAL bad cells \n");
201         
202         for (Int_t i=0; i<nSM; ++i) 
203         {
204           TH2I *hbm = fRecoUtils->GetEMCALChannelStatusMap(i);
205           
206           if (hbm)
207             delete hbm;
208           
209           hbm=(TH2I*)arrayBCpass->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
210           
211           if (!hbm) 
212           {
213             AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
214             continue;
215           }
216           
217           hbm->SetDirectory(0);
218           fRecoUtils->SetEMCALChannelStatusMap(i,hbm);
219           
220         } // loop
221       } else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT remove EMCAL bad channels 1\n"); // pass array
222     } else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT remove EMCAL bad channels 2\n"); // run array
223   }  // Remove bad
224   
225   // Energy Recalibration
226   if(fRecoUtils->IsRecalibrationOn())
227   {
228     AliOADBContainer *contRF=new AliOADBContainer("");
229     
230     contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePath.Data()),"AliEMCALRecalib");
231     
232     TObjArray *recal=(TObjArray*)contRF->GetObject(runnumber); 
233     
234     if(recal)
235     {
236       TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
237       
238       if(recalpass)
239       {
240         TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
241         
242         if(recalib)
243         {
244           printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Recalibrate EMCAL \n");
245           for (Int_t i=0; i<nSM; ++i) 
246           {
247             TH2F *h = fRecoUtils->GetEMCALChannelRecalibrationFactors(i);
248             
249             if (h)
250               delete h;
251             
252             h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
253             
254             if (!h) 
255             {
256               AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
257               continue;
258             }
259             
260             h->SetDirectory(0);
261             
262             fRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
263           } // SM loop
264         }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL 1\n"); // array ok
265       }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL 2\n"); // array pass ok
266     }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL 3\n");  // run number array ok
267     
268     // once set, apply run dependent corrections if requested
269     fRecoUtils->SetRunDependentCorrections(runnumber);
270     
271   } // Recalibration on
272   
273   
274   // Time Recalibration
275   if(fRecoUtils->IsTimeRecalibrationOn())
276   {
277     AliOADBContainer *contTRF=new AliOADBContainer("");
278     
279     contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePath.Data()),"AliEMCALTimeCalib");
280     
281     TObjArray *trecal=(TObjArray*)contTRF->GetObject(runnumber); 
282     
283     if(trecal)
284     {
285       TObjArray *trecalpass=(TObjArray*)trecal->FindObject(pass);
286       
287       if(trecalpass)
288       {
289         printf("AliCalorimeterUtils::SetOADBParameters() - Time Recalibrate EMCAL \n");
290         for (Int_t ibc = 0; ibc < 4; ++ibc) 
291         {
292           TH1F *h = fRecoUtils->GetEMCALChannelTimeRecalibrationFactors(ibc);
293           
294           if (h)
295             delete h;
296           
297           h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
298           
299           if (!h) 
300           {
301             AliError(Form("Could not load hAllTimeAvBC%d",ibc));
302             continue;
303           }
304           
305           h->SetDirectory(0);
306           
307           fRecoUtils->SetEMCALChannelTimeRecalibrationFactors(ibc,h);
308         } // bunch crossing loop
309       }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate time EMCAL 1\n"); // array pass ok
310     }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate time EMCAL 2\n");  // run number array ok
311     
312   } // Time recalibration on    
313   
314   // Parameters already set once, so do not it again
315   fOADBSet = kTRUE;
316   
317 }  
318
319 //_________________________________________________
320 Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
321 {
322   //Access to OCDB stuff
323   
324   fEvent = InputEvent();
325   if (!fEvent)
326   {
327     Warning("AccessOCDB","Event not available!!!");
328     return kFALSE;
329   }
330   
331   if (fEvent->GetRunNumber()==fRun)
332     return kTRUE;
333   fRun = fEvent->GetRunNumber();
334   
335   if(DebugLevel() > 1 )
336     printf("AliAnalysisTaksEMCALClusterize::AccessODCD() - Begin");
337     
338   AliCDBManager *cdb = AliCDBManager::Instance();
339   
340   
341   if (fOCDBpath.Length()){
342     cdb->SetDefaultStorage(fOCDBpath.Data());
343     printf("AliAnalysisTaksEMCALClusterize::AccessOCDB() - Default storage %s",fOCDBpath.Data());
344   }
345   
346   cdb->SetRun(fEvent->GetRunNumber());
347   
348   //
349   // EMCAL from RAW OCDB
350   if (fOCDBpath.Contains("alien:"))
351   {
352     cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
353     cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
354   }
355   
356   TString path = cdb->GetDefaultStorage()->GetBaseFolder();
357   
358   // init parameters:
359   
360   //Get calibration parameters  
361   if(!fCalibData)
362   {
363     AliCDBEntry *entry = (AliCDBEntry*) 
364     AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
365     
366     if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
367   }
368   
369   if(!fCalibData)
370     AliFatal("Calibration parameters not found in CDB!");
371   
372   //Get calibration parameters  
373   if(!fPedestalData)
374   {
375     AliCDBEntry *entry = (AliCDBEntry*) 
376     AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
377     
378     if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
379   }
380   
381   if(!fPedestalData)
382     AliFatal("Dead map not found in CDB!");
383   
384   return kTRUE;
385 }
386
387 //_____________________________________________________
388 void AliAnalysisTaskEMCALClusterize::CheckAndGetEvent()
389 {
390   // Get the input event, it can depend in embedded events what you want to get
391   // Also check if the quality of the event is good if not reject it
392   
393   fEvent = 0x0;
394   
395   AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
396   Int_t eventN = Entry();
397   if(aodIH) eventN = aodIH->GetReadEntry(); 
398   
399   if (eventN > fMaxEvent) 
400     return ;
401   
402   //printf("Clusterizer --- Event %d-- \n",eventN);
403
404   //Check if input event are embedded events
405   //If so, take output event
406   if (aodIH && aodIH->GetMergeEvents()) 
407   {
408     fEvent  = AODEvent();
409     
410     if(!aodIH->GetMergeEMCALCells()) 
411       AliFatal("Events merged but not EMCAL cells, check analysis settings!");
412     
413     if(DebugLevel() > 1)
414     {
415       printf("AliAnalysisTaksEMCALClusterize::UserExec() - Use embedded events\n");
416       
417       printf("\t InputEvent  N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
418              InputEvent()->GetEMCALCells()->GetNumberOfCells());
419       
420       printf("\t MergedEvent  N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(),
421              aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells());
422       
423       for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++) 
424       {
425         AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
426         if(sigCluster->IsEMCAL()) printf("\t \t Signal cluster: i %d, E  %f\n",icl,sigCluster->E());
427       }
428       
429       printf("\t OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
430              AODEvent()->GetEMCALCells()->GetNumberOfCells());
431     }
432   }
433   else 
434   {
435     fEvent =  InputEvent();
436     if(fFillAODCaloCells) FillAODCaloCells();   
437     if(fFillAODHeader)    FillAODHeader();
438   }
439   
440   if (!fEvent) 
441   {
442     Error("UserExec","Event not available");
443     return ;
444   }
445   
446   //-------------------------------------------------------------------------------------
447   // Reject events if LED was firing, use only for LHC11a data 
448   // Reject event if triggered by exotic cell and remove exotic cells if not triggered
449   //-------------------------------------------------------------------------------------
450   
451   if( IsLEDEvent( InputEvent()->GetRunNumber() ) ) { fEvent = 0x0 ; return ; }
452   
453   if( IsExoticEvent() )                            { fEvent = 0x0 ; return ; }
454   
455   //Magic line to write events to AOD filem put after event rejection
456   AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
457   
458 }
459
460 //____________________________________________________
461 void AliAnalysisTaskEMCALClusterize::ClusterizeCells()
462
463   // Recluster calocells, transform them into digits, 
464   // feed the clusterizer with them and get new list of clusters
465   
466   //In case of MC, first loop on the clusters and fill MC label to array  
467   Int_t nClusters     = fEvent->GetNumberOfCaloClusters();
468   Int_t nClustersOrg  = 0;
469
470   AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
471   if(aodIH && aodIH->GetEventToMerge())  //Embedding
472     nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
473   
474   for (Int_t i = 0; i < nClusters; i++)
475   {
476     AliVCluster *clus = 0;
477     if(aodIH && aodIH->GetEventToMerge()) //Embedding
478       clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
479     else      
480       clus = fEvent->GetCaloCluster(i);
481     
482     if(!clus) return;
483     
484     if(clus->IsEMCAL())
485     {   
486       Int_t label = clus->GetLabel();
487       Int_t label2 = -1 ;
488       //printf("Org cluster E %f, Time  %e, Id = ", clus->E(), clus->GetTOF() );
489       if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
490       UShort_t * index    = clus->GetCellsAbsId() ;
491       for(Int_t icell=0; icell < clus->GetNCells(); icell++ )
492       {
493         fCellLabels[index[icell]]       = label;
494         fCellSecondLabels[index[icell]] = label2;
495         fCellTime[index[icell]]         = clus->GetTOF();    
496         fCellMatchdEta[index[icell]]    = clus->GetTrackDz();
497         fCellMatchdPhi[index[icell]]    = clus->GetTrackDx();
498         //printf(" %d,", index[icell] );
499       }
500       nClustersOrg++;
501     }
502      // printf("\n");
503   } 
504   
505   // Transform CaloCells into Digits
506
507   Int_t    idigit =  0;
508   Int_t    id     = -1;
509   Float_t  amp    = -1; 
510   Double_t time   = -1; 
511   
512   AliVCaloCells *cells   = fEvent->GetEMCALCells();
513   
514   TFile* file = new TFile("tmpTree.root","recreate"); // Trick to avoid annoying messages FIXME
515   TTree *digitsTree = new TTree("digitstree","digitstree");
516   digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
517
518   Int_t bc = InputEvent()->GetBunchCrossNumber();
519   for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
520   {
521     
522     // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
523     id = cells->GetCellNumber(icell);
524     Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,cells);
525     
526     // Do not include cells with too low energy, nor exotic cell
527     if(amp < fRecParam->GetMinECut() ) accept = kFALSE;
528     
529     // In case of AOD analysis cell time is 0, approximate replacing by time of the cluster the digit belongs.
530     if (time*1e9 < 1.) 
531     { 
532       time = fCellTime[id];
533       //printf("cell %d time org %f - ",id, time*1.e9);
534       fRecoUtils->RecalibrateCellTime(id,bc,time);
535       //printf("recal %f\n",time*1.e9);
536     }
537     
538     if(  accept && fRecoUtils->IsExoticCell(id,cells,bc))
539     {
540       accept = kFALSE;
541     }
542
543     if( !accept )
544     {
545       fCellLabels[id]      =-1; //reset the entry in the array for next event
546       fCellSecondLabels[id]=-1; //reset the entry in the array for next event
547       fCellTime[id]        = 0.; 
548       fCellMatchdEta[id]   =-999;
549       fCellMatchdPhi[id]   =-999;
550       if( DebugLevel() > 2 )
551         printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - Remove channel absId %d, index %d of %d, amp %f, time %f\n",
552                                     id,icell, cells->GetNumberOfCells(), amp, time*1.e9);
553       continue;
554     }
555             
556     //Create the digit, put a fake primary deposited energy to trick the clusterizer when checking the most likely primary
557     new((*fDigitsArr)[idigit]) AliEMCALDigit( fCellLabels[id], fCellLabels[id],id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, 1); 
558  
559     fCellLabels[id]      =-1; //reset the entry in the array for next event
560     
561     idigit++;
562   }
563   
564   //Fill the tree with digits
565   digitsTree->Fill();
566   
567   //-------------------------------------------------------------------------------------
568   //Do the clusterization
569   //-------------------------------------------------------------------------------------        
570   TTree *clustersTree = new TTree("clustertree","clustertree");
571   
572   fClusterizer->SetInput(digitsTree);
573   fClusterizer->SetOutput(clustersTree);
574   fClusterizer->Digits2Clusters("");
575   
576   //-------------------------------------------------------------------------------------
577   //Transform the recpoints into AliVClusters
578   //-------------------------------------------------------------------------------------
579   
580   clustersTree->SetBranchStatus("*",0); //disable all branches
581   clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
582   
583   TBranch *branch = clustersTree->GetBranch("EMCALECARP");
584   branch->SetAddress(&fClusterArr);
585
586   branch->GetEntry(0);
587
588   RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
589   
590   if(!fCaloClusterArr)
591   { 
592     printf("AliAnalysisTaksEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast());
593     return;    
594   }
595   
596   if( DebugLevel() > 0 )
597   {
598     printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - N clusters: before recluster %d, after recluster %d\n",nClustersOrg, fCaloClusterArr->GetEntriesFast());
599
600     if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast()){
601       printf("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
602              fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
603              fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
604     }
605   }
606   
607   //Reset the array with second labels for this event
608   memset(fCellSecondLabels, -1, sizeof(fCellSecondLabels));
609   
610   //---CLEAN UP-----
611   fClusterizer->Clear();
612   fDigitsArr  ->Clear("C");
613   fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
614   
615   clustersTree->Delete("all");  
616   digitsTree  ->Delete("all");  
617   
618   // Trick to avoid annoying messages FIXME
619   file->Close();
620   delete file;
621 }
622
623 //_____________________________________________________
624 void AliAnalysisTaskEMCALClusterize::ClusterUnfolding()
625 {
626   // Take the event clusters and unfold them
627   
628   AliVCaloCells *cells   = fEvent->GetEMCALCells();
629   Double_t cellAmplitude = 0;
630   Double_t cellTime      = 0;
631   Short_t  cellNumber    = 0;
632   Int_t    nClustersOrg  = 0;
633   
634   // Fill the array with the EMCAL clusters, copy them
635   for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++)
636   {
637     AliVCluster *clus = fEvent->GetCaloCluster(i);
638     if(clus->IsEMCAL()){        
639       
640       //recalibrate/remove bad channels/etc if requested
641       if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells())){
642         continue;
643       } 
644       
645       if(fRecoUtils->IsRecalibrationOn()){
646         
647         //Calibrate cluster
648         fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
649         
650         //CalibrateCells
651         for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
652         {
653           if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
654             break;
655           
656           Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
657           fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta); 
658           fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);       
659           
660           //Do not include bad channels found in analysis?
661           if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() && 
662              fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
663             fCellLabels[cellNumber]      =-1; //reset the entry in the array for next event
664             fCellSecondLabels[cellNumber]=-1; //reset the entry in the array for next event
665             fCellTime[cellNumber]        = 0.;   
666             fCellMatchdEta[cellNumber]   =-999;
667             fCellMatchdPhi[cellNumber]   =-999;
668             continue;
669           }
670           
671           cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
672           
673         }// cells loop            
674       }// recalibrate
675       
676       //Cast to ESD or AOD, needed to create the cluster array
677       AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
678       AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
679       if     (esdCluster){
680         fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );   
681       }//ESD
682       else if(aodCluster){
683         fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );   
684       }//AOD
685       else 
686         Warning("UserExec()"," - Wrong CaloCluster type?");
687       nClustersOrg++;
688     }
689   }
690   
691   //Do the unfolding
692   fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
693   
694   //CLEAN-UP
695   fUnfolder->Clear();
696   
697 }
698
699 //_____________________________________________________
700 void AliAnalysisTaskEMCALClusterize::FillAODCaloCells() 
701 {
702   // Put calo cells in standard branch  
703   AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
704   Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
705   
706   AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
707   aodEMcells.CreateContainer(nEMcell);
708   aodEMcells.SetType(AliVCaloCells::kEMCALCell);
709   Double_t calibFactor = 1.;   
710   for (Int_t iCell = 0; iCell < nEMcell; iCell++) { 
711     Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
712     fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta); 
713     fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);     
714     
715     if(fRecoUtils->IsRecalibrationOn()){ 
716       calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
717     }
718     
719     if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ //Channel is not declared as bad
720       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor);
721     }
722     else {
723       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0);
724     }
725   }
726   aodEMcells.Sort();
727   
728 }
729
730 //__________________________________________________
731 void AliAnalysisTaskEMCALClusterize::FillAODHeader() 
732 {
733   //Put event header information in standard AOD branch
734   
735   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent);
736   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent);
737   
738   Double_t pos[3]   ;
739   Double_t covVtx[6];
740   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
741   
742   AliAODHeader* header = AODEvent()->GetHeader();
743   header->SetRunNumber(fEvent->GetRunNumber());
744   
745   if(esdevent){
746     TTree* tree = fInputHandler->GetTree();
747     if (tree) {
748       TFile* file = tree->GetCurrentFile();
749       if (file) header->SetESDFileName(file->GetName());
750     }
751   }
752   else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName());
753   
754   header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
755   header->SetOrbitNumber(fEvent->GetOrbitNumber());
756   header->SetPeriodNumber(fEvent->GetPeriodNumber());
757   header->SetEventType(fEvent->GetEventType());
758   
759   //Centrality
760   if(fEvent->GetCentrality()){
761     header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
762   }
763   else{
764     header->SetCentrality(0);
765   }
766   
767   //Trigger  
768   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
769   if      (esdevent) header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses());
770   else if (aodevent) header->SetFiredTriggerClasses(aodevent->GetFiredTriggerClasses());
771   header->SetTriggerMask(fEvent->GetTriggerMask()); 
772   header->SetTriggerCluster(fEvent->GetTriggerCluster());
773   if(esdevent){
774     header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());    
775     header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());    
776     header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());    
777   }
778   else if (aodevent){
779     header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());    
780     header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());    
781     header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());    
782   }
783   
784   header->SetMagneticField(fEvent->GetMagneticField());
785   //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.); 
786   
787   header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
788   header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
789   header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
790   header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
791   header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
792   
793   Float_t diamxy[2]={fEvent->GetDiamondX(),fEvent->GetDiamondY()};
794   Float_t diamcov[3];
795   fEvent->GetDiamondCovXY(diamcov);
796   header->SetDiamond(diamxy,diamcov);
797   if      (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
798   else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
799   
800   //
801   //
802   Int_t nVertices = 1 ;/* = prim. vtx*/;
803   Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
804   
805   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
806   
807   // Access to the AOD container of vertices
808   TClonesArray &vertices = *(AODEvent()->GetVertices());
809   Int_t jVertices=0;
810   
811   // Add primary vertex. The primary tracks will be defined
812   // after the loops on the composite objects (V0, cascades, kinks)
813   fEvent->GetPrimaryVertex()->GetXYZ(pos);
814   Float_t chi = 0;
815   if      (esdevent){
816     esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
817     chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
818   }
819   else if (aodevent){
820     aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
821     chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
822   }
823   
824   AliAODVertex * primary = new(vertices[jVertices++])
825   AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
826   primary->SetName(fEvent->GetPrimaryVertex()->GetName());
827   primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
828   
829 }
830
831 //_______________________________________________
832 TString AliAnalysisTaskEMCALClusterize::GetPass()
833 {
834   // Get passx from filename.
835   
836   if (!AliAnalysisManager::GetAnalysisManager()->GetTree()) 
837   {
838     AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Pointer to tree = 0, returning null\n");
839     return TString("");
840   }
841   
842   if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()) 
843   {
844     AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Null pointer input file, returning null\n");
845     return TString("");
846   }
847   
848   TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
849   if      (pass.Contains("ass1")) return TString("pass1");
850   else if (pass.Contains("ass2")) return TString("pass2");
851   else if (pass.Contains("ass3")) return TString("pass3");
852   
853   // No condition fullfilled, give a default value
854   printf("AliAnalysisTaskEMCALClusterize::GetPass() - Pass number string not found \n");
855   return TString("");            
856   
857 }
858
859 //_________________________________________
860 void AliAnalysisTaskEMCALClusterize::Init()
861 {
862   //Init analysis with configuration macro if available
863   
864   fOADBSet           = kFALSE;
865   if(fOADBFilePath == "") fOADBFilePath = "$ALICE_ROOT/OADB/EMCAL" ;          
866   
867   fDigitsArr         = new TClonesArray("AliEMCALDigit",12000);
868   fClusterArr        = new TObjArray(10000);
869   fCaloClusterArr    = new TObjArray(10000);
870   fBranchNames       = "ESD:AliESDHeader.,EMCALCells.";
871   
872   if(!fRecParam)     fRecParam  = new AliEMCALRecParam;
873   if(!fRecoUtils)    fRecoUtils = new AliEMCALRecoUtils();  
874   
875   if(fMaxEvent          <= 0) fMaxEvent          = 1000000000;
876   if(fSelectCellMinE    <= 0) fSelectCellMinE    = 0.005;     
877   if(fSelectCellMinFrac <= 0) fSelectCellMinFrac = 0.001;
878   
879   if (fOCDBpath            == "") fOCDBpath            = "raw://" ;
880   if (fOutputAODBranchName == "") fOutputAODBranchName = "newEMCALClusters" ;
881   
882   if(gROOT->LoadMacro(fConfigName) >=0)
883   {
884     printf("AliAnalysisTaksEMCALClusterize::Init() - Configure analysis with %s\n",fConfigName.Data());
885     AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
886     fGeomName         = clus->fGeomName; 
887     fLoadGeomMatrices = clus->fLoadGeomMatrices;
888     fOCDBpath         = clus->fOCDBpath;   
889     fAccessOCDB       = clus->fAccessOCDB;
890     fRecParam         = clus->fRecParam;
891     fJustUnfold       = clus->fJustUnfold;
892     fFillAODFile      = clus->fFillAODFile;
893     fRecoUtils        = clus->fRecoUtils; 
894     fConfigName       = clus->fConfigName;
895     fMaxEvent         = clus->fMaxEvent;
896     fDoTrackMatching  = clus->fDoTrackMatching;
897     fOutputAODBranchName = clus->fOutputAODBranchName;
898     for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
899   }
900   
901   // Init geometry, I do not like much to do it like this ...
902   if(fImportGeometryFromFile && !gGeoManager) 
903   {
904     if (fImportGeometryFilePath == "") fImportGeometryFilePath = "$ALICE_ROOT/PWGGA/EMCALTasks/macros/geometry.root" ; // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
905     printf("AliAnalysisTaskEMCALClusterize::Init() - Import %s\n",fImportGeometryFilePath.Data());
906     TGeoManager::Import(fImportGeometryFilePath) ; 
907   }
908
909 }  
910
911 //_______________________________________________________
912 void AliAnalysisTaskEMCALClusterize::InitClusterization()
913 {
914   //Select clusterization/unfolding algorithm and set all the needed parameters
915   
916   if (fJustUnfold){
917     // init the unfolding afterburner 
918     delete fUnfolder;
919     fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
920     return;
921   }
922   
923   //First init the clusterizer
924   delete fClusterizer;
925   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
926     fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
927   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2) 
928     fClusterizer = new AliEMCALClusterizerv2(fGeom, fCalibData, fPedestalData);
929   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN){ 
930     fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
931     fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
932     fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
933   } else {
934     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
935   }
936   
937   //Now set the parameters
938   fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
939   fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
940   fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );    
941   fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
942   fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
943   fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
944   fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
945   fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
946   fClusterizer->SetInputCalibrated       ( kTRUE                               );
947   fClusterizer->SetJustClusters          ( kTRUE                               );  
948   //In case of unfolding after clusterization is requested, set the corresponding parameters
949   if(fRecParam->GetUnfold()){
950     Int_t i=0;
951     for (i = 0; i < 8; i++) {
952       fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
953     }//end of loop over parameters
954     for (i = 0; i < 3; i++) {
955       fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
956       fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
957     }//end of loop over parameters
958     
959     fClusterizer->InitClusterUnfolding();
960     
961   }// to unfold
962 }
963
964 //_________________________________________________
965 void AliAnalysisTaskEMCALClusterize::InitGeometry()
966 {
967   // Init geometry and set the geometry matrix, for the first event, skip the rest
968   // Also set once the run dependent calibrations
969   
970   if(fGeomMatrixSet) return;
971   
972   Int_t runnumber = InputEvent()->GetRunNumber() ;
973   if (!fGeom)
974   {
975     if(fGeomName=="")
976     {
977       if     (runnumber < 140000) fGeomName = "EMCAL_FIRSTYEARV1";
978       else if(runnumber < 171000) fGeomName = "EMCAL_COMPLETEV1";
979       else                        fGeomName = "EMCAL_COMPLETE12SMV1";  
980       printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Set EMCAL geometry name to <%s> for run %d\n",fGeomName.Data(),runnumber);
981     }
982     
983                 fGeom = AliEMCALGeometry::GetInstance(fGeomName);
984     
985                 if(fDebug > 0)
986     {
987                         printf("AliAnalysisTaskEMCALClusterize::InitGeometry(run=%d)",runnumber);
988                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
989                         printf("\n");
990                 }
991         } // geometry pointer did not exist before
992   
993   if(fLoadGeomMatrices)
994   {
995     printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Load user defined EMCAL geometry matrices\n");
996     
997     // OADB if available
998     AliOADBContainer emcGeoMat("AliEMCALgeo");
999     emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
1000     TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
1001     
1002     
1003     for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1004     {
1005       
1006       if (!fGeomMatrix[mod]) // Get it from OADB
1007       {
1008         if(fDebug > 1 ) 
1009           printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - EMCAL matrices SM %d, %p\n",
1010                  mod,((TGeoHMatrix*) matEMCAL->At(mod)));
1011         //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
1012         
1013         fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
1014       }        
1015       
1016       if(fGeomMatrix[mod])
1017       {
1018         if(DebugLevel() > 1) 
1019           fGeomMatrix[mod]->Print();
1020         
1021         fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;  
1022       }
1023       
1024       fGeomMatrixSet=kTRUE;
1025       
1026     }//SM loop
1027   }//Load matrices
1028   else if(!gGeoManager)
1029   {
1030     printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - Get geo matrices from data");
1031     //Still not implemented in AOD, just a workaround to be able to work at least with ESDs     
1032     if(!strcmp(fEvent->GetName(),"AliAODEvent")) 
1033     {
1034       if(DebugLevel() > 1) 
1035         Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
1036     }//AOD
1037     else 
1038     {   
1039       if(DebugLevel() > 1) 
1040         printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
1041       
1042       AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ;
1043       
1044       if(!esd)
1045       {
1046         Error("InitGeometry"," - This event does not contain ESDs?");
1047         AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
1048         return;
1049       }
1050       
1051       for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1052       {
1053         if(DebugLevel() > 1) 
1054           esd->GetEMCALMatrix(mod)->Print();
1055         
1056         if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
1057         
1058       } 
1059       
1060       fGeomMatrixSet=kTRUE;
1061       
1062     }//ESD
1063   }//Load matrices from Data 
1064   
1065   
1066 }
1067
1068 //____________________________________________________
1069 Bool_t AliAnalysisTaskEMCALClusterize::IsExoticEvent()
1070 {
1071   
1072   // Check if event is exotic, get an exotic cell and compare with triggered patch
1073   // If there is a match remove event ... to be completed, filled with something provisional
1074   
1075   if(!fRemoveExoticEvents) return kFALSE;
1076   
1077   // Loop on cells
1078   AliVCaloCells * cells = fEvent->GetEMCALCells();
1079   Float_t totCellE = 0;
1080   Int_t bc = InputEvent()->GetBunchCrossNumber();
1081   for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++){
1082     
1083     Float_t  ecell = 0 ;
1084     Double_t tcell = 0 ;
1085     
1086     Int_t absID   = cells->GetCellNumber(icell);
1087     Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,cells);
1088     if(accept && !fRecoUtils->IsExoticCell(absID,cells,bc)) totCellE += ecell;
1089   }
1090   
1091   //  TString triggerclasses = "";
1092   //  if(esdevent) triggerclasses = esdevent             ->GetFiredTriggerClasses();
1093   //  else         triggerclasses = ((AliAODEvent*)event)->GetFiredTriggerClasses();
1094   //    //  
1095   //    printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster  - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
1096   //    AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
1097   //    return;
1098   //  
1099   
1100   //printf("TotE cell %f\n",totCellE);
1101   if(totCellE < 1) return kTRUE;
1102   
1103   return kFALSE;
1104   
1105
1106
1107 //________________________________________________________________
1108 Bool_t AliAnalysisTaskEMCALClusterize::IsLEDEvent(const Int_t run)
1109 {
1110   //Check if event is LED
1111   
1112   if(!fRemoveLEDEvents) return kFALSE;
1113   
1114   //check events of LHC11a period
1115   if(run < 140000  || run > 146860) return kFALSE ;
1116   
1117   // Count number of cells with energy larger than 0.1 in SM3, cut on this number
1118   Int_t ncellsSM3 = 0;
1119   Int_t ncellsSM4 = 0;
1120   AliVCaloCells * cells = fEvent->GetEMCALCells();
1121   for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1122   {
1123     if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
1124     if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==4) ncellsSM4++;      
1125   }
1126   
1127   TString triggerclasses = "";
1128   
1129   AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(fEvent);
1130   if(esdevent) triggerclasses = esdevent              ->GetFiredTriggerClasses();
1131   else         triggerclasses = ((AliAODEvent*)fEvent)->GetFiredTriggerClasses();
1132   
1133   Int_t ncellcut = 21;
1134   if(triggerclasses.Contains("EMC")) ncellcut = 35;
1135   
1136   if( ncellsSM3 >= ncellcut || ncellsSM4 >= 100 )
1137   {
1138     printf("AliAnalysisTaksEMCALClusterize::IsLEDEvent() - reject event %d with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
1139     AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
1140     return kTRUE;
1141   }
1142   
1143   return kFALSE;
1144   
1145
1146
1147 //______________________________________________________________________________
1148 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, 
1149                                                         TObjArray *recPoints, 
1150                                                         TObjArray *clusArray)
1151 {
1152   // Restore clusters from recPoints
1153   // Cluster energy, global position, cells and their amplitude fractions are restored
1154   Int_t j = 0;
1155   for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++)
1156   {
1157     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i);
1158     
1159     const Int_t ncells = recPoint->GetMultiplicity();
1160     Int_t ncellsTrue = 0;
1161     
1162     if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
1163     
1164     // cells and their amplitude fractions
1165     UShort_t   absIds[ncells];  
1166     Double32_t ratios[ncells];
1167     
1168     //For later check embedding
1169     AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1170     
1171     Float_t clusterE = 0; 
1172     for (Int_t c = 0; c < ncells; c++) 
1173     {
1174       AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]);
1175       
1176       absIds[ncellsTrue] = digit->GetId();
1177       ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
1178       
1179       // In case of unfolding, remove digits with unfolded energy too low      
1180       if(fSelectCell){
1181         if     (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)  
1182         {
1183           if(DebugLevel() > 1)  
1184           {
1185             printf("AliAnalysisTaksEMCALClusterize::RecPoints2Clusters() - Too small energy in cell of cluster: cluster cell %f, digit %f\n",
1186                    recPoint->GetEnergiesList()[c],digit->GetAmplitude());
1187           }
1188           
1189           continue;
1190           
1191         } // if cuts
1192       }// Select cells
1193       
1194       //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
1195       ncellsTrue++;
1196       clusterE  +=recPoint->GetEnergiesList()[c];
1197       
1198       // In case of embedding, fill ratio with amount of signal, 
1199       if (aodIH && aodIH->GetMergeEvents()) 
1200       {
1201         //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
1202         AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
1203         AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
1204         
1205         Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1206         //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1207         Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1208         //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
1209         
1210         if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
1211         //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
1212         
1213       }//Embedding
1214       
1215     }// cluster cell loop
1216     
1217     if (ncellsTrue < 1) 
1218     {
1219       if (DebugLevel() > 1) 
1220         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Skipping cluster with no cells avobe threshold E = %f, ncells %d\n",
1221                recPoint->GetEnergy(), ncells);
1222       continue;
1223     }
1224     
1225     //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
1226     
1227     if(clusterE <  fRecParam->GetClusteringThreshold()) 
1228     {
1229       if (DebugLevel()>1)
1230         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Remove cluster with energy below seed threshold %f\n",clusterE);
1231       continue;
1232     }
1233     
1234     TVector3 gpos;
1235     Float_t g[3];
1236     
1237     // calculate new cluster position
1238     recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr);
1239     recPoint->GetGlobalPosition(gpos);
1240     gpos.GetXYZ(g);
1241     
1242     // create a new cluster
1243     (*clusArray)[j] = new AliAODCaloCluster() ;
1244     AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( clusArray->At(j) ) ;
1245     j++;
1246     clus->SetType(AliVCluster::kEMCALClusterv1);
1247     clus->SetE(clusterE);
1248     clus->SetPosition(g);
1249     clus->SetNCells(ncellsTrue);
1250     clus->SetCellsAbsId(absIds);
1251     clus->SetCellsAmplitudeFraction(ratios);
1252     clus->SetChi2(-1); //not yet implemented
1253     clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
1254     clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
1255     clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower()); 
1256
1257     if(ncells == ncellsTrue)
1258     {
1259       Float_t elipAxis[2];
1260       recPoint->GetElipsAxis(elipAxis);
1261       clus->SetM02(elipAxis[0]*elipAxis[0]) ;
1262       clus->SetM20(elipAxis[1]*elipAxis[1]) ;
1263       clus->SetDispersion(recPoint->GetDispersion());
1264     }
1265     else if(fSelectCell)
1266     {
1267       // In case some cells rejected, in unfolding case, recalculate
1268       // shower shape parameters and position
1269       if(DebugLevel() > 1) 
1270         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape\n",ncells,ncellsTrue);
1271       AliVCaloCells* cells = 0x0; 
1272       if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent()  ->GetEMCALCells();
1273       else                                  cells = InputEvent()->GetEMCALCells();
1274       fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus);
1275       fRecoUtils->RecalculateClusterPID(clus);
1276       fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus); 
1277       
1278     }
1279     
1280     //MC
1281     Int_t  parentMult = 0;
1282     Int_t *parentList = recPoint->GetParents(parentMult);
1283     clus->SetLabel(parentList, parentMult); 
1284     
1285     //Write the second major contributor to each MC cluster.
1286     Int_t iNewLabel ;
1287     for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1288     {
1289       
1290       Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1291       if(idCell>=0){
1292         iNewLabel = 1 ; //iNewLabel makes sure we  don't write twice the same label.
1293         for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
1294         {
1295           if ( fCellSecondLabels[idCell] == -1 )  iNewLabel = 0;  // -1 is never a good second label.
1296           if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) )  iNewLabel = 0;
1297         }
1298         if (iNewLabel == 1) 
1299         {
1300           Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
1301           for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ )
1302           {
1303             newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
1304           }
1305           
1306           newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
1307           clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
1308           delete [] newLabelArray;
1309         }
1310       }//positive cell number
1311     }
1312     
1313   } // recPoints loop
1314   
1315 }
1316
1317
1318 //____________________________________________________________
1319 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
1320 {
1321   // Init geometry, create list of output clusters
1322   
1323   if(fOutputAODBranchName.Length()!=0)
1324   {
1325     fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1326     fOutputAODBranch->SetName(fOutputAODBranchName);
1327     //fOutputAODBranch->SetOwner(kFALSE);
1328     AddAODBranch("TClonesArray", &fOutputAODBranch);
1329   }
1330   else
1331   {
1332     AliFatal("fOutputAODBranchName not set\n");
1333   }
1334   
1335   //PostData(0,fOutputAODBranch);
1336   
1337 }
1338
1339 //_______________________________________________________
1340 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *) 
1341 {
1342   // Main loop
1343   // Called for each event
1344   
1345   //Remove the contents of output list set in the previous event 
1346   fOutputAODBranch->Clear("C");
1347   
1348   LoadBranches();
1349   
1350   //Get the event, do some checks and settings
1351   CheckAndGetEvent() ;
1352   
1353   if (!fEvent) 
1354   {
1355     if(DebugLevel() > 0 ) printf("AliAnalysisTaksEMCALClusterize::UserExec() - Skip Event %d", (Int_t) Entry());
1356     return ;
1357   }
1358   
1359   //Init pointers, geometry, clusterizer, ocdb, aodb
1360   
1361   InitGeometry(); // only once, must be done before OADB, geo OADB accessed here
1362   
1363   if(fAccessOCDB) AccessOCDB();
1364   if(fAccessOADB) AccessOADB(); // only once
1365
1366   InitClusterization();
1367   
1368   // Make clusters
1369   if (fJustUnfold) ClusterUnfolding();
1370   else             ClusterizeCells() ;
1371     
1372   //Recalculate track-matching for the new clusters
1373   if(fDoTrackMatching) 
1374   {
1375     fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom);
1376   }
1377   //Put the new clusters in the AOD list
1378   
1379   Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntriesFast();
1380   for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
1381     AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
1382     
1383     newCluster->SetID(i);
1384
1385     //Add matched track
1386     if(fDoTrackMatching)
1387     {
1388       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
1389       if(trackIndex >= 0)
1390       {
1391         newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex));
1392         if(DebugLevel() > 1) 
1393           printf("AliAnalysisTaksEMCALClusterize::UserExec() - Matched Track index %d to new cluster %d \n",trackIndex,i);
1394       }
1395       
1396       Float_t dR = 999., dZ = 999.;
1397       fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dR,dZ);
1398       newCluster->SetTrackDistance(dR,dZ);
1399       
1400     }
1401     else 
1402     {// Assign previously assigned matched track in reco, very very rough
1403       Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster
1404       newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]);
1405     }
1406
1407     //printf("New cluster E %f, Time  %e, Id = ", newCluster->E(), newCluster->GetTOF() );
1408     //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] );
1409     //printf("\n");
1410
1411     // Calculate distance to bad channel for new cluster. Make sure you give the list of bad channels.
1412     fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fEvent->GetEMCALCells(), newCluster);
1413     
1414     new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
1415         
1416     if(DebugLevel() > 1 )    
1417       printf("AliAnalysisTaksEMCALClusterize::UserExec() - New cluster %d of %d, energy %f\n",newCluster->GetID(), kNumberOfCaloClusters, newCluster->E());
1418     
1419   } // cluster loop
1420   
1421   fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
1422   
1423   // Clean up
1424   fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
1425   
1426 }      
1427