diff --git a/.vscode/settings.json b/.vscode/settings.json index 51aa676..daae08a 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -78,6 +78,14 @@ "variant": "cpp", "bitset": "cpp", "forward_list": "cpp", - "anasenMS.C": "cpp" + "anasenMS.C": "cpp", + "*.icc": "cpp", + "any": "cpp", + "cfenv": "cpp", + "csignal": "cpp", + "regex": "cpp", + "scoped_allocator": "cpp", + "shared_mutex": "cpp", + "valarray": "cpp" } } \ No newline at end of file diff --git a/Armory/ClassAnasen.h b/Armory/ClassAnasen.h index ee385ce..763001c 100644 --- a/Armory/ClassAnasen.h +++ b/Armory/ClassAnasen.h @@ -13,123 +13,43 @@ #include "TPolyLine3D.h" #include "TRandom.h" -struct SX3{ - short id = -1; // -1 when no hit - short chUp; - short chDown; - short chBack; - - double zFrac; // from +1 (downstream) to -1 (upstream) - - double eUp; - double eDown; - double eBack; - - TVector3 hitPos; - - void CalZFrac(){ - zFrac = (eUp - eDown)/(eUp + eDown); - } - - void Print(){ - if( id == -1 ){ - printf("Did not hit any SX3.\n"); - }else{ - printf("ID: %d, U,D,B: %d %d %d| zFrac : %.2f\n", id, chUp, chDown, chBack, zFrac); - printf("Hit Pos: %.2f, %.2f, %.2f\n", hitPos.X(), hitPos.Y(), hitPos.Z()); - } - } - -}; - -struct PW{ // proportional wire - // the nearest wire - short anode1 = -1; - short cathode1 = -1; - double anodeDis1 = 999999999; - double cathodeDis1 = 999999999; - - // the 2nd nearest wire - short anode2 = -1; - short cathode2 = -1; - double anodeDis2 = 999999999; - double cathodeDis2 = 999999999; - - void Print(){ - printf(" The nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", anode1, anodeDis1, cathode1, cathodeDis1); - printf(" The 2nd nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", anode2, anodeDis2, cathode2, cathodeDis2); - } - -}; +#include "ClassSX3.h" +#include "ClassPW.h" class ANASEN{ public: ANASEN(); ~ANASEN(); - void CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose = false); - TVector3 CalSX3Pos(unsigned short ID, unsigned short chUp, unsigned short chDown, unsigned short chBack, float eUp = 0, float eDown = 0 ); - - TVector3 GetTrackPos() const {return trackPos;} - TVector3 GetTrackVec() const {return trackVec;} - double GetTrackTheta() const {return trackVec.Theta();} - double GetTrackPhi() const {return trackVec.Phi();} - - void DrawAnasen(int anodeID1 = -1, int anodeID2 = -1, int cathodeID1 = -1, int cathodeID2 = -1, int sx3ID = -1, bool DrawQQQ = false ); - void DrawDeducedTrack(TVector3 sx3Pos, int anodeID, int cathodeID); - - //Simulation - SX3 FindSX3Pos(TVector3 pos, TVector3 direction, bool verbose = false); - PW FindWireID(TVector3 pos, TVector3 direction, bool verbose = false); void DrawTrack(TVector3 pos, TVector3 direction, bool drawEstimatedTrack = false); - - std::pair GetAnode(unsigned short id) const{return An[id];}; - std::pair GetCathode(unsigned short id) const{return Ca[id];}; + void DrawDeducedTrack(TVector3 sx3Pos, int anodeID, int cathodeID); + void DrawAnasen(int anodeID1 = -1, + int anodeID2 = -1, + int cathodeID1 = -1, + int cathodeID2 = -1, + int sx3ID = -1, + bool DrawQQQ = false ); private: - const int nWire = 24; - const int wireShift = 3; - const float zLen = 380; //mm - const float radiusA = 37; - const float radiusC = 43; - - std::vector> An; // the anode wire position vector in space - std::vector> Ca; // the cathode wire position vector in space - - std::vector> SDn; // coners of the SX3 0-11, z = mid point - std::vector> SUp; // coners of the SX3 12-23, z = mid point - std::vector SNorml; // normal of the SX3 (outward) - - void CalGeometry(); - - TVector3 trackPos; - TVector3 trackVec; - - double trackPosErrorZ; // mm - TVector3 tracePosErrorXY; // the mag is the size of the error - - TVector3 trackVecErrorA; // error vector prependicular to the Anode-Pos plan - TVector3 trackVecErrorC; // error vector prependicular to the Cathode-Pos plan - - const int nSX3 = 12; - const float sx3Radius = 88; - const float sx3Width = 40; - const float sx3Length = 75; - const float sx3Gap = 46; + PW pw; + SX3 sx3; const float qqqR1 = 50; const float qqqR2 = 100; - const float qqqZPos = sx3Gap/2 + sx3Length + 30; + const float qqqZPos = 23 + 75 + 30; + + void CalGeometry(); - // int geomID; TGeoManager *geom; TGeoVolume *worldBox; - void Construct3DModel(int anodeID1 = -1, int anodeID2 = -1, int cathodeID1 = -1, int cathodeID2 = -1, int sx3ID = -1, bool DrawQQQ = true); - - double Distance(TVector3 a1, TVector3 a2, TVector3 b1, TVector3 b2); - std::pair Intersect(TVector3 a1, TVector3 a2, TVector3 b1, TVector3 b2, bool verbose = false); + void Construct3DModel(int anodeID1 = -1, + int anodeID2 = -1, + int cathodeID1 = -1, + int cathodeID2 = -1, + int sx3ID = -1, + bool DrawQQQ = true); }; @@ -138,7 +58,6 @@ inline ANASEN::ANASEN(){ CalGeometry(); - // geomID = 0; geom = nullptr; worldBox = nullptr; @@ -152,61 +71,8 @@ inline ANASEN::~ANASEN(){ //!============================================== inline void ANASEN::CalGeometry(){ - std::pair p1; // anode - std::pair q1; // cathode - - //anode and cathode start at pos-Y axis and count in left-Hand - - //anode wire shift is right-hand. - //cathode wire shift is left-hand. - - for(int i = 0; i < nWire; i++ ){ - // Anode rotate right-hand - p1.first.SetXYZ( radiusA * TMath::Cos( TMath::TwoPi() / nWire * (-i) + TMath::PiOver2()), - radiusA * TMath::Sin( TMath::TwoPi() / nWire * (-i) + TMath::PiOver2()), - zLen/2); - p1.second.SetXYZ( radiusA * TMath::Cos( TMath::TwoPi() / nWire * (-i + wireShift) + TMath::PiOver2()), - radiusA * TMath::Sin( TMath::TwoPi() / nWire * (-i + wireShift) + TMath::PiOver2()), - -zLen/2); - An.push_back(p1); - - // Cathod rotate left-hand - q1.first.SetXYZ( radiusC * TMath::Cos( TMath::TwoPi() / nWire * (-i) + TMath::PiOver2()), - radiusC * TMath::Sin( TMath::TwoPi() / nWire * (-i) + TMath::PiOver2()), - zLen/2); - q1.second.SetXYZ( radiusC * TMath::Cos( TMath::TwoPi() / nWire * (- i - wireShift) + TMath::PiOver2()), - radiusC * TMath::Sin( TMath::TwoPi() / nWire * (- i - wireShift) + TMath::PiOver2()), - -zLen/2); - Ca.push_back(q1); - } - - // SX3 is couned in left-hand, started at neg-Y axis - - TVector3 sa, sb, sc, sn; - for(int i = 0; i < nSX3; i++){ - sa.SetXYZ( sx3Radius, -sx3Width/2, sx3Gap/2 + sx3Length/2 ); - sb.SetXYZ( sx3Radius, sx3Width/2, sx3Gap/2 + sx3Length/2 ); - - double rot = TMath::TwoPi() / nSX3 * (-i - 0.5) - TMath::PiOver2(); - - sa.RotateZ( rot ); - sb.RotateZ( rot ); - SDn.push_back(std::pair(sa,sb)); - - - sc.SetXYZ( sx3Radius, -sx3Width/2, sx3Gap/2 ); - sc.RotateZ( rot ); - - sn = ((sc-sa).Cross(sb-sa)).Unit(); - SNorml.push_back(sn); - - sa.SetXYZ( sx3Radius, -sx3Width/2, -sx3Gap/2 - sx3Length/2 ); - sb.SetXYZ( sx3Radius, sx3Width/2, -sx3Gap/2 - sx3Length/2 ); - - sa.RotateZ( rot ); - sb.RotateZ( rot ); - SUp.push_back(std::pair(sa,sb)); - } + sx3.ConstructGeo(); + pw.ConstructGeo(); } @@ -245,83 +111,73 @@ inline void ANASEN::Construct3DModel(int anodeID1, int anodeID2, int cathodeID1, worldBox->AddNode(axisZ, 1, new TGeoTranslation(0, 0, 5)); //.......... convert to wire center dimensions - double dAngle = wireShift * TMath::TwoPi() / nWire; - double wireALength = TMath::Sqrt( zLen*zLen + TMath::Power(2* radiusA * TMath::Sin(dAngle/2),2) ); - // double radiusAnew = radiusA * TMath::Cos( dAngle / 2.); - // double wireATheta = TMath::ATan2( 2* radiusA * TMath::Sin( dAngle / 2.), zLen); - - TGeoVolume *pcA = geom->MakeTube("tub1", Al, 0, 0.01, wireALength/2); + TGeoVolume *pcA = geom->MakeTube("tub1", Al, 0, 0.01, pw.GetAnodeLength()/2); pcA->SetLineColor(4); int startID = 0; - int endID = nWire - 1; + int endID = pw.GetNumWire() - 1; if( anodeID1 >= 0 && anodeID2 >= 0 ){ startID = anodeID1; endID = anodeID2; if( anodeID1 > anodeID2 ) { - endID = nWire + anodeID2; + endID = pw.GetNumWire() + anodeID2; } } for( int i = startID; i <= endID; i++){ - TVector3 a = (An[i].first + An[i].second) * 0.5; - double wireATheta = (An[i].first - An[i].second).Theta()* TMath::RadToDeg(); - double wireAPhi = (An[i].first - An[i].second).Phi() * TMath::RadToDeg() + 90; + TVector3 a = pw.GetAnodneMid(i); + double wireTheta = pw.GetAnodeTheta(i) * TMath::RadToDeg(); + double wirePhi = pw.GetAnodePhi(i) * TMath::RadToDeg() + 90; worldBox->AddNode(pcA, i+1, new TGeoCombiTrans( a.X(), a.Y(), a.Z(), - new TGeoRotation("rot1", wireAPhi , wireATheta, 0.))); + new TGeoRotation("rot1", wirePhi, wireTheta, 0.))); } - double wireCLength = TMath::Sqrt( zLen*zLen + TMath::Power(2* radiusC * TMath::Sin(dAngle/2),2) ); - // double radiusCnew = radiusC * TMath::Cos( dAngle / 2.); - // double wireCTheta = TMath::ATan2( 2* radiusC * TMath::Sin( dAngle / 2.), zLen); - - TGeoVolume *pcC = geom->MakeTube("tub2", Al, 0, 0.01, wireCLength/2); + TGeoVolume *pcC = geom->MakeTube("tub2", Al, 0, 0.01, pw.GetCathodeLength()/2); pcC->SetLineColor(6); startID = 0; - endID = nWire - 1; + endID = pw.GetNumWire() - 1; if( cathodeID1 >= 0 && cathodeID2 >= 0 ){ startID = cathodeID1; endID = cathodeID2; if( cathodeID1 > cathodeID2 ) { - endID = nWire + cathodeID2; + endID = pw.GetNumWire() + cathodeID2; } } for( int i = startID; i <= endID; i++){ - TVector3 a = (Ca[i].first + Ca[i].second) * 0.5; - double wireATheta = (Ca[i].first - Ca[i].second).Theta()* TMath::RadToDeg(); - double wireAPhi = (Ca[i].first - Ca[i].second).Phi() * TMath::RadToDeg() + 90; + TVector3 a = pw.GetCathodneMid(i); + double wireTheta = pw.GetCathodeTheta(i) * TMath::RadToDeg(); + double wirePhi = pw.GetCathodePhi(i) * TMath::RadToDeg() + 90; worldBox->AddNode(pcC, i+1, new TGeoCombiTrans( a.X(), a.Y(), a.Z(), - new TGeoRotation("rot1", wireAPhi , wireATheta, 0.))); + new TGeoRotation("rot1", wirePhi , wireTheta, 0.))); } - TGeoVolume * sx3 = geom->MakeBox("box", Al, 0.1, sx3Width/2, sx3Length/2); - sx3->SetLineColor(kGreen+3); + TGeoVolume * sx3Det = geom->MakeBox("box", Al, 0.1, sx3.GetWidth()/2, sx3.GetLength()/2); + sx3Det->SetLineColor(kGreen+3); - for( int i = 0; i < nSX3; i++){ + for( int i = 0; i < sx3.GetNumDet(); i++){ if( sx3ID != -1 && i != sx3ID ) continue; - TVector3 aUp = (SUp[i].first + SUp[i].second)*0.5; // center of the SX3 upstream - TVector3 aDn = (SDn[i].first + SDn[i].second)*0.5; // center of the SX3 Downstream - double phi = (SUp[i].second - SUp[i].first).Phi() * TMath::RadToDeg() + 90; - - worldBox->AddNode(sx3, 2*i+1., new TGeoCombiTrans( aUp.X(), - aUp.Y(), - aUp.Z(), - new TGeoRotation("rot1", phi, 0., 0.))); - worldBox->AddNode(sx3, 2*i+1., new TGeoCombiTrans( aDn.X(), - aDn.Y(), - aDn.Z(), - new TGeoRotation("rot1", phi, 0., 0.))); + TVector3 aUp = sx3.GetUpMid(i); // center of the SX3 upstream + TVector3 aDn = sx3.GetDnMid(i); // center of the SX3 Downstream + double phi = sx3.GetDetPhi(i) * TMath::RadToDeg() + 90; + worldBox->AddNode(sx3Det, 2*i+1., new TGeoCombiTrans( aUp.X(), + aUp.Y(), + aUp.Z(), + new TGeoRotation("rot1", phi, 0., 0.))); + worldBox->AddNode(sx3Det, 2*i+1., new TGeoCombiTrans( aDn.X(), + aDn.Y(), + aDn.Z(), + new TGeoRotation("rot1", phi, 0., 0.))); } if( DrawQQQ ){ @@ -337,167 +193,6 @@ inline void ANASEN::Construct3DModel(int anodeID1, int anodeID2, int cathodeID1, } -//!============================================== Aux Functions -inline double ANASEN::Distance(TVector3 a1, TVector3 a2, TVector3 b1, TVector3 b2){ - TVector3 na = a1 - a2; - TVector3 nb = b1 - b2; - TVector3 nd = (na.Cross(nb)).Unit(); - return TMath::Abs(nd.Dot(a1-b2)); -} - - -inline std::pair ANASEN::Intersect(TVector3 An, TVector3 p2, TVector3 q1, TVector3 q2, bool verbose){ - - //see https://nukephysik101.wordpress.com/2023/12/30/intersect-between-2-line-segments/ - //zero all z-component - TVector3 a0 = An; a0.SetZ(0); - TVector3 a1 = p2; a1.SetZ(0); - - TVector3 b0 = q1; b0.SetZ(0); - TVector3 b1 = q2; b1.SetZ(0); - - double A = ((b0-b1).Cross(a0-a1)).Mag(); - - double h = ((b0-a0).Cross(b1-a0)).Z()/ A; - double k = ((a1-b0).Cross(a0-b0)).Z()/ A; - - if( verbose ) printf(" ----h, k : %f, %f\n", h, k); - - return std::pair(h,k); -} - -//!============================================== Given a position and a direction, find wireID and SX3 position -inline PW ANASEN::FindWireID(TVector3 pos, TVector3 direction, bool verbose ){ - - PW pw; - - double phi = direction.Phi(); - - for( int i = 0; i < nWire; i++){ - - double disA = 99999999; - double disC = 99999999; - - double phiS = An[i].first.Phi() - TMath::PiOver4(); - double phiL = An[i].second.Phi() + TMath::PiOver4(); - - // printf("A%2d: %f %f | %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg(), phi * TMath::RadToDeg()); - - if( phi > 0 && phiS > phiL ) { - phiL = phiL + TMath::TwoPi(); - // printf("------ %f %f\n", phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg()); - } - if( phi < 0 && phiS > phiL ) { - phiS = phiS - TMath::TwoPi(); - // printf("------ %f %f\n", phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg()); - } - - if( phiS < phi && phi < phiL) { - disA = Distance( pos, pos + direction, An[i].first, An[i].second); - if( disA < pw.anodeDis1 ){ - pw.anodeDis1 = disA; - pw.anode1 = i; - } - } - - phiS = Ca[i].second.Phi()- TMath::PiOver4(); - phiL = Ca[i].first.Phi() + TMath::PiOver4(); - // printf("C%2d: %f %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg()); - if( phi > 0 && phiS > phiL ) { - phiL = phiL + TMath::TwoPi(); - // printf("------ %f %f\n", phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg()); - } - if( phi < 0 && phiS > phiL ) { - phiS = phiS - TMath::TwoPi(); - // printf("------ %f %f\n", phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg()); - } - - if(phiS < phi && phi < phiL) { - disC = Distance( pos, pos + direction, Ca[i].first, Ca[i].second); - - if( disC < pw.cathodeDis1 ){ - pw.cathodeDis1 = disC; - pw.cathode1 = i; - } - } - - if(verbose) printf(" %2d | %8.2f, %8.2f\n", i, disA, disC); - } - - //==== find the 2nd nearest wire - double haha1 = Distance( pos, pos + direction, An[pw.anode1-1].first, An[pw.anode1-1].second); - double haha2 = Distance( pos, pos + direction, An[pw.anode1+1].first, An[pw.anode1+1].second); - if( haha1 < haha2){ - pw.anode2 = pw.anode1-1; - pw.anodeDis2 = haha1; - }else{ - pw.anode2 = pw.anode1+1; - pw.anodeDis2 = haha2; - } - - haha1 = Distance( pos, pos + direction, Ca[pw.cathode1-1].first, Ca[pw.cathode1-1].second); - haha2 = Distance( pos, pos + direction, Ca[pw.cathode1+1].first, Ca[pw.cathode1+1].second); - if( haha1 < haha2){ - pw.cathode2 = pw.cathode1-1; - pw.cathodeDis2 = haha1; - }else{ - pw.cathode2 = pw.cathode1+1; - pw.cathodeDis2 = haha2; - } - - - if( verbose ) pw.Print(); - return pw; -} - -inline SX3 ANASEN::FindSX3Pos(TVector3 pos, TVector3 direction, bool verbose){ - - SX3 haha; - - haha.id = -1; - for( int i = 0 ; i < nSX3; i++){ - - if(verbose) printf(" %d ", i); - std::pair frac = Intersect( pos, pos + direction, SDn[i].first, SDn[i].second, verbose); - - - if( frac.second < 0 || frac.second > 1 ) continue; - haha.hitPos = pos + frac.first * direction; - - double dis = haha.hitPos.Dot(SNorml[i]); - - if(verbose) { - printf("reduced distance : %f\n", dis); - printf(" %d*", (i+1)%nSX3); - Intersect( pos, pos + direction, SDn[(i+1)%nSX3].first, SDn[(i+1)%nSX3].second, verbose); - } - - if( TMath::Abs(dis - sx3Radius) > 0.1 ) continue; - - haha.chDown = 2 * TMath::Floor(frac.second * 4); - haha.chUp = haha.chDown + 1; - - double zPos = haha.hitPos.Z(); - if( (sx3Gap/2 < zPos && zPos < sx3Gap/2 + sx3Length ) || (-sx3Gap/2 - sx3Length < zPos && zPos < -sx3Gap/2 ) ){ - - haha.id = zPos > 0 ? i : i + 12; - - haha.zFrac = zPos > 0 ? (zPos - sx3Gap/2. - sx3Length/2.)/sx3Length : (zPos - ( - sx3Gap/2. - sx3Length/2.) )/sx3Length ; - - haha.chBack = TMath::Floor( (haha.zFrac + 0.5) * 4 ) + 8; - - if( verbose) haha.Print(); - - return haha; - }else{ - if( verbose ) printf(" zPos out of sensitive region\n"); - } - } - - if( verbose) haha.Print(); - return haha; -} - //!============================================== Drawing functions inline void ANASEN::DrawAnasen(int anodeID1, int anodeID2, int cathodeID1, int cathodeID2, int sx3ID, bool DrawQQQ ){ @@ -511,21 +206,17 @@ inline void ANASEN::DrawAnasen(int anodeID1, int anodeID2, int cathodeID1, int c inline void ANASEN::DrawTrack(TVector3 pos, TVector3 direction, bool drawEstimatedTrack){ - PW pw = FindWireID(pos, direction); + pw.FindWireID(pos, direction); + sx3.FindSX3Pos(pos, direction); - SX3 sx3 = FindSX3Pos(pos, direction); - - int a1 = pw.anode1 - 1; if( a1 < 0 ) a1 += nWire; - int b1 = pw.cathode1 - 1; if( b1 < 0 ) b1 += nWire; - - //Construct3DModel(a1, id.first+1, b1, id.second+1, false); - Construct3DModel(pw.anode1, pw.anode1, pw.cathode1, pw.cathode1, -1, false); + std::pair wireID = pw.GetNearestID(); + + Construct3DModel(wireID.first, wireID.first, wireID.second, wireID.second, -1, false); double theta = direction.Theta() * TMath::RadToDeg(); double phi = direction.Phi() * TMath::RadToDeg(); // printf("Theta, Phi = %.2f %.2f \n", theta, phi); // pos.Print(); - TGeoVolume * Track = geom->MakeTube("track", 0, 0, 0.1, 150.); Track->SetLineColor(kRed); worldBox->AddNode(Track, 1, new TGeoCombiTrans( pos.X(), pos.Y(), pos.Z(), new TGeoRotation("rotA", phi + 90, theta, 0.))); @@ -534,20 +225,23 @@ inline void ANASEN::DrawTrack(TVector3 pos, TVector3 direction, bool drawEstima startPos->SetLineColor(kBlack); worldBox->AddNode(startPos, 3, new TGeoCombiTrans( pos.X(), pos.Y(), pos.Z(), new TGeoRotation("rotA", 0, 0, 0.))); - if( sx3.id >= 0 ){ + if( sx3.GetID() >= 0 ){ TGeoVolume * hit = geom->MakeSphere("hitpos", 0, 0, 3); hit->SetLineColor(kRed); - worldBox->AddNode(hit, 2, new TGeoCombiTrans( sx3.hitPos.X(), sx3.hitPos.Y(), sx3.hitPos.Z(), new TGeoRotation("rotA", 0, 0, 0.))); + + TVector3 hitPos = sx3.GetHitPos(); + + worldBox->AddNode(hit, 2, new TGeoCombiTrans( hitPos.X(), hitPos.Y(), hitPos.Z(), new TGeoRotation("rotA", 0, 0, 0.))); if( drawEstimatedTrack ){ - CalTrack(sx3.hitPos, pw.anode1, pw.cathode1, true); + pw.CalTrack(hitPos, wireID.first, wireID.second, true); - double thetaDeduce = trackVec.Theta() * TMath::RadToDeg(); - double phiDeduce = trackVec.Phi() * TMath::RadToDeg(); + double thetaDeduce = pw.GetTrackTheta() * TMath::RadToDeg(); + double phiDeduce = pw.GetTrackPhi() * TMath::RadToDeg(); TGeoVolume * trackDeduce = geom->MakeTube("trackDeduce", 0, 0, 0.1, 100.); trackDeduce->SetLineColor(kOrange); - worldBox->AddNode(trackDeduce, 1, new TGeoCombiTrans( sx3.hitPos.X(), sx3.hitPos.Y(), sx3.hitPos.Z(), new TGeoRotation("rotA", phiDeduce + 90, thetaDeduce, 0.))); + worldBox->AddNode(trackDeduce, 1, new TGeoCombiTrans( hitPos.X(), hitPos.Y(), hitPos.Z(), new TGeoRotation("rotA", phiDeduce + 90, thetaDeduce, 0.))); } @@ -561,12 +255,12 @@ inline void ANASEN::DrawTrack(TVector3 pos, TVector3 direction, bool drawEstima inline void ANASEN::DrawDeducedTrack(TVector3 sx3Pos, int anodeID, int cathodeID){ - CalTrack(sx3Pos, anodeID, cathodeID); + pw.CalTrack(sx3Pos, anodeID, cathodeID); Construct3DModel(anodeID, anodeID, cathodeID, cathodeID, -1, false); - double theta = trackVec.Theta() * TMath::RadToDeg(); - double phi = trackVec.Phi() * TMath::RadToDeg(); + double theta = pw.GetTrackTheta() * TMath::RadToDeg(); + double phi = pw.GetTrackPhi() * TMath::RadToDeg(); TGeoVolume * Track = geom->MakeTube("axisX", 0, 0, 0.1, 100.); Track->SetLineColor(kRed); @@ -576,63 +270,11 @@ inline void ANASEN::DrawDeducedTrack(TVector3 sx3Pos, int anodeID, int cathodeI hit->SetLineColor(kRed); worldBox->AddNode(hit, 2, new TGeoCombiTrans( sx3Pos.X(), sx3Pos.Y(), sx3Pos.Z(), new TGeoRotation("rotA", 0, 0, 0.))); - geom->CloseGeometry(); geom->SetVisLevel(4); worldBox->Draw("ogle"); } -//!============================================== Duduce trace from experiment -inline void ANASEN::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose){ - - trackPos = sx3Pos; - - TVector3 n1 = (An[anodeID].first - An[anodeID].second).Cross((sx3Pos - An[anodeID].second)).Unit(); - TVector3 n2 = (Ca[cathodeID].first - Ca[cathodeID].second).Cross((sx3Pos - Ca[cathodeID].second)).Unit(); - - // if the handiness of anode and cathode revered, it should be n2 cross n1 - trackVec = (n2.Cross(n1)).Unit(); - - if( verbose ) printf("Theta, Phi = %f, %f \n", trackVec.Theta() *TMath::RadToDeg(), trackVec.Phi()*TMath::RadToDeg()); - -} - -inline TVector3 ANASEN::CalSX3Pos(unsigned short ID, unsigned short chUp, unsigned short chDown, unsigned short chBack, float eUp, float eDown){ - - TVector3 haha; - - if( (chUp - chDown) != 1 || (chDown % 2) != 0) return haha; - - int reducedID = ID % nSX3; - - TVector3 sa, sb; - - if( ID < nSX3 ){ //down - - sa = SDn[reducedID].second; - sb = SDn[reducedID].first; - - }else{ - - sa = SUp[reducedID].second; - sb = SUp[reducedID].first; - - } - - haha.SetX( (sb.X() - sa.X()) * chUp/8 + sa.X()); - haha.SetY( (sb.Y() - sa.Y()) * chUp/8 + sa.Y()); - - if( eUp == 0 || eDown == 0 ){ - haha.SetZ( sa.Z() + (2*(chBack - 7)-1) * sx3Length / 8 ); - }else{ - double frac = (eUp - eDown)/(eUp + eDown); // from +1 (downstream) to -1 (upstream) - double zPos = sa.Z() + sx3Length * frac/2; - haha.SetZ( zPos ); - } - - return haha; - -} #endif \ No newline at end of file diff --git a/Armory/ClassPW.h b/Armory/ClassPW.h new file mode 100644 index 0000000..d0273ad --- /dev/null +++ b/Armory/ClassPW.h @@ -0,0 +1,214 @@ +#ifndef ClassPW_h +#define ClassPW_h + +#include +#include +#include + +class PW{ // proportional wire + +public: + PW(){ Clear(); }; + ~PW(){}; + + std::pair GetNearestID() const {return std::pair(anode1, cathode1);} + std::pair GetNearestDistance() const {return std::pair(anodeDis1,cathodeDis1);} + std::pair Get2ndNearestID() const {return std::pair(anode2,cathode2);} + std::pair Get2ndNearestDistance() const {return std::pair(anodeDis2,cathodeDis2);} + + TVector3 GetTrackPos() const {return trackPos;} + TVector3 GetTrackVec() const {return trackVec;} + double GetTrackTheta() const {return trackVec.Theta();} + double GetTrackPhi() const {return trackVec.Phi();} + + int GetNumWire() const {return nWire;} + double GetDeltaAngle() const {return dAngle;} + double GetAnodeLength() const {return anodeLength;} + double GetCathodeLength() const {return cathodeLength;} + TVector3 GetAnodeDn(short id) const {return An[id].first;} + TVector3 GetAnodeUp(short id) const {return An[id].second;} + TVector3 GetCathodeDn(short id) const {return Ca[id].first;} + TVector3 GetCathodeUp(short id) const {return Ca[id].second;} + + TVector3 GetAnodneMid(short id) const {return (An[id].first + An[id].second) * 0.5; } + double GetAnodeTheta(short id) const {return (An[id].first - An[id].second).Theta();} + double GetAnodePhi(short id) const {return (An[id].first - An[id].second).Phi();} + + TVector3 GetCathodneMid(short id) const {return (Ca[id].first + Ca[id].second) * 0.5; } + double GetCathodeTheta(short id) const {return (Ca[id].first - Ca[id].second).Theta();} + double GetCathodePhi(short id) const {return (Ca[id].first - Ca[id].second).Phi();} + + void Clear(); + void ConstructGeo(); + void FindWireID(TVector3 pos, TVector3 direction, bool verbose = false); + void CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose = false); + + void Print(){ + printf(" The nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", anode1, anodeDis1, cathode1, cathodeDis1); + printf(" The 2nd nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", anode2, anodeDis2, cathode2, cathodeDis2); + } + +private: + + // the nearest wire + short anode1; + short cathode1; + double anodeDis1; + double cathodeDis1; + + // the 2nd nearest wire + short anode2; + short cathode2; + double anodeDis2; + double cathodeDis2; + + TVector3 trackPos; + TVector3 trackVec; + + const int nWire = 24; + const int wireShift = 3; + const float zLen = 380; //mm + const float radiusA = 37; + const float radiusC = 43; + + double dAngle; + double anodeLength; + double cathodeLength; + + std::vector> An; // the anode wire position vector in space + std::vector> Ca; // the cathode wire position vector in space + + double Distance(TVector3 a1, TVector3 a2, TVector3 b1, TVector3 b2){ + TVector3 na = a1 - a2; + TVector3 nb = b1 - b2; + TVector3 nd = (na.Cross(nb)).Unit(); + return TMath::Abs(nd.Dot(a1-b2)); + } + +}; + +inline void PW::Clear(){ + anode1 = -1; + cathode1 = -1; + anodeDis1 = 999999999; + cathodeDis1 = 999999999; + anode2 = -1; + cathode2 = -1; + anodeDis2 = 999999999; + cathodeDis2 = 999999999; +} + +inline void PW::ConstructGeo(){ + std::pair p1; // anode + std::pair q1; // cathode + + //anode and cathode start at pos-Y axis and count in left-Hand + //anode wire shift is right-hand. + //cathode wire shift is left-hand. + + for(int i = 0; i < nWire; i++ ){ + // Anode rotate right-hand + p1.first.SetXYZ( radiusA * TMath::Cos( TMath::TwoPi() / nWire * (-i) + TMath::PiOver2()), + radiusA * TMath::Sin( TMath::TwoPi() / nWire * (-i) + TMath::PiOver2()), + zLen/2); + p1.second.SetXYZ( radiusA * TMath::Cos( TMath::TwoPi() / nWire * (-i + wireShift) + TMath::PiOver2()), + radiusA * TMath::Sin( TMath::TwoPi() / nWire * (-i + wireShift) + TMath::PiOver2()), + -zLen/2); + An.push_back(p1); + + // Cathod rotate left-hand + q1.first.SetXYZ( radiusC * TMath::Cos( TMath::TwoPi() / nWire * (-i) + TMath::PiOver2()), + radiusC * TMath::Sin( TMath::TwoPi() / nWire * (-i) + TMath::PiOver2()), + zLen/2); + q1.second.SetXYZ( radiusC * TMath::Cos( TMath::TwoPi() / nWire * (- i - wireShift) + TMath::PiOver2()), + radiusC * TMath::Sin( TMath::TwoPi() / nWire * (- i - wireShift) + TMath::PiOver2()), + -zLen/2); + Ca.push_back(q1); + } + + dAngle = wireShift * TMath::TwoPi() / nWire; + anodeLength = TMath::Sqrt( zLen*zLen + TMath::Power(2* radiusA * TMath::Sin(dAngle/2),2) ); + cathodeLength = TMath::Sqrt( zLen*zLen + TMath::Power(2* radiusC * TMath::Sin(dAngle/2),2) ); +} + +inline void PW::FindWireID(TVector3 pos, TVector3 direction, bool verbose ){ + + Clear(); + + double phi = direction.Phi(); + for( int i = 0; i < nWire; i++){ + + double disA = 99999999; + double disC = 99999999; + + double phiS = An[i].first.Phi() - TMath::PiOver4(); + double phiL = An[i].second.Phi() + TMath::PiOver4(); + + // printf("A%2d: %f %f | %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg(), phi * TMath::RadToDeg()); + if( phi > 0 && phiS > phiL ) phiL = phiL + TMath::TwoPi(); + if( phi < 0 && phiS > phiL ) phiS = phiS - TMath::TwoPi(); + + if( phiS < phi && phi < phiL) { + disA = Distance( pos, pos + direction, An[i].first, An[i].second); + if( disA < anodeDis1 ){ + anodeDis1 = disA; + anode1 = i; + } + } + + phiS = Ca[i].second.Phi()- TMath::PiOver4(); + phiL = Ca[i].first.Phi() + TMath::PiOver4(); + // printf("C%2d: %f %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg()); + if( phi > 0 && phiS > phiL ) phiL = phiL + TMath::TwoPi(); + if( phi < 0 && phiS > phiL ) phiS = phiS - TMath::TwoPi(); + + if(phiS < phi && phi < phiL) { + disC = Distance( pos, pos + direction, Ca[i].first, Ca[i].second); + if( disC < cathodeDis1 ){ + cathodeDis1 = disC; + cathode1 = i; + } + } + + if(verbose) printf(" %2d | %8.2f, %8.2f\n", i, disA, disC); + } + + //==== find the 2nd nearest wire + 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); + if( haha1 < haha2){ + anode2 = anode1-1; + anodeDis2 = haha1; + }else{ + anode2 = anode1+1; + anodeDis2 = haha2; + } + + 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); + if( haha1 < haha2){ + cathode2 = cathode1-1; + cathodeDis2 = haha1; + }else{ + cathode2 = cathode1+1; + cathodeDis2 = haha2; + } + + if( verbose ) Print(); +} + +inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose){ + + trackPos = sx3Pos; + + TVector3 n1 = (An[anodeID].first - An[anodeID].second).Cross((sx3Pos - An[anodeID].second)).Unit(); + TVector3 n2 = (Ca[cathodeID].first - Ca[cathodeID].second).Cross((sx3Pos - Ca[cathodeID].second)).Unit(); + + // if the handiness of anode and cathode revered, it should be n2 cross n1 + trackVec = (n2.Cross(n1)).Unit(); + + if( verbose ) printf("Theta, Phi = %f, %f \n", trackVec.Theta() *TMath::RadToDeg(), trackVec.Phi()*TMath::RadToDeg()); + +} + +#endif \ No newline at end of file diff --git a/Armory/ClassSX3.h b/Armory/ClassSX3.h new file mode 100644 index 0000000..48c2807 --- /dev/null +++ b/Armory/ClassSX3.h @@ -0,0 +1,219 @@ +#ifndef ClassSX3_h +#define ClassSX3_h + +#include +#include +#include + +class SX3{ +public: + SX3(){} + ~SX3(){} + + short GetID() const {return id;} + short GetChUp() const {return chUp;} + short GetChDn() const {return chDn;} + short GetChBk() const {return chBk;} + + TVector3 GetHitPos() const {return hitPos;} + + double GetZFrac() const {return zFrac;} // range from -0.5 to 0.5 + + void Clear(); + void ConstructGeo(); + void FindSX3Pos(TVector3 pos, TVector3 direction, bool verbose = false); + void CalSX3Pos(unsigned short ID, unsigned short chUp, unsigned short chDown, unsigned short chBack, float eUp, float eDown); + + double GetNumDet() const {return nSX3;} + double GetWidth() const {return sx3Width;} + double GetLength() const {return sx3Length;} + TVector3 GetDnL(short id) const {return SDn[id].first; } // lower strip ID + TVector3 GetDnH(short id) const {return SDn[id].second; } // higher strip ID + TVector3 GetUpL(short id) const {return SUp[id].first; } // lower strip ID + TVector3 GetUpH(short id) const {return SUp[id].second; } // higher strip ID + + TVector3 GetDnMid(short id) const { return (SDn[id].first + SDn[id].second)*0.5;} + TVector3 GetUpMid(short id) const { return (SUp[id].first + SUp[id].second)*0.5;} + + double GetDetPhi(short id) const { return (SUp[id].second - SUp[id].first).Phi();} + + void Print(){ + if( id == -1 ){ + printf("Did not hit any SX3.\n"); + }else{ + printf("ID: %d, U,D,B: %d %d %d| zFrac : %.2f\n", id, chUp, chDn, chBk, zFrac); + printf("Hit Pos: %.2f, %.2f, %.2f\n", hitPos.X(), hitPos.Y(), hitPos.Z()); + } + } + + // void CalZFrac(){ + // zFrac = (eUp - eDn)/(eUp + eDn); + // } + +private: + + short id; // -1 when no hit + short chUp; + short chDn; + short chBk; + + double zFrac; // from +1 (downstream) to -1 (upstream) + + double eUp; + double eDn; + double eBk; + + TVector3 hitPos; + + const int nSX3 = 12; + const float sx3Radius = 88; + const float sx3Width = 40; + const float sx3Length = 75; + const float sx3Gap = 46; + + std::vector> SDn; // coners of the SX3 0-11, z = mid point + std::vector> SUp; // coners of the SX3 12-23, z = mid point + std::vector SNorml; // normal of the SX3 (outward) + + + std::pair Intersect(TVector3 p1, TVector3 p2, TVector3 q1, TVector3 q2, bool verbose){ + + //see https://nukephysik101.wordpress.com/2023/12/30/intersect-between-2-line-segments/ + //zero all z-component + TVector3 a0 = p1; a0.SetZ(0); + TVector3 a1 = p2; a1.SetZ(0); + + TVector3 b0 = q1; b0.SetZ(0); + TVector3 b1 = q2; b1.SetZ(0); + + double A = ((b0-b1).Cross(a0-a1)).Mag(); + + double h = ((b0-a0).Cross(b1-a0)).Z()/ A; + double k = ((a1-b0).Cross(a0-b0)).Z()/ A; + + if( verbose ) printf(" ----h, k : %f, %f\n", h, k); + + return std::pair(h,k); + } + +}; + +inline void SX3::Clear(){ + id = -1; + chUp = -1; + chDn = -1; + chBk = -1; + zFrac = TMath::QuietNaN(); + + eUp = TMath::QuietNaN(); + eDn = TMath::QuietNaN(); + eBk = TMath::QuietNaN(); +} + +inline void SX3::ConstructGeo(){ + TVector3 sa, sb, sc, sn; + for(int i = 0; i < nSX3; i++){ + sa.SetXYZ( sx3Radius, -sx3Width/2, sx3Gap/2 + sx3Length/2 ); + sb.SetXYZ( sx3Radius, sx3Width/2, sx3Gap/2 + sx3Length/2 ); + + double rot = TMath::TwoPi() / nSX3 * (-i - 0.5) - TMath::PiOver2(); + + sa.RotateZ( rot ); + sb.RotateZ( rot ); + SDn.push_back(std::pair(sa,sb)); + + + sc.SetXYZ( sx3Radius, -sx3Width/2, sx3Gap/2 ); + sc.RotateZ( rot ); + + sn = ((sc-sa).Cross(sb-sa)).Unit(); + SNorml.push_back(sn); + + sa.SetXYZ( sx3Radius, -sx3Width/2, -sx3Gap/2 - sx3Length/2 ); + sb.SetXYZ( sx3Radius, sx3Width/2, -sx3Gap/2 - sx3Length/2 ); + + sa.RotateZ( rot ); + sb.RotateZ( rot ); + SUp.push_back(std::pair(sa,sb)); + } +} + + +inline void SX3::FindSX3Pos(TVector3 pos, TVector3 direction, bool verbose){ + + id = -1; + for( int i = 0 ; i < nSX3; i++){ + + if(verbose) printf(" %d ", i); + std::pair frac = Intersect( pos, pos + direction, SDn[i].first, SDn[i].second, verbose); + + + if( frac.second < 0 || frac.second > 1 ) continue; + hitPos = pos + frac.first * direction; + + double dis = hitPos.Dot(SNorml[i]); + + if(verbose) { + printf("reduced distance : %f\n", dis); + printf(" %d*", (i+1)%nSX3); + Intersect( pos, pos + direction, SDn[(i+1)%nSX3].first, SDn[(i+1)%nSX3].second, verbose); + } + + if( TMath::Abs(dis - sx3Radius) > 0.1 ) continue; + + chDn = 2 * TMath::Floor(frac.second * 4); + chUp = chDn + 1; + + double zPos = hitPos.Z(); + if( (sx3Gap/2 < zPos && zPos < sx3Gap/2 + sx3Length ) || (-sx3Gap/2 - sx3Length < zPos && zPos < -sx3Gap/2 ) ){ + + id = zPos > 0 ? i : i + 12; + + zFrac = zPos > 0 ? (zPos - sx3Gap/2. - sx3Length/2.)/sx3Length : (zPos - ( - sx3Gap/2. - sx3Length/2.) )/sx3Length ; + + chBk = TMath::Floor( (zFrac + 0.5) * 4 ) + 8; + + if( verbose) Print(); + return ; + + }else{ + if( verbose ) printf(" zPos out of sensitive region\n"); + } + } + + if( verbose) Print(); +} + + +inline void SX3::CalSX3Pos(unsigned short ID, unsigned short chUp, unsigned short chDown, unsigned short chBack, float eUp, float eDown){ + + hitPos.Clear(); + + if( (chUp - chDown) != 1 || (chDown % 2) != 0) return ; + + int reducedID = ID % nSX3; + + TVector3 sa, sb; + + if( ID < nSX3 ){ //down + sa = SDn[reducedID].second; + sb = SDn[reducedID].first; + }else{ + sa = SUp[reducedID].second; + sb = SUp[reducedID].first; + } + + hitPos.SetX( (sb.X() - sa.X()) * chUp/8 + sa.X()); + hitPos.SetY( (sb.Y() - sa.Y()) * chUp/8 + sa.Y()); + + if( eUp == 0 || eDown == 0 ){ + hitPos.SetZ( sa.Z() + (2*(chBk - 7)-1) * sx3Length / 8 ); + }else{ + double frac = (eUp - eDown)/(eUp + eDown); // from +1 (downstream) to -1 (upstream) + double zPos = sa.Z() + sx3Length * frac/2; + hitPos.SetZ( zPos ); + } + +} + +#endif \ No newline at end of file diff --git a/script.C b/script.C index b867f6f..f24508a 100644 --- a/script.C +++ b/script.C @@ -273,8 +273,13 @@ void script(TString fileName = "", int maxEvent = -1){ haha->DrawTrack(pos, dir, true); - PC pc = haha->FindWireID(pos, dir, true); - SX3 sx3 = haha->FindSX3Pos(pos, dir, true); + PW pc; + pc.ConstructGeo(); + pc.FindWireID(pos, dir, true); + + SX3 sx3; + sx3.ConstructGeo(); + sx3.FindSX3Pos(pos, dir, true); // haha->CalTrack(sx3.hitPos, wireID.first, wireID.second, true);