[Frost] FROST meeting follow up

[Michael Dugger]

Hi,

Durring the meeting I was asked about the use of GENBOD within ROOT.

A simple implementation of GENBOD within ROOT is shown on slides 26-30 of 
the pdf at
http://hadron.physics.fsu.edu/~skpark/document/ROOT/RootLecture/RootLecture120405.pdf

The code of interest is on slide 29. It does a phase space generation for 
the reaction "gamma p -> p pi0" with a photon energy of 400 MeV, and plots 
the pi0 lab frame polar angle.

The code is short (17 lines!) and is very simple to use :)

The code I have been playing with can be found at the bottom of this 
email. It assumes the initial reaction: "gamma C12 -> B11 + R", where R 
has a mass of 1535 MeV and the incident photon is set to 850 MeV. Then he 
reaction: "R -> p eta" takes place after the initial reaction. Once the 
reactions have been simulated, the overall reaction is analyzed as though 
it were "gamma p -> p X". The mass of X is then plotted versus cos(theta_cm).
Also, there is no fermi motion correction and all reactions are phase 
space generated.

Disclaimer:
The code I have is just a trivial hack of the code on slide 29, but it 
might be useful for others. If there are errors, please let me know.

It is easy to change the code for other reactions. For example:
To change from C12 to He4, and from B11 to tritium you just change 
the assignment of mTarg and mRecoil:

mTarg = mB11;
mRecoil = mT;

To change the resonance to the Delta1232 and the meson to pi0, you change 
the assignment of mR and mMeson:

mR = 1232;
mMeson = mPi0;

To change the photon energy from 0.85 to 1.0, you set
eGamma = 1.0;

The variable nToGen is the number of iterations and is currently set to 1 
million.

Take care,
Michael

CODE:
--------------------------------
{
   gROOT->Reset();
   gROOT->SetStyle("Plain");
   gStyle->SetOptStat(0);
   gStyle->SetPalette(1);

   gSystem.Load("libPhysics");

   int nToGen = 1000000;
   double mC12,mB11,mR,uToGeV,eGamma,mEta,mPro,mX,mXsq;
   double mTarg,mRecoil,mHe3,mHe4,mD,mT,mAl27,mMg26;
   double mO16,mN15,mMeson,mPi0;
   uToGeV = 0.931502;
   mC12 = 12.0*uToGeV;
   mB11 = 11.009305*uToGeV;
   mHe3 = 3.016029*uToGeV;
   mHe4 = 4.002603*uToGeV;
   mD = 2.014102*uToGeV;
   mT = 3.016049*uToGeV;
   mAl27 = 26.981539*uToGeV;
   mMg26 = 25.982594*uToGeV;
   mO16 = 15.994915*uToGeV;
   mN15 = 15.000109*uToGeV;

   mR = 1.535;
   mPro = 0.938;
   mEta = 0.547;
   mPi0 = 0.135;

   eGamma = 0.85;
   mMeson = mEta;

   mTarg = mC12;
   mRecoil = mB11;

   //mTarg = mHe4;
   //mRecoil = mT;

   //mTarg = mHe3;
   //mRecoil = mD;

   //mTarg = mAl27;
   //mRecoil = mMg26;

   //mTarg = mO16;
   //mRecoil = mN15;

   TLorentzVector target(0.0, 0.0, 0.0, mTarg);
   TLorentzVector proI(0.0, 0.0, 0.0, mPro);
   TLorentzVector beam(0.0, 0.0, eGamma, eGamma);
   TLorentzVector W = beam + target;
   Double_t masses1[2] = {mRecoil, mR} ;
   Double_t masses2[2] = {mPro,mMeson} ;
   TGenPhaseSpace event1;
   TGenPhaseSpace event2;
   event1.SetDecay(W, 2, masses1);
   TH1D *h100 = new TH1D("h100","",100,0.0,180.0);
   TH1D *h200 = new TH1D("h200","",100,-1.0,1.0);
   TH1D *h300 = new TH1D("h300","",100,-1.0,1.0);

   TH1D *h101 = new TH1D("h101","",100,0.0,180.0);
   TH1D *h201 = new TH1D("h201","",100,-1.0,1.0);
   TH1D *h301 = new TH1D("h301","",100,-1.0,1.0);

   TH1D *h400 = new TH1D("h400","",200,-1.0,1.0);
   TH1D *h500 = new TH1D("h500","",140,0.0,1.4);
   TH2D *h501 = new TH2D("h501","",140,0.0,1.4,10,-1.0,1.0);

   for (Int_t n=0;n<nToGen;n++) {

     //INITIAL REACTION
     event1.Generate();
     TLorentzVector *pB11 = event1.GetDecay(0);
     TLorentzVector pB11cm = *pB11;
     TLorentzVector *pR = event1.GetDecay(1);
     TLorentzVector pRcm = *pR;
     pRcm.Boost(-W.BoostVector());
     pB11cm.Boost(-W.BoostVector());

     h100->Fill(pR->Theta()*57.3);
     h200->Fill(pR->CosTheta());
     h300->Fill(pRcm.CosTheta());

     //RESONANCE DECAY
     TLorentzVector pRtoGo = *pR;
     event2.SetDecay(pRtoGo, 2, masses2);
     event2.Generate();
     TLorentzVector *pPro = event2.GetDecay(0);
     TLorentzVector *pMeson = event2.GetDecay(1);
     TLorentzVector pMesoncm = *pMeson;
     pMesoncm->Boost(-pR->BoostVector());

     h101->Fill(pMeson->Theta()*57.3);
     h201->Fill(pMeson->CosTheta());
     h301->Fill(pMesoncm.CosTheta());

     //ETA AS IF FROM GAMMA P REACTION
     TLorentzVector proF = *pPro;
     TLorentzVector pX = proI + beam - proF;
     TLorentzVector pXcm = pX;
     TLorentzVector cms = proI + beam;
     pXcm->Boost(-cms->BoostVector());

     mXsq = pX->Mag2();
     mX = pX->Mag();
     h400->Fill(mXsq);
     h500->Fill(mX);
     h501->Fill(mX,pXcm.CosTheta());
   }
   h501->Draw("colz");

}