This is the legacy ITS Confluence server. All content has been migrated to the new server located at:
For questions, please email atlassiansupport@virginia.edu
#if ROOT_VERSION_CODE >= ROOT_VERSION(6,00,0)
#include <TSystem.h>
#include <TStyle.h>
#include <TFile.h>
#include <TTree.h>
#include <TCanvas.h>
#include <TH1D.h>
#include <TLorentzVector.h>
#include "TFile.h"
#include "TTree.h"
#include "TH2D.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TMath.h"
#include <iostream>
#endif
using namespace std;
#define PI 3.14159265359
namespace {
int nfiles = 1;
const char* inputs [] = {
"trk_aug29_8M.root" //input file
};
}
void asymmetry() {
gStyle->SetOptFit();
TCanvas *c0 = new TCanvas("c0","c0"); c0->SetGrid();
TCanvas *c1 = new TCanvas("c1","c1"); c1->SetGrid();
TCanvas *c3 = new TCanvas("c3","c3"); c3->SetGrid();
TCanvas *c4 = new TCanvas("c4","c4"); c4->SetGrid();
TCanvas *c5 = new TCanvas("c5","c5"); c5->SetGrid();
TCanvas *c6 = new TCanvas("c6","c6"); c6->SetGrid();
TCanvas *c7 = new TCanvas("c7","c7"); c7->SetGrid();
TCanvas *c8 = new TCanvas("c8","c8"); c8->SetGrid();
TCanvas *c9 = new TCanvas("c9","c9"); c9->SetGrid();
for (int i=0; i<nfiles; ++i) {
TH1D *hx = new TH1D("hx","hx",20,0,1); //Ndimu vs Xtarget
TH1D *hnum = new TH1D("hnum","hnum",10, 0,10);//Ndimu vs mass
TH1D *hphi_s = new TH1D("hphi_s","hphi_s",30, -1*PI,PI);//Ndimu vs phi_target_rest_frame
TH1D *asym = new TH1D("Asym","Asym",5,0,0.5);//Sivers asymetry An
//Ndimu vs phi_target_rest_frame for each bin of xtarget:
TH1D *x_target[10];
char name1[100];
for (int ii=0; ii<10; ii++)
{
float xlow = 0.1*ii;
sprintf(name1,"xTarget_%f",xlow+0.05);
x_target[ii] = new TH1D(name1,name1, 20, -1*PI, PI);
x_target[ii]->Sumw2();
}
//
// Loop a TTree in a file
TFile *f = TFile::Open(inputs[i],"read");
TTree *T = (TTree*) f->Get("Truth");
#define MAX_PARTICLES 100
int gnhodo[MAX_PARTICLES];
float gpx[MAX_PARTICLES], gpy[MAX_PARTICLES], gpz[MAX_PARTICLES];
int n_particles, krecstat;
T->SetBranchAddress("gnhodo", &gnhodo);
T->SetBranchAddress("gpx", &gpx);
T->SetBranchAddress("gpy", &gpy);
T->SetBranchAddress("gpz", &gpz);
T->SetBranchAddress("krecstat", &krecstat);
float mu_mass = .105658;
float pro_mass = 0.93827;
float beamE = 120.0;
float gE[2];
float phi, phi_s;
for(int ientry=0; ientry<T->GetEntries(); ++ientry) {
T->GetEntry(ientry);
//Initial 4-vector definition
TLorentzVector mu0, mu1, dimu, hadron_b, hadron_t, initial,spin;
gE[0] = TMath::Sqrt(gpx[0]*gpx[0]+gpy[0]*gpy[0]+gpz[0]*gpz[0]+mu_mass*mu_mass);
gE[1] = TMath::Sqrt(gpx[1]*gpx[1]+gpy[1]*gpy[1]+gpz[1]*gpz[1]+mu_mass*mu_mass);
mu0.SetPxPyPzE(gpx[0], gpy[0], gpz[0], gE[0]);
mu1.SetPxPyPzE(gpx[1], gpy[1], gpz[1], gE[1]);
hadron_b.SetPxPyPzE(0.0,0.0,TMath::Sqrt(beamE*beamE-pro_mass*pro_mass),beamE);
hadron_t.SetPxPyPzE(0.0,0.0,0.0,pro_mass);
if (ientry % 2 ==0)
{
spin.SetPxPyPzE(0.0,1.0,0.0,0.0);
}
else
{
spin.SetPxPyPzE(0.0,-1.0,0.0,0.0);
}
initial = hadron_b + hadron_t;
dimu = mu0+mu1;
float dimu_gmass = dimu.Mag();
//boost to centre of mass
TLorentzRotation cm_r;
TVector3 cm_v = initial.BoostVector();
cm_r.Boost(-cm_v);
TLorentzVector hadron_b_cm = cm_r * hadron_b;
TLorentzVector hadron_t_cm = cm_r * hadron_t;
TLorentzVector mu0_cm = cm_r * mu0;
TLorentzVector mu1_cm = cm_r * mu1;
TLorentzVector dimu_cm = cm_r * dimu;
TLorentzVector spin_cm = cm_r * spin;
float s = 2.0*beamE*pro_mass + 2.0*pro_mass*pro_mass;
float sqrt_s = TMath::Sqrt(s);
float dimu_E = dimu_cm.E();
float dimu_pl = dimu_cm.Pz();
float dimu_pl_max = sqrt_s/2.0*(1-(dimu_gmass*dimu_gmass)/s);
float xf = dimu_pl/dimu_pl_max;
float y = 0.5*TMath::Log((dimu_E+dimu_pl)/(dimu_E-dimu_pl));
float tau = (dimu_gmass*dimu_gmass)/s;
float xb = 0.5*(TMath::Sqrt(xf*xf+4*tau)+xf);
float xt = 0.5*(TMath::Sqrt(xf*xf+4*tau)-xf);
//boost to dilepton restframe
TLorentzRotation cm_dimu;
TVector3 cm_v2 = dimu_cm.BoostVector();
cm_dimu.Boost(-cm_v2);
TLorentzVector hadron_b_cm_dimu = cm_dimu * hadron_b_cm;
TLorentzVector hadron_t_cm_dimu = cm_dimu * hadron_t_cm;
TLorentzVector mu0_cm_dimu = cm_dimu * mu0_cm;
TLorentzVector mu1_cm_dimu = cm_dimu * mu1_cm;
TLorentzVector dimu_cm_dimu = cm_dimu * dimu_cm;
TLorentzVector spin_cm_dimu = cm_dimu * spin_cm;
TVector3 zaxis = (hadron_b_cm_dimu.Vect()).Unit()-(hadron_t_cm_dimu.Vect()).Unit();
TVector3 zunit = zaxis.Unit();
TVector3 hadron_plane_norm = zunit.Cross((hadron_b_cm_dimu.Vect()).Unit());
TVector3 lepton_plane_norm = zunit.Cross((mu0_cm_dimu.Vect()).Unit());
TVector3 yunit = lepton_plane_norm.Unit();
TVector3 xunit = (yunit.Cross(zunit)).Unit();
//below are phi_s_tf from target rest frame. The frame for the Asymmetry calculation
TVector3 spin_vect_tf = spin.Vect();
TVector3 z_prime_unit = (hadron_b.Vect()).Unit();
TVector3 y_prime_unit = (z_prime_unit.Cross(dimu.Vect())).Unit();
TVector3 x_prime_unit = (y_prime_unit.Cross(z_prime_unit)).Unit();
float phi_s_tf;
float cosphi_s_tf = (x_prime_unit.Dot(spin_vect_tf))/((x_prime_unit.Mag())*(spin_vect_tf.Mag()));
if ( acos( spin_vect_tf.Dot(y_prime_unit) / ( (spin_vect_tf).Mag() * y_prime_unit.Mag() ) ) > (PI/2) )
{
phi_s_tf = -acos(cosphi_s_tf);
}
else
{
phi_s_tf = acos(cosphi_s_tf);
}
int bin_xtarget = xt/0.1;
cerr<<ientry<<" "<<T->GetEntries()<<" "<<dimu_gmass<<" "<<xt<<" "<<phi_s_tf<<" "<<endl;
//good_event status -> currently only for successfully reconstructed event. We can also add trigger condition later for the good event status:
bool good_event=false;
if (krecstat==0)
{
good_event=true;
}
if(good_event) //only for successfully reconstructed event
{
hnum->Fill(dimu_gmass);
hx->Fill(xt);
hphi_s->Fill(phi_s_tf);
if (bin_xtarget >=0 && bin_xtarget < 10)
{
x_target[bin_xtarget]->Fill(phi_s_tf);
}
}
}
// End of loop
cerr<<"end of loop"<<endl;
cerr<<endl;
cerr<<"get An"<<endl;
//Asymmetry calculation:
for (int j=0; j<5; j++)
{
float sum=0;
float sum_N=0;
float An=0;
float num_err=0;
float den_err=0;
float tot_err=0;
float new_err=0;
for (int k=1; k<=20; k++)
{
float phi_angle = -1.0*PI+(k-1)*2.0*PI/20.0+0.5*2.0*PI/20.0;
float a = x_target[j]->GetBinContent(k);
float b = x_target[j]->GetBinError(k);
sum = sum + a*sin(phi_angle);
sum_N = sum_N+a;
num_err = num_err + (b*sin(phi_angle))*(b*sin(phi_angle));
den_err = den_err + b*b;
//New error calculation from Kenichi. Kenichi pointed out that they are correlated:
new_err = new_err + pow(b * (sin(phi_angle) * sum_N - sum) / pow(sum_N, 2), 2);
}
num_err = TMath::Sqrt(num_err);
den_err = TMath::Sqrt(den_err);
An = 2*sum/sum_N;
tot_err = 2*An*TMath::Sqrt((num_err/sum)*(num_err/sum)+(den_err/sum_N)*(den_err/sum_N));
if(fabs(An)<=1)
{
asym->SetBinContent(j+1,An);
//asym->SetBinError(j+1,tot_err);
//applying new error calculation from kenichi:
new_err = 2*TMath::Sqrt(new_err);
asym->SetBinError(j+1,new_err);
}
cerr<<"Bin Xtarget = "<<j+1<<" An = "<<An<<endl;
}
//some cosmetics
hnum->SetMarkerStyle(20);
hx->SetMarkerStyle(20);
hphi_s->SetMarkerStyle(20);
asym->SetMarkerStyle(20);
for (int j=0; j<10; j++)
{
x_target[j]->SetMarkerStyle(20);
}
int color = kBlack;
if(i==0) {
color = kBlack;
c0->cd();
//c0->SetLogy();
hnum->SetTitle("nDimu-trig. vs. M_{#mu#mu}; M_{#mu#mu} [GeV/c^{2}]; nDimu-trig.");
hnum->SetMarkerColor(color);
hnum->SetLineColor(color);
hnum->SetMinimum(1);
hnum->Draw("e");
hnum->SetStats(0);
c1->cd();
//c1->SetLogy();
hx->SetTitle("nDimu-trig. vs. x_target; x_target; nDimu-trig.");
hx->SetMarkerColor(color);
hx->SetLineColor(color);
hx->SetMinimum(1);
hx->Draw("e");
hx->SetStats(0);
c3->cd();
hphi_s->SetTitle("nDimu-trig. vs. #phi_{s}; #phi_{s}; nDimu-trig.");
hphi_s->SetMarkerColor(color);
hphi_s->SetLineColor(color);
hphi_s->SetMinimum(0.);
//hphi->SetMaximum(1.1);
hphi_s->Draw("e");
hphi_s->SetStats(0);
c4->cd();
x_target[0]->SetTitle("nDimu-trig. vs. #phi_{s} (xt = 0.05) ; #phi_{s}; nDimu-trig.");
x_target[0]->SetMarkerColor(color);
x_target[0]->SetLineColor(color);
x_target[0]->SetMinimum(0.);
x_target[0]->Draw("e");
x_target[0]->SetStats(0);
c5->cd();
x_target[1]->SetTitle("nDimu-trig. vs. #phi_{s} (xt = 0.15) ; #phi_{s}; nDimu-trig.");
x_target[1]->SetMarkerColor(color);
x_target[1]->SetLineColor(color);
x_target[1]->SetMinimum(0.);
x_target[1]->Draw("e");
x_target[1]->SetStats(0);
c6->cd();
x_target[2]->SetTitle("nDimu-trig. vs. #phi_{s} (xt = 0.25) ; #phi_{s}; nDimu-trig.");
x_target[2]->SetMarkerColor(color);
x_target[2]->SetLineColor(color);
x_target[2]->SetMinimum(0.);
x_target[2]->Draw("e");
x_target[2]->SetStats(0);
c7->cd();
x_target[3]->SetTitle("nDimu-trig. vs. #phi_{s} (xt = 0.35) ; #phi_{s}; nDimu-trig.");
x_target[3]->SetMarkerColor(color);
x_target[3]->SetLineColor(color);
x_target[3]->SetMinimum(0.);
x_target[3]->Draw("e");
x_target[3]->SetStats(0);
c8->cd();
x_target[4]->SetTitle("nDimu-trig. vs. #phi_{s} (xt = 0.45) ; #phi_{s}; nDimu-trig.");
x_target[4]->SetMarkerColor(color);
x_target[4]->SetLineColor(color);
x_target[4]->SetMinimum(0.);
x_target[4]->Draw("e");
x_target[4]->SetStats(0);
c9->cd();
asym->SetTitle("Drell-Yan Target Single-Spin Asymmetry; X_{target}; A_{N}");
asym->SetMarkerColor(color);
asym->SetLineColor(color);
asym->SetMinimum(-0.5);
asym->SetMaximum(0.5);
asym->Draw("e");
asym->SetStats(0);
}
else {
color = i+1;
c0->cd();
hnum->SetMarkerColor(color);
hnum->SetLineColor(color);
hnum->Draw("esame");
c1->cd();
hx->SetMarkerColor(color);
hx->SetLineColor(color);
hx->Draw("esame");
c3->cd();
hphi_s->SetMarkerColor(color);
hphi_s->SetLineColor(color);
hphi_s->Draw("esame");
}
}
}