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