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

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

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

今回
  • x軸をkeV, yをchにした
  • Gaussian+linearもできた

  • 511のsigmaが期待よりかなりでかいのはなぜ

To Do 

  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_lin(){  //Last modified 2011/11/22
  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(filname,"22na20111104.001.root");
  24.   sprintf(filename,"22na-10.03.001-002.root");
  25.   //sprintf(filename,"22na-11.10.20-21.root");
  26.   //sprintf(filename,"22na-20111028.002.root");
  27.   //  sprintf(filename,"22na20111104.001.root");
  28.   //  sprintf(filename,"22na-20111111.001.root");
  29.   //  sprintf(filename,"22na-20111118.001.root");
  30.   //  sprintf(filename,"22na-20111201.001.root");  
  31. //filenameという変数を宣言して，そこにファイル名を代入してます．この変数はファイルを開くのと，あとでPDFファイルをつくるときに使います．sprintfはprintfと似た関数で，変数に文字列をいれることができます．ROOTのマクロを書くときにはけっこう便利です．あと、文字列をいれるときは，charというクラスでその文字数以上の配列として宣言するみたいです。
  32.  
  33.  
  34.   TFile *f1=new TFile(filename);//opening root file
  35. //ファイルを開いています。
  36.  
  37.   Double_t p[9][9][9];
  38. //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です。
  39.  
  40.  
  41.  
  42.   TH1F *h[99];
  43.   h[1]=h101;
  44.   h[2]=h102;
  45.   h[3]=h103;
  46.   h[4]=h104;
  47.   h[5]=h105;
  48.   h[6]=h106;
  49. //TH1Fは1次のFloatのHistgramのクラスです。ROOTファイルに入っているh101～h106までのヒストグラムをプログラム中で扱いやすいようにh[i]という配列を宣言して、そこに代入しています。
  50.  
  51.  
  52.  
  53.   // end of front matter
  54.  
  55.  
  56.   //gStyle->SetOptFit(0001);
  57.   gStyle->SetOptFit(0000);
  58. //ここはフィッティングのパラメータをどのようにキャンバスに出力するか決めています  
  59.  
  60.  
  61.  
  62.   TCanvas *c1 = new TCanvas("c1","Raw Data and Fitting");
  63.   c1->Divide(2,3);
  64.  
  65.   TCanvas *c2 = new TCanvas("c2","Enegy Calibration");
  66.   c2->Divide(2,3);
  67.  
  68.   TCanvas *c3 = new TCanvas("c3","Calibrated Hist");
  69.   c3->Divide(2,3);
  70.  
  71.   TCanvas *c4 = new TCanvas("c4","Fitting with Gaussian and Linear function");
  72.   c4->Divide(2,3);
  73.  
  74.   TCanvas *c5 = new TCanvas("c5","Energy Calibration with G+L");
  75.   c5->Divide(2,3);
  76.  
  77.   TCanvas *c6 = new TCanvas("c6","Calibrated Hist with G+L");
  78.   c6->Divide(2,3);
  79.  
  80.   TCanvas *c7 = new TCanvas("c7","Energy dependence of Sigma");
  81.   c7->Divide(2,3);
  82.  
  83.   TCanvas *c8 = new TCanvas("c8","E dependence of Sigma with G+L");
  84.   c8->Divide(2,3);
  85.  
  86.  
  87.  
  88. //TCanvasというクラスはヒストグラムとかを出力するためのクラスです．キャンバスをDivideするとキャンバスを複数に分けることができます
  89.  
  90.  
  91.  
  92.  
  93.   for(Int_t ch=1;ch<=6;ch++){
  94.  
  95.     char hclone[256]; 
  96.     sprintf(hclone,"hclone%d",ch);
  97.     TH1F *hoge=h[ch]->Clone(hclone);
  98.     h[ch+10]=hoge;
  99.  
  100.     sprintf(hclone,"hclone2%d",ch);
  101.     TH1F *hoge2=h[ch]->Clone(hclone);
  102.     h[ch+20]=hoge2;
  103.     //h101-h106を違う目的で使いたいので複製しています
  104.   }
  105.  
  106.   Double_t sigma1[3],sigma2[3];
  107.   Double_t sigmae1[3],sigmae2[3];
  108.  
  109.  
  110.   for(Int_t ch=1;ch<=6;ch++){
  111.  
  112.     c1->cd(ch)->SetLogy();    
  113.     h[ch]->Draw();
  114.     //まず，キャンバスc1の6つに分割したうちの1番目に移動し，ch=1のときはh[1]=h101をDrawします．
  115.  
  116.     //Preparation of Fitting function
  117.     //ここからはピークをフィッティングするための関数を準備しています．
  118.  
  119.     TF1 *fit1;    
  120.     char fname[256];
  121.     char mean1[256],mean2[256],mean3[256];
  122.     char resol1[256],resol2[256],resol3[256];
  123.  
  124.     sprintf(fname,"f[%d]",ch);
  125.     //fnameは関数の名前で，ひとつめのピークにはチャネルごとにf[1]-f[6]をフィッティングします
  126.     fit1 = new TF1(fname,"gaus",p1start[ch],p1end[ch]);
  127.     //この行ではフィッティングのための関数fit1をつくっています．名前はfnameに入っている値，種類はガウシアン，定義域はp1startとp1endに入っている値です
  128.     h[ch]->Fit(fit1,"","",p1start[ch],p1end[ch]); 
  129.     //h[ch]にfit1をフィッティングします
  130.     fit1->GetParameters(p[ch][1]);
  131.     sigmae1[0]=fit1->GetParError(2);
  132.     //フィットしたfit1について，そのフィッティングパラメータを得ます．これによってp[ch][1][1]にはMeanが，p[ch][1][2]にはSigmaが入ります．
  133.     sprintf(mean1," Mean =%lf",p[ch][1][1]);
  134.     sprintf(resol1,"Resol.=%lf",2*sqrt(2*log(2.0))*p[ch][1][2]/p[ch][1][1]);
  135.     //凡例に表示するための文字列を準備しています
  136.  
  137.  
  138.     //同じことを2つ目，3つ目のピークについても行います．
  139.     TF1 *fit2;
  140.     sprintf(fname,"f[1%d]",ch);
  141.     fit2 = new TF1(fname,"gaus",p2start[ch],p2end[ch]);
  142.     h[ch]->Fit(fit2,"+","",p2start[ch],p2end[ch]); 
  143.     fit2->GetParameters(p[ch][2]);
  144.     sigmae1[1]=fit2->GetParError(2);    
  145.     sprintf(mean2," Mean =%lf",p[ch][2][1]);
  146.     sprintf(resol2,"Resol.=%lf",2*sqrt(2*log(2.0))*p[ch][2][2]/p[ch][2][1]);
  147.  
  148.  
  149.     TF1 *fit3;
  150.     sprintf(fname,"f[10%d]",ch);
  151.     fit3 = new TF1(fname,"gaus",p3start[ch],p3end[ch]);
  152.     h[ch]->Fit(fname,"+","",p3start[ch],p3end[ch]); 
  153.     fit3->GetParameters(p[ch][3]);    
  154.     sigmae1[2]=fit3->GetParError(2);
  155.     sprintf(mean3," Mean =%lf",p[ch][3][1]);
  156.     sprintf(resol3,"Resol.=%lf",2*sqrt(2*log(2.0))*p[ch][3][2]/p[ch][3][1]);
  157.  
  158.     h[ch]->Draw();//もう１回Drawします．なぜかこうしないとうまくいかなかったきがする
  159.  
  160.     // fit1->Draw("sames");//fit1をDrawします．"sames"と書くと上書きされます．
  161.  
  162.     //    fit2->SetLineColor(3);
  163.     //fit2->Draw("sames");//fit2の線の色を設定したのち，上書きします
  164.  
  165.     //fit3->SetLineColor(7);
  166.     //  fit3->Draw("sames");
  167.  
  168.  
  169.     //ここからは凡例(Legend)の設定です    
  170.     TLegend *l = new TLegend(0.6255,0.75,0.875,0.45);
  171.     //TLegendクラスの変数 l を宣言しました
  172.     l->AddEntry(fit1,"511 keV","lp");
  173.     //lにfit1を追加し，その表示名を511 keVとしました
  174.     l->AddEntry((TObject*)0, mean1,"");
  175.     //文字列mean1をlに追加しました
  176.     l->AddEntry((TObject*)0, resol1,"");
  177.     l->AddEntry(fit2,"1274 keV","lp");
  178.     l->AddEntry((TObject*)0, mean2,"");  
  179.     l->AddEntry((TObject*)0, resol2,"");  
  180.     l->AddEntry(fit3,"511+1274 keV","lp");
  181.     l->AddEntry((TObject*)0, mean3,""); 
  182.     l->AddEntry((TObject*)0, resol3,"");   
  183.     l->SetFillColor(0);// l の背景を白にします
  184.     l->Draw();//lをキャンバスに描きます
  185.  
  186.     // Graph of calibration   
  187.     //エネルギーキャリブレーションの式を求めるグラフをつくります
  188.  
  189.     c2->cd(ch);
  190.  
  191.     //xとyの散布図をつくるために，値を準備しておきます  
  192.     Double_t x[3],y[3];
  193.     x[0]=510.9989;
  194.     x[1]=1274.58;
  195.     x[2]=510.9989+1274.58;
  196.     y[0]=p[ch][1][1];
  197.     y[1]=p[ch][2][1];
  198.     y[2]=p[ch][3][1];
  199.  
  200.  
  201.     TGraph *g=new TGraph(3,x,y);//gをxとyの散布図として宣言してつくります
  202.  
  203.     char gtitle[256];
  204.     sprintf(gtitle, "Calibration for ADC%d",ch);
  205.     //グラフのタイトルを準備します    
  206.  
  207.     g->SetTitle(gtitle);//グラフのタイトルをgtitleに入っている値にしました
  208.     g->SetMarkerStyle(5);//グラフのデータ点の設定
  209.     g->SetMarkerSize(3);//データ点の大きさ
  210.     g->SetMarkerColor(1);//データ点の色
  211.     g->GetYaxis()->SetTitle("ADC Channel");//x軸の名前
  212.     g->GetXaxis()->SetTitle("Energy (keV)");
  213.     g->Draw("AP");//グラフをDrawします．APってオプションはなんだか忘れました
  214.  
  215.  
  216.  
  217.     char fgname[256];
  218.  
  219.     sprintf(fgname,"fg[%d]",ch);
  220.     //    TF1 *fg = new TF1(fgname,"pol1",0,4000);
  221.     TF1 *fg = new TF1(fgname,"[1]*x+[0]",0,4000);
  222.     g->Fit(fg,"","",0,4000); 
  223.     //つくったグラフに対して1次元の1次関数fgをfitします 
  224.  
  225.     char func[256];
  226.  
  227.     Double_t pcal[9][5];
  228.     fg->GetParameters(pcal[ch]);    
  229.     sprintf(func,"y=%lfx+%lf",pcal[ch][1],pcal[ch][0]);
  230. //fitしたfgに対してパラメータを取得し，funcという文字列に式を書いておきます
  231.  
  232.     TLegend *lg = new TLegend(0.5255,0.40,0.875,0.15);
  233.     lg->AddEntry(g,"Data Point","p");
  234.     lg->AddEntry(fg,"Fitted line","l");
  235.     lg->AddEntry((TObject*)0, func, "");  
  236.     lg->SetFillColor(0); 
  237.     lg->Draw();
  238. //canvas2にもlgというlegendをつくりました．    
  239.  
  240.  
  241.     //次にヒストグラムの軸をkeVで書き直したものをつくります
  242.  
  243.     c3->cd(ch)->SetLogy();
  244.  
  245.     //hstartとhendは最後に横軸をkeVで書き換えたあと、各チャンネルのプロットする範囲を決める変数です。
  246.  
  247.     double hstart, hend;
  248.     hstart=pcal[ch][1]*0+pcal[ch][0];
  249.     hend=pcal[ch][1]*4000.0+pcal[ch][0];
  250.     h[ch+10]->SetBins(4000,hstart,hend);
  251.     h[ch+10]->Draw();
  252.     //x軸を書き直す作業はこれだけで終わりです
  253.  
  254.     //あとはさっきと同じことをしてる
  255.  
  256.     //Gaussian: Resolution=FWHM/Mean
  257.     //FWHM=2*Sigma*sqrt(2*log2)
  258.     // 分解能出すときにcalibrationの後もう一度フィッティングしてからのほうがいいのかも
  259.     sprintf(mean1," Mean =%lf keV",1/pcal[ch][1]*p[ch][1][1]-pcal[ch][0]/pcal[ch][1]);
  260.     sprintf(resol1,"Sigma/mean=%lf",(1/pcal[ch][1]*p[ch][1][2])/(1/pcal[ch][1]*p[ch][1][1]-pcal[ch][0]/pcal[ch][1]));
  261.     sprintf(mean2," Mean =%lf keV",1/pcal[ch][1]*p[ch][2][1]-pcal[ch][0]/pcal[ch][1]);
  262.     sprintf(resol2,"Sigma/mean=%lf",(1/pcal[ch][1]*p[ch][2][2])/(1/pcal[ch][1]*p[ch][2][1]-pcal[ch][0]/pcal[ch][1]));
  263.     sprintf(mean3," Mean =%lf keV",1/pcal[ch][1]*p[ch][3][1]-pcal[ch][0]/pcal[ch][1]);
  264.     sprintf(resol3,"Sigma/mean=%lf",(1/pcal[ch][1]*p[ch][3][2])/(1/pcal[ch][1]*p[ch][3][1]-pcal[ch][0]/pcal[ch][1]));
  265.  
  266.  
  267.     TLegend *lc = new TLegend(0.6255,0.75,0.875,0.45);
  268.     lc->AddEntry(fit1,"511 keV","lp");
  269.     lc->AddEntry((TObject*)0, mean1,"");
  270.     lc->AddEntry((TObject*)0, resol1,"");
  271.     lc->AddEntry(fit2,"1274 keV","lp");
  272.     lc->AddEntry((TObject*)0, mean2,"");  
  273.     lc->AddEntry((TObject*)0, resol2,"");  
  274.     lc->AddEntry(fit3,"511+1274 keV","lp");
  275.     lc->AddEntry((TObject*)0, mean3,""); 
  276.     lc->AddEntry((TObject*)0, resol3,"");   
  277.     lc->SetFillColor(0);
  278.     lc->Draw();
  279.  
  280.  
  281.  
  282.     //こっからは新機能  
  283.     //1次式＋ガウシアンでフィットし直す
  284.  
  285.     char fglname[99];
  286.     Double_t paratest[5];
  287.     Double_t ex[3];
  288.  
  289.     c4->cd(ch)->SetLogy();
  290.     //c4->cd(ch);
  291.  
  292.     h[ch+20]->Draw();
  293.     TLegend *l3 = new TLegend(0.6255,0.75,0.875,0.45);
  294.  
  295.  
  296.  
  297.     TF1 *fit4;
  298.     sprintf(fglname,"fgl1%d",ch);
  299.     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
  300.     fit4->SetParameters(p[ch][1][0],p[ch][1][2],p[ch][1][1],0,0);//パラメータの初期値の設定
  301.     h[ch+20]->Fit(fglname,"","",p1start[ch],p1end[ch]);
  302.     fit4 = h[ch+20]->GetFunction(fglname);
  303.     fit4->GetParameters(paratest);// parameterをgetしてこれを使えば直線の式が
  304.     x[0]=paratest[2];
  305.     ex[0]=fit4->GetParError(2);
  306.     sigma2[0]=paratest[1];
  307.     sigmae2[0]=fit4->GetParError(1);
  308.  
  309.  
  310.     TF1 *flin1= new TF1("flin1","[0]*x+[1]",p1start[ch],p1end[ch]);
  311.     flin1->FixParameter(0,paratest[3]);
  312.     flin1->FixParameter(1,paratest[4]);//"fix" is proper, not "set"
  313.     h[ch+20]->Fit(flin1,"Q0+","",p1start[ch],p1end[ch]);
  314.     flin1=h[ch+20]->GetFunction("flin1");
  315.  
  316.     char caption[32]; 
  317.     l3->AddEntry(fit4,"511 keV","lp");
  318.     sprintf(caption,"mean=%lf",paratest[2]);
  319.     l3->AddEntry((TObject*)0, caption,"");
  320.     sprintf(caption,"sigma=%lf",paratest[1]);
  321.     l3->AddEntry((TObject*)0, caption,"");
  322.  
  323.  
  324.     TF1 *fit5;
  325.     sprintf(fglname,"fgl2%d",ch);
  326.     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
  327.     fit5->SetParameters(p[ch][2][0],p[ch][2][2],p[ch][2][1],0,0);
  328.     h[ch+20]->Fit(fglname,"0+","",p2start[ch],p2end[ch]);
  329.     fit5 = h[ch+20]->GetFunction(fglname);
  330.     fit5->GetParameters(paratest);
  331.     x[1]=paratest[2];
  332.     ex[1]=fit5->GetParError(2);
  333.     sigma2[1]=paratest[1];
  334.     sigmae2[1]=fit5->GetParError(1);
  335.  
  336.     TF1 *flin2= new TF1("flin2","[0]*x+[1]",p2start[ch],p2end[ch]);
  337.     flin2->FixParameter(0,paratest[3]);
  338.     flin2->FixParameter(1,paratest[4]);//"fix" is proper, not "set"
  339.     h[ch+20]->Fit(flin2,"Q0+","",p2start[ch],p2end[ch]);
  340.     flin2=h[ch+20]->GetFunction("flin2");
  341.  
  342.     char caption[32]; 
  343.     l3->AddEntry(fit5,"1274.5 keV","lp");
  344.     sprintf(caption,"mean=%lf",paratest[2]);
  345.     l3->AddEntry((TObject*)0, caption,"");
  346.     sprintf(caption,"sigma=%lf",paratest[1]);
  347.     l3->AddEntry((TObject*)0, caption,"");
  348.  
  349.  
  350.  
  351.     TF1 *fit6;	 	 
  352.     sprintf(fglname,"fgl3%d",ch);
  353.     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
  354.     fit6->SetParameters(p[ch][3][0],p[ch][3][2],p[ch][3][1],0,0);
  355.     h[ch+20]->Fit(fglname,"0+","",p3start[ch],p3end[ch]);
  356.     fit6 = h[ch+20]->GetFunction(fglname);
  357.     fit6->GetParameters(paratest);
  358.     x[2]=paratest[2];
  359.     ex[2]=fit6->GetParError(2);
  360.     sigma2[2]=paratest[1];
  361.     sigmae2[2]=fit6->GetParError(1);
  362.     TF1 *flin3= new TF1("flin3","[0]*x+[1]",p3start[ch],p3end[ch]);
  363.     flin3->FixParameter(0,paratest[3]);
  364.     flin3->FixParameter(1,paratest[4]);//"fix" is proper, not "set"
  365.     h[ch+20]->Fit(flin3,"Q0+","",p3start[ch],p3end[ch]);
  366.     flin3=h[ch+20]->GetFunction("flin3");
  367.  
  368.     char caption[32]; 
  369.     l3->AddEntry(fit6,"511+1274 keV","lp");
  370.     sprintf(caption,"mean=%lf",paratest[2]);
  371.     l3->AddEntry((TObject*)0, caption,"");
  372.     sprintf(caption,"sigma=%lf",paratest[1]);
  373.     l3->AddEntry((TObject*)0, caption,"");
  374.  
  375.  
  376.  
  377.     h[ch+20]->Draw("same");
  378.  
  379.     fit4->Draw("same");  
  380.     flin1->Draw("same");		 
  381.  
  382.     fit5->SetLineColor(3);
  383.     fit5->Draw("same");  
  384.     flin2->SetLineColor(3); 
  385.     flin2->Draw("same");		 
  386.  
  387.     fit6->SetLineColor(7);
  388.     fit6->Draw("same"); 
  389.     flin3->SetLineColor(7);
  390.     flin3->Draw("same");	
  391.  
  392.     l3->SetFillColor(0);// l の背景を白にします
  393.     l3->Draw();//lをキャンバスに描きます
  394.  
  395.     /* Energy Calibration*/
  396.     c5->cd(ch);
  397.  
  398.  
  399.     //  TGraph *g=new TGraphErrors(3,x,y,ex,0);//なぜかxだけ誤差ありだとフィッティングできない
  400.     //    TGraph *g=new TGraph(3,x,y);
  401.     TGraph *g=new TGraphErrors(3,y,x,0,ex);//なぜかxだけ誤差ありだとフィッティングできない
  402.  
  403.     char gtitle[256];
  404.     sprintf(gtitle, "Calibration for ADC%d",ch);
  405.     //グラフのタイトルを準備します    
  406.  
  407.     g->SetTitle(gtitle);//グラフのタイトルをgtitleに入っている値にしました
  408.     g->SetMarkerStyle(5);//グラフのデータ点の設定
  409.     g->SetMarkerSize(2);//データ点の大きさ
  410.     g->SetMarkerColor(1);//データ点の色
  411.     g->GetYaxis()->SetTitle("ADC Channel");//x軸の名前
  412.     g->GetXaxis()->SetTitle("Energy (keV)");
  413.     g->Draw("AP");//グラフをDrawします．APってオプションはなんだか忘れました
  414.  
  415.     char fgname[256];
  416.     sprintf(fgname,"fg2",ch);
  417.     //    TF1 *fg2 = new TF1(fgname,"[2]*x**2+[1]*x+[0]",0,4000);
  418.     TF1 *fg2 = new TF1(fgname,"[1]*x+[0]",0,4000);
  419.     g->Fit(fg2,"","",0,4000);
  420.     fg=g->GetFunction("fg2");
  421.  
  422.     char func[256];
  423.  
  424.     Double_t pcal2[9][5];
  425.     fg->GetParameters(pcal2[ch]);    
  426.     //    sprintf(func,"y=%lfx^2+%lfx+%lf",pcal2[ch][2],pcal2[ch][1],pcal2[ch][0]);
  427.     sprintf(func,"y=%lfx+%lf",pcal2[ch][1],pcal2[ch][0]);
  428.  
  429.     TLegend *lg = new TLegend(0.5255,0.40,0.875,0.15);
  430.     lg->AddEntry(g,"Data Point","p");
  431.     lg->AddEntry(fg,"Fitted line","l");
  432.     lg->AddEntry((TObject*)0, func, "");  
  433.     lg->SetFillColor(0); 
  434.     lg->Draw();
  435.     fg->Draw("same");
  436.  
  437.     c6->cd(ch);
  438.  
  439.     sprintf(mean1," Mean =%lf keV",1/pcal2[ch][1]*p[ch][1][1]-pcal2[ch][0]/pcal2[ch][1]);
  440.     sprintf(resol1,"Sigma/mean=%lf",(1/pcal2[ch][1]*p[ch][1][2])/(1/pcal2[ch][1]*p[ch][1][1]-pcal2[ch][0]/pcal2[ch][1]));
  441.     sprintf(mean2," Mean =%lf keV",1/pcal2[ch][1]*p[ch][2][1]-pcal2[ch][0]/pcal2[ch][1]);
  442.     sprintf(resol2,"Sigma/mean=%lf",(1/pcal2[ch][1]*p[ch][2][2])/(1/pcal2[ch][1]*p[ch][2][1]-pcal2[ch][0]/pcal2[ch][1]));
  443.     sprintf(mean3," Mean =%lf keV",1/pcal2[ch][1]*p[ch][3][1]-pcal2[ch][0]/pcal2[ch][1]);
  444.     sprintf(resol3,"Sigma/mean=%lf",(1/pcal2[ch][1]*p[ch][3][2])/(1/pcal2[ch][1]*p[ch][3][1]-pcal2[ch][0]/pcal2[ch][1]));
  445.  
  446.  
  447.     TLegend *lc = new TLegend(0.6255,0.75,0.875,0.45);
  448.     lc->AddEntry(fit1,"511 keV","lp");
  449.     lc->AddEntry((TObject*)0, mean1,"");
  450.     lc->AddEntry((TObject*)0, resol1,"");
  451.     lc->AddEntry(fit2,"1274 keV","lp");
  452.     lc->AddEntry((TObject*)0, mean2,"");  
  453.     lc->AddEntry((TObject*)0, resol2,"");  
  454.     lc->AddEntry(fit3,"511+1274 keV","lp");
  455.     lc->AddEntry((TObject*)0, mean3,""); 
  456.     lc->AddEntry((TObject*)0, resol3,"");   
  457.     lc->SetFillColor(0);
  458.     lc->Draw();
  459.  
  460.  
  461.     c7->cd(ch);
  462.     //y is already set for energy
  463.  
  464.  
  465.  
  466.     sigma1[0]= p[ch][1][2] / pcal[ch][1];  
  467.     sigma1[1]= p[ch][2][2] / pcal[ch][1];
  468.     sigma1[2]= p[ch][3][2] / pcal[ch][1];
  469.     y[0]= p[ch][1][1] / pcal[ch][1] -pcal[ch][0] / pcal[ch][1];
  470.     y[1]= p[ch][2][1] / pcal[ch][1] -pcal[ch][0] / pcal[ch][1];
  471.     y[2]= p[ch][3][1] / pcal[ch][1] -pcal[ch][0] / pcal[ch][1];
  472.     sigmae1[0]= sigmae1[0] / pcal[ch][1];
  473.     sigmae1[1]= sigmae1[1] / pcal[ch][1];
  474.     sigmae1[2]= sigmae1[2] / pcal[ch][1];
  475.  
  476.     /*
  477.     sigma1[0]=4.611 / pcal[1][1];//using pedestal of ch1  
  478.     sigma1[0]=12.0 / pcal[2][1];//using pedestal of ch2
  479.     //  sigma1[0]=3.37 / pcal[3][1];//using pedestal of ch3
  480.       // sigma1[0]=3.88945 / pcal[4][1];//using pedestal of ch4
  481.     //   sigma1[0]=3.24 / pcal[5][1];//using pedestal of ch5
  482.     //    sigma1[0]=6 / pcal[6][1];//using pedestal of ch6
  483.     */
  484.  
  485.     sigma1[0]=39;//using Cs of ch1  
  486.     //    sigma1[0]=44;//using Cs of ch2
  487.     //  sigma1[0]=37;//using Cs of ch3
  488.       // sigma1[0]=46;//using Cs of ch4
  489.        sigma1[0]=40;//using Cs of ch5
  490.     //    sigma1[0]=57;//using Cs of ch6
  491.  
  492.     y[0]=662;
  493.     sigmae1[0]= sigmae1[1] ;
  494.  
  495.  
  496.  
  497.     TGraph *g=new TGraphErrors(3,y,sigma1,0,sigmae1);
  498.  
  499.  
  500.     sprintf(gtitle, "Sigma(E) of ADC%d (Gaussian Only)",ch);
  501.     //グラフのタイトルを準備します    
  502.  
  503.     g->SetTitle(gtitle);//グラフのタイトルをgtitleに入っている値にしました
  504.     g->SetMarkerStyle(1);//グラフのデータ点の設定
  505.     g->SetMarkerSize(1);//データ点の大きさ
  506.     g->SetMarkerColor(1);//データ点の色
  507.     g->GetXaxis()->SetTitle("Energy (keV)");//x軸の名前
  508.     g->GetYaxis()->SetTitle("Sigma (keV)");
  509.     g->Draw("AP");//グラフをDrawします．APってオプションはなんだか忘れました
  510.  
  511.     char fgname[256];
  512.     sprintf(fgname,"fg2",ch);
  513.     TF1 *fg2 = new TF1(fgname,"[1]*x**(0.5)+[0]",0,4000);
  514.     g->Fit(fg2,"","",0,4000);
  515.     fg=g->GetFunction("fg2");
  516.  
  517.     char func[256];
  518.  
  519.     Double_t psigma1[9][5];
  520.     fg->GetParameters(psigma1[ch]);    
  521.     sprintf(func,"y=+%lf*sqrt(x)+%lf",psigma1[ch][1],psigma1[ch][0]);
  522.  
  523.     TLegend *lg = new TLegend(0.5255,0.40,0.875,0.15);
  524.     lg->AddEntry(g,"Data Point","p");
  525.     lg->AddEntry(fg,"Fitted line","l");
  526.     lg->AddEntry((TObject*)0, func, "");  
  527.     lg->SetFillColor(0); 
  528.     lg->Draw();
  529.     fg->Draw("same");
  530.  
  531.     /*
  532.   
  533.           Canvas No. 8
  534.  
  535.      */
  536.     c8->cd(ch);
  537.  
  538.     sigma2[0] = sigma2[0] / pcal2[ch][1];  
  539.     sigma2[1] = sigma2[1] / pcal2[ch][1];
  540.     sigma2[2] = sigma2[2] / pcal2[ch][1];
  541.     x[0] = x[0] / pcal2[ch][1] - pcal2[ch][0] / pcal2[ch][1];
  542.     x[1] = x[1] / pcal2[ch][1] - pcal2[ch][0] / pcal2[ch][1];
  543.     x[2] = x[2] / pcal2[ch][1] - pcal2[ch][0] / pcal2[ch][1];
  544.     sigmae2[0] = sigmae2[0] / pcal2[ch][1];
  545.     sigmae2[1] = sigmae2[1] / pcal2[ch][1];
  546.     sigmae2[2] = sigmae2[2] / pcal2[ch][1];
  547.  
  548.     TGraph *g=new TGraphErrors(3,x,sigma2,0,sigmae2);
  549.  
  550.  
  551.     sprintf(gtitle, "Sigma(E) of ADC%d (G + Linear)",ch);
  552.     //グラフのタイトルを準備します    
  553.  
  554.     g->SetTitle(gtitle);//グラフのタイトルをgtitleに入っている値にしました
  555.     g->SetMarkerStyle(1);//グラフのデータ点の設定
  556.     g->SetMarkerSize(1);//データ点の大きさ
  557.     g->SetMarkerColor(1);//データ点の色
  558.     g->GetXaxis()->SetTitle("Energy (keV)");//x軸の名前
  559.     g->GetYaxis()->SetTitle("Sigma (keV)");
  560.     g->Draw("AP");//グラフをDrawします．APってオプションはなんだか忘れました
  561.  
  562.     char fgname[256];
  563.     sprintf(fgname,"fg2",ch);
  564.     TF1 *fg2 = new TF1(fgname,"[1]*x**(0.5)+[0]",0,4000);
  565.     g->Fit(fg2,"","",0,4000);
  566.     fg=g->GetFunction("fg2");
  567.  
  568.     char func[256];
  569.  
  570.     Double_t psigma1[9][5];
  571.     fg->GetParameters(psigma1[ch]);    
  572.     sprintf(func,"y=+%lf*sqrt(x)+%lf",psigma1[ch][1],psigma1[ch][0]);
  573.  
  574.     TLegend *lg = new TLegend(0.5255,0.40,0.875,0.15);
  575.     lg->AddEntry(g,"Data Point","p");
  576.     lg->AddEntry(fg,"Fitted line","l");
  577.     lg->AddEntry((TObject*)0, func, "");  
  578.     lg->SetFillColor(0); 
  579.     lg->Draw();
  580.     fg->Draw("same");
  581.  
  582.   }
  583.  
  584.  
  585.  
  586.   //キャンバスをPDFに出力します
  587.   //Print
  588.   char pdfstart[99];
  589.   sprintf(pdfstart,"%s.pdf(",filename);
  590.   char pdfmiddle[99];
  591.   sprintf(pdfmiddle,"%s.pdf",filename);
  592.   char pdfend[99];
  593.   sprintf(pdfend,"%s.pdf)",filename);
  594.  
  595.   c1->Print(pdfstart);
  596.   c2->Print(pdfmiddle);
  597.   c3->Print(pdfmiddle);
  598.   c4->Print(pdfend);
  599.  
  600.  
  601. //この下では別のプログラムで使うためのヘッダファイルをつくっています．あんまりきにしなくてもいいはず
  602.  
  603.   //フィッティングパラメーターをファイルに書きこむ
  604.   char pfile[99];
  605.   sprintf(pfile,"parameter%s.h",filename);
  606.   FILE *fp;
  607.   fp=fopen(pfile,"w");
  608.   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");
  609.   fprintf(fp,"//Energy = Naiadc[ch]*para1[ch]+para2[ch]\n");
  610.  
  611.   //Writing ADC Parameters
  612.   char para1[99];
  613.   char para2[99];
  614.   for(Int_t ch=1; ch<7; ch++){
  615.     sprintf(para1,"para1[%d]=%lf;\n",ch,pcal[ch-1][3]);
  616.     fprintf(fp,para1);
  617.     sprintf(para2,"para2[%d]=%lf;\n",ch,pcal[ch-1][2]);
  618.     fprintf(fp,para2);
  619.   }
  620.  
  621.   //Writing TDC Parameters
  622.   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");
  623.   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");  
  624.  
  625.   fclose(fp);
  626.  
  627. }
  628.  
  629.  
  630.  
「macro1228」をウィキ内検索
最終更新：2011年12月28日 01:02