From 6187dd6643378666d49f937e6ac0c38ac3246195 Mon Sep 17 00:00:00 2001 From: abilandz Date: Fri, 20 Oct 2023 10:42:32 +0200 Subject: [PATCH 1/4] preliminary support for Monte Carlo analysis + support for PROCESS_SWITCH + switch to templates in few function to trim down code bloat --- .../Core/MuPa-Configurables.h | 10 +- .../Core/MuPa-DataMembers.h | 13 +- .../Core/MuPa-Enums.h | 4 +- .../Core/MuPa-MemberFunctions.h | 566 ++++++++++++++++-- .../Tasks/multiparticle-correlations-ab.cxx | 276 +++++---- 5 files changed, 683 insertions(+), 186 deletions(-) diff --git a/PWGCF/MultiparticleCorrelations/Core/MuPa-Configurables.h b/PWGCF/MultiparticleCorrelations/Core/MuPa-Configurables.h index e2abe06f3b5..c99bf50c2e8 100644 --- a/PWGCF/MultiparticleCorrelations/Core/MuPa-Configurables.h +++ b/PWGCF/MultiparticleCorrelations/Core/MuPa-Configurables.h @@ -18,19 +18,21 @@ Configurable cfTaskName{ "cfTaskName", "Default task name", "set task name - use eventually to determine weights for this task"}; -Configurable cfRunNumber{ - "cfRunNumber", "Some run number", - "set run number - TBI temporarily here, this shall be eventually " - "automatically obtained from e.g. collision.bc().runNumber()"}; Configurable cfVerbose{ "cfVerbose", false, "run or not in verbose mode (but not for function calls per particle)"}; Configurable cfVerboseForEachParticle{ "cfVerboseForEachParticle", false, "run or not in verbose mode (also for function calls per particle)"}; +Configurable cfDoAdditionalInsanityChecks{ + "cfDoAdditionalInsanityChecks", false, + "do additional insanity checks at run time (this leads to small loss of performance)"}; Configurable cfUseCCDB{ "cfUseCCDB", true, "access personal files from CCDB or from home dir in AliEn"}; +Configurable cfWhatToProcess{ + "cfWhatToProcess", "Rec", + "Rec = process only reconstructed, Sim = process only simulated, RecSim = process both reconstructed and simulated"}; // Test0: Configurable cfCalculateTest0{"cfCalculateTest0", false, diff --git a/PWGCF/MultiparticleCorrelations/Core/MuPa-DataMembers.h b/PWGCF/MultiparticleCorrelations/Core/MuPa-DataMembers.h index e1dd0e9f0aa..017b04d5d2f 100644 --- a/PWGCF/MultiparticleCorrelations/Core/MuPa-DataMembers.h +++ b/PWGCF/MultiparticleCorrelations/Core/MuPa-DataMembers.h @@ -42,18 +42,24 @@ UInt_t fRandomSeed = 0; // argument to TRandom3 constructor. By default is 0, TString fTaskName = ""; // task name - this one is used to get the right weights // programatically for this analysis TString fRunNumber = ""; // over which run number this task is executed +Bool_t fRunNumberIsDetermined = + kFALSE; // ensures that run number is determined in process() and propagated to already booked objects only once Bool_t fVerbose = kFALSE; // print additional info like Green(__PRETTY_FUNCTION__); etc., to // be used during debugging, but not for function calls per particle Bool_t fVerboseForEachParticle = kFALSE; // print additional info like Green(__PRETTY_FUNCTION__); etc., to // be used during debugging, also for function calls per particle +Bool_t fDoAdditionalInsanityChecks = + kFALSE; // do additional insanity checks at run time, at the expense of losing a bit of performance + // For instance, check if the run number in the current 'collision' is the same as run number in the first 'collision', etc. Bool_t fUseCCDB = kFALSE; // access personal files from CCDB (kTRUE, this is set as default in // Configurables), or from home dir in AliEn (kFALSE, use with care, // as this is discouraged) -// TString fModusOperandi; // "Rec" = process only reconstructed, "Sim" = -// process only simulated, "RecSim" = process both reconstructed and simulated +Bool_t fProcessRemainingEvents = + kTRUE; // if certain criteria is reached, e.g. max number of processed events, ignore all subsequent events TBI 20231019 I need instead graceful exit, see preamble of MaxNumberOfEvents() +TString fWhatToProcess = "Rec"; // "Rec" = process only reconstructed, "Sim" = process only simulated, "RecSim" = process both reconstructed and simulated // UInt_t fRandomSeed; // argument to TRandom3 constructor. By default it is 0, // use SetRandomSeed(...) to change it Bool_t fUseFisherYates; // use // SetUseFisherYates(kTRUE); in the steering macro to randomize particle indices @@ -141,7 +147,8 @@ struct ParticleWeights_Arrays { } pw_a; // "pw_a" labels an instance of this group of histograms, e.g. // pw_a.fWeightsHist[0] TString fFileWithWeights = - ""; // path to external ROOT file which holds all particle weights + ""; // path to external ROOT file which holds all particle weights +Bool_t fParticleWeightsAreFetched = kFALSE; // ensures that particle weights are fetched only once // *) Nested loops: TList* fNestedLoopsList = NULL; // list to hold all nested loops objects diff --git a/PWGCF/MultiparticleCorrelations/Core/MuPa-Enums.h b/PWGCF/MultiparticleCorrelations/Core/MuPa-Enums.h index 623912b17ad..dc2e9658fc4 100644 --- a/PWGCF/MultiparticleCorrelations/Core/MuPa-Enums.h +++ b/PWGCF/MultiparticleCorrelations/Core/MuPa-Enums.h @@ -21,8 +21,8 @@ enum eConfiguration { eConfiguration_N }; -enum eRecoSim { eRec = 0, - eSim = 1 }; +enum eRecSim { eRec = 0, + eSim = 1 }; enum eBeforeAfter { eBefore = 0, eAfter = 1 }; diff --git a/PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h b/PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h index 843de91ca9a..d937629b85d 100644 --- a/PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h +++ b/PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h @@ -18,12 +18,13 @@ // *) Utility; // a) Called directly in init(...): -// void BookBaseList() +// void BookBaseList(); // void DefaultConfiguration(); // void DefaultBooking(); // void DefaultBinning(); // void DefaultCuts(); // Remark: has to be called after DefaultBinning(), since // some default cuts are defined through default binning, to ease bookeeping +// void WhatToProcess(); // void BookAndNestAllLists() // void BookEventHistograms() // void BookParticleHistograms() @@ -37,18 +38,23 @@ // void BookResultsHistograms() // b) Called directly in process(...): +// void Preprocess(aod::Collision const& collision, eRecSim rs); +// Bool_t MaxNumberOfEvents(const Int_t rs); +// void GetParticleWeights(); +// void DetermineAndPropagateRunNumber(aod::Collision const& collision); +// void CheckCurrentRunNumber(aod::Collision const& collision); +// template void FillEventHistograms(aod::Collision const& collision, T const& tracks, const Int_t rs, eBeforeAfter ba) +// template Bool_t EventCuts(aod::Collision const& collision, T const& tracks) +// template void MainLoopOverParticles(T const& tracks, eRecSim rs) +// template void FillParticleHistograms(T const& track, eRecSim rs, eBeforeAfter ba) +// template void ParticleCuts(T const& track, eRecSim rs) // void ResetEventByEventQuantities(); -// void FillEventHistograms(aod::Collision const& collision, aod::Tracks const& -// tracks, const Int_t rs, const Int_t ba); // reco or sim, before or after -// event cuts Bool_t EventCuts(aod::Collision const& collision, aod::Tracks -// const& tracks) void FillParticleHistograms(aod::Track const& track, const -// Int_t rs, const Int_t ba); // reco or sim, before or after particle cuts -// Bool_t ParticleCuts(aod::Track const& track) + // Double_t Weight(const Double_t &value, const char *variable) -// void CalculateCorrelations(); -// void CalculateNestedLoops(); // calculate all standard isotropic correlations -// with nested loops Double_t CalculateCustomNestedLoop(TArrayI *harmonics); // -// calculate nested loop for the specified harmonics +// void CalculateEverything(); +// void CalculateCorrelations(); +// void CalculateNestedLoops(); // calculate all standard isotropic correlations with nested loops +// Double_t CalculateCustomNestedLoop(TArrayI *harmonics); // calculate nested loop for the specified harmonics // *) Called after all events are processed (former "Terminate()"): // void ComparisonNestedLoopsVsCorrelations(); @@ -115,6 +121,40 @@ void BookBaseList() //============================================================ +void WhatToProcess() +{ + // Set here what to process: only rec, both rec and sim, only sim. + // Use in combination with configurable cfWhatToProcess. + // TBI 20231017 I call this function, but it still has no desired effect, until I can call PROCESS_SWITCH( ) by passing a variable, instead only literals 'true' or 'false', as it is now + + if (fVerbose) { + LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); + } + + if (fWhatToProcess.EqualTo("Rec")) { + gProcessRec = true; + } else if (fWhatToProcess.EqualTo("RecSim")) { + gProcessRecSim = true; + } else if (fWhatToProcess.EqualTo("Sim")) { + gProcessSim = true; + } else { + LOGF(info, "\033[1;32m This option is not supported! fWhatToProcess = %s \033[0m", fWhatToProcess.Data()); + LOGF(fatal, "in function \033[1;31m%s at line %d\033[0m", __PRETTY_FUNCTION__, __LINE__); + } + + // Make sure that only one of these flags is set to kTRUE. + // This is needed, becase these flags are used in PROCESS_SWITCH, + // and if 2 or more are kTRUE, then corresponding process function + // is executed over ALL data, then another process(...) function, etc. + if ((Int_t)gProcessRec + (Int_t)gProcessRecSim + (Int_t)gProcessSim > 1) { + LOGF(info, "\033[1;32m Only one flag can be kTRUE: gProcessRec = %d, gProcessRecSim = %d, gProcessSim = %d \033[0m", (Int_t)gProcessRec, (Int_t)gProcessRecSim, (Int_t)gProcessSim); + LOGF(fatal, "in function \033[1;31m%s at line %d\033[0m", __PRETTY_FUNCTION__, __LINE__); + } + +} // WhatToProcess() + +//============================================================ + void DefaultConfiguration() { // Default task configuration. @@ -139,19 +179,21 @@ void DefaultConfiguration() // Configurable cfTaskName{ ... } fTaskName = TString(cfTaskName); - // Configurable cfRunNumber{ ... } // TBI temporarily here, this shall - // be eventually automatically obtained from e.g. collision.bc().runNumber() - fRunNumber = TString(cfRunNumber); - // Configurable cfVerbose{ ... } fVerbose = cfVerbose; // Configurable cfVerboseForEachParticle{ ... } fVerboseForEachParticle = cfVerboseForEachParticle; + // Configurable cfDoAdditionalInsanityChecks{ ... } + fDoAdditionalInsanityChecks = cfDoAdditionalInsanityChecks; + // Configurable cfUseCCDB{ ... } fUseCCDB = cfUseCCDB; + // Configurable cfWhatToProcess{ ... ) + fWhatToProcess = TString(cfWhatToProcess); + // ... // Configurable cfCalculateTest0{ ... }; @@ -294,7 +336,19 @@ void DefaultBinning() cph_a.fParticleHistogramsBins[eEta][1] = -1.; cph_a.fParticleHistogramsBins[eEta][2] = 1.; - // ctd with etpcNClsCrossedRows, eDCA_xy, eDCA_z + cph_a.fParticleHistogramsBins[etpcNClsCrossedRows][0] = 200; + cph_a.fParticleHistogramsBins[etpcNClsCrossedRows][1] = 0.; + cph_a.fParticleHistogramsBins[etpcNClsCrossedRows][2] = 200.; + + cph_a.fParticleHistogramsBins[eDCA_xy][0] = 2000; + cph_a.fParticleHistogramsBins[eDCA_xy][1] = -10.; + cph_a.fParticleHistogramsBins[eDCA_xy][2] = 10.; + + cph_a.fParticleHistogramsBins[eDCA_z][0] = 2000; + cph_a.fParticleHistogramsBins[eDCA_z][1] = -10.; + cph_a.fParticleHistogramsBins[eDCA_z][2] = 10.; + + // ... } // void DefaultBinning() @@ -485,6 +539,9 @@ void BookEventHistograms() } for (Int_t rs = 0; rs < 2; rs++) // reco/sim { + if ((gProcessRec && rs == eSim) || (gProcessSim && rs == eRec)) { + continue; // if I am analyzing only reconstructed data, do not book histos for simulated, and vice versa. + } for (Int_t ba = 0; ba < 2; ba++) // before/after cuts { ceh_a.fEventHistograms[t][rs][ba] = new TH1D( @@ -545,6 +602,9 @@ void BookParticleHistograms() } for (Int_t rs = 0; rs < 2; rs++) // reco/sim { + if ((gProcessRec && rs == eSim) || (gProcessSim && rs == eRec)) { + continue; // if I am analyzing only reconstructed data, do not book histos for simulated, and vice versa. + } for (Int_t ba = 0; ba < 2; ba++) // before/after cuts { cph_a.fParticleHistograms[t][rs][ba] = new TH1D( @@ -881,7 +941,7 @@ void BookTest0Histograms() } // b) Book placeholder and make sure all labels are stored in the placeholder: - StoreLabelsInPlaceholder(); + this->StoreLabelsInPlaceholder(); if (fTest0LabelsPlaceholder) { fTest0List->Add(fTest0LabelsPlaceholder); } @@ -924,6 +984,95 @@ void BookResultsHistograms() //============================================================ +void Preprocess(aod::Collision const& collision, eRecSim rs) +{ + // Do all thingies before starting to process data (e.g. count number of events, fetch the run number, etc.). + + if (fVerbose) { + LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); + } + + // *) If I reached max number of events, ignore the remaining collisions: + if (MaxNumberOfEvents(rs)) { + fProcessRemainingEvents = kFALSE; + } + + // *) Determine and propagate run number info to already booked objects: + if (!fRunNumberIsDetermined) { + DetermineAndPropagateRunNumber(collision); + } + if (fDoAdditionalInsanityChecks && fRunNumberIsDetermined) { + CheckCurrentRunNumber(collision); + } + + // *) Fetch the weights for this particular run number. Do it only once. + // TBI 20231012 If eventualy I can access programatically run number in init(...) at run time, this shall go there. + if (!fParticleWeightsAreFetched) { + if (pw_a.fUseWeights[wPHI] || pw_a.fUseWeights[wPT] || pw_a.fUseWeights[wETA]) { + GetParticleWeights(); + fParticleWeightsAreFetched = kTRUE; + } + } + + // *) Check if gProcessRec, gProcessRecSim and gProcessSim global flags are in sync with eRecSim rs in this PROCESS_SWITCH: + // TBI 20231020 I need to perform this check only once + implement something also for gProcessRecSim + if (gProcessRec && rs != eRec) { + LOGF(error, "\033[1;33m%s gProcessRec && rs != eRec \033[0m", __PRETTY_FUNCTION__); + LOGF(fatal, "gProcessRec = %d, rs = %d", gProcessRec, eRec); + } else if (gProcessSim && rs != eSim) { + LOGF(error, "\033[1;33m%s gProcessSim && rs != eSim \033[0m", __PRETTY_FUNCTION__); + LOGF(fatal, "gProcessSim = %d, rs = %d", gProcessSim, eSim); + } + +} // void Preprocess(aod::Collision const& collision, eRecSim rs) + +//============================================================ + +void DetermineAndPropagateRunNumber(aod::Collision const& collision) +{ + // Determine and propagate run number info to already booked objects, wherever it's relevant. + // Make sure in process(...) that this function is called only once. + + // TBI 20231018 At the moment I can access run number info only in process(...) via collision.bc().runNumber(), but not in init(...) + // Once I can access run number info in init(...), this function shall be called in init(...), not in process(...) + + // a) Determine run number; + // b) Propagate run number to all booked objects, wherever that info is relevant. + + if (fVerbose) { + LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); + } + + // a) Determine run number: + fRunNumber = Form("%d", collision.bc().runNumber()); + if (fRunNumber.EqualTo("")) { + LOGF(error, "\033[1;33m%s fRunNumber is empty, collision.bc().runNumber() failed...\033[0m", __PRETTY_FUNCTION__); + LOGF(fatal, "collision.bc().runNumber() = %d", collision.bc().runNumber()); + } + fRunNumberIsDetermined = kTRUE; + + // b) Propagate run number to all booked objects, wherever that info is relevant: + fBasePro->GetXaxis()->SetBinLabel(eRunNumber, Form("fRunNumber = %s", fRunNumber.Data())); + // ... + ; + +} // void DetermineAndPropagateRunNumber(aod::Collision const& collision) + +//============================================================ + +void CheckCurrentRunNumber(aod::Collision const& collision) +{ + // Insanity check for the current run number. + + if (!fRunNumber.EqualTo(Form("%d", collision.bc().runNumber()))) { + LOGF(error, "\033[1;33m%s Run number changed within process(). This most likely indicates that a given masterjob is processing 2 or more different runs in one go.\033[0m", __PRETTY_FUNCTION__); + LOGF(fatal, "fRunNumber = %s, collision.bc().runNumber() = %d", fRunNumber.Data(), collision.bc().runNumber()); + } + +} // void CheckCurrentRunNumber(aod::Collision const& collision) + +//============================================================ + void ResetEventByEventQuantities() { // Reset all global event-by-event quantities here: @@ -972,11 +1121,67 @@ void ResetEventByEventQuantities() //============================================================ -// Bool_t EventCuts(aod::Collision const& collision, soa::Join const& tracks) -Bool_t EventCuts(aod::Collision const& collision, aod::Tracks const& tracks) +Bool_t EventCuts(aod::Collision const& collision, TracksRec const& tracksRec) { - // ... + // Event cuts on reconstructed data. + + if (fVerbose) { + LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); + } + + // NumberOfEvents: => cut directly in void process( ... ) + + // TotalMultiplicity: + if ((tracksRec.size() < ceh_a.fEventCuts[eTotalMultiplicity][eMin]) || + (tracksRec.size() > ceh_a.fEventCuts[eTotalMultiplicity][eMax])) { + return kFALSE; + } + + // SelectedTracks: => cut directly in void process( ... ) + + // Centrality: TBI + // if (( TBI < ceh_a.fEventCuts[eCentrality][eMin]) || ( TBI > + // ceh_a.fEventCuts[eCentrality][eMax])) { + // return kFALSE; + // } + + // Vertex_x: + if ((collision.posX() < ceh_a.fEventCuts[eVertex_x][eMin]) || + (collision.posX() > ceh_a.fEventCuts[eVertex_x][eMax])) { + return kFALSE; + } + + // Vertex_y: + if ((collision.posY() < ceh_a.fEventCuts[eVertex_y][eMin]) || + (collision.posY() > ceh_a.fEventCuts[eVertex_y][eMax])) { + return kFALSE; + } + + // Vertex_z: + if ((collision.posZ() < ceh_a.fEventCuts[eVertex_z][eMin]) || + (collision.posZ() > ceh_a.fEventCuts[eVertex_z][eMax])) { + return kFALSE; + } + + // NContributors: + if ((collision.numContrib() < ceh_a.fEventCuts[eNContributors][eMin]) || + (collision.numContrib() > ceh_a.fEventCuts[eNContributors][eMax])) { + return kFALSE; + } + + return kTRUE; + +} // void EventCuts(aod::Collision const& collision, TracksRec const& tracksRec) + +//============================================================ + +template +Bool_t EventCuts(aod::Collision const& collision, T const& tracks) +{ + + // Event cuts on reconstructed and simulated data. + + // TBI 20231019 Do I need an additional argument "eRecSim rs" if (fVerbose) { LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); @@ -1024,43 +1229,65 @@ Bool_t EventCuts(aod::Collision const& collision, aod::Tracks const& tracks) return kTRUE; -} // void EventCuts(aod::Collision const& collision, aod::Tracks const& tracks) +} // template Bool_t EventCuts(aod::Collision const& collision, T const& tracks) //============================================================ -// void FillEventHistograms(aod::Collision const& collision, -// soa::Join const& tracks, const -// Int_t rs, const Int_t ba) // reco or sim, before or after event cuts -void FillEventHistograms( - aod::Collision const& collision, aod::Tracks const& tracks, const Int_t rs, - const Int_t ba) // reco or sim, before or after event cuts +template +void FillEventHistograms(aod::Collision const& collision, T const& tracks, eRecSim rs, eBeforeAfter ba) { - // Fill all event histograms. + // Fill all event histograms for reconstructed or simulated data. + + // a) Common to both reconstructed and simulated; + // b) Only reconstructed; + // c) Only simulated. if (fVerbose) { LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); } + // a) Common to both reconstructed and simulated; + + cout << Form("ceh_a.fEventHistograms[eNumberOfEvents][%d][%d] = ", rs, ba) << ceh_a.fEventHistograms[eNumberOfEvents][rs][ba] << endl; + ceh_a.fEventHistograms[eNumberOfEvents][rs][ba]->Fill(0.5); - ceh_a.fEventHistograms[eTotalMultiplicity][rs][ba]->Fill(tracks.size()); - // SelectedTracks => filled directly in void process( ... ) - // ceh_a.fEventHistograms[eCentrality][rs][ba]->Fill( TBI ); => TBI - // 20231007 not ready yet - ceh_a.fEventHistograms[eVertex_x][rs][ba]->Fill(collision.posX()); - ceh_a.fEventHistograms[eVertex_y][rs][ba]->Fill(collision.posY()); - ceh_a.fEventHistograms[eVertex_z][rs][ba]->Fill(collision.posZ()); - ceh_a.fEventHistograms[eNContributors][rs][ba]->Fill(collision.numContrib()); - -} // void FillEventHistograms(aod::Collision const& collision, aod::Tracks -// const& tracks, const Int_t rs, const Int_t ba); // reco or sim, before or -// after event cuts + + switch (rs) { + // b) Only reconstructed: + case eRec: + ceh_a.fEventHistograms[eTotalMultiplicity][eRec][ba]->Fill(tracks.size()); + // SelectedTracks => filled directly in void process( ... ) + // ceh_a.fEventHistograms[eCentrality][eRec][ba]->Fill( TBI ); => TBI + // 20231007 not ready yet + ceh_a.fEventHistograms[eVertex_x][eRec][ba]->Fill(collision.posX()); + ceh_a.fEventHistograms[eVertex_y][eRec][ba]->Fill(collision.posY()); + ceh_a.fEventHistograms[eVertex_z][eRec][ba]->Fill(collision.posZ()); + ceh_a.fEventHistograms[eNContributors][eRec][ba]->Fill(collision.numContrib()); + break; + + // c) Only simulated: + case eSim: + ceh_a.fEventHistograms[eTotalMultiplicity][eSim][ba]->Fill(tracks.size()); + // SelectedTracks => filled directly in void process( ... ) + // ceh_a.fEventHistograms[eCentrality][eSim][ba]->Fill( TBI ); => TBI + // 20231007 not ready yet + /* + ceh_a.fEventHistograms[eVertex_x][eSim][ba]->Fill(collision.posX()); + ceh_a.fEventHistograms[eVertex_y][eSim][ba]->Fill(collision.posY()); + ceh_a.fEventHistograms[eVertex_z][eSim][ba]->Fill(collision.posZ()); + ceh_a.fEventHistograms[eNContributors][eSim][ba]->Fill(collision.numContrib()); + */ + break; + } // switch(rs) { + +} // template void FillEventHistograms(...) //============================================================ -Bool_t ParticleCuts(aod::Track const& track) -// Bool_t ParticleCuts(auto const& track) +template +Bool_t ParticleCuts(T const& track, eRecSim rs) { - // ... + // Particles cuts. if (fVerboseForEachParticle) { LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); @@ -1071,18 +1298,14 @@ Bool_t ParticleCuts(aod::Track const& track) } return kTRUE; - -} // void ParticleCuts(aod::Track const& tracks) +} //============================================================ -// void FillParticleHistograms(auto const& track, const Int_t rs, const Int_t -// ba) // reco or sim, before or after particle cuts -void FillParticleHistograms( - aod::Track const& track, const Int_t rs, - const Int_t ba) // reco or sim, before or after particle cuts +template +void FillParticleHistograms(T const& track, eRecSim rs, eBeforeAfter ba) { - // Fill all particle histograms. + // Fill all particle histograms for reconstructed and simulated data. if (fVerboseForEachParticle) { LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); @@ -1093,17 +1316,34 @@ void FillParticleHistograms( cph_a.fParticleHistograms[ePt][rs][ba]->Fill(track.pt()); cph_a.fParticleHistograms[eEta][rs][ba]->Fill(track.eta()); - /* + switch (rs) { + // a) Only reconstructed: + case eRec: + cph_a.fParticleHistograms[ePhi][rs][ba]->Fill(track.phi()); + cph_a.fParticleHistograms[ePt][rs][ba]->Fill(track.pt()); + cph_a.fParticleHistograms[eEta][rs][ba]->Fill(track.eta()); + break; + + // b) Only simulated: + case eSim: + /* TBI 20231019 not ready yet + cph_a.fParticleHistograms[ePhi][rs][ba]->Fill(track.phi()); + cph_a.fParticleHistograms[ePt][rs][ba]->Fill(track.pt()); + cph_a.fParticleHistograms[eEta][rs][ba]->Fill(track.eta()); +*/ + break; + } // switch(rs) { + + /* TBI 20231019 use also these + check further // From aod::TracksExtra cph_a.fParticleHistograms[etpcNClsCrossedRows][rs][ba]->Fill(track.tpcNClsCrossedRows()); // From aod::TracksDCA cph_a.fParticleHistograms[eDCA_xy][rs][ba]->Fill(track.dcaXY()); cph_a.fParticleHistograms[eDCA_z][rs][ba]->Fill(track.dcaZ()); - */ +*/ -} // void FillParticleHistograms(aod::Track const& track, const Int_t rs, const -// Int_t ba); // reco or sim, before or after particle cuts +} // template void FillParticleHistograms(T const& track, eRecSim rs, eBeforeAfter ba) //============================================================ @@ -1529,11 +1769,9 @@ TH1D* GetHistogramWithWeights(const char* filePath, const char* runNumber, // c) Determine from filePath if the file in on a local machine, or in home // dir AliEn, or in CCDB: // Algorithm: If filePath begins with "/alice/cern.ch/" then it's in home - // dir AliEn. - // If filePath begins with "/alice-ccdb.cern.ch/" then it's in - // CCDB. Therefore, files in AliEn and CCDB must be specified - // with abs path, for local files both abs and relative paths - // are just fine. + // dir AliEn. If filePath begins with "/alice-ccdb.cern.ch/" then it's in + // CCDB. Therefore, files in AliEn and CCDB must be specified with abs path, + // for local files both abs and relative paths are just fine. Bool_t bFileIsInAliEn = kFALSE; Bool_t bFileIsInCCDB = kFALSE; if (TString(filePath).BeginsWith("/alice/cern.ch/")) { @@ -2015,4 +2253,210 @@ Double_t Weight(const Double_t& value, //============================================================ +void GetParticleWeights() +{ + // Get the particle weights. Call this function only once. + + // TBI 20231012 Here the current working assumption is that: + // a) Corrections do not change within a given run; + // b) Hyperloop proceeses the dataset one masterjob per run number. + // If any of these 2 assumptions are violated, this code will have to be modified. + + if (fVerbose) { + LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); + } + + if (pw_a.fUseWeights[wPHI]) { + TH1D* phiWeights = GetHistogramWithWeights(fFileWithWeights.Data(), fRunNumber.Data(), "phi"); + if (!phiWeights) { + LOGF(fatal, "in function \033[1;31m%s at line %d, phiWeights is NULL. Check the external file %s with particle weights\033[0m", __PRETTY_FUNCTION__, __LINE__, fFileWithWeights.Data()); + } + SetWeightsHist(phiWeights, "phi"); + } + + if (pw_a.fUseWeights[wPT]) { + TH1D* ptWeights = GetHistogramWithWeights(fFileWithWeights.Data(), fRunNumber.Data(), "pt"); + if (!ptWeights) { + LOGF(fatal, "in function \033[1;31m%s at line %d, ptWeights is NULL. Check the external file %s with particle weights\033[0m", __PRETTY_FUNCTION__, __LINE__, fFileWithWeights.Data()); + } + SetWeightsHist(ptWeights, "pt"); + } + + if (pw_a.fUseWeights[wETA]) { + TH1D* etaWeights = GetHistogramWithWeights(fFileWithWeights.Data(), fRunNumber.Data(), "eta"); + if (!etaWeights) { + LOGF(fatal, "in function \033[1;31m%s at line %d, etaWeights is NULL. Check the external file %s with particle weights\033[0m", __PRETTY_FUNCTION__, __LINE__, fFileWithWeights.Data()); + } + SetWeightsHist(etaWeights, "eta"); + } + +} // void GetParticleWeights() + +//============================================================ + +Bool_t MaxNumberOfEvents(const Int_t rs) +{ + // Check if max number of events was reached. See also configurable cNumberOfEvents_max. + + if (fVerbose) { + LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); + } + + Bool_t reachedMaxNumberOfEvents = kFALSE; + + if (ceh_a.fEventHistograms[eNumberOfEvents][rs][eAfter]->GetBinContent(1) >= ceh_a.fEventCuts[eNumberOfEvents][eMax]) { + reachedMaxNumberOfEvents = kTRUE; + } + + return reachedMaxNumberOfEvents; + +} // void MaxNumberOfEvents() + +//============================================================ + +template +void MainLoopOverParticles(T const& tracks, eRecSim rs) +{ + // This is the main loop over particles, in which Q-vectors and particle histograms are filled, particle cuts applied, etc. + // If I want to access both rec and sim info in this loop, set gProcessRecSim = true via configurable "cfWhatToProcess" + + if (fVerbose) { + LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); + } + + Double_t dPhi = 0., wPhi = 1.; // azimuthal angle and corresponding phi weight + Double_t dPt = 0., wPt = 1.; // transverse momentum and corresponding pT weight + Double_t dEta = 0., wEta = 1.; // pseudorapidity and corresponding eta weight + Double_t wToPowerP = 1.; // weight raised to power p + fSelectedTracks = 0; // reset number of selected tracks + for (auto& track : tracks) { + + // *) Fill particle histograms for reconstructed particle before particle cuts: + FillParticleHistograms(track, rs, eBefore); + + // *) Fill particle histograms for corresponding simulated particle before particle cuts: + if (gProcessRecSim) { + if (track.has_mcParticle()) { + auto mcParticle = track.mcParticle(); + FillParticleHistograms(mcParticle, eSim, eBefore); // great, this works straight like this, thanks O2 people! + } else { + LOGF(debug, "in function \033[1;31m%s at line %d, track.has_mcParticle() is false. \033[0m", __PRETTY_FUNCTION__, __LINE__); + } + } // if (gProcessRecSim) { + + // *) Particle cuts: + if (!ParticleCuts(track, rs)) { + continue; + } + + // *) Fill particle histograms for reconstructed data after particle cuts: + FillParticleHistograms(track, rs, eAfter); + + // *) Fill particle histograms for corresponding simulated particle after particle cuts: + // TBI 20231020 there is a duplication of code here, see above. I should get 'mcParticle' only once for this 'track' in this loop + if (gProcessRecSim) { + if (track.has_mcParticle()) { + auto mcParticle = track.mcParticle(); + FillParticleHistograms(mcParticle, eSim, eAfter); // great, this works straight like this, thanks O2 people! + } else { + LOGF(debug, "in function \033[1;31m%s at line %d, track.has_mcParticle() is false. \033[0m", __PRETTY_FUNCTION__, __LINE__); + } + } // if (gProcessRecSim) { + + // *) Fill Q-vectors: + if (gProcessRec || gProcessRecSim) { + // Take kinematics from reconstructed particles: + dPhi = track.phi(); + dPt = track.pt(); + dEta = track.eta(); + } else if (gProcessSim) { + // Take kinematics from simulated particles: + dPhi = track.mcParticle().phi(); + dPt = track.mcParticle().pt(); + dEta = track.mcParticle().eta(); + } + + // Particle weights: + if (pw_a.fUseWeights[wPHI]) { + wPhi = Weight(dPhi, "phi"); // corresponding phi weight + if (!(wPhi > 0.)) { + LOGF(error, "\033[1;33m%s wPhi is not positive, skipping this particle for the time being...\033[0m", __PRETTY_FUNCTION__); + LOGF(error, "dPhi = %f\nwPhi = %f", dPhi, wPhi); + continue; + } + } // if(pw_a.fUseWeights[wPHI]) + if (pw_a.fUseWeights[wPT]) { + wPt = Weight(dPt, "pt"); // corresponding pt weight + if (!(wPt > 0.)) { + LOGF(error, "\033[1;33m%s wPt is not positive, skipping this particle for the time being...\033[0m", __PRETTY_FUNCTION__); + LOGF(error, "dPt = %f\nwPt = %f", dPt, wPt); + continue; + } + } // if(pw_a.fUseWeights[wPT]) + if (pw_a.fUseWeights[wETA]) { + wEta = Weight(dEta, "eta"); // corresponding eta weight + if (!(wEta > 0.)) { + LOGF(error, "\033[1;33m%s wEta is not positive, skipping this particle for the time being...\033[0m", __PRETTY_FUNCTION__); + LOGF(error, "dEta = %f\nwEta = %f", dEta, wEta); + continue; + } + } // if(pw_a.fUseWeights[wETA]) + + for (Int_t h = 0; h < gMaxHarmonic * gMaxCorrelator + 1; h++) { + for (Int_t wp = 0; wp < gMaxCorrelator + 1; wp++) { // weight power + if (pw_a.fUseWeights[wPHI] || pw_a.fUseWeights[wPT] || pw_a.fUseWeights[wETA]) { + wToPowerP = pow(wPhi * wPt * wEta, wp); + } + qv_a.fQvector[h][wp] += TComplex(wToPowerP * TMath::Cos(h * dPhi), wToPowerP * TMath::Sin(h * dPhi)); + } // for(Int_t wp=0;wpAddAt(dPhi, fSelectedTracks); + } // remember that the 2nd argument here must start from 0 + if (nl_a.ftaNestedLoops[1]) { + nl_a.ftaNestedLoops[1]->AddAt(wPhi * wPt * wEta, fSelectedTracks); + } // remember that the 2nd argument here must start from 0 + } // if(fCalculateNestedLoops||fCalculateCustomNestedLoop) + + // *) Counter of selected tracks in the current event: + fSelectedTracks++; + if (fSelectedTracks >= cSelectedTracks_max) { + break; + } + + } // for (auto& track : tracks) + +} // template void MainLoopOverParticles(T const& tracks) + +//============================================================ + +void CalculateEverything() +{ + // Calculate everything for selected events and particles. + // Remark: Data members for Q-vectors, containers for nested loops, etc., must all be filled when this function is called. + + if (fVerbose) { + LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); + } + + // *) Calculate multiparticle correlations (standard, isotropic, same harmonic): + if (fCalculateCorrelations) { + CalculateCorrelations(); + } + + // *) Calculate nested loops: + if (fCalculateNestedLoops) { + CalculateNestedLoops(); + + // TBI 20220823 this shall be called after all events are processed, only temporarily here called for each event: + ComparisonNestedLoopsVsCorrelations(); + } + +} // void CalculateEverything() + +//============================================================ + #endif // PWGCF_MULTIPARTICLECORRELATIONS_CORE_MUPA_MEMBERFUNCTIONS_H_ diff --git a/PWGCF/MultiparticleCorrelations/Tasks/multiparticle-correlations-ab.cxx b/PWGCF/MultiparticleCorrelations/Tasks/multiparticle-correlations-ab.cxx index cd8af99ec71..d474b30b975 100644 --- a/PWGCF/MultiparticleCorrelations/Tasks/multiparticle-correlations-ab.cxx +++ b/PWGCF/MultiparticleCorrelations/Tasks/multiparticle-correlations-ab.cxx @@ -13,10 +13,20 @@ #include #include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" -//#include "Common/DataModel/TrackSelectionTables.h" TBI 20231008 needed for aod::TracksDCA +#include "Framework/AnalysisDataModel.h" +#include "Common/DataModel/TrackSelectionTables.h" // needed for aod::TracksDCA table using namespace o2; using namespace o2::framework; +using TracksRec = soa::Join; +using TrackRec = soa::Join::iterator; + +using TracksRecSim = soa::Join; +using TrackRecSim = soa::Join::iterator; + +using TracksSim = soa::Join; // TBI 20231018 only temporarily defined this way, check still aod::McParticles +using TrackSim = soa::Join::iterator; // TBI 20231018 only temporarily defined this way check still aod::McParticles + // ROOT: #include "TList.h" #include "TSystem.h" @@ -34,6 +44,13 @@ using namespace std; // *) Global constants: #include "PWGCF/MultiparticleCorrelations/Core/MuPa-GlobalConstants.h" +// *) These are indended flags for PROCESS_SWITCH, have to be global, at least for the time being... +// TBI 20231017 check this further, it doesn't work yet this way. It seems I have to pass to PROCESS_SWITCH( ) only literals 'true' or 'false' +// TBI 20231020 I could as well re-define them as data members... +bool gProcessRec = false; +bool gProcessRecSim = false; +bool gProcessSim = false; + // *) Main task: struct MultiparticleCorrelationsAB // this name is used in lower-case format to name the TDirectoryFile in AnalysisResults.root { @@ -67,33 +84,22 @@ struct MultiparticleCorrelationsAB // this name is used in lower-case format to Bool_t oldHistAddStatus = TH1::AddDirectoryStatus(); TH1::AddDirectory(kFALSE); - // *) Book base list: - BookBaseList(); - // *) Default configuration, booking, binning and cuts: DefaultConfiguration(); DefaultBooking(); DefaultBinning(); DefaultCuts(); // Remark: has to be called after DefaultBinning(), since some default cuts are defined through default binning, to ease bookeeping - // *) Particle weights: - if (pw_a.fUseWeights[wPHI]) { - TH1D* phiWeights = GetHistogramWithWeights(fFileWithWeights.Data(), fRunNumber.Data(), "phi"); - SetWeightsHist(phiWeights, "phi"); - } - if (pw_a.fUseWeights[wPT]) { - TH1D* ptWeights = GetHistogramWithWeights(fFileWithWeights.Data(), fRunNumber.Data(), "pt"); - SetWeightsHist(ptWeights, "pt"); - } - if (pw_a.fUseWeights[wETA]) { - TH1D* etaWeights = GetHistogramWithWeights(fFileWithWeights.Data(), fRunNumber.Data(), "eta"); - SetWeightsHist(etaWeights, "eta"); - } + // *) Set what to process - only rec, both rec and sim, only sim: + WhatToProcess(); // yes, this can be called here, after calling all Default* member functions above, because this has an effect only on Book* members functions // *) Book random generator: delete gRandom; gRandom = new TRandom3(fRandomSeed); // if uiSeed is 0, the seed is determined uniquely in space and time via TUUID + // *) Book base list: + BookBaseList(); + // *) Book all remaining objects; BookAndNestAllLists(); BookEventHistograms(); @@ -112,129 +118,167 @@ struct MultiparticleCorrelationsAB // this name is used in lower-case format to // ------------------------------------------- - // *) Process the data: - // void process(aod::Collision const& collision, soa::Join const& tracks) // called once per collision found in the time frame - void process(aod::Collision const& collision, aod::Tracks const& tracks) // called once per collision found in the time frame + // Since I am subscribing to different tables in each case, there are 3 separate implementations of process(...) + // A) Process only reconstructed data; + // B) Process both reconstructed and simulated data; + // C) Process only simulated data. + // TBI 20231017 There is a (minor) code bloat here obviously. Fix this when it becomes possible to post the last argument to PROCESS_SWITCH via variable (ideally Configurable) + + // A) Process only reconstructed data: + void processRec(aod::Collision const& collision, aod::BCs const&, TracksRec const& tracksRec) // called once per collision found in the time frame { + // ... - // *) TBI 20231008 Temporarily here: If I reached max number of events, ignore the remaining collisions. - // But what I really need here is a graceful exit from subsequent processing (which will also dump the output file, etc.) - if (ceh_a.fEventHistograms[eNumberOfEvents][eRec][eAfter]->GetBinContent(1) >= ceh_a.fEventCuts[eNumberOfEvents][eMax]) { - return; + // *) If I reached max number of events, ignore the remaining collisions: + if (!fProcessRemainingEvents) { + return; // TBI 20231008 Temporarily implemented this way. But what I really need here is a graceful exit from subsequent processing (which will also dump the output file, etc.) + } + + // *) TBI 20231020 Temporary here (use configurable 'cfWhatToProcess' to set this flag corerectly): + if (!gProcessRec) { + LOGF(fatal, "in function \033[1;31m%s at line %d\033[0m", __PRETTY_FUNCTION__, __LINE__); } + // *) Do all thingies before starting to process data (e.g. count number of events, fetch the run number, etc.): + Preprocess(collision, eRec); + // *) Fill event histograms for reconstructed data before event cuts: - FillEventHistograms(collision, tracks, eRec, eBefore); + FillEventHistograms(collision, tracksRec, eRec, eBefore); // *) Event cuts: - if (!EventCuts(collision, tracks)) { + if (!EventCuts(collision, tracksRec)) { return; } // *) Fill event histograms for reconstructed data after event cuts: - FillEventHistograms(collision, tracks, eRec, eAfter); - - // *) Main loop over particles: - Double_t dPhi = 0., wPhi = 1.; // azimuthal angle and corresponding phi weight - Double_t dPt = 0., wPt = 1.; // transverse momentum and corresponding pT weight - Double_t dEta = 0., wEta = 1.; // pseudorapidity and corresponding eta weight - Double_t wToPowerP = 1.; // weight raised to power p - fSelectedTracks = 0; // reset number of selected tracks - for (auto& track : tracks) { - - // *) Fill particle histograms for reconstructed data before particle cuts: - FillParticleHistograms(track, eRec, eBefore); - - // *) Particle cuts: - if (!ParticleCuts(track)) { - continue; - } - - // *) Fill particle histograms for reconstructed data after particle cuts: - FillParticleHistograms(track, eRec, eAfter); - - // *) Fill Q-vectors: - dPhi = track.phi(); - dPt = track.pt(); - dEta = track.eta(); - // Particle weights: - if (pw_a.fUseWeights[wPHI]) { - wPhi = Weight(dPhi, "phi"); // corresponding phi weight - if (!(wPhi > 0.)) { - LOGF(error, "\033[1;33m%s wPhi is not positive, skipping this particle for the time being...\033[0m", __PRETTY_FUNCTION__); - LOGF(error, "dPhi = %f\nwPhi = %f", dPhi, wPhi); - continue; - } - } // if(pw_a.fUseWeights[wPHI]) - if (pw_a.fUseWeights[wPT]) { - wPt = Weight(dPt, "pt"); // corresponding pt weight - if (!(wPt > 0.)) { - LOGF(error, "\033[1;33m%s wPt is not positive, skipping this particle for the time being...\033[0m", __PRETTY_FUNCTION__); - LOGF(error, "dPt = %f\nwPt = %f", dPt, wPt); - continue; - } - } // if(pw_a.fUseWeights[wPT]) - if (pw_a.fUseWeights[wETA]) { - wEta = Weight(dEta, "eta"); // corresponding eta weight - if (!(wEta > 0.)) { - LOGF(error, "\033[1;33m%s wEta is not positive, skipping this particle for the time being...\033[0m", __PRETTY_FUNCTION__); - LOGF(error, "dEta = %f\nwEta = %f", dEta, wEta); - continue; - } - } // if(pw_a.fUseWeights[wETA]) - - for (Int_t h = 0; h < gMaxHarmonic * gMaxCorrelator + 1; h++) { - for (Int_t wp = 0; wp < gMaxCorrelator + 1; wp++) { // weight power - if (pw_a.fUseWeights[wPHI] || pw_a.fUseWeights[wPT] || pw_a.fUseWeights[wETA]) { - wToPowerP = pow(wPhi * wPt * wEta, wp); - } - qv_a.fQvector[h][wp] += TComplex(wToPowerP * TMath::Cos(h * dPhi), wToPowerP * TMath::Sin(h * dPhi)); - } // for(Int_t wp=0;wpAddAt(dPhi, fSelectedTracks); - } // remember that the 2nd argument here must start from 0 - if (nl_a.ftaNestedLoops[1]) { - nl_a.ftaNestedLoops[1]->AddAt(wPhi * wPt * wEta, fSelectedTracks); - } // remember that the 2nd argument here must start from 0 - } // if(fCalculateNestedLoops||fCalculateCustomNestedLoop) - - // *) Counter of selected tracks in the current event: - fSelectedTracks++; - if (fSelectedTracks >= cSelectedTracks_max) { - break; - } - - } // for (auto& track : tracks) - - // *) Fill remaining event histograms for reconstructed data after event (and particle) cuts: + FillEventHistograms(collision, tracksRec, eRec, eAfter); + + // *) Main loop over reconstructed particles: + //MainLoopOverParticles(tracksRec, eRec); TBI 20231020 yes, this is a bug, becase I call track.has_mcParticle() there, etc. Re-think how to re-implement this + + // *) Fill remaining event histograms for reconstructed data after event AND after particle cuts: ceh_a.fEventHistograms[eSelectedTracks][eRec][eAfter]->Fill(fSelectedTracks); - // *) Remaining event cuts: + // *) Remaining event cuts which can be applied only after the loop over particles is performed: if ((fSelectedTracks < ceh_a.fEventCuts[eSelectedTracks][eMin]) || (fSelectedTracks > ceh_a.fEventCuts[eSelectedTracks][eMax])) { return; } - // *) Calculate multiparticle correlations (standard, isotropic, same harmonic): - if (fCalculateCorrelations) { - CalculateCorrelations(); + // *) Calculate everything for selected events and particles: + CalculateEverything(); + + // *) Reset event-by-event quantities: + ResetEventByEventQuantities(); + + } // void processRec(...) + PROCESS_SWITCH(MultiparticleCorrelationsAB, processRec, "process only reconstructed information", false); + + // ------------------------------------------- + + // B) Process both reconstructed and simulated data: + void processRecSim(aod::Collision const& collision, aod::BCs const&, TracksRecSim const& tracksRecSim, aod::McParticles const&) // called once per collision found in the time frame + { + // Remark: Implemented separately to processRec(...), because here I need to subscribe also to Monte Carlo tables. + + // *) If I reached max number of events, ignore the remaining collisions: + if (!fProcessRemainingEvents) { + return; // TBI 20231008 Temporarily implemented this way. But what I really need here is a graceful exit from subsequent processing (which will also dump the output file, etc.) + } + + // *) TBI 20231020 Temporary here (use configurable 'cfWhatToProcess' to set this flag corerectly): + if (!gProcessRecSim) { + LOGF(fatal, "in function \033[1;31m%s at line %d\033[0m", __PRETTY_FUNCTION__, __LINE__); } - // *) Calculate nested loops: - if (fCalculateNestedLoops) { - CalculateNestedLoops(); + // *) Do all thingies before starting to process data (e.g. count number of events, fetch the run number, etc.): + Preprocess(collision, eRec); - // TBI 20220823 this shall be called after all events are processed, only temporarily here called for each event: - ComparisonNestedLoopsVsCorrelations(); + // *) Fill event histograms for reconstructed and simulated data before event cuts: + FillEventHistograms(collision, tracksRecSim, eRec, eBefore); + FillEventHistograms(collision, tracksRecSim, eSim, eBefore); + + // *) Event cuts: + if (!EventCuts(collision, tracksRecSim)) { + return; } - // *) Reset event-by-event objects: + // *) Fill event histograms for reconstructed and simulated data after event cuts: + FillEventHistograms(collision, tracksRecSim, eRec, eAfter); + FillEventHistograms(collision, tracksRecSim, eSim, eAfter); + + // *) Main loop over reconstructed particles: + // For simulated particles there is NO separate loop, if flag gProcessRecSim = true (via configurable "cfWhatToProcess"), MC info is accessed via label. + MainLoopOverParticles(tracksRecSim, eRec); + + // *) Fill remaining event histograms for reconstructed and simulated data after event AND after particle cuts: + ceh_a.fEventHistograms[eSelectedTracks][eRec][eAfter]->Fill(fSelectedTracks); + // ceh_a.fEventHistograms[eSelectedTracks][eSim][eAfter]->Fill(fSelectedTracks); // TBI 20231020 do I really need this one? + + // *) Remaining event cuts which can be applied only after the loop over particles is performed: + if ((fSelectedTracks < ceh_a.fEventCuts[eSelectedTracks][eMin]) || (fSelectedTracks > ceh_a.fEventCuts[eSelectedTracks][eMax])) { + return; + } + + // *) Calculate everything for selected events and particles: + CalculateEverything(); + + // *) Reset event-by-event quantities: + ResetEventByEventQuantities(); + + } // void processRecSim(...) + PROCESS_SWITCH(MultiparticleCorrelationsAB, processRecSim, "process both reconstructed and simulated information", true); + + // ------------------------------------------- + + // C) Process only simulated data: + void processSim(aod::Collision const& collision, aod::BCs const&, TracksSim const& tracksSim, aod::McParticles const&) // called once per collision found in the time frame + { + // ... + + // *) If I reached max number of events, ignore the remaining collisions: + if (!fProcessRemainingEvents) { + return; // TBI 20231008 Temporarily implemented this way. But what I really need here is a graceful exit from subsequent processing (which will also dump the output file, etc.) + } + + // *) TBI 20231020 Temporary here (use configurable 'cfWhatToProcess' to set this flag corerectly): + if (!gProcessSim) { + LOGF(fatal, "in function \033[1;31m%s at line %d\033[0m", __PRETTY_FUNCTION__, __LINE__); + } + + // *) Do all thingies before starting to process data (e.g. count number of events, fetch the run number, etc.): + Preprocess(collision, eSim); + + // *) Fill event histograms for simulated data before event cuts: + FillEventHistograms(collision, tracksSim, eRec, eBefore); + + // *) Event cuts: + if (!EventCuts(collision, tracksSim)) { + return; + } + + // *) Fill event histograms for simulated data after event cuts: + FillEventHistograms(collision, tracksSim, eSim, eAfter); + + // *) Main loop over reconstructed particles: + // For simulated particles there is NO separate loop, if flag gProcessRecSim = true (via configurable "cfWhatToProcess"), MC info is accessed via label. + // MainLoopOverParticles(tracksSim, eSim); // TBI 20231020 it will NOT work this way + re-think the definition of tracksSim in the preamble + + // *) Fill remaining event histograms for reconstructed and simulated data after event AND after particle cuts: + ceh_a.fEventHistograms[eSelectedTracks][eSim][eAfter]->Fill(fSelectedTracks); + + // *) Remaining event cuts which can be applied only after the loop over particles is performed: + if ((fSelectedTracks < ceh_a.fEventCuts[eSelectedTracks][eMin]) || (fSelectedTracks > ceh_a.fEventCuts[eSelectedTracks][eMax])) { + return; + } + + // *) Calculate everything for selected events and particles: + CalculateEverything(); + + // *) Reset event-by-event quantities: ResetEventByEventQuantities(); - } // void process(...) + } // void processSim(...) + PROCESS_SWITCH(MultiparticleCorrelationsAB, processSim, "process only simulated information", false); }; // struct MultiparticleCorrelationsAB From 7f4acb1f1ff00953eb9b6fbd1b3f7c6fc6450796 Mon Sep 17 00:00:00 2001 From: abilandz Date: Fri, 20 Oct 2023 10:49:16 +0200 Subject: [PATCH 2/4] fix for formatter #1 --- .../Tasks/multiparticle-correlations-ab.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGCF/MultiparticleCorrelations/Tasks/multiparticle-correlations-ab.cxx b/PWGCF/MultiparticleCorrelations/Tasks/multiparticle-correlations-ab.cxx index d474b30b975..55e536e2d6f 100644 --- a/PWGCF/MultiparticleCorrelations/Tasks/multiparticle-correlations-ab.cxx +++ b/PWGCF/MultiparticleCorrelations/Tasks/multiparticle-correlations-ab.cxx @@ -154,7 +154,7 @@ struct MultiparticleCorrelationsAB // this name is used in lower-case format to FillEventHistograms(collision, tracksRec, eRec, eAfter); // *) Main loop over reconstructed particles: - //MainLoopOverParticles(tracksRec, eRec); TBI 20231020 yes, this is a bug, becase I call track.has_mcParticle() there, etc. Re-think how to re-implement this + // MainLoopOverParticles(tracksRec, eRec); TBI 20231020 yes, this is a bug, becase I call track.has_mcParticle() there, etc. Re-think how to re-implement this // *) Fill remaining event histograms for reconstructed data after event AND after particle cuts: ceh_a.fEventHistograms[eSelectedTracks][eRec][eAfter]->Fill(fSelectedTracks); From 1d29fa498da658438927a326279b1277c4007459 Mon Sep 17 00:00:00 2001 From: abilandz Date: Sat, 21 Oct 2023 13:17:21 +0200 Subject: [PATCH 3/4] implementing Victor comments --- .../Core/MuPa-Enums.h | 3 +- .../Core/MuPa-MemberFunctions.h | 502 +++++++++++------- .../Tasks/multiparticle-correlations-ab.cxx | 138 ++--- 3 files changed, 353 insertions(+), 290 deletions(-) diff --git a/PWGCF/MultiparticleCorrelations/Core/MuPa-Enums.h b/PWGCF/MultiparticleCorrelations/Core/MuPa-Enums.h index dc2e9658fc4..1a9c8b715c5 100644 --- a/PWGCF/MultiparticleCorrelations/Core/MuPa-Enums.h +++ b/PWGCF/MultiparticleCorrelations/Core/MuPa-Enums.h @@ -22,7 +22,8 @@ enum eConfiguration { }; enum eRecSim { eRec = 0, - eSim = 1 }; + eSim = 1, + eRecAndSim = 2 }; // TBI 20231021 find a better name? enum eBeforeAfter { eBefore = 0, eAfter = 1 }; diff --git a/PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h b/PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h index d937629b85d..fafe9e0cb24 100644 --- a/PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h +++ b/PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h @@ -38,17 +38,22 @@ // void BookResultsHistograms() // b) Called directly in process(...): -// void Preprocess(aod::Collision const& collision, eRecSim rs); -// Bool_t MaxNumberOfEvents(const Int_t rs); -// void GetParticleWeights(); -// void DetermineAndPropagateRunNumber(aod::Collision const& collision); -// void CheckCurrentRunNumber(aod::Collision const& collision); -// template void FillEventHistograms(aod::Collision const& collision, T const& tracks, const Int_t rs, eBeforeAfter ba) -// template Bool_t EventCuts(aod::Collision const& collision, T const& tracks) -// template void MainLoopOverParticles(T const& tracks, eRecSim rs) -// template void FillParticleHistograms(T const& track, eRecSim rs, eBeforeAfter ba) -// template void ParticleCuts(T const& track, eRecSim rs) -// void ResetEventByEventQuantities(); +// void SteerRec(CollisionRec const& collision, TracksRec const& tracks) +// void SteerRecSim(CollisionRec const& collision, TracksRecSim const& tracks) +// void SteerSim(CollisionSim const& collision, TracksSim const& tracks) +// template void Preprocess(T const& collision) +// Bool_t MaxNumberOfEvents(const Int_t rs); +// void GetParticleWeights(); +// template void DetermineAndPropagateRunNumber(T const& collision) +// template void CheckCurrentRunNumber(T const& collision); +// template void FillEventHistograms(T1 const& collision, T2 const& tracks, eBeforeAfter ba) +// template Bool_t EventCuts(T1 const& collision, T2 const& tracks) +// void MainLoopOverParticlesRec(TracksRec const& tracks) +// void MainLoopOverParticlesRecSim(TracksRecSim const& tracks) +// void MainLoopOverParticlesSim(TracksSim const& tracks) +// template void FillParticleHistograms(T const& track, eBeforeAfter ba) +// template void ParticleCuts(T const& track, eRecSim rs) +// void ResetEventByEventQuantities(); // Double_t Weight(const Double_t &value, const char *variable) // void CalculateEverything(); @@ -984,7 +989,8 @@ void BookResultsHistograms() //============================================================ -void Preprocess(aod::Collision const& collision, eRecSim rs) +template +void Preprocess(T const& collision) { // Do all thingies before starting to process data (e.g. count number of events, fetch the run number, etc.). @@ -993,7 +999,7 @@ void Preprocess(aod::Collision const& collision, eRecSim rs) } // *) If I reached max number of events, ignore the remaining collisions: - if (MaxNumberOfEvents(rs)) { + if (MaxNumberOfEvents()) { fProcessRemainingEvents = kFALSE; } @@ -1014,21 +1020,12 @@ void Preprocess(aod::Collision const& collision, eRecSim rs) } } - // *) Check if gProcessRec, gProcessRecSim and gProcessSim global flags are in sync with eRecSim rs in this PROCESS_SWITCH: - // TBI 20231020 I need to perform this check only once + implement something also for gProcessRecSim - if (gProcessRec && rs != eRec) { - LOGF(error, "\033[1;33m%s gProcessRec && rs != eRec \033[0m", __PRETTY_FUNCTION__); - LOGF(fatal, "gProcessRec = %d, rs = %d", gProcessRec, eRec); - } else if (gProcessSim && rs != eSim) { - LOGF(error, "\033[1;33m%s gProcessSim && rs != eSim \033[0m", __PRETTY_FUNCTION__); - LOGF(fatal, "gProcessSim = %d, rs = %d", gProcessSim, eSim); - } - -} // void Preprocess(aod::Collision const& collision, eRecSim rs) +} // template void Preprocess(T const& collision) //============================================================ -void DetermineAndPropagateRunNumber(aod::Collision const& collision) +template +void DetermineAndPropagateRunNumber(T const& collision) { // Determine and propagate run number info to already booked objects, wherever it's relevant. // Make sure in process(...) that this function is called only once. @@ -1054,13 +1051,13 @@ void DetermineAndPropagateRunNumber(aod::Collision const& collision) // b) Propagate run number to all booked objects, wherever that info is relevant: fBasePro->GetXaxis()->SetBinLabel(eRunNumber, Form("fRunNumber = %s", fRunNumber.Data())); // ... - ; -} // void DetermineAndPropagateRunNumber(aod::Collision const& collision) +} // template void DetermineAndPropagateRunNumber(T const& collision) //============================================================ -void CheckCurrentRunNumber(aod::Collision const& collision) +template +void CheckCurrentRunNumber(T const& collision) { // Insanity check for the current run number. @@ -1069,7 +1066,7 @@ void CheckCurrentRunNumber(aod::Collision const& collision) LOGF(fatal, "fRunNumber = %s, collision.bc().runNumber() = %d", fRunNumber.Data(), collision.bc().runNumber()); } -} // void CheckCurrentRunNumber(aod::Collision const& collision) +} // template void CheckCurrentRunNumber(aod::Collision const& collision) //============================================================ @@ -1121,62 +1118,8 @@ void ResetEventByEventQuantities() //============================================================ -Bool_t EventCuts(aod::Collision const& collision, TracksRec const& tracksRec) -{ - // Event cuts on reconstructed data. - - if (fVerbose) { - LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); - } - - // NumberOfEvents: => cut directly in void process( ... ) - - // TotalMultiplicity: - if ((tracksRec.size() < ceh_a.fEventCuts[eTotalMultiplicity][eMin]) || - (tracksRec.size() > ceh_a.fEventCuts[eTotalMultiplicity][eMax])) { - return kFALSE; - } - - // SelectedTracks: => cut directly in void process( ... ) - - // Centrality: TBI - // if (( TBI < ceh_a.fEventCuts[eCentrality][eMin]) || ( TBI > - // ceh_a.fEventCuts[eCentrality][eMax])) { - // return kFALSE; - // } - - // Vertex_x: - if ((collision.posX() < ceh_a.fEventCuts[eVertex_x][eMin]) || - (collision.posX() > ceh_a.fEventCuts[eVertex_x][eMax])) { - return kFALSE; - } - - // Vertex_y: - if ((collision.posY() < ceh_a.fEventCuts[eVertex_y][eMin]) || - (collision.posY() > ceh_a.fEventCuts[eVertex_y][eMax])) { - return kFALSE; - } - - // Vertex_z: - if ((collision.posZ() < ceh_a.fEventCuts[eVertex_z][eMin]) || - (collision.posZ() > ceh_a.fEventCuts[eVertex_z][eMax])) { - return kFALSE; - } - - // NContributors: - if ((collision.numContrib() < ceh_a.fEventCuts[eNContributors][eMin]) || - (collision.numContrib() > ceh_a.fEventCuts[eNContributors][eMax])) { - return kFALSE; - } - - return kTRUE; - -} // void EventCuts(aod::Collision const& collision, TracksRec const& tracksRec) - -//============================================================ - -template -Bool_t EventCuts(aod::Collision const& collision, T const& tracks) +template +Bool_t EventCuts(T1 const& collision, T2 const& tracks) { // Event cuts on reconstructed and simulated data. @@ -1233,8 +1176,8 @@ Bool_t EventCuts(aod::Collision const& collision, T const& tracks) //============================================================ -template -void FillEventHistograms(aod::Collision const& collision, T const& tracks, eRecSim rs, eBeforeAfter ba) +template +void FillEventHistograms(T1 const& collision, T2 const& tracks, eBeforeAfter ba) { // Fill all event histograms for reconstructed or simulated data. @@ -1247,48 +1190,48 @@ void FillEventHistograms(aod::Collision const& collision, T const& tracks, eRecS } // a) Common to both reconstructed and simulated; - - cout << Form("ceh_a.fEventHistograms[eNumberOfEvents][%d][%d] = ", rs, ba) << ceh_a.fEventHistograms[eNumberOfEvents][rs][ba] << endl; - ceh_a.fEventHistograms[eNumberOfEvents][rs][ba]->Fill(0.5); - switch (rs) { - // b) Only reconstructed: - case eRec: - ceh_a.fEventHistograms[eTotalMultiplicity][eRec][ba]->Fill(tracks.size()); - // SelectedTracks => filled directly in void process( ... ) - // ceh_a.fEventHistograms[eCentrality][eRec][ba]->Fill( TBI ); => TBI - // 20231007 not ready yet - ceh_a.fEventHistograms[eVertex_x][eRec][ba]->Fill(collision.posX()); - ceh_a.fEventHistograms[eVertex_y][eRec][ba]->Fill(collision.posY()); - ceh_a.fEventHistograms[eVertex_z][eRec][ba]->Fill(collision.posZ()); - ceh_a.fEventHistograms[eNContributors][eRec][ba]->Fill(collision.numContrib()); - break; - - // c) Only simulated: - case eSim: - ceh_a.fEventHistograms[eTotalMultiplicity][eSim][ba]->Fill(tracks.size()); - // SelectedTracks => filled directly in void process( ... ) - // ceh_a.fEventHistograms[eCentrality][eSim][ba]->Fill( TBI ); => TBI - // 20231007 not ready yet - /* - ceh_a.fEventHistograms[eVertex_x][eSim][ba]->Fill(collision.posX()); - ceh_a.fEventHistograms[eVertex_y][eSim][ba]->Fill(collision.posY()); - ceh_a.fEventHistograms[eVertex_z][eSim][ba]->Fill(collision.posZ()); - ceh_a.fEventHistograms[eNContributors][eSim][ba]->Fill(collision.numContrib()); - */ - break; - } // switch(rs) { - -} // template void FillEventHistograms(...) + // b) Only reconstructed: + if constexpr (rs == eRec) { + ceh_a.fEventHistograms[eTotalMultiplicity][eRec][ba]->Fill(tracks.size()); + // SelectedTracks => filled directly in void process( ... ) + // ceh_a.fEventHistograms[eCentrality][eRec][ba]->Fill( TBI ); => TBI + // 20231007 not ready yet + ceh_a.fEventHistograms[eVertex_x][eRec][ba]->Fill(collision.posX()); + ceh_a.fEventHistograms[eVertex_y][eRec][ba]->Fill(collision.posY()); + ceh_a.fEventHistograms[eVertex_z][eRec][ba]->Fill(collision.posZ()); + ceh_a.fEventHistograms[eNContributors][eRec][ba]->Fill(collision.numContrib()); + } // if constexpr (rs == eRec) { + + // c) Only simulated: + if constexpr (rs == eSim) { + + // TBI 20231021 here I need to access info from CollisionSim = aod::McCollision; + + // ceh_a.fEventHistograms[eTotalMultiplicity][eSim][ba]->Fill(tracks.size()); + // SelectedTracks => filled directly in void process( ... ) + // ceh_a.fEventHistograms[eCentrality][eSim][ba]->Fill( TBI ); => TBI + // 20231007 not ready yet + /* + ceh_a.fEventHistograms[eVertex_x][eSim][ba]->Fill(collision.posX()); + ceh_a.fEventHistograms[eVertex_y][eSim][ba]->Fill(collision.posY()); + ceh_a.fEventHistograms[eVertex_z][eSim][ba]->Fill(collision.posZ()); + ceh_a.fEventHistograms[eNContributors][eSim][ba]->Fill(collision.numContrib()); + */ + } // if constexpr (rs == eSim) { + +} // template void FillEventHistograms(...) //============================================================ template -Bool_t ParticleCuts(T const& track, eRecSim rs) +Bool_t ParticleCuts(T const& track) { // Particles cuts. + // TBI 20231021 implement remaining cuts + if (fVerboseForEachParticle) { LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); } @@ -1298,12 +1241,13 @@ Bool_t ParticleCuts(T const& track, eRecSim rs) } return kTRUE; -} + +} // template Bool_t ParticleCuts(T const& track) //============================================================ -template -void FillParticleHistograms(T const& track, eRecSim rs, eBeforeAfter ba) +template +void FillParticleHistograms(T const& track, eBeforeAfter ba) { // Fill all particle histograms for reconstructed and simulated data. @@ -1311,28 +1255,17 @@ void FillParticleHistograms(T const& track, eRecSim rs, eBeforeAfter ba) LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); } - // From aod::Tracks - cph_a.fParticleHistograms[ePhi][rs][ba]->Fill(track.phi()); - cph_a.fParticleHistograms[ePt][rs][ba]->Fill(track.pt()); - cph_a.fParticleHistograms[eEta][rs][ba]->Fill(track.eta()); - - switch (rs) { - // a) Only reconstructed: - case eRec: - cph_a.fParticleHistograms[ePhi][rs][ba]->Fill(track.phi()); - cph_a.fParticleHistograms[ePt][rs][ba]->Fill(track.pt()); - cph_a.fParticleHistograms[eEta][rs][ba]->Fill(track.eta()); - break; + if constexpr (rs == eRec) { + cph_a.fParticleHistograms[ePhi][eRec][ba]->Fill(track.phi()); + cph_a.fParticleHistograms[ePt][eRec][ba]->Fill(track.pt()); + cph_a.fParticleHistograms[eEta][eRec][ba]->Fill(track.eta()); + } - // b) Only simulated: - case eSim: - /* TBI 20231019 not ready yet - cph_a.fParticleHistograms[ePhi][rs][ba]->Fill(track.phi()); - cph_a.fParticleHistograms[ePt][rs][ba]->Fill(track.pt()); - cph_a.fParticleHistograms[eEta][rs][ba]->Fill(track.eta()); -*/ - break; - } // switch(rs) { + if constexpr (rs == eSim) { + cph_a.fParticleHistograms[ePhi][eSim][ba]->Fill(track.phi()); + cph_a.fParticleHistograms[ePt][eSim][ba]->Fill(track.pt()); + cph_a.fParticleHistograms[eEta][eSim][ba]->Fill(track.eta()); + } /* TBI 20231019 use also these + check further // From aod::TracksExtra @@ -1341,9 +1274,9 @@ void FillParticleHistograms(T const& track, eRecSim rs, eBeforeAfter ba) // From aod::TracksDCA cph_a.fParticleHistograms[eDCA_xy][rs][ba]->Fill(track.dcaXY()); cph_a.fParticleHistograms[eDCA_z][rs][ba]->Fill(track.dcaZ()); -*/ + */ -} // template void FillParticleHistograms(T const& track, eRecSim rs, eBeforeAfter ba) +} // template void FillParticleHistograms(...) //============================================================ @@ -2294,7 +2227,7 @@ void GetParticleWeights() //============================================================ -Bool_t MaxNumberOfEvents(const Int_t rs) +Bool_t MaxNumberOfEvents() { // Check if max number of events was reached. See also configurable cNumberOfEvents_max. @@ -2302,23 +2235,35 @@ Bool_t MaxNumberOfEvents(const Int_t rs) LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); } + // *) Return value: Bool_t reachedMaxNumberOfEvents = kFALSE; + // *) Determine from which histogram the relevant info will be taken: + Int_t rs = -44; // reconstructed or simulated + if (gProcessRec || gProcessRecSim) { + rs = eRec; + } else if (gProcessSim) { + rs = eSim; + } else { + LOGF(fatal, "in function \033[1;31m%s at line %d, not a single flag gProcess* is true \033[0m", __PRETTY_FUNCTION__, __LINE__); + } + + // *) Okay, do the thing: if (ceh_a.fEventHistograms[eNumberOfEvents][rs][eAfter]->GetBinContent(1) >= ceh_a.fEventCuts[eNumberOfEvents][eMax]) { reachedMaxNumberOfEvents = kTRUE; } + // *) Hasta la vista: return reachedMaxNumberOfEvents; } // void MaxNumberOfEvents() //============================================================ -template -void MainLoopOverParticles(T const& tracks, eRecSim rs) +void MainLoopOverParticlesRec(TracksRec const& tracks) { // This is the main loop over particles, in which Q-vectors and particle histograms are filled, particle cuts applied, etc. - // If I want to access both rec and sim info in this loop, set gProcessRecSim = true via configurable "cfWhatToProcess" + // Only 'rec' info is processed in this loop, set gProcessRec = true via configurable "cfWhatToProcess" if (fVerbose) { LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); @@ -2332,49 +2277,129 @@ void MainLoopOverParticles(T const& tracks, eRecSim rs) for (auto& track : tracks) { // *) Fill particle histograms for reconstructed particle before particle cuts: - FillParticleHistograms(track, rs, eBefore); + FillParticleHistograms(track, eBefore); - // *) Fill particle histograms for corresponding simulated particle before particle cuts: - if (gProcessRecSim) { - if (track.has_mcParticle()) { - auto mcParticle = track.mcParticle(); - FillParticleHistograms(mcParticle, eSim, eBefore); // great, this works straight like this, thanks O2 people! - } else { - LOGF(debug, "in function \033[1;31m%s at line %d, track.has_mcParticle() is false. \033[0m", __PRETTY_FUNCTION__, __LINE__); + // *) Particle cuts: + if (!ParticleCuts(track)) { + continue; + } + + // *) Fill particle histograms for reconstructed data after particle cuts: + FillParticleHistograms(track, eAfter); + + // *) Fill Q-vectors: + // Take kinematics from reconstructed particles: + dPhi = track.phi(); + dPt = track.pt(); + dEta = track.eta(); + + // Particle weights: + if (pw_a.fUseWeights[wPHI]) { + wPhi = Weight(dPhi, "phi"); // corresponding phi weight + if (!(wPhi > 0.)) { + LOGF(error, "\033[1;33m%s wPhi is not positive, skipping this particle for the time being...\033[0m", __PRETTY_FUNCTION__); + LOGF(error, "dPhi = %f\nwPhi = %f", dPhi, wPhi); + continue; + } + } // if(pw_a.fUseWeights[wPHI]) + if (pw_a.fUseWeights[wPT]) { + wPt = Weight(dPt, "pt"); // corresponding pt weight + if (!(wPt > 0.)) { + LOGF(error, "\033[1;33m%s wPt is not positive, skipping this particle for the time being...\033[0m", __PRETTY_FUNCTION__); + LOGF(error, "dPt = %f\nwPt = %f", dPt, wPt); + continue; + } + } // if(pw_a.fUseWeights[wPT]) + if (pw_a.fUseWeights[wETA]) { + wEta = Weight(dEta, "eta"); // corresponding eta weight + if (!(wEta > 0.)) { + LOGF(error, "\033[1;33m%s wEta is not positive, skipping this particle for the time being...\033[0m", __PRETTY_FUNCTION__); + LOGF(error, "dEta = %f\nwEta = %f", dEta, wEta); + continue; } - } // if (gProcessRecSim) { + } // if(pw_a.fUseWeights[wETA]) + + for (Int_t h = 0; h < gMaxHarmonic * gMaxCorrelator + 1; h++) { + for (Int_t wp = 0; wp < gMaxCorrelator + 1; wp++) { // weight power + if (pw_a.fUseWeights[wPHI] || pw_a.fUseWeights[wPT] || pw_a.fUseWeights[wETA]) { + wToPowerP = pow(wPhi * wPt * wEta, wp); + } + qv_a.fQvector[h][wp] += TComplex(wToPowerP * TMath::Cos(h * dPhi), wToPowerP * TMath::Sin(h * dPhi)); + } // for(Int_t wp=0;wpAddAt(dPhi, fSelectedTracks); + } // remember that the 2nd argument here must start from 0 + if (nl_a.ftaNestedLoops[1]) { + nl_a.ftaNestedLoops[1]->AddAt(wPhi * wPt * wEta, fSelectedTracks); + } // remember that the 2nd argument here must start from 0 + } // if(fCalculateNestedLoops||fCalculateCustomNestedLoop) + + // *) Counter of selected tracks in the current event: + fSelectedTracks++; + if (fSelectedTracks >= cSelectedTracks_max) { + break; + } + + } // for (auto& track : tracks) + +} // void MainLoopOverParticlesRec(TracksRec const& tracks) + +//============================================================ + +void MainLoopOverParticlesRecSim(TracksRecSim const& tracks) +{ + // This is the main loop over particles, in which Q-vectors and particle histograms are filled, particle cuts applied, etc. + // Both 'rec' and 'sim' info are processed in this loop, set gProcessRecSim = true via configurable "cfWhatToProcess" + // It sufficcess to pass only rec info, because here I access sim info via labels. + + if (fVerbose) { + LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); + } + + Double_t dPhi = 0., wPhi = 1.; // azimuthal angle and corresponding phi weight + Double_t dPt = 0., wPt = 1.; // transverse momentum and corresponding pT weight + Double_t dEta = 0., wEta = 1.; // pseudorapidity and corresponding eta weight + Double_t wToPowerP = 1.; // weight raised to power p + fSelectedTracks = 0; // reset number of selected tracks + for (auto& track : tracks) { + + // *) Fill particle histograms for reconstructed particle before particle cuts: + FillParticleHistograms(track, eBefore); + + // *) Fill particle histograms for corresponding simulated particle before particle cuts: + if (track.has_mcParticle()) { + auto mcParticle = track.mcParticle(); + FillParticleHistograms(mcParticle, eBefore); // great, this works straight like this, thanks O2 people! + } else { + LOGF(debug, "in function \033[1;31m%s at line %d, track.has_mcParticle() is false. \033[0m", __PRETTY_FUNCTION__, __LINE__); + } // *) Particle cuts: - if (!ParticleCuts(track, rs)) { + if (!ParticleCuts(track)) { continue; } // *) Fill particle histograms for reconstructed data after particle cuts: - FillParticleHistograms(track, rs, eAfter); + FillParticleHistograms(track, eAfter); // *) Fill particle histograms for corresponding simulated particle after particle cuts: // TBI 20231020 there is a duplication of code here, see above. I should get 'mcParticle' only once for this 'track' in this loop - if (gProcessRecSim) { - if (track.has_mcParticle()) { - auto mcParticle = track.mcParticle(); - FillParticleHistograms(mcParticle, eSim, eAfter); // great, this works straight like this, thanks O2 people! - } else { - LOGF(debug, "in function \033[1;31m%s at line %d, track.has_mcParticle() is false. \033[0m", __PRETTY_FUNCTION__, __LINE__); - } - } // if (gProcessRecSim) { + if (track.has_mcParticle()) { + auto mcParticle = track.mcParticle(); + FillParticleHistograms(mcParticle, eAfter); // great, this works straight like this, thanks O2 people! + } else { + LOGF(debug, "in function \033[1;31m%s at line %d, track.has_mcParticle() is false. \033[0m", __PRETTY_FUNCTION__, __LINE__); + } // *) Fill Q-vectors: - if (gProcessRec || gProcessRecSim) { - // Take kinematics from reconstructed particles: - dPhi = track.phi(); - dPt = track.pt(); - dEta = track.eta(); - } else if (gProcessSim) { - // Take kinematics from simulated particles: - dPhi = track.mcParticle().phi(); - dPt = track.mcParticle().pt(); - dEta = track.mcParticle().eta(); - } + // Take kinematics from reconstructed particles: + dPhi = track.phi(); + dPt = track.pt(); + dEta = track.eta(); // Particle weights: if (pw_a.fUseWeights[wPHI]) { @@ -2429,7 +2454,22 @@ void MainLoopOverParticles(T const& tracks, eRecSim rs) } // for (auto& track : tracks) -} // template void MainLoopOverParticles(T const& tracks) +} // void MainLoopOverParticlesRecSim(TracksRecSim const& tracks) + +//============================================================ + +void MainLoopOverParticlesSim(TracksSim const& tracks) +{ + // This is the main loop over particles, in which Q-vectors and particle histograms are filled, particle cuts applied, etc. + // Only 'sim' info is processed in this loop, set gProcessSim = true via configurable "cfWhatToProcess" + + if (fVerbose) { + LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); + } + + // TBI 20231021 tbc... + +} // void MainLoopOverParticlesSim(TracksSim const& tracks) //============================================================ @@ -2459,4 +2499,104 @@ void CalculateEverything() //============================================================ +void SteerRec(CollisionRec const& collision, TracksRec const& tracks) +{ + // This is the only function to be called in processRec(...), processRecSim(...), and processSim(...). + // All analysis workflow is defined step-by-step here, via dedicated function calls. + // Order of function calls obviously does matter. + + if (fVerbose) { + LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); + } + + // *) Do all thingies before starting to process data from this collision (e.g. count number of events, fetch the run number, etc.): + Preprocess(collision); // see enum eRecSim for possible values of 'rsb' + + // *) Fill event histograms for reconstructed data before event cuts: + FillEventHistograms(collision, tracks, eBefore); + + // *) Event cuts: + if (!EventCuts(collision, tracks)) { + return; + } + + // *) Fill event histograms for reconstructed data after event cuts: + FillEventHistograms(collision, tracks, eAfter); + + // *) Main loop over particles: + MainLoopOverParticlesRec(tracks); + + // *) Fill remaining event histograms for reconstructed data after event AND after particle cuts: + ceh_a.fEventHistograms[eSelectedTracks][eRec][eAfter]->Fill(fSelectedTracks); + + // *) Remaining event cuts which can be applied only after the loop over particles is performed: + if ((fSelectedTracks < ceh_a.fEventCuts[eSelectedTracks][eMin]) || (fSelectedTracks > ceh_a.fEventCuts[eSelectedTracks][eMax])) { + return; + } + + // *) Calculate everything for selected events and particles: + CalculateEverything(); + + // *) Reset event-by-event quantities: + ResetEventByEventQuantities(); + +} // void SteerRec(CollisionRec const& collision, TracksRec const& tracks) + +//============================================================ + +void SteerRecSim(CollisionRec const& collision, TracksRecSim const& tracks) +{ + // This is the only function to be called in processRec(...), processRecSim(...), and processSim(...). + // All analysis workflow is defined step-by-step here, via dedicated function calls. + // Order of function calls obviously does matter. + + if (fVerbose) { + LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__); + } + + // *) Do all thingies before starting to process data from this collision (e.g. count number of events, fetch the run number, etc.): + Preprocess(collision); // see enum eRecSim for possible values of 'rsb' + + // *) Fill event histograms for reconstructed data before event cuts: + FillEventHistograms(collision, tracks, eBefore); + FillEventHistograms(collision, tracks, eBefore); + + // *) Event cuts: + if (!EventCuts(collision, tracks)) { + return; + } + + // *) Fill event histograms for reconstructed data after event cuts: + FillEventHistograms(collision, tracks, eAfter); + FillEventHistograms(collision, tracks, eAfter); + + // *) Main loop over particles: + MainLoopOverParticlesRecSim(tracks); + + // *) Fill remaining event histograms for reconstructed data after event AND after particle cuts: + ceh_a.fEventHistograms[eSelectedTracks][eRec][eAfter]->Fill(fSelectedTracks); + // ceh_a.fEventHistograms[eSelectedTracks][eSim][eAfter]->Fill(fSelectedTracks); // TBI 20231020 do I really need this one? + + // *) Remaining event cuts which can be applied only after the loop over particles is performed: + if ((fSelectedTracks < ceh_a.fEventCuts[eSelectedTracks][eMin]) || (fSelectedTracks > ceh_a.fEventCuts[eSelectedTracks][eMax])) { + return; + } + + // *) Calculate everything for selected events and particles: + CalculateEverything(); + + // *) Reset event-by-event quantities: + ResetEventByEventQuantities(); + +} // void SteerRecSim(CollisionRec const& collision, TracksRecSim const& tracks) + +//============================================================ + +void SteerSim(CollisionSim const& collision, TracksSim const& tracks) +{ + + // TBI 20231021 tbc... + +} // void SteerSim(CollisionSim const& collision, TracksSim const& tracks) + #endif // PWGCF_MULTIPARTICLECORRELATIONS_CORE_MUPA_MEMBERFUNCTIONS_H_ diff --git a/PWGCF/MultiparticleCorrelations/Tasks/multiparticle-correlations-ab.cxx b/PWGCF/MultiparticleCorrelations/Tasks/multiparticle-correlations-ab.cxx index 55e536e2d6f..27c5a683c10 100644 --- a/PWGCF/MultiparticleCorrelations/Tasks/multiparticle-correlations-ab.cxx +++ b/PWGCF/MultiparticleCorrelations/Tasks/multiparticle-correlations-ab.cxx @@ -18,14 +18,17 @@ using namespace o2; using namespace o2::framework; +using CollisionRec = aod::Collision; +using CollisionSim = aod::McCollision; + using TracksRec = soa::Join; using TrackRec = soa::Join::iterator; using TracksRecSim = soa::Join; using TrackRecSim = soa::Join::iterator; -using TracksSim = soa::Join; // TBI 20231018 only temporarily defined this way, check still aod::McParticles -using TrackSim = soa::Join::iterator; // TBI 20231018 only temporarily defined this way check still aod::McParticles +using TracksSim = aod::McParticles; +using TrackSim = aod::McParticles::iterator; // ROOT: #include "TList.h" @@ -122,53 +125,29 @@ struct MultiparticleCorrelationsAB // this name is used in lower-case format to // A) Process only reconstructed data; // B) Process both reconstructed and simulated data; // C) Process only simulated data. - // TBI 20231017 There is a (minor) code bloat here obviously. Fix this when it becomes possible to post the last argument to PROCESS_SWITCH via variable (ideally Configurable) + + // ------------------------------------------- // A) Process only reconstructed data: - void processRec(aod::Collision const& collision, aod::BCs const&, TracksRec const& tracksRec) // called once per collision found in the time frame + void processRec(CollisionRec const& collision, aod::BCs const&, TracksRec const& tracks) { // ... - // *) If I reached max number of events, ignore the remaining collisions: - if (!fProcessRemainingEvents) { - return; // TBI 20231008 Temporarily implemented this way. But what I really need here is a graceful exit from subsequent processing (which will also dump the output file, etc.) - } - // *) TBI 20231020 Temporary here (use configurable 'cfWhatToProcess' to set this flag corerectly): if (!gProcessRec) { LOGF(fatal, "in function \033[1;31m%s at line %d\033[0m", __PRETTY_FUNCTION__, __LINE__); } - // *) Do all thingies before starting to process data (e.g. count number of events, fetch the run number, etc.): - Preprocess(collision, eRec); - - // *) Fill event histograms for reconstructed data before event cuts: - FillEventHistograms(collision, tracksRec, eRec, eBefore); - - // *) Event cuts: - if (!EventCuts(collision, tracksRec)) { - return; - } - - // *) Fill event histograms for reconstructed data after event cuts: - FillEventHistograms(collision, tracksRec, eRec, eAfter); - - // *) Main loop over reconstructed particles: - // MainLoopOverParticles(tracksRec, eRec); TBI 20231020 yes, this is a bug, becase I call track.has_mcParticle() there, etc. Re-think how to re-implement this - - // *) Fill remaining event histograms for reconstructed data after event AND after particle cuts: - ceh_a.fEventHistograms[eSelectedTracks][eRec][eAfter]->Fill(fSelectedTracks); - - // *) Remaining event cuts which can be applied only after the loop over particles is performed: - if ((fSelectedTracks < ceh_a.fEventCuts[eSelectedTracks][eMin]) || (fSelectedTracks > ceh_a.fEventCuts[eSelectedTracks][eMax])) { - return; + // *) If I reached max number of events, ignore the remaining collisions: + if (!fProcessRemainingEvents) { + return; // TBI 20231008 Temporarily implemented this way. But what I really need here is a graceful exit + // from subsequent processing (which will also dump the output file, etc.). When that's possible, + // move this to a member function Steer*(...) } - // *) Calculate everything for selected events and particles: - CalculateEverything(); - - // *) Reset event-by-event quantities: - ResetEventByEventQuantities(); + // *) Steer all analysis steps: + SteerRec(collision, tracks); + // TBI 20231021 If I want to do some postprocessing after Steer(...), re-define Steer from void to bool, so that event cuts have effect, etc, } // void processRec(...) PROCESS_SWITCH(MultiparticleCorrelationsAB, processRec, "process only reconstructed information", false); @@ -176,54 +155,25 @@ struct MultiparticleCorrelationsAB // this name is used in lower-case format to // ------------------------------------------- // B) Process both reconstructed and simulated data: - void processRecSim(aod::Collision const& collision, aod::BCs const&, TracksRecSim const& tracksRecSim, aod::McParticles const&) // called once per collision found in the time frame + void processRecSim(CollisionRec const& collision, aod::BCs const&, TracksRecSim const& tracks, aod::McParticles const&) { - // Remark: Implemented separately to processRec(...), because here I need to subscribe also to Monte Carlo tables. - - // *) If I reached max number of events, ignore the remaining collisions: - if (!fProcessRemainingEvents) { - return; // TBI 20231008 Temporarily implemented this way. But what I really need here is a graceful exit from subsequent processing (which will also dump the output file, etc.) - } + // ... // *) TBI 20231020 Temporary here (use configurable 'cfWhatToProcess' to set this flag corerectly): if (!gProcessRecSim) { LOGF(fatal, "in function \033[1;31m%s at line %d\033[0m", __PRETTY_FUNCTION__, __LINE__); } - // *) Do all thingies before starting to process data (e.g. count number of events, fetch the run number, etc.): - Preprocess(collision, eRec); - - // *) Fill event histograms for reconstructed and simulated data before event cuts: - FillEventHistograms(collision, tracksRecSim, eRec, eBefore); - FillEventHistograms(collision, tracksRecSim, eSim, eBefore); - - // *) Event cuts: - if (!EventCuts(collision, tracksRecSim)) { - return; - } - - // *) Fill event histograms for reconstructed and simulated data after event cuts: - FillEventHistograms(collision, tracksRecSim, eRec, eAfter); - FillEventHistograms(collision, tracksRecSim, eSim, eAfter); - - // *) Main loop over reconstructed particles: - // For simulated particles there is NO separate loop, if flag gProcessRecSim = true (via configurable "cfWhatToProcess"), MC info is accessed via label. - MainLoopOverParticles(tracksRecSim, eRec); - - // *) Fill remaining event histograms for reconstructed and simulated data after event AND after particle cuts: - ceh_a.fEventHistograms[eSelectedTracks][eRec][eAfter]->Fill(fSelectedTracks); - // ceh_a.fEventHistograms[eSelectedTracks][eSim][eAfter]->Fill(fSelectedTracks); // TBI 20231020 do I really need this one? - - // *) Remaining event cuts which can be applied only after the loop over particles is performed: - if ((fSelectedTracks < ceh_a.fEventCuts[eSelectedTracks][eMin]) || (fSelectedTracks > ceh_a.fEventCuts[eSelectedTracks][eMax])) { - return; + // *) If I reached max number of events, ignore the remaining collisions: + if (!fProcessRemainingEvents) { + return; // TBI 20231008 Temporarily implemented this way. But what I really need here is a graceful exit + // from subsequent processing (which will also dump the output file, etc.). When that's possible, + // move this to a member function Steer(...) } - // *) Calculate everything for selected events and particles: - CalculateEverything(); - - // *) Reset event-by-event quantities: - ResetEventByEventQuantities(); + // *) Steer all analysis steps: + SteerRecSim(collision, tracks); + // TBI 20231021 If I want to do some postprocessing after Steer(...), re-define Steer from void to bool, so that event cuts have effect, etc, } // void processRecSim(...) PROCESS_SWITCH(MultiparticleCorrelationsAB, processRecSim, "process both reconstructed and simulated information", true); @@ -231,7 +181,7 @@ struct MultiparticleCorrelationsAB // this name is used in lower-case format to // ------------------------------------------- // C) Process only simulated data: - void processSim(aod::Collision const& collision, aod::BCs const&, TracksSim const& tracksSim, aod::McParticles const&) // called once per collision found in the time frame + void processSim(CollisionSim const& collision, aod::BCs const&, TracksSim const& tracks) { // ... @@ -245,37 +195,9 @@ struct MultiparticleCorrelationsAB // this name is used in lower-case format to LOGF(fatal, "in function \033[1;31m%s at line %d\033[0m", __PRETTY_FUNCTION__, __LINE__); } - // *) Do all thingies before starting to process data (e.g. count number of events, fetch the run number, etc.): - Preprocess(collision, eSim); - - // *) Fill event histograms for simulated data before event cuts: - FillEventHistograms(collision, tracksSim, eRec, eBefore); - - // *) Event cuts: - if (!EventCuts(collision, tracksSim)) { - return; - } - - // *) Fill event histograms for simulated data after event cuts: - FillEventHistograms(collision, tracksSim, eSim, eAfter); - - // *) Main loop over reconstructed particles: - // For simulated particles there is NO separate loop, if flag gProcessRecSim = true (via configurable "cfWhatToProcess"), MC info is accessed via label. - // MainLoopOverParticles(tracksSim, eSim); // TBI 20231020 it will NOT work this way + re-think the definition of tracksSim in the preamble - - // *) Fill remaining event histograms for reconstructed and simulated data after event AND after particle cuts: - ceh_a.fEventHistograms[eSelectedTracks][eSim][eAfter]->Fill(fSelectedTracks); - - // *) Remaining event cuts which can be applied only after the loop over particles is performed: - if ((fSelectedTracks < ceh_a.fEventCuts[eSelectedTracks][eMin]) || (fSelectedTracks > ceh_a.fEventCuts[eSelectedTracks][eMax])) { - return; - } - - // *) Calculate everything for selected events and particles: - CalculateEverything(); - - // *) Reset event-by-event quantities: - ResetEventByEventQuantities(); + // *) Steer all analysis steps: + SteerSim(collision, tracks); + // TBI 20231021 If I want to do some postprocessing after Steer(...), re-define Steer from void to bool, so that event cuts have effect, etc, } // void processSim(...) PROCESS_SWITCH(MultiparticleCorrelationsAB, processSim, "process only simulated information", false); From 8b83af4e990044244f7889b6340cd87caee74e6d Mon Sep 17 00:00:00 2001 From: abilandz Date: Sat, 21 Oct 2023 13:26:48 +0200 Subject: [PATCH 4/4] fix for formatter #2 --- PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h b/PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h index fafe9e0cb24..0d3e899509c 100644 --- a/PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h +++ b/PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h @@ -1213,7 +1213,7 @@ void FillEventHistograms(T1 const& collision, T2 const& tracks, eBeforeAfter ba) // SelectedTracks => filled directly in void process( ... ) // ceh_a.fEventHistograms[eCentrality][eSim][ba]->Fill( TBI ); => TBI // 20231007 not ready yet - /* + /* ceh_a.fEventHistograms[eVertex_x][eSim][ba]->Fill(collision.posX()); ceh_a.fEventHistograms[eVertex_y][eSim][ba]->Fill(collision.posY()); ceh_a.fEventHistograms[eVertex_z][eSim][ba]->Fill(collision.posZ());