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

Added root macros

Browse Source
This commit is contained in:
nemu 2014-08-15 17:18:30 +02:00
parent f60f4874d6
commit c6ebc44ccc
25 changed files with 5433 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1021
run/macros/Beam.old.C Normal file
View File

File diff suppressed because it is too large Load Diff

227
run/macros/BeamEnv.C Normal file
View File

@ -0,0 +1,227 @@
#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", 64, -40., 40., 64, -40., 40.,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,Zrms,Xavg,Yavg,Zavg,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",Cond);
// t1->Draw("save_z:save_y:save_x>>bspot_hist(64,-40,40,100,-600,6000,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);
Zrms=bspot_hist->GetRMS(3);
Xavg=bspot_hist->GetMean(1);
Yavg=bspot_hist->GetMean(2);
Zavg=bspot_hist->GetMean(3);
Nbins=bspot_hist->Integral();
printf("%d\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n",SaveID,zSlice,Xrms,Yrms,Zrms,Xavg,Yavg,Zavg,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]:muIniPosZ[0]>>ini_hist");
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(64,-40,40,64,-40,40)",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 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 Zrms \t Xavg \t Yavg \t Zavg \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(64,-40,40,64,-40,40,100,-2600,300)",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();
}
void ToF(char* fname,Int_t SaveID)
{
// For new musrSim simulation files
TFile* f1 = new TFile(fname);
TString Cond;
Cond.Form("save_detID== %03d&&save_particleID==-13",SaveID);
t1->Draw("save_time>>tof",Cond);
// t1->Draw("save_y:save_x>>bspot_hist(64,-40,40,64,-40,40)",Cond);
// f1->Close();
}

24
run/macros/BeamEnvLoop.C Normal file
View File

@ -0,0 +1,24 @@
// To run: root -b LoopFiles.C
{
TString fileDir,fileName,file;
Int_t fileNumber;
gStyle->SetPalette(1,0);
TCanvas* c1= new TCanvas("c1","canvas 1");
gROOT->ProcessLine(".L BeamEnv.C");
// gROOT->ProcessLine(".L sr_tr.C");
fileDir = "./data/";
printf("Run \t Z \t Xrms \t Yrms \t Xavg \t Yavg \t N\n");
for (fileNumber=52421; fileNumber<=52424; fileNumber++)
{
file.Form("musr_%04d.root",fileNumber);
fileName= fileDir + file;
EnvatZ(fileNumber,fileName,c1);
}
}

21
run/macros/Coincedence.C Normal file
View File

@ -0,0 +1,21 @@
{
TH1F* hEdeposited = new TH1F("hEdeposited","Energy spectrum; E [MeV]", 250,0.,2.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* hDettUp = new TH1F("hDettUp", "Muon arrival times Up (#mus)", 130,0.,13.);
TH1F* hDettDown = new TH1F("hDettDown", "Muon arrival times Down (#mus)", 130,0.,13.);
TH1F* hDettRight = new TH1F("hDettRight", "Muon arrival times Right (#mus)", 130,0.,13.);
TH1F* hDettLeft = new TH1F("hDettLeft", "Muon arrival times Left (#mus)", 130,0.,13.);
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]", 140, -140., 140.);
TH1F* hEventID = new TH1F("hEventID", "Event ID", 10001, -0.5, 10000.5);
for (Int_t i=0; i<100; i++) {
hDettLeft->Fill(det_time_start[j]);
hDetz->Fill(det_z[i]);
hEventID->Fill(eventID);
}
}

6
run/macros/Commands.C Normal file
View File

@ -0,0 +1,6 @@
{
TFile* f1 = new TFile("lem4_5000.root");
gROOT->ProcessLine(".L NewAnalysis.C");
NewSpec t ;
t.Coincidence(1,2);
}

39
run/macros/ExportAscii.C Normal file
View File

@ -0,0 +1,39 @@
/**
* \brief Export Single Histogram into ASCII file
*/
Bool_t SingleExportAscii(TH1* hist, TString &filename, TString folder="", TString separator="\t")
{
Int_t i,j;
Double_t xcenter, xwidth;
Bool_t success=kFALSE;
filename = folder + hist->GetName() + ".dat";
ofstream file_out(filename);
file_out << "# Output " << hist->ClassName() << ": " << hist->GetName() << " (" << hist->GetTitle() << ")\n";
if (hist->GetDimension()==1)
{
file_out << "# BinCenter" << separator << "Content" << separator << "BinHalfWidth" << separator << "Error\n";
for (i=1; i<=hist->GetNbinsX(); i++)
file_out << hist->GetBinCenter(i) << separator << hist->GetBinContent(i) << separator << hist->GetBinWidth(i)/2 << separator << hist->GetBinError(i) << endl;
if (i>1)
success=kTRUE;
}
else if (hist->GetDimension()==2)
{
file_out << "# xBinCenter" << separator << "yBinCenter" << separator << "Content" << separator << "xBinHalfWidth" << separator << "yBinHalfWidth" << separator << "Error" << endl;
for (i=1; i <= hist->GetNbinsX(); i++)
{
xcenter = hist->GetXaxis()->GetBinCenter(i);
xwidth = hist->GetXaxis()->GetBinWidth(i)/2;
for (j=1; j <= hist->GetNbinsY(); j++)
file_out << xcenter << separator << hist->GetYaxis()->GetBinCenter(j) << separator << hist->GetBinContent(i,j) << separator << xwidth << separator << hist->GetYaxis()->GetBinWidth(j)/2 << separator << hist->GetBinError(i,j) << endl;
if (j>1)
file_out << endl; // produce a blank line after each set of Ybins for a certain Xbin. Gnuplot likes this.
}
if (i>1)
success=kTRUE;
}
file_out.close();
if (success == kTRUE)
cout << "*** TRemiHistExport: Histogram " << hist->GetName() << " written to " << filename << endl;
return success;
}

253
run/macros/HiFiSim.C Normal file
View File

@ -0,0 +1,253 @@
#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();
}

24
run/macros/LoopFiles.C Normal file
View File

@ -0,0 +1,24 @@
// To run: root -b LoopFiles.C
{
TString fileDir,fileName,file;
Int_t fileNumber,L1,L3;
gROOT->ProcessLine(".L BeamEnv.C");
fileDir = "data/";
//printf("Run \t Asy \t dAsy \t Chi2 \t N\n");
fileNumber=0;
for (L3=0; L3<=16 ; L3++) {
for (L1=0; L1<=16 ; L1++) {
fileNumber++;
file.Form("musr_142%02d.root",fileNumber);
fileName= fileDir + file;
// printf("File is : %s\n", (char*) fileName);
NatSlice(fileName,901);
}
printf("\n");
}
}

24
run/macros/LoopFiles1.C Normal file
View File

@ -0,0 +1,24 @@
// To run: root -b LoopFiles.C
{
TString fileDir,fileName,file;
Int_t fileNumber,L1,L3;
gROOT->ProcessLine(".L BeamEnv.C");
fileDir = "data/";
//printf("Run \t Asy \t dAsy \t Chi2 \t N\n");
fileNumber=0;
for (L3=0; L3<=16 ; L3++) {
for (L1=0; L1<=16 ; L1++) {
fileNumber++;
file.Form("musr_142%02d.root",fileNumber);
fileName= fileDir + file;
// printf("File is : %s\n", (char*) fileName);
NatSlice(fileName,901);
}
printf("\n");
}
}

24
run/macros/LoopFilesOne.C Normal file
View File

@ -0,0 +1,24 @@
// To run: root -b LoopFiles.C
{
TString fileDir,fileName,file;
Int_t fileNumber,CR;
gROOT->ProcessLine(".L BeamEnv.C");
fileDir = "data/";
//printf("Run \t Asy \t dAsy \t Chi2 \t N\n");
fileNumber=0;
for (CR=0; CR<=16 ; CR++) {
fileNumber++;
file.Form("musr_101%02d.root",fileNumber);
fileName= fileDir + file;
// printf("File is : %s\n", (char*) fileName);
// NatSlice(fileName,901);
printf("%d\t",CR);
BeamSlice(fileName,901);
// printf("\n");
}
}

26
run/macros/LoopFilesThree.C Normal file
View File

@ -0,0 +1,26 @@
// To run: root -b LoopFiles.C
{
TString fileDir,fileName,file;
Int_t fileNumber,L1,L3;
gROOT->ProcessLine(".L BeamEnv.C");
fileDir = "data/";
//printf("Run \t Asy \t dAsy \t Chi2 \t N\n");
fileNumber=0;
for (L3=0; L3<=10 ; L3++) {
for (L1=0; L1<=10 ; L1++) {
for (CR=0; CR<=16; CR++) {
fileNumber++;
file.Form("musr_115%02d.root",fileNumber);
fileName= fileDir + file;
// printf("File is : %s\n", (char*) fileName);
NatSlice(fileName,901);
}
printf("\n");
}
}
}

169
run/macros/NewAnalysis.C Normal file
View File

@ -0,0 +1,169 @@
#define analysis_cxx
#include "analysis.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
void analysis::Loop()
{
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.,2.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* hDettUp = new TH1F("hDettUp", "Muon arrival times Up (#mus)", 130,0.,13.);
TH1F* hDettDown = new TH1F("hDettDown", "Muon arrival times Down (#mus)", 130,0.,13.);
TH1F* hDettRight = new TH1F("hDettRight", "Muon arrival times Right (#mus)", 130,0.,13.);
TH1F* hDettLeft = new TH1F("hDettLeft", "Muon arrival times Left (#mus)", 130,0.,13.);
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]", 140, -140., 140.);
TH1F* hEventID = new TH1F("hEventID", "Event ID", 10001, -0.5, 10000.5);
hEdeposited->Sumw2();
hEdepositCF->Sumw2();
hEdepTrig->Sumw2();
Long64_t nentries = fChain->GetEntriesFast();
//nentries=1000;
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.;
// Coincidence between 1 and 33 scintillators
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==1) {
hEdeposited->Fill(det_edep[i]);
for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==36) {
printf("Found coincedence\n");
if (det_edep[i]>eCut && det_edep[j]>eCut && tofFlag){
hDettUp->Fill(det_time_start[j]);
hDetz->Fill(det_z[i]);
hEventID->Fill(eventID);
break; //fill only once
}
}
}
}
}
// Coincidence between 17 and 49 scintillators
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==17) {
hEdeposited->Fill(det_edep[i]);
for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==49){
if (det_edep[i]>eCut && det_edep[j]>eCut && tofFlag){
hDettDown->Fill(det_time_start[j]);
hDetz->Fill(det_z[i]);
hEventID->Fill(eventID);
break;
}
}
}
}
}
// Coincidence between 9 and 42 scintillators
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==9) {
hEdeposited->Fill(det_edep[i]);
for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==42) {
if (det_edep[i]>eCut && det_edep[j]>eCut && tofFlag){
hDettRight->Fill(det_time_start[j]);
hDetz->Fill(det_z[i]);
hEventID->Fill(eventID);
break;
}
}
}
}
}
// Coincidence between 25 and 57 scintillators
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==25) {
hEdeposited->Fill(det_edep[i]);
for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==57){
if (det_edep[i]>eCut && det_edep[j]>eCut && tofFlag){
hDettLeft->Fill(det_time_start[j]);
hDetz->Fill(det_z[i]);
hEventID->Fill(eventID);
break;
}
}
}
}
}
}
// PLOT HISTOGRAMS
TCanvas* c1= new TCanvas("c1","canvas 1");
// hEdeposited->Scale(1/(number_of_generated_events*hEdeposited->GetBinWidth(1)));
TH1F *htemp1 = (TH1F*) hDettUp->Clone(); htemp1->SetName("htemp1");
TH1F *htemp2 = (TH1F*) hDettUp->Clone(); htemp2->SetName("htemp2");
TH1F *htemp3 = (TH1F*) hDettUp->Clone(); htemp3->SetName("htemp3");
TH1F *htemp4 = (TH1F*) hDettUp->Clone(); htemp4->SetName("htemp4");
TH1F *htemp5 = (TH1F*) hDettUp->Clone(); htemp5->SetName("htemp5");
htemp5->SetTitle("Muon decay asymmetry U-D; Time (#mus); Asymmetry");
htemp5->Sumw2();
TH1F *htemp6 = (TH1F*) hDettUp->Clone(); htemp6->SetName("htemp6");
htemp6->SetTitle("Muon decay asymmetry L-R; Time (#mus); Asymmetry");
htemp6->Sumw2();
htemp1->Add(hDettUp,hDettDown,1.,-1.);
htemp2->Add(hDettUp,hDettDown,1., 1.);
htemp5->Divide(htemp1,htemp2,1.,1.);
htemp3->Add(hDettLeft,hDettRight,1.,-1.);
htemp4->Add(hDettLeft,hDettRight,1., 1.);
htemp6->Divide(htemp3,htemp4,1.,1.);
htemp5->SetLineWidth(2);
htemp5->SetLineColor(2);
htemp5->SetAxisRange(-0.8,0.8,"Y");
htemp5->Draw("hist");
htemp6->SetLineWidth(2);
htemp6->SetLineColor(4);
htemp6->SetAxisRange(-0.8,0.8,"Y");
htemp6->Draw("same");
hEdepositCF->SetLineWidth(2);
hEdepositCF->SetLineColor(2);
hEdepositCF->Draw("hist");
//hEdepositCF->Fit("gaus");
TCanvas* c2= new TCanvas("c2","canvas 2");
hEdeposited->Draw();
// sprintf(filenamePrint,"musrSrClass_%s.ps",fileNr);
// c1->Print(filenamePrint);
}

1292
run/macros/NewSpec.C Normal file
View File

File diff suppressed because it is too large Load Diff

648
run/macros/NewSpeccp.C Normal file
View File

@ -0,0 +1,648 @@
#define NewSpec_cxx
#include "NewSpec.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>
void NewSpec::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.,2.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]", 140, -140., 140.);
TH1F* hEventID = new TH1F("hEventID", "Event ID", 10001, -0.5, 10000.5);
// Back inner histogram, i.e. all counts in segments 1-16
TH1F* hBackI = new TH1F("hBackI","Muon arrival times Back I (#mus)",1300,0.,13.);
// Forward inner histogram, i.e. all counts in segments 17-32
TH1F* hForwI = new TH1F("hForwI","Muon arrival times Forw I (#mus)",1300,0.,13.);
// Back outer histogram, i.e. all counts in segments 33-48
TH1F* hBackO = new TH1F("hBackO","Muon arrival times Back O (#mus)",1300,0.,13.);
// Forward outer histogram, i.e. all counts in segments 49-64
TH1F* hForwO = new TH1F("hForwO","Muon arrival times Forw O (#mus)",1300,0.,13.);
hEdeposited->Sumw2();
hEdepositCF->Sumw2();
hEdepTrig->Sumw2();
Long64_t nentries = fChain->GetEntriesFast();
//nentries=1000;
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 Back I detector (1-16)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]<=16 && det_ID[i] >=1) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
hBackI->Fill(det_time_start[i]);
hDetz->Fill(det_z[i]);
hEventID->Fill(eventID);
break; //fill only once
}
}
}
// Hist in Forw I detector (17-32)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]>=17 && det_ID[i]<=32) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
hForwI->Fill(det_time_start[i]);
hDetz->Fill(det_z[i]);
hEventID->Fill(eventID);
break; //fill only once
}
}
}
// Hist in Back O detector (33-48)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]<=48 && det_ID[i] >=33) {
// hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
hBackO->Fill(det_time_start[i]);
hDetz->Fill(det_z[i]);
hEventID->Fill(eventID);
break; //fill only once
}
}
}
// Hist in Forw O detector (49-64)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]>=49 && det_ID[i]<=64) {
// hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
hForwO->Fill(det_time_start[i]);
hDetz->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(2,3);
c1->cd(1);
hBackI->Draw();
c1->cd(2);
hForwI->Draw();
c1->cd(3);
hBackO->Draw();
c1->cd(4);
hForwO->Draw();
c1->cd(5);
hAsyFB->Draw();
hAsyFB -> Fit("pol0","Q");
gStyle->SetOptStat(1001111);
gStyle->SetOptFit(0001);
gStyle->SetLineColor(2);
c1->cd(6);
hDetz->Draw();
} else {
hAsyFB -> Fit("pol0","NQ");
}
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 NewSpec::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.,2.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]", 140, -140., 140.);
TH1F* hEventID = new TH1F("hEventID", "Event ID", 10001, -0.5, 10000.5);
// Back inner histogram, i.e. all counts in segments 1-16
TH1F* hBack = new TH1F("hBack","Muon arrival times Back (#mus)",1300,0.,13.);
// Forward inner histogram, i.e. all counts in segments 17-32
TH1F* hForw = new TH1F("hForw","Muon arrival times Forw (#mus)",1300,0.,13.);
hEdeposited->Sumw2();
hEdepositCF->Sumw2();
hEdepTrig->Sumw2();
Long64_t nentries = fChain->GetEntriesFast();
//nentries=1000;
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 Back detector (1-16) coincidence (33-48)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]<=16 && det_ID[i] >=1) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
for (Int_t j=0; j<det_n; j++){
if (det_ID[j]<=48 && det_ID[j] >=33) {
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 I detector (17-32) coincidence (49-64)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]>=17 && det_ID[i]<=32) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
for (Int_t j=0; j<det_n; j++){
if (det_ID[j]<=64 && det_ID[j] >=49) {
if (det_edep[j]>eCut){
hForw->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");
// 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();
// Asymmetry histograms!
TH1F *hAsyFB = (TH1F*) hForw->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(2,2);
c1->cd(1);
hBack->Draw();
c1->cd(2);
hForw->Draw();
c1->cd(3);
hAsyFB->Draw();
hAsyFB -> Fit("pol0","Q");
gStyle->SetOptStat(1001111);
gStyle->SetOptFit(0001);
gStyle->SetLineColor(2);
c1->cd(4);
hDetz->Draw();
} else {
hAsyFB -> Fit("pol0","NQ");
}
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);
}
void NewSpec::TBLRCoinIO( 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.,2.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]", 140, -140., 140.);
TH1F* hEventID = new TH1F("hEventID", "Event ID", 10001, -0.5, 10000.5);
// Top, Bottom, Left and Right histograms
TH1F* hT = new TH1F("hT","Muon arrival times Top (#mus)",1300,0.,13.);
TH1F* hB = new TH1F("hB","Muon arrival times Bottom (#mus)",1300,0.,13.);
TH1F* hL = new TH1F("hL","Muon arrival times Left (#mus)",1300,0.,13.);
TH1F* hR = new TH1F("hR","Muon arrival times Right (#mus)",1300,0.,13.);
hEdeposited->Sumw2();
hEdepositCF->Sumw2();
hEdepTrig->Sumw2();
Long64_t nentries = fChain->GetEntriesFast();
//nentries=1000;
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 Top detector (11-14)+(27-30) coincidence (43-46)+(59-62)
for (Int_t i=0; i<det_n; i++) {
if ((det_ID[i]<=14 && det_ID[i] >=11) || (det_ID[i]<=30 && det_ID[i] >=27)) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
for (Int_t j=0; j<det_n; j++){
if ((det_ID[j]<=46 && det_ID[j] >=43) || (det_ID[j]<=62 && det_ID[j] >=59)) {
if (det_edep[j]>eCut){
hT->Fill(det_time_start[j]);
hDetz->Fill(det_z[j]);
hEventID->Fill(eventID);
break; //fill only once
}
}
}
}
}
}
// Hist in Bottom detector (3-6)+(19-22) coincidence (35-38)+(51-54)
for (Int_t i=0; i<det_n; i++) {
if ((det_ID[i]<=6 && det_ID[i] >=3) || (det_ID[i]<=22 && det_ID[i] >=19)) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
for (Int_t j=0; j<det_n; j++){
if ((det_ID[j]<=38 && det_ID[j] >=35) || (det_ID[j]<=54 && det_ID[j] >=51)) {
if (det_edep[j]>eCut){
hB->Fill(det_time_start[j]);
hDetz->Fill(det_z[j]);
hEventID->Fill(eventID);
break; //fill only once
}
}
}
}
}
}
// Hist in Left detector (7-10)+(23-26) coincidence (39-42)+(55-58)
for (Int_t i=0; i<det_n; i++) {
if ((det_ID[i]<=10 && det_ID[i] >=7) || (det_ID[i]<=26 && det_ID[i] >=23)) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
for (Int_t j=0; j<det_n; j++){
if ((det_ID[j]<=42 && det_ID[j] >=39) || (det_ID[j]<=58 && det_ID[j] >=55)) {
if (det_edep[j]>eCut){
hL->Fill(det_time_start[j]);
hDetz->Fill(det_z[j]);
hEventID->Fill(eventID);
break; //fill only once
}
}
}
}
}
}
// Hist in Right detector (1,2,15,16)+(17,18,31,32) coincidence (33,34,47,48)+(49,50,63,64)
for (Int_t i=0; i<det_n; i++) {
if ((det_ID[i]<=18 && det_ID[i] >=15) || (det_ID[i]<=2 && det_ID[i] >=1) || (det_ID[i]<=32 && det_ID[i] >=31)) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
for (Int_t j=0; j<det_n; j++){
if ((det_ID[j]<=50 && det_ID[j] >=47) || (det_ID[j]<=34 && det_ID[j] >=33) || (det_ID[j]<=64 && det_ID[j] >=63)) {
if (det_edep[j]>eCut){
hR->Fill(det_time_start[j]);
hDetz->Fill(det_z[j]);
hEventID->Fill(eventID);
break; //fill only once
}
}
}
}
}
}
}
// Calculate T-B Asymmetry
// Temporary T and B histograms
TH1F *hTtemp = (TH1F*) hT->Clone(); hTtemp->SetName("hTtemp");
TH1F *hBtemp = (TH1F*) hB->Clone(); hBtemp->SetName("hBtemp");
// Sum and difference T and B
TH1F *hSumTB = (TH1F*) hT->Clone(); hSumTB->SetName("hSumTB");
hSumTB->Sumw2();
TH1F *hDifTB = (TH1F*) hT->Clone(); hDifTB->SetName("hDifTB");
hDifTB->Sumw2();
// Asymmetry histograms!
TH1F *hAsyTB = (TH1F*) hT->Clone(); hAsyTB->SetName("hAsyTB");
hAsyTB->SetTitle("Muon decay asymmetry T-B; Time (#mus); Asymmetry");
hAsyTB->Sumw2();
// Calculate difference and sum, then divide
hDifTB->Add(hTtemp,hBtemp,1.,-1.);
hSumTB->Add(hTtemp,hBtemp,1., 1.);
hAsyTB->Divide(hDifTB,hSumTB,1.,1.);
// Calculate L-R Asymmetry
// Temporary L and R histograms
TH1F *hLtemp = (TH1F*) hL->Clone(); hLtemp->SetName("hLtemp");
TH1F *hRtemp = (TH1F*) hR->Clone(); hRtemp->SetName("hRtemp");
// Sum and difference L and R
TH1F *hSumLR = (TH1F*) hL->Clone(); hSumLR->SetName("hSumLR");
hSumLR->Sumw2();
TH1F *hDifLR = (TH1F*) hL->Clone(); hDifLR->SetName("hDifLR");
hDifLR->Sumw2();
// Asymmetry histograms!
TH1F *hAsyLR = (TH1F*) hL->Clone(); hAsyLR->SetName("hAsyLR");
hAsyLR->SetTitle("Muon decay asymmetry L-R; Time (#mus); Asymmetry");
hAsyLR->Sumw2();
// Calculate difference and sum, then divide
hDifLR->Add(hLtemp,hRtemp,1.,-1.);
hSumLR->Add(hLtemp,hRtemp,1., 1.);
hAsyLR->Divide(hDifLR,hSumLR,1.,1.);
if (FigFlag) {
// Top - Bottom
TCanvas* c1= new TCanvas("c1","canvas 1");
c1->Divide(2,2);
c1->cd(1);
hT->Draw();
c1->cd(2);
hB->Draw();
c1->cd(3);
hAsyTB->Draw();
hAsyTB -> Fit("pol0","Q");
// Top - Bottom
Double_t TBchi2=pol0->GetChisquare();
Double_t TBp1=pol0->GetParameter(0);
Double_t TBe1=pol0->GetParError(0);
gStyle->SetOptStat(1001111);
gStyle->SetOptFit(0001);
gStyle->SetLineColor(2);
c1->cd(4);
hDetz->Draw();
// Left - Right
TCanvas* c2= new TCanvas("c2","canvas 2");
c2->Divide(2,2);
c2->cd(1);
hL->Draw();
c2->cd(2);
hR->Draw();
c2->cd(3);
hAsyLR->Draw();
hAsyLR -> Fit("pol0","Q");
// Left - Right
Double_t LRchi2=pol0->GetChisquare();
Double_t LRp1=pol0->GetParameter(0);
Double_t LRe1=pol0->GetParError(0);
gStyle->SetOptStat(1001111);
gStyle->SetOptFit(0001);
gStyle->SetLineColor(2);
c2->cd(4);
hDetz->Draw();
} else {
hAsyTB -> Fit("pol0","NQ");
// Top - Bottom
Double_t TBchi2=pol0->GetChisquare();
Double_t TBp1=pol0->GetParameter(0);
Double_t TBe1=pol0->GetParError(0);
hAsyLR -> Fit("pol0","NQ");
// Left - Right
Double_t LRchi2=pol0->GetChisquare();
Double_t LRp1=pol0->GetParameter(0);
Double_t LRe1=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\t%g\n",TBp1,TBe1,LRp1,LRe1,NDet);
}
void NewSpec::SCoinIO( 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.,2.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]", 140, -140., 140.);
TH1F* hEventID = new TH1F("hEventID", "Event ID", 10001, -0.5, 10000.5);
// Sigmenets histograms 1-32
TH1F* h1 = new TH1F("h1","Muon arrival times in 1 (#mus)",1300,0.,13.);
hEdeposited->Sumw2();
hEdepositCF->Sumw2();
hEdepTrig->Sumw2();
Long64_t nentries = fChain->GetEntriesFast();
//nentries=1000;
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.;
// Segments (1-16) coincidence (33-48)
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==1) {
hEdeposited->Fill(det_edep[i]);
if (det_edep[i]>eCut){
for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==33) {
if (det_edep[j]>eCut){
h1->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");
// 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();
// Asymmetry histograms!
TH1F *hAsyFB = (TH1F*) hForw->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(2,2);
c1->cd(1);
hBack->Draw();
c1->cd(2);
hForw->Draw();
c1->cd(3);
hAsyFB->Draw();
hAsyFB -> Fit("pol0","Q");
gStyle->SetOptStat(1001111);
gStyle->SetOptFit(0001);
gStyle->SetLineColor(2);
c1->cd(4);
hDetz->Draw();
} else {
hAsyFB -> Fit("pol0","NQ");
}
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);
}

366
run/macros/SpinRotAna.C Normal file
View File

@ -0,0 +1,366 @@
char popis[100];
void SpinRotAna(){
cout<<"Syntax: root[]> SpinRotAna.C(\"20021\")"<<endl;
exit(0);
}
void SpinRotAna(char* runID) {
cout<<"runID="<<runID<<endl;
char inputFileName[1000]="data/musr_"; strcat(inputFileName,runID); strcat(inputFileName,".root");
cout<<"Opening file "<<inputFileName<<endl;
TFile* f1=new TFile(inputFileName);
// TFile* f1=new TFile("data/musr_20021.root");
// PlotOneSaveVolume();
strcpy(popis,"undefined");
if (strcmp(runID,"20053")==0) {strcpy(popis,"#splitline{Run 20053: E_{const}=2.43 kV/mm; v04=0.03 T}{pie3\_HiField\_d05\_x30\_z-15111}");}
if (strcmp(runID,"20063")==0) {strcpy(popis,"#splitline{Run 20063: e01=2.55 kV/mm; v04=0.03 T}{pie3\_HiField\_d05\_x30\_z-15111}");}
if (strcmp(runID,"20064")==0) {strcpy(popis,"#splitline{Run 20064: e01=2.57 kV/mm; v04=0.03 T}{pie3\_HiField\_d05\_x30\_z-15111}");}
if (strcmp(runID,"20073")==0) {strcpy(popis,"#splitline{Run 20073: e01=2.90 kV/mm; v04=0.0339 T}{pie3\_HiField\_d05\_x30\_z-15111}");}
if (strcmp(runID,"20074")==0) {strcpy(popis,"#splitline{Run 20074: e01=2.95 kV/mm; v04=0.0339 T}{pie3\_HiField\_d05\_x30\_z-15111}");}
if (strcmp(runID,"20075")==0) {strcpy(popis,"#splitline{Run 20075: e01=3.00 kV/mm; v04=0.0339 T}{pie3\_HiField\_d05\_x30\_z-15111}");}
if (strcmp(runID,"20076")==0) {strcpy(popis,"#splitline{Run 20076: e01=3.05 kV/mm; v04=0.0339 T}{pie3\_HiField\_d05\_x30\_z-15111}");}
if (strcmp(runID,"20081")==0) {strcpy(popis,"#splitline{Run 20081: e01=3.00 kV/mm; v04=0.0339 T}{pie3\_HiField\_d05\_x30\_z-15111}");}
if (strcmp(runID,"20082")==0) {strcpy(popis,"#splitline{Run 20082: e01=3.00 kV/mm; v04=0.0339 T}{pie3\_HiField\_d05\_x30\_z-15111}");}
if (strcmp(runID,"20111")==0) {strcpy(popis,"#splitline{Run 20111: e01=3.00 kV/mm; v04=0.0339 T}{pie3\_HiField\_d05\_x30\_z-15111}");}
if (strcmp(runID,"20116")==0) {strcpy(popis,"#splitline{Run 20116: e01=3.00 kV/mm; v04=0.0339 T}{pie3\_HiField\_d05\_x30\_z-15111}");}
if (strcmp(runID,"20121")==0) {strcpy(popis,"#splitline{Run 20121: e01=3.00 kV/mm; b05c=0.0339 T}{pie3\_HiField\_d05\_x30\_z-15111}");}
if (strcmp(runID,"20131")==0) {strcpy(popis,"#splitline{Run 20131: e01=3.00 kV/mm; b05c=0.0339 T}{reggiani\_Jan2010\_z-12032}");}
if (strcmp(runID,"20132")==0) {strcpy(popis,"#splitline{Run 20132: e01=3.03 kV/mm; b05c=0.0339 T}{reggiani\_Jan2010\_z-12032}");}
if (strcmp(runID,"20133")==0) {strcpy(popis,"#splitline{Run 20133: e01=3.06 kV/mm; b05c=0.0339 T}{reggiani\_Jan2010\_z-12032}");}
if (strcmp(runID,"20134")==0) {strcpy(popis,"#splitline{Run 20134: e01=2.97 kV/mm; b05c=0.0339 T}{reggiani\_Jan2010\_z-12032}");}
if (strcmp(runID,"20135")==0) {strcpy(popis,"#splitline{Run 20135: e01=2.94 kV/mm; b05c=0.0339 T}{reggiani\_Jan2010\_z-12032}");}
if (strcmp(runID,"20141")==0) {strcpy(popis,"#splitline{Run 20141: e01=3.00 kV/mm; v04=0.0339 T}{reggiani\_Jan2010\_z-12032}");}
if (strcmp(runID,"20142")==0) {strcpy(popis,"#splitline{Run 20142: e01=3.03 kV/mm; v04=0.0339 T}{reggiani\_Jan2010\_z-12032}");}
if (strcmp(runID,"20143")==0) {strcpy(popis,"#splitline{Run 20143: e01=3.06 kV/mm; v04=0.0339 T}{reggiani\_Jan2010\_z-12032}");}
if (strcmp(runID,"20144")==0) {strcpy(popis,"#splitline{Run 20144: e01=2.97 kV/mm; v04=0.0339 T}{reggiani\_Jan2010\_z-12032}");}
if (strcmp(runID,"20145")==0) {strcpy(popis,"#splitline{Run 20145: e01=2.94 kV/mm; v04=0.0339 T}{reggiani\_Jan2010\_z-12032}");}
if (strcmp(runID,"20181")==0) {strcpy(popis,"#splitline{Run 20181: e01=3.09 kV/mm; b05c=0.0339 T}{reggiani\_Jan2010\_z-12032}");}
if (strcmp(runID,"20191")==0) {strcpy(popis,"#splitline{Run 20191: e01=3.09 kV/mm; v04=0.0339 T}{reggiani\_Jan2010\_z-12032}");}
if (strcmp(runID,"20201")==0) {strcpy(popis,"#splitline{Run 20201: e01=3.09 kV/mm; b05c=0.0339 T}{reggiani\_Jan2010\_z-12032, y-x}");}
if (strcmp(runID,"20201")==0) {strcpy(popis,"#splitline{Run 20201: e01=3.09 kV/mm; b05c=0.0339 T}{reggiani\_Jan2010\_z-12032, y-x}");}
if (strcmp(runID,"20202")==0) {strcpy(popis,"#splitline{Run 20202: e01=3.09 kV/mm; b05c=0.0339 T}{reggiani\_Jan2010\_z-12032, -yx}");}
if (strcmp(runID,"20203")==0) {strcpy(popis,"#splitline{Run 20203: e01=3.09 kV/mm; b05c=0.0339 T}{reggiani\_Jan2010\_z-12032}");}
if (strcmp(runID,"20211")==0) {strcpy(popis,"#splitline{Run 20211: e01=3.09 kV/mm; b05c=0.0339 T}{reggiani\_Jan2010\_z-12032, y-x}");}
if (strcmp(runID,"20212")==0) {strcpy(popis,"#splitline{Run 20212: e01=3.09 kV/mm; b05c=0.0339 T}{reggiani\_Jan2010\_z-12032, -yx}");}
if (strcmp(runID,"20213")==0) {strcpy(popis,"#splitline{Run 20213: e01=3.09 kV/mm; b05c=0.0339 T}{reggiani\_Jan2010\_z-12032}");}
if (strcmp(runID,"20221")==0) {strcpy(popis,"#splitline{Run 20221: e01=3.09 kV/mm; b05c=0.0339 T}{reggiani\_Jan2010\_z-12032\_NEW, y-x}");}
if (strcmp(runID,"20222")==0) {strcpy(popis,"#splitline{Run 20222: e01=3.09 kV/mm; b05c=0.0339 T}{reggiani\_Jan2010\_z-12032\_NEW, -yx}");}
if (strcmp(runID,"20223")==0) {strcpy(popis,"#splitline{Run 20223: e01=3.09 kV/mm; b05c=0.0339 T}{reggiani\_Jan2010\_z-12032\_NEW}");}
if (strcmp(runID,"20231")==0) {strcpy(popis,"#splitline{Run 20231: e01=3.09 kV/mm; b05c=0.0339 T}{reggiani\_Jan2010\_z-12032\_NEW}");}
if (strcmp(runID,"20301")==0) {strcpy(popis,"#splitline{Run 20301: e01=3.09 kV/mm; b05c=0.0339 T}{like 20181 with Geant4.9.3}");}
if (strcmp(runID,"20302")==0) {strcpy(popis,"#splitline{Run 20302: e01=3.09 kV/mm; b05c=0.0339 T}{like 20181}");}
if (strcmp(runID,"20303")==0) {strcpy(popis,"#splitline{Run 20303: e01=3.09 kV/mm; b05c=0.0339 T}{like 20181 with Geant4.9.3 new}");}
if (strcmp(runID,"20311")==0) {strcpy(popis,"#splitline{Run 20311: e01=3.09 kV/mm; b05c=0.0339 T}{Feb Turtle, Geant4.9.3}");}
if (strcmp(runID,"20312")==0) {strcpy(popis,"#splitline{Run 20312: e01=3.00 kV/mm; b05c=0.0339 T}{Feb Turtle, Geant4.9.3}");}
if (strcmp(runID,"20313")==0) {strcpy(popis,"#splitline{Run 20313: e01=3.03 kV/mm; b05c=0.0339 T}{Feb Turtle, Geant4.9.3}");}
if (strcmp(runID,"20314")==0) {strcpy(popis,"#splitline{Run 20314: e01=3.06 kV/mm; b05c=0.0339 T}{Feb Turtle, Geant4.9.3}");}
if (strcmp(runID,"20315")==0) {strcpy(popis,"#splitline{Run 20315: e01=3.12 kV/mm; b05c=0.0339 T}{Feb Turtle, Geant4.9.3}");}
if (strcmp(runID,"20316")==0) {strcpy(popis,"#splitline{Run 20316: e01=3.15 kV/mm; b05c=0.0339 T}{Feb Turtle, Geant4.9.3}");}
if (strcmp(runID,"20321")==0) {strcpy(popis,"#splitline{Run 20321: e01=3.03 kV/mm; b05c=0.0339 T}{Feb Turtle, Geant4.9.3}");}
if (strcmp(runID,"20331")==0) {strcpy(popis,"#splitline{Run 20331: E_{unif}=3.03 kV/mm; B_{unif}=0.0392 T}{Feb Turtle, Geant4.9.3}");}
if (strcmp(runID,"20332")==0) {strcpy(popis,"#splitline{Run 20332: E_{unif}=3.00 kV/mm; B_{unif}=0.0392 T}{Feb Turtle, Geant4.9.3}");}
if (strcmp(runID,"20333")==0) {strcpy(popis,"#splitline{Run 20333: E_{unif}=3.06 kV/mm; B_{unif}=0.0392 T}{Feb Turtle, Geant4.9.3}");}
if (strcmp(runID,"20334")==0) {strcpy(popis,"#splitline{Run 20334: E_{unif}=3.09 kV/mm; B_{unif}=0.0392 T}{Feb Turtle, Geant4.9.3}");}
if (strcmp(runID,"20335")==0) {strcpy(popis,"#splitline{Run 20335: E_{unif}=2.97 kV/mm; B_{unif}=0.0392 T}{Feb Turtle, Geant4.9.3}");}
if (strcmp(runID,"20351")==0) {strcpy(popis,"#splitline{Run 20351: E_{unif}=3.08544 kV/mm; B_{unif}=0.041 T}{Feb Turtle}");}
if (strcmp(runID,"20361")==0) {strcpy(popis,"#splitline{Run 20361: E_{unif}=3.08544 kV/mm; B_{unif}=0.041 T}{p=27.4 MeV/c, x=y=0}");}
if (strcmp(runID,"20362")==0) {strcpy(popis,"#splitline{Run 20362: E_{unif}=3.08544 kV/mm; B_{unif}=0.041 T}{p=27.4+/-0.81 MeV/c, x=y=0}");}
if (strcmp(runID,"20363")==0) {strcpy(popis,"#splitline{Run 20363: E_{unif}=3.08544 kV/mm; B_{unif}=0.041 T}{p=27.4 MeV/c, x,y=Gaussian}");}
if (strcmp(runID,"20371")==0) {strcpy(popis,"#splitline{Run 20371: e01=3.03 kV/mm; b05c=0.0339 T}{Feb Turtle, Geant4.9.3}");}
if (strcmp(runID,"20381")==0) {strcpy(popis,"#splitline{Run 20381: E_{unif}=3.08544 kV/mm; B_{unif}=0.041 T}{Feb Turtle, #mu does not decay}");}
if (strcmp(runID,"20391")==0) {strcpy(popis,"#splitline{Run 20391: E_{unif}=3.08544 kV/mm; B_{unif}=0.041 T}{p=27.4 MeV/c, x=y=0, #mu does not decay}");}
if (strcmp(runID,"20392")==0) {strcpy(popis,"#splitline{Run 20392: E_{unif}=3.08544 kV/mm; B_{unif}=0.041 T}{p=27.4+/-0.81 MeV/c, x=y=0, #mu does not decay}");}
if (strcmp(runID,"20393")==0) {strcpy(popis,"#splitline{Run 20393: E_{unif}=3.08544 kV/mm; B_{unif}=0.041 T}{p=27.4 MeV/c, x,y=Gaussian, #mu does not decay}");}
TCanvas* cR1IN = PlotOneSaveVolume ("R1IN", "save_detID==1010&&save_particleID==-13");
TCanvas* cR1IN2= PlotOneSaveVolumeTwo("R1IN2","save_detID==1010&&save_particleID==-13");
char outputFileName_cR1IN[1000]; sprintf(outputFileName_cR1IN,"pict/SpinRotAna_R1IN_%s.ps",runID);
cR1IN->Print(outputFileName_cR1IN);
char outputFileName_cR1IN2[1000]; sprintf(outputFileName_cR1IN2,"pict/SpinRotAna_R1IN_two_%s.ps",runID);
cR1IN2->Print(outputFileName_cR1IN2);
TCanvas* cR1OU = PlotOneSaveVolume ("R1OU", "save_detID==1020&&save_particleID==-13");
TCanvas* cR1OU2= PlotOneSaveVolumeTwo("R1OU2","save_detID==1020&&save_particleID==-13");
char outputFileName_cR1OU[1000]; sprintf(outputFileName_cR1OU,"pict/SpinRotAna_R1OU_%s.ps",runID);
cR1OU->Print(outputFileName_cR1OU);
char outputFileName_cR1OU2[1000]; sprintf(outputFileName_cR1OU2,"pict/SpinRotAna_R1OU_two_%s.ps",runID);
cR1OU2->Print(outputFileName_cR1OU2);
TCanvas* cR2IN = PlotOneSaveVolume ("R2IN", "save_detID==1050&&save_particleID==-13");
TCanvas* cR2IN2= PlotOneSaveVolumeTwo("R2IN2","save_detID==1050&&save_particleID==-13");
char outputFileName_cR2IN[1000]; sprintf(outputFileName_cR2IN,"pict/SpinRotAna_R2IN_%s.ps",runID);
cR2IN->Print(outputFileName_cR2IN);
char outputFileName_cR2IN2[1000]; sprintf(outputFileName_cR2IN2,"pict/SpinRotAna_R2IN_two_%s.ps",runID);
cR2IN2->Print(outputFileName_cR2IN2);
TCanvas* cR2OU = PlotOneSaveVolume ("R2OU", "save_detID==1060&&save_particleID==-13");
TCanvas* cR2OU2= PlotOneSaveVolumeTwo("R2OU2","save_detID==1060&&save_particleID==-13");
char outputFileName_cR2OU[1000]; sprintf(outputFileName_cR2OU,"pict/SpinRotAna_R2OU_%s.ps",runID);
cR2OU->Print(outputFileName_cR2OU);
char outputFileName_cR2OU2[1000]; sprintf(outputFileName_cR2OU2,"pict/SpinRotAna_R2OU_two_%s.ps",runID);
cR2OU2->Print(outputFileName_cR2OU2);
TCanvas* cEnvelope = new TCanvas("cEnvelope","cEnvelope");
// TH1D* hisEnvelopeX = new TH1D("hisEnvelopeX","Beam envelope; z (cm); 2*RMS of x and -y (cm)",1200,-1700,-500);
// TH1D* hisEnvelopeY = new TH1D("hisEnvelopeY","Beam envelope; z (cm); 2*RMS of x and -y (cm)",1200,-1700,-500);
TH1D* hisEnvelopeX = new TH1D("hisEnvelopeX","Beam envelope; z (cm); #pm 2*RMS of x and y (cm)",1200,-1400,-200);
TH1D* hisEnvelopeY = new TH1D("hisEnvelopeY","Beam envelope; z (cm); #pm 2*RMS of x and y (cm)",1200,-1400,-200);
Double_t fRunID=atof(runID);
std::cout<<"fRunID="<<fRunID<<std::endl;
if (fRunID<20130) {
FillEnvelope(hisEnvelopeX,-1511.15,"save_x/10>>hh","save_detID==1010&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-1241.0,"save_x/10>>hh","save_detID==1020&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-1200,"save_x/10>>hh","save_detID==1025&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-1160,"save_x/10>>hh","save_detID==1028&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-1146.6,"save_x/10>>hh","save_detID==1030&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-1086.6,"save_x/10>>hh","save_detID==1033&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-1026.6,"save_x/10>>hh","save_detID==1036&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-972.15,"save_x/10>>hh","save_detID==1040&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-962.15,"save_x/10>>hh","save_detID==1041&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-952.15,"save_x/10>>hh","save_detID==1042&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-942.15,"save_x/10>>hh","save_detID==1043&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-932.15,"save_x/10>>hh","save_detID==1044&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-922.15,"save_x/10>>hh","save_detID==1045&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-912.15,"save_x/10>>hh","save_detID==1046&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-902.15,"save_x/10>>hh","save_detID==1047&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-892.15,"save_x/10>>hh","save_detID==1048&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-882.15,"save_x/10>>hh","save_detID==1049&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-872.15,"save_x/10>>hh","save_detID==1050&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-602.0 ,"save_x/10>>hh","save_detID==1060&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-592.0 ,"save_x/10>>hh","save_detID==1061&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-582.0 ,"save_x/10>>hh","save_detID==1062&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-572.0 ,"save_x/10>>hh","save_detID==1063&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-562.0 ,"save_x/10>>hh","save_detID==1064&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-552.0 ,"save_x/10>>hh","save_detID==1065&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-542.0 ,"save_x/10>>hh","save_detID==1066&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-532.0 ,"save_x/10>>hh","save_detID==1067&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-522.0 ,"save_x/10>>hh","save_detID==1068&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-512.0 ,"save_x/10>>hh","save_detID==1069&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-1511.15,"save_y/10>>hh","save_detID==1010&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-1241.0,"save_y/10>>hh","save_detID==1020&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-1200,"save_y/10>>hh","save_detID==1025&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-1160,"save_y/10>>hh","save_detID==1028&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-1146.6,"save_y/10>>hh","save_detID==1030&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-1086.6,"save_y/10>>hh","save_detID==1033&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-1026.6,"save_y/10>>hh","save_detID==1036&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-972.15,"save_y/10>>hh","save_detID==1040&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-962.15,"save_y/10>>hh","save_detID==1041&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-952.15,"save_y/10>>hh","save_detID==1042&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-942.15,"save_y/10>>hh","save_detID==1043&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-932.15,"save_y/10>>hh","save_detID==1044&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-922.15,"save_y/10>>hh","save_detID==1045&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-912.15,"save_y/10>>hh","save_detID==1046&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-902.15,"save_y/10>>hh","save_detID==1047&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-892.15,"save_y/10>>hh","save_detID==1048&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-882.15,"save_y/10>>hh","save_detID==1049&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-872.15,"save_y/10>>hh","save_detID==1050&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-602.0 ,"save_y/10>>hh","save_detID==1060&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-592.0 ,"save_y/10>>hh","save_detID==1061&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-582.0 ,"save_y/10>>hh","save_detID==1062&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-572.0 ,"save_y/10>>hh","save_detID==1063&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-562.0 ,"save_y/10>>hh","save_detID==1064&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-552.0 ,"save_y/10>>hh","save_detID==1065&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-542.0 ,"save_y/10>>hh","save_detID==1066&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-532.0 ,"save_y/10>>hh","save_detID==1067&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-522.0 ,"save_y/10>>hh","save_detID==1068&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-512.0 ,"save_y/10>>hh","save_detID==1069&&save_particleID==-13");
}
else {
FillEnvelope(hisEnvelopeX,-1183.3,"save_x/10>>hh","save_detID==1010&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-913.1,"save_x/10>>hh","save_detID==1020&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-903.1,"save_x/10>>hh","save_detID==1025&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-893.1,"save_x/10>>hh","save_detID==1026&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-880.0,"save_x/10>>hh","save_detID==1028&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-870.0,"save_x/10>>hh","save_detID==1030&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-820.0,"save_x/10>>hh","save_detID==1033&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-759.0,"save_x/10>>hh","save_detID==1036&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-710.0,"save_x/10>>hh","save_detID==1045&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-707.0,"save_x/10>>hh","save_detID==1046&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-697.0,"save_x/10>>hh","save_detID==1047&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-687.0,"save_x/10>>hh","save_detID==1048&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-677.0,"save_x/10>>hh","save_detID==1049&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-667.0,"save_x/10>>hh","save_detID==1050&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-395.0,"save_x/10>>hh","save_detID==1060&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-385.0,"save_x/10>>hh","save_detID==1061&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-375.0,"save_x/10>>hh","save_detID==1062&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-365.0,"save_x/10>>hh","save_detID==1063&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-355.0,"save_x/10>>hh","save_detID==1064&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-345.0,"save_x/10>>hh","save_detID==1065&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-335.0,"save_x/10>>hh","save_detID==1066&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-325.0,"save_x/10>>hh","save_detID==1067&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-315.0,"save_x/10>>hh","save_detID==1068&&save_particleID==-13");
FillEnvelope(hisEnvelopeX,-305.0,"save_x/10>>hh","save_detID==1069&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-1183.3,"save_y/10>>hh","save_detID==1010&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-913.1,"save_y/10>>hh","save_detID==1020&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-903.1,"save_y/10>>hh","save_detID==1025&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-893.1,"save_y/10>>hh","save_detID==1026&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-880.0,"save_y/10>>hh","save_detID==1028&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-870.0,"save_y/10>>hh","save_detID==1030&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-820.0,"save_y/10>>hh","save_detID==1033&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-759.0,"save_y/10>>hh","save_detID==1036&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-710.0,"save_y/10>>hh","save_detID==1045&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-707.0,"save_y/10>>hh","save_detID==1046&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-697.0,"save_y/10>>hh","save_detID==1047&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-687.0,"save_y/10>>hh","save_detID==1048&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-677.0,"save_y/10>>hh","save_detID==1049&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-667.0,"save_y/10>>hh","save_detID==1050&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-395.0,"save_y/10>>hh","save_detID==1060&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-385.0,"save_y/10>>hh","save_detID==1061&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-375.0,"save_y/10>>hh","save_detID==1062&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-365.0,"save_y/10>>hh","save_detID==1063&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-355.0,"save_y/10>>hh","save_detID==1064&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-345.0,"save_y/10>>hh","save_detID==1065&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-335.0,"save_y/10>>hh","save_detID==1066&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-325.0,"save_y/10>>hh","save_detID==1067&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-315.0,"save_y/10>>hh","save_detID==1068&&save_particleID==-13");
FillEnvelope(hisEnvelopeY,-305.0,"save_y/10>>hh","save_detID==1069&&save_particleID==-13");
}
hisEnvelopeX->SetMaximum(19);
hisEnvelopeX->SetMinimum(-19);
hisEnvelopeX->SetMarkerStyle(20);
hisEnvelopeX->Draw("P");
hisEnvelopeY->SetMarkerStyle(20);
hisEnvelopeY->Draw("Psame");
Vysvetlivka(popis,0.15,0.14,90,0.025);
char outputFileName_cR2OU[1000]; sprintf(outputFileName_cR2OU,"pict/SpinRotAna_Envelope_%s.ps",runID);
gPad->SetGrid();
cEnvelope->Print(outputFileName_cR2OU);
}
TCanvas* PlotOneSaveVolume(char* name, char* condition) {
cout<<"PlotOneSaveVolume called"<<endl;
char cName[100]; strcpy(cName,"c"); strcat(cName,name);
TCanvas* c1= new TCanvas(cName,cName,700,900);
c1->Divide(2,3);
c1->cd(1);
char hNameX[200]; strcpy(hNameX,"h"); strcat(hNameX,name); strcat(hNameX,"x");
TH1D* hisX = new TH1D(hNameX,hNameX,100,-25,25);
hisX->GetXaxis()->SetTitle("x (cm)");
char whatToPlotX[200]; strcpy(whatToPlotX,"save_x/10>>"); strcat(whatToPlotX,hNameX);
t1->Draw(whatToPlotX,condition);
if (strcmp(popis,"undefined")) {
Vysvetlivka(popis,0.2,0.14,90);
}
char chTransmission[200];
Double_t transmission = (hisX->GetEntries())/(t1->GetEntries())*100;
sprintf(chTransmission,"Transmission: %.1f %%",transmission);
Vysvetlivka(chTransmission,0.8,0.14,90);
c1->cd(2);
char hNameY[200]; strcpy(hNameY,"h"); strcat(hNameY,name); strcat(hNameY,"y");
cout<<hNameY<<endl;
TH1D* hisY = new TH1D(hNameY,hNameY,100,-25,25);
hisY->GetXaxis()->SetTitle("y (cm)");
char whatToPlotY[200]; strcpy(whatToPlotY,"save_y/10>>"); strcat(whatToPlotY,hNameY);
t1->Draw(whatToPlotY,condition);
c1->cd(3);
char hNameXX[200]; strcpy(hNameXX,"h"); strcat(hNameXX,name); strcat(hNameXX,"Xprime");
cout<<hNameXX<<endl;
TH1D* hisXX = new TH1D(hNameXX,hNameXX,100,-100,100);
hisXX->GetXaxis()->SetTitle("x_{angle} (mrad)");
char whatToPlotXX[200]; strcpy(whatToPlotXX,"atan2(save_px,save_pz)*1000>>"); strcat(whatToPlotXX,hNameXX);
cout<<whatToPlotXX<<endl;
t1->Draw(whatToPlotXX,condition);
c1->cd(4);
char hNameYY[200]; strcpy(hNameYY,"h"); strcat(hNameYY,name); strcat(hNameYY,"Yprime");
cout<<hNameYY<<endl;
TH1D* hisYY = new TH1D(hNameYY,hNameYY,100,-100,100);
hisYY->GetXaxis()->SetTitle("y_{angle} (mrad)");
char whatToPlotYY[200]; strcpy(whatToPlotYY,"atan2(save_py,save_pz)*1000>>"); strcat(whatToPlotYY,hNameYY);
cout<<whatToPlotYY<<endl;
t1->Draw(whatToPlotYY,condition);
// c1->cd(5);
// char hNamePolT[200]; strcpy(hNamePolT,"h"); strcat(hNamePolT,name); strcat(hNamePolT,"PolT");
// cout<<hNamePolT<<endl;
// TH1D* hisPolT = new TH1D(hNamePolT,hNamePolT,102,-1.02,1.02);
// hisPolT->GetXaxis()->SetTitle("#sqrt{(P_{x}^{2}+P_{y}^{2})}");
// char whatToPlotPolT[200]; strcpy(whatToPlotPolT,"sqrt(save_polx*save_polx+save_poly*save_poly)>>"); strcat(whatToPlotPolT,hNamePolT);
// cout<<whatToPlotPolT<<endl;
// t1->Draw(whatToPlotPolT,condition);
c1->cd(5);
char hNameThetaPol[200]; strcpy(hNameThetaPol,"h"); strcat(hNameThetaPol,name); strcat(hNameThetaPol,"ThetaPol");
cout<<hNameThetaPol<<endl;
TH1D* hisThetaPol = new TH1D(hNameThetaPol,hNameThetaPol,90,-90.,90);
hisThetaPol->GetXaxis()->SetTitle("#theta(P)-90\,deg");
// char whatToPlotPolT[200]; strcpy(whatToPlotPolT,"180/3.14159*atan(sqrt(save_polx*save_polx+save_poly*save_poly)/save_polz)>>"); strcat(whatToPlotPolT,hNameThetaPol);
char whatToPlotPolT[200]; strcpy(whatToPlotPolT,"180/3.14159*atan(save_polz/sqrt(save_polx*save_polx+save_poly*save_poly))>>"); strcat(whatToPlotPolT,hNameThetaPol);
cout<<whatToPlotPolT<<endl;
t1->Draw(whatToPlotPolT,condition);
c1->cd(6);
char hNamePz[200]; strcpy(hNamePz,"h"); strcat(hNamePz,name); strcat(hNamePz,"pz");
TH1D* hisPz = new TH1D(hNamePz,hNamePz,100,-50,50);
hisPz->GetXaxis()->SetTitle("(p-p0)/p0 (%)");
char whatToPlotPz[200]; strcpy(whatToPlotPz,"(save_pz-27.4)/27.4*100>>"); strcat(whatToPlotPz,hNamePz);
t1->Draw(whatToPlotPz,condition);
return c1;
}
TCanvas* PlotOneSaveVolumeTwo(char* name, char* condition) {
cout<<"PlotOneSaveVolume called"<<endl;
char cName[100]; strcpy(cName,"c"); strcat(cName,name);
TCanvas* c1= new TCanvas(cName,cName,700,900);
c1->Divide(2,3);
c1->cd(1);
char hName1[200]; strcpy(hName1,"h"); strcat(hName1,name); strcat(hName1,"xXxprime");
TH2D* his1 = new TH2D(hName1,hName1,40,-25,25,40,-25,25);
his1->GetXaxis()->SetTitle("x (cm)");
his1->GetYaxis()->SetTitle("y (cm)");
char whatToPlotX[200]; strcpy(whatToPlotX,"save_y/10:save_x/10>>"); strcat(whatToPlotX,hName1);
t1->Draw(whatToPlotX,condition,"box");
if (strcmp(popis,"undefined")) {
Vysvetlivka(popis,0.2,0.14,90);
}
c1->cd(2);
char hName2[200]; strcpy(hName2,"h"); strcat(hName2,name); strcat(hName2,"xXxprime");
TH2D* his2 = new TH2D(hName2,hName2,40,-25,25,40,-100,100);
his2->GetXaxis()->SetTitle("x (cm)");
his2->GetYaxis()->SetTitle("x_{angle} (mrad)");
char whatToPlotX[200]; strcpy(whatToPlotX,"(atan2(save_px,save_pz)*1000):save_x/10>>"); strcat(whatToPlotX,hName2);
t1->Draw(whatToPlotX,condition,"box");
c1->cd(3);
char hName3[200]; strcpy(hName3,"h"); strcat(hName3,name); strcat(hName3,"xXp");
TH2D* his3 = new TH2D(hName3,hName3,40,-25,25,40,-50,50);
his3->GetXaxis()->SetTitle("x (cm)");
his3->GetYaxis()->SetTitle("(p-p0)/p0 (%)");
char whatToPlotX[200]; strcpy(whatToPlotX,"((save_pz-27.4)/27.4*100):(save_x/10)>>"); strcat(whatToPlotX,hName3);
t1->Draw(whatToPlotX,condition,"box");
return c1;
}
void FillEnvelope(TH1D* hist, const Double_t z, char* velicina, char* condition) {
TH1D hh = TH1D("hh","hh",1000,-50.,50.);
t1->Draw(velicina,condition);
Double_t RMS2 = 2*(hh.GetRMS());
if (strcmp(velicina,"save_y/10>>hh")==0) {RMS2 *=-1;}
std::cout<<"z="<<z<<" RMS2="<<RMS2<<std::endl;
hist->Fill(z,RMS2);
}
void Vysvetlivka(char* textik, Double_t x, Double_t y, Double_t uhel=270, Double_t textSize=0.04) {
TLatex latex;
latex.SetNDC();
latex.SetTextAngle(uhel);
latex.SetTextSize(textSize);
latex.DrawLatex(x,y,textik);
}

4
run/macros/TFrame.C Normal file
View File

@ -0,0 +1,4 @@
{
//========= Macro generated from object: TFrame/Pad graphics frame
//========= by ROOT version5.24/00
}

228
run/macros/analysis.C Normal file
View File

@ -0,0 +1,228 @@
#define analysis_cxx
#include "analysis.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
void analysis::Loop()
{
// In a ROOT session, you can do:
// Root > .L analysis.C
// Root > analysis t
// Root > t.GetEntry(12); // Fill t data members with entry number 12
// Root > t.Show(); // Show values of entry 12
// Root > t.Show(16); // Read and show values of entry 16
// Root > t.Loop(); // Loop on all entries <- Very Important!
// Root > .ls // Shows what histos are available for plotting
//
// This is the loop skeleton where:
// jentry is the global entry number in the chain
// ientry is the entry number in the current Tree
// Note that the argument to GetEntry must be:
// jentry for TChain::GetEntry
// ientry for TTree::GetEntry and TBranch::GetEntry
//
// To read only selected branches, insert statements as:
// METHOD1:
// fChain->SetBranchStatus("*",0); // disable all branches
// fChain->SetBranchStatus("branchname",1); // activate branchname
// METHOD2: replace line
// fChain->GetEntry(jentry); //read all branches
// b_branchname->GetEntry(ientry); //read only this branch
//
//
// To create TEMPLATE files for analysis proceed as follows:
// 1. Load the root file containing all the variables, e.g.:
// TFile* f=new TFile("data/lem4_1049.root");
// 2. Quickly check the Tree and Variable names:
// f->ls(); // Shows tree name (i.e. t1)
// t1->Print(); // Shows variable names
// 3. Create the template files .C and .h (t1 is the Tree name)
// t1->MakeClass("myAnalysis")
//
//
// How to write data to an EXTERNAL file (use cout instead of outdata for printing on screen):
// See also: www.bgsu.edu/departments/compsci/docs/index.html
//
// ofstream outdata;
// outdata.open("export.dat");
// for (Int_t i=0; i<hEdepositCF->GetNbinsX(); i++) {outdata<<hEdepositCF->GetBinContent(i)<<endl;}
// outdata.close();
//
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.,2.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* hDettUp = new TH1F("hDettUp", "Muon arrival times Up (#mus)", 1300,0.,13.);
TH1F* hDettDown = new TH1F("hDettDown", "Muon arrival times Down (#mus)", 1300,0.,13.);
TH1F* hDettRight = new TH1F("hDettRight", "Muon arrival times Right (#mus)", 1300,0.,13.);
TH1F* hDettLeft = new TH1F("hDettLeft", "Muon arrival times Left (#mus)", 1300,0.,13.);
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]", 140, -140., 140.);
TH1F* hEventID = new TH1F("hEventID", "Event ID", 10001, -0.5, 10000.5);
hEdeposited->Sumw2();
hEdepositCF->Sumw2();
hEdepTrig->Sumw2();
Long64_t nentries = fChain->GetEntriesFast();
//nentries=1000;
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.;
// Coincidence between Up scintillators
// for (Int_t i=0; i<det_n; i++) {
// if (det_ID[i]==21) {
// hEdeposited->Fill(det_edep[i]);
// for (Int_t j=0; j<det_n; j++){
// if (det_ID[j]==11) {
// if (det_edep[i]>eCut && det_edep[j]>eCut && tofFlag){
// hDettUp->Fill(det_time_start[j]);
// hDetz->Fill(det_z[i]);
// hEventID->Fill(eventID);
// break; //fill only once
// }
// }
// }
// }
}
// Coincidence between Down scintillators
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==23) {
hEdeposited->Fill(det_edep[i]);
for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==13){
if (det_edep[i]>eCut && det_edep[j]>eCut && tofFlag){
hDettDown->Fill(det_time_start[j]);
hDetz->Fill(det_z[i]);
hEventID->Fill(eventID);
break;
}
}
}
}
}
// Coincidence between Right scintillators
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==22) {
hEdeposited->Fill(det_edep[i]);
for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==12) {
if (det_edep[i]>eCut && det_edep[j]>eCut && tofFlag){
hDettRight->Fill(det_time_start[j]);
hDetz->Fill(det_z[i]);
hEventID->Fill(eventID);
break;
}
}
}
}
}
// Coincidence between Left scintillators
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==24) {
hEdeposited->Fill(det_edep[i]);
for (Int_t j=0; j<det_n; j++){
if (det_ID[j]==14){
if (det_edep[i]>eCut && det_edep[j]>eCut && tofFlag){
hDettLeft->Fill(det_time_start[j]);
hDetz->Fill(det_z[i]);
hEventID->Fill(eventID);
break;
}
}
}
}
}
// Testing event coincidence with SciFi
Bool_t twentyone= false;
Bool_t twozero7 = false;
Bool_t twozero8 = false;
Int_t ikk=0;
for (Int_t i=0; i<det_n; i++) {
if ((det_ID[i]== 21)&&(det_edep[i]>0.0)) {twentyone = true; ikk = i;}
if ((det_ID[i]==207)&&(det_edep[i]>0.0)) twozero7 = true;
if ((det_ID[i]==208)&&(det_edep[i]>0.0)) twozero8 = true;
}
//if (twentyone) hEdepoTest->Fill(det_edep[ikk]); // Fiber hits
//if (twozero7&& twozero8) hEdepoTest->Fill(det_edep[ikk]); // Coincidence
if (twentyone&&twozero7&& twozero8) hEdepoTest->Fill(det_edep[ikk]); // Fibre hits + coincidence
}
// PLOT HISTOGRAMS
TCanvas* c1= new TCanvas("c1","canvas 1");
// hEdeposited->Scale(1/(number_of_generated_events*hEdeposited->GetBinWidth(1)));
TH1F *htemp1 = (TH1F*) hDettUp->Clone(); htemp1->SetName("htemp1");
TH1F *htemp2 = (TH1F*) hDettUp->Clone(); htemp2->SetName("htemp2");
TH1F *htemp3 = (TH1F*) hDettUp->Clone(); htemp3->SetName("htemp3");
TH1F *htemp4 = (TH1F*) hDettUp->Clone(); htemp4->SetName("htemp4");
TH1F *htemp5 = (TH1F*) hDettUp->Clone(); htemp5->SetName("htemp5");
htemp5->SetTitle("Muon decay asymmetry U-D; Time (#mus); Asymmetry");
htemp5->Sumw2();
TH1F *htemp6 = (TH1F*) hDettUp->Clone(); htemp6->SetName("htemp6");
htemp6->SetTitle("Muon decay asymmetry L-R; Time (#mus); Asymmetry");
htemp6->Sumw2();
htemp1->Add(hDettUp,hDettDown,1.,-1.);
htemp2->Add(hDettUp,hDettDown,1., 1.);
htemp5->Divide(htemp1,htemp2,1.,1.);
htemp3->Add(hDettLeft,hDettRight,1.,-1.);
htemp4->Add(hDettLeft,hDettRight,1., 1.);
htemp6->Divide(htemp3,htemp4,1.,1.);
htemp5->SetLineWidth(2);
htemp5->SetLineColor(2);
htemp5->SetAxisRange(-0.8,0.8,"Y");
htemp5->Draw("hist");
htemp6->SetLineWidth(2);
htemp6->SetLineColor(4);
htemp6->SetAxisRange(-0.8,0.8,"Y");
htemp6->Draw("same");
hEdepositCF->SetLineWidth(2);
hEdepositCF->SetLineColor(2);
hEdepositCF->Draw("hist");
//hEdepositCF->Fit("gaus");
TCanvas* c2= new TCanvas("c2","canvas 2");
hEdeposited->Draw();
// sprintf(filenamePrint,"musrSrClass_%s.ps",fileNr);
// c1->Print(filenamePrint);
}

44
run/macros/invertZ.C Normal file
View File

@ -0,0 +1,44 @@
{
TString fileIn, fileOut;
Double_t *r, *z, *Er, *Ez;
Int_t fileSize, rValues, zValues, index;
fileIn = "NosePlate_2DE.map";
fileOut = "NosePlate_rawscaled_2DE.map";
TGraphErrors *fileGraph = new TGraphErrors(fileIn);
if ( fileGraph->IsZombie() ){
cout << endl;
cout << "File " << fileIn << " does not exist!!!" << endl << endl;
return;
}
r = fileGraph->GetX();
z = fileGraph->GetY();
Er = fileGraph->GetEX();
Ez = fileGraph->GetEY();
fileSize = fileGraph->GetN();
rValues = r[0];
zValues = z[0];
Double_t zScale[fileSize];
Double_t EzScale[fileSize];
zScale[0] = z[0];
EzScale[0] = Ez[0];
for (Int_t i=1; i<fileSize; i++){
zScale[i] = z[i]*-1;
EzScale[i] = Ez[i]*-1.;
}
FILE *fp = fopen(fileOut, "w");
fprintf(fp,"%4.0f\t%4.0f\t mm\t0.00001000\n",rValues, zValues);
for (Int_t j=0; j<rValues; j++){
for (i=zValues; i>0; i--){
index = j*zValues + i;
fprintf(fp,"%4.0f\t%4.0f\t%10.3f\t%10.3f\n",r[index], zScale[index], Er[index], EzScale[index]);
}
}
fclose(fp);
}

309
run/macros/sr_tr.C Normal file
View File

@ -0,0 +1,309 @@
// To run root without logon screen and script use: root -l -n
// To run a named macro in root use:
//.L sr_tr.C
// sr_transm("data/sr1_1100.root")
#include <string>
void sr_transm(char* fname)
{
TFile* f1 = new TFile(fname);
gStyle->SetPalette(1,0);
TCanvas* c1 = new TCanvas("c1","Polarization and Transmission",60,40,600,800);
c1->Divide(2,3);
c1->SetGridx();
c1->SetGridy();
// Resulting polarization in the x direction
c1->cd(1);
t1->Draw("save_polx >> pol_x(128, -1, 1)","save_detID==902&&save_particleID==-13"); //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");
// Resulting polarization in the z direction
c1->cd(2);
t1->Draw("save_polz >> pol_z(128, -1, 1)","save_detID==902&&save_particleID==-13"); //-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();
// Initial beam spot
c1->cd(3)->SetGridx();c1->cd(3)->SetGridy(); // c2->cd(3);
t1->Draw("muIniPosZ[0]:muIniPosY[0]>>ini_hist(64,-40,40,64,-1710,-1646)","muIniPolZ > -2");
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);
// Resulting beam spot
c1->cd(4)->SetGridx();c1->cd(4)->SetGridy(); // c1->cd(3);
// t1->Draw("save_y:save_x>>bspot_hist(64,-40,40,64,-40,40)","save_detID==902&&save_particleID==-13");
t1->Draw("save_y:save_x>>bspot_hist(64,-40,40,64,-40,40)","save_detID==902&&save_particleID==-13&&sqrt(save_x**2+save_y**2)<35");
bspot_hist->SetTitle("Muon beam cross section [mm]");
/* bspot_hist->GetXaxis()->SetTitle("x at z=16 [mm]")
bspot_hist->GetYaxis()->SetTitle("y at z=16 [mm]")*/
bspot_hist->Draw();
bspot_hist->Draw("cont0 same");
bspot_hist->ProjectionX("xfsect",31,33);
bspot_hist->ProjectionY("yfsect",31,33);
// Initial beam spot cross section
c1->cd(5);
xisect->SetTitle("x initial cross section [mm]");
xisect->SetFillStyle(1001);
xisect->SetFillColor(45);
xisect->Draw();
// Resulting beam spot cross section
c1->cd(6);
/* c1->cd(6)->SetGridx();c1->cd(6)->SetGridy();
t1->Draw("save_y:save_x>>bspot_hist1(64,-40,40,64,-40,40)","save_detID==901&&save_particleID==-13&&sqrt(save_x**2+save_y**2)<=15");
bspot_hist1->SetTitle("Muon beam cross section [mm]");
bspot_hist1->GetXaxis()->SetTitle("x at z=0 (mm)")
bspot_hist1->GetYaxis()->SetTitle("y at z=0 (mm)")
bspot_hist1->Draw();
bspot_hist1->Draw("cont0 same");*/
xfsect->SetTitle("x cross section [mm]");
xfsect->SetFillStyle(1001);
xfsect->SetFillColor(45);
xfsect->Draw();
}
void phase_space(char* fname)
{
TFile* f1 = new TFile(fname);
gStyle->SetPalette(1,0);
TCanvas* c1 = new TCanvas("c1","Phase Space",60,40,600,800);
c1->Divide(2,2);
c1->SetGridx();
c1->SetGridy();
// Resulting phase space px,py
c1->cd(1);
t1->Draw("muIniMomY[0]:muIniMomX[0]>>ini_hist(64,-40,40,64,-40,40)","muIniPolX > -2");
ini_hist->SetTitle("Muon initial beam phase space []");
ini_hist->Draw();
ini_hist->Draw("cont0 same");
ini_hist->ProjectionX("xisect",31,33);
ini_hist->ProjectionY("yisect",31,33);
// Resulting polarization in the z direction
c1->cd(2);
t1->Draw("muTargetPosY[0]:muTargetPosX[0]>>bspot_hist(64,-40,40,64,-40,40)","muTargetPolX > -2");
bspot_hist->SetTitle("Muon beam cross section [mm]");
bspot_hist->Draw();
bspot_hist->Draw("cont0 same");
bspot_hist->ProjectionX("xfsect",31,33);
bspot_hist->ProjectionY("yfsect",31,33);
// Initial beam spot
c1->cd(3)->SetGridx();c1->cd(3)->SetGridy(); // c2->cd(3);
t1->Draw("muIniPosY[0]:muIniPosX[0]>>ini_hist(64,-40,40,64,-40,40)","muIniPolX > -2");
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);
// Resulting beam spot
c1->cd(4)->SetGridx();c1->cd(4)->SetGridy(); // c1->cd(3);
t1->Draw("muTargetPosY[0]:muTargetPosX[0]>>bspot_hist(64,-40,40,64,-40,40)","muTargetPolX > -2");
bspot_hist->SetTitle("Muon beam cross section [mm]");
bspot_hist->Draw();
bspot_hist->Draw("cont0 same");
bspot_hist->ProjectionX("xfsect",31,33);
bspot_hist->ProjectionY("yfsect",31,33);
}
void Slice(char* fname, Int_t Slice)
{
TFile* f1 = new TFile(fname);
TString Cond,Title;
Double_t zSlice;
gStyle->SetPalette(1,0);
TCanvas* c1 = new TCanvas("c1","Polarization and Transmission",60,40,600,800);
c1->Divide(2,3);
c1->SetGridx();
c1->SetGridy();
// Resulting polarization in the x direction
c1->cd(1);
Cond.Form("save_detID== %04d&&save_particleID==-13",Slice);
t1->Draw("save_polx >> pol_x(64,-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");
// Resulting polarization in the z direction
c1->cd(2);
t1->Draw("save_polz >> pol_z(64,-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();
// Initial beam spot
c1->cd(3)->SetGridx();c1->cd(3)->SetGridy(); // c2->cd(3);
t1->Draw("muIniPosZ[0]:muIniPosY[0]>>ini_hist(64,-40,40,64,-1710,-1646)","muIniPolZ > -2");
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);
// Resulting beam spot
c1->cd(4)->SetGridx();c1->cd(4)->SetGridy(); // c1->cd(3);
t1->Draw("save_z>>zhist",Cond);
zSlice=zhist->GetMean(1);
t1->Draw("save_y:save_x>>bspot_hist(64,-40,40,64,-40,40)",Cond);
Title.Form("Muon beam cross section [mm], z=%g [mm]",zSlice);
bspot_hist->SetTitle(Title);
bspot_hist->Draw();
bspot_hist->Draw("cont0 same");
bspot_hist->ProjectionX("xfsect",31,33);
bspot_hist->ProjectionY("yfsect",31,33);
// Initial beam spot cross section
c1->cd(5);
xisect->SetTitle("x initial cross section [mm]");
xisect->SetFillStyle(1001);
xisect->SetFillColor(45);
xisect->Draw();
// Resulting beam spot cross section
c1->cd(6);
xfsect->SetTitle("x cross section [mm]");
xfsect->SetFillStyle(1001);
xfsect->SetFillColor(45);
xfsect->Draw();
}
void VarSlice(char* fname, char* Var, Int_t Slice, Int_t SubPlot)
{
TFile* f1 = new TFile(fname);
TString Cond,Title,Histo;
Double_t zSlice;
gStyle->SetPalette(1,0);
// TCanvas* c1 = new TCanvas("c1",Var);
// c1->Divide(3,5);
c1->cd(SubPlot);
Cond.Form("save_detID== %04d&&save_particleID==-13",Slice);
Histo.Form("%s >> histo",Var);
t1->Draw(Histo,Cond);
histo->SetTitle(Var);
histo->Draw("F");
histo->SetFillStyle(1001);
histo->SetFillColor(kGreen-5);
histo->Draw();
}
void VarxySlice(char* fname, char* Varx, char* Vary, Int_t Slice, Int_t SubPlot)
{
TFile* f1 = new TFile(fname);
TString Cond,Title,Histo;
Double_t zSlice;
gStyle->SetPalette(1,0);
// TCanvas* c1 = new TCanvas("c1",Var);
// c1->Divide(3,5);
c1->cd(SubPlot);
Cond.Form("save_detID== %04d&&save_particleID==-13",Slice);
Histo.Form("%s:%s >> histo",Vary,Varx);
t1->Draw(Histo,Cond);
histo->SetTitle(Histo);
histo->Draw();
histo->Draw("cont0 same");
}
void Beam_Env(char* fname, Int_t iSlice, Int_t fSlice)
{
TFile* f1 = new TFile(fname);
TString Cond,Title;
Double_t zSlice,Xrms,Yrms,Xavg,Yavg,Nbins;
Int_t i;
gStyle->SetPalette(1,0);
TCanvas* c1 = new TCanvas("c1","Polarization and Transmission");
// Initial beam spot
// t1->Draw("muIniPosZ[0]:muIniPosY[0]:muIniPosX[0]>>ini_hist(64,-40,40,64,-40,40)","muIniPolX > -2");
// 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);
printf("SliceID \t Z \t Xrms \t Yrms \t Xavg \t Yavg \t N\n");
for (i=iSlice; i<=fSlice; i++)
{
// Resulting beam spot
Cond.Form("save_detID==%d&&save_particleID==-13",i);
// Cond.Form("save_detID==%d&& (save_particleID!=22) ",i);
t1->Draw("save_z:save_y:save_x>>bspot_hist(64,-40,40,64,-40,40,100,-2600,300)",Cond);
// t1->Draw("save_z:save_y:save_x>>bspot_hist(64,-40,40,64,-40,40,100,-290,290)",Cond);
// Get properties of slice
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);
// delete bspot_hist;
}
// t1->Draw("save_z:save_y:save_x>>bspot_hist(64,-40,40,64,-40,40,100,-290,290)","save_particleID==-13");
// t1->Draw("save_z:save_y:save_x>>bspot_hist(64,-40,40,64,-40,40,100,-2600,300)","save_particleID==-13");
// bspot_hist->Draw("Iso");
}

301
run/macros/sr_tr_cp.C Normal file
View File

@ -0,0 +1,301 @@
// To run root without logon screen and script use: root -l -n
// To run a named macro in root use:
//.L sr_tr.C
// sr_transm("data/sr1_1100.root")
#include <string>
void sr_transm(char* fname)
{
TFile* f1 = new TFile(fname);
gStyle->SetPalette(1,0);
TCanvas* c1 = new TCanvas("c1","Polarization and Transmission",60,40,600,800);
c1->Divide(2,3);
c1->SetGridx();
c1->SetGridy();
// Resulting polarization in the x direction
c1->cd(1);
t1->Draw("save_polx >> pol_x(128, -1, 1)","save_detID==902&&save_particleID==-13"); //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");
// Resulting polarization in the z direction
c1->cd(2);
t1->Draw("save_polz >> pol_z(128, -1, 1)","save_detID==902&&save_particleID==-13"); //-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();
// Initial beam spot
c1->cd(3)->SetGridx();c1->cd(3)->SetGridy(); // c2->cd(3);
t1->Draw("muIniPosZ[0]:muIniPosY[0]>>ini_hist(64,-40,40,64,-1710,-1646)","muIniPolZ > -2");
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);
// Resulting beam spot
c1->cd(4)->SetGridx();c1->cd(4)->SetGridy(); // c1->cd(3);
// t1->Draw("save_y:save_x>>bspot_hist(64,-40,40,64,-40,40)","save_detID==902&&save_particleID==-13");
t1->Draw("save_y:save_x>>bspot_hist(64,-40,40,64,-40,40)","save_detID==902&&save_particleID==-13&&sqrt(save_x**2+save_y**2)<20");
bspot_hist->SetTitle("Muon beam cross section [mm]");
bspot_hist->Draw();
bspot_hist->Draw("cont0 same");
bspot_hist->ProjectionX("xfsect",31,33);
bspot_hist->ProjectionY("yfsect",31,33);
// Initial beam spot cross section
c1->cd(5);
xisect->SetTitle("x initial cross section [mm]");
xisect->SetFillStyle(1001);
xisect->SetFillColor(45);
xisect->Draw();
// Resulting beam spot cross section
c1->cd(6);
xfsect->SetTitle("x cross section [mm]");
xfsect->SetFillStyle(1001);
xfsect->SetFillColor(45);
xfsect->Draw();
}
void phase_space(char* fname)
{
TFile* f1 = new TFile(fname);
gStyle->SetPalette(1,0);
TCanvas* c1 = new TCanvas("c1","Phase Space",60,40,600,800);
c1->Divide(2,2);
c1->SetGridx();
c1->SetGridy();
// Resulting phase space px,py
c1->cd(1);
t1->Draw("muIniMomY[0]:muIniMomX[0]>>ini_hist(64,-40,40,64,-40,40)","muIniPolX > -2");
ini_hist->SetTitle("Muon initial beam phase space []");
ini_hist->Draw();
ini_hist->Draw("cont0 same");
ini_hist->ProjectionX("xisect",31,33);
ini_hist->ProjectionY("yisect",31,33);
// Resulting polarization in the z direction
c1->cd(2);
t1->Draw("muTargetPosY[0]:muTargetPosX[0]>>bspot_hist(64,-40,40,64,-40,40)","muTargetPolX > -2");
bspot_hist->SetTitle("Muon beam cross section [mm]");
bspot_hist->Draw();
bspot_hist->Draw("cont0 same");
bspot_hist->ProjectionX("xfsect",31,33);
bspot_hist->ProjectionY("yfsect",31,33);
// Initial beam spot
c1->cd(3)->SetGridx();c1->cd(3)->SetGridy(); // c2->cd(3);
t1->Draw("muIniPosY[0]:muIniPosX[0]>>ini_hist(64,-40,40,64,-40,40)","muIniPolX > -2");
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);
// Resulting beam spot
c1->cd(4)->SetGridx();c1->cd(4)->SetGridy(); // c1->cd(3);
t1->Draw("muTargetPosY[0]:muTargetPosX[0]>>bspot_hist(64,-40,40,64,-40,40)","muTargetPolX > -2");
bspot_hist->SetTitle("Muon beam cross section [mm]");
bspot_hist->Draw();
bspot_hist->Draw("cont0 same");
bspot_hist->ProjectionX("xfsect",31,33);
bspot_hist->ProjectionY("yfsect",31,33);
}
void Slice(char* fname, Int_t Slice)
{
TFile* f1 = new TFile(fname);
TString Cond,Title;
Double_t zSlice;
gStyle->SetPalette(1,0);
TCanvas* c1 = new TCanvas("c1","Polarization and Transmission",60,40,600,800);
c1->Divide(2,3);
c1->SetGridx();
c1->SetGridy();
// Resulting polarization in the x direction
c1->cd(1);
Cond.Form("save_detID== %04d&&save_particleID==-13",Slice);
t1->Draw("save_polx >> pol_x(64,-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");
// Resulting polarization in the z direction
c1->cd(2);
t1->Draw("save_polz >> pol_z(64,-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();
// Initial beam spot
c1->cd(3)->SetGridx();c1->cd(3)->SetGridy(); // c2->cd(3);
t1->Draw("muIniPosZ[0]:muIniPosY[0]>>ini_hist(64,-40,40,64,-1710,-1646)","muIniPolZ > -2");
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);
// Resulting beam spot
c1->cd(4)->SetGridx();c1->cd(4)->SetGridy(); // c1->cd(3);
t1->Draw("save_z>>zhist",Cond);
zSlice=zhist->GetMean(1);
t1->Draw("save_y:save_x>>bspot_hist(64,-40,40,64,-40,40)",Cond);
Title.Form("Muon beam cross section [mm], z=%g [mm]",zSlice);
bspot_hist->SetTitle(Title);
bspot_hist->Draw();
bspot_hist->Draw("cont0 same");
bspot_hist->ProjectionX("xfsect",31,33);
bspot_hist->ProjectionY("yfsect",31,33);
// Initial beam spot cross section
c1->cd(5);
xisect->SetTitle("x initial cross section [mm]");
xisect->SetFillStyle(1001);
xisect->SetFillColor(45);
xisect->Draw();
// Resulting beam spot cross section
c1->cd(6);
xfsect->SetTitle("x cross section [mm]");
xfsect->SetFillStyle(1001);
xfsect->SetFillColor(45);
xfsect->Draw();
}
void VarSlice(char* fname, char* Var, Int_t Slice, Int_t SubPlot)
{
TFile* f1 = new TFile(fname);
TString Cond,Title,Histo;
Double_t zSlice;
gStyle->SetPalette(1,0);
// TCanvas* c1 = new TCanvas("c1",Var);
// c1->Divide(3,5);
c1->cd(SubPlot);
Cond.Form("save_detID== %04d&&save_particleID==-13",Slice);
Histo.Form("%s >> histo",Var);
t1->Draw(Histo,Cond);
histo->SetTitle(Var);
histo->Draw("F");
histo->SetFillStyle(1001);
histo->SetFillColor(kGreen-5);
histo->Draw();
}
void VarxySlice(char* fname, char* Varx, char* Vary, Int_t Slice, Int_t SubPlot)
{
TFile* f1 = new TFile(fname);
TString Cond,Title,Histo;
Double_t zSlice;
gStyle->SetPalette(1,0);
// TCanvas* c1 = new TCanvas("c1",Var);
// c1->Divide(3,5);
c1->cd(SubPlot);
Cond.Form("save_detID== %04d&&save_particleID==-13",Slice);
Histo.Form("%s:%s >> histo",Vary,Varx);
t1->Draw(Histo,Cond);
histo->SetTitle(Histo);
histo->Draw();
histo->Draw("cont0 same");
}
void Beam_Env(char* fname, Int_t iSlice, Int_t fSlice)
{
TFile* f1 = new TFile(fname);
TString Cond,Title;
Double_t zSlice,Xrms,Yrms,Xavg,Yavg,Nbins;
Int_t i;
gStyle->SetPalette(1,0);
TCanvas* c1 = new TCanvas("c1","Polarization and Transmission");
// Initial beam spot
// t1->Draw("muIniPosZ[0]:muIniPosY[0]:muIniPosX[0]>>ini_hist(64,-40,40,64,-40,40)","muIniPolX > -2");
// 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);
printf("SliceID \t Z \t Xrms \t Yrms \t Xavg \t Yavg \t N\n");
for (i=iSlice; i<=fSlice; i++)
{
// Resulting beam spot
Cond.Form("save_detID==%d&&save_particleID==-13",i);
// Cond.Form("save_detID==%d&& (save_particleID!=22) ",i);
t1->Draw("save_z:save_y:save_x>>bspot_hist(64,-40,40,64,-40,40,100,-2600,300)",Cond);
// t1->Draw("save_z:save_y:save_x>>bspot_hist(64,-40,40,64,-40,40,100,-290,290)",Cond);
// Get properties of slice
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);
// delete bspot_hist;
}
// t1->Draw("save_z:save_y:save_x>>bspot_hist(64,-40,40,64,-40,40,100,-290,290)","save_particleID==-13");
// t1->Draw("save_z:save_y:save_x>>bspot_hist(64,-40,40,64,-40,40,100,-2600,300)","save_particleID==-13");
// bspot_hist->Draw("Iso");
}

197
run/macros/sr_transm.C Normal file
View File

@ -0,0 +1,197 @@
// To run root without logon screen and script use: root -l -n
// To run a named macro in root use:
//.L sr_transm.C
// sr_transm("data/sr1_1100.root")
#include <string>
void sr_transm(char* fname)
{
TFile* f1 = new TFile(fname);
//TCanvas* c1 = new TCanvas("c1","Polarization and Transmission",60,40,800,800);
gStyle->SetPalette(1,0);
TCanvas* c1 = new TCanvas("c1","Polarization and Transmission",60,40,800,800);
c1->Divide(2,2);
c1->SetGridx();
c1->SetGridy();
//c1->cd(1)->SetGridx();c1->cd(1)->SetGridy();
c1->cd(4)->SetGridx();c1->cd(4)->SetGridy(); // c1->cd(3);
t1->Draw("muTargetPosY[0]:muTargetPosX[0]>>bspot_hist(64,-40,40,64,-40,40)","muTargetPolX > -2");
//t1->Draw("det_y:det_x>>bspot_hist(64,-80,80,64,-80,80)"); //t1->Draw("mcpHit.positiony*10:mcpHit.positionx*10>>bspot_hist(64,-40,40,64,-40,40)","mcpHit.positionz==2&&mcpHit.ID==0");
bspot_hist->SetTitle("Muon beam cross section [mm]");
//bspot_hist->Draw("cont0");
bspot_hist->Draw();
bspot_hist->Draw("cont0 same");
// t1->Draw("mcpHit.positionx*10:mcpHit.positiony*10>>bspot_hist(128,-12,12,128,-12,12)","","goff");
bspot_hist->ProjectionX("xsect",31,33);
bspot_hist->ProjectionY("ysect",31,33);
// TText *pt = new TText(-36., 33., "L1 = 11.0 kV");
// pt->SetTextSize(0.05);
// pt->Draw();
//pl=new TPaveLabel(2, 155, 16, 175, "RA_Up = 11 kV");
//pl=new TPaveLabel(-75, 60, -15, 75, "#splitline{RA_Up = 8.6 kV}{L3 = 6.78 kV}");
//pl->SetBorderSize(1);
//pl->SetTextSize(0.43);
//pl->Draw();
c1->cd(3);
xsect->SetTitle("x cross section [mm]");
xsect->SetFillStyle(1001);
//xsect->Draw("F");
xsect->SetFillColor(45);
xsect->Draw();
c1->cd(2);
t1->Draw("muTargetPolZ >> pol_z(128, -1, -0.95)","(muTargetPolZ > -2)"); //-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();
c1->cd(1);
//Energy loss in CFoil
//t1->Draw("10-det_edep_mup*1000","(det_n==1)&&(det_edep<0.1)&&(det_edep_mup*1000>0.1)");
t1->Draw("muTargetPolX >> pol_x(128, -0.4, 0)","(muTargetPolX > -2)"); //0, 2
//t1->Draw("det_edep_mup*1000","(det_n==1)&&(det_edep<0.1)");
//t1->Draw("det_edep_mup*1000>>tmp_hist(64,-80,80)","(det_n==1)&&(det_edep<0.1)");
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(kBlue-5);
pol_x->SetFillColor(kGreen-5);
pol_x->Draw();
/*
TAxis* yax = ang_hist->GetYaxis();
yax->SetNdivisions(206); // N2*100 + N1 (N2 = secondary, N1 = primary) //yax->SetNdivisions(410);
yax->SetTickLength(0.018);
yax->SetLabelSize(0.04);
yax->SetTitleSize(0.045);
yax->SetTitleOffset(1.3);
yax->SetTitle("Counts (arb. units)");
*/
c1->SaveAs("data/sr_bspot_0.eps");
}
/*
/// Plot of muon beam angular dispersion after lens 1:
TCanvas* c2 = new TCanvas("c2","Angular Dispersion",60,40,800,800);
//gStyle->SetPalette(1,0);
//c1->Divide(2,2);
c2->SetGridx();
c2->SetGridy();
//c2->SetFillColor(kWhite); // Fill background with white for printing
//c2->cd(1);
t1->Draw("save_pz","(save_detID==800)"); //xsect->SetTitle("x cross section [mm]");
htemp->SetTitle("Position 1 [mm]");
htemp->Draw("F");
htemp->SetFillStyle(1001);
htemp->SetFillColor(38);
htemp->Draw();
//t1->Draw("atan2(mcpHit.posyini,mcpHit.posxini)*57.29577951308:atan2(mcpHit.positiony,mcpHit.positionx)*57.29577951308>
// Double_t Pi = TMath::Pi(); // 180/Pi
//c2->cd(2)->SetGridx(); c1->cd(2)->SetGridy();
//t1->Draw("atan2(save_py,save_px)*57.29577951308:atan2(mcpHit.positiony,mcpHit.positionx)*57.29577951308>>bspot_hist3(64,-180,180,64,-180,180)","mcpHit.ID==0");
///t1->Draw("atan2(save_py,save_px)*57.29577951308");
//t1->Draw("det_y[0]:det_x[0]>>bspot_hist(64,-40,40,64,-40,40)");
t1->Draw("90-atan2(save_pz,sqrt(save_px^2 + save_py^2))*57.29577951308>>ang_hist(64,0,5)","(save_detID==900)");
//t1->Draw("90-atan2(save_pz,sqrt(save_px^2 + save_py^2))*57.29577951308");
//htemp->SetTitle("Correlation angle for L1 (8.7 kV, 15 keV); Final angle [deg.]; Inital angle [deg.]");
ang_hist->GetXaxis()->SetRangeUser(-5., 5.);
ang_hist->Draw("F");
ang_hist->SetFillStyle(1001);
ang_hist->SetFillColor(805);
ang_hist->SetTitleSize(0.025);
ang_hist->SetTitle("Muon angular distribution after L1");
ang_hist->GetXaxis()->SetNdivisions(206);
ang_hist->GetXaxis()->SetTickLength(0.018);
ang_hist->GetXaxis()->SetLabelSize(0.04);
ang_hist->GetXaxis()->SetTitleSize(0.045);
ang_hist->GetYaxis()->SetTitleOffset(1.0);
ang_hist->GetXaxis()->SetTitle("Angle (deg.)");
TAxis* yax = ang_hist->GetYaxis();
yax->SetNdivisions(206); // N2*100 + N1 (N2 = secondary, N1 = primary) //yax->SetNdivisions(410);
yax->SetTickLength(0.018);
yax->SetLabelSize(0.04);
yax->SetTitleSize(0.045);
yax->SetTitleOffset(1.3);
yax->SetTitle("Counts (arb. units)");
ang_hist->Draw();
TText *pt = new TText(3.2, 200., "Plane at -300 mm");
pt->SetTextSize(0.03);
pt->Draw();
//TPaveText *pt = new TPaveText(0.6,0.4,0.89, 0.45,"trNDC");
//pt->SetTextSize(0.03);
//pt->SetFillColor(0); //pt->SetFillColor(390);
//pt->SetBorderSize(0);
//pt->SetTextAlign(12);
//pte = pt->AddText("Plane at z = -300 mm");
//pt->Draw();
//c2->SaveAs("data/sr1_ang_dist_300.eps");
}
*/
/// Insert simple text directly from command line
// TText text(3, 300, "Plane at -380 mm");
// text.SetTextSize(0.03);
// text.Draw();
/// Optional changes of statistics
// TPaveStats* sb2=(TPaveStats *)(ang_hist->GetListOfFunctions()->FindObject("stats"));
// sb2->SetTextSize(0.025);
// sb2->SetFillColor(390);
/// An easy way to plot functions
// (new TF1("fun1","(1-x)*sqrt(1-x*x)",-1,1))->Draw();
// (new TF2("fun2","(175/80)**2*(1-y)*(1-x**2)+(1+y)*(1-x)**2",-1,1,-1,1))->Draw("lego")
// From: www-pnp.physics.ox.ac.uk/~west/root/plotting_with_style.html#demo_background
/// For publishing (insert these lines in rootlogon.C):
// gStyle->SetLineWidth(2.);
// gStyle->SetTextSize(1.1);
// gStyle->SetLabelSize(0.06,"xy");
// gStyle->SetTitleSize(0.06,"xy");
// gStyle->SetTitleOffset(1.2,"x");
// gStyle->SetTitleOffset(1.0,"y");
// gStyle->SetPadTopMargin(0.1);
// gStyle->SetPadRightMargin(0.1);
// gStyle->SetPadBottomMargin(0.16);
// gStyle->SetPadLeftMargin(0.12);

96
run/macros/td_mcp2_tof.C Normal file
View File

@ -0,0 +1,96 @@
// To run root without logon screen and script use: root -l -n
// To run a named macro in root use:
//.L td_mcp2_tof.C
// td_mcp2_tof("data/sr1_1100.root")
#include <string>
#define NSIZE 100000
void td_mcp2_tof(char* fname)
{
Double_t *ptofTD, *ptofMCP2, *pEventIDTD, *pEventIDMCP2;
Double_t tofTD[NSIZE], tofMCP2[NSIZE], EventIDTD[NSIZE], EventIDMCP2[NSIZE];
Long64_t i, j, k, nEventsTD, nEventsMCP2;
TFile* f1 = new TFile(fname);
// muon time-of-flight to TD
t1->Draw("eventID:1000*save_time","save_detID==890&&save_particleID==-13"); //0, 2
nEventsTD = t1->GetSelectedRows();
ptofTD = t1->GetV2();
pEventIDTD = t1->GetV1();
for (i = 0; i<nEventsTD; i++){
EventIDTD[i] = pEventIDTD[i];
tofTD[i] = ptofTD[i];
}
// muon time-of-flight to MCP2, mu+
t1->Draw("eventID:1000*save_time","save_detID==902&&save_particleID==-13"); //0, 2
nEventsMCP2 = t1->GetSelectedRows();
ptofMCP2 = t1->GetV2();
pEventIDMCP2 = t1->GetV1();
for (i = 0; i<nEventsMCP2; i++){
EventIDMCP2[i] = pEventIDMCP2[i];
tofMCP2[i] = ptofMCP2[i];
}
// TOF TD-MCP2
TH1D *htof = new TH1D("Time-of-flight TD-MCP2", "Time-of-flight TD-MCP2", 1000, 0.25, 500.25);
k = 0;
for (i = 0; i<nEventsMCP2; i++){
for (j = k; j<nEventsTD; j++){
// cout << i << ", " << j << ", " << EventIDMCP2[i] << ", " << EventIDTD[j] << ", " << tofMCP2[i] << ", " << tofTD[j] << endl;
if (EventIDMCP2[i] == EventIDTD[j]){
htof->Fill(tofMCP2[i] - tofTD[j]);
k = j + 1;
break;
}
}
}
gStyle->SetPalette(1,0);
// gStyle->SetName(fname);
gStyle->SetOptTitle(1);
TCanvas* c1 = new TCanvas("Time of Flight","Time of Flight",60,40,600,800);
c1->Divide(1,3);
c1->SetGridx();
c1->SetGridy();
// muon time-of-flight to TD
c1->cd(1);
t1->Draw("1000*save_time >> tof_td(2000, 0.25, 1000.25)","save_detID==890&&save_particleID==-13"); //0, 2
tof_td->SetTitle("Time-of-flight to TD");
tof_td->GetXaxis()->SetTitle("TD time of flight (ns)");
// tof_td->GetXaxis()->SetNdivisions(405);
tof_td->GetYaxis()->SetNdivisions(406);
// tof_td->GetXaxis()->SetTickLength(0.018);
tof_td->GetYaxis()->SetTickLength(0.018);
tof_td->Draw("F");
tof_td->SetFillStyle(1001);
tof_td->SetFillColor(kGreen-5);
tof_td->Draw();
//c1->SaveAs("data/sr_bspot_0.eps");
// muon time-of-flight to MCP2
c1->cd(2);
t1->Draw("1000*save_time >> tof_mcp2(2000, 0.25, 1000.25)","save_detID==902&&((save_particleID==-13)||(save_particleID==-1313))"); //0, 2
tof_mcp2->SetTitle("Time-of-flight to MCP2");
tof_mcp2->GetXaxis()->SetTitle("MCP2 time of flight (ns)");
// tof_mcp2->GetXaxis()->SetNdivisions(405);
tof_mcp2->GetYaxis()->SetNdivisions(406);
// tof_mcp2->GetXaxis()->SetTickLength(0.018);
tof_mcp2->GetYaxis()->SetTickLength(0.018);
tof_mcp2->Draw("F");
tof_mcp2->SetFillStyle(1001);
tof_mcp2->Draw();
c1->cd(3);
htof->SetTitle("Time-of-flight TD-MCP2");
htof->GetXaxis()->SetTitle("TD-MCP2 time of flight (ns)");
htof->GetYaxis()->SetNdivisions(406);
htof->GetYaxis()->SetTickLength(0.018);
htof->Draw("F");
htof->Draw();
}

9
run/macros/tmp.C Normal file
View File

@ -0,0 +1,9 @@
// To run: root -b LoopFiles.C
{
TString fileDir,fileName,file;
Int_t fileNumber,L1,L3;
gROOT->ProcessLine(".L BeamEnv.C");
BeamSlices("data/musr_142180.root",851,904);
}

70
run/macros/treeviewer.C Normal file
View File

@ -0,0 +1,70 @@
void treeviewer() {
//=========Macro generated by ROOT version5.32/00
//=========for tree "t1" (Sat May 5 17:31:13 2012)
//===This macro can be opened from a TreeViewer session after loading
//===the corresponding tree, or by running root with the macro name argument
open_session();
}
open_session(void *p = 0) {
gSystem->Load("libTreeViewer");
TTreeViewer *treeview = (TTreeViewer *) p;
if (!treeview) treeview = new TTreeViewer();
TTree *tv_tree = (TTree*)gROOT->FindObject("t1");
TFile *tv_file = (TFile*)gROOT->GetListOfFiles()->FindObject("data/musr_50703.root");
if (!tv_tree) {
if (!tv_file) tv_file = new TFile("data/musr_50703.root");
if (tv_file) tv_tree = (TTree*)tv_file->Get("t1");
if(!tv_tree) {
printf("Tree %s not found", fTree->GetName());
return;
}
}
treeview->SetTreeName("t1");
treeview->SetNexpressions(10);
// Set expressions on axis and cut
TTVLVEntry *item;
// X expression
item = treeview->ExpressionItem(0);
item->SetExpression("save_z", "~save_z");
// Y expression
item = treeview->ExpressionItem(1);
item->SetExpression("save_ke", "~save_ke");
// Z expression
item = treeview->ExpressionItem(2);
item->SetExpression("", "-empty-");
// Cut expression
item = treeview->ExpressionItem(3);
item->SetExpression("save_detID>870&&save_particleID==-13", "~ScanAlongZ");
// Scan list
item = treeview->ExpressionItem(4);
item->SetExpression("", "Scan box");
// User defined expressions
item = treeview->ExpressionItem(5);
item->SetExpression("save_detID==902&&save_particleID==-13", "~SamplePosition", kTRUE);
item = treeview->ExpressionItem(6);
item->SetExpression("save_detID>870&&save_particleID==-13", "~ScanAlongZ", kTRUE);
item = treeview->ExpressionItem(7);
item->SetExpression("save_detID==890&&save_particleID==-13", "~TDAfterFoil", kTRUE);
item = treeview->ExpressionItem(8);
item->SetExpression("acos(save_pz/sqrt(save_px^2+save_py^2+save_pz^2))*180/3.14", "~AcosTheta", kFALSE);
item = treeview->ExpressionItem(9);
item->SetExpression("save_time-0.3152", "~TOFTDSample", kFALSE);
item = treeview->ExpressionItem(10);
item->SetExpression("", "-empty-", kFALSE);
item = treeview->ExpressionItem(11);
item->SetExpression("", "-empty-", kFALSE);
item = treeview->ExpressionItem(12);
item->SetExpression("", "-empty-", kFALSE);
item = treeview->ExpressionItem(13);
item->SetExpression("", "-empty-", kFALSE);
item = treeview->ExpressionItem(14);
item->SetExpression("", "-empty-", kFALSE);
//--- session object
tv_session = new TTVSession(treeview);
treeview->SetSession(tv_session);
//--- Connect first record
tv_session->First();
}

11
run/serialrun.sh Executable file
View File

@ -0,0 +1,11 @@
#!/bin/bash
# i=730044
# ../bin/Linux-g++/musrSim $i.mac
# i=730050
# ../bin/Linux-g++/musrSim $i.mac
i=720004
while [ "$i" -lt 720016 ]; do
../bin/Linux-g++/musrSim $i.mac
i=$(expr $i + 1)
done