From e442fc11025f67bc0d09beef2bbca41d6ebf3a24 Mon Sep 17 00:00:00 2001 From: Ryan Tang Date: Fri, 2 Feb 2024 15:40:43 -0500 Subject: [PATCH] modify anasenMS.cpp for the last change --- Armory/ClassAnasen.h | 76 +++++++++++++++++++++++++------------------- Armory/ClassPW.h | 3 ++ Armory/ClassSX3.h | 59 ++++++++++++++++++---------------- Armory/anasenMS.cpp | 59 +++++++++++++++++----------------- 4 files changed, 107 insertions(+), 90 deletions(-) diff --git a/Armory/ClassAnasen.h b/Armory/ClassAnasen.h index 763001c..3fb7d29 100644 --- a/Armory/ClassAnasen.h +++ b/Armory/ClassAnasen.h @@ -30,10 +30,15 @@ public: int sx3ID = -1, bool DrawQQQ = false ); + + + PW * GetPW() {return pw;} + SX3 * GetSX3() {return sx3;} + private: - PW pw; - SX3 sx3; + PW * pw; + SX3 * sx3; const float qqqR1 = 50; const float qqqR2 = 100; @@ -56,6 +61,9 @@ private: //!============================================== inline ANASEN::ANASEN(){ + pw = new PW(); + sx3 = new SX3(); + CalGeometry(); geom = nullptr; @@ -67,12 +75,16 @@ inline ANASEN::~ANASEN(){ delete geom; + delete pw; + delete sx3; + } + //!============================================== inline void ANASEN::CalGeometry(){ - sx3.ConstructGeo(); - pw.ConstructGeo(); + sx3->ConstructGeo(); + pw->ConstructGeo(); } @@ -111,24 +123,24 @@ inline void ANASEN::Construct3DModel(int anodeID1, int anodeID2, int cathodeID1, worldBox->AddNode(axisZ, 1, new TGeoTranslation(0, 0, 5)); //.......... convert to wire center dimensions - TGeoVolume *pcA = geom->MakeTube("tub1", Al, 0, 0.01, pw.GetAnodeLength()/2); + TGeoVolume *pcA = geom->MakeTube("tub1", Al, 0, 0.01, pw->GetAnodeLength()/2); pcA->SetLineColor(4); int startID = 0; - int endID = pw.GetNumWire() - 1; + int endID = pw->GetNumWire() - 1; if( anodeID1 >= 0 && anodeID2 >= 0 ){ startID = anodeID1; endID = anodeID2; if( anodeID1 > anodeID2 ) { - endID = pw.GetNumWire() + anodeID2; + endID = pw->GetNumWire() + anodeID2; } } for( int i = startID; i <= endID; i++){ - TVector3 a = pw.GetAnodneMid(i); - double wireTheta = pw.GetAnodeTheta(i) * TMath::RadToDeg(); - double wirePhi = pw.GetAnodePhi(i) * 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(), @@ -136,24 +148,24 @@ inline void ANASEN::Construct3DModel(int anodeID1, int anodeID2, int cathodeID1, new TGeoRotation("rot1", wirePhi, wireTheta, 0.))); } - TGeoVolume *pcC = geom->MakeTube("tub2", Al, 0, 0.01, pw.GetCathodeLength()/2); + TGeoVolume *pcC = geom->MakeTube("tub2", Al, 0, 0.01, pw->GetCathodeLength()/2); pcC->SetLineColor(6); startID = 0; - endID = pw.GetNumWire() - 1; + endID = pw->GetNumWire() - 1; if( cathodeID1 >= 0 && cathodeID2 >= 0 ){ startID = cathodeID1; endID = cathodeID2; if( cathodeID1 > cathodeID2 ) { - endID = pw.GetNumWire() + cathodeID2; + endID = pw->GetNumWire() + cathodeID2; } } for( int i = startID; i <= endID; i++){ - TVector3 a = pw.GetCathodneMid(i); - double wireTheta = pw.GetCathodeTheta(i) * TMath::RadToDeg(); - double wirePhi = pw.GetCathodePhi(i) * 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(), @@ -161,14 +173,14 @@ inline void ANASEN::Construct3DModel(int anodeID1, int anodeID2, int cathodeID1, new TGeoRotation("rot1", wirePhi , wireTheta, 0.))); } - TGeoVolume * sx3Det = geom->MakeBox("box", Al, 0.1, sx3.GetWidth()/2, sx3.GetLength()/2); + TGeoVolume * sx3Det = geom->MakeBox("box", Al, 0.1, sx3->GetWidth()/2, sx3->GetLength()/2); sx3Det->SetLineColor(kGreen+3); - for( int i = 0; i < sx3.GetNumDet(); i++){ + for( int i = 0; i < sx3->GetNumDet(); i++){ if( sx3ID != -1 && i != sx3ID ) continue; - 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; + 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(), @@ -206,10 +218,10 @@ inline void ANASEN::DrawAnasen(int anodeID1, int anodeID2, int cathodeID1, int c inline void ANASEN::DrawTrack(TVector3 pos, TVector3 direction, bool drawEstimatedTrack){ - pw.FindWireID(pos, direction); - sx3.FindSX3Pos(pos, direction); + pw->FindWireID(pos, direction); + sx3->FindSX3Pos(pos, direction); - std::pair wireID = pw.GetNearestID(); + std::pair wireID = pw->GetNearestID(); Construct3DModel(wireID.first, wireID.first, wireID.second, wireID.second, -1, false); @@ -225,19 +237,19 @@ 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.GetID() >= 0 ){ + if( sx3->GetID() >= 0 ){ TGeoVolume * hit = geom->MakeSphere("hitpos", 0, 0, 3); hit->SetLineColor(kRed); - TVector3 hitPos = sx3.GetHitPos(); + TVector3 hitPos = sx3->GetHitPos(); worldBox->AddNode(hit, 2, new TGeoCombiTrans( hitPos.X(), hitPos.Y(), hitPos.Z(), new TGeoRotation("rotA", 0, 0, 0.))); if( drawEstimatedTrack ){ - pw.CalTrack(hitPos, wireID.first, wireID.second, true); + pw->CalTrack(hitPos, wireID.first, wireID.second, true); - double thetaDeduce = pw.GetTrackTheta() * TMath::RadToDeg(); - double phiDeduce = pw.GetTrackPhi() * 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); @@ -255,12 +267,12 @@ inline void ANASEN::DrawTrack(TVector3 pos, TVector3 direction, bool drawEstima inline void ANASEN::DrawDeducedTrack(TVector3 sx3Pos, int anodeID, int cathodeID){ - pw.CalTrack(sx3Pos, anodeID, cathodeID); + pw->CalTrack(sx3Pos, anodeID, cathodeID); Construct3DModel(anodeID, anodeID, cathodeID, cathodeID, -1, false); - double theta = pw.GetTrackTheta() * TMath::RadToDeg(); - double phi = pw.GetTrackPhi() * 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); diff --git a/Armory/ClassPW.h b/Armory/ClassPW.h index d0273ad..55b88dc 100644 --- a/Armory/ClassPW.h +++ b/Armory/ClassPW.h @@ -96,6 +96,9 @@ inline void PW::Clear(){ cathode2 = -1; anodeDis2 = 999999999; cathodeDis2 = 999999999; + + An.clear(); + Ca.clear(); } inline void PW::ConstructGeo(){ diff --git a/Armory/ClassSX3.h b/Armory/ClassSX3.h index 48c2807..3f29e59 100644 --- a/Armory/ClassSX3.h +++ b/Armory/ClassSX3.h @@ -7,7 +7,7 @@ class SX3{ public: - SX3(){} + SX3(){Clear();}; ~SX3(){} short GetID() const {return id;} @@ -24,9 +24,9 @@ public: 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;} + double GetNumDet() const {return numDet;} + double GetWidth() const {return width;} + double GetLength() const {return length;} 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 @@ -52,6 +52,12 @@ public: private: + const int numDet = 12; + const float radius = 88; + const float width = 40; + const float length = 75; + const float gap = 46; + short id; // -1 when no hit short chUp; short chDn; @@ -65,12 +71,6 @@ private: 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) @@ -108,29 +108,32 @@ inline void SX3::Clear(){ eUp = TMath::QuietNaN(); eDn = TMath::QuietNaN(); eBk = TMath::QuietNaN(); + + SDn.clear(); + SUp.clear(); } 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(); + for(int i = 0; i < numDet; i++){ + sa.SetXYZ( radius, -width/2, gap/2 + length/2 ); + sb.SetXYZ( radius, width/2, gap/2 + length/2 ); + + double rot = TMath::TwoPi() / numDet * (-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.SetXYZ( radius, -width/2, gap/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.SetXYZ( radius, -width/2, -gap/2 - length/2 ); + sb.SetXYZ( radius, width/2, -gap/2 - length/2 ); sa.RotateZ( rot ); sb.RotateZ( rot ); @@ -142,7 +145,7 @@ inline void SX3::ConstructGeo(){ inline void SX3::FindSX3Pos(TVector3 pos, TVector3 direction, bool verbose){ id = -1; - for( int i = 0 ; i < nSX3; i++){ + for( int i = 0 ; i < numDet; i++){ if(verbose) printf(" %d ", i); std::pair frac = Intersect( pos, pos + direction, SDn[i].first, SDn[i].second, verbose); @@ -155,21 +158,21 @@ inline void SX3::FindSX3Pos(TVector3 pos, TVector3 direction, bool verbose){ 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); + printf(" %d*", (i+1)%numDet); + Intersect( pos, pos + direction, SDn[(i+1)%numDet].first, SDn[(i+1)%numDet].second, verbose); } - if( TMath::Abs(dis - sx3Radius) > 0.1 ) continue; + if( TMath::Abs(dis - radius) > 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 ) ){ + if( (gap/2 < zPos && zPos < gap/2 + length ) || (-gap/2 - length < zPos && zPos < -gap/2 ) ){ id = zPos > 0 ? i : i + 12; - zFrac = zPos > 0 ? (zPos - sx3Gap/2. - sx3Length/2.)/sx3Length : (zPos - ( - sx3Gap/2. - sx3Length/2.) )/sx3Length ; + zFrac = zPos > 0 ? (zPos - gap/2. - length/2.)/length : (zPos - ( - gap/2. - length/2.) )/length ; chBk = TMath::Floor( (zFrac + 0.5) * 4 ) + 8; @@ -191,11 +194,11 @@ inline void SX3::CalSX3Pos(unsigned short ID, unsigned short chUp, unsigned shor if( (chUp - chDown) != 1 || (chDown % 2) != 0) return ; - int reducedID = ID % nSX3; + int reducedID = ID % numDet; TVector3 sa, sb; - if( ID < nSX3 ){ //down + if( ID < numDet ){ //down sa = SDn[reducedID].second; sb = SDn[reducedID].first; }else{ @@ -207,10 +210,10 @@ inline void SX3::CalSX3Pos(unsigned short ID, unsigned short chUp, unsigned shor 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 ); + hitPos.SetZ( sa.Z() + (2*(chBk - 7)-1) * length / 8 ); }else{ double frac = (eUp - eDown)/(eUp + eDown); // from +1 (downstream) to -1 (upstream) - double zPos = sa.Z() + sx3Length * frac/2; + double zPos = sa.Z() + length * frac/2; hitPos.SetZ( zPos ); } diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index 125d82e..ecb3bce 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -45,7 +45,9 @@ int main(int argc, char **argv){ int nExA = ExAList.size(); int nEx = ExList.size(); - ANASEN anasen; + ANASEN * anasen = new ANASEN(); + SX3 * sx3 = anasen->GetSX3(); + PW * pw = anasen->GetPW(); TString saveFileName = "SimAnasen.root"; printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data()); @@ -92,12 +94,12 @@ int main(int argc, char **argv){ tree->Branch("aID", &anodeID, "anodeID/I"); tree->Branch("cID", &cathodeID, "cathodeID/I"); - int sx3ID, sx3Up, sx3Down, sx3Back; + int sx3ID, sx3Up, sx3Dn, sx3Bk; double sx3ZFrac; - tree->Branch("sx3ID", &sx3ID, "sx3ID/I"); - tree->Branch("sx3Up", &sx3Up, "sx3Up/I"); - tree->Branch("sx3Down", &sx3Down, "sx3Down/I"); - tree->Branch("sx3Back", &sx3Back, "sx3Back/I"); + tree->Branch("sx3ID", &sx3ID, "sx3ID/I"); + tree->Branch("sx3Up", &sx3Up, "sx3Up/I"); + tree->Branch("sx3Dn", &sx3Dn, "sx3Dn/I"); + tree->Branch("sx3Bk", &sx3Bk, "sx3Bk/I"); tree->Branch("sx3ZFrac", &sx3ZFrac, "sx3ZFrac/D"); double reTheta, rePhi; @@ -150,38 +152,35 @@ int main(int argc, char **argv){ dir.SetTheta(thetab * TMath::DegToRad()); dir.SetPhi(phib * TMath::DegToRad()); - std::pair wireID = anasen.FindWireID(vertex, dir, false); - SX3 sx3 = anasen.FindSX3Pos(vertex, dir, false); + + pw->FindWireID(vertex, dir, false); + sx3->FindSX3Pos(vertex, dir, false); + + std::pair wireID = pw->GetNearestID(); anodeID = wireID.first; cathodeID = wireID.second; - sx3ID = sx3.id; - if( sx3.id >= 0 ){ - sx3Up = sx3.chUp; - sx3Down = sx3.chDown; - sx3Back = sx3.chBack; - sx3ZFrac = sx3.zFrac; + sx3ID = sx3->GetID(); + if( sx3ID >= 0 ){ + sx3Up = sx3->GetChUp(); + sx3Dn = sx3->GetChDn(); + sx3Bk = sx3->GetChBk(); + sx3ZFrac = sx3->GetZFrac(); - sx3X = sx3.hitPos.X(); - sx3Y = sx3.hitPos.Y(); - sx3Z = sx3.hitPos.Z(); + sx3X = sx3->GetHitPos().X(); + sx3Y = sx3->GetHitPos().Y(); + sx3Z = sx3->GetHitPos().Z(); - // for( int i = 0; i < 12; i++){ - // sx3Index[i] = -1; - // if( i == sx3Up ) sx3Index[i] = sx3ID * 12 + sx3Up; - // if( i == sx3Down ) sx3Index[i] = sx3ID * 12 + sx3Down; - // if( i == sx3Back ) sx3Index[i] = sx3ID * 12 + sx3Back; - // } - anasen.CalTrack(sx3.hitPos, wireID.first, wireID.second, false); + pw->CalTrack(sx3->GetHitPos(), wireID.first, wireID.second, false); - reTheta = anasen.GetTrackTheta() * TMath::RadToDeg(); - rePhi = anasen.GetTrackPhi() * TMath::RadToDeg(); + reTheta = pw->GetTrackTheta() * TMath::RadToDeg(); + rePhi = pw->GetTrackPhi() * TMath::RadToDeg(); }else{ sx3Up = -1; - sx3Down = -1; - sx3Back = -1; + sx3Dn = -1; + sx3Bk = -1; sx3ZFrac = TMath::QuietNaN(); sx3X = TMath::QuietNaN(); @@ -196,8 +195,6 @@ int main(int argc, char **argv){ rePhi = TMath::QuietNaN(); } - - tree->Fill(); //#################################################################### Timer @@ -224,6 +221,8 @@ int main(int argc, char **argv){ printf("=============== done. saved as %s. count(hit==1) : %d\n", saveFileName.Data(), count); + delete anasen; + return 0; } \ No newline at end of file