FairRoot/PandaRoot
LimitsAlign.C
Go to the documentation of this file.
1 
3 #include <TFile.h>
4 #include <TH1.h>
5 #include <TH2.h>
6 #include <TString.h>
7 #include <TCanvas.h>
8 #include <TNtuple.h>
9 #include <TGraphErrors.h>
10 #include <TMultiGraph.h>
11 #include <TLegend.h>
12 #include <TMath.h>
13 #include <TStopwatch.h>
14 #include <sstream>
15 
16 #include <TFile.h>
17 #include <TH1.h>
18 #include <TH2.h>
19 #include <TString.h>
20 #include <TCanvas.h>
21 
22 #include <sstream>
23 using namespace std;
24 int LimitsAlign(TString pathG="/home/karavdin/datastorage/AlignmentLMDpixel/10data_20consts_SamplesStudyOct2013/BOX/mom_15/", double daGL=0,int iSc=0)
25 {
26 
27  const int nDt = 8;
28  double v_tr_sc[nDt]={0,50,100,200,300,400,500,600};//misalign, trans
29 
30  const int nDa = 1; //TEST
31  double v_rt_sc[nDa]={daGL};//misalign, trans
32  const int nSdGL=100;//number of data samples
33  const int nSdstep = 10; //step size for "number of used trks in Knossos"
34  const int nuSd = nSdGL/nSdstep;
35  const int nDS2 = 10;//number of samples simulated in stepC for the same case
36 
37  // for(int iSc=0;iSc<nSc;iSc++){//diff output files
38  for(int iDa=0;iDa<nDa;iDa++){//diff output files
39  double rt_sc1 = v_rt_sc[iDa];
40  double av_theta_mean_b[nuSd+1][nDt],av_theta_rms_b[nuSd+1][nDt];
41  double av_theta_mean_a[nuSd+1][nDt],av_theta_rms_a[nuSd+1][nDt];
42  double av_theta_mean_r[nuSd+1][nDt],av_theta_rms_r[nuSd+1][nDt];
43  double avEr_theta_mean_b[nuSd+1][nDt],avEr_theta_rms_b[nuSd+1][nDt];
44  double avEr_theta_mean_a[nuSd+1][nDt],avEr_theta_rms_a[nuSd+1][nDt];
45  double avEr_theta_mean_r[nuSd+1][nDt],avEr_theta_rms_r[nuSd+1][nDt];
46  double av_events[nuSd+1][nDt], avEr_events[nuSd+1][nDt];
47  double av_mis_a_mean[6][nuSd+1][nDt];
48  double av_mis_b_mean[6][nuSd+1][nDt];
49  double avEr_mis_a_mean[6][nuSd+1][nDt];
50  double avEr_mis_b_mean[6][nuSd+1][nDt];
51 
52 
53  TMultiGraph *mgr_mis[6];
54  TMultiGraph *mgr_theta = new TMultiGraph("mgr_theta_mean","");
55  TMultiGraph *mgr_theta_rms = new TMultiGraph("mgr_theta_rms","");
56  TMultiGraph *mgr_theta_b = new TMultiGraph("mgr_theta_mean_b","");
57  TMultiGraph *mgr_theta_rms_b = new TMultiGraph("mgr_theta_rms_b","");
58  TMultiGraph *mgr_events = new TMultiGraph("mgr_events","");
59 
60 
61  //How to save data
62  TString resname= pathG+"/AlignmentLimits_";
63  resname+=v_rt_sc[iDa];
64  resname+="mrad_";
65  resname+=iSc;
66  resname+="constsSample";
67  TString resname_pdf = resname+".pdf";
68  TString resname_pdf_o = resname_pdf+"(";
69  TString resname_pdf_c = resname_pdf+")";
70  for(int ijsum=1;ijsum<nuSd+1;ijsum++){ //loop over diff stat.samples (plot with diff markers)
71  int isum = 1+nSdstep*(ijsum-1);
72  for(int idt=0;idt<nDt;idt++){//loop over dt values (x-axis of plots)
73  // for(int ipart=0;ipart<nSd;ipart++){//loop over samples (take average)
74  double tr_sc = v_tr_sc[idt];
75  double rt_sc = v_rt_sc[iDa];
76  //Where search for files
77  TString path = pathG+"/";
78  path+=tr_sc;
79  path+="mkm_";
80  path+=rt_sc;
81  path+="mrad/";
82  path+=iSc;
83  path+="Sum";
84  // path+=nSd; //sample aligned on max number of tracks
85  path+=isum; //sample for alignment with sum of samples
86  path+="Samples/resultsSUM/";
87 
89  TString name1 = path + "/Lumi_QA_before.root";
90  TString name2 = path + "/Lumi_QA_after.root";
91  TString name3 = pathG+"/0mkm_0mrad/";
92  name3+=iSc;
93  name3+="Sum";
94  name3+=isum; //sample for alignment with sum of samples
95  name3+="Samples/resultsSUM/Lumi_QA_before.root";
96  TH1 *hbefore;
97  TH1 *hafter;
98  TH1 *hrefer;
99 
100  TFile *f1 = new TFile(name1,"READ");
101  if (f1->IsZombie()) {
102  std::cout << "!!! Error opening file " <<name1<< std::endl;
103  continue;
104  // return;
105  }
106 
107  hbefore = (TH1*)f1->Get("NearIP/hResTheta");
108  hbefore->SetName("hResThetaBefore");
109  hbefore->SetTitle("#theta resolution");
110  av_theta_mean_b[isum][idt] = 1e6*(hbefore->GetMean());
111  av_theta_rms_b[isum][idt] = 1e6*(hbefore->GetRMS());
112  avEr_theta_mean_b[isum][idt] = 1e6*(hbefore->GetMeanError());
113  avEr_theta_rms_b[isum][idt] = 1e6*(hbefore->GetRMSError());
114  // events[isum][idt] = 1e-3*(hbefore->GetEntries());
115  // cout<<"ipart = "<<ipart<<" for sum of "<<isum<<" samlpes, with iDa = "<<iDa<<" idt = "<<idt<<" has "<<events[isum][idt]<<" events"<<endl;
116  TFile *f2 = new TFile(name2,"READ");
117  if (f2->IsZombie()) {
118  std::cout << "!!! Error opening file " <<name2<< std::endl;
119  continue;
120  }
121  hafter = (TH1F*)f2->Get("NearIP/hResTheta");
122  hafter->SetName("hResThetaAfter");
123  hafter->SetTitle("#theta resolution");
124  av_theta_mean_a[isum][idt] = 1e6*(hafter->GetMean());
125  av_theta_rms_a[isum][idt] = 1e6*(hafter->GetRMS());
126  avEr_theta_mean_a[isum][idt] = 1e6*(hafter->GetMeanError());
127  avEr_theta_rms_a[isum][idt] = 1e6*(hafter->GetRMSError());
128  TFile *f3 = new TFile(name3,"READ");
129  if (f3->IsZombie()) {
130  std::cout << "!!! Error opening file " <<name2<< std::endl;
131  continue;
132  }
133  hrefer = (TH1F*)f3->Get("NearIP/hResTheta");
134  hrefer->SetName("hResThetaRef");
135  hrefer->SetTitle("#theta resolution");
136  av_theta_mean_r[isum][idt] = 1e6*(hrefer->GetMean());
137  av_theta_rms_r[isum][idt] = 1e6*(hrefer->GetRMS());
138  avEr_theta_mean_r[isum][idt] = 1e6*(hrefer->GetMeanError());
139  avEr_theta_rms_r[isum][idt] = 1e6*(hrefer->GetRMSError());
141  f1->Close();
142  f2->Close();
143  f3->Close();
144 
146  TString namemisc = pathG+"/";
147  namemisc+=tr_sc;
148  namemisc+="mkm_";
149  namemisc+=rt_sc;
150  namemisc+="mrad/";
151  namemisc+=iSc;
152  namemisc+="Sum";
153  // namemisc+=nSd; //sample aligned on max number of tracks [take result from Knossos with the largerst num of rec.trks]
154  namemisc+=isum; //sample for alignment with sum of samples
155  namemisc+="Samples/KnossosResults.root";
156  TFile *fmisc = new TFile(namemisc,"READ");
157  if (fmisc->IsZombie()) {
158  std::cout << "!!! Error opening file " <<namemisc<< std::endl;
159  }
160 
161  TH2F *hmis_b[6];
162  TH2F *hmis_a[6];
163  for(int imdf=0;imdf<6;imdf++){
164  TString nb = "mis_before_";
165  nb+=imdf;
166  hmis_b[imdf] = (TH2F*)fmisc->Get(nb);
167  av_mis_b_mean[imdf][isum][idt] = 1e4*(hmis_b[imdf]->ProjectionY()->GetRMS());
168  avEr_mis_b_mean[imdf][isum][idt] = 1e4*(hmis_b[imdf]->ProjectionY()->GetRMSError());
169 
170  TString na = "mis_diff_";
171  na+=imdf;
172  hmis_a[imdf] = (TH2F*)fmisc->Get(na);
173  av_mis_a_mean[imdf][isum][idt] = 1e4*(hmis_a[imdf]->ProjectionY()->GetRMS());
174  avEr_mis_a_mean[imdf][isum][idt] = 1e4*(hmis_a[imdf]->ProjectionY()->GetRMSError());
175  }
177  fmisc->Close();
178 
179  }//end //loop over dt values
180  }//end loop over diff sum values
181 
183  TCanvas cRes;
184  TGraphErrors *gr_theta_b[nuSd+1];
185  TGraphErrors *gr_theta_a[nuSd+1];
186  TGraphErrors *gr_theta_r[nuSd+1];
187  TGraphErrors *gr_events[nuSd+1];
188  TLegend *legThMb = new TLegend(0.1,0.2,0.38,0.58);
189  TLegend *legThMa = new TLegend(0.1,0.2,0.38,0.58);
190  for(int isum=1;isum<nuSd+1;isum++){//loop over diff stat.samples (plot with diff markers)
191  TString sname = " ";
192  sname += isum;
193  sname +=" samples";
194 
195  gr_events[isum] = new TGraphErrors(nDt,v_tr_sc, av_events[isum],0, avEr_events[isum]);
196  TString name_ev = "number_ev_sum";
197  name_ev += isum;
198  gr_events[isum]->SetName(name_ev);
199  gr_events[isum]->SetMarkerStyle(20+isum);
200  gr_events[isum]->SetMarkerColor(kGreen-3);
201  gr_events[isum]->SetMarkerSize(1.5);
202  mgr_events->Add(gr_events[isum]);
203 
204  gr_theta_b[isum] = new TGraphErrors(nDt,v_tr_sc, av_theta_mean_b[isum],0, avEr_theta_mean_b[isum]);
205  TString name_thm_b = "th_mean_b_sum";
206  name_thm_b += isum;
207  gr_theta_b[isum]->SetName(name_thm_b);
208  gr_theta_b[isum]->SetMarkerStyle(20+isum);
209  gr_theta_b[isum]->SetMarkerColor(kGreen-3);
210  gr_theta_b[isum]->SetMarkerSize(1.5);
211  legThMb->AddEntry(gr_theta_b[isum],sname,"P");
212 
213  gr_theta_a[isum] = new TGraphErrors(nDt,v_tr_sc, av_theta_mean_a[isum],0, avEr_theta_mean_a[isum]);
214  // gr_theta_a[isum]->SetMarkerStyle(30+isum);
215  TString name_thm_a = "th_mean_a_sum";
216  name_thm_a += isum;
217  gr_theta_a[isum]->SetName(name_thm_a);
218  gr_theta_a[isum]->SetMarkerStyle(20+isum);
219  gr_theta_a[isum]->SetMarkerColor(kOrange+7);
220  gr_theta_a[isum]->SetMarkerSize(1.5);
221  //
222  legThMa->AddEntry(gr_theta_a[isum],sname,"P");
223  gr_theta_r[isum] = new TGraphErrors(nDt,v_tr_sc, av_theta_mean_r[isum],0, avEr_theta_mean_r[isum]);
224  // gr_theta_r[isum]->SetMarkerStyle(26+isum);
225  TString name_thm_r = "th_mean_r_sum";
226  name_thm_r += isum;
227  gr_theta_r[isum]->SetName(name_thm_r);
228  gr_theta_r[isum]->SetMarkerStyle(20+isum);
229  gr_theta_r[isum]->SetMarkerColor(15);
230  gr_theta_r[isum]->SetMarkerSize(1.5);
231  mgr_theta_b->Add(gr_theta_b[isum]);
232  mgr_theta->Add(gr_theta_a[isum]);
233  mgr_theta->Add(gr_theta_r[isum]);
234 
235  }
236  mgr_theta->Draw("AP");
237  mgr_theta->GetXaxis()->SetTitle("#delta_{t}, #mum");
238  mgr_theta->GetYaxis()->SetTitle("#theta_{mean}, #murad");
239 
240  legThMa->SetFillColor(0);
241  legThMb->SetFillColor(0);
242  legThMa->Draw();
243  cRes.Print(resname_pdf_o); //write canvas and keep the ps file open
244  cRes.Clear();
245  mgr_theta_b->Draw("AP");
246  mgr_theta_b->GetXaxis()->SetTitle("#delta_{t}, #mum");
247  mgr_theta_b->GetYaxis()->SetTitle("#theta_{mean}, #murad");
248  legThMb->Draw();
249  cRes.Print(resname_pdf_o); //write canvas and keep the ps file open
250  cRes.Clear();
251  TGraphErrors *gr_theta_rms_b[nuSd+1];
252  TGraphErrors *gr_theta_rms_a[nuSd+1];
253  TGraphErrors *gr_theta_rms_r[nuSd+1];
254  for(int isum=1;isum<nuSd+1;isum++){//loop over diff stat.samples (plot with diff markers)
255  gr_theta_rms_b[isum] = new TGraphErrors(nDt,v_tr_sc, av_theta_rms_b[isum],0, avEr_theta_rms_b[isum]);
256  gr_theta_rms_b[isum]->SetMarkerStyle(20+isum);
257  gr_theta_rms_b[isum]->SetMarkerColor(kGreen-3);
258  gr_theta_rms_b[isum]->SetMarkerSize(1.5);
259  TString name_thm_b = "th_rms_b_sum";
260  name_thm_b += isum;
261  gr_theta_rms_b[isum]->SetName(name_thm_b);
262 
263  gr_theta_rms_a[isum] = new TGraphErrors(nDt,v_tr_sc, av_theta_rms_a[isum],0, avEr_theta_rms_a[isum]);
264  // gr_theta_rms_a[isum]->SetMarkerStyle(30+isum);
265  gr_theta_rms_a[isum]->SetMarkerStyle(20+isum);
266  gr_theta_rms_a[isum]->SetMarkerColor(kOrange+7);
267  gr_theta_rms_a[isum]->SetMarkerSize(1.5);
268  TString name_thm_a = "th_rms_a_sum";
269  name_thm_a += isum;
270  gr_theta_rms_a[isum]->SetName(name_thm_a);
271  gr_theta_rms_r[isum] = new TGraphErrors(nDt,v_tr_sc, av_theta_rms_r[isum],0, avEr_theta_rms_r[isum]);
272  // gr_theta_rms_r[isum]->SetMarkerStyle(26+isum);
273  gr_theta_rms_r[isum]->SetMarkerStyle(20+isum);
274  gr_theta_rms_r[isum]->SetMarkerColor(15);
275  gr_theta_rms_r[isum]->SetMarkerSize(1.5);
276  TString name_thm_r = "th_rms_r_sum";
277  name_thm_r += isum;
278  gr_theta_rms_r[isum]->SetName(name_thm_r);
279  mgr_theta_rms_b->Add(gr_theta_rms_b[isum]);
280  mgr_theta_rms->Add(gr_theta_rms_a[isum]);
281  mgr_theta_rms->Add(gr_theta_rms_r[isum]);
282 
283  }
284  mgr_theta_rms->Draw("AP");
285  mgr_theta_rms->GetXaxis()->SetTitle("#delta_{t}, #mum");
286  mgr_theta_rms->GetYaxis()->SetTitle("#theta_{#sigma}, #murad");
287  cRes.Print(resname_pdf_o); //write canvas and keep the ps file open
288  cRes.Clear();
289  mgr_theta_rms_b->Draw("AP");
290  mgr_theta_rms_b->GetXaxis()->SetTitle("#delta_{t}, #mum");
291  mgr_theta_rms_b->GetYaxis()->SetTitle("#theta_{#sigma}, #murad");
292  cRes.Print(resname_pdf_o); //write canvas and keep the ps file open
293  cRes.Clear();
294  mgr_events->Draw("APL");
295  mgr_events->GetXaxis()->SetTitle("#delta_{t}, #mum");
296  mgr_events->GetYaxis()->SetTitle("events, 10^{3}");
297  cRes.Print(resname_pdf_o); //write canvas and keep the ps file open
298  cRes.Clear();
299  //put together missalignment consts results
300  TGraphErrors *gr_mis_b[6][nuSd+1];
301  TGraphErrors *gr_mis_a[6][nuSd+1];
302 
303  for(int ialc=0;ialc<6;ialc++){
304  for(int isum=1;isum<nuSd+1;isum++){//loop over diff stat.samples (plot with diff markers)
305  gr_mis_b[ialc][isum] = new TGraphErrors(nDt,v_tr_sc,av_mis_b_mean[ialc][isum],0,avEr_mis_b_mean[ialc][isum]);
306  gr_mis_b[ialc][isum]->SetMarkerStyle(20+isum);
307  gr_mis_b[ialc][isum]->SetMarkerColor(kGreen-3);
308  gr_mis_b[ialc][isum]->SetMarkerSize(1.5);
309  gr_mis_a[ialc][isum] = new TGraphErrors(nDt,v_tr_sc,av_mis_a_mean[ialc][isum],0,avEr_mis_a_mean[ialc][isum]);
310  gr_mis_a[ialc][isum]->SetMarkerStyle(20+isum);
311  gr_mis_a[ialc][isum]->SetMarkerColor(kOrange+7);
312  gr_mis_a[ialc][isum]->SetMarkerSize(1.5);
313  TString name_miss_a = "mis_a_";
314  name_miss_a +=ialc;
315  name_miss_a +="_sum";
316  name_miss_a +=isum;
317  TString name_miss_b = "mis_b_";
318  name_miss_b +=ialc;
319  name_miss_b +="_sum";
320  name_miss_b +=isum;
321  gr_mis_b[ialc][isum]->SetName(name_miss_b);
322  gr_mis_a[ialc][isum]->SetName(name_miss_a);
323  mgr_mis[ialc]->Add(gr_mis_b[ialc][isum]);
324  mgr_mis[ialc]->Add(gr_mis_a[ialc][isum]);
325  }
326  mgr_mis[ialc]->Draw("AP");
327  mgr_mis[ialc]->GetXaxis()->SetTitle("#delta_{t}, #mum");
328  if(ialc==0) mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{x}, #mum");
329  if(ialc==1) mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{y}, #mum");
330  if(ialc==2) mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{z}, #mum");
331  if(ialc==3) mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{#alpha}, mrad");
332  if(ialc==4) mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{#beta}, mrad");
333  if(ialc==5) mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{#gamma}, mrad");
334  if(ialc<5){
335  if(ialc==0 || ialc==1)
336  cRes.Print(resname_pdf_o); //write canvas and keep the file open
337  }
338  else cRes.Print(resname_pdf_c); //write canvas and close pdf file
339  cRes.Clear();
340  }
341 
342 
343 
345  TString out = resname+".root";
346  TFile *f = new TFile(out,"RECREATE");
347  mgr_theta->Write();
348  mgr_theta_rms->Write();
349  mgr_theta_b->Write();
350  mgr_theta_rms_b->Write();
351  for(int ialc=0;ialc<6;ialc++){
352  mgr_mis[ialc]->Write();
353  }
354 
355  for(int isum=1;isum<nuSd+1;isum++){//loop over diff stat.samples (plot with diff markers)
356  gr_events[isum]->Write();
357  gr_theta_b[isum]->Write();
358  gr_theta_a[isum]->Write();
359  gr_theta_r[isum]->Write();
360  gr_theta_rms_b[isum]->Write();
361  gr_theta_rms_a[isum]->Write();
362  gr_theta_rms_r[isum]->Write();
363  for(int ialc=0;ialc<6;ialc++){
364  gr_mis_b[ialc][isum]->Write();
365  gr_mis_a[ialc][isum]->Write();
366  }
367  }
368 
369  f->Write();
370  f->Close();
371  }
372  return 0;
373 }
374 
375 int main(int __argc,char *__argv[]) {
376  TString storePath="/home/karavdin/datastorage/AlignmentLMDpixel/100data_10consts_SamplesStudyOct2013/BOX/mom_15/";
377  double daGL=0;
378  int iSc=0;
379  std::string pathStr="",daStr="",iScStr="";
380  // decode arguments
381  if( __argc>1 && ( strcmp( __argv[1], "-help" ) == 0
382  || strcmp( __argv[1], "--help" ) == 0 ) ){
383  std::cout << "This is script for alignment data analysis \n"
384  << "with parameters\n"
385  <<"-path path to the file(s) \n"
386  <<"-da scale for rotation misalignment \n"
387  <<"-iSc number of misalignment consts sample \n"
388  <<"Have fun! \n"
389  << std::endl;
390  return 0;
391  }
392  while ((optind < (__argc-1) ) && (__argv[optind][0]=='-')) {
393  bool found=false;
394  std::string sw = __argv[optind];
395  if (sw=="-path"){
396  optind++;
397  pathStr = __argv[optind];
398  found=true;
399  }
400  if (sw=="-da"){
401  optind++;
402  daStr = __argv[optind];
403  found=true;
404  }
405 
406  if (sw=="-iSc"){
407  optind++;
408  iScStr = __argv[optind];
409  found=true;
410  }
411  if (!found){
412  std::cout<< "Unknown switch: "
413  << __argv[optind] <<std::endl;
414  optind++;
415  }
416  while ( (optind < __argc ) && __argv[optind][0]!='-' ) optind++;
417  }
418 
419  std::stringstream pathSStr(pathStr), daSStr(daStr),iScSStr(iScStr);
420  pathSStr >> storePath;
421  daSStr >> daGL;
422  // isumSStr >> isum;
423  iScSStr >> iSc;
424  LimitsAlign(storePath,daGL,iSc);
425  return 0;
426 }
name3
TString name3
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:200
main
int main(int argc, char **argv)
Definition: MvdMQ/run/runFileSampler.cxx:29
name2
TString name2
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:200
f1
TF1 * f1
Definition: reco_analys2.C:50
storePath
TString storePath
Definition: runGeaneTestSimBox.cxx:50
f3
TFile * f3
Definition: ReadMesDecCorrRes.C:8
name1
TString name1
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:200
f
TFile * f
Definition: bump_analys.C:12
out
TFile * out
Definition: reco_muo.C:20
f2
TFile * f2
Definition: ReadMesDecCorrRes.C:5
LimitsAlign
int LimitsAlign(TString pathG="/home/karavdin/datastorage/AlignmentLMDpixel/10data_20consts_SamplesStudyOct2013/BOX/mom_15/", double daGL=0, int iSc=0)
Definition: LimitsAlign.C:24
Generated on Wed Apr 3 2019 10:02:47 for FairRoot/PandaRoot by   doxygen 1.8.5