[u/mrichter/AliRoot.git] / JETAN / AliJetFillUnitArrayEMCalDigits.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)

// --- Standard library ---
#include <Riostream.h>
#include <assert.h>

// --- ROOT system ---
#include <TSystem.h>
#include <TLorentzVector.h>
#include <TRefArray.h> 
#include "TTask.h"
#include <TGeoManager.h>
#include <TMath.h>
#include <TClonesArray.h>
#include <TArrayS.h>

// --- AliRoot header files ---
#include "AliJetFinder.h"
#include "AliJetReaderHeader.h"
#include "AliJetReader.h"
#include "AliJetESDReader.h"
#include "AliJetESDReaderHeader.h"
#include "AliESDEvent.h"
#include "AliESDVertex.h"
#include "AliJetDummyGeo.h"
#include "AliESDCaloCluster.h"
#include "AliESDCaloCells.h"
#include "AliJetUnitArray.h"
#include "AliJetFillUnitArrayEMCalDigits.h"
// Remove CDB dependence under construction
//#include "AliEMCALCalibData.h"
//#include "AliCDBManager.h"

//class AliCDBStorage;
//#include "AliCDBEntry.h"


ClassImp(AliJetFillUnitArrayEMCalDigits)

//_____________________________________________________________________________
AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits()
  : TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
    fESD(0),
    fNumUnits(0), 
    fEtaMinCal(0.),
    fEtaMaxCal(0.),
    fPhiMinCal(0.),
    fPhiMaxCal(0.),
    fNIn(0),
    fOpt(0),
    fCluster(0),
    fDebug(0),
    fNCEMCAL(0),
    fNCPHOS(0),
    fNCCalo(0),
    fTPCGrid(0x0),
    fEMCalGrid(0x0),
    fECorrection(0),
    fReaderHeader(0x0),
    fMomentumArray(0x0),
    fUnitArray(0x0),
    fRefArray(0x0),
    fProcId(kFALSE),
    fGeom(0x0),
    fClus(0x0),
    fNDigitEmcal(0),
    fNDigitEmcalCut(0),
    fCalibData(0x0),
    fADCchannelECA(0),
    fADCpedestalECA(0)
{
  // constructor
}

//_____________________________________________________________________________
AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(AliESDEvent *esd)
  : TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
    fESD(esd),
    fNumUnits(0), 
    fEtaMinCal(0.),
    fEtaMaxCal(0.),
    fPhiMinCal(0.),
    fPhiMaxCal(0.),
    fNIn(0),
    fOpt(0),
    fCluster(0),
    fDebug(0),
    fNCEMCAL(0),
    fNCPHOS(0),
    fNCCalo(0),
    fTPCGrid(0x0),
    fEMCalGrid(0x0),
    fECorrection(0),
    fReaderHeader(0x0),
    fMomentumArray(0x0),
    fUnitArray(0x0),
    fRefArray(0x0),
    fProcId(kFALSE),
    fGeom(0x0),
    fClus(0x0),
    fNDigitEmcal(0),
    fNDigitEmcalCut(0),
    fCalibData(0x0),
    fADCchannelECA(0),
    fADCpedestalECA(0)
{
  // constructor 2
}

//_____________________________________________________________________________
AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(const AliJetFillUnitArrayEMCalDigits &det)
  : TTask(det),//"AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
    fESD(det.fESD),
    fNumUnits(det.fNumUnits), 
    fEtaMinCal(det.fEtaMinCal),
    fEtaMaxCal(det.fEtaMaxCal),
    fPhiMinCal(det.fPhiMinCal),
    fPhiMaxCal(det.fPhiMaxCal),
    fNIn(det.fNIn),
    fOpt(det.fOpt),
    fCluster(det.fCluster),
    fDebug(det.fDebug),
    fNCEMCAL(det.fNCEMCAL),
    fNCPHOS(det.fNCPHOS),
    fNCCalo(det.fNCCalo),
    fTPCGrid(det.fTPCGrid),
    fEMCalGrid(det.fEMCalGrid),
    fECorrection(det.fECorrection),
    fReaderHeader(det.fReaderHeader),
    fMomentumArray(det.fMomentumArray),
    fUnitArray(det.fUnitArray),
    fRefArray(det.fRefArray),
    fProcId(det.fProcId),
    fGeom(det.fGeom),
    fClus(det.fClus),
    fNDigitEmcal(det.fNDigitEmcal),
    fNDigitEmcalCut(det.fNDigitEmcalCut),
    fCalibData(det.fCalibData),
    fADCchannelECA(det.fADCchannelECA),
    fADCpedestalECA(det.fADCpedestalECA)
{
  // Copy constructor
}

//_____________________________________________________________________________
AliJetFillUnitArrayEMCalDigits& AliJetFillUnitArrayEMCalDigits::operator=(const AliJetFillUnitArrayEMCalDigits& other)
{
  // Assignment
  
  fESD = other.fESD;
  fNumUnits = other.fNumUnits; 
  fEtaMinCal = other.fEtaMinCal;
  fEtaMaxCal = other.fEtaMaxCal;
  fPhiMinCal = other.fPhiMinCal;
  fPhiMaxCal = other.fPhiMaxCal;
  fNIn = other.fNIn;
  fOpt = other.fOpt;
  fCluster = other.fCluster;
  fDebug = other.fDebug;
  fNCEMCAL = other.fNCEMCAL;
  fNCPHOS = other.fNCPHOS;
  fNCCalo = other.fNCCalo;
  fTPCGrid = other.fTPCGrid;
  fEMCalGrid = other.fEMCalGrid;
  fECorrection = other.fECorrection;
  fReaderHeader = other.fReaderHeader;
  fMomentumArray = other.fMomentumArray;
  fUnitArray = other.fUnitArray;
  fRefArray = other.fRefArray;
  fProcId = other.fProcId;
  fGeom = other.fGeom;
  fClus = other.fClus;
  fNDigitEmcal = other.fNDigitEmcal;
  fNDigitEmcalCut = other.fNDigitEmcalCut;
  fCalibData = other.fCalibData;
  fADCchannelECA = other.fADCchannelECA;
  fADCpedestalECA = other.fADCpedestalECA;

  return (*this);
  
}

//____________________________________________________________________________
void AliJetFillUnitArrayEMCalDigits::InitParameters()
{
  fNumUnits = fGeom->GetNCells();      // Number of towers in EMCAL

  fEtaMinCal = fGeom->GetArm1EtaMin();
  fEtaMaxCal = fGeom->GetArm1EtaMax();
  fPhiMinCal = fGeom->GetArm1PhiMin();
  fPhiMaxCal = fGeom->GetArm1PhiMax(); 
  fClus      = 0;

  // Get calibration parameters from file or digitizer default values.
  // Under construction
  //  GetCalibrationParameters() ;

  if(fDebug>1) printf("\n EMCAL parameters initiated ! \n");

}

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

//_____________________________________________________________________________
void AliJetFillUnitArrayEMCalDigits::Exec(Option_t* /*option*/)
{
  //
  // Main method.
  // 
  
  fDebug   = fReaderHeader->GetDebug();
  fOpt     = fReaderHeader->GetDetector();
  fCluster = fReaderHeader->GetCluster();
  
  // Init parameters
  InitParameters();
  
  Int_t   goodDigit      = 0;
  Int_t   index          = 0;
  
  if(!fCluster) { // Keep all digit information
    // Loop over all cell information
    //------------------------------------------------------------------
    AliESDCaloCells &cells = *(fESD->GetEMCALCells());
    Int_t            ncell = cells.GetNumberOfCells() ;
    
    for (Int_t icell=  0; icell <  ncell; icell++) {
      Int_t     digitID   = cells.GetCellNumber(icell);
      Double_t  digitAmp  = cells.GetAmplitude(icell);
      Float_t digitEn   = digitAmp*0.0153; // Last correct
      //      Float_t digitEn = Calibrate((Int_t)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 vertex_position[3] ;
    fESD->GetVertex()->GetXYZ(vertex_position) ;

    // Get CaloCells
    AliESDCaloCells &cells= *(fESD->GetEMCALCells());
    
    for(Int_t j = beg; j < nclus; j++) { // loop over clusters
      // Retrieve cluster from esd
      AliESDCaloCluster *fClus = (AliESDCaloCluster *) caloClusters->At(j) ;

      // Get the cluster info
      Float_t energy       = fClus->E() ;
      Int_t iprim          = fClus->GetLabel();
      Int_t trackIndex     = fClus->GetTrackMatched();

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

      Int_t     digMult = fClus->GetNCells() ;
      UShort_t *digID   = fClus->GetCellsAbsId() ;
      //Print cluster info
      if(fDebug>2) cout<<"Cluster "<< j <<"; digits mult "<<digMult<<"; type "<<(Int_t )fClus->GetClusterType()
		       <<"; Energy "<<energy<< "; transverse energy:" <<  energy*TMath::Abs(TMath::Sin(EtaToTheta(vpos.Eta())))
		       <<"; Phi "<<vpos.Phi()<<"; Eta "<<vpos.Eta() <<"; label "<<iprim<<endl;

      // Do double-counted electron correction 
      if (fECorrection != 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) ;
	  Double_t  digitAmp    = cells.GetCellAmplitude(digitID) ;
	  
	  // Calibration for an energy in GeV
	  Float_t digitEn = digitAmp*0.0153;
	  //	  Float_t digitEn = Calibrate((Int_t)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);
	  
	} // 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 AliJetFillUnitArrayEMCalDigits::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  AliJetFillUnitArrayEMCalDigits::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 ***************");
    
//     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 ;

// }


//_____________________________________________________________________________
Float_t  AliJetFillUnitArrayEMCalDigits::EtaToTheta(Float_t arg)
{
  //  return (180./TMath::Pi())*2.*atan(exp(-arg));
  return 2.*atan(exp(-arg));


}
Commit	Line	Data
ee7de0dd	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	**************************************************************************/
8838ab7a	16
8838ab7a	17	// Fill Unit Array
ee7de0dd	18	// Called by ESD reader for jet analysis
ee7de0dd	19	// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
ee7de0dd	20
	21	// --- Standard library ---
	22	#include <Riostream.h>
8838ab7a	23	#include <assert.h>
ee7de0dd	24
	25	// --- ROOT system ---
	26	#include <TSystem.h>
	27	#include <TLorentzVector.h>
	28	#include <TRefArray.h>
	29	#include "TTask.h"
	30	#include <TGeoManager.h>
	31	#include <TMath.h>
	32	#include <TClonesArray.h>
	33	#include <TArrayS.h>
	34
	35	// --- AliRoot header files ---
	36	#include "AliJetFinder.h"
	37	#include "AliJetReaderHeader.h"
	38	#include "AliJetReader.h"
	39	#include "AliJetESDReader.h"
	40	#include "AliJetESDReaderHeader.h"
ee7de0dd	41	#include "AliESDEvent.h"
8838ab7a	42	#include "AliESDVertex.h"
ee7de0dd	43	#include "AliJetDummyGeo.h"
ee7de0dd	44	#include "AliESDCaloCluster.h"
8838ab7a	45	#include "AliESDCaloCells.h"
ee7de0dd	46	#include "AliJetUnitArray.h"
ee7de0dd	47	#include "AliJetFillUnitArrayEMCalDigits.h"
8838ab7a	48	// Remove CDB dependence under construction
	49	//#include "AliEMCALCalibData.h"
	50	//#include "AliCDBManager.h"
	51
	52	//class AliCDBStorage;
	53	//#include "AliCDBEntry.h"
	54
ee7de0dd	55
	56	ClassImp(AliJetFillUnitArrayEMCalDigits)
	57
	58	//_____________________________________________________________________________
	59	AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits()
	60	: TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
8838ab7a	61	fESD(0),
ee7de0dd	62	fNumUnits(0),
	63	fEtaMinCal(0.),
	64	fEtaMaxCal(0.),
	65	fPhiMinCal(0.),
	66	fPhiMaxCal(0.),
	67	fNIn(0),
	68	fOpt(0),
8838ab7a	69	fCluster(0),
ee7de0dd	70	fDebug(0),
	71	fNCEMCAL(0),
	72	fNCPHOS(0),
	73	fNCCalo(0),
	74	fTPCGrid(0x0),
	75	fEMCalGrid(0x0),
8838ab7a	76	fECorrection(0),
ee7de0dd	77	fReaderHeader(0x0),
	78	fMomentumArray(0x0),
	79	fUnitArray(0x0),
	80	fRefArray(0x0),
8838ab7a	81	fProcId(kFALSE),
ee7de0dd	82	fGeom(0x0),
	83	fClus(0x0),
	84	fNDigitEmcal(0),
8838ab7a	85	fNDigitEmcalCut(0),
	86	fCalibData(0x0),
	87	fADCchannelECA(0),
	88	fADCpedestalECA(0)
ee7de0dd	89	{
	90	// constructor
	91	}
	92
	93	//_____________________________________________________________________________
8838ab7a	94	AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(AliESDEvent *esd)
ee7de0dd	95	: TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
8838ab7a	96	fESD(esd),
ee7de0dd	97	fNumUnits(0),
	98	fEtaMinCal(0.),
	99	fEtaMaxCal(0.),
	100	fPhiMinCal(0.),
	101	fPhiMaxCal(0.),
	102	fNIn(0),
	103	fOpt(0),
8838ab7a	104	fCluster(0),
ee7de0dd	105	fDebug(0),
	106	fNCEMCAL(0),
	107	fNCPHOS(0),
	108	fNCCalo(0),
	109	fTPCGrid(0x0),
	110	fEMCalGrid(0x0),
8838ab7a	111	fECorrection(0),
ee7de0dd	112	fReaderHeader(0x0),
	113	fMomentumArray(0x0),
	114	fUnitArray(0x0),
	115	fRefArray(0x0),
8838ab7a	116	fProcId(kFALSE),
ee7de0dd	117	fGeom(0x0),
	118	fClus(0x0),
	119	fNDigitEmcal(0),
8838ab7a	120	fNDigitEmcalCut(0),
	121	fCalibData(0x0),
	122	fADCchannelECA(0),
	123	fADCpedestalECA(0)
ee7de0dd	124	{
8838ab7a	125	// constructor 2
	126	}
	127
	128	//_____________________________________________________________________________
	129	AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(const AliJetFillUnitArrayEMCalDigits &det)
	130	: TTask(det),//"AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
	131	fESD(det.fESD),
	132	fNumUnits(det.fNumUnits),
	133	fEtaMinCal(det.fEtaMinCal),
	134	fEtaMaxCal(det.fEtaMaxCal),
	135	fPhiMinCal(det.fPhiMinCal),
	136	fPhiMaxCal(det.fPhiMaxCal),
	137	fNIn(det.fNIn),
	138	fOpt(det.fOpt),
	139	fCluster(det.fCluster),
	140	fDebug(det.fDebug),
	141	fNCEMCAL(det.fNCEMCAL),
	142	fNCPHOS(det.fNCPHOS),
	143	fNCCalo(det.fNCCalo),
	144	fTPCGrid(det.fTPCGrid),
	145	fEMCalGrid(det.fEMCalGrid),
	146	fECorrection(det.fECorrection),
	147	fReaderHeader(det.fReaderHeader),
	148	fMomentumArray(det.fMomentumArray),
	149	fUnitArray(det.fUnitArray),
	150	fRefArray(det.fRefArray),
	151	fProcId(det.fProcId),
	152	fGeom(det.fGeom),
	153	fClus(det.fClus),
	154	fNDigitEmcal(det.fNDigitEmcal),
	155	fNDigitEmcalCut(det.fNDigitEmcalCut),
	156	fCalibData(det.fCalibData),
	157	fADCchannelECA(det.fADCchannelECA),
	158	fADCpedestalECA(det.fADCpedestalECA)
	159	{
	160	// Copy constructor
ee7de0dd	161	}
ee7de0dd	162
8838ab7a	163	//_____________________________________________________________________________
	164	AliJetFillUnitArrayEMCalDigits& AliJetFillUnitArrayEMCalDigits::operator=(const AliJetFillUnitArrayEMCalDigits& other)
	165	{
	166	// Assignment
	167
	168	fESD = other.fESD;
	169	fNumUnits = other.fNumUnits;
	170	fEtaMinCal = other.fEtaMinCal;
	171	fEtaMaxCal = other.fEtaMaxCal;
	172	fPhiMinCal = other.fPhiMinCal;
	173	fPhiMaxCal = other.fPhiMaxCal;
	174	fNIn = other.fNIn;
	175	fOpt = other.fOpt;
	176	fCluster = other.fCluster;
	177	fDebug = other.fDebug;
	178	fNCEMCAL = other.fNCEMCAL;
	179	fNCPHOS = other.fNCPHOS;
	180	fNCCalo = other.fNCCalo;
	181	fTPCGrid = other.fTPCGrid;
	182	fEMCalGrid = other.fEMCalGrid;
	183	fECorrection = other.fECorrection;
	184	fReaderHeader = other.fReaderHeader;
	185	fMomentumArray = other.fMomentumArray;
	186	fUnitArray = other.fUnitArray;
	187	fRefArray = other.fRefArray;
	188	fProcId = other.fProcId;
	189	fGeom = other.fGeom;
	190	fClus = other.fClus;
	191	fNDigitEmcal = other.fNDigitEmcal;
	192	fNDigitEmcalCut = other.fNDigitEmcalCut;
	193	fCalibData = other.fCalibData;
	194	fADCchannelECA = other.fADCchannelECA;
	195	fADCpedestalECA = other.fADCpedestalECA;
	196
	197	return (*this);
	198
	199	}
ee7de0dd	200
	201	//____________________________________________________________________________
	202	void AliJetFillUnitArrayEMCalDigits::InitParameters()
	203	{
	204	fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL
	205
	206	fEtaMinCal = fGeom->GetArm1EtaMin();
	207	fEtaMaxCal = fGeom->GetArm1EtaMax();
	208	fPhiMinCal = fGeom->GetArm1PhiMin();
	209	fPhiMaxCal = fGeom->GetArm1PhiMax();
	210	fClus = 0;
	211
8838ab7a	212	// Get calibration parameters from file or digitizer default values.
	213	// Under construction
	214	// GetCalibrationParameters() ;
	215
ee7de0dd	216	if(fDebug>1) printf("\n EMCAL parameters initiated ! \n");
	217
	218	}
	219
	220	//_____________________________________________________________________________
	221	AliJetFillUnitArrayEMCalDigits::~AliJetFillUnitArrayEMCalDigits()
	222	{
	223	// destructor
	224	}
	225
	226	//_____________________________________________________________________________
	227	void AliJetFillUnitArrayEMCalDigits::Exec(Option_t* /option/)
	228	{
	229	//
	230	// Main method.
8838ab7a	231	//
	232
	233	fDebug = fReaderHeader->GetDebug();
	234	fOpt = fReaderHeader->GetDetector();
	235	fCluster = fReaderHeader->GetCluster();
	236
ee7de0dd	237	// Init parameters
ee7de0dd	238	InitParameters();
ee7de0dd	239
ee7de0dd	240	Int_t goodDigit = 0;
8838ab7a	241	Int_t index = 0;
	242
	243	if(!fCluster) { // Keep all digit information
	244	// Loop over all cell information
	245	//------------------------------------------------------------------
	246	AliESDCaloCells &cells = *(fESD->GetEMCALCells());
	247	Int_t ncell = cells.GetNumberOfCells() ;
	248
	249	for (Int_t icell= 0; icell < ncell; icell++) {
	250	Int_t digitID = cells.GetCellNumber(icell);
	251	Double_t digitAmp = cells.GetAmplitude(icell);
	252	Float_t digitEn = digitAmp*0.0153; // Last correct
	253	// Float_t digitEn = Calibrate((Int_t)digitAmp,digitID);
	254
	255	Float_t etaD=-10., phiD=-10.;
	256	fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
	257	// fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
ee7de0dd	258
8838ab7a	259	phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
	260
	261	Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
ee7de0dd	262
8838ab7a	263	AliJetUnitArray uArray = (AliJetUnitArray)fUnitArray->At(digitID);
ee7de0dd	264
8838ab7a	265	if(uArray->GetUnitEnergy() == 0.) goodDigit++;
	266	uArray->SetUnitTrackID(digitID);
	267
	268	Float_t unitEnergy = 0.;
	269	Bool_t ok = kFALSE;
	270	unitEnergy = uArray->GetUnitEnergy();
	271
	272	if(unitEnergy==0){
	273	if(!fProcId){
	274	new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
	275	fProcId = kTRUE;
	276	}
	277	fRefArray->Add(uArray);
	278	fNDigitEmcal++;
	279	ok = kTRUE;
	280	}
	281
	282	// Detector flag
	283	if(unitEnergy>0.)
	284	uArray->SetUnitDetectorFlag(kAll);
	285	else uArray->SetUnitDetectorFlag(kEmcal);
ee7de0dd	286
8838ab7a	287	uArray->SetUnitEnergy(unitEnergy+digitEt);
ee7de0dd	288
8838ab7a	289	uArray->SetUnitCutFlag(kPtHigher);
	290
	291	// To be modified !!!
	292	uArray->SetUnitSignalFlag(kGood);
	293
	294	// This is for jet multiplicity
	295	uArray->SetUnitClusterID(index);
	296
	297	if(fDebug > 1) printf("goodDigit : %d\n", goodDigit);
	298
	299	} // End loop over cells
	300
	301	} // end if !fCluster
	302	else { // Keep digit information from clusterization
	303
	304	// Loop over calo clusters
	305	//------------------------------------------------------------------
	306
	307	//select EMCAL clusters only
	308	TRefArray * caloClusters = new TRefArray();
	309	fESD->GetEMCALClusters(caloClusters);
	310
	311	// Total number of EMCAL cluster
	312	Int_t nclus = caloClusters->GetEntries() ;
	313	Int_t beg = 0;
	314	Float_t pos[3] ;
	315
	316	// Get reconstructed vertex position
	317	Double_t vertex_position[3] ;
	318	fESD->GetVertex()->GetXYZ(vertex_position) ;
	319
	320	// Get CaloCells
	321	AliESDCaloCells &cells= *(fESD->GetEMCALCells());
	322
	323	for(Int_t j = beg; j < nclus; j++) { // loop over clusters
	324	// Retrieve cluster from esd
	325	AliESDCaloCluster fClus = (AliESDCaloCluster ) caloClusters->At(j) ;
	326
	327	// Get the cluster info
	328	Float_t energy = fClus->E() ;
	329	Int_t iprim = fClus->GetLabel();
	330	Int_t trackIndex = fClus->GetTrackMatched();
	331
	332	fClus->GetPosition(pos) ;
	333	TVector3 vpos(pos[0],pos[1],pos[2]) ;
	334	TLorentzVector p ;
	335	fClus->GetMomentum(p,vertex_position);
	336
	337	Int_t digMult = fClus->GetNCells() ;
	338	UShort_t *digID = fClus->GetCellsAbsId() ;
	339	//Print cluster info
	340	if(fDebug>2) cout<<"Cluster "<< j <<"; digits mult "<<digMult<<"; type "<<(Int_t )fClus->GetClusterType()
	341	<<"; Energy "<<energy<< "; transverse energy:" << energy*TMath::Abs(TMath::Sin(EtaToTheta(vpos.Eta())))
	342	<<"; Phi "<<vpos.Phi()<<"; Eta "<<vpos.Eta() <<"; label "<<iprim<<endl;
	343
	344	// Do double-counted electron correction
	345	if (fECorrection != 0 && trackIndex !=-1 )
	346	{
	347	// The electron correction go there
	348	// Under construction !!!!
	349
	350	} // End of Electron correction
	351
	352	// Get CaloCells of cluster and fill the unitArray
353	for(Int_t i = 0; i < digMult ; i++)
354	{
355	Int_t digitID = digID[i]; // or clus->GetCellNumber(i) ;
356	Double_t digitAmp = cells.GetCellAmplitude(digitID) ;
ee7de0dd	357
8838ab7a	358	// Calibration for an energy in GeV
	359	Float_t digitEn = digitAmp*0.0153;
	360	// Float_t digitEn = Calibrate((Int_t)digitAmp,digitID);
ee7de0dd	361
8838ab7a	362	Float_t etaD=-10., phiD=-10.;
	363	fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
	364	// fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
ee7de0dd	365
8838ab7a	366	phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
	367
	368	Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
	369
	370	AliJetUnitArray uArray = (AliJetUnitArray)fUnitArray->At(digitID);
	371	if(uArray->GetUnitEnergy() == 0.) goodDigit++;
	372	uArray->SetUnitTrackID(digitID);
	373
	374	Float_t unitEnergy = 0.;
	375	Bool_t ok = kFALSE;
	376	unitEnergy = uArray->GetUnitEnergy();
	377
	378	if(unitEnergy==0){
	379	if(!fProcId){
	380	new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
	381	fProcId = kTRUE;
	382	}
	383	fRefArray->Add(uArray);
	384	fNDigitEmcal++;
	385	ok = kTRUE;
	386	}
	387
	388	// Detector flag
	389	if(unitEnergy>0.)
	390	uArray->SetUnitDetectorFlag(kAll);
	391	else uArray->SetUnitDetectorFlag(kEmcal);
	392
	393	uArray->SetUnitEnergy(unitEnergy+digitEt);
ee7de0dd	394	uArray->SetUnitCutFlag(kPtHigher);
8838ab7a	395
	396	// To be modified !!!
	397	uArray->SetUnitSignalFlag(kGood);
	398
	399	// This is for jet multiplicity
	400	uArray->SetUnitClusterID(index);
	401
	402	} // End loop over cells
	403	} // End loop over clusters
	404	} // end else
	405
ee7de0dd	406	fNIn += goodDigit;
ee7de0dd	407
8838ab7a	408	if(fDebug>1)
	409	{
	410	printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
	411	printf("goodDigit : %d\n", goodDigit);
	412	}
	413
ee7de0dd	414	}
ee7de0dd	415
8838ab7a	416	// //____________________________________________________________________________
	417	// void AliJetFillUnitArrayEMCalDigits::GetCalibrationParameters()
	418	// {
	419	// // Set calibration parameters:
	420	// // if calibration database exists, they are read from database,
	421	// // otherwise, they are taken from digitizer.
	422	// //
	423	// // It is a user responsilibity to open CDB before reconstruction,
	424	// // for example:
	425	// // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
	426
	427	// //Check if calibration is stored in data base
	428
	429	// if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
	430	// {
	431	// AliCDBEntry entry = (AliCDBEntry)
	432	// AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
	433	// if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
	434	// }
	435
	436	// if(!fCalibData)
	437	// printf("*********** Calibration parameters not found in CDB! **************");
	438	// // AliFatal("Calibration parameters not found in CDB!");
	439
	440
	441	// }
	442
	443	// //____________________________________________________________________________
	444	// Float_t AliJetFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
	445	// {
	446
	447	// // Convert digitized amplitude into energy.
	448	// // Calibration parameters are taken from calibration data base for raw data,
	449	// // or from digitizer parameters for simulated data.
	450
	451	// if(fCalibData){
	452
	453	// if (fGeom==0)
	454	// printf("*********** Did not get geometry from EMCALLoader *************");
	455
	456	// Int_t iSupMod = -1;
	457	// Int_t nModule = -1;
	458	// Int_t nIphi = -1;
	459	// Int_t nIeta = -1;
	460	// Int_t iphi = -1;
	461	// Int_t ieta = -1;
	462
	463	// Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
	464	// if(!bCell) {
	465	// // fGeom->PrintGeometry();
	466	// Error("Calibrate()"," Wrong cell id number : %i", AbsId);
	467	// assert(0);
	468	// }
	469
	470	// fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
	471
	472	// fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
	473	// fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
	474
	475	// return -fADCpedestalECA + amp * fADCchannelECA ;
	476
	477	// }
	478	// else //Return energy with default parameters if calibration is not available
	479	// return -fADCpedestalECA + amp * fADCchannelECA ;
480
481	// }
482
483
ee7de0dd	484	//_____________________________________________________________________________
	485	Float_t AliJetFillUnitArrayEMCalDigits::EtaToTheta(Float_t arg)
	486	{
	487	// return (180./TMath::Pi())2.atan(exp(-arg));
	488	return 2.*atan(exp(-arg));
	489
	490
	491	}