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 _41_GEV_
3 //#define _100_GEV_
4 //#define _275_GEV_
5 
6 //#define _THETA_PLOT_
7 
8 // Assume 41 x 5 GeV case in table 3.3;
9 #ifdef _41_GEV_
10 // Beam emittance; converted to [mm];
11 #define _EX_ (45.9E-6)
12 #define _EY_ ( 2.6E-6)
13 
14 // Beam energy spread;
15 #define _DPP_ (10.3E-4)
16 #endif
17 // Assume 100 x 10 GeV case in table 3.3;
18 #ifdef _100_GEV_
19 #define _EX_ (33.8E-6)
20 #define _EY_ ( 1.8E-6)
21 
22 #define _DPP_ ( 9.1E-4)
23 #endif
24 // Assume 275 x 18 GeV case in table 3.3;
25 #ifdef _275_GEV_
26 #define _EX_ (12.4E-6)
27 #define _EY_ ( 5.5E-6)
28 
29 #define _DPP_ ( 4.6E-4)
30 #endif
31 
32 //
33 // Assume the rest is the same for all three cases;
34 //
35 
36 // Central trajectory offset; must be tuned by hand;
37 #define _X0_ (-0.73)
38 #define _Y0_ ( 0.00)
39 #define _SX0_ (38.76)
40 #define _SY0_ ( 0.00)
41 
42 // This cut will limit acceptance at small theta;
43 #define _SIGMA_CUT_ (10.0)
44 
45 // Dispersion at RP location in [mm/10^-3];
46 #define _DX_ (-0.19)
47 
48 // Beta functions around Z~30m for this beam tune; converted to [mm];
49 #define _BETA_X_ (475E3)
50 #define _BETA_Y_ (520E3)
51 
52 {
53  // Load basic libraries;
54  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
55 
56  //+gStyle->SetOptStat(0);
57  gStyle->SetPadLeftMargin(0.13);
58  gStyle->SetPadBottomMargin(0.13);
59 
60  // Input simulated & reconstructed files;
61  TFile *ff = new TFile("simulation.root");
62  TTree *cbmsim = ff->Get("cbmsim");
63  cbmsim->AddFriend("cbmsim", "digitization.root");
64 
65  TClonesArray *rpHitArray = new TClonesArray("EicTrackingDigiHitOrth2D");
66  cbmsim->SetBranchAddress("RpTrackingDigiHit", &rpHitArray);
67  TClonesArray *b0HitArray = new TClonesArray("EicTrackingDigiHitOrth2D");
68  cbmsim->SetBranchAddress("B0trackerTrackingDigiHit", &b0HitArray);
69  TClonesArray *mcTrackArray = new TClonesArray("PndMCTrack");
70  cbmsim->SetBranchAddress("MCTrack", &mcTrackArray);
71 
72  TCanvas *c1 = new TCanvas("c1", "c1", 0, 0, 800, 500);
73 
74  TH1D *th0 = new TH1D("th0", "", 250, 0., 25.);
75  th0->SetMinimum(0); th0->SetLineWidth(2); th0->SetLineColor(kBlack);
76  TH1D *th1 = new TH1D("th1", "", 250, 0., 25.); th1->SetLineWidth(2); th1->SetLineColor(kGreen);
77  TH1D *th2 = new TH1D("th2", "", 250, 0., 25.); th2->SetLineWidth(2); th2->SetLineColor(kRed);
78  TH1D *th3 = new TH1D("th3", "", 250, 0., 25.); th3->SetLineWidth(2); th3->SetLineColor(kBlue);
79 
80  //TH1D *pres = new TH1D("pres", "", 100, -10., 10.);
81  //TH1D *dpt = new TH1D("dpt", "dpt", 100, -300., 300.);
82  //TH1D *hxl = new TH1D("hxl", "" , 300, -300., 300.0);
83 
84  TH1D *pt0 = new TH1D("pt0", "" ,150, 0., 1.5);
85  pt0->SetLineColor(kBlack); pt0->SetMinimum(0); pt0->SetLineWidth(1);
86  TH1D *pt1 = new TH1D("pt1", "", 150, 0., 1.5); pt1->SetLineWidth(3); pt1->SetLineColor(kGreen);
87  TH1D *pt2 = new TH1D("pt2", "", 150, 0., 1.5); pt2->SetLineWidth(3); pt2->SetLineColor(kRed);
88  TH1D *pt3 = new TH1D("pt3", "", 150, 0., 1.5); pt3->SetLineWidth(3); pt3->SetLineColor(kBlue);
89 
90  // For test purposes only;
91  TH1D *tst = new TH1D("tst", "", 10000, -300., 300.);
92 
93  //TH1D *hq2d = new TH1D("hq2d", "hq2d", 100, 0., 100.);
94  //TH1D *phi= new TH1D("phi","phi",100, -3.14159, 3.14159);
95  //TH2D *acc= new TH2D("acc","acc",100, -10., 10., 100, -10., 10.);
96  //TH1D *x1 = new TH1D("x1", "x1", 3000, -30., 30.);
97  //TH1D *sx = new TH1D("sx", "sx", 1000, -10., 10.);
98  //TH1D *y1 = new TH1D("y1", "y1", 2000, -20., 20.);
99  //TH1D *sy = new TH1D("sy", "sy", 1000, -10., 10.);
100  //TH2D *xy = new TH2D("xy", "xy", 50, -120., 120., 50, -120., 120.);
101 
102  TGeoRotation *rw = new TGeoRotation();
103  double _22mrad_ = 0.022 * TMath::RadToDeg();
104  rw->RotateY(_22mrad_);
105 
106  // FIXME: these are the values tuned for the particular optics of course;
107  //double aX[2][2] = {{-1.810, 24.070}, {-0.165, 0.240}}, aXinv[2][2], aYinv = 1/7.725;
108  //double det = aX[0][0]*aX[1][1] - aX[0][1]*aX[1][0]; assert(det);
109  //aXinv[0][0] = aX[1][1]/det;
110  //aXinv[0][1] = -aX[0][1]/det;
111  //aXinv[1][0] = -aX[1][0]/det;
112  //aXinv[1][1] = aX[0][0]/det;
113 
114  // Beam 1 sigma size at the RP location;
115  double sigx = sqrt(_EX_*_BETA_X_), sigy = sqrt(_EY_*_BETA_Y_);
116  double sigp = _DX_ * _DPP_ * 1E3;
117  sigx = sqrt(sigx*sigx + sigp*sigp);
118 
119  // Loop through all events; NB: for box-generated events without secondaries
120  // could simply use cbmsim->Project() as well; in general EicEventAssembler
121  // should be used for "true" physics events;
122  int nEvents = cbmsim->GetEntries();
123  for(unsigned ev=0; ev<nEvents; ev++) {
124  cbmsim->GetEntry(ev);
125 
126  assert(mcTrackArray->GetEntriesFast() == 1);
127  PndMCTrack *mctrack = mcTrackArray->At(0);
128 
129  if (mctrack->GetPdgCode() == 2212 && mctrack->GetMotherID() == -1) {
130  TVector3 pin = mctrack->GetMomentum();
131  double arr[3], out[3]; arr[0] = pin.x(); arr[1] = pin.y(); arr[2] = pin.z();
132 
133  rw->MasterToLocalVect(arr, out);
134  TVector3 psim(out), prec;
135  double theta = 1000*psim.Theta(), ptsim = psim.Pt();
136  printf("%d %d -> %7.2f %7.2f\n", b0HitArray->GetEntriesFast(), rpHitArray->GetEntriesFast(), theta, psim.Mag());
137  th0->Fill(theta); /*pq3->Fill(psim.Mag());*/ pt0->Fill(ptsim);
138  if (!b0HitArray->GetEntriesFast() && rpHitArray->GetEntriesFast() == 2) {
139  th1->Fill(theta);
140  pt1->Fill(ptsim);
141 
142  // NB: will use the same _X0_ & _Y0_ for both; does not matter as long as [0] is
143  // used for slope calculation only;
144  //double XL[2], YL[2], base = 200.0;
145  for(unsigned iq=0; iq<2; iq++) {
146  EicTrackingDigiHitOrth2D *rphit = rpHitArray->At(iq);
147 
148  double xl = 10*rphit->_GetCoord(0) - _X0_, yl = 10*rphit->_GetCoord(1) - _Y0_;
149 
150  //XL[iq] = xl; YL[iq] = yl;
151 
152  // NB: prefer to work with hit#1, for historic reasons;
153  if (iq == 1) {
154  //printf("XL = %7.2f, YL = %7.2f [mm]\n", xl, yl);
155  double qx = xl/sigx, qy = yl/sigy, q2d = sqrt(qx*qx+qy*qy);
156 
157  // 10 sigma in 2D;
158  if (q2d > _SIGMA_CUT_) {
159  th3->Fill(theta);
160  pt3->Fill(ptsim);
161 
162  //thq3->Fill(theta);
163  //pq3->Fill(psim.Mag());
164  //pthq3->Fill(psim.Mag(), theta);
165  } //if
166 
167  //hxl->Fill(xl);
168  //x1->Fill(xl);
169  //y1->Fill(yl);
170  //if (fabs(theta)< 1.) xx->Fill(xl);
171  //if (fabs(theta)< 1.) yy->Fill(yl);
172  //xy->Fill(xl, yl);
173 
174  tst->Fill(10.*rphit->_GetCoord(0));
175  } //if
176  } //for iq
177 
178  //double Xrp[2], Xip[2]; Xrp[0] = XL[1]; Xrp[1] = 1000*(XL[1] - XL[0])/base - _SX0_;
179  //double Yrp = YL[1], Yip = aYinv * Yrp;
180  //Xip[0] = Xip[1] = 0.0;
181  //for(unsigned i0=0; i0<2; i0++)
182  //for(unsigned i1=0; i1<2; i1++)
183  // Xip[i0] += aXinv[i0][i1]*Xrp[i1];
184 
185  //printf("%7.3f %7.3f -> %7.3f\n", theta, Xip[1], Xip[1] - theta);
186  //double perc = 100*(psim.Mag()-275.)/275.;
187  //{
188  //double rsx = Xip[1]/1000., rsy = Yip/1000., p = 275.*(1 + 0.01*Xip[0]);
189  //double norm = sqrt(1.0 + rsx*rsx + rsy*rsy);
190  //double arr[3], out[3]; arr[0] = p*rsx/norm; arr[1] = p*rsy/norm; arr[2] = p/norm;
191 
192  //rw->MasterToLocalVect(arr, out);
193  //prec = TVector3(arr);//out);
194  //printf("%f\n", prec.Mag());
195  //}
196  //printf("%7.3f %7.3f -> %7.3f\n", perc, Xip[0], Xip[0] - perc);
197  //thres->Fill(Xip[1] - theta);
198  //thres->Fill(Yip - theta);
199  //pres->Fill( Xip[0] - perc);
200  //pres->Fill(100.*(prec.Mag() - psim.Mag())/psim.Mag());
201  //printf("%f\n", 1000*(prec.Pt() - psim.Pt()));
202  //dpt->Fill(1000*(prec.Pt() - psim.Pt()));
203 
204  //printf("%f\n", 1000*(XL[1] - XL[0])/base - _SX0_);
205  //sx->Fill(1000*(XL[1] - XL[0])/base - _SX0_);
206  //sy->Fill(1000*(YL[1] - YL[0])/base - _SY0_);
207  } //if
208 
209  if (b0HitArray->GetEntriesFast() == 4) {
210  th2->Fill(theta);
211  pt2->Fill(ptsim);
212  } //if
213  } //if
214  } //for ev
215 
216 #if _LATER_
217  pq3->GetXaxis()->SetTitle("Spectator proton momentum [GeV/c]");
218  pq3->GetXaxis()->SetTitleOffset(0.9);
219  pq3->GetXaxis()->SetLabelFont(52);
220  pq3->GetXaxis()->SetLabelSize(0.050);
221  pq3->GetXaxis()->SetTitleFont(52);
222  pq3->GetXaxis()->SetTitleSize(0.060);
223 
224  pq3->GetYaxis()->SetTitle("Events");
225  pq3->GetYaxis()->SetTitleOffset(1.0);
226  pq3->GetYaxis()->SetLabelFont(52);
227  pq3->GetYaxis()->SetLabelSize(0.050);
228  pq3->GetYaxis()->SetTitleFont(52);
229  pq3->GetYaxis()->SetTitleSize(0.060);
230 
231  thq3->GetXaxis()->SetTitle("Spectator proton theta [mrad]");
232  thq3->GetXaxis()->SetTitleOffset(0.9);
233  thq3->GetXaxis()->SetLabelFont(52);
234  thq3->GetXaxis()->SetLabelSize(0.050);
235  thq3->GetXaxis()->SetTitleFont(52);
236  thq3->GetXaxis()->SetTitleSize(0.060);
237 
238  thq3->GetYaxis()->SetTitle("Events");
239  thq3->GetYaxis()->SetTitleOffset(1.0);
240  thq3->GetYaxis()->SetLabelFont(52);
241  thq3->GetYaxis()->SetLabelSize(0.050);
242  thq3->GetYaxis()->SetTitleFont(52);
243  thq3->GetYaxis()->SetTitleSize(0.060);
244 
245  pthq3->GetXaxis()->SetTitle("Spectator proton momentum [GeV/c]");
246  pthq3->GetXaxis()->SetTitleOffset(0.9);
247  pthq3->GetXaxis()->SetLabelFont(52);
248  pthq3->GetXaxis()->SetLabelSize(0.050);
249  pthq3->GetXaxis()->SetTitleFont(52);
250  pthq3->GetXaxis()->SetTitleSize(0.060);
251 
252  //pthq3->GetYaxis()->SetTitle("Events");
253  pthq3->GetYaxis()->SetTitle("Spectator proton theta [mrad]");
254  pthq3->GetYaxis()->SetTitleOffset(0.8);
255  pthq3->GetYaxis()->SetLabelFont(52);
256  pthq3->GetYaxis()->SetLabelSize(0.050);
257  pthq3->GetYaxis()->SetTitleFont(52);
258  pthq3->GetYaxis()->SetTitleSize(0.060);
259 #endif
260 
261 #ifdef _THETA_PLOT_
262  th0->GetXaxis()->SetTitle("Leading proton theta [mrad]");
263  th0->GetXaxis()->SetTitleOffset(0.9);
264  th0->GetXaxis()->SetLabelFont(52);
265  th0->GetXaxis()->SetLabelSize(0.050);
266  th0->GetXaxis()->SetTitleFont(52);
267  th0->GetXaxis()->SetTitleSize(0.060);
268 
269  th0->GetYaxis()->SetTitle("Events");
270  th0->GetYaxis()->SetTitleOffset(0.8);
271  th0->GetYaxis()->SetLabelFont(52);
272  th0->GetYaxis()->SetLabelSize(0.050);
273  th0->GetYaxis()->SetTitleFont(52);
274  th0->GetYaxis()->SetTitleSize(0.060);
275  th0->Draw();
276  th2->Draw("SAME");
277  th3->Draw("SAME");
278 #else
279  pt0->GetXaxis()->SetTitle("DVCS proton P_{t} [GeV/c]");
280  pt0->GetXaxis()->SetTitleOffset(0.9);
281  pt0->GetXaxis()->SetLabelFont(52);
282  pt0->GetXaxis()->SetLabelSize(0.050);
283  pt0->GetXaxis()->SetTitleFont(52);
284  pt0->GetXaxis()->SetTitleSize(0.060);
285 
286  pt0->GetYaxis()->SetTitle("Events");
287  pt0->GetYaxis()->SetTitleOffset(0.8);
288  pt0->GetYaxis()->SetLabelFont(52);
289  pt0->GetYaxis()->SetLabelSize(0.050);
290  pt0->GetYaxis()->SetTitleFont(52);
291  pt0->GetYaxis()->SetTitleSize(0.060);
292 
293  pt0->Draw();
294  pt2->Draw("SAME");
295  pt3->Draw("SAME");
296 #endif
297 
298  tst->Draw();
299 
300  gPad->SetGrid();
301 } // acceptance()