DEMPgen/src/SigmaCalc.cxx

0001 #include <vector>
0002 #include "json/json.h"
0003
0004 #include "TMath.h"
0005
0006 #include "DEMPEvent.hxx"
0007 #include "Asymmetry.hxx"
0008
0009 #include "SigmaL.hxx"
0010 #include "SigmaT.hxx"
0011 #include "SigmaLT.hxx"
0012 #include "SigmaTT.hxx"
0013 #include "Constants.hxx"
0014
0015 #include "correction.h"
0016
0017 #include "SigmaCalc.hxx"
0018
0019 using namespace std;
0020 using namespace constants;
0021 using namespace TMath;
0022
0023 extern Json::Value obj;
0024
0025
0026 SigmaCalc::SigmaCalc(DEMPEvent* in_VertEvent,
0027                      DEMPEvent* in_CofMEvent,
0028                      DEMPEvent* in_RestEvent,
0029                      DEMPEvent* in_TConEvent):
0030   VertEvent(in_VertEvent),
0031   CofMEvent(in_CofMEvent),
0032   RestEvent(in_RestEvent),
0033   TConEvent(in_TConEvent)
0034 {
0035
0036
0037   vector<double> qsq;
0038   qsq.push_back(4.107);
0039   qsq.push_back(4.335);
0040   qsq.push_back(4.845);
0041   qsq.push_back(5.138);
0042   qsq.push_back(5.557);
0043   qsq.push_back(5.948);
0044   qsq.push_back(6.049);
0045   qsq.push_back(6.393);
0046   qsq.push_back(6.778);
0047   qsq.push_back(6.894);
0048   qsq.push_back(7.617);
0049
0050   Asyms = new vector<Asymmetry*>(5);
0051
0052   Asyms->at(0) =
0053     new Asymmetry("asy",
0054                   "[0]*exp([1]*x)+(-[2]-[0])*exp([3]*x)+[2]",
0055                   qsq, false);
0056
0057   Asyms->at(1) =
0058     new Asymmetry("asysfi",
0059                   "[0]*exp([1]*x)+[2]",
0060                   qsq, false);
0061
0062   Asyms->at(2) =
0063     new Asymmetry("asy2fi",
0064                   "[0]*exp([1]*x)+[2]",
0065                   qsq, false);
0066
0067   Asyms->at(3) =
0068     new Asymmetry("asyfpfs",
0069                   "[0]*exp([1]*x)+[2]",
0070                   qsq, false);
0071
0072   Asyms->at(4) =
0073     new Asymmetry("asy3f",
0074                   "[0]*exp([1]*x)+[2]",
0075                   qsq, false);
0076
0077 }
0078
0079 double SigmaCalc::sigma_l()
0080 {
0081   if(*VertEvent->t_GeV < -1.5)
0082     return 1e-6;
0083   else{
0084     return MySigmaL(*VertEvent->qsq_GeV,
0085                     -*VertEvent->t_GeV,
0086                     *VertEvent->w_GeV);
0087   }
0088   return MySigmaL(*VertEvent->qsq_GeV,
0089                   -*VertEvent->t_GeV,
0090                   *VertEvent->w_GeV);
0091 }
0092
0093 double SigmaCalc::sigma_t()
0094 {
0095   if (*VertEvent->t_GeV < -1.5)
0096     return 1e-6;
0097   else{
0098     return MySigmaT(*VertEvent->qsq_GeV,
0099                     -*VertEvent->t_GeV,
0100                     *VertEvent->w_GeV);
0101   }
0102   return MySigmaT(*VertEvent->qsq_GeV,
0103                   -*VertEvent->t_GeV,
0104                   *VertEvent->w_GeV);
0105 }
0106
0107 double SigmaCalc::sigma_tt()
0108 {
0109   if (*VertEvent->t_GeV < -1.5)
0110     return 1e-7;
0111   else{
0112     double sig_tt = MySigmaTT(*VertEvent->qsq_GeV,
0113                               -*VertEvent->t_GeV,
0114                               *VertEvent->w_GeV);
0115     if ((sig_tt<0) && (Abs(sig_tt)>abs(this->sigma_t())))
0116       sig_tt = correctionToSigTT(sig_tt,
0117                                  this->sigma_t(),
0118                                  *VertEvent->qsq_GeV);
0119     return sig_tt;
0120   }
0121   double sig_tt = MySigmaTT(*VertEvent->qsq_GeV,
0122                             -*VertEvent->t_GeV,
0123                             *VertEvent->w_GeV);
0124   if ((sig_tt<0) && (Abs(sig_tt)>abs(this->sigma_t())))
0125     sig_tt = correctionToSigTT(sig_tt,
0126                                this->sigma_t(),
0127                                *VertEvent->qsq_GeV);
0128   return sig_tt;
0129 }
0130
0131 double SigmaCalc::sigma_lt()
0132 {
0133   if (*VertEvent->t_GeV < -1.5)
0134     return 1e-7;
0135   else{
0136     double sig_lt = MySigmaLT(*VertEvent->qsq_GeV,
0137                               -*VertEvent->t_GeV,
0138                               *VertEvent->w_GeV);
0139     if ((sig_lt<0) && (Abs(sig_lt)>abs(this->sigma_t())))
0140       sig_lt = correctionToSigLT(sig_lt,
0141                                  this->sigma_t(),
0142                                  *VertEvent->qsq_GeV);
0143     return sig_lt;
0144   }
0145   double sig_lt = MySigmaLT(*VertEvent->qsq_GeV,
0146                             -*VertEvent->t_GeV,
0147                             *VertEvent->w_GeV);
0148   if ((sig_lt<0) && (Abs(sig_lt)>abs(this->sigma_t())))
0149     sig_lt = correctionToSigLT(sig_lt,
0150                                this->sigma_t(),
0151                                *VertEvent->qsq_GeV);
0152   return sig_lt;
0153 }
0154
0155 double SigmaCalc::epsilon()
0156 {
0157   double q = RestEvent->VirtPhot->P()/1000;
0158   double theta = RestEvent->ScatElec->Theta();
0159
0160   return 1.0/(1.0 + 2.0*(q*q)/(*VertEvent->qsq_GeV)*
0161               Power(Tan(theta/2.0),2));
0162
0163 }
0164
0165 double SigmaCalc::sigma_uu()
0166 {
0167
0168   double phi = *TConEvent->Phi;
0169   double eps = this->epsilon();
0170   double siguu = this->sigma_t();
0171   //cout << siguu << endl;
0172   siguu += eps*this->sigma_l();
0173   siguu += Sqrt(2*eps*(1+eps))*this->sigma_lt()*Cos(phi);
0174   siguu += eps*this->sigma_tt()*Cos(2*phi);
0175   //siguu /= 2*Pi();
0176   return siguu;
0177
0178   //double siguu = ( this->sigma_t() + eps * this->sigma_l() +
0179   //               eps * TMath::Cos( 2.0 * phi ) * this->sigma_tt() +
0180   //              TMath::Sqrt( 0.5 * eps * ( 1.0 + eps ) ) * 2.0 * TMath::Cos( phi ) * this->sigma_lt() );
0181   //return siguu;
0182 }
0183
0184 double SigmaCalc::Sigma_k(int k)
0185 {
0186   double sigk = this -> Asyms->at(k)
0187     ->GetAsyAmp(*VertEvent->qsq_GeV,
0188                 *VertEvent->t_prime_GeV);
0189   // cout << sigk << endl;
0190   return sigk;
0191 }
0192
0193 double SigmaCalc::sigma_ut()
0194 {
0195   double phi = *TConEvent->Phi;
0196   double phi_s = *TConEvent->Phi_s;
0197   double theta = *RestEvent->Theta;
0198   double pt = *TConEvent->P_T;
0199
0200   pt = 0.865; // This is an approximations used in Ahmed's code.
0201   //It assumes the the q vector is parallel to the lab z axis.
0202
0203   double sigut = Sin(phi - phi_s)*this->Sigma_k(0);
0204   //cout << "sigma_ut1\t"<<sigut << endl;
0205   sigut += Sin(phi_s)*this->Sigma_k(1);
0206   sigut += Sin(2*phi-phi_s)*this->Sigma_k(2);
0207   sigut += Sin(phi+phi_s)*this->Sigma_k(3);
0208   sigut += Sin(3*phi-phi_s)*this->Sigma_k(4);
0209
0210   //cout << "sigma_ut2\t"<< sigut << endl;
0211
0212
0213   sigut *= pt/(Sqrt(1-Power(Sin(theta)*Sin(phi_s),2)));
0214   //cout << "sigma_ut3\t"<< sigut << endl;
0215   sigut *= this->sigma_uu();
0216
0217   //cout <<  "sigma_u4\t"<<sigut << endl;
0218
0219
0220   return sigut;
0221 }
0222
0223 double SigmaCalc::sigma()
0224 {
0225   double sigma = (this->sigma_uu()+sigma_ut());
0226   sigma *= this->fluxfactor_col();
0227   sigma *= this->jacobian_cm() * CofMEvent->ProdPion->P() / (1000*Pi());
0228   //sigma *= this->jacobian_A();
0229   sigma *= this->jacobian_cm_col();
0230   return sigma;
0231 }
0232
0233 double SigmaCalc::weight(int nGen)
0234 {
0235   return (this->sigma())*PSF()*uBcm2*Lumi/nGen;
0236 }
0237
0238 double SigmaCalc::fluxfactor_col()
0239 {
0240   double ff = alpha/(2*Power(Pi(),2));
0241   ff *= (VertEvent->ScatElec->E())/(VertEvent->BeamElec->E());
0242   ff *= (Power((*VertEvent->w_GeV),2))-Power(neutron_mass_gev,2);
0243   ff /= 2*neutron_mass_gev;
0244   ff /= *(VertEvent->qsq_GeV);
0245   ff /= (1-this->epsilon());
0246   return ff;
0247 }
0248
0249
0250 double SigmaCalc::jacobian_cm()
0251 {
0252   double beta = VertEvent->VirtPhot->P()/(VertEvent->VirtPhot->E()+neutron_mass_mev);
0253   double gamma = (VertEvent->VirtPhot->E() + neutron_mass_mev)/(*VertEvent->w_GeV * 1000);
0254   double jcm = RestEvent->VirtPhot->P()/1000;
0255   jcm -= beta * RestEvent->VirtPhot->E()/1000;
0256   jcm /= gamma * (1-Power(beta,2));
0257   return jcm;
0258 }
0259
0260 double SigmaCalc::jacobian_cm_col()
0261 {
0262   double jcol = Power(VertEvent->ProdPion->P(),2)*(*VertEvent->w_GeV*1000);
0263   jcol /= (CofMEvent->ProdPion->P()
0264            * Abs((neutron_mass_mev+VertEvent->VirtPhot->E())
0265                  * VertEvent->ProdPion->P()
0266                  - (VertEvent->ProdPion->E()*VertEvent->VirtPhot->P()
0267                     * Cos(VertEvent->ProdPion->Theta())
0268                     )
0269                  )
0270            );
0271   return jcol;
0272 }
0273
0274 double SigmaCalc::PSF()
0275 {
0276   double psf = (obj["scat_elec_Emax"].asDouble()-obj["scat_elec_Emin"].asDouble())/1000;
0277   psf *= (Sin(obj["scat_elec_thetamax"].asDouble()/DEG)
0278           -Sin(obj["scat_elec_thetamin"].asDouble()/DEG));
0279   psf *= (Sin(obj["prod_pion_thetamax"].asDouble()/DEG)
0280           -Sin(obj["prod_pion_thetamin"].asDouble()/DEG));
0281   psf *= 4*Pi()*Pi();
0282
0283   return psf;
0284 }
0285
0286 double SigmaCalc::jacobian_A(){
0287   //return CofMEvent->VirtPhot->Vect().Dot(CofMEvent->ProdPion->Vect())/(1000000*Pi());
0288   return (this->jacobian_cm() * CofMEvent->ProdPion->P() / (1000*Pi()));
0289 }