musrsim/run/macros/HiFiSim.C
2025-03-27 22:05:05 +01:00

254 lines
7.1 KiB
C

#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <string>
void EnvatZ(Int_t fnumber, char* fname, TCanvas* c1)
{
// For old lem4 simulation files
TFile* f1 = new TFile(fname);
//DEFINE HISTOGRAMS
TH3F* hBeamSpot = new TH3F("hBeamSpot", " x,y,z", 32, -100., 100., 32, -100., 100.,2900,-900,2000);
// TH1F* hPosZ = new TH1F("hPosZ", "z",1600,-800.,800.);
// hBeamSpot->Fill(muTargetPosX[0], muTargetPosY[0]);
// hPosZ->Fill(muTargetPosZ[0]);
// TCanvas* c1= new TCanvas("c1","canvas 1");
// hPosZ->Draw();
// gStyle->SetPalette(1,0);
// TCanvas* c1= new TCanvas("c1","canvas 1");
// c1->cd(1);
t1->Draw("muTargetPosZ[0]:muTargetPosY[0]:muTargetPosX[0]>>hBeamSpot");
// hBeamSpot->Draw("lego");
// hBeamSpot->Draw();
hBeamSpot->DrawCopy("cont0 same");
Double_t Xrms=hBeamSpot->GetRMS(1);
Double_t Yrms=hBeamSpot->GetRMS(2);
Double_t Zrms=hBeamSpot->GetRMS(3);
Double_t Xavg=hBeamSpot->GetMean(1);
Double_t Yavg=hBeamSpot->GetMean(2);
Double_t Zavg=hBeamSpot->GetMean(3);
Double_t Nbins=hBeamSpot->Integral();
// Long64_t Nbins=t1.save_x->GetEntries();
if (Xrms != 0) {
printf("%d\t%g\t%g\t%g\t%g\t%g\t%g\n",fnumber,Zavg,Xrms,Yrms,Xavg,Yavg,Nbins);
}
// Clean up
// delete c1;
delete f1;
// hPosZ->Draw();
// sprintf(filenamePrint,"musrSrClass_%s.ps",fileNr);
// c1->Print(filenamePrint);
}
void BeamSlice(char* fname,Int_t SaveID)
{
Double_t zSlice,Xrms,Yrms,Xavg,Yavg,Nbins;
// For new musrSim simulation files
TFile* f1 = new TFile(fname);
TString Cond;
// Define an array with save volume IDs
Cond.Form("save_detID== %04d&&save_particleID==-13",SaveID);
t1->Draw("save_z:save_y:save_x>>bspot_hist(32,-100,100,32,-100,100,100,-2600,1300)",Cond);
// t1->Draw("save_y:save_x>>bspot_hist(64,-40,40,64,-40,40)",Cond);
gStyle->SetPalette(1,0);
bspot_hist->SetTitle("Beam Cross Section [mm]");
bspot_hist->Draw();
bspot_hist->Draw("cont0 same");
zSlice=bspot_hist->GetMean(3);
Xrms=bspot_hist->GetRMS(1);
Yrms=bspot_hist->GetRMS(2);
Xavg=bspot_hist->GetMean(1);
Yavg=bspot_hist->GetMean(2);
Nbins=bspot_hist->Integral();
printf("%d\t%g\t%g\t%g\t%g\t%g\t%g\n",SaveID,zSlice,Xrms,Yrms,Xavg,Yavg,Nbins);
f1->Close();
}
void Slice(char* fname,Int_t SaveID)
{
// For new musrSim simulation files
TFile* f1 = new TFile(fname);
TString Cond;
// Define an array with save volume IDs
Cond.Form("save_detID== %04d&&save_particleID==-13",SaveID);
gStyle->SetPalette(1,0);
TCanvas* c1 = new TCanvas("c1","Slice Properties",60,40,600,800);
c1->Divide(2,3);
// Initial beam spot
c1->cd(1);
t1->Draw("muIniPosY[0]:muIniPosX[0]>>ini_hist(32,-100,100,32,-100,100)");
ini_hist->SetTitle("Muon initial beam cross section [mm]");
ini_hist->Draw();
ini_hist->Draw("cont0 same");
ini_hist->ProjectionX("xisect",31,33);
ini_hist->ProjectionY("yisect",31,33);
// Beam spot at slice
c1->cd(2);
t1->Draw("save_y:save_x>>bspot_hist(32,-100,100,32,-100,100)",Cond);
// t1->Draw("save_y:save_x>>bspot_hist(64,-40,40,64,-40,40)",Cond);
bspot_hist->SetTitle("Beam Cross Section [mm]");
bspot_hist->Draw();
bspot_hist->Draw("cont0 same");
// Initial x polarization
c1->cd(3);
t1->Draw("muIniPolX >> ipol_x(128, -1.1, 1.1)"); //0, 2
ipol_x->SetTitle("Polarization X");
ipol_x->GetXaxis()->SetNdivisions(405);
ipol_x->GetYaxis()->SetNdivisions(406);
ipol_x->GetXaxis()->SetTickLength(0.018);
ipol_x->GetYaxis()->SetTickLength(0.018);
ipol_x->Draw("F");
ipol_x->SetFillStyle(1001);
ipol_x->SetFillColor(kGreen-5);
ipol_x->Draw();
// x polarization at slice
c1->cd(4);
t1->Draw("save_polx >> pol_x(128, -1.1, 1.1)",Cond); //0, 2
pol_x->SetTitle("Polarization X");
pol_x->GetXaxis()->SetNdivisions(405);
pol_x->GetYaxis()->SetNdivisions(406);
pol_x->GetXaxis()->SetTickLength(0.018);
pol_x->GetYaxis()->SetTickLength(0.018);
pol_x->Draw("F");
pol_x->SetFillStyle(1001);
pol_x->SetFillColor(kGreen-5);
pol_x->Draw();
//c1->SaveAs("data/sr_bspot_0.eps");
// Initial z polarization
c1->cd(5);
t1->Draw("muIniPolZ >> ipol_z(128, -1.1, 1.1)"); //0, 2
ipol_z->SetTitle("Polarization Z");
ipol_z->GetXaxis()->SetNdivisions(405);
ipol_z->GetYaxis()->SetNdivisions(406);
ipol_z->GetXaxis()->SetTickLength(0.018);
ipol_z->GetYaxis()->SetTickLength(0.018);
ipol_z->Draw("F");
ipol_z->SetFillStyle(1001);
ipol_z->SetFillColor(kGreen-5);
ipol_z->Draw();
// z polarization at slice
c1->cd(6);
t1->Draw("save_polz >> pol_z(128, -1.1, 1.1)",Cond); //-1., 1
pol_z->SetTitle("Polarization Z");
pol_z->GetXaxis()->SetNdivisions(405);
pol_z->GetYaxis()->SetNdivisions(406);
pol_z->GetXaxis()->SetTickLength(0.018);
pol_z->GetYaxis()->SetTickLength(0.018);
pol_z->Draw("F");
pol_z->SetFillStyle(1001);
pol_z->SetFillColor(38);
pol_z->Draw();
// f1->Close();
}
void CompSlice(Int_t ifnamei, Int_t ifnamef, Int_t SaveID)
{
TString fname;
TString title;
TString Cond;
Double_t B,Xrms,Yrms,Xavg,Yavg,Nbins;
Cond.Form("save_detID== %04d&&save_particleID==-13&&save_n==1",SaveID);
gStyle->SetPalette(1,0);
TCanvas* c1 = new TCanvas("c1","Slice Properties",60,40,600,800);
c1->Divide(4,4);
Int_t ifname;
Int_t counter=0;
printf("B \t Xrms \t Yrms \t Xavg \t Yavg \t N\n");
for (ifname=ifnamei; ifname<=ifnamef; ifname++) {
counter++;
fname.Form("data/musr_%d.root",ifname);
TFile* f1 = new TFile(fname);
// Initial beam spot
c1->cd(counter);
t1->Draw("save_y:save_x>>bspot_hist(32,-100,100,32,-100,100)",Cond);
B=(ifname-ifnamei)*0.05+0.6;
title.Form("B=%f Tesla",B);
bspot_hist->SetTitle(title);
bspot_hist->Draw();
bspot_hist->Draw("cont0 same");
Xrms=bspot_hist->GetRMS(1);
Yrms=bspot_hist->GetRMS(2);
Xavg=bspot_hist->GetMean(1);
Yavg=bspot_hist->GetMean(2);
Nbins=bspot_hist->Integral();
printf("%g\t%g\t%g\t%g\t%g\t%g\n",B,Xrms,Yrms,Xavg,Yavg,Nbins);
// f1->Close();
}
}
void BeamSlices(char* fname,Int_t SaveIDi, Int_t SaveIDf)
{
// For new musrSim simulation files
TString Cond;
Int_t i;
printf("SliceID \t Z \t Xrms \t Yrms \t Xavg \t Yavg \t N\n");
for (i=SaveIDi; i<=SaveIDf; i++) {
BeamSlice(fname,i);
}
}
void NatSlice(char* fname,Int_t SaveID)
{
// For new musrSim simulation files
TString Cond;
Double_t zSlice,Xrms,Yrms,Xavg,Yavg,Nbins;
Int_t i;
TFile* f1 = new TFile(fname);
TString Cond;
// Define an array with save volume IDs
Cond.Form("save_detID== %04d&&save_particleID==-13",SaveID);
t1->Draw("save_z:save_y:save_x>>bspot_hist(32,-100,100,32,-100,100,100,-2600,1300)",Cond);
// t1->Draw("save_y:save_x>>bspot_hist(64,-40,40,64,-40,40)",Cond);
gStyle->SetPalette(1,0);
bspot_hist->SetTitle("Beam Cross Section [mm]");
bspot_hist->Draw();
bspot_hist->Draw("cont0 same");
zSlice=bspot_hist->GetMean(3);
Xrms=bspot_hist->GetRMS(1);
Yrms=bspot_hist->GetRMS(2);
Xavg=bspot_hist->GetMean(1);
Yavg=bspot_hist->GetMean(2);
Nbins=bspot_hist->Integral();
// printf("%d\t%g\t%g\t%g\t%g\t%g\t%g\n",i,zSlice,Xrms,Yrms,Xavg,Yavg,Nbins);
printf("%d\t",Nbins);
f1->Close();
}