| 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>
|
|---|
| 40 | using namespace std;
|
|---|
| 41 |
|
|---|
| 42 | //-----------------------------------------------------------------------------
|
|---|
| 43 | // Decleration of variables
|
|---|
| 44 | //-----------------------------------------------------------------------------
|
|---|
| 45 | UInt_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;
|
|---|
| 58 | vector<proj> channel;
|
|---|
| 59 |
|
|---|
| 60 | typedef struct{
|
|---|
| 61 | Double_t constant;
|
|---|
| 62 | Double_t mean;
|
|---|
| 63 | Double_t sigma;
|
|---|
| 64 | } dist_t;
|
|---|
| 65 | dist_t distribution;
|
|---|
| 66 |
|
|---|
| 67 | //Verbosity for Histograms
|
|---|
| 68 | TString verbos[4]={"0Q","Q", "", "V"};
|
|---|
| 69 | //-----------------------------------------------------------------------------
|
|---|
| 70 | // Decleration of Histograms
|
|---|
| 71 | //-----------------------------------------------------------------------------
|
|---|
| 72 |
|
|---|
| 73 | //Amplification Distribution
|
|---|
| 74 | TH1F *hAmplDistribution = NULL;
|
|---|
| 75 | TH1F *hAmplDistribution2 = NULL;
|
|---|
| 76 | TH1F *hGainFitBsl = NULL;
|
|---|
| 77 |
|
|---|
| 78 | //Lists of Histograms that have to be saved
|
|---|
| 79 | TObjArray hList;
|
|---|
| 80 | TObjArray hListPeakSpektra;
|
|---|
| 81 | TObjArray hListCrazy;
|
|---|
| 82 |
|
|---|
| 83 | //-----------------------------------------------------------------------------
|
|---|
| 84 | // Functions
|
|---|
| 85 | //-----------------------------------------------------------------------------
|
|---|
| 86 | Double_t spek_fit(Double_t *, Double_t *);
|
|---|
| 87 | Double_t spek_fit2(Double_t *, Double_t *);
|
|---|
| 88 |
|
|---|
| 89 | void BookHistos();
|
|---|
| 90 | void SaveHistograms(const char * );
|
|---|
| 91 | bool SearchAnormal(UInt_t, Double_t, int verbosityLevel );
|
|---|
| 92 |
|
|---|
| 93 | //-----------------------------------------------------------------------------
|
|---|
| 94 | // Main function
|
|---|
| 95 | //-----------------------------------------------------------------------------
|
|---|
| 96 |
|
|---|
| 97 |
|
|---|
| 98 | int 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 | }
|
|---|
| 408 | if (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();
|
|---|
| 432 | SaveHistograms( RootOutFileName );
|
|---|
| 433 | return( 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 |
|
|---|
| 455 | void 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 |
|
|---|
| 484 | void 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 |
|
|---|
| 505 | Double_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 |
|
|---|
| 528 | Double_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 |
|
|---|
| 551 | bool 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 |
|
|---|
