correct compilation warning
[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / 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 /* $Id: $ */
16
17 //_________________________________________________________________________
18 // This analysis provides a new list of clusters to be used in other analysis
19 //
20 // Author: Gustavo Conesa Balbastre,
21 //         Adapted from analysis class from Deepa Thomas
22 //
23 //
24 //_________________________________________________________________________
25
26 // --- Root ---
27 #include "TString.h"
28 #include "TRefArray.h"
29 #include "TClonesArray.h"
30 #include "TTree.h"
31 #include "TGeoManager.h"
32 #include "TROOT.h"
33 #include "TInterpreter.h"
34 #include "TFile.h"
35
36 // --- AliRoot Analysis Steering
37 #include "AliAnalysisTask.h"
38 #include "AliAnalysisManager.h"
39 #include "AliESDEvent.h"
40 #include "AliGeomManager.h"
41 #include "AliVCaloCells.h"
42 #include "AliAODCaloCluster.h"
43 #include "AliCDBManager.h"
44 #include "AliCDBStorage.h"
45 #include "AliCDBEntry.h"
46 #include "AliLog.h"
47 #include "AliVEventHandler.h"
48 #include "AliAODInputHandler.h"
49
50 // --- EMCAL
51 #include "AliEMCALRecParam.h"
52 #include "AliEMCALAfterBurnerUF.h"
53 #include "AliEMCALGeometry.h"
54 #include "AliEMCALClusterizerNxN.h"
55 #include "AliEMCALClusterizerv1.h"
56 #include "AliEMCALRecPoint.h"
57 #include "AliEMCALDigit.h"
58 #include "AliCaloCalibPedestal.h"
59 #include "AliEMCALCalibData.h"
60 #include "AliEMCALRecoUtils.h"
61
62 #include "AliAnalysisTaskEMCALClusterize.h"
63
64 ClassImp(AliAnalysisTaskEMCALClusterize)
65
66 //________________________________________________________________________
67 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name) 
68   : AliAnalysisTaskSE(name)
69   , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"),  fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
70   , fCalibData(0),       fPedestalData(0),     fOCDBpath("raw://")
71   , fDigitsArr(0),       fClusterArr(0),       fCaloClusterArr(0)
72   , fRecParam(0),        fClusterizer(0),      fUnfolder(0),           fJustUnfold(kFALSE) 
73   , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters")
74   , fFillAODFile(kTRUE), fFillAODHeader(0),    fFillAODCaloCells(0)
75   , fRun(-1),            fRecoUtils(0),        fConfigName("")
76   , fCellLabels(),       fCellSecondLabels(),  fMaxEvent(1000000000)
77   
78   {
79   //ctor
80   for(Int_t i = 0; i < 10;    i++)  fGeomMatrix[i] =  0;
81     for(Int_t j = 0; j < 24*48*11; j++)  {
82       fCellLabels[j]       = -1;
83       fCellSecondLabels[j] = -1;
84     }  
85   fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
86   fClusterArr      = new TObjArray(100);
87   fCaloClusterArr  = new TObjArray(100);
88   fRecParam        = new AliEMCALRecParam;
89   fBranchNames     = "ESD:AliESDHeader.,EMCALCells.";
90   fRecoUtils       = new AliEMCALRecoUtils();
91 }
92
93 //________________________________________________________________________
94 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize() 
95   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
96   , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"),   fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
97   , fCalibData(0),        fPedestalData(0),     fOCDBpath("raw://")
98   , fDigitsArr(0),        fClusterArr(0),       fCaloClusterArr(0)
99   , fRecParam(0),         fClusterizer(0),      fUnfolder(0),           fJustUnfold(kFALSE)
100   , fOutputAODBranch(0),  fOutputAODBranchName("newEMCALClusters")
101   , fFillAODFile(kFALSE), fFillAODHeader(0),    fFillAODCaloCells(0)
102   , fRun(-1),             fRecoUtils(0),        fConfigName("") 
103   , fCellLabels(),        fCellSecondLabels(),  fMaxEvent(1000000000)
104 {
105   // Constructor
106   for(Int_t i = 0; i < 10;    i++)  fGeomMatrix[i] =  0;
107   for(Int_t j = 0; j < 24*48*11; j++)  {
108     fCellLabels[j]       = -1;
109     fCellSecondLabels[j] = -1;
110   }
111   fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
112   fClusterArr      = new TObjArray(100);
113   fCaloClusterArr  = new TObjArray(100);
114   fRecParam        = new AliEMCALRecParam;
115   fBranchNames     = "ESD:AliESDHeader.,EMCALCells.";
116   fRecoUtils       = new AliEMCALRecoUtils();
117 }
118
119
120 //________________________________________________________________________
121 AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
122 {
123   //dtor 
124   
125   if (fDigitsArr){
126     fDigitsArr->Clear("C");
127     delete fDigitsArr; 
128   }
129   
130   if (fClusterArr){
131     fClusterArr->Delete();
132     delete fClusterArr;
133   }
134   
135   if (fCaloClusterArr){
136     fCaloClusterArr->Delete();
137     delete fCaloClusterArr; 
138   }
139
140   if(fClusterizer) delete fClusterizer;
141   if(fUnfolder)    delete fUnfolder;   
142   if(fRecoUtils)   delete fRecoUtils;
143
144 }
145
146 //-------------------------------------------------------------------
147 void AliAnalysisTaskEMCALClusterize::Init()
148 {
149   //Init analysis with configuration macro if available
150   
151   if(gROOT->LoadMacro(fConfigName) >=0){
152     printf("Configure analysis with %s\n",fConfigName.Data());
153     AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
154     fGeomName         = clus->fGeomName; 
155     fLoadGeomMatrices = clus->fLoadGeomMatrices;
156     fOCDBpath         = clus->fOCDBpath;
157     fRecParam         = clus->fRecParam;
158     fJustUnfold       = clus->fJustUnfold;
159     fFillAODFile      = clus->fFillAODFile;
160     fRecoUtils        = clus->fRecoUtils; 
161     fConfigName       = clus->fConfigName;
162     fOutputAODBranchName = clus->fOutputAODBranchName;
163     for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
164
165   }
166
167 }  
168
169 //-------------------------------------------------------------------
170 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
171 {
172   // Init geometry, create list of output clusters
173
174   fGeom =  AliEMCALGeometry::GetInstance(fGeomName) ;   
175   if(fOutputAODBranchName.Length()!=0){
176     fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
177     fOutputAODBranch->SetName(fOutputAODBranchName);
178     AddAODBranch("TClonesArray", &fOutputAODBranch);
179   }
180   else {
181     AliFatal("fOutputAODBranchName not set\n");
182   }
183 }
184
185 //________________________________________________________________________
186 void AliAnalysisTaskEMCALClusterize::FillAODCaloCells() {
187   // Put calo cells in standard branch  
188   AliVEvent * event = InputEvent();
189   AliVCaloCells &eventEMcells = *(event->GetEMCALCells());
190   Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
191   
192   AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
193   aodEMcells.CreateContainer(nEMcell);
194   aodEMcells.SetType(AliVCaloCells::kEMCALCell);
195   Double_t calibFactor = 1.;   
196   for (Int_t iCell = 0; iCell < nEMcell; iCell++) { 
197     Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
198     fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta); 
199     fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);     
200     
201     if(fRecoUtils->IsRecalibrationOn()){ 
202       calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
203     }
204     
205     if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ //Channel is not declared as bad
206       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor);
207       //printf("GOOD channel\n");
208     }
209     else {
210       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0);
211       //printf("BAD channel\n");
212     }
213   }
214   aodEMcells.Sort();
215   
216 }
217
218 //________________________________________________________________________
219 void AliAnalysisTaskEMCALClusterize::FillAODHeader() {
220   //Put event header information in standard AOD branch
221   
222   AliVEvent* event      = InputEvent();
223   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
224   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
225
226   Double_t pos[3]   ;
227   Double_t covVtx[6];
228   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
229   
230   AliAODHeader* header = AODEvent()->GetHeader();
231   header->SetRunNumber(event->GetRunNumber());
232   
233   if(esdevent){
234     TTree* tree = fInputHandler->GetTree();
235     if (tree) {
236       TFile* file = tree->GetCurrentFile();
237       if (file) header->SetESDFileName(file->GetName());
238     }
239   }
240   else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName());
241   
242   header->SetBunchCrossNumber(event->GetBunchCrossNumber());
243   header->SetOrbitNumber(event->GetOrbitNumber());
244   header->SetPeriodNumber(event->GetPeriodNumber());
245   header->SetEventType(event->GetEventType());
246   
247   //Centrality
248   if(event->GetCentrality()){
249     header->SetCentrality(new AliCentrality(*(event->GetCentrality())));
250   }
251   else{
252     header->SetCentrality(0);
253   }
254   
255   //Trigger  
256   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
257   if      (esdevent) header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses());
258   else if (aodevent) header->SetFiredTriggerClasses(aodevent->GetFiredTriggerClasses());
259   header->SetTriggerMask(event->GetTriggerMask()); 
260   header->SetTriggerCluster(event->GetTriggerCluster());
261   if(esdevent){
262     header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());    
263     header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());    
264     header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());    
265   }
266   else if (aodevent){
267     header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());    
268     header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());    
269     header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());    
270   }
271   
272   header->SetMagneticField(event->GetMagneticField());
273   //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.); 
274   
275   header->SetZDCN1Energy(event->GetZDCN1Energy());
276   header->SetZDCP1Energy(event->GetZDCP1Energy());
277   header->SetZDCN2Energy(event->GetZDCN2Energy());
278   header->SetZDCP2Energy(event->GetZDCP2Energy());
279   header->SetZDCEMEnergy(event->GetZDCEMEnergy(0),event->GetZDCEMEnergy(1));
280   
281   Float_t diamxy[2]={event->GetDiamondX(),event->GetDiamondY()};
282   Float_t diamcov[3];
283   event->GetDiamondCovXY(diamcov);
284   header->SetDiamond(diamxy,diamcov);
285   if      (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
286   else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
287   
288   //
289   //
290   Int_t nVertices = 1 ;/* = prim. vtx*/;
291   Int_t nCaloClus = event->GetNumberOfCaloClusters();
292   
293   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
294   
295   // Access to the AOD container of vertices
296   TClonesArray &vertices = *(AODEvent()->GetVertices());
297   Int_t jVertices=0;
298   
299   // Add primary vertex. The primary tracks will be defined
300   // after the loops on the composite objects (V0, cascades, kinks)
301   event->GetPrimaryVertex()->GetXYZ(pos);
302   Float_t chi = 0;
303   if      (esdevent){
304     esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
305     chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
306   }
307   else if (aodevent){
308     aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
309     chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
310   }
311   
312   AliAODVertex * primary = new(vertices[jVertices++])
313   AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
314   primary->SetName(event->GetPrimaryVertex()->GetName());
315   primary->SetTitle(event->GetPrimaryVertex()->GetTitle());
316   
317 }
318
319 //________________________________________________________________________
320 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *) 
321 {
322   // Main loop
323   // Called for each event
324   
325   AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
326   Int_t eventN = Entry();
327   if(aodIH) 
328     eventN = aodIH->GetReadEntry();
329   //printf("Clusterizer --- Event %d-- \n",eventN);
330
331   if (eventN > fMaxEvent) return;
332
333   //Remove the contents of output list set in the previous event 
334   fOutputAODBranch->Clear("C");
335   
336 //  if(fOutputAODBranchName.Length()!=0){
337 //    //Remove the contents of output list set in the previous event 
338 //    fOutputAODBranch->Clear("C");
339 //  }
340 //  else{ 
341 //    //Reset and put new clusters in standard branch, also cells
342 //    AODEvent()->ResetStd(0, 1, 0, 0, 0, AODEvent()->GetNumberOfCaloClusters(), 0, 0);
343 //    fOutputAODBranch = AODEvent()->GetCaloClusters();// Put new clusters in standard branch
344 //  }
345   
346   //Magic line to write events to AOD file
347   AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
348   LoadBranches();
349   
350   AccessOCDB();
351   
352   //Get the event
353   AliVEvent   * event    = 0;
354   AliESDEvent * esdevent = 0;
355     
356   //Fill output event with headerø
357
358   //Check if input event are embedded events
359   //If so, take output event
360   if (aodIH && aodIH->GetMergeEvents()) {
361     //printf("AliAnalysisTaskEMCALClusterize::UserExec() - Use embedded events\øn");
362     event  = AODEvent();
363     
364     if(!aodIH->GetMergeEMCALCells()) 
365       AliFatal("Events merged but not EMCAL cells, check analysis settings!");
366     
367 //    printf("InputEvent  N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
368 //           InputEvent()->GetEMCALCells()->GetNumberOfCells());
369 //    printf("MergedEvent  N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(),
370 //           aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells());
371 //    printf("OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
372 //           AODEvent()->GetEMCALCells()->GetNumberOfCells());
373     
374   }
375   else {
376     event =  InputEvent();
377     esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
378     if(fFillAODCaloCells) FillAODCaloCells();   
379     if(fFillAODHeader)    FillAODHeader();
380   }
381   
382   if (!event) {
383     Error("UserExec","Event not available");
384     return;
385   }
386
387   //-------------------------------------------------------------------------------------
388   //Set the geometry matrix, for the first event, skip the rest
389   //-------------------------------------------------------------------------------------
390   if(!fGeomMatrixSet){
391     if(fLoadGeomMatrices){
392       for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
393         if(fGeomMatrix[mod]){
394           if(DebugLevel() > 1) 
395             fGeomMatrix[mod]->Print();
396           fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;  
397         }
398         fGeomMatrixSet=kTRUE;
399       }//SM loop
400     }//Load matrices
401     else if(!gGeoManager){
402       Info("UserExec","Get geo matrices from data");
403       //Still not implemented in AOD, just a workaround to be able to work at least with ESDs   
404       if(!strcmp(event->GetName(),"AliAODEvent")) {
405         if(DebugLevel() > 1) 
406           Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
407       }//AOD
408       else {    
409         if(DebugLevel() > 1) 
410           Info("UserExec","AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
411         AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
412         if(!esd) {
413           Error("UserExec","This event does not contain ESDs?");
414           return;
415         }
416         for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
417           if(DebugLevel() > 1) 
418             esd->GetEMCALMatrix(mod)->Print();
419           if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
420         } 
421         fGeomMatrixSet=kTRUE;
422       }//ESD
423     }//Load matrices from Data 
424     
425     //Recover time dependent corrections, put then in recalibration histograms. Do it once
426     fRecoUtils->SetTimeDependentCorrections(InputEvent()->GetRunNumber());
427
428   }//first event
429
430   //Get the list of cells needed for unfolding and reclustering
431   AliVCaloCells *cells= event->GetEMCALCells();
432   Int_t nClustersOrg = 0;
433   //-------------------------------------------
434   //---------Unfolding clusters----------------
435   //-------------------------------------------
436   if (fJustUnfold) {
437     
438     //Fill the array with the EMCAL clusters, copy them
439     for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
440     {
441       AliVCluster *clus = event->GetCaloCluster(i);
442       if(clus->IsEMCAL()){        
443         //printf("Org Cluster %d, E %f\n",i, clus->E());
444         
445         //recalibrate/remove bad channels/etc if requested
446         if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells())){
447           //printf("Remove cluster\n");
448           continue;
449         } 
450         
451         if(fRecoUtils->IsRecalibrationOn()){
452           //printf("Energy before %f ",clus->E());
453           fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
454           //printf("; after %f\n",clus->E());
455         }
456         //Cast to ESD or AOD, needed to create the cluster array
457         AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
458         AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
459         if     (esdCluster){
460           fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );   
461         }//ESD
462         else if(aodCluster){
463           fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );   
464         }//AOD
465         else 
466           Warning("UserExec()"," - Wrong CaloCluster type?");
467         nClustersOrg++;
468       }
469     }
470     
471     //Do the unfolding
472     fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
473     
474     //CLEAN-UP
475     fUnfolder->Clear();
476     
477     
478     //printf("nClustersOrg %d\n",nClustersOrg);
479   }
480   //-------------------------------------------
481   //---------- Recluster cells ----------------
482   //-------------------------------------------
483   
484   else{
485     //-------------------------------------------------------------------------------------
486     //Transform CaloCells into Digits
487     //-------------------------------------------------------------------------------------
488     
489     //In case of MC, first loop on the clusters and fill MC label to array
490     //.....................................................................
491     
492 //    for(Int_t j = 0; j < 24*48*11; j++)  {
493 //      if(fCellLabels[j]!=-1) printf("Not well initialized cell %d, label %d\n",j,fCellLabels[j]) ;
494 //    }
495     
496     Int_t nClusters = event->GetNumberOfCaloClusters();
497     if(aodIH && aodIH->GetEventToMerge())  //Embedding
498       nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
499     for (Int_t i = 0; i < nClusters; i++)
500     {
501       AliVCluster *clus = 0;
502       if(aodIH && aodIH->GetEventToMerge()) //Embedding
503         clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
504       else      
505         clus = event->GetCaloCluster(i);
506
507       if(!clus) {
508         printf("AliEMCALReclusterize::UserExec() - No Clusters\n");
509         return;
510       }
511       
512       if(clus->IsEMCAL()){   
513         //printf("Cluster Signal %d, energy %f\n",clus->GetID(),clus->E());
514         Int_t label = clus->GetLabel();
515         Int_t label2 = -1 ;
516         if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
517         UShort_t * index    = clus->GetCellsAbsId() ;
518         for(Int_t icell=0; icell < clus->GetNCells(); icell++ ){
519           fCellLabels[index[icell]]=label;
520           fCellSecondLabels[index[icell]]=label2;
521           //printf("1) absID %d, label[0] %d label[1] %d\n",index[icell], fCellLabels[index[icell]],fCellSecondLabels[index[icell]]);
522         }
523         nClustersOrg++;
524       }
525     } 
526
527     // Create digits 
528     //.................
529     Int_t idigit =  0;
530     Int_t id     = -1;
531     Float_t amp  = -1; 
532     Float_t time = -1; 
533     
534     Double_t cellAmplitude = 0;
535     Double_t cellTime      = 0;
536     Short_t  cellNumber    = 0;
537         
538     TTree *digitsTree = new TTree("digitstree","digitstree");
539     digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
540     
541     
542     for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
543     {
544       if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
545         break;
546       
547       time = cellTime;
548       amp  = cellAmplitude;
549       id   = cellNumber;
550       
551       Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
552       fGeom->GetCellIndex(id,imod,iTower,iIphi,iIeta); 
553       fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);   
554       
555       //Do not include bad channels found in analysis?
556       if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() && 
557           fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
558           fCellLabels[id]=-1; //reset the entry in the array for next event
559         //printf("Remove channel %d\n",id);
560         continue;
561       }
562              
563       //Recalibrate?
564       if(fRecoUtils->IsRecalibrationOn()){ 
565         //printf("CalibFactor %f times %f for id %d\n",fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),amp,id);
566         amp *=fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
567       }
568            
569       //Create the digit, put a fake primary deposited energy to trick the clusterizer when checking the most likely primary
570       new((*fDigitsArr)[idigit]) AliEMCALDigit( fCellLabels[id], fCellLabels[id],id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, 1); 
571       //if(fCellLabels[id]>=0)printf("2) Digit cell %d, label %d\n",id,fCellLabels[id]) ;
572       //else                  printf("2) Digit cell %d, no label, amp %f \n",id,amp) ;
573       fCellLabels[id]=-1; //reset the entry in the array for next event
574       
575       //AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
576       //digit->SetId(id);
577       //digit->SetAmplitude(amp);
578       //digit->SetTime(time);
579       //digit->SetTimeR(time);
580       //digit->SetIndexInList(idigit);
581       //digit->SetType(AliEMCALDigit::kHG);
582       
583       //printf("Digit: Id %d, amp %f, time %e, index %d\n",id, amp,time,idigit);
584       idigit++;
585     }
586     
587     //Fill the tree with digits
588     digitsTree->Fill();
589     
590     //-------------------------------------------------------------------------------------
591     //Do the clusterization
592     //-------------------------------------------------------------------------------------        
593     TTree *clustersTree = new TTree("clustertree","clustertree");
594
595     fClusterizer->SetInput(digitsTree);
596     fClusterizer->SetOutput(clustersTree);
597     fClusterizer->Digits2Clusters("");
598
599     //-------------------------------------------------------------------------------------
600     //Transform the recpoints into AliVClusters
601     //-------------------------------------------------------------------------------------
602     
603     clustersTree->SetBranchStatus("*",0); //disable all branches
604     clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
605     
606     TBranch *branch = clustersTree->GetBranch("EMCALECARP");
607     branch->SetAddress(&fClusterArr);
608     branch->GetEntry(0);
609
610     RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
611     
612     //---CLEAN UP-----
613     fClusterizer->Clear();
614     fDigitsArr  ->Clear("C");
615     fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
616
617     clustersTree->Delete("all");
618     digitsTree  ->Delete("all");
619   }
620   
621   //Recalculate track-matching for the new clusters, only with ESDs
622   if(esdevent)fRecoUtils->FindMatches(esdevent,fCaloClusterArr);
623
624   
625   //-------------------------------------------------------------------------------------
626   //Put the new clusters in the AOD list
627   //-------------------------------------------------------------------------------------
628   
629   Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntries();
630   //printf("New clusters %d, Org clusters %d\n",kNumberOfCaloClusters, nClustersOrg);
631   for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
632     AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
633     //if(Entry()==0) Info("UserExec","newCluster E %f\n", newCluster->E());
634     
635     //Add matched track, if any, only with ESDs
636     if(esdevent){
637       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
638       if(trackIndex >= 0){
639         newCluster->AddTrackMatched(event->GetTrack(trackIndex));
640         if(DebugLevel() > 1) 
641           Info("UserExec","Matched Track index %d to new cluster %d \n",trackIndex,i);
642       }
643     }
644     
645     //In case of new bad channels, recalculate distance to bad channels
646     if(fRecoUtils->IsBadChannelsRemovalSwitchedOn()){
647       //printf("DistToBad before %f ",newCluster->GetDistanceToBadChannel());
648       fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, event->GetEMCALCells(), newCluster);
649       //printf("; after %f \n",newCluster->GetDistanceToBadChannel());
650     }
651     
652 //    if(newCluster->GetNLabels()>0){
653 //      printf("3) MC: N labels %d, label %d ; ", newCluster->GetNLabels(), newCluster->GetLabel() );
654 //      UShort_t * newindex    = newCluster->GetCellsAbsId() ;
655 //      for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ){
656 //       printf("\t absID %d\n",newindex[icell]);
657 //     }
658 //    }
659     
660 //    printf("Cluster %d, energy %f\n",newCluster->GetID(),newCluster->E());
661 //    printf("labels: ");
662 //    for(UInt_t ilab = 0; ilab < newCluster->GetNLabels(); ilab++)
663 //      printf("Label[%d]: %d ",ilab, newCluster->GetLabelAt(ilab)); 
664 //     printf("\n ");
665     
666     newCluster->SetID(i);
667     new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
668   }
669   
670   //if(fOutputAODBranchName.Length()!=0)
671   fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
672   
673   //---CLEAN UP-----
674   fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
675 }      
676
677 //_____________________________________________________________________
678 Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
679 {
680   //Access to OCDB stuff
681
682   AliVEvent * event = InputEvent();
683   if (!event)
684   {
685     Warning("AccessOCDB","Event not available!!!");
686     return kFALSE;
687   }
688
689   if (event->GetRunNumber()==fRun)
690     return kTRUE;
691   fRun = event->GetRunNumber();
692
693   if(DebugLevel() > 1 )
694     Info("AccessODCD()"," Begin");
695
696   fGeom = AliEMCALGeometry::GetInstance(fGeomName);
697   
698   
699   AliCDBManager *cdb = AliCDBManager::Instance();
700     
701
702   if (fOCDBpath.Length()){
703     cdb->SetDefaultStorage(fOCDBpath.Data());
704     Info("AccessOCDB","Default storage %s",fOCDBpath.Data());
705   }
706   
707   cdb->SetRun(event->GetRunNumber());
708
709   //
710   // EMCAL from RAW OCDB
711   if (fOCDBpath.Contains("alien:"))
712   {
713     cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
714     cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
715   }
716
717   TString path = cdb->GetDefaultStorage()->GetBaseFolder();
718   
719   // init parameters:
720   
721   //Get calibration parameters  
722   if(!fCalibData)
723   {
724     AliCDBEntry *entry = (AliCDBEntry*) 
725       AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
726     
727     if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
728   }
729   
730   if(!fCalibData)
731     AliFatal("Calibration parameters not found in CDB!");
732     
733   //Get calibration parameters  
734   if(!fPedestalData)
735   {
736     AliCDBEntry *entry = (AliCDBEntry*) 
737       AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
738     
739     if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
740   }
741     
742   if(!fPedestalData)
743     AliFatal("Dead map not found in CDB!");
744
745   InitClusterization();
746   
747   return kTRUE;
748 }
749
750 //________________________________________________________________________________________
751 void AliAnalysisTaskEMCALClusterize::InitClusterization()
752 {
753   //Select clusterization/unfolding algorithm and set all the needed parameters
754   
755   if (fJustUnfold){
756     // init the unfolding afterburner 
757     delete fUnfolder;
758     fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
759     return;
760  }
761
762   //First init the clusterizer
763   delete fClusterizer;
764   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
765     fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
766   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) 
767     fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
768   else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
769    AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
770    clusterizer->SetNRowDiff(2);
771    clusterizer->SetNColDiff(2);
772    fClusterizer = clusterizer;
773   } else{
774     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
775   }
776     
777   //Now set the parameters
778   fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
779   fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
780   fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );    
781   fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
782   fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
783   fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
784   fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
785   fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
786   fClusterizer->SetInputCalibrated       ( kTRUE                               );
787   fClusterizer->SetJustClusters          ( kTRUE                               );  
788   //In case of unfolding after clusterization is requested, set the corresponding parameters
789   if(fRecParam->GetUnfold()){
790     Int_t i=0;
791     for (i = 0; i < 8; i++) {
792       fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
793     }//end of loop over parameters
794     for (i = 0; i < 3; i++) {
795       fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
796       fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
797     }//end of loop over parameters
798     
799     fClusterizer->InitClusterUnfolding();
800   }// to unfold
801 }
802
803 //________________________________________________________________________________________
804 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, TObjArray *recPoints, TObjArray *clusArray)
805 {
806   // Restore clusters from recPoints
807   // Cluster energy, global position, cells and their amplitude fractions are restored
808   
809   for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++)
810   {
811     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i);
812     
813     const Int_t ncells = recPoint->GetMultiplicity();
814     Int_t ncells_true = 0;
815     
816     // cells and their amplitude fractions
817     UShort_t   absIds[ncells];  
818     Double32_t ratios[ncells];
819     
820     for (Int_t c = 0; c < ncells; c++) {
821       AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]);
822       
823       absIds[ncells_true] = digit->GetId();
824       ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
825       
826       if (ratios[ncells_true] > 0.001) ncells_true++;
827     }
828     
829     if (ncells_true < 1) {
830       Warning("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters", "skipping cluster with no cells");
831       continue;
832     }
833     
834     TVector3 gpos;
835     Float_t g[3];
836     
837     // calculate new cluster position
838     recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr);
839     recPoint->GetGlobalPosition(gpos);
840     gpos.GetXYZ(g);
841     
842     // create a new cluster
843     AliAODCaloCluster *clus = new AliAODCaloCluster();
844     clus->SetType(AliVCluster::kEMCALClusterv1);
845     clus->SetE(recPoint->GetEnergy());
846     clus->SetPosition(g);
847     clus->SetNCells(ncells_true);
848     clus->SetCellsAbsId(absIds);
849     clus->SetCellsAmplitudeFraction(ratios);
850     clus->SetDispersion(recPoint->GetDispersion());
851     clus->SetChi2(-1); //not yet implemented
852     clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
853     clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
854     Float_t elipAxis[2];
855     recPoint->GetElipsAxis(elipAxis);
856     clus->SetM02(elipAxis[0]*elipAxis[0]) ;
857     clus->SetM20(elipAxis[1]*elipAxis[1]) ;
858     clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower()); 
859     
860     //MC
861     Int_t  parentMult  = 0;
862     Int_t *parentList = recPoint->GetParents(parentMult);
863     clus->SetLabel(parentList, parentMult); 
864     
865     //Write the second major contributor to each MC cluster.
866     Int_t iNewLabel ;
867     for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ ){
868       
869       Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
870       if(idCell>=0){
871         iNewLabel = 1 ; //iNewLabel makes sure we  don't write twice the same label.
872         for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
873         {
874           if ( fCellSecondLabels[idCell] == -1 )  iNewLabel = 0;  // -1 is never a good second label.
875           if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) )  iNewLabel = 0;
876         }
877         if (iNewLabel == 1) 
878         {
879           Int_t * newLabelArray = (Int_t*)malloc((clus->GetNLabels()+1)*sizeof(Int_t)) ;
880           for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ ){
881             newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
882           }
883
884           newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
885           clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
886           delete [] newLabelArray;
887         }
888       }//positive cell number
889     }
890     
891     clusArray->Add(clus);
892   } // recPoints loop
893 }