2006-02-25 19:42:00 +00:00

598 lines
13 KiB
C

gStyle->SetPalette(1);
gStyle->SetErrorX(0);
gStyle->SetOptStat(0);
TFile *trans = TFile::Open("SDInnerScint.root");
TTree *tree= (TTree*)trans->Get("tree");
TFile *foil = TFile::Open("SDCryoE0.root");
TTree *plane= (TTree*)foil->Get("tree");
TCanvas *ec1;
TCanvas *ec2;
double meanMU, meanmuplus, ie ;
double implE;
double Bfield=0;
TGraph *g1, *g2;
void esquisse(){
gStyle->SetPalette(1);
ec1= new TCanvas();
ec2= new TCanvas();
ec1->SetGrid();
ec1->Divide(2,2);
TString poscut="positionz==14";
ec1->cd(1);
plane->Draw("positionz>>muPOS(1500,14,14.00003)");
ec1->cd(2);
plane->Draw("muon",poscut);
ec1->cd(3);
plane->Draw("globaltime>>hTOF(150)",poscut,"E1");
ec1->cd(4);
plane->Draw("kenergy+edeposit>>hKE(150)",poscut,"E1");
ec2->SetGrid();
ec2->Divide(2,2);
TFile *foil = TFile::Open("SDCryoE0.root");
TTree *plane= (TTree*)foil->Get("tree");
ec2->cd(1);
plane->Draw("positiony:positionx>>mxy(100,100)",poscut+"&&muon==1","zcol");
ec2->cd(2);
plane->Draw("positiony:positionx>>Mxy(100,100)",poscut+"&&muonium==1","zcol");
ec2->cd(3);
plane->Draw("globaltime>>mpg(100)",poscut+"&&muon==1","E1");
ec2->cd(4);
plane->Draw("globaltime>>mg(100)",poscut+"&&muonium==1","E1");
}
void moyennes(int ID=-1)
{
gStyle->SetPalette(1);
gStyle->SetErrorX(0);
gStyle->SetOptStat(11);
TCanvas *ec1;
TCanvas *ec2;
double meanMU, meanmuplus, ie ;
double implE;
TString cut = "positionz==14";
if(ID!=-1)
{
cut=cut+"&&runID=="+FS(ID);
}
plane->Draw("globaltime>>hmuplus(1000,0,500)",cut+"&&muon==1");
meanmuplus = hmuplus->GetMean();
plane->Draw("globaltime>>hMU(1000,0,500)",cut+"&&muonium==1");
meanMU = hMU->GetMean();
// plane->Draw("Init_Energy>>hIE");
plane->Draw("kenergy+edeposit>>hIE(500,0,5)",cut+"&&muon==1");
ie=hIE->GetMean();
ofstream outf;
outf.open("means.dat");
if(outf.is_open())
{
std::cout<<meanmuplus << " " << meanMU<<" " << ie <<std::endl;
outf<<meanmuplus << " " << meanMU<<" " << ie <<std::endl;
outf.close();
std::cout<<"Output file means.dat created."<<std::endl;
}
else
{
cout << "Error opening file";
return 0;
}
}
void mtime(TString particle="mu+",TString grids="ON", TString guards="ON"){
gDirectory->DeleteAll();
int ID;
std::cout<<"Run? (-1 for all): " <<std::endl;
std::cin>>ID;
std::cout<<"Magnetic Field [Gauss]?" <<std::endl;
std::cin>>Bfield;
// std::cout<<"Implantation Energy?" <<std::endl;
// std::cin>>implE;
int b_f=1;
if (Bfield==0)b_f=0;
#define GMUON 0.01355
#define TAU_MUON 2.197
#define ASYM 0.28
Int_t i,j;
double r;
char func_str[256];
char text_str[256];
GetValues();
sprintf(func_str, "exp(-x/%f)", TAU_MUON);
TF1 *f1 = new TF1("f1", func_str, 0., 20.);
//sprintf(func_str, "exp(-x/%f)", TAU_MUON);
TF1 *f2 = new TF1("f2", func_str, 0., 20.);
TF1 *fitdecay = new TF1("fitdecay", TFdecayfit, meanmuplus, 10000., 2);
TF1 *fitasym1 = new TF1("fitasym1", TFasymuonfit, meanmuplus, 10000., 3);
TF1 *fitasym2 = new TF1("fitasym2", TFasymuonfit, meanmuplus, 10000., 3);
/*----------- initialise fit hforw and hback ----------------*/
fitdecay->SetParNames("N(0)","MuonLifeTime");
fitdecay->SetParameters(0, 450.);
fitdecay->SetParameters(1, 2.2);
// convert frequence from Hz 1/s to 1/ns (10^-9 Hz)
double fz = 13.*b_f;
double ampl = 0.2;
if(particle=="Mu")
{
fz = -1.39e3*b_f;
ampl = 0.16;
}
fitasym1->SetParNames("Amplitude","Frequency [kHz]","Offset/Pi");
fitasym1->SetParameters(ampl,fz,0*b_f);
fitasym1->SetLineColor(46);
fitasym2->SetParNames("Amplitude","Frequency [kHz]","Offset/Pi");
fitasym2->SetParameters(ampl,fz,0.5*b_f);
fitasym2->SetLineColor(50);
//===================================================================
//===================================================================
/*---------------------------------------------------------*/
TH1F *hRight = new TH1F("hRight", "hRight",100, 0., 2400.);
TH1F *hLeft = new TH1F("hLeft", "hLeft",100, 0., 2400.);
TH1F *LRdiff = new TH1F("LRdiff", "LRdiff",100, 0., 2400.);
TH1F *LRsum = new TH1F("LRsum", "LRsum",100, 0., 2400.);
TH1F *LRasym = new TH1F("LRasym", "LRasym",100, 0., 2400.);
/*---------------------------------------------------------*/
TH1F *hTop = new TH1F("hTop", "hTop",100, 0., 2400.);
TH1F *hBottom = new TH1F("hBottom", "hBottom",100, 0., 2400.);
TH1F *BTdiff = new TH1F("BTdiff", "BTdiff",100, 0., 2400.);
TH1F *BTsum = new TH1F("BTsum", "BTsum",100, 0., 2400.);
TH1F *BTasym = new TH1F("BTasym", "BTasym",100, 0., 2400.);
/*---------------------------------------------------------*/
if(particle=="Mu")
{
hRight->SetBins( 100, 0., meanMU+200.);
hLeft->SetBins( 100, 0., meanMU+200.);
LRdiff->SetBins( 100, 0., meanMU+200.);
LRsum->SetBins( 100, 0., meanMU+200.);
LRasym->SetBins( 100, 0., meanMU+200.);
/*---------------------------------------------------------*/
hTop->SetBins(100, 0., 500.);
hBottom->SetBins(100, 0., 500.);
BTdiff->SetBins(100, 0., 500.);
BTsum->SetBins(100, 0., 500.);
BTasym->SetBins(100, 0., 500.);
}
/*---------------------------------------------------------*/
hLeft->Sumw2();
hRight->Sumw2();
hTop->Sumw2();
hBottom->Sumw2();
LRdiff->Sumw2();
LRsum->Sumw2();
LRasym->Sumw2();
BTdiff->Sumw2();
BTsum->Sumw2();
BTasym->Sumw2();
//===================================================================
//===================================================================
// moyennes();
TString cut2;
cut2 = "";
if(ID==-1)
{
cut2 = "";
}
else
{
cut2 = "positron.runID=="+FS(ID)+"&&";
}
tree->Draw("positron.globaltime>>hLeft",cut2+"positron.Left==1","");
tree->Draw("positron.globaltime>>hRight",cut2+"positron.Right==1","");
tree->Draw("positron.globaltime>>hBottom",cut2+"positron.Bottom==1","");
tree->Draw("positron.globaltime>>hTop",cut2+"positron.Top==1","");
LRsum->Add(hLeft,hRight,1,1);
LRdiff->Add(hLeft,hRight,1,-1);
LRasym->Divide(LRdiff,LRsum);
hLeft->SetLineColor(1);
hTop->SetLineColor(2);
hRight->SetLineColor(3);
hBottom->SetLineColor(4);
//=============================================================================
// LEFT-RIGHT
//=============================================================================
c1=new TCanvas("c1","LR Counts");
c1->SetGrid();
c1->SetFillColor(0);
c1->SetFrameFillColor(0);
c1->SetFrameBorderMode(0);
hLeft->SetTitle("Counts in Left and Right Scintillators");
hLeft->GetYaxis()->SetTitle("Counts (arbitrary scale)");
hLeft->GetXaxis()->SetTitle("Time [ns]");
hLeft->SetMarkerStyle(3);
hLeft->SetMarkerColor(4);
hLeft->SetLineColor(4);
hLeft->Draw("E1lp");
pt = new TPaveText(0.01,0.91,0.6,0.99,"blNDC");
pt->SetName("title");
pt->SetTextSize(0.032);
pt->SetBorderSize(0);
pt->SetFillColor(0);
text = pt->AddText("");
pt->Draw();
hRight->SetMarkerStyle(26);
hRight->SetMarkerColor(2);
hRight->SetLineColor(2);
hRight->Draw("E1lpsame");
TOF_plot();
leg = new TLegend(0.2,0.75,0.43,0.88);
leg->SetFillColor(0);
leg->SetBorderSize(1);
leg->AddEntry(hLeft,"Left ","p");
leg->AddEntry(hRight,"Right ","p");
leg->AddEntry(g1,"#mu+ TOF ("+FS(meanmuplus)+" ns)","l");
leg->AddEntry(g2,"Mu TOF("+FS(meanMU)+" ns)","l");
leg->Draw();
//++++++++++++++++++++++++++++++++++
//++++++++++++++++++++++++++++++++++
c3=new TCanvas("c3","LR Asymmetry");
c3->SetGrid();
c3->SetFillColor(0);
c3->SetFrameFillColor(0);
c3->SetFrameBorderMode(0);
TString title = "Left/Right Asymmetry for "+FS(implE)+" KeV impl. energy. (B="+FS(Bfield)+"G)";
LRasym->SetTitle(title);
LRasym->GetXaxis()->SetTitle("Time [ns]");
LRasym->GetYaxis()->SetTitle("Asymmetry");
LRasym->GetYaxis()->SetRangeUser(-1,1);
LRasym->SetMarkerStyle(2);
LRasym->SetMarkerColor(38);
LRasym->SetLineColor(38);
LRasym->Fit("fitasym1");
// gStyle->SetOptStat(0);
gStyle->SetOptFit(111);
LRasym->Draw("e1");
pt = new TPaveText(0.01,0.91,0.6,0.99,"blNDC");
pt->SetName("title");
pt->SetBorderSize(0);
pt->SetTextSize(0.032);
pt->SetFillColor(0);
text = pt->AddText("");
pt->Draw();
TOF_plot();
leg = new TLegend(0.2,0.75,0.43,0.88);
leg->SetFillColor(0);
leg->SetBorderSize(1);
leg->AddEntry(LRasym,"L-R Asym","p");
leg->AddEntry(g1,"#mu+ TOF ("+FS(meanmuplus)+" ns)","l");
leg->AddEntry(g2,"Mu TOF("+FS(meanMU)+" ns)","l");
leg->Draw();
//=============================================================================
// Bottom-Top
//=============================================================================
gStyle->SetPalette(1);
gStyle->SetOptStat(11);
gStyle->SetErrorX(0);
BTsum->Add(hBottom,hTop,1,1);
BTdiff->Add(hBottom,hTop,1,-1);
BTasym->Divide(BTdiff,BTsum);
BTasym->SetMarkerStyle(3);
c4=new TCanvas("c4","BT Counts");
c4->SetGrid();
c4->SetFillColor(0);
c4->SetFillColor(0);
c4->SetFrameFillColor(0);
c4->SetFrameBorderMode(0);
hTop->SetTitle("Counts in Top and Bottom Scintillators");
hTop->SetMarkerStyle(26);
hTop->SetMarkerColor(2);
hTop->SetLineColor(2);
hTop->GetYaxis()->SetTitle("Counts (arbitrary scale)");
hTop->GetXaxis()->SetTitle("Time [ns]");
hTop->Draw("lp");
pt = new TPaveText(0.01,0.91,0.6,0.99,"blNDC");
pt->SetName("title");
pt->SetBorderSize(0);
pt->SetTextSize(0.032);
pt->SetFillColor(0);
text = pt->AddText("");
pt->Draw();
hBottom->SetMarkerStyle(3);
hBottom->SetMarkerColor(4);
hBottom->SetLineColor(4);
hBottom->Draw("lpsame");
TOF_plot();
leg = new TLegend(0.2,0.75,0.43,0.88);
leg->SetFillColor(0);
leg->SetBorderSize(1);
leg->AddEntry(hBottom,"Bottom ","p");
leg->AddEntry(hTop,"Top ","p");
leg->AddEntry(g1,"#mu+ TOF("+FS(meanmuplus)+" ns)","l");
leg->AddEntry(g2,"Mu TOF("+FS(meanMU)+" ns)","l");
leg->Draw();
//++++++++++++++++++++++++++++++++++
//++++++++++++++++++++++++++++++++++
c6=new TCanvas("c6","BT Asymmetry");
c6->SetGrid();
c6->SetFillColor(0);
c6->SetFrameFillColor(0);
c6->SetFrameBorderMode(0);
title = "Bottom/Top Asymmetry for "+FS(implE)+" KeV impl. energy. (B="+FS(Bfield)+"G)";
BTasym->SetTitle(title);
BTasym->GetXaxis()->SetTitle("Time [ns]");
BTasym->GetYaxis()->SetTitle("Asymmetry");
BTasym->GetYaxis()->SetRangeUser(-1,1);
BTasym->SetMarkerColor(31);
BTasym->SetLineColor(31);
// BTasym->GetXaxis()->SetRangeUser(0,5000);
BTasym->Fit("fitasym2");
gStyle->SetOptStat(11);
gStyle->SetOptFit(111);
BTasym->Draw("p");
pt = new TPaveText(0.01,0.91,0.6,0.99,"blNDC");
pt->SetName("title");
pt->SetTextSize(0.032);
pt->SetBorderSize(0);
pt->SetFillColor(0);
text = pt->AddText("");
pt->Draw();
TOF_plot();
leg = new TLegend(0.2,0.75,0.43,0.88);
leg->SetFillColor(0);
leg->SetBorderSize(1);
leg->AddEntry(BTasym,"B-T Asym","p");
leg->AddEntry(g1,"#mu+ TOF("+FS(meanmuplus)+" ns)","l");
leg->AddEntry(g2,"Mu TOF("+FS(meanMU)+" ns)","l");
leg->Draw();
// c6->SetEditable(0);
}
TString flt2string(const float& number)
{
ostringstream oss;
oss << number;
return oss.str();
}
TString FS(const float& number)
{
ostringstream oss;
oss << number;
return oss.str();
}
saveall( TString save_ans="")
{
c1->SaveAs("LRcounts"+save_ans+".eps");
c3->SaveAs("LRasym"+save_ans+".eps");
c4->SaveAs("BTcounts"+save_ans+".eps");
c6->SaveAs("BTasym"+save_ans+".eps");
/* c1->SaveAs("LRcounts"+FS(implE)+save_ans+".gif");
c3->SaveAs("LRasym"+FS(implE)+save_ans+".gif");
c4->SaveAs("BTcounts"+FS(implE)+save_ans+".gif");
c6->SaveAs("BTasym"+FS(implE)+save_ans+".gif");*/
}
void TOF_plot()
{
Double_t a[2],b[2];
a[0]=meanmuplus;a[1]=meanmuplus;
b[0]=-1;b[1]=2500;
g1= new TGraph(2,a,b);
a[0]=meanMU;a[1]=meanMU;
b[0]=-1;b[1]=2500;
g2= new TGraph(2,a,b);
g1->SetLineColor(1);
g1->SetLineWidth(2);
g1->Draw("");
g2->SetLineColor(2);
g2->SetLineWidth(2);
g2->Draw("");
}
void GetValues()
{
ifstream inf;
inf.open("means.dat");
while(!inf.is_open())
{
std::cout<<"No 'means.dat' file.\nAttempt to create it."<<std::endl;
moyennes();
inf.open("means.dat");
}
inf >> meanmuplus >> meanMU >> implE;// ie;
inf.close();
}
Double_t MuonLifeFunc(Double_t *x, Double_t *par) {
// two fit parameters: the normalisation
// the muon life time
return par[0]*TMath::Exp(-x[0]/par[1]);
}
Double_t TFasymuonfit(Double_t *x, Double_t *par) {
// one fit parameter: the asymmetry
Double_t fitval, arg;
arg=-x[0]*Bfield*6.28*par[1]*1.e-6+3.14*par[2];
fitval = par[0]*TMath::Cos(arg);
return fitval;
}
Double_t TFdecayfit(Double_t *x, Double_t *par) {
//
Double_t fitval;
fitval = MuonLifeFunc(x, par);
return fitval;
}
format()
{
for(Int_t i==0; i<3; i++)
{
CASYM->cd(i)->SetFrameBorderMode(0);
CASYM->cd(i)->SetFrameBorderSize(0);
CASYM->cd(i)->SetGrid();
CASYM->cd(i)->SetEditable(1);
}
format(TCanvas cv)
{
cv->SetFrameBorderMode(0);
cv->SetFrameBorderSize(0);
cv->SetGrid();
cv->SetEditable(0);
}
}