63 for (
int i = 0 ; i < 6 ; i++) {
64 h_tpc_rate[i] =
new TH1F(TString::Format(
"h_tpc_rate_%d", i),
"detector #", 8, 0., 8.);
67 h_mctpc_recoil[0] =
new TH3F(
"h_mctpc_recoil_He",
"Neutron recoil energy [MeV]", 13, -0.5, 12.5, 8, -0.5, 7.5, 1000, 0., 10.);
68 h_mctpc_recoilW[0] =
new TH3F(
"h_mctpc_recoil_w_He",
"Neutron recoil energy [MeV]", 13, -0.5, 12.5, 8, -0.5, 7.5, 1000, 0., 10.);
72 h_mctpc_recoil[1] =
new TH3F(
"h_mctpc_recoil_O",
"Neutron recoil energy [MeV]", 13, -0.5, 12.5, 8, -0.5, 7.5, 1000, 0., 10.);
73 h_mctpc_recoilW[1] =
new TH3F(
"h_mctpc_recoil_w_O",
"Neutron recoil energy [MeV]", 13, -0.5, 12.5, 8, -0.5, 7.5, 1000, 0., 10.);
77 h_mctpc_recoil[2] =
new TH3F(
"h_mctpc_recoil_C",
"Neutron recoil energy [MeV]", 13, -0.5, 12.5, 8, -0.5, 7.5, 1000, 0., 10.);
78 h_mctpc_recoilW[2] =
new TH3F(
"h_mctpc_recoil_w_C",
"Neutron recoil energy [MeV]", 13, -0.5, 12.5, 8, -0.5, 7.5, 1000, 0., 10.);
83 for (
int i = 0 ; i < 12 ; i++) {
84 h_mctpc_kinetic[i] =
new TH2F(TString::Format(
"h_mctpc_kinetic_%d", i),
"Neutron kin. energy [GeV]", 8, -0.5, 7.5, 1000, 0., 10.);
85 h_mctpc_kinetic_zoom[i] =
new TH2F(TString::Format(
"h_mctpc_kinetic_zoom_%d", i),
"Neutron kin. energy [MeV]", 8, -0.5, 7.5, 1000,
87 h_mctpc_tvp[i] =
new TH2F(TString::Format(
"h_mctpc_tvp_%d", i),
"theta v phi", 180, 0., 180., 360, -180., 180.);
88 h_mctpc_tvpW[i] =
new TH2F(TString::Format(
"h_mctpc_tvpW_%d", i),
"theta v phi weighted by kin", 180, 0., 180., 360, -180., 180.);
89 h_mctpc_zr[i] =
new TH2F(TString::Format(
"h_mctpc_zr_%d", i),
"r v z", 200, -400., 400., 200, 0., 400.);
96 for (
int i = 0 ; i < 8 ; i++) {
97 for (
int j = 0; j < 12; j++) {
98 h_Wtvp1[i][j] =
new TH2F(TString::Format(
"h_Wtvp1_%d_%d", i, j),
"Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
99 h_Wtvp2[i][j] =
new TH2F(TString::Format(
"h_Wtvp2_%d_%d", i, j),
"Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
100 h_Wevtrl1[i][j] =
new TH2F(TString::Format(
"h_Wevtrl1_%d_%d", i, j),
"Deposited energy [keV] v. track length [cm]", 2000, 0., 10000,
102 h_Wevtrl2[i][j] =
new TH2F(TString::Format(
"h_Wevtrl2_%d_%d", i, j),
"Deposited energy [keV] v. track length [cm]", 2000, 0., 10000,
107 for (
int i = 0 ; i < 8 ; i++) {
108 h_z[i] =
new TH1F(TString::Format(
"h_z_%d", i),
"Charged density per cm^2", 2000, 0.0, 20.0);
110 h_zr[i] =
new TH2F(TString::Format(
"h_zr_%d", i),
"Charged density vs z vs r", 100, 0, 20, 100, 0., 5.);
112 h_xy[i] =
new TH2F(TString::Format(
"h_xy_%d", i),
"Charged density vs y vs x", 100, -5., 5., 100, -5., 5.);
114 h_zx[i] =
new TH2F(TString::Format(
"h_zx_%d", i),
"Charged density vs x vs r", 100, 0, 20, 100, -5., 5.);
116 h_zy[i] =
new TH2F(TString::Format(
"h_zy_%d", i),
"Charged density vs y vs r", 100, 0, 20, 100, -5., 5.);
118 h_evtrl[i] =
new TH2F(TString::Format(
"h_evtrl_%d", i),
"Deposited energy [keV] v. track length [cm]", 2000, 0., 10000, 200, 0.,
120 h_evtrlb[i] =
new TH2F(TString::Format(
"h_evtrlb_%d", i),
"Deposited energy [keV] v. track length [cm]", 2000, 0., 10000, 200, 0.,
122 h_evtrlc[i] =
new TH2F(TString::Format(
"h_evtrlc_%d", i),
"Deposited energy [keV] v. track length [cm]", 2000, 0., 10000, 200, 0.,
124 h_evtrld[i] =
new TH2F(TString::Format(
"h_evtrld_%d", i),
"Deposited energy [keV] v. track length [cm]", 2000, 0., 10000, 200, 0.,
127 h_evtrl_p[i] =
new TH2F(TString::Format(
"h_evtrl_p_%d", i),
"Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200, 0.,
129 h_evtrl_x[i] =
new TH2F(TString::Format(
"h_evtrl_x_%d", i),
"Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200, 0.,
131 h_evtrl_Hex[i] =
new TH2F(TString::Format(
"h_evtrl_Hex_%d", i),
"Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200,
133 h_evtrl_He[i] =
new TH2F(TString::Format(
"h_evtrl_He_%d", i),
"Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200,
135 h_evtrl_C[i] =
new TH2F(TString::Format(
"h_evtrl_C_%d", i),
"Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200, 0.,
137 h_evtrl_O[i] =
new TH2F(TString::Format(
"h_evtrl_O_%d", i),
"Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200, 0.,
139 h_evtrl_He_pure[i] =
new TH2F(TString::Format(
"h_evtrl_He_pure_%d", i),
"Deposited energy [keV] v. track length [cm]", 2000, 0.,
143 h_tevtrl[i] =
new TH2F(TString::Format(
"h_tevtrl_%d", i),
"t: Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200, 0.,
145 h_tevtrl_p[i] =
new TH2F(TString::Format(
"h_tevtrl_p_%d", i),
"t: Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200,
148 h_tevtrl_x[i] =
new TH2F(TString::Format(
"h_tevtrl_x_%d", i),
"t: Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200,
151 h_tevtrl_Hex[i] =
new TH2F(TString::Format(
"h_tevtrl_Hex_%d", i),
"t: Deposited energy [keV] v. track length [cm]", 2000, 0., 2000,
154 h_tevtrl_He[i] =
new TH2F(TString::Format(
"h_tevtrl_He_%d", i),
"t: Deposited energy [keV] v. track length [cm]", 2000, 0., 2000,
157 h_tevtrl_C[i] =
new TH2F(TString::Format(
"h_tevtrl_C_%d", i),
"t: Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200,
160 h_tevtrl_O[i] =
new TH2F(TString::Format(
"h_tevtrl_O_%d", i),
"t: Deposited energy [keV] v. track length [cm]", 2000, 0., 2000, 200,
163 h_tevtrl_He_pure[i] =
new TH2F(TString::Format(
"h_tevtrl_He_pure_%d", i),
"t: Deposited energy [keV] v. track length [cm]", 2000,
167 h_tvp[i] =
new TH2F(TString::Format(
"h_tvp_%d", i),
"Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
168 h_wtvpb[i] =
new TH2F(TString::Format(
"h_wtvpb_%d", i),
"Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
169 h_wtvpc[i] =
new TH2F(TString::Format(
"h_wtvpc_%d", i),
"Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
170 h_wtvpd[i] =
new TH2F(TString::Format(
"h_wtvpd_%d", i),
"Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
172 h_tvpb[i] =
new TH2F(TString::Format(
"h_tvpb_%d", i),
"Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
173 h_tvpc[i] =
new TH2F(TString::Format(
"h_tvpc_%d", i),
"Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
174 h_tvpd[i] =
new TH2F(TString::Format(
"h_tvpd_%d", i),
"Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
175 h_ttvp[i] =
new TH2F(TString::Format(
"h_ttvp_%d", i),
"t: Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
176 h_wtvp[i] =
new TH2F(TString::Format(
"h_wtvp_%d", i),
"Phi [deg] v. theta [deg] - weighted", 180, 0., 180, 360, -180., 180.);
177 h_tvp_x[i] =
new TH2F(TString::Format(
"h_tvp_x_%d", i),
"Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
178 h_ttvp_x[i] =
new TH2F(TString::Format(
"h_ttvp_x_%d", i),
"t: Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
179 h_wtvp_x[i] =
new TH2F(TString::Format(
"h_wtvp_x_%d", i),
"Phi [deg] v. theta [deg] - weighted", 180, 0., 180, 360, -180., 180.);
180 h_tvp_p[i] =
new TH2F(TString::Format(
"h_tvp_p_%d", i),
"Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
181 h_ttvp_p[i] =
new TH2F(TString::Format(
"h_ttvp_p_%d", i),
"t: Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
182 h_wtvp_p[i] =
new TH2F(TString::Format(
"h_wtvp_p_%d", i),
"Phi [deg] v. theta [deg] - weighted", 180, 0., 180, 360, -180., 180.);
183 h_tvp_He[i] =
new TH2F(TString::Format(
"h_tvp_He_%d", i),
"Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
184 h_ttvp_He[i] =
new TH2F(TString::Format(
"h_ttvp_He_%d", i),
"t: Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
185 h_wtvp_He[i] =
new TH2F(TString::Format(
"h_wtvp_He_%d", i),
"Phi [deg] v. theta [deg] - weighted", 180, 0., 180, 360, -180., 180.);
186 h_tvp_Hex[i] =
new TH2F(TString::Format(
"h_tvp_Hex_%d", i),
"Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
187 h_ttvp_Hex[i] =
new TH2F(TString::Format(
"h_ttvp_Hex_%d", i),
"t: Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
188 h_wtvp_Hex[i] =
new TH2F(TString::Format(
"h_wtvp_Hex_%d", i),
"Phi [deg] v. theta [deg] - weighted", 180, 0., 180, 360, -180.,
190 h_tvp_He_pure[i] =
new TH2F(TString::Format(
"h_tvp_He_pure_%d", i),
"Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
191 h_ttvp_He_pure[i] =
new TH2F(TString::Format(
"h_ttvp_He_pure_%d", i),
"t: Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180.,
193 h_wtvp_He_pure[i] =
new TH2F(TString::Format(
"h_wtvp_He_pure_%d", i),
"Phi [deg] v. theta [deg] - weighted", 180, 0., 180, 360,
195 h_twtvp_He_pure[i] =
new TH2F(TString::Format(
"h_twtvp_He_pure_%d", i),
"t: Phi [deg] v. theta [deg] - weighted", 180, 0., 180,
198 h_tvp_C[i] =
new TH2F(TString::Format(
"h_tvp_C_%d", i),
"Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
199 h_ttvp_C[i] =
new TH2F(TString::Format(
"h_ttvp_C_%d", i),
"t: Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
200 h_wtvp_C[i] =
new TH2F(TString::Format(
"h_wtvp_C_%d", i),
"Phi [deg] v. theta [deg] - weighted", 180, 0., 180, 360, -180., 180.);
201 h_tvp_O[i] =
new TH2F(TString::Format(
"h_tvp_O_%d", i),
"Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
202 h_ttvp_O[i] =
new TH2F(TString::Format(
"h_ttvp_O_%d", i),
"t: Phi [deg] v. theta [deg]", 180, 0., 180, 360, -180., 180.);
203 h_wtvp_O[i] =
new TH2F(TString::Format(
"h_wtvp_O_%d", i),
"Phi [deg] v. theta [deg] - weighted", 180, 0., 180, 360, -180., 180.);
258 int ring_section = 0;
259 const int section_ordering[12] = {1, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2};
260 for (
const auto& sadMetaHit : sadMetaHits) {
261 rate = sadMetaHit.getrate();
262 double ss = sadMetaHit.getss() / 100.;
263 if (ss < 0) ss += 3000.;
264 int section = (int)(ss / 250.);
265 if (section >= 0 && section < 12) ring_section = section_ordering[section];
284 for (
int i = 0; i < 8; i++) {
299 auto phiArray =
new vector<double>[8]();
300 auto thetaArray =
new vector<double>[8]();
301 auto pidArray =
new vector<int>[8]();
303 auto esumArray =
new vector<double>[8]();
304 auto trlArray =
new vector<double>[8]();
306 ROOT::Math::XYZVector EndPoint;
308 for (
int i = 0; i < nSimHits; i++) {
311 ROOT::Math::XYZVector position = aHit->
gettkPos();
312 double xpos = position.X() / 100. -
TPCCenter[detNb].X();
313 double ypos = position.Y() / 100. -
TPCCenter[detNb].Y();
314 double zpos = position.Z() / 100. -
TPCCenter[detNb].Z() +
m_z_DG / 2.;
315 if (0. < zpos && zpos <
m_z_DG) {
318 if (PDGid == 1000080160)
ORec[detNb] =
true;
319 if (PDGid == 1000060120)
CRec[detNb] =
true;
320 if (PDGid == 1000020040)
HeRec[detNb] =
true;
326 double ioni = (edep - niel) * 1e3;
328 double r =
sqrt(xpos * xpos + ypos * ypos);
329 h_z[detNb]->Fill(zpos, ioni);
330 h_zr[detNb]->Fill(zpos, r, ioni);
331 h_zx[detNb]->Fill(zpos, xpos, ioni);
332 h_xy[detNb]->Fill(xpos, ypos, ioni);
333 h_zy[detNb]->Fill(zpos, ypos, ioni);
334 ROOT::Math::XYZVector direction = aHit->
gettkMomDir();
335 double theta = direction.Theta() * TMath::RadToDeg();
336 double phi = direction.Phi() * TMath::RadToDeg();
340 (0. < zpos && zpos <
m_z_DG)) {
348 if (esum[detNb] > 0) {
349 esumArray[detNb].push_back(esum[detNb]);
350 ROOT::Math::XYZVector BeginPoint;
351 BeginPoint.SetXYZ(xpos, ypos, zpos);
352 double trl0 = BeginPoint.Dot(direction.Unit());
353 double trl1 = EndPoint.Dot(direction.Unit());
354 trlArray[detNb].push_back(fabs(trl0 - trl1));
390 thetaArray[detNb].push_back(theta);
391 phiArray[detNb].push_back(phi);
392 pidArray[detNb].push_back(PDGid);
400 EndPoint.SetXYZ(xpos, ypos, zpos);
449 for (
const auto& mcpart : mcparts) {
450 const double energy = mcpart.getEnergy();
451 const double mass = mcpart.getMass();
452 double kin = energy - mass;
453 const double PDG = mcpart.getPDG();
454 const ROOT::Math::XYZVector vtx = mcpart.getProductionVertex();
455 const ROOT::Math::XYZVector mom = mcpart.getMomentum();
456 double theta = mom.Theta() * TMath::RadToDeg();
457 double phi = mom.Phi() * TMath::RadToDeg();
459 double r =
sqrt(vtx.X() * vtx.X() + vtx.Y() * vtx.Y());
460 int partID[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
461 if (trID == mcpart.getTrackID())
continue;
462 else trID = mcpart.getTrackID();
465 for (
const auto& sHit : SimHits) {
466 if (sHit.gettkID() == trID) {
467 detNb = sHit.getdetNb();
468 kin = sHit.gettkKEnergy() / 1000;
474 double trlen = abs(2. / TMath::Sin(mom.Theta()));
475 if (trlen > 10) trlen = 10.;
478 double recoil = gRandom->Uniform(fract) * kin * 1e3;
479 double weight =
m_intProb[irecoil]->Eval(kin * 1e3) * trlen;
480 if (weight < 0) weight = 0;
494 else if (PDG == 1000080160) partID[5] = 1;
495 else if (PDG == 1000060120) partID[6] = 1;
496 else if (PDG == 1000020040) partID[7] = 1;
517 for (
int i = 0; i < 9; i++) {
518 if (partID[i] == 1) {
532 for (
const auto& aTrack : Tracks) {
533 const int detNb = aTrack.getdetNb();
534 const float phi = aTrack.getphi();
535 const float theta = aTrack.gettheta();
536 const float trl = aTrack.gettrl();
537 const float tesum = aTrack.getesum();
538 const int pixnb = aTrack.getpixnb();
541 for (
int j = 0; j < 16; j++) {
543 side[j] = aTrack.getside()[j];
545 Bool_t EdgeCuts =
false;
546 if (side[0] == 0 && side[1] == 0 && side[2] == 0 && side[3] == 0) EdgeCuts =
true;
547 Bool_t Asource =
false;
548 if (side[4] == 2 && side[5] == 2) Asource =
true;
555 for (
int j = 0; j < 6; j++) partID[j + 1] = aTrack.getpartID()[j];
573 h_evtrl[detNb]->Fill(tesum, trl);
574 h_tvp[detNb]->Fill(theta, phi);
575 h_wtvp[detNb]->Fill(theta, phi, tesum);
576 h_Wtvp1[detNb][0]->Fill(theta, phi, rate);
577 h_Wevtrl1[detNb][0]->Fill(tesum, trl, rate);
578 h_Wtvp2[detNb][0]->Fill(theta, phi, rate * tesum);
580 if (EdgeCuts && pixnb > 10. && tesum > 10.) {
582 h_tvpb[detNb]->Fill(theta, phi);
583 h_wtvpb[detNb]->Fill(theta, phi, tesum);
584 h_Wtvp1[detNb][1]->Fill(theta, phi, rate);
585 h_Wevtrl1[detNb][1]->Fill(tesum, trl, rate);
586 h_Wtvp2[detNb][1]->Fill(theta, phi, rate * tesum);
589 for (
int j = 0; j < 7; j++) {
590 if (j == 3 && !EdgeCuts && (partID[1] == 1 || partID[2] == 1 || partID[4] == 1 || partID[5] == 1 || partID[6] == 1)) partID[j] = 0;
591 if ((j == 4 || j == 5) && !Asource) partID[j] = 0;
592 if (partID[j] == 1) {
593 h_Wtvp1[detNb][2 + j]->Fill(theta, phi, rate);
594 h_Wevtrl1[detNb][2 + j]->Fill(tesum, trl, rate);
595 h_Wtvp2[detNb][2 + j]->Fill(theta, phi, rate * tesum);
598 h_tvpc[detNb]->Fill(theta, phi);
599 h_wtvpc[detNb]->Fill(theta, phi, tesum);
603 h_tvpd[detNb]->Fill(theta, phi);
604 h_wtvpd[detNb]->Fill(theta, phi, tesum);
608 h_tvp_x[detNb]->Fill(theta, phi);
609 h_wtvp_x[detNb]->Fill(theta, phi, tesum);
613 h_tvp_p[detNb]->Fill(theta, phi);
614 h_wtvp_p[detNb]->Fill(theta, phi, tesum);
618 h_tvp_x[detNb]->Fill(theta, phi);
619 h_wtvp_x[detNb]->Fill(theta, phi, tesum);
624 h_wtvp_He[detNb]->Fill(theta, phi, tesum);
639 delete [] thetaArray;