ExampleAnalysis_KVEventMixerN_3Body.cpp

Build correlated & uncorrelated excitation energy spectra for \f${}^{7}Be\rightarrow p+d+\alpha\f$ 3-body decay using event mixing technique

Example of use of KVEventMixerN which is a generic implementation of event mixing for N-particle correlation studies. In this example, N=3

#ifndef __CORRELATOR_H
#define __CORRELATOR_H
 
#include "KVReconEventSelector.h"
#include "KVEventMixerN.h"
#include "KVEventClassifier.h"
 
class ExampleAnalysis_KVEventMixerN_3Body : public KVReconEventSelector {
 
   struct particle {
      int z, a, det_index; 
      TVector3 momentum;
      particle() = default;
      particle(particle&&) = default;
      particle(const particle&) = default;
      particle(const KVReconstructedNucleus& p) : z{p.GetZ()},
         a{p.GetA()},
         det_index{p.GetStoppingDetector()->GetIndex()},
         momentum{p.GetMomentum()}
      {}
   };
 
   KVEventMixerN<3, KVReconstructedNucleus, particle, 5> event_mixer{1};
 
   KVEventClassifier* mult_bin;
 
   TString get_cor_histo_name(int bin, const TString& quantity)
   {
      return Form("h_cor_%s_bin_%d", quantity.Data(), bin);
   }
   TString get_uncor_histo_name(int bin, const TString& quantity)
   {
      return Form("h_uncor_%s_bin_%d", quantity.Data(), bin);
   }
   int nevent = 0;
public:
   ExampleAnalysis_KVEventMixerN_3Body() {}
 
   virtual void InitRun();
   virtual void EndRun() {}
   virtual void InitAnalysis();
   virtual Bool_t Analysis();
   virtual void EndAnalysis() {}
 
   ClassDef(ExampleAnalysis_KVEventMixerN_3Body, 0) //User analysis class
};
 
#endif

#include "ExampleAnalysis_KVEventMixerN_3Body.h"
 
ClassImp(ExampleAnalysis_KVEventMixerN_3Body)
 
void ExampleAnalysis_KVEventMixerN_3Body::InitAnalysis(void)
{
   // Declaration of histograms, global variables, etc.
   // Called at the beginning of the analysis
   // The examples given are compatible with interactive, batch,
   // and PROOFLite analyses.
 
   // multiplicities of each type of particle
   auto gv = AddGV("KVMult", "m_proton");
   gv->SetSelection({"proton", [](const KVNucleus * _n)
   {
      const KVReconstructedNucleus* n = dynamic_cast<const KVReconstructedNucleus*>(_n);
      return n->IsOK() && n->IsAMeasured() && n->IsIsotope(1, 1);
   }});
   // select events with at least 1 proton
   gv->SetEventSelection([](const KVVarGlob * vg) {
      return vg->GetValue() > 0;
   });
 
   gv = AddGV("KVMult", "m_deuton");
   gv->SetSelection({"deuton", [](const KVNucleus * _n)
   {
      const KVReconstructedNucleus* n = dynamic_cast<const KVReconstructedNucleus*>(_n);
      return n->IsOK() && n->IsAMeasured() && n->IsIsotope(1, 2);
   }});
   // select events with at least 1 deuton
   gv->SetEventSelection([](const KVVarGlob * vg) {
      return vg->GetValue() > 0;
   });
 
   gv = AddGV("KVMult", "m_alpha");
   gv->SetSelection({"alpha", [](const KVNucleus * _n)
   {
      const KVReconstructedNucleus* n = dynamic_cast<const KVReconstructedNucleus*>(_n);
      return n->IsOK() && n->IsAMeasured() && n->IsIsotope(2, 4);
   }});
   // select events with at least 1 alpha
   gv->SetEventSelection([](const KVVarGlob * vg) {
      return vg->GetValue() > 0;
   });
 
   AddGV("KVMult", "mtot");   // total multiplicity for event sorting
   mult_bin = GetGVList()->AddEventClassifier("mtot");
   // Multiplicity cuts corresponding to 5 equal statistics bins of Xe+Sn 50AMeV
   std::vector<Double_t> mult_slices = {8.4195047, 14.449584, 21.173604, 28.044737};
   for (auto cut : mult_slices) mult_bin->AddCut(cut);
 
   for (size_t i = 0; i <= mult_slices.size(); ++i) {
      AddHisto<TH1F>(get_cor_histo_name(i, "Ex"), Form("Correlated spectrum E* bin %ld", i), 500, 0., 100.);
      AddHisto<TH1F>(get_uncor_histo_name(i, "Ex"), Form("Uncorrelated spectrum E* bin %ld", i), 500, 0., 100.);
   }
 
   auto t = AddTree("check_tree", "check event classifier");
   GetGVList()->MakeBranches(t);
 
   /*** DEFINE WHERE TO SAVE THE RESULTS ***/
   SetJobOutputFileName("ExampleAnalysis_KVEventMixerN_3Body_results.root");
}
 
void ExampleAnalysis_KVEventMixerN_3Body::InitRun(void)
{
   // Reject events with less identified particles than the acquisition multiplicity trigger
   SetTriggerConditionsForRun(GetCurrentRun()->GetNumber());
}
 
Bool_t ExampleAnalysis_KVEventMixerN_3Body::Analysis(void)
{
   auto bin = mult_bin->GetEventClassification();
 
   GetGVList()->FillBranches();
   FillTree();
 
   event_mixer.ProcessEvent(bin,
                            // correlated particles function, 'CorFunc'
   [ = ](int bin_num, KVRefVec<KVReconstructedNucleus> parts) {
      // calculate E* for p+d+a from same event
      KVNucleus BE;
      for (auto& n : parts) {
         BE += n;
      }
      FillHisto(get_cor_histo_name(bin_num, "Ex"), BE.GetExcitEnergy());
   },
   // uncorrelated particles function, 'UnCorFunc'
   [ = ](int bin_num, KVReconstructedNucleus & part, KVRefVec<particle> other_parts) {
      // calculate E* for 'part' from current event, other 2 particles from previous events
      KVNucleus BE;
      BE += part;
      for (auto& n : other_parts) {
         BE += KVNucleus(n.z, n.a, n.momentum);
      }
      FillHisto(get_uncor_histo_name(bin_num, "Ex"), BE.GetExcitEnergy());
   },
   // iterator for protons
   ReconEventIterator(GetEvent(), {
      "proton",
      [](const KVReconstructedNucleus * n)
      {
         return n->IsOK() && n->IsAMeasured() && n->IsIsotope(1, 1);
      }
   }),
   // iterator for deuterons
   ReconEventIterator(GetEvent(), {
      "deuton",
      [](const KVReconstructedNucleus * n)
      {
         return n->IsOK() && n->IsAMeasured() && n->IsIsotope(1, 2);
      }
   }),
   // iterator for alphas
   ReconEventIterator(GetEvent(), {
      "alpha",
      [](const KVReconstructedNucleus * n)
      {
         return n->IsOK() && n->IsAMeasured() && n->IsIsotope(2, 4);
      }
   })
 
                           );
 
   return kTRUE;
}