FMD/scripts/XSection.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 #include <TGraph.h>
  26 #include <TGraphErrors.h>
  27 #include <TAxis.h>
  28 #include <TFile.h>
  29 #include <TH2.h>
  30 #include <TLegend.h>
  31 #include <AliRun.h>
  32 #include <TCanvas.h>
  33 #include <cmath>
  34 #include <string>
  35 #include <vector>
  36 #include <sstream>
  37 #include <iomanip>
  38 /** @defgroup xsec_script X-section script
  39     @ingroup FMD_script
  40 */
  41 //____________________________________________________________________
  42 /** @ingroup xsec_script
  43  */
  44 struct Mech
  45 {
  46   char* fName;
  47   char* fTitle;
  48   char* fUnit;
  49   int   fColor;
  50   int   fStyle;
  51 };
  52
  53 Mech fgMechs[] =
  54   {{ "HADF","FLUKA total hadronic x-section",        "cm^{-1}",1,1},
  55    { "INEF","FLUKA hadronic inelastic x-section",    "cm^{-1}",2,1},
  56    { "ELAF","FLUKA hadronic elastic x-section",      "cm^{-1}",3,1},
  57    { "HADG","GHEISHA total hadronic x-section",      "cm^{-1}",4,1},
  58    { "INEG","GHEISHA hadronic inelastic x-section",  "cm^{-1}",6,1},
  59    { "ELAG","GHEISHA hadronic elastic x-section",    "cm^{-1}",7,1},
  60    { "FISG","GHEISHA nuclear fission x-section",     "cm^{-1}",1,2},
  61    { "CAPG","GHEISHA neutron capture x-section",     "cm^{-1}",2,2},
  62    { "LOSS","stopping power",                        "MeV/cm", 3,2},
  63    { "PHOT","photoelectric x-section",               "cm^{-1}",4,2},
  64    { "ANNI","positron annihilation x-section",       "cm^{-1}",6,2},
  65    { "COMP","Compton effect x-section",              "cm^{-1}",7,2},
  66    { "BREM","bremsstrahlung x-section",              "cm^{-1}",1,3},
  67    { "PAIR","photon and muon direct- pair x-section","cm^{-1}",2,3},
  68    { "DRAY","delta-rays x-section",                  "cm^{-1}",3,3},
  69    { "PFIS","photo-fission x-section",               "cm^{-1}",4,3},
  70    { "RAYL","Rayleigh scattering x-section",         "cm^{-1}",6,3},
  71    { "MUNU","muon-nuclear interaction x-section",    "cm^{-1}",7,3},
  72    { "RANG","range",                                 "cm",     1,4},
  73    { "STEP","maximum step",                          "cm",     2,4},
  74    { 0, 0, 0, 0, 0}};
  75
  76 //____________________________________________________________________
  77 /** @ingroup xsec_script
  78  */
  79 struct MechValue
  80 {
  81   MechValue(const Mech& m, size_t n)
  82     : fMech(m), fValues(n), fStatus(0), fELoss(0)
  83   {
  84     fGraph.SetName(fMech.fName);
  85     fGraph.SetTitle(fMech.fTitle);
  86     fGraph.GetXaxis()->SetTitle("#beta#gamma");
  87     fGraph.GetYaxis()->SetTitle(fMech.fUnit);
  88     fGraph.SetLineColor(fMech.fColor);
  89     fGraph.SetLineStyle(fMech.fStyle);
  90     fGraph.SetLineWidth(2);
  91     std::string name(fMech.fName);
  92     if (name != "LOSS") return;
  93
  94     fELoss = new TGraphErrors(n);
  95     fELoss->SetName("ELOSS");
  96     fELoss->SetTitle(fMech.fTitle);
  97     fELoss->GetXaxis()->SetTitle("#beta#gamma");
  98     fELoss->GetYaxis()->SetTitle(fMech.fUnit);
  99     fELoss->SetLineColor(fMech.fColor);
 100     fELoss->SetLineStyle(fMech.fStyle);
 101     fELoss->SetLineWidth(2);
 102     fELoss->SetFillColor(fMech.fColor+1);
 103     fELoss->SetFillStyle(3001);
 104   }
 105   bool Get(TGeant3* mc, std::vector<float>& t,
 106            std::vector<float>& c, int medNo, int pidNo,
 107            float mass)
 108   {
 109     int ixst;
 110     mc->Gftmat(medNo, pidNo, fMech.fName, t.size(),
 111                &(t[0]), &(fValues[0]), &(c[0]), ixst);
 112     fStatus = ixst;
 113     if (!fStatus) return false;
 114     fGraph.Set(t.size());
 115     for (size_t i = 0; i < t.size(); i++) {
 116       fGraph.SetPoint(i, t[i]/mass, fValues[i]);
 117       if (!fELoss) continue;
 118       fELoss->SetPoint(i, t[i]/mass, fValues[i]);
 119       // ~ 5 sigma
 120       fELoss->SetPointError(i, 0, .5 * fValues[i]);
 121     }
 122     fGraph.Write();
 123     if (fELoss) fELoss->Write();
 124     return true;
 125   }
 126   TGraph& Draw()
 127   {
 128     if (fELoss) fELoss->DrawClone("4 same");
 129     fGraph.DrawClone("l same");
 130     return fGraph;
 131   }
 132   const Mech&        fMech;
 133   std::vector<float> fValues;
 134   bool               fStatus;
 135   TGraph             fGraph;
 136   TGraphErrors*      fELoss;
 137 };
 138
 139 //____________________________________________________________________
 140 /** @ingroup xsec_script
 141  */
 142 struct XSections
 143 {
 144   XSections(size_t n=91, float emin=1e-5, float emax=1e4)
 145     : fTKine(n),
 146       fCuts(5)
 147   {
 148     float dp   = 1. / log10(emax/emin);
 149     float pmin = log10(emin);
 150     fTKine[0]  = emin;
 151     for (size_t i = 1; i < fTKine.size(); i++) {
 152       float el  = pmin + i * dp;
 153       fTKine[i] = pow(10, el);
 154     }
 155     for (float_array::iterator i = fCuts.begin(); i != fCuts.end(); ++i)
 156       *i = 1e-4;
 157     Mech* mech = &(fgMechs[0]);
 158     while (mech->fName) {
 159       fMechs.push_back(new MechValue(*mech, n));
 160       mech++;
 161     }
 162   }
 163   void Run(const std::string& medName, const std::string& pdgName)
 164   {
 165     TGeant3* mc = static_cast<TGeant3*>(gMC);
 166     if (!mc) {
 167       std::cerr << "Couldn't get VMC" << std::endl;
 168       return;
 169     }
 170     TGeoMedium* medium = gGeoManager->GetMedium(medName.c_str());
 171     if (!medium) {
 172       std::cerr << "Couldn't find medium " << medName << std::endl;
 173       return;
 174     }
 175     int medNo = medium->GetMaterial()->GetUniqueID();
 176     TDatabasePDG* pdgDb = TDatabasePDG::Instance();
 177     fPDG                = pdgDb->GetParticle(pdgName.c_str());
 178     if (!fPDG) {
 179       std::cerr << "Couldn't find particle " << pdgName << std::endl;
 180       return;
 181     }
 182     int pdgNo = fPDG->PdgCode();
 183     int pidNo = mc->IdFromPDG(pdgNo);
 184
 185     std::stringstream vars;
 186     vars << "betagamma/F";
 187
 188     size_t nOk   = 0;
 189     // Loop over defined mechanisms
 190     for (mech_array::iterator i = fMechs.begin(); i != fMechs.end(); ++i) {
 191       if (!(*i)->Get(mc, fTKine, fCuts, medNo, pidNo, fPDG->Mass()))continue;
 192       vars << ":" << (*i)->fMech.fName;
 193       nOk ++;
 194     }
 195
 196     std::stringstream tName;
 197     tName << medName << "_" << pdgName;
 198     TTree* tree = new TTree(tName.str().c_str(), tName.str().c_str());
 199
 200     float_array cache(nOk+1);
 201     tree->Branch("xsec", &(cache[0]), vars.str().c_str());
 202     for (size_t i = 0; i < fTKine.size(); i++) {
 203       cache[0] = fTKine[i] / fPDG->Mass();
 204       int k = 0;
 205       for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) {
 206         if (!(*j)->fStatus) continue;
 207         cache[++k] = (*j)->fValues[i];
 208       }
 209       tree->Fill();
 210     }
 211     tree->Write();
 212   }
 213   void Draw()
 214   {
 215     float min = 100000;
 216     float max = 0;
 217     std::vector<TGraph*> gs;
 218     float_array bg(fTKine.size());
 219     for (size_t i = 0; i < fTKine.size(); i++) bg[i] = fTKine[i]/fPDG->Mass();
 220     for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) {
 221       if (!(*j)->fStatus) continue;
 222       for (size_t i = 0; i < fTKine.size(); i++) {
 223         if ((*j)->fValues[i] == 0) continue;
 224         min = std::min(min, (*j)->fValues[i]);
 225         max = std::max(max, (*j)->fValues[i]);
 226       }
 227     }
 228
 229     TCanvas* c = new TCanvas("c", "C", 700, 700);
 230     c->SetFillColor(0);
 231     c->SetLogy();
 232     c->SetLogx();
 233     c->SetGridy();
 234     c->SetGridx();
 235     float_array y(101);
 236     float ymin = log10(min);
 237     float dy   = (log10(max)+.5 - log10(min)) / y.size();
 238     for (size_t i = 1; i < y.size(); i++)  y[i] = pow(10, ymin + i * dy);
 239     TH2* f = new TH2F("x", "X-sec",bg.size()-1,&(bg[0]),y.size()-1,&(y[0]));
 240     f->SetXTitle("#beta#gamma");
 241     f->SetDirectory(0);
 242     f->SetStats(kFALSE);
 243     f->Draw();
 244     TLegend* l = new TLegend(0.45, 0.125, 0.90, 0.45);
 245     l->SetFillColor(0);
 246     // l->SetFillStyle(0);
 247     for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) {
 248       if (!(*j)->fStatus) continue;
 249       TGraph& g = (*j)->Draw();
 250       l->AddEntry(&g, g.GetTitle(), "l");
 251     }
 252     l->Draw("same");
 253   }
 254 protected:
 255   typedef std::vector<float> float_array;
 256   typedef std::vector<MechValue*> mech_array;
 257   float_array   fTKine;
 258   float_array   fCuts;
 259   mech_array    fMechs;
 260   TParticlePDG* fPDG;
 261 };
 262
 263 bool init = false;
 264
 265 void XSection(const char* pdgName="pi-")
 266 {
 267   if (!init) {
 268     gAlice->InitMC("$(ALICE_ROOT)/FMD/Config.C");
 269     init = true;
 270   }
 271   TFile* file = TFile::Open("xsec.root", "RECREATE");
 272   XSections xs;
 273   xs.Run("FMD_Si$", pdgName);
 274   xs.Draw();
 275   file->Close();
 276 }
 277 //____________________________________________________________________
 278 //
 279 // EOF
 280 //