[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 with EMCal information 
// Called by ESD reader for jet analysis
// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
//------------------------------------------------------------- 

// --- Standard library ---
#include <Riostream.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 "AliESD.h"
#include "AliESDEvent.h"
#include "AliJetDummyGeo.h"
#include "AliESDCaloCluster.h"
#include "AliJetUnitArray.h"
#include "AliJetFillUnitArrayEMCalDigits.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),
    fDebug(0),
    fNCEMCAL(0),
    fNCPHOS(0),
    fNCCalo(0),
    fTPCGrid(0x0),
    fEMCalGrid(0x0),
    fReaderHeader(0x0),
    fMomentumArray(0x0),
    fUnitArray(0x0),
    fRefArray(0x0),
    fGeom(0x0),
    fClus(0x0),
    fNDigitEmcal(0),
    fNDigitEmcalCut(0)
{
  // constructor
}

//_____________________________________________________________________________
AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(AliESDEvent */*esd*/)
  : 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),
    fDebug(0),
    fNCEMCAL(0),
    fNCPHOS(0),
    fNCCalo(0),
    fTPCGrid(0x0),
    fEMCalGrid(0x0),
    fReaderHeader(0x0),
    fMomentumArray(0x0),
    fUnitArray(0x0),
    fRefArray(0x0),
    fGeom(0x0),
    fClus(0x0),
    fNDigitEmcal(0),
    fNDigitEmcalCut(0)
{
  // constructor
}


//____________________________________________________________________________
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;

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

}

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

//_____________________________________________________________________________
void AliJetFillUnitArrayEMCalDigits::Exec(Option_t* /*option*/)
{
  //
  // Main method.
  // Explain

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

  // Init parameters
  InitParameters();

  // Get number of clusters from EMCAL
  
  Int_t   nDigitTot      = 0;
  Int_t   goodDigit      = 0;
  Int_t   beg            = 0;
  Int_t   end            = 0;
  Float_t ptMin          = fReaderHeader->GetPtCut();

  // Loop over calo clusters 
  //------------------------------------------------------------------
  Int_t type = 0; 
  Int_t index = 0;

  // Total number of EMCAL cluster
  end =  fESD->GetNumberOfCaloClusters();

  for(Int_t j = beg; j < end; j++) {
      fClus = fESD->GetCaloCluster(j);
      if(!fClus->IsEMCAL()) continue;

      type = fClus->GetClusterType(); 
      index = fClus->GetID();
      nDigitTot = fClus->GetNumberOfDigits();
      
      // Keep clusters or pseudo clusters
      if (type != AliESDCaloCluster::kEMCALClusterv1) continue;
      //      if (type != AliESDCaloCluster::kPseudoCluster) continue;

      // Get the digit index and the digit information
      //============================================================
	  
      // Get number of digits in a cluster
      Int_t nD = fClus->GetNumberOfDigits();

      TArrayS *digID = fClus->GetDigitIndex();      
      TArrayS *digEnergy = fClus->GetDigitAmplitude();
      Float_t *digitEnergy = new Float_t[nD];      
      //      Float_t digitEn = 0.;
	  
      // Loop over digits
      for(Int_t k=0; k<nD; k++) {
	
	// Convert energy in GeV
	Int_t idF = (Int_t)digID->At(k);
	// Calibration for an energy in GeV
	digitEnergy[k] = (Float_t)digEnergy->At(k)/500.; 

	// Second method to extract eta, phi positions of a digit
	//=================================================================

	Float_t etaD=-10., phiD=-10.;
	fGeom->EtaPhiFromIndex(idF,etaD,phiD); 
	//	    fEMCalGrid->GetEtaPhiFromIndex2(idF,phiD,etaD);
	phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);

	Float_t etDigit = digitEnergy[k]*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));

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

	Float_t unitEnergy = 0.;
	Bool_t ok = kFALSE;
	unitEnergy = uArray->GetUnitEnergy();
	if(unitEnergy==0){
	  fRefArray->Add(uArray);
	  fNDigitEmcal++;
	  ok = kTRUE;
	}
	uArray->SetUnitEnergy(unitEnergy+etDigit);
	// Put a cut flag
	if(uArray->GetUnitEnergy()<ptMin)
	  uArray->SetUnitCutFlag(kPtSmaller);
	else {
	  uArray->SetUnitCutFlag(kPtHigher);
	  if(ok) fNDigitEmcalCut++;
	}
	// Detector flag
	if(unitEnergy>0.) 
	  uArray->SetUnitDetectorFlag(kAll);
	else uArray->SetUnitDetectorFlag(kEmcal);
	
	// This is for jet multiplicity
	uArray->SetUnitClusterID(index);
	
	if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
	
      } // End loop over digits
      
  } // End loop over clusters
  
  fNIn += goodDigit;

}

//_____________________________________________________________________________
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	**************************************************************************/
	16
	17	//-------------------------------------------------------------
	18	// Fill Unit Array with EMCal information
	19	// Called by ESD reader for jet analysis
	20	// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
	21	//-------------------------------------------------------------
	22
	23	// --- Standard library ---
	24	#include <Riostream.h>
	25
	26	// --- ROOT system ---
	27	#include <TSystem.h>
	28	#include <TLorentzVector.h>
	29	#include <TRefArray.h>
	30	#include "TTask.h"
	31	#include <TGeoManager.h>
	32	#include <TMath.h>
	33	#include <TClonesArray.h>
	34	#include <TArrayS.h>
	35
	36	// --- AliRoot header files ---
	37	#include "AliJetFinder.h"
	38	#include "AliJetReaderHeader.h"
	39	#include "AliJetReader.h"
	40	#include "AliJetESDReader.h"
	41	#include "AliJetESDReaderHeader.h"
	42	//#include "AliESD.h"
	43	#include "AliESDEvent.h"
	44	#include "AliJetDummyGeo.h"
	45	#include "AliESDCaloCluster.h"
	46	#include "AliJetUnitArray.h"
	47	#include "AliJetFillUnitArrayEMCalDigits.h"
ee7de0dd	48
	49	ClassImp(AliJetFillUnitArrayEMCalDigits)
	50
	51	//_____________________________________________________________________________
	52	AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits()
	53	: TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
	54	fESD(0),
	55	fNumUnits(0),
	56	fEtaMinCal(0.),
	57	fEtaMaxCal(0.),
	58	fPhiMinCal(0.),
	59	fPhiMaxCal(0.),
	60	fNIn(0),
	61	fOpt(0),
	62	fDebug(0),
	63	fNCEMCAL(0),
	64	fNCPHOS(0),
	65	fNCCalo(0),
	66	fTPCGrid(0x0),
	67	fEMCalGrid(0x0),
	68	fReaderHeader(0x0),
	69	fMomentumArray(0x0),
	70	fUnitArray(0x0),
	71	fRefArray(0x0),
	72	fGeom(0x0),
	73	fClus(0x0),
	74	fNDigitEmcal(0),
	75	fNDigitEmcalCut(0)
	76	{
	77	// constructor
	78	}
	79
	80	//_____________________________________________________________________________
	81	AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(AliESDEvent /esd*/)
	82	: TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
	83	fESD(0),
	84	fNumUnits(0),
	85	fEtaMinCal(0.),
	86	fEtaMaxCal(0.),
	87	fPhiMinCal(0.),
	88	fPhiMaxCal(0.),
	89	fNIn(0),
	90	fOpt(0),
	91	fDebug(0),
	92	fNCEMCAL(0),
	93	fNCPHOS(0),
	94	fNCCalo(0),
	95	fTPCGrid(0x0),
	96	fEMCalGrid(0x0),
	97	fReaderHeader(0x0),
	98	fMomentumArray(0x0),
	99	fUnitArray(0x0),
	100	fRefArray(0x0),
	101	fGeom(0x0),
	102	fClus(0x0),
	103	fNDigitEmcal(0),
	104	fNDigitEmcalCut(0)
	105	{
	106	// constructor
	107	}
	108
	109
	110	//____________________________________________________________________________
	111	void AliJetFillUnitArrayEMCalDigits::InitParameters()
112	{
113	fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL
114
115	fEtaMinCal = fGeom->GetArm1EtaMin();
116	fEtaMaxCal = fGeom->GetArm1EtaMax();
117	fPhiMinCal = fGeom->GetArm1PhiMin();
118	fPhiMaxCal = fGeom->GetArm1PhiMax();
119	fClus = 0;
120
121	if(fDebug>1) printf("\n EMCAL parameters initiated ! \n");
122
123	}
124
125	//_____________________________________________________________________________
126	AliJetFillUnitArrayEMCalDigits::~AliJetFillUnitArrayEMCalDigits()
127	{
128	// destructor
129	}
130
131	//_____________________________________________________________________________
132	void AliJetFillUnitArrayEMCalDigits::Exec(Option_t* /option/)
133	{
134	//
135	// Main method.
136	// Explain
137
138	fDebug = fReaderHeader->GetDebug();
139	fOpt = fReaderHeader->GetDetector();
140
141	// Init parameters
142	InitParameters();
143
144	// Get number of clusters from EMCAL
145
146	Int_t nDigitTot = 0;
147	Int_t goodDigit = 0;
148	Int_t beg = 0;
149	Int_t end = 0;
150	Float_t ptMin = fReaderHeader->GetPtCut();
151
152	// Loop over calo clusters
153	//------------------------------------------------------------------
154	Int_t type = 0;
155	Int_t index = 0;
156
157	// Total number of EMCAL cluster
158	end = fESD->GetNumberOfCaloClusters();
159
160	for(Int_t j = beg; j < end; j++) {
161	fClus = fESD->GetCaloCluster(j);
162	if(!fClus->IsEMCAL()) continue;
163
164	type = fClus->GetClusterType();
165	index = fClus->GetID();
166	nDigitTot = fClus->GetNumberOfDigits();
167
168	// Keep clusters or pseudo clusters
8ada0ffe	169	if (type != AliESDCaloCluster::kEMCALClusterv1) continue;
ee7de0dd	170	// if (type != AliESDCaloCluster::kPseudoCluster) continue;
	171
	172	// Get the digit index and the digit information
	173	//============================================================
	174
	175	// Get number of digits in a cluster
	176	Int_t nD = fClus->GetNumberOfDigits();
	177
90dc70f0	178	TArrayS *digID = fClus->GetDigitIndex();
	179	TArrayS *digEnergy = fClus->GetDigitAmplitude();
	180	Float_t *digitEnergy = new Float_t[nD];
ee7de0dd	181	// Float_t digitEn = 0.;
	182
	183	// Loop over digits
	184	for(Int_t k=0; k<nD; k++) {
	185
	186	// Convert energy in GeV
90dc70f0	187	Int_t idF = (Int_t)digID->At(k);
ee7de0dd	188	// Calibration for an energy in GeV
90dc70f0	189	digitEnergy[k] = (Float_t)digEnergy->At(k)/500.;
ee7de0dd	190
	191	// Second method to extract eta, phi positions of a digit
	192	//=================================================================
	193
	194	Float_t etaD=-10., phiD=-10.;
	195	fGeom->EtaPhiFromIndex(idF,etaD,phiD);
	196	// fEMCalGrid->GetEtaPhiFromIndex2(idF,phiD,etaD);
	197	phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
	198
	199	Float_t etDigit = digitEnergy[k]*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
	200
	201	AliJetUnitArray uArray = (AliJetUnitArray)fUnitArray->At(idF);
	202	if(uArray->GetUnitEnergy() == 0.) goodDigit++;
	203
	204	Float_t unitEnergy = 0.;
	205	Bool_t ok = kFALSE;
	206	unitEnergy = uArray->GetUnitEnergy();
	207	if(unitEnergy==0){
	208	fRefArray->Add(uArray);
	209	fNDigitEmcal++;
	210	ok = kTRUE;
	211	}
	212	uArray->SetUnitEnergy(unitEnergy+etDigit);
	213	// Put a cut flag
	214	if(uArray->GetUnitEnergy()<ptMin)
	215	uArray->SetUnitCutFlag(kPtSmaller);
	216	else {
	217	uArray->SetUnitCutFlag(kPtHigher);
	218	if(ok) fNDigitEmcalCut++;
	219	}
	220	// Detector flag
	221	if(unitEnergy>0.)
	222	uArray->SetUnitDetectorFlag(kAll);
	223	else uArray->SetUnitDetectorFlag(kEmcal);
	224
	225	// This is for jet multiplicity
	226	uArray->SetUnitClusterID(index);
	227
	228	if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
	229
	230	} // End loop over digits
	231
	232	} // End loop over clusters
	233
	234	fNIn += goodDigit;
	235
	236	}
	237
	238	//_____________________________________________________________________________
	239	Float_t AliJetFillUnitArrayEMCalDigits::EtaToTheta(Float_t arg)
	240	{
	241	// return (180./TMath::Pi())2.atan(exp(-arg));
	242	return 2.*atan(exp(-arg));
	243
	244
	245	}