BeamSummary.cxx

Go to the documentation of this file.

#include <iostream>
#include <cstdlib>
#include <cmath>

#include "BeamData/ana/Summary/BeamSummary.h"
#include "MessageService/MsgService.h"

#include "TTree.h"

CVSID("$Id: BeamSummary.cxx,v 1.7 2006/05/27 07:24:02 rhatcher Exp $");

using namespace std;

BeamSummary::BeamSummary(){
    clear();
    fChain=0;
}

BeamSummary::BeamSummary(TTree* btree)
{
    clear();
    SetTree(btree);
    fChain->SetBranchAddress("Beaminfo",this);
}

void BeamSummary::SetTree(TTree* btree)
{
    fChain = btree;
}

void BeamSummary::clear()
{
    beamIntensity = 0.0;    
    hBeamWidth = vBeamWidth = 0.0;       
    hPosatTargetPM = vPosatTargetPM = 0.0;
    hPosatTargetBPM = vPosatTargetBPM = 0.0;  
    tanHoriz = tanVert = 0.0;          
    hornPeakCurrent = nuTarZ = 0.0;           
    hadTotalPH = hadMeanX = hadMeanY = 0.0;         
    mu1TotalPH = mu1MeanX = mu1MeanY = 0.0;         
    mu2TotalPH = mu2MeanX = mu2MeanY = 0.0;         
    mu3TotalPH = mu3MeanX = mu3MeanY = 0.0;         
    timeStampD = 0.0;       
    timeStampSec = 0;
}

Int_t BeamSummary::MatchSpillWithTime(Int_t tsec, Int_t tnanosec=0)
{
    static Long64_t lastindex=0;
    Long64_t index=0;
    Int_t found=0;

   // Get the total number of entries in the ntuple
    Long64_t nentries = fChain->GetEntries();
    if (nentries<1) {
        MSG("BeamSummary", Msg::kWarning) << "No entries in the beam ntuple" << endl;
    }
    // First set the index to the position of the previous call
    index = lastindex;
    // Check if this index is not larger than the number of
    // entries in the ntuple. If so, put it to the last entry in the
    // tree
    while (index>=nentries) {
        --index;
    }        
    // the search process can put the lastindex to a negative value
    // if so, put the index to 0
    if (index<0) index=0;

    // The ntuple is supposed to be order in time, so it is assumed here...
    // Read in the previous entry and determine whether this occured earlier
    // in time, if not change position to the left, else proceed to the right.
    
    clear();
    Long64_t ientry = fChain->LoadTree(index);
    if (ientry < 0) {
        MSG("BeamSummary",Msg::kWarning) << "Error reading beam ntuple" << endl;
        return found;
    }
    Int_t nb = fChain->GetEntry(index);
    if (nb<0) {
        MSG("BeamSummary",Msg::kWarning) << "No byte read in" << endl;
        return found;
    }
    Int_t timedif = tsec-timeStampSec;
    if (abs(timedif)<2) {
        // hey, we are already at the right spot
        found =1;
    }
    else if (timedif<=-2) {
        // time to look up is earlier than current entry, proceed left
        index--;
        for (Long64_t ient=index; ient>=0; ient--) {
            index = ient;
            fChain->GetEntry(index);
            timedif = tsec-timeStampSec;
            if (abs(timedif)<2) {
                found = 1;
                break;
            }
        }
    }
    else if (timedif>=2) {
        // time to look up is later than current entry, proceed right
        index++;
        for (Long64_t ient=index; ient<nentries; ient++) {
            index = ient;
            fChain->GetEntry(index);
            timedif = tsec-timeStampSec;
            if (abs(timedif)<2) {
                found = 1;
                break;
            } 
        }
    }
    else {
        MSG("BeamSummary",Msg::kDebug) << "This is not supposed to happen!" << endl;
    }
    
    if (found==0) {
        clear();
        MSG("BeamSummary",Msg::kWarning)
            << "No entry found in Beam ntuple corresponding to time stamp:"
            << " " << tsec << " seconds " << tnanosec << " nanoseconds" << endl; 
    } else {
        // Matching with seconds does not guarantee it is the closest
        // spill. This can happen when the repitition rate of the
        // spills is smaller than 2 secs.
        Double_t dtime = tsec + (tnanosec/1e9);
        // Look up the actual time from this, previous and next spill
        Double_t time_this =  timeStampD;
        Double_t time_prev = 0.0; 
        Double_t time_next = 0.0;
        if (index<nentries-1) {
            fChain->GetEntry(index+1);
            time_next = timeStampD;
        }    
        if (index>0) {
            fChain->GetEntry(index-1);
            time_prev = timeStampD;
        }    
        Double_t Deltat_this = fabs(time_this - dtime);
        Double_t Deltat_prev = fabs(time_prev - dtime);
        Double_t Deltat_next = fabs(time_next - dtime);
        MSG("BeamSummary",Msg::kDebug) << "Time difference i = " << Deltat_this<< endl;
        MSG("BeamSummary",Msg::kDebug) << "Time difference i-1 = " << Deltat_prev<< endl;
        MSG("BeamSummary",Msg::kDebug) << "Time difference i+1 = " << Deltat_next<< endl;
        if (Deltat_this>Deltat_prev) index--;
        else if (Deltat_this>Deltat_next) index++;
        fChain->GetEntry(index);
    }
    lastindex=index;
    return found;
}

---