Merge pull request #12 from dglazier/detach

Detach

README.md

@@ -3,7 +3,7 @@
 Event Generator framework for incorporating Spectroscopy into electro / photoproduction reactions.
 
 # Prerequisites
-
+ 
 ROOT 6
 
 with MathMore RooFit GenVector EG

core/DecayModelQ2W.cpp

@@ -18,15 +18,6 @@ namespace elSpectro{
   }
   ///////////////////////////////////////////////////////
   ///constructor includes subseqent decay of Ngamma* system
-  /* DecayModelQ2W::DecayModelQ2W( double thresh, DecayModel* gNmodel) :
-    _threshold{thresh},
-    DecayModel{{ new DecayingParticle{-2211,primary} },{11}}
-  {
-    _name={"DecayModelQ2W_with_primary_decay"};
-    Init();
-    }*/
-  ///////////////////////////////////////////////////////
-  ///constructor includes subseqent decay of Ngamma* system
   DecayModelQ2W::DecayModelQ2W( double thresh,
 				DecayModel* gNmodel,DecayVectors* gNdecayer) :
     _threshold{thresh},
@@ -115,7 +106,6 @@ namespace elSpectro{
     auto delta = 2*escat::M2_el()/getQ2()*(1-epsilon);
 
     _photonPol.SetEpsilon(epsilon);
-    //_photonPol.SetEpsilon(1);
     _photonPol.SetDelta(delta);
 
     //Get envelope weight from integrated cross section
@@ -193,14 +183,7 @@ namespace elSpectro{
       std::cout<<"DecayModelQ2W::FindExcitationSpectra()  result   "<<hist.GetMaximum()<<" "<<hist.GetBinCenter(hist.GetMaximumBin())<<" "<<hist.GetNbinsX()<<std::endl;
       hist.SetName("Wdist");
 
-      /*
-      TFile ff("debug.root","recreate");
-      hist.Write();
-      histpeak.Write();
-      histlow.Write();
-      ff.Close();
-      */
-
+
      _Wrealphoto_Dist.reset( new DistTH1(hist) );
     }
     else{

core/DecayModelst.cpp

@@ -125,25 +125,14 @@ namespace elSpectro{
       {
 	_s = x[0]*x[0];
 	_W=x[0];
-	//if( _W < Wmin ) return 1E2;
-	//if( _W < M3+M4 ) return 1E2;
 	if( _W < Wmin ) return 0.;
 	if( _W < M3+M4 ) return 0.;
-	//	if( x[1] > kine::t0(_W,M1,M2,M3,M4) ) return 0.;
-	//	if( x[1]< kine::tmax(_W,M1,M2,M3,M4) ) return 0.;
 
 	auto myt0=kine::t0(_W,M1,M2,M3,M4);
 	auto mytmax=kine::tmax(_W,M1,M2,M3,M4);
 
 	auto currt=x[1];
 
-	/*if( TMath::IsNaN(val) ) {
-	  if(TMath::Abs(x[1]-myt0)<TMath::Abs(x[1]-myt0) )
-	    curr_t = 0;
-	  else x[1] = mytmax- myt0;
-	  }*/
-
-	//	if( TMath::IsNaN(x[1]) ) return 1E2;
 	if( TMath::IsNaN(x[1]) ) return 0.;
 
 	//make smooth function at boundary, but must be lower than max
@@ -152,23 +141,16 @@ namespace elSpectro{
 	if( x[1] >  myt0){
 	  _t=myt0;
 	  factor*=(1-1*(x[1]-myt0)); //lower result
-	  //return 0.;
-	  // std::cout<<"too low "<<x[1]<<" "<<_t<<" "<<myt0<<std::endl;
 	}
 	else if( x[1]< mytmax ){
 	  _t=mytmax;
 	  factor*=(1-1*(mytmax-x[1])); //lower result
-	  // std::cout<<"too high "<<x[1]<<" "<<_t<<" "<<mytmax<<std::endl;
-	  //return 0.;
 	}
 	else _t=x[1];
 	if(factor<0) factor = 0;
 
 	double val = DifferentialXSect()* (myt0-mytmax)*factor;
-	//std::cout<<factor<<" "<<val<<" "<<(myt0-mytmax)<<"t "<<_t<<" "<<myt0<<" "<<mytmax<<" W "<<_W<<" x1 "<<x[1]<<"         "<<DifferentialXSect()<<std::endl;
-      	//if(val==0)return 1E2;
 	if( TMath::IsNaN(val) ) return 0.;
-	//if( TMath::IsNaN(val) ) return 1E2;
 	return -(val); //using a minimiser!
       };
 

core/DecayingParticle.cpp

@@ -1,6 +1,7 @@
 #include "DecayingParticle.h"
 #include "TwoBodyFlat.h"
 #include "Manager.h"
+#include "Interface.h"
 #include <TRandom.h>
 #include <TDatabasePDG.h>
 
@@ -25,10 +26,13 @@ namespace elSpectro{
   void DecayingParticle::PostInit(ReactionInfo* info) {
     //Decay type depends on Lifetime
     if(TDatabasePDG::Instance()->GetParticle(Pdg())){
-      double meanFreePath=TDatabasePDG::Instance()->GetParticle(Pdg())->Lifetime();
-      meanFreePath*=TMath::C()*100; //in mm
-      if( meanFreePath>0.01 ){ //10umm
+      double lifetime=TDatabasePDG::Instance()->GetParticle(Pdg())->Lifetime();
+      double meanFreePath=lifetime*TMath::C()*1000; //in mm
+      if( meanFreePath>0.1 ){ //0.1mm
     	_decayType=DecayType::Detached;
+	_decVertexDist = new DistTF1(TF1("MFP","TMath::Exp(-x/[0])",0,25*lifetime));//in mm
+	_decVertexDist->GetTF1().SetParameter(0,lifetime);
+	_decVertexDist->GetTF1().SetNpx(500);
       }
       else _decayType=DecayType::Attached;
     }
@@ -141,5 +145,22 @@ namespace elSpectro{
 
   }
 
-
+   void DecayingParticle::GenerateVertexPosition()  noexcept{
+      //auto old=VertexPosition();
+      if( IsDecay()==DecayType::Detached){
+	Double_t t0=_decVertexDist->SampleSingle();//in s
+	//Need lab 4-vector
+	LorentzVector lab=P4();
+	generator().BoostToLab(lab);
+
+	Double_t r= t0 * lab.Gamma() * TMath::C() * lab.Beta() *1000; //Lorentz contraction , mm
+	Double_t labP=lab.P();
+	//Set in direction of particle momentum
+	//with length of decay
+	_decayVertex.SetXYZT(lab.X()/labP*r,lab.Y()/labP*r,lab.Z()/labP*r,r/1000/TMath::C());
+	//add production vertex
+	_decayVertex+=*VertexPosition(); 
+      }
+    }
+
 }

core/DecayingParticle.h

@@ -14,6 +14,7 @@
 #include "DecayVectors.h"
 #include "TwoBodyFlat.h"
 #include "ReactionInfo.h"
+#include "DistTF1.h"
 
 namespace elSpectro{
 
@@ -76,13 +77,8 @@ namespace elSpectro{
     virtual void PostInit(ReactionInfo* info);
 
     //temporary until deal with vertices properly i.e. non zero
-    virtual void GenerateVertexPosition()  noexcept{
-      //auto old=VertexPosition();
-      if( IsDecay()==DecayType::Detached){
-
-      }
-      // _decayVertex.SetXYZT(old.X(),old.Y(),old.Z(),old.T());
-    }
+    virtual void GenerateVertexPosition()  noexcept;
+
     const LorentzVector& DecayVertexPosition()const noexcept{return _decayVertex;}
     int DecayVertexID()const noexcept{return _decayVertexID;}
 
@@ -106,7 +102,8 @@ namespace elSpectro{
     DecayModel* _decay={nullptr}; //not owner
 
     std::unique_ptr<DecayVectors> _decayer={nullptr}; //owner
-
+    DistTF1* _decVertexDist=nullptr;//! needed if detached vertex
+
     double _minMass={0};
     LorentzVector _decayVertex;
     int _decayVertexID={0};

core/DistTF1.h

@@ -45,7 +45,7 @@ namespace elSpectro{
     double GetWeightFor(double valX) const {return _tf1.Eval(valX)/_max_val;}
     double GetValueFor(double valX,double valY=0) final {return _tf1.Eval(valX);}
 
-    const TF1& GetTF1() const noexcept {return _tf1;}
+    TF1& GetTF1()  noexcept {return _tf1;}
 
   private:
     //no one should use default constructor

core/DistVirtPhotFlux_xy.cpp

@@ -16,7 +16,7 @@ namespace elSpectro{
     _ebeam(eb),
     _mTar(mion)
    {
-     // SetWThreshold(Wmin);
+
    }
 
   void DistVirtPhotFlux_xy::SetWThreshold(double Wmin){
@@ -94,11 +94,6 @@ namespace elSpectro{
 
     std::cout<<"DistVirtPhotFlux_xy max "<<_max_val<<" within y limits "<<TMath::Exp(_lnymin)<<" "<<TMath::Exp(_lnymax)<<" minimum possible x "<<_maxPossiblexRange<<std::endl;
     std::cout<<"  Other limits : "<<std::endl;
-
-    // if(_maxPossiblexRange==0)
-    // _maxPossiblexRange=200;
-    //else
-    //  _maxPossiblexRange=-TMath::Log(_maxPossiblexRange);
 
     if(_requestQ2min!=0)std::cout<<"\tQ2min "<<_requestQ2min<<std::endl;
     if(_requestQ2max!=0)std::cout<<"\tQ2max "<<_requestQ2max<<std::endl;
@@ -117,8 +112,6 @@ namespace elSpectro{
     //Finally, Integrate over photon flux
     auto xvar = RooRealVar("x","x",-(_lnxmax-_lnxmin)/2,_lnxmin,_lnxmax,"");
     auto yvar = RooRealVar("y","y",-(_lnymax-_lnymin)/2,_lnymin,_lnymax,"");
-    // auto xvar = RooRealVar("x","x",TMath::Exp((_lnxmax-_lnxmin)/2),TMath::Exp(_lnxmin),TMath::Exp(_lnxmax),"");
-    //auto yvar = RooRealVar("y","y",TMath::Exp((_lnymax-_lnymin)/2),TMath::Exp(_lnymin),TMath::Exp(_lnymax),"");
     xvar.Print();
     yvar.Print();
 
@@ -179,7 +172,7 @@ namespace elSpectro{
 	    (_val=escat::flux_dlnxdlny(_ebeam,lnx,lny)) ) {
       if(_val>_max_val){
 	_max_val=_val;
-	std::cout<<" MAX REACHED "<<_val<<" "<<_max_val<<std::endl;
+	std::cout<<"DistVirtPhotFlux_xy::FindWithAcceptReject() MAX REACHED "<<_val<<" "<<_max_val<<std::endl;
 	exit(0);
       }
       getRandomXY(); //sample another pair
@@ -198,83 +191,5 @@ namespace elSpectro{
 
   }
 
-  /*
-   void DistVirtPhotFlux_xy::FindFlat(){
-
-     _val=0;
-     _xy=std::make_pair(-1,-1); //unphysical should never be used
-     double y = gRandom->Uniform( TMath::Exp(_lnymin),TMath::Exp(_lnymax) );
-     //double y = gRandom->Uniform( 1E-8,1 );
-     //if(y<0.1) return; 
-     //calculate the fraction of x-space available
-     //now calculate x limits
-     //y = r/2ME and x = Q2/r 
-     double avail_xmax = XMax(y);
-     double avail_xmin = XMin(y);
-  
-     //if(avail_xmin>=1){ return; }
-     //_maxPossiblexRange=0;
-     double x=-1;
-     if(_maxPossiblexRange)
-       //for efficiency we need not sample below lowest possible x value
-       x = gRandom->Uniform(_maxPossiblexRange,1);
-     else 
-       x = gRandom->Uniform(0,1);
-     
-     //check if we are within allowed x-range
-     if( x<avail_xmax && x>avail_xmin )
-       _val=escat::flux_dxdy(_ebeam,x,y);
-     
-     //return x and y values
-     _xy=std::make_pair(x,y);
-     
-       
-     return;
-     
-   }
-  */
-  /*
-   void DistVirtPhotFlux_xy::FindFlat(){
-
-     TF1 sampleX("sampleX","TMath::Exp(-x[0])",0,1);
-     sampleX.SetRange(0,1);
-     sampleX.SetNpx(1E4);
-     TF1 sampleY("sampleY","TMath::Exp(-x[0])",0,1);
-     sampleY.SetRange(TMath::Exp(_lnymin),TMath::Exp(_lnymax));
-     sampleY.SetNpx(1E4);
-    
-     _val=0;
-     _xy=std::make_pair(-1,-1); //unphysical should never be used
-    
-      double y = sampleY.GetRandom();
-      double sampleYVal=sampleY.Eval(y)/sampleYIntegral;
-     
-     //calculate the fraction of x-space available
-     //now calculate x limits
-     //y = r/2ME and x = Q2/r 
-     double avail_xmax = XMax(y);
-     double avail_xmin = XMin(y);
-  
-     if(avail_xmin>=1){ return; }
-
-     double x=-1;
-     while(x<_maxPossiblexRange)
-        x = sampleX.GetRandom();
-     
-     double sampleXVal=sampleX.Eval(x)/sampleXIntegral;  //value of PDF
-    
-     //check if we are within allowed x-range
-     if( x<avail_xmax && x>avail_xmin )
-       _val=escat::flux_dxdy(_ebeam,x,y)/sampleYVal/sampleXVal;
-     
-     //return x and y values
-     _xy=std::make_pair(x,y);
-     
-     
-      
-     return;
-     
-   }
 
-  */
 }//namespace

core/DistVirtPhotFlux_xy.h

@@ -19,15 +19,13 @@ namespace elSpectro{
 
   public :
 
-    //DistVirtPhotFlux_xy(double ebeam,float xmin,float xmax,float ymin,float ymax);
     DistVirtPhotFlux_xy(double eb, double mion, double Wmin);
 
     double SampleSingle()   noexcept final {
       return 0;
     }
 
     dist_pair SamplePair()   noexcept final {
-      //    _forIntegral==false ? FindWithAcceptReject() : FindFlat();
       FindWithAcceptReject();
       return _xy;
     }
@@ -48,7 +46,6 @@ namespace elSpectro{
     double GetWMin() const noexcept {return TMath::Sqrt(_Wthresh2);}
 
     void FindWithAcceptReject();
-    // void FindFlat();
 
     void SetElecE(double ee){_ebeam=ee;}
     void SetM(double m){_mTar=m;}

core/ElectronScattering.cpp

@@ -306,6 +306,7 @@ namespace elSpectro{
     //Boost into ion rest frame
     auto prBoost=_beamNucl.P4().BoostToCM();
     collision=boost(collision,prBoost);
+    SetBoostToLab(-prBoost);
     _nuclRestNucl=LorentzVector(0,0,0,_beamNucl.Mass());
     _nuclRestElec= boost(_beamElec.P4(),prBoost);
     //set decay parent for e -> e'g*

core/LundWriter.h

@@ -54,10 +54,11 @@ namespace elSpectro{
        auto p4=p->P4();
        auto ver=p->VertexPosition();
        //second entry 0. == lifetime Could add to Particle.h
+       //note convert vertex mm->cm
        _stream<<_id++<<" "<<0.<<" "<<status
 	      <<" "<<p->Pdg()<<" "<<0<<" "<<0<<" "
 	      <<p4.X()<<" "<<p4.Y()<<" "<<p4.Z()<<" "<<p4.T()<<" "
-	      <<p4.M()<<" "<<ver->X()<<" "<<ver->Y()<<" "<<ver->Z()<<"\n";
+	      <<p4.M()<<" "<<ver->X()/10<<" "<<ver->Y()/10<<" "<<ver->Z()/10<<"\n";
      }
 
      //data members

core/Manager.h

@@ -74,8 +74,10 @@ namespace elSpectro{
      void Reaction(ProductionProcess* prod){
        _process.reset(prod);
      }
-
-     ProductionProcess* Reaction(){return _process.get();}
+    void BoostToLab(LorentzVector& boostme){
+      _process->BoostToLab(boostme);
+    }
+    ProductionProcess* Reaction(){return _process.get();}
 
      void SetSeed(ULong_t seed = 0){gRandom->SetSeed(seed);}
 

core/ParticleManager.cpp

@@ -1,11 +1,14 @@
 #include "ParticleManager.h"
 #include <TDatabasePDG.h>
+#include <TSystem.h>
 
 namespace elSpectro{
 
   ParticleManager::ParticleManager(){
 
-    TDatabasePDG *pdgDB = TDatabasePDG::Instance();
+    TDatabasePDG *pdgDB = new TDatabasePDG();
+    pdgDB->ReadPDGTable(Form("%s/etc/el_pdg_table.txt",gSystem->Getenv("ELSPECTRO")));
+
     //name,title,mass,stable,width,charge,type.code 
     pdgDB->AddParticle("gamma_star","gamma_star", 0.0, kFALSE,
 		       0, 0, "virtual", -22);
@@ -32,6 +35,8 @@ namespace elSpectro{
 
     pdgDB->AddParticle("deuteron","deuteron", 1.875612, kTRUE,0, 1, "Baryon", 45); //Jlab CLAS numbering
     pdgDB->AddParticle("deuteron","deuteron", 1.875612, kTRUE,0, 1, "Baryon", 1000010020); //PDG code numbering
+
+    //Lambda des not have lifetime in pdg_table.txt
 
   }
   Particle*  ParticleManager::Take(Particle* p){