[u/mrichter/AliRoot.git] / JETAN / AliJetAODFillUnitArrayEMCalDigits.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 AOD reader for jet analysis
// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
//======================================================================

class TSystem;
class TLorentzVector;
class TGeoManager;
class TArrayS;

// --- AliRoot header files ---
#include "AliJetUnitArray.h"
#include "AliJetAODFillUnitArrayEMCalDigits.h"
class AliJetFinder;

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


ClassImp(AliJetAODFillUnitArrayEMCalDigits)

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

//_____________________________________________________________________________
AliJetAODFillUnitArrayEMCalDigits::AliJetAODFillUnitArrayEMCalDigits(AliAODEvent *aod)
  : AliJetFillUnitArray(),
    fAOD(aod),
    fNIn(0),
    fCluster(0),
    fNCEMCAL(0),
    fNCPHOS(0),
    fNCCalo(0),
    fApplyElectronCorrection(kFALSE),
    fApplyFractionHadronicCorrection(kFALSE),
    fFractionHadronicCorrection(0.3),
    fClus(0x0),
    fNDigitEmcal(0),
    fNDigitEmcalCut(0)
{
  // constructor
}

//_____________________________________________________________________________
AliJetAODFillUnitArrayEMCalDigits::AliJetAODFillUnitArrayEMCalDigits(const AliJetAODFillUnitArrayEMCalDigits &det)
  : AliJetFillUnitArray(det),
    fAOD(det.fAOD),
    fNIn(det.fNIn),
    fCluster(det.fCluster),
    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
}

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

  fAOD = other.fAOD;
  fNIn = other.fNIn;
  fCluster = other.fCluster;
  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);

}

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


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

  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
    //------------------------------------------------------------------
    AliAODCaloCells &cells= *(fAOD->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);
      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;
      }

      // Hadronic Correction
      double correction=0;
// Temporarily commented - will be updated removing the AliAODpid dependence
/*
      if (fApplyFractionHadronicCorrection) {
        TArrayS* matched = new TArrayS();
	GetTracksPointingToCell(matched, etaD, phiD,0.02);

	double ptot = 0;
	for(int itrack = 0; itrack <  matched->GetSize(); itrack++)
	  {
	    const short indexS = matched->At(itrack);
	    if (indexS>-1)
	      {	
		AliAODTrack *track = fAOD->GetTrack(indexS);
		ptot += track->P();
	      }
	  }
	
	correction = ptot * fFractionHadronicCorrection;

	if (ptot>0 && fDebug>1 ) {
	   printf("SUM of track momentum for this cell: p=%f \n", ptot);	
	   printf("fractional correction=%f \n", fFractionHadronicCorrection);
	   printf("TOTAL correction=%f \n", correction );
	}
	
	delete matched;
	
      }//end hadronic correction
*/
      double enerCorr = digitEn;
      if (correction >= enerCorr) enerCorr = 0;
      else enerCorr -= correction;
      if (correction>0 && fDebug>1) {
        printf("AliJetAODFillUnitArrayEMCALDigits---Hadronic Correction: uncorrected E=%f, corrected E=%f \n", digitEn, enerCorr);
      }

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

      // 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 if !fCluster
  else { // Keep digit information from clusterization

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

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

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

    // Get CaloCells
    AliAODCaloCells &cells= *(fAOD->GetEMCALCells());

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

      // Get the cluster info

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

      Int_t     digMult = fClus->GetNCells() ;
      UShort_t *digID   = fClus->GetCellsAbsId() ;
      Int_t     trackIndex = -1;
      if(fClus->GetNTracksMatched()!=0) 
	trackIndex = ((AliAODTrack*)fClus->GetTrackMatched(0))->GetID();

      // 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) ;
        Float_t  digitAmp     = cells.GetCellAmplitude(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);

	// Hadronic Correction
	double correction=0;
// Temporarily commented - will be updated removing the AliAODpid dependence
/*
	if (fApplyFractionHadronicCorrection) {
	  TArrayS* matched = new TArrayS();
	  GetTracksPointingToCell(matched, etaD, phiD,0.02);
	  
	  double ptot = 0;
	  for(int itrack = 0; itrack <  matched->GetSize(); itrack++)
	    {
	      const short indexS = matched->At(itrack);
	      if (indexS>-1)
		{	
		  AliAODTrack *track = fAOD->GetTrack(indexS);
		  ptot += track->P();
		}
	    }
	
	  correction = ptot * fFractionHadronicCorrection;

	  if (ptot>0 && fDebug>1 ) {
	    printf("SUM of track momentum for this cell: p=%f \n", ptot);	
	    printf("fractional correction=%f \n", fFractionHadronicCorrection);
	    printf("TOTAL correction=%f \n", correction );
	  }
	
	  delete matched;
	
	}//end hadronic correction
*/	
	double enerCorr = digitEn;
	if (correction >= enerCorr) enerCorr = 0;
	else enerCorr -= correction;
	if (correction>0 && fDebug>1) {
	  printf("AliJetAODFillUnitArrayEMCALDigits---Hadronic Correction: uncorrected E=%f, corrected E=%f \n", digitEn, enerCorr);
	}

	Float_t digitEt = enerCorr*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);

	// Signal or background
	// 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 AliJetAODFillUnitArrayEMCalDigits::GetTracksPointingToCell(TArrayS* array,Double_t eta, Double_t phi, Double_t cut)
//{ // Temporarily commented -> will be corrected removing the dependence to AliAODPid
//// Get all tracks pointing to cell 
//  int size=0;
//  
//  for (Int_t itrk =  0; itrk <  fAOD->GetNumberOfTracks() ; itrk++) { //track loop
//  
//    AliAODTrack * track = (AliAODTrack*) fAOD->GetTrack(itrk) ;  
//    AliAODPid*    pid   = (AliAODPid*) track->GetDetPid();
//    
//    if(pid) {
//      Double_t emcpos[3];
//      pid->GetEMCALPosition(emcpos);      
//      TVector3 tpos(emcpos[0],emcpos[1],emcpos[2]);
//
//      Double_t deta = tpos.Eta() - eta;
//      Double_t dphi = tpos.Phi() - phi;
//
//      if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
//      if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
//
//      Double_t res = sqrt(dphi*dphi + deta*deta);
//      
//      if (res< cut) {
//        //add this track-index
//        size++;
//        array->Set( size );
//        array->AddAt( itrk, (size-1) ); 
//	if(fDebug>1) printf("MTH:: track %d matched cell at eta=%f , phi=%f \n", itrk, eta, phi);
//		      
//      }
//    }
//  }
//
//}

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