卒研めも/Calibration/Cs-137/macro0109 - おかめも

わりと完成版。まず661kevをGaussianだけでフィットして、そのデータを利用してガウシアン＋1次式でフィット。出力されるヒストグラムのレジェンドにはmeanとsigmaをchとkeVの単位で書く。あとsigmaのerrorもkeV単位で。

#include<fstream>
#include<TH1.h>
#include<TH2.h>
#include<TFile.h>
#include<TTree.h>
#include<TCanvas.h>
#include<parameter_et.h>
#include<parameter_cs.h>
#include<math.h>
//#include<iostream>
 
#include<iostream> using namespace std;
 
//ここまでで使うライブラリをincludeしてます．parameter.hというのは自分でつくってヘッダファイルで，キャリブレーション用のパラメターが書いてあります．
 
 
void macro_calb_cs(){  //Last modified 2012/01/09
  //difference from macro_calb_oka:
  //fit by Gaussian but not draw it,
  //and fit again by "Gaussian + linear" and draw it
  //for 137Cs
  char filename[99];//root file name
 
  sprintf(filename,"137cs-20111201.001.root");  
 
 
//filenameという変数を宣言して，そこにファイル名を代入してます．この変数はファイルを開くのと，あとでPDFファイルをつくるときに使います．sprintfはprintfと似た関数で，変数に文字列をいれることができます．ROOTのマクロを書くときにはけっこう便利です．あと、文字列をいれるときは，charというクラスでその文字数以上の配列として宣言するみたいです。
 
 
  TFile *f1=new TFile(filename);//opening root file
//ファイルを開いています。
 
  Double_t p[9][9][9];
//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です。
 
 
  TH1F *h[99];
  h[1]=h101;
  h[2]=h102;
  h[3]=h103;
  h[4]=h104;
  h[5]=h105;
  h[6]=h106;
//TH1Fは1次のFloatのHistgramのクラスです。ROOTファイルに入っているh101～h106までのヒストグラムをプログラム中で扱いやすいようにh[i]という配列を宣言して、そこに代入しています。
 
  // end of front matter
 
 
  //gStyle->SetOptFit(0001);
  gStyle->SetOptFit(0000);
//ここはフィッティングのパラメータをどのようにキャンバスに出力するか決めています  
 
 
 
  TCanvas *c1 = new TCanvas("c1","Raw Data and Fitting");
  c1->Divide(2,3);
 
  TCanvas *ctest = new TCanvas("ctest","Fitting Only With Gaussian");
  ctest->Divide(2,3);
 
 
 
  for(Int_t ch=1;ch<=6;ch++){
 
    char hclone[256]; 
    sprintf(hclone,"hclone%d",ch);
    TH1F *hoge=h[ch]->Clone(hclone);
    h[ch+10]=hoge;
    /* 
    sprintf(hclone,"hclone2%d",ch);
    TH1F *hoge2=h[ch]->Clone(hclone);
    h[ch+20]=hoge2;
    //h101-h106を違う目的で使いたいので複製しています
    */  
  }
 
 
 
  Double_t sigma1[3],sigma2[3];
  Double_t sigmae1[3],sigmae2[3];
  Double_t ene[3];  
  ene[0]=510.9989;
  ene[1]=1274.58;
  ene[2]=ene[0]+ene[1];
 
 
  for(Int_t ch=1;ch<=6;ch++){
 
    ctest->cd(ch)->SetLogy();    
 
    char fglname[99];
    Double_t paratest[5];
    Double_t ex[3];
 
    h[ch+10]->Draw();
    TLegend *l3 = new TLegend(0.6255,0.75,0.875,0.45);
 
    TF1 *ftest=new TF1("ftest","gaus",p1start[ch],p1end[ch]);
    h[ch+10]->Fit("ftest","","",p1start[ch],p1end[ch]);    
 
    Double_t parag[9];        
    ftest->GetParameters(parag);
 
    c1->cd(ch)->SetLogy();         
 
    TF1 *fit4;
    sprintf(fglname,"fgl1%d",ch);
    fit4 = new TF1(fglname,"[0]/sqrt(2.0*3.141592)/[1]*exp(-(x-[2])*(x-[2])/2.0/[1]/[1])+[3]*x+[4]");//[0]:constant ,[1]:sigma, [2]:mean
    fit4->SetParameters(parag[0],parag[2],parag[1],0,0);//パラメータの初期値の設定
    h[ch]->Fit(fglname,"","",p1start[ch],p1end[ch]);
    fit4 = h[ch]->GetFunction(fglname);
    fit4->GetParameters(paratest);// parameterをgetしてこれを使えば直線の式が
 
    Double_t Esigma;   
    Esigma[0]=fit4->GetParError(1);
 
    TF1 *flin1= new TF1("flin1","[0]*x+[1]",p1start[ch],p1end[ch]);
    flin1->FixParameter(0,paratest[3]);
    flin1->FixParameter(1,paratest[4]);//"fix" is proper, not "set"
    h[ch]->Fit(flin1,"Q0+","",p1start[ch],p1end[ch]);
    flin1=h[ch]->GetFunction("flin1");
 
    char caption[32]; 
    sprintf(caption,"mean=%lf(ch)",paratest[2]);
    l3->AddEntry((TObject*)0, caption,"");
    sprintf(caption,"sigma=%lf(ch)",paratest[1]);
    l3->AddEntry((TObject*)0, caption,"");
    sprintf(caption,"mean=%lf(keV)",paratest[2]*para1[ch]+para2[ch]);
    l3->AddEntry((TObject*)0, caption,"");
    sprintf(caption,"sigma=%lf(keV)",paratest[1]*para1[ch]);
    l3->AddEntry((TObject*)0, caption,"");
    sprintf(caption,"Errsig=%lf(keV)",Esigma*para1[ch]);
    l3->AddEntry((TObject*)0, caption,"");
 
 
 
    h[ch]->Draw();
    fit4->Draw("same");  
    flin1->Draw("same");		 
 
 
    l3->SetFillColor(0);// l の背景を白にします
    l3->Draw();//lをキャンバスに描きます
 
 
  }
 
 
 
  //キャンバスをPDFに出力します
  //Print
  /*  char pdfstart[99];
  sprintf(pdfstart,"%s.pdf(",filename);
  char pdfmiddle[99];
  sprintf(pdfmiddle,"%s.pdf",filename);
  char pdfend[99];
  sprintf(pdfend,"%s.pdf)",filename);
  
  c1->Print(pdfstart);
  c1->Print(pdfmiddle);
  //  c3->Print(pdfmiddle);
  c1->Print(pdfend);
  */
 
}