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