rtang/ANASEN_analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

modify anasenMS.cpp for the last change

Browse Source
This commit is contained in:
Ryan Tang 2024-02-02 15:40:43 -05:00
parent 50b95b33ed
commit e442fc1102
4 changed files with 107 additions and 90 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

76
Armory/ClassAnasen.h
View File

@ -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<short, short> wireID = pw.GetNearestID();
std::pair<short, short> 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);

3
Armory/ClassPW.h
View File

@ -96,6 +96,9 @@ inline void PW::Clear(){
cathode2 = -1;
anodeDis2 = 999999999;
cathodeDis2 = 999999999;
An.clear();
Ca.clear();
}
inline void PW::ConstructGeo(){

59
Armory/ClassSX3.h
View File

@ -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<std::pair<TVector3,TVector3>> SDn; // coners of the SX3 0-11, z = mid point
std::vector<std::pair<TVector3,TVector3>> SUp; // coners of the SX3 12-23, z = mid point
std::vector<TVector3> 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<double, double> 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 );
}

55
Armory/anasenMS.cpp
View File

@ -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("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<int, int> 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<int, int> 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;
}