modified: .gitignore
modified: Armory/ClassPW.h new geo + crossover moved to PW.h modified: MakeVertex.C added plots new file: anasen_fem/README.md new file: anasen_fem/clean.sh new file: anasen_fem/junk/wires.py new file: anasen_fem/junk/wires2d_test.sif new file: anasen_fem/junk/wires_gmsh.py changed wire radius from Sudarsan's version new file: anasen_fem/junk/wires_gmsh_bc.py new file: anasen_fem/paraview_plotter.py new file: anasen_fem/run.py new file: anasen_fem/scalars.dat new file: anasen_fem/scalars.dat.names new file: anasen_fem/wires2d.sif new file: anasen_fem/wires_gmsh2d_bc.py
This commit is contained in:
parent
65af941b39
commit
ec9d25b048
1
.gitignore
vendored
1
.gitignore
vendored
|
|
@ -8,6 +8,7 @@ EventBuilder*
|
|||
*.err
|
||||
*.seq
|
||||
*.png
|
||||
*.gif
|
||||
Mapper
|
||||
AnasenMS
|
||||
|
||||
|
|
|
|||
|
|
@ -156,7 +156,6 @@ inline void PW::ConstructGeo()
|
|||
double k = TMath::TwoPi()/24.; //48 solder thru holes, wires in every other one
|
||||
double offset_a1 = -6*k-3*k;
|
||||
double offset_c1 = -4*k -2*k - TMath::TwoPi()/48; //correct for a half-turn
|
||||
#include "../scratch/testing.h"
|
||||
double offset_a2 = offset_a1+wireShift*k;
|
||||
double offset_c2 = offset_c1-wireShift*k;
|
||||
|
||||
|
|
|
|||
585
MakeVertex.C
585
MakeVertex.C
|
|
@ -1,5 +1,13 @@
|
|||
#define MakeVertex_cxx
|
||||
|
||||
Int_t colors[40] = {
|
||||
kBlack, kRed, kGreen, kBlue, kYellow, kMagenta, kCyan, kOrange,
|
||||
kSpring, kTeal, kAzure, kViolet, kPink, kGray, kWhite,
|
||||
kRed + 2, kGreen + 2, kBlue + 2, kYellow + 2, kMagenta + 2, kCyan + 2, kOrange + 2,
|
||||
kSpring + 2, kTeal + 2, kAzure + 2, kViolet + 2, kPink + 2,
|
||||
kRed - 7, kGreen - 7, kBlue - 7, kYellow - 7, kMagenta - 7, kCyan - 7, kOrange - 7,
|
||||
kSpring - 7, kTeal - 7, kAzure - 7, kViolet - 7, kPink - 7, kGray + 2};
|
||||
|
||||
#include "MakeVertex.h"
|
||||
#include "Armory/ClassPW.h"
|
||||
#include "Armory/HistPlotter.h"
|
||||
|
|
@ -10,7 +18,12 @@
|
|||
#include <TCanvas.h>
|
||||
#include <TMath.h>
|
||||
#include <TBranch.h>
|
||||
#include "TVector3.h"
|
||||
#include <TVector3.h>
|
||||
#include <TGraph2D.h>
|
||||
#include <TView.h>
|
||||
#include <TPolyLine3D.h>
|
||||
#include <TPolyMarker3D.h>
|
||||
#include <TH3D.h>
|
||||
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
|
|
@ -21,8 +34,32 @@
|
|||
#include <utility>
|
||||
#include <algorithm>
|
||||
|
||||
bool realtime = true;
|
||||
const double source_vertex = 53;
|
||||
const double qqq_z = 100.0;
|
||||
const double anode_gain = 1.5146e-5; // channels --> MeV
|
||||
|
||||
TApplication *app = NULL;
|
||||
TH1F *hha = NULL, *hhc = NULL;
|
||||
TH3D *frame = NULL;
|
||||
TCanvas *can1 = NULL, *can2 = NULL;
|
||||
|
||||
TPolyLine3D *pla[24] = {NULL};
|
||||
TPolyLine3D *plc[24] = {NULL};
|
||||
TPolyLine3D *qqqw[16][4] = {NULL};
|
||||
TGraph2D *qqqg = NULL, *crossoverg = NULL, *guessg = NULL;
|
||||
|
||||
double z_to_crossover_rho(double z)
|
||||
{
|
||||
return 9.20645e-5 * z * z + 34.1973;
|
||||
}
|
||||
|
||||
double z_to_crossover_rho_cathode(double z)
|
||||
{
|
||||
return 9.20645e-5 * z * z + 34.1973;
|
||||
}
|
||||
|
||||
// Global instances
|
||||
PW pw_contr;
|
||||
PW pwinstance;
|
||||
TVector3 hitPos;
|
||||
double qqqenergy, qqqtimestamp;
|
||||
|
|
@ -42,15 +79,13 @@ public:
|
|||
double Time2 = -1;
|
||||
int Anodech = -1;
|
||||
int Cathodech = -1;
|
||||
int multi1 = -1, multi2 = -1;
|
||||
};
|
||||
|
||||
// Calibration globals
|
||||
const int MAX_QQQ = 4;
|
||||
const int MAX_RING = 16;
|
||||
const int MAX_WEDGE = 16;
|
||||
const double qqqpos = 100.0;
|
||||
const double vertexpos = 53.4;
|
||||
const double pcrad = 37.0;
|
||||
double qqqGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
|
||||
bool qqqGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
|
||||
double qqqCalib[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
|
||||
|
|
@ -76,8 +111,9 @@ void MakeVertex::Begin(TTree * /*tree*/)
|
|||
{
|
||||
TString option = GetOption();
|
||||
plotter = new HistPlotter("Analyzer_SX3.root", "TFILE");
|
||||
pw_contr.ConstructGeo();
|
||||
pwinstance.ConstructGeo();
|
||||
// if (gROOT->IsBatch())
|
||||
realtime = false;
|
||||
|
||||
// ---------------------------------------------------------
|
||||
// 1. CRITICAL FIX: Initialize PC Arrays to Default (Raw)
|
||||
|
|
@ -88,40 +124,6 @@ void MakeVertex::Begin(TTree * /*tree*/)
|
|||
pcIntercept[i] = 0.0; // Default intercept = 0
|
||||
}
|
||||
|
||||
// Calculate Crossover Geometry ONCE
|
||||
TVector3 a, c, diff;
|
||||
double a2, ac, c2, adiff, cdiff, denom, alpha;
|
||||
|
||||
for (size_t i = 0; i < pwinstance.An.size(); i++)
|
||||
{
|
||||
a = pwinstance.An[i].first - pwinstance.An[i].second;
|
||||
|
||||
for (size_t j = 0; j < pwinstance.Ca.size(); j++)
|
||||
{
|
||||
c = pwinstance.Ca[j].first - pwinstance.Ca[j].second;
|
||||
diff = pwinstance.An[i].first - pwinstance.Ca[j].first;
|
||||
a2 = a.Dot(a);
|
||||
c2 = c.Dot(c);
|
||||
ac = a.Dot(c);
|
||||
adiff = a.Dot(diff);
|
||||
cdiff = c.Dot(diff);
|
||||
denom = a2 * c2 - ac * ac;
|
||||
alpha = (ac * cdiff - c2 * adiff) / denom;
|
||||
|
||||
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 || (i + j) % 24 == 12)
|
||||
{
|
||||
Crossover[i][j][0].z = 9999999;
|
||||
}
|
||||
|
||||
Crossover[i][j][1].x = alpha;
|
||||
Crossover[i][j][1].y = 0;
|
||||
}
|
||||
}
|
||||
|
||||
// Load PC Calibrations
|
||||
std::ifstream inputFile("slope_intercept_results.dat");
|
||||
if (inputFile.is_open())
|
||||
|
|
@ -223,6 +225,61 @@ void MakeVertex::Begin(TTree * /*tree*/)
|
|||
}
|
||||
infile.close();
|
||||
}
|
||||
if (realtime)
|
||||
{
|
||||
can1 = new TCanvas("wireindex", "c1", 0, 0, 640, 480);
|
||||
can2 = new TCanvas("3d", "c2", 650, 0, 640, 480);
|
||||
can1->cd();
|
||||
// can2->SetFillColor(30);
|
||||
frame = new TH3D("frame", "frame", 1000, -100, 100, 1000, -100, 100, 1000, -200, 200);
|
||||
hha = new TH1F("hha", "Anode Ecal vs wire#", 48, -12, 36);
|
||||
hhc = new TH1F("hhc", "Cathode Ecal vs wire#", 48, -12, 36);
|
||||
hha->SetLineColor(kRed);
|
||||
hha->GetYaxis()->SetRangeUser(0, 16384);
|
||||
hha->GetXaxis()->SetTitle("press any key, interrupt/refresh or double click to continue..");
|
||||
hha->Draw();
|
||||
hhc->Draw("SAME");
|
||||
can1->Modified();
|
||||
can1->Update();
|
||||
can1->BuildLegend();
|
||||
|
||||
can2->cd();
|
||||
frame->Draw();
|
||||
for (int i = 0; i < 24; i++)
|
||||
{
|
||||
plc[i] = new TPolyLine3D(2);
|
||||
pla[i] = new TPolyLine3D(2);
|
||||
pla[i]->SetPoint(0, pwinstance.An[i].first.X(), pwinstance.An[i].first.Y(), pwinstance.An[i].first.Z());
|
||||
pla[i]->SetPoint(1, pwinstance.An[i].second.X(), pwinstance.An[i].second.Y(), pwinstance.An[i].second.Z());
|
||||
plc[i]->SetPoint(0, pwinstance.Ca[i].first.X(), pwinstance.Ca[i].first.Y(), pwinstance.Ca[i].first.Z());
|
||||
plc[i]->SetPoint(1, pwinstance.Ca[i].second.X(), pwinstance.Ca[i].second.Y(), pwinstance.Ca[i].second.Z());
|
||||
plc[i]->SetLineStyle(kDotted);
|
||||
pla[i]->SetLineStyle(kDotted);
|
||||
pla[i]->SetLineWidth(1.);
|
||||
plc[i]->SetLineWidth(1.);
|
||||
plc[i]->Draw("same");
|
||||
pla[i]->Draw("same");
|
||||
plc[i]->SetLineColor(colors[i]);
|
||||
pla[i]->SetLineColor(colors[i]);
|
||||
}
|
||||
crossoverg = new TGraph2D(1);
|
||||
crossoverg->SetName("crossoverg");
|
||||
crossoverg->SetMarkerStyle(20);
|
||||
crossoverg->SetMarkerColor(kBlue + 3);
|
||||
qqqg = new TGraph2D(1);
|
||||
qqqg->SetName("qqqg");
|
||||
qqqg->SetMarkerColor(kRed);
|
||||
qqqg->SetMarkerStyle(42);
|
||||
|
||||
crossoverg->SetPoint(0, 0, 0, 0);
|
||||
qqqg->SetPoint(0, 0, 0, qqq_z);
|
||||
crossoverg->Draw("P same");
|
||||
qqqg->Draw("P same");
|
||||
|
||||
can2->Modified();
|
||||
can2->Update();
|
||||
}
|
||||
|
||||
std::cout << "aaa" << std::endl;
|
||||
}
|
||||
|
||||
|
|
@ -254,53 +311,76 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
pc.CalIndex();
|
||||
|
||||
std::vector<Event> sx3Events;
|
||||
// if(sx3.multi>1) {
|
||||
// std::array<sx3det,12> Fsx3;
|
||||
// //std::cout << "-----" << std::endl;
|
||||
// for(int i=0; i<sx3.multi; i++) {
|
||||
// if(sx3.id[i]>=12) continue;
|
||||
// int id = sx3.id[i];
|
||||
// if(sx3.ch[i]>=8) {
|
||||
// int sx3ch=sx3.ch[i]-8;
|
||||
// sx3ch=(sx3ch+3)%4;
|
||||
// if(sx3ch==0 || sx3ch==3) continue;
|
||||
// float value=sx3.e[i];
|
||||
// int gch = sx3.id[i]*4+(sx3.ch[i]-8);
|
||||
// Fsx3.at(id).fillevent("BACK",sx3ch,value);
|
||||
// Fsx3.at(id).ts = static_cast<double>(sx3.t[i]);
|
||||
// plotter->Fill2D("sx3backs_raw",100,0,100,800,0,4096,gch,sx3.e[i]);
|
||||
// } else {
|
||||
// int sx3ch=sx3.ch[i]/2;
|
||||
// double value=sx3.e[i];
|
||||
// if(sx3.ch[i]%2==0) {
|
||||
// Fsx3.at(id).fillevent("FRONT_L",sx3ch,value*sx3RightGain[id][sx3ch]);
|
||||
// } else {
|
||||
// Fsx3.at(id).fillevent("FRONT_R",sx3ch,value);
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// for(int id=0; id<12; id++) {
|
||||
// Fsx3.at(id).validate();
|
||||
// auto det = Fsx3.at(id);
|
||||
// bool no_charge_sharing_strict = det.valid_front_chans.size()==1 && det.valid_back_chans.size()==1;
|
||||
// if(det.valid) {
|
||||
// //std::cout << det.frontEL << " " << det.frontEL*sx3RightGain[id][det.stripF] << std::endl;
|
||||
// plotter->Fill2D("be_vs_x_sx3_id_"+std::to_string(id)+"_f"+std::to_string(det.stripF)+"_b"+std::to_string(det.stripB),200,-1,1,800,0,8192,
|
||||
// det.frontX,det.backE,"evsx");
|
||||
// //std::cout << sx3BackGain[id][det.stripF][det.stripB] << " " << sx3FrontGain[id][det.stripF] << std::endl;
|
||||
// plotter->Fill2D("matched_be_vs_x_sx3_id_"+std::to_string(id)+"_f"+std::to_string(det.stripF)+"_"+std::to_string(id*4+det.stripF),200,-30,30,800,0,8192,
|
||||
// det.frontX*sx3FrontGain[id][det.stripF]+sx3FrontOffset[id][det.stripF],det.backE*sx3BackGain[id][det.stripF][det.stripB],"evsx_matched");
|
||||
// //plotter->Fill2D("fe_vs_x_sx3_id_"+std::to_string(id)+"_f"+std::to_string(det.stripF)+"_"+std::to_string(det.stripB),200,-1,1,800,0,4096,det.frontX,det.backE,"evsx");
|
||||
// plotter->Fill2D("l_vs_r_sx3_id_"+std::to_string(id)+"_f"+std::to_string(det.stripF),800,0,4096,800,0,4096,det.frontEL,det.frontER,"l_vs_r");
|
||||
// }
|
||||
// if(det.valid && (id ==9 || id==7 || id == 1 || id==3) && det.stripF!=DEFAULT_NULL && det.stripB!=DEFAULT_NULL) {
|
||||
// double z = det.frontX*sx3FrontGain[id][det.stripF]+sx3FrontOffset[id][det.stripF];
|
||||
// double backE = det.backE*sx3BackGain[id][det.stripF][det.stripB];
|
||||
// Event sx3ev(TVector3(0,0,z),backE,-1,det.ts,-1,det.stripB+4*id,det.stripF+4*id);
|
||||
// sx3Events.push_back(sx3ev);
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
if (sx3.multi > 1)
|
||||
{
|
||||
std::array<sx3det, 24> Fsx3;
|
||||
// std::cout << "-----" << std::endl;
|
||||
for (int i = 0; i < sx3.multi; i++)
|
||||
{
|
||||
int id = sx3.id[i];
|
||||
// if(id>=12) continue;
|
||||
if (sx3.ch[i] >= 8)
|
||||
{
|
||||
int sx3ch = sx3.ch[i] - 8;
|
||||
sx3ch = (sx3ch + 3) % 4;
|
||||
if (sx3ch == 0 || sx3ch == 3)
|
||||
continue;
|
||||
float value = sx3.e[i];
|
||||
int gch = sx3.id[i] * 4 + (sx3.ch[i] - 8);
|
||||
Fsx3.at(id).fillevent("BACK", sx3ch, value);
|
||||
Fsx3.at(id).ts = static_cast<double>(sx3.t[i]);
|
||||
plotter->Fill2D("sx3backs_raw", 100, 0, 100, 800, 0, 4096, gch, sx3.e[i]);
|
||||
}
|
||||
else
|
||||
{
|
||||
int sx3ch = sx3.ch[i] / 2;
|
||||
double value = sx3.e[i];
|
||||
if (sx3.ch[i] % 2 == 0)
|
||||
{
|
||||
Fsx3.at(id).fillevent("FRONT_L", sx3ch, value * sx3RightGain[id][sx3ch]);
|
||||
}
|
||||
else
|
||||
{
|
||||
Fsx3.at(id).fillevent("FRONT_R", sx3ch, value);
|
||||
}
|
||||
}
|
||||
}
|
||||
for (int id = 0; id < 24; id++)
|
||||
{
|
||||
// std::cout << id << " " << Fsx3.at(id).valid_front_chans.size() << " " << Fsx3.at(id).valid_back_chans.size() << std::endl;;
|
||||
try
|
||||
{
|
||||
Fsx3.at(id).validate();
|
||||
}
|
||||
catch (std::exception exc)
|
||||
{
|
||||
std::cout << "oops! anyway" << std::endl;
|
||||
continue;
|
||||
}
|
||||
auto det = Fsx3.at(id);
|
||||
bool no_charge_sharing_strict = det.valid_front_chans.size() == 1 && det.valid_back_chans.size() == 1;
|
||||
if (det.valid)
|
||||
{
|
||||
// std::cout << det.frontEL << " " << det.frontEL*sx3RightGain[id][det.stripF] << std::endl;
|
||||
plotter->Fill2D("be_vs_x_sx3_id_" + std::to_string(id) + "_f" + std::to_string(det.stripF) + "_b" + std::to_string(det.stripB), 200, -1, 1, 800, 0, 8192,
|
||||
det.frontX, det.backE, "evsx");
|
||||
// std::cout << sx3BackGain[id][det.stripF][det.stripB] << " " << sx3FrontGain[id][det.stripF] << std::endl;
|
||||
plotter->Fill2D("matched_be_vs_x_sx3_id_" + std::to_string(id) + "_f" + std::to_string(det.stripF), 200, -30, 30, 800, 0, 8192,
|
||||
det.frontX * sx3FrontGain[id][det.stripF] + sx3FrontOffset[id][det.stripF], det.backE * sx3BackGain[id][det.stripF][det.stripB], "evsx_matched");
|
||||
// plotter->Fill2D("fe_vs_x_sx3_id_"+std::to_string(id)+"_f"+std::to_string(det.stripF)+"_"+std::to_string(det.stripB),200,-1,1,800,0,4096,det.frontX,det.backE,"evsx");
|
||||
plotter->Fill2D("l_vs_r_sx3_id_" + std::to_string(id) + "_f" + std::to_string(det.stripF), 800, 0, 4096, 800, 0, 4096, det.frontEL, det.frontER, "l_vs_r");
|
||||
}
|
||||
if (det.valid && (id == 9 || id == 7 || id == 1 || id == 3) && det.stripF != DEFAULT_NULL && det.stripB != DEFAULT_NULL)
|
||||
{
|
||||
double z = det.frontX * sx3FrontGain[id][det.stripF] + sx3FrontOffset[id][det.stripF];
|
||||
double backE = det.backE * sx3BackGain[id][det.stripF][det.stripB];
|
||||
double beta_n = 15.0 + TMath::ATan2((2 * det.stripF - 3) * 40.30, 8.0 * 88.0 * TMath::Cos(15.0 * M_PI / 180.0)) * 180. / M_PI; // how much to add per strip to the starting position
|
||||
double phi_n = ((-id + 0.5) * 30 + beta_n) * M_PI / 180.; // starting-position phi + strip contribution
|
||||
Event sx3ev(TVector3(88.0 * TMath::Cos(phi_n), 88.0 * TMath::Sin(phi_n), z), backE, -1, det.ts, -1, det.stripB + 4 * id, det.stripF + 4 * id);
|
||||
sx3Events.push_back(sx3ev);
|
||||
}
|
||||
}
|
||||
}
|
||||
// return kTRUE;
|
||||
// QQQ Processing
|
||||
|
||||
|
|
@ -333,7 +413,6 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
// if(qqq.multi>=3) std::cout << "-----" << std::endl;
|
||||
for (int i = 0; i < qqq.multi; i++)
|
||||
{
|
||||
// if(qqq.multi>=3) std::cout << std::setprecision(16) << "qqq"<< qqq.id[i] << " " << std::string(qqq.ch[i]/16?"ring":"wedge") << qqq.ch[i]%16 << " " << qqq.e[i] << " " << qqq.t[i] - qqq.t[0] << std::endl;
|
||||
if (qqq.ch[i] / 16)
|
||||
{
|
||||
if (qvecr[qqq.id[i]].find(qqq.ch[i]) != qvecr[qqq.id[i]].end())
|
||||
|
|
@ -350,6 +429,8 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
}
|
||||
|
||||
bool PCQQQTimeCut = false;
|
||||
bool PCAQQQTimeCut = false;
|
||||
bool PCCQQQTimeCut = false;
|
||||
for (int i = 0; i < qqq.multi; i++)
|
||||
{
|
||||
plotter->Fill2D("QQQ_Index_Vs_Energy", 16 * 8, 0, 16 * 8, 2000, 0, 16000, qqq.index[i], qqq.e[i], "hRawQQQ");
|
||||
|
|
@ -423,20 +504,21 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
continue;
|
||||
// if(eRingMeV<4.0 || eWedgeMeV<4.0) continue;
|
||||
|
||||
double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
|
||||
// double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5); old method
|
||||
double theta = 2 * TMath::Pi() * (-qqq.id[i] * 16 + (15 - chWedge) + 0.5) / (16 * 4);
|
||||
double rho = 50. + (50. / 16.) * (chRing + 0.5); //"?"
|
||||
// z used to be 75+30+23=128
|
||||
// we found a 12mm shift towards the vertex later --> 116
|
||||
Event qqqevent(TVector3(rho * TMath::Cos(theta), rho * TMath::Sin(theta), qqqpos), eRingMeV, eWedgeMeV, tRing, tWedge, chRing + qqq.id[i] * 16, chWedge + qqq.id[i] * 16);
|
||||
Event qqqeventr(TVector3(rho * TMath::Cos(theta), rho * TMath::Sin(theta), qqqpos), eRing, eWedge, tRing, tWedge, chRing + qqq.id[i] * 16, chWedge + qqq.id[i] * 16);
|
||||
QQQ_Events.push_back(qqqevent);
|
||||
QQQ_Events_Raw.push_back(qqqeventr);
|
||||
Event qqqevent(TVector3(rho * TMath::Cos(theta), rho * TMath::Sin(theta), qqq_z), eRingMeV, eWedgeMeV, tRing, tWedge, chRing + qqq.id[i] * 16, chWedge + qqq.id[i] * 16);
|
||||
Event qqqeventr(TVector3(rho * TMath::Cos(theta), rho * TMath::Sin(theta), qqq_z), eRing, eWedge, tRing, tWedge, chRing + qqq.id[i] * 16, chWedge + qqq.id[i] * 16);
|
||||
if (qqq.id[i] >= 1)
|
||||
{
|
||||
QQQ_Events.push_back(qqqevent);
|
||||
QQQ_Events_Raw.push_back(qqqeventr);
|
||||
}
|
||||
|
||||
plotter->Fill2D("QQQCartesianPlot", 200, -100, 100, 200, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hCalQQQ");
|
||||
plotter->Fill2D("QQQCartesianPlot" + std::to_string(qqq.id[i]), 200, -100, 100, 200, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hCalQQQ");
|
||||
if (PCQQQTimeCut)
|
||||
{
|
||||
plotter->Fill2D("PC_XY_Projection_QQQ_TimeCut" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hPCQQQ");
|
||||
}
|
||||
plotter->Fill2D("PC_XY_Projection_QQQ" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hPCQQQ");
|
||||
}
|
||||
else
|
||||
|
|
@ -463,23 +545,29 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
// if (tRing - static_cast<double>(pc.t[k]) < -150 && tRing - static_cast<double>(pc.t[k]) > -450) // 27Al
|
||||
// if (tRing - static_cast<double>(pc.t[k]) < -70 && tRing - static_cast<double>(pc.t[k]) > -150) // 17F
|
||||
{
|
||||
PCQQQTimeCut = true;
|
||||
PCAQQQTimeCut = true;
|
||||
}
|
||||
}
|
||||
|
||||
if (pc.index[k] >= 24 && pc.e[k] > 50)
|
||||
if (pc.index[k] >= 24 && pc.e[k] > 10)
|
||||
{
|
||||
if (tRing - static_cast<double>(pc.t[k]) < -200)
|
||||
PCCQQQTimeCut = true;
|
||||
plotter->Fill2D("Timing_Difference_QQQ_PC_Cathode", 500, -2000, 2000, 16, 0, 16, tRing - static_cast<double>(pc.t[k]), chRing, "hTiming");
|
||||
}
|
||||
} // end of pc k loop
|
||||
|
||||
if (!HitNonZero)
|
||||
{
|
||||
double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
|
||||
// double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5); old method
|
||||
double theta = 2 * TMath::Pi() * (-qqq.id[i] * 16 + (15 - chWedge) + 0.5) / (16 * 4);
|
||||
double rho = 50. + (50. / 16.) * (chRing + 0.5); //"?"
|
||||
|
||||
double x = rho * TMath::Cos(theta);
|
||||
double y = rho * TMath::Sin(theta);
|
||||
hitPos.SetXYZ(x, y, (qqqpos));
|
||||
hitPos.SetXYZ(x, y, (qqq_z));
|
||||
if (realtime)
|
||||
qqqg->SetPoint(0, hitPos.X(), hitPos.Y(), hitPos.Z());
|
||||
qqqenergy = eRingMeV;
|
||||
qqqtimestamp = tRing;
|
||||
HitNonZero = true;
|
||||
|
|
@ -488,19 +576,19 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
} // j loop end
|
||||
} // i loop end
|
||||
|
||||
plotter->Fill1D("QQQ_Multiplicity", 10, 0, 10, qqqCount, "hRawQQQ");
|
||||
PCQQQTimeCut = PCAQQQTimeCut && PCCQQQTimeCut;
|
||||
|
||||
/*if(QQQ_Events.size()>=1) {
|
||||
std::cout<< " ---->" << std::endl;
|
||||
for(auto qe: QQQ_Events) {
|
||||
std::cout << qe.ch1/16 << " " <<qe.ch2/16 << " " << qe.ch1%16 << " "<< qe.ch2%16 << " " << qe.Energy1 << " " << qe.Energy2 << " " << std::endl;
|
||||
}
|
||||
}*/
|
||||
plotter->Fill1D("QQQ_Multiplicity", 10, 0, 10, qqqCount, "hRawQQQ");
|
||||
|
||||
typedef std::unordered_map<int, std::tuple<int, double, double>> WireEvent; // this stores nearest neighbour wire events, or a 'cluster'
|
||||
WireEvent aWireEvents, cWireEvents; // naming for book keeping
|
||||
aWireEvents.clear();
|
||||
aWireEvents.reserve(24);
|
||||
if (realtime)
|
||||
{
|
||||
hha->Reset();
|
||||
hhc->Reset();
|
||||
}
|
||||
|
||||
// PC Gain Matching and Filling
|
||||
double anodeT = -99999;
|
||||
|
|
@ -509,7 +597,7 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
int cathodeIndex = -1;
|
||||
for (int i = 0; i < pc.multi; i++)
|
||||
{
|
||||
if (pc.e[i] > 50)
|
||||
if (pc.e[i] > 20)
|
||||
{
|
||||
plotter->Fill2D("PC_Index_Vs_Energy", 48, 0, 48, 2000, 0, 30000, pc.index[i], static_cast<double>(pc.e[i]), "hRawPC");
|
||||
}
|
||||
|
|
@ -529,12 +617,16 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
anodeT = static_cast<double>(pc.t[i]);
|
||||
anodeIndex = pc.index[i];
|
||||
aWireEvents[pc.index[i]] = std::tuple(pc.index[i], pc.e[i], static_cast<double>(pc.t[i]));
|
||||
if (realtime)
|
||||
hha->SetBinContent(hha->FindFixBin(anodeIndex), pc.e[i]);
|
||||
}
|
||||
else
|
||||
{
|
||||
cathodeT = static_cast<double>(pc.t[i]);
|
||||
cathodeIndex = pc.index[i] - 24;
|
||||
cWireEvents[pc.index[i] - 24] = std::tuple(pc.index[i] - 24, pc.e[i], static_cast<double>(pc.t[i]));
|
||||
if (realtime)
|
||||
hhc->SetBinContent(hhc->FindFixBin(cathodeIndex), pc.e[i]);
|
||||
}
|
||||
|
||||
if (anodeT != -99999 && cathodeT != 99999)
|
||||
|
|
@ -575,7 +667,9 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
double aESum = 0;
|
||||
double cESum = 0;
|
||||
double aEMax = 0;
|
||||
double cEMax = 0;
|
||||
int aIDMax = 0;
|
||||
int cIDMax = 0;
|
||||
|
||||
for (int i = 0; i < pc.multi; i++)
|
||||
{
|
||||
|
|
@ -592,10 +686,10 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
}
|
||||
}
|
||||
|
||||
std::sort(anodeHits.begin(), anodeHits.end(), [](std::pair<int, double> a, std::pair<int, double> b)
|
||||
{ return a.first < b.first; });
|
||||
std::sort(cathodeHits.begin(), cathodeHits.end(), [](std::pair<int, double> a, std::pair<int, double> b)
|
||||
{ return a.first < b.first; });
|
||||
// std::sort(anodeHits.begin(), anodeHits.end(), [](std::pair<int, double> a, std::pair<int, double> b)
|
||||
// { return a.first < b.first; });
|
||||
// std::sort(cathodeHits.begin(), cathodeHits.end(), [](std::pair<int, double> a, std::pair<int, double> b)
|
||||
// { return a.first < b.first; });
|
||||
|
||||
// clusters = collection of (collection of wires) where each wire is (index, energy, timestamp)
|
||||
std::vector<std::vector<std::tuple<int, double, double>>> aClusters = pwinstance.Make_Clusters(aWireEvents);
|
||||
|
|
@ -608,15 +702,17 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
{
|
||||
// if (aCluster.size() <= 1 && cCluster.size() <= 1)
|
||||
// continue;
|
||||
if (aCluster.size() <= 1 && cCluster.size() == 0)
|
||||
continue;
|
||||
auto [crossover, alpha, apSumE, cpSumE, apMaxE, cpMaxE, apTSMaxE, cpTSMaxE] = pwinstance.FindCrossoverProperties(aCluster, cCluster);
|
||||
if (alpha != 9999999 && apSumE != -1)
|
||||
{
|
||||
// Event PCEvent(crossover,apMaxE,cpMaxE,apTSMaxE,cpTSMaxE);
|
||||
// Event PCEvent(crossover,apSumE,cpSumE,apTSMaxE,cpTSMaxE);
|
||||
// Event PCEvent(crossover, apSumE, cpMaxE, apTSMaxE, cpTSMaxE); // run12 shows cathode-max and anode-sum provide best dE signals.
|
||||
Event PCEvent(crossover, apSumE, cpMaxE, apTSMaxE, cpTSMaxE, aCluster.size(), cCluster.size()); // changed to include cluster size info --VS
|
||||
// Event PCEvent(crossover, apSumE, cpMaxE, apTSMaxE, cpTSMaxE, std::get<0>(aCluster.back()), std::get<0>(cCluster.back()), aCluster.size(), cCluster.size()); // changed to include cluster size info and channel numbers--VS
|
||||
Event PCEvent(crossover, apSumE, cpMaxE, apTSMaxE, cpTSMaxE); // run12 shows cathode-max and anode-sum provide best dE signals.
|
||||
// std::cout << apSumE << " " << crossover.Perp() << " " << apMaxE << " " << apTSMaxE << std::endl;
|
||||
// PCEvent.multi1=aCluster.size();
|
||||
// PCEvent.multi2=cCluster.size();
|
||||
PC_Events.push_back(PCEvent);
|
||||
sumE_AC.push_back(std::pair(apSumE, cpSumE));
|
||||
}
|
||||
|
|
@ -625,8 +721,33 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
if (QQQ_Events.size() && PC_Events.size())
|
||||
plotter->Fill2D("PCEv_vs_QQQEv", 20, 0, 20, 20, 0, 20, QQQ_Events.size(), PC_Events.size());
|
||||
|
||||
plotter->Fill2D("ac_vs_cc", 20, 0, 20, 20, 0, 20, aClusters.size(), cClusters.size(), "wiremult");
|
||||
for (auto cluster : aClusters)
|
||||
{
|
||||
plotter->Fill1D("aClusters" + std::to_string(aClusters.size()), 20, -5, 15, cluster.size(), "wiremult");
|
||||
}
|
||||
for (auto cluster : cClusters)
|
||||
{
|
||||
plotter->Fill1D("cClusters" + std::to_string(cClusters.size()), 20, -5, 15, cluster.size(), "wiremult");
|
||||
}
|
||||
|
||||
if (cClusters.size() && aClusters.size())
|
||||
{
|
||||
plotter->Fill2D("ac_vs_cc_ign0", 20, 0, 20, 20, 0, 20, aClusters.size(), cClusters.size(), "wiremult");
|
||||
}
|
||||
|
||||
for (auto pcevent : PC_Events)
|
||||
{
|
||||
if (aClusters.size() == 1 && cClusters.size() == 1)
|
||||
{
|
||||
// plotter->Fill1D("pcz_a"+std::to_string(aClusters.at(0).size())+"_c"+std::to_string(cClusters.at(0).size()),800,-200,200,pcevent.pos.Z(),"wiremult");
|
||||
std::string detid = "_+_";
|
||||
if (sx3Events.size())
|
||||
detid = "+sx3";
|
||||
if (QQQ_Events.size())
|
||||
detid = "+qqq";
|
||||
plotter->Fill1D("pcz_a" + std::to_string(aClusters.at(0).size()) + "_c" + std::to_string(cClusters.at(0).size()) + detid, 800, -200, 200, pcevent.pos.Z(), "wiremult");
|
||||
}
|
||||
for (auto sx3event : sx3Events)
|
||||
{
|
||||
plotter->Fill1D("dt_pcA_sx3B" + std::to_string(sx3event.ch2), 640, -2000, 2000, sx3event.Time1 - pcevent.Time1);
|
||||
|
|
@ -636,11 +757,11 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
plotter->Fill2D("dE_E_Cathodesx3B", 400, 0, 10, 800, 0, 10000, sx3event.Energy1 * 0.001, pcevent.Energy2);
|
||||
double sx3z = sx3event.pos.Z() + (75.0 / 2.0) - 3.0; // w.r.t target origin at 90 for run12
|
||||
double sx3rho = 88.0; // approximate barrel radius
|
||||
double sx3theta = TMath::ATan2(sx3rho, sx3z - vertexpos);
|
||||
double pczguess = pcrad / TMath::Tan(sx3theta) + vertexpos;
|
||||
plotter->Fill2D("pcz_vs_sx3pczguess", 300, 0, 200, 150, 0, 200, pczguess, pcevent.pos.Z());
|
||||
plotter->Fill2D("pcz_vs_sx3pczguess" + std::to_string(sx3event.ch2), 300, 0, 200, 150, 0, 200, pczguess, pcevent.pos.Z());
|
||||
plotter->Fill2D("pcz_vs_sx3z", 300, 0, 200, 150, 0, 200, sx3z, pcevent.pos.Z());
|
||||
double sx3theta = TMath::ATan2(sx3rho, sx3z - source_vertex);
|
||||
double pczguess = 37.0 / TMath::Tan(sx3theta) + source_vertex;
|
||||
plotter->Fill2D("pcz_vs_sx3pczguess", 300, -178, 178, 150, 0, 200, pczguess, pcevent.pos.Z());
|
||||
plotter->Fill2D("pcz_vs_sx3pczguess" + std::to_string(sx3event.ch2), 300, -178, 178, 150, 0, 200, pczguess, pcevent.pos.Z());
|
||||
plotter->Fill2D("pcz_vs_sx3z", 300, 0, 178, 300, -200, 200, sx3z, pcevent.pos.Z());
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -652,18 +773,62 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
// continue;
|
||||
if (aCluster.size() == 1 && cCluster.size() == 1)
|
||||
{
|
||||
plotter->Fill2D("AnodeE_vs_CathodeE_TC" + std::to_string(PCQQQTimeCut) + "_a" + std::to_string(std::get<0>(aCluster.back())) + "c" + std::to_string(std::get<0>(cCluster.back())), 800, 0, 20000, 800, 0, 7000, std::get<1>(aCluster.back()), std::get<1>(cCluster.back()), "AvC");
|
||||
// plotter->Fill2D("AnodeE_vs_CathodeE_TC" + std::to_string(PCQQQTimeCut) + "_a" + std::to_string(std::get<0>(aCluster.back())) + "c" + std::to_string(std::get<0>(cCluster.back())), 800, 0, 20000, 800, 0, 7000, std::get<1>(aCluster.back()), std::get<1>(cCluster.back()), "AvC");
|
||||
plotter->Fill2D("AnodeE_vs_CathodeE_TC" + std::to_string(PCQQQTimeCut), 800, 0, 20000, 800, 0, 7000, std::get<1>(aCluster.back()), std::get<1>(cCluster.back()), "AvC");
|
||||
}
|
||||
else if (aCluster.size() == 1 && cCluster.size() == 2)
|
||||
{
|
||||
plotter->Fill2D("CCh1_vsCCh2", 24, 0, 24, 24, 0, 24, std::get<0>(cCluster.back()), std::get<0>(cCluster.front()), "AvC");
|
||||
if (std::get<1>(cCluster.back()) + std::get<1>(cCluster.front()) < 3400)
|
||||
{
|
||||
plotter->Fill2D("CCh1_vsCCh2_gated", 24, 0, 24, 24, 0, 24, std::get<0>(cCluster.back()), std::get<0>(cCluster.front()), "AvC");
|
||||
|
||||
if (std::get<1>(cCluster.back()) > std::get<1>(cCluster.front()))
|
||||
{
|
||||
plotter->Fill2D("C1vsC2_gated", 400, 0, 8000, 400, 0, 8000, std::get<1>(cCluster.back()), std::get<1>(cCluster.front()), "AvC");
|
||||
}
|
||||
else if (std::get<1>(cCluster.back()) < std::get<1>(cCluster.front()))
|
||||
{
|
||||
plotter->Fill2D("C1vsC2_gated", 400, 0, 8000, 400, 0, 8000, std::get<1>(cCluster.front()), std::get<1>(cCluster.back()), "AvC");
|
||||
}
|
||||
}
|
||||
plotter->Fill2D("AnodeE_vs_CathodeESum_TC" + std::to_string(PCQQQTimeCut), 800, 0, 20000, 800, 0, 14000, std::get<1>(aCluster.back()), std::get<1>(cCluster.back()) + std::get<1>(cCluster.front()), "AvC");
|
||||
plotter->Fill2D("C1vsC2", 800, 0, 8000, 800, 0, 8000, std::get<1>(cCluster.back()), std::get<1>(cCluster.front()), "AvC");
|
||||
// if (std::get<1>(cCluster.back()) > std::get<1>(cCluster.front()))
|
||||
|
||||
plotter->Fill2D("C1vsC2", 400, 0, 8000, 400, 0, 8000, std::get<1>(cCluster.front()), std::get<1>(cCluster.back()), "AvC");
|
||||
plotter->Fill2D("C1vsC2_normA", 1000, 0, 1, 1000, 0, 1, std::get<1>(cCluster.front()) / std::get<1>(aCluster.back()), std::get<1>(cCluster.back()) / std::get<1>(aCluster.back()), "AvC");
|
||||
plotter->Fill2D("C1vsC2_normCsum", 1000, 0, 1, 1000, 0, 1, std::get<1>(cCluster.front()) /( std::get<1>(cCluster.back()) + std::get<1>(cCluster.front())), std::get<1>(cCluster.back())/( std::get<1>(cCluster.back()) + std::get<1>(cCluster.front())), "AvC");
|
||||
plotter->Fill2D("C1vsC2_normA_TC" + std::to_string(PCQQQTimeCut), 1000, 0, 1, 1000, 0, 1, std::get<1>(cCluster.front()) / std::get<1>(aCluster.back()), std::get<1>(cCluster.back()) / std::get<1>(aCluster.back()), "AvC");
|
||||
plotter->Fill2D("C1vsC2_TC" + std::to_string(PCQQQTimeCut), 400, 0, 8000, 400, 0, 8000, std::get<1>(cCluster.front()), std::get<1>(cCluster.back()), "AvC");
|
||||
|
||||
for (auto qqqevent : QQQ_Events)
|
||||
{
|
||||
plotter->Fill2D("qqqER_2Cathode_dESum", 800, 0, 10, 800, 0, 14000, qqqevent.Energy1, std::get<1>(cCluster.back()) + std::get<1>(cCluster.front()), "AvC");
|
||||
plotter->Fill2D("qqqER_AnodeE", 800, 0, 10, 800, 0, 14000, qqqevent.Energy1, std::get<1>(aCluster.back()), "AvC");
|
||||
}
|
||||
}
|
||||
else if (aCluster.size() == 2 && cCluster.size() == 1)
|
||||
{
|
||||
plotter->Fill2D("ACh1_vsACh2", 24, 0, 24, 24, 0, 24, std::get<0>(aCluster.back()), std::get<0>(aCluster.front()), "AvC");
|
||||
if (std::get<1>(aCluster.back()) + std::get<1>(aCluster.front()) < 6800)
|
||||
{
|
||||
plotter->Fill2D("ACh1_vsACh2_gated", 24, 0, 24, 24, 0, 24, std::get<0>(aCluster.back()), std::get<0>(aCluster.front()), "AvC");
|
||||
// if (std::get<1>(aCluster.back()) > std::get<1>(aCluster.front()))
|
||||
{
|
||||
plotter->Fill2D("A1vsA2_gated", 400, 0, 20000, 400, 0, 20000, std::get<1>(aCluster.back()), std::get<1>(aCluster.front()), "AvC");
|
||||
}
|
||||
}
|
||||
plotter->Fill2D("AnodeESum_vs_CathodeE_TC" + std::to_string(PCQQQTimeCut) + "_a" + std::to_string(std::get<0>(aCluster.back())) + "c" + std::to_string(std::get<0>(cCluster.back())), 800, 0, 30000, 800, 0, 7000, std::get<1>(aCluster.back()) + std::get<1>(aCluster.front()), std::get<1>(cCluster.back()), "AvC");
|
||||
plotter->Fill2D("AnodeESum_vs_CathodeE_TC" + std::to_string(PCQQQTimeCut), 800, 0, 30000, 800, 0, 7000, std::get<1>(aCluster.back()) + std::get<1>(aCluster.front()), std::get<1>(cCluster.back()), "AvC");
|
||||
plotter->Fill2D("A1vsA2", 1600, 0, 20000, 1600, 0, 20000, std::get<1>(aCluster.back()), std::get<1>(aCluster.front()), "AvC");
|
||||
// if (std::get<1>(aCluster.back()) > std::get<1>(aCluster.front()))
|
||||
{
|
||||
plotter->Fill2D("A1vsA2", 400, 0, 20000, 400, 0, 20000, std::get<1>(aCluster.back()), std::get<1>(aCluster.front()), "AvC");
|
||||
plotter->Fill2D("A1vsA2_TC" + std::to_string(PCQQQTimeCut), 400, 0, 20000, 400, 0, 20000, std::get<1>(aCluster.back()), std::get<1>(aCluster.front()), "AvC");
|
||||
}
|
||||
for (auto qqqevent : QQQ_Events)
|
||||
{
|
||||
plotter->Fill2D("qqqER_2Anode_dESum", 800, 0, 10, 800, 0, 14000, qqqevent.Energy1, std::get<1>(cCluster.back()) + std::get<1>(cCluster.front()), "AvC");
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -707,17 +872,90 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
plotter->Fill1D("dt_pcC_qqqW", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2);
|
||||
plotter->Fill2D("dE_E_AnodeQQQR", 400, 0, 10, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1);
|
||||
plotter->Fill2D("dE_E_CathodeQQQR", 400, 0, 10, 800, 0, 10000, qqqevent.Energy2, pcevent.Energy2);
|
||||
double sinTheta = TMath::Sin((qqqevent.pos - TVector3(0, 0, vertexpos)).Theta()) / TMath::Sin((TVector3(51.5, 0, qqqpos) - TVector3(0, 0, vertexpos)).Theta());
|
||||
plotter->Fill2D("dE2_E_AnodeQQQR", 400, 0, 10, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1 * sinTheta);
|
||||
plotter->Fill2D("dE2_E_CathodeQQQR", 400, 0, 10, 800, 0, 10000, qqqevent.Energy2, pcevent.Energy2 * sinTheta);
|
||||
double sinTheta = TMath::Sin((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta());
|
||||
if ((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI > 52)
|
||||
{
|
||||
plotter->Fill2D("dE2_E_AnodeQQQR_outer", 400, 0, 10, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1 * sinTheta);
|
||||
plotter->Fill2D("dE2_E_CathodeQQQR_outer", 400, 0, 10, 800, 0, 10000, qqqevent.Energy2, pcevent.Energy2 * sinTheta);
|
||||
plotter->Fill2D("dE_E_AnodeQQQR_outer", 400, 0, 10, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1);
|
||||
plotter->Fill2D("dE_E_CathodeQQQR_outer", 400, 0, 10, 800, 0, 10000, qqqevent.Energy2, pcevent.Energy2);
|
||||
}
|
||||
else
|
||||
{
|
||||
plotter->Fill2D("dE2_E_AnodeQQQR_inner", 400, 0, 10, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1 * sinTheta);
|
||||
plotter->Fill2D("dE2_E_CathodeQQQR_inner", 400, 0, 10, 800, 0, 10000, qqqevent.Energy2, pcevent.Energy2 * sinTheta);
|
||||
plotter->Fill2D("dE_E_AnodeQQQR_inner", 400, 0, 10, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1);
|
||||
plotter->Fill2D("dE_E_CathodeQQQR_inner", 400, 0, 10, 800, 0, 10000, qqqevent.Energy2, pcevent.Energy2);
|
||||
}
|
||||
|
||||
bool timecut = (qqqevent.Time1 - pcevent.Time1 < -150);
|
||||
if (timecut)
|
||||
{ // && qqqevent.pos.Phi() <= pcevent.pos.Phi()+TMath::Pi()/4. && qqqevent.pos.Phi() >= pcevent.pos.Phi()-TMath::Pi()/4. ) {
|
||||
plotter->Fill2D("dE_theta_AnodeQQQR", 75, 0, 90, 400, 0, 20000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy1);
|
||||
plotter->Fill2D("dE2_theta_AnodeQQQR", 75, 0, 90, 400, 0, 20000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy1 * sinTheta);
|
||||
|
||||
plotter->Fill2D("E_theta_AnodeQQQR", 75, 0, 90, 300, 0, 15, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, qqqevent.Energy1);
|
||||
plotter->Fill2D("E2_theta_AnodeQQQR", 75, 0, 90, 300, 0, 15, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, qqqevent.Energy1);
|
||||
plotter->Fill2D("Etot2_theta_AnodeQQQR", 75, 0, 90, 300, 0, 15, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, qqqevent.Energy1 + pcevent.Energy1 * anode_gain * sinTheta);
|
||||
|
||||
plotter->Fill2D("dE_theta_CathodeQQQR", 75, 0, 90, 800, 0, 10000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy2);
|
||||
plotter->Fill2D("dE2_theta_CathodeQQQR", 75, 0, 90, 800, 0, 10000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, pcevent.Energy2 * sinTheta);
|
||||
|
||||
plotter->Fill2D("dE_phi_AnodeQQQR", 100, -180, 180, 800, 0, 40000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Phi() * 180 / M_PI, pcevent.Energy1);
|
||||
|
||||
plotter->Fill2D("dE_phi_CathodeQQQR", 100, -180, 180, 800, 0, 10000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Phi() * 180 / M_PI, pcevent.Energy2);
|
||||
plotter->Fill1D("PCZ", 800, -200, 200, pcevent.pos.Z(), "phicut");
|
||||
plotter->Fill1D("PCZ_phicut_a" + std::to_string(aClusters.at(0).size()) + "_c" + std::to_string(cClusters.at(0).size()), 800, -200, 200, pcevent.pos.Z(), "wiremult");
|
||||
|
||||
double pcz_guess_37 = 37. / TMath::Tan((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()) + source_vertex;
|
||||
plotter->Fill2D("pczguess_vs_pc_37", 180, 0, 200, 150, 0, 200, pcz_guess_37, pcevent.pos.Z(), "phicut");
|
||||
|
||||
double pcz_guess_42 = 42. / TMath::Tan((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()) + source_vertex;
|
||||
plotter->Fill2D("pczguess_vs_pc_42", 180, 0, 200, 150, 0, 200, pcz_guess_42, pcevent.pos.Z(), "phicut");
|
||||
|
||||
double pcz_guess_int = z_to_crossover_rho(pcevent.pos.Z()) / TMath::Tan((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()) + source_vertex;
|
||||
// plotter->Fill2D("pczguess_vs_pc_int",180,0,200,150,0,200,pcz_guess_int,pcevent.pos.Z(),"phicut");
|
||||
plotter->Fill2D("pczguess_vs_pc_int", 180, 0, 200, 600, -400, 400, pcz_guess_int, pcevent.pos.Z(), "phicut");
|
||||
|
||||
double qqqrho = qqqevent.pos.Perp();
|
||||
double qqqz = (qqqevent.pos - TVector3(0, 0, source_vertex)).Z();
|
||||
double tan_theta = qqqrho / qqqz;
|
||||
double pcz_guess_int2 = z_to_crossover_rho(pcevent.pos.Z()) / tan_theta + source_vertex;
|
||||
plotter->Fill2D("pczguess_vs_pc_int2", 180, 0, 200, 150, 0, 200, pcz_guess_int, pcevent.pos.Z(), "phicut");
|
||||
plotter->Fill2D("pczguess_vs_pc_int2_a" + std::to_string(pcevent.multi1) + "_c" + std::to_string(pcevent.multi2), 180, 0, 200, 150, 0, 200, pcz_guess_int, pcevent.pos.Z(), "phicut");
|
||||
|
||||
double pcz_guess = pcz_guess_int;
|
||||
plotter->Fill2D("pctheta_vs_qqqtheta", 320, 0, 160, 320, 0, 160, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, (pcevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, "phicut");
|
||||
|
||||
plotter->Fill2D("pczguess_vs_pc_phi=" + std::to_string(qqqevent.pos.Phi() * 180. / M_PI), 300, 0, 200, 150, 0, 200, pcz_guess, pcevent.pos.Z(), "phicut");
|
||||
|
||||
// plotter->Fill1D("PCZ",800,-200,200,pcevent.pos.Z(),"phicut");
|
||||
}
|
||||
|
||||
if (qqqevent.pos.Phi() <= pcevent.pos.Phi() + TMath::Pi() / 4. && qqqevent.pos.Phi() >= pcevent.pos.Phi() - TMath::Pi() / 4.)
|
||||
{
|
||||
plotter->Fill1D("PCZ", 800, -200, 200, pcevent.pos.Z(), "phicut");
|
||||
double pcz_guess = pcrad / TMath::Tan((qqqevent.pos - TVector3(0, 0, vertexpos)).Theta()) + vertexpos; // this is ideally kept to be all QQQ+userinput for calibration of pcz
|
||||
plotter->Fill2D("pczguess_vs_pc", 300, 0, 200, 150, 0, 200, pcz_guess, pcevent.pos.Z(), "phicut");
|
||||
double pcz_guess_37 = 37. / TMath::Tan((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()) + source_vertex;
|
||||
plotter->Fill2D("pczguess_vs_pc_37", 180, 0, 200, 150, 0, 200, pcz_guess_37, pcevent.pos.Z(), "phicut");
|
||||
|
||||
double pcz_guess_42 = 42. / TMath::Tan((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()) + source_vertex;
|
||||
plotter->Fill2D("pczguess_vs_pc_42", 180, 0, 200, 150, 0, 200, pcz_guess_42, pcevent.pos.Z(), "phicut");
|
||||
|
||||
double pcz_guess_int = z_to_crossover_rho(pcevent.pos.Z()) / TMath::Tan((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()) + source_vertex;
|
||||
// plotter->Fill2D("pczguess_vs_pc_int",180,0,200,150,0,200,pcz_guess_int,pcevent.pos.Z(),"phicut");
|
||||
plotter->Fill2D("pczguess_vs_pc_int", 180, 0, 200, 600, -400, 400, pcz_guess_int, pcevent.pos.Z(), "phicut");
|
||||
|
||||
double qqqrho = qqqevent.pos.Perp();
|
||||
double qqqz = (qqqevent.pos - TVector3(0, 0, source_vertex)).Z();
|
||||
double tan_theta = qqqrho / qqqz;
|
||||
double pcz_guess_int2 = z_to_crossover_rho(pcevent.pos.Z()) / tan_theta + source_vertex;
|
||||
plotter->Fill2D("pczguess_vs_pc_int2", 180, 0, 200, 150, 0, 200, pcz_guess_int, pcevent.pos.Z(), "phicut");
|
||||
plotter->Fill2D("pczguess_vs_pc_int2_a" + std::to_string(pcevent.multi1) + "_c" + std::to_string(pcevent.multi2), 180, 0, 200, 150, 0, 200, pcz_guess_int, pcevent.pos.Z(), "phicut");
|
||||
|
||||
double pcz_guess = pcz_guess_int;
|
||||
plotter->Fill2D("pctheta_vs_qqqtheta", 320, 0, 160, 320, 0, 160, (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, (pcevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, "phicut");
|
||||
|
||||
plotter->Fill2D("pczguess_vs_pc_phi=" + std::to_string(qqqevent.pos.Phi() * 180. / M_PI), 300, 0, 200, 150, 0, 200, pcz_guess, pcevent.pos.Z(), "phicut");
|
||||
// plotter->Fill1D("PCZ",800,-200,200,pcevent.pos.Z(),"phicut");
|
||||
}
|
||||
|
||||
if (cSize == 1)
|
||||
|
|
@ -786,6 +1024,11 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
{
|
||||
corrcatMax.push_back(std::pair<int, double>(cID, cE));
|
||||
cESum += cE;
|
||||
if (cE > cEMax)
|
||||
{
|
||||
cEMax = cE;
|
||||
cIDMax = cID;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -800,13 +1043,13 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
double x = 0, y = 0, z = 0;
|
||||
for (const auto &corr : corrcatMax)
|
||||
{
|
||||
if (Crossover[aIDMax][corr.first][0].z > 9000000)
|
||||
if (pwinstance.Crossover[aIDMax][corr.first][0].z > 9000000)
|
||||
continue;
|
||||
if (cESum > 0)
|
||||
{
|
||||
x += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].x;
|
||||
y += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].y;
|
||||
z += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].z;
|
||||
x += (corr.second) / cESum * pwinstance.Crossover[aIDMax][corr.first][0].x;
|
||||
y += (corr.second) / cESum * pwinstance.Crossover[aIDMax][corr.first][0].y;
|
||||
z += (corr.second) / cESum * pwinstance.Crossover[aIDMax][corr.first][0].z;
|
||||
}
|
||||
}
|
||||
if (x == 0 && y == 0 && z == 0)
|
||||
|
|
@ -823,13 +1066,35 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
PCQQQPhiCut = true;
|
||||
}
|
||||
|
||||
for (double Tz = 60; Tz <= 100; Tz += 1.0)
|
||||
if (anodeIndex != -1 && cathodeIndex != -1 && hitPos.Perp() != 0 && anodeIntersection.Perp() != 0 && realtime)
|
||||
{
|
||||
TVector3 TargetPos(0, 0, Tz);
|
||||
if (PCQQQPhiCut && anodeIntersection.Perp() > 0 && anodeIntersection.Z() != 0 && cathodeHits.size() >= 2)
|
||||
can1->Modified();
|
||||
can1->Update();
|
||||
for (auto cath : corrcatMax)
|
||||
{
|
||||
plotter->Fill2D("Inttheta_vs_QQQtheta_TC" + std::to_string(PCQQQTimeCut) + "_TZ" + std::to_string(Tz), 400, 0, 180, 90, 0, 90, (anodeIntersection - TargetPos).Theta() * 180. / TMath::Pi(), (hitPos - TargetPos).Theta() * 180. / TMath::Pi(), "TPosVariation");
|
||||
// plotter->Fill2D("R_ratio_to_Z_ratio" + std::to_string(PCQQQTimeCut) + "_TZ" + std::to_string(Tz), 100, -2, 2, 100, -2, 2, (anodeIntersection - TargetPos).Z()/(hitPos-TargetPos).Z(), ((anodeIntersection - TargetPos).Perp()+2.5)/(hitPos-TargetPos).Perp(), "TPosVariation");
|
||||
plc[cath.first]->SetLineWidth(3);
|
||||
// plc[cath.first]->SetLineStyle(kLine);
|
||||
}
|
||||
for (auto anodeW : anodeHits)
|
||||
{
|
||||
pla[anodeW.first]->SetLineWidth(3);
|
||||
// pla[anodeW.first]->SetLineStyle(kLine);
|
||||
}
|
||||
// can2->Modified();
|
||||
can2->Update();
|
||||
while (can1->WaitPrimitive())
|
||||
;
|
||||
// pla[anodeIndex]->SetLineWidth(1);
|
||||
// pla[anodeIndex]->SetLineStyle(kDotted);
|
||||
for (auto anodeW : anodeHits)
|
||||
{
|
||||
pla[anodeW.first]->SetLineWidth(1);
|
||||
pla[anodeW.first]->SetLineStyle(kDotted);
|
||||
}
|
||||
for (auto cath : corrcatMax)
|
||||
{
|
||||
plc[cathodeIndex]->SetLineStyle(kDotted);
|
||||
plc[cath.first]->SetLineWidth(1);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -918,12 +1183,12 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
|
||||
if (HitNonZero && anodeIntersection.Z() != 0)
|
||||
{
|
||||
pw_contr.CalTrack2(hitPos, anodeIntersection);
|
||||
plotter->Fill1D("VertexRecon", 600, -1300, 1300, pw_contr.GetZ0());
|
||||
plotter->Fill1D("VertexRecon_TC" + std::to_string(PCQQQTimeCut) + "_PhiC" + std::to_string(PCQQQPhiCut), 600, -1300, 1300, pw_contr.GetZ0());
|
||||
pwinstance.CalTrack2(hitPos, anodeIntersection);
|
||||
plotter->Fill1D("VertexRecon", 600, -1300, 1300, pwinstance.GetZ0());
|
||||
plotter->Fill1D("VertexRecon_TC" + std::to_string(PCQQQTimeCut) + "_PhiC" + std::to_string(PCQQQPhiCut), 600, -1300, 1300, pwinstance.GetZ0());
|
||||
|
||||
if (cathodeHits.size() == 2)
|
||||
plotter->Fill1D("VertexRecon_2c_TC" + std::to_string(PCQQQTimeCut) + "_PhiC" + std::to_string(PCQQQPhiCut), 600, -1300, 1300, pw_contr.GetZ0());
|
||||
plotter->Fill1D("VertexRecon_2c_TC" + std::to_string(PCQQQTimeCut) + "_PhiC" + std::to_string(PCQQQPhiCut), 600, -1300, 1300, pwinstance.GetZ0());
|
||||
|
||||
TVector3 x2(anodeIntersection), x1(hitPos);
|
||||
|
||||
|
|
@ -1005,6 +1270,7 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
// if (anodeIntersection.Z() != 0)
|
||||
{
|
||||
plotter->Fill2D("PC_Z_vs_QQQRing", 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCzQQQ");
|
||||
plotter->Fill2D("PC_Z_vs_QQQRho", 600, -300, 300, 40, 40, 110, anodeIntersection.Z(), hitPos.Perp(), "hPCzQQQ");
|
||||
}
|
||||
|
||||
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 2)
|
||||
|
|
@ -1046,15 +1312,15 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
} // j loop end
|
||||
} // qqq i loop end
|
||||
|
||||
TVector3 guessVertex(0, 0, 90.); // for run12, subtract anodeIntersection.Z() by ~74.0 seems to work
|
||||
TVector3 guessVertex(0, 0, source_vertex); // for run12, subtract anodeIntersection.Z() by ~74.0 seems to work
|
||||
// rho=40.0 mm is halfway between the cathodes(rho=42) and anodes(rho=37)
|
||||
double pcz_guess = 42.0 / TMath::Tan((hitPos - guessVertex).Theta()) + guessVertex.Z(); // this is ideally kept to be all QQQ+userinput for calibration of pcz
|
||||
double pcz_guess = z_to_crossover_rho(anodeIntersection.Z()) / TMath::Tan((hitPos - guessVertex).Theta()) + guessVertex.Z(); // this is ideally kept to be all QQQ+userinput for calibration of pcz
|
||||
if (PCQQQTimeCut && PCQQQPhiCut && hitPos.Perp() > 0 && anodeIntersection.Perp() > 0 && cathodeHits.size() >= 2)
|
||||
{
|
||||
plotter->Fill2D("pczguess_vs_qqqE", 100, 0, 200, 800, 0, 20, pcz_guess, qqqenergy, "pczguess");
|
||||
double pczoffset = 30.0;
|
||||
double pczoffset = 0.0;
|
||||
// plotter->Fill2D("pczguess_vs_pcz_rad="+std::to_string(hitPos.Perp()),100,0,200,150,0,200,pcz_guess,anodeIntersection.Z(),"pczguess"); //entirely qqq-derived position vs entirely PC derived position
|
||||
plotter->Fill2D("pczguess_vs_pcz_phi=" + std::to_string(hitPos.Phi() * 180. / M_PI), 100, 0, 200, 150, 0, 200, pcz_guess, anodeIntersection.Z() + pczoffset, "pczguess"); // entirely qqq-derived position vs entirely PC derived position
|
||||
plotter->Fill2D("pczguess_vs_pcz_phi=" + std::to_string(hitPos.Phi() * 180. / M_PI), 100, 0, 200, 150, 0, 200, pcz_guess, anodeIntersection.Z(), "pczguess"); // entirely qqq-derived position vs entirely PC derived position
|
||||
plotter->Fill2D("pczguess_vs_pcz", 100, 0, 200, 150, 0, 200, pcz_guess, anodeIntersection.Z() + pczoffset);
|
||||
plotter->Fill2D("pcz_vs_pcPhi_rad=" + std::to_string(hitPos.Perp()), 360, 0, 360, 150, 0, 200, anodeIntersection.Phi() * 180. / M_PI, anodeIntersection.Z() + pczoffset, "pczguess");
|
||||
}
|
||||
|
|
@ -1065,16 +1331,25 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
plotter->Fill2D("PCPhi_vs_SX3Strip", 100, -200, 200, 8 * 24, 0, 8 * 24, anodeIntersection.Phi() * 180. / TMath::Pi(), sx3.id[i] * 8 + sx3.ch[i]);
|
||||
}
|
||||
|
||||
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 3)
|
||||
if (anodeIntersection.Perp() != 0)
|
||||
{
|
||||
plotter->Fill1D("PC_Z_proj_3C", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
|
||||
}
|
||||
plotter->Fill2D("AnodeMaxE_Vs_Cathode_Sum_Energy", 2000, 0, 20000, 2000, 0, 10000, aEMax, cESum, "hGMPC");
|
||||
plotter->Fill2D("AnodeSumE_Vs_Cathode_Max_Energy", 800, 0, 20000, 800, 0, 10000, aESum, cEMax, "hGMPC");
|
||||
plotter->Fill2D("AnodeMaxE_Vs_Cathode_Max_Energy", 800, 0, 20000, 800, 0, 10000, aEMax, cEMax, "hGMPC");
|
||||
// double sinTheta = TMath::Sin((anodeIntersection - TVector3(0,0,source_vertex)).Theta());///TMath::Sin((TVector3(51.5,0,128.) - TVector3(0,0,85)).Theta());
|
||||
// plotter->Fill2D("AnodeMaxE_Vs_Cathode_Max_Energy_path_corrected", 800, 0, 20000, 800, 0, 10000, aEMax*sinTheta, cEMax*sinTheta, "hGMPC");
|
||||
plotter->Fill2D("AnodeSumE_Vs_Cathode_Sum_Energy", 800, 0, 20000, 800, 0, 10000, aESum, cESum, "hGMPC");
|
||||
plotter->Fill2D("AnodeSumE_Vs_Cathode_Max_Energy_TC" + std::to_string(PCQQQTimeCut) + "_PC" + std::to_string(PCQQQPhiCut), 800, 0, 20000, 800, 0, 10000, aESum, cEMax, "hGMPC");
|
||||
// plotter->Fill2D("AnodeSumE_Vs_Cathode_Max_Energy_path_corrected"+std::to_string(PCQQQTimeCut)+"_PC"+std::to_string(PCQQQPhiCut), 800, 0, 20000, 800, 0, 10000, aESum*sinTheta, cEMax*sinTheta, "hGMPC");
|
||||
// plotter->Fill2D("AnodeSumE_Vs_Cathode_Max_Energy_path_corrected", 800, 0, 20000, 800, 0, 10000, aESum*sinTheta, cEMax*sinTheta, "hGMPC");
|
||||
|
||||
plotter->Fill2D("AnodeMaxE_Vs_Cathode_Sum_Energy", 2000, 0, 30000, 2000, 0, 30000, aEMax, cESum, "hGMPC");
|
||||
plotter->Fill1D("Correlated_Cathode_MaxAnode", 6, 0, 5, corrcatMax.size(), "hGMPC");
|
||||
plotter->Fill2D("Correlated_Cathode_VS_MaxAnodeEnergy", 6, 0, 5, 2000, 0, 30000, corrcatMax.size(), aEMax, "hGMPC");
|
||||
plotter->Fill1D("AnodeHits", 12, 0, 11, anodeHits.size(), "hGMPC");
|
||||
plotter->Fill2D("AnodeMaxE_vs_AnodeHits", 12, 0, 11, 2000, 0, 30000, anodeHits.size(), aEMax, "hGMPC");
|
||||
if (PCQQQTimeCut && PCQQQPhiCut)
|
||||
{
|
||||
plotter->Fill2D("AnodeSumE_Vs_Cathode_Max_Energy_TC" + std::to_string(PCQQQTimeCut) + "_PC" + std::to_string(PCQQQPhiCut) + "_cMax" + std::to_string(cIDMax), 800, 0, 20000, 800, 0, 10000, aESum, cEMax, "hGMPC");
|
||||
}
|
||||
// plotter->Fill2D("AnodeSumE_Vs_CathodeSum_Energy_path_corrected", 800, 0, 20000, 800, 0, 10000, aESum*sinTheta, cESum*sinTheta, "hGMPC");
|
||||
// plotter->Fill2D("AnodeSumE_Vs_CathodeSum_Energy_path_corrected_TC"+std::to_string(PCQQQTimeCut)+"_PC"+std::to_string(PCQQQPhiCut), 800, 0, 20000, 800, 0, 10000, aESum*sinTheta, cESum*sinTheta, "hGMPC"); */
|
||||
}
|
||||
|
||||
if (anodeHits.size() < 1)
|
||||
{
|
||||
|
|
|
|||
27
anasen_fem/README.md
Executable file
27
anasen_fem/README.md
Executable file
|
|
@ -0,0 +1,27 @@
|
|||
### README for ANASEN fem simulations:
|
||||
|
||||
* There are a few iterations of these simulations that already exist. Be sure to also locate and refer to them if necessary.
|
||||
* Install gmsh and its python api by running (Ubuntu 22.04 LTS)
|
||||
|
||||
```
|
||||
sudo apt install gmsh python3-gmsh
|
||||
```
|
||||
|
||||
* Gmsh gives us the tools to create a meshgrid that samples the 2d space appropriately to plot the field/equipotential lines.
|
||||
* The output file typically has the .msh extension. This is read as input to Elmer, which is the FEM differential-equation solver.
|
||||
|
||||
* Install Elmer via the following steps:
|
||||
|
||||
```
|
||||
sudo add-apt-repository ppa:elmer-csc-ubuntu/elmer-csc-ppa
|
||||
sudo apt install elmerfem-csc-eg
|
||||
```
|
||||
|
||||
* Install ParaView for visualizations by downloading from the Linux .tar.gz link at https://www.paraview.org/download/
|
||||
- The current version is tested to work on Paraview 6.1.0. The default version in Ubuntu 22.04 repositories has some trouble with scripting
|
||||
* v0.0.1, March 10 2026
|
||||
- 2d simulations of fields only. gmsh for meshing, elmer for fem, paraview to plot
|
||||
- Before running, open `paraview_plotter.py` to make the bash shebang (#!) point to the location of `pvpython` or `pvbatch`
|
||||
- `python3 run.py` should run everything in order, and is hopefully all the files are self-documenting
|
||||
* v0.0.2, planned TODO
|
||||
- Garfield to take Elmer results and perform charge-transport
|
||||
2
anasen_fem/clean.sh
Executable file
2
anasen_fem/clean.sh
Executable file
|
|
@ -0,0 +1,2 @@
|
|||
rm wires2d.msh
|
||||
rm wires2d/*
|
||||
114
anasen_fem/junk/wires.py
Executable file
114
anasen_fem/junk/wires.py
Executable file
|
|
@ -0,0 +1,114 @@
|
|||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from mpl_toolkits.mplot3d import Axes3D
|
||||
|
||||
acolors = plt.get_cmap('tab20',24)
|
||||
ccolors = plt.get_cmap('tab20',24)
|
||||
|
||||
k=-2*np.pi/24.
|
||||
offset = 6*k + 3*k #-pi/2
|
||||
xarra_1 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)])
|
||||
yarra_1 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)])
|
||||
labelsa_1 = np.array([i for i in np.arange(0,24)])
|
||||
|
||||
fig,ax = plt.subplots(figsize=(10,10))
|
||||
ax.invert_yaxis()
|
||||
ax.plot(xarra_1,yarra_1,"x",label="anode, z=-L/2")
|
||||
|
||||
for x,y,label in zip(xarra_1,yarra_1,labelsa_1):
|
||||
ax.text(x,y,label)
|
||||
|
||||
kc=2*np.pi/24.
|
||||
offsetc = -4*kc + 2*kc - np.pi/24 #-pi/4
|
||||
xarrc_1 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,24)])
|
||||
yarrc_1 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,24)])
|
||||
labelsc_1 = np.array([i for i in np.arange(0,24)])
|
||||
|
||||
ax.plot(xarrc_1,yarrc_1,"o",label="cathode, z=-L/2, where they are picked up")
|
||||
|
||||
for x,y,label in zip(xarrc_1,yarrc_1,labelsc_1):
|
||||
ax.text(x,y,label)
|
||||
plt.title("z=-L/2 plane, beam going into the plane along +z, (+x right, +y down)")
|
||||
plt.grid()
|
||||
plt.legend()
|
||||
plt.savefig("plane1.png")
|
||||
plt.show()
|
||||
|
||||
fig,ax = plt.subplots(figsize=(10,10))
|
||||
ax.invert_yaxis()
|
||||
|
||||
offset = offset-3*k
|
||||
xarra_2 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)])
|
||||
yarra_2 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)])
|
||||
labelsa_2 = np.array([i for i in np.arange(0,24)])
|
||||
|
||||
ax.plot(xarra_2,yarra_2,"x",label="anode, z=+L/2, where they are picked up")
|
||||
|
||||
for x,y,label in zip(xarra_2,yarra_2,labelsa_2):
|
||||
ax.text(x,y,label)
|
||||
|
||||
offsetc = offsetc-3*kc
|
||||
xarrc_2 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,24)])
|
||||
yarrc_2 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,24)])
|
||||
labelsc_2 = np.array([i for i in np.arange(0,24)])
|
||||
|
||||
ax.plot(xarrc_2,yarrc_2,"o",label="cathode, z=+L/2")
|
||||
|
||||
for x,y,label in zip(xarrc_2,yarrc_2,labelsc_2):
|
||||
ax.text(x,y,label)
|
||||
plt.title("z=+L/2 plane, beam going into the plane along +z, (+x right, +y down)")
|
||||
plt.grid()
|
||||
plt.legend()
|
||||
plt.savefig("plane2.png")
|
||||
plt.show()
|
||||
|
||||
fig = plt.figure(figsize=(10,10))
|
||||
ax = fig.add_subplot(111,projection='3d')
|
||||
ax.set_xlabel("x")
|
||||
ax.set_ylabel("y")
|
||||
ax.set_zlabel("z")
|
||||
|
||||
xx_a = np.array([[x1,x2] for x1,x2 in zip(xarra_1,xarra_2)])
|
||||
yy_a = np.array([[y1,y2] for y1,y2 in zip(yarra_1,yarra_2)])
|
||||
zz_a = np.array([[-173.5,173.5] for x1,x2 in zip(xarra_1,xarra_2)])
|
||||
for i,[xx,yy,zz] in enumerate(zip(xx_a,yy_a,zz_a)):
|
||||
ax.plot(xx,yy,zz,'-',color=acolors(i/24))
|
||||
|
||||
for i,[x,y,label] in enumerate(zip(xarra_1,yarra_1,labelsa_1)):
|
||||
ax.text(x,y,-173,"a"+str(label),color=acolors(i/24))
|
||||
for i,[x,y,label] in enumerate(zip(xarra_2,yarra_2,labelsa_2)):
|
||||
ax.text(x,y,+173,"a"+str(label),color=acolors(i/24))
|
||||
|
||||
xx_c = np.array([[x1,x2] for x1,x2 in zip(xarrc_1,xarrc_2)])
|
||||
yy_c = np.array([[y1,y2] for y1,y2 in zip(yarrc_1,yarrc_2)])
|
||||
zz_c = np.array([[-173.5,173.5] for x1,x2 in zip(xarrc_1,xarrc_2)])
|
||||
for i,[xx,yy,zz] in enumerate(zip(xx_c,yy_c,zz_c)):
|
||||
ax.plot(xx,yy,zz,'--',color=ccolors(((47-i)%24)/24))
|
||||
|
||||
for i,[x,y,label] in enumerate(zip(xarrc_1,yarrc_1,labelsc_1)):
|
||||
ax.text(x,y,-173,"c"+str(label),color=ccolors(((25-i)%24)/24))
|
||||
for i,[x,y,label] in enumerate(zip(xarrc_2,yarrc_2,labelsc_2)):
|
||||
ax.text(x,y,+173,"c"+str(label),color=ccolors(((47-i)%24)/24))
|
||||
ax.view_init(elev=-53, azim=-106, roll=18)
|
||||
plt.tight_layout()
|
||||
plt.show()
|
||||
phi_qqq = np.array([[2*np.pi*(-i*16+j+0.5)/(16*4) for i in range(4)] for j in range(16)])
|
||||
print(phi_qqq)
|
||||
#'''
|
||||
for i,[phi1,phi2,phi3,phi4] in enumerate(phi_qqq):
|
||||
ax.plot([50*np.cos(phi1),100*np.cos(phi1)],[50*np.sin(phi1),100*np.sin(phi1)],[100,100],'-',color='red')
|
||||
ax.text(104*np.cos(phi1),104*np.sin(phi1),100,"0_%d"%(i),color="red")
|
||||
|
||||
ax.plot([50*np.cos(phi2),100*np.cos(phi2)],[50*np.sin(phi2),100*np.sin(phi2)],[100,100],'-',color='green')
|
||||
ax.text(104*np.cos(phi2),104*np.sin(phi2),100,"1_%d"%(i),color="green")
|
||||
|
||||
ax.plot([50*np.cos(phi3),100*np.cos(phi3)],[50*np.sin(phi3),100*np.sin(phi3)],[100,100],'-',color='blue')
|
||||
ax.text(104*np.cos(phi3),104*np.sin(phi3),100,"2_%d"%(i),color="blue")
|
||||
|
||||
ax.plot([50*np.cos(phi4),100*np.cos(phi4)],[50*np.sin(phi4),100*np.sin(phi4)],[100,100],'-',color='brown')
|
||||
ax.text(104*np.cos(phi4),104*np.sin(phi4),100,"3_%d"%(i),color="brown")
|
||||
#'''
|
||||
#coords_qqq = np.array([[50*np.cos(phi),100*np.sin(phi),100] for phi in phi_qqq]).T
|
||||
plt.tight_layout()
|
||||
plt.show()
|
||||
|
||||
71
anasen_fem/junk/wires2d_test.sif
Executable file
71
anasen_fem/junk/wires2d_test.sif
Executable file
|
|
@ -0,0 +1,71 @@
|
|||
Check Keywords Warn
|
||||
|
||||
Header
|
||||
Mesh DB "." "wires2d"
|
||||
End
|
||||
|
||||
Simulation
|
||||
Coordinate System = Cartesian 2D
|
||||
Simulation Type = Steady State
|
||||
Steady State Max Iterations = 1
|
||||
Output File = "elstatics.result"
|
||||
Post File = "elstatics.ep"
|
||||
End
|
||||
|
||||
Constants
|
||||
Permittivity Of Vacuum = 8.8542e-12
|
||||
End
|
||||
|
||||
|
||||
Body 1
|
||||
Target Bodies(1) = 13
|
||||
Equation = 1
|
||||
Material = 1
|
||||
End
|
||||
|
||||
|
||||
Equation 1
|
||||
Active Solvers(2) = 1 2
|
||||
End
|
||||
|
||||
|
||||
Material 1
|
||||
Relative Permittivity = 1
|
||||
End
|
||||
|
||||
|
||||
Solver 1
|
||||
Equation = Electrostatics
|
||||
Procedure = "StatElecSolve" "StatElecSolver"
|
||||
|
||||
Variable = Potential
|
||||
Variable DOFs = 1
|
||||
|
||||
Calculate Electric Field = True
|
||||
Calculate Electric Flux = False
|
||||
|
||||
Linear System Solver = Iterative
|
||||
Linear System Iterative Method = CG
|
||||
Linear System Max Iterations = 500
|
||||
Linear System Convergence Tolerance = 1.0e-8
|
||||
Linear System Preconditioning = ILU1
|
||||
End
|
||||
|
||||
|
||||
Solver 2
|
||||
Equation = Electric Force
|
||||
Procedure = "ElectricForce" "StatElecForce"
|
||||
End
|
||||
|
||||
|
||||
Boundary Condition 1
|
||||
Target Boundaries = 10
|
||||
Potential = 650
|
||||
Calculate Electric Force = True
|
||||
End
|
||||
|
||||
|
||||
Boundary Condition 2
|
||||
Target Boundaries = 20
|
||||
Potential = 0
|
||||
End
|
||||
91
anasen_fem/junk/wires_gmsh.py
Executable file
91
anasen_fem/junk/wires_gmsh.py
Executable file
|
|
@ -0,0 +1,91 @@
|
|||
import numpy as np
|
||||
import gmsh
|
||||
|
||||
gmsh.initialize()
|
||||
gmsh.model.add("adaptive_mesh")
|
||||
gmsh.option.setNumber('General.NumThreads', 4)
|
||||
#gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfElements", 200000)
|
||||
#gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfNodes", 200000)
|
||||
#gmsh.option.setNumber("Mesh.Adapt.MaxIter",5)
|
||||
#gmsh.option.setNumber("Mesh.MeshSizeMin", 5e-3)
|
||||
#gmsh.option.setNumber("Mesh.MeshSizeMax", 10.0)
|
||||
#gmsh.option.setNumber("Mesh.CharacteristicLengthFromCurvature", 0)
|
||||
|
||||
lc = 0.01
|
||||
|
||||
#anodes, plane 1 at -zmax/2
|
||||
k=-2*np.pi/24.
|
||||
offset = 6*k + 3*k #-pi/2
|
||||
xarra_1 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)])
|
||||
yarra_1 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)])
|
||||
|
||||
#cathodes, plane 1 at -zmax/2
|
||||
kc=2*np.pi/24.
|
||||
offsetc = -4*kc + 2*kc - np.pi/24 #-pi/4
|
||||
xarrc_1 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,24)])
|
||||
yarrc_1 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,24)])
|
||||
|
||||
#anodes, plane 2 at +zmax/2
|
||||
offset = offset-3*k
|
||||
xarra_2 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)])
|
||||
yarra_2 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)])
|
||||
|
||||
#cathodes, plane2 at +zmax/2
|
||||
offsetc = offsetc-3*kc
|
||||
xarrc_2 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,24)])
|
||||
yarrc_2 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,24)])
|
||||
|
||||
pa1 = []
|
||||
pa2 = []
|
||||
pc1 = []
|
||||
pc2 = []
|
||||
|
||||
wire_radius = 0.254 #mm
|
||||
anode_wires = []
|
||||
cathode_wires = []
|
||||
aw_tags = [(3,i) for i in range(24)]
|
||||
cw_tags = [(3,i+24) for i in range(24)]
|
||||
|
||||
for i,[xa,ya,xc,yc,xa2,ya2,xc2,yc2] in enumerate(zip(xarra_1,yarra_1,xarrc_1,yarrc_1,xarra_2,yarra_2,xarrc_2,yarrc_2)):
|
||||
print(i,xa,ya,-178.3,xc,yc,-178.3,xa2,ya2,178.3,xc2,yc2,178.3)
|
||||
anode_wires.append(gmsh.model.occ.addCylinder(xa,ya,-178.3,(xa2-xa),(ya2-ya),178.3*2,wire_radius,i)) #x,y,z of first face center, dx,dy,dz of the axis, then the wire radius
|
||||
cathode_wires.append(gmsh.model.occ.addCylinder(xc,yc,-178.3,(xc2-xc),(yc2-yc),178.3*2,wire_radius,i+24)) #cathode tags 24-47, anode 0-23
|
||||
|
||||
|
||||
anasen_barrel = gmsh.model.occ.addCylinder(0,0,-500,0,0,500+605,300,1234) #tag 1234
|
||||
#anasen_barrel = gmsh.model.occ.addCylinder(0,0,-500,0,0,500+605,300,1234) #tag 1234
|
||||
|
||||
gmsh.model.occ.synchronize()
|
||||
|
||||
all_wires = aw_tags+cw_tags
|
||||
gmsh.model.occ.fragment([(3,1234)],all_wires)
|
||||
gmsh.model.occ.removeAllDuplicates()
|
||||
gmsh.model.occ.synchronize()
|
||||
gmsh.option.setNumber("Geometry.Tolerance", 1e-6)
|
||||
gmsh.option.setNumber("Geometry.OCCFixDegenerated", 1)
|
||||
gmsh.option.setNumber("Geometry.OCCFixSmallEdges", 1)
|
||||
gmsh.option.setNumber("Geometry.OCCFixSmallFaces", 1)
|
||||
|
||||
wire_surfs = []
|
||||
for w in anode_wires + cathode_wires:
|
||||
wire_surfs += [s[1] for s in gmsh.model.getBoundary([(3,w)], oriented=False) if s[0] == 2]
|
||||
#'''
|
||||
f1 = gmsh.model.mesh.field.add("Distance")
|
||||
gmsh.model.mesh.field.setNumbers(f1, "FacesList", wire_surfs) # Example curves
|
||||
f2 = gmsh.model.mesh.field.add("Threshold")
|
||||
gmsh.model.mesh.field.setNumber(f2, "InField", f1)
|
||||
gmsh.model.mesh.field.setNumber(f2, "SizeMin", 0.05)
|
||||
gmsh.model.mesh.field.setNumber(f2, "SizeMax", 5.)
|
||||
gmsh.model.mesh.field.setNumber(f2, "DistMin", 1.)
|
||||
gmsh.model.mesh.field.setNumber(f2, "DistMax", 20.)
|
||||
gmsh.model.mesh.field.setAsBackgroundMesh(f2)
|
||||
#'''
|
||||
gmsh.option.setNumber("Mesh.Algorithm", 2)
|
||||
gmsh.option.setNumber("Mesh.Algorithm3D", 1) # For 3D meshes
|
||||
|
||||
gmsh.model.mesh.generate(dim=3)
|
||||
#gmsh.model.mesh.refine()
|
||||
#gmsh.model.mesh.refine()
|
||||
#gmsh.model.mesh.refine()
|
||||
gmsh.fltk.run()
|
||||
gmsh.finalize()
|
||||
87
anasen_fem/junk/wires_gmsh_bc.py
Executable file
87
anasen_fem/junk/wires_gmsh_bc.py
Executable file
|
|
@ -0,0 +1,87 @@
|
|||
import numpy as np
|
||||
import gmsh
|
||||
|
||||
gmsh.initialize()
|
||||
#gmsh.model.add("adaptive_mesh")
|
||||
gmsh.option.setNumber('General.NumThreads', 4)
|
||||
#gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfElements", 200000)
|
||||
#gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfNodes", 200000)
|
||||
#gmsh.option.setNumber("Mesh.Adapt.MaxIter",5)
|
||||
#gmsh.option.setNumber("Mesh.MeshSizeMin", 5e-3)
|
||||
#gmsh.option.setNumber("Mesh.MeshSizeMax", 10.0)
|
||||
gmsh.option.setNumber("Geometry.Tolerance", 1e-2)
|
||||
#gmsh.option.setNumber("Mesh.CharacteristicLengthFromCurvature", 0)
|
||||
|
||||
lc = 0.04
|
||||
#anodes, plane 1 at -zmax/2
|
||||
k=-2*np.pi/24.
|
||||
offset = 6*k + 3*k #-pi/2
|
||||
xarra_1 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)])
|
||||
yarra_1 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)])
|
||||
|
||||
#cathodes, plane 1 at -zmax/2
|
||||
kc=2*np.pi/24.
|
||||
offsetc = -4*kc + 2*kc - np.pi/24 #-pi/4
|
||||
xarrc_1 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,24)])
|
||||
yarrc_1 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,24)])
|
||||
|
||||
#anodes, plane 2 at +zmax/2
|
||||
offset = offset-3*k
|
||||
xarra_2 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)])
|
||||
yarra_2 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)])
|
||||
|
||||
#cathodes, plane2 at +zmax/2
|
||||
offsetc = offsetc-3*kc
|
||||
xarrc_2 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,24)])
|
||||
yarrc_2 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,24)])
|
||||
|
||||
wire_radius = 0.254 #mm
|
||||
anode_wires = []
|
||||
cathode_wires = []
|
||||
aw_tags = [(3,i) for i in range(24)]
|
||||
cw_tags = [(3,i+24) for i in range(24)]
|
||||
|
||||
for i,[xa,ya,xc,yc,xa2,ya2,xc2,yc2] in enumerate(zip(xarra_1,yarra_1,xarrc_1,yarrc_1,xarra_2,yarra_2,xarrc_2,yarrc_2)):
|
||||
print(i,xa,ya,-178.3,xc,yc,-178.3,xa2,ya2,178.3,xc2,yc2,178.3)
|
||||
pa1 = gmsh.model.occ.addPoint(xa,ya,-178.3,lc)
|
||||
pa2 = gmsh.model.occ.addPoint(xa2,ya2,178.3,lc)
|
||||
pc1 = gmsh.model.occ.addPoint(xc,yc,-178.3,lc)
|
||||
pc2 = gmsh.model.occ.addPoint(xc2,yc2,178.3,lc)
|
||||
linea = gmsh.model.occ.addLine(pa1,pa2)
|
||||
linec = gmsh.model.occ.addLine(pc1,pc2)
|
||||
anode_wires.append(linea)
|
||||
cathode_wires.append(linec)
|
||||
|
||||
#anasen_barrel = gmsh.model.occ.addCylinder(0,0,-500,0,0,500+605,300,1234) #tag 1234
|
||||
anasen_barrel = gmsh.model.occ.addCylinder(0,0,-200, 0,0,400, 300) #tag 1234
|
||||
|
||||
gmsh.model.occ.synchronize()
|
||||
gmsh.model.mesh.embed(1,anode_wires+cathode_wires,3,anasen_barrel)
|
||||
|
||||
f1 = gmsh.model.mesh.field.add("Distance")
|
||||
gmsh.model.mesh.field.setNumbers(f1,"CurvesList",anode_wires+cathode_wires)
|
||||
|
||||
f2 = gmsh.model.mesh.field.add("Threshold")
|
||||
gmsh.model.mesh.field.setNumber(f2,"InField",f1)
|
||||
gmsh.model.mesh.field.setNumber(f2,"SizeMin",0.08)
|
||||
gmsh.model.mesh.field.setNumber(f2,"SizeMax",5)
|
||||
gmsh.model.mesh.field.setNumber(f2,"DistMin",1)
|
||||
gmsh.model.mesh.field.setNumber(f2,"DistMax",20)
|
||||
|
||||
gmsh.model.mesh.field.setAsBackgroundMesh(f2)
|
||||
|
||||
gmsh.model.addPhysicalGroup(1, anode_wires, tag=10)
|
||||
gmsh.model.setPhysicalName(1,10,"anode_wires")
|
||||
|
||||
gmsh.model.addPhysicalGroup(1, cathode_wires, tag=20)
|
||||
gmsh.model.setPhysicalName(1,20,"cathode_wires")
|
||||
|
||||
gmsh.option.setNumber("Mesh.Algorithm",6)
|
||||
gmsh.option.setNumber("Mesh.Algorithm3D", 10) # For 3D meshes
|
||||
|
||||
gmsh.model.mesh.generate(dim=3)
|
||||
gmsh.model.mesh.refine()
|
||||
#gmsh.model.mesh.refine()
|
||||
#gmsh.model.mesh.refine()
|
||||
gmsh.fltk.run()
|
||||
gmsh.finalize()
|
||||
62
anasen_fem/paraview_plotter.py
Executable file
62
anasen_fem/paraview_plotter.py
Executable file
|
|
@ -0,0 +1,62 @@
|
|||
#!~/ParaView-6.1.0-RC1-MPI-Linux-Python3.12-x86_64/bin/pvbatch
|
||||
import numpy as np
|
||||
import sys
|
||||
from paraview.simple import *
|
||||
|
||||
reader = XMLUnstructuredGridReader(FileName=["wires2d/elfield_anasen_t0001.vtu"])
|
||||
|
||||
contour_filter = Contour(Input=reader,ContourBy = 'potential')
|
||||
contour_filter.Isosurfaces = [i for i in np.arange(0,660,650/24.)]
|
||||
|
||||
renderView = GetActiveViewOrCreate('RenderView')
|
||||
renderView.ViewSize = [800,800]
|
||||
renderView.OrientationAxesVisibility = 0 # Hide axis
|
||||
renderView.UseColorPaletteForBackground=0
|
||||
renderView.Background = [0.1, 0.1, 0.1] # Set background to dark gray (RGB 0-1)
|
||||
|
||||
renderView.MultiSamples = 8 # 0 disables it, 4-8 is usually sufficient
|
||||
|
||||
ResetCamera()
|
||||
|
||||
contour_display = Show(contour_filter, renderView)
|
||||
|
||||
#colorbar
|
||||
contour_display_potentialLUT = GetColorTransferFunction('potential', contour_display, separate=True)
|
||||
contour_display_potentialLUT.ApplyPreset('Cool to Warm', True)
|
||||
contour_display.SetScalarBarVisibility(renderView, True)
|
||||
|
||||
#axesGrid = renderView.AxesGrid
|
||||
#axesGrid.Visibility = 1
|
||||
#axesGrid.XTitle = "x (mm)"
|
||||
#axesGrid.YTitle = "y (mm)"
|
||||
|
||||
# 1. Get the active view
|
||||
view = GetActiveView()
|
||||
|
||||
# 2. Define your desired coordinate ranges (x_min, x_max, y_min, y_max, z_min, z_max)
|
||||
# Example: Look at a box from -10 to 10 in all dimensions
|
||||
x_min, x_max = -50.0, 50.0
|
||||
y_min, y_max = -50.0, 50.0
|
||||
z_min, z_max = -50.0, 50.0
|
||||
|
||||
# 3. Calculate Center, Position, and Parallel Scale
|
||||
center = [(x_min + x_max) / 2.0, (y_min + y_max) / 2.0, (z_min + z_max) / 2.0]
|
||||
# Position the camera far away along Z to look at the center
|
||||
position = [center[0], center[1], z_min - 30.0]
|
||||
# Parallel scale defines how much of the scene is visible.
|
||||
# It is usually half the height of the viewed area.
|
||||
view.CameraParallelScale = max((x_max - x_min), (y_max - y_min))/1.6
|
||||
|
||||
# 4. Apply settings
|
||||
view.CenterOfRotation = center
|
||||
view.CameraPosition = position
|
||||
view.CameraFocalPoint = center
|
||||
view.CameraViewUp = [0.0, 1.0, 0.0] # Y-axis is up
|
||||
|
||||
# 5. Enable Parallel Projection (optional, often better for exact mapping)
|
||||
view.CameraParallelProjection = 1
|
||||
|
||||
#ResetCamera()
|
||||
Render()
|
||||
|
||||
SaveScreenshot("contour_output.png")
|
||||
15
anasen_fem/run.py
Executable file
15
anasen_fem/run.py
Executable file
|
|
@ -0,0 +1,15 @@
|
|||
import os
|
||||
|
||||
#val=-178.3
|
||||
val=17.83
|
||||
count=11
|
||||
while val<178.3+0.1:
|
||||
print(val)
|
||||
os.system("python3 wires_gmsh2d_bc.py "+str(val))
|
||||
os.system("ElmerGrid 14 2 wires2d.msh")
|
||||
os.system("ElmerSolver wires2d.sif")
|
||||
os.system("./paraview_plotter.py")
|
||||
os.system("cp contour_output.png contour_output_z_%02d_%1.4f.png"%(count,val))
|
||||
val=val+17.83
|
||||
count = count + 1
|
||||
|
||||
1
anasen_fem/scalars.dat
Executable file
1
anasen_fem/scalars.dat
Executable file
|
|
@ -0,0 +1 @@
|
|||
6.500000000000E+002
|
||||
8
anasen_fem/scalars.dat.names
Executable file
8
anasen_fem/scalars.dat.names
Executable file
|
|
@ -0,0 +1,8 @@
|
|||
Metadata for SaveScalars file: ./scalars.dat
|
||||
Elmer version: 26.1
|
||||
Elmer compilation date: 2026-03-05
|
||||
Solver input file: wires2d.sif
|
||||
File started at: 2026/03/11 23:33:58
|
||||
|
||||
Variables in columns of matrix:
|
||||
1: res: potential difference
|
||||
103
anasen_fem/wires2d.sif
Executable file
103
anasen_fem/wires2d.sif
Executable file
|
|
@ -0,0 +1,103 @@
|
|||
Check Keywords Warn
|
||||
|
||||
Header
|
||||
Mesh DB "." "wires2d"
|
||||
End
|
||||
|
||||
Simulation
|
||||
Coordinate System = Cartesian 2D
|
||||
Simulation Type = Steady State
|
||||
Steady State Max Iterations = 1
|
||||
Output File = "elstatics.result"
|
||||
Post File = "elstatics.ep"
|
||||
End
|
||||
|
||||
Constants
|
||||
Permittivity Of Vacuum = 8.8542e-12
|
||||
End
|
||||
|
||||
|
||||
Body 1
|
||||
Target Bodies(1) = 13
|
||||
Equation = 1
|
||||
Material = 1
|
||||
End
|
||||
|
||||
|
||||
Equation 1
|
||||
Active Solvers(2) = 1 2
|
||||
End
|
||||
|
||||
|
||||
Material 1
|
||||
Relative Permittivity = 1
|
||||
End
|
||||
|
||||
|
||||
Solver 1
|
||||
Equation = Electrostatics
|
||||
Procedure = "StatElecSolve" "StatElecSolver"
|
||||
|
||||
Variable = Potential
|
||||
Variable DOFs = 1
|
||||
|
||||
Calculate Electric Field = True
|
||||
Calculate Electric Flux = False
|
||||
|
||||
Linear System Solver = Iterative
|
||||
Linear System Iterative Method = CG
|
||||
Linear System Max Iterations = 5000
|
||||
Linear System Convergence Tolerance = 1.0e-8
|
||||
Linear System Preconditioning = ILU1
|
||||
Calculate Vectors = Logical True
|
||||
End
|
||||
|
||||
Solver 2
|
||||
Equation = "Electric Field"
|
||||
Procedure = "FluxSolver" "FluxSolver"
|
||||
|
||||
! Calculate from the potential
|
||||
Target Variable = "Potential"
|
||||
! Name of the output vector field in VTU
|
||||
Flux Variable = String "Electric Field"
|
||||
|
||||
! Use 2D components (x, y)
|
||||
Flux Coefficient = String "Permittivity"
|
||||
|
||||
Calculate Vectors = Logical True
|
||||
End
|
||||
|
||||
Solver 3
|
||||
Equation = Result Output
|
||||
Procedure = "ResultOutputSolve" "ResultOutputSolver"
|
||||
Output File Name = elfield_anasen ! Sets prefix for output files
|
||||
Output Format = Vtu
|
||||
! Optional: Select specific variables to save
|
||||
Scalar Field 1 = Potential
|
||||
Vector Field 1 = Electric Field
|
||||
End
|
||||
|
||||
Solver 4
|
||||
Exec Solver = After All
|
||||
Equation = SaveScalars
|
||||
Procedure = "SaveData" "SaveScalars"
|
||||
Filename = "scalars.dat"
|
||||
End
|
||||
|
||||
|
||||
Boundary Condition 1
|
||||
Target Boundaries = 10
|
||||
Potential = 650
|
||||
Calculate Electric Force = True
|
||||
End
|
||||
|
||||
|
||||
Boundary Condition 2
|
||||
Target Boundaries = 20
|
||||
Potential = 0
|
||||
End
|
||||
|
||||
!Boundary Condition 2
|
||||
! Target Boundaries = 30
|
||||
! Potential = 0
|
||||
!End
|
||||
127
anasen_fem/wires_gmsh2d_bc.py
Executable file
127
anasen_fem/wires_gmsh2d_bc.py
Executable file
|
|
@ -0,0 +1,127 @@
|
|||
import numpy as np
|
||||
import gmsh,sys
|
||||
|
||||
gmsh.initialize()
|
||||
gmsh.model.add("adaptive_mesh")
|
||||
gmsh.option.setNumber('General.NumThreads', 4)
|
||||
#gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfElements", 200000)
|
||||
#gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfNodes", 200000)
|
||||
#gmsh.option.setNumber("Mesh.Adapt.MaxIter",5)
|
||||
#gmsh.option.setNumber("Mesh.MeshSizeMin", 5e-3)
|
||||
#gmsh.option.setNumber("Mesh.MeshSizeMax", 10.0)
|
||||
gmsh.option.setNumber("Geometry.Tolerance", 4e-2)
|
||||
#gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
|
||||
|
||||
lc = 0.04
|
||||
#z_loc = -178.3
|
||||
|
||||
if len(sys.argv) < 2:
|
||||
print("Usage: python3 wires_gmsh2d_bc.py <z_locus in mm>")
|
||||
quit()
|
||||
|
||||
z_loc = float(sys.argv[1])
|
||||
|
||||
k=(2*np.pi/24.)
|
||||
|
||||
#anodes, plane 1 at -zmax/2
|
||||
k=-2*np.pi/24.
|
||||
offset = 6*k + 3*k #-pi/2
|
||||
xarra_1 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)])
|
||||
yarra_1 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)])
|
||||
|
||||
#cathodes, plane 1 at -zmax/2
|
||||
kc=2*np.pi/24.
|
||||
offsetc = -4*kc + 2*kc - np.pi/24 #-pi/4
|
||||
xarrc_1 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,24)])
|
||||
yarrc_1 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,24)])
|
||||
|
||||
#anodes, plane 2 at +zmax/2
|
||||
offset = offset-3*k
|
||||
xarra_2 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)])
|
||||
yarra_2 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)])
|
||||
|
||||
#cathodes, plane2 at +zmax/2
|
||||
offsetc = offsetc-3*kc
|
||||
xarrc_2 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,24)])
|
||||
yarrc_2 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,24)])
|
||||
|
||||
direction_anodes_x = xarra_2 - xarra_1
|
||||
direction_anodes_y = yarra_2 - yarra_1
|
||||
|
||||
direction_cathodes_x = xarrc_2 - xarrc_1
|
||||
direction_cathodes_y = yarrc_2 - yarrc_1
|
||||
|
||||
t = (z_loc+178.3)/(2*178.3) #z=-178.3 is 0, z=+178.3 is 1
|
||||
xloc_a = xarra_1 + t*direction_anodes_x
|
||||
yloc_a = yarra_1 + t*direction_anodes_y
|
||||
xloc_c = xarrc_1 + t*direction_cathodes_x
|
||||
yloc_c = yarrc_1 + t*direction_cathodes_y
|
||||
|
||||
wire_radius_a = 0.018 #mm
|
||||
wire_radius_c = 0.0762 #mm
|
||||
anode_wires = []
|
||||
cathode_wires = []
|
||||
aw_tags = [(3,i) for i in range(24)]
|
||||
cw_tags = [(3,i+24) for i in range(24)]
|
||||
|
||||
#for i,[xa,ya,xc,yc] in enumerate(zip(xarra_1,yarra_1,xarrc_1,yarrc_1)):
|
||||
for i,[xa,ya,xc,yc] in enumerate(zip(xloc_a,yloc_a,xloc_c,yloc_c)):
|
||||
print(i,xa,ya,-178.3,xc,yc,-178.3)
|
||||
adisk = gmsh.model.occ.addDisk(xa,ya,0,wire_radius_a,wire_radius_a)
|
||||
cdisk = gmsh.model.occ.addDisk(xc,yc,0,wire_radius_c,wire_radius_c)
|
||||
anode_wires.append(adisk)
|
||||
cathode_wires.append(cdisk)
|
||||
|
||||
anasen_barrel = gmsh.model.occ.addDisk(0,0,0,500,500)
|
||||
#gmsh.model.occ.synchronize()
|
||||
#gmsh.model.mesh.embed(1,anode_wires+cathode_wires,2,anasen_barrel)
|
||||
|
||||
gmsh.option.setNumber("Geometry.Tolerance", 1e-6)
|
||||
gmsh.option.setNumber("Geometry.OCCFixDegenerated", 1)
|
||||
gmsh.model.occ.synchronize()
|
||||
|
||||
awire_surfs = []
|
||||
for w in anode_wires:
|
||||
awire_surfs += [s[1] for s in gmsh.model.getBoundary([(2,w)], oriented=False) if s[0] == 1]
|
||||
|
||||
cwire_surfs = []
|
||||
for w in cathode_wires:
|
||||
cwire_surfs += [s[1] for s in gmsh.model.getBoundary([(2,w)], oriented=False) if s[0] == 1]
|
||||
gmsh.model.mesh.embed(1,cwire_surfs+awire_surfs,2,anasen_barrel)
|
||||
|
||||
for s in gmsh.model.getBoundary([(2,w)],oriented=False):
|
||||
if s[0] == 1:
|
||||
anasen_bdry=s[1]
|
||||
|
||||
|
||||
f1 = gmsh.model.mesh.field.add("Distance")
|
||||
gmsh.model.mesh.field.setNumbers(f1,"CurvesList",cwire_surfs+awire_surfs)
|
||||
|
||||
f2 = gmsh.model.mesh.field.add("Threshold")
|
||||
gmsh.model.mesh.field.setNumber(f2,"InField",f1)
|
||||
gmsh.model.mesh.field.setNumber(f2,"SizeMin",0.1)
|
||||
gmsh.model.mesh.field.setNumber(f2,"SizeMax",5)
|
||||
gmsh.model.mesh.field.setNumber(f2,"DistMin",1)
|
||||
gmsh.model.mesh.field.setNumber(f2,"DistMax",20)
|
||||
|
||||
gmsh.model.mesh.field.setAsBackgroundMesh(f2)
|
||||
|
||||
gmsh.model.addPhysicalGroup(1, awire_surfs, tag=10)
|
||||
gmsh.model.setPhysicalName(1,10,"anode_wires")
|
||||
|
||||
gmsh.model.addPhysicalGroup(1, cwire_surfs, tag=20)
|
||||
gmsh.model.setPhysicalName(1,20,"cathode_wires")
|
||||
|
||||
#gmsh.model.addPhysicalGroup(1, [anasen_bdry], tag=30)
|
||||
#gmsh.model.setPhysicalName(1,30,"barrel_boundary")
|
||||
|
||||
gmsh.model.addPhysicalGroup(2,[anasen_barrel],tag=13)
|
||||
gmsh.model.setPhysicalName(1,13,"gas")
|
||||
|
||||
gmsh.option.setNumber("Mesh.Algorithm", 6)
|
||||
|
||||
gmsh.model.mesh.generate(dim=2)
|
||||
gmsh.model.mesh.refine()
|
||||
gmsh.write("wires2d.msh")
|
||||
#gmsh.fltk.run()
|
||||
gmsh.finalize()
|
||||
Loading…
Reference in New Issue
Block a user