アットウィキロゴ
おかめも
アットウィキロゴ
おかめも
ページ一覧

卒研めも > Calibration > Co-60 > macro0111

Co用だがいまのところNa-22とほとんど変わんない状態。コバルトの2つのピークはうまくフィッティングできないので、ダブルガウシアンによるフィッティングに変える予定。 

  1. #include<fstream>
  2. #include<TH1.h>
  3. #include<TH2.h>
  4. #include<TFile.h>
  5. #include<TTree.h>
  6. #include<TCanvas.h>
  7. #include<parameter_co.h>
  8. #include<math.h>
  9. //#include<iostream>
  10.  
  11. #include<iostream> using namespace std;
  12.  
  13. //ここまでで使うライブラリをincludeしてます．parameter.hというのは自分でつくってヘッダファイルで，キャリブレーション用のパラメターが書いてあります．
  14.  
  15.  
  16. void macro_calb_co(){  //Last modified 2012/01/11
  17.   //difference from macro_calb_oka:
  18.   //fit by Gaussian but not draw it,
  19.   //and fit again by "Gaussian + linear" and draw it
  20.  
  21.   char filename[99];//root file name
  22.  
  23.   sprintf(filename,"60co-20120111.001.root");  
  24.  
  25.  
  26.   TFile *f1=new TFile(filename);//opening root file
  27. //ファイルを開いています。
  28.  
  29.   Double_t p[9][9][9];
  30. //pというDouble_t クラスの配列は、フィッティングのときに使うパラメターです。p[i][j][k]はチャネルiのj番目のピークのk番目のパラメターを意味しています。j=1は511keVのピーク、j=2は1275のピーク、j=3はサムピークです。k=0はガウシアン全体にかかる定数（使わない）k=1はMean，k=2はSigmaです。
  31.  
  32.   TH1F *h[99];
  33.   h[1]=h101;
  34.   h[2]=h102;
  35.   h[3]=h103;
  36.   h[4]=h104;
  37.   h[5]=h105;
  38.   h[6]=h106;
  39. //TH1Fは1次のFloatのHistgramのクラスです。ROOTファイルに入っているh101～h106までのヒストグラムをプログラム中で扱いやすいようにh[i]という配列を宣言して、そこに代入しています。
  40.  
  41.  
  42.   // end of front matter
  43.  
  44.  
  45.   //gStyle->SetOptFit(0001);
  46.   gStyle->SetOptFit(0000);
  47. //ここはフィッティングのパラメータをどのようにキャンバスに出力するか決めています  
  48.  
  49.  
  50.  
  51.   TCanvas *c1 = new TCanvas("c1","Raw Data and Fitting");
  52.   c1->Divide(2,3);
  53.  
  54.  
  55.   TCanvas *c4 = new TCanvas("c4","Fitting with Gaussian and Linear function");
  56.   c4->Divide(2,3);
  57.  
  58.  
  59.  
  60. //TCanvasというクラスはヒストグラムとかを出力するためのクラスです．キャンバスをDivideするとキャンバスを複数に分けることができます
  61.  
  62.  
  63.   for(Int_t ch=1;ch<=6;ch++){
  64.  
  65.     char hclone[256]; 
  66.     sprintf(hclone,"hclone%d",ch);
  67.     TH1F *hoge=h[ch]->Clone(hclone);
  68.     h[ch+10]=hoge;
  69.  
  70.     sprintf(hclone,"hclone2%d",ch);
  71.     TH1F *hoge2=h[ch]->Clone(hclone);
  72.     h[ch+20]=hoge2;
  73.     //h101-h106を違う目的で使いたいので複製しています
  74.   }
  75.  
  76.  
  77.   for(Int_t ch=1;ch<=6;ch++){
  78.     /*
  79.            Canvas No. 1
  80.      */
  81.  
  82.     c1->cd(ch)->SetLogy();    
  83.     h[ch]->Draw();
  84.  
  85.     //Preparation of Fitting function
  86.     //ここからはピークをフィッティングするための関数を準備しています．
  87.  
  88.     TF1 *fit1;    
  89.     char fname[256];
  90.     char mean1[256],mean2[256],mean3[256];
  91.     char resol1[256],resol2[256],resol3[256];
  92.  
  93.     Double_t sigmae1[3];  
  94.  
  95.     sprintf(fname,"f[%d]",ch);
  96.  
  97.     fit1 = new TF1(fname,"gaus",p1start[ch],p1end[ch]);
  98.  
  99.     h[ch]->Fit(fit1,"","",p1start[ch],p1end[ch]); 
  100.     //h[ch]にfit1をフィッティングします
  101.     fit1->GetParameters(p[ch][1]);
  102.     sigmae1[0]=fit1->GetParError(2);
  103.     //フィットしたfit1について，そのフィッティングパラメータを得ます．これによってp[ch][1][1]にはMeanが，p[ch][1][2]にはSigmaが入ります．
  104.     sprintf(mean1," Mean =%lf",p[ch][1][1]);
  105.     sprintf(resol1,"Sigma =%lf",p[ch][1][2]);
  106.     //凡例に表示するための文字列を準備しています
  107.  
  108.  
  109.     TF1 *fit2;
  110.     sprintf(fname,"f[1%d]",ch);
  111.     fit2 = new TF1(fname,"gaus",p2start[ch],p2end[ch]);
  112.     h[ch]->Fit(fit2,"+","",p2start[ch],p2end[ch]); 
  113.     fit2->GetParameters(p[ch][2]);
  114.     sigmae1[1]=fit2->GetParError(2);    
  115.     sprintf(mean2," Mean =%lf",p[ch][2][1]);
  116.     sprintf(resol2,"Sigma =%lf",p[ch][2][2]);
  117.  
  118.     TF1 *fit3;
  119.     sprintf(fname,"f[1%d]",ch);
  120.     fit3 = new TF1(fname,"gaus",p3start[ch],p3end[ch]);
  121.     h[ch]->Fit(fit3,"+","",p3start[ch],p3end[ch]); 
  122.     fit3->GetParameters(p[ch][3]);
  123.     sigmae1[2]=fit3->GetParError(2);    
  124.     sprintf(mean3," Mean =%lf",p[ch][3][1]);
  125.     sprintf(resol3,"Sigma =%lf",p[ch][3][2]);
  126.  
  127.  
  128.     h[ch]->Draw();//もう１回Drawします．なぜかこうしないとうまくいかなかったきがする
  129.  
  130.  
  131.     //ここからは凡例(Legend)の設定です    
  132.     TLegend *l = new TLegend(0.6255,0.75,0.875,0.45);
  133.     //TLegendクラスの変数 l を宣言しました
  134.     l->AddEntry(fit1,"511 keV","lp");
  135.     //lにfit1を追加し，その表示名を511 keVとしました
  136.     l->AddEntry((TObject*)0, mean1,"");
  137.     //文字列mean1をlに追加しました
  138.     l->AddEntry((TObject*)0, resol1,"");
  139.     l->AddEntry(fit2,"1274 keV","lp");
  140.     l->AddEntry((TObject*)0, mean2,"");  
  141.     l->AddEntry((TObject*)0, resol2,"");  
  142.     l->AddEntry(fit2,"1274 keV","lp");
  143.     l->AddEntry((TObject*)0, mean3,"");  
  144.     l->AddEntry((TObject*)0, resol3,"");  
  145.  
  146.     l->SetFillColor(0);// l の背景を白にします
  147.     l->Draw();//lをキャンバスに描きます
  148.  
  149.  
  150.  
  151.  
  152.     //1次式＋ガウシアンでフィットし直す
  153.  
  154.     char fglname[99];
  155.     Double_t paratest[5];
  156.     Double_t ex[3];
  157.  
  158.     Double_t x[3],sigma2[3],sigmae2[3];
  159.  
  160.     /*
  161.                Canvas No.4
  162.  
  163.      */
  164.  
  165.  
  166.  
  167.     c4->cd(ch)->SetLogy();
  168.  
  169.     h[ch+20]->Draw();
  170.     TLegend *l3 = new TLegend(0.6255,0.75,0.875,0.45);
  171.  
  172.     TF1 *fit4;
  173.     sprintf(fglname,"fgl1%d",ch);
  174.     fit4 = new TF1(fglname,"[0]/sqrt(2.0*3.14)/[1]*exp(-(x-[2])*(x-[2])/2.0/[1]/[1])+[3]*x+[4]");//[0]:constant ,[1]:sigma, [2]:mean
  175.     fit4->SetParameters(p[ch][1][0],p[ch][1][2],p[ch][1][1],0,0);//パラメータの初期値の設定
  176.     h[ch+20]->Fit(fglname,"","",p1start[ch],p1end[ch]);
  177.     fit4 = h[ch+20]->GetFunction(fglname);
  178.     fit4->GetParameters(paratest);// parameterをgetしてこれを使えば直線の式が
  179.     x[0]=paratest[2];
  180.     ex[0]=fit4->GetParError(2);
  181.     sigma2[0]=paratest[1];
  182.     sigmae2[0]=fit4->GetParError(1);
  183.  
  184.  
  185.     TF1 *flin1= new TF1("flin1","[0]*x+[1]",p1start[ch],p1end[ch]);
  186.     flin1->FixParameter(0,paratest[3]);
  187.     flin1->FixParameter(1,paratest[4]);//"fix" is proper, not "set"
  188.     h[ch+20]->Fit(flin1,"Q0+","",p1start[ch],p1end[ch]);
  189.     flin1=h[ch+20]->GetFunction("flin1");
  190.  
  191.     char caption[32]; 
  192.     l3->AddEntry(fit4,"511 keV","lp");
  193.     sprintf(caption,"mean=%lf",paratest[2]);
  194.     l3->AddEntry((TObject*)0, caption,"");
  195.     sprintf(caption,"sigma=%lf",paratest[1]);
  196.     l3->AddEntry((TObject*)0, caption,"");
  197.  
  198.  
  199.     TF1 *fit5;
  200.     sprintf(fglname,"fgl2%d",ch);
  201.     fit5 = new TF1(fglname,"[0]/sqrt(2.0*3.14)/[1]*exp(-(x-[2])*(x-[2])/2.0/[1]/[1])+[3]*x+[4]");//[0]:constant ,[1]:sigma, [2]:mean
  202.     fit5->SetParameters(p[ch][2][0],p[ch][2][2],p[ch][2][1],0,0);
  203.     h[ch+20]->Fit(fglname,"0+","",p2start[ch],p2end[ch]);
  204.     fit5 = h[ch+20]->GetFunction(fglname);
  205.     fit5->GetParameters(paratest);
  206.     x[1]=paratest[2];
  207.     ex[1]=fit5->GetParError(2);
  208.     sigma2[1]=paratest[1];
  209.     sigmae2[1]=fit5->GetParError(1);
  210.  
  211.     TF1 *flin2= new TF1("flin2","[0]*x+[1]",p2start[ch],p2end[ch]);
  212.     flin2->FixParameter(0,paratest[3]);
  213.     flin2->FixParameter(1,paratest[4]);//"fix" is proper, not "set"
  214.     h[ch+20]->Fit(flin2,"Q0+","",p2start[ch],p2end[ch]);
  215.     flin2=h[ch+20]->GetFunction("flin2");
  216.  
  217.     char caption[32]; 
  218.     l3->AddEntry(fit5,"1274.5 keV","lp");
  219.     sprintf(caption,"mean=%lf",paratest[2]);
  220.     l3->AddEntry((TObject*)0, caption,"");
  221.     sprintf(caption,"sigma=%lf",paratest[1]);
  222.     l3->AddEntry((TObject*)0, caption,"");
  223.  
  224.  
  225.  
  226.     TF1 *fit6;	 	 
  227.     sprintf(fglname,"fgl3%d",ch);
  228.     fit6 = new TF1(fglname,"[0]/sqrt(2.0*3.14)/[1]*exp(-(x-[2])*(x-[2])/2.0/[1]/[1])+[3]*x+[4]");//[0]:constant ,[1]:sigma, [2]:mean
  229.     fit6->SetParameters(p[ch][3][0],p[ch][3][2],p[ch][3][1],0,0);
  230.     h[ch+20]->Fit(fglname,"0+","",p3start[ch],p3end[ch]);
  231.     fit6 = h[ch+20]->GetFunction(fglname);
  232.     fit6->GetParameters(paratest);
  233.     x[2]=paratest[2];
  234.     ex[2]=fit6->GetParError(2);
  235.     sigma2[2]=paratest[1];
  236.     sigmae2[2]=fit6->GetParError(1);
  237.     TF1 *flin3= new TF1("flin3","[0]*x+[1]",p3start[ch],p3end[ch]);
  238.     flin3->FixParameter(0,paratest[3]);
  239.     flin3->FixParameter(1,paratest[4]);//"fix" is proper, not "set"
  240.     h[ch+20]->Fit(flin3,"Q0+","",p3start[ch],p3end[ch]);
  241.     flin3=h[ch+20]->GetFunction("flin3");
  242.  
  243.     char caption[32]; 
  244.     l3->AddEntry(fit6,"511+1274 keV","lp");
  245.     sprintf(caption,"mean=%lf",paratest[2]);
  246.     l3->AddEntry((TObject*)0, caption,"");
  247.     sprintf(caption,"sigma=%lf",paratest[1]);
  248.     l3->AddEntry((TObject*)0, caption,"");
  249.  
  250.  
  251.  
  252.     h[ch+20]->Draw("same");
  253.  
  254.     fit4->Draw("same");  
  255.     flin1->Draw("same");		 
  256.  
  257.     fit5->SetLineColor(3);
  258.     fit5->Draw("same");  
  259.     flin2->SetLineColor(3); 
  260.     flin2->Draw("same");		 
  261.  
  262.     fit6->SetLineColor(7);
  263.     fit6->Draw("same"); 
  264.     flin3->SetLineColor(7);
  265.     flin3->Draw("same");	
  266.  
  267.     l3->SetFillColor(0);// l の背景を白にします
  268.     l3->Draw();//lをキャンバスに描きます
  269.  
  270.  
  271.  
  272.  
  273.   }
  274.  
  275.  
  276.   /* 
  277.   //キャンバスをPDFに出力します
  278.   //Print
  279.   char pdfstart[99];
  280.   sprintf(pdfstart,"%s.pdf(",filename);
  281.   char pdfmiddle[99];
  282.   sprintf(pdfmiddle,"%s.pdf",filename);
  283.   char pdfend[99];
  284.   sprintf(pdfend,"%s.pdf)",filename);
  285.   
  286.   c1->Print(pdfstart);
  287.   c2->Print(pdfmiddle);
  288.   c3->Print(pdfmiddle);
  289.   c4->Print(pdfend);
  290.   */  
  291.  
  292.  
  293. }
  294.  
  295.  
  296.  
  297.  
「macro0111」をウィキ内検索
最終更新：2012年01月11日 18:31