Page MenuHomeHEPForge
Files F8308725

No OneTemporary
Actions

View Options

Index: trunk/macros/tools/runCore.C
===================================================================
--- trunk/macros/tools/runCore.C (revision 0)
+++ trunk/macros/tools/runCore.C (revision 55)
@@ -0,0 +1,205 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "../runglauber.C"
+#endif
+
+void runQPS(const Int_t n,
+ const char *sysA,
+ const char *sysB,
+ const Double_t signn,
+ const Double_t mind,
+ const Int_t nquarks,
+ Double_t sigqq,
+ const char *fname)
+{
+ // Run Glauber and store ntuple with smeared eccentricities in file.
+
+ TGlauberMC *mcg = new TGlauberMC(sysA,sysB,signn);
+ mcg->SetMinDistance(mind);
+ TFile *out = TFile::Open(fname,"recreate",fname,9);
+ if (!out)
+ return;
+ TNtuple *nt = new TNtuple("nt","nt",
+ "Npart:Ncoll:B:Npqq:Ncqq:Ecc1P:Ecc2P:Ecc3P:Ecc4P:Ecc5P:Ecc1Q:Ecc2Q:Ecc3Q:Ecc4Q:Ecc5Q");
+ nt->SetDirectory(out);
+
+ const TGlauNucleus *nucA = mcg->GetNucleusA();
+ const Int_t AN = nucA->GetN();
+ const Int_t nAN = AN*nquarks;
+ Double_t *xA = new Double_t[nAN];
+ Double_t *yA = new Double_t[nAN];
+ Double_t *zA = new Double_t[nAN];
+ Int_t *nA = new Int_t[nAN];
+ const TGlauNucleus *nucB = mcg-> GetNucleusB();
+ const Int_t BN = nucB->GetN();
+ const Int_t nBN = BN*nquarks;
+ Double_t *xB = new Double_t[nBN];
+ Double_t *yB = new Double_t[nBN];
+ Double_t *zB = new Double_t[nBN];
+ Int_t *nB = new Int_t[nBN];
+
+ TF1 *rad = new TF1("rad","x*x*TMath::Exp(-[0]*x)",0,3);
+ rad->SetParameter(0,4.27);
+ Double_t d2 = (Double_t)sigqq/(TMath::Pi()*10);
+
+ Int_t totalEvs = 0;
+ Int_t accEvs = 0;
+ while (accEvs<n) {
+ if (1)
+ mcg->NextEvent();
+ else
+ while (!mcg->NextEvent()) {;}
+
+ const TObjArray *nucleonsA = nucA->GetNucleons();
+ for (Int_t i = 0,q=0; i<AN; ++i) {
+ TGlauNucleon *nucleonA=(TGlauNucleon*)(nucleonsA->At(i));
+ for (Int_t j=0; j<nquarks; ++j) {
+ Double_t sr = rad->GetRandom();
+ Double_t sp = gRandom->Uniform(-TMath::Pi(), +TMath::Pi());
+ Double_t ctheta = 2*gRandom->Rndm() - 1 ;
+ Double_t stheta = TMath::Sqrt(1-ctheta*ctheta);
+ xA[q] = nucleonA->GetX() + sr*stheta*TMath::Cos(sp);
+ yA[q] = nucleonA->GetY() + sr*stheta*TMath::Sin(sp);
+ zA[q] = nucleonA->GetZ() + sr*ctheta;
+ nA[q] = 0;
+ ++q;
+ }
+ }
+ const TObjArray *nucleonsB = nucB->GetNucleons();
+ for (Int_t i = 0,q=0; i<BN; ++i) {
+ TGlauNucleon *nucleonB=(TGlauNucleon*)(nucleonsB->At(i));
+ for (Int_t j=0; j<nquarks; ++j) {
+ Double_t sr = rad->GetRandom();
+ Double_t sp = gRandom->Uniform(-TMath::Pi(), +TMath::Pi());
+ Double_t ctheta = 2*gRandom->Rndm() - 1 ;
+ Double_t stheta = TMath::Sqrt(1-ctheta*ctheta);
+ xB[q] = nucleonB->GetX() + sr*stheta*TMath::Cos(sp);
+ yB[q] = nucleonB->GetY() + sr*stheta*TMath::Sin(sp);
+ zB[q] = nucleonB->GetZ() + sr*ctheta;
+ nB[q] = 0;
+ ++q;
+ }
+ }
+
+ Int_t ncollqq=0;
+ for (Int_t i = 0; i<nBN; ++i) {
+ for (Int_t j = 0; j<nAN; ++j) {
+ Double_t dx = xB[i]-xA[j];
+ Double_t dy = yB[i]-yA[j];
+ Double_t dij = dx*dx+dy*dy;
+ if (dij < d2) {
+ ++nA[j];
+ ++nB[i];
+ ++ncollqq;
+ }
+ }
+ }
+
+ Int_t npartqq = 0;
+ Double_t xqqvals[99999] = {0};
+ Double_t yqqvals[99999] = {0};
+ for (Int_t i = 0; i<nAN; ++i) {
+ if (nA[i]==0)
+ continue;
+ xqqvals[npartqq] = xA[i];
+ yqqvals[npartqq] = yA[i];
+ ++npartqq;
+
+ }
+ for (Int_t i = 0; i<nBN; ++i) {
+ if (nB[i]==0)
+ continue;
+ xqqvals[npartqq] = xB[i];
+ yqqvals[npartqq] = yB[i];
+ ++npartqq;
+ }
+
+ // *** event counting ***
+ ++totalEvs;
+ if (npartqq==0)
+ continue;
+ ++accEvs;
+ if (!(accEvs%50)) {
+ Double_t bm = mcg->GetBMax();
+ Double_t totxs = (1.*accEvs/totalEvs)*TMath::Pi()*bm*bm/100;
+ Double_t totes = totxs/TMath::Sqrt(accEvs*TMath::Sqrt(1.-accEvs/totalEvs));
+ cout << "Event # " << accEvs << " q-sect = " << totxs << " +- " << totes << " b \r" << flush;
+ }
+ // **********************
+
+ Double_t MeanXQQParts = 0;
+ Double_t MeanYQQParts = 0;
+ Double_t MeanXQQParts2 = 0;
+ Double_t MeanYQQParts2 = 0;
+ Double_t sx2qq = 0;
+ Double_t sy2qq = 0;
+
+ for (Int_t i = 0; i<npartqq; ++i) {
+ MeanXQQParts += xqqvals[i];
+ MeanYQQParts += yqqvals[i];
+ MeanXQQParts2 += xqqvals[i]*xqqvals[i];
+ MeanYQQParts2 += yqqvals[i]*yqqvals[i];
+ }
+ MeanXQQParts /= npartqq;
+ MeanYQQParts /= npartqq;
+ MeanXQQParts2 /= npartqq;
+ MeanYQQParts2 /= npartqq;
+ sx2qq += MeanXQQParts2-MeanXQQParts*MeanXQQParts;
+ sy2qq += MeanYQQParts2-MeanYQQParts*MeanYQQParts;
+
+ Double_t sinphi[6] = {0};
+ Double_t cosphi[6] = {0};
+ Double_t rn[6] = {0};
+ Double_t ecc[6] = {0};
+ Double_t psi[6] = {0};
+ for (Int_t j=1; j<5; ++j) {
+ for (Int_t i = 0; i<npartqq; ++i) {
+ Double_t x = xqqvals[i] - MeanXQQParts;
+ Double_t y = yqqvals[i] - MeanYQQParts;
+ Double_t r = TMath::Sqrt(x*x+y*y);
+ Double_t phi = TMath::ATan2(y,x);
+ Double_t w = j;
+ if (j==1)
+ w = 3; // use r^3 weighting for Ecc1/Psi1
+ cosphi[j] += TMath::Power(r,w)*TMath::Cos(j*phi);
+ sinphi[j] += TMath::Power(r,w)*TMath::Sin(j*phi);
+ rn[j] += TMath::Power(r,w);
+ }
+ }
+ for (Int_t j=1; j<5; ++j) {
+ psi[j] = (TMath::ATan2(sinphi[j],cosphi[j]) + TMath::Pi())/j;
+ ecc[j] = TMath::Sqrt(sinphi[j]*sinphi[j] + cosphi[j]*cosphi[j]) / rn[j];
+ }
+
+ Float_t v[15]; Int_t f=0;
+ v[f++] = mcg->GetNpart();
+ v[f++] = mcg->GetNcoll();
+ v[f++] = mcg->GetB();
+ v[f++] = npartqq;
+ v[f++] = ncollqq;
+ v[f++] = mcg->GetEcc(1);
+ v[f++] = mcg->GetEcc(2);
+ v[f++] = mcg->GetEcc(3);
+ v[f++] = mcg->GetEcc(4);
+ v[f++] = mcg->GetEcc(5);
+ v[f++] = ecc[1];
+ v[f++] = ecc[2];
+ v[f++] = ecc[3];
+ v[f++] = ecc[4];
+ v[f++] = ecc[5];
+
+ nt->Fill(v);
+ }
+ out->Write();
+ out->Close();
+ delete out;
+
+ cout << endl;
+ if (1) {
+ Double_t bm = mcg->GetBMax();
+ Double_t totxs = (1.*accEvs/totalEvs)*TMath::Pi()*bm*bm/100;
+ Double_t totes = totxs/TMath::Sqrt(accEvs*TMath::Sqrt(1.-accEvs/totalEvs));
+ cout << " q-sect = " << totxs << " +- " << totes << " b" << endl;
+ cout << " x-sect = " << mcg->GetTotXSect() << " +- " << mcg->GetTotXSectErr() << " b" << endl;
+ }
+}
+

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 12:58 PM (1 d, 1 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4022861
Default Alt Text
(6 KB)

Event Timeline