[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)

Bool_t AliEMCALClusterizer::fgkIsInputCalibrated = kFALSE;


//____________________________________________________________________________
AliEMCALClusterizer::AliEMCALClusterizer():
  fDigitsArr(NULL),
  fTreeR(NULL),
  fRecPoints(NULL),
  fGeom(NULL),
  fCalibData(NULL), 
  fCaloPed(NULL),
  fADCchannelECA(0.),fADCpedestalECA(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): 
  fDigitsArr(NULL),
  fTreeR(NULL),
  fRecPoints(NULL),
  fGeom(geometry),
  fCalibData(NULL), 
  fCaloPed(NULL),
  fADCchannelECA(0.),fADCpedestalECA(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): 
  fDigitsArr(NULL),
  fTreeR(NULL),
  fRecPoints(NULL),
  fGeom(geometry),
  fCalibData(calib),
  fCaloPed(caloped),
  fADCchannelECA(0.),fADCpedestalECA(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(fGeom)      delete fGeom;
//   if(fCalibData) delete fCalibData;
//   if(fCaloPed)   delete fCaloPed;

//   if (fDigitsArr) {
//     fDigitsArr->Clear("C");
//     delete fDigitsArr;
//   }
//   if (fRecPoints) {
//     fRecPoints->Delete();
//     delete fRecPoints;
//  }
  
  if(fClusterUnfolding) delete fClusterUnfolding;
}

//____________________________________________________________________________
Float_t  AliEMCALClusterizer::Calibrate(const Float_t amp, const Float_t time, const Int_t absId) 
{
  
  // Convert digitized amplitude into energy.
  // Calibration parameters are taken from calibration data base for raw data,
  // or from digitizer parameters for simulated data.
  if(fCalibData){
    
    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.
      return 0;
    }
    
    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
    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));
		  return 0;
    }
    //Check if time is too large or too small, indication of a noisy channel, remove in this case
    if(time > fTimeMax || time < fTimeMin) return 0;
    
    if (AliEMCALClusterizer::fgkIsInputCalibrated == kTRUE)
    {
      AliDebug(10, Form("Input already calibrated!"));
      return amp;
    }
	  
    fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
    fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
    
    
    return -fADCpedestalECA + amp * fADCchannelECA ;        
    
  }
  else{ //Return energy with default parameters if calibration is not available
    
    if (AliEMCALClusterizer::fgkIsInputCalibrated == kTRUE)
    {
      AliDebug(10, Form("Input already calibrated!"));
      return amp;
    }    
    
    return -fADCpedestalECA + amp * fADCchannelECA ; 
  }
  
}

//____________________________________________________________________________
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");
  
  //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");
  
  //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
  fNumberOfECAClusters = 0 ;
  fCalibData           = 0 ;
  fCaloPed             = 0 ;
	
  const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
  if(!recParam) {
    AliFatal("Reconstruction parameters for EMCAL not set!");
  } else {
    fECAClusteringThreshold = recParam->GetClusteringThreshold();
    fECAW0                  = recParam->GetW0();
    fMinECut                = recParam->GetMinECut();    
    fToUnfold               = recParam->GetUnfold();
    fECALocMaxCut           = recParam->GetLocMaxCut();
    fTimeCut                = recParam->GetTimeCut();
    fTimeMin                = recParam->GetTimeMin();
    fTimeMax                = recParam->GetTimeMax();
    
    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
  }// recparam not null
  
}


//____________________________________________________________________________
void AliEMCALClusterizer::Print(Option_t * /*option*/)const
{
  // Print clusterizer parameters
  
  TString message("\n") ; 
  
  if( strcmp(GetName(), "") !=0 ){
    
    // 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("------------------------------------------------------------------"); 
  }
  else
    printf("AliEMCALClusterizer not initialized ") ;
}

//____________________________________________________________________________
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){
        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
  fTreeR = clustersTree;
  
  AliDebug(9, "Making array for EMCAL clusters");
  fRecPoints = new TObjArray(100) ;
  Int_t split = 0;
  Int_t bufsize = 32000;
  fTreeR->Branch("EMCALECARP", "TObjArray", &fRecPoints, bufsize, split);
}

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