| 9129699a |
1 | #include <math.h> |
| 2 | #include <time.h> |
| 3 | #include <stdio.h> |
| 4 | #include <stdlib.h> |
| 5 | #include <Riostream.h> |
| 6 | #include <complex> |
| 7 | |
| 8 | #include "TObject.h" |
| 9 | #include "TTree.h" |
| 10 | #include "TBranch.h" |
| 11 | #include "TLeaf.h" |
| 12 | #include "TVector2.h" |
| 13 | #include "TVector3.h" |
| 14 | #include "TFile.h" |
| 15 | #include "TString.h" |
| 16 | #include "TF1.h" |
| 17 | #include "TH1.h" |
| 18 | #include "TH2.h" |
| 19 | #include "TH3.h" |
| 20 | #include "TProfile.h" |
| 21 | #include "TProfile2D.h" |
| 22 | #include "TMath.h" |
| 23 | #include "TText.h" |
| 24 | #include "TRandom3.h" |
| 25 | #include "TArray.h" |
| 26 | #include "TLegend.h" |
| 27 | #include "TStyle.h" |
| 28 | #include "TMinuit.h" |
| 29 | #include "TChain.h" |
| 30 | #include "TLorentzVector.h" |
| 31 | |
| 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 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame