// Plot the TF3 Integral function used for the fit of 3 GeV data // Fitting one particle species. // root -x 'BWModified_Proton1.C("proton_pt_data_3Gev.dat")' //#include "TF23.h" // class projectiong a TF3(x,y,z) in a TF2 (x,z) struct TF3Projector { TF3Projector(TF3 * f, double y) : fFunc(f), fX(), fP() { fX[0] = y;} double operator() (const double *x, const double *p) { fX[1] = x[0]; fX[2] = x[1]; fP[0] = p[0]; fP[1] = p[1]; fP[2] = p[2]; fP[3] = p[3]; fP[4] = p[4]; return (*fFunc)( fX,fP ); } TF3 * fFunc; double fX[3]; double fP[5]; }; double Gaussian(double y, double sigma,double mean){ double f1 = pow((y-mean),2); double f2 = 2*pow(sigma,2); double f3 = -f1/f2; return TMath::Exp(f3); } Double_t fRy(Double_t *x,Double_t *par){ //par[0] = Tkin //par[1] = betaS //par[2] = n //par[3] = sigma //par[4] = Normalization //x[0] = pT //x[1] = y //x[2] = r const double R = 10.; const double m = 0.938; const double ycm = 1.045; double beta = par[1] * pow(x[2]/R, par[2]); double rho = TMath::ATanH(beta); double mT = sqrt(x[0]*x[0] + m*m); return par[4] * (mT * x[2] * TMath::BesselI0(x[0] * TMath::SinH(rho) / par[0]) * TMath::CosH(x[1]) * TMath::Exp(-mT * TMath::CosH(x[1]) * TMath::CosH(rho) / par[0]) * Gaussian(x[1],par[3],ycm)); } double getYield(Double_t *x,Double_t *par){ TF3 * f3 = new TF3("fRy",fRy); /* TF3Projector * p = new TF3Projector(f3,par[0]); TF2 *fpT = new TF2("fpT",p,0.0,10.0,-2.04,2.04,4);//fpT should have only 4 pars. after integration. fpT->SetParameter(0,par[1]); fpT->SetParameter(1,par[2]); fpT->SetParameter(2,par[3]); fpT->SetParameter(3,par[4]); return fpT->Integral(0.,10.,-2.04,2.04); */ TF3Projector * p = new TF3Projector(f3,x[0]); TF2 *fpT = new TF2("fpT",p,-2.04,2.04,0,10,5);//fpT should have only 5 pars. after integration. fpT->SetParameter(0,par[0]); fpT->SetParameter(1,par[1]); fpT->SetParameter(2,par[2]); fpT->SetParameter(3,par[3]); fpT->SetParameter(4,par[4]); double f = fpT->Integral(-2.04,2.04,0,10); delete fpT; delete p; delete f3; return f; } void BWModified_Proton1(Char_t *inFile){ ifstream mFile; mFile.open(inFile); if (mFile.is_open() == 0) { cout << " file not found! " << endl; return 0; } Char_t buf[256]; const unsigned int MaxNBins = 100; Double_t x[MaxNBins]; Double_t y[MaxNBins]; Double_t ey[MaxNBins]; Int_t Bins; Int_t j=0; while (!mFile.eof()){ mFile >> x[j] >> y[j] >> ey[j] ; cout << x[j] << " " <GetFrame()->SetFillColor(21); c1->GetFrame()->SetBorderSize(12); //TH2F *fhist = new TH2F("fhist","",50,0.15,0.7,200,0.5e1,1e2); TH2F *fhist = new TH2F("fhist","",50,0.6,1.3,200,4,25); gStyle->SetOptStat(0); fhist->SetTitle("Blast Wave Fit for p;p_{t} (GeV/c);#frac{1}{2#pi N_{evt}} #frac{1}{p_{T}} #frac{d^{2}N}{dp_{T}dy} (GeV/c)^{-2};"); fhist->GetXaxis()->CenterTitle(); fhist->GetYaxis()->CenterTitle(); c1->SetLogy(); fhist->Draw(""); Double_t par0[9] = {0.09,0.1,0.1,0.1,0.1, 0.1,0.1,0.1, 0.1}; // Temperature Double_t par1[9] = {0.8,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5}; // Beta_s Double_t par2[9] = {1.0,1, 1,1.0,1,1,1,1,1}; // n Double_t par3[9] = {1e5,1e4,1e5,1e6,1e5,1e5,1e5,1e5,1e5}; // Normalization constant TVirtualFitter::SetDefaultFitter("Minuit"); TVirtualFitter::SetMaxIterations(500000); TF1 *fns = NULL; //initialize /* fns = new TF1("fns",getYield,.1,.5,4); fns->SetParNames("T","beta_s" , "n" , "norm"); fns->SetParameters(par0[0],par1[0],par2[0],par3[0]); */ fns = new TF1("fns",getYield,.1,2,5); fns->SetParNames("T","beta_s" , "n" , "sigma", "norm"); fns->SetParameters(0.073,0.63,0.5,0.5,1.1e7); fns->SetLineColor(kBlue); fns->SetParLimits(2,0,2); fns->FixParameter(3,0.5); grBW->Fit("fns","S"); grBW->Draw("P"); grBW->SetMarkerColor(kRed); grBW->SetMarkerStyle(21); }