13 #include <TDirectoryFile.h>
14 #include <TLorentzVector.h>
16 #include <JANA/JApplication.h>
17 #include <JANA/JEventLoop.h>
22 #include <PID/DParticle.h>
23 #include <PID/DPhoton.h>
53 mm_gp_to_pX =
new TH1D(
"mm_gp_to_pX",
"Missing mass from #gamma p -> p X",4001, 0.0, 4.0);
54 mm_gp_to_pX->SetXTitle(
"Missing mass (GeV)");
55 mm_gp_to_pX_thrown = (TH1D*)mm_gp_to_pX->Clone(
"mm_gp_to_pX_thrown");
57 t_pX =
new TH1D(
"t_pX",
"-t for #gammap->pX",20, 0.0, 6.0);
58 t_pX->SetXTitle(
"-t (GeV)");
60 sqrt_s =
new TH1D(
"sqrt_s",
"Center of mass energy #sqrt{s}",2001, 0.0, 6.0);
61 sqrt_s->SetXTitle(
"#sqrt{s} C.M. energy (GeV)");
65 tree =
new TTree(
"t",
"#rho p events");
66 tree->Branch(
"E",&evt);
68 pthread_mutex_init(&mutex, NULL);
79 vector<const DBeamPhoton*> beam_photons;
80 vector<const DParticle*> particles;
81 vector<const DParticle*> particles_thrown;
83 loop->Get(beam_photons);
85 loop->Get(particles_thrown,
"THROWN");
88 TLorentzVector target(0.0, 0.0, 0.0, 0.93827);
91 vector<TLorentzVector> rec_piplus;
92 vector<TLorentzVector> rec_piminus;
93 vector<TLorentzVector> rec_protons;
94 SortChargedParticles(particles, rec_piplus, rec_piminus, rec_protons);
97 vector<TLorentzVector> thrown_piplus;
98 vector<TLorentzVector> thrown_piminus;
99 vector<TLorentzVector> thrown_protons;
100 SortChargedParticles(particles_thrown, thrown_piplus, thrown_piminus, thrown_protons);
104 if(beam_photons.size()==0){
112 beam_photons.push_back(&beam);
119 japp->RootWriteLock();
122 evt->event = eventnumber;
125 for(
unsigned int i=0; i<beam_photons.size(); i++){
128 TLorentzVector beam_photon = MakeTLorentz(beam_photons[i], 0.0);
129 evt->beam = beam_photon;
132 sqrt_s->Fill((beam_photon+target).M());
135 if(rec_protons.size()==1){
136 TLorentzVector missing = beam_photon + target - rec_protons[0];
137 mm_gp_to_pX->Fill(missing.M());
141 if(thrown_protons.size()==1){
142 TLorentzVector missing = beam_photon + target - thrown_protons[0];
143 mm_gp_to_pX_thrown->Fill(missing.M());
155 if(thrown_piplus.size()!=1)
160 if(thrown_piminus.size()!=1)
167 evt->rho_thrown.isfiducial = IsFiducial(thrown_piplus[0]) && IsFiducial(thrown_piminus[0]);
170 evt->rho_thrown.pip = thrown_piplus[0];
171 evt->rho_thrown.pim = thrown_piminus[0];
172 evt->rho_thrown.m = (evt->rho_thrown.pip+evt->rho_thrown.pim).M();
175 for(
unsigned int j=0; j<rec_piplus.size(); j++){
176 for(
unsigned int k=0; k<rec_piminus.size(); k++){
177 evt->AddRho(rec_piplus[j], rec_piminus[k]);
182 if(thrown_protons.size()==1)evt->proton_thrown = thrown_protons[0];
185 if(rec_protons.size()==1){
186 TLorentzVector beam_photon = beam_photons.size()>0 ? MakeTLorentz(beam_photons[0], 0.0):TLorentzVector(0.0, 0.0, 0.0, 9.0);
187 TLorentzVector &proton = rec_protons[0];
188 TLorentzVector p3 = beam_photon + target - proton;
189 double t = (beam_photon - p3).Mag2();
210 for(
unsigned int j=0; j<particles.size(); j++){
211 const DParticle *part = particles[j];
217 int type = part->charge()<0.0 ? 9:8;
221 vector<const DMCThrown*> throwns;
224 if(throwns.size()>0)type = throwns[0]->type;
228 case 8: rec_piplus.push_back(MakeTLorentz(part, 0.13957));
break;
229 case 9: rec_piminus.push_back(MakeTLorentz(part, 0.13957));
break;
230 case 14: rec_protons.push_back(MakeTLorentz(part, 0.93827));
break;
246 double theta = kd->
momentum().Theta();
248 double px = p*
sin(theta)*cos(phi);
250 double pz = p*cos(theta);
251 double E =
sqrt(mass*mass + p*p);
253 return TLorentzVector(px,py,pz,E);
261 double theta = pion.Theta()*TMath::RadToDeg();
262 if(theta<2.0 || theta>110.0)
return false;
263 if(pion.P()<0.500)
return false;
jerror_t fini(void)
Invoked via DEventProcessor virtual method.
void setMomentum(const DVector3 &aMomentum)
jerror_t init(void)
Invoked via DEventProcessor virtual method.
jerror_t evnt(JEventLoop *loop, uint64_t eventnumber)
Invoked via DEventProcessor virtual method.
void SortChargedParticles(vector< const DParticle * > &particles, vector< TLorentzVector > &rec_piplus, vector< TLorentzVector > &rec_piminus, vector< TLorentzVector > &rec_protons)
const DVector3 & momentum(void) const
jerror_t erun(void)
Invoked via DEventProcessor virtual method.
void setPosition(const DVector3 &aPosition)
bool IsFiducial(TLorentzVector &pion)
TLorentzVector MakeTLorentz(const DKinematicData *track, double mass)