macros/diskinematics_analysis/makePlot_truth_xQ2_Studies.C

0001 /**
0002
0003  *\Brief Macro to compare EIC theoretical kinematics to Pythia simulation data
0004  *\Author: Thomas Krahulik <thomas.krahulik@stonybrook.edu>
0005
0006  **/
0007
0008
0009 int makePlot_truth_xQ2_Studies(
0010                    const TString infile = "../../data/Sample_DISReco_ep.root",
0011                    //const TString infile = "../../data/EventGenAna_Pythia6_DIS_10x250_100k.root",
0012                    const bool save_figures = true
0013                    )
0014
0015 {
0016   gStyle->SetOptStat(kFALSE);
0017
0018   /*--------------Get Input File--------------*/
0019
0020   TFile *f_truth = new TFile( infile, "OPEN");
0021   TTree* t_truth = (TTree*)f_truth->Get("event_truth");
0022   t_truth->Print();
0023
0024   /*------------------------------------------*/
0025   /*---------------Define Cuts----------------*/
0026
0027   TCut electron_cut = "em_evtgen_pid == 11";
0028
0029   /*------------------------------------------*/
0030   /*---------x-Q2 Space for Electrons---------*/
0031   /*------------------------------------------*/
0032
0033   /*------------------------------------------*/
0034   /*------------Create Histogram--------------*/
0035   /*------------------------------------------*/
0036   /**
0037    *\ The creation of this histogram includes a method for setting the bins to a logarithmic scale along with the axes
0038    **/
0039   /*------------------------------------------*/
0040
0041   TH2F* hQ2x_e = new TH2F("hQ2x_e", "", 40, -4, 0, 60, 0, 3);
0042   TAxis *xaxis = hQ2x_e->GetXaxis();
0043   TAxis *yaxis = hQ2x_e->GetYaxis();
0044   int xbins = xaxis->GetNbins();
0045   int ybins = yaxis->GetNbins();
0046
0047   Axis_t xmin = xaxis->GetXmin();
0048   Axis_t xmax = xaxis->GetXmax();
0049   Axis_t xwidth = (xmax - xmin ) / xbins;
0050   Axis_t *new_xbins = new Axis_t[xbins + 1];
0051
0052   for( int i =0; i <= xbins; i++)
0053     {
0054       new_xbins[i] = TMath::Power( 10, xmin + i*xwidth);
0055     }
0056   xaxis->Set(xbins, new_xbins);
0057
0058   Axis_t ymin = yaxis->GetXmin();
0059   Axis_t ymax = yaxis->GetXmax();
0060   Axis_t ywidth = (ymax - ymin) / ybins;
0061   Axis_t *new_ybins = new Axis_t[ybins + 1];
0062
0063   for( int i2 =0; i2 <= ybins; i2++)
0064     {
0065       new_ybins[i2] = TMath::Power( 10, ymin + i2*ywidth);
0066     }
0067
0068   yaxis->Set(ybins, new_ybins);
0069   /*------------------------------------------*/
0070   /*-------------Plot x-Q2 Space--------------*/
0071   /*------------------------------------------*/
0072
0073   /*------------------------------------------*/
0074   /*----------Scattered Lepton Angle----------*/
0075   /*------------------------------------------*/
0076
0077   TCanvas *cQ2x_e1 = new TCanvas("cQ2x_e1");
0078   cQ2x_e1->SetLogx();
0079   cQ2x_e1->SetLogy();
0080   t_truth->Draw("evtgen_Q2:evtgen_x>>hQ2x_e", "" , "colz");
0081
0082   hQ2x_e->GetXaxis()->SetTitle("x");
0083   hQ2x_e->GetYaxis()->SetTitle("Q^{2} [GeV^{2}]");
0084
0085   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);
0086   f_LepA_0->SetParameter( 0, 10);
0087   f_LepA_0->SetParameter( 1, 250);
0088   f_LepA_0->SetParameter( 2, (0 * TMath::Pi())/180 );
0089   TF1 *f_LepA_90 = (TF1*)f_LepA_0->Clone("f_LepA_90");
0090   f_LepA_90->SetParameter( 2, (90 * TMath::Pi())/180 );
0091   TF1 *f_LepA_150 = (TF1*)f_LepA_0->Clone("f_LepA_150");
0092   f_LepA_150->SetParameter( 2, (150 * TMath::Pi())/180 );
0093   TF1 *f_LepA_170 = (TF1*)f_LepA_0->Clone("f_LepA_170");
0094   f_LepA_170->SetParameter( 2, (170 * TMath::Pi())/180 );
0095
0096   f_LepA_0->Draw("SAME");
0097   f_LepA_0->SetLineColor(1);
0098   f_LepA_90->Draw("SAME");
0099   f_LepA_90->SetLineColor(1);
0100   f_LepA_150->Draw("SAME");
0101   f_LepA_150->SetLineColor(1);
0102   f_LepA_170->Draw("SAME");
0103   f_LepA_170->SetLineColor(1);
0104
0105   /*---------------Label Lines----------------*/
0106
0107   TPaveText *pt_LepA_0 = new TPaveText(0.4, 1.2, 0.6, 1.4);
0108   pt_LepA_0->AddText("0^{/circ}");
0109   pt_LepA_0->Draw();
0110
0111   /*------------------------------------------*/
0112   /*---------Scattered Lepton Energy----------*/
0113   /*------------------------------------------*/
0114
0115   TCanvas *cQ2x_e2 = new TCanvas("cQ2x_e2");
0116   cQ2x_e2->SetLogx();
0117   cQ2x_e2->SetLogy();
0118   t_truth->Draw("evtgen_Q2:evtgen_x>>hQ2x_e", "" , "colz");
0119
0120   hQ2x_e->GetXaxis()->SetTitle("x");
0121   hQ2x_e->GetYaxis()->SetTitle("Q^{2} [GeV^{2}]");
0122
0123   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);
0124   f_LepE_5->SetParameter( 0, 10);
0125   f_LepE_5->SetParameter( 1, 250);
0126   f_LepE_5->SetParameter( 2, 5);
0127   TF1 *f_LepE_9 = (TF1*)f_LepE_5->Clone("f_LepE_9");
0128   f_LepE_9->SetParameter( 2, 9);
0129   TF1 *f_LepE_11 = (TF1*)f_LepE_5->Clone("f_LepE_11");
0130   f_LepE_11->SetParameter( 2, 11);
0131   TF1 *f_LepE_15 = (TF1*)f_LepE_5->Clone("f_LepE_15");
0132   f_LepE_15->SetParameter( 2, 15);
0133
0134   f_LepE_5->Draw("SAME");
0135   f_LepE_5->SetLineColor(1);
0136   f_LepE_9->Draw("SAME");
0137   f_LepE_9->SetLineColor(1);
0138   f_LepE_11->Draw("SAME");
0139   f_LepE_11->SetLineColor(1);
0140   f_LepE_15->Draw("SAME");
0141   f_LepE_15->SetLineColor(1);
0142
0143   /*------------------------------------------*/
0144   /*------------Current Jet Angle-------------*/
0145   /*------------------------------------------*/
0146
0147   TCanvas *cQ2x_e3 = new TCanvas("cQ2x_e3");
0148   cQ2x_e3->SetLogx();
0149   cQ2x_e3->SetLogy();
0150   t_truth->Draw("evtgen_Q2:evtgen_x>>hQ2x_e", "" , "colz");
0151
0152   hQ2x_e->GetXaxis()->SetTitle("x");
0153   hQ2x_e->GetYaxis()->SetTitle("Q^{2} [GeV^{2}]");
0154
0155   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);
0156   f_JetA_180->SetParameter( 0 , 10);
0157   f_JetA_180->SetParameter( 1 , 250.);
0158   f_JetA_180->SetParameter( 2 , (180 * TMath::Pi())/180 );
0159   TF1 *f_JetA_150 = (TF1*)f_JetA_180->Clone("f_JetA_150");
0160   f_JetA_150->SetParameter( 2 , (150 * TMath::Pi())/180 );
0161   TF1 *f_JetA_30 = (TF1*)f_JetA_180->Clone("f_JetA_30");
0162   f_JetA_30->SetParameter( 2 , (30 * TMath::Pi())/180 );
0163   TF1 *f_JetA_10 = (TF1*)f_JetA_180->Clone("f_JetA_10");
0164   f_JetA_10->SetParameter( 2 , (10 * TMath::Pi())/180 );
0165
0166   f_JetA_180->Draw("SAME");
0167   f_JetA_180->SetLineColor(1);
0168   f_JetA_150->Draw("SAME");
0169   f_JetA_150->SetLineColor(1);
0170   f_JetA_30->Draw("SAME");
0171   f_JetA_30->SetLineColor(1);
0172   f_JetA_10->Draw("SAME");
0173   f_JetA_10->SetLineColor(1);
0174
0175   /*------------------------------------------*/
0176   /*------------Current Jet Energy------------*/
0177   /*------------------------------------------*/
0178
0179   TCanvas *cQ2x_e4 = new TCanvas("cQ2x_e4");
0180   cQ2x_e4->SetLogx();
0181   cQ2x_e4->SetLogy();
0182   t_truth->Draw("evtgen_Q2:evtgen_x>>hQ2x_e", "" , "colz");
0183
0184   hQ2x_e->GetXaxis()->SetTitle("x");
0185   hQ2x_e->GetYaxis()->SetTitle("Q^{2} [GeV^{2}]");
0186
0187   TF1 *f_JetE_5 = new TF1("f_JetE_5" , "(4*x*[0]*[1]*([2]-x*[1])) / ([0] - x*[1])", 10e-6, 1);
0188   f_JetE_5->SetParameter( 0 , 10);
0189   f_JetE_5->SetParameter( 1 , 250.);
0190   f_JetE_5->SetParameter( 2 , 5. );
0191   TF1 *f_JetE_7 = (TF1*)f_JetE_5->Clone("f_JetE_7");
0192   f_JetE_7->SetParameter( 2 , 7);
0193   TF1 *f_JetE_10 = (TF1*)f_JetE_5->Clone("f_JetE_10");
0194   f_JetE_10->SetParameter( 2 , 10);
0195   TF1 *f_JetE_20 = (TF1*)f_JetE_5->Clone("f_JetE_20");
0196   f_JetE_20->SetParameter( 2 , 20);
0197
0198   f_JetE_5->Draw("SAME");
0199   f_JetE_5->SetLineColor(1);
0200   f_JetE_7->Draw("SAME");
0201   f_JetE_7->SetLineColor(1);
0202   f_JetE_10->Draw("SAME");
0203   f_JetE_10->SetLineColor(1);
0204   f_JetE_20->Draw("SAME");
0205   f_JetE_20->SetLineColor(1);
0206
0207   /*------------------------------------------*/
0208   /*---------------Inelasticity---------------*/
0209   /*------------------------------------------*/
0210
0211   TCanvas *cQ2x_e5 = new TCanvas("cQ2x_e5");
0212   cQ2x_e5->SetLogx();
0213   cQ2x_e5->SetLogy();
0214   t_truth->Draw("evtgen_Q2:evtgen_x>>hQ2x_e", "" , "colz");
0215
0216   hQ2x_e->GetXaxis()->SetTitle("x");
0217   hQ2x_e->GetYaxis()->SetTitle("Q^{2} [GeV^{2}]");
0218
0219   TF1 *f_y1 = new TF1("f_y1", "4*x*[0]*[1]*[2]", 10e-5, 1);
0220   f_y1->SetParameter( 0, 10);
0221   f_y1->SetParameter( 1, 250);
0222   f_y1->SetParameter( 2, 1);
0223   TF1 *f_y01 = (TF1*)f_y1->Clone("f_y01");
0224   f_y01->SetParameter(2 , 0.1);
0225   TF1 *f_y001 = (TF1*)f_y1->Clone("f_y001");
0226   f_y001->SetParameter(2 , 0.01);
0227   TF1 *f_y0001 = (TF1*)f_y1->Clone("f_y0001");
0228   f_y0001->SetParameter(2 , 0.001);
0229
0230   f_y1->Draw("SAME");
0231   f_y1->SetLineColor(1);
0232   f_y01->Draw("SAME");
0233   f_y01->SetLineColor(1);
0234   f_y001->Draw("SAME");
0235   f_y001->SetLineColor(1);
0236   f_y0001->Draw("SAME");
0237   f_y0001->SetLineColor(1);
0238
0239   /*------------------------------------------*/
0240   /*-----------Final Hadronic Mass------------*/
0241   /*------------------------------------------*/
0242
0243   TCanvas *cQ2x_e6 = new TCanvas("cQ2x_e6");
0244   cQ2x_e6->SetLogx();
0245   cQ2x_e6->SetLogy();
0246   t_truth->Draw("evtgen_Q2:evtgen_x>>hQ2x_e", "" , "colz");
0247
0248   hQ2x_e->GetXaxis()->SetTitle("x");
0249   hQ2x_e->GetYaxis()->SetTitle("Q^{2} [GeV^{2}]");
0250
0251   TF1 *f_W10 = new TF1("f_W10" , "([1] - [0]**2) / ((1/x) - 1)", 10e-7, 1);
0252   f_W10->SetParameter( 0 , 0.938);
0253   f_W10->SetParameter( 1 , (10)**2 );
0254   TF1 *f_W20 = (TF1*)f_W10->Clone("f_W20");
0255   f_W20->SetParameter( 1 , (20)**2 );
0256   TF1 *f_W50 = (TF1*)f_W10->Clone("f_W50");
0257   f_W50->SetParameter( 1 , (50)**2 );
0258   TF1 *f_W100 = (TF1*)f_W10->Clone("f_W100");
0259   f_W100->SetParameter( 1 , (100)**2 );
0260   TF1 *f_W3 = (TF1*)f_W10->Clone("f_W3");
0261   f_W3->SetParameter( 1 , (3)**2 );
0262
0263   f_W10->Draw("SAME");
0264   f_W10->SetLineColor(1);
0265   f_W20->Draw("SAME");
0266   f_W20->SetLineColor(1);
0267   f_W50->Draw("SAME");
0268   f_W50->SetLineColor(1);
0269   f_W100->Draw("SAME");
0270   f_W100->SetLineColor(1);
0271   f_W3->Draw("SAME");
0272   f_W3->SetLineColor(1);
0273
0274   if ( save_figures )
0275     {
0276       cQ2x_e1->Print("plots/Pythia_Q2x_LepA_e.eps");
0277       cQ2x_e2->Print("plots/Pythia_Q2x_LepE_e.eps");
0278       cQ2x_e3->Print("plots/Pythia_Q2x_JetA_e.eps");
0279       cQ2x_e4->Print("plots/Pythia_Q2x_JetE_e.eps");
0280       cQ2x_e5->Print("plots/Pythia_Q2x_y_e.eps");
0281       cQ2x_e6->Print("plots/Pythia_Q2x_W_e.eps");
0282     }
0283
0284   return 0;
0285
0286 }