Ignore:
Timestamp:
04/19/13 01:37:40 (12 years ago)
Author:
Jens Buss
Message:
add functions
File:
1 edited

Legend:

Unmodified
Added
Removed
  • fact/tools/rootmacros/zerosearch.C

    r15119 r15368  
    1919    int VerbosityLevel
    2020){
    21 
    2221vector<Region> * ZeroPositions = new vector<Region>;
    23 Region currentRegion = {0, 0, 0, 0.0, 0, 0.0, 0.0, 0};
     22Region currentRegion = {0, 0, 0, 0.0, 0, 0.0, 0.0, 0, 0};
     23
    2424    if (step < 1){
    2525        if (VerbosityLevel > 0)
     
    6767            }
    6868        }
    69 
    7069    }
    7170    else // search for zero-x-ings on both esges
     
    8180        }
    8281    }
    83 
    8482
    8583    return ZeroPositions;
     
    172170                << reg->maxPos << " with Value:" << reg->maxVal << endl;
    173171        }
     172
     173        if (reg->maxPos <= 0)
     174        {
     175            reg = regions.erase( reg ) ;
     176            --reg;
     177            continue;
     178        }
    174179    }
    175180    return regions.size();
     
    411416*/
    412417
     418//size_t calcAmplitudeWeightedTime(
     419//        vector<Region>  &regions,
     420//        vector<float>   &data,
     421
     422//        ){
     423
     424//}
     425
     426size_t calcCFDPositions(
     427        vector<Region>  &regions,
     428        vector<float>   &cfd_data
     429        )
     430{
     431    vector<Region>::iterator reg;
     432    for (reg = regions.begin() ; reg < regions.end() ; reg++)
     433    {
     434        int     neg_zero_crossing = -1;
     435        int     min_in_cfd        = -1;
     436        int     pos_zero_crossing = reg->maxPos;
     437        float   minVal            = 10.0;
     438
     439        for (int i = pos_zero_crossing; i > pos_zero_crossing - 40; i--)
     440        {
     441            if (i <= 1 || i >= 1439) break;
     442
     443            if ( cfd_data[i] < minVal)
     444            {
     445                minVal      = cfd_data[i];
     446                min_in_cfd  = i;
     447            }
     448
     449            if (cfd_data[i] > 0 && cfd_data[i+1] <= 0 && cfd_data[i+2] < 0 )
     450            {
     451                neg_zero_crossing = i+1;
     452            }
     453        }
     454
     455        reg->cfdMinPos          = min_in_cfd;
     456        reg->cfdNegZerocrossing = neg_zero_crossing;
     457    }
     458    return regions.size();
     459}
     460
     461
    413462
    414463size_t findTimeOfHalfMaxLeft(
    415     vector<Region> &regions,
    416     vector<float> &data,
    417     float baseline,
    418     int beginRisingEdge,
    419     int endRisingEdge,
    420     int VerbosityLevel)
    421 {
    422     float pulse_size = 0;     //
    423     float thr = 0;
    424 
     464    vector<Region>  &regions,
     465    vector<float>   &data,
     466    float           baseline,
     467    int             beginRisingEdge,
     468    int             endRisingEdge,
     469    int             VerbosityLevel
     470        )
     471{
     472//    float   pulse_size      = 0;     //
     473    float   thr             = 0;
    425474    //int begin = beginRisingEdge;
    426     int end = endRisingEdge;
    427     int counter = 0;
     475    int     end             = 20;
     476    int     counter        = 0;
    428477    // local vars for line fitting
    429     float xmean = 0.0;
    430     float ymean = 0.0;
    431     float cov = 0.0;          // empiric covarianve between x and y
    432     float var = 0.0;          // empiric variance of x
     478    float   xmean          = 0.0;
     479    float   ymean          = 0.0;
     480    float   cov             = 0.0;          // empiric covariance between x and y
     481    float   var            = 0.0;          // empiric variance of x
    433482    // result of line fitting
    434     float slope = 0.0;
    435     float intercept = 0.0;
     483    float   slope          = 0.0;
     484    float   intercept      = 0.0;
    436485    // result of this function -- the position of the half rising edge
    437     float pos_of_thr_Xing = 0.0;
    438     int result = 0;
     486    float   pos_of_thr_Xing = 0.0;
     487    int     result          = 0;
     488
    439489
    440490    vector<Region>::iterator reg;
    441     for (reg = regions.begin() ; reg < regions.end() ; reg++){
    442         counter = 0;
     491    for (reg = regions.begin() ; reg < regions.end() ; reg++)
     492    {
    443493        // check if linear rising edge is completely in pipeline
    444494        // if not --> delete
     
    450500            continue;
    451501        }
    452 
    453         // The maximum was probably found using smoothed data,
    454         // so I read the value at the position of the maximum from the data
    455         // again. So I rely on a well defined maxPos.
    456         if (VerbosityLevel > 1) cout << "## baseline: " << baseline << endl;
    457         pulse_size = data[reg->maxPos] - baseline;
    458         if (VerbosityLevel > 1) cout << "## pulse_size: " << pulse_size << endl;
    459         thr = pulse_size / 2. + baseline;
    460         if (VerbosityLevel > 1) cout << "## thr: " << thr << endl;
    461502
    462503        // I think 5 slices left of the max there begins the rising edge
     
    469510        // ALARM check for out of range values...
    470511
     512        float average_of_max = -1;
     513        int stepwidth = 2;
     514
     515        counter = 0;
     516        for ( int i = reg->maxPos -stepwidth; i <= reg->maxPos + stepwidth; i++)
     517        {
     518            average_of_max += data[i];
     519            counter++;
     520        }
     521        average_of_max /= counter;
     522
     523        thr = average_of_max / 2.0;
     524
     525        float   max_slope       = 0;
     526        int     max_slope_pos   = -1;
     527        int     half_height_pos = -1;
     528
    471529        bool passed = false;
    472         for (int slice=reg->maxPos;
    473                     slice > 1; --slice)
    474         {
    475             if ( data[slice] < 0.9*data[reg->maxPos] && !passed)
    476             {
    477                 beginRisingEdge = reg->maxPos - slice;
    478                 passed = true;
    479             }
    480 
    481             if ( data[slice] < 0.1*data[reg->maxPos])
    482             {
    483                 endRisingEdge = reg->maxPos - slice;
    484                 passed = false;
    485                 reg->lengthOfRisingEdge=endRisingEdge - beginRisingEdge;
    486                 break;
    487             }
    488         }
    489 
    490 
     530        for (int slice=reg->maxPos; slice > reg->maxPos - end; --slice)
     531        {
     532            if (slice <= 0) break;
     533//            if ( data[slice] < 0.9*data[reg->maxPos] && !passed)
     534//            {
     535//                endRisingEdge = slice;
     536//                beginRisingEdge = endRisingEdge -3;
     537//                passed = true;
     538//            }
     539//            else if ( data[slice] < 0) //0.4*data[reg->maxPos] )
     540//            {
     541//                beginRisingEdge = slice;
     542//                reg->lengthOfRisingEdge=endRisingEdge - beginRisingEdge;
     543//                break;
     544//            }
     545
     546            //calculate max slope on leading edge
     547            if (data[slice] - data[slice-4] > max_slope && slice-4 > 0)
     548            {
     549                max_slope   = data[slice] - data[slice-4];
     550                max_slope_pos      = slice-2;
     551            }
     552
     553            if ( data[slice] <= thr && !passed){
     554                half_height_pos  = slice;
     555                passed  = true;
     556            }
     557        }
     558        reg->maxSlope           = max_slope/4;
     559        reg->maxSlopePos        = max_slope_pos;
     560        reg->halfRisingEdgePos  = half_height_pos;
     561        reg->halfHeight         = data[half_height_pos];
     562
     563        result                  = half_height_pos;
     564
     565        endRisingEdge   = result +2;
     566        beginRisingEdge = result -2;
     567
     568        counter = 0;
    491569        //calculate mean of x and y coordinate
    492         for (int slice=reg->maxPos - beginRisingEdge;
    493                     slice > reg->maxPos - endRisingEdge; --slice)
    494         {
     570        for (int slice= beginRisingEdge; slice < endRisingEdge; ++slice)
     571        {
     572            if (slice <= 0) break;
     573
    495574            xmean += slice;
    496575            ymean += data[slice];
    497576            counter++;
    498577        }
    499 //        xmean /= beginRisingEdge - endRisingEdge;
    500 //        ymean /= beginRisingEdge - endRisingEdge;
    501578        xmean /= counter;
    502579        ymean /= counter;
     
    508585        var = 0.0;
    509586
    510         for (int slice=reg->maxPos - beginRisingEdge;
    511                     slice > reg->maxPos - endRisingEdge; --slice)
    512         {
     587        //calculate covariance and variance
     588        for (int slice=beginRisingEdge; slice < endRisingEdge; ++slice)
     589        {
     590            if (slice <= 0) break;
     591
    513592            cov += (data[slice] - ymean) * (slice - xmean);
    514593            var += (slice - xmean) * (slice - xmean);
     
    517596        if (VerbosityLevel > 2) cout << "## var: " << var << endl;
    518597
    519         slope = cov / var;
    520         intercept = ymean - slope * xmean;
     598        slope       = cov / var;
     599        intercept   = ymean - slope * xmean;
     600
     601        // The maximum was probably found using smoothed data,
     602        // so I read the value at the position of the maximum from the data
     603        // again. So I rely on a well defined maxPos.
     604//        if (VerbosityLevel > 4) cout << "## baseline: " << baseline << endl;
     605//        pulse_size = data[reg->maxPos] - baseline;
     606//        if (VerbosityLevel > 4) cout << "## pulse_size: " << pulse_size << endl;
     607//        thr = pulse_size / 2. + baseline;
     608//        if (VerbosityLevel > 4) cout << "## thr: " << thr << endl;
    521609
    522610        // now calculate, where the fittet line crosses the threshold
    523         pos_of_thr_Xing = (thr - intercept) / slope;
    524         result = (int)(pos_of_thr_Xing + 0.5);
     611//        pos_of_thr_Xing = (thr - intercept) / slope;
     612//        result = (int)(pos_of_thr_Xing + 0.5);
    525613
    526614
     
    537625        }
    538626
    539         reg->halfRisingEdgePos = result;
    540 
    541         reg->distanceEdgeToMax = reg->maxPos - result;
    542         reg->interceptRisingEdge = intercept;
    543         reg->slopeOfRisingEdge = slope;
     627        reg->distanceEdgeToMax      = reg->maxPos - result;
     628        reg->interceptRisingEdge    = reg->maxPos - intercept;
     629        reg->slopeOfRisingEdge      = slope;
     630
     631        if (reg->halfRisingEdgePos <= 0 || reg->halfRisingEdgePos >= 1024)
     632        {
     633            reg = regions.erase( reg ) ;
     634            --reg;
     635            continue;
     636        }
    544637    }
    545638    return regions.size();
     
    570663size_t GetMaxPositions(
    571664    vector<Region>  &regions,
    572     vector<int>     &maxPositions,
    573     int VerbosityLevel)
    574 {
    575     maxPositions.clear();
     665    vector<int>     &positions,
     666    int VerbosityLevel)
     667{
     668    positions.clear();
    576669    vector<Region>::iterator it = regions.begin();
    577670    while( it != regions.end() )
    578671    {
    579         maxPositions.push_back(it->maxPos);
     672        positions.push_back(it->maxPos);
    580673        if (VerbosityLevel > 3){
    581674            cout << "getting maxPos@ " << it->maxPos << "\t";
     
    590683size_t GetEdgePositions(
    591684    vector<Region>  &regions,
    592     vector<int>     &edgePos,
    593     int VerbosityLevel)
    594 {
    595     maxPositions.clear();
     685    vector<int>     &positions,
     686    int VerbosityLevel)
     687{
     688    positions.clear();
    596689    vector<Region>::iterator it = regions.begin();
    597690    while( it != regions.end() )
    598691    {
    599         maxPositions.push_back(it->halfRisingEdgePos);
     692        positions.push_back(it->halfRisingEdgePos);
    600693        if (VerbosityLevel > 3){
    601694            cout << "getting maxPos@ " << it->halfRisingEdgePos << "\t";
Note: See TracChangeset for help on using the changeset viewer.