PWGCF/FEMTOSCOPY/macros/Therm.C

   1 #include <math.h>
   2 #include <time.h>
   3 #include <stdio.h>
   4 #include <stdlib.h>
   5 #include <Riostream.h>
   6 #include <complex>
   7
   8 #include "TObject.h"
   9 #include "TTree.h"
  10 #include "TBranch.h"
  11 #include "TLeaf.h"
  12 #include "TVector2.h"
  13 #include "TVector3.h"
  14 #include "TFile.h"
  15 #include "TString.h"
  16 #include "TF1.h"
  17 #include "TH1.h"
  18 #include "TH2.h"
  19 #include "TH3.h"
  20 #include "TProfile.h"
  21 #include "TProfile2D.h"
  22 #include "TMath.h"
  23 #include "TText.h"
  24 #include "TRandom3.h"
  25 #include "TArray.h"
  26 #include "TLegend.h"
  27 #include "TStyle.h"
  28 #include "TMinuit.h"
  29 #include "TChain.h"
  30 #include "TLorentzVector.h"
  31
  32 #define PI 3.1415926
  33 #define BohrR 1963.6885 // Mate's value(1963.6885) ~ 387.5/0.197327(1963.7455)
  34 #define FmToGeV 0.197327 // 0.197327
  35 #define fzero_aa 0.942597 // 0.186fm/FmToGeV (scattering length of pi+pi-) = 0.942597
  36 #define fzero_ab -0.89192 // -0.176fm/FmtoGeV (scattering length (pi+pi-)-->(pi0pi0) = -0.89192
  37 #define fzero_bb 0.3119 // fzero_aa + 1/sqrt(2)*fzero_ab = 0.0615/FmToGeV = 0.3119
  38 #define dzero -50.6773 // -10/FmToGeV = -50.6773 (effective range)
  39 #define EulerC 0.5772156649 // Euler constant
  40 #define masspi0 0.1349766 // pi0 mass
  41 #define masspiC 0.1395702 // pi+ mass
  42 #define massKch 0.493677 // K+- mass
  43 #define massK0 0.497614 // K0 mass
  44
  45 #define SF 1.0 // Scale Factor for spatial coordinates
  46 #define SF_Mom 1.0 // Scale Factor to reduce momenta
  47 #define RstarMax 405.4/SF // 405.4 (80 fm / FmToGeV), tried 253.4 (50 fm / FmToGeV) as a variation
  48 #define RstarMin 0.507/SF // 0.507 (0.1 fm / FmToGeV)
  49 using namespace std;
  50
  51 // g and p used in Gamma function
  52 const int g_gamma = 7;
  53 double p_gamma[9] = {0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7};
  54 // h coefficients for very small kstar (eta > 0.2)
  55 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};
  56 // Scattering length (in GeV) coefficients in isospin representation from Lednicky's code from Colangelo (2001)
  57 double fzero[2][6]={{.220, .268,-.0139,-.00139, 36.77, 15*0}, {-.0444,-.0857,-.00221,-.000129,-21.62, 15*0}};
  58
  59 double massPOI;
  60 int PIDPOI;
  61 double qCutOff;
  62
  63 float Gamov(float);
  64 float fact(float);
  65 double hFunction(double);
  66 double KinverseFunction(double, int);
  67 void Shuffle(int*, int, int);
  68 void BoostPRF(float [4], float [4], float [4], TVector3*);
  69
  70
  71 complex<double> gamma(complex<double> z)
  72 {// Lanczos approximation
  73
  74   if(floor(abs(z)) == z)
  75     return fact(abs(z)-1.0);
  76
  77   if(real(z)<(double)(0.5)) return PI/(sin(PI*z)*gamma((1.0-z)));
  78   z -= 1.0;
  79   complex<double> x = p_gamma[0];
  80   for(int i=1;i<g_gamma+2;i++) x = x+p_gamma[i]/(z+(double)i);
  81   complex<double> t = z+(double)g_gamma+0.5;
  82   return sqrt(2*PI)*pow(t,z+0.5)*exp(-t)*x;
  83 }
  84
  85 // confluent hypergeometric function
  86 complex<double> conf_hyp(complex<double> a, complex<double> b, complex<double> z)
  87 {// b is assumed to be 1.0 here
  88   complex<double> S(1.0, 0.0);
  89   complex<double> F=S;
  90   complex<double> alf1 = a - 1.0;
  91   double J=0;
  92   double errorCutOff=0.00001;// series truncation point (0.00001, 1e-6 for strict case)
  93   if(abs(z) > 40.) {// high k*r* case (was 20., now 40. Improved accuracy for qinv>60 MeV)
  94     double eta = -a.imag();
  95     double RKS = z.imag();
  96     double D1 = log(RKS)*eta;
  97     double D2 = pow(eta,2)/RKS;
  98     complex<double> arg(1.0, eta);
  99     complex<double> EIDC = gamma(arg);
 100     EIDC /= abs(EIDC);
 101     complex<double> ZZ1(cos(D1), sin(D1));
 102     ZZ1 /= EIDC;
 103     complex<double> value(1.0, -D2);
 104     value *= ZZ1;
 105     complex<double> value2(cos(RKS), sin(RKS));
 106     F = value - value2*eta/RKS/ZZ1;
 107     return F;
 108   }else {// low k*r* case
 109     while(fabs(S.real())+fabs(S.imag()) > errorCutOff){
 110       J++;
 111       complex<double> A = (alf1 + J)/(pow(J,2));
 112       S *= A*z;
 113       F += S;
 114     }
 115     return F;
 116   }
 117 }
 118
 119 void Therm(){
 120
 121   TRandom3 *RandomNumber = new TRandom3();
 122   RandomNumber->SetSeed(0);
 123
 124
 125   bool StrongFSI=kTRUE;
 126   bool LambdaDirect=kFALSE;
 127   bool RemoveXP=kFALSE;
 128   bool Gaussian=kFALSE;
 129   bool ThreeParticle=kFALSE;
 130   bool BoostPairsIntoPRF=kTRUE;
 131   massPOI = masspiC;
 132   PIDPOI = 211;// 211 for pi+-
 133   qCutOff = 0.1;// 0.1 or 0.4 for pp
 134   int qBins = int(qCutOff/0.002);
 135   //
 136   int Nfiles=1000;// 1000 for b2, 500 b3-b5, 1500 b7-b8, 3000 b9
 137   int bValue=2;
 138
 139   double RValue=8;// Gaussian radius
 140   double NsigmaGauss=2;// number of sigma to integrate over for a Gaussian source
 141   if(bValue==2) RValue=8;//8
 142   if(bValue==3) RValue=8;//8
 143   if(bValue==5) RValue=7;//7
 144   if(bValue==7) RValue=6;//6
 145   if(bValue==8) RValue=5;//5
 146   if(bValue==9) RValue=5;//5
 147   TF1 *SingleParticleGaussian = new TF1("SingleParticleGaussian","exp(-pow(x/(sqrt(2.)*[0]),2))",0,1000);
 148   SingleParticleGaussian->SetParameter(0,RValue);
 149
 150   TFile *outfile = new TFile("Therm_dist_temp.root","RECREATE");
 151   TH1D *r_dist=new TH1D("r_dist","",400,0,100);
 152   TH1D *x_dist=new TH1D("x_dist","",400,0,100);
 153   TH1D *Pt_dist=new TH1D("Pt_dist","",200,0,2);
 154   TH1D *P_dist=new TH1D("P_dist","",900,0.1,1);
 155   TH1D *P_dist_lowq=new TH1D("P_dist_lowq","",900,0.1,1);
 156   TProfile *AvgMult = new TProfile("AvgMult","",1,0.5,1.5, 0,2000,"");
 157   //
 158   TH2D *rstarPRF_dist_ss=new TH2D("rstarPRF_dist_ss","",20,0,1.0, 400,0,100);
 159   TH2D *tstar_dist_ss=new TH2D("tstar_dist_ss","",20,0,1.0, 400,0,100);
 160   //
 161   TH2D *rstarPRF_dist_os=new TH2D("rstarPRF_dist_os","",20,0,1.0, 400,0,100);
 162   TH2D *tstar_dist_os=new TH2D("tstar_dist_os","",20,0,1.0, 400,0,100);
 163
 164   TH2D *Num_Cos_ss=new TH2D("Num_Cos_ss","",20,0,1, 40,0,0.2);
 165   TH2D *Num_Cos_os=new TH2D("Num_Cos_os","",20,0,1, 40,0,0.2);
 166   TH2D *Num_CosFSI_ss=new TH2D("Num_CosFSI_ss","",20,0,1, 40,0,0.2);
 167   TH2D *Num_CosFSI_os=new TH2D("Num_CosFSI_os","",20,0,1, 40,0,0.2);
 168   TH2D *NumSq_CosFSI_ss=new TH2D("NumSq_CosFSI_ss","",20,0,1, 40,0,0.2);
 169   TH2D *NumSq_CosFSI_os=new TH2D("NumSq_CosFSI_os","",20,0,1, 40,0,0.2);
 170   TH2D *Num_PrimCosFSI_ss=new TH2D("Num_PrimCosFSI_ss","",20,0,1, 40,0,0.2);
 171   TH2D *Num_PrimCosFSI_os=new TH2D("Num_PrimCosFSI_os","",20,0,1, 40,0,0.2);
 172   //
 173   TH2D *Den_ss=new TH2D("Den_ss","",20,0,1, 40,0,0.2);
 174   TH2D *Den_os=new TH2D("Den_os","",20,0,1, 40,0,0.2);
 175   TH2D *DenEM_ss=new TH2D("DenEM_ss","",20,0,1, 40,0,0.2);
 176   TH2D *DenEM_os=new TH2D("DenEM_os","",20,0,1, 40,0,0.2);
 177   TH2D *LargeRpairs_ss=new TH2D("LargeRpairs_ss","",20,0,1, 40,0,0.2);
 178   TH2D *LargeRpairs_os=new TH2D("LargeRpairs_os","",20,0,1, 40,0,0.2);
 179   //
 180   TH2D *NumLambda1=new TH2D("NumLambda1","",40,0,0.2, 100,0,100);
 181   TH1D *DenLambda1=new TH1D("DenLambda1","",40,0,0.2);
 182   TH3D *NumLambda2_ss=new TH3D("NumLambda2_ss","",20,0,1, 40,0,0.2, 100,0,100);
 183   TH3D *NumLambda2_os=new TH3D("NumLambda2_os","",20,0,1, 40,0,0.2, 100,0,100);
 184   TH2D *DenLambda2_ss=new TH2D("DenLambda2_ss","",20,0,1, 40,0,0.2);
 185   TH2D *DenLambda2_os=new TH2D("DenLambda2_os","",20,0,1, 40,0,0.2);
 186   TH3D *NumLambda3_ss=new TH3D("NumLambda3_ss","",40,0,.2, 40,0,.2, 40,0,.2);
 187   TH3D *DenLambda3_ss=new TH3D("DenLambda3_ss","",40,0,.2, 40,0,.2, 40,0,.2);
 188   TH3D *NumLambda3_os=new TH3D("NumLambda3_os","",40,0,.2, 40,0,.2, 40,0,.2);
 189   TH3D *DenLambda3_os=new TH3D("DenLambda3_os","",40,0,.2, 40,0,.2, 40,0,.2);
 190   //
 191   TH3D *rstar3_dist_ss=new TH3D("rstar3_dist_ss","",200,0,100, 200,0,100, 200,0,100);
 192   TH3D *tstar3_dist_ss=new TH3D("tstar3_dist_ss","",200,0,100, 200,0,100, 200,0,100);
 193   //
 194   TH3D *r3num=new TH3D("r3num","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 195   TH3D *r3den1=new TH3D("r3den1","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 196   TH3D *r3den2=new TH3D("r3den2","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 197   TH3D *r3den3=new TH3D("r3den3","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 198   TH3D *r3numSq=new TH3D("r3numSq","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 199   TH3D *r3den1Sq=new TH3D("r3den1Sq","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 200   TH3D *r3den2Sq=new TH3D("r3den2Sq","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 201   TH3D *r3den3Sq=new TH3D("r3den3Sq","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 202   TH3D *r3numEn=new TH3D("r3numEn","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 203   TH3D *r3den1En=new TH3D("r3den1En","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 204   TH3D *r3den2En=new TH3D("r3den2En","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 205   TH3D *r3den3En=new TH3D("r3den3En","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 206   TH3D *r3numME=new TH3D("r3numME","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 207   TH3D *r3den1ME=new TH3D("r3den1ME","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 208   TH3D *r3den2ME=new TH3D("r3den2ME","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 209   TH3D *r3den3ME=new TH3D("r3den3ME","",qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 210   //
 211   TH1D *Num_Cos12=new TH1D("Num_Cos12","",40,0,0.2);
 212   TH1D *Num_Cos23=new TH1D("Num_Cos23","",40,0,0.2);
 213   TH1D *Num_Cos31=new TH1D("Num_Cos31","",40,0,0.2);
 214   TH1D *Den_Cos12=new TH1D("Den_Cos12","",40,0,0.2);
 215   TH1D *Den_Cos23=new TH1D("Den_Cos23","",40,0,0.2);
 216   TH1D *Den_Cos31=new TH1D("Den_Cos31","",40,0,0.2);
 217   //
 218   //
 219   TH2D *K2_ss = new TH2D("K2_ss","",20,0,1, qBins,0,qCutOff);
 220   TH2D *PlaneWF_ss = new TH2D("PlaneWF_ss","",20,0,1, qBins,0,qCutOff);
 221   TH2D *K2_os = new TH2D("K2_os","",20,0,1, qBins,0,qCutOff);
 222   TH2D *PlaneWF_os = new TH2D("PlaneWF_os","",20,0,1, qBins,0,qCutOff);
 223   //
 224   TH1D *PlaneWF3ss = new TH1D("PlaneWF3ss","",qBins,0,qCutOff);
 225   TH1D *PlaneWF3os = new TH1D("PlaneWF3os","",qBins,0,qCutOff);
 226   TH1D *K3ss = new TH1D("K3ss","",qBins,0,qCutOff);
 227   K3ss->GetXaxis()->SetTitle("Q3 (GeV/c)");
 228   K3ss->GetYaxis()->SetTitle("K_{3}");
 229   TH1D *K3os = new TH1D("K3os","",qBins,0,qCutOff);
 230   K3os->GetXaxis()->SetTitle("Q3 (GeV/c)");
 231   K3os->GetYaxis()->SetTitle("K_{3}");
 232   TH3D *K3ss_3D = new TH3D("K3ss_3D","", qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 233   K3ss_3D->GetXaxis()->SetTitle("q12");
 234   K3ss_3D->GetYaxis()->SetTitle("q23");
 235   K3ss_3D->GetZaxis()->SetTitle("q31");
 236   TH3D *PlaneWF3ss_3D = new TH3D("PlaneWF3ss_3D","", qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 237   TH3D *K3os_3D = new TH3D("K3os_3D","", qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 238   K3os_3D->GetXaxis()->SetTitle("q12");
 239   K3os_3D->GetYaxis()->SetTitle("q23");
 240   K3os_3D->GetZaxis()->SetTitle("q31");
 241   TH3D *PlaneWF3os_3D = new TH3D("PlaneWF3os_3D","", qBins,0,qCutOff, qBins,0,qCutOff, qBins,0,qCutOff);
 242
 243   TH1D *test=new TH1D("test","",200,-2*PI,2*PI);
 244
 245   ////////////////////////
 246
 247
 248   complex<double> binput(1.0,0.0);
 249   complex<double> CHG(0.,0.);// for Strong FSI in 3-body WF
 250   double kstar=0;
 251   double rstar=0;
 252   double kstar_1=0, kstar_2=0;
 253   double rstar_1=0, rstar_2=0;
 254   double phase_1=0, phase_2=0, phase_3=0, phase_4=0;
 255   double qinvMixedC_ss=0, qinvMixedC_os1=0, qinvMixedC_os2=0;
 256   double WeightWF_aa=1.0, WeightWF_bb=1.0;
 257   complex<double> CHG_1(0.,0.);
 258   complex<double> CHG_2(0.,0.);
 259   complex<double> CHG_3(0.,0.);
 260   complex<double> CHG_4(0.,0.);
 261   complex<double> AlphaWF[4];
 262   complex<double> BetaWF[4];
 263   double ThreePhase_1=0, ThreePhase_2=0, ThreePhase_3=0;
 264   double ThreePhase_4=0, ThreePhase_5=0, ThreePhase_6=0;
 265   double TwoPhase_12=0, TwoPhase_23=0, TwoPhase_31=0;
 266
 267   const int MAXPIONS = 2000;
 268   const int NevtBuff = 1;
 269   float P[NevtBuff][MAXPIONS][4]={{{0}}};
 270   float X[NevtBuff][MAXPIONS][4]={{{0}}};
 271   int PID[NevtBuff][MAXPIONS]={{0}};
 272   bool Primary[NevtBuff][MAXPIONS]={{0}};
 273   int Nparticles[NevtBuff]={0};
 274   Int_t ShuffledIndex[MAXPIONS]={0};// x-p decorrelation
 275   TVector3 P1_PRF;
 276   TVector3 P2_PRF;
 277   TVector3 P3_PRF;
 278   TVector3 X1_PRF;
 279   TVector3 X2_PRF;
 280   TVector3 X3_PRF;
 281   TVector3 X1_RF;
 282   TVector3 X2_RF;
 283   float P_dummy[4]={0};
 284
 285
 286
 287   // First Create B(rho,eta) and P(rho,eta) and h(eta) functions (for G hypergeometric functions (strong+Coulomb FSI))
 288   const int Btermslimit=500;// 500
 289   double rhoLimit = .05*100/FmToGeV;// 100 fm "limit"
 290   int rhoBins = rhoLimit/(0.001/FmToGeV);// this binning provides small changes between nearby bins
 291   double rhoStep = rhoLimit/rhoBins;
 292   double etaLimit = 1/(0.0005*BohrR);
 293   int etaBins = 2.0*etaLimit/(0.001);// this binning provides small changes between nearby bins
 294   double etaStep = (2*etaLimit)/etaBins;
 295   //cout<<rhoBins<<"  "<<etaBins<<"  "<<rhoLimit/rhoBins<<"  "<<etaLimit/etaBins<<endl;
 296   TH2D *B_histogram = new TH2D("B_histogram","",rhoBins,0,rhoLimit, etaBins,-etaLimit,etaLimit);
 297   TH2D *P_histogram = new TH2D("P_histogram","",rhoBins,0,rhoLimit, etaBins,-etaLimit,etaLimit);
 298   TH1D *h_histogram = new TH1D("h_histogram","",etaBins,-etaLimit,etaLimit);
 299
 300   for(int i=1; i<=rhoBins; i++){
 301     for(int j=1; j<=etaBins; j++){
 302       double rho = (i+0.5)*rhoStep;
 303       double eta = (j+0.5)*etaStep - etaLimit;// starts from negative values
 304       if(fabs(eta) < 0.0001) continue;
 305
 306       double Bfunc[Btermslimit]={0};
 307       double Pfunc[Btermslimit]={0};
 308       int StopPoint=Btermslimit;
 309       Bfunc[0]=1; Bfunc[1]=eta*rho;
 310       Pfunc[0]=1; Pfunc[1]=0;
 311       for(int ii=2; ii<Btermslimit; ii++) {
 312         Bfunc[ii] = (2*eta*rho*Bfunc[ii-1] - rho*rho*Bfunc[ii-2])/(double(ii*(ii+1.)));
 313         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));
 314         if(fabs(Bfunc[ii])+fabs(Bfunc[ii-1])+fabs(Pfunc[ii])+fabs(Pfunc[ii-1]) < 1.0e-7) {StopPoint=ii; break;}
 315       }
 316       //cout<<StopPoint+1<<endl;
 317       double Bsum=0, Psum=0;
 318       for(int ii=0; ii<StopPoint+1; ii++) {Bsum += Bfunc[ii]; Psum += Pfunc[ii];}
 319       B_histogram->SetBinContent(i,j, Bsum);
 320       P_histogram->SetBinContent(i,j, Psum);
 321     }
 322   }
 323   cout<<"Done preparing B and P histograms"<<endl;
 324   for(int j=1; j<=etaBins; j++){
 325     double eta = (j+0.5)*etaStep - etaLimit;// starts from negative values
 326     if(fabs(eta) < 0.0001) continue;
 327     h_histogram->SetBinContent(j, hFunction(eta));
 328   }
 329   cout<<"Done preparing h histogram"<<endl;
 330
 331
 332   int seconds = time(0);
 333
 334   for(int FileN=0; FileN<Nfiles; FileN++){
 335     cout<<"File #"<<FileN+1<<endl;
 336     //if(FileN<1260) continue;
 337     //if(FileN>230) continue;
 338     TString *fname=new TString("b");
 339     *fname += bValue;
 340     fname->Append("/event");
 341     *fname += FileN+1;
 342     fname->Append(".root");
 343     TFile *infile=new TFile(fname->Data(),"READ");
 344     if(!infile->IsOpen()) {cout<<fname->Data()<<" does not exist"<<endl; continue;}
 345     TTree *event_tree=(TTree*)infile->Get("events");
 346     TTree *particles_tree=(TTree*)infile->Get("particles");
 347
 348
 349     int Nevents = event_tree->GetEntries();
 350     int NpartList=0;
 351     int NpartList_past = 0;
 352     /////////////////////////////////////////
 353     // Create Pion list first
 354     for(int Evt=0; Evt<Nevents; Evt++){
 355
 356       //cout<<"Event # "<<Evt+1<<endl;
 357       event_tree->GetEntry(Evt);
 358
 359       TBranch *br=(TBranch*)event_tree->GetBranch("event");
 360       TLeaf *lf=(TLeaf*)br->GetLeaf("entries");
 361       NpartList_past += NpartList;
 362       NpartList = lf->GetValue();
 363
 364       /////////////////////////////////////
 365       // past event buffer
 366       for(int pastEvt=NevtBuff-1; pastEvt>0; pastEvt--){
 367
 368         for(int particle=0; particle<Nparticles[pastEvt-1]; particle++){
 369           P[pastEvt][particle][0] = P[pastEvt-1][particle][0];
 370           P[pastEvt][particle][1] = P[pastEvt-1][particle][1];
 371           P[pastEvt][particle][2] = P[pastEvt-1][particle][2];
 372           P[pastEvt][particle][3] = P[pastEvt-1][particle][3];
 373           //
 374           X[pastEvt][particle][0] = X[pastEvt-1][particle][0];
 375           X[pastEvt][particle][1] = X[pastEvt-1][particle][1];
 376           X[pastEvt][particle][2] = X[pastEvt-1][particle][2];
 377           X[pastEvt][particle][3] = X[pastEvt-1][particle][3];
 378           //
 379           PID[pastEvt][particle] = PID[pastEvt-1][particle];
 380           Primary[pastEvt][particle] = Primary[pastEvt-1][particle];
 381         }
 382
 383         Nparticles[pastEvt]=Nparticles[pastEvt-1];
 384       }
 385       Nparticles[0]=0;
 386       /////////////////////////////////////
 387
 388
 389       // Create Pion list first
 390       int Npions=0;
 391
 392       for(int index1=NpartList_past; index1<NpartList_past + NpartList; index1++){
 393
 394         if(Npions >= MAXPIONS) {cout<<"Too Many Pions!!!"<<endl; break;}
 395         particles_tree->GetEntry(index1);
 396         //
 397
 398         TLeaf *lf_pid=particles_tree->GetLeaf("pid");
 399         int pid = int(lf_pid->GetValue());
 400         if(abs(pid) != PIDPOI) continue;// pions only
 401         TLeaf *lf_decayed=(TLeaf*)particles_tree->GetLeaf("decayed");
 402         if(lf_decayed->GetValue() != 0) continue;
 403
 404
 405         TLeaf *lf_px=(TLeaf*)particles_tree->GetLeaf("px");
 406         TLeaf *lf_py=(TLeaf*)particles_tree->GetLeaf("py");
 407         TLeaf *lf_pz=(TLeaf*)particles_tree->GetLeaf("pz");
 408         TLeaf *lf_x=(TLeaf*)particles_tree->GetLeaf("x");
 409         TLeaf *lf_y=(TLeaf*)particles_tree->GetLeaf("y");
 410         TLeaf *lf_z=(TLeaf*)particles_tree->GetLeaf("z");
 411         TLeaf *lf_t=(TLeaf*)particles_tree->GetLeaf("t");
 412
 413         float px = lf_px->GetValue(); float py = lf_py->GetValue(); float pz = lf_pz->GetValue();
 414         float x = lf_x->GetValue(); float y = lf_y->GetValue(); float z = lf_z->GetValue();
 415         double t = lf_t->GetValue();
 416
 417         if(Gaussian){// Gaussian randomization
 418           bool goodchoice=kFALSE;
 419           while(goodchoice==kFALSE){
 420             x = NsigmaGauss*RValue*RandomNumber->Rndm();
 421             float height = RandomNumber->Rndm();
 422             if(SingleParticleGaussian->Eval(x) >= height) goodchoice=kTRUE;
 423           }
 424           if(RandomNumber->Rndm()<0.5) x = -x;
 425           //
 426           goodchoice=kFALSE;
 427           while(goodchoice==kFALSE){
 428             y = NsigmaGauss*RValue*RandomNumber->Rndm();
 429             float height = RandomNumber->Rndm();
 430             if(SingleParticleGaussian->Eval(y) >= height) goodchoice=kTRUE;
 431           }
 432           if(RandomNumber->Rndm()<0.5) y = -y;
 433           //
 434           goodchoice=kFALSE;
 435           while(goodchoice==kFALSE){
 436             z = NsigmaGauss*RValue*RandomNumber->Rndm();
 437             float height = RandomNumber->Rndm();
 438             if(SingleParticleGaussian->Eval(z) >= height) goodchoice=kTRUE;
 439           }
 440           if(RandomNumber->Rndm()<0.5) z = -z;
 441           //
 442           goodchoice=kFALSE;
 443           while(goodchoice==kFALSE){
 444             t = NsigmaGauss*RValue*RandomNumber->Rndm();
 445             float height = RandomNumber->Rndm();
 446             if(SingleParticleGaussian->Eval(t) >= height) goodchoice=kTRUE;
 447           }
 448           if(RandomNumber->Rndm()<0.5) t = -t;
 449           //t=0;
 450         }
 451
 452         float pt = sqrt(pow(px,2)+pow(py,2));
 453         float eta = asinh(pz/pt);
 454         float E = sqrt(pow(px,2)+pow(py,2)+pow(pz,2)+pow(massPOI,2));
 455
 456         if(pt<0.16 || pt>1.0) continue;
 457         if(fabs(eta)>0.8) continue;
 458         if(sqrt(pow(px,2)+pow(py,2)+pow(pz,2)) > 1.0) continue;
 459         //
 460         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;
 461         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;
 462         PID[0][Npions] = pid;
 463         float r = sqrt(pow(x,2)+pow(y,2)+pow(z,2));
 464         if(t < 1.e6) Primary[0][Npions]=kTRUE;// 1e6
 465         else {Primary[0][Npions]=kFALSE;}
 466
 467
 468         if(RemoveXP){
 469           if(r > 100./SF) continue;
 470         }
 471
 472         if(Primary[0][Npions]) P_dist->Fill(sqrt( pow(px,2) + pow(py,2) + pow(pz,2)));
 473
 474         if(Primary[0][Npions]){
 475           r_dist->Fill(r);
 476           x_dist->Fill(fabs(X[0][Npions][1])*FmToGeV);
 477           Pt_dist->Fill(pt);
 478         }
 479
 480         //if(Primary[0][Npions]==kFALSE) continue;
 481
 482         //
 483
 484
 485         Npions++;
 486       }
 487
 488       Nparticles[0]=Npions;
 489       AvgMult->Fill(1, Npions);
 490
 491       // Shuffle
 492       for (Int_t i = 0; i < Npions; i++) {
 493         // generate random index call
 494         int j = (Int_t) (gRandom->Rndm() * Npions);
 495         // temporarily store values from random index
 496         float t = X[0][j][0], x = X[0][j][1], y = X[0][j][2], z = X[0][j][3];
 497         float E = P[0][j][0], px = P[0][j][1], py = P[0][j][2], pz = P[0][j][3];
 498         int pid = PID[0][j];
 499         bool prim = Primary[0][j];
 500         //
 501         // Swap value locations
 502         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];
 503         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];
 504         PID[0][j] = PID[0][i];
 505         Primary[0][j] = Primary[0][i];
 506         //
 507         X[0][i][0] = t; X[0][i][1] = x; X[0][i][2] = y; X[0][i][3] = z;
 508         P[0][i][0] = E; P[0][i][1] = px; P[0][i][2] = py; P[0][i][3] = pz;
 509         PID[0][i] = pid;
 510         Primary[0][i] = prim;
 511       }
 512
 513       ///////////////////////////////////////////
 514       // de-correlate x-p
 515       if(RemoveXP){
 516         for (Int_t i = 0; i < Npions; i++) ShuffledIndex[i] = i;
 517         Shuffle(ShuffledIndex, 0, Npions-1);
 518         for (Int_t i = 0; i < Npions; i++) {// decorrelate x-p
 519           float t = X[0][i][0], x = X[0][i][1], y = X[0][i][2], z = X[0][i][3];
 520           X[0][i][0] = X[0][ShuffledIndex[i]][0];
 521           X[0][i][1] = X[0][ShuffledIndex[i]][1];
 522           X[0][i][2] = X[0][ShuffledIndex[i]][2];
 523           X[0][i][3] = X[0][ShuffledIndex[i]][3];
 524           //
 525           X[0][ShuffledIndex[i]][0] = t;
 526           X[0][ShuffledIndex[i]][1] = x;
 527           X[0][ShuffledIndex[i]][2] = y;
 528           X[0][ShuffledIndex[i]][3] = z;
 529         }
 530       }
 531
 532
 533       //////////////////////////////////////////
 534
 535
 536
 537       //////////////////////////////////////////////////
 538       // Start Correlation Analysis
 539
 540       for(int EN=0; EN<1; EN++){// current or 1st past event
 541
 542         for(int index1=0; index1<Nparticles[EN]; index1++){// current or past event pion loop
 543
 544           int start=index1+1;
 545           if(EN>0) start=0;
 546           for(int index2=start; index2<Npions; index2++){// current event pion loop (EN=0)
 547
 548             float kt = sqrt(pow(P[EN][index1][1]+P[0][index2][1],2)+pow(P[EN][index1][2]+P[0][index2][2],2))/2.;
 549             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));
 550             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));
 551
 552             if(qinv12 < 0.001) continue;
 553             if(qinv12 > qCutOff) continue;
 554
 555             if(EN>0){
 556               if(PID[EN][index1] == PID[0][index2]){
 557                 DenEM_ss->Fill(kt, qinv12);
 558               }else{
 559                 DenEM_os->Fill(kt, qinv12);
 560               }
 561               if(LambdaDirect) continue;
 562             }
 563
 564             BoostPRF(P[EN][index1], P_dummy, X[EN][index1], &X1_RF);
 565             DenLambda1->Fill(qinv12);
 566             if(X1_RF.Mag()<RstarMax) NumLambda1->Fill(qinv12, X1_RF.Mag()*FmToGeV);
 567             BoostPRF(P[0][index2], P_dummy, X[0][index2], &X2_RF);
 568             DenLambda1->Fill(qinv12);
 569             if(X2_RF.Mag()<RstarMax) NumLambda1->Fill(qinv12, X2_RF.Mag()*FmToGeV);
 570
 571
 572             P_dist_lowq->Fill(sqrt(pow(P[EN][index1][1],2)+pow(P[EN][index1][2],2)+pow(P[EN][index1][3],2)));
 573             P_dist_lowq->Fill(sqrt(pow(P[0][index2][1],2)+pow(P[0][index2][2],2)+pow(P[0][index2][3],2)));
 574
 575             // boost P into PRF
 576             BoostPRF(P[EN][index1], P[0][index2], P[EN][index1], &P1_PRF);
 577             BoostPRF(P[EN][index1], P[0][index2], P[0][index2], &P2_PRF);
 578             if(BoostPairsIntoPRF){// boost X into PRF
 579               BoostPRF(P[EN][index1], P[0][index2], X[EN][index1], &X1_PRF);
 580               BoostPRF(P[EN][index1], P[0][index2], X[0][index2], &X2_PRF);
 581             }else{
 582               X1_PRF[0]=X[EN][index1][1]; X1_PRF[1]=X[EN][index1][2]; X1_PRF[2]=X[EN][index1][3];
 583               X2_PRF[0]=X[0][index2][1]; X2_PRF[1]=X[0][index2][2]; X2_PRF[2]=X[0][index2][3];
 584             }
 585             //cout<<"out: "<<P1_PRF[0]<<"  "<<P2_PRF[0]<<endl;
 586             //cout<<"side: "<<P1_PRF[1]<<"  "<<P2_PRF[1]<<endl;
 587             //cout<<"long: "<<P1_PRF[2]<<"  "<<P2_PRF[2]<<" ++++++"<<endl;
 588             //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;
 589             //
 590             double rstar12_PRF = (X1_PRF-X2_PRF).Mag();
 591
 592             if(PID[EN][index1] == PID[0][index2]) DenLambda2_ss->Fill(kt, qinv12);
 593             else DenLambda2_os->Fill(kt, qinv12);
 594
 595             if(rstar12_PRF < RstarMin) continue;// removes most pairs from the same decay
 596
 597             if(rstar12_PRF < RstarMax){
 598               if(PID[EN][index1] == PID[0][index2]) NumLambda2_ss->Fill(kt, qinv12, rstar12_PRF*FmToGeV);
 599               else NumLambda2_os->Fill(kt, qinv12, rstar12_PRF*FmToGeV);
 600             }
 601
 602             if(rstar12_PRF > RstarMax) {
 603               if(PID[EN][index1]==PID[0][index2]) LargeRpairs_ss->Fill(kt, qinv12);
 604               else LargeRpairs_os->Fill(kt, qinv12);
 605               if(!LambdaDirect) continue;
 606             }
 607             float FVP12 = (P1_PRF-P2_PRF).Mag();
 608             FVP12 -= pow(P[EN][index1][0]-P[0][index2][0],2);
 609             double tstar12_PRF = sqrt(fabs(pow(rstar12_PRF,2)-FVP12));// fabs biases only 0.00001 % of the pairs (likely from same decay)
 610
 611
 612
 613             //
 614             double rho = rstar12_PRF * qinv12/2.;
 615             double phase = (P1_PRF-P2_PRF)*(X1_PRF-X2_PRF);
 616             double Cosphase = phase/(rstar12_PRF)/(qinv12);
 617
 618             double eta = 1/(qinv12/2.*BohrR);
 619             if(PID[EN][index1] != PID[0][index2]) eta = -eta;
 620             complex<double> ainput(0.0, -eta);
 621             complex<double> zinput12(0.0, rho + rho*Cosphase);
 622             complex<double> zinput21(0.0, rho - rho*Cosphase);
 623             complex<double> CHG12(1.0, 0);
 624             complex<double> CHG21(1.0, 0);
 625             double GamovFactor=1.0;
 626             double BEE = 0.;
 627             double WeightWF = 1.0;
 628             double WeightPlaneWF = 1.0;
 629             // Strong FSI
 630             complex<double> Psi_alpha(0., 0.);
 631             if(EN==0 && !LambdaDirect && qinv12 < qCutOff && rstar12_PRF < RstarMax){
 632               if(Primary[EN][index1]==kTRUE && Primary[0][index2]==kTRUE){
 633                 if(StrongFSI && PID[EN][index1]!=PID[0][index2]){
 634                   double kb = sqrt(pow(massPOI,2)+pow(qinv12/2.,2) - pow(masspi0,2));
 635                   double G2 = KinverseFunction(pow(qinv12/2.,2),1);
 636                   double G1 = KinverseFunction(pow(qinv12/2.,2),0);
 637                   double RK11 = 2/3.*G1 + 1/3.*G2;
 638                   double RK22 = 2/3.*G2 + 1/3.*G1;
 639                   double RK12 = -sqrt(2/9.)*(G1-G2);
 640                   complex<double> Chi(h_histogram->GetBinContent(h_histogram->GetXaxis()->FindBin(eta)), Gamov(eta)/(2*eta));
 641                   complex<double> C3 = RK11 - 2.0*Chi/BohrR;
 642                   complex<double> C5(RK22, -kb);
 643                   complex<double> C10 = C3*C5 - pow(RK12,2);
 644                   complex<double> fc_aa = C5/C10;
 645                   complex<double> G = P_histogram->GetBinContent(P_histogram->GetXaxis()->FindBin(rho), P_histogram->GetYaxis()->FindBin(eta));
 646                   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));
 647                   Psi_alpha = fc_aa*G/(rstar12_PRF);
 648                 }
 649
 650                 CHG12= conf_hyp(ainput,binput,zinput12);
 651                 CHG21= conf_hyp(ainput,binput,zinput21);
 652                 GamovFactor = Gamov(eta);
 653                 BEE = cos(phase);
 654
 655                 complex<double> planewave12(double(TMath::Cos(-rho*Cosphase)), double(TMath::Sin(-rho*Cosphase)));
 656                 complex<double> planewave21(planewave12.real(), -planewave12.imag());
 657                 //
 658                 complex<double> Full12(0.0,0.0);
 659                 Full12 += planewave12;
 660                 Full12 *= CHG12;
 661                 if(PID[EN][index1]!=PID[0][index2]) Full12 += Psi_alpha;
 662                 //
 663                 complex<double> Full21(0.0,0.0);
 664                 Full21 += planewave21;
 665                 Full21 *= CHG21;
 666                 //
 667                 complex<double> FullWF = Full12;
 668                 complex<double> FullPlaneWF = planewave12;
 669                 if(PID[EN][index1] == PID[0][index2]) {
 670                   FullWF += Full21;
 671                   FullWF /= sqrt(2.);
 672                   FullPlaneWF += planewave21;
 673                   FullPlaneWF /= sqrt(2.);
 674                 }
 675
 676                 WeightWF = GamovFactor*pow(abs(FullWF),2);
 677                 WeightPlaneWF = pow(abs(FullPlaneWF),2);
 678
 679                 // Fill undiluted histos
 680                 if(PID[EN][index1] == PID[0][index2]){
 681                   Num_PrimCosFSI_ss->Fill(kt, qinv12, WeightWF);
 682                 }else {
 683                   Num_PrimCosFSI_os->Fill(kt, qinv12, WeightWF);
 684                 }
 685
 686               }// primary
 687             }// same-event, LambdaDirect, rstar and qinv cut
 688
 689
 690             if(PID[EN][index1] == PID[0][index2]){// same-charge
 691               rstarPRF_dist_ss->Fill(kt, rstar12_PRF*FmToGeV);
 692               tstar_dist_ss->Fill(kt, fabs(X[EN][index1][0]-X[0][index2][0])*FmToGeV);
 693               Num_Cos_ss->Fill(kt, qinv12, BEE);
 694               Num_CosFSI_ss->Fill(kt, qinv12, WeightWF);
 695               NumSq_CosFSI_ss->Fill(kt, qinv12, pow(WeightWF,2));
 696               Den_ss->Fill(kt, qinv12);
 697               //
 698               K2_ss->Fill(kt, qinv12, WeightWF);
 699               PlaneWF_ss->Fill(kt, qinv12, WeightPlaneWF);
 700             }else {// mixed-charge
 701               rstarPRF_dist_os->Fill(kt, rstar12_PRF*FmToGeV);
 702               tstar_dist_os->Fill(kt, fabs(X[EN][index1][0]-X[0][index2][0])*FmToGeV);
 703               Num_Cos_os->Fill(kt, qinv12, BEE);
 704               Num_CosFSI_os->Fill(kt, qinv12, WeightWF);
 705               NumSq_CosFSI_os->Fill(kt, qinv12, pow(WeightWF,2));
 706               Den_os->Fill(kt, qinv12);
 707               //
 708               K2_os->Fill(kt, qinv12, WeightWF);
 709               PlaneWF_os->Fill(kt, qinv12, WeightPlaneWF);
 710             }
 711
 712
 713             if(!ThreeParticle) continue;
 714
 715
 716             TLorentzVector LP1(P[EN][index1][1], P[EN][index1][2], P[EN][index1][3], P[EN][index1][0]);
 717             TLorentzVector LX1(X[EN][index1][1], X[EN][index1][2], X[EN][index1][3], X[EN][index1][0]);
 718             TLorentzVector LP2(P[0][index2][1], P[0][index2][2], P[0][index2][3], P[0][index2][0]);
 719             TLorentzVector LX2(X[0][index2][1], X[0][index2][2], X[0][index2][3], X[0][index2][0]);
 720
 721
 722             if(!LambdaDirect) {if(!Primary[EN][index1] || !Primary[0][index2]) continue;}
 723             if(qinv12 > qCutOff) continue;
 724             //
 725             int EN2=0;
 726             int start3=index2+1;
 727             if(EN==1) {EN2=2; start3=0;}
 728
 729             for(int index3=start3; index3<Nparticles[EN2]; index3++){
 730
 731               if(!LambdaDirect && !Primary[EN2][index3]) continue;
 732
 733               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));
 734               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));
 735               if(qinv23 > qCutOff || qinv31 > qCutOff) continue;
 736               if(qinv23 < 0.001 || qinv31 < 0.001) continue;
 737
 738               if(PID[EN][index1] == PID[0][index2] && PID[EN][index1] == PID[EN2][index3]){
 739                 DenLambda3_ss->Fill(qinv12, qinv23, qinv31);
 740               }else{
 741                 DenLambda3_os->Fill(qinv12, qinv23, qinv31);
 742               }
 743
 744
 745               bool MixedCharge=kFALSE;
 746               if(PID[EN][index1] != PID[0][index2] || PID[EN][index1] != PID[EN2][index3] || PID[0][index2] != PID[EN2][index3]) MixedCharge=kTRUE;
 747
 748               if(EN2==2){
 749                 if(MixedCharge==kFALSE){
 750                   r3numME->Fill(qinv12, qinv23, qinv31);
 751                   r3den1ME->Fill(qinv12, qinv23, qinv31);
 752                   r3den2ME->Fill(qinv12, qinv23, qinv31);
 753                   r3den3ME->Fill(qinv12, qinv23, qinv31);
 754                 }
 755                 continue;
 756               }
 757
 758
 759
 760               // 12 boost
 761               BoostPRF(P[EN][index1], P[0][index2], P[EN2][index3], &P3_PRF);
 762               BoostPRF(P[EN][index1], P[0][index2], P[0][index2], &P2_PRF);
 763               BoostPRF(P[EN][index1], P[0][index2], P[EN][index1], &P1_PRF);
 764               X1_PRF[0]=X[EN][index1][1]; X1_PRF[1]=X[EN][index1][2]; X1_PRF[2]=X[EN][index1][3];// default CMS
 765               X2_PRF[0]=X[0][index2][1]; X2_PRF[1]=X[0][index2][2]; X2_PRF[2]=X[0][index2][3];// default CMS
 766               X3_PRF[0]=X[EN2][index3][1]; X3_PRF[1]=X[EN2][index3][2]; X3_PRF[2]=X[EN2][index3][3];// default CMS
 767
 768               //
 769               if(BoostPairsIntoPRF){
 770                 BoostPRF(P[EN][index1], P[0][index2], X[EN2][index3], &X3_PRF);
 771                 BoostPRF(P[EN][index1], P[0][index2], X[0][index2], &X2_PRF);
 772                 BoostPRF(P[EN][index1], P[0][index2], X[EN][index1], &X1_PRF);
 773               }
 774
 775
 776               rstar12_PRF = (X1_PRF-X2_PRF).Mag();
 777               if(rstar12_PRF > RstarMax) continue;
 778               if(rstar12_PRF < RstarMin) continue;
 779               double phase_p12_x12 = (P1_PRF-P2_PRF)*(X1_PRF-X2_PRF)/2.;
 780               double phase_p12_x32 = (P1_PRF-P2_PRF)*(X3_PRF-X2_PRF)/2.;
 781               double phase_p12_x13 = (P1_PRF-P2_PRF)*(X1_PRF-X3_PRF)/2.;
 782               complex<double> z_p12_x12(0.0, (P1_PRF-P2_PRF)*(X1_PRF-X2_PRF)/2. + qinv12/2. * (X1_PRF-X2_PRF).Mag());
 783               complex<double> z_p12_x21(0.0, (P1_PRF-P2_PRF)*(X2_PRF-X1_PRF)/2. + qinv12/2. * (X2_PRF-X1_PRF).Mag());
 784               complex<double> z_p12_x32(0.0, (P1_PRF-P2_PRF)*(X3_PRF-X2_PRF)/2. + qinv12/2. * (X3_PRF-X2_PRF).Mag());
 785               complex<double> z_p12_x13(0.0, (P1_PRF-P2_PRF)*(X1_PRF-X3_PRF)/2. + qinv12/2. * (X1_PRF-X3_PRF).Mag());
 786               complex<double> z_p12_x23(0.0, (P1_PRF-P2_PRF)*(X2_PRF-X3_PRF)/2. + qinv12/2. * (X2_PRF-X3_PRF).Mag());
 787               complex<double> z_p12_x31(0.0, (P1_PRF-P2_PRF)*(X3_PRF-X1_PRF)/2. + qinv12/2. * (X3_PRF-X1_PRF).Mag());
 788
 789
 790               // 23 boost
 791               BoostPRF(P[0][index2], P[EN2][index3], P[EN2][index3], &P3_PRF);
 792               BoostPRF(P[0][index2], P[EN2][index3], P[0][index2], &P2_PRF);
 793               BoostPRF(P[0][index2], P[EN2][index3], P[EN][index1], &P1_PRF);
 794               if(BoostPairsIntoPRF){
 795                 BoostPRF(P[0][index2], P[EN2][index3], X[EN2][index3], &X3_PRF);
 796                 BoostPRF(P[0][index2], P[EN2][index3], X[0][index2], &X2_PRF);
 797                 BoostPRF(P[0][index2], P[EN2][index3], X[EN][index1], &X1_PRF);
 798               }
 799               double rstar23_PRF = (X2_PRF-X3_PRF).Mag();
 800               if(rstar23_PRF > RstarMax) continue;
 801               if(rstar23_PRF < RstarMin) continue;
 802               double phase_p23_x13 = (P2_PRF-P3_PRF)*(X1_PRF-X3_PRF)/2.;
 803               double phase_p23_x23 = (P2_PRF-P3_PRF)*(X2_PRF-X3_PRF)/2.;
 804               double phase_p23_x21 = (P2_PRF-P3_PRF)*(X2_PRF-X1_PRF)/2.;
 805               double FVP23 = (P2_PRF-P3_PRF).Mag();
 806               FVP23 -= pow(P[0][index2][0]-P[EN2][index3][0],2);
 807               double tstar23_PRF = sqrt(fabs(pow(rstar23_PRF,2)-FVP23));
 808               //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));
 809               complex<double> z_p23_x23(0.0, (P2_PRF-P3_PRF)*(X2_PRF-X3_PRF)/2. + qinv23/2. * (X2_PRF-X3_PRF).Mag());
 810               complex<double> z_p23_x13(0.0, (P2_PRF-P3_PRF)*(X1_PRF-X3_PRF)/2. + qinv23/2. * (X1_PRF-X3_PRF).Mag());
 811               complex<double> z_p23_x21(0.0, (P2_PRF-P3_PRF)*(X2_PRF-X1_PRF)/2. + qinv23/2. * (X2_PRF-X1_PRF).Mag());
 812               complex<double> z_p23_x32(0.0, (P2_PRF-P3_PRF)*(X3_PRF-X2_PRF)/2. + qinv23/2. * (X3_PRF-X2_PRF).Mag());
 813               complex<double> z_p23_x31(0.0, (P2_PRF-P3_PRF)*(X3_PRF-X1_PRF)/2. + qinv23/2. * (X3_PRF-X1_PRF).Mag());
 814               complex<double> z_p23_x12(0.0, (P2_PRF-P3_PRF)*(X1_PRF-X2_PRF)/2. + qinv23/2. * (X1_PRF-X2_PRF).Mag());
 815
 816
 817               // 31 boost
 818               BoostPRF(P[EN2][index3], P[EN][index1], P[EN2][index3], &P3_PRF);
 819               BoostPRF(P[EN2][index3], P[EN][index1], P[0][index2], &P2_PRF);
 820               BoostPRF(P[EN2][index3], P[EN][index1], P[EN][index1], &P1_PRF);
 821               if(BoostPairsIntoPRF){
 822                 BoostPRF(P[EN2][index3], P[EN][index1], X[EN2][index3], &X3_PRF);
 823                 BoostPRF(P[EN2][index3], P[EN][index1], X[0][index2], &X2_PRF);
 824                 BoostPRF(P[EN2][index3], P[EN][index1], X[EN][index1], &X1_PRF);
 825               }
 826               double rstar31_PRF = (X3_PRF-X1_PRF).Mag();
 827               if(rstar31_PRF > RstarMax) continue;
 828               if(rstar31_PRF < RstarMin) continue;
 829               double phase_p31_x31 = (P3_PRF-P1_PRF)*(X3_PRF-X1_PRF)/2.;
 830               double phase_p31_x32 = (P3_PRF-P1_PRF)*(X3_PRF-X2_PRF)/2.;
 831               double phase_p31_x21 = (P3_PRF-P1_PRF)*(X2_PRF-X1_PRF)/2.;
 832               double FVP31 = (P3_PRF-P1_PRF).Mag();
 833               FVP31 -= pow(P[EN2][index3][0]-P[EN][index1][0],2);
 834               double tstar31_PRF = sqrt(fabs(pow(rstar31_PRF,2)-FVP31));
 835               //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));
 836               complex<double> z_p31_x31(0.0, (P3_PRF-P1_PRF)*(X3_PRF-X1_PRF)/2. + qinv31/2. * (X3_PRF-X1_PRF).Mag());
 837               complex<double> z_p31_x32(0.0, (P3_PRF-P1_PRF)*(X3_PRF-X2_PRF)/2. + qinv31/2. * (X3_PRF-X2_PRF).Mag());
 838               complex<double> z_p31_x13(0.0, (P3_PRF-P1_PRF)*(X1_PRF-X3_PRF)/2. + qinv31/2. * (X1_PRF-X3_PRF).Mag());
 839               complex<double> z_p31_x21(0.0, (P3_PRF-P1_PRF)*(X2_PRF-X1_PRF)/2. + qinv31/2. * (X2_PRF-X1_PRF).Mag());
 840               complex<double> z_p31_x12(0.0, (P3_PRF-P1_PRF)*(X1_PRF-X2_PRF)/2. + qinv31/2. * (X1_PRF-X2_PRF).Mag());
 841               complex<double> z_p31_x23(0.0, (P3_PRF-P1_PRF)*(X2_PRF-X3_PRF)/2. + qinv31/2. * (X2_PRF-X3_PRF).Mag());
 842
 843
 844
 845               ////////////////////////////////////////////
 846               ////////////////////////////////////////////
 847               double Q3 = sqrt(pow(qinv12,2) + pow(qinv23,2) + pow(qinv31,2));
 848
 849               if(!LambdaDirect){
 850                 int CP12=+1, CP23=+1, CP31=+1;
 851                 if(PID[EN][index1] != PID[0][index2]) CP12 = -1;
 852                 if(PID[0][index2] != PID[EN2][index3]) CP23 = -1;
 853                 if(PID[EN2][index3] != PID[EN][index1]) CP31 = -1;
 854
 855                 double Ac_12 = Gamov(CP12/(qinv12/2.*BohrR));
 856                 double Ac_23 = Gamov(CP23/(qinv23/2.*BohrR));
 857                 double Ac_31 = Gamov(CP31/(qinv31/2.*BohrR));
 858
 859                 //
 860                 complex<double> a_12(0.0, -CP12/(qinv12/2.*BohrR));
 861                 complex<double> a_23(0.0, -CP23/(qinv23/2.*BohrR));
 862                 complex<double> a_31(0.0, -CP31/(qinv31/2.*BohrR));
 863                 //
 864                 //
 865                 //
 866               TLorentzVector LP3(P[EN2][index3][1], P[EN2][index3][2], P[EN2][index3][3], P[EN2][index3][0]);
 867               TLorentzVector LX3(X[EN2][index3][1], X[EN2][index3][2], X[EN2][index3][3], X[EN2][index3][0]);
 868
 869               double phasePW = ( (LP1-LP2)*(LX1-LX2) + (LP1-LP3)*(LX1-LX3) + (LP2-LP3)*(LX2-LX3) )/3.;// base
 870               complex<double> planewave1(cos(phasePW), sin(phasePW));
 871               //
 872               phasePW = ( (LP1-LP2)*(LX2-LX1) + (LP1-LP3)*(LX2-LX3) + (LP2-LP3)*(LX1-LX3) )/3.;// 12 swap
 873               complex<double> planewave2(cos(phasePW), sin(phasePW));
 874               //
 875               phasePW = ( (LP1-LP2)*(LX3-LX2) + (LP1-LP3)*(LX3-LX1) + (LP2-LP3)*(LX2-LX1) )/3.;// 13 swap
 876               complex<double> planewave3(cos(phasePW), sin(phasePW));
 877               //
 878               phasePW = ( (LP1-LP2)*(LX1-LX3) + (LP1-LP3)*(LX1-LX2) + (LP2-LP3)*(LX3-LX2) )/3.;// 23 swap
 879               complex<double> planewave4(cos(phasePW), sin(phasePW));
 880               //
 881               phasePW = ( (LP1-LP2)*(LX2-LX3) + (LP1-LP3)*(LX2-LX1) + (LP2-LP3)*(LX3-LX1) )/3.;// 123 swap (1st type)
 882               complex<double> planewave5(cos(phasePW), sin(phasePW));
 883               //
 884               phasePW = ( (LP1-LP2)*(LX3-LX1) + (LP1-LP3)*(LX3-LX2) + (LP2-LP3)*(LX1-LX2) )/3.;// 123 swap (2nd type)
 885               complex<double> planewave6(cos(phasePW), sin(phasePW));
 886               //
 887               //
 888               //
 889               complex<double> CHG_p12_x12= conf_hyp(a_12, binput, z_p12_x12);
 890               complex<double> CHG_p12_x21= conf_hyp(a_12, binput, z_p12_x21);
 891               complex<double> CHG_p12_x32= conf_hyp(a_12, binput, z_p12_x32);
 892               complex<double> CHG_p12_x13= conf_hyp(a_12, binput, z_p12_x13);
 893               complex<double> CHG_p12_x23= conf_hyp(a_12, binput, z_p12_x23);
 894               complex<double> CHG_p12_x31= conf_hyp(a_12, binput, z_p12_x31);
 895               //
 896               complex<double> CHG_p31_x31= conf_hyp(a_31, binput, z_p31_x31);
 897               complex<double> CHG_p31_x32= conf_hyp(a_31, binput, z_p31_x32);
 898               complex<double> CHG_p31_x13= conf_hyp(a_31, binput, z_p31_x13);
 899               complex<double> CHG_p31_x21= conf_hyp(a_31, binput, z_p31_x21);
 900               complex<double> CHG_p31_x12= conf_hyp(a_31, binput, z_p31_x12);
 901               complex<double> CHG_p31_x23= conf_hyp(a_31, binput, z_p31_x23);
 902               //
 903               complex<double> CHG_p23_x23= conf_hyp(a_23, binput, z_p23_x23);
 904               complex<double> CHG_p23_x13= conf_hyp(a_23, binput, z_p23_x13);
 905               complex<double> CHG_p23_x21= conf_hyp(a_23, binput, z_p23_x21);
 906               complex<double> CHG_p23_x32= conf_hyp(a_23, binput, z_p23_x32);
 907               complex<double> CHG_p23_x31= conf_hyp(a_23, binput, z_p23_x31);
 908               complex<double> CHG_p23_x12= conf_hyp(a_23, binput, z_p23_x12);
 909
 910               ThreePhase_1 = X1_PRF*(P1_PRF-P2_PRF) + X2_PRF*(P2_PRF-P3_PRF) + X3_PRF*(P3_PRF-P1_PRF);
 911               ThreePhase_2 = X1_PRF*(P1_PRF-P2_PRF) + X3_PRF*(P2_PRF-P3_PRF) + X2_PRF*(P3_PRF-P1_PRF);
 912               ThreePhase_3 = X2_PRF*(P1_PRF-P2_PRF) + X1_PRF*(P2_PRF-P3_PRF) + X3_PRF*(P3_PRF-P1_PRF);
 913               ThreePhase_4 = X2_PRF*(P1_PRF-P2_PRF) + X3_PRF*(P2_PRF-P3_PRF) + X1_PRF*(P3_PRF-P1_PRF);
 914               ThreePhase_5 = X3_PRF*(P1_PRF-P2_PRF) + X1_PRF*(P2_PRF-P3_PRF) + X2_PRF*(P3_PRF-P1_PRF);
 915               ThreePhase_6 = X3_PRF*(P1_PRF-P2_PRF) + X2_PRF*(P2_PRF-P3_PRF) + X1_PRF*(P3_PRF-P1_PRF);
 916               TwoPhase_12 = (P1_PRF-P2_PRF)*(X1_PRF-X2_PRF);
 917               TwoPhase_23 = (P2_PRF-P3_PRF)*(X2_PRF-X3_PRF);
 918               TwoPhase_31 = (P3_PRF-P1_PRF)*(X3_PRF-X1_PRF);
 919               //
 920               //
 921
 922               // Strong FSI
 923               if(MixedCharge){
 924                 if(CP12==+1){// 1 and 2 are same-charge
 925                   kstar_1 = qinv31/2.; kstar_2 = qinv23/2.;
 926                   rstar_1 = rstar31_PRF; rstar_2 = rstar23_PRF;
 927                   phase_1 = phase_p31_x31; phase_2 = phase_p31_x32; phase_3 = phase_p23_x13; phase_4 = phase_p23_x23;
 928                   CHG_1 = CHG_p31_x31;
 929                   CHG_2 = CHG_p31_x32;
 930                   CHG_3 = CHG_p23_x13;
 931                   CHG_4 = CHG_p23_x23;
 932                   qinvMixedC_ss=qinv12; qinvMixedC_os1=qinv31; qinvMixedC_os2=qinv23;
 933                 }else if(CP31==+1){// 1 and 3 are same-charge
 934                   kstar_1 = qinv12/2.; kstar_2 = qinv23/2.;
 935                   rstar_1 = rstar12_PRF; rstar_2 = rstar23_PRF;
 936                   phase_1 = phase_p12_x12; phase_2 = phase_p12_x32; phase_3 = phase_p23_x21; phase_4 = phase_p23_x23;
 937                   CHG_1 = CHG_p12_x12;
 938                   CHG_2 = CHG_p12_x32;
 939                   CHG_3 = CHG_p23_x21;
 940                   CHG_4 = CHG_p23_x23;
 941                   qinvMixedC_ss=qinv31; qinvMixedC_os1=qinv12; qinvMixedC_os2=qinv23;
 942                 }else {// 2 and 3 are same-charge
 943                   kstar_1 = qinv12/2.; kstar_2 = qinv31/2.;
 944                   rstar_1 = rstar12_PRF; rstar_2 = rstar31_PRF;
 945                   phase_1 = phase_p12_x12; phase_2 = phase_p12_x13; phase_3 = phase_p31_x21; phase_4 = phase_p31_x31;
 946                   CHG_1 = CHG_p12_x12;
 947                   CHG_2 = CHG_p12_x13;
 948                   CHG_3 = CHG_p31_x21;
 949                   CHG_4 = CHG_p31_x31;
 950                   qinvMixedC_ss=qinv23; qinvMixedC_os1=qinv31; qinvMixedC_os2=qinv31;
 951                 }
 952
 953
 954                 for(int piece=0; piece<4; piece++){
 955                   if(piece==0) {// P31X31
 956                     kstar=kstar_1;
 957                     rstar=rstar_1;
 958                     phase=phase_1;
 959                     CHG=CHG_1;
 960                   }else if(piece==1) {// P31X32
 961                     kstar=kstar_1;
 962                     rstar=rstar_2;
 963                     phase=phase_2;
 964                     CHG=CHG_2;
 965                   }else if(piece==2) {// P23X13
 966                     kstar=kstar_2;
 967                     rstar=rstar_1;
 968                     phase=phase_3;
 969                     CHG=CHG_3;
 970                   }else {// P23X23
 971                     kstar=kstar_2;
 972                     rstar=rstar_2;
 973                     phase=phase_4;
 974                     CHG=CHG_4;
 975                   }
 976
 977                   rho = kstar*rstar;
 978                   eta = -1./(kstar*BohrR);
 979                   double kb = sqrt(pow(massPOI,2)+pow(kstar,2) - pow(masspi0,2));
 980                   double G2 = KinverseFunction(pow(kstar,2),1);
 981                   double G1 = KinverseFunction(pow(kstar,2),0);
 982                   double RK11 = 2/3.*G1 + 1/3.*G2;
 983                   double RK22 = 2/3.*G2 + 1/3.*G1;
 984                   double RK12 = -sqrt(2/9.)*(G1-G2);
 985                   complex<double> Chi(h_histogram->GetBinContent(h_histogram->GetXaxis()->FindBin(eta)), Gamov(eta)/(2*eta));
 986                   complex<double> C3 = RK11 - 2.0*Chi/BohrR;
 987                   complex<double> C5(RK22, -kb);
 988                   complex<double> C10 = C3*C5 - pow(RK12,2);
 989                   complex<double> fc_aa = C5/C10;
 990                   complex<double> fc_ab = -RK12/C10;
 991                   complex<double> G = P_histogram->GetBinContent(P_histogram->GetXaxis()->FindBin(rho), P_histogram->GetYaxis()->FindBin(eta));
 992                   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));
 993                   complex<double> PW(cos(phase), sin(phase));
 994                   AlphaWF[piece] = CHG + PW*fc_aa*G/rstar;
 995                   complex<double> spherewave(cos(rstar*kb), sin(rstar*kb));
 996                   BetaWF[piece] = PW*fc_ab*sqrt(masspi0/massPOI)*spherewave/rstar;
 997
 998                 }// piece
 999
1000               }// Mixed-charge case with Strong FSI
1001
1002               complex<double> FullFSI_WF[2];// aa, bb
1003               if(MixedCharge==kTRUE){
1004                 if(StrongFSI==kFALSE) {
1005                   FullFSI_WF[0] = planewave1*CHG_p12_x12*CHG_p31_x31*CHG_p23_x23;
1006                   if(CP12==+1){// 1 and 2 are same-charge
1007                     FullFSI_WF[0] += planewave2*CHG_p12_x21*CHG_p31_x32*CHG_p23_x13;// 12 symmetrization
1008                   }else if(CP31==+1){// 1 and 3 are same-charge
1009                     FullFSI_WF[0] += planewave3*CHG_p12_x32*CHG_p31_x13*CHG_p23_x21;// 13 symmetrization
1010                   }else{
1011                     FullFSI_WF[0] += planewave4*CHG_p12_x13*CHG_p31_x21*CHG_p23_x32;// 23 symmetrization
1012                   }
1013                   FullFSI_WF[1] = 0.;
1014                 }else {// strong FSI
1015                   if(CP12==+1){// 1 and 2 are same-charge
1016                     FullFSI_WF[0] = planewave1*CHG_p12_x12*AlphaWF[0]*AlphaWF[3];
1017                     FullFSI_WF[1] = planewave1*CHG_p12_x12*BetaWF[0]*BetaWF[3];
1018                     FullFSI_WF[0] += planewave2*CHG_p12_x21*AlphaWF[1]*AlphaWF[2];// 12 symmetrization
1019                     FullFSI_WF[1] += planewave2*CHG_p12_x21*BetaWF[1]*BetaWF[2];// 12 symmetrization
1020                   }else if(CP31==+1){// 1 and 3 are same-charge
1021                     FullFSI_WF[0] = planewave1*CHG_p31_x31*AlphaWF[0]*AlphaWF[3];
1022                     FullFSI_WF[1] = planewave1*CHG_p31_x31*BetaWF[0]*BetaWF[3];
1023                     FullFSI_WF[0] += planewave3*CHG_p31_x13*AlphaWF[1]*AlphaWF[2];// 13 symmetrization
1024                     FullFSI_WF[1] += planewave3*CHG_p31_x13*BetaWF[1]*BetaWF[2];// 13 symmetrization
1025                   }else{
1026                     FullFSI_WF[0] = planewave1*CHG_p23_x23*AlphaWF[0]*AlphaWF[3];
1027                     FullFSI_WF[1] = planewave1*CHG_p23_x23*BetaWF[0]*BetaWF[3];
1028                     FullFSI_WF[0] += planewave4*CHG_p23_x32*AlphaWF[1]*AlphaWF[2];// 13 symmetrization
1029                     FullFSI_WF[1] += planewave4*CHG_p23_x32*BetaWF[1]*BetaWF[2];// 13 symmetrization
1030                   }
1031                 }
1032                 FullFSI_WF[0] /= sqrt(2.); FullFSI_WF[1] /= sqrt(2.);
1033               }else {
1034                  FullFSI_WF[0] = planewave1*CHG_p12_x12*CHG_p31_x31*CHG_p23_x23;// base
1035                  FullFSI_WF[0] += planewave2*CHG_p12_x21*CHG_p31_x32*CHG_p23_x13;// 12 symmetrization
1036                  FullFSI_WF[0] += planewave3*CHG_p12_x32*CHG_p31_x13*CHG_p23_x21;// 13
1037                  FullFSI_WF[0] += planewave4*CHG_p12_x13*CHG_p31_x21*CHG_p23_x32;// 23
1038                  FullFSI_WF[0] += planewave5*CHG_p12_x23*CHG_p31_x12*CHG_p23_x31;// 123
1039                  FullFSI_WF[0] += planewave6*CHG_p12_x31*CHG_p31_x23*CHG_p23_x12;// 123
1040                  FullFSI_WF[1] = 0.;
1041                  FullFSI_WF[0] /= sqrt(6.); FullFSI_WF[1] /= sqrt(6.);
1042               }
1043
1044               complex<double> FullPlaneWF = planewave1;// base
1045               if(MixedCharge && CP12==+1) FullPlaneWF += planewave2;
1046               if(MixedCharge && CP31==+1) FullPlaneWF += planewave3;
1047               if(MixedCharge && CP23==+1) FullPlaneWF += planewave4;
1048               if(!MixedCharge) {
1049                 FullPlaneWF += planewave2 + planewave3 + planewave4 + planewave5 + planewave6;
1050                 FullPlaneWF /= sqrt(6.);
1051               }else FullPlaneWF /= sqrt(2.);
1052
1053               WeightWF_aa = pow(abs(FullFSI_WF[0]),2);
1054               WeightWF_bb = pow(abs(FullFSI_WF[1]),2);
1055               WeightWF_aa *= Ac_12*Ac_23*Ac_31;
1056               WeightWF_bb *= Ac_12*Ac_23*Ac_31;
1057               WeightPlaneWF = pow(abs(FullPlaneWF),2);
1058
1059               }// !LambdaDirect
1060               ////////////////////////////////////////////
1061               ////////////////////////////////////////////
1062               //
1063
1064
1065               if(MixedCharge){
1066                 K3os->Fill(Q3, WeightWF_aa + WeightWF_bb);
1067                 PlaneWF3os->Fill(Q3, WeightPlaneWF);
1068                 K3os_3D->Fill(qinvMixedC_ss, qinvMixedC_os1, qinvMixedC_os2, WeightWF_aa + WeightWF_bb);
1069                 PlaneWF3os_3D->Fill(qinvMixedC_ss, qinvMixedC_os1, qinvMixedC_os2, WeightPlaneWF);
1070                 NumLambda3_os->Fill(qinv12, qinv23, qinv31);
1071               }else{
1072                 K3ss->Fill(Q3, WeightWF_aa + WeightWF_bb);
1073                 PlaneWF3ss->Fill(Q3, WeightPlaneWF);
1074                 K3ss_3D->Fill(qinv12, qinv23, qinv31, WeightWF_aa + WeightWF_bb);
1075                 PlaneWF3ss_3D->Fill(qinv12, qinv23, qinv31, WeightPlaneWF);
1076                 NumLambda3_ss->Fill(qinv12, qinv23, qinv31);
1077                 rstar3_dist_ss->Fill(rstar12_PRF, rstar23_PRF, rstar31_PRF);
1078                 //rstar3_dist_ss->Fill(rstar12_CMS, rstar23_CMS, rstar31_CMS);
1079                 tstar3_dist_ss->Fill(tstar12_PRF, tstar23_PRF, tstar31_PRF);
1080                 //
1081                 double cosThreePhase = 1/6. * (cos(ThreePhase_1)+cos(ThreePhase_2)+cos(ThreePhase_3)+cos(ThreePhase_4)+cos(ThreePhase_5)+cos(ThreePhase_6));
1082                 r3num->Fill(qinv12, qinv23, qinv31, 2*cosThreePhase);
1083                 r3den1->Fill(qinv12, qinv23, qinv31, cos(TwoPhase_12));
1084                 r3den2->Fill(qinv12, qinv23, qinv31, cos(TwoPhase_23));
1085                 r3den3->Fill(qinv12, qinv23, qinv31, cos(TwoPhase_31));
1086                 r3numSq->Fill(qinv12, qinv23, qinv31, pow(2*cosThreePhase,2));
1087                 r3den1Sq->Fill(qinv12, qinv23, qinv31, pow(cos(TwoPhase_12),2));
1088                 r3den2Sq->Fill(qinv12, qinv23, qinv31, pow(cos(TwoPhase_23),2));
1089                 r3den3Sq->Fill(qinv12, qinv23, qinv31, pow(cos(TwoPhase_31),2));
1090                 r3numEn->Fill(qinv12, qinv23, qinv31);
1091                 r3den1En->Fill(qinv12, qinv23, qinv31);
1092                 r3den2En->Fill(qinv12, qinv23, qinv31);
1093                 r3den3En->Fill(qinv12, qinv23, qinv31);
1094                 //
1095                 Num_Cos12->Fill(qinv12,cos(TwoPhase_12));
1096                 Num_Cos23->Fill(qinv23,cos(TwoPhase_23));
1097                 Num_Cos31->Fill(qinv31,cos(TwoPhase_31));
1098                 Den_Cos12->Fill(qinv12);
1099                 Den_Cos23->Fill(qinv23);
1100                 Den_Cos31->Fill(qinv31);
1101               }
1102
1103
1104             }// index3
1105
1106           }// index2
1107
1108         }// index1
1109
1110       }// EN buffer
1111
1112     }// Nevents
1113
1114     infile->Close();
1115   }// Nfiles
1116
1117   outfile->cd();
1118   //
1119   r_dist->Write();
1120   x_dist->Write();
1121   Pt_dist->Write();
1122   P_dist->Write();
1123   P_dist_lowq->Write();
1124   AvgMult->Write();
1125
1126   rstarPRF_dist_ss->Write();
1127   tstar_dist_ss->Write();
1128   rstarPRF_dist_os->Write();
1129   tstar_dist_os->Write();
1130   NumLambda1->Write();
1131   DenLambda1->Write();
1132   NumLambda2_ss->Write();
1133   NumLambda2_os->Write();
1134   DenLambda2_ss->Write();
1135   DenLambda2_os->Write();
1136   NumLambda3_ss->Write();
1137   DenLambda3_ss->Write();
1138   NumLambda3_os->Write();
1139   DenLambda3_os->Write();
1140   Den_ss->Write();
1141   Den_os->Write();
1142   DenEM_ss->Write();
1143   DenEM_os->Write();
1144   Num_Cos_ss->Write();
1145   Num_Cos_os->Write();
1146   Num_CosFSI_ss->Write();
1147   Num_CosFSI_os->Write();
1148   NumSq_CosFSI_ss->Write();
1149   NumSq_CosFSI_os->Write();
1150   LargeRpairs_ss->Write();
1151   LargeRpairs_os->Write();
1152   //
1153   Num_PrimCosFSI_ss->Write();
1154   Num_PrimCosFSI_os->Write();
1155   //
1156   rstar3_dist_ss->Write();
1157   tstar3_dist_ss->Write();
1158   //
1159   K2_ss->Write();
1160   PlaneWF_ss->Write();
1161   K2_os->Write();
1162   PlaneWF_os->Write();
1163   //
1164   K3os->Write();
1165   PlaneWF3os->Write();
1166   K3os_3D->Write();
1167   PlaneWF3os_3D->Write();
1168   K3ss->Write();
1169   PlaneWF3ss->Write();
1170   K3ss_3D->Write();
1171   PlaneWF3ss_3D->Write();
1172   //
1173   r3num->Write();
1174   r3den1->Write();
1175   r3den2->Write();
1176   r3den3->Write();
1177   r3numSq->Write();
1178   r3den1Sq->Write();
1179   r3den2Sq->Write();
1180   r3den3Sq->Write();
1181   r3numEn->Write();
1182   r3den1En->Write();
1183   r3den2En->Write();
1184   r3den3En->Write();
1185   r3numME->Write();
1186   r3den1ME->Write();
1187   r3den2ME->Write();
1188   r3den3ME->Write();
1189   //
1190   Num_Cos12->Write();
1191   Num_Cos23->Write();
1192   Num_Cos31->Write();
1193   Den_Cos12->Write();
1194   Den_Cos23->Write();
1195   Den_Cos31->Write();
1196   //
1197   test->Write();
1198   //
1199   outfile->Close();
1200
1201   seconds = time(0) - seconds;
1202   cout<<"Minutes = "<<seconds/60.<<endl;
1203
1204
1205 }
1206 float Gamov(float eta){
1207   return (2*PI*eta/(exp(2*PI*eta)-1));
1208 }
1209 float fact(float n){
1210   return (n < 1.00001) ? 1 : fact(n - 1) * n;
1211 }
1212 void Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
1213 {
1214   Int_t j, k;
1215   Int_t a = i2 - i1;
1216   for (Int_t i = i1; i < i2+1; i++) {
1217     j = (Int_t) (gRandom->Rndm() * a);
1218     k = iarr[j];
1219     iarr[j] = iarr[i];
1220     iarr[i] = k;
1221   }
1222 }
1223 void BoostPRF(float P1[4], float P2[4], float V[4], TVector3 *T){
1224   // P1 is particle 1 momentum in reaction CMS (0=E,1=px,2=py,3=pz)
1225   // P2 is particle 2 momentum in reaction CMS (0=E,1=px,2=py,3=pz)
1226   // V is the vector to be transformed (V is in CMS)
1227   // T is the PRF vector
1228
1229   float MT=sqrt(pow(P1[0]+P2[0],2)-pow(P1[3]+P2[3],2));
1230   float PT=sqrt(pow(P1[1]+P2[1],2)+pow(P1[2]+P2[2],2));
1231   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));
1232   // LCMS
1233   T->SetX( ((P1[1]+P2[1])*V[1] + (P1[2]+P2[2])*V[2])/PT );
1234   T->SetY( ((P1[1]+P2[1])*V[2] - (P1[2]+P2[2])*V[1])/PT );
1235   T->SetZ( ((P1[0]+P2[0])*V[3] - (P1[3]+P2[3])*V[0])/MT );
1236   // PRF
1237   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]) );
1238 }
1239 // h function from Lednicky, phys. of part. and nuclei 40:307 (2009)
1240 double hFunction(double eta){
1241   double DS=1.0;
1242   double N=0;
1243   double S=0;
1244   double returnValue=0;
1245   if(eta < 3.0){
1246     while(DS > 1.0e-13){
1247       N++;
1248       DS = 1/(N*(pow(N/eta,2)+1.0));
1249       S += DS;
1250     }
1251     returnValue = (S - EulerC - log(fabs(eta)));
1252   }else{// small kstar, high eta
1253     double etaSquared=pow(eta,2);
1254     double etaPowered=etaSquared;
1255     for(int ii=0; ii<9; ii++){// 9 was maximum value for Lednicky's code
1256       returnValue += 1/etaPowered * hcoef[ii];
1257       etaPowered *= etaSquared;
1258     }
1259   }
1260   return returnValue;
1261 }
1262 double KinverseFunction(double kstarSq, int J){
1263   double ESq = kstarSq + pow(massPOI,2);
1264   double E = sqrt(ESq);
1265   double xSq = kstarSq/pow(massPOI,2);
1266   double Kinverse = E*(4*ESq-fzero[J][4])/(4*pow(massPOI,2)-fzero[J][4]);
1267   Kinverse /= fzero[J][0] + fzero[J][1]*xSq + fzero[J][2]*pow(xSq,2) + fzero[J][3]*pow(xSq,3);
1268   return Kinverse;
1269 }
1270