From f0f33bb3a676e3018e2d65bc04cda80926212f0b Mon Sep 17 00:00:00 2001 From: Benjamin Linnik Date: Mon, 24 Aug 2015 20:40:02 +0200 Subject: [PATCH] Run analyzer: Round up of average cluster calculation routine --- MABS_run_analyzer/ChargeSpektrum.c | 2 +- MABS_run_analyzer/MAPS.h | 28 +++---- MABS_run_analyzer/Run.c | 118 ++++++++++++++++++++++++++--- MABS_run_analyzer/Run.h | 40 +++++++--- 4 files changed, 153 insertions(+), 35 deletions(-) diff --git a/MABS_run_analyzer/ChargeSpektrum.c b/MABS_run_analyzer/ChargeSpektrum.c index 2375ddb..464e796 100644 --- a/MABS_run_analyzer/ChargeSpektrum.c +++ b/MABS_run_analyzer/ChargeSpektrum.c @@ -162,7 +162,7 @@ void ChargeSpektrum(TString runnumber = "") // runs[runi]->plotAllHistograms(); // runs[runi]->plotAllHistogramsThresholdCluster(); runs[runi]->plotClusterDistribution(&runs[runi]->histogramCalibrated); - runs[runi]->plotAllHistogramsThresholdClusterCalibrated(); +// runs[runi]->plotAllHistogramsThresholdClusterCalibrated(); // runs[runi]->plotAllHistogramsCalibrated(); // runs[runi]->writeAllHistogramsToFile(); } diff --git a/MABS_run_analyzer/MAPS.h b/MABS_run_analyzer/MAPS.h index 876c85c..9bbb504 100644 --- a/MABS_run_analyzer/MAPS.h +++ b/MABS_run_analyzer/MAPS.h @@ -69,7 +69,7 @@ private: Int_t fRunNumber; /// Number of Pixels found in Raw Data, is set im @c checkConf() Int_t fPixelsData; - /// Number of Pixels written to #fCdsmatrix and used for analysis + /// Number of Pixels written to #MAPS::fCdsmatrix and used for analysis Int_t fPixels; /// number of rows in the sensor, set and passed initMapsRun() to in the constructor, if Mimosa34 and system==USB then it is equal to 8, for PXI system it is equal to 16 Int_t fRows; @@ -94,12 +94,12 @@ private: int fNoiseDyn; Float_t HIT_tmp[29]; - /// total number of events, stated in the config file #fRootFile + /// total number of events, stated in the config file #MAPS::fRootFile UInt_t FileTotalEvNbInConfig = 0; - /// number of events per file, stated in the config file #fRootFile + /// number of events per file, stated in the config file #MAPS::fRootFile UInt_t FileEvNbInConfig = 0; - /// Average noise over all pixels in #fPixels in current frame to be binned into #fDynNoiseTree + /// Average noise over all pixels in #MAPS::fPixels in current frame to be binned into #MAPS::fDynNoiseTree Float_t fNoiseMean; /// Average pedestal to be binned into fDynNoiseTree Float_t fPedestalsMean; @@ -107,7 +107,7 @@ private: Int_t* fEvents; /// Array with f0 values of pixels in current frame Int_t* fF0matrix; - /// Array with f1 values of pixels in current frame, used to generate #fCdsmatrix = #fF0matrix - #fF1matrix + /// Array with f1 values of pixels in current frame, used to generate #MAPS::fCdsmatrix = #MAPS::fF0matrix - #MAPS::fF1matrix Int_t* fF1matrix; /// Array with CDS values of pixels in current frame Float_t* fCdsmatrix; @@ -223,7 +223,7 @@ private: * @brief Compares provided information with that found in the RAW data * @return true on success, false on failure * - * It also sets #FileTotalEvNbInConfig, #FileEvNbInConfig and @a fPixelsData accordingly + * It also sets #MAPS::FileTotalEvNbInConfig, #MAPS::FileEvNbInConfig and @a fPixelsData accordingly */ bool checkConf(Int_t&); @@ -286,7 +286,7 @@ private: void debugStream(const arraytype* (a), Int_t n=512, Int_t columns=8, Int_t precision=2, float highlightabove = 99999999); /** - * @brief Initialize histogram labels and histograms #hint1, #hint2, #fdiscriminatedhitmatrix, #fADCHitmatrix + * @brief Initialize histogram labels and histograms #MAPS::hint1, #MAPS::hint2, #MAPS::fdiscriminatedhitmatrix, #MAPS::fADCHitmatrix */ void initHistograms(); @@ -343,7 +343,7 @@ public: Bool_t initOldRootFile(); - /** @brief Refill the noise branch every #numberofframesfornoise frames or just use the initial initialization + /** @brief Refill the noise branch every #MAPS::numberofframesfornoise frames or just use the initial initialization * * turning this on will slow down the analysis, but allows time dependant noise investigation */ @@ -368,14 +368,14 @@ public: /** - * @brief Calculates a first estimate of the noise and pedestial of each pixel in #frames + * @brief Calculates a first estimate of the noise and pedestial of each pixel in #MAPS::frames * * Calculate noise and pedestial, trying to avoid pixels with hit charge * * Steps: * 1. Calculate noise and pedestial * 2. Recalculate noise and pedestial without pixels with higher CDS then 5 sigma - * 3. Recalculate noise and pedestial again using second estimate, without pixels with higher CDS then 3 sigma. This noise and pedestial is called #fPedestals and #fNoise. + * 3. Recalculate noise and pedestial again using second estimate, without pixels with higher CDS then 3 sigma. This noise and pedestial is called #MAPS::fPedestals and #MAPS::fNoise. * This new noise and pedestial will be used for hit analysis (@c hitana()) * * @return true if no errors occured @@ -387,7 +387,7 @@ public: /** * @brief Eliminates line structures in the chip * - * @Todo An average over the noise of a line is calculated and subtracted from each #fCdsmatrix value of every pixel in this line.
+ * @Todo An average over the noise of a line is calculated and subtracted from each #MAPS::fCdsmatrix value of every pixel in this line.
* Noise structure between different lines is eliminated this way. * * @see fCdsmatrix @@ -448,7 +448,7 @@ public: Long_t plus = 0; Long_t minus = 0; /** - * @brief Analysis current frame for hits, saves cluster to #fFrameInfo.p + * @brief Analysis current frame for hits, saves cluster to #MAPS::fFrameInfo.p * * loops over all pixel, compares * @@ -457,8 +457,8 @@ public: * * if true, saves pixel as seed pixel candidate, creates 5x5 cluster * around it. If cluster has biggest value in the middle, then count it as a valid - * cluster and indicate a hit in the discriminator matrix #fHittedPixel, TH2F variable #fdiscriminatedhitmatrix - * and TH2F ADC matrix #fADCHitmatrix + * cluster and indicate a hit in the discriminator matrix #MAPS::fHittedPixel, TH2F variable #MAPS::fdiscriminatedhitmatrix + * and TH2F ADC matrix #MAPS::fADCHitmatrix * */ void hitana (); diff --git a/MABS_run_analyzer/Run.c b/MABS_run_analyzer/Run.c index ecb94f0..568d079 100644 --- a/MABS_run_analyzer/Run.c +++ b/MABS_run_analyzer/Run.c @@ -784,8 +784,6 @@ Bool_t Run::binSeedSumVeto() Bool_t Run::binCluster() { - // histogram.histAvgCluster->SetBit(TH1::kCanRebin); - // histogramthreshold.histAvgCluster->SetBit(TH1::kCanRebin); Float_t rotateangle = 0; if (curpixelinfo.staggered) { rotateangle = TMath::ATan(curpixelinfo.pitchY/(curpixelinfo.pitchX/2))*180.0/TMath::Pi(); @@ -1209,8 +1207,7 @@ Bool_t Run::plotClusterDistribution(histogramstruct* histogramstructpointer) Canvas_1->Modified(); Canvas_1->cd(); Canvas_1->SetSelected(Canvas_1); - - + TImageDump *img = new TImageDump(savepathresults + "/" + canvastitle + ".png"); Canvas_1->Paint(); img->Close(); @@ -1225,7 +1222,54 @@ Bool_t Run::plotClusterDistribution(histogramstruct* histogramstructpointer) // gStyle->SetPaperSize(10.,10.); // canvas->Print(savepathresults + "/" + canvastitle + ".tex"); - return 0; + // create lorentz fit of a slice + TH1F* xslicetroughcluster = new TH1F(Form("%s_xslice", histogramstructpointer->histAvgCluster->GetTitle()),"X slice",histogramstructpointer->histAvgCluster->GetXaxis()->GetNbins(),histogramstructpointer->histAvgCluster->GetXaxis()->GetBinLowEdge(1),histogramstructpointer->histAvgCluster->GetXaxis()->GetBinUpEdge(histogramstructpointer->histAvgCluster->GetXaxis()->GetNbins())); + Double_t xslicetroughclusterpar[4]; + Int_t middlebin = (int)(histogramstructpointer->histAvgCluster->GetXaxis()->GetNbins()/2.0 + 0.5); + for (Int_t bini=1; bini <= histogramstructpointer->histAvgCluster->GetXaxis()->GetNbins(); bini++) + xslicetroughcluster->SetBinContent(bini,histogramstructpointer->histAvgCluster->GetBinContent(bini,middlebin)); + + TH1F* yslicetroughcluster = new TH1F(Form("%s_yslice", histogramstructpointer->histAvgCluster->GetTitle()),"Y slice",histogramstructpointer->histAvgCluster->GetYaxis()->GetNbins(),histogramstructpointer->histAvgCluster->GetYaxis()->GetBinLowEdge(1),histogramstructpointer->histAvgCluster->GetYaxis()->GetBinUpEdge(histogramstructpointer->histAvgCluster->GetYaxis()->GetNbins())); + Double_t yslicetroughclusterpar[4]; + middlebin = (int)(histogramstructpointer->histAvgCluster->GetYaxis()->GetNbins()/2.0 + 0.5); + for (Int_t bini=1; bini <= histogramstructpointer->histAvgCluster->GetYaxis()->GetNbins(); bini++) + yslicetroughcluster->SetBinContent(bini,histogramstructpointer->histAvgCluster->GetBinContent(middlebin,bini)); + + random = random1->Rndm()*100000; + canvastitle = Form("%s Cluster slices", runcode.Data()); + if (histogramstructpointer->calibrated) + canvastitle += "_el"; + if (histogramstructpointer->thresholdcluster) + canvastitle += "_trsh"; + canvasname = Form("%s%d",runcode.Data(),random); + TCanvas* canvas = new TCanvas(canvasname, canvastitle, 1200, 700); + canvas->Divide(2,1); + canvas->SetFillColor(0); + canvas->SetBorderSize(2); + + canvas->cd(1); + TPad* firstpart = (TPad*)canvas->GetPad(1); + firstpart->SetGrid(); + FitPerform(xslicetroughcluster,"lorentz", false, &xslicetroughclusterpar[0]); + + + canvas->cd(2); + TPad* secondpart = (TPad*)canvas->GetPad(2); + secondpart->SetGrid(); + FitPerform(yslicetroughcluster,"lorentz", false, &yslicetroughclusterpar[0]); + + canvas->Draw(); + + TImageDump *img2 = new TImageDump(savepathresults + "/" + canvastitle + ".png"); + canvas->Paint(); + img2->Close(); + + f->Append(canvas); + f->Append(img2); + f->Write(); + // .Data(), Form("Avg. cluster distribution; x; y %s", humanreadablestr.Data()),sizearrrotx, xcoord_abs_min,xcoord_abs_max,sizearrroty,ycoord_abs_min,ycoord_abs_max); + + return 0; } Bool_t Run::plotClusterDistribution() @@ -1346,7 +1390,7 @@ Bool_t Run::integrateSr90Spectra(TH1F* histogrampointer, Bool_t verbose) return 0; } -Float_t Run::FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verbose) +Float_t Run::FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verbose, Double_t* parameters) { Float_t posMax = 0; Float_t posMax2 = 0; @@ -1355,12 +1399,12 @@ Float_t Run::FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verb if (doFits) { - 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()); - TF1* fitFunc = new TF1("fitFunc",fitFuncType,noiseborder,posMaxValHist); - if (fitFuncType.Contains("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()); + TF1* fitFunc = new TF1("fitFunc",fitFuncType,noiseborder,posMaxValHist); + if (TString(histogrampointer->GetName()).Contains("Veto")) { Float_t peak1 = histogrampointer->GetMaximumBin(); @@ -1445,6 +1489,10 @@ Float_t Run::FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verb } else if (fitFuncType=="landau") { + 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()); + TF1* fitFunc = new TF1("fitFunc",fitFuncType,noiseborder,posMaxValHist); + Float_t fitMax1 = 1000; Float_t fitMax2 = 1000; Float_t fitMax3 = 1000; @@ -1472,12 +1520,59 @@ Float_t Run::FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verb //fitLandauErrorLeft.push_back(posMax - minFitMax); //fitLandauErrorRight.push_back(maxFitMax - posMax); + } else if (fitFuncType=="lorentz") + { + // create a TF1 with the range from 0 to 3 and 6 parameters + TF1 *fitFcn = new TF1("fitFcn",lorentzianPeak,0,160,4); + fitFcn->SetNpx(1000); + fitFcn->SetLineWidth(4); + fitFcn->SetLineColor(rootcolors[plotStyle]); + // set start values for some parameters + fitFcn->SetParName(0,"Area"); + fitFcn->SetParameter(0,21655); // Area + fitFcn->SetParName(1,"Width"); + fitFcn->SetParameter(1,34); // width + fitFcn->SetParName(2,"Peak pos."); + fitFcn->SetParameter(2,80); // peak position x + fitFcn->SetParName(3,"Y shift"); + fitFcn->SetParameter(3,-22); // y-shift + gStyle->SetOptFit(0111); + +// fitFcn->FixParameter(0,21655.90045); +// fitFcn->FixParameter(1,34.31885101); + fitFcn->FixParameter(2,histogrampointer->GetXaxis()->GetBinCenter((int)(histogrampointer->GetXaxis()->GetNbins()/2.0+0.5))); +// fitFcn->FixParameter(3,-22.43016284); + + //histogrampointer->Fit("fitFcn","V+","ep"); + histogrampointer->Fit("fitFcn","Q,W","ep"); + Double_t par[4]; + + // writes the fit results into the par array + fitFcn->GetParameters(par); + for (Int_t pari=0; pari<4; pari++) + parameters[pari]=fitFcn->GetParameter(pari); + if (verbose) + { + Printf("width: %.6f ",fitFcn->GetParameter(1)); + cout << "width: " << fitFcn->GetParameter(1) << endl; + cout << "peak: " << fitFcn->GetParameter(2) << endl; + cout << "y-shift:" << fitFcn->GetParameter(3) << endl; + } } + } return posMax; } +// Lorenzian Peak function +Double_t Run::lorentzianPeak(Double_t *x, Double_t *par) { + return par[3]+(0.5*par[0]*par[1]/TMath::Pi()) / + TMath::Max( 1.e-10,(x[0]-par[2])*(x[0]-par[2]) + + .25*par[1]*par[1]); +} + + void Run::plotVerticalLine(TH1F* histogrampointer, Float_t xVal) { if (xVal > 0) { @@ -1576,7 +1671,7 @@ Bool_t Run::write2DHistogramToFile(TH2F* onehistogram) for (Int_t biny_i = 1; biny_i<=onehistogram->GetXaxis()->GetNbins(); biny_i++) { for (Int_t binx_i = 1; binx_i<=onehistogram->GetXaxis()->GetNbins(); binx_i++) { outline=Form("%.2f\t%.2f\t%.2f", onehistogram->GetXaxis()->GetBinCenter(binx_i), onehistogram->GetYaxis()->GetBinCenter(biny_i), onehistogram->GetBinContent(binx_i, biny_i)); - cout << outline << endl; // debug line + //cout << outline << endl; // debug line *fout< #include +#include +#include + #include #include #include @@ -107,6 +110,10 @@ private: */ string removeNumbers(std::string x); + /** + * @brief Lorenzian Peak function */ + static Double_t lorentzianPeak(Double_t *x, Double_t *par); + /** * @brief removes all punctuation from given string */ @@ -144,13 +151,13 @@ private: char hostname[1024]; /** - * @brief fills noise #histogram */ + * @brief fills noise histograms in #Run::histogramstruct */ Bool_t binNoise(); /** - * @brief fills Seedm, Sum and Veto #histogram */ + * @brief fills Seedm, Sum and Veto #Run::histogramstruct */ Bool_t binSeedSumVeto(); /** - * @brief fills the average cluster distribution into a 2 dimension histogram #histogram*/ + * @brief fills the average cluster distribution into a 2 dimension histogram #Run::histogramstruct */ Bool_t binCluster(); /// noise quantiles: mean value, sigma in postive direction and sigma in negative direction Double_t noisequantiles[3]; @@ -174,7 +181,7 @@ private: * @return peak position of the fit * */ - Float_t FitPerform(TH1F*, TString fitFuncType="landau", Bool_t verbose=true); + Float_t FitPerform(TH1F*, TString fitFuncType="landau", Bool_t verbose=true, Double_t* parameters=NULL); /** * @brief find the border between the noise and the signal in Sr90 runs @@ -323,7 +330,7 @@ public: /** @brief Writes a given histogram into a file */ Bool_t write2DHistogramToFile(TH2F* onehistogram); - /** @brief Writes all histogram of @c Run::histogramsctruct into a file */ + /** @brief Writes all histogram of @c Run::histogramstruct into a file */ Bool_t writeAllHistogramsToFile(); /** @brief Makes a gnuplot file to plot the histogram data created in @c Run::writeAllHistogramsToFile() */ @@ -502,16 +509,31 @@ public: histogramstruct* plothistogramstructpointer; TH1F** plothistogrampointer; - /** @brief Plots all histograms from #histogram into one canvas */ + /** @brief Plots all histograms from #Run::histogramstruct into one canvas + * + * @param histogramstruct pointer to a histogram structure of type #Run::histogramstruct + */ Bool_t plotAllHistograms(histogramstruct*); - /** @brief Plots average cluster charge distribution */ + /** @brief Plots average cluster charge distribution + * + * Creates a graphical view of the average charge distribution in a + * cluster + * + * @param histogramstruct Pointer to a structure of type #Run::histogramstruct + * + */ Bool_t plotClusterDistribution(histogramstruct*); - /// init all histograms, set binwidth, bin amount and names + /// + /** @brief initializes all histograms, set binwidth, bin amount and names + * * + * @param histogramstruct Pointer to a structure of type #Run::histogramstruct + * @param suffix A string which will be appended to the generated title of each histograms + */ void initHistograms(histogramstruct*, TString suffix = ""); - /** @brief initializes the TH2F cluster for the binning for a specific structure of type #histogramstruct one points to*/ + /** @brief initializes the TH2F cluster for the binning for a specific structure of type #Run::histogramstruct one points to*/ void initCluster(histogramstruct*, TString suffix, Float_t xcoord_min_step, Float_t xcoord_abs_min, Float_t xcoord_abs_max, Float_t ycoord_min_step, Float_t ycoord_abs_min, Float_t ycoord_abs_max); pixelinfo pixelinfoMi34[32]; -- 2.43.0