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