diff --git a/PWGCF/MultiparticleCorrelations/Tasks/multiharmonicCorrelations.cxx b/PWGCF/MultiparticleCorrelations/Tasks/multiharmonicCorrelations.cxx index 5ea45627e26..d19e1b7233c 100644 --- a/PWGCF/MultiparticleCorrelations/Tasks/multiharmonicCorrelations.cxx +++ b/PWGCF/MultiparticleCorrelations/Tasks/multiharmonicCorrelations.cxx @@ -25,6 +25,7 @@ #include #include +#include #include #include #include @@ -88,20 +89,25 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to Configurable cfDryRun{"cfDryRun", false, "book all histos and run without filling and calculating anything"}; // example for built-in type (float, string, etc.) Configurable> cf_pt_bins{"cf_pt_bins", {1000, 0., 100.}, "nPtBins, ptMin, ptMax"}; // example for an array Configurable> cf_phi_bins{"cf_phi_bins", {100, 0., o2::constants::math::TwoPI}, "nPhiBins, phiMin, phiMax"}; - Configurable> cf_centr_bins{"cf_centr_bins", {1000, 0., 100.}, "nCentrBins, centrMin, centrMax"}; + Configurable> cf_centr_bins{"cf_centr_bins", {100, 0., 100.}, "nCentrBins, centrMin, centrMax"}; Configurable> cf_x_bins{"cf_x_bins", {1000, -100., 100.}, "nXBins, xMin, xMax"}; Configurable> cf_y_bins{"cf_y_bins", {1000, -100., 100.}, "nYBins, yMin, yMax"}; Configurable> cf_z_bins{"cf_z_bins", {1000, -100., 100.}, "nZBins, zMin, zMax"}; Configurable> cf_mult_bins{"cf_mult_bins", {50, 0, 3e3}, "nMultBins, multMin, multMax"}; + Configurable> cf_tpcncls_bins{"cf_tpcncls_bins", {100, 0., 1000.}, "ntpcnclsBins, tpnclsMin, tpcnclsMax"}; + Configurable> cf_dcaxy_bins{"cf_dcaxy_bins", {1000, -20., 20.}, "ndcaxyBins, dcaxyMin, dcaxyMax"}; + Configurable> cf_dcaz_bins{"cf_dcaz_bins", {1000, -10., 10.}, "ndcazBins, dcazMin, dcazMax"}; + Configurable> cf_ncontr_bins{"cf_ncontr_bins", {100, 0., 1000}, "nNContrBins, NContrMin, NContrMax"}; Configurable cfCent{"cfCent", 1, "centrality estimator"}; - Configurable cfMult{"cfMult", 1, "multiplicity"}; + Configurable cfMult{"cfMult", "TPC", "multiplicity"}; Configurable cfQA{"cfQA", true, "quality assurance"}; Configurable> cfVertexZ{"cfVertexZ", {-10, 10.}, "vertex z position range: {min, max}[cm], with convention: min <= Vz < max"}; Configurable> cfPt{"cfPt", {0.2, 5.0}, "transverse momentum range"}; + Configurable> cfEta{"cfEta", {-0.8, 0.8}, "eta range"}; - Configurable cfFileWithWeights{"cfFileWithWeights", "~/O2/weights.root", "path to external ROOT file which holds all particle weights"}; + Configurable cfFileWithWeights{"cfFileWithWeights", "/alice-ccdb.cern.ch/Users/p/pengchon/test04", "path to external ROOT file which holds all particle weights"}; // *) Define and initialize all data members to be called in the main process* functions: // **) Task configuration: @@ -115,6 +121,10 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to TList* fParticleHistogramsList = NULL; //!Getenv("USER"), "", ""); // do not forget to add #include to the preamble of your analysis task @@ -302,7 +329,7 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to float vertexZmin = static_cast(vertexZ[0]); float vertexZmax = static_cast(vertexZ[1]); float posZ = collision.posZ(); - if (posZ < vertexZmin || posZ > vertexZmax) + if (posZ < vertexZmin || posZ > vertexZmax || (!collision.sel8())) return false; return true; } @@ -314,11 +341,32 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to float ptcutmin = static_cast(Pt[0]); float ptcutmax = static_cast(Pt[1]); float pt = track.pt(); - if (pt < ptcutmin || pt > ptcutmax) + vector Eta = cfEta.value; + float etacutmin = static_cast(Eta[0]); + float etacutmax = static_cast(Eta[1]); + float eta = track.eta(); + if (pt < ptcutmin || pt > ptcutmax || eta < etacutmin || eta > etacutmax) return false; return true; } + TComplex Q(Int_t n) + { + // Using the fact that Q{-n,p} = Q{n,p}^*. + if (n >= 0) { + return cor.Qvector[n]; + } + return TComplex::Conjugate(cor.Qvector[-n]); + } + + TComplex Four(Int_t n1, Int_t n2, Int_t n3, Int_t n4) + { // Generic four-particle correlation . + TComplex four = Q(n1) * Q(n2) * Q(n3) * Q(n4) - Q(n1 + n2) * Q(n3) * Q(n4) - Q(n2) * Q(n1 + n3) * Q(n4) - Q(n1) * Q(n2 + n3) * Q(n4) + 2. * Q(n1 + n2 + n3) * Q(n4) - Q(n2) * Q(n3) * Q(n1 + n4) + Q(n2 + n3) * Q(n1 + n4) - Q(n1) * Q(n3) * Q(n2 + n4) + Q(n1 + n3) * Q(n2 + n4) + 2. * Q(n3) * Q(n1 + n2 + n4) - Q(n1) * Q(n2) * Q(n3 + n4) + Q(n1 + n2) * Q(n3 + n4) + 2. * Q(n2) * Q(n1 + n3 + n4) + 2. * Q(n1) * Q(n2 + n3 + n4) - 6. * Q(n1 + n2 + n3 + n4); + + return four; + + } // TComplex Four(Int_t n1, Int_t n2, Int_t n3, Int_t n4) + template void Steer(T1 const& collision, T2 const& tracks) { @@ -327,37 +375,40 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to return; } // Print current run number: - LOGF(info, "Run number: %d", collision.bc().runNumber()); - // Print centrality estimated with "FT0M" estimator: - LOGF(info, "Centrality: %f", collision.centFT0M()); - // Print vertex X position: - LOGF(info, "Vertex X position: %f", collision.posX()); + // LOGF(info, "Run number: %d", collision.bc().runNumber()); + + float zrec = 0., zsim = 0., centr = 0, M = 0.; + + float v22, v32, v42; + float four32, four42; - float zrec = 0., zsim = 0.; if constexpr (rs == eRec || rs == eRecAndSim) { event.fHistX[eRec]->Fill(collision.posX()); event.fHistY[eRec]->Fill(collision.posY()); event.fHistZ[eRec]->Fill(collision.posZ()); event.fEventHistograms[eVertexZ][eRec][0]->Fill(collision.posZ()); zrec = collision.posZ(); - float centr = 0; + event.fHistNContr->Fill(collision.numContrib()); if (cfCent == 1) centr = collision.centFT0C(); - if (cfCent == 2) + else if (cfCent == 2) centr = collision.centFT0M(); - if (cfCent == 3) + else if (cfCent == 3) centr = collision.centFT0A(); event.fHistCentr[eRec]->Fill(centr); - if (cfMult == 1) - event.fHistMult[eRec]->Fill(collision.multTPC()); - if (cfMult == 2) - event.fHistMult[eRec]->Fill(collision.multFV0M()); - if (cfMult == 3) - event.fHistMult[eRec]->Fill(collision.multFT0C()); - if (cfMult == 4) - event.fHistMult[eRec]->Fill(collision.multFT0M()); - if (cfMult == 5) - event.fHistMult[eRec]->Fill(collision.multNTracksPV()); + + std::string multType = "TPC"; + if (cfMult.value == "TPC") + M = collision.multTPC(); + else if (cfMult.value == "FV0M") + M = collision.multFV0M(); + else if (cfMult.value == "FT0C") + M = collision.multFT0C(); + else if (cfMult.value == "FT0M") + M = collision.multFT0M(); + else if (cfMult.value == "NTracksPV") + M = collision.multNTracksPV(); + event.fHistMult[eRec]->Fill(M); if constexpr (rs == eRecAndSim) { auto mccollision = collision.mcCollision(); @@ -369,8 +420,8 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to event.fEventHistograms[eVertexZ][eSim][0]->Fill(mccollision.posZ()); zsim = mccollision.posZ(); event.fHistCentr[eSim]->Fill(centrsim); - event.fQA->Fill(centrsim, centr); - // event.fHistMult[eSim]->Fill(tracks.size()); + qa.fQA->Fill(centrsim, centr); + centr = centrsim; } // *) Event cuts: @@ -382,12 +433,17 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to event.fEventHistograms[eVertexZ][eSim][1]->Fill(zsim); } + // before loop over particles + float eta = 0; + float phi = 0; + for (int ih = 0; ih < maxHarmonic; ih++) { + cor.Qvector[ih] = TComplex(0., 0.); + } + // Main loop over particles: auto track = tracks.iteratorAt(0); // set the type and scope from one instance for (int64_t i = 0; i < tracks.size(); i++) { - // Print track azimuthal angle: - // Fill reconstructed ...: float ptrec = 0., ptsim = 0.; if constexpr (rs == eRec || rs == eRecAndSim) { @@ -395,7 +451,13 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to pc.fHistPt[eRec]->Fill(track.pt()); event.fEventHistograms[ePt][eRec][0]->Fill(track.pt()); ptrec = track.pt(); + eta = track.eta(); + phi = track.phi(); pc.fHistPhi[eRec]->Fill(track.phi()); + pc.fHistCharge[eRec]->Fill(track.sign()); + pc.fHistTPCncls[eRec]->Fill(track.tpcNClsFindable()); + pc.fHistTracksdcaXY[eRec]->Fill(track.dcaXY()); + pc.fHistTracksdcaZ[eRec]->Fill(track.dcaZ()); // ... and corresponding MC truth simulated: // See https://github.com/AliceO2Group/O2Physics/blob/master/Tutorials/src/mcHistograms.cxx @@ -409,8 +471,14 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to pc.fHistPt[eSim]->Fill(mcparticle.pt()); event.fEventHistograms[ePt][eSim][0]->Fill(mcparticle.pt()); ptsim = mcparticle.pt(); + eta = mcparticle.eta(); + phi = mcparticle.phi(); pc.fHistPhi[eSim]->Fill(mcparticle.phi()); - } // end of if constexpr (rs == eRecAndSim) { + // pc.fHistCharge[eSim]->Fill(mcparticle.sign()); + // pc.fHistTPCncls[eSim]->Fill(mcparticle.tpcNClsFindable()); + // pc.fHistTracksdcaXY[eSim]->Fill(mcparticle.dcaXY()); + // pc.fHistTracksdcaZ[eSim]->Fill(mcparticle.dcaZ()); + } // end of if constexpr (rs == eRecAndSim) // *) Particle cuts: if (!ParticleCuts(track)) { // Main call for particle cuts. @@ -420,10 +488,27 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to if constexpr (rs == eRecAndSim) event.fEventHistograms[ePt][eSim][0]->Fill(ptsim); - } // if constexpr (rs == eRec || rs == eRecAndSim) { - } // end of for (int64_t i = 0; i < tracks.size(); i++) { + } // if constexpr (rs == eRec || rs == eRecAndSim) - } // end of template void Steer(T1 const& collision, T2 const& tracks) { + // analysis in the loop over particle + for (int ih = 0; ih < maxHarmonic; ih++) { + cor.Qvector[ih] += TComplex(TMath::Cos(ih * phi), TMath::Sin(ih * phi)); + } + } // end of for (int64_t i = 0; i < tracks.size(); i++) + // calculate correlations + four32 = Four(3, 2, -3, -2).Re() / Four(0, 0, 0, 0).Re(); + four42 = Four(4, 2, -4, -2).Re() / Four(0, 0, 0, 0).Re(); + v22 = (Q(2).Rho2() - M) / (M * (M - 1.)); + v32 = (Q(3).Rho2() - M) / (M * (M - 1.)); + v42 = (Q(4).Rho2() - M) / (M * (M - 1.)); + + cor.pv22_centr->Fill(centr, v22, M * (M - 1)); + cor.pv32_centr->Fill(centr, v32, M * (M - 1)); + cor.pv42_centr->Fill(centr, v42, M * (M - 1)); + cor.pfour32_centr->Fill(centr, four32, M * (M - 1)); + cor.pfour42_centr->Fill(centr, four42, M * (M - 1)); + + } // end of template void Steer(T1 const& collision, T2 const& tracks) // *) Initialize and book all objects: void init(InitContext&) @@ -454,6 +539,16 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to event.fEventHistogramsList->SetOwner(true); fBaseList->Add(event.fEventHistogramsList); + qa.fQAList = new TList(); + qa.fQAList->SetName("QAHistograms"); + qa.fQAList->SetOwner(true); + fBaseList->Add(qa.fQAList); + + cor.fCorrelationVariablesList = new TList(); + cor.fCorrelationVariablesList->SetName("CorrelationVariables"); + cor.fCorrelationVariablesList->SetOwner(true); + fBaseList->Add(cor.fCorrelationVariablesList); + // *) Book pt distribution with binning defined through configurables in the json file: vector l_pt_bins = cf_pt_bins.value; // define local array and initialize it from an array set in the configurables vector l_phi_bins = cf_phi_bins.value; @@ -462,6 +557,10 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to vector l_y_bins = cf_y_bins.value; vector l_z_bins = cf_z_bins.value; vector l_mult_bins = cf_mult_bins.value; + vector l_tpcncls_bins = cf_tpcncls_bins.value; + vector l_dcaxy_bins = cf_dcaxy_bins.value; + vector l_dcaz_bins = cf_dcaz_bins.value; + vector l_ncontr_bins = cf_ncontr_bins.value; int nBins = static_cast(l_pt_bins[0]); int nBinsphi = static_cast(l_phi_bins[0]); int nBinscentr = static_cast(l_centr_bins[0]); @@ -469,6 +568,11 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to int nBinsy = static_cast(l_y_bins[0]); int nBinsz = static_cast(l_z_bins[0]); int nBinsmult = static_cast(l_mult_bins[0]); + int nBinscharge = 2; + int nBinstpcncls = static_cast(l_tpcncls_bins[0]); + int nBinsdcaxy = static_cast(l_dcaxy_bins[0]); + int nBinsdcaz = static_cast(l_dcaz_bins[0]); + int nBinsncontr = static_cast(l_ncontr_bins[0]); float min = l_pt_bins[1]; float max = l_pt_bins[2]; @@ -484,20 +588,54 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to float maxz = l_z_bins[2]; float minmult = l_mult_bins[1]; float maxmult = l_mult_bins[2]; + float mincharge = -2.; + float maxcharge = 2.; + float mintpcncls = l_tpcncls_bins[1]; + float maxtpcncls = l_tpcncls_bins[2]; + float mindcaxy = l_dcaxy_bins[1]; + float maxdcaxy = l_dcaxy_bins[2]; + float mindcaz = l_dcaz_bins[1]; + float maxdcaz = l_dcaz_bins[2]; + float maxncontr = l_ncontr_bins[1]; + float minncontr = l_ncontr_bins[2]; pc.fHistPt[eRec] = new TH1F("fHistPt[eRec]", "pt distribution for reconstructed particles", nBins, min, max); - pc.fHistPt[eRec]->GetXaxis()->SetTitle("p_{T}"); - pc.fParticleHistogramsList->Add(pc.fHistPt[eRec]); pc.fHistPhi[eRec] = new TH1F("fHistPhi[eRec]", "phi distribution for reconstructed particles", nBinsphi, minphi, maxphi); + pc.fHistCharge[eRec] = new TH1F("fHistCharge[eRec]", "charge distribution for reconstructed particles", nBinscharge, mincharge, maxcharge); + pc.fHistTPCncls[eRec] = new TH1F("fHistTPCncls[eRec]", "tpcncls distribution for reconstructed particles", nBinstpcncls, mintpcncls, maxtpcncls); + pc.fHistTracksdcaXY[eRec] = new TH1F("fHistTracksdcaXY[eRec]", "dcaxy distribution for reconstructed particles", nBinsdcaxy, mindcaxy, maxdcaxy); + pc.fHistTracksdcaZ[eRec] = new TH1F("fHistTracksdcaZ[eRec]", "dcaz distribution for reconstructed particles", nBinsdcaz, mindcaz, maxdcaz); + pc.fHistPt[eRec]->GetXaxis()->SetTitle("p_{T}"); pc.fHistPhi[eRec]->GetXaxis()->SetTitle("phi"); + pc.fHistCharge[eRec]->GetXaxis()->SetTitle("charge"); + pc.fHistTPCncls[eRec]->GetXaxis()->SetTitle("TPCNClsFindable"); + pc.fHistTracksdcaXY[eRec]->GetXaxis()->SetTitle("DCA XY"); + pc.fHistTracksdcaZ[eRec]->GetXaxis()->SetTitle("DCA Z"); + pc.fParticleHistogramsList->Add(pc.fHistPt[eRec]); pc.fParticleHistogramsList->Add(pc.fHistPhi[eRec]); + pc.fParticleHistogramsList->Add(pc.fHistCharge[eRec]); + pc.fParticleHistogramsList->Add(pc.fHistTPCncls[eRec]); + pc.fParticleHistogramsList->Add(pc.fHistTracksdcaXY[eRec]); + pc.fParticleHistogramsList->Add(pc.fHistTracksdcaZ[eRec]); pc.fHistPt[eSim] = new TH1F("fHistPt[eSim]", "pt distribution for simulated particles", nBins, min, max); - pc.fHistPt[eSim]->GetXaxis()->SetTitle("p_{T}"); - pc.fParticleHistogramsList->Add(pc.fHistPt[eSim]); pc.fHistPhi[eSim] = new TH1F("fHistPhi[eSim]", "phi distribution for simulated particles", nBinsphi, minphi, maxphi); + pc.fHistCharge[eSim] = new TH1F("fHistCharge[eSim]", "charge distribution for simulated particles", nBinscharge, mincharge, maxcharge); + pc.fHistTPCncls[eSim] = new TH1F("fHistTPCncls[eSim]", "tpcncls distribution for simulated particles", nBinstpcncls, minphi, maxtpcncls); + pc.fHistTracksdcaXY[eSim] = new TH1F("fHistTracksdcaXY[eSim]", "dcaxy distribution for simulated particles", nBinsdcaxy, mindcaxy, maxdcaxy); + pc.fHistTracksdcaZ[eSim] = new TH1F("fHistTracksdcaZ[eSim]", "dcaz distribution for simulated particles", nBinsdcaz, mindcaz, maxdcaz); + pc.fHistPt[eSim]->GetXaxis()->SetTitle("p_{T}"); pc.fHistPhi[eSim]->GetXaxis()->SetTitle("phi"); + pc.fHistCharge[eSim]->GetXaxis()->SetTitle("charge"); + pc.fHistTPCncls[eSim]->GetXaxis()->SetTitle("TPCNClsFindable"); + pc.fHistTracksdcaXY[eSim]->GetXaxis()->SetTitle("DCA XY"); + pc.fHistTracksdcaZ[eSim]->GetXaxis()->SetTitle("DCA Z"); + pc.fParticleHistogramsList->Add(pc.fHistPt[eSim]); pc.fParticleHistogramsList->Add(pc.fHistPhi[eSim]); + pc.fParticleHistogramsList->Add(pc.fHistCharge[eSim]); + pc.fParticleHistogramsList->Add(pc.fHistTPCncls[eSim]); + pc.fParticleHistogramsList->Add(pc.fHistTracksdcaXY[eSim]); + pc.fParticleHistogramsList->Add(pc.fHistTracksdcaZ[eSim]); pc.histWeights = GetHistogramWithWeights(cfFileWithWeights.value.c_str(), "000123456"); pc.fParticleHistogramsList->Add(pc.histWeights); @@ -507,16 +645,19 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to event.fHistY[eRec] = new TH1F("fHistY[eRec]", "posY distribution for reconstructed particles", nBinsy, miny, maxy); event.fHistZ[eRec] = new TH1F("fHistZ[eRec]", "posZ distribution for reconstructed particles", nBinsz, minz, maxz); event.fHistMult[eRec] = new TH1I("fHistMult[eRec]", "mult distribution for reconstructed particles", nBinsmult, minmult, maxmult); + event.fHistNContr = new TH1I("fHistNContr", "NContr distribution", nBinsncontr, minncontr, maxncontr); event.fHistCentr[eRec]->GetXaxis()->SetTitle("centrality"); event.fHistX[eRec]->GetXaxis()->SetTitle("x"); event.fHistY[eRec]->GetXaxis()->SetTitle("y"); event.fHistZ[eRec]->GetXaxis()->SetTitle("z"); event.fHistMult[eRec]->GetXaxis()->SetTitle("multiplicity"); + event.fHistNContr->GetXaxis()->SetTitle("numContrib"); event.fEventHistogramsList->Add(event.fHistCentr[eRec]); event.fEventHistogramsList->Add(event.fHistX[eRec]); event.fEventHistogramsList->Add(event.fHistY[eRec]); event.fEventHistogramsList->Add(event.fHistZ[eRec]); event.fEventHistogramsList->Add(event.fHistMult[eRec]); + event.fEventHistogramsList->Add(event.fHistNContr); event.fHistCentr[eSim] = new TH1F("fHistCentr[eSim]", "centrality distribution for simulated particles", nBinscentr, mincentr, maxcentr); event.fHistX[eSim] = new TH1F("fHistX[eSim]", "posX distribution for simulated particles", nBinsx, minx, maxx); @@ -552,11 +693,29 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to } } - event.fQA = new TH2F("QA", "quality assurance", nBinscentr, mincentr, maxcentr, nBinscentr, mincentr, maxcentr); + qa.fQA = new TH2F("QA", "quality assurance", nBinscentr, mincentr, maxcentr, nBinscentr, mincentr, maxcentr); if (cfQA) { - event.fEventHistogramsList->Add(event.fQA); + qa.fQAList->Add(qa.fQA); } + // float quantiles[10] = {0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8}; + float quantiles[10] = {0, 5, 10, 20, 30, 40, 50, 60, 70, 80}; + cor.pv22_centr = new TProfile("pv22", "profile of v2^2", 9, quantiles); + cor.pv32_centr = new TProfile("pv32", "profile of v3^2", 9, quantiles); + cor.pv42_centr = new TProfile("pv42", "profile of v4^2", 9, quantiles); + cor.pfour32_centr = new TProfile("pfour32", "profile of v2^2*v3^2", 9, quantiles); + cor.pfour42_centr = new TProfile("pfour32", "profile of v2^2*v4^2", 9, quantiles); + cor.pv22_centr->GetXaxis()->SetTitle("v_{2}^{2}"); + cor.pv32_centr->GetXaxis()->SetTitle("v_{3}^{2}"); + cor.pv42_centr->GetXaxis()->SetTitle("v_{4}^{2}"); + cor.pfour32_centr->GetXaxis()->SetTitle("v_{2}^{2}v_{3}^{2}"); + cor.pfour42_centr->GetXaxis()->SetTitle("v_{2}^{2}v_{4}^{2}"); + cor.fCorrelationVariablesList->Add(cor.pv22_centr); + cor.fCorrelationVariablesList->Add(cor.pv32_centr); + cor.fCorrelationVariablesList->Add(cor.pv42_centr); + cor.fCorrelationVariablesList->Add(cor.pfour32_centr); + cor.fCorrelationVariablesList->Add(cor.pfour42_centr); + } // end of void init(InitContext&) { // A) Process only reconstructed data: