Changeset 15368 for fact

Timestamp:

04/19/13 01:37:40 (12 years ago)

Author:

Jens Buss

Message:

add functions

File:

• fact/tools/rootmacros/zerosearch.C (modified) (12 diffs)

fact/tools/rootmacros/zerosearch.C

@@ -19 +19 @@
     int VerbosityLevel
 ){
-
 vector<Region> * ZeroPositions = new vector<Region>;
-Region currentRegion = {0, 0, 0, 0.0, 0, 0.0, 0.0, 0};
+Region currentRegion = {0, 0, 0, 0.0, 0, 0.0, 0.0, 0, 0};
+
     if (step < 1){
         if (VerbosityLevel > 0)
@@ -67 +67 @@
             }
         }
-
     }
     else // search for zero-x-ings on both esges
@@ -81 +80 @@
         }
     }
-
 
     return ZeroPositions;
@@ -172 +170 @@
                 << reg->maxPos << " with Value:" << reg->maxVal << endl;
         }
+
+        if (reg->maxPos <= 0)
+        {
+            reg = regions.erase( reg ) ;
+            --reg;
+            continue;
+        }
     }
     return regions.size();
@@ -411 +416 @@
 */
 
+//size_t calcAmplitudeWeightedTime(
+//        vector<Region>  &regions,
+//        vector<float>   &data,
+
+//        ){
+
+//}
+
+size_t calcCFDPositions(
+        vector<Region>  &regions,
+        vector<float>   &cfd_data
+        )
+{
+    vector<Region>::iterator reg;
+    for (reg = regions.begin() ; reg < regions.end() ; reg++)
+    {
+        int     neg_zero_crossing = -1;
+        int     min_in_cfd        = -1;
+        int     pos_zero_crossing = reg->maxPos;
+        float   minVal            = 10.0;
+
+        for (int i = pos_zero_crossing; i > pos_zero_crossing - 40; i--)
+        {
+            if (i <= 1 || i >= 1439) break;
+
+            if ( cfd_data[i] < minVal)
+            {
+                minVal      = cfd_data[i];
+                min_in_cfd  = i;
+            }
+
+            if (cfd_data[i] > 0 && cfd_data[i+1] <= 0 && cfd_data[i+2] < 0 )
+            {
+                neg_zero_crossing = i+1;
+            }
+        }
+
+        reg->cfdMinPos          = min_in_cfd;
+        reg->cfdNegZerocrossing = neg_zero_crossing;
+    }
+    return regions.size();
+}
+
+
 
 size_t findTimeOfHalfMaxLeft(
-    vector<Region> &regions,
-    vector<float> &data,
-    float baseline,
-    int beginRisingEdge,
-    int endRisingEdge,
-    int VerbosityLevel)
-{
-    float pulse_size = 0;     //
-    float thr = 0;
-
+    vector<Region>  &regions,
+    vector<float>   &data,
+    float           baseline,
+    int             beginRisingEdge,
+    int             endRisingEdge,
+    int             VerbosityLevel
+        )
+{
+//    float   pulse_size      = 0;     //
+    float   thr             = 0;
     //int begin = beginRisingEdge;
-    int end = endRisingEdge;
-    int counter = 0;
+    int     end             = 20;
+    int     counter         = 0;
     // local vars for line fitting
-    float xmean = 0.0;
-    float ymean = 0.0;
-    float cov = 0.0;          // empiric covarianve between x and y
-    float var = 0.0;          // empiric variance of x
+    float   xmean           = 0.0;
+    float   ymean           = 0.0;
+    float   cov             = 0.0;          // empiric covariance between x and y
+    float   var             = 0.0;          // empiric variance of x
     // result of line fitting
-    float slope = 0.0;
-    float intercept = 0.0;
+    float   slope           = 0.0;
+    float   intercept       = 0.0;
     // result of this function -- the position of the half rising edge
-    float pos_of_thr_Xing = 0.0;
-    int result = 0;
+    float   pos_of_thr_Xing = 0.0;
+    int     result          = 0;
+
 
     vector<Region>::iterator reg;
-    for (reg = regions.begin() ; reg < regions.end() ; reg++){
-        counter = 0;
+    for (reg = regions.begin() ; reg < regions.end() ; reg++)
+    {
         // check if linear rising edge is completely in pipeline
         // if not --> delete
@@ -450 +500 @@
             continue;
         }
-
-        // The maximum was probably found using smoothed data,
-        // so I read the value at the position of the maximum from the data
-        // again. So I rely on a well defined maxPos.
-        if (VerbosityLevel > 1) cout << "## baseline: " << baseline << endl;
-        pulse_size = data[reg->maxPos] - baseline;
-        if (VerbosityLevel > 1) cout << "## pulse_size: " << pulse_size << endl;
-        thr = pulse_size / 2. + baseline;
-        if (VerbosityLevel > 1) cout << "## thr: " << thr << endl;
 
         // I think 5 slices left of the max there begins the rising edge
@@ -469 +510 @@
         // ALARM check for out of range values...
 
+        float average_of_max = -1;
+        int stepwidth = 2;
+
+        counter = 0;
+        for ( int i = reg->maxPos -stepwidth; i <= reg->maxPos + stepwidth; i++)
+        {
+            average_of_max += data[i];
+            counter++;
+        }
+        average_of_max /= counter;
+
+        thr = average_of_max / 2.0;
+
+        float   max_slope       = 0;
+        int     max_slope_pos   = -1;
+        int     half_height_pos = -1;
+
         bool passed = false;
-        for (int slice=reg->maxPos;
-                    slice > 1; --slice)
-        {
-            if ( data[slice] < 0.9*data[reg->maxPos] && !passed)
-            {
-                beginRisingEdge = reg->maxPos - slice;
-                passed = true;
-            }
-
-            if ( data[slice] < 0.1*data[reg->maxPos])
-            {
-                endRisingEdge = reg->maxPos - slice;
-                passed = false;
-                reg->lengthOfRisingEdge=endRisingEdge - beginRisingEdge;
-                break;
-            }
-        }
-
-
+        for (int slice=reg->maxPos; slice > reg->maxPos - end; --slice)
+        {
+            if (slice <= 0) break;
+//            if ( data[slice] < 0.9*data[reg->maxPos] && !passed)
+//            {
+//                endRisingEdge = slice;
+//                beginRisingEdge = endRisingEdge -3;
+//                passed = true;
+//            }
+//            else if ( data[slice] < 0) //0.4*data[reg->maxPos] )
+//            {
+//                beginRisingEdge = slice;
+//                reg->lengthOfRisingEdge=endRisingEdge - beginRisingEdge;
+//                break;
+//            }
+
+            //calculate max slope on leading edge
+            if (data[slice] - data[slice-4] > max_slope && slice-4 > 0)
+            {
+                max_slope   = data[slice] - data[slice-4];
+                max_slope_pos      = slice-2;
+            }
+
+            if ( data[slice] <= thr && !passed){
+                half_height_pos  = slice;
+                passed  = true;
+            }
+        }
+        reg->maxSlope           = max_slope/4;
+        reg->maxSlopePos        = max_slope_pos;
+        reg->halfRisingEdgePos  = half_height_pos;
+        reg->halfHeight         = data[half_height_pos];
+
+        result                  = half_height_pos;
+
+        endRisingEdge   = result +2;
+        beginRisingEdge = result -2;
+
+        counter = 0;
         //calculate mean of x and y coordinate
-        for (int slice=reg->maxPos - beginRisingEdge;
-                    slice > reg->maxPos - endRisingEdge; --slice)
-        {
+        for (int slice= beginRisingEdge; slice < endRisingEdge; ++slice)
+        {
+            if (slice <= 0) break;
+
             xmean += slice;
             ymean += data[slice];
             counter++;
         }
-//        xmean /= beginRisingEdge - endRisingEdge;
-//        ymean /= beginRisingEdge - endRisingEdge;
         xmean /= counter;
         ymean /= counter;
@@ -508 +585 @@
         var = 0.0;
 
-        for (int slice=reg->maxPos - beginRisingEdge;
-                    slice > reg->maxPos - endRisingEdge; --slice)
-        {
+        //calculate covariance and variance
+        for (int slice=beginRisingEdge; slice < endRisingEdge; ++slice)
+        {
+            if (slice <= 0) break;
+
             cov += (data[slice] - ymean) * (slice - xmean);
             var += (slice - xmean) * (slice - xmean);
@@ -517 +596 @@
         if (VerbosityLevel > 2) cout << "## var: " << var << endl;
 
-        slope = cov / var;
-        intercept = ymean - slope * xmean;
+        slope       = cov / var;
+        intercept   = ymean - slope * xmean;
+
+        // The maximum was probably found using smoothed data,
+        // so I read the value at the position of the maximum from the data
+        // again. So I rely on a well defined maxPos.
+//        if (VerbosityLevel > 4) cout << "## baseline: " << baseline << endl;
+//        pulse_size = data[reg->maxPos] - baseline;
+//        if (VerbosityLevel > 4) cout << "## pulse_size: " << pulse_size << endl;
+//        thr = pulse_size / 2. + baseline;
+//        if (VerbosityLevel > 4) cout << "## thr: " << thr << endl;
 
         // now calculate, where the fittet line crosses the threshold
-        pos_of_thr_Xing = (thr - intercept) / slope;
-        result = (int)(pos_of_thr_Xing + 0.5);
+//        pos_of_thr_Xing = (thr - intercept) / slope;
+//        result = (int)(pos_of_thr_Xing + 0.5);
 
 
@@ -537 +625 @@
         }
 
-        reg->halfRisingEdgePos = result;
-
-        reg->distanceEdgeToMax = reg->maxPos - result;
-        reg->interceptRisingEdge = intercept;
-        reg->slopeOfRisingEdge = slope;
+        reg->distanceEdgeToMax      = reg->maxPos - result;
+        reg->interceptRisingEdge    = reg->maxPos - intercept;
+        reg->slopeOfRisingEdge      = slope;
+
+        if (reg->halfRisingEdgePos <= 0 || reg->halfRisingEdgePos >= 1024)
+        {
+            reg = regions.erase( reg ) ;
+            --reg;
+            continue;
+        }
     }
     return regions.size();
@@ -570 +663 @@
 size_t GetMaxPositions(
     vector<Region>  &regions,
-    vector<int>     &maxPositions,
-    int VerbosityLevel)
-{
-    maxPositions.clear();
+    vector<int>     &positions,
+    int VerbosityLevel)
+{
+    positions.clear();
     vector<Region>::iterator it = regions.begin();
     while( it != regions.end() )
     {
-        maxPositions.push_back(it->maxPos);
+        positions.push_back(it->maxPos);
         if (VerbosityLevel > 3){
             cout << "getting maxPos@ " << it->maxPos << "\t";
@@ -590 +683 @@
 size_t GetEdgePositions(
     vector<Region>  &regions,
-    vector<int>     &edgePos,
-    int VerbosityLevel)
-{
-    maxPositions.clear();
+    vector<int>     &positions,
+    int VerbosityLevel)
+{
+    positions.clear();
     vector<Region>::iterator it = regions.begin();
     while( it != regions.end() )
     {
-        maxPositions.push_back(it->halfRisingEdgePos);
+        positions.push_back(it->halfRisingEdgePos);
         if (VerbosityLevel > 3){
             cout << "getting maxPos@ " << it->halfRisingEdgePos << "\t";