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 }
0053
0054 void ProductGen::SetInteraction(Particle * inInteraction)
0055 {
0056   Interaction = inInteraction;
0057   W_in_val = W_in();
0058 }
0059
0060 void ProductGen::SetTarget(Particle * inTarget)
0061 {
0062   Target = inTarget;
0063   W_in_val = W_in();
0064 }
0065
0066 double ProductGen::W_in()
0067 {
0068   return (*Interaction+*Target).Mag2();
0069 }
0070
0071 double ProductGen::W_out()
0072 {
0073   return (*Proton+*Pion).Mag2();
0074 }
0075
0076 double ProductGen::Qsq_in()
0077 {
0078   return - Interaction->Mag2();
0079 }
0080
0081 int ProductGen::Solve()
0082 {
0083   double theta = AngleGen->Theta();
0084   double phi = AngleGen->Phi();
0085
0086   return this->Solve(theta, phi);
0087 }
0088
0089 int ProductGen::Solve(double theta, double phi)
0090 {
0091   *Initial = *Interaction+*Target;
0092
0093   W_in_val = W_in();
0094
0095   if (W_in_val<0){
0096     return 1;
0097   }
0098
0099   UnitVect->SetTheta(theta);
0100   UnitVect->SetPhi(phi);
0101   UnitVect->SetMag(1);
0102
0103   double pars[9];
0104   pars[0] = UnitVect->X();
0105   pars[1] = UnitVect->Y();
0106   pars[2] = UnitVect->Z();
0107   pars[3] = Initial->Px();
0108   pars[4] = Initial->Py();
0109   pars[5] = Initial->Pz();
0110   pars[6] = Initial->E();
0111   pars[7] = pion_mass_mev;
0112   pars[8] = proton_mass_mev;
0113
0114   F->SetParameters(pars);
0115
0116   double P = F->GetX(0, 0, pars[6], 0.0001, 10000);
0117
0118   Particle * Pion1 = new Particle(pion_mass_mev,
0119                                   P*pars[0],
0120                                   P*pars[1],
0121                                   P*pars[2]);
0122   *Pion = *Pion1;
0123
0124   Particle * Proton1 = new Particle();
0125   *Proton1 = *Initial-*Pion;
0126   *Proton = *Proton1;
0127
0128   if (TMath::Abs(F->Eval(P)) > 1){
0129     delete Pion1;
0130     delete Proton1;
0131     return 1;
0132   }
0133
0134   if (!SolnCheck()){
0135     delete Pion1;
0136     delete Proton1;
0137     return 1;
0138   }
0139
0140   //Check for Second solution:
0141   double P2 = F->GetX(0, P+100, pars[6], 0.0001, 10000);
0142
0143   if (TMath::Abs(F->Eval(P2))> 1){
0144     //No second soln
0145     delete Pion1;
0146     delete Proton1;
0147     return 0;
0148   }
0149
0150   //Try second solution
0151   Particle * Pion2 = new Particle(pion_mass_mev,
0152                                   P*pars[0],
0153                                   P*pars[1],
0154                                   P*pars[2]);
0155   *Pion = *Pion2;
0156
0157   Particle * Proton2 = new Particle();
0158   *Proton2 = *Initial - * Pion;
0159   *Proton = *Proton2;
0160
0161   if (SolnCheck()){
0162     //Toss a coin
0163     if (CoinToss->Uniform(0,1)>0.5){
0164       delete Pion1;
0165       delete Pion2;
0166       delete Proton1;
0167       delete Proton2;
0168       return 0; // Keep second solution
0169     }
0170   }
0171
0172   //Either SolnCheck or coin toss failed
0173   //Revert to original solution
0174
0175   *Proton = *Proton1;
0176   *Pion = *Pion1;
0177
0178   delete Pion1;
0179   delete Pion2;
0180   delete Proton1;
0181   delete Proton2;
0182
0183 //  exit(0);
0184   return 0;
0185 }
0186
0187 bool ProductGen::SolnCheck()
0188 {
0189   // Double Checking for solution viability
0190   if (TMath::Abs(proton_mass_mev-Proton->M())>1){
0191     //cerr << "Mass Missmatch" << endl;
0192     //cerr << TMath::Abs(proton_mass_mev-Proton->M()) << endl;
0193     return false;
0194   }
0195   if (TMath::Abs(W_in()-W_out())>1){
0196     //cerr << "W Missmatch" << endl;
0197     //cerr << TMath::Abs(W_in()-W_out()) << endl;
0198     return false;
0199   }
0200   *Final = *Proton + *Pion;
0201
0202   if (TMath::Abs(Initial->Px()-Final->Px())>1){
0203     //cerr << "Px Missmatch" << endl;
0204     //cerr << TMath::Abs(Initial->Px()-Final->Px()) << endl;
0205     return false;
0206   }
0207
0208   if (TMath::Abs(Initial->Py()-Final->Py())>1){
0209     //cerr << "Py Missmatch" << endl;
0210     //cerr << TMath::Abs(Initial->Py()-Final->Py()) << endl;
0211     return false;
0212   }
0213
0214   if (TMath::Abs(Initial->Pz()-Final->Pz())>1){
0215     //cerr << "Pz Missmatch" << endl;
0216     //cerr << TMath::Abs(Initial->Pz()-Final->Pz()) << endl;
0217     return false;
0218   }
0219
0220   if (TMath::Abs(Initial->E()-Final->E())>1){
0221     return false;
0222   }
0223   return true;
0224 }
0225
0226 Particle * ProductGen::ProdPion()
0227 {
0228   return Pion;
0229 }
0230
0231 Particle * ProductGen::ProdProton()
0232 {
0233   return Proton;
0234 }