--- /dev/null
+/*******************************************************************************
+ * *
+ * HYDJET++ , event generator under the ROOT FRAMEWORK for simulation of *
+ * relativistic heavy ion AA collisions as the superposition of soft, *
+ * hydro-type state and hard, multi-parton state. *
+ * *
+ * The main routine is written in the object-oriented C++ language *
+ * under the ROOT environment. The hard, multi-partonic part of *
+ * HYDJET++ event is identical to the hard part of Fortran-written *
+ * HYDJET (PYTHIA6.4xx + PYQUEN1.5) and is included in the generator *
+ * structure as the separate directory. The soft part of HYDJET++ *
+ * event represents the "thermal" hadronic state obtained with the *
+ * parameterization Bjorken-like of freeze-out hypersurface and *
+ * includes longitudinal, radial and elliptic flow effects and *
+ * decays of hadronic resonances. The corresponding fast *
+ * Monte-Carlo simulation procedure (C++ code) FAST MC is adapted. *
+ * -------------------------------------------------------------- *
+ * Web-page: *
+ * http://cern.ch/lokhtin/hydjet++ *
+ * -------------------------------------------------------------- *
+ * * *
+ * *
+ * This program is a free software; you can use and redistribute it freely. *
+ * Any publication of results obtained using this code must reference *
+ * *
+ * *
+ * *
+ * Main reference for HYDJET++: *
+ * I.P. Lokhtin, L.V. Malinina, S.V. Petrushanko, A.M. Snigirev, *
+ * I. Arsene, K. Tywoniuk, submitted to Comp. Phys. Comm. *
+ * *
+ * Reference for HYDJET and PYQUEN: *
+ * I.P. Lokhtin, A.M. Snigirev, Eur. Phys. J. C 46 (2006) 211; *
+ * http://cern.ch/lokhtin/hydro/hydjet.html *
+ * http://cern.ch/lokhtin/pyquen. *
+ * *
+ * Reference for PYTHIA6.4: *
+ * T.Sjostrand, S. Mrenna and P. Skands, JHEP05 (2006) 026; *
+ * http://home.thep.lu.se/~torbjorn/Pythia.html. *
+ * *
+ * References for FAST MC: *
+ * N.S. Amelin, R. Lednicky, T.A. Pocheptsov, I.P. Lokhtin, *
+ * L.V. Malinina, A.M. Snigirev, Iu.A. Karpenko and Yu.M. Sinyukov, *
+ * Phys. Rev. C 74 (2006) 064901; *
+ * N.S. Amelin, I. Arsene, L. Bravina, Iu.A. Karpenko, R. Lednicky, *
+ * I.P. Lokhtin, L.V. Malinina, A.M. Snigirev and Yu.M. Sinyukov, *
+ * Phys. Rev. C 77 (2008) 014903; *
+ * http://uhkm.jinr.ru. *
+ * *
+ * Reference for nuclear shadowing model: *
+ * K. Tywoniuk, I.C. Arsene, L. Bravina, A. Kaidalov and *
+ * E. Zabrodin, Phys. Lett. B 657 (2007) 170. *
+ * *
+ * version 2.0: *
+ * *
+ * Igor Lokhtin, SINP MSU, Moscow, RU *
+ * e-mail: Igor.Lokhtin@cern.ch *
+ * *
+ * Ludmila Malinina, SINP MSU, Moscow, RU *
+ * e-mail: malinina@lav01.sinp.msu.ru *
+ * *
+ *******************************************************************************/
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <time.h>
+
+#include <TNtuple.h>
+#include <TError.h>
+#include <TTree.h>
+#include <TH1D.h>
+#include <TFile.h>
+
+#include "InitialState.h"
+#include "InitialStateHydjet.h"
+
+
+#include <TRandom.h>
+
+#include "Particle.h"
+//#include "HYJET_COMMONS.h"
+//extern SERVICECommon SERVICE;
+
+
+//Main program:
+//reads input parameters from file "RunInputBjorken" or "RunInputHubble";
+//calculates particle densities and average initial multiplicities and writes them
+//in output file "multiplicities.txt";
+//creates trees (tree with direct hadrons and hadrons after resonance decays)
+//with space-time and momentum-energy information of produced hadrons;
+//writes trees in file "RunOutput.root".
+
+Int_t main() {
+
+ clock_t start;
+ start = clock();
+
+
+//new
+ time_t now;
+ struct tm *ts;
+ char buf[80];
+
+ // Get the current time
+ time(&now);
+
+ // Format and print the time, "ddd yyyy-mm-dd hh:mm:ss zzz"
+ ts = localtime(&now);
+ strftime(buf, sizeof(buf), "%a %Y-%m-%d %H:%M:%S %Z", ts);
+ printf("%s\n", buf);
+
+
+
+
+ TFile *outputFile=new TFile("RunOutput.root", "RECREATE");
+
+ //SET MAXIMAl VALUE OF PARTICLE MULTIPLICITY!!!
+ const Int_t kMax = 500000;
+ //define hadron number
+ Int_t ntot;
+ //define event number
+ Int_t nev;
+ //define hadron characteristic vectors
+ std::vector<Int_t> pdg(kMax); //pdg encodings
+ std::vector<Int_t> Mpdg(kMax);//pdg encodings for mother hadrons
+ std::vector<Int_t> type(kMax);//type: 0-from hydro or decay, 1111 from jets
+ std::vector<Float_t> Px(kMax);//x-hadron momentum component,[GeV/c]
+ std::vector<Float_t> Py(kMax);//y-hadron momentum component,[GeV/c]
+ std::vector<Float_t> Pz(kMax);//z-hadron momentum component,[GeV/c]
+ std::vector<Float_t> E(kMax); //hadron total energy,[GeV]
+ std::vector<Float_t> X(kMax);//x-hadron coordinate component,[fm]
+ std::vector<Float_t> Y(kMax);//y-hadron coordinate component,[fm]
+ std::vector<Float_t> Z(kMax);//z-hadron coordinate component,[fm]
+ std::vector<Float_t> T(kMax);//hadron time,[fm/c]
+
+ TH1D *hpt1 = new TH1D("hpt1", "hpt1", 100, 0., 20.);
+ TH1D *hpt1j = new TH1D("hpt1j", "hpt1j", 100, 0., 20.);
+ TH1D *hpt1h = new TH1D("hpt1h", "hpt1h", 100, 0., 20.);
+
+ TH1D *hv2 = new TH1D("hv2", "hv2", 100, 0.0, 10.);
+ TH1D *hv0 = new TH1D("hv0", "hv0", 100, 0.0, 10.);
+
+ TH1D *hy = new TH1D("hy", "hy", 51, -5.1, 5.1);
+ TH1D *hyjets = new TH1D("hyjets", "hyjets", 51, -5.1, 5.1);
+ TH1D *hyhydro = new TH1D("hyhydro", "hyhydro", 51, -5.1, 5.1);
+
+
+ double pdg1, Mpdg1, Px1, Py1, E1, Z1, Pz1, pt, phi, v2, eta;
+ int type1;
+
+ InitialState *FASTMC;
+
+ FASTMC = new InitialStateHydjet();
+
+ if(!FASTMC->ReadParams()) {
+ Error("RunHadronSource::main", "No initial model parameters found!!\n");
+ return 0;
+ }
+
+
+ if(!FASTMC->MultIni()) {
+ Error("RunHadronSource::main", "Initial multiplicities are zero!!\n");
+ return 0;
+ }
+
+ ParticleAllocator allocator;
+ List_t source;
+ List_t secondaries;
+ std::cout << "Generating " << FASTMC->GetNev() << " events" << std::endl;
+ std::cout << "Starting the event loop" << std::endl;
+
+
+ // Loop over events
+ for(Int_t ev = 0; ev < FASTMC->GetNev(); ++ev) {
+ nev = ev;
+ // Initialize the source
+ FASTMC->Initialize(source, allocator);
+ if(source.empty()) {
+ Error("RunHadronSource::main", "Source is not initialized!!");
+ //return 0;
+ continue;
+ }
+
+ // Run the decays //fDecay
+ if(FASTMC->GetTime() >= 0.)
+ FASTMC->Evolve(source, secondaries, allocator, FASTMC->GetWeakDecayLimit());
+
+ std::cout << "event #" << ev << "\r" << std::flush;
+// npart = 0;
+ LPIT_t it;
+ LPIT_t e;
+
+ // Fill the decayed tree
+// npart = 0;
+
+ for(it = secondaries.begin(), e = secondaries.end(); it != e; ++it) {
+ TVector3 pos(it->Pos().Vect());
+ TVector3 mom(it->Mom().Vect());
+ Float_t m1 = it->TableMass();
+ pdg1 = it->Encoding();
+ Mpdg1 = it->GetLastMotherPdg();
+ Px1 = mom[0];
+ Py1 = mom[1];
+ Pz1 = mom[2];
+ E1 = TMath::Sqrt(mom.Mag2() + m1*m1);
+ type1 = it->GetType();
+ if(pdg1==211 && abs(0.5*log((E1+Pz1)/(E1-Pz1)))<1.) {
+ hpt1->Fill(sqrt(Px1*Px1+Py1*Py1),1./sqrt(Px1*Px1+Py1*Py1));
+ }
+
+ if(pdg1==211 && abs(0.5*log((E1+Pz1)/(E1-Pz1)))<1. && type1==0) hpt1h->Fill(sqrt(Px1*Px1+Py1*Py1),1./sqrt(Px1*Px1+Py1*Py1));
+ if(pdg1==211 && abs(0.5*log((E1+Pz1)/(E1-Pz1)))<1. && type1==1)hpt1j->Fill(sqrt(Px1*Px1+Py1*Py1),1./sqrt(Px1*Px1+Py1*Py1));
+
+ if(((abs(pdg1)==211)||(abs(pdg1)==321)||(abs(pdg1)==2212))
+ && (abs(0.5*log((E1+Pz1)/(E1-Pz1)))<1.0)){
+ pt = TMath::Sqrt(Px1*Px1+Py1*Py1);
+ phi = TMath::ATan2(Py1,Px1);
+ v2 = TMath::Cos(2*phi);
+ hv2->Fill(pt,v2);
+ hv0->Fill(pt,1.);
+ }
+
+ if((abs(pdg1)==211)||(abs(pdg1)==321)||(abs(pdg1)==2212)){
+ eta=0.5*TMath::Log((sqrt(Px1*Px1+Py1*Py1+Pz1*Pz1)+Pz1)/(sqrt(Px1*Px1+Py1*Py1+Pz1*Pz1)-Pz1));
+ if(type1==1)hyjets->Fill(eta);
+ if(type1==0)hyhydro->Fill(eta);
+ hy->Fill(eta);
+ }
+
+ // npar++;
+ // if(npart > kMax)
+ // Error("in main:", "npart is too large %d", npart);
+
+
+ }
+
+ allocator.FreeList(source);
+ allocator.FreeList(secondaries);
+
+
+ }
+
+ hpt1->Write();
+ hpt1h->Write();
+ hpt1j->Write();
+ hv2->Write();
+ hv0->Write();
+ hyhydro->Write();
+ hyjets->Write();
+ hy->Write();
+
+ clock_t stop;
+ stop = clock();
+ std::cout << "*********************************************" << std::endl;
+ std::cout << "Execution time: " << (stop - start)/CLOCKS_PER_SEC << " seconds" << std::endl;
+ std::cout << "*********************************************" << std::endl;
+
+
+//new
+ time_t now1;
+ struct tm *ts1;
+ char buf1[80];
+
+ // Get the current time
+ time(&now1);
+
+ // Format and print the time, "ddd yyyy-mm-dd hh:mm:ss zzz"
+ ts1 = localtime(&now1);
+ strftime(buf1, sizeof(buf1), "%a %Y-%m-%d %H:%M:%S %Z", ts1);
+ printf("%s\n", buf1);
+
+
+
+ return 0;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff