Analysis Examples

Histogramming ω Invariant Mass in b1pi

• Assumes the DReaction has been setup as shown in: DReaction

Minimal/Simplest

• Minimal: No check against double-counting, using kinematic fit data only.
• In your DEventProcessor class definition (header):

TH1D* dHist_omegaMass;

• In your DEventProcessor init() method:

dHist_omegaMass = new TH1D("OmegaInvariantMass", ";#pi^{+}#pi^{-}#pi^{0} Invariant Mass (GeV/c^{2})", 600, 0.2, 1.4);

• In your DEventProcessor evnt() method:

vector<const DParticleCombo*> locParticleCombos;
locEventLoop->Get(locParticleCombos);
 
for(size_t loc_i = 0; loc_i < locParticleCombos.size(); ++loc_i)
{
  if(locParticleCombos[loc_i]->Get_Reaction()->Get_ReactionName() != "b1pi")
    continue; //wrong DReaction
 
  // get omega combo step //0 is photoproduction, 1 is X(2000) decay, 2 is b1 decay, etc.
  const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(3);
 
  // calculate omega p4 with kinematic fit data
  TLorentzVector locP4(0.0, 0.0, 0.0, 0.0);
  locP4 += locParticleComboStep->Get_FinalParticle(0)->lorentzMomentum(); //pi+
  locP4 += locParticleComboStep->Get_FinalParticle(1)->lorentzMomentum(); //pi-
  locP4 += locParticleComboStep->Get_FinalParticle(2)->lorentzMomentum(); //pi0
 
  Get_Application()->RootWriteLock(); //ACQUIRE ROOT LOCK!!
  {
    dHist_omegaMass->Fill(locP4.M());
  }
  Get_Application()->RootUnLock(); //RELEASE ROOT LOCK!!
  break;
}

Using an Analysis Action

• Action class definition (put in a new header file in your plugin):

class DHistogramAction_omegaMass : public DAnalysisAction
{
  //assumes b1pi reaction
  public:
    DHistogramAction_omegaMass(const DReaction* locReaction, bool locUseKinFitResultsFlag, string locActionUniqueString = "") :
    DAnalysisAction(locReaction, "Hist_omegaMass", locUseKinFitResultsFlag, locActionUniqueString) {} //"Hist_omegaMass" should be unique to this class
 
  private:
    bool Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo);
    void Initialize(JEventLoop* locEventLoop);
 
    TH1D* dHist_omegaMass;
};

• Initialize() method: create histogram

void DHistogramAction_omegaMass::Initialize(JEventLoop* locEventLoop)
{
  //CREATE THE HISTOGRAMS
  Get_Application()->RootWriteLock(); //ACQUIRE ROOT LOCK!!
 
  CreateAndChangeTo_ActionDirectory(); //go to the correct directory
 
  string locHistName = "InvariantMass";
  if(gDirectory->Get(locHistName.c_str()) != NULL) //already created by another thread
    dHist_omegaMass = static_cast<TH1D*>(gDirectory->Get(locHistName.c_str()));
  else
    dHist_omegaMass = new TH1D(locHistName.c_str(), ";#pi^{+}#pi^{-}#pi^{0} Invariant Mass (GeV/c^{2})", 600, 0.2, 1.4);
 
  Get_Application()->RootUnLock(); //RELEASE ROOT LOCK!!
}

• Perform_Action() method: fill histogram

bool DHistogramAction_omegaMass::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
{
  // get omega combo step //0 is photoproduction, 1 is X(2000) decay, 2 is b1 decay, etc.
  size_t locOmegaStepIndex = 3;
  const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(locOmegaStepIndex);
 
  // check to make sure not double-counting
  if(!Get_UseKinFitResultsFlag()) //set in constructor
  {
    //measured data (kinematic fit results will be unique for each DParticleCombo)
    deque<pair<const DParticleCombo*, bool> > locPreviousParticleCombos;
    Get_PreviousParticleCombos(locPreviousParticleCombos); // DParticleCombos for which this action has been previously executed for this event: compare against these
    //can either compare manually or by using the DAnalysisUtilities methods
    if(Get_AnalysisUtilities()->Find_SimilarCombos(locParticleCombo, locOmegaStepIndex, locPreviousParticleCombos))
      return true; //dupe: already histed! //return true: this isn't a cut, the combo didn't fail
  }
 
  // two ways to calculate invariant mass: 
  TLorentzVector locP4_v1 = Get_AnalysisUtilities()->Calc_FinalStateP4(locParticleCombo, locOmegaStepIndex, Get_UseKinFitResultsFlag());
  // or manually:
  TLorentzVector locP4_v2(0.0, 0.0, 0.0, 0.0);
  if(!Get_UseKinFitResultsFlag())
  {
    // measured data
    locP4_v2 += locParticleComboStep->Get_FinalParticle_Measured(0)->lorentzMomentum(); //pi+
    locP4_v2 += locParticleComboStep->Get_FinalParticle_Measured(1)->lorentzMomentum(); //pi-
 
    // pi0 isn't measured: must get data from gamma, gamma
    const DParticleComboStep* locPi0Step = locParticleCombo->Get_ParticleComboStep(4); //pi0 decay
    locP4_v2 += locPi0Step->Get_FinalParticle_Measured(0)->lorentzMomentum(); //gamma
    locP4_v2 += locPi0Step->Get_FinalParticle_Measured(1)->lorentzMomentum(); //gamma
  }
  else
  {
    // kinematic fit data
    locP4_v2 += locParticleComboStep->Get_FinalParticle(0)->lorentzMomentum(); //pi+
    locP4_v2 += locParticleComboStep->Get_FinalParticle(1)->lorentzMomentum(); //pi-
    locP4_v2 += locParticleComboStep->Get_FinalParticle(2)->lorentzMomentum(); //pi0
  }
 
  // lock around histogram filling
  Get_Application()->RootWriteLock(); //ACQUIRE ROOT LOCK!!
  {
    dHist_omegaMass->Fill(locP4_v1.M());
  }
  Get_Application()->RootUnLock(); //RELEASE ROOT LOCK!!
 
  return true; //return true: this isn't a cut, the combo didn't fail
}

• Add to your DReaction:

locReaction->Add_AnalysisAction(new DHistogramAction_omegaMass(locReaction, false)); //false: measured data
locReaction->Add_AnalysisAction(new DHistogramAction_omegaMass(locReaction, true, "KinFit")); //true: kinfit data