From: Benjamin Linnik Date: Sat, 21 Mar 2015 04:01:37 +0000 (+0100) Subject: Run analyzer: Code cleanup X-Git-Url: https://jspc29.x-matter.uni-frankfurt.de/git/?a=commitdiff_plain;h=ec25c633674281b2f00e6f3893bb222e43783d9f;p=radhard.git Run analyzer: Code cleanup --- diff --git a/MABS_run_analyzer/ChargeSpektrum.c b/MABS_run_analyzer/ChargeSpektrum.c index 62d2b70..ffef131 100644 --- a/MABS_run_analyzer/ChargeSpektrum.c +++ b/MABS_run_analyzer/ChargeSpektrum.c @@ -28,6 +28,8 @@ Int_t* ReadRunList(Int_t*); +Run** runs; + void ChargeSpektrum() { cout << endl << endl; @@ -36,13 +38,14 @@ void ChargeSpektrum() Int_t numberRuns=0; ReadRunList(&numberRuns); /// array with run numbers - Int_t* runList=new Int_t[numberRuns]; - Run* runs[numberRuns]; + Int_t* runList=new Int_t[numberRuns]; + runList=ReadRunList(&numberRuns); + runs = new Run*[numberRuns]; cout << "Found " << numberRuns << " run(s) in 'runlist.txt'." << endl; for(Int_t runi=0;runiupperpart = 0; } } + runs[runi]->setResultsPath("/home/blinnik/Dropbox/Promotion/MABS_Analyse/HR18vs20/"); + runs[runi]->setPlotStyle(runi); + runs[runi]->useDynamicalNoise(true); runs[runi]->analyzeRun(true); // creates or opens .root file, can analyze the RAW data -// runs[runi]->analyzeFrame(100000); -// runs[runi]->analyzeFrame(86000); -// runs[runi]->analyzeFrame(0); -// runs[runi]->processed->plotFrame(39); + gROOT->SetBatch(kTRUE); + runs[runi]->plotSeed(); + runs[runi]->plotSum(); + runs[runi]->plotVeto(); + runs[runi]->plotNoise(); + gROOT->SetBatch(kFALSE); + runs[runi]->plotAllHistograms(); + runs[runi]->writeAllHistogramsToFile(); } + + TCanvas* canvas = new TCanvas("Summary", "Summary", 800, 600); + canvas->SetGridy(kTRUE); + canvas->SetGridx(kTRUE); + TLegend* leg = new TLegend(0.5,0.8,0.89,0.89); + leg->SetFillColor(0); + leg->SetBorderSize(0); + for(Int_t runi=0;runihistogram.Seed->Draw("SAME"); + leg->AddEntry(runs[runi]->histogram.Seed, runs[runi]->histogram.Seed->GetTitle(), "l"); + } + leg->SetTextSize(0.03); + leg->Draw(); + canvas->Modified(); + canvas->cd(); + canvas->SetSelected(canvas); + +// canvas -> SaveAs( savepathresults + "/" + runcode + " " + histogram->GetName() + ".eps"); +// +// TImage *img = TImage::Create(); +// img->FromPad(canvas); +// img->WriteImage(savepathresults + "/" + runcode + " " + histogram->GetName() + ".png"); +// +// TFile *f = new TFile(savepathresults + "/" + runcode + " " + histogram->GetName() + ".root","UPDATE"); +// f->WriteTObject(canvas); +// f->WriteTObject(img); + + + +// cout << "Bevor delete! " << endl; +// for(Int_t runi=0;runilabbook.system = Form("USB"); } - run->setSystemInfo(); + run->setSystemSpecificParameters(); initMapsRun( ); } sleep(1); // TODO test if removable @@ -42,55 +42,66 @@ MAPS::MAPS(Run* runp) { Bool_t MAPS::initNewRootFile() { + fSave = true; if ( checkConf(fPixelsData) ) // TODO FileEvNbInConfig { - //----------------------------------------------- - //Check and open Data Files - int MaxFiles = FileTotalEvNbInConfig/FileEvNbInConfig; - if( checkDataFiles(MaxFiles) ) - { - fOutputFile = new TFile(fRootFile,"RECREATE"); - // Hit TTree - fHitTree = new TTree("hit", "hit"); - fHitTree->Branch("frame" , &fFrameInfo.frame , "frame/i" , 32000); - fHitTree->Branch("hits" , &fFrameInfo.hits , "hits/i" , 32000); - fHitTree->Branch("pixel" , &fFrameInfo.pixel[0] , "pixel[hits]/i" , 32000); - for(int i=0; i<25; i++) { - fHitTree->Branch( Form("p%i",i+1) , &fFrameInfo.p[i][0] , Form("p%i [hits]/F",i+1) , 32000); - } - // Noise TTree - fNoiseTree = new TTree("noise", "noise"); - fDynNoiseTree = fNoiseTree; - cout<<"-----------------------"<Branch("frame" , &fFrameInfo.frame , "frame/i" , 32000); + fHitTree->Branch("hits" , &fFrameInfo.hits , "hits/i" , 32000); + fHitTree->Branch("pixel" , &fFrameInfo.pixel[0] , "pixel[hits]/i" , 32000); + for(int i=0; i<25; i++) { + fHitTree->Branch( Form("p%i",i+1) , &fFrameInfo.p[i][0] , Form("p%i [hits]/F",i+1) , 32000); + } + // Noise TTree + fNoiseTree = new TTree("noise", "noise"); + fDynNoiseTree = fNoiseTree; + initHistograms(); + cout<<"-----------------------"<IsZombie()) { - cout << "Error opening old ROOT file!" << endl; + cout << colorred << "Error opening old ROOT file!"<< endlr; return 1; } fHitTree = (TTree*) fOutputFile->Get("hit"); + if (!fHitTree) // got pointer + { + cout<< colorred << "Error finding hit TTree in old ROOT file!"<< endlr; + return 1; + } fHitTree->SetBranchAddress("frame" , &fFrameInfo.frame ); fHitTree->SetBranchAddress("hits" , &fFrameInfo.hits ); fHitTree->SetBranchAddress ("pixel", &fFrameInfo.pixel[0]); for(int i=0; i<25; i++) { fHitTree->SetBranchAddress( Form("p%i",i+1) , &fFrameInfo.p[i][0]); } - - fHitTree = (TTree*) fOutputFile->Get("noise"); + fNoiseTree = (TTree*) fOutputFile->Get("noise"); + if (!fNoiseTree) // got pointer + { + cout<< colorred << "Error finding noise TTree in old ROOT file!"<< endlr; + return 1; + } fDynNoiseTree = fNoiseTree; - + initHistograms(); cout<<"-----------------------"<cd(); @@ -138,15 +148,15 @@ MAPS::~MAPS(void) { { fInn[i].close(); } - + + cout << "Bevor delete MAPS Arrays! " << endl; delete[] fEvents; delete[] fF0matrix; delete[] fF1matrix; delete[] fCdsmatrix; delete[] fNoise; delete[] fPedestals; - delete[] fDynFrameArr; - delete[] fDynCounter; + cout << "Nach delete MAPS Arrays! " << endl; cout<<"================================================================="<storepathLinux; - fOutDir = run->storepathLinux; // default ouput directory is input directory + fInDir=run->storepathRAWLinux; + fOutDir = run->storepathRAWLinux; // default ouput directory is input directory fRunNumber = run->labbook.runnumber; - fRows = run->sensorinfocurrent.rows; - fColumns = run->sensorinfocurrent.columns; + fRows = run->sensorinfocur.rows; + fColumns = run->sensorinfocur.columns; fPixels = fRows*fColumns; fMatrix = run->labbook.matrix; fSystem = run->labbook.system; @@ -172,13 +182,20 @@ void MAPS::initMapsRun( ) { cout<<" == Initiate MAPS() ... == "< "<Reset(); fNoiseTree->Branch("pixel" , &PIXEL , "pixel/i" , 32000); fNoiseTree->Branch("noise" , &NOISE , "noise/F" , 32000); fNoiseTree->Branch("pedestal" , &PEDESTAL , "pedestal/F" , 32000); @@ -665,7 +681,7 @@ bool MAPS::getNoise(Int_t Start, Int_t Frames) { NOISE+=TMath::Power(ARR[j*fPixels+i]-PEDESTAL,2); } - NOISE = TMath::Sqrt(NOISE/(Frames)); + NOISE = TMath::Sqrt(NOISE/(Frames-1)); fNoise[i] = NOISE; fPedestals[i] = PEDESTAL; @@ -691,11 +707,6 @@ bool MAPS::getNoise(Int_t Start, Int_t Frames) { //#################################################################### // initial Dynamic Noise + Pedestals calculation bool MAPS::InitialDynNoise(Int_t startframe, Int_t frames) { - cout<<" Dynamic noise calc : "; - if (useexponentialdecayingnoisewindow) - cout << "Modified moving average (MMA)" << endl; - else - cout << "Simple moving average (SMA)" << endl; // temporal local variables /// temporal variable to save second pedestial estimate Float_t* pixelpedestal2 = new Float_t[fPixels](); @@ -703,6 +714,7 @@ bool MAPS::InitialDynNoise(Int_t startframe, Int_t frames) { Float_t* pixelnoise2 = new Float_t[fPixels](); noiselastframes = new std::queue[fPixels](); + CDSlastframes = new std::queue[fPixels](); /// number of frames cut away in second noise estimate (propably hit frames) Int_t* nframescutawaysecondnoiseestimate = new Int_t[fPixels](); @@ -714,9 +726,9 @@ bool MAPS::InitialDynNoise(Int_t startframe, Int_t frames) { Float_t curpixelpedestal; Float_t curpixelnoise; - fDynNoiseTree->Branch("frame" , &fFrameNumber , "frame/i" , 32000); - fDynNoiseTree->Branch("noise" , &fNoiseMean , "noise/F" , 32000); - fDynNoiseTree->Branch("pedestal" , &fPedestalsMean , "pedestal/F" , 32000); +// fDynNoiseTree->Branch("frame" , &fFrameNumber , "frame/i" , 32000); +// fDynNoiseTree->Branch("noise" , &fNoiseMean , "noise/F" , 32000); +// fDynNoiseTree->Branch("pedestal" , &fPedestalsMean , "pedestal/F" , 32000); // check if all RAW data files are OK if(fOk) @@ -835,7 +847,7 @@ bool MAPS::InitialDynNoise(Int_t startframe, Int_t frames) { { fPedestals[pixeli] += fCdsmatrix[pixeli]; if (!useexponentialdecayingnoisewindow) - noiselastframes[pixeli].push(fCdsmatrix[pixeli]); + CDSlastframes[pixeli].push(fCdsmatrix[pixeli]); } else { @@ -858,6 +870,8 @@ bool MAPS::InitialDynNoise(Int_t startframe, Int_t frames) { if (abs(1.*(fCdsmatrix[pixeli]-pixelpedestal2[pixeli])) < 5.*pixelnoise2[pixeli]) { fNoise[pixeli] += TMath::Power(1.*fCdsmatrix[pixeli]-fPedestals[pixeli],2); + if (!useexponentialdecayingnoisewindow) + noiselastframes[pixeli].push(TMath::Power(1.*fCdsmatrix[pixeli]-fPedestals[pixeli],2)); } } } @@ -869,20 +883,25 @@ bool MAPS::InitialDynNoise(Int_t startframe, Int_t frames) { npixelscutawaythirdnoiseestimate = SumOverArray(nframescutawaythirdnoiseestimate, fPixels); cout << "In third estimate of dynamic noise, cutted away " << npixelscutawaythirdnoiseestimate << " pixels in total in " << frames << " frames. (" << (npixelscutawaythirdnoiseestimate*100./fPixels/frames) << " %)" << endl; - if( fSave ) { - fFrameNumber = 0; - fNoiseMean = 0; - fPedestalsMean = 0; - - fDynNoiseTree->Fill(); + Float_t PEDESTAL; + Float_t NOISE; + Int_t PIXEL; - fFrameNumber = 0; - fNoiseMean = (float)TMath::Mean((const int)fPixels, fNoise); - fPedestalsMean = (float)TMath::Mean((const int)fPixels, fPedestals); + fNoiseTree->Branch("pixel" , &PIXEL , "pixel/i" , 32000); + fNoiseTree->Branch("noise" , &NOISE , "noise/F" , 32000); + fNoiseTree->Branch("pedestal" , &PEDESTAL , "pedestal/F" , 32000); - fDynNoiseTree->Fill(); + for(Int_t pixeli=0; pixeliFill(); + } +// getNoise(startframe, frames); } fNoiseOk = true; @@ -904,19 +923,21 @@ arraytype MAPS::SumOverArray(const arraytype* (a), Int_t n) { template -arraytype MAPS::debugStream(const arraytype* (a), Int_t n, Int_t columns, Int_t precision) { - cout << endl << endl; - cout << "\033[1;29m----------------------------------\033[0m" << endl; +arraytype MAPS::debugStream(const arraytype* (a), Int_t n, Int_t columns, Int_t precision, float highlightabove) { + cout << colorwhite << "----------------------------------" << colorreset << endl; for (int i=0; i highlightabove) + std::cout << colorred; std::cout.width(precision+4); std::cout << std::fixed; - std::cout << std::left << std::setprecision(precision) << a[i]; + std::cout << std::left << std::setprecision(precision); + std::cout << a[i] << colorreset; if ((i+1)%columns==0) cout << endl << endl; } - cout << "\033[1;29m----------------------------------\033[0m" << endl; - cout << endl << endl; + cout << colorwhite << "----------------------------------" << colorreset << endl; + cout << endl; } @@ -1078,7 +1099,8 @@ bool MAPS::getDynNoise2(Int_t Frames) { bool MAPS::regetDynNoise(Int_t Frames) { for(Int_t pixeli=0; pixeli20) + { + cout<<"\rFrame: "< Pdestial 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) + { + cout <<"\rFrame: "< Noise suspiciously high! "<(fCdsmatrix , fPixels, fColumns, 1, 20); +// cout<<" Noise of sensor: "; +// debugStream<>(fNoise , fPixels, fColumns, 1, 20); +// cout<<" Hitted Pixel discriminator matrix: "; +// debugStream<>(fHittedPixel, fPixels, fColumns, 1, 1); +// exit(1); + } } } else { // cout << "Skipped pixel " << pixeli << " for noise calc" << endl; - } - - if (fFrameNumber%numberofframesfornoise == 0) + } + } + if (fFrameNumber%Frames == 0 && RefillNoiseBranch) + { + if(fSave) { - if(fSave) + Float_t PEDESTAL; + Float_t NOISE; + Int_t PIXEL; + +// fNoiseTree->Reset(); + fNoiseTree->Branch("pixel" , &PIXEL , "pixel/i" , 32000); + fNoiseTree->Branch("noise" , &NOISE , "noise/F" , 32000); + fNoiseTree->Branch("pedestal" , &PEDESTAL , "pedestal/F" , 32000); + + for(Int_t pixeli=0; pixeliFill(); + + PEDESTAL = fPedestals[pixeli]; + NOISE = fNoise[pixeli]; + PIXEL = pixeli; + fNoiseTree->Fill(); } + RefilledNoiseBranch = true; +// for(Int_t pixeli=0; pixeliFill(); +// } +// +// cout << fNoiseMean << endl; +// fNoiseMean = TMath::Mean((const int)fPixels, fNoise); +// fPedestalsMean = TMath::Mean((const int)fPixels, fPedestals); +// fDynNoiseTree->Fill(); } } // debug // if ( fFrameNumber%10000==0) // { -// debugStream<>(fPedestals, fPixels, fColumns, 2); +// debugStream<>(fNoise, fPixels, fColumns, 2, 5); // } return true; } @@ -1210,6 +1282,7 @@ void MAPS::hitana() { Int_t DUMMY=0; Float_t CLUSTER[25]; Int_t CHANCE; + Bool_t bordercluster=false; for(Int_t i=0; i Remove hits where the seed pixel is within 2 pixels beside the edge // Hitlist[hit]%fColumn = x-coordinate of the seed pixel, Hitlist[hit]/fColumns = y-coordinate of the seed pixel - if (Hitlist[hit]%fColumns < 2 || Hitlist[hit]%fColumns > fColumns-3 || Hitlist[hit]/fColumns < 2 || Hitlist[hit]/fColumns > fRows-3) - continue; +// if (Hitlist[hit]%fColumns < 2 || Hitlist[hit]%fColumns > fColumns-3 || Hitlist[hit]/fColumns < 2 || Hitlist[hit]/fColumns > fRows-3) +// continue; + //Provide 5x5 clusters with CDS - content: for(Int_t row=0; row<5; row++) @@ -1258,25 +1332,25 @@ void MAPS::hitana() { // wenn erste Reihe des Clusters gefüllt werden soll // und das Seed Pixel am Rand - in erster oder zweiten Reihe des Matrix liegt if ( (row==0) && (A<2) ) { - cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; +// cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; CLUSTER[(row*5)+column]=DUMMY; } // wenn zweite Reihe des Clusters gefüllt werden soll // und das Seed Pixel am Rand - in erster Reihe des Matrix liegt else if ( (row==1) && (A<1) ) { - cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; +// cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; CLUSTER[(row*5)+column]=DUMMY; } // wenn vierte Reihe des Clusters gefüllt werden soll // und das Seed Pixel am Rand - in letzter Reihe des Matrix liegt else if ( (row==3) && (A>= (fRows-1)) ) { - cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; +// cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; CLUSTER[(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) && (A>= (fRows-2)) ) { - cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; +// cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; CLUSTER[(row*5)+column]=DUMMY; } else @@ -1284,25 +1358,25 @@ void MAPS::hitana() { // 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) && (B<2) ) { - cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; +// cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; CLUSTER[(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) && (B==0) ) { - cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; +// cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; CLUSTER[(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) && (B==fColumns-1)) { - cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; +// cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; CLUSTER[(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) && (B>fColumns-3)) { - cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; +// cout << "WARNING: Clusters are filled with dummy values, please check MAPS.c hitana() cluster fill procedure." << endl; CLUSTER[(row*5)+column]=DUMMY; } else { @@ -1338,7 +1412,7 @@ void MAPS::hitana() { CHANCE=0; //NOTE: potential error source break; } - else { + else { CHANCE=100; if(i%5==3) { i+=2; @@ -1349,19 +1423,12 @@ void MAPS::hitana() { //Begin: loop evaluate true seeds: if(CHANCE==100) { - //Fill hit TTree: - fFrameInfo.pixel[fHits] = Hitlist[hit]; - for(int clupos=0; clupos<25; clupos++) - { - fFrameInfo.p [clupos][fHits] = CLUSTER[clupos]; - } - - fHits++; - - if(fSave) { - hint1->Fill( Hitlist[hit]%fColumns, (int)(Hitlist[hit]/fColumns) ); - } - + if (B < 2 || B > fColumns-3 || A < 2 || A > fRows-3) + 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++) @@ -1377,12 +1444,33 @@ void MAPS::hitana() { else if ( (column==3) && (B==fColumns-1)) { } else if ( (column==4) && (B>fColumns-3)) { } else { - fHittedPixel[Hitlist[hit]+(row-2)*fColumns+(column-2)] = 1; + fHittedPixel[Hitlist[hit]+(row-2)*fColumns+(column-2)] = bordercluster?-1:1; } } } } - fHittedPixel[Hitlist[hit]] = 2; +// cout<<"Hitted pixel discriminator matrix:"<(fHittedPixel, fPixels, fColumns, 1, 1); + if (bordercluster) + fHittedPixel[Hitlist[hit]] = -2; + else + { + fHittedPixel[Hitlist[hit]] = 2; + + //Fill hit TTree: + fFrameInfo.pixel[fHits] = Hitlist[hit]; + for(int clupos=0; clupos<25; clupos++) + { + fFrameInfo.p [clupos][fHits] = CLUSTER[clupos]; + } + + fHits++; + + if(fSave) { + hint1->Fill( Hitlist[hit]%fColumns, (int)(Hitlist[hit]/fColumns) ); + } + } + } } //End: loop evaluate true seeds: @@ -1392,53 +1480,8 @@ void MAPS::hitana() { { for(Int_t i=0; i0) + if (fHittedPixel[i]!=0) { - if (fHittedPixel[i]>1) - { -// cout << "Frame: " << fFrameNumber << " Hit at pixel: " << i << " with ADC: " << fFrameInfo.p [(2*5)+2][fHits-1]+fPedestals[Hitlist[0]] << endl; -// for(Int_t row=0; row<5; row++) -// { -// for(Int_t column=0; column<5; column++) -// { -// std::cout.width(10); -// std::cout << std::fixed; -// std::cout << std::left << std::setprecision(2) << fFrameInfo.p [(row*5)+column][fHits-1]+fPedestals[Hitlist[0]+(row-2)*fColumns+(column-2)]; -// } -// cout << endl; -// } -// cout << endl; -// cout << "Noise : "<< endl; -// for(Int_t row=0; row<5; row++) -// { -// for(Int_t column=0; column<5; column++) -// { -// std::cout.width(10); -// std::cout << std::fixed; -// std::cout << std::left << std::setprecision(2) << fNoise[Hitlist[0]+(row-2)*fColumns+(column-2)]; -// } -// cout << endl; -// } -// cout << endl; -// cout << "Pedestial : "<< endl; -// for(Int_t row=0; row<5; row++) -// { -// for(Int_t column=0; column<5; column++) -// { -// std::cout.width(10); -// std::cout << std::fixed; -// std::cout << std::left << std::setprecision(2) << fPedestals[Hitlist[0]+(row-2)*fColumns+(column-2)]; -// } -// cout << endl; -// } -// cout << endl; -// -// cout << "Noise of my pixel #0 : "<< endl; -// std::cout.width(10); -// std::cout << std::fixed; -// std::cout << std::left << std::setprecision(5) << fNoise[0]; -// cout << endl; - } fdiscriminatedhitmatrix->SetBinContent(i%fColumns, (int)(i/fColumns), fHittedPixel[i]); fADCHitmatrix->SetBinContent(i%fColumns, (int)(i/fColumns), fCdsmatrix[i]); } @@ -1479,8 +1522,7 @@ void MAPS::filterCommonMode() { Arr[column] = fCdsmatrix[ row*fColumns+column ] - fPedestals[ row*fColumns+column ]; // Bugfix: Considered pedestals in common mode calculation } - CommonModeInRow = TMath::Mean(fColumns, Arr); -// CommonModeInRow = TMath::Median(fColumns, Arr); + CommonModeInRow = TMath::Median(fColumns, Arr); // Don't use mean, consider hitted pixel ! if (CommonModeInRow < 0) plus++; else @@ -1491,7 +1533,8 @@ void MAPS::filterCommonMode() { warning = true; } if(warning) { cout<<"\rFrame: "< Common Mode suspiciously high! "<(fPedestals, fPixels, fColumns, 2); + debugStream<>(fCdsmatrix, fPixels, fColumns, 1, 20); + debugStream<>(Arr, fColumns, fColumns, 1, 20); } for(int column=0; column(fCdsmatrix, fPixels, fColumns, 1, 19); + + cout<<"Hitted pixel discriminator matrix:"<(fHittedPixel, fPixels, fColumns, 1, 1); + TCanvas* cm1; cm1 = new TCanvas("cm1","Matrix",50,100,1200,800); cm1->Divide(2,3); diff --git a/MABS_run_analyzer/MAPS.h b/MABS_run_analyzer/MAPS.h index 23625b4..77070bf 100644 --- a/MABS_run_analyzer/MAPS.h +++ b/MABS_run_analyzer/MAPS.h @@ -118,6 +118,7 @@ private: Float_t* fPedestals; std::queue* noiselastframes; + std::queue* CDSlastframes; Int_t fHits; @@ -140,7 +141,7 @@ private: /// system of the run, "PXI" or "USB", set and passed initMapsRun() to in the constructor TString fSystem; /// default number of frames to use for noise calculation - const static int numberofframesfornoise = 100; + const static int numberofframesfornoise =100; /// full path in LINUX style to the .root file, all hit information is stored there @@ -266,13 +267,19 @@ private: arraytype SumOverArray(const arraytype* (a), Int_t n); template - arraytype debugStream(const arraytype* (a), Int_t n=512, Int_t columns=8, Int_t precision=2); + arraytype 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 */ void initHistograms(); + const TString colorwhite = "\033[1;29m"; + const TString colorred = "\033[1;31m"; + const TString coloryellow = "\033[1;33m"; + const TString colorreset = "\033[0m"; + /// resets color and prints a new line + const TString endlr = "\033[0m\n"; public: /** @@ -304,6 +311,16 @@ public: */ Bool_t initOldRootFile(); + + /** @brief Refill the noise branch every #numberofframesfornoise frames or just use the initial initialization + * + * turning this on will slow down the analysis, but allows time dependant noise investigation + */ + Bool_t RefillNoiseBranch = false; + + /** @brief indicates that the "noise" branch was filled with new values, needed for noise analysis */ + Bool_t RefilledNoiseBranch = false; + /** * @brief Reads in RAW data, offsets to specific frame, sets fF0matrix, fF1matrix and fCdsmatrix * diff --git a/MABS_run_analyzer/Run.c b/MABS_run_analyzer/Run.c index e172746..c331da2 100644 --- a/MABS_run_analyzer/Run.c +++ b/MABS_run_analyzer/Run.c @@ -18,9 +18,8 @@ Run::Run( void ) Run::~Run( void) { - delete processed; db->Close(); - rootfile->Close(); + cout << "gelöscht Run" << endl; } Run::Run(Int_t runnumber) @@ -31,8 +30,8 @@ Run::Run(Int_t runnumber) //db = TSQLServer::Connect("mysql://jspc29.x-matter.uni-frankfurt.de","radhard","mimosa88"); try { - string selectquery = "select System,TempCooling,COALESCE(TempChipStart,-10000) as TempChipStart,COALESCE(TempChipEnd,-10000) as TempChipEnd,ChipNum,RadiationSource,Matrix,Clock,StorePath,ChipGen,COALESCE(VetoPeak,-1) as VetoPeak from radhard.labbook WHERE runnumber=" + numberToString<>(labbook.runnumber); - res = db->Query(prepareSQLStatement(selectquery)); + string selectquery = prepareSQLStatement("select `System`, `TempCooling`, COALESCE(`TempChipStart`,-10000) as `TempChipStart`, COALESCE(`TempChipEnd`,-10000) as `TempChipEnd`, `ChipNum`, `RadiationSource`, `Matrix`, `Clock`, `StorePath`, `ChipGen`,COALESCE(`VetoPeak`,-1) as `VetoPeak` from `radhard`.`labbook` WHERE `runnumber`=" + numberToString<>(labbook.runnumber)); + res = db->Query(selectquery.c_str()); nrows = res->GetRowCount(); if (nrows > 0) { @@ -48,23 +47,33 @@ Run::Run(Int_t runnumber) // replace windows drive notation with linux style if (labbook.storepath.Length() > 0) { - storepathLinux = labbook.storepath; - storepathLinux = storepathLinux.ReplaceAll("H:","/jspc53_H"); - storepathLinux = storepathLinux.ReplaceAll("h:","/jspc53_H"); - storepathLinux = storepathLinux.ReplaceAll("U:","/jspc53_U"); - storepathLinux = storepathLinux.ReplaceAll("u:","/jspc53_U"); - storepathLinux = storepathLinux.ReplaceAll("F:","/jspc12_F"); - storepathLinux = storepathLinux.ReplaceAll("f:","/jspc12_F"); - storepathLinux = storepathLinux.ReplaceAll("O:","/d/garlic"); - storepathLinux = storepathLinux.ReplaceAll("o:","/d/garlic"); - storepathLinux = storepathLinux.ReplaceAll("\\","/"); + storepathRAWLinux = labbook.storepath; +// storepathRAWLinux = storepathRAWLinux.ReplaceAll("H:","/jspc53_H"); +// storepathRAWLinux = storepathRAWLinux.ReplaceAll("h:","/jspc53_H"); + storepathRAWLinux = storepathRAWLinux.ReplaceAll("H:","/d/garlic"); + storepathRAWLinux = storepathRAWLinux.ReplaceAll("h:","/d/garlic"); + storepathRAWLinux = storepathRAWLinux.ReplaceAll("U:","/jspc53_U"); + storepathRAWLinux = storepathRAWLinux.ReplaceAll("u:","/jspc53_U"); + storepathRAWLinux = storepathRAWLinux.ReplaceAll("F:","/jspc12_F"); + storepathRAWLinux = storepathRAWLinux.ReplaceAll("f:","/jspc12_F"); + storepathRAWLinux = storepathRAWLinux.ReplaceAll("O:","/d/garlic"); + storepathRAWLinux = storepathRAWLinux.ReplaceAll("o:","/d/garlic"); + storepathRAWLinux = storepathRAWLinux.ReplaceAll("S:","/d/garlic"); + storepathRAWLinux = storepathRAWLinux.ReplaceAll("s:","/d/garlic"); +// storepathRAWLinux = storepathRAWLinux.ReplaceAll("O:","/home/blinnik/garlic/local"); +// storepathRAWLinux = storepathRAWLinux.ReplaceAll("o:","/home/blinnik/garlic/local"); +// storepathRAWLinux = storepathRAWLinux.ReplaceAll("S:","/home/blinnik/garlic/local"); +// storepathRAWLinux = storepathRAWLinux.ReplaceAll("s:","/home/blinnik/garlic/local"); + storepathRAWLinux = storepathRAWLinux.ReplaceAll("\\","/"); } + savepathresults = storepathRAWLinux; labbook.chipGen = (rowsql->GetField(9) != NULL)?std::string(rowsql->GetField(9)):""; labbook.posVetoDB = (rowsql->GetField(10) != NULL)?atoi(rowsql->GetField(10)):-1; delete res; if (labbook.chip.Length() > 0 && labbook.chipGen.Length() > 0) // versuche infos zum Chip aus der ChipDatenbank zu bekommen { - res = db->Query(prepareSQLStatement("select epi_thickness, resistivity, `ChipRadiation Ion`, `ChipRadiation NonIon` from radhard.chips WHERE no='" + numberToString<>(labbook.chip) + "' AND chipgen='" + numberToString<>(labbook.chipGen) + "'")); + selectquery=prepareSQLStatement("select `epi_thickness`, `resistivity`, `ChipRadiation Ion`, `ChipRadiation NonIon` from `radhard`.`chips` WHERE `no`='" + numberToString<>(labbook.chip) + "' AND `chipgen`='" + numberToString<>(labbook.chipGen) + "'"); + res = db->Query(selectquery.c_str()); nrows = res->GetRowCount(); if (nrows > 0) { @@ -80,8 +89,14 @@ Run::Run(Int_t runnumber) { getVetoPeakPositionFromFe55Run(); } + cout << colorwhite << "setSystemSpecificParameters():" << endlr; setSystemSpecificParameters(); + cout << colorwhite << "generateReadableRunCode():" << endlr; generateReadableRunCode(); + cout << colorwhite << "initRootParameters():" << endlr; + initRootParameters(); + cout << colorwhite << "initHistograms():" << endlr; + initHistograms(); runexistsinDB = 1; } else @@ -93,13 +108,34 @@ Run::Run(Int_t runnumber) { cout << "\033[1;31mError while reading laboratory book (run number " << numberToString<>(labbook.runnumber) << ")!\033[0m\n"; } +} + +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; +} + +Bool_t Run::setResultsPath(TString path) +{ + cout << "Changing save directory to '" << path << "'" << endl; + if (system("mkdir '"+ path + "' -p")) + { + cout << "\033[1;31mError changing directory, save path is: '\033[1;37m" << savepathresults << "\033[1;31m'\033[0m" << endl; + return 1; + } + if (path.EndsWith("/")) + path.Remove(path.Length()-1); + savepathresults = path; + return 0; } Bool_t Run::analyzeRun(Bool_t force) { if (runexistsinDB) { + cout << colorwhite << "instatiate MAPS(this)" << endlr; processed = new MAPS(this); if (!runAllreadyAnalyzed() || force) { @@ -114,65 +150,56 @@ Bool_t Run::analyzeRun(Bool_t force) for(int i=0; igetFrame(i); - processed->filterCommonMode(); +// processed->filterCommonMode(); processed->hitana(); - processed->regetDynNoise(); + if (dynamicalNoise) + { + processed->RefillNoiseBranch=true; + processed->regetDynNoise(); + if (processed->RefillNoiseBranch && processed->RefilledNoiseBranch) + { + binNoise(); + plotNoise(); + processed->RefilledNoiseBranch = false; + } + } 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(storepathLinux + "/rootfilecreated",std::ios::out); + fstream* fout = new fstream(storepathRAWLinux + "/rootfilecreated",std::ios::out); *fout << "" << endl; fout->close(); } else { cout << "\033[1;33mSkipped analysis of run " << labbook.runnumber << ", I think the root file for this run allready exists.\033[0m" << endl; - - processed->initOldRootFile(); -// TString rootfilepath = storepathLinux + Form ("/RUN_%i_i.root", labbook.runnumber); -// -// -// rootfile = new TFile(rootfilepath); -// if (!(rootfile->IsZombie())) -// { -// runexistsAsRootFile = 1; -// } -// else -// { -// // if we cannot open the file, print an error messageForm("hist%i",nHistNtuple) and return immediatly -// cout << "\033[1;31mError: cannot open " << rootfilepath << "\033[0m\n"; -// exit(0); -// } + cout << colorwhite << "initOldRootFile():" << endlr; + error = processed->initOldRootFile(); } -// binNoise(); -// binSeedSumVeto(); - - -// if (!(rootfile->IsZombie())) -// { -// runexistsAsRootFile = 1; -// initNoise(); -// } -// else -// { -// // if we cannot open the file, print an error messageForm("hist%i",nHistNtuple) and return immediatly -// cout << "\033[1;31mError: cannot open " << rootfilepath << "\033[0m\n"; -// exit(0); -// } - return true; + if (!error) + { + cout << colorwhite << "binNoise():" << endlr; + binNoise(); + cout << colorwhite << "binSeedSumVeto():" << endlr; + binSeedSumVeto(); + cout << colorwhite << "delete processed:" << endlr; + delete processed; + return true; + } + delete processed; } return false; } - Bool_t Run::analyzeFrame(Int_t frame) { if (runexistsinDB) { processed = new MAPS(this); + processed->initNewRootFile(); int entries = processed->GetNumberFrames(); if (frame < entries) { @@ -181,15 +208,16 @@ Bool_t Run::analyzeFrame(Int_t frame) processed->filterCommonMode(); processed->hitana(); processed->plotFrame(frame); - binNoise(); - return true; + delete processed; + return 0; } else { cout << "\033[1;33mFrame number too big, max frame: " << entries << "\033[0m" << endl; } } - return false; + delete processed; + return 1; } // generateSeedSpectrum() @@ -200,7 +228,7 @@ Bool_t Run::analyzeFrame(Int_t frame) Bool_t Run::runAllreadyAnalyzed() { - return checkFileExists(storepathLinux + "/rootfilecreated"); + return checkFileExists(storepathRAWLinux + "/rootfilecreated"); } Bool_t Run::checkFileExists(TString file) @@ -284,34 +312,32 @@ Float_t Run::stringToNumber ( string Text ) void Run::setSystemSpecificParameters() -{ - // if USB system, threshold = 5, if PXI = 65 - vetothreshold=labbook.system=="USB"?5:75; - setSystemInfo(); +{ + if (labbook.system.EqualTo("USB")) // && labbook.chipGen.EqualTo("34") ) + systemparamcur = systemparamUSB; + else if (labbook.system.EqualTo("PXI")) // && labbook.chipGen.EqualTo("34") ) + systemparamcur = systemparamPXI; + + setSensorInSystemParam(); } -void Run::setSystemInfo() -{ - sensorinfostruct sensorinfoMi34USB; - sensorinfoMi34USB.columns=64; - sensorinfoMi34USB.rows=8; - sensorinfostruct sensorinfoMi34PXI; - sensorinfoMi34PXI.columns=64; - sensorinfoMi34PXI.rows=16; - +void Run::setSensorInSystemParam() +{ + sensorinfostruct sensorinfoMi34USB( 8, 64 ); + sensorinfostruct sensorinfoMi34PXI( 16, /* rows */ 64 /* columns */ ); if (labbook.system.EqualTo("USB")) // && labbook.chipGen.EqualTo("34") ) - sensorinfocurrent=sensorinfoMi34USB; + sensorinfocur=sensorinfoMi34USB; else if (labbook.system.EqualTo("PXI")) // && labbook.chipGen.EqualTo("34") ) - sensorinfocurrent=sensorinfoMi34PXI; + sensorinfocur=sensorinfoMi34PXI; } template -char const* Run::prepareSQLStatement( T statement ) +string Run::prepareSQLStatement( T statement ) { stringstream ss; ss << statement; cout << ss.str() << endl << endl; // for debugging - return ss.str().c_str(); + return ss.str(); } void Run::getVetoPeakPositionFromFe55Run() @@ -319,7 +345,8 @@ void Run::getVetoPeakPositionFromFe55Run() try { // query database and print results - res = db->Query(prepareSQLStatement("select COALESCE(VetoPeak,-1) as VetoPeak, runnumber from radhard.labbook WHERE ChipNum='" + numberToString<>(labbook.chip) + "' AND RadiationSource='Fe55' AND Matrix='" + numberToString<>(labbook.matrix) + "' AND TempCooling=" + to_str_w_prec(labbook.temp,3) + " AND System='" + numberToString<>(labbook.system) + "'")); + TString query=prepareSQLStatement("select COALESCE(VetoPeak,-1) as VetoPeak, runnumber from radhard.labbook WHERE ChipNum='" + numberToString<>(labbook.chip) + "' AND RadiationSource='Fe55' AND Matrix='" + numberToString<>(labbook.matrix) + "' AND TempCooling=" + to_str_w_prec(labbook.temp,3) + " AND System='" + numberToString<>(labbook.system) + "'"); + res = db->Query(query.Data()); nrows = res->GetRowCount(); if (nrows > 0) { @@ -336,7 +363,8 @@ void Run::getVetoPeakPositionFromFe55Run() } if (!(Fe55run.posVeto > 0)) // no calibration peak found { - res = db->Query(prepareSQLStatement("select COALESCE(VetoPeak,-1) as VetoPeak, runnumber, TempCooling from radhard.labbook WHERE ChipNum='" + numberToString<>(labbook.chip) + "' AND RadiationSource='Fe55' AND Matrix='" + numberToString<>(labbook.matrix) + "' AND System='" + numberToString<>(labbook.system) + "'")); + query=prepareSQLStatement("select COALESCE(VetoPeak,-1) as VetoPeak, runnumber, TempCooling from radhard.labbook WHERE ChipNum='" + numberToString<>(labbook.chip) + "' AND RadiationSource='Fe55' AND Matrix='" + numberToString<>(labbook.matrix) + "' AND System='" + numberToString<>(labbook.system) + "'"); + res = db->Query(query.Data()); nrows = res->GetRowCount(); if (nrows > 0) { @@ -392,8 +420,8 @@ void Run::updateDatabase() { { try { - sqlupdatequery = "UPDATE `labbook` SET " + sqlupdatequery + " WHERE `runnumber`=" + numberToString<>(labbook.runnumber); - Bool_t sucess = db->Exec(prepareSQLStatement(sqlupdatequery)); + sqlupdatequery = prepareSQLStatement("UPDATE `labbook` SET " + sqlupdatequery + " WHERE `runnumber`=" + numberToString<>(labbook.runnumber)); + Bool_t sucess = db->Exec(sqlupdatequery.c_str()); if (!sucess) { cout << "\033[1;31mError while writing laboratory book TO SQL database (run number " << labbook.runnumber << ")!\033[0m\n"; @@ -422,24 +450,19 @@ Bool_t Run::binNoise() 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->fDynNoiseTree->SetBranchAddress("noise", &noise, &noiseBranch); - // create histogram - if (labbook.system == "PXI") - histogram.Noise = new TH1F("Noise " + runcode, "Noise" + runcode, 150, 0, 150); - else - histogram.Noise = new TH1F("Noise " + runcode, "Noise" + runcode, 100, 0, 10); - for (Int_t cnt=0; cntfDynNoiseTree->GetEntries(); cnt++) { processed->fDynNoiseTree->GetEntry(cnt); histogram.Noise->Fill(noise); } histogram.Noise->GetQuantiles( 3, noisequantiles, probabilities); - if (labbook.system == "PXI") - for (int j=0; j<3; j++) - noisequantiles[j] /= 16.0; // TODO analyze PXI scales - +// if (labbook.system == "PXI") +// for (int j=0; j<3; j++) +// noisequantiles[j] /= 16.0; // TODO analyze PXI scales + return false; } @@ -457,7 +480,7 @@ Bool_t Run::binSeedSumVeto() /// collected charge in cluster Float_t pixelSum = 0; Float_t notSeedSum = 0; - + for (Int_t framei=0; frameifHitTree->GetEntries(); framei++) // loop over all frames { processed->fHitTree->GetEntry(framei); @@ -481,27 +504,312 @@ Bool_t Run::binSeedSumVeto() histogram.Sum->Fill(pixelSum); // veto spectrum - if (TMath::Abs(notSeedSum) < vetothreshold) + if (TMath::Abs(notSeedSum) < systemparamcur.vetothreshold && labbook.source.Contains("Fe")) histogram.Veto->Fill(processed->fFrameInfo.p[12][hiti]); // histogram with the single pixel } } } } +void Run::setPlotStyle(Int_t x){ + plotStyle = x%13; + histogram.Seed->SetLineColor(rootcolors[plotStyle]); + histogram.Sum->SetLineColor(rootcolors[plotStyle]); + histogram.Veto->SetLineColor(rootcolors[plotStyle]); + histogram.Noise->SetLineColor(rootcolors[plotStyle]); +} + Bool_t Run::plotNoise() { - return plot1DHistogram(histogram.Noise); + if (!error) + { + TString legendEntry = Form("Noise: %.2f + %.2f - %.2f",noisequantiles[1],noisequantiles[1] - noisequantiles[0],noisequantiles[2] - noisequantiles[1] ); + TCanvas* canvas = plot1DHistogram(histogram.Noise, "", legendEntry); + return 0; + } + return 0; } -Bool_t Run::plot1DHistogram(TH1F* histogram, TString titlestr, TString legendstr) +Bool_t Run::plotSeed() { - TCanvas* canvas = new TCanvas(); - canvas->SetTitle(runcode + " " + titlestr); - histogram->Draw(); - TLegend* leg = new TLegend(0.5,0.5,0.89,0.89);//(0.6,0.7,0.89,0.89); - leg->SetFillStyle(0); - leg->AddEntry((TObject*) 0, runcode + " " + legendstr, ""); + if (!error) + { + plot1DHistogram(histogram.Seed); + return 0; + } +} + + +Bool_t Run::plotSum() +{ + if (!error) + { + plot1DHistogram(histogram.Sum); + return 0; + } + return 0; +} + +Bool_t Run::plotVeto() +{ + if (!error) + { + plot1DHistogram(histogram.Veto); + return 0; + } + return 0; +} + +Bool_t Run::plotAllHistograms() +{ + if (!error) + { + TString canvastitle = Form("%s", runcode.Data()); + TCanvas* canvas = new TCanvas(canvastitle, canvastitle, 1200, 800); + canvas->Divide(2,2); + canvas->cd(1); + histogram.Seed->Draw(); + canvas->cd(2); + histogram.Sum->Draw(); + canvas->cd(3); + histogram.Veto->Draw(); + canvas->cd(4); + histogram.Noise->Draw(); + TString legendEntry = Form("Noise: %.2f + %.2f - %.2f",noisequantiles[1],noisequantiles[1] - noisequantiles[0],noisequantiles[2] - noisequantiles[1] ); + TLegend* leg = new TLegend(0.5,0.8,0.89,0.89);//(0.6,0.7,0.89,0.89); + leg->AddEntry(histogram.Veto, legendEntry, "l"); + leg->SetTextSize(0.03); + leg->Draw(); + return 0; + } +} + +TCanvas* Run::plot1DHistogram(TH1F* onehistogram, TString titlestr, TString legendstr) +{ + if (titlestr.Length() < 1) + titlestr = Form("%s",onehistogram->GetTitle()); + TString canvastitle = Form("%s %s", onehistogram->GetName(), runcode.Data()); + TCanvas* canvas = new TCanvas(canvastitle, canvastitle, 700, 500); + onehistogram->SetTitle(titlestr); + onehistogram->Draw(); + TLegend* leg = new TLegend(0.5,0.8,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->Draw(); + + canvas -> SaveAs( savepathresults + "/" + runcode + " " + onehistogram->GetName() + ".eps"); + + TImage *img = TImage::Create(); + img->FromPad(canvas); + img->WriteImage(savepathresults + "/" + runcode + " " + onehistogram->GetName() + ".png"); + + TFile *f = new TFile(savepathresults + "/" + runcode + " " + onehistogram->GetName() + ".root","UPDATE"); + f->WriteTObject(canvas); + f->WriteTObject(img); + + gStyle->SetPaperSize(10.,10.); + canvas->Print(savepathresults + "/" + runcode + " " + onehistogram->GetName() + ".tex"); + + writeHistogramToFile(onehistogram); + return canvas; +} + +Bool_t Run::writeHistogramToFile(TH1F* onehistogram) +{ + system("mkdir "+ savepathresults + " -p"); + 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"); + *fout << header << endl; + + TString outline; + for(Int_t bin=0;binGetNbinsX();bin++) + { + outline=Form("%.3f\t%.3f", onehistogram->GetBinCenter(bin), onehistogram->GetBinContent(bin)); + *fout<close(); + + fout = new fstream(filename,ios::out); + return 0; +} + +Bool_t Run::writeAllHistogramsToFile() +{ + system("mkdir "+ savepathresults + " -p"); + TString filename= savepathresults + "/" + runcode + " histograms.dat"; + fstream* fout = new fstream(filename,ios::out); + + TString header = Form("#bin [ADU]\tSeed\tSum\t\Veto\tbin noise\tnoise\tbin [e]\tbin [Veto%]\tSeed"); + *fout << header << endl; + + TString outline; + for(Int_t bin=0;binGetNbinsX();bin++) + { + outline=Form("%.1f\t%.1f\t%.1f\t%.1f", histogram.Seed->GetBinCenter(bin), histogram.Seed->GetBinContent(bin), histogram.Sum->GetBinContent(bin), histogram.Veto->GetBinContent(bin)); + if (bin < systemparamcur.nbinsnoise) + outline+=Form("\t%.1f\t%.1f", histogram.Noise->GetBinCenter(bin),histogram.Noise->GetBinContent(bin)); + *fout<close(); return 0; -} \ No newline at end of file +} + +void Run::MakeGnuplotFile() +{ + system("mkdir "+ savepathresults + " -p"); + TString filename= savepathresults + "/" + runcode+"_gnuplot.plt"; + fstream* fout = new fstream(filename,ios::out); + + *fout << "reset" << endl; + *fout << "cd \"" + savepathresults + "\"" << endl; + + *fout << "datafile = \"" + runcode.ReplaceAll("_", "\\_") + " histograms.dat" + "\"" << endl; + *fout << "outputfile = \"" + runcode.ReplaceAll("_", "\\_") + " plot.png" + "\"" << endl; + *fout << "title = \"\"" << endl; + + *fout << "#set terminal postscript eps enhanced color \"Helvetica\" 18" << endl; + *fout << "#set terminal epslatex color dl 4" << endl; + *fout << "set terminal png giant size 1024,768 font \"Arial\" 18" << endl; + + *fout << "#set style line 1 pt 6 ps 1 lw 2 lt 1 lc rgb \"#0052A3\"" << endl; + *fout << "set style line 1 pt 7 ps 2 lw 2 lt 1 lc rgb \"#000000\"" << endl; + *fout << "set style line 2 pt 6 ps 1 lw 2 lt 2 lc rgb \"#ff3e3e\"" << endl; + *fout << "#set style line 3 pt 6 ps 1 lw 2 lt 1 lc rgb \"#E6B800\"" << endl; + *fout << "set style line 3 pt 7 ps 2 lw 2 lt 1 lc rgb \"#0000ff\"" << endl; + *fout << "set style line 4 pt 6 ps 1 lw 2 lt 2 lc rgb \"#00ffff\"" << endl; + *fout << "set style line 5 pt 7 ps 2 lw 2 lt 1 lc rgb \"#E65C00\"" << endl; + *fout << "set style line 6 pt 6 ps 1 lw 2 lt 2 lc rgb \"#ffe553\"" << endl; + *fout << "set style line 7 pt 6 ps 1 lw 2 lt 1 lc rgb \"#B20000\"" << endl; + *fout << "set style line 8 pt 6 ps 1 lw 2 lt 2 lc rgb \"#3eff4c\"" << endl; + *fout << "set style line 9 pt 6 ps 1 lw 2 lt 1 lc rgb \"#660099\"" << endl; + *fout << "set style line 10 pt 6 ps 1 lw 2 lt 2 lc rgb \"#c049fc\"" << endl; + *fout << "set style line 11 pt 6 ps 1 lw 2 lt 1 lc rgb \"#ff0000\"" << endl; + *fout << "set style line 12 pt 6 ps 1 lw 2 lt 2 lc rgb \"#ffaaaa\"" << endl; + + *fout << "set style line 88 pt 6 ps 1 lw 2 lt 1 lc rgb \"#000000\"" << endl; + *fout << "set style line 85 pt 6 ps 1 lw 12 lt 1 lc rgb \"#aaaaaa\"" << endl; + *fout << "set style line 86 pt 6 ps 1 lw 0.5 lt 4 lc rgb \"#aaaaaa\" " << endl; + + *fout << "set style line 90 pt 6 ps 1 lw 20 lt 1 lc rgb \"#ffcccc\"" << endl; + *fout << "set style line 91 pt 6 ps 1 lw 20 lt 1 lc rgb \"#ccccff\"" << endl; + + *fout << "set style fill transparent solid 0.5" << endl; + *fout << "#set style fill solid 0.5" << endl; + + *fout << "#set key top right spacing 1.1 height +0.5 width -1 box lt -1 lw 0.4" << endl; + *fout << "set key top right spacing 1.1 height +0.5 width -2" << endl; + *fout << "#set key top left Right reverse samplen 2" << endl; + *fout << "set ytics font \"Arial, 12\"" << endl; + *fout << "set y2tics font \"Arial, 12\"" << endl; + *fout << "set ylabel font \"Arial, 12\"" << endl; + *fout << "set y2label font \"Arial, 12\"" << endl; + *fout << "set grid front" << endl; + *fout << "scale = 1" << endl; + + *fout << "############ All parameters ##################" << endl; + + *fout << "titlename = sprintf(\"%s\",title)" << endl; + *fout << "plotdata = sprintf(\"%s\",datafile)" << endl; + *fout << "outputname = sprintf(\"%s.png\",outputfile)" << endl; + *fout << "set output outputname" << endl; + + *fout << "set multiplot" << endl; + *fout << "set size 1.062,1" << endl; + *fout << "set origin 0,0" << endl; + *fout << "set grid xtics" << endl; + *fout << "set title titlename font \"Arial, 22\"" << endl; + *fout << "set xlabel \"Non ionizing radiation [10^13 neq/cm^2]\"" << endl; + *fout << "set xrange [0:*]" << endl; + *fout << "set yrange [*:100]" << endl; + *fout << "set ylabel \"CCE [%]\"" << endl; + *fout << "set y2label \"Gain [ADU/e]\"" << endl; + *fout << "set y2label offset -15,3" << endl; + *fout << "set y2label textcolor rgb \"blue\"" << endl; + *fout << "set y2tics auto" << endl; + *fout << "set y2tics offset -5,0" << endl; + *fout << "set y2tics textcolor rgb \"blue\" norotate" << endl; + *fout << "set y2range [*:*]" << endl; + *fout << "plot plotdata u ($7):($13) w lp ls 1 t \"CCE\",\\" << endl; + *fout << "plotdata u ($7):($14) w lp ls 3 axes x1y2 t \"Gain\"" << endl; + *fout << "set size 1,0.88" << endl; + *fout << "set origin 0.055,0.053" << endl; + *fout << "set key top right spacing 17 height +0.5 width -2" << endl; + *fout << "unset title" << endl; + *fout << "unset grid" << endl; + *fout << "unset xtics" << endl; + *fout << "set xlabel \" \"" << endl; + *fout << "set ytics \" \"" << endl; + *fout << "set ylabel \" \"" << endl; + *fout << "set y2range [*:*]" << endl; + *fout << "set y2tics auto" << endl; + *fout << "set y2tics out offset 0,0" << endl; + *fout << "set y2label \"Avg. Noise [e]\"" << endl; + *fout << "set y2label textcolor rgb \"#ff0000\"" << endl; + *fout << "set y2label offset -3,3" << endl; + *fout << "set y2tics textcolor rgb \"#ff0000\" norotate" << endl; + *fout << "plot plotdata u ($7):($14) w lp ls 5 axes x1y2 t \"Noise\"" << endl; + *fout << "unset multiplot" << endl; + *fout << "set output \"temp\"" << endl; + + fout->close(); + system("gnuplot "+ filename); + //system("okular "+ ergebnisfile + ".png &"); +} + +void Run::initHistograms() +{ + TString prefix = "Seed"; + TString humanreadablestr = Form("%s, %s spectrum, Mi%s, chip %s, %s, T=%.1f",prefix.Data(), labbook.source.Data(), labbook.chipGen.Data(), labbook.chip.Data(), labbook.matrix.Data(), labbook.temp); + histogram.Seed=new TH1F(prefix.Data(), humanreadablestr.Data(), systemparamcur.nbins, 0, systemparamcur.maxbin); + histogram.Seed->SetLineStyle(rootlinestyle[plotStyle]); + histogram.Seed->SetLineColor(rootcolors[plotStyle]); + histogram.Seed->SetStats(kTRUE); + histogram.Seed->SetStats(111111111); + histogram.Seed->SetLineWidth(3); + + prefix = "Sum"; + humanreadablestr = Form("%s, %s spectrum, Mi%s, chip %s, %s, T=%.1f",prefix.Data(), labbook.source.Data(), labbook.chipGen.Data(), labbook.chip.Data(), labbook.matrix.Data(), labbook.temp); + histogram.Sum=new TH1F(prefix.Data(), humanreadablestr.Data(), systemparamcur.nbins, 0, systemparamcur.maxbin); + histogram.Sum->SetLineStyle(rootlinestyle[plotStyle]); + histogram.Sum->SetLineColor(rootcolors[plotStyle]); + histogram.Sum->SetLineWidth(3); + + /// TH1F histogram with veto spectrum, used to identify the veto peak + prefix = "Veto"; + humanreadablestr = Form("%s, %s spectrum, Mi%s, chip %s, %s, T=%.1f",prefix.Data(), labbook.source.Data(), labbook.chipGen.Data(), labbook.chip.Data(), labbook.matrix.Data(), labbook.temp); + histogram.Veto=new TH1F(prefix.Data(), humanreadablestr.Data(), systemparamcur.nbins, 0, systemparamcur.maxbin); + histogram.Veto->SetLineStyle(rootlinestyle[plotStyle]); + histogram.Veto->SetLineColor(rootcolors[plotStyle]); + histogram.Veto->SetLineWidth(3); + histogram.Veto->SetAxisRange(300,550,"X"); + + prefix = "Noise"; + humanreadablestr = Form("%s, %s spectrum, Mi%s, chip %s, %s, T=%.1f",prefix.Data(), labbook.source.Data(), labbook.chipGen.Data(), labbook.chip.Data(), labbook.matrix.Data(), labbook.temp); + histogram.Noise=new TH1F(prefix.Data(), humanreadablestr.Data(), systemparamcur.nbinsnoise, 0, systemparamcur.maxbinnoise); + histogram.Noise->SetLineStyle(rootlinestyle[plotStyle]); + histogram.Noise->SetLineColor(rootcolors[plotStyle]); + histogram.Noise->SetLineWidth(3); +} + + +// void Run::initHistogram(TH1F* curhistogram, TString prefix) +// { +// TString humanreadablestr = Form("%s, %s spectrum, Mi%s, chip %s, %s, T=%.1f",prefix.Data(), labbook.source.Data(), labbook.chipGen.Data(), labbook.chip.Data(), labbook.matrix.Data(), labbook.temp); +// curhistogram=new TH1F(prefix.Data(), humanreadablestr.Data(), systemparamcur.maxbin, 0, systemparamcur.nbins); +// curhistogram->SetLineStyle(rootlinestyle[plotStyle]); +// curhistogram->SetLineColor(rootcolors[plotStyle]); +// curhistogram->SetLineWidth(3); +// } + +void Run::initRootParameters() +{ rootcolors = new Int_t[13]{1, 2, 4, 6, 8, 13, 46, 28, 32, 33, 12, 20, 40}; + rootlinestyle = new Int_t[13]{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; + +} diff --git a/MABS_run_analyzer/Run.h b/MABS_run_analyzer/Run.h index fbfee1d..a11303e 100644 --- a/MABS_run_analyzer/Run.h +++ b/MABS_run_analyzer/Run.h @@ -12,10 +12,10 @@ #include #include #include -#include "TLegend.h" - -#include "TH1F.h" -#include "TF1.h" +#include +#include +#include +#include #include "help.h" @@ -48,6 +48,8 @@ private: /// number of rows the SQL database found to a given request int nrows=0; + Int_t plotStyle=0; + /** * @brief Takes a number, returns a string, used to build SQL queries */ @@ -56,7 +58,7 @@ private: /** @brief formats and checks SQL statements */ template - char const* prepareSQLStatement( T ); + string prepareSQLStatement( T ); /** * @brief tries to get veto peak position from other, but close runs @@ -78,11 +80,11 @@ private: struct sensorinfostruct { /// number of rows read out in the system the run was taken with (USB or PXI) - Int_t rows=0; + Int_t rows; /// number of columns read out in the system the run was taken with (USB or PXI) - Int_t columns=0; - // sensorinfo():rows(0),columns(0){}; - // sensorinfo( int d, int m ) : rows( d ), columns( m ) {} + Int_t columns; + sensorinfostruct():rows(0),columns(0){}; + sensorinfostruct( Int_t r, Int_t c ) : rows( r ), columns( c ) {} }; /** @@ -126,6 +128,9 @@ private: /// is set to true if the RAW Data is consistent Bool_t runRAWok = 0; + /// is set to true if an error occured + Bool_t error = 0; + /** * @brief fills noise #histogram */ Bool_t binNoise(); @@ -135,12 +140,66 @@ private: /// noise quantiles: mean value, sigma in postive direction and sigma in negative direction Double_t noisequantiles[3]; - Bool_t plot1DHistogram(TH1F* histogram, TString titlestr = "", TString legendstr = ""); + TCanvas* plot1DHistogram(TH1F* onehistogram, TString titlestr = "", TString legendstr = ""); /** * @brief Checks if a file exists */ Bool_t checkFileExists(TString); + + + /** @brief information about the submatrices of the sensor + */ + struct pixelinfo + { + Int_t ndiods; + Float_t diodesize; + /// TODO: add more + }; + + /** @brief system specific constants + */ + struct systemparam + { + Int_t maxbin; + Int_t nbins; + Int_t vetothreshold; + Int_t maxbinnoise; + Int_t nbinsnoise; + }; + systemparam systemparamUSB { + 800, // maxbin; + 800/4,// nbins; + 5, //vetothreshold + 10, + 100 + }; + systemparam systemparamPXI { + 800*16, // maxbin + 800,// nbins; + 75, //vetothreshold + 150, + 150 + }; + + Int_t* rootcolors; + Int_t* rootlinestyle; + + /// init all histograms, set binwidth, bin amount and names + void initHistograms(); + + // init one specific histogram +// void initHistogram(TH1F* curhistogram, TString prefix); + + void initRootParameters(); + + Bool_t dynamicalNoise = 1; + + const TString colorwhite = "\033[1;29m"; + const TString colorred = "\033[1;31m"; + const TString coloryellow = "\033[1;33m"; + const TString colorreset = "\033[0m"; + const TString endlr = "\033[0m\n"; public: /** @brief empty constructor */ @@ -150,6 +209,22 @@ public: /** @brief writes values back to the database, destroys the #MAPS object */ ~Run (void); + systemparam systemparamcur; + + void setPlotStyle(Int_t); + + /** @brief Plots all histograms from #histogram into one canvas */ + Bool_t plotAllHistograms(); + + /** @brief Writes a given histogram into a file */ + Bool_t writeHistogramToFile(TH1F* onehistogram); + + /** @brief Writes all histogram of #histogramsctruct into a file */ + Bool_t writeAllHistogramsToFile(); + + /** @brief Makes a gnuplot file to plot the histogram data created in @c writeAllHistogramsToFile() */ + void MakeGnuplotFile(); + /** * @brief analysis the RAW data * @@ -172,6 +247,17 @@ public: Bool_t runAllreadyAnalyzed(); Bool_t plotNoise(); + Bool_t plotSeed(); + Bool_t plotSum(); + Bool_t plotVeto(); + + /** @brief Turn on or off the use of dynamical noise adjustment @c MAPS::regetDynNoise() */ + Bool_t useDynamicalNoise(Bool_t); + + /** + * @brief set pathwhere images and data to save to + */ + Bool_t setResultsPath(TString path); /** * @brief provides information from the SQL database @@ -213,12 +299,15 @@ public: labbooksctruct labbook; /// path to the RAW files on LINUX systems - TString storepathLinux = ""; + TString storepathRAWLinux = ""; + + /// path to the RAW files on LINUX systems + TString savepathresults = "."; TFile* rootfile; /** @brief sets #sensorinfocurrent after run data got from db, USB or PXI */ - void setSystemInfo(); + void setSensorInSystemParam(); /** @brief sets system dependant variables after run data got from db, USB or PXI */ void setSystemSpecificParameters(); @@ -243,19 +332,6 @@ public: }; histogramstruct histogram; - /// maximal bin in - Int_t maxbinusb; - Int_t maxbinPXI; - - - /** @brief information about the submatrices of the sensor - */ - struct pixelinfo - { - Int_t ndiods; - Float_t diodesize; - /// TODO: add more - }; pixelinfo pixelinfoMi34[32]; /// analyze only half of matrix @@ -269,10 +345,8 @@ public: Float_t posSum=0; /// fitted position of the calibration peak of Fe55-beta-photons in the seed spectrum, from a run best suited to the current Float_t posVeto=0; - - /// veto threshold, if sum in cluster small then given value, consider as direct diode hit and add to veto histogram - /// is set to 5 if USB system and to 75 if PXI - Float_t vetothreshold; + + Int_t maxbin; /** @brief related Fe55 run * @@ -322,7 +396,8 @@ public: MAPS *processed; /// sensor information to use in analysis, is the system read out by USB or PXI? Number of rows differ - sensorinfostruct sensorinfocurrent; + sensorinfostruct sensorinfocur; + }; #endif \ No newline at end of file diff --git a/MABS_run_analyzer/help.h b/MABS_run_analyzer/help.h index 6ce53da..138635c 100644 --- a/MABS_run_analyzer/help.h +++ b/MABS_run_analyzer/help.h @@ -10,8 +10,9 @@ #include #include "TStopwatch.h" -#define MAXHITS 1000000 +#include +#define MAXHITS 1000000 //#################################################################### Int_t print_progress(Float_t ProgressInPercent) { @@ -108,5 +109,35 @@ struct frameInfo{ Float_t p [25][MAXHITS]; }; + +void preparecanvas() { +// TCanvas* canvas = new TCanvas("Canvas1", 700, 500); +// histogram->SetTitle(titlestr); +// histogram->SetLineStyle(rootlinestyle[plotStyle]); +// histogram->SetLineColor(rootcolors[plotStyle]); +// histogram->SetLineWidth(3); +// histogram->Draw(); +// TLegend* leg = new TLegend(0.5,0.8,0.89,0.89);//(0.6,0.7,0.89,0.89); +// leg->SetFillColor(0); +// leg->SetBorderSize(0); +// if (legendstr.Length() < 1) +// legendstr = humanreadablestr; +// leg->AddEntry(histogram, legendstr, "l"); +// +// leg->SetTextSize(0.03); +// leg->Draw(); +// +// canvas -> SaveAs( savepathresults + "/" + runcode + " " + histogram->GetName() + ".eps"); +// +// TImage *img = TImage::Create(); +// img->FromPad(canvas); +// img->WriteImage(savepathresults + "/" + runcode + " " + histogram->GetName() + ".png"); +// +// TFile *f = new TFile(savepathresults + "/" + runcode + " " + histogram->GetName() + ".root","UPDATE"); +// f->WriteTObject(canvas); +// f->WriteTObject(img); +// return 0; +} + //#################################################################### #endif \ No newline at end of file diff --git a/MABS_run_analyzer/runlist.txt b/MABS_run_analyzer/runlist.txt index f9dc87f..d43bbc4 100644 --- a/MABS_run_analyzer/runlist.txt +++ b/MABS_run_analyzer/runlist.txt @@ -1 +1 @@ -342244 +342327 \ No newline at end of file