| 1 | #include <TFile.h>
|
|---|
| 2 | #include <TH2F.h>
|
|---|
| 3 | #include <TH1D.h>
|
|---|
| 4 | #include <TROOT.h>
|
|---|
| 5 | #include <TF1.h>
|
|---|
| 6 | #include <TTimer.h>
|
|---|
| 7 | #include <TString.h>
|
|---|
| 8 | #include <Getline.h>
|
|---|
| 9 | #include <TCanvas.h>
|
|---|
| 10 | #include <TGraph.h>
|
|---|
| 11 | #include <TMultiGraph.h>
|
|---|
| 12 |
|
|---|
| 13 | #include <iostream>
|
|---|
| 14 | #include <fstream>
|
|---|
| 15 | using namespace std;
|
|---|
| 16 |
|
|---|
| 17 | int gainfit(
|
|---|
| 18 | const char * InputRootFileName,
|
|---|
| 19 | const char * OutTextFileName,
|
|---|
| 20 | bool debug = false
|
|---|
| 21 | ) {
|
|---|
| 22 |
|
|---|
| 23 | TFile * tf = TFile::Open( InputRootFileName );
|
|---|
| 24 | TH2F *h = (TH2F*)gROOT->FindObject("hAmplSpek_cfd");
|
|---|
| 25 |
|
|---|
| 26 | TH1D *p1 = NULL;
|
|---|
| 27 | TH1D *p2 = NULL;
|
|---|
| 28 |
|
|---|
| 29 | TF1 *f_peaks = new TF1("total","[2]*TMath::Exp(-(x-[0])*(x-[0])/(2*[3]*[3])) "
|
|---|
| 30 | "+ [4]*TMath::Exp(-(x-[0]-[1])*(x-[0]-[1])/(2*[5]*[5])) "
|
|---|
| 31 | "+ [6]*TMath::Exp((-x-[7])/[8])", 3, 23);
|
|---|
| 32 | Double_t par_init[9] = {8, 8, 1000 , 1 , 190 , 2, 30 , 5 , 10};
|
|---|
| 33 | f_peaks->SetParameters(par_init);
|
|---|
| 34 |
|
|---|
| 35 | Int_t n = 1440;
|
|---|
| 36 | Double_t simple_gain[n];
|
|---|
| 37 | Double_t gain[n];
|
|---|
| 38 | Double_t pixelvec[n];
|
|---|
| 39 | for (int i=0; i<1440; i++){
|
|---|
| 40 | pixelvec[i]=i;
|
|---|
| 41 | simple_gain[i]=0;
|
|---|
| 42 | gain[i]=0;
|
|---|
| 43 | }
|
|---|
| 44 | TGraph *gr1 = new TGraph(n,pixelvec,simple_gain);
|
|---|
| 45 | TGraph *gr2 = new TGraph(n,pixelvec,gain);
|
|---|
| 46 |
|
|---|
| 47 | gr1->SetMarkerStyle(5);
|
|---|
| 48 | gr2->SetMarkerStyle(2);
|
|---|
| 49 | gr1->SetMarkerColor(kRed);
|
|---|
| 50 | gr2->SetMarkerColor(kBlack);
|
|---|
| 51 | gr1->SetMarkerSize(0.3);
|
|---|
| 52 | gr2->SetMarkerSize(0.3);
|
|---|
| 53 | TMultiGraph *mg = new TMultiGraph();
|
|---|
| 54 | mg->Add(gr1);
|
|---|
| 55 | mg->Add(gr2);
|
|---|
| 56 |
|
|---|
| 57 |
|
|---|
| 58 | TCanvas *cGainFit = new TCanvas("cGainFit", "Fit of Amplitude Spektra", 0,0, 1000,500);
|
|---|
| 59 | cGainFit->Divide(1,2);
|
|---|
| 60 |
|
|---|
| 61 | // bad conding ... 1440 shouldn be hardcoded here
|
|---|
| 62 | for ( int pixelid =0; pixelid < 1440; pixelid++){
|
|---|
| 63 | cout << "PID: " << pixelid << endl;
|
|---|
| 64 | p1 = h->ProjectionY("p1",pixelid+1, pixelid+1);
|
|---|
| 65 | p2 = (TH1D*) p1->Clone("p2");
|
|---|
| 66 |
|
|---|
| 67 |
|
|---|
| 68 | cGainFit->cd(1);
|
|---|
| 69 | p1->GetXaxis()->SetRangeUser( 5, 25);
|
|---|
| 70 | p1->Draw();
|
|---|
| 71 |
|
|---|
| 72 | Double_t peak[2];
|
|---|
| 73 | Double_t rms[2];
|
|---|
| 74 |
|
|---|
| 75 | peak[0]= p1->GetBinCenter(p1->GetMaximumBin());
|
|---|
| 76 | rms[0] = p1->GetRMS();
|
|---|
| 77 | cout << "1st position:" << peak[0] << endl;
|
|---|
| 78 | cout << "1st rms:" << rms[0] << endl;
|
|---|
| 79 | f_peaks->SetParameters(par_init);
|
|---|
| 80 | p1->Fit(f_peaks, "RQ");
|
|---|
| 81 |
|
|---|
| 82 | p2->GetXaxis()->SetRangeUser( peak[0] + 1.5*rms[0] , 25);
|
|---|
| 83 |
|
|---|
| 84 | peak[1] = p2->GetBinCenter(p2->GetMaximumBin());
|
|---|
| 85 | rms[1] = p2->GetRMS();
|
|---|
| 86 | cout << "2nd position:" << peak[1] << endl;
|
|---|
| 87 | cout << "2nd rms:" << rms[1] << endl;
|
|---|
| 88 | cout << "easy gain" << peak[1]-peak[0] << "\t\t";
|
|---|
| 89 | cout << "fittet gain" << f_peaks->GetParameter(1) << endl;
|
|---|
| 90 |
|
|---|
| 91 | simple_gain[pixelid]=peak[1]-peak[0];
|
|---|
| 92 | gain[pixelid] = f_peaks->GetParameter(1);
|
|---|
| 93 | gr1->SetPoint(pixelid, pixelvec[pixelid], simple_gain[pixelid]);
|
|---|
| 94 | gr2->SetPoint(pixelid, pixelvec[pixelid], gain[pixelid]);
|
|---|
| 95 |
|
|---|
| 96 | cGainFit->cd(2);
|
|---|
| 97 | mg->Draw("AP");
|
|---|
| 98 | if (pixelid % 50 ==0){
|
|---|
| 99 | cGainFit->Modified();
|
|---|
| 100 | cGainFit->Update();
|
|---|
| 101 | }
|
|---|
| 102 | if(debug){
|
|---|
| 103 | cGainFit->Modified();
|
|---|
| 104 | cGainFit->Update();
|
|---|
| 105 | TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
|
|---|
| 106 | timer.TurnOn();
|
|---|
| 107 | TString input = Getline("Type 'q' to exit, <return> to go on: ");
|
|---|
| 108 | timer.TurnOff();
|
|---|
| 109 | if (input=="q\n") {
|
|---|
| 110 | break;
|
|---|
| 111 | }
|
|---|
| 112 | }
|
|---|
| 113 | }
|
|---|
| 114 |
|
|---|
| 115 | //tf->Close();
|
|---|
| 116 | string filename = InputRootFileName;
|
|---|
| 117 | size_t lastslash = filename.find_last_of('/');
|
|---|
| 118 | filename = filename.substr(lastslash+1);
|
|---|
| 119 | size_t firstdot = filename.find_first_of('.');
|
|---|
| 120 | filename = filename.substr(0, firstdot);
|
|---|
| 121 |
|
|---|
| 122 | string::iterator it;
|
|---|
| 123 | for (it=filename.begin(); it < filename.end(); ){
|
|---|
| 124 |
|
|---|
| 125 | if (!isdigit(*it)){
|
|---|
| 126 | it=filename.erase(it);
|
|---|
| 127 | }
|
|---|
| 128 | else
|
|---|
| 129 | ++it;
|
|---|
| 130 | }
|
|---|
| 131 | cout << filename << endl;
|
|---|
| 132 |
|
|---|
| 133 |
|
|---|
| 134 | ofstream out;
|
|---|
| 135 | out.open( OutTextFileName );
|
|---|
| 136 | out << "pixelID/I:gain/F:fitgain:filename/l" << endl;
|
|---|
| 137 | out << "# ---- the 1st line helps to use TTree::ReadFile, for reading again this txt file with root. " << endl;
|
|---|
| 138 | out << "#" << InputRootFileName << endl;
|
|---|
| 139 | out << "# PixelID : simple_gain : gain " << endl;
|
|---|
| 140 | out << "# simple_gain is just the distance of the local maxima in the amplitude spektra " << endl;
|
|---|
| 141 | out << "# gain is derived from a fit of two gaussians " << endl;
|
|---|
| 142 | for (int i=0; i<1440; i++){
|
|---|
| 143 | out << pixelvec[i] << "\t" << simple_gain[i] << "\t" << gain[i] << "\t" << filename << endl;
|
|---|
| 144 |
|
|---|
| 145 | }
|
|---|
| 146 | out.close();
|
|---|
| 147 |
|
|---|
| 148 | return 0;
|
|---|
| 149 | }
|
|---|
