BOSS 7.0.7
BESIII Offline Software System
Loading...
Searching...
No Matches
wireFit.cxx
Go to the documentation of this file.
1const double m_norm = 500.734;
2const int WireNo = 6796;
3
4void wireFit()
5{
6 TFile *f=new TFile("wiregain.root");
7
8 double mean[WireNo];
9 double sigma[WireNo];
10 double chisq[WireNo];
11 double entryNo[WireNo];
12 double wireg[WireNo];
13 double wire[WireNo];
14 double fitmean[WireNo];
15 double fitmeanerr[WireNo];
16 double fitsigma[WireNo];
17 wiregcalib -> SetBranchAddress("wireg",wireg);
18 wiregcalib -> SetBranchAddress("mean", mean);
19 wiregcalib -> SetBranchAddress("sigma", sigma);
20 wiregcalib -> SetBranchAddress("chisq", chisq);
21 wiregcalib -> SetBranchAddress("entryNo", entryNo);
22 wiregcalib -> SetBranchAddress("wire", wire);
23 wiregcalib -> SetBranchAddress("fitmean", fitmean);
24 wiregcalib -> SetBranchAddress("fitmeanerr", fitmeanerr);
25 wiregcalib -> SetBranchAddress("fitsigma", fitsigma);
26 wiregcalib->GetEntry(0);
27
28 TH1F *h;
29 gStyle->SetOptFit(1111);
30 double gain;
31 ostringstream strout;
32
33 TFile* g = new TFile("wiregain_new.root","RECREATE");
34 TCanvas c1("c1", "canvas", 800, 642);
35 c1.Print("wiregain_new.root.ps[");
36 for(int i=0;i<WireNo;i++){
37// for(int i=241;i<242;i++){
38 strout.str("");
39 strout<<"dEdx_Wire_"<<i;
40 h=(TH1F*)f->Get(strout.str().c_str());
41 if(h->Integral()>100 && (chisq[i]>2 || chisq[i]<0.5)){
42 cout<<"Fit hisogram "<<strout.str().c_str()<<endl;
43 if(h->Integral((int)(fitmean[i]-1.0*fitsigma[i]-h->GetBinLowEdge(1))/h->GetBinWidth(1),(int)(fitmean[i]+1.0*fitsigma[i]-h->GetBinLowEdge(1))/h->GetBinWidth(1))>100) h->Fit("gaus","rq","",fitmean[i]-1.0*(fitsigma[i]),fitmean[i]+1.0*(fitsigma[i]));
44 chisq[i] = h->GetFunction("gaus")->GetChisquare()/(h->GetFunction("gaus")->GetNDF());
45 if(chisq[i]>2){
46 fitmean[i] = h->GetBinCenter(h->GetMaximumBin());
47 fitsigma[i]= 0;
48 }
49 else{
50 fitmean[i] = h->GetFunction("gaus")->GetParameter(1);
51 fitsigma[i] = h->GetFunction("gaus")->GetParameter(2);
52 }
53 wireg[i] = fitmean[i]/m_norm;
54 }
55 h->Draw();
56 c1.Print("wiregain_new.root.ps");
57 g->cd();
58 h->Write();
59 }
60 c1.Print("wiregain_new.root.ps]");
61
62 TTree* wiregain = new TTree("wiregcalib", "wiregcalib");
63 wiregain -> Branch("wireg", wireg, "wireg[6796]/D");
64 wiregain -> Branch("mean", mean, "mean[6796]/D");
65 wiregain -> Branch("sigma", sigma, "sigma[6796]/D");
66 wiregain -> Branch("chisq", chisq, "chisq[6796]/D");
67 wiregain -> Branch("entryNo", entryNo, "entryNo[6796]/D");
68 wiregain -> Branch("wire", wire, "wire[6796]/D");
69 wiregain -> Branch("fitmean", fitmean, "fitmean[6796]/D");
70 wiregain -> Branch("fitsigma", fitsigma, "fitsigma[6796]/D");
71
72 wiregain->Fill();
73 g->cd();
74 wiregain->Write();
75 g->Close();
76
77 TFile *f1 = new TFile("wiregain.root");
78 TTree *t1 = (TTree*)f1->Get("wiregcalib");
79 TFile *f2 = new TFile("wiregain_new.root");
80 TTree *t2 = (TTree*)f2->Get("wiregcalib");
81 Double_t x[6796], g_old[6796], g_new[6796], y[6796];
82 t1->SetBranchAddress("wire", x);
83 t1->SetBranchAddress("wireg", g_old);
84 t1->GetEntry(0);
85 t2->SetBranchAddress("wireg", g_new);
86 t2->GetEntry(0);
87 Double_t m_add = 0;
88 int m_count = 0;
89 for(int i=0; i<6796; i++){
90 if(g_old[i]>0.3){
91 y[i] = g_new[i]/g_old[i];
92 if(i<600) { m_add + = y[i]; m_count ++; }
93 }
94 else{
95 y[i] = 1;
96 cout << "i " << i << " wireg " << g_old[i] << endl;
97 }
98 }
99 cout << "ave diff before 600 wire: " << m_add/m_count << endl;
100
101 TGraph *gra = new TGraph(6796, x,y);
102 gra->SetTitle("ratio of wireg_new/wireg_old");
103 TCanvas *c2 = new TCanvas("wiregain", "ratio of wireg_new/wireg_old", 10,10,800,642);
104 gra->Draw("ap");
105 c2->SaveAs("compare_fit.eps");
106}
107
TTree * sigma
curve Branch("CurveSize",&CurveSize,"CurveSize/I")
double y[1000]
double x[1000]
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
curve SetBranchAddress("CurveSize",&CurveSize)
void wireFit()
Definition: wireFit.cxx:4
const int WireNo
Definition: wireFit.cxx:2
const double m_norm
Definition: wireFit.cxx:1