runMC_zdc.C 5.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191
  1. // Macro for running Fair with Geant3 or Geant4 (M. Al-Turany , D. Bertini)
  2. // Modified 22/06/2005 D.Bertini
  3. // Modified ... roleg (adopted for mpdroot)
  4. // Modified 27-Mar-2008 litvin (field definition fixed to try to store the field parameters)
  5. // Modified 26-Dec-2008 litvin (adopted to current mpdroot changes)
  6. // Modified 30-Jan-2009 litvin (new hostname mpd.jinr.ru)
  7. // Modified 17-Nov-2011 litvin (use macro find_path_to_URQMD_files)
  8. runMC_zdc (const char *datadir="")
  9. {
  10. Int_t nEvents = 20;
  11. TString the_datadir=datadir;
  12. if (the_datadir=="")
  13. the_datadir=".";
  14. TString outFile = the_datadir+ Form("/zdc_%devents.root",nEvents);
  15. TString parFile = outFile ;
  16. TStopwatch timer;
  17. timer.Start();
  18. gDebug=0;
  19. #define URQMD
  20. gROOT->LoadMacro("$VMCWORKDIR/macro/mpd/mpdloadlibs.C");
  21. mpdloadlibs(1,1); // load all libraries
  22. // gROOT->LoadMacro("$VMCWORKDIR/macro/mpd/geometry_v2_option.C");
  23. gROOT->LoadMacro("$VMCWORKDIR/macro/mpd/geometry_v2_minimal.C"); // macro for mini-set
  24. FairRunSim *fRun = new FairRunSim();
  25. // set the MC version used
  26. // ------------------------
  27. // fRun->SetName("TGeant4");
  28. fRun->SetName("TGeant3");
  29. fRun->SetOutputFile(outFile.Data());
  30. // Load mpd detectors libraries with changed ZDC geometry file:
  31. // ------------------------
  32. // geometry_v2_option (fRun, kTRUE, "ZDC","zdc_modules84_layers60_16_4.geo");
  33. // geometry_v2_minimal(fRun); // load minimal set of detector geometries
  34. // // (cave+pipe+magnet+tpc+tof)
  35. geometry_v2_minimal(fRun); // load minimal set of detector geometries
  36. // (cave+pipe+magnet+tpc+zdc)
  37. // Create and Set Event Generator
  38. //-------------------------------
  39. FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
  40. fRun->SetGenerator(primGen);
  41. #ifdef URQMD
  42. // Urqmd Generator
  43. TString hostname = gSystem->HostName();
  44. TString dataFile;
  45. dataFile=find_path_to_URQMD_files();
  46. if ((hostname=="lxmpd-ui.jinr.ru")||(hostname=="lxmpd-ui"))
  47. dataFile += "auau.09gev.mbias.10k.f14";
  48. else
  49. dataFile += "auau.09gev.mbias.98k.ftn14";
  50. FairUrqmdGenerator* urqmdGen = new FairUrqmdGenerator(dataFile);
  51. primGen->AddGenerator(urqmdGen);
  52. #else
  53. #ifdef PART
  54. // ------- Particle Generator
  55. FairParticleGenerator* partGen =
  56. new FairParticleGenerator(211, 10, 1, 0, 3, kTRUE);
  57. primGen->AddGenerator(partGen);
  58. #else
  59. #ifdef ION
  60. // ------- Ion Generator
  61. FairIonGenerator *fIongen= new FairIonGenerator(79, 197,79,1, 0.,0., 25, 0.,0.,-1.);
  62. primGen->AddGenerator(fIongen);
  63. #else
  64. #ifdef BOX
  65. // Box Generator
  66. FairBoxGenerator* boxGen = new
  67. FairBoxGenerator(2212, 10); // 2212 =proton ; 10 = multipl.
  68. // FairBoxGenerator* boxGen = new
  69. // FairBoxGenerator(13, 100); // 13 = muon; 1 = multipl.
  70. boxGen->SetPRange(0.25,2.5); // GeV/c //setPRange vs setPtRange
  71. boxGen->SetPhiRange(0, 360); // Azimuth angle range [degree]
  72. boxGen->SetThetaRange(0, 180); // Polar angle in lab system range [degree]
  73. boxGen->SetXYZ(0., 0., 0.); // mm o cm ??
  74. primGen->AddGenerator(boxGen);
  75. #endif
  76. #endif
  77. #endif
  78. #endif
  79. // Field Map Definition
  80. // --------------------
  81. // // 1- Reading the new field map in the old format
  82. // // Constant Field
  83. // PndConstField *fMagField=new PndConstField();
  84. // fMagField->SetField(0, 0 , 5. ); // values are in kG: 1T = 10kG
  85. // // MinX=-75, MinY=-40,MinZ=-12 ,MaxX=75, MaxY=40 ,MaxZ=124 ); // values are in cm
  86. // fMagField->SetFieldRegion(-205, 205, -205, 205, -261, 261); //cm
  87. // // 2- Reading the new field map in the new format
  88. // // FairField *fMagField= new FairFieldMapSym3("FieldActive");
  89. // // Active Shielding
  90. // fRun->SetField(fMagField);
  91. PndMultiField *fField= new PndMultiField();
  92. PndConstField *fMagField=new PndConstField();
  93. fMagField->SetField(0, 0 , 5. ); // values are in kG: 1T = 10kG
  94. // fMagField->SetFieldRegion(-205, 205, -205, 205, -270, 270); //cm
  95. #define magnet_v4
  96. #ifdef magnet_v4
  97. fMagField->SetFieldRegion(-230, 230, -230, 230, -375, 375); //cm // magnet_v4 or later
  98. #else
  99. fMagField->SetFieldRegion(-205, 205, -205, 205, -270, 270); //cm // old magnet
  100. #endif
  101. fField->AddField(fMagField);
  102. fRun->SetField(fField);
  103. fRun->SetStoreTraj(kTRUE);
  104. fRun->Init();
  105. // -Trajectories Visualization (TGeoManager Only )
  106. // -----------------------------------------------
  107. fRun->SetStoreTraj(kFALSE);
  108. ;
  109. // Set cuts for storing the trajectpries
  110. FairTrajFilter* trajFilter = FairTrajFilter::Instance();
  111. trajFilter->SetStepSizeCut(0.01); // 1 cm
  112. // trajFilter->SetVertexCut(-2000., -2000., 4., 2000., 2000., 100.);
  113. trajFilter->SetMomentumCutP(.50); // p_lab > 500 MeV
  114. // trajFilter->SetEnergyCut(.2, 3.02); // 0 < Etot < 1.04 GeV
  115. trajFilter->SetStorePrimaries(kTRUE);
  116. trajFilter->SetStoreSecondaries(kFALSE);
  117. // trajFilter->SetStoreSecondaries(kTRUE);
  118. // Fill the Parameter containers for this run
  119. //-------------------------------------------
  120. FairRuntimeDb *rtdb=fRun->GetRuntimeDb();
  121. Bool_t kParameterMerged=kTRUE;
  122. FairParRootFileIo* output=new FairParRootFileIo(kParameterMerged);
  123. // output->open(parFile.Data());
  124. output->open(gFile);
  125. rtdb->setOutput(output);
  126. rtdb->saveOutput();
  127. rtdb->print();
  128. // Transport nEvents
  129. // -----------------
  130. fRun->Run(nEvents);
  131. timer.Stop();
  132. Double_t rtime = timer.RealTime();
  133. Double_t ctime = timer.CpuTime();
  134. printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
  135. cout << " Test passed" << endl;
  136. cout << " All ok " << endl;
  137. exit(0);
  138. }