DEMPgen/src/DEMPEvent.cxx

0001 #include "DEMPEvent.hxx"
0002 #include "Particle.hxx"
0003 #include "Constants.hxx"
0004
0005 #include "TMath.h"
0006 #include "TVector3.h"
0007 #include "TLorentzVector.h"
0008
0009 #include <iostream>
0010 #include <stdio.h>
0011 #include <string.h>
0012 #include "json/json.h"
0013
0014 using namespace constants;
0015
0016 DEMPEvent::DEMPEvent()
0017 {
0018   qsq_GeV = new double(0);
0019   t_GeV = new double(0);
0020   w_GeV = new double(0);
0021   t_para_GeV = new double(0);
0022   t_prime_GeV = new double(0);
0023   negt = new double(0);
0024   x_B = new double(0);
0025   Phi = new double(0);
0026   Phi_s = new double(0);
0027   Theta = new double(0);
0028   Phi_deg = new double(0);
0029   Phi_s_deg = new double(0);
0030   Theta_deg = new double(0);
0031   P_T = new double(0);
0032   Vertex_x = new double(0);
0033   Vertex_y = new double(0);
0034   Vertex_z = new double(0);
0035
0036   STphi = TMath::Pi();
0037   //std::cout<<STphi<<std::endl;
0038
0039   extern Json::Value obj;
0040
0041   TargPol = new TLorentzVector(obj["targ_pol_x"].asDouble(),
0042                                obj["targ_pol_y"].asDouble(),
0043                                obj["targ_pol_z"].asDouble(),
0044                                0);
0045
0046 }
0047
0048 DEMPEvent::DEMPEvent(const char* prefix)
0049 {
0050   qsq_GeV = new double(0);
0051   t_GeV = new double(0);
0052   w_GeV = new double(0);
0053   t_para_GeV = new double(0);
0054   t_prime_GeV = new double(0);
0055   negt = new double(0);
0056   x_B = new double(0);
0057   Phi = new double(0);
0058   Phi_s = new double(0);
0059   Theta = new double(0);
0060   Phi_deg = new double(0);
0061   Phi_s_deg = new double(0);
0062   Theta_deg = new double(0);
0063   P_T = new double(0);
0064   Vertex_x = new double(0);
0065   Vertex_y = new double(0);
0066   Vertex_z = new double(0);
0067
0068   STphi = TMath::Pi();
0069   //std::cout<<STphi<<std::endl;
0070
0071   char be_name[100];
0072   strcpy(be_name, prefix);
0073   strcat(be_name, "BeamElec");
0074   BeamElec = new Particle(electron_mass_mev, be_name, pid_elec);
0075   BeamElec->SetCharge(-1);
0076
0077   char tn_name[100];
0078   strcpy(tn_name, prefix);
0079   strcat(tn_name, "TargNeut");
0080   TargNeut = new Particle(neutron_mass_mev, tn_name, pid_neut);
0081   TargNeut->SetCharge(0);
0082
0083   char se_name[100];
0084   strcpy(se_name, prefix);
0085   strcat(se_name, "ScatElec");
0086   ScatElec = new Particle(electron_mass_mev, se_name, pid_elec);
0087   ScatElec->SetCharge(-1);
0088
0089   char ppi_name[100];
0090   strcpy(ppi_name, prefix);
0091   strcat(ppi_name, "ProdPion");
0092   ProdPion = new Particle(pion_mass_mev, ppi_name, pid_pion);
0093   ProdPion->SetCharge(-1);
0094
0095   char ppr_name[100];
0096   strcpy(ppr_name, prefix);
0097   strcat(ppr_name, "ProdProt");
0098   ProdProt = new Particle(proton_mass_mev, ppr_name, pid_prot);
0099   ProdProt->SetCharge(1);
0100
0101   /*
0102   char vp_name[100];
0103   strcpy(vp_name, prefix);
0104   strcat(vp_name, "VirtPhot");
0105   ProdProt = new Particle(0, vp_name, pid_prot);
0106   ProdProt->SetCharge(0);
0107   */
0108   VirtPhot = new Particle();
0109
0110   extern Json::Value obj;
0111
0112   TargPol = new TLorentzVector(obj["targ_pol_x"].asDouble(),
0113                                obj["targ_pol_y"].asDouble(),
0114                                obj["targ_pol_z"].asDouble(),
0115                                0);
0116
0117 }
0118
0119 void DEMPEvent::Update()
0120 {
0121   *VirtPhot = *BeamElec - *ScatElec;
0122
0123   *qsq_GeV = -(VirtPhot->Mag2())/1000000;
0124
0125   *w_GeV = (*VirtPhot+*TargNeut).Mag()/1000;
0126
0127   *t_GeV = (*VirtPhot-*ProdPion).Mag2()/1000000;
0128
0129
0130   double E1 = VirtPhot->E();
0131   double P1 = VirtPhot->P();
0132
0133   double E2 = ProdPion->E();
0134   double P2 = ProdPion->P();
0135
0136   *t_para_GeV = (((E1-E2)*(E1-E2))-((P1-P2)*(P1-P2)))/1000000;
0137
0138   *t_prime_GeV = (*t_GeV) - (*t_para_GeV);
0139
0140   *negt = (-1.0)*(*t_GeV);
0141
0142   *x_B = (*qsq_GeV)*1000000/(2*proton_mass_mev *
0143                              this->VirtPhot->E());
0144
0145   *Phi_s = TargPol->Phi() - ScatElec->Phi();
0146   *Phi_s_deg = *Phi_s * DEG;
0147
0148   *Phi = ProdPion->Phi();
0149   *Phi_deg = *Phi * DEG;
0150
0151   *Theta = VirtPhot->Theta();
0152   *Theta_deg = *Theta * DEG;
0153
0154   *P_T = TMath::Abs(TargPol->Perp(VirtPhot->Vect()));
0155
0156 }
0157
0158 // Calculate and return the center of momentum 3-vector
0159 TVector3 DEMPEvent::CoM()
0160 {
0161   return (BeamElec->Vect() + TargNeut->Vect())*(1/(BeamElec->E()+TargNeut->E()));
0162 }
0163
0164 // Copy each particle object from another event.
0165 DEMPEvent DEMPEvent::operator = (const DEMPEvent& q)
0166 {
0167   *(this->BeamElec) = *(q.BeamElec);
0168   *(this->TargNeut) = *(q.TargNeut);
0169   *(this->ScatElec) = *(q.ScatElec);
0170   *(this->ProdPion) = *(q.ProdPion);
0171   *(this->ProdProt) = *(q.ProdProt);
0172   *(this->VirtPhot) = *(q.VirtPhot);
0173   *(this->TargPol) = *(q.TargPol);
0174
0175   this->STphi = q.STphi;
0176 }
0177
0178 void DEMPEvent::Boost(TVector3 boostvect)
0179 {
0180   //std::cout << boostvect.Z() << std::endl;
0181   //std::cout << BeamElec->T() << std::endl;
0182   BeamElec->Boost(boostvect);
0183   //std::cout << BeamElec->T() << std::endl;
0184   TargNeut->Boost(boostvect);
0185   ScatElec->Boost(boostvect);
0186   ProdProt->Boost(boostvect);
0187   ProdPion->Boost(boostvect);
0188   VirtPhot->Boost(boostvect);
0189   TargPol->Boost(boostvect);
0190 }
0191
0192 void DEMPEvent::Rotate(double rottheta, double rotphi)
0193 {
0194
0195   //std::cout << ScatElec->Phi() << std::endl;
0196
0197   BeamElec->RotateZ(rotphi);
0198   TargNeut->RotateZ(rotphi);
0199   ScatElec->RotateZ(rotphi);
0200   ProdProt->RotateZ(rotphi);
0201   ProdPion->RotateZ(rotphi);
0202   VirtPhot->RotateZ(rotphi);
0203
0204   TargPol->RotateZ(rotphi);
0205
0206   //std::cout << ScatElec->Phi() << std::endl;
0207
0208   //std::cout<< VirtPhot->Theta()*DEG << std::endl;
0209
0210   BeamElec->RotateY(rottheta);
0211   TargNeut->RotateY(rottheta);
0212   ScatElec->RotateY(rottheta);
0213   ProdProt->RotateY(rottheta);
0214   ProdPion->RotateY(rottheta);
0215   VirtPhot->RotateY(rottheta);
0216
0217   TargPol->RotateY(rottheta);
0218
0219   //std::cout<< VirtPhot->Theta()*DEG << std::endl;
0220
0221   //std::cout<<STphi*DEG<<std::endl;
0222
0223   STphi += rotphi;
0224   //std::cout<<STphi*DEG<<std::endl;
0225
0226 }