[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),
    fAttachClusters(0),
    fRecalibOnly(0),
    fSubBackground(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),
    fAttachClusters(0),
    fRecalibOnly(0),
    fSubBackground(0)
{ 
  // Constructor

  fBranchNames     = "ESD:AliESDHeader.,AliESDRun.,EMCALCells. AOD:header,emcalCells";
  for(Int_t i = 0; i < 12; ++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");
    return;
  }

  FillDigitsArray();

  if (fRecalibOnly) {
    UpdateCells();
    return; // not requested to run clusterizer
  }

  Clusterize();
  UpdateCells();
  UpdateClusters();

  if (fOutputAODBranch)
    RecPoints2Clusters(fOutputAODBranch);
}

//________________________________________________________________________
void AliAnalysisTaskEMCALClusterizeFast::Clusterize()
{
  // Clusterize

  if (fSubBackground) {
    fClusterizer->SetInputCalibrated(kTRUE);   
    fClusterizer->SetCalibrationParameters(0);
  }
  fClusterizer->Digits2Clusters("");
  if (fSubBackground) {
    fClusterizer->SetInputCalibrated(kFALSE);   
    fClusterizer->SetCalibrationParameters(fCalibData);
  }
}

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

  fDigitsArr->Clear("C");
  AliVCaloCells *cells = InputEvent()->GetEMCALCells();
  Double_t avgE = 0; // for background subtraction
  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->New(idigit));
    digit->SetId(cellNumber);
    digit->SetTime(cellTime);
    digit->SetTimeR(cellTime);
    digit->SetIndexInList(idigit);
    digit->SetType(AliEMCALDigit::kHG);
    if (fRecalibOnly||fSubBackground) {
      Double_t energy = fClusterizer->Calibrate(cellAmplitude,cellTime,cellNumber);
      digit->SetAmplitude(energy);
      avgE += energy;
    } else {
      digit->SetAmplitude(cellAmplitude);
    }
    idigit++;
  }

  if (fSubBackground) {
    avgE /= AliEMCALGeometry::GetInstance(fGeomName)->GetNumberOfSuperModules()*48*24;
    Int_t ndigis = fDigitsArr->GetEntries();
    for (Int_t i = 0; i < ndigis; ++i) {
      AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
      Double_t energy = digit->GetAmplitude() - avgE;
      if (energy<=0.001) {
        //fDigitsArr->RemoveAt(i);
        digit->SetAmplitude(0);
      } else {
        digit->SetAmplitude(energy);
      }
    }
  }
  fDigitsArr->Compress();
  fDigitsArr->Sort();
}

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

  Bool_t esdobjects = 0;
  if (strcmp(clus->GetClass()->GetName(),"AliESDCaloCluster")==0)
    esdobjects = 1;

  Int_t Ncls = fClusterArr->GetEntriesFast();
  for(Int_t i=0, nout=clus->GetEntries(); 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) 
        continue;
      ++ncells_true;
    }
    
    if (ncells_true < 1) {
      AliWarning("Skipping cluster with no cells");
      continue;
    }
    
    // calculate new cluster position
    TVector3 gpos;
    recpoint->GetGlobalPosition(gpos);
    Float_t g[3];
    gpos.GetXYZ(g);
    
    AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
    c->SetType(AliVCluster::kEMCALClusterv1);
    c->SetE(recpoint->GetEnergy());
    c->SetPosition(g);
    c->SetNCells(ncells_true);
    c->SetDispersion(recpoint->GetDispersion());
    c->SetEmcCpvDistance(-1);            //not yet implemented
    c->SetChi2(-1);                      //not yet implemented
    c->SetTOF(recpoint->GetTime()) ;     //time-of-flight
    c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
    Float_t elipAxis[2];
    recpoint->GetElipsAxis(elipAxis);
    c->SetM02(elipAxis[0]*elipAxis[0]) ;
    c->SetM20(elipAxis[1]*elipAxis[1]) ;
    c->SetDistanceToBadChannel(recpoint->GetDistanceToBadTower()); 
    if (esdobjects) {
      AliESDCaloCluster *cesd = static_cast<AliESDCaloCluster*>(c);
      cesd->SetCellsAbsId(absIds);
      cesd->SetCellsAmplitudeFraction(ratios);
    } else {
      AliAODCaloCluster *caod = static_cast<AliAODCaloCluster*>(c);
      caod->SetCellsAbsId(absIds);
      caod->SetCellsAmplitudeFraction(ratios);
    }
  }
 
  AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
  if (!esdevent)
    return;
  if (!fRecoUtils)
    return;

  AliAnalysisManager::GetAnalysisManager()->LoadBranch("Tracks");
  fRecoUtils->FindMatches(esdevent,clus);
  
  if (!esdobjects) {
    Int_t Nclus = clus->GetEntries();
    for(Int_t i=0; i < Nclus; ++i) {
      AliAODCaloCluster *c = static_cast<AliAODCaloCluster*>(clus->At(i));
      Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
      if(trackIndex >= 0) {
        Float_t dR, dZ;
        fRecoUtils->GetMatchedResiduals(i,dR, dZ);
        c->AddTrackMatched(0x0); //esdevent->GetTrack(trackIndex));
        c->SetTrackDistance(dR,dZ); // not implemented
        c->SetEmcCpvDistance(dR);
        c->SetChi2(dZ);
        if(DebugLevel() > 1) 
          AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
      }
    }
  } else {
    Int_t Nclus = clus->GetEntries();
    for(Int_t i=0; i < Nclus; ++i) {
      AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(clus->At(i));
      Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
      if(trackIndex >= 0) {
        Float_t dR, dZ;
        fRecoUtils->GetMatchedResiduals(i,dR, dZ);
        c->SetTrackDistance(dR,dZ);
        c->SetEmcCpvDistance(dR); //to be consistent with AODs
        c->SetChi2(dZ);           //to be consistent with AODs
        TArrayI tm(1,&trackIndex);
        c->AddTracksMatched(tm);
        if(DebugLevel() > 1) 
          AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
      }
    }
  }
}

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

  if (!fCalibData&&!fSubBackground)
    return;

  AliVCaloCells *cells = InputEvent()->GetEMCALCells();
  Int_t ncells = cells->GetNumberOfCells();
  Int_t ndigis = fDigitsArr->GetEntries();

  if (ncells!=ndigis) {
    cells->DeleteContainer();
    cells->CreateContainer(ndigis);
  }
  for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
    AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
    Double_t cellAmplitude = digit->GetCalibAmp();
    Short_t cellNumber = digit->GetId();
    Double_t cellTime = digit->GetTime();
    cells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
  }
}

//________________________________________________________________________
void AliAnalysisTaskEMCALClusterizeFast::UpdateClusters()
{
  // Update cells in case re-calibration was done.

  if (!fAttachClusters)
    return;

  TClonesArray *clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
  if (!clus)
    clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
  if(!clus)
    return;

  Int_t nents = clus->GetEntries();
  for (Int_t i=0;i<nents;++i) {
    AliVCluster *c = static_cast<AliVCluster*>(clus->At(i));
    if (!c)
      continue;
    if (c->IsEMCAL())
      clus->RemoveAt(i);
  }
  clus->Compress();
  RecPoints2Clusters(clus);
}

//________________________________________________________________________________________
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;
  }

  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 { // get matrix from file (work around bug in aliroot)
      for(Int_t mod=0; mod < geometry->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
        const TGeoHMatrix *gm = 0;
        AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(event->GetHeader());
        if (esdevent) {
          gm = esdevent->GetEMCALMatrix(mod);
        } else {
          AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(event->GetHeader());
          if (aodheader) {
            gm = aodheader->GetEMCALMatrix(mod);
          }
        }
        if (gm) {
          if(DebugLevel() > 2) 
            gm->Print();
          geometry->SetMisalMatrix(gm,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) || (!fPedestalData&&fLoadPed)) {
    AliCDBManager *cdb = AliCDBManager::Instance();
    if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
      cdb->SetDefaultStorage(fOCDBpath);
    if (fRun!=cdb->GetRun())
      cdb->SetRun(fRun);
  }

  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
2f7259cf	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
2431b20f	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),
4cf9204b	66	fLoadPed(0),
b0e354fe	67	fAttachClusters(0),
95f69c66	68	fRecalibOnly(0),
95f69c66	69	fSubBackground(0)
2431b20f	70	{
	71	// Constructor
	72	}
	73
2f7259cf	74	//________________________________________________________________________
	75	AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast(const char *name)
	76	: AliAnalysisTaskSE(name),
	77	fRun(0),
	78	fDigitsArr(0),
	79	fClusterArr(0),
	80	fRecParam(new AliEMCALRecParam),
	81	fClusterizer(0),
	82	fUnfolder(0),
	83	fJustUnfold(kFALSE),
	84	fGeomName("EMCAL_FIRSTYEARV1"),
	85	fGeomMatrixSet(kFALSE),
	86	fLoadGeomMatrices(kFALSE),
	87	fOCDBpath(),
	88	fCalibData(0),
	89	fPedestalData(0),
	90	fOutputAODBranch(0),
	91	fOutputAODBrName(),
10d33986	92	fRecoUtils(0),
10d33986	93	fLoadCalib(0),
4cf9204b	94	fLoadPed(0),
b0e354fe	95	fAttachClusters(0),
95f69c66	96	fRecalibOnly(0),
95f69c66	97	fSubBackground(0)
2f7259cf	98	{
	99	// Constructor
	100
01bc1ce8	101	fBranchNames = "ESD:AliESDHeader.,AliESDRun.,EMCALCells. AOD:header,emcalCells";
95f69c66	102	for(Int_t i = 0; i < 12; ++i)
2f7259cf	103	fGeomMatrix[i] = 0;
	104	}
	105
	106	//________________________________________________________________________
	107	AliAnalysisTaskEMCALClusterizeFast::~AliAnalysisTaskEMCALClusterizeFast()
	108	{
	109	// Destructor.
	110
	111	delete fDigitsArr;
	112	delete fClusterizer;
	113	delete fUnfolder;
	114	delete fRecoUtils;
	115	}
	116
	117	//-------------------------------------------------------------------
	118	void AliAnalysisTaskEMCALClusterizeFast::UserCreateOutputObjects()
	119	{
	120	// Create output objects.
	121
	122	if (!fOutputAODBrName.IsNull()) {
	123	fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
	124	fOutputAODBranch->SetName(fOutputAODBrName);
	125	AddAODBranch("TClonesArray", &fOutputAODBranch);
	126	AliInfo(Form("Created Branch: %s",fOutputAODBrName.Data()));
	127	}
	128	}
	129
	130	//________________________________________________________________________
	131	void AliAnalysisTaskEMCALClusterizeFast::UserExec(Option_t *)
	132	{
	133	// Main loop, called for each event
	134
	135	// remove the contents of output list set in the previous event
	136	if (fOutputAODBranch)
	137	fOutputAODBranch->Clear("C");
	138
	139	AliESDEvent esdevent = dynamic_cast<AliESDEvent>(InputEvent());
	140	AliAODEvent aodevent = dynamic_cast<AliAODEvent>(InputEvent());
	141
	142	if (!esdevent&&!aodevent) {
	143	Error("UserExec","Event not available");
	144	return;
	145	}
	146
	147	LoadBranches();
	148
	149	Init();
	150
	151	if (fJustUnfold) {
	152	AliWarning("Unfolding not implemented");
84ccad86	153	return;
	154	}
	155
95f69c66	156	FillDigitsArray();
84ccad86	157
95f69c66	158	if (fRecalibOnly) {
	159	UpdateCells();
	160	return; // not requested to run clusterizer
84ccad86	161	}
7099ec08	162
95f69c66	163	Clusterize();
	164	UpdateCells();
	165	UpdateClusters();
2431b20f	166
95f69c66	167	if (fOutputAODBranch)
95f69c66	168	RecPoints2Clusters(fOutputAODBranch);
2431b20f	169	}
	170
	171	//________________________________________________________________________
95f69c66	172	void AliAnalysisTaskEMCALClusterizeFast::Clusterize()
2431b20f	173	{
95f69c66	174	// Clusterize
7099ec08	175
95f69c66	176	if (fSubBackground) {
	177	fClusterizer->SetInputCalibrated(kTRUE);
	178	fClusterizer->SetCalibrationParameters(0);
2f7259cf	179	}
95f69c66	180	fClusterizer->Digits2Clusters("");
	181	if (fSubBackground) {
	182	fClusterizer->SetInputCalibrated(kFALSE);
	183	fClusterizer->SetCalibrationParameters(fCalibData);
4cf9204b	184	}
4cf9204b	185	}
4cf9204b	186
2f7259cf	187	//________________________________________________________________________
95f69c66	188	void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray()
2f7259cf	189	{
	190	// Fill digits from cells.
	191
	192	fDigitsArr->Clear("C");
95f69c66	193	AliVCaloCells *cells = InputEvent()->GetEMCALCells();
95f69c66	194	Double_t avgE = 0; // for background subtraction
2f7259cf	195	Int_t ncells = cells->GetNumberOfCells();
2f7259cf	196	for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
	197	Double_t cellAmplitude=0, cellTime=0;
	198	Short_t cellNumber=0;
	199	if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
	200	break;
95f69c66	201	AliEMCALDigit digit = static_cast<AliEMCALDigit>(fDigitsArr->New(idigit));
2f7259cf	202	digit->SetId(cellNumber);
2f7259cf	203	digit->SetTime(cellTime);
	204	digit->SetTimeR(cellTime);
	205	digit->SetIndexInList(idigit);
	206	digit->SetType(AliEMCALDigit::kHG);
95f69c66	207	if (fRecalibOnly\|\|fSubBackground) {
b0e354fe	208	Double_t energy = fClusterizer->Calibrate(cellAmplitude,cellTime,cellNumber);
95f69c66	209	digit->SetAmplitude(energy);
	210	avgE += energy;
	211	} else {
	212	digit->SetAmplitude(cellAmplitude);
b0e354fe	213	}
2f7259cf	214	idigit++;
2f7259cf	215	}
2f7259cf	216
95f69c66	217	if (fSubBackground) {
	218	avgE /= AliEMCALGeometry::GetInstance(fGeomName)->GetNumberOfSuperModules()4824;
	219	Int_t ndigis = fDigitsArr->GetEntries();
	220	for (Int_t i = 0; i < ndigis; ++i) {
	221	AliEMCALDigit digit = static_cast<AliEMCALDigit>(fDigitsArr->At(i));
	222	Double_t energy = digit->GetAmplitude() - avgE;
	223	if (energy<=0.001) {
	224	//fDigitsArr->RemoveAt(i);
	225	digit->SetAmplitude(0);
	226	} else {
	227	digit->SetAmplitude(energy);
	228	}
b0e354fe	229	}
2f7259cf	230	}
95f69c66	231	fDigitsArr->Compress();
84ccad86	232	fDigitsArr->Sort();
2f7259cf	233	}
	234
	235	//________________________________________________________________________________________
95f69c66	236	void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters(TClonesArray *clus)
4cf9204b	237	{
	238	// Cluster energy, global position, cells and their amplitude fractions are restored.
	239
95f69c66	240	Bool_t esdobjects = 0;
	241	if (strcmp(clus->GetClass()->GetName(),"AliESDCaloCluster")==0)
	242	esdobjects = 1;
	243
4cf9204b	244	Int_t Ncls = fClusterArr->GetEntriesFast();
	245	for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i) {
	246	AliEMCALRecPoint recpoint = static_cast<AliEMCALRecPoint>(fClusterArr->At(i));
	247	Int_t ncells_true = 0;
	248	const Int_t ncells = recpoint->GetMultiplicity();
	249	UShort_t absIds[ncells];
	250	Double32_t ratios[ncells];
	251	Int_t *dlist = recpoint->GetDigitsList();
	252	Float_t *elist = recpoint->GetEnergiesList();
	253	for (Int_t c = 0; c < ncells; ++c) {
	254	AliEMCALDigit digit = static_cast<AliEMCALDigit>(fDigitsArr->At(dlist[c]));
	255	absIds[ncells_true] = digit->GetId();
	256	ratios[ncells_true] = elist[c]/digit->GetAmplitude();
84ccad86	257	if (ratios[ncells_true] < 0.001)
	258	continue;
	259	++ncells_true;
4cf9204b	260	}
	261
	262	if (ncells_true < 1) {
	263	AliWarning("Skipping cluster with no cells");
	264	continue;
	265	}
	266
	267	// calculate new cluster position
	268	TVector3 gpos;
4cf9204b	269	recpoint->GetGlobalPosition(gpos);
84ccad86	270	Float_t g[3];
4cf9204b	271	gpos.GetXYZ(g);
4cf9204b	272
95f69c66	273	AliVCluster c = static_cast<AliVCluster>(clus->New(nout++));
4cf9204b	274	c->SetType(AliVCluster::kEMCALClusterv1);
	275	c->SetE(recpoint->GetEnergy());
	276	c->SetPosition(g);
	277	c->SetNCells(ncells_true);
4cf9204b	278	c->SetDispersion(recpoint->GetDispersion());
84ccad86	279	c->SetEmcCpvDistance(-1); //not yet implemented
4cf9204b	280	c->SetChi2(-1); //not yet implemented
	281	c->SetTOF(recpoint->GetTime()) ; //time-of-flight
	282	c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
	283	Float_t elipAxis[2];
	284	recpoint->GetElipsAxis(elipAxis);
	285	c->SetM02(elipAxis[0]*elipAxis[0]) ;
	286	c->SetM20(elipAxis[1]*elipAxis[1]) ;
17cff35e	287	c->SetDistanceToBadChannel(recpoint->GetDistanceToBadTower());
95f69c66	288	if (esdobjects) {
	289	AliESDCaloCluster cesd = static_cast<AliESDCaloCluster>(c);
	290	cesd->SetCellsAbsId(absIds);
	291	cesd->SetCellsAmplitudeFraction(ratios);
	292	} else {
	293	AliAODCaloCluster caod = static_cast<AliAODCaloCluster>(c);
	294	caod->SetCellsAbsId(absIds);
	295	caod->SetCellsAmplitudeFraction(ratios);
	296	}
84ccad86	297	}
	298
	299	AliESDEvent esdevent = dynamic_cast<AliESDEvent>(InputEvent());
	300	if (!esdevent)
	301	return;
	302	if (!fRecoUtils)
	303	return;
	304
95f69c66	305	AliAnalysisManager::GetAnalysisManager()->LoadBranch("Tracks");
84ccad86	306	fRecoUtils->FindMatches(esdevent,clus);
84ccad86	307
95f69c66	308	if (!esdobjects) {
	309	Int_t Nclus = clus->GetEntries();
	310	for(Int_t i=0; i < Nclus; ++i) {
	311	AliAODCaloCluster c = static_cast<AliAODCaloCluster>(clus->At(i));
	312	Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
	313	if(trackIndex >= 0) {
	314	Float_t dR, dZ;
	315	fRecoUtils->GetMatchedResiduals(i,dR, dZ);
	316	c->AddTrackMatched(0x0); //esdevent->GetTrack(trackIndex));
	317	c->SetTrackDistance(dR,dZ); // not implemented
	318	c->SetEmcCpvDistance(dR);
	319	c->SetChi2(dZ);
	320	if(DebugLevel() > 1)
	321	AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
	322	}
2f7259cf	323	}
95f69c66	324	} else {
84ccad86	325	Int_t Nclus = clus->GetEntries();
	326	for(Int_t i=0; i < Nclus; ++i) {
	327	AliESDCaloCluster c = static_cast<AliESDCaloCluster>(clus->At(i));
2f7259cf	328	Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
2f7259cf	329	if(trackIndex >= 0) {
84ccad86	330	Float_t dR, dZ;
	331	fRecoUtils->GetMatchedResiduals(i,dR, dZ);
	332	c->SetTrackDistance(dR,dZ);
	333	c->SetEmcCpvDistance(dR); //to be consistent with AODs
	334	c->SetChi2(dZ); //to be consistent with AODs
4cf9204b	335	TArrayI tm(1,&trackIndex);
4cf9204b	336	c->AddTracksMatched(tm);
2f7259cf	337	if(DebugLevel() > 1)
	338	AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
	339	}
	340	}
	341	}
	342	}
	343
95f69c66	344	//________________________________________________________________________
	345	void AliAnalysisTaskEMCALClusterizeFast::UpdateCells()
	346	{
	347	// Update cells in case re-calibration was done.
	348
	349	if (!fCalibData&&!fSubBackground)
	350	return;
	351
	352	AliVCaloCells *cells = InputEvent()->GetEMCALCells();
	353	Int_t ncells = cells->GetNumberOfCells();
	354	Int_t ndigis = fDigitsArr->GetEntries();
	355
	356	if (ncells!=ndigis) {
	357	cells->DeleteContainer();
	358	cells->CreateContainer(ndigis);
	359	}
	360	for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
	361	AliEMCALDigit digit = static_cast<AliEMCALDigit>(fDigitsArr->At(idigit));
	362	Double_t cellAmplitude = digit->GetCalibAmp();
	363	Short_t cellNumber = digit->GetId();
	364	Double_t cellTime = digit->GetTime();
	365	cells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
	366	}
	367	}
	368
	369	//________________________________________________________________________
	370	void AliAnalysisTaskEMCALClusterizeFast::UpdateClusters()
	371	{
	372	// Update cells in case re-calibration was done.
	373
	374	if (!fAttachClusters)
	375	return;
	376
	377	TClonesArray clus = dynamic_cast<TClonesArray>(InputEvent()->FindListObject("caloClusters"));
	378	if (!clus)
	379	clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
	380	if(!clus)
	381	return;
	382
	383	Int_t nents = clus->GetEntries();
	384	for (Int_t i=0;i<nents;++i) {
	385	AliVCluster c = static_cast<AliVCluster>(clus->At(i));
	386	if (!c)
	387	continue;
	388	if (c->IsEMCAL())
	389	clus->RemoveAt(i);
	390	}
	391	clus->Compress();
	392	RecPoints2Clusters(clus);
	393	}
	394
2f7259cf	395	//________________________________________________________________________________________
	396	void AliAnalysisTaskEMCALClusterizeFast::Init()
	397	{
	398	//Select clusterization/unfolding algorithm and set all the needed parameters
01bc1ce8	399
2f7259cf	400	AliVEvent * event = InputEvent();
	401	if (!event) {
	402	AliWarning("Event not available!!!");
	403	return;
	404	}
	405
	406	if (event->GetRunNumber()==fRun)
	407	return;
	408	fRun = event->GetRunNumber();
	409
	410	if (fJustUnfold){
	411	// init the unfolding afterburner
	412	delete fUnfolder;
	413	fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
	414	return;
	415	}
	416
2f7259cf	417	AliEMCALGeometry *geometry = AliEMCALGeometry::GetInstance(fGeomName);
	418	if (!geometry) {
	419	AliFatal("Geometry not available!!!");
	420	return;
	421	}
	422
	423	if (!fGeomMatrixSet) {
	424	if (fLoadGeomMatrices) {
	425	for(Int_t mod=0; mod < (geometry->GetEMCGeometry())->GetNumberOfSuperModules(); ++mod) {
	426	if(fGeomMatrix[mod]){
	427	if(DebugLevel() > 2)
	428	fGeomMatrix[mod]->Print();
	429	geometry->SetMisalMatrix(fGeomMatrix[mod],mod);
	430	}
	431	}
29b496d5	432	} else { // get matrix from file (work around bug in aliroot)
2f7259cf	433	for(Int_t mod=0; mod < geometry->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
29b496d5	434	const TGeoHMatrix *gm = 0;
	435	AliESDEvent esdevent = dynamic_cast<AliESDEvent>(event->GetHeader());
	436	if (esdevent) {
	437	gm = esdevent->GetEMCALMatrix(mod);
	438	} else {
	439	AliAODHeader aodheader = dynamic_cast<AliAODHeader>(event->GetHeader());
	440	if (aodheader) {
	441	gm = aodheader->GetEMCALMatrix(mod);
	442	}
	443	}
	444	if (gm) {
2f7259cf	445	if(DebugLevel() > 2)
29b496d5	446	gm->Print();
29b496d5	447	geometry->SetMisalMatrix(gm,mod);
2f7259cf	448	}
	449	}
	450	}
	451	fGeomMatrixSet=kTRUE;
	452	}
	453
	454	// setup digit array if needed
	455	if (!fDigitsArr) {
	456	fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
	457	fDigitsArr->SetOwner(1);
	458	}
	459
	460	// then setup clusterizer
	461	delete fClusterizer;
	462	if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
	463	fClusterizer = new AliEMCALClusterizerv1(geometry);
	464	else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
	465	fClusterizer = new AliEMCALClusterizerNxN(geometry);
	466	else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
	467	AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(geometry);
	468	clusterizer->SetNRowDiff(2);
	469	clusterizer->SetNColDiff(2);
	470	fClusterizer = clusterizer;
	471	} else{
	472	AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
	473	}
	474	fClusterizer->InitParameters(fRecParam);
01bc1ce8	475
	476	if ((!fCalibData&&fLoadCalib) \|\| (!fPedestalData&&fLoadPed)) {
	477	AliCDBManager *cdb = AliCDBManager::Instance();
	478	if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
	479	cdb->SetDefaultStorage(fOCDBpath);
	480	if (fRun!=cdb->GetRun())
	481	cdb->SetRun(fRun);
	482	}
	483
10d33986	484	if (!fCalibData&&fLoadCalib) {
2f7259cf	485	AliCDBEntry entry = static_cast<AliCDBEntry>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
	486	if (entry)
	487	fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
	488	if (!fCalibData)
	489	AliFatal("Calibration parameters not found in CDB!");
	490	}
10d33986	491	if (!fPedestalData&&fLoadPed) {
2f7259cf	492	AliCDBEntry entry = static_cast<AliCDBEntry>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
	493	if (entry)
	494	fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
	495	}
10d33986	496	if (fCalibData) {
	497	fClusterizer->SetInputCalibrated(kFALSE);
	498	fClusterizer->SetCalibrationParameters(fCalibData);
	499	fClusterizer->SetCaloCalibPedestal(fPedestalData);
	500	} else {
	501	fClusterizer->SetInputCalibrated(kTRUE);
	502	}
	503	fClusterizer->SetJustClusters(kTRUE);
2f7259cf	504	fClusterizer->SetDigitsArr(fDigitsArr);
	505	fClusterizer->SetOutput(0);
	506	fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
	507	}