X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWGCF%2FEBYE%2FNetParticle%2FAliAnalysisNetParticleDistribution.cxx;h=c54431df1bc110341ee3465579f3616eeb4636d3;hb=d85f6819269447083582dbf4349395fe419ff1f5;hp=ab36cc5c93347bbd172697540fed88c1d58f38e0;hpb=3aa68b926baeb1745c1c4be38b42ed58d958fed8;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWGCF/EBYE/NetParticle/AliAnalysisNetParticleDistribution.cxx b/PWGCF/EBYE/NetParticle/AliAnalysisNetParticleDistribution.cxx index ab36cc5c933..c54431df1bc 100644 --- a/PWGCF/EBYE/NetParticle/AliAnalysisNetParticleDistribution.cxx +++ b/PWGCF/EBYE/NetParticle/AliAnalysisNetParticleDistribution.cxx @@ -2,31 +2,34 @@ #include "TMath.h" #include "TAxis.h" -#include "TSystem.h" #include "TProfile.h" -#include "TH2F.h" -#include "TH3F.h" -#include "TFile.h" -#include "TPRegexp.h" +#include "TProfile2D.h" +#include "TH2D.h" +#include "TH3D.h" #include "AliStack.h" #include "AliMCEvent.h" #include "AliMCParticle.h" +#include "AliESDEvent.h" #include "AliESDtrackCuts.h" -#include "AliESDInputHandler.h" -#include "AliESDpid.h" + #include "AliCentrality.h" #include "AliTracker.h" -#include "AliAODInputHandler.h" + #include "AliAODEvent.h" #include "AliAODTrack.h" +#include "AliAODMCParticle.h" #include "AliAnalysisNetParticleDistribution.h" using namespace std; -// Task for NetParticle checks -// Author: Jochen Thaeder +/** + * Class for NetParticle Distributions + * -- Create input for distributions + * Authors: Jochen Thaeder + * Michael Weber + */ ClassImp(AliAnalysisNetParticleDistribution) @@ -37,31 +40,19 @@ ClassImp(AliAnalysisNetParticleDistribution) */ //________________________________________________________________________ -AliAnalysisNetParticleDistribution::AliAnalysisNetParticleDistribution() : - fHelper(NULL), +AliAnalysisNetParticleDistribution::AliAnalysisNetParticleDistribution() : + AliAnalysisNetParticleBase("Dist", "Dist"), fOutList(NULL), - fESDHandler(NULL), - fPIDResponse(NULL), - fESD(NULL), - fAODHandler(NULL), - fAOD(NULL), - fIsMC(kFALSE), - fMCEvent(NULL), - fStack(NULL), - fESDTrackCuts(NULL), - fEtaMax(0.9), - fPtRange(NULL), - fAODtrackCutBit(1024), + fOrder(8), + fNNp(6), fNp(NULL), - fNCorrNp(2), - fCorrNp(NULL), - fNMCNp(5), + fNpPt(NULL), + fNMCNp(7), fMCNp(NULL), - fNControlMCNp(5), - fControlMCNp(NULL), - + fMCNpPt(NULL), + fRedFactp(NULL), fHnTrackUnCorr(NULL) { // Constructor @@ -72,19 +63,33 @@ AliAnalysisNetParticleDistribution::AliAnalysisNetParticleDistribution() : AliAnalysisNetParticleDistribution::~AliAnalysisNetParticleDistribution() { // Destructor - if (fNp) delete[] fNp; + for (Int_t ii = 0; ii < fNNp; ++ii) + if (fNp[ii]) delete[] fNp[ii]; + if (fNp) delete[] fNp; - for (Int_t ii = 0; ii < fNCorrNp; ++ii) - if (fCorrNp[ii]) delete[] fCorrNp[ii]; - if (fCorrNp) delete[] fCorrNp; + for (Int_t ii = 0; ii < 2; ++ii) { + for (Int_t kk = 0; kk < 2; ++kk) + if (fNpPt[ii][kk]) delete[] fNpPt[ii][kk]; + if (fNpPt[ii]) delete[] fNpPt[ii]; + } + if (fNpPt) delete[] fNpPt; + // ------------------------------------------- for (Int_t ii = 0; ii < fNMCNp; ++ii) if (fMCNp[ii]) delete[] fMCNp[ii]; if (fMCNp) delete[] fMCNp; - for (Int_t ii = 0; ii < fNControlMCNp; ++ii) - if (fControlMCNp[ii]) delete[] fControlMCNp[ii]; - if (fControlMCNp) delete[] fControlMCNp; + for (Int_t ii = 0; ii < 2; ++ii) { + for (Int_t kk = 0; kk < 2; ++kk) + if (fMCNpPt[ii][kk]) delete[] fMCNpPt[ii][kk]; + if (fMCNpPt[ii]) delete[] fMCNpPt[ii]; + } + if (fMCNpPt) delete[] fMCNpPt; + + // ------------------------------------------- + for (Int_t ii = 0; ii <= fOrder; ++ii) + if (fRedFactp[ii]) delete[] fRedFactp[ii]; + if (fRedFactp) delete[] fRedFactp; return; } @@ -96,240 +101,291 @@ AliAnalysisNetParticleDistribution::~AliAnalysisNetParticleDistribution() { */ //________________________________________________________________________ -Int_t AliAnalysisNetParticleDistribution::Initialize(AliAnalysisNetParticleHelper* helper, AliESDtrackCuts* cuts, Bool_t isMC, Float_t *ptRange, Float_t etaMax, Int_t trackCutBit) { - // -- Initialize - - fHelper = helper; - fESDTrackCuts = cuts; - fIsMC = isMC; - fPtRange = ptRange; - fEtaMax = etaMax; - fAODtrackCutBit = trackCutBit; - - // ------------------------------------------------------------------ - // -- N particles / N anti-particles - // ------------------------------------------------------------------ - // Np : arr[particle] - // CorrNp : arr[corrSet][particle] - // MCNp : arr[corrSet][particle] - MC - // ControlMCNp : arr[corrSet][particle] - Control MC - - fNp = new Float_t[2]; - - fCorrNp = new Float_t*[fNCorrNp]; - for (Int_t ii = 0 ; ii < fNCorrNp; ++ii) - fCorrNp[ii] = new Float_t[2]; - - fMCNp = new Float_t*[fNMCNp]; - for (Int_t ii = 0 ; ii < fNMCNp; ++ii) - fMCNp[ii] = new Float_t[2]; - - fControlMCNp = new Float_t*[fNControlMCNp]; - for (Int_t ii = 0 ; ii < fNControlMCNp; ++ii) - fControlMCNp[ii] = new Float_t[2]; +void AliAnalysisNetParticleDistribution::Process() { + // -- Process NetParticle Distributions - ResetEvent(); + // -- Fill ESD/AOD tracks + ProcessTracks(); + + // -- Fill MC truth particles (missing for AOD MW - However AliVParticle already used) + if (fIsMC) + ProcessParticles(); - return 0; + return; } +/* + * --------------------------------------------------------------------------------- + * Methods - private + * --------------------------------------------------------------------------------- + */ + //________________________________________________________________________ -void AliAnalysisNetParticleDistribution::CreateHistograms(TList* outList) { - // -- Add histograms to outlist +void AliAnalysisNetParticleDistribution::Init() { + // -- Init eventwise - fOutList = outList; - // ------------------------------------------------------------------ - // -- Create net particle histograms + // -- N particles / N anti-particles // ------------------------------------------------------------------ + // Np : arr[set][particle] + // MCNp : arr[set][particle] - MC + // Factorials : arr[order][particle] - AddHistSet("fHDist", "Uncorrected"); - AddHistSet("fHDistCorr0", "Corrected [without cross section correction]"); - AddHistSet("fHDistCorr1", "Corrected [with cross section correction]"); + // NpPt : arr[phiBin][particle][ptBins] + // MCNpPt : arr[phiBin][particle][ptBins] - MC + // FactorialsPt : arr[order][particle][ptBins] - if (fIsMC) { - AddHistSet("fHMCrapidity", "MC primary in |y| < 0.5"); - AddHistSet("fHMCptMin", "MC primary in |y| + #it{p}_{T} > 0.1"); - AddHistSet("fHMCpt", Form("MC primary in |y| < 0.5 + #it{p}_{T} [%.1f,%.1f]", fPtRange[0], fPtRange[1])); - AddHistSet("fHMCeta", Form("MC primary in |y| < 0.5 + |#eta| < %.1f", fEtaMax)); - AddHistSet("fHMCetapt", Form("MC primary in |y| < 0.5 + |#eta| < %.1f + #it{p}_{T} [%.1f,%.1f]", fEtaMax, fPtRange[0], fPtRange[1])); - - AddHistSet("fHControlMCLambdarapidity", "Control MC Lambda primary in |y| < 0.5"); - AddHistSet("fHControlMCLambdaptMin", "Control MC Lambda primary in |y| + #it{p}_{T} > 0.1"); - AddHistSet("fHControlMCLambdapt", Form("Control MC primary in |y| < 0.5 + #it{p}_{T} [%.1f,%.1f]", fPtRange[0], fPtRange[1])); - AddHistSet("fHControlMCLambdaeta", Form("Control MC primary in |y| < 0.5 + |#eta| < %.1f", fEtaMax)); - AddHistSet("fHControlMCLambdaetapt", Form("Control MC primary in |y| < 0.5 + |#eta| < %.1f + #it{p}_{T} [%.1f,%.1f]", fEtaMax, fPtRange[0], fPtRange[1])); + // set -> fNNp / fNMCNp + // order -> fOrder + // particle -> 2 {0,1} + // phiBin -> 2 {0,1} + + fNp = new Int_t*[fNNp]; + for (Int_t ii = 0 ; ii < fNNp; ++ii) + fNp[ii] = new Int_t[2]; + + fNpPt = new Int_t**[2]; + for (Int_t ii = 0 ; ii < 2; ++ii) { + fNpPt[ii] = new Int_t*[2]; + for (Int_t kk = 0 ; kk < 2; ++kk) + fNpPt[ii][kk] = new Int_t[AliAnalysisNetParticleHelper::fgkfHistNBinsPt]; } - // ------------------------------------------------------------------ - // -- Get Probe Particle Container - // ------------------------------------------------------------------ - - Double_t dCent[2] = {-0.5, 8.5}; - Int_t iCent = 9; - - Double_t dEta[2] = {-0.9, 0.9}; // -> 37 bins - Int_t iEta = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ; + // ------------------------------------------- + fMCNp = new Int_t*[fNMCNp]; + for (Int_t ii = 0 ; ii < fNMCNp; ++ii) + fMCNp[ii] = new Int_t[2]; - Double_t dRap[2] = {-0.5, 0.5}; - Int_t iRap = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ; + fMCNpPt = new Int_t**[2]; + for (Int_t ii = 0 ; ii < 2; ++ii) { + fMCNpPt[ii] = new Int_t*[2]; + for (Int_t kk = 0 ; kk < 2; ++kk) + fMCNpPt[ii][kk] = new Int_t[AliAnalysisNetParticleHelper::fgkfHistNBinsPt]; + } - Double_t dPhi[2] = {0.0, TMath::TwoPi()}; - Int_t iPhi = 42; + // ------------------------------------------- + fRedFactp = new Double_t*[fOrder+1]; + for (Int_t ii = 0 ; ii <= fOrder; ++ii) + fRedFactp[ii] = new Double_t[2]; +} - Double_t dPt[2] = {0.1, 3.0}; - Int_t iPt = Int_t((dPt[1]-dPt[0]) / 0.075); +//________________________________________________________________________ +void AliAnalysisNetParticleDistribution::Reset() { + // -- Reset eventwise - Double_t dSign[2] = {-1.5, 1.5}; - Int_t iSign = 3; + // -- Reset N particles/anti-particles + for (Int_t ii = 0; ii < fNNp; ++ii) + for (Int_t jj = 0; jj < 2; ++jj) + fNp[ii][jj] = 0; - // cent: pt: sign: eta: phi: y - Int_t binShort[6] = { iCent, iPt, iSign, iEta, iPhi, iRap}; - Double_t minShort[6] = {dCent[0], dPt[0], dSign[0], dEta[0], dPhi[0], dRap[0]}; - Double_t maxShort[6] = {dCent[1], dPt[1], dSign[1], dEta[1], dPhi[1], dRap[1]}; + for (Int_t ii = 0; ii < 2; ++ii) + for (Int_t jj = 0; jj < 2; ++jj) + for (Int_t kk = 0; kk < AliAnalysisNetParticleHelper::fgkfHistNBinsPt; ++kk) + fNpPt[ii][jj][kk] = 0; - // -- UnCorrected - fOutList->Add(new THnSparseF("fHnTrackUnCorr", "cent:pt:sign:eta:phi:y", 6, binShort, minShort, maxShort)); - fHnTrackUnCorr = static_cast(fOutList->Last()); - fHnTrackUnCorr->GetAxis(0)->SetTitle("centrality"); - fHnTrackUnCorr->GetAxis(1)->SetTitle("#it{p}_{T} (GeV/#it{c})"); - fHnTrackUnCorr->GetAxis(2)->SetTitle("sign"); - fHnTrackUnCorr->GetAxis(3)->SetTitle("#eta"); - fHnTrackUnCorr->GetAxis(4)->SetTitle("#varphi"); - fHnTrackUnCorr->GetAxis(5)->SetTitle("#it{y}"); + // -- Reset N MC particles/anti-particles + for (Int_t ii = 0; ii < fNMCNp; ++ii) + for (Int_t jj = 0; jj < 2; ++jj) + fMCNp[ii][jj] = 0; - fHelper->BinLogAxis(fHnTrackUnCorr, 1); + for (Int_t ii = 0; ii < 2; ++ii) + for (Int_t jj = 0; jj < 2; ++jj) + for (Int_t kk = 0; kk < AliAnalysisNetParticleHelper::fgkfHistNBinsPt; ++kk) + fMCNpPt[ii][jj][kk] = 0; - // ------------------------------------------------------------------ + // -- Reset reduced factorials for particles/anti-particles + for (Int_t ii = 0; ii <= fOrder; ++ii) + for (Int_t jj = 0; jj < 2; ++jj) + fRedFactp[ii][jj] = 1.; - return; } //________________________________________________________________________ -Int_t AliAnalysisNetParticleDistribution::SetupEvent(AliESDInputHandler *esdHandler, AliAODInputHandler *aodHandler, AliMCEvent *mcEvent) { - // -- Setup Event +void AliAnalysisNetParticleDistribution::CreateHistograms() { + // -- Add histograms to outlist - ResetEvent(); + // ------------------------------------------------------------------ + // -- Get Probe Particle Container + // ------------------------------------------------------------------ - // -- Get ESD objects - if(esdHandler){ - fESDHandler = esdHandler; - fPIDResponse = esdHandler->GetPIDResponse(); - fESD = fESDHandler->GetEvent(); - } + Int_t binHnUnCorr[6] = {AliAnalysisNetParticleHelper::fgkfHistNBinsCent, AliAnalysisNetParticleHelper::fgkfHistNBinsEta, // cent | eta + AliAnalysisNetParticleHelper::fgkfHistNBinsRap, AliAnalysisNetParticleHelper::fgkfHistNBinsPhi, // y | phi + AliAnalysisNetParticleHelper::fgkfHistNBinsPt, AliAnalysisNetParticleHelper::fgkfHistNBinsSign}; // pt | sign + + Double_t minHnUnCorr[6] = {AliAnalysisNetParticleHelper::fgkfHistRangeCent[0], AliAnalysisNetParticleHelper::fgkfHistRangeEta[0], + AliAnalysisNetParticleHelper::fgkfHistRangeRap[0], AliAnalysisNetParticleHelper::fgkfHistRangePhi[0], + AliAnalysisNetParticleHelper::fgkfHistRangePt[0], AliAnalysisNetParticleHelper::fgkfHistRangeSign[0]}; + + Double_t maxHnUnCorr[6] = {AliAnalysisNetParticleHelper::fgkfHistRangeCent[1], AliAnalysisNetParticleHelper::fgkfHistRangeEta[1], + AliAnalysisNetParticleHelper::fgkfHistRangeRap[1], AliAnalysisNetParticleHelper::fgkfHistRangePhi[1], + AliAnalysisNetParticleHelper::fgkfHistRangePt[1], AliAnalysisNetParticleHelper::fgkfHistRangeSign[1]}; + + // -- UnCorrected + fOutList->Add(new THnSparseD("hnTrackUnCorr", "cent:eta:y:phi:pt:sign", 6, binHnUnCorr, minHnUnCorr, maxHnUnCorr)); + fHnTrackUnCorr = static_cast(fOutList->Last()); + fHnTrackUnCorr->Sumw2(); + fHnTrackUnCorr->GetAxis(0)->SetTitle("centrality"); + fHnTrackUnCorr->GetAxis(1)->SetTitle("#eta"); + fHnTrackUnCorr->GetAxis(2)->SetTitle("#it{y}"); + fHnTrackUnCorr->GetAxis(3)->SetTitle("#varphi"); + fHnTrackUnCorr->GetAxis(4)->SetTitle("#it{p}_{T} (GeV/#it{c})"); + fHnTrackUnCorr->GetAxis(5)->SetTitle("sign"); - // -- Get AOD objects - else if(aodHandler){ - fAODHandler = aodHandler; - fPIDResponse = aodHandler->GetPIDResponse(); - fAOD = fAODHandler->GetEvent(); - } + fHelper->BinLogAxis(fHnTrackUnCorr, 4, fESDTrackCuts); - // -- Get MC objects - fMCEvent = mcEvent; - if (fMCEvent) - fStack = fMCEvent->Stack(); + for (Int_t idx = 1 ; idx <= fHnTrackUnCorr->GetAxis(4)->GetNbins(); ++idx) + printf("%02d | %f > %f < %f\n", idx, fHnTrackUnCorr->GetAxis(4)->GetBinLowEdge(idx), fHnTrackUnCorr->GetAxis(4)->GetBinCenter(idx), fHnTrackUnCorr->GetAxis(4)->GetBinUpEdge(idx)); - return 0; -} + // ------------------------------------------------------------------ + // -- Create net particle histograms + // ------------------------------------------------------------------ -//________________________________________________________________________ -void AliAnalysisNetParticleDistribution::ResetEvent() { - // -- Reset event - - // -- Reset ESD Event - fESD = NULL; + // -- Get ranges for pt and eta + Float_t etaRange[2]; + fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]); - // -- Reset AOD Event - fAOD = NULL; + Float_t ptRange[2]; + fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]); - // -- Reset MC Event - if (fIsMC) - fMCEvent = NULL; + // ------------------------------------------------------------------ - // -- Reset N particles/anti-particles - for (Int_t jj = 0; jj < 2; ++jj) - fNp[jj] = 0.; + TString sTitle(""); + sTitle = (fHelper->GetUsePID()) ? Form("|y|<%.1f", fHelper->GetRapidityMax()) : Form("|#eta|<%.1f", etaRange[1]); + + // -- centrality dependent + AddHistSetCent("Dist", Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), ptRange[0], ptRange[1])); + AddHistSetCent("DistTPC", Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired())); + AddHistSetCent("DistTOF", Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1])); + +#if USE_PHI + AddHistSetCent("Distphi", Form("%s,#it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", + sTitle.Data(), ptRange[0], ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax())); + AddHistSetCent("DistTPCphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", + sTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired(), fHelper->GetPhiMin(), fHelper->GetPhiMax())); + AddHistSetCent("DistTOFphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", + sTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax())); +#endif + + // -- centrality + pt dependent + AddHistSetCentPt("PtDist", Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), ptRange[0], ptRange[1])); +#if USE_PHI + AddHistSetCentPt("PtDistphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", + sTitle.Data(), ptRange[0], ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax())); +#endif - // -- Reset N corrected particles/anti-particles - for (Int_t ii = 0; ii < fNCorrNp; ++ii) - for (Int_t jj = 0; jj < 2; ++jj) - fCorrNp[ii][jj] = 0.; + if (fIsMC) { + TString sMCTitle(""); + sMCTitle = (fHelper->GetUsePID()) ? Form("MC primary in |y| < %.1f", fHelper->GetRapidityMax()) : Form("MC primary in |#eta| < %.1f", etaRange[1]); + + // -- centrality dependent + AddHistSetCent("MC", Form("%s", sTitle.Data())); + + AddHistSetCent("MCpt", Form("%s, #it{p}_{T} [%.1f,%.1f]", sMCTitle.Data(), ptRange[0], ptRange[1])); + AddHistSetCent("MCptTPC", Form("%s, #it{p}_{T} [%.1f,%.1f]", sMCTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired())); + AddHistSetCent("MCptTOF", Form("%s, #it{p}_{T} [%.1f,%.1f]", sMCTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1])); + +#if USE_PHI + AddHistSetCent("MCptphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", + sMCTitle.Data(), ptRange[0], ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax())); + AddHistSetCent("MCptTPCphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", + sMCTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired(), fHelper->GetPhiMin(), fHelper->GetPhiMax())); + AddHistSetCent("MCptTOFphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", + sMCTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax())); +#endif + + // -- centrality + pt dependent + AddHistSetCentPt("PtMCpt", Form("%s, #it{p}_{T} [%.1f,%.1f]", sMCTitle.Data(), ptRange[0], ptRange[1])); +#if USE_PHI + AddHistSetCentPt("PtMCptphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", + sMCTitle.Data(), ptRange[0], ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax())); +#endif + } - // -- Reset N MC particles/anti-particles - for (Int_t ii = 0; ii < fNMCNp; ++ii) - for (Int_t jj = 0; jj < 2; ++jj) - fMCNp[ii][jj] = 0.; + // ------------------------------------------------------------------ - // -- Reset N control MC particles/anti-particles - for (Int_t ii = 0; ii < fNControlMCNp; ++ii) - for (Int_t jj = 0; jj < 2; ++jj) - fControlMCNp[ii][jj] = 0.; + return; } +/* + * --------------------------------------------------------------------------------- + * Methods - private + * --------------------------------------------------------------------------------- + */ + //________________________________________________________________________ -Int_t AliAnalysisNetParticleDistribution::Process() { - // -- Process NetParticle Distributions +Int_t AliAnalysisNetParticleDistribution::ProcessTracks() { + // -- Process ESD/AOD tracks and fill histograms - // -- Fill ESD tracks - if (fESD) ProcessESDTracks(); + // -- Get ranges for AOD particles + Float_t etaRange[2]; + fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]); - // -- Fill AOD tracks - else if (fAOD) ProcessAODTracks(); + Float_t ptRange[2]; + fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]); + + // -- Track Loop + for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) { - // -- Fill MC truth particles (missing for AOD XXX) - if (fIsMC) { - ProcessStackParticles(); - ProcessStackControlParticles(); - } - - return 0; -} - -//________________________________________________________________________ -Int_t AliAnalysisNetParticleDistribution::ProcessESDTracks() { - // -- Process ESD tracks and fill histograms - - for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) { - AliESDtrack *track = fESD->GetTrack(idxTrack); + AliVTrack *track = (fESD) ? static_cast(fESD->GetTrack(idxTrack)) : static_cast(fAOD->GetTrack(idxTrack)); // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- // -- Check track cuts // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- - - // -- Check if charged track is accepted for basic parameters + + // -- Check if track is accepted for basic parameters if (!fHelper->IsTrackAcceptedBasicCharged(track)) continue; - // -- Check if accepted - if (!fESDTrackCuts->AcceptTrack(track)) + // -- Check if accepted - ESD + if (fESD && !fESDTrackCuts->AcceptTrack(dynamic_cast(track))) continue; - - // -- Check if accepted in rapidity window + + // -- Check if accepted - AOD + if (fAOD){ + AliAODTrack * trackAOD = dynamic_cast(track); + + if (!trackAOD) { + AliError("Pointer to dynamic_cast(track) = ZERO"); + continue; + } + if (!trackAOD->TestFilterBit(fAODtrackCutBit)) + continue; + + // -- Check if in pT and eta range (is done in ESDTrackCuts for ESDs) + if(!(track->Pt() > ptRange[0] && track->Pt() <= ptRange[1] && TMath::Abs(track->Eta()) <= etaRange[1])) + continue; + } + + // -- Check if accepted in rapidity window -- for identified particles + // for charged - eta check is done in the trackcuts Double_t yP; - if (!fHelper->IsTrackAcceptedRapidity(track, yP)) + if (fHelper->GetUsePID() && !fHelper->IsTrackAcceptedRapidity(track, yP)) continue; - // -- Check if accepted bt PID from TPC or TPC+TOF - Double_t pid[2]; - if (!fHelper->IsTrackAcceptedPID(track, pid)) + // -- Check if accepted with thighter DCA cuts + // -- returns kTRUE for AODs for now : MW + if (!fHelper->IsTrackAcceptedDCA(track)) continue; - // -- Check if accepted with thighter DCA cuts - if (fHelper->IsTrackAcceptedDCA(track)) + // -- Check if accepted by PID from TPC or TPC+TOF + Double_t pid[3]; + if (!fHelper->IsTrackAcceptedPID(track, pid)) continue; - + // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- // -- Fill Probe Particle // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- + if (!fHelper->GetUsePID()) + yP = track->Eta(); + Double_t aTrack[6] = { - Double_t(fHelper->GetCentralityBin()), // 0 centrality - track->Pt(), // 1 pt - track->GetSign(), // 2 sign - track->Eta(), // 3 eta - track->Phi(), // 4 phi - yP // 5 rapidity + Double_t(fCentralityBin), // 0 centrality + track->Eta(), // 1 eta + yP, // 2 rapidity + track->Phi(), // 3 phi + track->Pt(), // 4 pt + static_cast(track->Charge()) // 5 sign }; fHnTrackUnCorr->Fill(aTrack); @@ -338,116 +394,85 @@ Int_t AliAnalysisNetParticleDistribution::ProcessESDTracks() { // ------------------------------------------------------------------ // idxPart = 0 -> anti particle // idxPart = 1 -> particle + Int_t idxPart = (track->Charge() < 0) ? 0 : 1; - Int_t idxPart = (track->GetSign() < 0) ? 0 : 1; - fNp[idxPart] += 1.; - - for (Int_t ii = 0; ii < fNCorrNp; ++ii) - fCorrNp[ii][idxPart] += fHelper->GetTrackbyTrackCorrectionFactor(aTrack, ii); - - } // for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) { - - // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- - // -- Fill Particle Fluctuation Histograms - // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- + // idxPhi = 0 -> Full Acceptance + // idxPhi = 1 -> Partial Acceptance + Int_t idxPhi = 0; - // -- Uncorrected - FillHistSet("fHDist", fNp); - - // -- Corrected - for (Int_t ii = 0 ; ii < fNCorrNp; ++ii) - FillHistSet(Form("fHDistCorr%d", ii), fCorrNp[ii]); + // idxPt -> [0, nBins-1] + Int_t idxPt = static_cast(fHnTrackUnCorr->GetAxis(4))->FindBin(track->Pt())-1; - return 0; -} + // -- using pt Bins + fNpPt[idxPhi][idxPart][idxPt] += 1; -//________________________________________________________________________ -Int_t AliAnalysisNetParticleDistribution::ProcessAODTracks() { - // -- Process ESD tracks and fill histograms + // -- in pt Range + fNp[0][idxPart] += 1; - for (Int_t idxTrack = 0; idxTrack < fAOD->GetNumberOfTracks(); ++idxTrack) { - AliAODTrack *track = (AliAODTrack*)fAOD->GetTrack(idxTrack); + // -- in TPC pt Range + if (track->Pt() <= fHelper->GetMinPtForTOFRequired()) + fNp[1][idxPart] += 1; - // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- - // -- Check track cuts - // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- + // -- in TPC+TOF pt Range + if (track->Pt() > fHelper->GetMinPtForTOFRequired()) + fNp[2][idxPart] += 1; - // -- Check if charged track is accepted for basic parameters - if (!fHelper->IsTrackAcceptedBasicCharged(track)) + // -- check phi range ---------------------------------------------------------- +#if USE_PHI + if(!fHelper->IsTrackAcceptedPhi(track)) continue; + idxPhi = 1; - // -- Check if accepted - if(!track->TestFilterBit(fAODtrackCutBit)) - continue; - - // -- Check if in pT and eta range (is done in ESDTrackCuts for ESDs) - if(!(track->Pt() > fPtRange[0] && track->Pt() < fPtRange[1] && TMath::Abs(track->Eta()) <= fEtaMax)) - continue; + // -- using pt Bins + fNpPt[idxPhi][idxPart][idxPt] += 1; - // -- Check if accepted in rapidity window - Double_t yP; - if (!fHelper->IsTrackAcceptedRapidity(track, yP)) - continue; - - // -- Check if accepted bt PID from TPC or TPC+TOF - Double_t pid[2]; - if (!fHelper->IsTrackAcceptedPID(track, pid)) - continue; - - // -- Check if accepted with thighter DCA cuts XXX - // if (fHelper->IsTrackAcceptedDCA(track)) - // continue; - - // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- - // -- Fill Probe Particle - // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- - - Double_t aTrack[6] = { - Double_t(fHelper->GetCentralityBin()), // 0 centrality - track->Pt(), // 1 pt - track->Charge(), // 2 sign - track->Eta(), // 3 eta - track->Phi(), // 4 phi - yP // 5 rapidity - }; - - fHnTrackUnCorr->Fill(aTrack); - - // -- Count particle / anti-particle - // ------------------------------------------------------------------ - // idxPart = 0 -> anti particle - // idxPart = 1 -> particle + // -- in pt Range + fNp[3][idxPart] += 1; - Int_t idxPart = (track->Charge() < 0) ? 0 : 1; - fNp[idxPart] += 1.; - - for (Int_t ii = 0; ii < fNCorrNp; ++ii) - fCorrNp[ii][idxPart] += fHelper->GetTrackbyTrackCorrectionFactor(aTrack, ii); + // -- in TPC pt Range + if (track->Pt() <= fHelper->GetMinPtForTOFRequired()) + fNp[4][idxPart] += 1; + // -- in TPC+TOF pt Range + if (track->Pt() > fHelper->GetMinPtForTOFRequired()) + fNp[5][idxPart] += 1; +#endif } // for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) { // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- // -- Fill Particle Fluctuation Histograms // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- - // -- Uncorrected - FillHistSet("fHDist", fNp); - - // -- Corrected - for (Int_t ii = 0 ; ii < fNCorrNp; ++ii) - FillHistSet(Form("fHDistCorr%d", ii), fCorrNp[ii]); + FillHistSetCent("Dist", 0, kFALSE); + FillHistSetCent("DistTPC", 1, kFALSE); + FillHistSetCent("DistTOF", 2, kFALSE); + + FillHistSetCentPt("PtDist", 0, kFALSE); +#if USE_PHI + FillHistSetCent("Distphi", 3, kFALSE); + FillHistSetCent("DistTPCphi", 4, kFALSE); + FillHistSetCent("DistTOFphi", 5, kFALSE); + + FillHistSetCentPt("PtDistphi", 1, kFALSE); +#endif + return 0; } //________________________________________________________________________ -Int_t AliAnalysisNetParticleDistribution::ProcessStackParticles() { +Int_t AliAnalysisNetParticleDistribution::ProcessParticles() { // -- Process primary particles from the stack and fill histograms - - Int_t pdgCode = AliPID::ParticleCode(fHelper->GetParticleSpecies()); + Float_t etaRange[2]; + fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]); + + Float_t ptRange[2]; + fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]); + for (Int_t idxMC = 0; idxMC < fStack->GetNprimary(); ++idxMC) { - TParticle* particle = fStack->Particle(idxMC); + AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(idxMC) : NULL; + if (!particle) continue; @@ -456,229 +481,514 @@ Int_t AliAnalysisNetParticleDistribution::ProcessStackParticles() { continue; // -- Check if particle / anti-particle - if (TMath::Abs(particle->GetPdgCode()) != pdgCode) + if (fHelper->GetUsePID() && TMath::Abs(particle->PdgCode()) != fPdgCode) continue; - // -- Get particle : 0 anti-particle / 1 particle - Int_t idxPart = (particle->GetPdgCode() < 0) ? 0 : 1; - - // >> NOW only anti-particle / particle - // >> With idxPart - - // -- Check rapidity window + // -- Check rapidity window -- for identfied particles Double_t yMC; - if (!fHelper->IsParticleAcceptedRapidity(particle, yMC)) + if (fHelper->GetUsePID() && !fHelper->IsParticleAcceptedRapidity(particle, yMC)) + continue; + + // -- Check eta window -- for charged particles + if (!fHelper->GetUsePID() && TMath::Abs(particle->Eta()) > etaRange[1]) continue; - fMCNp[0][idxPart] += 1.; // -> MCrapidity - - // -- Check acceptance - if (particle->Pt() > 0.1 ) - fMCNp[1][idxPart] += 1.; // -> MCptMin - - if (particle->Pt() > fPtRange[0] && particle->Pt() <= fPtRange[1]) - fMCNp[2][idxPart] += 1.; // -> MCpt - - if (TMath::Abs(particle->Eta()) <= fEtaMax) - fMCNp[3][idxPart] += 1.; // -> MCeta + // -- Count particle / anti-particle + // ------------------------------------------------------------------ + // idxPart = 0 -> anti particle + // idxPart = 1 -> particle + Int_t idxPart = (particle->PdgCode() < 0) ? 0 : 1; + + // idxPhi = 0 -> Full Acceptance + // idxPhi = 1 -> Partial Acceptance + Int_t idxPhi = 0; + + // idxPt -> [0, nBins-1] + Int_t idxPt = static_cast(fHnTrackUnCorr->GetAxis(4))->FindBin(particle->Pt()) - 1; + + // -- MCrapidity for identfied particles + // MCeta for charged particles + fMCNp[0][idxPart] += 1.; + + // -- Check main pt window + if (!(particle->Pt() > ptRange[0] && particle->Pt() <= ptRange[1])) + continue; + + // -- using pt Bin + fMCNpPt[idxPhi][idxPart][idxPt] += 1; + + // -- in pt Range + fMCNp[1][idxPart] += 1; - if (particle->Pt() > fPtRange[0] && particle->Pt() > fPtRange[1] && TMath::Abs(particle->Eta()) <= fEtaMax) - fMCNp[4][idxPart] += 1.; // -> MCetapt + // -- in TPC pt Range + if (particle->Pt() <= fHelper->GetMinPtForTOFRequired()) + fMCNp[2][idxPart] += 1; + + // -- in TPC+TOF pt Range + if (particle->Pt() > fHelper->GetMinPtForTOFRequired()) + fMCNp[3][idxPart] += 1; + + // -- check phi range ---------------------------------------------------------- +#if USE_PHI + if(!fHelper->IsParticleAcceptedPhi(particle)) + continue; + idxPhi = 1; + + // -- using pt Bin + fMCNpPt[idxPhi][idxPart][idxPt] += 1; + // -- in pt Range + fMCNp[4][idxPart] += 1; + + // -- in TPC pt Range + if (particle->Pt() <= fHelper->GetMinPtForTOFRequired()) + fMCNp[5][idxPart] += 1; + + // -- in TPC+TOF pt Range + if (particle->Pt() > fHelper->GetMinPtForTOFRequired()) + fMCNp[6][idxPart] += 1; +#endif } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) { // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- // -- Fill Particle Fluctuation Histograms // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- - FillHistSet("fHMCrapidity", fMCNp[0]); - FillHistSet("fHMCptMin", fMCNp[1]); - FillHistSet("fHMCpt", fMCNp[2]); - FillHistSet("fHMCeta", fMCNp[3]); - FillHistSet("fHMCetapt", fMCNp[4]); + FillHistSetCent("MC", 0, kTRUE); + FillHistSetCent("MCpt", 1, kTRUE); + FillHistSetCent("MCptTPC", 2, kTRUE); + FillHistSetCent("MCptTOF", 3, kTRUE); + + FillHistSetCentPt("PtMCpt", 0, kTRUE); + +#if USE_PHI + FillHistSetCent("MCptphi", 4, kTRUE); + FillHistSetCent("MCptTPCphi", 5, kTRUE); + FillHistSetCent("MCptTOFphi", 6, kTRUE); + + FillHistSetCentPt("PtMCptphi", 1, kTRUE); +#endif return 0; } +/* + * --------------------------------------------------------------------------------- + * Helper Methods - private + * --------------------------------------------------------------------------------- + */ + //________________________________________________________________________ -Int_t AliAnalysisNetParticleDistribution::ProcessStackControlParticles() { - // -- Process primary particles from the stack and fill histograms +void AliAnalysisNetParticleDistribution::AddHistSetCent(const Char_t *name, const Char_t *title) { + // -- Add histogram sets for particle and anti-particle + // dependence : centrality - Int_t pdgCode = fHelper->GetControlParticleCode(); - Bool_t isNeutral = fHelper->IsControlParticleNeutral(); - const Char_t* name = fHelper->GetControlParticleName().Data(); + TString sName(name); + TString sTitle(title); - for (Int_t idxMC = 0; idxMC < fStack->GetNprimary(); ++idxMC) { - TParticle* particle = fStack->Particle(idxMC); - if (!particle) - continue; - - // -- Check basic MC properties -> neutral or charged physical primary - if (isNeutral) { - if (!fHelper->IsParticleAcceptedBasicNeutral(particle, idxMC)) - continue; - } - else { - if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC)) - continue; - } - - // -- Check if particle / anti-particle - if (TMath::Abs(particle->GetPdgCode()) != pdgCode) - continue; + // -- Add List + fOutList->Add(new TList); + TList *list = static_cast(fOutList->Last()); + list->SetName(Form("f%s", name)); + list->SetOwner(kTRUE); - // -- Get particle : 0 anti-particle / 1 particle - Int_t idxPart = (particle->GetPdgCode() < 0) ? 0 : 1; + // -- Get Bin Ranges + Int_t nBinsCent = AliAnalysisNetParticleHelper::fgkfHistNBinsCent; + Double_t centBinRange[] = {AliAnalysisNetParticleHelper::fgkfHistRangeCent[0], AliAnalysisNetParticleHelper::fgkfHistRangeCent[1]}; - // >> NOW only anti-particle / particle - // >> With idxPart - - // -- Check rapidity window - Double_t yMC; - if (!fHelper->IsParticleAcceptedRapidity(particle, yMC)) - continue; + // -- Create Titles + TString sNetTitle(Form("N_{%s} - N_{%s}", fHelper->GetParticleTitleLatex(1).Data(), fHelper->GetParticleTitleLatex(0).Data())); + TString sSumTitle(Form("N_{%s} + N_{%s}", fHelper->GetParticleTitleLatex(1).Data(), fHelper->GetParticleTitleLatex(0).Data())); - fControlMCNp[0][idxPart] += 1.; // -> ControlMCrapidity + // -- Add Particle / Anti-Particle Distributions + for (Int_t idxPart = 0; idxPart < 2; ++idxPart) { + list->Add(new TH2D(Form("h%s%s", name, fHelper->GetParticleName(idxPart).Data()), + Form("N_{%s} : %s;Centrality;N_{%s}", fHelper->GetParticleTitleLatex(idxPart).Data(), sTitle.Data(), fHelper->GetParticleTitleLatex(idxPart).Data()), + nBinsCent, centBinRange[0], centBinRange[1], 2601, -0.5, 2600.49)); + } // for (Int_t idxPart = 0; idxPart < 2; ++idxPart) { + + // -- Add NetParticle Distributions + list->Add(new TH2D(Form("h%sNet%s", name, fHelper->GetParticleName(1).Data()), + Form("%s : %s;Centrality;%s", sNetTitle.Data(), sTitle.Data(), sNetTitle.Data()), + nBinsCent, centBinRange[0], centBinRange[1], 601, -300.5, 300.49)); - // -- Check acceptance - if (particle->Pt() > 0.1 ) - fControlMCNp[1][idxPart] += 1.; // -> ControlMCptMin - - if (particle->Pt() > fPtRange[0] && particle->Pt() <= fPtRange[1]) - fControlMCNp[2][idxPart] += 1.; // -> ControlMCpt - - if (TMath::Abs(particle->Eta()) <= fEtaMax) - fControlMCNp[3][idxPart] += 1.; // -> ControlMCeta - - if (particle->Pt() > fPtRange[0] && particle->Pt() > fPtRange[1] && TMath::Abs(particle->Eta()) <= fEtaMax) - fControlMCNp[4][idxPart] += 1.; // -> ControlMCetapt - } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) { + // -- Add NetParticle vs SumParticle + list->Add(new TH2D(Form("h%sNet%sOverSum", name, fHelper->GetParticleName(1).Data()), + Form("(%s)/(%s) : %s;Centrality;(%s)/(%s)", sNetTitle.Data(), sSumTitle.Data(), sTitle.Data(), sNetTitle.Data(), sSumTitle.Data()), + nBinsCent, centBinRange[0], centBinRange[1], 41, -2.5, 2.49)); - // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- - // -- Fill Particle Fluctuation Histograms - // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- + // ----------------------------------------------------------------------------------------------- - FillHistSet(Form("fHControlMC%srapidity", name), fControlMCNp[0], 0); - FillHistSet(Form("fHControlMC%sptMin", name), fControlMCNp[1], 1); - FillHistSet(Form("fHControlMC%spt", name), fControlMCNp[2], 2); - FillHistSet(Form("fHControlMC%seta", name), fControlMCNp[3], 3); - FillHistSet(Form("fHControlMC%setapt", name), fControlMCNp[4], 4); + // -- Add Particle / Anti-Particle Distributions + for (Int_t idxPart = 0; idxPart < 2; ++idxPart) { + list->Add(new TH2D(Form("h%s%sX", name, fHelper->GetParticleName(idxPart).Data()), + Form("N_{%s} : %s;Centrality;N_{%s}", fHelper->GetParticleTitleLatex(idxPart).Data(), sTitle.Data(), fHelper->GetParticleTitleLatex(idxPart).Data()), + nBinsCent, centBinRange[0], centBinRange[1], 2601, -0.5, 2600.49)); + } // for (Int_t idxPart = 0; idxPart < 2; ++idxPart) { + + // -- Add NetParticle Distributions + list->Add(new TH2D(Form("h%sNet%sX", name, fHelper->GetParticleName(1).Data()), + Form("%s : %s;Centrality;%s", sNetTitle.Data(), sTitle.Data(), sNetTitle.Data()), + nBinsCent, centBinRange[0], centBinRange[1], 601, -300.5, 300.49)); + + // -- Add NetParticle vs SumParticle + list->Add(new TH2D(Form("h%sNet%sOverSumX", name, fHelper->GetParticleName(1).Data()), + Form("(%s)/(%s) : %s;Centrality;(%s)/(%s)", sNetTitle.Data(), sSumTitle.Data(), sTitle.Data(), sNetTitle.Data(), sSumTitle.Data()), + nBinsCent, centBinRange[0], centBinRange[1], 41, -2.5, 2.49)); + + // ----------------------------------------------------------------------------------------------- + // -- Add TProfiles for + // ----------------------------------------------------------------------------------------------- + for (Int_t idx = 1; idx <= fOrder; ++idx) { + list->Add(new TProfile(Form("p%sNet%s%dM", name, fHelper->GetParticleName(1).Data(), idx), + Form("(%s)^{%d} : %s;Centrality;(%s)^{%d}", sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx), + nBinsCent, centBinRange[0], centBinRange[1])); + } - return 0; -} + // ----------------------------------------------------------------------------------------------- + // -- Add TProfiles for for every SubSample + // ----------------------------------------------------------------------------------------------- + for (Int_t idxSub = 0; idxSub < fHelper->GetNSubSamples(); ++ idxSub) { + for (Int_t idx = 1; idx <= fOrder; ++idx) { + list->Add(new TProfile(Form("p%sNet%s%dM_%02d", name, fHelper->GetParticleName(1).Data(), idx, idxSub), + Form("(%s)^{%d} : %s;Centrality;(%s)^{%d}", sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx), + nBinsCent, centBinRange[0], centBinRange[1])); + } + } + + // ----------------------------------------------------------------------------------------------- + // -- Add TProfiles for + // ----------------------------------------------------------------------------------------------- + list->Add(new TList); + TList *fikList = static_cast(list->Last()); + fikList->SetName(Form("f%sFik",name)); + fikList->SetOwner(kTRUE); + + for (Int_t ii = 0; ii <= fOrder; ++ii) { + for (Int_t kk = 0; kk <= fOrder; ++kk) { + fikList->Add(new TProfile(Form("p%sNet%sF%02d%02d", name, fHelper->GetParticleName(1).Data(), ii, kk), + Form("f_{%02d%02d} : %s;Centrality;f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk), + nBinsCent, centBinRange[0], centBinRange[1])); + } + } -//////////////////////////////////////////////////////////////////////////////////////////// + // -- Add counter for number of Non-zero entries + for (Int_t idxCent = 0; idxCent < nBinsCent; ++idxCent) + fikList->Add(new TH2D(Form("p%sNet%sFCounts_%02d", name, fHelper->GetParticleName(1).Data(), idxCent), + Form("f_ik counts : %s Cent %d;i;k;", sTitle.Data(), idxCent), + fOrder+1, -0.5, Double_t(fOrder)+0.49, fOrder+1, -0.5, Double_t(fOrder)+0.49)); + + // ----------------------------------------------------------------------------------------------- + // -- Add TProfiles for for every SubSample + // ----------------------------------------------------------------------------------------------- + for (Int_t idxSub = 0; idxSub < fHelper->GetNSubSamples(); ++ idxSub) { + list->Add(new TList); + TList *fikListSub = static_cast(list->Last()); + fikListSub->SetName(Form("f%sFik_%02d",name, idxSub)); + fikListSub->SetOwner(kTRUE); + + for (Int_t ii = 0; ii <= fOrder; ++ii) { + for (Int_t kk = 0; kk <= fOrder; ++kk) { + fikListSub->Add(new TProfile(Form("p%sNet%sF%02d%02d_%02d", name, fHelper->GetParticleName(1).Data(), ii, kk, idxSub), + Form("f_{%02d%02d} : %s;Centrality;f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk), + nBinsCent, centBinRange[0], centBinRange[1])); + } + } + // -- Add counter for number of Non-zero entries + for (Int_t idxCent = 0; idxCent < nBinsCent; ++idxCent) + fikListSub->Add(new TH2D(Form("p%sNet%sFCounts_%02d_%02d", name, fHelper->GetParticleName(1).Data(), idxCent, idxSub), + Form("f_ik counts : %s Cent %d;i;k;", sTitle.Data(), idxCent), + fOrder+1, -0.5, Double_t(fOrder)+0.49, fOrder+1, -0.5, Double_t(fOrder)+0.49)); + } + + return; +} //________________________________________________________________________ -void AliAnalysisNetParticleDistribution::AddHistSet(const Char_t *name, const Char_t *title) { +void AliAnalysisNetParticleDistribution::AddHistSetCentPt(const Char_t *name, const Char_t *title) { // -- Add histogram sets for particle and anti-particle + // dependence : centrality and pt + + TString sName(name); + TString sTitle(title); - Int_t pid = fHelper->GetParticleSpecies(); - if( pid < 0 || pid > AliPID::kSPECIES){ - AliError("Particle ID not in AliPID::kSPECIES"); - return; - } - + // -- Add List fOutList->Add(new TList); TList *list = static_cast(fOutList->Last()); - list->SetName(name) ; + list->SetName(Form("f%s", name)); list->SetOwner(kTRUE); - TString sName(name); + // -- Get Bin Ranges + Int_t nBinsPt = AliAnalysisNetParticleHelper::fgkfHistNBinsPt; + Double_t ptBinRange[] = {-0.5, AliAnalysisNetParticleHelper::fgkfHistNBinsPt - 0.5}; - TString sPtTitle(""); - if (!sName.Contains("fHMC") || !sName.Contains("fHControlMC")) - sPtTitle += Form("#it{p}_{T} [%.1f,%.1f]", fPtRange[0], fPtRange[1]); - + Int_t nBinsCent = AliAnalysisNetParticleHelper::fgkfHistNBinsCent; + Double_t centBinRange[] = {AliAnalysisNetParticleHelper::fgkfHistRangeCent[0], AliAnalysisNetParticleHelper::fgkfHistRangeCent[1]}; + + // -- Create Titles + TString sNetTitle(Form("N_{%s} - N_{%s}", fHelper->GetParticleTitleLatex(1).Data(), fHelper->GetParticleTitleLatex(0).Data())); + TString sSumTitle(Form("N_{%s} + N_{%s}", fHelper->GetParticleTitleLatex(1).Data(), fHelper->GetParticleTitleLatex(0).Data())); + + // -- Add Particle / Anti-Particle Distributions for (Int_t idxPart = 0; idxPart < 2; ++idxPart) { - list->Add(new TH2F(Form("%s%s", name, fHelper->GetParticleName(idxPart).Data()), - Form("%s %s Dist %s;Centrality;N_{%s}", title, fHelper->GetParticleTitle(idxPart).Data(), sPtTitle.Data(), fHelper->GetParticleTitleLatex(idxPart).Data()), - 24, -0.5, 11.5, 501, -0.5, 500.49)); + list->Add(new TH3D(Form("h%s%s", name, fHelper->GetParticleName(idxPart).Data()), + Form("N_{%s} : %s;Centrality;N_{%s}", fHelper->GetParticleTitleLatex(idxPart).Data(), sTitle.Data(), fHelper->GetParticleTitleLatex(idxPart).Data()), + nBinsCent, centBinRange[0], centBinRange[1], nBinsPt, ptBinRange[0], ptBinRange[1], 2601, -0.5, 2600.49)); } // for (Int_t idxPart = 0; idxPart < 2; ++idxPart) { - list->Add(new TH2F(Form("%sNet%s", name, fHelper->GetParticleName(1).Data()), - Form("%s Net %s Dist %s;Centrality;N_{%s} - N_{%s}", title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()), 24, -0.5, 11.5, 201, -100.5, 100.49)); - - list->Add(new TProfile(Form("%sNet%sM", name, fHelper->GetParticleName(1).Data()), - Form("%s Net %s Dist %s;Centrality;N_{%s} - N_{%s}", title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()), 24, -0.5, 11.5)); - list->Add(new TProfile(Form("%sNet%s2M", name, fHelper->GetParticleName(1).Data()), - Form("%s (Net %s)^{2} Dist %s;Centrality;(N_{%s} - N_{%s})^{2}", title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()), 24, -0.5, 11.5)); - list->Add(new TProfile(Form("%sNet%s3M", name, fHelper->GetParticleName(1).Data()), - Form("%s (Net %s)^{3} Dist %s;Centrality;(N_{%s} - N_{%s})^{3}", title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()), 24, -0.5, 11.5)); - list->Add(new TProfile(Form("%sNet%s4M", name, fHelper->GetParticleName(1).Data()), - Form("%s (Net %s)^{4} Dist %s;Centrality;(N_{%s} - N_{%s})^{4}", title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()), 24, -0.5, 11.5)); - list->Add(new TProfile(Form("%sNet%s5M", name, fHelper->GetParticleName(1).Data()), - Form("%s (Net %s)^{5} Dist %s;Centrality;(N_{%s} - N_{%s})^{5}", title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()), 24, -0.5, 11.5)); - list->Add(new TProfile(Form("%sNet%s6M", name, fHelper->GetParticleName(1).Data()), - Form("%s (Net %s)^{6} Dist %s;Centrality;(N_{%s} - N_{%s})^{6}", title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()), 24, -0.5, 11.5)); - - list->Add(new TH2F(Form("%sNet%sOverSum", name, fHelper->GetParticleName(1).Data()), - Form("%s (Net %s)/ Sum Dist %s;Centrality;(N_{%s} - N_{%s})/(N_{%s} + N_{%s})", title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()), - 24, -0.5, 11.5, 801, -50.5, 50.49)); - - if (sName.Contains("fHControlMC")) { - for (Int_t idxPart = 0; idxPart < 2; ++idxPart) { - list->Add(new TH2F(Form("%s%sOver%s", name, fHelper->GetParticleName(idxPart).Data(), fHelper->GetControlParticleName(idxPart).Data()), - Form("%s %s / %s Dist %s;Centrality;N_{%s}/N_{Tracks}", - title, fHelper->GetParticleTitle(idxPart).Data(), fHelper->GetControlParticleTitle(idxPart).Data(), sPtTitle.Data(), fHelper->GetParticleTitleLatex(idxPart).Data()), - 24, -0.5, 11.5, 101, 0., 1.)); - } // for (Int_t idxPart = 0; idxPart < 2; ++idxPart) { + // -- Add NetParticle Distributions + list->Add(new TH3D(Form("h%sNet%s", name, fHelper->GetParticleName(1).Data()), + Form("%s : %s;Centrality;%s", sNetTitle.Data(), sTitle.Data(), sNetTitle.Data()), + nBinsCent, centBinRange[0], centBinRange[1], nBinsPt, ptBinRange[0], ptBinRange[1], 601, -300.5, 300.49)); + + // -- Add NetParticle vs SumParticle + list->Add(new TH3D(Form("h%sNet%sOverSum", name, fHelper->GetParticleName(1).Data()), + Form("(%s)/(%s) : %s;Centrality;(%s)/(%s)", sNetTitle.Data(), sSumTitle.Data(), sTitle.Data(), sNetTitle.Data(), sSumTitle.Data()), + nBinsCent, centBinRange[0], centBinRange[1], nBinsPt, ptBinRange[0], ptBinRange[1], 41, -2.5, 2.49)); + + // ----------------------------------------------------------------------------------------------- + // -- Add TProfiles for + // ----------------------------------------------------------------------------------------------- + for (Int_t idx = 1; idx <= fOrder; ++idx) { + list->Add(new TProfile2D(Form("p%sNet%s%dM", name, fHelper->GetParticleName(1).Data(), idx), + Form("(%s)^{%d} : %s;Centrality;(%s)^{%d}", sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx), + nBinsCent, centBinRange[0], centBinRange[1], nBinsPt, ptBinRange[0], ptBinRange[1])); + } + + // ----------------------------------------------------------------------------------------------- + // -- Add TProfiles for for every SubSample + // ----------------------------------------------------------------------------------------------- + for (Int_t idxSub = 0; idxSub < fHelper->GetNSubSamples(); ++ idxSub) { + for (Int_t idx = 1; idx <= fOrder; ++idx) { + list->Add(new TProfile2D(Form("p%sNet%s%dM_%02d", name, fHelper->GetParticleName(1).Data(), idx, idxSub), + Form("(%s)^{%d} : %s;Centrality;(%s)^{%d}", sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx), + nBinsCent, centBinRange[0], centBinRange[1], nBinsPt, ptBinRange[0], ptBinRange[1])); + } } + + // ----------------------------------------------------------------------------------------------- + // -- Add TProfiles for + // ----------------------------------------------------------------------------------------------- + list->Add(new TList); + TList *fikListPt = static_cast(list->Last()); + fikListPt->SetName(Form("f%sPtFik",name)); + fikListPt->SetOwner(kTRUE); + + for (Int_t ii = 0; ii <= fOrder; ++ii) { + for (Int_t kk = 0; kk <= fOrder; ++kk) { + fikListPt->Add(new TProfile2D(Form("p%sNet%sF%02d%02d", name, fHelper->GetParticleName(1).Data(), ii, kk), + Form("f_{%02d%02d} : %s;Centrality;#it{p}_{T} (GeV/#it{c});f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk), + nBinsCent, centBinRange[0], centBinRange[1], nBinsPt, ptBinRange[0], ptBinRange[1])); + } + } + + // -- Add counter for number of Non-zero entries + for (Int_t idxCent = 0; idxCent < nBinsCent; ++idxCent) + fikListPt->Add(new TH3D(Form("p%sNet%sFCounts_%02d", name, fHelper->GetParticleName(1).Data(), idxCent), + Form("f_ik counts : %s Cent %d;i;k;#it{p}_{T} Bins", sTitle.Data(), idxCent), + fOrder+1, -0.5, Double_t(fOrder)+0.49, fOrder+1, -0.5, Double_t(fOrder)+0.49, nBinsPt+1, -0.5, Double_t(nBinsPt)+0.49)); - // if (sName.Contains("fHDist")) { - // list->Add(new TH3F(Form("%sNet%s_RecVsMC", name, fHelper->GetParticleName(1).Data()), - // Form("%s Net %s Dist - Rec Vs MC - %s;Centrality;(N_{%s} - N_{%s})_{rec};(N_{%s} - N_{%s})_{MC}", - // title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()), - // 24, -0.5, 11.5, 201, -100.5, 100.49, 201, -100.5, 100.49)); - // } - + // ----------------------------------------------------------------------------------------------- + // -- Add TProfiles for for every SubSample + // ----------------------------------------------------------------------------------------------- + for (Int_t idxSub = 0; idxSub < fHelper->GetNSubSamples(); ++ idxSub) { + list->Add(new TList); + TList *fikListPtSub = static_cast(list->Last()); + fikListPtSub->SetName(Form("f%sPtFik_%02d",name, idxSub)); + fikListPtSub->SetOwner(kTRUE); + + for (Int_t ii = 0; ii <= fOrder; ++ii) { + for (Int_t kk = 0; kk <= fOrder; ++kk) { + fikListPtSub->Add(new TProfile2D(Form("p%sNet%sF%02d%02d_%02d", name, fHelper->GetParticleName(1).Data(), ii, kk, idxSub), + Form("f_{%02d%02d} : %s;Centrality;#it{p}_{T} (GeV/#it{c});f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk), + nBinsCent, centBinRange[0], centBinRange[1], nBinsPt, ptBinRange[0], ptBinRange[1])); + } + } + + // -- Add counter for number of Non-zero entries + for (Int_t idxCent = 0; idxCent < nBinsCent; ++idxCent) + fikListPtSub->Add(new TH3D(Form("p%sNet%sFCounts_%02d_%02d", name, fHelper->GetParticleName(1).Data(), idxCent, idxSub), + Form("f_ik counts : %s Cent %d;i;k;#it{p}_{T} Bins", sTitle.Data(), idxCent), + fOrder+1, -0.5, Double_t(fOrder)+0.49, fOrder+1, -0.5, Double_t(fOrder)+0.49, nBinsPt+1, -0.5, Double_t(nBinsPt)+0.49)); + } + return; } //________________________________________________________________________ -void AliAnalysisNetParticleDistribution::FillHistSet(const Char_t *name, Float_t *np, Int_t controlIdx) { - // -- Add histogram sets for particle and anti-particle +void AliAnalysisNetParticleDistribution::FillHistSetCent(const Char_t *name, Int_t idx, Bool_t isMC) { + // -- Fill histogram sets for particle and anti-particle + // dependence : centrality + + // -- Get List + TList *list = static_cast(fOutList->FindObject(Form("f%s",name))); + + // -- Get Centrality Bin + Float_t centralityBin = fHelper->GetCentralityBin(); + + // -- Select MC or Data + Int_t **np = (isMC) ? fMCNp : fNp; - Int_t pid = fHelper->GetParticleSpecies(); - if( pid < 0 || pid > AliPID::kSPECIES){ - AliError("Particle ID not in AliPID::kSPECIES"); - return; + // ----------------------------------------------------------------------------------------------- + + Int_t sumNp = np[idx][1]+np[idx][0]; // p + pbar + Int_t deltaNp = np[idx][1]-np[idx][0]; // p - pbar + + // -- Fill Particle / Anti-Particle Distributions + (static_cast(list->FindObject(Form("h%s%s", name, fHelper->GetParticleName(0).Data()))))->Fill(centralityBin, np[idx][0]); + (static_cast(list->FindObject(Form("h%s%s", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, np[idx][1]); + + // -- Fill NetParticle Distributions + (static_cast(list->FindObject(Form("h%sNet%s", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp); + + // -- Fill NetParticle vs SumParticle + Double_t deltaNpOverSumNp = (sumNp == 0.) ? 0. : deltaNp/Double_t(sumNp); + (static_cast(list->FindObject(Form("h%sNet%sOverSum", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNpOverSumNp); + + // ----------------------------------------------------------------------------------------------- + + Double_t CENT[7] = {1.157990, 1.153218, 1.147382, 1.142218, 1.139568, 1.138152, 1.134517}; + + Double_t sumNpX = np[idx][1]+(np[idx][0]*CENT[Int_t(centralityBin)]); + Double_t deltaNpX = np[idx][1]-(np[idx][0]*CENT[Int_t(centralityBin)]); + + // -- Fill Particle / Anti-Particle Distributions + (static_cast(list->FindObject(Form("h%s%sX", name, fHelper->GetParticleName(0).Data()))))->Fill(centralityBin, np[idx][0]*CENT[Int_t(centralityBin)]); + (static_cast(list->FindObject(Form("h%s%sX", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, np[idx][1]); + + // -- Fill NetParticle Distributions + (static_cast(list->FindObject(Form("h%sNet%sX", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNpX); + + // -- Fill NetParticle vs SumParticle + Double_t deltaNpXOverSumNpX = (sumNpX == 0.) ? 0. : deltaNpX/sumNpX; + (static_cast(list->FindObject(Form("h%sNet%sOverSumX", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNpXOverSumNpX); + + // ----------------------------------------------------------------------------------------------- + + // -- Fill TProfile for + Double_t delta = 1.; + for (Int_t idxOrder = 1; idxOrder <= fOrder; ++idxOrder) { + delta *= deltaNp; + (static_cast(list->FindObject(Form("p%sNet%s%dM", name, fHelper->GetParticleName(1).Data(), idxOrder))))->Fill(centralityBin, delta); + (static_cast(list->FindObject(Form("p%sNet%s%dM_%02d", name, fHelper->GetParticleName(1).Data(), idxOrder, fHelper->GetSubSampleIdx()))))->Fill(centralityBin, delta); } - TList *list = static_cast(fOutList->FindObject(name)); + // -- Generate reduced factorials - explictly removing the factorials + // - Reset all to 1 + for (Int_t idxOrder = 0; idxOrder <= fOrder; ++ idxOrder) { + fRedFactp[idxOrder][0] = 1.; + fRedFactp[idxOrder][1] = 1.; + } + + // - start at idx 1, as idx 0 = 1 done by reset + for (Int_t idxOrder = 1; idxOrder <= fOrder; ++ idxOrder) { + fRedFactp[idxOrder][0] = fRedFactp[idxOrder-1][0] * Double_t(np[idx][0]-(idxOrder-1)); + fRedFactp[idxOrder][1] = fRedFactp[idxOrder-1][1] * Double_t(np[idx][1]-(idxOrder-1)); + } + // -- Fill TProfiles for + TList *fikList = static_cast(list->FindObject(Form("f%sFik",name))); + TList *fikListSub = static_cast(list->FindObject(Form("f%sFik_%02d",name, fHelper->GetSubSampleIdx()))); + TH2D *hCntik = static_cast(fikList->FindObject(Form("p%sNet%sFCounts_%02d", name, fHelper->GetParticleName(1).Data(), Int_t(centralityBin)))); + TH2D *hCntikSub = static_cast(fikListSub->FindObject(Form("p%sNet%sFCounts_%02d_%02d", name, fHelper->GetParticleName(1).Data(), Int_t(centralityBin), fHelper->GetSubSampleIdx()))); + + for (Int_t ii = 0; ii <= fOrder; ++ii) { // ii -> p -> n1 + for (Int_t kk = 0; kk <= fOrder; ++kk) { // kk -> pbar -> n2 + // -- use the reduced factorials only + Double_t fik = fRedFactp[ii][1] * fRedFactp[kk][0]; // n1 *n2 -> p * pbar + (static_cast(fikList->FindObject(Form("p%sNet%sF%02d%02d", name, fHelper->GetParticleName(1).Data(), ii, kk))))->Fill(centralityBin, fik); + (static_cast(fikListSub->FindObject(Form("p%sNet%sF%02d%02d_%02d", + name, fHelper->GetParticleName(1).Data(), ii, kk, fHelper->GetSubSampleIdx()))))->Fill(centralityBin, fik); + + if (fik != 0.) { + hCntik->Fill(ii, kk); + hCntikSub->Fill(ii, kk); + } + } + } + + return; +} + +//________________________________________________________________________ +void AliAnalysisNetParticleDistribution::FillHistSetCentPt(const Char_t *name, Int_t idx, Bool_t isMC) { + // -- Add histogram sets for particle and anti-particle + // dependence : centrality and pt + + // -- Get List + TList *list = static_cast(fOutList->FindObject(Form("f%s",name))); + + // -- Get Centrality Bin Float_t centralityBin = fHelper->GetCentralityBin(); - (static_cast(list->FindObject(Form("%s%s", name, fHelper->GetParticleName(0).Data()))))->Fill(centralityBin, np[0]); - (static_cast(list->FindObject(Form("%s%s", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, np[1]); + // -- Select MC or Data + Int_t ***npPt = (isMC) ? fMCNpPt : fNpPt; + + // ----------------------------------------------------------------------------------------------- - Float_t sumNp = np[1]+np[0]; - Float_t deltaNp = np[1]-np[0]; - Float_t deltaNp2 = deltaNp * deltaNp; - Float_t deltaNp3 = deltaNp2 * deltaNp; + // -- Loop over the pt bins + for (Int_t idxPt = 0; idxPt < AliAnalysisNetParticleHelper::fgkfHistNBinsPt; ++idxPt) { + + Int_t deltaNp = npPt[idx][1][idxPt]-npPt[idx][0][idxPt]; // p - pbar + Int_t sumNp = npPt[idx][1][idxPt]+npPt[idx][0][idxPt]; // p + pbar - (static_cast(list->FindObject(Form("%sNet%s", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp); + // -- Fill Particle / Anti-Particle Distributions + (static_cast(list->FindObject(Form("h%s%s", name, fHelper->GetParticleName(0).Data()))))->Fill(centralityBin, idxPt, npPt[idx][0][idxPt]); + (static_cast(list->FindObject(Form("h%s%s", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, idxPt, npPt[idx][1][idxPt]); + + // -- Fill NetParticle Distributions + (static_cast(list->FindObject(Form("h%sNet%s", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, idxPt, deltaNp); + + // -- Fill NetParticle vs SumParticle + Double_t deltaNpOverSumNp = (sumNp == 0.) ? 0. : deltaNp/Double_t(sumNp); + (static_cast(list->FindObject(Form("h%sNet%sOverSum", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, idxPt, deltaNpOverSumNp); - (static_cast(list->FindObject(Form("%sNet%sM", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp); - (static_cast(list->FindObject(Form("%sNet%s2M", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp2); - (static_cast(list->FindObject(Form("%sNet%s3M", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp2*deltaNp); - (static_cast(list->FindObject(Form("%sNet%s4M", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp2*deltaNp2); - (static_cast(list->FindObject(Form("%sNet%s5M", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp3*deltaNp2); - (static_cast(list->FindObject(Form("%sNet%s6M", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp3*deltaNp3); + // ----------------------------------------------------------------------------------------------- - (static_cast(list->FindObject(Form("%sNet%sOverSum", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp/sumNp); + // -- Fill TProfile for + Double_t delta = 1.; + for (Int_t idxOrder = 1; idxOrder <= fOrder; ++idxOrder) { + delta *= deltaNp; + (static_cast(list->FindObject(Form("p%sNet%s%dM", name, fHelper->GetParticleName(1).Data(), idxOrder))))->Fill(centralityBin, idxPt, delta); + (static_cast(list->FindObject(Form("p%sNet%s%dM_%02d", + name, fHelper->GetParticleName(1).Data(), idxOrder, fHelper->GetSubSampleIdx()))))->Fill(centralityBin, idxPt, delta); + } + + // -- Generate reduced factorials - explictly removing the factorials + // - Reset all to 1 + for (Int_t idxOrder = 0; idxOrder <= fOrder; ++ idxOrder) { + fRedFactp[idxOrder][0] = 1.; + fRedFactp[idxOrder][1] = 1.; + } - TString sName(name); - if (sName.Contains("fHControlMC") && controlIdx >= 0) { - (static_cast(list->FindObject(Form("%s%sOverlambdabar", name, fHelper->GetParticleName(0).Data()))))->Fill(centralityBin, np[0]/fControlMCNp[controlIdx][0]); - (static_cast(list->FindObject(Form("%s%sOverlambda", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, np[1]/fControlMCNp[controlIdx][1]); - } + // - start at idx 1, as idx 0 = 1 done by reset + for (Int_t idxOrder = 1; idxOrder <= fOrder; ++ idxOrder) { + fRedFactp[idxOrder][0] = fRedFactp[idxOrder-1][0] * Double_t(npPt[idx][0][idxPt]-(idxOrder-1)); + fRedFactp[idxOrder][1] = fRedFactp[idxOrder-1][1] * Double_t(npPt[idx][1][idxPt]-(idxOrder-1)); + } - // if (sName.Contains("fHDist")) { - // Float_t deltaMCNp = fMCNp[controlIdx][1]-fMCNp[controldx][0]; - // (static_cast(list->FindObject(Form("%sNetp_RecVsMC", name))))->Fill(centralityBin, deltaNp, deltaMCNp); - // } + // -- Fill TProfiles for + TList *fikListPt = static_cast(list->FindObject(Form("f%sPtFik",name))); + TList *fikListPtSub = static_cast(list->FindObject(Form("f%sPtFik_%02d",name, fHelper->GetSubSampleIdx()))); + TH3D *hCntikPt = static_cast(fikListPt->FindObject(Form("p%sNet%sFCounts_%02d", name, fHelper->GetParticleName(1).Data(), Int_t(centralityBin)))); + TH3D *hCntikPtSub = static_cast(fikListPtSub->FindObject(Form("p%sNet%sFCounts_%02d_%02d", name, fHelper->GetParticleName(1).Data(), + Int_t(centralityBin), fHelper->GetSubSampleIdx()))); + for (Int_t ii = 0; ii <= fOrder; ++ii) { // ii -> p -> n1 + for (Int_t kk = 0; kk <= fOrder; ++kk) { // kk -> pbar -> n2 + Double_t fik = fRedFactp[ii][1] * fRedFactp[kk][0]; // n1 *n2 -> p * pbar + (static_cast(fikListPt->FindObject(Form("p%sNet%sF%02d%02d", name, fHelper->GetParticleName(1).Data(), ii, kk))))->Fill(centralityBin, idxPt, fik); + (static_cast(fikListPtSub->FindObject(Form("p%sNet%sF%02d%02d_%02d", + name, fHelper->GetParticleName(1).Data(), ii, kk, fHelper->GetSubSampleIdx()))))->Fill(centralityBin, idxPt, fik); + + if (fik != 0.) { + hCntikPt->Fill(ii, kk, idxPt); + hCntikPtSub->Fill(ii, kk, idxPt); + } + } + } + + } // for (Int_t idxPt = 0; idxPt < AliAnalysisNetParticleHelper::fgkfHistNBinsPt; ++idxPt) { return; }