musrsim/run/SiPixel/SiPixel.C

343 lines
11 KiB
C

// To run root without logon screen and script use: root -l -n
// To run a named macro in root use:
//.L SiPixel_tr.C
// muSpot("data/musr_1100.root", 868)
#include <string>
void muSpot(const char* fname, int saveID)
{
TString cutStr;
TH1D *pol_x, *pol_z, *xisect, *yisect, *xfsect, *yfsect;
TH2D *ini_hist, *bspot_hist, *bspot_hist1;
TTree *t1;
TFile* f1 = new TFile(fname);
t1 = (TTree*) gDirectory->FindObjectAny("t1");
// Get muon (particleID == -13) hits in saveID
cutStr.Form("save_detID== %3i &&save_particleID==-13",saveID);
// Get positron (particleID == 22) hits in saveID
//sprintf(poscutStr,"save_detID== %3i &&save_particleID==22",saveID);
gStyle->SetPalette(1,0);
TCanvas* c1 = new TCanvas("c1","Beam Spot",60,60,600,600);
c1->Divide(2,2);
c1->SetGridx();
c1->SetGridy();
// Initial beam spot
c1->cd(1)->SetGridx();c1->cd(1)->SetGridy();
t1->Draw("muIniPosY[0]:muIniPosX[0]>>ini_hist(60,-20, 20,60,-20,20)");
ini_hist = (TH2D*) gDirectory->Get("ini_hist");
ini_hist->SetTitle("Muon initial beam cross section [mm]");
ini_hist->Draw();
ini_hist->Draw("cont0 same");
ini_hist->ProjectionX("xisect");
ini_hist->ProjectionY("yisect");
// Resulting beam spot
c1->cd(2)->SetGridx();c1->cd(2)->SetGridy();
// Muons decay position (on target)
//t1->Draw("muDecayPosX:muDecayPosY>>bspot_hist(60,-20,20,60,-20,20)");
// Muons passing through save_detID
t1->Draw("save_y:save_x>>bspot_hist(60,-20,20,60,-20,20)",cutStr);
bspot_hist = (TH2D*) gDirectory->Get("bspot_hist");
bspot_hist->SetTitle("Muon beam cross section [mm]");
bspot_hist->Draw();
bspot_hist->Draw("cont0 same");
bspot_hist->ProjectionX("xfsect");
bspot_hist->ProjectionY("yfsect");
// Initial beam spot cross section
c1->cd(3);
c1->cd(3)->SetGridx();c1->cd(3)->SetGridy();
xisect = (TH1D*) gDirectory->Get("xisect");
xisect->SetTitle("x initial cross section [mm]");
xisect->SetFillStyle(1001);
xisect->SetFillColor(45);
xisect->Draw();
// Resulting beam spot cross section
c1->cd(4);
c1->cd(4)->SetGridx();c1->cd(4)->SetGridy();
xfsect = (TH1D*) gDirectory->Get("xfsect");
xfsect->SetTitle("x muon cross section [mm]");
xfsect->SetFillStyle(1001);
xfsect->SetFillColor(45);
xfsect->Draw();
// yfsect = (TH1D*) gDirectory->Get("yfsect");
// yfsect->SetTitle("x muon cross section [mm]");
// yfsect->SetFillStyle(1001);
// yfsect->SetFillColor(46);
// yfsect->Draw();
// bspot_hist1 = (TH2D*) gDirectory->Get("bspot_hist1");
// bspot_hist1->SetTitle("Muon beam cross section [mm]");
// bspot_hist1->GetXaxis()->SetTitle("x (mm)");
// bspot_hist1->GetYaxis()->SetTitle("y (mm)");
// bspot_hist1->Draw();
// bspot_hist1->Draw("cont0 same");
// cout << "Mean x, Mean y, RMS x, RMS y: " << bspot_hist1->GetMean(1) << ", " << bspot_hist1->GetMean(2) << ", " << bspot_hist1->GetRMS(1) << ", " << bspot_hist1->GetRMS(2) << ", " << endl;
// cout << "Fraction in [-10.,10.]: " << bspot_hist1->Integral(21,60,21,60)/bspot_hist1->Integral(1,80,1,80) << endl;
// cout << "Fraction in [-12.5,12.5]: " << bspot_hist1->Integral(16,65,16,65)/bspot_hist1->Integral(1,80,1,80) << endl;
// cout << "Fraction in [-5.0,5.0]: " << bspot_hist1->Integral(31,50,31,50)/bspot_hist1->Integral(1,80,1,80) << endl;
// cout << "Fraction in [-2.5,2.5]: " << bspot_hist1->Integral(36,45,36,45)/bspot_hist1->Integral(1,80,1,80) << endl;
// cout << "N[-20,20]: " << bspot_hist1->Integral(1,80,1,80) << endl;
}
void phase_space(const char* fname)
{
TH1D *pol_x, *pol_z, *xisect, *yisect, *xfsect, *yfsect;
TH2D *ini_hist, *bspot_hist, *bspot_hist1;
TTree *t1;
TFile* f1 = new TFile(fname);
t1 = (TTree*) gDirectory->FindObjectAny("t1");
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(const char* fname, Int_t SliceID)
{
TH1D *pol_x, *pol_z, *xisect, *yisect, *xfsect, *yfsect, *zhist;
TH2D *ini_hist, *bspot_hist;
TTree *t1;
TFile* f1 = new TFile(fname);
TString Cond,Title;
Double_t zSlice;
t1 = (TTree*) gDirectory->FindObjectAny("t1");
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",SliceID);
Cond.Form("save_detID== %04d",SliceID);
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(const char* fname, char* Var, Int_t Slice, Int_t SubPlot)
{
TH1D *histo;
TTree *t1;
TFile* f1 = new TFile(fname);
t1 = (TTree*) gDirectory->FindObjectAny("t1");
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(const char* fname, char* Varx, char* Vary, Int_t Slice, Int_t SubPlot)
{
TH1D *histo;
TTree *t1;
TFile* f1 = new TFile(fname);
TString Cond,Title,Histo;
Double_t zSlice;
t1 = (TTree*) gDirectory->FindObjectAny("t1");
gStyle->SetPalette(1,0);
TCanvas* c1 = new TCanvas("c1", Varx);
// 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(const char* fname, Int_t iSlice, Int_t fSlice)
{
TH2D *bspot_hist;
TTree *t1;
TFile* f1 = new TFile(fname);
TString Cond,Title;
Double_t zSlice,Xrms,Yrms,Xavg,Yavg,Nbins;
Int_t i;
t1 = (TTree*) gDirectory->FindObjectAny("t1");
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(100,-100,100,100,-100,100,100,-2600,1500)",Cond);
// t1->Draw("save_z:save_y:save_x>>bspot_hist(64,-40,40,64,-40,40,100,-290,290)",Cond);
bspot_hist = (TH2D*) gDirectory->Get("bspot_hist");
// 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%10.1f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\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");
}
void ToF(const char* fname,Int_t SaveID)
{
TTree *t1;
TFile* f1 = new TFile(fname);
TString Cond;
t1 = (TTree*) gDirectory->FindObjectAny("t1");
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();
}