void analysis(Int_t run_i,Int_t run_f,Int_t df_run) { gROOT->Reset(); gStyle->SetOptStat(0); TH1F *hsthe = new TH1F("hstheta","HMS theta in degree",100,0,50); TH1F *hsdel = new TH1F("hsdelta","HMS \delta E",100,-12,12); TH1F *hsen = new TH1F("hse","HMS energy in GeV",100,1.,5.); TH1F *Q2 = new TH1F("Q2","Q^{2} in GeV^{2}",100,0.,3.); TH1F *W = new TH1F("W","W in GeV",100,0.8,3.); Float_t Wmin = 1.05; //To apply dilution factor in [0.0,3.0] with 200 bins, Wmin should be divided by 0.015 GeV. Int_t Wbin = 42; // For 15MeV resoultion deltaW=1.05 for Wbin=70 Float_t Wbinsize = 0.015; // 15MeV Float_t Wmax = Wmin + Wbin*Wbinsize; TH2F *Q2vsW = new TH2F("Q2vsW","Q^{2} vs W",Wbin,Wmin,Wmax,100,0.,3.); //TH1F *W_p = new TH1F("W_p","W of positive helicity",Wbin,Wmin,Wmax); //TH1F *W_n = new TH1F("W_n","W of negative helicity",Wbin,Wmin,Wmax); //TH1D *W_C_p = new TH1D("W_C_p","W of positive helicity, charge normalized, LT correction muliplied",Wbin,Wmin,Wmax); //TH1D *W_C_n = new TH1D("W_C_n","W of negative helicity, charge normalized, LT correction muliplied",Wbin,Wmin,Wmax); //TH1D *W_CPbPt_p = new TH1D("W_C_p","W of positive helicity, charge normalized, divided by Pb*Pt",Wbin,Wmin,Wmax); //TH1D *W_CPbPt_n = new TH1D("W_C_n","W of negative helicity, charge normalized, divided by Pb*Pt",Wbin,Wmin,Wmax); TH1D *W_asym_raw_all = new TH1D("W_asym_raw_all","(N_{+} - N_{-})/(N_{+} + N_{-}) along W for all runs",Wbin,Wmin,Wmax); TH1D *Error_W_asym_raw_all = new TH1D("Error_W_asym_all","Accumulated inverse of squared error of W_asym",Wbin,Wmin,Wmax); TH1D *W_asym_CPbPt_all = new TH1D("W_asym_CPbPt_all","#frac{1}{PbPt} (N_{+}/C_{+} - N_{-}/C_{-})/(N_{+}/C_{+} + N_{-}/C_{-}) along W for all runs",Wbin,Wmin,Wmax); TH1D *Error_W_asym_CPbPt_all = new TH1D("Error_W_asym_CPbPt_all","Accumulated inverse of squared error of W_asym_CPbPt",Wbin,Wmin,Wmax); TH1D *W_asym_CPbPtdf_all = new TH1D("W_asym_CPbPtdf_all","#frac{1}{PbPtdf} (N_{+}/C_{+} - N_{-}/C_{-})/(N_{+}/C_{+} + N_{-}/C_{-}) along W for all runs",Wbin,Wmin,Wmax); TChain chain("h9010"); Long_t *id,*size,*flags,*mt; for (Int_t nrun=run_i; nrun<=run_f; nrun++) { Int_t nfile = 0; chain.Reset(""); for (Int_t index=1; index<10; index++) { if (gSystem->GetPathInfo(filename(nrun,index),id,size,flags,mt)==0) { chain.Add(filename(nrun,index)); cout << filename(nrun,index) << " loaded : " << endl; nfile++; } } if (nfile!=0) { //chain.Add(filename(72957,2)); //chain.Add(filename(72957,3)); TH1F *W_p = new TH1F("W_p","W of positive helicity",Wbin,Wmin,Wmax); TH1F *W_n = new TH1F("W_n","W of negative helicity",Wbin,Wmin,Wmax); //Kinematic variables Run_info *r = new Run_info(nrun); r->summary(); Int_t hwplate = r->hwplate(); Float_t Pb = r->Pb()/100.; // Pb in [0,1] Float_t Pt = r->Pt()/100.; // Pt in [0,1] Int_t hel_flip = 0; if(hwplate==1) hel_flip = -1; if(hwplate==0) hel_flip = 1; // cout << "hel_flip = " << hel_flip << " hwplate = " << hwplate << endl; Float_t C_p_over_n = 1.; if(hel_flip > 0) C_p_over_n = r->Charge_help()/r->Charge_helm(); if(hel_flip < 0) C_p_over_n = r->Charge_helm()/r->Charge_help(); // cout << "Charge_help = " << r->Charge_help() << " Charge_helm = " << r->Charge_helm() << " Cp/Cn = " << C_p_over_n << endl ; Float_t LT_cor = r->LT_cor(); cout << "LT_correction = " << LT_cor << endl; Kin_var *k = new Kin_var(r->BeamE(),r->P_HMS(),r->Theta_HMS()); //k->set(4700,2200,16); cout << k->Q2()/1e6 << " " << sqrt(k->hsca2pe()) << endl; k->summary(); //delete r; //delete k; TObjArray Hlist(0); Int_t nentries = chain.GetEntries(); Float_t hsdelta, hstheta, hse, hsshtrk; //Float_t polariza; Int_t helicite; const Int_t r2d = 180./TMath::Pi(); chain.SetBranchAddress("hsdelta",&hsdelta); chain.SetBranchAddress("hstheta",&hstheta); chain.SetBranchAddress("hse",&hse); chain.SetBranchAddress("hsshtrk",&hsshtrk); chain.SetBranchAddress("helicite",&helicite); //chain.SetBranchAddress("polariza",&polariza); cout << "Number of entries :" << nentries << endl; for (Int_t i = 0; i0. && (hsshtrk/hse)>0.4) { k->set(r->BeamE(),hse*1e3,hstheta*r2d); //k->set(r->BeamE(),3200.*(1.+hsdelta/100.),hstheta*r2d); Float_t sW = sqrt(k->W2())/1e3; // W in GeV //Float_t Pt = polariza/100.; // Pt in [0,1] hsthe->Fill(k->theta()); hsdel->Fill(hsdelta); hsen->Fill(hse); Q2->Fill(k->Q2()/1e6); W->Fill(sW); Q2vsW->Fill(sW,k->Q2()/1e6); if(TMath::Abs(Pt)>0.2) { if(helicite*hel_flip>0) { W_p->Fill(sW); //W_C_p->Fill(sW,LT_cor); // charge, LT correction //W_CPbPt_p->Fill(sW,1./Pb/Pt); // charge, Pb, Pt } else if(helicite*hel_flip<0) { W_n->Fill(sW); //W_C_n->Fill(sW,LT_cor*C_p_over_n); // charge, LT correction //W_CPbPt_n->Fill(sW,C_p_over_n/Pb/Pt); // charge, Pb, Pt } } } } //Getting asymmetry TH1D *W_p_p_n = new TH1D("W_p_p_n","N_{+} + N_{-} along W",Wbin,Wmin,Wmax); TH1D *W_p_m_n = new TH1D("W_p_m_n","N_{+} - N_{-} along W",Wbin,Wmin,Wmax); TH1D *W_CPbPt_p_p_n = new TH1D("W_CPbPt_p_p_n","#frac{1}{PbPt} (N_{+}/C_{+} + N_{-}/C_{-}) along W",Wbin,Wmin,Wmax); TH1D *W_CPbPt_p_m_n = new TH1D("W_CPbPt_p_m_n","#frac{1}{PbPt} (N_{+}/C_{+} - N_{-}/C_{-}) along W",Wbin,Wmin,Wmax); TH1D *W_asym_raw = new TH1D("W_asym_raw","(N_{+} - N_{-})/(N_{+} + N_{-}) along W",Wbin,Wmin,Wmax); TH1D *W_asym_CPbPt = new TH1D("W_asym_CPbPt","#frac{1}{PbPt} (N_{+}/C_{+} - N_{-}/C_{-})/(N_{+}/C_{+} + N_{-}/C_{-}) along W",Wbin,Wmin,Wmax); W_p_p_n->Add(W_p,W_n,1,1); W_p_m_n->Add(W_p,W_n,1,-1); W_asym_raw->Divide(W_p_m_n,W_p_p_n,1,1,""); W_CPbPt_p_p_n->Add(W_p,W_n,1.,1.*C_p_over_n); W_CPbPt_p_m_n->Add(W_p,W_n,1.,-1.*C_p_over_n); W_CPbPt_p_m_n->Scale(LT_cor/Pb/Pt,""); W_asym_CPbPt->Divide(W_CPbPt_p_m_n,W_CPbPt_p_p_n,1,1,""); //W_CPbPt_p_m_n->Add(W_CPbPt_p,W_CPbPt_n,1,-1); //W_asym_CPbPt->Divide(W_CPbPt_p_m_n,W_p_p_n,1,1,""); //W_p->Add(W_n,-1); //W_asym_raw->Divide(W_p,W_p_p_n,1,1,""); //Error calculation TH1D *E2_W_asym_raw = new TH1D("E2_W_asym_raw","Inverse of squared errors of W_asym_raw",Wbin,Wmin,Wmax); TH1D *W_asym_raw_weighted = new TH1D("W_asym_raw_weighted","W_asym_raw divided by error^2",Wbin,Wmin,Wmax); TH1D *E2_W_asym_CPbPt_CPbPt = new TH1D("E2_W_asym_CPbPt","Inverse of squared errors of W_asym_CPbPt",Wbin,Wmin,Wmax); TH1D *W_asym_CPbPt_weighted = new TH1D("W_asym_CPbPt_weighted","W_asym_CPbPt divided by error^2",Wbin,Wmin,Wmax); for(Int_t ibin = 1; ibin<=Wbin; ibin++) { cout << ibin << " " << W_p->GetBinCenter(ibin) << " " << W_p->GetBinContent(ibin) << endl; Float_t N_p = W_p->GetBinContent(ibin); Float_t N_n = W_n->GetBinContent(ibin); Float_t sigma2_W_asym_raw = 4*N_p*N_n/pow(N_p+N_n,3); Float_t sigma2_W_asym_CPbPt = 4*pow(LT_cor/Pb/Pt,2)*pow(C_p_over_n,2)*N_p*N_n*(N_p+N_n)/pow(N_p+N_n*C_p_over_n,4); //cout << ibin << " " << W_asym_raw->GetBinContent(ibin) << " " << W_asym_raw->GetBinError(ibin) << endl; W_asym_raw->SetBinError(ibin,sqrt(sigma2_W_asym_raw)); W_asym_CPbPt->SetBinError(ibin,sqrt(sigma2_W_asym_CPbPt)); //cout << ibin << " " << W_asym_raw->GetBinContent(ibin) << " " << W_asym_raw->GetBinError(ibin) << endl; E2_W_asym_raw->SetBinContent(ibin,1./sigma2_W_asym_raw); E2_W_asym_CPbPt->SetBinContent(ibin,1./sigma2_W_asym_CPbPt); } //Accumulation of asymmetry/sigma^2 and 1/sigma^2 W_asym_raw_weighted->Multiply(W_asym_raw,E2_W_asym_raw,1,1,""); W_asym_raw_all->Add(W_asym_raw_weighted,1); Error_W_asym_raw_all->Add(E2_W_asym_raw,1); W_asym_CPbPt_weighted->Multiply(W_asym_CPbPt,E2_W_asym_CPbPt,1,1,""); W_asym_CPbPt_all->Add(W_asym_CPbPt_weighted,1); Error_W_asym_CPbPt_all->Add(E2_W_asym_CPbPt,1); delete k; delete r; } } //Calculation of final asymmetry of all runs for(Int_t ibin = 1; ibin<=Wbin; ibin++) { Float_t sum_asym_raw_of_sigma2 = W_asym_raw_all->GetBinContent(ibin); Float_t sum_raw_of_sigma2 = Error_W_asym_raw_all->GetBinContent(ibin); Float_t sum_asym_CPbPt_of_sigma2 = W_asym_CPbPt_all->GetBinContent(ibin); Float_t sum_CPbPt_of_sigma2 = Error_W_asym_CPbPt_all->GetBinContent(ibin); W_asym_raw_all->SetBinContent(ibin,sum_asym_raw_of_sigma2/sum_raw_of_sigma2); W_asym_raw_all->SetBinError(ibin,sqrt(1./sum_raw_of_sigma2)); W_asym_CPbPt_all->SetBinContent(ibin,sum_asym_CPbPt_of_sigma2/sum_CPbPt_of_sigma2); W_asym_CPbPt_all->SetBinError(ibin,sqrt(1./sum_CPbPt_of_sigma2)); } //Loading the dilution factor (If df_run = 1, df = 1.0) Long_t *id,*size,*flags,*mt; if (gSystem->GetPathInfo(df_filename(df_run),id,size,flags,mt)==0) { //TFile *df_file = new TFile(df_filename(df_run)); //TH1D *df = (TH1D*)df_file->Get("h300"); TH1D *df = new TH1D("df","Dilution factor without error",Wbin,Wmin,Wmax); ifstream df_text; df_text.open(df_filename(df_run)); Int_t iWmin = ceil(Wmin/Wbinsize); cout << "iWmin = " << iWmin << endl; cout << "index \t W \t df \t error" << endl; Int_t line_number = 1; Int_t ibin = 1; TString line; while(!df_text.eof()) { line.ReadLine(df_text); if(line_number>iWmin && ibin<=Wbin) { TObjArray *sublist = line.Tokenize("\t "); TObjString *sub = sublist->At(3); Float_t wprnh = atof(sub->GetString().Data()); df->SetBinContent(ibin,wprnh); df->SetBinError(ibin,0); cout << line.Data() << endl; cout << ibin << " " << df->GetBinCenter(ibin) << " " << wprnh << " " << df->GetBinError(ibin) << endl; ibin++; } line_number++; } df_text.close(); TCanvas *df_c = new TCanvas("df_c","DF canvas"); df_c->cd(1); df->Draw("E"); W_asym_CPbPtdf_all->Divide(W_asym_CPbPt_all,df,1,1,""); } else { cout << "No dilution factor exists!" << endl; W_asym_CPbPtdf_all->Add(W_asym_CPbPt_all,1); } Hlist.Add(hsthe); Hlist.Add(hsdel); Hlist.Add(hsen); Hlist.Add(Q2); Hlist.Add(W); Hlist.Add(Q2vsW); Hlist.Add(W_asym_raw_all); Hlist.Add(W_asym_CPbPt_all); Hlist.Add(W_asym_CPbPtdf_all); Hlist.Add(df); //Drawing for check TCanvas *c1 = new TCanvas("c1","canvas1"); c1->Divide(2,2); c1->cd(1); hsthe->Draw(); c1->cd(2); hsen->Draw(); c1->cd(3); W->Draw(); c1->cd(4); Q2->Draw(); TCanvas *c2 = new TCanvas("c2","canvas2"); c2->cd(1); Q2vsW->SetOption("colz"); Q2vsW->Draw(); TCanvas *c3 = new TCanvas("c3","canvas3"); c3->Divide(2,2); c3->cd(1); W_p->Draw("E"); c3->cd(2); W_n->Draw("E"); c3->cd(3); W_p_p_n->Draw("E"); c3->cd(4); W_p_m_n->Draw("E"); TCanvas *c4 = new TCanvas("c4","canvas4"); c4->cd(1); W_asym_raw->Draw("E"); TCanvas *c5 = new TCanvas("c5","canvas5"); c5->cd(1); W_asym_CPbPt->Draw("E"); TCanvas *c6 = new TCanvas("c6","canvas6"); W_asym_raw_all->SetMaximum(0.1); W_asym_raw_all->SetMinimum(-0.1); c6->cd(1); W_asym_raw_all->Draw("E"); TCanvas *c7 = new TCanvas("c7","canvas7"); W_asym_CPbPt_all->SetMaximum(0.2); W_asym_CPbPt_all->SetMinimum(-0.2); c7->cd(1); W_asym_CPbPt_all->Draw("E"); TCanvas *c8 = new TCanvas("c8","canvas8"); W_asym_CPbPtdf_all->SetMaximum(1.0); W_asym_CPbPtdf_all->SetMinimum(-1.0); c8->cd(1); W_asym_CPbPtdf_all->Draw("E"); //Writing to file TFile *result = new TFile("new.root","RECREATE"); Hlist.Write(); result->Close(); } TString filename(Int_t nrun, Int_t index) { TString dir, prefix, suffix, s_number; dir = "ntup/"; prefix = "hms"; suffix = ".root"; s_number = Form("%d.%d",nrun,index); return dir + prefix + s_number + suffix; } TString df_filename(Int_t nrun) { TString dir, prefix, suffix, s_number; dir = "df/"; prefix = "df_"; suffix = ".out"; //".root"; s_number = Form("%d",nrun); return dir + prefix + s_number + suffix; } // Kinematics from report files class Run_info { private: Int_t Run; Float_t BeamE; Float_t HMS_Deadtime; Float_t Current; Float_t Runtime; Float_t P_HMS; Float_t Theta_HMS; Float_t Charge_help; Float_t Charge_helm; Int_t hwplate; Float_t Pb; Float_t Pt; Float_t LT_cor; Float_t line_elem(TString*, Int_t); public: Run_info(Int_t); ~Run_info(); void read(Int_t); void summary(); Int_t Run() { return Run; } Float_t BeamE() { return BeamE; } Float_t HMS_Deadtime() { return HMS_Deadtime; } Float_t Current() { return Current; } Float_t Runtime() { return Runtime; } Float_t P_HMS() { return P_HMS; } Float_t Theta_HMS() { return Theta_HMS; } Int_t hwplate() { return hwplate; } Float_t Pb() { return Pb; } Float_t Pt() { return Pt; } Float_t LT_cor() { return LT_cor; } Float_t Charge_help() { return Charge_help; } Float_t Charge_helm() { return Charge_helm; } }; Run_info::Run_info(Int_t nrun) { read(nrun); } Run_info::~Run_info() { } Run_info::read(Int_t nrun) { // BeamE in MeV // HMS_Deadtime as ratio between 0 to 1 (BeamOn) // Current in microA (BCM1 BeamOn) // Runtime in sec (BCM1 BeamOn) // P_HMS in MeV/c // Theta_HMS in degrees // Charge_help/m in microC // hwplate 1/0 // Pb, Pt in percent TString genfile = Form("scalers/gen%d.txt", nrun); TString run = Form("%d", nrun); TString line; Run=nrun; // Run infomation from E_Pb file ifstream Beaminfo; Beaminfo.open("E_Pb.txt"); while(!Beaminfo.eof()) { line.ReadLine(Beaminfo); if(line.BeginsWith(run.Data())) { cout << "E_Pb : " << line.Data() << endl; BeamE = line_elem(line, 1); hwplate = (Int_t)line_elem(line, 5); Pb = line_elem(line,6); } } Beaminfo.close(); // Live time correction from livetime_list file ifstream LTinfo; LTinfo.open("livetime_list.txt"); while(!LTinfo.eof()) { line.ReadLine(LTinfo); if(line.BeginsWith(run.Data())) { cout << "LT_correction : " << line.Data() << endl; LT_cor = line_elem(line, 3); } } LTinfo.close(); // Target polarization using 2011/02/09 J. Maxwell charge averaged polarization with 50 nA cut or no cut ifstream Ptinfo; Ptinfo.open("charge_averaged_runs.nocut.txt"); while(!Ptinfo.eof()) { line.ReadLine(Ptinfo); if(line.BeginsWith(run.Data())) { cout << "Pt : " << line.Data() << endl; Pt = line_elem(line, 1); } } Ptinfo.close(); ifstream Report; Report.open(genfile.Data()); while(!Report.eof()) { line.ReadLine(Report); if(line.BeginsWith("HMS_Deadtime")) HMS_Deadtime = line_elem(line,2); if(line.BeginsWith("Current")) Current = line_elem(line,2); if(line.BeginsWith("Runtime")) Runtime = line_elem(line,2); if(line.BeginsWith("P_HMS")) P_HMS = line_elem(line,2)*1e3; if(line.BeginsWith("Theta_HMS")) Theta_HMS = line_elem(line,2); if(line.BeginsWith("Charge_help")) Charge_help = line_elem(line,2); if(line.BeginsWith("Charge_helm")) Charge_helm = line_elem(line,2); } Report.close(); } void Run_info::summary() { cout << "Summary of Run information" << endl; cout << "(Run#, BeamE, HMS_Deadtime, Current, Runtime, P_HMS, Theta_HMS, hwplate, Pb, Pt, LT_correction)" << endl; cout << Run << ", "; cout << BeamE << ", "; cout << HMS_Deadtime << ", "; cout << Current << ", "; cout << Runtime << ", "; cout << P_HMS << ", "; cout << Theta_HMS << ", "; cout << hwplate << ", "; cout << Pb << ", "; cout << Pt << ", "; cout << LT_cor << endl; } Float_t Run_info::line_elem(TString line, Int_t index) { TObjArray *sublist = line.Tokenize("\t "); TObjString *sub = sublist->At(index); return atof(sub->GetString().Data()); } // Kinematic variables class Kin_var { private: Float_t sEi; Float_t sEf; Float_t stheta; Float_t stheta_degree; Float_t sPhi; Float_t sQ2; Float_t M; public: Kin_var(Float_t, Float_t, Float_t); ~Kin_var(); void set(Float_t, Float_t, Float_t); void summary(); Float_t Ei() { return sEi; } Float_t Ef() { return sEf; } Float_t theta() { return stheta_degree; } Float_t nu() { return sEi-sEf; } Float_t Q2() { return sQ2; } Float_t W2(); Float_t x_bj(); Float_t gamma(); Float_t hsca1pa(); Float_t hsca1pe(); Float_t hsca2pe(); Float_t hsna1(); Float_t hsna2(); }; Kin_var::Kin_var(Float_t Ei,Float_t Ef, Float_t theta) { set(Ei,Ef,theta); } Kin_var::~Kin_var() { } void Kin_var::set(Float_t Ei,Float_t Ef, Float_t theta) { // Energy in MeV, theta in degree, but stheta in radian sEi = Ei; sEf = Ef; stheta = theta*TMath::Pi()/180; stheta_degree = theta; sPhi = 0.0; sQ2 = 4*sEi*sEf*TMath::Sin(stheta/2)*TMath::Sin(stheta/2); M = 955.74; } void Kin_var::summary() { cout << "Summary of Kinematic variables" << endl; cout << "Ei = " << sEi << endl; cout << "Ef = " << sEf << endl; cout << "theta = " << stheta << endl; cout << "nu = " << nu() << endl; cout << "Q2 = " << Q2() << endl; cout << "W2 = " << W2() << endl; cout << "x_bj = " << x_bj() << endl; cout << "gamma = " << gamma() << endl; } Float_t Kin_var::W2() { return M*M + 2*M*nu() - sQ2; } Float_t Kin_var::x_bj() { return sQ2/(2*M*nu()); } Float_t Kin_var::gamma() { return 2*M*x_bj()/sqrt(sQ2); } //# Col 55: HSCA1PA = E-EP*COS(TH) //# Col 56: HSCA1PE = EP*SIN(TH)/COS(PHI) //# Col 57: HSCA2PE = (E-EP*COS(TH) )/(EP*SIN(TH)COS(PHI) //# Col 58: HSNA1 = 1/(E+EP) //# Col 59: HSNA2 = SQRT(Q2)/(2E) Float_t Kin_var::hsca1pa() { return sEi - sEf*TMath::Cos(stheta); } Float_t Kin_var::hsca1pe() { return sEf*TMath::Sin(stheta)/TMath::Cos(sPhi); } Float_t Kin_var::hsca2pe() { return (sEi - sEf*TMath::Cos(stheta))/(sEf*TMath::Sin(stheta)*TMath::Cos(sPhi)); } Float_t Kin_var::hsna1() { return 1/(sEi + sEf); } Float_t Kin_var::hsna2() { return sqrt(sQ2)/(2*sEi); }