Dubna1mHBC ZakladneSystemySuradnic

Základné systémy súradníc a generovanie pružného rozptylu [math]\displaystyle{ dp \rightarrow dp }[/math]

Základné systémy súradníc a kinematické veličiny: modul hybnosti a uhol výletu (polárny uhol), pre názornosť a porovnanie v sústavách laboratorná, ťažisková a antilaboratorná pre produkty reakcie

Pre analýzu zrážky dvoch častíc a a b sa najčastejsie používajú tri systemy súradníc: laboratorný (L), ťažiskový (^*) a antilaboratorný (AL). Súradnicový system vyberieme tak, aby počiatok bol v bode zrážky a os z bola totožná s nenulovou hybnosťou jednej z častíc v počiatočnom stave. Označme štvorhybnosti častíc v počiatočnom stave [math]\displaystyle{ P_{a}=(E_{a}, \bf{p}_{a}) }[/math] a [math]\displaystyle{ P_{b}=(E_{b}, \bf{p}_{b}) }[/math]. Tým pádom systemy sú určené hybnosťami [math]\displaystyle{ \bf{p}_{a} }[/math] alebo [math]\displaystyle{ \bf{p}_{b} }[/math]:

1. Laboratorný system (tiež kľudová sústava terčíka), sa používa pri meraní a prvotné experimentálne údaje sú nabrané v tomto systeme a odtiaľ sa transformujú do iných. Keď zväzkovú časticu označíme [math]\displaystyle{ a }[/math] a terčíkovú [math]\displaystyle{ b }[/math] potom ich štvorhybnosti sú nasledovné: [math]\displaystyle{ P_{a}^{(L)}=(E_{a}^{(L)}, 0, 0, p_{a}^{(L)}) }[/math] a [math]\displaystyle{ P_{b}^{(L)}=(m_{b}, 0, 0, 0) }[/math].
2. Ťažiskový system, v ktorom suma trojhybností je nulová: [math]\displaystyle{ {\bf{p}}_{a}^* + {\bf{p}}_{b}^* = 0 }[/math].
3. Antilaboratorný (tiež kľudová sústava zväzku), [math]\displaystyle{ P_{a}^{AL}=(m_{a}, 0, 0, 0) }[/math] a [math]\displaystyle{ P_{b}^{(AL)}=(E_{b}^{(AL)}, 0, 0, p_{b}^{(AL)}) }[/math], kde [math]\displaystyle{ E_{a}^{(L)}= \sqrt{m_{a}^2 + p_{a}^{(L)2}} }[/math] a [math]\displaystyle{ E_{b}^{(AL)}= \sqrt{m_{b}^2 + p_{b}^{(AL)2}} }[/math] energia častice [math]\displaystyle{ a }[/math] a [math]\displaystyle{ b }[/math] a [math]\displaystyle{ m_a }[/math] a [math]\displaystyle{ m_b }[/math] sú ich hmotnosti.
4. Sústava protibežných zväzkov (colliding beams) obe častice [math]\displaystyle{ a }[/math] a [math]\displaystyle{ b }[/math] sa zrážajú pod malým uhlom [math]\displaystyle{ \pi - \theta }[/math], kde [math]\displaystyle{ \theta }[/math] je uhol medzi ich hybnosťami. Túto sústavu nebudeme ďalej rozoberať.

Prechod z jednej sústavy do druhej sa uskutoční Lorentzovou špeciálnou transformáciou [math]\displaystyle{ E^\prime = \gamma(E - v p_z) }[/math], [math]\displaystyle{ p_{x}^\prime = p_{x} }[/math], [math]\displaystyle{ p_{y}^\prime = p_{y} }[/math], [math]\displaystyle{ p_{z}^\prime = \gamma(E - v p_z) }[/math], kde [math]\displaystyle{ {\bf{v}} = (0, 0, v) }[/math] je rýchlosť čiarkovanej sústavy voči pôvodnej, [math]\displaystyle{ \gamma = (1-v^2)^{-1/2} }[/math]:

• Rýchlosť tažiskovej sústavy (TS) voči laboratornej sústave [math]\displaystyle{ {\bf v}_{TS} = (0, 0, v_{TS}) }[/math] je: [math]\displaystyle{ v_{TS} = \frac{p_{a}^{(L)}}{E_a^{L} + m_{b}} }[/math].
• Rýchlosť systemu AL voči laboratornej sústave [math]\displaystyle{ {\bf v}_{AL} = (0, 0, v_{AL}) }[/math] je: [math]\displaystyle{ v_{AL} = \frac{p_{a}^{(L)}}{E_a^{(L)}} }[/math].

Polárny a azimutálny uhol sekundárnych častíc v sférických súradniciach

Označme štvorhybnosť produkovanej častice [math]\displaystyle{ i }[/math] ako [math]\displaystyle{ P_i =(E_i, \bf{p}_i) }[/math]. Potom polárny uhol (uhol výletu) v jednotlivých sústavách je definovaný nasledovne:

• Laboratorná sústava: [math]\displaystyle{ cos\theta^{(L)} = \frac{\bf{p}_{i}^{(L)}\cdot \bf{p}_{a}^{(L)}}{p_{i}^{(L)}p_{a}^{(L)}} }[/math].
• Ťažisková sústava: [math]\displaystyle{ cos\theta^{*} = \frac{\bf{p}_{i}^{*} \cdot \bf{p}_{a}^{*}}{p_{i}^{*}p_{a}^{*}} }[/math].
• Antilaboratorná sústava: [math]\displaystyle{ cos\theta^{(AL)} = \frac{\bf{p}_{i}^{(AL)} \cdot \bf{p}_{b}^{(AL)}}{p_{i}^{(AL)}p_{b}^{(AL)}} }[/math].
• Azimutálny uhol je invariantný : [math]\displaystyle{ \tan\phi_{i} = \frac{p_{yi}}{p_{xi}} }[/math]

Generovanie prípadov pružného rozptylu dp pri hybnosti deuterona 3,35 GeV/c

Vstupné údaje:[math]\displaystyle{ P_{d}^{(L)}=(E_{d}^{(L)}, 0, 0, p_{d}^{(L)}) }[/math] a [math]\displaystyle{ P_{p}^{(L)}=(m_{p}, 0, 0, 0) }[/math], kde [math]\displaystyle{ E_{d}^{(L)} = \sqrt{m_{d}^{2} + p_{d}^{(L)2}} }[/math] energia deuterona, [math]\displaystyle{ p_{d} = 3.35 GeV/c }[/math] hybnosť deuterona, [math]\displaystyle{ m_{d} }[/math] a [math]\displaystyle{ m_{p} }[/math] hmotnosť deuterona a protona.

Generovanie je rozumné urobiť v tažiskovej sústave, kde sekundárne častice sú rozdelené izotropne a použiť polárne súradnice:

• v počiatočnom stave vypočítať hybnosť v ťažiskovom systeme [math]\displaystyle{ \bf{p}^* = {\bf{p}}_{d}^* = -{\bf{p}}_{p}^* }[/math]
• generovať rovnomerné rozdelenie [math]\displaystyle{ cos\theta^* }[/math],
• generovať rovnomerné rozdelenie azimutálneho uhla [math]\displaystyle{ \phi }[/math],
• vypočítať energiu rozptýleného deuterona [math]\displaystyle{ E_{d}^{\prime*} = \sqrt{m_{d}^2 + p^{*2}} }[/math],
• prejsť do kartézskych súradníc (x,y,z): [math]\displaystyle{ p_{dx}^{\prime*} = p^{*}sin\theta^{*} cos\phi }[/math], [math]\displaystyle{ p_{dy}^{\prime*} = p^{*}sin\theta^{*} sin\phi }[/math], [math]\displaystyle{ p_{dz}^{\prime*} = p^{*}cos\theta^{*} }[/math]
• odrazenému protonu priradiť hybnosť [math]\displaystyle{ \bf{p}_{p}^{\prime *} = -\bf{p}_{d}^{\prime *} }[/math] a energiu [math]\displaystyle{ E_{p}^{\prime*} = \sqrt{m_{p}^2 + p^{*2}} }[/math]
• postupne transformovať štvorhybnosti [math]\displaystyle{ P_{d}^{\prime*} =(E_{d}^{\prime*}, \bf{p}_{d}^{\prime *}) }[/math] a [math]\displaystyle{ P_{p}^{\prime*} =(E_{p}^{\prime*}, \bf{p}_{p}^{\prime *}) }[/math] do laboratornej a antilaboratornej sústavy a zostrojiť pre sekundárne častice rozdelenia modulov hybností, polárneho uhla, polárneho uhla versus modul hybnosti a/alebo cosinus polárneho uhla versus modul hybnosti
• obdržané rozdelenie veličiny kvadrátu prenesenej štvorhybnosti od počiatocného protonu k odrazenému protonu [math]\displaystyle{ t = (P_{p} - P_{p}^{\prime})^2 }[/math] je [math]\displaystyle{ t \lt 0 }[/math] a je rovnomerné. Tento výsledok sa nezhoduje s experimentálnymi výsledkami, kde [math]\displaystyle{ t }[/math] je rozdelené exponenciálne [math]\displaystyle{ \frac{dN}{dt} \sim exp(bt) }[/math] so sklonom [math]\displaystyle{ b \approx 20 (GeV)^{-2} }[/math]. Možné korekcie:
  • každému generovanému prípadu priradíme váhu [math]\displaystyle{ w = b \cdot exp(bt) }[/math]
  • generujeme [math]\displaystyle{ t }[/math] s rozdelením [math]\displaystyle{ exp(bt) }[/math] a polárny uhol v ťažiskovej sústave sa dá potom určiť zo vzťahu [math]\displaystyle{ cos\theta^* = 1 + \frac{t}{2p^{*2}} }[/math]

Realizácia v ROOT-e

Vzhľadom k demonštračnému charekteru úlohy pre generovanie rozdelení môžeme použiť ľubovoľný generátor náhodných čísiel napr. ten nastavený v gRandom. Pre Lorentzovu transformáciu štvorhyností TLorentzVector sú k dipozicii BoostVector() a Boost(). Pozor BoostVector() vráti rýchlosť pôvodného systemu voči čiarkovanému. Pre hmotnosti stačí použiť približné hodnoty: [math]\displaystyle{ m_{d} = 1.88 GeV }[/math], [math]\displaystyle{ m_{p} = 0.938 GeV }[/math], presnejšie hodnoty možno získať z databázy po spustení makra writeujPDG.C, ako je ilustrované v ďalších priložených makrách dp2dpByHand.C a dp3_35Elastic.C. Pred použitím týchto makier je potrebné doinštalovať do ROOT-u balik Monte Carlo nasledujúcim príkazom:

• "sudo dnf install root-montecarlo-eg" .

Ako kontrolu výsledkov generovania možno využiť TGenPhaseSpace (). Príklad jeho použitia https://root.cern/doc/v618/PhaseSpace_8C_source.html alebo priložené macra ROOT-u. Pre obdržanie realistického rozdelenia kvadrátu prenesenej štvorhybnosti od počiatocného protonu k odrazenému protonu generovanému prípadu priradíme váhu [math]\displaystyle{ w = b \cdot exp(bt) }[/math].

Prílohy

Makro na doplnenie databázy častíc

Ak nie je Monte Carlo balik v roote treba dat prikaz (Fedora Silverblue toolbox):

sudo dnf install root-montecarlo-eg

void writeujPDG(const char* filename="ujpdg.txt") {
  const Int_t kion=1000000000;
  const Double_t khSlash = 1.0545726663e-27;
  const Double_t kErg2Gev = 1/1.6021773349e-3;
  const Double_t khShGev = khSlash*kErg2Gev;
  const Double_t kYear2Sec = 3600*24*365.25;
  // IONS
  // charge in units of |e|/3
  TDatabasePDG *pdgDB = new TDatabasePDG();
  Int_t ionCode = kion+10020;
  if(!pdgDB->GetParticle(ionCode)){
      pdgDB->AddParticle("Deuteron","Deuteron", 1.875613, kTRUE,
                         0,3,"Ion",ionCode);
  }
  pdgDB->AddAntiParticle("AntiDeuteron", - ionCode);
  ionCode = kion+10030;
  if(!pdgDB->GetParticle(ionCode)){
    pdgDB->AddParticle("Triton","Triton", 2.80925, kFALSE,
                     khShGev/(12.33*kYear2Sec),3,"Ion",ionCode);
  }
  pdgDB->AddAntiParticle("AntiTriton", - ionCode);
  ionCode = kion+20030;
  if(!pdgDB->GetParticle(ionCode)){
    pdgDB->AddParticle("HE3","HE3", 2.80923,kFALSE,
                     0,6,"Ion",ionCode);
  }
  pdgDB->AddAntiParticle("AntiHE3", - ionCode);
  ionCode = kion+20040;
  if(!pdgDB->GetParticle(ionCode)){
    pdgDB->AddParticle("Alpha","Alpha", 3.727379, kTRUE,
                       khShGev/(12.33*kYear2Sec), 6, "Ion", ionCode);
  }
  pdgDB->AddAntiParticle("AntiAlpha", - ionCode);
  pdgDB->WritePDGTable(filename);
}

Makro na generovanie prípadov dp pružného rozptylu - ručne

void dp2dpByHand(){
// dp elastic @ d lab mom 3,35 GeVC/c   
   const char *tunit="(GeV)^{2}", *punit="GeV/c", *Eunit="GeV";
// Get masses
//Momentum, Energy units are Gev/c, GeV   
   TDatabasePDG *pdg = new TDatabasePDG();
   pdg -> ReadPDGTable("ujpdg.txt");
   TParticlePDG * pp = pdg->GetParticle("Deuteron");
   Double_t mdeuteron = pp->Mass();
   pp = pdg->GetParticle("proton");
   Double_t mproton = pp->Mass();
   Double_t mtarget = mproton;
   Double_t mbeam = mdeuteron ;
   Double_t pbeamMod = 3.35; // GeV/c beam pz
  cout <<  "mdeuteron = " << mdeuteron << "   mproton = " << mproton << endl;
   TVector3 pbeam(0, 0, pbeamMod); // Lab frame
   TLorentzVector target(0.0, 0.0, 0.0, mtarget);
   TLorentzVector beam(0.0, 0.0, pbeamMod, sqrt(pow(pbeamMod, 2) + pow(mbeam, 2)) );
   TLorentzVector W = beam + target;
   Double_t s =  W.Mag2() ;
   cout << "beam   Lab 4-mom = " ;
   beam.Print();
   cout << "target Lab 4-mom = " ;
   target.Print();
       cout << endl ;
   Double_t finalMasses[] = {mdeuteron, mproton};
   Int_t nfinalPart = end(finalMasses) - begin(finalMasses);
   cout << "nfinalPart = " << nfinalPart << endl;
   cout << "finalMasses = " << finalMasses[0] << "  " << finalMasses[1] << endl;
       cout << endl ;
 TH2D *hDLabmomvsPolAngle = new TH2D("hDLabmomvsPolAngle", "Deuteron Lab mom mod vs cos polar angle ; p [GeV/c] ;cos#theta", 20, 0, 4, 40, -1, 1);
   TH2D *hPLabmomvsPolAngle = new TH2D("hPLabmomvsPolAngle", "Proton Lab mom mod vs cos polar angle ; p [GeV/c] ;cos#theta", 20, 0, 4, 40, -1, 1);
   TH2D *hLabPmomvsDmom = new TH2D("hLabPmomvsDmom", "Proton Lab mom mod vs Deuteron Lab mom mod ;  p_{D} [GeV/c]; p_{p} [GeV/c]", 20, 0, 4, 20, 0, 4);
   TH1D *hchecksum = new TH1D("hchecksum", "Check sum ", 100, -1, 1);
   TH1D *hDCMScosth = new TH1D("hDCMScosth", "D CMS cos polar angle; cos#theta; N", 40, -1, 1 );
   hDCMScosth ->GetYaxis()->SetTitle(Form("N/%g", hDCMScosth->GetBinWidth(1))) ;
   TH1D *hPCMScosth = new TH1D("hPCMScosth", "P CMS cos polar angle; cos#theta; N", 40, -1, 1 );
   hPCMScosth ->GetYaxis()->SetTitle(Form("N/%g", hPCMScosth->GetBinWidth(1))) ;
   TH1D *hDALABcosth = new TH1D("hDALABcosth", "D ALAB cos polar angle; cos#theta; N", 40, -1, 1 );
   hDALABcosth ->GetYaxis()->SetTitle(Form("N/%g", hDALABcosth->GetBinWidth(1))) ;
   TH1D *hPALABcosth = new TH1D("hPALABcosth", "P ALAB cos polar angle; cos#theta; N", 40, -1, 1 );
   hPALABcosth ->GetYaxis()->SetTitle(Form("N/%g", hPALABcosth->GetBinWidth(1))) ;
   TH1D *htp_p = new TH1D("htp_p", "4-mom squared p-p; t [(GeV)^{2}]", 5000, -5, 0 );
   htp_p ->GetYaxis()->SetTitle(Form("N/(%g %s)", htp_p->GetBinWidth(1), tunit)) ;
   Double_t deltat = htp_p -> GetBinWidth(1);
   Double_t weight = 1 ;
//     Lorentz boost Lab --> CMS
       TVector3 v = (beam + target).BoostVector();   // Lab frame velocity w.r.t. CMS
       TLorentzVector pbeam4CMS = beam;
       pbeam4CMS.Boost(-v);  // Lab -->  CMS
       cout << "beam   CMS 4-mom = " ;
       pbeam4CMS.Print();
       TLorentzVector ptarget4CMS = target;
       ptarget4CMS.Boost(-v);
       cout << "target CMS 4-mom = " ;
       ptarget4CMS.Print();
       cout << endl ;
//    Lorentz boost Lab --> ALAB
       TVector3 vALAB2Lab = (beam).BoostVector();   // ALAB frame velocity w.r.t. Lab
       TLorentzVector pbeam4ALAB = beam;
       pbeam4ALAB.Boost(-vALAB2Lab);  // Lab -->  ALAB
       cout << "beam   ALAB 4-mom = " ;
       pbeam4ALAB.Print();
       TLorentzVector ptarget4ALAB = target;
       ptarget4ALAB.Boost(-vALAB2Lab);
       cout << "target ALAB 4-mom = " ;
       ptarget4ALAB.Print();
      Int_t option(0); // 1 w=1,  2 t generated exp(20t), 3 event reweigthed w = 20 exp(20t)   ;  t < 0 !!!!
      option = 3 ;
   Int_t ngen = 1000000 ;
   for (int n = 0 ; n < ngen ; n++){
//    Secondaries in CMS
       Double_t pCMS =  pbeam4CMS.P();
       Double_t costhCMS ;
       if (option == 1 || option == 3) {costhCMS = gRandom->Uniform(-1, 1) ;}
      Double_t t(0);
      if (option == 2 ){
// generate f(t) = 20 * exp(20*t) , t <0
       costhCMS = 0.5 * gRandom->Exp(-0.05)/pow(pCMS, 2) + 1 ; 
      }
       Double_t phi = gRandom->Uniform(0, 2*TMath::Pi()) ;
       Double_t sinthCMS = sqrt(1-costhCMS * costhCMS);
      TVector3 pD3CMS(pCMS*sinthCMS*cos(phi), pCMS*sinthCMS*sin(phi), pCMS*costhCMS);
      TLorentzVector pD4CMS(pD3CMS, sqrt(pCMS*pCMS + mdeuteron*mdeuteron) );
      TVector3 pP3CMS = -1*pD3CMS ;
      TLorentzVector pP4CMS(pP3CMS, sqrt(pCMS*pCMS + mproton*mproton) );
      TVector3 pP3CMS = -1*pD3CMS ;
      TLorentzVector pP4CMS(pP3CMS, sqrt(pCMS*pCMS + mproton*mproton) );
       TLorentzVector p4CMSSum = pD4CMS + pP4CMS ;
       Double_t sum = p4CMSSum.P();
       //cout << sum << endl;
       hchecksum->Fill(sum, weight);
// CMS --> Lab
      TLorentzVector pD4Lab = pD4CMS ;
      pD4Lab.Boost(v);
      TLorentzVector pP4Lab = pP4CMS ;
      pP4Lab.Boost(v);
//  Lab --> ALAB
      TLorentzVector pD4ALAB = pD4Lab;
      pD4ALAB.Boost(-vALAB2Lab);  // Lab -->  ALAB
      TLorentzVector pP4ALAB = pP4Lab;
      pP4ALAB.Boost(-vALAB2Lab);  // Lab -->  ALAB
       t = (pP4Lab - target).Mag2();
      if (option == 3 ) { 
// weighting  by f(t) = 20 * exp(20*t) , t <0
       weight = 20 * exp(20 * t) ;    
      }
        hDLabmomvsPolAngle->Fill( pD4Lab.P(), pD4Lab.CosTheta(), weight);
        hPLabmomvsPolAngle->Fill( pP4Lab.P(), pP4Lab.CosTheta(), weight);
        hLabPmomvsDmom->Fill( pD4Lab.P(), pP4Lab.P(), weight);
       hDCMScosth->Fill(pD4CMS.CosTheta(), weight);
       hPCMScosth->Fill(pP4CMS.CosTheta(), weight);
       hDALABcosth->Fill( ((pD4ALAB.Vect()).Unit()).Dot((ptarget4ALAB.Vect()).Unit()), weight);
       hPALABcosth->Fill( ((pP4ALAB.Vect().Unit())).Dot((ptarget4ALAB.Vect()).Unit()), weight);
       htp_p->Fill(t, weight);
   }
   gStyle->SetOptStat(111111111);
 TCanvas *c1 = new TCanvas("c1", "Deuteron - c1");
 hDLabmomvsPolAngle->Draw() ;
 TCanvas *c2 = new TCanvas("c2", "Proton - c3");
 hPLabmomvsPolAngle->Draw() ;
 TCanvas *c3 = new TCanvas("c3", "P vs D - c3");
 hLabPmomvsDmom->Draw() ;
 TCanvas *c4 = new TCanvas("c4", "PCMS+ DCMS mom modul - c4");
 hchecksum->Draw();
 TCanvas *c5 = new TCanvas("c5", "DCMS costh - c5");
 hDCMScosth->Draw();
 TCanvas *c6 = new TCanvas("c6", "PCMS costh - c6");
 hPCMScosth->Draw();
 TCanvas *c8 = new TCanvas("c8", "PCMS costh - c8");
 hDALABcosth->Draw();
 TCanvas *c9 = new TCanvas("c9", "PCMS costh - c9");
 hPALABcosth->Draw();
 TCanvas *c7 = new TCanvas("c7", "t - c7");
 htp_p->Draw();
 if(option == 1 ) {return ;}                            
 // htp_p normalization   20*exp(20*t)
 TH1D * htp_pnorm = (TH1D *) htp_p->Clone("htp_pnorm");
 htp_pnorm -> Scale(1/htp_p->GetSumOfWeights()/deltat);
 TCanvas *c71 = new TCanvas("c71", "t normalized - c71");
 htp_pnorm ->Draw();
 htp_pnorm ->Fit("expo", "", "", -5, 0);
 TF1 *f =  htp_pnorm->GetFunction("expo");
 cout << "p0 = " << f->GetParameter(0) << " +- " << f->GetParError(0) <<"   exp(p0) = " << exp(f->GetParameter(0)) <<  " +- " << exp(f->GetParameter(0))*f->GetParError(0) << endl;
 cout << "slope p1 = " << f->GetParameter(1) << " +- " << f->GetParError(1)<< endl;
}

Výsledky generovania prípadov dp pružného rozptylu

beam mass = 1.87561   target mass = 0.938272
beam   Lab 4-mom :
(x,y,z,t)=(0.000000,0.000000,3.350000,3.839325) 
(P,eta,phi,E)=(3.350000,100000000000.000000,0.000000,3.839325)

target Lab 4-mom : 
(x,y,z,t)=(0.000000,0.000000,0.000000,0.938272) 
(P,eta,phi,E)=(0.000000,0.000000,0.000000,0.938272)

s = 11.6029

nfinalPart = 2
finalMasses = 1.87561,  0.938272

Lab frame velocity w.r.t. CMS v = TVector3 A 3D physics vector:
(x,y,z)=(0.000000,0.000000,0.701189) 
(rho,theta,phi)=(0.701189,0.000000,0.000000)

beam   CMS 4-mom :
(x,y,z,t)=(0.000000,0.000000,0.922762,2.090312) 
(P,eta,phi,E)=(0.922762,100000000000.000000,0.000000,2.090312)

target CMS 4-mom : 
(x,y,z,t)=(0.000000,0.000000,-0.922762,1.315996) 
(P,eta,phi,E)=(0.922762,-100000000000.000000,0.000000,1.315996)

Lab frame velocity w.r.t. ALab vALAB2Lab TVector3 A 3D physics vector :
(x,y,z)=(0.000000,0.000000,0.872549) 
(rho,theta,phi)=(0.872549,0.000000,0.000000)

beam   ALAB 4-mom : 
(x,y,z,t)=(0.000000,0.000000,0.000000,1.875610) 
(P,eta,phi,E)=(0.000000,0.000000,0.000000,1.875610)

target ALAB 4-mom : 
(x,y,z,t)=(0.000000,0.000000,-1.675834,1.920618) 
(P,eta,phi,E)=(1.675834,-100000000000.000000,0.000000,1.920618)

Scattered deuteron CMS 4-mom : 
(x,y,z,t)=(0.015429,0.025324,0.922285,2.090312) 
(P,eta,phi,E)=(0.922762,4.130660,1.023593,2.090312)

Recoiled  proton   CMS 4-mom : 
(x,y,z,t)=(-0.015429,-0.025324,-0.922285,1.315996) (P,eta,phi,E)=(0.922762,-4.130660,-2.118000,1.315996)

Scattered deuteron LAB 4-mom : 
(x,y,z,t)=(0.015429,0.025324,3.349332,3.838856) 
(P,eta,phi,E)=(3.349463,5.420082,1.023593,3.838856)

Recoiled  proton   LAB 4-mom : 
(x,y,z,t)=(-0.015429,-0.025324,0.000668,0.938741) 
(P,eta,phi,E)=(0.029662,0.022541,-2.118000,0.938741)

Scattered deuteron ALAB 4-mom : 
(x,y,z,t)=(0.015429,0.025324,-0.000531,1.875844) 
(P,eta,phi,E)=(0.029659,-0.017911,1.023593,1.875844)

Recoiled  proton   ALAB 4-mom : 
(x,y,z,t)=(-0.015429,-0.025324,-1.675303,1.920384) (P,eta,phi,E)=(1.675565,-4.727374,-2.118000,1.920384)

4-mom squaqred form target proton to recoiled proton (Lab fr evaluation)  t = -0.000879593

Event weight = 1

Makro na generovanie prípadov dp pružného rozptylu - TGenPhaseSpace

void dp3_35Elastic(){
// dp elastic @ d lab mom 3,35 GeVC/c   
   const char *tunit="(GeV)^{2}", *punit="GeV/c", *Eunit="GeV";
// Get masses
//Momentum, Energy units are Gev/c, GeV   
   TDatabasePDG *pdg = new TDatabasePDG();
   pdg -> ReadPDGTable("ujpdg.txt");
   TParticlePDG * pp = pdg->GetParticle("Deuteron");
   Double_t mdeuteron = pp->Mass();
   pp = pdg->GetParticle("proton");
   Double_t mproton = pp->Mass();
   Double_t mtarget = mproton;
   Double_t mbeam = mdeuteron ;
   Double_t pbeamMod = 3.35; // GeV/c beam pz
  cout <<  "mbeam " << mbeam << "   mtarget = " << mtarget << endl;
   TVector3 pbeam(0, 0, pbeamMod); // Lab frame
   TLorentzVector target(0.0, 0.0, 0.0, mtarget);
   TLorentzVector beam(0.0, 0.0, pbeamMod, sqrt(pow(pbeamMod, 2) + pow(mbeam, 2)) );
   TLorentzVector W = beam + target;
   Double_t s =  W.Mag2() ;
   cout << "beam   Lab 4-mom = " ;
   beam.Print();
   cout << "target Lab 4-mom = " ;
   target.Print();
       cout << endl ;
//     Lorentz boost Lab --> CMS
       TVector3 vinit = (beam + target).BoostVector();   // Lab frame velocity w.r.t. CMS
       TLorentzVector beamCMS = beam;
       beamCMS.Boost(-vinit);  // Lab -->  CMS
       cout << "beam   CMS 4-mom = " ;
       beamCMS.Print();
       TLorentzVector targetCMS = target;
       targetCMS.Boost(-vinit);
       cout << "target CMS 4-mom = " ;
       targetCMS.Print();
       cout << endl ;
   Double_t finalMasses[] = {mdeuteron, mproton};
   Int_t nfinalPart = end(finalMasses) - begin(finalMasses);
   cout << "nfinalPart = " << nfinalPart << endl;
   cout << "finalMasses = " << finalMasses[0] << "  " << finalMasses[1] << endl;
       cout << endl ;
   TGenPhaseSpace dpelastic ;
//   dpelastic.SetDecay(W, nfinalPart, finalMasses, "Fermi");
   dpelastic.SetDecay(W, nfinalPart, finalMasses);
   TH2D *hDLabmomvsPolAngle = new TH2D("hDLabmomvsPolAngle", "Deuteron Lab mom mod vs cos polar angle ; p [GeV/c] ;cos#theta ", 20, 0, 4, 40, -1, 1);
   TH2D *hPLabmomvsPolAngle = new TH2D("hPLabmomvsPolAngle", "Proton Lab mom mod vs cos polar angle ; p [GeV/c] ;cos#theta ", 20, 0, 4, 40, -1, 1);
   TH2D *hLabPmomvsDmom = new TH2D("hLabPmomvsDmom", "Proton Lab mom mod vs Deuteron Lab mom mod ; p_{D} [GeV/c] ; p_{P} [GeV/c]", 20, 0, 4, 20, 0, 4);
   TH1D *hchecksum = new TH1D("hchecksum", "Check sum ", 100, -1, 1);
   TH1D *hDCMScosth = new TH1D("hDCMScosth", "D CMS cos polar angle; cos#theta^{*}; N", 40, -1, 1 );
   hDCMScosth ->GetYaxis()->SetTitle(Form("N/%g", hDCMScosth->GetBinWidth(1))) ;
   TH1D *hPCMScosth = new TH1D("hPCMScosth", "P CMS cos polar angle;  cos#theta^{*}; N", 40, -1, 1 );
   hPCMScosth ->GetYaxis()->SetTitle(Form("N/%g", hPCMScosth->GetBinWidth(1))) ;
   TH1D *hDphi = new TH1D("hDphi", "D CMS azimuthal angle; #phi [rad]; N", 40, -3.2, 3.2 );
   hDphi ->GetYaxis()->SetTitle(Form("N/%g", hDphi->GetBinWidth(1))) ;
   TH1D *hPphi = new TH1D("hPphi", "P CMS azimuthal angle; #phi [rad]; N", 40, -3.2, 3.2 );
   hPphi ->GetYaxis()->SetTitle(Form("N/%g", hPphi->GetBinWidth(1))) ;
   TH1D *htp_p = new TH1D("htp_p", "4-mom squared p-p; t [(GeV)^{2}]; N", 500, -5, 0 );
   htp_p ->GetYaxis()->SetTitle(Form("N/(%g %s)", htp_p->GetBinWidth(1), tunit)) ;
   TH1D *htd_d = new TH1D("htd_d", "4-mom squared d_d; t [(GeV)^{2}]; N", 500, -5, 0 );
   htd_d ->GetYaxis()->SetTitle(Form("N/(%g %s)", htd_d->GetBinWidth(1), tunit)) ;
   TH1D *htp_pCMS = new TH1D("htp_pCMS", "4-mom squared p-p CMS; t [(GeV)^{2}]; N", 500, -5, 0 );
   htp_pCMS ->GetYaxis()->SetTitle(Form("N/(%g %s)", htp_pCMS->GetBinWidth(1), tunit)) ;
   TH1D *htd_dCMS = new TH1D("htd_dCMS", "4-mom squared d_d CMS; t [(GeV)^{2}]; N", 500, -5, 0 );
   htd_dCMS ->GetYaxis()->SetTitle(Form("N/(%g %s)", htd_dCMS->GetBinWidth(1), tunit)) ;
   TH1D *htp_p_reweighted = new TH1D("htp_p_reweighted", "4-mom squared p-p reweighted; t [(GeV)^{2}]; N", 500, -5, 0 );
   Int_t ngen = 100000 ;
   for (int n = 0 ; n < ngen ; n++){
        Double_t weight = dpelastic.Generate();
        TLorentzVector pD4LAB(*dpelastic.GetDecay(0));
        TVector3 pD3Lab = pD4LAB.Vect();
        TLorentzVector pP4LAB(*dpelastic.GetDecay(1));
        TVector3 pP3Lab = pP4LAB.Vect();
        hDLabmomvsPolAngle->Fill( pD3Lab.Mag(), (pD3Lab.Unit()).Dot(pbeam.Unit()), weight);
        hPLabmomvsPolAngle->Fill( pP3Lab.Mag(), (pP3Lab.Unit()).Dot(pbeam.Unit()), weight);
        hLabPmomvsDmom->Fill( pD3Lab.Mag(), pP3Lab.Mag(), weight);
//     Lorentz boost to CMS, it is superflouos ... ony to demostrate how it works
       TVector3 v ;
       v = (pD4LAB + pP4LAB).BoostVector();
       TLorentzVector pD4CMS = pD4LAB;
       pD4CMS.Boost(-v);
       TLorentzVector pP4CMS = pP4LAB;
       pP4CMS.Boost(-v);
       Double_t tCMS = (pP4CMS - targetCMS).Mag2();
       Double_t td_dCMS = (pD4CMS - beamCMS).Mag2();
       TLorentzVector p4CMSSum = pD4CMS + pP4CMS ;
       Double_t sum = p4CMSSum.P();
       hchecksum->Fill(sum, weight);
       hDCMScosth->Fill(pD4CMS.CosTheta(), weight);
       hPCMScosth->Fill(pP4CMS.CosTheta(), weight);
       hDphi->Fill(pD4CMS.Phi(), weight);
       hPphi->Fill(pP4CMS.Phi(), weight);
       htp_p->Fill(t, weight);
       htd_d->Fill(td_d, weight);
       htp_pCMS->Fill(tCMS, weight);
       htd_dCMS->Fill(td_dCMS, weight);
// reweighting  by f(t) = 20 * exp(20*t) , t <0, only t reweigthed to show it works
       weight = 20 * exp(20 * t) ;
       htp_p_reweighted -> Fill(t, weight);
   }
   gStyle->SetOptStat(111111111);
 TCanvas *c1 = new TCanvas("c1", "Deuteron");
 hDLabmomvsPolAngle->Draw() ;
 TCanvas *c2 = new TCanvas("c2", "Proton");
 hPLabmomvsPolAngle->Draw() ;
 TCanvas *c3 = new TCanvas("c3", "P vs D");
 hLabPmomvsDmom->Draw() ;
 TCanvas *c4 = new TCanvas("c4", "PCMS+ DCMS mom modul");
 hchecksum->Draw();
 TCanvas *c5 = new TCanvas("c5", "DCMS costheta");
 hDCMScosth->Draw();
 TCanvas *c6 = new TCanvas("c6", "PCMS costheta");
 hPCMScosth->Draw();
 TCanvas *c7 = new TCanvas("c7", "P-P Lab t");
 htp_p->Draw();
 TCanvas *c71 = new TCanvas("c71", "D-D Lab t");
 htd_d->Draw();
 TCanvas *c72 = new TCanvas("c72", "P-P CMS t");
 htp_pCMS->Draw();
 TCanvas *c73 = new TCanvas("c73", "D-D CMS t");
 htd_dCMS->Draw();
 TCanvas *c8 = new TCanvas("c8", "D CMS phi");
 hDphi-> Draw();
 TCanvas *c9 = new TCanvas("c9", "P CMS phi");
 hPphi-> Draw();
 TCanvas *c7rew = new TCanvas("c7rew", "P-P Lab t");
 htp_p_reweighted -> Scale(1/htp_p_reweighted->GetSumOfWeights()/htp_p_reweighted->GetBinWidth(1));
 htp_p_reweighted ->Draw();
 htp_p_reweighted ->Fit("expo", "", "", -5, 0);
 TF1 *f =  htp_p_reweighted->GetFunction("expo");
 cout << "p0 = " << f->GetParameter(0) << " +- " << f->GetParError(0) <<"   exp(p0) = " << exp(f->GetParameter(0)) <<  " +- " << exp(f->GetParameter(0))*f->GetParError(0) << endl;
 cout << "slope p1 = " << f->GetParameter(1) << " +- " << f->GetParError(1)<< endl;
}