123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981 |
- // -------------------------------------------------------------------------
- // ----- MpdLAQGSMGenerator source file -----
- // -------------------------------------------------------------------------
- /** MpdLAQGSMGenerator
- *@author Elena Litvinenko <litvin@nf.jinr.ru>
- *@version 15.02.2016
- *
- ** The MpdLAQGSMGenerator uses the ASCII input
- ** provided by K.Gudima LAQGSM event generator.
- **/
- #include "MpdLAQGSMGenerator.h"
- #include "MpdHypYPtGenerator.h" //AZ
- #include "MpdStack.h"
- #include "FairMCEventHeader.h"
- #include "FairPrimaryGenerator.h"
- #include "FairIon.h"
- #include "FairRunSim.h"
- #include "FairParticle.h"
- #include "TDatabasePDG.h"
- #include "TParticlePDG.h"
- #include "TParticle.h"
- #include "TSystem.h"
- #include "TParticle.h"
- #include "TVirtualMC.h"
- #include <iostream>
- using std::cout;
- using std::endl;
- using std::map;
- const Double_t kProtonMass = 0.938272013;
- const Double_t kNeutronMass = 0.93957;
- // ----- Default constructor ------------------------------------------
- MpdLAQGSMGenerator::MpdLAQGSMGenerator() :
- fInputFile(NULL),
- fGZInputFile(NULL)
- {}
- // ------------------------------------------------------------------------
- // ----- Standard constructor -----------------------------------------
- MpdLAQGSMGenerator::MpdLAQGSMGenerator(const char* fileName, const Bool_t use_collider_system, Int_t QGSM_format_ID,Int_t Max_Event_Number ) :
- FairGenerator(),
- fIonMap(),
- fInputFile(NULL),
- fGZInputFile(NULL)
- {
- //AZ memset(la_tab,0,sizeof(la_tab));
- fPDG=TDatabasePDG::Instance();
- fFileName = fileName;
- cout << "-I- MpdLAQGSMGenerator: Opening input file " << fileName << endl;
- TString sFileName = fileName;
- if (sFileName.Contains(".gz"))
- fGZ_input=1;
- else
- fGZ_input=0;
- if (fGZ_input)
- fGZInputFile = gzopen(fFileName, "rb");
- else
- fInputFile = fopen(fFileName, "r");
- if (( ! fInputFile ) && ( ! fGZInputFile ))
- Fatal("MpdLAQGSMGenerator","Cannot open input file.");
- fUseColliderSystem = use_collider_system;
- fQGSM_format_ID=QGSM_format_ID;
- Init(); //AZ
- cout << "-I- MpdLAQGSMGenerator: Looking for ions..." << endl;
- //AZ Int_t nIons = RegisterIons();
- // Int_t nIons = RegisterIons1();
- // Int_t nIons = RegisterIons(); // MG
- Int_t nIons = RegisterIons(Max_Event_Number); // EL
- cout << "-I- MpdLAQGSMGenerator: " << nIons << " ions registered."
- << endl;
- CloseInput();
- cout << "-I- MpdLAQGSMGenerator: Reopening input file " << fileName
- << endl;
- if (fGZ_input)
- fGZInputFile= gzopen(fFileName, "rb");
- else
- fInputFile = fopen(fFileName, "r");
-
- }
- // ------------------------------------------------------------------------
- // ----- Destructor ---------------------------------------------------
- MpdLAQGSMGenerator::~MpdLAQGSMGenerator()
- {
- CloseInput();
- Int_t nPart = fLa_tab.size();
- for (Int_t i = 0; i < nPart; ++i) delete fLa_tab[i];
- fLa_tab.clear();
- }
- // ------------------------------------------------------------------------
- // --------------Public method Init (fill list of known light particles)---
- void MpdLAQGSMGenerator::Init (const char *light_particles_filename)
- {
- TString fname, dd;
- if (!light_particles_filename) {
- dd = getenv("VMCWORKDIR");
- if (gSystem->OpenDirectory(dd+"/mpdgenerators"))
- dd += "/mpdgenerators";
- else {
- if (gSystem->OpenDirectory(dd+"/generators"))
- dd += "/generators";
- else
- dd = gSystem->WorkingDirectory();
- }
- fname = dd+ "/MpdLAQGSM_light_particles.dat";
- }
- else
- fname = light_particles_filename;
- Int_t i=0, i_val, i_max=84;
- Float_t a_val;
- char ss[300];
- FILE *fInputFile_light = fopen(fname.Data(), "r");
- //AZ while (!(feof(fInputFile_light))&&(i<i_max)) {
- while (!(feof(fInputFile_light))) {
- la_tab_t *lt = new la_tab_t;
- /*AZ
- fscanf(fInputFile_light, "%d", &i_val); la_tab[i].pdg=i_val;
- fscanf(fInputFile_light, "%d", &i_val); la_tab[i].Z=i_val;
- fscanf(fInputFile_light, "%d", &i_val); la_tab[i].lepton=i_val;
- fscanf(fInputFile_light, "%d", &i_val); la_tab[i].strange=i_val;
- fscanf(fInputFile_light, "%d", &i_val); la_tab[i].A=i_val;
- fscanf(fInputFile_light, "%f", &a_val); la_tab[i].mass=a_val;
- fscanf(fInputFile_light, "%s", ss); strncpy(la_tab[i].name,ss,10);
- i++;
- */
- i = fscanf(fInputFile_light, "%d", &i_val); lt->pdg=i_val;
- i = fscanf(fInputFile_light, "%d", &i_val); lt->Z=i_val;
- i = fscanf(fInputFile_light, "%d", &i_val); lt->lepton=i_val;
- i = fscanf(fInputFile_light, "%d", &i_val); lt->strange=i_val;
- i = fscanf(fInputFile_light, "%d", &i_val); lt->A=i_val;
- i = fscanf(fInputFile_light, "%f", &a_val); lt->mass=a_val;
- i = fscanf(fInputFile_light, "%s", ss); strncpy(lt->name,ss,10);
- fLa_tab.push_back(lt);
- if(!i) cerr<<__FILE__<<__func__<<" something went wrong in reading input.";
- }
- fclose(fInputFile_light);
- }
- // ------------------------------------------------------------------------
- Bool_t MpdLAQGSMGenerator::general_fgets (char *ss, Int_t nn, FILE* p)
- {
- Bool_t finished=0;
- if (fGZ_input) {
- gzgets(fGZInputFile,ss,nn);
- finished = (gzeof(fGZInputFile));
- }
- else {
- if(!fgets(ss,nn,fInputFile))
- cerr<<__FILE__<<__func__<<" something went wrong in reading input.";
- finished = ( feof(fInputFile) );
- }
- return finished;
- }
- Bool_t MpdLAQGSMGenerator::GetEventHeader(char *ss)
- {
- Bool_t finished=0;
- // if (fGZ_input) {
- // gzgets(fGZInputFile,ss,250);
- // finished = (gzeof(fGZInputFile));
- // }
- // else {
- // fgets(ss,250,fInputFile);
- // finished = ( feof(fInputFile) );
- // }
- finished = general_fgets (ss,250,fInputFile);
- // If end of input file is reached : close it and abort run
- if ( finished ) {
- cout << "-I- MpdLAQGSMGenerator: End of input file reached " << endl;
- CloseInput();
- return kFALSE;
- }
- Bool_t wrong_file=kFALSE;
- TString tss(ss);
- if (tss.Contains("QGSM")) { // 0
- general_fgets(ss,250,fInputFile);
- tss=ss;
- Int_t lines=0;
- while (!(tss.Contains("particles, b,bx,by"))) {
- if (tss.Contains("QGSM format ID")) {
- sscanf(ss,"QGSM format ID=%d",&fQGSM_format_ID);
- cout << "QGSM format ID " << fQGSM_format_ID << endl;
- if (fQGSM_format_ID==2) { // correction of incorrect format_ID in some data files
- fpos_t file_loc;
- z_off_t gz_file_loc;
- if (fGZ_input)
- gz_file_loc = gztell(fGZInputFile);
- else
- fgetpos(fInputFile, &file_loc);
- for (int k=0;k<5;k++) general_fgets(ss,250,fInputFile);
- if (strlen(ss)>90)
- fQGSM_format_ID=3;
- cout << "QGSM format ID now " << fQGSM_format_ID << endl;
- if (fGZ_input)
- gzseek(fGZInputFile,gz_file_loc,SEEK_SET);
- else
- fsetpos(fInputFile, &file_loc);
- }
- }
- general_fgets(ss,250,fInputFile);
- tss=ss;
- lines++;
- if ((fQGSM_format_ID>=2)&&(lines>=4))
- return kTRUE;
- if (lines>25) {
- wrong_file=kTRUE;
- break;
- }
- }
- }
- else {
- if (fQGSM_format_ID>=2)
- return kTRUE;
- if (!tss.Contains("particles, b,bx,by")) {
- general_fgets(ss,250,fInputFile);
- tss=ss;
- }
- }
- if (wrong_file) {
- cout << "-E- MpdLAQGSMGenerator: Wrong file header! " << endl;
- CloseInput();
- return kFALSE;
- }
-
- if (!tss.Contains("event")) { // 3
- if (fQGSM_format_ID>=2)
- return kTRUE;
- else {
- cout << "-E- MpdLAQGSMGenerator: Wrong event header!" << endl;
- CloseInput();
- return kFALSE;
- }
- }
- char tmp[250];
- switch (fQGSM_format_ID) {
- case 0:
- case 1:
- general_fgets(tmp,250,fInputFile);
- break;
- case 2:
- case 3:
- general_fgets(ss,250,fInputFile);
- general_fgets(ss,250,fInputFile);
- break;
- case 4:
- general_fgets(ss,250,fInputFile);
- general_fgets(ss,250,fInputFile);
- break;
- default:
- break;
- }
- return kTRUE;
- }
- // ------------------------------------------------------------------------
- //AZ ----- Public method SkipEvents --------------------------------------
- Bool_t MpdLAQGSMGenerator::SkipEvents(Int_t count) {
- if (count<=0) return kTRUE;
- // Check for input file
- if (fGZ_input) {
- if ( ! fGZInputFile ) {
- cout << "-E- MpdLAQGSMGenerator: Input file not open!" << endl;
- return kFALSE;
- }
- }
- else {
- if ( ! fInputFile ) {
- cout << "-E- MpdLAQGSMGenerator: Input file not open!" << endl;
- return kFALSE;
- }
- }
- Int_t eventId = 0, nTracks = 0;
- for(Int_t ii=0; ii<count;ii++) {
- char ss[252];
- if (!GetEventHeader(ss))
- return kFALSE;
- sscanf(ss," %d %d", &eventId, &nTracks);
- //cout << ii << " " << eventId << endl;
- for (Int_t itrack=0; itrack<nTracks; itrack++) {
- general_fgets(ss,250,fInputFile);
- }
- }
- return kTRUE;
- }
- // ------------------------------------------------------------------------
- // ----- Public method ReadEvent --------------------------------------
- Bool_t MpdLAQGSMGenerator::ReadEvent(FairPrimaryGenerator* primGen) {
- // Check for input file
- if (fGZ_input) {
- if ( ! fGZInputFile ) {
- cout << "-E- MpdLAQGSMGenerator: Input file not open!" << endl;
- return kFALSE;
- }
- }
- else {
- if ( ! fInputFile ) {
- cout << "-E- MpdLAQGSMGenerator: Input file not open!" << endl;
- return kFALSE;
- }
- }
- // ---> Check for primary generator
- if ( ! primGen ) {
- cout << "-E- MpdLAQGSMGenerator::ReadEvent: "
- << "No PrimaryGenerator!" << endl;
- return kFALSE;
- }
- //AZ if (la_tab[0].pdg<=0)
- if (fLa_tab[0]->pdg<=0)
- Init();
- // Define event variable to be read from file
- Int_t eventId = 0;
- Int_t nTracks = 0;
- Float_t b = 0., bx=0, by=0;
- // Define track variables to be read from file
- Int_t
- iCharge = 0,
- iLeptonic = 0,
- iStrange = 0,
- iBarionic=0, iCode=0, io1=0, io2=0, io3=0 //,
- //A=0,
- //Z=0
- ;
- Float_t
- px = 0.,
- py = 0.,
- pz = 0.,
- pz1 = 0,
- pza1 = 0,
- // Ekin = 0,
- // Ekin1 = 0,
- mass = 0, pol0 = 0, pol1 = 0;
- Int_t PDG=0;
- char ss[252];
- char ionName[30];
- if (!GetEventHeader(ss))
- return kFALSE;
- //Double_t projA,projZ,projEex, targA, targZ, targEex;
- Float_t projA,projZ,projEex, targA, targZ, targEex;
- //cout <<"ss " <<ss <<endl;
- int bpos;
- sscanf(ss," %d %d %n", &eventId, &nTracks, &bpos);
- //cout <<"eventId " <<eventId <<" " <<nTracks <<" " <<fQGSM_format_ID <<endl;
- if (!fQGSM_format_ID)
- sscanf(&(ss[bpos]),"%g %g %g %g %g %g %g %g %g", &b, &bx, &by, &projA,&projZ,&projEex, &targA, &targZ, &targEex );
- else
- sscanf(&(ss[bpos]),"%g %g %g", &b, &bx, &by);
- cout << "-I- MpdLAQGSMGenerator::ReadEvent: Event " << eventId
- << ", b = " << b << " fm, multiplicity " << nTracks << endl;
- TVector2 bb; // MG
- bb.Set(bx,by); // MG
- // Set event id and impact parameter in MCEvent if not yet done
- FairMCEventHeader* event = primGen->GetEvent();
- if ( event && (! event->IsSet()) ) {
- event->SetEventID(eventId-1);
- event->SetB(b);
- event->SetRotZ(bb.Phi()); // MG
- event->MarkSet(kTRUE);
- }
- /*
- //AZ - Loop over tracks to find hypernuclei
- fpos_t position;
- fgetpos (fInputFile, &position); // save position
- Int_t iHyperOk = 0;
- for (Int_t itrack = 0; itrack < nTracks; itrack++) {
- fgets(ss,250,fInputFile);
- sscanf(ss," %d %d %d %d %d", &iCharge,&iLeptonic,&iStrange,&iBarionic, &iCode);
- if (iStrange == -1 && iCharge == 1 && iBarionic == 3) {
- // select only H3L
- iHyperOk = 1;
- break;
- }
- // Select H4LL
- //if (iStrange == -2 && iCharge == 1 && iBarionic == 4) {
- //iHyperOk = 1;
- //break;
- //}
- }
- if (iHyperOk) fsetpos (fInputFile, &position);
- //else nTracks = 0; // do not store any track
- else {
- // Add particle manually
- TParticlePDG* part = fPDG->GetParticle("H3L");
- Int_t pdgType = part->PdgCode();
- static MpdHypYPtGenerator *hypGen = 0x0;
- if (hypGen == 0x0) {
- hypGen = new MpdHypYPtGenerator(pdgType,1);
- hypGen->SetDistributionPt(0.163); // temperature - H3L: NICA @ 5 GeV
- hypGen->SetDistributionY(0.0, 0.373); // y0, Sigma_y - H3L: NICA @ 5 GeV
- hypGen->Init();
- }
- hypGen->ReadEvent(primGen);
- fsetpos (fInputFile, &position);
- }
- // AZ
- */
- // Loop over tracks in the current event
- for (Int_t itrack=0; itrack<nTracks; itrack++) {
- general_fgets(ss,250,fInputFile);
- if (fQGSM_format_ID<3)
- sscanf(ss," %d %d %d %d %d", &iCharge,&iLeptonic,&iStrange,&iBarionic, &iCode);
- else{
- if(fQGSM_format_ID==4){
- //cout << ss << endl;
- sscanf(ss," %d %d %d %d %d %g %g %g %g %g %g", &iCharge,&iLeptonic,&iStrange,&iBarionic, &PDG, &px, &py, &pz, &pz1, &pza1, &mass);
- //printf("!!!!!!!!!!!!!!!!!!!!!reading: %d %d %d %d %d %g %g %g %g %g %g\n", iCharge,iLeptonic,iStrange,iBarionic, PDG, px, py, pz, pz1, pza1, mass);
- //primGen->AddTrack(PDG, px, py, pz, 0., 0., 0.);
- if (fUseColliderSystem)
- primGen->AddTrack(PDG, px, py, pz, 0., 0., 0.);
- else
- primGen->AddTrack(PDG, px, py, pz1, 0., 0., 0.);
- }
- else
- sscanf(ss," %d %d %d %d %d %d %d", &iCharge,&iLeptonic,&iStrange,&iBarionic, &io1, &io2, &io3);
- }
- if(fQGSM_format_ID!=4){
- Int_t p_num = 21;
- if (fQGSM_format_ID==3)
- p_num=36;
-
- if (!fQGSM_format_ID)
- sscanf(&(ss[p_num]),"%g%g%g%g%g",&px,&py,&pz,&pz1,&mass);
- else if (strlen(ss) > 105)
- sscanf(&(ss[p_num]),"%g%g%g%g%g%g%g",&px,&py,&pz,&pz1,&pol0,&pol1,&mass); //AZ - polarization data
- else
- sscanf(&(ss[p_num]),"%g%g%g%g%g%g",&px,&py,&pz,&pz1,&pza1,&mass);
-
-
- if (FindParticle (iCharge, iStrange, iLeptonic, iBarionic, mass, PDG, ionName)) {
-
- if (!fPDG->GetParticle(PDG)) {
- cout << "-W- MpdLAQGSMGenerator::ReadEvent: Cannot find " << PDG << " "
- << ionName << " in database!" << endl;
- continue;
- }
-
- // Give track to PrimaryGenerator
- if (fUseColliderSystem)
- primGen->AddTrack(PDG, px, py, pz, 0., 0., 0.);
- else
- primGen->AddTrack(PDG, px, py, pz1, 0., 0., 0.);
- //AZ - add polarization for light hyperons (pdg == 31** and 32**)
- //*
- //if (PDG == 3122) {
- if (TMath::Abs(PDG/100) == 31 || TMath::Abs(PDG/100) == 32) {
- Int_t nTr = gMC->GetStack()->GetNtrack();
- TParticle *part = ((MpdStack*)gMC->GetStack())->GetParticle(nTr-1);
- Double_t pol[3] = {pol1,0,0}; // polarization after rescattering
- pol[1] = TMath::Sqrt (1.0 - pol[0] * pol[0]);
- part->SetPolarisation(pol[0], pol[1], pol[2]);
- }
- //*/
- }//if (FindParticle (iCharge
- }//if(fQGSM_format_ID!=4)
- }
-
- return kTRUE;
- }
- // ------------------------------------------------------------------------
- // ----- Private method CloseInput ------------------------------------
- void MpdLAQGSMGenerator::CloseInput() {
- if ( fInputFile ) {
- if (!fGZ_input)
- fclose(fInputFile);
- fInputFile = NULL;
- }
- if ( fGZInputFile ) {
- if (fGZ_input)
- gzclose(fGZInputFile);
- fGZInputFile = NULL;
- }
- }
- // ------------------------------------------------------------------------
- Int_t MpdLAQGSMGenerator::CreatePdgCode(Int_t Z, Int_t A, Int_t Strange,Int_t user)
- {
- // 10LZZZAAAI
- const Int_t kion =1000000000;
- const Int_t kion_strange = kion+10000000;
- const Int_t kspecial =50000000;
- if (Strange) {
- Int_t ss=abs(Strange);
- if (Z) {
- if (ss==1)
- return (kion_strange+Z*10000+A*10);
- else
- return (kion+10000000*ss+Z*10000+A*10);
- }
- else
- return (kspecial + A*10);
- }
- else {
- return (kion+Z*10000+A*10+user);
- }
- }
- Bool_t MpdLAQGSMGenerator::general_feof (void *p)
- {
- Bool_t finished=0;
- if (fGZ_input)
- finished = (gzeof(fGZInputFile));
- else
- finished = ( feof(fInputFile) );
- return finished;
- }
- // ----- Private method RegisterIons ----------------------------------
- Int_t MpdLAQGSMGenerator::RegisterIons(Int_t Max_Event_Number) {
- Int_t nIons = 0;
- // Int_t eventId, nTracks, iPid, iMass, iCharge;
- // Double_t pBeam, b, px, py, pz;
- Int_t eventId = 0;
- Int_t nTracks = 0;
- // Float_t b = 0., bx=0, by=0;
- // Define track variables to be read from file
- Int_t
- iCharge = 0,
- iLeptonic = 0,
- iStrange = 0,
- iBarionic=0,
- iCode =0;
- Float_t
- px = 0.,
- py = 0.,
- pz = 0.,
- pz1 = 0,
- pza1 = 0,
- // Ekin = 0,
- // Ekin1 = 0,
- mass = 0,
- excEnergy=0;
- char
- ss[252],
- buf_ionName[30];
- Int_t
- A=0, Z=0, Q=0;
- fIonMap.clear();
- if (!GetEventHeader(ss))
- return 0;
- Bool_t geant4_bug_case=1;
- Int_t PDG;
- TString ionName;
- while ( (!general_feof(fInputFile)) && ((!Max_Event_Number)||(eventId<Max_Event_Number)) ) {
- sscanf(ss," %d %d", &eventId, &nTracks);
- if ( general_feof(fInputFile) ) continue;
- for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
- general_fgets(ss,250,fInputFile);
- sscanf(ss," %d %d %d %d %d", &iCharge,&iLeptonic,&iStrange,&iBarionic,&iCode);
-
- // if (!fQGSM_format_ID)
- // sscanf(&(ss[21]),"%g%g%g%g%g",&px,&py,&pz,&pz1,&mass);
- // else
- // sscanf(&(ss[21]),"%g%g%g%g%g%g",&px,&py,&pz,&pz1,&pza1,&mass);
- Int_t p_num = 21;
- if (fQGSM_format_ID>2)
- p_num=36;
- if (!fQGSM_format_ID) {
- sscanf(&(ss[p_num]),"%g%g%g%g%g",&px,&py,&pz,&pz1,&mass);
- //if(iCharge>20)
- // cout << "format 1 " << mass << " " << fQGSM_format_ID << endl;
- }
- else{
- sscanf(&(ss[p_num]),"%g%g%g%g%g%g",&px,&py,&pz,&pz1,&pza1,&mass);
- //if(iCharge>20)
- // cout << "format 2 " << px << " " << py << " " << pz << " " << pz1 << " " << pza1 << " "<< mass << " " << fQGSM_format_ID << endl;
- }
- if (iBarionic>1){ // ion
- A = iBarionic;
- Z = iCharge;
- Q = iCharge;
- if (A>9)
- geant4_bug_case=0;
- //if(iCharge>20)
- // cout << "??????????????????1 " << iCharge << " " << iStrange << " " << iLeptonic << " " << iBarionic<< " " << mass << endl;
- if (FindParticle (iCharge, iStrange, iLeptonic, iBarionic, mass, PDG, buf_ionName)) {
- //if(iCharge>20)
- // cout << "??????????????????2 " << iCharge << " " << iStrange << " " << iLeptonic << " " << iBarionic<< " " << mass << endl;
- if (fPDG->GetParticle(PDG))
- continue;
-
- // cout << "!!!!!!!!!!!!!!!!!!!! " << PDG << endl;
- ionName = buf_ionName;
- if (!iStrange) { // normal ion
- if(Z>2 || (Z==2 && A>4)) { // MG
- if ( fIonMap.find(ionName) == fIonMap.end() ) { // new ion
- // excEnergy = fabs(mass - kProtonMass*Z -kNeutronMass*(iBarionic-Z));
- FairIon* ion = new FairIon(ionName, Z, A, Q) ; // , excEnergy,mass); // from MG
- fIonMap[ionName] = ion;
- nIons++;
- }
- }
- }
- else { // hyper
- if (CreateNucleus (Z, mass, PDG, buf_ionName)) {
- // FairParticle *pa=new FairParticle(buf_ionName,Z, iBarionic, abs(iStrange), mass,Z,0,2.6e-10);
- /*
- TParticle *tpa= new TParticle();
- if (tpa) {
- tpa->SetPdgCode(PDG);
- FairParticle *pa=new FairParticle(PDG,tpa);
- if (pa)
- (FairRunSim::Instance())->AddNewParticle(pa);
- }
- */
- if (Z) { //AZ
- FairParticle *pa = new FairParticle(PDG, buf_ionName, kPTHadron, mass, Z, 2.6e-10, "hadron");
- if (pa) FairRunSim::Instance()->AddNewParticle(pa);
- } else {
- FairParticle *pa = new FairParticle(PDG, buf_ionName, kPTNeutron, mass, Z, 2.6e-10, "special");
- if (pa) FairRunSim::Instance()->AddNewParticle(pa);
- } //AZ
- }
- }
- }
- } // ion
-
- } // track loop
- if (!GetEventHeader(ss))
- break;
- } // event loop
- if (geant4_bug_case) { // artificial heavy ion
- mass = 9.98;
- excEnergy = fabs(mass - (kProtonMass+kNeutronMass)*5);
- TString ionName_dummy="Ion_10_5";
- if ( fIonMap.find(ionName_dummy) == fIonMap.end() ) {
- FairIon* ion = new FairIon(ionName_dummy, 5,10, 5, excEnergy,mass);
- fIonMap[ionName_dummy] = ion;
- nIons++;
- }
- }
- FairRunSim* run = FairRunSim::Instance();
- map<TString, FairIon*>::iterator mapIt;
- for (mapIt=fIonMap.begin(); mapIt!=fIonMap.end(); ++mapIt) {
- FairIon* ion = (*mapIt).second;
- run->AddNewIon(ion);
- }
- return nIons;
- }
- // ----- Private method RegisterIons1 ----------------------------------
- Int_t MpdLAQGSMGenerator::RegisterIons1()
- {
- // Read file with particle data
- Int_t nIons = 0, nPart = fLa_tab.size();
- fIonMap.clear();
- Bool_t geant4_bug_case = 1;
- Int_t PDG;
- TString ionName;
- Float_t mass = 0.0, excEnergy = 0.0;
- char buf_ionName[30];
- for (Int_t i = 0; i < nPart; ++i) {
- Int_t Z = fLa_tab[i]->Z, iCharge = Z, Q = Z;
- Int_t iLeptonic = fLa_tab[i]->lepton;
- Int_t iStrange = fLa_tab[i]->strange;
- Int_t A = fLa_tab[i]->A, iBarionic = A;
- mass = fLa_tab[i]->mass;
- if (iBarionic > 1) {
- // ion
- //if (A > 9) geant4_bug_case = 0;
- if (FindParticle (iCharge, iStrange, iLeptonic, iBarionic, mass, PDG, buf_ionName)) {
- if (fPDG->GetParticle(PDG)) continue;
- ionName = buf_ionName;
- if (!iStrange) {
- // normal ion
- if ( fIonMap.find(ionName) == fIonMap.end() ) { // new ion
- excEnergy = fabs(mass - kProtonMass*Z -kNeutronMass*(iBarionic-Z));
- FairIon* ion = new FairIon(ionName, Z, A, Q, excEnergy,mass);
- fIonMap[ionName] = ion;
- nIons++;
- }
- } else {
- // hyper
- if (CreateNucleus (Z, mass, PDG, buf_ionName)) {
- // FairParticle *pa=new FairParticle(buf_ionName,Z, iBarionic, abs(iStrange), mass,Z,0,2.6e-10);
- if (Z) {
- FairParticle *pa = new FairParticle(PDG, buf_ionName, kPTHadron, mass, Z, 2.6e-10, "hadron");
- if (pa) FairRunSim::Instance()->AddNewParticle(pa);
- } else {
- FairParticle *pa = new FairParticle(PDG, buf_ionName, kPTNeutron, mass, Z, 2.6e-10, "special");
- if (pa) FairRunSim::Instance()->AddNewParticle(pa);
- }
- }
- }
- }
- } // if (iBarionic > 1)
-
- } // for (Int_t i = 0; i < nPart;
- if (geant4_bug_case) { // artificial heavy ion
- mass = 9.98;
- excEnergy = fabs(mass - (kProtonMass+kNeutronMass)*5);
- TString ionName_dummy="Ion_10_5";
- if ( fIonMap.find(ionName_dummy) == fIonMap.end() ) {
- FairIon* ion = new FairIon(ionName_dummy, 5,10, 5, excEnergy,mass);
- fIonMap[ionName_dummy] = ion;
- nIons++;
- }
- }
- FairRunSim* run = FairRunSim::Instance();
- map<TString, FairIon*>::iterator mapIt;
- for (mapIt=fIonMap.begin(); mapIt!=fIonMap.end(); ++mapIt) {
- FairIon* ion = (*mapIt).second;
- run->AddNewIon(ion);
- }
- return nIons;
- }
- // ------------------------------------------------------------------------
- Bool_t MpdLAQGSMGenerator::FindParticle (Int_t Z, Int_t strange, Int_t lepton, Int_t A, Float_t mass, Int_t &PDG, char name[11])
- {
- Int_t
- n_mass,
- result=0;
- if (mass<6) { // use table
- Int_t
- l_min=-1, // line number first
- l_max=-1; // line number last
- n_mass = (Int_t) (mass/0.1);
- switch (n_mass) {
- case 1: l_min=8; l_max=12; // MU+,MU-,PI0,PI+,PI-
- break;
- case 9: l_min=28; l_max=32; // P,AP,AN,N,ETAP
- break;
- case 4: l_min=13; l_max=18; // K+,K-,AK0,KL,KS,K0
- break;
- case 18: l_min=66; l_max=66; // d
- break;
- case 5: l_min=19; l_max=19; // ETA
- break;
- case 11: l_min=34; l_max=41; // AL,L,S+,AS+,AS0,S0,AS-,S-
- break;
- case 28: l_min=68; l_max=69; // He3,t
- break;
- case 37: l_min=71; l_max=71; // He4
- break;
- case 13: l_min=50; l_max=59; // AXI0,XI0,AXI-,XI-,AS*0,S*0,S*+,AS*+,S*-,AS*-
- break;
- case 29: l_min=70; l_max=70; // H3L
- break;
- case 22:
- if (!Z)
- { l_min=67; l_max=67;} // LL
- else
- { l_min=81; l_max=81;} // KsiNP
- break;
- case 39: l_min=72; l_max=73; // He4L,H4L
- break;
- case 48: l_min=75; l_max=75; // He5L
- break;
- case 16: l_min=64; l_max=65; // AOM-,OM-
- break;
- case 41: l_min=74; l_max=74; // H4LL
- break;
- case 50: l_min=76; l_max=77; // He5LL,H5LL
- break;
- case 59: l_min=78; l_max=78; // He6LL
- break;
- case 0: l_min=1; l_max=7; // GM,ANUM,ANUE,NUE,NUM,E+,E-
- break;
- case 7: l_min=20; l_max=23; // RHO+,RHO-,RHO0,OMEG
- break;
- case 8: l_min=24; l_max=27; // K*+,K*-,AK*0,K*0
- break;
- case 10: l_min=33; l_max=33; // PHI
- break;
- case 12: l_min=42; l_max=49; // DL++,ADL++,DL+,DL-,ADL0,DL0,ADL-,ADL+
- break;
- case 15: l_min=60; l_max=63; // AXI*0,XI*0,AXI*-,XI*-
- break;
- case 20: l_min=79; l_max=80; // LN, LP
- break;
- case 24: l_min=82; l_max=83; // KsiL0, KsiL1
- break;
- case 26: l_min=84; l_max=84; // KsiKsi
- break;
- default: // unknown particle
- break;
- }
- if (l_min>0){
- for (int i=l_min-1; i<l_max; i++ ) {
- /*AZ
- if (Z!=la_tab[i].Z) continue;
- if (strange!=la_tab[i].strange) continue;
- if (lepton!=la_tab[i].lepton) continue;
- if (A!=la_tab[i].A) continue;
- if (mass<(la_tab[i].mass*0.99)) continue;
- if (mass>(la_tab[i].mass*1.01)) continue;
- PDG=la_tab[i].pdg;
- strncpy(name,la_tab[i].name,10);
- */
- if (Z != fLa_tab[i]->Z) continue;
- if (strange != fLa_tab[i]->strange) continue;
- if (lepton != fLa_tab[i]->lepton) continue;
- if (A != fLa_tab[i]->A) continue;
- if (mass < (fLa_tab[i]->mass*0.99)) continue;
- if (mass > (fLa_tab[i]->mass*1.01)) continue;
- PDG = fLa_tab[i]->pdg;
- strncpy(name,fLa_tab[i]->name,10);
- result=1;
- break;
- }
- }
- }
- else { // must be ion
- if (A>1) {
- PDG = CreatePdgCode(Z,A,strange,0 ); // 1); // MG
- sprintf(name,"Ion_%d_%d", A,Z);
- result=1;
- }
- }
- return (result);
- }
-
- // ------------------------------------------------------------------------
- Bool_t MpdLAQGSMGenerator::CreateNucleus (Int_t Z, Float_t mass, Int_t pdgCode, char pdgName[11])
- {
- // const Double_t kAu2Gev=0.9314943228; // constants from TGeant3::AddParticlesToPdgDataBase code
- const Double_t khSlash = 1.0545726663e-27; // constants from TGeant3::AddParticlesToPdgDataBase code
- const Double_t kErg2Gev = 1/1.6021773349e-3;
- const Double_t khShGev = khSlash*kErg2Gev;
- const Double_t kYear2Sec = 3600*24*365.25;
- Double_t width = khShGev/(12.33*kYear2Sec); // ? EL
- cout << "-I- MpdLAQGSMGenerator::CreateNucleus " << pdgName << " " << pdgCode << endl;
- TParticlePDG* p=0;
- if (Z)
- p = fPDG->AddParticle(pdgName,pdgName,mass,kFALSE,
- width,Z*3,"nucleus",pdgCode,0,kPTHadron);
- else
- p = fPDG->AddParticle(pdgName,pdgName,mass,kFALSE,
- 6.58211889e-25/2.6e-10,0,
- // "special",pdgCode,0,kPTHadron);
- "special",pdgCode,0,kPTNeutron); //AZ!!!
- if (p)
- p->Print();
- return (p);
- }
- // ------------------------------------------------------------------------
- ClassImp(MpdLAQGSMGenerator)
|