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