Rho photonuclear STARLight MC

pair_pt.C

#include "TFile.h"
#include "TH1F.h"
#include "TH1D.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TVector3.h"
#include "TLegend.h"
#include "TF1.h"
#include "TRandom3.h"
#include "TStyle.h"
#include "TLorentzVector.h"


TF1 * fFFQ = 0;
TF1 * fgkt = 0;

TVector3 impact_parameter( 96, 0, 0 );

TH2 * hDeltaPhiPt;


// boost a daughter vector according to the parent (lab) vector boost
void applyBoost( TLorentzVector &_parent_lv, TLorentzVector &_d_lv ){

    float betaX = _parent_lv.Px() / _parent_lv.E();
    float betaY = _parent_lv.Py() / _parent_lv.E();
    float betaZ = _parent_lv.Pz() / _parent_lv.E();

    _d_lv.Boost(betaX,betaY,betaZ);
}

/* Two Body Decay in the rest frame of the parent
 *
 * The decay is computed and then the daughters are boosted into the frame of the parent
 */
vector<TLorentzVector> twoBodyDecay( TLorentzVector _parent_lv, double M_d1, double M_d2, double wint ){
    // cout<< "twoBodyDecay( lv, M_d1=" << M_d1 << ", M_d2=" << M_d2 << ")" << endl;

    double dm2 = pow(M_d1 - M_d2,2);
    double sm2 = pow(M_d1 + M_d2,2);

    double M_p = _parent_lv.M();
    double M_p2 = pow(_parent_lv.M(), 2);

    double E1 = (M_p2 + pow(M_d1,2) - pow(M_d2,2)) / ( 2 * M_p );
    double E2 = M_p - E1;

    // cout << "E1 = " << E1 << endl;
    // cout << "E2 = " << E2 << endl;

    double p1 = sqrt( pow(E1,2) - pow(M_d1,2) );
    double p2 = sqrt( pow(E2,2) - pow(M_d2,2) );

    

    double p = sqrt( (M_p2 - sm2) * (M_p2 - dm2) ) / ( 2. * _parent_lv.M() );

    // cout << "p1 = " << p1 << " and p2 = " << p2 << " vs. p = " << p << endl;


    double theta = TMath::ACos( gRandom->Uniform(-1, 1) );// gRandom->Uniform(0,TMath::Pi()/2.0);
    // costheta = gRandom->Uniform(0,1.);
    while ( 2*cos(2*theta) < wint ) {
        theta = gRandom->Uniform(0,TMath::Pi()/2.0);
    }
    double costheta = cos(theta);

    double phi = gRandom->Uniform(0,TMath::Pi()*2);

    double pz = p1*costheta;
    double px = p1*sqrt(1.-costheta*costheta)*cos(phi);
    double py = p1*sqrt(1.-costheta*costheta)*sin(phi);

    TLorentzVector daughter1( px, py, pz, sqrt( p1*p1 + M_d1*M_d1 ));
    TLorentzVector daughter2( -px, -py, -pz, sqrt( p2*p2 + M_d2*M_d2 ) );

    auto v = _parent_lv.Vect().Unit();
    daughter1.RotateUz( v );
    daughter2.RotateUz( v );

    applyBoost(_parent_lv,daughter1);
    applyBoost(_parent_lv,daughter2);

    vector<TLorentzVector> dlv;
    dlv.push_back( daughter1 );
    dlv.push_back( daughter2 );
    return dlv;
}

// compute deltaPhi from two input vectors
double deltaPhi( TLorentzVector lv1, TLorentzVector lv2 ) {
    TLorentzVector lvp = lv1 + lv2;
    TLorentzVector lvn = lv1 - lv2;

    if (gRandom->Rndm() < 0.5)
        lvn = lv2 - lv1;

    return lvn.DeltaPhi( lvp );
}

void analyzeDecay( double ipol, TLorentzVector lv ) {

    double ctheta = ipol/2.0;
    // if (gRandom->Rndm() > 0.5)
    //     ctheta *= -1;

    auto pions = twoBodyDecay( lv, 0.139, 0.139, ctheta );

    hDeltaPhiPt->Fill( lv.Pt(), deltaPhi( pions[0], pions[1] ) );

}




double PomeronIF( double *x, double *par ){
	double v = fFFQ->Eval( x[0] );
	double iff = x[0] / 0.4;
	if ( x[0] > 0.4 ) iff = 1.0;
	return v * ( iff );
	// TVector3 p(0,0,0);
	// p.SetPtEtaPhi( v, 0, gRandom->Rndm() * TMath::Pi() * 2 );
	
}

double PhotonKt( double *x, double *par ){
	double v = fgkt->Eval( x[0] / 1.0 );
	return v;
}


void pair_pt(){


	TFile * f = new TFile( "SnapshotFormFactor4.root" );

	fFFQ = (TF1*)f->Get( "fFFQ" );
	fgkt = (TF1*)f->Get( "fPhotonkt" );

	gRandom = new TRandom3();
	gRandom->SetSeed(0);

	TCanvas *can = new TCanvas( "c", "c", 1200, 900 );
	gStyle->SetOptStat(0);

	TFile * fOut = new TFile( "data.root", "RECREATE" );

	TVector3 p, g, v;

	TH1 * hPairPt = new TH1F( "hPairPt", "", 250, 0, 0.24 );
	TH1 * hPairPtVS = new TH1D( "hPairPtVS", "; p_{T} (GeV/c)", 100, 0, 0.24 );
	TH1 * hPairPtIF = new TH1D( "hPairPtIF", "; p_{T} (GeV/c)", 100, 0, 0.24 );
	TH1 * hPairT = new TH1D( "hPairT", "", 625, 0, 0.0625 );

	TH2 * hgp = new TH2F( "hgp", ";g;p", 100, 0, 0.24, 100, 0, 0.24 );
	TH1 * hdphi = new TH1F( "hdphi", ";#Delta #phi", 200, -TMath::Pi(), TMath::Pi() );
	TH2 * hdPhiPt = new TH2F( "hdPhiPt", ";p_{T} (GeV/c); #Delta #phi", 100, 0, 0.24, 200, -TMath::Pi(), TMath::Pi() );
	TH2 * hdPhiPtIF = new TH2F( "hdPhiPtIF", ";p_{T} (GeV/c); #Delta #phi", 100, 0, 0.24, 200, -TMath::Pi(), TMath::Pi() );
	TH2 * hdPomPhiPt = new TH2F( "hdPomPhiPt", ";p_{T} (GeV/c); #Delta #phi", 100, 0, 0.24, 200, -TMath::Pi(), TMath::Pi() );

	TH2 * hCosdPhiPt = new TH2F( "hCosdPhiPt", ";p_{T} (GeV/c); cos(#Delta #phi( VM, #gamma ))", 100, 0, 0.24, 200, -1.0, 1.0 );
	TH2 * hCosdPhiPtIF = new TH2F( "hCosdPhiPtIF", ";p_{T} (GeV/c); cos(#Delta #phi( VM, #gamma ))", 100, 0, 0.24, 200, -1.0, 1.0 );
	TH2 * hCosdPhiPtIFAbs = new TH2F( "hCosdPhiPtIFAbs", ";p_{T} (GeV/c); cos(#Delta #phi( VM, #gamma ))", 100, 0, 0.24, 200, 0.0, 1.0 );


	TH2 * hCosdPhiPtIso = new TH2F( "hCosdPhiPtIso", ";p_{T} (GeV/c); cos(#Delta #phi( VM, #gamma ))", 100, 0, 0.24, 200, -1.0, 1.0 );
	TH2 * hCosdPhiPtIsoIF = new TH2F( "hCosdPhiPtIsoIF", ";p_{T} (GeV/c); cos(#Delta #phi( VM, #gamma ))", 100, 0, 0.24, 200, -1.0, 1.0 );
	TH2 * hCosdPhiPtIsoIFAbs = new TH2F( "hCosdPhiPtIsoIFAbs", ";p_{T} (GeV/c); cos(#Delta #phi( VM, #gamma ))", 100, 0, 0.24, 200, 0.0, 1.0 );

	TH2 * hdPhiPtIso = new TH2F( "hdPhiPtIso", ";p_{T} (GeV/c); #Delta #phi", 100, 0, 0.24, 200, -TMath::Pi(), TMath::Pi() );
	TH2 * hdPhiPtIFIso = new TH2F( "hdPhiPtIFIso", ";p_{T} (GeV/c); #Delta #phi", 100, 0, 0.24, 200, -TMath::Pi(), TMath::Pi() );

	TH2 * hIF = new TH2F( "hIF", "hIF;p{_T} (GeV/c)", 100, 0, 0.24, 400, -2, 2 );
	TH2 * hdPomVsPho = new TH2F( "hdPomVsPho", ";p_{T} (GeV/c); #Delta #phi", 200, -TMath::Pi(), TMath::Pi(), 200, -TMath::Pi(), TMath::Pi() );
	TH2 * hdMesonVsPho = new TH2F( "hdMesonVsPho", ";p_{T} (GeV/c); #Delta #phi", 200, -TMath::Pi(), TMath::Pi(), 200, -TMath::Pi(), TMath::Pi() );
	TH3 * hdMesonVsPhoVsPt = new TH3F( "hdMesonVsPhoVsPt", ";p_{T} (GeV/c); #Delta #phi", 200, -TMath::Pi(), TMath::Pi(), 200, -TMath::Pi(), TMath::Pi(), 100, 0, 0.24 );
	TH2 * hdMesonVsPom = new TH2F( "hdMesonVsPom", ";p_{T} (GeV/c); #Delta #phi", 200, -TMath::Pi(), TMath::Pi(), 200, -TMath::Pi(), TMath::Pi() );
	TH1 * hImpactParameter = new TH1F( "hImpactParameter", ";b (fm)", 115, 0, 115 );

    hDeltaPhiPt = new TH2F( "hDeltaPhiPt", ";p_{T} GeV/c; #Delta #phi[(#pi^{+}+#pi^{-}), (#pi^{+}-#pi^{-})]", 200, 0, 0.200, 100, -3.14, 3.14 );

	TF1 * fPomeronIF = new TF1( "fPomeronIF", &PomeronIF, 0, 1, 0 );
	TF1 * fPhotonkt = new TF1( "fPhotonkt", &PhotonKt, 0, 1, 0 );
	fPhotonkt->SetNpx( 5000 );
	
	// 1/b^4 weight plus an efficiency term for the STAR trigger
	TF1 * fimpact = new TF1( "fimpact", "pow( x, -4 ) * pow(exp(-1/(x-2*6.36)),50)" );
	fimpact->SetRange( 6.38*2, 115 );

	TLorentzVector lvRho;
	// cout << "photon_dir angle = " << photon_dir << endl;
	// generate the pair pT
	for ( int i = 0; i < 1e5; i++ ){
		double pkt = fFFQ->GetRandom();
		// double pkt = fPomeronIF->GetRandom();
		if ( i % (int)1e4 == 0 ) cout << "." << std::flush;
		for ( int j = 0; j < 1e2; j++ ){
			impact_parameter.SetXYZ( fimpact->GetRandom(), 0, 0 );
			hImpactParameter->Fill( impact_parameter.X() );
			double photon_dir = TMath::ASin( 6.38 / impact_parameter.X() );
			
			// double gkt = fgkt->GetRandom();
			// if ( gkt > 0.055 ) continue;
			double gkt = fPhotonkt->GetRandom();

			p.SetPtEtaPhi( pkt, 0, gRandom->Rndm() * TMath::Pi() * 2 );
			// p.SetPtEtaPhi( pkt, 0, -TMath::Pi() + (TMath::Pi()/22.0) * gRandom->Rndm()
			g.SetPtEtaPhi( gkt, 0, gRandom->Rndm() * TMath::Pi() * 2 );
			// g.SetPtEtaPhi( gkt, 0, TMath::Pi() + (photon_dir) * (2*gRandom->Rndm() - 1.0) );
			// cout << "Px = " << g.Px() << ", Py = " << g.Py() << endl;






			v = p + g;
            lvRho.SetPtEtaPhiM( v.Pt(), gRandom->Uniform(-5.,5.), v.Phi(), gRandom->Uniform( 0.600, 0.950 ) );
            if ( fabs(lvRho.Rapidity()) > 1.0 ) continue;
            // cout << "Rap: " << lvRho.Rapidity() << endl;
			// cout << "v dot b = " << v.Dot( impact_parameter ) << ", pT = " << v.Pt() << endl;

			if ( v.Pt() > 0.1 && v.Pt() < 0.15 ){
				hgp->Fill( gkt, pkt );
			}

			double iff = (1 - cos( v.Dot(impact_parameter) ));

            analyzeDecay( iff, lvRho );



			hdphi->Fill( v.DeltaPhi( g )  );
			hdPhiPt->Fill( v.Pt(), v.DeltaPhi( g ) );
			hdPhiPtIF->Fill( v.Pt(), v.DeltaPhi( g ), iff );

			hCosdPhiPt->Fill( v.Pt(), cos(v.DeltaPhi( g )) );
			hCosdPhiPtIF->Fill( v.Pt(), cos(v.DeltaPhi( g )), iff );
			hCosdPhiPtIFAbs->Fill( v.Pt(), fabs(cos(v.DeltaPhi( g ))), iff );

			hdPomPhiPt->Fill( v.Pt(), v.DeltaPhi( p ) );

			hdMesonVsPom->Fill( v.Phi(), p.Phi() );
			hdMesonVsPho->Fill( v.Phi(), g.Phi() );
			hdMesonVsPhoVsPt->Fill( v.Phi(), g.Phi(), v.Pt() );
			hdPomVsPho->Fill( p.Phi(), g.Phi() );


			

			hPairPt -> Fill( pkt + gkt );
			hPairPtVS -> Fill( v.Pt() );
			hPairPtIF -> Fill( v.Pt(), iff );
			hIF->Fill( v.Pt(), iff );
			hPairT -> Fill( v.Pt()*v.Pt() );


			// isotropic gamma
			// g.SetPtEtaPhi( gkt, 0, gRandom->Rndm() * TMath::Pi() * 2 );
			// v = p + g;
			// iff = (1 - cos( v.Dot(impact_parameter) ));

			// hdPhiPtIso->Fill( v.Pt(), v.DeltaPhi( g ) );
			// hdPhiPtIFIso->Fill( v.Pt(), v.DeltaPhi( g ), iff );

			// hCosdPhiPtIso->Fill( v.Pt(), cos(v.DeltaPhi( g )) );
			// hCosdPhiPtIsoIF->Fill( v.Pt(), cos(v.DeltaPhi( g )), iff );
			// hCosdPhiPtIsoIFAbs->Fill( v.Pt(), fabs(cos(v.DeltaPhi( g ))), iff );

		}
	}

	hPairPt->Scale( 1e-1 / hPairPt->Integral( ), "width" );
	hPairPtVS->Scale( 1.0 / hPairPtVS->GetBinContent( 6 ) );
	// hPairPtIF->Scale( fFFQ->Eval( hPairPtIF->GetBinCenter(26) ) / hPairPtIF->GetBinContent( 26 ) );
	hPairPtIF->Scale( hPairPtVS->GetBinContent( 42 ) / hPairPtIF->GetBinContent( 42 ) );
	hPairPtVS->SetLineColor(kBlack);


	// hPairPt->Draw();
	hPairPtVS->SetMinimum( 4e-5 );
	hPairPtVS->Draw( "hist c" );
	hPairPtVS->SetLineWidth(2);

	hPairPtIF->Draw( "same hist c" );
	hPairPtIF->SetLineStyle( 7 );
	hPairPtIF->SetLineWidth(2);
	hPairPtIF->SetLineColor( kBlack);

	fFFQ->SetLineStyle( 9 );
	fFFQ->SetLineWidth(2);
	fFFQ->Draw("same");

	fPhotonkt->SetLineStyle( 2 );
	fPhotonkt->SetLineColor( kBlue );
	fPhotonkt->SetLineWidth(2);
	
	fPhotonkt->Draw("same");

	TLegend * leg = new TLegend( 0.6, 0.6, 0.89, 0.89 );
	leg->SetBorderSize(0);
	leg->AddEntry( fPhotonkt, "Photon", "l" );
	leg->AddEntry( fFFQ, "Pomeron", "l" );
	leg->AddEntry( hPairPtVS, "Vector Meson", "l" );
	leg->AddEntry( hPairPtIF, "w/Interference", "l" );
	leg->Draw();

	gPad->SetLogy(1);

	can->Print( "transverse_momenta.pdf" );
	can->Print( "transverse_momenta.png" );

    gPad->SetLogy(0);
    hDeltaPhiPt->Draw("colz");

	// make a profile of the cosphi
	hCosdPhiPt->ProfileX( "ProfileXhCosdPhiPt" );
	hCosdPhiPtIF->ProfileX( "ProfileXhCosdPhiPtIF" );
	hCosdPhiPtIFAbs->ProfileX( "ProfileXhCosdPhiPtIFAbs" );

	hCosdPhiPtIso->ProfileX( "ProfileXhCosdPhiPtIso" );
	hCosdPhiPtIsoIF->ProfileX( "ProfileXhCosdPhiPtIsoIF" );
	hCosdPhiPtIsoIFAbs->ProfileX( "ProfileXhCosdPhiPtIsoIFAbs" );


	fOut->Write();
	
}