MVariableOperation.cc

Go to the documentation of this file.

#include "MVariableOperation.hh"
#include "QBaseType.hh"
#include "QEvent.hh"
#include <algorithm> 
#include "QEventList.hh"
#include "QRawEvent.hh"
   
REGISTER_MODULE(MVariableOperation)
 
using namespace Diana;
using std::cout;
using std::endl;
using std::vector;
 
void MVariableOperation::Init(QEvent& ev)
{
    //Getting Arguments from cfg
    fInputVar= GetString("Variable","",true);
    fOutputVar = GetString("OutputVar","",true);
    fFunction = GetString("Function","",true);
        
        fParamNames=FindVars(fFunction);
        
        fOtherVars = FindQuantities(fFunction);
    
    fChannels = GetVectorInt("Channels",vector<int>(),true);
    std::sort(fChannels.begin(),fChannels.end());
    channelparams.resize(fChannels.size());
    for(size_t i =0; i<fChannels.size();i++){
        channelparams[i].channel_number=fChannels[i];
    }
    cout<<endl;
 
    StringParams = GetString("Parameters","",true);
    FindParameters(StringParams);
    cout<<"InputVar\t"<<fInputVar.GetStringLabel()<<endl;
    cout<<"OutputVar\t"<<fOutputVar<<endl;
    cout<<"Function\t"<<fFunction<<endl;
    cout<<"Parameter Names: ";
    for(size_t k=0;k<fParamNames.size();k++){
        cout<<"["<<fParamNames[k]<<"] ";
    }
    cout<<endl;
    cout<<"Parameters:"<<endl;
    for(size_t k=0;k<channelparams.size(); k++){
        if(channelparams[k].fParameters.size()!=fParamNames.size()){
                Panic("Number of parameters in formula don't match parameter values provided for channel %d",channelparams[k].channel_number);
        }
        cout<<"\tChannel "<< channelparams[k].channel_number<<" Parameters: ";
        for(size_t i =0;i<channelparams[k].fParameters.size();i++){
            cout<<channelparams[k].fParameters[i]<<" ";
        }
        cout<<endl;
    }
    
        fFunc = TF1("EvalFunc",fFunction,-1e7,1e7);
        cout<<"Found Quantities"<<endl;
        for(size_t i =0;i<fOtherVars.size();i++){
                cout<<"\t["<<fParamNames.size()+i<<"] "<<fOtherVars[i].GetStringLabel()<<endl;
                SetReq(fOtherVars[i],ev);
        }
        SetReq(fInputVar,ev);
    ev.Require<QPulseInfo>("DAQ","PulseInfo");
    ev.Add<QDouble>(fOutputVar.Data()); 
}
 
void MVariableOperation::SetReq(QEventLabel EvLab,QEvent& ev){
        TString lab = TString(EvLab.GetStringLabel());
        if(lab.Contains("BaselineData")){
                ev.Require<QBaselineData>(EvLab.GetOwner().c_str(),"BaselineData"); 
        }else{
                ev.RequireByLabel<QDouble>(EvLab);
        }
}               
        
double MVariableOperation::GetValue(QEventLabel EvLab,QEvent& ev){
        TString lab = TString(EvLab.GetStringLabel());
        double variable=0;
        if(lab.Contains("BaselineData")){
                const QBaselineData& bData=ev.Get<QBaselineData>(EvLab.GetOwner().c_str(),"BaselineData"); 
                TString varlab = lab(lab.First(".")+1,lab.Length());
                if(varlab.CompareTo("fBaseline")==0){
                variable = bData.GetBaseline();
  
                }else if(varlab.CompareTo("fBaselineFlatRMS")==0){
                variable = bData.GetBaselineFlatRMS();
  
                }else if(varlab.CompareTo("fBaselineIntercept")==0){
                variable = bData.GetBaselineIntercept();
  
                }else if(varlab.CompareTo("fBaselineSlope")==0){
                variable = bData.GetBaselineSlope();
  
                }else if(varlab.CompareTo("fBaselineSlopeRMSWindow")==0){
                variable = bData.GetBaselineSlopeRMSWindow();
  
                }else if(varlab.CompareTo("fBaselineRMS")==0){
                variable = bData.GetBaselineRMS();
  
                }else if(varlab.CompareTo("fRightBaseline")==0){
                variable = bData.GetRightBaseline();
  
                }else if(varlab.CompareTo("fRightBaselineRMS")==0){
                variable = bData.GetRightBaselineInRMS();
                }else if(varlab.CompareTo("fRightLeftBaseline")==0){
                variable = bData.GetRightLeftBaseline();
                }else if(varlab.CompareTo("fRightLeftBaselineInRMS")==0){
                variable = bData.GetRightLeftBaselineInRMS();
                }else{
                Panic(TString("'")+varlab+TString("' if not a valid choice for the 'Variable' field."));
                }
        }else{
                variable = ev.GetByLabel<QDouble>(EvLab);
        }
        return variable;
 
}
void MVariableOperation::Do(QEvent& ev)
{
    const QPulseInfo& pulseInfo = ev.Get<QPulseInfo>("DAQ","PulseInfo");
    const int chan = pulseInfo.GetChannelId();
    size_t chan_counter=channelparams.size();
    for(size_t i =0;i<channelparams.size();i++){
        if(channelparams[i].channel_number==chan){
            chan_counter=i;
            break;
        }
    }
    
    if(chan_counter==channelparams.size()){
        Panic("Channel %d not found in provided channels from cfg!",chan);
    }
    for(size_t i =0;i<fParamNames.size();i++){
                fFunc.SetParameter(fParamNames[i],channelparams[chan_counter].fParameters[i]);
        }
    
        const double variable=GetValue(fInputVar,ev);
        for(size_t i =0;i<fOtherVars.size();i++){
                const double tmp = GetValue(fOtherVars[i],ev);
                fFunc.SetParameter(fParamNames.size()+i,tmp);
        }
        ev.Get<QDouble>(fOutputVar.Data()) = fFunc.Eval(variable);
}  
 
vector<Int_t> MVariableOperation::FindVars(TString instring){
        vector<TString> out;
        size_t idxopen=0;
        size_t idxclose=0;
        bool condition = true;
        while(condition){
                idxclose=TString(instring(int(idxopen),instring.Length()-int(idxopen))).First(']') + idxopen;
                idxopen=TString(instring(int(idxopen),instring.Length()-int(idxopen))).First('[') + idxopen;
                if(idxclose<=idxopen){
                        condition=false;
                        continue;
                }
 
        TString substr = TString(instring(int(idxopen),int(idxclose+1-idxopen)));
                if(!substr.BeginsWith("[")){
                        condition=false;
                        continue;
                }
                
                bool found=false;
                for(size_t j=0;j<out.size();j++){
                        if(out[j].CompareTo(TString(substr(1,substr.Length()-2)))==0){
                                found=true;
                        }
                }
                if(!found){
                        out.push_back(TString(substr(1,substr.Length()-2)));
                }
                idxopen=idxclose+1;
        }
        vector<Int_t> out2;
    for(size_t j=0;j<out.size();j++){
                out2.push_back(out[j].Atoi());
        }
        std::sort (out2.begin(), out2.end()); 
        return out2;
 
        
}
 
vector<Diana::QEventLabel> MVariableOperation::FindQuantities(TString instring){
        vector<TString> out;
        size_t idxopen=0;
        size_t idxclose=0;
        bool condition=true;
        while(condition){
                idxclose=TString(instring(int(idxopen),instring.Length()-int(idxopen))).First('&') + idxopen;
                idxopen=TString(instring(int(idxopen),instring.Length()-int(idxopen))).First('^') + idxopen;
                if(idxclose<=idxopen){
                        condition=false;
                        continue;
                }
 
        TString substr = TString(instring(int(idxopen),int(idxclose+1-idxopen)));
                if(!substr.BeginsWith("^")){
                        condition=false;
                        continue;
                }
                
                bool found=false;
                for(size_t j=0;j<out.size();j++){
                        if(out[j].CompareTo(TString(substr(1,substr.Length()-2)))==0){
                                found=true;
                        }
                }
                if(!found){
                        out.push_back(TString(substr(1,substr.Length()-2)));
                }
                idxopen=idxclose+1;
        }
        TString modin = instring;
        bool stop=false;
        while(!stop){
                idxopen = modin.First("^");
                idxclose = modin.First("&");
                TString v = TString(modin(idxopen+1,idxclose-idxopen-1));
                if(idxclose<=idxopen){
                        stop=true;
                        continue;
                }
                size_t idx=0;
                condition=true;
                while(condition){
                        if(idx>=out.size()){
                                condition=false;
                                continue;
                        }
                        if (out[idx].CompareTo(v)!=0){
                                idx++;
                        }else{
                                condition=false;
                                continue;
                        }
                                
                }
                TString idxsub;
                idxsub.Form("%lu",fParamNames.size()+idx);
                modin = modin(0,int(idxopen))+"["+idxsub+"]"+modin(int(idxclose)+1,int(modin.Length()));
        }
        fFunction =modin;
        vector<Diana::QEventLabel> out2;
        for(size_t j=0;j<out.size();j++){
                out2.push_back(Diana::QEventLabel(out[j]));
        }
        return out2;
 
        
}
 
void MVariableOperation::FindParameters(TString instring){
    //Remove Spaces
    TString manip="";
    for(int i =0; i< instring.Length();i++){
        if(TString(instring(i,1)).CompareTo(" ")!=0){
            manip+=TString(instring(i,1));
        }
    }
 
    //Find all "}," indicating channel separation
    vector<int> separations={-1};
    bool found=true;
    while(found){
        TString tmpstr=TString(manip(separations[separations.size()-1]+1,manip.Length()));
        int tmp = -1;
        for(int kk=0;kk<tmpstr.Length()-1;kk++){
            if(TString(tmpstr(kk,2)).CompareTo("},")==0){
                tmp=kk;
                break;
            }
        }
        if(tmp==-1 or TString(manip(separations[separations.size()-1]+1,manip.Length())).Length()/2+TString(manip(separations[separations.size()-1]+1,manip.Length())).Length()%2==tmp){
            found=false;
        }else{
            separations.push_back(tmp+separations[separations.size()-1]+2);
        }
    }
    //separating channels
    size_t found_channels=0;
    for(size_t i =0;i<separations.size();i++){
        TString params;
        if(i < separations.size()-1){
            params = TString(manip(separations[i]+1,separations[i+1]-separations[i]-1));
        }else{
            params = TString(manip(separations[i]+1,manip.Length()));
        }
 
        //Channel identification
        int tmp = params.First(":");
        if(tmp==-1 or tmp ==params.Length()){
            Panic("Please provide channel ! String %s not valid! Format channel:{param1,param2,..}",params.Data());
        }
        tmp++;
        int ch = TString(params(0,tmp)).Atoi();
        params = TString(params(tmp,params.Length()));
        
        //Parameter identification
        int start=params.First("{");
        if(start==-1 or start ==params.Length()){
            Panic("Wrong format ! String %s not valid! Format channel:{param1,param2,..}",params.Data());
        }
        int stop=params.First("}");
        if(stop==-1 or stop ==params.Length()){
            Panic("Wrong format ! String %s not valid! Format channel:{param1,param2,..}",params.Data());
        }
        start++;
        params = TString(params(start,stop-start));
        found = true;
        vector<int> paramssep={-1};
        while(found){
            TString tmpstr2=TString(params(paramssep[paramssep.size()-1]+1,params.Length()));
            tmp = -1;
            for(int kk=0;kk<tmpstr2.Length();kk++){
                if(TString(tmpstr2(kk,1)).CompareTo(",")==0){
                    tmp=kk;
                    break;
                }
            }
            if(tmp==-1 or TString(params(paramssep[paramssep.size()-1],params.Length())).Length()==tmp){
                found=false;
            }else{
                paramssep.push_back(tmp+paramssep[paramssep.size()-1]+1);
            }
        }
        
        size_t chcounter=channelparams.size();
        
        for(size_t k=0;k<channelparams.size(); k++){
            if(channelparams[k].channel_number == ch){
                chcounter=k;
                
            }
        }
        if(chcounter == channelparams.size()){
            Panic("Selected Channel not found! Please include it in the cfg parameter!");
        }else{
            found_channels++;
        }
        
        for(size_t i =0;i<paramssep.size();i++){
            TString params2;
            if(i == 0){
                params2 = TString(params(0,paramssep[1]));
            }else if(i < paramssep.size()-1){
                params2 = TString(params(paramssep[i]+1,paramssep[i+1]-paramssep[i]));
            
            }else{
                params2 = TString(params(paramssep[i]+1,params.Length()));
            }
            channelparams[chcounter].fParameters.push_back(params2.Atof());
            
        }
    }
    if(found_channels <channelparams.size()){
        Panic("Please indicate parameters for all channels!");
    }else if(found_channels >channelparams.size()){
        Panic("Something weird happening... more channels found than present");
    }
}
 
void MVariableOperation::Done()
{
    /* This method is called at the end of the event loop. 
     * Here you can:
     *
     * 1) Operate on the sequence execution parameters(see QBaseModule.hh).
     *
     * 2) Read/Write global data.
     */
}