All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Attributes | List of all members
fluxr::FluxReader Class Reference

#include <FluxReader.h>

Public Member Functions

 FluxReader (fhicl::ParameterSet const &pset, art::ProductRegistryHelper &helper, art::SourceHelper const &pm)
 
void closeCurrentFile ()
 
void readFile (std::string const &name, art::FileBlock *&fb)
 
bool readNext (art::RunPrincipal *const &inR, art::SubRunPrincipal *const &inSR, art::RunPrincipal *&outR, art::SubRunPrincipal *&outSR, art::EventPrincipal *&outE)
 

Private Attributes

art::SourceHelper const & fSourceHelper
 
art::SubRunID fSubRunID
 
uint32_t fEventCounter
 
uint32_t fEntry
 
int fMaxEvents
 
uint32_t fSkipEvents
 
std::string fInputType
 
float fPOT
 
float fCurrentPOT
 
bool fSelfIncrementRuns
 
art::RunNumber_t fIncrement = 1
 
FluxInterfacefFluxDriver
 
TFile * fFluxInputFile
 
TH1D * fHFlux [4]
 
TH1D * fHFluxParent [4][4]
 
TH1D * fHFluxSec [4][5]
 
fhicl::ParameterSet fConfigPS
 
art::TypeLabel fTLmctruth
 
art::TypeLabel fTLmcflux
 
art::TypeLabel fTLdk2nu
 
art::TypeLabel fTLnuchoice
 

Detailed Description

Definition at line 20 of file FluxReader.h.

Constructor & Destructor Documentation

fluxr::FluxReader::FluxReader ( fhicl::ParameterSet const &  pset,
art::ProductRegistryHelper &  helper,
art::SourceHelper const &  pm 
)

Definition at line 47 of file FluxReader.cxx.

50  :
51  fSourceHelper(pm),
52  fSubRunID(),
53  fInputType(ps.get<std::string>("inputType")),
54  fSelfIncrementRuns(ps.get<bool>("SelfIncrementRun")),
55  fIncrement(1),
56  fTLmctruth{helper.reconstitutes<std::vector<simb::MCTruth>, art::InEvent>("flux")},
57  fTLmcflux{helper.reconstitutes<std::vector<simb::MCFlux>, art::InEvent>("flux")},
58  fTLdk2nu{fInputType=="dk2nu" ? helper.reconstitutes<std::vector<bsim::Dk2Nu>, art::InEvent>("flux"):fTLmctruth},
59  fTLnuchoice{fInputType=="dk2nu" ? helper.reconstitutes<std::vector<bsim::NuChoice>, art::InEvent>("flux"):fTLmctruth}
60  {
61 
62  fPOT = fCurrentPOT = 0;
63 
64  helper.reconstitutes<sumdata::POTSummary, art::InSubRun >("flux");
65 
66  if (fInputType=="dk2nu") {
67  helper.reconstitutes<art::Assns<simb::MCTruth, bsim::Dk2Nu>, art::InEvent>("flux");
68  helper.reconstitutes<art::Assns<simb::MCTruth, bsim::NuChoice>, art::InEvent>("flux");
69  helper.reconstitutes<art::Assns<simb::MCTruth, simb::MCFlux>, art::InEvent>("flux");
70 
71  fConfigPS=ps.get<fhicl::ParameterSet>(ps.get<std::string>("dk2nuConfig"));
72  }
73 
74  art::ServiceHandle<art::TFileService> tfs;
75  //initialize beam histograms specified in fhicl file
76  art::TFileDirectory tffluxdir = tfs->mkdir( "Flux" );
77  string nutype[] ={ "nue", "nuebar", "numu", "numubar"};
78  string nultx[] ={"#nu_{e}", "#bar{#nu}_{e}", "#nu_{#mu}", "#bar{#nu}_{#mu}"};
79  string pltx[] ={"#mu^{#pm}","#pi^{#pm}","K^{0}_{L}","K^{#pm}"};
80  string secltx[] ={"pBe->#pi^{#pm}->...->#mu^{#pm}",
81  "pBe->#pi^{#pm}->..(not #mu^{#pm})..",
82  "pBe->K^{0}_{L}->...",
83  "pBe->K^{#pm}->...",
84  "pBe->(p or n)->..."};
85 
86 
87  int nbins=ps.get<int>("nBins",0);
88  int Elow=ps.get<float>("Elow",0);
89  int Ehigh=ps.get<float>("Ehigh",0);
90 
91  for (int i=0;i<4;i++) {
92  fHFlux[i]=tffluxdir.make<TH1D>(Form("h50%i",i+1),
93  Form("%s (all);Energy %s (GeV);#phi(%s)/50MeV/POT",nultx[i].c_str(),nultx[i].c_str(),nultx[i].c_str()),
94  nbins,Elow,Ehigh);
95  fHFlux[i]->Sumw2();
96  }
97 
98  for (int inu=0;inu<4;inu++) {
99  for (int ipar=0;ipar<4;ipar++) {
100  fHFluxParent[inu][ipar]=tffluxdir.make<TH1D>(Form("h5%i%i",ipar+1,inu+1),
101  Form("...->%s->%s;Energy %s (GeV);#phi(%s)/50MeV/POT",pltx[ipar].c_str(),nultx[inu].c_str(),nultx[inu].c_str(),nultx[inu].c_str()),
102  nbins,Elow,Ehigh);
103 
104  fHFluxParent[inu][ipar]->Sumw2();
105  }
106  for (int isec=0;isec<5;isec++) {
107  fHFluxSec[inu][isec]=tffluxdir.make<TH1D>(Form("h7%i%i",isec+1,inu+1),
108  Form("%s->%s;Energy %s (GeV);#phi(%s)/50MeV/POT",secltx[isec].c_str(),nultx[inu].c_str(),nultx[inu].c_str(),nultx[inu].c_str()),
109  nbins,Elow,Ehigh);
110  fHFluxSec[inu][isec]->Sumw2();
111  }
112  }
113 
114 
115  // TH1D* h=tffluxdir.make<TH1D>("numu","numu",240,0,120);
116  fEventCounter=0;
117  fEntry=ps.get<uint32_t>("skipEvents", 0);
118  fMaxEvents=ps.get<int>("maxEvents", -1);
119 
120  }
TH1D * fHFluxSec[4][5]
Definition: FluxReader.h:55
uint32_t fEntry
Definition: FluxReader.h:42
TH1D * fHFlux[4]
Definition: FluxReader.h:53
art::TypeLabel fTLdk2nu
Definition: FluxReader.h:61
std::string fInputType
Definition: FluxReader.h:45
art::SourceHelper const & fSourceHelper
Definition: FluxReader.h:38
TH1D * fHFluxParent[4][4]
Definition: FluxReader.h:54
art::SubRunID fSubRunID
Definition: FluxReader.h:39
fhicl::ParameterSet fConfigPS
Definition: FluxReader.h:57
art::TypeLabel fTLnuchoice
Definition: FluxReader.h:62
art::TypeLabel fTLmcflux
Definition: FluxReader.h:60
art::ServiceHandle< art::TFileService > tfs
art::TypeLabel fTLmctruth
Definition: FluxReader.h:59
uint32_t fEventCounter
Definition: FluxReader.h:41
art::RunNumber_t fIncrement
Definition: FluxReader.h:49
bool fSelfIncrementRuns
Definition: FluxReader.h:48

Member Function Documentation

void fluxr::FluxReader::closeCurrentFile ( )

Definition at line 122 of file FluxReader.cxx.

123  {
124  //mf::LogInfo(__FUNCTION__)<<"File boundary (processed "<<fEventCounter<<" events)"<<std::endl;
125  fSubRunID.flushSubRun();
126  fEventCounter=0;
127  fFluxInputFile->Close();
128  delete fFluxInputFile;
129  fSkipEvents=0; //if crossing file boundary don't skip events in next file
130  fEntry=0;
131  }
uint32_t fEntry
Definition: FluxReader.h:42
art::SubRunID fSubRunID
Definition: FluxReader.h:39
TFile * fFluxInputFile
Definition: FluxReader.h:52
uint32_t fSkipEvents
Definition: FluxReader.h:44
uint32_t fEventCounter
Definition: FluxReader.h:41
void fluxr::FluxReader::readFile ( std::string const &  name,
art::FileBlock *&  fb 
)

Definition at line 133 of file FluxReader.cxx.

135  {
136  std::cout << "readFile " << name << std::endl;
137  // Fill and return a new Fileblock.
138  fb = new art::FileBlock(art::FileFormatVersion(1, "FluxReader"), name);
139 
140  fFluxInputFile=new TFile(name.c_str());
141  if (fFluxInputFile->IsZombie()) {
142  //throw cet::exception(__PRETTY_FUNCTION__) << "Failed to open "<<fFluxInputFile<<std::endl;
143  }
144 
145  if (fInputType=="gsimple") {
146  fFluxDriver=new GSimpleInterface();
147  ((GSimpleInterface*)fFluxDriver)->SetRootFile(fFluxInputFile);
148  } else if (fInputType=="dk2nu") {
149  fFluxDriver=new DK2NuInterface();
150  ((DK2NuInterface*)fFluxDriver)->SetRootFile(fFluxInputFile);
151  ((DK2NuInterface*)fFluxDriver)->Init(fConfigPS);
152  } else if (fInputType=="boone") {
153  fFluxDriver=new BooNEInterface();
154  ((BooNEInterface*)fFluxDriver)->SetRootFile(fFluxInputFile);
155  } else {
156  throw cet::exception(__PRETTY_FUNCTION__) << "Ntuple format "<<fInputType<<" not supported"<<std::endl;
157  }
158  std::cout << "Reading in file " << name << std::endl;
159  std::cout<<"File has "<<fFluxDriver->GetEntries()<<" entries"<<std::endl;
160  std::cout<<"POT = "<<fFluxDriver->GetPOT()<<std::endl;
161  std::cout<<"Run = "<<fFluxDriver->GetRun()<<std::endl;
162  fPOT+=fFluxDriver->GetPOT();
164 
165  if (fSelfIncrementRuns) {
166  fIncrement ++;
168  }
169  }
virtual const float GetPOT()=0
std::string fInputType
Definition: FluxReader.h:45
virtual const int GetRun()=0
TFile * fFluxInputFile
Definition: FluxReader.h:52
fhicl::ParameterSet fConfigPS
Definition: FluxReader.h:57
virtual const Long64_t GetEntries()=0
virtual const void SetRun(int)=0
then echo fcl name
FluxInterface * fFluxDriver
Definition: FluxReader.h:51
BEGIN_PROLOG could also be cout
art::RunNumber_t fIncrement
Definition: FluxReader.h:49
bool fSelfIncrementRuns
Definition: FluxReader.h:48
bool fluxr::FluxReader::readNext ( art::RunPrincipal *const &  inR,
art::SubRunPrincipal *const &  inSR,
art::RunPrincipal *&  outR,
art::SubRunPrincipal *&  outSR,
art::EventPrincipal *&  outE 
)

Definition at line 172 of file FluxReader.cxx.

177  {
178  if (fMaxEvents > 0 && fEventCounter == unsigned(fMaxEvents))
179  return false;
180 
181  //if (fEventCounter%10000==0)
182  //mf::LogInfo(__FUNCTION__)<<"Attempting to read event: "<<fEventCounter<<std::endl;
183  // Create empty result, then fill it from current file:
184  std::unique_ptr< std::vector<simb::MCFlux> > mcfluxvec(new std::vector<simb::MCFlux >);
185  std::unique_ptr< std::vector<simb::MCTruth> > mctruthvec(new std::vector<simb::MCTruth >);
186 
187  std::unique_ptr< std::vector<bsim::Dk2Nu> > dk2nuvec(new std::vector<bsim::Dk2Nu >);
188  std::unique_ptr< std::vector<bsim::NuChoice> > nuchoicevec(new std::vector<bsim::NuChoice >);
189  std::unique_ptr< art::Assns<simb::MCTruth, bsim::Dk2Nu> >
190  dk2nuassn(new art::Assns<simb::MCTruth, bsim::Dk2Nu>);
191  std::unique_ptr< art::Assns<simb::MCTruth, bsim::NuChoice> >
192  nuchoiceassn(new art::Assns<simb::MCTruth, bsim::NuChoice>);
193  std::unique_ptr< art::Assns<simb::MCTruth, simb::MCFlux> >
194  mcfluxassn(new art::Assns<simb::MCTruth, simb::MCFlux>);
195 
196  simb::MCFlux flux;
197  if (!fFluxDriver->FillMCFlux(fEntry,flux))
198  return false;
199 
200  //check if neutrino goes through volTPCActive
201  //geo::GeometryCore const* geom = lar::providerFrom<geo::Geometry>();
202 
203  //std::string tmp[]={"volTPCActive"};
204  //std::set<std::string> volName(tmp,tmp+sizeof(tmp)/sizeof(tmp[0]));
205 
206  //std::vector< TGeoNode const * > geoNodeVec=geom->FindAllVolumes(volName);
207  //std::cout <<"vector with nodes is this long "<<geoNodeVec.size()<<std::endl;
208 
209  // std::cout<<fEventCounter<<std::endl;
210  //fake mctruth product to cheat eventweight that gets neutrino energy from it
211  simb::MCTruth mctruth;
212  simb::MCParticle mcpnu(0,flux.fntype,"Flux");
213  mcpnu.AddTrajectoryPoint(fFluxDriver->GetNuPosition(), fFluxDriver->GetNuMomentum());
214 
215  mctruth.Add(mcpnu);
216  mctruth.SetNeutrino(0,0,0,0,0,0,0,0,0,0);
217  mctruthvec->push_back(mctruth);
218 
219  if (fInputType=="dk2nu"){
220  dk2nuvec->push_back(*((DK2NuInterface*)fFluxDriver)->GetDk2Nu());
221  nuchoicevec->push_back(*((DK2NuInterface*)fFluxDriver)->GetNuChoice());
222  }
223 
224  int ipdg=0;
225  //select index matching order in nutype defined in constructor
226  if (flux.fntype==12) ipdg=0;
227  else if (flux.fntype==-12) ipdg=1;
228  else if (flux.fntype==14) ipdg=2;
229  else if (flux.fntype==-14) ipdg=3;
230 
231  double enu=flux.fnenergyn;
232  double totwgh=flux.fnwtnear*flux.fnimpwt/fCurrentPOT;
233  if (fInputType=="gsimple") totwgh=1./fCurrentPOT;//for gsimple files weight is actually 1
234  else if (totwgh<0 || totwgh>100) {
235  std::cout<<" Bad weight "<<totwgh<<std::endl;
236  totwgh=0;
237  }
238 
239  // if (ipdg == 2) std::cout << "Filling histo for 14 with weight " << totwgh << ", fCurrentPOT is " << fCurrentPOT << std::endl;
240  fHFlux[ipdg]->Fill(enu,totwgh);
241 
242  if (flux.fptype==13 || flux.fptype==-13) //mu+-
243  fHFluxParent[ipdg][0]->Fill(enu,totwgh);
244  else if (flux.fptype==211 || flux.fptype==-211) //pi+-
245  fHFluxParent[ipdg][1]->Fill(enu,totwgh);
246  else if (flux.fptype==130) //K0L
247  fHFluxParent[ipdg][2]->Fill(enu,totwgh);
248  else if (flux.fptype==321 || flux.fptype==-321) //K+-
249  fHFluxParent[ipdg][3]->Fill(enu,totwgh);
250 
251  if (fabs(flux.ftptype==211) && fabs(flux.fptype)==13)
252  fHFluxSec[ipdg][0]->Fill(enu,totwgh);
253  else if (fabs(flux.ftptype)==211)
254  fHFluxSec[ipdg][1]->Fill(enu,totwgh);
255  else if (fabs(flux.ftptype)==130)
256  fHFluxSec[ipdg][2]->Fill(enu,totwgh);
257  else if (fabs(flux.ftptype)==321)
258  fHFluxSec[ipdg][3]->Fill(enu,totwgh);
259  else if (flux.ftptype==2212 || flux.ftptype==2112)
260  fHFluxSec[ipdg][4]->Fill(enu,totwgh);
261 
262  mcfluxvec->push_back(flux);
263  fEventCounter++;
264  fEntry++;
265 
266  art::RunNumber_t rn = fFluxDriver->GetRun();
267  if (rn==0) rn=999999;
268  art::Timestamp tstamp(time(0));
269 
270  art::SubRunID newID(rn, 0); //subrun not used in flux files, so set to 0
271  if (fSubRunID.runID() != newID.runID()) { // New Run
272  outR = fSourceHelper.makeRunPrincipal(rn, tstamp);
273  }
274  if (fSubRunID != newID) { // New SubRun
275  outSR = fSourceHelper.makeSubRunPrincipal(rn,0,tstamp);
276  std::unique_ptr<sumdata::POTSummary> pot(new sumdata::POTSummary);
277  pot->totpot = fPOT;
278  pot->totgoodpot = fPOT;
279 
280  art::put_product_in_principal(std::move(pot),
281  *outSR,
282  "flux");
283 
284  fSubRunID = newID;
285  }
286 
287 
288  outE = fSourceHelper.makeEventPrincipal(fSubRunID.run(),
289  fSubRunID.subRun(),
291  tstamp);
292 
293  // Put products in the event.
294  art::put_product_in_principal(std::move(mcfluxvec),
295  *outE,
296  "flux"); // Module label
297  art::put_product_in_principal(std::move(mctruthvec),
298  *outE,
299  "flux"); // Module label
300  if (fInputType=="dk2nu") {
301  art::put_product_in_principal(std::move(dk2nuvec),
302  *outE,
303  "flux"); // Module label
304  art::put_product_in_principal(std::move(nuchoicevec),
305  *outE,
306  "flux"); // Module label
307 
308  auto aptr = fSourceHelper.makePtr<simb::MCTruth>(fTLmctruth, *outE, 0);
309  auto bptr = fSourceHelper.makePtr<bsim::Dk2Nu>(fTLdk2nu,*outE, 0);
310  auto cptr = fSourceHelper.makePtr<bsim::NuChoice>(fTLnuchoice, *outE, 0);
311  auto dptr = fSourceHelper.makePtr<simb::MCFlux>(fTLmcflux, *outE, 0);
312 
313  dk2nuassn->addSingle(aptr,bptr);
314  nuchoiceassn->addSingle(aptr,cptr);
315  mcfluxassn->addSingle(aptr,dptr);
316 
317  art::put_product_in_principal(std::move(dk2nuassn),
318  *outE,
319  "flux"); // Module label
320  art::put_product_in_principal(std::move(nuchoiceassn),
321  *outE,
322  "flux"); // Module label
323  art::put_product_in_principal(std::move(mcfluxassn),
324  *outE,
325  "flux"); // Module label
326  }
327 
328  return true;
329  }
virtual const TLorentzVector GetNuMomentum()=0
TH1D * fHFluxSec[4][5]
Definition: FluxReader.h:55
virtual bool FillMCFlux(Long64_t ientry, simb::MCFlux &mclux)=0
uint32_t fEntry
Definition: FluxReader.h:42
TH1D * fHFlux[4]
Definition: FluxReader.h:53
art::TypeLabel fTLdk2nu
Definition: FluxReader.h:61
std::string fInputType
Definition: FluxReader.h:45
virtual const TLorentzVector GetNuPosition()=0
virtual const int GetRun()=0
art::SourceHelper const & fSourceHelper
Definition: FluxReader.h:38
TH1D * fHFluxParent[4][4]
Definition: FluxReader.h:54
art::SubRunID fSubRunID
Definition: FluxReader.h:39
art::TypeLabel fTLnuchoice
Definition: FluxReader.h:62
art::TypeLabel fTLmcflux
Definition: FluxReader.h:60
art::TypeLabel fTLmctruth
Definition: FluxReader.h:59
FluxInterface * fFluxDriver
Definition: FluxReader.h:51
uint32_t fEventCounter
Definition: FluxReader.h:41
BEGIN_PROLOG could also be cout

Member Data Documentation

fhicl::ParameterSet fluxr::FluxReader::fConfigPS
private

Definition at line 57 of file FluxReader.h.

float fluxr::FluxReader::fCurrentPOT
private

Definition at line 47 of file FluxReader.h.

uint32_t fluxr::FluxReader::fEntry
private

Definition at line 42 of file FluxReader.h.

uint32_t fluxr::FluxReader::fEventCounter
private

Definition at line 41 of file FluxReader.h.

FluxInterface* fluxr::FluxReader::fFluxDriver
private

Definition at line 51 of file FluxReader.h.

TFile* fluxr::FluxReader::fFluxInputFile
private

Definition at line 52 of file FluxReader.h.

TH1D* fluxr::FluxReader::fHFlux[4]
private

Definition at line 53 of file FluxReader.h.

TH1D* fluxr::FluxReader::fHFluxParent[4][4]
private

Definition at line 54 of file FluxReader.h.

TH1D* fluxr::FluxReader::fHFluxSec[4][5]
private

Definition at line 55 of file FluxReader.h.

art::RunNumber_t fluxr::FluxReader::fIncrement = 1
private

Definition at line 49 of file FluxReader.h.

std::string fluxr::FluxReader::fInputType
private

Definition at line 45 of file FluxReader.h.

int fluxr::FluxReader::fMaxEvents
private

Definition at line 43 of file FluxReader.h.

float fluxr::FluxReader::fPOT
private

Definition at line 46 of file FluxReader.h.

bool fluxr::FluxReader::fSelfIncrementRuns
private

Definition at line 48 of file FluxReader.h.

uint32_t fluxr::FluxReader::fSkipEvents
private

Definition at line 44 of file FluxReader.h.

art::SourceHelper const& fluxr::FluxReader::fSourceHelper
private

Definition at line 38 of file FluxReader.h.

art::SubRunID fluxr::FluxReader::fSubRunID
private

Definition at line 39 of file FluxReader.h.

art::TypeLabel fluxr::FluxReader::fTLdk2nu
private

Definition at line 61 of file FluxReader.h.

art::TypeLabel fluxr::FluxReader::fTLmcflux
private

Definition at line 60 of file FluxReader.h.

art::TypeLabel fluxr::FluxReader::fTLmctruth
private

Definition at line 59 of file FluxReader.h.

art::TypeLabel fluxr::FluxReader::fTLnuchoice
private

Definition at line 62 of file FluxReader.h.


The documentation for this class was generated from the following files: