777 std::vector<const anab::T0*> T0s;
786 std::cout <<
"No T0s found" <<
std::endl;
803 fTimeStamp = ts2.AsDouble();
807 if(!
evt.isRealData()){
808 for(
int k=0;
k < 6;
k++)
812 for(
int k=0;
k < 6;
k++)
816 bool beamTriggerEvent =
false;
818 auto mcTruths =
evt.getValidHandle<std::vector<simb::MCTruth>>(
fGeneratorTag);
819 if(!
evt.isRealData()){
821 auto beamHandle =
evt.getValidHandle<std::vector<beam::ProtoDUNEBeamEvent>>(
"generator");
822 std::vector<art::Ptr<beam::ProtoDUNEBeamEvent>> beamVec;
823 if( beamHandle.isValid()){
830 std::cout <<
"Number of reconstructed beam momenta from spec: " << momenta.size() <<
std::endl;
832 if( momenta.size() > 0 ) std::cout <<
"Measured Beam Momentum from spec: " << momenta.at(0) <<
std::endl;
835 for (
size_t i = 0; i<momenta.size(); ++i){
841 std::cout<<
"beam trk size:"<<btracks.size()<<
std::endl;
842 for (
size_t i = 0; i<btracks.size(); ++i){
843 std::cout<<
"beamPosx/beamPosy/beamPosz:"<<btracks[i].End().X()<<
"/"<<btracks[i].End().Y()<<
"/"<<btracks[i].End().Z()<<
std::endl;
844 std::cout<<
"beamDirx/beamDiry/beamDirz:"<<btracks[i].StartDirection().X()<<
"/"<<btracks[i].StartDirection().Y()<<
"/"<<btracks[i].StartDirection().Z()<<
std::endl;
862 auto mcTruths =
evt.getValidHandle<std::vector<simb::MCTruth>>(
fGeneratorTag);
879 if(geantGoodParticle != 0x0){
880 std::cout <<
"Found GEANT particle corresponding to the good particle with pdg = " << geantGoodParticle->PdgCode()
881 <<
" , track id = " << geantGoodParticle->TrackId()
882 <<
" , Vx/Vy/Vz = " << geantGoodParticle->Vx() <<
"/"<< geantGoodParticle->Vy() <<
"/" << geantGoodParticle->Vz()
888 beamTriggerEvent =
true;
903 for(
size_t i_s=0; i_s < geantGoodParticle->NumberTrajectoryPoints(); i_s++){
904 beamtrk_x.push_back(geantGoodParticle->Position(i_s).X());
905 beamtrk_y.push_back(geantGoodParticle->Position(i_s).Y());
906 beamtrk_z.push_back(geantGoodParticle->Position(i_s).Z());
908 beamtrk_Px.push_back(geantGoodParticle->Momentum(i_s).X());
909 beamtrk_Py.push_back(geantGoodParticle->Momentum(i_s).Y());
910 beamtrk_Pz.push_back(geantGoodParticle->Momentum(i_s).Z());
912 beamtrk_Eng.push_back(geantGoodParticle->Momentum(i_s).E()-geantGoodParticle->Mass());
916 int key_reach_tpc=-99;
917 for (
size_t kk=0; kk<
beamtrk_z.size(); ++kk) {
919 if ((zpos_beam+0.49375)<0.01&&(zpos_beam+0.49375)>-0.01) {
925 if (key_reach_tpc!=-99) {
936 std::cout<<
"This particle doesn't enter TPC!!!"<<
std::endl;
943 beamid = geantGoodParticle->TrackId();
982 std::cout <<
"Doing reweight" <<
std::endl;
988 if (created && theTraj.GetNSteps()) {
994 std::vector<double> weights_vec =
MultiRW.GetWeightFromAll1DThrows(
997 weights_vec.begin(), weights_vec.end());
999 for (
size_t i = 0; i <
ParSet.size(); ++i) {
1000 std::pair<double, double> pm_weights =
1001 MultiRW.GetPlusMinusSigmaParWeight(theTraj, i);
1010 for (
size_t i = 0; i < 300; ++i) {
1011 for (
size_t j = 0; j < 200; ++j) {
1012 double tmp_p1=(.9 + i*.001);
1013 double tmp_p2=(.9 + j*.001);
1015 std::vector<double> input_values = {tmp_p1,tmp_p2};
1017 bool set_values =
MultiRW.SetAllParameterValues(input_values);
1018 if (!set_values)
continue;
1023 MultiRW.GetWeightFromSetParameters(theTraj)
1045 std::vector< art::Ptr<beam::ProtoDUNEBeamEvent> > beaminfo;
1046 auto pdbeamHandle =
evt.getHandle< std::vector<beam::ProtoDUNEBeamEvent> >(
fBeamModuleLabel);
1050 for(
unsigned int i = 0; i < beaminfo.size(); ++i){
1056 if(beaminfo[i]->GetTOFChan() >= 0)
1057 ftof = beaminfo[i]->GetTOF();
1060 if(beaminfo[i]->GetBITrigger() == 1){
1066 auto &
tracks = beaminfo[i]->GetBeamTracks();
1077 auto & beammom = beaminfo[i]->GetRecoBeamMomenta();
1078 if(!beammom.empty())
1108 auto recoParticles =
evt.getValidHandle<std::vector<recob::PFParticle>>(
fPFParticleTag);
1119 std::vector<art::Ptr<recob::Track> > tracklist;
1120 std::vector<art::Ptr<recob::PFParticle> > pfplist;
1122 auto trackListHandle =
evt.getHandle< std::vector<recob::Track> >(
"pandoraTrack");
1126 auto PFPListHandle =
evt.getHandle< std::vector<recob::PFParticle> > (
"pandora");
1129 std::vector<art::Ptr<recob::Cluster>> clusterlist;
1130 auto clusterListHandle =
evt.getHandle< std::vector<recob::Cluster> >(
"pandora") ;
1133 art::FindManyP<recob::Cluster> fmcp(PFPListHandle,
evt,
"pandora");
1134 art::FindManyP<recob::Track> pftrack(PFPListHandle,
evt,
"pandoraTrack");
1136 std::cout<<
"number of pfp_particles "<<pfplist.size()<<
std::endl;
1137 std::cout<<
" size of pfParticles "<<pfParticles.size()<<
std::endl;
1138 art::FindManyP<recob::Hit, recob::TrackHitMeta> fmthm(trackListHandle,
evt,
"pandoraTrack");
1171 if(!pfT0vec.empty())
1183 if(thisTrack != 0x0) {
1186 std::cout<<
"inside the this track loop "<<
std::endl;
1187 if(mcparticle0!=0x0) {
1191 truthid=mcparticle0->
TrackId();
1247 if (hi->WireID().Plane==2) {
1251 if (hi->WireID().Plane==1) {
1255 if (hi->WireID().Plane==0) {
1289 if(calovector.size() != 3)
1290 std::cerr <<
"WARNING::Calorimetry vector size for primary is = " << calovector.size() <<
std::endl;
1293 std::vector<double> tmp_primtrk_dqdx;
1294 std::vector<double> tmp_primtrk_resrange;
1295 std::vector<double> tmp_primtrk_dedx;
1296 std::vector<double> tmp_primtrk_hitx;
1297 std::vector<double> tmp_primtrk_hity;
1298 std::vector<double> tmp_primtrk_hitz;
1299 std::vector<double> tmp_primtrk_pitch;
1300 std::vector<int> tmp_primtrk_wid;
1302 for (
auto &
calo : calovector) {
1303 if (
calo.PlaneID().Plane == 2){
1306 for (
size_t ihit = 0; ihit <
calo.dQdx().size(); ++ihit){
1307 tmp_primtrk_dqdx.push_back(
calo.dQdx()[ihit]);
1308 tmp_primtrk_resrange.push_back(
calo.ResidualRange()[ihit]);
1309 tmp_primtrk_dedx.push_back(
calo.dEdx()[ihit]);
1310 tmp_primtrk_pitch.push_back(
calo.TrkPitchVec()[ihit]);
1312 const auto &primtrk_pos=(
calo.XYZ())[ihit];
1313 tmp_primtrk_hitx.push_back(primtrk_pos.X());
1314 tmp_primtrk_hity.push_back(primtrk_pos.Y());
1315 tmp_primtrk_hitz.push_back(primtrk_pos.Z());
1317 double pos_reco[3]={primtrk_pos.X(),primtrk_pos.Y(),primtrk_pos.Z()};
1330 tmp_primtrk_wid.push_back(wireID.
Wire);
1333 tmp_primtrk_wid.push_back(-9999);
1354 std::vector<double> tmp_primtrk_hitx_true;
1355 std::vector<double> tmp_primtrk_hity_true;
1356 std::vector<double> tmp_primtrk_hitz_true;
1357 std::vector<double> tmp_primtrk_trkid_true;
1358 std::vector<double> tmp_primtrk_edept_true;
1360 std::vector<double> tmp_primtrk_true_x;
1361 std::vector<double> tmp_primtrk_true_y;
1362 std::vector<double> tmp_primtrk_true_z;
1363 std::vector<double> tmp_primtrk_true_trkid;
1364 std::vector<double> tmp_primtrk_true_edept;
1365 std::vector<int> tmp_primtrk_true_wid;
1369 std::vector<double> tmp_primtrj_true_x;
1370 std::vector<double> tmp_primtrj_true_y;
1371 std::vector<double> tmp_primtrj_true_z;
1372 std::vector<double> tmp_primtrj_true_edept;
1374 if(geantGoodParticle1!= 0x0 && geantGoodParticle1->
Process()==
"primary") {
1377 double ke_reco=
ke_ff;
1378 double ke_true=
ke_ff;
1386 if (thisTrajectoryProcessMap1.size()){
1387 for(
auto const& couple: thisTrajectoryProcessMap1){
1407 if (thisTrajectoryProcessMap1.size()) {
1408 for(
auto const& couple1: thisTrajectoryProcessMap1) {
1409 interactionX.push_back(((truetraj.
at(couple1.first)).first).X());
1410 interactionY.push_back(((truetraj.
at(couple1.first)).first).Y());
1411 interactionZ.push_back(((truetraj.
at(couple1.first)).first).Z());
1412 interactionE.push_back(((truetraj.
at(couple1.first)).first).E());
1416 double xval=((truetraj.
at(couple1.first)).first).X();
1417 double yval=((truetraj.
at(couple1.first)).first).Y();
1418 double zval=((truetraj.
at(couple1.first)).first).Z();
1419 unsigned int tpcno=1;
1420 if(xval<=0 && zval<232) tpcno=1;
1421 if(xval<=0 && zval>232 && zval<464) tpcno=5;
1422 if(xval<=0 && zval>=464) tpcno=9;
1423 if(xval>0 && zval<232) tpcno=2;
1424 if(xval>0 && zval>232 && zval<464) tpcno=6;
1425 if(xval>0 && zval>=464) tpcno=10;
1451 if ((truetraj.
KeyToProcess(couple1.second)).find(
"CoulombScat")!= std::string::npos)
continue;
1454 auto interactionPos4D = (truetraj.
at(couple1.first)).first ;
1456 if (interactionPos4D.Z() <
minZ || interactionPos4D.Z() >
maxZ )
continue;
1457 else if (interactionPos4D.X() <
minX || interactionPos4D.X() >
maxX )
continue;
1458 else if (interactionPos4D.Y() <
minY || interactionPos4D.Y() >
maxY )
continue;
1461 double interactionAngle = 999999.;
1464 auto prevInteractionPos4D = (truetraj.
at(couple1.first-1)).first ;
1465 auto prevInteractionPos3D = prevInteractionPos4D.Vect() ;
1466 auto interactionPos3D = interactionPos4D.Vect() ;
1467 auto distanceBtwPoint = interactionPos3D - prevInteractionPos3D;
1470 if (truetraj.
size() > couple1.first + 1) {
1472 auto nextInteractionPos4D = (truetraj.
at(couple1.first+1)).first ;
1473 auto nextInteractionPos3D = nextInteractionPos4D.Vect() ;
1474 auto distanceBtwPointNext = nextInteractionPos3D - interactionPos3D;
1475 interactionAngles.push_back(TMath::ACos(distanceBtwPointNext.Dot(distanceBtwPoint)/(distanceBtwPointNext.Mag()*distanceBtwPoint.Mag() ) ));
1482 double maxAngle = 0.;
1483 int numberOfTroubleDaugther = 0;
1485 const sim::ParticleList& plist=pi_serv->
ParticleList();
1487 for(
size_t dPart1=0;(dPart1<plist.size()) && (plist.begin()!=plist.end());++dPart1) {
1489 if (dPart->
Mother() != 1 )
continue;
1495 interactionAngle=999999;
1499 numberOfTroubleDaugther++;
1501 auto daughtFirstPt = ((daugthTraj[0]).first).Vect() ;
1502 auto daughtSecondPt = ((daugthTraj[1]).first).Vect() ;
1503 auto distanceBtwPointNext = daughtSecondPt - daughtFirstPt;
1504 interactionAngle = TMath::ACos(distanceBtwPointNext.Dot(distanceBtwPoint)/(distanceBtwPointNext.Mag()*distanceBtwPoint.Mag()));
1505 if ( maxAngle < interactionAngle ) maxAngle = interactionAngle;
1506 interactionAngle = maxAngle;
1508 if (!numberOfTroubleDaugther) interactionAngle = 9999999.;
1509 if (interactionAngle < 0.0001) std::cout<<
"-------------------------------------------------> \n";
1518 std::cout<<
"\n MCTrajectory of this prim. trk has:"<<truetraj.
size()<<
" hits!!"<<
std::endl;
1519 for (
size_t tt=0;
tt<truetraj.
size(); ++
tt) {
1520 tmp_primtrj_true_x.push_back(truetraj.
X(
tt));
1521 tmp_primtrj_true_y.push_back(truetraj.
Y(
tt));
1522 tmp_primtrj_true_z.push_back(truetraj.
Z(
tt));
1523 tmp_primtrj_true_edept.push_back(truetraj.
E(
tt));
1531 std::map<double, sim::IDE> orderedSimIDE;
1532 for (
auto& ide : simIDE_prim) orderedSimIDE[ide->z]= *ide;
1533 auto inTPCPoint = truetraj.
begin();
1534 auto Momentum0 = inTPCPoint->second;
1541 bool run_me_onetime=
true;
1542 double xi=0.0;
double yi=0.0;
double zi=0.0;
1545 if(tmp_primtrk_dqdx.size()!=0) {
1546 if(tmp_primtrk_hitz[0]>tmp_primtrk_hitz[tmp_primtrk_hitz.size()-1]) {
1547 std::reverse(tmp_primtrk_hitz.begin(), tmp_primtrk_hitz.end());
1548 std::reverse(tmp_primtrk_hity.begin(), tmp_primtrk_hity.end());
1549 std::reverse(tmp_primtrk_hitx.begin(), tmp_primtrk_hitx.end());
1550 std::reverse(tmp_primtrk_pitch.begin(),tmp_primtrk_pitch.end());
1551 std::reverse(tmp_primtrk_dedx.begin(), tmp_primtrk_dedx.end());
1552 std::reverse(tmp_primtrk_dqdx.begin(), tmp_primtrk_dqdx.end());
1553 std::reverse(tmp_primtrk_resrange.begin(), tmp_primtrk_resrange.end());
1557 for(
size_t idx1=0; idx1<tmp_primtrk_dqdx.size()-1; idx1++) {
1558 ke_reco-=tmp_primtrk_pitch[idx1]*tmp_primtrk_dedx[idx1];
1559 double currentDepEng = 0.;
1560 for (
auto iter= orderedSimIDE.begin(); iter!= orderedSimIDE.end(); iter++) {
1561 auto currentIde = iter->second;
1564 if(currentIde.z<
minZ || currentIde.z >
maxZ )
continue;
1565 else if (currentIde.x <
minX || currentIde.x >
maxX )
continue;
1566 else if (currentIde.y <
minY || currentIde.y >
maxY )
continue;
1569 if(cnt==0&¤tIde.trackID>=0) {
1575 tmp_primtrk_hitx_true.push_back(currentIde.x);
1576 tmp_primtrk_hity_true.push_back(currentIde.y);
1577 tmp_primtrk_hitz_true.push_back(currentIde.z);
1578 tmp_primtrk_trkid_true.push_back(currentIde.trackID);
1579 tmp_primtrk_edept_true.push_back(currentIde.energy);
1585 if (currentIde.trackID>=0) {
1586 if(cnt>0 && run_me_onetime) {
1597 tmp_primtrk_hitx_true.push_back(currentIde.x);
1598 tmp_primtrk_hity_true.push_back(currentIde.y);
1599 tmp_primtrk_hitz_true.push_back(currentIde.z);
1600 tmp_primtrk_trkid_true.push_back(currentIde.trackID);
1601 tmp_primtrk_edept_true.push_back(currentIde.energy);
1603 xi=currentIde.x; yi=currentIde.y; zi=currentIde.z;
1609 if ( currentIde.z <= tmp_primtrk_hitz[idx1])
continue;
1610 if ( currentIde.z > tmp_primtrk_hitz[idx1+1])
continue;
1611 currentDepEng += currentIde.energy;
1613 ke_true -= currentDepEng;
1614 run_me_onetime=
false;
1616 if(currentDepEng>0.001) {
1628 for (
auto iter2=orderedSimIDE.begin(); iter2!=orderedSimIDE.end(); iter2++) {
1629 auto currentIde2 = iter2->second;
1632 if(currentIde2.z<
minZ || currentIde2.z >
maxZ )
continue;
1633 else if (currentIde2.x <
minX || currentIde2.x >
maxX )
continue;
1634 else if (currentIde2.y <
minY || currentIde2.y >
maxY )
continue;
1638 tmp_primtrk_true_x.push_back(currentIde2.x);
1639 tmp_primtrk_true_y.push_back(currentIde2.y);
1640 tmp_primtrk_true_z.push_back(currentIde2.z);
1641 tmp_primtrk_true_trkid.push_back(currentIde2.trackID);
1642 tmp_primtrk_true_edept.push_back(currentIde2.energy);
1644 double pos_true[3] = {currentIde2.x, currentIde2.y, currentIde2.z};
1656 tmp_primtrk_true_wid.push_back(wireID.
Wire);
1659 tmp_primtrk_true_wid.push_back(-9999);
1731 tmp_primtrk_dqdx.clear();
1732 tmp_primtrk_resrange.clear();
1733 tmp_primtrk_dedx.clear();
1734 tmp_primtrk_hitx.clear();
1735 tmp_primtrk_hity.clear();
1736 tmp_primtrk_hitz.clear();
1737 tmp_primtrk_pitch.clear();
1738 tmp_primtrk_wid.clear();
1741 tmp_primtrk_hitx_true.clear();
1742 tmp_primtrk_hity_true.clear();
1743 tmp_primtrk_hitz_true.clear();
1744 tmp_primtrk_trkid_true.clear();
1745 tmp_primtrk_edept_true.clear();
1747 tmp_primtrk_true_x.clear();
1748 tmp_primtrk_true_y.clear();
1749 tmp_primtrk_true_z.clear();
1750 tmp_primtrk_true_trkid.clear();
1751 tmp_primtrk_true_edept.clear();
1752 tmp_primtrk_true_wid.clear();
1756 tmp_primtrj_true_x.clear();
1757 tmp_primtrj_true_y.clear();
1758 tmp_primtrj_true_z.clear();
1759 tmp_primtrj_true_edept.clear();
1761 for(
size_t k = 0;
k < calovector.size() &&
k<3;
k++){
1775 if(fmthm.isValid()){
1778 for (
size_t ii = 0; ii<vhit.size(); ++ii){
1779 bool fBadhit =
false;
1785 if (vmeta[ii]->
Index()>=tracklist[
fprimaryID]->NumberTrajectoryPoints()){
1786 throw cet::exception(
"Calorimetry_module.cc") <<
"Requested track trajectory index "<<vmeta[ii]->Index()<<
" exceeds the total number of trajectory points "<<tracklist[
fprimaryID]->NumberTrajectoryPoints()<<
" for track index "<<
fprimaryID<<
". Something is wrong with the track reconstruction. Please contact tjyang@fnal.gov!!";
1803 unsigned int tpc_no=1;
1804 if(xpos<=0 && zpos<232) tpc_no=1;
1805 if(xpos<=0 && zpos>232 && zpos<464) tpc_no=5;
1806 if(xpos<=0 && zpos>=464) tpc_no=9;
1807 if(xpos>0 && zpos<232) tpc_no=2;
1808 if(xpos>0 && zpos>232 && zpos<464) tpc_no=6;
1809 if(xpos>0 && zpos>=464) tpc_no=10;
1812 if (fBadhit)
continue;
1813 if (zpos<-100)
continue;
1814 planenum=vhit[ii]->WireID().Plane;
1827 x_c.push_back(xpos);
1828 y_c.push_back(ypos);
1829 z_c.push_back(zpos);
1834 tt_c.push_back(detProp.ConvertXToTicks(xpos, planenum, tpc_no, 0));
1837 ch_c.push_back(vhit[ii]->Channel());
1839 pt_c.push_back(vhit[ii]->PeakTime());
1840 q_c.push_back(vhit[ii]->Integral());
1841 a_c.push_back(vhit[ii]->PeakAmplitude());
1847 unsigned int wireno=vhit[ii]->WireID().Wire;
1849 wirez_c.push_back(xyzStart[2]);
1876 tt_v.push_back(detProp.ConvertXToTicks(xpos, planenum, tpc_no, 0));
1879 ch_v.push_back(vhit[ii]->Channel());
1881 pt_v.push_back(vhit[ii]->PeakTime());
1882 q_v.push_back(vhit[ii]->Integral());
1883 a_v.push_back(vhit[ii]->PeakAmplitude());
1889 unsigned int wireno=vhit[ii]->WireID().Wire;
1891 wirez_v.push_back(xyzStart[2]);
1899 tt_u.push_back(detProp.ConvertXToTicks(xpos, planenum, tpc_no, 0));
1902 ch_u.push_back(vhit[ii]->Channel());
1904 pt_u.push_back(vhit[ii]->PeakTime());
1905 q_u.push_back(vhit[ii]->Integral());
1906 a_u.push_back(vhit[ii]->PeakAmplitude());
1912 unsigned int wireno=vhit[ii]->WireID().Wire;
1914 wirez_u.push_back(xyzStart[2]);
1930 std::vector<float> Stw, Endw, Stt, Endt, Stwires, Endwires, Stticks, Endticks, TPCb, TPCcl;
1931 Stw.clear(); Endw.clear(); Stt.clear(); Endt.clear(); Stwires.clear(); Endwires.clear(); Stticks.clear(); Endticks.clear(); TPCb.clear(); TPCcl.clear();
1933 float numw, numt,wire_no,ticks_no;
1934 for(
size_t p1=0;p1<pfplist.size();p1++){
1935 std::vector<art::Ptr<recob::Track>>
trk=pftrack.at(p1);
1936 std::vector<art::Ptr<recob::Cluster>> allClusters=fmcp.at(p1);
1938 for(
size_t c1=0;
c1<allClusters.size();
c1++){
1941 Stw.push_back(allClusters[
c1]->StartWire());
1942 Endw.push_back(allClusters[
c1]->EndWire());
1943 Stt.push_back(allClusters[
c1]->StartTick());
1944 Endt.push_back(allClusters[
c1]->EndTick());
1945 TPCb.push_back(allClusters[
c1]->
Plane().
TPC);
1949 Stwires.push_back(allClusters[
c1]->StartWire());
1950 Endwires.push_back(allClusters[
c1]->EndWire());
1951 Stticks.push_back(allClusters[
c1]->StartTick());
1952 Endticks.push_back(allClusters[
c1]->EndTick());
1953 TPCcl.push_back(allClusters[
c1]->
Plane().
TPC);
1958 for(
size_t clt=0;clt<Stw.size();clt++){
1959 for(
size_t cl1=0;cl1<Stwires.size();cl1++){
1960 if(TPCcl[cl1]!=TPCb[clt])
continue;
1962 den=(Stw[clt]-Endw[clt])*(Stticks[cl1]-Endticks[cl1])-(Stt[clt]-Endt[clt])*(Stwires[cl1]-Endwires[cl1]);
1963 if(den==0)
continue;
1964 numw=(Stw[clt]*Endt[clt]-Stt[clt]*Endw[clt])*(Stwires[cl1]-Endwires[cl1])-(Stw[clt]-Endw[clt])*(Stwires[cl1]*Endticks[cl1]-Stticks[cl1]*Endwires[cl1]);
1965 numt=(Stw[clt]*Endt[clt]-Stt[clt]*Endw[clt])*(Stticks[cl1]-Endticks[cl1])-(Stt[clt]-Endt[clt])*(Stwires[cl1]*Endticks[cl1]-Stticks[cl1]*Endwires[cl1]);
1969 if(((Stw[clt]<wire_no && Endw[clt]>wire_no)||(Stw[clt]>wire_no && Endw[clt]<wire_no))&&((Stt[clt]<ticks_no && Endt[clt]>ticks_no)||(Stt[clt]>ticks_no && Endt[clt]<ticks_no)) && ((Stwires[cl1]<wire_no && Endwires[cl1]>wire_no)||(Stwires[cl1]>wire_no && Endwires[cl1]<wire_no)) && ((Stticks[cl1]<ticks_no && Endticks[cl1]>ticks_no)||(Stticks[cl1]>ticks_no && Endticks[cl1]<ticks_no))) {
1970 std::cout<<
"intersection wire and ticks are "<<std::round(wire_no)<<
" "<<ticks_no<<
" Stw Endw StT EndT "<<Stwires[cl1]<<
" "<<Endwires[cl1]<<
" "<<Stticks[cl1]<<
" "<<Endticks[cl1]<<
std::endl;
1973 unsigned int wireno=std::round(wire_no);
1976 std::cout<<
"Z position of intersection = "<<xyzStart[2]<<
" "<<xyzEnd[2]<<
" "<<wireno<<
std::endl;
1991 else if(thisShower != 0x0){
2009 if( (thisShower->
dEdx()).
size() > 0 )
2013 std::cout <<
"INFO::Primary pfParticle is not track or shower. Skip!" <<
std::endl;
2030 std::cout <<
"INFO::Number of daughters is " << particle->NumDaughters() <<
". Only the first NMAXDAUGTHERS are processed." <<
std::endl;
2033 for(
const int daughterID : particle->Daughters()){
2036 std::cout <<
"Daughter " << daughterID <<
" has " << daughterParticle->
NumDaughters() <<
" daughters" <<
std::endl;
2041 if(daughterTrack != 0x0){
2066 if(daughtercalovector.size() != 3)
2067 std::cerr <<
"WARNING::Calorimetry vector size for daughter is = " << daughtercalovector.size() <<
std::endl;
2069 for(
size_t k = 0;
k < daughtercalovector.size() &&
k<3;
k++){
2076 if(mcdaughterparticle != 0x0){
2101 std::cout <<
"Daughter Process = " << (mcdaughterparticle->
Process()).c_str()
2102 <<
" , mother = " << mcdaughterparticle->
Mother()
2106 else if(daughterShower != 0x0){
2122 if( (daughterShower->
dEdx()).
size() > 0 )
2126 std::cout <<
"INFO::Daughter pfParticle is not track or shower. Skip!" <<
std::endl;
2155 if(beamTriggerEvent)
double E(const int i=0) const
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
double fdaughterTheta[NMAXDAUGTHERS]
std::vector< double > q_c
std::vector< double > g4rw_p2
std::vector< double > beamtrk_Pz
std::vector< double > wid_c
double fprimaryMomentumByRangeMuon
int fdaughter_truth_Pdg[NMAXDAUGTHERS]
std::vector< double > Zintersection
std::vector< std::vector< double > > primtrk_hitz_true
std::vector< double > pt_v
unsigned int NumberTrajectoryPoints() const
std::vector< double > primtrk_range
const TVector3 & ShowerStart() const
double Z(const size_type i) const
double X(const size_type i) const
double EndMomentum() const
std::vector< double > wirez_u
std::vector< std::vector< int > > primtrk_wid
constexpr std::uint32_t timeLow() const
std::vector< double > beamPosz_spec
std::vector< std::string > interactionProcesslist
double Py(const int i=0) const
std::vector< double > beamDirz_spec
double g4rw_primary_singular_weight
std::vector< anab::T0 > GetPFParticleT0(const recob::PFParticle &particle, art::Event const &evt, std::string particleLabel) const
Get the T0(s) from a given PFParticle.
std::vector< double > pfphit_peaktime_v
int NumDaughters() const
Returns the number of daughter particles flowing from this one.
std::vector< double > wid_u
const simb::MCParticle * GetGeantGoodParticle(const simb::MCTruth &genTruth, const art::Event &evt) const
std::vector< int > wireno_u
int fprimaryShowerBestPlane
std::vector< std::vector< double > > primtrj_true_x
double E(const size_type i) const
int fdaughter_truth_Process[NMAXDAUGTHERS]
std::vector< std::vector< double > > primtrk_resrange
std::vector< double > g4rw_set_weights
double fprimaryShowerdEdx
double fprimary_truth_EndPosition[4]
double fdaughter_truth_EndMomentum[NMAXDAUGTHERS][4]
std::vector< double > a_c
double fprimaryShowerEnergy
const recob::Shower * GetPFParticleShower(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel, const std::string showerLabel) const
Get the shower associated to this particle. Returns a null pointer if not found.
std::vector< double > interaction_wid_c
std::vector< double > tt_c
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
const simb::MCTrajectory & Trajectory() const
std::string fGeneratorTag
std::string fPFParticleTag
std::vector< double > beamtrk_x
std::vector< std::vector< double > > primtrk_dedx
const recob::TrackTrajectory & Trajectory() const
Access to the stored recob::TrackTrajectory.
std::vector< double > timeintersection
std::vector< anab::Calorimetry > GetRecoTrackCalorimetry(const recob::Track &track, art::Event const &evt, const std::string trackModule, const std::string caloModule) const
Get the Calorimetry(s) from a given reco track.
const value_type & at(const size_type i) const
constexpr std::uint32_t timeHigh() const
std::vector< double > pfphit_wireid_c
double fdaughterEndMomentum[NMAXDAUGTHERS]
The data type to uniquely identify a Plane.
bool isValid
Whether this ID points to a valid element.
std::string KeyToProcess(unsigned char const &key) const
int fisdaughtershower[NMAXDAUGTHERS]
const std::vector< const recob::Hit * > GetPFParticleHits(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel) const
Get the hits associated to the PFParticle.
const std::vector< double > & Energy() const
std::vector< double > a_u
double Px(const int i=0) const
std::vector< double > pfphit_wireid_v
protoana::ProtoDUNEPFParticleUtils pfpUtil
double fdaughter_truth_P[NMAXDAUGTHERS]
const TVector3 GetPFParticleVertex(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel, const std::string trackLabel) const
Function to find the interaction vertex of a primary PFParticle.
std::vector< double > interactionAngles
double fprimary_truth_Mass
Plane(const Point_t &planePos, const Vector_t &planeDir)
Constructor from reference position on the plane and direction orthogonal to the plane.
std::vector< const sim::IDE * > TrackIdToSimIDEs_Ps(int const &id) const
std::vector< std::vector< double > > primtrk_hitz
std::vector< std::vector< double > > primtrk_edept_true
WireID_t Wire
Index of the wire within its plane.
std::vector< double > beamDirx_spec
std::vector< double > wirez_c
const art::InputTag fBeamModuleLabel
double fdaughter_truth_StartPosition[NMAXDAUGTHERS][4]
std::string Process() const
int fdaughterShowerBestPlane[NMAXDAUGTHERS]
unsigned char ProcessToKey(std::string const &process) const
std::vector< std::vector< double > > primtrk_hitx_true
std::vector< double > primtrk_ke_true
double fprimary_truth_EndMomentum[4]
double fprimaryStartDirection[3]
std::vector< double > a_v
int fprimary_truth_NDaughters
std::vector< std::vector< double > > primtrk_true_edept
std::vector< double > interaction_wid_v
std::vector< double > tt_u
double fprimary_truth_StartPosition[4]
int NumberDaughters() const
std::vector< double > q_u
std::vector< double > tt_v
std::vector< int > wireno_c
std::vector< double > beamtrk_z
double fdaughterPhi[NMAXDAUGTHERS]
std::vector< std::vector< double > > primtrk_true_x
std::vector< double > g4rw_p1
std::vector< double > beamtrk_Py
std::vector< double > interaction_tt_c
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
std::vector< double > beamtrk_y
protoana::ProtoDUNEDataUtils dataUtil
protoana::ProtoDUNETrackUtils trackUtil
std::vector< double > beamDiry_spec
std::vector< double > beamPosy_spec
double Pt(const int i=0) const
std::vector< double > pfphit_peaktime_c
int fdaughter_truth_TrackId[NMAXDAUGTHERS]
std::vector< std::vector< double > > primtrk_true_trkid
double fdaughterLength[NMAXDAUGTHERS]
double Length(size_t p=0) const
Access to various track properties.
int fprimary_truth_TrackId
double Y(const size_type i) const
const std::vector< double > & dEdx() const
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
const std::vector< double > & MIPEnergy() const
std::vector< fhicl::ParameterSet > ParSet
double fprimaryKineticEnergy[3]
std::vector< double > primtrk_ke_reco
std::vector< double > beamtrk_Px
double fdaughter_truth_EndPosition[NMAXDAUGTHERS][4]
const std::vector< double > & GetRecoBeamMomenta() const
std::vector< double > beamtrk_Eng
double Theta() const
Access to spherical or geographical angles at vertex or at any point.
double fdaughterMomentumByRangeProton[NMAXDAUGTHERS]
std::vector< double > Xintersection
geo::GeometryCore const * fGeometry
double P(const int i=0) const
const recob::Track * GetPFParticleTrack(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel, const std::string trackLabel) const
Get the track associated to this particle. Returns a null pointer if not found.
double fprimaryEndDirection[3]
double T(const int i=0) const
std::vector< double > interactionE
std::vector< double > q_v
const simb::MCParticle * GetMCParticleFromRecoTrack(detinfo::DetectorClocksData const &clockData, const recob::Track &track, art::Event const &evt, std::string trackModule) const
double fdaughterShowerMIPEnergy[NMAXDAUGTHERS]
const std::vector< const recob::PFParticle * > GetPFParticlesFromBeamSlice(art::Event const &evt, const std::string particleLabel) const
Return the pointers for the PFParticles in the beam slice. Returns an empty vector is no beam slice w...
const TVector3 & Direction() const
bool CreateRWTraj(const simb::MCParticle &part, const sim::ParticleList &plist, art::ServiceHandle< geo::Geometry > geo_serv, int event, G4ReweightTraj *theTraj)
double fdaughter_truth_Phi[NMAXDAUGTHERS]
double fdaughterOpeningAngle[NMAXDAUGTHERS]
std::vector< double > interactionZ
static int max(int a, int b)
The data type to uniquely identify a TPC.
std::vector< double > interaction_tt_v
std::vector< std::vector< double > > primtrk_true_y
double fdaughter_truth_Theta[NMAXDAUGTHERS]
double fdaughterStartDirection[NMAXDAUGTHERS][3]
int fprimary_truth_Isbeammatched
ProcessMap const & TrajectoryProcesses() const
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
std::string fprimary_truth_EndProcess
double fdaughter_truth_Momentum[NMAXDAUGTHERS][4]
std::vector< double > pfphit_peaktime_u
double fdaughterMomentumByRangeMuon[NMAXDAUGTHERS]
std::vector< std::string > g4rw_primary_var
std::vector< double > g4rw_primary_minus_sigma_weight
double fdaughterEndDirection[NMAXDAUGTHERS][3]
G4MultiReweighter MultiRW
const sim::ParticleList & ParticleList() const
Class storing the result of the Maximum Likelihood fit of Multiple Coulomb Scattering angles between ...
std::vector< std::vector< int > > primtrk_true_wid
std::vector< double > interactionX
Vector_t EndDirection() const
Returns the direction of the trajectory at the last point.
std::vector< int > wireno_v
double fdaughterRange[NMAXDAUGTHERS][3]
std::vector< std::vector< double > > primtrk_hitx
std::vector< std::vector< double > > primtrk_pitch
double fdaughterKineticEnergy[NMAXDAUGTHERS][3]
double Vx(const int i=0) const
double fdaughter_truth_Mass[NMAXDAUGTHERS]
std::vector< std::vector< double > > primtrk_trkid_true
std::vector< double > beamMomentum_spec
std::vector< double > primtrk_range_true
double StartMomentum() const
int fdaughterID[NMAXDAUGTHERS]
Hierarchical representation of particle flow.
std::vector< std::vector< double > > primtrj_true_edept
double fprimary_truth_StartPosition_MC[4]
std::vector< double > interaction_tt_u
std::vector< double > g4rw_primary_weights
int fdaughterNHits[NMAXDAUGTHERS]
std::vector< double > pfphit_wireid_u
double fprimaryEndPosition[3]
int GetNActiveFembsForAPA(art::Event const &evt, int apa) const
Get number of active fembs in an APA.
double fdaughterShowerEnergy[NMAXDAUGTHERS]
std::vector< double > Yintersection
static unsigned int reverse(QString &chars, unsigned char *level, unsigned int a, unsigned int b)
std::vector< std::vector< double > > primtrj_true_z
Point_t const & End() const
Returns the position of the last valid point of the trajectory [cm].
int fisdaughtertrack[NMAXDAUGTHERS]
const TLorentzVector & Momentum(const int i=0) const
std::vector< std::vector< double > > primtrk_hity
double Pz(const int i=0) const
art::ServiceHandle< geo::Geometry > fGeometryService_rw
double Vz(const int i=0) const
double fdaughterMomentum[NMAXDAUGTHERS]
Exception thrown on invalid wire number.
void WireEndPoints(geo::WireID const &wireid, double *xyzStart, double *xyzEnd) const
Fills two arrays with the coordinates of the wire end points.
const TVector3 GetPFParticleSecondaryVertex(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel, const std::string trackLabel) const
Function to find the secondary interaction vertex of a primary PFParticle.
geo::WireID NearestWireID(geo::Point_t const &point, geo::PlaneID const &planeid) const
Returns the ID of wire closest to position in the specified TPC.
const std::vector< recob::Track > & GetBeamTracks() const
std::vector< double > pos_ff
detail::Node< FrameID, bool > PlaneID
std::vector< double > interaction_wid_u
std::vector< double > wirez_v
std::vector< std::vector< double > > primtrk_hity_true
2D representation of charge deposited in the TDC/wire plane
std::vector< double > wid_v
std::vector< std::vector< double > > primtrk_true_z
std::vector< double > beamPosx_spec
std::string fCalorimetryTag
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
TPCID_t TPC
Index of the TPC within its cryostat.
double fprimary_truth_EndPosition_MC[4]
std::vector< std::vector< double > > primtrk_dqdx
unsigned int GetNumberPrimaryPFParticle(art::Event const &evt, const std::string particleLabel) const
Get the number of primary PFParticles.
std::vector< double > interactionY
double fprimary_truth_Momentum[4]
Vector_t StartDirection() const
Returns the direction of the trajectory at the first point.
second_as<> second
Type of time stored in seconds, in double precision.
geo::WireID suggestedWireID() const
Returns a better wire ID.
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
recob::tracking::Plane Plane
double fprimary_truth_Theta
Point_t const & Start() const
Returns the position of the first valid point of the trajectory [cm].
double fbeamtrackMomentum
double fprimaryStartPosition[3]
protoana::ProtoDUNETruthUtils truthUtil
double fprimaryEndMomentum
float GetBeamCosmicScore(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel) const
Access the BDT output used to decide if a slice is beam-like or cosmic-like.
std::vector< double > pt_u
double Vy(const int i=0) const
double fprimaryMomentumByRangeProton
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
double GetTrackMomentum(double trkrange, int pdg) const
double fprimary_truth_Phi
double fdaughterShowerdEdx[NMAXDAUGTHERS]
cet::coded_exception< error, detail::translate > exception
QTextStream & endl(QTextStream &s)
bool IsBeamTrigger(art::Event const &evt) const
double fdaughterEndPosition[NMAXDAUGTHERS][3]
std::vector< double > pt_c
double fprimaryOpeningAngle
double fdaughterStartPosition[NMAXDAUGTHERS][3]
std::string primtrk_end_g4process
std::vector< std::vector< double > > primtrj_true_y
std::vector< double > g4rw_primary_plus_sigma_weight
double fdaughter_truth_Pt[NMAXDAUGTHERS]
double fprimaryShowerMIPEnergy