Ignore:
Timestamp:
04/24/12 11:02:49 (13 years ago)
Author:
Jens Buss
Message:
added function findTimeOfHalfMaxLeft
File:
1 edited

Legend:

Unmodified
Added
Removed

  • fact/tools/rootmacros/zerosearch.C

    r12682 r13424  
    33
    44// searches for zero crossings in a given vector of floats
    5 // for zero crossings in falling edge 
     5// for zero crossings in falling edge
    66// give edge = -1
    7 // 
    8 // 
     7//
     8//
    99// returns pointer to vetctor of ints
    1010
    1111vector<Region> *zerosearch(
    12         vector<float> &input,
    13         int edge,                       // search for transitions on rising edge=1, -1:falling
    14         unsigned int step,                      // search in steps of step
    15         int VerbosityLevel
     12    vector<float> &input,
     13    int edge,               // search on rising edge=1, -1:falling
     14    unsigned int step,      // search in steps of step
     15    int VerbosityLevel
    1616){
    1717
    1818vector<Region> * ZeroPositions = new vector<Region>;
    19 Region currentRegion = {0, 0, 0, 0.0};
    20         if (step < 1){
    21                 if (VerbosityLevel > 0)
    22                         cout << "zerosearch - stepsize=" << step
    23                                 << " is smaller than 1. returning." << endl;
    24                 return ZeroPositions;
    25         }
    26         if (input.size() < step){
    27                 if (VerbosityLevel > 0)
    28                         cout << "zerosearch - input-vector.size()="<< input.size()
    29                         << " is smaller than stepsize=" << step
    30                         << "returning." << endl;
    31                 return ZeroPositions;
    32         }
    33        
    34         if (edge > 0) // search for zero-x-ings on the rising edge
    35         {
    36                 for (unsigned int sl = 0 ; sl < (input.size()-step) ; sl+=step)
    37                 {
    38                         if (input[sl] > 0.0) // can't be a rising 0-x-ing
    39                         {
    40                                 continue;
    41                         }
    42                         if (input[sl] * input[sl+step] <= 0.0)  // 0-x-ing or both are 0
    43                         {
    44                                 currentRegion.begin     = sl;
    45                                 currentRegion.end       = sl+step;
    46                                 ZeroPositions->push_back(currentRegion);
    47                         }
    48                 }
    49         }
    50         else if (edge < 0) // search for zero-x-ings on the falling edge
    51         {
    52                 for (unsigned int sl = 0 ; sl < (input.size()-step) ; sl+=step)
    53                 {
    54                         if (input[sl] < 0.0) // can't be a falling 0-x-ing
    55                         {
    56                                 continue;
    57                         }
    58                         if (input[sl] * input[sl+step] <= 0.0)  // 0-x-ing or both are 0
    59                         {
    60                                 currentRegion.begin     = sl;
    61                                 currentRegion.end       = sl+step;
    62                                 ZeroPositions->push_back(currentRegion);
    63                         }
    64                 }
    65 
    66         }
    67         else // search for zero-x-ings on both esges
    68         {
    69                 for (unsigned int sl = 0 ; sl < (input.size()-step) ; sl+=step)
    70                 {
    71                         if (input[sl] * input[sl+step] <= 0.0)  // 0-x-ing or both are 0
    72                         {
    73                                 currentRegion.begin     = sl;
    74                                 currentRegion.end       = sl+step;
    75                                 ZeroPositions->push_back(currentRegion);
    76                         }
    77                 }
    78         }
    79 
    80        
    81         return ZeroPositions;
     19Region currentRegion = {0, 0, 0, 0.0, 0};
     20    if (step < 1){
     21        if (VerbosityLevel > 0)
     22            cout << "zerosearch - stepsize=" << step
     23                << " is smaller than 1. returning." << endl;
     24        return ZeroPositions;
     25    }
     26    if (input.size() < step){
     27        if (VerbosityLevel > 0)
     28            cout << "zerosearch - input-vector.size()="<< input.size()
     29            << " is smaller than stepsize=" << step
     30            << "returning." << endl;
     31        return ZeroPositions;
     32    }
     33
     34    if (edge > 0) // search for zero-x-ings on the rising edge
     35    {
     36        for (unsigned int sl = 0 ; sl < (input.size()-step) ; sl+=step)
     37        {
     38            if (input[sl] > 0.0) // can't be a rising 0-x-ing
     39            {
     40                continue;
     41            }
     42            if (input[sl] * input[sl+step] <= 0.0)  // 0-x-ing or both are 0
     43            {
     44                currentRegion.begin     = sl;
     45                currentRegion.end       = sl+step;
     46                ZeroPositions->push_back(currentRegion);
     47            }
     48        }
     49    }
     50    else if (edge < 0) // search for zero-x-ings on the falling edge
     51    {
     52        for (unsigned int sl = 0 ; sl < (input.size()-step) ; sl+=step)
     53        {
     54            if (input[sl] < 0.0) // can't be a falling 0-x-ing
     55            {
     56                continue;
     57            }
     58            if (input[sl] * input[sl+step] <= 0.0)  // 0-x-ing or both are 0
     59            {
     60                currentRegion.begin     = sl;
     61                currentRegion.end       = sl+step;
     62                ZeroPositions->push_back(currentRegion);
     63            }
     64        }
     65
     66    }
     67    else // search for zero-x-ings on both esges
     68    {
     69        for (unsigned int sl = 0 ; sl < (input.size()-step) ; sl+=step)
     70        {
     71            if (input[sl] * input[sl+step] <= 0.0)  // 0-x-ing or both are 0
     72            {
     73                currentRegion.begin     = sl;
     74                currentRegion.end       = sl+step;
     75                ZeroPositions->push_back(currentRegion);
     76            }
     77        }
     78    }
     79
     80
     81    return ZeroPositions;
    8282}
    8383
    8484size_t ShiftRegionBy(vector<Region> &src,
    85         int Shift,
    86         int VerbosityLevel)
    87 {
    88 vector<Region>::iterator it;
    89         // I copy code here ... not good, but I hope nobody is giving VL > 10 ever...
    90         if (VerbosityLevel > 10){
    91                 for (it = src.begin() ; it < src.end() ; it++){
    92                         cout << " it->begin="<< it->begin;
    93                         it->begin += Shift;
    94                         cout << "--> shifted="<< it->begin << endl;
    95                         cout << " it->end="<< it->end;
    96                         it->end += Shift;
    97                         cout << "--> shifted="<< it->end << endl;
    98                         cout << " it->maxPos="<< it->maxPos;
    99                         it->maxPos += Shift;
    100                         cout << "--> shifted="<< it->maxPos << endl;
    101                 }
    102                 return src.size();
    103         }
    104 
    105         for (it = src.begin() ; it < src.end() ; it++){
    106                 it->begin += Shift;
    107                 it->end += Shift;
    108                 it->maxPos += Shift;
    109         }
    110         return src.size();
     85    int Shift,
     86    int VerbosityLevel)
     87{
     88    vector<Region>::iterator it;
     89    // I copy code here ... not good, but I hope nobody is giving VL > 10 ever...
     90    if (VerbosityLevel > 10){
     91        for (it = src.begin() ; it < src.end() ; it++){
     92            cout << " it->begin="<< it->begin;
     93            it->begin += Shift;
     94            cout << "--> shifted="<< it->begin << endl;
     95            cout << " it->end="<< it->end;
     96            it->end += Shift;
     97            cout << "--> shifted="<< it->end << endl;
     98            cout << " it->maxPos="<< it->maxPos;
     99            it->maxPos += Shift;
     100            cout << "--> shifted="<< it->maxPos << endl;
     101        }
     102    return src.size();
     103    }
     104
     105    for (it = src.begin() ; it < src.end() ; it++){
     106        it->begin += Shift;
     107        it->end += Shift;
     108        it->maxPos += Shift;
     109    }
     110    return src.size();
    111111}
    112112
    113113size_t EnlargeRegion(vector<Region> &src,
    114         int Left,
    115         int Right,
    116         int VerbosityLevel)
    117 {
    118 vector<Region>::iterator it;
    119         // I copy code here ... not good, but I hope nobody is giving VL > 10 ever...
    120         if (VerbosityLevel > 10){
    121                 for (it = src.begin() ; it < src.end() ; it++){
    122                         cout << " it->begin="<< it->begin;
    123                         it->begin -= Left;
    124                         cout << "--> enlarged="<< it->begin << endl;
    125                         cout << " it->end="<< it->end;
    126                         it->end += Right;
    127                         cout << "--> enlarged="<< it->end << endl;
    128                 }
    129                 return src.size();
    130         }
    131 
    132 
    133                 for (it = src.begin() ; it < src.end() ; it++){
    134                 it->begin -= Left;
    135                 it->end += Right;
    136         }
    137         return src.size();
     114    int Left,
     115    int Right,
     116    int VerbosityLevel)
     117{
     118    vector<Region>::iterator it;
     119    // I copy code here ... not good, but I hope nobody is giving VL > 10 ever...
     120    if (VerbosityLevel > 10){
     121        for (it = src.begin() ; it < src.end() ; it++){
     122            cout << " it->begin="<< it->begin;
     123            it->begin -= Left;
     124            cout << "--> enlarged="<< it->begin << endl;
     125            cout << " it->end="<< it->end;
     126            it->end += Right;
     127            cout << "--> enlarged="<< it->end << endl;
     128        }
     129        return src.size();
     130    }
     131
     132
     133    for (it = src.begin() ; it < src.end() ; it++){
     134        it->begin -= Left;
     135        it->end += Right;
     136    }
     137    return src.size();
    138138
    139139}
    … …  
    141141#include <limits>
    142142size_t findAbsMaxInRegions(
    143         vector<Region> &regions,
    144         vector<float> &data,
    145         int VerbosityLevel)
    146 {
    147 vector<Region>::iterator reg;
    148         for (reg = regions.begin() ; reg < regions.end() ; reg++){
    149                 reg->maxVal=-numeric_limits<float>::max();
    150                 // 1st check if both ends of the region are in the data-vector
    151                 if (reg->begin < 0) reg->begin=0;
    152                 if ((unsigned int)reg->end > data.size()-1) reg->end=data.size()-1;
    153 
    154                 // TODO or TOTHINKOF:
    155                 // Region might overlap ... I don't care .. I just search
    156                 // in a later step Maxima at the same Positions can be removed...
    157                 // of course this means, some searches were needless ...
    158                 //cout << "searching from " << reg->begin << " to " << reg->end << endl;
    159                 for (unsigned int pos=reg->begin; pos <= (unsigned int)reg->end; pos++){
    160                         if (data[pos] > reg->maxVal){
    161                                 reg->maxVal = data[pos];
    162                                 reg->maxPos = pos;
    163                         }
    164                 }
    165                 if (VerbosityLevel > 3) {
    166                         cout << "findAbsMaxInRegions - found Max at:"
    167                                 << reg->maxPos << " with Value:" << reg->maxVal << endl;
    168                 }
    169         }
    170         return regions.size();
     143    vector<Region> &regions,
     144    vector<float> &data,
     145    int VerbosityLevel)
     146{
     147    vector<Region>::iterator reg;
     148    for (reg = regions.begin() ; reg < regions.end() ; reg++){
     149        reg->maxVal=-numeric_limits<float>::max();
     150        // 1st check if both ends of the region are in the data-vector
     151        if (reg->begin < 0) reg->begin=0;
     152        if ((unsigned int)reg->end > data.size()-1) reg->end=data.size()-1;
     153
     154        // TODO or TOTHINKOF:
     155        // Region might overlap ... I don't care .. I just search
     156        // in a later step Maxima at the same Positions can be removed...
     157        // of course this means, some searches were needless ...
     158        //cout << "searching from " << reg->begin << " to " << reg->end << endl;
     159        for (unsigned int pos=reg->begin; pos <= (unsigned int)reg->end; pos++){
     160            if (data[pos] > reg->maxVal){
     161                reg->maxVal = data[pos];
     162                reg->maxPos = pos;
     163            }
     164        }
     165        if (VerbosityLevel > 3) {
     166            cout << "findAbsMaxInRegions - found Max at:"
     167                << reg->maxPos << " with Value:" << reg->maxVal << endl;
     168        }
     169    }
     170    return regions.size();
    171171}
    172172
    … …  
    177177
    178178size_t removeEqualMaxima(
    179         vector<Region> &regions,
    180         int VerbosityLevel)
    181 {
    182         if (VerbosityLevel > 2){
    183                 cout << "removeEqualMaxima:" << endl;
    184                 cout << "region before contains:" << endl;
    185                 for (unsigned int i=0; i<regions.size(); i++){
    186                         cout << i <<"\t";
    187                         cout << regions[i].begin <<"\t";
    188                         cout << regions[i].end <<"\t";
    189                         cout << regions[i].maxPos <<"\t";
    190                         cout << regions[i].maxVal <<"\t";
    191                         cout << endl;
    192                 }
    193         }
    194 
    195         vector<Region>::iterator it;
    196         it = unique (regions.begin(), regions.end() , compMaxPosOfRegions);
    197         regions.resize( it - regions.begin() );
    198 
    199         if (VerbosityLevel > 1){
    200                 cout << "region now contains:" << endl;
    201                 for (unsigned int i=0; i<regions.size(); i++){
    202                         cout << i <<"\t";
    203                         cout << regions[i].begin <<"\t";
    204                         cout << regions[i].end <<"\t";
    205                         cout << regions[i].maxPos <<"\t";
    206                         cout << regions[i].maxVal <<"\t";
    207                         cout << endl;
    208                 }
    209         }
    210         return regions.size();
     179    vector<Region> &regions,
     180    int VerbosityLevel)
     181{
     182    if (VerbosityLevel > 2){
     183        cout << "removeEqualMaxima:" << endl;
     184        cout << "region before contains:" << endl;
     185        for (unsigned int i=0; i<regions.size(); i++){
     186            cout << i <<"\t";
     187            cout << regions[i].begin <<"\t";
     188            cout << regions[i].end <<"\t";
     189            cout << regions[i].maxPos <<"\t";
     190            cout << regions[i].maxVal <<"\t";
     191            cout << endl;
     192        }
     193    }
     194
     195    vector<Region>::iterator it;
     196    it = unique (regions.begin(), regions.end() , compMaxPosOfRegions);
     197    regions.resize( it - regions.begin() );
     198
     199    if (VerbosityLevel > 1){
     200        cout << "region now contains:" << endl;
     201        for (unsigned int i=0; i<regions.size(); i++){
     202            cout << i <<"\t";
     203            cout << regions[i].begin <<"\t";
     204            cout << regions[i].end <<"\t";
     205            cout << regions[i].maxPos <<"\t";
     206            cout << regions[i].maxVal <<"\t";
     207            cout << endl;
     208        }
     209    }
     210    return regions.size();
    211211}
    212212
    213213size_t removeRegionOnFallingEdge(
    214         vector<Region> &regions,
    215         unsigned int FallingEdgeWidth,
    216         int VerbosityLevel)
     214    vector<Region> &regions,
     215    unsigned int FallingEdgeWidth,
     216    int VerbosityLevel)
    217217{
    218218if ( FallingEdgeWidth < 1){
    … …  
    224224    }
    225225  }
    226         if (regions.size() < 1){
    227                 if (VerbosityLevel > 3)
    228                         cout << "removeRegionOnFallingEdge: source vector empty." << endl;
    229                 return 0;
    230         }
     226    if (regions.size() < 1){
     227        if (VerbosityLevel > 3)
     228            cout << "removeRegionOnFallingEdge: source vector empty." << endl;
     229        return 0;
     230    }
    231231
    232232 vector<Region>::reverse_iterator it = regions.rbegin();
    … …  
    234234  {
    235235  if (VerbosityLevel >10)
    236                 cout << it->maxPos << "\t" << (it+1)->maxPos << "\t" << it->maxPos - (it+1)->maxPos << endl;
     236        cout << it->maxPos << "\t" << (it+1)->maxPos << "\t" << it->maxPos - (it+1)->maxPos << endl;
    237237
    238238    if ( it->maxPos - (it+1)->maxPos < (int)FallingEdgeWidth) {
    … …  
    250250
    251251size_t removeRegionWithMaxOnEdge(
    252         vector<Region> &regions,
    253         unsigned int EdgeWidth,
    254         int VerbosityLevel)
    255 {
    256         if (EdgeWidth < 1){
    257                 if (VerbosityLevel > 0){
    258                         cout << "removeReginWithMaximaOnEdge: EdgeWidth < 1" << endl;
    259                         cout << "EdgeWidth=" << EdgeWidth << endl;
    260                         cout << "returning." << endl;
    261                 return regions.size();
    262                 }
    263         }
    264 
    265         vector<Region>::iterator it = regions.begin();
    266         while( it != regions.end() )
    267         {
     252    vector<Region> &regions,
     253    unsigned int EdgeWidth,
     254    int VerbosityLevel)
     255{
     256    if (EdgeWidth < 1){
     257        if (VerbosityLevel > 0){
     258            cout << "removeReginWithMaximaOnEdge: EdgeWidth < 1" << endl;
     259            cout << "EdgeWidth=" << EdgeWidth << endl;
     260            cout << "returning." << endl;
     261        return regions.size();
     262        }
     263    }
     264
     265    vector<Region>::iterator it = regions.begin();
     266    while( it != regions.end() )
     267    {
    268268//      cout << it->begin << "\t";
    269269//      cout << it->maxPos << "\t";
    … …  
    271271//      cout << it->maxVal << endl;
    272272
    273                 if (it->maxPos < (int)(it->begin+EdgeWidth)) {
    274                         if (VerbosityLevel > 3)
    275                                 cout << "erasing Region(max@left edge) " << it->maxPos << endl;
    276                         it = regions.erase( it ) ;
    277                         //++it;
    278                 }else if (it->maxPos > (int)(it->end-EdgeWidth)) {
    279                         if (VerbosityLevel > 3)
    280                                 cout << "erasing Region(max@right edge) " << it->maxPos << endl;
    281                         it = regions.erase( it ) ;
    282                         //++it;
    283                 }else
    284                         ++it;
    285 
    286         }
    287 
    288         return regions.size();
     273        if (it->maxPos < (int)(it->begin+EdgeWidth)) {
     274            if (VerbosityLevel > 3)
     275                cout << "erasing Region(max@left edge) " << it->maxPos << endl;
     276            it = regions.erase( it ) ;
     277            //++it;
     278        }else if (it->maxPos > (int)(it->end-EdgeWidth)) {
     279            if (VerbosityLevel > 3)
     280                cout << "erasing Region(max@right edge) " << it->maxPos << endl;
     281            it = regions.erase( it ) ;
     282            //++it;
     283        }else
     284            ++it;
     285
     286    }
     287
     288    return regions.size();
    289289}
    290290
    291291size_t removeMaximaBelow(
    292         vector<Region> &regions,
    293         float threshold,
    294         int VerbosityLevel)
     292    vector<Region> &regions,
     293    float threshold,
     294    int VerbosityLevel)
    295295{
    296296//      if (threshold < 0){
    … …  
    302302//      }
    303303
    304         vector<Region>::iterator it = regions.begin();
    305         while( it != regions.end() )
    306         {
    307                 if (it->maxVal < threshold ) {
    308                         if (VerbosityLevel > 3){
    309                                 cout << "removing max " << it->maxVal << "\t";
    310                                 cout << " @ " << it->maxPos << "\t";
    311                                 cout << endl;
    312                         }
    313                         it = regions.erase( it ) ;
    314                 }else
    315                         ++it;
    316         }
    317 
    318         return regions.size();
     304    vector<Region>::iterator it = regions.begin();
     305    while( it != regions.end() )
     306    {
     307        if (it->maxVal < threshold ) {
     308            if (VerbosityLevel > 3){
     309                cout << "removing max " << it->maxVal << "\t";
     310                cout << " @ " << it->maxPos << "\t";
     311                cout << endl;
     312            }
     313            it = regions.erase( it ) ;
     314        }else
     315            ++it;
     316    }
     317
     318    return regions.size();
    319319}
    320320
    321321size_t removeMaximaAbove(
    322         vector<Region> &regions,
    323         float threshold,
    324         int VerbosityLevel)
    325 {
    326         if (threshold < 0){
    327                 if (VerbosityLevel > 0)
    328                         cout << "removeMaximaAbove: threshold < 0" << endl;
    329                         cout << "threshold=" << threshold << endl;
    330                         cout << "returning." << endl;
    331                 return regions.size();
    332         }
    333 
    334         vector<Region>::iterator it = regions.begin();
    335         while( it != regions.end() )
    336         {
    337                 if (it->maxVal > threshold ) {
    338                         if (VerbosityLevel > 3){
    339                                 cout << "removing max " << it->maxVal << "\t";
    340                                 cout << " @ " << it->maxPos << "\t";
    341                                 cout << endl;
    342                         }
    343                         it = regions.erase( it ) ;
    344                 }else
    345                         ++it;
    346         }
    347 
    348         return regions.size();
     322    vector<Region> &regions,
     323    float threshold,
     324    int VerbosityLevel)
     325{
     326    if (threshold < 0){
     327        if (VerbosityLevel > 0)
     328            cout << "removeMaximaAbove: threshold < 0" << endl;
     329            cout << "threshold=" << threshold << endl;
     330            cout << "returning." << endl;
     331        return regions.size();
     332    }
     333
     334    vector<Region>::iterator it = regions.begin();
     335    while( it != regions.end() )
     336    {
     337        if (it->maxVal > threshold ) {
     338            if (VerbosityLevel > 3){
     339                cout << "removing max " << it->maxVal << "\t";
     340                cout << " @ " << it->maxPos << "\t";
     341                cout << endl;
     342            }
     343            it = regions.erase( it ) ;
     344        }else
     345            ++it;
     346    }
     347
     348    return regions.size();
    349349}
    350350
    351351Region FindAbsMax(
    352         vector<Region> &regions,
    353         int VerbosityLevel)
    354 {
    355         Region AbsMax;
    356         AbsMax.begin = -1;
    357         AbsMax.maxPos = -1;
    358         AbsMax.end = -1;
    359         AbsMax.maxVal=-numeric_limits<float>::max();
    360         if (regions.size() < 1)
    361                 return AbsMax;
    362 
    363         for (vector<Region>::iterator it = regions.begin(); it < regions.end(); ++it)
    364         {
    365                 if (it->maxVal > AbsMax.maxVal ) {
    366                         AbsMax = *it;
    367                 }
    368         }
    369         if (VerbosityLevel > 5){
    370                 AbsMax.begin +=0;
    371         }
    372         return AbsMax;
     352    vector<Region> &regions,
     353    int VerbosityLevel)
     354{
     355    Region AbsMax;
     356    AbsMax.begin = -1;
     357    AbsMax.maxPos = -1;
     358    AbsMax.end = -1;
     359    AbsMax.maxVal=-numeric_limits<float>::max();
     360    if (regions.size() < 1)
     361        return AbsMax;
     362
     363    for (vector<Region>::iterator it = regions.begin(); it < regions.end(); ++it)
     364    {
     365        if (it->maxVal > AbsMax.maxVal ) {
     366            AbsMax = *it;
     367        }
     368    }
     369    if (VerbosityLevel > 5){
     370        AbsMax.begin +=0;
     371    }
     372    return AbsMax;
    373373}
    374374
    … …  
    376376////////////// old Version of the code
    377377/*
    378         for (unsigned int sl = step; sl < input.size(); sl+=step){
    379                 if (input[sl] * input[sl-] < 0 || input[sl]==0.0 ){ // sign change --> zero crossing OR really zero
    380 
    381                         if ( (input[sl-1] - input[sl]) * edge < 0){ // this is the crossing the user wanted
    382 
    383                                 // check if we go lower than a certain limit in the next few slices
    384                                 for ( int lala=0; lala<post; lala++){
    385                                         if (input[sl+lala] > 1.5) {
    386                                                 zeroPositions->push_back(sl);
    387                                                 sl += lala;
    388                                                 break;
    389                                         }
    390                                 }
    391 
    392                         } else if ( (input[sl-1] - input[sl]) * edge > 0){ // this is the zero x-ing the user did not want
    393 
    394                                 // do nothing
    395 
    396                         } else { // sl and sl-1 are equal .. don't know waht do to...
    397 
    398                         }
    399 
    400                 }
    401 
    402         } // end of loop over slices - between pre and post
     378    for (unsigned int sl = step; sl < input.size(); sl+=step){
     379        if (input[sl] * input[sl-] < 0 || input[sl]==0.0 ){ // sign change --> zero crossing OR really zero
     380
     381            if ( (input[sl-1] - input[sl]) * edge < 0){ // this is the crossing the user wanted
     382
     383                // check if we go lower than a certain limit in the next few slices
     384                for ( int lala=0; lala<post; lala++){
     385                    if (input[sl+lala] > 1.5) {
     386                        zeroPositions->push_back(sl);
     387                        sl += lala;
     388                        break;
     389                    }
     390                }
     391
     392            } else if ( (input[sl-1] - input[sl]) * edge > 0){ // this is the zero x-ing the user did not want
     393
     394                // do nothing
     395
     396            } else { // sl and sl-1 are equal .. don't know waht do to...
     397
     398            }
     399
     400        }
     401
     402    } // end of loop over slices - between pre and post
    403403
    404404*/
     405
     406
     407size_t findTimeOfHalfMaxLeft(
     408    vector<Region> &regions,
     409    vector<float> &data,
     410    float baseline,
     411    int beginRisingEdge,
     412    int endRisingEdge,
     413    int VerbosityLevel)
     414{
     415    float pulse_size = 0;     //
     416    float thr = 0;
     417
     418    //int begin = beginRisingEdge;
     419    int end = endRisingEdge;
     420    int counter = 0;
     421    // local vars for line fitting
     422    float xmean = 0.0;
     423    float ymean = 0.0;
     424    float cov = 0.0;          // empiric covarianve between x and y
     425    float var = 0.0;          // empiric variance of x
     426    // result of line fitting
     427    float slope = 0.0;
     428    float intercept = 0.0;
     429    // result of this function -- the position of the half rising edge
     430    float pos_of_thr_Xing = 0.0;
     431    int result = 0;
     432
     433    vector<Region>::iterator reg;
     434    for (reg = regions.begin() ; reg < regions.end() ; reg++){
     435
     436        // check if linear falling edge is completely in pipeline
     437        // if not --> delete
     438        if ( reg->maxPos - end < 0 )
     439        {
     440            // delete this region from vector<Region>
     441            reg = regions.erase( reg ) ;
     442            --reg;
     443            continue;
     444        }
     445
     446        // The maximum was probably found using smoothed data,
     447        // so I read the value at the position of the maximum from the data
     448        // again. So I rely on a well defined maxPos.
     449        if (VerbosityLevel > 1) cout << "## baseline: " << baseline << endl;
     450        pulse_size = data[reg->maxPos] - baseline;
     451        if (VerbosityLevel > 1) cout << "## pulse_size: " << pulse_size << endl;
     452        thr = pulse_size / 2. + baseline;
     453        if (VerbosityLevel > 1) cout << "## thr: " << thr << endl;
     454
     455        // I think 5 slices left of the max there begins the rising edge
     456        // and I think the rising edge is about 6 slices long.
     457        // but these numbers are input as parameters
     458        // beginRisingEdge &
     459        // endRisingEdge
     460
     461        // I fit a linear function to the falling edge
     462        // ALARM check for out of range values...
     463
     464
     465        for (int slice=reg->maxPos - beginRisingEdge;
     466                    slice > reg->maxPos - endRisingEdge; --slice)
     467        {
     468            xmean += slice;
     469            ymean += data[slice];
     470            counter++;
     471        }
     472//        xmean /= beginRisingEdge - endRisingEdge;
     473//        ymean /= beginRisingEdge - endRisingEdge;
     474        xmean /= counter;
     475        ymean /= counter;
     476        if (VerbosityLevel > 1) cout << "## xmean: " << xmean << endl;
     477        if (VerbosityLevel > 1) cout << "## ymean: " << ymean << endl;
     478
     479        for (int slice=reg->maxPos - beginRisingEdge;
     480                    slice > reg->maxPos - endRisingEdge; --slice)
     481        {
     482            cov = (data[slice] - ymean) * (slice - xmean);
     483            var = (slice - xmean) * (slice - xmean);
     484        }
     485        if (VerbosityLevel > 1) cout << "## cov: " << cov << endl;
     486        if (VerbosityLevel > 1) cout << "## var: " << var << endl;
     487
     488        slope = cov / var;
     489        intercept = ymean - slope * xmean;
     490        // now calculate, where the fittet line crosses the threshold
     491        pos_of_thr_Xing = (thr - intercept) / slope;
     492        result = (int)(pos_of_thr_Xing + 0.5);
     493
     494
     495        if (VerbosityLevel > 0)
     496        {
     497            cout << "findTimeOfHalfMaxLeft() is still in debugging phase:" << endl;
     498            cout << "please edit the code in oder to suppress this output:" << endl;
     499            cout << "threshold: " << thr << endl;
     500            cout << "slope: " << slope << " [mV/timeslice]" << endl;
     501            cout << "intercept: " << intercept << "[mV]" << endl;
     502            cout << "position where line crosses threshold " << pos_of_thr_Xing << endl;
     503            cout << "converted to slice number " << result << endl;
     504
     505        }
     506
     507        reg->halfRisingEdgePos = result;
     508        reg->slopeOfRisingEdge = slope;
     509    }
     510    return regions.size();
     511}
Note: See TracChangeset for help on using the changeset viewer.