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 #ifdef __CINT__
  16 XSection()
  17 {
  18   gROOT->ProcessLine(".x Compile.C(\"$(ALICE_ROOT)/FMD/scripts/XSection.C\"");
  19   gAlice->InitMC("$(ALICE_ROOT)/FMD/Config.C");
  20   TFile* file = TFile::Open("xsec.root", "RECREATE");
  21   GetXsection("FMD_Si$", "pi+");
  22   file->Close();
  23 }
  24 #else
  25 #include <TArrayF.h>
  26 #include <TTree.h>
  27 #include <TMath.h>
  28 #include <iostream>
  29 #include <TVirtualMC.h>
  30 #include <TDatabasePDG.h>
  31 #include <TString.h>
  32 #include <TGeoManager.h>
  33 #include <TGeoMedium.h>
  34 #include <TGeant3.h>
  35 #include <TGraph.h>
  36 #include <TAxis.h>
  37 #include <cmath>
  38 #include <string>
  39 #include <vector>
  40 #include <sstream>
  41 #include <iomanip>
  42 /** @defgroup xsec_script X-section script
  43     @ingroup FMD_script
  44 */
  45 //____________________________________________________________________
  46 /** @ingroup xsec_script
  47  */
  48 struct Mech
  49 {
  50   char* fName;
  51   char* fTitle;
  52   char* fUnit;
  53 };
  54
  55 Mech fgMechs[] =
  56   {{ "HADF","total hadronic x-section according to FLUKA",      "cm^{1}"},
  57    { "INEF","hadronic inelastic x-section according to FLUKA",  "cm^{1}"},
  58    { "ELAF","hadronic elastic x-section according to FLUKA",    "cm^{1}"},
  59    { "HADG","total hadronic x-section according to GHEISHA",    "cm^{1}"},
  60    { "INEG","hadronic inelastic x-section according to GHEISHA","cm^{1}"},
  61    { "ELAG","hadronic elastic x-section according to GHEISHA",  "cm^{1}"},
  62    { "FISG","nuclear fission x-section according to GHEISHA",   "cm^{1}"},
  63    { "CAPG","neutron capture x-section according to GHEISHA",   "cm^{1}"},
  64    { "LOSS","stopping power",                                   "cm^{1}"},
  65    { "PHOT","photoelectric x-section",                          "MeV/cm"},
  66    { "ANNI","positron annihilation x-section",                  "cm^{1}"},
  67    { "COMP","Compton effect x-section",                         "cm^{1}"},
  68    { "BREM","bremsstrahlung x-section",                         "cm^{1}"},
  69    { "PAIR","photon and muon direct- pair x-section",           "cm^{1}"},
  70    { "DRAY","delta-rays x-section",                             "cm^{1}"},
  71    { "PFIS","photo-fission x-section",                          "cm^{1}"},
  72    { "RAYL","Rayleigh scattering x-section",                    "cm^{1}"},
  73    { "MUNU","muon-nuclear interaction x-section",               "cm^{1}"},
  74    { "RANG","range",                                            "cm"},
  75    { "STEP","maximum step",                                     "cm"},
  76    { 0, 0, 0,}};
  77
  78 //____________________________________________________________________
  79 /** @ingroup xsec_script
  80  */
  81 struct MechValue
  82 {
  83   MechValue(const Mech& m, size_t n)
  84     : fMech(m), fValues(n), fStatus(0)
  85   {}
  86   bool Get(TGeant3* mc, std::vector<float>& t,
  87            std::vector<float>& c, int medNo, int pidNo)
  88   {
  89     TGraph g(t.size());
  90     g.SetName(fMech.fName);
  91     g.SetTitle(fMech.fTitle);
  92     g.GetXaxis()->SetTitle("p [GeV]");
  93     g.GetYaxis()->SetTitle(fMech.fUnit);
  94     int ixst;
  95     mc->Gftmat(medNo, pidNo, fMech.fName, t.size(),
  96                &(t[0]), &(fValues[0]), &(c[0]), ixst);
  97     fStatus = ixst;
  98     if (!fStatus) return false;
  99     for (size_t i = 0; i < t.size(); i++) g.SetPoint(i, t[i], fValues[i]);
 100     g.Write();
 101     return true;
 102   }
 103   const Mech&        fMech;
 104   std::vector<float> fValues;
 105   bool               fStatus;
 106 };
 107
 108 //____________________________________________________________________
 109 /** @ingroup xsec_script
 110  */
 111 struct XSection
 112 {
 113   XSection(size_t n=91, float emin=1e-5, float emax=1e4)
 114     : fTKine(n),
 115       fCuts(5)
 116   {
 117     float dp   = 1. / log10(emax/emin);
 118     float pmin = log10(emin);
 119     fTKine[0]  = emin;
 120     for (size_t i = 1; i < fTKine.size(); i++) {
 121       float el  = pmin + i * dp;
 122       fTKine[i] = pow(10, el);
 123     }
 124     for (float_array::iterator i = fCuts.begin(); i != fCuts.end(); ++i)
 125       *i = 1e-4;
 126     Mech* mech = &(fgMechs[0]);
 127     size_t i = 0;
 128     while (mech->fName) {
 129       fMechs.push_back(new MechValue(*mech, n));
 130       mech++;
 131     }
 132   }
 133   void Run(const std::string& medName, const std::string& pdgName)
 134   {
 135     TGeant3* mc = static_cast<TGeant3*>(gMC);
 136     if (!mc) {
 137       std::cerr << "Couldn't get VMC" << std::endl;
 138       return;
 139     }
 140     TGeoMedium* medium = gGeoManager->GetMedium(medName.c_str());
 141     if (!medium) {
 142       std::cerr << "Couldn't find medium " << medName << std::endl;
 143       return;
 144     }
 145     int medNo = medium->GetMaterial()->GetUniqueID();
 146     TDatabasePDG* pdgDb = TDatabasePDG::Instance();
 147     TParticlePDG* pdgP  = pdgDb->GetParticle(pdgName.c_str());
 148     if (!pdgP) {
 149       std::cerr << "Couldn't find particle " << pdgName << std::endl;
 150       return;
 151     }
 152     int pdgNo = pdgP->PdgCode();
 153     int pidNo = mc->IdFromPDG(pdgNo);
 154
 155     std::stringstream vars;
 156     vars << "T/F";
 157
 158     size_t nOk   = 0;
 159     // Loop over defined mechanisms
 160     for (mech_array::iterator i = fMechs.begin(); i != fMechs.end(); ++i) {
 161       if (!(*i)->Get(mc, fTKine, fCuts, medNo, pidNo))continue;
 162       vars << ":" << (*i)->fMech.fName;
 163       nOk ++;
 164     }
 165
 166     std::stringstream tName;
 167     tName << medName << "_" << pdgName;
 168     TTree* tree = new TTree(tName.str().c_str(), tName.str().c_str());
 169
 170     float_array cache(nOk+1);
 171     tree->Branch("xsec", &(cache[0]), vars.str().c_str());
 172     for (size_t i = 0; i < fTKine.size(); i++) {
 173       cache[0] = fTKine[i];
 174       int k = 0;
 175       for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) {
 176         if (!(*j)->fStatus) continue;
 177         cache[++k] = (*j)->fValues[i];
 178       }
 179       tree->Fill();
 180     }
 181     tree->Write();
 182   }
 183 protected:
 184   typedef std::vector<float> float_array;
 185   typedef std::vector<MechValue*> mech_array;
 186   float_array  fTKine;
 187   float_array  fCuts;
 188   mech_array   fMechs;
 189 };
 190
 191 #if 0
 192 //____________________________________________________________________
 193 /** @ingroup xsec_script
 194     @param medName
 195     @param pdgName
 196     @param n
 197     @param emin
 198     @param emax
 199 */
 200 void
 201 GetXsection(const char* medName, const char* pdgName,
 202             Int_t n=91, Float_t emin=1e-5, Float_t emax=1e4)
 203 {
 204   TArrayF tkine(n);
 205   Float_t dp   = 1/TMath::Log10(emax/emin);
 206   Float_t pmin = TMath::Log10(emin);
 207   tkine[0]     = emin;
 208   for (Int_t i=1; i < tkine.fN; i++) {
 209     Float_t el = pmin + i * dp;
 210     tkine[i]   = TMath::Power(10, el);
 211   }
 212   TArrayF cuts(5);
 213   cuts.Reset(1e-4);
 214
 215   Mech mechs[] =
 216     {{ "HADF","total hadronic x-section according to FLUKA","cm^{1}",n,0},
 217      { "INEF","hadronic inelastic x-section according to FLUKA","cm^{1}",n,0},
 218      { "ELAF","hadronic elastic x-section according to FLUKA","cm^{1}",n,0},
 219      { "HADG","total hadronic x-section according to GHEISHA","cm^{1}",n,0},
 220      { "INEG","hadronic inelastic x-section according to GHEISHA","cm^{1}",n,0},
 221      { "ELAG","hadronic elastic x-section according to GHEISHA","cm^{1}",n,0},
 222      { "FISG","nuclear fission x-section according to GHEISHA","cm^{1}",n,0},
 223      { "CAPG","neutron capture x-section according to GHEISHA","cm^{1}",n,0},
 224      { "LOSS","stopping power","cm^{1}",n,0},
 225      { "PHOT","photoelectric x-section","MeV/cm",n,0},
 226      { "ANNI","positron annihilation x-section","cm^{1}",n,0},
 227      { "COMP","Compton effect x-section","cm^{1}",n,0},
 228      { "BREM","bremsstrahlung x-section","cm^{1}",n,0},
 229      { "PAIR","photon and muon direct- pair x-section","cm^{1}",n,0},
 230      { "DRAY","delta-rays x-section","cm^{1}",n,0},
 231      { "PFIS","photo-fission x-section","cm^{1}",n,0},
 232      { "RAYL","Rayleigh scattering x-section","cm^{1}",n,0},
 233      { "MUNU","muon-nuclear interaction x-section","cm^{1}",n,0},
 234      { "RANG","range","cm",n,0},
 235      { "STEP","maximum step","cm",n,0},
 236      { 0, 0, 0, 0, 0}};
 237   TGeant3* mc = (TGeant3*)gMC;
 238   if (!mc) {
 239     std::cerr << "Couldn't get VMC" << std::endl;
 240     return;
 241   }
 242   TGeoMedium* medium = gGeoManager->GetMedium(medName);
 243   if (!medium) {
 244     std::cerr << "Couldn't find medium " << medName << std::endl;
 245     return;
 246   }
 247   Int_t medNo = medium->GetMaterial()->GetUniqueID();
 248   TDatabasePDG* pdgDb = TDatabasePDG::Instance();
 249   TParticlePDG* pdgP  = pdgDb->GetParticle(pdgName);
 250   if (!pdgP) {
 251     std::cerr << "Couldn't find particle " << pdgName << std::endl;
 252     return;
 253   }
 254   Int_t pdgNo = pdgP->PdgCode();
 255   Int_t pidNo = mc->IdFromPDG(pdgNo);
 256
 257   Mech* mech = &(mechs[0]);
 258   Int_t nMech = 0;
 259   Int_t nOk = 0;
 260   TString vars("T/F");
 261   while (mech->name) {
 262     cout << mech->name << ": " << mech->title << " ... " << std::flush;
 263     nMech++;
 264     Int_t ixst;
 265     mc->Gftmat(medNo, pidNo, mech->name, n,
 266                tkine.fArray, mech->values.fArray, cuts.fArray, ixst);
 267     mech->status = ixst;
 268     if (ixst) {
 269       nOk++;
 270       vars.Append(Form(":%s", mech->name));
 271       if (!strcmp("LOSS", mech->name)) {
 272         for (Int_t i = 0; i < n; i++)
 273           std::cout << i << "\t" << tkine[i] << "\t"
 274                     << mech->values[i] << std::endl;
 275       }
 276     }
 277     std::cout << (ixst ? "ok" : "failed") << std::endl;
 278     mech++;
 279   }
 280   // TFile* file = TFile::Open(Form("xsec-%d.root", pdgNo),
 281   // "RECREATE");
 282   TArrayF cache(nOk+1);
 283   TTree* tree = new TTree(Form("%s_%s", medName, pdgName),
 284                           Form("%s_%s", medName, pdgName));
 285   tree->Branch("xsec", cache.fArray, vars.Data());
 286   for (Int_t i = 0; i < n; i++) {
 287     cache[0] = tkine[i];
 288     Int_t k = 0;
 289     for (Int_t j = 0; j < nMech; j++) {
 290       if (mechs[j].status) {
 291         if (!strcmp(mechs[j].name, "LOSS"))
 292           std::cout << tkine[i] << "\t" << mechs[j].values[i] << std::endl;
 293         cache[k+1] = mechs[j].values[i];
 294         k++;
 295       }
 296     }
 297     std::cout << k << "\t" << (k == nOk) << std::endl;
 298     tree->Fill();
 299   }
 300   tree->Write();
 301 }
 302 #endif
 303 #endif
 304 //____________________________________________________________________
 305 //
 306 // EOF
 307 //