Adding new tasks and moving old ones for 'Femto3D' package by glromane · Pull Request #3675 · AliceO2Group/O2Physics

glromane · 2023-10-23T12:53:08Z

--Implementing small changes in "singletrackselector" task and moving them to "Femto3D" folder in PWGCF branch;

--adding a QA task that uses tables from "singletrackselector";

--adding a femto task that does mixing, pair creation for obtaining a femtoscopic correlation function;

Please consider the following formatting changes to AliceO2Group#3675

victor-gonzalez

Please, use the same sub-directories structure for Femto3D as the other functional areas, as for example FemtoUniverse

Have in mind that the mandate is that there should be only one derived data datamodel for Femto analysis. There has been extensive work in this area and such a datamodel already exists. We will further iterate on that (see my e-mail)

zchochul · 2023-11-02T10:37:42Z

There are only few comments like that, but please try to reduce the amount of commented code in the future PR. Like here:

O2Physics/PWGCF/Femto3D/femto3dPairTask.cxx

Line 45 in 08eb0e5

    
           // if(PDGcode == 2212) return TDatabasePDG::Instance()->GetParticle(2212)->Mass();

Please consider the following formatting changes to AliceO2Group#3675

glromane · 2023-11-03T18:00:43Z

Dear Victor,

Thanks a lot for today's discussion. As agreed, here we commit our tasks with the proper order.

Sincerely yours, Gleb.

victor-gonzalez

Thanks for reallocating the source files in the suggested structure

For the actual code, please, consider the suggested comments, specially the one encouraging the use of the standard O2Physics mixed events framework and not to replicate in memory what is received from tables

victor-gonzalez · 2023-11-06T09:51:39Z

PWGCF/Femto3D/Core/femto3dPairTask.h

+
+class FemtoParticle
+{
+ public:
+  FemtoParticle() {}
+  FemtoParticle(TLorentzVector* fourmomentum, const float& eta, const float& phi);
+  // FemtoParticle(TLorentzVector* fourmomentum);
+  FemtoParticle(const float& E, const float& px, const float& py, const float& pz, const float& eta, const float& phi);
+  // FemtoParticle(const float& E, const float& px, const float& py, const float& pz);
+  FemtoParticle(const FemtoParticle& obj);
+  FemtoParticle(const FemtoParticle* obj);
+  ~FemtoParticle();
+  FemtoParticle& operator=(const FemtoParticle& obj);
+
+  void Set4momentum(TLorentzVector* fourmomentum) { _fourmomentum = fourmomentum; }
+  void SetEta(const float& eta) { _eta = eta; }
+  void SetPhi(const float& phi) { _phi = phi; }
+  void SetSign(const float& sign) { _sign = sign; }
+  void SetMagField(const float& magField) { _magField = magField; }
+
+  TLorentzVector* Get4momentum() const { return _fourmomentum; }
+  float GetEta() const { return _eta; }
+  float GetPhi() const { return _phi; }
+  float GetPhiStar(const float& radius = 1.2) { return _phi + asin(-0.3 * _magField * _sign * radius / (2.0 * _fourmomentum->Pt())); }
+  float GetMagField() const { return _magField; }
+
+ private:
+  float _eta, _phi, _sign, _magField = 0.0;
+  TLorentzVector* _fourmomentum;
+};
+
+FemtoParticle::FemtoParticle(TLorentzVector* fourmomentum, const float& eta, const float& phi)
+{
+  _fourmomentum = fourmomentum;
+  _eta = eta;
+  _phi = phi;
+}
+
+// FemtoParticle::FemtoParticle(TLorentzVector* fourmomentum){
+//   _fourmomentum = fourmomentum;
+// }
+
+FemtoParticle::FemtoParticle(const float& E, const float& px, const float& py, const float& pz, const float& eta, const float& phi)
+{
+  _fourmomentum = new TLorentzVector(px, py, pz, E);
+  _eta = eta;
+  _phi = phi;
+}
+
+// FemtoParticle::FemtoParticle(const float& E, const float& px, const float& py, const float& pz){
+//   _fourmomentum = new TLorentzVector(*px, *py, *pz, *E);
+// }
+
+FemtoParticle::FemtoParticle(const FemtoParticle& obj)
+{
+  Set4momentum(obj.Get4momentum());
+  SetEta(obj.GetEta());
+  SetPhi(obj.GetPhi());
+}
+
+FemtoParticle::FemtoParticle(const FemtoParticle* obj)
+{
+  Set4momentum(obj->Get4momentum());
+  SetEta(obj->GetEta());
+  SetPhi(obj->GetPhi());
+}
+
+FemtoParticle::~FemtoParticle()
+{
+}
+
+FemtoParticle& FemtoParticle::operator=(const FemtoParticle& obj)
+{
+  if (this != &obj) {
+    Set4momentum(obj.Get4momentum());
+    SetEta(obj.GetEta());
+    SetPhi(obj.GetPhi());
+  }
+
+  return *this;
+}
+
+//====================================================================================
+
+class FemtoPair
+{
+ public:
+  FemtoPair(){};
+  FemtoPair(FemtoParticle* first, FemtoParticle* second)
+  {
+    _first = first;
+    _second = second;
+  }
+  FemtoPair(FemtoParticle* first, FemtoParticle* second, const bool& isidentical)
+  {
+    _first = first;
+    _second = second;
+    _isidentical = isidentical;
+  }
+
+  FemtoPair(const FemtoPair& obj)
+  {
+    SetFirstParticle(obj.GetFirstParticle());
+    SetSecondParticle(obj.GetSecondParticle());
+  }
+  FemtoPair(const FemtoPair* obj)
+  {
+    SetFirstParticle(obj->GetFirstParticle());
+    SetSecondParticle(obj->GetSecondParticle());
+  }
+  ~FemtoPair() {}
+  FemtoPair& operator=(const FemtoPair& obj)
+  {
+    if (this != &obj) {
+      SetFirstParticle(obj.GetFirstParticle());
+      SetSecondParticle(obj.GetSecondParticle());
+    }
+    return *this;
+  }
+
+  void SetFirstParticle(FemtoParticle* first) { _first = first; }
+  void SetSecondParticle(FemtoParticle* second) { _second = second; }
+  void SetIdentical(const bool& isidentical) { _isidentical = isidentical; }
+
+  FemtoParticle* GetFirstParticle() const { return _first; }
+  FemtoParticle* GetSecondParticle() const { return _second; }
+  bool IsIdentical() { return _isidentical; }
+
+  bool IsClosePair(const float& deta = 0.01, const float& dphi = 0.01, const float& radius = 1.2);
+  float GetEtaDiff() const { return _first->GetEta() - _second->GetEta(); }
+  float GetPhiStarDiff(const float& radius = 1.2) const { return _first->GetPhiStar(radius) - _second->GetPhiStar(radius); }
+  float GetKstar() const;
+
+ private:
+  FemtoParticle* _first;
+  FemtoParticle* _second;
+  bool _isidentical = true;
+};
+
+bool FemtoPair::IsClosePair(const float& deta, const float& dphi, const float& radius)
+{
+  if (_first == NULL || _second == NULL)
+    return true;
+  if (!(_first->GetMagField() * _second->GetMagField()))
+    return true;
+  if (abs(GetEtaDiff()) < deta || abs(GetPhiStarDiff(radius)) < dphi)
+    return true;
+
+  return false;
+}
+
+float FemtoPair::GetKstar() const
+{
+  if (_first == NULL || _second == NULL)
+    return -1000;
+  if (!(_first->GetMagField() * _second->GetMagField()))
+    return -1000;
+
+  TLorentzVector* first4momentum = new TLorentzVector(*(_first->Get4momentum()));
+  TLorentzVector* second4momentum = new TLorentzVector(*(_second->Get4momentum()));
+
+  if (_isidentical) {
+    TLorentzVector fourmomentadiff = *first4momentum - *second4momentum;
+    return 0.5 * abs(fourmomentadiff.Mag());
+  } else {
+    TLorentzVector fourmomentasum = *first4momentum + *second4momentum;
+
+    first4momentum->Boost((-1) * fourmomentasum.BoostVector());
+    second4momentum->Boost((-1) * fourmomentasum.BoostVector());
+
+    TVector3 qinv = first4momentum->Vect() - second4momentum->Vect();
+    return 0.5 * abs(qinv.Mag());
+  }
+}


If you need these structures they should be under an appropriate namespace

Fixed, put them in the "o2::aod::singletrackselector::" namespace.

victor-gonzalez · 2023-11-06T09:53:20Z

PWGCF/Femto3D/DataModel/singletrackselector.h

+} // namespace singletrackselector
+
+DECLARE_SOA_TABLE_FULL(SingleTrackSel, "SelTracks", "AOD", "STSEL", // Table of the variables for single track selection.
+                       o2::soa::Index<>,


Please, keep the plural for the table and the singular for the iterator

victor-gonzalez · 2023-11-06T09:58:17Z

PWGCF/Femto3D/DataModel/singletrackselector.h

+
+template <typename TrackType>
+inline bool TPCselection(TrackType const& track, std::pair<int, std::vector<float>> const& PIDcuts)
+{


Do you really need the pair structure?
It seems it is created for an ephemeral life and at a per particle basis.
This will definitely have an impact in your job execution.

This structure contains a PDG/PIDcuts pair and if I'm not mistaken the pair is created in only ones at the beginning of a task (QA of PairTask) and this function only reads it since it's by a "const" reference.

victor-gonzalez · 2023-11-06T10:00:42Z

PWGCF/Femto3D/TableProducer/singleTrackSelector.cxx

+        // if satisfied, want to continue in the upper loop (over tracks) -- skip the current track
+        // cannot use simple 'continue' since it will be applied to the current loop, so have to use 'goto'
+        if (TOFselection(track, std::make_pair(i, rejectWithinNsigmaTOF))) {
+          goto skip;


Please, no goto
Use meaningful name variables to control the execution flow

I put a flag instead. No 'goto' anymore.

victor-gonzalez · 2023-11-06T10:01:52Z

PWGCF/Femto3D/TableProducer/singleTrackSelector.cxx

+          break; // break the loop with particlesToKeep after the 'if' condition is satisfied -- don't want double entries
+        }
+
+    skip:; // here we send with 'goto'


Please, no labels
Use meaningful name variables to control the execution flow

victor-gonzalez · 2023-11-06T10:17:32Z

PWGCF/Femto3D/Tasks/femto3dPairTask.cxx

+          Particle = new FemtoParticle(track.energy(particle_mass(_particlePDG_1)), track.px(), track.py(), track.pz(), track.eta(), track.phi());
+          Particle->SetSign(_sign_1);
+          Particle->SetMagField(collision.magField());
+          selectedtracks_1[track.singleCollSelId()].push_back(Particle);
+
+          registry.fill(HIST("p_first"), track.p());
+          if (_particlePDG_1 == 2212) {
+            registry.fill(HIST("nsigmaTOF_first"), track.p(), track.tofNSigmaPr());
+            registry.fill(HIST("nsigmaTPC_first"), track.p(), track.tpcNSigmaPr());
+          }
+          if (_particlePDG_1 == 1000010020) {
+            registry.fill(HIST("nsigmaTOF_first"), track.p(), track.tofNSigmaDe());
+            registry.fill(HIST("nsigmaTPC_first"), track.p(), track.tpcNSigmaDe());
+          }


You are replicating in memory everything that you get from the tables with a much more inefficient way of accessing to the stored information.
This will have a huge penalty in execution time and in memory requirements
This goes against the O2 execution model (SoA vs AoS)
It seems that one of the reasons for storing the information is for extracting the mixed event results
You have to use the mixed events framework that is standard in O2Physics
It works without any kind of data replication

Maybe it doesn't seem right but I'll try to explain why I think it's more convenient to use it like this. The main reason is the mixing (I put a more detailed comment about the mixing in the main reply). Since we need to store some info per particle for the mixing later I decided to use this "FemtoParticle" object that contains the bare minimum: pT, Eta, Phi, magField and signedPDG (updated in the last commit wrt to the previous one -- now consumes less memory).

Actually I was trying to store iterators/pointers to the tracks but the ones that are created after slicing or filtering are transient (temporary) ones and cannot be used to access an object (for the mixing) after the selection. So, I don't see any other way to do it apart from creating another additional table which would lead to duplication anyway. But using the class "FemtoPatricle" is more convenient because I can put there needed functions, pass the object to the "FemtoPair" etc.

victor-gonzalez · 2023-11-06T10:18:08Z

PWGCF/Femto3D/Tasks/femto3dPairTask.cxx

+          Particle = new FemtoParticle(track.energy(particle_mass(_particlePDG_2)), track.px(), track.py(), track.pz(), track.eta(), track.phi());
+          Particle->SetSign(_sign_2);
+          Particle->SetMagField(collision.magField());
+          selectedtracks_2[track.singleCollSelId()].push_back(Particle);
+
+          registry.fill(HIST("p_second"), track.p());
+          if (_particlePDG_2 == 2212) {
+            registry.fill(HIST("nsigmaTOF_second"), track.p(), track.tofNSigmaPr());
+            registry.fill(HIST("nsigmaTPC_second"), track.p(), track.tpcNSigmaPr());
+          }
+          if (_particlePDG_2 == 1000010020) {
+            registry.fill(HIST("nsigmaTOF_second"), track.p(), track.tofNSigmaDe());
+            registry.fill(HIST("nsigmaTPC_second"), track.p(), track.tpcNSigmaDe());
+          }


Same comment as above

victor-gonzalez · 2023-11-06T10:20:34Z

PWGCF/Femto3D/Tasks/femto3dPairTask.cxx

+    if (IsIdentical) { // start of same event identical
+
+      for (auto i = selectedtracks_1.begin(); i != selectedtracks_1.end(); i++) { // start of same event identical
+        if ((i->second).size() < 2)
+          continue;
+
+        std::vector<std::vector<int>> comb_idx;
+        comb(comb_idx, (i->second).size(), 2);
+
+        for (auto indx : comb_idx) {
+          FemtoPair* Pair = new FemtoPair((i->second)[indx[0]], (i->second)[indx[1]], IsIdentical);
+          if (!Pair->IsClosePair(_deta, _dphi, _radiusTPC))
+            registry.fill(HIST("SE"), Pair->GetKstar()); // close pair separation doesn't work properly, need eta and phi at some position in the barrel
+        }
+      } // end of same event identical
+
+      for (auto i = mixbins.begin(); i != mixbins.end(); i++) { // start of mixed event identical
+        if ((i->second).size() < 2)
+          continue;
+
+        std::vector<std::vector<int>> comb_idx;
+        comb(comb_idx, (i->second).size(), 2);
+
+        for (auto indx : comb_idx) {
+          int colId_1 = (i->second)[indx[0]];
+          int colId_2 = (i->second)[indx[1]];
+
+          for (auto ii : selectedtracks_1[colId_1]) {
+            for (auto iii : selectedtracks_1[colId_2]) {
+              FemtoPair* Pair = new FemtoPair(ii, iii, IsIdentical);
+              if (!Pair->IsClosePair(_deta, _dphi, _radiusTPC))
+                registry.fill(HIST("ME"), Pair->GetKstar());
+            }
+          }
+        }
+      } // end of mixed event identical
+
+    } //====================================== end of mixing identical ======================================
+
+    else { //====================================== mixing non-identical ======================================
+
+      for (auto i = selectedtracks_1.begin(); i != selectedtracks_1.end(); i++) { // start of same event non-identical
+
+        auto j = selectedtracks_2.find(i->first);
+        if (j == selectedtracks_2.end())
+          continue;
+
+        for (auto ii : i->second) {
+          for (auto iii : j->second) {
+            FemtoPair* Pair = new FemtoPair(ii, iii, IsIdentical);
+            if (!Pair->IsClosePair(_deta, _dphi, _radiusTPC))
+              registry.fill(HIST("SE"), Pair->GetKstar()); // close pair separation doesn't work properly, need eta and phi at some position in the barrel
+          }
+        }
+      } // end of same event non-identical
+
+      for (auto i = mixbins.begin(); i != mixbins.end(); i++) { // start of mixed event non-identical
+        if ((i->second).size() < 2)
+          continue;
+
+        std::vector<std::vector<int>> comb_idx;
+        comb(comb_idx, (i->second).size(), 2);
+
+        for (auto indx : comb_idx) {
+          int colId_1 = (i->second)[indx[0]];
+          int colId_2 = (i->second)[indx[1]];
+
+          for (auto ii : selectedtracks_1[colId_1]) {
+            for (auto iii : selectedtracks_2[colId_2]) {
+              FemtoPair* Pair = new FemtoPair(ii, iii, IsIdentical);
+              if (!Pair->IsClosePair(_deta, _dphi, _radiusTPC))
+                registry.fill(HIST("ME"), Pair->GetKstar());
+            }
+          }
+        }
+      } // end of mixed event non-idetical
+
+    } //====================================== end of mixing non-identical ======================================


Same comment as above

Here you're right. I fixed this and now I use only one "FemtoPair" object and just change the stored objects inside it. Thanks for pointing it out to me.

Please consider the following formatting changes to AliceO2Group#3675

glromane · 2023-11-10T15:25:22Z

Thanks for reallocating the source files in the suggested structure

For the actual code, please, consider the suggested comments, specially the one encouraging the use of the standard O2Physics mixed events framework and not to replicate in memory what is received from tables

Dear Victor, thanks for your comments, here is my overall reply (I also commented on your separate comments).

I took into account your suggestions and applied them in the new commit.
Small things as namespace, goto etc. I fixed. For the "FemtoParticle", "FemtoPair" objects and mixing procedure I'm putting a more detailed comment:

First, for the tables I cannot use filtering for the PID selection which means that even in the filtered with the track cuts table I will have lots of unwanted entries. And if I pass the table with the "fake" tracks to the mixer (as in the FemtoDream) I'll have lots of "fake" pairs to process the PID selection for. The number of pairs scales as n*(n-1) (x2 because it's a pair) and I think it's very inefficient and instead I do the PID before the mixing and pairs I build doesn't need to be checked anymore. So, instead of 2n(n-1) PID checks, I do it only n-times. Yes, for that I need to store the info and allocate memory but I gain a lot in CPU efficiency. And for the use of O2 mixing algorithm after the PID I'd need to make another table and duplicate the objects anyway, so it's the same only use of the "FemtoParticle" class is much more convenient for me.

There are also some other for wrt to the previous commit made for the sake of efficiency (for example I got rid of the function that does the permutations and replaced it with just nested loops).

We will continue working on improvements and of course, if we don't meet the hyper loop requirements we will deal with it, but before we need to have our tasks there to test and see how to proceed.

If there is something unclear, please, tell me, I'll give you more details. Thank you.

Sincerely yours, Gleb.

victor-gonzalez · 2023-11-15T22:37:10Z

We cannot approve the code as it is
One of the rules we have is to approve code that could be use as reference for other people
Clearly this code is against the O2 rules

Consider this procedure. Once you have your filtered collisions/tracks, your derived data, consider using an additional column you can join to your filtered tracks table which will tag your filtered tracks according to your selection criteria; then use the standard O2 mixed events after filtering the filtered tracks over the additional column

If you plan to elaborate more on your code exclude everything which is not the filter
I will approve the filter. Then, running the filter on hyperloop will allow you to get enough derived data for your local tests in case you need them

jgrosseo · 2023-11-16T19:05:45Z

Let me comment as well. Indeed memory copy needs to be avoided here. For example instead of creating a FemtoPair instance (and then deleting it after use) make a templated function isClosePair which takes two track iterators and gives you the answer. An example can be seen here: https://github.com/jgrosseo/O2Physics/blob/master/PWGCF/Core/PairCuts.h#L114

Very important is also to use the O2 mixing event framework. This was created to avoid the buffers which we had in AliPhysics for storing the objects which slowed down the processing significantly. It works very well, you can define the mixing bins you need yourself and you get the events to mix with automatically.

jgrosseo

See: #3675 (comment)

jgrosseo · 2023-11-20T08:12:06Z

Just to add: if you think something is missing in the O2 event mixing framework, please let us know. We can organize a short slot in the WP4/14 meeting and it can be discussed if it cannot be resolved by simply replying here...

glromane requested review from alibuild, jgrosseo, saganatt, victor-gonzalez and zchochul as code owners October 23, 2023 12:53

glromane added a commit to glromane/O2Physics that referenced this pull request Oct 23, 2023

Merge pull request #1 from alibuild/alibot-cleanup-3675

bcfce38

Please consider the following formatting changes to AliceO2Group#3675

victor-gonzalez requested changes Oct 25, 2023

View reviewed changes

glromane force-pushed the master branch from 08eb0e5 to aa7a8fa Compare November 3, 2023 17:56

glromane added a commit to glromane/O2Physics that referenced this pull request Nov 3, 2023

Merge pull request #3 from alibuild/alibot-cleanup-3675

8faa771

Please consider the following formatting changes to AliceO2Group#3675

victor-gonzalez requested changes Nov 6, 2023

View reviewed changes

glromane force-pushed the master branch from 8faa771 to 7b84cdc Compare November 10, 2023 15:17

glromane added a commit to glromane/O2Physics that referenced this pull request Nov 10, 2023

Merge pull request #4 from alibuild/alibot-cleanup-3675

897144a

Please consider the following formatting changes to AliceO2Group#3675

jgrosseo requested changes Nov 16, 2023

View reviewed changes

glromane mentioned this pull request Nov 17, 2023

Femto 3D table producer #3893

Merged

glromane closed this Dec 21, 2023

glromane force-pushed the master branch from 114e56d to 2e1b15a Compare December 21, 2023 17:10

Conversation

glromane commented Oct 23, 2023

Uh oh!

victor-gonzalez left a comment

Choose a reason for hiding this comment

Uh oh!

zchochul commented Nov 2, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

glromane commented Nov 3, 2023

Uh oh!

victor-gonzalez left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

glromane commented Nov 10, 2023

Uh oh!

victor-gonzalez commented Nov 15, 2023

Uh oh!

jgrosseo commented Nov 16, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jgrosseo left a comment

Choose a reason for hiding this comment

Uh oh!

jgrosseo commented Nov 20, 2023

Uh oh!

Reviewers

Assignees

Labels

Milestone

Development

Uh oh!

4 participants

zchochul commented Nov 2, 2023 •

edited

Loading

jgrosseo commented Nov 16, 2023 •

edited

Loading