TUHKMgen/UHKM/RunHadronSourceHISTO.cxx

   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(Px1*Px1+Py1*Py1),1./sqrt(Px1*Px1+Py1*Py1));
 209          }
 210
 211       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));
 212       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));
 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(Px1*Px1+Py1*Py1);
 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.5*TMath::Log((sqrt(Px1*Px1+Py1*Py1+Pz1*Pz1)+Pz1)/(sqrt(Px1*Px1+Py1*Py1+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 }