X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=PWGJE%2FAliAnalysisTaskJetSpectrum2.cxx;h=158a8f347f2bbfedace7cb46274da43f8c4f80d9;hb=f15c1f69167973bb1eb7e04eb04c846e5529df1a;hp=c99f2a69bd535ae88a60866eb9f8f707650d9556;hpb=441afcb96e84fce4bd3498d32a7167a5c999f65e;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWGJE/AliAnalysisTaskJetSpectrum2.cxx b/PWGJE/AliAnalysisTaskJetSpectrum2.cxx index c99f2a69bd5..158a8f347f2 100644 --- a/PWGJE/AliAnalysisTaskJetSpectrum2.cxx +++ b/PWGJE/AliAnalysisTaskJetSpectrum2.cxx @@ -18,7 +18,6 @@ * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ - #include #include @@ -61,14 +60,14 @@ #include "AliJetKineReaderHeader.h" #include "AliGenCocktailEventHeader.h" #include "AliInputEventHandler.h" - +#include "AliCFContainer.h" #include "AliAnalysisHelperJetTasks.h" ClassImp(AliAnalysisTaskJetSpectrum2) AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): - AliAnalysisTaskSE(), +AliAnalysisTaskSE(), fJetHeaderRec(0x0), fJetHeaderGen(0x0), fAODIn(0x0), @@ -93,11 +92,17 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): fDoMatching(kFALSE), fNMatchJets(5), fNRPBins(3), + fTRP(0), + fDebug(0), fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered), + fJetTriggerBestMask(AliAODJet::kHighTrackPtBest), fFilterMask(0), fEventSelectionMask(0), fNTrigger(0), fTriggerBit(0x0), + fNAcceptance(0), + fNBinsLeadingTrackPt(10), + fNBinsMult(20), fAnalysisType(0), fTrackTypeRec(kTrackUndef), fTrackTypeGen(kTrackUndef), @@ -110,6 +115,10 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): fMinJetPt(0), fMinTrackPt(0.15), fDeltaPhiWindow(90./180.*TMath::Pi()), + fAcceptancePhiMin(0x0), + fAcceptancePhiMax(0x0), + fAcceptanceEtaMin(0x0), + fAcceptanceEtaMax(0x0), fCentrality(100), fRPAngle(0), fMultRec(0), @@ -117,6 +126,7 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): fTriggerName(0x0), fh1Xsec(0x0), fh1Trials(0x0), + fh1AvgTrials(0x0), fh1PtHard(0x0), fh1PtHardNoW(0x0), fh1PtHardTrials(0x0), @@ -128,7 +138,13 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): fh2MultGen(0x0), fh2RPCentrality(0x0), fh2PtFGen(0x0), + fh2deltaPt1Pt2(0x0), fh2RelPtFGen(0x0), + fh3RelPtFGenLeadTrkPt(0x0), + fh1EvtSelection(0), + fMaxVertexZ(100.), + fMinNcontributors(0), + fRejectPileup(0), fHistList(0x0) { @@ -138,6 +154,7 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): fh1SumPtTrack[ij] = 0; fh1PtJetsIn[ij] = 0; fh1PtJetsInRej[ij] = 0; + fh1PtJetsInBest[ij] = 0; fh1PtTracksIn[ij] = 0; fh1PtTracksInLow[ij] = 0; fh2NJetsPt[ij] = 0; @@ -145,6 +162,8 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): fp2MultRPPhiTrackPt[ij] = 0; fp2CentRPPhiTrackPt[ij] = 0; fhnJetPt[ij] = 0; + fhnJetPtBest[ij] = 0; + fhnJetPtRej[ij] = 0; fhnJetPtQA[ij] = 0; fhnTrackPt[ij] = 0; fhnTrackPtQA[ij] = 0; @@ -188,11 +207,17 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name): fDoMatching(kFALSE), fNMatchJets(5), fNRPBins(3), + fTRP(0), + fDebug(0), fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered), + fJetTriggerBestMask(AliAODJet::kHighTrackPtBest), fFilterMask(0), fEventSelectionMask(0), fNTrigger(0), fTriggerBit(0x0), + fNAcceptance(0), + fNBinsLeadingTrackPt(10), + fNBinsMult(20), fAnalysisType(0), fTrackTypeRec(kTrackUndef), fTrackTypeGen(kTrackUndef), @@ -205,6 +230,10 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name): fMinJetPt(0), fMinTrackPt(0.15), fDeltaPhiWindow(90./180.*TMath::Pi()), + fAcceptancePhiMin(0x0), + fAcceptancePhiMax(0x0), + fAcceptanceEtaMin(0x0), + fAcceptanceEtaMax(0x0), fCentrality(100), fRPAngle(0), fMultRec(0), @@ -212,6 +241,7 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name): fTriggerName(0x0), fh1Xsec(0x0), fh1Trials(0x0), + fh1AvgTrials(0x0), fh1PtHard(0x0), fh1PtHardNoW(0x0), fh1PtHardTrials(0x0), @@ -223,7 +253,13 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name): fh2MultGen(0x0), fh2RPCentrality(0x0), fh2PtFGen(0x0), + fh2deltaPt1Pt2(0x0), fh2RelPtFGen(0x0), + fh3RelPtFGenLeadTrkPt(0x0), + fh1EvtSelection(0), + fMaxVertexZ(100.), + fMinNcontributors(0), + fRejectPileup(0), fHistList(0x0) { @@ -233,6 +269,7 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name): fh1SumPtTrack[ij] = 0; fh1PtJetsIn[ij] = 0; fh1PtJetsInRej[ij] = 0; + fh1PtJetsInBest[ij] = 0; fh1PtTracksIn[ij] = 0; fh1PtTracksInLow[ij] = 0; fh2NJetsPt[ij] = 0; @@ -240,6 +277,8 @@ AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name): fp2MultRPPhiTrackPt[ij] = 0; fp2CentRPPhiTrackPt[ij] = 0; fhnJetPt[ij] = 0; + fhnJetPtBest[ij] = 0; + fhnJetPtRej[ij] = 0; fhnJetPtQA[ij] = 0; fhnTrackPt[ij] = 0; fhnTrackPtQA[ij] = 0; @@ -308,6 +347,7 @@ Bool_t AliAnalysisTaskJetSpectrum2::Notify() // construct a poor man average trials Float_t nEntries = (Float_t)tree->GetTree()->GetEntries(); if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries; + fh1Trials->Fill("#sum{ntrials}",ftrials); } if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName()); @@ -335,13 +375,13 @@ void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects() // event npsparse cent, mult - const Int_t nBinsSparse0 = 3; - const Int_t nBins0[nBinsSparse0] = { 100, 500,fNTrigger}; - const Double_t xmin0[nBinsSparse0] = { 0, 0, -0.5}; - const Double_t xmax0[nBinsSparse0] = { 100,5000,fNTrigger-0.5}; + const Int_t nBinsSparse0 = 4; + const Int_t nBins0[nBinsSparse0] = { 100, 500,fNTrigger,125}; + const Double_t xmin0[nBinsSparse0] = { 0, 0, -0.5,-2}; + const Double_t xmax0[nBinsSparse0] = { 100,5000,fNTrigger-0.5,248}; - fhnEvent = new THnSparseF("fhnEvent",";cent;mult",nBinsSparse0, + fhnEvent = new THnSparseF("fhnEvent",";cent;mult:trigger;#rho",nBinsSparse0, nBins0,xmin0,xmax0); fHistList->Add(fhnEvent); @@ -354,15 +394,25 @@ void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects() // Histogram + fh1EvtSelection = new TH1F("fh1EvtSelection", "Event Selection", 7, -0.5, 6.5); + fh1EvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED"); + fh1EvtSelection->GetXaxis()->SetBinLabel(2,"event selection: rejected"); + fh1EvtSelection->GetXaxis()->SetBinLabel(3,"event class: rejected"); + fh1EvtSelection->GetXaxis()->SetBinLabel(4,"vertex Ncontr: rejected"); + fh1EvtSelection->GetXaxis()->SetBinLabel(5,"vertex z: rejected"); + fh1EvtSelection->GetXaxis()->SetBinLabel(6,"vertex type: rejected"); + fh1EvtSelection->GetXaxis()->SetBinLabel(7,"pileup: rejected"); - const Int_t nBinPt = 120; + fHistList->Add(fh1EvtSelection); + + const Int_t nBinPt = 150; Double_t binLimitsPt[nBinPt+1]; for(Int_t iPt = 0;iPt <= nBinPt;iPt++){ if(iPt == 0){ binLimitsPt[iPt] = -50.0; } else {// 1.0 - binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.5; + binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.; } } const Int_t nBinPhi = 90; @@ -394,6 +444,9 @@ void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects() fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1); fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}"); fHistList->Add(fh1Trials); + fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1); + fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}"); + fHistList->Add(fh1AvgTrials); fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt); fHistList->Add(fh1PtHard); fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt); @@ -423,9 +476,30 @@ void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects() fh2PtFGen = new TH2F("fh2PtFGen",Form("%s vs. %s;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt); fHistList->Add(fh2PtFGen); + // fh2deltaPt1Pt2(0x0), + fh2deltaPt1Pt2 = new TH3F("fh2deltaPt1Pt2",Form("%s vs. %s;delta;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt); + fHistList->Add(fh2deltaPt1Pt2); + fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen}",nBinPt,binLimitsPt,241,-2.41,2.41); fHistList->Add(fh2RelPtFGen); - + + const Int_t nBinsLeadingTrackPt = 10; + const Double_t binArrayLeadingTrackPt[nBinsLeadingTrackPt+1] = {0.,1.,2.,3.,4.,5.,6.,8.,10.,12.,200.}; //store pT of leading track in jet + + const Int_t nBinDeltaPt = 241; + Double_t binLimitsDeltaPt[nBinDeltaPt+1]; + for(Int_t iDeltaPt = 0;iDeltaPt <= nBinDeltaPt;iDeltaPt++){ + if(iDeltaPt == 0){ + binLimitsDeltaPt[iDeltaPt] = -2.41; + } + else { + binLimitsDeltaPt[iDeltaPt] = binLimitsDeltaPt[iDeltaPt-1] + 0.02; + } + } + + fh3RelPtFGenLeadTrkPt = new TH3F("fh3RelPtFGenLeadTrkPt",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen};p_{T}^{leading track}",nBinPt,binLimitsPt,nBinDeltaPt,binLimitsDeltaPt,nBinsLeadingTrackPt,binArrayLeadingTrackPt); + fHistList->Add(fh3RelPtFGenLeadTrkPt); + for(int ij = 0;ij Add(fh1PtJetsInRej[ij]); + + fh1PtJetsInBest[ij] = new TH1F(Form("fh1PtJets%sInBest",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt); + fHistList->Add(fh1PtJetsInBest[ij]); fh1PtTracksIn[ij] = new TH1F(Form("fh1PtTracks%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt); fHistList->Add(fh1PtTracksIn[ij]); @@ -479,102 +556,131 @@ void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects() fp2CentRPPhiTrackPt[ij] = new TProfile2D(Form("fp2CentRPPhiTrackPt%s",cAdd.Data()),"RP phi vs cent;# cent;#Delta#phi_{RP}; ",10,0,100,181,-1./180.*TMath::Pi(),TMath::Pi(),"S"); fHistList->Add(fp2CentRPPhiTrackPt[ij]); - // Bins: Jet number: pTJet, cent, mult, RP, Area, trigger total bins = 4.5M - const Int_t nBinsSparse1 = 7; - Int_t nBins1[nBinsSparse1] = { kMaxJets+1,120, 10, 25, fNRPBins, 10,fNTrigger}; + // Bins: Jet number: pTJet, cent, mult, RP, Area, trigger, leading track pT total bins = 4.5M + const Int_t nBinsSparse1 = 9; + const Int_t nBinsArea = 10; + Int_t nbinr5=120; + Int_t nbinr5min=-50; + + if(cJetBranch.Contains("KT05")){ + nbinr5=132; + nbinr5min=-80;} + + Int_t nBins1[nBinsSparse1] = { kMaxJets+1,nbinr5, 10, fNBinsMult, fNRPBins, nBinsArea,fNTrigger,fNBinsLeadingTrackPt,fNAcceptance+1}; if(cJetBranch.Contains("RandomCone")){ nBins1[1] = 600; nBins1[5] = 1; } - const Double_t xmin1[nBinsSparse1] = { -0.5,-50, 0, 0, -0.5, 0.,-0.5}; - const Double_t xmax1[nBinsSparse1] = {kMaxJets+0.5,250,100,5000,fNRPBins-0.5,1.0,fNTrigger-0.5}; - - fhnJetPt[ij] = new THnSparseF(Form("fhnJetPt%s",cAdd.Data()),";jet number;p_{T,jet};cent;# tracks;RP;area;trigger",nBinsSparse1,nBins1,xmin1,xmax1); + const Double_t xmin1[nBinsSparse1] = { -0.5, static_cast(nbinr5min), 0, 0, -0.5, 0., -0.5, 0., -0.5,}; + const Double_t xmax1[nBinsSparse1] = {kMaxJets+0.5,250,100,4000,fNRPBins-0.5,1.0,fNTrigger-0.5,200.,fNAcceptance+0.5}; + + const Double_t binArrayArea[nBinsArea+1] = {xmin1[5],0.07,0.1,0.15,0.2,0.25,0.3,0.4,0.5,0.6,xmax1[5]}; + + + fhnJetPt[ij] = new THnSparseF(Form("fhnJetPt%s",cAdd.Data()),";jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leading track p_{T};acceptance bin",nBinsSparse1,nBins1,xmin1,xmax1); + fhnJetPt[ij]->SetBinEdges(5,binArrayArea); + if(fNBinsLeadingTrackPt>1) fhnJetPt[ij]->SetBinEdges(7,binArrayLeadingTrackPt); fHistList->Add(fhnJetPt[ij]); + + + // Bins: pTJet, cent, trigger, + const Int_t nBinsSparse1b = 3; + Int_t nBins1b[nBinsSparse1b] = {120, 10,fNTrigger}; + const Double_t xmin1b[nBinsSparse1b] = {-50, 0,-0.5}; + const Double_t xmax1b[nBinsSparse1b] = {250,100,fNTrigger-0.5}; + + fhnJetPtBest[ij] = new THnSparseF(Form("fhnJetPtBest%s",cAdd.Data()),";p_{T,jet};cent;trigger",nBinsSparse1b,nBins1b,xmin1b,xmax1b); + fHistList->Add(fhnJetPtBest[ij]); + + fhnJetPtRej[ij] = new THnSparseF(Form("fhnJetPtRej%s",cAdd.Data()),";p_{T,jet};cent;trigger",nBinsSparse1b,nBins1b,xmin1b,xmax1b); + fHistList->Add(fhnJetPtRej[ij]); - // Bins: Jet number: pTJet, cent, eta, phi, Area. total bins = 9.72 M - const Int_t nBinsSparse2 = 6; - Int_t nBins2[nBinsSparse2] = { kMaxJets+1, 25, 5, 18, 360, 10}; + // Bins: Jet number: pTJet, cent, eta, phi, Area, trigger, acceptance, signed pT leading + const Int_t nBinsSparse2 = 9; + Int_t nBins2[nBinsSparse2] = { kMaxJets+1, 60, 8, 6, 90, 100,fNTrigger,fNAcceptance+1,20}; if(cJetBranch.Contains("RandomCone")){ nBins2[5] = 1; } - const Double_t xmin2[nBinsSparse2] = { -0.5, 0, 0,-0.9, 0, 0.}; - const Double_t xmax2[nBinsSparse2] = {kMaxJets+0.5,250, 100, 0.9, 2.*TMath::Pi(),1.0}; - fhnJetPtQA[ij] = new THnSparseF(Form("fhnJetPtQA%s",cAdd.Data()),";jet number;p_{T,jet};cent;#eta;#phi;area",nBinsSparse2,nBins2,xmin2,xmax2); + const Double_t xmin2[nBinsSparse2] = { -0.5, -50, 0,-0.9, 0, 0., -0.5, -0.5, -100}; + const Double_t xmax2[nBinsSparse2] = {kMaxJets+0.5, 250, 80, 0.9, 2.*TMath::Pi(),1.0,fNTrigger-0.5,fNAcceptance+0.5, 100}; + fhnJetPtQA[ij] = new THnSparseF(Form("fhnJetPtQA%s",cAdd.Data()),";jet number;p_{T,jet};cent;#eta;#phi;area;trigger;acceptance bin;signed pt leading",nBinsSparse2,nBins2,xmin2,xmax2); + // fhnJetPt[ij]->SetBinEdges(5,binArrayArea); fHistList->Add(fhnJetPtQA[ij]); - // Bins:track number pTtrack, cent, mult, RP. total bins = 224 k - const Int_t nBinsSparse3 = 5; - const Int_t nBins3[nBinsSparse3] = { 2, 100, 10, 20, fNRPBins}; - const Double_t xmin3[nBinsSparse3] = { -0.5, 0, 0, 0, -0.5}; - const Double_t xmax3[nBinsSparse3] = { 1.5, 200, 100, 4000,fNRPBins-0.5}; + // Bins:track number pTtrack, cent, mult, RP. charge total bins = 224 k + const Int_t nBinsSparse3 = 7; + const Int_t nRPBinsSparse3 = 3; + const Int_t nBins3[nBinsSparse3] = { 2, 100, 10, 1, nRPBinsSparse3,fNTrigger,2}; + const Double_t xmin3[nBinsSparse3] = { -0.5, 0, 0, 0, -0.5,-0.5,-1.5}; + const Double_t xmax3[nBinsSparse3] = { 1.5, 200, 100, 4000,nRPBinsSparse3-0.5,fNTrigger-0.5,1.5}; - // change the binning ot the pT axis: + // change the binning of the pT axis: Double_t *xPt3 = new Double_t[nBins3[1]+1]; xPt3[0] = 0.; for(int i = 1; i<=nBins3[1];i++){ if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.05; // 1 - 40 else if(xPt3[i-1]<4)xPt3[i] = xPt3[i-1] + 0.2; // 41 - 50 else if(xPt3[i-1]<10)xPt3[i] = xPt3[i-1] + 0.5; // 50 - 62 - else if(xPt3[i-1]<20)xPt3[i] = xPt3[i-1] + 1.; // 62 - 72 - else if(xPt3[i-1]<30)xPt3[i] = xPt3[i-1] + 2.5; // 74 - 78 - else xPt3[i] = xPt3[i-1] + 5.; // 78 - 100 = 140 + else if(xPt3[i-1]<15)xPt3[i] = xPt3[i-1] + 1.; // 62 - 67 + else if(xPt3[i-1]<20)xPt3[i] = xPt3[i-1] + 2.; // 67 - 72 + else if(xPt3[i-1]<60)xPt3[i] = xPt3[i-1] + 5; // 72 - 76 + else xPt3[i] = xPt3[i-1] + 10; // 76 - 100 = 140 } - fhnTrackPt[ij] = new THnSparseF(Form("fhnTrackPt%s",cAdd.Data()),";track number;p_{T};cent;#tracks;RP",nBinsSparse3,nBins3,xmin3,xmax3); - fhnTrackPt[ij]->SetBinEdges(1,xPt3); - fHistList->Add(fhnTrackPt[ij]); - delete [] xPt3; - - // Track QA bins track nr, pTrack, cent, eta, phi bins 5.4 M - const Int_t nBinsSparse4 = 5; - const Int_t nBins4[nBinsSparse4] = { 2, 50, 10, 20, 360}; - const Double_t xmin4[nBinsSparse4] = { -0.5, 0, 0, -1.0, 0.}; - const Double_t xmax4[nBinsSparse4] = { 1.5,150, 100, 1.0,2.*TMath::Pi()}; - - // change the binning ot the pT axis: - Double_t *xPt4 = new Double_t[nBins4[1]+1]; - xPt4[0] = 0.; - for(int i = 1; i<=nBins4[1];i++){ - if(xPt4[i-1]<2)xPt4[i] = xPt4[i-1] + 0.1; - else if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 0.5; - else if(xPt4[i-1]<20)xPt4[i] = xPt4[i-1] + 1.; - else if(xPt4[i-1]<30)xPt4[i] = xPt4[i-1] + 2.5; - else xPt4[i] = xPt4[i-1] + 5.; - } - fhnTrackPtQA[ij] = new THnSparseF(Form("fhnTrackPtQA%s",cAdd.Data()),";track number;p_{T};cent;#eta;#phi",nBinsSparse4,nBins4,xmin4,xmax4); - fhnTrackPtQA[ij]->SetBinEdges(1,xPt4); - fHistList->Add(fhnTrackPtQA[ij]); - delete [] xPt4; - - for(int i = 0;i <= kMaxJets;++i){ - fh1PtIn[ij][i] = new TH1F(Form("fh1Pt%sIn_j%d",cAdd.Data(),i),Form("%s p_T input ;p_{T}",cAdd.Data()),nBinPt,binLimitsPt); - fHistList->Add(fh1PtIn[ij][i]); - - - if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2); - fh2LTrackPtJetPt[ij][i] = new TH2F(Form("fh2LTrackPtJetPt%s_j%d",cAdd.Data(),i), - Form("pt of leadin track within a jet vs jet %s;p_{T,lead in jet};p_{T.jet};", - cAdd.Data()), - 200,0.,200.,nBinPt,binLimitsPt); - fHistList->Add(fh2LTrackPtJetPt[ij][i]); - } + fhnTrackPt[ij] = new THnSparseF(Form("fhnTrackPt%s",cAdd.Data()),";track number;p_{T};cent;#tracks;RP;trigger",nBinsSparse3,nBins3,xmin3,xmax3); + fhnTrackPt[ij]->SetBinEdges(1,xPt3); + fHistList->Add(fhnTrackPt[ij]); + delete [] xPt3; + + // Track QA bins track nr, pTrack, cent, eta, phi bins 5.4 M + const Int_t nBinsSparse4 = 6; + const Int_t nBins4[nBinsSparse4] = { 2, 50, 10, 20, 180,2}; + const Double_t xmin4[nBinsSparse4] = { -0.5, 0, 0, -1.0, 0.,-1.5}; + const Double_t xmax4[nBinsSparse4] = { 1.5,150, 100, 1.0,2.*TMath::Pi(),1.5}; + + // change the binning ot the pT axis: + Double_t *xPt4 = new Double_t[nBins4[1]+1]; + xPt4[0] = 0.; + for(int i = 1; i<=nBins4[1];i++){ + if(xPt4[i-1]<2)xPt4[i] = xPt4[i-1] + 0.1; + else if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 0.5; + else if(xPt4[i-1]<20)xPt4[i] = xPt4[i-1] + 1.; + else if(xPt4[i-1]<30)xPt4[i] = xPt4[i-1] + 2.5; + else xPt4[i] = xPt4[i-1] + 5.; + } + fhnTrackPtQA[ij] = new THnSparseF(Form("fhnTrackPtQA%s",cAdd.Data()),";track number;p_{T};cent;#eta;#phi;sign",nBinsSparse4,nBins4,xmin4,xmax4); + fhnTrackPtQA[ij]->SetBinEdges(1,xPt4); + fHistList->Add(fhnTrackPtQA[ij]); + delete [] xPt4; + + for(int i = 0;i <= kMaxJets;++i){ + fh1PtIn[ij][i] = new TH1F(Form("fh1Pt%sIn_j%d",cAdd.Data(),i),Form("%s p_T input ;p_{T}",cAdd.Data()),nBinPt,binLimitsPt); + fHistList->Add(fh1PtIn[ij][i]); + + + if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2); + fh2LTrackPtJetPt[ij][i] = new TH2F(Form("fh2LTrackPtJetPt%s_j%d",cAdd.Data(),i), + Form("pt of leadin track within a jet vs jet %s;p_{T,lead in jet};p_{T.jet};", + cAdd.Data()), + 200,0.,200.,nBinPt,binLimitsPt); + fHistList->Add(fh2LTrackPtJetPt[ij][i]); + } - fh1DijetMinv[ij] = new TH1F(Form("fh1Dijet%sMinv",cAdd.Data()),"Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt); - fHistList->Add(fh1DijetMinv[ij]); + fh1DijetMinv[ij] = new TH1F(Form("fh1Dijet%sMinv",cAdd.Data()),"Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt); + fHistList->Add(fh1DijetMinv[ij]); - fh2DijetDeltaPhiPt[ij] = new TH2F(Form("fh2Dijet%sDeltaPhiPt",cAdd.Data()),"Difference in the azimuthal angle;#Delta#phi;p_{T,2};Entries",180,0.,TMath::Pi(),nBinPt,binLimitsPt); - fHistList->Add(fh2DijetDeltaPhiPt[ij]); + fh2DijetDeltaPhiPt[ij] = new TH2F(Form("fh2Dijet%sDeltaPhiPt",cAdd.Data()),"Difference in the azimuthal angle;#Delta#phi;p_{T,2};Entries",180,0.,TMath::Pi(),nBinPt,binLimitsPt); + fHistList->Add(fh2DijetDeltaPhiPt[ij]); - fh2DijetAsymPt[ij] = new TH2F(Form("fh2Dijet%sAsym",cAdd.Data()),"Pt asymmetry;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt); - fHistList->Add(fh2DijetAsymPt[ij]); + fh2DijetAsymPt[ij] = new TH2F(Form("fh2Dijet%sAsym",cAdd.Data()),"Pt asymmetry;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt); + fHistList->Add(fh2DijetAsymPt[ij]); - fh2DijetPt2vsPt1[ij] = new TH2F(Form("fh2Dijet%sPt2vsPt1",cAdd.Data()),"Pt2 versus Pt1;p_{T,1} (GeV/c);p_{T,2} (GeV/c)",250,0.,250.,250,0.,250.); - fHistList->Add(fh2DijetPt2vsPt1[ij]); - fh2DijetDifvsSum[ij] = new TH2F(Form("fh2Dijet%sDifvsSum",cAdd.Data()),"Pt difference vs Pt sum;p_{T,1}+p_{T,2} (GeV/c);#Deltap_{T} (GeV/c)",400,0.,400.,150,0.,150.); - fHistList->Add( fh2DijetDifvsSum[ij]); - } + fh2DijetPt2vsPt1[ij] = new TH2F(Form("fh2Dijet%sPt2vsPt1",cAdd.Data()),"Pt2 versus Pt1;p_{T,1} (GeV/c);p_{T,2} (GeV/c)",250,0.,250.,250,0.,250.); + fHistList->Add(fh2DijetPt2vsPt1[ij]); + fh2DijetDifvsSum[ij] = new TH2F(Form("fh2Dijet%sDifvsSum",cAdd.Data()),"Pt difference vs Pt sum;p_{T,1}+p_{T,2} (GeV/c);#Deltap_{T} (GeV/c)",400,0.,400.,150,0.,150.); + fHistList->Add( fh2DijetDifvsSum[ij]); + } // =========== Switch on Sumw2 for all histos =========== for (Int_t i=0; iGetEntries(); ++i) { TH1 *h1 = dynamic_cast(fHistList->At(i)); @@ -610,6 +716,14 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){ selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask); } + if(!selected){ + // no selection by the service task, we continue + if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass()); + fh1EvtSelection->Fill(1.); + PostData(1, fHistList); + return; + } + if(fEventClass>0){ selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass); } @@ -617,6 +731,7 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){ if(!selected){ // no selection by the service task, we continue if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass()); + fh1EvtSelection->Fill(2.); PostData(1, fHistList); return; } @@ -626,8 +741,8 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){ // take all other information from the aod we take the tracks from if(!aod){ - if(fUseAODTrackInput)aod = fAODIn; - else aod = fAODOut; + if(fUseAODTrackInput)aod = fAODIn; + else aod = fAODOut; } @@ -659,21 +774,21 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){ if(!aodRecJets){ if(fDebug){ - Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data()); - if(fAODIn){ - Printf("Input AOD >>>>"); - fAODIn->Print(); - } - if(fAODExtension){ - Printf("AOD Extension >>>>"); - fAODExtension->Print(); - } - if(fAODOut){ - Printf("Output AOD >>>>"); - fAODOut->Print(); - } + Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data()); + if(fAODIn){ + Printf("Input AOD >>>>"); + fAODIn->Print(); + } + if(fAODExtension){ + Printf("AOD Extension >>>>"); + fAODExtension->Print(); + } + if(fAODOut){ + Printf("Output AOD >>>>"); + fAODOut->Print(); + } } - return; + return; } TClonesArray *aodGenJets = 0; @@ -732,6 +847,44 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){ } } + AliAODVertex* primVtx = aod->GetPrimaryVertex(); + Int_t nTracksPrim = primVtx->GetNContributors(); + + + if (fDebug > 1) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim); + if(nTracksPrim < fMinNcontributors){ + if (fDebug > 1) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__); + fh1EvtSelection->Fill(3.); + PostData(1, fHistList); + return; + } + + if(TMath::Abs(primVtx->GetZ())>fMaxVertexZ){ + if (fDebug > 1) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ()); + fh1EvtSelection->Fill(4.); + PostData(1, fHistList); + return; + } + + TString primVtxName(primVtx->GetName()); + + if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){ + if (fDebug > 1) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__); + fh1EvtSelection->Fill(5.); + PostData(1, fHistList); + return; + } + + if(fRejectPileup && AliAnalysisHelperJetTasks::IsPileUp()){ + if (fDebug > 1) Printf("%s:%d SPD pileup: event REJECTED...",(char*)__FILE__,__LINE__); + fh1EvtSelection->Fill(6.); + PostData(1, fHistList); + return; + } + + if (fDebug > 1) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__); + fh1EvtSelection->Fill(0.); + // new Scheme // first fill all the pure histograms, i.e. full jets @@ -764,7 +917,7 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){ Double_t eventW = 1; Double_t ptHard = 0; Double_t nTrials = 1; // Trials for MC trigger - fh1Trials->Fill("#sum{ntrials}",fAvgTrials); + fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials); // Getting some global properties fCentrality = GetCentrality(); @@ -803,6 +956,7 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){ TList genParticles; Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec); + if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries()); nT = GetListOfTracks(&genParticles,fTrackTypeGen); if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries()); @@ -810,9 +964,11 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){ // CalculateReactionPlaneAngle(&recParticles); fRPAngle = 0; - if(fRPMethod==0)fRPAngle = aod->GetHeader()->GetEventplane(); + if(fRPMethod==0)fRPAngle = ((AliVAODHeader*)aod->GetHeader())->GetEventplane(); else if(fRPMethod==1||fRPMethod==2){ - fRPAngle = aod->GetHeader()->GetQTheta(fRPMethod); + AliAODHeader * aodheader = dynamic_cast(aod->GetHeader()); + if(!aodheader) AliFatal("Not a standard AOD"); + fRPAngle = aodheader->GetQTheta(fRPMethod); } fh1RP->Fill(fRPAngle); fh2RPCentrality->Fill(fCentrality,fRPAngle); @@ -841,12 +997,13 @@ void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){ fMultGen = genMult1; if(fMultGen<=0)fMultGen = genMult2; - Double_t var0[3] = {0,}; + Double_t var0[4] = {0,}; var0[0] = fCentrality; var0[1] = fMultRec; for(int it=0;itIsEventSelected()&fTriggerBit[it]){ var0[2] = it; + var0[3] = GetRho(recJetsList); fhnEvent->Fill(var0); } } @@ -893,39 +1050,56 @@ void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particles if(nJets<=0)return; - Float_t ptOld = 1.E+32; - Float_t pT0 = 0; Float_t pT1 = 0; Float_t phi0 = 0; Float_t phi1 = 0; Int_t ij0 = -1; Int_t ij1 = -1; - Double_t var1[7] = {0,}; // jet number;p_{T,jet};cent;# tracks;RP;area + Double_t var1[9] = {0,}; // jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leadingTrackPt + Double_t var1b[3] = {0,}; // p_{T,jet};cent;trigger; + var1[2] = fCentrality; + var1b[1] = fCentrality; var1[3] = refMult; - Double_t var2[6] = {0,}; // jet number;p_{T,jet};cent;#eta;#phi;area + + + Double_t var2[9] = {0,}; // jet number;p_{T,jet};cent;#eta;#phi;area;trigger var2[2] = fCentrality; for(int ij = 0;ij < nJets;ij++){ AliAODJet *jet = (AliAODJet*)jetsList.At(ij); Float_t ptJet = jet->Pt(); + + if(ptJet<0.150)ptJet = jet->GetPtSubtracted(0); + + var1b[0] = ptJet; + if(jet->Trigger()&fJetTriggerBestMask){ + fh1PtJetsInBest[iType]->Fill(ptJet); + for(int it = 0;it IsEventSelected()&fTriggerBit[it]){ + var1b[2] = it; + fhnJetPtBest[iType]->Fill(var1b); + } + } + } if(jet->Trigger()&fJetTriggerExcludeMask){ fh1PtJetsInRej[iType]->Fill(ptJet); + for(int it = 0;it IsEventSelected()&fTriggerBit[it]){ + var1b[2] = it; + fhnJetPtRej[iType]->Fill(var1b); + } + } continue; } fh1PtJetsIn[iType]->Fill(ptJet); - if(ptJet>ptOld){ - Printf("%s:%d Jets Type %d Not Sorted !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,iType,ij,ptJet,ij-1,ptOld); - } - ptOld = ptJet; // find the dijets assume sorting and acceptance cut... if(ij==0){ ij0 = ij; - pT0 = ptJet; phi0 = jet->Phi(); if(phi0<0)phi0+=TMath::Pi()*2.; } @@ -942,60 +1116,84 @@ void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particles } } // fill jet histos for kmax jets - - Float_t phiJet = jet->Phi(); - Float_t etaJet = jet->Eta(); - if(phiJet<0)phiJet+=TMath::Pi()*2.; - fh1TmpRho->Reset(); - if(ijFill(ptJet); - - fh1PtIn[iType][kMaxJets]->Fill(ptJet); - // fill leading jets... - AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet); - // AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList); - Int_t phiBin = GetPhiBin(phiJet-fRPAngle); - + Float_t phiJet = jet->Phi(); + Float_t etaJet = jet->Eta(); + if(phiJet<0)phiJet+=TMath::Pi()*2.; + fh1TmpRho->Reset(); + if(ijFill(ptJet); + + fh1PtIn[iType][kMaxJets]->Fill(ptJet); + // fill leading jets... + AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet); + //AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList); + Int_t phiBin = GetPhiBin(phiJet-fRPAngle); + Double_t ptLead = jet->GetPtLeading(); //pT of leading jet + AliVParticle *tt = (AliVParticle*)particlesList.At(0); + Float_t ttphi=tt->Phi(); + if(ttphi<0)ttphi+=TMath::Pi()*2.; + Float_t ttpt=tt->Pt(); + Int_t phiBintt = GetPhiBin(ttphi-fRPAngle); + Double_t dphitrigjet=RelativePhi(ttphi,phiJet); + if(fTRP==1){ + if(TMath::Abs(dphitrigjet)EffectiveAreaCharged(); + var1[7] = ttpt; + var1[8] = CheckAcceptance(phiJet,etaJet); + } + + if(fTRP==0){ var1[1] = ptJet; var1[4] = phiBin; var1[5] = jet->EffectiveAreaCharged(); - - var2[1] = ptJet; - var2[3] = etaJet; - var2[4] = phiJet; - var2[5] = jet->EffectiveAreaCharged(); - if(ijFill(leadTrack->Pt(),ptJet); - var1[0] = ij; - var2[0] = ij; - for(int it = 0;it IsEventSelected()&fTriggerBit[it]){ - var1[6] = it; - fhnJetPt[iType]->Fill(var1); - } - } - fhnJetPtQA[iType]->Fill(var2); - } - var1[0] = kMaxJets;// fill for all jets - var2[0] = kMaxJets;// fill for all jets + var1[7] = ptLead; + var1[8] = CheckAcceptance(phiJet,etaJet); + } + + //jet number;p_{T,jet};cent;#eta;#phi;area;trigger;acceptance bin;signed pt leading + var2[1] = ptJet; + var2[3] = etaJet; + var2[4] = phiJet; + var2[5] = jet->EffectiveAreaCharged(); + var2[7] = var1[8]; + var2[8] = (leadTrack?leadTrack->Charge()*ptLead:0);//pT of leading jet x charge + + if(ijFill(jet->GetPtLeading(),ptJet); + var1[0] = ij; + var2[0] = ij; for(int it = 0;it IsEventSelected()&fTriggerBit[it]){ var1[6] = it; + var2[6] = it; fhnJetPt[iType]->Fill(var1); + fhnJetPtQA[iType]->Fill(var2); } } + } + var1[0] = kMaxJets;// fill for all jets + var2[0] = kMaxJets;// fill for all jets + for(int it = 0;it IsEventSelected()&fTriggerBit[it]){ + var1[6] = it; + var2[6] = it; + fhnJetPt[iType]->Fill(var1); + fhnJetPtQA[iType]->Fill(var2); + } + } - fhnJetPtQA[iType]->Fill(var2); - if(leadTrack)fh2LTrackPtJetPt[iType][kMaxJets]->Fill(leadTrack->Pt(),ptJet); + fh2LTrackPtJetPt[iType][kMaxJets]->Fill(jet->GetPtLeading(),ptJet); - if(particlesList.GetSize()&&ijDeltaR(part); - if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet); - } - // fill the jet shapes - }// if we have particles + if(particlesList.GetSize()&&ijDeltaR(part); + if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet); + } + // fill the jet shapes + }// if we have particles }// Jet Loop @@ -1019,16 +1217,16 @@ void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particles Float_t asym = 9999; if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1); - fh2DijetAsymPt[iType]->Fill(asym,ptJet0); - fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1); - fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1); - Float_t minv = 2.*(jet0->P()*jet1->P()- - jet0->Px()*jet1->Px()- - jet0->Py()*jet1->Py()- - jet0->Pz()*jet1->Pz()); // assume mass == 0; - if(minv<0)minv=0; // prevent numerical instabilities - minv = TMath::Sqrt(minv); - fh1DijetMinv[iType]->Fill(minv); + fh2DijetAsymPt[iType]->Fill(asym,ptJet0); + fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1); + fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1); + Float_t minv = 2.*(jet0->P()*jet1->P()- + jet0->Px()*jet1->Px()- + jet0->Py()*jet1->Py()- + jet0->Pz()*jet1->Pz()); // assume mass == 0; + if(minv<0)minv=0; // prevent numerical instabilities + minv = TMath::Sqrt(minv); + fh1DijetMinv[iType]->Fill(minv); } @@ -1067,10 +1265,10 @@ void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType } // - Double_t var3[5]; // track number;p_{T};cent;#tracks;RP + Double_t var3[7]; // track number;p_{T};cent;#tracks;RP;cahrge var3[2] = fCentrality; var3[3] = refMult; - Double_t var4[5]; // track number;p_{T};cent;#eta;#phi + Double_t var4[6]; // track number;p_{T};cent;#eta;#phi;sign var4[2] = fCentrality; Int_t nTrackOver = particlesList.GetSize(); // do the same for tracks and jets @@ -1123,23 +1321,41 @@ void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType var3[0] = 1; var3[1] = tmpPt; var3[4] = phiBin; + var3[6] = tmpTrack->Charge(); var4[0] = 1; var4[1] = tmpPt; var4[3] = tmpTrack->Eta(); var4[4] = tmpPhi; + var4[5] = tmpTrack->Charge(); + + + for(int it = 0;it IsEventSelected()&fTriggerBit[it]){ + var3[0] = 1; + var3[5] = it; + fhnTrackPt[iType]->Fill(var3); + if(tmpTrack==leading){ + var3[0] = 0; + fhnTrackPt[iType]->Fill(var3); + } + } + } - fhnTrackPt[iType]->Fill(var3); - fhnTrackPtQA[iType]->Fill(var4); + + fhnTrackPtQA[iType]->Fill(var4); if(tmpTrack==leading){ - var3[0] = 0; var4[0] = 0; - fhnTrackPt[iType]->Fill(var3); fhnTrackPtQA[iType]->Fill(var4); continue; } + + + + + } fh1SumPtTrack[iType]->Fill(sumPt); delete trackIter; @@ -1179,7 +1395,9 @@ void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJ } } - Double_t container[6]; + Double_t container[8]; + Double_t containerRec[4]; + Double_t containerGen[4]; // loop over generated jets // consider the @@ -1189,18 +1407,31 @@ void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJ Double_t phiGen = genJet->Phi(); if(phiGen<0)phiGen+=TMath::Pi()*2.; Double_t etaGen = genJet->Eta(); - container[3] = ptGen; - container[4] = etaGen; - container[5] = phiGen; - fhnJetContainer->Fill(&container[3],kStep0); - if(JetSelected(genJet)){ - fhnJetContainer->Fill(&container[3],kStep1); - Int_t ir = aRecIndex[ig]; - if(ir>=0&&irFill(&container[3],kStep2); - AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir); - if(JetSelected(recJet))fhnJetContainer->Fill(&container[3],kStep4); - if(JetSelected(recJet))fhnJetContainer->Fill(&container[3],kStep3); + containerGen[0] = ptGen; + containerGen[1] = etaGen; + containerGen[2] = phiGen; + containerGen[3] = genJet->GetPtLeading(); + + fhnJetContainer->Fill(containerGen,kStep0); //all generated jets + + if(JetSelected(genJet)) + fhnJetContainer->Fill(containerGen,kStep1); // all generated jets in eta window + + Int_t ir = aRecIndex[ig]; + if(ir>=0&&irFill(containerGen,kStep2); // all generated jets with reconstructed partner + + if(JetSelected(genJet)){ + fhnJetContainer->Fill(containerGen,kStep3); // all generated jets in eta window with reconstructed partner + + if(JetSelected(recJet)) { + fhnJetContainer->Fill(containerGen,kStep4); // all generated jets in eta window with reconstructed partner in eta window + + } + containerGen[3] = recJet->GetPtLeading(); + fhnJetContainer->Fill(containerGen,kStep5); // all generated jets in eta window with reconstructed partner with leading track on reconstructed level } } }// loop over generated jets used for matching... @@ -1214,42 +1445,52 @@ void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJ Double_t ptRec = recJet->Pt(); Double_t phiRec = recJet->Phi(); if(phiRec<0)phiRec+=TMath::Pi()*2.; - // do something with dijets... + containerRec[0] = ptRec; + containerRec[1] = etaRec; + containerRec[2] = phiRec; + containerRec[3] = recJet->GetPtLeading(); + container[0] = ptRec; container[1] = etaRec; container[2] = phiRec; + container[3] = recJet->GetPtLeading(); - fhnJetContainer->Fill(container,kStep0+kMaxStep); if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__); if(JetSelected(recJet)){ - fhnJetContainer->Fill(container,kStep1+kMaxStep); + fhnJetContainer->Fill(containerRec,kStep7); //all rec jets in eta window // Fill Correlation Int_t ig = aGenIndex[ir]; if(ig>=0 && igFill(container,kStep2+kMaxStep); if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir); AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig); Double_t ptGen = genJet->Pt(); Double_t phiGen = genJet->Phi(); if(phiGen<0)phiGen+=TMath::Pi()*2.; Double_t etaGen = genJet->Eta(); - - container[3] = ptGen; - container[4] = etaGen; - container[5] = phiGen; - // - // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance - // - if(JetSelected(genJet))fhnJetContainer->Fill(container,kStep4+kMaxStep); - fhnJetContainer->Fill(container,kStep3+kMaxStep); - fhnCorrelation->Fill(container,0); - if(ptGen>0){ - Float_t delta = (ptRec-ptGen)/ptGen; - fh2RelPtFGen->Fill(ptGen,delta); - fh2PtFGen->Fill(ptGen,ptRec); - } + + containerGen[0] = ptGen; + containerGen[1] = etaGen; + containerGen[2] = phiGen; + containerGen[3] = genJet->GetPtLeading(); + + container[4] = ptGen; + container[5] = etaGen; + container[6] = phiGen; + container[7] = genJet->GetPtLeading(); + + fhnJetContainer->Fill(containerGen,kStep6); // all rec jets in eta window with generated partner + + fhnCorrelation->Fill(container); + + Float_t delta = (ptRec-ptGen)/ptGen; + fh2RelPtFGen->Fill(ptGen,delta); + fh3RelPtFGenLeadTrkPt->Fill(ptGen,delta,recJet->GetPtLeading()); + fh2PtFGen->Fill(ptGen,ptRec); + + fh2deltaPt1Pt2->Fill(ptRec-ptGen,ptRec,ptGen); + } }// loop over reconstructed jets } @@ -1262,10 +1503,11 @@ void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){ // Create the particle container for the correction framework manager and // link it // - const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi + const Int_t kNvar = 4 ; //number of variables on the grid:pt,eta, phi const Double_t kPtmin = 0.0, kPtmax = 250.; // we do not want to have empty bins at the beginning... - const Double_t kEtamin = -3.0, kEtamax = 3.0; + const Double_t kEtamin = -1.0, kEtamax = 1.0; const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi(); + const Double_t kPtLeadingTrackPtMin = 0., kPtLeadingTrackPtMax = 200.; // can we neglect migration in eta and phi? // phi should be no problem since we cover full phi and are phi symmetric @@ -1275,8 +1517,9 @@ void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){ //arrays for the number of bins in each dimension Int_t iBin[kNvar]; iBin[0] = 125; //bins in pt - iBin[1] = 1; //bins in eta - iBin[2] = 1; // bins in phi + iBin[1] = 4; //bins in eta + iBin[2] = 1; // bins in phi + iBin[3] = 10; // bins in leading track Pt //arrays for lower bounds : Double_t* binEdges[kNvar]; @@ -1288,9 +1531,19 @@ void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){ for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i; for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i; for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i; - - - fhnJetContainer = new AliTHn(Form("fahnJetContainer"),Form("AliTHn jet info"),kMaxStep*2,kNvar,iBin); + binEdges[3][0]= kPtLeadingTrackPtMin; + binEdges[3][1]= 1.; + binEdges[3][2]= 2.; + binEdges[3][3]= 3.; + binEdges[3][4]= 4.; + binEdges[3][5]= 5.; + binEdges[3][6]= 6.; + binEdges[3][7]= 8.; + binEdges[3][8]= 10.; + binEdges[3][9]= 12.; + binEdges[3][10]= kPtLeadingTrackPtMax; + + fhnJetContainer = new AliCFContainer(Form("fhnJetContainer"),Form("AliCFContainer jet info"),kMaxStep,kNvar,iBin); for (int k=0; kSetBinLimits(k,binEdges[k]); } @@ -1304,10 +1557,10 @@ void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){ thnDim[k+kNvar] = iBin[k]; } - fhnCorrelation = new AliTHn("fahnCorrelation","AliTHn with correlations",1,2*kNvar,thnDim); + fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparseF with correlations",2*kNvar,thnDim); for (int k=0; kSetBinLimits(k,binEdges[k]); - fhnCorrelation->SetBinLimits(k+kNvar,binEdges[k]); + fhnCorrelation->SetBinEdges(k,binEdges[k]); + fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]); } for(Int_t ivar = 0; ivar < kNvar; ivar++) @@ -1318,9 +1571,9 @@ void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){ void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/) { -// Terminate analysis -// - if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n"); + // Terminate analysis + // + if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n"); } @@ -1337,7 +1590,8 @@ Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){ return iCount; } for(int it = 0;it < aod->GetNumberOfTracks();++it){ - AliAODTrack *tr = aod->GetTrack(it); + AliAODTrack *tr = dynamic_cast(aod->GetTrack(it)); + if(!tr) AliFatal("Not a standard AOD"); if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue; if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue; if(tr->Pt()(InputEvent()); - else aod = AODEvent(); - if(!aod){ - return 101; - } - return aod->GetHeader()->GetCentrality(); + AliAODEvent *aod = 0; + if(fUseAODTrackInput)aod = dynamic_cast(InputEvent()); + else aod = AODEvent(); + if(!aod){ + return 101; + } + return ((AliVAODHeader*)aod->GetHeader())->GetCentrality(); } @@ -1522,38 +1776,143 @@ AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TLi return leading; } + +Double_t AliAnalysisTaskJetSpectrum2::RelativePhi(Double_t mphi,Double_t vphi){ + + if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi()); + else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi()); + if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi()); + else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi()); + double dphi = mphi-vphi; + if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi()); + else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi()); + return dphi;//dphi in [-Pi, Pi] +} + + + + + + Int_t AliAnalysisTaskJetSpectrum2::GetPhiBin(Double_t phi) { - Int_t phibin=-1; - if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;} - Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi))); - phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi())); - if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");} - return phibin; + Int_t phibin=-1; + if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;} + Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi))); + phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi())); + if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");} + return phibin; } void AliAnalysisTaskJetSpectrum2::SetNTrigger(Int_t n){ + // + // set number of trigger bins + // if(n>0){ - fNTrigger = n; - delete [] fTriggerName; - fTriggerName = new TString [fNTrigger]; - delete [] fTriggerBit;fTriggerBit = 0; - fTriggerBit = new UInt_t [fNTrigger]; + fNTrigger = n; + delete [] fTriggerName; + fTriggerName = new TString [fNTrigger]; + delete [] fTriggerBit;fTriggerBit = 0; + fTriggerBit = new UInt_t [fNTrigger]; } else{ fNTrigger = 0; } } + void AliAnalysisTaskJetSpectrum2::SetTrigger(Int_t i,UInt_t it,const char* c){ + // + // set trigger bin + // if(i0){ + fNAcceptance = n; + delete [] fAcceptancePhiMin; + delete [] fAcceptancePhiMax; + delete [] fAcceptanceEtaMin; + delete [] fAcceptanceEtaMax; + + fAcceptancePhiMin = new Float_t[fNAcceptance]; + fAcceptancePhiMax = new Float_t[fNAcceptance]; + fAcceptanceEtaMin = new Float_t[fNAcceptance]; + fAcceptanceEtaMax = new Float_t[fNAcceptance]; + } + else{ + fNTrigger = 0; + } +} + +void AliAnalysisTaskJetSpectrum2::SetAcceptance(Int_t i,Float_t phiMin,Float_t phiMax,Float_t etaMin,Float_t etaMax){ + // + // Set acceptance bins + // + + if(ifAcceptancePhiMax[i])continue; + if(etafAcceptanceEtaMax[i])continue; + return i; + } + // catch the rest + return fNAcceptance; +} + +Float_t AliAnalysisTaskJetSpectrum2::GetRho(TList &list){ + + // invert the correction + AliAODJet *jet = (AliAODJet*)list.At(0); // highest pt jet + if(!jet)return -1; + if(jet->EffectiveAreaCharged()<=0)return -1; + Float_t rho = jet->ChargedBgEnergy()/jet->EffectiveAreaCharged(); + return rho; + +} + + + AliAnalysisTaskJetSpectrum2::~AliAnalysisTaskJetSpectrum2(){ // delete [] fTriggerBit; delete [] fTriggerName; + + delete [] fAcceptancePhiMin; + delete [] fAcceptancePhiMax; + delete [] fAcceptanceEtaMin; + delete [] fAcceptanceEtaMax; + + + }