From: Benjamin Linnik Date: Sun, 6 Dec 2015 09:36:54 +0000 (+0100) Subject: Run analyzer: added improved summary plot, expanded Seed Percentage histogram scala... X-Git-Url: https://jspc29.x-matter.uni-frankfurt.de/git/?a=commitdiff_plain;h=1d4008fedc1b1b97d83cdc25e71a5c3af028506d;p=radhard.git Run analyzer: added improved summary plot, expanded Seed Percentage histogram scala, added ChargeSpectrum.c due to important changes in it, changed overall style of plots --- diff --git a/.gitignore b/.gitignore index aa483a7..cdc4701 100644 --- a/.gitignore +++ b/.gitignore @@ -10,6 +10,6 @@ *.MOV */results/* *runlist.txt -ChargeSpektrum.c +#ChargeSpektrum.c MABS_run_analyzer/html/* MABS_run_analyzer/latex/* \ No newline at end of file diff --git a/MABS_run_analyzer/ChargeSpektrum.c b/MABS_run_analyzer/ChargeSpektrum.c new file mode 100644 index 0000000..a52a2ae --- /dev/null +++ b/MABS_run_analyzer/ChargeSpektrum.c @@ -0,0 +1,180 @@ +/** + * @file ChargeSpektrum.c + * @brief Use brief, otherwise the index won't have a brief explanation. + * + * Detailed explanation. + * + * + */ + +#include "MAPS.c" +#include "Run.c" +#include "CSVRow.h" +#include "CSVRow.C" +#include +#include "help.h" + +Run** runs; +Int_t numberRuns; +Bool_t isBatch = kFALSE; +TString savepathresults = "./results/"; + +#include "ChargeSpektrumFunctions.c" + +void ChargeSpektrum(TString runnumber = "") +{ + #include "SetStyle.c" + cout << endl << endl; + if (gROOT->IsBatch()) + isBatch = kTRUE; + + numberRuns=0; + std::vector runList; + std::vector sumuprunList; + if (runnumber.Length() > 0) + { + if (runnumber.Contains("-")) { + TObjArray* runarray = runnumber.Tokenize("-"); + if (runarray->GetEntries()==2) + { + TObjString* tempstrobj; + tempstrobj=static_cast(runarray->At(0)); + if (tempstrobj->GetString().Length()>0) + { + Int_t tempintstart = tempstrobj->GetString().Atoi(); + tempstrobj=static_cast(runarray->At(1)); + if (tempstrobj->GetString().Length()>0) + { + Int_t tempintend = tempstrobj->GetString().Atoi(); + if (tempintend-tempintstart > 0) + { + for (Int_t i=tempintstart; i <= tempintend; i++) + { + runList.push_back(i); + numberRuns++; + } + } + else + cout << coloryellow << "Invalid run number range "<< colorreset << colorwhite << tempintstart << " till " << tempintend << endlr; + } + else + cout << coloryellow << "Invalid run number "<< colorreset << colorwhite << tempstrobj->GetString() << endlr; + + } + else + cout << coloryellow << "Invalid run number "<< colorreset << colorwhite << tempstrobj->GetString() << endlr; + } + else + { + cout << endl << colorred << "Given parameters have wrong format, please enter integer run numbers devided by '" << colorreset << colorwhite << "," << colorreset << colorred << " or 2 runnumbers devided by " << colorreset << colorwhite << "-" << colorreset << colorred << "' and surrounded by '" << colorreset << colorwhite << "\"" << colorreset << colorred << "'." << endlr; + cout << "For example run: " << colorwhite << "root -l 'ChargeSpektrum.c+(\"342815,342816\")'" << endl; + exit(1); + } + + } else { + TObjArray* runarray = runnumber.Tokenize(","); + if (runarray->GetEntries()>0) + { + TObjString* tempstrobj; + for (Int_t i=0; i < runarray->GetEntries(); i++) + { + tempstrobj=static_cast(runarray->At(i)); + if (tempstrobj->GetString().Length()>0) + { + TString sumrunsstr = ""; + TObjArray* sumrunarray = tempstrobj->GetString().Tokenize("+"); + TObjString* sumuprunsstrobj; + for (Int_t j=0; j < sumrunarray->GetEntries(); j++) + { + sumuprunsstrobj=static_cast(sumrunarray->At(j)); + Int_t tempint = sumuprunsstrobj->GetString().Atoi(); + if (tempint > 0) + { + runList.push_back(tempint); + sumrunsstr.Append(Form("%d,",tempint)); + numberRuns++; + } + else + cout << coloryellow << "Invalid run number "<< colorreset << colorwhite << tempint << endlr; + } + sumuprunList.push_back(sumrunsstr); + } + else + cout << coloryellow << "Invalid run number "<< colorreset << colorwhite << tempstrobj->GetString() << endlr; + } + } + else + { + cout << endl << colorred << "Given parameters have wrong format, please enter integer run numbers devided by '" << colorreset << colorwhite << "," << colorreset << colorred << "' and surrounded by '" << colorreset << colorwhite << "\"" << colorreset << colorred << "'." << endlr; + cout << "For example run: " << colorwhite << "root -l 'ChargeSpektrum.c+(\"342815,342816\")'" << endl; + exit(1); + } + } + } + else + { + /// number of runs to be analyzed, number of lines read by @c ReadRunList() + ReadRunList(&runList); + } + runs = new Run*[numberRuns]; + + cout << "Found " << numberRuns << " run(s) in 'runlist.txt' or given as parameters." << endl; + for(Int_t runi=0;runi0) + { + runs[runi] = new Run(runList[runi], runi); + if (!runs[runi]->error) + { + runs[runi]->setResultsPath(savepathresults); + runs[runi]->error=runs[runi]->initRun(); + runs[runi]->setDynamicalNoiseMode("simple"); + runs[runi]->useDynamicalNoise(true); + runs[runi]->useCommonModeFilter(false); + runs[runi]->setFixedThresholdValueElectrons(1440); + +// runs[runi]->analyzeFrame(57826); +// runs[runi]->analyzeFrame(57827); +// runs[runi]->analyzeFrame(983603); + + // creates or opens .root file, can analyze the RAW data + runs[runi]->error=runs[runi]->analyzeRun(false); // creates or opens .root file, can analyze the RAW data + if (runs[runi]->error) // first time analysis without forcing new ROOT file failed, retry with forcing new analysis + { + runs[runi]->error = false; + runs[runi]->error=runs[runi]->analyzeRun(true); + } + if (!runs[runi]->error) + { + // gROOT->SetBatch(kTRUE); + if (!isBatch) + gROOT->SetBatch(kFALSE); +// runs[runi]->plot1DHistogram(runs[runi]->histogram,runs[runi]->histogram->Seed,"landau","Test"); +// runs[runi]->compareHistogramClassVector.push_back(runs[runi]->histogram->calibrated); +// runs[runi]->compareHistogramClassVector.push_back(runs[runi]->histogramthreshold->calibrated); +// runs[runi]->compareHistogramClassVector.push_back(runs[runi]->histogramfixedthreshold); +// runs[runi]->plotCompareHistograms(); + runs[runi]->plotAllHistograms(runs[runi]->histogramthreshold); + compareHistogramClassVector.push_back(runs[runi]->histogramthreshold->calibrated); +// compareHistogramClassVector2.push_back(runs[runi]->histogram->calibrated); + // runs[runi]->plot1DHistogram(runs[runi]->histogramthreshold, runs[runi]->histogramthreshold->Veto, "gaus"); +// runs[runi]->plot1DHistogram(runs[runi]->histogramthreshold, runs[runi]->histogramthreshold->SeedPerc, "landau", 0); + compareHistogramVector.push_back(runs[runi]->histogramthreshold->calibrated->Seed); + runs[runi]->writeAllHistogramsToFile(); + } + //cout << runs[runi]->histogram + } + } + } + plotAllRuns(&compareHistogramClassVector); + printSummaryTable(&compareHistogramClassVector); +// plotAllRuns(&compareHistogramClassVector2); +// plotAllRuns(&compareHistogramClassVector3); + writeObservableToFile(); +// plotAllRuns("seed threshold calibrated"); +// setCustomPath("Excel/"); +// writeObservableToFile("seed threshold calibrated"); +// writeObservableToFile("seed threshold"); +// writeObservableToFile("seed threshold"); +// writeObservableToFile("sum threshold"); +} diff --git a/MABS_run_analyzer/ChargeSpektrumFunctions.c b/MABS_run_analyzer/ChargeSpektrumFunctions.c index 752677e..e2872f6 100644 --- a/MABS_run_analyzer/ChargeSpektrumFunctions.c +++ b/MABS_run_analyzer/ChargeSpektrumFunctions.c @@ -8,8 +8,11 @@ #include #include #include +#include +#include using namespace std; +TString trimRunnumberAtBegin(TString str); Int_t* ReadRunList(); Int_t ReadRunList(std::vector*); Bool_t plotAllRuns(); @@ -167,19 +170,37 @@ Bool_t plotAllRuns(vector* ptCompareHistogramClassVector) if (testifMixingCalibration(ptCompareHistogramClassVector)) return 1; // legend entries - Float_t height = ptCompareHistogramClassVector->size() * 0.04; - TLegend* leg1 = new TLegend(0.4,0.89-height,0.89,0.89);//(0.6,0.7,0.89,0.89); - leg1->SetTextSize(0.02); + Float_t height = ptCompareHistogramClassVector->size() * 0.045; + TLegend* leg1 = new TLegend(0.35,0.89-height,0.89,0.89);//(0.6,0.7,0.89,0.89); + leg1->SetTextSize(0.025); + leg1->SetFillStyle(0); + leg1->SetTextFont(12); leg1->SetFillColor(0); leg1->SetBorderSize(0); TLegend* leg2 = (TLegend*) leg1->Clone(); TString legendEntry; + TPaveText *owntitle = new TPaveText(0.15,0.9,0.930401,0.995,"nbNDC"); + owntitle->SetFillStyle(0); + owntitle->SetBorderSize(0); + owntitle->SetTextAlign(22); + owntitle->SetTextSize(0.05); + owntitle->SetTextFont(42); + TPaveText *owntitle2 = (TPaveText*)owntitle->Clone(); + TPaveText *owntitle3 = (TPaveText*)owntitle->Clone(); + TPaveText *owntitle4 = (TPaveText*)owntitle->Clone(); + gStyle->SetOptTitle(0); + Float_t lastbin1=0; Float_t lastbin2=0; Float_t lastbin3=0; - Float_t lastbin4=0; - TString canvastitle = Form("Summary"); - TTimeStamp* time = new TTimeStamp(); + Float_t lastbin4=0; + TString canvastitle = Form("Summary"); + for (UInt_t histogrami=0; histogrami < ptCompareHistogramClassVector->size(); histogrami++) + { + HistogramType* curhistogramclassp = ptCompareHistogramClassVector->at(histogrami); + canvastitle+= Form("_%d",curhistogramclassp->labbook->runnumber); + } + TTimeStamp* time = new TTimeStamp(); TString canvasname = Form("%d",time->GetNanoSec()); TCanvas* canvas = new TCanvas(canvasname, canvastitle, 1200, 800); canvas->Divide(2,2); @@ -187,9 +208,15 @@ Bool_t plotAllRuns(vector* ptCompareHistogramClassVector) { HistogramType* curhistogramclassp = ptCompareHistogramClassVector->at(histogrami); TH1F* curhistogramclone; + canvas->cd(1); curhistogramclone = (TH1F*) curhistogramclassp->Seed->Clone(); +// TLatex Tl; +// Tl.SetNDC(); +// Tl.DrawText(0.3, 0.2, "Titel titel"); + //TPaveText *t = new TPaveText(0.0, 0.9, 0.3, 1.0, "brNDC"); // left-up // curhistogramclone->SetLineColor(rootcolors[histogrami]); + curhistogramclone->Draw("SAME"); legendEntry = Form("%d %s", curhistogramclassp->labbook->runnumber, curhistogramclone->GetTitle()); leg1->AddEntry(curhistogramclone, legendEntry, "l"); @@ -197,6 +224,10 @@ Bool_t plotAllRuns(vector* ptCompareHistogramClassVector) lastbin1 = (curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,1))>lastbin1)?curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,1)):lastbin1; curhistogramclone->SetAxisRange(0,lastbin1*1.1,"X"); gPad->SetLogy(1); + owntitle->Clear(); + owntitle->AddText(trimRunnumberAtBegin(curhistogramclone->GetName())); + owntitle->Draw("SAME"); + canvas->cd(2); curhistogramclone = (TH1F*) curhistogramclassp->Sum->Clone(); // curhistogramclone->SetLineColor(rootcolors[histogrami]); @@ -205,6 +236,10 @@ Bool_t plotAllRuns(vector* ptCompareHistogramClassVector) lastbin2 = (curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,1))>lastbin2)?curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,1)):lastbin2; curhistogramclone->SetAxisRange(0,lastbin2*1.1,"X"); gPad->SetLogy(1); + owntitle2->Clear(); + owntitle2->AddText(trimRunnumberAtBegin(curhistogramclone->GetName())); + owntitle2->Draw("SAME"); + canvas->cd(3); curhistogramclone = (TH1F*) curhistogramclassp->Veto->Clone(); // curhistogramclone->SetLineColor(rootcolors[histogrami]); @@ -212,6 +247,10 @@ Bool_t plotAllRuns(vector* ptCompareHistogramClassVector) leg1->Draw("SAME"); lastbin3 = (curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,1))>lastbin3)?curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,1)):lastbin3; curhistogramclone->SetAxisRange(0,lastbin3*1.1,"X"); + owntitle3->Clear(); + owntitle3->AddText(trimRunnumberAtBegin(curhistogramclone->GetName())); + owntitle3->Draw("SAME"); + canvas->cd(4); curhistogramclone = (TH1F*) curhistogramclassp->Noise->Clone(); // curhistogramclone->SetLineColor(rootcolors[histogrami]); @@ -219,8 +258,24 @@ Bool_t plotAllRuns(vector* ptCompareHistogramClassVector) legendEntry = Form("%d Noise: %.2f + %.2f - %.2f",curhistogramclassp->labbook->runnumber, curhistogramclassp->avgNoise, curhistogramclassp->avgNoisePlus, curhistogramclassp->avgNoiseMinus); leg2->AddEntry(curhistogramclone, legendEntry, "l"); leg2->Draw(); + owntitle4->Clear(); + owntitle4->AddText(trimRunnumberAtBegin(curhistogramclone->GetName())); + owntitle4->Draw("SAME"); } + canvas->Update(); + //gStyle->SetOptTitle(1); + + TImageDump *img = new TImageDump(savepathresults + "/" + canvastitle + ".png"); + canvas->Paint(); + img->Close(); + + TFile *f = new TFile(savepathresults + "/" + canvastitle + ".root","RECREATE"); + f->cd(); + f->Append(canvas); + f->Append(img); + f->Write(); + return 0; } return 1; @@ -236,6 +291,7 @@ Bool_t printSummaryTable(vector* ptCompareHistogramClassVector) cout << left << setw(width*1) << setfill(' ') << "Cut type: "<< left << setw(width*8) << setfill(' ') << ptCompareHistogramClassVector->at(0)->histogramdescription << endl; cout << left << setw(width-5) << setfill(' ') << "Runnumber"; + cout << left << setw(width-10) << setfill(' ') << "Chp#"; cout << left << setw(width-10) << setfill(' ') << "Mtrx"; Float_t printradion = 0; for (vector::iterator curHistogramClass = ptCompareHistogramClassVector->begin(); curHistogramClass != ptCompareHistogramClassVector->end(); curHistogramClass++) @@ -260,12 +316,16 @@ Bool_t printSummaryTable(vector* ptCompareHistogramClassVector) printIntegral += (*curHistogramClass)->sr90IntegralVal; if (printIntegral != 0) cout << left << setw(width) << setfill(' ') << "Integral"; cout << left << setw(width) << setfill(' ') << "Seed Peak"; + cout << left << setw(width) << setfill(' ') << "Veto Peak"; cout << left << setw(width) << setfill(' ') << "Noise"; + cout << left << setw(width) << setfill(' ') << "Noise threshold"; + //cout << left << setw(width) << setfill(' ') << "Frames"; cout << endl; for (UInt_t histogrami=0; histogrami < ptCompareHistogramClassVector->size(); histogrami++) { HistogramType* curhistogramclassp = ptCompareHistogramClassVector->at(histogrami); cout << left << setw(width-5) << setfill(' ') << curhistogramclassp->labbook->runnumber; + cout << left << setw(width-10) << setfill(' ') << curhistogramclassp->labbook->chip; cout << left << setw(width-10) << setfill(' ') << curhistogramclassp->labbook->matrix; if (printradion != 0) cout << left << setw(width-5) << setfill(' ') << curhistogramclassp->labbook->radDoseIon; if (printradnonion != 0) cout << left << setw(width-5) << setfill(' ') << curhistogramclassp->labbook->radDoseNonIon; @@ -274,7 +334,10 @@ Bool_t printSummaryTable(vector* ptCompareHistogramClassVector) if (printStoN != 0) cout << left << setw(width) << setfill(' ') << printTableElement(curhistogramclassp->StoN,ptCompareHistogramClassVector->at(0)->StoN); if (printIntegral != 0) cout << left << setw(width) << setfill(' ') << printTableElement(curhistogramclassp->sr90IntegralVal,ptCompareHistogramClassVector->at(0)->sr90IntegralVal); cout << left << setw(width) << setfill(' ') << printTableElement(curhistogramclassp->posSeed,ptCompareHistogramClassVector->at(0)->posSeed); + cout << left << setw(width) << setfill(' ') << printTableElement(curhistogramclassp->posVeto,ptCompareHistogramClassVector->at(0)->posVeto); cout << left << setw(width) << setfill(' ') << printTableElement(curhistogramclassp->avgNoise,ptCompareHistogramClassVector->at(0)->avgNoise); + cout << left << setw(width) << setfill(' ') << printTableElement(curhistogramclassp->noisethresholdborder,ptCompareHistogramClassVector->at(0)->noisethresholdborder); + //cout << left << setw(width) << setfill(' ') << printTableElement(curhistogramclassp->frames_found,ptCompareHistogramClassVector->at(0)->frames_found); cout << "" << endl; } cout << right << setw(width*9) << setfill('=') << " " << endl; @@ -292,7 +355,7 @@ template string printTableElement(varType t1, varType t2, cons string to_str_w_prec(const Float_t a_value, const int precision) { std::ostringstream out; - if (abs(a_value) > 1) { + if (abs(a_value) > 1 && abs(a_value) < 10000) { out << std::fixed << std::setprecision(precision) << a_value; } else { out << std::scientific << std::setprecision(precision) << a_value << std::fixed; out.setf(ios_base::fixed); @@ -300,3 +363,9 @@ string to_str_w_prec(const Float_t a_value, const int precision) return out.str(); } +TString trimRunnumberAtBegin(TString str) { + size_t i = 0; + while(isdigit(str(i++))); + return str.Remove(0, i); +} + diff --git a/MABS_run_analyzer/HistogramType.c b/MABS_run_analyzer/HistogramType.c index c5d8520..78a8f70 100644 --- a/MABS_run_analyzer/HistogramType.c +++ b/MABS_run_analyzer/HistogramType.c @@ -8,6 +8,7 @@ #define __HISTOGRAMTYPE__C #include "HistogramType.h" +#include using namespace std; @@ -44,7 +45,12 @@ void HistogramType::initHistograms(Int_t gotcolor, Int_t gotstyle) { initHistogram(Sum, "Sum" + histogramdescription, color, style); initHistogram(Veto, "Veto" + histogramdescription, color, style); initHistogram(Noise, "Noise" + histogramdescription, color, style); + initHistogram(SeedPerc, "Seed Percentage" + histogramdescription, color, style); + SeedPerc->SetBins(120, 0, 120); + SeedPerc->GetXaxis()->SetRangeUser(0,50); + SeedPerc->GetYaxis()->SetTitle(Form("Entries [1/%%]")); + SeedPerc->GetXaxis()->SetTitle("Q_coll [%]"); Noise->SetBins(cursystempar->nbinsnoise, 0, cursystempar->maxbinnoise); } @@ -88,11 +94,14 @@ Bool_t HistogramType::calibrateHistograms( Float_t gotgain ) { if (Sum != 0) calibrateHistogram(calibrated->Sum, Sum); if (Veto != 0) calibrateHistogram(calibrated->Veto, Veto); if (Noise != 0) calibrateHistogram(calibrated->Noise, Noise); + if (SeedPerc != 0) calibrated->SeedPerc = SeedPerc; if (histAvgCluster != 0) calibrate2DHistogramCounts(calibrated->histAvgCluster, histAvgCluster); calibrated->posSeed = posSeed * gain; calibrated->posSum = posSum * gain; calibrated->posVeto = posVeto * gain; + calibrated->posSeedPerc = posSeedPerc; + calibrated->sigmaSeedPerc = sigmaSeedPerc; calibrated->avgNoise = avgNoise * gain; calibrated->avgNoisePlus = avgNoisePlus * gain; calibrated->avgNoiseMinus = avgNoiseMinus * gain; @@ -108,8 +117,8 @@ Bool_t HistogramType::calibrateHistograms( Float_t gotgain ) { void HistogramType::calibrateHistogram(TH1F* &histogrampointernew, TH1F* &histogrampointerold) { histogrampointernew = (TH1F*)histogrampointerold->Clone(); - histogrampointernew->SetName(Form("%sC", histogrampointerold->GetName())); - histogrampointernew->SetTitle(Form("%sC", histogrampointerold->GetTitle())); + histogrampointernew->SetName(Form("%s Calibr.", histogrampointerold->GetName())); + histogrampointernew->SetTitle(Form("%s Calibr.", histogrampointerold->GetTitle())); histogrampointernew->GetXaxis()->SetTitle("Q_coll [e]"); int nbins = histogrampointernew->GetXaxis()->GetNbins(); double new_bins[nbins+1]; @@ -123,8 +132,8 @@ void HistogramType::calibrateHistogram(TH1F* &histogrampointernew, TH1F* &histog void HistogramType::calibrate2DHistogramCounts(TH2F* &histogrampointernew, TH2F* &histogrampointerold ) { histogrampointernew = (TH2F*)histogrampointerold->Clone(); - histogrampointernew->SetName(Form("%sC", histogrampointerold->GetName())); - histogrampointernew->SetTitle(Form("%sC", histogrampointerold->GetTitle())); + histogrampointernew->SetName(Form("%s Calibr.", histogrampointerold->GetName())); + histogrampointernew->SetTitle(Form("%s Calibr.", histogrampointerold->GetTitle())); histogrampointernew->GetZaxis()->SetTitle("Q_coll [e]"); int nbinsx = histogrampointernew->GetXaxis()->GetNbins(); int nbinsy = histogrampointernew->GetYaxis()->GetNbins(); @@ -136,16 +145,19 @@ void HistogramType::calibrate2DHistogramCounts(TH2F* &histogrampointernew, TH2F* } } -Float_t HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verbose, Double_t* parameters) { +Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verbose, Float_t fitstart) { + Double_t* parameters = 0; Float_t posMax = 0; Float_t posMax2 = 0; Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax(); // posMaxValHist/10 for USB system, the value is 90 Float_t noiseborder = 90; - if (noisethresholdborder>=0) + if (fitstart >= 0) + noiseborder = fitstart; + else if (noisethresholdborder>=0 ) noiseborder = noisethresholdborder; - if (fitFuncType.Contains("gaus")) + if (fitFuncType.EqualTo("gaus")) { histogrampointer->GetXaxis()->SetRange(histogrampointer->GetXaxis()->FindBin(noiseborder),histogrampointer->GetXaxis()->FindBin(posMaxValHist)); // look only for maxima with x greater than noiseborder, cut away noise // Int_t xValMax = histogrampointer->GetBinCenter(histogrampointer->GetMaximumBin()); @@ -213,6 +225,11 @@ Float_t HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, B histogrampointer->Fit(fitFunc, "N,Q,W", "", noiseborder, posMaxValHist); posMax = fitFunc->GetMaximumX(); // Methode 2 } + parameters = (Double_t *)calloc(3, sizeof(Double_t)); + for (Int_t pari=0; pari<3; pari++) { + //cout << colorcyan << fitFunc->GetParameter(pari) << endlr; + parameters[pari]=fitFunc->GetParameter(pari); + } Float_t sigma = fitFunc->GetParameter(2); posMax2 = fitFunc->GetMaximumX(); // Methode 2 // if (verbose) @@ -233,7 +250,8 @@ Float_t HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, B cout << colorred << " Could not find " << histogrampointer->GetName() << " peak" << endlr; } } - return 0; + parameters[1] = 0; + return parameters; } else if (sigma > 40 || peakposdifperc>0.1) { @@ -273,6 +291,10 @@ Float_t HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, B Float_t minFitMax = 1000; Float_t maxFitMax = 1000; histogrampointer->Fit(fitFunc, "N,Q,W", "", noiseborder, posMaxValHist); + parameters = (Double_t *)calloc(3, sizeof(Double_t)); + for (Int_t pari=0; pari<3; pari++) + parameters[pari]=fitFunc->GetParameter(pari); + posMax = fitFunc->GetParameter(1); // new method, just use MPV of fit, as noise border is well known. Compare with old method anyway and warn if something suspecious. // cout << "fitFunc->GetParameter(0): " << fitFunc->GetParameter(0) << endl; // cout << "fitFunc->GetParameter(1): " << fitFunc->GetParameter(1) << endl; @@ -298,7 +320,7 @@ Float_t HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, B posMax2 = fitMax1 + fitMax2 + fitMax3 - minFitMax - maxFitMax; if (abs(posMax-posMax2)/posMax > 0.11) { - cout << coloryellow << " Please check the Landau MPV value fit, as two alternate methods revealed results wich differ " << abs(posMax-posMax2)/posMax*100 << " %." << endlr; + cout << coloryellow << " Please check the Landau MPV value fit of " << histogrampointer->GetName() << ", as two alternate methods revealed results wich differ " << abs(posMax-posMax2)/posMax*100 << " %." << endlr; cout << " Primary, new MPV method: " << posMax << endl; cout << " Secondary, mean maxima mathod: " << posMax2 << endl; } @@ -329,7 +351,8 @@ Float_t HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, B // fitFcn->FixParameter(3,-22.43016284); //histogrampointer->Fit("fitFcn","V+","ep"); - histogrampointer->Fit("fitFcn","Q,W","ep"); + histogrampointer->Fit("fitFcn","Q,W","ep"); + parameters = (Double_t *)calloc(4, sizeof(Double_t)); for (Int_t pari=0; pari<4; pari++) parameters[pari]=fitFcn->GetParameter(pari); if (verbose) @@ -341,7 +364,7 @@ Float_t HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, B } } - return posMax; + return parameters; } Bool_t HistogramType::calculteCCE(Bool_t verbose) { diff --git a/MABS_run_analyzer/HistogramType.h b/MABS_run_analyzer/HistogramType.h index 700cc99..76c7a17 100644 --- a/MABS_run_analyzer/HistogramType.h +++ b/MABS_run_analyzer/HistogramType.h @@ -107,6 +107,8 @@ public: TH1F* Noise = 0; /// histogram with cluster data, rotated if staggered TH2F* histAvgCluster = 0; + /// Percentage of charge gathered by seed pixel (Seed to Sum ratio) + TH1F* SeedPerc = 0; //***************** // FITTING RESULTS OF TH1F HISTOGRAMS @@ -130,6 +132,11 @@ public: Double_t sr90IntegralVal = -1; Double_t sr90IntegralErr = -1; + /// Peak and sigma of Seed to Sum ratio + Float_t posSeedPerc = 0; + Float_t sigmaSeedPerc = 0; + + /// threshold for integrating the Sr90 spectrum, is set dynamically in @c integrateSr90Spectra() Float_t noisethresholdborder = -1; @@ -186,7 +193,7 @@ public: * @return peak position of the fit * */ - Float_t FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verbose = 0, Double_t* parameters = 0); + Double_t* FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verbose = 0, Float_t fitstart = -1); /** @brief calculates Charge Collection Efficiency if Fe55 run */ Bool_t calculteCCE(Bool_t verbose=false); /** diff --git a/MABS_run_analyzer/MAPS.c b/MABS_run_analyzer/MAPS.c index f9c1582..2cc1632 100644 --- a/MABS_run_analyzer/MAPS.c +++ b/MABS_run_analyzer/MAPS.c @@ -463,7 +463,7 @@ bool MAPS::checkDataFiles( Int_t MaxFiles ) { fOk = true; } cout<<"-----------------------"<posSum, 4); constructUpdateString(&sqlupdatequery, "SeedPeak", histogram->posSeed, 4); constructUpdateString(&sqlupdatequery, "VetoPeak", histogram->posVeto, 4); - constructUpdateString(&sqlupdatequery, "Avg.Noise", histogram->calibrated->avgNoise); - constructUpdateString(&sqlupdatequery, "Avg.Noise+", histogram->calibrated->avgNoisePlus, 2); - constructUpdateString(&sqlupdatequery, "Avg.Noise-", histogram->calibrated->avgNoiseMinus, 2); + if (histogram->calibrated != 0) + constructUpdateString(&sqlupdatequery, "Avg.Noise", histogram->calibrated->avgNoise); + if (histogram->calibrated != 0) + constructUpdateString(&sqlupdatequery, "Avg.Noise+", histogram->calibrated->avgNoisePlus, 2); + if (histogram->calibrated != 0) + constructUpdateString(&sqlupdatequery, "Avg.Noise-", histogram->calibrated->avgNoiseMinus, 2); constructUpdateString(&sqlupdatequery, "CCE_1", histogram->CCE_in_Perc_1); constructUpdateString(&sqlupdatequery, "CCE_25", histogram->CCE_in_Perc_25); constructUpdateString(&sqlupdatequery, "StoN", histogram->StoN, 3); constructUpdateString(&sqlupdatequery, "Avg.NoiseADC", histogram->avgNoise); constructUpdateString(&sqlupdatequery, "Frames_found", frames_found,100000000); - constructUpdateString(&sqlupdatequery, "Sr90IntegralVal", histogramfixedthreshold->calibrated->sr90IntegralVal,1000000000); + if (histogram->calibrated != 0) + constructUpdateString(&sqlupdatequery, "Sr90IntegralVal", histogramfixedthreshold->calibrated->sr90IntegralVal,1000000000); if (sqlupdatequery.length()>0) { @@ -801,7 +805,10 @@ Bool_t Run::binSeedSumVeto() notSeedSum += processed->fFrameInfo.p[clusteri][hiti]; } } - histogram->Sum->Fill(pixelSum); + histogram->Sum->Fill(pixelSum); + + // seed percentage spectrum + histogram->SeedPerc->Fill(processed->fFrameInfo.p[12][hiti]/pixelSum*100); // veto spectrum if (TMath::Abs(notSeedSum) < cursystemparam.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))) @@ -813,6 +820,7 @@ Bool_t Run::binSeedSumVeto() histogramthreshold->Seed->Fill(processed->fFrameInfo.p[12][hiti]); histogramthreshold->Sum->Fill(pixelSum); + histogramthreshold->SeedPerc->Fill(processed->fFrameInfo.p[12][hiti]/pixelSum*100); if (TMath::Abs(notSeedSum) < cursystemparam.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))) histogramthreshold->Veto->Fill(processed->fFrameInfo.p[12][hiti]); // histogram with the single pixel } @@ -823,6 +831,7 @@ Bool_t Run::binSeedSumVeto() histogramfixedthreshold->Seed->Fill(processed->fFrameInfo.p[12][hiti]); histogramfixedthreshold->Sum->Fill(pixelSum); + histogramfixedthreshold->SeedPerc->Fill(processed->fFrameInfo.p[12][hiti]/pixelSum*100); if (TMath::Abs(notSeedSum) < cursystemparam.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))) histogramfixedthreshold->Veto->Fill(processed->fFrameInfo.p[12][hiti]); // histogram with the single pixel } @@ -834,15 +843,27 @@ Bool_t Run::binSeedSumVeto() // gROOT->SetBatch(kTRUE); for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) { + Double_t* parameters = (Double_t *)calloc(3, sizeof(Double_t)); (*curHistogramClass)->FindNoisethresholdborder((*curHistogramClass)->Seed, false); if (labbook.source.Contains("Fe")||labbook.source.Contains("Cd")) - (*curHistogramClass)->posVeto=(*curHistogramClass)->FitPerform((*curHistogramClass)->Veto, "gaus", true); - (*curHistogramClass)->posSeed=(*curHistogramClass)->FitPerform((*curHistogramClass)->Seed, "landau", true); - (*curHistogramClass)->posSum=(*curHistogramClass)->FitPerform((*curHistogramClass)->Sum, "gaus", true); + parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Veto, "gaus", true); + (*curHistogramClass)->posVeto = parameters[1]; + parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Seed, "landau", true); + (*curHistogramClass)->posSeed = parameters[1]; + parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Sum, "gaus", true); + (*curHistogramClass)->posSum = parameters[1]; + + parameters =(*curHistogramClass)->FitPerform((*curHistogramClass)->SeedPerc, "landau", true); + (*curHistogramClass)->posSeedPerc = parameters[1]; + (*curHistogramClass)->sigmaSeedPerc = parameters[2]; + for (Int_t bini=1; bini <= cursystemparam.nbins; bini++) { // TODO: rescaled histogram to number of frames found, remove the inner for loop later // (*curHistogramClass)->Seed->SetBinContent(bini,(*curHistogramClass)->Seed->GetBinContent(bini)/(frames_found*1.0)); // (*curHistogramClass)->Sum->SetBinContent(bini,(*curHistogramClass)->Sum->GetBinContent(bini)/(frames_found*1.0)); } + for (Int_t bini=1; bini <= 100; bini++) { + (*curHistogramClass)->SeedPerc->SetBinContent(bini,(*curHistogramClass)->SeedPerc->GetBinContent(bini)/0.01/((*curHistogramClass)->SeedPerc->GetEntries())); + } } histogramfixedthreshold->noisethresholdborder=0; // gROOT->SetBatch(kFALSE); @@ -1325,13 +1346,13 @@ Bool_t Run::plotClusterDistribution(HistogramType* HistogramTypepointer) // create lorentz fit of a slice TH1F* xslicetroughcluster = new TH1F(Form("%s_xslice", HistogramTypepointer->histAvgCluster->GetTitle()),"X slice",HistogramTypepointer->histAvgCluster->GetXaxis()->GetNbins(),HistogramTypepointer->histAvgCluster->GetXaxis()->GetBinLowEdge(1),HistogramTypepointer->histAvgCluster->GetXaxis()->GetBinUpEdge(HistogramTypepointer->histAvgCluster->GetXaxis()->GetNbins())); - Double_t xslicetroughclusterpar[4]; + Double_t* xslicetroughclusterpar = (Double_t *)calloc(4, sizeof(Double_t)); // allocate 6 parameters for safety, maximum 4 used at the moment Int_t middlebin = (int)(HistogramTypepointer->histAvgCluster->GetXaxis()->GetNbins()/2.0 + 0.5); for (Int_t bini=1; bini <= HistogramTypepointer->histAvgCluster->GetXaxis()->GetNbins(); bini++) xslicetroughcluster->SetBinContent(bini,HistogramTypepointer->histAvgCluster->GetBinContent(bini,middlebin)); TH1F* yslicetroughcluster = new TH1F(Form("%s_yslice", HistogramTypepointer->histAvgCluster->GetTitle()),"Y slice",HistogramTypepointer->histAvgCluster->GetYaxis()->GetNbins(),HistogramTypepointer->histAvgCluster->GetYaxis()->GetBinLowEdge(1),HistogramTypepointer->histAvgCluster->GetYaxis()->GetBinUpEdge(HistogramTypepointer->histAvgCluster->GetYaxis()->GetNbins())); - Double_t yslicetroughclusterpar[4]; + Double_t* yslicetroughclusterpar = (Double_t *)calloc(4, sizeof(Double_t)); // allocate 6 parameters for safety, maximum 4 used at the moment middlebin = (int)(HistogramTypepointer->histAvgCluster->GetYaxis()->GetNbins()/2.0 + 0.5); for (Int_t bini=1; bini <= HistogramTypepointer->histAvgCluster->GetYaxis()->GetNbins(); bini++) yslicetroughcluster->SetBinContent(bini,HistogramTypepointer->histAvgCluster->GetBinContent(middlebin,bini)); @@ -1351,13 +1372,12 @@ Bool_t Run::plotClusterDistribution(HistogramType* HistogramTypepointer) canvas->cd(1); TPad* firstpart = (TPad*)canvas->GetPad(1); firstpart->SetGrid(); - HistogramTypepointer->FitPerform(xslicetroughcluster,"lorentz", false, &xslicetroughclusterpar[0]); - + xslicetroughclusterpar = HistogramTypepointer->FitPerform(xslicetroughcluster,"lorentz"); canvas->cd(2); TPad* secondpart = (TPad*)canvas->GetPad(2); secondpart->SetGrid(); - HistogramTypepointer->FitPerform(yslicetroughcluster,"lorentz", false, &yslicetroughclusterpar[0]); + yslicetroughclusterpar = HistogramTypepointer->FitPerform(yslicetroughcluster,"lorentz"); canvas->Draw(); @@ -1385,17 +1405,22 @@ void Run::plotVerticalLine(TH1F* histogrampointer, Float_t xVal) { l->SetLineStyle(2); // dashed l->Draw("same"); - TString legendEntry = TString(Form("peak position: %.1f",xVal )); + TString legendEntry; + if (abs(xVal) < 1) + legendEntry = TString(Form("peak position: %.2f",xVal )); + else + legendEntry = TString(Form("peak position: %.1f",xVal )); TLegend* leg = new TLegend(0.65,0.65,0.80,0.80);//(0.6,0.7,0.89,0.89); leg->SetFillColor(0); leg->SetBorderSize(0); + leg->SetFillStyle(0); leg->AddEntry((TObject*) 0, legendEntry, ""); leg->SetTextSize(0.025); leg->Draw(); } } -TCanvas* Run::plot1DHistogram(HistogramType* HistogramTypepointer, TH1F* onehistogram, TString fitFuncType, TString titlestr, TString legendstr) +TCanvas* Run::plot1DHistogram(HistogramType* HistogramTypepointer, TH1F* onehistogram, TString fitFuncType, Float_t fitstart, TString titlestr, TString legendstr) { if (onehistogram != nullptr) { @@ -1407,16 +1432,21 @@ TCanvas* Run::plot1DHistogram(HistogramType* HistogramTypepointer, TH1F* onehist TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700); onehistogram->SetTitle(titlestr); onehistogram->Draw(); - Float_t maxValue = HistogramTypepointer->FitPerform(onehistogram, fitFuncType, true); + Double_t* parameters = (Double_t *)calloc(6, sizeof(Double_t)); // allocate 6 parameters for safety, maximum 4 used at the moment + parameters = HistogramTypepointer->FitPerform(onehistogram, fitFuncType, true, fitstart); + Float_t maxValue = parameters[1]; plotVerticalLine(onehistogram, maxValue); - TLegend* leg = new TLegend(0.4,0.8,0.89,0.89);//(0.6,0.7,0.89,0.89); + Float_t startvalue = onehistogram->GetBinLowEdge(onehistogram->FindFirstBinAbove(0)); + Float_t endrange = onehistogram->GetBinCenter(onehistogram->FindLastBinAbove(0)); + onehistogram->GetXaxis()->SetRangeUser(startvalue,endrange); + TLegend* leg = new TLegend(0.2,0.85,0.89,0.89);//(0.6,0.7,0.89,0.89); leg->SetFillColor(0); leg->SetBorderSize(0); if (legendstr.Length() < 1) legendstr = Form("%s",onehistogram->GetTitle()); leg->AddEntry(onehistogram, legendstr, "l"); - leg->SetTextSize(0.03); + leg->SetTextSize(0.025); leg->Draw(); canvas -> SaveAs( savepathresults + "/" + runcode + " " + onehistogram->GetName() + ".eps"); diff --git a/MABS_run_analyzer/Run.h b/MABS_run_analyzer/Run.h index 7ff0f13..c0dd624 100644 --- a/MABS_run_analyzer/Run.h +++ b/MABS_run_analyzer/Run.h @@ -370,7 +370,7 @@ public: */ Bool_t plotClusterDistribution(HistogramType*); - TCanvas* plot1DHistogram(HistogramType* histogramtypepointer, TH1F* onehistogram, TString fitFuncType = "landau", TString titlestr = "", TString legendstr = ""); + TCanvas* plot1DHistogram(HistogramType* histogramtypepointer, TH1F* onehistogram, TString fitFuncType = "landau", Float_t fitstart=-1, TString titlestr = "", TString legendstr = ""); pixelinfo pixelinfoMi34[32]; diff --git a/MABS_run_analyzer/SetStyle.c b/MABS_run_analyzer/SetStyle.c new file mode 100644 index 0000000..d177108 --- /dev/null +++ b/MABS_run_analyzer/SetStyle.c @@ -0,0 +1,80 @@ +{ + //////////////////////////////////////////////////////////////////////////// + // // + // who: Pete Alonzi (lpa2a@virginia.edu) // + // what: rootlogon file (~/.rootlogon.C) // + // when: last update: 2009_June_11 // + // where: galileo.phys.virginia.edu // + // why: change the default root canvas to acceptable configuration // + // // + // This file is run when you open root. // + // Just place it in your home directory and it works automatically. // + // If you make a change to this file simply start // + // a new root session and the changes will be implemented. // + // // + // If you did not install roofit then there will be an error // + // message when you login, it can be ignored. // + // // + //////////////////////////////////////////////////////////////////////////// + + // Definition of new styles + TStyle *Nab_collaboration = new TStyle("Nab_Collaboration","How Nab_Collaboration likes it"); + TStyle *Nab_collaboration_work = new TStyle("Nab_Collaboration_Work","How Nab_Collaboration likes it with stats"); + + // Set the default parameters for style Nab_collaboration + Nab_collaboration->SetCanvasBorderMode(0); // Removes Canvas Border + Nab_collaboration->SetPadBorderMode(0); // Removes Pad Border + Nab_collaboration->SetCanvasColor(0); // Changes Canvas color to white + Nab_collaboration->SetPadColor(0); // Changes Pad color to white + Nab_collaboration->SetStatColor(0); // Changes Stats bg color to white + Nab_collaboration->SetTitleFillColor(0); // Changes Title bg color to white + Nab_collaboration->SetLabelFont(42,"xyz"); // Changes Label Font + Nab_collaboration->SetTitleFont(42,"xyz"); // Changes Axis Title Font + Nab_collaboration->SetPadTickX(1); // Sets tic marks on both horizontal axes + Nab_collaboration->SetPadTickY(1); // Sets tic marks on both vertical axes + Nab_collaboration->SetTickLength(0.018,"xyz"); // Sets tic length + Nab_collaboration->SetOptStat(0); // Turns off Statistics display + Nab_collaboration->SetOptTitle(0); // Turns off Title display + Nab_collaboration->SetHistLineWidth(2); // Changes Histogram Line width + Nab_collaboration->SetFrameBorderMode(0); // Removes the Frame Border + Nab_collaboration->SetFrameFillStyle(0); + Nab_collaboration->SetCanvasDefH(494); // Sets Default Canva Height + Nab_collaboration->SetCanvasDefW(800); // Sets Default Canvas Width + Nab_collaboration->SetPalette(1); + + + // The next for commands set the default margin size + // n.b. the margins do not take axis labels into account! grr!! + Nab_collaboration->SetPadTopMargin(0.1); // Set Margin Top + Nab_collaboration->SetPadRightMargin(0.1); // Set Margin Right + Nab_collaboration->SetPadBottomMargin(0.1); // Set Margin Bottom + Nab_collaboration->SetPadLeftMargin(0.1); // Set Margin Left + + //make a nice palette + int NRGBs = 7, NCont = 999; + Nab_collaboration->SetNumberContours(NCont); + Double_t stops[7] = { 0.00, 0.10, 0.25, 0.45, 0.60, 0.75, 1.00 }; + Double_t red[7] = { 1.00, 0.00, 0.00, 0.00, 0.97, 0.97, 0.10 }; + Double_t green[7] = { 1.00, 0.97, 0.30, 0.40, 0.97, 0.00, 0.00 }; + Double_t blue[7] = { 1.00, 0.97, 0.97, 0.00, 0.00, 0.00, 0.00 }; + TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont); + + // Make Nab_collaboration_work a copy of Nab_collaboration and remove stats + Nab_collaboration->Copy(*Nab_collaboration_work); + Nab_collaboration_work->SetOptStat(1); + Nab_collaboration_work->SetOptTitle(1); + Nab_collaboration_work->SetPadTopMargin(0.1); + Nab_collaboration_work->SetPadRightMargin(0.1); + + // Select Which Style to + gROOT->SetStyle("Nab_Collaboration_Work"); + // gROOT->SetStyle("Nab_Collaboration"); + + + // Use the RooFit functions + // using namespace RooFit; + + // Force Root to run Nab_collaboration/Nab_collaboration_work style + gROOT->ForceStyle(); + +}