diff --git a/PWGCF/MultiparticleCorrelations/Core/MuPa-DataMembers.h b/PWGCF/MultiparticleCorrelations/Core/MuPa-DataMembers.h index 1114066d5d3..3f97f56e7c5 100644 --- a/PWGCF/MultiparticleCorrelations/Core/MuPa-DataMembers.h +++ b/PWGCF/MultiparticleCorrelations/Core/MuPa-DataMembers.h @@ -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 @@ -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; @@ -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] diff --git a/PWGCF/MultiparticleCorrelations/Core/MuPa-GlobalConstants.h b/PWGCF/MultiparticleCorrelations/Core/MuPa-GlobalConstants.h index 714ff2c5971..d400ce14185 100644 --- a/PWGCF/MultiparticleCorrelations/Core/MuPa-GlobalConstants.h +++ b/PWGCF/MultiparticleCorrelations/Core/MuPa-GlobalConstants.h @@ -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_ diff --git a/PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h b/PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h index eb0cbd5fca2..4840cff7bc2 100644 --- a/PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h +++ b/PWGCF/MultiparticleCorrelations/Core/MuPa-MemberFunctions.h @@ -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); @@ -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); @@ -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__); } @@ -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. @@ -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__); @@ -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++) { @@ -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(); @@ -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() @@ -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__); @@ -2985,6 +2996,19 @@ void StoreLabelsInPlaceholder() // cout<GetEntries()<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() //============================================================ @@ -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. @@ -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()); + 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(track, eBefore);