// **************************************
-// Task used for the correction of determiantion of reconstructed jet spectra
+// used for the correction of determiantion of reconstructed jet spectra
// Compares input (gen) and output (rec) jets
// *******************************************
#include <TROOT.h>
#include <TRandom.h>
+#include <TRandom3.h>
#include <TSystem.h>
#include <TInterpreter.h>
#include <TChain.h>
fAODOut(0x0),
fAODExtension(0x0),
fhnCorrelation(0x0),
- fhnCorrelationPhiZRec(0x0),
+ fhnEvent(0x0),
f1PtScale(0x0),
fBranchRec("jets"),
fBranchGen(""),
fBranchBkgRec(""),
fBranchBkgGen(""),
fNonStdFile(""),
+ fRandomizer(0x0),
fUseAODJetInput(kFALSE),
fUseAODTrackInput(kFALSE),
fUseAODMCInput(kFALSE),
fUseExternalWeightOnly(kFALSE),
fLimitGenJetEta(kFALSE),
fNMatchJets(5),
+ fNRPBins(3),
fFilterMask(0),
fEventSelectionMask(0),
fAnalysisType(0),
fTrackTypeRec(kTrackUndef),
fTrackTypeGen(kTrackUndef),
fEventClass(0),
+ fRPSubeventMethod(0),
fAvgTrials(1),
fExternalWeight(1),
fJetRecEtaWindow(0.5),
fMinJetPt(0),
fMinTrackPt(0.15),
fDeltaPhiWindow(90./180.*TMath::Pi()),
+ fCentrality(100),
+ fRPAngle(0),
fMultRec(0),
fMultGen(0),
fh1Xsec(0x0),
fh1PtHardNoW(0x0),
fh1PtHardTrials(0x0),
fh1ZVtx(0x0),
+ fh1RP(0x0),
+ fh1Centrality(0x0),
fh1TmpRho(0x0),
fh2MultRec(0x0),
fh2MultGen(0x0),
+ fh2RPSubevents(0x0),
+ fh2RPCentrality(0x0),
+ fh2RPDeltaRP(0x0),
+ fh2RPQxQy(0x0),
+ fh2RPCosDeltaRP(0x0),
fh2PtFGen(0x0),
fh2RelPtFGen(0x0),
+ fh3PhiWeights(0x0),
+ fh3RPPhiTracks(0x0),
fHistList(0x0)
{
+
+ fFlatA[0] = fFlatA[1] = 0;
+ fFlatB[0] = fFlatB[1] = 0;
+ fDeltaQxy[0] = fDeltaQxy[1] = 0;
+
for(int i = 0;i < kMaxStep*2;++i){
fhnJetContainer[i] = 0;
}
for(int ij = 0;ij <kJetTypes;++ij){
+ fFlagJetType[ij] = 1; // default = on
fh1NJets[ij] = 0;
fh1SumPtTrack[ij] = 0;
fh1PtJetsIn[ij] = 0;
fh1PtTracksIn[ij] = 0;
- fh1PtTracksLeadingIn[ij] = 0;
- fh2MultJetPt[ij] = 0;
+ fh1PtTracksInLow[ij] = 0;
fh2NJetsPt[ij] = 0;
fh2NTracksPt[ij] = 0;
- fh2LeadingTrackPtTrackPhi[ij] = 0;
+ fhnJetPt[ij] = 0;
+ fhnJetPtQA[ij] = 0;
+ fhnTrackPt[ij] = 0;
+ fhnTrackPtQA[ij] = 0;
for(int i = 0;i <= kMaxJets;++i){
- fh2PhiPt[ij][i] = 0;
- fh2EtaPt[ij][i] = 0;
- fh2AreaPt[ij][i] = 0;
- fh2EtaArea[ij][i] = 0;
- fh2PhiEta[ij][i] = 0;
fh2LTrackPtJetPt[ij][i] = 0;
fh1PtIn[ij][i] = 0;
- fh2RhoPt[ij][i] = 0;
- fh2PsiPt[ij][i] = 0;
}
fh1DijetMinv[ij] = 0;
fAODOut(0x0),
fAODExtension(0x0),
fhnCorrelation(0x0),
- fhnCorrelationPhiZRec(0x0),
+ fhnEvent(0x0),
f1PtScale(0x0),
fBranchRec("jets"),
fBranchGen(""),
fBranchBkgRec(""),
fBranchBkgGen(""),
fNonStdFile(""),
+ fRandomizer(0x0),
fUseAODJetInput(kFALSE),
fUseAODTrackInput(kFALSE),
fUseAODMCInput(kFALSE),
fUseExternalWeightOnly(kFALSE),
fLimitGenJetEta(kFALSE),
fNMatchJets(5),
+ fNRPBins(3),
fFilterMask(0),
fEventSelectionMask(0),
fAnalysisType(0),
fTrackTypeRec(kTrackUndef),
fTrackTypeGen(kTrackUndef),
fEventClass(0),
+ fRPSubeventMethod(0),
fAvgTrials(1),
fExternalWeight(1),
fJetRecEtaWindow(0.5),
fMinJetPt(0),
fMinTrackPt(0.15),
fDeltaPhiWindow(90./180.*TMath::Pi()),
+ fCentrality(100),
+ fRPAngle(0),
fMultRec(0),
fMultGen(0),
fh1Xsec(0x0),
fh1PtHardNoW(0x0),
fh1PtHardTrials(0x0),
fh1ZVtx(0x0),
+ fh1RP(0x0),
+ fh1Centrality(0x0),
fh1TmpRho(0x0),
fh2MultRec(0x0),
fh2MultGen(0x0),
+ fh2RPSubevents(0x0),
+ fh2RPCentrality(0x0),
+ fh2RPDeltaRP(0x0),
+ fh2RPQxQy(0x0),
+ fh2RPCosDeltaRP(0x0),
fh2PtFGen(0x0),
fh2RelPtFGen(0x0),
+ fh3PhiWeights(0x0),
+ fh3RPPhiTracks(0x0),
fHistList(0x0)
{
+ fFlatA[0] = fFlatA[1] = 0;
+ fFlatB[0] = fFlatB[1] = 0;
+ fDeltaQxy[0] = fDeltaQxy[1] = 0;
+
for(int i = 0;i < kMaxStep*2;++i){
fhnJetContainer[i] = 0;
}
for(int ij = 0;ij <kJetTypes;++ij){
+ fFlagJetType[ij] = 1; // default = on
fh1NJets[ij] = 0;
fh1SumPtTrack[ij] = 0;
fh1PtJetsIn[ij] = 0;
fh1PtTracksIn[ij] = 0;
- fh1PtTracksLeadingIn[ij] = 0;
- fh2MultJetPt[ij] = 0;
+ fh1PtTracksInLow[ij] = 0;
fh2NJetsPt[ij] = 0;
fh2NTracksPt[ij] = 0;
- fh2LeadingTrackPtTrackPhi[ij] = 0;
+ fhnJetPt[ij] = 0;
+ fhnJetPtQA[ij] = 0;
+ fhnTrackPt[ij] = 0;
+ fhnTrackPtQA[ij] = 0;
for(int i = 0;i <= kMaxJets;++i){
- fh2PhiPt[ij][i] = 0;
- fh2EtaPt[ij][i] = 0;
- fh2AreaPt[ij][i] = 0;
- fh2EtaArea[ij][i] = 0;
- fh2PhiEta[ij][i] = 0;
fh2LTrackPtJetPt[ij][i] = 0;
fh1PtIn[ij][i] = 0;
- fh2RhoPt[ij][i] = 0;
- fh2PsiPt[ij][i] = 0;
}
fh1DijetMinv[ij] = 0;
void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
{
-
+
// Connect the AOD
if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
OpenFile(1);
if(!fHistList)fHistList = new TList();
+ PostData(1, fHistList); // post data in any case once
+
+ if(!fRandomizer)fRandomizer = new TRandom3(0);
+
fHistList->SetOwner(kTRUE);
Bool_t oldStatus = TH1::AddDirectoryStatus();
TH1::AddDirectory(kFALSE);
+
+
+ // event npsparse cent, mult
+ const Int_t nBinsSparse0 = 2;
+ const Int_t nBins0[nBinsSparse0] = { 100, 500};
+ const Double_t xmin0[nBinsSparse0] = { 0, 0};
+ const Double_t xmax0[nBinsSparse0] = { 100,5000};
+
+
+ fhnEvent = new THnSparseF("fhnEvent",";cent;mult",nBinsSparse0,
+ nBins0,xmin0,xmax0);
+ fHistList->Add(fhnEvent);
+
+ /*
MakeJetContainer();
fHistList->Add(fhnCorrelation);
- fHistList->Add(fhnCorrelationPhiZRec);
for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
-
+ */
//
// Histogram
- const Int_t nBinPt = 320;
+
+
+ const Int_t nBinPt = 100;
Double_t binLimitsPt[nBinPt+1];
for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
if(iPt == 0){
binLimitsPt[iPt] = 0.0;
}
else {// 1.0
- binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
+ binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
}
}
const Int_t nBinPhi = 90;
fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
- fHistList->Add(fh1Xsec);
+ fHistList->Add(fh1Xsec);
fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
fHistList->Add(fh1Trials);
fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
fHistList->Add(fh1ZVtx);
- fh2MultRec = new TH2F("fh2MultRec","multiplicity rec;# tracks;# jetrefs",400,-0.5,4000,400,0.,4000);
+
+ fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
+ fHistList->Add(fh1RP);
+
+ fh1Centrality = new TH1F("fh1Centrality","cent;cent (%)",101,-0.5,100.5);
+ fHistList->Add(fh1Centrality);
+
+ fh2MultRec = new TH2F("fh2MultRec","multiplicity rec;# tracks;# jetrefs",500,0,5000,500,0.,5000);
fHistList->Add(fh2MultRec);
- fh2MultGen = new TH2F("fh2MultGen","multiplicity gen;# tracks;# jetrefs",400,-0.5,4000,400,0.,4000);
+ fh2MultGen = new TH2F("fh2MultGen","multiplicity gen;# tracks;# jetrefs",400,0,5000,500,0.,5000);
fHistList->Add(fh2MultGen);
- 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);
+ fh2RPSubevents = new TH2F("fh2RPSubevents" ,"Reaction Plane Angle" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
+ fHistList->Add( fh2RPSubevents);
- 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);
+ fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle" , 20, 0.,100., 180, 0, TMath::Pi());
+ fHistList->Add(fh2RPCentrality);
- for(int ij = 0;ij <kJetTypes;++ij){
- TString cAdd = "";
- TString cJetBranch = "";
- if(ij==kJetRec){
- cAdd = "Rec";
- cJetBranch = fBranchRec.Data();
- }
- else if (ij==kJetGen){
- cAdd = "Gen";
- cJetBranch = fBranchGen.Data();
- }
- else if (ij==kJetRecFull){
- cAdd = "RecFull";
- cJetBranch = fBranchRec.Data();
- }
- else if (ij==kJetGenFull){
- cAdd = "GenFull";
- cJetBranch = fBranchGen.Data();
- }
+ fh2RPDeltaRP = new TH2F("fh2DeltaRP" ,"Delta Reaction Plane Angle" , 100, -TMath::Pi()/2, TMath::Pi()/2,20,0.,100.0);
+ fHistList->Add(fh2RPDeltaRP);
- fh1NJets[ij] =new TH1F(Form("fh1N%sJets",cAdd.Data()),Form("N %s jets",cAdd.Data()),50,-0.5,49.5);
- fHistList->Add(fh1NJets[ij]);
+ fh2RPQxQy = new TH2F("fh2RPQxQy" ,"" , 100, -100,100,100,-100,100);
+ fHistList->Add(fh2RPQxQy);
- fh1PtJetsIn[ij] = new TH1F(Form("fh1PtJets%sIn",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
- fHistList->Add(fh1PtJetsIn[ij]);
+ fh2RPCosDeltaRP = new TH2F("fh2RPCosDeltaRP" ,"" , 20, 0.001,100.001,100,-1,1);
+ fHistList->Add(fh2RPCosDeltaRP);
- 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]);
- fh1PtTracksLeadingIn[ij] = new TH1F(Form("fh1PtTracksLeading%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
- fHistList->Add(fh1PtTracksLeadingIn[ij]);
+ 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);
- fh1SumPtTrack[ij] = new TH1F(Form("fh1SumPtTrack%s",cAdd.Data()),Form("Sum %s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
- fHistList->Add(fh1SumPtTrack[ij]);
+ 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);
- fh2MultJetPt[ij] = new TH2F(Form("fh2MultJetPt%s",cAdd.Data()),Form("%s jets p_T;# tracks;;p_{T} (GeV/c)",cAdd.Data()),400,0,4000,nBinPt,binLimitsPt);
- fHistList->Add(fh2MultJetPt[ij]);
+ fh3RPPhiTracks = new TH3F("fh3RPPhiTracks","Phi Tracks Pt Centrality", 10, 0.,100.,20,-5,5,200, 0, 2*TMath::Pi());
+ fHistList->Add(fh3RPPhiTracks);
+
+ for(int ij = 0;ij <kJetTypes;++ij){
+ TString cAdd = "";
+ TString cJetBranch = "";
+ if(ij==kJetRec){
+ cAdd = "Rec";
+ cJetBranch = fBranchRec.Data();
+ }
+ else if (ij==kJetGen){
+ cAdd = "Gen";
+ cJetBranch = fBranchGen.Data();
+ }
+ else if (ij==kJetRecFull){
+ cAdd = "RecFull";
+ cJetBranch = fBranchRec.Data();
+ }
+ else if (ij==kJetGenFull){
+ cAdd = "GenFull";
+ cJetBranch = fBranchGen.Data();
+ }
- fh2NJetsPt[ij] = new TH2F(Form("fh2N%sJetsPt",cAdd.Data()),Form("Number of %s jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",cAdd.Data()),nBinPt,binLimitsPt,50,-0.5,49.5);
- fHistList->Add(fh2NJetsPt[ij]);
+ if(cJetBranch.Length()==0)fFlagJetType[ij] = 0;
+ if(!fFlagJetType[ij])continue;
- fh2NTracksPt[ij] = new TH2F(Form("fh2N%sTracksPt",cAdd.Data()),Form("Number of %s tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",cAdd.Data()),nBinPt,binLimitsPt,1000,0.,4000);
- fHistList->Add(fh2NTracksPt[ij]);
+ fh1NJets[ij] =new TH1F(Form("fh1N%sJets",cAdd.Data()),Form("N %s jets",cAdd.Data()),50,-0.5,49.5);
+ fHistList->Add(fh1NJets[ij]);
+
+ fh1PtJetsIn[ij] = new TH1F(Form("fh1PtJets%sIn",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
+ fHistList->Add(fh1PtJetsIn[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]);
+
+ fh1PtTracksInLow[ij] = new TH1F(Form("fh1PtTracks%sInLow",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),50,0.,5.);
+ fHistList->Add(fh1PtTracksInLow[ij]);
+
+ fh1SumPtTrack[ij] = new TH1F(Form("fh1SumPtTrack%s",cAdd.Data()),Form("Sum %s track p_T;p_{T} (GeV/c)",cAdd.Data()),1000,0.,3000.);
+ fHistList->Add(fh1SumPtTrack[ij]);
+
+ fh2NJetsPt[ij] = new TH2F(Form("fh2N%sJetsPt",cAdd.Data()),Form("Number of %s jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",cAdd.Data()),nBinPt,binLimitsPt,50,-0.5,49.5);
+ fHistList->Add(fh2NJetsPt[ij]);
+
+ fh2NTracksPt[ij] = new TH2F(Form("fh2N%sTracksPt",cAdd.Data()),Form("Number of %s tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",cAdd.Data()),nBinPt,binLimitsPt,1000,0.,5000);
+ fHistList->Add(fh2NTracksPt[ij]);
+
- fh2LeadingTrackPtTrackPhi[ij] = new TH2F(Form("fh2Leading%sTrackPtTrackPhi",cAdd.Data()),Form("phi of leading %s track;p_{T};#phi;",cAdd.Data()),
- nBinPt,binLimitsPt,nBinPhi,binLimitsPhi);
- fHistList->Add(fh2LeadingTrackPtTrackPhi[ij]);
+ // Bins: Jet number: pTJet, cent, mult, RP, Area. total bins = 4.5M
+ const Int_t nBinsSparse1 = 6;
+ const Int_t nBins1[nBinsSparse1] = { kMaxJets+1,100, 10, 25, fNRPBins, 10};
+ const Double_t xmin1[nBinsSparse1] = { -0.5, 0, 0, 0, -0.5, 0.};
+ const Double_t xmax1[nBinsSparse1] = {kMaxJets+0.5,200,100,5000,fNRPBins-0.5,1.0};
+
+ fhnJetPt[ij] = new THnSparseF(Form("fhnJetPt%s",cAdd.Data()),";jet number;p_{T,jet};cent;# tracks;RP;area",nBinsSparse1,nBins1,xmin1,xmax1);
+ fHistList->Add(fhnJetPt[ij]);
+
+ // Bins: Jet number: pTJet, cent, eta, phi, Area. total bins = 9.72 M
+ const Int_t nBinsSparse2 = 6;
+ const Int_t nBins2[nBinsSparse2] = { kMaxJets+1, 20, 10, 18, 90, 10};
+ const Double_t xmin2[nBinsSparse2] = { -0.5, 0, 0,-0.9, 0, 0.};
+ const Double_t xmax2[nBinsSparse2] = {kMaxJets+0.5,100, 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);
+ 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, 75, 10, 25, fNRPBins};
+ const Double_t xmin3[nBinsSparse3] = { -0.5, 0, 0, 0, -0.5};
+ const Double_t xmax3[nBinsSparse3] = { 1.5, 200, 100, 5000,fNRPBins-0.5};
+
+ // change the binning ot 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.1;
+ else if(xPt3[i-1]<10)xPt3[i] = xPt3[i-1] + 0.5;
+ else if(xPt3[i-1]<20)xPt3[i] = xPt3[i-1] + 1.;
+ else if(xPt3[i-1]<40)xPt3[i] = xPt3[i-1] + 2.;
+ else xPt3[i] = xPt3[i-1] + 5.;
+ }
+
+ 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,75, 10, 20, 180};
+ const Double_t xmin4[nBinsSparse4] = { -0.5, 0, 0, -1.0, 0.};
+ const Double_t xmax4[nBinsSparse4] = { 1.5,200, 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]<40)xPt4[i] = xPt4[i-1] + 2.;
+ 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]);
- fh2RhoPt[ij][i] = new TH2F(Form("fh2RhoPt%s_j%d",cAdd.Data(),i),Form("jet shape rho for %s jets;r;p_{T};",cAdd.Data()),
- 40,0.,2.,nBinPt,binLimitsPt);
- fHistList->Add(fh2RhoPt[ij][i]);
- if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
- fh2PsiPt[ij][i] = new TH2F(Form("fh2PsiPt%s_j%d",cAdd.Data(),i),Form("jet shape #psi for %s jets;r;p_{T};",cAdd.Data()),
- 40,0.,2.,nBinPt,binLimitsPt);
- fHistList->Add(fh2PsiPt[ij][i]);
- fh2PhiPt[ij][i] = new TH2F(Form("fh2PhiPt%s_j%d",cAdd.Data(),i),Form("pt vs phi %s;#phi;p_{T};",cAdd.Data()),
- nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
- fHistList->Add(fh2PhiPt[ij][i]);
- fh2EtaPt[ij][i] = new TH2F(Form("fh2EtaPt%s_j%d",cAdd.Data(),i),Form("pt vs eta %s;#eta;p_{T};",cAdd.Data()),
- 50,-1.,1.,nBinPt,binLimitsPt);
- fHistList->Add(fh2EtaPt[ij][i]);
- fh2AreaPt[ij][i] = new TH2F(Form("fh2AreaPt%s_j%d",cAdd.Data(),i),
- Form("pt vs area %s;area;p_{T};",
- cAdd.Data()),
- 50,0.,1.,nBinPt,binLimitsPt);
- fHistList->Add(fh2AreaPt[ij][i]);
- fh2EtaArea[ij][i] = new TH2F(Form("fh2EtaArea%s_j%d",cAdd.Data(),i),
- Form("area vs eta %s;#eta;area;",
- cAdd.Data()),
- 50,-1.,1.,50,0,1.);
- fHistList->Add(fh2EtaArea[ij][i]);
-
- fh2PhiEta[ij][i] = new TH2F(Form("fh2PhiEta%s_j%d",cAdd.Data(),i),Form("phi vs eta %s ;phi;p_{T};",cAdd.Data()),
- nBinPhi,binLimitsPhi,20,-1.,1.);
- fHistList->Add(fh2PhiEta[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()),
if(!aodRecJets){
- Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
- return;
+ 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();
+ }
+ }
+ return;
}
TClonesArray *aodGenJets = 0;
// second fill the correlation histograms, here we limit to the n-th leading jet in case of the reconstructed
+
TList genJetsList; // full acceptance
TList genJetsListCut; // acceptance cut
TList recJetsList; // full acceptance
Double_t nTrials = 1; // Trials for MC trigger
fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
+ // Getting some global properties
+ fCentrality = GetCentrality();
+ fh1Centrality->Fill(fCentrality);
+
+
if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
// this is the part we only use when we have MC information
AliMCEvent* mcEvent = MCEvent();
nT = GetListOfTracks(&genParticles,fTrackTypeGen);
if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
+ CalculateReactionPlaneAngle(&recParticles);
+
// Event control and counting ...
// MC
fh1PtHard->Fill(ptHard,eventW);
fMultGen = genMult1;
if(fMultGen<=0)fMultGen = genMult2;
+ Double_t var0[2] = {0,};
+ var0[0] = fCentrality;
+ var0[1] = fMultRec;
+ fhnEvent->Fill(var0);
+
// the loops for rec and gen should be indentical... pass it to a separate
// function ...
// Jet Loop
// Here follows the jet matching part
// delegate to a separated method?
- FillMatchHistos(recJetsList,genJetsList);
-
+ // FillMatchHistos(recJetsList,genJetsList);
+
if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
PostData(1, fHistList);
}
if(iType>=kJetTypes){
return;
}
+ if(!fFlagJetType[iType])return;
Int_t refMult = fMultRec;
if(iType==kJetGen||iType==kJetGenFull){
Int_t ij0 = -1;
Int_t ij1 = -1;
+ Double_t var1[6] = {0,}; // jet number;p_{T,jet};cent;# tracks;RP;area
+ var1[2] = fCentrality;
+ var1[3] = refMult;
+
+ Double_t var2[6] = {0,}; // jet number;p_{T,jet};cent;#eta;#phi;area
+ var2[2] = fCentrality;
for(int ij = 0;ij < nJets;ij++){
AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
Float_t ptJet = jet->Pt();
fh1PtJetsIn[iType]->Fill(ptJet);
- fh2MultJetPt[iType]->Fill(refMult,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);
}
if(phiJet<0)phiJet+=TMath::Pi()*2.;
fh1TmpRho->Reset();
if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
+
fh1PtIn[iType][kMaxJets]->Fill(ptJet);
// fill leading jets...
AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
// AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
- if(ptJet>10){
- if(ij<kMaxJets){
- fh2PhiEta[iType][ij]->Fill(phiJet,etaJet);
- fh2AreaPt[iType][ij]->Fill(jet->EffectiveAreaCharged(),ptJet);
- fh2EtaArea[iType][ij]->Fill(etaJet,jet->EffectiveAreaCharged());
- }
- fh2PhiEta[iType][kMaxJets]->Fill(phiJet,etaJet);
- fh2AreaPt[iType][kMaxJets]->Fill(jet->EffectiveAreaCharged(),ptJet);
- fh2EtaArea[iType][kMaxJets]->Fill(etaJet,jet->EffectiveAreaCharged());
- }
+ Int_t phiBin = GetPhiBin(phiJet-fRPAngle);
+
+ var1[1] = ptJet;
+ var1[4] = phiBin;
+ var1[5] = jet->EffectiveAreaCharged();
+
+ var2[1] = ptJet;
+ var2[3] = etaJet;
+ var2[4] = phiJet;
+ var2[5] = jet->EffectiveAreaCharged();
if(ij<kMaxJets){
- fh2PhiPt[iType][ij]->Fill(phiJet,ptJet);
- fh2EtaPt[iType][ij]->Fill(etaJet,ptJet);
if(leadTrack)fh2LTrackPtJetPt[iType][ij]->Fill(leadTrack->Pt(),ptJet);
+ var1[0] = ij;
+ var2[0] = ij;
+ fhnJetPt[iType]->Fill(var1);
+ fhnJetPtQA[iType]->Fill(var2);
}
- fh2PhiPt[iType][kMaxJets]->Fill(phiJet,ptJet);
- fh2EtaPt[iType][kMaxJets]->Fill(etaJet,ptJet);
+ var1[0] = kMaxJets;// fill for all jets
+ var2[0] = kMaxJets;// fill for all jets
+ fhnJetPt[iType]->Fill(var1);
+ fhnJetPtQA[iType]->Fill(var2);
if(leadTrack)fh2LTrackPtJetPt[iType][kMaxJets]->Fill(leadTrack->Pt(),ptJet);
if(particlesList.GetSize()&&ij<kMaxJets){
if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
}
// fill the jet shapes
- Float_t rhoSum = 0;
- for(int ibx = 1;ibx<=fh2RhoPt[iType][ij]->GetNbinsX();ibx++){
- Float_t r = fh2RhoPt[iType][ij]->GetXaxis()->GetBinCenter(ibx);
- Float_t rho = fh1TmpRho->GetBinContent(ibx);
- rhoSum += rho;
- fh2RhoPt[iType][ij]->Fill(r,ptJet,rho);
- fh2PsiPt[iType][ij]->Fill(r,ptJet,rhoSum);
- }
}// if we have particles
}// Jet Loop
void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
+ if(fFlagJetType[iType]<=0)return;
+ Int_t refMult = fMultRec;
+ if(iType==kJetGen||iType==kJetGenFull){
+ refMult = fMultGen;
+
+ }
+
+ //
+ Double_t var3[5]; // track number;p_{T};cent;#tracks;RP
+ var3[2] = fCentrality;
+ var3[3] = refMult;
+ Double_t var4[5]; // track number;p_{T};cent;#eta;#phi
+ var4[2] = fCentrality;
Int_t nTrackOver = particlesList.GetSize();
// do the same for tracks and jets
if(nTrackOver>0){
while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
Float_t tmpPt = tmpTrack->Pt();
fh1PtTracksIn[iType]->Fill(tmpPt);
+ fh1PtTracksInLow[iType]->Fill(tmpPt);
+
sumPt += tmpPt;
Float_t tmpPhi = tmpTrack->Phi();
if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
+ Int_t phiBin = GetPhiBin(tmpPhi-fRPAngle);
+ var3[0] = 1;
+ var3[1] = tmpPt;
+ var3[4] = phiBin;
+
+ var4[0] = 1;
+ var4[1] = tmpPt;
+ var4[3] = tmpTrack->Eta();
+ var4[4] = tmpPhi;
+
+
+ fhnTrackPt[iType]->Fill(var3);
+ fhnTrackPtQA[iType]->Fill(var4);
+
if(tmpTrack==leading){
- fh1PtTracksLeadingIn[iType]->Fill(tmpPt);
- fh2LeadingTrackPtTrackPhi[iType]->Fill(tmpPt,tmpPhi);
+ var3[0] = 0;
+ var4[0] = 0;
+ fhnTrackPt[iType]->Fill(var3);
+ fhnTrackPtQA[iType]->Fill(var4);
continue;
}
}
}
Double_t container[6];
- Double_t containerPhiZ[6];
// loop over generated jets
// consider the
container[0] = ptRec;
container[1] = etaRec;
container[2] = phiRec;
- containerPhiZ[0] = ptRec;
- containerPhiZ[1] = phiRec;
fhnJetContainer[kStep0+kMaxStep]->Fill(container);
if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
container[3] = ptGen;
container[4] = etaGen;
container[5] = phiGen;
- containerPhiZ[3] = ptGen;
//
// we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
//
fh2RelPtFGen->Fill(ptGen,delta);
fh2PtFGen->Fill(ptGen,ptRec);
}
- if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
}
- else{
- containerPhiZ[3] = 0;
- if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
- }
}// loop over reconstructed jets
}
if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
// link it
//
const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
- const Double_t kPtmin = 0.0, kPtmax = 320.; // we do not want to have empty bins at the beginning...
+ 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 kPhimin = 0., kPhimax = 2. * TMath::Pi();
- const Double_t kZmin = 0., kZmax = 1;
// can we neglect migration in eta and phi?
// phi should be no problem since we cover full phi and are phi symmetric
//arrays for the number of bins in each dimension
Int_t iBin[kNvar];
- iBin[0] = 320; //bins in pt
+ iBin[0] = 125; //bins in pt
iBin[1] = 1; //bins in eta
iBin[2] = 1; // bins in phi
}
fhnCorrelation->Sumw2();
- // for second correlation histogram
-
-
- const Int_t kNvarPhiZ = 4;
- //arrays for the number of bins in each dimension
- Int_t iBinPhiZ[kNvarPhiZ];
- iBinPhiZ[0] = 80; //bins in pt
- iBinPhiZ[1] = 72; //bins in phi
- iBinPhiZ[2] = 20; // bins in Z
- iBinPhiZ[3] = 80; //bins in ptgen
-
- //arrays for lower bounds :
- Double_t* binEdgesPhiZ[kNvarPhiZ];
- for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
- binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
-
- for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
- for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
- for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
- for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
-
- fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
- for (int k=0; k<kNvarPhiZ; k++) {
- fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
- }
- fhnCorrelationPhiZRec->Sumw2();
-
-
- // Add a histogram for Fake jets
for(Int_t ivar = 0; ivar < kNvar; ivar++)
delete [] binEdges[ivar];
- for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
- delete [] binEdgesPhiZ[ivar];
}
}
+Float_t AliAnalysisTaskJetSpectrum2::GetCentrality(){
+ AliAODEvent *aod = 0;
+ if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
+ else aod = AODEvent();
+ if(!aod){
+ return 101;
+ }
+ return aod->GetHeader()->GetCentrality();
+}
+
+
Bool_t AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
Bool_t selected = false;
}
return leading;
}
+
+Bool_t AliAnalysisTaskJetSpectrum2::CalculateReactionPlaneAngle(const TList *trackList)
+{
+
+ if(!trackList)return kFALSE;
+ fRPAngle=0;
+
+ // need to get this info from elsewhere??
+
+ Double_t fPsiRP =0,fDeltaPsiRP = 0;
+
+
+
+ TVector2 mQ,mQ1,mQ2;
+ Float_t mQx= fDeltaQxy[0], mQy=fDeltaQxy[1];
+
+ Float_t mQx1=fDeltaQxy[0], mQy1=fDeltaQxy[1];
+ Float_t mQx2=fDeltaQxy[0], mQy2=fDeltaQxy[1];
+
+ AliVParticle *track=0x0;
+ Int_t count[3]={0,0,0};
+
+
+ for (Int_t iter=0;iter<trackList->GetEntries();iter++){
+
+ track=(AliVParticle*)trackList->At(iter);
+
+ //cuts already applied before
+ // Comment DCA not correctly implemented yet for AOD tracks
+
+ Double_t momentum;
+ if(track->Charge()>0){momentum=track->Pt();}
+ else{momentum=-track->Pt();}
+
+
+
+ // For Weighting
+ fh3RPPhiTracks->Fill(fCentrality,momentum,track->Phi());
+ count[0]++;
+
+ Double_t phiweight=GetPhiWeight(track->Phi(),momentum);
+ // Double_t phiweight=1;
+ Double_t weight=2;
+ if(track->Pt()<2){weight=track->Pt();}
+
+
+ mQx += (cos(2*track->Phi()))*weight*phiweight;
+ mQy += (sin(2*track->Phi()))*weight*phiweight;
+
+ // Make random Subevents
+
+ if(fRPSubeventMethod==0){
+ if(fRandomizer->Binomial(1,0.5)){
+ mQx1 += (cos(2*track->Phi()))*weight*phiweight;
+ mQy1 += (sin(2*track->Phi()))*weight*phiweight;
+ count[1]++;}
+ else{
+ mQx2 += (cos(2*track->Phi()))*weight*phiweight;
+ mQy2 += (sin(2*track->Phi()))*weight*phiweight;
+ count[2]++;}
+ }
+ else if(fRPSubeventMethod==1){
+ // Make eta dependent subevents
+ if(track->Eta()>0){
+ mQx1 += (cos(2*track->Phi()))*weight*phiweight;
+ mQy1 += (sin(2*track->Phi()))*weight*phiweight;
+ count[1]++;}
+ else{
+ mQx2 += (cos(2*track->Phi()))*weight*phiweight;
+ mQy2 += (sin(2*track->Phi()))*weight*phiweight;
+ count[2]++;}
+ }
+
+ }
+
+
+
+ //If no track passes the cuts, the ,Q.Phi() will return Pi and a peak at Pi/2 in the RP Angular Distribution will appear
+ if(count[0]==0||count[1]==0||count[2]==0){
+ return kFALSE;
+ }
+
+ mQ.Set(mQx,mQy);
+ mQ1.Set(mQx1,mQy1);
+ mQ2.Set(mQx2,mQy2);
+
+ // cout<<"MQ"<<mQx<<" " <<mQy<<" psi"<<endl;
+
+ fPsiRP=mQ.Phi()/2;
+
+ //Correction
+ fPsiRP+=fFlatA[0]*TMath::Cos(2*fPsiRP)+fFlatB[0]*TMath::Sin(2*fPsiRP)+fFlatA[1]*TMath::Cos(4*fPsiRP)+fFlatB[1]*TMath::Sin(4*fPsiRP);
+
+ Double_t fPsiRP1=mQ1.Phi()/2;
+ fPsiRP1+=fFlatA[0]*TMath::Cos(2*fPsiRP1)+fFlatB[0]*TMath::Sin(2*fPsiRP1)+fFlatA[1]*TMath::Cos(4*fPsiRP1)+fFlatB[1]*TMath::Sin(4*fPsiRP1);
+ Double_t fPsiRP2=mQ2.Phi()/2;
+ fPsiRP2+=fFlatA[0]*TMath::Cos(2*fPsiRP2)+fFlatB[0]*TMath::Sin(2*fPsiRP2)+fFlatA[1]*TMath::Cos(4*fPsiRP2)+fFlatB[1]*TMath::Sin(4*fPsiRP2);
+ fDeltaPsiRP=fPsiRP1-fPsiRP2;
+
+ if(fPsiRP>TMath::Pi()){fPsiRP-=TMath::Pi();}
+ if(fPsiRP<0){fPsiRP+=TMath::Pi();}
+
+ // reactionplaneangle + Pi() is the same angle
+ if(TMath::Abs(fDeltaPsiRP)>TMath::Pi()/2){
+ if(fDeltaPsiRP>0)fDeltaPsiRP-=TMath::Pi();
+ else fDeltaPsiRP+=TMath::Pi();
+ }
+
+ Double_t cos2deltaRP=TMath::Cos(2*fDeltaPsiRP);
+
+ // FillHistograms
+ fh2RPSubevents->Fill(fPsiRP1,fPsiRP2);
+ fh1RP->Fill(fPsiRP);
+ fh2RPCentrality->Fill(fCentrality,fPsiRP);
+ fh2RPDeltaRP->Fill(fDeltaPsiRP,fCentrality);
+ fh2RPQxQy->Fill(mQx,mQy);
+ fh2RPCosDeltaRP->Fill(fCentrality,cos2deltaRP);
+
+ fRPAngle=fPsiRP;
+
+ return kTRUE;
+}
+
+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;
+}
+
+
+
+Double_t AliAnalysisTaskJetSpectrum2::GetPhiWeight(Double_t phi,Double_t signedpt){
+ if(!fh3PhiWeights)return 1;
+ else return fh3PhiWeights->GetBinContent(fh3PhiWeights->GetXaxis()->FindBin(fCentrality),fh3PhiWeights->GetYaxis()->FindBin(signedpt),fh3PhiWeights->GetZaxis()->FindBin(phi));
+}
+
+ //________________________________________________________________________
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff