PWG2/EVCHAR/macros/cent/testHijingGlauber.C

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