LBNEAnalysis.cc
Go to the documentation of this file.
1 //----------------------------------------------------------------------
2 // LBNEAnalysis.cc
3 //
4 // $Id: LBNEAnalysis.cc,v 1.3.2.6 2013/12/03 21:29:05 lebrun Exp $
5 //----------------------------------------------------------------------
6 #include <cmath>
7 #include <climits>
8 
9 #include <vector>
10 #include <fstream>
11 #include <iomanip>
12 #include <stdlib.h>
13 #include <map>
14 
15 //Root
16 #include <TSystem.h> // ROOT head file for a generic interface to the OS
17 #include <TStopwatch.h> // ROOT head file for stopwatch for real and cpu time
18 #include <TFile.h>
19 #include <TTree.h>
20 
21 //GEANT4
22 #include "globals.hh"
23 #include "G4ios.hh"
24 #include "G4Track.hh"
25 #include "G4SteppingManager.hh"
26 #include "G4ThreeVector.hh"
27 #include "G4TrajectoryContainer.hh"
28 #include "G4EventManager.hh"
29 #include "LBNERunManager.hh"
30 #include "G4ParticleDefinition.hh"
31 #include "G4ParticleTypes.hh"
32 #include "G4Navigator.hh"
33 #include "G4TransportationManager.hh"
34 #include "G4Run.hh"
35 #include "G4Version.hh"
36 #include "G4VProcess.hh"
37 #include "VersionAndContext.hh"
38 
39 //g4numi -> g4lbne..
40 #include "LBNEDataNtp_t.hh"
41 #include "LBNEParticleCode.hh"
42 #include "LBNEAnalysis.hh"
43 #include "LBNETrackInformation.hh"
44 #include "LBNEPrimaryGeneratorAction.hh"
45 #include "LBNENuWeight.hh"
46 #include "LBNEVolumePlacements.hh"
47 #include "LBNEDetectorConstruction.hh"
48 #include "LBNEQuickPiToNu.hh"
49 #include "dk2nu/tree/readWeightLocations.h"
50 #include "dk2nu/tree/calcLocationWeights.h"
51 
52 using namespace std;
53 
54 LBNEAnalysis* LBNEAnalysis::instance = 0;
55 
56 //------------------------------------------------------------------------------------
57 LBNEAnalysis::LBNEAnalysis() :
58  fDk2NuDetectorFileRead(false),
59  fHorn1TrackingTree(0),
60  fHorn2TrackingTree(0),
61  fDPTrackingTree(0),
62  fTargetOutputTree(0),
63  fOutFile(0),
64  fOutTree(0),
65  nuNtuple(0),
66  fOutFileDk2Nu(0),
67  fOutTreeDk2Nu(0),
68  fOutTreeDk2NuMeta(0),
69  fTrackingTree(0),
70  fDk2Nu(0),
71  fDkMeta(0),
72  fLBNEOutNtpData(0),
73  fTrackingPlaneData(0)
74 {
75  LBNERunManager *pRunManager=
76  dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
77 
78  if (pRunManager->GetVerboseLevel() > 0)
79  {
80  std::cout << "LBNEAnalysis Constructor Called." << std::endl;
81  }
82 
83 #ifdef G4ANALYSIS_USE
84 #endif
85 
86 
87  //g4data = new data_t();
88  fLBNEOutNtpData = new LBNEDataNtp_t();
89  fTrackingPlaneData = new LBNEDataNtp_t(); // Not sure what this does?
90 
91 
92  fcount = 0;
93  fentry = 0;
94 
95  //
96  //need code map defined but initialized
97  //or get a seg fault don't know why
98  //
99  /*code[-13] = 10;
100  code[13] = 11;
101  code[111] = 23;
102  code[211] = 13;
103  code[-211] = 14;
104  code[130] = 12;
105  code[321] = 15;
106  code[-321] = 16;
107  code[2112] = 8;
108  code[2212] = 1;
109  code[-2212] = 2;
110  code[310] = 19;
111  code[3122] = 17;
112  code[3222] = 21;
113  code[3212] = 22;
114  code[3112] = 20;*/
115 
116  setDetectorPositions();
117 
118 }
119 //------------------------------------------------------------------------------------
120 LBNEAnalysis::~LBNEAnalysis()
121 {
122  LBNERunManager *pRunManager=
123  dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
124 
125  if (pRunManager->GetVerboseLevel() > 0) {
126  std::cout << "LBNEAnalysis Destructor Called." << std::endl;
127  }
128 
129 #ifdef G4ANALYSIS_USE
130  // delete things
131 #endif
132 }
133 //------------------------------------------------------------------------------------
134 LBNEAnalysis* LBNEAnalysis::getInstance()
135 {
136  if (instance == 0) instance = new LBNEAnalysis;
137  return instance;
138 }
139 //------------------------------------------------------------------------------------
140 G4bool LBNEAnalysis::CreateOutput()
141 {
142 
143  LBNERunManager *pRunManager=
144  dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
145 
146  if (pRunManager->GetCreateOutput())
147  {
148  G4String spaces = " ";
149  std::cout << " LBNEAnalysis::CreateOutput() - Creating output ntuple..." << std::endl;
150 
151 
152  {
153  //std::cout << spaces << "Creating G4LBNEData Ntuple..." << std::endl;
154  //return LBNEAnalysis::CreateG4LBNEDataNtp();
155 
156  snprintf(nuNtupleFileName, 1023, "%s_%03d.root",
157  (pRunManager->GetOutputNtpFileName()).c_str(),pRunManager->GetCurrentRun()->GetRunID()); //April 29, 2015->Fixed the dk2nu extension (line 176)
158  fOutFile = new TFile(nuNtupleFileName,"RECREATE","root ntuple"); //which was giving the error due to conflict in name.
159  std::cerr << "Creating neutrino ntuple: " << nuNtupleFileName << std::endl;
160  //fOutTree = new TTree("nudata","g4numi Neutrino ntuple");
161  //fOutTree->Branch("data","data_t",&g4data,32000,1);
162  fOutTree = new TTree("nudata","g4lbne Neutrino ntuple");
163  fOutTree->Branch("data","LBNEDataNtp_t",&fLBNEOutNtpData,32000,1);
164 
165  if (pRunManager->GetCreateTrkPlaneOutput()) {
166  // create tree for tracking planes
167  std::cout << spaces << "Creating Tracking Plane Ntuple..." << std::endl;
168  fTrackingTree = new TTree("trackingdata", "Particles in Tracking Plane");
169  fTrackingTree->Branch("num", &fNParticles, "fNParticles/I");
170  fTrackingTree->Branch("trackID", &fTrackID, "fTrackID[fNParticles]/I]");
171  fTrackingTree->Branch("pdgCode", &fParticlePDG, "fParticlePDG[fNParticles]/I");
172  fTrackingTree->Branch("ImpWeight", &fNImpWt, "fNImpWt[fNParticles]/D");
173  fTrackingTree->Branch("positionX", &fParticleX, "fParticleX[fNParticles]/D");
174  fTrackingTree->Branch("positionY", &fParticleY, "fParticleY[fNParticles]/D");
175  fTrackingTree->Branch("positionZ", &fParticleZ, "fParticleZ[fNParticles]/D");
176  fTrackingTree->Branch("momentumX", &fParticleDX, "fParticleDX[fNParticles]/D");
177  fTrackingTree->Branch("momentumY", &fParticleDY, "fParticleDY[fNParticles]/D");
178  fTrackingTree->Branch("momentumZ", &fParticleDZ, "fParticleDZ[fNParticles]/D");
179  fTrackingTree->Branch("mass", &fParticleMass, "fParticleMass[fNParticles]/D");
180  fTrackingTree->Branch("energy", &fParticleEnergy, "fParticleEnergy[fNParticles]/D");
181  }
182  if (pRunManager->GetCreateTrkPlaneH1Output()) {
183  // create tree for Horn 1 tracking planes Amit Bashyal
184  std::cout << spaces << "Creating Horn1 Tracking Plane Ntuple..." << std::endl;
185  fHorn1TrackingTree = new TTree("H1trackingdata", "Particles in Horn1 Tracking Plane");
186  fHorn1TrackingTree->Branch("nH1P", &fH1NParticles, "fH1NParticles/I");
187  fHorn1TrackingTree->Branch("trackID", &fH1TrackID, "fH1TrackID/I]");
188  fHorn1TrackingTree->Branch("pdgCode", &fH1ParticlePDG, "fH1ParticlePDG/I");
189  fHorn1TrackingTree->Branch("H1ImpWeight", &fH1NImpWt, "fH1NImpWt/D");
190  fHorn1TrackingTree->Branch("H1positionX", &fH1ParticleX, "fH1ParticleX/D");
191  fHorn1TrackingTree->Branch("H1positionY", &fH1ParticleY, "fH1ParticleY/D");
192  fHorn1TrackingTree->Branch("H1positionZ", &fH1ParticleZ, "fH1ParticleZ/D");
193  fHorn1TrackingTree->Branch("H1momentumdirX", &fH1ParticleDX, "fH1ParticleDX/D");
194  fHorn1TrackingTree->Branch("H1momentumdirY", &fH1ParticleDY, "fH1ParticleDY/D");
195  fHorn1TrackingTree->Branch("H1momentumdirZ", &fH1ParticleDZ, "fH1ParticleDZ/D");
196  fHorn1TrackingTree->Branch("mass", &fH1ParticleMass, "fH1ParticleMass/D");
197  fHorn1TrackingTree->Branch("energy", &fH1ParticleEnergy, "fH1ParticleEnergy/D");
198  fHorn1TrackingTree->Branch("H1momentumX",&fH1ParticlePX,"fH1ParticlePX/D");
199  fHorn1TrackingTree->Branch("H1momentumY",&fH1ParticlePY,"fH1ParticlePY/D");
200  fHorn1TrackingTree->Branch("H1momentumZ",&fH1ParticlePZ,"fH1ParticlePZ/D");
201  fHorn1TrackingTree->Branch("H1PXPZ",&fH1ParticlePXPZ,"fH1ParticlePXPZ/D");
202  fHorn1TrackingTree->Branch("H1PYPZ",&fH1ParticlePYPZ,"fH1ParticlePYPZ/D");
203  fHorn1TrackingTree->Branch("H1PProductionX",&fH1PProductionX,"fH1PProductionX/D");
204  fHorn1TrackingTree->Branch("H1PProductionY",&fH1PProductionY,"fH1PProductionY/D");
205  fHorn1TrackingTree->Branch("H1PProductionZ",&fH1PProductionZ,"fH1PProductionZ/D");
206  fHorn1TrackingTree->Branch("H1PmomentumdirX",&fH1PProductionDX,"fH1PProductionDX/D");
207  fHorn1TrackingTree->Branch("H1PmomentmdirY",&fH1PProductionDY,"fH1PProductionDY/D");
208  fHorn1TrackingTree->Branch("H1PmomentumdirZ",&fH1PProductionDZ,"fH1PProductionDZ/D");
209  }
210 
211 
212  if (pRunManager->GetCreateTrkPlaneDPOutput()) {
213  // Create Planes for Decay Pipe End Tracking Plane Amit Bashyal
214  fDPTrackingTree = new TTree("DecayPipetrackingdata", "Particles Making to the Decay Pipe");
215  std::cerr << spaces << "Creating DecaypipeOutput Ntuples"<<std::endl;
216  fDPTrackingTree->Branch("DPImpWeight", &fDPNImpWt, "fDPNImpWt/D");
217  fDPTrackingTree->Branch("DPPProductionX",&fDPPProductionX,"fDPPProductionX/D");
218  fDPTrackingTree->Branch("DPPProductionY",&fDPPProductionY,"fDPPProductionY/D");
219  fDPTrackingTree->Branch("DPPProductionZ",&fDPPProductionZ,"fDPPProductionZ/D");
220  fDPTrackingTree->Branch("DPPtrackID", &fDPTrackID, "fDPTrackID/I]");
221  fDPTrackingTree->Branch("DPParentID", &fDPParentID, "fDPParentID/I]");
222  fDPTrackingTree->Branch("DPenergy", &fDPParticleEnergy, "fDPParticleEnergy/D");
223  fDPTrackingTree->Branch("DPpositionX", &fDPParticleX, "fDPParticleX/D");
224  fDPTrackingTree->Branch("DPpositionY", &fDPParticleY, "fDPParticleY/D");
225  fDPTrackingTree->Branch("DPpositionZ", &fDPParticleZ, "fDPParticleZ/D");
226  fDPTrackingTree->Branch("DPmomentumX",&fDPParticlePX,"fDPParticlePX/D");
227  fDPTrackingTree->Branch("DPmomentumY",&fDPParticlePY,"fDPParticlePY/D");
228  fDPTrackingTree->Branch("DPmomentumZ",&fDPParticlePZ,"fDPParticlePZ/D");
229  fDPTrackingTree->Branch("DPpdgCode", &fDPParticlePDG, "fDPParticlePDG/I");
230  }
231  if (pRunManager->GetCreateTargetOutput()) {
232  std::cout<<spaces<<"Creating TargetOutPut Ntuples"<<std::endl;
233  fTargetOutputTree = new TTree("TargetOutputData", "Particles Produced in the Target");
234  fTargetOutputTree->Branch("TtrackID", &fTTrackID, "fTTrackID/I]");
235  fTargetOutputTree->Branch("TpdgCode", &fTParticlePDG, "fTParticlePDG/I");
236  fTargetOutputTree->Branch("TpositionX", &fTParticleX, "fTParticleX/D");
237  fTargetOutputTree->Branch("TpositionY", &fTParticleY, "fTParticleY/D");
238  fTargetOutputTree->Branch("TpositionZ", &fTParticleZ, "fTParticleZ/D");
239  fTargetOutputTree->Branch("Tenergy", &fTParticleEnergy, "fTParticleEnergy/D");
240  fTargetOutputTree->Branch("TmomentumX",&fTParticlePX,"fTParticlePX/D");
241  fTargetOutputTree->Branch("TmomentumY",&fTParticlePY,"fTParticlePY/D");
242  fTargetOutputTree->Branch("TmomentumZ",&fTParticlePZ,"fTParticlePZ/D");
243 
244  }
245 
246  // fOutTree->Branch("horn2data","LBNEDataNtp_t",&fLBNEOutNtpData,32000,1);
247  if (pRunManager->GetCreateTrkPlaneH2Output()) {
248  // create tree for Horn 2 tracking planes
249  std::cout << spaces << "Creating Horn2 Tracking Plane Ntuple..." << std::endl;
250  fHorn2TrackingTree = new TTree("H2trackingdata", "Particles in Horn2 Tracking Plane");
251  fHorn2TrackingTree->Branch("nH2P", &fH2NParticles, "fH2NParticles/I");
252  fHorn2TrackingTree->Branch("trackID", &fH2TrackID, "fH2TrackID/I]");
253  fHorn2TrackingTree->Branch("pdgCode", &fH2ParticlePDG, "fH2ParticlePDG/I");
254  fHorn2TrackingTree->Branch("H2ImpWeight", &fH2NImpWt, "fH2NImpWt/D");
255  fHorn2TrackingTree->Branch("H2positionX", &fH2ParticleX, "fH2ParticleX/D");
256  fHorn2TrackingTree->Branch("H2positionY", &fH2ParticleY, "fH2ParticleY/D");
257  fHorn2TrackingTree->Branch("H2positionZ", &fH2ParticleZ, "fH2ParticleZ/D");
258  fHorn2TrackingTree->Branch("H2momentumdirX", &fH2ParticleDX, "fH2ParticleDX/D");
259  fHorn2TrackingTree->Branch("H2momentumdirY", &fH2ParticleDY, "fH2ParticleDY/D");
260  fHorn2TrackingTree->Branch("H2momentumdirZ", &fH2ParticleDZ, "fH2ParticleDZ/D");
261  fHorn2TrackingTree->Branch("mass", &fH2ParticleMass, "fH2ParticleMass/D");
262  fHorn2TrackingTree->Branch("energy", &fH2ParticleEnergy, "fH2ParticleEnergy/D");
263  fHorn2TrackingTree->Branch("H2momentumX",&fH2ParticlePX,"fH2ParticlePX/D");
264  fHorn2TrackingTree->Branch("H2momentumY",&fH2ParticlePY,"fH2ParticlePY/D");
265  fHorn2TrackingTree->Branch("H2momentumZ",&fH2ParticlePZ,"fH2ParticlePZ/D");
266  fHorn2TrackingTree->Branch("H2PXPZ",&fH2ParticlePXPZ,"fH2ParticlePXPZ/D");
267  fHorn2TrackingTree->Branch("H2PYPZ",&fH2ParticlePYPZ,"fH2ParticlePYPZ/D");
268  fHorn2TrackingTree->Branch("H2PProductionX",&fH2PProductionX,"fH2PProductionX/D");
269  fHorn2TrackingTree->Branch("H2PProductionY",&fH2PProductionY,"fH2PProductionY/D");
270  fHorn2TrackingTree->Branch("H2PProductionZ",&fH2PProductionZ,"fH2PProductionZ/D");
271  fHorn2TrackingTree->Branch("H2PmomentumdirX",&fH2PProductionDX,"fH2PProductionDX/D");
272  fHorn2TrackingTree->Branch("H2PmomentumdirY",&fH2PProductionDY,"fH2PProductionDY/D");
273  fHorn2TrackingTree->Branch("H2PmomentumdirZ",&fH2PProductionDZ,"fH2PProductionDZ/D");
274  fHorn2TrackingTree->Branch("H2CNProcess",&fH2CNProcess,"fH2CNProcess/I");
275 
276  }
277 
278  if (pRunManager->GetCreateAlcoveTrackingOutput()){
279  std::cout << spaces << "Creating Sculpted Absorber Alcove Tracking Plane Ntuple..." << std::endl;
280  fAlcoveTrackingTree = new TTree("AlcoveTracks","Alcove Tracking");
281 
282  fAlcoveTrackingTree->Branch("run",&fMuRunNo,"run/I");
283  fAlcoveTrackingTree->Branch("event",&fMuEvtNo,"event/I");
284  fAlcoveTrackingTree->Branch("np",&fMuNParticles,"np/I");
285  fAlcoveTrackingTree->Branch("ID",&fMuTrackID,"ID[np]/I");
286  fAlcoveTrackingTree->Branch("ParID",&fMuParentID,"ParID[np]/I");
287  fAlcoveTrackingTree->Branch("PDG",&fMuPDG,"PDG[np]/I");
288  fAlcoveTrackingTree->Branch("impwt",&fMuNimpWt,"impwt[np]/D");
289  fAlcoveTrackingTree->Branch("x",&fMuX,"x[np]/D");
290  fAlcoveTrackingTree->Branch("y",&fMuY,"y[np]/D");
291  fAlcoveTrackingTree->Branch("z",&fMuZ,"z[np]/D");
292  fAlcoveTrackingTree->Branch("startx",&fMuStartX,"startx[np]/D");
293  fAlcoveTrackingTree->Branch("starty",&fMuStartY,"starty[np]/D");
294  fAlcoveTrackingTree->Branch("startz",&fMuStartZ,"startz[np]/D");
295  fAlcoveTrackingTree->Branch("startE",&fMuStartE,"startE[np]/D");
296  fAlcoveTrackingTree->Branch("px",&fMuPX,"px[np]/D");
297  fAlcoveTrackingTree->Branch("py",&fMuPY,"py[np]/D");
298  fAlcoveTrackingTree->Branch("pz",&fMuPZ,"pz[np]/D");
299  fAlcoveTrackingTree->Branch("m",&fMuMass,"m[np]/D");
300  fAlcoveTrackingTree->Branch("E",&fMuEnergy,"E[np]/D");
301  fAlcoveTrackingTree->Branch("dEdx",&fMudEdx,"dEdx[np]/D");
302  fAlcoveTrackingTree->Branch("dEdx_ion",&fMudEdx_ion,"dEdx_ion[np]/D");
303  fAlcoveTrackingTree->Branch("cosTheta",&fMuTheta,"cosTheta[np]/D");
304  fAlcoveTrackingTree->Branch("edep",&fMuDE,"edep[np]/D");
305  fAlcoveTrackingTree->Branch("edep_ion",&fMuDEion,"edep_ion[np]/D");
306  fAlcoveTrackingTree->Branch("nsteps",&fMuNSteps,"nsteps[np]/I");
307  fAlcoveTrackingTree->Branch("parPDG",&fMuParPDG,"parPDG[np]/I");
308  fAlcoveTrackingTree->Branch("parE",&fMuParE,"parE[np]/D");
309  fAlcoveTrackingTree->Branch("parX",&fMuParX,"parX[np]/D");
310  fAlcoveTrackingTree->Branch("parY",&fMuParY,"parY[np]/D");
311  fAlcoveTrackingTree->Branch("parZ",&fMuParZ,"parZ[np]/D");
312  fAlcoveTrackingTree->Branch("parPX",&fMuParPX,"parPX[np]/D");
313  fAlcoveTrackingTree->Branch("parPY",&fMuParPY,"parPY[np]/D");
314  fAlcoveTrackingTree->Branch("parPZ",&fMuParPZ,"parPZ[np]/D");
315 
316 
317  fMuNParticles = 0;//ensure that we start with no particles
318  }//Sculpted absorber tracking plane
319 
320  return true;
321 
322  }
323 
324  return false;
325 
326  }
327 
328  return false;
329 
330 
331 }
332 // Clone from the above, but for Dk2Nu type of output
333 
334 G4bool LBNEAnalysis::CreateDk2NuOutput()
335 {
336 
337  LBNERunManager *pRunManager=
338  dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
339 
340  if (pRunManager->GetCreateDk2NuOutput())
341  {
342  G4String spaces = " ";
343  std::cout << " LBNEAnalysis::CreateOutput() - Creating output ntuple..." << std::endl;
344 
345 
346  {
347  //std::cout << spaces << "Creating G4LBNEData Ntuple..." << std::endl;
348  //return LBNEAnalysis::CreateG4LBNEDataNtp();
349 
350  pRunManager = dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
351  fDk2Nu = new bsim::Dk2Nu;
352  fDkMeta = new bsim::DkMeta;
353  snprintf(nuNtupleFileName, 1023, "%s_%05d.dk2nu.root",
354  (pRunManager->GetOutputDk2NuFileName()).c_str(),pRunManager->GetCurrentRun()->GetRunID());
355  fOutFileDk2Nu = new TFile(nuNtupleFileName,"RECREATE","root ntuple");
356  std::cerr << "Creating neutrino ntuple, Dk2Nu : " << nuNtupleFileName << std::endl;
357 
358  fOutTreeDk2Nu = new TTree("dk2nuTree", "g4lbne neutrino ntuple, dk2Nu format");
359  fOutTreeDk2Nu->Branch("dk2nu", "bsim::Dk2Nu", &fDk2Nu, 32000, 99);
360  fOutTreeDk2NuMeta = new TTree("dkmetaTree", "g4lbne neutrino ntuple, dk2Nu format, metadata");
361  fOutTreeDk2NuMeta->Branch("dkmeta", "bsim::DkMeta", &fDkMeta, 32000,99);
362  std::cerr << "Done Creating neutrino ntuple, Dk2Nu : " << nuNtupleFileName << std::endl;
363 
364  return true;
365 
366  }
367 
368  return false;
369 
370  }
371 
372  return false;
373 
374 
375 }
376 
377 void LBNEAnalysis::fillDkMeta() {
378 
379  VersionAndContext* versionContext = VersionAndContext::Instance();
380 
381  LBNERunManager* theRunManager = dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
382  if (!theRunManager->GetCreateDk2NuOutput()) return;
383 
384  fDkMeta->job = theRunManager->GetCurrentRun()->GetRunID();
385  fDkMeta->pots = theRunManager->GetNumberOfEventsLBNE();
386  fDkMeta->beamsim = std::string("g4lbnf-");
387  fDkMeta->beamsim += versionContext->GetGitDescription();
388  fDkMeta->beamsim += std::string (", builddir=") +
389  versionContext->GetG4lbnfDir();
390  fDkMeta->beamsim += std::string (", on") +
391  versionContext->GetBuildTimeStamp();
392  fDkMeta->physics = G4Version;
393  fDkMeta->physcuts = theRunManager->GetPhysicsListName();
394  fDkMeta->tgtcfg = "le00zdu";
395 
396  {
397  std::ostringstream oStrStr;
398  oStrStr << " --Kine Cut " << theRunManager->GetLBNESteppingManager()->GetKillTrackingThreshold()/CLHEP::GeV <<"CLHEP::GeV";
399  fDkMeta->physcuts += oStrStr.str();
400  }
401 
402  LBNEVolumePlacements *plManager = LBNEVolumePlacements::Instance();
403  {
404  //if (plManager->GetUse1p2MW()) fDkMeta->tgtcfg += std::string("1p2MW");
405  //else fDkMeta->tgtcfg += std::string("700kW");
406  std::ostringstream oStrStr;
407  oStrStr.precision(4);
408  oStrStr << "-length" << (plManager->GetTargetSLengthGraphite())/CLHEP::cm;
409  oStrStr << "-Pos" << (plManager->GetTargetSLengthGraphite() -
410  plManager->GetTargetLengthIntoHorn() + 25.3*CLHEP::mm)/CLHEP::cm;
411  oStrStr << "cm";
412  fDkMeta->tgtcfg = "le00zdu";//versionContext->GetMacroFileName();
413  }
414  {
415  std::ostringstream oStrStr;
416  oStrStr.precision(4);
417  const LBNEDetectorConstruction *pDet =
418  static_cast<const LBNEDetectorConstruction*> (theRunManager->GetUserDetectorConstruction());
419  oStrStr << "Current-" << pDet->GetHornCurrent()/(CLHEP::ampere*1000.); //in kAmp
420  // fDkMeta->horncfg = oStrStr.str();
421  oStrStr << "Current-" << pDet->GetHornCurrent()/(CLHEP::ampere*1000.); //in kAmp
422  // fDkMeta->horncfg = std::string("230i");
423  fDkMeta->horncfg = pDet->GetHornCurrent()/(CLHEP::ampere*1000.);
424  }
425  // fDkMeta->dkvolcfg = plManager->GetDecayPipeGas();
426  std::cout<<"GetDecayPipeGas() "<<plManager->GetDecayPipeGas()<<std::endl;
427  fDkMeta->dkvolcfg = std::string("helium");
428  {
429  std::ostringstream oStrStr;
430  oStrStr.precision(4);
431  oStrStr << "length-" << plManager->GetDecayPipeLength()/CLHEP::cm << "cm";
432  // fDkMeta->dkvolcfg += oStrStr.str();
433  }
434  //
435  // Beam parameters
436  //
437  const LBNEPrimaryGeneratorAction* NPGA=
438  static_cast<const LBNEPrimaryGeneratorAction*> (theRunManager->GetUserPrimaryGeneratorAction());
439  fDkMeta->beam0x = NPGA->GetBeamOffsetX();
440  fDkMeta->beam0y = NPGA->GetBeamOffsetY();
441  fDkMeta->beam0z = NPGA->GetBeamOffsetZ();
442  fDkMeta->beamhwidth = NPGA->GetBeamSigmaX();
443  fDkMeta->beamvwidth = NPGA->GetBeamSigmaY();
444  fDkMeta->beamdxdz = 0.;
445  fDkMeta->beamdydz = NPGA->GetBeamAngleTheta();
446 
447  //A. Bashyal stuff for the ppfx:
448  G4String name_gen[] = {"parent","granparent","greatgranparent"};
449  G4String name_vol[] = {"PHorn1IC","PHorn2IC","DPIP","DVOL"}; //First we try to verify if this will work on ppfx.
450  for(G4int ii=0;ii<nGenAbs;ii++){
451  GenAbsName.push_back(name_gen[ii]);
452  }
453  for(G4int ii=0;ii<nVolAbs;ii++){
454  VolAbsName.push_back(name_vol[ii]);
455  }
456 // char namevar[50];
457  for(int ii=0; ii<nVolAbs; ii++){
458  for(int jj=0; jj<nGenAbs; jj++){
459 // sprintf(namevar,"Material_%s_%s",VolAbsName[ii].c_str(),GenAbsName[jj].c_str());
460 // VolVdblName.push_back(G4String(namevar));
461 // recoding this in C++, potential overwrite with fix C-style arrays.
462  std::ostringstream nameVarStrStr; nameVarStrStr << "Material_" << VolAbsName[ii] << "_" << GenAbsName[jj];
463  const std::string nameVarStr(nameVarStrStr.str());
464  VolVdblName.push_back(nameVarStr);
465  }
466  }
467  nVdblTot = int(VolVdblName.size());
468 
469  //For vint branch
470  nVintTot = 2;
471  VolVintName.clear();
472  G4String name_vint[] = {"Index_Tar_In_Ancestry","playlistID"};
473  for(G4int ii=0;ii<nVintTot;ii++){
474  VolVintName.push_back(name_vint[ii]);
475  }
476  fDkMeta->vintnames.clear();
477  for(int ii=0;ii<nVintTot;ii++){
478  (fDkMeta->vintnames).push_back(VolVintName[ii]);
479  std::cout<< " LBNEAnalysis:: nVintnames are : "<<VolVintName[ii]<<std::endl;
480  }
481  fDkMeta->vdblnames.clear();
482  for(int ii=0;ii<nVdblTot;ii++){
483  (fDkMeta->vdblnames).push_back(VolVdblName[ii]);
484  std::cout<<" LBNEAnalysis: nVdblTot are : "<<VolVdblName[ii]<<std::endl;
485  }
486 
487  // Detector positions
488  if (!fDk2NuDetectorFileRead) {
489  fDkMeta->location.clear();
490  fDk2NuDetectorFileRead=true;
491  std::string aFileName(theRunManager->GetDetectorLocationFileName());
492  if (aFileName == G4String("?")) {
493  char *G4WORKDIRC = getenv("G4LBNE_DIR");
494  if (G4WORKDIRC == 0) {
495  G4WORKDIRC = getenv("G4LBNE_DIR");
496  if (G4WORKDIRC == 0) {
497  G4Exception("LBNEAnalysis::LBNEAnalysis::fillDkMeta", "InvalidSetup",
498  FatalErrorInArgument,
499  "Environmental variable G4LBNE_DIR or G4LBNE not set, can't find default Det. loc. file ");
500  }
501  }
502  aFileName = std::string(G4WORKDIRC) + std::string("/locations/etc/LBNElocations.txt");
503  std::cerr << " ... And this file name is ... " << aFileName << std::endl;
504  }
505  // Is this open?
506  std::ifstream fStr;
507  fStr.open(aFileName.c_str());
508  if (!fStr.is_open()) {
509  std::cerr << " LBNEAnalysis::LBNEAnalysis::fillDkMeta, with DK2NU set to " << std::string(getenv("G4LBNE_DIR")) << std::endl;
510  std::cerr << " can't open detector location file ... " << aFileName << std::endl;
511  G4Exception("LBNEAnalysis::LBNEAnalysis::fillDkMeta", "InvalidSetup",
512  FatalErrorInArgument, "Invalid Detector FileName (not found or corrupted)");
513  } else {
514  fStr.close();
515  }
516  bsim::readWeightLocations(aFileName, fDkMeta);
517  }
518  fOutTreeDk2NuMeta->Fill();
519 }
520 
521 //------------------------------------------------------------------------------------
522 void LBNEAnalysis::CloseOutput()
523 {
524  LBNERunManager *pRunManager=
525  dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
526  if (pRunManager->GetCreateOutput()) {
527  //if(fOutFile && fOutFile->IsOpen())
528  //{
529  fOutFile->cd();
530  fOutTree->Write();
531 
532  if (pRunManager->GetCreateTrkPlaneOutput()){
533  //---tracking plane data
534  std::cout << " Writing out Tracking Plane Ntuple to " << fOutFile -> GetName() << std::endl;
535  fTrackingTree->Write();
536  }
537 
538 
539 
540  if(pRunManager->GetCreateTrkPlaneH1Output()){
541  //---tracking plane data
542  std::cout << " Writing out Tracking Plane Ntuple to " << fOutFile -> GetName() << std::endl;
543  if (fHorn1TrackingTree != 0) fHorn1TrackingTree->Write();
544  }
545 
546  if (pRunManager->GetCreateTrkPlaneH2Output()){
547  //---tracking plane data
548  std::cout << " Writing out Tracking Plane Ntuple to " << fOutFile -> GetName() << std::endl;
549  if (fHorn2TrackingTree != 0) fHorn2TrackingTree->Write();
550  }
551  if
552  (pRunManager->GetCreateTrkPlaneDPOutput()){
553  //---tracking plane data
554  std::cout << " Writing out Tracking Plane Ntuple to " << fOutFile -> GetName() << std::endl;
555  if (fDPTrackingTree != 0) fDPTrackingTree->Write();
556  }
557 
558  if (pRunManager->GetCreateTargetOutput()){
559  //---tracking plane data
560  std::cout << " Writing out Tracking Plane Ntuple to " << fOutFile -> GetName() << std::endl;
561  if (fTargetOutputTree != 0) fTargetOutputTree->Write();
562  }
563 
564  if (pRunManager->GetCreateAlcoveTrackingOutput())
565  {
566  std::cout << " Writeing out sculpted absorber tracking plane to "<<fOutFile->GetName()<<std::endl;
567  if (fAlcoveTrackingTree != 0 ) fAlcoveTrackingTree->Write();
568 
569  }
570 
571 
572  fOutFile->Close();
573 
574 
575 /*
576  if(fOutFile->IsOpen())
577  {
578  std::cout << " PROBLEM: Failed to close Output Ntuple " << fOutFile -> GetName() << std::endl;
579  }
580  else
581  {*/
582  std::cout << " Sucessfully closed Output Ntuple " << fOutFile -> GetName() << std::endl;
583 /* }
584 */
585  delete fOutFile;
586  //}
587  }
588 
589 }
590 
591 void LBNEAnalysis::CloseDk2NuOutput()
592 {
593  LBNERunManager *pRunManager=
594  dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
595  if (pRunManager->GetCreateDk2NuOutput()) {
596  fOutFileDk2Nu->cd();
597  std::cerr << " Writing the Dk2Nu Ntuple ... " << std::endl;
598  fOutTreeDk2Nu->Write();
599  fOutTreeDk2NuMeta->Write();
600  std::cerr << " Written the Dk2Nu Ntuple and meta data.. ... " << std::endl;
601 
602  fOutFileDk2Nu->Close();
603  delete fOutFileDk2Nu;
604  //}
605  }
606 
607 }
608 
609 
610 
611 //------------------------------------------------------------------------------------
612 void LBNEAnalysis::SetCount(G4int count)
613 {
614  fcount = count;
615 }
616 //------------------------------------------------------------------------------------
617 G4int LBNEAnalysis::GetCount()
618 {
619  return fcount;
620 }
621 //------------------------------------------------------------------------------------
622 void LBNEAnalysis::SetEntry(G4int entry)
623 {
624  fentry = entry;
625 }
626 //------------------------------------------------------------------------------------
627 G4int LBNEAnalysis::GetEntry()
628 {
629  return fentry;
630 }
631 
632 
633 //------------------------------------------------------------------------------------
634 void LBNEAnalysis::FillNeutrinoNtuple(const G4Track& track,const
635 std::vector<G4VTrajectory*>& history)
636 {
637 
638  LBNERunManager *pRunManager=
639  dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
640  const LBNERunAction *runUserAction = reinterpret_cast<const LBNERunAction *>(pRunManager->GetUserRunAction());
641 
642  if (runUserAction->GetDoComputeEDepInHorns()) return;
643 
644  if (pRunManager->GetVerboseLevel() > 3)
645  { G4cout << "LBNEAnalysis::FillNeutrinoNtuple() called." << G4endl;}
646 
647  if ((!pRunManager->GetCreateOutput()) && (!pRunManager->GetCreateDk2NuOutput())) return;
648 
649  if (pRunManager->GetCreateOutput()) fLBNEOutNtpData -> Clear();
650 
651 
652  LBNEPrimaryGeneratorAction *NPGA = (LBNEPrimaryGeneratorAction*)(pRunManager)->GetUserPrimaryGeneratorAction();
653  G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
654 
655 
656  //Neutrino vertex position and momentum
657  G4ThreeVector pos = track.GetPosition()/CLHEP::mm;
658  bool validTrack = true;
659  if (std::isnan(pos.x()) || std::isnan(pos.y()) || std::isnan(pos.z())) {
660  std::cerr << " LBNEAnalysis::FillNeutrinoNtuple Bogus track position for track ID "
661  << track.GetTrackID() << " parentID " << track.GetParentID() << std::endl;
662  validTrack = false;
663  return;
664  }
665  x = pos.x();
666  y = pos.y();
667  z = pos.z();
668  NuMomentum = track.GetMomentum();
669  G4int parentID = track.GetParentID();
670 
671  if (pRunManager->GetCreateOutput()) fLBNEOutNtpData->ptrkid = parentID;
672 
673  if (parentID == 0)
674  { std::cout << "Event no: " << evtno << " LBNEAnalysis::FillNeutrinoNtuple - "
675  << "PROBLEM: Direct Nu Parent has track ID = 0." << std::endl;
676  return; //I have to make some changes so that neutrinos in fluka/mars ntuples don't crash
677  }
678 
679  // Debugging anomalous abosrption here...
680  if(LBNEQuickPiToNuVect::Instance()->doIt()) {
681  LBNEQuickPiToNuVect::Instance()->blessPion(parentID, track.GetTotalEnergy(), track.GetMomentum());
682  }
683 
684  LBNETrajectory* NuParentTrack = GetParentTrajectory(parentID);
685  G4int point_no = NuParentTrack->GetPointEntries();
686  G4ThreeVector ParentMomentumFinal = NuParentTrack->GetMomentum(point_no-1);
687  G4ThreeVector vertex_r = NuParentTrack->GetPoint(point_no-1)->GetPosition(); //Should be the same as Neutrino vertex
688  G4String parent_name = NuParentTrack->GetParticleName();
689  G4double Parent_mass = NuParentTrack->GetMass();
690  G4double gamma = sqrt(ParentMomentumFinal*ParentMomentumFinal+Parent_mass*Parent_mass)/Parent_mass;
691  G4double Parent_energy = gamma*Parent_mass;
692  G4ThreeVector beta_vec = ParentMomentumFinal/Parent_energy;
693  G4double partial = gamma*(beta_vec*NuMomentum);
694  G4double enuzr = gamma*(track.GetTotalEnergy())-partial; //neutrino energy in parent rest frame
695  G4double enuzrInGeV = enuzr/CLHEP::GeV; //neutrino energy in parent rest frame, in CLHEP::GeV
696  //fill histograms, ntuples,...
697  if (!validTrack) {
698 
699  std::cerr << " More info about bogus track... Parent name " << parent_name << " Mom "
700  << ParentMomentumFinal.x() << " / " << ParentMomentumFinal.y() << " / " << ParentMomentumFinal.z()
701  << std::endl;
702 
703  }
704  if (pRunManager->GetCreateOutput()) {
705  fLBNEOutNtpData->protonN = NPGA->GetProtonNumber();
706 
707  fLBNEOutNtpData->run = pRunManager->GetCurrentRun()->GetRunID();
708  fLBNEOutNtpData->evtno = pRunManager->GetCurrentEvent()->GetEventID();
709  fLBNEOutNtpData->beamHWidth = NPGA->GetBeamSigmaX()/CLHEP::cm;
710  fLBNEOutNtpData->beamVWidth = NPGA->GetBeamSigmaY()/CLHEP::cm;
711  fLBNEOutNtpData->beamX = NPGA->GetBeamOffsetX()/CLHEP::cm;
712  fLBNEOutNtpData->beamY = NPGA->GetBeamOffsetY()/CLHEP::cm;
713  }
714  //G4int particleID = track.GetParentID();
715  G4ThreeVector protonOrigin = NPGA->GetProtonOrigin();
716  if (pRunManager->GetCreateOutput()) {
717  fLBNEOutNtpData->protonX = protonOrigin[0];
718  fLBNEOutNtpData->protonY = protonOrigin[1];
719  fLBNEOutNtpData->protonZ = protonOrigin[2];
720  }
721 
722  G4ThreeVector protonMomentum = NPGA->GetProtonMomentum();
723  if (pRunManager->GetCreateOutput()) {
724  fLBNEOutNtpData->protonPx = protonMomentum[0];
725  fLBNEOutNtpData->protonPy = protonMomentum[1];
726  fLBNEOutNtpData->protonPz = protonMomentum[2];
727  }
728 // std::cerr << " protonOrigin " << protonOrigin << " momentum " << protonMomentum << std::endl;
729 
730  LBNEVolumePlacements *volDB = LBNEVolumePlacements::Instance();
731  if (pRunManager->GetCreateOutput()) {
732  fLBNEOutNtpData->nuTarZ = volDB->GetTargetLengthIntoHorn(); // A better info that the somewhat meaningless Z0
733  }
734  const LBNEDetectorConstruction *detDB =
735  dynamic_cast<const LBNEDetectorConstruction*>(pRunManager->GetUserDetectorConstruction());
736  //other info
737  // Neutrino origin:
738  // 3 From muon decay
739  // 1 From particle from target
740  // 2 From scraping
741  //check if nu is from muon decay or from a particle from target, otherwise Norig = 2
742  G4int Norig = 2;
743  if ((parent_name=="mu+") || (parent_name=="mu-")) Norig = 3;
744  G4String firstvolname = NuParentTrack->GetPreStepVolumeName(0);
745  if (firstvolname.contains("Baffle") || firstvolname.contains("TGT")) Norig = 1;
746 
747  if (pRunManager->GetCreateOutput()) {
748  fLBNEOutNtpData->hornCurrent = detDB->GetHornCurrent()/CLHEP::ampere/1000.;
749 
750  // Random decay - these neutrinos rarely hit any of the detectors
751  fLBNEOutNtpData->Ndxdz = NuMomentum[0]/NuMomentum[2];
752  fLBNEOutNtpData->Ndydz = NuMomentum[1]/NuMomentum[2];
753  fLBNEOutNtpData->Npz = NuMomentum[2]/CLHEP::GeV;
754  fLBNEOutNtpData->Nenergy = track.GetTotalEnergy()/CLHEP::GeV;
755 
756 
757  fLBNEOutNtpData->Norig = Norig;
758  fLBNEOutNtpData->Ndecay = NuParentTrack->GetDecayCode();
759 
760  G4ParticleDefinition * particleType = track.GetDefinition();
761  G4int ntype = LBNEParticleCode::AsInt(LBNEParticleCode::StringToEnum(particleType->GetParticleName()));
762  fLBNEOutNtpData->Ntype = ntype;
763  fLBNEOutNtpData->Vx = x/CLHEP::cm;
764  fLBNEOutNtpData->Vy = y/CLHEP::cm;
765  fLBNEOutNtpData->Vz = z/CLHEP::cm;
766  fLBNEOutNtpData->pdPx = ParentMomentumFinal[0]/CLHEP::GeV;
767  fLBNEOutNtpData->pdPy = ParentMomentumFinal[1]/CLHEP::GeV;
768  fLBNEOutNtpData->pdPz = ParentMomentumFinal[2]/CLHEP::GeV;
769 
770  G4ThreeVector ParentMomentumProduction = NuParentTrack->GetMomentum(0);
771  fLBNEOutNtpData->ppdxdz = ParentMomentumProduction[0]/ParentMomentumProduction[2];
772  fLBNEOutNtpData->ppdydz = ParentMomentumProduction[1]/ParentMomentumProduction[2];
773  fLBNEOutNtpData->pppz = ParentMomentumProduction[2]/CLHEP::GeV;
774 
775  G4double parentp = sqrt(ParentMomentumProduction*ParentMomentumProduction);
776 
777  fLBNEOutNtpData->ppenergy = sqrt((parentp*parentp+Parent_mass*Parent_mass))/CLHEP::GeV;
778 
779  fLBNEOutNtpData->ppmedium = 0.; //this is still empty
780 
781  fLBNEOutNtpData->ptype = LBNEParticleCode::AsInt(LBNEParticleCode::StringToEnum(parent_name));
782 
783  G4ThreeVector production_vertex = NuParentTrack->GetPoint(0)->GetPosition();
784  fLBNEOutNtpData->ppvx = production_vertex[0]/CLHEP::cm;
785  fLBNEOutNtpData->ppvy = production_vertex[1]/CLHEP::cm;
786  fLBNEOutNtpData->ppvz = production_vertex[2]/CLHEP::cm;
787 
788  //if nu parent is a muon then find muon parent info
789  if ((parent_name=="mu+" || parent_name=="mu-") && NuParentTrack->GetParentID()!=0)
790  {
791  G4int mupar = NuParentTrack->GetParentID();
792  LBNETrajectory* MuParentTrack = GetParentTrajectory(mupar);
793  G4int nopoint_mupar = MuParentTrack->GetPointEntries();
794  G4ThreeVector muparp = MuParentTrack->GetMomentum(nopoint_mupar-1);
795  G4double muparm = MuParentTrack->GetMass();
796  fLBNEOutNtpData->muparpx = muparp[0]/CLHEP::GeV; // vector of hadron parent of muon
797  fLBNEOutNtpData->muparpy = muparp[1]/CLHEP::GeV; //
798  fLBNEOutNtpData->muparpz = muparp[2]/CLHEP::GeV;
799  fLBNEOutNtpData->mupare = (sqrt(muparp*muparp+muparm*muparm))/CLHEP::GeV;
800  }
801  else
802  {
803  fLBNEOutNtpData->muparpx = -999999.;
804  fLBNEOutNtpData->muparpy = -999999.;
805  fLBNEOutNtpData->muparpz = -999999.;
806  fLBNEOutNtpData->mupare = -999999.;
807  }
808 
809  fLBNEOutNtpData->Necm = enuzr/CLHEP::GeV; // Neutrino energy in parent rest frame
810  LBNETrackInformation* info = (LBNETrackInformation*)(track.GetUserInformation());
811  fLBNEOutNtpData->Nimpwt = info->GetNImpWt(); // Importance weight
812  fLBNEOutNtpData->tgen = info->GetTgen()-1;
813 
814  fLBNEOutNtpData->xpoint = 0.; // x, y, z of parent at user selected vol
815  fLBNEOutNtpData->xpoint = 0.;
816  fLBNEOutNtpData->xpoint = 0.;
817 
818  {
819  G4ThreeVector ParticlePosition = NPGA->GetParticlePosition();
820  fLBNEOutNtpData->tvx = ParticlePosition[0]/CLHEP::cm;
821  fLBNEOutNtpData->tvy = ParticlePosition[1]/CLHEP::cm;
822  fLBNEOutNtpData->tvz = ParticlePosition[2]/CLHEP::cm;
823  G4ThreeVector ParticleMomentum = NPGA->GetParticleMomentum();
824  fLBNEOutNtpData->tpx = ParticleMomentum[0]/CLHEP::GeV;
825  fLBNEOutNtpData->tpy = ParticleMomentum[1]/CLHEP::GeV;
826  fLBNEOutNtpData->tpz = ParticleMomentum[2]/CLHEP::GeV;
827  fLBNEOutNtpData->tptype = NPGA->GetParticleType();
828  }
829  //end find particle exiting target
830 
831  } // Partial sequence of info fill for fLBNEOutNtpData
832 
833  //
834  // Dk2Nu filling
835  //
836  //
837  if (pRunManager->GetCreateDk2NuOutput()) {
838  fDk2Nu->nuray.clear();
839  fDk2Nu->ancestor.clear();
840  fDk2Nu->vint.clear();
841  LBNERunManager* theRunManager = dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
842  fDk2Nu->job = theRunManager->GetCurrentRun()->GetRunID();
843  fDk2Nu->potnum = G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID();
844  LBNETrackInformation* info = (LBNETrackInformation*)(track.GetUserInformation());
845 
846 // bsim::NuRay myNu(NuMomentum[0]/CLHEP::GeV, NuMomentum[1]/CLHEP::GeV, NuMomentum[2]/CLHEP::GeV,
847 // track.GetTotalEnergy()/CLHEP::GeV, info->GetNImpWt());
848  bsim::NuRay myNu(NuMomentum[0]/CLHEP::GeV, NuMomentum[1]/CLHEP::GeV, NuMomentum[2]/CLHEP::GeV,
849  track.GetTotalEnergy()/CLHEP::GeV, 1); // Sept 24, we do not place the importance weight here..
850 // std::cerr << " NuMomentum, Pz in CLHEP::GeV " << NuMomentum[2]/CLHEP::GeV << " in nuray " << myNu.pz << std::endl;
851  fDk2Nu->nuray.push_back(myNu);
852  fDk2Nu->decay.norig=Norig; // Same convention as in G4LBNE old Ntuple, i..e, as above
853  fDk2Nu->decay.ndecay=NuParentTrack->GetDecayCode(); // See details in LBNESteppingAction::CheckDecay
854  fDk2Nu->decay.ntype=track.GetDefinition()->GetPDGEncoding();
855 
856  fDk2Nu->decay.vx = x/CLHEP::cm;
857  fDk2Nu->decay.vy = y/CLHEP::cm;
858  fDk2Nu->decay.vz = z/CLHEP::cm;
859  fDk2Nu->decay.pdpx = ParentMomentumFinal[0]/CLHEP::GeV;
860  fDk2Nu->decay.pdpy = ParentMomentumFinal[1]/CLHEP::GeV;
861  fDk2Nu->decay.pdpz = ParentMomentumFinal[2]/CLHEP::GeV;
862  G4ThreeVector ParentMomentumProduction = NuParentTrack->GetMomentum(0);
863  fDk2Nu->decay.ppdxdz = ParentMomentumProduction[0]/ParentMomentumProduction[2];
864  fDk2Nu->decay.ppdydz = ParentMomentumProduction[1]/ParentMomentumProduction[2];
865  fDk2Nu->decay.pppz = ParentMomentumProduction[2]/CLHEP::GeV;
866  if (std::isnan(fDk2Nu->decay.pdpy) || (std::abs(fDk2Nu->decay.pdpy) < 1.0e-20) || std::isinf(fDk2Nu->decay.pdpy) ||
867  std::isnan(fDk2Nu->decay.pdpx) || (std::abs(fDk2Nu->decay.pdpx) < 1.0e-20) || std::isinf(fDk2Nu->decay.pdpx)){
868  std::cerr << " Anomalous value for the transverse momentum of parent of neutrino " <<
869  " Px " << fDk2Nu->decay.pdpy << "Py " << fDk2Nu->decay.pdpy << " Pz " << fDk2Nu->decay.pdpz << std::endl;
870  }
871  G4double parentp = sqrt(ParentMomentumProduction*ParentMomentumProduction);
872  fDk2Nu->decay.ppenergy = sqrt((parentp*parentp+Parent_mass*Parent_mass))/CLHEP::GeV;
873  fDk2Nu->decay.ppmedium = NuParentTrack->GetMaterialNumber1rst();
874  fDk2Nu->decay.ptype = NuParentTrack->GetPDGEncoding();
875  if ((parent_name=="mu+" || parent_name=="mu-") && NuParentTrack->GetParentID()!=0)
876  {
877  G4int mupar = NuParentTrack->GetParentID();
878  LBNETrajectory* MuParentTrack = GetParentTrajectory(mupar);
879  G4int nopoint_mupar = MuParentTrack->GetPointEntries();
880  G4ThreeVector muparp = MuParentTrack->GetMomentum(nopoint_mupar-1);
881  G4double muparm = MuParentTrack->GetMass();
882  fDk2Nu->decay.muparpx = muparp[0]/CLHEP::GeV; // vector of hadron parent of muon
883  fDk2Nu->decay.muparpy = muparp[1]/CLHEP::GeV; //
884  fDk2Nu->decay.muparpz = muparp[2]/CLHEP::GeV;
885  fDk2Nu->decay.mupare = (sqrt(muparp*muparp+muparm*muparm))/CLHEP::GeV;
886  }
887  else
888  {
889  fDk2Nu->decay.muparpx = -9999999.;
890  fDk2Nu->decay.muparpy = -9999999.;
891  fDk2Nu->decay.muparpz = -9999999.;
892  fDk2Nu->decay.mupare = -9999999.;
893  }
894  fDk2Nu->decay.necm = enuzrInGeV; // Now in CLHEP::GeV...
895  fDk2Nu->decay.nimpwt = info->GetNImpWt(); // duplicate info, to ease analysis. Per Robert Hatcher directive, April 24 2014.
896 
897  //For PPfx stuff..... A. Bashyal
898 
899  //So here we write down the tgtexit ntuples for dk2nu. A. Bashyal 12/22/2015
900  bsim::TgtExit tgt;
901  G4int tar_pdg = 0;
902  int Tptype = 99;
903  //LBNEParticleCode::AsInt(LBNEParticleCode::StringToEnum(PParentTrack->GetParticleName()));
904 
905  //int ptrkid = PParentTrack->GetTrackID();
906  G4ThreeVector TParticleMomentum = G4ThreeVector(-999999,-999999,-999999);
907  G4ThreeVector TParticlePosition = G4ThreeVector(-999999,-999999,-999999);
908 
909  //G4int tar_trackID = -1;
910  //G4bool findtarget = false;
911  //Trying to use an alternative method to find the target exit location for neutrino parents and see if it works:
912  size_t ntraj = history.size();
913  std::vector<G4VTrajectory*>TempHistory(history);
914  std::vector<int>catchercontainer;
915  std::vector<int>trackcontainer;
916  for (size_t iA=0; iA != ntraj; iA++) {
917  //Trying to see if I can inlcude the neutrino info in the trajectory
918  LBNETrajectory *tgtraj = dynamic_cast<LBNETrajectory*>(TempHistory.at(iA));
919  //As of june 20, 2016, it will record any particle exiting the target...this will be handled by ppfx later
920  int Points = tgtraj->GetPointEntries();
921  if((tgtraj->GetMaterialName(0)=="Target")&&(tgtraj->GetMaterialName(Points-1) != "Target")){
922  int transit = Points-1;
923  while (transit != 0){
924  G4String volmat = tgtraj->GetMaterialName(transit);
925  if(volmat.contains("Target")){
926  catchercontainer.push_back(transit);
927  trackcontainer.push_back(iA);
928  // std::cout<<" looking for the loop..."<<std::endl;
929  break;
930 
931  }
932 
933  else --transit;
934  // std::cout<< "transit is "<<transit<<std::endl;
935  }
936  //Now we know the last point at which the parent of the neutrino saw the target.
937  //if(transit<Points-1)std::cout<<" position and momentum "<<tgtraj->GetPoint(transit)->GetPosition()/CLHEP::cm<<std::endl;
938 
939  }
940  }
941  if(trackcontainer.size() != 0){
942  int trackno = trackcontainer.at(0);
943  int pointno = catchercontainer.at(0);
944  LBNETrajectory *tgtraj = dynamic_cast<LBNETrajectory*>(TempHistory.at(trackno));
945  TParticleMomentum = tgtraj->GetMomentum(pointno)/CLHEP::GeV;
946  TParticlePosition = tgtraj->GetPoint(pointno)->GetPosition()/CLHEP::cm;
947  Tptype = LBNEParticleCode::AsInt(LBNEParticleCode::StringToEnum(tgtraj->GetParticleName()));
948  }
949  tgt.tvx = TParticlePosition[0];
950  tgt.tvy = TParticlePosition[1];
951  tgt.tvz = TParticlePosition[2];
952  tgt.tpx = TParticleMomentum[0];
953  tgt.tpy = TParticleMomentum[1];
954  tgt.tpz = TParticleMomentum[2];
955  tgt.tptype = Tptype;
956  // std::cout<<"The tar_pdg is "<<Tptype<<" and momentum is "<<TParticleMomentum[2]/CLHEP::GeV<<" "<<tgt.tvz<<std::endl;
957  //std::cout<<"tgtexit : tvz "<<tgt.tvz<<std::endl;
958  //tgt.tgen = info->GetTgen()-1;
959  fDk2Nu->tgtexit = tgt;
960 
961  // ancestry info.
962  // start by counting the depth of the tree.
963  size_t numAncestors = 0;
964  G4int trackIDTmp = track.GetParentID();
965  LBNETrajectory *tmpTraj = GetTrajectory(trackIDTmp);
966  const LBNEVolumePlacements *aPlacementHandler = LBNEVolumePlacements::Instance();
967  std::vector<G4String> Horn2ICList = aPlacementHandler->GetHorn2ICList();
968  std::vector<G4String>Horn1ICList = aPlacementHandler->GetHorn1ICList();
969  while (trackIDTmp > 0) {
970  numAncestors++;
971  trackIDTmp = tmpTraj->GetParentID();
972  if(trackIDTmp > 0) tmpTraj = GetTrajectory(trackIDTmp);
973  }
974  std::vector<LBNETrajectory *> trajs(numAncestors, 0);
975  // Revert order, starting with the proton...
976  size_t nAnces = numAncestors;
977  trackIDTmp = track.GetParentID();
978  tmpTraj = GetTrajectory(trackIDTmp);
979  while (trackIDTmp > 0) {
980  nAnces--;
981  trajs[nAnces] = tmpTraj;
982  trackIDTmp = tmpTraj->GetParentID();
983  if(trackIDTmp > 0) tmpTraj = GetTrajectory(trackIDTmp);
984  }
985  // Now fill..
986  fDk2Nu->ancestor.clear();
987  fDk2Nu->vint.clear();
988  //A. Bashyal...few more additions for ppfx.
989  size_t nmax = 0;
990  int idx_tar_in_chain = -1;
991  size_t ntrajectory = history.size();
992  std::vector<G4VTrajectory*>TmpHistory(history);
993  std::reverse(TmpHistory.begin(),TmpHistory.end());
994  // for (size_t iA=0; iA != numAncestors; iA++) {
995  //Trying to see if I can inlcude the neutrino info in the trajectory
996  for (size_t iA=0;iA<ntrajectory;++iA){
997  bsim::Ancestor a;
998  //LBNETrajectory *t = trajs[iA];
999  LBNETrajectory *t = dynamic_cast<LBNETrajectory*>(TmpHistory.at(iA));
1000  a.pdg = t->GetPDGEncoding();
1001  //if(t->GetPDGEncoding()==14)std::cout<<"Muon neutrino found "<<std::endl;
1002  G4ThreeVector x1rst = t->GetPoint(0)->GetPosition();
1003  a.startx = x1rst[0]/CLHEP::cm;
1004  a.starty = x1rst[1]/CLHEP::cm;
1005  a.startz = x1rst[2]/CLHEP::cm;
1006  a.startt = t->GetTimeStart(); // ? Units?
1007  G4ThreeVector p1rst = t->GetMomentum(0);
1008  a.startpx = p1rst[0]/CLHEP::GeV;
1009  a.startpy = p1rst[1]/CLHEP::GeV;
1010  a.startpz = p1rst[2]/CLHEP::GeV;
1011  G4int np = t->GetPointEntries();
1012  G4ThreeVector pLast = t->GetMomentum(np-1);
1013  a.stoppx = pLast[0]/CLHEP::GeV;
1014  a.stoppy = pLast[1]/CLHEP::GeV;
1015  a.stoppz = pLast[2]/CLHEP::GeV;
1016  G4ThreeVector pol = t->GetPolarization();
1017  a.polx = pol[0];
1018  a.poly = pol[1];
1019  a.polz = pol[2];
1020  // We do not fill pprodpx, y,z, duplicate info
1021  //a.pprodpx = 0.; a.pprodpy = 0.; a.pprodpz = 0.;
1022  // prodpxyz ntuples are filled as per the requirement in ppfx.....A. Bashyal
1023  //As per the NumiAnalysis, the prodpx saves the momentum of the
1024  //particle before it produces the secondaries. I am assuming this is
1025  //equilavent to the particle momentum before the particle decays into
1026  //whatever it has to.
1027  G4ThreeVector momatProd = t->GetParentMomentumAtThisProduction();
1028  a.pprodpx = momatProd[0]/CLHEP::GeV;
1029  a.pprodpy = momatProd[1]/CLHEP::GeV;
1030  a.pprodpz = momatProd[2]/CLHEP::GeV;
1031  //std::cout<<"pprodpz "<<momatProd/CLHEP::GeV<<std::endl;
1032  //if((t->GetPDGEncoding() != 211) && (momatProd[2]/CLHEP::GeV == 0)std::cout<<" 0 momentum pion "<<t->GetProcessName()<<std::endl;
1033  a.nucleus = t->GetPDGNucleus();
1034  a.proc = t->GetProcessName();
1035  //a.ivol = t->GetVolName1rst();
1036  G4String volname = t->GetVolName1rst();
1037  for(size_t j = 0;j != Horn1ICList.size();j++)
1038  {
1039  if(volname==Horn1ICList.at(j))volname="Horn1IC";
1040  }
1041  for(size_t j = 0;j!=Horn2ICList.size();j++){
1042  if(volname==Horn2ICList.at(j))volname="Horn2IC";
1043  }
1044  //if(volname=="TargetNoSplitSegmentLeft"||volname=="TargetNoSplitSegmentRight")volname="TargetNoSplitSegment";
1045  //below changes are made because of the naming convention in the ppfx-ral
1046  if(volname=="TargetUpstrDownstrSegmentRight"||volname=="TargetUpstrDownstrSegmentLeft"||volname=="RAL-Target"){
1047  volname="TargetNoSplitSegment"; //A. bashyal at leaset upto v3r5c, reference design
1048  }
1049  a.ivol = volname;
1050  // std::cout<<"volname "<<volname<<std::endl;
1051  a.imat = t->GetMaterialName1rst();
1052  //More ppfx stuff.....A. Bashyal
1053  //for(int ap = 0;ap<np-1;ap++){if(t->GetMaterialName(ap).contains("Target"))idx_tar_in_chain = int(iA);}
1054  //Okay I think I did this part wrong. Should only look at the end point of the trajectory to determine this.
1055  if(t->GetMaterialName(np-1).contains("Target"))idx_tar_in_chain = int (iA);
1056  fDk2Nu->ancestor.push_back(a);
1057  }
1058  //
1059  fDk2Nu->vint.push_back(idx_tar_in_chain); //ppfx stuff A. Bashyal...gotta see this one more time..
1060  //Just for the sake of consistency with g4numi I am going to feed in a dummy variable. A. Bashyal
1061  fDk2Nu->vint.push_back(-1); //-1 since we are not shifting the target (see TargetAttenuationReweighter.cpp
1062  //Few stuff on ppfx distance travelled through out volume of interest A. Bashyal
1063  fDk2Nu->vdbl.clear();
1064 
1065 
1066  const int Nanc = 3; //parent, granparent and greatgranparent for ancestry
1067  const double fact_Al = (26.98/2.7);
1068  const double fact_steel316 = (52.73/8.0);
1069  const double fact_concrete = (33.63/2.03);
1070  const double fact_water = (18.01/1.0);
1071  const double fact_iron = (207.19/11.35);
1072  const double fact_helium = (4.003/0.000145);//This from G4numi
1073 
1074  for(int ii = 0;ii<Nanc;ii++){
1075  dist_IC1[ii] = -1*fact_Al;
1076  dist_IC2[ii] = -1*fact_Al;
1077  dist_DPIP[ii] = -1*fact_steel316;
1078  dist_DVOL[ii] = -1*fact_helium;
1079  }
1080 
1081  if(history.size()!=0){
1082 
1083  for(int ii = 0;ii<Nanc;ii++){
1084  dist_IC1[ii] = 0;
1085  dist_IC2[ii] = 0;
1086  dist_DPIP[ii] =0;
1087  dist_DVOL[ii] = 0;
1088  }
1089 
1090  }
1091 
1092  LBNETrajectory* temp_traj;
1093  //we need two different kinds of initializations of these things
1094  //for (size_t ii = 0;ii<std::min(size_t(3),trajs.size());++ii){ //This basically gives parent at 0, greatgranparent at 1 and so on.
1095  for(int ii = 1; ii<4;ii++){
1096  if(history.size()-ii<1)continue;
1097  temp_traj = dynamic_cast<LBNETrajectory*>(TmpHistory.at(ii)); //parsing in the distance travelled info
1098 
1099  double horn1IC = 0;
1100  double horn2IC = 0;
1101  //implementation...yeah!!
1102  for(size_t i = 0;i != Horn1ICList.size();i++)horn1IC += GetDistanceInVolume(temp_traj,Horn1ICList.at(i));
1103  for(size_t i = 0;i != Horn2ICList.size();i++)horn2IC += GetDistanceInVolume(temp_traj,Horn2ICList.at(i));
1104 
1105  //std::cout<<"Total IC1 and IC2 distance "<<horn1IC/fact_Al<<" "<<horn2IC/fact_Al<<std::endl;
1106 
1107  if(horn1IC != 0.)dist_IC1[ii]=horn1IC/fact_Al;
1108 
1109  if(horn2IC != 0.)dist_IC2[ii]=horn2IC/fact_Al;
1110 
1111  dist_DPIP[ii] = GetDistanceInVolume(temp_traj,"DecayPipeWall")/fact_steel316;
1112  dist_DVOL[ii] = GetDistanceInVolume(temp_traj,"DecayPipeVolume")/fact_helium;
1113 
1114  }
1115 
1116  for(int ii=0;ii<3;ii++){
1117  fDk2Nu->vdbl.push_back(dist_IC1[ii]);
1118 
1119  fDk2Nu->vdbl.push_back(dist_IC2[ii]);
1120 
1121  fDk2Nu->vdbl.push_back(dist_DPIP[ii]);
1122 
1123  fDk2Nu->vdbl.push_back(dist_DVOL[ii]);
1124  }
1125 
1126 
1127 
1128 
1129 
1130  //if(dist_DVOL[ii] != 0)std::cout<<dist_DVOL[ii]<<" and level is "<<ii<<std::endl;
1131 
1132 
1133  //Now we write the traj branches....A. Bashyal December 17 2015
1134  bsim::Traj tmp_traj; //Declaration outside the loop??
1135  double trkx[10];
1136  double trky[10];
1137  double trkz[10];
1138  double trkpx[10];
1139  double trkpy[10];
1140  double trkpz[10];
1141 
1142  for (int ii = 0;ii<10;++ii)
1143  {
1144  //bsim::Traj tmp_traj; //Declaration outside the loop??
1145  trkx[ii]=-99999;
1146  trky[ii]=-99999;
1147  trkz[ii]=-99999;
1148  trkpx[ii]=-99999;
1149  trkpy[ii]=-99999;
1150  trkpz[ii]=-99999;
1151  }
1152  G4int npoint = NuParentTrack->GetPointEntries();
1153  //std::cout<<"Number of iterations are :"<<npoint<<std::endl;
1154  G4ThreeVector ParentMomentum;
1155  G4ThreeVector ParentPosition;
1156  // if(h2entryindex>2){std::cout<<" For Particles in horn2 "<<h2entryindex<<" "<<NuParentTrack->GetPreStepVolumeName(h2entryindex)<<
1157  // h2entryindex+1<<" "<<NuParentTrack->GetPreStepVolumeName(h2entryindex+1)<<
1158  // " "<<h2entryindex-1<<NuParentTrack->GetPreStepVolumeName(h2entryindex-1)<<std::endl;}
1159  std::vector<int> h2exitindex;
1160  size_t points = 0;
1161  points = size_t(npoint);
1162  h2exitindex.clear();
1163  for(size_t ii = 0;ii != points-1;ii++){
1164  int iii = int(ii);
1165  G4String h2volname = NuParentTrack->GetPreStepVolumeName(iii);
1166  if(h2volname.contains("Horn2")&&(NuParentTrack->GetMaterialName(iii).contains("Aluminum"))){
1167  h2exitindex.push_back(iii);
1168 
1169 
1170  }
1171  }
1172  int h2exit = -1;
1173  int h2entry = -1;
1174  if(!h2exitindex.empty())h2exit = h2exitindex.back();
1175  if(!h2exitindex.empty())h2entry = h2exitindex.front();
1176  //if(h2exit>2&&h2exit<npoint){std::cout<<" For Particles in horn2 "<<h2exit<<" "<<NuParentTrack->GetPreStepVolumeName(h2exit)<<
1177  // h2exit+1<<" "<<NuParentTrack->GetPreStepVolumeName(h2exit+1)<<
1178  // " "<<h2exit-1<<NuParentTrack->GetPreStepVolumeName(h2exit-1)<<std::endl;}
1179  for(G4int ii = 0; ii<npoint-1;ii++){
1180  ParentMomentum = NuParentTrack->GetMomentum(ii);
1181  ParentPosition = NuParentTrack->GetPoint(ii)->GetPosition();
1182  G4String postvolname = " ";
1183  G4String prevolname = NuParentTrack->GetPreStepVolumeName(ii);
1184  if(ii<npoint-2) postvolname = NuParentTrack->GetPreStepVolumeName(ii+1);
1185 
1186  if((prevolname.contains("TargetNoSplitSegment")||prevolname.contains("TargetFinHorizontal"))&&ii==0){
1187  trkx[0] = ParentPosition[0]/CLHEP::cm;
1188  trky[0] = ParentPosition[1]/CLHEP::cm;
1189  trkz[0] = ParentPosition[2]/CLHEP::cm;
1190  trkpx[0] = ParentMomentum[0]/CLHEP::GeV;
1191  trkpy[0] = ParentMomentum[1]/CLHEP::GeV;
1192  trkpz[0] = ParentMomentum[2]/CLHEP::GeV;
1193 
1194  }
1195  if(prevolname.contains("TargetNoSplitM1")&&postvolname.contains("TargetHallAndHorn1")){
1196  // if(ii==Points){ //I meant to keep it exiting the target but requires more detailed arguments
1197  //std::cout<<"trkx1 "<<ParentPosition[0]/CLHEP::cm<<" from "<<prevolname<<" to "<<postvolname<<std::endl;
1198  trkx[1] = ParentPosition[0]/CLHEP::cm;
1199  trky[1] = ParentPosition[1]/CLHEP::cm;
1200  trkz[1] = ParentPosition[2]/CLHEP::cm;
1201  trkpx[1] = ParentMomentum[0]/CLHEP::GeV;
1202  trkpy[1] = ParentMomentum[1]/CLHEP::GeV;
1203  trkpz[1] = ParentMomentum[2]/CLHEP::GeV;
1204 
1205  }
1206  //if(prevolname.contains("Target"))std::cout<<" Exiting target Prevolname is "<<prevolname<<" and postvol name is "<<
1207  //postvolname<<" at xyz = "<<ParentPosition[0]/CLHEP::cm<<" "<<ParentPosition[1]/CLHEP::cm<<" "<<
1208  //ParentPosition[2]/CLHEP::cm<<std::endl;
1209  if(prevolname.contains("TargetHallAndHorn1") && postvolname.contains("Horn1PolyM1"))
1210  // && !(postvolname.contains("TargetHall")))
1211  {
1212  // std::cout<<"trkx2 "<<ParentPosition[0]/CLHEP::cm<<" from "<<prevolname<<" to "<<postvolname<<std::endl;
1213  trkx[2] = ParentPosition[0]/CLHEP::cm;
1214  trky[2] = ParentPosition[1]/CLHEP::cm;
1215  trkz[2] = ParentPosition[2]/CLHEP::cm;
1216  trkpx[2] = ParentMomentum[0]/CLHEP::GeV;
1217  trkpy[2] = ParentMomentum[1]/CLHEP::GeV;
1218  trkpz[2] = ParentMomentum[2]/CLHEP::GeV;
1219  }
1220  //if(postvolname.contains("Horn2"))std::cout<<" Postvolname contains Horn2 "<< postvolname<< " and prevol name is "<<
1221  // prevolname<<" at z = "<<ParentPosition[2]/CLHEP::cm<<std::endl;
1222  if(prevolname.contains("Horn1PolyM1")&& postvolname.contains("TargetHallAndHorn1"))
1223  {
1224  //std::cout<<"trkx3 "<<ParentPosition[0]/CLHEP::cm<<" from "<<prevolname<<" to "<<postvolname<<" and "<<ParentPosition[2]/CLHEP::cm<<std::endl;
1225  trkx[3] = ParentPosition[0]/CLHEP::cm;
1226  trky[3] = ParentPosition[1]/CLHEP::cm;
1227  trkz[3] = ParentPosition[2]/CLHEP::cm;
1228  trkpx[3] = ParentMomentum[0]/CLHEP::GeV;
1229  trkpy[3] = ParentMomentum[1]/CLHEP::GeV;
1230  trkpz[3] = ParentMomentum[2]/CLHEP::GeV;
1231  }
1232  if(ii==h2entry)
1233  { //need to make sure we dont include the particle exiting horn2 condition.
1234  //std::cout<<"trkx4 "<<ParentPosition[0]/CLHEP::cm<<" from "<<prevolname<<" to "<<postvolname<<std::endl;
1235  trkx[4] = ParentPosition[0]/CLHEP::cm;
1236  trky[4] = ParentPosition[1]/CLHEP::cm;
1237  trkz[4] = ParentPosition[2]/CLHEP::cm;
1238  trkpx[4] = ParentMomentum[0]/CLHEP::GeV;
1239  trkpy[4] = ParentMomentum[1]/CLHEP::GeV;
1240  trkpz[4] = ParentMomentum[2]/CLHEP::GeV;
1241  //std::cout<<" Entering from "<<prevolname<<" Going to "<<postvolname<<" at "<<
1242  //ParentPosition/CLHEP::cm<<std::endl;
1243  }
1244  if(ii==h2exit)
1245  {
1246  //std::cout<<"trkx5 "<<ParentPosition[0]/CLHEP::cm<<" from "<<prevolname<<" to "<<postvolname<<std::endl;
1247  trkx[5] = ParentPosition[0]/CLHEP::cm;
1248  trky[5] = ParentPosition[1]/CLHEP::cm;
1249  trkz[5] = ParentPosition[2]/CLHEP::cm;
1250  trkpx[5] = ParentMomentum[0]/CLHEP::GeV;
1251  trkpy[5] = ParentMomentum[1]/CLHEP::GeV;
1252  trkpz[5] = ParentMomentum[2]/CLHEP::GeV;
1253  //std::cout<<"Entering from"<<prevolname<<" Going to "<<postvolname<<" at "<<ParentPosition/CLHEP::cm<<std::endl;
1254  }
1255  if(prevolname.contains("DecayPipeHall")&&postvolname.contains("DecayPipeVolume"))
1256  {
1257  //std::cout<<"trkx6 "<<ParentPosition[0]/CLHEP::cm<<" from "<<prevolname<<" to "<<postvolname<<std::endl;
1258  trkx[6] = ParentPosition[0]/CLHEP::cm;
1259  trky[6] = ParentPosition[1]/CLHEP::cm;
1260  trkz[6] = ParentPosition[2]/CLHEP::cm;
1261  trkpx[6] = ParentMomentum[0]/CLHEP::GeV;
1262  trkpy[6] = ParentMomentum[1]/CLHEP::GeV;
1263  trkpz[6] = ParentMomentum[2]/CLHEP::GeV;
1264  }
1265 
1266  const double NomH1InnerRad = aPlacementHandler->GetHorn1NeckInnerRadius()/CLHEP::cm;
1267  const double NomH1NeckZloc = aPlacementHandler->GetHorn1NeckZPosition()/CLHEP::cm;
1268  double OptH1InnerRad = 0.0;
1269  double OptH1NeckZLoc = 0.0;
1270  if(aPlacementHandler->GetUseHorn1Polycone()||aPlacementHandler->GetNumberOfHornsPolycone() >0){
1271  OptH1InnerRad =
1272  aPlacementHandler->GetHorn1PolyInnerConductorRadius(1)/CLHEP::cm;
1273  OptH1NeckZLoc = aPlacementHandler->GetHorn1PolyInnerConductorZCoord(3)/CLHEP::cm; //This way the target interactions are not ignored
1274  } //1 for the inner neck radius...see DUNECDR_Optimized.mac in macros.
1275  double H1NeckInnerRad = 0.0;
1276  double H1NeckZLoc = 0.0;
1277  if(NomH1InnerRad != 0){
1278  H1NeckInnerRad = NomH1InnerRad;
1279  H1NeckZLoc = NomH1NeckZloc;
1280  }
1281  else {
1282  H1NeckInnerRad = OptH1InnerRad;
1283  H1NeckZLoc = OptH1NeckZLoc;
1284 
1285 
1286  }
1287  //std::cout<<H1NeckZLoc<< " is neck NLocation"<<std::endl;
1288  //const double H1NeckInnerRad = (aPlacementHandler->GetHorn1NeckInnerRadius())/CLHEP::cm;
1289  const double H2NeckInnerRad = (aPlacementHandler->GetHorn2NeckInnerRadius())/CLHEP::cm;
1290  //std::cout<<"Horn1neck and horn2 neck inner rad are" << H1NeckInnerRad <<" "<<H2NeckInnerRad<<std::endl;
1291 //Exclude the particles through the neck of the horns.
1292  if((sqrt(ParentPosition[0]*ParentPosition[0]+ParentPosition[1]*ParentPosition[1])<H1NeckInnerRad)&&ParentPosition[2]>H1NeckZLoc)
1293  {
1294  trkx[2]=99999;
1295  trky[2]=99999;
1296  trkz[2]=99999;
1297  }
1298  if(sqrt(ParentPosition[0]*ParentPosition[0]+ParentPosition[1]*ParentPosition[1])<H2NeckInnerRad)
1299  {
1300  trkx[4]=99999;
1301  trky[4]=99999;
1302  trkz[4]=99999;
1303  }
1304  }
1305 
1306  ParentMomentum = NuParentTrack->GetMomentum(npoint-1);
1307  ParentPosition = (NuParentTrack->GetPoint(npoint-1)->GetPosition()/CLHEP::m)*CLHEP::m;
1308  trkx[7] = ParentPosition[0]/CLHEP::cm;
1309  trky[7] = ParentPosition[1]/CLHEP::cm;
1310  trkz[7] = ParentPosition[2]/CLHEP::cm;
1311  trkpx[7] = ParentMomentum[0]/CLHEP::GeV;
1312  trkpy[7] = ParentMomentum[1]/CLHEP::GeV;
1313  trkpz[7] = ParentMomentum[2]/CLHEP::GeV;
1314 
1315 
1316  fDk2Nu->traj.clear();
1317 
1318  for(G4int ii = 0;ii<10;++ii)
1319  {
1320  //bsim::Traj tmp_traj;
1321  tmp_traj.trkx = trkx[ii];
1322  tmp_traj.trky = trky[ii];
1323  tmp_traj.trkz = trkz[ii];
1324  tmp_traj.trkpx = trkpx[ii];
1325  tmp_traj.trkpy = trkpy[ii];
1326  tmp_traj.trkpz = trkpz[ii];
1327  fDk2Nu->traj.push_back(tmp_traj);
1328  }
1329 
1330  // Last step : compute the so-called "location weights" for 3 detectors.
1331  //
1332 
1333 
1334  CalcLocationWeights(fDkMeta, fDk2Nu,pRunManager->GetUseRealisticNearDetectorVolume());
1335  //
1336  // Check the weights..
1337  //
1338  for (size_t kLoc = 0; kLoc != fDkMeta->location.size(); kLoc++) {
1339  if ((!std::isnan(fDk2Nu->nuray[kLoc].E)) && (!std::isnan(fDk2Nu->nuray[kLoc].wgt))) continue;
1340  std::cerr << " LBNEAnalysis::FillNeutrinoNtuple Problem with detector weight calculation for location "
1341  << kLoc << std::endl;
1342 
1343  std::cerr << " Decay Vertex : X = " << fDk2Nu->decay.vx << " Y = "
1344  << fDk2Nu->decay.vy << " Z = " << fDk2Nu->decay.vz << std::endl;
1345  std::cerr << " Original track position " << pos.x() << " / " << pos.y() << " / " << pos.z() << std::endl;
1346  std::cerr << " Decay Momentum : " << fDk2Nu->decay.pdpx
1347  << " / " << fDk2Nu->decay.pdpy << " / " << fDk2Nu->decay.pdpy << std::endl;
1348  std::cerr << " Neutrino from G4 decay momentum " << NuMomentum[0]/CLHEP::GeV
1349  << "/ " << NuMomentum[1]/CLHEP::GeV << " / " << NuMomentum[2]/CLHEP::GeV << std::endl;
1350  std::cerr << " decay type " << fDk2Nu->decay.ptype << " parent name " << parent_name <<std::endl;
1351 
1352  }
1353 // std::cerr << " Weight and energy for E [G4] = " << track.GetTotalEnergy()/CLHEP::GeV << " CLHEP::GeV " << std::endl;
1354 // for (size_t kDet = 0; kDet != 3; kDet++ ) {
1355 // std::cerr << " E [CLHEP::GeV] " << fDk2Nu->nuray[kDet].E << " Weight " << fDk2Nu->nuray[kDet].wgt << std::endl;
1356 // }
1357  } // on filling the Dk2nu Ntuple.
1358 
1359  unsigned int ndets = fXdet_near.size();
1360  for(unsigned int idet = 0; idet < ndets; ++idet)
1361  if (pRunManager->GetCreateOutput()) {
1362  fLBNEOutNtpData->NdxdzNear[idet] = (x-fXdet_near[idet])/(z-fZdet_near[idet]);
1363  fLBNEOutNtpData->NdydzNear[idet] = (y-fYdet_near[idet])/(z-fZdet_near[idet]);
1364 
1365  LBNENuWeight nuwgh;
1366  G4double nu_wght;
1367  G4double nu_energy;
1368  std::vector<double> r_det;
1369  r_det.push_back(fXdet_near[idet]/CLHEP::cm);
1370  r_det.push_back(fYdet_near[idet]/CLHEP::cm);
1371  r_det.push_back(fZdet_near[idet]/CLHEP::cm);
1372  nuwgh.GetWeight(fLBNEOutNtpData, r_det,nu_wght,nu_energy);
1373  fLBNEOutNtpData->NenergyN[idet] = nu_energy; //in CLHEP::GeV
1374  fLBNEOutNtpData->NWtNear[idet] = nu_wght;
1375  }
1376  //end near det
1377 
1378  //
1379  //Far Detector
1380 
1381  if(fXdet_far.size() != fYdet_far.size() ||
1382  fXdet_far.size() != fYdet_far.size() ||
1383  fYdet_far.size() != fZdet_far.size())
1384  {
1385  G4cout << "LBNEAnalysis::FillNeutrinoNtuple - "
1386  << "Far Detector vectors are not the same size." << G4endl;
1387  }
1388 
1389 // const int ntypeBefWeight = fLBNEOutNtpData->Ntype; Old statement, lost history... Paul Lebrun, Feb. 2015.
1390 
1391  ndets = fXdet_far.size();
1392  if (pRunManager->GetCreateOutput()) {
1393  for(unsigned int idet = 0; idet < ndets; ++idet)
1394  {
1395  fLBNEOutNtpData->NdxdzFar[idet] = (x-fXdet_far[idet])/(z-fZdet_far[idet]);
1396  fLBNEOutNtpData->NdydzFar[idet] = (y-fYdet_far[idet])/(z-fZdet_far[idet]);
1397 
1398  LBNENuWeight nuwgh;
1399  G4double nu_wght;
1400  G4double nu_energy;
1401  std::vector<double> r_det;
1402  r_det.push_back(fXdet_far[idet]/CLHEP::cm);
1403  r_det.push_back(fYdet_far[idet]/CLHEP::cm);
1404  r_det.push_back(fZdet_far[idet]/CLHEP::cm);
1405  nuwgh.GetWeight(fLBNEOutNtpData, r_det,nu_wght,nu_energy);
1406  fLBNEOutNtpData->NenergyF[idet] = nu_energy; //in CLHEP::GeV
1407  fLBNEOutNtpData->NWtFar[idet] = nu_wght;
1408  }
1409  }
1410 
1411 
1412  //Horn info starts here Amit Bashyal
1413 if(pRunManager->GetCreateOutput()){
1414 
1415  G4int numberOfPoints = NuParentTrack->GetPointEntries();
1416  // std::cout<<NuParentTrack->GetPDGEncoding()<<std::endl;
1417  G4ThreeVector ParticlePosition1 = NuParentTrack->GetPoint(0)->GetPosition();
1418  // std::cout<<ParticlePosition1/CLHEP::cm<<std::endl;
1419  for (G4int ii=numberOfPoints-1; ii > -1; --ii)
1420  {
1421 
1422  G4String prevolname = NuParentTrack->GetPreStepVolumeName(ii);
1423 
1424  G4String postvolname = "";
1425  if(ii < numberOfPoints-1) postvolname = NuParentTrack->GetPreStepVolumeName(ii+1);
1426  //std::cout<<postvolname<<std::endl;
1427  G4ThreeVector ParticleMomentum = NuParentTrack->GetMomentum(ii);
1428  //std::cout<<ParticleMomentum<<" "<<"eventid"<<fLBNEOutNtpData->evtno<<" "<<"trackid "<<NuParentTrack->GetTrackID()<<std::endl;
1429  G4ThreeVector ParticlePosition = NuParentTrack->GetPoint(ii)->GetPosition();
1430  if (((numberOfPoints - ii) < 3) || (ii < 2))
1431  {
1432  //++count;
1433  //if(count > 1)
1434 // std::cout << "eventid: " << fLBNEOutNtpData->evtno << " trackid: " << trackID << " ii = " << ii
1435 // << " particle = " << TrackTrajectory->GetParticleName()
1436 // << " prestepvolname = " << prevolname << " postvolname = " << postvolname
1437 // << " steplength = " << steplength/mm << " mm" << std::endl;
1438  }
1439 
1440  //
1441  //particle created in the target
1442  // For the neutrino ancestors passing through Horn 1 end plane Amit Bashyal
1443  if (prevolname.contains("Horn1TrackingPlane"))
1444  {
1445  fLBNEOutNtpData -> h1posx = ParticlePosition[0]/CLHEP::cm;
1446  fLBNEOutNtpData -> h1posy = ParticlePosition[1]/CLHEP::cm;
1447  fLBNEOutNtpData -> h1posz = ParticlePosition[2]/CLHEP::cm;
1448  fLBNEOutNtpData -> h1momx = ParticleMomentum[0]/CLHEP::GeV;
1449  fLBNEOutNtpData -> h1momy = ParticleMomentum[1]/CLHEP::GeV;
1450  fLBNEOutNtpData -> h1momz = ParticleMomentum[2]/CLHEP::GeV;
1451  fLBNEOutNtpData ->h1trackid = NuParentTrack->GetTrackID();
1452  }
1453 
1454  if(prevolname.contains("Horn2TrackingPlane"))
1455  { //For the Neutrino Ancestors passing through Horn 2 end plane Amit Bashyal
1456  fLBNEOutNtpData -> h2posx = ParticlePosition[0]/CLHEP::cm;
1457  fLBNEOutNtpData -> h2posy = ParticlePosition[1]/CLHEP::cm;
1458  fLBNEOutNtpData -> h2posz = ParticlePosition[2]/CLHEP::cm;
1459  fLBNEOutNtpData -> h2momx = ParticleMomentum[0]/CLHEP::GeV;
1460  fLBNEOutNtpData -> h2momy = ParticleMomentum[1]/CLHEP::GeV;
1461  fLBNEOutNtpData -> h2momz = ParticleMomentum[2]/CLHEP::GeV;
1462  fLBNEOutNtpData ->h2trackid = NuParentTrack->GetTrackID();
1463  }
1464 
1465 
1466  }
1467  }
1468 
1469  G4int trackID = track.GetParentID();
1470  LBNETrajectory* TrackTrajectory = GetTrajectory(trackID);
1471 
1472 
1473  while (trackID > 0)
1474  {
1475  LBNEAnalysis::TrackThroughGeometry(TrackTrajectory);
1476 
1477  trackID = TrackTrajectory->GetParentID();
1478  if(trackID > 0) TrackTrajectory = GetTrajectory(trackID);
1479 
1480 // std::cerr << " Fill Ntuple, Trajectory access for track ancestry " << trackID
1481 // << " From trajectory " << TrackTrajectory->GetParticleDefinition()->GetParticleName()
1482 // << " Parent Track ID " << TrackTrajectory->GetParentID()
1483 // << " num Points " << TrackTrajectory->GetPointEntries() << std::endl;
1484 
1485 
1486  }//end while trackID > 0
1487  //end tracking through geometry
1488 
1489  //
1490  //Done. Fill Tree.
1491  //
1492  if (pRunManager->GetCreateOutput()) fOutTree->Fill();
1493  if (pRunManager->GetCreateDk2NuOutput()) {
1494  fOutTreeDk2Nu->Fill();
1495  }
1496  //
1497  // Write to ascii file
1498  //
1499  if (pRunManager->GetCreateASCIIOutput())
1500  {
1501  std::string asciiFileName = pRunManager->GetOutputASCIIFileName();
1502  std::ofstream asciiFile(asciiFileName.c_str(), std::ios::app);
1503  if(asciiFile.is_open())
1504  {
1505  asciiFile << fLBNEOutNtpData->Ntype<< " " << fLBNEOutNtpData->Nenergy << " "
1506  << fLBNEOutNtpData->NenergyN[0] << " " << fLBNEOutNtpData->NWtNear[0];
1507  asciiFile << " " << fLBNEOutNtpData->NenergyF[0] << " " << fLBNEOutNtpData->NWtFar[0]
1508  <<" "<<fLBNEOutNtpData->Nimpwt<< G4endl;
1509  asciiFile.close();
1510  }
1511  }
1512 
1513 }
1514 
1515 //-----------------------------------------------------------------------------------------
1516 //
1517 void LBNEAnalysis::FillTrackingNtuple(const G4Track& track, LBNETrajectory* currTrajectory)
1518 {
1519 
1520  LBNERunManager* pRunManager =
1521  dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
1522 
1523  if (!pRunManager->GetCreateOutput()) return;
1524 
1525  fLBNEOutNtpData -> Clear();
1526 
1527  LBNEPrimaryGeneratorAction *NPGA =
1528  (LBNEPrimaryGeneratorAction*)(pRunManager)->GetUserPrimaryGeneratorAction();
1529 
1530 
1531 
1532 
1533  fLBNEOutNtpData->protonN = NPGA->GetProtonNumber();
1534 
1535  fLBNEOutNtpData->run = pRunManager->GetCurrentRun()->GetRunID();
1536  fLBNEOutNtpData->evtno = pRunManager->GetCurrentEvent()->GetEventID();
1537  fLBNEOutNtpData->beamHWidth = NPGA->GetBeamSigmaX()/CLHEP::cm;
1538  fLBNEOutNtpData->beamVWidth = NPGA->GetBeamSigmaY()/CLHEP::cm;
1539  fLBNEOutNtpData->beamX = NPGA->GetBeamOffsetX()/CLHEP::cm;
1540  fLBNEOutNtpData->beamY = NPGA->GetBeamOffsetY()/CLHEP::cm;
1541 
1542  //G4int particleID = track.GetParentID();
1543  G4ThreeVector protonOrigin = NPGA->GetProtonOrigin();
1544  fLBNEOutNtpData->protonX = protonOrigin[0];
1545  fLBNEOutNtpData->protonY = protonOrigin[1];
1546  fLBNEOutNtpData->protonZ = protonOrigin[2];
1547 
1548  G4ThreeVector protonMomentum = NPGA->GetProtonMomentum();
1549  fLBNEOutNtpData->protonPx = protonMomentum[0];
1550  fLBNEOutNtpData->protonPy = protonMomentum[1];
1551  fLBNEOutNtpData->protonPz = protonMomentum[2];
1552 
1553  LBNEVolumePlacements *volDB = LBNEVolumePlacements::Instance();
1554 
1555  fLBNEOutNtpData->nuTarZ = volDB->GetTargetLengthIntoHorn(); // A better info that the somewhat meaningless Z0
1556  const LBNEDetectorConstruction *detDB =
1557  dynamic_cast<const LBNEDetectorConstruction*>(pRunManager->GetUserDetectorConstruction());
1558  fLBNEOutNtpData->hornCurrent = detDB->GetHornCurrent()/CLHEP::ampere/1000.;
1559 
1560  G4int trackID = track.GetTrackID();
1561  G4String currentTrackVolName = track.GetVolume() -> GetName();
1562  LBNETrajectory* TrackTrajectory = currTrajectory;
1563 
1564  //G4int eventID = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
1565  //std::cout << "Event: " << eventID << " TrackID: " << trackID << " curr vol = " << currentTrackVolName << std::endl;
1566 
1567 
1568  while (trackID > 0)
1569  {
1570 
1571  LBNEAnalysis::TrackThroughGeometry(TrackTrajectory);
1572 
1573  trackID = TrackTrajectory->GetParentID();
1574  if(trackID > 0) TrackTrajectory = GetTrajectory(trackID);
1575 
1576  }
1577 
1578  fOutTree->Fill();
1579 
1580 }
1581 
1582 //-----------------------------------------------------------------------------------------
1583 void LBNEAnalysis::TrackThroughGeometry(const LBNETrajectory* TrackTrajectory)
1584 {
1585  G4int trackID = TrackTrajectory-> GetTrackID();
1586 
1587  G4int numberOfPoints = TrackTrajectory->GetPointEntries();
1588 
1589  for (G4int ii=numberOfPoints-1; ii > -1; --ii)
1590  {
1591 
1592  G4String prevolname = TrackTrajectory->GetPreStepVolumeName(ii);
1593  G4String postvolname = "";
1594  if(ii < numberOfPoints-1) postvolname = TrackTrajectory->GetPreStepVolumeName(ii+1);
1595 
1596  G4ThreeVector ParticleMomentum = TrackTrajectory->GetMomentum(ii);
1597  G4ThreeVector ParticlePosition = TrackTrajectory->GetPoint(ii)->GetPosition();
1598 
1599  if (((numberOfPoints - ii) < 3) || (ii < 2))
1600  {
1601  //++count;
1602 // G4double steplength = TrackTrajectory->GetStepLength(ii);
1603  //if(count > 1)
1604 // std::cout << "eventid: " << fLBNEOutNtpData->evtno << " trackid: " << trackID << " ii = " << ii
1605 // << " particle = " << TrackTrajectory->GetParticleName()
1606 // << " prestepvolname = " << prevolname << " postvolname = " << postvolname
1607 // << " steplength = " << steplength/mm << " mm" << std::endl;
1608  }
1609 
1610  TrackPoint_t TrkPt;
1611  TrkPt.type = LBNEParticleCode::AsInt(LBNEParticleCode::StringToEnum(TrackTrajectory->GetParticleName()));
1612  TrkPt.x = ParticlePosition[0]/CLHEP::cm;
1613  TrkPt.y = ParticlePosition[1]/CLHEP::cm;
1614  TrkPt.z = ParticlePosition[2]/CLHEP::cm;
1615  TrkPt.px = ParticleMomentum[0]/CLHEP::GeV;
1616  TrkPt.py = ParticleMomentum[1]/CLHEP::GeV;
1617  TrkPt.pz = ParticleMomentum[2]/CLHEP::GeV;
1618  TrkPt.trkid = trackID;
1619  TrkPt.impwt = TrackTrajectory->GetNImpWt();
1620  TrkPt.gen = (TrackTrajectory->GetTgen())-1;
1621 
1622  //
1623  //particle created in the target
1624  //
1625  if (prevolname.contains("TargetFin") && ii==0)
1626  {
1627  fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Target"), TrkPt);
1628  }
1629  //
1630  //particle as exits target
1631  //
1632  if ((prevolname.contains("Target")) &&
1633  postvolname.contains("Horn1") && (!postvolname.contains("Target")))
1634  {
1635  /*if(fLBNEOutNtpData->evtno == 1450 || fLBNEOutNtpData->evtno == 1679)
1636  {
1637  std::cout << "EVENT ID: " << fLBNEOutNtpData->evtno << std::endl;
1638  std::cout << " LBNEAnalysis particle exiting target..." << std::endl;
1639  std::cout << " track id = " << TrkPt.trkid << std::endl;
1640  std::cout << " type = " << TrkPt.type << std::endl;
1641  std::cout << " x pos = " << TrkPt.x << std::endl;
1642  std::cout << " y pos = " << TrkPt.y << std::endl;
1643  std::cout << " z pos = " << TrkPt.z << std::endl;
1644  std::cout << " px mom = " << TrkPt.px << std::endl;
1645  std::cout << " py mom = " << TrkPt.py << std::endl;
1646  std::cout << " pz mom = " << TrkPt.pz << std::endl;
1647  std::cout << " LBNEAnalysis Nu Parent found exiting target..." << std::endl;
1648  std::cout << " track id = " << tptrkid << std::endl;
1649  std::cout << " type = " << fLBNEOutNtpData->tptype << std::endl;
1650  std::cout << " x pos = " << fLBNEOutNtpData->tvx << std::endl;
1651  std::cout << " y pos = " << fLBNEOutNtpData->tvy << std::endl;
1652  std::cout << " z pos = " << fLBNEOutNtpData->tvz << std::endl;
1653  std::cout << " px mom = " << fLBNEOutNtpData->tpx << std::endl;
1654  std::cout << " py mom = " << fLBNEOutNtpData->tpy << std::endl;
1655  std::cout << " pz mom = " << fLBNEOutNtpData->tpz << std::endl;
1656  }*/
1657 
1658  fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("TargetExit"), TrkPt);
1659  }
1660  //
1661  //particle at plane of end of tgt. This only makes sense
1662  //if only running with the target and target hall constructed.
1663  //otherwise this will
1664  //
1665  if (prevolname.contains("TgtEndPlane"))
1666  {
1667  fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("TargetEndPlane"), TrkPt);
1668  }
1669 
1670 
1671 // if(LBNEData->GetSimulation() == "Target Tracking") continue;
1672  //
1673  // The following volume names need updating.....
1674  //
1675 
1676  //particle enters horn 1
1677  if ((prevolname.contains("Horn1TopLevelDownstr")) && postvolname.contains("Horn1TopLevelDownstr"))
1678  {
1679  fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Horn1Enter"), TrkPt);
1680  }
1681  //particle at neck of horn 1 the last entry will be just before it leaves the neck
1682  //if (prevolname.contains("PH01-02"))
1683  if (prevolname.contains("Horn1DownstrPart1"))
1684  {
1685  fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Horn1NeckPlane"), TrkPt);
1686  }
1687  //particle exists horn 1
1688  if ((prevolname.contains("Horn1")) && postvolname.contains("Tunnel"))
1689  {
1690  fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Horn1Exit"), TrkPt);
1691  }
1692  //particle at end plane of horn 1
1693  if (prevolname.contains("PH01EndPlane"))
1694  {
1695  fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Horn1EndPlane"), TrkPt);
1696  }
1697  //particle enters horn 2
1698  if ((prevolname.contains("Tunnel")) && postvolname.contains("Horn2"))
1699  {
1700  fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Horn2Enter"), TrkPt);
1701  }
1702  //particle at neck of horn 2 the last entry will be just before it leaves the neck
1703  if (prevolname.contains("Horn2Part2"))
1704  {
1705  fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Horn2NeckPlane"), TrkPt);
1706  }
1707  //particle exits horn 2
1708  if ((prevolname.contains("Horn2")) && postvolname.contains("Tunnel"))
1709  {
1710  fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Horn2Exit"), TrkPt);
1711  }
1712  //particle at end plane of horn 2
1713  if (prevolname.contains("PH02EndPlane"))
1714  {
1715  fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Horn2EndPlane"), TrkPt);
1716  }
1717  //Particle enters DP
1718  if ( (prevolname.contains("DecayPipeUsptrWindow")) &&
1719  (postvolname.contains("DecayPipeVolume"))) {
1720 
1721  fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("DPipeEnter"), TrkPt);
1722  }
1723  //Particle exits DP
1724  if ( (prevolname.contains("DecayPipe")) && (postvolname.contains("Tunnel")))
1725 
1726  {
1727  fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("DPipeExit"), TrkPt);
1728  }
1729 
1730  }//end loop over npts in Track
1731 
1732 }
1733 
1734 
1735 //------------------------------------------------------------------------------------
1736 LBNETrajectory* LBNEAnalysis::GetParentTrajectory(G4int parentID)
1737 {
1738 
1739  LBNERunManager *pRunManager=
1740  dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
1741 
1742  if (pRunManager->GetVerboseLevel() == 10) {
1743  G4cout << "LBNEAnalysis::GetParentTrajectory() called." << G4endl;
1744  }
1745 
1746  G4TrajectoryContainer* container =
1747  G4RunManager::GetRunManager()->GetCurrentEvent()->GetTrajectoryContainer();
1748  if(container==0) return 0;
1749 
1750  TrajectoryVector* vect = container->GetVector();
1751  G4VTrajectory* tr;
1752  G4int ii = 0;
1753  while (ii<G4int(vect->size())){
1754  tr = (*vect)[ii];
1755  LBNETrajectory* tr1 = (LBNETrajectory*)(tr);
1756  if(tr1->GetTrackID() == parentID) return tr1;
1757  ++ii;
1758  }
1759 
1760  return 0;
1761 }
1762 
1763 //----------------------------------------------------------------------------------------------------------
1764 LBNETrajectory* LBNEAnalysis::GetTrajectory(G4int trackID)
1765 {
1766  G4TrajectoryContainer* container = G4RunManager::GetRunManager()->GetCurrentEvent()->GetTrajectoryContainer();
1767 
1768  if(container==0)
1769  {
1770  G4cout << "LBNEAnalysis::GetTrajectory - PROBLEM: No Trajectory Container for track ID = " << trackID << endl;
1771  return 0;
1772  }
1773 
1774  TrajectoryVector* vect = container->GetVector();
1775  G4VTrajectory* tr;
1776  G4int ii = 0;
1777  while (ii < G4int(vect->size()))
1778  {
1779  tr = (*vect)[ii];
1780  LBNETrajectory* tr1 = (LBNETrajectory*)(tr);
1781  if(tr1->GetTrackID() == trackID) return tr1;
1782  ++ii;
1783  }
1784 
1785  G4cout << "LBNEAnalysis::GetTrajectory - PROBLEM: Failed to find track with ID = " << trackID << endl;
1786 
1787  return 0;
1788 }
1789 void LBNEAnalysis::setDetectorPositions() {
1790 
1791  // Info passed along in the Ntuple files, old or new (Dk2Nu), so we place these data statements in one place.
1792 
1793  const int nNear=5;
1794  const int nFar=3;
1795  G4double xdetNear[] = { 0. , 0. , 11.50, 0. , 25.84 };
1796  G4double ydetNear[] = { 0. , 0. , -3.08, 0. , 78.42 };
1797  G4double zdetNear[] = {574. ,1036.49 , 1000.97, 1030.99 , 745.25 };
1798  fDetNameNear.clear();
1799  fDetNameNear.push_back(std::string("LBNE"));
1800 // fDetNameNear.push_back(std::string("Minos"));
1801 // fDetNameNear.push_back(std::string("Nova"));
1802 // fDetNameNear.push_back(std::string("Minerva"));
1803 // fDetNameNear.push_back(std::string("MiniBooNE"));
1804  G4double xdetFar[] = { 0. , 0. , 11040. };
1805  G4double ydetFar[] = { 0. , 0. , -4200. };
1806  G4double zdetFar[] = { 1297000. , 735340. , 810450. };
1807  fDetNameFar.clear();
1808  fDetNameFar.push_back(std::string("LBNE"));
1809 // fDetNameFar.push_back(std::string("Minos"));
1810 // fDetNameFar.push_back(std::string("Nova"));
1811 
1812 // get(xdetNear[0], "FluxAreaX0near");
1813 // get(ydetNear[0], "FluxAreaY0near");
1814 // get(zdetNear[0], "FluxAreaZ0near");
1815 
1816 // get(xdetFar[0], "FluxAreaX0far");
1817 // get(ydetFar[0], "FluxAreaY0far");
1818 // get(zdetFar[0], "FluxAreaZ0far");
1819  fXdet_near.clear();
1820  fYdet_near.clear();
1821  fZdet_near.clear();
1822 //
1823  fXdet_near.resize(nNear);
1824  fYdet_near.resize(nNear);
1825  fZdet_near.resize(nNear);
1826 
1827  for(G4int ii=0;ii<nNear;ii++){
1828  fXdet_near[ii] = xdetNear[ii]*CLHEP::m;
1829  fYdet_near[ii] = ydetNear[ii]*CLHEP::m;
1830  fZdet_near[ii] = zdetNear[ii]*CLHEP::m;
1831  }
1832  //Far Detector
1833  // Neutrino data at different points
1834  // need neutrino parent info to be filled in fLBNEOutNtpData by this point
1835  //
1836  fXdet_far.clear();
1837  fYdet_far.clear();
1838  fZdet_far.clear();
1839  fXdet_far.resize(nFar);
1840  fYdet_far.resize(nFar);
1841  fZdet_far.resize(nFar);
1842 
1843  for(G4int ii=0;ii<nFar;ii++){
1844  fXdet_far[ii] = xdetFar[ii]*CLHEP::m;
1845  fYdet_far[ii] = ydetFar[ii]*CLHEP::m;
1846  fZdet_far[ii] = zdetFar[ii]*CLHEP::m;
1847  }
1848 
1849 
1850 }
1851 //----------------------------------------------------------------------------
1852 //---- for tracking planes
1853 void LBNEAnalysis::FillTrackingPlaneData(const G4Step& step){
1854 
1855  LBNERunManager *pRunManager=
1856  dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
1857 
1858  if (pRunManager->GetCreateTrkPlaneOutput()){
1859 
1860  if(fNParticles < kMaxP){
1861  G4Track *track = step.GetTrack();
1862  LBNETrackInformation *info =
1863  (LBNETrackInformation*)(track->GetUserInformation());
1864  if(info != 0){
1865  fNImpWt[fNParticles] = info->GetNImpWt();
1866  }
1867  G4StepPoint *point = step.GetPostStepPoint();
1868  G4ThreeVector pos = point->GetPosition();
1869  G4ThreeVector dir = point->GetMomentumDirection();
1870  G4ParticleDefinition *def = track->GetDefinition();
1871  fParticleX[fNParticles] = pos.x()/CLHEP::m;
1872  fParticleY[fNParticles] = pos.y()/CLHEP::m;
1873  fParticleZ[fNParticles] = pos.z()/CLHEP::m;
1874  fParticleDX[fNParticles] = dir.x();
1875  fParticleDY[fNParticles] = dir.y();
1876  fParticleDZ[fNParticles] = dir.z();
1877  fParticlePDG[fNParticles] = def->GetPDGEncoding();
1878  fParticleMass[fNParticles] = def->GetPDGMass();
1879  fParticleEnergy[fNParticles] = point->GetKineticEnergy()/CLHEP::GeV;
1880  fTrackID[fNParticles] = track->GetTrackID();
1881  fNParticles++;
1882  }
1883  // std::cout << "I stored muon data" << std::endl;
1884 
1885  }// end if track planes is on
1886 
1887 }
1888 
1889 // The ntuples for the Horn 1 end plane are generated here. Amit Bashyal
1890 void LBNEAnalysis::FillTrackingPlaneH1Data(const G4Step& step){
1891 
1892  LBNERunManager *pRunManager=
1893  dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
1894 
1895  if (pRunManager->GetCreateTrkPlaneH1Output()){
1896 
1897  G4Track *track = step.GetTrack();
1898 
1899 
1900  LBNETrackInformation *info =
1901  (LBNETrackInformation*)(track->GetUserInformation());
1902  if(info != 0){
1903  fH1NImpWt = info->GetNImpWt();
1904  }
1905  G4ParticleDefinition *def = track->GetDefinition();
1906  fH1ParticlePDG = def->GetPDGEncoding();
1907 
1908  G4StepPoint *point = step.GetPostStepPoint();
1909  G4ThreeVector pos = point->GetPosition();
1910  G4ThreeVector mom = point->GetMomentum();
1911 
1912  const G4ThreeVector ppoint = track->GetVertexPosition();
1913  const G4ThreeVector pmom = track->GetVertexMomentumDirection();
1914  fH1PProductionDX = pmom.x();
1915  fH1PProductionDY = pmom.y();
1916  fH1PProductionDZ = pmom.z();
1917  fH1PProductionX = ppoint.x()/CLHEP::cm;
1918  fH1PProductionY = ppoint.y()/CLHEP::cm;
1919  fH1PProductionZ = ppoint.z()/CLHEP::cm;
1920  fH1ParticlePX = mom.x()/CLHEP::GeV;
1921  fH1ParticlePY = mom.y()/CLHEP::GeV;
1922  fH1ParticlePZ = mom.z()/CLHEP::GeV;
1923  fH1ParticlePXPZ = mom[0]/mom[2];
1924  fH1ParticlePYPZ = mom[1]/mom[2];
1925  G4ThreeVector dir = point->GetMomentumDirection();
1926  fH1ParticleX = pos.x()/CLHEP::cm;
1927  fH1ParticleY = pos.y()/CLHEP::cm;
1928  fH1ParticleZ = pos.z()/CLHEP::cm;
1929  fH1ParticleDX = dir.x();
1930  fH1ParticleDY = dir.y();
1931  fH1ParticleDZ = dir.z();
1932  fH1ParticleMass = def->GetPDGMass();
1933  fH1ParticleEnergy = point->GetTotalEnergy()/CLHEP::GeV;
1934  fH1TrackID = track->GetTrackID();
1935 
1936 
1937 
1938  fH1NParticles++;
1939  if (fHorn1TrackingTree != 0) fHorn1TrackingTree->Fill();
1940  // }
1941  // std::cout << "I stored muon data" << std::endl;
1942 
1943  }// end if track planes is on
1944 
1945 }
1946 //The ntuples for Horn2 End plane are generated here. Amit Bashyal
1947 void LBNEAnalysis::FillTrackingPlaneH2Data(const G4Step& step){
1948 
1949  LBNERunManager *pRunManager=
1950  dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
1951 
1952  if (pRunManager->GetCreateTrkPlaneH2Output()){
1953 
1954 // if(fH2NParticles < HkMaxP){
1955  G4Track *track = step.GetTrack();
1956  LBNETrackInformation *info =
1957  (LBNETrackInformation*)(track->GetUserInformation());
1958  if(info != 0){
1959  fH2NImpWt = info->GetNImpWt();
1960  }
1961 
1962  if(track->GetCreatorProcess()!=0){
1963 
1964  fH2CProcess = track->GetCreatorProcess()->GetProcessName();
1965  if( fH2CProcess.contains("alphaInelastic"))fH2CNProcess = 1;
1966  else if(fH2CProcess.contains("AntiLambdaInelastic"))fH2CNProcess = 2;
1967  else if(fH2CProcess.contains("AntiNeutronInelastic"))fH2CNProcess = 3;
1968  else if(fH2CProcess.contains("AntiProtonInelastic"))fH2CNProcess = 4;
1969  else if(fH2CProcess.contains("AntiSigmaMinusInelastic"))fH2CNProcess = 5;
1970  else if(fH2CProcess.contains("AntiSigmaPlusInelastic"))fH2CNProcess = 6;
1971  else if(fH2CProcess.contains("AntiXiMinusInelastic"))fH2CNProcess = 7;
1972  else if(fH2CProcess.contains("Decay"))fH2CNProcess = 8;
1973  else if(fH2CProcess.contains("hadElastic"))fH2CNProcess = 9;
1974  else if(fH2CProcess.contains("KaonMinusInelastic"))fH2CNProcess = 10;
1975  else if(fH2CProcess.contains("KaonZeroLInelastic"))fH2CNProcess = 11;
1976  else if(fH2CProcess.contains("KaonZeroSInelastic"))fH2CNProcess = 12;
1977  else if(fH2CProcess.contains("LambdaInelastic"))fH2CNProcess = 13;
1978  else if(fH2CProcess.contains("NeutronInelastic"))fH2CNProcess = 14;
1979  else if(fH2CProcess.contains("PionMinusInelastic"))fH2CNProcess = 15;
1980  else if(fH2CProcess.contains("PionPlusInelastic"))fH2CNProcess = 16;
1981  else if(fH2CProcess.contains("ProtonInelastic"))fH2CNProcess = 17;
1982  else if(fH2CProcess.contains("SigmaMinusInelastic"))fH2CNProcess = 18;
1983  else if(fH2CProcess.contains("SigmaPlusInelastic"))fH2CNProcess = 19;
1984  else if(fH2CProcess.contains("XiMinusInelastic"))fH2CNProcess = 20;
1985  else if(fH2CProcess.contains("XiZeroInelastic"))fH2CNProcess = 21;
1986  else if(fH2CProcess.contains("tInelastic"))fH2CNProcess = 22;
1987 
1988  else fH2CNProcess = 25;
1989 
1990  }
1991  else{
1992  fH2CProcess = "NA";
1993  fH2CNProcess = 0;
1994  }
1995 
1996 
1997  const G4ThreeVector ppoint = track->GetVertexPosition();
1998  const G4ThreeVector pmom = track->GetVertexMomentumDirection();
1999  fH2PProductionDX = pmom.x();
2000  fH2PProductionDY = pmom.y();
2001  fH2PProductionDZ = pmom.z();
2002  G4StepPoint *point = step.GetPostStepPoint();
2003  G4ThreeVector pos = point->GetPosition();
2004  G4ThreeVector dir = point->GetMomentumDirection();
2005  G4ParticleDefinition *def = track->GetDefinition();
2006  G4ThreeVector mom = point->GetMomentum();
2007  fH2PProductionX = ppoint.x()/CLHEP::cm;
2008  fH2PProductionY = ppoint.y()/CLHEP::cm;
2009  fH2PProductionZ = ppoint.z()/CLHEP::cm;
2010  fH2ParticlePX = mom.x()/CLHEP::GeV;
2011  fH2ParticlePY = mom.y()/CLHEP::GeV;
2012  fH2ParticlePZ = mom.z()/CLHEP::GeV;
2013  fH2ParticlePXPZ = mom[0]/mom[2];
2014  fH2ParticlePYPZ = mom[1]/mom[2];
2015  fH2ParticleX = pos.x()/CLHEP::cm;
2016  fH2ParticleY = pos.y()/CLHEP::cm;
2017  fH2ParticleZ = pos.z()/CLHEP::cm;
2018  fH2ParticleDX = dir.x();
2019  fH2ParticleDY = dir.y();
2020  fH2ParticleDZ = dir.z();
2021  fH2ParticlePDG = def->GetPDGEncoding();
2022  fH2ParticleMass = def->GetPDGMass();
2023  fH2ParticleEnergy = point->GetTotalEnergy()/CLHEP::GeV;
2024  fH2TrackID = track->GetTrackID();
2025 
2026  // std::cout<<fH2CProcess<<" "<<fH2ParticlePDG<<" "<<ppoint<<" "<<mom<<std::endl;
2027 
2028 
2029  fH2NParticles++;
2030  if (fHorn2TrackingTree != 0) fHorn2TrackingTree->Fill();
2031 
2032  //}
2033 
2034  }// end if track planes is on
2035 
2036 }
2037 
2038 //The ntuples for Decay Pipe End plane are created here. Amit Bashyal
2039 void LBNEAnalysis::FillTrackingPlaneDPData(const G4Step& step){
2040 
2041  LBNERunManager *pRunManager=
2042  dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
2043 
2044  if (pRunManager->GetCreateTrkPlaneDPOutput()){
2045 
2046  G4Track *track = step.GetTrack();
2047  // LBNETrajectory *ptrack = (LBNETrajectory*)(track->GetTrajectory());
2048  LBNETrackInformation *info =
2049  (LBNETrackInformation*)(track->GetUserInformation());
2050  if(info != 0){
2051  fDPNImpWt = info->GetNImpWt();
2052  }
2053  G4StepPoint *point = step.GetPostStepPoint();
2054  G4ThreeVector pos = point->GetPosition();
2055  G4ThreeVector mom = point->GetMomentum();
2056  const G4ThreeVector ppoint = track->GetVertexPosition();
2057  const G4ThreeVector pmom = track->GetVertexMomentumDirection();
2058  fDPPProductionDX = pmom.x();
2059  fDPPProductionDY = pmom.y();
2060  fDPPProductionDZ = pmom.z();
2061  fDPPProductionX = ppoint.x()/CLHEP::cm;
2062  fDPPProductionY = ppoint.y()/CLHEP::cm;
2063  fDPPProductionZ = ppoint.z()/CLHEP::cm;
2064  fDPParticlePX = mom.x()/CLHEP::GeV;
2065  fDPParticlePY = mom.y()/CLHEP::GeV;
2066  fDPParticlePZ = mom.z()/CLHEP::GeV;
2067  fDPParticlePXPZ = mom[0]/mom[2];
2068  fDPParticlePYPZ = mom[1]/mom[2];
2069  G4ThreeVector dir = point->GetMomentumDirection();
2070  G4ParticleDefinition *def = track->GetDefinition();
2071  fDPParticleX = pos.x()/CLHEP::cm;
2072  fDPParticleY = pos.y()/CLHEP::cm;
2073  fDPParticleZ = pos.z()/CLHEP::cm;
2074  fDPParticleDX = dir.x();
2075  fDPParticleDY = dir.y();
2076  fDPParticleDZ = dir.z();
2077  fDPParticlePDG = def->GetPDGEncoding();
2078  fDPParticleMass = def->GetPDGMass();
2079  fDPParticleEnergy = point->GetTotalEnergy()/CLHEP::GeV;
2080  fDPTrackID = track->GetTrackID();
2081  fDPParentID = track->GetParentID();
2082  fDPNParticles++;
2083  if (fDPTrackingTree!= 0) fDPTrackingTree->Fill();
2084 
2085  }// end if track planes is on
2086 
2087 }
2088 // I tried to write the ntuples for the particles coming out of the target without implementing any target plane. However, because the LBNE target is made up of
2089 // small segments, the info of same particle is stored multiple times as it enters from one segment of the target to the another segment.
2090 // At this point, it is unreliable. Amit Bashyal
2091 void LBNEAnalysis::FillTargetOutputData(const G4Step& step){
2092  LBNERunManager *pRunManager=
2093  dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
2094  if (pRunManager->GetCreateTargetOutput()){
2095 
2096  if(fTNParticles < HkMaxP){
2097  G4Track *track = step.GetTrack();
2098  // We are ignoring the Importance weight here... Why.. ? Paul L. G. Lebrun.
2099  // LBNETrajectory *ptrack = (LBNETrajectory*)(track->GetTrajectory());
2100 // LBNETrackInformation *info =
2101 // (LBNETrackInformation*)(track->GetUserInformation());
2102 // if(info != 0){
2103  // fTNImpWt = info->GetNImpWt();
2104  // }
2105  G4String prematerial = step.GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial()->GetName();
2106  G4String posmaterial = step.GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial()->GetName();
2107  G4StepPoint *point = step.GetPostStepPoint();
2108  G4ThreeVector pos = point->GetPosition();
2109  G4ThreeVector mom = point->GetMomentum();
2110  fTParticlePX = mom.x()/CLHEP::GeV;
2111  fTParticlePY = mom.y()/CLHEP::GeV;
2112  fTParticlePZ = mom.z()/CLHEP::GeV;
2113  fTParticleX = pos.x()/CLHEP::cm;
2114  fTParticleY = pos.y()/CLHEP::cm;
2115  fTParticleZ = pos.z()/CLHEP::cm;
2116  G4ParticleDefinition *def = track->GetDefinition();
2117  fTParticlePDG = def->GetPDGEncoding();
2118  fTParticleEnergy = point->GetTotalEnergy()/CLHEP::GeV;
2119  fTTrackID = track->GetTrackID();
2120  fTParentID = track->GetParentID();
2121  fTNParticles++;
2122  }
2123  // fTargetOutputTree->Fill();
2124 
2125 
2126  }// end if track planes is on
2127 
2128 }
2129 
2130 void LBNEAnalysis::FillAlcoveTrackingPlaneData(const G4Step& aStep)
2131 {
2132 
2133  LBNERunManager *pRunManager=
2134  dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
2135  if (pRunManager->GetCreateAlcoveTrackingOutput())
2136  {
2137  G4Track *track = aStep.GetTrack();
2138  LBNETrackInformation *info =
2139  (LBNETrackInformation*)(track->GetUserInformation());
2140  G4StepPoint *point = aStep.GetPostStepPoint();
2141  G4ThreeVector mom = point->GetMomentum();
2142  G4ThreeVector pos = point->GetPosition();
2143 
2144  int theTrackID = aStep.GetTrack()->GetTrackID();
2145  int index = -1;
2146 
2147 
2148  //First look for this track
2149  for (int i = 0 ; i < fMuNParticles; i++)
2150  {
2151  if (fMuTrackID[i] == theTrackID){
2152  index = i;
2153  }
2154  }
2155 
2156  if (index > -1){//Found track
2157 
2158  fMuNSteps[index]++;
2159  fMuDE[index] += aStep.GetTotalEnergyDeposit()/CLHEP::MeV;
2160  fMuDEion[index] += (aStep.GetTotalEnergyDeposit()- aStep.GetNonIonizingEnergyDeposit())/CLHEP::MeV;
2161  fMudEdx[index] += aStep.GetStepLength()/CLHEP::cm;
2162  fMudEdx_ion[index] += aStep.GetStepLength()/CLHEP::cm;
2163  }else{//Track not already saved
2164 
2165  index = fMuNParticles;
2166  if (index >= kMaxP)
2167  {
2168  G4cout << "Alcove Tracking: Max Number of Particles Reached. Not saving ... " <<G4endl;
2169  return;
2170  }
2171  fMuNParticles++;
2172  fMuEvtNo = pRunManager->GetCurrentEvent()->GetEventID();
2173  fMuRunNo = pRunManager->GetCurrentRun()->GetRunID();
2174 
2175  G4ThreeVector vpos = track->GetVertexPosition();
2176  if(info != 0){
2177  fMuNimpWt[index] = info->GetNImpWt();
2178  }
2179 
2180  G4ParticleDefinition* theParticle = track->GetDefinition();
2181 
2182  fMuTrackID[index] = theTrackID;
2183  fMuParentID[index] = track->GetParentID();
2184 
2185  fMuPDG[index] = theParticle->GetPDGEncoding();
2186  fMuMass[index] = theParticle->GetPDGMass()/CLHEP::GeV;
2187 
2188  fMuT[index] = point->GetGlobalTime()/CLHEP::ns;
2189  fMuPX[index] = mom.x()/CLHEP::GeV;
2190  fMuPY[index] = mom.y()/CLHEP::GeV;
2191  fMuPZ[index] = mom.z()/CLHEP::GeV;
2192  fMuX[index] = pos.x()/CLHEP::cm;
2193  fMuY[index] = pos.y()/CLHEP::cm;
2194  fMuZ[index] = pos.z()/CLHEP::cm;
2195  fMuTheta[index] = mom.z()/sqrt(mom.z()*mom.z()+mom.y()*mom.y()+mom.x()*mom.x());
2196  fMuEnergy[index] = point->GetKineticEnergy()/CLHEP::GeV;
2197  fMuStartX[index] = vpos.x()/CLHEP::cm;
2198  fMuStartY[index] = vpos.y()/CLHEP::cm;
2199  fMuStartZ[index] = vpos.z()/CLHEP::cm;
2200  fMuStartE[index] = track->GetVertexKineticEnergy()/CLHEP::GeV;
2201 
2202  fMuDE[index] = aStep.GetTotalEnergyDeposit()/CLHEP::MeV;
2203  fMuDEion[index] = fMuDE[index] - aStep.GetNonIonizingEnergyDeposit()/CLHEP::MeV;
2204  fMudEdx[index] = aStep.GetStepLength()/CLHEP::cm;
2205  fMudEdx_ion[index] = aStep.GetStepLength()/CLHEP::cm;
2206  fMuNSteps[index] = 1;
2207  //Parent info
2208  LBNETrajectory* parTraj = GetParentTrajectory(fMuParentID[index]);
2209 
2210  if (parTraj){
2211  int npt = parTraj->GetPointEntries();
2212  fMuParPDG[index] = parTraj->GetPDGEncoding();
2213  fMuParX[index] = parTraj->GetPoint(npt-1)->GetPosition().x()/CLHEP::cm;
2214  fMuParY[index] = parTraj->GetPoint(npt-1)->GetPosition().y()/CLHEP::cm;
2215  fMuParZ[index] = parTraj->GetPoint(npt-1)->GetPosition().z()/CLHEP::cm;
2216  fMuParPX[index] = parTraj->GetMomentum(npt-1).x()/CLHEP::GeV;
2217  fMuParPY[index] = parTraj->GetMomentum(npt-1).y()/CLHEP::GeV;
2218  fMuParPZ[index] = parTraj->GetMomentum(npt-1).z()/CLHEP::GeV;
2219  fMuParE[index] = sqrt(fMuParPX[index]*fMuParPX[index]+fMuParPY[index]*fMuParPY[index]+fMuParZ[index]*fMuParZ[index] + parTraj->GetMass()*parTraj->GetMass() / (CLHEP::GeV*CLHEP::GeV) )-parTraj->GetMass()/ CLHEP::GeV;
2220  }else{
2221  fMuParPDG[index] = 0;
2222  fMuParE[index] = -1;
2223  }
2224  }
2225  }//end if sculpted muon alcove plane is on
2226 
2227 }
2228 //-------------------------------------------------------------------------------------------------
2229 //PPFx stuff again!!! A. Bashyal
2230 G4double LBNEAnalysis::GetDistanceInVolume(LBNETrajectory* wanted_traj, G4String wanted_vol){
2231  double dist_vol = 0;
2232  if(wanted_traj==0)return -1.;
2233 
2234  // G4double fact_Al = (26.98/2.7);
2235 // G4double fact_steel316 = (52.73/8.0);
2236 // G4double fact_concrete = (33.63/2.03);
2237 // G4double fact_water = (18.01/1.0);
2238 // G4double fact_iron = (207.19/11.35);
2239 // G4double fact_helium = (4.003/0.000145);//This from G4numi
2240 
2241  G4ThreeVector ParticlePos;
2242  G4int npoints = wanted_traj->GetPointEntries();
2243 
2244  G4ThreeVector tmp_ipos,tmp_fpos;
2245  G4ThreeVector tmp_disp;
2246 
2247  G4double tmp_dist = 0.0;
2248  G4bool enter_vol = false;
2249  G4bool exit_vol = false;
2250  for(G4int ii=0; ii<npoints; ++ii){
2251  ParticlePos = (wanted_traj->GetPoint(ii)->GetPosition()/CLHEP::m)*CLHEP::m;
2252  G4String postvol = "";
2253  G4String prevol = wanted_traj->GetPreStepVolumeName(ii);
2254  if(ii<npoints-1) postvol = wanted_traj->GetPreStepVolumeName(ii+1);
2255 
2256  G4bool vol_in = (!(prevol.contains(wanted_vol)) && (postvol.contains(wanted_vol)) ) || ( ii==0 && prevol.contains(wanted_vol));
2257  G4bool vol_out = (prevol.contains(wanted_vol)) && (!postvol.contains(wanted_vol));
2258  G4bool vol_through = (prevol.contains(wanted_vol))&&(postvol.contains(wanted_vol));
2259  //Get rid of the vol through and instead create an array of volume of interests and create new ntuples.
2260  // if(prevol.contains("Water"))vol_through = false; //We dont want water to be part of Inner Conductor.
2261  // if(vol_in) std::cout<<"Particle entered at "<<prevol<<" and exited at "<<postvol<<std::endl;
2262  // if((vol_through)&&(prevol != "DecayPipeVolume_P")) std::cout<<"Particle through the volume "<<prevol<<" and "<<postvol<<std::endl;
2263  if(vol_in){
2264  enter_vol = true;
2265  exit_vol = false;
2266  tmp_ipos = G4ThreeVector(ParticlePos[0]/CLHEP::cm,ParticlePos[1]/CLHEP::cm,ParticlePos[2]/CLHEP::cm);
2267  tmp_dist = 0.0;
2268  }
2269  if(enter_vol &&vol_through&& !exit_vol){
2270  tmp_fpos = G4ThreeVector(ParticlePos[0]/CLHEP::cm,ParticlePos[1]/CLHEP::cm,ParticlePos[2]/CLHEP::cm);
2271  tmp_disp = tmp_fpos - tmp_ipos;
2272  tmp_dist += tmp_disp.mag();
2273  tmp_ipos = tmp_fpos;
2274  }
2275  if(enter_vol && vol_out){
2276  tmp_fpos = G4ThreeVector(ParticlePos[0]/CLHEP::cm,ParticlePos[1]/CLHEP::cm,ParticlePos[2]/CLHEP::cm);
2277  tmp_disp = tmp_fpos - tmp_ipos;
2278  tmp_dist += tmp_disp.mag();
2279  tmp_ipos = tmp_fpos;
2280  exit_vol = true;
2281  enter_vol = false;
2282  dist_vol += tmp_dist;
2283  }
2284 
2285  }
2286 
2287 
2288  return dist_vol;
2289 
2290 }
2291 //-------------------------------------------------------------------------------------------
2292 void LBNEAnalysis::ResetEvent()
2293 {
2294  fH1NParticles = 0;
2295  fH2NParticles = 0;
2296  fNParticles = 0;
2297  fDPNParticles = 0;
2298  fTNParticles = 0;
2299  fMuNParticles = 0;
2300  // std::cout << "Resetting event" << std::endl;
2301 }
2302 
2303 
2304 void LBNEAnalysis::FillEvent()
2305 {
2306  if(fNParticles>0){
2307  fTrackingTree->Fill();
2308 }
2309 
2310  if(fTNParticles>0){
2311  fTargetOutputTree->Fill();
2312 }
2313 
2314 
2315  for (int i = 0 ; i < fMuNParticles; i++)
2316  {
2317  fMudEdx[i] = fMuDE[i]/fMudEdx[i];
2318  fMudEdx_ion[i] = fMuDEion[i]/fMudEdx_ion[i];
2319  }
2320  fAlcoveTrackingTree->Fill();
2321 
2322 
2323 }
2324 
2325 void LBNEAnalysis::CalcLocationWeights(const bsim::DkMeta* dkmeta, bsim::Dk2Nu* dk2nu, bool useRealisticNearDetectorVolume)
2326 {
2327 
2328  size_t nloc = dkmeta->location.size();
2329  for (size_t iloc = 0; iloc < nloc; ++iloc ) {
2330  // skip calculation for random location ... should already be filled
2331  const std::string rkey = "random decay";
2332  if ( dkmeta->location[iloc].name == rkey ) {
2333  if ( iloc != 0 ) {
2334  std::cerr << "calcLocationWeights \"" << rkey << "\""
2335  << " isn't the 0-th entry" << std::endl;
2336  assert(0);
2337  }
2338  if ( dk2nu->nuray.size() != 1 ) {
2339  std::cerr << "calcLocationWeights \"" << rkey << "\""
2340  << " nuenergy[" << iloc << "] not filled" << std::endl;
2341  assert(0);
2342  }
2343  continue;
2344  }
2345  double xdet = dkmeta->location[iloc].x;
2346  double ydet = dkmeta->location[iloc].y;
2347  double zdet = dkmeta->location[iloc].z;
2348 
2349  if(useRealisticNearDetectorVolume) {
2350  // randomize location within a volume
2351  // for now, assume a 2 m x 2m x 6 m near detector
2352  xdet += 200*(CLHEP::HepRandom::getTheEngine()->flat()-0.5);
2353  ydet += 200*(CLHEP::HepRandom::getTheEngine()->flat()-0.5);
2354  zdet += 600*(CLHEP::HepRandom::getTheEngine()->flat()-0.5);
2355  }
2356 
2357  TVector3 xyzDet(xdet,
2358  ydet,
2359  zdet); // position to evaluate
2360  double enu_xy = 0; // give a default value
2361  double wgt_xy = 0; // give a default value
2362  int status = bsim::calcEnuWgt(dk2nu,xyzDet,enu_xy,wgt_xy);
2363  if ( status != 0 ) {
2364  std::cerr << "bsim::calcEnuWgt returned " << status << " for "
2365  << dkmeta->location[iloc].name << std::endl;
2366  }
2367  // with the recalculated energy compute the momentum components
2368  TVector3 xyzDk(dk2nu->decay.vx,dk2nu->decay.vy,dk2nu->decay.vz); // origin of decay
2369  TVector3 p3 = enu_xy * (xyzDet - xyzDk).Unit();
2370  bsim::NuRay anuray(p3.x(), p3.y(), p3.z(), enu_xy, wgt_xy);
2371  dk2nu->nuray.push_back(anuray);
2372  }
2373 }
LBNETrajectory::GetVolName1rst
G4String GetVolName1rst() const
Definition: LBNETrajectory.hh:89
LBNEAnalysis::fH2CProcess
G4String fH2CProcess
Definition: LBNEAnalysis.hh:156
LBNEAnalysis::FillTrackingPlaneH1Data
void FillTrackingPlaneH1Data(const G4Step &aStep)
Definition: LBNEAnalysis.cc:1890
LBNEAnalysis::GetCount
G4int GetCount()
Definition: LBNEAnalysis.cc:617
LBNETrajectory::GetTimeStart
G4double GetTimeStart() const
Definition: LBNETrajectory.hh:92
LBNEAnalysis::fParticlePDG
int fParticlePDG[kMaxP]
Definition: LBNEAnalysis.hh:91
LBNEAnalysis::fHorn2TrackingTree
TTree * fHorn2TrackingTree
Definition: LBNEAnalysis.hh:243
LBNEDataNtp_t::muparpy
Float_t muparpy
Definition: LBNEDataNtp_t.hh:153
LBNEAnalysis::fH2ParticleDZ
double fH2ParticleDZ
Definition: LBNEAnalysis.hh:144
LBNEAnalysis::fDPParticleDY
double fDPParticleDY
Definition: LBNEAnalysis.hh:172
LBNEAnalysis::fillDkMeta
void fillDkMeta()
Definition: LBNEAnalysis.cc:377
LBNEAnalysis::dist_DPIP
G4double dist_DPIP[3]
Definition: LBNEAnalysis.hh:299
LBNEDataNtp_t::hornCurrent
Float_t hornCurrent
Definition: LBNEDataNtp_t.hh:120
LBNEAnalysis::fH2PProductionX
double fH2PProductionX
Definition: LBNEAnalysis.hh:150
LBNEAnalysis::fMudEdx
double fMudEdx[kMaxP]
Definition: LBNEAnalysis.hh:225
LBNEDataNtp_t::beamY
Float_t beamY
Definition: LBNEDataNtp_t.hh:112
LBNEAnalysis::fDPPProductionX
double fDPPProductionX
Definition: LBNEAnalysis.hh:179
LBNEAnalysis::GetParentTrajectory
LBNETrajectory * GetParentTrajectory(G4int parentID)
Definition: LBNEAnalysis.cc:1736
LBNEAnalysis::fMuParE
double fMuParE[kMaxP]
Definition: LBNEAnalysis.hh:235
LBNEDataNtp_t::ppvz
Float_t ppvz
Definition: LBNEDataNtp_t.hh:151
LBNEAnalysis::fOutFileDk2Nu
TFile * fOutFileDk2Nu
Definition: LBNEAnalysis.hh:263
LBNEAnalysis::FillNeutrinoNtuple
void FillNeutrinoNtuple(const G4Track &track, const std::vector< G4VTrajectory * > &nuHistory)
Definition: LBNEAnalysis.cc:634
LBNETrajectory::GetPolarization
G4ThreeVector GetPolarization() const
Definition: LBNETrajectory.hh:95
LBNEDataNtp_t::pdPz
Float_t pdPz
Definition: LBNEDataNtp_t.hh:141
entry
QList< Entry > entry
Definition: tclscanner.cpp:1150
LBNEAnalysis::fMuParX
double fMuParX[kMaxP]
Definition: LBNEAnalysis.hh:232
LBNEAnalysis::fH2ParticleX
double fH2ParticleX
Definition: LBNEAnalysis.hh:137
LBNEAnalysis::fParticleZ
double fParticleZ[kMaxP]
Definition: LBNEAnalysis.hh:95
LBNEAnalysis::fDPParticleX
double fDPParticleX
Definition: LBNEAnalysis.hh:166
LBNEDataNtp_t::xpoint
Float_t xpoint
Definition: LBNEDataNtp_t.hh:158
LBNEAnalysis::fH2ParticlePX
double fH2ParticlePX
Definition: LBNEAnalysis.hh:145
LBNENuWeight::GetWeight
double GetWeight(const LBNEDataNtp_t *nudata, const std::vector< double > xdet, double &nu_wght, double &nu_energy)
Definition: LBNENuWeight.cc:25
LBNEAnalysis::fDkMeta
bsim::DkMeta * fDkMeta
Definition: LBNEAnalysis.hh:271
LBNEDataNtp_t::Vy
Float_t Vy
Definition: LBNEDataNtp_t.hh:137
LBNEDataNtp_t::NdxdzNear
Float_t NdxdzNear[5]
Definition: LBNEDataNtp_t.hh:125
LBNERunManager::GetCreateTrkPlaneH1Output
bool GetCreateTrkPlaneH1Output() const
Definition: LBNERunManager.hh:122
LBNEAnalysis::fMuParPDG
int fMuParPDG[kMaxP]
Definition: LBNEAnalysis.hh:231
LBNETrajectory::GetMaterialName
virtual G4String GetMaterialName(G4int i) const
Definition: LBNETrajectory.hh:106
LBNEAnalysis::fMuParPZ
double fMuParPZ[kMaxP]
Definition: LBNEAnalysis.hh:238
LBNEAnalysis::fH1ParticleMass
double fH1ParticleMass
Definition: LBNEAnalysis.hh:112
LBNETrajectory::GetMass
G4double GetMass() const
Definition: LBNETrajectory.hh:54
LBNEDataNtp_t::tgen
Int_t tgen
Definition: LBNEDataNtp_t.hh:168
TrackPoint_t::z
Float_t z
Definition: TrackPoint_t.hh:54
string
std::string string
Definition: nybbler.cc:12
LBNEDataNtp_t::beamVWidth
Float_t beamVWidth
Definition: LBNEDataNtp_t.hh:110
LBNEDataNtp_t::Ntype
Int_t Ntype
Definition: LBNEDataNtp_t.hh:135
LBNEAnalysis::fDPPProductionY
double fDPPProductionY
Definition: LBNEAnalysis.hh:180
LBNEVolumePlacements::Instance
static LBNEVolumePlacements * Instance()
Definition: LBNEVolumePlacements.cc:57
LBNERunManager::GetCreateAlcoveTrackingOutput
bool GetCreateAlcoveTrackingOutput() const
Definition: LBNERunManager.hh:128
LBNEAnalysis::fDk2Nu
bsim::Dk2Nu * fDk2Nu
Definition: LBNEAnalysis.hh:270
LBNEAnalysis::fDPParticleDX
double fDPParticleDX
Definition: LBNEAnalysis.hh:171
LBNEAnalysis::fH1PProductionX
double fH1PProductionX
Definition: LBNEAnalysis.hh:121
LBNEAnalysis::fDPParticleEnergy
double fDPParticleEnergy
Definition: LBNEAnalysis.hh:169
LBNEAnalysis::fXdet_far
std::vector< G4double > fXdet_far
Definition: LBNEAnalysis.hh:286
VersionAndContext::GetGitDescription
std::string GetGitDescription() const
Definition: VersionAndContext.hh:34
LBNEAnalysis::fMudEdx_ion
double fMudEdx_ion[kMaxP]
Definition: LBNEAnalysis.hh:226
TrackPoint_t::type
Int_t type
Definition: TrackPoint_t.hh:48
LBNERunManager::GetDetectorLocationFileName
G4String GetDetectorLocationFileName() const
Definition: LBNERunManager.hh:118
LBNEAnalysis::fDetNameNear
std::vector< G4String > fDetNameNear
Definition: LBNEAnalysis.hh:289
LBNEDataNtp_t::ppenergy
Float_t ppenergy
Definition: LBNEDataNtp_t.hh:145
TrackPoint_t::trkid
Int_t trkid
Definition: TrackPoint_t.hh:49
LBNEDataNtp_t::protonPz
Float_t protonPz
Definition: LBNEDataNtp_t.hh:118
LBNERunManager::GetPhysicsListName
G4String GetPhysicsListName() const
Definition: LBNERunManager.hh:116
LBNEAnalysis::CloseDk2NuOutput
void CloseDk2NuOutput()
Definition: LBNEAnalysis.cc:591
LBNEAnalysis::NuMomentum
G4ThreeVector NuMomentum
Definition: LBNEAnalysis.hh:128
LBNEAnalysis::fMuT
double fMuT[kMaxP]
Definition: LBNEAnalysis.hh:214
LBNERunManager::GetCreateASCIIOutput
bool GetCreateASCIIOutput() const
Definition: LBNERunManager.hh:95
LBNEAnalysis::fParticleMass
double fParticleMass[kMaxP]
Definition: LBNEAnalysis.hh:97
LBNEAnalysis::FillTrackingPlaneH2Data
void FillTrackingPlaneH2Data(const G4Step &aStep)
Definition: LBNEAnalysis.cc:1947
LBNEDataNtp_t::mupare
Float_t mupare
Definition: LBNEDataNtp_t.hh:155
LBNEDataNtp_t::pdPy
Float_t pdPy
Definition: LBNEDataNtp_t.hh:140
LBNEAnalysis::fH1ParticlePX
double fH1ParticlePX
Definition: LBNEAnalysis.hh:116
LBNEDataNtp_t::nuTarZ
Float_t nuTarZ
Definition: LBNEDataNtp_t.hh:119
LBNEVolumePlacements::GetHorn2ICList
std::vector< G4String > GetHorn2ICList() const
Definition: LBNEVolumePlacements.hh:345
LBNEAnalysis::GetTrajectory
LBNETrajectory * GetTrajectory(G4int trackID)
Definition: LBNEAnalysis.cc:1764
LBNEDataNtp_t::tpy
Float_t tpy
Definition: LBNEDataNtp_t.hh:165
LBNEAnalysis::fDPParticleMass
double fDPParticleMass
Definition: LBNEAnalysis.hh:170
LBNEAnalysis::fDPParticleDZ
double fDPParticleDZ
Definition: LBNEAnalysis.hh:173
LBNEAnalysis::fDPParticleY
double fDPParticleY
Definition: LBNEAnalysis.hh:167
std
STL namespace.
LBNEAnalysis::fH2ParticleMass
double fH2ParticleMass
Definition: LBNEAnalysis.hh:141
LBNERunManager::GetCreateTrkPlaneDPOutput
bool GetCreateTrkPlaneDPOutput() const
Definition: LBNERunManager.hh:124
LBNEAnalysis::fMuPX
double fMuPX[kMaxP]
Definition: LBNEAnalysis.hh:221
LBNEAnalysis::dist_DVOL
G4double dist_DVOL[3]
Definition: LBNEAnalysis.hh:300
LBNEAnalysis::fH1ParticlePXPZ
double fH1ParticlePXPZ
Definition: LBNEAnalysis.hh:119
LBNEAnalysis::fH1ParticleZ
double fH1ParticleZ
Definition: LBNEAnalysis.hh:110
LBNEDataNtp_t::pdPx
Float_t pdPx
Definition: LBNEDataNtp_t.hh:139
LBNEAnalysis::SetEntry
void SetEntry(G4int entry)
Definition: LBNEAnalysis.cc:622
LBNEVolumePlacements::GetTargetSLengthGraphite
double GetTargetSLengthGraphite() const
Definition: LBNEVolumePlacements.hh:255
LBNEAnalysis::fOutTreeDk2NuMeta
TTree * fOutTreeDk2NuMeta
Definition: LBNEAnalysis.hh:265
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:45
LBNEAnalysis::y
G4double y
Definition: LBNEAnalysis.hh:85
LBNEAnalysis::fMuStartZ
double fMuStartZ[kMaxP]
Definition: LBNEAnalysis.hh:217
LBNEDataNtp_t::Nimpwt
Double_t Nimpwt
Definition: LBNEDataNtp_t.hh:157
LBNEAnalysis::SetCount
void SetCount(G4int count)
Definition: LBNEAnalysis.cc:612
LBNEQuickPiToNuVect::blessPion
void blessPion(int trNum1, double e, G4ThreeVector p)
Definition: LBNEQuickPiToNu.hh:100
LBNETrajectory::GetDecayCode
G4int GetDecayCode() const
Definition: LBNETrajectory.hh:72
LBNEAnalysis::fParticleEnergy
double fParticleEnergy[kMaxP]
Definition: LBNEAnalysis.hh:96
LBNERunManager::GetCreateTargetOutput
bool GetCreateTargetOutput() const
Definition: LBNERunManager.hh:120
VersionAndContext::Instance
static VersionAndContext * Instance()
Definition: VersionAndContext.cc:17
LBNEAnalysis::fYdet_far
std::vector< G4double > fYdet_far
Definition: LBNEAnalysis.hh:287
LBNEAnalysis::fH2PProductionDX
double fH2PProductionDX
Definition: LBNEAnalysis.hh:153
LBNEParticleCode::StringToEnum
LBNEParticleCode_t StringToEnum(G4String particleName)
Definition: LBNEParticleCode.cc:90
LBNEAnalysis::x
G4double x
Definition: LBNEAnalysis.hh:84
LBNEAnalysis::fH1PProductionZ
double fH1PProductionZ
Definition: LBNEAnalysis.hh:123
LBNEVolumePlacements::GetDecayPipeLength
double GetDecayPipeLength() const
Definition: LBNEVolumePlacements.hh:165
LBNEAnalysis::fMuPDG
int fMuPDG[kMaxP]
Definition: LBNEAnalysis.hh:209
LBNEAnalysis::fDk2NuDetectorFileRead
bool fDk2NuDetectorFileRead
Definition: LBNEAnalysis.hh:82
LBNEDataNtp_t::protonY
Float_t protonY
Definition: LBNEDataNtp_t.hh:114
LBNEAnalysis::FillTrackingNtuple
void FillTrackingNtuple(const G4Track &track, LBNETrajectory *currTrajectory)
Definition: LBNEAnalysis.cc:1517
LBNEAnalysis::fMuX
double fMuX[kMaxP]
Definition: LBNEAnalysis.hh:211
LBNEAnalysis::TrackThroughGeometry
void TrackThroughGeometry(const LBNETrajectory *TrackTrajectory)
Definition: LBNEAnalysis.cc:1583
LBNEAnalysis::fMuParPX
double fMuParPX[kMaxP]
Definition: LBNEAnalysis.hh:236
LBNEAnalysis::FillAlcoveTrackingPlaneData
void FillAlcoveTrackingPlaneData(const G4Step &aStep)
Definition: LBNEAnalysis.cc:2130
LBNETrajectory::GetTrackID
G4int GetTrackID() const
Definition: LBNETrajectory.hh:46
LBNEVolumePlacements::GetHorn1NeckInnerRadius
double GetHorn1NeckInnerRadius() const
Definition: LBNEVolumePlacements.hh:320
LBNEAnalysis::fTParticlePY
double fTParticlePY
Definition: LBNEAnalysis.hh:196
LBNEDataNtp_t::protonPx
Float_t protonPx
Definition: LBNEDataNtp_t.hh:116
LBNEDataNtp_t::ppdydz
Float_t ppdydz
Definition: LBNEDataNtp_t.hh:143
LBNEAnalysis::VolVintName
vstring_t VolVintName
Definition: LBNEAnalysis.hh:303
sensitivity_comps.i
int i
Definition: sensitivity_comps.py:15
LBNEAnalysis::fH1PProductionDX
double fH1PProductionDX
Definition: LBNEAnalysis.hh:124
LBNEAnalysis::fLBNEOutNtpData
LBNEDataNtp_t * fLBNEOutNtpData
Definition: LBNEAnalysis.hh:273
LBNEAnalysis::fDPParticlePYPZ
double fDPParticlePYPZ
Definition: LBNEAnalysis.hh:178
LBNETrajectory::GetPDGNucleus
G4int GetPDGNucleus() const
Definition: LBNETrajectory.hh:98
LBNEAnalysis::fXdet_near
std::vector< G4double > fXdet_near
Definition: LBNEAnalysis.hh:283
LBNEAnalysis::fDPPProductionDX
double fDPPProductionDX
Definition: LBNEAnalysis.hh:182
LBNEAnalysis::fDPPProductionZ
double fDPPProductionZ
Definition: LBNEAnalysis.hh:181
LBNEAnalysis::fMuMass
double fMuMass[kMaxP]
Definition: LBNEAnalysis.hh:219
LBNEVolumePlacements::GetDecayPipeGas
G4String GetDecayPipeGas() const
Definition: LBNEVolumePlacements.hh:169
LBNEAnalysis::fH2PProductionZ
double fH2PProductionZ
Definition: LBNEAnalysis.hh:152
LBNEAnalysis::getInstance
static LBNEAnalysis * getInstance()
Definition: LBNEAnalysis.cc:134
LBNEAnalysis::instance
static LBNEAnalysis * instance
Definition: LBNEAnalysis.hh:80
LBNEAnalysis::fH1PProductionDZ
double fH1PProductionDZ
Definition: LBNEAnalysis.hh:126
VersionAndContext::GetG4lbnfDir
std::string GetG4lbnfDir() const
Definition: VersionAndContext.hh:37
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
LBNEAnalysis::setDetectorPositions
void setDetectorPositions()
Definition: LBNEAnalysis.cc:1789
LBNEDataNtp_t::NWtFar
Double_t NWtFar[3]
Definition: LBNEDataNtp_t.hh:132
LBNEAnalysis::CalcLocationWeights
void CalcLocationWeights(const bsim::DkMeta *dkmeta, bsim::Dk2Nu *dk2nu, bool useRealisticNearDetectorVolume)
Definition: LBNEAnalysis.cc:2325
LBNEAnalysis::fDPTrackingTree
TTree * fDPTrackingTree
Definition: LBNEAnalysis.hh:244
LBNETrajectory::GetNImpWt
virtual G4double GetNImpWt() const
Definition: LBNETrajectory.hh:74
LBNERunManager::GetNumberOfEventsLBNE
int GetNumberOfEventsLBNE()
Definition: LBNERunManager.hh:133
LBNEAnalysis::fMuPY
double fMuPY[kMaxP]
Definition: LBNEAnalysis.hh:222
LBNERunManager::GetOutputASCIIFileName
G4String GetOutputASCIIFileName() const
Definition: LBNERunManager.hh:89
LBNEAnalysis::fTParticleZ
double fTParticleZ
Definition: LBNEAnalysis.hh:192
LBNEAnalysis::fMuEnergy
double fMuEnergy[kMaxP]
Definition: LBNEAnalysis.hh:220
LBNEDataNtp_t::tpz
Float_t tpz
Definition: LBNEDataNtp_t.hh:166
LBNEAnalysis::fMuPZ
double fMuPZ[kMaxP]
Definition: LBNEAnalysis.hh:223
LBNEAnalysis::fDPParticlePX
double fDPParticlePX
Definition: LBNEAnalysis.hh:174
LBNEAnalysis::GetEntry
G4int GetEntry()
Definition: LBNEAnalysis.cc:627
LBNEDataNtp_t::protonPy
Float_t protonPy
Definition: LBNEDataNtp_t.hh:117
LBNEDataNtp_t::tvx
Float_t tvx
Definition: LBNEDataNtp_t.hh:161
LBNERunAction::GetDoComputeEDepInHorns
bool GetDoComputeEDepInHorns() const
Definition: LBNERunAction.hh:62
LBNEVolumePlacements::GetHorn2NeckInnerRadius
double GetHorn2NeckInnerRadius() const
Definition: LBNEVolumePlacements.hh:331
LBNERunManager::GetCreateTrkPlaneOutput
bool GetCreateTrkPlaneOutput() const
Definition: LBNERunManager.hh:97
LBNETrajectory::GetPreStepVolumeName
virtual G4String GetPreStepVolumeName(G4int i) const
Definition: LBNETrajectory.hh:66
LBNEAnalysis::fMuTrackID
int fMuTrackID[kMaxP]
Definition: LBNEAnalysis.hh:207
LBNEAnalysis::FillTrackingPlaneData
void FillTrackingPlaneData(const G4Step &aStep)
Definition: LBNEAnalysis.cc:1853
LBNEAnalysis::fMuZ
double fMuZ[kMaxP]
Definition: LBNEAnalysis.hh:213
LBNEAnalysis::fParticleY
double fParticleY[kMaxP]
Definition: LBNEAnalysis.hh:94
LBNEAnalysis::fH2PProductionDZ
double fH2PProductionDZ
Definition: LBNEAnalysis.hh:155
cet::getenv
std::string getenv(std::string const &name)
Definition: getenv.cc:15
TrackPoint_t::impwt
Double_t impwt
Definition: TrackPoint_t.hh:51
LBNEAnalysis::fH1ParticleDX
double fH1ParticleDX
Definition: LBNEAnalysis.hh:113
LBNEAnalysis::fMuParY
double fMuParY[kMaxP]
Definition: LBNEAnalysis.hh:233
LBNEDataNtp_t::Ndxdz
Float_t Ndxdz
Definition: LBNEDataNtp_t.hh:121
LBNEDataNtp_t::protonX
Float_t protonX
Definition: LBNEDataNtp_t.hh:113
LBNEAnalysis::fMuY
double fMuY[kMaxP]
Definition: LBNEAnalysis.hh:212
generate_datataset.info
info
Definition: generate_datataset.py:142
LBNETrajectory::GetParentMomentumAtThisProduction
const G4ThreeVector GetParentMomentumAtThisProduction() const
Definition: LBNETrajectory.hh:108
LBNEAnalysis::fOutTreeDk2Nu
TTree * fOutTreeDk2Nu
Definition: LBNEAnalysis.hh:264
LBNEAnalysis::fH1ParticlePY
double fH1ParticlePY
Definition: LBNEAnalysis.hh:117
LBNEAnalysis::fZdet_far
std::vector< G4double > fZdet_far
Definition: LBNEAnalysis.hh:288
LBNEAnalysis::nuNtupleFileName
char nuNtupleFileName[1024]
Definition: LBNEAnalysis.hh:250
LBNEAnalysis::fDPParticleZ
double fDPParticleZ
Definition: LBNEAnalysis.hh:168
LBNEDataNtp_t::Norig
Int_t Norig
Definition: LBNEDataNtp_t.hh:133
LBNEDataNtp_t::ppvy
Float_t ppvy
Definition: LBNEDataNtp_t.hh:150
LBNEAnalysis::fMuDE
double fMuDE[kMaxP]
Definition: LBNEAnalysis.hh:227
LBNEQuickPiToNuVect::Instance
static LBNEQuickPiToNuVect * Instance()
Definition: LBNEQuickPiToNu.cc:21
TrkPoint::StringToEnum
TrkPoint_t StringToEnum(const std::string &trkpt)
Definition: LBNEDataNtp_t.cc:465
LBNEAnalysis::fH1ParticlePZ
double fH1ParticlePZ
Definition: LBNEAnalysis.hh:118
LBNETrajectory::GetPoint
virtual G4VTrajectoryPoint * GetPoint(G4int i) const
Definition: LBNETrajectory.hh:62
LBNEAnalysis::fDPParticlePY
double fDPParticlePY
Definition: LBNEAnalysis.hh:175
LBNEAnalysis::fMuStartY
double fMuStartY[kMaxP]
Definition: LBNEAnalysis.hh:216
LBNEDataNtp_t::h1trackid
Int_t h1trackid
Definition: LBNEDataNtp_t.hh:182
LBNETrajectory::GetTgen
virtual G4int GetTgen() const
Definition: LBNETrajectory.hh:70
LBNEDataNtp_t::NdxdzFar
Float_t NdxdzFar[3]
Definition: LBNEDataNtp_t.hh:129
LBNEDataNtp_t::muparpx
Float_t muparpx
Definition: LBNEDataNtp_t.hh:152
LBNEAnalysis::fH2PProductionDY
double fH2PProductionDY
Definition: LBNEAnalysis.hh:154
LBNEDataNtp_t::tpx
Float_t tpx
Definition: LBNEDataNtp_t.hh:164
LBNEAnalysis::fH1ParticleDZ
double fH1ParticleDZ
Definition: LBNEAnalysis.hh:115
LBNEAnalysis::fH1NImpWt
double fH1NImpWt
Definition: LBNEAnalysis.hh:105
LBNEAnalysis::fH2ParticlePXPZ
double fH2ParticlePXPZ
Definition: LBNEAnalysis.hh:148
LBNEDataNtp_t::beamX
Float_t beamX
Definition: LBNEDataNtp_t.hh:111
LBNEAnalysis::z
G4double z
Definition: LBNEAnalysis.hh:86
LBNEDataNtp_t::Necm
Float_t Necm
Definition: LBNEDataNtp_t.hh:156
LBNEAnalysis::fAlcoveTrackingTree
TTree * fAlcoveTrackingTree
Definition: LBNEAnalysis.hh:246
LBNEDataNtp_t::NdydzFar
Float_t NdydzFar[3]
Definition: LBNEDataNtp_t.hh:130
LBNEDataNtp_t::tvy
Float_t tvy
Definition: LBNEDataNtp_t.hh:162
LBNEAnalysis::fDPPProductionDZ
double fDPPProductionDZ
Definition: LBNEAnalysis.hh:184
LBNEDataNtp_t::NenergyF
Float_t NenergyF[3]
Definition: LBNEDataNtp_t.hh:131
LBNEAnalysis::fParticleDZ
double fParticleDZ[kMaxP]
Definition: LBNEAnalysis.hh:100
LBNERunManager::GetLBNESteppingManager
LBNESteppingAction * GetLBNESteppingManager()
Definition: LBNERunManager.hh:74
LBNEAnalysis::fTParticleEnergy
double fTParticleEnergy
Definition: LBNEAnalysis.hh:194
LBNERunManager::GetOutputNtpFileName
G4String GetOutputNtpFileName() const
Definition: LBNERunManager.hh:85
LBNETrajectory::GetProcessName
G4String GetProcessName() const
Definition: LBNETrajectory.hh:101
LBNEAnalysis::fMuStartE
double fMuStartE[kMaxP]
Definition: LBNEAnalysis.hh:218
LBNEAnalysis::fH1PProductionDY
double fH1PProductionDY
Definition: LBNEAnalysis.hh:125
LBNETrajectory::GetMaterialNumber1rst
G4int GetMaterialNumber1rst() const
Definition: LBNETrajectory.hh:79
LBNEAnalysis::fH2ParticleDY
double fH2ParticleDY
Definition: LBNEAnalysis.hh:143
LBNEAnalysis::fH2NImpWt
double fH2NImpWt
Definition: LBNEAnalysis.hh:136
LBNEAnalysis::fH1ParticleY
double fH1ParticleY
Definition: LBNEAnalysis.hh:109
LBNEAnalysis::fTrackID
int fTrackID[kMaxP]
Definition: LBNEAnalysis.hh:90
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
LBNEDataNtp_t
Definition: LBNEDataNtp_t.hh:58
TrackPoint_t::px
Float_t px
Definition: TrackPoint_t.hh:55
LBNEDataNtp_t::Vz
Float_t Vz
Definition: LBNEDataNtp_t.hh:138
a
const GenericPointer< typename T::ValueType > T2 T::AllocatorType & a
Definition: pointer.h:1124
LBNEVolumePlacements::GetHorn1NeckZPosition
double GetHorn1NeckZPosition() const
Definition: LBNEVolumePlacements.hh:317
LBNEVolumePlacements::GetHorn1PolyInnerConductorRadius
double GetHorn1PolyInnerConductorRadius(size_t numPts) const
Definition: LBNEVolumePlacements.hh:454
LBNEAnalysis::fH1ParticleX
double fH1ParticleX
Definition: LBNEAnalysis.hh:108
LBNEAnalysis::fTParticleY
double fTParticleY
Definition: LBNEAnalysis.hh:191
LBNERunManager::GetCreateOutput
bool GetCreateOutput() const
Definition: LBNERunManager.hh:91
LBNEAnalysis::fH2ParticlePYPZ
double fH2ParticlePYPZ
Definition: LBNEAnalysis.hh:149
LBNEAnalysis::fH1ParticleDY
double fH1ParticleDY
Definition: LBNEAnalysis.hh:114
LBNEAnalysis::fDPPProductionDY
double fDPPProductionDY
Definition: LBNEAnalysis.hh:183
LBNEAnalysis::fH2ParticlePZ
double fH2ParticlePZ
Definition: LBNEAnalysis.hh:147
LBNEAnalysis::fTParticlePZ
double fTParticlePZ
Definition: LBNEAnalysis.hh:197
LBNEAnalysis::VolAbsName
vstring_t VolAbsName
Definition: LBNEAnalysis.hh:302
LBNEDataNtp_t::tvz
Float_t tvz
Definition: LBNEDataNtp_t.hh:163
LBNEAnalysis::fParticleDY
double fParticleDY[kMaxP]
Definition: LBNEAnalysis.hh:99
LBNEAnalysis::fMuParPY
double fMuParPY[kMaxP]
Definition: LBNEAnalysis.hh:237
LBNEDataNtp_t::muparpz
Float_t muparpz
Definition: LBNEDataNtp_t.hh:154
LBNEAnalysis::FillTrackingPlaneDPData
void FillTrackingPlaneDPData(const G4Step &aStep)
Definition: LBNEAnalysis.cc:2039
LBNEDataNtp_t::Ndecay
Int_t Ndecay
Definition: LBNEDataNtp_t.hh:134
TrackPoint_t::py
Float_t py
Definition: TrackPoint_t.hh:56
LBNEAnalysis::fH1PProductionY
double fH1PProductionY
Definition: LBNEAnalysis.hh:122
VersionAndContext::GetBuildTimeStamp
std::string GetBuildTimeStamp() const
Definition: VersionAndContext.hh:40
LBNEAnalysis::GenAbsName
vstring_t GenAbsName
Definition: LBNEAnalysis.hh:304
reverse
static unsigned int reverse(QString &chars, unsigned char *level, unsigned int a, unsigned int b)
Definition: qstring.cpp:11649
LBNEVolumePlacements::GetHorn1PolyInnerConductorZCoord
double GetHorn1PolyInnerConductorZCoord(size_t numPts) const
Definition: LBNEVolumePlacements.hh:468
LBNEAnalysis::fYdet_near
std::vector< G4double > fYdet_near
Definition: LBNEAnalysis.hh:284
LBNEAnalysis::fTrackingTree
TTree * fTrackingTree
Definition: LBNEAnalysis.hh:267
LBNEVolumePlacements::GetTargetLengthIntoHorn
double GetTargetLengthIntoHorn() const
Definition: LBNEVolumePlacements.hh:282
LBNEDataNtp_t::Nenergy
Float_t Nenergy
Definition: LBNEDataNtp_t.hh:124
LBNETrajectory::GetPDGEncoding
G4int GetPDGEncoding() const
Definition: LBNETrajectory.hh:56
LBNEAnalysis::fMuNSteps
int fMuNSteps[kMaxP]
Definition: LBNEAnalysis.hh:230
LBNETrajectory::GetParentID
G4int GetParentID() const
Definition: LBNETrajectory.hh:48
LBNEAnalysis::fH2ParticlePY
double fH2ParticlePY
Definition: LBNEAnalysis.hh:146
LBNEDataNtp_t::protonZ
Float_t protonZ
Definition: LBNEDataNtp_t.hh:115
LBNEDataNtp_t::ptype
Int_t ptype
Definition: LBNEDataNtp_t.hh:147
LBNETrackInformation::GetNImpWt
G4double GetNImpWt() const
Definition: LBNETrackInformation.hh:39
LBNERunManager::GetOutputDk2NuFileName
G4String GetOutputDk2NuFileName() const
Definition: LBNERunManager.hh:87
LBNEAnalysis::fH1ParticlePYPZ
double fH1ParticlePYPZ
Definition: LBNEAnalysis.hh:120
run_cnn_1class.m
m
Definition: run_cnn_1class.py:75
LBNERunManager::GetVerboseLevel
int GetVerboseLevel() const
Definition: LBNERunManager.hh:135
LBNEAnalysis::fDPParticlePZ
double fDPParticlePZ
Definition: LBNEAnalysis.hh:176
LBNEDataNtp_t::NdydzNear
Float_t NdydzNear[5]
Definition: LBNEDataNtp_t.hh:126
LBNEDataNtp_t::tptype
Int_t tptype
Definition: LBNEDataNtp_t.hh:167
LBNEAnalysis::CreateDk2NuOutput
G4bool CreateDk2NuOutput()
Definition: LBNEAnalysis.cc:334
LBNEDataNtp_t::evtno
Int_t evtno
Definition: LBNEDataNtp_t.hh:107
LBNEAnalysis::fTargetOutputTree
TTree * fTargetOutputTree
Definition: LBNEAnalysis.hh:245
HkMaxP
const int HkMaxP
Definition: LBNEAnalysis.hh:33
LBNEAnalysis::fH2ParticleZ
double fH2ParticleZ
Definition: LBNEAnalysis.hh:139
LBNEDataNtp_t::run
Int_t run
Definition: LBNEDataNtp_t.hh:106
LBNEAnalysis::fDPNImpWt
double fDPNImpWt
Definition: LBNEAnalysis.hh:165
LBNEAnalysis::fTParticlePX
double fTParticlePX
Definition: LBNEAnalysis.hh:195
LBNEAnalysis::VolVdblName
vstring_t VolVdblName
Definition: LBNEAnalysis.hh:301
LBNETrajectory::GetMaterialName1rst
G4String GetMaterialName1rst() const
Definition: LBNETrajectory.hh:104
LBNEAnalysis::dist_IC2
G4double dist_IC2[3]
Definition: LBNEAnalysis.hh:298
LBNEAnalysis::fOutTree
TTree * fOutTree
Definition: LBNEAnalysis.hh:259
LBNETrajectory::GetParticleName
G4String GetParticleName() const
Definition: LBNETrajectory.hh:50
LBNEAnalysis::dist_IC1
G4double dist_IC1[3]
Definition: LBNEAnalysis.hh:297
LBNEDataNtp_t::NenergyN
Float_t NenergyN[5]
Definition: LBNEDataNtp_t.hh:127
LBNERunManager::GetUseRealisticNearDetectorVolume
bool GetUseRealisticNearDetectorVolume() const
Definition: LBNERunManager.hh:99
LBNEDataNtp_t::ppvx
Float_t ppvx
Definition: LBNEDataNtp_t.hh:149
LBNEAnalysis::fHorn1TrackingTree
TTree * fHorn1TrackingTree
Definition: LBNEAnalysis.hh:242
LBNEDataNtp_t::NWtNear
Double_t NWtNear[5]
Definition: LBNEDataNtp_t.hh:128
LBNEAnalysis::CreateOutput
G4bool CreateOutput()
Definition: LBNEAnalysis.cc:140
LBNEAnalysis::fMuStartX
double fMuStartX[kMaxP]
Definition: LBNEAnalysis.hh:215
LBNEAnalysis::fMuParZ
double fMuParZ[kMaxP]
Definition: LBNEAnalysis.hh:234
LBNEAnalysis::GetDistanceInVolume
G4double GetDistanceInVolume(LBNETrajectory *wanted_traj, G4String wanted_vol)
Definition: LBNEAnalysis.cc:2230
TrackPoint_t
Definition: TrackPoint_t.hh:23
LBNEDataNtp_t::Vx
Float_t Vx
Definition: LBNEDataNtp_t.hh:136
LBNEAnalysis::fH1ParticleEnergy
double fH1ParticleEnergy
Definition: LBNEAnalysis.hh:111
LBNEVolumePlacements::GetNumberOfHornsPolycone
int GetNumberOfHornsPolycone() const
Definition: LBNEVolumePlacements.hh:518
LBNEAnalysis::FillTargetOutputData
void FillTargetOutputData(const G4Step &aStep)
Definition: LBNEAnalysis.cc:2091
LBNEAnalysis::fDPParticlePXPZ
double fDPParticlePXPZ
Definition: LBNEAnalysis.hh:177
LBNEAnalysis::fMuNimpWt
double fMuNimpWt[kMaxP]
Definition: LBNEAnalysis.hh:210
LBNEAnalysis::fH2PProductionY
double fH2PProductionY
Definition: LBNEAnalysis.hh:151
LBNETrajectory::GetPointEntries
virtual int GetPointEntries() const
Definition: LBNETrajectory.hh:60
if
if(!yymsg) yymsg
LBNEAnalysis::fParticleX
double fParticleX[kMaxP]
Definition: LBNEAnalysis.hh:93
TrackPoint_t::gen
Int_t gen
Definition: TrackPoint_t.hh:50
LBNERunManager::GetCreateDk2NuOutput
bool GetCreateDk2NuOutput() const
Definition: LBNERunManager.hh:93
LBNEDataNtp_t::protonN
Int_t protonN
Definition: LBNEDataNtp_t.hh:108
ns
QAsciiDict< Entry > ns
Definition: tclscanner.cpp:1147
LBNEDataNtp_t::Npz
Float_t Npz
Definition: LBNEDataNtp_t.hh:123
TrackPoint_t::pz
Float_t pz
Definition: TrackPoint_t.hh:57
TrackPoint_t::x
Float_t x
Definition: TrackPoint_t.hh:52
LBNEAnalysis::fH2ParticleY
double fH2ParticleY
Definition: LBNEAnalysis.hh:138
LBNETrajectory::GetMomentum
virtual G4ThreeVector GetMomentum(G4int i) const
Definition: LBNETrajectory.hh:64
LBNEDataNtp_t::Ndydz
Float_t Ndydz
Definition: LBNEDataNtp_t.hh:122
LBNEParticleCode::AsInt
G4int AsInt(LBNEParticleCode_t pCode)
Definition: LBNEParticleCode.cc:50
LBNEAnalysis::fMuDEion
double fMuDEion[kMaxP]
Definition: LBNEAnalysis.hh:228
LBNEAnalysis::fNImpWt
double fNImpWt[kMaxP]
Definition: LBNEAnalysis.hh:92
kMaxP
const int kMaxP
Definition: LBNEAnalysis.hh:32
LBNEAnalysis::fH2ParticleDX
double fH2ParticleDX
Definition: LBNEAnalysis.hh:142
LBNEDataNtp_t::ptrkid
Int_t ptrkid
Definition: LBNEDataNtp_t.hh:148
LBNEAnalysis::fMuTheta
double fMuTheta[kMaxP]
Definition: LBNEAnalysis.hh:224
LBNEAnalysis::fTrackingPlaneData
LBNEDataNtp_t * fTrackingPlaneData
Definition: LBNEAnalysis.hh:274
LBNEDataNtp_t.hh
LBNEAnalysis::fDetNameFar
std::vector< G4String > fDetNameFar
Definition: LBNEAnalysis.hh:290
TrackPoint_t::y
Float_t y
Definition: TrackPoint_t.hh:53
LBNEDataNtp_t::beamHWidth
Float_t beamHWidth
Definition: LBNEDataNtp_t.hh:109
LBNEVolumePlacements::GetHorn1ICList
std::vector< G4String > GetHorn1ICList() const
Definition: LBNEVolumePlacements.hh:344
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
LBNEDataNtp_t::h2trackid
Int_t h2trackid
Definition: LBNEDataNtp_t.hh:183
LBNEAnalysis::fTParticleX
double fTParticleX
Definition: LBNEAnalysis.hh:190
LBNEAnalysis::fParticleDX
double fParticleDX[kMaxP]
Definition: LBNEAnalysis.hh:98
LBNESteppingAction::GetKillTrackingThreshold
double GetKillTrackingThreshold() const
Definition: LBNESteppingAction.hh:144
LBNEDataNtp_t::ppdxdz
Float_t ppdxdz
Definition: LBNEDataNtp_t.hh:142
LBNEAnalysis::CloseOutput
void CloseOutput()
Definition: LBNEAnalysis.cc:522
LBNEDataNtp_t::pppz
Float_t pppz
Definition: LBNEDataNtp_t.hh:144
LBNEAnalysis::fZdet_near
std::vector< G4double > fZdet_near
Definition: LBNEAnalysis.hh:285
LBNEAnalysis::fOutFile
TFile * fOutFile
Definition: LBNEAnalysis.hh:258
LBNEAnalysis::fMuParentID
int fMuParentID[kMaxP]
Definition: LBNEAnalysis.hh:208
LBNEAnalysis::fH2ParticleEnergy
double fH2ParticleEnergy
Definition: LBNEAnalysis.hh:140
LBNEDataNtp_t::ppmedium
Float_t ppmedium
Definition: LBNEDataNtp_t.hh:146
LBNERunManager::GetCreateTrkPlaneH2Output
bool GetCreateTrkPlaneH2Output() const
Definition: LBNERunManager.hh:123
Generated by   doxygen 1.8.11