add nested loops

Author: Pengchong Hu

PWGCF/MultiparticleCorrelations/Tasks/multiharmonicCorrelations.cxx

@@ -116,6 +116,7 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to
   Configurable<bool> cfInitsim{"cfInitsim", false, "init histograms of sim"};
   Configurable<bool> cfUseWeights{"cfUseWeights", true, "use weights"};
   Configurable<bool> cfToyModel{"cfToyModel", true, "phi-distribution from toy model"};
+  Configurable<bool> cfNest{"cfNest", true, "nested loops"};
 
   Configurable<std::vector<float>> cfVertexZ{"cfVertexZ", {-10, 10.}, "vertex z position range: {min, max}[cm], with convention: min <= Vz < max"};
   Configurable<std::vector<float>> cfPt{"cfPt", {0.2, 5.0}, "transverse momentum range"};
@@ -171,6 +172,12 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to
     TProfile* pfour32_centr = NULL;
     TProfile* pfour42_centr = NULL;
     TComplex Qvector[maxHarmonic][maxPower];
+    std::vector<float> vecphi;
+    std::vector<float> vecwei;
+    TProfile* nestedLoops[maxHarmonic] = {NULL};
+    TProfile* pv2_nest = NULL;
+    TProfile* pv3_nest = NULL;
+    TProfile* pv4_nest = NULL;
   } cor;
 
   struct PhiHist {
@@ -486,12 +493,24 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to
 
     // before loop over particles
     float phi = 0, weight = 1.;
+    vector<float>().swap(cor.vecphi);
+    vector<float>().swap(cor.vecwei);
     for (int ih = 0; ih < maxHarmonic; ih++) {
       for (int ip = 0; ip < maxPower; ip++) {
         cor.Qvector[ih][ip] = TComplex(0., 0.);
       }
     }
 
+    /*
+    for (int ih = 0; ih < maxHarmonic; ih++) {
+      for (int ip = 0; ip < maxPower; ip++) {
+        LOGF(info, "Qvector[%d][%d]=%f, ", ih, ip, cor.Qvector[ih][ip].Rho2());
+      }
+      LOGF(info, "\n");
+    }
+    LOGF(info, "\n\n");
+    */
+
     // Main loop over particles:
     for (const auto& track : tracks) {
 
@@ -527,6 +546,7 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to
         if (cfToyModel) {
           phi = f->GetRandom();
         }
+        cor.vecphi.push_back(phi);
         if (it != phih.histMap.end()) {
           it->second->Fill(phi);
         }
@@ -541,6 +561,7 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to
         } else {
           weight = 1;
         }
+        cor.vecwei.push_back(weight);
 
         // ... and corresponding MC truth simulated:
         // See https://github.com/AliceO2Group/O2Physics/blob/master/Tutorials/src/mcHistograms.cxx
@@ -570,6 +591,30 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to
       }
     } // end of for (auto track: tracks)
     event.fHistMsel[eAfter]->Fill(msel);
+
+    if (cfNest) {
+      float phi1 = 0., phi2 = 0., weight1 = 1., weight2 = 2.;
+      for (Int_t c = 0; c < maxHarmonic; c++) {
+        delete cor.nestedLoops[c];
+        cor.nestedLoops[c] = new TProfile("", "", 1, 0., 1.);
+        cor.nestedLoops[c]->Sumw2();
+      }
+      for (int i1 = 0; i1 < static_cast<int>(cor.vecphi.size()); i1++) { // nested loop of particles
+        phi1 = cor.vecphi[i1];
+        weight1 = cor.vecwei[i1];
+        for (int i2 = 0; i2 < static_cast<int>(cor.vecphi.size()); i2++) {
+          if (i2 == i1) {
+            continue;
+          }
+          phi2 = cor.vecphi[i2];
+          weight2 = cor.vecwei[i2];
+          cor.nestedLoops[0]->Fill(0.5, TMath::Cos(2 * phi1 - 2 * phi2), weight1 * weight2);
+          cor.nestedLoops[1]->Fill(0.5, TMath::Cos(3 * phi1 - 3 * phi2), weight1 * weight2);
+          cor.nestedLoops[2]->Fill(0.5, TMath::Cos(4 * phi1 - 4 * phi2), weight1 * weight2);
+        }
+      } // end of two nested loop
+    }
+
     // calculate correlations
     float Mmin = 4.;
     if (msel < Mmin)
@@ -593,6 +638,13 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to
     cor.pfour32_centr->Fill(centr, four32, wFour);
     cor.pfour42_centr->Fill(centr, four42, wFour);
 
+    if (cfNest) {
+      cor.pv2_nest->Fill(centr, cor.nestedLoops[0]->GetBinContent(1), wTwo);
+      cor.pv3_nest->Fill(centr, cor.nestedLoops[1]->GetBinContent(1), wTwo);
+      cor.pv4_nest->Fill(centr, cor.nestedLoops[2]->GetBinContent(1), wTwo);
+
+      LOGF(info, "v22=%f, v22_nest=%f", v22, cor.nestedLoops[0]->GetBinContent(1));
+    }
   } // end of template <eRecSim rs, typename T1, typename T2> void Steer(T1 const& collision, T2 const& tracks)
 
   // *) Initialize and book all objects:
@@ -827,6 +879,9 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to
     cor.pv42_centr = new TProfile("pv42", "profile of v_{4}^{2}", 9, quantiles);
     cor.pfour32_centr = new TProfile("pfour32", "profile of v_{2}^{2}*v_{3}^{2}", 9, quantiles);
     cor.pfour42_centr = new TProfile("pfour42", "profile of v_{2}^{2}*v_{4}^{2}", 9, quantiles);
+    cor.pv2_nest = new TProfile("pv2_nest", "profile of v_{2} from nest", 9, quantiles);
+    cor.pv3_nest = new TProfile("pv3_nest", "profile of v_{3} from nest", 9, quantiles);
+    cor.pv4_nest = new TProfile("pv4_nest", "profile of v_{4} from nest", 9, quantiles);
     cor.pv22_centr->GetYaxis()->SetTitle("v_{2}^{2}");
     cor.pv32_centr->GetYaxis()->SetTitle("v_{3}^{2}");
     cor.pv42_centr->GetYaxis()->SetTitle("v_{4}^{2}");
@@ -837,11 +892,22 @@ struct MultiharmonicCorrelations { // this name is used in lower-case format to
     cor.pfour42_centr->GetYaxis()->SetTitle("v_{2}^{2}v_{4}^{2}");
     cor.pfour32_centr->GetXaxis()->SetTitle("centrality");
     cor.pfour42_centr->GetXaxis()->SetTitle("centrality");
+    cor.pv2_nest->GetYaxis()->SetTitle("v_{2}");
+    cor.pv3_nest->GetYaxis()->SetTitle("v_{3}");
+    cor.pv4_nest->GetYaxis()->SetTitle("v_{4}");
+    cor.pv2_nest->GetXaxis()->SetTitle("centrality");
+    cor.pv3_nest->GetXaxis()->SetTitle("centrality");
+    cor.pv4_nest->GetXaxis()->SetTitle("centrality");
     cor.fCorrelationVariablesList->Add(cor.pv22_centr);
     cor.fCorrelationVariablesList->Add(cor.pv32_centr);
     cor.fCorrelationVariablesList->Add(cor.pv42_centr);
     cor.fCorrelationVariablesList->Add(cor.pfour32_centr);
     cor.fCorrelationVariablesList->Add(cor.pfour42_centr);
+    if (cfNest) {
+      cor.fCorrelationVariablesList->Add(cor.pv2_nest);
+      cor.fCorrelationVariablesList->Add(cor.pv3_nest);
+      cor.fCorrelationVariablesList->Add(cor.pv4_nest);
+    }
 
     // init of phi hist for different runs
     for (const int& run : targetRuns) {