ROOT logo
//____________________________________________________________________
//
// $Id$
//
// Script to get the various cross sections, energy loss, and ranges
// of a particular particle type in a particular medium. 
//
// This script should be compiled to speed it up.  
// 
// It creates a tree on the current output file, with the relevant
// information. 
//
// Note, that VMC _must_ be the TGeant3TGeo VMC. 
//
#ifdef __CINT__
XSection()
{
  gROOT->ProcessLine(".x Compile.C(\"$(ALICE_ROOT)/FMD/scripts/XSection.C\"");
  gAlice->InitMC("$(ALICE_ROOT)/FMD/Config.C");
  TFile* file = TFile::Open("xsec.root", "RECREATE");
  GetXsection("FMD_Si$", "pi+");
  file->Close();
}
#else  
#include <TArrayF.h>
#include <TTree.h>
#include <TMath.h>
#include <iostream>
#include <TVirtualMC.h>
#include <TDatabasePDG.h>
#include <TString.h>
#include <TGeoManager.h>
#include <TGeoMedium.h>
#include <TGeant3.h>
#include <TGraph.h>
#include <TAxis.h>
#include <cmath>
#include <string>
#include <vector>
#include <sstream>
#include <iomanip>
/** @defgroup xsec_script X-section script 
    @ingroup FMD_script 
*/
//____________________________________________________________________
/** @ingroup xsec_script
 */
struct Mech 
{
  char* fName;
  char* fTitle;
  char* fUnit;
};

Mech fgMechs[] = 
  {{ "HADF","total hadronic x-section according to FLUKA",      "cm^{1}"},
   { "INEF","hadronic inelastic x-section according to FLUKA",  "cm^{1}"},
   { "ELAF","hadronic elastic x-section according to FLUKA",    "cm^{1}"},
   { "HADG","total hadronic x-section according to GHEISHA",    "cm^{1}"},
   { "INEG","hadronic inelastic x-section according to GHEISHA","cm^{1}"},
   { "ELAG","hadronic elastic x-section according to GHEISHA",  "cm^{1}"},
   { "FISG","nuclear fission x-section according to GHEISHA",   "cm^{1}"},
   { "CAPG","neutron capture x-section according to GHEISHA",   "cm^{1}"},
   { "LOSS","stopping power",                                   "cm^{1}"},
   { "PHOT","photoelectric x-section",                          "MeV/cm"},
   { "ANNI","positron annihilation x-section",                  "cm^{1}"},
   { "COMP","Compton effect x-section",                         "cm^{1}"},
   { "BREM","bremsstrahlung x-section",                         "cm^{1}"},
   { "PAIR","photon and muon direct- pair x-section",           "cm^{1}"},
   { "DRAY","delta-rays x-section",                             "cm^{1}"},
   { "PFIS","photo-fission x-section",                          "cm^{1}"},
   { "RAYL","Rayleigh scattering x-section",                    "cm^{1}"},
   { "MUNU","muon-nuclear interaction x-section",               "cm^{1}"},
   { "RANG","range",                                            "cm"},
   { "STEP","maximum step",                                     "cm"},
   { 0, 0, 0,}};

//____________________________________________________________________
/** @ingroup xsec_script
 */
struct MechValue 
{
  MechValue(const Mech& m, size_t n) 
    : fMech(m), fValues(n), fStatus(0)
  {}
  bool Get(TGeant3* mc, std::vector<float>& t, 
	   std::vector<float>& c, int medNo, int pidNo) 
  {
    TGraph g(t.size());
    g.SetName(fMech.fName);
    g.SetTitle(fMech.fTitle);
    g.GetXaxis()->SetTitle("p [GeV]");
    g.GetYaxis()->SetTitle(fMech.fUnit);
    int ixst;
    mc->Gftmat(medNo, pidNo, fMech.fName, t.size(), 
	       &(t[0]), &(fValues[0]), &(c[0]), ixst);
    fStatus = ixst;
    if (!fStatus) return false;
    for (size_t i = 0; i < t.size(); i++) g.SetPoint(i, t[i], fValues[i]);
    g.Write();
    return true;
  }
  const Mech&        fMech;
  std::vector<float> fValues;
  bool               fStatus;
};

//____________________________________________________________________
/** @ingroup xsec_script
 */
struct XSection 
{
  XSection(size_t n=91, float emin=1e-5, float emax=1e4) 
    : fTKine(n), 
      fCuts(5)
  {
    float dp   = 1. / log10(emax/emin);
    float pmin = log10(emin);
    fTKine[0]  = emin;
    for (size_t i = 1; i < fTKine.size(); i++) {
      float el  = pmin + i * dp;
      fTKine[i] = pow(10, el);
    }
    for (float_array::iterator i = fCuts.begin(); i != fCuts.end(); ++i) 
      *i = 1e-4;
    Mech* mech = &(fgMechs[0]);
    size_t i = 0;
    while (mech->fName) {
      fMechs.push_back(new MechValue(*mech, n));
      mech++;
    }
  }
  void Run(const std::string& medName, const std::string& pdgName) 
  {
    TGeant3* mc = static_cast<TGeant3*>(gMC);
    if (!mc) {
      std::cerr << "Couldn't get VMC" << std::endl;
      return;
    }
    TGeoMedium* medium = gGeoManager->GetMedium(medName.c_str());
    if (!medium) {
      std::cerr << "Couldn't find medium " << medName << std::endl;
      return;
    }
    int medNo = medium->GetMaterial()->GetUniqueID();
    TDatabasePDG* pdgDb = TDatabasePDG::Instance();
    TParticlePDG* pdgP  = pdgDb->GetParticle(pdgName.c_str());
    if (!pdgP) {
      std::cerr << "Couldn't find particle " << pdgName << std::endl;
      return;
    }
    int pdgNo = pdgP->PdgCode();
    int pidNo = mc->IdFromPDG(pdgNo);

    std::stringstream vars;
    vars << "T/F";

    size_t nOk   = 0;
    // Loop over defined mechanisms 
    for (mech_array::iterator i = fMechs.begin(); i != fMechs.end(); ++i) {
      if (!(*i)->Get(mc, fTKine, fCuts, medNo, pidNo))continue;
      vars << ":" << (*i)->fMech.fName;
      nOk ++;
    }
    
    std::stringstream tName;
    tName << medName << "_" << pdgName;
    TTree* tree = new TTree(tName.str().c_str(), tName.str().c_str());

    float_array cache(nOk+1);
    tree->Branch("xsec", &(cache[0]), vars.str().c_str());
    for (size_t i = 0; i < fTKine.size(); i++) {
      cache[0] = fTKine[i];
      int k = 0;
      for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) {
	if (!(*j)->fStatus) continue;
	cache[++k] = (*j)->fValues[i];
      }
      tree->Fill();
    }
    tree->Write();
  }
protected: 
  typedef std::vector<float> float_array;
  typedef std::vector<MechValue*> mech_array;
  float_array  fTKine;
  float_array  fCuts;
  mech_array   fMechs;
};

#if 0
//____________________________________________________________________
/** @ingroup xsec_script
    @param medName 
    @param pdgName 
    @param n 
    @param emin 
    @param emax 
*/
void
GetXsection(const char* medName, const char* pdgName,
	    Int_t n=91, Float_t emin=1e-5, Float_t emax=1e4)
{
  TArrayF tkine(n);
  Float_t dp   = 1/TMath::Log10(emax/emin);
  Float_t pmin = TMath::Log10(emin);
  tkine[0]     = emin;
  for (Int_t i=1; i < tkine.fN; i++) {
    Float_t el = pmin + i * dp;
    tkine[i]   = TMath::Power(10, el);
  }
  TArrayF cuts(5);
  cuts.Reset(1e-4);

  Mech mechs[] = 
    {{ "HADF","total hadronic x-section according to FLUKA","cm^{1}",n,0},
     { "INEF","hadronic inelastic x-section according to FLUKA","cm^{1}",n,0},
     { "ELAF","hadronic elastic x-section according to FLUKA","cm^{1}",n,0},
     { "HADG","total hadronic x-section according to GHEISHA","cm^{1}",n,0},
     { "INEG","hadronic inelastic x-section according to GHEISHA","cm^{1}",n,0},
     { "ELAG","hadronic elastic x-section according to GHEISHA","cm^{1}",n,0},
     { "FISG","nuclear fission x-section according to GHEISHA","cm^{1}",n,0},
     { "CAPG","neutron capture x-section according to GHEISHA","cm^{1}",n,0},
     { "LOSS","stopping power","cm^{1}",n,0},
     { "PHOT","photoelectric x-section","MeV/cm",n,0},
     { "ANNI","positron annihilation x-section","cm^{1}",n,0},
     { "COMP","Compton effect x-section","cm^{1}",n,0},
     { "BREM","bremsstrahlung x-section","cm^{1}",n,0},
     { "PAIR","photon and muon direct- pair x-section","cm^{1}",n,0},
     { "DRAY","delta-rays x-section","cm^{1}",n,0},
     { "PFIS","photo-fission x-section","cm^{1}",n,0},
     { "RAYL","Rayleigh scattering x-section","cm^{1}",n,0},
     { "MUNU","muon-nuclear interaction x-section","cm^{1}",n,0},
     { "RANG","range","cm",n,0},
     { "STEP","maximum step","cm",n,0},
     { 0, 0, 0, 0, 0}};
  TGeant3* mc = (TGeant3*)gMC;
  if (!mc) {
    std::cerr << "Couldn't get VMC" << std::endl;
    return;
  }
  TGeoMedium* medium = gGeoManager->GetMedium(medName);
  if (!medium) {
    std::cerr << "Couldn't find medium " << medName << std::endl;
    return;
  }
  Int_t medNo = medium->GetMaterial()->GetUniqueID();
  TDatabasePDG* pdgDb = TDatabasePDG::Instance();
  TParticlePDG* pdgP  = pdgDb->GetParticle(pdgName);
  if (!pdgP) {
    std::cerr << "Couldn't find particle " << pdgName << std::endl;
    return;
  }
  Int_t pdgNo = pdgP->PdgCode();
  Int_t pidNo = mc->IdFromPDG(pdgNo);
    
  Mech* mech = &(mechs[0]);
  Int_t nMech = 0;
  Int_t nOk = 0;
  TString vars("T/F");
  while (mech->name) {
    cout << mech->name << ": " << mech->title << " ... " << std::flush;
    nMech++;
    Int_t ixst;
    mc->Gftmat(medNo, pidNo, mech->name, n, 
	       tkine.fArray, mech->values.fArray, cuts.fArray, ixst);
    mech->status = ixst;
    if (ixst) {
      nOk++;
      vars.Append(Form(":%s", mech->name));
      if (!strcmp("LOSS", mech->name)) {
	for (Int_t i = 0; i < n; i++) 
	  std::cout << i << "\t" << tkine[i] << "\t" 
		    << mech->values[i] << std::endl;
      }
    }
    std::cout << (ixst ? "ok" : "failed") << std::endl;
    mech++;
  }
  // TFile* file = TFile::Open(Form("xsec-%d.root", pdgNo),
  // "RECREATE");
  TArrayF cache(nOk+1);
  TTree* tree = new TTree(Form("%s_%s", medName, pdgName), 
			  Form("%s_%s", medName, pdgName));
  tree->Branch("xsec", cache.fArray, vars.Data());
  for (Int_t i = 0; i < n; i++) {
    cache[0] = tkine[i];
    Int_t k = 0;
    for (Int_t j = 0; j < nMech; j++) {
      if (mechs[j].status) {
	if (!strcmp(mechs[j].name, "LOSS")) 
	  std::cout << tkine[i] << "\t" << mechs[j].values[i] << std::endl;
	cache[k+1] = mechs[j].values[i];
	k++;
      }
    }
    std::cout << k << "\t" << (k == nOk) << std::endl;
    tree->Fill();
  }
  tree->Write();
}
#endif
#endif
//____________________________________________________________________
//
// EOF
//
 GetXsection.C:1
 GetXsection.C:2
 GetXsection.C:3
 GetXsection.C:4
 GetXsection.C:5
 GetXsection.C:6
 GetXsection.C:7
 GetXsection.C:8
 GetXsection.C:9
 GetXsection.C:10
 GetXsection.C:11
 GetXsection.C:12
 GetXsection.C:13
 GetXsection.C:14
 GetXsection.C:15
 GetXsection.C:16
 GetXsection.C:17
 GetXsection.C:18
 GetXsection.C:19
 GetXsection.C:20
 GetXsection.C:21
 GetXsection.C:22
 GetXsection.C:23
 GetXsection.C:24
 GetXsection.C:25
 GetXsection.C:26
 GetXsection.C:27
 GetXsection.C:28
 GetXsection.C:29
 GetXsection.C:30
 GetXsection.C:31
 GetXsection.C:32
 GetXsection.C:33
 GetXsection.C:34
 GetXsection.C:35
 GetXsection.C:36
 GetXsection.C:37
 GetXsection.C:38
 GetXsection.C:39
 GetXsection.C:40
 GetXsection.C:41
 GetXsection.C:42
 GetXsection.C:43
 GetXsection.C:44
 GetXsection.C:45
 GetXsection.C:46
 GetXsection.C:47
 GetXsection.C:48
 GetXsection.C:49
 GetXsection.C:50
 GetXsection.C:51
 GetXsection.C:52
 GetXsection.C:53
 GetXsection.C:54
 GetXsection.C:55
 GetXsection.C:56
 GetXsection.C:57
 GetXsection.C:58
 GetXsection.C:59
 GetXsection.C:60
 GetXsection.C:61
 GetXsection.C:62
 GetXsection.C:63
 GetXsection.C:64
 GetXsection.C:65
 GetXsection.C:66
 GetXsection.C:67
 GetXsection.C:68
 GetXsection.C:69
 GetXsection.C:70
 GetXsection.C:71
 GetXsection.C:72
 GetXsection.C:73
 GetXsection.C:74
 GetXsection.C:75
 GetXsection.C:76
 GetXsection.C:77
 GetXsection.C:78
 GetXsection.C:79
 GetXsection.C:80
 GetXsection.C:81
 GetXsection.C:82
 GetXsection.C:83
 GetXsection.C:84
 GetXsection.C:85
 GetXsection.C:86
 GetXsection.C:87
 GetXsection.C:88
 GetXsection.C:89
 GetXsection.C:90
 GetXsection.C:91
 GetXsection.C:92
 GetXsection.C:93
 GetXsection.C:94
 GetXsection.C:95
 GetXsection.C:96
 GetXsection.C:97
 GetXsection.C:98
 GetXsection.C:99
 GetXsection.C:100
 GetXsection.C:101
 GetXsection.C:102
 GetXsection.C:103
 GetXsection.C:104
 GetXsection.C:105
 GetXsection.C:106
 GetXsection.C:107
 GetXsection.C:108
 GetXsection.C:109
 GetXsection.C:110
 GetXsection.C:111
 GetXsection.C:112
 GetXsection.C:113
 GetXsection.C:114
 GetXsection.C:115
 GetXsection.C:116
 GetXsection.C:117
 GetXsection.C:118
 GetXsection.C:119
 GetXsection.C:120
 GetXsection.C:121
 GetXsection.C:122
 GetXsection.C:123
 GetXsection.C:124
 GetXsection.C:125
 GetXsection.C:126
 GetXsection.C:127
 GetXsection.C:128
 GetXsection.C:129
 GetXsection.C:130
 GetXsection.C:131
 GetXsection.C:132
 GetXsection.C:133
 GetXsection.C:134
 GetXsection.C:135
 GetXsection.C:136
 GetXsection.C:137
 GetXsection.C:138
 GetXsection.C:139
 GetXsection.C:140
 GetXsection.C:141
 GetXsection.C:142
 GetXsection.C:143
 GetXsection.C:144
 GetXsection.C:145
 GetXsection.C:146
 GetXsection.C:147
 GetXsection.C:148
 GetXsection.C:149
 GetXsection.C:150
 GetXsection.C:151
 GetXsection.C:152
 GetXsection.C:153
 GetXsection.C:154
 GetXsection.C:155
 GetXsection.C:156
 GetXsection.C:157
 GetXsection.C:158
 GetXsection.C:159
 GetXsection.C:160
 GetXsection.C:161
 GetXsection.C:162
 GetXsection.C:163
 GetXsection.C:164
 GetXsection.C:165
 GetXsection.C:166
 GetXsection.C:167
 GetXsection.C:168
 GetXsection.C:169
 GetXsection.C:170
 GetXsection.C:171
 GetXsection.C:172
 GetXsection.C:173
 GetXsection.C:174
 GetXsection.C:175
 GetXsection.C:176
 GetXsection.C:177
 GetXsection.C:178
 GetXsection.C:179
 GetXsection.C:180
 GetXsection.C:181
 GetXsection.C:182
 GetXsection.C:183
 GetXsection.C:184
 GetXsection.C:185
 GetXsection.C:186
 GetXsection.C:187
 GetXsection.C:188
 GetXsection.C:189
 GetXsection.C:190
 GetXsection.C:191
 GetXsection.C:192
 GetXsection.C:193
 GetXsection.C:194
 GetXsection.C:195
 GetXsection.C:196
 GetXsection.C:197
 GetXsection.C:198
 GetXsection.C:199
 GetXsection.C:200
 GetXsection.C:201
 GetXsection.C:202
 GetXsection.C:203
 GetXsection.C:204
 GetXsection.C:205
 GetXsection.C:206
 GetXsection.C:207
 GetXsection.C:208
 GetXsection.C:209
 GetXsection.C:210
 GetXsection.C:211
 GetXsection.C:212
 GetXsection.C:213
 GetXsection.C:214
 GetXsection.C:215
 GetXsection.C:216
 GetXsection.C:217
 GetXsection.C:218
 GetXsection.C:219
 GetXsection.C:220
 GetXsection.C:221
 GetXsection.C:222
 GetXsection.C:223
 GetXsection.C:224
 GetXsection.C:225
 GetXsection.C:226
 GetXsection.C:227
 GetXsection.C:228
 GetXsection.C:229
 GetXsection.C:230
 GetXsection.C:231
 GetXsection.C:232
 GetXsection.C:233
 GetXsection.C:234
 GetXsection.C:235
 GetXsection.C:236
 GetXsection.C:237
 GetXsection.C:238
 GetXsection.C:239
 GetXsection.C:240
 GetXsection.C:241
 GetXsection.C:242
 GetXsection.C:243
 GetXsection.C:244
 GetXsection.C:245
 GetXsection.C:246
 GetXsection.C:247
 GetXsection.C:248
 GetXsection.C:249
 GetXsection.C:250
 GetXsection.C:251
 GetXsection.C:252
 GetXsection.C:253
 GetXsection.C:254
 GetXsection.C:255
 GetXsection.C:256
 GetXsection.C:257
 GetXsection.C:258
 GetXsection.C:259
 GetXsection.C:260
 GetXsection.C:261
 GetXsection.C:262
 GetXsection.C:263
 GetXsection.C:264
 GetXsection.C:265
 GetXsection.C:266
 GetXsection.C:267
 GetXsection.C:268
 GetXsection.C:269
 GetXsection.C:270
 GetXsection.C:271
 GetXsection.C:272
 GetXsection.C:273
 GetXsection.C:274
 GetXsection.C:275
 GetXsection.C:276
 GetXsection.C:277
 GetXsection.C:278
 GetXsection.C:279
 GetXsection.C:280
 GetXsection.C:281
 GetXsection.C:282
 GetXsection.C:283
 GetXsection.C:284
 GetXsection.C:285
 GetXsection.C:286
 GetXsection.C:287
 GetXsection.C:288
 GetXsection.C:289
 GetXsection.C:290
 GetXsection.C:291
 GetXsection.C:292
 GetXsection.C:293
 GetXsection.C:294
 GetXsection.C:295
 GetXsection.C:296
 GetXsection.C:297
 GetXsection.C:298
 GetXsection.C:299
 GetXsection.C:300
 GetXsection.C:301
 GetXsection.C:302
 GetXsection.C:303
 GetXsection.C:304
 GetXsection.C:305
 GetXsection.C:306
 GetXsection.C:307
 GetXsection.C:308