FairRoot/PandaRoot
online_monitoring_studies.C
Go to the documentation of this file.
1 #include <iostream>
2 
3 #include<TCanvas.h>
4 #include<TApplication.h>
5 #include<TROOT.h>
6 #include<TChain.h>
7 #include<TStyle.h>
8 #include<TGaxis.h>
9 #include<TColor.h>
10 #include<TClonesArray.h>
11 #include<TH1.h>
12 #include<TH2.h>
13 #include<vector>
14 #include<TTree.h>
15 #include<sstream>
16 
17 #include<PndLmdDim.h>
18 #include<PndMCTrack.h>
19 #include<PndSdsMCPoint.h>
20 #include<TGraphErrors.h>
21 
22 using namespace std;
23 
24 // change some stuff in the ROOT appearance
25 void Root_Appearance();
26 
27 double last_percent(-1.);
28 // draw a progress bar only when the length changes significantly
29 void DrawProgressBar(int len, double percent);
30 
31 // the application
32 // Assuming to find an mc file called
33 // ./Lumi_MC_0.root
34 int online_monitoring_studies();
35 
36 int main() {
37  online_monitoring_studies();
38  return 0;
39 }
40 
41 int online_monitoring_studies(){
42  Root_Appearance();
43  TApplication myapp("myapp", 0, 0);
44  cout << " lmd online monitoring studies " << endl;
45 
46  // get the lmd dim help class
47  PndLmdDim& lmddim = PndLmdDim::Get_instance();
48 
49 
50  const int nmomenta = 1;//3;
51  double momenta[nmomenta] = {1.5};//{1.5, 4.06, 15.00};
52  double pixelsize = 80; // um
53  double ip_size = 1.; // mm
54  double ip_div = 1.; // mrad
55  double ip_width_z = 5.; // mm
56  // z positions of planes in m
57  vector < double > z; z.push_back(11.2); z.push_back(11.4); z.push_back(11.5); z.push_back(11.6);
58  // acceptance in the plane (already the power of 2)
59  const double Rmin = 0.035*0.035;
60  const double Rmax = (0.035+0.05)*(0.035+0.05);
61  // cut for a cone starting from ip
62  // angles for a cone cut
63  const double theta_cut_low = 4.; // mrad
64  const double theta_cut_high = 8.; // mrad
65  vector < double > Rmin_cuts;
66  vector < double > Rmax_cuts;
67  for (int iplane = 0; iplane < 4; iplane++ ){
68  Rmin_cuts.push_back(theta_cut_low*1e-3*z[iplane]*theta_cut_low*1e-3*z[iplane]);
69  Rmax_cuts.push_back(theta_cut_high*1e-3*z[iplane]*theta_cut_high*1e-3*z[iplane]);
70  }
71 
72  TCanvas canvas("canvas", "canvas", 600, 1000);
73  canvas.Divide(2,1);
74  canvas.cd(1);
75  TCanvas canvas_distributions("canvas_distributions", "distributions", 800, 800);
76  canvas_distributions.Divide(2,2);
77 
78  //TFile treefile("asymmetries_results.root", "RECREATE");
79  TTree results("asymmetry_results", "Asymmetry results");
80  double offsetx, offsety, tiltx, tilty, momentum,
81  offsetx_lmd, offsety_lmd, tiltx_lmd, tilty_lmd,
82  offsetx_fit, offsety_fit, tiltx_fit, tilty_fit,
83  meanx_0, meany_0, rmsx_0, rmsy_0,
84  meanx_1, meany_1, rmsx_1, rmsy_1,
85  meanx_2, meany_2, rmsx_2, rmsy_2,
86  meanx_3, meany_3, rmsx_3, rmsy_3;
87  int count_0[10] = {0,0,0,0,0,0,0,0,0,0};
88  int count_1[10] = {0,0,0,0,0,0,0,0,0,0};
89  int count_2[10] = {0,0,0,0,0,0,0,0,0,0};
90  int count_3[10] = {0,0,0,0,0,0,0,0,0,0};
91  TLine line;
92 
93  results.Branch("offsetx", &offsetx);
94  results.Branch("offsety", &offsety);
95  results.Branch("tiltx", &tiltx);
96  results.Branch("tilty", &tilty);
97  results.Branch("offsetx_lmd", &offsetx_lmd);
98  results.Branch("offsety_lmd", &offsety_lmd);
99  results.Branch("tiltx_lmd", &tiltx_lmd);
100  results.Branch("tilty_lmd", &tilty_lmd);
101  results.Branch("momentum", &momentum);
102  results.Branch("offsetx_fit", &offsetx_fit);
103  results.Branch("offsety_fit", &offsety_fit);
104  results.Branch("tiltx_fit", &tiltx_fit);
105  results.Branch("tilty_fit", &tilty_fit);
106  results.Branch("meanx_0", &meanx_0);
107  results.Branch("meany_0", &meany_0);
108  results.Branch("rmsx_0", &rmsx_0);
109  results.Branch("rmsy_0", &rmsy_0);
110  results.Branch("meanx_1", &meanx_1);
111  results.Branch("meany_1", &meany_1);
112  results.Branch("rmsx_1", &rmsx_1);
113  results.Branch("rmsy_1", &rmsy_1);
114  results.Branch("meanx_2", &meanx_2);
115  results.Branch("meany_2", &meany_2);
116  results.Branch("rmsx_2", &rmsx_2);
117  results.Branch("rmsy_2", &rmsy_2);
118  results.Branch("meanx_3", &meanx_3);
119  results.Branch("meany_3", &meany_3);
120  results.Branch("rmsx_3", &rmsx_3);
121  results.Branch("rmsy_3", &rmsy_3);
122  results.Branch("count_0", count_0, "count_0[10]/I");
123  results.Branch("count_1", count_1, "count_1[10]/I");
124  results.Branch("count_2", count_2, "count_2[10]/I");
125  results.Branch("count_3", count_3, "count_3[10]/I");
126 
127 
128  vector< string > files;
129 
130  // 15 GeV/c
131  double pbeam = 1.5;
132  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.1_0.1_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
133  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_-0.5_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
134  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.001_0.001_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
135  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_-0.0002_-0.0002_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
136  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_-0.0002_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
137  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_-0.0002_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
138  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_-0.0001_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
139  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
140  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_-0.0001_0.0001_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
141  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0002_-0.0002_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
142  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0001_-0.0001_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
143  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0001_0.0001_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
144  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_-0.0001_-0.0001_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
145  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_-0.0002_0.0002_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
146  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.001_0.001/200000/mc_data/asymmetries_results.root");
147  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.001_0.001/200000/mc_data/Lumi_MC_200000.root");
148  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_-0.0001_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
149  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0001_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
150  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0001_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
151  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0002_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
152  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0002_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
153  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.5_0.5_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
154  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.1_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
155  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.5_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
156  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.1_-0.1_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/asymmetries_results.root");
157  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.1_-0.1_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
158  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_-0.1_-0.1_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
159  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_-0.1_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
160  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.5_-0.5_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
161  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_-0.5_0.5_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
162  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_-0.5_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
163  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_-0.1_0.1_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
164  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_-0.1_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
165  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.1_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
166  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_-0.5_-0.5_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
167  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.5_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data/Lumi_MC_200000.root");
168 /*
169 
170  // 1.5 GeV/c
171  double pbeam = 1.5;
172  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.1_0.1_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data");
173  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_-0.5_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data");
174  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.001_0.001_0.0_0.0/200000/mc_data");
175  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data");
176  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.001_0.001/200000/mc_data");
177  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.5_0.5_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data");
178  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.1_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data");
179  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.5_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data");
180  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.1_-0.1_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data");
181  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_-0.1_-0.1_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data");
182  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_-0.1_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data");
183  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.5_-0.5_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data");
184  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_-0.5_0.5_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data");
185  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_-0.5_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data");
186  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_-0.1_0.1_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data");
187  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_-0.1_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data");
188  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.0_0.1_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data");
189  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_-0.5_-0.5_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data");
190  files.push_back("./dpm_elastic_thmin_0.06deg/ip_offset_XYZDXDYDZ_0.5_0.0_0.0_0.08_0.08_0.35/beam_grad_XYDXDY_0.0_0.0_0.0_0.0/200000/mc_data");
191 */
192  int nfiles = files.size();
193  for (int ifile = 0; ifile < nfiles; ifile++){
194  cout << endl << " file " << ifile << " of " << nfiles << " = " << files[ifile] << endl;
195  // load the MC File
196  TChain tMC("pndsim");
197  tMC.Add((files[ifile]+"/Lumi_MC*.root").c_str());
198 
199  //--- assign MC info -----------------------------------------------------
200  TClonesArray* true_tracks = new TClonesArray("PndMCTrack");
201  tMC.SetBranchAddress("MCTrack", &true_tracks); //True Track to compare
202 
203  TClonesArray* true_points = new TClonesArray("PndSdsMCPoint");
204  tMC.SetBranchAddress("LMDPoint", &true_points); //True Points to compare
205  stringstream filename;
206  filename << "asymmetry_studies_"<< momentum << "_GeVc_x" << offsetx << "_mm_y_" << offsety << "_mm_dxdz_" << tiltx << "_mrad_dydz_" << tilty << "_mrad";
207  //TFile file((filename.str()+".root").c_str(), "RECREATE");
208 
209  TH2F hist_gen_theta_phi ("hist_gen_theta_phi", "theta phi generated distribution", 100, 2e-3, 9e-3, 100, -3.141, 3.141);
210  TH2F hist_gen_x_y ("hist_gen_x_y", "x y generated distribution", 100, -5e-3, 5e-3, 100, -5e-3, 5e-3);
211  TH2F hist_gen_dx_dy ("hist_gen_dx_dy", "dx/dz dy/dz generated distribution", 100, -15e-3, 15e-3, 100, -15e-3, 15e-3);
212  TH2F hist_XY_plane_0 ("hist_XY_plane_0", "xy distribution at plane 0", 100, -0.100, 0.100, 100, -0.100, 0.100);
213  TH2F hist_XY_plane_1 ("hist_XY_plane_1", "xy distribution at plane 1", 100, -0.100, 0.100, 100, -0.100, 0.100);
214  TH2F hist_XY_plane_2 ("hist_XY_plane_2", "xy distribution at plane 2", 100, -0.100, 0.100, 100, -0.100, 0.100);
215  TH2F hist_XY_plane_3 ("hist_XY_plane_3", "xy distribution at plane 3", 100, -0.100, 0.100, 100, -0.100, 0.100);
216  TH2F* hist_XY_plane[4] = {&hist_XY_plane_0,&hist_XY_plane_1,&hist_XY_plane_2,&hist_XY_plane_3};
217  //TH2F hist_phi_count("hist_phi_count", "counts over plane number", 4, -0.5, 3.5, 10, -3.141, 3.141);
218  //TH2F hist_center_of_gravity("hist_center_of_gravity", "center of gravity over plane number", 4, -0.5, 3.5, 10, -3.141, 3.141);
219 
220  // calculate some multiple scattering constants
221 
222 // double beta = 1/sqrt(1+pow(0.938/momentum,2));
223 // cout << "the beta factor for " << momentum << " GeV/c antiprotons is " << beta << endl;
224 // double x_over_x0 = 2.329/21.82*350/1000/10;
225 // cout << "x/X0 is" << x_over_x0 << endl;
226 // double theta_0_slicon = 13.6/beta/(momentum*1e3)*sqrt(x_over_x0)*(1+0.038*log(x_over_x0));
227 // cout << "scattering angle of 350 mum silicon is " << theta_0_slicon << endl;
228  // assuming a multiple scattering effect in the transition region of a silicon plane of half thickness
229 // double theta_0_kapton = 13.6/beta/(momentum*1e3)*sqrt(x_over_x0/2.)*(1+0.038*log(x_over_x0/2.));
230 
231  int nEvents = tMC.GetEntries();
232  int nEvents_counted(0);
233  int nhitstotal(0);
234 
235  tiltx = 0.;
236  tilty = 0.;
237  offsetx = 0.;
238  offsety = 0.;
239  momentum = 0.;
240 
241  int npbar = 0; // to normalize the mean for the above values to
242 
243  cout << " reading " << nEvents << " Events " << endl;
244  for (Int_t j = 0; j < nEvents ; j++) {
245  nEvents_counted++;
246  DrawProgressBar(50, (j + 1) / ((double) nEvents));
247  tMC.GetEntry(j);
248 
249  int nHits = true_points->GetEntriesFast();
250  int nSdsHits = 0;
251 
252  const int nParticles = true_tracks->GetEntriesFast();
253 
254  // get the beam properties of the IP
255 
256  for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
257  PndMCTrack *mctrk = (PndMCTrack*) true_tracks->At(iParticle);
258  Int_t mcID = mctrk->GetPdgCode();
259  if (mctrk->IsGeneratorCreated() && mcID == -2212) {
260  npbar++;
261  TVector3 momMC = mctrk->GetMomentum();
262  tiltx += momMC.X()/momMC.Z();
263  tilty += momMC.Y()/momMC.Z();
264  momentum += momMC.Mag();
265  TVector3 posMC = mctrk->GetStartVertex();
266  offsetx += posMC.X();
267  offsety += posMC.Y();
268  hist_gen_theta_phi.Fill(momMC.Theta()*1e3, momMC.Phi());
269  hist_gen_x_y.Fill(posMC.X()*10., posMC.X()*10.);
270  hist_gen_dx_dy.Fill(momMC.X()/momMC.Z()*1e3, momMC.Y()/momMC.Z()*1e3);
271  }
272  }
273 
274  for (Int_t iHit = 0; iHit < nHits; iHit++) {
275  PndSdsMCPoint* mcpoint = (PndSdsMCPoint*) true_points->At(iHit);
276  double theta_ip = 0.;
277  //if (mcpoint->GetTrackID() >= 0) {
278  // (true_tracks->At(mcpoint->GetTrackID))->;
279  //}
280  // cout << mcpoint->GetEnergyLoss() << endl;
281  //hist_eloss_total->Fill(mcpoint->GetEnergyLoss()*1.e6); // in keV
282  //hist_eloss->Fill(mcpoint->GetEnergyLoss()*1.e9/50.); // in eV
283 
284  int sensID = mcpoint->GetSensorID();
285  int ihalf, iplane, imodule, iside, idie, isensor;
286  lmddim.Get_sensor_by_id(sensID, ihalf, iplane, imodule, iside, idie, isensor);
287  TVector3 hit_pos = mcpoint->GetPosition();
288  TVector3 hit_pos_lmd = lmddim.Transform_global_to_lmd_local(hit_pos, false);
289 
290 
291  nhitstotal++;
292 
293  double X = hit_pos_lmd.X()*0.01;
294  double Y = hit_pos_lmd.Y()*0.01;
295  //cout << hit_pos.Z() << endl;
296  double RR = X*X+Y*Y;
297 
298  /*
299  if (Rmin_cuts[iplane] < RR && RR < Rmax_cuts[iplane]){//(Rmin < RR && RR < Rmax){
300  hist_XY_plane[iplane]->Fill(X, Y);
301  }*/
302  }
303  }
304 
305  offsetx /= (double) npbar;
306  offsety /= (double) npbar;
307  TVector3 offset_ip(offsetx, offsety, 0);
308  offsetx *= 10.; // mm
309  offsety *= 10.; // mm
310  tiltx /= (double) npbar;
311  tilty /= (double) npbar;
312  TVector3 tilt_ip(tiltx, tilty, 1.);
313  tiltx *= 1000.; // mrad
314  tilty *= 1000.; // mrad
315  momentum /= (double) npbar;
316 
317  // propagate properties to lmd
318  TVector3 tilt_lmd = tilt_ip;
319  TVector3 offset_lmd = offset_ip;
320  cout << " before propagation " << endl;
321  offset_lmd.Print();
322  tilt_lmd.Print();
323  cout << endl;
324 
325  lmddim.Propagate_fast_ip_to_lmd(offset_lmd, tilt_lmd, pbeam);
326  cout << " after propagation " << endl;
327  offset_lmd.Print();
328  tilt_lmd.Print();
329  cout << endl;
330 
331  offset_lmd = lmddim.Transform_global_to_lmd_local(offset_lmd, false, true);
332  tilt_lmd = lmddim.Transform_global_to_lmd_local(tilt_lmd, true, true);
333  cout << " in lmd frame " << endl;
334  offset_lmd.Print();
335  tilt_lmd.Print();
336  cout << endl;
337 
338  offsetx_lmd = offset_lmd.X()/10.;
339  offsety_lmd = offset_lmd.Y()/10.;
340  tiltx_lmd = tilt_lmd.X()/tilt_lmd.Z() * 1000.;
341  tilty_lmd = tilt_lmd.Y()/tilt_lmd.Z() * 1000.;
342 
343 
344  cout << " analyzing histograms for momentum " << momentum << endl;
345  cout << " mean pos of ip was (generated)" << endl;
346  cout << " x = " << hist_gen_x_y.GetMean(1) << " mm ("<< offsetx <<" mm)" << endl;
347  cout << " y = " << hist_gen_x_y.GetMean(2) << " mm ("<< offsety <<" mm)" << endl;
348 
349  cout << " mean dir of beam was (generated)" << endl;
350  cout << " theta = " << hist_gen_theta_phi.GetMean(1) << " mrad ("<< tiltx <<" mrad)" << endl;
351  cout << " phi = " << hist_gen_theta_phi.GetMean(2) << " rad ("<< tilty <<" mrad)" << endl;
352  cout << " dx/dz = " << hist_gen_dx_dy.GetMean(1) << " mrad ("<< tiltx <<" mrad)" << endl;
353  cout << " dy/dz = " << hist_gen_dx_dy.GetMean(2) << " mrad ("<< tilty <<" mrad)" << endl;
354 
355  cout << " mean dir of beam was (propagated to lmd)" << endl;
356  //cout << " theta = " << " mrad ("<< theta_lmd <<" mrad)" << endl;
357  //cout << " phi = " << " rad ("<< phi_lmd <<" rad)" << endl;
358  cout << " dx/dz = " << " mrad ("<< tiltx_lmd <<" mrad)" << endl;
359  cout << " dy/dz = " << " mrad ("<< tilty_lmd <<" mrad)" << endl;
360 
361  cout << " mean in plane 0,1,2,3 is " << endl;
362  vector < double > means_x;
363  vector < double > widths_x;
364  vector < double > means_y;
365  vector < double > widths_y;
366  vector < double > z_err; z_err.push_back(0); z_err.push_back(0); z_err.push_back(0); z_err.push_back(0);
367  means_x.push_back(hist_XY_plane_0.GetMean(1)); meanx_0 = means_x[0];
368  means_x.push_back(hist_XY_plane_1.GetMean(1)); meanx_1 = means_x[1];
369  means_x.push_back(hist_XY_plane_2.GetMean(1)); meanx_2 = means_x[2];
370  means_x.push_back(hist_XY_plane_3.GetMean(1)); meanx_3 = means_x[3];
371  means_y.push_back(hist_XY_plane_0.GetMean(2)); meany_0 = means_x[0];
372  means_y.push_back(hist_XY_plane_1.GetMean(2)); meany_1 = means_x[1];
373  means_y.push_back(hist_XY_plane_2.GetMean(2)); meany_2 = means_x[2];
374  means_y.push_back(hist_XY_plane_3.GetMean(2)); meany_3 = means_x[3];
375 
376  widths_x.push_back(hist_XY_plane_0.GetRMS(1)/sqrt(hist_XY_plane_0.GetEntries())); rmsx_0 = hist_XY_plane_0.GetRMS(1);
377  widths_x.push_back(hist_XY_plane_1.GetRMS(1)/sqrt(hist_XY_plane_1.GetEntries())); rmsx_1 = hist_XY_plane_1.GetRMS(1);
378  widths_x.push_back(hist_XY_plane_2.GetRMS(1)/sqrt(hist_XY_plane_2.GetEntries())); rmsx_2 = hist_XY_plane_2.GetRMS(1);
379  widths_x.push_back(hist_XY_plane_3.GetRMS(1)/sqrt(hist_XY_plane_3.GetEntries())); rmsx_3 = hist_XY_plane_3.GetRMS(1);
380  widths_y.push_back(hist_XY_plane_0.GetRMS(2)/sqrt(hist_XY_plane_0.GetEntries())); rmsy_0 = hist_XY_plane_0.GetRMS(2);
381  widths_y.push_back(hist_XY_plane_1.GetRMS(2)/sqrt(hist_XY_plane_1.GetEntries())); rmsy_1 = hist_XY_plane_1.GetRMS(2);
382  widths_y.push_back(hist_XY_plane_2.GetRMS(2)/sqrt(hist_XY_plane_2.GetEntries())); rmsy_2 = hist_XY_plane_2.GetRMS(2);
383  widths_y.push_back(hist_XY_plane_3.GetRMS(2)/sqrt(hist_XY_plane_3.GetEntries())); rmsy_3 = hist_XY_plane_3.GetRMS(2);
384 
385  //TMarker markeroffset(offsetx, offsety, 4);
386  //markeroffset.SetMarkerColor(2);
387  canvas_distributions.cd(1);
388  hist_XY_plane_0.Draw("COLZ");
389  //markeroffset.Draw();
390  gPad->Update();
391  canvas_distributions.cd(2);
392  hist_XY_plane_1.Draw("COLZ");
393  //markeroffset.Draw();
394  gPad->Update();
395  canvas_distributions.cd(3);
396  hist_XY_plane_2.Draw("COLZ");
397  //markeroffset.Draw();
398  gPad->Update();
399  canvas_distributions.cd(4);
400  hist_XY_plane_3.Draw("COLZ");
401  //markeroffset.Draw();
402  gPad->Update();
403 
404  // save values to tree
405 
406 
407  cout << " x = " << meanx_0 << ", "<< meanx_1 << ", "<< meanx_2 << ", "<< meanx_3 << endl;
408  cout << " y = " << meany_0 << ", "<< meany_1 << ", "<< meany_2 << ", "<< meany_3 << endl;
409  cout << " rms x = " << rmsx_0 << ", "<< rmsx_1 << ", "<< rmsx_2 << ", "<< rmsx_3 << endl;
410  cout << " rms y = " << rmsy_0 << ", "<< rmsy_1 << ", "<< rmsy_2 << ", "<< rmsy_3 << endl;
411 
412  TGraphErrors graph_average_x(means_x.size(), &z[0], &means_x[0], NULL, &widths_x[0]);
413  graph_average_x.SetNameTitle("graph_average_x", "average x over z");
414  canvas.cd(1);
415  gPad->Clear();
416  graph_average_x.Draw("ALP");
417  graph_average_x.GetYaxis()->SetRangeUser(-0.02, 0.02);
418  TF1 funcx("funcx", "pol1", -0.1, 0.5);
419  graph_average_x.Fit(&funcx);
420  //line.DrawLine(z[0], offsetx*1e-3, z[3], offsetx*1e-3);
421  //line.DrawLine(z[0], tiltx*1e-3*z[0]+offsetx*1e-3, z[3], tiltx*1e-3*z[3]+offsetx*1e-3);
422  //line.DrawLine(z[0], offsetx_lmd/100., z[3], offsetx_lmd/100.);
423  line.DrawLine(z[0], offsetx_lmd/100., z[3], offsetx_lmd/100.+tiltx_lmd/1000.*(z[3]-z[0]));
424  gPad->Update();
425  //sleep(2);
426  //graph_average_x.Write();
427  canvas.cd(2);
428  TGraphErrors graph_average_y(means_y.size(), &z[0], &means_y[0], NULL, &widths_y[0]);
429  graph_average_y.SetNameTitle("graph_average_y", "average y over z");
430  gPad->Clear();
431  graph_average_y.Draw("ALP");
432  graph_average_y.GetYaxis()->SetRangeUser(-0.02, 0.02);
433  TF1 funcy("funcy", "pol1", -0.1, 0.5);
434  graph_average_y.Fit(&funcy);
435  //line.DrawLine(z[0], offsety*1e-3, z[3], offsety*1e-3);
436  //line.DrawLine(z[0], tilty*1e-3*z[0]+offsety*1e-3, z[3], tilty*1e-3*z[3]+offsety*1e-3);
437  //line.DrawLine(z[0], offsety_lmd/100., z[3], offsety_lmd/100.);
438  line.DrawLine(z[0], offsety_lmd/100., z[3], offsety_lmd/100.+tilty_lmd/1000.*(z[3]-z[0]));
439  //graph_average_y.Write();
440  gPad->Update();
441 
442  cout << " mean pos of ip fit (generated and propagated)" << endl;
443  offsetx_fit = funcx.Eval(z[0])*1e3;// GetParameter(0)*1e3;
444  offsety_fit = funcy.Eval(z[0])*1e3;// GetParameter(0)*1e3;
445  cout << " x = " << offsetx_fit << " mm ("<< offsetx_lmd <<" mm)" << endl;
446  cout << " y = " << offsety_fit << " mm ("<< offsety_lmd <<" mm)" << endl;
447 
448  tiltx_fit = funcx.GetParameter(1)*1e3;
449  tilty_fit = funcy.GetParameter(1)*1e3;
450  cout << " mean dir of beam fit (generated and propagated)" << endl;
451  cout << " dx/dz = " << tiltx_fit << " mrad ("<< tiltx_lmd <<" mrad)" << endl;
452  cout << " dy/dz = " << tilty_fit << " mrad ("<< tilty_lmd <<" mrad)" << endl;
453 
454  cout << endl;
455  results.Fill();
456  //file.Write();
457  //file.Close();
458  }
459 
460  TFile treefile("asymmetries_results.root", "RECREATE");
461  results.SetDirectory(&treefile);
462  results.Write();
463  treefile.Close();
464  cout << " done " << endl;
465  myapp.Run();
466  return 0;
467 }
468 
469 void Root_Appearance(){
470  TPad foo; // never remove this line :-)))
471  if(1){
472  gROOT->SetStyle("Plain");
473  const Int_t NRGBs = 5;
474  const Int_t NCont = 255;
475  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
476  Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
477  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
478  Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
479  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
480  gStyle->SetNumberContours(NCont);
481  gStyle->SetTitleFont(10*13+2,"xyz");
482  gStyle->SetTitleSize(0.06, "xyz");
483  gStyle->SetTitleOffset(1,"xy");
484  gStyle->SetTitleOffset(1.3,"z");
485  gStyle->SetLabelFont(10*13+2,"xyz");
486  gStyle->SetLabelSize(0.06,"xyz");
487  gStyle->SetLabelOffset(0.009,"xyz");
488  gStyle->SetPadBottomMargin(0.2);
489  gStyle->SetPadTopMargin(0.15);
490  gStyle->SetPadLeftMargin(0.15);
491  gStyle->SetPadRightMargin(0.20);
492  gStyle->SetOptTitle(1);
493  gStyle->SetOptStat(1);
494  gROOT->ForceStyle();
495  gStyle->SetFrameFillColor(0);
496  gStyle->SetFrameFillStyle(0);
497  TGaxis::SetMaxDigits(3);
498  }
499 }
500 
501 void DrawProgressBar(int len, double percent) {
502  if ((int) (last_percent * 100) == (int) (percent * 100))
503  return;
504  //cout << " drawing " << endl;
505  cout << "\x1B[2K"; // Erase the entire current line.
506  cout << "\x1B[0E"; // Move to the beginning of the current line.
507  string progress;
508  for (int i = 0; i < len; ++i) {
509  if (i < static_cast<int> (len * percent)) {
510  progress += "=";
511  } else {
512  progress += " ";
513  }
514  }
515  cout << "[" << progress << "] " << (static_cast<int> (100 * percent))
516  << "%";
517  flush( cout); // Required.
518  last_percent = percent;
519 }
PndLmdDim::Propagate_fast_ip_to_lmd
void Propagate_fast_ip_to_lmd(TVector3 &pos, TVector3 &mom, double pbeam)
Definition: PndLmdDim.cxx:2547
i
Int_t i
Definition: run_full.C:25
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
PndMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
last_percent
double last_percent(-1.)
PndLmdDim::Transform_global_to_lmd_local
TVector3 Transform_global_to_lmd_local(const TVector3 &point, bool isvector=false, bool aligned=true)
Definition: PndLmdDim.cxx:3113
Root_Appearance
void Root_Appearance()
Definition: online_monitoring_studies.C:469
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
online_monitoring_studies
int online_monitoring_studies()
Definition: online_monitoring_studies.C:41
main
int main()
Definition: online_monitoring_studies.C:36
DrawProgressBar
void DrawProgressBar(int len, double percent)
Definition: online_monitoring_studies.C:501
files
std::map< int, TString > files
Definition: simubg.C:28
Y
double Y
Definition: anaLmdDigi.C:68
PndSdsMCPoint::GetSensorID
Int_t GetSensorID() const
Definition: PndSdsMCPoint.h:89
PndMCTrack::IsGeneratorCreated
Bool_t IsGeneratorCreated(void) const
Definition: PndMCTrack.h:84
nHits
int nHits
Definition: RiemannTest.C:16
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
PndLmdDim::Get_sensor_by_id
void Get_sensor_by_id(const int sensor_id, int &ihalf, int &iplane, int &imodule, int &iside, int &idie, int &isensor)
Definition: PndLmdDim.h:526
nEvents
Int_t nEvents
Definition: hit_dirc.C:11
z
Double_t z
Definition: createRootFscGeometryFile.C:146
PndSdsMCPoint::GetPosition
TVector3 GetPosition() const
Definition: PndSdsMCPoint.h:90
X
double X
Definition: anaLmdDigi.C:68
PndLmdDim::Get_instance
static PndLmdDim & Get_instance()
Definition: PndLmdDim.cxx:242
PndMCTrack::GetStartVertex
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
filename
const string filename
Definition: Modify_root_geometry.C:23
Generated on Sun Aug 18 2019 03:01:16 for FairRoot/PandaRoot by   doxygen 1.8.8