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

Last change on this file since 12762 was 12754, checked in by Jens Buss, 13 years ago
deleted bug in Search for crazy pixels
File size: 19.8 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/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 = false,
104 bool debug = false,
105 int verbosityLevel = 4)
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 if( showHistos ){
328 if (pixel%60==0){
329 cGraphs->cd(1);
330 gr2->Draw("APL");
331 cGraphs->cd(2);
332 gr1->Draw("APL");
333 }
334
335 cGainFit->Modified();
336 cGainFit->Update();
337 cGraphs->Modified();
338 cGraphs->Update();
339 }
340
341 //-----------------------------------------------------------------------------
342 // Debug Modus
343 //-----------------------------------------------------------------------------
344 if(debug){
345
346 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
347 timer.TurnOn();
348 TString input = Getline("Type 'q' to exit, <return> to go on: ");
349 timer.TurnOff();
350 if (input=="q\n") {
351 break;
352 }
353 }
354 }
355//-----------------------------------------------------------------------------
356// End of Loop
357//-----------------------------------------------------------------------------
358
359
360// ( cGainFit->cd(4);
361// hPixelAmplSpek[ NumberOfPixels ]->Draw();
362
363 cGainFit->cd(2);
364 hAmplDistribution->Draw();
365
366 cGainFit->cd(5);
367 hAmplDistribution2->Draw();
368
369 cGainFit->cd(6);
370 hGainFitBsl->Draw();
371 hbaseline->Draw("same");
372
373 cGainFit->Modified();
374 cGainFit->Update();
375
376 cGraphs->cd(1);
377 gr2->Draw("APL");
378 cGraphs->cd(2);
379 gr1->Draw("APL");
380 hList.Add(gr1);
381 hList.Add(gr2);
382 cGraphs->Modified();
383 cGraphs->Update();
384
385//-----------------------------------------------------------------------------
386// Write into logfile
387//-----------------------------------------------------------------------------
388 string filename = InputRootFileName;
389 size_t lastslash = filename.find_last_of('/');
390 filename = filename.substr(lastslash+1);
391 size_t firstdot = filename.find_first_of('.');
392 filename = filename.substr(0, firstdot);
393
394 string::iterator it;
395 for (it=filename.begin(); it < filename.end(); ){
396
397 if (!isdigit(*it)){
398 it=filename.erase(it);
399 }
400 else
401 ++it;
402 }
403if (verbosityLevel > 0) cout << filename << endl;
404
405 ofstream out;
406 out.open( OutTextFileName );
407 out << "pixelID/I:gainSD/F:gainDT:filename/l" << endl;
408 out << "# ---- the 1st line helps to use TTree::ReadFile, for reading again this txt file with root. " << endl;
409 out << "#" << InputRootFileName << endl;
410 out << "# PixelID : gain between Single and Double Peak : gain between Double and Tripple Peak " << endl;
411 out << "# gain is derived from the parameters of a three gaussian fit function " << endl;
412 for (UInt_t pixel=0; pixel<NumberOfPixels; pixel++){
413 out << channel[ pixel ].PID << ";\t" ;
414 out << channel[ pixel ].ampl_fspek_fit << ";\t" ;
415 // out << channel[ pixel ].amplDT << ";\t" ;
416 out << channel[ pixel ].ampl_fspek_fit2 << "\t" ;
417
418 //check if pixel gain is compareable to gain of the others and save result to logfile
419 channel[ pixel ].crazy = SearchAnormal(pixel, MAX_SIGMA, verbosityLevel);
420 if (channel[ pixel ].crazy) hListCrazy.Add(hPixelAmplSpek[ pixel ]);
421 out << channel[ pixel ].crazy << "\t" ;
422
423 out << filename << endl;
424
425 }
426 out.close();
427SaveHistograms( RootOutFileName );
428return( 0 );
429}
430
431//-----------------------------------------------------------------------------
432// End of Main
433//-----------------------------------------------------------------------------
434
435
436
437
438//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
439//+ +
440//+ Function definitions +
441//+ +
442//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
443
444
445
446//-----------------------------------------------------------------------------
447// Function: Book Histograms
448//-----------------------------------------------------------------------------
449
450void BookHistos(){
451
452 //Amplification Distribution - Single2Double
453 hAmplDistribution = new TH1F("hAmplDistribution","Distribution of G-APD Amplification (with Baseline Parameter)", 500, 0 , 20 );
454 hAmplDistribution->SetXTitle("Amplification [mV]");
455 hAmplDistribution->SetYTitle("Counts");
456 hAmplDistribution->SetAxisRange(7.5, 12.5, "X");
457 hList.Add( hAmplDistribution );
458
459 //Amplification Distribution2 - Double2Tripple
460 hAmplDistribution2 = new TH1F("hAmplDistribution2","Distribution of G-APD Amplification(without Baseline Parameter)", 500, 0 , 20 );
461 hAmplDistribution2->SetXTitle("Amplification [mV]");
462 hAmplDistribution2->SetYTitle("Counts");
463 hAmplDistribution2->SetAxisRange(7.5, 12.5, "X");
464 hList.Add( hAmplDistribution2 );
465
466 //cross vs Channels
467 hGainFitBsl = new TH1F("hGainFitBsl","BaselineDistribution from plot", 400,-99.5,100.5);
468 hGainFitBsl->SetYTitle("Counts");
469 hGainFitBsl->SetXTitle("Voltage [mV]");
470 hGainFitBsl->SetLineColor(2);
471 hList.Add( hGainFitBsl );
472
473}
474
475//-----------------------------------------------------------------------------
476// Save Histograms
477//-----------------------------------------------------------------------------
478
479void SaveHistograms( const char * loc_fname){
480 // create the histogram file (replace if already existing)
481 TFile tf( loc_fname, "RECREATE");
482
483 hList.Write(); // write the major histograms into the top level directory
484 tf.mkdir("PeakSpektra"); tf.cd("PeakSpektra"); // go to new subdirectory
485 hListPeakSpektra.Write(); // write histos into subdirectory
486 tf.mkdir("CrazyPixels"); tf.cd("CrazyPixels"); // go to new subdirectory
487 hListCrazy.Write(); // write histos into subdirector
488 tf.Close(); // close the file
489}
490
491//-----------------------------------------------------------------------------
492// Fit-Function: Sum of several Gauß-Function with falling Maximum Amplitude
493//-----------------------------------------------------------------------------
494
495Double_t spek_fit(Double_t *x, Double_t *par){
496 Double_t arg = 0;
497 Double_t arg2 = 0;
498 Double_t arg3 = 0;
499 Double_t arg4 = 0;
500 Double_t cross = 0;
501 if (par[0] != 0) cross = par[3]/par[0];
502 if (par[2] != 0) arg = (x[0] - par[1]-par[4])/par[2];
503 if (par[2] != 0) arg2 = (x[0] - 2*par[1]-par[4])/par[2];
504 if (par[2] != 0) arg3 = (x[0] - 3*par[1]-par[4])/par[2];
505 if (par[2] != 0) arg4 = (x[0] - 4*par[1]-par[4])/par[2];
506 Double_t peak1 = par[0]*TMath::Exp(-0.5 * arg * arg);
507 Double_t peak2 = cross*par[0]*TMath::Exp(-0.5 * arg2 * arg2 /2);
508 Double_t peak3 = cross*cross*par[0]*TMath::Exp(-0.5 * arg3 * arg3 /3);
509 Double_t peak4 = cross*cross*cross*par[0]*TMath::Exp(-0.5 * arg4 * arg4 /4);
510 par[5]=cross;
511 return peak1 + peak2 + peak3 + peak4;
512}
513
514//-----------------------------------------------------------------------------
515// Fit-Function: Sum of several Gauß-Function with overlying potential decay
516//-----------------------------------------------------------------------------
517
518Double_t spek_fit2(Double_t *x, Double_t *par2){
519
520 Double_t arg = 0;
521 Double_t arg2 = 0;
522 Double_t arg3 = 0;
523 Double_t arg4 = 0;
524 Double_t cross = 0;
525 if (par2[0] != 0) cross = par2[3]/par2[0];
526 if (par2[2] != 0) arg = (x[0] - par2[1])/par2[2];
527 if (par2[2] != 0) arg2 = (x[0] - 2*par2[1])/par2[2];
528 if (par2[2] != 0) arg3 = (x[0] - 3*par2[1])/par2[2];
529 if (par2[2] != 0) arg4 = (x[0] - 4*par2[1])/par2[2];
530 Double_t peak1 = par2[0]*TMath::Exp(-0.5 * arg * arg);
531 Double_t peak2 = cross*par2[0]*TMath::Exp(-0.5 * arg2 * arg2 /2);
532 Double_t peak3 = cross*cross*par2[0]*TMath::Exp(-0.5 * arg3 * arg3 /3);
533 Double_t peak4 = cross*cross*cross*par2[0]*TMath::Exp(-0.5 * arg4 * arg4 /4);
534 return peak1 + peak2 + peak3 + peak4;
535}
536
537//-----------------------------------------------------------------------------
538// Function: Search for Pixel that do not have a compareable gain to the Other
539//-----------------------------------------------------------------------------
540
541bool SearchAnormal(UInt_t pixel, Double_t distance, int verbosityLevel){
542 Double_t check = channel[ pixel ].ampl_fspek_fit;
543 // check = TMath::Sqrt(check*check);
544 if ( check > ( distribution.mean - distance * distribution.sigma )
545 && check < ( distribution.mean + distance * distribution.sigma )){
546 return false;
547 }
548 else{
549 if (verbosityLevel > 2) cout << "Not normal behaviour in Pixel: " << pixel << endl;
550 return true;
551 }
552}
553
Note: See TracBrowser for help on using the repository browser.