Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Re-add get1DVarHist back into samplePDFFDBase #172

Merged
merged 3 commits into from
Oct 22, 2024
Merged

Re-add get1DVarHist back into samplePDFFDBase #172

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
653304e
Add get1DVarHist into core code for ability to plot variables defined…
dbarrow257 Oct 17, 2024
13e1e36
Merge remote-tracking branch 'origin/develop' into feature_FDVariable…
dbarrow257 Oct 17, 2024
ef6d901
Merge branch 'develop' into feature_FDVariablePlotting
dbarrow257 Oct 21, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions samplePDF/FDMCStruct.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,8 @@ struct fdmc_base {
int nupdg;
int nupdgUnosc;
bool signal; // true if signue
int nEvents; // how many MC events are there
int nEvents; // how many MC events are there
double ChannelIndex;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is ChannelIndex?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In FD code, there is the concept of oscillation channels. Each oscillation channel has an associated FDMCStruct object. These channels' indexing depends entirely on the ordering defined in the sample YAML.

If you want to loop over each oscillation channel and make an event rate breakdown by oscChan, you need some way to select events in oscillationChannel j (j=6, for example). To do this from ReturnKinematicParameter, we need to store the index of the oscillation channel somewhere in it's memory, thus ChannelIndex.

You can get name of the j'th oscillationChannel index from FDMCStruct.flavourName

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In DUNE code for example:

const double* samplePDFDUNEAtm::GetPointerToKinematicParameter(KinematicTypes KinPar, int iSample, int iEvent) {
  double* KinematicValue;


  switch (KinPar) {
  case kTrueNeutrinoEnergy:
    KinematicValue = &(dunemcSamples[iSample].rw_etru[iEvent]);
    break;
  case kRecoNeutrinoEnergy:
    KinematicValue = &(dunemcSamples[iSample].rw_erec[iEvent]);
    break;
  case kTrueCosZ:
    KinematicValue = &(dunemcSamples[iSample].rw_truecz[iEvent]);
    break;
  case kRecoCosZ:
    KinematicValue = &(dunemcSamples[iSample].rw_theta[iEvent]);
    break;
  case kOscChannel:
    KinematicValue = &(MCSamples[iSample].ChannelIndex);
    break;
  case kMode:
    KinematicValue = MCSamples[iSample].mode[iEvent];
    break;
  default:
    MACH3LOG_ERROR("Unknown KinPar: {}",KinPar);
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  
  return KinematicValue;
}

std::string flavourName;

int **Target; // target the interaction was on
Expand Down Expand Up @@ -52,7 +53,7 @@ struct fdmc_base {
double *rw_upper_xbinedge; // upper to check if Eb has moved the erec bin
double *rw_upper_upper_xbinedge; // upper to check if Eb has moved the erec bin

int **mode;
double **mode;

const _float_ **osc_w_pointer;
double *xsec_w;
Expand Down
63 changes: 61 additions & 2 deletions samplePDF/samplePDFFDBase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1537,7 +1537,7 @@ void samplePDFFDBase::InitialiseSingleFDMCObject(int iSample, int nEvents_) {
fdobj->rw_lower_lower_xbinedge = new double [fdobj->nEvents];
fdobj->rw_upper_xbinedge = new double [fdobj->nEvents];
fdobj->rw_upper_upper_xbinedge = new double [fdobj->nEvents];
fdobj->mode = new int*[fdobj->nEvents];
fdobj->mode = new double*[fdobj->nEvents];
fdobj->nxsec_norm_pointers = new int[fdobj->nEvents];
fdobj->xsec_norm_pointers = new const double**[fdobj->nEvents];
fdobj->xsec_norms_bins = new std::list< int >[fdobj->nEvents];
Expand All @@ -1554,7 +1554,7 @@ void samplePDFFDBase::InitialiseSingleFDMCObject(int iSample, int nEvents_) {
for(int iEvent = 0 ;iEvent < fdobj->nEvents ; ++iEvent){
fdobj->rw_etru[iEvent] = &fdobj->Unity;
//fdobj->rw_truecz[iEvent] = &fdobj->Unity;
fdobj->mode[iEvent] = &fdobj->Unity_Int;
fdobj->mode[iEvent] = &fdobj->Unity;
fdobj->Target[iEvent] = 0;
fdobj->NomXBin[iEvent] = -1;
fdobj->NomYBin[iEvent] = -1;
Expand Down Expand Up @@ -1604,3 +1604,62 @@ void samplePDFFDBase::InitialiseSplineObject() {

fillSplineBins();
}

TH1* samplePDFFDBase::get1DVarHist(std::string ProjectionVar_Str, std::vector< std::vector<double> > SelectionVec, int WeightStyle, TAxis* Axis) {
//DB Grab the associated enum with the argument string
int ProjectionVar_Int = ReturnKinematicParameterFromString(ProjectionVar_Str);

//DB Need to overwrite the Selection member variable so that IsEventSelected function operates correctly.
// Consequently, store the selection cuts already saved in the sample, overwrite the Selection variable, then reset
std::vector< std::vector<double> > tmp_Selection = Selection;

std::vector< std::vector<double> > SelectionVecToApply;

//DB Add all the predefined selections to the selection vector which will be applied
for (int iSelec=0;iSelec<Selection.size();iSelec++) {
SelectionVecToApply.emplace_back(Selection[iSelec]);
}

//DB Add all requested cuts from the arguement to the selection vector which will be applied
for (int iSelec=0;iSelec<SelectionVec.size();iSelec++) {
SelectionVecToApply.emplace_back(SelectionVec[iSelec]);
}

//DB Check the formatting of all requested cuts, should be [cutPar,lBound,uBound]
for (int iSelec=0;iSelec<SelectionVecToApply.size();iSelec++) {
if (SelectionVecToApply[iSelec].size()!=3) {
MACH3LOG_ERROR("Selection Vector[{}] is not formed correctly. Expect size == 3, given: {}",iSelec,SelectionVecToApply[iSelec].size());
throw MaCh3Exception(__FILE__, __LINE__);
}
}

//DB Set the member variable to be the cuts to apply
Selection = SelectionVecToApply;

//DB Define the histogram which will be returned
TH1D* _h1DVar;
if (Axis) {
_h1DVar = new TH1D("","",Axis->GetNbins(),Axis->GetXbins()->GetArray());
} else {
std::vector<double> xBinEdges = ReturnKinematicParameterBinning(ProjectionVar_Str);
double xbin_edges[xBinEdges.size()];
for (unsigned int i=0;i<xBinEdges.size();i++) {xbin_edges[i] = xBinEdges[i];}
_h1DVar = new TH1D("", "", xBinEdges.size()-1, xbin_edges);
}

//DB Loop over all events
for (int iSample=0;iSample<getNMCSamples();iSample++) {
for (int iEvent=0;iEvent<getNEventsInSample(iSample);iEvent++) {
if (IsEventSelected(iSample,iEvent)) {
double Weight = GetEventWeight(iSample,iEvent);
double Var = ReturnKinematicParameter(ProjectionVar_Int,iSample,iEvent);
_h1DVar->Fill(Var,Weight);
}
}
}

//DB Reset the saved selection
Selection = tmp_Selection;

return _h1DVar;
}
28 changes: 24 additions & 4 deletions samplePDF/samplePDFFDBase.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,30 @@ class samplePDFFDBase : public samplePDFBase
virtual void setupSplines(fdmc_base *FDObj, const char *SplineFileName, int nutype, int signal){};
void ReadSampleConfig();

int getNMCSamples() {return MCSamples.size();}

int getNEventsInSample(int iSample) {
if (iSample < 0 || iSample > getNMCSamples()) {
MACH3LOG_ERROR("Invalid Sample Requested: {}",iSample);
throw MaCh3Exception(__FILE__ , __LINE__);
}
return MCSamples[iSample].nEvents;
}

std::string getFlavourName(int iSample) {
if (iSample < 0 || iSample > getNMCSamples()) {
MACH3LOG_ERROR("Invalid Sample Requested: {}",iSample);
throw MaCh3Exception(__FILE__ , __LINE__);
}
return MCSamples[iSample].flavourName;
}

TH1* get1DVarHist(std::string ProjectionVar, std::vector< std::vector<double> > SelectionVec = std::vector< std::vector<double> >(), int WeightStyle=0, TAxis* Axis=nullptr);

//ETA - new function to generically convert a string from xsec cov to a kinematic type
virtual inline int ReturnKinematicParameterFromString(std::string KinematicStr) = 0;
virtual inline std::string ReturnStringFromKinematicParameter(int KinematicVariable) = 0;

protected:
/// @brief DB Function to determine which weights apply to which types of samples pure virtual!!
virtual void SetupWeightPointers() = 0;
Expand Down Expand Up @@ -128,10 +152,6 @@ class samplePDFFDBase : public samplePDFBase
virtual const double* GetPointerToKinematicParameter(std::string KinematicParamter, int iSample, int iEvent) = 0;
virtual const double* GetPointerToKinematicParameter(double KinematicVariable, int iSample, int iEvent) = 0;

//ETA - new function to generically convert a string from xsec cov to a kinematic type
virtual inline int ReturnKinematicParameterFromString(std::string KinematicStr) = 0;
virtual inline std::string ReturnStringFromKinematicParameter(int KinematicVariable) = 0;

// Function to setup Functional and shift parameters. This isn't idea but
// do this in your experiment specific code for now as we don't have a
// generic treatment for func and shift pars yet
Expand Down