731 fTimeStamp = ts2.AsDouble();
734 bool beamTriggerEvent =
false;
736 auto mcTruths =
evt.getValidHandle<std::vector<simb::MCTruth>>(
fGeneratorTag);
737 if(!
evt.isRealData()){
740 auto beamHandle =
evt.getValidHandle<std::vector<beam::ProtoDUNEBeamEvent>>(
"generator");
741 std::vector<art::Ptr<beam::ProtoDUNEBeamEvent>> beamVec;
742 if( beamHandle.isValid()){
749 std::cout <<
"Number of reconstructed beam momenta from spec: " << momenta.size() <<
std::endl;
751 if( momenta.size() > 0 ) std::cout <<
"Measured Beam Momentum from spec: " << momenta.at(0) <<
std::endl;
754 for (
size_t i = 0; i<momenta.size(); ++i){
760 std::cout<<
"beam trk size:"<<btracks.size()<<
std::endl;
761 for (
size_t i = 0; i<btracks.size(); ++i){
762 std::cout<<
"beamPosx/beamPosy/beamPosz:"<<btracks[i].End().X()<<
"/"<<btracks[i].End().Y()<<
"/"<<btracks[i].End().Z()<<
std::endl;
763 std::cout<<
"beamDirx/beamDiry/beamDirz:"<<btracks[i].StartDirection().X()<<
"/"<<btracks[i].StartDirection().Y()<<
"/"<<btracks[i].StartDirection().Z()<<
std::endl;
777 if(geantGoodParticle != 0x0){
778 std::cout<<
"geant good particle loop "<<
std::endl;
779 std::cout <<
"Found GEANT particle corresponding to the good particle with pdg = " << geantGoodParticle->
PdgCode()
780 <<
" , track id = " << geantGoodParticle->
TrackId()
781 <<
" , Vx/Vy/Vz = " << geantGoodParticle->
Vx() <<
"/"<< geantGoodParticle->
Vy() <<
"/" << geantGoodParticle->
Vz()
785 std::vector<double> tmp_primtrk_truth_Z;
786 std::vector<double> tmp_primtrk_truth_Eng;
787 std::vector<double> tmp_primtrk_truth_trkide;
788 std::vector<int> tmp_wire;
789 std::vector<int> tmp_tpc;
791 beamTriggerEvent =
true;
803 beamid = geantGoodParticle->
TrackId();
811 std::cout <<
"Doing reweight" <<
std::endl;
823 std::cout<<
"G4reweighting loop step 1"<<
std::endl;
825 for (
int i = 0; i < n_rw_array; ++i) {
826 for (
int j = 0; j < n_rw_array; ++j) {
827 std::cout<<
"loop i, j "<<i<<
" "<<j<<
std::endl;
828 double tmp_p1=(react_st + i*d_react);
829 double tmp_p2=(elast_st + j*d_elast);
830 std::vector<double> input_values = {tmp_p1,tmp_p2};
832 std::cout<<
"loop temp1, 2 "<<tmp_p1<<
" "<<tmp_p2<<
std::endl;
833 bool set_values =
MultiRW.SetAllParameterValues(input_values);
834 std::cout<<
"setvalues "<<set_values<<
std::endl;
835 if (!set_values)
continue;
864 double pos_true[3]={geantGoodParticle->
Position(i_s).X(),geantGoodParticle->
Position(i_s).Y(),geantGoodParticle->
Position(i_s).Z()};
934 if (thisTrajectoryProcessMap1.size()){
935 for(
auto const& couple: thisTrajectoryProcessMap1){
949 std::cout<<
"interaction map size "<<thisTrajectoryProcessMap1.size()<<
std::endl;
950 if (thisTrajectoryProcessMap1.size()){
951 for(
auto const& couple1: thisTrajectoryProcessMap1){
953 if ((truetraj.
KeyToProcess(couple1.second)).find(
"CoulombScat")!= std::string::npos)
continue;
955 auto interactionPos4D = (truetraj.
at(couple1.first)).first ;
956 if (interactionPos4D.Z() <
minZ || interactionPos4D.Z() >
maxZ )
continue;
957 else if (interactionPos4D.X() <
minX || interactionPos4D.X() >
maxX )
continue;
958 else if (interactionPos4D.Y() <
minY || interactionPos4D.Y() >
maxY )
continue;
962 double xval=((truetraj.
at(couple1.first)).first).X();
963 double zval=((truetraj.
at(couple1.first)).first).Z();
964 unsigned int tpcno=1;
965 if(xval<=0 && zval<232) tpcno=1;
966 if(xval<=0 && zval>232 && zval<464) tpcno=5;
967 if(xval<=0 && zval>=464) tpcno=9;
968 if(xval>0 && zval<232) tpcno=2;
969 if(xval>0 && zval>232 && zval<464) tpcno=6;
970 if(xval>0 && zval>=464) tpcno=10;
972 interactionT.push_back(detProp.ConvertXToTicks(((truetraj.
at(couple1.first)).first).X(), 2, tpcno, 0));
977 std::cout<<
"number of interactions "<<thisTrajectoryProcessMap1.size()<<
std::endl;
978 std::cout<<
"int X, Y, Z and process "<<((truetraj.
at(couple1.first)).first).X()<<
" "<<((truetraj.
at(couple1.first)).first).Y()<<
" "<<((truetraj.
at(couple1.first)).first).Z()<<
" "<<truetraj.
KeyToProcess(couple1.second)<<
std::endl;
980 double interactionAngle = 999999.;
986 auto prevInteractionPos4D = (truetraj.
at(couple1.first-1)).first ;
987 auto prevInteractionPos3D = prevInteractionPos4D.Vect() ;
988 auto interactionPos3D = interactionPos4D.Vect() ;
989 auto distanceBtwPoint = interactionPos3D - prevInteractionPos3D;
991 if (truetraj.
size() > couple1.first + 1) {
993 auto nextInteractionPos4D = (truetraj.
at(couple1.first+1)).first ;
994 auto nextInteractionPos3D = nextInteractionPos4D.Vect() ;
995 auto distanceBtwPointNext = nextInteractionPos3D - interactionPos3D;
996 std::cout<<
"distance between points values "<<distanceBtwPoint.X()<<
" "<<distanceBtwPoint.Z()<<
" "<<distanceBtwPointNext.X()<<
" "<<distanceBtwPointNext.Z()<<
std::endl;
998 interactionAngles.push_back(TMath::ACos(distanceBtwPointNext.Dot(distanceBtwPoint)/(distanceBtwPointNext.Mag()*distanceBtwPoint.Mag() ) ));
1003 double z0=0.4792*prevInteractionPos3D.Z();
1004 double x0=prevInteractionPos3D.X();
1008 double z1=0.4792*interactionPos3D.Z();
1009 double x1=interactionPos3D.Z();
1013 double z2=0.4792*nextInteractionPos3D.Z();
1014 double x2=nextInteractionPos3D.Z();
1032 std::map<double, sim::IDE> orderedSimIDE;
1033 for (
auto& ide : simIDE_prim) orderedSimIDE[ide->z]= *ide;
1034 auto inTPCPoint = truetraj.
begin();
1035 auto Momentum0 = inTPCPoint->second;
1036 auto old_iter = orderedSimIDE.begin();
1038 double xi=0.0;
double yi=0.0;
double zi=0.0;
1042 for (
auto iter= orderedSimIDE.begin(); iter!= orderedSimIDE.end(); iter++,old_iter++){
1043 auto currentIde = iter->second;
1044 if(currentIde.z<
minZ)
continue;
1045 else if (currentIde.x <
minX || currentIde.x >
maxX )
continue;
1046 else if (currentIde.y <
minY || currentIde.y >
maxY )
continue;
1047 tmp_primtrk_truth_Z.push_back(currentIde.z);
1048 tmp_primtrk_truth_Eng.push_back(currentIde.energy);
1049 tmp_primtrk_truth_trkide.push_back(currentIde.trackID);
1050 double pos_true[3]={currentIde.x,currentIde.y,currentIde.z};
1063 int tpc_no=tpc1.
TPC;
1073 tmp_wire.push_back(wireID.
Wire);
1074 tmp_tpc.push_back(wireID.
TPC);
1077 tmp_wire.push_back(-9999);
1078 tmp_tpc.push_back(-9999);
1088 if(currentIde.trackID>=0){
1092 xi=currentIde.x;yi=currentIde.y;zi=currentIde.z;
1104 tmp_primtrk_truth_Z.clear();
1105 tmp_primtrk_truth_Eng.clear();
1106 tmp_primtrk_truth_trkide.clear();
1135 auto recoParticles =
evt.getValidHandle<std::vector<recob::PFParticle>>(
fPFParticleTag);
1149 std::vector<art::Ptr<recob::Hit> > hitlist;
1150 auto hitListHandle =
evt.getHandle< std::vector<recob::Hit> >(
"hitpdune");
1154 std::vector<art::Ptr<recob::Track> > tracklist;
1155 auto trackListHandle =
evt.getHandle< std::vector<recob::Track> >(
"pandoraTrack");
1158 art::FindManyP<recob::Track> thass(hitListHandle,
evt,
"pandoraTrack");
1160 std::vector<art::Ptr<recob::PFParticle> > pfplist;
1161 auto PFPListHandle =
evt.getHandle< std::vector<recob::PFParticle> >(
"pandora");
1164 std::vector<art::Ptr<recob::Cluster>> clusterlist;
1165 auto clusterListHandle =
evt.getHandle< std::vector<recob::Cluster> >(
"pandora");
1168 art::FindManyP<recob::Cluster> fmcp(PFPListHandle,
evt,
"pandora");
1169 art::FindManyP<recob::Track> pftrack(PFPListHandle,
evt,
"pandoraTrack");
1170 art::FindManyP<recob::Hit> clhit(clusterListHandle,
evt,
"pandora");
1172 std::cout<<
"number of pfp_particles "<<pfplist.size()<<
std::endl;
1173 std::cout<<
" size of pfParticles testing size "<<pfParticles.size()<<
std::endl;
1174 art::FindManyP<recob::Hit, recob::TrackHitMeta> fmthm(trackListHandle,
evt,
"pandoraTrack");
1177 std::cout<<
"outside pfparticle loop"<<
std::endl;
1194 double beamendx=-30;
1196 double beamendy=420;
1198 double beamendz=100;
1200 if(thisTrack != 0x0){
1202 beamstx=thisTrack->
Start().X();
1203 beamsty=thisTrack->
Start().Y();
1204 beamstz=thisTrack->
Start().Z();
1205 beamendx=thisTrack->
End().X();
1206 beamendy=thisTrack->
End().Y();
1207 beamendz=thisTrack->
End().Z();
1208 std::cout<<
"beamstx "<<beamstx<<
std::endl;
1211 std::vector<double> secondarystartx1;
1212 std::vector<double> secondarystarty1;
1213 std::vector<double> secondarystartz1;
1214 std::vector<double> secondaryendx1;
1215 std::vector<double> secondaryendy1;
1216 std::vector<double> secondaryendz1;
1217 std::vector<double> dQmichel1;
1218 std::vector<double> dQtrackbegin1;
1219 std::vector<double> dQtrackend1;
1220 std::vector<double> tracklengthsecondary1;
1221 std::vector<double> primsectheta1;
1222 std::vector<int> endhitssecondary1;
1223 std::vector<int> MtrackID1;
1224 std::vector<double> trks1;
1225 std::vector<double> ems1;
1226 std::vector<double> michels1;
1227 std::vector<double> nones1;
1228 size_t NTracks = tracklist.size();
1229 std::cout<<
"number of tracks "<<NTracks<<
std::endl;
1230 for(
size_t i=0;i<NTracks;i++){
1236 double startx=
pos.X();
1237 double starty=
pos.Y();
1238 double startz=
pos.Z();
1239 double endx=
end.X();
1240 double endy=
end.Y();
1241 double endz=
end.Z();
1242 if(track.
Length()<5)
continue;
1243 if(TMath::Max(endy,starty)>500 || TMath::Min(endy, starty)<200 || TMath::Max(startx, endx)>0||TMath::Min(startx,endx)<-200||TMath::Max(startz,endz)<30)
continue;
1246 std::vector<int> wirenos;
1247 std::vector<float> peakts,dqbuff1;
1248 std::vector<float> dQstart,dQend;
1249 std::vector<double> micheldq;
1250 wirenos.clear();peakts.clear();dqbuff1.clear();
1254 float zlast0=-99999;
1256 std::vector<std::tuple<double,double,double,double,int,double>> buff_ZYXTWQ;
1257 buff_ZYXTWQ.clear();
1258 double thetavalue=
theta12(beamstx,beamendx,beamsty,beamendy,beamstz,beamendz,startx,endx,starty,endy,startz,endz);
1259 if(fmthm.isValid()){
1260 auto vhit=fmthm.at(i);
1261 auto vmeta=fmthm.data(i);
1262 for (
size_t ii = 0; ii<vhit.size(); ++ii){
1263 bool fBadhit =
false;
1269 if (vmeta[ii]->
Index()>=tracklist[i]->NumberTrajectoryPoints()){
1270 throw cet::exception(
"Calorimetry_module.cc") <<
"Requested track trajectory index "<<vmeta[ii]->Index()<<
" exceeds the total number of trajectory points "<<tracklist[i]->NumberTrajectoryPoints()<<
" for track index "<<i<<
". Something is wrong with the track reconstruction. Please contact tjyang@fnal.gov!!";
1272 if (!tracklist[i]->HasValidPoint(vmeta[ii]->
Index())){
1278 auto loc = tracklist[i]->LocationAtPoint(vmeta[ii]->
Index());
1279 if (fBadhit)
continue;
1280 if (
loc.Z()<-100)
continue;
1281 if(vhit[ii]->
WireID().Plane==2){
1282 buff_ZYXTWQ.push_back(std::make_tuple(
loc.Z(),
loc.Y(),
loc.X(),vhit[ii]->PeakTime(),vhit[ii]->WireID().Wire,vhit[ii]->Integral()));
1284 peakts.push_back(vhit[ii]->PeakTime());
1288 wireno=vhit[ii]->WireID().Wire;
1289 peaktime=vhit[ii]->PeakTime();
1290 tpcno=vhit[ii]->WireID().TPC;
1305 double trk_score=0.0;
1307 double michel_score=0;
1308 double none_score=0;
1309 for(
size_t hitl=0;hitl<hitlist.size();hitl++){
1310 std::array<float,4> cnn_out=hitResults.getOutput(hitlist[hitl]);
1311 auto &
tracks = thass.at(hitlist[hitl].
key());
1312 if (!
tracks.empty() &&
tracks[0].key()!=ptrack.key() && tracklist[
tracks[0].key()]->Length()>25)
continue;
1315 float peakth1=hitlist[hitl]->PeakTime();
1316 int wireh1=hitlist[hitl]->WireID().Wire;
1317 for(
size_t m=0;
m<wirenos.size();
m++){
1318 if(wireh1==wirenos[
m] && peakth1==peakts[
m]){
1324 int planeid=hitlist[hitl]->WireID().Plane;
1325 int tpcid=hitlist[hitl]->WireID().TPC;
1326 if(
abs(wireh1-wireno)<15 &&
abs(peakth1-peaktime)<100 && planeid==2 && tpcid==tpcno){
1330 micheldq.push_back(hitlist[hitl]->Integral());
1331 trk_score+=cnn_out[hitResults.getIndex(
"track")];
1332 em_score+=cnn_out[hitResults.getIndex(
"em")];
1333 michel_score+=cnn_out[hitResults.getIndex(
"michel")];
1334 none_score+=cnn_out[hitResults.getIndex(
"none")];
1339 if(buff_ZYXTWQ.size()<10)
continue;
1340 sort(buff_ZYXTWQ.begin(),buff_ZYXTWQ.end());
1341 dQstart.clear(); dQend.clear();
1342 int qi11=buff_ZYXTWQ.size();
1343 for(
int qi=5;
qi<TMath::Min(15,qi11);
qi++){
1344 dQstart.push_back(std::get<5>(buff_ZYXTWQ[
qi]));
1347 for(
int qi=qi11-5;
qi<qi11;
qi++){
1348 dQend.push_back(std::get<5>(buff_ZYXTWQ[
qi]));
1350 secondarystartx1.push_back(startx);
1351 secondarystarty1.push_back(starty);
1352 secondarystartz1.push_back(startz);
1353 secondaryendx1.push_back(endx);
1354 secondaryendy1.push_back(endy);
1355 secondaryendz1.push_back(endz);
1356 endhitssecondary1.push_back(counter1);
1357 tracklengthsecondary1.push_back(track.
Length());
1358 dQmichel1.push_back(TMath::Median(micheldq.size(),&micheldq[0]));
1359 dQtrackbegin1.push_back(TMath::Median(dQstart.size(),&dQstart[0]));
1360 dQtrackend1.push_back(TMath::Median(dQend.size(),&dQend[0]));
1361 primsectheta1.push_back(thetavalue);
1362 MtrackID1.push_back(track.
ID());
1363 trks1.push_back(trk_score);
1364 ems1.push_back(em_score);
1365 michels1.push_back(michel_score);
1366 nones1.push_back(none_score);
1367 std::cout<<
"avg, trackscore, emscore, michelscore, nonescore "<<trk_score<<
" "<<em_score<<
" "<<michel_score<<
" "<<none_score<<
std::endl;
1388 endhitssecondary1.clear();
1389 secondarystartx1.clear();
1390 secondaryendx1.clear();
1391 secondarystarty1.clear();
1392 secondaryendy1.clear();
1393 secondarystartz1.clear();
1394 secondaryendz1.clear();
1396 primsectheta1.clear();
1397 dQtrackbegin1.clear();
1398 dQtrackend1.clear();
1399 tracklengthsecondary1.clear();
1409 if(thisTrack != 0x0){
1413 if(mcparticle!=0x0){
1416 truthid=mcparticle->
TrackId();
1462 if(calovector.size() != 3)
1463 std::cerr <<
"WARNING::Calorimetry vector size for primary is = " << calovector.size() <<
std::endl;
1464 std::vector<double> tmp_primtrk_dqdx;
1465 std::vector<double> tmp_primtrk_resrange;
1466 std::vector<double> tmp_primtrk_dedx;
1467 std::vector<double> tmp_primtrk_hitx;
1468 std::vector<double> tmp_primtrk_hity;
1469 std::vector<double> tmp_primtrk_hitz;
1470 std::vector<double> tmp_primtrk_pitch;
1471 std::vector<int> tmp_zwire;
1472 std::vector<int> tmp_ztpc;
1475 for (
auto &
calo : calovector) {
1476 if (
calo.PlaneID().Plane == 2){
1478 for (
size_t ihit = 0; ihit <
calo.dQdx().size(); ++ihit){
1479 tmp_primtrk_dqdx.push_back(
calo.dQdx()[ihit]);
1480 tmp_primtrk_resrange.push_back(
calo.ResidualRange()[ihit]);
1481 tmp_primtrk_dedx.push_back(
calo.dEdx()[ihit]);
1482 tmp_primtrk_pitch.push_back(
calo.TrkPitchVec()[ihit]);
1484 const auto &primtrk_pos=(
calo.XYZ())[ihit];
1485 tmp_primtrk_hitx.push_back(primtrk_pos.X());
1486 tmp_primtrk_hity.push_back(primtrk_pos.Y());
1487 tmp_primtrk_hitz.push_back(primtrk_pos.Z());
1489 double pos_true[3]={primtrk_pos.X(),primtrk_pos.Y(),primtrk_pos.Z()};
1503 int tpc_no=tpc2.
TPC;
1513 tmp_zwire.push_back(wireID.
Wire);
1514 tmp_ztpc.push_back(wireID.
TPC);
1517 tmp_zwire.push_back(-9999);
1518 tmp_ztpc.push_back(-9999);
1524 if(tmp_primtrk_dqdx.size()!=0){
1525 if(tmp_primtrk_hitz[0]>tmp_primtrk_hitz[tmp_primtrk_hitz.size()-1]){
1526 std::reverse(tmp_primtrk_hitz.begin(),tmp_primtrk_hitz.end());
1527 std::reverse(tmp_primtrk_hity.begin(),tmp_primtrk_hity.end());
1528 std::reverse(tmp_primtrk_hitx.begin(),tmp_primtrk_hitx.end());
1529 std::reverse(tmp_primtrk_pitch.begin(),tmp_primtrk_pitch.end());
1530 std::reverse(tmp_primtrk_dedx.begin(),tmp_primtrk_dedx.end());
1531 std::reverse(tmp_primtrk_dqdx.begin(),tmp_primtrk_dqdx.end());
1532 std::reverse(tmp_primtrk_resrange.begin(),tmp_primtrk_resrange.end());
1549 tmp_primtrk_dqdx.clear();
1550 tmp_primtrk_resrange.clear();
1551 tmp_primtrk_dedx.clear();
1552 tmp_primtrk_hitx.clear();
1553 tmp_primtrk_hity.clear();
1554 tmp_primtrk_hitz.clear();
1555 tmp_primtrk_pitch.clear();
1559 for(
size_t k = 0;
k < calovector.size() &&
k<3;
k++){
1572 if(fmthm.isValid()){
1575 for (
size_t ii = 0; ii<vhit.size(); ++ii){
1576 bool fBadhit =
false;
1582 if (vmeta[ii]->
Index()>=tracklist[
fprimaryID]->NumberTrajectoryPoints()){
1583 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!!";
1597 if (fBadhit)
continue;
1598 if (zpos<-100)
continue;
1599 planenum=vhit[ii]->WireID().Plane;
1603 peakT_2.push_back(vhit[ii]->PeakTime());
1604 int_2.push_back(vhit[ii]->Integral());
1614 int_1.push_back(vhit[ii]->Integral());
1618 int_0.push_back(vhit[ii]->Integral());
1628 std::vector<float> Stw, Endw, Stt, Endt, Stwires, Endwires, Stticks, Endticks, TPCb, TPCcl, primwireb, primtickb, wirebuf, tickbuf, primdqb;
1629 std::vector<std::vector<float> > primwire;
1630 std::vector<std::vector<float> > clwire;
1631 std::vector<std::vector<float> > primtick;
1632 std::vector<std::vector<float> > cltick;
1633 Stw.clear(); Endw.clear(); Stt.clear(); Endt.clear(); Stwires.clear(); Endwires.clear(); Stticks.clear(); Endticks.clear(); TPCb.clear(); TPCcl.clear();
1635 float numw, numt,wire_no,ticks_no;
1636 for(
size_t p1=0;p1<pfplist.size();p1++){
1637 std::vector<art::Ptr<recob::Track>>
trk=pftrack.at(p1);
1638 std::vector<art::Ptr<recob::Cluster>> allClusters=fmcp.at(p1);
1639 for(
size_t c1=0;
c1<allClusters.size();
c1++){
1640 if(trk.size() &&
int(trk[0].
key())==
fprimaryID) std::cout<<
"cluster tpc and plane "<<allClusters[
c1]->
Plane().TPC<<
" "<<allClusters[
c1]->Plane().Plane<<
" size "<<clhit.at(allClusters[
c1].
key()).size()<<
std::endl;
1641 if(allClusters[
c1]->
Plane().
TPC!=1)
continue;
1643 std::vector<art::Ptr<recob::Hit>> allHits=clhit.at(allClusters[
c1].
key());
1644 for(
size_t h1=0;
h1<allHits.size();
h1++){
1646 primtickb.push_back(allHits[
h1]->PeakTime());
1647 primdqb.push_back(allHits[
h1]->Integral());
1651 dq_0.push_back(primdqb);
1657 std::vector<art::Ptr<recob::Hit>> allHits=clhit.at(allClusters[
c1].
key());
1658 for(
size_t h1=0;
h1<allHits.size();
h1++){
1660 primtickb.push_back(allHits[
h1]->PeakTime());
1661 primdqb.push_back(allHits[
h1]->Integral());
1665 dq_1.push_back(primdqb);
1673 std::vector<art::Ptr<recob::Hit>> allHits=clhit.at(allClusters[
c1].
key());
1674 Stw.push_back(allClusters[
c1]->StartWire());
1675 Endw.push_back(allClusters[
c1]->EndWire());
1676 Stt.push_back(allClusters[
c1]->StartTick());
1677 Endt.push_back(allClusters[
c1]->EndTick());
1678 TPCb.push_back(allClusters[
c1]->
Plane().
TPC);
1679 for(
size_t h1=0;
h1<allHits.size();
h1++){
1681 primtickb.push_back(allHits[
h1]->PeakTime());
1682 primdqb.push_back(allHits[
h1]->Integral());
1686 dq_2.push_back(primdqb);
1687 primwire.push_back(primwireb);
1688 primtick.push_back(primtickb);
1693 else if(trk.size()){
1694 std::vector<art::Ptr<recob::Hit>> allHits=clhit.at(allClusters[
c1].
key());
1696 Stwires.push_back(allClusters[
c1]->StartWire());
1697 Endwires.push_back(allClusters[
c1]->EndWire());
1698 Stticks.push_back(allClusters[
c1]->StartTick());
1699 Endticks.push_back(allClusters[
c1]->EndTick());
1700 TPCcl.push_back(allClusters[
c1]->
Plane().
TPC);
1701 for(
size_t h2=0;
h2<allHits.size();
h2++){
1703 tickbuf.push_back(allHits[
h2]->PeakTime());
1705 clwire.push_back(wirebuf);
1706 cltick.push_back(tickbuf);
1713 for(
size_t clt=0;clt<Stw.size();clt++){
1714 for(
size_t cl1=0;cl1<Stwires.size();cl1++){
1715 if(TPCcl[cl1]!=TPCb[clt])
continue;
1716 den=(Stw[clt]-Endw[clt])*(Stticks[cl1]-Endticks[cl1])-(Stt[clt]-Endt[clt])*(Stwires[cl1]-Endwires[cl1]);
1718 if(Stwires[cl1]==Endwires[cl1]){
1719 wire_no=Stwires[cl1];
1720 ticks_no=Stticks[cl1];
1723 if(den==0 && !testv)
continue;
1725 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]);
1726 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]);
1730 if(((Stw[clt]-20<wire_no && Endw[clt]+20>wire_no)||(Stw[clt]+20>wire_no && Endw[clt]<wire_no-20))&&((Stt[clt]-100<ticks_no && Endt[clt]+100>ticks_no)||(Stt[clt]+100>ticks_no && Endt[clt]-100<ticks_no)) && ((Stwires[cl1]-20<wire_no && Endwires[cl1]+20>wire_no)||(Stwires[cl1]+20>wire_no && Endwires[cl1]-20<wire_no)) && ((Stticks[cl1]-100<ticks_no && Endticks[cl1]+100>ticks_no)||(Stticks[cl1]+100>ticks_no && Endticks[cl1]-100<ticks_no))){
1733 unsigned int wireno=std::round(wire_no);
1735 if(wireno>=0 && wireno<=479){
1750 float smallest=wireno;
1751 float biggest=wireno;
1754 for(
size_t hit1=0;hit1<primwire[clt].size();hit1++){
1755 if(primwire[clt][hit1]>=wireno-30 && primwire[clt][hit1]<=wireno+30){
1756 if(primwire[clt][hit1]>=wireno-30 && primwire[clt][hit1]<=wireno){
1757 small=primwire[clt][hit1];
1759 pws=primwire[clt][hit1];
1760 pts=primtick[clt][hit1];
1762 smallest=primwire[clt][hit1];
1764 if(primwire[clt][hit1]<=wireno+30 && primwire[clt][hit1]>=wireno){
1765 big=primwire[clt][hit1];
1767 pwb=primwire[clt][hit1];
1768 ptb=primtick[clt][hit1];
1770 biggest=primwire[clt][hit1];
1784 for(
size_t hit1=0;hit1<clwire[cl1].size();hit1++){
1785 if(clwire[cl1][hit1]>=wireno-30 && clwire[cl1][hit1]<=wireno+30){
1786 if(clwire[cl1][hit1]>=wireno-30 && clwire[cl1][hit1]<=wireno){
1787 csmall=clwire[cl1][hit1];
1788 if(csmall<smallest){
1789 cws=clwire[cl1][hit1];
1790 cts=cltick[cl1][hit1];
1792 smallest=clwire[cl1][hit1];
1794 if(clwire[cl1][hit1]<=wireno+30 && clwire[cl1][hit1]>=wireno){
1795 cbig=clwire[cl1][hit1];
1797 cwb=clwire[cl1][hit1];
1798 ctb=cltick[cl1][hit1];
1800 biggest=clwire[cl1][hit1];
1806 soln(pws, pwb, cws, cwb, pts, ptb, cts, ctb,ans);
1809 int wir=round(ans[0]);
1810 double xyzStart1[3];
1818 if(((Stw[clt]-5<ans[0] && Endw[clt]+5>ans[0])||(Stw[clt]+5>ans[0] && Endw[clt]<ans[0]-5))&&((Stt[clt]-50<ans[1] && Endt[clt]+50>ans[1])||(Stt[clt]+50>ans[1] && Endt[clt]-50<ans[1])) && ((Stwires[cl1]-5<ans[0] && Endwires[cl1]+5>ans[0])||(Stwires[cl1]+5>ans[0] && Endwires[cl1]-5<ans[0])) && ((Stticks[cl1]-50<ans[1] && Endticks[cl1]+50>ans[1])||(Stticks[cl1]+50>ans[1] && Endticks[cl1]-50<ans[1]))){
1819 if(wir>=0 && wir<=479){
1823 std::cout<<
"wire X position is here "<<xyzStart1[0]<<
std::endl;
1824 std::cout<<
"wire X position End is here "<<xyzEnd1[0]<<
std::endl;
1839 if(((Stw[clt]<=ans[0] && Endw[clt]>=ans[0])||(Stw[clt]>=ans[0] && Endw[clt]<=ans[0]))&&((Stt[clt]<=ans[1] && Endt[clt]>ans[1])||(Stt[clt]>ans[1] && Endt[clt]<ans[1])) && ((Stwires[cl1]<ans[0] && Endwires[cl1]>ans[0])||(Stwires[cl1]>ans[0] && Endwires[cl1]<ans[0])) && ((Stticks[cl1]<ans[1] && Endticks[cl1]>ans[1])||(Stticks[cl1]>ans[1] && Endticks[cl1]<ans[1]))){
1843 else if((Stw[clt]<=ans[0] && Endw[clt]>=ans[0])||(Stw[clt]>=ans[0] && Endw[clt]<=ans[0])){
1847 else if((Stwires[cl1]<=ans[0] && Endwires[cl1]>=ans[0])||(Stwires[cl1]>=ans[0] && Endwires[cl1]<=ans[0])){
1876 else if(thisShower != 0x0){
1894 if( (thisShower->
dEdx()).
size() > 0 )
1898 std::cout <<
"INFO::Primary pfParticle is not track or shower. Skip!" <<
std::endl;
1915 std::cout <<
"INFO::Number of daughters is " << particle->NumDaughters() <<
". Only the first NMAXDAUGTHERS are processed." <<
std::endl;
1918 for(
const int daughterID : particle->Daughters()){
1925 if(daughterTrack != 0x0){
1950 if(daughtercalovector.size() != 3)
1951 std::cerr <<
"WARNING::Calorimetry vector size for daughter is = " << daughtercalovector.size() <<
std::endl;
1953 for(
size_t k = 0;
k < daughtercalovector.size() &&
k<3;
k++){
1960 if(mcdaughterparticle != 0x0){
1985 std::cout <<
"Daughter Process = " << (mcdaughterparticle->
Process()).c_str()
1986 <<
" , mother = " << mcdaughterparticle->
Mother()
1990 else if(daughterShower != 0x0){
2006 if( (daughterShower->
dEdx()).
size() > 0 )
2010 std::cout <<
"INFO::Daughter pfParticle is not track or shower. Skip!" <<
std::endl;
2033 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 fprimaryStartPosition[3]
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
std::vector< std::vector< int > > MtrackID
std::vector< std::vector< float > > wireno_1
std::vector< std::vector< float > > wireno_0
std::vector< std::vector< double > > tracklengthsecondary
double fdaughterMomentum[NMAXDAUGTHERS]
std::vector< double > g4rw_react
unsigned int NumberTrajectoryPoints() const
std::vector< float > beamtrk_Pz
const TVector3 & ShowerStart() const
const TLorentzVector & Position(const int i=0) const
int fdaughter_truth_TrackId[NMAXDAUGTHERS]
double EndMomentum() const
constexpr std::uint32_t timeLow() const
int fdaughter_truth_Process[NMAXDAUGTHERS]
std::vector< std::vector< float > > peakTime_2
std::vector< double > interactionAngles
double fprimary_truth_Momentum[4]
double Py(const int i=0) const
std::vector< float > peakT_2
std::vector< std::vector< float > > peakTime_1
std::vector< std::vector< int > > endhitssecondary
std::vector< double > beamDirx_spec
std::vector< fhicl::ParameterSet > ParSet
const simb::MCParticle * GetGeantGoodParticle(const simb::MCTruth &genTruth, const art::Event &evt) const
std::vector< double > beamDirz_spec
int fdaughterID[NMAXDAUGTHERS]
std::string fCalorimetryTag
std::vector< double > beamMomentum_spec
std::vector< std::vector< double > > dQmichel
std::vector< std::vector< double > > trackscore
std::vector< std::vector< int > > primtrk_hitz_wire
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.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
std::vector< std::vector< double > > primtrk_dqdx
std::vector< double > beamPosx_spec
const simb::MCTrajectory & Trajectory() const
double fdaughterEndMomentum[NMAXDAUGTHERS]
const recob::TrackTrajectory & Trajectory() const
Access to the stored recob::TrackTrajectory.
double fprimary_truth_Theta
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.
double angle2d(double x0, double y0, double x1, double y1, double x2, double y2)
std::vector< std::vector< double > > primtrk_truth_Eng
std::vector< std::vector< float > > peakTime_0
const value_type & at(const size_type i) const
constexpr std::uint32_t timeHigh() const
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
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
double Px(const int i=0) const
std::vector< float > beamtrk_x
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 > Zintersection1
double fdaughterEndPosition[NMAXDAUGTHERS][3]
std::vector< std::vector< double > > primtrk_hitx
std::vector< std::vector< double > > emscore
int fprimary_truth_Isbeammatched
std::vector< const sim::IDE * > TrackIdToSimIDEs_Ps(int const &id) const
std::vector< std::vector< int > > primtrk_truth_Z_wire
double GetNTrajWeightFromSetPars(const std::vector< G4ReweightTraj * > &trajs, G4MultiReweighter &rw)
WireID_t Wire
Index of the wire within its plane.
std::vector< double > deltaZint
std::vector< float > beamtrk_Px
std::vector< std::vector< double > > dQtrackbegin
std::vector< std::vector< double > > primtrk_truth_trkide
std::string Process() const
unsigned char ProcessToKey(std::string const &process) const
double fdaughterMomentumByRangeProton[NMAXDAUGTHERS]
int fisdaughtershower[NMAXDAUGTHERS]
double fdaughterRange[NMAXDAUGTHERS][3]
int fisdaughtertrack[NMAXDAUGTHERS]
bool IsBeamlike(const recob::Track &, const art::Event &, std::string)
std::vector< std::vector< float > > wireno_2
double fprimary_truth_EndMomentum[4]
int NumberDaughters() const
double fprimaryOpeningAngle
std::vector< double > interactionV
double fdaughterOpeningAngle[NMAXDAUGTHERS]
double fdaughter_truth_Pt[NMAXDAUGTHERS]
std::vector< double > timeintersection
double fdaughter_truth_StartPosition[NMAXDAUGTHERS][4]
std::vector< std::vector< double > > secondaryendx
double fdaughterTheta[NMAXDAUGTHERS]
std::vector< double > primtrk_range
int fprimaryShowerBestPlane
std::vector< double > interactionY
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 > beamPosy_spec
std::vector< std::vector< double > > secondaryendz
std::vector< float > hitz_1
std::vector< int > beamtrk_z_tpc
std::vector< float > hitz_0
double Pt(const int i=0) const
double fprimaryKineticEnergy[3]
std::vector< double > interactionZ
double Length(size_t p=0) const
Access to various track properties.
std::vector< std::vector< double > > michelscore
const std::vector< double > & dEdx() const
double fdaughterStartPosition[NMAXDAUGTHERS][3]
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
const std::vector< double > & MIPEnergy() const
std::string EndProcess() const
double fdaughterLength[NMAXDAUGTHERS]
int fdaughterShowerBestPlane[NMAXDAUGTHERS]
double fprimaryShowerEnergy
Point_t const & Start() const
Access to track position at different points.
std::vector< std::vector< int > > primtrk_truth_Z_tpc
int fprimary_truth_TrackId
std::vector< std::vector< double > > primtrk_hity
const std::vector< double > & GetRecoBeamMomenta() const
std::vector< float > beamtrk_z
std::vector< std::vector< double > > secondarystarty
double Theta() const
Access to spherical or geographical angles at vertex or at any point.
double fdaughter_truth_Momentum[NMAXDAUGTHERS][4]
double fbeamtrackMomentum
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.
Point_t const & Vertex() const
std::vector< std::vector< double > > primtrk_dedx
double fprimary_truth_StartPosition[4]
double fdaughterShowerMIPEnergy[NMAXDAUGTHERS]
std::vector< std::vector< float > > dq_2
double T(const int i=0) const
double fprimary_truth_Phi
double fprimary_truth_EndPosition[4]
std::vector< std::vector< double > > primsectheta
std::vector< double > Zintersection
geo::GeometryCore const * fGeometry
std::vector< double > interactionAnglesVT
std::vector< float > int_2
const simb::MCParticle * GetMCParticleFromRecoTrack(detinfo::DetectorClocksData const &clockData, const recob::Track &track, art::Event const &evt, std::string trackModule) const
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
std::vector< std::vector< double > > primtrk_resrange
std::vector< int > beamtrk_z_wire
static int max(int a, int b)
std::vector< double > interactionU
The data type to uniquely identify a TPC.
std::vector< float > beamtrk_y
double fdaughter_truth_Mass[NMAXDAUGTHERS]
double fprimaryEndPosition[3]
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::vector< double > beamDiry_spec
std::vector< float > int_0
double fprimary_truth_tracklength
float soln(float x1, float x2, float x3, float x4, float y1, float y2, float y3, float y4, float result[2])
const sim::ParticleList & ParticleList() const
double fdaughter_truth_EndPosition[NMAXDAUGTHERS][4]
double fdaughterKineticEnergy[NMAXDAUGTHERS][3]
std::vector< double > g4rw_elast
double theta12(double x1, double x2, double y1, double y2, double z1, double z2, double x1p, double x2p, double y1p, double y2p, double z1p, double z2p)
std::vector< std::vector< double > > primtrk_truth_Z
Vector_t EndDirection() const
Returns the direction of the trajectory at the last point.
std::string fPFParticleTag
std::string fprimary_truth_EndProcess
double fprimaryMomentumByRangeMuon
std::vector< std::vector< int > > primtrk_hitz_tpc
double fprimaryMomentumByRangeProton
double Vx(const int i=0) const
double fdaughterPhi[NMAXDAUGTHERS]
std::vector< float > hitz_2
protoana::ProtoDUNEBeamCuts beam_cuts
double StartMomentum() const
Hierarchical representation of particle flow.
std::vector< double > deltatimeint
double fdaughter_truth_P[NMAXDAUGTHERS]
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
int fprimary_truth_Process
int fdaughterNHits[NMAXDAUGTHERS]
double fdaughterShowerdEdx[NMAXDAUGTHERS]
static unsigned int reverse(QString &chars, unsigned char *level, unsigned int a, unsigned int b)
std::vector< std::string > interactionProcesslist
std::vector< double > timeintersection1
int fdaughter_truth_Pdg[NMAXDAUGTHERS]
Point_t const & End() const
Returns the position of the last valid point of the trajectory [cm].
int fprimary_truth_NDaughters
std::vector< double > g4rw_set_weights
const TLorentzVector & Momentum(const int i=0) const
std::string fGeneratorTag
double Pz(const int i=0) const
std::vector< std::vector< double > > primtrk_hitz
double fprimaryShowerMIPEnergy
double Vz(const int i=0) const
std::vector< float > beamtrk_Eng
std::vector< float > beamtrk_Py
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.
double fdaughterMomentumByRangeMuon[NMAXDAUGTHERS]
double fprimaryEndMomentum
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.
Point_t const & End() const
std::vector< std::vector< float > > dq_1
std::string truth_last_process
std::vector< std::vector< double > > primtrk_pitch
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
detail::Node< FrameID, bool > PlaneID
double fprimaryStartDirection[3]
std::vector< double > interactionAnglesZT
std::vector< double > interactionT
std::vector< double > beamPosz_spec
std::vector< double > interactionW
G4MultiReweighter MultiRW
std::vector< std::vector< double > > dQtrackend
double fdaughterStartDirection[NMAXDAUGTHERS][3]
double fdaughter_truth_Theta[NMAXDAUGTHERS]
std::vector< double > interactionX
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
TPCID_t TPC
Index of the TPC within its cryostat.
unsigned int GetNumberPrimaryPFParticle(art::Event const &evt, const std::string particleLabel) const
Get the number of primary PFParticles.
protoana::ProtoDUNETruthUtils truthUtil
std::vector< G4ReweightTraj * > CreateNRWTrajs(const simb::MCParticle &part, const sim::ParticleList &plist, art::ServiceHandle< geo::Geometry > geo_serv, int event, std::string material_name="LAr", bool fVerbose=false)
Vector_t StartDirection() const
Returns the direction of the trajectory at the first point.
std::vector< std::vector< double > > secondarystartx
second_as<> second
Type of time stored in seconds, in double precision.
geo::WireID suggestedWireID() const
Returns a better wire ID.
recob::tracking::Plane Plane
double fdaughterEndDirection[NMAXDAUGTHERS][3]
Point_t const & Start() const
Returns the position of the first valid point of the trajectory [cm].
std::vector< std::vector< double > > secondaryendy
std::vector< float > int_1
double fprimary_truth_Mass
std::vector< std::vector< double > > nonescore
protoana::ProtoDUNEPFParticleUtils pfpUtil
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 > interactionAnglesUT
double Vy(const int i=0) const
double fprimaryEndDirection[3]
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
cet::coded_exception< error, detail::translate > exception
double fdaughterShowerEnergy[NMAXDAUGTHERS]
QTextStream & endl(QTextStream &s)
std::vector< std::vector< double > > secondarystartz
double fdaughter_truth_EndMomentum[NMAXDAUGTHERS][4]
double fdaughter_truth_Phi[NMAXDAUGTHERS]
double fprimaryShowerdEdx
std::vector< std::vector< float > > dq_0
protoana::ProtoDUNETrackUtils trackUtil