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

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