Skip to content

Commit 81f3728

Browse files
SwatiSaha-1997SwatiSaha-1997
committed
Calculate observables as a func. of Ngen in MCgen and rec
1 parent 6f9193a commit 81f3728

File tree

1 file changed

+80
-6
lines changed

1 file changed

+80
-6
lines changed

‎PWGCF/EbyEFluctuations/Tasks/MeanptFluctuations.cxx‎

Lines changed: 80 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@
2727
#include "Framework/runDataProcessing.h"
2828
#include <CCDB/BasicCCDBManager.h>
2929

30+
#include "TDatabasePDG.h"
3031
#include <TDirectory.h>
3132
#include <TF1.h>
3233
#include <TFile.h>
@@ -149,6 +150,8 @@ struct MeanptFluctuationsAnalysis {
149150
HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};
150151
std::vector<std::vector<std::shared_ptr<TProfile2D>>> subSampleMcGen;
151152
std::vector<std::vector<std::shared_ptr<TProfile2D>>> subSample;
153+
std::vector<std::vector<std::shared_ptr<TProfile>>> subSampleMcGenV2;
154+
std::vector<std::vector<std::shared_ptr<TProfile2D>>> subSampleV2;
152155
TRandom3* fRndm = new TRandom3(0);
153156

154157
// filtering collisions and tracks***********
@@ -246,6 +249,7 @@ struct MeanptFluctuationsAnalysis {
246249
histos.add("MCGenerated/hPtParticleVsTrack", "", kTH2F, {ptAxis, ptAxis});
247250
histos.add("MCGenerated/hEtaParticleVsTrack", "", kTH2F, {{100, -2.01, 2.01}, {100, -2.01, 2.01}});
248251
histos.add("MCGenerated/hPhiParticleVsTrack", "", kTH2F, {{100, 0., o2::constants::math::TwoPI}, {100, 0., o2::constants::math::TwoPI}});
252+
histos.add("MCGenerated/hNgenVsNrec", "Gen multiplicity vs. Rec multiplicity in MC", kTH2F, {multAxis, multAxis});
249253

250254
// Analysis Profiles for central val
251255
histos.add("MCGenerated/AnalysisProfiles/Prof_mean_t1", "", {HistType::kTProfile2D, {centAxis, multAxis}});
@@ -255,16 +259,39 @@ struct MeanptFluctuationsAnalysis {
255259
histos.add("MCGenerated/AnalysisProfiles/Hist2D_Nch_centrality", "", {HistType::kTH2D, {centAxis, multAxis}});
256260
histos.add("MCGenerated/AnalysisProfiles/Hist2D_meanpt_centrality", "", {HistType::kTH2D, {centAxis, meanpTAxis}});
257261

262+
histos.add("MCGenerated/AnalysisProfilesV2/Prof_mean_t1", "", {HistType::kTProfile, {multAxis}});
263+
histos.add("MCGenerated/AnalysisProfilesV2/Prof_var_t1", "", {HistType::kTProfile, {multAxis}});
264+
histos.add("MCGenerated/AnalysisProfilesV2/Prof_skew_t1", "", {HistType::kTProfile, {multAxis}});
265+
histos.add("MCGenerated/AnalysisProfilesV2/Prof_kurt_t1", "", {HistType::kTProfile, {multAxis}});
266+
histos.add("AnalysisProfilesV2/Prof_mean_t1", "", {HistType::kTProfile2D, {multAxis, multAxis}});
267+
histos.add("AnalysisProfilesV2/Prof_var_t1", "", {HistType::kTProfile2D, {multAxis, multAxis}});
268+
histos.add("AnalysisProfilesV2/Prof_skew_t1", "", {HistType::kTProfile2D, {multAxis, multAxis}});
269+
histos.add("AnalysisProfilesV2/Prof_kurt_t1", "", {HistType::kTProfile2D, {multAxis, multAxis}});
270+
258271
// Analysis Profiles for error
259272
subSampleMcGen.resize(cfgNsubSample);
273+
subSampleMcGenV2.resize(cfgNsubSample);
274+
subSampleV2.resize(cfgNsubSample);
260275
for (int i = 0; i < cfgNsubSample; i++) {
261276
subSampleMcGen[i].resize(4);
277+
subSampleMcGenV2[i].resize(4);
278+
subSampleV2[i].resize(4);
262279
}
263280
for (int i = 0; i < cfgNsubSample; i++) {
264281
subSampleMcGen[i][0] = std::get<std::shared_ptr<TProfile2D>>(histos.add(Form("MCGenerated/AnalysisProfiles/subSample_%d/Prof_mean_t1", i), "", {HistType::kTProfile2D, {centAxis, multAxis}}));
265282
subSampleMcGen[i][1] = std::get<std::shared_ptr<TProfile2D>>(histos.add(Form("MCGenerated/AnalysisProfiles/subSample_%d/Prof_var_t1", i), "", {HistType::kTProfile2D, {centAxis, multAxis}}));
266283
subSampleMcGen[i][2] = std::get<std::shared_ptr<TProfile2D>>(histos.add(Form("MCGenerated/AnalysisProfiles/subSample_%d/Prof_skew_t1", i), "", {HistType::kTProfile2D, {centAxis, multAxis}}));
267284
subSampleMcGen[i][3] = std::get<std::shared_ptr<TProfile2D>>(histos.add(Form("MCGenerated/AnalysisProfiles/subSample_%d/Prof_kurt_t1", i), "", {HistType::kTProfile2D, {centAxis, multAxis}}));
285+
286+
subSampleMcGenV2[i][0] = std::get<std::shared_ptr<TProfile>>(histos.add(Form("MCGenerated/AnalysisProfilesV2/subSample_%d/Prof_mean_t1", i), "", {HistType::kTProfile, {multAxis}}));
287+
subSampleMcGenV2[i][1] = std::get<std::shared_ptr<TProfile>>(histos.add(Form("MCGenerated/AnalysisProfilesV2/subSample_%d/Prof_var_t1", i), "", {HistType::kTProfile, {multAxis}}));
288+
subSampleMcGenV2[i][2] = std::get<std::shared_ptr<TProfile>>(histos.add(Form("MCGenerated/AnalysisProfilesV2/subSample_%d/Prof_skew_t1", i), "", {HistType::kTProfile, {multAxis}}));
289+
subSampleMcGenV2[i][3] = std::get<std::shared_ptr<TProfile>>(histos.add(Form("MCGenerated/AnalysisProfilesV2/subSample_%d/Prof_kurt_t1", i), "", {HistType::kTProfile, {multAxis}}));
290+
291+
subSampleV2[i][0] = std::get<std::shared_ptr<TProfile2D>>(histos.add(Form("AnalysisProfilesV2/subSample_%d/Prof_mean_t1", i), "", {HistType::kTProfile2D, {multAxis, multAxis}}));
292+
subSampleV2[i][1] = std::get<std::shared_ptr<TProfile2D>>(histos.add(Form("AnalysisProfilesV2/subSample_%d/Prof_var_t1", i), "", {HistType::kTProfile2D, {multAxis, multAxis}}));
293+
subSampleV2[i][2] = std::get<std::shared_ptr<TProfile2D>>(histos.add(Form("AnalysisProfilesV2/subSample_%d/Prof_skew_t1", i), "", {HistType::kTProfile2D, {multAxis, multAxis}}));
294+
subSampleV2[i][3] = std::get<std::shared_ptr<TProfile2D>>(histos.add(Form("AnalysisProfilesV2/subSample_%d/Prof_kurt_t1", i), "", {HistType::kTProfile2D, {multAxis, multAxis}}));
268295
}
269296
}
270297

@@ -530,9 +557,11 @@ struct MeanptFluctuationsAnalysis {
530557
if (!mcParticle.has_mcCollision())
531558
continue;
532559

533-
int pdgCode = std::abs(mcParticle.pdgCode());
534-
bool extraPDGType = (pdgCode != PDG_t::kK0Short && pdgCode != PDG_t::kLambda0);
535-
if (extraPDGType && pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != kKPlus && pdgCode != PDG_t::kProton)
560+
// charged check
561+
auto pdgEntry = TDatabasePDG::Instance()->GetParticle(mcParticle.pdgCode());
562+
if (!pdgEntry)
563+
continue;
564+
if (pdgEntry->Charge() == 0)
536565
continue;
537566

538567
if (mcParticle.isPhysicalPrimary()) {
@@ -575,19 +604,29 @@ struct MeanptFluctuationsAnalysis {
575604
histos.fill(HIST("MCGenerated/AnalysisProfiles/Hist2D_Nch_centrality"), cent, nChgen);
576605
histos.fill(HIST("MCGenerated/AnalysisProfiles/Hist2D_meanpt_centrality"), cent, meanTerm1gen);
577606

607+
histos.get<TProfile>(HIST("MCGenerated/AnalysisProfilesV2/Prof_mean_t1"))->Fill(nChgen, meanTerm1gen);
608+
histos.get<TProfile>(HIST("MCGenerated/AnalysisProfilesV2/Prof_var_t1"))->Fill(nChgen, varianceTerm1gen);
609+
histos.get<TProfile>(HIST("MCGenerated/AnalysisProfilesV2/Prof_skew_t1"))->Fill(nChgen, skewnessTerm1gen);
610+
histos.get<TProfile>(HIST("MCGenerated/AnalysisProfilesV2/Prof_kurt_t1"))->Fill(nChgen, kurtosisTerm1gen);
611+
578612
// selecting subsample and filling profiles
579613
float lRandomMc = fRndm->Rndm();
580614
int sampleIndexMc = static_cast<int>(cfgNsubSample * lRandomMc);
581615
subSampleMcGen[sampleIndexMc][0]->Fill(cent, nChgen, meanTerm1gen);
582616
subSampleMcGen[sampleIndexMc][1]->Fill(cent, nChgen, varianceTerm1gen);
583617
subSampleMcGen[sampleIndexMc][2]->Fill(cent, nChgen, skewnessTerm1gen);
584618
subSampleMcGen[sampleIndexMc][3]->Fill(cent, nChgen, kurtosisTerm1gen);
619+
620+
subSampleMcGenV2[sampleIndexMc][0]->Fill(nChgen, meanTerm1gen);
621+
subSampleMcGenV2[sampleIndexMc][1]->Fill(nChgen, varianceTerm1gen);
622+
subSampleMcGenV2[sampleIndexMc][2]->Fill(nChgen, skewnessTerm1gen);
623+
subSampleMcGenV2[sampleIndexMc][3]->Fill(nChgen, kurtosisTerm1gen);
585624
}
586625
//-------------------------------------------------------------------------------------------
587626
}
588627
PROCESS_SWITCH(MeanptFluctuationsAnalysis, processMCGen, "Process Generated MC data", true);
589628

590-
void processMCRec(MyMCRecCollisions::iterator const& collision, MyMCTracks const& tracks, aod::McCollisions const&, aod::McParticles const&)
629+
void processMCRec(MyMCRecCollisions::iterator const& collision, MyMCTracks const& tracks, aod::McCollisions const&, aod::McParticles const& mcParticles)
591630
{
592631
histos.fill(HIST("MCGenerated/hMC"), 5.5);
593632

@@ -635,6 +674,29 @@ struct MeanptFluctuationsAnalysis {
635674
histos.fill(HIST("Hist2D_globalTracks_PVTracks"), collision.multNTracksPV(), tracks.size());
636675
histos.fill(HIST("Hist2D_cent_nch"), tracks.size(), centralityFT0C);
637676

677+
// Calculating generated no of particles for the collision event
678+
double noGen = 0.0;
679+
auto mcColl = collision.mcCollision();
680+
for (const auto& mcParticle : mcParticles) {
681+
if (!mcParticle.has_mcCollision())
682+
continue;
683+
684+
// charged check
685+
auto pdgEntry = TDatabasePDG::Instance()->GetParticle(mcParticle.pdgCode());
686+
if (!pdgEntry)
687+
continue;
688+
if (pdgEntry->Charge() == 0)
689+
continue;
690+
691+
if (mcParticle.isPhysicalPrimary()) {
692+
if ((mcParticle.pt() > cfgCutPtLower) && (mcParticle.pt() < cfgCutPreSelPt) && (std::abs(mcParticle.eta()) < cfgCutPreSelEta)) {
693+
if (mcParticle.pt() > cfgCutPtLower && mcParticle.pt() < cfgCutPtUpper) {
694+
noGen = noGen + 1.0;
695+
}
696+
}
697+
}
698+
} //! end particle loop
699+
638700
// variables
639701
double pTsum = 0.0;
640702
double nN = 0.0;
@@ -667,7 +729,7 @@ struct MeanptFluctuationsAnalysis {
667729
}
668730

669731
histos.fill(HIST("his2DdcaXYvsPtBeforePtDepSel"), track.pt(), track.dcaXY());
670-
histos.fill(HIST("his2DdcaZvsPtBeforePtDepSel"), track.pt(), track.dcaZ());
732+
histos.fill(HIST("his2DdcaZvsPtPtBeforePtDepSel"), track.pt(), track.dcaZ());
671733

672734
if (cfgUsePtDepDCAxy && !(std::abs(track.dcaXY()) < fPtDepDCAxy->Eval(track.pt()))) {
673735
continue;
@@ -706,6 +768,8 @@ struct MeanptFluctuationsAnalysis {
706768
}
707769
} // end track loop
708770

771+
histos.fill(HIST("MCGenerated/hNgenVsNrec"), noGen, nCh);
772+
709773
// MeanPt
710774
if (nN > 0.0f)
711775
histos.fill(HIST("hMeanPt"), cent, pTsum / nN);
@@ -725,13 +789,23 @@ struct MeanptFluctuationsAnalysis {
725789
histos.fill(HIST("AnalysisProfiles/Hist2D_Nch_centrality"), cent, nCh);
726790
histos.fill(HIST("AnalysisProfiles/Hist2D_meanpt_centrality"), cent, meanTerm1);
727791

792+
histos.get<TProfile2D>(HIST("AnalysisProfilesV2/Prof_mean_t1"))->Fill(noGen, nCh, meanTerm1);
793+
histos.get<TProfile2D>(HIST("AnalysisProfilesV2/Prof_var_t1"))->Fill(noGen, nCh, varianceTerm1);
794+
histos.get<TProfile2D>(HIST("AnalysisProfilesV2/Prof_skew_t1"))->Fill(noGen, nCh, skewnessTerm1);
795+
histos.get<TProfile2D>(HIST("AnalysisProfilesV2/Prof_kurt_t1"))->Fill(noGen, nCh, kurtosisTerm1);
796+
728797
// selecting subsample and filling profiles
729798
float lRandom = fRndm->Rndm();
730799
int sampleIndex = static_cast<int>(cfgNsubSample * lRandom);
731800
subSample[sampleIndex][0]->Fill(cent, nCh, meanTerm1);
732801
subSample[sampleIndex][1]->Fill(cent, nCh, varianceTerm1);
733802
subSample[sampleIndex][2]->Fill(cent, nCh, skewnessTerm1);
734803
subSample[sampleIndex][3]->Fill(cent, nCh, kurtosisTerm1);
804+
805+
subSampleV2[sampleIndex][0]->Fill(noGen, nCh, meanTerm1);
806+
subSampleV2[sampleIndex][1]->Fill(noGen, nCh, varianceTerm1);
807+
subSampleV2[sampleIndex][2]->Fill(noGen, nCh, skewnessTerm1);
808+
subSampleV2[sampleIndex][3]->Fill(noGen, nCh, kurtosisTerm1);
735809
}
736810
//-------------------------------------------------------------------------------------------
737811
}
@@ -804,7 +878,7 @@ struct MeanptFluctuationsAnalysis {
804878
}
805879

806880
histos.fill(HIST("his2DdcaXYvsPtBeforePtDepSel"), track.pt(), track.dcaXY());
807-
histos.fill(HIST("his2DdcaZvsPtBeforePtDepSel"), track.pt(), track.dcaZ());
881+
histos.fill(HIST("his2DdcaZvsPtPtBeforePtDepSel"), track.pt(), track.dcaZ());
808882

809883
if (cfgUsePtDepDCAxy && !(std::abs(track.dcaXY()) < fPtDepDCAxy->Eval(track.pt()))) {
810884
continue;

0 commit comments

Comments
 (0)