FMD/scripts/GetXsection.C

   1 //____________________________________________________________________
   2 //
   3 // $Id$
   4 //
   5 // Script to get the various cross sections, energy loss, and ranges
   6 // of a particular particle type in a particular medium.
   7 //
   8 // This script should be compiled to speed it up.
   9 //
  10 // It creates a tree on the current output file, with the relevant
  11 // information.
  12 //
  13 // Note, that VMC _must_ be the TGeant3TGeo VMC.
  14 //
  15 #include <TArrayF.h>
  16 #include <TTree.h>
  17 #include <TMath.h>
  18 #include <iostream>
  19 #include <TVirtualMC.h>
  20 #include <TDatabasePDG.h>
  21 #include <TString.h>
  22 #include <TGeoManager.h>
  23 #include <TGeoMedium.h>
  24 #include <TGeant3.h>
  25 /** @defgroup xsec_script X-section script
  26     @ingroup FMD_script
  27 */
  28 //____________________________________________________________________
  29 /** @ingroup xsec_script
  30  */
  31 struct Mech
  32 {
  33   char*    name;
  34   char*    title;
  35   char*    unit;
  36   TArrayF  values;
  37   int      status;
  38 };
  39
  40
  41 //____________________________________________________________________
  42 /** @ingroup xsec_script
  43     @param medName
  44     @param pdgName
  45     @param n
  46     @param emin
  47     @param emax
  48 */
  49 void
  50 GetXsection(const char* medName, const char* pdgName,
  51             Int_t n=91, Float_t emin=1e-5, Float_t emax=1e4)
  52 {
  53   TArrayF tkine(n);
  54   Float_t dp   = 1/TMath::Log10(emax/emin);
  55   Float_t pmin = TMath::Log10(emin);
  56   tkine[0]     = emin;
  57   for (Int_t i=1; i < tkine.fN; i++) {
  58     Float_t el = pmin + i * dp;
  59     tkine[i]   = TMath::Power(10, el);
  60   }
  61   TArrayF cuts(5);
  62   cuts.Reset(1e-4);
  63
  64   Mech mechs[] =
  65     {{ "HADF","total hadronic x-section according to FLUKA","cm^{1}",n,0},
  66      { "INEF","hadronic inelastic x-section according to FLUKA","cm^{1}",n,0},
  67      { "ELAF","hadronic elastic x-section according to FLUKA","cm^{1}",n,0},
  68      { "HADG","total hadronic x-section according to GHEISHA","cm^{1}",n,0},
  69      { "INEG","hadronic inelastic x-section according to GHEISHA","cm^{1}",n,0},
  70      { "ELAG","hadronic elastic x-section according to GHEISHA","cm^{1}",n,0},
  71      { "FISG","nuclear fission x-section according to GHEISHA","cm^{1}",n,0},
  72      { "CAPG","neutron capture x-section according to GHEISHA","cm^{1}",n,0},
  73      { "LOSS","stopping power","cm^{1}",n,0},
  74      { "PHOT","photoelectric x-section","MeV/cm",n,0},
  75      { "ANNI","positron annihilation x-section","cm^{1}",n,0},
  76      { "COMP","Compton effect x-section","cm^{1}",n,0},
  77      { "BREM","bremsstrahlung x-section","cm^{1}",n,0},
  78      { "PAIR","photon and muon direct- pair x-section","cm^{1}",n,0},
  79      { "DRAY","delta-rays x-section","cm^{1}",n,0},
  80      { "PFIS","photo-fission x-section","cm^{1}",n,0},
  81      { "RAYL","Rayleigh scattering x-section","cm^{1}",n,0},
  82      { "MUNU","muon-nuclear interaction x-section","cm^{1}",n,0},
  83      { "RANG","range","cm",n,0},
  84      { "STEP","maximum step","cm",n,0},
  85      { 0, 0, 0, 0, 0}};
  86   TGeant3* mc = (TGeant3*)gMC;
  87   if (!mc) {
  88     std::cerr << "Couldn't get VMC" << std::endl;
  89     return;
  90   }
  91   TGeoMedium* medium = gGeoManager->GetMedium(medName);
  92   if (!medium) {
  93     std::cerr << "Couldn't find medium " << medName << std::endl;
  94     return;
  95   }
  96   Int_t medNo = medium->GetMaterial()->GetUniqueID();
  97   TDatabasePDG* pdgDb = TDatabasePDG::Instance();
  98   TParticlePDG* pdgP  = pdgDb->GetParticle(pdgName);
  99   if (!pdgP) {
 100     std::cerr << "Couldn't find particle " << pdgName << std::endl;
 101     return;
 102   }
 103   Int_t pdgNo = pdgP->PdgCode();
 104   Int_t pidNo = mc->IdFromPDG(pdgNo);
 105
 106   Mech* mech = &(mechs[0]);
 107   Int_t nMech = 0;
 108   Int_t nOk = 0;
 109   TString vars("T/F");
 110   while (mech->name) {
 111     cout << mech->name << ": " << mech->title << " ... " << std::flush;
 112     nMech++;
 113     Int_t ixst;
 114     mc->Gftmat(medNo, pidNo, mech->name, n,
 115                tkine.fArray, mech->values.fArray, cuts.fArray, ixst);
 116     mech->status = ixst;
 117     if (ixst) {
 118       nOk++;
 119       vars.Append(Form(":%s", mech->name));
 120       if (!strcmp("LOSS", mech->name)) {
 121         for (Int_t i = 0; i < n; i++)
 122           std::cout << i << "\t" << tkine[i] << "\t"
 123                     << mech->values[i] << std::endl;
 124       }
 125     }
 126     std::cout << (ixst ? "ok" : "failed") << std::endl;
 127     mech++;
 128   }
 129   // TFile* file = TFile::Open(Form("xsec-%d.root", pdgNo),
 130   // "RECREATE");
 131   TArrayF cache(nOk+1);
 132   TTree* tree = new TTree(Form("%s_%s", medName, pdgName),
 133                           Form("%s_%s", medName, pdgName));
 134   tree->Branch("xsec", cache.fArray, vars.Data());
 135   for (Int_t i = 0; i < n; i++) {
 136     cache[0] = tkine[i];
 137     Int_t k = 0;
 138     for (Int_t j = 0; j < nMech; j++) {
 139       if (mechs[j].status) {
 140         if (!strcmp(mechs[j].name, "LOSS"))
 141           std::cout << tkine[i] << "\t" << mechs[j].values[i] << std::endl;
 142         cache[k+1] = mechs[j].values[i];
 143         k++;
 144       }
 145     }
 146     std::cout << k << "\t" << (k == nOk) << std::endl;
 147     tree->Fill();
 148   }
 149   tree->Write();
 150 }
 151 //____________________________________________________________________
 152 //
 153 // EOF
 154 //