--- /dev/null
+#include <math.h>
+#include <time.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <Riostream.h>
+#include <complex>
+
+#include "TVector2.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"
+
+using namespace std;
+
+void MakeKFileTherm(){
+ // 2-particles
+
+ const int kbVALUES=6;// 6 b values (2,3,5,7,8,9)
+
+ TFile *probeFile=new TFile("Therm_FSI_b2.root","READ");// all files have same binning.
+ TH2D *probeHisto = (TH2D*)probeFile->Get("K2_ss");// kt x qinv
+ int binsY= probeHisto->GetNbinsY();
+ double lowY = probeHisto->GetYaxis()->GetBinLowEdge(1);
+ double highY = probeHisto->GetYaxis()->GetBinUpEdge(binsY);
+ probeFile->Close();
+
+ TFile *OutFile=new TFile("KFile_temp.root","RECREATE");
+ TH2D *K2ssT = new TH2D("K2ssT","",kbVALUES,0.5,kbVALUES+0.5, binsY,lowY,highY);// kt integrated
+ TH2D *K2osT = new TH2D("K2osT","",kbVALUES,0.5,kbVALUES+0.5, binsY,lowY,highY);// kt integrated
+ TH3D *K2ssT_kt = new TH3D("K2ssT_kt","",kbVALUES,0.5,kbVALUES+0.5, 6,0.5,6+0.5, binsY,lowY,highY);// kt differential
+ TH3D *K2osT_kt = new TH3D("K2osT_kt","",kbVALUES,0.5,kbVALUES+0.5, 6,0.5,6+0.5, binsY,lowY,highY);// kt differential
+ K2ssT_kt->GetXaxis()->SetTitle("b bin"); K2osT_kt->GetXaxis()->SetTitle("b bin");
+ K2ssT_kt->GetYaxis()->SetTitle("kt bin"); K2osT_kt->GetYaxis()->SetTitle("kt bin");
+ K2ssT_kt->GetZaxis()->SetTitle("qinv"); K2osT_kt->GetZaxis()->SetTitle("qinv");
+
+ for(int r=0; r<kbVALUES; r++){
+
+ TString *nameT=new TString("Therm_FSI_b");
+
+
+ if(r==0) {*nameT += 2;}
+ if(r==1) {*nameT += 3;}
+ if(r==2) {*nameT += 5;}
+ if(r==3) {*nameT += 7;}
+ if(r==4) {*nameT += 8;}
+ if(r==5) {*nameT += 9;}
+
+ nameT->Append(".root");
+
+
+
+ //
+ TFile *file_T=new TFile(nameT->Data(),"READ");
+ TH2D *Num2_ss = (TH2D*)file_T->Get("K2_ss");
+ TH2D *Den2_ss = (TH2D*)file_T->Get("PlaneWF_ss");
+ TH1D *Num_ss = (TH1D*)Num2_ss->ProjectionY();// kt integrated
+ TH1D *Den_ss = (TH1D*)Den2_ss->ProjectionY();// kt integrated
+ TH2D *Num2_os = (TH2D*)file_T->Get("K2_os");
+ TH2D *Den2_os = (TH2D*)file_T->Get("PlaneWF_os");
+ TH1D *Num_os = (TH1D*)Num2_os->ProjectionY();// kt integrated
+ TH1D *Den_os = (TH1D*)Den2_os->ProjectionY();// kt integrated
+ Num_ss->Divide(Den_ss);
+ Num_os->Divide(Den_os);
+ for(int i=1; i<=binsY; i++){
+ K2ssT->SetBinContent(r+1, i, Num_ss->GetBinContent(i));
+ K2osT->SetBinContent(r+1, i, Num_os->GetBinContent(i));
+ }
+
+ // kt differential
+ for(int ktbin=1; ktbin<=6; ktbin++){
+ TString *name=new TString("pro_");
+ *name += ktbin;
+ *name += 1;
+ TH1D *Num_ss = (TH1D*)Num2_ss->ProjectionY(name->Data(), ktbin*2+3, ktbin*2+4);
+ *name += 2;
+ TH1D *Den_ss = (TH1D*)Den2_ss->ProjectionY(name->Data(), ktbin*2+3, ktbin*2+4);
+ *name += 3;
+ TH1D *Num_os = (TH1D*)Num2_os->ProjectionY(name->Data(), ktbin*2+3, ktbin*2+4);
+ *name += 4;
+ TH1D *Den_os = (TH1D*)Den2_os->ProjectionY(name->Data(), ktbin*2+3, ktbin*2+4);
+ Num_ss->Divide(Den_ss);
+ Num_os->Divide(Den_os);
+ for(int i=1; i<=binsY; i++){
+ K2ssT_kt->SetBinContent(r+1, ktbin, i, Num_ss->GetBinContent(i));
+ K2osT_kt->SetBinContent(r+1, ktbin, i, Num_os->GetBinContent(i));
+ }
+ }
+
+ // 3-particles
+
+ TH3D *Num3_ss=(TH3D*)file_T->Get("K3ss_3D");
+ TH3D *Den3_ss=(TH3D*)file_T->Get("PlaneWF3ss_3D");
+ Num3_ss->Divide(Den3_ss);
+ TString *OutNameSS = new TString("K3ss_");
+ *OutNameSS += r;
+ OutFile->cd();
+ TH3D *K3ss = (TH3D*)Num3_ss->Clone();
+ TH3D *temp3ss = (TH3D*)Num3_ss->Clone();
+ for(int i=1; i<=K3ss->GetNbinsX(); i++){
+ for(int j=1; j<=K3ss->GetNbinsY(); j++){
+ for(int k=1; k<=K3ss->GetNbinsZ(); k++){
+
+ // GRS
+ //double GRS = K2ssT->GetBinContent(r+1, i) * K2ssT->GetBinContent(r+1, j) * K2ssT->GetBinContent(r+1, k);
+ //K3ss->SetBinContent(i,j,k, GRS);
+
+ // Omega0
+ if(temp3ss->GetBinContent(i,j,k) > 1.0) K3ss->SetBinContent(i,j,k, 1.0);
+ else{
+ double mean=0;
+ double terms=0;
+ if(temp3ss->GetBinContent(i,j,k) < 1.0) {mean += temp3ss->GetBinContent(i,j,k); terms++;}
+ if(temp3ss->GetBinContent(i,k,j) < 1.0) {mean += temp3ss->GetBinContent(i,k,j); terms++;}
+ if(temp3ss->GetBinContent(j,i,k) < 1.0) {mean += temp3ss->GetBinContent(j,i,k); terms++;}
+ if(temp3ss->GetBinContent(j,k,i) < 1.0) {mean += temp3ss->GetBinContent(j,k,i); terms++;}
+ if(temp3ss->GetBinContent(k,i,j) < 1.0) {mean += temp3ss->GetBinContent(k,i,j); terms++;}
+ if(temp3ss->GetBinContent(k,j,i) < 1.0) {mean += temp3ss->GetBinContent(k,j,i); terms++;}
+
+ if(terms > 0) {mean /= terms; K3ss->SetBinContent(i,j,k, mean);}
+ else K3ss->SetBinContent(i,j,k, 0);
+ }
+
+ }
+ }
+ }
+ K3ss->Write(OutNameSS->Data());
+ //
+ TH3D *Num3_os=(TH3D*)file_T->Get("K3os_3D");
+ TH3D *Den3_os=(TH3D*)file_T->Get("PlaneWF3os_3D");
+ Num3_os->Divide(Den3_os);
+ //
+ TString *OutNameOS = new TString("K3os_");
+ *OutNameOS += r;
+ OutFile->cd();
+ TH3D *K3os = (TH3D*)Num3_os->Clone();
+ TH3D *temp3os = (TH3D*)Num3_os->Clone();
+ for(int i=1; i<=K3os->GetNbinsX(); i++){
+ for(int j=1; j<=K3os->GetNbinsY(); j++){
+ for(int k=1; k<=K3os->GetNbinsZ(); k++){
+
+ // GRS
+ //double GRS = K2ssT->GetBinContent(r+1, i) * K2osT->GetBinContent(r+1, j) * K2osT->GetBinContent(r+1, k);
+ //K3os->SetBinContent(i,j,k, GRS);
+
+ // Omega0
+ if(temp3os->GetBinContent(i,j,k) > 3.0) K3os->SetBinContent(i,j,k, 1.0);
+ else {
+ double mean=0;
+ double terms=0;
+ if(temp3os->GetBinContent(i,j,k) < 3.0) {mean += temp3os->GetBinContent(i,j,k); terms++;}
+ if(temp3os->GetBinContent(i,k,j) < 3.0) {mean += temp3os->GetBinContent(i,k,j); terms++;}
+
+ if(terms > 0) {mean /= terms; K3os->SetBinContent(i,j,k, mean);}
+ else K3os->SetBinContent(i,j,k, 0);
+ }
+
+ }
+ }
+ }
+
+ K3os->Write(OutNameOS->Data());
+
+
+ file_T->Close();
+
+ }// r loop
+
+ //
+ OutFile->cd();
+ //
+ K2ssT->Write("K2ssT");
+ K2osT->Write("K2osT");
+ K2ssT_kt->Write("K2ssT_kt");
+ K2osT_kt->Write("K2osT_kt");
+ //
+
+ OutFile->Close();
+
+}
--- /dev/null
+#include <math.h>
+#include <time.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <Riostream.h>
+
+#include "TVector2.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"
+
+using namespace std;
+
+void MakeMomResHisto_Merged(){
+
+ int MBmerged = 0;
+ //
+
+ //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_myRun_L0p68R11.root","READ");
+ //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_myRun_L0p68R10.root","READ");
+ //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_lego_L0p68R9andR8.root","READ");
+ //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_myRun_L0p68R7.root","READ");
+ //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_myRun_L0p68R6.root","READ");
+ //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_lego_L0p68R11_PID1p5_and_FB5.root","READ");
+ //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_myRun_L0p68R11_chi3p1.root","READ");
+ TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_myRun_L0p68R11_highKt.root","READ");
+ //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17d_fix_genSignal_Rinv11.root","READ");
+ //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_myRun_L0p68R11_ncls.root","READ");
+ //TFile *infile1 = new TFile("MyOutput.root","READ");
+ //
+ //TFile *infile2 = new TFile("Results/PDC_HIJING_12a17b_12a17e_noPadCut_lego.root","READ");
+ //TFile *infile2 = new TFile("Results/PDC_HIJING_12a17b_newr3_lego.root","READ");
+ //TFile *infile2 = new TFile("Results/PDC_HIJING_12a17b_GRS_lego.root","READ");
+ //TFile *infile2 = new TFile("Results/PDC_HIJING_12a17b_myRun_L0p4R8.root","READ");
+
+ //TDirectoryFile *tdir1 = (TDirectoryFile*)infile1->Get("PWGCF.outputChaoticityAnalysis.root");
+ //TList *list1=(TList*)tdir1->Get("ChaoticityOutput_2");
+ TList *list1=(TList*)infile1->Get("MyList");
+ infile1->Close();
+ //TDirectoryFile *tdir2 = (TDirectoryFile*)infile2->Get("PWGCF.outputChaoticityAnalysis.root");
+ //TList *list2=(TList*)tdir2->Get("ChaoticityOutput");
+ //TList *list2=(TList*)infile2->Get("MyList");
+ //infile2->Close();
+
+ TList *list;
+
+ const int Mbins=2;
+
+ TH2D *num_i_2[2];// 2-particle
+ TH2D *den_i_2[2];// 2-particle
+ TH2D *num_s_2[2];// 2-particle
+ TH2D *den_s_2[2];// 2-particle
+ //
+ TH2D *merged_num_i_2[2];
+ TH2D *merged_den_i_2[2];
+ TH2D *merged_num_s_2[2];
+ TH2D *merged_den_s_2[2];
+ //
+
+
+ TString *names[2][4];
+ // 3d
+ TString *names_terms[2][5][2];// charge-combo (ss or os), terms, ideal/smeared, Mbins
+ TH3D *merged_terms_i[2][5];// charge-combo (ss or os), terms, Mbins
+ TH3D *merged_terms_s[2][5];// charge-combo (ss or os), terms, Mbins
+ TH1D *merged_terms1D_i[2][5];
+ TH1D *merged_terms1D_s[2][5];
+ //
+ TString *names_termsK3[2][2];// charge-combo (ss or os), Sum/En, Mbins
+ TH3D *merged_terms_SumK3[2];// charge-combo (ss or os), Mbins
+ TH3D *merged_terms_EnK3[2];// charge-combo (ss or os), Mbins
+ //
+ TString *names_terms_FVP[5][2][2];// terms, ideal/smeared, Projection 1 or 2, Mbins
+ TH3D *merged_terms_FVP_i[5][2];// terms, Projection 1 or 2, Mbins
+ TH3D *merged_terms_FVP_s[5][2];// terms, Projection 1 or 2, Mbins
+ //
+ TString *names_terms_FVP_K[5][2][2];// terms, Sum/En, Projection 1 or 2, Mbins
+ TH3D *merged_terms_FVP_SumK[5][2];// terms, Projection 1 or 2, Mbins
+ TH3D *merged_terms_FVP_EnK[5][2];// terms, Projection 1 or 2, Mbins
+ //
+ TString *names_terms_FVPTPN[2][2];// ideal/smeared, Projection 1 or 2, Mbins
+ TH3D *merged_terms_FVPTPN_i[2];// Projection 1 or 2, Mbins
+ TH3D *merged_terms_FVPTPN_s[2];// Projection 1 or 2, Mbins
+
+
+ for(int mb=0; mb<Mbins; mb++){// Mbin loop
+ cout<<"Cent bin "<<mb<<endl;
+ for(int ch=0; ch<2; ch++){// + or - loop
+
+ int MixedCbin1=0, MixedCbin2=0, MixedCbin3=0;
+ // 3d histos
+ for(int term=0; term<5; term++){// term loop
+ names_terms[0][term][0]=new TString("PairCut3_Charge1_"); *names_terms[0][term][0] += ch;
+ names_terms[0][term][0]->Append("_Charge2_"); *names_terms[0][term][0] += ch;
+ names_terms[0][term][0]->Append("_Charge3_"); *names_terms[0][term][0] += ch;
+ names_terms[0][term][0]->Append("_SC_0_M_");
+ *names_terms[0][term][0] += mb;
+ names_terms[0][term][0]->Append("_ED_0_");
+ TString *TPNnameBase = new TString(names_terms[0][term][0]->Data());
+ names_terms[0][term][0]->Append("Term_"); *names_terms[0][term][0] += term+1;
+ names_terms[0][term][1]=new TString(names_terms[0][term][0]->Data());
+ if(term==0) {
+ names_termsK3[0][0]=new TString(names_terms[0][term][0]->Data());
+ names_termsK3[0][1]=new TString(names_terms[0][term][0]->Data());
+ names_termsK3[0][0]->Append("_SumK3");
+ names_termsK3[0][1]->Append("_EnK3");
+ }
+ //
+
+ names_terms_FVP[term][0][0] = new TString(names_terms[0][term][0]->Data());
+ names_terms_FVP[term][1][0] = new TString(names_terms[0][term][0]->Data());
+ names_terms_FVP[term][0][0]->Append("_4VectProd1Ideal");
+ names_terms_FVP[term][1][0]->Append("_4VectProd1Smeared");
+ //
+ names_terms_FVP[term][0][1] = new TString(names_terms[0][term][1]->Data());
+ names_terms_FVP[term][1][1] = new TString(names_terms[0][term][1]->Data());
+ names_terms_FVP[term][0][1]->Append("_4VectProd2Ideal");
+ names_terms_FVP[term][1][1]->Append("_4VectProd2Smeared");
+ //
+ names_terms_FVP_K[term][0][0] = new TString(names_terms[0][term][1]->Data());// Prod1 SumK
+ names_terms_FVP_K[term][1][0] = new TString(names_terms[0][term][1]->Data());// Prod1 EnK
+ names_terms_FVP_K[term][0][1] = new TString(names_terms[0][term][1]->Data());// Prod2 SumK
+ names_terms_FVP_K[term][1][1] = new TString(names_terms[0][term][1]->Data());// Prod2 EnK
+ if(term==0) {
+ names_terms_FVP_K[term][0][0]->Append("_4VectProd1SumK3");
+ names_terms_FVP_K[term][1][0]->Append("_4VectProd1EnK3");
+ names_terms_FVP_K[term][0][1]->Append("_4VectProd2SumK3");
+ names_terms_FVP_K[term][1][1]->Append("_4VectProd2EnK3");
+ names_terms_FVPTPN[0][0] = new TString(TPNnameBase->Data());
+ names_terms_FVPTPN[0][1] = new TString(TPNnameBase->Data());
+ names_terms_FVPTPN[1][0] = new TString(TPNnameBase->Data());
+ names_terms_FVPTPN[1][1] = new TString(TPNnameBase->Data());
+ names_terms_FVPTPN[0][0]->Append("TPN_0_4VectProd1Ideal");
+ names_terms_FVPTPN[0][1]->Append("TPN_0_4VectProd2Ideal");
+ names_terms_FVPTPN[1][0]->Append("TPN_0_4VectProd1Smeared");
+ names_terms_FVPTPN[1][1]->Append("TPN_0_4VectProd2Smeared");
+ }
+ if(term>0 && term<4){
+ names_terms_FVP_K[term][0][0]->Append("_4VectProd1SumK2");
+ names_terms_FVP_K[term][1][0]->Append("_4VectProd1EnK2");
+ names_terms_FVP_K[term][0][1]->Append("_4VectProd2SumK2");
+ names_terms_FVP_K[term][1][1]->Append("_4VectProd2EnK2");
+ }
+
+ names_terms[0][term][0]->Append("_Ideal");
+ names_terms[0][term][1]->Append("_Smeared");
+
+ //
+ if(ch==0) {MixedCbin1=0; MixedCbin2=0; MixedCbin3=1;}
+ else {MixedCbin1=0; MixedCbin2=1; MixedCbin3=1;}
+ names_terms[1][term][0]=new TString("PairCut3_Charge1_"); *names_terms[1][term][0] += MixedCbin1;
+ names_terms[1][term][0]->Append("_Charge2_"); *names_terms[1][term][0] += MixedCbin2;
+ names_terms[1][term][0]->Append("_Charge3_"); *names_terms[1][term][0] += MixedCbin3;
+ names_terms[1][term][0]->Append("_SC_0_M_");
+ *names_terms[1][term][0] += mb;
+ names_terms[1][term][0]->Append("_ED_0_Term_"); *names_terms[1][term][0] += term+1;
+ names_terms[1][term][1]=new TString(names_terms[1][term][0]->Data());
+ if(term==0) {
+ names_termsK3[1][0]=new TString(names_terms[1][term][0]->Data());
+ names_termsK3[1][1]=new TString(names_terms[1][term][0]->Data());
+ names_termsK3[1][0]->Append("_SumK3");
+ names_termsK3[1][1]->Append("_EnK3");
+ }
+ names_terms[1][term][0]->Append("_Ideal");
+ names_terms[1][term][1]->Append("_Smeared");
+ }// term loop
+
+
+ TString *base1 = new TString("Explicit2_Charge1_"); *base1 += ch;
+ base1->Append("_Charge2_"); *base1 += ch; base1->Append("_SC_0_M_");
+ names[0][0]=new TString(base1->Data());
+ *names[0][0] += mb;
+ names[0][0]->Append("_ED_0_Term_1_Ideal");
+ names[0][1]=new TString(base1->Data());
+ *names[0][1] += mb;
+ names[0][1]->Append("_ED_0_Term_2_Ideal");
+ names[0][2]=new TString(base1->Data());
+ *names[0][2] += mb;
+ names[0][2]->Append("_ED_0_Term_1_Smeared");
+ names[0][3]=new TString(base1->Data());
+ *names[0][3] += mb;
+ names[0][3]->Append("_ED_0_Term_2_Smeared");
+ // mixed charge 2-particle case is double counted here. It's ok (only the stat errors are wrong...not used anyway).
+ TString *base2 = new TString("Explicit2_Charge1_0_Charge2_1"); base2->Append("_SC_0_M_");
+ names[1][0]=new TString(base2->Data());
+ *names[1][0] += mb;
+ names[1][0]->Append("_ED_0_Term_1_Ideal");
+ names[1][1]=new TString(base2->Data());
+ *names[1][1] += mb;
+ names[1][1]->Append("_ED_0_Term_2_Ideal");
+ names[1][2]=new TString(base2->Data());
+ *names[1][2] += mb;
+ names[1][2]->Append("_ED_0_Term_1_Smeared");
+ names[1][3]=new TString(base2->Data());
+ *names[1][3] += mb;
+ names[1][3]->Append("_ED_0_Term_2_Smeared");
+ //
+ if(mb<2) list=(TList*)list1->Clone();
+ //else list=(TList*)list2->Clone();
+
+
+ for(int ChProd=0; ChProd<2; ChProd++){// charge product loop
+ //
+ num_i_2[ChProd]=(TH2D*)list->FindObject(names[ChProd][0]->Data());
+ den_i_2[ChProd]=(TH2D*)list->FindObject(names[ChProd][1]->Data());
+ num_s_2[ChProd]=(TH2D*)list->FindObject(names[ChProd][2]->Data());
+ den_s_2[ChProd]=(TH2D*)list->FindObject(names[ChProd][3]->Data());
+ //
+
+ if(ch==0 && mb==0){
+ for(int term=0; term<5; term++){
+ merged_terms_i[ChProd][term] = (TH3D*)(list->FindObject(names_terms[ChProd][term][0]->Data()))->Clone();
+ merged_terms_s[ChProd][term] = (TH3D*)(list->FindObject(names_terms[ChProd][term][1]->Data()))->Clone();
+
+ if(term==0){
+ merged_terms_SumK3[ChProd] = (TH3D*)(list->FindObject(names_termsK3[ChProd][0]->Data()))->Clone();
+ merged_terms_EnK3[ChProd] = (TH3D*)(list->FindObject(names_termsK3[ChProd][1]->Data()))->Clone();
+ }
+ if(ChProd==0) {
+ //merged_terms_FVP_i[term][0] = (TH3D*)(list->FindObject(names_terms_FVP[term][0][0]->Data()))->Clone();// ideal, Prod 1
+ //merged_terms_FVP_s[term][0] = (TH3D*)(list->FindObject(names_terms_FVP[term][1][0]->Data()))->Clone();// smeared, Prod 1
+ //merged_terms_FVP_i[term][1] = (TH3D*)(list->FindObject(names_terms_FVP[term][0][1]->Data()))->Clone();// ideal, Prod 2
+ //merged_terms_FVP_s[term][1] = (TH3D*)(list->FindObject(names_terms_FVP[term][1][1]->Data()))->Clone();// smeared, Prod 2
+ if(term !=4){
+ //merged_terms_FVP_SumK[term][0] = (TH3D*)(list->FindObject(names_terms_FVP_K[term][0][0]->Data()))->Clone();// Prod 1 Sum K
+ //merged_terms_FVP_EnK[term][0] = (TH3D*)(list->FindObject(names_terms_FVP_K[term][1][0]->Data()))->Clone();// Prod 1 En K
+ //merged_terms_FVP_SumK[term][1] = (TH3D*)(list->FindObject(names_terms_FVP_K[term][0][1]->Data()))->Clone();// Prod 2 Sum K
+ //merged_terms_FVP_EnK[term][1] = (TH3D*)(list->FindObject(names_terms_FVP_K[term][1][1]->Data()))->Clone();// Prod 2 En K
+ }
+ if(term==0){
+ //merged_terms_FVPTPN_i[0] = (TH3D*)(list->FindObject(names_terms_FVPTPN[0][0]->Data()))->Clone();// ideal, Prod 1
+ //merged_terms_FVPTPN_s[0] = (TH3D*)(list->FindObject(names_terms_FVPTPN[1][0]->Data()))->Clone();// smeared, Prod 1
+ //merged_terms_FVPTPN_i[1] = (TH3D*)(list->FindObject(names_terms_FVPTPN[0][1]->Data()))->Clone();// ideal, Prod 2
+ //merged_terms_FVPTPN_s[1] = (TH3D*)(list->FindObject(names_terms_FVPTPN[1][1]->Data()))->Clone();// smeared, Prod 2
+ }
+ }
+ }
+ }else{// ch==1 || mb!=0
+ for(int term=0; term<5; term++){
+ if(ChProd==0) {
+ merged_terms_i[0][term]->Add((TH3D*)(list->FindObject(names_terms[0][term][0]->Data())));
+ merged_terms_s[0][term]->Add((TH3D*)(list->FindObject(names_terms[0][term][1]->Data())));
+ if(term==0) {
+ merged_terms_SumK3[0]->Add((TH3D*)(list->FindObject(names_termsK3[0][0]->Data())));
+ merged_terms_EnK3[0]->Add((TH3D*)(list->FindObject(names_termsK3[0][1]->Data())));
+ /*merged_terms_FVPTPN_i[0]->Add((TH3D*)(list->FindObject(names_terms_FVPTPN[0][0]->Data())));
+ merged_terms_FVPTPN_i[1]->Add((TH3D*)(list->FindObject(names_terms_FVPTPN[0][1]->Data())));
+ merged_terms_FVPTPN_s[0]->Add((TH3D*)(list->FindObject(names_terms_FVPTPN[1][0]->Data())));
+ merged_terms_FVPTPN_s[1]->Add((TH3D*)(list->FindObject(names_terms_FVPTPN[1][1]->Data())));*/
+ }
+ /* merged_terms_FVP_i[term][0]->Add((TH3D*)(list->FindObject(names_terms_FVP[term][0][0]->Data())));
+ merged_terms_FVP_s[term][0]->Add((TH3D*)(list->FindObject(names_terms_FVP[term][1][0]->Data())));
+ merged_terms_FVP_i[term][1]->Add((TH3D*)(list->FindObject(names_terms_FVP[term][0][1]->Data())));
+ merged_terms_FVP_s[term][1]->Add((TH3D*)(list->FindObject(names_terms_FVP[term][1][1]->Data())));*/
+
+ if(term<4){
+ /*merged_terms_FVP_SumK[term][0]->Add((TH3D*)(list->FindObject(names_terms_FVP_K[term][0][0]->Data())));
+ merged_terms_FVP_EnK[term][0]->Add((TH3D*)(list->FindObject(names_terms_FVP_K[term][1][0]->Data())));
+ merged_terms_FVP_SumK[term][1]->Add((TH3D*)(list->FindObject(names_terms_FVP_K[term][0][1]->Data())));
+ merged_terms_FVP_EnK[term][1]->Add((TH3D*)(list->FindObject(names_terms_FVP_K[term][1][1]->Data())));*/
+ }
+ }else {// different placement of same-charge pair
+ int termOther=term;
+ if(ch==1){
+ if(term==0) termOther=0;
+ else if(term==1) termOther=3;
+ else if(term==2) termOther=2;
+ else if(term==3) termOther=1;
+ else termOther=4;
+ }
+ TH3D *temphisto_i=((TH3D*)(list->FindObject(names_terms[1][termOther][0]->Data())));
+ TH3D *temphisto_s=((TH3D*)(list->FindObject(names_terms[1][termOther][1]->Data())));
+ TH3D *temphisto_SumK3=((TH3D*)(list->FindObject(names_termsK3[1][0]->Data())));
+ TH3D *temphisto_EnK3=((TH3D*)(list->FindObject(names_termsK3[1][1]->Data())));
+ for(int bin1=1; bin1<merged_terms_i[1][term]->GetNbinsX(); bin1++){
+ for(int bin2=1; bin2<merged_terms_i[1][term]->GetNbinsY(); bin2++){
+ for(int bin3=1; bin3<merged_terms_i[1][term]->GetNbinsZ(); bin3++){
+
+ if(ch==0){
+ merged_terms_i[1][term]->SetBinContent(bin1,bin2,bin3, merged_terms_i[1][term]->GetBinContent(bin1,bin2,bin3)+temphisto_i->GetBinContent(bin1,bin2,bin3));
+ merged_terms_s[1][term]->SetBinContent(bin1,bin2,bin3, merged_terms_s[1][term]->GetBinContent(bin1,bin2,bin3)+temphisto_s->GetBinContent(bin1,bin2,bin3));
+ }else {
+ merged_terms_i[1][term]->SetBinContent(bin1,bin2,bin3, merged_terms_i[1][term]->GetBinContent(bin1,bin2,bin3)+temphisto_i->GetBinContent(bin3,bin2,bin1));
+ merged_terms_s[1][term]->SetBinContent(bin1,bin2,bin3, merged_terms_s[1][term]->GetBinContent(bin1,bin2,bin3)+temphisto_s->GetBinContent(bin3,bin2,bin1));
+ }
+
+ //
+ if(term==0){
+ //merged_terms_SumK3[1]->SetBinContent(bin1,bin2,bin3, merged_terms_SumK3[1]->GetBinContent(bin1,bin2,bin3)+temphisto_SumK3->GetBinContent(bin3,bin2,bin1));
+ //merged_terms_EnK3[1]->SetBinContent(bin1,bin2,bin3, merged_terms_EnK3[1]->GetBinContent(bin1,bin2,bin3)+temphisto_EnK3->GetBinContent(bin3,bin2,bin1));
+ }
+ }
+ }
+ }
+ }
+ }
+ }// else for 3-particle part
+
+ // 2-particles
+ if(ch==0 && mb==0) {
+ merged_num_i_2[ChProd] = (TH2D*)num_i_2[ChProd]->Clone();
+ merged_den_i_2[ChProd] = (TH2D*)den_i_2[ChProd]->Clone();
+ merged_num_s_2[ChProd] = (TH2D*)num_s_2[ChProd]->Clone();
+ merged_den_s_2[ChProd] = (TH2D*)den_s_2[ChProd]->Clone();
+ }else {
+ merged_num_i_2[ChProd]->Add(num_i_2[ChProd]);
+ merged_den_i_2[ChProd]->Add(den_i_2[ChProd]);
+ merged_num_s_2[ChProd]->Add(num_s_2[ChProd]);
+ merged_den_s_2[ChProd]->Add(den_s_2[ChProd]);
+
+ }
+ }// ch loop
+
+
+ }// mb loop
+
+ }// ChProd loop
+
+
+ cout<<"Start Writing to Output File"<<endl;
+ TFile *outfile = new TFile("MomResFile_temp.root","RECREATE");
+ TFile *outfileOffline = new TFile("MomResFile_Offline_temp.root","RECREATE");
+ int NbinsX = merged_num_i_2[0]->GetNbinsX();
+ int NbinsY = merged_num_i_2[0]->GetNbinsY();
+ float Xlow = merged_num_i_2[0]->GetXaxis()->GetBinLowEdge(1);
+ float Xhigh = merged_num_i_2[0]->GetXaxis()->GetBinUpEdge(NbinsX);
+ float Ylow = merged_num_i_2[0]->GetYaxis()->GetBinLowEdge(1);
+ float Yhigh = merged_num_i_2[0]->GetYaxis()->GetBinUpEdge(NbinsY);
+ TH2D *MomResHisto_pp = new TH2D("MomResHisto_pp","",NbinsX,Xlow,Xhigh, NbinsY,Ylow,Yhigh);
+ TH2D *MomResHisto_mp = new TH2D("MomResHisto_mp","",NbinsX,Xlow,Xhigh, NbinsY,Ylow,Yhigh);
+ //
+ TH2D *MomResHisto_pp_term1 = new TH2D("MomResHisto_pp_term1","",NbinsX,Xlow,Xhigh, NbinsY,Ylow,Yhigh);
+ TH2D *MomResHisto_pp_term2 = new TH2D("MomResHisto_pp_term2","",NbinsX,Xlow,Xhigh, NbinsY,Ylow,Yhigh);
+ TH2D *MomResHisto_mp_term1 = new TH2D("MomResHisto_mp_term1","",NbinsX,Xlow,Xhigh, NbinsY,Ylow,Yhigh);
+ TH2D *MomResHisto_mp_term2 = new TH2D("MomResHisto_mp_term2","",NbinsX,Xlow,Xhigh, NbinsY,Ylow,Yhigh);
+ //
+ // 3d part
+ int NbinsX_3d = merged_terms_i[0][0]->GetNbinsX();
+ float Xlow_3d = merged_terms_i[0][0]->GetXaxis()->GetBinLowEdge(1);
+ float Xhigh_3d = merged_terms_i[0][0]->GetXaxis()->GetBinUpEdge(NbinsX_3d);
+ TH3D *MomResHisto_3d_SC[5];
+ TH3D *MomResHisto_3d_MC[5];
+ TH1D *MomResHisto_1d_SC[5];
+ TH1D *MomResHisto_1d_MC[5];
+ TH3D *MomResHisto_3d_SC_K3;
+ TH3D *MomResHisto_3d_MC_K3;
+ // 4vectProd
+ /*const int NEdges_FVP_temp = merged_terms_FVP_i[0][0]->GetNbinsX() + 1;
+ const int NEdges_FVP = NEdges_FVP_temp;
+ double Edges[NEdges_FVP]={0};
+ for(int edg=0; edg<NEdges_FVP; edg++){
+ Edges[edg] = merged_terms_FVP_i[0][0]->GetXaxis()->GetBinLowEdge(edg+1);
+ //cout<<Edges[edg]<<endl;
+ }
+ TH3D *MomResHisto_3d_FVP1[5];
+ TH3D *MomResHisto_3d_FVP2[5];
+ TH3D *MomResHisto_3d_FVP1K[5];
+ TH3D *MomResHisto_3d_FVP2K[5];
+ TH3D *MomResHisto_TPN_FVP1;
+ TH3D *MomResHisto_TPN_FVP2;*/
+
+
+ for(int term=0; term<5; term++){
+ TString *nameSC = new TString("MomResHisto3_SC_term");
+ TString *nameMC = new TString("MomResHisto3_MC_term");
+ TString *nameSC1D = new TString("MomResHisto1D_SC_term");
+ TString *nameMC1D = new TString("MomResHisto1D_MC_term");
+ TString *nameFVP1 = new TString("MomResHisto3_FVP1_term");
+ TString *nameFVP2 = new TString("MomResHisto3_FVP2_term");
+ TString *nameFVP1K = new TString("MomResHisto3_FVP1K_term");
+ TString *nameFVP2K = new TString("MomResHisto3_FVP2K_term");
+ *nameSC += term+1;
+ *nameMC += term+1;
+ *nameSC1D += term+1;
+ *nameMC1D += term+1;
+ *nameFVP1 += term+1;
+ *nameFVP2 += term+1;
+ *nameFVP1K += term+1;
+ *nameFVP2K += term+1;
+ nameSC->Append("_M");
+ nameMC->Append("_M");
+ nameSC1D->Append("_M");
+ nameMC1D->Append("_M");
+ nameFVP1->Append("_M");
+ nameFVP2->Append("_M");
+ nameFVP1K->Append("_M");
+ nameFVP2K->Append("_M");
+ *nameSC += MBmerged;
+ *nameMC += MBmerged;
+ *nameSC1D += MBmerged;
+ *nameMC1D += MBmerged;
+ *nameFVP1 += MBmerged;
+ *nameFVP2 += MBmerged;
+ *nameFVP1K += MBmerged;
+ *nameFVP2K += MBmerged;
+ MomResHisto_3d_SC[term] = new TH3D(nameSC->Data(),"",NbinsX_3d,Xlow_3d,Xhigh_3d, NbinsX_3d,Xlow_3d,Xhigh_3d, NbinsX_3d,Xlow_3d,Xhigh_3d);
+ MomResHisto_3d_MC[term] = new TH3D(nameMC->Data(),"",NbinsX_3d,Xlow_3d,Xhigh_3d, NbinsX_3d,Xlow_3d,Xhigh_3d, NbinsX_3d,Xlow_3d,Xhigh_3d);
+ MomResHisto_1d_SC[term] = new TH1D(nameSC1D->Data(),"",20,0,0.2);
+ MomResHisto_1d_MC[term] = new TH1D(nameMC1D->Data(),"",20,0,0.2);
+ /*MomResHisto_3d_FVP1[term] = new TH3D(nameFVP1->Data(),"",NEdges_FVP-1,Edges, NEdges_FVP-1,Edges, NEdges_FVP-1,Edges);
+ MomResHisto_3d_FVP2[term] = new TH3D(nameFVP2->Data(),"",NEdges_FVP-1,Edges, NEdges_FVP-1,Edges, NEdges_FVP-1,Edges);
+ MomResHisto_3d_FVP1K[term] = new TH3D(nameFVP1K->Data(),"",NEdges_FVP-1,Edges, NEdges_FVP-1,Edges, NEdges_FVP-1,Edges);
+ MomResHisto_3d_FVP2K[term] = new TH3D(nameFVP2K->Data(),"",NEdges_FVP-1,Edges, NEdges_FVP-1,Edges, NEdges_FVP-1,Edges);*/
+ }
+ //TString *nameTPN_FVP1 = new TString("MomResHisto3_TPN_FVP1_M");
+ //TString *nameTPN_FVP2 = new TString("MomResHisto3_TPN_FVP2_M");
+ TString *nameSC_K3 = new TString("AvgK3_SC_M");
+ TString *nameMC_K3 = new TString("AvgK3_MC_M");
+ *nameSC_K3 += MBmerged;
+ *nameMC_K3 += MBmerged;
+ //*nameTPN_FVP1 += MBmerged;
+ //*nameTPN_FVP2 += MBmerged;
+ MomResHisto_3d_SC_K3 = new TH3D(nameSC_K3->Data(),"",NbinsX_3d,Xlow_3d,Xhigh_3d, NbinsX_3d,Xlow_3d,Xhigh_3d, NbinsX_3d,Xlow_3d,Xhigh_3d);
+ MomResHisto_3d_MC_K3 = new TH3D(nameMC_K3->Data(),"",NbinsX_3d,Xlow_3d,Xhigh_3d, NbinsX_3d,Xlow_3d,Xhigh_3d, NbinsX_3d,Xlow_3d,Xhigh_3d);
+ //MomResHisto_TPN_FVP1 = new TH3D(nameTPN_FVP1->Data(),"",NEdges_FVP-1,Edges, NEdges_FVP-1,Edges, NEdges_FVP-1,Edges);
+ //MomResHisto_TPN_FVP2 = new TH3D(nameTPN_FVP2->Data(),"",NEdges_FVP-1,Edges, NEdges_FVP-1,Edges, NEdges_FVP-1,Edges);
+
+
+
+ for(int term=0; term<5; term++){
+
+ // 1D Q3 projection
+ double Sum_i_SC[20]={0};
+ double Sum_i_MC[20]={0};
+ double Sum_s_SC[20]={0};
+ double Sum_s_MC[20]={0};
+
+ for(int i=1; i<NbinsX_3d; i++){
+ for(int j=1; j<NbinsX_3d; j++){
+ for(int l=1; l<NbinsX_3d; l++){
+ //if((i+j)==2 || (i+l)==2 || (j+l)==2) continue;
+ double q3 = sqrt(pow(merged_terms_i[0][term]->GetXaxis()->GetBinCenter(i+1),2)+pow(merged_terms_i[0][term]->GetYaxis()->GetBinCenter(j+1),2)+pow(merged_terms_i[0][term]->GetZaxis()->GetBinCenter(l+1),2));
+ int q3bin = MomResHisto_1d_SC[term]->GetXaxis()->FindBin(q3);
+ if(q3bin>20) continue;
+ Sum_i_SC[q3bin-1] += merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1);
+ Sum_i_MC[q3bin-1] += merged_terms_i[1][term]->GetBinContent(i+1,j+1,l+1);
+ Sum_s_SC[q3bin-1] += merged_terms_s[0][term]->GetBinContent(i+1,j+1,l+1);
+ Sum_s_MC[q3bin-1] += merged_terms_s[1][term]->GetBinContent(i+1,j+1,l+1);
+ }
+ }
+ }
+ for(int i=0; i<20; i++){
+
+ if(Sum_s_SC[i]>0) MomResHisto_1d_SC[term]->SetBinContent(i+1, Sum_i_SC[i]/Sum_s_SC[i]);
+ if(Sum_s_MC[i]>0) MomResHisto_1d_MC[term]->SetBinContent(i+1, Sum_i_MC[i]/Sum_s_MC[i]);
+ }
+
+ // full 3d method
+ merged_terms_i[0][term]->Divide(merged_terms_s[0][term]);
+ merged_terms_i[1][term]->Divide(merged_terms_s[1][term]);
+ //merged_terms_FVP_i[term][0]->Divide(merged_terms_FVP_s[term][0]);// Prod 1
+ //merged_terms_FVP_i[term][1]->Divide(merged_terms_FVP_s[term][1]);// Prod 2
+ if(term < 4){
+ //merged_terms_FVP_SumK[term][0]->Divide(merged_terms_FVP_EnK[term][0]);// Prod 1
+ //merged_terms_FVP_SumK[term][1]->Divide(merged_terms_FVP_EnK[term][1]);// Prod 2
+ }
+ }
+ //merged_terms_FVPTPN_i[0]->Divide(merged_terms_FVPTPN_s[0]);// Prod 1
+ //merged_terms_FVPTPN_i[1]->Divide(merged_terms_FVPTPN_s[1]);// Prod 2
+ merged_terms_SumK3[0]->Divide(merged_terms_EnK3[0]);// Avg SC K3 for Q12,Q13,Q23
+ merged_terms_SumK3[1]->Divide(merged_terms_EnK3[1]);// Avg MC K3 for Q12,Q13,Q23
+
+
+
+ for(int term=0; term<5; term++){
+ for(int i=0; i<NbinsX_3d; i++){
+ for(int j=0; j<NbinsX_3d; j++){
+ for(int l=0; l<NbinsX_3d; l++){
+
+ double MomResCorr=0; double terms=0;
+ if(term==0 || term==4){
+ if(merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1); terms++;}
+ if(merged_terms_i[0][term]->GetBinContent(i+1,l+1,j+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(i+1,l+1,j+1); terms++;}
+ if(merged_terms_i[0][term]->GetBinContent(j+1,i+1,l+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(j+1,i+1,l+1); terms++;}
+ if(merged_terms_i[0][term]->GetBinContent(j+1,l+1,i+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(j+1,l+1,i+1); terms++;}
+ if(merged_terms_i[0][term]->GetBinContent(l+1,i+1,j+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(l+1,i+1,j+1); terms++;}
+ if(merged_terms_i[0][term]->GetBinContent(l+1,j+1,i+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(l+1,j+1,i+1); terms++;}
+ }else if(term==1){
+ if(merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1); terms++;}
+ if(merged_terms_i[0][term]->GetBinContent(i+1,l+1,j+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(i+1,l+1,j+1); terms++;}
+ }else if(term==2){
+ if(merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1); terms++;}
+ if(merged_terms_i[0][term]->GetBinContent(l+1,j+1,i+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(l+1,j+1,i+1); terms++;}
+ }else {
+ if(merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1); terms++;}
+ if(merged_terms_i[0][term]->GetBinContent(j+1,i+1,l+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(j+1,i+1,l+1); terms++;}
+ }
+
+ if(terms > 0) MomResCorr /= terms;
+ MomResHisto_3d_SC[term]->SetBinContent(i+1,j+1,l+1, MomResCorr);
+
+ // Mixed-Charge
+ MomResCorr=0; terms=0;
+ if(term==0 || term==1 || term==4){
+ if(merged_terms_i[1][term]->GetBinContent(i+1,j+1,l+1)>0) {MomResCorr += merged_terms_i[1][term]->GetBinContent(i+1,j+1,l+1); terms++;}
+ if(merged_terms_i[1][term]->GetBinContent(i+1,l+1,j+1)>0) {MomResCorr += merged_terms_i[1][term]->GetBinContent(i+1,l+1,j+1); terms++;}
+ }else{
+ if(merged_terms_i[1][term]->GetBinContent(i+1,j+1,l+1)>0) {MomResCorr += merged_terms_i[1][term]->GetBinContent(i+1,j+1,l+1); terms++;}
+ }
+
+ if(terms > 0) MomResCorr /= terms;
+ MomResHisto_3d_MC[term]->SetBinContent(i+1,j+1,l+1, MomResCorr);
+
+ if(term==0){
+ MomResHisto_3d_SC_K3->SetBinContent(i+1,j+1,l+1, merged_terms_SumK3[0]->GetBinContent(i+1,j+1,l+1));
+ MomResHisto_3d_MC_K3->SetBinContent(i+1,j+1,l+1, merged_terms_SumK3[1]->GetBinContent(i+1,j+1,l+1));
+ }
+ // condition edge effects
+ if(MomResHisto_3d_SC[term]->GetBinContent(i+1,j+1,l+1) > 10) MomResHisto_3d_SC[term]->SetBinContent(i+1,j+1,l+1, 1.0);
+ if(MomResHisto_3d_MC[term]->GetBinContent(i+1,j+1,l+1) > 10) MomResHisto_3d_MC[term]->SetBinContent(i+1,j+1,l+1, 1.0);
+ }
+ }
+ }
+ }// term
+
+ // FVP
+ /*for(int mb=0; mb<1; mb++){
+ for(int term=0; term<5; term++){
+ for(int i=0; i<NEdges_FVP-1; i++){
+ for(int j=0; j<NEdges_FVP-1; j++){
+ for(int l=0; l<NEdges_FVP-1; l++){
+ MomResHisto_3d_FVP1[term]->SetBinContent(i+1,j+1,l+1, merged_terms_FVP_i[term][0]->GetBinContent(i+1,j+1,l+1));
+ MomResHisto_3d_FVP2[term]->SetBinContent(i+1,j+1,l+1, merged_terms_FVP_i[term][1]->GetBinContent(i+1,j+1,l+1));
+ if(term < 4){
+ MomResHisto_3d_FVP1K[term]->SetBinContent(i+1,j+1,l+1, merged_terms_FVP_SumK[term][0]->GetBinContent(i+1,j+1,l+1));
+ MomResHisto_3d_FVP2K[term]->SetBinContent(i+1,j+1,l+1, merged_terms_FVP_SumK[term][1]->GetBinContent(i+1,j+1,l+1));
+ }
+ if(term==0){
+ MomResHisto_TPN_FVP1->SetBinContent(i+1,j+1,l+1, merged_terms_FVPTPN_i[0]->GetBinContent(i+1,j+1,l+1));
+ MomResHisto_TPN_FVP2->SetBinContent(i+1,j+1,l+1, merged_terms_FVPTPN_i[1]->GetBinContent(i+1,j+1,l+1));
+ }
+ // condition edge effects
+ if(MomResHisto_3d_FVP1[term]->GetBinContent(i+1,j+1,l+1) > 10) MomResHisto_3d_FVP1[term]->SetBinContent(i+1,j+1,l+1, 1.0);
+ if(MomResHisto_3d_FVP2[term]->GetBinContent(i+1,j+1,l+1) > 10) MomResHisto_3d_FVP2[term]->SetBinContent(i+1,j+1,l+1, 1.0);
+ }
+ }
+ }
+ }
+ }*/
+
+ // 2-particles
+ for(int ChProd=0; ChProd<2; ChProd++){
+ for(int i=0; i<NbinsX; i++){
+ for(int j=0; j<NbinsY; j++){
+
+ double weight = 1.0, weight_term1=1.0, weight_term2=1.0;
+ double weight_e = 0, weight_term1_e=0, weight_term2_e=0;
+ bool goodbin=kTRUE;
+ if(merged_num_i_2[ChProd]->GetBinContent(i+1,j+1) == 0) goodbin=kFALSE;
+ if(merged_den_i_2[ChProd]->GetBinContent(i+1,j+1) == 0) goodbin=kFALSE;
+ if(merged_num_s_2[ChProd]->GetBinContent(i+1,j+1) == 0) goodbin=kFALSE;
+ if(merged_den_s_2[ChProd]->GetBinContent(i+1,j+1) == 0) goodbin=kFALSE;
+
+ if(goodbin==kTRUE){
+ double NUM_i=merged_num_i_2[ChProd]->GetBinContent(i+1,j+1);
+ double NUM_s=merged_num_s_2[ChProd]->GetBinContent(i+1,j+1);
+ double DEN_i=merged_den_i_2[ChProd]->GetBinContent(i+1,j+1);
+ double DEN_s=merged_den_s_2[ChProd]->GetBinContent(i+1,j+1);
+ weight = (NUM_i/DEN_i)/(NUM_s/DEN_s);
+ weight_e = pow((sqrt(NUM_s)*NUM_i/DEN_i)/(pow(NUM_s,2)/DEN_s),2);// approximate. does not take lambda into account!!
+ weight_e += pow((sqrt(DEN_s)*NUM_i/DEN_i)/(NUM_s),2);
+ weight_e = sqrt(weight_e);
+ //
+ weight_term1 = NUM_i/NUM_s;
+ weight_term1_e = pow(sqrt(NUM_i)/NUM_s,2);
+ weight_term1_e += pow(sqrt(NUM_s)*NUM_i/pow(NUM_s,2),2);
+ weight_term1_e = sqrt(weight_term1_e);
+ weight_term2 = DEN_i/DEN_s;
+ weight_term2_e = pow(sqrt(DEN_i)/DEN_s,2);
+ weight_term2_e += pow(sqrt(DEN_s)*DEN_i/pow(DEN_s,2),2);
+ weight_term2_e = sqrt(weight_term2_e);
+ }
+
+ if(ChProd==0) {
+ MomResHisto_pp->SetBinContent(i+1,j+1, weight);
+ MomResHisto_pp->SetBinError(i+1,j+1, weight_e);
+ MomResHisto_pp_term1->SetBinContent(i+1,j+1, weight_term1);
+ MomResHisto_pp_term2->SetBinContent(i+1,j+1, weight_term2);
+ MomResHisto_pp_term1->SetBinError(i+1,j+1, weight_term1_e);
+ MomResHisto_pp_term2->SetBinError(i+1,j+1, weight_term2_e);
+ }else {
+ MomResHisto_mp->SetBinContent(i+1,j+1, weight);
+ MomResHisto_mp->SetBinError(i+1,j+1, weight_e);
+ MomResHisto_mp_term1->SetBinContent(i+1,j+1, weight_term1);
+ MomResHisto_mp_term2->SetBinContent(i+1,j+1, weight_term2);
+ MomResHisto_mp_term1->SetBinError(i+1,j+1, weight_term1_e);
+ MomResHisto_mp_term2->SetBinError(i+1,j+1, weight_term2_e);
+ }
+ }
+ }
+ }
+
+
+
+ outfile->cd();
+ MomResHisto_pp->Write("MomResHisto_pp");
+ MomResHisto_mp->Write("MomResHisto_mp");
+ MomResHisto_pp_term1->Write("MomResHisto_pp_term1");
+ MomResHisto_pp_term2->Write("MomResHisto_pp_term2");
+ MomResHisto_mp_term1->Write("MomResHisto_mp_term1");
+ MomResHisto_mp_term2->Write("MomResHisto_mp_term2");
+ outfile->Close();
+
+ outfileOffline->cd();
+ if(MBmerged==0){
+ MomResHisto_pp->Write("MomResHisto_pp");
+ MomResHisto_mp->Write("MomResHisto_mp");
+ MomResHisto_pp_term1->Write("MomResHisto_pp_term1");
+ MomResHisto_pp_term2->Write("MomResHisto_pp_term2");
+ MomResHisto_mp_term1->Write("MomResHisto_mp_term1");
+ MomResHisto_mp_term2->Write("MomResHisto_mp_term2");
+ }
+ //
+
+ for(int term=0; term<5; term++){
+ MomResHisto_3d_SC[term]->Write();
+ MomResHisto_3d_MC[term]->Write();
+ MomResHisto_1d_SC[term]->Write();
+ MomResHisto_1d_MC[term]->Write();
+ //MomResHisto_3d_FVP1[term]->Write();
+ //MomResHisto_3d_FVP2[term]->Write();
+ if(term < 4){
+ //MomResHisto_3d_FVP1K[term]->Write();
+ //MomResHisto_3d_FVP2K[term]->Write();
+ }
+ }
+ MomResHisto_3d_SC_K3->Write();
+ MomResHisto_3d_MC_K3->Write();
+ //MomResHisto_TPN_FVP1->Write();
+ //MomResHisto_TPN_FVP2->Write();
+
+
+ outfileOffline->Close();
+
+}
--- /dev/null
+#!/bin/bash
+
+# Save, MCcase, IncludeEWfromTherm, SC, G, EW, GRS, EDBin, CH, Mbin, Ktbin
+# without G
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, kTRUE, 0, +1, '$1', 1)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, kTRUE, 0, +1, '$1', 2)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, kTRUE, 0, +1, '$1', 3)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, kTRUE, 0, +1, '$1', 4)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, kTRUE, 0, +1, '$1', 5)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kFALSE, kFALSE, kTRUE, 0, +1, '$1', 6)'
+# with G
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kTRUE, kFALSE, kTRUE, 0, +1, '$1', 1)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kTRUE, kFALSE, kTRUE, 0, +1, '$1', 2)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kTRUE, kFALSE, kTRUE, 0, +1, '$1', 3)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kTRUE, kFALSE, kTRUE, 0, +1, '$1', 4)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kTRUE, kFALSE, kTRUE, 0, +1, '$1', 5)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kTRUE, kFALSE, kTRUE, 0, +1, '$1', 6)'
+#
+# Include EW from Therminator fits
+# without G
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kTRUE, kTRUE, kFALSE, kFALSE, kTRUE, 0, +1, '$1', 1)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kTRUE, kTRUE, kFALSE, kFALSE, kTRUE, 0, +1, '$1', 2)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kTRUE, kTRUE, kFALSE, kFALSE, kTRUE, 0, +1, '$1', 3)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kTRUE, kTRUE, kFALSE, kFALSE, kTRUE, 0, +1, '$1', 4)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kTRUE, kTRUE, kFALSE, kFALSE, kTRUE, 0, +1, '$1', 5)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kTRUE, kTRUE, kFALSE, kFALSE, kTRUE, 0, +1, '$1', 6)'
+# with G
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kTRUE, kTRUE, kTRUE, kFALSE, kTRUE, 0, +1, '$1', 1)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kTRUE, kTRUE, kTRUE, kFALSE, kTRUE, 0, +1, '$1', 2)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kTRUE, kTRUE, kTRUE, kFALSE, kTRUE, 0, +1, '$1', 3)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kTRUE, kTRUE, kTRUE, kFALSE, kTRUE, 0, +1, '$1', 4)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kTRUE, kTRUE, kTRUE, kFALSE, kTRUE, 0, +1, '$1', 5)'
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kTRUE, kTRUE, kTRUE, kFALSE, kTRUE, 0, +1, '$1', 6)'
+#
+# 3-particle
+# EDbin 1
+# GRS
+# SC EW, +
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kFALSE, kTRUE, kTRUE, 0, +1, '$1', 10)'
+# SC EW, -
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kFALSE, kTRUE, kTRUE, 0, -1, '$1', 10)'
+# MC EW, +
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kTRUE, 0, +1, '$1', 10)'
+# MC EW, -
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kTRUE, 0, -1, '$1', 10)'
+# Omega0
+# SC EW, +
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kFALSE, kTRUE, kFALSE, 0, +1, '$1', 10)'
+# SC EW, -
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kFALSE, kTRUE, kFALSE, 0, -1, '$1', 10)'
+# MC EW, +
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kFALSE, 0, +1, '$1', 10)'
+# MC EW, -
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kFALSE, 0, -1, '$1', 10)'
+#
+# EDbin 2
+# GRS
+# SC EW, +
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kFALSE, kTRUE, kTRUE, 1, +1, '$1', 10)'
+# SC EW, -
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kFALSE, kTRUE, kTRUE, 1, -1, '$1', 10)'
+# MC EW, +
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kTRUE, 1, +1, '$1', 10)'
+# MC EW, -
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kTRUE, 1, -1, '$1', 10)'
+# Omega0
+# SC EW, +
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kFALSE, kTRUE, kFALSE, 1, +1, '$1', 10)'
+# SC EW, -
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kTRUE, kFALSE, kTRUE, kFALSE, 1, -1, '$1', 10)'
+# MC EW, +
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kFALSE, 1, +1, '$1', 10)'
+# MC EW, -
+root -b -q 'Plot_PDCumulants.C++(kTRUE, kFALSE, kFALSE, kFALSE, kFALSE, kTRUE, kFALSE, 1, -1, '$1', 10)'
\ No newline at end of file
--- /dev/null
+#!/bin/bash
+
+source MakeOutFiles.sh 0
+source MakeOutFiles.sh 1
+source MakeOutFiles.sh 2
+source MakeOutFiles.sh 3
+source MakeOutFiles.sh 4
+source MakeOutFiles.sh 5
+source MakeOutFiles.sh 6
+source MakeOutFiles.sh 7
+source MakeOutFiles.sh 8
+source MakeOutFiles.sh 9
--- /dev/null
+#include <math.h>
+#include <time.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <Riostream.h>
+
+#include "TVector2.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"
+
+using namespace std;
+
+
+void MakeWeightFile()
+{
+
+ gStyle->SetOptStat(0);
+ gStyle->SetOptDate(0);
+
+
+ //TFile *InputFile = new TFile("Results/RawWeightFile_11h.root","READ");
+ //TFile *InputFile = new TFile("Results/RawWeightFile_genSignal_Rinv11.root","READ");
+ //TFile *InputFile = new TFile("Results/RawWeightFile_genSignal_Rinv11_Smeared.root","READ");
+ //TFile *InputFile = new TFile("Results/RawWeightFile_11h_PID1p5.root","READ");
+ //TFile *InputFile = new TFile("Results/RawWeightFile_FB5.root","READ");
+ //TFile *InputFile = new TFile("Results/RawWeightFile_11h_Chi3p1.root","READ");
+ //TFile *InputFile = new TFile("Results/PDC_11h_standard_and_Raw0p04TTC.root","READ");// standard weights
+ //TFile *InputFile = new TFile("Results/RawWeightFile_11h_4ktbins_new.root","READ");
+ TFile *InputFile = new TFile("Results/RawWeightFile_11h_4ktbins_2sigma_3sigmaTTC.root","READ");// _1(2sigmaTTC), _2(3sigmaTTC)
+
+ //TFile *InputFile = new TFile("MyOutput.root","READ");
+ TDirectoryFile *tdir = (TDirectoryFile*)InputFile->Get("PWGCF.outputChaoticityAnalysis.root");
+ TList *MyList=(TList*)tdir->Get("ChaoticityOutput_2");
+ //TList *MyList=(TList*)InputFile->Get("MyList");
+ InputFile->Close();
+
+ const int KtBins=4;
+ const int KyBins=1;
+ const int EDBins=1;
+ const int MBins=10;//10, make sure MB assignment below is correct!!!!!!!!!!!!!!!
+
+ TFile *OutFile = new TFile("WeightFile_temp.root","RECREATE");
+ TH3F *WeightHistos[KtBins][MBins];
+
+ for(int ktB=1; ktB<=KtBins; ktB++){
+ for(int MB=1; MB<=MBins; MB++){
+ TString *InNameNum = new TString("TwoPart_num_Kt_");
+ *InNameNum += ktB-1;
+ InNameNum->Append("_Ky_0_M_");
+ if(MB<=10) *InNameNum += MB-1;// make sure MB assignment is correct!!!!!!!!!!!!!!!
+ else *InNameNum += 1;
+ InNameNum->Append("_ED_0");
+ //
+ TString *InNameDen = new TString("TwoPart_den_Kt_");
+ *InNameDen += ktB-1;
+ InNameDen->Append("_Ky_0_M_");
+ if(MB<=10) *InNameDen += MB-1;// make sure MB assignment is correct!!!!!!!!!!!!!!!
+ else *InNameDen += 1;
+ InNameDen->Append("_ED_0");
+ //
+ TH3D *tempNum = (TH3D*)MyList->FindObject(InNameNum->Data());
+ TH3D *tempDen = (TH3D*)MyList->FindObject(InNameDen->Data());
+ //
+ int NormBinStart = tempNum->GetXaxis()->FindBin(0.135);//0.135
+ int NormBinEnd = tempNum->GetXaxis()->FindBin(0.2);//0.2
+ double Norm = tempNum->Integral(NormBinStart,NormBinEnd, NormBinStart,NormBinEnd, NormBinStart,NormBinEnd);
+ Norm /= tempDen->Integral(NormBinStart,NormBinEnd, NormBinStart,NormBinEnd, NormBinStart,NormBinEnd);
+ cout<<"Normalization = "<<Norm<<endl;
+ cout<<tempNum->Integral(NormBinStart,NormBinEnd, NormBinStart,NormBinEnd, NormBinStart,NormBinEnd)<<" "<<tempDen->Integral(NormBinStart,NormBinEnd, NormBinStart,NormBinEnd, NormBinStart,NormBinEnd)<<endl;
+ //
+ TString *OutNameWeight = new TString("Weight_Kt_");
+ *OutNameWeight += ktB-1;
+ OutNameWeight->Append("_Ky_0_M_");
+ *OutNameWeight += MB-1;
+ OutNameWeight->Append("_ED_0");
+
+ int Nbins=tempNum->GetNbinsX();
+ double QLimit = tempNum->GetXaxis()->GetBinUpEdge(Nbins);
+ WeightHistos[ktB-1][MB-1] = new TH3F(OutNameWeight->Data(),"r3 Weights", Nbins,0,QLimit, Nbins,0,QLimit, Nbins,0,QLimit);
+ WeightHistos[ktB-1][MB-1]->GetXaxis()->SetTitle("out");
+ WeightHistos[ktB-1][MB-1]->GetYaxis()->SetTitle("side");
+ WeightHistos[ktB-1][MB-1]->GetZaxis()->SetTitle("long");
+ WeightHistos[ktB-1][MB-1]->GetXaxis()->SetTitleOffset(1.8);
+ WeightHistos[ktB-1][MB-1]->GetYaxis()->SetTitleOffset(1.8);
+ WeightHistos[ktB-1][MB-1]->GetZaxis()->SetTitleOffset(1.8);
+ for(int outB=1; outB<=Nbins; outB++){
+ for(int sideB=1; sideB<=Nbins; sideB++){
+ for(int longB=1; longB<=Nbins; longB++){
+ double weight=1, weight_e=0;
+ if(tempDen->GetBinContent(outB,sideB,longB) > 0 && tempNum->GetBinContent(outB,sideB,longB) > 0) {
+ //if(outB==sideB && outB==longB) cout<<outB<<" "<<tempNum->GetBinContent(outB,sideB,longB)<<" "<<tempDen->GetBinContent(outB,sideB,longB)<<endl;
+ weight = double(tempNum->GetBinContent(outB,sideB,longB))/double(tempDen->GetBinContent(outB,sideB,longB)) / Norm;
+ weight_e = pow(sqrt(double(tempNum->GetBinContent(outB,sideB,longB)))/double(tempDen->GetBinContent(outB,sideB,longB)) / Norm,2);
+ weight_e += pow(sqrt(double(tempDen->GetBinContent(outB,sideB,longB)))*double(tempNum->GetBinContent(outB,sideB,longB))/pow(double(tempDen->GetBinContent(outB,sideB,longB)),2) / Norm,2);
+ weight_e = sqrt(weight_e);
+ //if(weight < 0.8 && weight > 0.1) cout<<outB<<" "<<sideB<<" "<<longB<<" "<<tempNum->GetBinContent(outB,sideB,longB)<<" "<<tempDen->GetBinContent(outB,sideB,longB)<<endl;
+ }
+ if(weight==0){
+ WeightHistos[ktB-1][MB-1]->SetBinContent(outB,sideB,longB, 0);
+ WeightHistos[ktB-1][MB-1]->SetBinError(outB,sideB,longB, 0);
+ }else {
+ WeightHistos[ktB-1][MB-1]->SetBinContent(outB,sideB,longB, weight-1.0);// difference from unity
+ WeightHistos[ktB-1][MB-1]->SetBinError(outB,sideB,longB, weight_e);
+ }
+
+ }
+ }
+ }
+
+ WeightHistos[ktB-1][MB-1]->Write();
+
+ }// MB
+ }// ktB
+
+ int BOI1=2;
+ int BOI2=2;
+ TH3D *histoOI = WeightHistos[0][0]->Clone();
+ histoOI->SetDirectory(0);
+ //OutFile->Close();
+
+ TH1D *pro = WeightHistos[0][0]->ProjectionY("pro",BOI1,BOI1,BOI2,BOI2);
+ pro->Draw();
+
+ /*TFile *projections=new TFile("projections.root","RECREATE");
+ histoOI->GetXaxis()->SetRange(BOI1,BOI1);
+ TH2D *proYZ = (TH2D*)histoOI->Project3D("yz");
+ histoOI->GetXaxis()->SetRange(1,40);
+ histoOI->GetYaxis()->SetRange(BOI1,BOI1);
+ TH2D *proXZ = (TH2D*)histoOI->Project3D("xz");
+ histoOI->GetYaxis()->SetRange(1,40);
+ histoOI->GetZaxis()->SetRange(BOI1,BOI1);
+ TH2D *proXY = (TH2D*)histoOI->Project3D("xy");
+ proYZ->GetXaxis()->SetRangeUser(0,0.1);
+ proYZ->GetYaxis()->SetRangeUser(0,0.1);
+ proXZ->GetXaxis()->SetRangeUser(0,0.1);
+ proXZ->GetYaxis()->SetRangeUser(0,0.1);
+ proXY->GetXaxis()->SetRangeUser(0,0.1);
+ proXY->GetYaxis()->SetRangeUser(0,0.1);
+ proXY->Draw("lego");
+
+ proYZ->Write();
+ proXZ->Write();
+ proXY->Write();
+ */
+
+}
--- /dev/null
+#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");
+
+}
--- /dev/null
+#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"
+
+#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
+#define fzero_aa 0.942597 // 0.186fm/FmToGeV (scattering length of pi+pi-) = 0.942597
+#define fzero_ab -0.89192 // -0.176fm/FmtoGeV (scattering length (pi+pi-)-->(pi0pi0) = -0.89192
+#define fzero_bb 0.3119 // fzero_aa + 1/sqrt(2)*fzero_ab = 0.0615/FmToGeV = 0.3119
+#define dzero -50.6773 // -10/FmToGeV = -50.6773 (effective range)
+#define EulerC 0.5772156649 // Euler constant
+#define masspi0 0.1349766 // pi0 mass
+#define masspiC 0.1395702 // pi+ mass
+#define massKch 0.493677 // K+- mass
+#define massK0 0.497614 // K0 mass
+
+#define SF 1.0 // Scale Factor for spatial coordinates
+#define SF_Mom 1.0 // Scale Factor to reduce momenta
+#define RstarMax 405.4/SF // 405.4 (80 fm / FmToGeV), tried 253.4 (50 fm / FmToGeV) as a variation
+#define RstarMin 0.507/SF // 0.507 (0.1 fm / FmToGeV)
+using namespace std;
+
+// g and p used in Gamma function
+const int g_gamma = 7;
+double p_gamma[9] = {0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7};
+// h coefficients for very small kstar (eta > 0.2)
+double hcoef[15]={.8333333333e-1, .8333333333e-2, .396825396825e-2, .4166666667e-2, .7575757576e-2, .2109279609e-1, .8333333333e-1, .4432598039e0 , .305395433e1, .2645621212e2, .2814601449e3, .3607510546e4, .5482758333e5, .9749368235e6, .200526958e8};
+// Scattering length (in GeV) coefficients in isospin representation from Lednicky's code from Colangelo (2001)
+double fzero[2][6]={{.220, .268,-.0139,-.00139, 36.77, 15*0}, {-.0444,-.0857,-.00221,-.000129,-21.62, 15*0}};
+
+double massPOI;
+int PIDPOI;
+double qCutOff;
+
+float Gamov(float);
+float fact(float);
+double hFunction(double);
+double KinverseFunction(double, int);
+void Shuffle(int*, int, int);
+void BoostPRF(float [4], float [4], float [4], TVector3*);
+
+
+complex<double> gamma(complex<double> z)
+{// Lanczos approximation
+
+ if(floor(abs(z)) == z)
+ return fact(abs(z)-1.0);
+
+ if(real(z)<(double)(0.5)) return PI/(sin(PI*z)*gamma((1.0-z)));
+ z -= 1.0;
+ complex<double> x = p_gamma[0];
+ for(int i=1;i<g_gamma+2;i++) x = x+p_gamma[i]/(z+(double)i);
+ complex<double> t = z+(double)g_gamma+0.5;
+ return sqrt(2*PI)*pow(t,z+0.5)*exp(-t)*x;
+}
+
+// confluent hypergeometric function
+complex<double> conf_hyp(complex<double> a, complex<double> b, complex<double> z)
+{// b is assumed to be 1.0 here
+ complex<double> S(1.0, 0.0);
+ complex<double> F=S;
+ complex<double> alf1 = a - 1.0;
+ double J=0;
+ double errorCutOff=0.00001;// series truncation point (0.00001, 1e-6 for strict case)
+ if(abs(z) > 40.) {// high k*r* case (was 20., now 40. Improved accuracy for qinv>60 MeV)
+ double eta = -a.imag();
+ double RKS = z.imag();
+ double D1 = log(RKS)*eta;
+ double D2 = pow(eta,2)/RKS;
+ complex<double> arg(1.0, eta);
+ complex<double> EIDC = gamma(arg);
+ EIDC /= abs(EIDC);
+ complex<double> ZZ1(cos(D1), sin(D1));
+ ZZ1 /= EIDC;
+ complex<double> value(1.0, -D2);
+ value *= ZZ1;
+ complex<double> value2(cos(RKS), sin(RKS));
+ F = value - value2*eta/RKS/ZZ1;
+ return F;
+ }else {// low k*r* case
+ while(fabs(S.real())+fabs(S.imag()) > errorCutOff){
+ J++;
+ complex<double> A = (alf1 + J)/(pow(J,2));
+ S *= A*z;
+ F += S;
+ }
+ return F;
+ }
+}
+
+void Therm(){
+
+ TRandom3 *RandomNumber = new TRandom3();
+ RandomNumber->SetSeed(0);
+
+
+ bool StrongFSI=kTRUE;
+ bool LambdaDirect=kFALSE;
+ bool RemoveXP=kFALSE;
+ bool Gaussian=kFALSE;
+ bool ThreeParticle=kFALSE;
+ bool BoostPairsIntoPRF=kTRUE;
+ massPOI = masspiC;
+ PIDPOI = 211;// 211 for pi+-
+ qCutOff = 0.1;// 0.1 or 0.4 for pp
+ int qBins = int(qCutOff/0.002);
+ //
+ int Nfiles=1000;// 1000 for b2, 500 b3-b5, 1500 b7-b8, 3000 b9
+ int bValue=2;
+
+ double RValue=8;// Gaussian radius
+ double NsigmaGauss=2;// number of sigma to integrate over for a Gaussian source
+ if(bValue==2) RValue=8;//8
+ if(bValue==3) RValue=8;//8
+ if(bValue==5) RValue=7;//7
+ if(bValue==7) RValue=6;//6
+ if(bValue==8) RValue=5;//5
+ if(bValue==9) RValue=5;//5
+ TF1 *SingleParticleGaussian = new TF1("SingleParticleGaussian","exp(-pow(x/(sqrt(2.)*[0]),2))",0,1000);
+ SingleParticleGaussian->SetParameter(0,RValue);
+
+ TFile *outfile = new TFile("Therm_dist_temp.root","RECREATE");
+ TH1D *r_dist=new TH1D("r_dist","",400,0,100);
+ TH1D *x_dist=new TH1D("x_dist","",400,0,100);
+ TH1D *Pt_dist=new TH1D("Pt_dist","",200,0,2);
+ TH1D *P_dist=new TH1D("P_dist","",900,0.1,1);
+ TH1D *P_dist_lowq=new TH1D("P_dist_lowq","",900,0.1,1);
+ TProfile *AvgMult = new TProfile("AvgMult","",1,0.5,1.5, 0,2000,"");
+ //
+ TH2D *rstarPRF_dist_ss=new TH2D("rstarPRF_dist_ss","",20,0,1.0, 400,0,100);
+ TH2D *tstar_dist_ss=new TH2D("tstar_dist_ss","",20,0,1.0, 400,0,100);
+ //
+ TH2D *rstarPRF_dist_os=new TH2D("rstarPRF_dist_os","",20,0,1.0, 400,0,100);
+ TH2D *tstar_dist_os=new TH2D("tstar_dist_os","",20,0,1.0, 400,0,100);
+
+ TH2D *Num_Cos_ss=new TH2D("Num_Cos_ss","",20,0,1, 40,0,0.2);
+ TH2D *Num_Cos_os=new TH2D("Num_Cos_os","",20,0,1, 40,0,0.2);
+ TH2D *Num_CosFSI_ss=new TH2D("Num_CosFSI_ss","",20,0,1, 40,0,0.2);
+ TH2D *Num_CosFSI_os=new TH2D("Num_CosFSI_os","",20,0,1, 40,0,0.2);
+ TH2D *NumSq_CosFSI_ss=new TH2D("NumSq_CosFSI_ss","",20,0,1, 40,0,0.2);
+ TH2D *NumSq_CosFSI_os=new TH2D("NumSq_CosFSI_os","",20,0,1, 40,0,0.2);
+ TH2D *Num_PrimCosFSI_ss=new TH2D("Num_PrimCosFSI_ss","",20,0,1, 40,0,0.2);
+ TH2D *Num_PrimCosFSI_os=new TH2D("Num_PrimCosFSI_os","",20,0,1, 40,0,0.2);
+ //
+ TH2D *Den_ss=new TH2D("Den_ss","",20,0,1, 40,0,0.2);
+ TH2D *Den_os=new TH2D("Den_os","",20,0,1, 40,0,0.2);
+ TH2D *DenEM_ss=new TH2D("DenEM_ss","",20,0,1, 40,0,0.2);
+ TH2D *DenEM_os=new TH2D("DenEM_os","",20,0,1, 40,0,0.2);
+ TH2D *LargeRpairs_ss=new TH2D("LargeRpairs_ss","",20,0,1, 40,0,0.2);
+ TH2D *LargeRpairs_os=new TH2D("LargeRpairs_os","",20,0,1, 40,0,0.2);
+ //
+ TH2D *NumLambda1=new TH2D("NumLambda1","",40,0,0.2, 100,0,100);
+ TH1D *DenLambda1=new TH1D("DenLambda1","",40,0,0.2);
+ TH3D *NumLambda2_ss=new TH3D("NumLambda2_ss","",20,0,1, 40,0,0.2, 100,0,100);
+ TH3D *NumLambda2_os=new TH3D("NumLambda2_os","",20,0,1, 40,0,0.2, 100,0,100);
+ TH2D *DenLambda2_ss=new TH2D("DenLambda2_ss","",20,0,1, 40,0,0.2);
+ TH2D *DenLambda2_os=new TH2D("DenLambda2_os","",20,0,1, 40,0,0.2);
+ TH3D *NumLambda3_ss=new TH3D("NumLambda3_ss","",40,0,.2, 40,0,.2, 40,0,.2);
+ TH3D *DenLambda3_ss=new TH3D("DenLambda3_ss","",40,0,.2, 40,0,.2, 40,0,.2);
+ TH3D *NumLambda3_os=new TH3D("NumLambda3_os","",40,0,.2, 40,0,.2, 40,0,.2);
+ TH3D *DenLambda3_os=new TH3D("DenLambda3_os","",40,0,.2, 40,0,.2, 40,0,.2);
+ //
+ TH3D *rstar3_dist_ss=new TH3D("rstar3_dist_ss","",200,0,100, 200,0,100, 200,0,100);
+ TH3D *tstar3_dist_ss=new TH3D("tstar3_dist_ss","",200,0,100, 200,0,100, 200,0,100);
+ //
+ TH3D *r3num=new TH3D("r3num","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ TH3D *r3den1=new TH3D("r3den1","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ TH3D *r3den2=new TH3D("r3den2","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ TH3D *r3den3=new TH3D("r3den3","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ TH3D *r3numSq=new TH3D("r3numSq","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ TH3D *r3den1Sq=new TH3D("r3den1Sq","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ TH3D *r3den2Sq=new TH3D("r3den2Sq","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ TH3D *r3den3Sq=new TH3D("r3den3Sq","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ TH3D *r3numEn=new TH3D("r3numEn","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ TH3D *r3den1En=new TH3D("r3den1En","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ TH3D *r3den2En=new TH3D("r3den2En","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ TH3D *r3den3En=new TH3D("r3den3En","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ TH3D *r3numME=new TH3D("r3numME","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ TH3D *r3den1ME=new TH3D("r3den1ME","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ TH3D *r3den2ME=new TH3D("r3den2ME","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ TH3D *r3den3ME=new TH3D("r3den3ME","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ //
+ TH1D *Num_Cos12=new TH1D("Num_Cos12","",40,0,0.2);
+ TH1D *Num_Cos23=new TH1D("Num_Cos23","",40,0,0.2);
+ TH1D *Num_Cos31=new TH1D("Num_Cos31","",40,0,0.2);
+ TH1D *Den_Cos12=new TH1D("Den_Cos12","",40,0,0.2);
+ TH1D *Den_Cos23=new TH1D("Den_Cos23","",40,0,0.2);
+ TH1D *Den_Cos31=new TH1D("Den_Cos31","",40,0,0.2);
+ //
+ //
+ TH2D *K2_ss = new TH2D("K2_ss","",20,0,1, qBins,0,qCutOff);
+ TH2D *PlaneWF_ss = new TH2D("PlaneWF_ss","",20,0,1, qBins,0,qCutOff);
+ TH2D *K2_os = new TH2D("K2_os","",20,0,1, qBins,0,qCutOff);
+ TH2D *PlaneWF_os = new TH2D("PlaneWF_os","",20,0,1, qBins,0,qCutOff);
+ //
+ TH1D *PlaneWF3ss = new TH1D("PlaneWF3ss","",qBins,0,qCutOff);
+ TH1D *PlaneWF3os = new TH1D("PlaneWF3os","",qBins,0,qCutOff);
+ TH1D *K3ss = new TH1D("K3ss","",qBins,0,qCutOff);
+ K3ss->GetXaxis()->SetTitle("Q3 (GeV/c)");
+ K3ss->GetYaxis()->SetTitle("K_{3}");
+ TH1D *K3os = new TH1D("K3os","",qBins,0,qCutOff);
+ K3os->GetXaxis()->SetTitle("Q3 (GeV/c)");
+ K3os->GetYaxis()->SetTitle("K_{3}");
+ TH3D *K3ss_3D = new TH3D("K3ss_3D","", qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ K3ss_3D->GetXaxis()->SetTitle("q12");
+ K3ss_3D->GetYaxis()->SetTitle("q23");
+ K3ss_3D->GetZaxis()->SetTitle("q31");
+ TH3D *PlaneWF3ss_3D = new TH3D("PlaneWF3ss_3D","", qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ TH3D *K3os_3D = new TH3D("K3os_3D","", qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+ K3os_3D->GetXaxis()->SetTitle("q12");
+ K3os_3D->GetYaxis()->SetTitle("q23");
+ K3os_3D->GetZaxis()->SetTitle("q31");
+ TH3D *PlaneWF3os_3D = new TH3D("PlaneWF3os_3D","", qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
+
+ TH1D *test=new TH1D("test","",200,-2*PI,2*PI);
+
+ ////////////////////////
+
+
+ complex<double> binput(1.0,0.0);
+ complex<double> CHG(0.,0.);// for Strong FSI in 3-body WF
+ double kstar=0;
+ double rstar=0;
+ double kstar_1=0, kstar_2=0;
+ double rstar_1=0, rstar_2=0;
+ double phase_1=0, phase_2=0, phase_3=0, phase_4=0;
+ double qinvMixedC_ss=0, qinvMixedC_os1=0, qinvMixedC_os2=0;
+ double WeightWF_aa=1.0, WeightWF_bb=1.0;
+ complex<double> CHG_1(0.,0.);
+ complex<double> CHG_2(0.,0.);
+ complex<double> CHG_3(0.,0.);
+ complex<double> CHG_4(0.,0.);
+ complex<double> AlphaWF[4];
+ complex<double> BetaWF[4];
+ double ThreePhase_1=0, ThreePhase_2=0, ThreePhase_3=0;
+ double ThreePhase_4=0, ThreePhase_5=0, ThreePhase_6=0;
+ double TwoPhase_12=0, TwoPhase_23=0, TwoPhase_31=0;
+
+ const int MAXPIONS = 2000;
+ const int NevtBuff = 1;
+ float P[NevtBuff][MAXPIONS][4]={{{0}}};
+ float X[NevtBuff][MAXPIONS][4]={{{0}}};
+ int PID[NevtBuff][MAXPIONS]={{0}};
+ bool Primary[NevtBuff][MAXPIONS]={{0}};
+ int Nparticles[NevtBuff]={0};
+ Int_t ShuffledIndex[MAXPIONS]={0};// x-p decorrelation
+ TVector3 P1_PRF;
+ TVector3 P2_PRF;
+ TVector3 P3_PRF;
+ TVector3 X1_PRF;
+ TVector3 X2_PRF;
+ TVector3 X3_PRF;
+ TVector3 X1_RF;
+ TVector3 X2_RF;
+ float P_dummy[4]={0};
+
+
+
+ // First Create B(rho,eta) and P(rho,eta) and h(eta) functions (for G hypergeometric functions (strong+Coulomb FSI))
+ const int Btermslimit=500;// 500
+ double rhoLimit = .05*100/FmToGeV;// 100 fm "limit"
+ int rhoBins = rhoLimit/(0.001/FmToGeV);// this binning provides small changes between nearby bins
+ double rhoStep = rhoLimit/rhoBins;
+ double etaLimit = 1/(0.0005*BohrR);
+ int etaBins = 2.0*etaLimit/(0.001);// this binning provides small changes between nearby bins
+ double etaStep = (2*etaLimit)/etaBins;
+ //cout<<rhoBins<<" "<<etaBins<<" "<<rhoLimit/rhoBins<<" "<<etaLimit/etaBins<<endl;
+ TH2D *B_histogram = new TH2D("B_histogram","",rhoBins,0,rhoLimit, etaBins,-etaLimit,etaLimit);
+ TH2D *P_histogram = new TH2D("P_histogram","",rhoBins,0,rhoLimit, etaBins,-etaLimit,etaLimit);
+ TH1D *h_histogram = new TH1D("h_histogram","",etaBins,-etaLimit,etaLimit);
+
+ for(int i=1; i<=rhoBins; i++){
+ for(int j=1; j<=etaBins; j++){
+ double rho = (i+0.5)*rhoStep;
+ double eta = (j+0.5)*etaStep - etaLimit;// starts from negative values
+ if(fabs(eta) < 0.0001) continue;
+
+ double Bfunc[Btermslimit]={0};
+ double Pfunc[Btermslimit]={0};
+ int StopPoint=Btermslimit;
+ Bfunc[0]=1; Bfunc[1]=eta*rho;
+ Pfunc[0]=1; Pfunc[1]=0;
+ for(int ii=2; ii<Btermslimit; ii++) {
+ Bfunc[ii] = (2*eta*rho*Bfunc[ii-1] - rho*rho*Bfunc[ii-2])/(double(ii*(ii+1.)));
+ Pfunc[ii] = (2*eta*rho*Pfunc[ii-1] - rho*rho*Pfunc[ii-2] - (2*(ii-1)+1)*2*eta*rho*Bfunc[ii-1])/(double((ii-1)*ii));
+ if(fabs(Bfunc[ii])+fabs(Bfunc[ii-1])+fabs(Pfunc[ii])+fabs(Pfunc[ii-1]) < 1.0e-7) {StopPoint=ii; break;}
+ }
+ //cout<<StopPoint+1<<endl;
+ double Bsum=0, Psum=0;
+ for(int ii=0; ii<StopPoint+1; ii++) {Bsum += Bfunc[ii]; Psum += Pfunc[ii];}
+ B_histogram->SetBinContent(i,j, Bsum);
+ P_histogram->SetBinContent(i,j, Psum);
+ }
+ }
+ cout<<"Done preparing B and P histograms"<<endl;
+ for(int j=1; j<=etaBins; j++){
+ double eta = (j+0.5)*etaStep - etaLimit;// starts from negative values
+ if(fabs(eta) < 0.0001) continue;
+ h_histogram->SetBinContent(j, hFunction(eta));
+ }
+ cout<<"Done preparing h histogram"<<endl;
+
+
+ int seconds = time(0);
+
+ for(int FileN=0; FileN<Nfiles; FileN++){
+ cout<<"File #"<<FileN+1<<endl;
+ //if(FileN<1260) continue;
+ //if(FileN>230) continue;
+ TString *fname=new TString("b");
+ *fname += bValue;
+ fname->Append("/event");
+ *fname += FileN+1;
+ fname->Append(".root");
+ TFile *infile=new TFile(fname->Data(),"READ");
+ if(!infile->IsOpen()) {cout<<fname->Data()<<" does not exist"<<endl; continue;}
+ TTree *event_tree=(TTree*)infile->Get("events");
+ TTree *particles_tree=(TTree*)infile->Get("particles");
+
+
+ int Nevents = event_tree->GetEntries();
+ int NpartList=0;
+ int NpartList_past = 0;
+ /////////////////////////////////////////
+ // Create Pion list first
+ for(int Evt=0; Evt<Nevents; Evt++){
+
+ //cout<<"Event # "<<Evt+1<<endl;
+ event_tree->GetEntry(Evt);
+
+ TBranch *br=(TBranch*)event_tree->GetBranch("event");
+ TLeaf *lf=(TLeaf*)br->GetLeaf("entries");
+ NpartList_past += NpartList;
+ NpartList = lf->GetValue();
+
+ /////////////////////////////////////
+ // past event buffer
+ for(int pastEvt=NevtBuff-1; pastEvt>0; pastEvt--){
+
+ for(int particle=0; particle<Nparticles[pastEvt-1]; particle++){
+ P[pastEvt][particle][0] = P[pastEvt-1][particle][0];
+ P[pastEvt][particle][1] = P[pastEvt-1][particle][1];
+ P[pastEvt][particle][2] = P[pastEvt-1][particle][2];
+ P[pastEvt][particle][3] = P[pastEvt-1][particle][3];
+ //
+ X[pastEvt][particle][0] = X[pastEvt-1][particle][0];
+ X[pastEvt][particle][1] = X[pastEvt-1][particle][1];
+ X[pastEvt][particle][2] = X[pastEvt-1][particle][2];
+ X[pastEvt][particle][3] = X[pastEvt-1][particle][3];
+ //
+ PID[pastEvt][particle] = PID[pastEvt-1][particle];
+ Primary[pastEvt][particle] = Primary[pastEvt-1][particle];
+ }
+
+ Nparticles[pastEvt]=Nparticles[pastEvt-1];
+ }
+ Nparticles[0]=0;
+ /////////////////////////////////////
+
+
+ // Create Pion list first
+ int Npions=0;
+
+ for(int index1=NpartList_past; index1<NpartList_past + NpartList; index1++){
+
+ if(Npions >= MAXPIONS) {cout<<"Too Many Pions!!!"<<endl; break;}
+ particles_tree->GetEntry(index1);
+ //
+
+ TLeaf *lf_pid=particles_tree->GetLeaf("pid");
+ int pid = int(lf_pid->GetValue());
+ if(abs(pid) != PIDPOI) continue;// pions only
+ TLeaf *lf_decayed=(TLeaf*)particles_tree->GetLeaf("decayed");
+ if(lf_decayed->GetValue() != 0) continue;
+
+
+ TLeaf *lf_px=(TLeaf*)particles_tree->GetLeaf("px");
+ TLeaf *lf_py=(TLeaf*)particles_tree->GetLeaf("py");
+ TLeaf *lf_pz=(TLeaf*)particles_tree->GetLeaf("pz");
+ TLeaf *lf_x=(TLeaf*)particles_tree->GetLeaf("x");
+ TLeaf *lf_y=(TLeaf*)particles_tree->GetLeaf("y");
+ TLeaf *lf_z=(TLeaf*)particles_tree->GetLeaf("z");
+ TLeaf *lf_t=(TLeaf*)particles_tree->GetLeaf("t");
+
+ float px = lf_px->GetValue(); float py = lf_py->GetValue(); float pz = lf_pz->GetValue();
+ float x = lf_x->GetValue(); float y = lf_y->GetValue(); float z = lf_z->GetValue();
+ double t = lf_t->GetValue();
+
+ if(Gaussian){// Gaussian randomization
+ bool goodchoice=kFALSE;
+ while(goodchoice==kFALSE){
+ x = NsigmaGauss*RValue*RandomNumber->Rndm();
+ float height = RandomNumber->Rndm();
+ if(SingleParticleGaussian->Eval(x) >= height) goodchoice=kTRUE;
+ }
+ if(RandomNumber->Rndm()<0.5) x = -x;
+ //
+ goodchoice=kFALSE;
+ while(goodchoice==kFALSE){
+ y = NsigmaGauss*RValue*RandomNumber->Rndm();
+ float height = RandomNumber->Rndm();
+ if(SingleParticleGaussian->Eval(y) >= height) goodchoice=kTRUE;
+ }
+ if(RandomNumber->Rndm()<0.5) y = -y;
+ //
+ goodchoice=kFALSE;
+ while(goodchoice==kFALSE){
+ z = NsigmaGauss*RValue*RandomNumber->Rndm();
+ float height = RandomNumber->Rndm();
+ if(SingleParticleGaussian->Eval(z) >= height) goodchoice=kTRUE;
+ }
+ if(RandomNumber->Rndm()<0.5) z = -z;
+ //
+ goodchoice=kFALSE;
+ while(goodchoice==kFALSE){
+ t = NsigmaGauss*RValue*RandomNumber->Rndm();
+ float height = RandomNumber->Rndm();
+ if(SingleParticleGaussian->Eval(t) >= height) goodchoice=kTRUE;
+ }
+ if(RandomNumber->Rndm()<0.5) t = -t;
+ //t=0;
+ }
+
+ float pt = sqrt(pow(px,2)+pow(py,2));
+ float eta = asinh(pz/pt);
+ float E = sqrt(pow(px,2)+pow(py,2)+pow(pz,2)+pow(massPOI,2));
+
+ if(pt<0.16 || pt>1.0) continue;
+ if(fabs(eta)>0.8) continue;
+ if(sqrt(pow(px,2)+pow(py,2)+pow(pz,2)) > 1.0) continue;
+ //
+ P[0][Npions][0] = E/SF_Mom; P[0][Npions][1] = px/SF_Mom; P[0][Npions][2] = py/SF_Mom; P[0][Npions][3] = pz/SF_Mom;
+ X[0][Npions][0] = t/FmToGeV/SF; X[0][Npions][1] = x/FmToGeV/SF; X[0][Npions][2] = y/FmToGeV/SF; X[0][Npions][3] = z/FmToGeV/SF;
+ PID[0][Npions] = pid;
+ float r = sqrt(pow(x,2)+pow(y,2)+pow(z,2));
+ if(t < 1.e6) Primary[0][Npions]=kTRUE;// 1e6
+ else {Primary[0][Npions]=kFALSE;}
+
+
+ if(RemoveXP){
+ if(r > 100./SF) continue;
+ }
+
+ if(Primary[0][Npions]) P_dist->Fill(sqrt( pow(px,2) + pow(py,2) + pow(pz,2)));
+
+ if(Primary[0][Npions]){
+ r_dist->Fill(r);
+ x_dist->Fill(fabs(X[0][Npions][1])*FmToGeV);
+ Pt_dist->Fill(pt);
+ }
+
+ //if(Primary[0][Npions]==kFALSE) continue;
+
+ //
+
+
+ Npions++;
+ }
+
+ Nparticles[0]=Npions;
+ AvgMult->Fill(1, Npions);
+
+ // Shuffle
+ for (Int_t i = 0; i < Npions; i++) {
+ // generate random index call
+ int j = (Int_t) (gRandom->Rndm() * Npions);
+ // temporarily store values from random index
+ float t = X[0][j][0], x = X[0][j][1], y = X[0][j][2], z = X[0][j][3];
+ float E = P[0][j][0], px = P[0][j][1], py = P[0][j][2], pz = P[0][j][3];
+ int pid = PID[0][j];
+ bool prim = Primary[0][j];
+ //
+ // Swap value locations
+ X[0][j][0] = X[0][i][0]; X[0][j][1] = X[0][i][1]; X[0][j][2] = X[0][i][2]; X[0][j][3] = X[0][i][3];
+ P[0][j][0] = P[0][i][0]; P[0][j][1] = P[0][i][1]; P[0][j][2] = P[0][i][2]; P[0][j][3] = P[0][i][3];
+ PID[0][j] = PID[0][i];
+ Primary[0][j] = Primary[0][i];
+ //
+ X[0][i][0] = t; X[0][i][1] = x; X[0][i][2] = y; X[0][i][3] = z;
+ P[0][i][0] = E; P[0][i][1] = px; P[0][i][2] = py; P[0][i][3] = pz;
+ PID[0][i] = pid;
+ Primary[0][i] = prim;
+ }
+
+ ///////////////////////////////////////////
+ // de-correlate x-p
+ if(RemoveXP){
+ for (Int_t i = 0; i < Npions; i++) ShuffledIndex[i] = i;
+ Shuffle(ShuffledIndex, 0, Npions-1);
+ for (Int_t i = 0; i < Npions; i++) {// decorrelate x-p
+ float t = X[0][i][0], x = X[0][i][1], y = X[0][i][2], z = X[0][i][3];
+ X[0][i][0] = X[0][ShuffledIndex[i]][0];
+ X[0][i][1] = X[0][ShuffledIndex[i]][1];
+ X[0][i][2] = X[0][ShuffledIndex[i]][2];
+ X[0][i][3] = X[0][ShuffledIndex[i]][3];
+ //
+ X[0][ShuffledIndex[i]][0] = t;
+ X[0][ShuffledIndex[i]][1] = x;
+ X[0][ShuffledIndex[i]][2] = y;
+ X[0][ShuffledIndex[i]][3] = z;
+ }
+ }
+
+
+ //////////////////////////////////////////
+
+
+
+ //////////////////////////////////////////////////
+ // Start Correlation Analysis
+
+ for(int EN=0; EN<1; EN++){// current or 1st past event
+
+ for(int index1=0; index1<Nparticles[EN]; index1++){// current or past event pion loop
+
+ int start=index1+1;
+ if(EN>0) start=0;
+ for(int index2=start; index2<Npions; index2++){// current event pion loop (EN=0)
+
+ float kt = sqrt(pow(P[EN][index1][1]+P[0][index2][1],2)+pow(P[EN][index1][2]+P[0][index2][2],2))/2.;
+ float qinv12 = sqrt(pow(P[EN][index1][1]-P[0][index2][1],2)+pow(P[EN][index1][2]-P[0][index2][2],2)+pow(P[EN][index1][3]-P[0][index2][3],2)-pow(P[EN][index1][0]-P[0][index2][0],2));
+ float rstar12_CMS = sqrt(pow(X[EN][index1][1]-X[0][index2][1],2)+pow(X[EN][index1][2]-X[0][index2][2],2)+pow(X[EN][index1][3]-X[0][index2][3],2));
+
+ if(qinv12 < 0.001) continue;
+ if(qinv12 > qCutOff) continue;
+
+ if(EN>0){
+ if(PID[EN][index1] == PID[0][index2]){
+ DenEM_ss->Fill(kt, qinv12);
+ }else{
+ DenEM_os->Fill(kt, qinv12);
+ }
+ if(LambdaDirect) continue;
+ }
+
+ BoostPRF(P[EN][index1], P_dummy, X[EN][index1], &X1_RF);
+ DenLambda1->Fill(qinv12);
+ if(X1_RF.Mag()<RstarMax) NumLambda1->Fill(qinv12, X1_RF.Mag()*FmToGeV);
+ BoostPRF(P[0][index2], P_dummy, X[0][index2], &X2_RF);
+ DenLambda1->Fill(qinv12);
+ if(X2_RF.Mag()<RstarMax) NumLambda1->Fill(qinv12, X2_RF.Mag()*FmToGeV);
+
+
+ P_dist_lowq->Fill(sqrt(pow(P[EN][index1][1],2)+pow(P[EN][index1][2],2)+pow(P[EN][index1][3],2)));
+ P_dist_lowq->Fill(sqrt(pow(P[0][index2][1],2)+pow(P[0][index2][2],2)+pow(P[0][index2][3],2)));
+
+ // boost P into PRF
+ BoostPRF(P[EN][index1], P[0][index2], P[EN][index1], &P1_PRF);
+ BoostPRF(P[EN][index1], P[0][index2], P[0][index2], &P2_PRF);
+ if(BoostPairsIntoPRF){// boost X into PRF
+ BoostPRF(P[EN][index1], P[0][index2], X[EN][index1], &X1_PRF);
+ BoostPRF(P[EN][index1], P[0][index2], X[0][index2], &X2_PRF);
+ }else{
+ X1_PRF[0]=X[EN][index1][1]; X1_PRF[1]=X[EN][index1][2]; X1_PRF[2]=X[EN][index1][3];
+ X2_PRF[0]=X[0][index2][1]; X2_PRF[1]=X[0][index2][2]; X2_PRF[2]=X[0][index2][3];
+ }
+ //cout<<"out: "<<P1_PRF[0]<<" "<<P2_PRF[0]<<endl;
+ //cout<<"side: "<<P1_PRF[1]<<" "<<P2_PRF[1]<<endl;
+ //cout<<"long: "<<P1_PRF[2]<<" "<<P2_PRF[2]<<" ++++++"<<endl;
+ //cout<<qinv12<<" "<<sqrt(pow(P1_PRF[0]-P2_PRF[0],2)+pow(P1_PRF[1]-P2_PRF[1],2)+pow(P1_PRF[2]-P2_PRF[2],2))<<endl;
+ //
+ double rstar12_PRF = (X1_PRF-X2_PRF).Mag();
+
+ if(PID[EN][index1] == PID[0][index2]) DenLambda2_ss->Fill(kt, qinv12);
+ else DenLambda2_os->Fill(kt, qinv12);
+
+ if(rstar12_PRF < RstarMin) continue;// removes most pairs from the same decay
+
+ if(rstar12_PRF < RstarMax){
+ if(PID[EN][index1] == PID[0][index2]) NumLambda2_ss->Fill(kt, qinv12, rstar12_PRF*FmToGeV);
+ else NumLambda2_os->Fill(kt, qinv12, rstar12_PRF*FmToGeV);
+ }
+
+ if(rstar12_PRF > RstarMax) {
+ if(PID[EN][index1]==PID[0][index2]) LargeRpairs_ss->Fill(kt, qinv12);
+ else LargeRpairs_os->Fill(kt, qinv12);
+ if(!LambdaDirect) continue;
+ }
+ float FVP12 = (P1_PRF-P2_PRF).Mag();
+ FVP12 -= pow(P[EN][index1][0]-P[0][index2][0],2);
+ double tstar12_PRF = sqrt(fabs(pow(rstar12_PRF,2)-FVP12));// fabs biases only 0.00001 % of the pairs (likely from same decay)
+
+
+
+ //
+ double rho = rstar12_PRF * qinv12/2.;
+ double phase = (P1_PRF-P2_PRF)*(X1_PRF-X2_PRF);
+ double Cosphase = phase/(rstar12_PRF)/(qinv12);
+
+ double eta = 1/(qinv12/2.*BohrR);
+ if(PID[EN][index1] != PID[0][index2]) eta = -eta;
+ complex<double> ainput(0.0, -eta);
+ complex<double> zinput12(0.0, rho + rho*Cosphase);
+ complex<double> zinput21(0.0, rho - rho*Cosphase);
+ complex<double> CHG12(1.0, 0);
+ complex<double> CHG21(1.0, 0);
+ double GamovFactor=1.0;
+ double BEE = 0.;
+ double WeightWF = 1.0;
+ double WeightPlaneWF = 1.0;
+ // Strong FSI
+ complex<double> Psi_alpha(0., 0.);
+ if(EN==0 && !LambdaDirect && qinv12 < qCutOff && rstar12_PRF < RstarMax){
+ if(Primary[EN][index1]==kTRUE && Primary[0][index2]==kTRUE){
+ if(StrongFSI && PID[EN][index1]!=PID[0][index2]){
+ double kb = sqrt(pow(massPOI,2)+pow(qinv12/2.,2) - pow(masspi0,2));
+ double G2 = KinverseFunction(pow(qinv12/2.,2),1);
+ double G1 = KinverseFunction(pow(qinv12/2.,2),0);
+ double RK11 = 2/3.*G1 + 1/3.*G2;
+ double RK22 = 2/3.*G2 + 1/3.*G1;
+ double RK12 = -sqrt(2/9.)*(G1-G2);
+ complex<double> Chi(h_histogram->GetBinContent(h_histogram->GetXaxis()->FindBin(eta)), Gamov(eta)/(2*eta));
+ complex<double> C3 = RK11 - 2.0*Chi/BohrR;
+ complex<double> C5(RK22, -kb);
+ complex<double> C10 = C3*C5 - pow(RK12,2);
+ complex<double> fc_aa = C5/C10;
+ complex<double> G = P_histogram->GetBinContent(P_histogram->GetXaxis()->FindBin(rho), P_histogram->GetYaxis()->FindBin(eta));
+ G += 2.0*eta*rho*( log(fabs(2.0*eta*rho)) + 2.0*EulerC - 1.0 + Chi )*B_histogram->GetBinContent(B_histogram->GetXaxis()->FindBin(rho), B_histogram->GetYaxis()->FindBin(eta));
+ Psi_alpha = fc_aa*G/(rstar12_PRF);
+ }
+
+ CHG12= conf_hyp(ainput,binput,zinput12);
+ CHG21= conf_hyp(ainput,binput,zinput21);
+ GamovFactor = Gamov(eta);
+ BEE = cos(phase);
+
+ complex<double> planewave12(double(TMath::Cos(-rho*Cosphase)), double(TMath::Sin(-rho*Cosphase)));
+ complex<double> planewave21(planewave12.real(), -planewave12.imag());
+ //
+ complex<double> Full12(0.0,0.0);
+ Full12 += planewave12;
+ Full12 *= CHG12;
+ if(PID[EN][index1]!=PID[0][index2]) Full12 += Psi_alpha;
+ //
+ complex<double> Full21(0.0,0.0);
+ Full21 += planewave21;
+ Full21 *= CHG21;
+ //
+ complex<double> FullWF = Full12;
+ complex<double> FullPlaneWF = planewave12;
+ if(PID[EN][index1] == PID[0][index2]) {
+ FullWF += Full21;
+ FullWF /= sqrt(2.);
+ FullPlaneWF += planewave21;
+ FullPlaneWF /= sqrt(2.);
+ }
+
+ WeightWF = GamovFactor*pow(abs(FullWF),2);
+ WeightPlaneWF = pow(abs(FullPlaneWF),2);
+
+ // Fill undiluted histos
+ if(PID[EN][index1] == PID[0][index2]){
+ Num_PrimCosFSI_ss->Fill(kt, qinv12, WeightWF);
+ }else {
+ Num_PrimCosFSI_os->Fill(kt, qinv12, WeightWF);
+ }
+
+ }// primary
+ }// same-event, LambdaDirect, rstar and qinv cut
+
+
+ if(PID[EN][index1] == PID[0][index2]){// same-charge
+ rstarPRF_dist_ss->Fill(kt, rstar12_PRF*FmToGeV);
+ tstar_dist_ss->Fill(kt, fabs(X[EN][index1][0]-X[0][index2][0])*FmToGeV);
+ Num_Cos_ss->Fill(kt, qinv12, BEE);
+ Num_CosFSI_ss->Fill(kt, qinv12, WeightWF);
+ NumSq_CosFSI_ss->Fill(kt, qinv12, pow(WeightWF,2));
+ Den_ss->Fill(kt, qinv12);
+ //
+ K2_ss->Fill(kt, qinv12, WeightWF);
+ PlaneWF_ss->Fill(kt, qinv12, WeightPlaneWF);
+ }else {// mixed-charge
+ rstarPRF_dist_os->Fill(kt, rstar12_PRF*FmToGeV);
+ tstar_dist_os->Fill(kt, fabs(X[EN][index1][0]-X[0][index2][0])*FmToGeV);
+ Num_Cos_os->Fill(kt, qinv12, BEE);
+ Num_CosFSI_os->Fill(kt, qinv12, WeightWF);
+ NumSq_CosFSI_os->Fill(kt, qinv12, pow(WeightWF,2));
+ Den_os->Fill(kt, qinv12);
+ //
+ K2_os->Fill(kt, qinv12, WeightWF);
+ PlaneWF_os->Fill(kt, qinv12, WeightPlaneWF);
+ }
+
+
+ if(!ThreeParticle) continue;
+
+
+ TLorentzVector LP1(P[EN][index1][1], P[EN][index1][2], P[EN][index1][3], P[EN][index1][0]);
+ TLorentzVector LX1(X[EN][index1][1], X[EN][index1][2], X[EN][index1][3], X[EN][index1][0]);
+ TLorentzVector LP2(P[0][index2][1], P[0][index2][2], P[0][index2][3], P[0][index2][0]);
+ TLorentzVector LX2(X[0][index2][1], X[0][index2][2], X[0][index2][3], X[0][index2][0]);
+
+
+ if(!LambdaDirect) {if(!Primary[EN][index1] || !Primary[0][index2]) continue;}
+ if(qinv12 > qCutOff) continue;
+ //
+ int EN2=0;
+ int start3=index2+1;
+ if(EN==1) {EN2=2; start3=0;}
+
+ for(int index3=start3; index3<Nparticles[EN2]; index3++){
+
+ if(!LambdaDirect && !Primary[EN2][index3]) continue;
+
+ double qinv23 = sqrt(pow(P[0][index2][1]-P[EN2][index3][1],2)+pow(P[0][index2][2]-P[EN2][index3][2],2)+pow(P[0][index2][3]-P[EN2][index3][3],2)-pow(P[0][index2][0]-P[EN2][index3][0],2));
+ double qinv31 = sqrt(pow(P[EN2][index3][1]-P[EN][index1][1],2)+pow(P[EN2][index3][2]-P[EN][index1][2],2)+pow(P[EN2][index3][3]-P[EN][index1][3],2)-pow(P[EN2][index3][0]-P[EN][index1][0],2));
+ if(qinv23 > qCutOff || qinv31 > qCutOff) continue;
+ if(qinv23 < 0.001 || qinv31 < 0.001) continue;
+
+ if(PID[EN][index1] == PID[0][index2] && PID[EN][index1] == PID[EN2][index3]){
+ DenLambda3_ss->Fill(qinv12, qinv23, qinv31);
+ }else{
+ DenLambda3_os->Fill(qinv12, qinv23, qinv31);
+ }
+
+
+ bool MixedCharge=kFALSE;
+ if(PID[EN][index1] != PID[0][index2] || PID[EN][index1] != PID[EN2][index3] || PID[0][index2] != PID[EN2][index3]) MixedCharge=kTRUE;
+
+ if(EN2==2){
+ if(MixedCharge==kFALSE){
+ r3numME->Fill(qinv12, qinv23, qinv31);
+ r3den1ME->Fill(qinv12, qinv23, qinv31);
+ r3den2ME->Fill(qinv12, qinv23, qinv31);
+ r3den3ME->Fill(qinv12, qinv23, qinv31);
+ }
+ continue;
+ }
+
+
+
+ // 12 boost
+ BoostPRF(P[EN][index1], P[0][index2], P[EN2][index3], &P3_PRF);
+ BoostPRF(P[EN][index1], P[0][index2], P[0][index2], &P2_PRF);
+ BoostPRF(P[EN][index1], P[0][index2], P[EN][index1], &P1_PRF);
+ X1_PRF[0]=X[EN][index1][1]; X1_PRF[1]=X[EN][index1][2]; X1_PRF[2]=X[EN][index1][3];// default CMS
+ X2_PRF[0]=X[0][index2][1]; X2_PRF[1]=X[0][index2][2]; X2_PRF[2]=X[0][index2][3];// default CMS
+ X3_PRF[0]=X[EN2][index3][1]; X3_PRF[1]=X[EN2][index3][2]; X3_PRF[2]=X[EN2][index3][3];// default CMS
+
+ //
+ if(BoostPairsIntoPRF){
+ BoostPRF(P[EN][index1], P[0][index2], X[EN2][index3], &X3_PRF);
+ BoostPRF(P[EN][index1], P[0][index2], X[0][index2], &X2_PRF);
+ BoostPRF(P[EN][index1], P[0][index2], X[EN][index1], &X1_PRF);
+ }
+
+
+ rstar12_PRF = (X1_PRF-X2_PRF).Mag();
+ if(rstar12_PRF > RstarMax) continue;
+ if(rstar12_PRF < RstarMin) continue;
+ double phase_p12_x12 = (P1_PRF-P2_PRF)*(X1_PRF-X2_PRF)/2.;
+ double phase_p12_x32 = (P1_PRF-P2_PRF)*(X3_PRF-X2_PRF)/2.;
+ double phase_p12_x13 = (P1_PRF-P2_PRF)*(X1_PRF-X3_PRF)/2.;
+ complex<double> z_p12_x12(0.0, (P1_PRF-P2_PRF)*(X1_PRF-X2_PRF)/2. + qinv12/2. * (X1_PRF-X2_PRF).Mag());
+ complex<double> z_p12_x21(0.0, (P1_PRF-P2_PRF)*(X2_PRF-X1_PRF)/2. + qinv12/2. * (X2_PRF-X1_PRF).Mag());
+ complex<double> z_p12_x32(0.0, (P1_PRF-P2_PRF)*(X3_PRF-X2_PRF)/2. + qinv12/2. * (X3_PRF-X2_PRF).Mag());
+ complex<double> z_p12_x13(0.0, (P1_PRF-P2_PRF)*(X1_PRF-X3_PRF)/2. + qinv12/2. * (X1_PRF-X3_PRF).Mag());
+ complex<double> z_p12_x23(0.0, (P1_PRF-P2_PRF)*(X2_PRF-X3_PRF)/2. + qinv12/2. * (X2_PRF-X3_PRF).Mag());
+ complex<double> z_p12_x31(0.0, (P1_PRF-P2_PRF)*(X3_PRF-X1_PRF)/2. + qinv12/2. * (X3_PRF-X1_PRF).Mag());
+
+
+ // 23 boost
+ BoostPRF(P[0][index2], P[EN2][index3], P[EN2][index3], &P3_PRF);
+ BoostPRF(P[0][index2], P[EN2][index3], P[0][index2], &P2_PRF);
+ BoostPRF(P[0][index2], P[EN2][index3], P[EN][index1], &P1_PRF);
+ if(BoostPairsIntoPRF){
+ BoostPRF(P[0][index2], P[EN2][index3], X[EN2][index3], &X3_PRF);
+ BoostPRF(P[0][index2], P[EN2][index3], X[0][index2], &X2_PRF);
+ BoostPRF(P[0][index2], P[EN2][index3], X[EN][index1], &X1_PRF);
+ }
+ double rstar23_PRF = (X2_PRF-X3_PRF).Mag();
+ if(rstar23_PRF > RstarMax) continue;
+ if(rstar23_PRF < RstarMin) continue;
+ double phase_p23_x13 = (P2_PRF-P3_PRF)*(X1_PRF-X3_PRF)/2.;
+ double phase_p23_x23 = (P2_PRF-P3_PRF)*(X2_PRF-X3_PRF)/2.;
+ double phase_p23_x21 = (P2_PRF-P3_PRF)*(X2_PRF-X1_PRF)/2.;
+ double FVP23 = (P2_PRF-P3_PRF).Mag();
+ FVP23 -= pow(P[0][index2][0]-P[EN2][index3][0],2);
+ double tstar23_PRF = sqrt(fabs(pow(rstar23_PRF,2)-FVP23));
+ //double rstar23_CMS = sqrt(pow(X[0][index2][1]-X[EN2][index3][1],2)+pow(X[0][index2][2]-X[EN2][index3][2],2)+pow(X[0][index2][3]-X[EN2][index3][3],2));
+ complex<double> z_p23_x23(0.0, (P2_PRF-P3_PRF)*(X2_PRF-X3_PRF)/2. + qinv23/2. * (X2_PRF-X3_PRF).Mag());
+ complex<double> z_p23_x13(0.0, (P2_PRF-P3_PRF)*(X1_PRF-X3_PRF)/2. + qinv23/2. * (X1_PRF-X3_PRF).Mag());
+ complex<double> z_p23_x21(0.0, (P2_PRF-P3_PRF)*(X2_PRF-X1_PRF)/2. + qinv23/2. * (X2_PRF-X1_PRF).Mag());
+ complex<double> z_p23_x32(0.0, (P2_PRF-P3_PRF)*(X3_PRF-X2_PRF)/2. + qinv23/2. * (X3_PRF-X2_PRF).Mag());
+ complex<double> z_p23_x31(0.0, (P2_PRF-P3_PRF)*(X3_PRF-X1_PRF)/2. + qinv23/2. * (X3_PRF-X1_PRF).Mag());
+ complex<double> z_p23_x12(0.0, (P2_PRF-P3_PRF)*(X1_PRF-X2_PRF)/2. + qinv23/2. * (X1_PRF-X2_PRF).Mag());
+
+
+ // 31 boost
+ BoostPRF(P[EN2][index3], P[EN][index1], P[EN2][index3], &P3_PRF);
+ BoostPRF(P[EN2][index3], P[EN][index1], P[0][index2], &P2_PRF);
+ BoostPRF(P[EN2][index3], P[EN][index1], P[EN][index1], &P1_PRF);
+ if(BoostPairsIntoPRF){
+ BoostPRF(P[EN2][index3], P[EN][index1], X[EN2][index3], &X3_PRF);
+ BoostPRF(P[EN2][index3], P[EN][index1], X[0][index2], &X2_PRF);
+ BoostPRF(P[EN2][index3], P[EN][index1], X[EN][index1], &X1_PRF);
+ }
+ double rstar31_PRF = (X3_PRF-X1_PRF).Mag();
+ if(rstar31_PRF > RstarMax) continue;
+ if(rstar31_PRF < RstarMin) continue;
+ double phase_p31_x31 = (P3_PRF-P1_PRF)*(X3_PRF-X1_PRF)/2.;
+ double phase_p31_x32 = (P3_PRF-P1_PRF)*(X3_PRF-X2_PRF)/2.;
+ double phase_p31_x21 = (P3_PRF-P1_PRF)*(X2_PRF-X1_PRF)/2.;
+ double FVP31 = (P3_PRF-P1_PRF).Mag();
+ FVP31 -= pow(P[EN2][index3][0]-P[EN][index1][0],2);
+ double tstar31_PRF = sqrt(fabs(pow(rstar31_PRF,2)-FVP31));
+ //double rstar31_CMS = sqrt(pow(X[EN2][index3][1]-X[EN][index1][1],2)+pow(X[EN2][index3][2]-X[EN][index1][2],2)+pow(X[EN2][index3][3]-X[EN][index1][3],2));
+ complex<double> z_p31_x31(0.0, (P3_PRF-P1_PRF)*(X3_PRF-X1_PRF)/2. + qinv31/2. * (X3_PRF-X1_PRF).Mag());
+ complex<double> z_p31_x32(0.0, (P3_PRF-P1_PRF)*(X3_PRF-X2_PRF)/2. + qinv31/2. * (X3_PRF-X2_PRF).Mag());
+ complex<double> z_p31_x13(0.0, (P3_PRF-P1_PRF)*(X1_PRF-X3_PRF)/2. + qinv31/2. * (X1_PRF-X3_PRF).Mag());
+ complex<double> z_p31_x21(0.0, (P3_PRF-P1_PRF)*(X2_PRF-X1_PRF)/2. + qinv31/2. * (X2_PRF-X1_PRF).Mag());
+ complex<double> z_p31_x12(0.0, (P3_PRF-P1_PRF)*(X1_PRF-X2_PRF)/2. + qinv31/2. * (X1_PRF-X2_PRF).Mag());
+ complex<double> z_p31_x23(0.0, (P3_PRF-P1_PRF)*(X2_PRF-X3_PRF)/2. + qinv31/2. * (X2_PRF-X3_PRF).Mag());
+
+
+
+ ////////////////////////////////////////////
+ ////////////////////////////////////////////
+ double Q3 = sqrt(pow(qinv12,2) + pow(qinv23,2) + pow(qinv31,2));
+
+ if(!LambdaDirect){
+ int CP12=+1, CP23=+1, CP31=+1;
+ if(PID[EN][index1] != PID[0][index2]) CP12 = -1;
+ if(PID[0][index2] != PID[EN2][index3]) CP23 = -1;
+ if(PID[EN2][index3] != PID[EN][index1]) CP31 = -1;
+
+ double Ac_12 = Gamov(CP12/(qinv12/2.*BohrR));
+ double Ac_23 = Gamov(CP23/(qinv23/2.*BohrR));
+ double Ac_31 = Gamov(CP31/(qinv31/2.*BohrR));
+
+ //
+ complex<double> a_12(0.0, -CP12/(qinv12/2.*BohrR));
+ complex<double> a_23(0.0, -CP23/(qinv23/2.*BohrR));
+ complex<double> a_31(0.0, -CP31/(qinv31/2.*BohrR));
+ //
+ //
+ //
+ TLorentzVector LP3(P[EN2][index3][1], P[EN2][index3][2], P[EN2][index3][3], P[EN2][index3][0]);
+ TLorentzVector LX3(X[EN2][index3][1], X[EN2][index3][2], X[EN2][index3][3], X[EN2][index3][0]);
+
+ double phasePW = ( (LP1-LP2)*(LX1-LX2) + (LP1-LP3)*(LX1-LX3) + (LP2-LP3)*(LX2-LX3) )/3.;// base
+ complex<double> planewave1(cos(phasePW), sin(phasePW));
+ //
+ phasePW = ( (LP1-LP2)*(LX2-LX1) + (LP1-LP3)*(LX2-LX3) + (LP2-LP3)*(LX1-LX3) )/3.;// 12 swap
+ complex<double> planewave2(cos(phasePW), sin(phasePW));
+ //
+ phasePW = ( (LP1-LP2)*(LX3-LX2) + (LP1-LP3)*(LX3-LX1) + (LP2-LP3)*(LX2-LX1) )/3.;// 13 swap
+ complex<double> planewave3(cos(phasePW), sin(phasePW));
+ //
+ phasePW = ( (LP1-LP2)*(LX1-LX3) + (LP1-LP3)*(LX1-LX2) + (LP2-LP3)*(LX3-LX2) )/3.;// 23 swap
+ complex<double> planewave4(cos(phasePW), sin(phasePW));
+ //
+ phasePW = ( (LP1-LP2)*(LX2-LX3) + (LP1-LP3)*(LX2-LX1) + (LP2-LP3)*(LX3-LX1) )/3.;// 123 swap (1st type)
+ complex<double> planewave5(cos(phasePW), sin(phasePW));
+ //
+ phasePW = ( (LP1-LP2)*(LX3-LX1) + (LP1-LP3)*(LX3-LX2) + (LP2-LP3)*(LX1-LX2) )/3.;// 123 swap (2nd type)
+ complex<double> planewave6(cos(phasePW), sin(phasePW));
+ //
+ //
+ //
+ complex<double> CHG_p12_x12= conf_hyp(a_12, binput, z_p12_x12);
+ complex<double> CHG_p12_x21= conf_hyp(a_12, binput, z_p12_x21);
+ complex<double> CHG_p12_x32= conf_hyp(a_12, binput, z_p12_x32);
+ complex<double> CHG_p12_x13= conf_hyp(a_12, binput, z_p12_x13);
+ complex<double> CHG_p12_x23= conf_hyp(a_12, binput, z_p12_x23);
+ complex<double> CHG_p12_x31= conf_hyp(a_12, binput, z_p12_x31);
+ //
+ complex<double> CHG_p31_x31= conf_hyp(a_31, binput, z_p31_x31);
+ complex<double> CHG_p31_x32= conf_hyp(a_31, binput, z_p31_x32);
+ complex<double> CHG_p31_x13= conf_hyp(a_31, binput, z_p31_x13);
+ complex<double> CHG_p31_x21= conf_hyp(a_31, binput, z_p31_x21);
+ complex<double> CHG_p31_x12= conf_hyp(a_31, binput, z_p31_x12);
+ complex<double> CHG_p31_x23= conf_hyp(a_31, binput, z_p31_x23);
+ //
+ complex<double> CHG_p23_x23= conf_hyp(a_23, binput, z_p23_x23);
+ complex<double> CHG_p23_x13= conf_hyp(a_23, binput, z_p23_x13);
+ complex<double> CHG_p23_x21= conf_hyp(a_23, binput, z_p23_x21);
+ complex<double> CHG_p23_x32= conf_hyp(a_23, binput, z_p23_x32);
+ complex<double> CHG_p23_x31= conf_hyp(a_23, binput, z_p23_x31);
+ complex<double> CHG_p23_x12= conf_hyp(a_23, binput, z_p23_x12);
+
+ ThreePhase_1 = X1_PRF*(P1_PRF-P2_PRF) + X2_PRF*(P2_PRF-P3_PRF) + X3_PRF*(P3_PRF-P1_PRF);
+ ThreePhase_2 = X1_PRF*(P1_PRF-P2_PRF) + X3_PRF*(P2_PRF-P3_PRF) + X2_PRF*(P3_PRF-P1_PRF);
+ ThreePhase_3 = X2_PRF*(P1_PRF-P2_PRF) + X1_PRF*(P2_PRF-P3_PRF) + X3_PRF*(P3_PRF-P1_PRF);
+ ThreePhase_4 = X2_PRF*(P1_PRF-P2_PRF) + X3_PRF*(P2_PRF-P3_PRF) + X1_PRF*(P3_PRF-P1_PRF);
+ ThreePhase_5 = X3_PRF*(P1_PRF-P2_PRF) + X1_PRF*(P2_PRF-P3_PRF) + X2_PRF*(P3_PRF-P1_PRF);
+ ThreePhase_6 = X3_PRF*(P1_PRF-P2_PRF) + X2_PRF*(P2_PRF-P3_PRF) + X1_PRF*(P3_PRF-P1_PRF);
+ TwoPhase_12 = (P1_PRF-P2_PRF)*(X1_PRF-X2_PRF);
+ TwoPhase_23 = (P2_PRF-P3_PRF)*(X2_PRF-X3_PRF);
+ TwoPhase_31 = (P3_PRF-P1_PRF)*(X3_PRF-X1_PRF);
+ //
+ //
+
+ // Strong FSI
+ if(MixedCharge){
+ if(CP12==+1){// 1 and 2 are same-charge
+ kstar_1 = qinv31/2.; kstar_2 = qinv23/2.;
+ rstar_1 = rstar31_PRF; rstar_2 = rstar23_PRF;
+ phase_1 = phase_p31_x31; phase_2 = phase_p31_x32; phase_3 = phase_p23_x13; phase_4 = phase_p23_x23;
+ CHG_1 = CHG_p31_x31;
+ CHG_2 = CHG_p31_x32;
+ CHG_3 = CHG_p23_x13;
+ CHG_4 = CHG_p23_x23;
+ qinvMixedC_ss=qinv12; qinvMixedC_os1=qinv31; qinvMixedC_os2=qinv23;
+ }else if(CP31==+1){// 1 and 3 are same-charge
+ kstar_1 = qinv12/2.; kstar_2 = qinv23/2.;
+ rstar_1 = rstar12_PRF; rstar_2 = rstar23_PRF;
+ phase_1 = phase_p12_x12; phase_2 = phase_p12_x32; phase_3 = phase_p23_x21; phase_4 = phase_p23_x23;
+ CHG_1 = CHG_p12_x12;
+ CHG_2 = CHG_p12_x32;
+ CHG_3 = CHG_p23_x21;
+ CHG_4 = CHG_p23_x23;
+ qinvMixedC_ss=qinv31; qinvMixedC_os1=qinv12; qinvMixedC_os2=qinv23;
+ }else {// 2 and 3 are same-charge
+ kstar_1 = qinv12/2.; kstar_2 = qinv31/2.;
+ rstar_1 = rstar12_PRF; rstar_2 = rstar31_PRF;
+ phase_1 = phase_p12_x12; phase_2 = phase_p12_x13; phase_3 = phase_p31_x21; phase_4 = phase_p31_x31;
+ CHG_1 = CHG_p12_x12;
+ CHG_2 = CHG_p12_x13;
+ CHG_3 = CHG_p31_x21;
+ CHG_4 = CHG_p31_x31;
+ qinvMixedC_ss=qinv23; qinvMixedC_os1=qinv31; qinvMixedC_os2=qinv31;
+ }
+
+
+ for(int piece=0; piece<4; piece++){
+ if(piece==0) {// P31X31
+ kstar=kstar_1;
+ rstar=rstar_1;
+ phase=phase_1;
+ CHG=CHG_1;
+ }else if(piece==1) {// P31X32
+ kstar=kstar_1;
+ rstar=rstar_2;
+ phase=phase_2;
+ CHG=CHG_2;
+ }else if(piece==2) {// P23X13
+ kstar=kstar_2;
+ rstar=rstar_1;
+ phase=phase_3;
+ CHG=CHG_3;
+ }else {// P23X23
+ kstar=kstar_2;
+ rstar=rstar_2;
+ phase=phase_4;
+ CHG=CHG_4;
+ }
+
+ rho = kstar*rstar;
+ eta = -1./(kstar*BohrR);
+ double kb = sqrt(pow(massPOI,2)+pow(kstar,2) - pow(masspi0,2));
+ double G2 = KinverseFunction(pow(kstar,2),1);
+ double G1 = KinverseFunction(pow(kstar,2),0);
+ double RK11 = 2/3.*G1 + 1/3.*G2;
+ double RK22 = 2/3.*G2 + 1/3.*G1;
+ double RK12 = -sqrt(2/9.)*(G1-G2);
+ complex<double> Chi(h_histogram->GetBinContent(h_histogram->GetXaxis()->FindBin(eta)), Gamov(eta)/(2*eta));
+ complex<double> C3 = RK11 - 2.0*Chi/BohrR;
+ complex<double> C5(RK22, -kb);
+ complex<double> C10 = C3*C5 - pow(RK12,2);
+ complex<double> fc_aa = C5/C10;
+ complex<double> fc_ab = -RK12/C10;
+ complex<double> G = P_histogram->GetBinContent(P_histogram->GetXaxis()->FindBin(rho), P_histogram->GetYaxis()->FindBin(eta));
+ G += 2.0*eta*rho*( log(fabs(2.0*eta*rho)) + 2.0*EulerC - 1.0 + Chi )*B_histogram->GetBinContent(B_histogram->GetXaxis()->FindBin(rho), B_histogram->GetYaxis()->FindBin(eta));
+ complex<double> PW(cos(phase), sin(phase));
+ AlphaWF[piece] = CHG + PW*fc_aa*G/rstar;
+ complex<double> spherewave(cos(rstar*kb), sin(rstar*kb));
+ BetaWF[piece] = PW*fc_ab*sqrt(masspi0/massPOI)*spherewave/rstar;
+
+ }// piece
+
+ }// Mixed-charge case with Strong FSI
+
+ complex<double> FullFSI_WF[2];// aa, bb
+ if(MixedCharge==kTRUE){
+ if(StrongFSI==kFALSE) {
+ FullFSI_WF[0] = planewave1*CHG_p12_x12*CHG_p31_x31*CHG_p23_x23;
+ if(CP12==+1){// 1 and 2 are same-charge
+ FullFSI_WF[0] += planewave2*CHG_p12_x21*CHG_p31_x32*CHG_p23_x13;// 12 symmetrization
+ }else if(CP31==+1){// 1 and 3 are same-charge
+ FullFSI_WF[0] += planewave3*CHG_p12_x32*CHG_p31_x13*CHG_p23_x21;// 13 symmetrization
+ }else{
+ FullFSI_WF[0] += planewave4*CHG_p12_x13*CHG_p31_x21*CHG_p23_x32;// 23 symmetrization
+ }
+ FullFSI_WF[1] = 0.;
+ }else {// strong FSI
+ if(CP12==+1){// 1 and 2 are same-charge
+ FullFSI_WF[0] = planewave1*CHG_p12_x12*AlphaWF[0]*AlphaWF[3];
+ FullFSI_WF[1] = planewave1*CHG_p12_x12*BetaWF[0]*BetaWF[3];
+ FullFSI_WF[0] += planewave2*CHG_p12_x21*AlphaWF[1]*AlphaWF[2];// 12 symmetrization
+ FullFSI_WF[1] += planewave2*CHG_p12_x21*BetaWF[1]*BetaWF[2];// 12 symmetrization
+ }else if(CP31==+1){// 1 and 3 are same-charge
+ FullFSI_WF[0] = planewave1*CHG_p31_x31*AlphaWF[0]*AlphaWF[3];
+ FullFSI_WF[1] = planewave1*CHG_p31_x31*BetaWF[0]*BetaWF[3];
+ FullFSI_WF[0] += planewave3*CHG_p31_x13*AlphaWF[1]*AlphaWF[2];// 13 symmetrization
+ FullFSI_WF[1] += planewave3*CHG_p31_x13*BetaWF[1]*BetaWF[2];// 13 symmetrization
+ }else{
+ FullFSI_WF[0] = planewave1*CHG_p23_x23*AlphaWF[0]*AlphaWF[3];
+ FullFSI_WF[1] = planewave1*CHG_p23_x23*BetaWF[0]*BetaWF[3];
+ FullFSI_WF[0] += planewave4*CHG_p23_x32*AlphaWF[1]*AlphaWF[2];// 13 symmetrization
+ FullFSI_WF[1] += planewave4*CHG_p23_x32*BetaWF[1]*BetaWF[2];// 13 symmetrization
+ }
+ }
+ FullFSI_WF[0] /= sqrt(2.); FullFSI_WF[1] /= sqrt(2.);
+ }else {
+ FullFSI_WF[0] = planewave1*CHG_p12_x12*CHG_p31_x31*CHG_p23_x23;// base
+ FullFSI_WF[0] += planewave2*CHG_p12_x21*CHG_p31_x32*CHG_p23_x13;// 12 symmetrization
+ FullFSI_WF[0] += planewave3*CHG_p12_x32*CHG_p31_x13*CHG_p23_x21;// 13
+ FullFSI_WF[0] += planewave4*CHG_p12_x13*CHG_p31_x21*CHG_p23_x32;// 23
+ FullFSI_WF[0] += planewave5*CHG_p12_x23*CHG_p31_x12*CHG_p23_x31;// 123
+ FullFSI_WF[0] += planewave6*CHG_p12_x31*CHG_p31_x23*CHG_p23_x12;// 123
+ FullFSI_WF[1] = 0.;
+ FullFSI_WF[0] /= sqrt(6.); FullFSI_WF[1] /= sqrt(6.);
+ }
+
+ complex<double> FullPlaneWF = planewave1;// base
+ if(MixedCharge && CP12==+1) FullPlaneWF += planewave2;
+ if(MixedCharge && CP31==+1) FullPlaneWF += planewave3;
+ if(MixedCharge && CP23==+1) FullPlaneWF += planewave4;
+ if(!MixedCharge) {
+ FullPlaneWF += planewave2 + planewave3 + planewave4 + planewave5 + planewave6;
+ FullPlaneWF /= sqrt(6.);
+ }else FullPlaneWF /= sqrt(2.);
+
+ WeightWF_aa = pow(abs(FullFSI_WF[0]),2);
+ WeightWF_bb = pow(abs(FullFSI_WF[1]),2);
+ WeightWF_aa *= Ac_12*Ac_23*Ac_31;
+ WeightWF_bb *= Ac_12*Ac_23*Ac_31;
+ WeightPlaneWF = pow(abs(FullPlaneWF),2);
+
+ }// !LambdaDirect
+ ////////////////////////////////////////////
+ ////////////////////////////////////////////
+ //
+
+
+ if(MixedCharge){
+ K3os->Fill(Q3, WeightWF_aa + WeightWF_bb);
+ PlaneWF3os->Fill(Q3, WeightPlaneWF);
+ K3os_3D->Fill(qinvMixedC_ss, qinvMixedC_os1, qinvMixedC_os2, WeightWF_aa + WeightWF_bb);
+ PlaneWF3os_3D->Fill(qinvMixedC_ss, qinvMixedC_os1, qinvMixedC_os2, WeightPlaneWF);
+ NumLambda3_os->Fill(qinv12, qinv23, qinv31);
+ }else{
+ K3ss->Fill(Q3, WeightWF_aa + WeightWF_bb);
+ PlaneWF3ss->Fill(Q3, WeightPlaneWF);
+ K3ss_3D->Fill(qinv12, qinv23, qinv31, WeightWF_aa + WeightWF_bb);
+ PlaneWF3ss_3D->Fill(qinv12, qinv23, qinv31, WeightPlaneWF);
+ NumLambda3_ss->Fill(qinv12, qinv23, qinv31);
+ rstar3_dist_ss->Fill(rstar12_PRF, rstar23_PRF, rstar31_PRF);
+ //rstar3_dist_ss->Fill(rstar12_CMS, rstar23_CMS, rstar31_CMS);
+ tstar3_dist_ss->Fill(tstar12_PRF, tstar23_PRF, tstar31_PRF);
+ //
+ double cosThreePhase = 1/6. * (cos(ThreePhase_1)+cos(ThreePhase_2)+cos(ThreePhase_3)+cos(ThreePhase_4)+cos(ThreePhase_5)+cos(ThreePhase_6));
+ r3num->Fill(qinv12, qinv23, qinv31, 2*cosThreePhase);
+ r3den1->Fill(qinv12, qinv23, qinv31, cos(TwoPhase_12));
+ r3den2->Fill(qinv12, qinv23, qinv31, cos(TwoPhase_23));
+ r3den3->Fill(qinv12, qinv23, qinv31, cos(TwoPhase_31));
+ r3numSq->Fill(qinv12, qinv23, qinv31, pow(2*cosThreePhase,2));
+ r3den1Sq->Fill(qinv12, qinv23, qinv31, pow(cos(TwoPhase_12),2));
+ r3den2Sq->Fill(qinv12, qinv23, qinv31, pow(cos(TwoPhase_23),2));
+ r3den3Sq->Fill(qinv12, qinv23, qinv31, pow(cos(TwoPhase_31),2));
+ r3numEn->Fill(qinv12, qinv23, qinv31);
+ r3den1En->Fill(qinv12, qinv23, qinv31);
+ r3den2En->Fill(qinv12, qinv23, qinv31);
+ r3den3En->Fill(qinv12, qinv23, qinv31);
+ //
+ Num_Cos12->Fill(qinv12,cos(TwoPhase_12));
+ Num_Cos23->Fill(qinv23,cos(TwoPhase_23));
+ Num_Cos31->Fill(qinv31,cos(TwoPhase_31));
+ Den_Cos12->Fill(qinv12);
+ Den_Cos23->Fill(qinv23);
+ Den_Cos31->Fill(qinv31);
+ }
+
+
+ }// index3
+
+ }// index2
+
+ }// index1
+
+ }// EN buffer
+
+ }// Nevents
+
+ infile->Close();
+ }// Nfiles
+
+ outfile->cd();
+ //
+ r_dist->Write();
+ x_dist->Write();
+ Pt_dist->Write();
+ P_dist->Write();
+ P_dist_lowq->Write();
+ AvgMult->Write();
+
+ rstarPRF_dist_ss->Write();
+ tstar_dist_ss->Write();
+ rstarPRF_dist_os->Write();
+ tstar_dist_os->Write();
+ NumLambda1->Write();
+ DenLambda1->Write();
+ NumLambda2_ss->Write();
+ NumLambda2_os->Write();
+ DenLambda2_ss->Write();
+ DenLambda2_os->Write();
+ NumLambda3_ss->Write();
+ DenLambda3_ss->Write();
+ NumLambda3_os->Write();
+ DenLambda3_os->Write();
+ Den_ss->Write();
+ Den_os->Write();
+ DenEM_ss->Write();
+ DenEM_os->Write();
+ Num_Cos_ss->Write();
+ Num_Cos_os->Write();
+ Num_CosFSI_ss->Write();
+ Num_CosFSI_os->Write();
+ NumSq_CosFSI_ss->Write();
+ NumSq_CosFSI_os->Write();
+ LargeRpairs_ss->Write();
+ LargeRpairs_os->Write();
+ //
+ Num_PrimCosFSI_ss->Write();
+ Num_PrimCosFSI_os->Write();
+ //
+ rstar3_dist_ss->Write();
+ tstar3_dist_ss->Write();
+ //
+ K2_ss->Write();
+ PlaneWF_ss->Write();
+ K2_os->Write();
+ PlaneWF_os->Write();
+ //
+ K3os->Write();
+ PlaneWF3os->Write();
+ K3os_3D->Write();
+ PlaneWF3os_3D->Write();
+ K3ss->Write();
+ PlaneWF3ss->Write();
+ K3ss_3D->Write();
+ PlaneWF3ss_3D->Write();
+ //
+ r3num->Write();
+ r3den1->Write();
+ r3den2->Write();
+ r3den3->Write();
+ r3numSq->Write();
+ r3den1Sq->Write();
+ r3den2Sq->Write();
+ r3den3Sq->Write();
+ r3numEn->Write();
+ r3den1En->Write();
+ r3den2En->Write();
+ r3den3En->Write();
+ r3numME->Write();
+ r3den1ME->Write();
+ r3den2ME->Write();
+ r3den3ME->Write();
+ //
+ Num_Cos12->Write();
+ Num_Cos23->Write();
+ Num_Cos31->Write();
+ Den_Cos12->Write();
+ Den_Cos23->Write();
+ Den_Cos31->Write();
+ //
+ test->Write();
+ //
+ outfile->Close();
+
+ seconds = time(0) - seconds;
+ cout<<"Minutes = "<<seconds/60.<<endl;
+
+
+}
+float Gamov(float eta){
+ return (2*PI*eta/(exp(2*PI*eta)-1));
+}
+float fact(float n){
+ return (n < 1.00001) ? 1 : fact(n - 1) * n;
+}
+void Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
+{
+ Int_t j, k;
+ Int_t a = i2 - i1;
+ for (Int_t i = i1; i < i2+1; i++) {
+ j = (Int_t) (gRandom->Rndm() * a);
+ k = iarr[j];
+ iarr[j] = iarr[i];
+ iarr[i] = k;
+ }
+}
+void BoostPRF(float P1[4], float P2[4], float V[4], TVector3 *T){
+ // P1 is particle 1 momentum in reaction CMS (0=E,1=px,2=py,3=pz)
+ // P2 is particle 2 momentum in reaction CMS (0=E,1=px,2=py,3=pz)
+ // V is the vector to be transformed (V is in CMS)
+ // T is the PRF vector
+
+ float MT=sqrt(pow(P1[0]+P2[0],2)-pow(P1[3]+P2[3],2));
+ float PT=sqrt(pow(P1[1]+P2[1],2)+pow(P1[2]+P2[2],2));
+ float MINV=sqrt(pow(P1[0]+P2[0],2) - pow(P1[1]+P2[1],2) - pow(P1[2]+P2[2],2) - pow(P1[3]+P2[3],2));
+ // LCMS
+ T->SetX( ((P1[1]+P2[1])*V[1] + (P1[2]+P2[2])*V[2])/PT );
+ T->SetY( ((P1[1]+P2[1])*V[2] - (P1[2]+P2[2])*V[1])/PT );
+ T->SetZ( ((P1[0]+P2[0])*V[3] - (P1[3]+P2[3])*V[0])/MT );
+ // PRF
+ T->SetX( T->X()*MINV/MT - PT/(MT*MINV)*((P1[0]+P2[0])*V[0] - (P1[1]+P2[1])*V[1] - (P1[2]+P2[2])*V[2] - (P1[3]+P2[3])*V[3]) );
+}
+// h function from Lednicky, phys. of part. and nuclei 40:307 (2009)
+double hFunction(double eta){
+ double DS=1.0;
+ double N=0;
+ double S=0;
+ double returnValue=0;
+ if(eta < 3.0){
+ while(DS > 1.0e-13){
+ N++;
+ DS = 1/(N*(pow(N/eta,2)+1.0));
+ S += DS;
+ }
+ returnValue = (S - EulerC - log(fabs(eta)));
+ }else{// small kstar, high eta
+ double etaSquared=pow(eta,2);
+ double etaPowered=etaSquared;
+ for(int ii=0; ii<9; ii++){// 9 was maximum value for Lednicky's code
+ returnValue += 1/etaPowered * hcoef[ii];
+ etaPowered *= etaSquared;
+ }
+ }
+ return returnValue;
+}
+double KinverseFunction(double kstarSq, int J){
+ double ESq = kstarSq + pow(massPOI,2);
+ double E = sqrt(ESq);
+ double xSq = kstarSq/pow(massPOI,2);
+ double Kinverse = E*(4*ESq-fzero[J][4])/(4*pow(massPOI,2)-fzero[J][4]);
+ Kinverse /= fzero[J][0] + fzero[J][1]*xSq + fzero[J][2]*pow(xSq,2) + fzero[J][3]*pow(xSq,3);
+ return Kinverse;
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff