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