getPixelWithFakehits

diff --git a/MABS_run_analyzer/MAPS.c b/MABS_run_analyzer/MAPS.c
index 3490fde462a7fce3ffd7586332d2f8b9897c6853..8d3e6d628c5287d39f0529b9e34458f4a8d4f02a 100644 (file)
--- a/MABS_run_analyzer/MAPS.c
+++ b/MABS_run_analyzer/MAPS.c
@@ -77,36 +77,45 @@ Bool_t MAPS::initNewRootFile() {
 
 Bool_t MAPS::initOldRootFile() {
     fSave = false;
-    fOutputFile     = new TFile(fRootFile,"READ");
-    if (fOutputFile->IsZombie())
-    {
-        cout << colorred << "Error opening old ROOT file! -> Zombie"<< 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]);
-    fHitTree->SetBranchAddress ("pixelthreshold", &fFrameInfo.pixelthreshold[0]);
-    for(int i=0; i<25; i++) {
-        fHitTree->SetBranchAddress( Form("p%i",i+1)   , &fFrameInfo.p[i][0]);
+       if ( !checkConf(fPixelsData) ) { // TODO FileEvNbInConfig
+        defaultConf();
     }
-    fNoiseTree =  (TTree*) fOutputFile->Get("noise");
-    if (!fNoiseTree) // got pointer
+    //-----------------------------------------------
+    //Check and open Data Files
+    int MaxFiles = TMath::Ceil((Float_t) FileTotalEvNbInConfig/FileEvNbInConfig);
+    if( checkDataFiles(MaxFiles) )
     {
-        cout<< colorred << "Error finding noise TTree in old ROOT file!"<< endlr;
-        return 1;
-    }
-    fDynNoiseTree   = fNoiseTree;
-    initHistograms();  
-    cout<<"-----------------------"<<endl;
-    cout<<fRootFile<<" opened!"<<endl;
-    cout<<"-----------------------"<<endl;
+               fOutputFile     = new TFile(fRootFile,"READ");
+               if (fOutputFile->IsZombie())
+               {
+                       cout << colorred << "Error opening old ROOT file! -> Zombie"<< 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]);
+               fHitTree->SetBranchAddress ("pixelthreshold", &fFrameInfo.pixelthreshold[0]);
+               for(int i=0; i<25; i++) {
+                       fHitTree->SetBranchAddress( Form("p%i",i+1)   , &fFrameInfo.p[i][0]);
+               }
+               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<<"-----------------------"<<endl;
+               cout<<fRootFile<<" opened!"<<endl;
+               cout<<"-----------------------"<<endl;
+       }
     return 0;    
 }
 //####################################################################
@@ -816,12 +825,13 @@ bool MAPS::InitialDynNoise(Int_t startframe, Int_t frames) {
         
         // first estimate
         // sum up over all frames every CDS value, pixelwise
+               
         for(Int_t framei=startframe; framei<frames+startframe; framei++)
         {
             getFrame(framei);
-//             cout << "RAW CDS signal, frame: "<< framei;;
-//             debugStream(fCdsmatrix, fPixels, fColumns,3);
-            
+             // cout << "RAW CDS signal, frame: "<< framei;;
+             // debugStream(fCdsmatrix, fPixels, fColumns,3);
+           
             for(Int_t pixeli=0; pixeli<fPixels; pixeli++)
             {
                 fPedestals[pixeli] += fCdsmatrix[pixeli];
@@ -846,7 +856,41 @@ bool MAPS::InitialDynNoise(Int_t startframe, Int_t frames) {
         {
             fNoise[pixeli] = TMath::Sqrt(fNoise[pixeli]/(frames-1)); // Standard deviation is defined with n-1 
         }
-        
+        // ----------------------
+            // debugStream<>(fPedestals , fPixels, 1664, 1, 14);
+                       // cout<<"MedianAll"<<TMath::Median(fPixels,fPedestals)<<endl;
+                       // Float_t a0a[208];
+                       // Float_t a0b[208];
+                       // for(int i=0;i<208;i++)
+                       // {
+                       // a0a[i]=fPedestals[i];
+                       // a0b[i]=fPedestals[208+i];
+
+                       // }
+                       // for(int ichannel=0;ichannel<4;ichannel++)
+                       // {
+                               // cout<<"Pedestal A0[ "<<ichannel<<"]"<<endl;
+                               // for(int i=0;i<208;i++)
+                               // {
+                               // cout<<fPedestals[i+ichannel*208]<<",";
+                               // }
+                               // cout<<endl;
+                       // }
+                       // for(int ichannel=0;ichannel<4;ichannel++)
+                       // {
+                               // cout<<"Pedestal A1[ "<<ichannel<<"]"<<endl;
+                               // for(int i=0;i<208;i++)
+                               // {
+                               // cout<<fPedestals[i+ichannel*208+4*208]<<",";
+                               // }
+                               // cout<<endl;
+                       // }
+                       // cout<<"MedianAll: "<<TMath::Median(fPixels,fPedestals)<<endl;
+                       // cout<<"Mediana0a: "<<TMath::Median(208,a0a)<<endl;
+                       // cout<<"Mediana0b: "<<TMath::Median(208,a0b)<<endl;
+                       // exit(0);
+                       
+                       // -----------
         // second estimate
         // cut away pixels with cds - Pedestals > 5 * sigma
         // sum up over all frames every CDS value, pixelwise
@@ -964,7 +1008,7 @@ bool MAPS::InitialDynNoise(Int_t startframe, Int_t frames) {
             Float_t PEDESTAL;
             Float_t NOISE;
             Int_t PIXEL;
-            
+                       
             fNoiseTree->Branch("pixel"      , &PIXEL    , "pixel/i"     , 32000);
             fNoiseTree->Branch("noise"      , &NOISE    , "noise/F"     , 32000);
             fNoiseTree->Branch("pedestal"   , &PEDESTAL , "pedestal/F"  , 32000);
@@ -1256,7 +1300,19 @@ void MAPS::hitana() {
                        cdsmatrixCorrPed[i]=(float)(1.*fCdsmatrix[i]-fPedestals[i]);
             fHittedPixel[i] = 0;
         //    if( (float)(1.*fCdsmatrix[i]-fPedestals[i]) > (50.) )
-            if( (float)(1.*fCdsmatrix[i]-fPedestals[i]) > (5.*fNoise[i]) )
+                       if(fOrderCode=="FSBB")
+                       {
+                               if( (float)(1.*fCdsmatrix[i]-fPedestals[i]) > (5.*fNoise[i]) )
+                               {
+                HITNR++;
+                HITS.Set(HITNR);
+                HITS.AddAt(i,(HITNR-1));                
+                if(fSave)   {
+                    hint2->Fill(i%fColumns+0.1, (int)(i/fColumns)+0.1); // counts up in 2dimensional pixel matrix
+                }
+                               }
+                       }
+                       else if( (float)(1.*fCdsmatrix[i]-fPedestals[i]) > (5.*fNoise[i]) )
             {
                 HITNR++;
                 HITS.Set(HITNR);
@@ -1600,7 +1656,8 @@ void MAPS::hitana() {
             //Begin: loop evaluate true seeds:
             if(CHANCE==100)
             {
-                       //cout<<"CLUSTER ACCEPTED"<<Hitlist[hit]<< " with ADC: " << pixelchargeincluster[12]<<endl;
+               //      debugStream<>(fCdsmatrix, fPixels, fColumns, 1, 100);
+               //      cout<<"CLUSTER ACCEPTED"<<Hitlist[hit]<< " with ADC: " << pixelchargeincluster[12]<<" A: "<<A<<" B: "<<B<<endl;
                        if(fOrderCode=="FSBB")
                        {
                                if (B < 2 || B > fColumns-3 || A < 1 || A > fRows-2)
@@ -1642,7 +1699,8 @@ void MAPS::hitana() {
                 else
                 {
                     fHittedPixel[Hitlist[hit]] = 2;
-                
+                               //      if(Hitlist[hit]<415&&Hitlist[hit]>830)
+                               //      cout<<fFrameNumber<<","<<Hitlist[hit]<<endl;
                     //Fill hit TTree:
                     fFrameInfo.pixel[fHits] = Hitlist[hit];
                     for(int clupos=0; clupos<25; clupos++)
@@ -1786,16 +1844,17 @@ void MAPS::plotFrame(Int_t FrameNumber) {
        
         cout <<endl <<endl << colorwhite << "---------------- FRAME " << fFrameNumber << " ----------------" << colorreset << endl << endl;
         cout<<"CDS matrix:"<<endl;
-       // debugStream<>(fCdsmatrix, fPixels, fColumns, 0, 39);
+        debugStream<>(fCdsmatrix, fPixels, fColumns, 0, 39);
+
                Float_t cdsmatrixCorrPed[fPixels];
                for(Int_t i=0; i<fPixels; i++)
         {
                        cdsmatrixCorrPed[i]=(float)(1.*fCdsmatrix[i]-fPedestals[i]);
                }
-       //      debugStream<>(cdsmatrixCorrPed, fPixels, fColumns, 0, 10);
+               debugStream<>(cdsmatrixCorrPed, fPixels, fColumns, 0, 10);
                // cout<<fCdsmatrix[0]<<","<<fCdsmatrix[1]<<","<<fCdsmatrix[2]<<","<<fCdsmatrix[3]<<","<<fCdsmatrix[4]<<","<<fCdsmatrix[5]<<","<<fCdsmatrix[6]<<endl;
         //cout<<"Hitted pixel discriminator matrix:"<<endl;
-        //debugStream<>(fHittedPixel, fPixels, fColumns, 1, 1);
+    //    debugStream<>(fHittedPixel, fPixels, fColumns, 1, 1);
                
         cout << colorred <<"Accepted hit cluster:"<<endlr;
                for (Int_t hiti = 0; hiti < fHits; hiti++)
@@ -2714,7 +2773,7 @@ void MAPS::reorderFSBB3() {
 //####################################################################
 
 //####################################################################
-void MAPS::reorderFSBB() {
+void MAPS::reorderFSBB_B() {
                        
                Float_t CDSMATRIX       [fPixels];
         Int_t F0MATRIX [fPixels];
@@ -2762,9 +2821,9 @@ void MAPS::reorderFSBB() {
                           a1[iOutChannel][iPixelInZeile]=CDSMATRIX     [(4+iOutChannel)*subColumns+iPixelInZeile];//4+iOutChannel, because this is A1, which is shifted by four channels of A0
                        }
                        
-                       
+               /*      
                        // Start Filter CommonMode
-                       /*Float_t       CommonModeInRow = 0;
+                       Float_t CommonModeInRow = 0;
                        bool warning=false;
                        CommonModeInRow = TMath::Median(subColumns, a0[iOutChannel]); // Don't use mean, consider hitted pixel !
                        if (CommonModeInRow < 0)
@@ -2772,7 +2831,7 @@ void MAPS::reorderFSBB() {
                        else
                                minus++;
                        
-                       if(TMath::Abs(CommonModeInRow)>0) // TODO better warning criteria
+                       if(TMath::Abs(CommonModeInRow)>10000) // TODO better warning criteria
                        {
                                warning = true;
                        }
@@ -2785,6 +2844,27 @@ void MAPS::reorderFSBB() {
                        {
                                a0[iOutChannel][column ] -= CommonModeInRow;
                        }
+                       CommonModeInRow = 0;
+                       warning=false;
+                       CommonModeInRow = TMath::Median(subColumns, a1[iOutChannel]); // Don't use mean, consider hitted pixel !
+                       if (CommonModeInRow < 0)
+                               plus++;
+                       else
+                               minus++;
+                       
+                       if(TMath::Abs(CommonModeInRow)>10000) // TODO better warning criteria
+                       {
+                               warning = true;
+                       }
+                       if(warning) { cout<<"\rFrame: "<<fFrameNumber<< " row: " << iOutChannel << " Median of row: " << CommonModeInRow <<  " --> Common Mode suspiciously high!               "<<endl; 
+                                debugStream<>(fCdsmatrix, fPixels, fColumns, 0, 20);
+                               debugStream<>(a1[iOutChannel], subColumns, subColumns, 1, 20);
+                       }
+
+                       for(int column=0; column<subColumns; column++ )
+                       {
+                               a1[iOutChannel][column ] -= CommonModeInRow;
+                       }
                        // End Filter CommonMode*/
                }
                
@@ -2826,38 +2906,38 @@ void MAPS::reorderFSBB() {
                                
                        }
                }
-               //debugStream<>(fCdsmatrix, fPixels, fColumns, 0, 20);
+       //      debugStream<>(fCdsmatrix, fPixels, fColumns, 0, 20);
+               Float_t tMatrix [fPixels];
                for(Int_t iPixel=0; iPixel<fPixels;iPixel++)
                        { 
                //              cout<<fCdsmatrix[iPixel]<<","; 
-                               
+                               tMatrix[iPixel]=iPixel;
+                               //cout<<tMatrix[iPixel]<<","; 
                        }
+                       //debugStream<>(tMatrix, fPixels, fColumns, 0, 20);
                //exit(0);
 }
 
 //####################################################################
 //####################################################################
-void MAPS::reorderFSBB_B() {
-                       // War: A0-0, A1-0, A0-1, A1-1
-                       //0,1,...,415
-                       //416,...,
-                       //Ist Blödsinn, deshalb wird zuerst A0 ausgelesen, und dann A1
+void MAPS::reorderFSBB() {
+                       
                Float_t CDSMATRIX       [fPixels];
         Int_t F0MATRIX [fPixels];
         Int_t F1MATRIX [fPixels];
                for(Int_t i=0; i<fPixels; i++)
         {
-            CDSMATRIX  [i] = i;
+            CDSMATRIX  [i] = fCdsmatrix[i];
             F0MATRIX   [i] = fF0matrix [i];
             F1MATRIX   [i] = fF1matrix [i];
         }
                fRows           = 4;
         fColumns       = 416;
         fPixels                = fRows*fColumns;
-       //      debugStream(CDSMATRIX, fPixels, 416, 0, 2000);
-       
+       //      debugStream(CDSMATRIX, fPixels, fPixels, 0, 14);
+       //exit(0);
                Int_t subRows = 4;
-               Int_t subColumns = 2*208;
+               Int_t subColumns = 208;
                Int_t a0[subRows][subColumns];
                Int_t a1[subRows][subColumns];
                Int_t aMerge[fRows][fColumns];
@@ -2881,16 +2961,16 @@ void MAPS::reorderFSBB_B() {
                //Divide CDSMATRIX in the channels. First 208 pixels are from A0-0, then A0-1, then A0-2, then A0-3, then A1-0, then A1-1, then A1-2, then A1-3
                 for(Int_t iOutChannel=0; iOutChannel<subRows; iOutChannel++)
                 {
-                       for(Int_t iPixelInZeile=0; iPixelInZeile<fColumns; iPixelInZeile++)
+                       for(Int_t iPixelInZeile=0; iPixelInZeile<subColumns; iPixelInZeile++)
                        {
                        //      cout<<iPixelInZeile<<":"<<(4+iOutChannel)*subColumns+iPixelInZeile<<"-";
-                          a0[iOutChannel][iPixelInZeile]=CDSMATRIX     [iOutChannel*fColumns+iPixelInZeile];
-                       //   a1[iOutChannel][iPixelInZeile]=CDSMATRIX   [(4+iOutChannel)*subColumns+iPixelInZeile];//4+iOutChannel, because this is A1, which is shifted by four channels of A0
+                          a0[iOutChannel][iPixelInZeile]=CDSMATRIX     [iOutChannel*subColumns+iPixelInZeile];
+                          a1[iOutChannel][iPixelInZeile]=CDSMATRIX     [(4+iOutChannel)*subColumns+iPixelInZeile];//4+iOutChannel, because this is A1, which is shifted by four channels of A0
                        }
                        
-                       
+               /*      
                        // Start Filter CommonMode
-                       /*Float_t       CommonModeInRow = 0;
+                       Float_t CommonModeInRow = 0;
                        bool warning=false;
                        CommonModeInRow = TMath::Median(subColumns, a0[iOutChannel]); // Don't use mean, consider hitted pixel !
                        if (CommonModeInRow < 0)
@@ -2898,7 +2978,7 @@ void MAPS::reorderFSBB_B() {
                        else
                                minus++;
                        
-                       if(TMath::Abs(CommonModeInRow)>0) // TODO better warning criteria
+                       if(TMath::Abs(CommonModeInRow)>10000) // TODO better warning criteria
                        {
                                warning = true;
                        }
@@ -2911,10 +2991,31 @@ void MAPS::reorderFSBB_B() {
                        {
                                a0[iOutChannel][column ] -= CommonModeInRow;
                        }
+                       CommonModeInRow = 0;
+                       warning=false;
+                       CommonModeInRow = TMath::Median(subColumns, a1[iOutChannel]); // Don't use mean, consider hitted pixel !
+                       if (CommonModeInRow < 0)
+                               plus++;
+                       else
+                               minus++;
+                       
+                       if(TMath::Abs(CommonModeInRow)>10000) // TODO better warning criteria
+                       {
+                               warning = true;
+                       }
+                       if(warning) { cout<<"\rFrame: "<<fFrameNumber<< " row: " << iOutChannel << " Median of row: " << CommonModeInRow <<  " --> Common Mode suspiciously high!               "<<endl; 
+                                debugStream<>(fCdsmatrix, fPixels, fColumns, 0, 20);
+                               debugStream<>(a1[iOutChannel], subColumns, subColumns, 1, 20);
+                       }
+
+                       for(int column=0; column<subColumns; column++ )
+                       {
+                               a1[iOutChannel][column ] -= CommonModeInRow;
+                       }
                        // End Filter CommonMode*/
                }
                
-                // for(Int_t i=0; i<208*2; i++)
+                // for(Int_t i=0; i<208; i++)
         // {
           // cout<< a0[0][i]<<",";
         // }
@@ -2936,11 +3037,11 @@ void MAPS::reorderFSBB_B() {
                // }
                for(Int_t iOutChannel=0; iOutChannel<subRows; iOutChannel++)
         {
-                       for(Int_t iMerge2channels=0; iMerge2channels<fColumns;iMerge2channels++)
+                       for(Int_t iMerge2channels=0; iMerge2channels<subColumns;iMerge2channels++)
                        {
-                       //Merge A0-0[0],A1-0[0],A0-0[1],A1-0[1],A0-0[2],A1-0[2],... output is an array of 4 geometrically columns
+                       //Merge A0-0, then A1-0 in one line, then A0-1 and A1-1 in one line...
                                aMerge[iOutChannel][iMerge2channels]=a0[iOutChannel][iMerge2channels];//A0-iOutChannel[iMerge2channels]
-                       //      aMerge[iOutChannel][iMerge2channels+subColumns]=a1[iOutChannel][iMerge2channels];//A1-iOutChannel[iMerge2channels]
+                               aMerge[iOutChannel][iMerge2channels+subColumns]=a1[iOutChannel][iMerge2channels];//A1-iOutChannel[iMerge2channels]
                        }
                }
                
@@ -2952,12 +3053,15 @@ void MAPS::reorderFSBB_B() {
                                
                        }
                }
-               debugStream<>(fCdsmatrix, fPixels, fColumns, 0, 20);
+       //      debugStream<>(fCdsmatrix, fPixels, fColumns, 0, 20);
+               Float_t tMatrix [fPixels];
                for(Int_t iPixel=0; iPixel<fPixels;iPixel++)
                        { 
                //              cout<<fCdsmatrix[iPixel]<<","; 
-                               
+                               tMatrix[iPixel]=iPixel;
+                               //cout<<tMatrix[iPixel]<<","; 
                        }
+                       //debugStream<>(tMatrix, fPixels, fColumns, 0, 20);
                //exit(0);
 }
 

diff --git a/MABS_run_analyzer/MAPS.h b/MABS_run_analyzer/MAPS.h
index 05032f47886a8127dc1a3f48f5f1ee1d28cf31ae..ecdc529749c5f580cd519da921a155c5e5c9cf75 100644 (file)
--- a/MABS_run_analyzer/MAPS.h
+++ b/MABS_run_analyzer/MAPS.h
@@ -83,9 +83,7 @@ private:
     UInt_t       fFrameNumber; /**< enum value 1 */
     /// true if all data files OK
     bool        fOk; 
-    /// if set to true, a root file for the analyzed run will be created, set and passed initMapsRun() to  in the constructor
-    bool        fSave = true;
-    /// if set to true, the frame retrieved with @c getframe() is OK
+        /// if set to true, the frame retrieved with @c getframe() is OK
     bool        fFrameOk;
     /// if set to true, the file from where the noise is extracted is OK
     bool        fNoiseOk;
@@ -322,7 +320,8 @@ public:
     
     /// is set to true if somewhere a critical error occurs
     Bool_t error = false;
-    
+    /// if set to true, a root file for the analyzed run will be created, set and passed initMapsRun() to  in the constructor
+    bool        fSave = true;
     
     /// is set to true if somewhere a critical error occurs
     Bool_t used_default_config = false;

diff --git a/MABS_run_analyzer/Run.c b/MABS_run_analyzer/Run.c
index 9bd22c86bd858a438977fbe051abec3889b5decc..3c1b8f178584111651ef45b650d6a7c144fb6cfc 100644 (file)
--- a/MABS_run_analyzer/Run.c
+++ b/MABS_run_analyzer/Run.c
@@ -269,8 +269,8 @@ Bool_t Run::analyzeRun(Bool_t force)
             processed->InitialDynNoise();
             int start   = 0;
             int nframes = processed->GetNumberFrames();
-             for(int i=0; i<1000000;i++) // TODO remove 100000 run 342272
-         //  for(int i=0; i<nframes;i++) // TODO remove 100000 run 342272
+         //    for(int i=0; i<1000000;i++) // TODO remove 100000 run 342272
+           for(int i=0; i<nframes;i++) // TODO remove 100000 run 342272
             {
 //                 cout << "getframe " << i << endl;
                 processed->getFrame(i);
@@ -443,25 +443,28 @@ Bool_t Run::analyzeFrame(Int_t frame)
     if (!error)
     {
         processed = new MAPS(this);
-        processed->initNewRootFile();
-        int entries = processed->GetNumberFrames();
-        if (frame < entries)
-        {      
+               
+            processed->initOldRootFile();
+        //processed->plotPixSignal  (0, 100000, 350);
+               
+               int entries = processed->GetNumberFrames();
+        // if (frame < entries)
+        // {   
                        if (dynamicalNoise)
             processed->InitialDynNoise(100);
             processed->getFrame(frame);
             // processed->filterCommonMode();
             processed->hitana();
             processed->plotFrame(frame);
-            delete processed;
-            return 0;
-        }
-        else
-        {
-            cout << "\033[1;33mFrame number too big, max frame: " << entries << "\033[0m" << endl;
-        }
+          //  delete processed;
+            // return 0;
+        // }
+        // else
+        // {
+            // cout << "\033[1;33mFrame number too big, max frame: " << entries << "\033[0m" << endl;
+        // }
     }
-    delete processed;
+  //  delete processed;
     return 1;
 }
 
@@ -755,11 +758,21 @@ 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->fNoiseTree->SetBranchAddress("noise", &noise, &noiseBranch);    
-    for (Int_t cnt=0; cnt<processed->fNoiseTree->GetEntries(); cnt++) {
-        processed->fNoiseTree->GetEntry(cnt);
-        histogram.Noise->Fill(noise);
-    }    
+    processed->fNoiseTree->SetBranchAddress("noise", &noise, &noiseBranch); 
+       for (Int_t cnt=0; cnt<processed->fNoiseTree->GetEntries(); cnt++) {     
+               uint pixel_column_x = cnt%sensorinfocur.columns;    //  column of seed
+               uint pixel_row_y = cnt/sensorinfocur.columns;       //  row of seed
+               if (false) {
+                    cout << "submatrix_x_start: " << submatrix_x_start << ", submatrix_x_end: " << submatrix_x_end << endl;
+                    cout << "submatrix_y_start: " << submatrix_y_start << ", submatrix_y_end: " << submatrix_y_end << endl;
+                    cout << "pixel_column_x: " << pixel_column_x << ", pixel_row_y: " << pixel_row_y << endl;
+               }
+               if (pixel_column_x >= submatrix_x_start && pixel_column_x < submatrix_x_end && pixel_row_y >= submatrix_y_start && pixel_row_y < submatrix_y_end) // Diode sitzt oben im SeedPixel, da nach PitchY angeordnet
+               {
+                               processed->fNoiseTree->GetEntry(cnt);
+                               histogram.Noise->Fill(noise);
+               }
+       }    
     histogram.Noise->GetQuantiles( 3, noisequantiles, probabilities);
     histogram.avgNoise = noisequantiles[1];
     histogram.avgNoisePlus = noisequantiles[2] - noisequantiles[1];
@@ -783,6 +796,8 @@ Bool_t Run::binSeedSumVeto()
     /// collected charge in cluster
     Float_t pixelSum = 0;
     Float_t notSeedSum = 0;
+       getPixelWithFakehits();
+       
    //  for (Int_t framei=0; framei<10000; framei++) // loop over all frames   
     for (Int_t framei=0; framei<processed->fHitTree->GetEntries(); framei++) // loop over all frames
     {
@@ -820,7 +835,7 @@ Bool_t Run::binSeedSumVeto()
                                 notSeedSum += processed->fFrameInfo.p[clusteri][hiti];
                         }
                         TMath::Sort(processed->clustersize*processed->clustersize,clusterArray,index,1);
-                        for (Int_t clusteri=0; clusteri<4; clusteri++)
+                        for (Int_t clusteri=0; clusteri<21; clusteri++)
                         {
                             pixelSum += clusterArray[index[clusteri]];
                             
@@ -843,6 +858,22 @@ Bool_t Run::binSeedSumVeto()
                     
                     if (processed->fFrameInfo.pixelthreshold[hiti]>0)
                     {
+                                                Bool_t pixelfound=false;
+                                                for(int col=0;col<fcountPixelFakeHits;col++)
+                                                {
+                                                if(fPixelFakeHits[col]==processed->fFrameInfo.pixel[hiti])
+                                                       {
+                                                        pixelfound=true;
+                                                        }
+                                                }
+                                                if(pixelfound==true)
+                                                {
+                                                       continue;//If pixel is a fake hit pixel with a rate >1e-5, do not count it
+                                                }
+                                                 // if(pixelfound==false)
+                                                // {
+                                                       // cout<<"Pixel"<<processed->fFrameInfo.pixel[hiti]<<" NOT found in array"<<endl;
+                                                // }
                         histogramthreshold.numberofhits++;
                         
                         histogramthreshold.Seed->Fill(processed->fFrameInfo.p[12][hiti]);
@@ -868,20 +899,20 @@ Bool_t Run::binSeedSumVeto()
                         if (TMath::Abs(notSeedSum) < systemparamcur.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd")))
                             histogramthreshold.Veto->Fill(processed->fFrameInfo.p[12][hiti]);    // histogram with the single pixel
                                                        
-                                                if (TMath::Abs(pixelSum) > 200 && (labbook.source.Contains("Fe")))
-                                                { 
-                            histogramthreshold.VetoSum->Fill(processed->fFrameInfo.p[12][hiti]);    // histogram with the single pixel
+                                                //if (TMath::Abs(pixelSum) > 200 && (labbook.source.Contains("Fe")))
+                                                //{ 
+                          //  histoPixel->Fill(processed->fFrameInfo.pixel[hiti]);    // histogram with the single pixel
+                                                       
                                                        //histogramthreshold.VetoSum->Fill(pixelSum);
                                                        // cout<<"Hitnummer"<<framei<<":"<<hiti<<": ";
-                                                       // for(int col=0;col<25;col++)
-                                                       // cout<<processed->fFrameInfo.p[col][hiti]<<",";
-                                                       // cout<<endl;
-                                                }
+                                                        
+                                               // }
                     }
                 }
             }
         }
     }
+
 //     gROOT->SetBatch(kTRUE);
     if (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))
         histogram.posVeto=FitPerform(histogram.Veto, "gaus", true);
@@ -2004,5 +2035,43 @@ void Run::initRootParameters()
     rootlinestyle = new Int_t[13]{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
     
 }
-
+void Run::getPixelWithFakehits()
+{      
+       Int_t pixels=sensorinfocur.columns*sensorinfocur.rows;
+       histoPixel = new TH1F("pixel histo","pixel histo",pixels,0,pixels);
+       for (Int_t framei=0; framei<processed->fHitTree->GetEntries(); framei++) // loop over all frames
+    {
+        processed->fHitTree->GetEntry(framei);
+               for(Int_t hiti=0; (unsigned int)hiti<processed->fFrameInfo.hits;hiti++)
+            {
+                               histoPixel->Fill(processed->fFrameInfo.pixel[hiti]); 
+                       }
+       }
+       fPixelFakeHits  = new Float_t [pixels]();
+       fcountPixelFakeHits=0;
+       for(int col=0;col<pixels;col++)
+       {
+               Float_t fakehitrate=0.00001;
+               Float_t sourceIntensity=1;//Adjust it depending on the intensity of the used source
+               if(labbook.source.Contains("Fe"))
+               {
+                       sourceIntensity=30;
+               }
+               else if(labbook.source.Contains("Sr"))
+               {
+                       sourceIntensity=2;
+               }
+               else if(labbook.source.Contains("Cd"))
+               {
+                       sourceIntensity=2;
+               }
+               if(histoPixel->GetBinContent(col)>fakehitrate*sourceIntensity*labbook.frames_foundDB)//show only pixels with a fake hit rate above 1e-5
+               {       
+                       fPixelFakeHits[fcountPixelFakeHits]=col;
+                               fcountPixelFakeHits++;
+                        cout<<col<<":"<<histoPixel->GetBinContent(col)<<endl;
+               }
+       }
+return;        
+}
 #endif

diff --git a/MABS_run_analyzer/Run.h b/MABS_run_analyzer/Run.h
index dbafc9d70ada554e0426be02f269eeb97a74be57..bee99901a1cbf29c39548745de7ab84e8014b6fd 100644 (file)
--- a/MABS_run_analyzer/Run.h
+++ b/MABS_run_analyzer/Run.h
@@ -279,8 +279,12 @@ private:
     };
     Int_t* rootcolors;
     Int_t* rootlinestyle;
-        
-    /// init one specific histogram
+       
+        void getPixelWithFakehits();
+       TH1F* histoPixel; 
+    Float_t* fPixelFakeHits;
+       Int_t fcountPixelFakeHits;
+       /// init one specific histogram
     void initHistogram(TH1F* &histogrampointer, TString prefix);
     
     /** @brief calls #initCluster for each structure of type @c Run::histogram
@@ -319,7 +323,7 @@ public:
     Run    ( Int_t, Int_t );
     /** @brief writes values back to the database, destroys the #MAPS object */
     ~Run   (void);
-
+       
     systemparam systemparamcur; 
     
     void setPlotStyle(Int_t);