Difference between revisions of "Analysis Examples"

From GlueXWiki
Jump to: navigation, search
(→‎Histogramming ω Invariant Mass in b1pi)
(→‎Histogramming b1+ Invariant Mass in b1pi)
Line 1: Line 1:
−
== Histogramming b<sub>1</sub><sup>+</sup> Invariant Mass in b1pi ==
+
== Histogramming b<sub>1</sub>(1235)<sup>+</sup> Invariant Mass in b1pi ==
 
* Assumes the DReaction has been setup as shown in: [[Analysis_DReaction | DReaction]]
 
* Assumes the DReaction has been setup as shown in: [[Analysis_DReaction | DReaction]]
  
Line 12: Line 12:
 
* In your DEventProcessor init() method:
 
* In your DEventProcessor init() method:
 
<syntaxhighlight>
 
<syntaxhighlight>
−
dHist_omegaMass = new TH1D("OmegaInvariantMass", ";#pi^{+}#pi^{-}#pi^{0} Invariant Mass (GeV/c^{2})", 600, 0.2, 1.4);
+
dHist_omegaMass = new TH1D("B1InvariantMass", ";#pi^{+}#omega Invariant Mass (GeV/c^{2})", 600, 0.6, 1.8);
 
</syntaxhighlight>
 
</syntaxhighlight>
  
Line 26: Line 26:
  
 
   // get omega combo step //0 is photoproduction, 1 is omega decay, 2 is pi0 decay
 
   // get omega combo step //0 is photoproduction, 1 is omega decay, 2 is pi0 decay
−
   const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(1);
+
   const DParticleComboStep* locParticleComboStep0 = locParticleCombo->Get_ParticleComboStep(0);
 +
  const DParticleComboStep* locParticleComboStep1 = locParticleCombo->Get_ParticleComboStep(1);
  
−
   // calculate omega p4 with kinematic fit data
+
   // calculate b1+ p4 with kinematic fit data
−
   TLorentzVector locP4(0.0, 0.0, 0.0, 0.0);
+
   TLorentzVector locP4;
−
   locP4 += locParticleComboStep->Get_FinalParticle(0)->lorentzMomentum(); //pi+
+
  locP4 += locParticleComboStep0->Get_FinalParticle(2)->lorentzMomentum(); //pi+
−
   locP4 += locParticleComboStep->Get_FinalParticle(1)->lorentzMomentum(); //pi-
+
   locP4 += locParticleComboStep1->Get_FinalParticle(0)->lorentzMomentum(); //pi+
−
   locP4 += locParticleComboStep->Get_FinalParticle(2)->lorentzMomentum(); //pi0
+
   locP4 += locParticleComboStep1->Get_FinalParticle(1)->lorentzMomentum(); //pi-
 +
   locP4 += locParticleComboStep1->Get_FinalParticle(2)->lorentzMomentum(); //pi0
  
 
   Get_Application()->RootWriteLock(); //ACQUIRE ROOT LOCK!!
 
   Get_Application()->RootWriteLock(); //ACQUIRE ROOT LOCK!!
Line 45: Line 47:
 
=== Thorough: Using an Analysis Action ===
 
=== Thorough: Using an Analysis Action ===
 
* Checks against double-counting, shows/allows histogramming measured or kinfit data.  
 
* Checks against double-counting, shows/allows histogramming measured or kinfit data.  
−
* Action class definition (put in a new header file in your plugin):
+
* See DCustomAction_HistMass_b1_1235 in the [https://halldsvn.jlab.org/repos/trunk/sim-recon/src/programs/Analysis/plugins/b1pi_hists/ b1pi Plugin]
−
<syntaxhighlight>
+
−
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;
+
−
};
+
−
</syntaxhighlight>
+
−
 
+
−
* Initialize() method: create histogram
+
−
<syntaxhighlight>
+
−
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!!
+
−
}
+
−
</syntaxhighlight>
+
−
 
+
−
* Perform_Action() method: fill histogram
+
−
<syntaxhighlight>
+
−
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);
+
−
 
+
−
  // 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
+
−
}
+
−
</syntaxhighlight>
+
  
 
* Add to your DReaction anywhere in the analysis chain:
 
* Add to your DReaction anywhere in the analysis chain:

Revision as of 16:38, 5 January 2014

Histogramming b1(1235)+ 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.
    • Not best practice: intended only as a simple illustration.
  • In your DEventProcessor class definition (header):
TH1D* dHist_omegaMass;
  • In your DEventProcessor init() method:
dHist_omegaMass = new TH1D("B1InvariantMass", ";#pi^{+}#omega Invariant Mass (GeV/c^{2})", 600, 0.6, 1.8);
  • 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 omega decay, 2 is pi0 decay
  const DParticleComboStep* locParticleComboStep0 = locParticleCombo->Get_ParticleComboStep(0);
  const DParticleComboStep* locParticleComboStep1 = locParticleCombo->Get_ParticleComboStep(1);
 
  // calculate b1+ p4 with kinematic fit data
  TLorentzVector locP4;
  locP4 += locParticleComboStep0->Get_FinalParticle(2)->lorentzMomentum(); //pi+
  locP4 += locParticleComboStep1->Get_FinalParticle(0)->lorentzMomentum(); //pi+
  locP4 += locParticleComboStep1->Get_FinalParticle(1)->lorentzMomentum(); //pi-
  locP4 += locParticleComboStep1->Get_FinalParticle(2)->lorentzMomentum(); //pi0
 
  Get_Application()->RootWriteLock(); //ACQUIRE ROOT LOCK!!
  {
    dHist_omegaMass->Fill(locP4.M());
  }
  Get_Application()->RootUnLock(); //RELEASE ROOT LOCK!!
  break;
}

Thorough: Using an Analysis Action

  • Checks against double-counting, shows/allows histogramming measured or kinfit data.
  • See DCustomAction_HistMass_b1_1235 in the b1pi Plugin
  • Add to your DReaction anywhere in the analysis chain:
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)

Kinematic Fit Confidence Level

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

Getting the Thrown Track Topology

  • There are several ways of getting the thrown track information. The most direct way of getting it is to ask JANA for the DMCThrown particles (see below). The DMCThrown::parentid member variable matches the DMCThrown::myid member variable that the particle decayed from. If it did not decay from any particle (e.g. was directly photoproduced by the generator) then the DMCThrown::parentid member will equal 0.
vector<const DMCThrown*> locMCThrowns;
locEventLoop->Get(locMCThrowns);
  • An easier way of getting the information is by asking JANA for the thrown DReaction and DParticleCombo objects. These are accessible by requesting these objects with the "Thrown" tag (see below). These contain the PIDs (DReaction) and track kinematics (DParticleCombo) of the thrown particles, organized by the event's decay chain.
    • NOTE: when producing these, any decays of kaons, pions, muons, and neutrons are ignored.
    • NOTE: if a particle type is not defined in sim-recon/src/libraries/include/particleType.h, then it is ignored when generating these objects. 
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];
Retrieved from "https://halldweb.jlab.org/wiki/index.php?title=Analysis_Examples&oldid=53584"