[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 "AliEMCALRecoUtils.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)
  , fRun(-1),            fRecoUtils(0)
  
  {
  //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;
  fBranchNames     = "ESD:AliESDHeader.,EMCALCells.";
  fRecoUtils       = new AliEMCALRecoUtils();
}

//________________________________________________________________________
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)
  , fRun(-1),            fRecoUtils(0)
{
  // 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;
  fBranchNames     = "ESD:AliESDHeader.,EMCALCells.";
  fRecoUtils       = new AliEMCALRecoUtils();
}

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

  if(fClusterizer) delete fClusterizer;
  if(fUnfolder)    delete fUnfolder;   
  if(fRecoUtils)   delete fRecoUtils;

}

//-------------------------------------------------------------------
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);
  Info("UserCreateOutputObjects","Create Branch: %s",fOutputAODBranchName.Data());
}

//________________________________________________________________________
void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event
  
  //Remove the contents of output list set in the previous event 
  fOutputAODBranch->Clear("C");
  
  AliVEvent   * event    = InputEvent();
  AliESDEvent * esdevent = dynamic_cast<AliESDEvent*>(InputEvent());

  if (!event) {
    Error("UserExec","Event not available");
    return;
  }
  
  //Magic line to write events to AOD file
  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
  LoadBranches();
  
  AccessOCDB();

  //-------------------------------------------------------------------------------------
  //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){
      Info("UserExec","Get geo matrices from data");
      //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) 
          Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
      }//AOD
      else {	
        if(DebugLevel() > 1) 
          Info("UserExec","AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
        AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
        if(!esd) {
          Error("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());
        AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
        AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
        if     (esdCluster){
          fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );   
        }//ESD
        else if(aodCluster){
          fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );   
        }//AOD
        else 
          Warning("UserExec()"," - Wrong CaloCluster type?");
        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");
  }
  
  //Recalculate track-matching for the new clusters, only with ESDs
  if(esdevent)fRecoUtils->FindMatches(esdevent,fCaloClusterArr);

  
  //-------------------------------------------------------------------------------------
  //Put the new clusters in the AOD list
  //-------------------------------------------------------------------------------------
  
  Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntries();
  //Info("UserExec","New clusters %d",kNumberOfCaloClusters);
  //if(nClustersOrg!=kNumberOfCaloClusters) Info("UserExec","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) Info("UserExec","newCluster E %f\n", newCluster->E());
    
    //Add matched track, if any, only with ESDs
    if(esdevent){
      Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
      if(trackIndex >= 0){
        newCluster->AddTrackMatched(event->GetTrack(trackIndex));
        if(DebugLevel() > 1) 
          Info("UserExec","Matched Track index %d to new cluster %d \n",trackIndex,i);
      }
    }
    new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
  }
  
  //---CLEAN UP-----
  fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
}      

//_____________________________________________________________________
Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
{
  //Access to OCDB stuff

  AliVEvent * event = InputEvent();
  if (!event)
  {
    Warning("AccessOCDB","Event not available!!!");
    return kFALSE;
  }

  if (event->GetRunNumber()==fRun)
    return kTRUE;
  fRun = event->GetRunNumber();

  if(DebugLevel() > 1 )
    Info("AccessODCD()"," Begin");

  fGeom = AliEMCALGeometry::GetInstance(fGeomName);
  
  
  AliCDBManager *cdb = AliCDBManager::Instance();
    

  if (fOCDBpath.Length()){
    cdb->SetDefaultStorage(fOCDBpath.Data());
    Info("AccessOCDB","Default storage %s",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!");

  InitClusterization();
  
  return kTRUE;
}

//________________________________________________________________________________________
void AliAnalysisTaskEMCALClusterize::InitClusterization()
{
  //Select clusterization/unfolding algorithm and set all the needed parameters
  
  if (fJustUnfold){
    // init the unfolding afterburner 
    delete fUnfolder;
    fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
    return;
 }

  //First init the clusterizer
  delete fClusterizer;
  if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
    fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) 
    fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
  else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
   AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
   clusterizer->SetNRowDiff(2);
   clusterizer->SetNColDiff(2);
   fClusterizer = clusterizer;
  } else{
    AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
  }
    
  //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
}

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