+// **************************************
+// 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 <TRandom.h>
#include <TSystem.h>
#include <TInterpreter.h>
#include <TChain.h>
#include <TFile.h>
+#include <TKey.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include "AliJetReader.h"
#include "AliJetReaderHeader.h"
#include "AliUA1JetHeaderV1.h"
-#include "AliJet.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliAODHandler.h"
#include "AliGenCocktailEventHeader.h"
#include "AliInputEventHandler.h"
+
#include "AliAnalysisHelperJetTasks.h"
ClassImp(AliAnalysisTaskJetSpectrum)
fJetHeaderGen(0x0),
fAOD(0x0),
fBranchRec("jets"),
- fConfigRec("ConfigJets.C"),
fBranchGen(""),
- fConfigGen(""),
fUseAODInput(kFALSE),
fUseExternalWeightOnly(kFALSE),
fLimitGenJetEta(kFALSE),
fAnalysisType(0),
fExternalWeight(1),
+ fRecEtaWindow(0.5),
fh1Xsec(0x0),
fh1Trials(0x0),
fh1PtHard(0x0),
- fh1PtHard_NoW(0x0),
- fh1PtHard_Trials(0x0),
+ fh1PtHardNoW(0x0),
+ fh1PtHardTrials(0x0),
fh1NGenJets(0x0),
fh1NRecJets(0x0),
- fHistList(0x0) ,
+ fHistList(0x0) ,
////////////////
fh1JetMultiplicity(0x0) ,
fh2ERecZRec(0x0) ,
// Default constructor
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;
+ fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] = fh2Frag[ij] = fh2FragLn[ij] = fh2PtRecDeltaR[ij] = fh2PtGenDeltaR[ij] = fh2PtGenDeltaPhi[ij] = fh2PtGenDeltaEta[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;
fJetHeaderGen(0x0),
fAOD(0x0),
fBranchRec("jets"),
- fConfigRec("ConfigJets.C"),
fBranchGen(""),
- fConfigGen(""),
fUseAODInput(kFALSE),
fUseExternalWeightOnly(kFALSE),
fLimitGenJetEta(kFALSE),
fAnalysisType(0),
fExternalWeight(1),
+ fRecEtaWindow(0.5),
fh1Xsec(0x0),
fh1Trials(0x0),
fh1PtHard(0x0),
- fh1PtHard_NoW(0x0),
- fh1PtHard_Trials(0x0),
+ fh1PtHardNoW(0x0),
+ fh1PtHardTrials(0x0),
fh1NGenJets(0x0),
fh1NRecJets(0x0),
fHistList(0x0) ,
// Default constructor
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;
+ fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] = fh2Frag[ij] = fh2FragLn[ij] = fh2PtGenDeltaPhi[ij] = fh2PtGenDeltaEta[ij] = fh2PtRecDeltaR[ij] = fh2PtGenDeltaR[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;
+ Float_t ftrials = 0;
if(tree){
TFile *curfile = tree->GetCurrentFile();
if (!curfile) {
TString fileName(curfile->GetName());
if(fileName.Contains("AliESDs.root")){
- fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
+ fileName.ReplaceAll("AliESDs.root", "");
}
else if(fileName.Contains("AliAOD.root")){
- fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
+ fileName.ReplaceAll("AliAOD.root", "");
+ }
+ else if(fileName.Contains("AliAODs.root")){
+ fileName.ReplaceAll("AliAODs.root", "");
}
else if(fileName.Contains("galice.root")){
// for running with galice and kinematics alone...
- fileName.ReplaceAll("galice.root", "pyxsec.root");
+ fileName.ReplaceAll("galice.root", "");
}
- TFile *fxsec = TFile::Open(fileName.Data());
+ TFile *fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec.root"));
if(!fxsec){
- Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data());
- // no a severe condition
- return kTRUE;
+ if(fDebug>0)Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,Form("%s%s",fileName.Data(),"pyxsec.root"));
+ // next trial fetch the histgram file
+ fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec_hists.root"));
+ if(!fxsec){
+ // not a severe condition
+ if(fDebug>0)Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,Form("%s%s",fileName.Data(),"pyxsec_hists.root"));
+ return kTRUE;
+ }
+ else{
+ // find the tlist we want to be independtent of the name so use the Tkey
+ TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
+ if(!key){
+ if(fDebug>0)Printf("%s:%d key not found in the file",(char*)__FILE__,__LINE__);
+ return kTRUE;
+ }
+ TList *list = dynamic_cast<TList*>(key->ReadObj());
+ if(!list){
+ if(fDebug>0)Printf("%s:%d key is not a tlist",(char*)__FILE__,__LINE__);
+ return kTRUE;
+ }
+ xsection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
+ ftrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
+ }
}
- TTree *xtree = (TTree*)fxsec->Get("Xsection");
- if(!xtree){
- Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
+ else{
+ TTree *xtree = (TTree*)fxsec->Get("Xsection");
+ if(!xtree){
+ Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
+ }
+ xtree->SetBranchAddress("xsection",&xsection);
+ xtree->SetBranchAddress("ntrials",&ntrials);
+ ftrials = ntrials;
+ xtree->GetEntry(0);
}
- xtree->SetBranchAddress("xsection",&xsection);
- xtree->SetBranchAddress("ntrials",&ntrials);
- xtree->GetEntry(0);
fh1Xsec->Fill("<#sigma>",xsection);
- fh1Trials->Fill("#sum{ntrials}",ntrials);
+ fh1Trials->Fill("#sum{ntrials}",ftrials);
}
return kTRUE;
binLimitsPt[iPt] = 0.0;
}
else {// 1.0
- binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2;
- }
- }
- const Int_t nBinEta = 22;
- Double_t binLimitsEta[nBinEta+1];
- for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
- if(iEta==0){
- binLimitsEta[iEta] = -2.2;
- }
- else{
- binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.2;
+ binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.5;
}
}
+
+ const Int_t nBinEta = 26;
+ Double_t binLimitsEta[nBinEta+1] = {
+ -1.6,-1.4,-1.2,-1.0,
+ -0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,
+ 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
+ 1.0, 1.2, 1.4, 1.6
+ };
+
const Int_t nBinPhi = 30;
Double_t binLimitsPhi[nBinPhi+1];
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);
fh1PtGenOut[ij] = new TH1F(Form("fh1PtGenOut_j%d",ij),"found p_T output jets;p_{T,gen}",nBinPt,binLimitsPt);
+
fh2PtFGen[ij] = new TH2F(Form("fh2PtFGen_j%d",ij),"Pt Found vs. gen;p_{T,rec} (GeV/c);p_{T,gen} (GeV/c)",
nBinPt,binLimitsPt,nBinPt,binLimitsPt);
- fh2PhiFGen[ij] = new TH2F(Form("fh2PhiFGen_j%d",ij),"#phi Found vs. gen;#phi_{rec};phi_{gen}",
+ fh2PhiFGen[ij] = new TH2F(Form("fh2PhiFGen_j%d",ij),"#phi Found vs. gen;#phi_{rec};#phi_{gen}",
nBinPhi,binLimitsPhi,nBinPhi,binLimitsPhi);
- fh2EtaFGen[ij] = new TH2F(Form("fh2EtaFGen_j%d",ij),"#eta Found vs. gen;#eta_{rec};eta_{gen}",
+ fh2EtaFGen[ij] = new TH2F(Form("fh2EtaFGen_j%d",ij),"#eta Found vs. gen;#eta_{rec};#eta_{gen}",
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);
+ fh2PtRecDeltaR[ij] = new TH2F(Form("fh2PtRecDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,rec,j};#Delta R",
+ nBinPt,binLimitsPt,60,0,6.0);
+ fh2PtGenDeltaR[ij] = new TH2F(Form("fh2PtGenDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,gen,j};#Delta R",
+ nBinPt,binLimitsPt,60,0,6.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);
- fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
+ fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,} (GeV/c)",
nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
}
+
+ // tmp histos do not add to the header
+ TH2F *hCorrPt = new TH2F("fh2PtRecPhiCorrPt","#Delta#phi correlation pt weighted",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
+ fHistList->Add(hCorrPt);
+ TH2F *hCorrRanPt = new TH2F("fh2PtRecPhiCorrPtRan","#Delta#phi Random correlation pt weighted",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
+ fHistList->Add(hCorrRanPt);
+
+ TH2F *hCorr = new TH2F("fh2PtRecPhiCorr","#Delta#phi correlation",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
+ fHistList->Add(hCorr);
+ TH2F *hCorrRan = new TH2F("fh2PtRecPhiCorrRan","#Delta#phi Random correlation",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
+ fHistList->Add(hCorrRan);
+
+
/////////////////////////////////////////////////////////////////
fh1JetMultiplicity = new TH1F("fh1JetMultiplicity", "Jet Multiplicity", 51, 0., 50.);
// 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(fh2EtaFGen[ij]);
fHistList->Add(fh2PtGenDeltaEta[ij]);
fHistList->Add(fh2PtGenDeltaPhi[ij]);
+ fHistList->Add(fh2PtRecDeltaR[ij]);
+ fHistList->Add(fh2PtGenDeltaR[ij]);
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]);
}
}
}
return;
}
-
- // aodH->SelectEvent(kTRUE);
-
- // ========= These pointers need to be valid in any case =======
-
-
- /*
- AliUA1JetHeaderV1 *jhRec = dynamic_cast<AliUA1JetHeaderV1*>(fJetFinderRec->GetHeader());
- if(!jhRec){
- Printf("%s:%d No Jet Header found",(char*)__FILE__,__LINE__);
- return;
- }
- */
if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+
TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
if(!aodRecJets){
Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
nTrials = pythiaGenHeader->Trials();
ptHard = pythiaGenHeader->GetPtHard();
-
+ int iProcessType = pythiaGenHeader->ProcessType();
+ // 11 f+f -> f+f
+ // 12 f+barf -> f+barf
+ // 13 f+barf -> g+g
+ // 28 f+g -> f+g
+ // 53 g+g -> f+barf
+ // 68 g+g -> g+g
+ /*
+ if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
+ // if(iProcessType != 13 && iProcessType != 68){ // allow only glue
+ if(iProcessType != 11 && iProcessType != 12 && iProcessType != 53){ // allow only quark
+ // if(iProcessType != 28){ // allow only -> f+g
+ PostData(1, fHistList);
+ return;
+ }
+ */
+ if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
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){
}
- GetClosestJets(genJets,nGenJets,recJets,nRecJets,
+ AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
iGenIndex,iRecIndex,fDebug);
if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
for(int ir = 0;ir < nRecJets;++ir){
Double_t ptRec = recJets[ir].Pt();
Double_t phiRec = recJets[ir].Phi();
+ Double_t phiRecRan = TMath::Pi()*gRandom->Rndm(); // better take real jet axis from previous events (TPC acceptance in 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);
+ for(int irr = ir+1;irr<nRecJets;irr++){
+ fh2PtRecDeltaR[ir]->Fill(recJets[irr].Pt(),recJets[ir].DeltaR(&recJets[irr]));
+ }
// 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.;
Double_t etaGen = genJets[ig].Eta();
+
+ //
+ // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
+ //
+
+ if(TMath::Abs(etaRec)<fRecEtaWindow){
+
fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
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);
+ fh2PtGenDeltaPhi[ir]->Fill(ptGen,phiGen-phiRec,eventW);
fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
- fh3PtRecGenHard_NoW[ir]->Fill(ptRec,ptGen,ptHard,1);
+ fh3PtRecGenHardNoW[ir]->Fill(ptRec,ptGen,ptHard,1);
/////////////////////////////////////////////////////
- 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();
+ // CKB use p_T not Energy
+ // TODO recname variabeles and histos
+ 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 ptleading = -1;
Int_t nPart=0;
// loop over tracks
for (Int_t it = 0; it< nTracks; it++){
- 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);
+ // if (fAOD->GetTrack(it)->E() > eGen) continue; // CKB. Not allowed! cannot cut on gen properties in real events!
// 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){
+ // TODO implement esd filter flag to be the same as in the jet finder
+ // allow also for MC particles...
+ Float_t r = recJets[ir].DeltaR(fAOD->GetTrack(it));
+ if (r<0.4 && fAOD->GetTrack(it)->Pt()>ptleading){
+ 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) // CKB cannot cut on gen properties
+ if (fAOD->GetTrack(it)->Pt()>0.03*eRec && fAOD->GetTrack(it)->Pt()<=eRec && r<0.7)
nPart++;
+
+
+ // correlate jet axis of leading jet with particles
+ if(ir==0){
+ Float_t phi = fAOD->GetTrack(it)->Phi();
+ Float_t dPhi = phi - phiRec;
+ if(dPhi>TMath::Pi()/1.5)dPhi = dPhi - 2.*TMath::Pi();
+ if(dPhi<(-0.5*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
+ Float_t dPhiRan = phi - phiRecRan;
+ if(dPhiRan>TMath::Pi()/1.5)dPhiRan = dPhiRan - 2.*TMath::Pi();
+ if(dPhiRan<(-0.5*TMath::Pi()))dPhiRan = dPhiRan + 2.*TMath::Pi();
+ ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorr"))->Fill(ptRec,dPhi);
+ ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrRan"))->Fill(ptRec,dPhiRan);
+ ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrPt"))->Fill(ptRec,dPhi,fAOD->GetTrack(it)->Pt());
+ ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrPtRan"))->Fill(ptRec,dPhiRan,fAOD->GetTrack(it)->Pt());
+
+ }
}
+ 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]){ // is this corrected for CKB
+ fhnCorrelation[ic]->Fill(var);
+ break;
+ }
}
}
- }
+
+ }// if etarec in window
+
+ }
////////////////////////////////////////////////////
else{
- fh3RecEtaPhiPt_NoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
+ fh3RecEtaPhiPtNoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
}
}// loop over reconstructed jets
Double_t etaGen = genJets[ig].Eta();
fh3GenEtaPhiPt[ig]->Fill(etaGen,phiGen,ptGen,eventW);
fh1PtGenIn[ig]->Fill(ptGen,eventW);
+ for(int igg = ig+1;igg<nGenJets;igg++){
+ fh2PtGenDeltaR[ig]->Fill(genJets[igg].Pt(),genJets[ig].DeltaR(&genJets[igg]));
+ }
Int_t ir = iRecIndex[ig];
if(ir>=0&&ir<nRecJets){
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
//
if (fDebug > 1) printf("AnalysisJetSpectrum: Terminate() \n");
}
-
-
-void AliAnalysisTaskJetSpectrum::GetClosestJets(AliAODJet *genJets,Int_t &nGenJets,
- AliAODJet *recJets,Int_t &nRecJets,
- Int_t *iGenIndex,Int_t *iRecIndex,Int_t iDebug){
-
- //
- // Relate the two input jet Arrays
- //
-
- //
- // The association has to be unique
- // So check in two directions
- // find the closest rec to a gen
- // and check if there is no other rec which is closer
- // Caveat: Close low energy/split jets may disturb this correlation
-
- // Idea: search in two directions generated e.g (a--e) and rec (1--3)
- // Fill a matrix with Flags (1 for closest rec jet, 2 for closest rec jet
- // in the end we have something like this
- // 1 2 3
- // ------------
- // a| 3 2 0
- // b| 0 1 0
- // c| 0 0 3
- // d| 0 0 1
- // e| 0 0 1
- // Topology
- // 1 2
- // a b
- //
- // d c
- // 3 e
- // Only entries with "3" match from both sides
-
- const int kMode = 3;
-
- Int_t iFlag[kMaxJets][kMaxJets];
-
-
-
- for(int i = 0;i < kMaxJets;++i){
- iRecIndex[i] = -1;
- iGenIndex[i] = -1;
- for(int j = 0;j < kMaxJets;++j)iFlag[i][j] = 0;
- }
-
- if(nRecJets==0)return;
- if(nGenJets==0)return;
-
- const Float_t maxDist = 1.0;
- // find the closest distance to the generated
- for(int ig = 0;ig<nGenJets;++ig){
- Float_t dist = maxDist;
- if(iDebug>1)Printf("Gen (%d) p_T %3.3f eta %3.3f ph %3.3f ",ig,genJets[ig].Pt(),genJets[ig].Eta(),genJets[ig].Phi());
- for(int ir = 0;ir<nRecJets;++ir){
- Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
- if(iDebug>1)Printf("Rec (%d) p_T %3.3f eta %3.3f ph %3.3f ",ir,recJets[ir].Pt(),recJets[ir].Eta(),recJets[ir].Phi());
- if(iDebug>1)Printf("Distance (%d)--(%d) %3.3f ",ig,ir,dR);
- if(dR<dist){
- iRecIndex[ig] = ir;
- dist = dR;
- }
- }
- if(iRecIndex[ig]>=0)iFlag[ig][iRecIndex[ig]]+=1;
- // reset...
- iRecIndex[ig] = -1;
- }
- // other way around
- for(int ir = 0;ir<nRecJets;++ir){
- Float_t dist = maxDist;
- for(int ig = 0;ig<nGenJets;++ig){
- Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
- if(dR<dist){
- iGenIndex[ir] = ig;
- dist = dR;
- }
- }
- if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]][ir]+=2;
- // reset...
- iGenIndex[ir] = -1;
- }
-
- // check for "true" correlations
-
- if(iDebug>1)Printf(">>>>>> Matrix");
-
- for(int ig = 0;ig<nGenJets;++ig){
- for(int ir = 0;ir<nRecJets;++ir){
- // Print
- if(iDebug>1)printf("XFL %d ",iFlag[ig][ir]);
-
- if(kMode==3){
- // we have a uniqie correlation
- if(iFlag[ig][ir]==3){
- iGenIndex[ir] = ig;
- iRecIndex[ig] = ir;
- }
- }
- else{
- // we just take the correlation from on side
- if((iFlag[ig][ir]&2)==2){
- iGenIndex[ir] = ig;
- }
- if((iFlag[ig][ir]&1)==1){
- iRecIndex[ig] = ir;
- }
- }
- }
- if(iDebug>1)printf("\n");
- }
-}
-
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff