FairRoot/PandaRoot
trafo_matrix_fit.C
Go to the documentation of this file.
1 #include <iostream>
2 #include <sstream>
3 
4 #include<TCanvas.h>
5 #include<TApplication.h>
6 #include<TROOT.h>
7 #include<TChain.h>
8 #include<TStyle.h>
9 #include<TGaxis.h>
10 #include<TColor.h>
11 #include<TClonesArray.h>
12 #include<TH1.h>
13 #include<TH2.h>
14 #include<vector>
15 #include<TTree.h>
16 #include<sstream>
17 #include<string>
18 
19 #include<PndLmdDim.h>
20 #include<PndMCTrack.h>
21 #include<PndSdsMCPoint.h>
22 #include<TGraphErrors.h>
23 
24 #include "Math/Minimizer.h"
25 #include "Math/Factory.h"
26 #include "Math/Functor.h"
27 #include "TError.h"
28 #include "TLine.h"
29 
30 #include <iomanip>
31 
32 using namespace std;
33 
34 // change some stuff in the ROOT appearance
35 void Root_Appearance();
36 
37 double last_percent(-1.);
38 // draw a progress bar only when the length changes significantly
39 void DrawProgressBar(int len, double percent);
40 
41 // the application
42 // Assuming to find an mc file called
43 // ./Lumi_MC_0.root
44 // if invert, then the solution for a transformation
45 // from lmd to ip is searched
46 int trafo_matrix_fit(bool invert = false);
47 
48 // calling all vectors to minimize the parameters xx
49 // (Minimizer function)
50 double Chi2(const double *xx);
51 
52 // load vectors from a given MC file
53 // vectors from hits at the first LMD plane
54 // with corresponding hits at the IP will be loaded
55 // if invert, then the solution for a transformation
56 // from lmd to ip is searched
57 void load_vectors(string filename, bool invert = false);
58 
59 // test loaded vectors on a transformation matrix
60 // see code of Chi2() for the structure explanation of the matrix xx
61 // provide an output file name for the results
62 void trafo_matrix_test(const double *xx, string filename);
63 
64 // definition of the vector of parameters
65 // containing the homogeneous space coordinates
66 // and the 3 momenta of particles
67 // x,y,z,1,px/pz,py/pz,1
68 class TRepVect{
69 public:
70  TRepVect(){
71  for (int i = 0; i < 7; i++){
72  entry[i] = 0;
73  }
74  }
75  double entry[7];
76  double Mag2(){
77  double result = 0.;
78  for (int i = 0; i < 7; i++){
79  if (i < 4) result += entry[i]*entry[i];
80  // put a penalty on the normalized momentum direction vectors
81  // to have a good measure of the track direction
82  else result += entry[i]*entry[i];
83  }
84  return result;
85  }
86  double& operator[] (int i) {return entry[i];}
87  TRepVect operator- (const TRepVect vect){
88  TRepVect result;
89  for (int i = 0; i < 7; i++){
90  result.entry[i] = this->entry[i]-vect.entry[i];
91  }
92  return result;
93  }
94 };
95 
96 class TRepVectPair{
97 public:
98  TRepVect IP;
99  TRepVect LMD;
100  bool valid;
101  void invert(){
102  TRepVect _temp = IP;
103  IP = LMD;
104  LMD = _temp;
105  }
106 };
107 
108 // containing a 7 x 7 transformation matrix for TRepVect
109 /* Assuming the following structure of the matrix
110  * which is used to transform homogeneous coordinates
111  * as well as momenta directions
112  * (x,y,z,1,px/pz,py/pz,1)^T = xi
113  */
114 // The corresponding transformation would be
115 // xi' = Mij * xj
116 // where xi' are the vector parameters at the LMD
117 // and xj the vector parameters at the IP
118 // xi' will be written as xi
119 // Mij * xj - xi' should be ideally 0
120 // what will be used as an minimization criterion
121 class TRepMat{
122 public:
123  double M[7][7];
124  TRepMat(){
125  for (int i = 0; i < 7; i++){
126  for (int j = 0; j < 7; j++){
127  M[i][j] = 0;
128  }
129  }
130  }
131  // constructor for 49 entries in a row
132  // to be matched on the matrix
133  TRepMat(const double* mat){
134  for (int i = 0; i < 7*7; i++){
135  (*this)[i] = mat[i];
136  }
137  }
138  double& operator[](int const ij) {
139  return M[ij/7][ij%7];
140  }
141  double& GetSet(int col, int row){
142  return M[col][row];
143  }
144  double* Get(){
145  double* result = new double[7*7];
146  for (int i = 0; i < 7*7; i++){
147  result[i] = (*this)[i];
148  }
149  return result;
150  }
151  TRepVect operator*(const TRepVect vect){
152  TRepVect result;
153  for (int i = 0; i < 7; i++){
154  for (int j = 0; j < 7; j++){
155  result.entry[i] += this->M[i][j] * vect.entry[j];
156  }
157  }
158  return result;
159  }
160  TRepMat operator*(const TRepMat mat){
161  TRepMat result;
162  for (int i = 0; i < 7; i++){
163  for (int j = 0; j < 7; j++){
164  for (int k = 0; k < 7; k++){
165  result.M[i][j] += this->M[i][k] * mat.M[k][j];
166  }
167  }
168  }
169  return result;
170  }
171  void Print(){
172  cout << " The numbers are " << endl;
173  for (int i = 0; i < 7; i++){
174  for (int j = 0; j < 7; j++){
175  cout << std::setprecision(4) << "\t" << M[i][j] << ", ";
176  }
177  cout << endl;
178  }
179  }
180 };
181 
182 vector<TRepVectPair> repvectpairs;
183 
184 // see what is going on during minimization
185 TGraph convergence_graph;
186 TH1D Chi2_dist("Chi2_dist", "chi2 distribution", 100, 0, 0.002);
187 TLine meanline;
188 int ncalls(0);
189 
190 int main() {
191  trafo_matrix_fit(false);
192  return 0;
193 }
194 
195 double Chi2(const double *xx)
196 {
197  Chi2_dist.Reset();
198 /*
199  double Mij[7][7] = {
200  {xx[0],xx[1],xx[2],xx[3],xx[4],xx[5],xx[6]},
201  {xx[7],xx[8],xx[9],xx[10],xx[11],xx[12],xx[13]},
202  {xx[14],xx[15],xx[16],xx[17],xx[18],xx[19],xx[20]},
203  {xx[21],xx[22],xx[23],xx[24],xx[25],xx[26],xx[27]},
204  {xx[28],xx[29],xx[30],xx[31],xx[32],xx[33],xx[34]},
205  {xx[35],xx[36],xx[37],xx[38],xx[39],xx[40],xx[41]},
206  {xx[42],xx[43],xx[44],xx[45],xx[46],xx[47],xx[48]}
207  };
208  */
209  TRepMat Mij(xx);
210  double mean_diff = 0;
211  int nvectors = repvectpairs.size();
212  for (int ivector = 0; ivector < nvectors; ivector++){
213  TRepVect x0 = Mij*repvectpairs[ivector].IP;
214  x0 = x0 - repvectpairs[ivector].LMD;
215  //cout << endl;
216  //for (int i = 0; i < 7; i++){
217  // cout << " " << repvectpairs[ivector].IP[i] << endl;
218  //}
219  //for (int i = 0; i < 7; i++){
220  // for (int j = 0; j < 7; j++){
221  // cout << "\t" << Mij[i][j];
222  // }
223  // cout << endl;
224  //}
225  /*
226  for (int i = 0; i < 7; i++){
227  //cout << endl << i << " : " << repvectpairs[ivector].IP[i] << " becomes ";
228  for (int j = 0; j < 7; j++){
229  x0[i] += Mij[i][j] * repvectpairs[ivector].IP[j];
230  }
231  //cout << x0[i] << " and should be ";
232  x0[i] -= repvectpairs[ivector].LMD[i];
233  //cout << endl << repvectpairs[ivector].LMD[i] << endl;
234  }*/
235  //for (int i = 0; i < 7; i++){
236  // cout << " " << x0[i] << "\t" << repvectpairs[ivector].LMD[i] << endl;
237  //}
238  // calculate the length of the vector
239  Chi2_dist.Fill(x0.Mag2());
240  //cout << x0.Mag2() << endl;
241  mean_diff += x0.Mag2()/(double) nvectors;
242 
243  //break;
244  }
245  // normalize
246  //cout << " the difference before normalization is " << mean_diff << endl;
247  //cout << mean_diff << " " << endl;
248  //mean_diff /= (double) repvectpairs.size();
249  //cout << " " << mean_diff << endl;
250  /*cout << " The numbers are " << endl;
251  for (int i = 0; i < 7; i++){
252  for (int j = 0; j < 7; j++){
253  cout << "\t" << Mij[i][j];
254  }
255  cout << endl;
256  }*/
257  //cout << " the difference is " << mean_diff << endl;
258 
259  ncalls++;
260  if (ncalls > convergence_graph.GetN()){
261  convergence_graph.Set(ncalls+100);
262  for (int i = ncalls; i < ncalls + 100; i++)
263  convergence_graph.SetPoint(i, 0, convergence_graph.GetY()[0]);
264  }
265  convergence_graph.SetPoint(ncalls-1, ncalls, mean_diff);
266  if (ncalls%100 == 0){
267  Chi2_dist.Draw();
268  //convergence_graph.Draw("A*P");
269  gPad->SetLogy();
270  meanline.SetX1(mean_diff);
271  meanline.SetX2(mean_diff);
272  meanline.SetY1(1e-4);
273  meanline.SetY2(1000);
274  meanline.SetLineColor(2);
275  meanline.SetLineWidth(2);
276  meanline.Draw();
277  gPad->Update();
278  //sleep(1);
279  }
280  //mean_diff /= (double) repvectpairs.size();
281 
282  return Chi2_dist.GetMean();
283  return mean_diff; // the mean diff seems not a good measure without a cut such as the histogram is performing
284 }
285 
286 void trafo_matrix_test(const double *xx, string filename){
287  TCanvas canvas("canvas", "canvas", 600, 600);
288  TH1F hist_res_x("hist_res_x", "resolution x", 500, -0.005, 0.005);
289  TH1F hist_res_y("hist_res_y", "resolution y", 500, -0.005, 0.005);
290  TH1F hist_res_z("hist_res_z", "resolution z", 500, -0.005, 0.005);
291  TH1F hist_res_w("hist_res_w", "resolution w", 500, -0.005, 0.005);
292  TH1F hist_res_pxpz("hist_res_pxpz", "resolution px/pz", 500, -0.005, 0.005);
293  TH1F hist_res_pypz("hist_res_pypz", "resolution py/pz", 500, -0.005, 0.005);
294  TH1F hist_res_pw("hist_res_pw", "resolution pw", 500, -0.005, 0.005);
295  //TH1F hist_res_theta("hist_res_theta", "resolution #theta", 1000, -0.1, 0.1);
296  //TH1F hist_res_phi("hist_res_phi", "resolution #phi", 1000, -0.1, 0.1);
297  TH1F* hists[7] = {&hist_res_x, &hist_res_y, &hist_res_z, &hist_res_w, &hist_res_pxpz, &hist_res_pypz, &hist_res_pw};
298  TRepMat Mij(xx);
299  /*
300  double Mij[7][7] = {
301  {xx[0],xx[1],xx[2],xx[3],xx[4],xx[5],xx[6]},
302  {xx[7],xx[8],xx[9],xx[10],xx[11],xx[12],xx[13]},
303  {xx[14],xx[15],xx[16],xx[17],xx[18],xx[19],xx[20]},
304  {xx[21],xx[22],xx[23],xx[24],xx[25],xx[26],xx[27]},
305  {xx[28],xx[29],xx[30],xx[31],xx[32],xx[33],xx[34]},
306  {xx[35],xx[36],xx[37],xx[38],xx[39],xx[40],xx[41]},
307  {xx[42],xx[43],xx[44],xx[45],xx[46],xx[47],xx[48]}
308  };*/
309  int nvectors = repvectpairs.size();
310  for (int ivector = 0; ivector < nvectors; ivector++){
311  //TRepVect x0;
312  //TVector3 dir_lmd(repvectpairs[ivector].LMD[4], repvectpairs[ivector].LMD[5], repvectpairs[ivector].LMD[6]);
313  //TVector3 dir_after;
314  TRepVect x0 = Mij*repvectpairs[ivector].IP;
315  x0 = x0 - repvectpairs[ivector].LMD;
316  for (int i = 0; i < 7; i++){
317  //for (int j = 0; j < 7; j++){
318  // x0[i] += Mij[i][j] * repvectpairs[ivector].IP[j];
319  //}
320  //if (i == 4) dir_after.SetX(x0[i]);
321  //if (i == 5) dir_after.SetY(x0[i]);
322  //if (i == 6) dir_after.SetZ(x0[i]);
323  //x0[i] -= repvectpairs[ivector].LMD[i];
324  hists[i]->Fill(x0[i]);
325  }
326  //hist_res_theta.Fill(dir_after.Theta()-dir_lmd.Theta());
327  //hist_res_phi.Fill(dir_after.Phi()-dir_lmd.Phi());
328  }
329  //hist_res_theta.Draw("hist");
330  //canvas.Print((filename+"(").c_str());
331  //hist_res_phi.Draw("hist");
332  //canvas.Print((filename+"(").c_str());
333  for (int i = 0; i < 7; i++){
334  hists[i]->Draw("hist");
335  if (i < 6) canvas.Print((filename+"(").c_str());
336  else canvas.Print((filename+")").c_str());
337  }
338 // canvas.Print((filename.str()+"]").c_str());
339 }
340 
341 void load_vectors(string filename, bool invert){
342  PndLmdDim& lmddim = PndLmdDim::Get_instance();
343  TChain tMC("pndsim");
344  tMC.Add(filename.c_str());
345 
346  //--- assign MC info -----------------------------------------------------
347  TClonesArray* true_tracks = new TClonesArray("PndMCTrack");
348  tMC.SetBranchAddress("MCTrack", &true_tracks); //True Track to compare
349 
350  TClonesArray* true_points = new TClonesArray("PndSdsMCPoint");
351  tMC.SetBranchAddress("LMDPoint", &true_points); //True Points to compare
352 
353  // ******** testing
354  TCanvas canvas("canvas_test", "canvas_test", 600, 600);
355  TH1F hist_res_x("hist_res_x", "resolution x", 500, -0.1, 0.1);
356  TH1F hist_res_y("hist_res_y", "resolution y", 500, -0.1, 0.1);
357  TH1F hist_res_z("hist_res_z", "resolution z", 500, -0.5, 0.5);
358  TH1F hist_res_px("hist_res_px", "resolution px", 500, -0.001, 0.001);
359  TH1F hist_res_py("hist_res_py", "resolution py", 500, -0.001, 0.001);
360  TH1F hist_res_pz("hist_res_pz", "resolution pz", 500, -0.001, 0.001);
361  TH1F hist_res_theta("hist_res_theta", "resolution theta", 500, -0.001, 0.001);
362  TH1F hist_res_phi("hist_res_phi", "resolution phi", 500, -0.1, 0.1);
363  // ******** end testing
364 
365  int nEvents = tMC.GetEntries();
366  cout << " reading " << nEvents << " Events " << endl;
367  for (Int_t j = 0; j < nEvents ; j++) {
368  DrawProgressBar(50, (j + 1) / ((double) nEvents));
369  tMC.GetEntry(j);
370 
371  // For testing purposes only
372  TVector3 pos_ip_test; // for testing only
373  TVector3 mom_ip_test; // for testing only
374  TVector3 pos_lmd_test; // for testing only
375  TVector3 mom_lmd_test; // for testing only
376  // end of "testing purposes only"
377 
378  int nHits = true_points->GetEntriesFast();
379  int nSdsHits = 0;
380 
381  const int nParticles = true_tracks->GetEntriesFast();
382 
383  TRepVectPair repvectpair;
384  repvectpair.valid = false;
385 
386  // get the beam properties of the IP
387  for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
388  PndMCTrack *mctrk = (PndMCTrack*) true_tracks->At(iParticle);
389  Int_t mcID = mctrk->GetPdgCode();
390  if (mctrk->IsGeneratorCreated() && mcID == -2212) {
391  //npbar++;
392  TVector3 momMC = mctrk->GetMomentum();
393  mom_ip_test = momMC;
394  //tiltx += momMC.X()/momMC.Z();
395  //tilty += momMC.Y()/momMC.Z();
396  //momentum += momMC.Mag();
397  TVector3 posMC = mctrk->GetStartVertex();
398  pos_ip_test = posMC;
399  //offsetx += posMC.X();
400  //offsety += posMC.Y();
401  repvectpair.IP[0] = posMC.X()/100.;
402  repvectpair.IP[1] = posMC.Y()/100.;
403  repvectpair.IP[2] = posMC.Z()/100.;
404  repvectpair.IP[3] = 1.;
405  repvectpair.IP[4] = momMC.X()/momMC.Z()*10.; // scale by 10 in order to weight its importance
406  repvectpair.IP[5] = momMC.Y()/momMC.Z()*10.;
407  repvectpair.IP[6] = 1.;//momMC.Z()/momMC.Mag();
408  break;
409  }
410  }
411 
412  for (Int_t iHit = 0; iHit < nHits; iHit++) {
413  PndSdsMCPoint* mcpoint = (PndSdsMCPoint*) true_points->At(iHit);
414  if (mcpoint->GetTrackID() < 0) continue;
415  PndMCTrack *mctrk = (PndMCTrack*) true_tracks->At(mcpoint->GetTrackID());
416  if (mctrk->GetPdgCode() == -2212){
417 
418  }
419  // cout << mcpoint->GetEnergyLoss() << endl;
420  //hist_eloss_total->Fill(mcpoint->GetEnergyLoss()*1.e6); // in keV
421  //hist_eloss->Fill(mcpoint->GetEnergyLoss()*1.e9/50.); // in eV
422 
423  int sensID = mcpoint->GetSensorID();
424  int ihalf, iplane, imodule, iside, idie, isensor;
425  lmddim.Get_sensor_by_id(sensID, ihalf, iplane, imodule, iside, idie, isensor);
426  if (iplane > 0) continue; // take only hits on the very first plane
427  TVector3 posMC = mcpoint->GetPosition();
428  TVector3 momMC(mcpoint->GetPx(), mcpoint->GetPy(), mcpoint->GetPz());// = mctrk->GetMomentum();
429  pos_lmd_test = posMC;
430  mom_lmd_test = momMC;
431  repvectpair.LMD[0] = posMC.X()/100.;
432  repvectpair.LMD[1] = posMC.Y()/100.;
433  repvectpair.LMD[2] = posMC.Z()/100.;
434  repvectpair.LMD[3] = 1.;
435  repvectpair.LMD[4] = momMC.X()/momMC.Z()*10.;
436  repvectpair.LMD[5] = momMC.Y()/momMC.Z()*10.;
437  repvectpair.LMD[6] = 1.;//momMC.Z()/momMC.Mag();
438  repvectpair.valid = true;
439  break;
440  //TVector3 hit_pos_lmd = lmddim.Transform_global_to_lmd_local(hit_pos, false);
441  //nhitstotal++;
442  }
443  if (repvectpair.valid){
444  if (invert) repvectpair.invert();
445  repvectpairs.push_back(repvectpair);
446  /*
447  cout << endl;
448  for (int i = 0; i < 7; i++){
449  cout << i << " " << repvectpair.IP[i] << " ";
450  cout << i << " " << repvectpair.LMD[i] << " ";
451  }
452  cout << endl;*/
453  // testing fast tracking of lmd dim which is based on the results from this fit
454  //cout << endl;
455  //pos_ip_test.Print();
456  //mom_ip_test.Print();
457  lmddim.Propagate_fast_ip_to_lmd(pos_ip_test, mom_ip_test, 1.5); // modfy last number according to beam momentum
458  //cout << endl;
459  //pos_ip_test.Print();
460  //pos_lmd_test.Print();
461  //mom_ip_test.Print();
462  //mom_lmd_test.Print();
463 
464  double theta_test = mom_ip_test.Theta() - mom_lmd_test.Theta();
465  double phi_test = mom_ip_test.Phi() - mom_lmd_test.Phi();
466  pos_lmd_test = pos_ip_test - pos_lmd_test;
467  mom_lmd_test = mom_ip_test - mom_lmd_test;
468  hist_res_x.Fill(pos_lmd_test.X());
469  hist_res_y.Fill(pos_lmd_test.Y());
470  hist_res_z.Fill(pos_lmd_test.Z());
471  hist_res_px.Fill(mom_lmd_test.X());
472  hist_res_py.Fill(mom_lmd_test.Y());
473  hist_res_pz.Fill(mom_lmd_test.Z());
474  hist_res_theta.Fill(theta_test);
475  hist_res_phi.Fill(phi_test);
476  }
477  }
478  cout << endl << " found " << repvectpairs.size() << " vector pairs for analysis " << endl;
479 
480  // testing
481  static int datacounter (0);
482  datacounter++;
483  stringstream outfile;
484  outfile << "test_lmd_" << datacounter << ".pdf";
485 
486  hist_res_x.Draw("hist");
487  canvas.Print((outfile.str()+"(").c_str());
488  hist_res_y.Draw("hist");
489  canvas.Print((outfile.str()+"(").c_str());
490  hist_res_z.Draw("hist");
491  canvas.Print((outfile.str()+"(").c_str());
492  hist_res_px.Draw("hist");
493  canvas.Print((outfile.str()+"(").c_str());
494  hist_res_py.Draw("hist");
495  canvas.Print((outfile.str()+"(").c_str());
496  hist_res_pz.Draw("hist");
497  canvas.Print((outfile.str()+"(").c_str());
498  hist_res_theta.Draw("hist");
499  canvas.Print((outfile.str()+"(").c_str());
500  hist_res_phi.Draw("hist");
501  canvas.Print((outfile.str()+")").c_str());
502 }
503 
504 int trafo_matrix_fit(bool invert){
505  TApplication myapp("myapp", 0, 0);
506  // all 1.5 GeV/c files by stefan
507 
508 load_vectors("./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", invert);
509 load_vectors("./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", invert);
510 load_vectors("./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", invert);
511 load_vectors("./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", invert);
512 load_vectors("./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", invert);
513 load_vectors("./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", invert);
514 load_vectors("./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", invert);
515 load_vectors("./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", invert);
516 load_vectors("./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", invert);
517 load_vectors("./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", invert);
518 load_vectors("./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", invert);
519 load_vectors("./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", invert);
520 load_vectors("./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", invert);
521 load_vectors("./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", invert);
522 load_vectors("./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", invert);
523 load_vectors("./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", invert);
524 load_vectors("./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", invert);
525 load_vectors("./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", invert);
526 load_vectors("./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", invert);
527 load_vectors("./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", invert);
528 load_vectors("./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", invert);
529 load_vectors("./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", invert);
530 load_vectors("./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", invert);
531 load_vectors("./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", invert);
532 load_vectors("./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", invert);
533 load_vectors("./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", invert);
534 load_vectors("./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", invert);
535 load_vectors("./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", invert);
536 load_vectors("./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", invert);
537 load_vectors("./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", invert);
538 load_vectors("./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", invert);
539 load_vectors("./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", invert);
540 load_vectors("./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", invert);
541 load_vectors("./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", invert);
542 load_vectors("./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", invert);
543 load_vectors("./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", invert);
544 //load_vectors("./box_th_0.0-15.0mrad_recoil_corrected/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", invert);
545 //load_vectors("./box_th_0.0-15.0mrad_recoil_corrected/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", invert);
546 //load_vectors("./box_th_0.0-15.0mrad_recoil_corrected/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", invert);
547 
548 // load_vectors("./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*.root", invert);
549 
550  // 15 GeV/c files by stefan
551  //load_vectors("./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*.root", invert);
552 
553  ROOT::Math::Minimizer* min =
554  ROOT::Math::Factory::CreateMinimizer("Minuit2", "");
555  min->SetMaxFunctionCalls(100000);
556  min->SetMaxIterations(100000);
557  min->SetTolerance(0.00001);
558  min->SetPrintLevel(1);
559 
560  ROOT::Math::Functor f(&Chi2, 49);
561  /*
562  double Mtranssol[49] = {
563  1 ,0 ,0 ,0 ,0 ,0 ,0 , // angles are scaled by 10
564  0 ,1 ,0 ,0 ,0 ,0 ,0 , // 2 equals 2m
565  0 ,0 ,1 ,6.5 ,0 ,0 ,0 ,
566  0 ,0 ,0 ,1 ,0 ,0 ,0 ,
567  0 ,0 ,0 ,0 ,1 ,0 ,0 , // 40 mrad times 10
568  0 ,0 ,0 ,0 ,0 ,1 ,0 ,
569  0 ,0 ,0 ,0 ,0 ,0 ,1
570  };
571  double Mrotsol[49] = {
572  cos(.21) ,-sin(.21) ,0 ,0 ,0 ,0 ,0 , // angles are scaled by 10
573  sin(.21) ,cos(.21) ,0 ,0 ,0 ,0 ,0 , // 1.124 equals 11.24m
574  0 ,0 ,1 ,0 ,0 ,0 ,0 ,
575  0 ,0 ,0 ,1 ,0 ,0 ,0 ,
576  0 ,0 ,0 ,0 ,1 ,0 ,0 , // 40 mrad times 10
577  0 ,0 ,0 ,0 ,0 ,1 ,0 ,
578  0 ,0 ,0 ,0 ,0 ,0 ,1
579  };
580  double Mrotdip[49] = {
581  1 ,0 ,40e-3 ,0 ,0 ,0 ,0 , // angles are scaled by 10
582  0 ,1 ,0 ,0 ,0 ,0 ,0 , // 1.124 equals 11.24m
583  -40e-3 ,0 ,1 ,0 ,0 ,0 ,0 ,
584  0 ,0 ,0 ,1 ,0 ,0 ,0 ,
585  0 ,0 ,0 ,0 ,1 ,0 ,400e-3 , // 40 mrad times 10
586  0 ,0 ,0 ,0 ,0 ,1 ,0 ,
587  0 ,0 ,0 ,0 ,0 ,0 ,1
588  };
589  double Mtransdiplmd[49] = {
590  1 ,0 ,0 ,0 ,1.124 ,0 ,0 , // angles are scaled by 10
591  0 ,1 ,0 ,0 ,0 ,1.124 ,0 , // 2 equals 2m
592  0 ,0 ,1 ,4.74 ,0 ,0 ,0 ,
593  0 ,0 ,0 ,1 ,0 ,0 ,0 ,
594  0 ,0 ,0 ,0 ,1 ,0 ,0 , // 40 mrad times 10
595  0 ,0 ,0 ,0 ,0 ,1 ,0 ,
596  0 ,0 ,0 ,0 ,0 ,0 ,1
597  };*/
598 
599  //TRepMat M1(Mtranssol), M2(Mrotsol), M3(Mrotdip), M4(Mtransdiplmd);
600  //TRepMat MSum = M4* M3 * M2 * M1;
601  //M4.Print();
602  //M3.Print();
603  //M2.Print();
604  //M1.Print();
605  // hacking
606  //MSum.M[0][6] = 0;
607  //MSum.M[1][3] = 0;
608  //MSum.Print();
609  //double* Mij = MSum.Get();
610  //double phi = 0.21;
611  /*
612  double Mij[49] = {
613  1 ,0 ,0.04 ,6.5*0.04 ,1.0992 ,-0.234 ,0 , // angles are scaled by 10
614  0 ,1 ,0 ,0 ,+0.234 ,1.0992 ,0 , // 1.124 equals 11.24m
615  -0.04 ,0 ,1 ,11.24 ,0 ,0 ,0 ,
616  0 ,0 ,0 ,1 ,0 ,0 ,0 ,
617  0 ,0 ,0 ,0 ,0.978 ,-0.208 ,0.4 , // 40 mrad times 10
618  0 ,0 ,0 ,0 ,0.208 ,0.978 ,0.0 ,
619  0 ,0 ,0 ,0 ,0 ,0 ,1
620  };*/
621  double Mij[49] = {
622  1 ,0 ,40e-3 ,6.5*40e-3 ,1.124 ,0. ,0 , // angles are scaled by 10
623  0 ,1 ,0 ,0 ,0 ,1.124 ,0 , // 1.124 equals 11.24m
624  -40e-3 ,0 ,1 ,11.24 ,0 ,0 ,0 ,
625  0 ,0 ,0 ,1 ,0 ,0 ,0 ,
626  0 ,0 ,0 ,0 ,1 ,0 ,400e-3 , // 40 mrad times 10
627  0 ,0 ,0 ,0 ,0 ,1 ,0 ,
628  0 ,0 ,0 ,0 ,0 ,0 ,1
629  };
630  if (invert){
631  double Mij_inverted[49] = {
632  /* 1 ,0 ,-0.04 ,4.75*0.04 ,-1.124 ,0 ,0.4*1.124 , // angles are scaled by 10
633  0 ,1 ,0 ,0 ,0 ,-1.124 ,0 , // 1.124 equals 11.24m
634  0.04 ,0 ,1 ,-11.24 ,0 ,0 ,0 ,
635  0 ,0 ,0 ,1 ,0 ,0 ,0 ,
636  0 ,0 ,0 ,0 ,1 ,0 ,-400e-3 , // 40 mrad times 10
637  0 ,0 ,0 ,0 ,0 ,1 ,0 ,
638  0 ,0 ,0 ,0 ,0 ,0 ,1 */
639  0.979707, 0.00855083, -0.957085, 10.5033, -1.14336, -0.00568183,
640  0.45782, 0.0379627, 0.955083, 0.999504, -11.2443, 0.00606571,
641  -1.07678, -0.00239097, 9.70718, -4.38641, 244.35, -2749.01,
642  -0.0322177, 4.76301, 0.0134341, 0.0000190301, 0.000359596,
643  0.000503436, 0.994336, 0.0000179324, -0.00038629, -6.77702e-6,
644  0.147446, -0.229086, 3.19516, -35.9521, 0.974497, 0.28175,
645  -0.389989, 0.105476, 0.0543932, 0.267799, -3.03747, -0.167528,
646  0.943272, 0.0670344, 0.0000190301, 0.000359596, 0.000503436,
647  -0.00566317, 0.0000179324, -0.00038629, 0.999993
648  };
649  for (int i = 0; i < 49; i++)
650  Mij[i] = Mij_inverted[i];
651  }
652  TRepMat MSum(Mij);
653  cout << " the starting condition is " << endl;
654  MSum.Print();
655 /*
656  double Mij[49] = {
657  1 ,0 ,0 ,0 ,0 ,0 ,0 ,
658  0 ,1 ,0 ,0 ,0 ,0 ,0 ,
659  0 ,0 ,1 ,0 ,0 ,0 ,0 ,
660  0 ,0 ,0 ,1 ,0 ,0 ,0 ,
661  0 ,0 ,0 ,0 ,1 ,0 ,0 ,
662  0 ,0 ,0 ,0 ,0 ,1 ,0 ,
663  0 ,0 ,0 ,0 ,0 ,0 ,1
664  };*/
665  //Chi2(Mij);
666  //Chi2(Mij);
667  min->SetFunction(f);
668  double stepsize = 0.001;
669  for (int i = 0; i < 49; i++){
670  stringstream varname;
671  varname << "M" << i/7 << i%7;
672  if (i > 3) stepsize = 0.001;
673  else stepsize = 0.001;
674  //if (Mij[i] == 0) Mij[i] = 0.0001;
675  min->SetVariable(i, varname.str().c_str(), Mij[i], stepsize);
676  }
677  //trafo_matrix_test(0, "residuals_test.pdf");
678 
679  trafo_matrix_test(Mij, "residuals_before.pdf");
680  min->Minimize();
681  TRepMat Mij_min(min->X());
682  TRepMat Mij_err(min->Errors());
683  cout << " The minimized matrix is " << endl;
684  Mij_min.Print();
685  cout << " with the corresponding errors " << endl;
686  Mij_err.Print();
687 
688  trafo_matrix_test(Mij_min.Get(), "residuals_after.pdf");
689  //min->PrintResults();
690  myapp.Run();
691  return 0;
692 }
693 
694 void Root_Appearance(){
695  TPad foo; // never remove this line :-)))
696  if(1){
697  gROOT->SetStyle("Plain");
698  const Int_t NRGBs = 5;
699  const Int_t NCont = 255;
700  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
701  Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
702  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
703  Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
704  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
705  gStyle->SetNumberContours(NCont);
706  gStyle->SetTitleFont(10*13+2,"xyz");
707  gStyle->SetTitleSize(0.06, "xyz");
708  gStyle->SetTitleOffset(1,"xy");
709  gStyle->SetTitleOffset(1.3,"z");
710  gStyle->SetLabelFont(10*13+2,"xyz");
711  gStyle->SetLabelSize(0.06,"xyz");
712  gStyle->SetLabelOffset(0.009,"xyz");
713  gStyle->SetPadBottomMargin(0.2);
714  gStyle->SetPadTopMargin(0.15);
715  gStyle->SetPadLeftMargin(0.15);
716  gStyle->SetPadRightMargin(0.20);
717  gStyle->SetOptTitle(1);
718  gStyle->SetOptStat(1);
719  gROOT->ForceStyle();
720  gStyle->SetFrameFillColor(0);
721  gStyle->SetFrameFillStyle(0);
722  TGaxis::SetMaxDigits(3);
723  }
724 }
725 
726 void DrawProgressBar(int len, double percent) {
727  if ((int) (last_percent * 100) == (int) (percent * 100))
728  return;
729  //cout << " drawing " << endl;
730  cout << "\x1B[2K"; // Erase the entire current line.
731  cout << "\x1B[0E"; // Move to the beginning of the current line.
732  string progress;
733  for (int i = 0; i < len; ++i) {
734  if (i < static_cast<int> (len * percent)) {
735  progress += "=";
736  } else {
737  progress += " ";
738  }
739  }
740  cout << "[" << progress << "] " << (static_cast<int> (100 * percent))
741  << "%";
742  flush( cout); // Required.
743  last_percent = percent;
744 }
row
int row
Definition: anaLmdDigi.C:67
repvectpairs
vector< TRepVectPair > repvectpairs
Definition: trafo_matrix_fit.C:182
x0
Double_t x0
Definition: checkhelixhit.C:70
operator-
RhoError operator-(const RhoError &m1, const RhoError &m2)
Definition: RhoError.cxx:255
PndLmdDim::Propagate_fast_ip_to_lmd
void Propagate_fast_ip_to_lmd(TVector3 &pos, TVector3 &mom, double pbeam)
Definition: PndLmdDim.cxx:2547
ncalls
int ncalls(0)
i
Int_t i
Definition: run_full.C:25
PndMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
col
int col
Definition: anaLmdDigi.C:67
TRepMat::operator*
TRepVect operator*(const TRepVect vect)
Definition: trafo_matrix_fit.C:151
TRepMat::Print
void Print()
Definition: trafo_matrix_fit.C:171
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
meanline
TLine meanline
Definition: trafo_matrix_fit.C:187
PndSdsMCPoint::GetSensorID
Int_t GetSensorID() const
Definition: PndSdsMCPoint.h:89
TRepVect::entry
double entry[7]
Definition: trafo_matrix_fit.C:75
PndMCTrack::IsGeneratorCreated
Bool_t IsGeneratorCreated(void) const
Definition: PndMCTrack.h:84
nHits
int nHits
Definition: RiemannTest.C:16
convergence_graph
TGraph convergence_graph
Definition: trafo_matrix_fit.C:185
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
TRepVect::Mag2
double Mag2()
Definition: trafo_matrix_fit.C:76
load_vectors
void load_vectors(string filename, bool invert=false)
Definition: trafo_matrix_fit.C:341
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
trafo_matrix_test
void trafo_matrix_test(const double *xx, string filename)
Definition: trafo_matrix_fit.C:286
f
TFile * f
Definition: bump_analys.C:12
TRepMat::TRepMat
TRepMat(const double *mat)
Definition: trafo_matrix_fit.C:133
min
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:25
Root_Appearance
void Root_Appearance()
Definition: trafo_matrix_fit.C:694
operator[]
float & operator[](int i)
Definition: P4_F32vec4.h:6
DrawProgressBar
void DrawProgressBar(int len, double percent)
Definition: trafo_matrix_fit.C:726
PndSdsMCPoint::GetPosition
TVector3 GetPosition() const
Definition: PndSdsMCPoint.h:90
TRepMat::M
double M[7][7]
Definition: trafo_matrix_fit.C:123
TRepMat::operator*
TRepMat operator*(const TRepMat mat)
Definition: trafo_matrix_fit.C:160
Chi2_dist
TH1D Chi2_dist("Chi2_dist","chi2 distribution", 100, 0, 0.002)
last_percent
double last_percent(-1.)
PndLmdDim::Get_instance
static PndLmdDim & Get_instance()
Definition: PndLmdDim.cxx:242
Chi2
double Chi2(const double *xx)
Definition: trafo_matrix_fit.C:195
TRepMat::GetSet
double & GetSet(int col, int row)
Definition: trafo_matrix_fit.C:141
trafo_matrix_fit
int trafo_matrix_fit(bool invert=false)
Definition: trafo_matrix_fit.C:504
main
int main()
Definition: trafo_matrix_fit.C:190
PndMCTrack::GetStartVertex
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
TRepMat::Get
double * Get()
Definition: trafo_matrix_fit.C:144
outfile
TString outfile
Definition: createRootFscGeometryFile.C:57
filename
const string filename
Definition: Modify_root_geometry.C:23
TRepMat::operator[]
double & operator[](int const ij)
Definition: trafo_matrix_fit.C:138
Generated on Sun Aug 18 2019 03:02:19 for FairRoot/PandaRoot by   doxygen 1.8.8