EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
acceptance.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file acceptance.C
1 
2 #define _THETA_PLOT_
3 
4 // Assume 10 GeV energy; 2-d column in pCDR table 3.3;
5 
6 // Beam emittance; converted to [mm];
7 //#define _EX_ (20.0E-6)
8 //#define _EY_ ( 1.0E-6)
9 #define _EX_ (22.0E-6)
10 #define _EY_ ( 3.3E-6)
11 
12 // Beam energy spread;
13 #define _DPP_ (10E-4)
14 
15 // Central trajectory offset; must be tuned by hand;
16 //#define _X0_ ( 0.0)
17 #define _X0_ ( 207.24)
18 #define _Y0_ ( 0.00)
19 
20 // This cut will limit acceptance at small theta;
21 #define _SIGMA_CUT_ (10.0)
22 
23 // Dispersion at LQST location in [mm/10^-3];
24 #define _DX_ (-0.22)
25 
26 // Beta functions around Z ~ -30m for this beam tune; converted to [mm];
27 #define _BETA_X_ (87E3)
28 #define _BETA_Y_ (56E3)
29 
30 // FIXME: hardcode to 10 GeV/c;
31 #define _P0_ (18.0)
32 // Approximately :-);
33 #define _ME_ (0.5E-4)
34 
35 //void acceptance()
36 {
37  // Load basic libraries;
38  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
39 
40  gStyle->SetOptStat(0);
41  gStyle->SetPadBottomMargin(0.13);
42 
43  // Input simulated & reconstructed files;
44  TFile *ff = new TFile("simulation.root");
45  TTree *cbmsim = ff->Get("cbmsim");
46  cbmsim->AddFriend("cbmsim", "digitization.root");
47 
48  TClonesArray *lqstHitArray = new TClonesArray("EicTrackingDigiHitOrth2D");
49  cbmsim->SetBranchAddress("LqstTrackingDigiHit", &lqstHitArray);
50  TClonesArray *mcTrackArray = new TClonesArray("PndMCTrack");
51  cbmsim->SetBranchAddress("MCTrack", &mcTrackArray);
52 
53  TCanvas *c1 = new TCanvas("c1", "c1", 0, 0, 800, 500);
54 
55  TH1D *th0 = new TH1D("th0", "", 250, 0., 25.);
56  th0->SetMinimum(0); th0->SetLineWidth(2); th0->SetLineColor(kBlack);
57  TH1D *th1 = new TH1D("th1", "", 250, 0., 25.); th1->SetLineWidth(2); th1->SetLineColor(kGreen);
58  TH1D *th2 = new TH1D("th2", "", 250, 0., 25.); th2->SetLineWidth(2); th2->SetLineColor(kRed);
59  TH1D *th3 = new TH1D("th3", "", 250, 0., 25.); th3->SetLineWidth(2); th3->SetLineColor(kBlue);
60 
61 
62 #if _TODAY_
63  TH1D *thq3 = new TH1D("thq3", "", 40, 0., 4.);
64  thq3->SetMinimum(0); thq3->SetLineWidth(2); thq3->SetLineColor(kBlue);
65  TH1D *pq3 = new TH1D("pq3", "" , 80, 60., 140.0);
66  pq3->SetLineColor(kBlue); pq3->SetMinimum(0); pq3->SetLineWidth(2);
67  TH2D *pthq3 = new TH2D("pthq3", "" , 80, 60., 140.0, 40, 0., 4.);
68 
69 
70  //TH1D *thres = new TH1D("thres", "", 100, -1., 1.);
71  //TH1D *thyres = new TH1D("thyres", "", 100, -1., 1.);
72  TH1D *pres = new TH1D("pres", "", 100, -10., 10.);
73  TH1D *dpt = new TH1D("dpt", "dpt", 100, -300., 300.);
74 #endif
75  //TH1D *hxl = new TH1D("hxl", "" , 200, 206., 208.0);
76  TH1D *hxl = new TH1D("hxl", "" , 200, -1.0, 1.0);
77  //TH1D *q2h0 = new TH1D("q2h0", "" , 1000, 0.0, 0.03);
78  //TH1D *q2h3 = new TH1D("q2h3", "" , 1000, 0.0, 0.03);
79  TH1D *q2h0 = new TH1D("q2h0", "" , 100, -20.0, 5.0);
80  q2h0->SetMinimum(0); q2h0->SetLineWidth(3); q2h0->SetLineColor(kBlack);
81  TH1D *q2h1 = new TH1D("q2h1", "" , 100, -20.0, 5.0);
82  q2h1->SetLineWidth(2); q2h1->SetLineColor(kGreen);
83  TH1D *q2h3 = new TH1D("q2h3", "" , 100, -20.0, 5.0);
84  q2h3->SetLineWidth(1); q2h3->SetLineColor(kBlue);
85 #if _TODAY_
86  TH1D *pt0 = new TH1D("pt0", "" ,150, 0., 1.5);
87  pt0->SetLineColor(kBlack); pt0->SetMinimum(0); pt0->SetLineWidth(1);
88  TH1D *pt1 = new TH1D("pt1", "", 150, 0., 1.5); pt1->SetLineWidth(3); pt1->SetLineColor(kGreen);
89  TH1D *pt2 = new TH1D("pt2", "", 150, 0., 1.5); pt2->SetLineWidth(3); pt2->SetLineColor(kRed);
90  TH1D *pt3 = new TH1D("pt3", "", 150, 0., 1.5); pt3->SetLineWidth(3); pt3->SetLineColor(kBlue);
91 
92  TH1D *hq2d = new TH1D("hq2d", "hq2d", 100, 0., 100.);
93  TH1D *phi= new TH1D("phi","phi",100, -3.14159, 3.14159);
94  TH2D *acc= new TH2D("acc","acc",100, -10., 10., 100, -10., 10.);
95  //TH1D *x0 = new TH1D("x0", "x0", 1000, -10., 10.);
96  TH1D *x1 = new TH1D("x1", "x1", 3000, -30., 30.);
97  TH1D *sx = new TH1D("sx", "sx", 1000, -10., 10.);
98  //TH1D *xx = new TH1D("xx", "xx", 100, -650., -450.);
99  TH1D *y1 = new TH1D("y1", "y1", 2000, -20., 20.);
100  TH1D *sy = new TH1D("sy", "sy", 1000, -10., 10.);
101  TH2D *xy = new TH2D("xy", "xy", 50, -120., 120., 50, -120., 120.);
102  //TH2D *xy = new TH2D("xy", "xy", 50, -650., -450., 50, -100., 100.);
103 
104  //TGeoRotation *rw = new TGeoRotation();
105  //double _22mrad_ = 0.022 * TMath::RadToDeg();
106  //rw->RotateY(_22mrad_);
107 
108  // FIXME: these are the values tuned for the particular optics of course;
109  double aX[2][2] = {{-1.810, 24.070}, {-0.165, 0.240}}, aXinv[2][2], aYinv = 1/7.725;
110  double det = aX[0][0]*aX[1][1] - aX[0][1]*aX[1][0]; assert(det);
111  aXinv[0][0] = aX[1][1]/det;
112  aXinv[0][1] = -aX[0][1]/det;
113  aXinv[1][0] = -aX[1][0]/det;
114  aXinv[1][1] = aX[0][0]/det;
115 #endif
116 
117  // Beam 1 sigma size at the LQST location;
118  double sigx = sqrt(_EX_*_BETA_X_), sigy = sqrt(_EY_*_BETA_Y_); //printf("%f %f\n", sigx, sigy); exit(0);
119  double sigp = _DX_ * _DPP_ * 1E3; //printf("%f\n", sigp); exit(0);
120  sigx = sqrt(sigx*sigx + sigp*sigp);
121  //printf("%f %f\n", sqrt(_EX_*_BETA_X_), sqrt(_EY_*_BETA_Y_)); exit(0);
122  //printf("%f %f\n", _EX_*_BETA_X_, _EY_*_BETA_Y_); exit(0);
123 
124  // Loop through all events; NB: for box-generated events without secondaries
125  // could simply use cbmsim->Project() as well; in general EicEventAssembler
126  // should be used for "true" physics events;
127  int nEvents = cbmsim->GetEntries();
128  for(unsigned ev=0; ev<nEvents; ev++) {
129  cbmsim->GetEntry(ev);
130 
131  assert(mcTrackArray->GetEntriesFast() == 1);
132  PndMCTrack *mctrack = (PndMCTrack *)mcTrackArray->At(0);
133 
134  if (mctrack->GetPdgCode() == 11 && mctrack->GetMotherID() == -1) {
135  TVector3 pin = mctrack->GetMomentum();
136 
137  double out[3]; out[0] = pin.x(); out[1] = pin.y(); out[2] = pin.z();
138 
139  //rw->MasterToLocalVect(arr, out);
140  TVector3 psim(out);//, prec;
141  // FIXME: do by hand for now;
142  double p40[4], p41[4], dp4[4];
143  p40[0] = sqrt(_ME_*_ME_ + _P0_*_P0_); p40[1] = p40[2] = 0.0; p40[3] = -_P0_;
144  p41[0] = sqrt(_ME_*_ME_ + pin.Mag2()); p41[1] = pin.x(); p41[2] = pin.y(); p41[3] = pin.z();
145  //printf("%f\n", pin.Z());
146  for(unsigned iq=0; iq<4; iq++)
147  dp4[iq] = p41[iq] - p40[iq];
148  double q2 = dp4[0]*dp4[0];
149  for(unsigned iq=1; iq<4; iq++)
150  q2 -= dp4[iq]*dp4[iq];
151  double theta = 1000*(180.0 - psim.Theta() * TMath::RadToDeg())*TMath::DegToRad();//, ptsim = psim.Pt();
152  printf("%d -> %7.2f %7.2f %f\n", lqstHitArray->GetEntriesFast(), theta, psim.Mag(), q2);
153  //#if _TODAY_
154  //printf("%d %d -> %7.2f\n", b0HitArray->GetEntriesFast(), lqstHitArray->GetEntriesFast(), theta);
155  th0->Fill(theta); //pq3->Fill(psim.Mag()); pt0->Fill(ptsim);//pt0->Fill(_P0_*psim.Theta());
156  //q2h0->Fill(-q2);
157  q2h0->Fill(log(-q2));
158  //#endif
159  if (lqstHitArray->GetEntriesFast() == 2) {
160  //#if _TODAY_
161  th1->Fill(theta);
162  q2h1->Fill(log(-q2));
163  //pt1->Fill(_P0_*psim.Theta());
164  //pt1->Fill(ptsim);
165  //#endif
166 
167  // NB: will use the same _X0_ & _Y0_ for both; does not matter as long as [0] is
168  // used for slope calculation only;
169  //double XL[2], YL[2], base = 200.0;
170  //for(unsigned iq=0; iq<2; iq++) {
171  EicTrackingDigiHitOrth2D *lqsthit = lqstHitArray->At(0);
172 
173  double xl = 10*lqsthit->GetCoord(0) - _X0_, yl = 10*lqsthit->GetCoord(1) - _Y0_;
174  hxl->Fill(xl);
175 
176  //printf("XL = %7.2f, YL = %7.2f [mm]\n", xl, yl);
177  double qx = xl/sigx, qy = yl/sigy, q2d = sqrt(qx*qx+qy*qy);
178 
179  // 10 sigma in 2D;
180  if (q2d > _SIGMA_CUT_) {
181  th3->Fill(theta);
182  //q2h3->Fill(-q2);
183  q2h3->Fill(log(-q2));
184 #if _TODAY_
185  //pt3->Fill(_P0_*psim.Theta());
186  pt3->Fill(ptsim);
187 
188  thq3->Fill(theta);
189  pq3->Fill(psim.Mag());
190  //printf("--> %f %f\n",
191  pthq3->Fill(psim.Mag(), theta);
192 #endif
193  } //if
194 
195  hxl->Fill(xl);
196  //x1->Fill(xl);
197  //y1->Fill(yl);
198  //if (fabs(theta)< 1.) xx->Fill(xl);
199  //if (fabs(theta)< 1.) yy->Fill(yl);
200  //xy->Fill(xl, yl);
201 
202 #if _TODAY_
203  double Xlqst[2], Xip[2]; Xlqst[0] = XL[1]; Xlqst[1] = 1000*(XL[1] - XL[0])/base - _SX0_;
204  double Ylqst = YL[1], Yip = aYinv * Ylqst;
205  Xip[0] = Xip[1] = 0.0;
206  for(unsigned i0=0; i0<2; i0++)
207  for(unsigned i1=0; i1<2; i1++)
208  Xip[i0] += aXinv[i0][i1]*Xlqst[i1];
209 
210  //printf("%7.3f %7.3f -> %7.3f\n", theta, Xip[1], Xip[1] - theta);
211  double perc = 100*(psim.Mag()-275.)/275.;
212  {
213  //double sx = htctrack->BeamSX(), sy = htctrack->BeamSY(), p = htctrack->mMomentum;
214  double rsx = Xip[1]/1000., rsy = Yip/1000., p = 275.*(1 + 0.01*Xip[0]);
215  double norm = sqrt(1.0 + rsx*rsx + rsy*rsy);
216  double arr[3], out[3]; arr[0] = p*rsx/norm; arr[1] = p*rsy/norm; arr[2] = p/norm;
217 
218  //rw->MasterToLocalVect(arr, out);
219  prec = TVector3(arr);//out);
220  //printf("%f\n", prec.Mag());
221  }
222  //printf("%7.3f %7.3f -> %7.3f\n", perc, Xip[0], Xip[0] - perc);
223  //thres->Fill(Xip[1] - theta);
224  //thres->Fill(Yip - theta);
225  //pres->Fill( Xip[0] - perc);
226  pres->Fill(100.*(prec.Mag() - psim.Mag())/psim.Mag());
227  //printf("%f\n", 1000*(prec.Pt() - psim.Pt()));
228  dpt->Fill(1000*(prec.Pt() - psim.Pt()));
229 
230  //printf("%f\n", 1000*(XL[1] - XL[0])/base - _SX0_);
231  sx->Fill(1000*(XL[1] - XL[0])/base - _SX0_);
232  sy->Fill(1000*(YL[1] - YL[0])/base - _SY0_);
233 #endif
234  } //if
235  } //if
236  } //for ev
237 
238  q2h0->GetXaxis()->SetTitle("log(Q^{2})");
239  q2h0->GetXaxis()->SetTitleOffset(0.9);
240  q2h0->GetXaxis()->SetLabelFont(52);
241  q2h0->GetXaxis()->SetLabelSize(0.050);
242  q2h0->GetXaxis()->SetTitleFont(52);
243  q2h0->GetXaxis()->SetTitleSize(0.060);
244 
245  q2h0->GetYaxis()->SetTitle("Events");
246  q2h0->GetYaxis()->SetTitleOffset(0.8);
247  q2h0->GetYaxis()->SetLabelFont(52);
248  q2h0->GetYaxis()->SetLabelSize(0.050);
249  q2h0->GetYaxis()->SetTitleFont(52);
250  q2h0->GetYaxis()->SetTitleSize(0.060);
251 
252  q2h0->Draw();
253  q2h1->Draw("SAME");
254  q2h3->Draw("SAME");
255  gPad->SetGrid();
256 
257 #if _TODAY_
258  //phi->Draw();
259  //acc->Draw("COLZ");
260  //xy->Draw("COLZ");
261  pq3->GetXaxis()->SetTitle("Leading proton momentum [GeV/c]");
262  pq3->GetXaxis()->SetTitleOffset(0.9);
263  pq3->GetXaxis()->SetLabelFont(52);
264  pq3->GetXaxis()->SetLabelSize(0.050);
265  pq3->GetXaxis()->SetTitleFont(52);
266  pq3->GetXaxis()->SetTitleSize(0.060);
267 
268  pq3->GetYaxis()->SetTitle("Events");
269  pq3->GetYaxis()->SetTitleOffset(0.8);
270  pq3->GetYaxis()->SetLabelFont(52);
271  pq3->GetYaxis()->SetLabelSize(0.050);
272  pq3->GetYaxis()->SetTitleFont(52);
273  pq3->GetYaxis()->SetTitleSize(0.060);
274 
275  thq3->GetXaxis()->SetTitle("Spectator proton theta [mrad]");
276  thq3->GetXaxis()->SetTitleOffset(0.9);
277  thq3->GetXaxis()->SetLabelFont(52);
278  thq3->GetXaxis()->SetLabelSize(0.050);
279  thq3->GetXaxis()->SetTitleFont(52);
280  thq3->GetXaxis()->SetTitleSize(0.060);
281 
282  thq3->GetYaxis()->SetTitle("Events");
283  thq3->GetYaxis()->SetTitleOffset(0.8);
284  thq3->GetYaxis()->SetLabelFont(52);
285  thq3->GetYaxis()->SetLabelSize(0.050);
286  thq3->GetYaxis()->SetTitleFont(52);
287  thq3->GetYaxis()->SetTitleSize(0.060);
288 
289  pthq3->GetXaxis()->SetTitle("Spectator proton momentum [GeV/c]");
290  pthq3->GetXaxis()->SetTitleOffset(0.9);
291  pthq3->GetXaxis()->SetLabelFont(52);
292  pthq3->GetXaxis()->SetLabelSize(0.050);
293  pthq3->GetXaxis()->SetTitleFont(52);
294  pthq3->GetXaxis()->SetTitleSize(0.060);
295 
296  //pthq3->GetYaxis()->SetTitle("Events");
297  pthq3->GetYaxis()->SetTitle("Spectator proton theta [mrad]");
298  pthq3->GetYaxis()->SetTitleOffset(0.8);
299  pthq3->GetYaxis()->SetLabelFont(52);
300  pthq3->GetYaxis()->SetLabelSize(0.050);
301  pthq3->GetYaxis()->SetTitleFont(52);
302  pthq3->GetYaxis()->SetTitleSize(0.060);
303 
304 
305 #ifdef _THETA_PLOT_
306  th0->GetXaxis()->SetTitle("Leading proton theta [mrad]");
307  th0->GetXaxis()->SetTitleOffset(0.9);
308  th0->GetXaxis()->SetLabelFont(52);
309  th0->GetXaxis()->SetLabelSize(0.050);
310  th0->GetXaxis()->SetTitleFont(52);
311  th0->GetXaxis()->SetTitleSize(0.060);
312 
313  th0->GetYaxis()->SetTitle("Events");
314  th0->GetYaxis()->SetTitleOffset(0.8);
315  th0->GetYaxis()->SetLabelFont(52);
316  th0->GetYaxis()->SetLabelSize(0.050);
317  th0->GetYaxis()->SetTitleFont(52);
318  th0->GetYaxis()->SetTitleSize(0.060);
319  th0->Draw();
320  th2->Draw("SAME");
321  th3->Draw("SAME");
322 #else
323  pt0->GetXaxis()->SetTitle("DVCS proton P_{t} [GeV/c]");
324  pt0->GetXaxis()->SetTitleOffset(0.9);
325  pt0->GetXaxis()->SetLabelFont(52);
326  pt0->GetXaxis()->SetLabelSize(0.050);
327  pt0->GetXaxis()->SetTitleFont(52);
328  pt0->GetXaxis()->SetTitleSize(0.060);
329 
330  pt0->GetYaxis()->SetTitle("Events");
331  pt0->GetYaxis()->SetTitleOffset(0.8);
332  pt0->GetYaxis()->SetLabelFont(52);
333  pt0->GetYaxis()->SetLabelSize(0.050);
334  pt0->GetYaxis()->SetTitleFont(52);
335  pt0->GetYaxis()->SetTitleSize(0.060);
336 
337  pt0->Draw();
338  pt2->Draw("SAME");
339  pt3->Draw("SAME");
340 #endif
341 
342 #endif
343 } // acceptance()