FairRoot/PandaRoot
SamplesSummaryAlign.C
Go to the documentation of this file.
1 #include <TFile.h>
2 #include <TH1.h>
3 #include <TH2.h>
4 #include <TString.h>
5 #include <TCanvas.h>
6 
7 #include <sstream>
8 using namespace std;
9 //void TOTALSummaryAlign(TString pathG="/panda/pandaroot/macro/lmd/tmpOutputAlignTMP/", double tr_sc=200, double rt_sc=9)
10 int SamplesSummaryAlign(TString pathG="/home/karavdin/datastorage/AlignmentLMDpixel/LargeSamplesStudyMay/BOX/mom_15/", double tr_sc=300, double rt_sc=3)
11 {
12  const int ntrksSample = 2e4;
13  const double i_TrksSimi = 100./ntrksSample; //relative to simulated in %
14  const int nParDt=8;
15  double v_Dt[nParDt]={0,50,100,200,300,400,500,600};
16  const int nParDa=5;
17  double v_Da[nParDa]={0,1,3,6,9};
18  // const int nParDt=1;
19  // double v_Dt[nParDt]={600};
20  // const int nParDa=1;
21  // double v_Da[nParDa]={0};
22  for(int iDt=0;iDt<nParDt;iDt++){
23  tr_sc = v_Dt[iDt];
24  for(int iDa=0;iDa<nParDa;iDa++){
25  rt_sc= v_Da[iDa];
26  cout<<"work with dt~"<<tr_sc<<" da~"<<rt_sc<<endl;
27  const int nS=10;
28  // const int nS=3;
29  double TrksSim[nS];
30  //double TrksSim[nS]={10000, 100000, 200000, 400000, 600000, 800000, 1000000};
31  // const int nS=3;
32  // double TrksSim[nS]={1000,5000,10000};
33 
34  // double TrksSim[nS]={10, 100, 200, 300, 400, 500, 700, 1000};
35  double TrksBef[nS][nS];
36  double TrksAft[nS][nS];
37  double TrksDiff[nS][nS];
38  //How to save data
39  TString resname= pathG+"/ResultsSamplesSummary_";
40  resname+=tr_sc;
41  resname+="um_";
42  resname+=rt_sc;
43  resname+="mrad";
44 
45  TString resname_pdf = resname+".pdf";
46  TString resname_pdf_o = resname_pdf+"(";
47  TString resname_pdf_c = resname_pdf+")";
48 
49  double mis_b[nS][6][nS];
50  double mis_a[nS][6][nS];
51 
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];
54  // double theta_mean_b_sec[nS][10][nS], theta_mean_a_sec[nS][10][nS], theta_mean_r_sec[nS][10][nS];
55  // double theta_rms_b_sec[nS][10][nS], theta_rms_a_sec[nS][10][nS], theta_rms_r_sec[nS][10][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);
59 
60  for(int ipart=0;ipart<nS;ipart++){
61  for(int i=0;i<nS;i++){
62  TrksSim[i]=(i+1)*ntrksSample ; //each sample contains ntrksSample events with 1trk/event
63 
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;
70  // for(int isector=0;isector<10;isector++){
71  // theta_mean_b_sec[ipart][isector][i] = 0;
72  // theta_mean_a_sec[ipart][isector][i] = 0;
73  // theta_mean_r_sec[ipart][isector][i] = 0;
74  // theta_rms_b_sec[ipart][isector][i] = 0;
75  // theta_rms_a_sec[ipart][isector][i] =0;
76  // theta_rms_r_sec[ipart][isector][i] = 0;
77  // }
78  }
79  }
80 
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];
86  //Where search for files
87  TString path = pathG+"/";
88  path+=tr_sc;
89  path+="mkm_";
90  path+=rt_sc;
91  path+="mrad/Sample";
92  path +=ipart;
93  path+="/results/";
94 
95 
96 
98  TString name1 = path + "/Lumi_QA_before.root";
99  TString name2 = path + "/Lumi_QA_after.root";
100  // TString name2 = path + "/Lumi_QA_after_";
101  // name2 +=i;
102  // name2 +="SumAlign.root";
103 
104 
105  TString name3 = pathG+"/0mkm_0mrad/Sample";
106  name3+= ipart;
107  name3+="/results/Lumi_QA_before.root";
108 
109  // TH1F *hbeforetmp; TH1F *haftertmp;
110  TH1 *hbefore;
111  TH1 *hafter;
112  TH1 *hrefer;
113  // TNtuple *nsectors_b;
114  // TNtuple *nsectors_a;
115  // TNtuple *nsectors_r;
116 
117  TFile *f1 = new TFile(name1,"READ");
118  if (f1->IsZombie()) {
119  std::cout << "!!! Error opening file " <<name1<< std::endl;
120  continue;
121  // return;
122  }
123 
124  hbefore = (TH1*)f1->Get("NearIP/hResTheta");
125  //hbefore->Print();
126  hbefore->SetName("hResThetaBefore");
127  hbefore->SetTitle("#theta resolution");
128  // hbefore->Print();
129  TrksBef[ipart][i] = hbefore->GetEntries();
130  TrksBef[ipart][i] *=i_TrksSimi;
131  // // if(i>0) TrksBef[ipart][i] += TrksBef[ipart][i-1];
132 
133  // // cout<<" Trks["<<ipart<<"]["<<i<<"] = "<<Trks[ipart][i]<<endl;
134  // TF1 *funrth_b = new TF1("fitrth_b","gaus",-0.01,0.01);
135  // funrth_b->SetParameters(100,0,3e-3);
136  // funrth_b->SetParNames("Constant","Mean","Sigma");
137  // hbefore->Fit(funrth_b);
138  // theta_mean_b[ipart][i] = 1e6*(funrth_b->GetParameter("Mean"));
139  // theta_rms_b[ipart][i] = 1e6*(funrth_b->GetParameter("Sigma"));
140  theta_mean_b[ipart][i] = 1e6*(hbefore->GetMean());
141  theta_rms_b[ipart][i] = 1e6*(hbefore->GetRMS());
142  // nsectors_b = (TNtuple*)f1->Get("nsectors");
143  TFile *f2 = new TFile(name2,"READ");
144  if (f2->IsZombie()) {
145  std::cout << "!!! Error opening file " <<name2<< std::endl;
146  continue;
147  // return;
148  }
149  hafter = (TH1F*)f2->Get("NearIP/hResTheta");
150  hafter->SetName("hResThetaAfter");
151  hafter->SetTitle("#theta resolution");
152  // Trks[ipart][i] = hafter->GetEntries();
153  // if(i>0) Trks[ipart][i] += Trks[ipart][i-1];
154  TrksAft[ipart][i] = hafter->GetEntries();
155  TrksAft[ipart][i] *=i_TrksSimi;
156  // if(i>0) TrksAft[ipart][i] += TrksAft[ipart][i-1];
157 
158  TrksDiff[ipart][i] = (TrksAft[ipart][i] - TrksBef[ipart][i]);
159 
160  // TF1 *funrth_a = new TF1("fitrth_a","gaus",-0.01,0.01);
161  // funrth_a->SetParameters(100,0,3e-3);
162  // funrth_a->SetParNames("Constant","Mean","Sigma");
163  // hafter->Fit(funrth_a);
164 
165  // theta_mean_a[ipart][i] = 1e6*(funrth_a->GetParameter("Mean"));
166  // theta_rms_a[ipart][i] = 1e6*(funrth_a->GetParameter("Sigma"));
167  // // nsectors_a = (TNtuple*)f2->Get("nsectors");
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;
173  continue;
174  // return;
175  }
176  hrefer = (TH1F*)f3->Get("NearIP/hResTheta");
177  hrefer->SetName("hResThetaRef");
178  hrefer->SetTitle("#theta resolution");
179 
180  // TF1 *funrth_r = new TF1("fitrth_r","gaus",-0.01,0.01);
181  // funrth_r->SetParameters(100,0,3e-3);
182  // funrth_r->SetParNames("Constant","Mean","Sigma");
183  // hrefer->Fit(funrth_r);
184  // theta_mean_r[ipart][i] = 1e6*(funrth_r->GetParameter("Mean"));
185  // theta_rms_r[ipart][i] = 1e6*(funrth_r->GetParameter("Sigma"));
186  // nsectors_r = (TNtuple*)f3->Get("nsectors");
187  theta_mean_r[ipart][i] = 1e6*(hrefer->GetMean());
188  theta_rms_r[ipart][i] = 1e6*(hrefer->GetRMS());
189 
190  // for(int isector=0;isector<10;isector++){
191  // TString cond = "sector==";
192  // cond +=isector;
193  // nsectors_b->Project("hthetabefore","thetares",cond.Data());
194  // nsectors_a->Project("hthetaafter","thetares",cond.Data());
195  // nsectors_r->Project("hthetaref","thetares",cond.Data());
196  // hthetabefore->Fit(funrth_b);
197  // hthetaafter->Fit(funrth_a);
198  // hthetaref->Fit(funrth_r);
199  // theta_mean_b_sec[ipart][isector][i] = 1e6*(funrth_b->GetParameter("Mean"));
200  // theta_rms_b_sec[ipart][isector][i] = 1e6*(funrth_b->GetParameter("Sigma"));
201  // theta_mean_a_sec[ipart][isector][i] = 1e6*(funrth_a->GetParameter("Mean"));
202  // theta_rms_a_sec[ipart][isector][i] = 1e6*(funrth_a->GetParameter("Sigma"));
203  // theta_mean_r_sec[ipart][isector][i] = 1e6*(funrth_r->GetParameter("Mean"));
204  // theta_rms_r_sec[ipart][isector][i] = 1e6*(funrth_r->GetParameter("Sigma"));
205  // }
206 
208  f1->Close();
209  f2->Close();
210  f3->Close();
212  TString namemisc = pathG+"/";
213  namemisc+=tr_sc;
214  namemisc+="mkm_";
215  namemisc+=rt_sc;
216  namemisc+="mrad/Sum";
217  namemisc +=i;
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;
222  // return;
223  }
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());
236 
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());
250  fmisc->Close();
251  }
252  }
253  //return;
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];
262 
263  // double av_theta_mean_b_sec[10][nS],av_theta_rms_b_sec[10][nS];
264  // double av_theta_mean_a_sec[10][nS],av_theta_rms_a_sec[10][nS];
265  // double av_theta_mean_r_sec[10][nS],av_theta_rms_r_sec[10][nS];
266  // double avEr_theta_mean_b_sec[10][nS],avEr_theta_rms_b_sec[10][nS];
267  // double avEr_theta_mean_a_sec[10][nS],avEr_theta_rms_a_sec[10][nS];
268  // double avEr_theta_mean_r_sec[10][nS],avEr_theta_rms_r_sec[10][nS];
269 
270 
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];
275 
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;
283  av_TrksBef[i]=0;
284  av_TrksAft[i]=0;
285  av_TrksDiff[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;
292  avEr_TrksBef[i] = 0;
293  avEr_TrksAft[i] = 0;
294  // for(int isector=0;isector<10;isector++){
295  // av_theta_mean_b_sec[isector][i] = 0;
296  // av_theta_rms_b_sec[isector][i] = 0;
297  // av_theta_mean_a_sec[isector][i] = 0;
298  // av_theta_rms_a_sec[isector][i] = 0;
299  // av_theta_mean_r_sec[isector][i] = 0;
300  // av_theta_rms_r_sec[isector][i] = 0;
301  // avEr_theta_mean_b_sec[isector][i] = 0;
302  // avEr_theta_rms_b_sec[isector][i] = 0;
303  // avEr_theta_mean_a_sec[isector][i] = 0;
304  // avEr_theta_rms_a_sec[isector][i] = 0;
305  // avEr_theta_mean_r_sec[isector][i] = 0;
306  // avEr_theta_rms_r_sec[isector][i] = 0;
307  // }
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;
313  }
314  }
315  double i_nS = 1./nS;
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;
320  // cout<<"theta_mean_b["<<ipart<<"]["<<i<<"] = "<<theta_mean_b[ipart][i]<<endl;
321  // cout<<"theta_mean_a["<<ipart<<"]["<<i<<"] = "<<theta_mean_a[ipart][i]<<endl;
322  // cout<<"theta_mean_r["<<ipart<<"]["<<i<<"] = "<<theta_mean_r[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];
332  // cout<<" av_Trks["<<i<<"]="<<av_Trks[i]<<endl;
333  // cout<<"av_theta_mean_r["<<i<<"]="<<av_theta_mean_r[i]<<endl;
334  // for(int isector=0;isector<10;isector++){
335  // av_theta_mean_b_sec[isector][i]+=theta_mean_b_sec[ipart][isector][i];
336  // av_theta_mean_a_sec[isector][i]+=theta_mean_a_sec[ipart][isector][i];
337  // av_theta_mean_r_sec[isector][i]+=theta_mean_r_sec[ipart][isector][i];
338  // av_theta_rms_b_sec[isector][i]+=theta_rms_b_sec[ipart][isector][i];
339  // av_theta_rms_a_sec[isector][i]+=theta_rms_a_sec[ipart][isector][i];
340  // av_theta_rms_r_sec[isector][i]+=theta_rms_r_sec[ipart][isector][i];
341  // }
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];
345  }
346  }
347 
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;
354  av_TrksBef[i]*=i_nS;
355  av_TrksAft[i]*=i_nS;
356  av_TrksDiff[i]*=i_nS;
357  // for(int isector=0;isector<10;isector++){
358  // av_theta_mean_b_sec[isector][i] *=i_nS;
359  // av_theta_mean_a_sec[isector][i] *=i_nS;
360  // av_theta_mean_r_sec[isector][i] *=i_nS;
361  // av_theta_rms_b_sec[isector][i] *=i_nS;
362  // av_theta_rms_a_sec[isector][i] *=i_nS;
363  // av_theta_rms_r_sec[isector][i] *=i_nS;
364  // }
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;
368  }
369  }
370 
371 
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]);
384  // for(int isector=0;isector<10;isector++){
385  // avEr_theta_mean_b_sec[isector][i]+=(theta_mean_b_sec[ipart][isector][i]-av_theta_mean_b_sec[isector][i])*(theta_mean_b_sec[ipart][isector][i]-av_theta_mean_b_sec[isector][i]);
386  // avEr_theta_mean_a_sec[isector][i]+=(theta_mean_a_sec[ipart][isector][i]-av_theta_mean_a_sec[isector][i])*(theta_mean_a_sec[ipart][isector][i]-av_theta_mean_a_sec[isector][i]);
387  // avEr_theta_mean_r_sec[isector][i]+=(theta_mean_r_sec[ipart][isector][i]-av_theta_mean_r_sec[isector][i])*(theta_mean_r_sec[ipart][isector][i]-av_theta_mean_r_sec[isector][i]);
388  // avEr_theta_rms_b_sec[isector][i]+=(theta_rms_b_sec[ipart][isector][i]-av_theta_rms_b_sec[isector][i])*(theta_rms_b_sec[ipart][isector][i]-av_theta_rms_b_sec[isector][i]);
389  // avEr_theta_rms_a_sec[isector][i]+=(theta_rms_a_sec[ipart][isector][i]-av_theta_rms_a_sec[isector][i])*(theta_rms_a_sec[ipart][isector][i]-av_theta_rms_a_sec[isector][i]);
390  // avEr_theta_rms_r_sec[isector][i]+=(theta_rms_r_sec[ipart][isector][i]-av_theta_rms_r_sec[isector][i])*(theta_rms_r_sec[ipart][isector][i]-av_theta_rms_r_sec[isector][i]);
391  // }
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]);
395  }
396  }
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]);
415  // for(int ipart=0;ipart<nS;ipart++){
416  // for(int isector=0;isector<10;isector++){
417  // avEr_theta_mean_b_sec[isector][i] *=i_nS_r;
418  // avEr_theta_mean_a_sec[isector][i] *=i_nS_r;
419  // avEr_theta_mean_r_sec[isector][i] *=i_nS_r;
420  // avEr_theta_rms_b_sec[isector][i] *=i_nS_r;
421  // avEr_theta_rms_a_sec[isector][i] *=i_nS_r;
422  // avEr_theta_rms_r_sec[isector][i] *=i_nS_r;
423  // avEr_theta_mean_b_sec[isector][i] = sqrt(avEr_theta_mean_b_sec[isector][i]);
424  // avEr_theta_mean_a_sec[isector][i] = sqrt(avEr_theta_mean_a_sec[isector][i]);
425  // avEr_theta_mean_r_sec[isector][i] = sqrt(avEr_theta_mean_r_sec[isector][i]);
426  // avEr_theta_rms_b_sec[isector][i] = sqrt(avEr_theta_rms_b_sec[isector][i]);
427  // avEr_theta_rms_a_sec[isector][i] = sqrt(avEr_theta_rms_a_sec[isector][i]);
428  // avEr_theta_rms_r_sec[isector][i] = sqrt(avEr_theta_rms_r_sec[isector][i]);
429  // }
430  // cout<<"av_Trks["<<i<<"]="<<av_Trks[i]<<" avEr_Trks["<<i<<"]="<<avEr_Trks[i]<<endl;
431  // cout<<"av_theta_mean_r["<<i<<"]="<<av_theta_mean_r[i]<<endl;
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]);
437  }
438  }
439 
440 
441  TCanvas c1;
442  // cout<<"av_Trks = "<<av_Trks<<endl;
443  //put together theta resolution results
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);
456 
457 
458  TMultiGraph *mgr_theta = new TMultiGraph();
459  // mgr_theta->Add(gr_theta_b);
460  mgr_theta->Add(gr_theta_r);
461  mgr_theta->Add(gr_theta_a);
462 
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); //write canvas and keep the ps file open
467  c1.Clear();
468 
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();
482  // mgr_theta_rms->Add(gr_theta_rms_b);
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); //write canvas and keep the ps file open
489  c1.Clear();
490  // for(int isector=0;isector<10;isector++){
491  // //put together theta resolution results
492  // TGraphErrors *gr_theta_b_sec = new TGraphErrors(nS,TrksSim, av_theta_mean_b_sec[isector],0, avEr_theta_mean_b_sec[isector]);
493  // gr_theta_b_sec->SetMarkerStyle(20);
494  // gr_theta_b_sec->SetMarkerColor(kGreen-3);
495  // gr_theta_b_sec->SetMarkerSize(2.5);
496  // TGraphErrors *gr_theta_a_sec = new TGraphErrors(nS,TrksSim, av_theta_mean_a_sec[isector],0, avEr_theta_mean_a_sec[isector]);
497  // gr_theta_a_sec->SetMarkerStyle(29);
498  // gr_theta_a_sec->SetMarkerColor(kOrange+7);
499  // gr_theta_a_sec->SetMarkerSize(2.5);
500  // TGraphErrors *gr_theta_r_sec = new TGraphErrors(nS,TrksSim, av_theta_mean_r_sec[isector],0, avEr_theta_mean_r_sec[isector]);
501  // gr_theta_r_sec->SetMarkerStyle(21);
502  // gr_theta_r_sec->SetMarkerColor(15);
503  // gr_theta_r_sec->SetMarkerSize(2.5);
504 
505  // TMultiGraph *mgr_theta_sec = new TMultiGraph();
506  // // mgr_theta_sec->Add(gr_theta_b_sec);
507  // mgr_theta_sec->Add(gr_theta_r_sec);
508  // mgr_theta_sec->Add(gr_theta_a_sec);
509  // TString titlemean = "#theta_{mean}, #murad [sector No.";
510  // titlemean+=isector;
511  // titlemean +="]";
512  // mgr_theta_sec->Draw("AP");
513  // mgr_theta_sec->GetXaxis()->SetTitle("Average number of rec. trks ");
514  // mgr_theta_sec->GetYaxis()->SetTitle(titlemean.Data());
515  // c1.Print(resname_pdf_o); //write canvas and keep the ps file open
516  // c1.Clear();
517 
518  // TGraphErrors *gr_theta_rms_b_sec = new TGraphErrors(nS,TrksSim,av_theta_rms_b_sec[isector],0,avEr_theta_rms_b_sec[isector]);
519  // gr_theta_rms_b_sec->SetMarkerStyle(20);
520  // gr_theta_rms_b_sec->SetMarkerColor(kGreen-3);
521  // gr_theta_rms_b_sec->SetMarkerSize(2.5);
522  // TGraphErrors *gr_theta_rms_a_sec = new TGraphErrors(nS,TrksSim,av_theta_rms_a_sec[isector],0,avEr_theta_rms_a_sec[isector]);
523  // gr_theta_rms_a_sec->SetMarkerStyle(29);
524  // gr_theta_rms_a_sec->SetMarkerColor(kOrange+7);
525  // gr_theta_rms_a_sec->SetMarkerSize(2.5);
526  // TGraphErrors *gr_theta_rms_r_sec = new TGraphErrors(nS,TrksSim,av_theta_rms_r_sec[isector],0,avEr_theta_rms_a_sec[isector]);
527  // gr_theta_rms_r_sec->SetMarkerStyle(21);
528  // gr_theta_rms_r_sec->SetMarkerColor(15);
529  // gr_theta_rms_r_sec->SetMarkerSize(2.5);
530  // TMultiGraph *mgr_theta_rms_sec = new TMultiGraph();
531  // // mgr_theta_rms_sec->Add(gr_theta_rms_b_sec);
532  // mgr_theta_rms_sec->Add(gr_theta_rms_r_sec);
533  // mgr_theta_rms_sec->Add(gr_theta_rms_a_sec);
534  // mgr_theta_rms_sec->Draw("AP");
535  // mgr_theta_rms_sec->GetXaxis()->SetTitle("Average number of rec. trks ");
536  // TString titlerms = "#theta_{#sigma}, #murad [sector No.";
537  // titlerms+=isector;
538  // titlerms +="]";
539  // mgr_theta_rms_sec->GetYaxis()->SetTitle(titlerms.Data());
540  // c1.Print(resname_pdf_o); //write canvas and keep the ps file open
541  // c1.Clear();
542  // }
543 
544  //put together missalignment consts results
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); //write canvas and keep the file open
569  c1.Clear();
570  }
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}]");
584  //mgr_stat->SetTitle("Average number of rec. trks (empty = before Knossos; filled = after Knossos)");
585  c1.Print(resname_pdf_o); //write canvas and keep the file open
586  c1.Clear();
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}, %");
593  // grstatDiff->SetTitle("Diff number of rec. trks (after - before), %");
594  c1.Print(resname_pdf_c); //write canvas and close ps file
595  c1.Close();
596  TString out = resname+".root";
597  TFile *f = new TFile(out,"RECREATE");
598  // for (int g=0;g<6;g++){
599  // gr_mis_b[g]->Write();
600  // gr_mis_a[g]->Write();
601  // mgr_mis[g]->Write();
602  // }
603  gr_theta_b->Write();
604  gr_theta_r->Write();
605  gr_theta_a->Write();
606  mgr_theta->Write();
607  gr_theta_rms_b->Write();
608  gr_theta_rms_r->Write();
609  gr_theta_rms_a->Write();
610  mgr_theta_rms->Write();
611 
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();
616  }
617  mgr_stat->Write();
618  grstatDiff->Write();
619  f->Write();
620  f->Close();
621  }
622  }
623  return 0;
624 }
Int_t i
Definition: run_full.C:25
TF1 * f1
Definition: reco_analys2.C:50
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
int SamplesSummaryAlign(TString pathG="/home/karavdin/datastorage/AlignmentLMDpixel/LargeSamplesStudyMay/BOX/mom_15/", double tr_sc=300, double rt_sc=3)
TFile * f3
TFile * f
Definition: bump_analys.C:12
TFile * out
Definition: reco_muo.C:20
c1
Definition: plot_dirc.C:35
TFile * f2