makeCAF.cxx
Go to the documentation of this file.
1 #include "CAF.C"
2 #include "TRandom3.h"
3 #include "TFile.h"
4 #include "TTree.h"
5 #include "TVector3.h"
6 #include "TLorentzVector.h"
8 #include "EVGCore/EventRecord.h"
9 #include "nusystematics/artless/response_helper.hh"
10 #include <stdio.h>
11 
12 TRandom3 * rando;
13 const double mmu = 0.1056583745;
14 TF1 * tsmear; // angular resolution function
15 
16 // params will be extracted from command line, and passed to the reconstruction
17 struct params {
18  double OA_xcoord;
19  bool fhc, grid, IsGasTPC;
20  int seed, run, subrun, first, n, nfiles;
22  double em_const, em_sqrtE;
23  double michelEff;
24  double CC_trk_length;
27 };
28 
29 // Fill reco variables for muon reconstructed in magnetized tracker
30 void recoMuonTracker( CAF &caf, params &par )
31 {
32  // smear momentum by resolution
33  double p = sqrt(caf.LepE*caf.LepE - mmu*mmu);
34  double reco_p = rando->Gaus( p, p*par.trk_muRes );
35  caf.Elep_reco = sqrt(reco_p*reco_p + mmu*mmu);
36 
37  double true_tx = 1000.*atan(caf.LepMomX / caf.LepMomZ);
38  double true_ty = 1000.*atan(caf.LepMomY / caf.LepMomZ);
39  double evalTsmear = tsmear->Eval(caf.Elep_reco - mmu);
40  if( evalTsmear < 0. ) evalTsmear = 0.;
41  double reco_tx = true_tx + rando->Gaus(0., evalTsmear/sqrt(2.));
42  double reco_ty = true_ty + rando->Gaus(0., evalTsmear/sqrt(2.));
43  caf.theta_reco = 0.001*sqrt( reco_tx*reco_tx + reco_ty*reco_ty );
44 
45  // assume perfect charge reconstruction
46  caf.reco_q = (caf.LepPDG > 0 ? -1 : 1);
47 
48  // assume always muon for tracker-matched
49  caf.reco_numu = 1; caf.reco_nue = 0; caf.reco_nc = 0;
50  caf.muon_contained = 0; caf.muon_tracker = 1; caf.muon_ecal = 0; caf.muon_exit = 0;
51  caf.Ev_reco = caf.Elep_reco;
52 }
53 
54 // Fill reco muon variables for muon contained in LAr
55 void recoMuonLAr( CAF &caf, params &par )
56 {
57  // range-based, smear kinetic energy
58  double ke = caf.LepE - mmu;
59  double reco_ke = rando->Gaus( ke, ke*par.LAr_muRes );
60  caf.Elep_reco = reco_ke + mmu;
61 
62  double true_tx = 1000.*atan(caf.LepMomX / caf.LepMomZ);
63  double true_ty = 1000.*atan(caf.LepMomY / caf.LepMomZ);
64  double evalTsmear = tsmear->Eval(caf.Elep_reco - mmu);
65  if( evalTsmear < 0. ) evalTsmear = 0.;
66  double reco_tx = true_tx + rando->Gaus(0., evalTsmear/sqrt(2.));
67  double reco_ty = true_ty + rando->Gaus(0., evalTsmear/sqrt(2.));
68  caf.theta_reco = 0.001*sqrt( reco_tx*reco_tx + reco_ty*reco_ty );
69 
70 
71  // assume negative for FHC, require Michel for RHC
72  if( par.fhc ) caf.reco_q = -1;
73  else {
74  double michel = rando->Rndm();
75  if( caf.LepPDG == -13 && michel < par.michelEff ) caf.reco_q = 1; // correct mu+
76  else if( caf.LepPDG == 13 && michel < par.michelEff*0.25 ) caf.reco_q = 1; // incorrect mu-
77  else caf.reco_q = -1; // no reco Michel
78  }
79 
80  caf.reco_numu = 1; caf.reco_nue = 0; caf.reco_nc = 0;
81  caf.muon_contained = 1; caf.muon_tracker = 0; caf.muon_ecal = 0; caf.muon_exit = 0;
82  caf.Ev_reco = caf.Elep_reco;
83 }
84 
85 // Fill reco variables for muon reconstructed in magnetized tracker
86 void recoMuonECAL( CAF &caf, params &par )
87 {
88  // range-based KE
89  double ke = caf.LepE - mmu;
90  double reco_ke = rando->Gaus( ke, ke*par.ECAL_muRes );
91  caf.Elep_reco = reco_ke + mmu;
92 
93  double true_tx = 1000.*atan(caf.LepMomX / caf.LepMomZ);
94  double true_ty = 1000.*atan(caf.LepMomY / caf.LepMomZ);
95  double evalTsmear = tsmear->Eval(caf.Elep_reco - mmu);
96  if( evalTsmear < 0. ) evalTsmear = 0.;
97  double reco_tx = true_tx + rando->Gaus(0., evalTsmear/sqrt(2.));
98  double reco_ty = true_ty + rando->Gaus(0., evalTsmear/sqrt(2.));
99  caf.theta_reco = 0.001*sqrt( reco_tx*reco_tx + reco_ty*reco_ty );
100 
101  // assume perfect charge reconstruction -- these are fairly soft and should curve a lot in short distance
102  caf.reco_q = (caf.LepPDG > 0 ? -1 : 1);
103 
104  // assume always muon for ecal-matched
105  caf.reco_numu = 1; caf.reco_nue = 0; caf.reco_nc = 0;
106  caf.muon_contained = 0; caf.muon_tracker = 0; caf.muon_ecal = 1; caf.muon_exit = 0;
107  caf.Ev_reco = caf.Elep_reco;
108 }
109 
110 // Fill reco variables for true electron
111 void recoElectron( CAF &caf, params &par )
112 {
113  caf.reco_q = 0; // never know charge
114  caf.reco_numu = 0;
115  caf.muon_contained = 0; caf.muon_tracker = 1; caf.muon_ecal = 0; caf.muon_exit = 0;
116 
117  // fake efficiency...threshold of 300 MeV, eff rising to 100% by 700 MeV
118  if( rando->Rndm() > (caf.LepE-0.3)*2.5 ) { // reco as NC
119  caf.Elep_reco = 0.;
120  caf.reco_nue = 0; caf.reco_nc = 1;
121  caf.Ev_reco = caf.LepE; // include electron energy in Ev anyway, since it won't show up in reco hadronic energy
122  } else { // reco as CC
123  caf.Elep_reco = rando->Gaus( caf.LepE, caf.LepE*(par.em_const + par.em_sqrtE/sqrt(caf.LepE)) );
124  caf.reco_nue = 1; caf.reco_nc = 0;
125  caf.Ev_reco = caf.Elep_reco;
126  }
127 
128  double true_tx = 1000.*atan(caf.LepMomX / caf.LepMomZ);
129  double true_ty = 1000.*atan(caf.LepMomY / caf.LepMomZ);
130  double evalTsmear = 3. + tsmear->Eval(caf.Elep_reco - mmu);
131  if( evalTsmear < 0. ) evalTsmear = 0.;
132  double reco_tx = true_tx + rando->Gaus(0., evalTsmear/sqrt(2.));
133  double reco_ty = true_ty + rando->Gaus(0., evalTsmear/sqrt(2.));
134  caf.theta_reco = 0.001*sqrt( reco_tx*reco_tx + reco_ty*reco_ty );
135 
136 }
137 
138 void decayPi0( TLorentzVector pi0, TVector3 &gamma1, TVector3 &gamma2 )
139 {
140  double e = pi0.E();
141  double mp = 134.9766; // pi0 mass
142 
143  double beta = sqrt( 1. - (mp*mp)/(e*e) ); // velocity of pi0
144  double theta = 3.1416*rando->Rndm(); // theta of gamma1 w.r.t. pi0 direction
145  double phi = 2.*3.1416*rando->Rndm(); // phi of gamma1 w.r.t. pi0 direction
146 
147  double p = mp/2.; // photon momentum in pi0 rest frame
148  TLorentzVector g1( 0., 0., p, p ); // pre-rotation photon 1
149  TLorentzVector g2( 0., 0., -p, p ); // pre-rotation photon 2 is opposite
150 
151  // rotate to the random decay axis in pi0 rest frame. choice of rotation about x instead of y is arbitrary
152  g1.RotateX( theta );
153  g2.RotateX( theta );
154  g1.RotateZ( phi );
155  g2.RotateZ( phi );
156 
157  // boost to lab frame with pi0 velocity. pi0 direction is z axis for this
158  g1.Boost( 0., 0., beta );
159  g2.Boost( 0., 0., beta );
160 
161  // make gamma1 the more energetic one
162  if( g1.E() > g2.E() ) {
163  gamma1 = g1.Vect();
164  gamma2 = g2.Vect();
165  } else {
166  gamma1 = g2.Vect();
167  gamma2 = g1.Vect();
168  }
169 
170  // rotate from frame where pi0 is z' direction into neutrino frame
171  TVector3 pi0dir = pi0.Vect().Unit(); // actually w.r.t. neutrino direction
172  gamma1.RotateUz( pi0dir );
173  gamma2.RotateUz( pi0dir );
174 }
175 
176 // main loop function
177 void loop( CAF &caf, params &par, TTree * tree, TTree * gtree, std::string fhicl_filename )
178 {
179  // read in dumpTree output file
180  int ievt, lepPdg, muonReco, nFS;
181  float lepKE, muGArLen, muECalLen, hadTot, hadCollar;
182  float hadP, hadN, hadPip, hadPim, hadPi0, hadOther;
183  float p3lep[3], vtx[3], muonExitPt[3], muonExitMom[3];
184  int fsPdg[100];
185  float fsPx[100], fsPy[100], fsPz[100], fsE[100], fsTrkLen[100], fsTrkLenPerp[100];
186  tree->SetBranchAddress( "ievt", &ievt );
187  tree->SetBranchAddress( "lepPdg", &lepPdg );
188  tree->SetBranchAddress( "muonReco", &muonReco );
189  tree->SetBranchAddress( "lepKE", &lepKE );
190  tree->SetBranchAddress( "muGArLen", &muGArLen );
191  tree->SetBranchAddress( "muECalLen", &muECalLen );
192  tree->SetBranchAddress( "hadTot", &hadTot );
193  tree->SetBranchAddress( "hadCollar", &hadCollar );
194  tree->SetBranchAddress( "hadP", &hadP );
195  tree->SetBranchAddress( "hadN", &hadN );
196  tree->SetBranchAddress( "hadPip", &hadPip );
197  tree->SetBranchAddress( "hadPim", &hadPim );
198  tree->SetBranchAddress( "hadPi0", &hadPi0 );
199  tree->SetBranchAddress( "hadOther", &hadOther );
200  tree->SetBranchAddress( "p3lep", p3lep );
201  tree->SetBranchAddress( "vtx", vtx );
202  tree->SetBranchAddress( "muonExitPt", muonExitPt );
203  tree->SetBranchAddress( "muonExitMom", muonExitMom );
204  tree->SetBranchAddress( "lepDeath", &caf.muon_endpoint );
205  tree->SetBranchAddress( "muon_endVolName", &caf.muon_endVolName );
206  tree->SetBranchAddress( "nFS", &nFS );
207  tree->SetBranchAddress( "fsPdg", fsPdg );
208  tree->SetBranchAddress( "fsPx", fsPx );
209  tree->SetBranchAddress( "fsPy", fsPy );
210  tree->SetBranchAddress( "fsPz", fsPz );
211  tree->SetBranchAddress( "fsE", fsE );
212  tree->SetBranchAddress( "fsTrkLen", fsTrkLen );
213  tree->SetBranchAddress( "fsTrkLenPerp", fsTrkLenPerp );
214 
215  tree->SetBranchAddress( "geoEffThrowResults", &caf.geoEffThrowResults );
216 
217  // DUNE reweight getter
218  nusyst::response_helper rh( fhicl_filename );
219  // Get list of variations, and make CAF branch for each one
220  std::vector<unsigned int> parIds = rh.GetParameters();
221  for( unsigned int i = 0; i < parIds.size(); ++i ) {
222  systtools::SystParamHeader head = rh.GetHeader(parIds[i]);
223  printf( "Adding reweight branch %u for %s with %lu shifts\n", parIds[i], head.prettyName.c_str(), head.paramVariations.size() );
224  bool is_wgt = head.isWeightSystematicVariation;
225  std::string wgt_var = ( is_wgt ? "wgt" : "var" );
226  caf.addRWbranch( parIds[i], head.prettyName, wgt_var, head.paramVariations );
227  caf.iswgt[parIds[i]] = is_wgt;
228  }
229 
230  caf.pot = gtree->GetWeight();
231  gtree->SetBranchAddress( "gmcrec", &caf.mcrec );
232 
233  // Main event loop
234  int N = tree->GetEntries();
235  for( int ii = par.first; ii < N; ++ii ) {
236 
237  tree->GetEntry(ii);
238  if( ii % 100 == 0 ) printf( "Event %d of %d...\n", ii, N );
239 
240  caf.setToBS();
241 
242  // set defaults
243  for( int j = 0; j < 100; ++j ) {
244  caf.nwgt[j] = 7;
245  caf.cvwgt[j] = ( caf.iswgt[j] ? 1. : 0. );
246  for( unsigned int k = 0; k < 100; ++k ) {
247  caf.wgt[j][k] = ( caf.iswgt[j] ? 1. : 0. );
248  }
249  }
250 
251  caf.vtx_x = vtx[0];
252  caf.vtx_y = vtx[1];
253  caf.vtx_z = vtx[2];
254  caf.det_x = -100.*par.OA_xcoord;
255 
256  // configuration variables in CAF file; we don't use mvaresult so just set it to zero
257  caf.run = par.run;
258  caf.subrun = par.subrun;
259  caf.event = ii;
260  caf.isFD = 0;
261  caf.isFHC = par.fhc;
262 
263  // get GENIE event record
264  gtree->GetEntry( ievt );
265  genie::EventRecord * event = caf.mcrec->event;
266  genie::Interaction * in = event->Summary();
267 
268  // Get truth stuff out of GENIE ghep record
269  caf.neutrinoPDG = in->InitState().ProbePdg();
270  caf.neutrinoPDGunosc = in->InitState().ProbePdg(); // fill this for similarity with FD, but no oscillations
271  caf.mode = in->ProcInfo().ScatteringTypeId();
272  caf.Ev = in->InitState().ProbeE(genie::kRfLab);
273  caf.LepPDG = in->FSPrimLeptonPdg();
274  caf.isCC = (abs(caf.LepPDG) == 13 || abs(caf.LepPDG) == 11);
275 
276  TLorentzVector lepP4;
277  TLorentzVector nuP4nuc = *(in->InitState().GetProbeP4(genie::kRfHitNucRest));
278  TLorentzVector nuP4 = *(in->InitState().GetProbeP4(genie::kRfLab));
279 
280  caf.nP = 0;
281  caf.nN = 0;
282  caf.nipip = 0;
283  caf.nipim = 0;
284  caf.nipi0 = 0;
285  caf.nikp = 0;
286  caf.nikm = 0;
287  caf.nik0 = 0;
288  caf.niem = 0;
289  caf.niother = 0;
290  caf.nNucleus = 0;
291  caf.nUNKNOWN = 0; // there is an "other" category so this never gets used
292  caf.eP = 0.;
293  caf.eN = 0.;
294  caf.ePip = 0.;
295  caf.ePim = 0.;
296  caf.ePi0 = 0.;
297  caf.eOther = 0.;
298  caf.eRecoP = 0.;
299  caf.eRecoN = 0.;
300  caf.eRecoPip = 0.;
301  caf.eRecoPim = 0.;
302  caf.eRecoPi0 = 0.;
303  caf.eOther = 0.;
304  for( int i = 0; i < nFS; ++i ) {
305  double ke = 0.001*(fsE[i] - sqrt(fsE[i]*fsE[i] - fsPx[i]*fsPx[i] - fsPy[i]*fsPy[i] - fsPz[i]*fsPz[i]));
306  if( fsPdg[i] == caf.LepPDG ) {
307  lepP4.SetPxPyPzE( fsPx[i]*0.001, fsPy[i]*0.001, fsPz[i]*0.001, fsE[i]*0.001 );
308  caf.LepE = fsE[i]*0.001;
309  }
310  else if( fsPdg[i] == 2212 ) {caf.nP++; caf.eP += ke;}
311  else if( fsPdg[i] == 2112 ) {caf.nN++; caf.eN += ke;}
312  else if( fsPdg[i] == 211 ) {caf.nipip++; caf.ePip += ke;}
313  else if( fsPdg[i] == -211 ) {caf.nipim++; caf.ePim += ke;}
314  else if( fsPdg[i] == 111 ) {caf.nipi0++; caf.ePi0 += ke;}
315  else if( fsPdg[i] == 321 ) {caf.nikp++; caf.eOther += ke;}
316  else if( fsPdg[i] == -321 ) {caf.nikm++; caf.eOther += ke;}
317  else if( fsPdg[i] == 311 || fsPdg[i] == -311 || fsPdg[i] == 130 || fsPdg[i] == 310 ) {caf.nik0++; caf.eOther += ke;}
318  else if( fsPdg[i] == 22 ) {caf.niem++; caf.eOther += ke;}
319  else if( fsPdg[i] > 1000000000 ) caf.nNucleus++;
320  else {caf.niother++; caf.eOther += ke;}
321  }
322 
323  // true 4-momentum transfer
324  TLorentzVector q = nuP4-lepP4;
325 
326  // Q2, W, x, y frequently do not get filled in GENIE Kinematics object, so calculate manually
327  caf.Q2 = -q.Mag2();
328  caf.W = sqrt(0.939*0.939 + 2.*q.E()*0.939 + q.Mag2()); // "Wexp"
329  caf.X = -q.Mag2()/(2*0.939*q.E());
330  caf.Y = q.E()/caf.Ev;
331 
332  caf.theta_reco = -1.; // default value
333 
334  caf.NuMomX = nuP4.X();
335  caf.NuMomY = nuP4.Y();
336  caf.NuMomZ = nuP4.Z();
337  caf.LepMomX = lepP4.X();
338  caf.LepMomY = lepP4.Y();
339  caf.LepMomZ = lepP4.Z();
340  caf.LepE = lepP4.E();
341  caf.LepNuAngle = nuP4.Angle( lepP4.Vect() );
342 
343  // Add DUNErw weights to the CAF
344 
345  systtools::event_unit_response_w_cv_t resp = rh.GetEventVariationAndCVResponse(*event);
346  for( systtools::event_unit_response_w_cv_t::iterator it = resp.begin(); it != resp.end(); ++it ) {
347  caf.nwgt[(*it).pid] = (*it).responses.size();
348  caf.cvwgt[(*it).pid] = (*it).CV_response;
349  for( unsigned int i = 0; i < (*it).responses.size(); ++i ) {
350  caf.wgt[(*it).pid][i] = (*it).responses[i];
351  }
352  }
353 
354  //--------------------------------------------------------------------------
355  // Parameterized reconstruction
356  //--------------------------------------------------------------------------
357  if( !par.IsGasTPC ) {
358  // Loop over final-state particles
359  double longest_mip = 0.;
360  double longest_mip_KE = 0.;
361  int longest_mip_charge = 0;
362  caf.reco_lepton_pdg = 0;
363  int electrons = 0;
364  double electron_energy = 0.;
365  int reco_electron_pdg = 0;
366  for( int i = 0; i < nFS; ++i ) {
367  int pdg = fsPdg[i];
368  double p = sqrt(fsPx[i]*fsPx[i] + fsPy[i]*fsPy[i] + fsPz[i]*fsPz[i]);
369  double KE = fsE[i] - sqrt(fsE[i]*fsE[i] - p*p);
370 
371  if( (abs(pdg) == 13 || abs(pdg) == 211) && fsTrkLen[i] > longest_mip ) {
372  longest_mip = fsTrkLen[i];
373  longest_mip_KE = KE;
374  caf.reco_lepton_pdg = pdg;
375  if( pdg == 13 || pdg == -211 ) longest_mip_charge = -1;
376  else longest_mip_charge = 1;
377  }
378 
379  // pi0 as nu_e
380  if( pdg == 111 ) {
381  TVector3 g1, g2;
382  TLorentzVector pi0( fsPx[i], fsPy[i], fsPz[i], fsE[i] );
383  decayPi0( pi0, g1, g2 );
384  double g1conv = rando->Exp( 14. ); // conversion distance
385  bool compton = (rando->Rndm() < 0.15); // dE/dX misID probability for photon
386  // if energetic gamma converts in first wire, and other gamma is either too soft or too colinear
387  if( g1conv < 2.0 && compton && (g2.Mag() < 50. || g1.Angle(g2) < 0.01) ) electrons++;
388  electron_energy = g1.Mag();
389  reco_electron_pdg = 111;
390  }
391  }
392 
393  // True CC reconstruction
394  if( abs(lepPdg) == 11 ) { // true nu_e
395  recoElectron( caf, par );
396  electrons++;
397  reco_electron_pdg = lepPdg;
398  } else if( abs(lepPdg) == 13 ) { // true nu_mu
399  if ( muGArLen > 50. ) recoMuonTracker( caf, par ); // gas TPC match
400  else if( muonReco == 1 ) recoMuonLAr( caf, par ); // LAr-contained muon, this might get updated to NC...
401  else if( muonReco == 3 && muECalLen > 5. ) recoMuonECAL( caf, par ); // ECAL-stopper
402  else { // exiting but poorly-reconstructed muon
403  caf.Elep_reco = longest_mip * 0.0022;
404  caf.reco_q = 0;
405  caf.reco_numu = 1; caf.reco_nue = 0; caf.reco_nc = 0;
406  caf.muon_contained = 0; caf.muon_tracker = 0; caf.muon_ecal = 0; caf.muon_exit = 1;
407 
408  double true_tx = 1000.*atan(caf.LepMomX / caf.LepMomZ);
409  double true_ty = 1000.*atan(caf.LepMomY / caf.LepMomZ);
410  double evalTsmear = tsmear->Eval(caf.Elep_reco - mmu);
411  if( evalTsmear < 0. ) evalTsmear = 0.;
412  double reco_tx = true_tx + rando->Gaus(0., evalTsmear/sqrt(2.));
413  double reco_ty = true_ty + rando->Gaus(0., evalTsmear/sqrt(2.));
414  caf.theta_reco = 0.001*sqrt( reco_tx*reco_tx + reco_ty*reco_ty );
415  }
416  } else { // NC -- set PID variables, will get updated later if fake CC
417  caf.Elep_reco = 0.;
418  caf.reco_q = 0;
419  caf.reco_numu = 0; caf.reco_nue = 0; caf.reco_nc = 1;
420  caf.muon_contained = 0; caf.muon_tracker = 0; caf.muon_ecal = 0; caf.muon_exit = 0;
421  }
422 
423  // CC/NC confusion
424  if( electrons == 1 && muonReco <= 1 ) { // NC or numuCC reco as nueCC
425  caf.Elep_reco = electron_energy*0.001;
426  caf.reco_q = 0;
427  caf.reco_numu = 0; caf.reco_nue = 1; caf.reco_nc = 0;
428  caf.muon_contained = 0; caf.muon_tracker = 0; caf.muon_ecal = 0; caf.muon_exit = 0;
429  caf.reco_lepton_pdg = reco_electron_pdg;
430  } else if( muonReco <= 1 && !(abs(lepPdg) == 11 && caf.Elep_reco > 0.) && (longest_mip < par.CC_trk_length || longest_mip_KE/longest_mip > 3.) ) {
431  // reco as NC
432  caf.Elep_reco = 0.;
433  caf.reco_q = 0;
434  caf.reco_numu = 0; caf.reco_nue = 0; caf.reco_nc = 1;
435  caf.muon_contained = 0; caf.muon_tracker = 0; caf.muon_ecal = 0; caf.muon_exit = 0;
436  caf.reco_lepton_pdg = 0;
437  } else if( (abs(lepPdg) == 12 || abs(lepPdg) == 14) && longest_mip > par.CC_trk_length && longest_mip_KE/longest_mip < 3. ) { // true NC reco as CC numu
438  caf.Elep_reco = longest_mip_KE*0.001 + mmu;
439  if( par.fhc ) caf.reco_q = -1;
440  else {
441  double michel = rando->Rndm();
442  if( longest_mip_charge == 1 && michel < par.michelEff ) caf.reco_q = 1; // correct mu+
443  else if( michel < par.michelEff*0.25 ) caf.reco_q = 1; // incorrect mu-
444  else caf.reco_q = -1; // no reco Michel
445  }
446  caf.reco_numu = 1; caf.reco_nue = 0; caf.reco_nc = 0;
447  caf.muon_contained = 1; caf.muon_tracker = 0; caf.muon_ecal = 0; caf.muon_exit = 0;
448  }
449 
450  // Hadronic energy calorimetrically
451  caf.Ev_reco = caf.Elep_reco + hadTot*0.001;
452  caf.Ehad_veto = hadCollar;
453  caf.eRecoP = hadP*0.001;
454  caf.eRecoN = hadN*0.001;
455  caf.eRecoPip = hadPip*0.001;
456  caf.eRecoPim = hadPim*0.001;
457  caf.eRecoPi0 = hadPi0*0.001;
458  caf.eRecoOther = hadOther*0.001;
459 
460  caf.pileup_energy = 0.;
461  if( rando->Rndm() < par.pileup_frac ) caf.pileup_energy = rando->Rndm() * par.pileup_max;
462  caf.Ev_reco += caf.pileup_energy;
463  } else {
464  // gas TPC: FS particle loop look for long enough tracks and smear momenta
465  caf.Ev_reco = 0.;
466  caf.nFSP = nFS;
467  for( int i = 0; i < nFS; ++i ) {
468  double ptrue = 0.001*sqrt(fsPx[i]*fsPx[i] + fsPy[i]*fsPy[i] + fsPz[i]*fsPz[i]);
469  double mass = 0.001*sqrt(fsE[i]*fsE[i] - fsPx[i]*fsPx[i] - fsPy[i]*fsPy[i] - fsPz[i]*fsPz[i]);
470  caf.pdg[i] = fsPdg[i];
471  caf.ptrue[i] = ptrue;
472  caf.trkLen[i] = fsTrkLen[i];
473  caf.trkLenPerp[i] = fsTrkLenPerp[i];
474  // track length cut 6cm according to T Junk
475  if( fsTrkLen[i] > 0. && fsPdg[i] != 2112 ) { // basically select charged particles; somehow neutrons ocasionally get nonzero track length
476  double pT = 0.001*sqrt(fsPy[i]*fsPy[i] + fsPz[i]*fsPz[i]); // transverse to B field, in GeV
477  double nHits = fsTrkLen[i] / par.gastpc_padPitch; // doesn't matter if not integer as only used in eq
478  // Gluckstern formula, sigmapT/pT, with sigmaX and L in meters
479  double fracSig_meas = sqrt(720./(nHits+4)) * (0.01*par.gastpc_padPitch/sqrt(12.)) * pT / (0.3 * par.gastpc_B * 0.0001 * fsTrkLenPerp[i]*fsTrkLenPerp[i]);
480  // multiple scattering term
481  double fracSig_MCS = 0.052 / (par.gastpc_B * sqrt(par.gastpc_X0*fsTrkLenPerp[i]*0.0001));
482 
483  double sigmaP = ptrue * sqrt( fracSig_meas*fracSig_meas + fracSig_MCS*fracSig_MCS );
484  double preco = rando->Gaus( ptrue, sigmaP );
485  double ereco = sqrt( preco*preco + mass*mass ) - mass; // kinetic energy
486  if( abs(fsPdg[i]) == 211 ) ereco += mass; // add pion mass
487  else if( fsPdg[i] == 2212 && preco > 1.5 ) ereco += 0.1395; // mistake pion mass for high-energy proton
488  caf.partEvReco[i] = ereco;
489 
490  // threshold cut
491  if( fsTrkLen[i] > par.gastpc_len ) {
492  caf.Ev_reco += ereco;
493  if( fsPdg[i] == 211 || (fsPdg[i] == 2212 && preco > 1.5) ) caf.gastpc_pi_pl_mult++;
494  else if( fsPdg[i] == -211 ) caf.gastpc_pi_min_mult++;
495  }
496 
497  if( (fsPdg[i] == 13 || fsPdg[i] == -13) && fsTrkLen[i] > 100. ) { // muon, don't really care about nu_e CC for now
498  caf.partEvReco[i] += mass;
499  caf.Elep_reco = sqrt(preco*preco + mass*mass);
500  // angle reconstruction
501  double true_tx = 1000.*atan(caf.LepMomX / caf.LepMomZ);
502  double true_ty = 1000.*atan(caf.LepMomY / caf.LepMomZ);
503  double evalTsmear = tsmear->Eval(caf.Elep_reco - mmu);
504  if( evalTsmear < 0. ) evalTsmear = 0.;
505  double reco_tx = true_tx + rando->Gaus(0., evalTsmear/sqrt(2.));
506  double reco_ty = true_ty + rando->Gaus(0., evalTsmear/sqrt(2.));
507  caf.theta_reco = 0.001*sqrt( reco_tx*reco_tx + reco_ty*reco_ty );
508  // assume perfect charge reconstruction
509  caf.reco_q = (fsPdg[i] > 0 ? -1 : 1);
510  caf.reco_numu = 1; caf.reco_nue = 0; caf.reco_nc = 0;
511  caf.muon_tracker = 1;
512  }
513  } else if( fsPdg[i] == 111 || fsPdg[i] == 22 ) {
514  double ereco = 0.001 * rando->Gaus( fsE[i], 0.1*fsE[i] );
515  caf.partEvReco[i] = ereco;
516  caf.Ev_reco += ereco;
517  }
518  }
519  }
520 
521  caf.fill();
522  }
523 
524  // set POT
525  caf.meta_run = par.run;
526  caf.meta_subrun = par.subrun;
527 
528 }
529 
530 int main( int argc, char const *argv[] )
531 {
532 
533  if( (argc == 2) && ((std::string("--help") == argv[1]) || (std::string("-h") == argv[1])) ) {
534  std::cout << "Help yourself by looking at the source code to see what the options are." << std::endl;
535  return 0;
536  }
537 
538  // Need this to store event-by-event geometric efficiency
539  gInterpreter->GenerateDictionary("vector<vector<vector<uint64_t> > >", "vector");
540 
541  // get command line options
542  std::string gfile;
545  std::string fhicl_filename;
546 
547  // Make parameter object and set defaults
548  params par;
549  par.IsGasTPC = false;
550  par.OA_xcoord = 0.; // on-axis by default
551  par.fhc = true;
552  par.grid = false;
553  par.seed = 7; // a very random number
554  par.run = 1; // CAFAna doesn't like run number 0
555  par.subrun = 0;
556  par.first = 0;
557  par.trk_muRes = 0.02; // fractional muon energy resolution of HP GAr TPC
558  par.LAr_muRes = 0.05; // fractional muon energy resolution of muons contained in LAr
559  par.ECAL_muRes = 0.1; // fractional muon energy resolution of muons ending in ECAL
560  par.em_const = 0.03; // EM energy resolution constant term: A + B/sqrt(E) (GeV)
561  par.em_sqrtE = 0.1; // EM energy resolution 1/sqrt(E) term: A + B/sqrt(E) (GeV)
562  par.michelEff = 0.75; // Michel finder efficiency
563  par.CC_trk_length = 100.; // minimum track length for CC in cm
564  par.pileup_frac = 0.1; // fraction of events with non-zero pile-up
565  par.pileup_max = 0.5; // GeV
566  par.gastpc_len = 6.; // track length cut in cm
567  par.gastpc_B = 0.4; // B field strength in Tesla
568  par.gastpc_padPitch = 0.1; // 1 mm. Actual pad pitch varies, which is going to be impossible to implement
569  par.gastpc_X0 = 1300.; // cm = 13m radiation length
570 
571  int i = 0;
572  while( i < argc ) {
573  if( argv[i] == std::string("--infile") ) {
574  infile = argv[i+1];
575  i += 2;
576  } else if( argv[i] == std::string("--gfile") ) {
577  gfile = argv[i+1];
578  i += 2;
579  } else if( argv[i] == std::string("--outfile") ) {
580  outfile = argv[i+1];
581  i += 2;
582  } else if( argv[i] == std::string("--fhicl") ) {
583  fhicl_filename = argv[i+1];
584  i += 2;
585  } else if( argv[i] == std::string("--seed") ) {
586  par.seed = atoi(argv[i+1]);
587  par.run = par.seed;
588  i += 2;
589  } else if( argv[i] == std::string("--oa") ) {
590  par.OA_xcoord = atof(argv[i+1]);
591  i += 2;
592  } else if( argv[i] == std::string("--rhc") ) {
593  par.fhc = false;
594  i += 1;
595  } else if( argv[i] == std::string("--gastpc") ) {
596  par.IsGasTPC = true;
597  i += 1;
598  } else i += 1; // look for next thing
599  }
600 
601  rando = new TRandom3( par.seed );
602  CAF caf( outfile, par.IsGasTPC );
603 
604  // LAr driven smearing, maybe we want to change for gas?
605  tsmear = new TF1( "tsmear", "0.162 + 3.407*pow(x,-1.) + 3.129*pow(x,-0.5)", 0., 999.9 );
606 
607  TFile * tf = new TFile( infile.c_str() );
608  TTree * tree = (TTree*) tf->Get( "tree" );
609 
610  TFile * gf = new TFile( gfile.c_str() );
611  TTree * gtree = (TTree*) gf->Get( "gtree" );
612 
613  loop( caf, par, tree, gtree, fhicl_filename );
614 
615  caf.version = 4;
616  printf( "Run %d POT %g\n", caf.meta_run, caf.pot );
617  caf.fillPOT();
618 
619  // Copy geometric efficiency throws TTree to CAF file
620  std::cout << "Copying geometric efficiency throws TTree to output file" << std::endl;
621  TTree *tGeoEfficiencyThrowsOut = (TTree*) tf->Get("geoEffThrows");
622  caf.cafFile->cd();
623  tGeoEfficiencyThrowsOut->CloneTree()->Write();
624 
625  std::cout << "Writing CAF" << std::endl;
626  caf.write();
627 
628 }
intermediate_table::iterator iterator
void addRWbranch(int parId, std::string name, std::string wgt_var, std::vector< double > &vars)
int first
Definition: makeNUSYS.cxx:17
double ePip
Definition: CAF.h:31
double vtx_y
Definition: CAF.h:35
double eOther
Definition: CAF.h:31
int meta_subrun
Definition: NUSYS.h:42
double beta(double KE, const simb::MCParticle *part)
int neutrinoPDG
Definition: CAF.h:27
int muon_contained
Definition: CAF.h:41
void recoMuonLAr(CAF &caf, params &par)
Definition: makeCAF.cxx:55
int nNucleus
Definition: CAF.h:30
double michelEff
Definition: makeCAF.cxx:23
double LepNuAngle
Definition: CAF.h:28
double NuMomY
Definition: CAF.h:28
double eRecoP
Definition: CAF.h:32
int nP
Definition: CAF.h:30
int LepPDG
Definition: CAF.h:27
int pdg[100]
Definition: CAF.h:50
double eP
Definition: CAF.h:31
int seed
Definition: makeNUSYS.cxx:17
std::string string
Definition: nybbler.cc:12
bool fhc
Definition: makeNUSYS.cxx:16
bool iswgt[100]
Definition: NUSYS.h:35
int niem
Definition: CAF.h:30
double partEvReco[100]
Definition: CAF.h:51
void fill()
double Ev_reco
Definition: CAF.h:39
double eRecoPip
Definition: CAF.h:32
genie::NtpMCEventRecord * mcrec
Definition: NUSYS.h:38
double cvwgt[100]
Definition: NUSYS.h:33
std::vector< double > trkLen
Definition: CAF.h:78
std::vector< double > trkLenPerp
Definition: CAF.h:78
double em_const
Definition: makeCAF.cxx:22
void recoElectron(CAF &caf, params &par)
Definition: makeCAF.cxx:111
double gastpc_len
Definition: makeCAF.cxx:26
int nN
Definition: CAF.h:30
int nFSP
Definition: CAF.h:49
double OA_xcoord
Definition: makeNUSYS.cxx:15
int reco_q
Definition: CAF.h:40
double LepMomY
Definition: CAF.h:28
const double mmu
Definition: makeCAF.cxx:13
void fillPOT()
int nikp
Definition: CAF.h:30
double gastpc_padPitch
Definition: makeCAF.cxx:26
Definition: tf_graph.h:23
int subrun
Definition: NUSYS.h:25
int isFD
Definition: NUSYS.h:23
int main(int argc, char const *argv[])
Definition: makeCAF.cxx:530
int run
Definition: NUSYS.h:25
double NuMomZ
Definition: CAF.h:28
int gastpc_pi_pl_mult
Definition: CAF.h:48
double W
Definition: CAF.h:28
double LepE
Definition: CAF.h:28
int FSPrimLeptonPdg(void) const
final state primary lepton pdg
double det_x
Definition: CAF.h:36
int isCC
Definition: CAF.h:27
double pileup_max
Definition: makeCAF.cxx:25
double ePim
Definition: CAF.h:31
int subrun
Definition: makeNUSYS.cxx:17
int isFHC
Definition: NUSYS.h:23
int nfiles
Definition: makeNUSYS.cxx:17
int reco_lepton_pdg
Definition: CAF.h:41
void write()
int meta_run
Definition: NUSYS.h:42
Summary information for an interaction.
Definition: Interaction.h:56
T abs(T value)
TFile * cafFile
The output TFile pointer */.
Definition: CAF.h:53
double Q2
Definition: CAF.h:28
int neutrinoPDGunosc
Definition: CAF.h:27
double pileup_frac
Definition: makeCAF.cxx:25
const double e
std::string * muon_endVolName
Definition: CAF.h:43
string infile
double gastpc_X0
Definition: makeCAF.cxx:26
double gastpc_B
Definition: makeCAF.cxx:26
double vtx_x
Definition: CAF.h:35
int nipim
Definition: CAF.h:30
void recoMuonECAL(CAF &caf, params &par)
Definition: makeCAF.cxx:86
int ProbePdg(void) const
Definition: InitialState.h:64
double NuMomX
Definition: CAF.h:28
int run
Definition: makeNUSYS.cxx:17
p
Definition: test.py:223
double LepMomZ
Definition: CAF.h:28
TRandom3 * rando
Definition: makeCAF.cxx:12
int event
Definition: NUSYS.h:25
int nikm
Definition: CAF.h:30
double ECAL_muRes
Definition: makeCAF.cxx:21
int muon_exit
Definition: CAF.h:41
double eN
Definition: CAF.h:31
bool grid
Definition: makeNUSYS.cxx:16
ScatteringType_t ScatteringTypeId(void) const
double Y
Definition: CAF.h:28
int version
Definition: NUSYS.h:43
double theta_reco
Definition: CAF.h:39
tGeoEfficiencyThrowsOut
Definition: dumpTree.py:470
Definition: CAF.h:12
double em_sqrtE
Definition: makeCAF.cxx:22
bool IsGasTPC
Definition: makeCAF.cxx:19
double eRecoOther
Definition: CAF.h:32
double CC_trk_length
Definition: makeCAF.cxx:24
double Ehad_veto
Definition: CAF.h:44
int n
Definition: makeNUSYS.cxx:17
double X
Definition: CAF.h:28
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:37
int nwgt[100]
Definition: NUSYS.h:32
int reco_nue
Definition: CAF.h:40
int gastpc_pi_min_mult
Definition: CAF.h:48
double ptrue[100]
Definition: CAF.h:51
double Ev
Definition: CAF.h:28
double eRecoPim
Definition: CAF.h:32
double Elep_reco
Definition: CAF.h:39
genFinder * gf
int reco_numu
Definition: CAF.h:40
int reco_nc
Definition: CAF.h:40
int muon_ecal
Definition: CAF.h:41
double pileup_energy
Definition: CAF.h:45
void loop(CAF &caf, params &par, TTree *tree, TTree *gtree, std::string fhicl_filename)
Definition: makeCAF.cxx:177
int nipi0
Definition: CAF.h:30
int niother
Definition: CAF.h:30
const InitialState & InitState(void) const
Definition: Interaction.h:69
double trk_muRes
Definition: makeCAF.cxx:21
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double eRecoN
Definition: CAF.h:32
std::vector< std::vector< std::vector< uint64_t > > > * geoEffThrowResults
Definition: CAF.h:64
std::vector< int > mode
Definition: CAF.h:69
void setToBS()
int nik0
Definition: CAF.h:30
Common Analysis Files.
Definition: SRGAr.h:13
double vtx_z
Definition: CAF.h:35
double ePi0
Definition: CAF.h:31
double wgt[100][100]
Definition: NUSYS.h:34
void decayPi0(TLorentzVector pi0, TVector3 &gamma1, TVector3 &gamma2)
Definition: makeCAF.cxx:138
int muon_tracker
Definition: CAF.h:41
double LAr_muRes
Definition: makeCAF.cxx:21
float muon_endpoint[3]
Definition: CAF.h:42
TF1 * tsmear
Definition: makeCAF.cxx:14
double eRecoPi0
Definition: CAF.h:32
double ProbeE(RefFrame_t rf) const
int nipip
Definition: CAF.h:30
EventRecord * event
event
int nUNKNOWN
Definition: CAF.h:30
QTextStream & endl(QTextStream &s)
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
double LepMomX
Definition: CAF.h:28
Event finding and building.
double pot
Definition: NUSYS.h:41
void recoMuonTracker(CAF &caf, params &par)
Definition: makeCAF.cxx:30