FairRoot/PandaRoot
ModulesStatFast.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 using namespace std;
16 int ModulesStatFast(TString pathG="/home/karavdin/datastorage/AlignmentLMDpixel/100data_10consts_SamplesStudyOct2013/BOX/mom_15/", double daGL=0, int MisConstSamp=0)
17 {
18  TStopwatch timerB;
19  timerB.Start();
20  int iSc = MisConstSamp;
21  const int nDt = 8;
22  double v_tr_sc[nDt]={0,50,100,200,300,400,500,600};//misalign, trans
23  const int nMod = 10;
24  const int nDa = 1; //TEST
25  double v_rt_sc[nDa]={daGL};//misalign, trans
26 
27  const int nSdGL=100;//number of data samples
28  const int nSdstep = 10; //step size for "number of used trks in Knossos"
29  const int nuSd = nSdGL/nSdstep;
30  const int nDS2 = 10;//number of samples simulated in stepC for the same case
31  // cout<<"nuSd = "<<nuSd<<endl;
32  // const int nSc=10;//number of align.consts samples
33  // const int nSc=1;//number of align.consts samples TEST
34 
35  // for(int iSc=0;iSc<nSc;iSc++){//diff output files
36  for(int iDa=0;iDa<nDa;iDa++){//diff output files
37 
38  double av_evMod_b[nMod][nuSd+1][nDt];
39  double av_evMod_a[nMod][nuSd+1][nDt];
40  double err_av_evMod_b[nMod][nuSd+1][nDt];
41  double err_av_evMod_a[nMod][nuSd+1][nDt];
42  double av_resxMod_b[nMod][nuSd+1][nDt];
43  double av_resxMod_a[nMod][nuSd+1][nDt];
44  double err_av_resxMod_b[nMod][nuSd+1][nDt];
45  double err_av_resxMod_a[nMod][nuSd+1][nDt];
46  double av_resyMod_b[nMod][nuSd+1][nDt];
47  double av_resyMod_a[nMod][nuSd+1][nDt];
48  double err_av_resyMod_b[nMod][nuSd+1][nDt];
49  double err_av_resyMod_a[nMod][nuSd+1][nDt];
50  //How to save data
51  TString resname= pathG+"/AlignmentLimits_EvModStat_";
52  resname+=v_rt_sc[iDa];
53  resname+="mrad_";
54  resname+=iSc;
55  resname+="constsSample";
56  TString resname_pdf = resname+".pdf";
57  TString resname_pdf_o = resname_pdf+"(";
58  TString resname_pdf_c = resname_pdf+")";
59 
60  for(int idt=0;idt<nDt;idt++){//loop over dt values (x-axis of plots)
61  for(int imod=0;imod<nMod;imod++){// loop over modules (diff plots)
62  for(int ijsum=1;ijsum<nuSd+1;ijsum++){ //loop over diff stat.samples (plot with diff markers)
63  int isum = 1+nSdstep*(ijsum-1);
64  double tr_sc = v_tr_sc[idt];
65  double rt_sc = v_rt_sc[iDa];
66  //Where search for files
67  TString path = pathG+"/";
68  path+=tr_sc;
69  path+="mkm_";
70  path+=rt_sc;
71  path+="mrad/";
72  path+=iSc;
73  path+="Sum";
74  path+=isum; //sample for alignment with sum of samples
75  path+="Samples/resultsSUM/";
76 
78  TString name1 = path + "/hitsRes_before_alignment_sector";
79  name1+=imod;
80  name1+=".root";
81  TString name2 = path + "/hitsRes_after_alignment_sector";
82  name2+=imod;
83  name2+=".root";
84  TFile *f1 = new TFile(name1,"READ");
85  TNtuple *nb = (TNtuple*)f1->Get("nhits");
86  // nb->Print();
87  TH1I *ntrksb = (TH1I*)f1->Get("trks");
88 
89 
90  av_resxMod_b[imod][ijsum][idt]=0;
91  av_resxMod_a[imod][ijsum][idt]=0;
92  err_av_resxMod_b[imod][ijsum][idt]=0;
93  err_av_resxMod_a[imod][ijsum][idt]=0;
94  av_resyMod_b[imod][ijsum][idt]=0;
95  av_resyMod_a[imod][ijsum][idt]=0;
96  err_av_resyMod_b[imod][ijsum][idt]=0;
97  err_av_resyMod_a[imod][ijsum][idt]=0;
98 
99  av_evMod_b[imod][ijsum][idt]= 1e-3*double(ntrksb->GetEntries())/nDS2;
100  err_av_evMod_b[imod][ijsum][idt]= 0;
101  TH1D *hresx_b = new TH1D("hresx_b","",1e3,-1.,1.);
102  nb->Project("hresx_b","xrec-xmc");
103  av_resxMod_b[imod][ijsum][idt]=1e4*(hresx_b->GetMean());
104  err_av_resxMod_b[imod][ijsum][idt]=1e4*(hresx_b->GetRMS());
105  TH1D *hresy_b = new TH1D("hresy_b","",1e3,-1.,1.);
106  nb->Project("hresy_b","yrec-ymc");
107  av_resyMod_b[imod][ijsum][idt]=1e4*(hresy_b->GetMean());
108  err_av_resyMod_b[imod][ijsum][idt]=1e4*(hresy_b->GetRMS());
109  TFile *f2 = new TFile(name2,"READ");
110  TNtuple *na = (TNtuple*)f2->Get("nhits");
111  // ev_mod_a[ipart]= 1e-3*double(na->GetEntries());
112  TH1I *ntrksa = (TH1I*)f2->Get("trks");
113  av_evMod_a[imod][ijsum][idt]= 1e-3*double(ntrksa->GetEntries())/nDS2;
114  err_av_evMod_a[imod][ijsum][idt]= 0;
115  TH1D *hresx_a = new TH1D("hresx_a","",1e3,-1.,1.);
116  na->Project("hresx_a","xrec-xmc");
117  av_resxMod_a[imod][ijsum][idt]=1e4*(hresx_a->GetMean());
118  err_av_resxMod_a[imod][ijsum][idt]=1e4*(hresx_a->GetRMS());
119  TH1D *hresy_a = new TH1D("hresy_a","",1e3,-1.,1.);
120  na->Project("hresy_a","yrec-ymc");
121  // resy_mod_a[ipart]= 1e4*(hresy_a->GetMean());
122  av_resyMod_a[imod][ijsum][idt]=1e4*(hresy_a->GetMean());
123  err_av_resyMod_a[imod][ijsum][idt]=1e4*(hresy_a->GetRMS());
124  f1->Close();
125  f2->Close();
126 
127  cout<<"iS_AC:"<<iSc<<" ida:"<<iDa<<" dt:"<<idt<<" imod:"<<imod<<" isum:"<<isum<<endl;
128  // timerB.Stop();
129  // Double_t rtime = timerB.RealTime();
130  // Double_t ctime = timerB.CpuTime();
131  // printf("for the LOOP: RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
132  }//end loop over diff sum values
133  }// end loop over modules (diff plots)
134  }//end //loop over dt values
135 
136 
137  TString out = resname+".root";
138  TFile *f = new TFile(out,"RECREATE");
139 
141  TCanvas cRes;
142  for(int imod=0;imod<nMod;imod++){// loop over modules (diff plots)
143  TGraphErrors *gr_av_b[nuSd+1];
144  TGraphErrors *gr_av_a[nuSd+1];
145  TMultiGraph *mgr_av_b = new TMultiGraph();
146  TMultiGraph *mgr_av_a = new TMultiGraph();
147  TGraphErrors *gr_avResx_b[nuSd+1];
148  TGraphErrors *gr_avResx_a[nuSd+1];
149  TMultiGraph *mgr_avResx_b = new TMultiGraph();
150  TMultiGraph *mgr_avResx_a = new TMultiGraph();
151  TGraphErrors *gr_avResy_b[nuSd+1];
152  TGraphErrors *gr_avResy_a[nuSd+1];
153  TMultiGraph *mgr_avResy_b = new TMultiGraph();
154  TMultiGraph *mgr_avResy_a = new TMultiGraph();
155  TLegend *leg_b = new TLegend(0.91,0.65,0.99,0.98);
156  leg_b->SetHeader("BEFORE");
157  TLegend *leg_a = new TLegend(0.91,0.65,0.99,0.98);
158  leg_a->SetHeader("AFTER");
159  // for(int isum=nSd;isum>0;isum--){//loop over diff stat.samples (plot with diff markers)
160  // for(int isum=nSd;isum>0;isum=isum-10){//loop over diff stat.samples (plot with diff markers)
161  for(int ijsum=nuSd;ijsum>0;ijsum--){ //loop over diff stat.samples (plot with diff markers)
162  // int isum = nSdstep*ijsum;
163  int isum = 1+nSdstep*(ijsum-1);
164  TString sname = " ";
165  sname +=" sum of ";
166  sname += isum;
167  gr_av_b[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_evMod_b[imod][ijsum],0,err_av_evMod_b[imod][ijsum]);
168  gr_av_a[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_evMod_a[imod][ijsum],0,err_av_evMod_a[imod][ijsum]);
169  gr_avResx_b[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_resxMod_b[imod][ijsum],0,err_av_resxMod_b[imod][ijsum]);
170  gr_avResx_a[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_resxMod_a[imod][ijsum],0,err_av_resxMod_a[imod][ijsum]);
171  gr_avResx_b[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_resxMod_b[imod][ijsum],0,err_av_resxMod_b[imod][ijsum]);
172  gr_avResx_a[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_resxMod_a[imod][ijsum],0,err_av_resxMod_a[imod][ijsum]);
173  gr_avResy_b[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_resyMod_b[imod][ijsum],0,err_av_resyMod_b[imod][ijsum]);
174  gr_avResy_a[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_resyMod_a[imod][ijsum],0,err_av_resyMod_a[imod][ijsum]);
175  gr_avResy_b[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_resyMod_b[imod][ijsum],0,err_av_resyMod_b[imod][ijsum]);
176  gr_avResy_a[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_resyMod_a[imod][ijsum],0,err_av_resyMod_a[imod][ijsum]);
177  TString name_ev = "Ev_";
178  name_ev +=imod;
179  name_ev +="mod_sum";
180  name_ev += isum;
181  gr_av_b[ijsum]->SetName(name_ev+"_b");
182  gr_av_b[ijsum]->SetMarkerStyle(20+ijsum);
183  gr_av_b[ijsum]->SetMarkerColor(kGreen-3);
184  gr_av_b[ijsum]->SetMarkerSize(1.5);
185  gr_av_a[ijsum]->SetName(name_ev+"_a");
186  gr_av_a[ijsum]->SetMarkerStyle(20+ijsum);
187  gr_av_a[ijsum]->SetMarkerColor(kOrange+7);
188  gr_av_a[ijsum]->SetMarkerSize(1.5);
189  mgr_av_b->Add(gr_av_b[ijsum]);
190  mgr_av_a->Add(gr_av_a[ijsum]);
191  gr_av_b[ijsum]->Write();
192  gr_av_a[ijsum]->Write();
193  leg_b->AddEntry(gr_av_b[ijsum],sname,"P");
194  leg_a->AddEntry(gr_av_a[ijsum],sname,"P");
195 
196  gr_avResx_b[ijsum]->SetName(name_ev+"_b");
197  gr_avResx_b[ijsum]->SetMarkerStyle(20+ijsum);
198  gr_avResx_b[ijsum]->SetMarkerColor(kGreen-3);
199  gr_avResx_b[ijsum]->SetMarkerSize(1.5);
200  gr_avResx_a[ijsum]->SetName(name_ev+"_a");
201  gr_avResx_a[ijsum]->SetMarkerStyle(20+ijsum);
202  gr_avResx_a[ijsum]->SetMarkerColor(kOrange+7);
203  gr_avResx_a[ijsum]->SetMarkerSize(1.5);
204  mgr_avResx_b->Add(gr_avResx_b[ijsum]);
205  mgr_avResx_a->Add(gr_avResx_a[ijsum]);
206  gr_avResx_b[ijsum]->Write();
207  gr_avResx_a[ijsum]->Write();
208  gr_avResy_b[ijsum]->SetName(name_ev+"_b");
209  gr_avResy_b[ijsum]->SetMarkerStyle(20+ijsum);
210  gr_avResy_b[ijsum]->SetMarkerColor(kGreen-3);
211  gr_avResy_b[ijsum]->SetMarkerSize(1.5);
212  gr_avResy_a[ijsum]->SetName(name_ev+"_a");
213  gr_avResy_a[ijsum]->SetMarkerStyle(20+ijsum);
214  gr_avResy_a[ijsum]->SetMarkerColor(kOrange+7);
215  gr_avResy_a[ijsum]->SetMarkerSize(1.5);
216  mgr_avResy_b->Add(gr_avResy_b[ijsum]);
217  mgr_avResy_a->Add(gr_avResy_a[ijsum]);
218  gr_avResy_b[ijsum]->Write();
219  gr_avResy_a[ijsum]->Write();
220 
221  }//end over diff stat.samples (plot with diff markers)
222 
223  TString smname = "Module_";
224  smname += imod;
225  mgr_av_b->SetTitle(smname+"_b");
226  mgr_av_a->SetTitle(smname+"_a");
227  leg_a->SetFillColor(0);
228  leg_b->SetFillColor(0);
229  mgr_av_b->SetName(smname+"_b");
230  mgr_av_a->SetName(smname+"_a");
231  TString smnamex,smnamey;
232  smnamex =smname;
233  smnamex+="_resX";
234  smnamey =smname;
235  smnamey+="_resY";
236  mgr_avResx_b->SetName(smnamex+"_b");
237  mgr_avResx_a->SetName(smnamex+"_a");
238  mgr_avResy_b->SetName(smnamey+"_b");
239  mgr_avResy_a->SetName(smnamey+"_a");
240  mgr_av_b->Draw("AP");
241  mgr_av_b->GetXaxis()->SetTitle("#delta_{t}, #mum");
242  mgr_av_b->GetYaxis()->SetTitle("events, 10^{3}");
243  leg_b->Draw();
244  cRes.Print(resname_pdf_o); //write canvas and keep the ps file open
245  cRes.Clear();
246  mgr_av_a->Draw("AP");
247  mgr_av_a->GetXaxis()->SetTitle("#delta_{t}, #mum");
248  mgr_av_a->GetYaxis()->SetTitle("events, 10^{3}");
249  leg_a->Draw();
250  cRes.Print(resname_pdf_o); //write canvas and keep the ps file open
251  cRes.Clear();
252  mgr_av_a->Write();
253  mgr_av_b->Write();
254  if(imod<1){
255  leg_b->SetName("LegBefore");
256  leg_a->SetName("LegAfter");
257  leg_b->Write();
258  leg_a->Write();
259  }
260 
261  mgr_avResx_b->Draw("AP");
262  mgr_avResx_b->GetXaxis()->SetTitle("#delta_{t}, #mum");
263  mgr_avResx_b->GetYaxis()->SetTitle("res_{X}, mean [#mum]");
264  leg_b->Draw();
265  cRes.Print(resname_pdf_o); //write canvas and keep the ps file open
266  cRes.Clear();
267  mgr_avResx_a->Draw("AP");
268  mgr_avResx_a->GetXaxis()->SetTitle("#delta_{t}, #mum");
269  mgr_avResx_a->GetYaxis()->SetTitle("res_{X}, mean [#mum]");
270  leg_a->Draw();
271  cRes.Print(resname_pdf_o); //write canvas and keep the ps file open
272  cRes.Clear();
273  mgr_avResx_a->Write();
274  mgr_avResx_b->Write();
275 
276  mgr_avResy_b->Draw("AP");
277  mgr_avResy_b->GetXaxis()->SetTitle("#delta_{t}, #mum");
278  mgr_avResy_b->GetYaxis()->SetTitle("res_{Y}, mean [#mum]");
279  leg_b->Draw();
280  cRes.Print(resname_pdf_o); //write canvas and keep the ps file open
281  cRes.Clear();
282  mgr_avResy_a->Draw("AP");
283  mgr_avResy_a->GetXaxis()->SetTitle("#delta_{t}, #mum");
284  mgr_avResy_a->GetYaxis()->SetTitle("res_{Y}, mean [#mum]");
285  leg_a->Draw();
286  cRes.Print(resname_pdf_o); //write canvas and keep the ps file open
287  cRes.Clear();
288  mgr_avResy_a->Write();
289  mgr_avResy_b->Write();
290 
291  }//end loop over modules (diff plots)
292  cRes.Print(resname_pdf_c); //write canvas and close the ps file
294  f->Write();
295  f->Close();
296  }
297  timerB.Stop();
298  Double_t rtime = timerB.RealTime();
299  Double_t ctime = timerB.CpuTime();
300  printf("On HIMSTER it took: RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
301  return 0;
302 }
303 
304 int main(int __argc,char *__argv[]) {
305  TString storePath="/home/karavdin/datastorage/AlignmentLMDpixel/100data_10consts_SamplesStudyOct2013/BOX/mom_15/";
306  double daGL=0;
307  int iSc=0;
308  std::string pathStr="",daStr="",iScStr="";
309  // decode arguments
310  if( __argc>1 && ( strcmp( __argv[1], "-help" ) == 0
311  || strcmp( __argv[1], "--help" ) == 0 ) ){
312  std::cout << "This is script for alignment data analysis \n"
313  << "with parameters\n"
314  <<"-path path to the file(s) \n"
315  <<"-da scale for rotation misalignment \n"
316  <<"-iSc number of misalignment consts sample \n"
317  <<"Have fun! \n"
318  << std::endl;
319  return 0;
320  }
321  while ((optind < (__argc-1) ) && (__argv[optind][0]=='-')) {
322  bool found=false;
323  std::string sw = __argv[optind];
324  if (sw=="-path"){
325  optind++;
326  pathStr = __argv[optind];
327  found=true;
328  }
329  if (sw=="-da"){
330  optind++;
331  daStr = __argv[optind];
332  found=true;
333  }
334 
335  if (sw=="-iSc"){
336  optind++;
337  iScStr = __argv[optind];
338  found=true;
339  }
340  if (!found){
341  std::cout<< "Unknown switch: "
342  << __argv[optind] <<std::endl;
343  optind++;
344  }
345  while ( (optind < __argc ) && __argv[optind][0]!='-' ) optind++;
346  }
347 
348  std::stringstream pathSStr(pathStr), daSStr(daStr),iScSStr(iScStr);
349  pathSStr >> storePath;
350  daSStr >> daGL;
351  // isumSStr >> isum;
352  iScSStr >> iSc;
353  ModulesStatFast(storePath,daGL,iSc);
354  return 0;
355 }
main
int main(int argc, char **argv)
Definition: MvdMQ/run/runFileSampler.cxx:29
printf
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
name2
TString name2
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:200
f1
TF1 * f1
Definition: reco_analys2.C:50
storePath
TString storePath
Definition: runGeaneTestSimBox.cxx:50
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
name1
TString name1
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:200
f
TFile * f
Definition: bump_analys.C:12
ModulesStatFast
int ModulesStatFast(TString pathG="/home/karavdin/datastorage/AlignmentLMDpixel/100data_10consts_SamplesStudyOct2013/BOX/mom_15/", double daGL=0, int MisConstSamp=0)
Definition: ModulesStatFast.C:16
out
TFile * out
Definition: reco_muo.C:20
ctime
Double_t ctime
Definition: hit_dirc.C:114
f2
TFile * f2
Definition: ReadMesDecCorrRes.C:5
rtime
Double_t rtime
Definition: hit_dirc.C:113
Generated on Wed Apr 3 2019 10:02:47 for FairRoot/PandaRoot by   doxygen 1.8.5