[u/mrichter/AliRoot.git] / PWGLF / ThermalFits / FitNPartDependence.C

#if !defined (__CINT__) || (defined(__MAKECINT__))
#include <iostream>
#include <map>
#include <fstream>
#include "TString.h"
#include "TObjArray.h"
#include "TObjString.h"
#include "AliParticleYield.h"
#include "TGraphErrors.h"
#include "TF1.h"
#include "TROOT.h"
#include "TVirtualFitter.h"
#include "TMath.h"
#include "TFitResult.h"
#include "TMinuit.h"
#include "TLatex.h"
#include "TDatabasePDG.h"
#include "TPad.h"
#include "TCanvas.h"

#endif


std::map<TString,Float_t> npartPbPb;
std::map<TString,Float_t> npartPbPbErr;

void ReadCentralityFromFile() ;
Double_t FitShiftedGraphAndExtrapolate(TGraphErrors * gr, Int_t shift, TF1 * f1, const char * centr, Color_t color) ;
Double_t ErrorFunction (Double_t *xx, Double_t *p ) ;
Double_t FuncPlusErr (Double_t *xx, Double_t *p) ;

enum {kNoShift, kShiftUp, kShiftDown};


TF1 * fForErrors = 0;
const int npar = 3;
Double_t matrix[npar][npar];
TString centrFile;
TString systemAndEnergy;
Double_t maxy = 0;
//TString label = "#frac{d#it{N}}{d#it{y}} = #it{a} + #it{b} #times (#it{N}_{part})^{#it{c}}";
TString label = "#frac{d#it{N}}{d#it{y}} = #it{a} + #it{b} #times (d#it{N}/d#it{#eta})^{#it{c}}";
Int_t collSystem  = 2;
Float_t energy = 2760;

void FitNPartDependence() {

  // KStar
  //  centrFile = "npart_PbPb.txt";
  //  centrFile = "dndeta_PbPb.txt";
  // maxy = 50;
  // systemAndEnergy = "Pb-Pb #sqrt{#it{s}}_{NN} = 2.76 TeV";
  // const char * centralityToPlot[] = {   "V0M0020" ,  "V0M2040" ,  "V0M4060" ,  "V0M6080" , 0};
  // const char * centrToExtrapolate = "V0M0010";
  // Int_t pdg = 313;
  // TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("./PbPb_2760_Kstar892.txt");  
  // Deuteron
  // gROOT->ProcessLine(".x figTemplate.C(0,0.00001,400,0.2)");  
  // const char * centralityToPlot[] = {   "V0M0010" ,  "V0M1020" ,  "V0M2040" ,  "V0M4060" ,  "V0M6080",0 };
  // const char * centrToExtrapolate = "V0M0010";
  // Int_t pdg = 1000010020;
  // TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt");
  // Deuteron pPb
  centrFile = "dndeta_pPb.txt";
  const char * centralityToPlot[] = {   "V0A0010", "V0A1020", "V0A2040", "V0A4060", "V0A6000" ,0};
  const char * centrToExtrapolate = "V0A0005";
  Int_t pdg = 1000010020;
  TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("pPb_5020_deuteron.txt");
  maxy = 0.01;
  systemAndEnergy = "p-Pb #sqrt{#it{s}}_{NN} = 5.02 TeV";
  energy = 5020;
  collSystem = 1;
  // Helium3
  // gROOT->ProcessLine(".x figTemplate.C(0,0.00001,400,0.2)");  
  // const char * centralityToPlot[] = {   "V0M0020" ,  "V0M2080" ,0};
  // const char * centrToExtrapolate = "V0M0010";
  // Int_t pdg = 1000020030;
  // TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt");


  //  gPad->GetCanvas()->SetTitle(Form("c%s", TDatabasePDG::Instance()->GetParticle(pdg)->GetName()));

  ReadCentralityFromFile();
  TGraphErrors * grSyst = new TGraphErrors; 
  TGraphErrors * grStat = new TGraphErrors; 
  Double_t maxx = 1.1*npartPbPb[centrToExtrapolate];
  // Function
  // TF1 * f1 = new TF1 ("f1", "[0] + [1]*x", 0, maxx);
  // f1->SetParameters(1,1);
  TF1 * f1 = new TF1 ("f1", "[0] + [1]*x^[2]", 0, maxx);
  f1->SetParameters(0, 6.3266e-04, 8.99883e-01);

  
  gROOT->ProcessLine(Form(".x figTemplate.C(0,0.1,%f,%f)", maxx, maxy));
  gPad->GetCanvas()->SetTitle("cKstar");

  //  f1->FixParameter(0, 0);
  Int_t icentr = 0;
  AliParticleYield * part = 0;
  while (centralityToPlot[icentr]) {
    part =  AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centralityToPlot[icentr]);
    if(part) {
      grSyst->SetPoint     (icentr , npartPbPb[centralityToPlot[icentr]]    , part->GetYield());
      grSyst->SetPointError(icentr , npartPbPbErr[centralityToPlot[icentr]] , part->GetSystError());
      grStat->SetPoint     (icentr , npartPbPb[centralityToPlot[icentr]]    , part->GetYield());
      grStat->SetPointError(icentr , npartPbPbErr[centralityToPlot[icentr]] , part->GetStatError());
      
      //      part->Print();
    }
    icentr++;
  }
  grStat->Draw("P");
  grStat->SetMarkerStyle(24);
  grSyst->Draw("PE2");
  grSyst->SetFillStyle(0);
  Double_t yield          = FitShiftedGraphAndExtrapolate(grStat, kNoShift  , f1, centrToExtrapolate, kRed  );

  // We need to cache the covariance matrix and the function here, otherwise things get messed up
  fForErrors = (TF1*) f1->Clone();
  fForErrors->SetLineColor(kRed);
  fForErrors->SetLineStyle(0);
  fForErrors->Draw("same");
  gMinuit->mnemat(&matrix[0][0],npar);
  // The statistical error is computed propagating the uncertainty on the parameters to the function via the covariance matrix
  TF1 * fShiftedPlus   = new TF1("fShiftedPlus"  , FuncPlusErr, 0, maxx,1); fShiftedPlus->SetParameter(0,+1);
  TF1 * fShiftedMinus  = new TF1("fShiftedMinus" , FuncPlusErr, 0, maxx,1); fShiftedMinus->SetParameter(0,-1);
  fShiftedMinus->SetLineStyle(kDotted);
  fShiftedPlus->SetLineStyle (kDotted);
  fShiftedPlus ->Draw("same");
  fShiftedMinus->Draw("same");
  TF1 * fError = new TF1("fError", ErrorFunction, 0,maxx, 0);

  // The uncertainty on the systematics is computed shifting the graph up and down + refitting
  Double_t errorSystPlus  = FitShiftedGraphAndExtrapolate(grSyst, kShiftUp  , f1, centrToExtrapolate, kRed)-yield;
  Double_t errorSystMinus = FitShiftedGraphAndExtrapolate(grSyst, kShiftDown, f1, centrToExtrapolate, kRed)-yield;

  // Double_t errorStatPlus  = FitShiftedGraphAndExtrapolate(grStat, kShiftUp  , f1, centrToExtrapolate, kBlue) -yield;
  // Double_t errorStatMinus = FitShiftedGraphAndExtrapolate(grStat, kShiftDown, f1, centrToExtrapolate, kBlue) -yield;

  Double_t errorStat = fError->Eval(npartPbPb[centrToExtrapolate]);


  std::cout << "Yield from fit: ("<<centrToExtrapolate<< "="<<npartPbPb[centrToExtrapolate]<<")" 
            << yield 
            << "   +" << errorStat << " " 
            << "   +" << errorSystPlus << " " << errorSystMinus//<< std::endl;
            << "   +" << fError->Eval(npartPbPb[centrToExtrapolate])  << std::endl;
  std::cout << yield << "    " 
            << errorStat << "     " 
            << (errorSystPlus-errorSystMinus)/2 << std::endl;// here we need - errorneg because errorneg is negative (!)
  

  TGraphErrors * gExtrap = new TGraphErrors();
  gExtrap->SetMarkerStyle(20);
  gExtrap->SetPoint(0, npartPbPb[centrToExtrapolate], yield);
  gExtrap->SetPointError(0, 0, (errorSystPlus-errorSystMinus)/2);
  gExtrap->Draw("P");

  TLatex * text = new TLatex (0.2,0.81,systemAndEnergy);
  text->SetNDC();
  text->Draw();
  TLatex * text2 = new TLatex (0.2,0.72,part->GetLatexName());
  text2->SetNDC();
  text2->Draw(); 

  TLatex * func = new TLatex(0.53,0.23,label);
  func->SetNDC();
  func->Draw();


}

Double_t ErrorFunction (Double_t *xx, Double_t *p ) {


  //  Double_t x = xx[0];
  Double_t func[npar];
  // func[0] = 1;
  // func[1] = TMath::Power(x, fForErrors->GetParameter(2));
  // func[2] = fForErrors->GetParameter(1)*TMath::Power(x, fForErrors->GetParameter(2))*TMath::Log(x);
  //  In general, one can compute the derivative numerically using root.
  for(Int_t ipar = 0; ipar < npar; ipar++){
        func[ipar] =  fForErrors->GradientPar(ipar, xx);
  }

  
  Double_t variance = 0;
  
  for(Int_t ipar = 0; ipar < npar; ipar++){
    for(Int_t jpar = 0; jpar < npar; jpar++){
      variance = variance + func[ipar]*func[jpar] * matrix[ipar][jpar];
    }
    
  }
  Double_t error = TMath::Sqrt(variance);
  return error;
}

Double_t FuncPlusErr (Double_t *xx, Double_t *p) {

  Double_t value = fForErrors->Eval(xx[0]);
  Double_t error = ErrorFunction(xx, p);

  if(p[0] > 0) {
    return value + error;
  }
  return value -error;
}


Double_t FitShiftedGraphAndExtrapolate(TGraphErrors * gr, Int_t shift, TF1 * f1, const char * centr, Color_t color) {
  // Shift graph
  // 1 = up
  // 2 = down

  TGraphErrors * grLocal = (TGraphErrors*) gr->Clone();
  Int_t npoint = grLocal->GetN() ;
  if(shift != kNoShift) {
    for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
      Double_t relError = grLocal->GetEY()[ipoint]/grLocal->GetY()[ipoint];
      Double_t value    = grLocal->GetY() [ipoint];
      if(shift == kShiftUp  ) value += grLocal->GetEY()[ipoint];        
      if(shift == kShiftDown) value -= grLocal->GetEY()[ipoint];        
      grLocal->SetPoint     (ipoint, grLocal->GetX() [ipoint], value);
      grLocal->SetPointError(ipoint, grLocal->GetEX()[ipoint], relError*value);
    }
  }
  

  grLocal->Fit(f1, "", "Q");

  Double_t yield = f1->Eval(npartPbPb[centr]);
  TF1 * clone = (TF1*) f1->Clone();
  clone->SetLineWidth(1);
  clone->SetLineStyle(kDashed);
  clone->SetLineColor(color);
  clone->Draw("Same");
  delete grLocal;
  return yield;

}

void ReadCentralityFromFile() {
  ifstream filein (centrFile);
  TString line;
  while (line.ReadLine((filein))) {    
    TObjArray * tokens = line.Tokenize(" ");
    if(tokens->GetEntries() != 3) {
      delete tokens;
      continue;
    }
    TString str        = ((TObjString*)tokens->At(0))->String();
    Float_t npart     = ((TObjString*)tokens->At(1))->String().Atof();
    Float_t npartErr  = ((TObjString*)tokens->At(2))->String().Atof();
    //    std::cout << "["<<str.Data() <<"]" << npart << " " << npartErr << std::endl;
    if(npart) {
      npartPbPb   [str] = npart;
      npartPbPbErr[str]  = npartErr;
    }
    delete tokens;
  }

}
Commit	Line	Data
368b4921	1	#if !defined (__CINT__) \|\| (defined(__MAKECINT__))
	2	#include <iostream>
	3	#include <map>
	4	#include <fstream>
	5	#include "TString.h"
	6	#include "TObjArray.h"
	7	#include "TObjString.h"
	8	#include "AliParticleYield.h"
	9	#include "TGraphErrors.h"
	10	#include "TF1.h"
	11	#include "TROOT.h"
	12	#include "TVirtualFitter.h"
	13	#include "TMath.h"
	14	#include "TFitResult.h"
	15	#include "TMinuit.h"
	16	#include "TLatex.h"
	17	#include "TDatabasePDG.h"
	18	#include "TPad.h"
	19	#include "TCanvas.h"
	20
	21	#endif
	22
	23
	24
	25
	26	std::map<TString,Float_t> npartPbPb;
	27	std::map<TString,Float_t> npartPbPbErr;
	28
	29	void ReadCentralityFromFile() ;
	30	Double_t FitShiftedGraphAndExtrapolate(TGraphErrors * gr, Int_t shift, TF1 * f1, const char * centr, Color_t color) ;
	31	Double_t ErrorFunction (Double_t xx, Double_t p ) ;
	32	Double_t FuncPlusErr (Double_t xx, Double_t p) ;
	33
	34	enum {kNoShift, kShiftUp, kShiftDown};
	35
	36
	37	TF1 * fForErrors = 0;
	38	const int npar = 3;
	39	Double_t matrix[npar][npar];
	40	TString centrFile;
	41	TString systemAndEnergy;
	42	Double_t maxy = 0;
	43	//TString label = "#frac{d#it{N}}{d#it{y}} = #it{a} + #it{b} #times (#it{N}_{part})^{#it{c}}";
	44	TString label = "#frac{d#it{N}}{d#it{y}} = #it{a} + #it{b} #times (d#it{N}/d#it{#eta})^{#it{c}}";
	45	Int_t collSystem = 2;
	46	Float_t energy = 2760;
	47
	48	void FitNPartDependence() {
	49
	50	// KStar
	51	// centrFile = "npart_PbPb.txt";
	52	// centrFile = "dndeta_PbPb.txt";
	53	// maxy = 50;
	54	// systemAndEnergy = "Pb-Pb #sqrt{#it{s}}_{NN} = 2.76 TeV";
	55	// const char * centralityToPlot[] = { "V0M0020" , "V0M2040" , "V0M4060" , "V0M6080" , 0};
	56	// const char * centrToExtrapolate = "V0M0010";
	57	// Int_t pdg = 313;
	58	// TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("./PbPb_2760_Kstar892.txt");
	59	// Deuteron
	60	// gROOT->ProcessLine(".x figTemplate.C(0,0.00001,400,0.2)");
	61	// const char * centralityToPlot[] = { "V0M0010" , "V0M1020" , "V0M2040" , "V0M4060" , "V0M6080",0 };
	62	// const char * centrToExtrapolate = "V0M0010";
	63	// Int_t pdg = 1000010020;
	64	// TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt");
65	// Deuteron pPb
66	centrFile = "dndeta_pPb.txt";
67	const char * centralityToPlot[] = { "V0A0010", "V0A1020", "V0A2040", "V0A4060", "V0A6000" ,0};
68	const char * centrToExtrapolate = "V0A0005";
69	Int_t pdg = 1000010020;
70	TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("pPb_5020_deuteron.txt");
71	maxy = 0.01;
72	systemAndEnergy = "p-Pb #sqrt{#it{s}}_{NN} = 5.02 TeV";
73	energy = 5020;
74	collSystem = 1;
75	// Helium3
76	// gROOT->ProcessLine(".x figTemplate.C(0,0.00001,400,0.2)");
77	// const char * centralityToPlot[] = { "V0M0020" , "V0M2080" ,0};
78	// const char * centrToExtrapolate = "V0M0010";
79	// Int_t pdg = 1000020030;
80	// TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt");
81
82
83
84
85	// gPad->GetCanvas()->SetTitle(Form("c%s", TDatabasePDG::Instance()->GetParticle(pdg)->GetName()));
86
87	ReadCentralityFromFile();
88	TGraphErrors * grSyst = new TGraphErrors;
89	TGraphErrors * grStat = new TGraphErrors;
90	Double_t maxx = 1.1*npartPbPb[centrToExtrapolate];
91	// Function
92	// TF1 * f1 = new TF1 ("f1", "[0] + [1]*x", 0, maxx);
93	// f1->SetParameters(1,1);
94	TF1 * f1 = new TF1 ("f1", "[0] + [1]*x^[2]", 0, maxx);
95	f1->SetParameters(0, 6.3266e-04, 8.99883e-01);
96
97
98	gROOT->ProcessLine(Form(".x figTemplate.C(0,0.1,%f,%f)", maxx, maxy));
99	gPad->GetCanvas()->SetTitle("cKstar");
100
101	// f1->FixParameter(0, 0);
102	Int_t icentr = 0;
103	AliParticleYield * part = 0;
104	while (centralityToPlot[icentr]) {
105	part = AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centralityToPlot[icentr]);
106	if(part) {
107	grSyst->SetPoint (icentr , npartPbPb[centralityToPlot[icentr]] , part->GetYield());
108	grSyst->SetPointError(icentr , npartPbPbErr[centralityToPlot[icentr]] , part->GetSystError());
109	grStat->SetPoint (icentr , npartPbPb[centralityToPlot[icentr]] , part->GetYield());
110	grStat->SetPointError(icentr , npartPbPbErr[centralityToPlot[icentr]] , part->GetStatError());
111
112	// part->Print();
113	}
114	icentr++;
115	}
116	grStat->Draw("P");
117	grStat->SetMarkerStyle(24);
118	grSyst->Draw("PE2");
119	grSyst->SetFillStyle(0);
120	Double_t yield = FitShiftedGraphAndExtrapolate(grStat, kNoShift , f1, centrToExtrapolate, kRed );
121
122	// We need to cache the covariance matrix and the function here, otherwise things get messed up
123	fForErrors = (TF1*) f1->Clone();
124	fForErrors->SetLineColor(kRed);
125	fForErrors->SetLineStyle(0);
126	fForErrors->Draw("same");
127	gMinuit->mnemat(&matrix[0][0],npar);
128	// The statistical error is computed propagating the uncertainty on the parameters to the function via the covariance matrix
129	TF1 * fShiftedPlus = new TF1("fShiftedPlus" , FuncPlusErr, 0, maxx,1); fShiftedPlus->SetParameter(0,+1);
130	TF1 * fShiftedMinus = new TF1("fShiftedMinus" , FuncPlusErr, 0, maxx,1); fShiftedMinus->SetParameter(0,-1);
131	fShiftedMinus->SetLineStyle(kDotted);
132	fShiftedPlus->SetLineStyle (kDotted);
133	fShiftedPlus ->Draw("same");
134	fShiftedMinus->Draw("same");
135	TF1 * fError = new TF1("fError", ErrorFunction, 0,maxx, 0);
136
137	// The uncertainty on the systematics is computed shifting the graph up and down + refitting
138	Double_t errorSystPlus = FitShiftedGraphAndExtrapolate(grSyst, kShiftUp , f1, centrToExtrapolate, kRed)-yield;
139	Double_t errorSystMinus = FitShiftedGraphAndExtrapolate(grSyst, kShiftDown, f1, centrToExtrapolate, kRed)-yield;
140
141	// Double_t errorStatPlus = FitShiftedGraphAndExtrapolate(grStat, kShiftUp , f1, centrToExtrapolate, kBlue) -yield;
142	// Double_t errorStatMinus = FitShiftedGraphAndExtrapolate(grStat, kShiftDown, f1, centrToExtrapolate, kBlue) -yield;
143
144	Double_t errorStat = fError->Eval(npartPbPb[centrToExtrapolate]);
145
146
147	std::cout << "Yield from fit: ("<<centrToExtrapolate<< "="<<npartPbPb[centrToExtrapolate]<<")"
148	<< yield
149	<< " +" << errorStat << " "
150	<< " +" << errorSystPlus << " " << errorSystMinus//<< std::endl;
151	<< " +" << fError->Eval(npartPbPb[centrToExtrapolate]) << std::endl;
152	std::cout << yield << " "
153	<< errorStat << " "
154	<< (errorSystPlus-errorSystMinus)/2 << std::endl;// here we need - errorneg because errorneg is negative (!)
155
156
157	TGraphErrors * gExtrap = new TGraphErrors();
158	gExtrap->SetMarkerStyle(20);
159	gExtrap->SetPoint(0, npartPbPb[centrToExtrapolate], yield);
160	gExtrap->SetPointError(0, 0, (errorSystPlus-errorSystMinus)/2);
161	gExtrap->Draw("P");
162
163	TLatex * text = new TLatex (0.2,0.81,systemAndEnergy);
164	text->SetNDC();
165	text->Draw();
166	TLatex * text2 = new TLatex (0.2,0.72,part->GetLatexName());
167	text2->SetNDC();
168	text2->Draw();
169
170	TLatex * func = new TLatex(0.53,0.23,label);
171	func->SetNDC();
172	func->Draw();
173
174
175	}
176
177	Double_t ErrorFunction (Double_t xx, Double_t p ) {
178
179
180
181	// Double_t x = xx[0];
182	Double_t func[npar];
183	// func[0] = 1;
184	// func[1] = TMath::Power(x, fForErrors->GetParameter(2));
185	// func[2] = fForErrors->GetParameter(1)TMath::Power(x, fForErrors->GetParameter(2))TMath::Log(x);
186	// In general, one can compute the derivative numerically using root.
187	for(Int_t ipar = 0; ipar < npar; ipar++){
188	func[ipar] = fForErrors->GradientPar(ipar, xx);
189	}
190
191
192
193	Double_t variance = 0;
194
195	for(Int_t ipar = 0; ipar < npar; ipar++){
196	for(Int_t jpar = 0; jpar < npar; jpar++){
197	variance = variance + func[ipar]func[jpar] matrix[ipar][jpar];
198	}
199
200	}
201	Double_t error = TMath::Sqrt(variance);
202	return error;
203	}
204
205	Double_t FuncPlusErr (Double_t xx, Double_t p) {
206
207	Double_t value = fForErrors->Eval(xx[0]);
208	Double_t error = ErrorFunction(xx, p);
209
210	if(p[0] > 0) {
211	return value + error;
212	}
213	return value -error;
214	}
215
216
217	Double_t FitShiftedGraphAndExtrapolate(TGraphErrors * gr, Int_t shift, TF1 * f1, const char * centr, Color_t color) {
218	// Shift graph
219	// 1 = up
220	// 2 = down
221
222	TGraphErrors * grLocal = (TGraphErrors*) gr->Clone();
223	Int_t npoint = grLocal->GetN() ;
224	if(shift != kNoShift) {
225	for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
226	Double_t relError = grLocal->GetEY()[ipoint]/grLocal->GetY()[ipoint];
227	Double_t value = grLocal->GetY() [ipoint];
228	if(shift == kShiftUp ) value += grLocal->GetEY()[ipoint];
229	if(shift == kShiftDown) value -= grLocal->GetEY()[ipoint];
230	grLocal->SetPoint (ipoint, grLocal->GetX() [ipoint], value);
231	grLocal->SetPointError(ipoint, grLocal->GetEX()[ipoint], relError*value);
232	}
233	}
234
235
236	grLocal->Fit(f1, "", "Q");
237
238	Double_t yield = f1->Eval(npartPbPb[centr]);
239	TF1 * clone = (TF1*) f1->Clone();
240	clone->SetLineWidth(1);
241	clone->SetLineStyle(kDashed);
242	clone->SetLineColor(color);
243	clone->Draw("Same");
244	delete grLocal;
245	return yield;
246
247	}
248
249	void ReadCentralityFromFile() {
250	ifstream filein (centrFile);
251	TString line;
252	while (line.ReadLine((filein))) {
253	TObjArray * tokens = line.Tokenize(" ");
254	if(tokens->GetEntries() != 3) {
255	delete tokens;
256	continue;
257	}
258	TString str = ((TObjString*)tokens->At(0))->String();
259	Float_t npart = ((TObjString*)tokens->At(1))->String().Atof();
260	Float_t npartErr = ((TObjString*)tokens->At(2))->String().Atof();
261	// std::cout << "["<<str.Data() <<"]" << npart << " " << npartErr << std::endl;
262	if(npart) {
263	npartPbPb [str] = npart;
264	npartPbPbErr[str] = npartErr;
265	}
266	delete tokens;
267	}
268
269	}