diff --git a/.gitignore b/.gitignore index edeef8c..3f2c5b7 100644 --- a/.gitignore +++ b/.gitignore @@ -5,4 +5,4 @@ EventBuilder* *.root Mapper -anasenMS \ No newline at end of file +AnasenMS \ No newline at end of file diff --git a/Armory/ClassPW.h b/Armory/ClassPW.h index e430a7e..947f28b 100644 --- a/Armory/ClassPW.h +++ b/Armory/ClassPW.h @@ -188,26 +188,30 @@ inline void PW::FindWireID(TVector3 pos, TVector3 direction, bool verbose ){ //==== find the 2nd nearest wire short anode1 = hitInfo.nearestWire.first; + short aaa1 = anode1 - 1; if( aaa1 < 0 ) aaa1 += nWire; + short aaa2 = (anode1 + 1) % nWire; - double haha1 = Distance( pos, pos + direction, An[anode1-1].first, An[anode1-1].second); - double haha2 = Distance( pos, pos + direction, An[anode1+1].first, An[anode1+1].second); + double haha1 = Distance( pos, pos + direction, An[aaa1].first, An[aaa1].second); + double haha2 = Distance( pos, pos + direction, An[aaa2].first, An[aaa2].second); if( haha1 < haha2){ - hitInfo.nextNearestWire.first = anode1-1; + hitInfo.nextNearestWire.first = aaa1; hitInfo.nextNearestDist.first = haha1; }else{ - hitInfo.nextNearestWire.first = anode1+1; + hitInfo.nextNearestWire.first = aaa2; hitInfo.nextNearestDist.first = haha2; } short cathode1 = hitInfo.nearestWire.second; + short ccc1 = cathode1 - 1; if( ccc1 < 0 ) ccc1 += nWire; + short ccc2 = (cathode1 + 1) % nWire; - haha1 = Distance( pos, pos + direction, Ca[cathode1-1].first, Ca[cathode1-1].second); - haha2 = Distance( pos, pos + direction, Ca[cathode1+1].first, Ca[cathode1+1].second); + haha1 = Distance( pos, pos + direction, Ca[ccc1].first, Ca[ccc1].second); + haha2 = Distance( pos, pos + direction, Ca[ccc2].first, Ca[ccc2].second); if( haha1 < haha2){ - hitInfo.nextNearestWire.second = cathode1-1; + hitInfo.nextNearestWire.second = ccc1; hitInfo.nextNearestDist.second = haha1; }else{ - hitInfo.nextNearestWire.second = cathode1+1; + hitInfo.nextNearestWire.second = ccc2; hitInfo.nextNearestDist.second = haha2; } @@ -241,11 +245,11 @@ inline void PW::CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, bool verbose){ TVector3 shiftA2 = (An[anodeID2].second - An[anodeID1].second) * fracA; double totDistC = hitInfo.nearestDist.second + hitInfo.nextNearestDist.second; - double fracC = hitInfo.nearestDist.second / totDistA; + double fracC = hitInfo.nearestDist.second / totDistC; short cathodeID1 = hitInfo.nearestWire.second; short cathodeID2 = hitInfo.nextNearestWire.second; - TVector3 shiftC1 = (Ca[anodeID2].first - Ca[anodeID1].first) * fracC; - TVector3 shiftC2 = (Ca[anodeID2].second - Ca[anodeID1].second) * fracC; + TVector3 shiftC1 = (Ca[cathodeID2].first - Ca[cathodeID1].first) * fracC; + TVector3 shiftC2 = (Ca[cathodeID2].second - Ca[cathodeID1].second) * fracC; TVector3 a1 = An[anodeID1].first + shiftA1; TVector3 a2 = An[anodeID1].second + shiftA2; diff --git a/Armory/Makefile b/Armory/Makefile index 538e49e..ba76824 100644 --- a/Armory/Makefile +++ b/Armory/Makefile @@ -23,6 +23,6 @@ Mapper : Mapper.cpp ../mapping.h ClassDet.h @echo "--------- making Mapper" $(CC) $(COPTS) -o Mapper Mapper.cpp $(ROOTLIBS) -AnasenMS : ClassAnasen.h anasenMS.cpp Isotope.h constant.h +AnasenMS : constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassPW.h ClassAnasen.h anasenMS.cpp @echo "--------- making ANASEN Monte Carlo" - $(CC) $(COPTS) -o anasenMS anasenMS.cpp $(ROOTLIBS) \ No newline at end of file + $(CC) $(COPTS) -o AnasenMS anasenMS.cpp $(ROOTLIBS) diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index ecb3bce..238968d 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -90,9 +90,13 @@ int main(int argc, char **argv){ tree->Branch("sx3Y", &sx3Y, "sx3Y/D"); tree->Branch("sx3Z", &sx3Z, "sx3Z/D"); - int anodeID, cathodeID; - tree->Branch("aID", &anodeID, "anodeID/I"); - tree->Branch("cID", &cathodeID, "cathodeID/I"); + int anodeID[2], cathodeID[2]; + tree->Branch("aID", anodeID, "anodeID/I"); + tree->Branch("cID", cathodeID, "cathodeID/I"); + + double anodeDist[2], cathodeDist[2]; + tree->Branch("aDist", anodeDist, "anodeDist/D"); + tree->Branch("cDist", cathodeDist, "cathodeDist/D"); int sx3ID, sx3Up, sx3Dn, sx3Bk; double sx3ZFrac; @@ -106,6 +110,10 @@ int main(int argc, char **argv){ tree->Branch("reTheta", &reTheta, "reconstucted_theta/D"); tree->Branch("rePhi", &rePhi, "reconstucted_phi/D"); + double reTheta1, rePhi1; + tree->Branch("reTheta1", &reTheta1, "reconstucted_theta1/D"); + tree->Branch("rePhi1", &rePhi1, "reconstucted_phi1/D"); + //========timer TBenchmark clock; @@ -156,10 +164,17 @@ int main(int argc, char **argv){ pw->FindWireID(vertex, dir, false); sx3->FindSX3Pos(vertex, dir, false); - std::pair wireID = pw->GetNearestID(); + PWHitInfo hitInfo = pw->GetHitInfo(); - anodeID = wireID.first; - cathodeID = wireID.second; + anodeID[0] = hitInfo.nearestWire.first; + cathodeID[0] = hitInfo.nearestWire.second; + anodeID[1] = hitInfo.nextNearestWire.first; + cathodeID[1] = hitInfo.nextNearestWire.second; + + anodeDist[0] = hitInfo.nearestDist.first; + cathodeDist[0] = hitInfo.nearestDist.second; + anodeDist[1] = hitInfo.nextNearestDist.first; + cathodeDist[1] = hitInfo.nextNearestDist.second; sx3ID = sx3->GetID(); if( sx3ID >= 0 ){ @@ -172,11 +187,14 @@ int main(int argc, char **argv){ sx3Y = sx3->GetHitPos().Y(); sx3Z = sx3->GetHitPos().Z(); - pw->CalTrack(sx3->GetHitPos(), wireID.first, wireID.second, false); - + pw->CalTrack(sx3->GetHitPos(), anodeID[0], cathodeID[0], false); reTheta = pw->GetTrackTheta() * TMath::RadToDeg(); rePhi = pw->GetTrackPhi() * TMath::RadToDeg(); + pw->CalTrack2(sx3->GetHitPos(), hitInfo, false); + reTheta1 = pw->GetTrackTheta() * TMath::RadToDeg(); + rePhi1 = pw->GetTrackPhi() * TMath::RadToDeg(); + }else{ sx3Up = -1; sx3Dn = -1; @@ -192,7 +210,11 @@ int main(int argc, char **argv){ // } reTheta = TMath::QuietNaN(); - rePhi = TMath::QuietNaN(); + rePhi = TMath::QuietNaN(); + + reTheta1 = TMath::QuietNaN(); + rePhi1 = TMath::QuietNaN(); + } tree->Fill(); diff --git a/Armory/constant.h b/Armory/constant.h index 22cbd2f..db83ec5 100644 --- a/Armory/constant.h +++ b/Armory/constant.h @@ -7,7 +7,6 @@ * email: goluckyryan@gmail.com * ********************************************************************/ - #ifndef constant #define constant #include @@ -69,7 +68,7 @@ const double hbar = 197.326979; //====================================================================== -std::vector SplitStr(std::string tempLine, std::string splitter, int shift = 0){ +inline std::vector SplitStr(std::string tempLine, std::string splitter, int shift = 0){ std::vector output;