Theta resolution--------------------------------------------------------------------------————
Misalignment constants -----------------------------------------------------------——
END (Misalignment constants) --------------------------------------------------------—
12 const int ntrksSample = 2e4;
13 const double i_TrksSimi = 100./ntrksSample;
15 double v_Dt[nParDt]={0,50,100,200,300,400,500,600};
17 double v_Da[nParDa]={0,1,3,6,9};
22 for(
int iDt=0;iDt<nParDt;iDt++){
24 for(
int iDa=0;iDa<nParDa;iDa++){
26 cout<<
"work with dt~"<<tr_sc<<
" da~"<<rt_sc<<endl;
35 double TrksBef[nS][nS];
36 double TrksAft[nS][nS];
37 double TrksDiff[nS][nS];
39 TString resname= pathG+
"/ResultsSamplesSummary_";
45 TString resname_pdf = resname+
".pdf";
46 TString resname_pdf_o = resname_pdf+
"(";
47 TString resname_pdf_c = resname_pdf+
")";
49 double mis_b[nS][6][nS];
50 double mis_a[nS][6][nS];
52 double theta_mean_b[nS][nS], theta_mean_a[nS][nS], theta_mean_r[nS][nS];
53 double theta_rms_b[nS][nS], theta_rms_a[nS][nS], theta_rms_r[nS][nS];
56 TH1D *hthetabefore =
new TH1D(
"hthetabefore",
"#theta_{MC}-#theta_{rec};#delta#theta,rad",1e2,-1e-3,1e-3);
57 TH1D *hthetaafter =
new TH1D(
"hthetaafter",
"#theta_{MC}-#theta_{rec};#delta#theta,rad",1e2,-1e-3,1e-3);
58 TH1D *hthetaref =
new TH1D(
"hthetaref",
"#theta_{MC}-#theta_{rec};#delta#theta,rad",1e2,-1e-3,1e-3);
60 for(
int ipart=0;ipart<nS;ipart++){
61 for(
int i=0;
i<nS;
i++){
62 TrksSim[
i]=(
i+1)*ntrksSample ;
64 theta_mean_b[ipart][
i] = 0;
65 theta_mean_a[ipart][
i] = 0;
66 theta_mean_r[ipart][
i] = 0;
67 theta_rms_b[ipart][
i] = 0;
68 theta_rms_a[ipart][
i] = 0;
69 theta_rms_r[ipart][
i] = 0;
81 for(
int i=0;
i<nS;
i++){
82 cout<<
"Sum of "<<
i+1<<
"sample(s)"<<endl;
83 for(
int ipart=0;ipart<nS;ipart++){
84 cout<<
"part #"<<ipart<<endl;
85 int nTrks = TrksSim[
i];
107 name3+=
"/results/Lumi_QA_before.root";
117 TFile *
f1 =
new TFile(name1,
"READ");
118 if (f1->IsZombie()) {
119 std::cout <<
"!!! Error opening file " <<name1<< std::endl;
124 hbefore = (TH1*)f1->Get(
"NearIP/hResTheta");
126 hbefore->SetName(
"hResThetaBefore");
127 hbefore->SetTitle(
"#theta resolution");
129 TrksBef[ipart][
i] = hbefore->GetEntries();
130 TrksBef[ipart][
i] *=i_TrksSimi;
140 theta_mean_b[ipart][
i] = 1e6*(hbefore->GetMean());
141 theta_rms_b[ipart][
i] = 1e6*(hbefore->GetRMS());
143 TFile *
f2 =
new TFile(name2,
"READ");
144 if (f2->IsZombie()) {
145 std::cout <<
"!!! Error opening file " <<name2<< std::endl;
149 hafter = (TH1F*)f2->Get(
"NearIP/hResTheta");
150 hafter->SetName(
"hResThetaAfter");
151 hafter->SetTitle(
"#theta resolution");
154 TrksAft[ipart][
i] = hafter->GetEntries();
155 TrksAft[ipart][
i] *=i_TrksSimi;
158 TrksDiff[ipart][
i] = (TrksAft[ipart][
i] - TrksBef[ipart][
i]);
168 theta_mean_a[ipart][
i] = 1e6*(hafter->GetMean());
169 theta_rms_a[ipart][
i] = 1e6*(hafter->GetRMS());
170 TFile *
f3 =
new TFile(name3,
"READ");
171 if (f3->IsZombie()) {
172 std::cout <<
"!!! Error opening file " <<name2<< std::endl;
176 hrefer = (TH1F*)f3->Get(
"NearIP/hResTheta");
177 hrefer->SetName(
"hResThetaRef");
178 hrefer->SetTitle(
"#theta resolution");
187 theta_mean_r[ipart][
i] = 1e6*(hrefer->GetMean());
188 theta_rms_r[ipart][
i] = 1e6*(hrefer->GetRMS());
216 namemisc+=
"mrad/Sum";
218 namemisc+=
"Samples/KnossosResults.root";
219 TFile *fmisc =
new TFile(namemisc,
"READ");
220 if (fmisc->IsZombie()) {
221 std::cout <<
"!!! Error opening file " <<namemisc<< std::endl;
224 TH2F *hmis_b_0 = (TH2F*)fmisc->Get(
"mis_before_0");
225 mis_b[ipart][0][
i] = 1e4*(hmis_b_0->ProjectionY()->GetRMS());
226 TH2F *hmis_b_1 = (TH2F*)fmisc->Get(
"mis_before_1");
227 mis_b[ipart][1][
i] = 1e4*(hmis_b_1->ProjectionY()->GetRMS());
228 TH2F *hmis_b_2 = (TH2F*)fmisc->Get(
"mis_before_2");
229 mis_b[ipart][2][
i] = 1e4*(hmis_b_2->ProjectionY()->GetRMS());
230 TH2F *hmis_b_3 = (TH2F*)fmisc->Get(
"mis_before_3");
231 mis_b[ipart][3][
i] = 1e3*(hmis_b_3->ProjectionY()->GetRMS());
232 TH2F *hmis_b_4 = (TH2F*)fmisc->Get(
"mis_before_4");
233 mis_b[ipart][4][
i] = 1e3*(hmis_b_4->ProjectionY()->GetRMS());
234 TH2F *hmis_b_5 = (TH2F*)fmisc->Get(
"mis_before_5");
235 mis_b[ipart][5][
i] = 1e3*(hmis_b_5->ProjectionY()->GetRMS());
237 TH2F *hmis_a_0 = (TH2F*)fmisc->Get(
"mis_diff_0");
238 mis_a[ipart][0][
i] = 1e4*(hmis_a_0->ProjectionY()->GetRMS());
239 TH2F *hmis_a_1 = (TH2F*)fmisc->Get(
"mis_diff_1");
240 mis_a[ipart][1][
i] = 1e4*(hmis_a_1->ProjectionY()->GetRMS());
241 TH2F *hmis_a_2 = (TH2F*)fmisc->Get(
"mis_diff_2");
242 mis_a[ipart][2][
i] = 1e4*(hmis_a_2->ProjectionY()->GetRMS());
243 TH2F *hmis_a_3 = (TH2F*)fmisc->Get(
"mis_diff_3");
244 mis_a[ipart][3][
i] = 1e3*(hmis_a_3->ProjectionY()->GetRMS());
245 TH2F *hmis_a_4 = (TH2F*)fmisc->Get(
"mis_diff_4");
246 mis_a[ipart][4][
i] = 1e3*(hmis_a_4->ProjectionY()->GetRMS());
247 TH2F *hmis_a_5 = (TH2F*)fmisc->Get(
"mis_diff_5");
248 mis_a[ipart][5][
i] = 1e3*(hmis_a_5->ProjectionY()->GetRMS());
254 double av_theta_mean_b[nS],av_theta_rms_b[nS];
255 double av_theta_mean_a[nS],av_theta_rms_a[nS];
256 double av_theta_mean_r[nS],av_theta_rms_r[nS];
257 double av_TrksBef[nS];
double av_TrksAft[nS];
double av_TrksDiff[nS];
258 double avEr_theta_mean_b[nS],avEr_theta_rms_b[nS];
259 double avEr_theta_mean_a[nS],avEr_theta_rms_a[nS];
260 double avEr_theta_mean_r[nS],avEr_theta_rms_r[nS];
261 double avEr_TrksBef[nS];
double avEr_TrksAft[nS];
double avEr_TrksDiff[nS];
271 double av_mis_a_mean[6][nS];
272 double av_mis_b_mean[6][nS];
273 double avEr_mis_a_mean[6][nS];
274 double avEr_mis_b_mean[6][nS];
276 for(
int i=0;
i<nS;
i++){
277 av_theta_mean_b[
i] = 0;
278 av_theta_rms_b[
i] = 0;
279 av_theta_mean_a[
i] = 0;
280 av_theta_rms_a[
i] = 0;
281 av_theta_mean_r[
i] = 0;
282 av_theta_rms_r[
i] = 0;
286 avEr_theta_mean_b[
i] = 0;
287 avEr_theta_rms_b[
i] =0;
288 avEr_theta_mean_a[
i] = 0;
289 avEr_theta_rms_a[
i] = 0;
290 avEr_theta_mean_r[
i] = 0;
291 avEr_theta_rms_r[
i] = 0;
308 for(
int ialc=0;ialc<6;ialc++){
309 av_mis_a_mean[ialc][
i] = 0;
310 av_mis_b_mean[ialc][
i] = 0;
311 avEr_mis_a_mean[ialc][
i] = 0;
312 avEr_mis_b_mean[ialc][
i] = 0;
316 for(
int i=0;
i<nS;
i++){
317 for(
int ipart=0;ipart<nS;ipart++){
318 cout<<
"TrksAft["<<ipart<<
"]["<<
i<<
"] = "<<TrksAft[ipart][
i]<<
" TrksBef["<<ipart<<
"]["<<
i<<
"] = "<<TrksBef[ipart][
i]
319 <<
" TrksDiff["<<ipart<<
"]["<<
i<<
"] = "<<TrksDiff[ipart][
i]<<endl;
323 av_theta_mean_b[
i]+=theta_mean_b[ipart][
i];
324 av_theta_mean_a[
i]+=theta_mean_a[ipart][
i];
325 av_theta_mean_r[
i]+=theta_mean_r[ipart][
i];
326 av_theta_rms_b[
i]+=theta_rms_b[ipart][
i];
327 av_theta_rms_a[
i]+=theta_rms_a[ipart][
i];
328 av_theta_rms_r[
i]+=theta_rms_r[ipart][
i];
329 av_TrksBef[
i]+=TrksBef[ipart][
i];
330 av_TrksAft[
i]+=TrksAft[ipart][
i];
331 av_TrksDiff[
i]+=TrksDiff[ipart][
i];
342 for(
int ialc=0;ialc<6;ialc++){
343 av_mis_a_mean[ialc][
i] += mis_a[ipart][ialc][
i];
344 av_mis_b_mean[ialc][
i] += mis_b[ipart][ialc][
i];
348 av_theta_mean_b[
i] *=i_nS;
349 av_theta_mean_a[
i] *=i_nS;
350 av_theta_mean_r[
i] *=i_nS;
351 av_theta_rms_b[
i] *=i_nS;
352 av_theta_rms_a[
i] *=i_nS;
353 av_theta_rms_r[
i] *=i_nS;
356 av_TrksDiff[
i]*=i_nS;
365 for(
int ialc=0;ialc<6;ialc++){
366 av_mis_a_mean[ialc][
i] *= i_nS;
367 av_mis_b_mean[ialc][
i] *= i_nS;
372 double i_nS_r = 1./(nS-1);
373 for(
int i=0;
i<nS;
i++){
374 for(
int ipart=0;ipart<nS;ipart++){
375 avEr_theta_mean_b[
i]+=(theta_mean_b[ipart][
i]-av_theta_mean_b[
i])*(theta_mean_b[ipart][
i]-av_theta_mean_b[
i]);
376 avEr_theta_mean_a[
i]+=(theta_mean_a[ipart][
i]-av_theta_mean_a[
i])*(theta_mean_a[ipart][i]-av_theta_mean_a[i]);
377 avEr_theta_mean_r[
i]+=(theta_mean_r[ipart][
i]-av_theta_mean_r[
i])*(theta_mean_r[ipart][i]-av_theta_mean_r[i]);
378 avEr_theta_rms_b[
i]+=(theta_rms_b[ipart][
i]-av_theta_rms_b[
i])*(theta_rms_b[ipart][i]-av_theta_rms_b[i]);
379 avEr_theta_rms_a[
i]+=(theta_rms_a[ipart][
i]-av_theta_rms_a[
i])*(theta_rms_a[ipart][i]-av_theta_rms_a[i]);
380 avEr_theta_rms_r[
i]+=(theta_rms_r[ipart][
i]-av_theta_rms_r[
i])*(theta_rms_r[ipart][i]-av_theta_rms_r[i]);
381 avEr_TrksBef[
i]+=(TrksBef[ipart][
i]-av_TrksBef[
i])*(TrksBef[ipart][i]-av_TrksBef[i]);
382 avEr_TrksAft[
i]+=(TrksAft[ipart][
i]-av_TrksAft[
i])*(TrksAft[ipart][i]-av_TrksAft[i]);
383 avEr_TrksDiff[
i]+=(TrksDiff[ipart][
i]-av_TrksDiff[
i])*(TrksDiff[ipart][i]-av_TrksDiff[i]);
392 for(
int ialc=0;ialc<6;ialc++){
393 avEr_mis_a_mean[ialc][
i] +=(mis_a[ipart][ialc][
i]-av_mis_a_mean[ialc][
i])*(mis_a[ipart][ialc][i]-av_mis_a_mean[ialc][i]);
394 avEr_mis_b_mean[ialc][
i] +=(mis_b[ipart][ialc][
i]-av_mis_b_mean[ialc][
i])*(mis_b[ipart][ialc][i]-av_mis_b_mean[ialc][i]);
397 avEr_theta_mean_b[
i] *=i_nS_r;
398 avEr_theta_mean_a[
i] *=i_nS_r;
399 avEr_theta_mean_r[
i] *=i_nS_r;
400 avEr_theta_rms_b[
i] *=i_nS_r;
401 avEr_theta_rms_a[
i] *=i_nS_r;
402 avEr_theta_rms_r[
i] *=i_nS_r;
403 avEr_TrksBef[
i] *= i_nS_r;
404 avEr_TrksAft[
i] *= i_nS_r;
405 avEr_TrksDiff[
i] *= i_nS_r;
406 avEr_theta_mean_b[
i] =
sqrt(avEr_theta_mean_b[i]);
407 avEr_theta_mean_a[
i] =
sqrt(avEr_theta_mean_a[i]);
408 avEr_theta_mean_r[
i] =
sqrt(avEr_theta_mean_r[i]);
409 avEr_theta_rms_b[
i] =
sqrt(avEr_theta_rms_b[i]);
410 avEr_theta_rms_a[
i] =
sqrt(avEr_theta_rms_a[i]);
411 avEr_theta_rms_r[
i] =
sqrt(avEr_theta_rms_r[i]);
412 avEr_TrksBef[
i] =
sqrt(avEr_TrksBef[i]);
413 avEr_TrksAft[
i] =
sqrt(avEr_TrksAft[i]);
414 avEr_TrksDiff[
i] =
sqrt(avEr_TrksDiff[i]);
432 for(
int ialc=0;ialc<6;ialc++){
433 avEr_mis_a_mean[ialc][
i] *=i_nS_r;
434 avEr_mis_b_mean[ialc][
i] *=i_nS_r;
435 avEr_mis_a_mean[ialc][
i] =
sqrt(avEr_mis_a_mean[ialc][i]);
436 avEr_mis_b_mean[ialc][
i] =
sqrt(avEr_mis_b_mean[ialc][i]);
444 TGraphErrors *gr_theta_b =
new TGraphErrors(nS,TrksSim, av_theta_mean_b,0, avEr_theta_mean_b);
445 gr_theta_b->SetMarkerStyle(20);
446 gr_theta_b->SetMarkerColor(kGreen-3);
447 gr_theta_b->SetMarkerSize(2.5);
448 TGraphErrors *gr_theta_a =
new TGraphErrors(nS,TrksSim, av_theta_mean_a,0, avEr_theta_mean_a);
449 gr_theta_a->SetMarkerStyle(29);
450 gr_theta_a->SetMarkerColor(kOrange+7);
451 gr_theta_a->SetMarkerSize(2.5);
452 TGraphErrors *gr_theta_r =
new TGraphErrors(nS,TrksSim, av_theta_mean_r,0, avEr_theta_mean_r);
453 gr_theta_r->SetMarkerStyle(21);
454 gr_theta_r->SetMarkerColor(15);
455 gr_theta_r->SetMarkerSize(2.5);
458 TMultiGraph *mgr_theta =
new TMultiGraph();
460 mgr_theta->Add(gr_theta_r);
461 mgr_theta->Add(gr_theta_a);
463 mgr_theta->Draw(
"AP");
464 mgr_theta->GetXaxis()->SetTitle(
"Average number of rec. trks");
465 mgr_theta->GetYaxis()->SetTitle(
"#theta_{mean}, #murad");
466 c1.Print(resname_pdf_o);
469 TGraphErrors *gr_theta_rms_b =
new TGraphErrors(nS,TrksSim,av_theta_rms_b,0,avEr_theta_rms_b);
470 gr_theta_rms_b->SetMarkerStyle(20);
471 gr_theta_rms_b->SetMarkerColor(kGreen-3);
472 gr_theta_rms_b->SetMarkerSize(2.5);
473 TGraphErrors *gr_theta_rms_a =
new TGraphErrors(nS,TrksSim,av_theta_rms_a,0,avEr_theta_rms_a);
474 gr_theta_rms_a->SetMarkerStyle(29);
475 gr_theta_rms_a->SetMarkerColor(kOrange+7);
476 gr_theta_rms_a->SetMarkerSize(2.5);
477 TGraphErrors *gr_theta_rms_r =
new TGraphErrors(nS,TrksSim,av_theta_rms_r,0,avEr_theta_rms_a);
478 gr_theta_rms_r->SetMarkerStyle(21);
479 gr_theta_rms_r->SetMarkerColor(15);
480 gr_theta_rms_r->SetMarkerSize(2.5);
481 TMultiGraph *mgr_theta_rms =
new TMultiGraph();
483 mgr_theta_rms->Add(gr_theta_rms_r);
484 mgr_theta_rms->Add(gr_theta_rms_a);
485 mgr_theta_rms->Draw(
"AP");
486 mgr_theta_rms->GetXaxis()->SetTitle(
"Average number of rec. trks per sector");
487 mgr_theta_rms->GetYaxis()->SetTitle(
"#theta_{#sigma}, #murad");
488 c1.Print(resname_pdf_o);
545 TGraphErrors *gr_mis_b[6];
546 TGraphErrors *gr_mis_a[6];
547 TMultiGraph *mgr_mis[6];
548 for(
int ialc=0;ialc<6;ialc++){
549 gr_mis_b[ialc] =
new TGraphErrors(nS,TrksSim,av_mis_b_mean[ialc],0,avEr_mis_b_mean[ialc]);
550 gr_mis_b[ialc]->SetMarkerStyle(20);
551 gr_mis_b[ialc]->SetMarkerColor(kGreen-3);
552 gr_mis_b[ialc]->SetMarkerSize(2.5);
553 gr_mis_a[ialc] =
new TGraphErrors(nS,TrksSim,av_mis_a_mean[ialc],0,avEr_mis_a_mean[ialc]);
554 gr_mis_a[ialc]->SetMarkerStyle(29);
555 gr_mis_a[ialc]->SetMarkerColor(kOrange+7);
556 gr_mis_a[ialc]->SetMarkerSize(2.5);
557 mgr_mis[ialc] =
new TMultiGraph();
558 mgr_mis[ialc]->Add(gr_mis_b[ialc]);
559 mgr_mis[ialc]->Add(gr_mis_a[ialc]);
560 mgr_mis[ialc]->Draw(
"AP");
561 mgr_mis[ialc]->GetXaxis()->SetTitle(
"Average number of rec. trks ");
562 if(ialc==0) mgr_mis[ialc]->GetYaxis()->SetTitle(
"#Delta_{x}, #mum");
563 if(ialc==1) mgr_mis[ialc]->GetYaxis()->SetTitle(
"#Delta_{y}, #mum");
564 if(ialc==2) mgr_mis[ialc]->GetYaxis()->SetTitle(
"#Delta_{z}, #mum");
565 if(ialc==3) mgr_mis[ialc]->GetYaxis()->SetTitle(
"#Delta_{#alpha}, mrad");
566 if(ialc==4) mgr_mis[ialc]->GetYaxis()->SetTitle(
"#Delta_{#beta}, mrad");
567 if(ialc==5) mgr_mis[ialc]->GetYaxis()->SetTitle(
"#Delta_{#gamma}, mrad");
568 if(ialc<2 || ialc==5) c1.Print(resname_pdf_o);
571 TMultiGraph *mgr_stat =
new TMultiGraph();
572 TGraphErrors *grstatBef =
new TGraphErrors(nS,TrksSim,av_TrksBef,0,avEr_TrksBef);
573 grstatBef->SetMarkerStyle(24);
574 grstatBef->SetMarkerSize(2.5);
575 TGraphErrors *grstatAft =
new TGraphErrors(nS,TrksSim,av_TrksAft,0,avEr_TrksAft);
576 grstatAft->SetMarkerStyle(22);
577 grstatAft->SetMarkerSize(2.5);
578 grstatAft->SetTitle(
"Average number of rec. trks (empty = before Knossos; filled = after Knossos)");
579 mgr_stat->Add(grstatBef);
580 mgr_stat->Add(grstatAft);
581 mgr_stat->Draw(
"AP");
582 mgr_stat->GetXaxis()->SetTitle(
"Number of sim. trks");
583 mgr_stat->GetYaxis()->SetTitle(
"N_{REC}/N_{SIM}, % [empty = N^{ before}_{REC}; filled = N^{ after}_{REC}]");
585 c1.Print(resname_pdf_o);
587 TGraphErrors *grstatDiff =
new TGraphErrors(nS,TrksSim,av_TrksDiff,0,avEr_TrksDiff);
588 grstatDiff->SetMarkerStyle(22);
589 grstatDiff->SetMarkerSize(2.5);
590 grstatDiff->Draw(
"ALP");
591 grstatDiff->GetXaxis()->SetTitle(
"Number of sim. trks");
592 grstatDiff->GetYaxis()->SetTitle(
"(N^{ after}_{REC} - N^{ before}_{REC})/N_{SIM}, %");
594 c1.Print(resname_pdf_c);
597 TFile *
f =
new TFile(out,
"RECREATE");
607 gr_theta_rms_b->Write();
608 gr_theta_rms_r->Write();
609 gr_theta_rms_a->Write();
610 mgr_theta_rms->Write();
612 for(
int ialc=0;ialc<6;ialc++){
613 gr_mis_b[ialc]->Write();
614 gr_mis_a[ialc]->Write();
615 mgr_mis[ialc]->Write();
friend F32vec4 sqrt(const F32vec4 &a)