Synthetic flow generator added (#1721)

* Synthetic flow generator added

* Add tester for synthetic flow MC

(cherry picked from commit c2b733f)

Author: ddobrigk

MC/config/PWGLF/ini/GeneratorLF_SyntheFlow.ini

@@ -0,0 +1,6 @@
+[GeneratorExternal]
+fileName=${O2DPG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_syntheFlow.C
+funcName=generator_syntheFlow()
+
+[GeneratorPythia8]
+config=${O2DPG_ROOT}/MC/config/common/pythia8/generator/pythia8_hi.cfg

MC/config/PWGLF/ini/tests/GeneratorLF_SyntheFlow.C

@@ -0,0 +1,28 @@
+int External()
+{
+    std::string path{"o2sim_Kine.root"};
+
+    TFile file(path.c_str(), "READ");
+    if (file.IsZombie())
+    {
+        std::cerr << "Cannot open ROOT file " << path << "\n";
+        return 1;
+    }
+
+    auto tree = (TTree *)file.Get("o2sim");
+    if (!tree)
+    {
+        std::cerr << "Cannot find tree o2sim in file " << path << "\n";
+        return 1;
+    }
+    std::vector<o2::MCTrack> *tracks{};
+    tree->SetBranchAddress("MCTrack", &tracks);
+
+    auto nEvents = tree->GetEntries();
+    if (nEvents < 1)
+    {
+        std::cerr << "No events actually generated: not OK!";
+        return 1;
+    }
+    return 0;
+}

MC/config/PWGLF/pythia8/generator_pythia8_syntheFlow.C

@@ -0,0 +1,90 @@
+
+#include "Pythia8/Pythia.h"
+#include "Pythia8/HeavyIons.h"
+#include "FairGenerator.h"
+#include "FairPrimaryGenerator.h"
+#include "Generators/GeneratorPythia8.h"
+#include "TRandom3.h"
+#include "TParticlePDG.h"
+#include "TDatabasePDG.h"
+
+#include <map>
+#include <unordered_set>
+
+class GeneratorPythia8SyntheFlow : public o2::eventgen::GeneratorPythia8
+{
+public:
+  /// Constructor
+  GeneratorPythia8SyntheFlow() {
+    lutGen = new o2::eventgen::FlowMapper();
+    
+    // -------- CONFIGURE SYNTHETIC FLOW ------------
+    // specify a v2 vs pT here
+    TFile *filehep = new TFile("/Users/daviddc/Downloads/HEPData-ins1116150-v1-Table_1.root", "READ");
+    TH1D *hv = (TH1D*) filehep->Get("Table 1/Hist1D_y6");
+    
+    TFile *fileEcc = new TFile("/Users/daviddc/Downloads/eccentricityvsb.root", "READ");
+    TH1D *hEccentricities = (TH1D*) fileEcc->Get("hEccentricities");
+    
+    cout<<"Generating LUT for flow test"<<endl;
+    lutGen->CreateLUT(hv, hEccentricities);
+    cout<<"Finished creating LUT!"<<endl;
+    // -------- END CONFIGURE SYNTHETIC FLOW ------------
+  }
+
+  ///  Destructor
+  ~GeneratorPythia8SyntheFlow() = default;
+  
+  //__________________________________________________________________
+  Bool_t generateEvent() override {
+    
+    // Generate PYTHIA event
+    Bool_t lPythiaOK = kFALSE;
+    while (!lPythiaOK){
+      lPythiaOK = mPythia.next();
+    }
+    
+    //+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
+    // loop over the entire event record and rotate all particles
+    // synthetic flow exercise
+    // first: get event plane
+    float eventPlaneAngle = mPythia.info.hiInfo->phi();
+    float impactParameter = mPythia.info.hiInfo->b();
+
+    for ( Long_t j=0; j < mPythia.event.size(); j++ ) {
+      float pyphi = mPythia.event[j].phi();
+      float pypT = mPythia.event[j].pT();
+
+      // calculate delta with EP
+      float deltaPhiEP = pyphi - eventPlaneAngle;
+      float shift = 0.0;
+      while(deltaPhiEP<0.0){
+        deltaPhiEP += 2*TMath::Pi();
+        shift += 2*TMath::Pi();
+      }
+      while(deltaPhiEP>2*TMath::Pi()){
+        deltaPhiEP -= 2*TMath::Pi();
+        shift -= 2*TMath::Pi();
+      }
+      float newDeltaPhiEP = lutGen->MapPhi(deltaPhiEP, impactParameter, pypT);
+      float pyphiNew = newDeltaPhiEP - shift + eventPlaneAngle;
+
+      if(pyphiNew>TMath::Pi())
+        pyphiNew -= 2.0*TMath::Pi();
+      if(pyphiNew<-TMath::Pi())
+        pyphiNew += 2.0*TMath::Pi();
+      mPythia.event[j].rot(0.0, pyphiNew-pyphi);
+    }
+    //+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
+    
+    return true;
+  }
+  
+private:
+  o2::eventgen::FlowMapper *lutGen; // for mapping phi angles
+};
+
+ FairGenerator *generator_syntheFlow()
+ {
+   return new GeneratorPythia8SyntheFlow();
+ }