+// **************************************
+// Task used for the correction of determiantion of reconstructed jet spectra
+// Compares input (gen) and output (rec) jets
+// *******************************************
+
+
/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* *
* provided "as is" without express or implied warranty. *
**************************************************************************/
-
#include <TROOT.h>
#include <TSystem.h>
fh1Xsec(0x0),
fh1Trials(0x0),
fh1PtHard(0x0),
- fh1PtHard_NoW(0x0),
- fh1PtHard_Trials(0x0),
+ fh1PtHardNoW(0x0),
+ fh1PtHardTrials(0x0),
fh1NGenJets(0x0),
fh1NRecJets(0x0),
fHistList(0x0) ,
for(int ij = 0;ij<kMaxJets;++ij){
fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] = fh2Frag[ij] = fh2FragLn[ij] = fh2PtGenDeltaPhi[ij] = fh2PtGenDeltaEta[ij] = 0;
- fh3PtRecGenHard[ij] = fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3GenEtaPhiPt_NoFound[ij] = fh3GenEtaPhiPt[ij] = 0;
+ fh3PtRecGenHard[ij] = fh3PtRecGenHardNoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPtNoGen[ij] =fh3GenEtaPhiPtNoFound[ij] = fh3GenEtaPhiPt[ij] = 0;
}
for(int ic = 0;ic < kMaxCorrelation;ic++){
fhnCorrelation[ic] = 0;
fh1Xsec(0x0),
fh1Trials(0x0),
fh1PtHard(0x0),
- fh1PtHard_NoW(0x0),
- fh1PtHard_Trials(0x0),
+ fh1PtHardNoW(0x0),
+ fh1PtHardTrials(0x0),
fh1NGenJets(0x0),
fh1NRecJets(0x0),
fHistList(0x0) ,
fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] = fh2Frag[ij] = fh2FragLn[ij] = fh2PtGenDeltaPhi[ij] = fh2PtGenDeltaEta[ij] = 0;
- fh3PtRecGenHard[ij] = fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3GenEtaPhiPt_NoFound[ij] = fh3GenEtaPhiPt[ij] = 0;
+ fh3PtRecGenHard[ij] = fh3PtRecGenHardNoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPtNoGen[ij] =fh3GenEtaPhiPtNoFound[ij] = fh3GenEtaPhiPt[ij] = 0;
}
for(int ic = 0;ic < kMaxCorrelation;ic++){
Bool_t AliAnalysisTaskJetSpectrum::Notify()
{
+ //
+ // Implemented Notify() to read the cross sections
+ // and number of trials from pyxsec.root
+ //
TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
Double_t xsection = 0;
UInt_t ntrials = 0;
fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
- fh1PtHard_NoW = new TH1F("fh1PtHard_NoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
+ fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
- fh1PtHard_Trials = new TH1F("fh1PtHard_Trials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
+ fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
fh1NGenJets = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
nBinEta,binLimitsEta,nBinEta,binLimitsEta);
- fh2PtGenDeltaPhi[ij] = new TH2F(Form("fh2PtGenDeltaPhi_j%d",ij),"delta phi vs. P_tgen;p_{T,gen} (GeV/c);#phi_{gen}-phi_{rec}",
- 100,-1.0,1.0,nBinPt,binLimitsPt);
- fh2PtGenDeltaPhi[ij] = new TH2F(Form("fh2PtGenDeltaEta_j%d",ij),"delta eta vs. P_tgen;p_{T,gen} (GeV/c);#eta_{gen}-eta_{rec}",
- 100,-1.0,1.0,nBinPt,binLimitsPt);
+ fh2PtGenDeltaPhi[ij] = new TH2F(Form("fh2PtGenDeltaPhi_j%d",ij),"delta phi vs. P_{T,gen};p_{T,gen} (GeV/c);#phi_{gen}-phi_{rec}",
+ nBinPt,binLimitsPt,100,-1.0,1.0);
+ fh2PtGenDeltaEta[ij] = new TH2F(Form("fh2PtGenDeltaEta_j%d",ij),"delta eta vs. p_{T,gen};p_{T,gen} (GeV/c);#eta_{gen}-eta_{rec}",
+ nBinPt,binLimitsPt,100,-1.0,1.0);
- fh3PtRecGenHard_NoW[ij] = new TH3F(Form("fh3PtRecGenHard_NoW_j%d",ij), "Pt hard vs. pt gen vs. pt rec no weight;p_{T,rec};p_{T,gen} (GeV/c);p_{T,hard} (GeV/c)",nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
+ fh3PtRecGenHardNoW[ij] = new TH3F(Form("fh3PtRecGenHardNoW_j%d",ij), "Pt hard vs. pt gen vs. pt rec no weight;p_{T,rec};p_{T,gen} (GeV/c);p_{T,hard} (GeV/c)",nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
fh2Frag[ij] = new TH2F(Form("fh2Frag_j%d",ij),"Jet Fragmentation;x=E_{i}/E_{jet};E_{jet};1/N_{jet}dN_{ch}/dx",
- fh3RecEtaPhiPt_NoGen[ij] = new TH3F(Form("fh3RecEtaPhiPt_NoGen_j%d",ij),"No generated for found jet Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
+ fh3RecEtaPhiPtNoGen[ij] = new TH3F(Form("fh3RecEtaPhiPtNoGen_j%d",ij),"No generated for found jet Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
- fh3GenEtaPhiPt_NoFound[ij] = new TH3F(Form("fh3GenEtaPhiPt_NoFound_j%d",ij),"No found for generated jet eta, phi, pt; #eta; #phi; p_{T,gen} (GeV/c)",
+ fh3GenEtaPhiPtNoFound[ij] = new TH3F(Form("fh3GenEtaPhiPtNoFound_j%d",ij),"No found for generated jet eta, phi, pt; #eta; #phi; p_{T,gen} (GeV/c)",
nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
// Response map
//arrays for bin limits
const Int_t nbin[4] = {100, 100, 100, 100};
- Double_t LowEdge[4] = {0.,0.,0.,0.};
- Double_t UpEdge[4] = {250., 250., 1., 1.};
+ Double_t vLowEdge[4] = {0.,0.,0.,0.};
+ Double_t vUpEdge[4] = {250., 250., 1., 1.};
for(int ic = 0;ic < kMaxCorrelation;ic++){
- fhnCorrelation[ic] = new THnSparseF(Form("fhnCorrelation_%d",ic), "Response Map", 4, nbin, LowEdge, UpEdge);
+ fhnCorrelation[ic] = new THnSparseF(Form("fhnCorrelation_%d",ic), "Response Map", 4, nbin, vLowEdge, vUpEdge);
if(ic==0) fhnCorrelation[ic]->SetTitle(Form("ResponseMap 0 <= npart <= %.0E",fgkJetNpartCut[ic]));
else fhnCorrelation[ic]->SetTitle(Form("ResponseMap %.0E < npart <= %.0E",fgkJetNpartCut[ic-1],fgkJetNpartCut[ic]));
}
fHistList->Add(fh1Xsec);
fHistList->Add(fh1Trials);
fHistList->Add(fh1PtHard);
- fHistList->Add(fh1PtHard_NoW);
- fHistList->Add(fh1PtHard_Trials);
+ fHistList->Add(fh1PtHardNoW);
+ fHistList->Add(fh1PtHardTrials);
fHistList->Add(fh1NGenJets);
fHistList->Add(fh1NRecJets);
////////////////////////
fHistList->Add(fh3RecEtaPhiPt[ij]);
fHistList->Add(fh3GenEtaPhiPt[ij]);
if(saveLevel>2){
- fHistList->Add(fh3RecEtaPhiPt_NoGen[ij]);
- fHistList->Add(fh3GenEtaPhiPt_NoFound[ij]);
+ fHistList->Add(fh3RecEtaPhiPtNoGen[ij]);
+ fHistList->Add(fh3GenEtaPhiPtNoFound[ij]);
}
}
}
if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
fh1PtHard->Fill(ptHard,eventW);
- fh1PtHard_NoW->Fill(ptHard,1);
- fh1PtHard_Trials->Fill(ptHard,nTrials);
+ fh1PtHardNoW->Fill(ptHard,1);
+ fh1PtHardTrials->Fill(ptHard,nTrials);
// If we set a second branch for the input jets fetch this
if(fBranchGen.Length()>0){
Double_t phiRec = recJets[ir].Phi();
if(phiRec<0)phiRec+=TMath::Pi()*2.;
Double_t etaRec = recJets[ir].Eta();
-
+ if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
fh1E[ir]->Fill(recJets[ir].E(),eventW);
fh1PtRecIn[ir]->Fill(ptRec,eventW);
fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
// Fill Correlation
Int_t ig = iGenIndex[ir];
if(ig>=0 && ig<nGenJets){
+ if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+ if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
fh1PtRecOut[ir]->Fill(ptRec,eventW);
Double_t ptGen = genJets[ig].Pt();
Double_t phiGen = genJets[ig].Phi();
if(phiGen<0)phiGen+=TMath::Pi()*2.;
+ if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
Double_t etaGen = genJets[ig].Eta();
+ if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
+ if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
fh2PhiFGen[ir]->Fill(phiRec,phiGen,eventW);
fh2EtaFGen[ir]->Fill(etaRec,etaGen,eventW);
- fh2PtGenDeltaEta[ir]->Fill(etaGen-etaRec,eventW);
- fh2PtGenDeltaEta[ir]->Fill(phiGen-phiRec,eventW);
+ fh2PtGenDeltaEta[ir]->Fill(ptGen,etaGen-etaRec,eventW);
+ if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+ fh2PtGenDeltaPhi[ir]->Fill(ptGen,phiGen-phiRec,eventW);
+ if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
- fh3PtRecGenHard_NoW[ir]->Fill(ptRec,ptGen,ptHard,1);
+ if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+ fh3PtRecGenHardNoW[ir]->Fill(ptRec,ptGen,ptHard,1);
+ if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
/////////////////////////////////////////////////////
- Double_t ERec = recJets[ir].E();
- Double_t EGen = genJets[ig].E();
- fh2Efficiency->Fill(EGen, ERec/EGen);
- if (EGen>=0. && EGen<=250.){
- Double_t Eleading = -1;
+ Double_t eRec = recJets[ir].E();
+ Double_t eGen = genJets[ig].E();
+ fh2Efficiency->Fill(eGen, eRec/eGen);
+ if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+
+ if (eGen>=0. && eGen<=250.){
+ Double_t eLeading = -1;
Double_t ptleading = -1;
Int_t nPart=0;
// loop over tracks
for (Int_t it = 0; it< nTracks; it++){
- if (fAOD->GetTrack(it)->E() > EGen) continue;
+ if (fAOD->GetTrack(it)->E() > eGen) continue;
Double_t phiTrack = fAOD->GetTrack(it)->Phi();
if (phiTrack<0) phiTrack+=TMath::Pi()*2.;
Double_t etaTrack = fAOD->GetTrack(it)->Eta();
Float_t deta = etaRec - etaTrack;
Float_t dphi = TMath::Abs(phiRec - phiTrack);
- Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);
+ Float_t r = TMath::Sqrt(deta*deta + dphi*dphi);
// find leading particle
- if (R<0.4 && fAOD->GetTrack(it)->E()>Eleading){
- Eleading = fAOD->GetTrack(it)->E();
+ if (r<0.4 && fAOD->GetTrack(it)->E()>eLeading){
+ eLeading = fAOD->GetTrack(it)->E();
ptleading = fAOD->GetTrack(it)->Pt();
}
- if (fAOD->GetTrack(it)->Pt()>0.03*EGen && fAOD->GetTrack(it)->E()<=EGen && R<0.7)
+ if (fAOD->GetTrack(it)->Pt()>0.03*eGen && fAOD->GetTrack(it)->E()<=eGen && r<0.7)
nPart++;
}
+ if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+
// fill Response Map (4D histogram) and Energy vs z distributions
- Double_t var[4] = {EGen, ERec, ptleading/EGen, ptleading/ERec};
+ Double_t var[4] = {eGen, eRec, ptleading/eGen, ptleading/eRec};
fh2ERecZRec->Fill(var[1],var[3]); // this has to be filled always in the real case...
fh2EGenZGen->Fill(var[0],var[2]);
fh1JetMultiplicity->Fill(nPart);
- fh3EGenERecN->Fill(EGen, ERec, nPart);
+ fh3EGenERecN->Fill(eGen, eRec, nPart);
for(int ic = 0;ic <kMaxCorrelation;ic++){
- if (nPart<=fgkJetNpartCut[ic]){
- fhnCorrelation[ic]->Fill(var);
- break;
- }
+ if (nPart<=fgkJetNpartCut[ic]){
+ fhnCorrelation[ic]->Fill(var);
+ break;
+ }
}
}
}
////////////////////////////////////////////////////
else{
- fh3RecEtaPhiPt_NoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
+ fh3RecEtaPhiPtNoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
}
}// loop over reconstructed jets
fh1PtGenOut[ig]->Fill(ptGen,eventW);
}
else{
- fh3GenEtaPhiPt_NoFound[ig]->Fill(etaGen,phiGen,ptGen,eventW);
+ fh3GenEtaPhiPtNoFound[ig]->Fill(etaGen,phiGen,ptGen,eventW);
}
}// loop over reconstructed jets
}
-void AliAnalysisTaskJetSpectrum::GetClosestJets(AliAODJet *genJets,Int_t &nGenJets,
- AliAODJet *recJets,Int_t &nRecJets,
+void AliAnalysisTaskJetSpectrum::GetClosestJets(AliAODJet *genJets,const Int_t &nGenJets,
+ AliAODJet *recJets,const Int_t &nRecJets,
Int_t *iGenIndex,Int_t *iRecIndex,Int_t iDebug){
//
/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* See cxx source for full Copyright notice */
-
+
+// **************************************
+// Task used for the correction of determiantion of reconstructed jet spectra
+// Compares input (gen) and output (rec) jets
+// *******************************************
+
#include "AliAnalysisTaskSE.h"
-#include "THnSparse.h"
+#include "THnSparse.h" // cannot forward declare ThnSparseF
+
////////////////
class AliJetHeader;
class AliESDEvent;
class TProfile;
+
class AliAnalysisTaskJetSpectrum : public AliAnalysisTaskSE
{
public:
virtual void SetAODInput(Bool_t b){fUseAODInput = b;}
virtual void SetLimitGenJetEta(Bool_t b){fLimitGenJetEta = b;}
virtual void SetAnalysisType(Int_t i){fAnalysisType = i;}
- virtual void SetBranchGen(char* c){fBranchGen = c;}
- virtual void SetBranchRec(char* c){fBranchRec = c;}
+ virtual void SetBranchGen(const char* c){fBranchGen = c;}
+ virtual void SetBranchRec(const char* c){fBranchRec = c;}
// Helper
- static void GetClosestJets(AliAODJet *genJets,Int_t &nGenJets,
- AliAODJet *recJets,Int_t &nRecJets,
+ static void GetClosestJets(AliAODJet *genJets,const Int_t &nGenJets,
+ AliAODJet *recJets,const Int_t &nRecJets,
Int_t *iGenIndex,Int_t *iRecIndex,Int_t iDebug = 0);
//
TProfile* fh1Xsec; // pythia cross section and trials
TH1F* fh1Trials; // trials are added
TH1F* fh1PtHard; // Pt har of the event...
- TH1F* fh1PtHard_NoW; // Pt har of the event...
- TH1F* fh1PtHard_Trials; // Number of trials
+ TH1F* fh1PtHardNoW; // Pt har of the event...
+ TH1F* fh1PtHardTrials; // Number of trials
TH1F* fh1NGenJets;
TH1F* fh1NRecJets;
TH1F* fh1E[kMaxJets]; // Jet Energy
TH1F* fh1PtRecIn[kMaxJets]; // Jet pt for all
TH1F* fh1PtRecOut[kMaxJets]; // Jet pt with corellated generated jet
TH1F* fh1PtGenIn[kMaxJets]; // Detection efficiency for given p_T.gen
- TH1F* fh1PtGenOut[kMaxJets]; //
+ TH1F* fh1PtGenOut[kMaxJets]; // gen pT of found jets
- TH2F* fh2PtFGen[kMaxJets]; //
- TH2F* fh2PhiFGen[kMaxJets]; //
- TH2F* fh2EtaFGen[kMaxJets]; //
+ TH2F* fh2PtFGen[kMaxJets]; // correlation betwen genreated and found jet pT
+ TH2F* fh2PhiFGen[kMaxJets]; // correlation betwen genreated and found jet phi
+ TH2F* fh2EtaFGen[kMaxJets]; // correlation betwen genreated and found jet eta
TH2F* fh2Frag[kMaxJets]; // fragmentation function
- TH2F* fh2FragLn[kMaxJets]; //
- TH2F* fh2PtGenDeltaPhi[kMaxJets];
- TH2F* fh2PtGenDeltaEta[kMaxJets];
-
- TH3F* fh3PtRecGenHard[kMaxJets]; //
- TH3F* fh3PtRecGenHard_NoW[kMaxJets]; //
- TH3F* fh3RecEtaPhiPt[kMaxJets]; //
- TH3F* fh3RecEtaPhiPt_NoGen[kMaxJets]; //
- TH3F* fh3GenEtaPhiPt_NoFound[kMaxJets]; //
- TH3F* fh3GenEtaPhiPt[kMaxJets]; //
+ TH2F* fh2FragLn[kMaxJets]; // fragmetation in xi
+ TH2F* fh2PtGenDeltaPhi[kMaxJets]; // difference between generated and found jet phi
+ TH2F* fh2PtGenDeltaEta[kMaxJets]; // difference between generated and found jet eta
+
+ TH3F* fh3PtRecGenHard[kMaxJets]; // correlation beetwen pt hard, rec and gen
+ TH3F* fh3PtRecGenHardNoW[kMaxJets]; // correlation beetwen pt hard, rec and gen no weight
+ TH3F* fh3RecEtaPhiPt[kMaxJets]; // correlation between eta phi and rec pt
+ TH3F* fh3RecEtaPhiPtNoGen[kMaxJets]; // correlation between eta phi and rec pt of jets without a partner
+ TH3F* fh3GenEtaPhiPtNoFound[kMaxJets]; // correlation between eta phi and gen pt of jets without a partner
+ TH3F* fh3GenEtaPhiPt[kMaxJets]; // correlation between eta phi and gen pt of jets without a partner
// ========= Multiplicity dependence ======
// ==========TODO , flavaor dependence ========
TList *fHistList; // Output list
///////// For 2 dimensional unfolding //////////////////
- TH1F* fh1JetMultiplicity;
- TH2F* fh2ERecZRec;
- TH2F* fh2EGenZGen;
- TH2F* fh2Efficiency;
- TH3F* fh3EGenERecN;
- THnSparseF* fhnCorrelation[kMaxCorrelation];
+ TH1F* fh1JetMultiplicity; // JetMultiplicity
+ TH2F* fh2ERecZRec; // rec e and z
+ TH2F* fh2EGenZGen; // gen e and z
+ TH2F* fh2Efficiency; // 2 dimensional efficiency
+ TH3F* fh3EGenERecN; // correlation rec energy, gen energy N particles
+ THnSparseF* fhnCorrelation[kMaxCorrelation]; // 4d Correllation for unfolding
///////////////////////////////////////////////////////
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff