1 #if !defined(__CINT__) || defined(__MAKECINT__)
8 #include <TDatabasePDG.h>
13 #include "AliRunLoader.h"
14 #include "AliHeader.h"
17 #include "AliGenAmpt.h"
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() {;}
37 ClassDef(MyHeader,1) // My header class
40 class MyNuc : public TObject
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; }
58 ClassDef(MyNuc,1) // My nucleon class in cartesian coordinates
61 class MyPart : public TObject
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; }
82 ClassDef(MyPart,1) // My particle class in cylindrical coordinates
85 void createAmptMC(Int_t nEvents,
86 const char *outFileName = "amptout.root");
88 void anaAmptMC(Int_t nEvents,
89 const char *inFileNames = "/eliza6/alice/loizides/ampt/run_1/ampt*.root");
91 //-----------------------------------------------------------------------------------------------------
93 void createAmptMC(Int_t nEvents,
94 const char *outFileName)
96 TClass::GetClass("MyNuc")->IgnoreTObjectStreamer();
97 TClass::GetClass("MyPart")->IgnoreTObjectStreamer();
99 AliPDG::AddParticlesToPdgDataBase();
100 TDatabasePDG::Instance();
103 TFolder *folder = new TFolder("myfolder","myfolder");
104 AliRunLoader* rl = new AliRunLoader(folder);
107 AliStack* stack = rl->Stack();
108 //AliHeader* rheader = rl->GetHeader();
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);
117 genHi->SetImpactParameterRange(0,30);
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
121 genHi->SetSpectators(0); // track spectators
122 genHi->SetIsoft(4); // 4=string melting, 1=standard AMPT
123 genHi->SetXmu(3.2264); // parton xsection
124 genHi->SetNtMax(150); // time bins
125 if (0) { //RHIC settings
126 genHi->SetAlpha(0.47140452);
127 genHi->SetStringFrag(2.2,0.5);
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;
136 TFile *outFile = TFile::Open(outFileName, "RECREATE");
137 outFile->SetCompressionLevel(5);
138 TDirectory::TContext context(outFile);
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);
145 // create events and fill them
146 for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
148 cout << "Event " << iEvent+1 << "/" << nEvents << endl;;
152 const TClonesArray *iarr = genHi->GetParticles();
154 for(Int_t i=0;i<iarr->GetEntriesFast();++i) {
155 TParticle *p = (TParticle*)iarr->At(i);
156 new((*parts)[i]) MyPart(p->GetPdgCode(),
164 TAmpt *ampt = genHi->Ampt();
166 ampt->ImportNucleons(inucs);
168 for(Int_t i=0;i<inucs->GetEntriesFast();++i) {
169 TParticle *p = (TParticle*)inucs->At(i);
170 new((*nucs)[i]) MyNuc(p->GetPdgCode(),
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();
191 } // end of event loop
196 //-----------------------------------------------------------------------------------------------------
198 void anaAmptMC(Int_t nEvents,
199 const char *inFileNames)
202 TChain *c = new TChain("ampt");
206 TClonesArray *parts = 0;
207 TClonesArray *nucs = 0;
210 c->SetBranchAddress("info",&info);
211 c->SetBranchAddress("nucs",&nucs);
212 c->SetBranchAddress("parts",&parts);
214 Int_t nRead = nEvents;
216 nRead = c->GetEntries();
217 else if (0 && (nRead>c->GetEntries()))
218 nRead = c->GetEntries();
220 for (Int_t ev=0;ev<nRead;++ev) {
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);
240 // Glauber calculation
241 Double_t fXSect = 65; //mbarn
242 Double_t d2 = (Double_t)fXSect/(TMath::Pi()*10); // in fm^2
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 = dx*dx+dy*dy;
258 for (Int_t i = 0; i<fAN; i++) {
263 for (Int_t i = 0; i<fBN; i++) {
268 cout << ev << " : " << info->fBB << " np " << npart << " vs " << info->fNP+info->fNT << endl;