[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / AliAnalysisTaskEMCALClusterize.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
/* $Id: $ */

//_________________________________________________________________________
// This analysis provides a new list of clusters to be used in other analysis
//
// Author: Gustavo Conesa Balbastre,
//         Adapted from analysis class from Deepa Thomas
//
//
//_________________________________________________________________________

// --- Root ---
#include "TString.h"
#include "TRefArray.h"
#include "TClonesArray.h"
#include "TTree.h"
#include "TGeoManager.h"

// --- AliRoot Analysis Steering
#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"
#include "AliESDEvent.h"
#include "AliGeomManager.h"
#include "AliVCaloCells.h"
#include "AliAODCaloCluster.h"
#include "AliCDBManager.h"
#include "AliCDBStorage.h"
#include "AliCDBEntry.h"
#include "AliLog.h"
#include "AliVEventHandler.h"

// --- EMCAL
#include "AliEMCALRecParam.h"
#include "AliEMCALAfterBurnerUF.h"
#include "AliEMCALGeometry.h"
#include "AliEMCALClusterizerNxN.h"
#include "AliEMCALClusterizerv1.h"
#include "AliEMCALRecPoint.h"
#include "AliEMCALDigit.h"
#include "AliCaloCalibPedestal.h"
#include "AliEMCALCalibData.h"

#include "AliAnalysisTaskEMCALClusterize.h"

ClassImp(AliAnalysisTaskEMCALClusterize)

//________________________________________________________________________
AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name) 
  : AliAnalysisTaskSE(name)
  , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"),  fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
  , fCalibData(0),       fPedestalData(0),     fOCDBpath("raw://")
  , fDigitsArr(0),       fClusterArr(0),       fCaloClusterArr(0)
  , fRecParam(0),        fClusterizer(0),      fUnfolder(0),           fJustUnfold(kFALSE) 
  , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"),     fFillAODFile(kTRUE)
  
  {
  //ctor
  for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0 ;
  fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
  fClusterArr      = new TObjArray(100);
  fCaloClusterArr  = new TObjArray(100);
  fRecParam        = new AliEMCALRecParam;
}

//________________________________________________________________________
AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize() 
  : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
  , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"),  fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
  , fCalibData(0),       fPedestalData(0),     fOCDBpath("raw://")
  , fDigitsArr(0),       fClusterArr(0),       fCaloClusterArr(0)
  , fRecParam(0),        fClusterizer(0),      fUnfolder(0),           fJustUnfold(kFALSE) 
  , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"),     fFillAODFile(kFALSE)
{
  // Constructor
  for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0 ;
  fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
  fClusterArr      = new TObjArray(100);
  fCaloClusterArr  = new TObjArray(100);
  fRecParam        = new AliEMCALRecParam;
}

//________________________________________________________________________
AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
{
  //dtor 
  
  if (fDigitsArr){
    fDigitsArr->Clear("C");
    delete fDigitsArr; fDigitsArr=0;
  }
  
  if (fClusterArr){
    fClusterArr->Delete();
    delete fClusterArr; fClusterArr=0;
  }
  
  if (fCaloClusterArr){
    fCaloClusterArr->Delete();
    delete fCaloClusterArr; fCaloClusterArr=0;
  }

  if(fClusterizer) {delete fClusterizer; fClusterizer = 0;}
  if(fGeom)        {delete fGeom;        fGeom        = 0;}
  if(fUnfolder)    {delete fUnfolder;    fUnfolder    = 0;}

}

//-------------------------------------------------------------------
void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
{
  // Init geometry, create list of output clusters
  
  fGeom =  AliEMCALGeometry::GetInstance(fGeomName) ;	

  fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
  fOutputAODBranch->SetName(fOutputAODBranchName);
  AddAODBranch("TClonesArray", &fOutputAODBranch);

}

//________________________________________________________________________
void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event
  
  //printf("___ Event __ %d __\n",(Int_t)Entry());
  
  //Remove the contents of output list set in the previous event 
  fOutputAODBranch->Clear("C");
  
  AliVEvent * event = InputEvent();
  if (!event) {
    printf("ERROR: event not available\n");
    return;
  }
  
  //Magic line to write events to AOD file
  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
  
  //-------------------------------------------------------------------------------------
  //Set the geometry matrix, for the first event, skip the rest
  //-------------------------------------------------------------------------------------
  if(!fGeomMatrixSet){
    if(fLoadGeomMatrices){
      for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
        if(fGeomMatrix[mod]){
          if(DebugLevel() > 1) 
            fGeomMatrix[mod]->Print();
          fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;  
        }
        fGeomMatrixSet=kTRUE;
      }//SM loop
    }//Load matrices
    else if(!gGeoManager){
      printf("AliAnalysisTaskEMCALClusterize::UserExec() - Get geo matrices from data\n");
      //Still not implemented in AOD, just a workaround to be able to work at least with ESDs	
      if(!strcmp(event->GetName(),"AliAODEvent")) {
        if(DebugLevel() > 1) 
          printf("AliAnalysisTaskCaloFilter Use ideal geometry, values geometry matrix not kept in AODs.\n");
      }//AOD
      else {	
        if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALClusterize Load Misaligned matrices. \n");
        AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
        if(!esd) {
          printf("AliAnalysisTaskCaloFilter::UserExec() - This event does not contain ESDs?");
          return;
        }
        for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
          if(DebugLevel() > 1) 
            esd->GetEMCALMatrix(mod)->Print();
          if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
        } 
        fGeomMatrixSet=kTRUE;
      }//ESD
    }//Load matrices from Data 
  }//first event

  //Get the list of cells needed for unfolding and reclustering
  AliVCaloCells *cells= event->GetEMCALCells();
  Int_t nClustersOrg = 0;
  //-------------------------------------------
  //---------Unfolding clusters----------------
  //-------------------------------------------
  if (fJustUnfold) {
    
    //Fill the array with the EMCAL clusters, copy them
    for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
    {
      AliVCluster *clus = event->GetCaloCluster(i);
      if(clus->IsEMCAL()){        
        //printf("Org Cluster %d, E %f \n",i, clus->E());
        if(dynamic_cast<AliESDCaloCluster*> (clus)){
          fCaloClusterArr->Add( new AliESDCaloCluster(*(dynamic_cast<AliESDCaloCluster*> (clus))) );   
        }//ESD
        else{
          fCaloClusterArr->Add( new AliAODCaloCluster(*(dynamic_cast<AliAODCaloCluster*> (clus))) );   
        }//AOD
        nClustersOrg++;
      }
    }
    
    //Do the unfolding
    fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
    
    //CLEAN-UP
    fUnfolder->Clear();
    
    //printf("nClustersOrg %d\n",nClustersOrg);
  }
  //-------------------------------------------
  //---------- Recluster cells ----------------
  //-------------------------------------------
  
  else{
    //-------------------------------------------------------------------------------------
    //Transform CaloCells into Digits
    //-------------------------------------------------------------------------------------
    Int_t idigit =  0;
    Int_t id     = -1;
    Float_t amp  = -1; 
    Float_t time = -1; 
    
    Double_t cellAmplitude = 0;
    Double_t cellTime      = 0;
    Short_t  cellNumber    = 0;
        
    TTree *digitsTree = new TTree("digitstree","digitstree");
    digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
    
    for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
    {
      if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
        break;
      
      time = cellTime;
      amp  = cellAmplitude;
      id   = cellNumber;
      
      AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
      digit->SetId(id);
      digit->SetAmplitude(amp);
      digit->SetTime(time);
      digit->SetTimeR(time);
      digit->SetIndexInList(idigit);
      digit->SetType(AliEMCALDigit::kHG);
      //if(Entry()==0) printf("Digit: Id %d, amp %f, time %e, index %d\n",id, amp,time,idigit);
      idigit++;
    }
    
    //Fill the tree with digits
    digitsTree->Fill();
    
    //-------------------------------------------------------------------------------------
    //Do the clusterization
    //-------------------------------------------------------------------------------------        
    TTree *clustersTree = new TTree("clustertree","clustertree");
    
    fClusterizer->SetInput(digitsTree);
    fClusterizer->SetOutput(clustersTree);
    fClusterizer->Digits2Clusters("");
    
    //-------------------------------------------------------------------------------------
    //Transform the recpoints into AliVClusters
    //-------------------------------------------------------------------------------------
    
    clustersTree->SetBranchStatus("*",0); //disable all branches
    clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
    
    TBranch *branch = clustersTree->GetBranch("EMCALECARP");
    branch->SetAddress(&fClusterArr);
    branch->GetEntry(0);
    
    RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
    
    //---CLEAN UP-----
    fClusterizer->Clear();
    fDigitsArr  ->Clear("C");
    fClusterArr ->Delete(); // Do not Clear(), it leaks, why?

    clustersTree->Delete("all");
    digitsTree  ->Delete("all");
    
  }
  
  //-------------------------------------------------------------------------------------
  //Put the new clusters in the AOD list
  //-------------------------------------------------------------------------------------
  
  Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntries();
  //printf("New clusters %d\n",kNumberOfCaloClusters);
  //if(nClustersOrg!=kNumberOfCaloClusters) printf("Different number: Org %d, New %d\n",nClustersOrg,kNumberOfCaloClusters);
  for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
    AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
    //if(Entry()==0)printf("newCluster E %f\n", newCluster->E());
    new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
  }
  
  //---CLEAN UP-----
  fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
  
 }      


//_____________________________________________________________________
Bool_t AliAnalysisTaskEMCALClusterize::UserNotify()
{
  //Access to OCDB stuff
  if(DebugLevel() > 1 )printf(" AliAnalysisTaskEMCALClusterize::UserNotify() - Begin \n");
  AliVEvent * event = InputEvent();
  if (event)
  {
    fGeom = AliEMCALGeometry::GetInstance(fGeomName);
    AliCDBManager *cdb = AliCDBManager::Instance();
    cdb->SetDefaultStorage(fOCDBpath.Data());
    cdb->SetRun(event->GetRunNumber());
    //
    // EMCAL from RAW OCDB
    if (fOCDBpath.Contains("alien:"))
    {
      cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
      cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
    }
    TString path = cdb->GetDefaultStorage()->GetBaseFolder();
    
    // init parameters:
    //Get calibration parameters	
    if(!fCalibData)
    {
      AliCDBEntry *entry = (AliCDBEntry*) 
	    AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
      if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
    }
    
    if(!fCalibData)
      AliFatal("Calibration parameters not found in CDB!");
    
    //Get calibration parameters	
    if(!fPedestalData)
    {
      AliCDBEntry *entry = (AliCDBEntry*) 
	    AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
      if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
    }
    
    if(!fPedestalData)
      AliFatal("Dead map not found in CDB!");
    
    //      cout << "[i] Change of run number: " << fAOD->GetRunNumber() << endl;
    InitClusterization();
    
  }
  else
  {
    printf(" AliAnalysisTaskEMCALClusterize::UserNotify() - Event not available!!! \n");
  }

  return kTRUE;
  
}

//________________________________________________________________________________________
void AliAnalysisTaskEMCALClusterize::InitClusterization()
{
  //Select clusterization/unfolding algorithm and set all the needed parameters
  
  if(!fJustUnfold){
    
    //First init the clusterizer
    if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
      fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
    else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
      fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
    else{
      printf("AliAnalysisTaskEMCALClusterize::InitClusterization() - Clusterizer < %d > not available",
             fRecParam->GetClusterizerFlag());
      abort();
    }
    
    //Now set the parameters
    fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
    fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
    fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );    
    fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
    fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
    fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
    fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
    fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
    fClusterizer->SetInputCalibrated       ( kTRUE                               );
    
    //In case of unfolding after clusterization is requested, set the corresponding parameters
    if(fRecParam->GetUnfold()){
      
      Int_t i=0;
      for (i = 0; i < 8; i++) {
        fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
      }//end of loop over parameters
      for (i = 0; i < 3; i++) {
        fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
        fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
      }//end of loop over parameters
      
      fClusterizer->InitClusterUnfolding();
            
    }// to unfold
    
    
  }else{
    //Now init the unfolding afterburner 
    fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
 }
}
//________________________________________________________________________________________
void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, TObjArray *recPoints, TObjArray *clusArray)
{
  // Restore clusters from recPoints
  // Cluster energy, global position, cells and their amplitude fractions are restored
  
  for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++)
  {
    AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i);
    
    const Int_t ncells = recPoint->GetMultiplicity();
    Int_t ncells_true = 0;
    
    // cells and their amplitude fractions
    UShort_t   absIds[ncells];  
    Double32_t ratios[ncells];
    
    for (Int_t c = 0; c < ncells; c++) {
      AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]);
      
      absIds[ncells_true] = digit->GetId();
      ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
      
      if (ratios[ncells_true] > 0.001) ncells_true++;
    }
    
    if (ncells_true < 1) {
      Warning("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters", "skipping cluster with no cells");
      continue;
    }
    
    TVector3 gpos;
    Float_t g[3];
    
    // calculate new cluster position
    recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr);
    recPoint->GetGlobalPosition(gpos);
    gpos.GetXYZ(g);
    
    // create a new cluster
    AliAODCaloCluster *clus = new AliAODCaloCluster();
    clus->SetType(AliVCluster::kEMCALClusterv1);
    clus->SetE(recPoint->GetEnergy());
    clus->SetPosition(g);
    clus->SetNCells(ncells_true);
    clus->SetCellsAbsId(absIds);
    clus->SetCellsAmplitudeFraction(ratios);
    clus->SetDispersion(recPoint->GetDispersion());
    clus->SetChi2(-1); //not yet implemented
    clus->SetTOF(recPoint->GetTime()) ; //time-of-fligh
    clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
    Float_t elipAxis[2];
    recPoint->GetElipsAxis(elipAxis);
    clus->SetM02(elipAxis[0]*elipAxis[0]) ;
    clus->SetM20(elipAxis[1]*elipAxis[1]) ;

    clusArray->Add(clus);

  } // recPoints loop
  
}
Commit	Line	Data
6fdd30c4	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