DEMPgen/src/FSI.cxx

0001 #include "FSI.hxx"
0002
0003 #include "TMath.h"
0004
0005 #include "CustomRand.hxx"
0006 #include "Particle.hxx"
0007 #include "Constants.hxx"
0008 #include "PhaseShift.h"
0009
0010 #include <stdio.h>
0011 #include "json/json.h"
0012
0013 using namespace std;
0014 using namespace constants;
0015 using namespace TMath;
0016
0017 FSI::FSI()
0018 {
0019   double ERange[2] = {0,1};
0020   double ThetaRange[2] = {0,Pi()};
0021   double PhiRange[2] = {0,2*Pi()};
0022   char fname[100] = "FSI";
0023   SpherePicker = new CustomRand(fname, ERange, ThetaRange, PhiRange);
0024
0025   extern Json::Value obj;
0026   ProtGen = new TargetGen(proton_mass_mev, obj["fermi_momentum"].asBool());
0027
0028   CoP = new TVector3();
0029
0030   VertInPion = new Particle(pion_mass_mev, "", pid_pion);
0031   VertTargProt = new Particle(proton_mass_mev, "", pid_prot);
0032
0033   CMInPion = new Particle(pion_mass_mev, "", pid_pion);
0034   CMTargProt = new Particle(proton_mass_mev, "", pid_prot);
0035
0036   VertOutPion = new Particle(pion_mass_mev, "", pid_pion);
0037   VertOutProt = new Particle(proton_mass_mev, "", pid_prot);
0038
0039   CMOutPion = new Particle(pion_mass_mev, "", pid_pion);
0040   CMOutProt = new Particle(proton_mass_mev, "", pid_prot);
0041
0042   PhaseShiftWeight = new double(0);
0043
0044   WilliamsWeight = new double(0);
0045   DedrickWeight = new double(0);
0046   CatchenWeight = new double(0);
0047 }
0048
0049 int FSI::Generate()
0050 {
0051
0052   *VertTargProt = *ProtGen->GetParticle();
0053
0054   *CMInPion = *VertInPion;
0055   *CMTargProt = *VertTargProt;
0056
0057   //cout << "Pion Momentum Before:\t"<<CMInPion->P() << endl;
0058   //cout << "Proton Momentum Before:\t"<<CMTargProt->P()<<endl;
0059
0060   *CoP = ((VertInPion->Vect()+VertTargProt->Vect())*
0061           (1.0/(VertInPion->E()+VertTargProt->E())));
0062
0063
0064   CMInPion->Boost(-(*CoP));
0065   //cout << "Pion Momentum After:\t"<<CMInPion->P() << endl;
0066   CMTargProt->Boost(-(*CoP));
0067   //cout << "Proton Momentum After:\t"<<CMTargProt->P()<<endl;
0068
0069   theta_pion = SpherePicker->Theta();
0070   phi_pion = SpherePicker->Phi();
0071
0072   // Some factors to simplify the formula
0073   a = Sqrt(Power(pion_mass_mev, 2)+Power(CMInPion->P(), 2))
0074     + Sqrt(Power(proton_mass_mev, 2)+Power(CMTargProt->P(), 2));
0075   //cout<<a<<endl;
0076   b = Power(pion_mass_mev, 2);
0077   //cout<<b<<endl;
0078   c = Power(proton_mass_mev, 2);
0079   //cout<<c<<endl;
0080
0081   // Solution to E_i = E_f, taking advantage of the fact that the momenta of
0082   // outgoing particles are equal ond opposite in the CoM frame, leaving
0083   // pion momentum the only unknown.
0084   x = (a*a*a*a+b*b+c*c-2*a*a*b-2*a*a*c-2*b*c)/(4*a*a);
0085
0086   //cout<<x<<endl;
0087
0088   p_pion = Sqrt(x);
0089
0090   //cout<<p_pion<<endl;
0091
0092   CMOutPion->SetThetaPhiP(theta_pion, phi_pion, p_pion);
0093   //cout << "Pion final P:\t" << CMOutPion->P() << endl;
0094   CMOutProt->SetVectM(-(CMOutPion->Vect()), proton_mass_mev);
0095   //cout << "Prot final P:\t" << CMOutProt->P() << endl;
0096
0097   *VertOutPion = *CMOutPion;
0098   *VertOutProt = *CMOutProt;
0099
0100   VertOutProt->Boost((*CoP));
0101   VertOutPion->Boost((*CoP));
0102
0103  //cout << "Boost x:\t" << CoP->X() <<endl;
0104  //cout << "Prot x: \t" << CMOutProt->X()<<endl;
0105  //cout << "Boost y:\t" << CoP->Y() <<endl;
0106  //cout << "Prot y: \t" << CMOutProt->Y()<<endl;
0107  //cout << "Boost z:\t" << CoP->Z() <<endl;
0108  //cout << "Prot z: \t" << CMOutProt->Z()<<endl;
0109  //cout << "CM E:\t" << CMOutProt->E()<<endl;
0110  //cout << "Vert E:\t" << VertOutProt->E()<<endl;
0111
0112   //Check cons laws:
0113
0114   if (((*CMInPion+*CMTargProt)-(*CMOutProt+*CMOutPion)).Px() > 1.0){
0115     cout << "Px: " << ((*CMInPion+*CMTargProt)-(*CMOutProt+*CMOutPion)).Px() << endl;
0116     return 1;
0117   }
0118   if (((*CMInPion+*CMTargProt)-(*CMOutProt+*CMOutPion)).Py() > 1.0){
0119     cout << "Py: " << ((*CMInPion+*CMTargProt)-(*CMOutProt+*CMOutPion)).Py() << endl;
0120     return 1;
0121   }
0122   if (((*CMInPion+*CMTargProt)-(*CMOutProt+*CMOutPion)).Pz() > 1.0){
0123     cout << "Pz: " << ((*CMInPion+*CMTargProt)-(*CMOutProt+*CMOutPion)).Pz() << endl;
0124     return 1;
0125   }
0126   if (((*CMInPion+*CMTargProt)-(*CMOutProt+*CMOutPion)).E() > 1.0){
0127     cout << "E: " << ((*CMInPion+*CMTargProt)-(*CMOutProt+*CMOutPion)).E() << endl;
0128       return 1;
0129     }
0130     else
0131       return 0;
0132
0133 }
0134
0135 int FSI::CalculateWeights()
0136 {
0137   phaseshifts(2, CMOutPion->P()/1000, (*CMInPion+*CMTargProt).Mag2()/1000000);
0138
0139   //cout << CMOutPion->P() << endl;
0140   //cout << CMOutPion->Px() << endl;
0141
0142   //cout << "Check 1" << endl;
0143
0144   Z0 = getZ0();
0145   //cout << Z0 << endl;
0146   Z1 = getZ1();
0147   //cout << Z1 << endl;
0148   Z2 = getZ2();
0149   //cout << Z2 << endl;
0150
0151   //cout << "Check 2" << endl;
0152
0153   //cout << (Z0 +
0154   //        (Z1 * CMOutPion->Px()/CMOutPion->P()) +
0155   //        (Z2 * Power(CMOutPion->Px()/CMOutPion->P(), 2))
0156   //         ) << endl;
0157
0158   *PhaseShiftWeight = (Z0 +
0159                       (Z1 * CMOutPion->Px()/CMOutPion->P()) +
0160                       (Z2 * Power(CMOutPion->Px()/CMOutPion->P(), 2))
0161                       );
0162   //cout << *PhaseShiftWeight << endl;
0163   *PhaseShiftWeight *= 0.012; //.012 nucleons per half of He_3 nucleus area in milli barns
0164
0165   //cout << "Check 3" << endl;
0166
0167   beta = (VertInPion->P()+VertTargProt->P())/(VertInPion->E()+VertTargProt->E());
0168   gamma = (VertInPion->E()+VertTargProt->E()) / (*VertInPion+*VertTargProt).Mag();
0169
0170   //cout << "Check 4" << endl;
0171
0172   // Jacobian by W. S. C. Williams
0173   *WilliamsWeight = *PhaseShiftWeight * (VertInPion->Vect().Mag2()/
0174                                         (gamma*CMOutPion->P()*
0175                                          (VertInPion->P()-(beta*VertInPion->E()*
0176                                                            VertInPion->CosTheta())
0177                                           )
0178                                          )
0179                                         );
0180
0181   //cout << "Check 5" << endl;
0182
0183
0184   beta_pion = CMOutPion->P() / CMOutPion->E();
0185   g = beta / beta_pion;
0186
0187   //cout << "Check 6" << endl;
0188
0189
0190   // Jacobian by K. G. Dedrick. Rev. Mod. Phys. 34, 429 (1962)
0191   *DedrickWeight = *PhaseShiftWeight * ((Power
0192                                         (Power(g+CMOutPion->CosTheta(), 2)+
0193                                          (1 - beta * beta)*
0194                                          (1 - Power(CMOutPion->CosTheta(), 2)),
0195                                          1.5))/
0196                                        ((1-beta*beta)*Abs(1+g*CMOutPion->CosTheta())
0197                                         )
0198                                        );
0199
0200   //cout << "Check 8" << endl;
0201
0202
0203   // Jacobian by Gary L. Catchen J. Chem. Phys. 69(4), 15 Aug 1978
0204   *CatchenWeight = *PhaseShiftWeight * (Power(VertInPion->P(),2)*CMOutPion->E()/
0205                                        (Power(CMOutPion->P(),2)*VertInPion->E())
0206                                        );
0207
0208   //cout << "Check 9" << endl;
0209
0210
0211   return 0;
0212
0213 }
0214
0215 int FSI::GenerateNoRandom(){
0216   //Debug only
0217   //For this function, VertInPion and VertOutPion must be specified
0218   *VertTargProt = *ProtGen->GetParticle();
0219
0220   *CMInPion = *VertInPion;
0221   *CMTargProt = *VertTargProt;
0222
0223   *CoP = ((VertInPion->Vect()+VertTargProt->Vect())*
0224           (1.0/(VertInPion->E()+VertTargProt->E())));
0225
0226   CMInPion->Boost(-(*CoP));
0227   CMTargProt->Boost(-(*CoP));
0228
0229   *CMOutPion = *VertOutPion;
0230   CMOutPion->Boost(-(*CoP));
0231
0232   CMOutProt->SetVectM(-(CMOutPion->Vect()), proton_mass_mev);
0233
0234   *VertOutProt = *CMOutProt;
0235   VertOutProt->Boost(*CoP);
0236
0237   return 0;
0238
0239 }