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