Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
plot_amatrix.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file plot_amatrix.C
1 #include "TraceBox.h"
2 #include <cmath>
3 int
4 plot_amatrix()
5 {
6  /**** Chose input files ****/
7 
8  /* Open iput file with trajectories from GEANT4 */
9  TFile *fin_0 = new TFile( "data_amatrix/G4EICIR_DSTReader_p250_vx0_vy0_theta22mrad.root", "OPEN" );
10  TFile *fin_a11 = new TFile( "data_amatrix/G4EICIR_DSTReader_p250_vx0_vy1_theta0mrad.root", "OPEN" );
11  TFile *fin_a21 = new TFile( "data_amatrix/G4EICIR_DSTReader_p250_vx0_vy0_theta10mrad.root", "OPEN" );
12  TFile *fin_a12 = new TFile( "data_amatrix/G4EICIR_DSTReader_p250_vx0_vy1_theta0mrad.root", "OPEN" );
13  TFile *fin_a22 = new TFile( "data_amatrix/G4EICIR_DSTReader_p250_vx0_vy0_theta10mrad.root", "OPEN" );
14 
15  /* Get tree from file */
16  TTree *tin_0 = (TTree*)fin_0->Get("T");
17  TTree *tin_a11 = (TTree*)fin_a11->Get("T");
18  TTree *tin_a21 = (TTree*)fin_a21->Get("T");
19  TTree *tin_a12 = (TTree*)fin_a12->Get("T");
20  TTree *tin_a22 = (TTree*)fin_a22->Get("T");
21 
22  /* define initial values for read-in trajectories */
23  float tin_a11_vy = 1.0; // this is y_0
24  float tin_a21_theta = 10.0; // this is theta_y*
25  float tin_a12_vy = 1.0; // this is y_0
26  float tin_a22_theta = 10.0; // theta_y*
27 
28  /* get number of hits */
29  int nhits = 0;
30  tin_a11->SetBranchAddress("n_G4HIT_FWDDISC",&nhits);
31  tin_a11->GetEntry(0);
32  cout << "hits: " << nhits << endl;
33 
34  /* create graph of nominal beam position (z) */
35  tin_0->Draw( "0.5*(G4HIT_FWDDISC.y[][0]+G4HIT_FWDDISC.y[][1]) : 0.5*(G4HIT_FWDDISC.z[][0]+G4HIT_FWDDISC.z[][1])","Entry$==0","");
36  TGraph* g_0 = new TGraph(nhits, &(tin_0->GetV2()[0]), &(tin_0->GetV1()[0]));
37 
38  /* create graph of a11(z) */
39  // tin_a11->Draw(TString::Format( "0.5*(G4HIT_FWDDISC.y[][0]+G4HIT_FWDDISC.y[][1]) / %f : 0.5*(G4HIT_FWDDISC.z[][0]+G4HIT_FWDDISC.z[][1])", tin_a11_vy ),"Entry$==0","");
40 
41  /* a11 = y/y_0 */
42  tin_a11->Draw("0.5*(G4HIT_FWDDISC.y[][0]+G4HIT_FWDDISC.y[][1]) / 1 : 0.5*(G4HIT_FWDDISC.z[][0]+G4HIT_FWDDISC.z[][1])","Entry$==0","" );
43  TGraph* g_a11 = new TGraph(nhits, &(tin_a11->GetV2()[0]), &(tin_a11->GetV1()[0]));
44  g_a11->SetMarkerStyle(7);
45  g_a11->SetMarkerSize(1);
46  g_a11->SetMarkerColor(kRed);
47 
48  /* create graph of a21(z) */
49  //tin_a21->Draw(TString::Format( "0.5*(G4HIT_FWDDISC.y[][0]+G4HIT_FWDDISC.y[][1]) / %f : 0.5*(G4HIT_FWDDISC.z[][0]+G4HIT_FWDDISC.z[][1])", tin_a21_theta ),"Entry$==0","");
50 
51 /* a21 = y/theta* */
52  tin_a21->Draw( "0.5*(G4HIT_FWDDISC.y[][0]+G4HIT_FWDDISC.y[][1]) / 10 : 0.5*(G4HIT_FWDDISC.z[][0]+G4HIT_FWDDISC.z[][1])","Entry$==0","");
53  TGraph* g_a21 = new TGraph(nhits, &(tin_a21->GetV2()[0]), &(tin_a21->GetV1()[0]));
54  g_a21->SetMarkerStyle(7);
55  g_a21->SetMarkerSize(1);
56  g_a21->SetMarkerColor(kBlue);
57 
58  /* create graph of a12(z) */
59  // tin_a12->Draw(TString::Format( "0.5*(G4HIT_FWDDISC.y[][0]+G4HIT_FWDDISC.y[][1]) / %f : 0.5*(G4HIT_FWDDISC.z[][0]+G4HIT_FWDDISC.z[][1])", tin_a12_vy ),"Entry$==0","");
60 
61  /* a12 = theta_y/y_0, theta*=atan( |exit-entry| / thickness ), here thickness = 1 */
62  tin_a12->Draw("atan( abs(G4HIT_FWDDISC.y[][0]-G4HIT_FWDDISC.y[][1]) ) / 1 : 0.5*(G4HIT_FWDDISC.z[][0]+G4HIT_FWDDISC.z[][1])","Entry$==0","");
63  TGraph* g_a12 = new TGraph(nhits, &(tin_a12->GetV2()[0]), &(tin_a12->GetV1()[0]));
64  g_a12->SetMarkerStyle(7);
65  g_a12->SetMarkerSize(1);
66  g_a12->SetMarkerColor(kGreen);
67 
68  /* create graph of a22(z) */
69  // tin_a22->Draw(TString::Format( "0.5*(G4HIT_FWDDISC.y[][0]+G4HIT_FWDDISC.y[][1]) / %f : 0.5*(G4HIT_FWDDISC.z[][0]+G4HIT_FWDDISC.z[][1])", tin_a22_theta ),"Entry$==0","");
70 
71  /* a22 = thetay/thetay*, theta*=atan( |exit-entry| / thickness ), here thickness = 1 */
72  tin_a22->Draw("atan( abs(G4HIT_FWDDISC.y[][0]-G4HIT_FWDDISC.y[][1]) ) / 10 : 0.5*(G4HIT_FWDDISC.z[][0]+G4HIT_FWDDISC.z[][1])","Entry$==0","");
73  TGraph* g_a22 = new TGraph(nhits, &(tin_a22->GetV2()[0]), &(tin_a22->GetV1()[0]));
74  g_a22->SetMarkerStyle(7);
75  g_a22->SetMarkerSize(1);
76  g_a22->SetMarkerColor(kBlue);
77  g_a22->SetMarkerColor(kOrange);
78 
79  /* Create Legend */
80  TLegend* leg = new TLegend(0.3,0.7,0.6,0.9);
81  leg->AddEntry(g_0,"Nominal beam","l");
82  leg->AddEntry(g_a11,"a11","p");
83  leg->AddEntry(g_a21,"a21","p");
84  leg->AddEntry(g_a12,"a12","p");
85  leg->AddEntry(g_a22,"a22","p");
86 
87  /* Create frame histogram for plot */
88  TH1F *h1 = new TH1F("h1","",10,0,15000);
89  h1->GetXaxis()->SetRangeUser(0,12000);
90  h1->GetYaxis()->SetRangeUser(-50,200);
91  h1->GetXaxis()->SetTitle("Z(cm)");
92  h1->GetYaxis()->SetTitle("a_{xx}");
93 
94  /* Plot frame histogram */
95  TCanvas *c1 = new TCanvas();
96  h1->Draw("AXIS");
97 
98  leg->Draw();
99 
100  g_0->Draw("Lsame");
101  g_a11->Draw("Psame");
102  g_a21->Draw("Psame");
103  g_a12->Draw("Psame");
104  g_a22->Draw("Psame");
105 
106  c1->Print("amatrix_new.eps");
107 
108  return 0;
109 }