BMSpillAna.cxx

Go to the documentation of this file.

#include "BeamDataUtil/BMSpillAna.h"
#include "Validity/VldTimeStamp.h"

#include <MessageService/Msg.h>
#include <MessageService/MsgService.h>
CVSID("$Id: BMSpillAna.cxx,v 1.18 2008/08/15 18:28:17 loiacono Exp $");

#include "Conventions/Munits.h"

#include "TMath.h"
#include <cmath>

#include <cstdio>


BMSpillAna::BMSpillAna()
    :fSpill(0),fUserCuts(),fUseDBCuts(true),fCutsSet(0),fResPtr(),
     fResID(-1),fTimeDiff(-99999)
{
    this->ChangeCutValues(this->DefaultConfig());
}


BMSpillAna::~BMSpillAna()
{}


const Registry& BMSpillAna::DefaultConfig() const
{
    static Registry r;
    if (r.Size() == 0) {
        r.UnLockValues();
        r.Set("TimeDiffMax", 1.0*Munits::s);

        r.Set("PosTgtXMin",-2.0*Munits::mm);
        r.Set("PosTgtXMax",-0.01*Munits::mm);

        r.Set("PosTgtYMin",0.01*Munits::mm);
        r.Set("PosTgtYMax",2.0*Munits::mm);

        r.Set("WidXMin",0.1*Munits::mm);
        r.Set("WidXMax",1.5*Munits::mm);

        r.Set("WidYMin",0.1*Munits::mm);
        r.Set("WidYMax",2.0*Munits::mm);

        r.Set("UseSpotSizeCut",0);
        
        r.Set("UseProfMonOut",1);

        r.Set("TorIntMin",0.50);
        r.Set("TorIntMax",50.0);

        r.Set("HornCurMin",-2.0e5*Munits::ampere);
        r.Set("HornCurMax",-1.55e5*Munits::ampere);

        r.Set("TargetIn",1);
        r.Set("BeamType",-1);

        r.Set("FracOnTargetMin",0.0);
        r.Set("FracOnTargetMax",1.0);

        r.LockValues();
    }
    return r;
}

void BMSpillAna::Config(const Registry& r)
{
    fUserCuts=r;
    this->ChangeCutValues(r);
}

void BMSpillAna::ChangeCutValues(const Registry& r)
{
    if (r.Size()>0){
        r.Get("TimeDiffMax",fTimeDiffMax);
        
        r.Get("PosTgtXMin",fPosTgtXMin);
        r.Get("PosTgtXMax",fPosTgtXMax);
        
        r.Get("PosTgtYMin",fPosTgtYMin);
        r.Get("PosTgtYMax",fPosTgtYMax);
        
        r.Get("WidXMin",fWidXMin);
        r.Get("WidXMax",fWidXMax);
        
        r.Get("WidYMin",fWidYMin);
        r.Get("WidYMax",fWidYMax);

        r.Get("UseSpotSizeCut",fUseSpotSizeCut);
        
        r.Get("UseProfMonOut",fUseProfMonOut);
        
        r.Get("TorIntMin",fTorIntMin);
        r.Get("TorIntMax",fTorIntMax);
        
        r.Get("HornCurMin",fHornCurMin);
        r.Get("HornCurMax",fHornCurMax);
        
        r.Get("TargetIn",fTargetIn);
        r.Get("BeamType",fBeamType);
        
        r.Get("FracOnTargetMin",fFracOnTargetMin);
        r.Get("FracOnTargetMax",fFracOnTargetMax);
    }        
}
    
void BMSpillAna::ApplyUserCuts()
{
    this->ChangeCutValues(fUserCuts);
}


void BMSpillAna::SetSpill(const BeamMonSpill& spill)
{
    fSpill = &spill;
    // Reset the time difference every time to make sure that the user
    // updates it.
    this->SetTimeDiff(-99999);

}

void BMSpillAna::SetSpill(const NtpBDLiteRecord& ntpbdr, BeamMonSpill& spill)
{
    fSpill = &spill;

    this->SetTimeDiff(ntpbdr.GetHeader().GetTimeDiffStreamSpill());
    
    spill.fDaeTime = ntpbdr.GetHeader().GetEarliestTimeStamp();
    spill.fVmeTime = ntpbdr.GetHeader().GetEarliestTimeStamp();
    spill.fTor101  = ntpbdr.tor101;
    spill.fTr101d  = ntpbdr.tr101d;
    spill.fTortgt  = ntpbdr.tortgt;
    spill.fTrtgtd  = ntpbdr.trtgtd;
    spill.fHornCur = ntpbdr.horncur;
    for (int i=0;i<6;++i){
        spill.fTargBpmX[i] = ntpbdr.bposx[i];
        spill.fTargBpmY[i] = ntpbdr.bposy[i];
        spill.fBpmInt[i]   = ntpbdr.bpmint[i];
    }
    spill.fProfWidX = ntpbdr.bwidx;
    spill.fProfWidY = ntpbdr.bwidy;
    spill.fHadInt = ntpbdr.hadint;
    spill.fMuInt1 = ntpbdr.muint1;
    spill.fMuInt2 = ntpbdr.muint2;
    spill.fMuInt3 = ntpbdr.muint3;
    union {int integer; BeamMonSpill::StatusBits bits; } status;
    status.integer = ntpbdr.GetHeader().GetStatus();
    spill.SetStatusBits(status.bits);

}

void BMSpillAna::SetSnarlTime(const VldTimeStamp& vs_snarl)
{
    if (fSpill==0){
        MSG("BMSpillAna", Msg::kError)<< "Set the spill object before setting this time" << endl;
        return;
    }
    this->SetTimeDiff(vs_snarl-fSpill->SpillTime());
}

Bool_t BMSpillAna::SelectSpill()
{
    if (fTimeDiff==-99999)
        MSG("BMSpillAna", Msg::kWarning) <<
            "Time difference seems not to be set correctly" << endl;
    
    
    // If using the database, check whether the cuts need to be
    // updated
        
    if (fUseDBCuts){
        MAXMSG("BMSpillAna", Msg::kDebug,5) <<
            "Using database cuts" << endl;


        VldContext vc(Detector::kNear,SimFlag::kData,fSpill->SpillTime());
        Int_t nrows = fResPtr.NewQuery(vc,fCutsSet);
        if (nrows==0){
            MAXMSG("BMSpillAna",Msg::kWarning,20)
                << "No cuts found in database. This should not happen!"
                << endl;
        }
        else {
            if (nrows>1) {
                MAXMSG("BMSpillAna",Msg::kWarning,20)
                    << "More than one row found for VldContext " 
                    << vc << endl;
                MAXMSG("BMSpillAna",Msg::kWarning,20)
                    << "   --> Will only use first row! " << endl;
            }
            
            MAXMSG("BMSpillAna", Msg::kDebug,5) <<
                "Rows found in BEAMMONCUTS table" << endl;
            // If the cuts are still the same (i.e. the data in the
            // DbiResultPointer is the same, do nothing, just keep the
            // current cut values, else change the cut values and
            // apply user cuts

            Int_t newid = fResPtr.GetResultID();
            if (newid != fResID){
                MSG("BMSpillAna", Msg::kDebug) <<
                    "Cuts need to be updated for spill at " << fSpill->SpillTime()  << endl;
                fResID = newid;
                const BeamMonCuts* bmc = fResPtr.GetRow(0);
                Registry newreg;
                bmc->FillRegistry(&newreg);
                this->ChangeCutValues(newreg);

                Int_t loglevel = MsgService::Instance()->GetStream("BMSpillAna")->GetLogLevel();
                if (loglevel == Msg::kDebug){
                    cout << "Database cuts" << endl;
                    this->Print();
                }

                this->ApplyUserCuts();
                if (loglevel == Msg::kDebug){
                    cout << "After user cuts" << endl;
                    this->Print();
                }
            }
        }        
    }

    if (fabs(fTimeDiff)>fTimeDiffMax) return false;
    
    Double_t xmean=0;
    Double_t ymean=0;
    Double_t xrms=0;
    Double_t yrms=0;
    fSpill->BpmAtTarget(xmean,ymean,xrms,yrms);
    if (xmean < fPosTgtXMin || xmean > fPosTgtXMax) return false;
    if (ymean < fPosTgtYMin || ymean > fPosTgtYMax) return false;


    // Never select fit failures
    if (fSpill->fProfWidX < -0.1*Munits::mm
        || fSpill->fProfWidY < -0.1*Munits::mm) return false;
    //
    if (fUseSpotSizeCut){
        // The widths should never be negative, unless it's due to a
        // fit failure, but these are already excluded above. The
        // values will be zero if the profile monitor is out.
        Double_t spot_size=0;
        if (fSpill->fProfWidX>0 && fSpill->fProfWidY>0)
            spot_size = fSpill->fProfWidX*fSpill->fProfWidY;

        // put this defualt to a small negative value, so that the
        // pribility exists to select the spill if the profile monitor
        // is out, i.e. widths are zero.
        Double_t spot_size_min = -0.01*Munits::mm*Munits::mm;
        if (fWidXMin>0 && fWidYMin>0)
            spot_size_min = fWidXMin*fWidYMin;

        // Assume people are smart enough not to put negative values
        // for the upper limit on the beam widths. If they do, they
        // will loose all spills...
        Double_t spot_size_max = fWidXMax*fWidYMax;

        if (!(fUseProfMonOut) && spot_size < spot_size_min ) return false;
        else if (spot_size > spot_size_max) return false;

    }
    else {       
        if (!(fUseProfMonOut) && fSpill->fProfWidX < fWidXMin) return false;
        else if (fSpill->fProfWidX > fWidXMax) return false;
        
        if (!(fUseProfMonOut) && fSpill->fProfWidY < fWidYMin) return false;
        else if (fSpill->fProfWidY > fWidYMax) return false;
    }


    //changed to Trtgtd from Tortgt on 20080814
    if (fSpill->fTrtgtd < fTorIntMin || fSpill->fTrtgtd > fTorIntMax){
        return false;
    }
    
    // FIXME: values in the database are in kAmps, not in
    // Munits. This will get fixed eventually, but for now, there
    // need to be an explicit factor of 1e3.
    if (fSpill->fHornCur*1e3 < fHornCurMin || fSpill->fHornCur*1e3 > fHornCurMax){
        return false;
    }

    //added  fTargetIn < 0 || on 20080814
    if ( fTargetIn < 0 || (fTargetIn >= 0 && (Bool_t)fSpill->GetStatusBits().target_in != (Bool_t)fTargetIn) ){ 
        return false;
    }

    if (fBeamType >=0 && (fSpill->BeamType() != fBeamType)){
        return false;
    }

    if (fFracOnTargetMin > 0 &&  this->FractionOnTarget() < fFracOnTargetMin) {
        return false;
    }
    if (fFracOnTargetMax < 1 &&  this->FractionOnTarget() > fFracOnTargetMax) {
        return false;
    }
    
    return true;
}


Double_t BMSpillAna::FractionOnTarget()
{
    Double_t frac = -1;
//
    Double_t bposx = 0;
    Double_t bwidx = 0;
//
    Double_t bposy = 0;
    Double_t bwidy = 0;
    
    // Get the position at target 
    fSpill->BpmAtTarget(bposx,bposy,bwidx,bwidy);
    
    // Set the values of the widths to prof mon fits
    bwidx = fSpill->fProfWidX; 
    bwidy = fSpill->fProfWidY; 

    // information of profile monitors not available or fits failed
    if (bwidx<=0 || bwidy<=0) return frac;

    // large values also indicate fit failure
    if (bwidx>5*Munits::mm || bwidy>05*Munits::mm) return frac;

// Set the target center and edges in BPM coordinates in x and y
    Double_t tedgx = 3.2*Munits::mm;
    Double_t toffx = -1.2*Munits::mm;
    Double_t tedgy = 5.5*Munits::mm;
    Double_t toffy = 0.9*Munits::mm;


    frac = CalcFracOnTarget(bposx,bwidx,tedgx,toffx);
    frac *= CalcFracOnTarget(bposy,bwidy,tedgy,toffy);
    
    return frac;    
}

Double_t BMSpillAna::CalcFracOnTarget(Double_t bpos, Double_t bwid, Double_t tedg, Double_t toff)
{
    Double_t powl = tedg+bpos-toff;
    powl /= bwid;
    powl /= TMath::Sqrt(2.0);
    Double_t powu = tedg-bpos+toff;
    powu /= bwid;
    powu /= TMath::Sqrt(2.0);
    Double_t frac = TMath::Erf(powl)/2 +  TMath::Erf(powu)/2; 
    return frac;
}

void BMSpillAna::Print()
{
    cout << endl << "Beam Monitoring Cut Values:" << endl;
    cout         << "===========================" << endl;

    printf(" > Maximum time diffrence (s):         %5.3f\n",fTimeDiffMax);
    printf(" > Spill intensity (1e12 pot):     [%5.2f,%5.2f]\n",
           fTorIntMin, fTorIntMax);
    printf(" > Horn Current (kA):             [%+5.1f,%+5.1f]\n",
           fHornCurMin/Munits::ampere/1e3, fHornCurMax/Munits::ampere/1e3);
    printf(" > Target in/out:                      %3d\n",fTargetIn);
    printf(" > fBeamType:                          %3d\n",fBeamType);    
    printf(" > Horizontal beam position (mm): [%+5.3f,%+5.3f]\n",
           fPosTgtXMin/Munits::mm, fPosTgtXMax/Munits::mm);
    printf(" > Vertical beam position (mm):   [%+5.3f,%+5.3f]\n",
           fPosTgtYMin/Munits::mm, fPosTgtYMax/Munits::mm);
    printf(" > Horizontal beam width a (mm):  [%+5.3f,%+5.3f]\n",
           fWidXMin/Munits::mm, fWidXMax/Munits::mm);
    printf(" > Vertical beam width a (mm):    [%+5.3f,%+5.3f]\n",
           fWidYMin/Munits::mm, fWidYMax/Munits::mm);
    printf(" > Use spill when prof mon out:        %3d\n",fUseProfMonOut);
    printf(" > Beam fraction on target:       [%5.3f,%5.3f]\n",
           fFracOnTargetMin, fFracOnTargetMax);

}

---