GSI Forum
GSI Helmholtzzentrum für Schwerionenforschung

Home » Hades » Pluto » [SOLVED] Beam Smear and Cross Sections
Re: Bear Smear and Cross Sections [message #13877 is a reply to message #13873] Thu, 23 August 2012 05:05 Go to previous messageGo to previous message
Michael Kunkel is currently offline  Michael Kunkel
Messages: 53
Registered: June 2011
continuous participant
From: *hr.hr.cox.net
Greetings,

I am not sure if I am doing something wrong, or the model technique is flawed, or both.

I have a distribution cos(theta) vs. differential cross section, see plot below
index.php?t=getfile&id=7092&private=0

I then take this TGraph, and by using TSpline5, draw a histogram for PScatterCrossSection, see figure below.

index.php?t=getfile&id=7094&private=0

I then use this histogram as seen in the sample macro
Input: _x is c.m. cos(theta), _y is the dsig/dcos(theta)
Output: _f: cross section
    model->AddHistogram(example,"value = Eval(_x); _f = _y * value");


The output of PLUTO generated looks nothing like the inputted, as it can be seen from the first 2 figure to that of the figure below

index.php?t=getfile&id=7095&private=0


Moreover, the macro I use stipulates the c.m energy range as

model->SetRange(2.2,2.22); //in GeV


which corresponds to 2.11 to 2.16 GeV beam energy in lab frame for the reaction g p -> p eta.
As can be seen from the figures below, from using the functions

my_reaction.Do(histo1,"_x = [g+p]->M()");

index.php?t=getfile&id=7098&private=0
for c.m. energy, and

my_reaction.Do(beamProfile,"_x = ([g+p]->GetBeam())->E();");

index.php?t=getfile&id=7097&private=0
for lab energy, neither of these plots correspond to the input given in the macro.

I can no longer attach any more files, so I am adding my macro in the text body.
Another note, when using ROOTs ALCiC to compile the macro, I had to add these lines into the PScatterCrossSection.h

#include "/Users/Mike/Pluto/pluto_v5.40.5/src/PAngularDistribution.h"
#include "/Users/Mike/Pluto/pluto_v5.40.5/src/PProjector.h"
#include "/Users/Mike/Pluto/pluto_v5.40.5/src/PF2EvalBatch.h"



Here is my macro



#include "TGraphAsymmErrors.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TChain.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TGraphAsymmErrors.h"
#include "/Users/Mike/Pluto/pluto_v5.40.5/src/PParticle.h"
#include "/Users/Mike/Pluto/pluto_v5.40.5/src/PReaction.h"
#include "/Users/Mike/Pluto/pluto_v5.40.5/src/PBeamSmearing.h"
#include "/Users/Mike/Pluto/pluto_v5.40.5/src/PAnyDistribution.h"
#include "/Users/Mike/Pluto/pluto_v5.40.5/plugins/scatter_mod/PScatterCrossSection.h"


// 2.2 - 2.22 

void eta_XSection()
{
	
  
  //From data on eta differential cross section as function of c.m. cos(theta) for c.m. energy 2.2 -> 2.22 GeV
  //Making a TGraph then fit the TGraph to a spline
  const int p7694_d47x1y1_numpoints = 17;

  double p7694_d47x1y1_xval[p7694_d47x1y1_numpoints] = { -0.84, -0.75, -0.65, -0.55, -0.45, -0.35, -0.25, -0.15, -0.05, 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75 };
	double p7694_d47x1y1_xerrminus[p7694_d47x1y1_numpoints] = { 0.040000000000000036, 0.050000000000000044, 0.04999999999999993, 0.04999999999999993, 0.04999999999999999, 0.050000000000000044, 0.04999999999999999, 0.05000000000000002, 0.05, 0.05, 0.04999999999999999, 0.04999999999999999, 0.04999999999999999, 0.04999999999999999, 0.050000000000000044, 0.050000000000000044, 0.050000000000000044 };
	double p7694_d47x1y1_xerrplus[p7694_d47x1y1_numpoints] = { 0.039999999999999925, 0.050000000000000044, 0.050000000000000044, 0.050000000000000044, 0.04999999999999999, 0.04999999999999999, 0.04999999999999999, 0.04999999999999999, 0.05, 0.05, 0.05000000000000002, 0.04999999999999999, 0.050000000000000044, 0.04999999999999999, 0.04999999999999993, 0.04999999999999993, 0.050000000000000044 };
	double p7694_d47x1y1_yval[p7694_d47x1y1_numpoints] = { 0.0324, 0.0374, 0.0276, 0.0233, 0.0173, 0.0122, 0.0122, 0.0169, 0.0277, 0.0333, 0.0359, 0.0375, 0.0439, 0.0625, 0.093, 0.132, 0.1634 };
	double p7694_d47x1y1_yerrminus[p7694_d47x1y1_numpoints] = { 0.005, 0.00291, 0.00185, 0.00166, 0.00139, 0.00102, 9.7E-4, 0.00117, 0.00178, 0.00205, 0.00215, 0.00217, 0.00258, 0.00386, 0.005, 0.0069, 0.0093 };
	double p7694_d47x1y1_yerrplus[p7694_d47x1y1_numpoints] = { 0.005, 0.00291, 0.00185, 0.00166, 0.00139, 0.00102, 9.7E-4, 0.00117, 0.00178, 0.00205, 0.00215, 0.00217, 0.00258, 0.00386, 0.005, 0.0069, 0.0093 };
	double p7694_d47x1y1_ystatminus[p7694_d47x1y1_numpoints] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
	double p7694_d47x1y1_ystatplus[p7694_d47x1y1_numpoints] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
	TGraphAsymmErrors *p7694_d47x1y1 = new TGraphAsymmErrors(p7694_d47x1y1_numpoints, p7694_d47x1y1_xval, p7694_d47x1y1_yval, p7694_d47x1y1_xerrminus, p7694_d47x1y1_xerrplus, p7694_d47x1y1_yerrminus, p7694_d47x1y1_yerrplus);
	p7694_d47x1y1->SetTitle("Differential Cross Section for #eta c.m. energy 2.2 GeV#rightarrow 2.22 GeV"); ///HepData/7694/d47x1y1
  p7694_d47x1y1->GetXaxis()->SetTitle("Cos(#theta)"); 
	p7694_d47x1y1->GetYaxis()->SetTitle("#frac{d#sigma}{d#Omega}"); 
  TCanvas *TgCan = new TCanvas("TgCan","TgCan");
  TgCan->cd();
	p7694_d47x1y1->Draw("AP");
  
  TSpline5 *s3 = new TSpline5("s3",p7694_d47x1y1);
  s3->Draw("same");

  
  //Now from the spline, create a TH1D histogram for PLUTO to model the Xsection with
  TH1D *example = new TH1D("example","try this",2001,-1,1);
  example->SetTitle("Extrapolated Differential Cross Section for #eta c.m. energy 2.2 GeV#rightarrow 2.22 GeV"); ///HepData/7694/d47x1y1
  example->GetXaxis()->SetTitle("Cos(#theta)"); 
	example->GetYaxis()->SetTitle("#frac{d#sigma}{d#Omega}"); 
  example->Sumw2();
  int hold = 1; //place value for histogram bin filling
//############Filling Histogram#########################  
  for (double j = -1.; j<=1.; j = j +0.001) {
    double bin_content = s3->Eval(double(j));
    if(bin_content < 0.0){bin_content = 0.0;}
    //cout<<bin_content<<endl;
    example->SetBinContent(hold,bin_content);
    hold++;
  }
//################################################  
  TCanvas *canf = new TCanvas("canf","canf");
  canf->cd();
  example->Draw("ep");
  
  
//PLUTO TIME
  
  PScatterCrossSection * model = new PScatterCrossSection("mymodel","My cross section");
  model->Add("g,grandparent,beam");
  model->Add("p,grandparent,target");
  model->Add("q,parent");
  model->Add("p,daughter");
  model->Add("eta,daughter,primary");
  
  //Define the range of the c.m. sampling:
  model->SetRange(2.2,2.22);  //CLAS G12 c.m. model->SetRange(1.739,3.3268);
  

  //model->AddEquation("my_x = s3->Eval(_x);;_f = Eval(my_x)");
  model->AddHistogram(example,"value = Eval(_x); _f =_y * value");
  makeDistributionManager()->Add(model);
  
  
  
  
  PReaction my_reaction("_P1 = 2.113","g","p","p eta [dilepton [e+ e-] g]");
  
  TH1F * histo1 = new TH1F ("histo1","c.m.",100,2.1,2.5);
  histo1->Sumw2();
  histo1->GetXaxis()->SetTitle("c.m. energy [GeV]"); 

  my_reaction.Do(histo1,"_x = [g+p]->M()");
  
  TH2F * histo2 = new TH2F ("histo2","c.m. vs cos_theta",100,1.0,3.0,20,-1,1);
  my_reaction.Do(histo2,"_x = [g+p]->M(); myeta = [eta]; myeta->Boost([g+p]); _y = myeta->CosTheta()");
  
  //This histogram shows the beam profile:
  TH1F * beamProfile = new TH1F ("beamProfile","Beam Spectrum",100,2.0,2.7);
  beamProfile->Sumw2();
  beamProfile->GetXaxis()->SetTitle("E_{#gamma} GeV"); 

  my_reaction.Do(beamProfile,"_x = ([g+p]->GetBeam())->E();");
  
  
  TH1F * histo4 = new TH1F ("histo4","PLUTO generated c.m. Cos#theta",100,-1,1);
  my_reaction.Do(histo4,"myeta->Boost([g+p]); _x = myeta->CosTheta()");
  histo4->GetXaxis()->SetTitle("Cos(#theta)"); 
  histo4->Sumw2();
  
  
  my_reaction.Print();
  
  
  my_reaction.Loop(100000);
  
  
  TCanvas *can1 = new TCanvas("can1","can1");
  can1->cd();
  histo1->Draw();
  
  TCanvas *can2 = new TCanvas("can2","can2");
  can2->cd();
  beamProfile->Draw("ep");
  TCanvas *can3 = new TCanvas("can3","can3");

  can3->cd();
  histo4->Draw("ep");
  
  //TgCan->Print("XSection_Data.jpeg");
  //canf->Print("XSection_Interpolated.jpeg");
  //can1->Print("PLUTO_generated_cm_energy.jpeg");
  //can2->Print("Beam_Profile.jpeg");
  //can3->Print("PLUTO_generated_cos_theta.jpeg");
}






Thanks for all the help
Michael

[Updated on: Thu, 23 August 2012 05:09]

Report message to a moderator

 
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Read Message
Previous Topic: New Pluto web page
Next Topic: [SOLVED] PAnyDistribution broke in v5.40.5
Goto Forum:
  


Current Time: Tue Oct 26 14:34:04 CEST 2021

Total time taken to generate the page: 0.02834 seconds