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

modified: Analyzer.C implemented basic trackreconstruction

Browse Source 
modified:   Armory/ANASEN_model.C changed qqq radii
	modified:   Armory/ClassPW.h implemented basic trackreconstruction
This commit is contained in:
Vignesh Sitaraman 2025-02-27 10:34:41 -05:00
parent 9225620426
commit 39e8f41ab1
4 changed files with 154 additions and 185 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.gitignore vendored
View File

@ -10,3 +10,6 @@ AnasenMS
data/
data_proton/
Analyzer_C_ACLiC_dict0713aaa966_dictContent.h
.gitignore
Analyzer_C_ACLiC_dict5411fecd5c_dictUmbrella.h

103
Analyzer.C
View File

@ -39,7 +39,6 @@ TH1F *hCat4An;
TH1F *hCat0An;
TH1F *hAnodehits;
TH2F *hNosvAe;
TH2F *hUncorrCAN;
int padID = 0;
@ -92,7 +91,7 @@ void Analyzer::Begin(TTree * /*tree*/)
hsx3VpcE->SetNdivisions(-612, "x");
hsx3VpcE->SetNdivisions(-12, "y");
// hZProj = new TH1F("hZProj", "Z Projection", 1200, -600, 600);
hZProj = new TH1F("hZProj", "Z Projection", 1200, -600, 600);
hPCZProj = new TH1F("hPCZProj", "PC Z Projection", 600, -300, 300);
hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400, 0, 16000, 400, 0, 20000);
@ -100,7 +99,6 @@ void Analyzer::Begin(TTree * /*tree*/)
hCat0An = new TH1F("hCat0An", "Number of Cathodes without Anode", 24, 0, 24);
hAnodehits = new TH1F("hAnodehits", "Number of Anode hits", 24, 0, 24);
hNosvAe = new TH2F("hnosvAe", "Number of Cathodes/Anode vs Anode Energy", 20, 0, 20, 400, 0, 16000);
hUncorrCAN = new TH2F("hUncorrCAn", "Uncorrelated Cathodes/Anode", 24, 0, 24, 24, 0, 24);
// for (int i = 0; i < 24; i++)
// {
// TString histName = Form("hAnodeVsCathode_%d", i);
@ -182,9 +180,9 @@ Bool_t Analyzer::Process(Long64_t entry)
for (int i = 0; i < sx3.multi; i++)
{
ID.push_back(std::pair<int, int>(sx3.id[i], i));
hsx3IndexVE->Fill(sx3.index[i], sx3.e[i]);
if (sx3.e[i] > 1000 && sx3.e[i] < 2200)
if (sx3.e[i] > 100)
{
sx3ecut = true;
}
@ -279,15 +277,18 @@ Bool_t Analyzer::Process(Long64_t entry)
// hitPos.Print();
}
}
qqqEcut = false;
// //======================= QQQ
qqqEcut = false;
for (int i = 0; i < qqq.multi; i++)
{
// for( int j = 0; j < pc.multi; j++){
// if(pc.index[j]==4){
hqqqIndexVE->Fill(qqq.index[i], qqq.e[i]);
// }
if (qqq.e[i] > 1600 && qqq.e[i] < 3000)
// printf("QQQ ID : %d, ch : %d, e : %d \n", qqq.id[i], qqq.ch[i], qqq.e[i]);
if (qqq.e[i] > 100)
{
qqqEcut = true;
}
@ -307,7 +308,7 @@ Bool_t Analyzer::Process(Long64_t entry)
{
hqqqVpcE->Fill(qqq.e[i], pc.e[k]);
// hpcIndexVE->Fill( pc.index[i], pc.e[i] );
hqqqVpcIndex->Fill(qqq.index[i], pc.index[j]);
hqqqVpcIndex->Fill(qqq.index[i], pc.index[k]);
}
// }
}
@ -329,7 +330,6 @@ Bool_t Analyzer::Process(Long64_t entry)
chRing = qqq.ch[i];
chWedge = qqq.ch[j] - 16;
}
// printf(" ID : %d , chWedge : %d, chRing : %d \n", qqq.id[i], chWedge, chRing);
double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
@ -356,8 +356,8 @@ Bool_t Analyzer::Process(Long64_t entry)
pwinstance.ConstructGeo();
Coord Crossover[24][24][2];
TVector3 a, c, diff, an, cn;
double a2, an2, ac, c2, cn2, adiff, cdiff, denom, alpha, beta;
TVector3 a, c, diff;
double a2, ac, c2, adiff, cdiff, denom, alpha, beta;
int index = 0;
for (int i = 0; i < pwinstance.An.size(); i++)
{
@ -377,8 +377,6 @@ Bool_t Analyzer::Process(Long64_t entry)
ac = a.Dot(c);
adiff = a.Dot(diff);
cdiff = c.Dot(diff);
an2 = an.Dot(an);
cn2 = cn.Dot(cn);
denom = a2 * c2 - ac * ac;
alpha = (ac * cdiff - c2 * adiff) / denom;
beta = (a2 * cdiff - ac * adiff) / denom;
@ -386,7 +384,8 @@ Bool_t Analyzer::Process(Long64_t entry)
Crossover[i][j][0].x = pwinstance.An[i].first.X() + alpha * a.X();
Crossover[i][j][0].y = pwinstance.An[i].first.Y() + alpha * a.Y();
Crossover[i][j][0].z = pwinstance.An[i].first.Z() + alpha * a.Z();
if(Crossover[i][j][0].z <-190 || Crossover[i][j][0].z > 190) Crossover[i][j][0].z = 9999999;
if (Crossover[i][j][0].z < -190 || Crossover[i][j][0].z > 190)
Crossover[i][j][0].z = 9999999;
// placeholder variable Crossover[i][j][2].x has nothing to do with the geometry of the crossover and is being used to store the alpha value-
//-so that it can be used to sort "good" hits later
@ -502,7 +501,12 @@ Bool_t Analyzer::Process(Long64_t entry)
std::sort(cathodeHits.begin(), cathodeHits.end(), [](const std::pair<int, double> &a, const std::pair<int, double> &b)
{ return a.second > b.second; });
if (anodeHits.size() >= 1 && cathodeHits.size() >1)
bool SiPCflag;
corrcatMax.clear();
if (anodeHits.size() >= 1 && cathodeHits.size() > 1)
{
// if (((TMath::TanH(hitPos.Y() / hitPos.X())) > (TMath::TanH(a.Y() / a.X()) - TMath::PiOver4())) || ((TMath::TanH(hitPos.Y() / hitPos.X())) < (TMath::TanH(a.Y() / a.X()) + TMath::PiOver4())))
{
for (const auto &anode : anodeHits)
@ -528,58 +532,29 @@ Bool_t Analyzer::Process(Long64_t entry)
// std::cout << " Anode iD : " << aIDMax << " Energy : " << aEMax << std::endl;
// printf("aID : %d, aE : %f, cE : %f\n", aID, aE, cE);
corrcatMax.clear();
for (const auto &cathode : cathodeHits)
{
cID = cathode.first;
cE = cathode.second;
// std::cout << "Cathode ID : " << cID << " Energy : " << cE << std::endl;
// if (cE > cEMax)
// {
// cIDnextMax = cIDMax;
// cEnextMax = cEMax;
// cEMax = cE;
// cIDMax = cID;
// }
// if (cE > cEnextMax && cE < cEMax)
// {
// cEnextMax = cE;
// cIDnextMax = cID;
// }
hAVCcoin->Fill(aIDMax, cID);
// This section of code is used to find the cathodes are correlated with the max and next max anodes, as well as to figure out if there are any common cathodes
// the anodes are correlated with the cathodes +/-3 from the anode number in the reverse order
// bool corr = false;
hAVCcoin->Fill(aIDMax, cID);
for (int j = -4; j < 3; j++)
{
if ((aIDMax + 24 + j) % 24 == 23 - cID) /* the 23-cID is used to accomodate for the fact that the order of the cathodes was reversed relative top the physical geometry */
// if(Crossover[aIDMax][cID][0].z != 9999999)
// for (int j = -4; j < 3; j++)
// if ((aIDMax + 24 + j) % 24 == 23 - cID) /* the 23-cID is used to accomodate for the fact that the order of the cathodes was reversed relative top the physical geometry */
if (Crossover[aIDMax][cID][0].z != 9999999)
{
corrcatMax.push_back(std::pair<int, double>(cID, cE));
cESum += cE;
// printf("Max Anode : %d Correlated Cathode : %d Anode Energy : %f z value : %f \n", aIDMax, cID, cESum, Crossover[aIDMax][cID][1].z /*prints alpha*/);
// std::cout << " Cathode iD : " << cID << " Energy : " << cE << std::endl;
cESum += cE;
// corr = true;
}
// if ((aIDnextMax + 24 + j) % 24 == cID)
// {
// corrcatnextMax.push_back(std::pair<int, double>(cathode.first, cathode.second));
// std::sort(corrcatMax.begin(), corrcatMax.end(), [](const std::pair<int, double> &a, const std::pair<int, double> &b)S
// { return a.second > b.second; });
// }
// for (int k = 0; k < 5; k++)
// {
// if ((aIDMax + 24 + j) % 24 == (aIDnextMax + 24 + k) % 24)
// {
// commcat.push_back(std::pair<int, double>(cathode.first, cathode.second));
// }
// }
}
// if (!corr)
// {
// hUncorrCAN->Fill(aIDMax, cID);
// }
}
}
TVector3 anodeIntersection;
@ -601,8 +576,19 @@ Bool_t Analyzer::Process(Long64_t entry)
{
printf("Warning: No valid cathode hits to correlate with anode %d! \n", aIDMax);
}
// printf("EventID : %llu, Max Anode : %d Cathode: %d PC X and Y : (%f, %f) \n", entry, aIDMax, cID, Crossover[aIDMax][cID][0].x, Crossover[aIDMax][cID][0].y);
// for (int i = 0; i < sx3.multi; i++)
// {
// printf("EventID : %llu, HitPos X, Y, Z: %f %f %f SX3ID : %d %d \n", entry, hitPos.X(), hitPos.Y(), hitPos.Z(), sx3.id[i], sx3.ch[i]);
// }
// for (int i = 0; i < qqq.multi; i++)
// {
// printf("Max Anode : %d Cathode: %d PC X and Y : %f %f \n", aIDMax, cID, Crossover[aIDMax][cID][0].x, Crossover[aIDMax][cID][0].y);
// printf("HitPos X, Y, Z, QQQID : %f %f %f %d \n", hitPos.X(), hitPos.Y(), hitPos.Z(), qqq.id[i]);
// }
}
if(qqqEcut) anodeIntersection = TVector3(x, y, z);
anodeIntersection = TVector3(x, y, z);
// std::cout << "Anode Intersection " << anodeIntersection.Z() << " " << x << " " << y << " " << z << std::endl;
}
@ -659,11 +645,10 @@ Bool_t Analyzer::Process(Long64_t entry)
hCat0An->Fill(cathodeHits.size());
}
// if (HitNonZero)
// {
// pw_contr.CalTrack(hitPos, aID, cID);
// hZProj->Fill(pw_contr.GetZ0());
// }
if (HitNonZero)
{
pw_contr.CalTrack2(hitPos, anodeIntersection);
hZProj->Fill(pw_contr.GetZ0());
}
// ########################################################### Track constrcution

4
Armory/ANASEN_model.C
View File

@ -103,8 +103,8 @@ void ANASEN_model(int anodeID1 = -1, int anodeID2 = -1, int cathodeID1 = -1, int
new TGeoRotation("rot1", 360/nSX3 * (i + 0.5), 0., 0.)));
}
const int qqqR1 = 10;
const int qqqR2 = 50;
const int qqqR1 = 50;
const int qqqR2 = 100;
TGeoVolume *qqq = geom->MakeTubs("qqq", Al, qqqR1, qqqR2, 0.5, 5, 85);
qqq->SetLineColor(7);
for( int i = 0; i < 4; i++){

47
Armory/ClassPW.h
View File

@ -82,7 +82,7 @@ public:
void ConstructGeo();
void FindWireID(TVector3 pos, TVector3 direction, bool verbose = false);
void CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose = false);
void CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA = 0, double sigmaC = 0, bool verbose = false);
void CalTrack2(TVector3 sx3Pos, TVector3 anodeInt, bool verbose = false);
void Print()
{
@ -283,41 +283,22 @@ inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbo
printf("Theta, Phi = %f, %f \n", trackVec.Theta() * TMath::RadToDeg(), trackVec.Phi() * TMath::RadToDeg());
}
inline void PW::CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA, double sigmaC, bool verbose)
inline void PW::CalTrack2(TVector3 siPos, TVector3 anodeInt, bool verbose)
{
trackPos = sx3Pos;
double p1 = TMath::Abs(hitInfo.nearestDist.first + gRandom->Gaus(0, sigmaA));
double p2 = TMath::Abs(hitInfo.nextNearestDist.first + gRandom->Gaus(0, sigmaA));
double fracA = p1 / (p1 + p2);
short anodeID1 = hitInfo.nearestWire.first;
short anodeID2 = hitInfo.nextNearestWire.first;
TVector3 shiftA1 = (An[anodeID2].first - An[anodeID1].first) * fracA;
TVector3 shiftA2 = (An[anodeID2].second - An[anodeID1].second) * fracA;
double q1 = TMath::Abs(hitInfo.nearestDist.second + gRandom->Gaus(0, sigmaC));
double q2 = TMath::Abs(hitInfo.nextNearestDist.second + gRandom->Gaus(0, sigmaC));
double fracC = q1 / (q1 + q2);
short cathodeID1 = hitInfo.nearestWire.second;
short cathodeID2 = hitInfo.nextNearestWire.second;
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;
TVector3 c1 = Ca[cathodeID1].first + shiftC1;
TVector3 c2 = Ca[cathodeID1].second + shiftC2;
TVector3 n1 = (a1 - a2).Cross((sx3Pos - a2)).Unit();
TVector3 n2 = (c1 - c2).Cross((sx3Pos - c2)).Unit();
// if the handiness of anode and cathode revered, it should be n2 cross n1
trackVec = (n2.Cross(n1)).Unit();
float mx, my;
double z;
mx = siPos.X() / (siPos.X() - anodeInt.X());
my = siPos.Y() / (siPos.Y() - anodeInt.Y());
z=siPos.Z() + mx * (anodeInt.Z() - siPos.Z());
// if (mx == my)
{
trackVec=TVector3(0,0,z);
}
if (verbose)
printf("Theta, Phi = %f, %f \n", trackVec.Theta() * TMath::RadToDeg(), trackVec.Phi() * TMath::RadToDeg());
printf("X slope = %f and Y slope = %f \n", mx, my);
}
inline double PW::GetZ0()
@ -328,7 +309,7 @@ inline double PW::GetZ0()
double rho = TMath::Sqrt(x * x + y * y);
double theta = trackVec.Theta();
return trackPos.Z() - rho / TMath::Tan(theta);
return trackVec.Z();
}
#endif