[u/mrichter/AliRoot.git] / TUHKMgen / UHKM / RunHadronSourceHISTO.cxx

/*******************************************************************************
 *                                                                             *
 *    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;
}
Commit	Line	Data
b1c2e580	1	/*******************************************************************************
	2	* *
	3	* HYDJET++ , event generator under the ROOT FRAMEWORK for simulation of *
	4	* relativistic heavy ion AA collisions as the superposition of soft, *
	5	* hydro-type state and hard, multi-parton state. *
	6	* *
	7	* The main routine is written in the object-oriented C++ language *
	8	* under the ROOT environment. The hard, multi-partonic part of *
	9	* HYDJET++ event is identical to the hard part of Fortran-written *
	10	* HYDJET (PYTHIA6.4xx + PYQUEN1.5) and is included in the generator *
	11	* structure as the separate directory. The soft part of HYDJET++ *
	12	* event represents the "thermal" hadronic state obtained with the *
	13	* parameterization Bjorken-like of freeze-out hypersurface and *
	14	* includes longitudinal, radial and elliptic flow effects and *
	15	* decays of hadronic resonances. The corresponding fast *
	16	* Monte-Carlo simulation procedure (C++ code) FAST MC is adapted. *
	17	* -------------------------------------------------------------- *
	18	* Web-page: *
	19	* http://cern.ch/lokhtin/hydjet++ *
	20	* -------------------------------------------------------------- *
	21	* * *
	22	* *
	23	* This program is a free software; you can use and redistribute it freely. *
	24	* Any publication of results obtained using this code must reference *
	25	* *
	26	* *
	27	* *
	28	* Main reference for HYDJET++: *
	29	* I.P. Lokhtin, L.V. Malinina, S.V. Petrushanko, A.M. Snigirev, *
	30	* I. Arsene, K. Tywoniuk, submitted to Comp. Phys. Comm. *
	31	* *
	32	* Reference for HYDJET and PYQUEN: *
	33	* I.P. Lokhtin, A.M. Snigirev, Eur. Phys. J. C 46 (2006) 211; *
	34	* http://cern.ch/lokhtin/hydro/hydjet.html *
	35	* http://cern.ch/lokhtin/pyquen. *
	36	* *
	37	* Reference for PYTHIA6.4: *
	38	* T.Sjostrand, S. Mrenna and P. Skands, JHEP05 (2006) 026; *
	39	* http://home.thep.lu.se/~torbjorn/Pythia.html. *
	40	* *
	41	* References for FAST MC: *
	42	* N.S. Amelin, R. Lednicky, T.A. Pocheptsov, I.P. Lokhtin, *
	43	* L.V. Malinina, A.M. Snigirev, Iu.A. Karpenko and Yu.M. Sinyukov, *
	44	* Phys. Rev. C 74 (2006) 064901; *
	45	* N.S. Amelin, I. Arsene, L. Bravina, Iu.A. Karpenko, R. Lednicky, *
	46	* I.P. Lokhtin, L.V. Malinina, A.M. Snigirev and Yu.M. Sinyukov, *
	47	* Phys. Rev. C 77 (2008) 014903; *
	48	* http://uhkm.jinr.ru. *
	49	* *
	50	* Reference for nuclear shadowing model: *
	51	* K. Tywoniuk, I.C. Arsene, L. Bravina, A. Kaidalov and *
	52	* E. Zabrodin, Phys. Lett. B 657 (2007) 170. *
	53	* *
	54	* version 2.0: *
	55	* *
	56	* Igor Lokhtin, SINP MSU, Moscow, RU *
	57	* e-mail: Igor.Lokhtin@cern.ch *
	58	* *
	59	* Ludmila Malinina, SINP MSU, Moscow, RU *
	60	* e-mail: malinina@lav01.sinp.msu.ru *
	61	* *
	62	*******************************************************************************/
	63	#include <iostream>
	64	#include <fstream>
65	#include <vector>
66	#include <time.h>
67
68	#include <TNtuple.h>
69	#include <TError.h>
70	#include <TTree.h>
71	#include <TH1D.h>
72	#include <TFile.h>
73
74	#include "InitialState.h"
75	#include "InitialStateHydjet.h"
76
77
78	#include <TRandom.h>
79
80	#include "Particle.h"
81	//#include "HYJET_COMMONS.h"
82	//extern SERVICECommon SERVICE;
83
84
85	//Main program:
86	//reads input parameters from file "RunInputBjorken" or "RunInputHubble";
87	//calculates particle densities and average initial multiplicities and writes them
88	//in output file "multiplicities.txt";
89	//creates trees (tree with direct hadrons and hadrons after resonance decays)
90	//with space-time and momentum-energy information of produced hadrons;
91	//writes trees in file "RunOutput.root".
92
93	Int_t main() {
94
95	clock_t start;
96	start = clock();
97
98
99	//new
100	time_t now;
101	struct tm *ts;
102	char buf[80];
103
104	// Get the current time
105	time(&now);
106
107	// Format and print the time, "ddd yyyy-mm-dd hh:mm:ss zzz"
108	ts = localtime(&now);
109	strftime(buf, sizeof(buf), "%a %Y-%m-%d %H:%M:%S %Z", ts);
110	printf("%s\n", buf);
111
112
113
114
115	TFile *outputFile=new TFile("RunOutput.root", "RECREATE");
116
117	//SET MAXIMAl VALUE OF PARTICLE MULTIPLICITY!!!
118	const Int_t kMax = 500000;
119	//define hadron number
120	Int_t ntot;
121	//define event number
122	Int_t nev;
123	//define hadron characteristic vectors
124	std::vector<Int_t> pdg(kMax); //pdg encodings
125	std::vector<Int_t> Mpdg(kMax);//pdg encodings for mother hadrons
126	std::vector<Int_t> type(kMax);//type: 0-from hydro or decay, 1111 from jets
127	std::vector<Float_t> Px(kMax);//x-hadron momentum component,[GeV/c]
128	std::vector<Float_t> Py(kMax);//y-hadron momentum component,[GeV/c]
129	std::vector<Float_t> Pz(kMax);//z-hadron momentum component,[GeV/c]
130	std::vector<Float_t> E(kMax); //hadron total energy,[GeV]
131	std::vector<Float_t> X(kMax);//x-hadron coordinate component,[fm]
132	std::vector<Float_t> Y(kMax);//y-hadron coordinate component,[fm]
133	std::vector<Float_t> Z(kMax);//z-hadron coordinate component,[fm]
134	std::vector<Float_t> T(kMax);//hadron time,[fm/c]
135
136	TH1D *hpt1 = new TH1D("hpt1", "hpt1", 100, 0., 20.);
137	TH1D *hpt1j = new TH1D("hpt1j", "hpt1j", 100, 0., 20.);
138	TH1D *hpt1h = new TH1D("hpt1h", "hpt1h", 100, 0., 20.);
139
140	TH1D *hv2 = new TH1D("hv2", "hv2", 100, 0.0, 10.);
141	TH1D *hv0 = new TH1D("hv0", "hv0", 100, 0.0, 10.);
142
143	TH1D *hy = new TH1D("hy", "hy", 51, -5.1, 5.1);
144	TH1D *hyjets = new TH1D("hyjets", "hyjets", 51, -5.1, 5.1);
145	TH1D *hyhydro = new TH1D("hyhydro", "hyhydro", 51, -5.1, 5.1);
146
147
148	double pdg1, Mpdg1, Px1, Py1, E1, Z1, Pz1, pt, phi, v2, eta;
149	int type1;
150
151	InitialState *FASTMC;
152
153	FASTMC = new InitialStateHydjet();
154
155	if(!FASTMC->ReadParams()) {
156	Error("RunHadronSource::main", "No initial model parameters found!!\n");
157	return 0;
158	}
159
160
161	if(!FASTMC->MultIni()) {
162	Error("RunHadronSource::main", "Initial multiplicities are zero!!\n");
163	return 0;
164	}
165
166	ParticleAllocator allocator;
167	List_t source;
168	List_t secondaries;
169	std::cout << "Generating " << FASTMC->GetNev() << " events" << std::endl;
170	std::cout << "Starting the event loop" << std::endl;
171
172
173	// Loop over events
174	for(Int_t ev = 0; ev < FASTMC->GetNev(); ++ev) {
175	nev = ev;
176	// Initialize the source
177	FASTMC->Initialize(source, allocator);
178	if(source.empty()) {
179	Error("RunHadronSource::main", "Source is not initialized!!");
180	//return 0;
181	continue;
182	}
183
184	// Run the decays //fDecay
185	if(FASTMC->GetTime() >= 0.)
186	FASTMC->Evolve(source, secondaries, allocator, FASTMC->GetWeakDecayLimit());
187
188	std::cout << "event #" << ev << "\r" << std::flush;
189	// npart = 0;
190	LPIT_t it;
191	LPIT_t e;
192
193	// Fill the decayed tree
194	// npart = 0;
195
196	for(it = secondaries.begin(), e = secondaries.end(); it != e; ++it) {
197	TVector3 pos(it->Pos().Vect());
198	TVector3 mom(it->Mom().Vect());
199	Float_t m1 = it->TableMass();
200	pdg1 = it->Encoding();
201	Mpdg1 = it->GetLastMotherPdg();
202	Px1 = mom[0];
203	Py1 = mom[1];
204	Pz1 = mom[2];
205	E1 = TMath::Sqrt(mom.Mag2() + m1*m1);
206	type1 = it->GetType();
207	if(pdg1==211 && abs(0.5*log((E1+Pz1)/(E1-Pz1)))<1.) {
208	hpt1->Fill(sqrt(Px1Px1+Py1Py1),1./sqrt(Px1Px1+Py1Py1));
209	}
210
211	if(pdg1==211 && abs(0.5log((E1+Pz1)/(E1-Pz1)))<1. && type1==0) hpt1h->Fill(sqrt(Px1Px1+Py1Py1),1./sqrt(Px1Px1+Py1*Py1));
212	if(pdg1==211 && abs(0.5log((E1+Pz1)/(E1-Pz1)))<1. && type1==1)hpt1j->Fill(sqrt(Px1Px1+Py1Py1),1./sqrt(Px1Px1+Py1*Py1));
213
214	if(((abs(pdg1)==211)\|\|(abs(pdg1)==321)\|\|(abs(pdg1)==2212))
215	&& (abs(0.5*log((E1+Pz1)/(E1-Pz1)))<1.0)){
216	pt = TMath::Sqrt(Px1Px1+Py1Py1);
217	phi = TMath::ATan2(Py1,Px1);
218	v2 = TMath::Cos(2*phi);
219	hv2->Fill(pt,v2);
220	hv0->Fill(pt,1.);
221	}
222
223	if((abs(pdg1)==211)\|\|(abs(pdg1)==321)\|\|(abs(pdg1)==2212)){
224	eta=0.5TMath::Log((sqrt(Px1Px1+Py1Py1+Pz1Pz1)+Pz1)/(sqrt(Px1Px1+Py1Py1+Pz1*Pz1)-Pz1));
225	if(type1==1)hyjets->Fill(eta);
226	if(type1==0)hyhydro->Fill(eta);
227	hy->Fill(eta);
228	}
229
230	// npar++;
231	// if(npart > kMax)
232	// Error("in main:", "npart is too large %d", npart);
233
234
235	}
236
237	allocator.FreeList(source);
238	allocator.FreeList(secondaries);
239
240
241	}
242
243	hpt1->Write();
244	hpt1h->Write();
245	hpt1j->Write();
246	hv2->Write();
247	hv0->Write();
248	hyhydro->Write();
249	hyjets->Write();
250	hy->Write();
251
252	clock_t stop;
253	stop = clock();
254	std::cout << "*********************************************" << std::endl;
255	std::cout << "Execution time: " << (stop - start)/CLOCKS_PER_SEC << " seconds" << std::endl;
256	std::cout << "*********************************************" << std::endl;
257
258
259	//new
260	time_t now1;
261	struct tm *ts1;
262	char buf1[80];
263
264	// Get the current time
265	time(&now1);
266
267	// Format and print the time, "ddd yyyy-mm-dd hh:mm:ss zzz"
268	ts1 = localtime(&now1);
269	strftime(buf1, sizeof(buf1), "%a %Y-%m-%d %H:%M:%S %Z", ts1);
270	printf("%s\n", buf1);
271
272
273
274	return 0;
275	}