[u/mrichter/AliRoot.git] / PWGPP / EVCHAR / macros / cent / glauberfit.C

// $Id:$

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <Riostream.h>
#include <TH1D.h>
#include <TFile.h>
#include <TTree.h>
#include <TMath.h>
#include <TParticle.h>
#include <TPDGCode.h>
#include <TDatabasePDG.h>
#include <TRandom3.h>
#include <TChain.h>
#include <TROOT.h>
#include <TH2F.h>
#include <TCanvas.h>
#include <TNtuple.h>
#include <TRandom3.h>
#include <TProfile.h>
#include <TColor.h>
#include <TLegend.h>
#include <TF1.h>
#include <TStyle.h>
#include <TGraphErrors.h>
#include <TTimeStamp.h>
#include <createHijingGlauberTestTree.C>
#endif

TObjArray *gList = 0;
TTree     *gGTree = 0;

// from PPR
const Int_t nclasses = 6;
Double_t bmin[nclasses] = {0,3,6,9,12,15};
Double_t bmax[nclasses] = {3,6,9,12,15,18};
Double_t fxs[nclasses];

// analysis cuts
const int nclassesan = 11;
Double_t bminan[nclassesan];
Double_t bmaxan[nclassesan];
Double_t fxsan[nclassesan];
Double_t fxsan1[nclassesan] = {5,5,10,10,10,10,10,10,10,10,10};
Double_t npminan[nclassesan];
Double_t npmaxan[nclassesan];
Double_t npfxsan[nclassesan];
const char *lan[nclassesan] = 
  {"0-5","5-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80","80-90","90-100"};
Int_t colorcl[nclassesan] = 
  {kYellow-9,kYellow,kOrange-4,kOrange+6,kOrange+8,kRed,kRed+1,kRed+2,kMagenta+3,kBlue+3,kBlue+4};
Double_t npmean[nclassesan];
Double_t nprms[nclassesan];

TCanvas *Canvas(const char *name, const char *title=0, Int_t ww=600, Int_t wh=400);
void Classes(TH1 *h, Double_t *resmin, Double_t *resmax, Double_t *fxs, Int_t pos=1, Int_t verbose=0);
TObjArray *Draw(const char *expr, const char *sel=0, const char *opt=0, Int_t type=1,
                const char *scl=0, Double_t *smin=0, Double_t *smax=0, TTree *t=0);
TCanvas *FitNpartDists(const char *name, TObjArray *arr, Int_t verbose=0);
TH1 *Hist(const char *hname, const char *name, const char *title, const char *xtitle, 
          const char *ytitle, Bool_t stats=0, Int_t lc=0, Int_t ls=0);
TH1 *Hist(TH1 *h, const char *name, const char *title, const char *xtitle, 
          const char *ytitle, Bool_t stats=0, Int_t lc=0, Int_t ls=0);
TH1 *Hist(const char *name, const char *title, const char *xtitle, const char *ytitle, 
          Bool_t stats=0, Int_t lc=1, Int_t ls=1);
Double_t NBD(Int_t n, Double_t mu, Double_t k);
TH1F *NBDhist(Double_t mu, Double_t k);
void Store(TObject *o, const char *name=0, const char *fname="gres");

// macro starts here

void glauber() 
{
  Bool_t pImDist               = 0;
  Bool_t pNpDist               = 0;
  Bool_t pNpDistSelWithImp     = 0;
  Bool_t pNpDistSelWithImpFits = 0;
  Bool_t pMidRecResStudy       = 0;
  Bool_t pGlaubFitStudy        = 1;

  Double_t fwdres=0.00;
  Double_t midres=0.00;

  gStyle->SetOptFit(1);
  gROOT->ForceStyle();
  TH1::SetDefaultSumw2(1);

  if (gList) 
    delete gList;
  gList = new TObjArray;
  gList->SetOwner(1);

  TFile *f = TFile::Open("hj-unquenched.root");
  TTree *t = (TTree*)f->Get("glaubertree");
  if (!t) {
    cerr << " not find glaubertree" <<endl;
    return;
  }
  if (gGTree)
    delete gGTree;
  gGTree = t;

  if (1) {
    TNtuple *nt = new TNtuple("nt","nt","g1:g2:g3");
    nt->SetDirectory(0);
    Int_t nents = t->GetEntries();
    for (Int_t i=0;i<nents;++i) {
      nt->Fill(gRandom->Gaus(),gRandom->Gaus(),gRandom->Gaus());
    }
    t->AddFriend(nt,"nt");
  }

  t->SetAlias("Etmidn","response.fEtch0n");
  t->SetAlias("Etmidp","response.fEtch0p");
  t->SetAlias("Etfwdn","response.fEtch3n+response.fEtch4n");
  t->SetAlias("Etfwdp","response.fEtch3p+response.fEtch4p");
  t->SetAlias("Nmidn","response.fNch0n");
  t->SetAlias("Nmidp","response.fNch0p");
  t->SetAlias("Nfwdn","response.fNch3n+response.fNch4n");
  t->SetAlias("Nfwdp","response.fNch3p+response.fNch4p");
  t->SetAlias("Nmid","(Nmidn+Nmidp)/2.");

  t->SetAlias("Etfwdnres",Form("Etfwdn*(1+%f*nt.g1)",fwdres));
  t->SetAlias("Etfwdpres",Form("Etfwdp*(1+%f*nt.g2)",fwdres));
  t->SetAlias("Nmidrec",Form("Nmid*(1+%f*nt.g3)",midres));
  t->SetAlias("npart","header.fNT+header.fNP");
  t->SetAlias("ncoll","header.fN00+header.fN01+header.fN10+header.fN11");
  t->SetAlias("bb","header.fBB");

  t->SetAlias("tresp","1+npart*0.");
  //t->SetAlias("tresp","1-exp(-npart/2.)");
  //t->SetAlias("trig","1+Nmid*0.");
  t->SetAlias("trig","rndm<tresp&&rndm<tresp&&Etfwdnres>5&&Etfwdpres>5");

  TString name;
  if (1) {
    name="ImpactDistFine";
    t->Draw("bb>>htemp(3000,0,30)","1","goff");
    TH1F *hbb=(TH1F*)Hist(name,"","impact parameter [fm]","counts per bin");

    Double_t totxs = 0;
    for (Int_t i=0;i<nclasses;++i) {
      fxs[i] = 100. * hbb->Integral(hbb->FindBin(bmin[i]), hbb->FindBin(bmax[i])) / hbb->Integral();
      totxs += fxs[i];
      printf("Class PPR %d: %.1f - %.1f -> %.1f\n",i+1,bmin[i],bmax[i],fxs[i]);
    }
    cout << "Total from PPR: " << totxs << endl;

    Classes(hbb, bminan, bmaxan, fxsan, 1, 0);
    totxs = 0;
    for (Int_t i=0;i<nclassesan;++i) {
      printf("Centrality Class %d: %.2ffm - %.2ffm -> %.2f (%s)\n",i+1,bminan[i],bmaxan[i],fxsan[i],lan[i]);
      totxs+=fxsan[i];
    }
    cout << "Total: " << totxs << endl;

    if (pImDist) {
      name="ImpactDist";
      TCanvas *c = Canvas(name);
      t->Draw("bb>>htemp(440,0,22)");
      Hist(name,"","impact parameter [fm]","counts per bin");
      TObjArray *arr=Draw("bb>>htemp(440,0,22)",0,"hist");
      TLegend *l = dynamic_cast<TLegend*>(arr->At(0));
      if (l) {
        l->SetX1(0.15); l->SetX2(0.35); l->Draw();
      }
      Store(c);
    }
  }
  if (1) {
    name="NpartDistFine";
    t->Draw("npart>>htemp(440,0,440)","1","goff");
    TH1F *hnp=(TH1F*)Hist(name,"","number participants","counts per bin");

    Classes(hnp, npminan, npmaxan, npfxsan, -1, 0);
    Double_t totxs = 0;
    for (Int_t i=0;i<nclassesan;++i) {
      printf("Npart Class %d: %.1f - %.1f -> %.1f (%s)\n",i+1,npminan[i],npmaxan[i],npfxsan[i],lan[i]);
      totxs+=npfxsan[i];
    }
    cout << "Total: " << totxs << endl;
    if (pNpDist) {
      name="NpartDist";
      TCanvas *c1=Canvas(name);
      c1->SetLogy(1);
      t->Draw("npart>>htemp(440,0,440)","1","goff");
      TH1 *h=Hist(name,"","number participants","counts per bin");
      h->SetMinimum(1);
      h->Draw();
      /*TObjArray *arr=*/Draw("npart>>htemp(440,0,440)",0,"hist",-2);
      Store(c1);
    }
    if (pNpDistSelWithImp) {
      name="NpartDistsWithImpact";
      TCanvas *c1=Canvas(name);
      c1->SetLogy(1);
      t->Draw("npart>>htemp(440,0,440)","1","goff");
      TH1 *h=Hist(name,"","number participants","counts per bin");
      h->SetMinimum(1);
      h->Draw();
      TObjArray *arr=Draw("npart>>htemp(440,0,440)",0,"hist",-1);
      TGraph *ge1 = new TGraph(nclassesan);
      TGraph *ge2 = new TGraph(nclassesan);
      ge1->SetMarkerSize(1.2);
      ge1->SetMarkerStyle(20);
      ge2->SetMarkerSize(1.2);
      ge2->SetMarkerStyle(20);
      for (Int_t i=1;i<arr->GetEntries();++i) {
        TH1F *h = (TH1F*)arr->At(i);
        Int_t N=nclassesan-i;
        Double_t mean  = h->GetMean();
        Double_t width = h->GetRMS();
        npmean[N] = mean;
        nprms[N]  = width;
        ge1->SetPoint(N,N,mean);
        ge2->SetPoint(N,N,width);
      }
      Store(c1);
      Store(ge1,"gNpartMean");
      Store(ge2,"gNpartRms");
      if (pNpDistSelWithImpFits) {
        TCanvas *c2 = FitNpartDists("NpartDistsWithImpactFits",arr,1);
        Store(c2);
      }
    }
  }
  if (pMidRecResStudy) {
    Int_t resint = 100*midres;
    name=Form("MidRecDistribution_res%d",resint);
    t->Draw("Nmidrec>>htemp(3000,0,3000)","1","goff");
    TH1 *hNm = Hist(name,"","Nch in -0.5<#eta<0.5","counts per bin");
    Double_t resmin[nclassesan];
    Double_t resmax[nclassesan];
    Double_t fxs[nclassesan];
    Classes(hNm,resmin,resmax,fxs,-1,1);
    if (1) {
      TCanvas *c = Canvas(name);
      hNm->SetMinimum(1);
      hNm->Draw();
      Draw("Nmidrec>>htemp(3000,0,3000)",0,"hist",-3,"Nmidrec",resmin,resmax);
      Store(c);
    } else {
      delete hNm;
    }

    name=Form("NpartDistsWithTracks_res%d",resint);
    TCanvas *c1 = Canvas(name);
    t->Draw("npart>>htemp(440,0,440)","1","goff");
    TH1 *hNpart = Hist(name,"","number participants","counts per bin");
    hNpart->SetMinimum(1);
    hNpart->Draw();
    TObjArray *arr=Draw("npart>>htemp(440,0,440)",0,"hist",-3,"Nmidrec",resmin,resmax);
    Store(c1);
    if (1) {
      TCanvas *c2 = FitNpartDists(Form("NpartDistsWithTracksFits_res%d",resint),arr,1);
      Store(c2);
    }
    TGraph *ge1 = new TGraph(nclassesan);
    TGraph *ge2 = new TGraph(nclassesan);
    ge1->SetMarkerSize(1.2);
    ge1->SetMarkerStyle(20);
    ge2->SetMarkerSize(1.2);
    ge2->SetMarkerStyle(20);
    for (Int_t i=1;i<arr->GetEntries();++i) {
      Int_t N=nclassesan-i;
      TH1F *h = (TH1F*)arr->At(i);
      Double_t mean  = h->GetMean();
      Double_t rms = h->GetRMS();
      ge1->SetPoint(N,N,mean/npmean[N]-1);
      ge2->SetPoint(N,N,rms/mean*npmean[N]/nprms[N]-1);
      cout << i << " " << mean << " " << rms << " " << npmean[N] << " " << nprms[N] << endl;
    }
    Store(ge1,Form("gMidrecMean_res%d",resint));
    Store(ge2,Form("gMidrecRms_res%d",resint));
  }

  if (pGlaubFitStudy) {
    name="DataDistribution";
    TCanvas *c = Canvas(name);
    t->Draw("Nmid>>htemp(350,0,3500)","1","groff");
    TH1 *h = Hist(name,"","Nch in -0.5<#eta<0.5","counts per bin");
    h->Scale(1./h->Integral());
    h->Draw();
    Store(h);
    //t->Draw("Etfwdnres+Etfwdpres","trig>0");
    //TH1 *h1=Hist(name,"","sum p_{t} [GeV] in 3<|#eta|<5","counts per bin");

    MyHeader *header = 0;
    TBranch *hbr = t->GetBranch("header");
    hbr->SetAddress(&header);
    Int_t nents = t->GetEntries();
    Int_t col=1;
    TLegend *leg = new TLegend(0.5,0.5,0.9,0.9);

    for (Double_t mu=0.5;mu<5;mu+=0.5) {
      for (Double_t k=0.5;k<5;k+=0.5) {
        TH1F *hSample = NBDhist(mu,k);
        TH1F *h1 = new TH1F(Form("fit_%.2f_%2.f",mu,k),"",350,0,3500);
        h1->SetDirectory(0);
        for (Int_t i=0;i<nents;++i) {
          t->GetEntry(i);
          hbr->GetEntry(i);
          Int_t np = header->fNT+header->fNP;
          Double_t alpha=1.1;
          Int_t n = TMath::Nint(TMath::Power(np,alpha));
          Int_t ntot=0;
          for(Int_t j = 0; j<n; ++j)
            ntot+=hSample->GetRandom();
          h1->Fill(ntot);
        }
        h1->SetLineColor(col++);
        h1->Scale(1./h1->Integral());
        leg->AddEntry(h1,Form("mu=%.2f,k=%.2f",mu,k),"l");
        h1->Draw("same");
        Store(h1);
      }
    }
    leg->Draw();
    Store(c);
  }
  return;

  if (0) {
    name="FwdSumCorr";
    Canvas(name);
    t->Draw("Etfwdn:Etfwdp");
    Hist(name,"","sum p_{t} [GeV] in 3<#eta<5","sum p_{t} [GeV] in -3<#eta<-5");
  }
  if (0) {
    name="FwdNvsNpart";
    Canvas(name);
    t->Draw("Etfwdn:npart","1","prof");
    Hist(name,"","Npart","sum p_{t} [GeV] in -3<#eta<-5");
  }
  if (0) {
    name="FwdPvsNpart";
    Canvas(name);
    t->Draw("Etfwdp:npart","1","prof");
    Hist(name,"","Npart","sum p_{t} [GeV] in 3<#eta<5");
  }
  if (0) {
    name="FwdSumvsNpart";
    Canvas(name);
    t->Draw("Etfwdn+Etfwdp:npart","1","prof");
    Hist(name,"","Npart","sum p_{t} [GeV] in 3<|#eta|<5");
  }
  if (0) {
    name="FwdSumDistribution";
    Canvas(name);
    t->Draw("Etfwdn+Etfwdp","1","");
    TH1 *h1=Hist(name,"","sum p_{t} [GeV] in 3<|#eta|<5","counts per bin");
    name="FwdSumTriggered";
    t->Draw("Etfwdn+Etfwdp","trig>0","same");
    TH1 *h2=Hist(name,"","sum p_{t} [GeV] in 3<|#eta|<5","counts per bin",0,2,1);
    Double_t percent = h2->Integral()/h1->Integral()*100.;
    cout << "Recorded " << percent << "% of tot. cross section" << endl;
  }
  if (0) {
    name="Mid2Distribution";
    Canvas(name);
    t->Draw("Nmidn+Nmidp","1","");
    Hist(name,"","Nch in -1<#eta<1","counts per bin");
    name="Mid2Triggered";
    t->Draw("Nmidn+Nmidp","trig>0","same");
    Hist(name,"","Nch in -1<#eta<1","counts per bin",0,2,1);
  }
  if (0) {
    name="MidDistribution";
    Canvas(name);
    t->Draw("Nmid","1","");
    Hist(name,"","Nch in -0.5<#eta<0.5","counts per bin");
    name="MidTriggered";
    t->Draw("Nmid","trig>0","same");
    Hist(name,"","Nch in -0.5<#eta<0.5","counts per bin",0,2,1);
  }
}

//--------------------------------------------------------------------------------------------------

TCanvas *Canvas(const char *name, const char *title, Int_t ww, Int_t wh)
{
  if (!name)
    return 0;
  TString hname(Form("c%s",name));
  if (!title)
    title = name;
  TCanvas *c = new TCanvas(hname,title,ww,wh);
  return c;
}

//--------------------------------------------------------------------------------------------------

void Classes(TH1 *h, Double_t *resmin, Double_t *resmax, Double_t *fxs, Int_t pos, Int_t verbose)
{
  Int_t bfrom = 0;
  Int_t cbin  = h->GetNbinsX();
  if (pos<0) {
    pos   = -1;
    bfrom = h->GetNbinsX()+1;
    cbin  = 1;
  } else {
    pos = 1;
  }

  Double_t totxs=0;
  for (Int_t i=0,lbin=bfrom,bin=lbin;i<nclassesan;++i) {
    lbin = bin+pos;
    Int_t lxs = 0;
    Double_t norm = h->Integral();
    while (1) {
      bin += pos;
      lxs += h->GetBinContent(bin);
      Double_t pxs = lxs/norm*100;
      Double_t tdiff = (lxs+h->GetBinContent(bin+pos))/norm*100;
      //cout << pos << " " << bin << " " << pxs << " " << tdiff << " " << lbin << endl;
      if ((pxs>1&&TMath::Abs(pxs-fxsan1[i])<=TMath::Abs(tdiff-fxsan1[i]))||(bin==cbin)) {
        if (pos>0) {
          resmin[i] = h->GetBinLowEdge(lbin);
          resmax[i] = h->GetBinLowEdge(bin+1);
        } else {
          resmin[i] = h->GetBinLowEdge(bin);
          resmax[i] = h->GetBinLowEdge(lbin+1);
        }
        fxs[i] = pxs;
        if (verbose)
          printf("Class %d: %.1f - %.1f -> %.1f (%s)\n",i+1,resmin[i],resmax[i],pxs,lan[i]);
        totxs += pxs;
        break;
      }
    }
  }
  if (verbose)
    cout << "Total: " << totxs << endl;
}

//--------------------------------------------------------------------------------------------------

TObjArray *Draw(const char *expr, const char *sel, const char *opt, Int_t type,
                const char *scl, Double_t *smin, Double_t *smax, TTree *t)
{
  if (!t) 
    t = gGTree;
  
  TObjArray *oarr = new TObjArray;
  oarr->SetOwner();

  TLegend *l = new TLegend(0.2,0.4,0.4,0.9);
  l->SetBorderSize(0);
  l->SetFillStyle(0);
  oarr->Add(l);

  TString tsel("1");
  if (sel)
    tsel=sel;
  TString doopt("same");
  if (opt)
    doopt=Form("%s,same",opt);

  Int_t beg=0;
  Int_t end=nclassesan;
  if (type<0) {
    beg=nclassesan-1;
    end=-1;
  }
  Double_t *cmin = bminan;
  Double_t *cmax = bmaxan;
  const char *varname="bb";
  if (TMath::Abs(type)==2) {
    cmin = npminan;
    cmax = npmaxan;
    varname="npart";
  } else if (TMath::Abs(type)>2) {
    cmin    = smin;
    cmax    = smax;
    varname = scl;
  }

  while(beg!=end) {
    Int_t i=beg;
    TString dosel(Form("(%s>%.2f&&%s<%.2f)&&(%s)",varname,cmin[i],varname,cmax[i],tsel.Data()));
    t->Draw(expr,dosel,doopt);
    TH1 *h=dynamic_cast<TH1*>(gROOT->FindObject("htemp"));
    if (h) {
      h->SetName(Form("cl_%s_%s",expr,dosel.Data()));
      h->SetLineColor(colorcl[i]);
      h->SetMarkerColor(colorcl[i]);
      h->SetFillColor(colorcl[i]);
      h->SetFillStyle(1000);
      l->AddEntry(h,Form("%s%%",lan[i]),"f");
      oarr->Add(h);
    } else {
      cerr << "Could not obtain htemp for: " << expr << " " << dosel << " " << doopt << endl;
    }
    if(type<0)
      --beg;
    else
      ++beg;
  }
  return oarr;
}

//--------------------------------------------------------------------------------------------------

TCanvas *FitNpartDists(const char *name, TObjArray *arr, Int_t verbose)
{
  TCanvas *c = Canvas(name,0,1200,900);
  c->Divide(4,3);
  TGraphErrors *ge = new TGraphErrors(nclassesan);
  for (Int_t i=1;i<arr->GetEntries();++i) {
    c->cd(i);
    Int_t N=nclassesan-i;
    TH1F *h = (TH1F*)arr->At(i);
    Double_t mean  = h->GetMean();
    Double_t width = h->GetRMS();
    TH1 *h1 = h->DrawCopy("hist");
    h1->SetName(lan[N]);
    h1->SetTitle(Form("%s%% most central",lan[N]));
    h1->SetStats(1);
    h1->SetXTitle("Npart");
    h1->SetYTitle("counts per bin");
    TF1 *fit = new TF1(Form("fit%d",i),"gaus(0)",0,440);
    fit->SetParameters(1,mean,width);
    fit->SetLineWidth(3);
    h1->Fit(fit,"QM0");
    fit->Draw("same");
    ge->SetPoint(N,mean,width);
    //ge->SetPointError(N,0,width);
    if (verbose) 
      cout << i << " hist: " << mean << " " << width 
           << " fit:  " << fit->GetParameter(1) << " " << fit->GetParameter(2) << endl;
  }
  c->cd(12);
  TH2F *h2f = new TH2F(Form("%s_summary",name),";#LTNpart#GT;RMS",1,0,440,1,0,49);
  h2f->SetStats(0);
  h2f->Draw();
  ge->SetMarkerSize(1.2);
  ge->SetMarkerStyle(20);
  ge->Draw("P");
  return c;
}

//--------------------------------------------------------------------------------------------------

TH1 *Hist(const char *name, const char *title, const char *xtitle, const char *ytitle, 
          Bool_t stats, Int_t lc, Int_t ls)
{
  TH1 *h=dynamic_cast<TH1*>(gROOT->FindObject("htemp"));
  return Hist(h,name,title,xtitle,ytitle,stats,lc,ls);
}

//--------------------------------------------------------------------------------------------------

TH1 *Hist(const char *hname, const char *name, const char *title, const char *xtitle, 
          const char *ytitle, Bool_t stats, Int_t lc, Int_t ls)
{
  TH1 *h=dynamic_cast<TH1*>(gROOT->FindObject(hname));
  return Hist(h,name,title,xtitle,ytitle,stats,lc,ls);
}

//--------------------------------------------------------------------------------------------------

TH1 *Hist(TH1 *h, const char *name, const char *title, const char *xtitle, 
          const char *ytitle, Bool_t stats, Int_t lc, Int_t ls)
{
  if (!h && !name) 
    return 0;
  TString hname(Form("h%s",name));
  h->SetName(hname);
  if (!title)
    title = name;
  h->SetTitle(title);
  if (xtitle)
    h->SetXTitle(xtitle);
  if (ytitle)
    h->SetYTitle(ytitle);
  h->GetXaxis()->SetTitleOffset(1.1);
  h->GetYaxis()->SetTitleOffset(1.2);
  h->SetStats(stats);
  h->SetLineColor(lc);
  h->SetLineStyle(ls);
  gList->Add(h);
  return h;
}

//--------------------------------------------------------------------------------------------------

void Store(TObject *o, const char *name, const char *fname)
{
  if (!o || !fname)
    return;

  TTimeStamp t;
  TString filename(Form("%s_%d.root",fname, t.GetDate()));
  TFile *f = new TFile(filename,"update");
  if (!name)
    name = o->GetName();
  if (o) {
    f->Delete(Form("%s;1",name));
    o->Write(name);
  }

  f->Close();
  delete f;
}

//--------------------------------------------------------------------------------------------------

Double_t NBD(Int_t n, Double_t mu, Double_t k)
{
  Double_t mudk = mu/k;
  Double_t ret = TMath::Gamma(n+k) / TMath::Gamma(k) / TMath::Factorial(n) *
                 TMath::Power(mudk,n) / TMath::Power(1+mudk,n+k);
  return ret;
}

//--------------------------------------------------------------------------------------------------

TH1F *NBDhist(Double_t mu, Double_t k)
{
  TH1F *h = new TH1F("htemp","",100,0,100);
  h->SetName(Form("nbd_%f_%f",mu,k));
  h->SetDirectory(0);
  for (Int_t i=0;i<100;++i) {
    Double_t val = NBD(i,mu,k);
    h->Fill(i,val);
  }
  return h;
}
Commit	Line	Data
8faabf6b	1	// $Id:$
	2
	3	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	4	#include <Riostream.h>
	5	#include <TH1D.h>
	6	#include <TFile.h>
	7	#include <TTree.h>
	8	#include <TMath.h>
	9	#include <TParticle.h>
	10	#include <TPDGCode.h>
	11	#include <TDatabasePDG.h>
	12	#include <TRandom3.h>
	13	#include <TChain.h>
	14	#include <TROOT.h>
	15	#include <TH2F.h>
	16	#include <TCanvas.h>
	17	#include <TNtuple.h>
	18	#include <TRandom3.h>
	19	#include <TProfile.h>
	20	#include <TColor.h>
	21	#include <TLegend.h>
	22	#include <TF1.h>
	23	#include <TStyle.h>
	24	#include <TGraphErrors.h>
	25	#include <TTimeStamp.h>
fc57ea63	26	#include <createHijingGlauberTestTree.C>
8faabf6b	27	#endif
	28
	29	TObjArray *gList = 0;
	30	TTree *gGTree = 0;
	31
	32	// from PPR
	33	const Int_t nclasses = 6;
	34	Double_t bmin[nclasses] = {0,3,6,9,12,15};
	35	Double_t bmax[nclasses] = {3,6,9,12,15,18};
	36	Double_t fxs[nclasses];
	37
	38	// analysis cuts
	39	const int nclassesan = 11;
	40	Double_t bminan[nclassesan];
	41	Double_t bmaxan[nclassesan];
	42	Double_t fxsan[nclassesan];
	43	Double_t fxsan1[nclassesan] = {5,5,10,10,10,10,10,10,10,10,10};
	44	Double_t npminan[nclassesan];
	45	Double_t npmaxan[nclassesan];
	46	Double_t npfxsan[nclassesan];
	47	const char *lan[nclassesan] =
	48	{"0-5","5-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80","80-90","90-100"};
	49	Int_t colorcl[nclassesan] =
	50	{kYellow-9,kYellow,kOrange-4,kOrange+6,kOrange+8,kRed,kRed+1,kRed+2,kMagenta+3,kBlue+3,kBlue+4};
	51	Double_t npmean[nclassesan];
	52	Double_t nprms[nclassesan];
	53
	54	TCanvas Canvas(const char name, const char *title=0, Int_t ww=600, Int_t wh=400);
	55	void Classes(TH1 h, Double_t resmin, Double_t resmax, Double_t fxs, Int_t pos=1, Int_t verbose=0);
	56	TObjArray Draw(const char expr, const char sel=0, const char opt=0, Int_t type=1,
	57	const char scl=0, Double_t smin=0, Double_t smax=0, TTree t=0);
	58	TCanvas FitNpartDists(const char name, TObjArray *arr, Int_t verbose=0);
	59	TH1 Hist(const char hname, const char name, const char title, const char *xtitle,
	60	const char *ytitle, Bool_t stats=0, Int_t lc=0, Int_t ls=0);
	61	TH1 Hist(TH1 h, const char name, const char title, const char *xtitle,
	62	const char *ytitle, Bool_t stats=0, Int_t lc=0, Int_t ls=0);
	63	TH1 Hist(const char name, const char title, const char xtitle, const char *ytitle,
	64	Bool_t stats=0, Int_t lc=1, Int_t ls=1);
	65	Double_t NBD(Int_t n, Double_t mu, Double_t k);
	66	TH1F *NBDhist(Double_t mu, Double_t k);
	67	void Store(TObject o, const char name=0, const char *fname="gres");
	68
	69	// macro starts here
	70
	71	void glauber()
	72	{
	73	Bool_t pImDist = 0;
	74	Bool_t pNpDist = 0;
	75	Bool_t pNpDistSelWithImp = 0;
	76	Bool_t pNpDistSelWithImpFits = 0;
	77	Bool_t pMidRecResStudy = 0;
	78	Bool_t pGlaubFitStudy = 1;
	79
	80	Double_t fwdres=0.00;
	81	Double_t midres=0.00;
	82
	83	gStyle->SetOptFit(1);
	84	gROOT->ForceStyle();
	85	TH1::SetDefaultSumw2(1);
	86
	87	if (gList)
	88	delete gList;
	89	gList = new TObjArray;
	90	gList->SetOwner(1);
91
92	TFile *f = TFile::Open("hj-unquenched.root");
93	TTree t = (TTree)f->Get("glaubertree");
94	if (!t) {
95	cerr << " not find glaubertree" <<endl;
96	return;
97	}
98	if (gGTree)
99	delete gGTree;
100	gGTree = t;
101
102	if (1) {
103	TNtuple *nt = new TNtuple("nt","nt","g1:g2:g3");
104	nt->SetDirectory(0);
105	Int_t nents = t->GetEntries();
106	for (Int_t i=0;i<nents;++i) {
107	nt->Fill(gRandom->Gaus(),gRandom->Gaus(),gRandom->Gaus());
108	}
109	t->AddFriend(nt,"nt");
110	}
111
112	t->SetAlias("Etmidn","response.fEtch0n");
113	t->SetAlias("Etmidp","response.fEtch0p");
114	t->SetAlias("Etfwdn","response.fEtch3n+response.fEtch4n");
115	t->SetAlias("Etfwdp","response.fEtch3p+response.fEtch4p");
116	t->SetAlias("Nmidn","response.fNch0n");
117	t->SetAlias("Nmidp","response.fNch0p");
118	t->SetAlias("Nfwdn","response.fNch3n+response.fNch4n");
119	t->SetAlias("Nfwdp","response.fNch3p+response.fNch4p");
120	t->SetAlias("Nmid","(Nmidn+Nmidp)/2.");
121
122	t->SetAlias("Etfwdnres",Form("Etfwdn(1+%fnt.g1)",fwdres));
123	t->SetAlias("Etfwdpres",Form("Etfwdp(1+%fnt.g2)",fwdres));
124	t->SetAlias("Nmidrec",Form("Nmid(1+%fnt.g3)",midres));
125	t->SetAlias("npart","header.fNT+header.fNP");
126	t->SetAlias("ncoll","header.fN00+header.fN01+header.fN10+header.fN11");
127	t->SetAlias("bb","header.fBB");
128
129	t->SetAlias("tresp","1+npart*0.");
130	//t->SetAlias("tresp","1-exp(-npart/2.)");
131	//t->SetAlias("trig","1+Nmid*0.");
132	t->SetAlias("trig","rndm<tresp&&rndm<tresp&&Etfwdnres>5&&Etfwdpres>5");
133
134	TString name;
135	if (1) {
136	name="ImpactDistFine";
137	t->Draw("bb>>htemp(3000,0,30)","1","goff");
138	TH1F hbb=(TH1F)Hist(name,"","impact parameter [fm]","counts per bin");
139
140	Double_t totxs = 0;
141	for (Int_t i=0;i<nclasses;++i) {
142	fxs[i] = 100. * hbb->Integral(hbb->FindBin(bmin[i]), hbb->FindBin(bmax[i])) / hbb->Integral();
143	totxs += fxs[i];
144	printf("Class PPR %d: %.1f - %.1f -> %.1f\n",i+1,bmin[i],bmax[i],fxs[i]);
145	}
146	cout << "Total from PPR: " << totxs << endl;
147
148	Classes(hbb, bminan, bmaxan, fxsan, 1, 0);
149	totxs = 0;
150	for (Int_t i=0;i<nclassesan;++i) {
151	printf("Centrality Class %d: %.2ffm - %.2ffm -> %.2f (%s)\n",i+1,bminan[i],bmaxan[i],fxsan[i],lan[i]);
152	totxs+=fxsan[i];
153	}
154	cout << "Total: " << totxs << endl;
155
156	if (pImDist) {
157	name="ImpactDist";
158	TCanvas *c = Canvas(name);
159	t->Draw("bb>>htemp(440,0,22)");
160	Hist(name,"","impact parameter [fm]","counts per bin");
161	TObjArray *arr=Draw("bb>>htemp(440,0,22)",0,"hist");
162	TLegend l = dynamic_cast<TLegend>(arr->At(0));
163	if (l) {
164	l->SetX1(0.15); l->SetX2(0.35); l->Draw();
165	}
166	Store(c);
167	}
168	}
169	if (1) {
170	name="NpartDistFine";
171	t->Draw("npart>>htemp(440,0,440)","1","goff");
172	TH1F hnp=(TH1F)Hist(name,"","number participants","counts per bin");
173
174	Classes(hnp, npminan, npmaxan, npfxsan, -1, 0);
175	Double_t totxs = 0;
176	for (Int_t i=0;i<nclassesan;++i) {
177	printf("Npart Class %d: %.1f - %.1f -> %.1f (%s)\n",i+1,npminan[i],npmaxan[i],npfxsan[i],lan[i]);
178	totxs+=npfxsan[i];
179	}
180	cout << "Total: " << totxs << endl;
181	if (pNpDist) {
182	name="NpartDist";
183	TCanvas *c1=Canvas(name);
184	c1->SetLogy(1);
185	t->Draw("npart>>htemp(440,0,440)","1","goff");
186	TH1 *h=Hist(name,"","number participants","counts per bin");
187	h->SetMinimum(1);
188	h->Draw();
189	/TObjArray arr=*/Draw("npart>>htemp(440,0,440)",0,"hist",-2);
190	Store(c1);
191	}
192	if (pNpDistSelWithImp) {
193	name="NpartDistsWithImpact";
194	TCanvas *c1=Canvas(name);
195	c1->SetLogy(1);
196	t->Draw("npart>>htemp(440,0,440)","1","goff");
197	TH1 *h=Hist(name,"","number participants","counts per bin");
198	h->SetMinimum(1);
199	h->Draw();
200	TObjArray *arr=Draw("npart>>htemp(440,0,440)",0,"hist",-1);
201	TGraph *ge1 = new TGraph(nclassesan);
202	TGraph *ge2 = new TGraph(nclassesan);
203	ge1->SetMarkerSize(1.2);
204	ge1->SetMarkerStyle(20);
205	ge2->SetMarkerSize(1.2);
206	ge2->SetMarkerStyle(20);
207	for (Int_t i=1;i<arr->GetEntries();++i) {
208	TH1F h = (TH1F)arr->At(i);
209	Int_t N=nclassesan-i;
210	Double_t mean = h->GetMean();
211	Double_t width = h->GetRMS();
212	npmean[N] = mean;
213	nprms[N] = width;
214	ge1->SetPoint(N,N,mean);
215	ge2->SetPoint(N,N,width);
216	}
217	Store(c1);
218	Store(ge1,"gNpartMean");
219	Store(ge2,"gNpartRms");
220	if (pNpDistSelWithImpFits) {
221	TCanvas *c2 = FitNpartDists("NpartDistsWithImpactFits",arr,1);
222	Store(c2);
223	}
224	}
225	}
226	if (pMidRecResStudy) {
227	Int_t resint = 100*midres;
228	name=Form("MidRecDistribution_res%d",resint);
229	t->Draw("Nmidrec>>htemp(3000,0,3000)","1","goff");
230	TH1 *hNm = Hist(name,"","Nch in -0.5<#eta<0.5","counts per bin");
231	Double_t resmin[nclassesan];
232	Double_t resmax[nclassesan];
233	Double_t fxs[nclassesan];
234	Classes(hNm,resmin,resmax,fxs,-1,1);
235	if (1) {
236	TCanvas *c = Canvas(name);
237	hNm->SetMinimum(1);
238	hNm->Draw();
239	Draw("Nmidrec>>htemp(3000,0,3000)",0,"hist",-3,"Nmidrec",resmin,resmax);
240	Store(c);
241	} else {
242	delete hNm;
243	}
244
245	name=Form("NpartDistsWithTracks_res%d",resint);
246	TCanvas *c1 = Canvas(name);
247	t->Draw("npart>>htemp(440,0,440)","1","goff");
248	TH1 *hNpart = Hist(name,"","number participants","counts per bin");
249	hNpart->SetMinimum(1);
250	hNpart->Draw();
251	TObjArray *arr=Draw("npart>>htemp(440,0,440)",0,"hist",-3,"Nmidrec",resmin,resmax);
252	Store(c1);
253	if (1) {
254	TCanvas *c2 = FitNpartDists(Form("NpartDistsWithTracksFits_res%d",resint),arr,1);
255	Store(c2);
256	}
257	TGraph *ge1 = new TGraph(nclassesan);
258	TGraph *ge2 = new TGraph(nclassesan);
259	ge1->SetMarkerSize(1.2);
260	ge1->SetMarkerStyle(20);
261	ge2->SetMarkerSize(1.2);
262	ge2->SetMarkerStyle(20);
263	for (Int_t i=1;i<arr->GetEntries();++i) {
264	Int_t N=nclassesan-i;
265	TH1F h = (TH1F)arr->At(i);
266	Double_t mean = h->GetMean();
267	Double_t rms = h->GetRMS();
268	ge1->SetPoint(N,N,mean/npmean[N]-1);
269	ge2->SetPoint(N,N,rms/mean*npmean[N]/nprms[N]-1);
270	cout << i << " " << mean << " " << rms << " " << npmean[N] << " " << nprms[N] << endl;
271	}
272	Store(ge1,Form("gMidrecMean_res%d",resint));
273	Store(ge2,Form("gMidrecRms_res%d",resint));
274	}
275
276	if (pGlaubFitStudy) {
277	name="DataDistribution";
278	TCanvas *c = Canvas(name);
279	t->Draw("Nmid>>htemp(350,0,3500)","1","groff");
280	TH1 *h = Hist(name,"","Nch in -0.5<#eta<0.5","counts per bin");
281	h->Scale(1./h->Integral());
282	h->Draw();
283	Store(h);
284	//t->Draw("Etfwdnres+Etfwdpres","trig>0");
285	//TH1 *h1=Hist(name,"","sum p_{t} [GeV] in 3<\|#eta\|<5","counts per bin");
286
287	MyHeader *header = 0;
288	TBranch *hbr = t->GetBranch("header");
289	hbr->SetAddress(&header);
290	Int_t nents = t->GetEntries();
291	Int_t col=1;
292	TLegend *leg = new TLegend(0.5,0.5,0.9,0.9);
293
294	for (Double_t mu=0.5;mu<5;mu+=0.5) {
295	for (Double_t k=0.5;k<5;k+=0.5) {
296	TH1F *hSample = NBDhist(mu,k);
297	TH1F *h1 = new TH1F(Form("fit_%.2f_%2.f",mu,k),"",350,0,3500);
298	h1->SetDirectory(0);
299	for (Int_t i=0;i<nents;++i) {
300	t->GetEntry(i);
301	hbr->GetEntry(i);
302	Int_t np = header->fNT+header->fNP;
303	Double_t alpha=1.1;
304	Int_t n = TMath::Nint(TMath::Power(np,alpha));
305	Int_t ntot=0;
306	for(Int_t j = 0; j<n; ++j)
307	ntot+=hSample->GetRandom();
308	h1->Fill(ntot);
309	}
310	h1->SetLineColor(col++);
311	h1->Scale(1./h1->Integral());
312	leg->AddEntry(h1,Form("mu=%.2f,k=%.2f",mu,k),"l");
313	h1->Draw("same");
314	Store(h1);
315	}
316	}
317	leg->Draw();
318	Store(c);
319	}
320	return;
321
322	if (0) {
323	name="FwdSumCorr";
324	Canvas(name);
325	t->Draw("Etfwdn:Etfwdp");
326	Hist(name,"","sum p_{t} [GeV] in 3<#eta<5","sum p_{t} [GeV] in -3<#eta<-5");
327	}
328	if (0) {
329	name="FwdNvsNpart";
330	Canvas(name);
331	t->Draw("Etfwdn:npart","1","prof");
332	Hist(name,"","Npart","sum p_{t} [GeV] in -3<#eta<-5");
333	}
334	if (0) {
335	name="FwdPvsNpart";
336	Canvas(name);
337	t->Draw("Etfwdp:npart","1","prof");
338	Hist(name,"","Npart","sum p_{t} [GeV] in 3<#eta<5");
339	}
340	if (0) {
341	name="FwdSumvsNpart";
342	Canvas(name);
343	t->Draw("Etfwdn+Etfwdp:npart","1","prof");
344	Hist(name,"","Npart","sum p_{t} [GeV] in 3<\|#eta\|<5");
345	}
346	if (0) {
347	name="FwdSumDistribution";
348	Canvas(name);
349	t->Draw("Etfwdn+Etfwdp","1","");
350	TH1 *h1=Hist(name,"","sum p_{t} [GeV] in 3<\|#eta\|<5","counts per bin");
351	name="FwdSumTriggered";
352	t->Draw("Etfwdn+Etfwdp","trig>0","same");
353	TH1 *h2=Hist(name,"","sum p_{t} [GeV] in 3<\|#eta\|<5","counts per bin",0,2,1);
354	Double_t percent = h2->Integral()/h1->Integral()*100.;
355	cout << "Recorded " << percent << "% of tot. cross section" << endl;
356	}
357	if (0) {
358	name="Mid2Distribution";
359	Canvas(name);
360	t->Draw("Nmidn+Nmidp","1","");
361	Hist(name,"","Nch in -1<#eta<1","counts per bin");
362	name="Mid2Triggered";
363	t->Draw("Nmidn+Nmidp","trig>0","same");
364	Hist(name,"","Nch in -1<#eta<1","counts per bin",0,2,1);
365	}
366	if (0) {
367	name="MidDistribution";
368	Canvas(name);
369	t->Draw("Nmid","1","");
370	Hist(name,"","Nch in -0.5<#eta<0.5","counts per bin");
371	name="MidTriggered";
372	t->Draw("Nmid","trig>0","same");
373	Hist(name,"","Nch in -0.5<#eta<0.5","counts per bin",0,2,1);
374	}
375	}
376
377	//--------------------------------------------------------------------------------------------------
378
379	TCanvas Canvas(const char name, const char *title, Int_t ww, Int_t wh)
380	{
381	if (!name)
382	return 0;
383	TString hname(Form("c%s",name));
384	if (!title)
385	title = name;
386	TCanvas *c = new TCanvas(hname,title,ww,wh);
387	return c;
388	}
389
390	//--------------------------------------------------------------------------------------------------
391
392	void Classes(TH1 h, Double_t resmin, Double_t resmax, Double_t fxs, Int_t pos, Int_t verbose)
393	{
394	Int_t bfrom = 0;
395	Int_t cbin = h->GetNbinsX();
396	if (pos<0) {
397	pos = -1;
398	bfrom = h->GetNbinsX()+1;
399	cbin = 1;
400	} else {
401	pos = 1;
402	}
403
404	Double_t totxs=0;
405	for (Int_t i=0,lbin=bfrom,bin=lbin;i<nclassesan;++i) {
406	lbin = bin+pos;
407	Int_t lxs = 0;
408	Double_t norm = h->Integral();
409	while (1) {
410	bin += pos;
411	lxs += h->GetBinContent(bin);
412	Double_t pxs = lxs/norm*100;
413	Double_t tdiff = (lxs+h->GetBinContent(bin+pos))/norm*100;
414	//cout << pos << " " << bin << " " << pxs << " " << tdiff << " " << lbin << endl;
415	if ((pxs>1&&TMath::Abs(pxs-fxsan1[i])<=TMath::Abs(tdiff-fxsan1[i]))\|\|(bin==cbin)) {
416	if (pos>0) {
417	resmin[i] = h->GetBinLowEdge(lbin);
418	resmax[i] = h->GetBinLowEdge(bin+1);
419	} else {
420	resmin[i] = h->GetBinLowEdge(bin);
421	resmax[i] = h->GetBinLowEdge(lbin+1);
422	}
423	fxs[i] = pxs;
424	if (verbose)
425	printf("Class %d: %.1f - %.1f -> %.1f (%s)\n",i+1,resmin[i],resmax[i],pxs,lan[i]);
426	totxs += pxs;
427	break;
428	}
429	}
430	}
431	if (verbose)
432	cout << "Total: " << totxs << endl;
433	}
434
435	//--------------------------------------------------------------------------------------------------
436
437	TObjArray Draw(const char expr, const char sel, const char opt, Int_t type,
438	const char scl, Double_t smin, Double_t smax, TTree t)
439	{
440	if (!t)
441	t = gGTree;
442
443	TObjArray *oarr = new TObjArray;
444	oarr->SetOwner();
445
446	TLegend *l = new TLegend(0.2,0.4,0.4,0.9);
447	l->SetBorderSize(0);
448	l->SetFillStyle(0);
449	oarr->Add(l);
450
451	TString tsel("1");
452	if (sel)
453	tsel=sel;
454	TString doopt("same");
455	if (opt)
456	doopt=Form("%s,same",opt);
457
458	Int_t beg=0;
459	Int_t end=nclassesan;
460	if (type<0) {
461	beg=nclassesan-1;
462	end=-1;
463	}
464	Double_t *cmin = bminan;
465	Double_t *cmax = bmaxan;
466	const char *varname="bb";
467	if (TMath::Abs(type)==2) {
468	cmin = npminan;
469	cmax = npmaxan;
470	varname="npart";
471	} else if (TMath::Abs(type)>2) {
472	cmin = smin;
473	cmax = smax;
474	varname = scl;
475	}
476
477	while(beg!=end) {
478	Int_t i=beg;
479	TString dosel(Form("(%s>%.2f&&%s<%.2f)&&(%s)",varname,cmin[i],varname,cmax[i],tsel.Data()));
480	t->Draw(expr,dosel,doopt);
481	TH1 h=dynamic_cast<TH1>(gROOT->FindObject("htemp"));
482	if (h) {
483	h->SetName(Form("cl_%s_%s",expr,dosel.Data()));
484	h->SetLineColor(colorcl[i]);
485	h->SetMarkerColor(colorcl[i]);
486	h->SetFillColor(colorcl[i]);
487	h->SetFillStyle(1000);
488	l->AddEntry(h,Form("%s%%",lan[i]),"f");
489	oarr->Add(h);
490	} else {
491	cerr << "Could not obtain htemp for: " << expr << " " << dosel << " " << doopt << endl;
492	}
493	if(type<0)
494	--beg;
495	else
496	++beg;
497	}
498	return oarr;
499	}
500
501	//--------------------------------------------------------------------------------------------------
502
503	TCanvas FitNpartDists(const char name, TObjArray *arr, Int_t verbose)
504	{
505	TCanvas *c = Canvas(name,0,1200,900);
506	c->Divide(4,3);
507	TGraphErrors *ge = new TGraphErrors(nclassesan);
508	for (Int_t i=1;i<arr->GetEntries();++i) {
509	c->cd(i);
510	Int_t N=nclassesan-i;
511	TH1F h = (TH1F)arr->At(i);
512	Double_t mean = h->GetMean();
513	Double_t width = h->GetRMS();
514	TH1 *h1 = h->DrawCopy("hist");
515	h1->SetName(lan[N]);
516	h1->SetTitle(Form("%s%% most central",lan[N]));
517	h1->SetStats(1);
518	h1->SetXTitle("Npart");
519	h1->SetYTitle("counts per bin");
520	TF1 *fit = new TF1(Form("fit%d",i),"gaus(0)",0,440);
521	fit->SetParameters(1,mean,width);
522	fit->SetLineWidth(3);
523	h1->Fit(fit,"QM0");
524	fit->Draw("same");
525	ge->SetPoint(N,mean,width);
526	//ge->SetPointError(N,0,width);
527	if (verbose)
528	cout << i << " hist: " << mean << " " << width
529	<< " fit: " << fit->GetParameter(1) << " " << fit->GetParameter(2) << endl;
530	}
531	c->cd(12);
532	TH2F *h2f = new TH2F(Form("%s_summary",name),";#LTNpart#GT;RMS",1,0,440,1,0,49);
533	h2f->SetStats(0);
534	h2f->Draw();
535	ge->SetMarkerSize(1.2);
536	ge->SetMarkerStyle(20);
537	ge->Draw("P");
538	return c;
539	}
540
541	//--------------------------------------------------------------------------------------------------
542
543	TH1 Hist(const char name, const char title, const char xtitle, const char *ytitle,
544	Bool_t stats, Int_t lc, Int_t ls)
545	{
546	TH1 h=dynamic_cast<TH1>(gROOT->FindObject("htemp"));
547	return Hist(h,name,title,xtitle,ytitle,stats,lc,ls);
548	}
549
550	//--------------------------------------------------------------------------------------------------
551
552	TH1 Hist(const char hname, const char name, const char title, const char *xtitle,
553	const char *ytitle, Bool_t stats, Int_t lc, Int_t ls)
554	{
555	TH1 h=dynamic_cast<TH1>(gROOT->FindObject(hname));
556	return Hist(h,name,title,xtitle,ytitle,stats,lc,ls);
557	}
558
559	//--------------------------------------------------------------------------------------------------
560
561	TH1 Hist(TH1 h, const char name, const char title, const char *xtitle,
562	const char *ytitle, Bool_t stats, Int_t lc, Int_t ls)
563	{
564	if (!h && !name)
565	return 0;
566	TString hname(Form("h%s",name));
567	h->SetName(hname);
568	if (!title)
569	title = name;
570	h->SetTitle(title);
571	if (xtitle)
572	h->SetXTitle(xtitle);
573	if (ytitle)
574	h->SetYTitle(ytitle);
575	h->GetXaxis()->SetTitleOffset(1.1);
576	h->GetYaxis()->SetTitleOffset(1.2);
577	h->SetStats(stats);
578	h->SetLineColor(lc);
579	h->SetLineStyle(ls);
580	gList->Add(h);
581	return h;
582	}
583
584	//--------------------------------------------------------------------------------------------------
585
586	void Store(TObject o, const char name, const char *fname)
587	{
588	if (!o \|\| !fname)
589	return;
590
591	TTimeStamp t;
592	TString filename(Form("%s_%d.root",fname, t.GetDate()));
593	TFile *f = new TFile(filename,"update");
594	if (!name)
595	name = o->GetName();
596	if (o) {
597	f->Delete(Form("%s;1",name));
598	o->Write(name);
599	}
600
601	f->Close();
602	delete f;
603	}
604
605	//--------------------------------------------------------------------------------------------------
606
607	Double_t NBD(Int_t n, Double_t mu, Double_t k)
608	{
609	Double_t mudk = mu/k;
610	Double_t ret = TMath::Gamma(n+k) / TMath::Gamma(k) / TMath::Factorial(n) *
611	TMath::Power(mudk,n) / TMath::Power(1+mudk,n+k);
612	return ret;
613	}
614
615	//--------------------------------------------------------------------------------------------------
616
617	TH1F *NBDhist(Double_t mu, Double_t k)
618	{
619	TH1F *h = new TH1F("htemp","",100,0,100);
620	h->SetName(Form("nbd_%f_%f",mu,k));
621	h->SetDirectory(0);
622	for (Int_t i=0;i<100;++i) {
623	Double_t val = NBD(i,mu,k);
624	h->Fill(i,val);
625	}
626	return h;
627	}