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

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