EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
analysis.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file analysis.C
1 
2 //
3 // This script is specifically tuned to be run on the PC farm
4 // production output; may want to optimize all the logic for speed later; for
5 // now just loop through directories one by one;
6 //
7 // If run scripts by hand (and produce single reconstruction.root file), then just:
8 //
9 // root -l 'analysis.C(".", 1)'
10 //
11 
12 // For production mode running; otherwise expect to find files right in
13 // the "dir" directory (see parameters);
14 //#define _EXPECT_SUBDIRECTORIES_
15 
16 //#define _SHOW_PURITY_
17 //#define _SHOW_KINEMATICS_
18 
19 #define _X_MIN_ (1E-5)
20 #define _X_MAX_ (1)
21 #if defined(_SHOW_PURITY_) || defined(_SHOW_KINEMATICS_)
22 #define _X_BNUM_ (5*5*1)
23 #else
24 #define _X_BNUM_ (200)
25 #endif
26 #define _LOGX_BWID_ ((log(_X_MAX_)-log(_X_MIN_))/_X_BNUM_)
27 
28 #define _Q2_MIN_ (1)
29 #define _Q2_MAX_ (1E3)
30 #if defined(_SHOW_PURITY_) || defined(_SHOW_KINEMATICS_)
31 #define _Q2_BNUM_ (3*4)
32 #else
33 #define _Q2_BNUM_ (200)
34 #endif
35 #define _LOGQ2_BWID_ ((log(_Q2_MAX_)-log(_Q2_MIN_))/_Q2_BNUM_)
36 
37 #define _Y_BNUM_ (200)
38 
39 void analysis(const char *dir, unsigned nn = 1)
40 {
41 #if defined(_SHOW_PURITY_) && defined(_SHOW_KINEMATICS_)
42  printf("defined(_SHOW_PURITY_) && defined(_SHOW_KINEMATICS_): one at a time, please!\n");
43  exit(0);
44 #endif
45 
46  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
47  gStyle->SetOptStat(0);
48  gStyle->SetPalette(55, 0);
49  gStyle->SetLabelSize(0.04, "xy");
50  gStyle->SetPadLeftMargin(0.12);
51  gStyle->SetPadBottomMargin(0.13);
52  gStyle->SetTextFont(52);
53 
54  unsigned ok [_X_BNUM_][_Q2_BNUM_];
55  unsigned in [_X_BNUM_][_Q2_BNUM_];
56  memset(ok, 0x00, sizeof(ok));
57  memset(in, 0x00, sizeof(in));
58 
59  double xb[_X_BNUM_+1], q2b[_Q2_BNUM_+1];
60 
61  for(int ip=0; ip<=_X_BNUM_; ip++)
62  {
63  xb [ip] = exp(log(_X_MIN_) + ip*_LOGX_BWID_);
64  //printf("%2d -> %f\n", ip, xb[ip]);
65  }
66  for(int iq=0; iq<=_Q2_BNUM_; iq++)
67  {
68  q2b[iq] = exp(log(_Q2_MIN_) + iq*_LOGQ2_BWID_);
69  //printf("%2d -> %f\n", iq, q2b[iq]);
70  }
71 
72  TH2D *xx = new TH2D("xx", "", _X_BNUM_, xb, _X_BNUM_, xb);
73  TH2D *yy = new TH2D("yy", "", _Y_BNUM_, 0.01, 1.0, _Y_BNUM_, 0.01, 1.0);
74  TH1D *y1d = new TH1D("y1d", "", _Y_BNUM_, 0.01, 1.0);
75  TH2D *qq = new TH2D("qq", "", _Q2_BNUM_, q2b, _Q2_BNUM_, q2b);
76  TH1D *q1d = new TH1D("q1d", "q1d", _Q2_BNUM_, q2b);
77  TH2D *xq = new TH2D("xq", "xq", _X_BNUM_, xb, _Q2_BNUM_, q2b);
78  TH2D *mg = new TH2D("mg", "", _X_BNUM_, xb, _Q2_BNUM_, q2b);
79 
80 #ifndef _EXPECT_SUBDIRECTORIES_
81  nn = 1;
82 #endif
83 
84  // Input simulated & reconstructed files;
85  for(unsigned iqn=0; iqn<nn; iqn++) {
86  char qdir[128] = "/";
87 #ifdef _EXPECT_SUBDIRECTORIES_
88  snprintf(qdir, 128-1, "/%05d/", iqn);
89 #endif
90  TString mcInFile = TString(dir) + TString(qdir) + "simulation.root";
91  TString rcInFile = TString(dir) + TString(qdir) + "reconstruction.root";
92 
93  EicRootManager *io = new EicRootManager(mcInFile, 0, rcInFile);
94  if (io->GetEicRcEventBranch()) {
95  printf("chunk# %4d, %6d entries ...\n", iqn, io->GetEicRcTreeEntries());
96 
97  for(unsigned ev=0; ev<io->GetEicRcTreeEntries(); ev++) {
98  io->GetEicRcTreeEntry(ev);
99 
100  const EicRcEvent *rcEvent = io->GetEicRcEvent();
101  erhic::EventMC *mcEvent = rcEvent->GetGenMcEvent();
102 
103  if (!rcEvent->GetNTracks()) continue;
104 
105  {
106  erhic::LeptonKinematicsComputer rcDis(*rcEvent);
107  erhic::DisKinematics *rcKin = rcDis.Calculate();
108  //printf("X: %7.5f; Q^2: %8.1f\n", rcKin->mX, rcKin->mQ2);
109  unsigned rcIx = mg->GetXaxis()->FindBin(rcKin->mX) - 1;
110  unsigned rcIq = mg->GetYaxis()->FindBin(rcKin->mQ2) - 1;
111 
112  erhic::LeptonKinematicsComputer mcDis(*mcEvent);
113  erhic::DisKinematics *mcKin = mcDis.Calculate();
114  //printf("X: %7.5f; Q^2: %8.1f\n\n", mcKin->mX, mcKin->mQ2);
115  unsigned mcIx = mg->GetXaxis()->FindBin(mcKin->mX) - 1;
116  unsigned mcIq = mg->GetYaxis()->FindBin(mcKin->mQ2) - 1;
117 
118  if (mcIx >=0 && mcIx < _X_BNUM_ &&
119  mcIq >=0 && mcIq < _Q2_BNUM_) {
120  if (mcIx == rcIx && mcIq == rcIq)
121  ok[mcIx][mcIq]++;
122  else {
123  if (rcIx >=0 && rcIx < _X_BNUM_ &&
124  rcIq >=0 && rcIq < _Q2_BNUM_)
125  in[rcIx][rcIq]++;
126  } //if
127 
128  // Yes, in particular want Y-smearing plot to have Q^2 < 1 cut applied;
129  xx->Fill (mcKin->mX, rcKin->mX);
130  yy->Fill (mcKin->mY, rcKin->mY);
131  y1d->Fill(mcKin->mY);
132  qq->Fill (mcKin->mQ2, rcKin->mQ2);
133 
134  q1d->Fill(mcKin->mQ2);
135  xq->Fill (mcKin->mX, mcKin->mQ2);
136  } //if
137  }
138  } //for ev
139  }
140  else
141  printf("chunk# %4d: corrupted (or missing)\n", iqn);
142 
143  delete io;
144  } //for iqn
145 
146  for(int ix=0; ix<_X_BNUM_; ix++)
147  for(int iq=0; iq<_Q2_BNUM_; iq++) {
148  double value = ok[ix][iq] ? 1.0*ok[ix][iq]/(ok[ix][iq]+in[ix][iq]) : 0.0;
149  mg->SetBinContent(ix+1, iq+1, value);
150  } //for ix..iq
151 
152 #if !defined(_SHOW_PURITY_) && !defined(_SHOW_KINEMATICS_)
153  TCanvas *c11 = new TCanvas("c11", "c11", 0, 0, 1400, 400);
154  c11->Divide(3,1);
155 
156  c11->cd(1);
157  gPad->SetLogz();
158 
159  yy->SetMinimum(2);
160  yy->GetXaxis()->SetTitle("Simulated Y");
161  yy->GetXaxis()->SetTitleSize(0.05);
162  yy->GetXaxis()->SetTitleFont(52);
163  yy->GetXaxis()->SetLabelSize(0.04);
164  yy->GetXaxis()->SetLabelFont(52);
165  yy->GetXaxis()->SetTitleOffset(1.2);
166  yy->GetYaxis()->SetTitle("Reconstructed Y");
167  yy->GetYaxis()->SetTitleSize(0.05);
168  yy->GetYaxis()->SetTitleFont(52);
169  yy->GetYaxis()->SetLabelSize(0.04);
170  yy->GetYaxis()->SetLabelFont(52);
171  yy->Draw("COLZ");
172 
173  c11->cd(2);
174  gPad->SetLogx();
175  gPad->SetLogy();
176  gPad->SetLogz();
177  xx->SetMinimum(2);
178  xx->GetXaxis()->SetTitle("Simulated X_{Bj}");
179  xx->GetXaxis()->SetTitleSize(0.05);
180  xx->GetXaxis()->SetTitleFont(52);
181  xx->GetXaxis()->SetLabelSize(0.04);
182  xx->GetXaxis()->SetLabelFont(52);
183  xx->GetXaxis()->SetTitleOffset(1.2);
184  xx->GetYaxis()->SetTitle("Reconstructed X_{Bj}");
185  xx->GetYaxis()->SetTitleSize(0.05);
186  xx->GetYaxis()->SetTitleFont(52);
187  xx->GetYaxis()->SetLabelSize(0.04);
188  xx->GetYaxis()->SetLabelFont(52);
189  xx->Draw("COLZ");
190 
191  c11->cd(3);
192  gPad->SetLogx();
193  gPad->SetLogy();
194  gPad->SetLogz();
195  qq->SetMinimum(2);
196  qq->GetXaxis()->SetTitle("Simulated Q^{2}, [GeV^{2}]");
197  qq->GetXaxis()->SetTitleSize(0.05);
198  qq->GetXaxis()->SetTitleFont(52);
199  qq->GetXaxis()->SetLabelSize(0.04);
200  qq->GetXaxis()->SetLabelFont(52);
201  qq->GetXaxis()->SetTitleOffset(1.2);
202  qq->GetYaxis()->SetTitle("Reconstructed Q^{2}, [GeV^{2}]");
203  qq->GetYaxis()->SetTitleSize(0.05);
204  qq->GetYaxis()->SetTitleFont(52);
205  qq->GetYaxis()->SetLabelSize(0.04);
206  qq->GetYaxis()->SetLabelFont(52);
207  qq->Draw("COLZ");
208 #else
209  TCanvas *c11 = new TCanvas("c11", "c11", 0, 0, 700, 450);
210 #ifdef _SHOW_PURITY_
211  TH2D *hh = mg;
212 #else
213  TH2D *hh = xq;
214  c11->SetLogz();
215 #endif
216  c11->SetLogx();
217  c11->SetLogy();
218 
219  hh->GetXaxis()->SetTitle("Simulated X_{Bj}");
220  hh->GetXaxis()->SetTitleSize(0.05);
221  hh->GetXaxis()->SetTitleFont(52);
222  hh->GetXaxis()->SetLabelSize(0.04);
223  hh->GetXaxis()->SetLabelFont(52);
224  hh->GetXaxis()->SetTitleOffset(1.2);
225  hh->GetYaxis()->SetTitle("Simulated Q^{2}, [GeV^{2}]");
226  hh->GetYaxis()->SetTitleSize(0.05);
227  hh->GetYaxis()->SetTitleFont(52);
228  hh->GetYaxis()->SetLabelSize(0.04);
229  hh->GetYaxis()->SetLabelFont(52);
230  hh->SetMaximum(1.0);
231  hh->Draw("COLZ");
232  //q1d->Draw();
233 
234  TF1 *y01 = new TF1("y01", "0.01*140*140*x");
235  y01->SetLineColor(kBlack);
236  y01->Draw("SAME");
237  TText *t01 = new TText(0.025, 2.5, "y = 0.01");
238  t01->SetTextAngle(45);
239  t01->Draw();
240 
241  TF1 *y10 = new TF1("y10", "0.10*140*140*x");
242  y10->SetLineColor(kBlack);
243  y10->Draw("SAME");
244  TText *t10 = new TText(0.0025, 2.5, "y = 0.10");
245  t10->SetTextAngle(45);
246  t10->Draw();
247 
248  TF1 *y95 = new TF1("y95", "0.95*140*140*x");
249  y95->SetLineColor(kBlack);
250  y95->Draw("SAME");
251  TText *t95 = new TText(0.00025, 2.5, "y = 0.95");
252  t95->SetTextAngle(45);
253  t95->Draw();
254 #endif
255 } // analysis()