177 #include "libxml/xmlmemory.h" 178 #include "libxml/parser.h" 185 #include <TChainElement.h> 187 #include <TStopwatch.h> 219 #ifdef GNUMI_TEST_XY_WGT 220 static genie::flux::xypartials gpartials;
223 using namespace genie;
237 std::vector<long int> GetIntVector(
std::string str);
241 void ParseParamSet(xmlDocPtr&, xmlNodePtr&);
242 void ParseBeamDir(xmlDocPtr&, xmlNodePtr&);
244 void ParseRotSeries(xmlDocPtr&, xmlNodePtr&);
245 void ParseWindowSeries(xmlDocPtr&, xmlNodePtr&);
257 TVector3 fFluxWindowPtXML[3];
265 const TLorentzVector kPosCenterNearBeam(0.,0., 1039.35,0.);
304 LOG(
"Flux",
pNOTICE) <<
"GenerateNext signaled End() ";
309 bool nextok = this->GenerateNext_weighted();
310 if ( fGenWeighted )
return nextok;
311 if ( ! nextok )
continue;
326 double f = this->
Weight() / fMaxWeight;
330 fMaxWeight = this->
Weight() * fMaxWgtFudge;
332 <<
"** Fractional weight = " << f
333 <<
" > 1 !! Bump fMaxWeight estimate to " << fMaxWeight
334 << PassThroughInfo();
336 double r = (f < 1.) ? rnd->
RndFlux().Rndm() : 0;
337 bool accept = ( r <
f );
340 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 342 <<
"Generated beam neutrino: " 343 <<
"\n pdg-code: " << fCurEntry->fgPdgC
366 if ( ! fG3NuMI && ! fG4NuMI && ! fFlugg ) {
368 <<
"The flux driver has not been properly configured";
377 if ( fIUse < fNUse && fIEntry >= 0 ) {
382 this->ResetCurrent();
385 if ( fIEntry >= fNEntries ) {
388 if ( fICycle < fNCycles || fNCycles == 0 ) {
393 <<
"No more entries in input flux neutrino ntuple, cycle " 394 << fICycle <<
" of " << fNCycles;
402 fG3NuMI->GetEntry(fIEntry);
403 fCurEntry->MakeCopy(fG3NuMI);
404 }
else if ( fG4NuMI ) {
405 fG4NuMI->GetEntry(fIEntry);
406 fCurEntry->MakeCopy(fG4NuMI);
407 }
else if ( fFlugg ) {
408 fFlugg->GetEntry(fIEntry);
409 fCurEntry->MakeCopy(fFlugg);
411 LOG(
"Flux",
pERROR) <<
"No ntuple configured";
416 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 418 <<
"got " << fNNeutrinos <<
" new fIEntry " << fIEntry
419 <<
" evtno " << fCurEntry->evtno;
423 fCurEntry->pcodes = 0;
424 fCurEntry->units = 0;
428 fCurEntry->ConvertPartCodes();
430 fCurEntry->fgPdgC = fCurEntry->ntype;
445 fAccumPOTs += fEffPOTsPerNu / fMaxWeight;
448 if ( ! fPdgCList->ExistsInPDGCodeList(fCurEntry->fgPdgC) ) {
452 int badpdg = fCurEntry->fgPdgC;
453 if ( ! fPdgCListRej->ExistsInPDGCodeList(badpdg) ) {
454 fPdgCListRej->push_back(badpdg);
456 <<
"Encountered neutrino specie (" << badpdg
457 <<
" pcodes=" << fCurEntry->pcodes <<
")" 458 <<
" that wasn't in SetFluxParticles() list, " 459 <<
"\nDeclared list of neutrino species: " << *fPdgCList;
474 fCurEntry->fgX4 = fFluxWindowBase;
477 double& wgt_xy = fCurEntry->fgXYWgt;
478 switch ( fUseFluxAtDetCenter ) {
480 wgt_xy = fCurEntry->nwtnear;
481 Ev = fCurEntry->nenergyn;
484 wgt_xy = fCurEntry->nwtfar;
485 Ev = fCurEntry->nenergyf;
489 fCurEntry->fgX4 += ( rnd->
RndFlux().Rndm()*fFluxWindowDir1 +
490 rnd->
RndFlux().Rndm()*fFluxWindowDir2 );
491 fCurEntry->CalcEnuWgt(fCurEntry->fgX4,Ev,wgt_xy);
497 <<
"Flux neutrino energy exceeds declared maximum neutrino energy" 498 <<
"\nEv = " << Ev <<
"(> Ev{max} = " << fMaxEv <<
")";
503 fgX4dkvtx = TLorentzVector( fCurEntry->vx,
508 TVector3 dirNu = (fCurEntry->fgX4.Vect() - fgX4dkvtx.Vect()).
Unit();
509 fCurEntry->fgP4.SetPxPyPzE( Ev*dirNu.X(),
515 fWeight = fCurEntry->nimpwt * fCurEntry->fgXYWgt;
516 if ( fApplyTiltWeight ) {
517 double tiltwgt = dirNu.Dot( FluxWindowNormal() );
518 fWeight *= TMath::Abs( tiltwgt );
522 fSumWeight += this->
Weight();
525 Beam2UserP4(fCurEntry->fgP4,fCurEntry->fgP4User);
526 Beam2UserPos(fCurEntry->fgX4,fCurEntry->fgX4User);
529 if ( TMath::Abs(fZ0) < 1.0e30 ) this->MoveToZ0(fZ0);
531 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 533 <<
"Generated neutrino: " << fIEntry <<
" " << fCurEntry->evtno
534 <<
" nenergyn " << fCurEntry->nenergyn
535 <<
"\n pdg-code: " << fCurEntry->fgPdgC
543 <<
"Generated neutrino had E_nu = " << Ev <<
" > " << fMaxEv
555 TVector3 x3diff = fCurEntry->fgX4.Vect() - fgX4dkvtx.Vect();
556 return x3diff.Mag() * fLengthScaleB2U;
564 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 566 <<
"MoveToZ0 (z0usr=" << z0usr <<
") before:" 571 double pzusr = fCurEntry->fgP4User.Pz();
572 if ( TMath::Abs(pzusr) < 1.0
e-30 ) {
575 <<
"MoveToZ0(" << z0usr <<
") not possible due to pz_usr (" << pzusr <<
")";
579 double scale = (z0usr - fCurEntry->fgX4User.Z()) / pzusr;
580 fCurEntry->fgX4User += (scale*fCurEntry->fgP4User);
581 fCurEntry->fgX4 += ((fLengthScaleU2B*
scale)*fCurEntry->fgP4);
583 fCurEntry->fgX4.SetT(0);
584 fCurEntry->fgX4User.SetT(0);
586 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 588 <<
"MoveToZ0 (" << z0usr <<
") after:" 599 if (!fNuFluxTree)
return;
605 const double kRDET = 100.0;
606 const double kRDET2 = kRDET * kRDET;
607 double flux_area = fFluxWindowDir1.Vect().Cross(fFluxWindowDir2.Vect()).
Mag();
608 LOG(
"Flux",
pNOTICE) <<
"in CalcEffPOTsPerNu, area = " << flux_area;
610 if ( flux_area < 1.0
e-30 ) {
612 <<
"CalcEffPOTsPerNu called with flux window area effectively zero";
615 double area_ratio = TMath::Pi() * kRDET2 / flux_area;
616 fEffPOTsPerNu = area_ratio * ( (double)fFilePOTs / (
double)fNEntries );
626 <<
"The flux driver has not been properly configured";
646 bool found_cfg = this->LoadConfig(config);
649 <<
"LoadBeamSimData could not find XML config \"" << config <<
"\"\n";
653 fNuFluxFilePatterns = patterns;
654 std::vector<int> nfiles_from_pattern;
658 std::set<std::string> fnames;
660 LOG(
"Flux",
pINFO) <<
"LoadBeamSimData was passed " << patterns.size()
663 for (
size_t ipatt = 0; ipatt < patterns.size(); ++ipatt ) {
664 string pattern = patterns[ipatt];
665 nfiles_from_pattern.push_back(0);
668 <<
"Loading gnumi flux tree from ROOT file pattern [" 669 << std::setw(3) << ipatt <<
"] \"" << pattern <<
"\"";
672 string dirname = gSystem->UnixPathName(gSystem->WorkingDirectory());
673 size_t slashpos = pattern.find_last_of(
"/");
675 if ( slashpos != std::string::npos ) {
676 dirname = pattern.substr(0,slashpos);
677 LOG(
"Flux",
pINFO) <<
"Look for flux using directory " << dirname;
678 fbegin = slashpos + 1;
679 }
else { fbegin = 0; }
681 void* dirp = gSystem->OpenDirectory(gSystem->ExpandPathName(dirname.c_str()));
684 pattern.substr(fbegin,pattern.size()-fbegin);
685 TRegexp re(basename.c_str(),kTRUE);
687 while ( ( onefile = gSystem->GetDirEntry(dirp) ) ) {
688 TString afile = onefile;
689 if ( afile==
"." || afile==
".." )
continue;
690 if ( basename!=afile && afile.Index(re) == kNPOS )
continue;
691 std::string fullname = dirname +
"/" + afile.Data();
692 fnames.insert(fullname);
693 nfiles_from_pattern[ipatt]++;
695 gSystem->FreeDirectory(dirp);
700 std::set<string>::const_iterator sitr = fnames.begin();
701 for ( ; sitr != fnames.end(); ++sitr, ++indx ) {
714 TFile*
tf = TFile::Open(filename.c_str(),
"READ");
715 int isflugg = ( filename.find(
"flugg") != string::npos ) ? 1 : 0;
717 const std::string gnames[] = {
"g3numi",
"g4numi",
"flugg",
"g4flugg"};
718 for (
int j = 0;
j < 2 ; ++
j ) {
719 TTree* atree = (TTree*)tf->Get(tnames[
j].c_str());
724 if ( ! fNuFluxTree ) {
725 this->SetTreeName(tname_this);
726 fNuFluxTree =
new TChain(fNuFluxTreeName.c_str());
727 fNuFluxGen = gname_this;
732 if ( fNuFluxTreeName != tname_this ||
733 fNuFluxGen != gname_this ) {
735 <<
"Inconsistent flux file types\n" 736 <<
"The input gnumi flux file \"" << filename
737 <<
"\"\ncontains a '" << tname_this <<
"' " << gname_this
739 <<
"but a '" << fNuFluxTreeName <<
"' " << fNuFluxGen
740 <<
" numi ntuple has alread been seen in the chain";
752 if ( fNuFluxTreeName ==
"" ) {
754 <<
"The input gnumi flux file doesn't exist! Initialization failed!";
757 if ( fNuFluxGen ==
"g3numi" ) fG3NuMI =
new g3numi(fNuFluxTree);
758 if ( fNuFluxGen ==
"g4numi" ) fG4NuMI =
new g4numi(fNuFluxTree);
759 if ( fNuFluxGen ==
"flugg" ) fFlugg =
new flugg(fNuFluxTree);
762 fNEntries = fNuFluxTree->GetEntries();
764 if ( fNEntries == 0 ) {
766 <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
768 <<
"Loaded flux tree contains " << fNEntries <<
" entries";
770 <<
"Was passed the file patterns: ";
771 for (
size_t ipatt = 0; ipatt < patterns.size(); ++ipatt ) {
772 string pattern = patterns[ipatt];
774 <<
" [" << std::setw(3) << ipatt <<
"] " <<
pattern;
777 <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
780 <<
"Loaded flux tree contains " << fNEntries <<
" entries" 781 <<
" from " << fnames.size() <<
" unique files";
782 for (
size_t ipatt = 0; ipatt < patterns.size(); ++ipatt ) {
783 string pattern = patterns[ipatt];
785 <<
" pattern: " << pattern <<
" yielded " 786 << nfiles_from_pattern[ipatt] <<
" files";
793 <<
"LoadBeamSimData left detector location unset";
797 <<
"Run ScanForMaxWeight() as part of LoadBeamSimData";
798 this->ScanForMaxWeight();
807 fIEntry = rnd->
RndFlux().Integer(fNEntries) - 1;
814 LOG(
"Flux",
pNOTICE) <<
"about to CalcEffPOTsPerNu";
815 this->CalcEffPOTsPerNu();
820 std::vector<std::string>& branchClassNames,
821 std::vector<void**>& branchObjPointers)
828 branchNames.push_back(
"flux");
829 branchClassNames.push_back(
"genie::flux::GNuMIFluxPassThroughInfo");
830 branchObjPointers.push_back((
void**)&fCurEntry);
840 <<
"Specify a detector location before scanning for max weight";
845 int ipos_estimator = fUseFluxAtDetCenter;
846 if ( ipos_estimator == 0 ) {
848 double zbase = fFluxWindowBase.Z();
849 if ( TMath::Abs(zbase-103648.837) < 10000. ) ipos_estimator = -1;
850 if ( TMath::Abs(zbase-73534000. ) < 10000. ) ipos_estimator = +1;
852 if ( ipos_estimator != 0 ) {
856 const char* ntwgtstrv[4] = {
"Nimpwt*Nwtnear",
859 "Nimpwt*NWtFar[0]" };
861 if ( ipos_estimator > 0 ) strindx = 1;
862 if ( fG4NuMI ) strindx += 2;
864 Long64_t nscan = TMath::Min(fNEntries,200000LL);
866 fNuFluxTree->Draw(ntwgtstrv[strindx],
"",
"goff",nscan);
871 Long64_t
idx = TMath::LocMax(fNuFluxTree->GetSelectedRows(),
872 fNuFluxTree->GetV1());
874 fMaxWeight = fNuFluxTree->GetV1()[
idx];
875 LOG(
"Flux",
pNOTICE) <<
"Maximum flux weight from Nwt in ntuple = " 877 if ( fMaxWeight <= 0 ) {
878 LOG(
"Flux",
pFATAL) <<
"Non-positive maximum flux weight!";
884 double wgtgenmx = 0, enumx = 0;
887 for (
int itry=0; itry < fMaxWgtEntries; ++itry) {
888 this->GenerateNext_weighted();
890 if ( wgt > wgtgenmx ) wgtgenmx =
wgt;
891 double enu = fCurEntry->fgP4.Energy();
892 if ( enu > enumx ) enumx =
enu;
896 LOG(
"Flux",
pNOTICE) <<
"Maximum flux weight for spin = " 897 << wgtgenmx <<
", energy = " << enumx
898 <<
" (" << fMaxWgtEntries <<
")";
900 if (wgtgenmx > fMaxWeight ) fMaxWeight = wgtgenmx;
902 fMaxWeight *= fMaxWgtFudge;
904 if ( enumx*fMaxEFudge > fMaxEv ) {
905 LOG(
"Flux",
pNOTICE) <<
"Adjust max: was=" << fMaxEv
906 <<
" now " << enumx <<
"*" << fMaxEFudge
907 <<
" = " << enumx*fMaxEFudge;
908 fMaxEv = enumx * fMaxEFudge;
911 LOG(
"Flux",
pNOTICE) <<
"Maximum flux weight = " << fMaxWeight
912 <<
", energy = " << fMaxEv;
918 fMaxEv = TMath::Max(0.,Ev);
921 <<
"Declared maximum flux neutrino energy: " << fMaxEv;
929 fNUse = TMath::Max(1
L, nuse);
934 fNuFluxTreeName =
name;
940 fFluxWindowBase = kPosCenterNearBeam;
941 fFluxWindowDir1 = TLorentzVector();
942 fFluxWindowDir2 = TLorentzVector();
943 fUseFluxAtDetCenter = -1;
951 fFluxWindowDir1 = TLorentzVector();
952 fFluxWindowDir2 = TLorentzVector();
953 fUseFluxAtDetCenter = +1;
962 double padbeam = padding / fLengthScaleB2U;
964 <<
"SetBeamFluxWindow " << (
int)stdwindow <<
" padding " << padbeam <<
" cm";
967 switch ( stdwindow ) {
968 #ifdef THIS_IS_NOT_YET_IMPLEMENTED 970 SetBeamFluxWindow(103648.837);
973 SetBeamFluxWindow(73534000.);
982 case kMinosNearCenter:
985 fFluxWindowBase = kPosCenterNearBeam;
986 fFluxWindowDir1 = TLorentzVector();
987 fFluxWindowDir2 = TLorentzVector();
988 TLorentzVector usrpos;
989 Beam2UserPos(fFluxWindowBase, usrpos);
990 fFluxWindowPtUser[0] = usrpos.Vect();
991 fFluxWindowPtUser[1] = fFluxWindowPtUser[0];
992 fFluxWindowPtUser[2] = fFluxWindowPtUser[0];
997 case kMinosFarCenter:
1001 fFluxWindowDir1 = TLorentzVector();
1002 fFluxWindowDir2 = TLorentzVector();
1003 TLorentzVector usrpos;
1004 Beam2UserPos(fFluxWindowBase, usrpos);
1005 fFluxWindowPtUser[0] = usrpos.Vect();
1006 fFluxWindowPtUser[1] = fFluxWindowPtUser[0];
1007 fFluxWindowPtUser[2] = fFluxWindowPtUser[0];
1008 fFluxWindowLen1 = 0;
1009 fFluxWindowLen2 = 0;
1014 <<
"SetBeamFluxWindow - StdFluxWindow " << stdwindow
1015 <<
" not yet implemented";
1018 LOG(
"Flux",
pNOTICE) <<
"about to CalcEffPOTsPerNu";
1019 this->CalcEffPOTsPerNu();
1028 fUseFluxAtDetCenter = 0;
1029 fDetLocIsSet =
true;
1031 fFluxWindowPtUser[0] =
p0;
1032 fFluxWindowPtUser[1] =
p1;
1033 fFluxWindowPtUser[2] =
p2;
1037 TLorentzVector ptbm0, ptbm1, ptbm2;
1038 User2BeamPos(TLorentzVector(fFluxWindowPtUser[0],0),ptbm0);
1039 User2BeamPos(TLorentzVector(fFluxWindowPtUser[1],0),ptbm1);
1040 User2BeamPos(TLorentzVector(fFluxWindowPtUser[2],0),ptbm2);
1042 fFluxWindowBase = ptbm0;
1043 fFluxWindowDir1 = ptbm1 - ptbm0;
1044 fFluxWindowDir2 = ptbm2 - ptbm0;
1046 fFluxWindowLen1 = fFluxWindowDir1.Mag();
1047 fFluxWindowLen2 = fFluxWindowDir2.Mag();
1048 fWindowNormal = fFluxWindowDir1.Vect().Cross(fFluxWindowDir2.Vect()).
Unit();
1050 double dot = fFluxWindowDir1.Dot(fFluxWindowDir2);
1051 if ( TMath::Abs(dot) > 1.0
e-8 )
1052 LOG(
"Flux",
pWARN) <<
"Dot product between window direction vectors was " 1053 << dot <<
"; please check for orthoganality";
1055 LOG(
"Flux",
pNOTICE) <<
"about to CalcEffPOTsPerNu";
1056 this->CalcEffPOTsPerNu();
1063 p0 = fFluxWindowPtUser[0];
1064 p1 = fFluxWindowPtUser[1];
1065 p2 = fFluxWindowPtUser[2];
1072 fBeamRot = TLorentzRotation(beamrot);
1073 fBeamRotInv = fBeamRot.Inverse();
1080 fBeamZero = TLorentzVector(beam0,0);
1090 const TLorentzRotation& rot4 = fBeamRot;
1091 TVector3 newX(rot4.XX(),rot4.XY(),rot4.XZ());
1092 TVector3 newY(rot4.YX(),rot4.YY(),rot4.YZ());
1093 TVector3 newZ(rot4.ZX(),rot4.ZY(),rot4.ZZ());
1094 rot3.RotateAxes(newX,newY,newZ);
1095 return rot3.Inverse();
1099 TVector3 beam0 = fBeamZero.Vect();
1128 TLorentzVector& usrxyz)
const 1130 usrxyz = fLengthScaleB2U*(fBeamRot*beamxyz) + fBeamZero;
1133 TLorentzVector& usrdir)
const 1135 usrdir = fLengthScaleB2U*(fBeamRot*beamdir);
1138 TLorentzVector& usrp4 )
const 1140 usrp4 = fBeamRot*beamp4;
1144 TLorentzVector& beamxyz)
const 1146 beamxyz = fLengthScaleU2B*(fBeamRotInv*(usrxyz-fBeamZero));
1149 TLorentzVector& beamdir)
const 1151 beamdir = fLengthScaleU2B*(fBeamRotInv*usrdir);
1154 TLorentzVector& beamp4)
const 1156 beamp4 = fBeamRotInv*usrp4;
1162 LOG(
"Flux",
pNOTICE) <<
"CurrentEntry:" << *fCurEntry;
1169 LOG(
"Flux",
pWARN) <<
"GSimpleNtpFlux::Clear(" << opt <<
") called";
1183 fGenWeighted = gen_weighted;
1188 LOG(
"Flux",
pNOTICE) <<
"Initializing GNuMIFlux driver";
1199 fNuFluxTreeName =
"";
1213 fMaxWgtFudge = 1.05;
1214 fMaxWgtEntries = 2500000;
1222 fGenWeighted =
false;
1223 fApplyTiltWeight =
true;
1224 fUseFluxAtDetCenter = 0;
1225 fDetLocIsSet =
false;
1229 this->SetDefaults();
1230 this->ResetCurrent();
1245 LOG(
"Flux",
pNOTICE) <<
"Setting default GNuMIFlux driver options";
1253 this->SetFluxParticles (particles);
1254 this->SetMaxEnergy (120.);
1255 this->SetUpstreamZ (-3.4e38);
1256 this->SetNumOfCycles (0);
1257 this->SetEntryReuse (1);
1259 this->SetXMLFileBase(
"GNuMIFlux.xml");
1267 fCurEntry->ResetCurrent();
1268 fCurEntry->ResetCopy();
1275 if (fPdgCList)
delete fPdgCList;
1276 if (fPdgCListRej)
delete fPdgCListRej;
1277 if (fCurEntry)
delete fCurEntry;
1279 if ( fG3NuMI )
delete fG3NuMI;
1280 if ( fG4NuMI )
delete fG4NuMI;
1281 if ( fFlugg )
delete fFlugg;
1284 <<
" flux file cycles: " << fICycle <<
" of " << fNCycles
1285 <<
", entry " << fIEntry <<
" use: " << fIUse <<
" of " << fNUse;
1293 ULong64_t
nentries = thetree->GetEntries();
1304 thetree->SetMakeClass(1);
1310 TBranch* br_evtno = 0;
1311 thetree->SetBranchAddress(
"evtno",&evtno, &br_evtno);
1312 Int_t evt_1 = 0x7FFFFFFF;
1316 for (
int j=0;
j<50; ++
j) {
1317 thetree->GetEntry(
j);
1318 if (evtno != 0) evt_1 = TMath::Min(evtno,evt_1);
1319 thetree->GetEntry(nentries-1 -
j );
1320 if (evtno != 0) evt_N = TMath::Max(evtno,evt_N);
1328 const Int_t
nquant = 1000;
1329 const Double_t rquant =
nquant;
1331 for (
int j=0;
j<50; ++
j) {
1332 thetree->GetEntry(
j);
1333 if (evtno != 0) evt_1 = TMath::Min(evtno,evt_1);
1336 for (
int j=0;
j<50; ++
j) {
1337 thetree->GetEntry(nentries-1 -
j );
1338 if (evtno != 0) evt_N = TMath::Max(evtno,evt_N);
1339 std::cout <<
"[" << (nentries-1-
j) <<
"] evtno=" << evtno <<
" evt_N=" << evt_N <<
std::endl;
1342 Int_t nquant = 1000;
1343 Double_t rquant =
nquant;
1346 Int_t est_1 = (TMath::FloorNint(evt_1/rquant))*nquant + 1;
1347 Int_t est_N = (TMath::FloorNint((evt_N-1)/rquant)+1)*nquant;
1348 ULong64_t npots = est_N - est_1 + 1;
1350 << fNuFluxTreeName <<
"->AddFile() of " << nentries <<
" entries [" 1351 << evt_1 <<
":" << evt_N <<
"%" << nquant
1352 <<
"(" << est_1 <<
":" << est_N <<
")=" 1353 << npots <<
" POTs] in {" << fNuFluxGen <<
"} file: " <<
fname;
1358 fNuFluxTree->AddFile(fname.c_str());
1374 double rescale = fLengthUnits / user_units;
1375 fLengthUnits = user_units;
1377 fLengthScaleB2U = cm / user_units;
1378 fLengthScaleU2B = user_units /
cm;
1382 fCurEntry->fgX4User *= rescale;
1383 fBeamZero *= rescale;
1384 fFluxWindowPtUser[0] *= rescale;
1385 fFluxWindowPtUser[1] *= rescale;
1386 fFluxWindowPtUser[2] *= rescale;
1400 return fLengthScaleB2U *
cm ;
1480 #ifndef SKIP_MINERVA_MODS 1487 for (
unsigned int itclear = 0; itclear <
MAX_N_TRAJ; itclear++ ) {
1520 fgP4.SetPxPyPzE(0.,0.,0.,0.);
1521 fgX4.SetXYZT(0.,0.,0.,0.);
1612 <<
"ConvertPartCodes saw ntype " <<
ntype <<
" -- unknown ";
1617 }
else if (
pcodes != 1 ) {
1620 <<
"ConvertPartCodes saw pcodes flag as " <<
pcodes;
1705 const int kNearIndx = 0;
1706 const int kFarIndx = 0;
1773 #ifndef SKIP_MINERVA_MODS
1784 for ( Int_t ic = 0; ic < ntraj; ++ic ) {
1886 double&
enu,
double& wgt_xy)
const 1931 const double kPIMASS = 0.13957;
1932 const double kKMASS = 0.49368;
1933 const double kK0MASS = 0.49767;
1934 const double kMUMASS = 0.105658389;
1935 const double kOMEGAMASS = 1.67245;
1937 const int kpdg_nue = 12;
1938 const int kpdg_nuebar = -12;
1939 const int kpdg_numu = 14;
1940 const int kpdg_numubar = -14;
1942 const int kpdg_muplus = -13;
1943 const int kpdg_muminus = 13;
1944 const int kpdg_pionplus = 211;
1945 const int kpdg_pionminus = -211;
1946 const int kpdg_k0long = 130;
1947 const int kpdg_k0short = 310;
1948 const int kpdg_k0mix = 311;
1949 const int kpdg_kaonplus = 321;
1950 const int kpdg_kaonminus = -321;
1951 const int kpdg_omegaminus = 3334;
1952 const int kpdg_omegaplus = -3334;
1954 const double kRDET = 100.0;
1956 double xpos = xyz.X();
1957 double ypos = xyz.Y();
1958 double zpos = xyz.Z();
1966 double parent_mass = kPIMASS;
1967 switch ( this->
ptype ) {
1969 case kpdg_pionminus:
1970 parent_mass = kPIMASS;
1973 case kpdg_kaonminus:
1974 parent_mass = kKMASS;
1979 parent_mass = kK0MASS;
1983 parent_mass = kMUMASS;
1985 case kpdg_omegaminus:
1986 case kpdg_omegaplus:
1987 parent_mass = kOMEGAMASS;
1992 LOG(
"Flux",
pFATAL) <<
"NU_REWGT unknown particle type " << this->
ptype;
1997 double parentp2 = ( this->
pdpx*this->
pdpx +
2001 parent_mass*parent_mass);
2004 double gamma = parent_energy / parent_mass;
2005 double gamma_sqr = gamma * gamma;
2006 double beta_mag =
TMath::Sqrt( ( gamma_sqr - 1.0 )/gamma_sqr );
2009 double enuzr = this->
necm;
2012 (ypos-this->
vy)*(ypos-this->
vy) +
2013 (zpos-this->
vz)*(zpos-this->
vz) );
2016 double costh_pardet = -999., theta_pardet = -999.;
2019 if ( parentp > 0. ) {
2020 costh_pardet = ( this->
pdpx*(xpos-this->
vx) +
2021 this->
pdpy*(ypos-this->
vy) +
2022 this->
pdpz*(zpos-this->
vz) )
2024 if ( costh_pardet > 1.0 ) costh_pardet = 1.0;
2025 if ( costh_pardet < -1.0 ) costh_pardet = -1.0;
2026 theta_pardet = TMath::ACos(costh_pardet);
2029 emrat = 1.0 / ( gamma * ( 1.0 - beta_mag * costh_pardet ));
2032 enu = emrat * enuzr;
2040 <<
" p " << parentp <<
" e " << parent_energy <<
" gamma " << gamma
2043 std::cout <<
" enuzr " << enuzr <<
" rad " << rad <<
" costh " << costh_pardet
2044 <<
" theta " << theta_pardet <<
" emrat " << emrat
2048 #ifdef GNUMI_TEST_XY_WGT 2049 gpartials.xdet = xpos;
2050 gpartials.ydet = ypos;
2051 gpartials.zdet = zpos;
2052 gpartials.parent_mass = parent_mass;
2053 gpartials.parentp = parentp;
2054 gpartials.parent_energy = parent_energy;
2055 gpartials.gamma = gamma;
2056 gpartials.beta_mag = beta_mag;
2057 gpartials.enuzr = enuzr;
2058 gpartials.rad =
rad;
2059 gpartials.costh_pardet = costh_pardet;
2060 gpartials.theta_pardet = theta_pardet;
2061 gpartials.emrat = emrat;
2062 gpartials.eneu =
enu;
2069 double sanddetcomp =
TMath::Sqrt(( (xpos-this->
vx)*(xpos-this->
vx) ) +
2070 ( (ypos-this->
vy)*(ypos-this->
vy) ) +
2071 ( (zpos-this->
vz)*(zpos-this->
vz) ) );
2072 double sangdet = ( 1.0 - TMath::Cos(TMath::ATan( kRDET / sanddetcomp)))/2.0;
2075 wgt_xy = sangdet * ( emrat * emrat );
2077 #ifdef GNUMI_TEST_XY_WGT 2078 gpartials.sangdet = sangdet;
2079 gpartials.wgt = wgt_xy;
2080 gpartials.ptype = this->
ptype;
2086 if ( this->
ptype == kpdg_muplus || this->
ptype == kpdg_muminus) {
2087 double beta[3], p_dcm_nu[4], p_nu[3], p_pcm_mp[3],
partial;
2090 beta[0] = this->
pdpx / parent_energy;
2091 beta[1] = this->
pdpy / parent_energy;
2092 beta[2] = this->
pdpz / parent_energy;
2093 p_nu[0] = (xpos-this->
vx)*enu/rad;
2094 p_nu[1] = (ypos-this->
vy)*enu/rad;
2095 p_nu[2] = (zpos-this->
vz)*enu/rad;
2097 (beta[0]*p_nu[0] + beta[1]*p_nu[1] + beta[2]*p_nu[2] );
2098 partial = enu - partial/(gamma+1.0);
2103 p_dcm_nu[0] = p_nu[0] - beta[0]*gamma*
partial;
2104 p_dcm_nu[1] = p_nu[1] - beta[1]*gamma*
partial;
2105 p_dcm_nu[2] = p_nu[2] - beta[2]*gamma*
partial;
2106 p_dcm_nu[3] =
TMath::Sqrt( p_dcm_nu[0]*p_dcm_nu[0] +
2107 p_dcm_nu[1]*p_dcm_nu[1] +
2108 p_dcm_nu[2]*p_dcm_nu[2] );
2110 #ifdef GNUMI_TEST_XY_WGT 2111 gpartials.betanu[0] = beta[0];
2112 gpartials.betanu[1] = beta[1];
2113 gpartials.betanu[2] = beta[2];
2114 gpartials.p_nu[0] = p_nu[0];
2115 gpartials.p_nu[1] = p_nu[1];
2116 gpartials.p_nu[2] = p_nu[2];
2118 gpartials.p_dcm_nu[0] = p_dcm_nu[0];
2119 gpartials.p_dcm_nu[1] = p_dcm_nu[1];
2120 gpartials.p_dcm_nu[2] = p_dcm_nu[2];
2121 gpartials.p_dcm_nu[3] = p_dcm_nu[3];
2125 double particle_energy = this->
ppenergy;
2126 gamma = particle_energy/parent_mass;
2127 beta[0] = this->
ppdxdz * this->
pppz / particle_energy;
2128 beta[1] = this->
ppdydz * this->
pppz / particle_energy;
2129 beta[2] = this->
pppz / particle_energy;
2130 partial = gamma * ( beta[0]*this->
muparpx +
2133 partial = this->
mupare - partial/(gamma+1.0);
2137 double p_pcm =
TMath::Sqrt ( p_pcm_mp[0]*p_pcm_mp[0] +
2138 p_pcm_mp[1]*p_pcm_mp[1] +
2139 p_pcm_mp[2]*p_pcm_mp[2] );
2148 #ifdef GNUMI_TEST_XY_WGT 2149 gpartials.muparent_px = this->
muparpx;
2150 gpartials.muparent_py = this->
muparpy;
2151 gpartials.muparent_pz = this->
muparpz;
2152 gpartials.gammamp = gamma;
2153 gpartials.betamp[0] = beta[0];
2154 gpartials.betamp[1] = beta[1];
2155 gpartials.betamp[2] = beta[2];
2157 gpartials.p_pcm_mp[0] = p_pcm_mp[0];
2158 gpartials.p_pcm_mp[1] = p_pcm_mp[1];
2159 gpartials.p_pcm_mp[2] = p_pcm_mp[2];
2160 gpartials.p_pcm = p_pcm;
2163 const double eps = 1.0e-30;
2164 if ( p_pcm < eps || p_dcm_nu[3] < eps ) {
2168 double costh = ( p_dcm_nu[0]*p_pcm_mp[0] +
2169 p_dcm_nu[1]*p_pcm_mp[1] +
2170 p_dcm_nu[2]*p_pcm_mp[2] ) /
2171 ( p_dcm_nu[3]*p_pcm );
2172 if ( costh > 1.0 ) costh = 1.0;
2173 if ( costh < -1.0 ) costh = -1.0;
2175 double wgt_ratio = 0.0;
2176 switch ( this->
ntype ) {
2179 wgt_ratio = 1.0 - costh;
2184 double xnu = 2.0 * enuzr / kMUMASS;
2185 wgt_ratio = ( (3.0-2.0*xnu ) - (1.0-2.0*xnu)*costh ) / (3.0-2.0*xnu);
2191 wgt_xy = wgt_xy * wgt_ratio;
2193 #ifdef GNUMI_TEST_XY_WGT 2194 gpartials.ntype = this->
ntype;
2195 gpartials.costhmu = costh;
2196 gpartials.wgt_ratio = wgt_ratio;
2213 stream <<
"\nGNuMIFlux run " << info.
run <<
" evtno " << info.
evtno 2214 <<
" (pcodes " << info.
pcodes <<
" units " << info.
units <<
")" 2215 <<
"\n random dk: dx/dz " << info.
ndxdz 2216 <<
" dy/dz " << info.
ndydz 2217 <<
" pz " << info.
npz <<
" E " << info.
nenergy 2218 <<
"\n near00 dk: dx/dz " << info.
ndxdznea 2221 <<
"\n far00 dk: dx/dz " << info.
ndxdzfar 2224 <<
"\n norig " << info.
norig <<
" ndecay " << info.
ndecay 2225 <<
" ntype " << info.
ntype 2226 <<
"\n had vtx " << info.
vx <<
" " << info.
vy <<
" " << info.
vz 2227 <<
"\n parent p3 @ dk " << info.
pdpx <<
" " << info.
pdpy <<
" " << info.
pdpz 2228 <<
"\n parent prod: dx/dz " << info.
ppdxdz 2229 <<
" dy/dz " << info.
ppdydz 2231 <<
"\n ppmedium " << info.
ppmedium <<
" ptype " << info.
ptype 2232 <<
" ppvtx " << info.
ppvx <<
" " << info.
ppvy <<
" " << info.
ppvz 2235 <<
"\n necm " << info.
necm <<
" nimpwt " << info.
nimpwt 2236 <<
"\n point x,y,z " << info.
xpoint <<
" " << info.
ypoint 2238 <<
"\n tv x,y,z " << info.
tvx <<
" " << info.
tvy <<
" " << info.
tvz 2239 <<
"\n tptype " << info.
tptype <<
" tgen " << info.
tgen 2240 <<
" tgptype " << info.
tgptype 2241 <<
"\n tgp px,py,pz " << info.
tgppx <<
" " << info.
tgppy 2242 <<
" " << info.
tgppz 2243 <<
"\n tpriv x,y,z " << info.
tprivx <<
" " << info.
tprivy 2245 <<
"\n beam x,y,z " << info.
beamx <<
" " << info.
beamy 2246 <<
" " << info.
beamz 2247 <<
"\n beam px,py,pz " << info.
beampx <<
" " << info.
beampy 2251 #ifndef SKIP_MINERVA_MODS 2255 stream <<
"\nNeutrino History : ntrajectories " << info.
ntrajectory 2256 <<
"\n (trkID, pdg) of nu parents: [";
2260 for ( Int_t itt = 0; itt < ntraj; ++itt )
2261 stream <<
"(" << info.
trackId[itt-1] <<
", " << info.
pdgcode[itt] <<
") ";
2266 stream <<
"\nCurrent: pdg " << info.
fgPdgC 2273 #ifdef GNUMI_TEST_XY_WGT 2274 stream <<
"\n" << xypartials::GetStaticInstance();
2305 #ifdef GNUMI_TEST_XY_WGT 2307 double fabserr(
double a,
double b)
2308 {
return TMath::Abs(a-b)/TMath::Max(TMath::Abs(b),1.0
e-30); }
2310 int outdiff(
double a,
double b,
double eps,
const char*
label)
2312 double err = fabserr(a,b);
2314 std::cout << std::setw(15) << label <<
" err " << err
2315 <<
" vals " << a <<
" " << b <<
std::endl;
2321 int gnumi2pdg(
int igeant)
2324 case 52:
return -12;
2326 case 55:
return -14;
2328 case 58:
return -16;
2335 void xypartials::ReadStream(std::ifstream& myfile)
2337 myfile >> parent_mass >> parentp >> parent_energy;
2338 myfile >> gamma >> beta_mag >> enuzr >>
rad;
2339 myfile >> costh_pardet >> theta_pardet >> emrat >> eneu;
2340 myfile >> sangdet >>
wgt;
2343 ptype = gnumi2pdg(ptype_g);
2346 myfile >> betanu[0] >> betanu[1] >> betanu[2]
2347 >> p_nu[0] >> p_nu[1] >> p_nu[2];
2349 >> p_dcm_nu[0] >> p_dcm_nu[1] >> p_dcm_nu[2] >> p_dcm_nu[3];
2351 myfile >> muparent_px >> muparent_py >> muparent_pz;
2352 myfile >> gammamp >> betamp[0] >> betamp[1] >> betamp[2];
2354 >> p_pcm_mp[0] >> p_pcm_mp[1] >> p_pcm_mp[2] >> p_pcm;
2356 if ( p_pcm != 0.0 && p_dcm_nu[3] != 0.0 ) {
2358 myfile >> costhmu >> ntype_g >> wgt_ratio;
2359 ntype = gnumi2pdg(ntype_g);
2366 double eps1 = 2.5e-5;
2367 double eps2 = 2.5e-5;
2368 double epsX = 2.5e-5;
2370 np += outdiff(parent_mass ,other.parent_mass ,eps1,
"parent_mass");
2371 np += outdiff(parentp ,other.parentp ,eps1,
"parentp");
2372 np += outdiff(parent_energy,other.parent_energy,eps1,
"parent_energy");
2373 np += outdiff(gamma ,other.gamma ,eps1,
"gamma");
2374 np += outdiff(beta_mag ,other.beta_mag ,eps1,
"beta_mag");
2375 np += outdiff(enuzr ,other.enuzr ,eps1,
"enuzr");
2376 np += outdiff(
rad ,other.rad ,eps1,
"rad");
2377 np += outdiff(costh_pardet ,other.costh_pardet ,eps1,
"costh_pardet");
2379 np += outdiff(emrat ,other.emrat ,eps1,
"emrat");
2380 np += outdiff(eneu ,other.eneu ,epsX,
"eneu");
2381 np += outdiff(sangdet ,other.sangdet ,eps1,
"sangdet");
2382 np += outdiff(
wgt ,other.wgt ,epsX,
"wgt");
2383 if (
ptype != other.ptype ) {
2387 if ( TMath::Abs(
ptype)==13 || TMath::Abs(other.ptype)==13 ) {
2389 np += outdiff(betanu[0] ,other.betanu[0] ,eps2,
"betanu[0]");
2390 np += outdiff(betanu[1] ,other.betanu[1] ,eps2,
"betanu[1]");
2391 np += outdiff(betanu[2] ,other.betanu[2] ,eps2,
"betanu[2]");
2392 np += outdiff(p_nu[0] ,other.p_nu[0] ,eps2,
"p_nu[0]");
2393 np += outdiff(p_nu[1] ,other.p_nu[1] ,eps2,
"p_nu[1]");
2394 np += outdiff(p_nu[2] ,other.p_nu[2] ,eps2,
"p_nu[2]");
2395 np += outdiff(partial1 ,other.partial1 ,eps2,
"partial1");
2396 np += outdiff(p_dcm_nu[0] ,other.p_dcm_nu[0] ,eps2,
"p_dcm_nu[0]");
2397 np += outdiff(p_dcm_nu[1] ,other.p_dcm_nu[1] ,eps2,
"p_dcm_nu[1]");
2398 np += outdiff(p_dcm_nu[2] ,other.p_dcm_nu[2] ,eps2,
"p_dcm_nu[2]");
2399 np += outdiff(p_dcm_nu[3] ,other.p_dcm_nu[3] ,eps2,
"p_dcm_nu[3]");
2401 np += outdiff(muparent_px ,other.muparent_px ,eps2,
"muparent_px");
2402 np += outdiff(muparent_py ,other.muparent_py ,eps2,
"muparent_py");
2403 np += outdiff(muparent_pz ,other.muparent_pz ,eps2,
"muparent_pz");
2404 np += outdiff(gammamp ,other.gammamp ,eps1,
"gammamp");
2405 np += outdiff(betamp[0] ,other.betamp[0] ,eps1,
"betamp[0]");
2406 np += outdiff(betamp[1] ,other.betamp[1] ,eps1,
"betamp[1]");
2407 np += outdiff(betamp[2] ,other.betamp[2] ,eps1,
"betamp[2]");
2408 np += outdiff(partial2 ,other.partial2 ,eps1,
"partial2");
2409 np += outdiff(p_pcm_mp[0] ,other.p_pcm_mp[0] ,eps1,
"p_pcm_mp[0]");
2410 np += outdiff(p_pcm_mp[1] ,other.p_pcm_mp[1] ,eps1,
"p_pcm_mp[1]");
2411 np += outdiff(p_pcm_mp[2] ,other.p_pcm_mp[2] ,eps1,
"p_pcm_mp[2]");
2412 np += outdiff(p_pcm ,other.p_pcm ,eps1,
"p_pcm");
2414 if (
ntype != other.ntype ) {
2418 np += outdiff(costhmu ,other.costhmu ,eps1,
"costhmu");
2419 np += outdiff(wgt_ratio ,other.wgt_ratio ,eps1,
"wgt_ratio");
2434 stream <<
"GNuMIFlux xypartials " <<
std::endl;
2435 stream <<
" xyzdet (" << xyp.xdet <<
"," << xyp.ydet <<
"," 2437 stream <<
" parent: mass=" << xyp.parent_mass <<
" p=" << xyp.parentp
2438 <<
" e=" << xyp.parent_energy <<
" gamma=" << xyp.gamma
2439 <<
" beta_mag=" << xyp.beta_mag <<
std::endl;
2440 stream <<
" enuzr=" << xyp.enuzr <<
" rad=" << xyp.rad
2441 <<
" costh_pardet=" << xyp.costh_pardet <<
std::endl;
2442 stream <<
" emrat=" << xyp.emrat <<
" sangdet=" << xyp.sangdet
2444 stream <<
" ptype=" << xyp.ptype <<
" " 2445 << ((TMath::Abs(xyp.ptype) == 13)?
"is-muon":
"not-muon")
2448 if ( TMath::Abs(xyp.ptype)==13 ) {
2449 stream <<
" betanu: [" << xyp.betanu[0] <<
"," << xyp.betanu[1]
2450 <<
"," << xyp.betanu[2] <<
"]" <<
std::endl;
2451 stream <<
" p_nu: [" << xyp.p_nu[0] <<
"," << xyp.p_nu[1]
2452 <<
"," << xyp.p_nu[2] <<
"]" <<
std::endl;
2453 stream <<
" partial1=" << xyp.partial1 <<
std::endl;
2454 stream <<
" p_dcm_nu: [" << xyp.p_dcm_nu[0] <<
"," << xyp.p_dcm_nu[1]
2455 <<
"," << xyp.p_dcm_nu[2]
2456 <<
"," << xyp.p_dcm_nu[3] <<
"]" <<
std::endl;
2457 stream <<
" muparent_p: [" << xyp.muparent_px <<
"," << xyp.muparent_py
2458 <<
"," << xyp.muparent_pz <<
"]" <<
std::endl;
2459 stream <<
" gammamp=" << xyp.gammamp <<
std::endl;
2460 stream <<
" betamp: [" << xyp.betamp[0] <<
"," << xyp.betamp[1] <<
"," 2462 stream <<
" partial2=" << xyp.partial2 <<
std::endl;
2463 stream <<
" p_pcm_mp: [" << xyp.p_pcm_mp[0] <<
"," << xyp.p_pcm_mp[1]
2464 <<
"," << xyp.p_pcm_mp[2] <<
"] p_pcm=" 2466 stream <<
" ntype=" << xyp.ntype
2467 <<
" costhmu=" << xyp.costhmu
2468 <<
" wgt_ratio=" << xyp.wgt_ratio <<
std::endl;
2475 xypartials& xypartials::GetStaticInstance()
2476 {
return gpartials; }
2482 const char* altxml = gSystem->Getenv(
"GNUMIFLUXXML");
2484 SetXMLFileBase(altxml);
2495 std::ostringstream
s;
2496 PDGCodeList::const_iterator itr = fPdgCList->begin();
2497 for ( ; itr != fPdgCList->end(); ++itr) s << (*itr) <<
" ";
2499 itr = fPdgCListRej->begin();
2500 for ( ; itr != fPdgCListRej->end(); ++itr) s << (*itr) <<
" ";
2503 std::ostringstream fpattout;
2504 for (
size_t i = 0;
i < fNuFluxFilePatterns.size(); ++
i)
2505 fpattout <<
"\n [" << std::setw(3) <<
i <<
"] " << fNuFluxFilePatterns[
i];
2507 std::ostringstream flistout;
2508 std::vector<std::string>
flist = GetFileList();
2509 for (
size_t i = 0;
i < flist.size(); ++
i)
2510 flistout <<
"\n [" << std::setw(3) <<
i <<
"] " << flist[
i];
2512 TLorentzVector usr0(0,0,0,0);
2513 TLorentzVector usr0asbeam;
2514 User2BeamPos(usr0,usr0asbeam);
2516 const int w=10,
p=6;
2517 std::ostringstream beamrot_str, beamrotinv_str;
2519 <<
"fBeamRot: " << std::setprecision(
p) <<
"\n" 2521 << std::setw(w) << fBeamRot.XX() <<
" " 2522 << std::setw(w) << fBeamRot.XY() <<
" " 2523 << std::setw(w) << fBeamRot.XZ() <<
" ]\n" 2525 << std::setw(w) << fBeamRot.YX() <<
" " 2526 << std::setw(w) << fBeamRot.YY() <<
" " 2527 << std::setw(w) << fBeamRot.YZ() <<
" ]\n" 2529 << std::setw(w) << fBeamRot.ZX() <<
" " 2530 << std::setw(w) << fBeamRot.ZY() <<
" " 2531 << std::setw(w) << fBeamRot.ZZ() <<
" ]";
2533 <<
"fBeamRotInv: " << std::setprecision(
p) <<
"\n" 2535 << std::setw(w) << fBeamRotInv.XX() <<
" " 2536 << std::setw(w) << fBeamRotInv.XY() <<
" " 2537 << std::setw(w) << fBeamRotInv.XZ() <<
" ]\n" 2539 << std::setw(w) << fBeamRotInv.YX() <<
" " 2540 << std::setw(w) << fBeamRotInv.YY() <<
" " 2541 << std::setw(w) << fBeamRotInv.YZ() <<
" ]\n" 2543 << std::setw(w) << fBeamRotInv.ZX() <<
" " 2544 << std::setw(w) << fBeamRotInv.ZY() <<
" " 2545 << std::setw(w) << fBeamRotInv.ZZ() <<
" ]";
2548 <<
"GNuMIFlux Config:" 2549 <<
"\n Enu_max " << fMaxEv
2550 <<
"\n pdg-codes: " << s.str() <<
"\n " 2551 << (fG3NuMI?
"g3numi":
"")
2552 << (fG4NuMI?
"g4numi":
"")
2553 << (fFlugg?
"flugg":
"")
2554 <<
"/" << fNuFluxGen <<
" " 2555 <<
"(" << fNuFluxTreeName <<
"), " << fNEntries <<
" entries" 2556 <<
" (FilePOTs " << fFilePOTs <<
") " 2557 <<
"in " << fNFiles <<
" files: " 2559 <<
"\n from file patterns:" 2561 <<
"\n wgt max=" << fMaxWeight <<
" fudge=" << fMaxWgtFudge <<
" using " 2562 << fMaxWgtEntries <<
" entries" 2563 <<
"\n Z0 pushback " << fZ0
2564 <<
"\n used entry " << fIEntry <<
" " << fIUse <<
"/" << fNUse
2565 <<
" times, in " << fICycle <<
"/" << fNCycles <<
" cycles" 2566 <<
"\n SumWeight " << fSumWeight <<
" for " << fNNeutrinos <<
" neutrinos" 2567 <<
"\n EffPOTsPerNu " << fEffPOTsPerNu <<
" AccumPOTs " << fAccumPOTs
2568 <<
"\n GenWeighted \"" << (fGenWeighted?
"true":
"false") <<
", " 2569 <<
"\", Detector location set \"" << (fDetLocIsSet?
"true":
"false") <<
"\"" 2570 <<
"\n LengthUnits " << fLengthUnits <<
", scale b2u " << fLengthScaleB2U
2571 <<
", scale u2b " << fLengthScaleU2B
2572 <<
"\n User Flux Window: " 2576 <<
"\n Flux Window (cm, beam coord): " 2581 <<
"\n User Beam Origin: " 2583 <<
"\n " << beamrot_str.str() <<
" " 2584 <<
"\n Detector Origin (beam coord): " 2586 <<
"\n " << beamrotinv_str.str() <<
" " 2587 <<
"\n UseFluxAtDetCenter " << fUseFluxAtDetCenter;
2594 std::vector<std::string>
flist;
2595 TObjArray *fileElements=fNuFluxTree->GetListOfFiles();
2596 TIter
next(fileElements);
2597 TChainElement *chEl=0;
2598 while (( chEl=(TChainElement*)
next() )) {
2599 flist.push_back(chEl->GetTitle());
2611 std::vector<double> vect;
2612 size_t ntok = strtokens.size();
2617 for (
size_t i=0;
i < ntok; ++
i) {
2619 if (
" " == trimmed ||
"" == trimmed )
continue;
2620 double val = strtod(trimmed.c_str(), (
char**)NULL);
2623 vect.push_back(val);
2634 std::vector<long int> vect;
2635 size_t ntok = strtokens.size();
2640 for (
size_t i=0;
i < ntok; ++
i) {
2642 if (
" " == trimmed ||
"" == trimmed )
continue;
2643 long int val = strtol(trimmed.c_str(),(
char**)NULL,10);
2646 vect.push_back(val);
2656 bool is_accessible = ! (gSystem->AccessPathName(fname.c_str()));
2657 if (!is_accessible) {
2659 <<
"The XML doc doesn't exist! (filename: " << fname <<
")";
2663 xmlDocPtr xml_doc = xmlParseFile( fname.c_str() );
2664 if ( xml_doc == NULL) {
2666 <<
"The XML doc can't be parsed! (filename: " << fname <<
")";
2670 xmlNodePtr xml_root = xmlDocGetRootElement( xml_doc );
2671 if ( xml_root == NULL ) {
2673 <<
"The XML doc is empty! (filename: " << fname <<
")";
2677 string rootele =
"gnumi_config";
2678 if ( xmlStrcmp(xml_root->name, (
const xmlChar*)rootele.c_str() ) ) {
2680 <<
"The XML doc has invalid root element! (filename: " << fname <<
")" 2681 <<
" expected \"" << rootele <<
"\", saw \"" << xml_root->name <<
"\"";
2685 SLOG(
"GNuMIFlux",
pINFO) <<
"Attempt to load config \"" << cfg
2686 <<
"\" from file: " <<
fname;
2688 bool found = this->LoadParamSet(xml_doc,cfg);
2698 xmlNodePtr xml_root = xmlDocGetRootElement( xml_doc );
2705 xmlNodePtr xml_pset = xml_root->xmlChildrenNode;
2706 for ( ; xml_pset != NULL ; xml_pset = xml_pset->next ) {
2707 if ( ! xmlStrEqual(xml_pset->name, (
const xmlChar*)
"param_set") )
continue;
2709 string param_set_name =
2712 if ( param_set_name != cfg )
continue;
2716 this->ParseParamSet(xml_doc,xml_pset);
2727 xmlNodePtr xml_child = xml_pset->xmlChildrenNode;
2728 for ( ; xml_child != NULL ; xml_child = xml_child->next ) {
2733 if ( pname ==
"text" || pname ==
"comment" )
continue;
2736 xmlNodeListGetString(xml_doc, xml_child->xmlChildrenNode, 1));
2740 <<
" pname \"" << pname <<
"\", string value \"" << pval <<
"\"";
2742 if ( pname ==
"verbose" ) {
2743 fVerbose = atoi(pval.c_str());
2745 }
else if ( pname ==
"using_param_set" ) {
2746 SLOG(
"GNuMIFlux",
pWARN) <<
"start using_param_set: \"" << pval <<
"\"";
2747 bool found = this->LoadParamSet(xml_doc,pval);
2749 SLOG(
"GNuMIFlux",
pFATAL) <<
"using_param_set: \"" << pval <<
"\" NOT FOUND";
2752 SLOG(
"GNuMIFlux",
pWARN) <<
"done using_param_set: \"" << pval <<
"\"";
2753 }
else if ( pname ==
"units" ) {
2755 fGNuMI->SetLengthUnits(scale);
2756 SLOG(
"GNuMIFlux",
pINFO) <<
"set user units to \"" << pval <<
"\"";
2758 }
else if ( pname ==
"beamdir" ) {
2759 ParseBeamDir(xml_doc,xml_child);
2760 fGNuMI->SetBeamRotation(fBeamRotXML);
2762 }
else if ( pname ==
"beampos" ) {
2764 fGNuMI->SetBeamCenter(fBeamPosXML);
2766 }
else if ( pname ==
"window" ) {
2767 ParseWindowSeries(xml_doc,xml_child);
2774 fGNuMI->SetFluxWindow(fFluxWindowPtXML[0],
2775 fFluxWindowPtXML[1],
2776 fFluxWindowPtXML[2]);
2778 }
else if ( pname ==
"enumax" ) {
2781 }
else if ( pname ==
"upstreamz" ) {
2782 double z0usr = -3.4e38;
2783 std::vector<double>
v = GetDoubleVector(pval);
2784 if ( v.size() > 0 ) z0usr = v[0];
2785 fGNuMI->SetUpstreamZ(z0usr);
2786 SLOG(
"GNuMIFlux",
pINFO) <<
"set upstreamz = " << z0usr;
2788 }
else if ( pname ==
"reuse" ) {
2789 long int nreuse = 1;
2790 std::vector<long int>
v = GetIntVector(pval);
2791 if ( v.size() > 0 ) nreuse = v[0];
2792 fGNuMI->SetEntryReuse(nreuse);
2793 SLOG(
"GNuMIFlux",
pINFO) <<
"set entry reuse = " << nreuse;
2797 <<
" NOT HANDLED: pname \"" << pname
2798 <<
"\", string value \"" << pval <<
"\"";
2808 fBeamRotXML.SetToIdentity();
2816 xmlNodeListGetString(xml_doc, xml_beamdir->xmlChildrenNode, 1));
2818 if ( dirtype ==
"series" ) {
2820 ParseRotSeries(xml_doc,xml_beamdir);
2822 }
else if ( dirtype ==
"thetaphi3") {
2824 std::vector<double> thetaphi3 = GetDoubleVector(pval);
2827 if ( thetaphi3.size() == 6 ) {
2829 TVector3 newX = AnglesToAxis(thetaphi3[0],thetaphi3[1],units);
2830 TVector3 newY = AnglesToAxis(thetaphi3[2],thetaphi3[3],units);
2831 TVector3 newZ = AnglesToAxis(thetaphi3[4],thetaphi3[5],units);
2832 fTempRot.RotateAxes(newX,newY,newZ);
2833 fBeamRotXML = fTempRot;
2836 <<
" type=\"" << dirtype <<
"\" within <beamdir> needs 6 values";
2839 }
else if ( dirtype ==
"newxyz" ) {
2841 std::vector<double> newdir = GetDoubleVector(pval);
2842 if ( newdir.size() == 9 ) {
2844 TVector3 newX = TVector3(newdir[0],newdir[1],newdir[2]).Unit();
2845 TVector3 newY = TVector3(newdir[3],newdir[4],newdir[5]).Unit();
2846 TVector3 newZ = TVector3(newdir[6],newdir[7],newdir[8]).Unit();
2847 fTempRot.RotateAxes(newX,newY,newZ);
2848 fBeamRotXML = fTempRot.Inverse();
2851 <<
" type=\"" << dirtype <<
"\" within <beamdir> needs 9 values";
2857 <<
" UNHANDLED type=\"" << dirtype <<
"\" within <beamdir>";
2860 if ( fVerbose > 1 ) {
2864 << std::setw(w) << fBeamRotXML.XX() <<
" " 2865 << std::setw(w) << fBeamRotXML.XY() <<
" " 2866 << std::setw(w) << fBeamRotXML.XZ() <<
endl 2868 << std::setw(w) << fBeamRotXML.YX() <<
" " 2869 << std::setw(w) << fBeamRotXML.YY() <<
" " 2870 << std::setw(w) << fBeamRotXML.YZ() <<
endl 2872 << std::setw(w) << fBeamRotXML.ZX() <<
" " 2873 << std::setw(w) << fBeamRotXML.ZY() <<
" " 2874 << std::setw(w) << fBeamRotXML.ZZ() <<
" ] " <<
std::endl;
2882 std::vector<double> xyz = GetDoubleVector(str);
2883 if ( xyz.size() == 3 ) {
2884 fBeamPosXML = TVector3(xyz[0],xyz[1],xyz[2]);
2885 }
else if ( xyz.size() == 6 ) {
2888 TVector3 userpos(xyz[0],xyz[1],xyz[2]);
2889 TVector3 beampos(xyz[3],xyz[4],xyz[5]);
2890 fBeamPosXML = userpos - fBeamRotXML*beampos;
2893 <<
"Unable to parse " << xyz.size() <<
" values in <beampos>";
2896 if ( fVerbose > 1 ) {
2898 std::cout <<
" fBeamPosXML: [ " << std::setprecision(
p)
2899 << std::setw(w) << fBeamPosXML.X() <<
" , " 2900 << std::setw(w) << fBeamPosXML.Y() <<
" , " 2901 << std::setw(w) << fBeamPosXML.Z() <<
" ] " 2908 std::vector<double>
v = GetDoubleVector(str);
2909 size_t n = v.size();
2911 fGNuMI->SetMaxEnergy(v[0]);
2916 fGNuMI->SetMaxEFudge(v[1]);
2922 fGNuMI->SetMaxWgtScan(v[2]);
2927 fGNuMI->SetMaxWgtScan(v[2],nentries);
2929 std::cout <<
"ParseEnuMax SetMaxWgtScan(" << v[2] <<
"," << nentries <<
")" <<
std::endl;
2938 xmlNodePtr xml_child = xml_pset->xmlChildrenNode;
2939 for ( ; xml_child != NULL ; xml_child = xml_child->next ) {
2944 if ( name ==
"text" || name ==
"comment" )
continue;
2946 if ( name ==
"rotation" ) {
2948 xmlNodeListGetString(xml_doc, xml_child->xmlChildrenNode, 1));
2955 double rot = atof(val.c_str());
2957 if (
'd' == units[0] ||
'D' == units[0] ) rot *= TMath::DegToRad();
2961 <<
" rotate " << rot <<
" radians around " << axis <<
" axis";
2963 if ( axis[0] ==
'x' || axis[0] ==
'X' ) fTempRot.RotateX(rot);
2964 else if ( axis[0] ==
'y' || axis[0] ==
'Y' ) fTempRot.RotateY(rot);
2965 else if ( axis[0] ==
'z' || axis[0] ==
'Z' ) fTempRot.RotateZ(rot);
2968 <<
" no " << axis <<
" to rotate around";
2973 <<
" found <" << name <<
"> within <beamdir type=\"series\">";
2977 fBeamRotXML = fTempRot.Inverse();
2985 xmlNodePtr xml_child = xml_pset->xmlChildrenNode;
2986 for ( ; xml_child != NULL ; xml_child = xml_child->next ) {
2991 if ( name ==
"text" || name ==
"comment" )
continue;
2993 if ( name ==
"point" ) {
2996 xmlNodeListGetString(xml_doc, xml_child->xmlChildrenNode, 1));
3000 std::vector<double> xyz = GetDoubleVector(val);
3001 if ( xyz.size() != 3 || coord !=
"det" ) {
3003 <<
"parsing <window> found <point> but size=" << xyz.size()
3004 <<
" (expect 3) and coord=\"" << coord <<
"\" (expect \"det\")" 3005 <<
" IGNORE problem";
3008 if ( ientry < 3 && ientry >= 0 ) {
3009 TVector3
pt(xyz[0],xyz[1],xyz[2]);
3010 if ( fVerbose > 0 ) {
3012 std::cout <<
" point[" << ientry <<
"] = [ " << std::setprecision(
p)
3013 << std::setw(w) << pt.X() <<
" , " 3014 << std::setw(w) << pt.Y() <<
" , " 3015 << std::setw(w) << pt.Z() <<
" ] " 3018 fFluxWindowPtXML[ientry] =
pt;
3021 <<
" <window><point> ientry " << ientry <<
" out of range (0-2)";
3026 <<
" found <" << name <<
"> within <window>";
3036 double scale = (
'd'==units[0]||
'D'==units[0]) ? TMath::DegToRad() : 1.0 ;
3038 xyz[0] = TMath::Cos(scale*phi)*TMath::Sin(scale*theta);
3039 xyz[1] = TMath::Sin(scale*phi)*TMath::Sin(scale*theta);
3040 xyz[2] = TMath::Cos(scale*theta);
3042 for (
int i=0;
i<3; ++
i) {
3043 const double eps = 1.0e-15;
3044 if (TMath::Abs(xyz[
i]) < eps ) xyz[
i] = 0;
3045 if (TMath::Abs(xyz[i]-1) < eps ) xyz[
i] = 1;
3046 if (TMath::Abs(xyz[i]+1) < eps ) xyz[
i] = -1;
3048 return TVector3(xyz[0],xyz[1],xyz[2]);
3053 std::vector<double> xyz = GetDoubleVector(str);
3054 if ( xyz.size() != 3 ) {
3057 <<
" ParseTV3 \"" << str <<
"\" had " << xyz.size() <<
" elements ";
3059 return TVector3(xyz[0],xyz[1],xyz[2]);
3064 #ifndef SKIP_MINERVA_MODS 3070 if(sval==
"AddedLV")ival=1;
3071 else if(sval==
"AlTube1LV")ival=2;
3072 else if(sval==
"AlTube2LV")ival=3;
3073 else if(sval==
"Al_BLK1")ival=4;
3074 else if(sval==
"Al_BLK2")ival=5;
3075 else if(sval==
"Al_BLK3")ival=6;
3076 else if(sval==
"Al_BLK4")ival=7;
3077 else if(sval==
"Al_BLK5")ival=8;
3078 else if(sval==
"Al_BLK6")ival=9;
3079 else if(sval==
"Al_BLK7")ival=10;
3080 else if(sval==
"Al_BLK8")ival=11;
3081 else if(sval==
"AlholeL")ival=12;
3082 else if(sval==
"AlholeR")ival=13;
3083 else if(sval==
"BEndLV")ival=14;
3084 else if(sval==
"BFrontLV")ival=15;
3085 else if(sval==
"BeDWLV")ival=16;
3086 else if(sval==
"BeUp1LV")ival=17;
3087 else if(sval==
"BeUp2LV")ival=18;
3088 else if(sval==
"BeUp3LV")ival=19;
3089 else if(sval==
"BodyLV")ival=20;
3090 else if(sval==
"BudalMonitor")ival=21;
3091 else if(sval==
"CLid1LV")ival=22;
3092 else if(sval==
"CLid2LV")ival=23;
3093 else if(sval==
"CShld_BLK11")ival=24;
3094 else if(sval==
"CShld_BLK12")ival=25;
3095 else if(sval==
"CShld_BLK2")ival=26;
3096 else if(sval==
"CShld_BLK3")ival=27;
3097 else if(sval==
"CShld_BLK4")ival=28;
3098 else if(sval==
"CShld_BLK7")ival=29;
3099 else if(sval==
"CShld_BLK8")ival=30;
3100 else if(sval==
"CShld_stl,BLK")ival=31;
3101 else if(sval==
"CerTubeLV")ival=32;
3102 else if(sval==
"CeramicRod")ival=33;
3103 else if(sval==
"ConcShield")ival=34;
3104 else if(sval==
"Concrete Chase Section")ival=35;
3105 else if(sval==
"Conn1LV")ival=36;
3106 else if(sval==
"Conn2LV")ival=37;
3107 else if(sval==
"Conn3LV")ival=38;
3108 else if(sval==
"DNWN")ival=39;
3109 else if(sval==
"DPIP")ival=40;
3110 else if(sval==
"DVOL")ival=41;
3111 else if(sval==
"DuratekBlock")ival=42;
3112 else if(sval==
"DuratekBlockCovering")ival=43;
3113 else if(sval==
"HadCell")ival=44;
3114 else if(sval==
"HadronAbsorber")ival=45;
3115 else if(sval==
"MuMonAlcvFill_0")ival=46;
3116 else if(sval==
"MuMonAlcv_0")ival=47;
3117 else if(sval==
"MuMonAlcv_1")ival=48;
3118 else if(sval==
"MuMonAlcv_2")ival=49;
3119 else if(sval==
"MuMon_0")ival=50;
3120 else if(sval==
"MuMon_1")ival=51;
3121 else if(sval==
"MuMon_2")ival=52;
3122 else if(sval==
"PHorn1CPB1slv")ival=53;
3123 else if(sval==
"PHorn1CPB2slv")ival=54;
3124 else if(sval==
"PHorn1F")ival=55;
3125 else if(sval==
"PHorn1Front")ival=56;
3126 else if(sval==
"PHorn1IC")ival=57;
3127 else if(sval==
"PHorn1InsRingslv")ival=58;
3128 else if(sval==
"PHorn1OC")ival=59;
3129 else if(sval==
"PHorn2CPB1slv")ival=60;
3130 else if(sval==
"PHorn2CPB2slv")ival=61;
3131 else if(sval==
"PHorn2F")ival=62;
3132 else if(sval==
"PHorn2Front")ival=63;
3133 else if(sval==
"PHorn2IC")ival=64;
3134 else if(sval==
"PHorn2InsRingslv")ival=65;
3135 else if(sval==
"PHorn2OC")ival=66;
3136 else if(sval==
"PVHadMon")ival=67;
3137 else if(sval==
"Pipe1")ival=68;
3138 else if(sval==
"Pipe1_water")ival=69;
3139 else if(sval==
"Pipe2")ival=70;
3140 else if(sval==
"Pipe2_water")ival=71;
3141 else if(sval==
"Pipe3")ival=72;
3142 else if(sval==
"Pipe3_water")ival=73;
3143 else if(sval==
"Pipe4")ival=74;
3144 else if(sval==
"Pipe4_water")ival=75;
3145 else if(sval==
"Pipe5")ival=76;
3146 else if(sval==
"Pipe5_water")ival=77;
3147 else if(sval==
"Pipe6")ival=78;
3148 else if(sval==
"Pipe6_water")ival=79;
3149 else if(sval==
"Pipe7")ival=80;
3150 else if(sval==
"Pipe8")ival=81;
3151 else if(sval==
"Pipe8_water")ival=82;
3152 else if(sval==
"Pipe9")ival=83;
3153 else if(sval==
"PipeAdapter1")ival=84;
3154 else if(sval==
"PipeAdapter1_water")ival=85;
3155 else if(sval==
"PipeAdapter2")ival=86;
3156 else if(sval==
"PipeAdapter2_water")ival=87;
3157 else if(sval==
"PipeBellowB")ival=88;
3158 else if(sval==
"PipeBellowB_water")ival=89;
3159 else if(sval==
"PipeBellowT")ival=90;
3160 else if(sval==
"PipeBellowT_water")ival=91;
3161 else if(sval==
"PipeC1")ival=92;
3162 else if(sval==
"PipeC1_water")ival=93;
3163 else if(sval==
"PipeC2")ival=94;
3164 else if(sval==
"PipeC2_water")ival=95;
3165 else if(sval==
"ROCK")ival=96;
3166 else if(sval==
"Ring1LV")ival=97;
3167 else if(sval==
"Ring2LV")ival=98;
3168 else if(sval==
"Ring3LV")ival=99;
3169 else if(sval==
"Ring4LV")ival=100;
3170 else if(sval==
"Ring5LV")ival=101;
3171 else if(sval==
"SC01")ival=102;
3172 else if(sval==
"SpiderSupport")ival=103;
3173 else if(sval==
"Stl_BLK1")ival=104;
3174 else if(sval==
"Stl_BLK10")ival=105;
3175 else if(sval==
"Stl_BLK2")ival=106;
3176 else if(sval==
"Stl_BLK3")ival=107;
3177 else if(sval==
"Stl_BLK4")ival=108;
3178 else if(sval==
"Stl_BLK5")ival=109;
3179 else if(sval==
"Stl_BLK6")ival=110;
3180 else if(sval==
"Stl_BLK7")ival=111;
3181 else if(sval==
"Stl_BLK8")ival=112;
3182 else if(sval==
"Stl_BLK9")ival=113;
3183 else if(sval==
"Stlhole")ival=114;
3184 else if(sval==
"TGAR")ival=115;
3185 else if(sval==
"TGT1")ival=116;
3186 else if(sval==
"TGTExitCyl2LV")ival=117;
3187 else if(sval==
"TUNE")ival=118;
3188 else if(sval==
"Tube1aLV")ival=119;
3189 else if(sval==
"Tube1bLV")ival=120;
3190 else if(sval==
"UpWn1")ival=121;
3191 else if(sval==
"UpWn2")ival=122;
3192 else if(sval==
"UpWnAl1SLV")ival=123;
3193 else if(sval==
"UpWnAl2SLV")ival=124;
3194 else if(sval==
"UpWnAl3SLV")ival=125;
3195 else if(sval==
"UpWnFe1SLV")ival=126;
3196 else if(sval==
"UpWnFe2SLV")ival=127;
3197 else if(sval==
"UpWnPolyCone")ival=128;
3198 else if(sval==
"blu_BLK25")ival=129;
3199 else if(sval==
"blu_BLK26")ival=130;
3200 else if(sval==
"blu_BLK27")ival=131;
3201 else if(sval==
"blu_BLK28")ival=132;
3202 else if(sval==
"blu_BLK29")ival=133;
3203 else if(sval==
"blu_BLK32")ival=134;
3204 else if(sval==
"blu_BLK37")ival=135;
3205 else if(sval==
"blu_BLK38")ival=136;
3206 else if(sval==
"blu_BLK39")ival=137;
3207 else if(sval==
"blu_BLK40")ival=138;
3208 else if(sval==
"blu_BLK45")ival=139;
3209 else if(sval==
"blu_BLK46")ival=140;
3210 else if(sval==
"blu_BLK47")ival=141;
3211 else if(sval==
"blu_BLK48")ival=142;
3212 else if(sval==
"blu_BLK49")ival=143;
3213 else if(sval==
"blu_BLK50")ival=144;
3214 else if(sval==
"blu_BLK51")ival=145;
3215 else if(sval==
"blu_BLK53")ival=146;
3216 else if(sval==
"blu_BLK55")ival=147;
3217 else if(sval==
"blu_BLK57")ival=148;
3218 else if(sval==
"blu_BLK59")ival=149;
3219 else if(sval==
"blu_BLK61")ival=150;
3220 else if(sval==
"blu_BLK63")ival=151;
3221 else if(sval==
"blu_BLK64")ival=152;
3222 else if(sval==
"blu_BLK65")ival=153;
3223 else if(sval==
"blu_BLK66")ival=154;
3224 else if(sval==
"blu_BLK67")ival=155;
3225 else if(sval==
"blu_BLK68")ival=156;
3226 else if(sval==
"blu_BLK69")ival=157;
3227 else if(sval==
"blu_BLK70")ival=158;
3228 else if(sval==
"blu_BLK72")ival=159;
3229 else if(sval==
"blu_BLK73")ival=160;
3230 else if(sval==
"blu_BLK75")ival=161;
3231 else if(sval==
"blu_BLK77")ival=162;
3232 else if(sval==
"blu_BLK78")ival=163;
3233 else if(sval==
"blu_BLK79")ival=164;
3234 else if(sval==
"blu_BLK81")ival=165;
3235 else if(sval==
"conc_BLK")ival=166;
3236 else if(sval==
"pvBaffleMother")ival=167;
3237 else if(sval==
"pvDPInnerTrackerTube")ival=168;
3238 else if(sval==
"pvMHorn1Mother")ival=169;
3239 else if(sval==
"pvMHorn2Mother")ival=170;
3240 else if(sval==
"pvTargetMother")ival=171;
3241 else if(sval==
"stl_slab1")ival=172;
3242 else if(sval==
"stl_slab4")ival=173;
3243 else if(sval==
"stl_slab5")ival=174;
3244 else if(sval==
"stl_slabL")ival=175;
3245 else if(sval==
"stl_slabR")ival=176;
3251 if(sval==
"AntiLambdaInelastic")ival=1;
3252 else if(sval==
"AntiNeutronInelastic")ival=2;
3253 else if(sval==
"AntiOmegaMinusInelastic")ival=3;
3254 else if(sval==
"AntiProtonInelastic")ival=4;
3255 else if(sval==
"AntiSigmaMinusInelastic")ival=5;
3256 else if(sval==
"AntiSigmaPlusInelastic")ival=6;
3257 else if(sval==
"AntiXiMinusInelastic")ival=7;
3258 else if(sval==
"AntiXiZeroInelastic")ival=8;
3259 else if(sval==
"Decay")ival=9;
3260 else if(sval==
"KaonMinusInelastic")ival=10;
3261 else if(sval==
"KaonPlusInelastic")ival=11;
3262 else if(sval==
"KaonZeroLInelastic")ival=12;
3263 else if(sval==
"KaonZeroSInelastic")ival=13;
3264 else if(sval==
"LambdaInelastic")ival=14;
3265 else if(sval==
"NeutronInelastic")ival=15;
3266 else if(sval==
"OmegaMinusInelastic")ival=16;
3267 else if(sval==
"PionMinusInelastic")ival=17;
3268 else if(sval==
"PionPlusInelastic")ival=18;
3269 else if(sval==
"Primary")ival=19;
3270 else if(sval==
"ProtonInelastic")ival=20;
3271 else if(sval==
"SigmaMinusInelastic")ival=21;
3272 else if(sval==
"SigmaPlusInelastic")ival=22;
3273 else if(sval==
"XiMinusInelastic")ival=23;
3274 else if(sval==
"XiZeroInelastic")ival=24;
3275 else if(sval==
"hElastic")ival=25;
bool LoadParamSet(xmlDocPtr &, std::string cfg)
TVector3 AnglesToAxis(double theta, double phi, std::string units="deg")
const XML_Char XML_Encoding * info
A class for generating concrete GFluxI derived classes based on the factory pattern. This code supplies a CPP macro which allows the classes to self-register and thus no modification of this class is needed in order to expand the list of classes it knows about.
std::vector< long int > GetIntVector(std::string str)
void User2BeamP4(const TLorentzVector &usrp4, TLorentzVector &beamp4) const
TVector3 GetBeamCenter() const
beam origin in user frame
virtual void GetBranchInfo(std::vector< std::string > &branchNames, std::vector< std::string > &branchClassNames, std::vector< void ** > &branchObjPointers)
TLorentzVector fgP4
generated nu 4-momentum beam coord
string P4AsShortString(const TLorentzVector *p)
void GetFluxWindow(TVector3 &p1, TVector3 &p2, TVector3 &p3) const
3 points define a plane in beam coordinate
void SetEntryReuse(long int nuse=1)
of times to use entry before moving to next
static const unsigned int MAX_N_TRAJ
Maximum number of trajectories to store.
THE MAIN GENIE PROJECT NAMESPACE
double POT_curr(void)
current average POT (RWH?)
Int_t run
current Tree number in a TChain
static RandomGen * Instance()
Access instance.
void SetLengthUnits(double user_units)
Set units assumed by user.
void UseFluxAtNearDetCenter(void)
force weights at MINOS detector "center" as found in ntuple
double pprodpx[MAX_N_TRAJ]
def AddFile(fout, deets, f)
string TrimSpaces(xmlChar *xmls)
TLorentzVector fgX4User
generated nu 4-position user coord
bool LoadConfig(string cfg)
load a named configuration
void ParseParamSet(xmlDocPtr &, xmlNodePtr &)
double stoppz[MAX_N_TRAJ]
double stoppy[MAX_N_TRAJ]
void ParseBeamDir(xmlDocPtr &, xmlNodePtr &)
int fVerbose
how noisy to be when parsing XML
Int_t run
current Tree number in a TChain
void SetTreeName(string name)
set input tree name (default: "h10")
virtual double GetTotalExposure() const
GNuMIFluxXMLHelper(GNuMIFlux *gnumi)
void Beam2UserP4(const TLorentzVector &beamp4, TLorentzVector &usrp4) const
void Compare(const std::vector< TH1 * > &hs)
TRotation GetBeamRotation() const
rotation to apply from beam->user
double startx[MAX_N_TRAJ]
double pprodpy[MAX_N_TRAJ]
A singleton holding random number generator classes. All random number generation in GENIE should tak...
void Initialize(Bool_t useTMVAStyle=kTRUE)
void MakeCopy(const g3numi *)
pull in from g3 ntuple
friend ostream & operator<<(ostream &stream, const GNuMIFluxPassThroughInfo &info)
A GENIE flux driver encapsulating the NuMI neutrino flux. It reads-in the official GNUMI neutrino flu...
void User2BeamPos(const TLorentzVector &usrxyz, TLorentzVector &beamxyz) const
int GeantToPdg(int geant_code)
TLorentzVector fgX4
generated nu 4-position beam coord
virtual void LoadBeamSimData(const std::vector< std::string > &filenames, const std::string &det_loc)
void MoveToZ0(double z0)
move ray origin to user coord Z0
void ParseBeamPos(std::string)
void ParseEnuMax(std::string)
std::vector< std::string > GetFileList()
list of files currently part of chain
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Int_t run
current Tree number in a TChain
static constexpr double L
void SetBeamRotation(TRotation beamrot)
< beam (0,0,0) relative to user frame, beam direction in user frame
double UnitFromString(string u)
GENIE interface for uniform flux exposure iterface.
string GetXMLFilePath(string basename)
void ParseRotSeries(xmlDocPtr &, xmlNodePtr &)
double pprodpz[MAX_N_TRAJ]
int CalcEnuWgt(const TLorentzVector &xyz, double &enu, double &wgt_xy) const
void SetBeamCenter(TVector3 beam0)
bool GenerateNext_weighted(void)
int getVolID(TString sval)
GNuMIFluxPassThroughInfo()
void PrintCurrent(void)
print current entry from leaves
bool SetFluxWindow(StdFluxWindow_t stdwindow, double padding=0)
return false if unhandled
bool GenerateNext(void)
generate the next flux neutrino (return false in err)
static constexpr Double_t rad
bool LoadConfig(std::string cfg)
double startpz[MAX_N_TRAJ]
double LengthUnits(void) const
Return user units.
double startz[MAX_N_TRAJ]
void Print(const Option_t *opt="") const
ClassImp(GNuMIFluxPassThroughInfo) const TLorentzVector kPosCenterNearBeam(0.
string TrimSpaces(string input)
void Clear(Option_t *opt)
reset state variables based on opt
void Print(std::string prefix, std::string name, std::string suffix="")
TVector3 ParseTV3(const std::string &)
double GetDecayDist() const
dist (user units) from dk to current pos
void Beam2UserPos(const TLorentzVector &beamxyz, TLorentzVector &usrxyz) const
double dot(const std::vector< double > &x, const std::vector< double > &y)
const TLorentzVector kPosCenterFarBeam(0., 0., 735340.00, 0.)
double startpx[MAX_N_TRAJ]
double startpy[MAX_N_TRAJ]
vector< string > Split(string input, string delim)
Var Sqrt(const Var &v)
Use to take sqrt of a var.
int fgPdgC
generated nu pdg-code
enum genie::flux::GNuMIFlux::EStdFluxWindow StdFluxWindow_t
void ParseWindowSeries(xmlDocPtr &, xmlNodePtr &)
virtual TTree * GetMetaDataTree()
double stoppx[MAX_N_TRAJ]
assert(nhit_max >=nhit_nbins)
void User2BeamDir(const TLorentzVector &usrdir, TLorentzVector &beamdir) const
void PrintConfig()
print the current configuration
void ScanForMaxWeight(void)
scan for max flux weight (before generating unweighted flux neutrinos)
int getProcessID(TString sval)
void SetMaxEnergy(double Ev)
specify maximum flx neutrino energy
TLorentzVector fgP4User
generated nu 4-momentum user coord
std::vector< double > GetDoubleVector(std::string str)
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
string X4AsString(const TLorentzVector *x)
double starty[MAX_N_TRAJ]
virtual long int NFluxNeutrinos() const
of rays generated
void Beam2UserDir(const TLorentzVector &beamdir, TLorentzVector &usrdir) const
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
string GetAttribute(xmlNodePtr xml_cur, string attr_name)
static constexpr Double_t cm
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
#define FLUXDRIVERREG4(_nsa, _nsb, _name, _fqname)
string Vec3AsString(const TVector3 *vec)
void AddFile(TTree *tree, string fname)
void UseFluxAtFarDetCenter(void)
void GenerateWeighted(bool gen_weighted)
set whether to generate weighted or unweighted neutrinos
void CalcEffPOTsPerNu(void)
void push_back(int pdg_code)
double UsedPOTs(void) const
of protons-on-target used