[u/mrichter/AliRoot.git] / PWG2 / EVCHAR / macros / cent / createHijingGlauberTestTree.C

// $Id$

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <Riostream.h>
#include <TH1D.h>
#include <TFile.h>
#include <TTree.h>
#include <TNtuple.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.3);
  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);

  TNtuple *ntuple = new TNtuple("gnt", "Glauber ntuple", "npart:ncoll:b");

  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();

    if (ntuple) {
      Int_t np = h->TargetParticipants() + h->ProjectileParticipants();
      Int_t nc = h->NwNw() + h->NwN() + h->NNw() + h->NN();
      Double_t b = h->ImpactParameter();
      ntuple->Fill(np,nc,b);
    }

  } // end of event loop

  tree->Write();
  ntuple->Write();
  outFile->Close();
}
Commit	Line	Data
5957f262	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>
b24d57aa	8	#include <TNtuple.h>
5957f262	9	#include <TMath.h>
	10	#include <TParticle.h>
	11	#include <TPDGCode.h>
	12	#include <TDatabasePDG.h>
	13	#include <TRandom3.h>
	14	#include <TChain.h>
	15	#include "AliRun.h"
	16	#include "AliRunLoader.h"
	17	#include "AliHeader.h"
	18	#include "AliStack.h"
	19	#include "AliPDG.h"
	20	#include "THijing/AliGenHijing.h"
	21	#endif
	22
	23	class MyHeader
	24	{
	25	public:
	26	MyHeader() : fNATT(0), fEATT(0),fJATT(0),fNT(0),fNP(0),
	27	fN00(0),fN01(0),fN10(0),fN11(0),fBB(0),fRP(0),
	28	fPSn(0), fPSp(0), fTSn(0), fTSp(0) {;}
	29	virtual ~MyHeader() {;}
	30	Int_t fNATT;
	31	Double32_t fEATT;
	32	Int_t fJATT;
	33	Int_t fNT;
	34	Int_t fNP;
	35	Int_t fN00;
	36	Int_t fN01;
	37	Int_t fN10;
	38	Int_t fN11;
	39	Double32_t fBB;
	40	Double32_t fRP;
	41	Int_t fPSn;
	42	Int_t fPSp;
	43	Int_t fTSn;
	44	Int_t fTSp;
	45	ClassDef(MyHeader,1) // My header class
	46	};
	47
	48	class MyResponse
	49	{
	50	public:
	51	MyResponse() : fEtch0p(0), fEtch1p(0), fEtch2p(0), fEtch3p(0), fEtch4p(0), fEtch5p(0), fEtchrp(0),
	52	fEtch0n(0), fEtch1n(0), fEtch2n(0), fEtch3n(0), fEtch4n(0), fEtch5n(0), fEtchrn(0),
	53	fNch0p(0), fNch1p(0), fNch2p(0), fNch3p(0), fNch4p(0), fNch5p(0), fNchrp(0),
	54	fNch0n(0), fNch1n(0), fNch2n(0), fNch3n(0), fNch4n(0), fNch5n(0), fNchrn(0) {;}
	55	virtual ~MyResponse() {;}
	56	Double32_t fEtch0p;
	57	Double32_t fEtch1p;
	58	Double32_t fEtch2p;
	59	Double32_t fEtch3p;
	60	Double32_t fEtch4p;
	61	Double32_t fEtch5p;
	62	Double32_t fEtchrp;
	63	Double32_t fEtch0n;
	64	Double32_t fEtch1n;
	65	Double32_t fEtch2n;
	66	Double32_t fEtch3n;
	67	Double32_t fEtch4n;
	68	Double32_t fEtch5n;
	69	Double32_t fEtchrn;
	70	Int_t fNch0p;
	71	Int_t fNch1p;
	72	Int_t fNch2p;
73	Int_t fNch3p;
74	Int_t fNch4p;
75	Int_t fNch5p;
76	Int_t fNchrp;
77	Int_t fNch0n;
78	Int_t fNch1n;
79	Int_t fNch2n;
80	Int_t fNch3n;
81	Int_t fNch4n;
82	Int_t fNch5n;
83	Int_t fNchrn;
84	ClassDef(MyResponse,1) // My reponse class
85	};
86
87	void createGlauberTree(Int_t nEvents,
88	const char *outFileName = "treeout.root");
89
90	//-----------------------------------------------------------------------------------------------------
91
92	void createGlauberTree(Int_t nEvents,
93	const char *outFileName)
94	{
95	AliPDG::AddParticlesToPdgDataBase();
96	TDatabasePDG::Instance();
97
98	// Run loader
99	TFolder *folder = new TFolder("myfolder","myfolder");
100	AliRunLoader* rl = new AliRunLoader(folder);
101	rl->MakeHeader();
102	rl->MakeStack();
103	AliStack* stack = rl->Stack();
104	//AliHeader* rheader = rl->GetHeader();
105
106	AliGenHijing *genHi = new AliGenHijing(-1);
107	genHi->SetStack(stack);
108	genHi->SetEnergyCMS(2760);
109	genHi->SetReferenceFrame("CMS");
110	genHi->SetProjectile("A", 208, 82);
111	genHi->SetTarget ("A", 208, 82);
4508fe8b	112	genHi->SetPtHardMin (2.3);
5957f262	113	genHi->SetImpactParameterRange(0.,30);
	114	genHi->SetJetQuenching(0); // enable jet quenching
	115	genHi->SetShadowing(1); // enable shadowing
	116	genHi->SetDecaysOff(1); // neutral pion and heavy particle decays switched off
	117	genHi->Init();
	118
	119	MyHeader *myheader = new MyHeader;
	120	MyResponse *myresp = new MyResponse;
	121
	122	TFile *outFile = TFile::Open(outFileName, "RECREATE");
	123	outFile->SetCompressionLevel(5);
	124	TDirectory::TContext context(outFile);
	125
	126	TTree *tree = new TTree("glaubertree", "Glauber tree");
	127	tree->Branch("header",&myheader, 32*1024, 99);
	128	tree->Branch("response",&myresp, 32*1024, 99);
	129
b24d57aa	130	TNtuple *ntuple = new TNtuple("gnt", "Glauber ntuple", "npart:ncoll:b");
b24d57aa	131
5957f262	132	Double_t etas[] = {-10,-5,-4,-3,-2,-1,0,1,2,3,4,5,10};
	133	TH1D *hNEta = new TH1D("hNeta","",12,etas);
	134	TH1D *hEtEta = new TH1D("hEteta","",12,etas);
	135
	136	// create events and fill them
	137	for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
	138
	139	cout << "Event " << iEvent+1 << "/" << nEvents << endl;;
	140	stack->Reset();
	141	hNEta->Reset();
	142	hEtEta->Reset();
	143	genHi->Generate();
	144
	145	AliStack *s = genHi->GetStack();
	146	const TObjArray *parts = s->Particles();
	147	Int_t nents = parts->GetEntries();
	148	for (Int_t i = 0; i<nents; ++i) {
	149	TParticle p = (TParticle)parts->At(i);
	150	//p->Print();
	151	TParticlePDG *pdg = p->GetPDG(1);
	152	Int_t c = (Int_t)(TMath::Abs(pdg->Charge()));
	153	if (c!=0) {
	154	hNEta->Fill(p->Eta());
	155	hEtEta->Fill(p->Eta(),p->Pt());
	156	}
	157	}
	158
	159	AliGenHijingEventHeader h = (AliGenHijingEventHeader)genHi->CollisionGeometry();
	160	myheader->fNATT = nents;
	161	myheader->fEATT = h->TotalEnergy();
	162	myheader->fJATT = h->HardScatters();
	163	myheader->fNT = h->TargetParticipants();
	164	myheader->fNP = h->ProjectileParticipants();
	165	myheader->fN00 = h->NwNw();
	166	myheader->fN01 = h->NwN();
	167	myheader->fN10 = h->NNw();
	168	myheader->fN11 = h->NN();
	169	myheader->fBB = h->ImpactParameter();
	170	myheader->fRP = h->ReactionPlaneAngle();
	171	myheader->fPSn = h->ProjSpectatorsn();
	172	myheader->fPSp = h->ProjSpectatorsp();
	173	myheader->fTSn = h->TargSpectatorsn();
	174	myheader->fTSp = h->TargSpectatorsn();
	175
	176	myresp->fEtch0p = hEtEta->GetBinContent(hEtEta->FindBin(0.5));
	177	myresp->fEtch1p = hEtEta->GetBinContent(hEtEta->FindBin(1.5));
	178	myresp->fEtch2p = hEtEta->GetBinContent(hEtEta->FindBin(2.5));
	179	myresp->fEtch3p = hEtEta->GetBinContent(hEtEta->FindBin(3.5));
	180	myresp->fEtch4p = hEtEta->GetBinContent(hEtEta->FindBin(4.5));
	181	myresp->fEtch5p = hEtEta->GetBinContent(hEtEta->FindBin(5.5));
	182	myresp->fEtchrp = hEtEta->GetBinContent(hEtEta->FindBin(10.5));
	183	myresp->fEtch0n = hEtEta->GetBinContent(hEtEta->FindBin(-0.5));
	184	myresp->fEtch1n = hEtEta->GetBinContent(hEtEta->FindBin(-1.5));
	185	myresp->fEtch2n = hEtEta->GetBinContent(hEtEta->FindBin(-2.5));
	186	myresp->fEtch3n = hEtEta->GetBinContent(hEtEta->FindBin(-3.5));
	187	myresp->fEtch4n = hEtEta->GetBinContent(hEtEta->FindBin(-4.5));
	188	myresp->fEtch5n = hEtEta->GetBinContent(hEtEta->FindBin(-5.5));
	189	myresp->fEtchrn = hEtEta->GetBinContent(hEtEta->FindBin(-10.5));
	190	myresp->fNch0p = hNEta->GetBinContent(hNEta->FindBin(0.5));
	191	myresp->fNch1p = hNEta->GetBinContent(hNEta->FindBin(1.5));
	192	myresp->fNch2p = hNEta->GetBinContent(hNEta->FindBin(2.5));
	193	myresp->fNch3p = hNEta->GetBinContent(hNEta->FindBin(3.5));
	194	myresp->fNch4p = hNEta->GetBinContent(hNEta->FindBin(4.5));
	195	myresp->fNch5p = hNEta->GetBinContent(hNEta->FindBin(5.5));
196	myresp->fNchrp = hNEta->GetBinContent(hNEta->FindBin(10.5));
197	myresp->fNch0n = hNEta->GetBinContent(hNEta->FindBin(-0.5));
198	myresp->fNch1n = hNEta->GetBinContent(hNEta->FindBin(-1.5));
199	myresp->fNch2n = hNEta->GetBinContent(hNEta->FindBin(-2.5));
200	myresp->fNch3n = hNEta->GetBinContent(hNEta->FindBin(-3.5));
201	myresp->fNch4n = hNEta->GetBinContent(hNEta->FindBin(-4.5));
202	myresp->fNch5n = hNEta->GetBinContent(hNEta->FindBin(-5.5));
203	myresp->fNchrn = hNEta->GetBinContent(hNEta->FindBin(-10.5));
204
205	tree->Fill();
b24d57aa	206
	207	if (ntuple) {
	208	Int_t np = h->TargetParticipants() + h->ProjectileParticipants();
	209	Int_t nc = h->NwNw() + h->NwN() + h->NNw() + h->NN();
	210	Double_t b = h->ImpactParameter();
	211	ntuple->Fill(np,nc,b);
	212	}
	213
5957f262	214	} // end of event loop
b24d57aa	215
5957f262	216	tree->Write();
b24d57aa	217	ntuple->Write();
5957f262	218	outFile->Close();
5957f262	219	}