#include "AliJet.h"
ClassImp(AliJet)
-
-
-////////////////////////////////////////////////////////////////////////
+
AliJet::AliJet()
{
- //
// Constructor
- //
fJets = new TClonesArray("TLorentzVector",1000);
fNInput=0;
fNJets=0;
+ fEtAvg = 0;
fInJet=TArrayI();
+ fPtIn=TArrayF();
+ fEtaIn=TArrayF();
+ fPhiIn=TArrayF();
fPtFromSignal=TArrayF();
fMultiplicities=TArrayI();
+ fNCells=TArrayI();
}
////////////////////////////////////////////////////////////////////////
Bool_t AliJet::OutOfRange(Int_t i, const char *s) const
{
- // checks if i is a valid index. s= name of calling method
+ // checks if i is a valid index. s = name of calling method
if (i >= fNJets || i < 0) {
cout << s << " Index " << i << " out of range" << endl;
return kTRUE;
{
// returns i-jet
if (OutOfRange(i, "AliJet::GetJet:")) return 0;
-
TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
return lv;
}
////////////////////////////////////////////////////////////////////////
-Int_t AliJet::GetMultiplicity(Int_t i)
+Int_t AliJet::GetMultiplicity(Int_t i) const
{
// gets multiplicity of i-jet
if (OutOfRange(i, "AliJet::GetMultiplicity:")) return 0;
////////////////////////////////////////////////////////////////////////
+Int_t AliJet::GetNCell(Int_t i) const
+{
+ // gets number of cell of i-jet
+ if (OutOfRange(i, "AliJet::GetNCell:")) return 0;
+ return fNCells[i];
+}
+
+////////////////////////////////////////////////////////////////////////
+
Double_t AliJet::GetPx(Int_t i)
{
// Get Px component of jet i
if (OutOfRange(i, "AliJet::GetPx:")) return -1e30;
-
TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
return lv->Px();
}
{
// Get Py component of jet i
if (OutOfRange(i, "AliJet::GetPy:")) return -1e30;
-
TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
return lv->Py();
}
{
// Get Pz component of jet i
if (OutOfRange(i, "AliJet::GetPz:")) return -1e30;
-
TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
return lv->Pz();
}
{
// Get momentum of jet i
if (OutOfRange(i, "AliJet::GetP:")) return -1e30;
-
TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
return lv->P();
}
{
// Get energy of jet i
if (OutOfRange(i, "AliJet::GetE:")) return -1e30;
-
TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
return lv->E();
}
{
// Get transverse momentum of jet i
if (OutOfRange(i, "AliJet::GetPt:")) return -1e30;
-
TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
return lv->Pt();
}
{
// Get eta of jet i
if (OutOfRange(i, "AliJet::GetEta:")) return -1e30;
-
TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
return lv->Eta();
}
{
// Get phi of jet i
if (OutOfRange(i, "AliJet::GetPhi:")) return -1e30;
-
TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
return ( (lv->Phi() < 0) ? (lv->Phi()) + 2. * TMath::Pi() : lv->Phi());
}
{
// Get theta of jet i
if (OutOfRange(i, "AliJet::GetTheta:")) return -1e30;
-
TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
return lv->Theta();
}
{
// Get invariant mass of jet i
if (OutOfRange(i, "AliJet::GetMass:")) return -1e30;
-
TLorentzVector *lv = (TLorentzVector*) fJets->At(i);
return lv->M();
}
void AliJet::SetInJet(Int_t* j)
{
// set information of which input object belongs
- // to each jet. filled in by AliJetFinder
+ // to each jet. If zero, object was not assigned to
+ // a jet, if n,positive, it was assiged to jet n
+ // if n, negative, it is within cone of jet n, but
+ // it did not passed the user cuts. filled in by AliJetFinder
+
if (fNInput>0) fInJet.Set(fNInput, j);
}
////////////////////////////////////////////////////////////////////////
+void AliJet::SetEtaIn(Float_t* r)
+{
+ if (fNInput>0) fEtaIn.Set(fNInput, r);
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJet::SetPtIn(Float_t* pt)
+{
+ if (fNInput>0) fPtIn.Set(fNInput, pt);
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJet::SetPhiIn(Float_t* x)
+{
+ if (fNInput>0) fPhiIn.Set(fNInput, x);
+}
+
+////////////////////////////////////////////////////////////////////////
+
void AliJet::SetPtFromSignal(Float_t* p)
{
// set information of percentage of pt of jets
////////////////////////////////////////////////////////////////////////
+void AliJet::SetNCells(Int_t* n)
+{
+ if (fNJets>0) fNCells.Set(fNJets, n);
+}
+
+////////////////////////////////////////////////////////////////////////
+
void AliJet::ClearJets(Option_t *option)
{
// reset all values
fNJets=0;
fMultiplicities.Set(0);
fInJet.Set(0);
+ fPtFromSignal.Set(0);
+ fPhiIn.Set(0);
+ fEtaIn.Set(0);
+ fPtIn.Set(0);
+ fNCells.Set(0);
}
////////////////////////////////////////////////////////////////////////
TClonesArray* GetJets() const {return fJets;}
TArrayI GetInJet() const {return fInJet;}
TArrayI GetMultiplicities() const {return fMultiplicities;}
+ TArrayI GetNCells() const {return fNCells;}
TArrayF GetPtFromSignal() const {return fPtFromSignal;}
+ TArrayF GetEtaIn() const { return fEtaIn; }
+ TArrayF GetPhiIn() const { return fPhiIn; }
+ TArrayF GetPtIn() const { return fPtIn; }
+ Double_t GetEtAvg() const { return fEtAvg; }
TLorentzVector* GetJet(Int_t i);
- Int_t GetMultiplicity(Int_t i);
+ Int_t GetMultiplicity(Int_t i) const;
+ Int_t GetNCell(Int_t i) const;
Double_t GetPx(Int_t i);
Double_t GetPy(Int_t i);
Double_t GetPz(Int_t i);
// Setters
void SetNinput(Int_t i) {fNInput = i;}
void AddJet(Double_t px, Double_t py, Double_t pz, Double_t e);
- void SetInJet(Int_t* j);
void SetMultiplicities(Int_t* m);
+ void SetNCells(Int_t* n);
void SetPtFromSignal(Float_t* p);
-
+ void SetEtaIn(Float_t* eta);
+ void SetPhiIn(Float_t* phi);
+ void SetPtIn(Float_t* pt);
+ void SetInJet(Int_t* idx);
+ void SetEtAvg(Double_t et) { fEtAvg = et; }
+
// others
Bool_t OutOfRange(Int_t i, const char *s) const;
void ClearJets(Option_t *option="");
Int_t fNInput; // number of input objects
Int_t fNJets; // number of jets found
+ Double_t fEtAvg; // average background et per cell
+
TArrayI fInJet; // i-input object belongs to k-jet
TArrayI fMultiplicities; // Multiplicity of each jet
+ TArrayI fNCells; // Number of cells in jet
TArrayF fPtFromSignal; // percentage of pt from signal
TClonesArray* fJets; // 4-momenta of jets
+ TArrayF fEtaIn; // Arrays of input particles kine:Eta
+ TArrayF fPhiIn; // Arrays of input particles kine:Phi
+ TArrayF fPtIn; // Arrays of input particles kine:Pt
+
ClassDef(AliJet,1)
};
//---------------------------------------------------------------------
// JetAnalysis class
-// Analyse Jets
-// Author: andreas.morsch@cern.ch
+// Analyse Jets (already found jets)
+// Authors: andreas.morsch@cern.ch, jgcn@mail.cern.ch
+// mercedes.lopez.noriega@cern.ch
//---------------------------------------------------------------------
+#include <Riostream.h>
#include "AliJetAnalysis.h"
ClassImp(AliJetAnalysis)
-
-
-////////////////////////////////////////////////////////////////////////
+
+// root
+#include <Riostream.h>
#include <TH1F.h>
+#include <TH2F.h>
#include <TProfile.h>
-#include <TCanvas.h>
#include <TFile.h>
#include <TTree.h>
#include <TStyle.h>
#include <TSystem.h>
#include <TLorentzVector.h>
-
+#include <TMath.h>
+// aliroot
#include "AliJetProductionDataPDC2004.h"
#include "AliJet.h"
+#include "AliJetKineReaderHeader.h"
#include "AliJetESDReaderHeader.h"
#include "AliUA1JetHeader.h"
#include "AliLeading.h"
+
+
+AliJetAnalysis::AliJetAnalysis()
+{
+ // Constructor
+ fDirectory = 0x0;
+ fBkgdDirectory = 0x0;
+ fFile = "anajets.root";
+ fEventMin = 0;
+ fEventMax = -1;
+ fRunMin = 0;
+ fRunMax = 11;
+ fminMult = 0;
+ fPercentage = -1.0;
+ fEfactor = 1.0;
+ SetReaderHeader();
+
+ // for Analize()
+ fPythia = kFALSE;
+ fDoPart = kTRUE;
+ fDoGenJ = kTRUE;
+ fDoRecJ = kTRUE;
+ fDoKine = kTRUE;
+ fDoCorr = kTRUE;
+ fDoCorr50 = kFALSE;
+ fDoShap = kTRUE;
+ fDoFrag = kTRUE;
+ fDoTrig = kTRUE;
+ fDoJt = kTRUE;
+ fDodNdxi = kTRUE;
+ fDoBkgd = kFALSE;
+
+ fPart = 0;
+ fGenJ = 0;
+ fRecJ = 0;
+ fRecB = 0;
+ fWeight = 1.0;
+ fWdEdr = 0.0;
+ fPartPtCut = 0.0;
+ fWJt = 0.0;
+ fWdNdxi = 0.0;
+
+ // for bkgd
+ fp0 = 0.0;
+ fPtJ = 0.0;
+ fEJ = 0.0;
+ fEtaJ = 0.0;
+ fPhiJ = 0.0;
+ fjv3X = 0.0;
+ fjv3Y = 0.0;
+ fjv3Z = 0.0;
+
+ // kine histos
+ fRKineEneH = 0;
+ fRKinePtH = 0;
+ fRKinePhiH = 0;
+ fRKineEtaH = 0;
+ fGKineEneH = 0;
+ fGKinePtH = 0;
+ fGKinePhiH = 0;
+ fGKineEtaH = 0;
+ fPKineEneH = 0;
+ fPKinePtH = 0;
+ fPKinePhiH = 0;
+ fPKineEtaH = 0;
+
+ // correlation histograms
+ fPGCorrEneH = 0;
+ fPGCorrPtH = 0;
+ fPGCorrEtaH = 0;
+ fPGCorrPhiH = 0;
+ fPRCorrEneH = 0;
+ fPRCorrPtH = 0;
+ fPRCorrEtaH = 0;
+ fPRCorrPhiH = 0;
+ fRGCorrEneH = 0;
+ fRGCorrPtH = 0;
+ fRGCorrEtaH = 0;
+ fRGCorrPhiH = 0;
+
+ // correlation histograms when one particle
+ // has more than 50% of the energy of the jet
+ fPRCorr50EneH = 0;
+ fPRCorr50PtH = 0;
+ fPRCorr50EtaH = 0;
+ fPRCorr50PhiH = 0;
+ fRGCorr50EneH = 0;
+ fRGCorr50PtH = 0;
+ fRGCorr50EtaH = 0;
+ fRGCorr50PhiH = 0;
+
+ // shape histos
+ fWShapR = 0.0;
+ fRShapSelH = 0;
+ fRShapRejH = 0;
+ fRShapAllH = 0;
+
+ // fragmentation function histos
+ fWFragR = 0.0;
+ fRFragSelH = 0;
+ fRFragRejH = 0;
+ fRFragAllH = 0;
+
+ // trigger bias histos
+ fGTriggerEneH = 0;
+ fRTriggerEneH = 0;
+ fGPTriggerEneH = 0;
+ fPTriggerEneH = 0;
+
+ // dE/dr histo and its dr
+ fdEdrH = 0;
+ fdEdrB = 0;
+ fPtEneH2 = 0;
+ fdEdrW = 0;
+ fdrdEdr = 1.0;
+
+ // jt histo and its dr
+ fJtH = 0;
+ fJtB = 0;
+ fJtW = 0;
+ fdrJt = 1.0;
+
+ // dN/dxi histo and its dr
+ fdNdxiH = 0;
+ fdNdxiB = 0;
+ fdNdxiW = 0;
+ fPtEneH = 0;
+ fdrdNdxi = 0.7;
+
+ // initialize weight for dE/dr histo
+ SetdEdrWeight();
+}
+
+AliJetAnalysis::~AliJetAnalysis()
+{
+ // Destructor
+}
+
+
+////////////////////////////////////////////////////////////////////////
+// define histogrames
+
+void AliJetAnalysis::DefineHistograms()
+{
+ // Define the histograms to be filled
+ if (fDoKine) DefineKineH();
+ if (fDoCorr) DefineCorrH();
+ if (fDoCorr50) DefineCorr50H();
+ if (fDoShap) DefineShapH();
+ if (fDoFrag) DefineFragH();
+ if (fDoTrig) DefineTrigH();
+ if (fDoJt) DefineJtH();
+ if (fDodNdxi) DefinedNdxiH();
+}
+
+void AliJetAnalysis::DefineKineH()
+{
+ // leading particle
+ if (fDoPart) {
+ fPKineEneH = new TH1F("PKineEne","Energy of leading particle",50,0,200);
+ SetProperties(fPKineEneH,"Energy (GeV)","Entries");
+ fPKinePtH = new TH1F("PKinePt","Pt of leading particle",50,0,200);
+ SetProperties(fPKinePtH,"P_{T} (GeV)","Entries");
+ fPKinePhiH = new TH1F("PKinePhiH","Azimuthal angle of leading particle",
+ 90,0.,2.0*TMath::Pi());
+ SetProperties(fPKinePhiH,"#phi","Entries");
+ fPKineEtaH = new TH1F("PKineEtaH","Pseudorapidity of leading particle",
+ 40,-1.0,1.0);
+ SetProperties(fPKineEtaH,"#eta","Entries");
+ }
+ // leading generated jet
+ if (fDoGenJ) {
+ fGKineEneH = new TH1F("GKineEne","Energy of generated jet",50,0,200);
+ SetProperties(fGKineEneH,"Energy (GeV)","Entries");
+ fGKinePtH = new TH1F("GKinePt","Pt of generated jet",50,0,200);
+ SetProperties(fGKinePtH,"P_{T} (GeV)","Entries");
+ fGKinePhiH = new TH1F("GKinePhiH","Azimuthal angle of generated jet",
+ 90,0.,2.0*TMath::Pi());
+ SetProperties(fGKinePhiH,"#phi","Entries");
+ fGKineEtaH = new TH1F("GKineEtaH","Pseudorapidity of generated jet",
+ 40,-1.0,1.0);
+ SetProperties(fGKineEtaH,"#eta","Entries");
+ }
+ // leading reconstructed jet
+ if (fDoRecJ) {
+ fRKineEneH = new TH1F("RKineEne","Energy of reconstructed jet",50,0,200);
+ SetProperties(fRKineEneH,"Energy (GeV)","Entries");
+ fRKinePtH = new TH1F("RKinePt","Pt of reconstructed jet",50,0,200);
+ SetProperties(fRKinePtH,"P_{T} (GeV)","Entries");
+ fRKinePhiH = new TH1F("RKinePhiH","Azimuthal angle of reconstructed jet",
+ 90,0.,2.0*TMath::Pi());
+ SetProperties(fRKinePhiH,"#phi","Entries");
+ fRKineEtaH = new TH1F("RKineEtaH","Pseudorapidity of reconstructed jet",
+ 40,-1.0,1.0);
+ SetProperties(fRKineEtaH,"#eta","Entries");
+ }
+}
- AliJetAnalysis::AliJetAnalysis()
+void AliJetAnalysis::DefineCorrH()
{
- // Constructor
- fDirectory = 0x0;
- fEventMin = 0;
- fEventMax = -1;
- fRunMin = 0;
- fRunMax = 11;
+ // correlation
+ if (fDoPart && fDoGenJ) {
+ fPGCorrEneH = new TH2F("PGCorrEne","Energy correlation part-gen jet",
+ 40,0,200,40,0,200);
+ SetProperties(fPGCorrEneH,"Part Energy (GeV)","Gen Jet Energy (GeV)");
+ fPGCorrPtH = new TH2F("PGCorrPt","Pt correlation part-gen jet",
+ 40,0,200,40,0,200);
+ SetProperties(fPGCorrPtH,"Part P_{T} (GeV)","Gen Jet P_{T} (GeV)");
+ fPGCorrEtaH = new TH2F("PGCorrEta","Pseudorapidity correlation part-gen jet",
+ 40,-1.0,1.0,40,-1.0,1.0);
+ SetProperties(fPGCorrEtaH,"Part #eta","Gen Jet #eta");
+ fPGCorrPhiH = new TH2F("PGCorrPhi","Azimuthal angle correlation part-gen jet",
+ 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
+ SetProperties(fPGCorrPhiH,"Part #phi","Gen Jet #phi");
+ }
+ if (fDoPart && fDoRecJ) {
+ fPRCorrEneH = new TH2F("PRCorrEne","Energy correlation part-rec jet",
+ 40,0,200,40,0,200);
+ SetProperties(fPRCorrEneH,"Part Jet Energy (GeV)","Rec Jet Energy (GeV)");
+ fPRCorrPtH = new TH2F("PRCorrPt","Pt correlation part-rec jet",
+ 40,0,200,40,0,200);
+ SetProperties(fPRCorrPtH,"Part Jet P_{T} (GeV)","Rec Jet P_{T} (GeV)");
+ fPRCorrEtaH = new TH2F("PRCorrEta","Pseudorapidity correlation part-rec jet",
+ 40,-1.0,1.0,40,-1.0,1.0);
+ SetProperties(fPRCorrEtaH,"part #eta","Rec Jet #eta");
+ fPRCorrPhiH = new TH2F("PRCorrPhi","Azimuthal angle correlation part-rec jet",
+ 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
+ SetProperties(fPRCorrPhiH,"Part #phi","Rec Jet #phi");
+ }
+ if (fDoGenJ && fDoRecJ) {
+ fRGCorrEneH = new TH2F("RGCorrEne","Energy correlation rec jet-gen jet",
+ 40,0,200,40,0,200);
+ SetProperties(fRGCorrEneH,"Rec Jet Energy (GeV)","Gen Jet Energy (GeV)");
+ fRGCorrPtH = new TH2F("RGCorrPt","Pt correlation rec jet-gen jet",
+ 40,0,200,40,0,200);
+ SetProperties(fRGCorrPtH,"Rec Jet P_{T} (GeV)","Gen Jet P_{T} (GeV)");
+ fRGCorrEtaH = new TH2F("RGCorrEta","Pseudorapidity correlation rec jet-gen jet",
+ 40,-1.0,1.0,40,-1.0,1.0);
+ SetProperties(fRGCorrEtaH,"Rec Jet #eta","Gen Jet #eta");
+ fRGCorrPhiH = new TH2F("RGCorrPhi","Azimuthal angle correlation rec jet-gen jet",
+ 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
+ SetProperties(fRGCorrPhiH,"Rec Jet #phi","Gen Jet #phi");
+ }
}
-void AliJetAnalysis::Analyse()
+void AliJetAnalysis::DefineCorr50H()
{
- //
- // Some histos
- //
- TH1F::AddDirectory(0);
- TProfile::AddDirectory(0);
+ // correlation
+ if (fDoPart && fDoRecJ) {
+ fPRCorr50EneH = new TH2F("PRCorr50Ene","Energy correlation part-rec jet",
+ 40,0,200,40,0,200);
+ SetProperties(fPRCorr50EneH,"Part Jet Energy (GeV)","Rec Jet Energy (GeV)");
+ fPRCorr50PtH = new TH2F("PRCorr50Pt","Pt correlation part-rec jet",
+ 40,0,200,40,0,200);
+ SetProperties(fPRCorr50PtH,"Part Jet P_{T} (GeV)","Rec Jet P_{T} (GeV)");
+ fPRCorr50EtaH = new TH2F("PRCorr50Eta","Pseudorapidity correlation part-rec jet",
+ 40,-1.0,1.0,40,-1.0,1.0);
+ SetProperties(fPRCorr50EtaH,"part #eta","Rec Jet #eta");
+ fPRCorr50PhiH = new TH2F("PRCorr50Phi","Azimuthal angle correlation part-rec jet",
+ 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
+ SetProperties(fPRCorr50PhiH,"Part #phi","Rec Jet #phi");
+ }
+
+ if (fDoGenJ && fDoRecJ) {
+ fRGCorr50EneH = new TH2F("RGCorr50Ene","Energy correlation rec jet-gen jet",
+ 40,0,200,40,0,200);
+ SetProperties(fRGCorr50EneH,"Rec Jet Energy (GeV)","Gen Jet Energy (GeV)");
+ fRGCorr50PtH = new TH2F("RGCorr50Pt","Pt correlation rec jet-gen jet",
+ 40,0,200,40,0,200);
+ SetProperties(fRGCorr50PtH,"Rec Jet P_{T} (GeV)","Gen Jet P_{T} (GeV)");
+ fRGCorr50EtaH = new TH2F("RGCorr50Eta","Pseudorapidity correlation rec jet-gen jet",
+ 40,-1.0,1.0,40,-1.0,1.0);
+ SetProperties(fRGCorr50EtaH,"Rec Jet #eta","Gen Jet #eta");
+ fRGCorr50PhiH = new TH2F("RGCorr50Phi","Azimuthal angle correlation rec jet-gen jet",
+ 90,0.,2.0*TMath::Pi(),90,0.,2.0*TMath::Pi());
+ SetProperties(fRGCorr50PhiH,"Rec Jet #phi","Gen Jet #phi");
+ }
+}
+
+void AliJetAnalysis::DefineShapH()
+{
+ // leading reconstructed jet
+ if (fDoRecJ) {
+ fdEdrH = new TH2F("fdEdrH","dE/dr histo",20,0,1,40,0,200);
+ SetProperties(fdEdrH,"r","Rec Jet P_{T}");
+ fdEdrB = new TH2F("fdEdrB","dE/dr bkgdhisto",20,0,1,40,0,200);
+ SetProperties(fdEdrB,"r","Rec P_{T}");
+ fPtEneH2 = new TH2F("fPtEneH2","fPtEneH2",40,0,200,40,0,200);
+ SetProperties(fPtEneH2,"Rec Jet E","Rec Jet P_{T}");
+ fdEdrW = new TH1F("fdEdrW","weights for dE/dr",40,0,200);
+ SetProperties(fdEdrW,"Rec Jet P_{T}","weights");
- TH1F* e0H = new TH1F("e0H" ,"Jet Energy (reconstructed)", 40, 0., 200.);
- TH1F* e1H = new TH1F("e1H" , "Jet Energy (generated)", 40, 0., 200.);
- TH1F* e2H = new TH1F("e2H" , "Jet Energy (generated, nrec = 0", 40, 0., 200.);
- TH1F* e3H = new TH1F("e3H" , "Jet Energy (leading)", 40, 0., 200.);
- TH1F* e4H = new TH1F("e4H" , "Jet Energy (reconstructed: 105 < Egen < 125", 40, 0., 200.);
+ fRShapSelH = new TH1F("RShapSel","Shape of generated jets (sel part)",20,0.,1.);
+ SetProperties(fRShapSelH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
+ fRShapRejH = new TH1F("RShapRej","Shape of generated jets (rej part)",20,0.,1.);
+ SetProperties(fRShapRejH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
+ fRShapAllH = new TH1F("RShapAll","Shape of generated jets (all part)",20,0.,1.);
+ SetProperties(fRShapAllH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R_{JET})");
+ }
+}
- TH1F* e5H = new TH1F("e5H" , "Jet Energy (generated)", 40, 0., 200.);
- TH1F* e6H = new TH1F("e6H" , "Jet Energy (generated)", 40, 0., 200.);
- TH1F* e7H = new TH1F("e7H" , "Jet Energy (generated)", 40, 0., 200.);
- TH1F* e8H = new TH1F("e8H" , "Jet Energy (generated)", 40, 0., 200.);
+void AliJetAnalysis::DefineJtH()
+{
+ // Define the histogram for J_T
+ if (fDoRecJ) {
+ fJtH = new TH2F("fJtH","J_{T} histogram",80,0.,4.,40,0.,200.);
+ SetProperties(fJtH,"J_{T}","Rec Jet P_{T}");
+ fJtB = new TH2F("fJtB","J_{T} bkgd histogram",80,0.,4.,40,0.,200.);
+ SetProperties(fJtB,"J_{T}","Rec P_{T}");
+ fJtW = new TH1F("fJtW","J_{T} weight",40,0,200);
+ SetProperties(fJtW,"J_{T}W","weight");
+ }
+}
- TProfile* r5H = new TProfile("r5H" , "rec/generated", 20, 0., 200, 0., 1., "S");
- TProfile* r6H = new TProfile("r6H" , "rec/generated", 20, 0., 200, 0., 1., "S");
+void AliJetAnalysis::DefinedNdxiH()
+{
+ // Define the histogram for dN/dxi
+ if (fDoRecJ) {
+ fdNdxiH = new TH2F("fdNdxiH","dN/d#xi histo",200,0,10,40,0,200);
+ SetProperties(fdNdxiH,"xi","Rec Jet P_{T}");
+ fdNdxiB = new TH2F("fdNdxiB","dN/d#xi bkgd histo",200,0,10,40,0,200);
+ SetProperties(fdNdxiB,"xi","Rec P_{T}");
+ fdNdxiW = new TH1F("fdNdxiW","dN/d#xi histo",40,0,200);
+ SetProperties(fdNdxiW,"Rec Jet P_{T}","weights");
+ fPtEneH = new TH2F("fPtEneH","fPtEneH",40,0,200,40,0,200);
+ SetProperties(fPtEneH,"Rec Jet E","Rec Jet P_{T}");
+ }
+}
- TProfile* r7H = new TProfile("r7H" , "rec/generated", 20, 0., 200, 0., 1., "S");
- TProfile* r8H = new TProfile("r8H" , "rec/generated", 20, 0., 200, 0., 1., "S");
+void AliJetAnalysis::DefineFragH()
+{
+ // leading reconstructed jet
+ if (fDoRecJ) {
+ fRFragSelH = new TH1F("RFragSel","Frag Fun of reconstructed jets (sel part)",20,0.,10.);
+ SetProperties(fRFragSelH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
+ fRFragRejH = new TH1F("RFragRej","Frag Fun of reconstructed jets (rej part)",20,0.,10.);
+ SetProperties(fRFragRejH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
+ fRFragAllH = new TH1F("RFragAll","Frag Fun of reconstructed jets (all part)",20,0.,10.);
+ SetProperties(fRFragAllH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
+ }
+}
+void AliJetAnalysis::DefineTrigH()
+{
+ // generated energy
+ fGTriggerEneH = new TProfile("GTriggerEne","Generated energy (trigger bias)",
+ 20, 0., 200., 0., 1., "S");
+ fGTriggerEneH->SetXTitle("E_{gen}");
+ fGTriggerEneH->SetYTitle("E_{rec}/E_{gen}");
+ // reconstructed energy
+ fRTriggerEneH = new TProfile("RTriggerEne","Reconstructed energy (trigger bias)",
+ 20, 0., 200., 0., 1., "S");
+ fRTriggerEneH->SetXTitle("E_{rec}");
+ fRTriggerEneH->SetYTitle("E_{rec}/E_{gen}");
+ // generated energy vs generated/leading
+ fGPTriggerEneH = new TProfile("GPTriggerEne","Generated energy (trigger bias)",
+ 20, 0., 200., 0., 1., "S");
+ fGPTriggerEneH->SetXTitle("E_{gen}");
+ fGPTriggerEneH->SetYTitle("E_{L}/E_{gen}");
+ // leading particle energy
+ fPTriggerEneH = new TProfile("PTriggerEne","Leading particle energy (trigger bias)",
+ 20, 0., 200., 0., 1., "S");
+ fPTriggerEneH->SetXTitle("E_{L}/0.2");
+ fPTriggerEneH->SetYTitle("E_{L}/E_{gen}");
- TH1F* dr1H = new TH1F("dr1H", "delta R", 160, 0., 2.);
- TH1F* dr2H = new TH1F("dr2H", "delta R", 160, 0., 2.);
- TH1F* dr3H = new TH1F("dr4H", "delta R", 160, 0., 2.);
+}
- TH1F* etaH = new TH1F("etaH", "eta", 160, -2., 2.);
- TH1F* eta1H = new TH1F("eta1H", "eta", 160, -2., 2.);
- TH1F* eta2H = new TH1F("eta2H", "eta", 160, -2., 2.);
+void AliJetAnalysis::SetProperties(TH1* h,const char* x, const char* y) const
+{
+ //Set properties of histos (x and y title and error propagation)
+ h->SetXTitle(x);
+ h->SetYTitle(y);
+ h->Sumw2();
+}
+
+////////////////////////////////////////////////////////////////////////
+// compute weights for dE/dr histogram
- TH1F* phiH = new TH1F("phiH", "phi", 160, -3., 3.);
- TH1F* dphiH = new TH1F("dphiH", "phi", 160, 0., 3.1415);
- TH1F* phi1H = new TH1F("phi1H", "phi", 160, 0., 6.28);
- TH1F* phi2H = new TH1F("phi2H", "phi", 160, 0., 6.28);
-
+void AliJetAnalysis::SetdEdrWeight()
+{
+ // Due to the limited acceptance, each bin in the dE/dr has a different
+ // weight given by the ratio of the area of a disk dR to the area
+ // within the eta acceptance. The weight depends on the eta position
+ // of the jet. Here a look up table for the weights is computed. It
+ // assumes |etaJet|<0.5 and |eta_lego|<0.9. It makes bins of 0.05
+ // units in eta. Note that this is fixed by the bin width chosen for
+ // the histogram --> this weights are tailored for the specific
+ // histogram definition used here!
+
+ // two dimensional table: first index, bin in eta of jet, second
+ // index bin of dE/dr histo
+
+ Int_t nBins = 20;
+ Float_t xMin = 0.0;
+ Float_t xMax = 1.0;
+ Float_t binSize = (xMax-xMin)/nBins;
- TProfile* drP1 = new TProfile("drP1" , "Delta_R", 20, 0., 200, -1., 1., "S");
- TProfile* drP2 = new TProfile("drP2" , "Delta_R", 20, 0., 200, -1., 1., "S");
+ Float_t da,ds,r1,r2,h1,h2,a1,a2,theta1,theta2,area1,area2;
+ Int_t ji;
+ for (Int_t i=0;i<(nBins/2);i++) {
+ // index of first histo bin needing a scale factor
+ ji=(nBins-2)-i;
+ for(Int_t j=0;j<nBins;j++) {
+ // area of ring.
+ da = TMath::Pi()*(binSize*binSize)*(1.0+2.0*j);
+ ds = 1.0;
+ if (j>=ji) { // compute missing area using segments of circle
+ r1=j*binSize;
+ r2=(j+1)*binSize;
+ h1=(j-ji)*binSize;
+ h2=(j+1-ji)*binSize;
+ a1=2.0*TMath::Sqrt(2.0*h1*r1-h1*h1);
+ a2=2.0*TMath::Sqrt(2.0*h2*r2-h2*h2);
+ theta1=2*TMath::ACos((r1-h1)/r1);
+ theta2=2*TMath::ACos((r2-h2)/r2);
+ area1=binSize*(r1*r1*theta1-a1*(r1-h1));
+ area2=binSize*(r2*r2*theta2-a2*(r2-h2));
+ ds = (da-(area2-area1))/da;
+ }
+ fWeightdEdr[i][j]=ds/da;
+ }
+ }
+}
+
+// get weight for dE/dr histogram
+Float_t AliJetAnalysis::GetdEdrWeight(Float_t etaJet, Float_t r)
+{
+ // Return the correponding weight for the dE/dr plot
+ Int_t nBins = 20;
+ Float_t xMin = 0.0;
+ Float_t xMax = 1.0;
+ Float_t binSize = (xMax-xMin)/nBins;
+
+ Float_t eta = TMath::Abs(etaJet);
+ if ((eta > 0.5) || (r > fdrdEdr)) return 0.0;
+ Int_t binJet = (Int_t) TMath::Floor(eta/binSize);
+ Int_t binR = (Int_t) TMath::Floor(r/binSize);
+ Float_t w = fWeightdEdr[binJet][binR];
+ return w;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+// fill histograms
+
+void AliJetAnalysis::FillHistograms()
+{
+ // fill histograms
// Run data
- AliJetProductionDataPDC2004* runData = new AliJetProductionDataPDC2004();
+ AliJetProductionDataPDC2004* runData = new AliJetProductionDataPDC2004();
+
+ // Loop over runs
+ TFile* jFile = 0x0;
+
+ for (Int_t iRun = fRunMin; iRun <= fRunMax; iRun++) {
+ // Open file
+ char fn[256];
+ sprintf(fn, "%s/%s.root", fDirectory, (runData->GetRunTitle(iRun)).Data());
+ jFile = new TFile(fn);
+ printf(" Analyzing run: %d %s\n", iRun,fn);
+
+ // Get reader header and events to be looped over
+ AliJetReaderHeader *jReaderH = (AliJetReaderHeader*)(jFile->Get(fReaderHeader));
+ if (fEventMin == -1) fEventMin = jReaderH->GetFirstEvent();
+ if (fEventMax == -1) {
+ fEventMax = jReaderH->GetLastEvent();
+ } else {
+ fEventMax = TMath::Min(fEventMax, jReaderH->GetLastEvent());
+ }
+
+ AliUA1JetHeader *jH = (AliUA1JetHeader *) (jFile->Get("AliUA1JetHeader"));
+ // Calculate weight
+ fWeight = runData->GetWeight(iRun)/ ( (Float_t) (fEventMax - fEventMin + 1));
- // Loop over runs
- TFile* jFile = 0x0;
-
- for (Int_t iRun = fRunMin; iRun <= fRunMax; iRun++)
- {
-
- // Open file
- char fn[20];
- sprintf(fn, "%s/%s.root", fDirectory, (runData->GetRunTitle(iRun)).Data());
-
-
- jFile = new TFile(fn);
-
- printf(" Analyzing run: %d %s\n", iRun,fn);
- // Get jet header and display parameters
- AliUA1JetHeader* jHeader =
- (AliUA1JetHeader*) (jFile->Get("AliUA1JetHeader"));
- // jHeader->PrintParameters();
-
- // Get reader header and events to be looped over
- AliJetESDReaderHeader *jReaderH =
- (AliJetESDReaderHeader*)(jFile->Get("AliJetKineReaderHeader"));
-
- if (fEventMin == -1) fEventMin = jReaderH->GetFirstEvent();
- if (fEventMax == -1) {
- fEventMax = jReaderH->GetLastEvent();
- } else {
- fEventMax = TMath::Min(fEventMax, jReaderH->GetLastEvent());
+ // Loop over events
+ for (Int_t i = fEventMin; i < fEventMax; i++) {
+ if (i%100 == 0) printf(" Analyzing event %d / %d \n",i,fEventMax);
+
+ // Get next tree with AliJet
+ char nameT[100];
+ sprintf(nameT, "TreeJ%d",i);
+ TTree *jetT =(TTree *)(jFile->Get(nameT));
+ if (fDoRecJ) jetT->SetBranchAddress("FoundJet", &fRecJ);
+ if (fDoGenJ) jetT->SetBranchAddress("GenJet", &fGenJ);
+ if (fDoPart) jetT->SetBranchAddress("LeadingPart", &fPart);
+ jetT->GetEntry(0);
+
+ int nJets = fRecJ->GetNJets();
+
+ TArrayI inJet = fRecJ->GetInJet();
+ if(inJet.GetSize()>fminMult){ // removes events with low multiplicity
+ if (fDoKine) FillKineH();
+ if (fDoCorr) FillCorrH();
+ if (fDoCorr50) FillCorr50H();
+ if (fDoShap) FillShapH(jH->GetRadius());
+ if (fDoFrag) FillFragH();
+ if (fDoTrig) FillTrigH();
+ if (fDoJt) FillJtH();
+ if (fDodNdxi) FilldNdxiH();
+ }
+ delete jetT; // jet should be deleted before creating a new one
+ if(inJet.GetSize()>fminMult){ // removes events with low multiplicity
+ if (fDoBkgd && nJets>0) FillBkgd(i,iRun);
+ }
+ } // end loop over events in one file
+ if (jFile) jFile->Close();
+ delete jFile;
+ } // end loop over files
+}
+
+void AliJetAnalysis::FillKineH()
+{
+ // leading particle
+ if (fDoPart && fPart->LeadingFound()){
+ fPKineEneH->Fill(fPart->GetE(),fWeight);
+ fPKinePtH->Fill(fPart->GetPt(),fWeight);
+ fPKinePhiH->Fill(fPart->GetPhi(),fWeight);
+ fPKineEtaH->Fill(fPart->GetEta(),fWeight);
+ }
+ // leading generated jet
+ if (fDoGenJ && fGenJ->GetNJets()> 0){
+ fGKineEneH->Fill(fGenJ->GetE(0),fWeight);
+ fGKinePtH->Fill(fGenJ->GetPt(0),fWeight);
+ fGKinePhiH->Fill(fGenJ->GetPhi(0),fWeight);
+ fGKineEtaH->Fill(fGenJ->GetEta(0),fWeight);
+ }
+ // leading reconstructed jet
+ if (fDoRecJ && fRecJ->GetNJets()> 0) {
+ TArrayF p=fRecJ->GetPtFromSignal();
+ if (p[0]>fPercentage) {
+ fRKineEneH->Fill(fRecJ->GetE(0)/fEfactor,fWeight);
+ fRKinePtH->Fill(fRecJ->GetPt(0),fWeight);
+ fRKinePhiH->Fill(fRecJ->GetPhi(0),fWeight);
+ fRKineEtaH->Fill(fRecJ->GetEta(0),fWeight);
+ }
+ }
+}
+
+void AliJetAnalysis::FillCorrH()
+{
+ // Fill correlation histograms
+ if (fDoPart && fPart->LeadingFound() && fDoGenJ && fGenJ->GetNJets()> 0)
+ Correlation(fPart->GetLeading(),fGenJ->GetJet(0),
+ fPGCorrEneH,fPGCorrPtH,fPGCorrEtaH,fPGCorrPhiH);
+ if (fDoPart && fPart->LeadingFound() && fDoRecJ && fRecJ->GetNJets()> 0) {
+ TArrayF p=fRecJ->GetPtFromSignal();
+ if (p[0]>fPercentage)
+ Correlation(fPart->GetLeading(),fRecJ->GetJet(0),
+ fPRCorrEneH,fPRCorrPtH,fPRCorrEtaH,fPRCorrPhiH);
+ }
+ if (fDoGenJ && fGenJ->GetNJets()> 0 && fDoRecJ && fRecJ->GetNJets()> 0) {
+ TArrayF p=fRecJ->GetPtFromSignal();
+ if (p[0]>fPercentage)
+ Correlation(fRecJ->GetJet(0), fGenJ->GetJet(0),
+ fRGCorrEneH,fRGCorrPtH,fRGCorrEtaH,fRGCorrPhiH);
+ }
+}
+
+void AliJetAnalysis::Correlation(TLorentzVector *lv1,TLorentzVector *lv2,
+ TH2F *h1, TH2F *h2, TH2F *h3, TH2F *h4)
+{
+ // Correlation histograms
+ h1->Fill(lv1->E(),lv2->E(),fWeight);
+ h2->Fill(lv1->Pt(),lv2->Pt(),fWeight);
+ h3->Fill(lv1->Eta(),lv2->Eta(),fWeight);
+ Float_t p1, p2;
+ p1 = ((lv1->Phi() < 0) ? (lv1->Phi()) + 2. * TMath::Pi() : lv1->Phi());
+ p2 = ((lv2->Phi() < 0) ? (lv2->Phi()) + 2. * TMath::Pi() : lv2->Phi());
+ h4->Fill(p1,p2,fWeight);
+}
+
+void AliJetAnalysis::FillCorr50H()
+{
+ // Fill correlation histograms when one particle has > 50% of jet energy
+ if (fDoRecJ && fRecJ->GetNJets()> 0) {
+ TArrayF p = fRecJ->GetPtFromSignal();
+ if (p[0]>fPercentage) {
+ if (fDoPart && fPart->LeadingFound())
+ Correlation50(fRecJ, fPart->GetLeading(),fRecJ->GetJet(0),
+ fPRCorr50EneH,fPRCorr50PtH,fPRCorr50EtaH,fPRCorr50PhiH);
+ if (fDoGenJ && fGenJ->GetNJets()> 0)
+ Correlation50(fRecJ, fRecJ->GetJet(0), fGenJ->GetJet(0),
+ fRGCorr50EneH,fRGCorr50PtH,fRGCorr50EtaH,fRGCorr50PhiH);
+ }
+ }
+}
+
+void AliJetAnalysis::Correlation50(AliJet *j,TLorentzVector *lv1,TLorentzVector *lv2,
+ TH2F *h1, TH2F *h2, TH2F *h3, TH2F *h4)
+{
+ // Correlation histograms when one particle has > 50% of jet energy
+ TArrayF ptin = j->GetPtIn();
+ TArrayF etain = j->GetEtaIn();
+ TArrayI inJet = j->GetInJet();
+
+ Int_t flag = 0;
+ for(Int_t i=0;i<(inJet.GetSize());i++) {
+ if (inJet[i] == 1) {
+ Float_t t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
+ Float_t x = (ptin[i]*TMath::Sqrt(1.+1./(t1*t1)))/(j->GetE(0));
+ if (x>0.5) flag++;
+ }
+ }
+ if (flag>1) cout << " Flag = " << flag << endl;
+ if (flag == 1) {
+ h1->Fill(lv1->E(),lv2->E(),fWeight);
+ h2->Fill(lv1->Pt(),lv2->Pt(),fWeight);
+ h3->Fill(lv1->Eta(),lv2->Eta(),fWeight);
+ Float_t p1, p2;
+ p1 = ((lv1->Phi() < 0) ?
+ (lv1->Phi()) + 2. * TMath::Pi() : lv1->Phi());
+ p2 = ((lv2->Phi() < 0) ?
+ (lv2->Phi()) + 2. * TMath::Pi() : lv2->Phi());
+ h4->Fill(p1,p2,fWeight);
+ }
+}
+
+void AliJetAnalysis::FillJtH()
+{
+ // Fill J_T histogram
+ if (fRecJ->GetNJets()> 0) {
+ fjv3X = 0.0; fjv3Y = 0.0; fjv3Z = 0.0;
+ TArrayF p=fRecJ->GetPtFromSignal();
+ if (p[0]>fPercentage) {
+ // initialize
+ const TVector3 kJv3 = fRecJ->GetJet(0)->Vect();
+ fjv3X = kJv3.X(); fjv3Y = kJv3.Y(); fjv3Z = kJv3.Z();
+ TVector3 trk;
+ TArrayI inJet = fRecJ->GetInJet();
+ TArrayF etain = fRecJ->GetEtaIn();
+ TArrayF ptin = fRecJ->GetPtIn();
+ TArrayF phiin = fRecJ->GetPhiIn();
+ Float_t deta, dphi,jt, dr;
+ for(Int_t i=0;i<inJet.GetSize();i++) {
+ deta = etain[i] - fRecJ->GetEta(0);
+ dphi = phiin[i] - fRecJ->GetPhi(0);
+ if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+ if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+ dr = TMath::Sqrt(deta * deta + dphi * dphi);
+ if ((dr<fdrJt) && (ptin[i] > fPartPtCut)) {
+ trk.SetPtEtaPhi(ptin[i],etain[i],phiin[i]);
+ jt = TMath::Sin(trk.Angle(kJv3))*trk.Mag();
+ fJtH->Fill(jt,fRecJ->GetPt(0),fWeight);
}
-
-
- // Calculate weight
- Float_t wgt = runData->GetWeight(iRun) / Float_t(fEventMax - fEventMin + 1);
- Float_t ptmin, ptmax;
- runData->GetPtHardLimits(iRun, ptmin, ptmax);
-
-
- // Loop over events
- AliJet *jets = 0x0;
- AliJet *gjets = 0x0;
- AliLeading *leading = 0x0;
- Float_t egen, etag, econ, erec;
-
-
- for (Int_t i = fEventMin; i < fEventMax; i++) {
- printf(" Analyzing run: %d Event %d / %d \n",
- iRun, i, fEventMax);
-
- // Het next tree with AliJet
- char nameT[100];
- sprintf(nameT, "TreeJ%d",i);
- TTree *jetT =(TTree *)(jFile->Get(nameT));
- jetT->SetBranchAddress("FoundJet", &jets);
- jetT->SetBranchAddress("GenJet", &gjets);
- jetT->SetBranchAddress("LeadingPart", &leading);
- jetT->GetEntry(0);
-
-//
-// Find the jet with the highest E_T within fiducial region
-//
- Int_t njets = jets->GetNJets();
- Int_t imax = -1;
- Int_t jmax = -1;
- Float_t emax = 0.;
-
- for (Int_t ij = 0; ij < njets; ij++) {
- if (jets->GetPt(ij) > emax &&
- TMath::Abs(jets->GetEta(ij)) < 0.60) {
- emax = jets->GetPt(ij);
- jmax = imax;
- imax = ij;
- }
- }
-
-
- if (imax == -1) {
- Int_t ngen = gjets->GetNJets();
- if(ngen>0) e2H->Fill(gjets->GetPt(0), wgt);
- } else {
-// printf("Reconstructed Jet %5d %13.3f %13.3f %13.3f\n",
-// imax, emax, jets->GetEta(imax), jets->GetPhi(imax));
- econ = jets->GetPt(imax);
- erec = jets->GetPt(imax) / 0.65;
- dr1H->Fill(jets->GetEta(imax) - gjets->GetEta(0));
-//
-// Find the generated jet closest to the reconstructed
-//
- Float_t rmin;
- Float_t etamin = 1.e6;
-
- Int_t igen;
- Float_t etaj = jets->GetEta(imax);
- Float_t phij = jets->GetPhi(imax);
-
- Int_t ngen = gjets->GetNJets();
-
- if (ngen != 0) {
- rmin = 1.e6;
- igen = -1;
- for (Int_t j = 0; j < ngen; j++) {
- etag = gjets->GetEta(j);
- Float_t phig = gjets->GetPhi(j);
- Float_t deta = etag - etaj;
- Float_t dphi = TMath::Abs(phig - phij);
- if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
- Float_t r = TMath::Sqrt(deta * deta + dphi * dphi);
- if (r < rmin) {
- rmin = r;
- etamin = deta;
- igen = j;
- }
- }
-
- egen = gjets->GetPt(igen);
- e1H->Fill(gjets->GetPt(igen), wgt);
- etag = gjets->GetEta(igen);
- Float_t phig = gjets->GetPhi(igen);
- Float_t dphi = phig - phij;
-
-// if (econ < ptmax) {
- e0H->Fill(erec, wgt);
-// } else {
-// e0H->Fill(erec, 6.7e-6);
-// }
-
-
-
- if (egen > 20. && egen < 40.) {
- phiH->Fill((dphi));
- etaH->Fill(etag - etaj, wgt);
- phi1H->Fill(phig);
- phi2H->Fill(phij);
- e4H->Fill(jets->GetPt(imax));
- }
-
- if (erec > 90. && erec < 110. && rmin < 0.1) {
- e5H->Fill(egen, wgt);
- dr2H->Fill(rmin);
- if (egen < 30.) {
- printf("Strange jet %6d %13.3f %13.3f %13.3f \n",
- imax, etaj, phij, erec);
- for (Int_t j = 0; j < ngen; j++) {
- printf("Generated %6d %13.3f %13.3f %13.3f\n",
- j, gjets->GetEta(j),
- gjets->GetPhi(j), gjets->GetPt(j));
-
- }
- }
- }
-
- if (rmin < 0.1) {
- r5H->Fill(egen, jets->GetPt(imax)/egen, wgt);
- r6H->Fill(jets->GetPt(imax)
- / 0.4, jets->GetPt(imax)/egen, wgt);
- e8H->Fill(erec, wgt);
- }
-
- if (rmin < 0.1) {
- drP1->Fill(egen, etamin, wgt);
- }
- } // ngen !=0
- } // has reconstructed jet
-
-//
-// Leading particle
-//
- if (leading->LeadingFound()) {
- Float_t etal = leading->GetLeading()->Eta();
- Float_t phil = leading->GetLeading()->Phi();
- Float_t el = leading->GetLeading()->E();
-// printf("Leading %f %f %f \n", etal, phil, el);
-
- e3H->Fill(el, wgt);
-
- Float_t rmin = 1.e6;
- Float_t etamin = 1.e6;
-
- Int_t igen = -1;
- Int_t ngen = gjets->GetNJets();
- for (Int_t j = 0; j < ngen; j++) {
- etag = gjets->GetEta(j);
- Float_t phig = gjets->GetPhi(j);
- Float_t deta = etag-etal;
- Float_t dphi = TMath::Abs(phig - phil);
- if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
-
- Float_t r = TMath::Sqrt(deta * deta + dphi * dphi);
-
- if (r < rmin) {
- rmin = r;
- igen = j;
- etamin = deta;
- }
- }
- if (egen > 20. && egen < 40.) {
- dr3H->Fill(rmin);
- eta1H->Fill(etag-etal, wgt);
-
- }
-
- if (el > 54. && el < 66.) e6H->Fill(egen, wgt);
- if (rmin < 0.3) {
- r7H->Fill(egen, el/egen, wgt);
- r8H->Fill(el / 0.2, el/egen, wgt);
- }
-
- if (rmin < 0.1) {
- drP2->Fill(egen, etamin, wgt);
- }
- } // if leading particle
-
-
-
-//
-// Generated Jet
-//
- Int_t ngen = gjets->GetNJets();
- emax = 0.;
- imax = -1;
-
- for (Int_t j = 0; j < ngen; j++) {
- if (gjets->GetPt(j) >
- emax && TMath::Abs(gjets->GetEta(j)) < 0.5) {
- emax = gjets->GetPt(j);
- imax = j;
- }
- }
- if (imax != -1) e7H->Fill(emax, wgt);
-
-
- delete jetT;
-
- } // events
- if (jFile) jFile->Close();
- delete jFile;
-
- } // runs
+ }
+ fJtW->Fill(fRecJ->GetPt(0),fWeight);
+ fWJt+=fWeight;
+ }
+ }
+}
+
+void AliJetAnalysis::FilldNdxiH()
+{
+ // Fill dN/d#xi histograms
+ if (fRecJ->GetNJets()> 0) {
+ TArrayF p=fRecJ->GetPtFromSignal();
+ if (p[0]>fPercentage) {
+ fp0 = p[0];
+ TArrayI inJet = fRecJ->GetInJet();
+ TArrayF etain = fRecJ->GetEtaIn();
+ TArrayF ptin = fRecJ->GetPtIn();
+ TArrayF phiin = fRecJ->GetPhiIn();
+ Float_t xi,t1,ene,dr,deta,dphi;
+ for(Int_t i=0;i<inJet.GetSize();i++) {
+ deta = etain[i] - fRecJ->GetEta(0);
+ dphi = phiin[i] - fRecJ->GetPhi(0);
+ if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+ if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+ dr = TMath::Sqrt(deta * deta + dphi * dphi);
+ if ((dr<fdrdNdxi) && (ptin[i] > fPartPtCut)) {
+ t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
+ ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
+ xi = (Float_t) TMath::Log((fRecJ->GetE(0)/fEfactor)/ene);
+ fdNdxiH->Fill(xi,fRecJ->GetPt(0),fWeight);
+ }
+ }
+ fdNdxiW->Fill(fRecJ->GetPt(0),fWeight);
+ fPtEneH->Fill(fRecJ->GetE(0)/fEfactor,fRecJ->GetPt(0),fWeight);
+ fWdNdxi+=fWeight;
+ }
+ }
+}
+
+void AliJetAnalysis::FillBkgd(Int_t eventN, Int_t runN)
+{
+ // Background calculated with hijing events (no pythia)
+ if (fp0>fPercentage) {
+ fPtJ=0.,fEJ=0.,fEtaJ=0.,fPhiJ=0.;
+ // Background calculated with hijing events (no pythia)
+ AliJetProductionDataPDC2004* runDataB = new AliJetProductionDataPDC2004();
+ TFile* jFileB =0x0;;
+ char fnB[256];
-// delete jFile;
-// if (jFile) jFile->Close();
-/*
- TFile* f = new TFile("j.root", "recreate");
- e0H->Write();
- e1H->Write();
- e2H->Write();
- e3H->Write();
- e4H->Write();
- e7H->Write();
- e8H->Write();
- f->Close();
-*/
-
- // Get Control Plots
-// gStyle->SetOptStat(0);
+ sprintf(fnB, "%s/%s.root", fBkgdDirectory, (runDataB->GetRunTitle(runN)).Data());
+ jFileB = new TFile(fnB);
+
+ char nameB[100];
+ sprintf(nameB, "TreeJ%d",eventN);
+ TTree *jetB =(TTree *)(jFileB->Get(nameB));
+ jetB->SetBranchAddress("FoundJet", &fRecB);
+ jetB->GetEntry(0);
- TCanvas* c1 = new TCanvas("c1");
- e0H->Draw();
- e1H->SetLineColor(2);
- e2H->SetLineColor(4);
- e3H->SetLineColor(5);
- e1H->Draw("same");
- e3H->Draw("same");
-
-
- TCanvas* c2 = new TCanvas("c2");
-// dr1H->Draw();
- dr2H->SetLineColor(2);
- dr2H->Draw("");
-
- TCanvas* c3 = new TCanvas("c3");
- dr2H->Draw();
- dr3H->Draw("same");
+ TArrayI inJetB = fRecB->GetInJet();
+ TArrayF etainB = fRecB->GetEtaIn();
+ TArrayF ptinB = fRecB->GetPtIn();
+ TArrayF phiinB = fRecB->GetPhiIn();
+ fPtJ = fRecJ->GetPt(0);
+ fEJ = fRecJ->GetE(0);
+ fEtaJ = fRecJ->GetEta(0);
+ fPhiJ = fRecJ->GetPhi(0);
+ Float_t t1,ene,xi,detaB,dphiB,drB,jt,wB;
+ TVector3 trkB;
+ TVector3 jv3B;
+ jv3B.SetX(fjv3X); jv3B.SetY(fjv3Y); jv3B.SetZ(fjv3Z);
+
+ for(Int_t k=0;k<inJetB.GetSize();k++){
+ if(ptinB[k] > fPartPtCut){
+ detaB = etainB[k] - fEtaJ;
+ dphiB = phiinB[k] - fPhiJ;
+ if (dphiB < -TMath::Pi()) dphiB= -dphiB - 2.0 * TMath::Pi();
+ if (dphiB > TMath::Pi()) dphiB = 2.0 * TMath::Pi() - dphiB;
+ drB = TMath::Sqrt(detaB * detaB + dphiB * dphiB);
+ t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etainB[k])));
+ ene = ptinB[k]*TMath::Sqrt(1.+1./(t1*t1));
+ trkB.SetPtEtaPhi(ptinB[k],etainB[k],phiinB[k]);
+ // --- dN/dxi
+ if (drB<fdrdNdxi) {
+ xi = (Float_t) TMath::Log(fEJ/ene);
+ fdNdxiB->Fill(xi,fPtJ,fWeight);
+ }
+ // --- Jt
+ if (drB<fdrJt) {
+ jt = TMath::Sin(trkB.Angle(jv3B))*(trkB.Mag());
+ fJtB->Fill(jt,fPtJ,fWeight);
+ }
+ // --- dE/dr
+ if (drB<fdrdEdr) {
+ wB = GetdEdrWeight(fEtaJ,drB)*fWeight*ene;
+ fdEdrB->Fill(drB,fPtJ,wB);
+ }
+ }
+ }
+ delete jetB;
+ if (jFileB) jFileB->Close();
+ delete jFileB;
+ }
+}
- TCanvas* c4 = new TCanvas("c4");
- e0H->Draw();
+void AliJetAnalysis::FillShapH(Float_t r)
+{
+ // Fill jet shape histograms
+ if (fDoRecJ && fRecJ->GetNJets()> 0) {
+ TArrayF p=fRecJ->GetPtFromSignal();
+ if (p[0]>fPercentage) {
+ Shape(fRecJ,fRShapSelH,fRShapRejH,fRShapAllH,fdEdrH,fPtEneH2,fdEdrW,r);
+ fWdEdr+=fWeight;
+ fWShapR+=fWeight;
+ }
+ }
+}
- TCanvas* c5 = new TCanvas("c5");
- etaH->SetXTitle("#eta_{rec} - #eta_{gen}");
-
- etaH->Draw();
- eta1H->SetLineColor(2);
+void AliJetAnalysis::Shape(AliJet *j,TH1F* hs, TH1F* hr, TH1F* ha,
+ TH2F *hdedr, TH2F *hptene, TH1F* wdedr, Float_t radius)
+{
+ // initialize
+ TArrayI inJet = j->GetInJet();
+ TArrayF etain = j->GetEtaIn();
+ TArrayF ptin = j->GetPtIn();
+ TArrayF phiin = j->GetPhiIn();
- eta1H->Draw("same");
+ // first compute dE/dr histo
+ Float_t etaj = j->GetEta(0);
+ Float_t ene,w,deta,dphi,dr,t1;
+ for(Int_t i=0;i<inJet.GetSize();i++) {
+ deta = etain[i] - j->GetEta(0);
+ dphi = phiin[i] - j->GetPhi(0);
+ if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+ if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+ dr = TMath::Sqrt(deta * deta + dphi * dphi);
+ if ((dr<fdrdEdr) && (ptin[i] > fPartPtCut)) {
+ t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
+ ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
+ w = GetdEdrWeight(etaj,dr)*fWeight*ene;
+ hdedr->Fill(dr,j->GetPt(0),w);
+ }
+ }
+ hptene->Fill(fRecJ->GetE(0),fRecJ->GetPt(0),fWeight);
+ wdedr->Fill(j->GetPt(0),fWeight);
- TCanvas* c5a = new TCanvas("c5a");
- eta1H->Draw();
-
- TCanvas* c5b = new TCanvas("c5b");
- eta2H->Draw();
-
- TCanvas* c6 = new TCanvas("c6");
- e4H->Draw();
- TCanvas* c7 = new TCanvas("c7");
- phiH->Draw();
-
- TCanvas* c7a = new TCanvas("c7a");
- phi1H->Draw();
- TCanvas* c7b = new TCanvas("c7b");
- phi2H->Draw();
-
- TCanvas* c8 = new TCanvas("c8");
- e5H->SetXTitle("E_{gen} (GeV)");
- e5H->Draw();
- e6H->SetLineColor(2);
- e6H->Draw("same");
+ // now compute shape histos
+ Int_t nBins = ha->GetNbinsX();
+ Float_t xMin = ha->GetXaxis()->GetXmin();
+ Float_t xMax = ha->GetXaxis()->GetXmax();
+ Float_t binSize,halfBin;
+ binSize = (xMax-xMin)/nBins;
+ halfBin = binSize/2;
+ Float_t rptS[20], rptR[20], rptA[20];
+ for(Int_t i=0;i<nBins;i++) rptS[i]=rptR[i]=rptA[i]=0.0;
+ // fill bins in r for leading jet
+ for(Int_t i=0;i<inJet.GetSize();i++) {
+ if (inJet[i] == 1 || (inJet[i] == -1 && fPythia)) {
+ deta = etain[i] - j->GetEta(0);
+ dphi = phiin[i] - j->GetPhi(0);
+ if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+ if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+ dr = TMath::Sqrt(deta * deta + dphi * dphi);
+ Float_t rR = dr/radius;
+ Int_t bin = (Int_t) TMath::Floor(rR/binSize);
+ rptA[bin]+=ptin[i]/(j->GetPt(0));
+ if (inJet[i] == 1) rptS[bin]+=ptin[i]/(j->GetPt(0));
+ if (fPythia && inJet[i] == -1)
+ rptR[bin]+=ptin[i]/(j->GetPt(0));
+ }
+ }
- TCanvas* c9 = new TCanvas("c9");
- e6H->Draw();
+ // compute shape and fill histogram
+ Float_t ptS,ptR,ptA,r;
+ ptS=ptR=ptA=0.0;
+ for (Int_t i=0;i<nBins;i++) {
+ ptS+=rptS[i];
+ ptR+=rptR[i];
+ ptA+=rptA[i];
+ r=(i+1)*binSize-halfBin;
+ hs->Fill(r,ptS*fWeight);
+ if(fPythia) {
+ hr->Fill(r,ptR*fWeight);
+ ha->Fill(r,ptA*fWeight);
+ }
+ }
+}
+
+void AliJetAnalysis::FillFragH()
+{
+ // Fill fragmentation histogram
+ if (fDoRecJ && fRecJ->GetNJets()> 0) {
+ TArrayF p=fRecJ->GetPtFromSignal();
+ if (p[0]>fPercentage) {
+ FragFun(fRecJ,fRFragSelH,fRFragRejH,fRFragAllH);
+ fWFragR+=fWeight;
+ }
+ }
+}
+void AliJetAnalysis::FragFun(AliJet *j,TH1F* hs, TH1F* hr, TH1F* ha)
+{
+ // Calculate de fragmentation function
+ TArrayI inJet = j->GetInJet();
+ TArrayF etain = j->GetEtaIn();
+ TArrayF ptin = j->GetPtIn();
+ Float_t xi,t1,ene;
- TCanvas* c10 = new TCanvas("c10");
+ for(Int_t i=0;i<(inJet.GetSize());i++) {
+ if (inJet[i] == 1 || (inJet[i] == -1 && fPythia)) {
+ t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
+ ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
+ xi = (Float_t) TMath::Log((j->GetE(0))/ene);
+ if (fPythia) ha->Fill(xi,fWeight);
+ if (inJet[i] == 1) hs->Fill(xi,fWeight);
+ if (fPythia && inJet[i] == -1) hr->Fill(xi,fWeight);
+ }
+ }
+}
+void AliJetAnalysis:: FillTrigH()
+{
+ // Fill trigger bias histograms
+ if(fDoTrig && fGenJ->GetNJets()>0 && fRecJ->GetNJets()>0){
+ TArrayF p=fRecJ->GetPtFromSignal();
+ if (p[0]>fPercentage) {
+ float genEne = fGenJ->GetE(0);
+ float recEne = fRecJ->GetE(0);
+ float eneRatio = (recEne/fEfactor)/genEne;
+
+ fGTriggerEneH->Fill(genEne, eneRatio, fWeight);
+ fRTriggerEneH->Fill(recEne/fEfactor, eneRatio, fWeight);
+
+ if (fPart->LeadingFound()){
+ float leaEne = fPart->GetE();
+ float eneRatio2 = leaEne/genEne;
+ fGPTriggerEneH->Fill(genEne, eneRatio2, fWeight);
+ fPTriggerEneH->Fill(leaEne/0.2, eneRatio2, fWeight);
+ }
+ }
+ }
+}
+
+
+////////////////////////////////////////////////////////////////////////
+// Normalize histogrames
+
+void AliJetAnalysis::NormHistograms()
+{
+ // normalize shape histograms
+ if (fDoShap) {
+ if (fDoRecJ && fWShapR>0.) { // leading reconstructed jet
+ fRShapSelH->Scale(1.0/fWShapR);
+ fRShapRejH->Scale(1.0/fWShapR);
+ fRShapAllH->Scale(1.0/fWShapR);
+ }
+ }
- r5H->SetMaximum(1.);
- r5H->Draw();
- r5H->SetXTitle("E_{gen} (GeV)");
- r5H->SetYTitle("E_{leading} / E_{gen}");
- r6H->SetLineColor(2);
- r6H->Draw("same");
-
- TCanvas* c11 = new TCanvas("c11");
-//
-// Leading particle rec/gen
-//
- r7H->SetMaximum(1.);
-
- r7H->Draw();
- r7H->SetXTitle("E_{gen} (GeV)");
- r7H->SetYTitle("E_{leading} / E_{gen}");
+ // normalize fragmentation function histograms
+ if (fDoFrag) {
+ if (fDoRecJ && fWFragR>0.) { // leading reconstructed jet
+ fRFragSelH->Scale(2.0/fWFragR);
+ fRFragRejH->Scale(2.0/fWFragR);
+ fRFragAllH->Scale(2.0/fWFragR);
+ }
+ }
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJetAnalysis::PlotHistograms()
+{
+ // Plot histogramas (to be done...)
+ if (fDoKine) PlotKineH();
+ if (fDoCorr) PlotCorrH();
+ if (fDoCorr50) PlotCorr50H();
+ if (fDoShap) PlotShapH();
+ if (fDoFrag) PlotFragH();
+ if (fDoTrig) PlotTrigH();
+}
+
+void AliJetAnalysis::PlotKineH() const
+{
+ // missing
+ if (fDoPart) ;
+ if (fDoGenJ) ;
+ if (fDoRecJ) ;
+}
+
+void AliJetAnalysis::PlotCorrH() const
+{
+ // missing
+ if (fDoPart && fDoGenJ) ;
+ if (fDoPart && fDoRecJ) ;
+ if (fDoGenJ && fDoRecJ) ;
+}
+void AliJetAnalysis::PlotCorr50H() const
+{
+ // missing
+ if (fDoPart && fDoGenJ) ;
+ if (fDoPart && fDoRecJ) ;
+ if (fDoGenJ && fDoRecJ) ;
+}
+
+void AliJetAnalysis::PlotShapH() const
+{
+ // missing
+ if (fDoGenJ) ;
+ if (fDoRecJ) ;
+}
+
+void AliJetAnalysis::PlotFragH() const
+{
+ // missing
+ if (fDoGenJ) ;
+ if (fDoRecJ) ;
+}
+
+void AliJetAnalysis::PlotTrigH()
+{
+ // missing
+
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliJetAnalysis::SaveHistograms()
+{
+ // Save histograms
+ TFile *fOut = new TFile(fFile,"recreate");
+ fOut->cd();
+ if (fDoKine) SaveKineH();
+ if (fDoCorr) SaveCorrH();
+ if (fDoCorr50) SaveCorr50H();
+ if (fDoShap) SaveShapH();
+ if (fDoFrag) SaveFragH();
+ if (fDoTrig) SaveTrigH();
+ if (fDoJt) SaveJtH();
+ if (fDodNdxi) SavedNdxiH();
+ fOut->Close();
+}
+
+void AliJetAnalysis::SaveKineH()
+{
+ // Save kinematic histograms
+ if (fDoPart) {
+ fPKineEneH->Write();
+ fPKinePtH->Write();
+ fPKinePhiH->Write();
+ fPKineEtaH->Write();
+ }
- r8H->SetLineColor(2);
- r8H->Draw("same");
+ if (fDoGenJ) {
+ fGKineEneH->Write();
+ fGKinePtH->Write();
+ fGKinePhiH->Write();
+ fGKineEtaH->Write();
+ }
- TCanvas* c12 = new TCanvas("c12");
- drP1->SetXTitle("E_{gen} (GeV)");
- drP1->SetYTitle("#eta_{rec} - #eta_{gen}");
- drP1->Draw();
-
- TCanvas* c13 = new TCanvas("c13");
- drP2->SetXTitle("E_{gen} (GeV)");
- drP2->SetYTitle("#eta_{leading} - #eta_{gen}");
+ if (fDoRecJ) {
+ fRKineEneH->Write();
+ fRKinePtH->Write();
+ fRKinePhiH->Write();
+ fRKineEtaH->Write();
+ }
+}
+
+void AliJetAnalysis::SaveCorrH()
+{
+ // Save correlation histograms
+ if (fDoPart && fDoGenJ) {
+ fPGCorrEneH->Write();
+ fPGCorrPtH->Write();
+ fPGCorrEtaH->Write();
+ fPGCorrPhiH->Write();
+ }
- drP2->Draw();
+ if (fDoPart && fDoRecJ) {
+ fPRCorrEneH->Write();
+ fPRCorrPtH->Write();
+ fPRCorrEtaH->Write();
+ fPRCorrPhiH->Write();
+ }
- TCanvas* c14 = new TCanvas("c14");
- dphiH->Draw();
-
-/*
- e1H->Write();
- e2H->Write();
- e3H->Write();
- e4H->Write();
-*/
-}
+ if (fDoGenJ && fDoRecJ) {
+ fRGCorrEneH->Write();
+ fRGCorrPtH->Write();
+ fRGCorrEtaH->Write();
+ fRGCorrPhiH->Write();
+ }
+}
+
+void AliJetAnalysis::SaveCorr50H()
+{
+ // Save correlation histograms (special case)
+ if (fDoPart && fDoRecJ) {
+ fPRCorr50EneH->Write();
+ fPRCorr50PtH->Write();
+ fPRCorr50EtaH->Write();
+ fPRCorr50PhiH->Write();
+ }
+ if (fDoGenJ && fDoRecJ) {
+ fRGCorr50EneH->Write();
+ fRGCorr50PtH->Write();
+ fRGCorr50EtaH->Write();
+ fRGCorr50PhiH->Write();
+ }
+}
+
+void AliJetAnalysis::SaveShapH()
+{
+ // Save jet shape histograms
+ if (fDoRecJ) {
+ fRShapSelH->Write();
+ fdEdrH->Write();
+ if(fDoBkgd) fdEdrB->Write();
+ fPtEneH2->Write();
+ fdEdrW->Write();
+ if (fPythia){
+ fRShapRejH->Write();
+ fRShapAllH->Write();
+ }
+ }
+}
+
+void AliJetAnalysis::SaveJtH()
+{
+ // Save J_T histograms
+ if (fDoRecJ) {
+ fJtH->Write();
+ if(fDoBkgd) fJtB->Write();
+ fJtW->Write();
+ }
+}
+
+void AliJetAnalysis::SavedNdxiH()
+{
+ // Save dN/d#xi histograms
+ if (fDoRecJ) {
+ fdNdxiH->Write();
+ if(fDoBkgd) fdNdxiB->Write();
+ fPtEneH->Write();
+ fdNdxiW->Write();
+ }
+}
+
+void AliJetAnalysis::SaveFragH()
+{
+ // Save fragmentation histograms
+ if (fDoRecJ) {
+ fRFragSelH->Write();
+ if(fPythia){
+ fRFragRejH->Write();
+ fRFragAllH->Write();
+ }
+ }
+}
+
+void AliJetAnalysis::SaveTrigH()
+{
+ // Save trigger bias histograms
+ if(fDoTrig){
+ fGTriggerEneH->Write();
+ fRTriggerEneH->Write();
+ fGPTriggerEneH->Write();
+ fPTriggerEneH->Write();
+ }
+}
+
+////////////////////////////////////////////////////////////////////////
+// main Analysis function
+
+void AliJetAnalysis::Analyze()
+
+{
+ // Kinematics for
+ // leading particle
+ // leading generated jet
+ // leading reconstructed jet
+
+ // Correlations amd resolutions
+ // a) correlations in energy, pt, phi, eta
+ // b) resolutions in energy, pt, phi, eta, r
+ // leading particle and leading generated jet
+ // leading particle and leading reconstructed jet
+ // leading generated jet and leading reconstructed jet
+
+ // Fragmentation functions and Shapes
+ // a) integrated over all pt
+ // b) in 3 flavors:
+ // b.1) only for user selected particles in jet
+ // b.2) only for user rejected particles in jet
+ // b.3) for all particles in jet
+
+ DefineHistograms();
+ FillHistograms();
+ NormHistograms();
+ PlotHistograms();
+ SaveHistograms();
+}
/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* See cxx source for full Copyright notice */
-
-
+
//---------------------------------------------------------------------
// JetAnalysis class
-// Perform Jet Analysis
-// Author: amdreas.morsch@cern.ch
+// Perform Jet Analysis on already found jets
+// Author: andreas.morsch@cern.ch, jgcn@mail.cern.ch
+// mercedes.lopez.noriega@cern.ch
//---------------------------------------------------------------------
+
#include <TObject.h>
+class AliLeading;
+class AliJet;
+class TH1;
+class TH1F;
+class TH2F;
+class TProfile;
+class TLorentzVector;
+
class AliJetAnalysis : public TObject
{
public:
AliJetAnalysis();
- ~AliJetAnalysis(){;}
+ virtual ~AliJetAnalysis();
- void Analyse();
- // Setter
- void SetDirectory(char* directory) {fDirectory = directory;}
- void SetEventRange(Int_t imin, Int_t imax) {fEventMin = imin; fEventMax = imax;}
- void SetRunRange(Int_t imin, Int_t imax) {fRunMin = imin; fRunMax = imax;}
+ void Analyze();
+ // define histograms
+ void DefineHistograms();
+ void DefineKineH();
+ void DefineCorrH();
+ void DefineCorr50H();
+ void DefineShapH();
+ void DefineFragH();
+ void DefineTrigH();
+ void DefineJtH();
+ void DefinedNdxiH();
+ // fill histograms
+ void FillHistograms();
+ void FillKineH();
+ void FillCorrH();
+ void FillCorr50H();
+ void FillShapH(Float_t r);
+ void FillFragH();
+ void FillTrigH();
+ void FillJtH();
+ void FilldNdxiH();
+ void FillBkgd(Int_t eventN, Int_t runN);
+ // normalize histograms
+ void NormHistograms();
+ // plot histograms
+ void PlotHistograms();
+ void PlotKineH() const;
+ void PlotCorrH() const;
+ void PlotCorr50H() const;
+ void PlotShapH() const;
+ void PlotFragH() const;
+ void PlotTrigH();
+ // save histograms
+ void SaveHistograms();
+ void SaveKineH();
+ void SaveCorrH();
+ void SaveCorr50H();
+ void SaveShapH();
+ void SaveFragH();
+ void SaveTrigH();
+ void SaveJtH();
+ void SavedNdxiH();
+ // other functions
+ void Shape(AliJet *j,TH1F* hs, TH1F* hr, TH1F* ha, TH2F* hd, TH2F* hp, TH1F* wd, Float_t r);
+ void FragFun(AliJet *j,TH1F* hs, TH1F* hr, TH1F* ha);
+ void Correlation(TLorentzVector *lv1,TLorentzVector *lv2,TH2F *h1, TH2F *h2, TH2F *h3, TH2F *h4);
+ void Correlation50(AliJet *j,TLorentzVector *lv1,TLorentzVector *lv2,TH2F *h1, TH2F *h2, TH2F *h3, TH2F *h4);
+ // setters
+ void SetDirectory(char* directory)
+ {fDirectory = directory;} // directory where file with jets is
+ void SetBkgdDirectory(char* directory)
+ {fBkgdDirectory = directory;} // directory where file with background is
+ void SetOutputFile(char* file) {fFile = file;} // file where plots will be saved
+ void SetPercentage(Float_t p) {fPercentage = p;} // minimum percentage of tracks coming from pythia (very aprox.)
+ void SetEventRange(Int_t imin, Int_t imax)
+ {fEventMin = imin; fEventMax = imax;} // first and last event
+ void SetRunRange(Int_t imin, Int_t imax)
+ {fRunMin = imin; fRunMax = imax;} // first and last run
+ void SetMinimumMult(Int_t m){fminMult = m;} // minimum multiplicity cut
+ void SetPythia(Bool_t f = kFALSE){fPythia = f;} // If only pythia, to save everything...
+ void SetDoJt(Bool_t f = kTRUE){fDoJt = f;} // To get j_T distribution
+ void SetDodNdxi(Bool_t f = kTRUE){fDodNdxi = f;} // To get #xi distribution
+ void SetDoBkgd(Bool_t f = kTRUE) {fDoBkgd = f;} // To get the bkgd for j_T, xi and dEdr in a hijing event
+ void SetDoLeadPart(Bool_t f = kTRUE){fDoPart = f;}// To make plots for leading particle
+ void SetDoGenJet(Bool_t f = kTRUE){fDoGenJ = f;} // To make plots for generated jets
+ void SetDoRecJet(Bool_t f = kTRUE){fDoRecJ = f;} // To make plots for reconstructed jets
+ void SetDoKinematics(Bool_t f = kTRUE){fDoKine = f;} // To make the kine plots
+ void SetDoCorrelations(Bool_t f = kTRUE){fDoCorr = f;} // Correlation histograms
+ void SetDoCorr50(Bool_t f = kFALSE){fDoCorr50 = f;} // Correlation histograms when one particle has more than 50% E
+ void SetDoShape(Bool_t f = kTRUE){fDoShap = f;} // Shape plots
+ void SetDoFragmentations(Bool_t f = kTRUE){fDoFrag = f;} // Fragmentation
+ void SetDoTriggerBias(Bool_t f = kTRUE){fDoTrig = f;} // Trigger bias plots
+ void SetDivideEnergy(Float_t Efactor){fEfactor = Efactor;} // Divides E of rec.jet by Efactor
+ void SetProperties(TH1* h,const char* x, const char* y) const;
+ void SetReaderHeader(char *s="AliJetKineReaderHeader"){fReaderHeader = s;}
+ void SetdEdrWeight();
+ void SetPartPtCut(Float_t c){fPartPtCut = c;}
+ void SetdrJt(Float_t r){fdrJt = r;}
+ void SetdrdNdxi(Float_t r){fdrdNdxi = r;}
+ void SetdrdEdr(Float_t r){fdrdEdr = r;}
+ // getters
+ Float_t GetdEdrWeight(Float_t eta, Float_t r);
+
private:
- char* fDirectory; // Directory
- Int_t fEventMin; // Minimum event number
- Int_t fEventMax; // Maximum event number
- Int_t fRunMin; // Minimum run number
- Int_t fRunMax; // Maximum run number
+ char* fReaderHeader; // Reader header
+ char* fDirectory; // Directory
+ char* fBkgdDirectory; // Directory for background
+ char* fFile; // Output file name
+ Int_t fEventMin; // Minimum event number
+ Int_t fEventMax; // Maximum event number
+ Int_t fRunMin; // Minimum run number
+ Int_t fRunMax; // Maximum run number
+ Int_t fminMult; // Minimum multiplicity for events
+ Float_t fPercentage; // percentage of pt from signal particles to accept a jet
+ Float_t fPartPtCut; // cut in the pt of particles in histos
+ Float_t fdrJt; // maximum dr for Jt plot
+ Float_t fdrdNdxi; // maximum dr for dN/dxi plot
+ Float_t fdrdEdr; // maximum dr for dE/dr plot
+ Float_t fEfactor; // factor by which energy the reconstructed jet will be divided
+
+ Float_t fp0; // percentage of tracks in reconstructed jet coming from pythia
+ // so far calculated in aprox. way, it needs to be improved!
+ // for background from hijing events:
+ Float_t fPtJ; // P_T of the pythia jet
+ Float_t fEJ; // Energy of the pythia jet
+ Float_t fEtaJ; // Eta of the pythia jet
+ Float_t fPhiJ; // Phi of the pythia jet
+ Float_t fjv3X, fjv3Y, fjv3Z; // x,y,z of the pythia jet
+
+ // user options
+ Bool_t fPythia; // if pythia events
+ Bool_t fDoPart; // do analysis of leading particle
+ Bool_t fDoGenJ; // do analysis of leading generated jet
+ Bool_t fDoRecJ; // do analysis of leading rec jet
+ Bool_t fDoKine; // do kinematic plots
+ Bool_t fDoCorr; // do correlation plots
+ Bool_t fDoCorr50; // do correlation plots when one track more than 50% of jet energy
+ Bool_t fDoShap; // do shape plots
+ Bool_t fDoFrag; // do fragmentation plots
+ Bool_t fDoTrig; // do trigger bias plots
+ Bool_t fDoJt; // do jt histo
+ Bool_t fDodNdxi; // do dN/dxi histo
+ Bool_t fDoBkgd; // get dN/dxi bkgd using hijing tracks only
-
+ // weights
+ Float_t fWeight; // event weight
+ Float_t fWShapR; // weighted number of jets
+ Float_t fWFragR; // weighted number of jets
+ Float_t fWeightdEdr[10][20]; // weight for acceptance of dE/dr histo
+ Float_t fWdEdr; // weighted number of events for dEdr histo
+ Float_t fWJt; // weight for Jt
+ Float_t fWdNdxi; // weight fro dNd#xi
+
+ // leading hets and particles
+ AliLeading* fPart; // pointer to leading particle
+ AliJet* fGenJ; // pointer to leading generated jet
+ AliJet* fRecJ; // pointer to leading reconstructed jet
+ AliJet* fRecB; // pointer to leading reconstructed jet for background
+
+ // kine histos
+ TH1F *fRKineEneH; // Reconstructed energy histo
+ TH1F *fRKinePtH; // Reconstructed Pt histo
+ TH1F *fRKinePhiH; // Reconstructed phi histo
+ TH1F *fRKineEtaH; // Reconstructed eta histo
+ TH1F *fGKineEneH; // Generated energy histo
+ TH1F *fGKinePtH; // Generated Pt histo
+ TH1F *fGKinePhiH; // Generated phi histo
+ TH1F *fGKineEtaH; // Generated eta histo
+ TH1F *fPKineEneH; // Pythia energy histo
+ TH1F *fPKinePtH; // Pythia Pt histo
+ TH1F *fPKinePhiH; // Pythia phi histo
+ TH1F *fPKineEtaH; // Pythia eta histo
+
+ // correlation histograms
+ TH2F *fPGCorrEneH; // Energy correlation part-gen jet
+ TH2F *fPGCorrPtH; // Pt correlation part-gen jet
+ TH2F *fPGCorrEtaH; // Pseudorapidity correlation part-gen jet
+ TH2F *fPGCorrPhiH; // Azimuthal angle correlation part-gen jet
+ TH2F *fPRCorrEneH; // Energy correlation part-rec jet
+ TH2F *fPRCorrPtH; // Pt correlation part-rec jet
+ TH2F *fPRCorrEtaH; // Pseudorapidity correlation part-rec jet
+ TH2F *fPRCorrPhiH; // Azimuthal angle correlation part-rec jet
+ TH2F *fRGCorrEneH; // Energy correlation rec jet-gen jet
+ TH2F *fRGCorrPtH; // Pt correlation rec jet-gen jet
+ TH2F *fRGCorrEtaH; // Pseudorapidity correlation rec jet-gen jet
+ TH2F *fRGCorrPhiH; // Azimuthal angle correlation rec jet-gen jet
+
+ // correlation histograms
+ TH2F *fPRCorr50EneH; // Energy correlation part-rec jet
+ TH2F *fPRCorr50PtH; // Pt correlation part-rec jet
+ TH2F *fPRCorr50EtaH; // Pseudorapidity correlation part-rec jet
+ TH2F *fPRCorr50PhiH; // Azimuthal angle correlation part-rec jet
+ TH2F *fRGCorr50EneH; // Energy correlation rec jet-gen jet
+ TH2F *fRGCorr50PtH; // Pt correlation rec jet-gen jet
+ TH2F *fRGCorr50EtaH; // Pseudorapidity correlation rec jet-gen jet
+ TH2F *fRGCorr50PhiH; // Azimuthal angle correlation rec jet-gen jet
+
+ // fragmentation function and shape histos
+ TH1F *fRFragSelH; // Frag Fun of reconstructed jets (sel part)
+ TH1F *fRFragRejH; // Frag Fun of reconstructed jets (rej part)
+ TH1F *fRFragAllH; // Frag Fun of reconstructed jets (all part)
+ TH1F *fRShapSelH; // Shape of generated jets (sel part)
+ TH1F *fRShapRejH; // Shape of generated jets (rej part)
+ TH1F *fRShapAllH; // Shape of generated jets (all part)
+
+ // trigger bias histos
+ TProfile *fGTriggerEneH; // Generated energy (trigger bias)
+ TProfile *fRTriggerEneH; // Reconstructed energy (trigger bias)
+ TProfile *fGPTriggerEneH; // Generated energy (trigger bias)
+ TProfile *fPTriggerEneH; // Leading particle energy (trigger bias)
+
+ // dE/dr histo
+ TH2F* fdEdrH; // dE/dr histo
+ TH2F* fdEdrB; // dE/dr bkgdhisto
+ TH2F* fPtEneH2;// fPtEneH2
+ TH1F* fdEdrW; // weights for dE/dr
+
+ // Jt histo
+ TH2F* fJtH; // J_{T} histogram
+ TH2F* fJtB; // J_{T} bkgd histogram
+ TH1F* fJtW; // J_{T} weight
+
+ // dN/dxi histo
+ TH2F* fdNdxiH; // dN/d#xi histo
+ TH2F* fdNdxiB; // dN/d#xi bkgd histo
+ TH1F* fdNdxiW; // dN/d#xi weight histo
+ TH2F* fPtEneH; // fPtEneH
+
ClassDef(AliJetAnalysis,1)
};
#include <TStyle.h>
#include <TCanvas.h>
-#include <TLorentzVector.h>
-#include <TFile.h>
-#include <TClonesArray.h>
#include <TH1I.h>
#include <TH1D.h>
#include "AliJetReader.h"
AliJetControlPlots::AliJetControlPlots()
{
- //
// Constructor
- //
fNJetT=0;
+ // general properties
fNJetsH = new TH1I("fNJetsH","Number of Jets",12,0,11);
SetProperties(fNJetsH,"Number of jets","Entries");
fMultH = new TH1I("fMultH","Multiplicity of Jets",22,0,21);
SetProperties(fMultH,"Multiplicity of jets","Entries");
+ fInJetH = new TH1D("fInJetH","Percentage of particles in jets",50,0,1);
+ SetProperties(fInJetH,"percentage of particles in jets","Entries");
+
+ // kinematics
fPtH = new TH1D("fPtH","Pt of Jets",50,0.,200.);
SetProperties(fPtH,"P_{t} [GeV]","Entries");
fEneH = new TH1D("fEneH","Energy of Jets",50,0.,200.);
SetProperties(fEneH,"Energy [GeV]","Entries");
- fFragH = new TH1D("fFragH","Jet Fragmentation",20,0.,1.);
- SetProperties(fFragH,"x=E_{i}/E_{JET}","1/N_{JET}dN_{ch}/dx");
-
- fFragLnH = new TH1D("fFragLnH","Jet Fragmentation",20,0.,10);
- SetProperties(fFragLnH,"#xi=ln(E_{JET}/E_{i}","1/N_{JET}dN_{ch}/d#xi");
-
fPhiH = new TH1D("fPhiH","Azimuthal angle of Jets",
- 60,-TMath::Pi(),TMath::Pi());
+ 60,0.,2.0*TMath::Pi());
SetProperties(fPhiH,"#phi","Entries");
-
- fInJetH = new TH1D("fInJetH","Percentage of particles in jets",
- 50,0,1);
- SetProperties(fInJetH,"percentage of particles in jets","Entries");
+
+ // fragmentation
+ fFragH = new TH1D("fFragH","Leading Jet Fragmentation (selected tracks)",
+ 20,0.,1.);
+ SetProperties(fFragH,"x=E_{i}/E_{JET}","1/N_{JET}dN_{ch}/dx");
+
+ fFragrH = new TH1D("fFragrH","Leading Jet Fragmentation (rejected tracks)",
+ 20,0.,1.);
+ SetProperties(fFragrH,"x=E_{i}/E_{JET}","1/N_{JET}dN_{ch}/dx");
+
+ // fragmentation log
+ fFragLnH = new TH1D("fFragLnH","Leading Jet Fragmentation (selected tracks)",
+ 20,0.,10);
+ SetProperties(fFragLnH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
+
+ fFragLnrH = new TH1D("fFragLnrH",
+ "Leading Jet Fragmentation (rejected tracks)",
+ 20,0.,10);
+ SetProperties(fFragLnrH,"#xi=ln(E_{JET}/E_{i})","1/N_{JET}dN_{ch}/d#xi");
+
+ // jet shape
+ fShapeH = new TH1D("fShapeH","Leading Jet Shape (selected tracks)",
+ 20,0.,1.);
+ SetProperties(fShapeH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R)");
+
+ fShaperH = new TH1D("fShaperH","Leading Jet Shape (rejected tracks)",
+ 20,0.,1.);
+ SetProperties(fShaperH,"r","1/N_{JET}#Sigma P_{T}(0,r)/P_{T}(0,R)");
}
////////////////////////////////////////////////////////////////////////
AliJetControlPlots::~AliJetControlPlots()
{
- //
// Destructor
- //
delete fNJetsH;
delete fMultH;
delete fPtH;
delete fEtaH;
delete fEneH;
- delete fFragH;
- delete fFragLnH;
delete fPhiH;
delete fInJetH;
+ delete fFragH;
+ delete fFragLnH;
+ delete fFragrH;
+ delete fFragLnrH;
+ delete fShapeH;
+ delete fShaperH;
}
////////////////////////////////////////////////////////////////////////
-void AliJetControlPlots::FillHistos(AliJet *j, AliJetReader *r)
+void AliJetControlPlots::FillHistos(AliJet *j)
{
-// Fills the histograms
-
+ // Fills the histograms
Int_t nj = j->GetNJets();
- fNJetsH->Fill(nj);
- fNJetT+=nj;
-
+ fNJetsH->Fill(nj,1);
+ if (nj == 0) return;
+
// kinematics, occupancy and multiplicities
TArrayI mj = j->GetMultiplicities();
Int_t mjTot=0;
for (Int_t i=0;i<nj;i++) {
mjTot+=mj[i];
- fMultH->Fill(mj[i]);
- fPtH->Fill(j->GetPt(i));
- fEtaH->Fill(j->GetEta(i));
- fEneH->Fill(j->GetE(i));
- fPhiH->Fill(j->GetPhi(i));
+ fMultH->Fill(mj[i],1);
+ fPtH->Fill(j->GetPt(i),1);
+ fEtaH->Fill(j->GetEta(i),1);
+ fEneH->Fill(j->GetE(i),1);
+ fPhiH->Fill(j->GetPhi(i),1);
}
- if (nj>0)
- fInJetH->Fill(((Double_t) mjTot)/((Double_t) j->GetNinput()));
-
- // fragmentation
- TClonesArray *lvArray = r->GetMomentumArray();
+ fInJetH->Fill(((Double_t) mjTot)/((Double_t) j->GetNinput()),1);
+
+ // fragmentation of leading jet
TArrayI inJet = j->GetInJet();
- Int_t nIn = j->GetNinput();
- for(Int_t i=0;i<nIn;i++) {
- if (inJet[i] == -1) continue;
- TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
- Double_t xe = (lv->E())/(j->GetE(inJet[i]));
- fFragH->Fill(xe);
- fFragLnH->Fill(TMath::Log(1.0/xe));
+ TArrayF etain = j->GetEtaIn();
+ TArrayF ptin = j->GetPtIn();
+ for(Int_t i=0;i<(inJet.GetSize());i++) {
+ Float_t t1 = TMath::Tan(2.0*TMath::ATan(TMath::Exp(-etain[i])));
+ Float_t ene = ptin[i]*TMath::Sqrt(1.+1./(t1*t1));
+ if (inJet[i] == 1) {
+ fFragH->Fill((Float_t) ene/(j->GetE(0)),1);
+ fFragLnH->Fill((Float_t) TMath::Log((j->GetE(0))/ene),1);
+ }
+ if (inJet[i] == -1) {
+ fFragrH->Fill((Float_t) ene/(j->GetE(0)),1);
+ fFragLnrH->Fill((Float_t) TMath::Log((j->GetE(0))/ene),1);
+ }
+ }
+
+ // shape of leading jet
+ // * CAREFUL: depends on number of bins and bin sizes
+ // of shape histo. HARDWIRED at the moment *
+ // Plot against dr NOT dr/R_jet!!!
+ TArrayF phiin = j->GetPhiIn();
+ Float_t rptS[20], rptR[20];
+ for(Int_t i=0;i<20;i++) rptS[i]=rptR[i]=0.0;
+
+ for(Int_t i=0;i<inJet.GetSize();i++) {
+ if (inJet[i] == 1 || inJet[i] == -1) {
+ Float_t deta = etain[i] - j->GetEta(0);
+ Float_t dphi = phiin[i] - j->GetPhi(0);
+ if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+ if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+ Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
+ if (dr>1) continue;
+ Int_t bin = (Int_t) TMath::Floor(dr/0.05);
+ if (inJet[i] == 1) rptS[bin]+=ptin[i]/(j->GetPt(0));
+ if (inJet[i] == -1) rptR[bin]+=ptin[i]/(j->GetPt(0));
+ }
+ }
+
+ Float_t ptS,ptR,r;
+ ptS=ptR=0.0;
+ for (Int_t i=0;i<20;i++) {
+ ptS+=rptS[i];
+ ptR+=rptR[i];
+ r=(i+1)*0.05-0.025;
+ fShapeH->Fill(r,ptS);
+ fShaperH->Fill(r,ptR);
}
+ fNJetT++;
}
////////////////////////////////////////////////////////////////////////
void AliJetControlPlots::Normalize()
{
+// *CAREFUL: depends on histogram number of bins
if (fNJetT == 0) return;
- fFragH->Sumw2();
fFragH->Scale(20.0/((Double_t) fNJetT));
- fFragLnH->Sumw2();
fFragLnH->Scale(2.0/((Double_t) fNJetT));
+ fFragrH->Scale(20.0/((Double_t) fNJetT));
+ fFragLnrH->Scale(2.0/((Double_t) fNJetT));
+ fShapeH->Scale(1.0/((Double_t) fNJetT));
+ fShaperH->Scale(1.0/((Double_t) fNJetT));
}
////////////////////////////////////////////////////////////////////////
gStyle->SetOptStat(kFALSE);
gStyle->SetOptTitle(kFALSE);
- TCanvas *c = new TCanvas("c","AliJetControlPlots",50,200,900,700);
- c->Divide(3,3);
+ TCanvas *c = new TCanvas("c","AliJetControlPlots",50,200,1200,700);
+ c->Divide(4,3);
c->cd(1); fNJetsH->Draw("e1");
c->cd(2); fMultH->Draw("e1");
c->cd(3); fInJetH->Draw("e1");
c->cd(5); fEtaH->Draw("e1");
c->cd(6); fPhiH->Draw("e1");
c->cd(7); fEneH->Draw("e1");
- c->cd(8); fFragH->Draw("e1");
- c->cd(9); fFragLnH->Draw("e1");
+ c->cd(8); fFragLnH->Draw("e1");
+ c->cd(9); fFragLnrH->Draw("e1");
+ c->cd(10); fShapeH->Draw("e1");
+ c->cd(11); fShaperH->Draw("e1");
}
////////////////////////////////////////////////////////////////////////
void AliJetControlPlots::SetProperties(TH1* h,const char* x, const char* y) const
{
-//
// Sets the histogram style properties
h->SetMarkerStyle(20);
h->SetMarkerSize(.5);
h->SetMarkerColor(2);
h->SetXTitle(x);
h->SetYTitle(y);
+ h->Sumw2();
}
// getters
TH1I *GetNJetsH() {return fNJetsH;}
TH1I *GetMultH() {return fMultH;}
+ TH1D *GetPhiH() {return fPhiH;}
+ TH1D *GetFractionInJetH() {return fInJetH;}
TH1D *GetEneH() {return fEneH;}
TH1D *GetPtH() {return fPtH;}
TH1D *GetEtaH() {return fEtaH;}
TH1D *GetFragH() {return fFragH;}
TH1D *GetFragLnH() {return fFragLnH;}
- TH1D *GetPhiH() {return fPhiH;}
- TH1D *GetFractionInJetH() {return fInJetH;}
+ TH1D *GetFragrH() {return fFragrH;}
+ TH1D *GetFragLnrH() {return fFragLnrH;}
+ TH1D *GetShapeH() {return fShapeH;}
+ TH1D *GetShaperH() {return fShaperH;}
// others
- void FillHistos(AliJet *j,AliJetReader *r);
+ void FillHistos(AliJet *j);
void PlotHistos();
void SetProperties(TH1* h,const char* x, const char* y) const;
void Normalize();
protected:
- TH1I *fNJetsH; // distribution of number of jets
+ TH1I *fNJetsH; // distribution of number of jets
TH1I *fMultH; // jet multiplicity
- TH1D *fPtH; // pt spectra
- TH1D *fEtaH; // eta distribution
- TH1D *fEneH; // energy distribution
- TH1D *fFragH; // jet fragmentation
- TH1D *fFragLnH; // jet fragmentation in ln scale
- TH1D *fPhiH; // phi distribution
- TH1D *fInJetH; // percentage of input particles in a jet
- Int_t fNJetT; // total number of jets for normalization
+ TH1D *fPtH; // pt spectra
+ TH1D *fEtaH; // eta distribution
+ TH1D *fEneH; // energy distribution
+ TH1D *fFragH; // leading jet fragmentation (selected part)
+ TH1D *fFragLnH; // leading jet fragmentation in ln scale
+ TH1D *fFragrH; // leading jet fragmentation (rejected part)
+ TH1D *fFragLnrH; // leading jet fragmentation in ln scale
+ TH1D *fShapeH; // leading jet shape (selected part)
+ TH1D *fShaperH; // leading jet shape (rejected part)
+ TH1D *fPhiH; // phi distribution
+ TH1D *fInJetH; // percentage of input particles in a jet
+ Int_t fNJetT; // total number of jets for normalization
ClassDef(AliJetControlPlots,1)
};
* provided "as is" without express or implied warranty. *
**************************************************************************/
-//
// Jet ESD Reader
// ESD reader for jet analysis
-// Author: Mercedes Lopez Noriega
-// mercedes.lopez.noriega@cern.ch
-//
-
+// Author: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
#include <Riostream.h>
#include <TSystem.h>
#include <TLorentzVector.h>
+#include <TVector3.h>
#include "AliJetESDReader.h"
#include "AliJetESDReaderHeader.h"
#include "AliESD.h"
AliJetESDReader::AliJetESDReader()
{
// Constructor
-
+ printf("\nIn reader constructor\n");
fReaderHeader = 0x0;
fMass = 0;
fSign = 0;
AliJetESDReader::~AliJetESDReader()
{
// Destructor
- // delete fReaderHeader;
}
//____________________________________________________________________________
void AliJetESDReader::OpenInputFiles()
-
{
+ // Open input files
+ printf("\nOpening files\n");
// chain for the ESDs
fChain = new TChain("esdTree");
fChainMC = new TChain("mcStackTree");
printf("Adding %s\n",name);
char path[256];
sprintf(path,"%s/%s",dirName,name);
- fChain->AddFile(path,-1,"esdTree");
- fChainMC->AddFile(path,-1,"mcStackTree");
+ fChain->AddFile(path,-1);
+ fChainMC->AddFile(path,-1);
}
}
void AliJetESDReader::FillMomentumArray(Int_t event)
{
-// Fills the momentum array
+ // Fill the momentum array for each track
Int_t goodTrack = 0;
Int_t nt = 0;
- Float_t pt;
- Double_t p[3]; // track 3 momentum vector
+ Float_t pt, eta;
+ TVector3 p3;
// clear momentum array
ClearArray();
+
// get event from chain
fChain->GetEntry(event);
fChainMC->GetEntry(event);
+
// get number of tracks in event (for the loop)
nt = fESD->GetNumberOfTracks();
- // tmporary storage of signal flags
- Int_t* flag = new Int_t[nt];
+
+ // temporary storage of signal and pt cut flag
+ Int_t* sflag = new Int_t[nt];
+ Int_t* cflag = new Int_t[nt];
// get cuts set by user
- Float_t ptMin = ((AliJetESDReaderHeader*) fReaderHeader)->GetPtCut();
+ Float_t ptMin = fReaderHeader->GetPtCut();
+ Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
+ Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
//loop over tracks
for (Int_t it = 0; it < nt; it++) {
AliESDtrack *track = fESD->GetTrack(it);
UInt_t status = track->GetStatus();
- if ((status & AliESDtrack::kITSrefit) == 0) continue; // quality check
-
+ p3 = track->P3();
+ pt = p3.Pt();
+ if (((status & AliESDtrack::kITSrefit) == 0) ||
+ ((status & AliESDtrack::kTPCrefit) == 0)) continue; // quality check
if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
- && TMath::Abs(track->GetLabel()) > 10000) continue;
-
- track->GetPxPyPz(p);
- pt = TMath::Sqrt(p[0] * p[0] + p[1] * p[1]); // pt of the track
-
- if (pt < ptMin) continue; //check cuts
-
- new ((*fMomentumArray)[goodTrack])
- TLorentzVector(p[0], p[1], p[2],
- TMath::Sqrt(pt * pt +p[2] * p[2]));
- flag[goodTrack]=0;
- if (TMath::Abs(track->GetLabel()) < 10000) flag[goodTrack]=1;
+ && TMath::Abs(track->GetLabel()) > 10000) continue; // quality check
+ if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
+ && TMath::Abs(track->GetLabel()) < 10000) continue; // quality check
+ eta = p3.Eta();
+ if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
+
+ new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
+ sflag[goodTrack]=0;
+ if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
+ cflag[goodTrack]=0;
+ if (pt > ptMin) cflag[goodTrack]=1; // pt cut
goodTrack++;
}
// set the signal flags
- fSignalFlag.Set(goodTrack,flag);
+ fSignalFlag.Set(goodTrack,sflag);
+ fCutFlag.Set(goodTrack,cflag);
- printf("\nIn event %d, number of good tracks %d \n", event, goodTrack);
+ //printf("\nIn event %d, number of good tracks %d \n", event, goodTrack);
}
/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* See cxx source for full Copyright notice */
+//---------------------------------------------------------------------
// Jet ESD Reader
// ESD reader for jet analysis
// Author: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
+//---------------------------------------------------------------------
#include "AliJetReader.h"
class AliJetESDReaderHeader;
virtual ~AliJetESDReader();
// Getters
- Float_t GetTrackMass() const {return fMass;} // returns mass of the track
- Int_t GetTrackSign() const {return fSign;} // returns sign of the track
+ Float_t GetTrackMass() const {return fMass;} // returns mass of the track
+ Int_t GetTrackSign() const {return fSign;} // returns sign of the track
// Setters
void FillMomentumArray(Int_t event);
* about the suitability of this software for any purpose. It is *
* provided "as is" without express or implied warranty. *
**************************************************************************/
-
-//
+
+//---------------------------------------------------------------------
// Jet ESD Reader Header
// Header for the ESD reader in the jet analysis
-// Author: Mercedes Lopez Noriega
-// mercedes.lopez.noriega@cern.ch
-//
+// Author: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
+//---------------------------------------------------------------------
#include "AliJetESDReaderHeader.h"
AliJetReaderHeader("AliJetESDReaderHeader")
{
// Constructor
- SetPtCut();
- SetDCA();
- SetTLength();
SetReadSignalOnly(kFALSE);
+ SetReadBkgdOnly(kFALSE);
}
virtual ~AliJetESDReaderHeader();
// Getters
- Float_t GetPtCut() const {return fPtCut;}
- Float_t GetDCA() const {return fDCA;} // not working so far..(always 0)
- Float_t GetTLength() const {return fTLength;} // not working so far.. (always 0)
Bool_t ReadSignalOnly() const {return fReadSignalOnly;}
-
+ Bool_t ReadBkgdOnly() const {return fReadBkgdOnly;}
// Setters
- virtual void SetPtCut(Float_t par = 2.0) {fPtCut = par;}
- virtual void SetDCA(Float_t dca = 3.0) {fDCA = dca;}
- virtual void SetTLength(Float_t length = 0.0) {fTLength = length;}
virtual void SetReadSignalOnly(Bool_t flag = kTRUE) {fReadSignalOnly = flag;}
+ virtual void SetReadBkgdOnly(Bool_t flag = kTRUE) {fReadBkgdOnly = flag;}
protected:
//parameters set by user
- Float_t fPtCut; // pt cut
- Float_t fDCA; // dca cut
- Float_t fTLength; // track length cut
Bool_t fReadSignalOnly; // read particles from signal event only
- ClassDef(AliJetESDReaderHeader,1);
+ Bool_t fReadBkgdOnly; // read particles from bkgd event only
+ ClassDef(AliJetESDReaderHeader,2);
};
#endif
* about the suitability of this software for any purpose. It is *
* provided "as is" without express or implied warranty. *
**************************************************************************/
-
-
+
//---------------------------------------------------------------------
// Jet finder base class
// manages the search for jets
#include "AliLeading.h"
#include "AliHeader.h"
-
ClassImp(AliJetFinder)
-////////////////////////////////////////////////////////////////////////
AliJetFinder::AliJetFinder()
{
- //
// Constructor
- //
fOut = 0x0;
fJets = new AliJet();
fGenJets = new AliJet();
SetPlotMode(kFALSE);
}
-
////////////////////////////////////////////////////////////////////////
AliJetFinder::~AliJetFinder()
{
- //
// destructor
- //
-
// here reset and delete jets
fJets->ClearJets();
delete fJets;
delete fOut;
// reset and delete control plots
if (fPlots) delete fPlots;
- // delete fLeading;
}
-
////////////////////////////////////////////////////////////////////////
void AliJetFinder::SetOutputFile(const char *name)
fOut = new TFile(name,"recreate");
}
-
////////////////////////////////////////////////////////////////////////
void AliJetFinder::PrintJets()
{
-//
-// Print jet information
+ // Print jet information
cout << " Jets found with jet algorithm:" << endl;
fJets->PrintJets();
cout << " Jets found by pythia:" << endl;
fGenJets->PrintJets();
}
-
////////////////////////////////////////////////////////////////////////
void AliJetFinder::SetPlotMode(Bool_t b)
{
-// Sets the plotting mode
+ // Sets the plotting mode
fPlotMode=b;
if (b && !fPlots) fPlots = new AliJetControlPlots();
}
void AliJetFinder::WriteJetsToFile(Int_t i)
{
-// Writes the jets to file
+ // Writes the jets to file
fOut->cd();
char hname[30];
sprintf(hname,"TreeJ%d",i);
rh->Write();
}
-
////////////////////////////////////////////////////////////////////////
void AliJetFinder::GetGenJets()
{
-// Get the generated jet information from mc header
+ // Get the generated jet information from mc header
AliHeader* alih = fReader->GetAliHeader();
if (alih == 0) return;
AliGenEventHeader * genh = alih->GenEventHeader();
WriteRHeaderToFile();
WriteJHeaderToFile();
}
-
// loop over events
Int_t nFirst,nLast;
nFirst = fReader->GetReaderHeader()->GetFirstEvent();
GetGenJets();
FindJets();
if (fOut) WriteJetsToFile(i);
- if (fPlots) fPlots->FillHistos(fJets,fReader);
+ if (fPlots) fPlots->FillHistos(fJets);
fLeading->Reset();
fGenJets->ClearJets();
Reset();
fOut->Close();
}
}
-
/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* See cxx source for full Copyright notice */
-
//---------------------------------------------------------------------
// Jet finder base class
// manages the search for jets
class AliJetControlPlots;
class AliLeading;
-
class AliJetFinder : public TObject
{
public:
virtual ~AliJetFinder();
// getters
- virtual AliJet *GetJets() {return fJets;}
- virtual Bool_t GetPlotMode() const {return fPlotMode;}
- virtual TFile* GetOutputFile() { return fOut; }
+ virtual AliJet *GetJets() {return fJets;}
+ virtual Bool_t GetPlotMode() const {return fPlotMode;}
+ virtual TFile* GetOutputFile() {return fOut;}
// setters
virtual void SetPlotMode(Bool_t b);
virtual void SetOutputFile(const char *name="jets.root");
protected:
Bool_t fPlotMode; // do you want control plots?
AliJet* fJets; // pointer to jet class
- AliJet* fGenJets; // pointer to generated jets
+ AliJet* fGenJets; // pointer to generated jets
AliLeading* fLeading; // pointer to leading particle data
AliJetReader* fReader; // pointer to reader
AliJetControlPlots* fPlots; // pointer to control plots
//
// Constructor
//
+ fJetEtaMax = 0.5;
+ fJetEtaMin = -0.5;
}
////////////////////////////////////////////////////////////////////////
//
// Constructor
//
+ fJetEtaMax = 0.5;
+ fJetEtaMin = -0.5;
}
////////////////////////////////////////////////////////////////////////
virtual ~AliJetHeader() { }
// Getters
- virtual TString GetComment() {return fComment;}
-
+ virtual TString GetComment() const {return fComment;}
+ virtual Float_t GetJetEtaMax() const {return fJetEtaMax;}
+ virtual Float_t GetJetEtaMin() const {return fJetEtaMin;}
+
// Setters
- virtual void SetComment(const char* com) {fComment=TString(com);}
+ virtual void SetComment(const char* com) {fComment=TString(com);}
+ virtual void SetJetEtaMax(Float_t eta= 0.5) {fJetEtaMax = eta;}
+ virtual void SetJetEtaMin(Float_t eta= -0.5) {fJetEtaMin = eta;}
+
// others
protected:
TString fComment; // a comment
-
+ Float_t fJetEtaMax; // maximum eta for the jet
+ Float_t fJetEtaMin; // minimum eta for the jet
+
ClassDef(AliJetHeader,1)
};
* provided "as is" without express or implied warranty. *
**************************************************************************/
-//
+//-------------------------------------------------------------------------
// Jet kinematics reader
// MC reader for jet analysis
-// Author: Andreas Morsch
-// andreas.morsch@cern.ch
-//
+// Author: Andreas Morsch (andreas.morsch@cern.ch)
+//-------------------------------------------------------------------------
// From root ...
#include <TClonesArray.h>
#include <TPDGCode.h>
#include <TParticle.h>
#include <TParticlePDG.h>
-#include <TVector3.h>
+//#include <TVector3.h>
#include <TLorentzVector.h>
#include <TSystem.h>
#include <TDatabasePDG.h>
void AliJetKineReader::OpenInputFiles()
{
- // Opens the input file using the run loader
- const char* dirName = fReaderHeader->GetDirectory();
- char path[256];
- sprintf(path, "%s/galice.root",dirName);
- fRunLoader = AliRunLoader::Open(path);
- fRunLoader->LoadKinematics();
- fRunLoader->LoadHeader();
-
- Int_t nMax = fRunLoader->GetNumberOfEvents();
- printf("\nTotal number of events = %d", nMax);
-
+ // Opens the input file using the run loader
+ const char* dirName = fReaderHeader->GetDirectory();
+ char path[256];
+ sprintf(path, "%s/galice.root",dirName);
+ fRunLoader = AliRunLoader::Open(path);
+ fRunLoader->LoadKinematics();
+ fRunLoader->LoadHeader();
+
+ // get number of events
+ Int_t nMax = fRunLoader->GetNumberOfEvents();
+ printf("\nTotal number of events = %d", nMax);
+
// set number of events in header
- if (fReaderHeader->GetLastEvent() == -1)
- fReaderHeader->SetLastEvent(nMax);
- else {
- Int_t nUsr = fReaderHeader->GetLastEvent();
- fReaderHeader->SetLastEvent(TMath::Min(nMax, nUsr));
- }
+ if (fReaderHeader->GetLastEvent() == -1)
+ fReaderHeader->SetLastEvent(nMax);
+ else {
+ Int_t nUsr = fReaderHeader->GetLastEvent();
+ fReaderHeader->SetLastEvent(TMath::Min(nMax, nUsr));
+ }
}
//____________________________________________________________________________
void AliJetKineReader::FillMomentumArray(Int_t event)
{
-//
-// Fill momentum array for event
-//
- Int_t goodTrack = 0;
- // Clear array
- ClearArray();
- // Get event from runloader
- fRunLoader->GetEvent(event);
- // Get the stack
- AliStack* stack = fRunLoader->Stack();
- // Number of primaries
- Int_t nt = stack->GetNprimary();
- // Get cuts set by user and header
- Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
- fAliHeader = fRunLoader->GetHeader();
-
- // Loop over particles
- Int_t* flag = new Int_t[nt];
- for (Int_t it = 0; it < nt; it++) {
- TParticle *part = stack->Particle(it);
- Int_t status = part->GetStatusCode();
- Int_t pdg = TMath::Abs(part->GetPdgCode());
- Float_t pt = part->Pt();
-
- // Skip non-final state particles, neutrinos and particles with pt < pt_min
-
- if (
- (status != 1)
- || (pdg == 12 || pdg == 14 || pdg == 16)
- || (pt < ptMin)
- ) continue;
-
-
- Float_t p = part->P();
- Float_t p0 = p;
- Float_t eta = part->Eta();
- Float_t phi = part->Phi();
-
- // Fast simulation of TPC if requested
- if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
- // Charged particles only
- Float_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
- if (charge == 0) continue;
- // Simulate efficiency
- if (!Efficiency(p0, eta, phi)) continue;
- // Simulate resolution
- p = SmearMomentum(4, p0);
- } // Fast TPC
-
- // Fast simulation of EMCAL if requested
- if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
- Float_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
- // Charged particles
- if (charge != 0) {
- // Simulate efficiency
- if (!Efficiency(p0, eta, phi)) continue;
- // Simulate resolution
- p = SmearMomentum(4, p0);
- } // charged
- // Neutral particles
- // Exclude K0L, n, nbar
- if (pdg == kNeutron || pdg == kK0Long) continue;
- } // Fast EMCAL
-
- fMass = part->GetCalcMass();
- fPdgC = part->GetPdgCode();
- // Fill momentum array
- Float_t r = p/p0;
-// printf("smeared %13.3f %13.3f %13.3f\n", p0, p, r);
-
- Float_t px = r * part->Px();
- Float_t py = r * part->Py();
- Float_t pz = r * part->Pz();
- Float_t m = part->GetMass();
- Float_t e = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
-
- new ((*fMomentumArray)[goodTrack])
- TLorentzVector(px, py, pz, e);
- goodTrack++;
+ // Fill momentum array for event
+
+ Int_t goodTrack = 0;
+ // Clear array
+ ClearArray();
+ // Get event from runloader
+ fRunLoader->GetEvent(event);
+ // Get the stack
+ AliStack* stack = fRunLoader->Stack();
+ // Number of primaries
+ Int_t nt = stack->GetNprimary();
+
+ // Get cuts set by user and header
+ Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
+ Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
+ Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
+ fAliHeader = fRunLoader->GetHeader();
+
+ // temporary storage of signal and cut flags
+ Int_t* sflag = new Int_t[nt];
+ Int_t* cflag = new Int_t[nt];
+
+ TLorentzVector p4;
+ // Loop over particles
+ for (Int_t it = 0; it < nt; it++) {
+ TParticle *part = stack->Particle(it);
+ Int_t status = part->GetStatusCode();
+ Int_t pdg = TMath::Abs(part->GetPdgCode());
+ Float_t pt = part->Pt();
+
+ // Skip non-final state particles, neutrinos
+ // and particles with pt < pt_min
+ if ((status != 1)
+ || (pdg == 12 || pdg == 14 || pdg == 16)) continue;
+
+ Float_t p = part->P();
+ Float_t p0 = p;
+ Float_t eta = part->Eta();
+ Float_t phi = part->Phi();
+
+ // Fast simulation of TPC if requested
+ if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
+ // Charged particles only
+ Float_t charge =
+ TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
+ if (charge == 0) continue;
+ // Simulate efficiency
+ if (!Efficiency(p0, eta, phi)) continue;
+ // Simulate resolution
+ p = SmearMomentum(4, p0);
+ } // End fast TPC
+
+ // Fast simulation of EMCAL if requested
+ if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
+ Float_t charge =
+ TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
+ // Charged particles only
+ if (charge != 0){
+ // Simulate efficiency
+ if (!Efficiency(p0, eta, phi)) continue;
+ // Simulate resolution
+ p = SmearMomentum(4, p0);
+ } // end "if" charged particles
+ // Neutral particles (exclude K0L, n, nbar)
+ if (pdg == kNeutron || pdg == kK0Long) continue;
+ } // End fast EMCAL
+
+ fMass = part->GetCalcMass();
+ fPdgC = part->GetPdgCode();
+
+ // Fill momentum array
+ Float_t r = p/p0;
+ Float_t px = r * part->Px();
+ Float_t py = r * part->Py();
+ Float_t pz = r * part->Pz();
+ Float_t m = part->GetMass();
+ Float_t e = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
+ p4 = TLorentzVector(px, py, pz, e);
+ if ( (p4.Eta()>etaMax) || (p4.Eta()<etaMin)) continue;
+ new ((*fMomentumArray)[goodTrack]) TLorentzVector(p4);
+
+ // all tracks are signal ... ?
+ sflag[goodTrack]=1;
+ cflag[goodTrack]=0;
+ if (pt > ptMin) cflag[goodTrack]=1; // track surviving pt cut
+ goodTrack++;
}
- fSignalFlag.Set(goodTrack,flag);
- printf("\nNumber of good tracks in event %d= %d \n",event,goodTrack);
+
+ // set the signal flags
+ fSignalFlag.Set(goodTrack,sflag);
+ fCutFlag.Set(goodTrack,cflag);
+ printf("\nIn event %d, number of good tracks %d \n", event, goodTrack);
+
}
Float_t AliJetKineReader::SmearMomentum(Int_t ind, Float_t p)
{
-//
-// The relative momentum error, i.e.
-// (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2,
-// where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity
-// (the lower value is for dn/d(eta) about 2000, and the higher one for 8000)
-//
-// If we include information from TRD b will be by a factor 2/3 smaller.
-//
-// ind = 1: high multiplicity
-// ind = 2: low multiplicity
-// ind = 3: high multiplicity + TRD
-// ind = 4: low multiplicity + TRD
-
- Float_t pSmeared;
- Float_t a = 0.75;
- Float_t b = 0.12;
-
- if (ind == 1) b = 0.12;
- if (ind == 2) b = 0.08;
- if (ind == 3) b = 0.12;
- if (ind == 4) b = 0.08;
-
- Float_t sigma = p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
- pSmeared = p + gRandom->Gaus(0., sigma);
- return pSmeared;
+ // The relative momentum error, i.e.
+ // (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2,
+ // where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity
+ // (the lower value is for dn/d(eta) about 2000, and the higher one for 8000)
+ //
+ // If we include information from TRD, b will be a factor 2/3 smaller.
+ //
+ // ind = 1: high multiplicity
+ // ind = 2: low multiplicity
+ // ind = 3: high multiplicity + TRD
+ // ind = 4: low multiplicity + TRD
+
+ Float_t pSmeared;
+ Float_t a = 0.75;
+ Float_t b = 0.12;
+
+ if (ind == 1) b = 0.12;
+ if (ind == 2) b = 0.08;
+ if (ind == 3) b = 0.12;
+ if (ind == 4) b = 0.08;
+
+ Float_t sigma = p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
+ pSmeared = p + gRandom->Gaus(0., sigma);
+ return pSmeared;
}
+
+
Bool_t AliJetKineReader::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi)
{
-//
-// Fast simulation of geometrical acceptance and tracking efficiency
-//
-// Tracking
- Float_t eff = 0.99;
- if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
-// Geometry
- if (p > 0.5) {
- phi *= 180. / TMath::Pi();
- // Sector number 0 - 17
- Int_t isec = Int_t(phi / 20.);
- // Sector centre
- Float_t phi0 = isec * 20. + 10.;
- Float_t phir = TMath::Abs(phi-phi0);
- // 2 deg of dead space
- if (phir > 9.) eff = 0.;
- }
- if (gRandom->Rndm() > eff) {
- return kFALSE;
- } else {
- return kTRUE;
- }
+ // Fast simulation of geometrical acceptance and tracking efficiency
+
+ // Tracking
+ Float_t eff = 0.99;
+ if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
+ // Geometry
+ if (p > 0.5) {
+ phi *= 180. / TMath::Pi();
+ // Sector number 0 - 17
+ Int_t isec = Int_t(phi / 20.);
+ // Sector centre
+ Float_t phi0 = isec * 20. + 10.;
+ Float_t phir = TMath::Abs(phi-phi0);
+ // 2 deg of dead space
+ if (phir > 9.) eff = 0.;
+ }
+
+ if (gRandom->Rndm() > eff) {
+ return kFALSE;
+ } else {
+ return kTRUE;
+ }
+
}
// Getters
Float_t GetParticleMass() const {return fMass;} // returns mass of the Track
Int_t GetParticlePdgCode() const {return fPdgC;} // returns Pdg code
-
// Setters
void FillMomentumArray(Int_t event);
void OpenInputFiles();
* provided "as is" without express or implied warranty. *
**************************************************************************/
+//---------------------------------------------------------------------
// Jet Kine Reader Header
// Header for the MC kinematics reader in the jet analysis
// Author: Andreas Morsch (andreas.morsch@cern.ch)
+//---------------------------------------------------------------------
#include "AliJetKineReaderHeader.h"
AliJetReaderHeader("AliJetKineReaderHeader")
{
// Constructor
- SetPtCut();
SetFastSimTPC(kFALSE);
+ SetFastSimEMCAL(kFALSE);
}
//____________________________________________________________________________
virtual ~AliJetKineReaderHeader();
// Setters
- void SetFastSimTPC(Bool_t flag = kTRUE) {fFastSimTPC = flag;}
- void SetFastSimEMCAL(Bool_t flag = kTRUE) {fFastSimEMCAL = flag;}
- virtual void SetPtCut(Float_t par = 2.0) {fPtCut = par;}
+ void SetFastSimTPC(Bool_t flag = kTRUE) {fFastSimTPC = flag;} // if TPC fast simulation
+ void SetFastSimEMCAL(Bool_t flag = kTRUE) {fFastSimEMCAL = flag;} // if EMCAL fast simulation
+
// Getter
Bool_t FastSimTPC() const {return fFastSimTPC;}
Bool_t FastSimEMCAL() const {return fFastSimEMCAL;}
- Float_t GetPtCut() const {return fPtCut;}
+
protected:
//parameters set by user
Bool_t fFastSimTPC;
Bool_t fFastSimEMCAL;
- Float_t fPtCut;
+
ClassDef(AliJetKineReaderHeader,1);
};
void GetPtHardLimits(Int_t bin, Float_t& ptmin, Float_t& ptmax);
TString GetRunTitle(Int_t bin);
Float_t GetWeight(Int_t bin);
- public:
+
Int_t fNbins; // Number of pt_hard bins used in the production
Float_t* fPtHardLimits; //[fNbins+1]
Float_t* fWeights; //[fNbins]
* provided "as is" without express or implied warranty. *
**************************************************************************/
-//---------------------------------------------------------------------
+//------------------------------------------------------------------------
// Jet reader base class
// manages the reading of input for jet algorithms
// Author: jgcn@mda.cinvestav.mx
-//---------------------------------------------------------------------
+//------------------------------------------------------------------------
+// root
#include <TClonesArray.h>
-
+//AliRoot
#include "AliJetReader.h"
#include "AliJetReaderHeader.h"
#include "AliESD.h"
fArrayMC = 0;
fAliHeader = 0;
fSignalFlag = TArrayI();
+ fCutFlag = TArrayI();
}
////////////////////////////////////////////////////////////////////////
/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* See cxx source for full Copyright notice */
-//---------------------------------------------------------------------
// Jet reader base class
// manages the reading of input for jet algorithms
// Author: jgcn@mda.cinvestav.mx
-//---------------------------------------------------------------------
#include <TObject.h>
#include <TChain.h>
class AliESD;
class AliHeader;
-
class AliJetReader : public TObject
{
public:
virtual TClonesArray *GetMomentumArray() {return fMomentumArray;}
virtual Int_t GetChainEntries() {return fChain->GetEntries();}
virtual AliJetReaderHeader* GetReaderHeader() { return fReaderHeader;}
- virtual AliHeader* GetAliHeader() { return fAliHeader;}
+ virtual AliHeader* GetAliHeader() {return fAliHeader;}
virtual Int_t GetSignalFlag(Int_t i) const {return fSignalFlag[i];}
+ virtual Int_t GetCutFlag(Int_t i) const {return fCutFlag[i];}
// Setters
virtual void FillMomentumArray(Int_t) {}
AliJetReaderHeader *fReaderHeader; // pointer to header
AliHeader *fAliHeader; // pointer to event header
TArrayI fSignalFlag; // to flag if a particle comes from pythia or
- // from the underlying event
+ // from the underlying event
+ TArrayI fCutFlag; // to flag if a particle passed the pt cut or not
ClassDef(AliJetReader,1)
};
* provided "as is" without express or implied warranty. *
**************************************************************************/
-//---------------------------------------------------------------------
+//-------------------------------------------------------------------------
+//-------------------------------------------------------------------------
// base class for Jet Reader Header
// Author: jgcn@mda.cinvestav.mx
-//---------------------------------------------------------------------
+//-------------------------------------------------------------------------
-
#include "AliJetReaderHeader.h"
ClassImp(AliJetReaderHeader)
TNamed("AliJetReaderHeader", "Jet Reader Header"),
fComment("No comment"),fDir(""),fPattern("")
{
- //
// Constructor
- //
fFirst = 0;
fLast = -1;
- fFiducialEtaMin = -0.5;
- fFiducialEtaMax = 0.5;
+ fFiducialEtaMin = -0.9;
+ fFiducialEtaMax = 0.9;
+ SetPtCut();
}
////////////////////////////////////////////////////////////////////////
TNamed(name, "Jet Reader Header"),
fComment("No comment"),fDir(""),fPattern("")
{
- //
// Constructor
- //
fFirst = 0;
fLast = -1;
- fFiducialEtaMin = -0.5;
- fFiducialEtaMax = 0.5;
+ fFiducialEtaMin = -0.9;
+ fFiducialEtaMax = 0.9;
+ SetPtCut();
}
////////////////////////////////////////////////////////////////////////
AliJetReaderHeader::~AliJetReaderHeader()
{
// destructor
- // delete fComment;
- //delete fDir;
- //delete fPattern;
+
}
/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* See cxx source for full Copyright notice */
-//---------------------------------------------------------------------
+//-------------------------------------------------------------------------
// base class for Jet Reader Header
// Author: jgcn@mda.cinvestav.mx
-//---------------------------------------------------------------------
+//-------------------------------------------------------------------------
#include <TNamed.h>
#include <TString.h>
virtual ~AliJetReaderHeader();
// Getters
- virtual TString GetComment() {return fComment;}
+ virtual TString GetComment() const {return fComment;}
virtual const char* GetDirectory() {return fDir.Data();}
virtual const char* GetPattern() {return fPattern.Data();}
virtual Float_t GetFiducialEtaMin() const {return fFiducialEtaMin;}
virtual Float_t GetFiducialEtaMax() const {return fFiducialEtaMax;}
-
+ virtual Float_t GetPtCut() const {return fPtCut;}
+ Float_t GetDCA() const {return fDCA;} // not working so far..(always 0)
+ Float_t GetTLength() const {return fTLength;} // not working so far.. (always 0)
Int_t GetNEvents() const {return fLast-fFirst;}
Int_t GetLastEvent() const {return fLast;}
Int_t GetFirstEvent() const {return fFirst;}
virtual void SetFirstEvent(Int_t i=0) {fFirst=i;}
virtual void SetLastEvent(Int_t i=-1) {fLast=i;}
virtual void SetFiducialEta(Float_t etamin, Float_t etamax)
- { fFiducialEtaMin = etamin; fFiducialEtaMax = etamax;}
+ { fFiducialEtaMin = etamin; fFiducialEtaMax = etamax;}
+ virtual void SetPtCut(Float_t par = 2.0) {fPtCut = par;}
+ virtual void SetDCA(Float_t dca = 0.0) {fDCA = dca;}
+ virtual void SetTLength(Float_t length = 0.0) {fTLength = length;}
+
protected:
Int_t fFirst; // First and last events analyzed
Int_t fLast; // in current set of files
Float_t fFiducialEtaMin; // Fiducial minimum eta
- Float_t fFiducialEtaMax; // Fiducial maximum eta
+ Float_t fFiducialEtaMax; // Fiducial maximum eta
+ Float_t fPtCut; // pt cut
+ Float_t fDCA; // dca cut
+ Float_t fTLength; // track length cut
TString fComment; // a comment
TString fDir; // directory with input files
TString fPattern; // pattern to look for input files
- ClassDef(AliJetReaderHeader,1);
+ ClassDef(AliJetReaderHeader,2);
};
#endif
AliLeading::AliLeading()
{
- //
// Constructor
- //
fNassoc = 0;
fLeading = new TLorentzVector(0.,0.,0.,0.);
fLow = -TMath::Pi()/2.0;
AliLeading::~AliLeading()
{
- //
// Destructor
- //
delete fLeading;
}
////////////////////////////////////////////////////////////////////////
void AliLeading::FindLeading(AliJetReader *reader)
-
{
- //
// find leading particle in the array of lorentz vectors
// lvArray and fill the correlation histogram
- //
-
-
+
AliJetReaderHeader* header = reader->GetReaderHeader();
-
+
TClonesArray* lvArray = reader->GetMomentumArray();
Int_t nIn = lvArray->GetEntries();
- fNassoc = nIn-1;
-
- if (fNassoc < 0) return;
-
+
// find max
Double_t ptMax = 0.0;
Int_t idxMax = -1;
for (Int_t i = 0; i < nIn; i++){
- TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
- if (lv->Pt() > ptMax &&
- lv->Eta() > header->GetFiducialEtaMin() &&
- lv->Eta() < header->GetFiducialEtaMax())
- {
- ptMax = lv->Pt();
- idxMax = i;
- }
+ TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
+ if ((reader->GetCutFlag(i) == 1) &&
+ lv->Pt() > ptMax &&
+ lv->Eta() > header->GetFiducialEtaMin() &&
+ lv->Eta() < header->GetFiducialEtaMax()){
+ ptMax = lv->Pt();
+ idxMax = i;
+ }
}
if (idxMax == -1) {
- fFound = kFALSE;
- Reset();
- return;
+ fFound = kFALSE;
+ Reset();
+ return;
}
// fill correlation array
fLeading = (TLorentzVector*) lvArray->At(idxMax);
fFound = kTRUE;
+ fNassoc = 0;
for (Int_t i = 0; i < nIn; i++) {
- if (i == idxMax) continue;
- TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
+ if (i == idxMax) continue;
+ TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
+ if ( (reader->GetCutFlag(i) == 1) &&
+ lv->Eta() > header->GetFiducialEtaMin() &&
+ lv->Eta() < header->GetFiducialEtaMax()) {
Double_t dphi = fLeading->DeltaPhi(*lv);
if (dphi < fLow) dphi = 2.0 * TMath::Pi() + dphi;
// find bin and fill array
-
Int_t iBin = (Int_t)
- TMath::Floor((dphi - fLow)
- *((Double_t) fnBin) / (2.0 * TMath::Pi()));
+ TMath::Floor((dphi - fLow)
+ *((Double_t) fnBin) / (2.0 * TMath::Pi()));
fCorr.AddAt(fCorr.At(iBin)+1,iBin);
+ fNassoc++;
+ }
}
}
////////////////////////////////////////////////////////////////////////
void AliLeading::PrintLeading()
-
{
// Print leading particle information
if (fNassoc<0) {
<< fLeading->Eta() << "," << fLeading->Phi() << ")" << endl;
cout << " " << fNassoc << " associated particles." << endl;
}
+
+////////////////////////////////////////////////////////////////////////
+
+Double_t AliLeading::GetE()
+{
+ return fLeading->E();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+Double_t AliLeading::GetPt()
+{
+ return fLeading->Pt();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+Double_t AliLeading::GetEta()
+{
+ return fLeading->Eta();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+Double_t AliLeading::GetPhi()
+{
+ // get phi of leading
+ return ( (fLeading->Phi() < 0) ?
+ (fLeading->Phi()) + 2. * TMath::Pi() :
+ fLeading->Phi());
+}
+
/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* See cxx source for full Copyright notice */
-
//---------------------------------------------------------------------
// Class to find and store the leading particle in event and
// store its correlation to associated particles
class TLorentzVector;
class AliJetReader;
-
class AliLeading : public TObject
{
public:
~AliLeading();
// Getters
- Int_t GetNassoc() const { return fNassoc;}
- TLorentzVector* GetLeading() const { return fLeading;}
- TArrayI GetCorr() const { return fCorr;}
- Bool_t LeadingFound() const { return fFound;}
-
+ Int_t GetNassoc() const {return fNassoc;}
+ TLorentzVector* GetLeading() const {return fLeading;}
+ TArrayI GetCorr() const {return fCorr;}
+ Bool_t LeadingFound() const {return fFound;}
+ Double_t GetE();
+ Double_t GetPt();
+ Double_t GetEta();
+ Double_t GetPhi();
+
// Setters
// Others
void FindLeading(AliJetReader* reader);
protected:
- Int_t fNassoc; // number of associated particles
+ Int_t fNassoc; // number of associated particles
TLorentzVector* fLeading; // leading particle
- TArrayI fCorr; // array to store azimuthal correlation
- // between leading and assoc particles
- Int_t fnBin; // number of bins in array
- Double_t fLow; // value corresponding to lower bound of bin 0
- Bool_t fFound;
+ TArrayI fCorr; // array to store azimuthal correlation
+ // between leading and assoc particles
+ Int_t fnBin; // number of bins in array
+ Double_t fLow; // value corresponding to lower bound of bin 0
+ Bool_t fFound; // leading found
ClassDef(AliLeading,1);
};
#include <Riostream.h>
#include <TLorentzVector.h>
+#include <TFile.h>
#include "AliPxconeJetFinder.h"
#include "AliPxconeJetHeader.h"
#include "AliJetReader.h"
AliPxconeJetHeader::AliPxconeJetHeader():
AliJetHeader("AliPxconeJetHeader")
{
- //
// Constructor
- //
SetMode();
SetRadius();
SetMinPt();
////////////////////////////////////////////////////////////////////////
-void AliPxconeJetHeader::PrintParameters()
+void AliPxconeJetHeader::PrintParameters() const
{
- //
// prints out parameters of jet algorithm
- //
-
cout << " PXCONE jet algorithm " << endl;
cout << " Running mode: " << fMode << endl;
cout << " Cone size: " << fRadius << endl;
void SetOverlap(Double_t o=.75) {fOverlap=o;}
// others
- void PrintParameters();
+ void PrintParameters() const;
protected:
Int_t fMode; // ee or pp mode
* about the suitability of this software for any purpose. It is *
* provided "as is" without express or implied warranty. *
**************************************************************************/
-
-
//---------------------------------------------------------------------
// UA1 Jet finder
#include <Riostream.h>
#include <TLorentzVector.h>
+#include <TFile.h>
#include <TH2F.h>
#include "AliUA1JetFinder.h"
#include "AliUA1JetHeader.h"
AliUA1JetFinder::AliUA1JetFinder()
{
- //
// Constructor
- //
fHeader = 0x0;
fLego = 0x0;
}
AliUA1JetFinder::~AliUA1JetFinder()
{
- //
// destructor
- //
-
- // delete fLego;
-
- // reset and delete header
}
////////////////////////////////////////////////////////////////////////
{
// initialize event, look for jets, download jet info
-
// initialization in 2 steps
// 1) transform input to pt,eta,phi plus lego
// 2) dump lego
if (nIn == 0) return;
// local arrays for input
+ Float_t* enT = new Float_t[nIn];
Float_t* ptT = new Float_t[nIn];
Float_t* etaT = new Float_t[nIn];
Float_t* phiT = new Float_t[nIn];
// load input vectors
for (Int_t i = 0; i < nIn; i++){
- TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
- ptT[i] = lv->Pt();
- etaT[i] = lv->Eta();
- phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
+ TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
+ enT[i] = lv->Energy();
+ ptT[i] = lv->Pt();
+ etaT[i] = lv->Eta();
+ phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
+ if (fReader->GetCutFlag(i) == 1)
fLego->Fill(etaT[i], phiT[i], ptT[i]);
}
fJets->SetNinput(nIn);
TAxis* xaxis = fLego->GetXaxis();
TAxis* yaxis = fLego->GetYaxis();
Float_t e = 0.0;
- for (Int_t i = 1; i <= fHeader->GetNbinEta(); i++) {
- for (Int_t j = 1; j <= fHeader->GetNbinPhi(); j++) {
+ for (Int_t i = 1; i <= fHeader->GetLegoNbinEta(); i++) {
+ for (Int_t j = 1; j <= fHeader->GetLegoNbinPhi(); j++) {
e = fLego->GetBinContent(i,j);
if (e > 0.0) e -= fHeader->GetMinCellEt();
if (e < 0.0) e = 0.;
}
// run the algo. Parameters from header
- Int_t nTot = (fHeader->GetNbinEta())*(fHeader->GetNbinPhi());
+ Int_t nTot = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
Float_t minmove = fHeader->GetMinMove();
Float_t maxmove = fHeader->GetMaxMove();
Int_t mode = fHeader->GetMode();
ua1_jet_finder(nCell, nTot, etCell, etaCell, phiCell,
minmove, maxmove, mode, precbg, ierror);
-
// sort jets
Int_t * idx = new Int_t[UA1JETS.njet];
TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
// download jet info.
- AliJetReaderHeader* rheader = fReader->GetReaderHeader();
for(Int_t i = 0; i < UA1JETS.njet; i++) {
// reject events outside acceptable eta range
- if (((UA1JETS.etaj[1][idx[i]])> (rheader->GetFiducialEtaMax()))
- || ((UA1JETS.etaj[1][idx[i]]) < (rheader->GetFiducialEtaMin())))
- continue;
- Float_t px, py,pz,en; // convert to 4-vector
- px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
- py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
- pz = UA1JETS.etj[idx[i]] /
- TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
- en = TMath::Sqrt(px * px + py * py + pz * pz);
- fJets->AddJet(px, py, pz, en);
+ if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
+ || ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
+ continue;
+ Float_t px, py,pz,en; // convert to 4-vector
+ px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
+ py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
+ pz = UA1JETS.etj[idx[i]] /
+ TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
+ en = TMath::Sqrt(px * px + py * py + pz * pz);
+ fJets->AddJet(px, py, pz, en);
}
// find multiplicities and relationship jet-particle
// find also percentage of pt from pythia
Int_t* injet = new Int_t[nIn];
- for (Int_t i = 0; i < nIn; i++) injet[i]= -1;
+ Int_t* sflag = new Int_t[nIn];
+ for (Int_t i = 0; i < nIn; i++) {injet[i]= 0;sflag[i]=0;}
Int_t* mult = new Int_t[fJets->GetNJets()];
+ Int_t* ncell = new Int_t[fJets->GetNJets()];
Float_t* percentage = new Float_t[fJets->GetNJets()];
for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
- Float_t pt_sig = 0.0;
- mult[i] = 0;
- for (Int_t j = 0; j < nIn; j++) {
- Float_t deta = etaT[j] - fJets->GetEta(i);
- Float_t dphi = phiT[j] - fJets->GetPhi(i);
- if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
- if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
- Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
- if (dr < fHeader->GetRadius() && injet[j] == -1) {
- if (fReader->GetSignalFlag(j) == 1) pt_sig+=ptT[j];
- injet[j] = i;
- mult[i]++;
+ Float_t pt_sig = 0.0;
+ mult[i] = 0;
+ ncell[i] = UA1JETS.ncellj[i];
+ for (Int_t j = 0; j < nIn; j++) {
+ Float_t deta = etaT[j] - fJets->GetEta(i);
+ Float_t dphi = phiT[j] - fJets->GetPhi(i);
+ if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
+ if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
+ Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
+ if (dr < fHeader->GetRadius() && injet[j] == 0) {
+ injet[j] = -(i+1);
+ if ((fReader->GetCutFlag(j)) == 1 &&
+ (etaT[j] < fHeader->GetLegoEtaMax()) &&
+ (etaT[j] > fHeader->GetLegoEtaMin())) {
+ injet[j] = i+1;
+ mult[i]++;
+ if (fReader->GetSignalFlag(j) == 1) {
+ pt_sig+=ptT[j];
+ sflag[j]=1;
+ }
+ }
+ }
}
- }
- percentage[i] = pt_sig/((Double_t) fJets->GetPt(i));
+ percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
+ ((Double_t) fJets->GetPt(i));
}
+ fJets->SetNCells(ncell);
fJets->SetPtFromSignal(percentage);
fJets->SetMultiplicities(mult);
fJets->SetInJet(injet);
+ fJets->SetEtaIn(etaT);
+ fJets->SetPhiIn(phiT);
+ fJets->SetPtIn(ptT);
+ fJets->SetEtAvg(UA1JETS.etavg);
+ delete ncell;
+ delete enT;
+ delete ptT;
+ delete etaT;
+ delete phiT;
+ delete injet;
+ delete idx;
+ delete mult;
+ delete percentage;
}
////////////////////////////////////////////////////////////////////////
void AliUA1JetFinder::Reset()
-
{
fLego->Reset();
fJets->ClearJets();
fHeader->Write();
}
-
////////////////////////////////////////////////////////////////////////
void AliUA1JetFinder::Init()
{
// initializes some variables
Float_t dEta, dPhi;
- dEta = (fHeader->GetEtaMax()-fHeader->GetEtaMin())
- /((Float_t) fHeader->GetNbinEta());
- dPhi = (fHeader->GetPhiMax()-fHeader->GetPhiMin())
- /((Float_t) fHeader->GetNbinPhi());
+ dEta = (fHeader->GetLegoEtaMax()-fHeader->GetLegoEtaMin())
+ /((Float_t) fHeader->GetLegoNbinEta());
+ dPhi = (fHeader->GetLegoPhiMax()-fHeader->GetLegoPhiMin())
+ /((Float_t) fHeader->GetLegoNbinPhi());
UA1CELL.etaCellSize = dEta;
UA1CELL.phiCellSize = dPhi;
// book lego
fLego = new
TH2F("legoH","eta-phi",
- fHeader->GetNbinEta(), fHeader->GetEtaMin(), fHeader->GetEtaMax(),
- fHeader->GetNbinPhi(), fHeader->GetPhiMin(), fHeader->GetPhiMax());
+ fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(),
+ fHeader->GetLegoEtaMax(), fHeader->GetLegoNbinPhi(),
+ fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
}
#include "AliUA1JetHeader.h"
ClassImp(AliUA1JetHeader)
-
////////////////////////////////////////////////////////////////////////
AliUA1JetHeader::AliUA1JetHeader():
AliJetHeader("AliUA1JetHeader")
{
- //
// Constructor
- //
fConeRadius = 0.3;
fEtSeed = 3.0;
fMinJetEt = 10.0;
fMinMove = 0.05;
fMaxMove = 0.15;
fPrecBg = 0.035;
- fNbinEta = 36;
- fNbinPhi = 124;
- fPhiMin = 0.;
- fPhiMax = 2. * TMath::Pi();
- fEtaMin = -0.9;
- fEtaMax = 0.9;
+ fLegoNbinEta = 36;
+ fLegoNbinPhi = 124;
+ fLegoPhiMin = 0.;
+ fLegoPhiMax = 2. * TMath::Pi();
+ fLegoEtaMin = -0.9;
+ fLegoEtaMax = 0.9;
+ fOnlySignal = kFALSE;
+ fOnlyBkgd = kFALSE;
}
-
////////////////////////////////////////////////////////////////////////
void AliUA1JetHeader::PrintParameters() const
-
{
- //
// prints out parameters of jet algorithm
- //
cout << " UA1 jet algorithm " << endl;
cout << " * Jet parameters: " << endl;
cout << " Maximum allowed move: " << fMaxMove << endl;
cout << " Precision for background: " << fPrecBg << endl;
cout << " * Lego parameters: " << endl;
- cout << " Number of bins in eta: " << fNbinEta<< endl;
- cout << " Number of bins in phi: " << fNbinPhi<< endl;
- cout << " Minimum azimuthal angle: " << fPhiMin<< endl;
- cout << " Maximum azimuthal angle: " << fPhiMax<< endl;
- cout << " Minimum rapidity angle: " << fEtaMin<< endl;
- cout << " Maximum rapidity angle: " << fEtaMax<< endl;
+ cout << " Number of bins in eta: " << fLegoNbinEta<< endl;
+ cout << " Number of bins in phi: " << fLegoNbinPhi<< endl;
+ cout << " Minimum azimuthal angle: " << fLegoPhiMin<< endl;
+ cout << " Maximum azimuthal angle: " << fLegoPhiMax<< endl;
+ cout << " Minimum rapidity angle: " << fLegoEtaMin<< endl;
+ cout << " Maximum rapidity angle: " << fLegoEtaMax<< endl;
}
/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* See cxx source for full Copyright notice */
-
//---------------------------------------------------------------------
// Jet header class for UA1 algorithm
// Stores the parameters of the UA1 jet algorithm
Float_t GetMinMove() const {return fMinMove;}
Float_t GetMaxMove() const {return fMaxMove;}
Float_t GetPrecBg() const {return fPrecBg;}
- Int_t GetNbinEta() const {return fNbinEta;}
- Int_t GetNbinPhi() const {return fNbinPhi;}
- Float_t GetEtaMin() const {return fEtaMin;}
- Float_t GetEtaMax() const {return fEtaMax;}
- Float_t GetPhiMin() const {return fPhiMin;}
- Float_t GetPhiMax() const {return fPhiMax;}
+ Int_t GetLegoNbinEta() const {return fLegoNbinEta;}
+ Int_t GetLegoNbinPhi() const {return fLegoNbinPhi;}
+ Float_t GetLegoEtaMin() const {return fLegoEtaMin;}
+ Float_t GetLegoEtaMax() const {return fLegoEtaMax;}
+ Float_t GetLegoPhiMin() const {return fLegoPhiMin;}
+ Float_t GetLegoPhiMax() const {return fLegoPhiMax;}
+ Bool_t GetOnlySignal() const {return fOnlySignal;}
+ Bool_t GetOnlyBkgd() const {return fOnlyBkgd;}
// Setters
void SetRadius(Float_t f) {fConeRadius=f;}
void SetMinMove(Float_t f) {fMinMove=f;}
void SetMaxMove(Float_t f) {fMaxMove=f;}
void SetPrecBg(Float_t f) {fPrecBg=f;}
- void SetNbinEta(Int_t f) {fNbinEta=f;}
- void SetNbinPhi(Int_t f) {fNbinPhi=f;}
- void SetEtaMin(Float_t f) {fEtaMin=f;}
- void SetEtaMax(Float_t f) {fEtaMax=f;}
- void SetPhiMin(Float_t f) {fPhiMin=f;}
- void SetPhiMax(Float_t f) {fPhiMax=f;}
+ void SetLegoNbinEta(Int_t f) {fLegoNbinEta=f;}
+ void SetLegoNbinPhi(Int_t f) {fLegoNbinPhi=f;}
+ void SetLegoEtaMin(Float_t f) {fLegoEtaMin=f;}
+ void SetLegoEtaMax(Float_t f) {fLegoEtaMax=f;}
+ void SetLegoPhiMin(Float_t f) {fLegoPhiMin=f;}
+ void SetLegoPhiMax(Float_t f) {fLegoPhiMax=f;}
+ void SetOnlySignal(Bool_t b) {fOnlySignal= b;}
+ void SetOnlyBkgd(Bool_t b) {fOnlyBkgd= b;}
// others
void PrintParameters() const;
Float_t fMaxMove; // max cone move
Float_t fPrecBg; // max value of change for BG (in %)
// parameters for legos
- Int_t fNbinEta; //! number of cells in eta
- Int_t fNbinPhi; //! number of cells in phi
- Float_t fEtaMin; //! minimum eta
- Float_t fEtaMax; //! maximum eta
- Float_t fPhiMin; //! minimun phi
- Float_t fPhiMax; //! maximum phi
+ Int_t fLegoNbinEta; //! number of cells in eta
+ Int_t fLegoNbinPhi; //! number of cells in phi
+ Float_t fLegoEtaMin; //! minimum eta
+ Float_t fLegoEtaMax; //! maximum eta
+ Float_t fLegoPhiMin; //! minimun phi
+ Float_t fLegoPhiMax; //! maximum phi
+ // parameters for Jet search
+ Bool_t fOnlySignal; // use only signal
+ Bool_t fOnlyBkgd; // use only background
ClassDef(AliUA1JetHeader,1)
};
#pragma link C++ class AliUA1JetFinder+;
#pragma link C++ class AliJetReaderHeader+;
#pragma link C++ class AliJetESDReaderHeader+;
-#pragma link C++ class AliJetMCReaderHeader+;
#pragma link C++ class AliJetReader+;
#pragma link C++ class AliJetESDReader+;
-#pragma link C++ class AliJetMCReader+;
#pragma link C++ class AliJetKineReader+;
#pragma link C++ class AliJetKineReaderHeader+;
#pragma link C++ class AliJetControlPlots+;
#pragma link C++ class AliJetProductionData+;
#pragma link C++ class AliJetProductionDataPDC2004+;
#pragma link C++ class AliJetAnalysis+;
+#pragma link C++ class AliJetDistributions+;
#endif
Float_t etaj[2][100];
Float_t phij[2][100];
Int_t ncellj[100];
+ Float_t etavg;
} UA1JetsCommon;
#define UA1JETS COMMON_BLOCK(UA1JETS,ua1jets)
AliJetFinder.cxx AliPxconeJetFinder.cxx AliJetReaderHeader.cxx \
AliJetESDReaderHeader.cxx AliJetReader.cxx AliJetESDReader.cxx \
AliJetControlPlots.cxx AliUA1JetHeader.cxx AliUA1JetFinder.cxx \
- AliJetMCReaderHeader.cxx AliJetMCReader.cxx AliLeading.cxx \
+ AliLeading.cxx \
AliJetKineReader.cxx AliJetKineReaderHeader.cxx \
AliJetProductionData.cxx AliJetProductionDataPDC2004.cxx \
- AliJetAnalysis.cxx
+ AliJetAnalysis.cxx AliJetDistributions.cxx
FSRCS = pxcone.F ua1_jet_finder.F
TH1F* dr1H = new TH1F("dr1H", "delta R", 160., 0., 2.);
TH1F* dr2H = new TH1F("dr2H", "delta R", 160., 0., 2.);
TH1F* dr3H = new TH1F("dr4H", "delta R", 160., 0., 2.);
- TH1F* etaH = new TH1F("etaH", "eta", 160., -2., 2.);
+ TH1F* etaH = new TH1F("etaH", "eta", 160., -2., 2.);
+ TH1F* eta1H = new TH1F("eta1H", "eta", 160., -2., 2.);
+ TH1F* eta2H = new TH1F("eta2H", "eta", 160., -2., 2.);
+
+ TH1F* phiH = new TH1F("phiH", "phi", 160., -3., 3.);
+ TH1F* phi1H = new TH1F("phi1H", "phi", 160., 0., 6.28);
+ TH1F* phi2H = new TH1F("phi2H", "phi", 160., 0., 6.28);
// load jet library
Float_t egen = gjets->GetPt(igen);
e1H->Fill(gjets->GetPt(igen));
-
- if (egen > 105. && egen < 125.) {
+ Float_t etag = gjets->GetEta(igen);
+ Float_t phig = gjets->GetPhi(igen);
+ Float_t dphi = phig - phij;
+
+ if (egen > 125. && egen < 150.) {
+ phiH->Fill((dphi));
+ etaH->Fill(etag - etaj);
+ phi1H->Fill(phig);
+ phi2H->Fill(phij);
+ eta1H->Fill(etag);
+ eta2H->Fill(etaj);
e4H->Fill(emax);
- if (rmin > 0.3) etaH->Fill(etaj);
dr2H->Fill(rmin);
}
}
Float_t deta = etag-etal;
Float_t dphi = TMath::Abs(phig - phil);
if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
+
Float_t r = TMath::Sqrt(deta * deta + dphi * dphi);
if (r < rmin) {
igen = j;
}
}
-
+ if (egen > 125. && egen < 150.)
dr3H->Fill(rmin);
// cout << " Generated Jets:" << endl;
dr2H->Draw("");
TCanvas* c3 = new TCanvas("c3");
- dr3H->Draw();
- dr2H->Draw("same");
+ dr2H->Draw();
+ dr3H->Draw("same");
TCanvas* c4 = new TCanvas("c4");
eH->Draw();
TCanvas* c5 = new TCanvas("c5");
etaH->Draw();
+ TCanvas* c5a = new TCanvas("c5a");
+ eta1H->Draw();
+
+ TCanvas* c5b = new TCanvas("c5b");
+ eta2H->Draw();
+
TCanvas* c6 = new TCanvas("c6");
e4H->Draw();
+ TCanvas* c7 = new TCanvas("c7");
+ phiH->Draw();
+
+ TCanvas* c7a = new TCanvas("c7a");
+ phi1H->Draw();
+ TCanvas* c7b = new TCanvas("c7b");
+ phi2H->Draw();
}
integer kpri
data kpri /0/
integer njet, ncellj
- real etj, etaj, phij
+ real etj, etaj, phij, etavg
* Results
COMMON /UA1JETS/ NJET, ETJ(100), ETAJ(100,2), PHIJ(100,2),
- + NCELLJ(100)
+ + NCELLJ(100), ETAVG
* Cell Geometry
COMMON /UA1CELL/ etaCellSize, phiCellSize
* Parameters
c*-reject cluster below minimum Ej_min
c* protection (am)
- etas=eta+etas/ets
- arg = 0.
- if (ets .ne. 0.) then
- if (abs(etas/ets) .lt. 23.719) then
- arg = ets * cosh(etas/ets)
- else
- arg = 1.e10
- endif
- endif
+
+c arg = 0.
+c if (ets .ne. 0.) then
+c if (abs(etas/ets) .lt. 23.719) then
+c arg = ets * cosh(etas/ets)
+c else
+c arg = 1.e10
+c endif
+c endif
- if(arg .lt. ej_min) then
+ if(ets .lt. ej_min) then
do k=1,ncell
if(flag(k).le.0) flag(k)=0
enddo
do k=1,ncell
if(flag(k).eq.-1) flag(k)=1
enddo
+ etas=eta+etas/ets
phi=phi+phis/ets
do while(phi .ge. C_2PI)
phi=phi-C_2PI
etaj(njet,2)=etas
phij(njet,2)=phi
ncellj(njet)=nc
+ etavg = et_ave
endif
endif
i=i+1
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff