ExhumeInterface/src/TwoSpace.cpp

0001 #include "GeneratorInterface/ExhumeInterface/interface/TwoSpace.h"
0002 #include "CLHEP/Random/RandomEngine.h"
0003
0004 Exhume::TwoSpace::TwoSpace(const edm::ParameterSet& pset) : CrossSection(pset) {
0005   Partons.resize(2);
0006
0007   dirty_weighting = true;
0008   PhiMax = 0.0;
0009   PartonMass = 0.0;  //default, override in constructor of derived class
0010   MaximumSubProcessValue = 0.0;
0011   SetMassAtThetaScan(2000);
0012   Fudge = 1.2;
0013 }
0014
0015 Exhume::TwoSpace::~TwoSpace() {}
0016
0017 int Exhume::TwoSpace::GetNumberOfSubParameters() { return (2); }
0018
0019 void Exhume::TwoSpace::SetSubParameters() {
0020   double cos_th = randomEngine->flat();
0021
0022   CosTheta = GetValue(cos_th);
0023   Phi = 2 * M_PI * randomEngine->flat();
0024
0025   //std::cout<<CosTheta<<std::endl;
0026   LIPS2Amp();
0027 }
0028
0029 void Exhume::TwoSpace::MaximiseSubParameters() {
0030   if (dirty_weighting == true) {
0031     //  std::cout<<CosThetaMax<<std::endl;
0032     WeightInit(CosThetaMin, CosThetaMax);
0033     dirty_weighting = false;
0034   }
0035   CosTheta = MaximumSubProcessCosTheta;
0036   //  std::cout<<CosTheta<<std::endl;
0037   Phi = PhiMax;
0038   LIPS2Amp();
0039   return;
0040 }
0041
0042 void Exhume::TwoSpace::SetPartons() {
0043   Amp2LIPS();  //set lips parameters, not needed?
0044   double _SinTheta = sin(acos(CosTheta));
0045   double _E = 0.5 * SqrtsHat;
0046   double _Pmod = sqrt(_E * _E - PartonMass * PartonMass);
0047   double _Px = _Pmod * _SinTheta * cos(Phi);
0048   double _Py = _Pmod * _SinTheta * sin(Phi);
0049   double _Pz = _Pmod * CosTheta;
0050
0051   Partons[0].p.setPx(_Px);
0052   Partons[0].p.setPy(_Py);
0053   Partons[0].p.setPz(_Pz);
0054   Partons[0].p.setE(_E);
0055
0056   Partons[1].p.setPx(-_Px);
0057   Partons[1].p.setPy(-_Py);
0058   Partons[1].p.setPz(-_Pz);
0059   Partons[1].p.setE(_E);
0060
0061   Partons[0].p.boost(CentralVector.boostVector());
0062   Partons[1].p.boost(CentralVector.boostVector());
0063
0064   return;
0065 }
0066
0067 double Exhume::TwoSpace::SubParameterWeight() {
0068   //bit of a misnomer SubParameterWeight tells the Event class the difference between the subprocess Maximum and this specific chance
0069   Amp2LIPS();  //not needed?
0070   double _SubWgt;
0071   _SubWgt = GetFunc(CosTheta);
0072
0073   // std::cout<<_SubWgt/MaximumSubProcessValue<<std::endl;
0074   return (Fudge * _SubWgt / (MaximumSubProcessValue));
0075 }
0076
0077 double Exhume::TwoSpace::SubParameterRange() { return (Fudge * 2 * M_PI * TotalIntegral / MaximumSubProcessValue); }
0078 //
0079 //double Exhume::TwoSpace::SubProcess(){
0080 //return(AmplitudeSq*TwoParticleLIPS*Amp2CrossSection);
0081 //define overall factor in font of Msq + DLIPS factor
0082 //}
0083 double Exhume::TwoSpace::WeightFunc(const double& _CosTheta) {
0084   //arbitrary weighting at low mass
0085   //std::cout<<_CosTheta<<std::endl;
0086   CosTheta = _CosTheta;
0087   //std::cout<<CosTheta<<std::endl;
0088   Phi = PhiMax;  //arbitrary angle but amp might use it but shouldn'
0089
0090   //double wgt_hold=0.0;
0091   //for(int ii=0;ii<100;ii++){
0092   //double _SqrtsHat = (MassAtThetaScanHigh - MassAtThetaScanLow)*double(ii)/10
0093   // + MassAtThetaScanLow;
0094   //double _SqrtsHat=sqrt(4*PartonMass*PartonMass +
0095   //        MassAtThetaScan*MassAtThetaScan);
0096   double _SqrtsHat;
0097   if (PartonMass > 2.0) {
0098     _SqrtsHat = 2 * PartonMass + 0.001;
0099   } else {
0100     _SqrtsHat = sqrt(4 * PartonMass * PartonMass + 10.0);
0101   }
0102
0103   SetKinematics(_SqrtsHat, 0.0, 0.0, 0.0, 0.0, 0.0);
0104   //std::cout<<"here"<<std::endl;
0105   LIPS2Amp();                  //
0106                                //std::cout<<"here"<<std::endl;
0107   double _wgt = SubProcess();  //AmplitudeSq();
0108   //if(_wgt > wgt_hold){
0109   //wgt_hold=wgt;
0110   //}
0111
0112   if (_wgt > MaximumSubProcessValue) {
0113     MaximumSubProcessCosTheta = _CosTheta;
0114     MaximumSubProcessValue = _wgt;
0115     std::cout << MaximumSubProcessCosTheta << std::endl;
0116     std::cout << MaximumSubProcessValue << std::endl;
0117   }
0118   //std::cout<<_wgt<<"\t"<<CosTheta<<std::endl;
0119   return (_wgt);
0120 }
0121
0122 void Exhume::TwoSpace::SetThetaMin(const double& ThetaMin_) {
0123   ThetaMin = ThetaMin_;
0124   ThetaMax = M_PI - ThetaMin_;
0125   // SinThetaMin = sin(ThetaMin);
0126   CosThetaMin = cos(ThetaMax);
0127   CosThetaMax = cos(ThetaMin);
0128   MaximumSubProcessCosTheta = CosThetaMax;
0129
0130   //std::cout<<"called\t"<<CosThetaMin<<std::endl;
0131
0132   return;
0133 }
0134
0135 //Params outside in event are M,y,t1,t2,phi1,phi2