runMC_ndet.C 5.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196
  1. // Macro for running Fair with Geant3 or Geant4 (M. Al-Turany , D. Bertini)
  2. // Modified 22/06/2005 D.Bertini
  3. runMC_ndet (const char *datadir="")
  4. {
  5. TString the_datadir=datadir;
  6. if (the_datadir=="")
  7. the_datadir=".";
  8. TString outFile = the_datadir+ "/evetest.root";
  9. //TString parFile = the_datadir+ "/testparams.root";
  10. TString parFile = outFile;
  11. TStopwatch timer;
  12. timer.Start();
  13. gDebug=0;
  14. #define URQMD
  15. gROOT->LoadMacro("$VMCWORKDIR/macro/mpd/mpdloadlibs.C");
  16. mpdloadlibs(0,1); // load main and detectors libraries
  17. gROOT->LoadMacro("$VMCWORKDIR/macro/mpd/geometry_v2_minimal.C");
  18. FairRunSim *fRun = new FairRunSim();
  19. // set the MC version used
  20. // ------------------------
  21. // fRun->SetName("TGeant4");
  22. fRun->SetName("TGeant3");
  23. // Choose the Geant Navigation System
  24. // fRun->SetGeoModel("G3Native");
  25. fRun->SetOutputFile(outFile.Data());
  26. geometry_v2_minimal(fRun); // load mpd "standard" + "ndet" geometry
  27. FairDetector *NDet= new MpdNDet("Ndet", kTRUE , "$VMCWORKDIR/geometry/ndet.geo");
  28. fRun->AddModule(NDet);
  29. // Create and Set Event Generator
  30. //-------------------------------
  31. FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
  32. fRun->SetGenerator(primGen);
  33. #ifdef URQMD
  34. // Urqmd Generator
  35. TString hostname = gSystem->HostName();
  36. TString dataFile;
  37. // if ((hostname=="nc2")||(hostname=="nc3")||
  38. // (hostname=="nc2.jinr.ru")||(hostname=="nc3.jinr.ru")) { // nc2, nc3
  39. // // dataFile = "/1/data4mpd/";
  40. // dataFile = "/1/data4mpd/UrQMD/1.3/";
  41. // dataFile += "auau.09gev.mbias.98k.ftn14";
  42. // }
  43. // else {
  44. // if ((hostname=="lxmpd-ui.jinr.ru")||(hostname=="lxmpd-ui")) { // linux farm
  45. // dataFile = "/opt/exp_soft/mpd/urqmd/";
  46. // dataFile += "auau.09gev.mbias.10k.f14";
  47. // }
  48. // else {
  49. // if ((hostname=="mpd")||(hostname=="mpd.jinr.ru")||(hostname=="nc12.jinr.ru")||
  50. // (hostname=="nc11.jinr.ru")||(hostname=="nc13.jinr.ru")||(hostname=="nc14.jinr.ru"))
  51. // dataFile = "/opt/data/"; // mpd, nc11
  52. // else{
  53. // //------------------- kanske.itep.ru ----------------// // Moscow
  54. // dataFile ="/scratch2/kmikhail/data4mpd/UrQMD/1.3/"; //
  55. // dataFile += "auau.03.8gev.0_3fm.10k.ftn14"; //
  56. // //---------------------------------------------------//
  57. // }
  58. // if (!((hostname=="kanske")||(hostname=="kanske.itep.ru"))) // Dubna
  59. // dataFile += "auau.09gev.mbias.98k.ftn14";
  60. // }
  61. // }
  62. dataFile=find_path_to_URQMD_files();
  63. if ((hostname=="lxmpd-ui.jinr.ru")||(hostname=="lxmpd-ui"))
  64. dataFile += "auau.09gev.mbias.10k.f14";
  65. else {
  66. if ((hostname=="kanske")||(hostname=="kanske.itep.ru")) // Moscow
  67. dataFile += "auau.03.8gev.0_3fm.10k.ftn14";
  68. else
  69. dataFile += "auau.09gev.mbias.98k.ftn14";
  70. }
  71. FairUrqmdGenerator* urqmdGen = new FairUrqmdGenerator(dataFile);
  72. primGen->AddGenerator(urqmdGen);
  73. #else
  74. #ifdef PART
  75. // ------- Particle Generator
  76. FairParticleGenerator* partGen =
  77. new FairParticleGenerator(211, 10, 1, 0, 3, kTRUE);
  78. primGen->AddGenerator(partGen);
  79. #else
  80. #ifdef ION
  81. // ------- Ion Generator
  82. FairIonGenerator *fIongen= new FairIonGenerator(79, 197,79,1, 0.,0., 25, 0.,0.,-1.);
  83. primGen->AddGenerator(fIongen);
  84. #else
  85. #ifdef BOX
  86. // Box Generator
  87. FairBoxGenerator* boxGen = new
  88. FairBoxGenerator(13, 1); // 13 = muon; 1 = multipl.
  89. boxGen->SetPRange(0.25,2.5); // GeV/c //setPRange vs setPtRange
  90. boxGen->SetPhiRange(0, 360); // Azimuth angle range [degree]
  91. boxGen->SetThetaRange(0, 180); // Polar angle in lab system range [degree]
  92. boxGen->SetXYZ(0., 0., 0.); // mm o cm ??
  93. primGen->AddGenerator(boxGen);
  94. #endif
  95. #endif
  96. #endif
  97. #endif
  98. // Field Map Definition
  99. // --------------------
  100. // 1- Reading the new field map in the old format
  101. // Constant Field
  102. PndConstField *fMagField = new PndConstField();
  103. fMagField->SetField(0, 0 , 5. ); // values are in kG: 1T = 10kG
  104. // MinX=-75, MinY=-40,MinZ=-12 ,MaxX=75, MaxY=40 ,MaxZ=124 ); // values are in cm
  105. fMagField->SetFieldRegion(-205, 205, -205, 205, -261, 261); //cm
  106. fRun->SetField(fMagField);
  107. fRun->SetStoreTraj(kTRUE);
  108. fRun->Init();
  109. // -Trajectories Visualization (TGeoManager Only )
  110. // -----------------------------------------------
  111. //fRun->SetStoreTraj(kFALSE);
  112. ;
  113. // Set cuts for storing the trajectpries
  114. FairTrajFilter* trajFilter = FairTrajFilter::Instance();
  115. trajFilter->SetStepSizeCut(0.01); // 1 cm
  116. // trajFilter->SetVertexCut(-2000., -2000., 4., 2000., 2000., 100.);
  117. trajFilter->SetMomentumCutP(.50); // p_lab > 500 MeV
  118. // trajFilter->SetEnergyCut(.2, 3.02); // 0 < Etot < 1.04 GeV
  119. trajFilter->SetStorePrimaries(kTRUE);
  120. trajFilter->SetStoreSecondaries(kFALSE);
  121. // trajFilter->SetStoreSecondaries(kTRUE);
  122. // Fill the Parameter containers for this run
  123. //-------------------------------------------
  124. FairRuntimeDb *rtdb=fRun->GetRuntimeDb();
  125. Bool_t kParameterMerged=kTRUE;
  126. FairParRootFileIo* output=new FairParRootFileIo(kParameterMerged);
  127. // output->open(parFile.Data());
  128. output->open(gFile);
  129. rtdb->setOutput(output);
  130. rtdb->saveOutput();
  131. rtdb->print();
  132. // Transport nEvents
  133. // -----------------
  134. Int_t nEvents = 1; // 10;//1000;
  135. fRun->Run(nEvents);
  136. timer.Stop();
  137. Double_t rtime = timer.RealTime();
  138. Double_t ctime = timer.CpuTime();
  139. printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
  140. cout << " Test passed" << endl;
  141. cout << " All ok " << endl;
  142. exit(0);
  143. }