changes from gsi. Using mult if no centrality. testfilterbit 128
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / macros / Extract_IntegratedPi0Yield.C
CommitLineData
5fdef45d 1// provided by Gamma Conversion Group, PWGGA, Kathrin Koch, kkoch@physi.uni-heidelberg.de
baf06eb6 2
3#include <fstream>
4#include <Riostream.h>
5//#include <PlottingGammaHistograms.h>
6
7extern TRandom *gRandom;
8extern TBenchmark *gBenchmark;
9extern TSystem *gSystem;
10
11void 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
5fdef45d 47 TDirectory *pwg4dir =(TDirectory*)f.Get(Form("PWGGA_GammaConversion_%s",cutSelection));
baf06eb6 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}