Difference between revisions of "Analysis Examples"

TH1D* dHist_omegaMass;

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

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;
}

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;
};

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!!
}

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
}

locReaction->Add_AnalysisAction(new DHistogramAction_omegaMass(locReaction, false)); //false: measured data
locReaction->Add_AnalysisAction(new DHistogramAction_omegaMass(locReaction, true, "KinFit")); //true: kinfit data //"KinFit:" a unique, distinguishing string (could have been anything)

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
  const DKinFitResults* locKinFitResults = locParticleCombos[loc_i]->Get_KinFitResults();
  double locConfidenceLevel = locKinFitResults->Get_ConfidenceLevel();
}

vector<const DMCThrown*> locMCThrowns;
locEventLoop->Get(locMCThrowns);

vector<const DReaction*> locThrownReactions;
locEventLoop->Get(locThrownReactions, "Thrown");
const DReaction* locThrownReaction = locThrownReactions[0];
deque<Particle_t> locThrownPIDs;
locThrownReaction->Get_DetectedFinalPIDs(locThrownPIDs, true); //true tells it to include duplicates (so if two pi-'s, list it twice).  
 
vector<const DParticleCombo*> locThrownParticleCombos;
locEventLoop->Get(locThrownParticleCombos, "Thrown");
const DParticleCombo* locThrownParticleCombo = locThrownParticleCombos[0];