#include "GENIE/Framework/Numerical/RandomGen.h"
#include "GENIE/Tools/Flux/GSimpleNtpFlux.h"
#include "GENIE/Framework/Utils/UnitUtils.h"
#include "TSystem.h"
#include "TStopwatch.h"
#include "TLorentzVector.h"
#include "TNtuple.h"
#include "TFile.h"
#include "TChain.h"
#include "TF1.h"
#include "TH1F.h"
#include "TRandom.h"
#include <iostream>
#include <iomanip>
#include <string>
#include <sstream>
#include <set>
#include <algorithm>

Functions
double	pick_energy (double, double)

void	warp_entry (GSimpleNtpEntry entry, GSimpleNtpAux aux, GSimpleNtpNuMI *numi)

void	warp_meta (GSimpleNtpMeta *meta, string fnamein)

double	fit_fcn (double x, double par)

void	update_flavors (GSimpleNtpMeta *meta, std::set< int > &gFlavors)

void	warp_gsimple (string fnameout="gsimple_output.root", string fnamein="gsimple_input.root", bool dowarp=true)

Variables
TF1 *	fcn

bool	gDoWarp = true

Function Documentation

double fit_fcn	(	double *	x,
		double *	par
	)

Definition at line 109 of file warp_gsimple.C.

References energy, update_flavors(), and submit_syst::y.

Referenced by warp_gsimple().

                                       {
   double energy = x[0];
   auto profile = [](double y, double* p){
     return p[0]+p[1]*TMath::Power(y, -p[2])+p[3]*TMath::Power(y, -p[4])+p[5]*y;
   };
   double thresholdE = 0.1;
   if(energy > thresholdE)
     return profile(energy, par);
   else
     return profile(100., par);
 }

double pick_energy	(	double	emin,
		double	emax
	)

Definition at line 93 of file warp_gsimple.C.

References fcn.

Referenced by warp_entry().

 {
   // pick a distribution of shape 1/E
   // between the values [emin, emax]
   // must have low energy cut to avoid infrared divergence
   // double r = gRandom->Rndm(); // flat variate between (0,1]
   // // double c = log(emin);
   // // double A = 1. / (log(emax) - c);
   // // return exp(r/A + c);
   // double c = emin;
   // double A = 1./(emax - c);
   // return r/A + c;
   return fcn->GetRandom(emin, emax);
 
 }

void update_flavors	(	GSimpleNtpMeta *	meta,
		std::set< int > &	gFlavors
	)

Definition at line 301 of file warp_gsimple.C.

References genie::flux::GSimpleNtpMeta::pdglist.

Referenced by fit_fcn(), and warp_gsimple().

 {
   // add any new flavors here, if necessary
   std::set<int>::const_iterator flitr = gFlavors.begin();
   std::vector<Int_t>& flvlist = meta->pdglist;
   for ( ; flitr != gFlavors.end(); ++flitr ) {
     int  seen_pdg = *flitr;
     std::vector<Int_t>::iterator knwitr = 
       std::find(flvlist.begin(),flvlist.end(),seen_pdg);
     bool known_pdg = ( knwitr != flvlist.end() );
     if ( ! known_pdg ) meta->pdglist.push_back(seen_pdg);
   }
 }

void warp_entry	(	GSimpleNtpEntry *	entry,
		GSimpleNtpAux *	aux,
		GSimpleNtpNuMI *	numi
	)

Definition at line 44 of file warp_gsimple.C.

References genie::flux::GSimpleNtpAux::auxint, om::cout, genie::flux::GSimpleNtpEntry::E, allTimeWatchdog::endl, genie::flux::GSimpleNtpNuMI::entryno, gDoWarp, pick_energy(), genie::flux::GSimpleNtpEntry::px, genie::flux::GSimpleNtpEntry::py, genie::flux::GSimpleNtpEntry::pz, and std::sqrt().

Referenced by warp_gsimple().

 {
   if ( ! gDoWarp ) return;
   // this is a dumb weighting scheme
   //if ( entry->E > 1.5 ) {
     //// don't set absolute weight in case we started with weighted entries
     //// scale whatever is there from GSimpleNtpFlux driver
     //entry->wgt *= 0.5;
   //}
   //// change some flavors
   //if ( gRandom->Rndm() < 0.25 ) {
     //entry->pdg = ( entry->pdg > 0 ) ? 16 : -16;
   //}
   
   // record the 4 momentum before warp
   TLorentzVector fmb4 = TLorentzVector(entry->px, entry->py, entry->pz,
                                        entry->E);
   // sample an energy value from the 1/E distribution
   // min energy 0.01 GeV, where GENIE spline starts
   // max energy 100 GeV, by eye
   entry->E = pick_energy(0.01, 50);
   // scale momentum so that the new momentum is in the same direction as the
   // old one but obeys the energy-momentum relation
   double k = sqrt(1+(entry->E*entry->E-fmb4.E()*fmb4.E())/fmb4.Vect().Mag2());
   entry->px *= k;
   entry->py *= k;
   entry->pz *= k;
 
   // use all branches so there are no complaints from compiler
   static double trash = 0;
   if ( numi ) trash = numi->entryno;
   if ( aux  ) trash = aux->auxint.size();
   if ( trash == -1 ) cout << "trash was -1" << endl;
 }

void warp_gsimple	(	string	fnameout = `"gsimple_output.root"`,
		string	fnamein = `"gsimple_input.root"`,
		bool	dowarp = `true`
	)

Definition at line 130 of file warp_gsimple.C.

 {
   
   gRandom->SetSeed(0);
   gDoWarp = dowarp;
 
   fcn = new TF1("fcn", fit_fcn, 0., 50, 6);
   // fcn->SetParameters(-8.10907e+01,-2.83834e+02,-3.21121e+00,1.24471e+04,9.42167e-01);
   fcn->SetParameters(-2.43864e+02, 4.42211e+02, 3.00555e+00, 1.23288e+04, 8.95128e-01, 2.38448e+00);
   fcn->SetNpx(100000);
 
   cout << "-----------" << endl;
   string fnameinlong(gSystem->ExpandPathName(fnamein.c_str()));
   cout << fnameinlong << endl;
   TFile* fin = TFile::Open(fnameinlong.c_str(),"READONLY");
   TTree* etreein = (TTree*)fin->Get("flux");
   TTree* mtreein = (TTree*)fin->Get("meta");
   genie::flux::GSimpleNtpEntry* entry_in = new genie::flux::GSimpleNtpEntry;
   genie::flux::GSimpleNtpNuMI*  numi_in  = new genie::flux::GSimpleNtpNuMI;
   genie::flux::GSimpleNtpAux*   aux_in   = new genie::flux::GSimpleNtpAux;
   genie::flux::GSimpleNtpMeta*  meta_in  = new genie::flux::GSimpleNtpMeta;
 
   long int nentries = etreein->GetEntries();
 
   int sba_status[4] = { -999, -999, -999, -999 };
   sba_status[0] = etreein->SetBranchAddress("entry",&entry_in);
   sba_status[1] = etreein->SetBranchAddress("numi",&numi_in);
   sba_status[2] = etreein->SetBranchAddress("aux",&aux_in);
   sba_status[3] = mtreein->SetBranchAddress("meta",&meta_in);
   cout << "SetBranchAddress results "
        << sba_status[0] << "," 
        << sba_status[1] << "," 
        << sba_status[2] << ","
        << sba_status[3] 
        << endl;
   bool donumi = ( sba_status[1] == 0 );
   bool doaux  = ( sba_status[2] == 0 );
 
   int nindices = mtreein->BuildIndex("metakey"); // tie metadata to entry
   cout << "saw " << nindices << " metakey indices" << endl;
 
   cout << "Creating:    " << fnameout << endl;
   cout << "Input file:  " << fnamein << endl;
   cout << "Branches:    entry " 
        << ( (doaux) ? "aux ":"" ) 
        << ( (donumi) ? "numi ":"" ) 
        << endl;
   if ( ! gDoWarp ) cout << "++++++++ NO ACTUAL WARP APPLIED" << endl;
 
   TFile* fout = TFile::Open(fnameout.c_str(),"RECREATE");
   TTree* fluxntp = new TTree("flux","a simple flux n-tuple");
   TTree* metantp = new TTree("meta","metadata for flux n-tuple");
   genie::flux::GSimpleNtpEntry* fentry = new genie::flux::GSimpleNtpEntry;
   genie::flux::GSimpleNtpAux*   faux   = new genie::flux::GSimpleNtpAux;
   genie::flux::GSimpleNtpNuMI*  fnumi  = new genie::flux::GSimpleNtpNuMI;
   genie::flux::GSimpleNtpMeta*  fmeta  = new genie::flux::GSimpleNtpMeta;
 
   fluxntp->Branch("entry",&fentry);
   if ( doaux  ) fluxntp->Branch("aux",&faux);
   if ( donumi ) fluxntp->Branch("numi",&fnumi);
   metantp->Branch("meta",&fmeta);
 
   cout << "=========================== Start " << endl;
 
   TStopwatch sw;
   sw.Start();
 
   const double large = 1.0e300;
   double minwgt = large, maxwgt = -large, maxenergy = -large;
   std::set<int> gFlavors;
 
   for ( long int ientry = 0; ientry < nentries; ++ientry ) {
 
     // reset what's been read in anticipation of a new entry
     entry_in->Reset();
     aux_in->Reset();
     numi_in->Reset();
     
     // read the next entry, get metadata if it's different
     //int nbytes = 
     etreein->GetEntry(ientry);
     UInt_t metakey = entry_in->metakey;
     if ( fmeta->metakey != metakey ) {
       // UInt_t oldkey = meta_in->metakey;
       int nbmeta = mtreein->GetEntryWithIndex(metakey);
       cout << "on entry " << ientry << " fetched metadata w/ key "
            << metakey << " read " << nbmeta << " bytes" << endl;
     }
 
     // reset what we're writing
     fentry->Reset();
     fnumi->Reset();
     faux->Reset();
 
     // copy read in objects to output objects
     *fentry = *entry_in;
     *fnumi  = *numi_in;
     *faux   = *aux_in;
 
     // mess with the values to your hearts content
     warp_entry(fentry,faux,fnumi);
     // keep track of any weights for the meta data
     if ( fentry->wgt < minwgt ) minwgt = fentry->wgt;
     if ( fentry->wgt > maxwgt ) maxwgt = fentry->wgt;
     // user might have changed an energy larger than any in original file
     if ( fentry->E > maxenergy ) maxenergy = fentry->E;
     // user might have changed a flavor keep track of all we see
     gFlavors.insert(fentry->pdg);
 
     // process currently held metadata after transition to metadata key
     if ( fmeta->metakey != meta_in->metakey ) {
       // new meta data found
       cout << "new meta found " << fmeta->metakey << "/" << meta_in->metakey << endl;
       if ( fmeta->metakey != 0 ) {
         // have metadata that needs adjustment and writing
         // before processing new metadata
         warp_meta(fmeta,fnamein);
         fmeta->minWgt = minwgt;
         fmeta->maxWgt = maxwgt;
         fmeta->maxEnergy = maxenergy;
         update_flavors(fmeta,gFlavors);
         gFlavors.clear();
         cout << "metantp->Fill() " << *fmeta << endl;
         metantp->Fill();
         // next meta-data restarts weight range (and energy max)
         minwgt    =  large;
         maxwgt    = -large;
         maxenergy = -large;
       }
       // get a copy of metadata for accumulating adjustments
       *fmeta = *meta_in;
       //cout << " copy meta_in " << *meta_in << endl << *fmeta << endl;
     }
  
     fluxntp->Fill();
   }
 
   // write last set of meta-data after all is done
   if ( fmeta->metakey != 0 ) {
     cout << "final meta flush " << fmeta->metakey << endl;
     // have metadata that needs adjustment and writing
     warp_meta(fmeta,fnamein);
     fmeta->minWgt = minwgt;
     fmeta->maxWgt = maxwgt;
     fmeta->maxEnergy = maxenergy;
     update_flavors(fmeta,gFlavors);
     gFlavors.clear();
     cout << "metantp->Fill() " << *fmeta << endl;
     metantp->Fill();
   }
 
   cout << "=========================== Complete " << endl;
 
   sw.Stop();
   cout << "Generated " << nentries 
        << endl
        << "Time to generate: " << endl;
   sw.Print();
   
   fout->cd();
 
   // write ntuples out
   fluxntp->Write();
   metantp->Write();
   fout->Close();
 
   cout << endl << endl;
 }

void warp_meta	(	GSimpleNtpMeta *	meta,
		string	fnamein
	)

Definition at line 81 of file warp_gsimple.C.

References gDoWarp, genie::flux::GSimpleNtpMeta::infiles, and lem_server::msg.

Referenced by warp_gsimple().

 {
   meta->infiles.push_back("THIS IS MDC WARPED FLUX FROM:");
   meta->infiles.push_back(fnamein);
   if ( ! gDoWarp ) {
     meta->infiles.push_back("NO ACTUAL WARPING APPLIED");
   } else {
     string msg ="USING WARP FUNCTION CODENAMED WHISKEY-TANGO-FOXTROT";
     meta->infiles.push_back(msg);
   }
 } 

Variable Documentation

TF1* fcn

Definition at line 35 of file warp_gsimple.C.

Referenced by art::EndPathExecutor::doForAllEnabledOutputWorkers_(), art::EndPathExecutor::doForAllEnabledWorkers_(), format_product_node(), pick_energy(), and warp_gsimple().

bool gDoWarp = true

Definition at line 39 of file warp_gsimple.C.

Referenced by warp_entry(), warp_gsimple(), and warp_meta().

Functions

Variables

Function Documentation

Variable Documentation