Page MenuHomeHEPForge
Files F8723619

cent.C
No OneTemporary
Actions

cent.C
View Options

#include <TH2.h>
#include <TCanvas.h>
#include <TRandom.h>
#include <TMath.h>
#include <TSpline.h>
#include <TF1.h>
#include <TFile.h>
#include <TNtuple.h>
#include <TKey.h>
Double_t NBD(Int_t n, Double_t mu, Double_t k);
TH1F *NBDhist(Double_t mu, Double_t k);
void cent(const char *fname, Double_t k=1.5, Double_t mu=46.6, Double_t f=0.801, Int_t lpe=1)
{
Bool_t draw=0;
TFile *file=TFile::Open(fname);
if (!file)
return;
TIter next(file->GetListOfKeys());
TKey *key;
while ((key=(TKey*)next())) {
TString className(key->GetClassName());
TString keyName(key->GetName());
if (className=="TNtuple")
break;
}
TNtuple *nt = dynamic_cast<TNtuple*>(key->ReadObj());
if (!nt) {
file->ls();
return;
}
cout << nt->GetTitle() << endl;
//PbPb at 5 TeV
Double_t paramf[3] = {1.5, 46.4, 0.801};
paramf[0]=k; paramf[1]=mu; paramf[2]=f;
Int_t nEvents = nt->GetEntries();
cout << "Number of events: " << nEvents << endl;
Float_t impp = 0;
nt->SetBranchAddress("B", &impp);
Float_t Npart = 0;
nt->SetBranchAddress("Npart", &Npart);
Float_t Ncoll = 0;
nt->SetBranchAddress("Ncoll", &Ncoll);
TH1F *nbdh = NBDhist(paramf[1], paramf[0]);
TH1D* hMultV0 = new TH1D("hMultV0", "; Multiplicity V0; Counts", 40000, 0, 40000);
TH1D* hImp = new TH1D("hImpact", "; Impact parameter; Counts", 2000, 0, 20);
for (Int_t iev = 0; iev < nEvents; ++iev) {
nt->GetEntry(iev);
hImp->Fill(impp);
for (Int_t j=0;j<lpe;j++) {
Double_t mult = 0;
Double_t alpha=paramf[2];
Int_t n = (Int_t)(alpha *Npart + (1-alpha) * Ncoll);
for(Int_t m = 0; m<n; ++m) {
mult += (Int_t)nbdh->GetRandom();
}
hMultV0->Fill(mult);
}
}
if (draw) {
TCanvas* cMultV0 = new TCanvas("cMultV0", "cMultV0");
cMultV0->cd();
hMultV0->Draw();
}
Int_t nBm = hMultV0->GetNbinsX();
TH1D* hMultV0Int = new TH1D("hMultV0Int", "; Multiplicity V0; Integral", nBm, 0, 40000);
Double_t intMultV0 = hMultV0->Integral(0, nBm);
Double_t sumMultV0 = 0;
for (Int_t ibQ = 0; ibQ <= nBm; ibQ++){
sumMultV0 += hMultV0->GetBinContent(ibQ);
hMultV0Int->SetBinContent(ibQ, sumMultV0/intMultV0);
}
TSpline* sp3MultV0 = new TSpline3(hMultV0Int);
sp3MultV0->SetName("spMultV0");
if (draw) {
TCanvas* cSp = new TCanvas("cSp", "cSp");
cSp->cd();
sp3MultV0->Draw();
}
Int_t nGm = hImp->GetNbinsX();
TH1D* hImpInt = new TH1D("hImpInt", "; Impact parameter; Integral", nGm, 0, 20.);
Double_t intImp = hImp->Integral(0, nGm);
Double_t sumImp = 0;
for (Int_t ibQ = 0; ibQ <= nGm; ibQ++){
sumImp += hImp->GetBinContent(ibQ);
hImpInt->SetBinContent(ibQ, sumImp/intImp);
}
TSpline* simp = new TSpline3(hImpInt);
simp->SetName("sImp");
if (draw) {
TCanvas* cImp = new TCanvas("cImp", "cImp");
cImp->cd();
simp->Draw();
}
Float_t bcen;
TBranch *bpc = nt->Branch("bcen",&bcen,"bcen/F");
Int_t v0mult;
TBranch *v0m = nt->Branch("v0m",&v0mult,"v0m/I");
Float_t mcen;
TBranch *mpc = nt->Branch("mcen",&mcen,"mcen/F");
for (Int_t iev = 0; iev < nEvents; ++iev) {
nt->GetEntry(iev);
Double_t mult = 0;
Double_t alpha=paramf[2];
Int_t n = (Int_t)(alpha *Npart + (1-alpha) * Ncoll);
for(Int_t m = 0; m<n; ++m) {
mult += (Int_t)nbdh->GetRandom();
}
v0mult=mult;
mcen = 100*(1-sp3MultV0->Eval(mult));
if (mcen<0)
mcen=0;
bcen = 100.*(simp->Eval(impp));
if (bcen>100)
bcen=100;
v0m->Fill();
bpc->Fill();
mpc->Fill();
//cout << bcen << " " << impp << " " << mcen << " " << mult << endl;
}
nt->SetDirectory(0);
TFile *outFile = TFile::Open("out.root","recreate");
nt->Write();
outFile->Close();
}
Double_t NBD(Int_t n, Double_t mu, Double_t k)
{
// Compute NBD.
Double_t F;
Double_t f;
if (n+k > 100.0) {
// log method for handling large numbers
F = TMath::LnGamma(n + k)- TMath::LnGamma(n + 1.)- TMath::LnGamma(k);
f = n * TMath::Log(mu/k) - (n + k) * TMath::Log(1.0 + mu/k);
F = F+f;
F = TMath::Exp(F);
} else {
F = TMath::Gamma(n + k) / ( TMath::Gamma(n + 1.) * TMath::Gamma(k) );
f = n * TMath::Log(mu/k) - (n + k) * TMath::Log(1.0 + mu/k);
f = TMath::Exp(f);
F *= f;
}
return F;
}
TH1F *NBDhist(Double_t mu, Double_t k)
{
// Interface for TH1F.
TH1F *h = new TH1F("hnbd","",100,-0.5,299.5);
h->SetName(Form("nbd_%f_%f",mu,k));
h->SetDirectory(0);
for (Int_t i=0;i<300;++i) {
Double_t val = NBD(i,mu,k);
if (val>1e-20) h->Fill(i,val);
}
return h;
}

File Metadata

Mime Type
text/x-c
Expires
Mon, Jan 20, 8:57 PM (23 h, 7 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242357
Default Alt Text
cent.C (4 KB)

Event Timeline