new task for reclusterization in EMCAL
[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
33 // --- AliRoot Analysis Steering
34 #include "AliAnalysisTask.h"
35 #include "AliAnalysisManager.h"
36 #include "AliESDEvent.h"
37 #include "AliGeomManager.h"
38 #include "AliVCaloCells.h"
39 #include "AliAODCaloCluster.h"
40 #include "AliCDBManager.h"
41 #include "AliCDBStorage.h"
42 #include "AliCDBEntry.h"
43 #include "AliLog.h"
44 #include "AliVEventHandler.h"
45
46 // --- EMCAL
47 #include "AliEMCALRecParam.h"
48 #include "AliEMCALAfterBurnerUF.h"
49 #include "AliEMCALGeometry.h"
50 #include "AliEMCALClusterizerNxN.h"
51 #include "AliEMCALClusterizerv1.h"
52 #include "AliEMCALRecPoint.h"
53 #include "AliEMCALDigit.h"
54 #include "AliCaloCalibPedestal.h"
55 #include "AliEMCALCalibData.h"
56
57 #include "AliAnalysisTaskEMCALClusterize.h"
58
59 ClassImp(AliAnalysisTaskEMCALClusterize)
60
61 //________________________________________________________________________
62 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name) 
63   : AliAnalysisTaskSE(name)
64   , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"),  fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
65   , fCalibData(0),       fPedestalData(0),     fOCDBpath("raw://")
66   , fDigitsArr(0),       fClusterArr(0),       fCaloClusterArr(0)
67   , fRecParam(0),        fClusterizer(0),      fUnfolder(0),           fJustUnfold(kFALSE) 
68   , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"),     fFillAODFile(kTRUE)
69   
70   {
71   //ctor
72   for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0 ;
73   fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
74   fClusterArr      = new TObjArray(100);
75   fCaloClusterArr  = new TObjArray(100);
76   fRecParam        = new AliEMCALRecParam;
77 }
78
79 //________________________________________________________________________
80 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize() 
81   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
82   , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"),  fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
83   , fCalibData(0),       fPedestalData(0),     fOCDBpath("raw://")
84   , fDigitsArr(0),       fClusterArr(0),       fCaloClusterArr(0)
85   , fRecParam(0),        fClusterizer(0),      fUnfolder(0),           fJustUnfold(kFALSE) 
86   , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"),     fFillAODFile(kFALSE)
87 {
88   // Constructor
89   for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0 ;
90   fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
91   fClusterArr      = new TObjArray(100);
92   fCaloClusterArr  = new TObjArray(100);
93   fRecParam        = new AliEMCALRecParam;
94 }
95
96 //________________________________________________________________________
97 AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
98 {
99   //dtor 
100   
101   if (fDigitsArr){
102     fDigitsArr->Clear("C");
103     delete fDigitsArr; fDigitsArr=0;
104   }
105   
106   if (fClusterArr){
107     fClusterArr->Delete();
108     delete fClusterArr; fClusterArr=0;
109   }
110   
111   if (fCaloClusterArr){
112     fCaloClusterArr->Delete();
113     delete fCaloClusterArr; fCaloClusterArr=0;
114   }
115
116   if(fClusterizer) {delete fClusterizer; fClusterizer = 0;}
117   if(fGeom)        {delete fGeom;        fGeom        = 0;}
118   if(fUnfolder)    {delete fUnfolder;    fUnfolder    = 0;}
119
120 }
121
122 //-------------------------------------------------------------------
123 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
124 {
125   // Init geometry, create list of output clusters
126   
127   fGeom =  AliEMCALGeometry::GetInstance(fGeomName) ;   
128
129   fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
130   fOutputAODBranch->SetName(fOutputAODBranchName);
131   AddAODBranch("TClonesArray", &fOutputAODBranch);
132
133 }
134
135 //________________________________________________________________________
136 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *) 
137 {
138   // Main loop
139   // Called for each event
140   
141   //printf("___ Event __ %d __\n",(Int_t)Entry());
142   
143   //Remove the contents of output list set in the previous event 
144   fOutputAODBranch->Clear("C");
145   
146   AliVEvent * event = InputEvent();
147   if (!event) {
148     printf("ERROR: event not available\n");
149     return;
150   }
151   
152   //Magic line to write events to AOD file
153   AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
154   
155   //-------------------------------------------------------------------------------------
156   //Set the geometry matrix, for the first event, skip the rest
157   //-------------------------------------------------------------------------------------
158   if(!fGeomMatrixSet){
159     if(fLoadGeomMatrices){
160       for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
161         if(fGeomMatrix[mod]){
162           if(DebugLevel() > 1) 
163             fGeomMatrix[mod]->Print();
164           fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;  
165         }
166         fGeomMatrixSet=kTRUE;
167       }//SM loop
168     }//Load matrices
169     else if(!gGeoManager){
170       printf("AliAnalysisTaskEMCALClusterize::UserExec() - Get geo matrices from data\n");
171       //Still not implemented in AOD, just a workaround to be able to work at least with ESDs   
172       if(!strcmp(event->GetName(),"AliAODEvent")) {
173         if(DebugLevel() > 1) 
174           printf("AliAnalysisTaskCaloFilter Use ideal geometry, values geometry matrix not kept in AODs.\n");
175       }//AOD
176       else {    
177         if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALClusterize Load Misaligned matrices. \n");
178         AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
179         if(!esd) {
180           printf("AliAnalysisTaskCaloFilter::UserExec() - This event does not contain ESDs?");
181           return;
182         }
183         for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
184           if(DebugLevel() > 1) 
185             esd->GetEMCALMatrix(mod)->Print();
186           if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
187         } 
188         fGeomMatrixSet=kTRUE;
189       }//ESD
190     }//Load matrices from Data 
191   }//first event
192
193   //Get the list of cells needed for unfolding and reclustering
194   AliVCaloCells *cells= event->GetEMCALCells();
195   Int_t nClustersOrg = 0;
196   //-------------------------------------------
197   //---------Unfolding clusters----------------
198   //-------------------------------------------
199   if (fJustUnfold) {
200     
201     //Fill the array with the EMCAL clusters, copy them
202     for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
203     {
204       AliVCluster *clus = event->GetCaloCluster(i);
205       if(clus->IsEMCAL()){        
206         //printf("Org Cluster %d, E %f \n",i, clus->E());
207         if(dynamic_cast<AliESDCaloCluster*> (clus)){
208           fCaloClusterArr->Add( new AliESDCaloCluster(*(dynamic_cast<AliESDCaloCluster*> (clus))) );   
209         }//ESD
210         else{
211           fCaloClusterArr->Add( new AliAODCaloCluster(*(dynamic_cast<AliAODCaloCluster*> (clus))) );   
212         }//AOD
213         nClustersOrg++;
214       }
215     }
216     
217     //Do the unfolding
218     fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
219     
220     //CLEAN-UP
221     fUnfolder->Clear();
222     
223     //printf("nClustersOrg %d\n",nClustersOrg);
224   }
225   //-------------------------------------------
226   //---------- Recluster cells ----------------
227   //-------------------------------------------
228   
229   else{
230     //-------------------------------------------------------------------------------------
231     //Transform CaloCells into Digits
232     //-------------------------------------------------------------------------------------
233     Int_t idigit =  0;
234     Int_t id     = -1;
235     Float_t amp  = -1; 
236     Float_t time = -1; 
237     
238     Double_t cellAmplitude = 0;
239     Double_t cellTime      = 0;
240     Short_t  cellNumber    = 0;
241         
242     TTree *digitsTree = new TTree("digitstree","digitstree");
243     digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
244     
245     for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
246     {
247       if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
248         break;
249       
250       time = cellTime;
251       amp  = cellAmplitude;
252       id   = cellNumber;
253       
254       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
255       digit->SetId(id);
256       digit->SetAmplitude(amp);
257       digit->SetTime(time);
258       digit->SetTimeR(time);
259       digit->SetIndexInList(idigit);
260       digit->SetType(AliEMCALDigit::kHG);
261       //if(Entry()==0) printf("Digit: Id %d, amp %f, time %e, index %d\n",id, amp,time,idigit);
262       idigit++;
263     }
264     
265     //Fill the tree with digits
266     digitsTree->Fill();
267     
268     //-------------------------------------------------------------------------------------
269     //Do the clusterization
270     //-------------------------------------------------------------------------------------        
271     TTree *clustersTree = new TTree("clustertree","clustertree");
272     
273     fClusterizer->SetInput(digitsTree);
274     fClusterizer->SetOutput(clustersTree);
275     fClusterizer->Digits2Clusters("");
276     
277     //-------------------------------------------------------------------------------------
278     //Transform the recpoints into AliVClusters
279     //-------------------------------------------------------------------------------------
280     
281     clustersTree->SetBranchStatus("*",0); //disable all branches
282     clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
283     
284     TBranch *branch = clustersTree->GetBranch("EMCALECARP");
285     branch->SetAddress(&fClusterArr);
286     branch->GetEntry(0);
287     
288     RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
289     
290     //---CLEAN UP-----
291     fClusterizer->Clear();
292     fDigitsArr  ->Clear("C");
293     fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
294
295     clustersTree->Delete("all");
296     digitsTree  ->Delete("all");
297     
298   }
299   
300   //-------------------------------------------------------------------------------------
301   //Put the new clusters in the AOD list
302   //-------------------------------------------------------------------------------------
303   
304   Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntries();
305   //printf("New clusters %d\n",kNumberOfCaloClusters);
306   //if(nClustersOrg!=kNumberOfCaloClusters) printf("Different number: Org %d, New %d\n",nClustersOrg,kNumberOfCaloClusters);
307   for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
308     AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
309     //if(Entry()==0)printf("newCluster E %f\n", newCluster->E());
310     new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
311   }
312   
313   //---CLEAN UP-----
314   fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
315   
316  }      
317
318
319 //_____________________________________________________________________
320 Bool_t AliAnalysisTaskEMCALClusterize::UserNotify()
321 {
322   //Access to OCDB stuff
323   if(DebugLevel() > 1 )printf(" AliAnalysisTaskEMCALClusterize::UserNotify() - Begin \n");
324   AliVEvent * event = InputEvent();
325   if (event)
326   {
327     fGeom = AliEMCALGeometry::GetInstance(fGeomName);
328     AliCDBManager *cdb = AliCDBManager::Instance();
329     cdb->SetDefaultStorage(fOCDBpath.Data());
330     cdb->SetRun(event->GetRunNumber());
331     //
332     // EMCAL from RAW OCDB
333     if (fOCDBpath.Contains("alien:"))
334     {
335       cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
336       cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
337     }
338     TString path = cdb->GetDefaultStorage()->GetBaseFolder();
339     
340     // init parameters:
341     //Get calibration parameters        
342     if(!fCalibData)
343     {
344       AliCDBEntry *entry = (AliCDBEntry*) 
345             AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
346       if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
347     }
348     
349     if(!fCalibData)
350       AliFatal("Calibration parameters not found in CDB!");
351     
352     //Get calibration parameters        
353     if(!fPedestalData)
354     {
355       AliCDBEntry *entry = (AliCDBEntry*) 
356             AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
357       if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
358     }
359     
360     if(!fPedestalData)
361       AliFatal("Dead map not found in CDB!");
362     
363     //      cout << "[i] Change of run number: " << fAOD->GetRunNumber() << endl;
364     InitClusterization();
365     
366   }
367   else
368   {
369     printf(" AliAnalysisTaskEMCALClusterize::UserNotify() - Event not available!!! \n");
370   }
371
372   return kTRUE;
373   
374 }
375
376 //________________________________________________________________________________________
377 void AliAnalysisTaskEMCALClusterize::InitClusterization()
378 {
379   //Select clusterization/unfolding algorithm and set all the needed parameters
380   
381   if(!fJustUnfold){
382     
383     //First init the clusterizer
384     if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
385       fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
386     else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
387       fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
388     else{
389       printf("AliAnalysisTaskEMCALClusterize::InitClusterization() - Clusterizer < %d > not available",
390              fRecParam->GetClusterizerFlag());
391       abort();
392     }
393     
394     //Now set the parameters
395     fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
396     fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
397     fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );    
398     fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
399     fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
400     fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
401     fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
402     fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
403     fClusterizer->SetInputCalibrated       ( kTRUE                               );
404     
405     //In case of unfolding after clusterization is requested, set the corresponding parameters
406     if(fRecParam->GetUnfold()){
407       
408       Int_t i=0;
409       for (i = 0; i < 8; i++) {
410         fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
411       }//end of loop over parameters
412       for (i = 0; i < 3; i++) {
413         fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
414         fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
415       }//end of loop over parameters
416       
417       fClusterizer->InitClusterUnfolding();
418             
419     }// to unfold
420     
421     
422   }else{
423     //Now init the unfolding afterburner 
424     fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
425  }
426 }
427 //________________________________________________________________________________________
428 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, TObjArray *recPoints, TObjArray *clusArray)
429 {
430   // Restore clusters from recPoints
431   // Cluster energy, global position, cells and their amplitude fractions are restored
432   
433   for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++)
434   {
435     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i);
436     
437     const Int_t ncells = recPoint->GetMultiplicity();
438     Int_t ncells_true = 0;
439     
440     // cells and their amplitude fractions
441     UShort_t   absIds[ncells];  
442     Double32_t ratios[ncells];
443     
444     for (Int_t c = 0; c < ncells; c++) {
445       AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]);
446       
447       absIds[ncells_true] = digit->GetId();
448       ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
449       
450       if (ratios[ncells_true] > 0.001) ncells_true++;
451     }
452     
453     if (ncells_true < 1) {
454       Warning("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters", "skipping cluster with no cells");
455       continue;
456     }
457     
458     TVector3 gpos;
459     Float_t g[3];
460     
461     // calculate new cluster position
462     recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr);
463     recPoint->GetGlobalPosition(gpos);
464     gpos.GetXYZ(g);
465     
466     // create a new cluster
467     AliAODCaloCluster *clus = new AliAODCaloCluster();
468     clus->SetType(AliVCluster::kEMCALClusterv1);
469     clus->SetE(recPoint->GetEnergy());
470     clus->SetPosition(g);
471     clus->SetNCells(ncells_true);
472     clus->SetCellsAbsId(absIds);
473     clus->SetCellsAmplitudeFraction(ratios);
474     clus->SetDispersion(recPoint->GetDispersion());
475     clus->SetChi2(-1); //not yet implemented
476     clus->SetTOF(recPoint->GetTime()) ; //time-of-fligh
477     clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
478     Float_t elipAxis[2];
479     recPoint->GetElipsAxis(elipAxis);
480     clus->SetM02(elipAxis[0]*elipAxis[0]) ;
481     clus->SetM20(elipAxis[1]*elipAxis[1]) ;
482
483     clusArray->Add(clus);
484
485   } // recPoints loop
486   
487 }
488
489
490