Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
makePlot_truth_xQ2_Studies.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file makePlot_truth_xQ2_Studies.C
1 
9 int makePlot_truth_xQ2_Studies(
10  const TString infile = "../../data/Sample_DISReco_ep.root",
11  //const TString infile = "../../data/EventGenAna_Pythia6_DIS_10x250_100k.root",
12  const bool save_figures = true
13  )
14 
15 {
16  gStyle->SetOptStat(kFALSE);
17 
18  /*--------------Get Input File--------------*/
19 
20  TFile *f_truth = new TFile( infile, "OPEN");
21  TTree* t_truth = (TTree*)f_truth->Get("event_truth");
22  t_truth->Print();
23 
24  /*------------------------------------------*/
25  /*---------------Define Cuts----------------*/
26 
27  TCut electron_cut = "em_evtgen_pid == 11";
28 
29  /*------------------------------------------*/
30  /*---------x-Q2 Space for Electrons---------*/
31  /*------------------------------------------*/
32 
33  /*------------------------------------------*/
34  /*------------Create Histogram--------------*/
35  /*------------------------------------------*/
39  /*------------------------------------------*/
40 
41  TH2F* hQ2x_e = new TH2F("hQ2x_e", "", 40, -4, 0, 60, 0, 3);
42  TAxis *xaxis = hQ2x_e->GetXaxis();
43  TAxis *yaxis = hQ2x_e->GetYaxis();
44  int xbins = xaxis->GetNbins();
45  int ybins = yaxis->GetNbins();
46 
47  Axis_t xmin = xaxis->GetXmin();
48  Axis_t xmax = xaxis->GetXmax();
49  Axis_t xwidth = (xmax - xmin ) / xbins;
50  Axis_t *new_xbins = new Axis_t[xbins + 1];
51 
52  for( int i =0; i <= xbins; i++)
53  {
54  new_xbins[i] = TMath::Power( 10, xmin + i*xwidth);
55  }
56  xaxis->Set(xbins, new_xbins);
57 
58  Axis_t ymin = yaxis->GetXmin();
59  Axis_t ymax = yaxis->GetXmax();
60  Axis_t ywidth = (ymax - ymin) / ybins;
61  Axis_t *new_ybins = new Axis_t[ybins + 1];
62 
63  for( int i2 =0; i2 <= ybins; i2++)
64  {
65  new_ybins[i2] = TMath::Power( 10, ymin + i2*ywidth);
66  }
67 
68  yaxis->Set(ybins, new_ybins);
69  /*------------------------------------------*/
70  /*-------------Plot x-Q2 Space--------------*/
71  /*------------------------------------------*/
72 
73  /*------------------------------------------*/
74  /*----------Scattered Lepton Angle----------*/
75  /*------------------------------------------*/
76 
77  TCanvas *cQ2x_e1 = new TCanvas("cQ2x_e1");
78  cQ2x_e1->SetLogx();
79  cQ2x_e1->SetLogy();
80  t_truth->Draw("evtgen_Q2:evtgen_x>>hQ2x_e", "" , "colz");
81 
82  hQ2x_e->GetXaxis()->SetTitle("x");
83  hQ2x_e->GetYaxis()->SetTitle("Q^{2} [GeV^{2}]");
84 
85  TF1 *f_LepA_0 = new TF1("f_LepA_0", "2*[0]*[1]*x*( 1 + cos([2])) / (1 + ((x*[1])/(2*[0]))*(1-cos([2])) - 0.5 * (1 - cos([2])))", 10e-5, 1);
86  f_LepA_0->SetParameter( 0, 10);
87  f_LepA_0->SetParameter( 1, 250);
88  f_LepA_0->SetParameter( 2, (0 * TMath::Pi())/180 );
89  TF1 *f_LepA_90 = (TF1*)f_LepA_0->Clone("f_LepA_90");
90  f_LepA_90->SetParameter( 2, (90 * TMath::Pi())/180 );
91  TF1 *f_LepA_150 = (TF1*)f_LepA_0->Clone("f_LepA_150");
92  f_LepA_150->SetParameter( 2, (150 * TMath::Pi())/180 );
93  TF1 *f_LepA_170 = (TF1*)f_LepA_0->Clone("f_LepA_170");
94  f_LepA_170->SetParameter( 2, (170 * TMath::Pi())/180 );
95 
96  f_LepA_0->Draw("SAME");
97  f_LepA_0->SetLineColor(1);
98  f_LepA_90->Draw("SAME");
99  f_LepA_90->SetLineColor(1);
100  f_LepA_150->Draw("SAME");
101  f_LepA_150->SetLineColor(1);
102  f_LepA_170->Draw("SAME");
103  f_LepA_170->SetLineColor(1);
104 
105  /*---------------Label Lines----------------*/
106 
107  TPaveText *pt_LepA_0 = new TPaveText(0.4, 1.2, 0.6, 1.4);
108  pt_LepA_0->AddText("0^{/circ}");
109  pt_LepA_0->Draw();
110 
111  /*------------------------------------------*/
112  /*---------Scattered Lepton Energy----------*/
113  /*------------------------------------------*/
114 
115  TCanvas *cQ2x_e2 = new TCanvas("cQ2x_e2");
116  cQ2x_e2->SetLogx();
117  cQ2x_e2->SetLogy();
118  t_truth->Draw("evtgen_Q2:evtgen_x>>hQ2x_e", "" , "colz");
119 
120  hQ2x_e->GetXaxis()->SetTitle("x");
121  hQ2x_e->GetYaxis()->SetTitle("Q^{2} [GeV^{2}]");
122 
123  TF1 *f_LepE_5 = new TF1("f_LepE_5", "2*[0]*[2]*( 1 + ((((x*[1]*([2]-[0]))/(x*[1]-[0])) - ([0] - (([0]*([2]-[0]))/(x*[1]-[0])))) / (((x*[1]*([2]-[0]))/(x*[1]-[0])) + ([0] - (([0]*([2]-[0]))/(x*[1]-[0]))))))", 10e-5, 1);
124  f_LepE_5->SetParameter( 0, 10);
125  f_LepE_5->SetParameter( 1, 250);
126  f_LepE_5->SetParameter( 2, 5);
127  TF1 *f_LepE_9 = (TF1*)f_LepE_5->Clone("f_LepE_9");
128  f_LepE_9->SetParameter( 2, 9);
129  TF1 *f_LepE_11 = (TF1*)f_LepE_5->Clone("f_LepE_11");
130  f_LepE_11->SetParameter( 2, 11);
131  TF1 *f_LepE_15 = (TF1*)f_LepE_5->Clone("f_LepE_15");
132  f_LepE_15->SetParameter( 2, 15);
133 
134  f_LepE_5->Draw("SAME");
135  f_LepE_5->SetLineColor(1);
136  f_LepE_9->Draw("SAME");
137  f_LepE_9->SetLineColor(1);
138  f_LepE_11->Draw("SAME");
139  f_LepE_11->SetLineColor(1);
140  f_LepE_15->Draw("SAME");
141  f_LepE_15->SetLineColor(1);
142 
143  /*------------------------------------------*/
144  /*------------Current Jet Angle-------------*/
145  /*------------------------------------------*/
146 
147  TCanvas *cQ2x_e3 = new TCanvas("cQ2x_e3");
148  cQ2x_e3->SetLogx();
149  cQ2x_e3->SetLogy();
150  t_truth->Draw("evtgen_Q2:evtgen_x>>hQ2x_e", "" , "colz");
151 
152  hQ2x_e->GetXaxis()->SetTitle("x");
153  hQ2x_e->GetYaxis()->SetTitle("Q^{2} [GeV^{2}]");
154 
155  TF1 *f_JetA_180 = new TF1("f_JetA_180" ,"(4 * (x**2)*([1]**2)*[0]*(1-cos([2]))) / (cos([2]) * ([0] -x*[1]) + [0] + x*[1])" , 10e-6, 1);
156  f_JetA_180->SetParameter( 0 , 10);
157  f_JetA_180->SetParameter( 1 , 250.);
158  f_JetA_180->SetParameter( 2 , (180 * TMath::Pi())/180 );
159  TF1 *f_JetA_150 = (TF1*)f_JetA_180->Clone("f_JetA_150");
160  f_JetA_150->SetParameter( 2 , (150 * TMath::Pi())/180 );
161  TF1 *f_JetA_30 = (TF1*)f_JetA_180->Clone("f_JetA_30");
162  f_JetA_30->SetParameter( 2 , (30 * TMath::Pi())/180 );
163  TF1 *f_JetA_10 = (TF1*)f_JetA_180->Clone("f_JetA_10");
164  f_JetA_10->SetParameter( 2 , (10 * TMath::Pi())/180 );
165 
166  f_JetA_180->Draw("SAME");
167  f_JetA_180->SetLineColor(1);
168  f_JetA_150->Draw("SAME");
169  f_JetA_150->SetLineColor(1);
170  f_JetA_30->Draw("SAME");
171  f_JetA_30->SetLineColor(1);
172  f_JetA_10->Draw("SAME");
173  f_JetA_10->SetLineColor(1);
174 
175  /*------------------------------------------*/
176  /*------------Current Jet Energy------------*/
177  /*------------------------------------------*/
178 
179  TCanvas *cQ2x_e4 = new TCanvas("cQ2x_e4");
180  cQ2x_e4->SetLogx();
181  cQ2x_e4->SetLogy();
182  t_truth->Draw("evtgen_Q2:evtgen_x>>hQ2x_e", "" , "colz");
183 
184  hQ2x_e->GetXaxis()->SetTitle("x");
185  hQ2x_e->GetYaxis()->SetTitle("Q^{2} [GeV^{2}]");
186 
187  TF1 *f_JetE_5 = new TF1("f_JetE_5" , "(4*x*[0]*[1]*([2]-x*[1])) / ([0] - x*[1])", 10e-6, 1);
188  f_JetE_5->SetParameter( 0 , 10);
189  f_JetE_5->SetParameter( 1 , 250.);
190  f_JetE_5->SetParameter( 2 , 5. );
191  TF1 *f_JetE_7 = (TF1*)f_JetE_5->Clone("f_JetE_7");
192  f_JetE_7->SetParameter( 2 , 7);
193  TF1 *f_JetE_10 = (TF1*)f_JetE_5->Clone("f_JetE_10");
194  f_JetE_10->SetParameter( 2 , 10);
195  TF1 *f_JetE_20 = (TF1*)f_JetE_5->Clone("f_JetE_20");
196  f_JetE_20->SetParameter( 2 , 20);
197 
198  f_JetE_5->Draw("SAME");
199  f_JetE_5->SetLineColor(1);
200  f_JetE_7->Draw("SAME");
201  f_JetE_7->SetLineColor(1);
202  f_JetE_10->Draw("SAME");
203  f_JetE_10->SetLineColor(1);
204  f_JetE_20->Draw("SAME");
205  f_JetE_20->SetLineColor(1);
206 
207  /*------------------------------------------*/
208  /*---------------Inelasticity---------------*/
209  /*------------------------------------------*/
210 
211  TCanvas *cQ2x_e5 = new TCanvas("cQ2x_e5");
212  cQ2x_e5->SetLogx();
213  cQ2x_e5->SetLogy();
214  t_truth->Draw("evtgen_Q2:evtgen_x>>hQ2x_e", "" , "colz");
215 
216  hQ2x_e->GetXaxis()->SetTitle("x");
217  hQ2x_e->GetYaxis()->SetTitle("Q^{2} [GeV^{2}]");
218 
219  TF1 *f_y1 = new TF1("f_y1", "4*x*[0]*[1]*[2]", 10e-5, 1);
220  f_y1->SetParameter( 0, 10);
221  f_y1->SetParameter( 1, 250);
222  f_y1->SetParameter( 2, 1);
223  TF1 *f_y01 = (TF1*)f_y1->Clone("f_y01");
224  f_y01->SetParameter(2 , 0.1);
225  TF1 *f_y001 = (TF1*)f_y1->Clone("f_y001");
226  f_y001->SetParameter(2 , 0.01);
227  TF1 *f_y0001 = (TF1*)f_y1->Clone("f_y0001");
228  f_y0001->SetParameter(2 , 0.001);
229 
230  f_y1->Draw("SAME");
231  f_y1->SetLineColor(1);
232  f_y01->Draw("SAME");
233  f_y01->SetLineColor(1);
234  f_y001->Draw("SAME");
235  f_y001->SetLineColor(1);
236  f_y0001->Draw("SAME");
237  f_y0001->SetLineColor(1);
238 
239  /*------------------------------------------*/
240  /*-----------Final Hadronic Mass------------*/
241  /*------------------------------------------*/
242 
243  TCanvas *cQ2x_e6 = new TCanvas("cQ2x_e6");
244  cQ2x_e6->SetLogx();
245  cQ2x_e6->SetLogy();
246  t_truth->Draw("evtgen_Q2:evtgen_x>>hQ2x_e", "" , "colz");
247 
248  hQ2x_e->GetXaxis()->SetTitle("x");
249  hQ2x_e->GetYaxis()->SetTitle("Q^{2} [GeV^{2}]");
250 
251  TF1 *f_W10 = new TF1("f_W10" , "([1] - [0]**2) / ((1/x) - 1)", 10e-7, 1);
252  f_W10->SetParameter( 0 , 0.938);
253  f_W10->SetParameter( 1 , (10)**2 );
254  TF1 *f_W20 = (TF1*)f_W10->Clone("f_W20");
255  f_W20->SetParameter( 1 , (20)**2 );
256  TF1 *f_W50 = (TF1*)f_W10->Clone("f_W50");
257  f_W50->SetParameter( 1 , (50)**2 );
258  TF1 *f_W100 = (TF1*)f_W10->Clone("f_W100");
259  f_W100->SetParameter( 1 , (100)**2 );
260  TF1 *f_W3 = (TF1*)f_W10->Clone("f_W3");
261  f_W3->SetParameter( 1 , (3)**2 );
262 
263  f_W10->Draw("SAME");
264  f_W10->SetLineColor(1);
265  f_W20->Draw("SAME");
266  f_W20->SetLineColor(1);
267  f_W50->Draw("SAME");
268  f_W50->SetLineColor(1);
269  f_W100->Draw("SAME");
270  f_W100->SetLineColor(1);
271  f_W3->Draw("SAME");
272  f_W3->SetLineColor(1);
273 
274  if ( save_figures )
275  {
276  cQ2x_e1->Print("plots/Pythia_Q2x_LepA_e.eps");
277  cQ2x_e2->Print("plots/Pythia_Q2x_LepE_e.eps");
278  cQ2x_e3->Print("plots/Pythia_Q2x_JetA_e.eps");
279  cQ2x_e4->Print("plots/Pythia_Q2x_JetE_e.eps");
280  cQ2x_e5->Print("plots/Pythia_Q2x_y_e.eps");
281  cQ2x_e6->Print("plots/Pythia_Q2x_W_e.eps");
282  }
283 
284  return 0;
285 
286 }