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

卒研めも > Calibration > Na-22 > macro1221

課題
  • 誤差付きのプロット（165行目附近）。ただガウシアンのフィッティングによる誤差がかなり小さいのであまり関係ないかも
    • eyていうのはADC Channelの誤差なので実はこの場合だとx軸の誤差
  • ガウシアン＋線形項によるフィッティング


  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.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_oka(){  //Last modified 2011/11/22
  17.  
  18.   char filename[99];//root file name
  19.  
  20.   // sprintf(filname,"22na20111104.001.root");
  21.   //sprintf(filename,"22na-10.03.001-002.root");
  22.   //spintf(filename,"22na-11.10.20-21.root");
  23.   sprintf(filename,"22na-20111028.002.root");
  24.   //  sprintf(filename,"22na20111104.001.root");
  25.   //  sprintf(filename,"22na-20111111.001.root");
  26.   //  sprintf(filename,"22na-20111118.001.root");
  27.   //  sprintf(filename,"22na-20111201.001.root");  
  28. //filenameという変数を宣言して，そこにファイル名を代入してます．この変数はファイルを開くのと，あとでPDFファイルをつくるときに使います．sprintfはprintfと似た関数で，変数に文字列をいれることができます．ROOTのマクロを書くときにはけっこう便利です．あと、文字列をいれるときは，charというクラスでその文字数以上の配列として宣言するみたいです。
  29.  
  30.  
  31.   TFile *f1=new TFile(filename);//opening root file
  32. //ファイルを開いています。
  33.  
  34.   Double_t p[9][9][9];
  35. //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です。
  36.  
  37.   Double_t ey[3];
  38.  
  39.   TH1F *h[99];
  40.   h[1]=h101;
  41.   h[2]=h102;
  42.   h[3]=h103;
  43.   h[4]=h104;
  44.   h[5]=h105;
  45.   h[6]=h106;
  46. //TH1Fは1次のFloatのHistgramのクラスです。ROOTファイルに入っているh101～h106までのヒストグラムをプログラム中で扱いやすいようにh[i]という配列を宣言して、そこに代入しています。
  47.  
  48.  
  49.  
  50.   // end of front matter
  51.  
  52.  
  53.   //  gStyle->SetOptFit(0001);
  54.     gStyle->SetOptFit(0000);
  55. //ここはフィッティングのパラメータをどのようにキャンバスに出力するか決めています  
  56.  
  57.  
  58.  
  59.   TCanvas *c1 = new TCanvas("c1","Raw Data and Fitting");
  60.   c1->Divide(2,3);
  61.  
  62.   TCanvas *c2 = new TCanvas("c2","Enegy Calibration");
  63.   c2->Divide(2,3);
  64.  
  65.   TCanvas *c3 = new TCanvas("c3","Calibrated Hist");
  66.   c3->Divide(2,3);
  67.  
  68. //TCanvasというクラスはヒストグラムとかを出力するためのクラスです．キャンバスをDivideするとキャンバスを複数に分けることができます
  69.  
  70.  
  71.  
  72.   for(Int_t ch=1;ch<=6;ch++){
  73.  
  74.     /*    char hclone[256]; 
  75.     sprintf(hclone,"hclone%d",ch);
  76.     hoge = h[ch]->Clone(hclone);
  77.     h[ch+10]=hoge;
  78.     //h101〜h106を違う目的で使いたいので複製しています*/
  79.     c1->cd(ch)->SetLogy();    
  80.     h[ch]->Draw();
  81. //まず，キャンバスc1の6つに分割したうちの1番目に移動し，ch=1のときはh[1]=h101をDrawします．
  82.  
  83.     //Preparation of Fitting function
  84. //ここからはピークをフィッティングするための関数を準備しています．
  85.  
  86.     TF1 *fit1;    
  87.     char fname[256];
  88.     char mean1[256],mean2[256],mean3[256];
  89.     char resol1[256],resol2[256],resol3[256];
  90.  
  91.     sprintf(fname,"f[%d]",ch);
  92. //fnameは関数の名前で，ひとつめのピークにはチャネルごとにf[1]?f[6]をフィッティングします
  93.     fit1 = new TF1(fname,"gaus",p1start[ch],p1end[ch]);
  94. //この行ではフィッティングのための関数fit1をつくっています．名前はfnameに入っている値，種類はガウシアン，定義域はp1startとp1endに入っている値です
  95.     h[ch]->Fit(fit1,"","",p1start[ch],p1end[ch]); 
  96. //h[ch]にfit1をフィッティングします
  97.     fit1->GetParameters(p[ch][1]); 
  98. //フィットしたfit1について，そのフィッティングパラメータを得ます．これによってp[ch][1][1]にはMeanが，p[ch][1][2]にはSigmaが入ります．
  99.     ey[0] = fit1->GetParError(1); 
  100.     sprintf(mean1," Mean =%lf",p[ch][1][1]);
  101.     sprintf(resol1,"Resol.=%lf",2*sqrt(2*log(2.0))*p[ch][1][2]/p[ch][1][1]);
  102. //凡例に表示するための文字列を準備しています
  103.  
  104. //同じことを2つ目，3つ目のピークについても行います．
  105.     TF1 *fit2;
  106.     sprintf(fname,"f[1%d]",ch);
  107.     fit2 = new TF1(fname,"gaus",p2start[ch],p2end[ch]);
  108.     h[ch]->Fit(fit2,"","",p2start[ch],p2end[ch]); 
  109.     fit2->GetParameters(p[ch][2]);    
  110.     ey[1] = fit2->GetParError(1); 
  111.     sprintf(mean2," Mean =%lf",p[ch][2][1]);
  112.     sprintf(resol2,"Resol.=%lf",2*sqrt(2*log(2.0))*p[ch][2][2]/p[ch][2][1]);
  113.  
  114.     TF1 *fit3;
  115.     sprintf(fname,"f[10%d]",ch);
  116.     fit3 = new TF1(fname,"gaus",p3start[ch],p3end[ch]);
  117.     h[ch]->Fit(fname,"","",p3start[ch],p3end[ch]); 
  118.     fit3->GetParameters(p[ch][3]);    
  119.     ey[2] = fit3->GetParError(1);
  120.     sprintf(mean3," Mean =%lf",p[ch][3][1]);
  121.     sprintf(resol3,"Resol.=%lf",2*sqrt(2*log(2.0))*p[ch][3][2]/p[ch][3][1]);
  122.  
  123.     h[ch]->Draw();//もう１回Drawします．なぜかこうしないとうまくいかなかったきがする
  124.  
  125.     fit1->Draw("sames");//fit1をDrawします．"sames"と書くと上書きされます．
  126.  
  127.     fit2->SetLineColor(3);
  128.     fit2->Draw("sames");//fit2の線の色を設定したのち，上書きします
  129.  
  130.     fit3->SetLineColor(7);
  131.     fit3->Draw("sames");
  132.  
  133.  
  134. //ここからは凡例(Legend)の設定です    
  135.     TLegend *l = new TLegend(0.6255,0.75,0.875,0.45);
  136. //TLegendクラスの変数 l を宣言しました
  137.     l->AddEntry(fit1,"511 keV","lp");
  138.     //lにfit1を追加し，その表示名を511 keVとしました
  139.     l->AddEntry((TObject*)0, mean1,"");
  140. //文字列mean1をlに追加しました
  141.     l->AddEntry((TObject*)0, resol1,"");
  142.     l->AddEntry(fit2,"1274 keV","lp");
  143.     l->AddEntry((TObject*)0, mean2,"");  
  144.     l->AddEntry((TObject*)0, resol2,"");  
  145.     l->AddEntry(fit3,"511+1274 keV","lp");
  146.     l->AddEntry((TObject*)0, mean3,""); 
  147.     l->AddEntry((TObject*)0, resol3,"");   
  148.     l->SetFillColor(0);// l の背景を白にします
  149.     l->Draw();//lをキャンバスに描きます
  150.  
  151.     // Graph of calibration   
  152.     //エネルギーキャリブレーションの式を求めるグラフをつくります
  153.  
  154.     c2->cd(ch);
  155.  
  156. //xとyの散布図をつくるために，値を準備しておきます  
  157.     Double_t x[3],y[3];
  158.     y[0]=511;
  159.     y[1]=1274;
  160.     y[2]=511+1274;
  161.     x[0]=p[ch][1][1];
  162.     x[1]=p[ch][2][1];
  163.     x[2]=p[ch][3][1];
  164.  
  165.     //    TGraph *g=new TGraphErrors(3,x,y,ey,0);//gをxとyの散布図として宣言してつくります
  166.     TGraph *g=new TGraph(3,x,y);//gをxとyの散布図として宣言してつくります
  167.  
  168.     char gtitle[256];
  169.     sprintf(gtitle, "Calibration for ADC%d",ch);
  170. //グラフのタイトルを準備します    
  171.  
  172.     g->SetTitle(gtitle);//グラフのタイトルをgtitleに入っている値にしました
  173.     g->SetMarkerStyle(5);//グラフのデータ点の設定
  174.     g->SetMarkerSize(1);//データ点の大きさ
  175.     g->SetMarkerColor(1);//データ点の色
  176.     g->GetXaxis()->SetTitle("ADC Channel");//x軸の名前
  177.     g->GetYaxis()->SetTitle("Energy (keV)");
  178.     g->Draw("AP");//グラフをDrawします．APってオプションはなんだか忘れました
  179.  
  180.  
  181.  
  182.     char fgname[256];
  183.  
  184.     sprintf(fgname,"fg[%d]",ch);
  185.     TF1 *fg = new TF1(fgname,"pol1",0,4000);
  186.     g->Fit(fg,"","",0,4000); 
  187. //つくったグラフに対して1次元の1次関数fgをfitします 
  188.  
  189.     char func[256];
  190.  
  191.     Double_t pcal[9][2];
  192.     fg->GetParameters(pcal[ch]);    
  193.     sprintf(func,"y=%lfx+%lf",pcal[ch][1],pcal[ch][0]);
  194. //fitしたfgに対してパラメータを取得し，funcという文字列に式を書いておきます
  195.  
  196.     TLegend *lg = new TLegend(0.5255,0.40,0.875,0.15);
  197.     lg->AddEntry(g,"Data Point","p");
  198.     lg->AddEntry(fg,"Fitted line","l");
  199.     lg->AddEntry((TObject*)0, func, "");  
  200.     lg->SetFillColor(0); 
  201.     lg->Draw();
  202. //canvas2にもlgというlegendをつくりました．    
  203.  
  204.  
  205.  
  206.  
  207.     //次にヒストグラムの軸をkeVで書き直したものをつくります
  208.  
  209.     c3->cd(ch)->SetLogy();
  210.  
  211. //hstartとhendは最後に横軸をkeVで書き換えたあと、各チャンネルのプロットする範囲を決める変数です。
  212.  
  213.     double hstart, hend;
  214.     hstart=pcal[ch][1]*0+pcal[ch][0];
  215.     hend=pcal[ch][1]*4000.0+pcal[ch][0];
  216.     /* h[ch+10]->SetBins(4000,hstart,hend);
  217.        h[ch+10]->Draw();*/
  218. //x軸を書き直す作業はこれだけで終わりです
  219.  
  220. //あとはさっきと同じことをしてる
  221.  
  222.     //Gaussian: Resolution=FWHM/Mean
  223.     //FWHM=2*Sigma*sqrt(2*log2)
  224.     // 分解能出すときにcalibrationの後もう一度フィッティングしてからのほうがいいのかも
  225.     sprintf(mean1," Mean =%lf keV",pcal[ch][1]*p[ch][1][1]+pcal[ch][0]);
  226.     sprintf(resol1,"Resol.=%lf",2*sqrt(2*log(2.0))*(pcal[ch][1]*p[ch][1][2])/(pcal[ch][1]*p[ch][1][1]+pcal[ch][0]));
  227.     sprintf(mean2," Mean =%lf keV",pcal[ch][1]*p[ch][2][1]+pcal[ch][0]);
  228.     sprintf(resol2,"Resol.=%lf",2*sqrt(2*log(2.0))*(pcal[ch][1]*p[ch][2][2])/(pcal[ch][1]*p[ch][2][1]+pcal[ch][0]));
  229.     sprintf(mean3," Mean =%lf keV",pcal[ch][1]*p[ch][3][1]+pcal[ch][0]);
  230.     sprintf(resol3,"Resol.=%lf",2*sqrt(2*log(2.0))*(pcal[ch][1]*p[ch][3][2])/(pcal[ch][1]*p[ch][3][1]+pcal[ch][0]));
  231.  
  232.     TLegend *lc = new TLegend(0.6255,0.75,0.875,0.45);
  233.     //    lc->AddEntry(fit1,"511 keV","lp");
  234.     lc->AddEntry((TObject*)0, mean1,"");
  235.     lc->AddEntry((TObject*)0, resol1,"");
  236.     //    lc->AddEntry(fit2,"1274 keV","lp");
  237.     lc->AddEntry((TObject*)0, mean2,"");  
  238.     lc->AddEntry((TObject*)0, resol2,"");  
  239.     //    lc->AddEntry(fit3,"511+1274 keV","lp");
  240.     lc->AddEntry((TObject*)0, mean3,""); 
  241.     lc->AddEntry((TObject*)0, resol3,"");   
  242.     lc->SetFillColor(0);
  243.     lc->Draw();
  244.   }
  245.  
  246.  
  247. //こっからは新機能  
  248.   //1次式＋ガウシアンでフィットし直す
  249.  
  250.   char fglname[99];
  251.  
  252.   TF1 fgl[7][4];//function gauss + linear [ch No.][peak No.]
  253.  
  254.   TCanvas *c4 = new TCanvas("c4","Fitting with Gaussian and Linear function");
  255.   c4->Divide(2,3);
  256.  
  257.   Double_t paratest[7][5];
  258.  
  259.  
  260.   for (Int_t ch=1; ch<7; ch++){
  261.  
  262.     c4->cd(ch)->SetLogy();
  263.     h[ch]->Draw();
  264.  
  265.     sprintf(fglname,"fgl1%d",ch);
  266.     fgl[ch][1] = 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
  267.     fgl[ch][1]->SetParameters(p[ch][1][0],p[ch][1][2],p[ch][1][1],0,0);//パラメータの初期値の設定
  268.     h[ch]->Fit(fglname,"","",p1start[ch],p1end[ch]);
  269.  
  270.     fgl[ch][1]->GetParameters(paratest[ch]);// parameterをgetしてこれを使えば直線の式が出せるはずだがなんかできない
  271.  
  272.  
  273.     //ftest= new TF1("ftest","[0]*x+[1]",p1start[ch],p1end[ch]);
  274.     TF1 *ftest = new TF1("ftest","[0]*x+[1]");
  275.     ftest->SetParameter(0,paratest[ch][3]);
  276.     ftest->SetParameter(1,paratest[ch][4]);
  277.  
  278.     ftest->SetLineColor(7);
  279.  
  280.  
  281.     fgl[ch][1]->Draw("sames");
  282.     ftest->Draw("sames");	
  283.  
  284.  
  285.     /*   //2番目3番目のピークを重ね書きするのはまだできてない
  286. 	 sprintf(fglname,"fgl2%d",ch);
  287. 	 fgl[ch][2] = 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
  288. 	 fgl[ch][2]->SetParameters(p[ch][2][0],p[ch][2][2],p[ch][2][1],0,0);
  289. 	 h[ch]->Fit(fglname,"","",p2start[ch],p2end[ch]);
  290. 	 fgl[ch][2]->Draw("sames");
  291.  
  292. 	 sprintf(fglname,"fgl3%d",ch);
  293. 	 fgl[ch][3] = 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
  294. 	 fgl[ch][3]->SetParameters(p[ch][3][0],p[ch][3][2],p[ch][3][1],0,0);
  295. 	 h[ch]->Fit(fglname,"","",p3start[ch],p3end[ch]);
  296. 	 fgl[ch][3]->Draw("sames");
  297.     */
  298.   }
  299.  
  300.  
  301.  
  302.   //キャンバスをPDFに出力します
  303.   //Print
  304.   char pdfstart[99];
  305.   sprintf(pdfstart,"%s.pdf(",filename);
  306.   char pdfmiddle[99];
  307.   sprintf(pdfmiddle,"%s.pdf",filename);
  308.   char pdfend[99];
  309.   sprintf(pdfend,"%s.pdf)",filename);
  310.  
  311.   c1->Print(pdfstart);
  312.   c2->Print(pdfmiddle);
  313.   c3->Print(pdfmiddle);
  314.   c4->Print(pdfend);
  315.  
  316.  
  317. //この下では別のプログラムで使うためのヘッダファイルをつくっています．あんまりきにしなくてもいいはず
  318.  
  319.   //フィッティングパラメーターをファイルに書きこむ
  320.   char pfile[99];
  321.   sprintf(pfile,"parameter%s.h",filename);
  322.   FILE *fp;
  323.   fp=fopen(pfile,"w");
  324.   fprintf(fp,"Double_t para1[7], para2[7],para3[7],para4[7];\nDouble_t para12,para13,para14,para15,para16, para23,para24,para25,para26,para34,para35,para36,para45,para46,para56;\n\n");
  325.   fprintf(fp,"//Energy = Naiadc[ch]*para1[ch]+para2[ch]\n");
  326.  
  327.   //Writing ADC Parameters
  328.   char para1[99];
  329.   char para2[99];
  330.   for(Int_t ch=1; ch<7; ch++){
  331.     sprintf(para1,"para1[%d]=%lf;\n",ch,pcal[ch-1][3]);
  332.     fprintf(fp,para1);
  333.     sprintf(para2,"para2[%d]=%lf;\n",ch,pcal[ch-1][2]);
  334.     fprintf(fp,para2);
  335.   }
  336.  
  337.   //Writing TDC Parameters
  338.   fprintf(fp,"para3[1]=0.5707;\npara4[1]=42.063;\npara3[2]=0.5679;\npara4[2]=47.576;\npara3[3]=0.5361;\npara4[3]=67.592;\npara3[4]=0.6024;\npara4[4]=12.809;\npara3[5]=0.5663;\npara4[5]=50.119;\npara3[6]=0.56559;\npara4[6]=45.407;\n\n");
  339.   fprintf(fp,"para12=0;\npara13=0;\npara14=0;\npara15=0;\npara16=0;\npara23=0;\npara24=0;\npara25=0;\npara26=0;\npara34=0;\npara35=0;\npara36=0;\npara45=0;\npara46=0;\npara56=0;\n");  
  340.  
  341.   fclose(fp);
  342.  
  343. }
  344.  
  345.  
「macro1221」をウィキ内検索
最終更新：2011年12月21日 23:40