From 89d2eb23bbf26b5e77145c128380aa9404e6167a Mon Sep 17 00:00:00 2001 From: Benjamin Linnik Date: Tue, 2 Jun 2015 14:28:52 +0200 Subject: [PATCH] Anylyzer: Added pegasus support, ignoring cluster around pixel with highly negative CDS value in dynamical noise calculation, added switching between SMA and MMA dynamical noise calculation, code cleanup and messing up --- MABS_run_analyzer/ChargeSpektrum.c | 30 +++- MABS_run_analyzer/MAPS.c | 224 +++++++++++++++++++++++------ MABS_run_analyzer/Run.c | 145 +++++++++++-------- MABS_run_analyzer/Run.h | 27 +++- 4 files changed, 309 insertions(+), 117 deletions(-) diff --git a/MABS_run_analyzer/ChargeSpektrum.c b/MABS_run_analyzer/ChargeSpektrum.c index 8aaa8ac..d43a9c2 100644 --- a/MABS_run_analyzer/ChargeSpektrum.c +++ b/MABS_run_analyzer/ChargeSpektrum.c @@ -132,6 +132,8 @@ void ChargeSpektrum(TString runnumber = "") /// number of runs to be analyzed, number of lines read by @c ReadRunList() ReadRunList(&runList); } + numberRuns=2; + runList.push_back(runList.front()); runs = new Run*[numberRuns]; cout << "Found " << numberRuns << " run(s) in 'runlist.txt' or given as parameters." << endl; @@ -162,8 +164,22 @@ void ChargeSpektrum(TString runnumber = "") } } runs[runi]->setResultsPath("./results/"); + runs[runi]->error=runs[runi]->initRun(); +// if (runi%2) +// { + if(runs[runi]->setDynamicalNoiseMode("simple")) cout << coloryellow << "Passed noise mode not recognized, using simple moving average" << endlr; +// } +// else +// { +// if(runs[runi]->setDynamicalNoiseMode("exponential")) cout << coloryellow << "Passed noise mode not recognized, using simple moving average" << endlr; +// } runs[runi]->useDynamicalNoise(true); - // runs[runi]->analyzeFrame(542875); // creates or opens .root file, can analyze the RAW data + +// runs[runi]->analyzeFrame(15684); +// runs[runi]->analyzeFrame(803316); +// runs[runi]->analyzeFrame(8319); + + // creates or opens .root file, can analyze the RAW data runs[runi]->error=runs[runi]->analyzeRun(true); // creates or opens .root file, can analyze the RAW data if (!runs[runi]->error) { @@ -171,15 +187,15 @@ void ChargeSpektrum(TString runnumber = "") // runs[runi]->plotSeed(); // runs[runi]->plotSeedThresholdCalibrated(); // runs[runi]->plotSeedThreshold(); - // runs[runi]->plotSeed(); + runs[runi]->plotSeed(); // runs[runi]->plotSum(); // runs[runi]->plotVeto(); // runs[runi]->plotNoise(); if (!isBatch) gROOT->SetBatch(kFALSE); -// runs[runi]->plotAllHistograms(); - runs[runi]->plotAllHistogramsThresholdCluster(); - runs[runi]->plotAllHistogramsThresholdClusterCalibrated(); + runs[runi]->plotAllHistograms(); +// runs[runi]->plotAllHistogramsThresholdCluster(); +// runs[runi]->plotAllHistogramsThresholdClusterCalibrated(); // runs[runi]->plotAllHistogramsCalibrated(); runs[runi]->writeAllHistogramsToFile(); } @@ -187,9 +203,9 @@ void ChargeSpektrum(TString runnumber = "") } } // plotAllRuns(""); -plotAllRuns("threshold"); + plotAllRuns("threshold calibrated"); // writeObservableToFile("seed threshold calibrated"); - writeObservableToFile("seed threshold"); + writeObservableToFile("seed threshold calibrated"); } Int_t ReadRunList(std::vector* runlist) diff --git a/MABS_run_analyzer/MAPS.c b/MABS_run_analyzer/MAPS.c index 8ee2102..acdbeb1 100644 --- a/MABS_run_analyzer/MAPS.c +++ b/MABS_run_analyzer/MAPS.c @@ -217,11 +217,6 @@ Bool_t MAPS::initMapsRun( ) { cout<<" Ordering according to : "<= 20000) cout << colorred << "Run could be corrupted, number of events in file doesn't match configuration! ("<< Frames << " != " << FileEvNbInConfig << ")" << endlr; // cout< highlightabove) std::cout << colorred; - std::cout.width(precision+4); + std::cout.width(precision+6); std::cout << std::fixed; std::cout << std::left << std::setprecision(precision); std::cout << a[i] << colorreset; @@ -989,7 +984,7 @@ bool MAPS::regetDynNoise(Int_t Frames) { { if (useexponentialdecayingnoisewindow) { - fNoise[pixeli] = ( (Frames-1)*fNoise[pixeli] + TMath::Power(1.*fCdsmatrix[pixeli]-fPedestals[pixeli],2) )/(Frames); + fNoise[pixeli] = TMath::Sqrt(( (Frames-2)*TMath::Power(fNoise[pixeli],2) + TMath::Power(1.*fCdsmatrix[pixeli]-fPedestals[pixeli],2) )/(Frames-1)); fPedestals[pixeli] = ( (Frames-1)*fPedestals[pixeli] + fCdsmatrix[pixeli] )/Frames; } else @@ -1000,43 +995,46 @@ bool MAPS::regetDynNoise(Int_t Frames) { fPedestals[pixeli] -= (CDSlastframes[pixeli].front()-fCdsmatrix[pixeli])/CDSlastframes[pixeli].size(); CDSlastframes[pixeli].pop(); CDSlastframes[pixeli].push(fCdsmatrix[pixeli]); - if (abs(fPedestals[pixeli])>20) - { - fPedestalhighPixel++; - if (!pedestalhighinthisframe) { - fPedestalhighFrames++; - pedestalhighinthisframe=true; - } -// cout<<"\rFrame: "< Pedestal suspiciously high!"<(fCdsmatrix , fPixels, fColumns, 1, 20); -// cout<<" Pedestial of sensor: "; -// debugStream<>(fPedestals , fPixels, fColumns, 1, 20); -// cout<<" Hitted Pixel discriminator matrix: "; -// debugStream<>(fHittedPixel, fPixels, fColumns, 1, 1); -// exit(1); + } + if (abs(fPedestals[pixeli])>20) + { + fPedestalhighPixel++; + if (!pedestalhighinthisframe) { + fPedestalhighFrames++; + pedestalhighinthisframe=true; } - if (abs(fNoise[pixeli])>20) - { - fNoiseHighPixel++; - if (!noisehighinthisframe) { - fNoiseHighFrames++; - noisehighinthisframe=true; - } + // cout<<"\rFrame: "< Pedestal suspiciously high!"<(fCdsmatrix , fPixels, fColumns, 1, 20); + // cout<<" Pedestial of sensor: "; + // debugStream<>(fPedestals , fPixels, fColumns, 1, 20); + // cout<<" Hitted Pixel discriminator matrix: "; + // debugStream<>(fHittedPixel, fPixels, fColumns, 1, 1); + // exit(1); + } + if (abs(fNoise[pixeli])>20) + { + fNoiseHighPixel++; + if (!noisehighinthisframe) { + fNoiseHighFrames++; + noisehighinthisframe=true; + } + // cout << colorred << fFrameNumber << endlr; +// if (fFrameNumber==15683 || fFrameNumber==15684) +// { // cout <<"\rFrame: "< Noise suspiciously high! "<(fCdsmatrix , fPixels, fColumns, 1, 20); +// debugStream<>(fCdsmatrix , fPixels, fColumns, 0, 20); // cout<<" Noise of sensor: "; -// debugStream<>(fNoise , fPixels, fColumns, 1, 20); +// debugStream<>(fNoise , fPixels, fColumns, 0, 20); // cout<<" Hitted Pixel discriminator matrix: "; -// debugStream<>(fHittedPixel, fPixels, fColumns, 1, 1); -// exit(1); - } +// debugStream<>(fHittedPixel, fPixels, fColumns, 0, 1); +// } } } else { -// cout << "Skipped pixel " << pixeli << " for noise calc" << endl; +// cout << "Skipped pixel " << pixeli << " for noise calc in frame " << fFrameNumber << ", because fHittedPixel: " << fHittedPixel[pixeli] << endl; } } if (pedestalhighinthisframe) @@ -1045,6 +1043,7 @@ bool MAPS::regetDynNoise(Int_t Frames) { cout<<"\rFrame: "< Noise suspiciously high!"< rechargePixellist; Int_t A; Int_t B; + Int_t seedrow; + Int_t seedcolumn; Int_t DUMMY=0; /// F0 - F1 of specific pixel in cluster Float_t pixelchargeincluster[25]={ 0 }; @@ -1190,7 +1193,7 @@ void MAPS::hitana() { Bool_t bordercluster=false; Float_t noisesumincluster=0; Float_t chargesumincluster=0; - + for(Int_t i=0; iFill(i%fColumns+0.1, (int)(i/fColumns)+0.1); // counts up in 2dimensional pixel matrix } } + + if( (float)(1.*fCdsmatrix[i]-fPedestals[i]) < (-5.*fNoise[i]) ) // loop to find pixel with highly negative CDS + { + rechargePixellist.push_back(i); + } } //Rewrite HITS to fHitlist array (why?) @@ -1216,6 +1222,131 @@ void MAPS::hitana() { Hitlist[i]=HITS.At(i); } +// Begin loop over all negative CDS pixel +// Determine clusters around them + for(Int_t rechargingpixeli=0; rechargingpixeli= (fRows-1)) ) { + // cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; + pixelchargeincluster[(row*5)+column]=DUMMY; + } + // wenn vierte Reihe des Clusters gefüllt werden soll + // und das Seed Pixel am Rand - in vorletzter Reihe des Matrix liegt + else if ( (row==4) && (seedrow>= (fRows-2)) ) { + // cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; + pixelchargeincluster[(row*5)+column]=DUMMY; + } + else + { + // wenn erste Spalte des Clusters gefüllt werden soll + // und das Seed Pixel am Rand - in erster oder zweiten Spalte des Matrix liegt + if ( (column==0) && (seedcolumn<2) ) { + // cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; + pixelchargeincluster[(row*5)+column]=DUMMY; + } + // wenn zweite Spalte des Clusters gefüllt werden soll + // und das Seed Pixel am Rand - in erster Spalte der Matrix liegt + else if ( (column==1) && (seedcolumn==0) ) { + // cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; + pixelchargeincluster[(row*5)+column]=DUMMY; + } + // wenn dritte Spalte des Clusters gefüllt werden soll + // und das Seed Pixel am Rand - in letzter Spalte der Matrix liegt + else if ( (column==3) && (seedcolumn==fColumns-1)) { + // cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; + pixelchargeincluster[(row*5)+column]=DUMMY; + } + // wenn letzte Spalte des Clusters gefüllt werden soll + // und das Seed Pixel am Rand - in vorletzter Spalte der Matrix liegt + else if ( (column==4) && (seedcolumn>fColumns-3)) { + // cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; + pixelchargeincluster[(row*5)+column]=DUMMY; + } + else { + pixelchargeincluster[(row*5)+column] = 1.*fCdsmatrix [rechargePixellist[rechargingpixeli]+(row-2)*fColumns+(column-2)] - fPedestals [rechargePixellist[rechargingpixeli]+(row-2)*fColumns+(column-2)]; + } + } + } + } + + //Check seeds (whether lowest entry in cluster): + for(Int_t i=6; i<19; i++) + { + if ( (i!=12) && (pixelchargeincluster[i] < pixelchargeincluster[12]) ) { + CHANCE=0; + break; + } + // maybe this is unnecessary, if upper if clause is a >= comparison + else if ( (i!=12) && (pixelchargeincluster[i] == pixelchargeincluster[12]) && i>12 ) { + cout << "WARNING: next pixel value identical to precessor" << endl; + CHANCE=0; //NOTE: potential error source + break; + } + else { + CHANCE=100; + if(i%5==3) { + i+=2; + }; + } + } + + //Begin: loop evaluate true seeds: + if(CHANCE==100) + { + if (seedcolumn < 2 || seedcolumn > fColumns-3 || seedrow < 2 || seedrow > fRows-3) + bordercluster = true; + else + bordercluster = false; + for(Int_t row=0; row<5; row++) + { + for(Int_t column=0; column<5; column++) + { + if ( (row==0) && (seedrow<2) ) { } + else if ( (row==1) && (seedrow<1) ) { } + else if ( (row==3) && (seedrow>= (fRows-1)) ) { } + else if ( (row==4) && (seedrow>= (fRows-2)) ) { } + else + { + if ( (column==0) && (seedcolumn<2) ) { } + else if ( (column==1) && (seedcolumn==0) ) { } + else if ( (column==3) && (seedcolumn==fColumns-1)) { } + else if ( (column==4) && (seedcolumn>fColumns-3)) { } + else { + fHittedPixel[rechargePixellist[rechargingpixeli]+(row-2)*fColumns+(column-2)] = -3; + } + } + } + } + // cout<<"Hitted pixel discriminator matrix:"<(fHittedPixel, fPixels, fColumns, 1, 1); + fHittedPixel[rechargePixellist[rechargingpixeli]] = -4; + + } + + } + //Begin: loop over all potential seed pixels: //Determine 'hit-vicinity': fHits = 0; @@ -1344,6 +1475,10 @@ void MAPS::hitana() { } } + + + + //Begin: loop evaluate true seeds: if(CHANCE==100) { @@ -1491,8 +1626,7 @@ void MAPS::plotFrame(Int_t FrameNumber) { cout<<"Hitted pixel discriminator matrix:"<(fHittedPixel, fPixels, fColumns, 1, 1); - TCanvas* cm1; - cm1 = new TCanvas("cm1","Matrix",50,100,1200,800); + TCanvas* cm1 = new TCanvas(TString(FrameNumber),TString(FrameNumber),50,100,1200,800); cm1->Divide(2,3); gStyle->SetOptFit(1011); gStyle->SetPalette(1); diff --git a/MABS_run_analyzer/Run.c b/MABS_run_analyzer/Run.c index 3fb01cd..d3ca0c1 100644 --- a/MABS_run_analyzer/Run.c +++ b/MABS_run_analyzer/Run.c @@ -61,8 +61,8 @@ Run::Run(Int_t runnumber, Int_t loopi) storepathRAWLinux = labbook.storepath; // storepathRAWLinux = storepathRAWLinux.ReplaceAll("H:","/jspc53_H"); // storepathRAWLinux = storepathRAWLinux.ReplaceAll("h:","/jspc53_H"); - storepathRAWLinux = storepathRAWLinux.ReplaceAll("H:","/d/jspc53_H"); - storepathRAWLinux = storepathRAWLinux.ReplaceAll("h:","/d/jspc53_H"); + storepathRAWLinux = storepathRAWLinux.ReplaceAll("H:","/d/garlic"); + storepathRAWLinux = storepathRAWLinux.ReplaceAll("h:","/d/garlic"); if (TString(hostname).EqualTo("jspc48")) { storepathRAWLinux = storepathRAWLinux.ReplaceAll("U:","/d/jspc28/jspc53_U"); @@ -126,7 +126,6 @@ Run::Run(Int_t runnumber, Int_t loopi) cout << colorwhite << "initHistograms():" << endlr; initHistograms(&histogram); initHistograms(&histogramthreshold, " threshold"); - runexistsinDB = 1; debugDBreadout(); } else @@ -195,72 +194,95 @@ Bool_t Run::setResultsPath(TString path) return 0; } -Bool_t Run::analyzeRun(Bool_t force) +Bool_t Run::setDynamicalNoiseMode(TString mode) +{ + cout<<" Dynamic noise calc : "; + if (mode.Contains("exponential", TString::kIgnoreCase) || mode.Contains("MMA", TString::kIgnoreCase)) + { + cout << "Modified moving average (MMA)" << endl; + processed->useexponentialdecayingnoisewindow = 1; + return 0; + } + else if (mode.Contains("simple", TString::kIgnoreCase) || mode.Contains("SMA", TString::kIgnoreCase)) + { + cout << "Simple moving average (SMA)" << endl; + processed->useexponentialdecayingnoisewindow = 0; + return 0; + } + return 1; +} + +Bool_t Run::initRun() { - if (runexistsinDB) + if (!error) { cout << colorwhite << "instatiate MAPS(this)" << endlr; processed = new MAPS(this); error = processed->error; - if (!error) - { - if (!runAllreadyAnalyzed() || force) - { - if (processed->initNewRootFile()) return 1; - /// progress meter, temporal variable - ULong_t progress_tmp=-1; - /// progress meter - Float_t progress; - processed->InitialDynNoise(); - int start = 0; - int nframes = processed->GetNumberFrames(); - for(int i=0; igetFrame(i); - // processed->filterCommonMode(); - // cout << "hitana " << endl; - processed->hitana(); - // cout << "regetDynNoise " << endl; - if (dynamicalNoise) - processed->regetDynNoise(); - 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; + return 0; + } + return 1; +} - // 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(); - processed->save(); - } - else + +Bool_t Run::analyzeRun(Bool_t force) +{ + if (!error) + { + if (!runAllreadyAnalyzed() || force) + { + if (processed->initNewRootFile()) return 1; + /// progress meter, temporal variable + ULong_t progress_tmp=-1; + /// progress meter + Float_t progress; + processed->InitialDynNoise(); + int start = 0; + int nframes = processed->GetNumberFrames(); + // for(int i=15580; i<15685;i++) // TODO remove 100000 run 342272 + for(int i=0; igetFrame(i); + // processed->filterCommonMode(); +// cout << "hitana " << i << endl; + processed->hitana(); +// cout << "regetDynNoise " << endl; + if (dynamicalNoise) + processed->regetDynNoise(); + 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; + + // 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(); + processed->save(); } - if (!error) + else { - if (processed->initOldRootFile()) return 1; - cout << colorwhite << "binNoise():" << endlr; - binNoise(); - cout << colorwhite << "binSeedSumVeto():" << endlr; - binSeedSumVeto(); - cout << colorwhite << "rescaleHistograms():" << endlr; - histogramCalibrated.calibrated = rescaleHistograms(); - histogramthresholdCalibrated.calibrated = histogramCalibrated.calibrated; - cout << colorwhite << "calculateCCE():" << endlr; - calculteCCE(); - cout << colorwhite << "updateDatabase():" << endlr; - updateDatabase(); - cout << colorwhite << "delete MAPS class:" << endlr; - delete processed; - return 0; + cout << "Skipped analysis of run " << labbook.runnumber << ", I think the root file for this run allready exists." << endl; + cout << colorwhite << "initOldRootFile():" << endlr; } + if (processed->initOldRootFile()) return 1; + cout << colorwhite << "binNoise():" << endlr; + binNoise(); + cout << colorwhite << "binSeedSumVeto():" << endlr; + binSeedSumVeto(); + cout << colorwhite << "rescaleHistograms():" << endlr; + histogramCalibrated.calibrated = rescaleHistograms(); + histogramthresholdCalibrated.calibrated = histogramCalibrated.calibrated; + cout << colorwhite << "calculateCCE():" << endlr; + calculteCCE(); + cout << colorwhite << "updateDatabase():" << endlr; + updateDatabase(); + cout << colorwhite << "delete MAPS class:" << endlr; delete processed; + return 0; } + delete processed; return 1; } @@ -347,7 +369,7 @@ void Run::rescaleHistogram(TH1F* &histogrampointernew, TH1F* &histogrampointerol Bool_t Run::analyzeFrame(Int_t frame) { - if (runexistsinDB) + if (!error) { processed = new MAPS(this); processed->initNewRootFile(); @@ -356,7 +378,7 @@ Bool_t Run::analyzeFrame(Int_t frame) { processed->InitialDynNoise(100); processed->getFrame(frame); - processed->filterCommonMode(); +// processed->filterCommonMode(); processed->hitana(); processed->plotFrame(frame); delete processed; @@ -469,6 +491,8 @@ void Run::setSystemSpecificParameters() systemparamcur = systemparamUSB; else if (labbook.system.EqualTo("PXI")) // && labbook.chipGen.EqualTo("34") ) systemparamcur = systemparamPXI; + else if (labbook.system.EqualTo("Pegasus")) // && labbook.chipGen.EqualTo("34") ) + systemparamcur = systemparamPegasus; setSensorInSystemParam(); } @@ -477,10 +501,13 @@ void Run::setSensorInSystemParam() { sensorinfostruct sensorinfoMi34USB( 8, 64 ); sensorinfostruct sensorinfoMi34PXI( 16, /* rows */ 64 /* columns */ ); + sensorinfostruct sensorinfoPegasus( 8, 56 ); if (labbook.system.EqualTo("USB")) // && labbook.chipGen.EqualTo("34") ) sensorinfocur=sensorinfoMi34USB; else if (labbook.system.EqualTo("PXI")) // && labbook.chipGen.EqualTo("34") ) sensorinfocur=sensorinfoMi34PXI; + else if (labbook.system.EqualTo("Pegasus")) // && labbook.chipGen.EqualTo("34") ) + sensorinfocur=sensorinfoPegasus; } template @@ -877,7 +904,7 @@ Float_t Run::FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verb Float_t posMax = 0; Float_t posMax2 = 0; Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax(); - Float_t noiseborder = posMaxValHist/10; // for USB system, the value is 90 + Float_t noiseborder = 90; // posMaxValHist/10 for USB system, the value is 90 if (doFits) { diff --git a/MABS_run_analyzer/Run.h b/MABS_run_analyzer/Run.h index 76d5402..a31fe6c 100644 --- a/MABS_run_analyzer/Run.h +++ b/MABS_run_analyzer/Run.h @@ -123,10 +123,7 @@ private: * The labbook data is used to make a string with useful information about the run */ Bool_t generateReadableRunCode(); - - /// is set to true if an entry was found in the SQL database - Bool_t runexistsinDB = 0; - + /// is set to true if an entry was found in the SQL database Bool_t runexistsAsRootFile = 0; @@ -208,8 +205,15 @@ private: Int_t nbinsnoise; }; systemparam systemparamUSB { - 900, // maxbin; - 900/5,// nbins; + 2800, // maxbin; + 2800/16,// nbins; + 25, //vetothreshold + 10, + 100 + }; + systemparam systemparamPegasus { + 2800, // maxbin; + 2800/2,// nbins; 25, //vetothreshold 10, 100 @@ -285,6 +289,12 @@ public: */ Bool_t analyzeRun(Bool_t force = 0); + + /** + * @brief initialize the RAW root file + */ + Bool_t initRun(); + /** * @brief analysis the RAW data of a given frame * @@ -315,6 +325,11 @@ public: */ Bool_t setResultsPath(TString path); + /** + * @brief set dynamical noise mode, exponential decaying noise windows (modified moving average), simple moving average (SMA) + */ + Bool_t setDynamicalNoiseMode(TString mode); + /** * @brief provides information from the SQL database */ -- 2.43.0