IAP GITLAB

Commit c8b05187 authored by Colin Baus's avatar Colin Baus

option 2 (as discussed via mail) for sibyll boosting of final particle stack....

option 2 (as discussed via mail) for sibyll boosting of final particle stack. off-shell particles cannot use true particle masses in utlob5

git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/trunk@4707 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent 5e1d0505
......@@ -43,23 +43,23 @@ int main (int argc, char **argv)
{
//-------------------SET UP DATA
TFile* theOutFile;
string outFileName("new_histogram_file.root");
string outFileName ("new_histogram_file.root");
cout << " ! Opening output file: " << outFileName << endl;
theOutFile = new TFile(outFileName.c_str(),"RECREATE");
theOutFile = new TFile (outFileName.c_str(),"RECREATE");
vector<string> filesModel1;
filesModel1.push_back("files/test.hepmc"); //add your files here with additional push_back()
filesModel1.push_back ("files/test.hepmc"); //add your files here with additional push_back()
DataManager data;
data.SetFiles(filesModel1); //for more models, loop over models and call SetFiles each time
theOutFile->mkdir("model1");
data.SetFiles (filesModel1); //for more models, loop over models and call SetFiles each time
theOutFile->mkdir ("model1");
//------------------SET UP HISTOGRAMS
TH1D* examplehist = new TH1D("dNdeta",";#eta;dN/d#eta",21,-10,10);
TH1D* exampleHist = new TH1D ("dNdeta",";#eta;dN/d#eta",21,-10,10);
//-------------------EVENT LOOP
int nEvts = 0;
while (data.GetNextEvent())
while (data.GetNextEvent ())
{
++nEvts;
//HepMC::HeavyIon* heavyIonInfo = NULL; //helpful to get cross section: heavyIonInfo->sigma_inel_NN ()
......@@ -71,7 +71,7 @@ int main (int argc, char **argv)
{
HepMC::GenParticle* p = (*par);
if (p->status () != 1) continue; //get final state particles status==2 would be decayed particles, status == 4 is beam particles
if (p->status () != 1) continue; //get final state particles. status == 2 are decayed particles, status == 4 is beam particles
//HepMC::GenVertex* parent_vertex = p->production_vertex();
//const int id = p->pdg_id ();
......@@ -83,17 +83,21 @@ int main (int argc, char **argv)
//-------------------EVENT SELECTION
//-------------------FILL HISTOGRAMS WITH PER PARTICLE VARIABLES
examplehist->Fill(eta);
exampleHist->Fill (eta);
}//PARTICLE LOOP
//-------------------FILL HISTOGRAMS WITH PER EVENT VARIABLES
}//EVENT LOOP
//---------------FINALISE HISTOGRAMS
examplehist->Scale(1./nEvts,"width");
exampleHist->Scale (1. / nEvts, "width");
std::cout << " ! Writing output file " << outFileName << std::endl;
//----------------Closing Files
std::cout << " ! Writing output file: " << outFileName << std::endl;
theOutFile->Write();
std::cout << " ! Closing output file " << outFileName << std::endl;
delete exampleHist;
exampleHist = 0;
std::cout << " ! Closing output file: " << outFileName << std::endl;
theOutFile->Close();
delete theOutFile;
theOutFile = 0;
return 0;
}
......@@ -34,7 +34,8 @@ c Input values
double precision decms,e1,e2
double precision ycm2det
common/boostvars/ycm2det
logical doBoost
common/boostvars/ycm2det,doBoost
common/producetab/ producetables !used to link CRMC
logical producetables !with EPOS and QII
......@@ -61,6 +62,7 @@ c Calculations of energy of the center-of-mass in the detector frame
e2=dsqrt(dble(m2)**2+ptarg**2)
decms=dsqrt((e1+e2)**2-(pproj+ptarg)**2)
c Later a rapidity boost back into the detector system will be performed
doBoost = .true.
c ycm2det defines this rapidity
if (((e1+e2)-(pproj+ptarg)) .le. 0d0) then
ycm2det=1d99
......@@ -136,8 +138,10 @@ c Output quantities
integer outstat(*)
double precision boostvec1,boostvec2,boostvec3,boostvec4,boostvec5
double precision mass
double precision ycm2det
common/boostvars/ycm2det !is this common block needed?
logical doBoost
common/boostvars/ycm2det,doBoost
integer i!,k
......@@ -165,17 +169,24 @@ c stop
c define vec to boost from cm. to cms frame
boostvec1=0d0
boostvec2=0d0
boostvec3=dsinh(dble(ycm2det))
boostvec4=dcosh(dble(ycm2det))
boostvec3=dsinh(ycm2det)
boostvec4=dcosh(ycm2det)
boostvec5=1d0
c write(*,*)nevhep,nhep,boostvec3,boostvec4,ycm2det
do i=1,nhep
c boost output to cms frame
if(doBoost.eqv..true..and.ycm2det.ne.0d0) then
if(model.eq.6) then ! sibyll needs calculated mass because itlob5 does not work with off-shell particles
mass = dsqrt(phep(4,i)**2
+ - phep(1,i)**2 - phep(2,i)**2 - phep(3,i)**2)
else
mass = phep(5,i)
endif
call utlob2(-1,boostvec1,boostvec2,boostvec3
+ ,boostvec4,boostvec5
+ ,vhep(1,i),vhep(2,i),vhep(3,i),vhep(4,i),-99)
call utlob5dbl(-ycm2det
+ ,phep(1,i), phep(2,i), phep(3,i), phep(4,i), phep(5,i))
+ ,phep(1,i), phep(2,i), phep(3,i), phep(4,i), mass)
outpart(i)=idhep(i)
outpx(i)=phep(1,i)
outpy(i)=phep(2,i)
......@@ -183,6 +194,7 @@ c boost output to cms frame
oute(i)=phep(4,i)
outm(i)=phep(5,i)
outstat(i)=isthep(i)
endif
c write(*,'(4x,i6,1x,4(e12.6,1x))')idhep(i),(vhep(k,i),k=1,4)
c write(*,'(i5,3x,i2,2x,2i5,2x,2i5)')i,isthep(i)
c * ,jmohep(1,i),jmohep(2,i),jdahep(1,i),jdahep(2,i)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment