Skip to content

support for randomization via Fisher-Yates algorithm #4360

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
victor-gonzalez merged 1 commit into from
Jan 18, 2024
Merged

support for randomization via Fisher-Yates algorithm #4360

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 4 additions & 3 deletions PWGCF/MultiparticleCorrelations/Core/MuPa-DataMembers.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@ struct TaskConfiguration {
TString fWhatToProcess = "Rec"; // "Rec" = process only reconstructed, "Sim" = process only simulated, "RecSim" = process both reconstructed and simulated
UInt_t fRandomSeed = 0; // argument to TRandom3 constructor. By default it is 0 (i.e. seed is guaranteed to be unique in time and space), use SetRandomSeed(...) to change it
Bool_t fUseFisherYates = kFALSE; // algorithm used to randomize particle indices, set via configurable
TArrayI* fRandomIndices = NULL; // array to store random indices obtained from Fisher-Yates algorithm
Int_t fFixedNumberOfRandomlySelectedTracks = -1; // use a fixed number of randomly selected particles in each event. It is set and applied, if > 0. Set to <=0 to ignore.

// Bool_t fRescaleWithTheoreticalInput; // if kTRUE, all measured correlators are
Expand Down Expand Up @@ -126,7 +127,7 @@ Bool_t fCalculateCorrelations =
struct Correlations_Arrays {
TProfile* fCorrelationsPro[4][gMaxHarmonic][eAsFunctionOf_N] = {
{{NULL}}}; //! multiparticle correlations
//! [2p=0,4p=1,6p=2,8p=3][n=1,n=2,...,n=6][0=integrated,1=vs.
//! [2p=0,4p=1,6p=2,8p=3][n=1,n=2,...,n=gMaxHarmonic][0=integrated,1=vs.
//! multiplicity,2=vs. centrality,3=pT,4=eta]
} c_a;

Expand All @@ -151,9 +152,9 @@ Bool_t fCalculateNestedLoops = kTRUE; // calculate and store correlations with
Bool_t fCalculateCustomNestedLoop =
kFALSE; // validate e-b-e all correlations with custom nested loop
struct NestedLoops_Arrays {
TProfile* fNestedLoopsPro[4][6][eAsFunctionOf_N] = {
TProfile* fNestedLoopsPro[4][gMaxHarmonic][eAsFunctionOf_N] = {
{{NULL}}}; //! multiparticle correlations from nested loops
//! [2p=0,4p=1,6p=2,8p=3][n=1,n=2,...,n=6][0=integrated,1=vs.
//! [2p=0,4p=1,6p=2,8p=3][n=1,n=2,...,n=gMaxHarmonic][0=integrated,1=vs.
//! multiplicity,2=vs. centrality,3=pT,4=eta]
TArrayD* ftaNestedLoops[2] = {NULL}; //! e-b-e container for nested loops
//! [0=angles;1=product of all weights]
Expand Down
2 changes: 1 addition & 1 deletion PWGCF/MultiparticleCorrelations/Core/MuPa-GlobalConstants.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
#define PWGCF_MULTIPARTICLECORRELATIONS_CORE_MUPA_GLOBALCONSTANTS_H_

const Int_t gMaxCorrelator = 12;
const Int_t gMaxHarmonic = 6;
const Int_t gMaxHarmonic = 9;
const Int_t gMaxIndex = 300; // per order, used only in Test0

#endif // PWGCF_MULTIPARTICLECORRELATIONS_CORE_MUPA_GLOBALCONSTANTS_H_
92 changes: 83 additions & 9 deletions PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ void BookBaseList()
Form("fTaskName = %s", tc.fTaskName.Data()));

fBasePro->GetXaxis()->SetBinLabel(eRunNumber,
Form("tc.fRunNumber = %s", tc.fRunNumber.Data()));
Form("fRunNumber = %s", tc.fRunNumber.Data()));

fBasePro->GetXaxis()->SetBinLabel(eVerbose, "fVerbose");
fBasePro->Fill(eVerbose - 0.5, (Int_t)tc.fVerbose);
Expand All @@ -57,7 +57,7 @@ void BookBaseList()
fBasePro->Fill(eProcessRemainingEvents - 0.5, (Int_t)tc.fProcessRemainingEvents);

fBasePro->GetXaxis()->SetBinLabel(eWhatToProcess,
Form("tc.WhatToProcess = %s", tc.fWhatToProcess.Data()));
Form("WhatToProcess = %s", tc.fWhatToProcess.Data()));

fBasePro->GetXaxis()->SetBinLabel(eRandomSeed, "fRandomSeed");
fBasePro->Fill(eRandomSeed - 0.5, (Int_t)tc.fRandomSeed);
Expand Down Expand Up @@ -126,6 +126,10 @@ void DefaultConfiguration()
// identical to the internal definitions in MuPa-Configurables.h, the
// settings in json file are silently ignored.

// c) Scientific notation is NOT supported in json file. E.g. if you have
// "cSelectedTracks_max": "1e3",
// that setting and ALL other ones in json are silently ignored.

if (tc.fVerbose) {
LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__);
}
Expand Down Expand Up @@ -1010,7 +1014,7 @@ void BookTest0Histograms()

// a) Book the profile holding flags;
// b) Book placeholder and make sure all labels are stored in the placeholder;
// c) Retreive labels from placeholder;
// c) Retrieve labels from placeholder;
// d) Book what needs to be booked;
// e) Few quick insanity checks on booking.

Expand All @@ -1036,7 +1040,7 @@ void BookTest0Histograms()
fTest0List->Add(fTest0LabelsPlaceholder);
}

// c) Retreive labels from placeholder:
// c) Retrieve labels from placeholder:
if (!(this->RetrieveCorrelationsLabels())) {
LOGF(fatal, "in function \033[1;31m%s at line %d\033[0m",
__PRETTY_FUNCTION__, __LINE__);
Expand Down Expand Up @@ -1221,7 +1225,7 @@ void ResetEventByEventQuantities()
fSelectedTracks = 0;
fCentrality = 0;

// c) Q-vectors:
// b) Q-vectors:
if (fCalculateQvector) {
ResetQ(); // generic Q-vector
for (Int_t h = 0; h < gMaxHarmonic * gMaxCorrelator + 1; h++) {
Expand All @@ -1232,7 +1236,7 @@ void ResetEventByEventQuantities()
}
} // if(fCalculateQvector)

// d) Reset ebe containers for nested loops:
// c) Reset ebe containers for nested loops:
if (fCalculateNestedLoops || fCalculateCustomNestedLoop) {
if (nl_a.ftaNestedLoops[0]) {
nl_a.ftaNestedLoops[0]->Reset();
Expand All @@ -1246,7 +1250,13 @@ void ResetEventByEventQuantities()

} // if(fCalculateNestedLoops||fCalculateCustomNestedLoop)

// ... TBI 20220809 port the rest ...
// d) Fisher-Yates algorithm:
if (tc.fUseFisherYates) {
delete tc.fRandomIndices;
tc.fRandomIndices = NULL;
}

// ... TBI 20240117 port the rest ...

} // void ResetEventByEventQuantities()

Expand Down Expand Up @@ -2930,7 +2940,8 @@ void StoreLabelsInPlaceholder()
// a) Initialize all counters;
// b) Fetch TObjArray with labels from an external file;
// c) Book the placeholder fTest0LabelsPlaceholder for all labels;
// d) Finally, store the labels from external source into placeholder.
// d) Finally, store the labels from external source into placeholder;
// e) Insantity check on labels.

if (tc.fVerbose) {
LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__);
Expand Down Expand Up @@ -2985,6 +2996,19 @@ void StoreLabelsInPlaceholder()
// cout<<oa->GetEntries()<<endl;
} // for(Int_t e=0; e<nLabels; e++)

// e) Insantity check on labels:
// Here I am merely checking that harmonic larget than gMaxHarmonic was not requested.
for (Int_t b = 1; b <= fTest0LabelsPlaceholder->GetXaxis()->GetNbins(); b++) {
TObjArray* temp = TString(fTest0LabelsPlaceholder->GetXaxis()->GetBinLabel(b)).Tokenize(" ");
for (Int_t h = 0; h < temp->GetEntries(); h++) {
if (TMath::Abs(TString(temp->At(h)->GetName()).Atoi()) > gMaxHarmonic) {
LOGF(info, "\033[1;31m bin = %d, label = %s, gMaxHarmonic = %d\033[0m", b, fTest0LabelsPlaceholder->GetXaxis()->GetBinLabel(b), (Int_t)gMaxHarmonic);
LOGF(fatal, "in function \033[1;31m%s at line %d\033[0m", __PRETTY_FUNCTION__, __LINE__);
} // if(TString(temp->At(h)->GetName()).Atoi() > gMaxHarmonic) {
} // for(Int_t h = 0; h < temp->GetEntries(); h++) {
delete temp; // yes, otherwise it's a memory leak
} // for(Int_t b = 1; b <= fTest0LabelsPlaceholder->GetXaxis()->GetNbins(); b++) {

} // void StoreLabelsInPlaceholder()

//============================================================
Expand Down Expand Up @@ -3426,6 +3450,35 @@ void DetermineCentrality()

//============================================================

void RandomIndices(Int_t nTracks)
{
// Randomize indices using Fisher-Yates algorithm.

if (tc.fVerbose) {
LOGF(info, "\033[1;32m%s\033[0m", __PRETTY_FUNCTION__);
}

if (nTracks < 1) {
return;
}

// Fisher-Yates algorithm:
tc.fRandomIndices = new TArrayI(nTracks);
tc.fRandomIndices->Reset(); // just in case there is some random garbage in memory at init
for (Int_t i = 0; i < nTracks; i++) {
tc.fRandomIndices->AddAt(i, i);
}
for (Int_t i = nTracks - 1; i >= 1; i--) {
Int_t j = gRandom->Integer(i + 1);
Int_t temp = tc.fRandomIndices->GetAt(j);
tc.fRandomIndices->AddAt(tc.fRandomIndices->GetAt(i), j);
tc.fRandomIndices->AddAt(temp, i);
} // end of for(Int_t i=nTracks-1;i>=1;i--)

} // void RandomIndices(Int_t nTracks)

//============================================================

void CalculateEverything()
{
// Calculate everything for selected events and particles.
Expand Down Expand Up @@ -3531,7 +3584,28 @@ void MainLoopOverParticles(T const& tracks)
Double_t wToPowerP = 1.; // weight raised to power p
fSelectedTracks = 0; // reset number of selected tracks

for (auto& track : tracks) {
// *) If random access of tracks from collection is requested, use Fisher-Yates algorithm to generate random indices:
if (tc.fUseFisherYates) {
if (tc.fRandomIndices) {
LOGF(fatal, "in function \033[1;31m%s at line %d\033[0m", __PRETTY_FUNCTION__, __LINE__);
}
this->RandomIndices(tracks.size());

@victor-gonzalez victor-gonzalez Jan 18, 2024

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is redundant, right?

@abilandz abilandz Jan 18, 2024

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, but it's a matter of taste - I like to use it to indicate that I am calling a member function directly implemented in a given class, and not inherited. It's not a very strong argument, I know...

if (!tc.fRandomIndices) {
LOGF(fatal, "in function \033[1;31m%s at line %d\033[0m", __PRETTY_FUNCTION__, __LINE__);
}
}

// *) Main loop over particles:
// for (auto& track : tracks) { // default standard way of looping of tracks
auto track = tracks.iteratorAt(0); // set the type and scope from one instance
for (int64_t i = 0; i < tracks.size(); i++) {

// *) Access track sequentially from collection of tracks (default), or randomly using Fisher-Yates algorithm:
if (!tc.fUseFisherYates) {
track = tracks.iteratorAt(i);
} else {
track = tracks.iteratorAt((int64_t)tc.fRandomIndices->GetAt(i));
}

// *) Fill particle histograms before particle cuts:
FillParticleHistograms<rs>(track, eBefore);
Expand Down