[u/mrichter/AliRoot.git] / FMD / scripts / XSection.C

//____________________________________________________________________
//
// $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. 
//
#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 <TGraphErrors.h>
#include <TAxis.h>
#include <TFile.h>
#include <TH2.h>
#include <TLegend.h>
#include <AliRun.h>
#include <TCanvas.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;
  int   fColor;
  int   fStyle;
};

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

//____________________________________________________________________
/** @ingroup xsec_script
 */
struct MechValue 
{
  MechValue(const Mech& m, size_t n) 
    : fMech(m), fValues(n), fStatus(0), fELoss(0)
  {
    fGraph.SetName(fMech.fName);
    fGraph.SetTitle(fMech.fTitle);
    fGraph.GetXaxis()->SetTitle("#beta#gamma");
    fGraph.GetYaxis()->SetTitle(fMech.fUnit);
    fGraph.SetLineColor(fMech.fColor);
    fGraph.SetLineStyle(fMech.fStyle);
    fGraph.SetLineWidth(2);
    std::string name(fMech.fName);
    if (name != "LOSS") return;
    
    fELoss = new TGraphErrors(n);
    fELoss->SetName("ELOSS");
    fELoss->SetTitle(fMech.fTitle);
    fELoss->GetXaxis()->SetTitle("#beta#gamma");
    fELoss->GetYaxis()->SetTitle(fMech.fUnit);
    fELoss->SetLineColor(fMech.fColor);
    fELoss->SetLineStyle(fMech.fStyle);
    fELoss->SetLineWidth(2);
    fELoss->SetFillColor(fMech.fColor+1);
    fELoss->SetFillStyle(3001);
  }
  bool Get(TGeant3* mc, std::vector<float>& t, 
	   std::vector<float>& c, int medNo, int pidNo, 
	   float mass) 
  {
    int ixst;
    mc->Gftmat(medNo, pidNo, fMech.fName, t.size(), 
	       &(t[0]), &(fValues[0]), &(c[0]), ixst);
    fStatus = ixst;
    if (!fStatus) return false;
    fGraph.Set(t.size());
    for (size_t i = 0; i < t.size(); i++) {
      fGraph.SetPoint(i, t[i]/mass, fValues[i]);
      if (!fELoss) continue;
      fELoss->SetPoint(i, t[i]/mass, fValues[i]);
      // ~ 5 sigma
      fELoss->SetPointError(i, 0, .5 * fValues[i]);
    }
    fGraph.Write();
    if (fELoss) fELoss->Write();
    return true;
  }
  TGraph& Draw()
  {
    if (fELoss) fELoss->DrawClone("4 same");
    fGraph.DrawClone("l same");
    return fGraph;
  }
  const Mech&        fMech;
  std::vector<float> fValues;
  bool               fStatus;
  TGraph             fGraph;
  TGraphErrors*      fELoss;
};

//____________________________________________________________________
/** @ingroup xsec_script
 */
struct XSections 
{
  XSections(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]);
    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();
    fPDG                = pdgDb->GetParticle(pdgName.c_str());
    if (!fPDG) {
      std::cerr << "Couldn't find particle " << pdgName << std::endl;
      return;
    }
    int pdgNo = fPDG->PdgCode();
    int pidNo = mc->IdFromPDG(pdgNo);

    std::stringstream vars;
    vars << "betagamma/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, fPDG->Mass()))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] / fPDG->Mass();
      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();
  }
  void Draw() 
  {
    float min = 100000;
    float max = 0;
    std::vector<TGraph*> gs;
    float_array bg(fTKine.size());
    for (size_t i = 0; i < fTKine.size(); i++) bg[i] = fTKine[i]/fPDG->Mass();
    for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) {
      if (!(*j)->fStatus) continue;
      for (size_t i = 0; i < fTKine.size(); i++) {
	if ((*j)->fValues[i] == 0) continue;
	min = std::min(min, (*j)->fValues[i]);
	max = std::max(max, (*j)->fValues[i]);
      }
    }

    TCanvas* c = new TCanvas("c", "C", 700, 700);
    c->SetFillColor(0);
    c->SetLogy();
    c->SetLogx();
    c->SetGridy();
    c->SetGridx();
    float_array y(101);
    float ymin = log10(min);
    float dy   = (log10(max)+.5 - log10(min)) / y.size();
    for (size_t i = 1; i < y.size(); i++)  y[i] = pow(10, ymin + i * dy);
    TH2* f = new TH2F("x", "X-sec",bg.size()-1,&(bg[0]),y.size()-1,&(y[0]));
    f->SetXTitle("#beta#gamma");
    f->SetDirectory(0);
    f->SetStats(kFALSE);
    f->Draw();
    TLegend* l = new TLegend(0.45, 0.125, 0.90, 0.45);
    l->SetFillColor(0);
    // l->SetFillStyle(0);
    for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) {
      if (!(*j)->fStatus) continue;
      TGraph& g = (*j)->Draw();
      l->AddEntry(&g, g.GetTitle(), "l");
    }
    l->Draw("same");
  }
protected: 
  typedef std::vector<float> float_array;
  typedef std::vector<MechValue*> mech_array;
  float_array   fTKine;
  float_array   fCuts;
  mech_array    fMechs;
  TParticlePDG* fPDG;
};

bool init = false;

void XSection(const char* pdgName="pi-")
{
  if (!init) {
    gAlice->InitMC("$(ALICE_ROOT)/FMD/Config.C");
    init = true;
  }
  TFile* file = TFile::Open("xsec.root", "RECREATE");
  XSections xs;
  xs.Run("FMD_Si$", pdgName);
  xs.Draw();
  file->Close();
}
//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
2b893216	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	//