[u/mrichter/AliRoot.git] / JETAN / AliJetESDFillUnitArrayEMCalDigits.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.                  *
 **************************************************************************/
 
//---------------------------------------------------------------------
// Fill Unit Array 
// Called by ESD reader for jet analysis
// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
//---------------------------------------------------------------------

// --- ROOT system ---
#include <TLorentzVector.h>
#include <TMath.h>

// --- AliRoot header files ---
#include "AliESDCaloCluster.h"
#include "AliESDCaloCells.h"
#include "AliJetUnitArray.h"
#include "AliJetESDFillUnitArrayEMCalDigits.h"
// #include "AliEMCALCalibData.h"
// #include "AliCDBManager.h"
// // class AliCDBStorage;
// #include "AliCDBEntry.h"

// --- ROOT system ---
class TSystem;
class TGeoManager;
class TArrayS;

// --- AliRoot header files ---
class AliJetFinder;
class AliJetReader;
class AliJetESDReader;
// class AliEMCALCalibData;
// class AliCDBManager;
// class AliCDBStorage;
// class AliCDBEntry;

using std::cout;
using std::endl;
ClassImp(AliJetESDFillUnitArrayEMCalDigits)

//_____________________________________________________________________________
AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits()
  : AliJetFillUnitArray(),
    fESD(0),
    fNIn(0),
    fOpt(0),
    fCluster(0),
    fDebug(0),
    fNCEMCAL(0),
    fNCPHOS(0),
    fNCCalo(0),
    fApplyElectronCorrection(kFALSE),
    fApplyFractionHadronicCorrection(kFALSE),
    fFractionHadronicCorrection(0.3),
    fClus(0x0),
    fNDigitEmcal(0),
    fNDigitEmcalCut(0)
{
  // constructor
}

//_____________________________________________________________________________
AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits(AliESDEvent *esd)
  : AliJetFillUnitArray(),
    fESD(esd),
    fNIn(0),
    fOpt(0),
    fCluster(0),
    fDebug(0),
    fNCEMCAL(0),
    fNCPHOS(0),
    fNCCalo(0),
    fApplyElectronCorrection(kFALSE),
    fApplyFractionHadronicCorrection(kFALSE),
    fFractionHadronicCorrection(0.3),
    fClus(0x0),
    fNDigitEmcal(0),
    fNDigitEmcalCut(0)
{
  // constructor
}

//_____________________________________________________________________________
AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits(const AliJetESDFillUnitArrayEMCalDigits &det)
  : AliJetFillUnitArray(det),
    fESD(det.fESD),
    fNIn(det.fNIn),
    fOpt(det.fOpt),
    fCluster(det.fCluster),
    fDebug(det.fDebug),
    fNCEMCAL(det.fNCEMCAL),
    fNCPHOS(det.fNCPHOS),
    fNCCalo(det.fNCCalo),
    fApplyElectronCorrection(det.fApplyElectronCorrection),
    fApplyFractionHadronicCorrection(det.fApplyFractionHadronicCorrection),
    fFractionHadronicCorrection(det.fFractionHadronicCorrection),
    fClus(det.fClus),
    fNDigitEmcal(det.fNDigitEmcal),
    fNDigitEmcalCut(det.fNDigitEmcalCut)
{
  // Copy constructor
}

//_____________________________________________________________________________
AliJetESDFillUnitArrayEMCalDigits& AliJetESDFillUnitArrayEMCalDigits::operator=(const AliJetESDFillUnitArrayEMCalDigits& other)
{
  // Assignment

  fESD = other.fESD;
  fNIn = other.fNIn;
  fOpt = other.fOpt;
  fCluster = other.fCluster;
  fDebug = other.fDebug;
  fNCEMCAL = other.fNCEMCAL;
  fNCPHOS = other.fNCPHOS;
  fNCCalo = other.fNCCalo;
  fApplyElectronCorrection = other.fApplyElectronCorrection;
  fApplyFractionHadronicCorrection = other.fApplyFractionHadronicCorrection;
  fFractionHadronicCorrection = other.fFractionHadronicCorrection;
  fClus = other.fClus;
  fNDigitEmcal = other.fNDigitEmcal;
  fNDigitEmcalCut = other.fNDigitEmcalCut;

  return (*this);

}

//_____________________________________________________________________________
AliJetESDFillUnitArrayEMCalDigits::~AliJetESDFillUnitArrayEMCalDigits()
{
  // destructor
}

//_____________________________________________________________________________
void AliJetESDFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/)
{
  //
  // Main method.
  // Fill the unit array with the neutral particle information from the EMCal cells in ESD
  //

  fDebug = fReaderHeader->GetDebug();
  fOpt = fReaderHeader->GetDetector();
  fCluster = fReaderHeader->GetCluster();

  // Init parameters
  InitParameters();

//(not used ?)  Int_t   nDigitTot      = 0;
  Int_t   goodDigit      = 0;
//(not used ?)  Int_t   beg            = 0;
//(not used ?)  Int_t   end            = 0; 
  Int_t   index          = 0;
//(not used ?)  Int_t   count          = 0;

//(not used ?)  Int_t nRefEnt = fRefArray->GetEntries();

  if(!fCluster) { // Keep all digit information
    // Loop over all cell information
    //------------------------------------------------------------------
    AliESDCaloCells &cells= *(fESD->GetEMCALCells());
    Int_t ncell = cells.GetNumberOfCells() ;
//(not used ?)    Int_t type = cells.GetType();

    for (Int_t icell=  0; icell <  ncell; icell++) {
      Int_t   digitID   = cells.GetCellNumber(icell);
      Float_t digitAmp  = cells.GetAmplitude(icell);
//(not used ?)      Float_t digitTime = cells.GetTime(icell);
      Float_t digitEn   = digitAmp*0.0153; // Last correct
      //  Float_t digitEn = Calibrate(digitAmp,digitID);

      Float_t etaD=-10., phiD=-10.;
      fGeom->EtaPhiFromIndex(digitID,etaD,phiD); 
      //  fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);

      phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
      
      Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));

      AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);

      if(uArray->GetUnitEnergy() == 0.) goodDigit++;
      uArray->SetUnitTrackID(digitID);

      Float_t unitEnergy = 0.;
      Bool_t ok = kFALSE;
      unitEnergy = uArray->GetUnitEnergy();

      if(unitEnergy==0){
	if(!fProcId){
	  new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
	  fProcId = kTRUE;
	}
	fRefArray->Add(uArray);
	fNDigitEmcal++;
	ok = kTRUE;
      }

      // Detector flag
      if(unitEnergy>0.)
	uArray->SetUnitDetectorFlag(kAll);
      else uArray->SetUnitDetectorFlag(kEmcal);
      
      uArray->SetUnitEnergy(unitEnergy+digitEt);

      uArray->SetUnitCutFlag(kPtHigher);

      // To be modified !!!
      uArray->SetUnitSignalFlag(kGood);

      // This is for jet multiplicity
      uArray->SetUnitClusterID(index);
	
      if(fDebug > 1) printf("goodDigit : %d\n", goodDigit);

    } // End loop over cells
  } // end if !fCluster
  else { // Keep digit information from clusterization

    // Loop over calo clusters
    //------------------------------------------------------------------

    //select EMCAL clusters only
    TRefArray * caloClusters  = new TRefArray();
    fESD->GetEMCALClusters(caloClusters);

    // Total number of EMCAL cluster
    Int_t nclus = caloClusters->GetEntries() ;
    Int_t beg   = 0;
    Float_t pos[3] ;

    // Get reconstructed vertex position
    Double_t vertexPosition[3] ;
    fESD->GetVertex()->GetXYZ(vertexPosition) ;

    // Get CaloCells
    AliESDCaloCells &cells= *(fESD->GetEMCALCells());
//(not used ?)    Int_t ncell = cells.GetNumberOfCells() ;
//(not used ?)    Int_t type = cells.GetType();

    for(Int_t j = beg; j < nclus; j++) { // loop over clusters
      // Retrieve cluster from esd
      fClus = (AliESDCaloCluster *) caloClusters->At(j) ;

      // Get the cluster info
//(not used ?)      Float_t energy         = fClus->E() ;
//(not used ?)      Int_t   iprim          = fClus->GetLabel();

      fClus->GetPosition(pos) ;
      TVector3 vpos(pos[0],pos[1],pos[2]) ;
      TLorentzVector p ;
      fClus->GetMomentum(p,vertexPosition);

      Int_t     digMult = fClus->GetNCells() ;
      UShort_t *digID   = fClus->GetCellsAbsId() ;
//(not used ?)      Double_t *digAmpFrac = fClus->GetCellsAmplitudeFraction() ;
      Int_t     trackIndex = fClus->GetTrackMatchedIndex();

      // Do double-counted electron correction 
      if (fApplyElectronCorrection != 0 && trackIndex !=-1 )
	{
          // The electron correction go there
          // Under construction !!!!
	}  // End of Electron correction

      // Get CaloCells of cluster and fill the unitArray
      for(Int_t i = 0; i < digMult ; i++){
        Int_t    digitID   = digID[i]; // or clus->GetCellNumber(i) ;
//(not used ?)        Double_t digitAmpFrac = digAmpFrac[i];
        Float_t  digitAmp     = cells.GetCellAmplitude(digitID) ;
//(not used ?)        Float_t  digitTime    = cells.GetCellTime(digitID);

        // Calibration for an energy in GeV
        Float_t digitEn = digitAmp*0.0153;

        Float_t etaD=-10., phiD=-10.;
        fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
        //  fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);

        phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);

        Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));

	cout << "Digit " << i << ", eta: " << etaD << ", phi: " << phiD << endl;

	AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
	if(uArray->GetUnitEnergy() == 0.) goodDigit++;
	uArray->SetUnitTrackID(digitID);

	Float_t unitEnergy = 0.;
	Bool_t ok = kFALSE;
	unitEnergy = uArray->GetUnitEnergy();

	if(unitEnergy==0){
	  if(!fProcId){
	    new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
	    fProcId = kTRUE;
	  }
	  fRefArray->Add(uArray);
	  fNDigitEmcal++;
	  ok = kTRUE;
	}

	// Detector flag
	if(unitEnergy>0.)
	  uArray->SetUnitDetectorFlag(kAll);
	else uArray->SetUnitDetectorFlag(kEmcal);
      
	uArray->SetUnitEnergy(unitEnergy+digitEt);

	uArray->SetUnitCutFlag(kPtHigher);

	// To be modified !!!
	uArray->SetUnitSignalFlag(kGood);

	// This is for jet multiplicity
	uArray->SetUnitClusterID(index);
	
	if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
      
      } // End loop over cells
    } // End loop over clusters
  } // end else

  fNIn += goodDigit;

  if(fDebug>1)
    {
      printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
      printf("goodDigit : %d\n", goodDigit);
    }
}

/*
//____________________________________________________________________________
void AliJetESDFillUnitArrayEMCalDigits::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 && (AliCDBManager::Instance()->IsDefaultStorageSet()))
    {
      AliCDBEntry *entry = (AliCDBEntry*)
        AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
      if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
    }

  if(!fCalibData)
    printf("************* Calibration parameters not found in CDB! ****************");
//    AliFatal("Calibration parameters not found in CDB!");


}

//____________________________________________________________________________
Float_t  AliJetESDFillUnitArrayEMCalDigits::Calibrate(Int_t amp, 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)
      printf("************* Did not get geometry from EMCALLoader ***************");
//      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();
      Error("Calibrate()"," Wrong cell id number : %i", AbsId);
      assert(0);
    }

    fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);

    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
    return -fADCpedestalECA + amp * fADCchannelECA ;

}
*/
Commit	Line	Data
be6e5811	1
	2	/**************************************************************************
	3	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	4	* *
	5	* Author: The ALICE Off-line Project. *
	6	* Contributors are mentioned in the code where appropriate. *
	7	* *
	8	* Permission to use, copy, modify and distribute this software and its *
	9	* documentation strictly for non-commercial purposes is hereby granted *
	10	* without fee, provided that the above copyright notice appears in all *
	11	* copies and that both the copyright notice and this permission notice *
	12	* appear in the supporting documentation. The authors make no claims *
	13	* about the suitability of this software for any purpose. It is *
	14	* provided "as is" without express or implied warranty. *
	15	**************************************************************************/
	16
	17	//---------------------------------------------------------------------
	18	// Fill Unit Array
	19	// Called by ESD reader for jet analysis
	20	// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
	21	//---------------------------------------------------------------------
	22
	23	// --- ROOT system ---
	24	#include <TLorentzVector.h>
	25	#include <TMath.h>
	26
	27	// --- AliRoot header files ---
	28	#include "AliESDCaloCluster.h"
	29	#include "AliESDCaloCells.h"
	30	#include "AliJetUnitArray.h"
	31	#include "AliJetESDFillUnitArrayEMCalDigits.h"
	32	// #include "AliEMCALCalibData.h"
	33	// #include "AliCDBManager.h"
	34	// // class AliCDBStorage;
	35	// #include "AliCDBEntry.h"
	36
	37	// --- ROOT system ---
	38	class TSystem;
	39	class TGeoManager;
	40	class TArrayS;
	41
	42	// --- AliRoot header files ---
	43	class AliJetFinder;
	44	class AliJetReader;
	45	class AliJetESDReader;
	46	// class AliEMCALCalibData;
	47	// class AliCDBManager;
	48	// class AliCDBStorage;
	49	// class AliCDBEntry;
	50
3a7af7bd	51	using std::cout;
3a7af7bd	52	using std::endl;
be6e5811	53	ClassImp(AliJetESDFillUnitArrayEMCalDigits)
	54
	55	//_____________________________________________________________________________
	56	AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits()
	57	: AliJetFillUnitArray(),
	58	fESD(0),
	59	fNIn(0),
	60	fOpt(0),
	61	fCluster(0),
	62	fDebug(0),
	63	fNCEMCAL(0),
	64	fNCPHOS(0),
	65	fNCCalo(0),
	66	fApplyElectronCorrection(kFALSE),
	67	fApplyFractionHadronicCorrection(kFALSE),
	68	fFractionHadronicCorrection(0.3),
	69	fClus(0x0),
	70	fNDigitEmcal(0),
	71	fNDigitEmcalCut(0)
	72	{
	73	// constructor
	74	}
	75
	76	//_____________________________________________________________________________
	77	AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits(AliESDEvent *esd)
	78	: AliJetFillUnitArray(),
	79	fESD(esd),
	80	fNIn(0),
	81	fOpt(0),
	82	fCluster(0),
	83	fDebug(0),
	84	fNCEMCAL(0),
	85	fNCPHOS(0),
	86	fNCCalo(0),
	87	fApplyElectronCorrection(kFALSE),
	88	fApplyFractionHadronicCorrection(kFALSE),
	89	fFractionHadronicCorrection(0.3),
	90	fClus(0x0),
	91	fNDigitEmcal(0),
	92	fNDigitEmcalCut(0)
	93	{
	94	// constructor
	95	}
	96
	97	//_____________________________________________________________________________
	98	AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits(const AliJetESDFillUnitArrayEMCalDigits &det)
	99	: AliJetFillUnitArray(det),
	100	fESD(det.fESD),
	101	fNIn(det.fNIn),
	102	fOpt(det.fOpt),
	103	fCluster(det.fCluster),
	104	fDebug(det.fDebug),
	105	fNCEMCAL(det.fNCEMCAL),
	106	fNCPHOS(det.fNCPHOS),
	107	fNCCalo(det.fNCCalo),
	108	fApplyElectronCorrection(det.fApplyElectronCorrection),
	109	fApplyFractionHadronicCorrection(det.fApplyFractionHadronicCorrection),
	110	fFractionHadronicCorrection(det.fFractionHadronicCorrection),
	111	fClus(det.fClus),
	112	fNDigitEmcal(det.fNDigitEmcal),
	113	fNDigitEmcalCut(det.fNDigitEmcalCut)
	114	{
	115	// Copy constructor
	116	}
117
118	//_____________________________________________________________________________
119	AliJetESDFillUnitArrayEMCalDigits& AliJetESDFillUnitArrayEMCalDigits::operator=(const AliJetESDFillUnitArrayEMCalDigits& other)
120	{
121	// Assignment
122
123	fESD = other.fESD;
124	fNIn = other.fNIn;
125	fOpt = other.fOpt;
126	fCluster = other.fCluster;
127	fDebug = other.fDebug;
128	fNCEMCAL = other.fNCEMCAL;
129	fNCPHOS = other.fNCPHOS;
130	fNCCalo = other.fNCCalo;
131	fApplyElectronCorrection = other.fApplyElectronCorrection;
132	fApplyFractionHadronicCorrection = other.fApplyFractionHadronicCorrection;
133	fFractionHadronicCorrection = other.fFractionHadronicCorrection;
134	fClus = other.fClus;
135	fNDigitEmcal = other.fNDigitEmcal;
136	fNDigitEmcalCut = other.fNDigitEmcalCut;
137
138	return (*this);
139
140	}
141
142	//_____________________________________________________________________________
143	AliJetESDFillUnitArrayEMCalDigits::~AliJetESDFillUnitArrayEMCalDigits()
144	{
145	// destructor
146	}
147
148	//_____________________________________________________________________________
149	void AliJetESDFillUnitArrayEMCalDigits::Exec(Option_t* const /option/)
150	{
151	//
152	// Main method.
153	// Fill the unit array with the neutral particle information from the EMCal cells in ESD
154	//
155
156	fDebug = fReaderHeader->GetDebug();
157	fOpt = fReaderHeader->GetDetector();
158	fCluster = fReaderHeader->GetCluster();
159
160	// Init parameters
161	InitParameters();
162
163	//(not used ?) Int_t nDigitTot = 0;
164	Int_t goodDigit = 0;
165	//(not used ?) Int_t beg = 0;
166	//(not used ?) Int_t end = 0;
167	Int_t index = 0;
168	//(not used ?) Int_t count = 0;
169
170	//(not used ?) Int_t nRefEnt = fRefArray->GetEntries();
171
172	if(!fCluster) { // Keep all digit information
173	// Loop over all cell information
174	//------------------------------------------------------------------
175	AliESDCaloCells &cells= *(fESD->GetEMCALCells());
176	Int_t ncell = cells.GetNumberOfCells() ;
177	//(not used ?) Int_t type = cells.GetType();
178
179	for (Int_t icell= 0; icell < ncell; icell++) {
180	Int_t digitID = cells.GetCellNumber(icell);
181	Float_t digitAmp = cells.GetAmplitude(icell);
182	//(not used ?) Float_t digitTime = cells.GetTime(icell);
183	Float_t digitEn = digitAmp*0.0153; // Last correct
184	// Float_t digitEn = Calibrate(digitAmp,digitID);
185
186	Float_t etaD=-10., phiD=-10.;
187	fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
188	// fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
189
190	phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
191
192	Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
193
194	AliJetUnitArray uArray = (AliJetUnitArray)fUnitArray->At(digitID);
195
196	if(uArray->GetUnitEnergy() == 0.) goodDigit++;
197	uArray->SetUnitTrackID(digitID);
198
199	Float_t unitEnergy = 0.;
200	Bool_t ok = kFALSE;
201	unitEnergy = uArray->GetUnitEnergy();
202
203	if(unitEnergy==0){
204	if(!fProcId){
205	new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
206	fProcId = kTRUE;
207	}
208	fRefArray->Add(uArray);
209	fNDigitEmcal++;
210	ok = kTRUE;
211	}
212
213	// Detector flag
214	if(unitEnergy>0.)
215	uArray->SetUnitDetectorFlag(kAll);
216	else uArray->SetUnitDetectorFlag(kEmcal);
217
218	uArray->SetUnitEnergy(unitEnergy+digitEt);
219
220	uArray->SetUnitCutFlag(kPtHigher);
221
222	// To be modified !!!
223	uArray->SetUnitSignalFlag(kGood);
224
225	// This is for jet multiplicity
226	uArray->SetUnitClusterID(index);
227
228	if(fDebug > 1) printf("goodDigit : %d\n", goodDigit);
229
230	} // End loop over cells
231	} // end if !fCluster
232	else { // Keep digit information from clusterization
233
234	// Loop over calo clusters
235	//------------------------------------------------------------------
236
237	//select EMCAL clusters only
238	TRefArray * caloClusters = new TRefArray();
239	fESD->GetEMCALClusters(caloClusters);
240
241	// Total number of EMCAL cluster
242	Int_t nclus = caloClusters->GetEntries() ;
243	Int_t beg = 0;
244	Float_t pos[3] ;
245
246	// Get reconstructed vertex position
247	Double_t vertexPosition[3] ;
248	fESD->GetVertex()->GetXYZ(vertexPosition) ;
249
250	// Get CaloCells
251	AliESDCaloCells &cells= *(fESD->GetEMCALCells());
252	//(not used ?) Int_t ncell = cells.GetNumberOfCells() ;
253	//(not used ?) Int_t type = cells.GetType();
254
255	for(Int_t j = beg; j < nclus; j++) { // loop over clusters
256	// Retrieve cluster from esd
08c28025	257	fClus = (AliESDCaloCluster *) caloClusters->At(j) ;
be6e5811	258
	259	// Get the cluster info
	260	//(not used ?) Float_t energy = fClus->E() ;
	261	//(not used ?) Int_t iprim = fClus->GetLabel();
	262
	263	fClus->GetPosition(pos) ;
	264	TVector3 vpos(pos[0],pos[1],pos[2]) ;
	265	TLorentzVector p ;
	266	fClus->GetMomentum(p,vertexPosition);
	267
	268	Int_t digMult = fClus->GetNCells() ;
	269	UShort_t *digID = fClus->GetCellsAbsId() ;
	270	//(not used ?) Double_t *digAmpFrac = fClus->GetCellsAmplitudeFraction() ;
c8fe2783	271	Int_t trackIndex = fClus->GetTrackMatchedIndex();
be6e5811	272
	273	// Do double-counted electron correction
	274	if (fApplyElectronCorrection != 0 && trackIndex !=-1 )
	275	{
	276	// The electron correction go there
	277	// Under construction !!!!
	278	} // End of Electron correction
	279
	280	// Get CaloCells of cluster and fill the unitArray
	281	for(Int_t i = 0; i < digMult ; i++){
	282	Int_t digitID = digID[i]; // or clus->GetCellNumber(i) ;
	283	//(not used ?) Double_t digitAmpFrac = digAmpFrac[i];
	284	Float_t digitAmp = cells.GetCellAmplitude(digitID) ;
	285	//(not used ?) Float_t digitTime = cells.GetCellTime(digitID);
	286
	287	// Calibration for an energy in GeV
	288	Float_t digitEn = digitAmp*0.0153;
	289
	290	Float_t etaD=-10., phiD=-10.;
	291	fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
	292	// fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
	293
	294	phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
	295
	296	Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
	297
	298	cout << "Digit " << i << ", eta: " << etaD << ", phi: " << phiD << endl;
	299
	300	AliJetUnitArray uArray = (AliJetUnitArray)fUnitArray->At(digitID);
	301	if(uArray->GetUnitEnergy() == 0.) goodDigit++;
	302	uArray->SetUnitTrackID(digitID);
	303
	304	Float_t unitEnergy = 0.;
	305	Bool_t ok = kFALSE;
	306	unitEnergy = uArray->GetUnitEnergy();
	307
	308	if(unitEnergy==0){
	309	if(!fProcId){
	310	new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
	311	fProcId = kTRUE;
	312	}
	313	fRefArray->Add(uArray);
	314	fNDigitEmcal++;
	315	ok = kTRUE;
	316	}
	317
	318	// Detector flag
	319	if(unitEnergy>0.)
	320	uArray->SetUnitDetectorFlag(kAll);
	321	else uArray->SetUnitDetectorFlag(kEmcal);
	322
	323	uArray->SetUnitEnergy(unitEnergy+digitEt);
	324
	325	uArray->SetUnitCutFlag(kPtHigher);
	326
	327	// To be modified !!!
	328	uArray->SetUnitSignalFlag(kGood);
	329
	330	// This is for jet multiplicity
	331	uArray->SetUnitClusterID(index);
	332
	333	if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
	334
	335	} // End loop over cells
336	} // End loop over clusters
337	} // end else
338
339	fNIn += goodDigit;
340
341	if(fDebug>1)
342	{
343	printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
344	printf("goodDigit : %d\n", goodDigit);
345	}
346	}
347
348	/*
349	//____________________________________________________________________________
350	void AliJetESDFillUnitArrayEMCalDigits::GetCalibrationParameters()
351	{
352	// Set calibration parameters:
353	// if calibration database exists, they are read from database,
354	// otherwise, they are taken from digitizer.
355	//
356	// It is a user responsilibity to open CDB before reconstruction,
357	// for example:
358	// AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
359
360	//Check if calibration is stored in data base
361
362	if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
363	{
364	AliCDBEntry entry = (AliCDBEntry)
365	AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
366	if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
367	}
368
369	if(!fCalibData)
370	printf("*********** Calibration parameters not found in CDB! **************");
371	// AliFatal("Calibration parameters not found in CDB!");
372
373
374	}
375
376	//____________________________________________________________________________
377	Float_t AliJetESDFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
378	{
379
380	// Convert digitized amplitude into energy.
381	// Calibration parameters are taken from calibration data base for raw data,
382	// or from digitizer parameters for simulated data.
383
384	if(fCalibData){
385
386	if (fGeom==0)
387	printf("*********** Did not get geometry from EMCALLoader *************");
388	// AliFatal("Did not get geometry from EMCALLoader") ;
389
390	Int_t iSupMod = -1;
391	Int_t nModule = -1;
392	Int_t nIphi = -1;
393	Int_t nIeta = -1;
394	Int_t iphi = -1;
395	Int_t ieta = -1;
396
397	Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
398	if(!bCell) {
399	// fGeom->PrintGeometry();
400	Error("Calibrate()"," Wrong cell id number : %i", AbsId);
401	assert(0);
402	}
403
404	fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
405
406	fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
407	fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
408
409	return -fADCpedestalECA + amp * fADCchannelECA ;
410
411	}
412	else //Return energy with default parameters if calibration is not available
413	return -fADCpedestalECA + amp * fADCchannelECA ;
414
415	}
416	*/
417