FairRoot/PandaRoot
Functions
error_matrix_fit.C File Reference

Go to the source code of this file.

Functions

int error_matrix_fit (Int_t emcGeometry=17, TString fileVersion="1")
 
Double_t fit1 (Double_t *x, Double_t *par)
 
Double_t fit2 (Double_t *x, Double_t *par)
 

Function Documentation

int error_matrix_fit ( Int_t  emcGeometry = 17,
TString  fileVersion = "1" 
)

Definition at line 7 of file error_matrix_fit.C.

References basiclibs(), Bool_t, c1, c2, c3, c4, cluster_array, cluster_energy, cluster_phi, cluster_theta, cos(), Double_t, energy, PndEmcCluster::energy(), f1, f10, f2, f3, f4, fit1(), fit2(), PndMCTrack::GetMomentum(), h2phi, i, max_energy, mctrack, mctrack_array, outFile, phi, phi_diff, photon_momentum, Pi, rootlogon(), PndEmcErrorMatrixParObject::SetErrorMatrix(), sin(), t, theta, TString, PndEmcCluster::where(), x, y, and z.

8 {
9  gROOT->SetStyle("Plain");
10  gROOT->LoadMacro("$VMCWORKDIR/gconfig/rootlogon.C");
11  gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
12  rootlogon();
13  basiclibs();
14 
15  TH2F *h2phi= new TH2F("h2phi","Phi difference",70,0.,7.,100,-2.,2.);
16  TH2F *h2x= new TH2F("h2x","X difference",70,0.,7.,100,-2.,2.);
17  TH2F *h2y= new TH2F("h2y","Y difference",70,0.,7.,100,-2.,2.);
18  TH2F *h2z= new TH2F("h2z","Z difference",70,0.,7.,100,-2.,2.);
19  TH2F *h2energy= new TH2F("h2energy","Energy difference",70,0.,7.,100,-0.2,0.2);
20 
21  TH1F *h1energy_1= new TH1F("h1energy_1","Cluster energy",70,0.,7.);
22  TH1F *h1energy_2= new TH1F("h1energy_2","Cluster energy",70,0.,7.);
23  TH1F *h1energy_3= new TH1F("h1energy_3","Cluster energy",70,0.,7.);
24  TH1F *h1energy_4= new TH1F("h1energy_4","Cluster energy",70,0.,7.);
25 
26  // Barrel
27  TH1F *h1phi= new TH1F("hphi","Phi difference",70,0.,7.);
28  TH1F *h1z= new TH1F("hz","Z difference",70,0.,7.);
29 
30  // Fwd endcap, bwd endcap, shashlyk
31  TH1F *h1x_2= new TH1F("hx_2","X difference",70,0.,7.);
32  TH1F *h1y_2= new TH1F("hy_2","Y difference",70,0.,7.);
33  TH1F *h1x_3= new TH1F("hx_3","X difference",70,0.,7.);
34  TH1F *h1y_3= new TH1F("hy_3","Y difference",70,0.,7.);
35  TH1F *h1x_4= new TH1F("hx_4","X difference",70,0.,7.);
36  TH1F *h1y_4= new TH1F("hy_4","Y difference",70,0.,7.);
37 
38 
39  TVector3 photon_momentum;
40  double cluster_energy, energy;
41  double cluster_theta, cluster_phi, cluster_x, cluster_y, cluster_z; //position of the cluster
42  double theta, phi, x, y, z; // position of the initial particle
43  double x_diff, y_diff, z_diff, phi_diff, e_diff;
44  double max_energy=0;
45  Int_t emc_component;
46  double sigma_e, sigma_x, sigma_y, sigma_z, sigma_phi, sigma_e_e;
47  double sigma_e_error, sigma_x_error, sigma_y_error, sigma_z_error, sigma_phi_error, sigma_e_e_error;
48  TH1D *h_phi, *h_x, *h_y, *h_z, *h_energy;
49 
50  // Flags set if data files for EMC component is present
51  Bool_t barrel_flag=false, fwcap_flag=false, bwcap_flag=false, fsc_flag=false;
52  // Fit parametrs
53  Double_t pars1[10]={0.}, pars2[10]={0.}, pars3[10]={0.}, pars4[10]={0.};
54 
58  emc_component=1;
59  TString outFile="emc_data_";
60  outFile+=emcGeometry;
61  outFile+="_";
62  outFile+=emc_component;
63  outFile+=".root";
64 
65  TFile* file1 = new TFile(outFile); //file you want to analyse
66  if (!file1->IsZombie())
67  {
68  barrel_flag=true;
69  h2phi->Reset();h2x->Reset();h2y->Reset();h2z->Reset();h2energy->Reset();
70  TTree *t=(TTree *) file1->Get("pndsim") ;
71  TClonesArray* cluster_array=new TClonesArray("PndEmcCluster");
72  t->SetBranchAddress("EmcCluster",&cluster_array);
73  TClonesArray* mctrack_array=new TClonesArray("PndMCTrack");
74  t->SetBranchAddress("MCTrack",&mctrack_array);
75 
76  // Cluster energy
77  for (Int_t j=0; j< t->GetEntriesFast(); j++)
78  {
79  t->GetEntry(j);
80 
81  PndMCTrack *mctrack=(PndMCTrack *) mctrack_array->At(0);
82  photon_momentum=mctrack->GetMomentum();
83  energy=photon_momentum.Mag();
84  theta=photon_momentum.Theta();
85  phi=photon_momentum.Phi();
86 
87  if (emc_component==1)
88  {
89  // Z is calculated at R=100 cm
90  z=100*cos(theta);
91  } else
92  {
93  // X and Y are calculated at Z=100 cm
94  x=100*sin(theta)*cos(phi);
95  y=100*sin(theta)*sin(phi);
96  }
97 
98  max_energy=0;
99 
100  for (Int_t i=0; i<cluster_array->GetEntriesFast(); i++)
101  {
102  PndEmcCluster *cluster=(PndEmcCluster*)cluster_array->At(i);
103  cluster_energy=cluster->energy();
104  if (cluster_energy>max_energy)
105  {
106  max_energy=cluster_energy;
107  TVector3 cluster_pos=cluster->where();
108  cluster_phi=cluster_pos.Phi();
109  cluster_theta=cluster_pos.Theta();
110  if (emc_component==1)
111  {
112  cluster_z=100*cos(cluster_theta);
113  } else {
114  cluster_x=100*sin(cluster_theta)*cos(cluster_phi);
115  cluster_y=100*sin(cluster_theta)*sin(cluster_phi);
116  }
117  }
118  }
119 
120  e_diff=energy-max_energy;
121  h2energy->Fill(energy,e_diff);
122 
123  if (emc_component==1)
124  {
125  phi_diff=(cluster_phi-phi)*180./TMath::Pi();
126  h2phi->Fill(energy,phi_diff);
127  z_diff=cluster_z-z;
128  h2z->Fill(energy,z_diff);
129  } else {
130  x_diff=cluster_x-x;
131  h2x->Fill(energy,x_diff);
132  y_diff=cluster_y-y;
133  h2y->Fill(energy,y_diff);
134  }
135  }
136 
137  for (int i=0; i<h2energy->GetNbinsX();i++)
138  {
139  h_energy = h2energy->ProjectionY("hp_energy_1",i,i);
140  sigma_e=h_energy->GetRMS();
141  sigma_e_error=h_energy->GetRMSError();
142  energy=h2energy->GetXaxis()->GetBinCenter(i);
143  if (energy!=0) {
144  sigma_e_e=sigma_e/energy;
145  sigma_e_e_error=sigma_e_error/energy;
146  }
147  else {
148  sigma_e_e=0;
149  sigma_e_e_error=0;
150  }
151  h1energy_1->SetBinContent(i,sigma_e_e);
152  h1energy_1->SetBinError(i,sigma_e_e_error);
153  }
154 
155  for (int i=0; i<h2phi->GetNbinsX();i++)
156  {
157  h_phi = h2phi->ProjectionY("hp_phi",i,i);
158  sigma_phi=h_phi->GetRMS()*TMath::Pi()/180.;
159  sigma_phi_error=h_phi->GetRMSError()*TMath::Pi()/180.;
160  h1phi->SetBinContent(i,sigma_phi);
161  h1phi->SetBinError(i,sigma_phi_error);
162  }
163  for (int i=0; i<h2z->GetNbinsX();i++)
164  {
165  h_z = h2z->ProjectionY("hp_z",i,i);
166  sigma_z=h_z->GetRMS();
167  sigma_z_error=h_z->GetRMSError();
168  h1z->SetBinContent(i,sigma_z);
169  h1z->SetBinError(i,sigma_z_error);
170  }
171  std::cout<<"emc_component "<<emc_component<<" done."<<std::endl;
172  } else {
173  std::cout<<"EMC data file for component "<<emc_component<<" for gemetry "<<emcGeometry<<" does not exist"<<std::endl;
174  }
175 
179  emc_component=2;
180  TString outFile="emc_data_";
181  outFile+=emcGeometry;
182  outFile+="_";
183  outFile+=emc_component;
184  outFile+=".root";
185 
186  TFile* file2 = new TFile(outFile); //file you want to analyse
187  if (!file2->IsZombie())
188  {
189  fwcap_flag=true;
190  h2phi->Reset();h2x->Reset();h2y->Reset();h2z->Reset();h2energy->Reset();
191  TTree *t=(TTree *) file2->Get("pndsim") ;
192  TClonesArray* cluster_array=new TClonesArray("PndEmcCluster");
193  t->SetBranchAddress("EmcCluster",&cluster_array);
194  TClonesArray* mctrack_array=new TClonesArray("PndMCTrack");
195  t->SetBranchAddress("MCTrack",&mctrack_array);
196 
197  // Cluster energy
198  for (Int_t j=0; j< t->GetEntriesFast(); j++)
199  {
200  t->GetEntry(j);
201 
202  PndMCTrack *mctrack=(PndMCTrack *) mctrack_array->At(0);
203  photon_momentum=mctrack->GetMomentum();
204  energy=photon_momentum.Mag();
205  theta=photon_momentum.Theta();
206  phi=photon_momentum.Phi();
207 
208  if (emc_component==1)
209  {
210  // Z is calculated at R=100 cm
211  z=100*cos(theta);
212  } else
213  {
214  // X and Y are calculated at Z=100 cm
215  x=100*sin(theta)*cos(phi);
216  y=100*sin(theta)*sin(phi);
217  }
218 
219  max_energy=0;
220 
221  for (Int_t i=0; i<cluster_array->GetEntriesFast(); i++)
222  {
223  PndEmcCluster *cluster=(PndEmcCluster*)cluster_array->At(i);
224  cluster_energy=cluster->energy();
225  if (cluster_energy>max_energy)
226  {
227  max_energy=cluster_energy;
228  TVector3 cluster_pos=cluster->where();
229  cluster_phi=cluster_pos.Phi();
230  cluster_theta=cluster_pos.Theta();
231  if (emc_component==1)
232  {
233  cluster_z=100*cos(cluster_theta);
234  } else {
235  cluster_x=100*sin(cluster_theta)*cos(cluster_phi);
236  cluster_y=100*sin(cluster_theta)*sin(cluster_phi);
237  }
238  }
239  }
240 
241  e_diff=energy-max_energy;
242  h2energy->Fill(energy,e_diff);
243 
244  if (emc_component==1)
245  {
246  phi_diff=(cluster_phi-phi)*180./TMath::Pi();
247  h2phi->Fill(energy,phi_diff);
248  z_diff=cluster_z-z;
249  h2z->Fill(energy,z_diff);
250  } else {
251  x_diff=cluster_x-x;
252  h2x->Fill(energy,x_diff);
253  y_diff=cluster_y-y;
254  h2y->Fill(energy,y_diff);
255  }
256  }
257 
258 
259  for (int i=0; i<h2energy->GetNbinsX();i++)
260  {
261  h_energy = h2energy->ProjectionY("hp_energy_2",i,i);
262  sigma_e=h_energy->GetRMS();
263  sigma_e_error=h_energy->GetRMSError();
264  energy=h2energy->GetXaxis()->GetBinCenter(i);
265  if (energy!=0) {
266  sigma_e_e=sigma_e/energy;
267  sigma_e_e_error=sigma_e_error/energy;
268  }
269  else {
270  sigma_e_e=0;
271  sigma_e_e_error=0;
272  }
273  h1energy_2->SetBinContent(i,sigma_e_e);
274  h1energy_2->SetBinError(i,sigma_e_e_error);
275  }
276 
277  for (int i=0; i<h2x->GetNbinsX();i++)
278  {
279  h_x = h2x->ProjectionY("hp_x_2",i,i);
280  sigma_x=h_x->GetRMS();
281  sigma_x_error=h_x->GetRMSError();
282  h1x_2->SetBinContent(i,sigma_x);
283  h1x_2->SetBinError(i,sigma_x_error);
284  }
285  for (int i=0; i<h2y->GetNbinsX();i++)
286  {
287  h_y = h2y->ProjectionY("hp_y_2",i,i);
288  sigma_y=h_y->GetRMS();
289  sigma_y_error=h_y->GetRMSError();
290  h1y_2->SetBinContent(i,sigma_y);
291  h1y_2->SetBinError(i,sigma_y_error);
292  }
293 
294  std::cout<<"emc_component "<<emc_component<<" done."<<std::endl;
295  } else {
296  std::cout<<"EMC data file for component "<<emc_component<<" for gemetry "<<emcGeometry<<" does not exist"<<std::endl;
297  }
298 
302  emc_component=3;
303  TString outFile="emc_data_";
304  outFile+=emcGeometry;
305  outFile+="_";
306  outFile+=emc_component;
307  outFile+=".root";
308 
309  TFile* file3 = new TFile(outFile); //file you want to analyse
310  if (!file3->IsZombie())
311  {
312  bwcap_flag=true;
313  h2phi->Reset();h2x->Reset();h2y->Reset();h2z->Reset();h2energy->Reset();
314  TTree *t=(TTree *) file3->Get("pndsim") ;
315  TClonesArray* cluster_array=new TClonesArray("PndEmcCluster");
316  t->SetBranchAddress("EmcCluster",&cluster_array);
317  TClonesArray* mctrack_array=new TClonesArray("PndMCTrack");
318  t->SetBranchAddress("MCTrack",&mctrack_array);
319 
320  // Cluster energy
321  for (Int_t j=0; j< t->GetEntriesFast(); j++)
322  {
323  t->GetEntry(j);
324 
325  PndMCTrack *mctrack=(PndMCTrack *) mctrack_array->At(0);
326  photon_momentum=mctrack->GetMomentum();
327  energy=photon_momentum.Mag();
328  theta=photon_momentum.Theta();
329  phi=photon_momentum.Phi();
330 
331  if (emc_component==1)
332  {
333  // Z is calculated at R=100 cm
334  z=100*cos(theta);
335  } else
336  {
337  // X and Y are calculated at Z=100 cm
338  x=100*sin(theta)*cos(phi);
339  y=100*sin(theta)*sin(phi);
340  }
341 
342  max_energy=0;
343 
344  for (Int_t i=0; i<cluster_array->GetEntriesFast(); i++)
345  {
346  PndEmcCluster *cluster=(PndEmcCluster*)cluster_array->At(i);
347  cluster_energy=cluster->energy();
348  if (cluster_energy>max_energy)
349  {
350  max_energy=cluster_energy;
351  TVector3 cluster_pos=cluster->where();
352  cluster_phi=cluster_pos.Phi();
353  cluster_theta=cluster_pos.Theta();
354  if (emc_component==1)
355  {
356  cluster_z=100*cos(cluster_theta);
357  } else {
358  cluster_x=100*sin(cluster_theta)*cos(cluster_phi);
359  cluster_y=100*sin(cluster_theta)*sin(cluster_phi);
360  }
361  }
362  }
363 
364  e_diff=energy-max_energy;
365  h2energy->Fill(energy,e_diff);
366 
367  if (emc_component==1)
368  {
369  phi_diff=(cluster_phi-phi)*180./TMath::Pi();
370  h2phi->Fill(energy,phi_diff);
371  z_diff=cluster_z-z;
372  h2z->Fill(energy,z_diff);
373  } else {
374  x_diff=cluster_x-x;
375  h2x->Fill(energy,x_diff);
376  y_diff=cluster_y-y;
377  h2y->Fill(energy,y_diff);
378  }
379  }
380 
381 
382  for (int i=0; i<h2energy->GetNbinsX();i++)
383  {
384  h_energy = h2energy->ProjectionY("hp_energy_3",i,i);
385  sigma_e=h_energy->GetRMS();
386  sigma_e_error=h_energy->GetRMSError();
387  energy=h2energy->GetXaxis()->GetBinCenter(i);
388  if (energy!=0) {
389  sigma_e_e=sigma_e/energy;
390  sigma_e_e_error=sigma_e_error/energy;
391  }
392  else {
393  sigma_e_e=0;
394  sigma_e_e_error=0;
395  }
396  h1energy_3->SetBinContent(i,sigma_e_e);
397  h1energy_3->SetBinError(i,sigma_e_e_error);
398  }
399 
400  for (int i=0; i<h2x->GetNbinsX();i++)
401  {
402  h_x = h2x->ProjectionY("hp_x_3",i,i);
403  sigma_x=h_x->GetRMS();
404  sigma_x_error=h_x->GetRMSError();
405  h1x_3->SetBinContent(i,sigma_x);
406  h1x_3->SetBinError(i,sigma_x_error);
407  }
408  for (int i=0; i<h2y->GetNbinsX();i++)
409  {
410  h_y = h2y->ProjectionY("hp_y_3",i,i);
411  sigma_y=h_y->GetRMS();
412  sigma_y_error=h_y->GetRMSError();
413  h1y_3->SetBinContent(i,sigma_y);
414  h1y_3->SetBinError(i,sigma_y_error);
415  }
416  std::cout<<"emc_component "<<emc_component<<" done."<<std::endl;
417  } else {
418  std::cout<<"EMC data file for component "<<emc_component<<" for gemetry "<<emcGeometry<<" does not exist"<<std::endl;
419  }
420 
424  emc_component=4;
425  TString outFile="emc_data_";
426  outFile+=emcGeometry;
427  outFile+="_";
428  outFile+=emc_component;
429  outFile+=".root";
430 
431  TFile* file4 = new TFile(outFile); //file you want to analyse
432  if (!file4->IsZombie())
433  {
434  fsc_flag=true;
435  h2phi->Reset();h2x->Reset();h2y->Reset();h2z->Reset();h2energy->Reset();
436  TTree *t=(TTree *) file4->Get("pndsim") ;
437  TClonesArray* cluster_array=new TClonesArray("PndEmcCluster");
438  t->SetBranchAddress("EmcCluster",&cluster_array);
439  TClonesArray* mctrack_array=new TClonesArray("PndMCTrack");
440  t->SetBranchAddress("MCTrack",&mctrack_array);
441 
442  // Cluster energy
443  for (Int_t j=0; j< t->GetEntriesFast(); j++)
444  {
445  t->GetEntry(j);
446 
447  PndMCTrack *mctrack=(PndMCTrack *) mctrack_array->At(0);
448  photon_momentum=mctrack->GetMomentum();
449  energy=photon_momentum.Mag();
450  theta=photon_momentum.Theta();
451  phi=photon_momentum.Phi();
452 
453  if (emc_component==1)
454  {
455  // Z is calculated at R=100 cm
456  z=100*cos(theta);
457  } else
458  {
459  // X and Y are calculated at Z=100 cm
460  x=100*sin(theta)*cos(phi);
461  y=100*sin(theta)*sin(phi);
462  }
463 
464  max_energy=0;
465 
466  for (Int_t i=0; i<cluster_array->GetEntriesFast(); i++)
467  {
468  PndEmcCluster *cluster=(PndEmcCluster*)cluster_array->At(i);
469  cluster_energy=cluster->energy();
470  if (cluster_energy>max_energy)
471  {
472  max_energy=cluster_energy;
473  TVector3 cluster_pos=cluster->where();
474  cluster_phi=cluster_pos.Phi();
475  cluster_theta=cluster_pos.Theta();
476  if (emc_component==1)
477  {
478  cluster_z=100*cos(cluster_theta);
479  } else {
480  cluster_x=100*sin(cluster_theta)*cos(cluster_phi);
481  cluster_y=100*sin(cluster_theta)*sin(cluster_phi);
482  }
483  }
484  }
485 
486  e_diff=energy-max_energy;
487  h2energy->Fill(energy,e_diff);
488 
489  if (emc_component==1)
490  {
491  phi_diff=(cluster_phi-phi)*180./TMath::Pi();
492  h2phi->Fill(energy,phi_diff);
493  z_diff=cluster_z-z;
494  h2z->Fill(energy,z_diff);
495  } else {
496  x_diff=cluster_x-x;
497  h2x->Fill(energy,x_diff);
498  y_diff=cluster_y-y;
499  h2y->Fill(energy,y_diff);
500  }
501  }
502 
503 
504  for (int i=0; i<h2energy->GetNbinsX();i++)
505  {
506  h_energy = h2energy->ProjectionY("hp_energy_4",i,i);
507  sigma_e=h_energy->GetRMS();
508  sigma_e_error=h_energy->GetRMSError();
509  energy=h2energy->GetXaxis()->GetBinCenter(i);
510  if (energy!=0) {
511  sigma_e_e=sigma_e/energy;
512  sigma_e_e_error=sigma_e_error/energy;
513  }
514  else {
515  sigma_e_e=0;
516  sigma_e_e_error=0;
517  }
518  h1energy_4->SetBinContent(i,sigma_e_e);
519  h1energy_4->SetBinError(i,sigma_e_e_error);
520  }
521 
522  for (int i=0; i<h2x->GetNbinsX();i++)
523  {
524  h_x = h2x->ProjectionY("hp_x_4",i,i);
525  sigma_x=h_x->GetRMS();
526  sigma_x_error=h_x->GetRMSError();
527  h1x_4->SetBinContent(i,sigma_x);
528  h1x_4->SetBinError(i,sigma_x_error);
529  }
530  for (int i=0; i<h2y->GetNbinsX();i++)
531  {
532  h_y = h2y->ProjectionY("hp_y_4",i,i);
533  sigma_y=h_y->GetRMS();
534  sigma_y_error=h_y->GetRMSError();
535  h1y_4->SetBinContent(i,sigma_y);
536  h1y_4->SetBinError(i,sigma_y_error);
537  }
538  std::cout<<"emc_component "<<emc_component<<" done."<<std::endl;
539  } else {
540  pars4[0]=0.202292;
541  pars4[1]=0.717711;
542  pars4[2]=0.083642;
543  pars4[3]=0.000000;
544  pars4[4]=0.333771;
545  pars4[5]=0.715637;
546  pars4[6]=0.341344;
547  pars4[7]=0.393149;
548  pars4[8]=0.397055;
549  pars4[9]=0.321306;
550  std::cout<<"EMC data file for component "<<emc_component<<" for gemetry "<<emcGeometry<<" does not exist"<<std::endl;
551  }
552 
556 
557 
559  if (barrel_flag){
560 
561  TCanvas* c1 = new TCanvas("c1", "Energy-position resolution vs Energy", 100, 100, 800, 800);
562  c1->Divide(2,2);
563  c1->cd(1);
564  h1energy_1->Draw("E1");
565 
566  TF1 *f1 = new TF1("f1", fit1, 0., 7.,4);
567  f1->SetLineColor(2);
568  h1energy_1->Fit("f1");
569  f1->GetParameters(pars1);
570 
571  c1->cd(2);
572  h1z->Draw("E1");
573  TF1 *f2 = new TF1("f2", fit2, 0., 7.,3);
574  f2->SetLineColor(2);
575  h1z->Fit("f2");
576  f2->GetParameters(pars1+4);
577 
578  c1->cd(3);
579  h1phi->Draw("E1");
580  TF1 *f3 = new TF1("f3", fit2, 0., 7.,3);
581  f3->SetLineColor(2);
582  h1phi->Fit("f3");
583  f3->GetParameters(pars1+7);
584  }
585 
587  if (fwcap_flag)
588  {
589  TCanvas* c2 = new TCanvas("c2", "Energy-position resolution vs Energy", 100, 100, 800, 800);
590  c2->Divide(2,2);
591  c2->cd(1);
592  h1energy_2->Draw("E1");
593 
594  TF1 *f4 = new TF1("f4", fit1, 0., 7.,4);
595  f4->SetLineColor(2);
596  h1energy_2->Fit("f4");
597  f4->GetParameters(pars2);
598 
599  c2->cd(2);
600  h1x_2->Draw("E1");
601  TF1 *f5 = new TF1("f5", fit2, 0., 7.,3);
602  f5->SetLineColor(2);
603  h1x_2->Fit("f5");
604  f5->GetParameters(pars2+4);
605 
606  c2->cd(3);
607  h1y_2->Draw("E1");
608  TF1 *f6 = new TF1("f6", fit2, 0., 7.,3);
609  f6->SetLineColor(2);
610  h1y_2->Fit("f6");
611  f6->GetParameters(pars2+7);
612  }
613 
615  if (bwcap_flag)
616  {
617  TCanvas* c3 = new TCanvas("c3", "Energy-position resolution vs Energy", 100, 100, 800, 800);
618  c3->Divide(2,2);
619  c3->cd(1);
620  h1energy_3->Draw("E1");
621 
622  TF1 *f7 = new TF1("f7", fit1, 0., 7.,4);
623  f7->SetLineColor(2);
624  h1energy_3->Fit("f7");
625  f7->GetParameters(pars3);
626 
627  c3->cd(2);
628  h1x_3->Draw("E1");
629  TF1 *f8 = new TF1("f8", fit2, 0., 7.,3);
630  f8->SetLineColor(2);
631  h1x_3->Fit("f8");
632  f8->GetParameters(pars3+4);
633 
634  c3->cd(3);
635  h1y_3->Draw("E1");
636  TF1 *f9 = new TF1("f9", fit2, 0., 7.,3);
637  f9->SetLineColor(2);
638  h1y_3->Fit("f9");
639  f9->GetParameters(pars3+7);
640  }
641 
643  if (fsc_flag)
644  {
645  TCanvas* c4 = new TCanvas("c4", "Energy-position resolution vs Energy", 100, 100, 800, 800);
646  c4->Divide(2,2);
647  c4->cd(1);
648  h1energy_4->Draw("E1");
649 
650  TF1 *f10 = new TF1("f10", fit1, 0., 7.,4);
651  f10->SetLineColor(2);
652  h1energy_4->Fit("f10");
653  f10->GetParameters(pars4);
654 
655  c4->cd(2);
656  h1x_4->Draw("E1");
657  TF1 *f11 = new TF1("f11", fit2, 0., 7.,3);
658  f11->SetLineColor(2);
659  h1x_4->Fit("f11");
660  f11->GetParameters(pars4+4);
661 
662  c4->cd(3);
663  h1y_4->Draw("E1");
664  TF1 *f12 = new TF1("f12", fit2, 0., 7.,3);
665  f12->SetLineColor(2);
666  h1y_4->Fit("f12");
667  f12->GetParameters(pars4+7);
668  }
669 
671  TString fileName="emc_error_matrix_";
672  fileName+=fileVersion;
673  fileName+=".root";
674  TFile *errorfile = new TFile(fileName, "RECREATE");
675 
676  PndEmcErrorMatrixParObject *parObject=new PndEmcErrorMatrixParObject();
677  enum {barrel, fwcap, bwcap, fsc};
678 
679  parObject->SetErrorMatrix(barrel, pars1);
680  parObject->SetErrorMatrix(fwcap, pars2);
681  parObject->SetErrorMatrix(bwcap, pars3);
682  parObject->SetErrorMatrix(fsc, pars4);
683 
684  parObject->Write();
685  errorfile->Close();
686 
687  return 0;
688 }
mctrack
PndMCTrack * mctrack
Definition: analysis_point_hit_fwendcap.C:250
phi_diff
double phi_diff
Definition: macro/detectors/emc/dedicated/reco_analys.C:25
theta
double theta
Definition: macro/detectors/emc/dedicated/reco_analys.C:24
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
PndEmcCluster::where
TVector3 where() const
Definition: PndEmcCluster.cxx:111
basiclibs
basiclibs()
f10
TF1 * f10
Definition: plot_acceptance_stt.C:259
c4
c4
Definition: plot_dirc.C:71
i
Int_t i
Definition: run_full.C:25
photon_momentum
TVector3 photon_momentum
Definition: anal_hit_digi_cluster_fwendcap.C:45
f1
TF1 * f1
Definition: reco_analys2.C:50
outFile
TString outFile
Definition: hit_dirc.C:17
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
PndEmcErrorMatrixParObject
Container class for EMC error matrix parameter class is inherited from FairParGenericSet.
Definition: PndEmcErrorMatrixPar.h:38
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
c2
c2
Definition: plot_dirc.C:39
f4
TFile * f4
Definition: ReadMesDecCorrRes.C:9
h2phi
TH2F * h2phi
Definition: macro/detectors/emc/dedicated/reco_analys.C:33
cluster_theta
Double_t cluster_theta
Definition: anal_hit_digi_cluster_fwendcap.C:31
f3
TFile * f3
Definition: ReadMesDecCorrRes.C:8
cluster_array
TClonesArray * cluster_array
Definition: anal_hit_digi_cluster_fwendcap.C:27
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
cluster_energy
Double_t cluster_energy
Definition: anal_hit_digi_cluster_fwendcap.C:31
fit2
Double_t fit2(Double_t *x, Double_t *par)
Definition: error_matrix_fit.C:706
cluster_phi
Double_t cluster_phi
Definition: anal_hit_digi_cluster_fwendcap.C:31
z
Double_t z
Definition: createRootFscGeometryFile.C:146
rootlogon
rootlogon()
Definition: outdated/mpiTools/macros/emc/rootlogon.C:1
PndEmcCluster
a cluster (group of neighboring crystals) of hit emc crystals
Definition: PndEmcCluster.h:29
c1
c1
Definition: plot_dirc.C:35
c3
c3
Definition: plot_dirc.C:50
max_energy
double max_energy
Definition: macro/detectors/emc/dedicated/reco_analys.C:27
phi
double phi
Definition: macro/detectors/emc/dedicated/reco_analys.C:24
x
Double_t x
Definition: createRootFscGeometryFile.C:159
f2
TFile * f2
Definition: ReadMesDecCorrRes.C:5
fit1
Double_t fit1(Double_t *x, Double_t *par)
Definition: error_matrix_fit.C:690
PndEmcCluster::energy
virtual Double_t energy() const
Definition: PndEmcCluster.cxx:79
t
TTree * t
Definition: bump_analys.C:13
y
Double_t y
Definition: createRootFscGeometryFile.C:159
mctrack_array
TClonesArray * mctrack_array
Definition: anal_hit_digi_cluster_fwendcap.C:11
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
PndEmcErrorMatrixParObject::SetErrorMatrix
void SetErrorMatrix(Int_t detectorComponent, Double_t *pars)
Definition: PndEmcErrorMatrixPar.cxx:22
energy
Double_t energy
Definition: plot_dirc.C:15
Double_t fit1 ( Double_t *  x,
Double_t *  par 
)

Definition at line 690 of file error_matrix_fit.C.

References a, c, and res.

Referenced by error_matrix_fit().

691 {
692  double xx=x[0];
693 
694  double a=par[0];
695  double power=par[1];
696  double c=par[2];
697  double q=par[3];
698 
699  double res;
700 
701  res=a*a/pow(xx,power)+c*c+q*q/(xx*xx);
702 
703  return res;
704 }
res
Int_t res
Definition: anadigi.C:166
par
Double_t par[3]
Definition: anaLmdCluster.C:146
a
Int_t a
Definition: anaLmdDigi.C:126
c
GeV c
Definition: ReadMesDecCorrRes_PiT.C:28
x
Double_t x
Definition: createRootFscGeometryFile.C:159
Double_t fit2 ( Double_t *  x,
Double_t *  par 
)

Definition at line 706 of file error_matrix_fit.C.

References a, c, and res.

Referenced by error_matrix_fit().

707 {
708  double xx=x[0];
709 
710  double a=par[0];
711  double power=par[1];
712  double c=par[2];
713 
714  double res;
715 
716  res=a*a/pow(xx,power)+c*c;
717 
718  return res;
719 }
res
Int_t res
Definition: anadigi.C:166
par
Double_t par[3]
Definition: anaLmdCluster.C:146
a
Int_t a
Definition: anaLmdDigi.C:126
c
GeV c
Definition: ReadMesDecCorrRes_PiT.C:28
x
Double_t x
Definition: createRootFscGeometryFile.C:159
Generated on Wed Apr 3 2019 10:03:13 for FairRoot/PandaRoot by   doxygen 1.8.5