Series->Histogram for SPHERES

pub/src/commands.cpp

@@ -526,6 +526,7 @@ bool fridaCommand( string cmd )
             "  ryc  raw counts (no normalization)\n"
             "  rs   load series of energy spectra\n"
             "  rs2  load also spectra for open state\n"
+            "  rsh  load series and make count histograms\n"
             ;
 
     } else if (cmd == "ry" ) {
@@ -546,6 +547,8 @@ bool fridaCommand( string cmd )
         NRSSM::ReadSeries( 0 );
     } else if (cmd == "rs2") {
         NRSSM::ReadSeries( 1 );
+    } else if (cmd == "rsh") {
+        NRSSM::ReadSeries2Histo( 0 );
 
 
     } else if (cmd == "!") {

pub/src/rssm.cpp

@@ -212,6 +212,8 @@ void CRawfileSpheres::RdRawYam( ifstream& F_in )
     }
 }    
 
+//! Read one raw data file.
+
 void NRSSM::ReadSpec( int flag )
 
 // flag values (to be OR'ed):
@@ -368,6 +370,7 @@ void NRSSM::ReadSpec( int flag )
 
 }
 
+//! Read a series of raw data files.
 
 void NRSSM::ReadSeries( int flag )
 
@@ -564,3 +567,143 @@ void NRSSM::ReadSeries( int flag )
     }
 
 }
+
+//! Read a series of raw data files to produce count histograms.
+
+void NRSSM::ReadSeries2Histo( int flag )
+{
+    char tstrg[30];
+    string fser, fnam, name;
+    time_t t;
+    uint isub;
+    vector<CRawfileSpheres> RR;
+
+    fser = wask("Read SPHERES data from series");
+    static uint isub0 = 0;
+    isub0 = iask( "Start with subscan", isub0 );
+    static uint mE = 32;
+    mE = iask( "E binning factor", mE);
+    static uint mR = 8;
+    mR = iask( "Raw binning factor", mR);
+
+    // Read consolidated files:
+    for( isub=isub0; ; ++isub ){
+        fnam = fser+"a"+S(isub);
+        ifstream F_in;
+        F_in.open(fnam.c_str(),fstream::in);
+        if( !(F_in.is_open()) )
+            break;
+        // cout << "successfully opened "<<fnam<<"\n";
+        CRawfileSpheres rf; 
+        rf.RdRawYam( F_in );
+        RR.push_back( rf );
+    }
+    int nsub = RR.size();
+    if( nsub==0 )
+        throw "could not open file " + fnam;
+    printf( "successfully read %d files\n", nsub );
+    nsub = nsub/mR*mR;
+
+    // Correct time:
+    double mean_time =
+        ( RR[nsub-1].measured_until - RR[0].measured_until ) / nsub;
+    printf( "mean measuring time %g seconds per file\n", mean_time );
+    t = (time_t) (RR[0].measured_until-mean_time);
+    strftime( tstrg, 30, "%F %b %H:%M:%S", localtime( &t ) );
+    printf( "measurement times relative to start time %s\n", tstrg );
+    for( isub=0; isub<nsub; ++isub )
+        RR[isub].measured_at = RR[isub].measured_until - 
+            RR[0].measured_until + (int) (mean_time/2);
+
+    // Get parameters, check consistency:
+    int ndet = RR[0].ndet;
+    printf( "used %i detectors\n", ndet );
+    double tstep = RR[0].daq_time_step;
+    bool incremental_in = RR[0].incremental;
+
+    for( isub=1; isub<nsub; ++isub ){
+        if( RR[isub].ndet!=ndet )
+            throw "inconsistent no of det: file " + S(isub)
+                + " has " + S(RR[isub].ndet) + " instead of " + S(ndet);
+        if( RR[isub].daq_time_step!=tstep )
+            throw "inconsistent time step";
+        if( RR[isub].incremental!=incremental_in )
+            throw "inconsistent increment status";
+    }
+
+    // Set file header:
+    POld olf( new COld );
+    olf->lDoc.push_back( "acquire histograms from series"+fser );
+    olf->xco = CCoord( "cts", "" );
+    olf->yco = CCoord( "n", "" );
+    olf->ZCo.push_back( CCoord( "det", "") );
+    olf->ZCo.push_back( CCoord( "E", "ueV") );
+    triv::fname_divide( fser, 0, &name, 0 );
+    olf->name = name;
+
+    if( RR[0].maj_outform<3 )
+        throw "old format not supported";
+
+    uint nE = RR[0].rawdata[0].size()/(1+ndet*2);
+    nE = nE/mE*mE;
+    vector<double> *raw;
+    raw = RR[0].rawdata;
+    uint j, iE;
+
+    for( j=0; j<ndet; ++j ){
+        for( iE=0; iE<nE; iE += mE ){
+            PSpec s( new CSpec );
+            s->z.resize( 2 );
+            s->z[0] = j;
+            s->z[1] = ( raw[0][iE*(1+ndet*2)]+raw[0][(iE+mE-1)*(1+ndet*2)] )/2;
+            s->resize( 24, false );
+            olf->V.push_back( s );
+        }
+    }
+
+    // Reverse incrementation?
+    if( incremental_in ){
+        printf( "reversing incrementation\n" );
+        for( isub=nsub-1; isub>=1; --isub ){
+            for( uint iarr=0; iarr<4; ++iarr ){
+                for( iE=0; iE<nE; ++iE ){
+                    for( j=1; j<1+ndet*2; ++j ){
+                        RR[isub].rawdata[iarr][iE*(1+ndet*2)+j] -= 
+                            RR[isub-1].rawdata[iarr][iE*(1+ndet*2)+j];
+                    }
+                }
+            }
+        }
+    }
+
+    // Transfer data from raw arrays, sum over halfspaces:
+    double count;
+    vector<double> x( 24 );
+    for ( uint cou= 0; cou<24; ++cou )
+        x[cou] = cou;
+
+    for( j=0; j<ndet; ++j ){
+        for( iE=0; iE<nE; iE+=mE ){
+            vector<double> y( 24, 0. );
+            for( isub=0; isub<nsub; isub += mR ){
+                for( uint iiE=0; iiE<mE; ++iiE ){
+                    raw = RR[isub].rawdata;
+                    uint iraw = (iE+iiE)*(1+ndet*2);
+                    count = raw[0][iraw+1+j] + raw[1][iraw+1+j];
+                    for( uint iiR=1; iiR<mR; ++iiR ){
+                        raw = RR[isub+iiR].rawdata;
+                        count += raw[0][iraw+1+j] + raw[1][iraw+1+j];
+                    }
+                    uint cou = (uint) (count+0.5);
+                    if( cou>= 24 )
+                        cou = 24-1;
+                    ++y[cou];
+                }
+            }
+            olf->VS( (j*nE+iE)/mE )->x = x;
+            olf->VS( (j*nE+iE)/mE )->y = y;
+        }
+    }
+    
+    NOlm::mem_store( olf );
+}

pub/src/rssm.h

@@ -12,4 +12,5 @@
 namespace NRSSM {
     void ReadSpec( int flag );
     void ReadSeries( int flag );
+    void ReadSeries2Histo( int flag );
 };