[u/mrichter/AliRoot.git] / TAmpt / macros / createAmptMC.C

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <Riostream.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 <TFolder.h>
#include "AliRun.h"
#include "AliRunLoader.h"
#include "AliHeader.h"
#include "AliStack.h"
#include "AliPDG.h"
#include "AliGenAmpt.h"
#include "TAmpt.h"
#endif

class MyHeader
{
public:
  MyHeader() : fNATT(0), fEATT(0),fJATT(0),fNT(0),fNP(0),fN0(0),fN01(0),fN10(0),fN11(0),fBB(0),fRP(0) {;}
  virtual ~MyHeader() {;}
  Int_t       fNATT;
  Float_t     fEATT;
  Int_t       fJATT;
  Int_t       fNT;
  Int_t       fNP;
  Int_t       fN0;
  Int_t       fN01;
  Int_t       fN10;
  Int_t       fN11;
  Float_t     fBB;
  Float_t     fRP;
  ClassDef(MyHeader,1) // My header class
};

class MyNuc : public TObject
{
public:
  MyNuc(Int_t pid=0, Int_t st=0, Int_t type=0, Double_t x=0, Double_t y=0, Double_t z=0) :
    TObject(), fPid(pid), fSt(st), fType(type), fX(x), fY(y), fZ(z) {;}
  Double_t    X()    const { return fX; }
  Double_t    Y()    const { return fY; }
  Double_t    Z()    const { return fZ; }
  Int_t       Pid()  const { return fPid; }
  Int_t       St()   const { return fSt; }
  Int_t       Type() const { return fType; }
protected:
  Int_t       fPid;
  Int_t       fSt;
  Int_t       fType;
  Double32_t  fX;
  Double32_t  fY;
  Double32_t  fZ;
  ClassDef(MyNuc,1) // My nucleon class in cartesian coordinates
};

class MyPart : public TObject
{
public:
  MyPart(Int_t pid=0, Int_t st=0, Int_t type=0, Double_t pt=0, Double_t eta=0, Double_t phi=0) :
    TObject(), fPid(pid), fSt(st), fType(type), fPt(pt), fEta(eta), fPhi(phi) {;}
  Double_t    Px()   const { return fPt*TMath::Cos(fPhi);  }
  Double_t    Py()   const { return fPt*TMath::Sin(fPhi);  }
  Double_t    Pz()   const { return fPt*TMath::SinH(fEta); }
  Double_t    Pt()   const { return fPt;  }
  Double_t    Eta()  const { return fEta; }
  Double_t    Phi()  const { return fPhi; }
  Int_t       Pid()  const { return fPid; }
  Int_t       St()   const { return fSt; }
  Int_t       Type() const { return fType; }
protected:
  Int_t       fPid;
  Int_t       fSt;
  Int_t       fType;
  Double32_t  fPt;
  Double32_t  fEta;
  Double32_t  fPhi;
  ClassDef(MyPart,1) // My particle class in cylindrical coordinates
};

void createAmptMC(Int_t nEvents,
                  const char *outFileName = "amptout.root");

void anaAmptMC(Int_t nEvents,
               const char *inFileNames = "/eliza6/alice/loizides/ampt/run_1/ampt*.root");

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

void createAmptMC(Int_t nEvents,
                  const char *outFileName) 
{
  TClass::GetClass("MyNuc")->IgnoreTObjectStreamer();
  TClass::GetClass("MyPart")->IgnoreTObjectStreamer();

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

  AliGenAmpt *genHi = new AliGenAmpt(-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->SetSpectators(0);   // track spectators 
  genHi->SetIsoft(4);        // 4=string melting, 1=standard AMPT
  genHi->SetXmu(3.2264);     // parton xsection
  genHi->SetNtMax(150);      // time bins
  if (0) { //RHIC settings
    genHi->SetAlpha(0.47140452);
    genHi->SetStringFrag(2.2,0.5);
  }
  genHi->Init();

  TClonesArray *inucs = new TClonesArray("TParticle",10000);
  TClonesArray *parts = new TClonesArray("MyPart",10000);
  TClonesArray *nucs  = new TClonesArray("MyNuc",500);
  MyHeader *myheader = new MyHeader;

  TFile *outFile = TFile::Open(outFileName, "RECREATE");
  outFile->SetCompressionLevel(5);
  TDirectory::TContext context(outFile);

  TTree *tree = new TTree("ampt", "AMPT tree");
  tree->Branch("parts",&parts, 256*1024, 99);
  tree->Branch("nucs",&nucs, 64*1024, 99);
  tree->Branch("info",&myheader, 32*1024, 99);

  // create events and fill them
  for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {

    cout << "Event " << iEvent+1 << "/" << nEvents << endl;;
    stack->Reset();
    genHi->Generate();
    parts->Clear();
    const TClonesArray *iarr = genHi->GetParticles();
    if (iarr) {
      for(Int_t i=0;i<iarr->GetEntriesFast();++i) {
        TParticle *p = (TParticle*)iarr->At(i);
        new((*parts)[i]) MyPart(p->GetPdgCode(),
                                p->GetStatusCode(),
                                p->GetUniqueID(),
                                p->Pt(),
                                p->Eta(),
                                p->Phi());
      }
    }
    TAmpt *ampt = genHi->Ampt();
    if (ampt) {
      ampt->ImportNucleons(inucs);
      nucs->Clear();
      for(Int_t i=0;i<inucs->GetEntriesFast();++i) {
        TParticle *p = (TParticle*)inucs->At(i);
        new((*nucs)[i]) MyNuc(p->GetPdgCode(),
                              p->GetStatusCode(),
                              p->GetUniqueID(),
                              p->Vx(),
                              p->Vy(),
                              p->Vz());
      }

      myheader->fNATT=ampt->GetNATT();
      myheader->fEATT=ampt->GetEATT();
      myheader->fJATT=ampt->GetJATT();
      myheader->fNT=ampt->GetNT();
      myheader->fNP=ampt->GetNP();
      myheader->fN0=ampt->GetN0();
      myheader->fN01=ampt->GetN01();
      myheader->fN10=ampt->GetN10();
      myheader->fN11=ampt->GetN11();
      myheader->fBB=ampt->GetBB();
      myheader->fRP=ampt->GetPhi();
    }
    tree->Fill();
  } // end of event loop
  tree->Write();
  outFile->Close();
}

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

void anaAmptMC(Int_t nEvents,
               const char *inFileNames)
{

  TChain *c = new TChain("ampt");
  c->Add(inFileNames);


  TClonesArray *parts = 0;
  TClonesArray *nucs  = 0;
  MyHeader *info = 0;

  c->SetBranchAddress("info",&info);
  c->SetBranchAddress("nucs",&nucs);
  c->SetBranchAddress("parts",&parts);

  Int_t nRead = nEvents;
  if (nRead<0)
    nRead = c->GetEntries();
  else if (0 && (nRead>c->GetEntries()))
    nRead = c->GetEntries();

   for (Int_t ev=0;ev<nRead;++ev) {
      c->GetEntry(ev);

      Int_t fAN=0, fBN=0, fAStat[250], fBStat[250];
      Double_t fXA[250], fXB[250], fYA[250], fYB[250];
      for (Int_t k=0;k<nucs->GetEntries();++k) {
        MyNuc *n = (MyNuc*)nucs->At(k);
        if (n->Type()>0) {
          fAStat[fAN] =0;
          fXA[fAN] = n->X();
          fYA[fAN] = n->Y();
          ++fAN;
        } else {
          fBStat[fBN] =0;
          fXB[fBN] = n->X();
          fYB[fBN] = n->Y();
          ++fBN;
        }
      }

      // Glauber calculation
      Double_t fXSect = 65; //mbarn
      Double_t d2 = (Double_t)fXSect/(TMath::Pi()*10); // in fm^2

      // For each of the A nucleons in nucleus B
      for (Int_t i = 0; i<fBN; i++) {
        for (Int_t j = 0 ; j < fAN ;j++) {
          Double_t dx = fXB[i]-fXA[j];
          Double_t dy = fYB[i]-fYA[j];
          Double_t dij = dx*dx+dy*dy;
          if (dij < d2) {
            fBStat[i]++;
            fAStat[j]++;
          }
        }
      }
      // Calculate npart
      Int_t npart=0;
      for (Int_t i = 0; i<fAN; i++) {
        if (fAStat[i]>0) {
          npart++;
        }
      }
      for (Int_t i = 0; i<fBN; i++) {
        if (fBStat[i]>0) {
          npart++;
        }
      }
      cout << ev << " : "  << info->fBB << " np " << npart << " vs " << info->fNP+info->fNT << endl;
   }
}
Commit	Line	Data
0119ef9a	1	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	2	#include <Riostream.h>
	3	#include <TFile.h>
	4	#include <TTree.h>
	5	#include <TMath.h>
	6	#include <TParticle.h>
	7	#include <TPDGCode.h>
	8	#include <TDatabasePDG.h>
	9	#include <TRandom3.h>
	10	#include <TChain.h>
b72f0cd4	11	#include <TFolder.h>
0119ef9a	12	#include "AliRun.h"
	13	#include "AliRunLoader.h"
	14	#include "AliHeader.h"
	15	#include "AliStack.h"
	16	#include "AliPDG.h"
	17	#include "AliGenAmpt.h"
	18	#include "TAmpt.h"
	19	#endif
	20
	21	class MyHeader
	22	{
	23	public:
	24	MyHeader() : fNATT(0), fEATT(0),fJATT(0),fNT(0),fNP(0),fN0(0),fN01(0),fN10(0),fN11(0),fBB(0),fRP(0) {;}
	25	virtual ~MyHeader() {;}
	26	Int_t fNATT;
	27	Float_t fEATT;
	28	Int_t fJATT;
	29	Int_t fNT;
	30	Int_t fNP;
	31	Int_t fN0;
	32	Int_t fN01;
	33	Int_t fN10;
	34	Int_t fN11;
	35	Float_t fBB;
	36	Float_t fRP;
	37	ClassDef(MyHeader,1) // My header class
	38	};
	39
	40	class MyNuc : public TObject
	41	{
	42	public:
	43	MyNuc(Int_t pid=0, Int_t st=0, Int_t type=0, Double_t x=0, Double_t y=0, Double_t z=0) :
	44	TObject(), fPid(pid), fSt(st), fType(type), fX(x), fY(y), fZ(z) {;}
	45	Double_t X() const { return fX; }
	46	Double_t Y() const { return fY; }
	47	Double_t Z() const { return fZ; }
	48	Int_t Pid() const { return fPid; }
	49	Int_t St() const { return fSt; }
	50	Int_t Type() const { return fType; }
	51	protected:
	52	Int_t fPid;
	53	Int_t fSt;
	54	Int_t fType;
	55	Double32_t fX;
	56	Double32_t fY;
	57	Double32_t fZ;
	58	ClassDef(MyNuc,1) // My nucleon class in cartesian coordinates
	59	};
	60
	61	class MyPart : public TObject
	62	{
	63	public:
	64	MyPart(Int_t pid=0, Int_t st=0, Int_t type=0, Double_t pt=0, Double_t eta=0, Double_t phi=0) :
	65	TObject(), fPid(pid), fSt(st), fType(type), fPt(pt), fEta(eta), fPhi(phi) {;}
	66	Double_t Px() const { return fPt*TMath::Cos(fPhi); }
	67	Double_t Py() const { return fPt*TMath::Sin(fPhi); }
	68	Double_t Pz() const { return fPt*TMath::SinH(fEta); }
	69	Double_t Pt() const { return fPt; }
	70	Double_t Eta() const { return fEta; }
	71	Double_t Phi() const { return fPhi; }
	72	Int_t Pid() const { return fPid; }
	73	Int_t St() const { return fSt; }
	74	Int_t Type() const { return fType; }
	75	protected:
76	Int_t fPid;
77	Int_t fSt;
78	Int_t fType;
79	Double32_t fPt;
80	Double32_t fEta;
81	Double32_t fPhi;
82	ClassDef(MyPart,1) // My particle class in cylindrical coordinates
83	};
84
85	void createAmptMC(Int_t nEvents,
86	const char *outFileName = "amptout.root");
87
88	void anaAmptMC(Int_t nEvents,
89	const char inFileNames = "/eliza6/alice/loizides/ampt/run_1/ampt.root");
90
91	//-----------------------------------------------------------------------------------------------------
92
93	void createAmptMC(Int_t nEvents,
94	const char *outFileName)
95	{
96	TClass::GetClass("MyNuc")->IgnoreTObjectStreamer();
97	TClass::GetClass("MyPart")->IgnoreTObjectStreamer();
98
99	AliPDG::AddParticlesToPdgDataBase();
100	TDatabasePDG::Instance();
101
102	// Run loader
103	TFolder *folder = new TFolder("myfolder","myfolder");
104	AliRunLoader* rl = new AliRunLoader(folder);
105	rl->MakeHeader();
106	rl->MakeStack();
107	AliStack* stack = rl->Stack();
108	//AliHeader* rheader = rl->GetHeader();
109
110	AliGenAmpt *genHi = new AliGenAmpt(-1);
111	genHi->SetStack(stack);
112	genHi->SetEnergyCMS(2760);
113	genHi->SetReferenceFrame("CMS");
114	genHi->SetProjectile("A", 208, 82);
115	genHi->SetTarget ("A", 208, 82);
116	genHi->SetPtHardMin (2);
b72f0cd4	117	genHi->SetImpactParameterRange(0,30);
0119ef9a	118	genHi->SetJetQuenching(0); // enable jet quenching
	119	genHi->SetShadowing(1); // enable shadowing
	120	genHi->SetDecaysOff(1); // neutral pion and heavy particle decays switched off
b72f0cd4	121	genHi->SetSpectators(0); // track spectators
b72f0cd4	122	genHi->SetIsoft(4); // 4=string melting, 1=standard AMPT
0119ef9a	123	genHi->SetXmu(3.2264); // parton xsection
b72f0cd4	124	genHi->SetNtMax(150); // time bins
a004b331	125	if (0) { //RHIC settings
	126	genHi->SetAlpha(0.47140452);
	127	genHi->SetStringFrag(2.2,0.5);
	128	}
0119ef9a	129	genHi->Init();
	130
	131	TClonesArray *inucs = new TClonesArray("TParticle",10000);
	132	TClonesArray *parts = new TClonesArray("MyPart",10000);
	133	TClonesArray *nucs = new TClonesArray("MyNuc",500);
	134	MyHeader *myheader = new MyHeader;
	135
	136	TFile *outFile = TFile::Open(outFileName, "RECREATE");
	137	outFile->SetCompressionLevel(5);
	138	TDirectory::TContext context(outFile);
	139
	140	TTree *tree = new TTree("ampt", "AMPT tree");
	141	tree->Branch("parts",&parts, 256*1024, 99);
	142	tree->Branch("nucs",&nucs, 64*1024, 99);
	143	tree->Branch("info",&myheader, 32*1024, 99);
	144
	145	// create events and fill them
	146	for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
	147
	148	cout << "Event " << iEvent+1 << "/" << nEvents << endl;;
	149	stack->Reset();
	150	genHi->Generate();
	151	parts->Clear();
	152	const TClonesArray *iarr = genHi->GetParticles();
	153	if (iarr) {
	154	for(Int_t i=0;i<iarr->GetEntriesFast();++i) {
	155	TParticle p = (TParticle)iarr->At(i);
	156	new((*parts)[i]) MyPart(p->GetPdgCode(),
	157	p->GetStatusCode(),
	158	p->GetUniqueID(),
	159	p->Pt(),
	160	p->Eta(),
	161	p->Phi());
	162	}
	163	}
	164	TAmpt *ampt = genHi->Ampt();
	165	if (ampt) {
	166	ampt->ImportNucleons(inucs);
	167	nucs->Clear();
	168	for(Int_t i=0;i<inucs->GetEntriesFast();++i) {
	169	TParticle p = (TParticle)inucs->At(i);
	170	new((*nucs)[i]) MyNuc(p->GetPdgCode(),
	171	p->GetStatusCode(),
	172	p->GetUniqueID(),
	173	p->Vx(),
	174	p->Vy(),
	175	p->Vz());
	176	}
	177
	178	myheader->fNATT=ampt->GetNATT();
	179	myheader->fEATT=ampt->GetEATT();
	180	myheader->fJATT=ampt->GetJATT();
	181	myheader->fNT=ampt->GetNT();
	182	myheader->fNP=ampt->GetNP();
	183	myheader->fN0=ampt->GetN0();
	184	myheader->fN01=ampt->GetN01();
	185	myheader->fN10=ampt->GetN10();
	186	myheader->fN11=ampt->GetN11();
	187	myheader->fBB=ampt->GetBB();
	188	myheader->fRP=ampt->GetPhi();
	189	}
	190	tree->Fill();
	191	} // end of event loop
	192	tree->Write();
193	outFile->Close();
194	}
195
196	//-----------------------------------------------------------------------------------------------------
197
198	void anaAmptMC(Int_t nEvents,
199	const char *inFileNames)
200	{
201
202	TChain *c = new TChain("ampt");
203	c->Add(inFileNames);
204
205
206	TClonesArray *parts = 0;
207	TClonesArray *nucs = 0;
208	MyHeader *info = 0;
209
210	c->SetBranchAddress("info",&info);
211	c->SetBranchAddress("nucs",&nucs);
212	c->SetBranchAddress("parts",&parts);
213
214	Int_t nRead = nEvents;
215	if (nRead<0)
216	nRead = c->GetEntries();
217	else if (0 && (nRead>c->GetEntries()))
218	nRead = c->GetEntries();
219
220	for (Int_t ev=0;ev<nRead;++ev) {
221	c->GetEntry(ev);
222
223	Int_t fAN=0, fBN=0, fAStat[250], fBStat[250];
224	Double_t fXA[250], fXB[250], fYA[250], fYB[250];
225	for (Int_t k=0;k<nucs->GetEntries();++k) {
226	MyNuc n = (MyNuc)nucs->At(k);
227	if (n->Type()>0) {
228	fAStat[fAN] =0;
229	fXA[fAN] = n->X();
230	fYA[fAN] = n->Y();
231	++fAN;
232	} else {
233	fBStat[fBN] =0;
234	fXB[fBN] = n->X();
235	fYB[fBN] = n->Y();
236	++fBN;
237	}
238	}
239
240	// Glauber calculation
241	Double_t fXSect = 65; //mbarn
242	Double_t d2 = (Double_t)fXSect/(TMath::Pi()*10); // in fm^2
243
244	// For each of the A nucleons in nucleus B
245	for (Int_t i = 0; i<fBN; i++) {
246	for (Int_t j = 0 ; j < fAN ;j++) {
247	Double_t dx = fXB[i]-fXA[j];
248	Double_t dy = fYB[i]-fYA[j];
249	Double_t dij = dxdx+dydy;
250	if (dij < d2) {
251	fBStat[i]++;
252	fAStat[j]++;
253	}
254	}
255	}
256	// Calculate npart
257	Int_t npart=0;
258	for (Int_t i = 0; i<fAN; i++) {
259	if (fAStat[i]>0) {
260	npart++;
261	}
262	}
263	for (Int_t i = 0; i<fBN; i++) {
264	if (fBStat[i]>0) {
265	npart++;
266	}
267	}
268	cout << ev << " : " << info->fBB << " np " << npart << " vs " << info->fNP+info->fNT << endl;
269	}
270	}