Home » Hades » Pluto » Sampling dilepton mass
Sampling dilepton mass [message #15055] |
Tue, 30 July 2013 20:58 |
Michael Kunkel
Messages: 53 Registered: June 2011
|
continuous participant |
From: *odu.edu
|
|
Greetings,
I am trying to sample the dilepton mass of the pi0 dalitz decay.
I am able to use the line
makeStaticData()->SetParticleMass("dilepton",0.06);
and this samples the dilepton mas from 0.6 -> Mass(pi0) using the code below.
{
TH1F * histo2 = new TH1F ("histo2","IM e+ e-",80,0.58,0.146);
makeStaticData()->SetParticleMass("dilepton",0.06);
double ebeam_min = 1.1725;
double ebeam_max = 2.9;
PBeamSmearing *beam_smear = new PBeamSmearing("beam_smear", "Beam smearing");
TF1* beam_smear_fn = new TF1("beam_smear_fn", "(1./x)", ebeam_min, ebeam_max);
beam_smear->SetReaction("g + p");
beam_smear->SetMomentumFunction(beam_smear_fn);
makeDistributionManager()->Add(beam_smear);
PReaction *rr=new PReaction("_P1 = 2.2","g","p","p pi0 [dilepton [e+ e-] g]","ppi0");
rr->Do(histo2," _x = ([dilepton])->M() ");
rr->Print();
rr->loop(125000);
histo2->Draw();
}
however, if I try to use the line
makeStaticData()->SetParticleMass("dilepton",0.06);
When using a model i.e.
#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 "TFile.h"
#include "/w/hallb/clasg12/mkunkel/PLUTO/pluto_v5.42/src/PParticle.h"
#include "/w/hallb/clasg12/mkunkel/PLUTO/pluto_v5.42/src/PReaction.h"
#include "/w/hallb/clasg12/mkunkel/PLUTO/pluto_v5.42/src/PBeamSmearing.h"
#include "/w/hallb/clasg12/mkunkel/PLUTO/pluto_v5.42/src/PAnyDistribution.h"
#include "/w/hallb/clasg12/mkunkel/PLUTO/pluto_v5.42/plugins/scatter_mod/PScatterCrossSection.h"
void pi0_XSection(int input, int events) //int input, int events
{
//Open TFile and get histograms
TFile *fhist = new TFile("PI0_Hists_SAID.root");
//TFile *f2hist = new TFile("/Volumes/Mac_Storage/Work_Data/DALITZ_SIMULATION/crosssectionData/Eta_CrossSection/diff_Xsection/TriggerHist.root");
//##################LOAD HISTOGRAMS HERE###################
TH1D *lab_profile = (TH1D*)fhist->Get("Blab_smear");
TH1D *cm_profile = (TH1D*)fhist->Get("B_smear");
//TH1D *trigger_profile = (TH1D*)f2hist->Get("trigger4_beam");
//TH1D *example0 = (TH1D*)fhist->Get("hist0");
TH1D *example1 = (TH1D*)fhist->Get("hist1");
TH1D *example2 = (TH1D*)fhist->Get("hist2");
TH1D *example3 = (TH1D*)fhist->Get("hist3");
TH1D *example4 = (TH1D*)fhist->Get("hist4");
TH1D *example5 = (TH1D*)fhist->Get("hist5");
TH1D *example6 = (TH1D*)fhist->Get("hist6");
TH1D *example7 = (TH1D*)fhist->Get("hist7");
TH1D *example8 = (TH1D*)fhist->Get("hist8");
TH1D *example9 = (TH1D*)fhist->Get("hist9");
TH1D *example10 = (TH1D*)fhist->Get("hist10");
TH1D *example11 = (TH1D*)fhist->Get("hist11");
TH1D *example12 = (TH1D*)fhist->Get("hist12");
TH1D *example13 = (TH1D*)fhist->Get("hist13");
TH1D *example14 = (TH1D*)fhist->Get("hist14");
TH1D *example15 = (TH1D*)fhist->Get("hist15");
TH1D *example16 = (TH1D*)fhist->Get("hist16");
TH1D *example17 = (TH1D*)fhist->Get("hist17");
TH1D *example18 = (TH1D*)fhist->Get("hist18");
TH1D *example19 = (TH1D*)fhist->Get("hist19");
TH1D *example20 = (TH1D*)fhist->Get("hist20");
TH1D *example21 = (TH1D*)fhist->Get("hist21");
TH1D *example22 = (TH1D*)fhist->Get("hist22");
TH1D *example23 = (TH1D*)fhist->Get("hist23");
TH1D *example24 = (TH1D*)fhist->Get("hist24");
TH1D *example25 = (TH1D*)fhist->Get("hist25");
TH1D *example26 = (TH1D*)fhist->Get("hist26");
TH1D *example27 = (TH1D*)fhist->Get("hist27");
TH1D *example28 = (TH1D*)fhist->Get("hist28");
TH1D *example29 = (TH1D*)fhist->Get("hist29");
TH1D *example30 = (TH1D*)fhist->Get("hist30");
TH1D *example31 = (TH1D*)fhist->Get("hist31");
TH1D *example32 = (TH1D*)fhist->Get("hist32");
TH1D *example33 = (TH1D*)fhist->Get("hist33");
TH1D *example34 = (TH1D*)fhist->Get("hist34");
TH1D *example35 = (TH1D*)fhist->Get("hist35");
TH1D *example36 = (TH1D*)fhist->Get("hist36");
TH1D *example37 = (TH1D*)fhist->Get("hist37");
TH1D *example38 = (TH1D*)fhist->Get("hist38");
TH1D *example39 = (TH1D*)fhist->Get("hist39");
TH1D *example40 = (TH1D*)fhist->Get("hist40");
TH1D *example41 = (TH1D*)fhist->Get("hist41");
TH1D *example42 = (TH1D*)fhist->Get("hist42");
TH1D *example43 = (TH1D*)fhist->Get("hist43");
TH1D *example44 = (TH1D*)fhist->Get("hist44");
TH1D *example45 = (TH1D*)fhist->Get("hist45");
TH1D *example46 = (TH1D*)fhist->Get("hist46");
TH1D *example47 = (TH1D*)fhist->Get("hist47");
TH1D *example48 = (TH1D*)fhist->Get("hist48");
TH1D *example49 = (TH1D*)fhist->Get("hist49");
TH1D *example50 = (TH1D*)fhist->Get("hist50");
TH1D *example51 = (TH1D*)fhist->Get("hist51");
TH1D *example52 = (TH1D*)fhist->Get("hist52");
TH1D *example53 = (TH1D*)fhist->Get("hist53");
TH1D *example54 = (TH1D*)fhist->Get("hist54");
TH1D *example55 = (TH1D*)fhist->Get("hist55");
TH1D *example56 = (TH1D*)fhist->Get("hist56");
TH1D *example57 = (TH1D*)fhist->Get("hist57");
TH1D *example58 = (TH1D*)fhist->Get("hist58");
TH1D *example59 = (TH1D*)fhist->Get("hist59");
TH1D *example60 = (TH1D*)fhist->Get("hist60");
TH1D *example61 = (TH1D*)fhist->Get("hist61");
TH1D *example62 = (TH1D*)fhist->Get("hist62");
TH1D *example63 = (TH1D*)fhist->Get("hist63");
TH1D *example64 = (TH1D*)fhist->Get("hist64");
TH1D *example65 = (TH1D*)fhist->Get("hist65");
TH1D *example66 = (TH1D*)fhist->Get("hist66");
TH1D *example67 = (TH1D*)fhist->Get("hist67");
TH1D *example68 = (TH1D*)fhist->Get("hist68");
TH1D *example69 = (TH1D*)fhist->Get("hist69");
//##################END HISTOGRAMS HERE###################//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("pi0,daughter,primary");
//Define the range of the c.m. sampling: G12 c.m. beam energy 1.71565, 2.55105 is 1.1 - 3.0 GeV in lab
double q_min = 1.71565;//1.71565; //;
double q_max = 2.55105;// 2.7//;
model->SetRange(q_min,q_max);
//These Histograms are for the model
//These Histograms are for the model
//model->AddHistogram(example0,"if (_y > 1.69534 && _y <= 1.70912) value = Eval(_x); _f =value");
model->AddHistogram(example1,"if (_y > 1.70912 && _y <= 1.72279) value = Eval(_x); _f =value");
model->AddHistogram(example2,"if (_y > 1.72279 && _y <= 1.73635) value = Eval(_x); _f =value");
model->AddHistogram(example3,"if (_y > 1.73635 && _y <= 1.74981) value = Eval(_x); _f =value");
model->AddHistogram(example4,"if (_y > 1.74981 && _y <= 1.76316) value = Eval(_x); _f =value");
model->AddHistogram(example5,"if (_y > 1.76316 && _y <= 1.77642) value = Eval(_x); _f =value");
model->AddHistogram(example6,"if (_y > 1.77642 && _y <= 1.78957) value = Eval(_x); _f =value");
model->AddHistogram(example7,"if (_y > 1.78957 && _y <= 1.80263) value = Eval(_x); _f =value");
model->AddHistogram(example8,"if (_y > 1.80263 && _y <= 1.8156) value = Eval(_x); _f =value");
model->AddHistogram(example9,"if (_y > 1.8156 && _y <= 1.82847) value = Eval(_x); _f =value");
model->AddHistogram(example10,"if (_y > 1.82847 && _y <= 1.84126) value = Eval(_x); _f =value");
model->AddHistogram(example11,"if (_y > 1.84126 && _y <= 1.85395) value = Eval(_x); _f =value");
model->AddHistogram(example12,"if (_y > 1.85395 && _y <= 1.86656) value = Eval(_x); _f =value");
model->AddHistogram(example13,"if (_y > 1.86656 && _y <= 1.87909) value = Eval(_x); _f =value");
model->AddHistogram(example14,"if (_y > 1.87909 && _y <= 1.89153) value = Eval(_x); _f =value");
model->AddHistogram(example15,"if (_y > 1.89153 && _y <= 1.90389) value = Eval(_x); _f =value");
model->AddHistogram(example16,"if (_y > 1.90389 && _y <= 1.91617) value = Eval(_x); _f =value");
model->AddHistogram(example17,"if (_y > 1.91617 && _y <= 1.92837) value = Eval(_x); _f =value");
model->AddHistogram(example18,"if (_y > 1.92837 && _y <= 1.9405) value = Eval(_x); _f =value");
model->AddHistogram(example19,"if (_y > 1.9405 && _y <= 1.95255) value = Eval(_x); _f =value");
model->AddHistogram(example20,"if (_y > 1.95255 && _y <= 1.96453) value = Eval(_x); _f =value");
model->AddHistogram(example21,"if (_y > 1.96453 && _y <= 1.97643) value = Eval(_x); _f =value");
model->AddHistogram(example22,"if (_y > 1.97643 && _y <= 1.98826) value = Eval(_x); _f =value");
model->AddHistogram(example23,"if (_y > 1.98826 && _y <= 2.00003) value = Eval(_x); _f =value");
model->AddHistogram(example24,"if (_y > 2.00003 && _y <= 2.01172) value = Eval(_x); _f =value");
model->AddHistogram(example25,"if (_y > 2.01172 && _y <= 2.02335) value = Eval(_x); _f =value");
model->AddHistogram(example26,"if (_y > 2.02335 && _y <= 2.03491) value = Eval(_x); _f =value");
model->AddHistogram(example27,"if (_y > 2.03491 && _y <= 2.0464) value = Eval(_x); _f =value");
model->AddHistogram(example28,"if (_y > 2.0464 && _y <= 2.05783) value = Eval(_x); _f =value");
model->AddHistogram(example29,"if (_y > 2.05783 && _y <= 2.0692) value = Eval(_x); _f =value");
model->AddHistogram(example30,"if (_y > 2.0692 && _y <= 2.08051) value = Eval(_x); _f =value");
model->AddHistogram(example31,"if (_y > 2.08051 && _y <= 2.09175) value = Eval(_x); _f =value");
model->AddHistogram(example32,"if (_y > 2.09175 && _y <= 2.10293) value = Eval(_x); _f =value");
model->AddHistogram(example33,"if (_y > 2.10293 && _y <= 2.11406) value = Eval(_x); _f =value");
model->AddHistogram(example34,"if (_y > 2.11406 && _y <= 2.12513) value = Eval(_x); _f =value");
model->AddHistogram(example35,"if (_y > 2.12513 && _y <= 2.13613) value = Eval(_x); _f =value");
model->AddHistogram(example36,"if (_y > 2.13613 && _y <= 2.14709) value = Eval(_x); _f =value");
model->AddHistogram(example37,"if (_y > 2.14709 && _y <= 2.15798) value = Eval(_x); _f =value");
model->AddHistogram(example38,"if (_y > 2.15798 && _y <= 2.16883) value = Eval(_x); _f =value");
model->AddHistogram(example39,"if (_y > 2.16883 && _y <= 2.17962) value = Eval(_x); _f =value");
model->AddHistogram(example40,"if (_y > 2.17962 && _y <= 2.19035) value = Eval(_x); _f =value");
model->AddHistogram(example41,"if (_y > 2.19035 && _y <= 2.20103) value = Eval(_x); _f =value");
model->AddHistogram(example42,"if (_y > 2.20103 && _y <= 2.21167) value = Eval(_x); _f =value");
model->AddHistogram(example43,"if (_y > 2.21167 && _y <= 2.22225) value = Eval(_x); _f =value");
model->AddHistogram(example44,"if (_y > 2.22225 && _y <= 2.23278) value = Eval(_x); _f =value");
model->AddHistogram(example45,"if (_y > 2.23278 && _y <= 2.24326) value = Eval(_x); _f =value");
model->AddHistogram(example46,"if (_y > 2.24326 && _y <= 2.25369) value = Eval(_x); _f =value");
model->AddHistogram(example47,"if (_y > 2.25369 && _y <= 2.26407) value = Eval(_x); _f =value");
model->AddHistogram(example48,"if (_y > 2.26407 && _y <= 2.27441) value = Eval(_x); _f =value");
model->AddHistogram(example49,"if (_y > 2.27441 && _y <= 2.2847) value = Eval(_x); _f =value");
model->AddHistogram(example50,"if (_y > 2.2847 && _y <= 2.29495) value = Eval(_x); _f =value");
model->AddHistogram(example51,"if (_y > 2.29495 && _y <= 2.30514) value = Eval(_x); _f =value");
model->AddHistogram(example52,"if (_y > 2.30514 && _y <= 2.3153) value = Eval(_x); _f =value");
model->AddHistogram(example53,"if (_y > 2.3153 && _y <= 2.32541) value = Eval(_x); _f =value");
model->AddHistogram(example54,"if (_y > 2.32541 && _y <= 2.33547) value = Eval(_x); _f =value");
model->AddHistogram(example55,"if (_y > 2.33547 && _y <= 2.34549) value = Eval(_x); _f =value");
model->AddHistogram(example56,"if (_y > 2.34549 && _y <= 2.35547) value = Eval(_x); _f =value");
model->AddHistogram(example57,"if (_y > 2.35547 && _y <= 2.36541) value = Eval(_x); _f =value");
model->AddHistogram(example58,"if (_y > 2.36541 && _y <= 2.37531) value = Eval(_x); _f =value");
model->AddHistogram(example59,"if (_y > 2.37531 && _y <= 2.38516) value = Eval(_x); _f =value");
model->AddHistogram(example60,"if (_y > 2.38516 && _y <= 2.39498) value = Eval(_x); _f =value");
model->AddHistogram(example61,"if (_y > 2.39498 && _y <= 2.40475) value = Eval(_x); _f =value");
model->AddHistogram(example62,"if (_y > 2.40475 && _y <= 2.41449) value = Eval(_x); _f =value");
model->AddHistogram(example63,"if (_y > 2.41449 && _y <= 2.42418) value = Eval(_x); _f =value");
model->AddHistogram(example64,"if (_y > 2.42418 && _y <= 2.43384) value = Eval(_x); _f =value");
model->AddHistogram(example65,"if (_y > 2.43384 && _y <= 2.44346) value = Eval(_x); _f =value");
model->AddHistogram(example66,"if (_y > 2.44346 && _y <= 2.45304) value = Eval(_x); _f =value");
model->AddHistogram(example67,"if (_y > 2.45304 && _y <= 2.46258) value = Eval(_x); _f =value");
model->AddHistogram(example68,"if (_y > 2.46258 && _y <= 2.47209) value = Eval(_x); _f =value");
//model->AddHistogram(example69,"if (_y > 2.47209 && _y <= 2.48156) value = Eval(_x); _f =value");
model->AddHistogram(example69,"if (_y > 2.47209 && _y <= 2.7) value = Eval(_x); _f =value");
//These Histograms are for the beam smear. There are 2 methods, 1 method for each histogram
//Use only one method
//Method 1: Beam profile is in c.m. and distributed as 1/E
//Just evaluate histogram
//model->AddEquation("b_value = (_y);");
//model->AddHistogram(cm_profile,"beam_value = Eval(b_value);");
//Method 2: Beam profile is in lab frame (t_lab, as in Ingo's notation from his example macro)
//Must convert this profile to c.m. ie. t_lab = (_y*_y - g.mass*g.mass - p.mass*p.mass)/(2*p.mass*p.mass) - g.mass;
model->AddEquation("t_lab = (_y*_y - g.mass*g.mass - p.mass*p.mass)/(2*p.mass*p.mass) - g.mass;");
model->AddHistogram(lab_profile,"_f = _f * Eval(t_lab);");
//model->SetNpy(300);
//model->SetNpx(300); //Use this after SetNpy analysis
makeDistributionManager()->Add(model);
// TCanvas *q_cos = new TCanvas("q_cos","q_cos");
// q_cos->cd();
// model->GetFunction()->Draw("surf2");
//Set up reaction as usual now
((PDalitzDecay * )makeDistributionManager()->GetDistribution("pi0_dalitz"))->SetUseQED(1);
//TH1F * histo2 = new TH1F ("histo2","IM e+ e-",80,0.58,0.146);
//gROOT->Reset();
PUtils::SetSeed(input*1000); //input
makeStaticData()->SetParticleMass("dilepton",0.06);
PReaction my_reaction("_P1 = 2.2","g","p","p pi0 [dilepton [e+ e-] g]","pi0_dalitz_QEDFF",1,0,0,0);
//my_reaction.Do(histo2," _x = ([dilepton])->M() ");
my_reaction.Print();
my_reaction.Loop(events);
//histo2->Draw();
}
This does not sample the dilepton mass as the previous code.
Could I please be informed on how to sample the dilepton mass when using a model.
Thanks
Michael
[Updated on: Tue, 30 July 2013 20:58] Report message to a moderator
|
|
|
Goto Forum:
Current Time: Sun Sep 01 19:05:50 CEST 2024
Total time taken to generate the page: 0.00881 seconds
|