[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizer.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$ */

//_________________________________________________________________________
//  Base class for the clusterization algorithm (pure abstract)
//*--
//*-- Author: Yves Schutz  SUBATECH 
// 
//   Clusterization mother class. Contains common methods/data members of different 
//   clusterizers. GCB 2010
//////////////////////////////////////////////////////////////////////////////

// --- ROOT system ---
#include "TClonesArray.h"
#include "TTree.h"
#include <TFile.h> 
class TFolder;
#include <TMath.h> 
#include <TMinuit.h>
#include <TTree.h> 
class TSystem; 
#include <TBenchmark.h>
#include <TBrowser.h>
#include <TROOT.h>

// --- Standard library ---
#include <cassert>

// --- AliRoot header files ---
#include "AliEMCALClusterizer.h"
#include "AliEMCALReconstructor.h"
#include "AliRunLoader.h"
#include "AliRun.h"
#include "AliLog.h"
#include "AliEMCAL.h"
#include "AliEMCALRecPoint.h"
#include "AliEMCALRecParam.h"
#include "AliEMCALGeometry.h"
#include "AliEMCALRecParam.h"
#include "AliCDBManager.h"
#include "AliCaloCalibPedestal.h"
#include "AliEMCALCalibData.h"
class AliCDBStorage;
#include "AliCDBEntry.h"

ClassImp(AliEMCALClusterizer)

//____________________________________________________________________________
AliEMCALClusterizer::AliEMCALClusterizer():
  fIsInputCalibrated(kFALSE),
  fJustClusters(kFALSE),
  fDigitsArr(NULL),
  fTreeR(NULL),
  fRecPoints(NULL),
  fGeom(NULL),
  fCalibData(NULL), 
  fCaloPed(NULL),
  fADCchannelECA(0.),fADCpedestalECA(0.), fTimeECA(0.),
  fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
  fDefaultInit(kFALSE),fToUnfold(kFALSE),
  fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
  fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
  fClusterUnfolding(NULL)
{
  // ctor
  
  Init();
}

//____________________________________________________________________________
AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry): 
  fIsInputCalibrated(kFALSE),
  fJustClusters(kFALSE),
  fDigitsArr(NULL),
  fTreeR(NULL),
  fRecPoints(NULL),
  fGeom(geometry),
  fCalibData(NULL), 
  fCaloPed(NULL),
  fADCchannelECA(0.),fADCpedestalECA(0.), fTimeECA(0.),
  fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
  fDefaultInit(kFALSE),fToUnfold(kFALSE),
  fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
  fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
  fClusterUnfolding(NULL)
{
  // Ctor with the indication of the file where header Tree and digits Tree are stored.
  // Use this contructor to avoid usage of Init() which uses runloader.
  // Change needed by HLT - MP.
  // Note for the future: the use on runloader should be avoided or optional at least.
  // Another way is to make Init virtual and protected at least 
  // such that the deriving classes can overload Init();
  
  if (!fGeom)
  {
    AliFatal("Geometry not initialized.");
  }
  Int_t i=0;
  for (i = 0; i < 8; i++)
    fSSPars[i] = 0.;
  for (i = 0; i < 3; i++) {
    fPar5[i] = 0.;
    fPar6[i] = 0.;
  }
}

//____________________________________________________________________________
AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry *geometry, 
                                         AliEMCALCalibData *calib, 
                                         AliCaloCalibPedestal *caloped): 
  fIsInputCalibrated(kFALSE),
  fJustClusters(kFALSE),
  fDigitsArr(NULL),
  fTreeR(NULL),
  fRecPoints(NULL),
  fGeom(geometry),
  fCalibData(calib),
  fCaloPed(caloped),
  fADCchannelECA(0.),fADCpedestalECA(0.), fTimeECA(0.),
  fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
  fDefaultInit(kFALSE),fToUnfold(kFALSE),
  fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
  fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
  fClusterUnfolding(NULL)
{
  // ctor, geometry and calibration are initialized elsewhere.
  
  if (!fGeom)
    AliFatal("Geometry not initialized.");
  
  Int_t i=0;
  for (i = 0; i < 8; i++)
    fSSPars[i] = 0.;
  for (i = 0; i < 3; i++) {
    fPar5[i] = 0.;
    fPar6[i] = 0.;
  }
}

//____________________________________________________________________________
AliEMCALClusterizer::~AliEMCALClusterizer()
{
  // dtor
  //Already deleted in AliEMCALReconstructor.

  if(fClusterUnfolding) delete fClusterUnfolding;

  // make sure we delete the rec points array
  DeleteRecPoints();

  //Delete digits array
  DeleteDigits();

}

//____________________________________________________________________________
void AliEMCALClusterizer::DeleteRecPoints()
{
  // free the cluster array
  if (fRecPoints) 
    {
      AliDebug(2, "Deleting fRecPoints.");
      fRecPoints->Delete();
      delete fRecPoints;
    }
}

//____________________________________________________________________________
void AliEMCALClusterizer::DeleteDigits()
{
  // free the digits array
  if (fDigitsArr) 
    {
      AliDebug(2, "Deleting fDigitsArr.");
      fDigitsArr->Clear("C");
      delete fDigitsArr;
    }
}

//____________________________________________________________________________
void AliEMCALClusterizer::Calibrate(Float_t & amp, Float_t & time, const Int_t absId) 
{
  // Convert digitized amplitude into energy, calibrate time
  // Calibration parameters are taken from OCDB : OCDB/EMCAL/Calib/Data

  //Check if time is too large or too small, indication of a noisy channel, remove in this case
  if(time > fTimeMax || time < fTimeMin) {
    amp  = 0 ;
    time = 0 ;
    return ;
  }  
  
  //Return energy with default parameters if calibration is not available
  if (!fCalibData && !fCaloPed) {
    if (fIsInputCalibrated == kTRUE)
    {
      AliDebug(10, Form("Input already calibrated!"));
      return ;
    }
    else{
      AliFatal("OCDB calibration and bad map parameters are not available");
      return;
    }   
  }
  
  if (fGeom==0)
    AliFatal("Did not get geometry from EMCALLoader") ;
    
  Int_t iSupMod = -1;
  Int_t nModule = -1;
  Int_t nIphi   = -1;
  Int_t nIeta   = -1;
  Int_t iphi    = -1;
  Int_t ieta    = -1;
    
  Bool_t bCell = fGeom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
  if(!bCell) {
    fGeom->PrintGeometry();
    AliError(Form("Wrong cell id number : %i", absId));
    //assert(0); // GCB: This aborts reconstruction of raw simulations 
    //where simulation had more SM than default geometry, 
    //change to return 0, to avoid aborting good generations.
    amp  = 0;
    time = 0;
    return ;
  }
    
  fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
	  
  // Check if channel is bad (dead or hot), in this case return 0.	
  // Gustavo: 15-12-09 In case of RAW data this selection is already done, but not in simulation.
  // for the moment keep it here but remember to do the selection at the sdigitizer level 
  // and remove it from here
  if (fCaloPed) {
    Int_t channelStatus = (Int_t)(fCaloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
    if(channelStatus == AliCaloCalibPedestal::kHot || channelStatus == AliCaloCalibPedestal::kDead) {
      AliDebug(2,Form("Tower from SM %d, ieta %d, iphi %d is BAD : status %d !!!",iSupMod,ieta,iphi, channelStatus));
      amp  = 0 ;
      time = 0 ;
      return ;
    }
  }
    
  if (fIsInputCalibrated || !fCalibData)
  {
    AliDebug(10, Form("Input already calibrated!"));
    return ;
  }
	  
  fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
  fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
  fTimeECA        = fCalibData->GetTimeChannel(iSupMod,ieta,iphi);
  
  time -= fTimeECA ;
  amp   = amp * fADCchannelECA - fADCpedestalECA ;  
  
}

//____________________________________________________________________________
void AliEMCALClusterizer::GetCalibrationParameters() 
{
  // Set calibration parameters:
  // If calibration database exists, they are read from database,
  // otherwise, they are taken from digitizer.
  // It is a user responsilibity to open CDB before reconstruction, 
  // for example: 
  // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");

  if (fIsInputCalibrated)
    return;
  
  //Check if calibration is stored in data base
  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!");
}

//____________________________________________________________________________
void AliEMCALClusterizer::GetCaloCalibPedestal() 
{
  // Set calibration parameters:
  // if calibration database exists, they are read from database,
  // otherwise, they are taken from digitizer.
  //
  // It is a user responsilibity to open CDB before reconstruction, 
  // for example: 
  // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");

  if (fIsInputCalibrated)
    return;
  
  // Check if calibration is stored in data base
  if(!fCaloPed)
    {
      AliCDBEntry *entry = (AliCDBEntry*) 
	AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
      if (entry) fCaloPed =  (AliCaloCalibPedestal*) entry->GetObject();
    }
  
  if(!fCaloPed)
    AliFatal("Pedestal info not found in CDB!");
}

//____________________________________________________________________________
void AliEMCALClusterizer::Init()
{
  // Make all memory allocations which can not be done in default constructor.
  // Attach the Clusterizer task to the list of EMCAL tasks
  
  AliRunLoader *rl = AliRunLoader::Instance();
  if (rl->GetAliRun()){
    AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
    if(emcal)fGeom = emcal->GetGeometry();
  }
  
  if(!fGeom){ 
    fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
  }
  
  AliDebug(1,Form("geom %p",fGeom));
  
  if(!gMinuit) 
    gMinuit = new TMinuit(100) ;
  
  Int_t i=0;
  for (i = 0; i < 8; i++)
    fSSPars[i] = 0.;
  for (i = 0; i < 3; i++) {
    fPar5[i] = 0.;
    fPar6[i] = 0.;
  }
}

//____________________________________________________________________________
void AliEMCALClusterizer::InitParameters()
{ 
  // Initializes the parameters for the Clusterizer from AliEMCALReconstructor::GetRecParam().

  return InitParameters(AliEMCALReconstructor::GetRecParam());
}

//____________________________________________________________________________
void AliEMCALClusterizer::InitParameters(const AliEMCALRecParam* recParam)
{ 
  // Initializes the parameters for the Clusterizer

  fNumberOfECAClusters = 0 ;
  fCalibData           = 0 ;
  fCaloPed             = 0 ;
	
  if(!recParam) {
    AliFatal("Reconstruction parameters for EMCAL not set!");
  } 

  fECAClusteringThreshold = recParam->GetClusteringThreshold();
  fECAW0                  = recParam->GetW0();
  fMinECut                = recParam->GetMinECut();    
  fToUnfold               = recParam->GetUnfold();
  fECALocMaxCut           = recParam->GetLocMaxCut();
  fTimeCut                = recParam->GetTimeCut();
  fTimeMin                = recParam->GetTimeMin();
  fTimeMax                = recParam->GetTimeMax();
  
  //For NxN
  SetNRowDiff(recParam->GetNRowDiff());
  SetNColDiff(recParam->GetNColDiff());
  
  AliDebug(1,Form("Reconstruction parameters: fECAClusteringThreshold=%.3f GeV, fECAW=%.3f, fMinECut=%.3f GeV, "
                  "fToUnfold=%d, fECALocMaxCut=%.3f GeV, fTimeCut=%e s,fTimeMin=%e s,fTimeMax=%e s",
                  fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax));

  if (fToUnfold) {
    Int_t i=0;
    for (i = 0; i < 8; i++) {
      fSSPars[i] = recParam->GetSSPars(i);
    } //end of loop over parameters
    for (i = 0; i < 3; i++) {
      fPar5[i] = recParam->GetPar5(i);
      fPar6[i] = recParam->GetPar6(i);
    } //end of loop over parameters
      
    InitClusterUnfolding();
      
    for (i = 0; i < 8; i++) {
      AliDebug(1,Form("unfolding shower shape parameters: fSSPars=%f \n",fSSPars[i]));
    }
    for (i = 0; i < 3; i++) {
      AliDebug(1,Form("unfolding parameter 5: fPar5=%f \n",fPar5[i]));
      AliDebug(1,Form("unfolding parameter 6: fPar6=%f \n",fPar6[i]));
    }
  } // to unfold
}

//____________________________________________________________________________
void AliEMCALClusterizer::Print(Option_t * /*option*/)const
{
  // Print clusterizer parameters
  
  TString message("\n") ; 
  
  if (strcmp(GetName(),"") == 0) {
    printf("AliEMCALClusterizer not initialized\n");
    return;
  }
    
  // Print parameters
  TString taskName(Version()) ;
    
  printf("--------------- "); 
  printf("%s",taskName.Data()) ; 
  printf(" "); 
  printf("Clusterizing digits: "); 
  printf("\n                       ECA Local Maximum cut    = %f", fECALocMaxCut); 
  printf("\n                       ECA Logarithmic weight   = %f", fECAW0); 
  if (fToUnfold) {
    printf("\nUnfolding on\n");
    printf("Unfolding parameters: fSSpars: \n");
    Int_t i=0;
    for (i = 0; i < 8; i++) {
      printf("fSSPars[%d] = %f \n", i, fSSPars[i]);
    }
    printf("Unfolding parameter 5 and 6: fPar5 and fPar6: \n");
    for (i = 0; i < 3; i++) {
      printf("fPar5[%d] = %f \n", i, fPar5[i]);
      printf("fPar6[%d] = %f \n", i, fPar6[i]);
    }
  }
  else
    printf("\nUnfolding off\n");
    
  printf("------------------------------------------------------------------"); 
}

//____________________________________________________________________________
void AliEMCALClusterizer::PrintRecPoints(Option_t * option)
{
  // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer

  if (strstr(option,"deb")) {
    printf("PrintRecPoints: Clusterization result:") ; 
    printf("           Found %d ECA Rec Points\n ", 
           fRecPoints->GetEntriesFast()) ; 
  }
  
  if (strstr(option,"all")) {
    if (strstr(option,"deb")) {
      printf("\n-----------------------------------------------------------------------\n") ;
      printf("Clusters in ECAL section\n") ;
      printf("Index    Ene(GeV) Multi Module     GX    GY   GZ  lX    lY   lZ   Dispersion Lambda 1   Lambda 2  # of prim  Primaries list\n") ;
    }
    Int_t index; 
    for (index =  0 ; index < fRecPoints->GetEntries() ; index++) {
      AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(fRecPoints->At(index)) ; 
      if (!rp) 
        continue;

      TVector3  globalpos;  
      //rp->GetGlobalPosition(globalpos);
      TVector3  localpos;  
      rp->GetLocalPosition(localpos);
      Float_t lambda[2]; 
      rp->GetElipsAxis(lambda);
        
      Int_t nprimaries=0;
      Int_t * primaries = rp->GetPrimaries(nprimaries);
      if(strstr(option,"deb")) 
        printf("\n%6d  %8.4f  %3d     %4.1f    %4.1f %4.1f  %4.1f %4.1f %4.1f    %4.1f   %4f  %4f    %2d     : ", 
               rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
               globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), 
               rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ; 
      if(strstr(option,"deb")){ 
        for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
          printf("%d ", primaries[iprimary] ) ; 
        }
      }
    }
    
    if(strstr(option,"deb"))
      printf("\n-----------------------------------------------------------------------\n");
  }
}

//___________________________________________________________________
void  AliEMCALClusterizer::PrintRecoInfo()
{
  // Print reco version

  printf(" AliEMCALClusterizer::PrintRecoInfo() : version %s \n", Version() );
}

//____________________________________________________________________________
void AliEMCALClusterizer::SetInput(TTree *digitsTree)
{
  // Read the digits from the input tree

  TBranch *branch = digitsTree->GetBranch("EMCAL");
  if (!branch) { 
    AliError("can't get the branch with the EMCAL digits !");
    return;
  }
  if (!fDigitsArr)
    fDigitsArr = new TClonesArray("AliEMCALDigit",100);
  branch->SetAddress(&fDigitsArr);
  branch->GetEntry(0);
}

//____________________________________________________________________________
void AliEMCALClusterizer::SetOutput(TTree *clustersTree)
{
  // Read the digits from the input tree

  AliDebug(9, "Making array for EMCAL clusters");
  fRecPoints = new TObjArray(1000);
  if (clustersTree) {
    fTreeR = clustersTree;
    Int_t split = 0;
    Int_t bufsize = 32000;
    fTreeR->Branch("EMCALECARP", "TObjArray", &fRecPoints, bufsize, split);
  }
}

//___________________________________________________________________
void AliEMCALClusterizer::SetInputCalibrated(Bool_t val)
{
  // Flag to indicate that input is calibrated - the case when we run already on ESD

  fIsInputCalibrated = val;
}

//___________________________________________________________________
void AliEMCALClusterizer::SetJustClusters(Bool_t val)
{
  // Flag to indicate that we are running on ESDs, when calling 
  // rp->EvalAll(fECAW0,fDigitsArr,fJustClusters); in derived classes

  fJustClusters = val;
}
Commit	Line	Data
483b0559	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	//_________________________________________________________________________
	19	// Base class for the clusterization algorithm (pure abstract)
	20	//*--
	21	//*-- Author: Yves Schutz SUBATECH
ee08edde	22	//
	23	// Clusterization mother class. Contains common methods/data members of different
	24	// clusterizers. GCB 2010
483b0559	25	//////////////////////////////////////////////////////////////////////////////
	26
	27	// --- ROOT system ---
c47157cd	28	#include "TClonesArray.h"
c47157cd	29	#include "TTree.h"
ee08edde	30	#include <TFile.h>
	31	class TFolder;
	32	#include <TMath.h>
	33	#include <TMinuit.h>
	34	#include <TTree.h>
	35	class TSystem;
	36	#include <TBenchmark.h>
	37	#include <TBrowser.h>
	38	#include <TROOT.h>
483b0559	39
483b0559	40	// --- Standard library ---
ee08edde	41	#include <cassert>
483b0559	42
483b0559	43	// --- AliRoot header files ---
483b0559	44	#include "AliEMCALClusterizer.h"
ee08edde	45	#include "AliEMCALReconstructor.h"
	46	#include "AliRunLoader.h"
	47	#include "AliRun.h"
c47157cd	48	#include "AliLog.h"
ee08edde	49	#include "AliEMCAL.h"
	50	#include "AliEMCALRecPoint.h"
	51	#include "AliEMCALRecParam.h"
	52	#include "AliEMCALGeometry.h"
	53	#include "AliEMCALRecParam.h"
	54	#include "AliCDBManager.h"
	55	#include "AliCaloCalibPedestal.h"
	56	#include "AliEMCALCalibData.h"
	57	class AliCDBStorage;
	58	#include "AliCDBEntry.h"
483b0559	59
	60	ClassImp(AliEMCALClusterizer)
	61
	62	//____________________________________________________________________________
c47157cd	63	AliEMCALClusterizer::AliEMCALClusterizer():
ac084be7	64	fIsInputCalibrated(kFALSE),
294553d1	65	fJustClusters(kFALSE),
c47157cd	66	fDigitsArr(NULL),
c47157cd	67	fTreeR(NULL),
ee08edde	68	fRecPoints(NULL),
	69	fGeom(NULL),
	70	fCalibData(NULL),
	71	fCaloPed(NULL),
e93968a8	72	fADCchannelECA(0.),fADCpedestalECA(0.), fTimeECA(0.),
ee08edde	73	fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
	74	fDefaultInit(kFALSE),fToUnfold(kFALSE),
	75	fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
65bec413	76	fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
65bec413	77	fClusterUnfolding(NULL)
483b0559	78	{
483b0559	79	// ctor
ee08edde	80
ee08edde	81	Init();
ee08edde	82	}
	83
	84	//____________________________________________________________________________
	85	AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry):
ac084be7	86	fIsInputCalibrated(kFALSE),
294553d1	87	fJustClusters(kFALSE),
ee08edde	88	fDigitsArr(NULL),
	89	fTreeR(NULL),
	90	fRecPoints(NULL),
	91	fGeom(geometry),
	92	fCalibData(NULL),
	93	fCaloPed(NULL),
e93968a8	94	fADCchannelECA(0.),fADCpedestalECA(0.), fTimeECA(0.),
ee08edde	95	fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
	96	fDefaultInit(kFALSE),fToUnfold(kFALSE),
	97	fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
65bec413	98	fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
65bec413	99	fClusterUnfolding(NULL)
ee08edde	100	{
ac084be7	101	// Ctor with the indication of the file where header Tree and digits Tree are stored.
	102	// Use this contructor to avoid usage of Init() which uses runloader.
	103	// Change needed by HLT - MP.
	104	// Note for the future: the use on runloader should be avoided or optional at least.
	105	// Another way is to make Init virtual and protected at least
	106	// such that the deriving classes can overload Init();
ee08edde	107
	108	if (!fGeom)
	109	{
	110	AliFatal("Geometry not initialized.");
	111	}
65bec413	112	Int_t i=0;
	113	for (i = 0; i < 8; i++)
	114	fSSPars[i] = 0.;
	115	for (i = 0; i < 3; i++) {
	116	fPar5[i] = 0.;
	117	fPar6[i] = 0.;
	118	}
483b0559	119	}
839828a6	120
ee08edde	121	//____________________________________________________________________________
ac084be7	122	AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry *geometry,
	123	AliEMCALCalibData *calib,
	124	AliCaloCalibPedestal *caloped):
	125	fIsInputCalibrated(kFALSE),
294553d1	126	fJustClusters(kFALSE),
ee08edde	127	fDigitsArr(NULL),
	128	fTreeR(NULL),
	129	fRecPoints(NULL),
	130	fGeom(geometry),
	131	fCalibData(calib),
	132	fCaloPed(caloped),
e93968a8	133	fADCchannelECA(0.),fADCpedestalECA(0.), fTimeECA(0.),
ee08edde	134	fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
	135	fDefaultInit(kFALSE),fToUnfold(kFALSE),
	136	fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
65bec413	137	fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
65bec413	138	fClusterUnfolding(NULL)
ee08edde	139	{
c526514b	140	// ctor, geometry and calibration are initialized elsewhere.
	141
	142	if (!fGeom)
	143	AliFatal("Geometry not initialized.");
ee08edde	144
65bec413	145	Int_t i=0;
	146	for (i = 0; i < 8; i++)
	147	fSSPars[i] = 0.;
	148	for (i = 0; i < 3; i++) {
	149	fPar5[i] = 0.;
	150	fPar6[i] = 0.;
	151	}
ee08edde	152	}
ee08edde	153
483b0559	154	//____________________________________________________________________________
c47157cd	155	AliEMCALClusterizer::~AliEMCALClusterizer()
483b0559	156	{
c47157cd	157	// dtor
23ca956b	158	//Already deleted in AliEMCALReconstructor.
23ca956b	159
65bec413	160	if(fClusterUnfolding) delete fClusterUnfolding;
c828bc97	161
	162	// make sure we delete the rec points array
	163	DeleteRecPoints();
8d0210ea	164
	165	//Delete digits array
	166	DeleteDigits();
	167
c828bc97	168	}
	169
	170	//____________________________________________________________________________
	171	void AliEMCALClusterizer::DeleteRecPoints()
	172	{
	173	// free the cluster array
	174	if (fRecPoints)
	175	{
	176	AliDebug(2, "Deleting fRecPoints.");
	177	fRecPoints->Delete();
	178	delete fRecPoints;
8d0210ea	179	}
	180	}
	181
	182	//____________________________________________________________________________
	183	void AliEMCALClusterizer::DeleteDigits()
	184	{
	185	// free the digits array
	186	if (fDigitsArr)
	187	{
	188	AliDebug(2, "Deleting fDigitsArr.");
	189	fDigitsArr->Clear("C");
	190	delete fDigitsArr;
c828bc97	191	}
18a21c7c	192	}
18a21c7c	193
ee08edde	194	//____________________________________________________________________________
783153ff	195	void AliEMCALClusterizer::Calibrate(Float_t & amp, Float_t & time, const Int_t absId)
ee08edde	196	{
783153ff	197	// Convert digitized amplitude into energy, calibrate time
783153ff	198	// Calibration parameters are taken from OCDB : OCDB/EMCAL/Calib/Data
ac084be7	199
db8df68e	200	//Check if time is too large or too small, indication of a noisy channel, remove in this case
	201	if(time > fTimeMax \|\| time < fTimeMin) {
	202	amp = 0 ;
	203	time = 0 ;
	204	return ;
	205	}
	206
ac084be7	207	//Return energy with default parameters if calibration is not available
e93968a8	208	if (!fCalibData && !fCaloPed) {
ac084be7	209	if (fIsInputCalibrated == kTRUE)
b1324a01	210	{
b1324a01	211	AliDebug(10, Form("Input already calibrated!"));
783153ff	212	return ;
e93968a8	213	}
	214	else{
	215	AliFatal("OCDB calibration and bad map parameters are not available");
	216	return;
	217	}
ac084be7	218	}
	219
	220	if (fGeom==0)
	221	AliFatal("Did not get geometry from EMCALLoader") ;
b1324a01	222
ac084be7	223	Int_t iSupMod = -1;
	224	Int_t nModule = -1;
	225	Int_t nIphi = -1;
	226	Int_t nIeta = -1;
	227	Int_t iphi = -1;
	228	Int_t ieta = -1;
ee08edde	229
ac084be7	230	Bool_t bCell = fGeom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
	231	if(!bCell) {
	232	fGeom->PrintGeometry();
	233	AliError(Form("Wrong cell id number : %i", absId));
	234	//assert(0); // GCB: This aborts reconstruction of raw simulations
	235	//where simulation had more SM than default geometry,
	236	//change to return 0, to avoid aborting good generations.
783153ff	237	amp = 0;
	238	time = 0;
	239	return ;
ee08edde	240	}
b1324a01	241
ac084be7	242	fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
	243
	244	// Check if channel is bad (dead or hot), in this case return 0.
	245	// Gustavo: 15-12-09 In case of RAW data this selection is already done, but not in simulation.
	246	// for the moment keep it here but remember to do the selection at the sdigitizer level
	247	// and remove it from here
294553d1	248	if (fCaloPed) {
	249	Int_t channelStatus = (Int_t)(fCaloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
	250	if(channelStatus == AliCaloCalibPedestal::kHot \|\| channelStatus == AliCaloCalibPedestal::kDead) {
	251	AliDebug(2,Form("Tower from SM %d, ieta %d, iphi %d is BAD : status %d !!!",iSupMod,ieta,iphi, channelStatus));
783153ff	252	amp = 0 ;
	253	time = 0 ;
	254	return ;
294553d1	255	}
ac084be7	256	}
b1324a01	257
db8df68e	258	if (fIsInputCalibrated \|\| !fCalibData)
ac084be7	259	{
ac084be7	260	AliDebug(10, Form("Input already calibrated!"));
783153ff	261	return ;
b1324a01	262	}
ac084be7	263
	264	fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
	265	fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
783153ff	266	fTimeECA = fCalibData->GetTimeChannel(iSupMod,ieta,iphi);
	267
	268	time -= fTimeECA ;
	269	amp = amp * fADCchannelECA - fADCpedestalECA ;
	270
ee08edde	271	}
	272
	273	//____________________________________________________________________________
	274	void AliEMCALClusterizer::GetCalibrationParameters()
	275	{
	276	// Set calibration parameters:
ac084be7	277	// If calibration database exists, they are read from database,
ee08edde	278	// otherwise, they are taken from digitizer.
ee08edde	279	// It is a user responsilibity to open CDB before reconstruction,
	280	// for example:
	281	// AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
ac084be7	282
	283	if (fIsInputCalibrated)
	284	return;
ee08edde	285
ee08edde	286	//Check if calibration is stored in data base
ee08edde	287	if(!fCalibData)
	288	{
	289	AliCDBEntry entry = (AliCDBEntry)
	290	AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
	291	if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
	292	}
	293
	294	if(!fCalibData)
	295	AliFatal("Calibration parameters not found in CDB!");
ee08edde	296	}
	297
	298	//____________________________________________________________________________
	299	void AliEMCALClusterizer::GetCaloCalibPedestal()
	300	{
c526514b	301	// Set calibration parameters:
	302	// if calibration database exists, they are read from database,
	303	// otherwise, they are taken from digitizer.
	304	//
	305	// It is a user responsilibity to open CDB before reconstruction,
	306	// for example:
	307	// AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
ac084be7	308
	309	if (fIsInputCalibrated)
	310	return;
c526514b	311
ac084be7	312	// Check if calibration is stored in data base
c526514b	313	if(!fCaloPed)
	314	{
	315	AliCDBEntry entry = (AliCDBEntry)
	316	AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
	317	if (entry) fCaloPed = (AliCaloCalibPedestal*) entry->GetObject();
	318	}
	319
	320	if(!fCaloPed)
	321	AliFatal("Pedestal info not found in CDB!");
ee08edde	322	}
	323
	324	//____________________________________________________________________________
	325	void AliEMCALClusterizer::Init()
	326	{
	327	// Make all memory allocations which can not be done in default constructor.
	328	// Attach the Clusterizer task to the list of EMCAL tasks
	329
	330	AliRunLoader *rl = AliRunLoader::Instance();
a51e676d	331	if (rl->GetAliRun()){
	332	AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
	333	if(emcal)fGeom = emcal->GetGeometry();
	334	}
	335
	336	if(!fGeom){
ac084be7	337	fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
a51e676d	338	}
ee08edde	339
	340	AliDebug(1,Form("geom %p",fGeom));
	341
	342	if(!gMinuit)
	343	gMinuit = new TMinuit(100) ;
	344
65bec413	345	Int_t i=0;
	346	for (i = 0; i < 8; i++)
	347	fSSPars[i] = 0.;
	348	for (i = 0; i < 3; i++) {
	349	fPar5[i] = 0.;
	350	fPar6[i] = 0.;
	351	}
ee08edde	352	}
	353
	354	//____________________________________________________________________________
	355	void AliEMCALClusterizer::InitParameters()
ac084be7	356	{
	357	// Initializes the parameters for the Clusterizer from AliEMCALReconstructor::GetRecParam().
	358
	359	return InitParameters(AliEMCALReconstructor::GetRecParam());
	360	}
	361
	362	//____________________________________________________________________________
	363	void AliEMCALClusterizer::InitParameters(const AliEMCALRecParam* recParam)
ee08edde	364	{
ee08edde	365	// Initializes the parameters for the Clusterizer
ac084be7	366
ee08edde	367	fNumberOfECAClusters = 0 ;
	368	fCalibData = 0 ;
	369	fCaloPed = 0 ;
	370
ee08edde	371	if(!recParam) {
ee08edde	372	AliFatal("Reconstruction parameters for EMCAL not set!");
ac084be7	373	}
	374
	375	fECAClusteringThreshold = recParam->GetClusteringThreshold();
	376	fECAW0 = recParam->GetW0();
	377	fMinECut = recParam->GetMinECut();
	378	fToUnfold = recParam->GetUnfold();
	379	fECALocMaxCut = recParam->GetLocMaxCut();
	380	fTimeCut = recParam->GetTimeCut();
	381	fTimeMin = recParam->GetTimeMin();
	382	fTimeMax = recParam->GetTimeMax();
d464daf4	383
	384	//For NxN
	385	SetNRowDiff(recParam->GetNRowDiff());
	386	SetNColDiff(recParam->GetNColDiff());
	387
ac084be7	388	AliDebug(1,Form("Reconstruction parameters: fECAClusteringThreshold=%.3f GeV, fECAW=%.3f, fMinECut=%.3f GeV, "
	389	"fToUnfold=%d, fECALocMaxCut=%.3f GeV, fTimeCut=%e s,fTimeMin=%e s,fTimeMax=%e s",
	390	fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax));
	391
	392	if (fToUnfold) {
	393	Int_t i=0;
	394	for (i = 0; i < 8; i++) {
	395	fSSPars[i] = recParam->GetSSPars(i);
	396	} //end of loop over parameters
	397	for (i = 0; i < 3; i++) {
	398	fPar5[i] = recParam->GetPar5(i);
	399	fPar6[i] = recParam->GetPar6(i);
	400	} //end of loop over parameters
622e10be	401
ac084be7	402	InitClusterUnfolding();
622e10be	403
ac084be7	404	for (i = 0; i < 8; i++) {
	405	AliDebug(1,Form("unfolding shower shape parameters: fSSPars=%f \n",fSSPars[i]));
	406	}
	407	for (i = 0; i < 3; i++) {
	408	AliDebug(1,Form("unfolding parameter 5: fPar5=%f \n",fPar5[i]));
	409	AliDebug(1,Form("unfolding parameter 6: fPar6=%f \n",fPar6[i]));
	410	}
	411	} // to unfold
ee08edde	412	}
ee08edde	413
ee08edde	414	//____________________________________________________________________________
	415	void AliEMCALClusterizer::Print(Option_t * /option/)const
	416	{
	417	// Print clusterizer parameters
	418
	419	TString message("\n") ;
	420
ac084be7	421	if (strcmp(GetName(),"") == 0) {
	422	printf("AliEMCALClusterizer not initialized\n");
	423	return;
	424	}
ee08edde	425
ac084be7	426	// Print parameters
ac084be7	427	TString taskName(Version()) ;
ee08edde	428
ac084be7	429	printf("--------------- ");
	430	printf("%s",taskName.Data()) ;
	431	printf(" ");
	432	printf("Clusterizing digits: ");
	433	printf("\n ECA Local Maximum cut = %f", fECALocMaxCut);
	434	printf("\n ECA Logarithmic weight = %f", fECAW0);
	435	if (fToUnfold) {
	436	printf("\nUnfolding on\n");
	437	printf("Unfolding parameters: fSSpars: \n");
	438	Int_t i=0;
	439	for (i = 0; i < 8; i++) {
	440	printf("fSSPars[%d] = %f \n", i, fSSPars[i]);
	441	}
	442	printf("Unfolding parameter 5 and 6: fPar5 and fPar6: \n");
	443	for (i = 0; i < 3; i++) {
	444	printf("fPar5[%d] = %f \n", i, fPar5[i]);
	445	printf("fPar6[%d] = %f \n", i, fPar6[i]);
c526514b	446	}
ee08edde	447	}
ee08edde	448	else
ac084be7	449	printf("\nUnfolding off\n");
	450
	451	printf("------------------------------------------------------------------");
ee08edde	452	}
	453
	454	//____________________________________________________________________________
	455	void AliEMCALClusterizer::PrintRecPoints(Option_t * option)
	456	{
	457	// Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
ac084be7	458
ac084be7	459	if (strstr(option,"deb")) {
ee08edde	460	printf("PrintRecPoints: Clusterization result:") ;
ee08edde	461	printf(" Found %d ECA Rec Points\n ",
	462	fRecPoints->GetEntriesFast()) ;
	463	}
	464
ac084be7	465	if (strstr(option,"all")) {
ac084be7	466	if (strstr(option,"deb")) {
ee08edde	467	printf("\n-----------------------------------------------------------------------\n") ;
	468	printf("Clusters in ECAL section\n") ;
	469	printf("Index Ene(GeV) Multi Module GX GY GZ lX lY lZ Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ;
	470	}
	471	Int_t index;
	472	for (index = 0 ; index < fRecPoints->GetEntries() ; index++) {
	473	AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(fRecPoints->At(index)) ;
ac084be7	474	if (!rp)
	475	continue;
	476
	477	TVector3 globalpos;
	478	//rp->GetGlobalPosition(globalpos);
	479	TVector3 localpos;
	480	rp->GetLocalPosition(localpos);
	481	Float_t lambda[2];
	482	rp->GetElipsAxis(lambda);
51e4198e	483
ac084be7	484	Int_t nprimaries=0;
	485	Int_t * primaries = rp->GetPrimaries(nprimaries);
	486	if(strstr(option,"deb"))
	487	printf("\n%6d %8.4f %3d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ",
	488	rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
	489	globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
	490	rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
	491	if(strstr(option,"deb")){
	492	for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
	493	printf("%d ", primaries[iprimary] ) ;
ee08edde	494	}
	495	}
	496	}
	497
	498	if(strstr(option,"deb"))
	499	printf("\n-----------------------------------------------------------------------\n");
	500	}
	501	}
	502
	503	//___________________________________________________________________
	504	void AliEMCALClusterizer::PrintRecoInfo()
	505	{
ac084be7	506	// Print reco version
ac084be7	507
ee08edde	508	printf(" AliEMCALClusterizer::PrintRecoInfo() : version %s \n", Version() );
ee08edde	509	}
ee08edde	510
18a21c7c	511	//____________________________________________________________________________
c47157cd	512	void AliEMCALClusterizer::SetInput(TTree *digitsTree)
18a21c7c	513	{
c47157cd	514	// Read the digits from the input tree
ac084be7	515
c47157cd	516	TBranch *branch = digitsTree->GetBranch("EMCAL");
	517	if (!branch) {
	518	AliError("can't get the branch with the EMCAL digits !");
	519	return;
	520	}
3f9ad99f	521	if (!fDigitsArr)
3f9ad99f	522	fDigitsArr = new TClonesArray("AliEMCALDigit",100);
c47157cd	523	branch->SetAddress(&fDigitsArr);
c47157cd	524	branch->GetEntry(0);
483b0559	525	}
	526
	527	//____________________________________________________________________________
c47157cd	528	void AliEMCALClusterizer::SetOutput(TTree *clustersTree)
483b0559	529	{
c47157cd	530	// Read the digits from the input tree
ac084be7	531
c47157cd	532	AliDebug(9, "Making array for EMCAL clusters");
ac084be7	533	fRecPoints = new TObjArray(1000);
	534	if (clustersTree) {
	535	fTreeR = clustersTree;
	536	Int_t split = 0;
	537	Int_t bufsize = 32000;
	538	fTreeR->Branch("EMCALECARP", "TObjArray", &fRecPoints, bufsize, split);
	539	}
88cb7938	540	}
ee08edde	541
b1324a01	542	//___________________________________________________________________
ac084be7	543	void AliEMCALClusterizer::SetInputCalibrated(Bool_t val)
b1324a01	544	{
ac084be7	545	// Flag to indicate that input is calibrated - the case when we run already on ESD
ee08edde	546
ac084be7	547	fIsInputCalibrated = val;
ac084be7	548	}
294553d1	549
	550	//___________________________________________________________________
	551	void AliEMCALClusterizer::SetJustClusters(Bool_t val)
	552	{
	553	// Flag to indicate that we are running on ESDs, when calling
	554	// rp->EvalAll(fECAW0,fDigitsArr,fJustClusters); in derived classes
	555
	556	fJustClusters = val;
	557	}