// $Id$ // // File: JEventProcessor_cpp_analysis.cc // Created: Thu Feb 28 17:03:05 EST 2013 // Creator: davidl (on Darwin harriet.jlab.org 11.4.2 i386) // #include #include "JEventProcessor_cpp_analysis.h" using namespace jana; #include #include //------------------ // JEventProcessor_cpp_analysis (Constructor) //------------------ JEventProcessor_cpp_analysis::JEventProcessor_cpp_analysis() { } //------------------ // ~JEventProcessor_cpp_analysis (Destructor) //------------------ JEventProcessor_cpp_analysis::~JEventProcessor_cpp_analysis() { } //------------------ // init //------------------ jerror_t JEventProcessor_cpp_analysis::init(void) { phi_pipi = new TH1D("phi_pipi", "#phi of #pi#pi system", 300, -M_PI, M_PI); phi_pipi_thrown = new TH1D("phi_pipi_thrown", "#phi of #pi#pi system (generated)", 300, -M_PI, M_PI); W_pipi = new TH1D("W_pipi", "W_{#pi#pi}", 500, 250.0, 550.0); W_pipi_thrown = new TH1D("W_pipi_thrown", "W_{#pi#pi}", 500, 250.0, 550.0); phi_pipi_res = new TH1D("phi_pipi_res", "#phi resolution of #pi#pi system", 300, -1000.0, 1000.0); theta_pipi_res = new TH1D("theta_pipi_res", "#theta resolution of #pi#pi system", 5000, -30.0, 500.0); pt_pipi_res = new TH1D("pt_pipi_res", "p_{t} resolution of #pi#pi system", 300, -1.0, 1.0); W_pipi_res = new TH1D("W_pipi_res", "W_{#pi#pi} resolution", 1000, -100.0, 100.0); missing_mass = new TH1D("missing_mass", "Missing Mass", 4000.0, 0.0, 400.0); theta_pip_cm = new TH1D("theta_pip_cm", "Reconstructed #pi^{+} #theta in CM", 300, 0.0, 180.0); theta_pip_cm_thrown = new TH1D("theta_pip_cm_thrown", "Generated #pi^{+} #theta in CM", 300, 0.0, 180.0); theta_pip_cm_res = new TH1D("theta_pip_cm_res", "Reconstructed #pi^{+} #theta resolution in CM", 1000, -500.0, +500.0); phi_pipi_res->SetXTitle("#delta#phi (mrad)"); theta_pipi_res->SetXTitle("#delta#phi (mrad)"); pt_pipi_res->SetXTitle("#deltap_{t}/p_{t}"); missing_mass->SetXTitle("missing mass (GeV/c^2)"); W_pipi_res->SetXTitle("#DeltaW_{#pi#pi} (MeV/c^2)"); W_pipi->SetXTitle("W_{#pi#pi} (MeV/c^2)"); W_pipi_thrown->SetXTitle("W_{#pi#pi} (MeV/c^2)"); theta_pip_cm->SetXTitle("CM angle (degrees)"); theta_pip_cm_thrown->SetXTitle("CM angle (degrees)"); theta_pip_cm_res->SetXTitle("#delta#theta_{CM} (mrad)"); return NOERROR; } //------------------ // brun //------------------ jerror_t JEventProcessor_cpp_analysis::brun(JEventLoop *eventLoop, int runnumber) { // This is called whenever the run number changes return NOERROR; } //------------------ // evnt //------------------ jerror_t JEventProcessor_cpp_analysis::evnt(JEventLoop *loop, int eventnumber) { // Get generated values for pi+ pi- DLorentzVector pip(0.0, 0.0, 0.0, 0.0); DLorentzVector pim(0.0, 0.0, 0.0, 0.0); double pip_cm_theta_degrees_thrown = -1000.0; vector throwns; loop->Get(throwns); for(unsigned int i=0; itype == PiPlus )pip = thrown->lorentzMomentum(); if(thrown->type == PiMinus)pim = thrown->lorentzMomentum(); } DLorentzVector pipi_thrown = pip + pim; if(pip.P()!=0.0 && pim.P()!=0.0){ phi_pipi_thrown->Fill(pipi_thrown.Phi()); W_pipi_thrown->Fill(pipi_thrown.M()*1000.0); // in MeV TLorentzVector pip_cm = pip; pip_cm.Boost(-pipi_thrown.BoostVector()); pip_cm_theta_degrees_thrown = pip_cm.Theta()*57.3; theta_pip_cm_thrown->Fill(pip_cm_theta_degrees_thrown); } // Have ANALYSIS library process all defined reactions vector locAnalysisResultsVector; loop->Get(locAnalysisResultsVector); // Loop over analysis results for(size_t loc_i = 0; loc_i < locAnalysisResultsVector.size(); ++loc_i) { // Check if this result has any particle combinations at all const DAnalysisResults* locAnalysisResults = locAnalysisResultsVector[loc_i]; if(locAnalysisResults->Get_NumPassedParticleCombos() == 0) continue; // pi+ pi- reaction if(locAnalysisResults->Get_Reaction()->Get_ReactionName() == "pi+pi-"){ // Get pointers to all particle combos passing our criteria deque locPassedParticleCombos; locAnalysisResults->Get_PassedParticleCombos(locPassedParticleCombos); if(locPassedParticleCombos.size()>0){ const DParticleComboStep *locStep = locPassedParticleCombos[0]->Get_ParticleComboStep(0); if(locStep->Get_NumFinalParticles() >= 2){ DLorentzVector pip = locStep->Get_FinalParticle(0)->lorentzMomentum(); DLorentzVector pim = locStep->Get_FinalParticle(1)->lorentzMomentum(); // Pb208? const DKinematicData *missing_kd = locStep->Get_MissingParticle(); if(missing_kd){ DLorentzVector missing = missing_kd->lorentzMomentum(); missing_mass->Fill(missing_kd->mass()); } // Histogram results DLorentzVector pipi = pip + pim; phi_pipi->Fill(pipi.Phi()); W_pipi->Fill(pipi.M()*1000.0); // in MeV double delta_phi = pipi.Phi() - pipi_thrown.Phi(); // in radians double delta_theta = pipi.Theta() - pipi_thrown.Theta(); // in radians double delta_pt = pipi.Pt() - pipi_thrown.Pt(); // in GeV/c double delta_W = pipi.M() - pipi_thrown.M(); phi_pipi_res->Fill(delta_phi * 1000.0); // in mrad theta_pipi_res->Fill(delta_theta * 1000.0); // in mrad pt_pipi_res->Fill(delta_pt/pipi_thrown.Pt()); // in GeV W_pipi_res->Fill(delta_W*1000.0); // in MeV TLorentzVector pip_cm = pip; pip_cm.Boost(-pipi.BoostVector()); double pip_cm_theta_degrees = pip_cm.Theta()*57.3; theta_pip_cm->Fill(pip_cm_theta_degrees); double delta_cm_theta_degrees = pip_cm_theta_degrees - pip_cm_theta_degrees_thrown; theta_pip_cm_res->Fill(delta_cm_theta_degrees/57.3*1000.0); // convert to milliradians } } } // pi0 reaction if(locAnalysisResults->Get_Reaction()->Get_ReactionName() == "pi0"){ // Get pointers to all particle combos passing our criteria deque locPassedParticleCombos; locAnalysisResults->Get_PassedParticleCombos(locPassedParticleCombos); } } return NOERROR; } //------------------ // erun //------------------ jerror_t JEventProcessor_cpp_analysis::erun(void) { // This is called whenever the run number changes, before it is // changed to give you a chance to clean up before processing // events from the next run number. return NOERROR; } //------------------ // fini //------------------ jerror_t JEventProcessor_cpp_analysis::fini(void) { // Called before program exit after event processing is finished. return NOERROR; }