From: Benjamin Linnik Date: Thu, 7 May 2015 10:04:13 +0000 (+0200) Subject: Analysis: few bugfixes, less warnings, new veto peak search algorithm, ROOT file... X-Git-Url: https://jspc29.x-matter.uni-frankfurt.de/git/?a=commitdiff_plain;h=2268c3460aa6709817db5a5acb6321ada1166236;p=radhard.git Analysis: few bugfixes, less warnings, new veto peak search algorithm, ROOT file is saved/wrote to disk prior to analysis --- diff --git a/MABS_run_analyzer/ChargeSpektrum.c b/MABS_run_analyzer/ChargeSpektrum.c index f397d89..51c756f 100644 --- a/MABS_run_analyzer/ChargeSpektrum.c +++ b/MABS_run_analyzer/ChargeSpektrum.c @@ -27,7 +27,7 @@ #include "CSVRow.C" #include -Int_t* ReadRunList(Int_t*); +Int_t* ReadRunList(); void plotAllRuns(); Run** runs; @@ -52,10 +52,10 @@ void ChargeSpektrum(Int_t runnumber = -1) { /// number of runs to be analyzed, number of lines read by @c ReadRunList() numberRuns=0; - ReadRunList(&numberRuns); + ReadRunList(); /// array with run numbers runList=new Int_t[numberRuns]; - runList=ReadRunList(&numberRuns); + runList=ReadRunList(); } runs = new Run*[numberRuns]; @@ -132,8 +132,7 @@ void ChargeSpektrum(Int_t runnumber = -1) } - -Int_t* ReadRunList(Int_t* numberRuns) +Int_t* ReadRunList() { Int_t* runList=new Int_t[1000]; @@ -152,7 +151,7 @@ Int_t* ReadRunList(Int_t* numberRuns) cout << "File ended"; } } - *numberRuns=runLauf; + numberRuns=runLauf; return runList; } diff --git a/MABS_run_analyzer/MAPS.c b/MABS_run_analyzer/MAPS.c index c39143a..5f4118a 100644 --- a/MABS_run_analyzer/MAPS.c +++ b/MABS_run_analyzer/MAPS.c @@ -19,9 +19,6 @@ MAPS::MAPS() { MAPS::MAPS(Run* runp) { run = runp; initMapsRun(); - - UInt_t FileEvNb; // Event number per file - UInt_t TotEvNb; // Total events in RUN - should be less if RUN was stopped before endl // check if sepcified system is correct, if not, switch and reinitialize if ( switchsystem() ) @@ -160,6 +157,26 @@ MAPS::~MAPS(void) { }; + +void MAPS::save() { + if(fSave && fOk) + { + fOutputFile->cd(); + + fHitTree->Write("",TObject::kOverwrite); + fNoiseTree->Write("",TObject::kOverwrite); + fDynNoiseTree->Write("",TObject::kOverwrite); + + hint1->Write("",TObject::kOverwrite); + hint2->Write("",TObject::kOverwrite); + + fOutputFile->Save(); + fOutputFile->Close(); + + cout<< fRootFile <<" saved!"<Branch("frame" , &fFrameNumber , "frame/i" , 32000); // fDynNoiseTree->Branch("noise" , &fNoiseMean , "noise/F" , 32000); // fDynNoiseTree->Branch("pedestal" , &fPedestalsMean , "pedestal/F" , 32000); @@ -1427,8 +1439,6 @@ void MAPS::hitana() { bordercluster = true; else bordercluster = false; - Int_t posmatrixrow; - Int_t posmatrixcol; for(Int_t row=0; row<5; row++) { for(Int_t column=0; column<5; column++) diff --git a/MABS_run_analyzer/MAPS.h b/MABS_run_analyzer/MAPS.h index 77070bf..0efd512 100644 --- a/MABS_run_analyzer/MAPS.h +++ b/MABS_run_analyzer/MAPS.h @@ -78,7 +78,7 @@ private: Int_t fColumns; /// Variable to be able to subdivide the Matrix according to fOrderCode TString fMatrix; - Int_t fEventsSum; + UInt_t fEventsSum; Int_t fFile; Int_t fFrameNumber; /**< enum value 1 */ @@ -301,7 +301,15 @@ public: * */ Bool_t initNewRootFile(); + + /** + * @brief Saves the TTrees into ROOT file + * + * Saves the ROOT file + * + */ + void save(); /** * @brief Reopens an allready existing ROOT file diff --git a/MABS_run_analyzer/Run.c b/MABS_run_analyzer/Run.c index 737fe95..532dc9b 100644 --- a/MABS_run_analyzer/Run.c +++ b/MABS_run_analyzer/Run.c @@ -150,12 +150,14 @@ Bool_t Run::debugDBreadout() cout << "| CCE_in_Perc_1DB: " << std::right << colorwhite << labbook.CCE_in_Perc_1DB << endlr; cout << "|______________________________ " << endlr; cout << endlr; + return 0; } Bool_t Run::useDynamicalNoise(Bool_t var) { cout<<" Dynamical Noise calc. : " << colorwhite << (var?"1":"0") << colorreset << " <-- only used if RAW files are analyzed, force analysis to make sure" << endl; dynamicalNoise = var; + return 0; } @@ -199,22 +201,21 @@ Bool_t Run::analyzeRun(Bool_t force) progress = (Int_t)(((i-start)*100)/(nframes-1)*10); if (progress!=progress_tmp) { print_progress( (((i-start)*100.)/(nframes-1)) ); progress_tmp=progress;} } - cout << processed->plus << " vs. " << processed->minus << " : " << (processed->plus*1.0)/processed->minus<< endl; +// cout << processed->plus << " vs. " << processed->minus << " : " << (processed->plus*1.0)/processed->minus<< endl; // print a dummy file to indicate, that a root file was created once fstream* fout = new fstream(storepathRAWLinux + "/rootfilecreated",std::ios::out); *fout << "" << endl; fout->close(); - - // TODO - + processed->save(); } else { cout << "Skipped analysis of run " << labbook.runnumber << ", I think the root file for this run allready exists." << endl; cout << colorwhite << "initOldRootFile():" << endlr; - error = processed->initOldRootFile(); } + error = processed->initOldRootFile(); + if (!error) { cout << colorwhite << "binNoise():" << endlr; @@ -284,7 +285,6 @@ Bool_t Run::rescaleHistograms() histogramCalibrated.avgNoise = histogram.avgNoise * gain; histogramCalibrated.avgNoisePlus = histogram.avgNoisePlus * gain; histogramCalibrated.avgNoiseMinus = histogram.avgNoiseMinus * gain; - return 1; } @@ -381,6 +381,7 @@ Bool_t Run::generateReadableRunCode() else runcode+= Form("_m%.0f-%.0f_", devided_lower_matr_min, devided_lower_matr_max); } + return 0; } string Run::removeNumbers(std::string x) @@ -561,15 +562,12 @@ Bool_t Run::binNoise() Float_t noise; TBranch* noiseBranch; Double_t const probabilities[] = {0.3415/2, 0.5, 1-0.3415/2}; // sigma/2 from gaus to the left and to the right //{0.17, 0.5, 1-0.17}; - histogram.Noise->Reset(); - processed->fNoiseTree->SetBranchAddress("noise", &noise, &noiseBranch); - + processed->fNoiseTree->SetBranchAddress("noise", &noise, &noiseBranch); for (Int_t cnt=0; cntfNoiseTree->GetEntries(); cnt++) { processed->fNoiseTree->GetEntry(cnt); histogram.Noise->Fill(noise); - } - + } histogram.Noise->GetQuantiles( 3, noisequantiles, probabilities); histogram.avgNoise = noisequantiles[1]; histogram.avgNoisePlus = noisequantiles[2] - noisequantiles[1]; @@ -577,20 +575,15 @@ Bool_t Run::binNoise() // if (labbook.system == "PXI") // for (int j=0; j<3; j++) // noisequantiles[j] /= 16.0; // TODO analyze PXI scales - - return false; + return 0; } Bool_t Run::binSeedSumVeto() { - /// framenumber - UInt_t framen; - /// number of hits in current frame - UInt_t hitsInframe; - /// pixel number of seed pixel, position on sensor - UInt_t seedPixel[10000]; - /// Array of CLUSTERSIZE * CLUSTERSIZE clusters, seed pixel in the middle - Float_t pixelcluster[processed->clustersize*processed->clustersize][10000]; +// /// pixel number of seed pixel, position on sensor +// UInt_t seedPixel[10000]; +// /// Array of CLUSTERSIZE * CLUSTERSIZE clusters, seed pixel in the middle +// Float_t pixelcluster[processed->clustersize*processed->clustersize][10000]; /// collected charge in cluster Float_t pixelSum = 0; @@ -600,9 +593,9 @@ Bool_t Run::binSeedSumVeto() { processed->fHitTree->GetEntry(framei); // account only frames with less then 10 hits - if (processed->fFrameInfo.hits<10) + if (processed->fFrameInfo.hits<(unsigned int)10) { - for(Int_t hiti=0; hitifFrameInfo.hits;hiti++) + for(Int_t hiti=0; (unsigned int)hitifFrameInfo.hits;hiti++) { // histogram with the single pixel histogram.Seed->Fill(processed->fFrameInfo.p[12][hiti]); @@ -630,6 +623,7 @@ Bool_t Run::binSeedSumVeto() histogram.posSeed=FitPerform(histogram.Seed); histogram.posSum=FitPerform(histogram.Sum, "gaus"); gROOT->SetBatch(kFALSE); + return 0; } void Run::setPlotStyle(Int_t x){ @@ -648,7 +642,7 @@ Bool_t Run::plotNoise() noisecanvas = plot1DHistogram(histogram.Noise, "landau", "", legendEntry); return 0; } - return 0; + return 1; } Bool_t Run::plotSeed() @@ -668,7 +662,7 @@ Bool_t Run::plotSum() plot1DHistogram(histogram.Sum, "gaus"); return 0; } - return 0; + return 1; } Bool_t Run::plotVeto() @@ -678,7 +672,7 @@ Bool_t Run::plotVeto() plot1DHistogram(histogram.Veto, "gaus"); return 0; } - return 0; + return 1; } // Int_t findMaxBin(TH1F* th1fpointer) { @@ -739,6 +733,7 @@ Bool_t Run::plotAllHistograms(histogramstruct* histogramstructpointer) return 0; } + return 1; } Bool_t Run::plotAllHistograms() @@ -748,6 +743,7 @@ Bool_t Run::plotAllHistograms() plotAllHistograms(&histogram); return 0; } + return 1; } Bool_t Run::plotAllHistogramsCalibrated() @@ -765,6 +761,7 @@ Bool_t Run::plotAllHistogramsCalibrated() Float_t Run::FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verbose) { Float_t posMax = 0; + Float_t posMax2 = 0; Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax(); Float_t noiseborder = posMaxValHist/15; @@ -780,55 +777,74 @@ Float_t Run::FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verb { Float_t peak1 = histogrampointer->GetMaximumBin(); Float_t peak2; - if (histogrampointer->GetBinCenter(peak1)<200) + if (histogrampointer->GetBinCenter(peak1)<3.3*noiseborder) // vermutlich ist veto peak nicht am höchsten { Float_t avg = 0; for(Int_t bin=histogrampointer->GetXaxis()->FindBin(noiseborder);binFindLastBinAbove(0);bin++) avg += histogrampointer->GetBinContent(bin); avg /= histogrampointer->FindLastBinAbove(0) - histogrampointer->GetXaxis()->FindBin(noiseborder); - peak2 = histogrampointer->FindLastBinAbove(avg/4); + peak2 = histogrampointer->FindLastBinAbove(avg/3); } else { peak2 = peak1; peak1 = histogrampointer->GetXaxis()->FindBin(noiseborder); - } + } histogrampointer->GetXaxis()->SetRange(peak1,peak2-1); // look only for maxima with x greater than noiseborder, cut away noise Float_t min = histogrampointer->GetMinimumBin(); histogrampointer->GetXaxis()->SetRange(min,histogrampointer->FindLastBinAbove(0)); // look only for maxima with x greater than noiseborder, cut away noise // DEBUG - if (verbose) - { - cout << coloryellow << "peak1 " << histogrampointer->GetXaxis()->GetBinCenter(peak1) << endlr; - cout << coloryellow << "peak1val " << histogrampointer->GetBinContent(peak1) << endlr; - cout << coloryellow << "peak2 " << histogrampointer->GetXaxis()->GetBinCenter(peak2) << endlr; - cout << coloryellow << "peak2val " << histogrampointer->GetBinContent(peak2) << endlr; - cout << coloryellow << "min " << histogrampointer->GetXaxis()->GetBinCenter( min ) << endlr; - } +// if (verbose) +// { +// cout << coloryellow << "peak1 " << histogrampointer->GetXaxis()->GetBinCenter(peak1) << endlr; +// cout << coloryellow << "peak1val " << histogrampointer->GetBinContent(peak1) << endlr; +// cout << coloryellow << "peak2 " << histogrampointer->GetXaxis()->GetBinCenter(peak2) << endlr; +// cout << coloryellow << "peak2val " << histogrampointer->GetBinContent(peak2) << endlr; +// cout << coloryellow << "min " << histogrampointer->GetXaxis()->GetBinCenter( min ) << endlr; +// } histogrampointer->Fit(fitFunc, "N,Q,W", "", histogrampointer->GetXaxis()->GetBinCenter(min), posMaxValHist); + posMax = histogrampointer->GetXaxis()->GetBinCenter(histogrampointer->GetMaximumBin()); // Methode 1 histogrampointer->GetXaxis()->SetRange(histogrampointer->GetXaxis()->FindBin(noiseborder),histogrampointer->GetXaxis()->FindBin(posMaxValHist)); // look only for maxima with x greater than noiseborder, cut away noise } else { histogrampointer->Fit(fitFunc, "N,Q,W", "", noiseborder, posMaxValHist); + posMax = fitFunc->GetMaximumX(); // Methode 2 } Float_t sigma = fitFunc->GetParameter(2); - if (sigma > 260) + posMax2 = fitFunc->GetMaximumX(); // Methode 2 + Float_t peakposdifperc = abs(posMax-posMax2)/min(posMax,posMax2); + if (sigma > 260 || peakposdifperc>0.2) { if (verbose) { - cout << "Sigma suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma << endl; - cout << colorred << " Could not find " << histogrampointer->GetName() << " peak" << endlr; + if (sigma > 260) + { + cout << "Sigma suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma << endl; + cout << colorred << " Could not find " << histogrampointer->GetName() << " peak" << endlr; + } + if (peakposdifperc>0.2) + { + cout << "Maximum peak position and fit gaussian peak position doesn't match in " << histogrampointer->GetName() << " spectrum. Difference: " << peakposdifperc*100 <<" % " << endl; + cout << colorred << " Could not find " << histogrampointer->GetName() << " peak" << endlr; + } } return posMax; } - else if (sigma > 40) + else if (sigma > 40 || peakposdifperc>0.1) { if (verbose) { - cout << "Sigma suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma << endl; - cout << coloryellow << " Please check " << histogrampointer->GetName() << " peak position: " << colorreset << colorwhite << fitFunc->GetMaximumX() << endlr; + if (sigma > 260) + { + cout << "Sigma or suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma << endl; + cout << coloryellow << " Please check " << histogrampointer->GetName() << " peak position: " << colorreset << colorwhite << posMax << endlr; + } + if (peakposdifperc>0.1) + { + cout << "Maximum peak position and fit gaussian peak position in " << histogrampointer->GetName() << " have difference of " << peakposdifperc*100 <<" % " << endl; + cout << coloryellow << " Please check " << histogrampointer->GetName() << " peak position: " << colorreset << colorwhite << posMax << endlr; + } } - } - posMax = fitFunc->GetMaximumX(); + } fitFunc->DrawCopy("same"); TString legendEntry = TString(Form("%s",runcode.Data())); TLegend* leg = new TLegend(0.5,0.5,0.89,0.89);//(0.6,0.7,0.89,0.89); @@ -942,7 +958,7 @@ Bool_t Run::writeHistogramToFile(TH1F* onehistogram) TString filename= savepathresults + "/" + runcode + " " + onehistogram->GetName() + " hist.dat"; fstream* fout = new fstream(filename,ios::out); - TString header = Form("#bin [ADU]\tCounts\tbin [e]\tbin [Veto%]\tSeed"); + TString header = Form("#bin [ADU]\tCounts"); *fout << header << endl; TString outline; @@ -963,7 +979,7 @@ Bool_t Run::writeAllHistogramsToFile() TString filename= savepathresults + "/" + runcode + " histograms.dat"; fstream* fout = new fstream(filename,ios::out); - TString header = Form("#bin [ADU]\tbin [e]\tSeed\tSum\t\Veto\tbin noise [ADU]\tbin noise [e]\tnoise\tbin [Veto%]\tSeed\n"); + TString header = Form("#bin [ADU]\tbin [e]\tSeed\tSum\t\Veto\tbin noise [ADU]\tbin noise [e]\tnoise\n"); header += Form("#posVeto, run: %.1f, DB: %.1f, Fe55 DB (%d, %.1f): %.1f\n", histogram.posVeto, labbook.posVetoDB, Fe55run.posVetorunnumber, Fe55run.temperature, Fe55run.posVeto); header += Form("#posSeed, run: %.1f, DB: %.1f\n", histogram.posSeed, labbook.posSeedDB); header += Form("#posSum, run: %.1f, DB: %.1f\n", histogram.posSum, labbook.posSumDB);