EIC-IR/macros_plots/plot_amatrix.C

0001 #include "TraceBox.h"
0002 #include <cmath>
0003 int
0004 plot_amatrix()
0005 {
0006   /**** Chose input files ****/
0007
0008   /* Open iput file with trajectories from GEANT4 */
0009   TFile *fin_0 =   new TFile( "data_amatrix/G4EICIR_DSTReader_p250_vx0_vy0_theta22mrad.root", "OPEN" );
0010   TFile *fin_a11 = new TFile( "data_amatrix/G4EICIR_DSTReader_p250_vx0_vy1_theta0mrad.root", "OPEN" );
0011   TFile *fin_a21 = new TFile( "data_amatrix/G4EICIR_DSTReader_p250_vx0_vy0_theta10mrad.root", "OPEN" );
0012   TFile *fin_a12 = new TFile( "data_amatrix/G4EICIR_DSTReader_p250_vx0_vy1_theta0mrad.root", "OPEN" );
0013   TFile *fin_a22 = new TFile( "data_amatrix/G4EICIR_DSTReader_p250_vx0_vy0_theta10mrad.root", "OPEN" );
0014
0015   /* Get tree from file */
0016   TTree *tin_0 = (TTree*)fin_0->Get("T");
0017   TTree *tin_a11 = (TTree*)fin_a11->Get("T");
0018   TTree *tin_a21 = (TTree*)fin_a21->Get("T");
0019   TTree *tin_a12 = (TTree*)fin_a12->Get("T");
0020   TTree *tin_a22 = (TTree*)fin_a22->Get("T");
0021
0022   /* define initial values for read-in trajectories */
0023   float tin_a11_vy = 1.0; // this is y_0
0024   float tin_a21_theta = 10.0; // this is theta_y*
0025   float tin_a12_vy = 1.0; // this is y_0
0026   float tin_a22_theta = 10.0; // theta_y*
0027
0028   /* get number of hits */
0029   int nhits = 0;
0030   tin_a11->SetBranchAddress("n_G4HIT_FWDDISC",&nhits);
0031   tin_a11->GetEntry(0);
0032   cout << "hits: " << nhits << endl;
0033
0034   /* create graph of nominal beam position (z) */
0035   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","");
0036   TGraph* g_0 = new TGraph(nhits, &(tin_0->GetV2()[0]), &(tin_0->GetV1()[0]));
0037
0038   /* create graph of a11(z) */
0039   //   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","");
0040
0041   /* a11 = y/y_0 */
0042   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","" );
0043   TGraph* g_a11 = new TGraph(nhits, &(tin_a11->GetV2()[0]), &(tin_a11->GetV1()[0]));
0044   g_a11->SetMarkerStyle(7);
0045   g_a11->SetMarkerSize(1);
0046   g_a11->SetMarkerColor(kRed);
0047
0048  /* create graph of a21(z) */
0049   //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","");
0050
0051 /* a21 = y/theta* */
0052   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","");
0053   TGraph* g_a21 = new TGraph(nhits, &(tin_a21->GetV2()[0]), &(tin_a21->GetV1()[0]));
0054   g_a21->SetMarkerStyle(7);
0055   g_a21->SetMarkerSize(1);
0056   g_a21->SetMarkerColor(kBlue);
0057
0058   /* create graph of a12(z) */
0059   //  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","");
0060
0061   /* a12 = theta_y/y_0, theta*=atan( |exit-entry| / thickness ), here thickness = 1  */
0062   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","");
0063   TGraph* g_a12 = new TGraph(nhits, &(tin_a12->GetV2()[0]), &(tin_a12->GetV1()[0]));
0064   g_a12->SetMarkerStyle(7);
0065   g_a12->SetMarkerSize(1);
0066   g_a12->SetMarkerColor(kGreen);
0067
0068  /* create graph of a22(z) */
0069   //  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","");
0070
0071   /* a22 = thetay/thetay*, theta*=atan( |exit-entry| / thickness ), here thickness = 1 */
0072   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","");
0073   TGraph* g_a22 = new TGraph(nhits, &(tin_a22->GetV2()[0]), &(tin_a22->GetV1()[0]));
0074   g_a22->SetMarkerStyle(7);
0075   g_a22->SetMarkerSize(1);
0076   g_a22->SetMarkerColor(kBlue);
0077   g_a22->SetMarkerColor(kOrange);
0078
0079   /* Create Legend */
0080   TLegend* leg = new TLegend(0.3,0.7,0.6,0.9);
0081   leg->AddEntry(g_0,"Nominal beam","l");
0082   leg->AddEntry(g_a11,"a11","p");
0083   leg->AddEntry(g_a21,"a21","p");
0084   leg->AddEntry(g_a12,"a12","p");
0085   leg->AddEntry(g_a22,"a22","p");
0086
0087   /* Create frame histogram for plot */
0088   TH1F *h1 = new TH1F("h1","",10,0,15000);
0089   h1->GetXaxis()->SetRangeUser(0,12000);
0090   h1->GetYaxis()->SetRangeUser(-50,200);
0091   h1->GetXaxis()->SetTitle("Z(cm)");
0092   h1->GetYaxis()->SetTitle("a_{xx}");
0093
0094   /* Plot frame histogram */
0095   TCanvas *c1 = new TCanvas();
0096   h1->Draw("AXIS");
0097
0098   leg->Draw();
0099
0100    g_0->Draw("Lsame");
0101    g_a11->Draw("Psame");
0102    g_a21->Draw("Psame");
0103    g_a12->Draw("Psame");
0104    g_a22->Draw("Psame");
0105
0106   c1->Print("amatrix_new.eps");
0107
0108   return 0;
0109 }