Analysis: few bugfixes, less warnings, new veto peak search algorithm, ROOT file...

diff --git a/MABS_run_analyzer/ChargeSpektrum.c b/MABS_run_analyzer/ChargeSpektrum.c
index f397d89c77a08fb96af29a5161c888915e5589ca..51c756f124f994310d8fd2e1d623cea0801b0b89 100644 (file)
--- a/MABS_run_analyzer/ChargeSpektrum.c
+++ b/MABS_run_analyzer/ChargeSpektrum.c
@@ -27,7 +27,7 @@
 #include "CSVRow.C"
 #include <TTimeStamp.h>
 
-Int_t* ReadRunList(Int_t*);
+Int_t* ReadRunList();
 void plotAllRuns();
 
 Run** runs;
@@ -52,10 +52,10 @@ void ChargeSpektrum(Int_t runnumber = -1)
     {
         /// number of runs to be analyzed, number of lines read by @c ReadRunList() 
         numberRuns=0;
-        ReadRunList(&numberRuns);
+        ReadRunList();
         /// array with run numbers
         runList=new Int_t[numberRuns];
-        runList=ReadRunList(&numberRuns);
+        runList=ReadRunList();
     }
     runs = new Run*[numberRuns];
     
@@ -132,8 +132,7 @@ void ChargeSpektrum(Int_t runnumber = -1)
     
 }
 
-
-Int_t* ReadRunList(Int_t* numberRuns)
+Int_t* ReadRunList()
 {
     Int_t* runList=new Int_t[1000];
     
@@ -152,7 +151,7 @@ Int_t* ReadRunList(Int_t* numberRuns)
             cout << "File ended";
         }
     }
-    *numberRuns=runLauf;
+    numberRuns=runLauf;
     return runList;
 }
 

diff --git a/MABS_run_analyzer/MAPS.c b/MABS_run_analyzer/MAPS.c
index c39143ad1a867c091196e5f418dd18cffc52bf1e..5f4118a0d487652e789ec5701f7fa23ecf9d6035 100644 (file)
--- a/MABS_run_analyzer/MAPS.c
+++ b/MABS_run_analyzer/MAPS.c
@@ -19,9 +19,6 @@ MAPS::MAPS() {
 MAPS::MAPS(Run* runp) {
        run = runp;
     initMapsRun();
-       
-    UInt_t     FileEvNb;               // Event number per file
-    UInt_t     TotEvNb;                // Total events in RUN - should be less if RUN was stopped before endl
 
     // check if sepcified system is correct, if not, switch and reinitialize
     if ( switchsystem() )
@@ -160,6 +157,26 @@ MAPS::~MAPS(void) {
 
 };
 
+
+void MAPS::save() {
+    if(fSave && fOk)
+    {
+        fOutputFile->cd();
+        
+        fHitTree->Write("",TObject::kOverwrite);
+        fNoiseTree->Write("",TObject::kOverwrite);
+        fDynNoiseTree->Write("",TObject::kOverwrite);
+        
+        hint1->Write("",TObject::kOverwrite);
+        hint2->Write("",TObject::kOverwrite);
+        
+        fOutputFile->Save();
+        fOutputFile->Close();
+        
+        cout<< fRootFile <<" saved!"<<endl;
+    }
+}
+
 //####################################################################
 
 void MAPS::initMapsRun( ) {
@@ -227,8 +244,6 @@ bool MAPS::switchsystem( )
 
         NrAdcBoards = littleEndian32( RAWDATA, 4*9);    // Number of Adc boards installed in the system
 
-
-        TString System[5] = {"","","USB","","PXI"};
         if( fSystem == "" )
         {
             if( NrAdcBoards == 1)   {
@@ -721,10 +736,7 @@ bool MAPS::InitialDynNoise(Int_t startframe, Int_t frames) {
     /// number of frames cut away in third noise estimate (hopefully hit frames)
     Int_t* nframescutawaythirdnoiseestimate = new Int_t[fPixels]();   
     Int_t npixelscutawaythirdnoiseestimate=0; 
-    
-    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);
@@ -1427,8 +1439,6 @@ void MAPS::hitana() {
                     bordercluster = true;
                 else
                     bordercluster = false;
-                Int_t posmatrixrow;
-                Int_t posmatrixcol;
                 for(Int_t row=0; row<5; row++)
                 {
                     for(Int_t column=0; column<5; column++)

diff --git a/MABS_run_analyzer/MAPS.h b/MABS_run_analyzer/MAPS.h
index 77070bf6be1123b981e33cc5c4823ef4dbe20e78..0efd512781da7c9f091ac179782cf6c379c4601d 100644 (file)
--- a/MABS_run_analyzer/MAPS.h
+++ b/MABS_run_analyzer/MAPS.h
@@ -78,7 +78,7 @@ private:
     Int_t       fColumns;
     /// Variable to be able to subdivide the Matrix according to fOrderCode
     TString       fMatrix;            
-    Int_t       fEventsSum;
+    UInt_t       fEventsSum;
     Int_t       fFile;
     
     Int_t       fFrameNumber; /**< enum value 1 */
@@ -301,7 +301,15 @@ public:
      * 
      */
     Bool_t initNewRootFile();
+
     
+    /**
+     * @brief Saves the TTrees into ROOT file
+     *
+     * Saves the ROOT file
+     * 
+     */
+    void save();
     
     /**
      * @brief Reopens an allready existing ROOT file

diff --git a/MABS_run_analyzer/Run.c b/MABS_run_analyzer/Run.c
index 737fe955f484bf85265bbc245c0916cf687ef427..532dc9bd4efcdfa22de2eb63d62e97ebe55e0535 100644 (file)
--- a/MABS_run_analyzer/Run.c
+++ b/MABS_run_analyzer/Run.c
@@ -150,12 +150,14 @@ Bool_t Run::debugDBreadout()
     cout << "| CCE_in_Perc_1DB: " << std::right << colorwhite << labbook.CCE_in_Perc_1DB  <<  endlr;
     cout << "|______________________________ " << endlr;
     cout << endlr;
+    return 0;
 }
 
 Bool_t Run::useDynamicalNoise(Bool_t var)
 {
     cout<<" Dynamical Noise calc. : " << colorwhite << (var?"1":"0") << colorreset << "  <-- only used if RAW files are analyzed, force analysis to make sure" << endl;    
     dynamicalNoise = var;
+    return 0;
 }
     
 
@@ -199,22 +201,21 @@ Bool_t Run::analyzeRun(Bool_t force)
                 progress = (Int_t)(((i-start)*100)/(nframes-1)*10);
                 if (progress!=progress_tmp)  { print_progress( (((i-start)*100.)/(nframes-1)) ); progress_tmp=progress;}
             }
-            cout << processed->plus << " vs. " << processed->minus << " : " << (processed->plus*1.0)/processed->minus<< endl;
+//             cout << processed->plus << " vs. " << processed->minus << " : " << (processed->plus*1.0)/processed->minus<< endl;
 
             // print a dummy file to indicate, that a root file was created once
             fstream* fout = new fstream(storepathRAWLinux + "/rootfilecreated",std::ios::out);
             *fout << "" << endl;
             fout->close();
-            
-            // TODO
-            
+            processed->save();
         }
         else
         {
             cout << "Skipped analysis of run " << labbook.runnumber << ", I think the root file for this run allready exists." << endl;
             cout << colorwhite << "initOldRootFile():" << endlr;
-            error = processed->initOldRootFile();
         }
+        error = processed->initOldRootFile();
+        
         if (!error)
         {
             cout << colorwhite << "binNoise():" << endlr;
@@ -284,7 +285,6 @@ Bool_t Run::rescaleHistograms()
     histogramCalibrated.avgNoise = histogram.avgNoise * gain;
     histogramCalibrated.avgNoisePlus = histogram.avgNoisePlus * gain;
     histogramCalibrated.avgNoiseMinus = histogram.avgNoiseMinus * gain;
-    
     return 1;    
 }
 
@@ -381,6 +381,7 @@ Bool_t Run::generateReadableRunCode()
         else
             runcode+= Form("_m%.0f-%.0f_", devided_lower_matr_min, devided_lower_matr_max);
     }
+    return 0;
 }
 
 string Run::removeNumbers(std::string x)
@@ -561,15 +562,12 @@ Bool_t Run::binNoise()
     Float_t noise;
     TBranch* noiseBranch;
     Double_t const probabilities[] = {0.3415/2, 0.5, 1-0.3415/2}; // sigma/2 from gaus to the left and to the right //{0.17, 0.5, 1-0.17};
-    
     histogram.Noise->Reset();
-    processed->fNoiseTree->SetBranchAddress("noise", &noise, &noiseBranch);
-    
+    processed->fNoiseTree->SetBranchAddress("noise", &noise, &noiseBranch);    
     for (Int_t cnt=0; cnt<processed->fNoiseTree->GetEntries(); cnt++) {
         processed->fNoiseTree->GetEntry(cnt);
         histogram.Noise->Fill(noise);
-    }
-    
+    }    
     histogram.Noise->GetQuantiles( 3, noisequantiles, probabilities);
     histogram.avgNoise = noisequantiles[1];
     histogram.avgNoisePlus = noisequantiles[2] - noisequantiles[1];
@@ -577,20 +575,15 @@ Bool_t Run::binNoise()
 //     if (labbook.system == "PXI")
 //         for (int j=0; j<3; j++)
 //             noisequantiles[j] /= 16.0; // TODO analyze PXI scales
-        
-    return false;
+    return 0;
 }
 
 Bool_t Run::binSeedSumVeto()
 {    
-    /// framenumber
-    UInt_t framen;
-    /// number of hits in current frame
-    UInt_t hitsInframe;
-    /// pixel number of seed pixel, position on sensor
-    UInt_t seedPixel[10000];
-    /// Array of CLUSTERSIZE * CLUSTERSIZE clusters, seed pixel in the middle
-    Float_t pixelcluster[processed->clustersize*processed->clustersize][10000];
+//     /// pixel number of seed pixel, position on sensor
+//     UInt_t seedPixel[10000];
+//     /// Array of CLUSTERSIZE * CLUSTERSIZE clusters, seed pixel in the middle
+//     Float_t pixelcluster[processed->clustersize*processed->clustersize][10000];
     
     /// collected charge in cluster
     Float_t pixelSum = 0;
@@ -600,9 +593,9 @@ Bool_t Run::binSeedSumVeto()
     {
         processed->fHitTree->GetEntry(framei);
         // account only frames with less then 10 hits
-        if (processed->fFrameInfo.hits<10)
+        if (processed->fFrameInfo.hits<(unsigned int)10)
         {
-            for(Int_t hiti=0; hiti<processed->fFrameInfo.hits;hiti++)
+            for(Int_t hiti=0; (unsigned int)hiti<processed->fFrameInfo.hits;hiti++)
             {
                 // histogram with the single pixel
                 histogram.Seed->Fill(processed->fFrameInfo.p[12][hiti]);
@@ -630,6 +623,7 @@ Bool_t Run::binSeedSumVeto()
     histogram.posSeed=FitPerform(histogram.Seed);
     histogram.posSum=FitPerform(histogram.Sum, "gaus");
     gROOT->SetBatch(kFALSE);
+    return 0;
 }
 
 void Run::setPlotStyle(Int_t x){
@@ -648,7 +642,7 @@ Bool_t Run::plotNoise()
         noisecanvas = plot1DHistogram(histogram.Noise, "landau", "", legendEntry);
         return 0;
     }
-    return 0;
+    return 1;
 }
 
 Bool_t Run::plotSeed()
@@ -668,7 +662,7 @@ Bool_t Run::plotSum()
         plot1DHistogram(histogram.Sum, "gaus");
         return 0;
     }
-    return 0;
+    return 1;
 }
 
 Bool_t Run::plotVeto()
@@ -678,7 +672,7 @@ Bool_t Run::plotVeto()
         plot1DHistogram(histogram.Veto, "gaus");
         return 0;
     }
-    return 0;
+    return 1;
 }
 
 // Int_t findMaxBin(TH1F* th1fpointer) {
@@ -739,6 +733,7 @@ Bool_t Run::plotAllHistograms(histogramstruct* histogramstructpointer)
         
         return 0;
     }
+    return 1;
 }
 
 Bool_t Run::plotAllHistograms()
@@ -748,6 +743,7 @@ Bool_t Run::plotAllHistograms()
         plotAllHistograms(&histogram);
         return 0;
     }
+    return 1;
 }
 
 Bool_t Run::plotAllHistogramsCalibrated()
@@ -765,6 +761,7 @@ Bool_t Run::plotAllHistogramsCalibrated()
 Float_t Run::FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verbose)
 {
     Float_t posMax = 0;
+    Float_t posMax2 = 0;
     Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax();
     Float_t noiseborder = posMaxValHist/15;
     
@@ -780,55 +777,74 @@ Float_t Run::FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verb
             { 
                 Float_t peak1 = histogrampointer->GetMaximumBin();
                 Float_t peak2;
-                if (histogrampointer->GetBinCenter(peak1)<200)
+                if (histogrampointer->GetBinCenter(peak1)<3.3*noiseborder) // vermutlich ist veto peak nicht am höchsten
                 {
                     Float_t avg = 0;
                     for(Int_t bin=histogrampointer->GetXaxis()->FindBin(noiseborder);bin<histogrampointer->FindLastBinAbove(0);bin++)
                         avg += histogrampointer->GetBinContent(bin);
                     avg /= histogrampointer->FindLastBinAbove(0) - histogrampointer->GetXaxis()->FindBin(noiseborder);
-                    peak2 = histogrampointer->FindLastBinAbove(avg/4);
+                    peak2 = histogrampointer->FindLastBinAbove(avg/3);
                 }
                 else
                 {
                     peak2 = peak1;
                     peak1 = histogrampointer->GetXaxis()->FindBin(noiseborder);
-                }             
+                }            
                 histogrampointer->GetXaxis()->SetRange(peak1,peak2-1);   // look only for maxima with x greater than noiseborder, cut away noise   
                 Float_t min = histogrampointer->GetMinimumBin();
                 histogrampointer->GetXaxis()->SetRange(min,histogrampointer->FindLastBinAbove(0));   // look only for maxima with x greater than noiseborder, cut away noise
                 // DEBUG 
-                if (verbose)
-                {
-                cout << coloryellow << "peak1 " << histogrampointer->GetXaxis()->GetBinCenter(peak1) << endlr;
-                cout << coloryellow << "peak1val " << histogrampointer->GetBinContent(peak1) << endlr;
-                cout << coloryellow << "peak2 " << histogrampointer->GetXaxis()->GetBinCenter(peak2) << endlr;
-                cout << coloryellow << "peak2val " << histogrampointer->GetBinContent(peak2) << endlr;
-                cout << coloryellow << "min " << histogrampointer->GetXaxis()->GetBinCenter( min ) << endlr;
-                }
+//                 if (verbose)
+//                 {
+//                 cout << coloryellow << "peak1 " << histogrampointer->GetXaxis()->GetBinCenter(peak1) << endlr;
+//                 cout << coloryellow << "peak1val " << histogrampointer->GetBinContent(peak1) << endlr;
+//                 cout << coloryellow << "peak2 " << histogrampointer->GetXaxis()->GetBinCenter(peak2) << endlr;
+//                 cout << coloryellow << "peak2val " << histogrampointer->GetBinContent(peak2) << endlr;
+//                 cout << coloryellow << "min " << histogrampointer->GetXaxis()->GetBinCenter( min ) << endlr;
+//                 }
                 histogrampointer->Fit(fitFunc, "N,Q,W", "", histogrampointer->GetXaxis()->GetBinCenter(min), posMaxValHist);
+                posMax = histogrampointer->GetXaxis()->GetBinCenter(histogrampointer->GetMaximumBin()); // Methode 1
                 histogrampointer->GetXaxis()->SetRange(histogrampointer->GetXaxis()->FindBin(noiseborder),histogrampointer->GetXaxis()->FindBin(posMaxValHist));   // look only for maxima with x greater than noiseborder, cut away noise
             } else {
                 histogrampointer->Fit(fitFunc, "N,Q,W", "", noiseborder, posMaxValHist);
+                posMax = fitFunc->GetMaximumX(); // Methode 2
             }
             Float_t sigma = fitFunc->GetParameter(2);
-            if (sigma > 260)
+            posMax2 = fitFunc->GetMaximumX(); // Methode 2
+            Float_t peakposdifperc = abs(posMax-posMax2)/min(posMax,posMax2);
+            if (sigma > 260 || peakposdifperc>0.2)
             {
                 if (verbose)
                 {
-                    cout << "Sigma suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma  << endl;
-                    cout << colorred << "  Could not find " << histogrampointer->GetName() << " peak" << endlr;
+                    if (sigma > 260)
+                    {
+                        cout << "Sigma suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma  << endl;
+                        cout << colorred << "  Could not find " << histogrampointer->GetName() << " peak" << endlr;
+                    }
+                    if (peakposdifperc>0.2)
+                    {
+                        cout << "Maximum peak position and fit gaussian peak position doesn't match in " << histogrampointer->GetName() << " spectrum. Difference: " <<  peakposdifperc*100 <<" % " << endl;
+                        cout << colorred << "  Could not find " << histogrampointer->GetName() << " peak" << endlr;
+                    }
                 }
                 return posMax;
             }
-            else if (sigma > 40)
+            else if (sigma > 40 || peakposdifperc>0.1)
             {
                 if (verbose)
                 {
-                    cout << "Sigma suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma  << endl;
-                    cout << coloryellow << "  Please check " << histogrampointer->GetName() << " peak position: " << colorreset << colorwhite << fitFunc->GetMaximumX()  << endlr;
+                    if (sigma > 260)
+                    {
+                        cout << "Sigma or  suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma  << endl;
+                        cout << coloryellow << "  Please check " << histogrampointer->GetName() << " peak position: " << colorreset << colorwhite << posMax  << endlr;
+                    }
+                    if (peakposdifperc>0.1)
+                    {
+                        cout << "Maximum peak position and fit gaussian peak position in " << histogrampointer->GetName() << " have difference of " <<  peakposdifperc*100 <<" % " << endl;
+                        cout << coloryellow << "  Please check " << histogrampointer->GetName() << " peak position: " << colorreset << colorwhite << posMax  << endlr;
+                    }
                 }
-            }            
-            posMax = fitFunc->GetMaximumX();
+            }       
             fitFunc->DrawCopy("same");
             TString  legendEntry = TString(Form("%s",runcode.Data()));
             TLegend* leg = new TLegend(0.5,0.5,0.89,0.89);//(0.6,0.7,0.89,0.89);
@@ -942,7 +958,7 @@ Bool_t Run::writeHistogramToFile(TH1F* onehistogram)
     TString filename= savepathresults + "/" + runcode + " " + onehistogram->GetName() + " hist.dat";
     fstream* fout = new fstream(filename,ios::out);        
     
-    TString header = Form("#bin [ADU]\tCounts\tbin [e]\tbin [Veto%]\tSeed");
+    TString header = Form("#bin [ADU]\tCounts");
     *fout << header << endl;
     
     TString outline;
@@ -963,7 +979,7 @@ Bool_t Run::writeAllHistogramsToFile()
     TString filename= savepathresults + "/" + runcode + " histograms.dat";
     fstream* fout = new fstream(filename,ios::out);
     
-    TString header = Form("#bin [ADU]\tbin [e]\tSeed\tSum\t\Veto\tbin noise [ADU]\tbin noise [e]\tnoise\tbin [Veto%]\tSeed\n");
+    TString header = Form("#bin [ADU]\tbin [e]\tSeed\tSum\t\Veto\tbin noise [ADU]\tbin noise [e]\tnoise\n");
     header += Form("#posVeto, run: %.1f, DB: %.1f, Fe55 DB (%d, %.1f): %.1f\n", histogram.posVeto, labbook.posVetoDB, Fe55run.posVetorunnumber, Fe55run.temperature, Fe55run.posVeto);
     header += Form("#posSeed, run: %.1f, DB: %.1f\n", histogram.posSeed, labbook.posSeedDB);
     header += Form("#posSum, run: %.1f, DB: %.1f\n", histogram.posSum, labbook.posSumDB);