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