[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetContainer.cxx

//
// Jet Container
//
// Author: M. Verweij

#include <TROOT.h>
#include <TSystem.h>
#include <TInterpreter.h>

#include <TChain.h>
#include <TClonesArray.h>
#include <TList.h>
#include <TObject.h>
#include "AliEmcalJet.h"
#include "AliVEvent.h"
#include "AliLog.h"
#include "AliRhoParameter.h"
#include "AliEMCALGeometry.h"

#include "AliJetContainer.h"

ClassImp(AliJetContainer)

//________________________________________________________________________
AliJetContainer::AliJetContainer():
  AliEmcalContainer("AliJetContainer"),
  fJetAcceptanceType(kUser),
  fJetRadius(0),
  fRhoName(),
  fPtBiasJetTrack(0),
  fPtBiasJetClus(0),
  fJetPtCut(1),
  fJetAreaCut(-1),
  fAreaEmcCut(0),
  fJetMinEta(-0.9),
  fJetMaxEta(0.9),
  fJetMinPhi(-10),
  fJetMaxPhi(10),
  fMaxClusterPt(1000),
  fMaxTrackPt(100),
  fLeadingHadronType(0),
  fNLeadingJets(1),
  fJetBitMap(0),
  fRho(0),
  fGeom(0),
  fRunNumber(0)
{
  // Default constructor.

}

//________________________________________________________________________
AliJetContainer::AliJetContainer(const char *name):
  AliEmcalContainer(name),
  fJetAcceptanceType(kUser),
  fJetRadius(0),
  fRhoName(),
  fPtBiasJetTrack(0),
  fPtBiasJetClus(0),
  fJetPtCut(1),
  fJetAreaCut(-1),
  fAreaEmcCut(0),
  fJetMinEta(-0.9),
  fJetMaxEta(0.9),
  fJetMinPhi(-10),
  fJetMaxPhi(10),
  fMaxClusterPt(1000),
  fMaxTrackPt(100),
  fLeadingHadronType(0),
  fNLeadingJets(1),
  fJetBitMap(0),
  fRho(0),
  fGeom(0),
  fRunNumber(0)
{
  // Standard constructor.

}

//________________________________________________________________________
AliJetContainer::~AliJetContainer()
{
  // Destructor.
}

//________________________________________________________________________
void AliJetContainer::SetJetArray(AliVEvent *event) 
{
  // Set jet array

  SetArray(event, "AliEmcalJet");

  if(fJetAcceptanceType==kTPC) {
    AliDebug(2,Form("%s: set TPC acceptance cuts",GetName()));
    SetJetEtaPhiTPC();
  }
  else if(fJetAcceptanceType==kEMCAL) {
    AliDebug(2,Form("%s: set EMCAL acceptance cuts",GetName()));
    SetJetEtaPhiEMCAL();
 }
}


//________________________________________________________________________
void AliJetContainer::SetEMCALGeometry() {
  fGeom = AliEMCALGeometry::GetInstance();
  if (!fGeom) {
    AliError(Form("%s: Can not create geometry", GetName()));
    return;
  }
}

//________________________________________________________________________
void AliJetContainer::LoadRho(AliVEvent *event)
{
  //Load rho

  if (!fRhoName.IsNull() && !fRho) {
    fRho = dynamic_cast<AliRhoParameter*>(event->FindListObject(fRhoName));
    if (!fRho) {
      AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
      return;
    }
  }
}

//________________________________________________________________________
AliEmcalJet* AliJetContainer::GetJet(Int_t i) const {

  //Get i^th jet in array

  if(i<0 || i>fClArray->GetEntriesFast()) return 0;
  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fClArray->At(i));
  return jet;

}


//________________________________________________________________________
Double_t AliJetContainer::GetJetPtCorr(Int_t i) const {
  AliEmcalJet *jet = GetJet(i);

  return jet->Pt() - fRho->GetVal()*jet->Area();
}

//________________________________________________________________________
AliEmcalJet* AliJetContainer::GetAcceptJet(Int_t i) const {

  //Only return jet if is accepted

  AliEmcalJet *jet = GetJet(i);
  if(!AcceptJet(jet)) return 0;

  return jet;
}

//________________________________________________________________________
Bool_t AliJetContainer::AcceptBiasJet(AliEmcalJet *jet) const
{ 
  // Accept jet with a bias.

  if (fLeadingHadronType == 0) {
    if (jet->MaxTrackPt() < fPtBiasJetTrack) return kFALSE;
  }
  else if (fLeadingHadronType == 1) {
    if (jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
  }
  else {
    if (jet->MaxTrackPt() < fPtBiasJetTrack && jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
  }

  return kTRUE;


}

//________________________________________________________________________
Bool_t AliJetContainer::AcceptJet(AliEmcalJet *jet) const
{   
  // Return true if jet is accepted.
  if (!jet)
    return kFALSE;

  if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap)
    return kFALSE;

  if (jet->Pt() <= fJetPtCut) 
    return kFALSE;
 
  if (jet->Area() <= fJetAreaCut) 
    return kFALSE;

  if (jet->AreaEmc()<fAreaEmcCut)
    return kFALSE;

  if (!AcceptBiasJet(jet))
    return kFALSE;
  
  if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
    return kFALSE;

  Double_t jetPhi = jet->Phi();
  Double_t jetEta = jet->Eta();
  
  if (fJetMinPhi < 0) // if limits are given in (-pi, pi) range
    jetPhi -= TMath::Pi() * 2;

  return (Bool_t)(jetEta > fJetMinEta && jetEta < fJetMaxEta && jetPhi > fJetMinPhi && jetPhi < fJetMaxPhi);
}

//________________________________________________________________________
Double_t AliJetContainer::GetLeadingHadronPt(AliEmcalJet *jet) const
{
  if (fLeadingHadronType == 0)       // charged leading hadron
    return jet->MaxTrackPt();
  else if (fLeadingHadronType == 1)  // neutral leading hadron
    return jet->MaxClusterPt();
  else                               // charged or neutral
    return jet->MaxPartPt();
}


//________________________________________________________________________
void AliJetContainer::SetJetEtaPhiEMCAL()
{
  //Set default cuts for full jets

  if(!fGeom) SetEMCALGeometry();
  if(fGeom) {
    SetJetEtaLimits(fGeom->GetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius);

    if(fRunNumber>=177295 && fRunNumber<=197470) //small SM masked in 2012 and 2013
      SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
    else
      SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius);

  }
  else {
    AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
    SetJetEtaLimits(-0.7+fJetRadius,0.7-fJetRadius);
    SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
  }

}

//________________________________________________________________________
void AliJetContainer::SetJetEtaPhiTPC()
{
  //Set default cuts for charged jets

  SetJetEtaLimits(-0.9+fJetRadius, 0.9-fJetRadius);
  SetJetPhiLimits(-10, 10);

}

//________________________________________________________________________
void AliJetContainer::ResetCuts() 
{
  // Reset cuts to default values

  fPtBiasJetTrack = 0;
  fPtBiasJetClus = 0;
  fJetPtCut = 1;
  fJetAreaCut = -1;
  fAreaEmcCut = 0;
  fJetMinEta = -0.9;
  fJetMaxEta = 0.9;
  fJetMinPhi = -10;
  fJetMaxPhi = 10;
  fMaxClusterPt = 1000;
  fMaxTrackPt = 100;
  fLeadingHadronType = 0;

}
Commit	Line	Data
421b12a3	1	//
	2	// Jet Container
	3	//
	4	// Author: M. Verweij
	5
	6	#include <TROOT.h>
	7	#include <TSystem.h>
	8	#include <TInterpreter.h>
	9
	10	#include <TChain.h>
	11	#include <TClonesArray.h>
	12	#include <TList.h>
	13	#include <TObject.h>
	14	#include "AliEmcalJet.h"
	15	#include "AliVEvent.h"
	16	#include "AliLog.h"
	17	#include "AliRhoParameter.h"
	18	#include "AliEMCALGeometry.h"
	19
	20	#include "AliJetContainer.h"
	21
	22	ClassImp(AliJetContainer)
	23
	24	//________________________________________________________________________
	25	AliJetContainer::AliJetContainer():
	26	AliEmcalContainer("AliJetContainer"),
	27	fJetAcceptanceType(kUser),
	28	fJetRadius(0),
	29	fRhoName(),
	30	fPtBiasJetTrack(0),
	31	fPtBiasJetClus(0),
	32	fJetPtCut(1),
	33	fJetAreaCut(-1),
	34	fAreaEmcCut(0),
	35	fJetMinEta(-0.9),
	36	fJetMaxEta(0.9),
	37	fJetMinPhi(-10),
	38	fJetMaxPhi(10),
	39	fMaxClusterPt(1000),
	40	fMaxTrackPt(100),
	41	fLeadingHadronType(0),
	42	fNLeadingJets(1),
	43	fJetBitMap(0),
	44	fRho(0),
56f370b0	45	fGeom(0),
56f370b0	46	fRunNumber(0)
421b12a3	47	{
	48	// Default constructor.
	49
	50	}
	51
	52	//________________________________________________________________________
	53	AliJetContainer::AliJetContainer(const char *name):
	54	AliEmcalContainer(name),
	55	fJetAcceptanceType(kUser),
	56	fJetRadius(0),
	57	fRhoName(),
	58	fPtBiasJetTrack(0),
	59	fPtBiasJetClus(0),
	60	fJetPtCut(1),
	61	fJetAreaCut(-1),
	62	fAreaEmcCut(0),
	63	fJetMinEta(-0.9),
	64	fJetMaxEta(0.9),
	65	fJetMinPhi(-10),
	66	fJetMaxPhi(10),
	67	fMaxClusterPt(1000),
	68	fMaxTrackPt(100),
	69	fLeadingHadronType(0),
	70	fNLeadingJets(1),
	71	fJetBitMap(0),
	72	fRho(0),
56f370b0	73	fGeom(0),
56f370b0	74	fRunNumber(0)
421b12a3	75	{
	76	// Standard constructor.
	77
	78	}
	79
	80	//________________________________________________________________________
	81	AliJetContainer::~AliJetContainer()
	82	{
	83	// Destructor.
	84	}
	85
	86	//________________________________________________________________________
	87	void AliJetContainer::SetJetArray(AliVEvent *event)
	88	{
	89	// Set jet array
	90
421b12a3	91	SetArray(event, "AliEmcalJet");
421b12a3	92
56f370b0	93	if(fJetAcceptanceType==kTPC) {
56f370b0	94	AliDebug(2,Form("%s: set TPC acceptance cuts",GetName()));
421b12a3	95	SetJetEtaPhiTPC();
56f370b0	96	}
	97	else if(fJetAcceptanceType==kEMCAL) {
	98	AliDebug(2,Form("%s: set EMCAL acceptance cuts",GetName()));
421b12a3	99	SetJetEtaPhiEMCAL();
56f370b0	100	}
421b12a3	101	}
421b12a3	102
56f370b0	103
421b12a3	104	//________________________________________________________________________
	105	void AliJetContainer::SetEMCALGeometry() {
	106	fGeom = AliEMCALGeometry::GetInstance();
	107	if (!fGeom) {
	108	AliError(Form("%s: Can not create geometry", GetName()));
	109	return;
	110	}
	111	}
	112
	113	//________________________________________________________________________
	114	void AliJetContainer::LoadRho(AliVEvent *event)
	115	{
	116	//Load rho
	117
	118	if (!fRhoName.IsNull() && !fRho) {
	119	fRho = dynamic_cast<AliRhoParameter*>(event->FindListObject(fRhoName));
	120	if (!fRho) {
	121	AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
	122	return;
	123	}
	124	}
	125	}
	126
	127	//________________________________________________________________________
	128	AliEmcalJet* AliJetContainer::GetJet(Int_t i) const {
	129
	130	//Get i^th jet in array
	131
	132	if(i<0 \|\| i>fClArray->GetEntriesFast()) return 0;
	133	AliEmcalJet jet = static_cast<AliEmcalJet>(fClArray->At(i));
	134	return jet;
	135
	136	}
	137
	138
	139	//________________________________________________________________________
	140	Double_t AliJetContainer::GetJetPtCorr(Int_t i) const {
	141	AliEmcalJet *jet = GetJet(i);
	142
	143	return jet->Pt() - fRho->GetVal()*jet->Area();
	144	}
	145
	146	//________________________________________________________________________
	147	AliEmcalJet* AliJetContainer::GetAcceptJet(Int_t i) const {
	148
	149	//Only return jet if is accepted
	150
	151	AliEmcalJet *jet = GetJet(i);
	152	if(!AcceptJet(jet)) return 0;
	153
	154	return jet;
	155	}
	156
	157	//________________________________________________________________________
	158	Bool_t AliJetContainer::AcceptBiasJet(AliEmcalJet *jet) const
	159	{
	160	// Accept jet with a bias.
	161
	162	if (fLeadingHadronType == 0) {
	163	if (jet->MaxTrackPt() < fPtBiasJetTrack) return kFALSE;
	164	}
	165	else if (fLeadingHadronType == 1) {
	166	if (jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
	167	}
168	else {
169	if (jet->MaxTrackPt() < fPtBiasJetTrack && jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
170	}
171
172	return kTRUE;
173
174
175	}
176
177	//________________________________________________________________________
178	Bool_t AliJetContainer::AcceptJet(AliEmcalJet *jet) const
179	{
180	// Return true if jet is accepted.
181	if (!jet)
182	return kFALSE;
183
184	if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap)
185	return kFALSE;
186
187	if (jet->Pt() <= fJetPtCut)
188	return kFALSE;
189
190	if (jet->Area() <= fJetAreaCut)
191	return kFALSE;
192
193	if (jet->AreaEmc()<fAreaEmcCut)
194	return kFALSE;
195
196	if (!AcceptBiasJet(jet))
197	return kFALSE;
198
199	if (jet->MaxTrackPt() > fMaxTrackPt \|\| jet->MaxClusterPt() > fMaxClusterPt)
200	return kFALSE;
201
202	Double_t jetPhi = jet->Phi();
203	Double_t jetEta = jet->Eta();
204
205	if (fJetMinPhi < 0) // if limits are given in (-pi, pi) range
206	jetPhi -= TMath::Pi() * 2;
207
208	return (Bool_t)(jetEta > fJetMinEta && jetEta < fJetMaxEta && jetPhi > fJetMinPhi && jetPhi < fJetMaxPhi);
209	}
210
211	//________________________________________________________________________
212	Double_t AliJetContainer::GetLeadingHadronPt(AliEmcalJet *jet) const
213	{
214	if (fLeadingHadronType == 0) // charged leading hadron
215	return jet->MaxTrackPt();
216	else if (fLeadingHadronType == 1) // neutral leading hadron
217	return jet->MaxClusterPt();
218	else // charged or neutral
219	return jet->MaxPartPt();
220	}
221
222
223	//________________________________________________________________________
224	void AliJetContainer::SetJetEtaPhiEMCAL()
225	{
226	//Set default cuts for full jets
227
228	if(!fGeom) SetEMCALGeometry();
229	if(fGeom) {
230	SetJetEtaLimits(fGeom->GetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius);
56f370b0	231
	232	if(fRunNumber>=177295 && fRunNumber<=197470) //small SM masked in 2012 and 2013
	233	SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
	234	else
	235	SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius);
	236
421b12a3	237	}
421b12a3	238	else {
56f370b0	239	AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
421b12a3	240	SetJetEtaLimits(-0.7+fJetRadius,0.7-fJetRadius);
	241	SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
	242	}
	243
	244	}
	245
	246	//________________________________________________________________________
	247	void AliJetContainer::SetJetEtaPhiTPC()
	248	{
56f370b0	249	//Set default cuts for charged jets
421b12a3	250
56f370b0	251	SetJetEtaLimits(-0.9+fJetRadius, 0.9-fJetRadius);
421b12a3	252	SetJetPhiLimits(-10, 10);
	253
	254	}
	255
	256	//________________________________________________________________________
	257	void AliJetContainer::ResetCuts()
	258	{
	259	// Reset cuts to default values
	260
	261	fPtBiasJetTrack = 0;
	262	fPtBiasJetClus = 0;
	263	fJetPtCut = 1;
	264	fJetAreaCut = -1;
	265	fAreaEmcCut = 0;
	266	fJetMinEta = -0.9;
	267	fJetMaxEta = 0.9;
	268	fJetMinPhi = -10;
	269	fJetMaxPhi = 10;
	270	fMaxClusterPt = 1000;
	271	fMaxTrackPt = 100;
	272	fLeadingHadronType = 0;
	273
	274	}