LMU/musrsim
LMU/musrsim
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fixed some bugs with the code and included better handeling for Si Pixel detectors.

Browse Source
This commit is contained in:
salman 2023-02-02 09:01:37 +01:00
parent 91f145b2fb
commit 5399b0d6a6
2 changed files with 56 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

92
run/SiPixel/SiSpect.C
View File

@ -1,3 +1,17 @@
/*
Particle codes:
MuonPlus= -13;
MuonMinus= 13;
Electron= 11;
Positron= -11;
Gamma= 22;
Proton= 2212;
Neutron= 2112;
PionMinus= -211;
PionPlus= 211;
*/
#define SiSpect_cxx #define SiSpect_cxx
#include "SiSpect.h" #include "SiSpect.h"
#include <TH2.h> #include <TH2.h>
@ -85,20 +99,19 @@ void SiSpect::CreateIO( Bool_t FigFlag, Double_t eCut )
// } // }
hTof->Fill(muTargetTime); hTof->Fill(muTargetTime);
if (muTargetTime>0.23) tofFlag = 1; //if (muTargetTime>0.23) tofFlag = 1;
//tofFlag = 1.; //tofFlag = 1.;
for (Int_t i=0; i<det_n; i++) { for (Int_t i=0; i<det_n; i++) {
// Only positrons // Only muons
if (det_ID[i] == 501 && save_particleID[i]==-13) { if (det_ID[i] == 501 && det_VrtxParticleID[i]==-13) {
hEdepositMuI->Fill(det_edep[i]); hEdepositMuI->Fill(det_edep[i]);
break; //fill only once break; //fill only once
} }
} }
for (Int_t i=0; i<det_n; i++) { for (Int_t i=0; i<det_n; i++) {
// Only positrons // Only muons
if (det_ID[i] == 601 && save_particleID[i]==-13) { if (det_ID[i] == 601 && det_VrtxParticleID[i]==-13) {
hEdepositMuO->Fill(det_edep[i]); hEdepositMuO->Fill(det_edep[i]);
break; //fill only once break; //fill only once
} }
@ -108,7 +121,10 @@ void SiSpect::CreateIO( Bool_t FigFlag, Double_t eCut )
// Hist in Back I detector (501) // Hist in Back I detector (501)
for (Int_t i=0; i<det_n; i++) { for (Int_t i=0; i<det_n; i++) {
// Only positrons // Only positrons
if (det_ID[i] == 501 && save_particleID[i]==22) { if (jentry > 1000 && jentry < 1100) {
printf("entry = %lld det_ID = %d det_Vrtx = %d\n",jentry,det_ID[i],det_VrtxParticleID[i]);
}
if (det_ID[i] == 501 && det_VrtxParticleID[i]==-11) {
hEdeposited->Fill(det_edep[i]); hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){ if (det_edep[i]>eCut){
hBackI->Fill(det_time_start[i]); hBackI->Fill(det_time_start[i]);
@ -122,7 +138,7 @@ void SiSpect::CreateIO( Bool_t FigFlag, Double_t eCut )
// Hist in Forw I detector (502) // Hist in Forw I detector (502)
for (Int_t i=0; i<det_n; i++) { for (Int_t i=0; i<det_n; i++) {
if (det_ID[i] == 502 && save_particleID[i]==22) { if (det_ID[i] == 502 && det_VrtxParticleID[i]==-11) {
hEdeposited->Fill(det_edep[i]); hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){ if (det_edep[i]>eCut){
hForwI->Fill(det_time_start[i]); hForwI->Fill(det_time_start[i]);
@ -136,7 +152,7 @@ void SiSpect::CreateIO( Bool_t FigFlag, Double_t eCut )
// Hist in Top I detector (503) // Hist in Top I detector (503)
for (Int_t i=0; i<det_n; i++) { for (Int_t i=0; i<det_n; i++) {
if (det_ID[i] == 503 && save_particleID[i]==22) { if (det_ID[i] == 503 && det_VrtxParticleID[i]==-11) {
hEdeposited->Fill(det_edep[i]); hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){ if (det_edep[i]>eCut){
hTopI->Fill(det_time_start[i]); hTopI->Fill(det_time_start[i]);
@ -150,7 +166,7 @@ void SiSpect::CreateIO( Bool_t FigFlag, Double_t eCut )
// Hist in Down I detector (504) // Hist in Down I detector (504)
for (Int_t i=0; i<det_n; i++) { for (Int_t i=0; i<det_n; i++) {
if (det_ID[i] == 504 && save_particleID[i]==22) { if (det_ID[i] == 504 && det_VrtxParticleID[i]==-11) {
hEdeposited->Fill(det_edep[i]); hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){ if (det_edep[i]>eCut){
hDownI->Fill(det_time_start[i]); hDownI->Fill(det_time_start[i]);
@ -164,7 +180,7 @@ void SiSpect::CreateIO( Bool_t FigFlag, Double_t eCut )
// Hist in Back O detector (601) // Hist in Back O detector (601)
for (Int_t i=0; i<det_n; i++) { for (Int_t i=0; i<det_n; i++) {
if (det_ID[i] == 601 && save_particleID[i]==22) { if (det_ID[i] == 601 && det_VrtxParticleID[i]==-11) {
hEdeposited->Fill(det_edep[i]); hEdeposited->Fill(det_edep[i]);
// Only positrons // Only positrons
if (det_edep[i]>eCut){ if (det_edep[i]>eCut){
@ -179,7 +195,7 @@ void SiSpect::CreateIO( Bool_t FigFlag, Double_t eCut )
// Hist in Forw O detector (602) // Hist in Forw O detector (602)
for (Int_t i=0; i<det_n; i++) { for (Int_t i=0; i<det_n; i++) {
if (det_ID[i] == 602 && save_particleID[i]==22) { if (det_ID[i] == 602 && det_VrtxParticleID[i]==-11) {
hEdeposited->Fill(det_edep[i]); hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){ if (det_edep[i]>eCut){
hForwO->Fill(det_time_start[i]); hForwO->Fill(det_time_start[i]);
@ -193,7 +209,7 @@ void SiSpect::CreateIO( Bool_t FigFlag, Double_t eCut )
// Hist in Top O detector (603) // Hist in Top O detector (603)
for (Int_t i=0; i<det_n; i++) { for (Int_t i=0; i<det_n; i++) {
if (det_ID[i] == 603 && save_particleID[i]==22) { if (det_ID[i] == 603 && det_VrtxParticleID[i]==-11) {
hEdeposited->Fill(det_edep[i]); hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){ if (det_edep[i]>eCut){
hTopO->Fill(det_time_start[i]); hTopO->Fill(det_time_start[i]);
@ -207,7 +223,7 @@ void SiSpect::CreateIO( Bool_t FigFlag, Double_t eCut )
// Hist in Down O detector (604) // Hist in Down O detector (604)
for (Int_t i=0; i<det_n; i++) { for (Int_t i=0; i<det_n; i++) {
if (det_ID[i] == 604 && save_particleID[i]==22) { if (det_ID[i] == 604 && det_VrtxParticleID[i]==-11) {
hEdeposited->Fill(det_edep[i]); hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){ if (det_edep[i]>eCut){
hDownO->Fill(det_time_start[i]); hDownO->Fill(det_time_start[i]);
@ -312,13 +328,13 @@ void SiSpect::CoinIO( Bool_t FigFlag, Double_t eCut )
TH1F* hEdepTrig = new TH1F("hEdepTrig", "Radioactive electron kin. energy",250,0.,2.5); TH1F* hEdepTrig = new TH1F("hEdepTrig", "Radioactive electron kin. energy",250,0.,2.5);
TH1F* hEdepoTest = new TH1F("hEdepoTest", "Number of events in coincidence", 100,0.,1.); TH1F* hEdepoTest = new TH1F("hEdepoTest", "Number of events in coincidence", 100,0.,1.);
TH1F* hTof = new TH1F("hTof", "time-of-flight (#mus)", 1000, 0., 1.); TH1F* hTof = new TH1F("hTof", "time-of-flight (#mus)", 1000, 0., 1.);
TH2F* hBeamSpot = new TH2F("hBeamSpot", " x,y", 40, -40., 40., 40, -40., 40.); TH2F* hBeamSpot = new TH2F("hBeamSpot", "Beam on sample x,y", 40, -40., 40., 40, -40., 40.);
TH1F* hDetz = new TH1F("hDetz", "z detector [mm]", 100, -40., 40.); TH1F* hDetz = new TH1F("hDetz", "z detector [mm]", 100, -40., 40.);
TH1F* hEventID = new TH1F("hEventID", "Event ID", 10001, -0.5, 10000.5); TH1F* hEventID = new TH1F("hEventID", "Event ID", 10001, -0.5, 10000.5);
TH2F* hMuIxy = new TH2F("hMuIxy", " x,y", 40, -40., 40., 40, -40., 40.); TH2F* hMuIxy = new TH2F("hMuIxy", "Beam in inner x,y", 40, -40., 40., 40, -40., 40.);
TH2F* hMuOxy = new TH2F("hMuOxy", " x,y", 40, -40., 40., 40, -40., 40.); TH2F* hMuOxy = new TH2F("hMuOxy", "Beam on outer x,y", 40, -40., 40., 40, -40., 40.);
TH2F* hMuIOxy = new TH2F("hMuIOxy", " x,y", 40, -40., 40., 40, -40., 40.); TH2F* hMuIOxy = new TH2F("hMuIOxy", "Diff. between reconstucted-actual Mu position", 40, -40., 40., 40, -40., 40.);
TH2F* hMuSamxy = new TH2F("hMuIOxy", " x,y", 40, -40., 40., 40, -40., 40.); TH2F* hMuSamxy = new TH2F("hMuIOxy", "Beam extrapolated on sample x,y", 40, -40., 40., 40, -40., 40.);
// Back histogram, i.e. all counts in coincedence 501-601 // Back histogram, i.e. all counts in coincedence 501-601
TH1F* hBack = new TH1F("hBack","Back (#mus)",130,0.,13.); TH1F* hBack = new TH1F("hBack","Back (#mus)",130,0.,13.);
@ -344,7 +360,7 @@ void SiSpect::CoinIO( Bool_t FigFlag, Double_t eCut )
Long64_t nbytes = 0, nb = 0; Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) { for (Long64_t jentry=0; jentry<nentries;jentry++) {
tofFlag = 0; tofFlag = 1;
/* Long64_t ientry = LoadTree(jentry); /* Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;*/ if (ientry < 0) break;*/
nb = fChain->GetEntry(jentry); nbytes += nb; nb = fChain->GetEntry(jentry); nbytes += nb;
@ -358,23 +374,20 @@ void SiSpect::CoinIO( Bool_t FigFlag, Double_t eCut )
// } // }
hTof->Fill(muTargetTime); hTof->Fill(muTargetTime);
if (muTargetTime>0.23) tofFlag = 1; //if (muTargetTime>0.23) tofFlag = 1;
//tofFlag = 1.; //tofFlag = 1.;
// Hist in Muon detector (501) coincidence (601) // Hist in Muon detector (501) coincidence (601)
for (Int_t i=0; i<det_n; i++) { for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==501 && save_particleID[i]==-13) { if (det_ID[i]==601 && det_VrtxParticleID[i]==-13) {
// Mu outer hit // Mu outer hit
xo=det_x[i]; xo=det_x[i];
yo=det_y[i]; yo=det_y[i];
zo=det_z[i]; zo=det_z[i];
hMuOxy->Fill(xo, yo); hMuOxy->Fill(xo, yo);
// printf("Outer: %6.2f\t%6.2f\n",xo,yo);
hEdepositMu->Fill(det_edep[i]); hEdepositMu->Fill(det_edep[i]);
for (Int_t j=0; j<det_n; j++){ for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==601 && save_particleID[j]==-13) { if (det_ID[j]==501 && det_VrtxParticleID[j]==-13) {
printf("Inner %d\t%d\n",det_ID[j] ,save_particleID[j]);
//printf("Inner\n");
// Mu inner hit // Mu inner hit
xi=det_x[j]; xi=det_x[j];
yi=det_y[j]; yi=det_y[j];
@ -387,23 +400,22 @@ void SiSpect::CoinIO( Bool_t FigFlag, Double_t eCut )
// Propagate Mu trajectory on sample // Propagate Mu trajectory on sample
xs=xi+((xi-xo)/(zi-zo))*(zs-zi); xs=xi+((xi-xo)/(zi-zo))*(zs-zi);
ys=yi+((yi-yo)/(zi-zo))*(zs-zi); ys=yi+((yi-yo)/(zi-zo))*(zs-zi);
printf("Outer: (%6.2g,%6.2g) - Inner: (%6.2g,%6.2g) - Extr. sample: (%6.2g,%6.2g) - Real sample: (%6.2g,%6.2g)\n",xo,yo,xi,yi,xs,ys,save_x[j],save_y[j]);
hMuIOxy->Fill(xs-save_x[j], ys-save_y[j]); hMuIOxy->Fill(xs-save_x[j], ys-save_y[j]);
hMuSamxy->Fill(save_x[j], save_y[j]); hMuSamxy->Fill(xs, ys);
// printf("%6.2f\t%6.2f\n",xs,ys);
break; //fill only once break; //fill only once
} }
} }
} }
//printf("\n");
} }
// Hist in Back detector (501) coincidence (601) // Hist in Back detector (501) coincidence (601)
for (Int_t i=0; i<det_n; i++) { for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==501 && save_particleID[i]==22) { if (det_ID[i]==601 && det_VrtxParticleID[i]==-11) {
hEdeposited->Fill(det_edep[i]); hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){ if (det_edep[i]>eCut){
for (Int_t j=0; j<det_n; j++){ for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==601 && save_particleID[j]==22) { if (det_ID[j]==501 && det_VrtxParticleID[j]==-11) {
hEdeposited->Fill(det_edep[j]); hEdeposited->Fill(det_edep[j]);
if (det_edep[j]>eCut){ if (det_edep[j]>eCut){
hBack->Fill(det_time_start[j]); hBack->Fill(det_time_start[j]);
@ -419,11 +431,11 @@ void SiSpect::CoinIO( Bool_t FigFlag, Double_t eCut )
// Hist in Forw detector (502) coincidence (602) // Hist in Forw detector (502) coincidence (602)
for (Int_t i=0; i<det_n; i++) { for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==502 && save_particleID[i]==22) { if (det_ID[i]==602 && det_VrtxParticleID[i]==-11) {
hEdeposited->Fill(det_edep[i]); hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){ if (det_edep[i]>eCut){
for (Int_t j=0; j<det_n; j++){ for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==602 && save_particleID[j]==22) { if (det_ID[j]==502 && det_VrtxParticleID[j]==-11) {
if (det_edep[j]>eCut){ if (det_edep[j]>eCut){
hForw->Fill(det_time_start[j]); hForw->Fill(det_time_start[j]);
hDetz->Fill(det_z[j]); hDetz->Fill(det_z[j]);
@ -438,11 +450,11 @@ void SiSpect::CoinIO( Bool_t FigFlag, Double_t eCut )
// Hist in Top detector (503) coincidence (603) // Hist in Top detector (503) coincidence (603)
for (Int_t i=0; i<det_n; i++) { for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==503 && save_particleID[i]==22) { if (det_ID[i]==603 && det_VrtxParticleID[i]==-11) {
hEdeposited->Fill(det_edep[i]); hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){ if (det_edep[i]>eCut){
for (Int_t j=0; j<det_n; j++){ for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==603 && save_particleID[j]==22) { if (det_ID[j]==503 && det_VrtxParticleID[j]==-11) {
if (det_edep[j]>eCut){ if (det_edep[j]>eCut){
hTop->Fill(det_time_start[j]); hTop->Fill(det_time_start[j]);
hDetz->Fill(det_z[j]); hDetz->Fill(det_z[j]);
@ -457,11 +469,11 @@ void SiSpect::CoinIO( Bool_t FigFlag, Double_t eCut )
// Hist in Down detector (504) coincidence (604) // Hist in Down detector (504) coincidence (604)
for (Int_t i=0; i<det_n; i++) { for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==504 && save_particleID[i]==22) { if (det_ID[i]==604 && det_VrtxParticleID[i]==-11) {
hEdeposited->Fill(det_edep[i]); hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){ if (det_edep[i]>eCut){
for (Int_t j=0; j<det_n; j++){ for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==604 && save_particleID[j]==22) { if (det_ID[j]==504 && det_VrtxParticleID[j]==-11) {
if (det_edep[j]>eCut){ if (det_edep[j]>eCut){
hDown->Fill(det_time_start[j]); hDown->Fill(det_time_start[j]);
hDetz->Fill(det_z[j]); hDetz->Fill(det_z[j]);
@ -538,7 +550,8 @@ void SiSpect::CoinIO( Bool_t FigFlag, Double_t eCut )
c1->cd(3); c1->cd(3);
hDetz->Draw(); hDetz->Draw();
c1->cd(4); c1->cd(4);
hMu->Draw(); hBeamSpot->Draw();
hBeamSpot->Draw("cont0 same");
c1->cd(5); c1->cd(5);
hEdeposited->Draw(); hEdeposited->Draw();
hEdepositMu->Draw("sames"); hEdepositMu->Draw("sames");
@ -553,6 +566,9 @@ void SiSpect::CoinIO( Bool_t FigFlag, Double_t eCut )
c1->cd(7); c1->cd(7);
hMuSamxy->Draw(); hMuSamxy->Draw();
hMuSamxy->Draw("cont0 same"); hMuSamxy->Draw("cont0 same");
c1->cd(8);
c1->cd(9);
} else { } else {
hAsyFB -> Fit("pol0","NQ","",0.1, 13.); hAsyFB -> Fit("pol0","NQ","",0.1, 13.);
} }

2
run/SiPixel/SiSpect.h
View File

@ -58,6 +58,7 @@ public :
Double_t det_x[100]; //[det_n] Double_t det_x[100]; //[det_n]
Double_t det_y[100]; //[det_n] Double_t det_y[100]; //[det_n]
Double_t det_z[100]; //[det_n] Double_t det_z[100]; //[det_n]
Int_t det_VrtxParticleID[100]; //[det_n]
Int_t save_n; Int_t save_n;
Int_t save_detID[100]; //[save_n] Int_t save_detID[100]; //[save_n]
Int_t save_particleID[100]; //[save_n] Int_t save_particleID[100]; //[save_n]
@ -226,6 +227,7 @@ void SiSpect::Init(TTree *tree)
fChain->SetBranchAddress("det_x",det_x); fChain->SetBranchAddress("det_x",det_x);
fChain->SetBranchAddress("det_y",det_y); fChain->SetBranchAddress("det_y",det_y);
fChain->SetBranchAddress("det_z",det_z); fChain->SetBranchAddress("det_z",det_z);
fChain->SetBranchAddress("det_VrtxParticleID",det_VrtxParticleID);
fChain->SetBranchAddress("save_n", &save_n, &b_save_n); fChain->SetBranchAddress("save_n", &save_n, &b_save_n);
fChain->SetBranchAddress("save_detID", save_detID, &b_save_detID); fChain->SetBranchAddress("save_detID", save_detID, &b_save_detID);
fChain->SetBranchAddress("save_particleID", save_particleID, &b_save_particleID); fChain->SetBranchAddress("save_particleID", save_particleID, &b_save_particleID);