[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / macros / PlotTherm_r3.C

#include <math.h>
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <Riostream.h>
#include <complex>

#include "TObject.h"
#include "TTree.h"
#include "TBranch.h"
#include "TLeaf.h"
#include "TVector2.h"
#include "TVector3.h"
#include "TFile.h"
#include "TString.h"
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TProfile.h"
#include "TProfile2D.h"
#include "TMath.h"
#include "TText.h"
#include "TRandom3.h"
#include "TArray.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TMinuit.h"
#include "TChain.h"
#include "TLorentzVector.h"
#include "TPaveStats.h"
#include "TLatex.h"
#include "TCanvas.h"
#include "TPad.h"

#define PI 3.1415926
#define BohrR 1963.6885 // Mate's value(1963.6885) ~ 387.5/0.197327(1963.7455)
#define FmToGeV 0.197327 // 0.197327

using namespace std;

void PlotTherm_r3(){

  double epsilon=1.0;// chaoticity

  gStyle->SetOptStat(0);
  //gStyle->SetOptFit(111);
 
  //TFile *infile=new TFile("Therm_r3_b9.root","READ");// For v2 paper (r<100 and r*<80)
  TFile *infile=new TFile("Therm_FSI_b2.root","READ");
  //TFile *infile=new TFile("Therm_dist_temp.root","READ");
 
  //
  TH3D *r3num_init=(TH3D*)infile->Get("r3num");
  TH3D *r3den1_init=(TH3D*)infile->Get("r3den1");
  TH3D *r3den2_init=(TH3D*)infile->Get("r3den2");
  TH3D *r3den3_init=(TH3D*)infile->Get("r3den3");
  TH3D *r3numSq_init=(TH3D*)infile->Get("r3numSq");
  TH3D *r3den1Sq_init=(TH3D*)infile->Get("r3den1Sq");
  TH3D *r3den2Sq_init=(TH3D*)infile->Get("r3den2Sq");
  TH3D *r3den3Sq_init=(TH3D*)infile->Get("r3den3Sq");
  TH3D *r3numEn_init=(TH3D*)infile->Get("r3numEn");
  TH3D *r3den1En_init=(TH3D*)infile->Get("r3den1En");
  TH3D *r3den2En_init=(TH3D*)infile->Get("r3den2En");
  TH3D *r3den3En_init=(TH3D*)infile->Get("r3den3En");
  TH3D *r3numME_init=(TH3D*)infile->Get("r3numME");
  TH3D *r3den1ME_init=(TH3D*)infile->Get("r3den1ME");
  TH3D *r3den2ME_init=(TH3D*)infile->Get("r3den2ME");
  TH3D *r3den3ME_init=(TH3D*)infile->Get("r3den3ME");
  
  TH2D *R2_2D = (TH2D*)infile->Get("Num_Cos_ss");
  //TH2D *R2Sq_2D = (TH2D*)infile->Get("NumSq_Cos_ss");
  TH2D *R2En_2D = (TH2D*)infile->Get("Den_ss");
  TH1D *R2=(TH1D*)R2_2D->ProjectionY();
  //TH1D *R2Sq=(TH1D*)R2Sq_2D->ProjectionY();
  TH1D *R2En=(TH1D*)R2En_2D->ProjectionY();
  R2->Divide(R2En);
  //R2Sq->Divide(R2En);

  TH3D *r3num=(TH3D*)r3num_init->Clone();
  TH3D *r3den1=(TH3D*)r3den1_init->Clone();
  TH3D *r3den2=(TH3D*)r3den2_init->Clone();
  TH3D *r3den3=(TH3D*)r3den3_init->Clone();
  TH3D *r3numSq=(TH3D*)r3numSq_init->Clone();
  TH3D *r3den1Sq=(TH3D*)r3den1Sq_init->Clone();
  TH3D *r3den2Sq=(TH3D*)r3den2Sq_init->Clone();
  TH3D *r3den3Sq=(TH3D*)r3den3Sq_init->Clone();
  TH3D *r3numEn=(TH3D*)r3numEn_init->Clone();
  TH3D *r3den1En=(TH3D*)r3den1En_init->Clone();
  TH3D *r3den2En=(TH3D*)r3den2En_init->Clone();
  TH3D *r3den3En=(TH3D*)r3den3En_init->Clone();
  TH3D *r3numME=(TH3D*)r3numME_init->Clone();
  TH3D *r3den1ME=(TH3D*)r3den1ME_init->Clone();
  TH3D *r3den2ME=(TH3D*)r3den2ME_init->Clone();
  TH3D *r3den3ME=(TH3D*)r3den3ME_init->Clone();


  int FB=40;
  int LB=50;
  /*double SF = r3numEn->Integral(FB,LB,FB,LB,FB,LB)/r3numME->Integral(FB,LB,FB,LB,FB,LB);
  cout<<SF<<endl;
  r3numME->Scale(SF);
  SF = r3den1En->Integral(FB,LB,FB,LB,FB,LB)/r3den1ME->Integral(FB,LB,FB,LB,FB,LB);
  cout<<SF<<endl;
  r3den1ME->Scale(SF);
  SF = r3den2En->Integral(FB,LB,FB,LB,FB,LB)/r3den2ME->Integral(FB,LB,FB,LB,FB,LB);
  cout<<SF<<endl;
  r3den2ME->Scale(SF);
  SF = r3den3En->Integral(FB,LB,FB,LB,FB,LB)/r3den3ME->Integral(FB,LB,FB,LB,FB,LB);
  cout<<SF<<endl;
  r3den3ME->Scale(SF);*/
  //
  /*
  //////////////////////
  // Symmetrize bin contents
  for(int i=1; i<=50; i++){
    for(int j=1; j<=50; j++){
      for(int k=1; k<=50; k++){
	// 3-pion phase
	double num_sum = r3num_init->GetBinContent(i,j,k) + r3num_init->GetBinContent(i,k,j) + r3num_init->GetBinContent(j,i,k);
	num_sum += r3num_init->GetBinContent(j,k,i) + r3num_init->GetBinContent(k,i,j) + r3num_init->GetBinContent(k,j,i);
	double Sq_sum = r3numSq_init->GetBinContent(i,j,k) + r3numSq_init->GetBinContent(i,k,j) + r3numSq_init->GetBinContent(j,i,k);
	Sq_sum += r3numSq_init->GetBinContent(j,k,i) + r3numSq_init->GetBinContent(k,i,j) + r3numSq_init->GetBinContent(k,j,i);
	double den_sum = r3numEn_init->GetBinContent(i,j,k) + r3numEn_init->GetBinContent(i,k,j) + r3numEn_init->GetBinContent(j,i,k);
	den_sum += r3numEn_init->GetBinContent(j,k,i) + r3numEn_init->GetBinContent(k,i,j) + r3numEn_init->GetBinContent(k,j,i);
	if(den_sum>0) {
	  r3num->SetBinContent(i,j,k, num_sum/den_sum); r3num->SetBinContent(i,k,j, num_sum/den_sum); 
	  r3num->SetBinContent(j,i,k, num_sum/den_sum); r3num->SetBinContent(j,k,i, num_sum/den_sum);
	  r3num->SetBinContent(k,i,j, num_sum/den_sum); r3num->SetBinContent(k,j,i, num_sum/den_sum);
	  r3numSq->SetBinContent(i,j,k, Sq_sum/den_sum); r3numSq->SetBinContent(i,k,j, Sq_sum/den_sum); 
	  r3numSq->SetBinContent(j,i,k, Sq_sum/den_sum); r3numSq->SetBinContent(j,k,i, Sq_sum/den_sum);
	  r3numSq->SetBinContent(k,i,j, Sq_sum/den_sum); r3numSq->SetBinContent(k,j,i, Sq_sum/den_sum);
	  r3numEn->SetBinContent(i,j,k, den_sum); r3numEn->SetBinContent(i,k,j, den_sum); 
	  r3numEn->SetBinContent(j,i,k, den_sum); r3numEn->SetBinContent(j,k,i, den_sum);
	  r3numEn->SetBinContent(k,i,j, den_sum); r3numEn->SetBinContent(k,j,i, den_sum);
	}else {
	  r3num->SetBinContent(i,j,k, 0.); r3num->SetBinContent(i,k,j, 0.);
	  r3num->SetBinContent(j,i,k, 0.); r3num->SetBinContent(j,k,i, 0.);
	  r3num->SetBinContent(k,i,j, 0.); r3num->SetBinContent(k,j,i, 0.);
	  r3numEn->SetBinContent(i,j,k, 0.); r3numEn->SetBinContent(i,k,j, 0.); 
	  r3numEn->SetBinContent(j,i,k, 0.); r3numEn->SetBinContent(j,k,i, 0.);
	  r3numEn->SetBinContent(k,i,j, 0.); r3numEn->SetBinContent(k,j,i, 0.);
	}
	// 2-pion phases
	// 12
	num_sum = r3den1_init->GetBinContent(i,j,k) + r3den1_init->GetBinContent(i,k,j);
	Sq_sum = r3den1Sq_init->GetBinContent(i,j,k) + r3den1Sq_init->GetBinContent(i,k,j);
	den_sum = r3numEn_init->GetBinContent(i,j,k) + r3numEn_init->GetBinContent(i,k,j);
	if(den_sum>0) {
	  r3den1->SetBinContent(i,j,k, num_sum/den_sum); r3den1->SetBinContent(i,k,j, num_sum/den_sum); 
	  r3den1Sq->SetBinContent(i,j,k, Sq_sum/den_sum); r3den1Sq->SetBinContent(i,k,j, Sq_sum/den_sum);
	  r3den1En->SetBinContent(i,j,k, den_sum); r3den1En->SetBinContent(i,k,j, den_sum); 
	}else{
	  r3den1->SetBinContent(i,j,k, 0.); r3den1->SetBinContent(i,k,j, 0.); 
	  r3den1Sq->SetBinContent(i,j,k, 0.); r3den1Sq->SetBinContent(i,k,j, 0.);
	  r3den1En->SetBinContent(i,j,k, 0.); r3den1En->SetBinContent(i,k,j, 0.);
	}
	// 23
	num_sum = r3den2_init->GetBinContent(j,i,k) + r3den2_init->GetBinContent(k,i,j);
	Sq_sum = r3den2Sq_init->GetBinContent(j,i,k) + r3den2Sq_init->GetBinContent(k,i,j);
	den_sum = r3numEn_init->GetBinContent(j,i,k) + r3numEn_init->GetBinContent(k,i,j);
	if(den_sum>0) {
	  r3den2->SetBinContent(j,i,k, num_sum/den_sum); r3den2->SetBinContent(k,i,j, num_sum/den_sum); 
	  r3den2Sq->SetBinContent(j,i,k, Sq_sum/den_sum); r3den2Sq->SetBinContent(k,i,j, Sq_sum/den_sum);
	  r3den2En->SetBinContent(j,i,k, den_sum); r3den2En->SetBinContent(k,i,j, den_sum);
	}else{
	  r3den2->SetBinContent(j,i,k, 0.); r3den2->SetBinContent(k,i,j, 0.); 
	  r3den2Sq->SetBinContent(j,i,k, 0.); r3den2Sq->SetBinContent(k,i,j, 0.);
	  r3den2En->SetBinContent(j,i,k, 0.); r3den2En->SetBinContent(k,i,j, 0.);
	}
	// 12
	num_sum = r3den3_init->GetBinContent(j,k,i) + r3den3_init->GetBinContent(k,j,i);
	Sq_sum = r3den3Sq_init->GetBinContent(j,k,i) + r3den3Sq_init->GetBinContent(k,j,i);
	den_sum = r3numEn_init->GetBinContent(j,k,i) + r3numEn_init->GetBinContent(k,j,i);
	if(den_sum>0) {
	  r3den3->SetBinContent(j,k,i, num_sum/den_sum); r3den3->SetBinContent(k,j,i, num_sum/den_sum); 
	  r3den3Sq->SetBinContent(j,k,i, Sq_sum/den_sum); r3den3Sq->SetBinContent(k,j,i, Sq_sum/den_sum);
	  r3den3En->SetBinContent(j,k,i, den_sum); r3den3En->SetBinContent(k,j,i, den_sum);
	}else{
	  r3den3->SetBinContent(j,k,i, 0.); r3den3->SetBinContent(k,j,i, 0.); 
	  r3den3Sq->SetBinContent(j,k,i, 0.); r3den3Sq->SetBinContent(k,j,i, 0.);
	  r3den3En->SetBinContent(j,k,i, 0.); r3den3En->SetBinContent(k,j,i, 0.);
	}
	
      }
    }
  }
  cout<<"Done Symmetrizing"<<endl;
  */

  r3num->Sumw2();// r3numEn->Sumw2();
  r3den1->Sumw2();// r3den1En->Sumw2();
  r3den2->Sumw2();// r3den2En->Sumw2();
  r3den3->Sumw2();// r3den3En->Sumw2();


  r3num->Divide(r3numEn);
  r3den1->Divide(r3den1En);
  r3den2->Divide(r3den2En);
  r3den3->Divide(r3den3En);
  /*r3num->Divide(r3num,r3numEn,1,1,"B");
  r3den1->Divide(r3den1,r3den1En,1,1,"B");
  r3den2->Divide(r3den2,r3den2En,1,1,"B");
  r3den3->Divide(r3den3,r3den3En,1,1,"B");*/
  r3numSq->Divide(r3numEn);
  r3den1Sq->Divide(r3den1En);
  r3den2Sq->Divide(r3den2En);
  r3den3Sq->Divide(r3den3En);
  

  TH1D *r3=new TH1D("r3","",10,0,0.1);
  TH1D *r3den=new TH1D("r3den","",10,0,0.1);
  double r3num_e[20]={0};
  double r3den_e[20]={0};
  double NegativeDenCount[20]={0};
  double TotalDenCount[20]={0};
  double r3den_neg[20]={0};
  double r3den_neg_Sq[20]={0};

  for(int i=1; i<=50; i++){
    for(int j=1; j<=50; j++){
      for(int k=1; k<=50; k++){

	double q3 =sqrt(pow(0.002*(i-0.5),2)+pow(0.002*(j-0.5),2)+pow(0.002*(k-0.5),2));
	if(q3>=0.2) continue;

	if(r3numEn->GetBinContent(i,j,k)==0) continue;
	if(r3den1En->GetBinContent(i,j,k)==0) continue;
	if(r3den2En->GetBinContent(i,j,k)==0) continue;
	if(r3den3En->GetBinContent(i,j,k)==0) continue;

	
	double num = r3num->GetBinContent(i,j,k);
	num *= pow(epsilon,0.5) * (3-2*epsilon)/pow(2-epsilon,1.5);
	double den1 = r3den1->GetBinContent(i,j,k);
	double den2 = r3den2->GetBinContent(i,j,k);
	double den3 = r3den3->GetBinContent(i,j,k);
	double den = den1*den2*den3;
	
	TotalDenCount[int(q3/0.01)]++;
	if(den <= 0) {
	  NegativeDenCount[int(q3/0.01)]++;
	  r3den_neg[int(q3/0.01)] += sqrt(fabs(den));
	  r3den_neg_Sq[int(q3/0.01)] += fabs(den);
	  continue;
	}
	den = sqrt(den);
	double weightWidth = sqrt( (r3den1Sq->GetBinContent(i,j,k) - pow(r3den1->GetBinContent(i,j,k),2))/r3den1En->GetBinContent(i,j,k));
	double den_e = pow(0.5*weightWidth*den2*den3/den,2);
	weightWidth = sqrt( fabs(r3den2Sq->GetBinContent(i,j,k) - pow(r3den2->GetBinContent(i,j,k),2))/r3den2En->GetBinContent(i,j,k));
	den_e += pow(0.5*weightWidth*den1*den3/den,2);
	weightWidth = sqrt( fabs(r3den3Sq->GetBinContent(i,j,k) - pow(r3den3->GetBinContent(i,j,k),2))/r3den3En->GetBinContent(i,j,k));
	den_e += pow(0.5*weightWidth*den1*den2/den,2);
	//
	r3->Fill(q3,num);
	r3den->Fill(q3,den);
	//
	double num_e = (r3numSq->GetBinContent(i,j,k)-pow(r3num->GetBinContent(i,j,k),2))/r3numEn->GetBinContent(i,j,k);
	r3num_e[int(q3/0.01)] += num_e;
	r3den_e[int(q3/0.01)] += den_e;
	

	//if(q3<0.01) cout<<i<<j<<k<<"   "<<r3num_e[0]<<"  "<<r3den_e[0]<<"  "<<r3numSq->GetBinContent(i,j,k)<<"  "<<pow(r3num->GetBinContent(i,j,k),2)<<"  "<<r3numEn->GetBinContent(i,j,k)<<endl;
	//if(q3<0.02) cout<<r3den1->GetBinContent(i,j,k)<<"  "<<r3den1->GetBinContent(i,k,j)<<endl;
      }
    }
  }

  cout<<"Done calculating r3"<<endl;

  TCanvas *can1 = new TCanvas("can1", "can1",10,0,600,600);// 11,53,700,500
  can1->SetHighLightColor(2);
  //can1->Range(-0.7838115,-1.033258,0.7838115,1.033258);
  can1->SetFillColor(10);//10
  can1->SetBorderMode(0);
  can1->SetBorderSize(2);
  can1->SetFrameFillColor(0);
  can1->SetFrameBorderMode(0);
  can1->SetFrameBorderMode(0);
  
  
  TPad *pad1 = new TPad("pad1","pad1",0,0,1.,1.);
  gPad->SetGridx(0);
  gPad->SetGridy(0);
  gPad->SetTickx();
  gPad->SetTicky();
  pad1->SetTopMargin(0.02);
  pad1->SetRightMargin(0.03);
  pad1->SetBottomMargin(0.06);//0.12
  pad1->Draw();
  
  TLegend *legend1 = new TLegend(.15,.25,.65,.45,NULL,"brNDC");
  legend1->SetBorderSize(1);
  legend1->SetTextSize(.04);// small .03; large .036 
  //legend1->SetLineColor(0);
  legend1->SetFillColor(0);

  TF1 *ChaoticLimit=new TF1("ChaoticLimit","2",0,0.2);
  ChaoticLimit->SetLineStyle(3);
  TF1 *Quartic=new TF1("Quartic","[0]*(1-[1]*pow(x,4))",0,0.1);
  Quartic->SetLineColor(4);
  TF1 *Quadratic=new TF1("Quadratic","[0]*(1-[1]*pow(x,2))",0,0.1);
  //Quadratic->SetLineStyle(2);
  Quadratic->SetLineColor(6);

  
  for(int i=0; i<10; i++){
    cout<<"Lost Triplet Fraction: "<<NegativeDenCount[i]/TotalDenCount[i]<<endl;
    r3->SetBinError(i+1,sqrt(r3num_e[i]));
    double ExtraDenError = 0;
    if(NegativeDenCount[i]>0) ExtraDenError = sqrt( fabs(r3den_neg_Sq[i]/NegativeDenCount[i] - pow(r3den_neg[i]/NegativeDenCount[i],2))/NegativeDenCount[i]) * NegativeDenCount[i];
    r3den->SetBinContent(i+1, r3den->GetBinContent(i+1) - r3den_neg[i]);
    r3den->SetBinError(i+1,sqrt(r3den_e[i])+r3den_neg[i]);
  }
  //r3->Sumw2(); r3den->Sumw2();
  //r3->Divide(r3, r3den, 1,1,"B");
  r3->Divide(r3den);
  

  // first bin sometimes only has 1 entry (bad error bar)
  if(r3->GetBinError(1)<0.001) r3->SetBinError(1, r3->GetBinError(2));
  
  gPad->SetLeftMargin(0.12); gPad->SetRightMargin(0.03); gPad->SetBottomMargin(0.1); gPad->SetTopMargin(0.01);
  r3->SetMarkerStyle(24);
  r3->SetMarkerColor(1);
  r3->SetLineColor(1);
  r3->GetXaxis()->SetLabelFont(63); r3->GetYaxis()->SetLabelFont(63);
  r3->GetXaxis()->SetLabelSize(18); r3->GetYaxis()->SetLabelSize(18);
  r3->GetXaxis()->SetTitleFont(63); r3->GetYaxis()->SetTitleFont(63);
  r3->GetXaxis()->SetTitleSize(32); r3->GetYaxis()->SetTitleSize(32);
  r3->GetXaxis()->SetTitleOffset(0.75); r3->GetYaxis()->SetTitleOffset(1.);
  r3->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
  r3->GetYaxis()->SetTitle("r_{3}");
  r3->SetTitle("");
  r3->SetMinimum(1.1); r3->SetMaximum(2.5);
  r3->GetXaxis()->SetRangeUser(0,0.069);
  r3->Draw();
 
  //r3->Fit(Quartic,"IME","",0,0.09);
  //r3->Fit(Quadratic,"IME","",0,0.09);
  ChaoticLimit->Draw("same");
  //Quartic->Draw("same");
  //Quadratic->Draw("same");
  //legend1->AddEntry(Quartic,"Quartic: #chi^{2}/NDF=1.2","l");
  //legend1->AddEntry(Quadratic,"Quadratic: #chi^{2}/NDF=5.4","l");
  //legend1->Draw("same");
  
  //cout<<"Quartic Chi2/NDF = "<<Quartic->GetChisquare()/Quartic->GetNDF()<<endl;
  //cout<<"Quadratic Chi2/NDF = "<<Quadratic->GetChisquare()/Quadratic->GetNDF()<<endl;
  
  TLatex *Specif = new TLatex(0.005,1.4,"Therminator");// Therminator
  Specif->SetTextSize(.075);
  Specif->Draw("same");
  
}
Commit	Line	Data
9129699a	1	#include <math.h>
	2	#include <time.h>
	3	#include <stdio.h>
	4	#include <stdlib.h>
	5	#include <Riostream.h>
	6	#include <complex>
	7
	8	#include "TObject.h"
	9	#include "TTree.h"
	10	#include "TBranch.h"
	11	#include "TLeaf.h"
	12	#include "TVector2.h"
	13	#include "TVector3.h"
	14	#include "TFile.h"
	15	#include "TString.h"
	16	#include "TF1.h"
	17	#include "TH1.h"
	18	#include "TH2.h"
	19	#include "TH3.h"
	20	#include "TProfile.h"
	21	#include "TProfile2D.h"
	22	#include "TMath.h"
	23	#include "TText.h"
	24	#include "TRandom3.h"
	25	#include "TArray.h"
	26	#include "TLegend.h"
	27	#include "TStyle.h"
	28	#include "TMinuit.h"
	29	#include "TChain.h"
	30	#include "TLorentzVector.h"
	31	#include "TPaveStats.h"
	32	#include "TLatex.h"
	33	#include "TCanvas.h"
	34	#include "TPad.h"
	35
	36	#define PI 3.1415926
	37	#define BohrR 1963.6885 // Mate's value(1963.6885) ~ 387.5/0.197327(1963.7455)
	38	#define FmToGeV 0.197327 // 0.197327
	39
	40	using namespace std;
	41
	42	void PlotTherm_r3(){
	43
	44	double epsilon=1.0;// chaoticity
	45
	46	gStyle->SetOptStat(0);
	47	//gStyle->SetOptFit(111);
	48
	49	//TFile infile=new TFile("Therm_r3_b9.root","READ");// For v2 paper (r<100 and r<80)
	50	TFile *infile=new TFile("Therm_FSI_b2.root","READ");
	51	//TFile *infile=new TFile("Therm_dist_temp.root","READ");
	52
	53	//
	54	TH3D r3num_init=(TH3D)infile->Get("r3num");
	55	TH3D r3den1_init=(TH3D)infile->Get("r3den1");
	56	TH3D r3den2_init=(TH3D)infile->Get("r3den2");
	57	TH3D r3den3_init=(TH3D)infile->Get("r3den3");
	58	TH3D r3numSq_init=(TH3D)infile->Get("r3numSq");
	59	TH3D r3den1Sq_init=(TH3D)infile->Get("r3den1Sq");
	60	TH3D r3den2Sq_init=(TH3D)infile->Get("r3den2Sq");
	61	TH3D r3den3Sq_init=(TH3D)infile->Get("r3den3Sq");
	62	TH3D r3numEn_init=(TH3D)infile->Get("r3numEn");
	63	TH3D r3den1En_init=(TH3D)infile->Get("r3den1En");
	64	TH3D r3den2En_init=(TH3D)infile->Get("r3den2En");
65	TH3D r3den3En_init=(TH3D)infile->Get("r3den3En");
66	TH3D r3numME_init=(TH3D)infile->Get("r3numME");
67	TH3D r3den1ME_init=(TH3D)infile->Get("r3den1ME");
68	TH3D r3den2ME_init=(TH3D)infile->Get("r3den2ME");
69	TH3D r3den3ME_init=(TH3D)infile->Get("r3den3ME");
70
71	TH2D R2_2D = (TH2D)infile->Get("Num_Cos_ss");
72	//TH2D R2Sq_2D = (TH2D)infile->Get("NumSq_Cos_ss");
73	TH2D R2En_2D = (TH2D)infile->Get("Den_ss");
74	TH1D R2=(TH1D)R2_2D->ProjectionY();
75	//TH1D R2Sq=(TH1D)R2Sq_2D->ProjectionY();
76	TH1D R2En=(TH1D)R2En_2D->ProjectionY();
77	R2->Divide(R2En);
78	//R2Sq->Divide(R2En);
79
80	TH3D r3num=(TH3D)r3num_init->Clone();
81	TH3D r3den1=(TH3D)r3den1_init->Clone();
82	TH3D r3den2=(TH3D)r3den2_init->Clone();
83	TH3D r3den3=(TH3D)r3den3_init->Clone();
84	TH3D r3numSq=(TH3D)r3numSq_init->Clone();
85	TH3D r3den1Sq=(TH3D)r3den1Sq_init->Clone();
86	TH3D r3den2Sq=(TH3D)r3den2Sq_init->Clone();
87	TH3D r3den3Sq=(TH3D)r3den3Sq_init->Clone();
88	TH3D r3numEn=(TH3D)r3numEn_init->Clone();
89	TH3D r3den1En=(TH3D)r3den1En_init->Clone();
90	TH3D r3den2En=(TH3D)r3den2En_init->Clone();
91	TH3D r3den3En=(TH3D)r3den3En_init->Clone();
92	TH3D r3numME=(TH3D)r3numME_init->Clone();
93	TH3D r3den1ME=(TH3D)r3den1ME_init->Clone();
94	TH3D r3den2ME=(TH3D)r3den2ME_init->Clone();
95	TH3D r3den3ME=(TH3D)r3den3ME_init->Clone();
96
97
98	int FB=40;
99	int LB=50;
100	/*double SF = r3numEn->Integral(FB,LB,FB,LB,FB,LB)/r3numME->Integral(FB,LB,FB,LB,FB,LB);
101	cout<<SF<<endl;
102	r3numME->Scale(SF);
103	SF = r3den1En->Integral(FB,LB,FB,LB,FB,LB)/r3den1ME->Integral(FB,LB,FB,LB,FB,LB);
104	cout<<SF<<endl;
105	r3den1ME->Scale(SF);
106	SF = r3den2En->Integral(FB,LB,FB,LB,FB,LB)/r3den2ME->Integral(FB,LB,FB,LB,FB,LB);
107	cout<<SF<<endl;
108	r3den2ME->Scale(SF);
109	SF = r3den3En->Integral(FB,LB,FB,LB,FB,LB)/r3den3ME->Integral(FB,LB,FB,LB,FB,LB);
110	cout<<SF<<endl;
111	r3den3ME->Scale(SF);*/
112	//
113	/*
114	//////////////////////
115	// Symmetrize bin contents
116	for(int i=1; i<=50; i++){
117	for(int j=1; j<=50; j++){
118	for(int k=1; k<=50; k++){
119	// 3-pion phase
120	double num_sum = r3num_init->GetBinContent(i,j,k) + r3num_init->GetBinContent(i,k,j) + r3num_init->GetBinContent(j,i,k);
121	num_sum += r3num_init->GetBinContent(j,k,i) + r3num_init->GetBinContent(k,i,j) + r3num_init->GetBinContent(k,j,i);
122	double Sq_sum = r3numSq_init->GetBinContent(i,j,k) + r3numSq_init->GetBinContent(i,k,j) + r3numSq_init->GetBinContent(j,i,k);
123	Sq_sum += r3numSq_init->GetBinContent(j,k,i) + r3numSq_init->GetBinContent(k,i,j) + r3numSq_init->GetBinContent(k,j,i);
124	double den_sum = r3numEn_init->GetBinContent(i,j,k) + r3numEn_init->GetBinContent(i,k,j) + r3numEn_init->GetBinContent(j,i,k);
125	den_sum += r3numEn_init->GetBinContent(j,k,i) + r3numEn_init->GetBinContent(k,i,j) + r3numEn_init->GetBinContent(k,j,i);
126	if(den_sum>0) {
127	r3num->SetBinContent(i,j,k, num_sum/den_sum); r3num->SetBinContent(i,k,j, num_sum/den_sum);
128	r3num->SetBinContent(j,i,k, num_sum/den_sum); r3num->SetBinContent(j,k,i, num_sum/den_sum);
129	r3num->SetBinContent(k,i,j, num_sum/den_sum); r3num->SetBinContent(k,j,i, num_sum/den_sum);
130	r3numSq->SetBinContent(i,j,k, Sq_sum/den_sum); r3numSq->SetBinContent(i,k,j, Sq_sum/den_sum);
131	r3numSq->SetBinContent(j,i,k, Sq_sum/den_sum); r3numSq->SetBinContent(j,k,i, Sq_sum/den_sum);
132	r3numSq->SetBinContent(k,i,j, Sq_sum/den_sum); r3numSq->SetBinContent(k,j,i, Sq_sum/den_sum);
133	r3numEn->SetBinContent(i,j,k, den_sum); r3numEn->SetBinContent(i,k,j, den_sum);
134	r3numEn->SetBinContent(j,i,k, den_sum); r3numEn->SetBinContent(j,k,i, den_sum);
135	r3numEn->SetBinContent(k,i,j, den_sum); r3numEn->SetBinContent(k,j,i, den_sum);
136	}else {
137	r3num->SetBinContent(i,j,k, 0.); r3num->SetBinContent(i,k,j, 0.);
138	r3num->SetBinContent(j,i,k, 0.); r3num->SetBinContent(j,k,i, 0.);
139	r3num->SetBinContent(k,i,j, 0.); r3num->SetBinContent(k,j,i, 0.);
140	r3numEn->SetBinContent(i,j,k, 0.); r3numEn->SetBinContent(i,k,j, 0.);
141	r3numEn->SetBinContent(j,i,k, 0.); r3numEn->SetBinContent(j,k,i, 0.);
142	r3numEn->SetBinContent(k,i,j, 0.); r3numEn->SetBinContent(k,j,i, 0.);
143	}
144	// 2-pion phases
145	// 12
146	num_sum = r3den1_init->GetBinContent(i,j,k) + r3den1_init->GetBinContent(i,k,j);
147	Sq_sum = r3den1Sq_init->GetBinContent(i,j,k) + r3den1Sq_init->GetBinContent(i,k,j);
148	den_sum = r3numEn_init->GetBinContent(i,j,k) + r3numEn_init->GetBinContent(i,k,j);
149	if(den_sum>0) {
150	r3den1->SetBinContent(i,j,k, num_sum/den_sum); r3den1->SetBinContent(i,k,j, num_sum/den_sum);
151	r3den1Sq->SetBinContent(i,j,k, Sq_sum/den_sum); r3den1Sq->SetBinContent(i,k,j, Sq_sum/den_sum);
152	r3den1En->SetBinContent(i,j,k, den_sum); r3den1En->SetBinContent(i,k,j, den_sum);
153	}else{
154	r3den1->SetBinContent(i,j,k, 0.); r3den1->SetBinContent(i,k,j, 0.);
155	r3den1Sq->SetBinContent(i,j,k, 0.); r3den1Sq->SetBinContent(i,k,j, 0.);
156	r3den1En->SetBinContent(i,j,k, 0.); r3den1En->SetBinContent(i,k,j, 0.);
157	}
158	// 23
159	num_sum = r3den2_init->GetBinContent(j,i,k) + r3den2_init->GetBinContent(k,i,j);
160	Sq_sum = r3den2Sq_init->GetBinContent(j,i,k) + r3den2Sq_init->GetBinContent(k,i,j);
161	den_sum = r3numEn_init->GetBinContent(j,i,k) + r3numEn_init->GetBinContent(k,i,j);
162	if(den_sum>0) {
163	r3den2->SetBinContent(j,i,k, num_sum/den_sum); r3den2->SetBinContent(k,i,j, num_sum/den_sum);
164	r3den2Sq->SetBinContent(j,i,k, Sq_sum/den_sum); r3den2Sq->SetBinContent(k,i,j, Sq_sum/den_sum);
165	r3den2En->SetBinContent(j,i,k, den_sum); r3den2En->SetBinContent(k,i,j, den_sum);
166	}else{
167	r3den2->SetBinContent(j,i,k, 0.); r3den2->SetBinContent(k,i,j, 0.);
168	r3den2Sq->SetBinContent(j,i,k, 0.); r3den2Sq->SetBinContent(k,i,j, 0.);
169	r3den2En->SetBinContent(j,i,k, 0.); r3den2En->SetBinContent(k,i,j, 0.);
170	}
171	// 12
172	num_sum = r3den3_init->GetBinContent(j,k,i) + r3den3_init->GetBinContent(k,j,i);
173	Sq_sum = r3den3Sq_init->GetBinContent(j,k,i) + r3den3Sq_init->GetBinContent(k,j,i);
174	den_sum = r3numEn_init->GetBinContent(j,k,i) + r3numEn_init->GetBinContent(k,j,i);
175	if(den_sum>0) {
176	r3den3->SetBinContent(j,k,i, num_sum/den_sum); r3den3->SetBinContent(k,j,i, num_sum/den_sum);
177	r3den3Sq->SetBinContent(j,k,i, Sq_sum/den_sum); r3den3Sq->SetBinContent(k,j,i, Sq_sum/den_sum);
178	r3den3En->SetBinContent(j,k,i, den_sum); r3den3En->SetBinContent(k,j,i, den_sum);
179	}else{
180	r3den3->SetBinContent(j,k,i, 0.); r3den3->SetBinContent(k,j,i, 0.);
181	r3den3Sq->SetBinContent(j,k,i, 0.); r3den3Sq->SetBinContent(k,j,i, 0.);
182	r3den3En->SetBinContent(j,k,i, 0.); r3den3En->SetBinContent(k,j,i, 0.);
183	}
184
185	}
186	}
187	}
188	cout<<"Done Symmetrizing"<<endl;
189	*/
190
191	r3num->Sumw2();// r3numEn->Sumw2();
192	r3den1->Sumw2();// r3den1En->Sumw2();
193	r3den2->Sumw2();// r3den2En->Sumw2();
194	r3den3->Sumw2();// r3den3En->Sumw2();
195
196
197
198
199	r3num->Divide(r3numEn);
200	r3den1->Divide(r3den1En);
201	r3den2->Divide(r3den2En);
202	r3den3->Divide(r3den3En);
203	/*r3num->Divide(r3num,r3numEn,1,1,"B");
204	r3den1->Divide(r3den1,r3den1En,1,1,"B");
205	r3den2->Divide(r3den2,r3den2En,1,1,"B");
206	r3den3->Divide(r3den3,r3den3En,1,1,"B");*/
207	r3numSq->Divide(r3numEn);
208	r3den1Sq->Divide(r3den1En);
209	r3den2Sq->Divide(r3den2En);
210	r3den3Sq->Divide(r3den3En);
211
212
213
214
215	TH1D *r3=new TH1D("r3","",10,0,0.1);
216	TH1D *r3den=new TH1D("r3den","",10,0,0.1);
217	double r3num_e[20]={0};
218	double r3den_e[20]={0};
219	double NegativeDenCount[20]={0};
220	double TotalDenCount[20]={0};
221	double r3den_neg[20]={0};
222	double r3den_neg_Sq[20]={0};
223
224	for(int i=1; i<=50; i++){
225	for(int j=1; j<=50; j++){
226	for(int k=1; k<=50; k++){
227
228	double q3 =sqrt(pow(0.002(i-0.5),2)+pow(0.002(j-0.5),2)+pow(0.002*(k-0.5),2));
229	if(q3>=0.2) continue;
230
231	if(r3numEn->GetBinContent(i,j,k)==0) continue;
232	if(r3den1En->GetBinContent(i,j,k)==0) continue;
233	if(r3den2En->GetBinContent(i,j,k)==0) continue;
234	if(r3den3En->GetBinContent(i,j,k)==0) continue;
235
236
237	double num = r3num->GetBinContent(i,j,k);
238	num = pow(epsilon,0.5) (3-2*epsilon)/pow(2-epsilon,1.5);
239	double den1 = r3den1->GetBinContent(i,j,k);
240	double den2 = r3den2->GetBinContent(i,j,k);
241	double den3 = r3den3->GetBinContent(i,j,k);
242	double den = den1den2den3;
243
244	TotalDenCount[int(q3/0.01)]++;
245	if(den <= 0) {
246	NegativeDenCount[int(q3/0.01)]++;
247	r3den_neg[int(q3/0.01)] += sqrt(fabs(den));
248	r3den_neg_Sq[int(q3/0.01)] += fabs(den);
249	continue;
250	}
251	den = sqrt(den);
252	double weightWidth = sqrt( (r3den1Sq->GetBinContent(i,j,k) - pow(r3den1->GetBinContent(i,j,k),2))/r3den1En->GetBinContent(i,j,k));
253	double den_e = pow(0.5weightWidthden2*den3/den,2);
254	weightWidth = sqrt( fabs(r3den2Sq->GetBinContent(i,j,k) - pow(r3den2->GetBinContent(i,j,k),2))/r3den2En->GetBinContent(i,j,k));
255	den_e += pow(0.5weightWidthden1*den3/den,2);
256	weightWidth = sqrt( fabs(r3den3Sq->GetBinContent(i,j,k) - pow(r3den3->GetBinContent(i,j,k),2))/r3den3En->GetBinContent(i,j,k));
257	den_e += pow(0.5weightWidthden1*den2/den,2);
258	//
259	r3->Fill(q3,num);
260	r3den->Fill(q3,den);
261	//
262	double num_e = (r3numSq->GetBinContent(i,j,k)-pow(r3num->GetBinContent(i,j,k),2))/r3numEn->GetBinContent(i,j,k);
263	r3num_e[int(q3/0.01)] += num_e;
264	r3den_e[int(q3/0.01)] += den_e;
265
266
267	//if(q3<0.01) cout<<i<<j<<k<<" "<<r3num_e[0]<<" "<<r3den_e[0]<<" "<<r3numSq->GetBinContent(i,j,k)<<" "<<pow(r3num->GetBinContent(i,j,k),2)<<" "<<r3numEn->GetBinContent(i,j,k)<<endl;
268	//if(q3<0.02) cout<<r3den1->GetBinContent(i,j,k)<<" "<<r3den1->GetBinContent(i,k,j)<<endl;
269	}
270	}
271	}
272
273	cout<<"Done calculating r3"<<endl;
274
275	TCanvas *can1 = new TCanvas("can1", "can1",10,0,600,600);// 11,53,700,500
276	can1->SetHighLightColor(2);
277	//can1->Range(-0.7838115,-1.033258,0.7838115,1.033258);
278	can1->SetFillColor(10);//10
279	can1->SetBorderMode(0);
280	can1->SetBorderSize(2);
281	can1->SetFrameFillColor(0);
282	can1->SetFrameBorderMode(0);
283	can1->SetFrameBorderMode(0);
284
285
286	TPad *pad1 = new TPad("pad1","pad1",0,0,1.,1.);
287	gPad->SetGridx(0);
288	gPad->SetGridy(0);
289	gPad->SetTickx();
290	gPad->SetTicky();
291	pad1->SetTopMargin(0.02);
292	pad1->SetRightMargin(0.03);
293	pad1->SetBottomMargin(0.06);//0.12
294	pad1->Draw();
295
296	TLegend *legend1 = new TLegend(.15,.25,.65,.45,NULL,"brNDC");
297	legend1->SetBorderSize(1);
298	legend1->SetTextSize(.04);// small .03; large .036
299	//legend1->SetLineColor(0);
300	legend1->SetFillColor(0);
301
302	TF1 *ChaoticLimit=new TF1("ChaoticLimit","2",0,0.2);
303	ChaoticLimit->SetLineStyle(3);
304	TF1 Quartic=new TF1("Quartic","[0](1-[1]*pow(x,4))",0,0.1);
305	Quartic->SetLineColor(4);
306	TF1 Quadratic=new TF1("Quadratic","[0](1-[1]*pow(x,2))",0,0.1);
307	//Quadratic->SetLineStyle(2);
308	Quadratic->SetLineColor(6);
309
310
311
312	for(int i=0; i<10; i++){
313	cout<<"Lost Triplet Fraction: "<<NegativeDenCount[i]/TotalDenCount[i]<<endl;
314	r3->SetBinError(i+1,sqrt(r3num_e[i]));
315	double ExtraDenError = 0;
316	if(NegativeDenCount[i]>0) ExtraDenError = sqrt( fabs(r3den_neg_Sq[i]/NegativeDenCount[i] - pow(r3den_neg[i]/NegativeDenCount[i],2))/NegativeDenCount[i]) * NegativeDenCount[i];
317	r3den->SetBinContent(i+1, r3den->GetBinContent(i+1) - r3den_neg[i]);
318	r3den->SetBinError(i+1,sqrt(r3den_e[i])+r3den_neg[i]);
319	}
320	//r3->Sumw2(); r3den->Sumw2();
321	//r3->Divide(r3, r3den, 1,1,"B");
322	r3->Divide(r3den);
323
324
325	// first bin sometimes only has 1 entry (bad error bar)
326	if(r3->GetBinError(1)<0.001) r3->SetBinError(1, r3->GetBinError(2));
327
328	gPad->SetLeftMargin(0.12); gPad->SetRightMargin(0.03); gPad->SetBottomMargin(0.1); gPad->SetTopMargin(0.01);
329	r3->SetMarkerStyle(24);
330	r3->SetMarkerColor(1);
331	r3->SetLineColor(1);
332	r3->GetXaxis()->SetLabelFont(63); r3->GetYaxis()->SetLabelFont(63);
333	r3->GetXaxis()->SetLabelSize(18); r3->GetYaxis()->SetLabelSize(18);
334	r3->GetXaxis()->SetTitleFont(63); r3->GetYaxis()->SetTitleFont(63);
335	r3->GetXaxis()->SetTitleSize(32); r3->GetYaxis()->SetTitleSize(32);
336	r3->GetXaxis()->SetTitleOffset(0.75); r3->GetYaxis()->SetTitleOffset(1.);
337	r3->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
338	r3->GetYaxis()->SetTitle("r_{3}");
339	r3->SetTitle("");
340	r3->SetMinimum(1.1); r3->SetMaximum(2.5);
341	r3->GetXaxis()->SetRangeUser(0,0.069);
342	r3->Draw();
343
344	//r3->Fit(Quartic,"IME","",0,0.09);
345	//r3->Fit(Quadratic,"IME","",0,0.09);
346	ChaoticLimit->Draw("same");
347	//Quartic->Draw("same");
348	//Quadratic->Draw("same");
349	//legend1->AddEntry(Quartic,"Quartic: #chi^{2}/NDF=1.2","l");
350	//legend1->AddEntry(Quadratic,"Quadratic: #chi^{2}/NDF=5.4","l");
351	//legend1->Draw("same");
352
353	//cout<<"Quartic Chi2/NDF = "<<Quartic->GetChisquare()/Quartic->GetNDF()<<endl;
354	//cout<<"Quadratic Chi2/NDF = "<<Quadratic->GetChisquare()/Quadratic->GetNDF()<<endl;
355
356	TLatex *Specif = new TLatex(0.005,1.4,"Therminator");// Therminator
357	Specif->SetTextSize(.075);
358	Specif->Draw("same");
359
360	}