source: fact/tools/rootmacros/gainfit2.C@ 12712

Last change on this file since 12712 was 12712, checked in by Jens Buss, 13 years ago
added differnt typee of fit function to see which fits best; detect pixel with not normal gain
File size: 17.4 KB
Line 
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 <TMath.h>
12#include <TMultiGraph.h>
13#include <TH1.h>
14#include <TProfile.h>
15#include <TStyle.h>
16
17#include <iostream>
18#include <fstream>
19using namespace std;
20
21//-----------------------------------------------------------------------------
22// Decleration of variables
23//-----------------------------------------------------------------------------
24UInt_t NumberOfPixels;
25
26 typedef struct{
27 UInt_t PID;
28 vector<float> fitParameters[9];
29 Double_t amplSD;
30 Double_t amplDT;
31 Double_t amplMean;
32 Double_t ampl_fspek_fit;
33 bool crazy;
34 } proj;
35vector<proj> channel;
36
37typedef struct{
38 Double_t constant;
39 Double_t mean;
40 Double_t sigma;
41} dist_t;
42dist_t distribution;
43
44//-----------------------------------------------------------------------------
45// Decleration of Histograms
46//-----------------------------------------------------------------------------
47
48//Amplification Distribution
49TH1F *hAmplDistribution = NULL;
50TH1F *hAmplDistribution2 = NULL;
51TH1F *hPixelGainFit = NULL;
52
53//Lists of Histograms that have to be saved
54TObjArray hList;
55TObjArray hListPeakSpektra;
56
57//-----------------------------------------------------------------------------
58// Functions
59//-----------------------------------------------------------------------------
60Double_t spek_fit(Double_t *, Double_t *);
61Double_t spek_fit2(Double_t *, Double_t *);
62
63void BookHistos();
64void SaveHistograms(const char * );
65bool SearchAnormal(UInt_t, Double_t );
66
67//-----------------------------------------------------------------------------
68// Main function
69//-----------------------------------------------------------------------------
70
71
72int gainfit2(
73 const char * InputRootFileName = "../analysis/fpeak/20111109_006/20111109_006_fpeak.root",
74 const char * InputBaselineFileName = "../analysis/fpeak/20111109_006/20111109_006_fbsl.root",
75 const char * OutTextFileName = "../analysis/fpeak/20111109_006_gainfit_test.txt",
76 const char * RootOutFileName = "../analysis/fpeak/20111109_006_gainfit_test.root",
77 bool showHistos = false,
78 bool debug = false)
79{
80
81//-----------------------------------------------------------------------------
82// Histograms, Canvases and Fit Functions
83//-----------------------------------------------------------------------------
84
85//Open Rootfiles Files
86 TFile * baseline = TFile::Open( InputBaselineFileName );
87 TFile * file = TFile::Open( InputRootFileName );
88
89//Baseline
90 TH2F * hbaseline = (TH2F*)baseline->Get("histo_mean");
91//Amplitude spectrum
92 TH2F *hAmplSpek = (TH2F*)file->Get("hAmplSpek_cfd");
93 NumberOfPixels = hAmplSpek->GetNbinsX();
94 channel.resize(hAmplSpek->GetNbinsX());
95
96//Amplitude Spectrum of a Pixel
97 TH1D* hPixelAmplSpek[1440]; ///++++ TODO: Warning hardcoded pixelnumber, change that +++
98
99//Book Histograms
100 BookHistos();
101
102//Define Fit Functions
103 Double_t par[9];
104 TF1 *g1 = new TF1("g1","gaus",5,15);
105 TF1 *g2 = new TF1("g2","gaus",15,25);
106 TF1 *g3 = new TF1("g3","gaus",25,35);
107 TF1 *g4 = new TF1("g4","gaus");
108 TF1 *total = new TF1("total","gaus(0)+gaus(3)+gaus(6)",7,32);
109 TF1 *fspek_fit = new TF1("fspek_fit", spek_fit, 5, 60, 5);
110 TF1 *fspek_fit2 = new TF1("fspek_fit2", spek_fit2, 5, 60, 5);
111
112//Create Canvases
113 TCanvas *cGainFit = new TCanvas("cGainFit", "Fit of Amplitude Spektra", 0,0, 1200,400);
114 cGainFit->Clear();
115 cGainFit->Divide(3, 2);
116
117 TCanvas *cGraphs = new TCanvas("cGraphs", "gain vs. channel", 0,401, 1200,400);
118 cGraphs->Clear();
119 cGraphs->Divide(1, 2);
120
121//Draw imported Histograms
122 cGainFit->cd(1);
123 gStyle->SetPalette(1,0);
124 gROOT->SetStyle("Plain");
125 hAmplSpek->SetTitle("Amplitude Spectrum of all Pixels");
126 hAmplSpek->SetAxisRange(0, 100, "Y");
127 hAmplSpek->Draw("COLZ");
128
129
130 cGainFit->cd(2);
131 hbaseline->SetTitle("Baselinedistribution of all Pixels");
132 hbaseline->SetAxisRange(-5, 5, "X");
133 //hbaseline->SetAxisRange(0, 100, "Y");
134 hbaseline->Fit("gaus");
135 hbaseline->Draw();
136
137 //-----------------------------------------------------------------------------
138 // Graphs
139 //-----------------------------------------------------------------------------
140
141 Double_t simple_gain[NumberOfPixels];
142 Double_t gain[NumberOfPixels];
143 Double_t pixelvec[NumberOfPixels];
144 for (UInt_t i=0; i<NumberOfPixels; i++){
145 pixelvec[i]=i;
146 simple_gain[i]=0;
147 gain[i]=0;
148 }
149
150 TGraph *gr1 = new TGraph(NumberOfPixels,pixelvec,simple_gain);
151 TGraph *gr2 = new TGraph(NumberOfPixels,pixelvec,gain);
152
153 gr1->SetMarkerStyle(5);
154 gr2->SetMarkerStyle(2);
155 gr1->SetMarkerColor(kRed);
156 gr2->SetMarkerColor(kBlack);
157 gr1->SetMarkerSize(0.5);
158 gr2->SetMarkerSize(0.5);
159 TMultiGraph *mg = new TMultiGraph();
160 mg->Add(gr1);
161 mg->Add(gr2);
162
163//-----------------------------------------------------------------------------
164// Loop over all Pixels and fit Distribution of singles, doubles and tripples
165//-----------------------------------------------------------------------------
166
167 cout << "Number of Pixels: " << NumberOfPixels << endl;
168 // Begin of Loop
169 for (UInt_t pixel = 0; pixel < NumberOfPixels; pixel++){
170
171 //Title
172 TString histotitle;
173 histotitle += "hPixelPro_";
174 histotitle += pixel + 1 ;
175
176 cout << "-----------------------------------------------------------" << endl;
177 cout << " PID: " << pixel+1 << endl;
178 cout << "-----------------------------------------------------------" << endl;
179 //Amplitude Spectrum of a Pixel
180 hPixelAmplSpek[ pixel ] = new TH1D;
181 hPixelAmplSpek[ pixel ] = hAmplSpek->ProjectionY( histotitle , pixel+1, pixel+1); //Projectipon of each Pixel out of Ampl.Spectrum
182 hPixelAmplSpek[ pixel ]->SetAxisRange(0, 40, "X");
183
184 histotitle = "Amplitude Spectrum of Pixel #";
185 histotitle += pixel + 1 ;
186
187 hPixelAmplSpek[ pixel ]->SetTitle(histotitle);
188 hPixelAmplSpek[ pixel ]->SetYTitle("Counts [a.u.]");
189 hListPeakSpektra.Add(hPixelAmplSpek[ pixel ]);
190
191
192//-----------------------------------------------------------------------------
193// Single Gaussian Fits and Tripple-Gaussian Fit Funciton
194//-----------------------------------------------------------------------------
195
196 total->SetLineColor(2);
197 fspek_fit->SetLineColor(3);
198 fspek_fit2->SetLineColor(kBlue);
199 fspek_fit2->SetLineStyle(2);
200
201 cGainFit->cd(3);
202 gPad->SetLogy(1);
203
204 hPixelAmplSpek[ pixel ]->Fit(g1,"R");
205 hPixelAmplSpek[ pixel ]->Fit(g2,"R+");
206 hPixelAmplSpek[ pixel ]->Fit(g3,"R+");
207 g1->GetParameters(&par[0]);
208 g2->GetParameters(&par[3]);
209 g3->GetParameters(&par[6]);
210
211 fspek_fit->SetParameters(par);
212 fspek_fit2->SetParameters(par);
213 total->SetParameters(par);
214 hPixelAmplSpek[ pixel ]->Fit(total,"R+");
215 hPixelAmplSpek[ pixel ]->Fit(fspek_fit,"R+");
216 hPixelAmplSpek[ pixel ]->Fit(fspek_fit2,"R+");
217
218 if( showHistos ){
219 hPixelAmplSpek[ pixel ]->Draw();
220 }
221
222
223
224//-----------------------------------------------------------------------------
225// Save Parameters into vetor
226//-----------------------------------------------------------------------------
227
228// total->GetParameters(channel[ pixel ].fitParameters);
229 channel[ pixel ].PID = pixel+1;
230 channel[ pixel ].amplSD = total->GetParameter(4) - total->GetParameter(1);
231 channel[ pixel ].amplDT = total->GetParameter(7) - total->GetParameter(4);
232 channel[ pixel ].amplMean = fspek_fit->GetParameter(1);
233 channel[ pixel ].ampl_fspek_fit = fspek_fit->GetParameter(1);
234
235// channel[ pixel ].tripple_peakSigm = total->GetParameter(8);
236// channel[ pixel ].ampl = total->GetParameter(1) - total->GetParameter(4);
237
238//-----------------------------------------------------------------------------
239// Draw Histograms
240//-----------------------------------------------------------------------------
241
242 cGainFit->cd(4);
243 gPad->SetLogy(1);
244 g4->SetLineColor(2);
245 hAmplDistribution->Fill( channel[ pixel ].amplSD ) ;
246 hAmplDistribution->Fit(g4);
247 distribution.sigma = g4->GetParameter(2);
248 distribution.mean = g4->GetParameter(1);
249 distribution.constant = g4->GetParameter(0);
250
251 cGainFit->cd(5);
252 gPad->SetLogy(1);
253 g4->SetLineColor(2);
254 hAmplDistribution2->Fill(channel[ pixel ].amplMean) ;
255 hAmplDistribution2->Fit(g4);
256
257 cGainFit->cd(6);
258 hPixelGainFit->SetBinContent(pixel, channel[ pixel ].amplSD);
259
260 if( showHistos ){
261 if (pixel%40==0){
262 cGainFit->cd(4);
263 hAmplDistribution->Draw();
264 cGainFit->cd(5);
265 hAmplDistribution2->Draw();
266 }
267
268 cGainFit->cd(6);
269 hPixelGainFit->Draw();
270 }
271
272//-----------------------------------------------------------------------------
273// Draw Graphs of Distributions
274//-----------------------------------------------------------------------------
275
276 gr1->SetPoint(pixel, channel[pixel].PID, channel[pixel].amplSD);
277 gr2->SetPoint(pixel, channel[pixel].PID, channel[pixel].amplMean);
278
279
280 if( showHistos ){
281 if (pixel%60==0){
282 cGraphs->cd(2);
283 gr2->Draw("APL");
284 cGraphs->cd(1);
285 gr1->Draw("APL");
286 }
287
288 cGainFit->Modified();
289 cGainFit->Update();
290 cGraphs->Modified();
291 cGraphs->Update();
292 }
293
294 //-----------------------------------------------------------------------------
295 // Debug Modus
296 //-----------------------------------------------------------------------------
297 if(debug){
298
299 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
300 timer.TurnOn();
301 TString input = Getline("Type 'q' to exit, <return> to go on: ");
302 timer.TurnOff();
303 if (input=="q\n") {
304 break;
305 }
306 }
307 }
308//-----------------------------------------------------------------------------
309// End of Loop
310//-----------------------------------------------------------------------------
311
312
313// ( cGainFit->cd(3);
314// hPixelAmplSpek[ NumberOfPixels ]->Draw();
315
316 cGainFit->cd(4);
317 hAmplDistribution->Draw();
318
319 cGainFit->cd(5);
320 hAmplDistribution2->Draw();
321
322 cGainFit->cd(6);
323 hPixelGainFit->Draw();
324
325 cGainFit->Modified();
326 cGainFit->Update();
327
328 cGraphs->cd(2);
329 gr2->Draw("APL");
330 cGraphs->cd(1);
331 gr1->Draw("APL");
332 cGraphs->Modified();
333 cGraphs->Update();
334
335//-----------------------------------------------------------------------------
336// Write into logfile
337//-----------------------------------------------------------------------------
338 string filename = InputRootFileName;
339 size_t lastslash = filename.find_last_of('/');
340 filename = filename.substr(lastslash+1);
341 size_t firstdot = filename.find_first_of('.');
342 filename = filename.substr(0, firstdot);
343
344 string::iterator it;
345 for (it=filename.begin(); it < filename.end(); ){
346
347 if (!isdigit(*it)){
348 it=filename.erase(it);
349 }
350 else
351 ++it;
352 }
353 cout << filename << endl;
354
355 ofstream out;
356 out.open( OutTextFileName );
357 out << "pixelID/I:gainSD/F:gainDT:filename/l" << endl;
358 out << "# ---- the 1st line helps to use TTree::ReadFile, for reading again this txt file with root. " << endl;
359 out << "#" << InputRootFileName << endl;
360 out << "# PixelID : gain between Single and Double Peak : gain between Double and Tripple Peak " << endl;
361 out << "# gain is derived from the parameters of a three gaussian fit function " << endl;
362 for (UInt_t pixel=0; pixel<NumberOfPixels; pixel++){
363 out << channel[ pixel ].PID << ";\t" ;
364 out << channel[ pixel ].amplSD << ";\t" ;
365 out << channel[ pixel ].amplDT << ";\t" ;
366 out << channel[ pixel ].ampl_fspek_fit << "\t" ;
367
368 channel[ pixel ].crazy = SearchAnormal(pixel, 4); //check if pixel had a strange gain
369 out << channel[ pixel ].crazy << "\t" ;
370
371 out << filename << endl;
372
373 }
374 out.close();
375SaveHistograms( RootOutFileName );
376return( 0 );
377}
378
379//-----------------------------------------------------------------------------
380// End of Main
381//-----------------------------------------------------------------------------
382
383
384
385
386//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
387//+ +
388//+ Function definitions +
389//+ +
390//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
391
392
393
394//-----------------------------------------------------------------------------
395// Function: Book Histograms
396//-----------------------------------------------------------------------------
397
398void BookHistos(){
399
400 //Amplification Distribution - Single2Double
401 hAmplDistribution = new TH1F("hAmplDistribution","Distribution of G-APD Amplification (1+2)", 500, 0 , 20 );
402 hAmplDistribution->SetXTitle("Amplification [mV]");
403 hAmplDistribution->SetYTitle("Counts");
404 hAmplDistribution->SetAxisRange(7.5, 12.5, "X");
405 hList.Add( hAmplDistribution );
406
407 //Amplification Distribution2 - Double2Tripple
408 hAmplDistribution2 = new TH1F("hAmplDistribution2","Distribution of G-APD Amplification(2+3)", 500, 0 , 20 );
409 hAmplDistribution2->SetXTitle("Amplification [mV]");
410 hAmplDistribution2->SetYTitle("Counts");
411 hAmplDistribution2->SetAxisRange(7.5, 12.5, "X");
412 hList.Add( hAmplDistribution2 );
413
414 //Amplification vs Channels
415 hPixelGainFit = new TH1F("hPixelGainFit","Amplifications (SD) of each Channel", 1440, -0.5 , 1439.5 );
416 hPixelGainFit->SetYTitle("Amplification [mV]");
417 hPixelGainFit->SetXTitle("Channels");
418 hPixelGainFit->SetLineColor(2);
419 hList.Add( hPixelGainFit );
420
421}
422
423//-----------------------------------------------------------------------------
424// Save Histograms
425//-----------------------------------------------------------------------------
426
427void SaveHistograms( const char * loc_fname){
428 // create the histogram file (replace if already existing)
429 TFile tf( loc_fname, "RECREATE");
430
431 hList.Write(); // write the major histograms into the top level directory
432 tf.mkdir("PeakSpektra"); tf.cd("PeakSpektra"); // go to new subdirectory
433 hListPeakSpektra.Write(); // write histos into subdirectory
434
435 tf.Close(); // close the file
436}
437
438//-----------------------------------------------------------------------------
439// MulitOrder MultiGauß Function
440//-----------------------------------------------------------------------------
441
442Double_t spek_fit(Double_t *x, Double_t *par){
443
444 Double_t arg = 0;
445 Double_t arg2 = 0;
446 Double_t arg3 = 0;
447 Double_t arg4 = 0;
448 Double_t cross = 1;
449 if (par[0] != 0) cross = par[4]/par[0];
450 if (par[2] != 0) arg = (x[0] - par[1])/par[2];
451 if (par[2] != 0) arg2 = (x[0] - 2*par[1])/par[2];
452 if (par[2] != 0) arg3 = (x[0] - 3*par[1])/par[2];
453 if (par[2] != 0) arg4 = (x[0] - 4*par[1])/par[2];
454 cross = par[4]/par[0];
455 Double_t peak1 = par[0]*TMath::Exp(-0.5 * arg * arg);
456 Double_t peak2 = cross*par[0]*TMath::Exp(-0.5 * arg2 * arg2 /2);
457 Double_t peak3 = cross*cross*par[0]*TMath::Exp(-0.5 * arg3 * arg3 /3);
458 Double_t peak4 = cross*cross*cross*par[0]*TMath::Exp(-0.5 * arg4 * arg4 /4);
459 //Double_t decayoverlay = TMath::Exp(-x[0]);
460 return peak1 + peak2 + peak3 + peak4;
461}
462
463Double_t spek_fit2(Double_t *x, Double_t *par){
464
465 Double_t arg = 0;
466 Double_t arg2 = 0;
467 Double_t arg3 = 0;
468 Double_t arg4 = 0;
469 Double_t cross = 1;
470 if (par[0] != 0) cross = par[4]/par[0];
471 if (par[2] != 0) arg = (x[0] - par[1])/par[2];
472 if (par[2] != 0) arg2 = (x[0] - 2*par[1])/par[2];
473 if (par[2] != 0) arg3 = (x[0] - 3*par[1])/par[2];
474 if (par[2] != 0) arg4 = (x[0] - 4*par[1])/par[2];
475 cross = par[4]/par[0];
476 Double_t peak1 = par[0]*TMath::Exp(-0.5 * arg * arg);
477 Double_t peak2 = cross*par[0]*TMath::Exp(-0.5 * arg2 * arg2 /2);
478 Double_t peak3 = cross*cross*par[0]*TMath::Exp(-0.5 * arg3 * arg3 /3);
479 Double_t peak4 = cross*cross*cross*par[0]*TMath::Exp(-0.5 * arg4 * arg4 /4);
480 Double_t decayoverlay = TMath::Exp(-x[0]*TMath::Log(1.0001));
481 return peak1 + peak2 + peak3 + peak4 + decayoverlay;
482}
483
484bool SearchAnormal(UInt_t pixel, Double_t distance){ //// MACHT ALLES KEINEN SINN WENN MAN KEINE BETRÄGE BETRACHTET
485 Double_t check = channel[ pixel ].ampl_fspek_fit;
486 check = TMath::Sqrt(check*check);
487 if (check < (distance * distribution.sigma + distribution.mean)){
488 return false;
489 }
490 else{
491 cout << "Not normal behaviotruer in Pixel: " << pixel +1 << endl;
492 return true;
493 }
494}
495
496/*
497
498
499
500
501
502
503 */
504
505
506/*
507Double_t spek_fit(Double_t *x, Double_t *par, Uint_t *maxorder)
508 {
509 Double_t exp = 0;
510
511 //Double_t arg2 = 0;
512 //if (par[2] != 0) arg = (v[0] - par[1])/par[2];
513 //if (par[5] != 0) arg2 = (v[0] - par[4])/par[5];
514
515 //Double_t fitval = par[0]*TMath::Exp(-0.5*arg*arg) + par[3]*TMath::Exp(-0.5*arg2*arg2);
516
517 for(int order = 0; order < maxorder; order++){
518 peak[order]=par[0]*par[4]/(order+1)*((x[0]-(order+1)*par[1])*(x[0]-(order+1)*par[1]))/(par[2]*par[2]);
519 exp = exp + peak[order];
520 }
521
522 Double_t fitval = TMath::Exp(-0.5*exp);
523 return fitval;
524 }
525*/
Note: See TracBrowser for help on using the repository browser.