[u/mrichter/AliRoot.git] / FMD / scripts / GetXsection.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. 
//
#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
//
Commit	Line	Data
a1f80595	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	//
2b893216	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
a1f80595	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>
2b893216	35	#include <TGraph.h>
	36	#include <TAxis.h>
	37	#include <cmath>
	38	#include <string>
	39	#include <vector>
	40	#include <sstream>
	41	#include <iomanip>
9b48326f	42	/** @defgroup xsec_script X-section script
	43	@ingroup FMD_script
	44	*/
bf000c32	45	//____________________________________________________________________
9b48326f	46	/** @ingroup xsec_script
9b48326f	47	*/
a1f80595	48	struct Mech
a1f80595	49	{
2b893216	50	char* fName;
	51	char* fTitle;
	52	char* fUnit;
a1f80595	53	};
a1f80595	54
2b893216	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	};
a1f80595	190
2b893216	191	#if 0
bf000c32	192	//____________________________________________________________________
9b48326f	193	/** @ingroup xsec_script
	194	@param medName
	195	@param pdgName
	196	@param n
	197	@param emin
	198	@param emax
	199	*/
a1f80595	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	}
2b893216	302	#endif
2b893216	303	#endif
a1f80595	304	//____________________________________________________________________
	305	//
	306	// EOF
	307	//