LMU/musrsim
LMU/musrsim
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
musrsim/run/SiPixel/SiSpect.C

587 lines
18 KiB
C
Raw Blame History

#define SiSpect_cxx
#include "SiSpect.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <string>
#include <TROOT.h>
#include <TKey.h>
#include <TFile.h>
#include <TSystem.h>
#include <TTree.h>
/*
To use
TFile* fin=new TFile("data/musr_1021.root");
.L SiSpect.C
SiSpect t
t.CreateIO(1,0.0)
*/
void SiSpect::Loop()
{
}
void SiSpect::CreateIO( Bool_t FigFlag, Double_t eCut )
{
// Double_t eCut = 0.7; //minimum energy deposition in scintillators [MeV]
Bool_t tofFlag = 0;
if (fChain == 0) return;
//DEFINE HISTOGRAMS
TH1F* hEdeposited = new TH1F("hEdeposited","Energy spectrum; E [MeV]", 250,0.,0.5);
TH1F* hEdepositMuI = new TH1F("hEdepositMuI","Energy spectrum Mu/Pos; E [MeV]", 250,0.,0.5);
TH1F* hEdepositMuO = new TH1F("hEdepositMuO","Energy spectrum Mu/Pos; E [MeV]", 250,0.,0.5);
TH1F* hEdepositCF = new TH1F("hEdepositCF","Energy spectrum; E [keV]", 500,0.,20.0);
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* hTof = new TH1F("hTof", "time-of-flight (#mus)", 1000, 0., 1.);
TH2F* hBeamSpot = new TH2F("hBeamSpot", " x,y", 40, -40., 40., 40, -40., 40.);
TH1F* hDetz = new TH1F("hDetz", "z detector [mm]", 100, -40., 40.);
TH1F* hDetzI = new TH1F("hDetzI", "z I detector [mm]", 100, -40., 40.);
TH1F* hDetzO = new TH1F("hDetzO", "z O detector [mm]", 100, -40., 40.);
TH1F* hEventID = new TH1F("hEventID", "Event ID", 10001, -0.5, 10000.5);
// Back inner histogram (501)
TH1F* hBackI = new TH1F("hBackI","Back I (#mus)",1300,0.,13.);
// Forward inner histogram (502)
TH1F* hForwI = new TH1F("hForwI","Forw I (#mus)",1300,0.,13.);
// Top inner histogram (503)
TH1F* hTopI = new TH1F("hTopI","Top I (#mus)",1300,0.,13.);
// Down inner histogram (504)
TH1F* hDownI = new TH1F("hDownI","Down I (#mus)",1300,0.,13.);
// Back outer histogram (601)
TH1F* hBackO = new TH1F("hBackO","Back O (#mus)",1300,0.,13.);
// Forward outer histogram (602)
TH1F* hForwO = new TH1F("hForwO","Forw O (#mus)",1300,0.,13.);
// Top outer histogram (603)
TH1F* hTopO = new TH1F("hTopO","Top O (#mus)",1300,0.,13.);
// Down outer histogram (604)
TH1F* hDownO = new TH1F("hDownO","Down O (#mus)",1300,0.,13.);
Long64_t NFI,NFO,NBI,NBO,NTI,NTO,NDI,NDO;
hEdeposited->Sumw2();
hEdepositCF->Sumw2();
hEdepTrig->Sumw2();
Long64_t nentries = fChain->GetEntriesFast();
//nentries=1000;
printf("nentries = %lld\n",nentries);
Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
tofFlag = 0;
/* Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;*/
nb = fChain->GetEntry(jentry); nbytes += nb;
// FILL IN HISTOGRAMS
hEdepositCF->Fill(save_ke[0]/1000.);
hBeamSpot->Fill(save_x[0], save_y[0]);
//for (Int_t i=0; i<det_n; i++)
// { if (det_ID[i]==623) { hEdeposited->Fill(save_ke[i]);}
// }
hTof->Fill(muTargetTime);
if (muTargetTime>0.23) tofFlag = 1;
//tofFlag = 1.;
for (Int_t i=0; i<det_n; i++) {
// Only positrons
if (det_ID[i] == 501 && save_particleID[i]==-13) {
hEdepositMuI->Fill(det_edep[i]);
break; //fill only once
}
}
for (Int_t i=0; i<det_n; i++) {
// Only positrons
if (det_ID[i] == 601 && save_particleID[i]==-13) {
hEdepositMuO->Fill(det_edep[i]);
break; //fill only once
}
}
// Hist in Back I detector (501)
for (Int_t i=0; i<det_n; i++) {
// Only positrons
if (det_ID[i] == 501 && save_particleID[i]==22) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
hBackI->Fill(det_time_start[i]);
hDetz->Fill(det_z[i]);
hDetzI->Fill(det_z[i]);
hEventID->Fill(eventID);
break; //fill only once
}
}
}
// Hist in Forw I detector (502)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i] == 502 && save_particleID[i]==22) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
hForwI->Fill(det_time_start[i]);
hDetz->Fill(det_z[i]);
hDetzI->Fill(det_z[i]);
hEventID->Fill(eventID);
break; //fill only once
}
}
}
// Hist in Top I detector (503)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i] == 503 && save_particleID[i]==22) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
hTopI->Fill(det_time_start[i]);
hDetz->Fill(det_z[i]);
hDetzI->Fill(det_z[i]);
hEventID->Fill(eventID);
break; //fill only once
}
}
}
// Hist in Down I detector (504)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i] == 504 && save_particleID[i]==22) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
hDownI->Fill(det_time_start[i]);
hDetz->Fill(det_z[i]);
hDetzI->Fill(det_z[i]);
hEventID->Fill(eventID);
break; //fill only once
}
}
}
// Hist in Back O detector (601)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i] == 601 && save_particleID[i]==22) {
hEdeposited->Fill(det_edep[i]);
// Only positrons
if (det_edep[i]>eCut){
hBackO->Fill(det_time_start[i]);
hDetz->Fill(det_z[i]);
hDetzO->Fill(det_z[i]);
hEventID->Fill(eventID);
break; //fill only once
}
}
}
// Hist in Forw O detector (602)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i] == 602 && save_particleID[i]==22) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
hForwO->Fill(det_time_start[i]);
hDetz->Fill(det_z[i]);
hDetzO->Fill(det_z[i]);
hEventID->Fill(eventID);
break; //fill only once
}
}
}
// Hist in Top O detector (603)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i] == 603 && save_particleID[i]==22) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
hTopO->Fill(det_time_start[i]);
hDetz->Fill(det_z[i]);
hDetzO->Fill(det_z[i]);
hEventID->Fill(eventID);
break; //fill only once
}
}
}
// Hist in Down O detector (604)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i] == 604 && save_particleID[i]==22) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
hDownO->Fill(det_time_start[i]);
hDetz->Fill(det_z[i]);
hDetzO->Fill(det_z[i]);
hEventID->Fill(eventID);
break; //fill only once
}
}
}
}
// Calculate Asymmetry
// Temporary F and B histograms
TH1F *hFtemp = (TH1F*) hForwI->Clone(); hFtemp->SetName("hFtemp");
TH1F *hBtemp = (TH1F*) hBackI->Clone(); hBtemp->SetName("hBtemp");
// Sum and difference F and B
TH1F *hSumFB = (TH1F*) hForwI->Clone(); hSumFB->SetName("hSumFB");
hSumFB->Sumw2();
TH1F *hDifFB = (TH1F*) hForwI->Clone(); hDifFB->SetName("hDifFB");
hDifFB->Sumw2();
// Asymmetry histograms!
TH1F *hAsyFB = (TH1F*) hForwI->Clone(); hAsyFB->SetName("hAsyFB");
hAsyFB->SetTitle("Muon decay asymmetry F-B; Time (#mus); Asymmetry");
hAsyFB->Sumw2();
// Calculate difference and sum, then divide
hDifFB->Add(hFtemp,hBtemp,1.,-1.);
hSumFB->Add(hFtemp,hBtemp,1., 1.);
hAsyFB->Divide(hDifFB,hSumFB,1.,1.);
if (FigFlag) {
TCanvas* c1= new TCanvas("c1","canvas 1");
c1->Divide(3,2);
c1->cd(1);
hForwI->Draw();
hForwI->SetLineColor(kBlack);
hBackI->Draw("sames");
hBackI->SetLineColor(kRed);
// hBackI->SetX1NDC(10);
// hBackI->SetX2NDC(20);
hTopI->Draw("sames");
hTopI->SetLineColor(kGreen);
hDownI->Draw("sames");
hDownI->SetLineColor(kBlue);
gStyle->SetOptStat("ne");
c1->cd(2);
hForwO->Draw();
hForwO->SetLineColor(kBlack);
hBackO->Draw("sames");
hBackO->SetLineColor(kRed);
hTopO->Draw("sames");
hTopO->SetLineColor(kGreen);
hDownO->Draw("sames");
hDownO->SetLineColor(kBlue);
gStyle->SetOptStat("ne");
c1->cd(3);
hAsyFB->Draw();
hAsyFB -> Fit("pol0","Q","",0.6, 13.);
gStyle->SetOptStat(1001111);
gStyle->SetOptFit(0001);
gStyle->SetLineColor(2);
c1->cd(4);
hEdeposited->Draw();
hEdepositMuO->Draw("sames");
hEdepositMuO->SetLineColor(kRed);
hEdepositMuI->Draw("sames");
hEdeposited->SetLineColor(kGreen);
//gStyle->SetOptStat("nemr");
c1->cd(5);
hDetz->Draw();
hDetzI->Draw("sames");
hDetzI->SetLineColor(kRed);
hDetzO->Draw("sames");
hDetzO->SetLineColor(kGreen);
} else {
hAsyFB -> Fit("pol0","NQ","",0.6, 13.);
}
TF1 *pol0;
pol0 = (TF1*)gROOT->GetListOfFunctions()->FindObject("pol0");
Double_t chi2=pol0->GetChisquare();
Double_t p1=pol0->GetParameter(0);
Double_t e1=pol0->GetParError(0);
// printf("Chi=%g\tP1=%g +/- %g\n",chi2,p1,e1);
printf("%g\t%g\t%g\n",p1,e1,chi2);
}
void SiSpect::CoinIO( Bool_t FigFlag, Double_t eCut )
{
// Double_t eCut = 0.7; //minimum energy deposition in scintillators [MeV]
Bool_t tofFlag = 0;
if (fChain == 0) return;
//DEFINE HISTOGRAMS
TH1F* hEdeposited = new TH1F("hEdeposited","Energy spectrum; E [MeV]", 250,0.,0.5);
TH1F* hEdepositMu = new TH1F("hEdepositMuI","Energy spectrum Mu/Pos; E [MeV]", 250,0.,0.5);
TH1F* hEdepositCF = new TH1F("hEdepositCF","Energy spectrum; E [keV]", 500,0.,20.0);
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* hTof = new TH1F("hTof", "time-of-flight (#mus)", 1000, 0., 1.);
TH2F* hBeamSpot = new TH2F("hBeamSpot", " x,y", 40, -40., 40., 40, -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);
TH2F* hMuIxy = new TH2F("hMuIxy", " x,y", 40, -40., 40., 40, -40., 40.);
TH2F* hMuOxy = new TH2F("hMuOxy", " x,y", 40, -40., 40., 40, -40., 40.);
TH2F* hMuIOxy = new TH2F("hMuIOxy", " x,y", 40, -40., 40., 40, -40., 40.);
TH2F* hMuSamxy = new TH2F("hMuIOxy", " x,y", 40, -40., 40., 40, -40., 40.);
// Back histogram, i.e. all counts in coincedence 501-601
TH1F* hBack = new TH1F("hBack","Back (#mus)",130,0.,13.);
// Forward histogram, i.e. all counts in coincedence 502-602
TH1F* hForw = new TH1F("hForw","Forw (#mus)",130,0.,13.);
// Top histogram, i.e. all counts in coincedence 503-603
TH1F* hTop = new TH1F("hTop","Top (#mus)",130,0.,13.);
// Down histogram, i.e. all counts in coincedence 504-604
TH1F* hDown = new TH1F("hDown","Down (#mus)",130,0.,13.);
// Mu counter, i.e. all counts in coincedence 501-601
TH1F* hMu = new TH1F("hMu","Muons (#mus)",130,0.,13.);
Double_t xo,yo,zo,xi,yi,zi,xs,ys,zs;
// Sample position
zs = 0.0;
hEdeposited->Sumw2();
hEdepositCF->Sumw2();
hEdepTrig->Sumw2();
Long64_t nentries = fChain->GetEntriesFast();
/* nentries=50000;*/
Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
tofFlag = 0;
/* Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;*/
nb = fChain->GetEntry(jentry); nbytes += nb;
// FILL IN HISTOGRAMS
hEdepositCF->Fill(save_ke[0]/1000.);
hBeamSpot->Fill(save_x[0], save_y[0]);
//for (Int_t i=0; i<det_n; i++)
// { if (det_ID[i]==623) { hEdeposited->Fill(save_ke[i]);}
// }
hTof->Fill(muTargetTime);
if (muTargetTime>0.23) tofFlag = 1;
//tofFlag = 1.;
// Hist in Muon detector (501) coincidence (601)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==501 && save_particleID[i]==-13) {
// Mu outer hit
xo=det_x[i];
yo=det_y[i];
zo=det_z[i];
hMuOxy->Fill(xo, yo);
// printf("Outer: %6.2f\t%6.2f\n",xo,yo);
hEdepositMu->Fill(det_edep[i]);
for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==601 && save_particleID[j]==-13) {
printf("Inner %d\t%d\n",det_ID[j] ,save_particleID[j]);
//printf("Inner\n");
// Mu inner hit
xi=det_x[j];
yi=det_y[j];
zi=det_z[j];
hMuIxy->Fill(xi, yi);
hEdepositMu->Fill(det_edep[j]);
hMu->Fill(det_time_start[j]);
hDetz->Fill(det_z[j]);
hEventID->Fill(eventID);
// Propagate Mu trajectory on sample
xs=xi+((xi-xo)/(zi-zo))*(zs-zi);
ys=yi+((yi-yo)/(zi-zo))*(zs-zi);
hMuIOxy->Fill(xs-save_x[j], ys-save_y[j]);
hMuSamxy->Fill(save_x[j], save_y[j]);
// printf("%6.2f\t%6.2f\n",xs,ys);
break; //fill only once
}
}
}
//printf("\n");
}
// Hist in Back detector (501) coincidence (601)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==501 && save_particleID[i]==22) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==601 && save_particleID[j]==22) {
hEdeposited->Fill(det_edep[j]);
if (det_edep[j]>eCut){
hBack->Fill(det_time_start[j]);
hDetz->Fill(det_z[j]);
hEventID->Fill(eventID);
}
break; //fill only once
}
}
}
}
}
// Hist in Forw detector (502) coincidence (602)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==502 && save_particleID[i]==22) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==602 && save_particleID[j]==22) {
if (det_edep[j]>eCut){
hForw->Fill(det_time_start[j]);
hDetz->Fill(det_z[j]);
hEventID->Fill(eventID);
}
break; //fill only once
}
}
}
}
}
// Hist in Top detector (503) coincidence (603)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==503 && save_particleID[i]==22) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==603 && save_particleID[j]==22) {
if (det_edep[j]>eCut){
hTop->Fill(det_time_start[j]);
hDetz->Fill(det_z[j]);
hEventID->Fill(eventID);
}
break; //fill only once
}
}
}
}
}
// Hist in Down detector (504) coincidence (604)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==504 && save_particleID[i]==22) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==604 && save_particleID[j]==22) {
if (det_edep[j]>eCut){
hDown->Fill(det_time_start[j]);
hDetz->Fill(det_z[j]);
hEventID->Fill(eventID);
}
break; //fill only once
}
}
}
}
}
}
// Calculate Asymmetry
// Temporary F and B histograms
TH1F *hFtemp = (TH1F*) hForw->Clone(); hFtemp->SetName("hFtemp");
TH1F *hBtemp = (TH1F*) hBack->Clone(); hBtemp->SetName("hBtemp");
// Temporary T and D histograms
TH1F *hTtemp = (TH1F*) hTop->Clone(); hTtemp->SetName("hTtemp");
TH1F *hDtemp = (TH1F*) hDown->Clone(); hDtemp->SetName("hDtemp");
// Sum and difference F and B
TH1F *hSumFB = (TH1F*) hForw->Clone(); hSumFB->SetName("hSumFB");
hSumFB->Sumw2();
TH1F *hDifFB = (TH1F*) hForw->Clone(); hDifFB->SetName("hDifFB");
hDifFB->Sumw2();
// Sum and difference T and D
TH1F *hSumTD = (TH1F*) hTop->Clone(); hSumTD->SetName("hSumTD");
hSumTD->Sumw2();
TH1F *hDifTD = (TH1F*) hDown->Clone(); hDifTD->SetName("hDifTD");
hDifTD->Sumw2();
// Asymmetry histograms!
TH1F *hAsyFB = (TH1F*) hForw->Clone(); hAsyFB->SetName("hAsyFB");
hAsyFB->SetTitle("Muon decay asymmetry F-B; Time (#mus); Asymmetry");
hAsyFB->Sumw2();
TH1F *hAsyTD = (TH1F*) hTop->Clone(); hAsyTD->SetName("hAsyTD");
hAsyTD->SetTitle("Muon decay asymmetry T-D; Time (#mus); Asymmetry");
hAsyTD->Sumw2();
// Calculate difference and sum, then divide
hDifFB->Add(hFtemp,hBtemp,1.,-1.);
hSumFB->Add(hFtemp,hBtemp,1., 1.);
hAsyFB->Divide(hDifFB,hSumFB,1.,1.);
hDifTD->Add(hTtemp,hDtemp,1.,-1.);
hSumTD->Add(hTtemp,hDtemp,1., 1.);
hAsyTD->Divide(hDifTD,hSumTD,1.,1.);
if (FigFlag) {
TCanvas* c1= new TCanvas("c1","canvas 1");
c1->Divide(3,3);
c1->cd(1);
hBack->Draw();
hBack->SetLineColor(kBlack);
hForw->Draw("sames");
hForw->SetLineColor(kRed);
hTop->Draw("sames");
hTop->SetLineColor(kGreen);
hDown->Draw("sames");
hDown->SetLineColor(kBlue);
c1->cd(2);
hAsyFB->Draw();
hAsyFB -> Fit("pol0","Q","",0.1, 13.);
gStyle->SetOptStat(1001111);
gStyle->SetOptFit(0001);
gStyle->SetLineColor(2);
hAsyTD->Draw("sames");
hAsyTD -> Fit("pol0","Q","",0.1, 13.);
gStyle->SetOptStat(1001111);
gStyle->SetOptFit(0001);
gStyle->SetLineColor(2);
c1->cd(3);
hDetz->Draw();
c1->cd(4);
hMu->Draw();
c1->cd(5);
hEdeposited->Draw();
hEdepositMu->Draw("sames");
hEdepositMu->SetLineColor(kRed);
/* det_edep_pos->Draw("sames");
det_edep_pos->SetLineColor(kGreen);
det_edep_mup->Draw("sames");
det_edep_mup->SetLineColor(kPink);*/
c1->cd(6);
hMuIOxy->Draw();
hMuIOxy->Draw("cont0 same");
c1->cd(7);
hMuSamxy->Draw();
hMuSamxy->Draw("cont0 same");
} else {
hAsyFB -> Fit("pol0","NQ","",0.1, 13.);
}
TF1 *pol0;
pol0 = (TF1*)gROOT->GetListOfFunctions()->FindObject("pol0");
Double_t chi2=pol0->GetChisquare();
Double_t p1=pol0->GetParameter(0);
Double_t e1=pol0->GetParError(0);
Double_t NDet=hDetz->GetSum();
// printf("Chi=%g\tP1=%g +/- %g\n",chi2,p1,e1);
//printf("%g\t%g\t%g\t%g\n",p1,e1,chi2,NDet);
Double_t NhF = hBack->GetEntries();
Double_t NhB = hForw->GetEntries();
Double_t asymFB = (NhB-NhF)/(NhB+NhF);
Double_t NhT = hTop->GetEntries();
Double_t NhD = hDown->GetEntries();
Double_t asymTD = (NhT-NhD)/(NhT+NhD);
cout << "FB Asymmetry = " << asymFB << endl;
cout << "TD Asymmetry = " << asymTD << endl;
}
/*
ax=(x2-x1)/(z2-z1)
ay=(y2-y1)/(z2-z1)
x3=x2+ax*(z3-z2)
y3=y2+ay*(x3-z2)
*/
Reference in New Issue View Git Blame Copy Permalink