KVTemplateEvent.h

/***************************************************************************
                          kvevent.h  -  description
                             -------------------
    begin                : Sun May 19 2002
    copyright            : (C) 2002 by J.D. Frankland
    email                : frankland@ganil.fr
 
$Id: KVEvent.h,v 1.29 2008/12/17 11:23:12 ebonnet Exp $
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
 
#ifndef KVEVENT_H
#define KVEVENT_H
 
#define KVEVENT_PART_INDEX_OOB "Particle index %d out of bounds [1,%d]"
 
#include "TTree.h"
#include "TVector3.h"
#include "TClonesArray.h"
#include "KVEvent.h"
#include "KVConfig.h"
#include "TRotation.h"
#include "TLorentzRotation.h"
#include "KVTemplateParticleCondition.h"
#include "KVNameValueList.h"
#include "TMethodCall.h"
 
#include <KVFrameTransform.h>
#include <KVIntegerList.h>
#include <TH1.h>
#include <TObjString.h>
#include <iterator>
 
class KVIntegerList;
 
template <typename Particle>
class KVTemplateEvent: public KVEvent {
 
public:
 
   class Iterator {
   public:
      typedef std::forward_iterator_tag iterator_category;
      typedef Particle value_type;
      typedef std::ptrdiff_t difference_type;
      typedef Particle* pointer;
      typedef Particle& reference;
 
      enum Type {    // type of iterator
         Null,      // null value
         All,       // include all particles
         OK,        // include particles which are "OK"
         Group,      // include particles belonging to group
         Bad      // type mismatch: iterator goes straight to end
      };
 
   private:
      KVTemplateParticleCondition<Particle> fSelection;//condition for selecting particles
      TIter   fIter;//iterator over TClonesArray
      Type    fType;//iterator type
      mutable Bool_t  fIterating;//=kTRUE when iteration in progress
      Bool_t AcceptableIteration()
      {
 
         if (fType == Type::Bad)
            return kFALSE;
         return fSelection.Test(current());
      }
      Particle* current() const
      {
         return static_cast<Particle*>(*fIter);
      }
   public:
      Iterator()
         : fSelection(),
           fIter(static_cast<TIterator*>(nullptr)),
           fType(Type::Null),
           fIterating(kFALSE)
      {}
      Iterator(const Iterator& i)
         : fSelection(i.fSelection),
           fIter(i.fIter),
           fType(i.fType),
           fIterating(i.fIterating)
      {}
 
      Iterator(const KVEvent* e, const KVTemplateParticleCondition<Particle>& selection)
         : fSelection(selection), fIter(e->GetParticleArray()), fIterating(kTRUE)
      {
 
         if (!e->GetParticleArray()->GetClass()->InheritsFrom(Particle::Class())) {
            ::Warning("KVTemplateEvent::Iterator", "KVTemplateEvent<%s>::Iterator for %s particles requested for event containing %s particles. Iteration is aborted.",
                      Particle::Class()->GetName(), Particle::Class()->GetName(), e->GetParticleArray()->GetClass()->GetName());
            fType = Bad;
         }
 
         fIter.Begin();// set iterator to first particle of event corresponding to selection
         while ((current() != nullptr) && !AcceptableIteration()) ++fIter;
         if (current() == nullptr) fIterating = kFALSE;
      }
 
      Iterator(const KVEvent& e, const KVTemplateParticleCondition<Particle>& selection)
         : fSelection(selection), fIter(e.GetParticleArray()), fIterating(kTRUE)
      {
 
         if (!e.GetParticleArray()->GetClass()->InheritsFrom(Particle::Class())) {
            ::Warning("KVTemplateEvent::Iterator", "KVTemplateEvent<%s>::Iterator for %s particles requested for event containing %s particles. Iteration is aborted.",
                      Particle::Class()->GetName(), Particle::Class()->GetName(), e.GetParticleArray()->GetClass()->GetName());
            fType = Bad;
         }
 
         fIter.Begin();// set iterator to first particle of event corresponding to selection
         while ((current() != nullptr) && !AcceptableIteration()) ++fIter;
         if (current() == nullptr) fIterating = kFALSE;
      }
 
      Iterator(const KVEvent* e, Type t = Type::All, TString grp = "")
         : fSelection(), fIter(e->GetParticleArray()), fType(t), fIterating(kTRUE)
      {
 
         if (!e->GetParticleArray()->GetClass()->InheritsFrom(Particle::Class())) {
            ::Warning("KVTemplateEvent::Iterator", "KVTemplateEvent<%s>::Iterator for %s particles requested for event containing %s particles. Iteration is aborted.",
                      Particle::Class()->GetName(), Particle::Class()->GetName(), e->GetParticleArray()->GetClass()->GetName());
            fType = Bad;
         }
         if (fType == Type::OK) {
            fSelection.Set("ok", [](const Particle * n) {
               return n->IsOK();
            });
         }
         else if (fType == Type::Group) {
            fSelection.Set("group", [grp](const Particle * n) {
               return n->BelongsToGroup(grp);
            });
         }
 
         fIter.Begin();// set iterator to first particle of event corresponding to selection
         while ((current() != nullptr) && !AcceptableIteration()) ++fIter;
         if (current() == nullptr) fIterating = kFALSE;
      }
 
      Iterator(const KVEvent& e, Type t = Type::All, TString grp = "")
         : fSelection(), fIter(e.GetParticleArray()), fType(t), fIterating(kTRUE)
      {
 
         if (!e.GetParticleArray()->GetClass()->InheritsFrom(Particle::Class())) {
            ::Warning("KVTemplateEvent::Iterator", "KVTemplateEvent<%s>::Iterator for %s particles requested for event containing %s particles. Iteration is aborted.",
                      Particle::Class()->GetName(), Particle::Class()->GetName(), e.GetParticleArray()->GetClass()->GetName());
            fType = Bad;
         }
         if (fType == Type::OK) {
            fSelection.Set("ok", [](const Particle * n) {
               return n->IsOK();
            });
         }
         else if (fType == Type::Group) {
            fSelection.Set("group", [grp](const Particle * n) {
               return n->BelongsToGroup(grp);
            });
         }
 
         fIter.Begin();// set iterator to first particle of event corresponding to selection
         while ((current() != nullptr) && !AcceptableIteration()) ++fIter;
         if (current() == nullptr) fIterating = kFALSE;
      }
 
      Particle& operator* () const
      {
 
         return *(current());
      }
      template<typename PointerType = Particle>
      PointerType * get_pointer() const
      {
         return dynamic_cast<PointerType*>(current());
      }
      template<typename ReferenceType = Particle>
      ReferenceType & get_reference() const
      {
         return dynamic_cast<ReferenceType&>(*current());
      }
      template<typename PointerType = Particle>
      const PointerType * get_const_pointer() const
      {
         return dynamic_cast<const PointerType*>(current());
      }
      template<typename ReferenceType = Particle>
      const ReferenceType & get_const_reference() const
      {
         return dynamic_cast<const ReferenceType&>(*current());
      }
      Bool_t operator!= (const Iterator& it) const
      {
         return current() != it.current();
      }
      Bool_t operator== (const Iterator& it) const
      {
         return current() == it.current();
      }
      const Iterator& operator++ ()
      {
         if (current() == nullptr) fIterating = kFALSE;
         ++fIter;
         while (current() != nullptr && !AcceptableIteration()) ++fIter;
         return *this;
      }
      Iterator operator++ (int)
      {
         Iterator tmp(*this);
         operator++();
         return tmp;
      }
      Iterator& operator= (const Iterator& rhs)
      {
         if (this != &rhs) { // check self-assignment based on address of object
            fSelection = rhs.fSelection;
            fIter = rhs.fIter;
            fType = rhs.fType;
            fIterating = rhs.fIterating;
         }
         return *this;
      }
 
      static Iterator End()
      {
         return Iterator();
      }
 
      void Reset(Type t = Type::Null, TString grp = "")
      {
         if (t != Type::Null) {
            fType = t;
            if (fType == Type::OK) {
               fSelection.Set("ok", [](const Particle * n) {
                  return n->IsOK();
               });
            }
            else if (fType == Type::Group) {
               fSelection.Set("group", [grp](const Particle * n) {
                  return n->BelongsToGroup(grp);
               });
            }
         }
         fIter.Begin();
         fIterating = kTRUE;
         while ((current() != nullptr) && !AcceptableIteration()) ++fIter;
      }
      Bool_t IsIterating() const
      {
         return fIterating;
      }
      void SetIsIterating(Bool_t on = kTRUE)
      {
         fIterating = on;
      }
   };
protected:
   mutable Iterator fIter;
 
public:
   KVTemplateEvent(Int_t mult = 50)
      : KVEvent(Particle::Class(), mult)
   {}
 
   Particle* AddParticle()
   {
 
      Int_t mult = GetMult();
#ifdef __WITHOUT_TCA_CONSTRUCTED_AT
      Particle* tmp = (Particle*) ConstructedAt(mult, "C");
#else
      Particle* tmp = (Particle*) fParticles->ConstructedAt(mult, "C");
#endif
      if (!tmp) {
         Error("AddParticle", "Allocation failure, Mult=%d", mult);
         return 0;
      }
      return tmp;
   }
   Particle* GetParticle(Int_t npart) const
   {
 
      if (npart < 1 || npart > fParticles->GetEntriesFast()) {
         Error("GetParticle", KVEVENT_PART_INDEX_OOB, npart,
               fParticles->GetEntriesFast());
         return 0;
      }
 
      return (Particle*)((*fParticles)[npart - 1]);
   }
   virtual Int_t GetMult(Option_t* opt = "") const
   {
 
      if (TString(opt) == "") return KVEvent::GetMult();
      Int_t fMultOK = 0;
      for (Iterator it = GetNextParticleIterator(opt); it != end(); ++it) ++fMultOK;
      return fMultOK;
   }
   Int_t GetMultiplicity(Int_t Z, Int_t A = 0, Option_t* opt = "")
   {
 
      if (A > 0) return (Int_t)GetSum("IsIsotope", "int,int", Form("%d,%d", Z, A), opt);
      return (Int_t)GetSum("IsElement", "int", Form("%d", Z), opt);
   }
   void GetMultiplicities(Int_t mult[], const TString& species, Option_t* opt = "")
   {
 
      std::unique_ptr<TObjArray> spec(species.Tokenize(", "));// remove any spaces
      Int_t nspec = spec->GetEntries();
      memset(mult, 0, nspec * sizeof(Int_t)); // set multiplicities to zero
      for (Iterator it = GetNextParticleIterator(opt); it != end(); ++it) {
         for (int i = 0; i < nspec; i++) {
            if (((TObjString*)(*spec)[i])->String() == (*it).GetSymbol()) mult[i] += 1;
         }
      }
   }
   Double_t GetSum(const Char_t* Nucleus_method, Option_t* opt = "")
   {
 
      Double_t fSum = 0;
      TMethodCall mt;
      mt.InitWithPrototype(Particle::Class(), Nucleus_method, "");
 
      if (mt.IsValid()) {
         Iterator it = GetNextParticleIterator(opt);
         if (mt.ReturnType() == TMethodCall::kLong) {
            Long_t ret;
            for (; it != end(); ++it) {
               Particle* tmp = it.get_pointer();
               mt.Execute(tmp, "", ret);
               fSum += ret;
            }
         }
         else if (mt.ReturnType() == TMethodCall::kDouble) {
            Double_t ret;
            for (; it != end(); ++it) {
               Particle* tmp = it.get_pointer();
               mt.Execute(tmp, "", ret);
               fSum += ret;
            }
         }
      }
 
      return fSum;
   }
   Double_t GetSum(const Char_t* Nucleus_method, const Char_t* method_prototype, const Char_t* args, Option_t* opt = "")
   {
 
      Double_t fSum = 0;
      TMethodCall mt;
      mt.InitWithPrototype(Particle::Class(), Nucleus_method, method_prototype);
 
      if (mt.IsValid()) {
         Iterator it = GetNextParticleIterator(opt);
         if (mt.ReturnType() == TMethodCall::kLong) {
            Long_t ret;
            for (; it != end(); ++it) {
               Particle* tmp = it.get_pointer();
               mt.Execute(tmp, args, ret);
               fSum += ret;
            }
         }
         else if (mt.ReturnType() == TMethodCall::kDouble) {
            Double_t ret;
            for (; it != end(); ++it) {
               Particle* tmp = it.get_pointer();
               mt.Execute(tmp, args, ret);
               fSum += ret;
            }
         }
      }
 
      return fSum;
   }
   void FillHisto(TH1* h, const Char_t* Nucleus_method, Option_t* opt = "")
   {
 
      TMethodCall mt;
      mt.InitWithPrototype(Particle::Class(), Nucleus_method, "");
 
      if (mt.IsValid()) {
         Iterator it = GetNextParticleIterator(opt);
         if (mt.ReturnType() == TMethodCall::kLong) {
            Long_t ret;
            for (; it != end(); ++it) {
               Particle* tmp = it.get_pointer();
               mt.Execute(tmp, "", ret);
               h->Fill((Double_t)ret);
            }
         }
         else if (mt.ReturnType() == TMethodCall::kDouble) {
            Double_t ret;
            for (; it != end(); ++it) {
               Particle* tmp = it.get_pointer();
               mt.Execute(tmp, "", ret);
               h->Fill(ret);
            }
         }
      }
   }
   void FillHisto(TH1* h, const Char_t* Nucleus_method, const Char_t* method_prototype, const Char_t* args, Option_t* opt = "")
   {
 
      TMethodCall mt;
      mt.InitWithPrototype(Particle::Class(), Nucleus_method, method_prototype);
 
      if (mt.IsValid()) {
         Iterator it = GetNextParticleIterator(opt);
         if (mt.ReturnType() == TMethodCall::kLong) {
            Long_t ret;
            for (; it != end(); ++it) {
               Particle* tmp = it.get_pointer();
               mt.Execute(tmp, args, ret);
               h->Fill((Double_t)ret);
            }
         }
         else if (mt.ReturnType() == TMethodCall::kDouble) {
            Double_t ret;
            for (; it != end(); ++it) {
               Particle* tmp = it.get_pointer();
               mt.Execute(tmp, args, ret);
               h->Fill(ret);
            }
         }
      }
   }
   virtual void Print(Option_t* t = "") const
   {
 
      std::cout << "\nKVEvent with " << GetMult(t) << " particles :" << std::endl;
      std::cout << "------------------------------------" << std::endl;
      fParameters.Print();
      for (Iterator it = GetNextParticleIterator(t); it != end(); ++it)(*it).Print();
   }
   virtual void ls(Option_t* t = "") const
   {
      Print(t);
   }
   Particle* GetParticleWithName(const Char_t* name) const
   {
 
      Particle* tmp = (Particle*)fParticles->FindObject(name);
      return tmp;
   }
   Particle* GetParticle(const Char_t* group_name) const
   {
 
      Iterator it = GetNextParticleIterator(group_name);
      Particle* tmp = it.get_pointer();
      if (!tmp) Warning("GetParticle", "Particle not found: %s", group_name);
      return tmp;
   }
 
   Iterator begin() const
   {
      return Iterator(this);
   }
   Iterator end() const
   {
      return Iterator::End();
   }
 
   Particle* GetNextParticle(Option_t* opt = "") const
   {
 
      TString Opt(opt);
      Opt.ToUpper();
 
      if (fIter.IsIterating()) return &(*(fIter++));
 
      Bool_t ok_iter = (Opt == "OK");
      Bool_t grp_iter = (!ok_iter && Opt.Length());
 
      if (ok_iter) fIter = Iterator(this, Iterator::Type::OK);
      else if (grp_iter) fIter = Iterator(this, Iterator::Type::Group, Opt);
      else fIter = Iterator(this);
      return &(*(fIter++));
   }
   void ResetGetNextParticle() const
   {
 
      fIter.SetIsIterating(kFALSE);
   }
 
   void ResetEnergies()
   {
 
      for (Iterator it = begin(); it != end(); ++it) {
         (*it).ResetEnergy();
      }
   }
 
   void DefineGroup(const Char_t* groupname, const Char_t* from = "")
   {
 
      for (Iterator it = GetNextParticleIterator(from); it != end(); ++it) {
         (*it).AddGroup(groupname);
      }
   }
   void DefineGroup(const Char_t* groupname, KVTemplateParticleCondition<Particle>* cond, const Char_t* from = "")
   {
 
      for (Iterator it = GetNextParticleIterator(from); it != end(); ++it) {
         if (cond->Test(it.get_const_pointer()))(*it).AddGroup(groupname, from);
      }
   }
 
   void SetFrame(const Char_t* frame, const KVFrameTransform& ft)
   {
 
      for (auto& part : *this) {
         part.SetFrame(frame, ft);
      }
   }
   void SetFrame(const Char_t* newframe, const Char_t* oldframe, const KVFrameTransform& ft)
   {
 
      for (auto& part : *this) {
         part.SetFrame(newframe, oldframe, ft);
      }
   }
   void ChangeFrame(const KVFrameTransform& ft, const KVString& name = "")
   {
 
 
      for (auto& part : *this) {
         part.ChangeFrame(ft, name);
      }
      if (name != "") SetParameter("defaultFrame", name);
   }
   void ChangeDefaultFrame(const Char_t* newdef, const Char_t* defname = "")
   {
 
      for (auto& part : *this) {
         part.ChangeDefaultFrame(newdef, defname);
      }
      SetParameter("defaultFrame", newdef);
   }
   void UpdateAllFrames()
   {
 
      for (auto& part : *this) {
         part.UpdateAllFrames();
      }
   }
 
   template<typename U = Particle>
   typename std::enable_if<std::is_base_of<KVNucleus, U>::value>::type
   FillIntegerList(KVIntegerList* IL, Option_t* opt)
   {
 
      IL->Clear();
      for (Iterator it = GetNextParticleIterator(opt); it != end(); ++it) IL->Add((*it).GetZ());
      IL->SetPopulation(1);
      IL->CheckForUpdate();
   }
 
   void GetMasses(std::vector<Double_t>& mass)
   {
 
      mass.clear();
      mass.reserve(GetMult());
      int i = 0;
      for (Iterator it = begin(); it != end(); ++it) mass[i++] = (*it).GetMass();
   }
 
   template<typename U = Particle>
   typename std::enable_if<std::is_base_of<KVNucleus, U>::value>::type
   GetGSMasses(std::vector<Double_t>& mass)
   {
 
      mass.clear();
      mass.reserve(GetMult());
      int i = 0;
      for (Iterator it = begin(); it != end(); ++it) mass[i++] = (*it).GetMassGS();
   }
 
   template<typename U = Particle>
   typename std::enable_if<std::is_base_of<KVNucleus, U>::value, Double_t>::type
   get_channel_qvalue() const
   {
 
      Double_t sumM = 0;
      Particle CN;
      Int_t M = GetMult();
      for (int i = 1; i <= M; i++) {
         sumM += GetParticle(i)->GetMass();
         CN += *(GetParticle(i));
      }
      return CN.GetMassGS() - sumM;
   }
   template<typename U = Particle>
   typename std::enable_if < !std::is_base_of<KVNucleus, U>::value, Double_t >::type
   get_channel_qvalue() const
   {
      return 0;
   }
   Double_t GetChannelQValue() const
   {
      return get_channel_qvalue();
   }
   template<typename U = Particle>
   typename std::enable_if<std::is_base_of<KVNucleus, U>::value, Double_t>::type
   GetGSChannelQValue() const
   {
 
      Double_t sumM = 0;
      Particle CN;
      Int_t M = GetMult();
      for (int i = 1; i <= M; i++) {
         sumM += GetParticle(i)->GetMassGS();
         CN += *(GetParticle(i));
      }
      return CN.GetMassGS() - sumM;
   }
   template<typename U = Particle>
   typename std::enable_if<std::is_base_of<KVNucleus, U>::value, KVString>::type
   get_partition_name()
   {
      fParticles->Sort();
      KVString partition;
 
      KVNameValueList nvl;
      partition = "";
      for (Iterator it = begin(); it != end(); ++it) {
         TString st = (*it).GetSymbol();
         Int_t pop = TMath::Max(nvl.GetIntValue(st.Data()), 0);
         pop += 1;
         nvl.SetValue(st.Data(), pop);
      }
      for (Int_t ii = 0; ii < nvl.GetEntries(); ii += 1) {
         Int_t pop = nvl.GetIntValue(ii);
         if (pop == 1) partition += nvl.GetNameAt(ii);
         else        partition += Form("%s(%d)", nvl.GetNameAt(ii), pop);
         if (ii < nvl.GetEntries() - 1) partition += " ";
      }
      return partition;
   }
 
   template<typename U = Particle>
   typename std::enable_if < !std::is_base_of<KVNucleus, U>::value, KVString >::type
   get_partition_name()
   {
      return "";
   }
 
   KVString GetPartitionName()
   {
      return get_partition_name();
   }
 
   void SetFrameName(const KVString& name)
   {
 
      for (Iterator it = begin(); it != end(); ++it) {
         (*it).SetFrameName(name);
      }
      SetParameter("defaultFrame", name);
   }
 
   struct EventIterator {
      Iterator it;
      EventIterator(const KVEvent& event, typename Iterator::Type t = Iterator::Type::All, const TString& grp = "")
         : it(event, t, grp)
      {}
      EventIterator(const KVEvent* event, typename Iterator::Type t = Iterator::Type::All, const TString& grp = "")
         : it(event, t, grp)
      {}
      EventIterator(const KVEvent& event, const KVTemplateParticleCondition<Particle>& selection)
         : it(event, selection)
      {}
      EventIterator(const KVEvent* event, const KVTemplateParticleCondition<Particle>& selection)
         : it(event, selection)
      {}
      Iterator begin() const
      {
         return it;
      }
      Iterator end() const
      {
         return Iterator::End();
      }
   };
 
   struct EventOKIterator : public EventIterator {
      EventOKIterator(const KVEvent& event) :
         EventIterator(event, Iterator::Type::OK)
      {}
      EventOKIterator(const KVEvent* event) :
         EventIterator(event, Iterator::Type::OK)
      {}
   };
 
   struct EventGroupIterator : public EventIterator {
      EventGroupIterator(const KVEvent& event, const TString& grp) :
         EventIterator(event, Iterator::Type::Group, grp)
      {}
      EventGroupIterator(const KVEvent* event, const TString& grp) :
         EventIterator(event, Iterator::Type::Group, grp)
      {}
   };
 
   Iterator GetNextParticleIterator(Option_t* opt) const
   {
 
      TString Opt(opt);
      Opt.ToUpper();
      Bool_t ok_iter = (Opt == "OK");
      Bool_t grp_iter = (!ok_iter && Opt.Length());
 
      if (ok_iter) return EventOKIterator(this).begin();
      else if (grp_iter) return EventGroupIterator(this, Opt).begin();
      return EventIterator(this).begin();
   }
 
   EventIterator ConditionalIterator(const KVTemplateParticleCondition<Particle>& c)
   {
      return EventIterator(this, c);
   }
 
   ClassDef(KVTemplateEvent, 0)         //Base class for all types of multiparticle event
};
 
#endif