[u/mrichter/AliRoot.git] / PWGGA / GammaConv / macros / Extract_IntegratedPi0Yield.C

// provided by Gamma Conversion Group, PWG4, Kathrin Koch, kkoch@physi.uni-heidelberg.de

#include <fstream>
#include <Riostream.h>
//#include <PlottingGammaHistograms.h>

extern TRandom *gRandom;
extern TBenchmark *gBenchmark;
extern TSystem *gSystem;

void Extract_IntegratedPi0Yield(const char* cutSelection="", const char *inputRootFile = "AnalysisResults",const char *path = "./",const char* outputDir="./Output/",const char* outputDatFileBase= "RB-data",const char* suffix=""){	

  Bool_t writePlotFiles=kFALSE; // if one wants to write f.ex .gif files from the plots
  if(suffix!=""){
    writePlotFiles=kTRUE;
  }

  gROOT->Reset();	
  gROOT->SetStyle("Plain");
	
  // rebinning
  Int_t RebinMass = 2;
	
  Double_t BGFit_range[2] = {0.17,0.3};   	
	
  Double_t Pi0Mass = 0.135;
  Double_t Pi0Mass_range[2] = {0.12,0.15};
  Double_t Pi0Width = 0.003;
  Double_t Pi0Width_range[2] = {0.001,0.007};
  Double_t Pi0Slope = 0.007;
  Double_t Pi0Slope_range[2] = {0.001,0.016};
	
  Double_t Pi0FitRange[2] = {0.1,0.16};
	
  Bool_t MC = kFALSE;

	
  // --------------------------------- end of self definitions --------------------------------

  TString filename = Form("%s%s.root",path,inputRootFile);	
  TFile f(filename.Data());  

  TCanvas* c_mass = new TCanvas("c_mass","",200,10,600,600);  // gives the page size
  c_mass->SetFillColor(0);

      
  TDirectory *pwg4dir =(TDirectory*)f.Get(Form("PWG4_GammaConversion_%s",cutSelection));
  TList *fHistosGammaConversion = (TList*)pwg4dir->Get(Form("histogramsAliGammaConversion_%s",cutSelection));
  TList *fESDContainer = (TList*)fHistosGammaConversion->FindObject("ESD histograms");
  TList *fMCContainer = (TList*)fHistosGammaConversion->FindObject("MC histograms");
  TList *fBackgroundContainer = (TList*)fHistosGammaConversion->FindObject("Background histograms");
  TH1F *ESD_Mother_Mass= (TH1F*)fESDContainer->FindObject("ESD_Mother_InvMass");
	    
  TH1F *ESD_Background_Mass=(TH1F*)fBackgroundContainer->FindObject("ESD_Background_InvMass");
	    
  TH1F * ESD_NumberOfContributorsVtx=(TH1F*)fESDContainer->FindObject("ESD_NumberOfContributorsVtx");
  ESD_NumberOfContributorsVtx->SetAxisRange(1.,100.);
	    
  Float_t nGoodEvents = ESD_NumberOfContributorsVtx->Integral();

  ESD_Mother_Mass->Sumw2();
  ESD_Background_Mass->Sumw2();

  // for Pi0 reco efficiency
  Double_t nPi0_MC = 0;
  Double_t nPi0_MC_error = 0;
  if(MC){
    TH1D *MC_Pi0_Pt=fMCContainer->FindObject("MC_Pi0_Pt");
    nPi0_MC = MC_Pi0_Pt->GetEntries();
    nPi0_MC_error = sqrt(nPi0_MC);
  }

  // get counts for reco and background within pi0 range
  TAxis *xaxis_reco = ESD_Mother_Mass->GetXaxis();							
  Int_t r_1 = xaxis_reco->FindBin(BGFit_range[0]);
  Int_t r_2 = xaxis_reco->FindBin(BGFit_range[1]);  
  Double_t r =  ESD_Mother_Mass->Integral(r_1,r_2); // Integral(75,125)
  TAxis *xaxis_back =ESD_Background_Mass->GetXaxis();							
  Int_t b_1 = xaxis_back->FindBin(BGFit_range[0]);
  Int_t b_2 = xaxis_back->FindBin(BGFit_range[1]);   
  Double_t b =  ESD_Background_Mass->Integral(b_1,b_2);
  Double_t norm = 1;
  if(b != 0) norm = r/b;
  ESD_Background_Mass->Sumw2();		
  ESD_Background_Mass->Scale(norm);
	    
  ESD_Mother_Mass->Rebin(RebinMass);
  ESD_Background_Mass->Rebin(RebinMass);

  ESD_Mother_Mass->SetTitle(Form("%s - cut: %s",inputRootFile,cutSelection));
  ESD_Mother_Mass->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
  ESD_Mother_Mass->GetXaxis()->SetRangeUser(0.0,0.3);	
  ESD_Mother_Mass->SetLineColor(2);
  ESD_Mother_Mass->Draw("");
	    
  ESD_Background_Mass->GetXaxis()->SetRangeUser(0.0,0.3);	
  ESD_Background_Mass->SetLineColor(4);
  ESD_Background_Mass->Draw("same");
	    
  if(writePlotFiles){
    c_mass->Print(Form("%sSpectra-%s-%s.%s",outputDir,inputRootFile,cutSelection,suffix));
  }
  c_mass->Clear();

  TAxis *xaxis_Reco = ESD_Mother_Mass->GetXaxis();	
  Double_t ilowmassPi0_bin = xaxis_Reco->FindBin(0.1);  
  Double_t ihighmassPi0_bin = xaxis_Reco->FindBin(0.16);	
  Double_t Reco_error;
  Double_t Reco = ESD_Mother_Mass->IntegralAndError(ilowmassPi0_bin,ihighmassPi0_bin,Reco_error);
  //	    cout << "Reco:   " << Reco << "      error:  " << Reco_error << endl;
	    
  TAxis *xaxis_Back = ESD_Background_Mass->GetXaxis();	
  Double_t ilowmassPi0_bin = xaxis_Back->FindBin(0.1);  
  Double_t ihighmassPi0_bin = xaxis_Back->FindBin(0.16);	
  Double_t Back_error;
  Double_t Back = ESD_Background_Mass->IntegralAndError(ilowmassPi0_bin,ihighmassPi0_bin,Back_error);
  //	    cout << "Back:   " << Back << "      error:  " << Back_error << endl;
	    
		
  TH1D* Signal = (TH1D*)ESD_Mother_Mass->Clone();
  Signal->Sumw2();	
  Signal->Add(ESD_Background_Mass,-1.);
  Signal->Draw();
	    
  Double_t Pi0Amplitude = Signal->GetMaximum();
  Double_t Pi0Amplitude_Min = Pi0Amplitude*80/100;
  Double_t Pi0Amplitude_Max = Pi0Amplitude*100/100;

  TF1 *fReco = new TF1("Exp+Gauss","(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))",Pi0FitRange[0],Pi0FitRange[1]); 
	    
  // fit the peak around pi0 mass and draw		
  fReco->SetParameter(0,Pi0Amplitude);
  fReco->SetParameter(1,Pi0Mass);
  fReco->SetParameter(2,Pi0Width);
  fReco->SetParameter(3,Pi0Slope);
	    
  fReco->SetParLimits(0,Pi0Amplitude_Min,Pi0Amplitude_Max);
  fReco->SetParLimits(1,Pi0Mass_range[0],Pi0Mass_range[1]);
  fReco->SetParLimits(2,Pi0Width_range[0],Pi0Width_range[1]);
  fReco->SetParLimits(3,Pi0Slope_range[0],Pi0Slope_range[1]);
	    
  Signal->Fit(fReco,"RME");
  fReco->SetLineColor(3);
  fReco->SetLineWidth(0.7);
  fReco->SetLineStyle(4);		
  fReco->Draw("same");
	    
  Double_t parameter[4];
  fReco->GetParameters(parameter);
	    
  Double_t Mass = parameter[1];
  Double_t Mass_error = fReco->GetParError(1);
	    
  // calculation of FWHM + error
  Int_t Pi0Amplitude_fit = fReco->GetMaximum(Pi0FitRange[0],Pi0FitRange[1]);
  Int_t Pi0Amplitude_bin = fReco->GetXaxis()->FindBin(Pi0Amplitude_fit);
  Double_t MassFWHM1 = fReco->GetX(Pi0Amplitude_fit*0.5,0.1,0.135);
  Double_t MassFWHM2 = fReco->GetX(Pi0Amplitude_fit*0.5,0.135,0.16);		
  Double_t FWHM = MassFWHM2 - MassFWHM1;

  // + error
  TF1 *fReco_EPlus = new TF1("Exp+Gauss","(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))",Pi0FitRange[0],Pi0FitRange[1]); 
  fReco_EPlus->SetParameter(0,parameter[0] + fReco->GetParError(0));
  fReco_EPlus->SetParameter(1,parameter[1] + fReco->GetParError(1));
  fReco_EPlus->SetParameter(2,parameter[2] + fReco->GetParError(2));
  fReco_EPlus->SetParameter(3,parameter[3] + fReco->GetParError(3));

  Int_t Pi0Amplitude_EPlus_fit = fReco_EPlus->GetMaximum(Pi0FitRange[0],Pi0FitRange[1]);
  Int_t Pi0Amplitude_EPlus_bin = fReco_EPlus->GetXaxis()->FindBin(Pi0Amplitude_EPlus_fit);
  Double_t MassFWHM1_EPlus = fReco_EPlus->GetX(Pi0Amplitude_EPlus_fit*0.5,0.1,0.135);
  Double_t MassFWHM2_EPlus = fReco_EPlus->GetX(Pi0Amplitude_EPlus_fit*0.5,0.135,0.16);		
  Double_t FWHM_EPlus = MassFWHM2_EPlus - MassFWHM1_EPlus;
    
  // - error
  TF1 *fReco_EMinus = new TF1("Exp+Gauss","(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))",Pi0FitRange[0],Pi0FitRange[1]); 
  fReco_EMinus->SetParameter(0,parameter[0] - fReco->GetParError(0));
  fReco_EMinus->SetParameter(1,parameter[1] - fReco->GetParError(1));
  fReco_EMinus->SetParameter(2,parameter[2] - fReco->GetParError(2));
  fReco_EMinus->SetParameter(3,parameter[3] - fReco->GetParError(3));
	    
  Int_t Pi0Amplitude_EMinus_fit = fReco_EMinus->GetMaximum(Pi0FitRange[0],Pi0FitRange[1]);
  Int_t Pi0Amplitude_EMinus_bin = fReco_EMinus->GetXaxis()->FindBin(Pi0Amplitude_EMinus_fit);
  Double_t MassFWHM1_EMinus = fReco_EMinus->GetX(Pi0Amplitude_EMinus_fit*0.5,0.1,0.135);
  Double_t MassFWHM2_EMinus = fReco_EMinus->GetX(Pi0Amplitude_EMinus_fit*0.5,0.135,0.16);		
  Double_t FWHM_EMinus = MassFWHM2_EMinus - MassFWHM1_EMinus;		 
	    
  Double_t Error1 = TMath::Abs(FWHM-FWHM_EPlus);
  Double_t Error2 = TMath::Abs(FWHM-FWHM_EMinus);
	    
  Double_t FWHM_error;
  if(Error1>=Error2) FWHM_error = Error1;
  if(Error1 < Error2) FWHM_error = Error2;
	    
  fstream file;
  file.open(Form("%s%s-%s.dat",outputDir,outputDatFileBase,inputRootFile), ios::out|ios::app);
  file << inputRootFile << "	" << cutSelection << "    " << nGoodEvents<< "	" << Reco << "	" << Reco_error  << "	" << Back << "	" << Back_error << "	" << Mass << "	" << Mass_error << "	" << FWHM << "	" << FWHM_error << "    " << nPi0_MC<< "    " << nPi0_MC_error << endl;
  file.close();

  if(writePlotFiles){
    c_mass->Print(Form("%sInvMass-%s-%s.%s",outputDir,inputRootFile,cutSelection,suffix));
  }
  c_mass->Update();

  delete fReco;
  delete fReco_EPlus;
  delete fReco_EMinus;
  delete c_mass;
}
Commit	Line	Data
baf06eb6	1	// provided by Gamma Conversion Group, PWG4, Kathrin Koch, kkoch@physi.uni-heidelberg.de
	2
	3	#include <fstream>
	4	#include <Riostream.h>
	5	//#include <PlottingGammaHistograms.h>
	6
	7	extern TRandom *gRandom;
	8	extern TBenchmark *gBenchmark;
	9	extern TSystem *gSystem;
	10
	11	void Extract_IntegratedPi0Yield(const char* cutSelection="", const char inputRootFile = "AnalysisResults",const char path = "./",const char* outputDir="./Output/",const char* outputDatFileBase= "RB-data",const char* suffix=""){
	12
	13	Bool_t writePlotFiles=kFALSE; // if one wants to write f.ex .gif files from the plots
	14	if(suffix!=""){
	15	writePlotFiles=kTRUE;
	16	}
	17
	18	gROOT->Reset();
	19	gROOT->SetStyle("Plain");
	20
	21	// rebinning
	22	Int_t RebinMass = 2;
	23
	24	Double_t BGFit_range[2] = {0.17,0.3};
	25
	26	Double_t Pi0Mass = 0.135;
	27	Double_t Pi0Mass_range[2] = {0.12,0.15};
	28	Double_t Pi0Width = 0.003;
	29	Double_t Pi0Width_range[2] = {0.001,0.007};
	30	Double_t Pi0Slope = 0.007;
	31	Double_t Pi0Slope_range[2] = {0.001,0.016};
	32
	33	Double_t Pi0FitRange[2] = {0.1,0.16};
	34
	35	Bool_t MC = kFALSE;
	36
	37
	38	// --------------------------------- end of self definitions --------------------------------
	39
	40	TString filename = Form("%s%s.root",path,inputRootFile);
	41	TFile f(filename.Data());
	42
	43	TCanvas* c_mass = new TCanvas("c_mass","",200,10,600,600); // gives the page size
	44	c_mass->SetFillColor(0);
	45
	46
	47	TDirectory pwg4dir =(TDirectory)f.Get(Form("PWG4_GammaConversion_%s",cutSelection));
	48	TList fHistosGammaConversion = (TList)pwg4dir->Get(Form("histogramsAliGammaConversion_%s",cutSelection));
	49	TList fESDContainer = (TList)fHistosGammaConversion->FindObject("ESD histograms");
	50	TList fMCContainer = (TList)fHistosGammaConversion->FindObject("MC histograms");
	51	TList fBackgroundContainer = (TList)fHistosGammaConversion->FindObject("Background histograms");
	52	TH1F ESD_Mother_Mass= (TH1F)fESDContainer->FindObject("ESD_Mother_InvMass");
	53
	54	TH1F ESD_Background_Mass=(TH1F)fBackgroundContainer->FindObject("ESD_Background_InvMass");
	55
	56	TH1F * ESD_NumberOfContributorsVtx=(TH1F*)fESDContainer->FindObject("ESD_NumberOfContributorsVtx");
	57	ESD_NumberOfContributorsVtx->SetAxisRange(1.,100.);
	58
	59	Float_t nGoodEvents = ESD_NumberOfContributorsVtx->Integral();
	60
	61	ESD_Mother_Mass->Sumw2();
	62	ESD_Background_Mass->Sumw2();
	63
	64	// for Pi0 reco efficiency
65	Double_t nPi0_MC = 0;
66	Double_t nPi0_MC_error = 0;
67	if(MC){
68	TH1D *MC_Pi0_Pt=fMCContainer->FindObject("MC_Pi0_Pt");
69	nPi0_MC = MC_Pi0_Pt->GetEntries();
70	nPi0_MC_error = sqrt(nPi0_MC);
71	}
72
73	// get counts for reco and background within pi0 range
74	TAxis *xaxis_reco = ESD_Mother_Mass->GetXaxis();
75	Int_t r_1 = xaxis_reco->FindBin(BGFit_range[0]);
76	Int_t r_2 = xaxis_reco->FindBin(BGFit_range[1]);
77	Double_t r = ESD_Mother_Mass->Integral(r_1,r_2); // Integral(75,125)
78	TAxis *xaxis_back =ESD_Background_Mass->GetXaxis();
79	Int_t b_1 = xaxis_back->FindBin(BGFit_range[0]);
80	Int_t b_2 = xaxis_back->FindBin(BGFit_range[1]);
81	Double_t b = ESD_Background_Mass->Integral(b_1,b_2);
82	Double_t norm = 1;
83	if(b != 0) norm = r/b;
84	ESD_Background_Mass->Sumw2();
85	ESD_Background_Mass->Scale(norm);
86
87	ESD_Mother_Mass->Rebin(RebinMass);
88	ESD_Background_Mass->Rebin(RebinMass);
89
90	ESD_Mother_Mass->SetTitle(Form("%s - cut: %s",inputRootFile,cutSelection));
91	ESD_Mother_Mass->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
92	ESD_Mother_Mass->GetXaxis()->SetRangeUser(0.0,0.3);
93	ESD_Mother_Mass->SetLineColor(2);
94	ESD_Mother_Mass->Draw("");
95
96	ESD_Background_Mass->GetXaxis()->SetRangeUser(0.0,0.3);
97	ESD_Background_Mass->SetLineColor(4);
98	ESD_Background_Mass->Draw("same");
99
100	if(writePlotFiles){
101	c_mass->Print(Form("%sSpectra-%s-%s.%s",outputDir,inputRootFile,cutSelection,suffix));
102	}
103	c_mass->Clear();
104
105	TAxis *xaxis_Reco = ESD_Mother_Mass->GetXaxis();
106	Double_t ilowmassPi0_bin = xaxis_Reco->FindBin(0.1);
107	Double_t ihighmassPi0_bin = xaxis_Reco->FindBin(0.16);
108	Double_t Reco_error;
109	Double_t Reco = ESD_Mother_Mass->IntegralAndError(ilowmassPi0_bin,ihighmassPi0_bin,Reco_error);
110	// cout << "Reco: " << Reco << " error: " << Reco_error << endl;
111
112	TAxis *xaxis_Back = ESD_Background_Mass->GetXaxis();
113	Double_t ilowmassPi0_bin = xaxis_Back->FindBin(0.1);
114	Double_t ihighmassPi0_bin = xaxis_Back->FindBin(0.16);
115	Double_t Back_error;
116	Double_t Back = ESD_Background_Mass->IntegralAndError(ilowmassPi0_bin,ihighmassPi0_bin,Back_error);
117	// cout << "Back: " << Back << " error: " << Back_error << endl;
118
119
120	TH1D* Signal = (TH1D*)ESD_Mother_Mass->Clone();
121	Signal->Sumw2();
122	Signal->Add(ESD_Background_Mass,-1.);
123	Signal->Draw();
124
125	Double_t Pi0Amplitude = Signal->GetMaximum();
126	Double_t Pi0Amplitude_Min = Pi0Amplitude*80/100;
127	Double_t Pi0Amplitude_Max = Pi0Amplitude*100/100;
128
129	TF1 fReco = new TF1("Exp+Gauss","(x<[1])([0](exp(-0.5((x-[1])/[2])^2)+exp((x-[1])/[3])(1-exp(-0.5((x-[1])/[2])^2))))+(x>=[1])([0]exp(-0.5*((x-[1])/[2])^2))",Pi0FitRange[0],Pi0FitRange[1]);
130
131	// fit the peak around pi0 mass and draw
132	fReco->SetParameter(0,Pi0Amplitude);
133	fReco->SetParameter(1,Pi0Mass);
134	fReco->SetParameter(2,Pi0Width);
135	fReco->SetParameter(3,Pi0Slope);
136
137	fReco->SetParLimits(0,Pi0Amplitude_Min,Pi0Amplitude_Max);
138	fReco->SetParLimits(1,Pi0Mass_range[0],Pi0Mass_range[1]);
139	fReco->SetParLimits(2,Pi0Width_range[0],Pi0Width_range[1]);
140	fReco->SetParLimits(3,Pi0Slope_range[0],Pi0Slope_range[1]);
141
142	Signal->Fit(fReco,"RME");
143	fReco->SetLineColor(3);
144	fReco->SetLineWidth(0.7);
145	fReco->SetLineStyle(4);
146	fReco->Draw("same");
147
148	Double_t parameter[4];
149	fReco->GetParameters(parameter);
150
151	Double_t Mass = parameter[1];
152	Double_t Mass_error = fReco->GetParError(1);
153
154	// calculation of FWHM + error
155	Int_t Pi0Amplitude_fit = fReco->GetMaximum(Pi0FitRange[0],Pi0FitRange[1]);
156	Int_t Pi0Amplitude_bin = fReco->GetXaxis()->FindBin(Pi0Amplitude_fit);
157	Double_t MassFWHM1 = fReco->GetX(Pi0Amplitude_fit*0.5,0.1,0.135);
158	Double_t MassFWHM2 = fReco->GetX(Pi0Amplitude_fit*0.5,0.135,0.16);
159	Double_t FWHM = MassFWHM2 - MassFWHM1;
160
161	// + error
162	TF1 fReco_EPlus = new TF1("Exp+Gauss","(x<[1])([0](exp(-0.5((x-[1])/[2])^2)+exp((x-[1])/[3])(1-exp(-0.5((x-[1])/[2])^2))))+(x>=[1])([0]exp(-0.5*((x-[1])/[2])^2))",Pi0FitRange[0],Pi0FitRange[1]);
163	fReco_EPlus->SetParameter(0,parameter[0] + fReco->GetParError(0));
164	fReco_EPlus->SetParameter(1,parameter[1] + fReco->GetParError(1));
165	fReco_EPlus->SetParameter(2,parameter[2] + fReco->GetParError(2));
166	fReco_EPlus->SetParameter(3,parameter[3] + fReco->GetParError(3));
167
168	Int_t Pi0Amplitude_EPlus_fit = fReco_EPlus->GetMaximum(Pi0FitRange[0],Pi0FitRange[1]);
169	Int_t Pi0Amplitude_EPlus_bin = fReco_EPlus->GetXaxis()->FindBin(Pi0Amplitude_EPlus_fit);
170	Double_t MassFWHM1_EPlus = fReco_EPlus->GetX(Pi0Amplitude_EPlus_fit*0.5,0.1,0.135);
171	Double_t MassFWHM2_EPlus = fReco_EPlus->GetX(Pi0Amplitude_EPlus_fit*0.5,0.135,0.16);
172	Double_t FWHM_EPlus = MassFWHM2_EPlus - MassFWHM1_EPlus;
173
174	// - error
175	TF1 fReco_EMinus = new TF1("Exp+Gauss","(x<[1])([0](exp(-0.5((x-[1])/[2])^2)+exp((x-[1])/[3])(1-exp(-0.5((x-[1])/[2])^2))))+(x>=[1])([0]exp(-0.5*((x-[1])/[2])^2))",Pi0FitRange[0],Pi0FitRange[1]);
176	fReco_EMinus->SetParameter(0,parameter[0] - fReco->GetParError(0));
177	fReco_EMinus->SetParameter(1,parameter[1] - fReco->GetParError(1));
178	fReco_EMinus->SetParameter(2,parameter[2] - fReco->GetParError(2));
179	fReco_EMinus->SetParameter(3,parameter[3] - fReco->GetParError(3));
180
181	Int_t Pi0Amplitude_EMinus_fit = fReco_EMinus->GetMaximum(Pi0FitRange[0],Pi0FitRange[1]);
182	Int_t Pi0Amplitude_EMinus_bin = fReco_EMinus->GetXaxis()->FindBin(Pi0Amplitude_EMinus_fit);
183	Double_t MassFWHM1_EMinus = fReco_EMinus->GetX(Pi0Amplitude_EMinus_fit*0.5,0.1,0.135);
184	Double_t MassFWHM2_EMinus = fReco_EMinus->GetX(Pi0Amplitude_EMinus_fit*0.5,0.135,0.16);
185	Double_t FWHM_EMinus = MassFWHM2_EMinus - MassFWHM1_EMinus;
186
187	Double_t Error1 = TMath::Abs(FWHM-FWHM_EPlus);
188	Double_t Error2 = TMath::Abs(FWHM-FWHM_EMinus);
189
190	Double_t FWHM_error;
191	if(Error1>=Error2) FWHM_error = Error1;
192	if(Error1 < Error2) FWHM_error = Error2;
193
194	fstream file;
195	file.open(Form("%s%s-%s.dat",outputDir,outputDatFileBase,inputRootFile), ios::out\|ios::app);
196	file << inputRootFile << " " << cutSelection << " " << nGoodEvents<< " " << Reco << " " << Reco_error << " " << Back << " " << Back_error << " " << Mass << " " << Mass_error << " " << FWHM << " " << FWHM_error << " " << nPi0_MC<< " " << nPi0_MC_error << endl;
197	file.close();
198
199	if(writePlotFiles){
200	c_mass->Print(Form("%sInvMass-%s-%s.%s",outputDir,inputRootFile,cutSelection,suffix));
201	}
202	c_mass->Update();
203
204	delete fReco;
205	delete fReco_EPlus;
206	delete fReco_EMinus;
207	delete c_mass;
208	}