EIC-IR/macros_plots/plot_track_deviations.C

0001 #include "TraceBox.h"
0002 #include <cmath>
0003
0004 TGraph*
0005 graph_delta( TGraph* g1, TGraph* g2 )
0006 {
0007   /* make sure both graph have smae number of entries- return NULL pointer if they don't */
0008   if ( g1->GetN() != g2->GetN() )
0009     {
0010       cerr << "Mismatch in number of points for graphs in calculating delta!" << endl;
0011       return NULL;
0012     }
0013
0014   TGraph* gdelta = new TGraph( g1->GetN() );// (TGraph*)g1->Clone("gdelta");
0015
0016   for ( unsigned i = 0; i < g1->GetN(); i++ )
0017     {
0018       Double_t x1 = 0;
0019       Double_t y1 = 0;
0020       Double_t x2 = 0;
0021       Double_t y2 = 0;
0022
0023       g1->GetPoint( i , x1 , y1 );
0024       g2->GetPoint( i , x2 , y2 );
0025
0026       Double_t xdelta = x1;
0027       Double_t ydelta = y1-y2;
0028
0029       gdelta->SetPoint(i, xdelta, ydelta);
0030     }
0031
0032   return gdelta;
0033 }
0034
0035 int
0036 plot_track_deviations()
0037 {
0038   /**** Chose input files ****/
0039
0040   /* Open iput file with trajectories from GEANT4 */
0041   TFile *file_nominal = new TFile( "data_amatrix/G4EICIR_DSTReader_p250_vx0_vy0_theta22mrad.root", "OPEN" );
0042   TFile *file_nominal_upper = new TFile( "data_amatrix/G4EICIR_DSTReader_p250_vx0_vy1_theta0mrad.root", "OPEN" ); // +200 MeV pT
0043   TFile *file_nominal_lower = new TFile( "data_amatrix/G4EICIR_DSTReader_p250_vx0_vy0_theta10mrad.root", "OPEN" ); // -200 MeV pT
0044
0045   TFile *file_trial = new TFile( "data_amatrix/G4EICIR_DSTReader_p250_vx0_vy0_theta10mrad.root", "OPEN" );
0046
0047   /* Get tree from file */
0048   TTree *tree_nominal = (TTree*)file_nominal->Get("T");
0049   TTree *tree_nominal_upper = (TTree*)file_nominal_upper->Get("T");
0050   TTree *tree_nominal_lower = (TTree*)file_nominal_lower->Get("T");
0051
0052   TTree *tree_trial = (TTree*)file_trial->Get("T");
0053
0054   /* get number of hits */
0055   int nhits = 0;
0056   tree_nominal->SetBranchAddress("n_G4HIT_FWDDISC",&nhits);
0057   tree_nominal->GetEntry(0);
0058   cout << "hits: " << nhits << endl;
0059
0060   /* create graph of nominal beam position (z) */
0061   tree_nominal->Draw( "0.5*(G4HIT_FWDDISC.y[][0]+G4HIT_FWDDISC.y[][1]) : 0.5*(G4HIT_FWDDISC.z[][0]+G4HIT_FWDDISC.z[][1])","Entry$==0","");
0062   TGraph* g_nominal = new TGraph(nhits, &(tree_nominal->GetV2()[0]), &(tree_nominal->GetV1()[0]));
0063   g_nominal->SetName("g_nominal");
0064
0065   /* create graph of beam positions relative to nominal */
0066   TGraph* g_nominal_delta = graph_delta(g_nominal, g_nominal);
0067   if ( ! g_nominal_delta )
0068     return 1;
0069   g_nominal_delta->SetName("g_nominal_delta");
0070   g_nominal_delta->SetLineWidth(2);
0071
0072   /* create graph of nominal_upper beam position (z) */
0073   tree_nominal_upper->Draw( "0.5*(G4HIT_FWDDISC.y[][0]+G4HIT_FWDDISC.y[][1]) : 0.5*(G4HIT_FWDDISC.z[][0]+G4HIT_FWDDISC.z[][1])","Entry$==0","");
0074   TGraph* g_nominal_upper = new TGraph(nhits, &(tree_nominal_upper->GetV2()[0]), &(tree_nominal_upper->GetV1()[0]));
0075   g_nominal_upper->SetName("g_nominal_upper");
0076
0077   /* create graph of beam positions relative to nominal_upper */
0078   TGraph* g_nominal_upper_delta = graph_delta(g_nominal_upper, g_nominal);
0079   if ( ! g_nominal_upper_delta )
0080     return 1;
0081   g_nominal_upper_delta->SetName("g_nominal_upper_delta");
0082   g_nominal_upper_delta->SetLineWidth(2);
0083   g_nominal_upper_delta->SetLineColor(kBlue);
0084
0085   /* create graph of nominal_lower beam position (z) */
0086   tree_nominal_lower->Draw( "0.5*(G4HIT_FWDDISC.y[][0]+G4HIT_FWDDISC.y[][1]) : 0.5*(G4HIT_FWDDISC.z[][0]+G4HIT_FWDDISC.z[][1])","Entry$==0","");
0087   TGraph* g_nominal_lower = new TGraph(nhits, &(tree_nominal_lower->GetV2()[0]), &(tree_nominal_lower->GetV1()[0]));
0088   g_nominal_lower->SetName("g_nominal_lower");
0089
0090   /* create graph of beam positions relative to nominal_lower */
0091   TGraph* g_nominal_lower_delta = graph_delta(g_nominal_lower, g_nominal);
0092   if ( ! g_nominal_lower_delta )
0093     return 1;
0094   g_nominal_lower_delta->SetName("g_nominal_lower_delta");
0095   g_nominal_lower_delta->SetLineWidth(2);
0096   g_nominal_lower_delta->SetLineColor(kRed);
0097
0098   /* Create Legend */
0099 //  TLegend* leg = new TLegend(0.3,0.7,0.6,0.9);
0100 //  leg->AddEntry(g_0,"Nominal beam","l");
0101 //  leg->AddEntry(g_a11,"a11","p");
0102 //  leg->AddEntry(g_a21,"a21","p");
0103 //  leg->AddEntry(g_a12,"a12","p");
0104 //  leg->AddEntry(g_a22,"a22","p");
0105
0106   /* Create frame histogram for plot */
0107   TH1F *h1 = new TH1F("h1","",10,0,15000);
0108   h1->GetXaxis()->SetRangeUser(0,12000);
0109   h1->GetYaxis()->SetRangeUser(-10,10);
0110   h1->GetXaxis()->SetTitle("z (cm)");
0111   h1->GetYaxis()->SetTitle("#Delta_{x} (mm)");
0112
0113   /* Plot frame histogram */
0114   TCanvas *c1 = new TCanvas();
0115   h1->Draw("AXIS");
0116
0117   //leg->Draw();
0118
0119   g_nominal_delta->Draw("Lsame");
0120   g_nominal_upper_delta->Draw("Lsame");
0121   g_nominal_lower_delta->Draw("Lsame");
0122
0123   c1->Print("track_deviations_new.eps");
0124
0125   return 0;
0126 }