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)