first set of macros
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 30 Sep 2010 13:52:23 +0000 (13:52 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 30 Sep 2010 13:52:23 +0000 (13:52 +0000)
PWG2/EVCHAR/macros/cent/createHijingGlauberTestTree.C [new file with mode: 0644]
PWG2/EVCHAR/macros/cent/rootlogon.C [new file with mode: 0644]
PWG2/EVCHAR/macros/cent/testHijingGlauber.C [new file with mode: 0644]

diff --git a/PWG2/EVCHAR/macros/cent/createHijingGlauberTestTree.C b/PWG2/EVCHAR/macros/cent/createHijingGlauberTestTree.C
new file mode 100644 (file)
index 0000000..1878847
--- /dev/null
@@ -0,0 +1,206 @@
+// $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 "AliRun.h"
+#include "AliRunLoader.h"
+#include "AliHeader.h"
+#include "AliStack.h"
+#include "AliPDG.h"
+#include "THijing/AliGenHijing.h"
+#endif
+
+class MyHeader
+{
+public:
+  MyHeader() : fNATT(0), fEATT(0),fJATT(0),fNT(0),fNP(0),
+               fN00(0),fN01(0),fN10(0),fN11(0),fBB(0),fRP(0),
+               fPSn(0), fPSp(0), fTSn(0), fTSp(0) {;}
+  virtual ~MyHeader() {;}
+  Int_t       fNATT;
+  Double32_t  fEATT;
+  Int_t       fJATT;
+  Int_t       fNT;
+  Int_t       fNP;
+  Int_t       fN00;
+  Int_t       fN01;
+  Int_t       fN10;
+  Int_t       fN11;
+  Double32_t  fBB;
+  Double32_t  fRP;
+  Int_t       fPSn;
+  Int_t       fPSp;
+  Int_t       fTSn;
+  Int_t       fTSp;
+  ClassDef(MyHeader,1) // My header class
+};
+
+class MyResponse
+{
+public:
+  MyResponse() : fEtch0p(0), fEtch1p(0), fEtch2p(0), fEtch3p(0), fEtch4p(0), fEtch5p(0), fEtchrp(0), 
+                 fEtch0n(0), fEtch1n(0), fEtch2n(0), fEtch3n(0), fEtch4n(0), fEtch5n(0), fEtchrn(0),
+                 fNch0p(0), fNch1p(0), fNch2p(0), fNch3p(0), fNch4p(0), fNch5p(0), fNchrp(0),
+                 fNch0n(0), fNch1n(0), fNch2n(0), fNch3n(0), fNch4n(0), fNch5n(0), fNchrn(0) {;}
+  virtual ~MyResponse() {;}
+  Double32_t  fEtch0p;
+  Double32_t  fEtch1p;
+  Double32_t  fEtch2p;
+  Double32_t  fEtch3p;
+  Double32_t  fEtch4p;
+  Double32_t  fEtch5p;
+  Double32_t  fEtchrp;
+  Double32_t  fEtch0n;
+  Double32_t  fEtch1n;
+  Double32_t  fEtch2n;
+  Double32_t  fEtch3n;
+  Double32_t  fEtch4n;
+  Double32_t  fEtch5n;
+  Double32_t  fEtchrn;
+  Int_t       fNch0p;
+  Int_t       fNch1p;
+  Int_t       fNch2p;
+  Int_t       fNch3p;
+  Int_t       fNch4p;
+  Int_t       fNch5p;
+  Int_t       fNchrp;
+  Int_t       fNch0n;
+  Int_t       fNch1n;
+  Int_t       fNch2n;
+  Int_t       fNch3n;
+  Int_t       fNch4n;
+  Int_t       fNch5n;
+  Int_t       fNchrn;
+  ClassDef(MyResponse,1) // My reponse class
+};
+
+void createGlauberTree(Int_t nEvents,
+                       const char *outFileName = "treeout.root");
+
+//-----------------------------------------------------------------------------------------------------
+
+void createGlauberTree(Int_t nEvents,
+                       const char *outFileName) 
+{
+  AliPDG::AddParticlesToPdgDataBase();
+  TDatabasePDG::Instance();
+
+  // Run loader
+  TFolder *folder = new TFolder("myfolder","myfolder");
+  AliRunLoader* rl = new AliRunLoader(folder);
+  rl->MakeHeader();
+  rl->MakeStack();
+  AliStack* stack = rl->Stack();
+  //AliHeader* rheader = rl->GetHeader();
+
+  AliGenHijing *genHi = new AliGenHijing(-1);
+  genHi->SetStack(stack);
+  genHi->SetEnergyCMS(2760);
+  genHi->SetReferenceFrame("CMS");
+  genHi->SetProjectile("A", 208, 82);
+  genHi->SetTarget    ("A", 208, 82);
+  genHi->SetPtHardMin (2);
+  genHi->SetImpactParameterRange(0.,30);
+  genHi->SetJetQuenching(0); // enable jet quenching
+  genHi->SetShadowing(1);    // enable shadowing
+  genHi->SetDecaysOff(1);    // neutral pion and heavy particle decays switched off
+  genHi->Init();
+
+  MyHeader  *myheader = new MyHeader;
+  MyResponse *myresp  = new MyResponse;
+
+  TFile *outFile = TFile::Open(outFileName, "RECREATE");
+  outFile->SetCompressionLevel(5);
+  TDirectory::TContext context(outFile);
+
+  TTree *tree = new TTree("glaubertree", "Glauber tree");
+  tree->Branch("header",&myheader, 32*1024, 99);
+  tree->Branch("response",&myresp, 32*1024, 99);
+
+  Double_t etas[] = {-10,-5,-4,-3,-2,-1,0,1,2,3,4,5,10};
+  TH1D *hNEta = new TH1D("hNeta","",12,etas);
+  TH1D *hEtEta = new TH1D("hEteta","",12,etas);
+
+  // create events and fill them
+  for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
+
+    cout << "Event " << iEvent+1 << "/" << nEvents << endl;;
+    stack->Reset();
+    hNEta->Reset();
+    hEtEta->Reset();
+    genHi->Generate();
+  
+    AliStack *s = genHi->GetStack();
+    const TObjArray *parts = s->Particles();
+    Int_t nents = parts->GetEntries();
+    for (Int_t i = 0; i<nents; ++i) {
+      TParticle *p = (TParticle*)parts->At(i);
+      //p->Print();
+      TParticlePDG *pdg = p->GetPDG(1);
+      Int_t c = (Int_t)(TMath::Abs(pdg->Charge()));
+      if (c!=0) {
+        hNEta->Fill(p->Eta());
+        hEtEta->Fill(p->Eta(),p->Pt());
+      }
+    }
+
+    AliGenHijingEventHeader *h = (AliGenHijingEventHeader*)genHi->CollisionGeometry();
+    myheader->fNATT = nents;
+    myheader->fEATT = h->TotalEnergy();
+    myheader->fJATT = h->HardScatters();
+    myheader->fNT   = h->TargetParticipants();
+    myheader->fNP   = h->ProjectileParticipants();
+    myheader->fN00  = h->NwNw();
+    myheader->fN01  = h->NwN();
+    myheader->fN10  = h->NNw();
+    myheader->fN11  = h->NN();
+    myheader->fBB   = h->ImpactParameter();
+    myheader->fRP   = h->ReactionPlaneAngle();
+    myheader->fPSn  = h->ProjSpectatorsn();
+    myheader->fPSp  = h->ProjSpectatorsp();
+    myheader->fTSn  = h->TargSpectatorsn();
+    myheader->fTSp  = h->TargSpectatorsn();
+
+    myresp->fEtch0p = hEtEta->GetBinContent(hEtEta->FindBin(0.5));
+    myresp->fEtch1p = hEtEta->GetBinContent(hEtEta->FindBin(1.5));
+    myresp->fEtch2p = hEtEta->GetBinContent(hEtEta->FindBin(2.5));
+    myresp->fEtch3p = hEtEta->GetBinContent(hEtEta->FindBin(3.5));
+    myresp->fEtch4p = hEtEta->GetBinContent(hEtEta->FindBin(4.5));
+    myresp->fEtch5p = hEtEta->GetBinContent(hEtEta->FindBin(5.5));
+    myresp->fEtchrp = hEtEta->GetBinContent(hEtEta->FindBin(10.5));
+    myresp->fEtch0n = hEtEta->GetBinContent(hEtEta->FindBin(-0.5));
+    myresp->fEtch1n = hEtEta->GetBinContent(hEtEta->FindBin(-1.5));
+    myresp->fEtch2n = hEtEta->GetBinContent(hEtEta->FindBin(-2.5));
+    myresp->fEtch3n = hEtEta->GetBinContent(hEtEta->FindBin(-3.5));
+    myresp->fEtch4n = hEtEta->GetBinContent(hEtEta->FindBin(-4.5));
+    myresp->fEtch5n = hEtEta->GetBinContent(hEtEta->FindBin(-5.5));
+    myresp->fEtchrn = hEtEta->GetBinContent(hEtEta->FindBin(-10.5));
+    myresp->fNch0p  = hNEta->GetBinContent(hNEta->FindBin(0.5));
+    myresp->fNch1p  = hNEta->GetBinContent(hNEta->FindBin(1.5));
+    myresp->fNch2p  = hNEta->GetBinContent(hNEta->FindBin(2.5));
+    myresp->fNch3p  = hNEta->GetBinContent(hNEta->FindBin(3.5));
+    myresp->fNch4p  = hNEta->GetBinContent(hNEta->FindBin(4.5));
+    myresp->fNch5p  = hNEta->GetBinContent(hNEta->FindBin(5.5));
+    myresp->fNchrp  = hNEta->GetBinContent(hNEta->FindBin(10.5));
+    myresp->fNch0n  = hNEta->GetBinContent(hNEta->FindBin(-0.5));
+    myresp->fNch1n  = hNEta->GetBinContent(hNEta->FindBin(-1.5));
+    myresp->fNch2n  = hNEta->GetBinContent(hNEta->FindBin(-2.5));
+    myresp->fNch3n  = hNEta->GetBinContent(hNEta->FindBin(-3.5));
+    myresp->fNch4n  = hNEta->GetBinContent(hNEta->FindBin(-4.5));
+    myresp->fNch5n  = hNEta->GetBinContent(hNEta->FindBin(-5.5));
+    myresp->fNchrn  = hNEta->GetBinContent(hNEta->FindBin(-10.5));
+
+    tree->Fill();
+  } // end of event loop
+  tree->Write();
+  outFile->Close();
+}
diff --git a/PWG2/EVCHAR/macros/cent/rootlogon.C b/PWG2/EVCHAR/macros/cent/rootlogon.C
new file mode 100644 (file)
index 0000000..bcbd50e
--- /dev/null
@@ -0,0 +1,39 @@
+// $Id$
+{
+  gSystem->Load("libTree.so");
+  gSystem->Load("libGui.so");
+  gSystem->Load("libVMC.so");
+  gSystem->Load("libSTEERBase.so");
+  gSystem->Load("libESD.so");
+  gSystem->Load("libAOD.so"); 
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libANALYSISalice");
+  gSystem->Load("libCDB.so");
+  gSystem->Load("libMinuit.so");
+  gSystem->Load("libMinuit2.so");
+  gSystem->Load("libProof.so");
+  gSystem->Load("libRAWDatabase.so");
+  gSystem->Load("libSTEER.so");
+  gSystem->Load("libPhysics.so");
+  gSystem->Load("libEVGEN.so");
+  gSystem->Load("libFASTSIM.so");
+  if (1) {
+    gSystem->Load("libhijing.so");
+    gSystem->Load("libTHijing.so");
+  } else {
+    gSystem->Load("libampt.so");
+    gSystem->Load("libTAmpt.so");
+  }
+  gSystem->Load("libEGPythia6.so");
+  gSystem->Load("libpythia6.so");
+  gSystem->Load("libAliPythia6.so");
+  gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include -I$ALICE_ROOT");
+  if (gSystem->Getenv("TMPDIR")) 
+    gSystem->SetBuildDir(gSystem->Getenv("TMPDIR"));
+  if (!gAlice)
+    new AliRun("gAlice","The ALICE Off-line Simulation Framework");
+  if (1) {
+    delete gRandom;
+    gRandom = new TRandom3(0);
+  }
+}
diff --git a/PWG2/EVCHAR/macros/cent/testHijingGlauber.C b/PWG2/EVCHAR/macros/cent/testHijingGlauber.C
new file mode 100644 (file)
index 0000000..7e2151a
--- /dev/null
@@ -0,0 +1,578 @@
+// $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>
+#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 testHijingGlauber() 
+{
+  Bool_t pImDist               = 0;
+  Bool_t pNpDist               = 0;
+  Bool_t pNpDistSelWithImp     = 0;
+  Bool_t pNpDistSelWithImpFits = 0;
+  Bool_t pMidRecResStudy       = 0;
+
+  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 (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;
+}