Page MenuHomeHEPForge
Files F8310040

No OneTemporary
Actions

View Options

Index: trunk/macros/runAndCalcAbsorptionX.C
===================================================================
--- trunk/macros/runAndCalcAbsorptionX.C (revision 31)
+++ trunk/macros/runAndCalcAbsorptionX.C (revision 32)
@@ -1,131 +1,182 @@
#if !defined(__CINT__) || defined(__MAKECINT__)
#include "runglauber.C"
#endif
//---------------------------------------------------------------------------------
void runAndCalcAbsorptionX(const Int_t n,
const Double_t sigX = 1,
const char *sysA = "Pb",
const char *sysB = "Pb",
const Double_t signn = 64,
const Double_t mind = 0.4,
const Int_t nqA = -1,
const Int_t nqB = -1,
const char *fname = "glau_abs_ntuple.root");
//---------------------------------------------------------------------------------
void runAndCalcAbsorptionX(const Int_t n,
const Double_t sigX,
const char *sysA,
const char *sysB,
const Double_t signn,
const Double_t mind,
const Int_t nqA,
const Int_t nqB,
const char *fname)
{
// Run Glauber and store ntuple with smeared eccentricities in file.
TGlauberMC *mcg = new TGlauberMC(sysA,sysB,signn);
mcg->SetMinDistance(mind);
if (nqA>0)
mcg->GetNucleusA()->SetN(nqA);
if (nqB>0)
mcg->GetNucleusB()->SetN(nqB);
TFile *out = TFile::Open(fname,"recreate",fname,9);
if (!out)
return;
TNtuple *nt = new TNtuple("nt","nt",
- "Npart:Ncoll:B:Psi2:Ecc2:Psi3:Ecc3:Nxa:Nxb:S");
+ "Npart:Ncoll:B:Psi2:Ecc2:Psi3:Ecc3:Nxa:Nxb:Sabs:Nxp:Smelt");
nt->SetDirectory(out);
for (Int_t ievent=0; ievent<n; ++ievent) {
while (!mcg->NextEvent()) {}
const TGlauNucleus *nucA = mcg->GetNucleusA();
const TObjArray *nucleonsA = nucA->GetNucleons();
const Int_t AN = nucA->GetN();
const TGlauNucleus *nucB = mcg-> GetNucleusB();
const TObjArray *nucleonsB = nucB->GetNucleons();
const Int_t BN = nucB->GetN();
// find production point
TGlauNucleon *nhA=0,*nhB=0;
Double_t hX=-1;
Double_t hY=-1;
Double_t hZ=-1;
// "ball" diameter = distance at which two balls interact
Double_t d2 = (Double_t)signn/(TMath::Pi()*10); // in fm^2
// for each of the A nucleons in nucleus B
- for (Int_t i = 0; i<BN; i++) {
+ for (Int_t i = 0; i<BN; ++i) {
TGlauNucleon *nucleonB=(TGlauNucleon*)(nucleonsB->At(i));
- for (Int_t j = 0; j < AN ;j++) {
+ for (Int_t j = 0; j < AN ; ++j) {
TGlauNucleon *nucleonA=(TGlauNucleon*)(nucleonsA->At(j));
Double_t dx = nucleonB->GetX()-nucleonA->GetX();
Double_t dy = nucleonB->GetY()-nucleonA->GetY();
Double_t dij = dx*dx+dy*dy;
if (dij < d2) {
nhA=nucleonA;
nhB=nucleonB;
hX = (nucleonA->GetX()+nucleonB->GetX())/2;
hY = (nucleonA->GetY()+nucleonB->GetY())/2;
hZ = (nucleonA->GetZ()+nucleonB->GetZ())/2;
break;
}
}
}
Int_t nXa = 0;
Int_t nXb = 0;
Double_t dX2 = (Double_t)sigX/(TMath::Pi()*10); // in fm^2
- for (Int_t i = 0; i<AN; i++) {
+ for (Int_t i = 0; i<AN; ++i) {
TGlauNucleon *nucleonA=(TGlauNucleon*)(nucleonsA->At(i));
- //cout << "A: " << nucleonA->GetX() << " , " << nucleonA->GetY() << " , " << nucleonA->GetZ() << endl;
if (nhA==nucleonA)
continue;
- //if (nucleonA->GetZ()>hZ)
- // continue;
- Double_t dx = nucleonA->GetX()-hX;
- Double_t dy = nucleonA->GetY()-hY;
+ Double_t dx = nucleonA->GetX()-hX;
+ Double_t dy = nucleonA->GetY()-hY;
Double_t dij = dx*dx+dy*dy;
if (dij < dX2) {
- nXa++;
+ ++nXa;
}
}
- for (Int_t i = 0; i<BN; i++) {
+ for (Int_t i = 0; i<BN; ++i) {
TGlauNucleon *nucleonB=(TGlauNucleon*)(nucleonsB->At(i));
- //cout << "B: " << nucleonB->GetX() << " , " << nucleonB->GetY() << " , " << nucleonB->GetZ() << endl;
if (nhB==nucleonB)
continue;
- //if (nucleonB->GetZ()<hZ)
- // continue;
- Double_t dx = nucleonB->GetX()-hX;
- Double_t dy = nucleonB->GetY()-hY;
+ Double_t dx = nucleonB->GetX()-hX;
+ Double_t dy = nucleonB->GetY()-hY;
Double_t dij = dx*dx+dy*dy;
if (dij < dX2) {
- nXb++;
+ ++nXb;
}
}
- Float_t v[10];
+ Int_t npx = 0;
+ Double_t x0 = hX;
+ Double_t y0 = hY;
+ Double_t phi0 = gRandom->Uniform(0,TMath::TwoPi());
+ Double_t kcphi0 = TMath::Cos(phi0);
+ Double_t ksphi0 = TMath::Sin(phi0);
+ Double_t ksize2 = -1;
+ Double_t kDL = dX2/2.;
+
+ TObjArray parts(1000);
+ parts.SetOwner(0);
+ for (Int_t i = 0; i<AN; ++i) {
+ TGlauNucleon *nucleonA=(TGlauNucleon*)(nucleonsA->At(i));
+ if (nhA==nucleonA)
+ continue;
+ if (!nucleonA->IsWounded())
+ continue;
+ parts.Add(nucleonA);
+ Double_t x = nucleonA->GetX();
+ Double_t y = nucleonA->GetY();
+ Double_t d2 = x*x+y*y;
+ if (d2>ksize2)
+ ksize2 = d2*1.1;
+ }
+ for (Int_t j = 0; j<BN; ++j) {
+ TGlauNucleon *nucleonB=(TGlauNucleon*)(nucleonsB->At(j));
+ if (nhB==nucleonB)
+ continue;
+ if (!nucleonB->IsWounded())
+ continue;
+ parts.Add(nucleonB);
+ Double_t x = nucleonB->GetX();
+ Double_t y = nucleonB->GetY();
+ Double_t d2 = x*x+y*y;
+ if (d2>ksize2)
+ ksize2 = d2*1.1;
+ }
+
+ while(x0*x0+y0*y0<ksize2) {
+ x0 += kDL * kcphi0;
+ y0 += kDL * ksphi0;
+ for (Int_t i = 0; i<parts.GetEntries(); ++i) {
+ TGlauNucleon *part=(TGlauNucleon*)(parts.At(i));
+ if (!part)
+ continue;
+ Double_t dx = part->GetX()-x0;
+ Double_t dy = part->GetY()-y0;
+ Double_t dij = dx*dx+dy*dy;
+ if (dij < dX2) {
+ ++npx;
+ parts.Remove(part);
+ }
+ }
+ }
+
+ Float_t v[12];
v[0] = mcg->GetNpart();
v[1] = mcg->GetNcoll();
v[2] = mcg->GetB();
v[3] = mcg->GetPsi(2);
v[4] = mcg->GetEcc(2);
v[5] = mcg->GetPsi(3);
v[6] = mcg->GetEcc(3);
v[7] = nXa;
v[8] = nXb;
v[9] = TMath::Exp(-nXa-nXb);
+ v[10] = npx;
+ v[11] = TMath::Exp(-npx);
nt->Fill(v);
if (!(ievent%1000))
cout << "Event # " << ievent << " x-sect = " << mcg->GetTotXSect() << " +- " << mcg->GetTotXSectErr() << " b \r" << flush;
}
cout << endl;
out->Write();
out->Close();
delete out;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 5:13 PM (12 h, 31 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023593
Default Alt Text
(6 KB)

Event Timeline