wireFit.cxx

Go to the documentation of this file.

const double m_norm = 500.734;
const int WireNo = 6796;
 
void wireFit()
{
    TFile *f=new TFile("wiregain.root");
 
    double mean[WireNo];
    double sigma[WireNo];
    double chisq[WireNo];
    double entryNo[WireNo];
    double wireg[WireNo];
    double wire[WireNo];
    double fitmean[WireNo];
    double fitmeanerr[WireNo];
    double fitsigma[WireNo];
    wiregcalib -> SetBranchAddress("wireg",wireg);
    wiregcalib -> SetBranchAddress("mean", mean);
    wiregcalib -> SetBranchAddress("sigma", sigma);
    wiregcalib -> SetBranchAddress("chisq", chisq);
    wiregcalib -> SetBranchAddress("entryNo", entryNo);
    wiregcalib -> SetBranchAddress("wire", wire);
    wiregcalib -> SetBranchAddress("fitmean", fitmean);
    wiregcalib -> SetBranchAddress("fitmeanerr", fitmeanerr);   
    wiregcalib -> SetBranchAddress("fitsigma", fitsigma); 
    wiregcalib->GetEntry(0);
 
    TH1F *h;
    gStyle->SetOptFit(1111);
    double gain;
    ostringstream strout;
 
    TFile* g = new TFile("wiregain_new.root","RECREATE");
    TCanvas c1("c1", "canvas", 800, 642);
    c1.Print("wiregain_new.root.ps[");
    for(int i=0;i<WireNo;i++){
//    for(int i=241;i<242;i++){
      strout.str("");
      strout<<"dEdx_Wire_"<<i;
      h=(TH1F*)f->Get(strout.str().c_str());
      if(h->Integral()>100 && (chisq[i]>2 || chisq[i]<0.5)){
    cout<<"Fit hisogram "<<strout.str().c_str()<<endl;
    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]));
    chisq[i] = h->GetFunction("gaus")->GetChisquare()/(h->GetFunction("gaus")->GetNDF());
    if(chisq[i]>2){
      fitmean[i] = h->GetBinCenter(h->GetMaximumBin());
      fitsigma[i]= 0;
    }
    else{
      fitmean[i] = h->GetFunction("gaus")->GetParameter(1);
      fitsigma[i] = h->GetFunction("gaus")->GetParameter(2);
    }
    wireg[i] = fitmean[i]/m_norm;
      }
      h->Draw();
      c1.Print("wiregain_new.root.ps");                                                                     
      g->cd();
      h->Write();
    }                      
    c1.Print("wiregain_new.root.ps]");
    
    TTree* wiregain = new TTree("wiregcalib", "wiregcalib");
    wiregain -> Branch("wireg", wireg, "wireg[6796]/D");
    wiregain -> Branch("mean", mean, "mean[6796]/D");
    wiregain -> Branch("sigma", sigma, "sigma[6796]/D");
    wiregain -> Branch("chisq", chisq, "chisq[6796]/D");
    wiregain -> Branch("entryNo", entryNo, "entryNo[6796]/D");
    wiregain -> Branch("wire", wire, "wire[6796]/D");
    wiregain -> Branch("fitmean", fitmean, "fitmean[6796]/D");
    wiregain -> Branch("fitsigma", fitsigma, "fitsigma[6796]/D");
 
    wiregain->Fill();
    g->cd();
    wiregain->Write();
    g->Close();
 
    TFile *f1 = new TFile("wiregain.root");
    TTree *t1 = (TTree*)f1->Get("wiregcalib");
    TFile *f2 = new TFile("wiregain_new.root");
    TTree *t2 = (TTree*)f2->Get("wiregcalib");
    Double_t x[6796], g_old[6796], g_new[6796], y[6796];
    t1->SetBranchAddress("wire", x);
    t1->SetBranchAddress("wireg", g_old);
    t1->GetEntry(0);
    t2->SetBranchAddress("wireg", g_new);
    t2->GetEntry(0);
    Double_t m_add = 0;
    int m_count = 0;
    for(int i=0; i<6796; i++){
      if(g_old[i]>0.3){
    y[i] = g_new[i]/g_old[i];
    if(i<600) { m_add + = y[i]; m_count ++; }
      }
      else{
      y[i] = 1;
      cout << "i " << i << "  wireg " << g_old[i] << endl;
      }
    }
    cout << "ave diff before 600 wire: " << m_add/m_count << endl;
   
    TGraph *gra = new TGraph(6796, x,y);
    gra->SetTitle("ratio of wireg_new/wireg_old");
    TCanvas *c2 = new TCanvas("wiregain", "ratio of wireg_new/wireg_old", 10,10,800,642);
    gra->Draw("ap");
    c2->SaveAs("compare_fit.eps");
}