[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / AliAnalysisTaskEMCALClusterizeFast.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$ */

// --- Root ---
#include <TClonesArray.h>
#include <TGeoManager.h>
#include <TObjArray.h>
#include <TString.h>
#include <TTree.h>
#include "AliAODCaloCluster.h"
#include "AliAODEvent.h"
#include "AliAnalysisManager.h"
#include "AliCDBEntry.h"
#include "AliCDBManager.h"
#include "AliCaloCalibPedestal.h"
#include "AliEMCALAfterBurnerUF.h"
#include "AliEMCALCalibData.h"
#include "AliEMCALClusterizerNxN.h"
#include "AliEMCALClusterizerv1.h"
#include "AliEMCALDigit.h"
#include "AliEMCALGeometry.h"
#include "AliEMCALRecParam.h"
#include "AliEMCALRecPoint.h"
#include "AliEMCALRecoUtils.h"
#include "AliESDEvent.h"
#include "AliLog.h"

#include "AliAnalysisTaskEMCALClusterizeFast.h"

ClassImp(AliAnalysisTaskEMCALClusterizeFast)

//________________________________________________________________________
AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast() 
  : AliAnalysisTaskSE(), 
    fRun(0),
    fDigitsArr(0),       
    fClusterArr(0),       
    fRecParam(0),
    fClusterizer(0),
    fUnfolder(0),
    fJustUnfold(kFALSE),
    fGeomName(),
    fGeomMatrixSet(kFALSE), 
    fLoadGeomMatrices(kFALSE),
    fOCDBpath(),
    fCalibData(0),
    fPedestalData(0),
    fOutputAODBranch(0),
    fOutputAODBrName(),
    fRecoUtils(0),
    fLoadCalib(0),
    fLoadPed(0)
{ 
  // Constructor
}

//________________________________________________________________________
AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast(const char *name) 
  : AliAnalysisTaskSE(name), 
    fRun(0),
    fDigitsArr(0),       
    fClusterArr(0),       
    fRecParam(new AliEMCALRecParam),
    fClusterizer(0),
    fUnfolder(0),
    fJustUnfold(kFALSE),
    fGeomName("EMCAL_FIRSTYEARV1"),
    fGeomMatrixSet(kFALSE), 
    fLoadGeomMatrices(kFALSE),
    fOCDBpath(),
    fCalibData(0),
    fPedestalData(0),
    fOutputAODBranch(0),
    fOutputAODBrName(),
    fRecoUtils(0),
    fLoadCalib(0),
    fLoadPed(0)
{ 
  // Constructor

  fBranchNames     = "ESD:AliESDHeader.,EMCALCells. AOD:header,emcalCells";
  for(Int_t i = 0; i < 10; ++i) 
    fGeomMatrix[i] = 0;
}

//________________________________________________________________________
AliAnalysisTaskEMCALClusterizeFast::~AliAnalysisTaskEMCALClusterizeFast()
{
  // Destructor.

  delete fDigitsArr; 
  delete fClusterizer;
  delete fUnfolder;   
  delete fRecoUtils;
}

//-------------------------------------------------------------------
void AliAnalysisTaskEMCALClusterizeFast::UserCreateOutputObjects()
{
  // Create output objects.

  if (!fOutputAODBrName.IsNull()) {
    fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
    fOutputAODBranch->SetName(fOutputAODBrName);
    AddAODBranch("TClonesArray", &fOutputAODBranch);
    AliInfo(Form("Created Branch: %s",fOutputAODBrName.Data()));
  }
}

//________________________________________________________________________
void AliAnalysisTaskEMCALClusterizeFast::UserExec(Option_t *) 
{
  // Main loop, called for each event

  // remove the contents of output list set in the previous event 
  if (fOutputAODBranch)
    fOutputAODBranch->Clear("C");

  AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
  AliAODEvent *aodevent = dynamic_cast<AliAODEvent*>(InputEvent());

  if (!esdevent&&!aodevent) {
    Error("UserExec","Event not available");
    return;
  }

  LoadBranches();
  
  Init();

  if (fJustUnfold) {
    AliWarning("Unfolding not implemented");
  } else {
    if (esdevent) 
      FillDigitsArray(esdevent);
    else 
      FillDigitsArray(aodevent);
    fClusterizer->Digits2Clusters("");
    if (esdevent && fRecoUtils)
      fRecoUtils->FindMatches(esdevent,fClusterArr);
    if (fOutputAODBranch) {
      RecPoints2Clusters();
    }
    if (esdevent) {
      UpdateCells(esdevent);
    } else {
      UpdateCells(aodevent);
    }
  }
}

//________________________________________________________________________
void AliAnalysisTaskEMCALClusterizeFast::UpdateCells(AliAODEvent *event)
{
  // Update cells in case re-calibration was done.

  if (!fCalibData)
    return;

  AliAODCaloCells *cells = event->GetEMCALCells();
  Int_t ncells = cells->GetNumberOfCells();
  for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
    Double_t cellAmplitude=0, cellTime=0;
    Short_t cellNumber=0;
    if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
      break;
    AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
    cellAmplitude = digit->GetCalibAmp();
    cells->SetCell(icell, cellNumber, cellAmplitude, cellTime);
    idigit++;
  }
  cells->Sort();
}

//________________________________________________________________________
void AliAnalysisTaskEMCALClusterizeFast::UpdateCells(AliESDEvent *event)
{
  // Update cells in case re-calibration was done.

  if (!fCalibData)
    return;

  AliESDCaloCells *cells = event->GetEMCALCells();
  Int_t ncells = cells->GetNumberOfCells();
  for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
    Double_t cellAmplitude=0, cellTime=0;
    Short_t cellNumber=0;
    if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
      break;
    AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
    cellAmplitude = digit->GetCalibAmp();
    cells->SetCell(icell, cellNumber, cellAmplitude, cellTime);
    idigit++;
  }
  cells->Sort();
}

//________________________________________________________________________
void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray(AliAODEvent *event)
{
  // Fill digits from cells.

  fDigitsArr->Clear("C");
  AliAODCaloCells *cells = event->GetEMCALCells();
  Int_t ncells = cells->GetNumberOfCells();
  if (ncells>fDigitsArr->GetSize())
    fDigitsArr->Expand(2*ncells);
  for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
    Double_t cellAmplitude=0, cellTime=0;
    Short_t cellNumber=0;
    if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
      break;
    AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
    digit->SetId(cellNumber);
    digit->SetAmplitude(cellAmplitude);
    digit->SetTime(cellTime);
    digit->SetTimeR(cellTime);
    digit->SetIndexInList(idigit);
    digit->SetType(AliEMCALDigit::kHG);
    idigit++;
  }
}

//________________________________________________________________________
void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray(AliESDEvent *event)
{
  // Fill digits from cells.

  fDigitsArr->Clear("C");
  AliESDCaloCells *cells = event->GetEMCALCells();
  Int_t ncells = cells->GetNumberOfCells();
  if (ncells>fDigitsArr->GetSize())
    fDigitsArr->Expand(2*ncells);
  for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
    Double_t cellAmplitude=0, cellTime=0;
    Short_t cellNumber=0;
    if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
      break;
    AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
    digit->SetId(cellNumber);
    digit->SetAmplitude(cellAmplitude);
    digit->SetTime(cellTime);
    digit->SetTimeR(cellTime);
    digit->SetIndexInList(idigit);
    digit->SetType(AliEMCALDigit::kHG);
    idigit++;
  }
}

//________________________________________________________________________________________
void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters()
{
  // Cluster energy, global position, cells and their amplitude fractions are restored.

  AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());

  Int_t Ncls = fClusterArr->GetEntriesFast();
  for(Int_t i=0, nout=0; i < Ncls; ++i) {
    AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
    Int_t ncells_true = 0;
    const Int_t ncells = recpoint->GetMultiplicity();
    UShort_t   absIds[ncells];  
    Double32_t ratios[ncells];
    Int_t *dlist = recpoint->GetDigitsList();
    Float_t *elist = recpoint->GetEnergiesList();
    for (Int_t c = 0; c < ncells; ++c) {
      AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
      absIds[ncells_true] = digit->GetId();
      ratios[ncells_true] = elist[c]/digit->GetAmplitude();
      if (ratios[ncells_true] > 0.001) 
        ++ncells_true;
    }
    
    if (ncells_true < 1) {
      AliWarning("Skipping cluster with no cells");
      continue;
    }
    
    // calculate new cluster position
    TVector3 gpos;
    Float_t g[3];

    recpoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
    recpoint->GetGlobalPosition(gpos);
    gpos.GetXYZ(g);
    
    AliAODCaloCluster *clus = static_cast<AliAODCaloCluster*>(fOutputAODBranch->New(nout++));
    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()); 
    if (esdevent && fRecoUtils) {
      Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
      if(trackIndex >= 0) {
        clus->AddTrackMatched(esdevent->GetTrack(trackIndex));
        if(DebugLevel() > 1) 
          AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
      }
    }
  }
}

//________________________________________________________________________________________
void AliAnalysisTaskEMCALClusterizeFast::Init()
{
  //Select clusterization/unfolding algorithm and set all the needed parameters
  
  AliVEvent * event = InputEvent();
  if (!event) {
    AliWarning("Event not available!!!");
    return;
  }

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

  if (fJustUnfold){
    // init the unfolding afterburner 
    delete fUnfolder;
    fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
    return;
  }

  AliCDBManager *cdb = AliCDBManager::Instance();
  if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
    cdb->SetDefaultStorage(fOCDBpath);
  if (fRun!=cdb->GetRun())
    cdb->SetRun(fRun);

  AliEMCALGeometry *geometry = AliEMCALGeometry::GetInstance(fGeomName);
  if (!geometry) {
    AliFatal("Geometry not available!!!");
    return;
  }

  if (!fGeomMatrixSet) {
    if (fLoadGeomMatrices) {
      for(Int_t mod=0; mod < (geometry->GetEMCGeometry())->GetNumberOfSuperModules(); ++mod) {
        if(fGeomMatrix[mod]){
          if(DebugLevel() > 2) 
            fGeomMatrix[mod]->Print();
          geometry->SetMisalMatrix(fGeomMatrix[mod],mod);  
        }
      }
    } else {
      for(Int_t mod=0; mod < geometry->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
        if(event->GetEMCALMatrix(mod)) {
          if(DebugLevel() > 2) 
            event->GetEMCALMatrix(mod)->Print();
          geometry->SetMisalMatrix(event->GetEMCALMatrix(mod),mod);
        }
      }
    }
    fGeomMatrixSet=kTRUE;
  }
  
  // setup digit array if needed
  if (!fDigitsArr) {
    fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
    fDigitsArr->SetOwner(1);
  }

  // then setup clusterizer
  delete fClusterizer;
  if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
    fClusterizer = new AliEMCALClusterizerv1(geometry);
  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) 
    fClusterizer = new AliEMCALClusterizerNxN(geometry);
  else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
   AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(geometry);
   clusterizer->SetNRowDiff(2);
   clusterizer->SetNColDiff(2);
   fClusterizer = clusterizer;
  } else{
    AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
  }
  fClusterizer->InitParameters(fRecParam);
  if (!fCalibData&&fLoadCalib) {
    AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
    if (entry) 
      fCalibData =  static_cast<AliEMCALCalibData*>(entry->GetObject());
    if (!fCalibData)
      AliFatal("Calibration parameters not found in CDB!");
  }
  if (!fPedestalData&&fLoadPed) {
    AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
    if (entry) 
      fPedestalData =  static_cast<AliCaloCalibPedestal*>(entry->GetObject());
  }
  if (fCalibData) {
    fClusterizer->SetInputCalibrated(kFALSE);   
    fClusterizer->SetCalibrationParameters(fCalibData);
    fClusterizer->SetCaloCalibPedestal(fPedestalData);
  } else {
    fClusterizer->SetInputCalibrated(kTRUE);   
  }
  fClusterizer->SetJustClusters(kTRUE);
  fClusterizer->SetDigitsArr(fDigitsArr);
  fClusterizer->SetOutput(0);
  fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
}
Commit	Line	Data
	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
	18	// --- Root ---
	19	#include <TClonesArray.h>
	20	#include <TGeoManager.h>
	21	#include <TObjArray.h>
	22	#include <TString.h>
	23	#include <TTree.h>
	24	#include "AliAODCaloCluster.h"
	25	#include "AliAODEvent.h"
	26	#include "AliAnalysisManager.h"
	27	#include "AliCDBEntry.h"
	28	#include "AliCDBManager.h"
	29	#include "AliCaloCalibPedestal.h"
	30	#include "AliEMCALAfterBurnerUF.h"
	31	#include "AliEMCALCalibData.h"
	32	#include "AliEMCALClusterizerNxN.h"
	33	#include "AliEMCALClusterizerv1.h"
	34	#include "AliEMCALDigit.h"
	35	#include "AliEMCALGeometry.h"
	36	#include "AliEMCALRecParam.h"
	37	#include "AliEMCALRecPoint.h"
	38	#include "AliEMCALRecoUtils.h"
	39	#include "AliESDEvent.h"
	40	#include "AliLog.h"
	41
	42	#include "AliAnalysisTaskEMCALClusterizeFast.h"
	43
	44	ClassImp(AliAnalysisTaskEMCALClusterizeFast)
	45
	46	//________________________________________________________________________
	47	AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast()
	48	: AliAnalysisTaskSE(),
	49	fRun(0),
	50	fDigitsArr(0),
	51	fClusterArr(0),
	52	fRecParam(0),
	53	fClusterizer(0),
	54	fUnfolder(0),
	55	fJustUnfold(kFALSE),
	56	fGeomName(),
	57	fGeomMatrixSet(kFALSE),
	58	fLoadGeomMatrices(kFALSE),
	59	fOCDBpath(),
	60	fCalibData(0),
	61	fPedestalData(0),
	62	fOutputAODBranch(0),
	63	fOutputAODBrName(),
	64	fRecoUtils(0),
	65	fLoadCalib(0),
	66	fLoadPed(0)
	67	{
	68	// Constructor
	69	}
	70
	71	//________________________________________________________________________
	72	AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast(const char *name)
	73	: AliAnalysisTaskSE(name),
	74	fRun(0),
	75	fDigitsArr(0),
	76	fClusterArr(0),
	77	fRecParam(new AliEMCALRecParam),
	78	fClusterizer(0),
	79	fUnfolder(0),
	80	fJustUnfold(kFALSE),
	81	fGeomName("EMCAL_FIRSTYEARV1"),
	82	fGeomMatrixSet(kFALSE),
	83	fLoadGeomMatrices(kFALSE),
	84	fOCDBpath(),
	85	fCalibData(0),
	86	fPedestalData(0),
	87	fOutputAODBranch(0),
	88	fOutputAODBrName(),
	89	fRecoUtils(0),
	90	fLoadCalib(0),
	91	fLoadPed(0)
	92	{
	93	// Constructor
	94
	95	fBranchNames = "ESD:AliESDHeader.,EMCALCells. AOD:header,emcalCells";
	96	for(Int_t i = 0; i < 10; ++i)
	97	fGeomMatrix[i] = 0;
	98	}
	99
	100	//________________________________________________________________________
	101	AliAnalysisTaskEMCALClusterizeFast::~AliAnalysisTaskEMCALClusterizeFast()
	102	{
	103	// Destructor.
	104
	105	delete fDigitsArr;
	106	delete fClusterizer;
	107	delete fUnfolder;
	108	delete fRecoUtils;
	109	}
	110
	111	//-------------------------------------------------------------------
	112	void AliAnalysisTaskEMCALClusterizeFast::UserCreateOutputObjects()
	113	{
	114	// Create output objects.
	115
	116	if (!fOutputAODBrName.IsNull()) {
	117	fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
	118	fOutputAODBranch->SetName(fOutputAODBrName);
	119	AddAODBranch("TClonesArray", &fOutputAODBranch);
	120	AliInfo(Form("Created Branch: %s",fOutputAODBrName.Data()));
	121	}
	122	}
	123
	124	//________________________________________________________________________
	125	void AliAnalysisTaskEMCALClusterizeFast::UserExec(Option_t *)
	126	{
	127	// Main loop, called for each event
	128
	129	// remove the contents of output list set in the previous event
	130	if (fOutputAODBranch)
	131	fOutputAODBranch->Clear("C");
	132
	133	AliESDEvent esdevent = dynamic_cast<AliESDEvent>(InputEvent());
	134	AliAODEvent aodevent = dynamic_cast<AliAODEvent>(InputEvent());
	135
	136	if (!esdevent&&!aodevent) {
	137	Error("UserExec","Event not available");
	138	return;
	139	}
	140
	141	LoadBranches();
	142
	143	Init();
	144
	145	if (fJustUnfold) {
	146	AliWarning("Unfolding not implemented");
	147	} else {
	148	if (esdevent)
	149	FillDigitsArray(esdevent);
	150	else
	151	FillDigitsArray(aodevent);
	152	fClusterizer->Digits2Clusters("");
	153	if (esdevent && fRecoUtils)
	154	fRecoUtils->FindMatches(esdevent,fClusterArr);
	155	if (fOutputAODBranch) {
	156	RecPoints2Clusters();
	157	}
	158	if (esdevent) {
	159	UpdateCells(esdevent);
	160	} else {
	161	UpdateCells(aodevent);
	162	}
	163	}
	164	}
	165
	166	//________________________________________________________________________
	167	void AliAnalysisTaskEMCALClusterizeFast::UpdateCells(AliAODEvent *event)
	168	{
	169	// Update cells in case re-calibration was done.
	170
	171	if (!fCalibData)
	172	return;
	173
	174	AliAODCaloCells *cells = event->GetEMCALCells();
	175	Int_t ncells = cells->GetNumberOfCells();
	176	for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
	177	Double_t cellAmplitude=0, cellTime=0;
	178	Short_t cellNumber=0;
	179	if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
	180	break;
	181	AliEMCALDigit digit = static_cast<AliEMCALDigit>(fDigitsArr->At(idigit));
	182	cellAmplitude = digit->GetCalibAmp();
	183	cells->SetCell(icell, cellNumber, cellAmplitude, cellTime);
	184	idigit++;
	185	}
	186	cells->Sort();
	187	}
	188
	189	//________________________________________________________________________
	190	void AliAnalysisTaskEMCALClusterizeFast::UpdateCells(AliESDEvent *event)
	191	{
	192	// Update cells in case re-calibration was done.
	193
	194	if (!fCalibData)
	195	return;
	196
	197	AliESDCaloCells *cells = event->GetEMCALCells();
	198	Int_t ncells = cells->GetNumberOfCells();
	199	for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
	200	Double_t cellAmplitude=0, cellTime=0;
	201	Short_t cellNumber=0;
	202	if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
	203	break;
	204	AliEMCALDigit digit = static_cast<AliEMCALDigit>(fDigitsArr->At(idigit));
	205	cellAmplitude = digit->GetCalibAmp();
	206	cells->SetCell(icell, cellNumber, cellAmplitude, cellTime);
	207	idigit++;
	208	}
	209	cells->Sort();
	210	}
	211
	212	//________________________________________________________________________
	213	void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray(AliAODEvent *event)
	214	{
	215	// Fill digits from cells.
	216
	217	fDigitsArr->Clear("C");
	218	AliAODCaloCells *cells = event->GetEMCALCells();
	219	Int_t ncells = cells->GetNumberOfCells();
	220	if (ncells>fDigitsArr->GetSize())
	221	fDigitsArr->Expand(2*ncells);
	222	for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
	223	Double_t cellAmplitude=0, cellTime=0;
	224	Short_t cellNumber=0;
	225	if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
	226	break;
	227	AliEMCALDigit digit = (AliEMCALDigit) fDigitsArr->New(idigit);
	228	digit->SetId(cellNumber);
	229	digit->SetAmplitude(cellAmplitude);
	230	digit->SetTime(cellTime);
	231	digit->SetTimeR(cellTime);
	232	digit->SetIndexInList(idigit);
	233	digit->SetType(AliEMCALDigit::kHG);
	234	idigit++;
	235	}
	236	}
	237
	238	//________________________________________________________________________
	239	void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray(AliESDEvent *event)
	240	{
	241	// Fill digits from cells.
	242
	243	fDigitsArr->Clear("C");
	244	AliESDCaloCells *cells = event->GetEMCALCells();
	245	Int_t ncells = cells->GetNumberOfCells();
	246	if (ncells>fDigitsArr->GetSize())
	247	fDigitsArr->Expand(2*ncells);
	248	for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
	249	Double_t cellAmplitude=0, cellTime=0;
	250	Short_t cellNumber=0;
	251	if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
	252	break;
	253	AliEMCALDigit digit = (AliEMCALDigit) fDigitsArr->New(idigit);
	254	digit->SetId(cellNumber);
	255	digit->SetAmplitude(cellAmplitude);
	256	digit->SetTime(cellTime);
	257	digit->SetTimeR(cellTime);
	258	digit->SetIndexInList(idigit);
	259	digit->SetType(AliEMCALDigit::kHG);
	260	idigit++;
	261	}
	262	}
	263
	264	//________________________________________________________________________________________
	265	void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters()
	266	{
	267	// Cluster energy, global position, cells and their amplitude fractions are restored.
	268
	269	AliESDEvent esdevent = dynamic_cast<AliESDEvent>(InputEvent());
	270
	271	Int_t Ncls = fClusterArr->GetEntriesFast();
	272	for(Int_t i=0, nout=0; i < Ncls; ++i) {
	273	AliEMCALRecPoint recpoint = static_cast<AliEMCALRecPoint>(fClusterArr->At(i));
	274	Int_t ncells_true = 0;
	275	const Int_t ncells = recpoint->GetMultiplicity();
	276	UShort_t absIds[ncells];
	277	Double32_t ratios[ncells];
	278	Int_t *dlist = recpoint->GetDigitsList();
	279	Float_t *elist = recpoint->GetEnergiesList();
	280	for (Int_t c = 0; c < ncells; ++c) {
	281	AliEMCALDigit digit = static_cast<AliEMCALDigit>(fDigitsArr->At(dlist[c]));
	282	absIds[ncells_true] = digit->GetId();
	283	ratios[ncells_true] = elist[c]/digit->GetAmplitude();
	284	if (ratios[ncells_true] > 0.001)
	285	++ncells_true;
	286	}
	287
	288	if (ncells_true < 1) {
	289	AliWarning("Skipping cluster with no cells");
	290	continue;
	291	}
	292
	293	// calculate new cluster position
	294	TVector3 gpos;
	295	Float_t g[3];
	296
	297	recpoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
	298	recpoint->GetGlobalPosition(gpos);
	299	gpos.GetXYZ(g);
	300
	301	AliAODCaloCluster clus = static_cast<AliAODCaloCluster>(fOutputAODBranch->New(nout++));
	302	clus->SetType(AliVCluster::kEMCALClusterv1);
	303	clus->SetE(recpoint->GetEnergy());
	304	clus->SetPosition(g);
	305	clus->SetNCells(ncells_true);
	306	clus->SetCellsAbsId(absIds);
	307	clus->SetCellsAmplitudeFraction(ratios);
	308	clus->SetDispersion(recpoint->GetDispersion());
	309	clus->SetChi2(-1); //not yet implemented
	310	clus->SetTOF(recpoint->GetTime()) ; //time-of-flight
	311	clus->SetNExMax(recpoint->GetNExMax()); //number of local maxima
	312	Float_t elipAxis[2];
	313	recpoint->GetElipsAxis(elipAxis);
	314	clus->SetM02(elipAxis[0]*elipAxis[0]) ;
	315	clus->SetM20(elipAxis[1]*elipAxis[1]) ;
	316	clus->SetDistToBadChannel(recpoint->GetDistanceToBadTower());
	317	if (esdevent && fRecoUtils) {
	318	Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
	319	if(trackIndex >= 0) {
	320	clus->AddTrackMatched(esdevent->GetTrack(trackIndex));
	321	if(DebugLevel() > 1)
	322	AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
	323	}
	324	}
	325	}
	326	}
	327
	328	//________________________________________________________________________________________
	329	void AliAnalysisTaskEMCALClusterizeFast::Init()
	330	{
	331	//Select clusterization/unfolding algorithm and set all the needed parameters
	332
	333	AliVEvent * event = InputEvent();
	334	if (!event) {
	335	AliWarning("Event not available!!!");
	336	return;
	337	}
	338
	339	if (event->GetRunNumber()==fRun)
	340	return;
	341	fRun = event->GetRunNumber();
	342
	343	if (fJustUnfold){
	344	// init the unfolding afterburner
	345	delete fUnfolder;
	346	fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
	347	return;
	348	}
	349
	350	AliCDBManager *cdb = AliCDBManager::Instance();
	351	if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
	352	cdb->SetDefaultStorage(fOCDBpath);
	353	if (fRun!=cdb->GetRun())
	354	cdb->SetRun(fRun);
	355
	356	AliEMCALGeometry *geometry = AliEMCALGeometry::GetInstance(fGeomName);
	357	if (!geometry) {
	358	AliFatal("Geometry not available!!!");
	359	return;
	360	}
	361
	362	if (!fGeomMatrixSet) {
	363	if (fLoadGeomMatrices) {
	364	for(Int_t mod=0; mod < (geometry->GetEMCGeometry())->GetNumberOfSuperModules(); ++mod) {
	365	if(fGeomMatrix[mod]){
	366	if(DebugLevel() > 2)
	367	fGeomMatrix[mod]->Print();
	368	geometry->SetMisalMatrix(fGeomMatrix[mod],mod);
	369	}
	370	}
	371	} else {
	372	for(Int_t mod=0; mod < geometry->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
	373	if(event->GetEMCALMatrix(mod)) {
	374	if(DebugLevel() > 2)
	375	event->GetEMCALMatrix(mod)->Print();
	376	geometry->SetMisalMatrix(event->GetEMCALMatrix(mod),mod);
	377	}
	378	}
	379	}
	380	fGeomMatrixSet=kTRUE;
	381	}
	382
	383	// setup digit array if needed
	384	if (!fDigitsArr) {
	385	fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
	386	fDigitsArr->SetOwner(1);
	387	}
	388
	389	// then setup clusterizer
	390	delete fClusterizer;
	391	if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
	392	fClusterizer = new AliEMCALClusterizerv1(geometry);
	393	else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
	394	fClusterizer = new AliEMCALClusterizerNxN(geometry);
	395	else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
	396	AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(geometry);
	397	clusterizer->SetNRowDiff(2);
	398	clusterizer->SetNColDiff(2);
	399	fClusterizer = clusterizer;
	400	} else{
	401	AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
	402	}
	403	fClusterizer->InitParameters(fRecParam);
	404	if (!fCalibData&&fLoadCalib) {
	405	AliCDBEntry entry = static_cast<AliCDBEntry>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
	406	if (entry)
	407	fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
	408	if (!fCalibData)
	409	AliFatal("Calibration parameters not found in CDB!");
	410	}
	411	if (!fPedestalData&&fLoadPed) {
	412	AliCDBEntry entry = static_cast<AliCDBEntry>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
	413	if (entry)
	414	fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
	415	}
	416	if (fCalibData) {
	417	fClusterizer->SetInputCalibrated(kFALSE);
	418	fClusterizer->SetCalibrationParameters(fCalibData);
	419	fClusterizer->SetCaloCalibPedestal(fPedestalData);
	420	} else {
	421	fClusterizer->SetInputCalibrated(kTRUE);
	422	}
	423	fClusterizer->SetJustClusters(kTRUE);
	424	fClusterizer->SetDigitsArr(fDigitsArr);
	425	fClusterizer->SetOutput(0);
	426	fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
	427	}