From 83a444b167f04e7494dde6f6a71da3923497958a Mon Sep 17 00:00:00 2001 From: hristov Date: Thu, 19 Jan 2006 09:15:05 +0000 Subject: [PATCH] PPR version of the JETAN code (M. Lopez Noriega) --- JETAN/AliJet.cxx | 74 +- JETAN/AliJet.h | 24 +- JETAN/AliJetAnalysis.cxx | 1507 ++++++++++++++++++++++-------- JETAN/AliJetAnalysis.h | 236 ++++- JETAN/AliJetControlPlots.cxx | 161 +++- JETAN/AliJetControlPlots.h | 32 +- JETAN/AliJetESDReader.cxx | 69 +- JETAN/AliJetESDReader.h | 6 +- JETAN/AliJetESDReaderHeader.cxx | 13 +- JETAN/AliJetESDReaderHeader.h | 15 +- JETAN/AliJetFinder.cxx | 29 +- JETAN/AliJetFinder.h | 10 +- JETAN/AliJetHeader.cxx | 4 + JETAN/AliJetHeader.h | 15 +- JETAN/AliJetKineReader.cxx | 313 ++++--- JETAN/AliJetKineReader.h | 1 - JETAN/AliJetKineReaderHeader.cxx | 4 +- JETAN/AliJetKineReaderHeader.h | 10 +- JETAN/AliJetProductionData.h | 2 +- JETAN/AliJetReader.cxx | 8 +- JETAN/AliJetReader.h | 9 +- JETAN/AliJetReaderHeader.cxx | 24 +- JETAN/AliJetReaderHeader.h | 23 +- JETAN/AliLeading.cxx | 87 +- JETAN/AliLeading.h | 28 +- JETAN/AliPxconeJetFinder.cxx | 1 + JETAN/AliPxconeJetHeader.cxx | 7 +- JETAN/AliPxconeJetHeader.h | 2 +- JETAN/AliUA1JetFinder.cxx | 122 +-- JETAN/AliUA1JetHeader.cxx | 33 +- JETAN/AliUA1JetHeader.h | 44 +- JETAN/JetAnalysisLinkDef.h | 3 +- JETAN/UA1Common.h | 1 + JETAN/libJETAN.pkg | 4 +- JETAN/read_jets.C | 42 +- JETAN/ua1_jet_finder.F | 26 +- 36 files changed, 2066 insertions(+), 923 deletions(-) diff --git a/JETAN/AliJet.cxx b/JETAN/AliJet.cxx index bf6d17727a9..450f5a074d5 100644 --- a/JETAN/AliJet.cxx +++ b/JETAN/AliJet.cxx @@ -25,21 +25,22 @@ #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(); } //////////////////////////////////////////////////////////////////////// @@ -57,7 +58,7 @@ AliJet::~AliJet() 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; @@ -71,14 +72,13 @@ TLorentzVector* AliJet::GetJet(Int_t i) { // 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; @@ -87,11 +87,19 @@ Int_t AliJet::GetMultiplicity(Int_t i) //////////////////////////////////////////////////////////////////////// +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(); } @@ -102,7 +110,6 @@ Double_t AliJet::GetPy(Int_t i) { // Get Py component of jet i if (OutOfRange(i, "AliJet::GetPy:")) return -1e30; - TLorentzVector *lv = (TLorentzVector*) fJets->At(i); return lv->Py(); } @@ -113,7 +120,6 @@ Double_t AliJet::GetPz(Int_t i) { // Get Pz component of jet i if (OutOfRange(i, "AliJet::GetPz:")) return -1e30; - TLorentzVector *lv = (TLorentzVector*) fJets->At(i); return lv->Pz(); } @@ -124,7 +130,6 @@ Double_t AliJet::GetP(Int_t i) { // Get momentum of jet i if (OutOfRange(i, "AliJet::GetP:")) return -1e30; - TLorentzVector *lv = (TLorentzVector*) fJets->At(i); return lv->P(); } @@ -135,7 +140,6 @@ Double_t AliJet::GetE(Int_t i) { // Get energy of jet i if (OutOfRange(i, "AliJet::GetE:")) return -1e30; - TLorentzVector *lv = (TLorentzVector*) fJets->At(i); return lv->E(); } @@ -146,7 +150,6 @@ Double_t AliJet::GetPt(Int_t i) { // Get transverse momentum of jet i if (OutOfRange(i, "AliJet::GetPt:")) return -1e30; - TLorentzVector *lv = (TLorentzVector*) fJets->At(i); return lv->Pt(); } @@ -157,7 +160,6 @@ Double_t AliJet::GetEta(Int_t i) { // Get eta of jet i if (OutOfRange(i, "AliJet::GetEta:")) return -1e30; - TLorentzVector *lv = (TLorentzVector*) fJets->At(i); return lv->Eta(); } @@ -168,7 +170,6 @@ Double_t AliJet::GetPhi(Int_t i) { // 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()); } @@ -179,7 +180,6 @@ Double_t AliJet::GetTheta(Int_t i) { // Get theta of jet i if (OutOfRange(i, "AliJet::GetTheta:")) return -1e30; - TLorentzVector *lv = (TLorentzVector*) fJets->At(i); return lv->Theta(); } @@ -190,7 +190,6 @@ Double_t AliJet::GetMass(Int_t i) { // Get invariant mass of jet i if (OutOfRange(i, "AliJet::GetMass:")) return -1e30; - TLorentzVector *lv = (TLorentzVector*) fJets->At(i); return lv->M(); } @@ -209,12 +208,37 @@ void AliJet::AddJet(Double_t px, Double_t py, Double_t pz, Double_t e) 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 @@ -233,6 +257,13 @@ void AliJet::SetMultiplicities(Int_t* m) //////////////////////////////////////////////////////////////////////// +void AliJet::SetNCells(Int_t* n) +{ + if (fNJets>0) fNCells.Set(fNJets, n); +} + +//////////////////////////////////////////////////////////////////////// + void AliJet::ClearJets(Option_t *option) { // reset all values @@ -241,6 +272,11 @@ void AliJet::ClearJets(Option_t *option) fNJets=0; fMultiplicities.Set(0); fInJet.Set(0); + fPtFromSignal.Set(0); + fPhiIn.Set(0); + fEtaIn.Set(0); + fPtIn.Set(0); + fNCells.Set(0); } //////////////////////////////////////////////////////////////////////// diff --git a/JETAN/AliJet.h b/JETAN/AliJet.h index d9b4389b77f..c9c5da71a5c 100644 --- a/JETAN/AliJet.h +++ b/JETAN/AliJet.h @@ -31,10 +31,16 @@ class AliJet : public TObject 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); @@ -49,10 +55,15 @@ class AliJet : public TObject // 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=""); @@ -62,11 +73,18 @@ class AliJet : public TObject 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) }; diff --git a/JETAN/AliJetAnalysis.cxx b/JETAN/AliJetAnalysis.cxx index 2c7e0cca05e..0f9787d0bb9 100644 --- a/JETAN/AliJetAnalysis.cxx +++ b/JETAN/AliJetAnalysis.cxx @@ -15,440 +15,1177 @@ //--------------------------------------------------------------------- // 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 #include "AliJetAnalysis.h" ClassImp(AliJetAnalysis) - - -//////////////////////////////////////////////////////////////////////// + +// root +#include #include +#include #include -#include #include #include #include #include #include - +#include +// 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=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;iGetEta(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 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;iGetEta(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 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 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 (drBFill(xi,fPtJ,fWeight); + } + // --- Jt + if (drBFill(jt,fPtJ,fWeight); + } + // --- dE/dr + if (drBFill(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;iGetEta(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 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;iGetEta(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;iFill(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(); +} diff --git a/JETAN/AliJetAnalysis.h b/JETAN/AliJetAnalysis.h index cba11095c76..9c94dd4e383 100644 --- a/JETAN/AliJetAnalysis.h +++ b/JETAN/AliJetAnalysis.h @@ -3,34 +3,238 @@ /* 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 +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) }; diff --git a/JETAN/AliJetControlPlots.cxx b/JETAN/AliJetControlPlots.cxx index 6e37ef5b8bf..a9c1030d327 100755 --- a/JETAN/AliJetControlPlots.cxx +++ b/JETAN/AliJetControlPlots.cxx @@ -22,9 +22,6 @@ #include #include -#include -#include -#include #include #include #include "AliJetReader.h" @@ -37,17 +34,20 @@ ClassImp(AliJetControlPlots) 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"); @@ -57,85 +57,144 @@ AliJetControlPlots::AliJetControlPlots() 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;iFill(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;iAt(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;iGetEta(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)); } //////////////////////////////////////////////////////////////////////// @@ -146,8 +205,8 @@ void AliJetControlPlots::PlotHistos() 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"); @@ -155,21 +214,23 @@ void AliJetControlPlots::PlotHistos() 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(); } diff --git a/JETAN/AliJetControlPlots.h b/JETAN/AliJetControlPlots.h index ff372db74c6..d3b19defb5e 100755 --- a/JETAN/AliJetControlPlots.h +++ b/JETAN/AliJetControlPlots.h @@ -31,31 +31,39 @@ class AliJetControlPlots : public TObject // 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) }; diff --git a/JETAN/AliJetESDReader.cxx b/JETAN/AliJetESDReader.cxx index 4a2ce343831..a2988d26ee1 100755 --- a/JETAN/AliJetESDReader.cxx +++ b/JETAN/AliJetESDReader.cxx @@ -13,17 +13,14 @@ * 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 #include #include +#include #include "AliJetESDReader.h" #include "AliJetESDReaderHeader.h" #include "AliESD.h" @@ -34,7 +31,7 @@ ClassImp(AliJetESDReader) AliJetESDReader::AliJetESDReader() { // Constructor - + printf("\nIn reader constructor\n"); fReaderHeader = 0x0; fMass = 0; fSign = 0; @@ -45,14 +42,14 @@ AliJetESDReader::AliJetESDReader() 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"); @@ -69,8 +66,8 @@ void AliJetESDReader::OpenInputFiles() 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); } } @@ -95,50 +92,58 @@ void AliJetESDReader::OpenInputFiles() 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); } diff --git a/JETAN/AliJetESDReader.h b/JETAN/AliJetESDReader.h index 790c83a06f7..e425f9d7bcb 100755 --- a/JETAN/AliJetESDReader.h +++ b/JETAN/AliJetESDReader.h @@ -4,9 +4,11 @@ /* 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; @@ -19,8 +21,8 @@ class AliJetESDReader : public AliJetReader 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); diff --git a/JETAN/AliJetESDReaderHeader.cxx b/JETAN/AliJetESDReaderHeader.cxx index 41fb5aed4b1..bdd606d049d 100755 --- a/JETAN/AliJetESDReaderHeader.cxx +++ b/JETAN/AliJetESDReaderHeader.cxx @@ -12,13 +12,12 @@ * 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" @@ -29,10 +28,8 @@ AliJetESDReaderHeader::AliJetESDReaderHeader(): AliJetReaderHeader("AliJetESDReaderHeader") { // Constructor - SetPtCut(); - SetDCA(); - SetTLength(); SetReadSignalOnly(kFALSE); + SetReadBkgdOnly(kFALSE); } diff --git a/JETAN/AliJetESDReaderHeader.h b/JETAN/AliJetESDReaderHeader.h index 56ac1e6b002..a0c36b848b2 100755 --- a/JETAN/AliJetESDReaderHeader.h +++ b/JETAN/AliJetESDReaderHeader.h @@ -18,25 +18,18 @@ class AliJetESDReaderHeader : public AliJetReaderHeader 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 diff --git a/JETAN/AliJetFinder.cxx b/JETAN/AliJetFinder.cxx index 7b866dd2f97..37bfa564142 100644 --- a/JETAN/AliJetFinder.cxx +++ b/JETAN/AliJetFinder.cxx @@ -12,8 +12,7 @@ * 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 @@ -31,16 +30,12 @@ #include "AliLeading.h" #include "AliHeader.h" - ClassImp(AliJetFinder) -//////////////////////////////////////////////////////////////////////// AliJetFinder::AliJetFinder() { - // // Constructor - // fOut = 0x0; fJets = new AliJet(); fGenJets = new AliJet(); @@ -50,15 +45,11 @@ AliJetFinder::AliJetFinder() SetPlotMode(kFALSE); } - //////////////////////////////////////////////////////////////////////// AliJetFinder::~AliJetFinder() { - // // destructor - // - // here reset and delete jets fJets->ClearJets(); delete fJets; @@ -72,10 +63,8 @@ AliJetFinder::~AliJetFinder() delete fOut; // reset and delete control plots if (fPlots) delete fPlots; - // delete fLeading; } - //////////////////////////////////////////////////////////////////////// void AliJetFinder::SetOutputFile(const char *name) @@ -84,25 +73,22 @@ 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(); } @@ -111,7 +97,7 @@ void AliJetFinder::SetPlotMode(Bool_t b) 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); @@ -134,12 +120,11 @@ void AliJetFinder::WriteRHeaderToFile() 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(); @@ -175,7 +160,6 @@ void AliJetFinder::Run() WriteRHeaderToFile(); WriteJHeaderToFile(); } - // loop over events Int_t nFirst,nLast; nFirst = fReader->GetReaderHeader()->GetFirstEvent(); @@ -187,7 +171,7 @@ void AliJetFinder::Run() GetGenJets(); FindJets(); if (fOut) WriteJetsToFile(i); - if (fPlots) fPlots->FillHistos(fJets,fReader); + if (fPlots) fPlots->FillHistos(fJets); fLeading->Reset(); fGenJets->ClearJets(); Reset(); @@ -203,4 +187,3 @@ void AliJetFinder::Run() fOut->Close(); } } - diff --git a/JETAN/AliJetFinder.h b/JETAN/AliJetFinder.h index efefbcff44d..c3feee45740 100755 --- a/JETAN/AliJetFinder.h +++ b/JETAN/AliJetFinder.h @@ -4,7 +4,6 @@ /* 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 @@ -19,7 +18,6 @@ class AliJetReader; class AliJetControlPlots; class AliLeading; - class AliJetFinder : public TObject { public: @@ -28,9 +26,9 @@ class AliJetFinder : public TObject 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"); @@ -51,7 +49,7 @@ class AliJetFinder : public TObject 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 diff --git a/JETAN/AliJetHeader.cxx b/JETAN/AliJetHeader.cxx index 0c39befed80..bd50406b3c5 100755 --- a/JETAN/AliJetHeader.cxx +++ b/JETAN/AliJetHeader.cxx @@ -32,6 +32,8 @@ AliJetHeader::AliJetHeader(): // // Constructor // + fJetEtaMax = 0.5; + fJetEtaMin = -0.5; } //////////////////////////////////////////////////////////////////////// @@ -43,6 +45,8 @@ AliJetHeader::AliJetHeader(const char * name): // // Constructor // + fJetEtaMax = 0.5; + fJetEtaMin = -0.5; } //////////////////////////////////////////////////////////////////////// diff --git a/JETAN/AliJetHeader.h b/JETAN/AliJetHeader.h index d85de2f7460..55e45a5fca7 100755 --- a/JETAN/AliJetHeader.h +++ b/JETAN/AliJetHeader.h @@ -23,16 +23,23 @@ class AliJetHeader : public TNamed 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) }; diff --git a/JETAN/AliJetKineReader.cxx b/JETAN/AliJetKineReader.cxx index 2b96d2c301a..363650c33c0 100644 --- a/JETAN/AliJetKineReader.cxx +++ b/JETAN/AliJetKineReader.cxx @@ -13,19 +13,18 @@ * 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 #include #include #include -#include +//#include #include #include #include @@ -59,166 +58,180 @@ AliJetKineReader::~AliJetKineReader() 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() 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; + } + } diff --git a/JETAN/AliJetKineReader.h b/JETAN/AliJetKineReader.h index ee26dbbbf38..c507aca54e8 100644 --- a/JETAN/AliJetKineReader.h +++ b/JETAN/AliJetKineReader.h @@ -21,7 +21,6 @@ class AliJetKineReader : public AliJetReader // 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(); diff --git a/JETAN/AliJetKineReaderHeader.cxx b/JETAN/AliJetKineReaderHeader.cxx index 6d2682f1de4..e338e780fd5 100644 --- a/JETAN/AliJetKineReaderHeader.cxx +++ b/JETAN/AliJetKineReaderHeader.cxx @@ -13,9 +13,11 @@ * 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" @@ -28,8 +30,8 @@ AliJetKineReaderHeader::AliJetKineReaderHeader(): AliJetReaderHeader("AliJetKineReaderHeader") { // Constructor - SetPtCut(); SetFastSimTPC(kFALSE); + SetFastSimEMCAL(kFALSE); } //____________________________________________________________________________ diff --git a/JETAN/AliJetKineReaderHeader.h b/JETAN/AliJetKineReaderHeader.h index 49ef2e161e1..e190359f140 100644 --- a/JETAN/AliJetKineReaderHeader.h +++ b/JETAN/AliJetKineReaderHeader.h @@ -18,19 +18,19 @@ class AliJetKineReaderHeader : public AliJetReaderHeader 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); }; diff --git a/JETAN/AliJetProductionData.h b/JETAN/AliJetProductionData.h index ee58a84de3f..eb5f24d5162 100644 --- a/JETAN/AliJetProductionData.h +++ b/JETAN/AliJetProductionData.h @@ -23,7 +23,7 @@ class AliJetProductionData : public TObject 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] diff --git a/JETAN/AliJetReader.cxx b/JETAN/AliJetReader.cxx index 3835a8a81ca..3ea3e465f27 100755 --- a/JETAN/AliJetReader.cxx +++ b/JETAN/AliJetReader.cxx @@ -13,14 +13,15 @@ * 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 - +//AliRoot #include "AliJetReader.h" #include "AliJetReaderHeader.h" #include "AliESD.h" @@ -40,6 +41,7 @@ AliJetReader::AliJetReader() fArrayMC = 0; fAliHeader = 0; fSignalFlag = TArrayI(); + fCutFlag = TArrayI(); } //////////////////////////////////////////////////////////////////////// diff --git a/JETAN/AliJetReader.h b/JETAN/AliJetReader.h index cb8979d8bff..1c71365470b 100755 --- a/JETAN/AliJetReader.h +++ b/JETAN/AliJetReader.h @@ -4,11 +4,9 @@ /* 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 #include @@ -19,7 +17,6 @@ class AliJetReaderHeader; class AliESD; class AliHeader; - class AliJetReader : public TObject { public: @@ -30,8 +27,9 @@ class AliJetReader : public TObject 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) {} @@ -51,7 +49,8 @@ class AliJetReader : public TObject 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) }; diff --git a/JETAN/AliJetReaderHeader.cxx b/JETAN/AliJetReaderHeader.cxx index cf39a0ae09e..269c67a0622 100755 --- a/JETAN/AliJetReaderHeader.cxx +++ b/JETAN/AliJetReaderHeader.cxx @@ -13,12 +13,12 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -//--------------------------------------------------------------------- +//------------------------------------------------------------------------- +//------------------------------------------------------------------------- // base class for Jet Reader Header // Author: jgcn@mda.cinvestav.mx -//--------------------------------------------------------------------- +//------------------------------------------------------------------------- - #include "AliJetReaderHeader.h" ClassImp(AliJetReaderHeader) @@ -29,13 +29,12 @@ AliJetReaderHeader::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(); } //////////////////////////////////////////////////////////////////////// @@ -44,13 +43,12 @@ AliJetReaderHeader::AliJetReaderHeader(const char * name): 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(); } //////////////////////////////////////////////////////////////////////// @@ -58,7 +56,5 @@ AliJetReaderHeader::AliJetReaderHeader(const char * name): AliJetReaderHeader::~AliJetReaderHeader() { // destructor - // delete fComment; - //delete fDir; - //delete fPattern; + } diff --git a/JETAN/AliJetReaderHeader.h b/JETAN/AliJetReaderHeader.h index fd9a78060c9..7466b423a94 100755 --- a/JETAN/AliJetReaderHeader.h +++ b/JETAN/AliJetReaderHeader.h @@ -4,10 +4,10 @@ /* 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 #include @@ -21,12 +21,14 @@ class AliJetReaderHeader : public TNamed 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;} @@ -38,18 +40,25 @@ class AliJetReaderHeader : public TNamed 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 diff --git a/JETAN/AliLeading.cxx b/JETAN/AliLeading.cxx index c0d2e89ca1a..f7f07ec17a7 100644 --- a/JETAN/AliLeading.cxx +++ b/JETAN/AliLeading.cxx @@ -34,9 +34,7 @@ ClassImp(AliLeading) AliLeading::AliLeading() { - // // Constructor - // fNassoc = 0; fLeading = new TLorentzVector(0.,0.,0.,0.); fLow = -TMath::Pi()/2.0; @@ -49,66 +47,62 @@ AliLeading::AliLeading() 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++; + } } } @@ -126,7 +120,6 @@ void AliLeading::Reset() //////////////////////////////////////////////////////////////////////// void AliLeading::PrintLeading() - { // Print leading particle information if (fNassoc<0) { @@ -141,3 +134,35 @@ void AliLeading::PrintLeading() << 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()); +} + diff --git a/JETAN/AliLeading.h b/JETAN/AliLeading.h index b01c5b4d670..cef8512f668 100644 --- a/JETAN/AliLeading.h +++ b/JETAN/AliLeading.h @@ -4,7 +4,6 @@ /* 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 @@ -17,7 +16,6 @@ class TLorentzVector; class AliJetReader; - class AliLeading : public TObject { public: @@ -26,11 +24,15 @@ class AliLeading : public TObject ~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); @@ -39,13 +41,13 @@ class AliLeading : public TObject 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); }; diff --git a/JETAN/AliPxconeJetFinder.cxx b/JETAN/AliPxconeJetFinder.cxx index becf0072c29..4e92d8eaa10 100755 --- a/JETAN/AliPxconeJetFinder.cxx +++ b/JETAN/AliPxconeJetFinder.cxx @@ -23,6 +23,7 @@ #include #include +#include #include "AliPxconeJetFinder.h" #include "AliPxconeJetHeader.h" #include "AliJetReader.h" diff --git a/JETAN/AliPxconeJetHeader.cxx b/JETAN/AliPxconeJetHeader.cxx index eb996840f0a..32fb1f7aadd 100755 --- a/JETAN/AliPxconeJetHeader.cxx +++ b/JETAN/AliPxconeJetHeader.cxx @@ -29,9 +29,7 @@ ClassImp(AliPxconeJetHeader) AliPxconeJetHeader::AliPxconeJetHeader(): AliJetHeader("AliPxconeJetHeader") { - // // Constructor - // SetMode(); SetRadius(); SetMinPt(); @@ -41,13 +39,10 @@ AliPxconeJetHeader::AliPxconeJetHeader(): //////////////////////////////////////////////////////////////////////// -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; diff --git a/JETAN/AliPxconeJetHeader.h b/JETAN/AliPxconeJetHeader.h index 007b0e62185..1e1acb5470b 100755 --- a/JETAN/AliPxconeJetHeader.h +++ b/JETAN/AliPxconeJetHeader.h @@ -34,7 +34,7 @@ class AliPxconeJetHeader : public AliJetHeader void SetOverlap(Double_t o=.75) {fOverlap=o;} // others - void PrintParameters(); + void PrintParameters() const; protected: Int_t fMode; // ee or pp mode diff --git a/JETAN/AliUA1JetFinder.cxx b/JETAN/AliUA1JetFinder.cxx index 2e291deaba3..c9c1e953a00 100755 --- a/JETAN/AliUA1JetFinder.cxx +++ b/JETAN/AliUA1JetFinder.cxx @@ -12,8 +12,6 @@ * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ - - //--------------------------------------------------------------------- // UA1 Jet finder @@ -24,6 +22,7 @@ #include #include +#include #include #include "AliUA1JetFinder.h" #include "AliUA1JetHeader.h" @@ -40,9 +39,7 @@ ClassImp(AliUA1JetFinder) AliUA1JetFinder::AliUA1JetFinder() { - // // Constructor - // fHeader = 0x0; fLego = 0x0; } @@ -52,13 +49,7 @@ AliUA1JetFinder::AliUA1JetFinder() AliUA1JetFinder::~AliUA1JetFinder() { - // // destructor - // - - // delete fLego; - - // reset and delete header } //////////////////////////////////////////////////////////////////////// @@ -90,7 +81,6 @@ void AliUA1JetFinder::FindJets() { // initialize event, look for jets, download jet info - // initialization in 2 steps // 1) transform input to pt,eta,phi plus lego // 2) dump lego @@ -101,16 +91,19 @@ void AliUA1JetFinder::FindJets() 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); @@ -125,8 +118,8 @@ void AliUA1JetFinder::FindJets() 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.; @@ -140,7 +133,7 @@ void AliUA1JetFinder::FindJets() } // 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(); @@ -149,60 +142,83 @@ void AliUA1JetFinder::FindJets() 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(); @@ -216,17 +232,16 @@ void AliUA1JetFinder::WriteJHeaderToFile() 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; @@ -240,6 +255,7 @@ void AliUA1JetFinder::Init() // 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()); } diff --git a/JETAN/AliUA1JetHeader.cxx b/JETAN/AliUA1JetHeader.cxx index 03566ca7039..6d4aed8b0ef 100755 --- a/JETAN/AliUA1JetHeader.cxx +++ b/JETAN/AliUA1JetHeader.cxx @@ -24,15 +24,12 @@ #include "AliUA1JetHeader.h" ClassImp(AliUA1JetHeader) - //////////////////////////////////////////////////////////////////////// AliUA1JetHeader::AliUA1JetHeader(): AliJetHeader("AliUA1JetHeader") { - // // Constructor - // fConeRadius = 0.3; fEtSeed = 3.0; fMinJetEt = 10.0; @@ -41,23 +38,21 @@ AliUA1JetHeader::AliUA1JetHeader(): 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; @@ -71,10 +66,10 @@ void AliUA1JetHeader::PrintParameters() const 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; } diff --git a/JETAN/AliUA1JetHeader.h b/JETAN/AliUA1JetHeader.h index 62eb5aa3902..f22c963666e 100755 --- a/JETAN/AliUA1JetHeader.h +++ b/JETAN/AliUA1JetHeader.h @@ -4,7 +4,6 @@ /* 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 @@ -29,12 +28,14 @@ class AliUA1JetHeader : public AliJetHeader 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;} @@ -45,12 +46,14 @@ class AliUA1JetHeader : public AliJetHeader 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; @@ -68,12 +71,15 @@ protected: 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) }; diff --git a/JETAN/JetAnalysisLinkDef.h b/JETAN/JetAnalysisLinkDef.h index e575817e491..d8abac04ec5 100644 --- a/JETAN/JetAnalysisLinkDef.h +++ b/JETAN/JetAnalysisLinkDef.h @@ -13,10 +13,8 @@ #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+; @@ -24,4 +22,5 @@ #pragma link C++ class AliJetProductionData+; #pragma link C++ class AliJetProductionDataPDC2004+; #pragma link C++ class AliJetAnalysis+; +#pragma link C++ class AliJetDistributions+; #endif diff --git a/JETAN/UA1Common.h b/JETAN/UA1Common.h index e7f332d34fe..4e1e6115031 100755 --- a/JETAN/UA1Common.h +++ b/JETAN/UA1Common.h @@ -14,6 +14,7 @@ extern "C" { Float_t etaj[2][100]; Float_t phij[2][100]; Int_t ncellj[100]; + Float_t etavg; } UA1JetsCommon; #define UA1JETS COMMON_BLOCK(UA1JETS,ua1jets) diff --git a/JETAN/libJETAN.pkg b/JETAN/libJETAN.pkg index 80508afd2e6..6d0d1a0f8e8 100644 --- a/JETAN/libJETAN.pkg +++ b/JETAN/libJETAN.pkg @@ -4,10 +4,10 @@ SRCS = AliJet.cxx AliJetHeader.cxx AliPxconeJetHeader.cxx \ 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 diff --git a/JETAN/read_jets.C b/JETAN/read_jets.C index f7df5e7fbc9..bf11ac4023d 100644 --- a/JETAN/read_jets.C +++ b/JETAN/read_jets.C @@ -13,7 +13,13 @@ void read_jets(const char* fn = "jets.root") 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 @@ -97,10 +103,18 @@ void read_jets(const char* fn = "jets.root") 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); } } @@ -124,6 +138,7 @@ void read_jets(const char* fn = "jets.root") 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) { @@ -131,7 +146,7 @@ void read_jets(const char* fn = "jets.root") igen = j; } } - + if (egen > 125. && egen < 150.) dr3H->Fill(rmin); // cout << " Generated Jets:" << endl; @@ -156,8 +171,8 @@ void read_jets(const char* fn = "jets.root") dr2H->Draw(""); TCanvas* c3 = new TCanvas("c3"); - dr3H->Draw(); - dr2H->Draw("same"); + dr2H->Draw(); + dr3H->Draw("same"); TCanvas* c4 = new TCanvas("c4"); eH->Draw(); @@ -165,6 +180,19 @@ void read_jets(const char* fn = "jets.root") 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(); } diff --git a/JETAN/ua1_jet_finder.F b/JETAN/ua1_jet_finder.F index 1bf3c02b6aa..2b0a74fd3be 100755 --- a/JETAN/ua1_jet_finder.F +++ b/JETAN/ua1_jet_finder.F @@ -38,10 +38,10 @@ c:>------------------------------------------------------------------ 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 @@ -227,17 +227,17 @@ c*-sum up unused cells within required distance of given eta/phi 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 @@ -246,6 +246,7 @@ c*-eles, store flags and jet variables 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 @@ -260,6 +261,7 @@ c*-eles, store flags and jet variables etaj(njet,2)=etas phij(njet,2)=phi ncellj(njet)=nc + etavg = et_ave endif endif i=i+1 -- 2.43.0