[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;
      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;
	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)
{
// 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
	188	double correction=0;
	189	if (fApplyFractionHadronicCorrection) {
	190
	191	TArrayS* matched = new TArrayS();
	192	GetTracksPointingToCell(matched, etaD, phiD,0.02);
	193
	194	double ptot = 0;
	195	for(int itrack = 0; itrack < matched->GetSize(); itrack++)
	196	{
08c28025	197	const short indexS = matched->At(itrack);
08c28025	198	if (indexS>-1)
e36a3f22	199	{
08c28025	200	AliAODTrack *track = fAOD->GetTrack(indexS);
e36a3f22	201	ptot += track->P();
e36a3f22	202	}
	203	}
	204
	205	correction = ptot * fFractionHadronicCorrection;
	206
	207	if (ptot>0 && fDebug>1 ) {
	208	printf("SUM of track momentum for this cell: p=%f \n", ptot);
	209	printf("fractional correction=%f \n", fFractionHadronicCorrection);
	210	printf("TOTAL correction=%f \n", correction );
	211	}
	212
	213	delete matched;
	214
	215	}//end hadronic correction
	216
	217	double enerCorr = digitEn;
	218	if (correction >= enerCorr) enerCorr = 0;
	219	else enerCorr -= correction;
	220	if (correction>0 && fDebug>1) {
	221	printf("AliJetAODFillUnitArrayEMCALDigits---Hadronic Correction: uncorrected E=%f, corrected E=%f \n", digitEn, enerCorr);
	222	}
	223
	224	Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
	225
be6e5811	226	// Detector flag
	227	if(unitEnergy>0.)
	228	uArray->SetUnitDetectorFlag(kAll);
	229	else uArray->SetUnitDetectorFlag(kEmcal);
	230
	231	uArray->SetUnitEnergy(unitEnergy+digitEt);
	232
	233	uArray->SetUnitCutFlag(kPtHigher);
	234
	235	// To be modified !!!
	236	uArray->SetUnitSignalFlag(kGood);
	237
	238	// This is for jet multiplicity
	239	uArray->SetUnitClusterID(index);
	240
	241	if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
	242
	243	} // End loop over cells
	244	} // end if !fCluster
	245	else { // Keep digit information from clusterization
	246
	247	// Loop over calo clusters
	248	//------------------------------------------------------------------
	249
	250	// select EMCAL clusters only
	251	TRefArray * caloClusters = new TRefArray();
	252	fAOD->GetEMCALClusters(caloClusters);
	253
	254	// Total number of EMCAL cluster
	255	Int_t nclus = caloClusters->GetEntries() ;
	256	Int_t beg = 0;
	257	Float_t pos[3] ;
	258
	259	// Get CaloCells
	260	AliAODCaloCells &cells= *(fAOD->GetEMCALCells());
be6e5811	261
	262	for(Int_t j = beg; j < nclus; j++) { // loop over clusters
	263	// Retrieve cluster from aod
08c28025	264	fClus = (AliAODCaloCluster *) caloClusters->At(j) ;
be6e5811	265
be6e5811	266	// Get the cluster info
be6e5811	267
	268	fClus->GetPosition(pos) ;
	269	TVector3 vpos(pos[0],pos[1],pos[2]) ;
	270
	271	Int_t digMult = fClus->GetNCells() ;
	272	UShort_t *digID = fClus->GetCellsAbsId() ;
be6e5811	273	Int_t trackIndex = -1;
	274	if(fClus->GetNTracksMatched()!=0)
	275	trackIndex = ((AliAODTrack*)fClus->GetTrackMatched(0))->GetID();
e36a3f22	276
be6e5811	277	// Do double-counted electron correction
	278	if (fApplyElectronCorrection != 0 && trackIndex !=-1 )
	279	{
	280	// The electron correction go there
	281	// Under construction !!!!
	282	} // End of Electron correction
	283
	284	// Get CaloCells of cluster and fill the unitArray
	285	for(Int_t i = 0; i < digMult ; i++) {
	286	Int_t digitID = digID[i]; // or clus->GetCellNumber(i) ;
be6e5811	287	Float_t digitAmp = cells.GetCellAmplitude(digitID) ;
	288
	289	// Calibration for an energy in GeV
	290	Float_t digitEn = digitAmp*0.0153;
	291
	292	Float_t etaD=-10., phiD=-10.;
	293	fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
	294	// fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
	295
	296	phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
	297
e36a3f22	298	// Hadronic Correction
	299	double correction=0;
	300	if (fApplyFractionHadronicCorrection) {
	301
	302	TArrayS* matched = new TArrayS();
	303	GetTracksPointingToCell(matched, etaD, phiD,0.02);
	304
	305	double ptot = 0;
	306	for(int itrack = 0; itrack < matched->GetSize(); itrack++)
	307	{
08c28025	308	const short indexS = matched->At(itrack);
08c28025	309	if (indexS>-1)
e36a3f22	310	{
08c28025	311	AliAODTrack *track = fAOD->GetTrack(indexS);
e36a3f22	312	ptot += track->P();
e36a3f22	313	}
	314	}
	315
	316	correction = ptot * fFractionHadronicCorrection;
	317
	318	if (ptot>0 && fDebug>1 ) {
	319	printf("SUM of track momentum for this cell: p=%f \n", ptot);
	320	printf("fractional correction=%f \n", fFractionHadronicCorrection);
	321	printf("TOTAL correction=%f \n", correction );
	322	}
	323
	324	delete matched;
	325
	326	}//end hadronic correction
	327
	328	double enerCorr = digitEn;
	329	if (correction >= enerCorr) enerCorr = 0;
	330	else enerCorr -= correction;
	331	if (correction>0 && fDebug>1) {
	332	printf("AliJetAODFillUnitArrayEMCALDigits---Hadronic Correction: uncorrected E=%f, corrected E=%f \n", digitEn, enerCorr);
	333	}
	334
	335	Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
be6e5811	336
	337	AliJetUnitArray uArray = (AliJetUnitArray)fUnitArray->At(digitID);
	338	if(uArray->GetUnitEnergy() == 0.) goodDigit++;
	339	uArray->SetUnitTrackID(digitID);
	340
	341	Float_t unitEnergy = 0.;
	342	Bool_t ok = kFALSE;
	343	unitEnergy = uArray->GetUnitEnergy();
	344
	345	if(unitEnergy==0){
	346	if(!fProcId){
	347	new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
	348	fProcId = kTRUE;
	349	}
	350	fRefArray->Add(uArray);
	351	fNDigitEmcal++;
	352	ok = kTRUE;
	353	}
	354
	355	// Detector flag
	356	if(unitEnergy>0.)
	357	uArray->SetUnitDetectorFlag(kAll);
	358	else uArray->SetUnitDetectorFlag(kEmcal);
	359
	360	uArray->SetUnitEnergy(unitEnergy+digitEt);
	361
	362	uArray->SetUnitCutFlag(kPtHigher);
	363
	364	// Signal or background
	365	// To be modified !!!
	366	uArray->SetUnitSignalFlag(kGood);
	367
	368	// This is for jet multiplicity
	369	uArray->SetUnitClusterID(index);
	370
	371	if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
	372
	373	}// End loop over cells
	374	} // End loop over clusters
	375	} // end else
	376
	377	fNIn += goodDigit;
	378
	379	if(fDebug>1)
	380	{
	381	printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
	382	printf("goodDigit : %d\n", goodDigit);
	383	}
	384	}
	385
e36a3f22	386	//____________________________________________________________________________
	387	void AliJetAODFillUnitArrayEMCalDigits::GetTracksPointingToCell(TArrayS* array,Double_t eta, Double_t phi, Double_t cut)
	388	{
1240edf5	389	// Get all tracks pointing to cell
e36a3f22	390	int size=0;
	391
	392	for (Int_t itrk = 0; itrk < fAOD->GetNumberOfTracks() ; itrk++) { //track loop
	393
	394	AliAODTrack * track = (AliAODTrack*) fAOD->GetTrack(itrk) ;
	395	AliAODPid* pid = (AliAODPid*) track->GetDetPid();
	396
	397	if(pid) {
	398	Double_t emcpos[3];
	399	pid->GetEMCALPosition(emcpos);
	400	TVector3 tpos(emcpos[0],emcpos[1],emcpos[2]);
	401
	402	Double_t deta = tpos.Eta() - eta;
	403	Double_t dphi = tpos.Phi() - phi;
	404
	405	if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
	406	if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
	407
	408	Double_t res = sqrt(dphidphi + detadeta);
	409
	410	if (res< cut) {
	411	//add this track-index
	412	size++;
	413	array->Set( size );
	414	array->AddAt( itrk, (size-1) );
	415	if(fDebug>1) printf("MTH:: track %d matched cell at eta=%f , phi=%f \n", itrk, eta, phi);
	416
	417	}
	418	}
	419	}
	420
	421	}
	422
be6e5811	423	/*
	424	//____________________________________________________________________________
	425	void AliJetAODFillUnitArrayEMCalDigits::GetCalibrationParameters()
	426	{
	427	// Set calibration parameters:
	428	// if calibration database exists, they are read from database,
	429	// otherwise, they are taken from digitizer.
	430	//
	431	// It is a user responsilibity to open CDB before reconstruction,
	432	// for example:
	433	// AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
	434
	435	//Check if calibration is stored in data base
	436
	437	if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
	438	{
	439	AliCDBEntry entry = (AliCDBEntry)
	440	AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
	441	if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
	442	}
	443
	444	if(!fCalibData)
	445	printf("*********** Calibration parameters not found in CDB! **************");
	446	// AliFatal("Calibration parameters not found in CDB!");
	447
	448
	449	}
	450
	451	//____________________________________________________________________________
	452	Float_t AliJetAODFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
	453	{
	454
	455	// Convert digitized amplitude into energy.
	456	// Calibration parameters are taken from calibration data base for raw data,
	457	// or from digitizer parameters for simulated data.
	458
	459	if(fCalibData){
	460
	461	if (fGeom==0)
	462	printf("*********** Did not get geometry from EMCALLoader *************");
	463	// AliFatal("Did not get geometry from EMCALLoader") ;
	464
	465	Int_t iSupMod = -1;
	466	Int_t nModule = -1;
	467	Int_t nIphi = -1;
	468	Int_t nIeta = -1;
	469	Int_t iphi = -1;
	470	Int_t ieta = -1;
	471
	472	Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
	473	if(!bCell) {
	474	// fGeom->PrintGeometry();
	475	Error("Calibrate()"," Wrong cell id number : %i", AbsId);
	476	assert(0);
	477	}
	478
	479	fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
	480
	481	fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
	482	fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
	483
	484	return -fADCpedestalECA + amp * fADCchannelECA ;
	485
	486	}
487	else //Return energy with default parameters if calibration is not available
488	return -fADCpedestalECA + amp * fADCchannelECA ;
489
490	}
491	*/
492
493