DEMPgen/src/ProductGen.cxx

0001 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0002   This class uses GSL MathMore's MultiRootFinder to find the
0003   momenta of the produced particles. In this case, a proton,
0004   and a negative pion.
0005   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
0006
0007 #include <vector>
0008 #include <stdio.h>
0009 #include <iostream>
0010
0011 #include "ProductGen.hxx"
0012 #include "CustomRand.hxx"
0013 #include "Constants.hxx"
0014
0015 #include "TMath.h"
0016 #include "TRandom3.h"
0017
0018 #include "json/json.h"
0019
0020 using namespace std;
0021 using namespace constants;
0022
0023 ProductGen::ProductGen(Particle* inInteraction, Particle* inTarget):
0024   Interaction(inInteraction), Target(inTarget)
0025 {
0026
0027   extern Json::Value obj;
0028
0029   char AngleGenName[100] = "AngleGen";
0030   double dummy[2] = {0,1};
0031   double ThetaRange[2] = {obj["prod_pion_thetamin"].asDouble()*TMath::DegToRad(),
0032                           obj["prod_pion_thetamax"].asDouble()*TMath::DegToRad()};
0033   double PhiRange[2] = {0, 360*TMath::DegToRad()};
0034   AngleGen = new CustomRand(AngleGenName, dummy,
0035                             ThetaRange, PhiRange);
0036   CoinToss = new TRandom3();
0037
0038   W_in_val = W_in();
0039
0040   Proton = new Particle();
0041   Pion = new Particle();
0042   Final = new Particle();
0043
0044   Initial = new Particle();
0045   *Initial = *Interaction+*Target;
0046
0047   UnitVect = new TVector3(0,0,1);
0048   F = new TF1("F",
0049               "[6]-sqrt([7]**2+x**2)-sqrt([8]**2+([3]-[0]*x)**2+([4]-[1]*x)**2+([5]-[2]*x)**2)",
0050               0, 12000);
0051
0052 //  cout << "CCCC "<< inInteraction->Px() << "  " << inInteraction->Py() << "  " << inInteraction->Pz() << "  " << inInteraction->E() << "  " << inInteraction->GetMass() << endl;
0053
0054
0055 }
0056
0057 void ProductGen::SetInteraction(Particle * inInteraction)
0058 {
0059   Interaction = inInteraction;
0060   W_in_val = W_in();
0061 }
0062
0063 void ProductGen::SetTarget(Particle * inTarget)
0064 {
0065   Target = inTarget;
0066   W_in_val = W_in();
0067 }
0068
0069 double ProductGen::W_in()
0070 {
0071   return (*Interaction+*Target).Mag2();
0072 }
0073
0074 double ProductGen::W_out()
0075 {
0076   return (*Proton+*Pion).Mag2();
0077 }
0078
0079 double ProductGen::Qsq_in()
0080 {
0081   return - Interaction->Mag2();
0082 }
0083
0084 int ProductGen::Solve()
0085 {
0086   double theta = AngleGen->Theta();
0087   double phi = AngleGen->Phi();
0088
0089 //  theta = 0.282478;
0090 //  phi = 3.49651;
0091
0092 //  cout << " Theta Phi: "<< theta << "   " << phi << endl;
0093
0094   return this->Solve(theta, phi);
0095 }
0096
0097 int ProductGen::Solve(double theta, double phi)
0098 {
0099   *Initial = *Interaction+*Target;
0100
0101   W_in_val = W_in();
0102
0103   if (W_in_val<0){
0104     //cout << "W < 0 " << endl;
0105     return 1;
0106   }
0107
0108   //cout << proton_mass_mev << endl;
0109
0110   UnitVect->SetTheta(theta);
0111   UnitVect->SetPhi(phi);
0112   UnitVect->SetMag(1);
0113
0114   double pars[9];
0115   pars[0] = UnitVect->X();
0116   pars[1] = UnitVect->Y();
0117   pars[2] = UnitVect->Z();
0118   pars[3] = Initial->Px();
0119   pars[4] = Initial->Py();
0120   pars[5] = Initial->Pz();
0121   pars[6] = Initial->E();
0122   pars[7] = pion_mass_mev;
0123   pars[8] = proton_mass_mev;
0124
0125 //  cout << Interaction->Px() << "  " << Interaction->Py() << "  " << Interaction->Pz() << "  " << Interaction->E() << "  " << Interaction->GetMass() << endl;
0126 //
0127 //  cout << Target->Px() << "  " << Target->Py() << "  " << Target->Pz() << "  " << Target->E() << "  " << Target->GetMass() << endl;
0128 //
0129 //  cout << endl;
0130 //  cout << pars[0] << "  " << pars[1]  << "  " << pars[2] << endl;
0131 //  cout << pars[3] << "  " << pars[4]  << "  " << pars[5] << endl;
0132 //  cout << pars[6] << "  " << pars[7]  << "  " << pars[8] << endl;
0133
0134
0135   F->SetParameters(pars);
0136
0137   double P = F->GetX(0, 0, pars[6], 0.0001, 10000);
0138
0139   //std::cout << "Zero: " << F->Eval(P) << std::endl;
0140
0141
0142   Particle * Pion1 = new Particle(pion_mass_mev,
0143                                   P*pars[0],
0144                                   P*pars[1],
0145                                   P*pars[2]);
0146   *Pion = *Pion1;
0147
0148   Particle * Proton1 = new Particle();
0149   *Proton1 = *Initial-*Pion;
0150   *Proton = *Proton1;
0151
0152 //  cout << Pion->E() << "   " << Proton->E() << endl;
0153
0154   if (TMath::Abs(F->Eval(P)) > 1){
0155     delete Pion1;
0156     delete Proton1;
0157     return 1;
0158   }
0159 //  cout << "AAAAAAAA " <<  SolnCheck() <<  endl;
0160
0161   if (!SolnCheck()){
0162     delete Pion1;
0163     delete Proton1;
0164     return 1;
0165   }
0166
0167
0168   //Check for Second solution:
0169   double P2 = F->GetX(0, P+100, pars[6], 0.0001, 10000);
0170
0171   if (TMath::Abs(F->Eval(P2))> 1){
0172     //No second soln
0173     delete Pion1;
0174     delete Proton1;
0175     return 0;
0176   }
0177
0178   //Try second solution
0179   Particle * Pion2 = new Particle(pion_mass_mev,
0180                                   P*pars[0],
0181                                   P*pars[1],
0182                                   P*pars[2]);
0183   *Pion = *Pion2;
0184
0185   Particle * Proton2 = new Particle();
0186   *Proton2 = *Initial - * Pion;
0187   *Proton = *Proton2;
0188
0189   if (SolnCheck()){
0190     //Toss a coin
0191     if (CoinToss->Uniform(0,1)>0.5){
0192       delete Pion1;
0193       delete Pion2;
0194       delete Proton1;
0195       delete Proton2;
0196       return 0; // Keep second solution
0197     }
0198   }
0199
0200   //Either SolnCheck or coin toss failed
0201   //Revert to original solution
0202
0203   *Proton = *Proton1;
0204   *Pion = *Pion1;
0205
0206   delete Pion1;
0207   delete Pion2;
0208   delete Proton1;
0209   delete Proton2;
0210
0211 //  exit(0);
0212   return 0;
0213 }
0214
0215 bool ProductGen::SolnCheck()
0216 {
0217   // Double Checking for solution viability
0218   if (TMath::Abs(proton_mass_mev-Proton->M())>1){
0219     //cerr << "Mass Missmatch" << endl;
0220     //cerr << TMath::Abs(proton_mass_mev-Proton->M()) << endl;
0221     return false;
0222   }
0223   if (TMath::Abs(W_in()-W_out())>1){
0224     //cerr << "W Missmatch" << endl;
0225     //cerr << TMath::Abs(W_in()-W_out()) << endl;
0226     return false;
0227   }
0228   *Final = *Proton + *Pion;
0229
0230   if (TMath::Abs(Initial->Px()-Final->Px())>1){
0231     //cerr << "Px Missmatch" << endl;
0232     //cerr << TMath::Abs(Initial->Px()-Final->Px()) << endl;
0233     return false;
0234   }
0235
0236   if (TMath::Abs(Initial->Py()-Final->Py())>1){
0237     //cerr << "Py Missmatch" << endl;
0238     //cerr << TMath::Abs(Initial->Py()-Final->Py()) << endl;
0239     return false;
0240   }
0241
0242   if (TMath::Abs(Initial->Pz()-Final->Pz())>1){
0243     //cerr << "Pz Missmatch" << endl;
0244     //cerr << TMath::Abs(Initial->Pz()-Final->Pz()) << endl;
0245     return false;
0246   }
0247
0248   if (TMath::Abs(Initial->E()-Final->E())>1){
0249     return false;
0250   }
0251   return true;
0252 }
0253
0254 Particle * ProductGen::ProdPion()
0255 {
0256   return Pion;
0257 }
0258
0259 Particle * ProductGen::ProdProton()
0260 {
0261   return Proton;
0262 }