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

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