| 47747ea2 |
1 | /* $Id$ */ |
| 2 | |
| 3 | #include "TFile.h" |
| 6aa36480 |
4 | #include "TCanvas.h" |
| 47747ea2 |
5 | #include "TDirectory.h" |
| 6 | #include "TNamed.h" |
| 7 | #include "TF1.h" |
| 8 | #include "TAttLine.h" |
| 9 | #include "TAttMarker.h" |
| 20fb7b55 |
10 | #include "TH1.h" |
| 47747ea2 |
11 | #include "TH1F.h" |
| 12 | #include "TH2F.h" |
| 13 | #include "TSystem.h" |
| 14 | #include "TStyle.h" |
| 47747ea2 |
15 | #include "TMath.h" |
| 16 | #include "TList.h" |
| 17 | #include <TString.h> |
| 18 | #include <TObject.h> |
| 19 | #include <TROOT.h> |
| 20 | #include <THashList.h> |
| 21 | #include <Rtypes.h> |
| 22 | #include <TPRegexp.h> |
| 23 | #include "TFitResult.h" |
| d1f44d99 |
24 | #include <TMap.h> |
| 25 | #include <TObjString.h> |
| 6aa36480 |
26 | #include "TPad.h" |
| 20fb7b55 |
27 | #include "TAxis.h" |
| 9583e896 |
28 | #include <../cint/cint/include/constants.h> |
| 47747ea2 |
29 | |
| 30 | namespace RawProduction { |
| 31 | |
| a9e9e511 |
32 | bool ignoreErrors = false; |
| 04f85d20 |
33 | const int centBinVersion = 2; |
| 47747ea2 |
34 | |
| 35 | // Fit range |
| 36 | Double_t rangeMin=0.05 ; |
| 37 | Double_t rangeMax=0.3 ; |
| a9e9e511 |
38 | |
| 47747ea2 |
39 | // Fit limits |
| 40 | double lowerMass = 0.110; double upperMass = 0.145; |
| 41 | double lowerWidth = 0.001; double upperWidth = 0.012; |
| 42 | |
| 04df14a2 |
43 | // Scale integral range |
| 44 | double scalarEMBSFrom = 0.2, scalarEMBSTo = 0.4; |
| 9583e896 |
45 | double AFrom = 0.125, ATo = 0.16; |
| 46 | double B1From = 0.1, B1To = 0.12; |
| 47 | double B2From = 0.2, B2To = 0.25; |
| 48 | double GetA(TH1* hist, double from=AFrom, double to=ATo); |
| 49 | double GetADiffB(TH1* hist=0, TH1* scaledMixHist=0); |
| 50 | double GetP0(TH1* hist=0); |
| 51 | double GetP1(TH1* hist=0); |
| 04df14a2 |
52 | |
| 47747ea2 |
53 | Double_t PeakPosition(Double_t pt){ |
| 54 | //Fit to the measured peak position |
| 55 | return 4.99292e-003*exp(-pt*9.32300e-001)+1.34944e-001 ; |
| 56 | } |
| 57 | //----------------------------------------------------------------------------- |
| 58 | Double_t PeakWidth(Double_t pt){ |
| 59 | //fit to the measured peak width |
| 60 | Double_t a=0.0068, b=0.0025, c=0.000319 ; |
| 61 | return TMath::Sqrt(a*a+b*b/pt/pt+c*c*pt*pt) ; |
| 62 | } |
| 63 | |
| 499f0073 |
64 | const char* GetCentString(int centrality); |
| 65 | |
| 47747ea2 |
66 | // Pt bin parameters |
| 67 | Int_t nPtBins=0; |
| 68 | Double_t ptBinEdges[1000] = {0}; |
| 69 | int fool; |
| 70 | double GetPtBin(int bin); |
| 71 | void MakePtBins(); |
| a9e9e511 |
72 | |
| 73 | // various parameters |
| 47747ea2 |
74 | const Double_t kMean=0.136 ; //Approximate peak position to facilitate error estimate |
| 75 | const char format[] = ".pdf"; // say if you want .pdf |
| a9e9e511 |
76 | |
| 47747ea2 |
77 | class Input; |
| 78 | |
| 79 | TH1* GetHistogram_cent(Input& input, const TString& name, int centrality); |
| 80 | TH1* MergeHistogram_cent(Input& input, const TString& name, int newCentIndex, int fromCentIndex, int toCentIndex); |
| 81 | int kNCentralityBins = 3; |
| 82 | |
| 83 | Double_t CGausPol1(Double_t * x, Double_t * par); |
| 84 | Double_t CGausPol2(Double_t * x, Double_t * par); |
| 85 | Double_t CGausPol0(Double_t * x, Double_t * par); |
| 86 | Double_t CPol1(Double_t * x, Double_t * par); |
| 87 | Double_t CPol2(Double_t * x, Double_t * par); |
| a9e9e511 |
88 | |
| 89 | |
| d1f44d99 |
90 | // Output Bin |
| a53ea5d3 |
91 | class TriggerBin { |
| d1f44d99 |
92 | public: |
| a53ea5d3 |
93 | TriggerBin(const TString& trigger = "kCentral" ); |
| 94 | TriggerBin(const char* trigger); |
| 3374aa5b |
95 | const TString& Dir() const {return fDir;} |
| d1f44d99 |
96 | const TString& Key() const {return fKey;} |
| 97 | const TString& Trigger() const {return fTrigger;} |
| a9e9e511 |
98 | protected: |
| d1f44d99 |
99 | TString fTrigger; |
| 3374aa5b |
100 | TString fDir; |
| d1f44d99 |
101 | TString fKey; |
| 102 | }; |
| 47747ea2 |
103 | // Object describing the input for the macros |
| 104 | class Input { |
| 105 | public: |
| a53ea5d3 |
106 | Input(const TString& fileName, const TriggerBin& inputBin, TString listPath = ""); |
| 47747ea2 |
107 | TH1 * GetHistogram(const char* name=""); |
| a53ea5d3 |
108 | const TriggerBin& Bin() const {return fBin;} |
| 47747ea2 |
109 | private: |
| 110 | static TFile* fFile; |
| 111 | TList* fList; |
| a53ea5d3 |
112 | TriggerBin fBin; |
| 47747ea2 |
113 | }; |
| a9e9e511 |
114 | |
| d1f44d99 |
115 | // Output Bin |
| a53ea5d3 |
116 | class TriCenPidBin : public TriggerBin { |
| d1f44d99 |
117 | public: |
| a53ea5d3 |
118 | TriCenPidBin(Int_t centrality, const TString& pid, const TString& trigger); |
| d1f44d99 |
119 | Int_t Centrality() const {return fCentrality;} |
| 120 | const TString& PID() const {return fPID;} |
| 121 | private: |
| 122 | Int_t fCentrality; |
| 123 | TString fPID; |
| 124 | }; |
| 125 | // Object describing the output of the macros |
| 126 | class Output { |
| 127 | public: |
| 6b5d1b53 |
128 | Output(const TString& fileName = "RawProduction.root", const char* options = "UPDATE"); |
| 8375919f |
129 | ~Output(); |
| a53ea5d3 |
130 | TH1* GetHistogram(const TString& name, const TriggerBin& inBin); |
| baacba22 |
131 | TH1* GetHistogram(const TString& name); |
| a53ea5d3 |
132 | void SetDir(const TriggerBin& inBin); |
| d1f44d99 |
133 | void Write(); |
| 134 | private: |
| d1f44d99 |
135 | TFile* fFile; |
| d1f44d99 |
136 | }; |
| 47747ea2 |
137 | |
| 5e8d81ec |
138 | void MakePi0Fit(Input& input, const TriggerBin& outBin, Output& output); |
| 139 | void MakePi0FitTCP(Input& input, const TriCenPidBin& outBin, Output& output); |
| 140 | |
| 141 | void MakePi0Fit(Input& input, const TriggerBin& outBin, Output& output) { |
| 142 | MakePtBins(); |
| 143 | Printf("\nMakePi0Fit(%s)", outBin.Key().Data()); |
| 144 | output.SetDir(outBin); |
| 145 | |
| 146 | for(int m1i = 1; m1i <=3; ++m1i) { |
| 147 | for(int m2i = m1i; m2i <=3; ++m2i) { |
| 148 | TStringToken names("A;M;W;p0;p1;p2", ";"); |
| 149 | TStringToken titles("Amplitude;Mass;Width;p0;p1;p2", ";"); |
| 150 | while( names.NextToken() && titles.NextToken() ) { |
| 151 | TString name(Form("%s_M%i%i", names.Data(), m1i, m2i)); |
| 152 | TString title(Form("%s_M%i%i", titles.Data(), m1i, m2i)); |
| 153 | new TH1D(name.Data(), title.Data(), nPtBins,ptBinEdges); |
| 154 | new TH1D(Form("%s_error", name.Data()), title.Data(), nPtBins,ptBinEdges); |
| 155 | } |
| 156 | } |
| 157 | } |
| 158 | |
| 159 | TF1 * pol2Gaus = new TF1("pol2Gaus", CGausPol2, 0., 1., 6); |
| 160 | pol2Gaus->SetParNames("A", "m_{0}", "#sigma", "a_{0}", "a_{1}", "a_{2}"); |
| 161 | for(int m1i = 1; m1i <=3; ++m1i) { |
| 162 | for(int m2i = m1i; m2i <=3; ++m2i) { |
| 163 | TH2F* hPi0MNN = (TH2F*) input.GetHistogram(Form("hPi0M%i%i", m1i, m2i)); |
| 164 | for(Int_t ptBin=1; ptBin<=nPtBins; ptBin++){ |
| 165 | Int_t imin = hPi0MNN->GetYaxis()->FindBin(ptBinEdges[ptBin-1]+0.0001); |
| 166 | Int_t imax = hPi0MNN->GetYaxis()->FindBin(ptBinEdges[ptBin]+0.0001) -1; |
| 167 | TH1D* hPi0MNNProj = hPi0MNN->ProjectionX(Form("pt%03i_hPi0M%i%i",ptBin, m1i, m2i), imin, imax); |
| 168 | hPi0MNNProj->SetTitle(Form("M_{#gamma#gamma}, M%i%i, %.1f<p_{T}<%.1f, %s", m1i, m2i, ptBinEdges[ptBin-1], ptBinEdges[ptBin], outBin.Trigger().Data())); |
| 169 | hPi0MNNProj->Sumw2(); |
| 170 | int entries = hPi0MNNProj->Integral(hPi0MNNProj->FindBin(lowerMass), hPi0MNNProj->FindBin(upperMass)); |
| 171 | Printf("pt bin %i, %3.1f to %3.1f, entries: %i", ptBin, ptBinEdges[ptBin-1], ptBinEdges[ptBin], entries); |
| 172 | if( entries < 10 ) continue; |
| 173 | |
| 174 | |
| 175 | // Error Fix |
| 176 | for(Int_t ib=1; ib<=hPi0MNNProj->GetNbinsX();ib++){if(hPi0MNNProj ->GetBinContent(ib)==0)hPi0MNNProj ->SetBinError(ib,1.);} |
| 177 | |
| 178 | pol2Gaus->SetParameters(entries, 0.134, 0.006, entries, 0, 0); |
| 179 | pol2Gaus->SetParLimits(1, lowerMass, upperMass); |
| 180 | pol2Gaus->SetParLimits(2, lowerWidth, upperWidth); |
| 181 | TFitResultPtr resPtr = hPi0MNNProj->Fit(pol2Gaus, "MSQ", "", rangeMin, rangeMax); |
| 182 | Int_t error = int(resPtr) != 4000; |
| 183 | error = error || TMath::Abs(pol2Gaus->GetParameter(1) - lowerMass) < 0.0001 || TMath::Abs(pol2Gaus->GetParameter(1) - upperMass) <0.0001; |
| 184 | error = error || TMath::Abs(pol2Gaus->GetParameter(2) - lowerWidth) < 0.0001 || TMath::Abs(pol2Gaus->GetParameter(2) - upperWidth) <0.0001; |
| 185 | |
| 186 | if(error) { |
| 187 | ((TH1D*) output.GetHistogram(Form("A_M%i%i_error", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(0)); |
| 188 | ((TH1D*) output.GetHistogram(Form("M_M%i%i_error", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(1)); |
| 189 | ((TH1D*) output.GetHistogram(Form("W_M%i%i_error", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(2)); |
| 190 | ((TH1D*) output.GetHistogram(Form("p0_M%i%i_error", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(3)); |
| 191 | ((TH1D*) output.GetHistogram(Form("p1_M%i%i_error", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(4)); |
| 192 | ((TH1D*) output.GetHistogram(Form("p2_M%i%i_error", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(5)); |
| 193 | ((TH1D*) output.GetHistogram(Form("A_M%i%i_error", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(0)); |
| 194 | ((TH1D*) output.GetHistogram(Form("M_M%i%i_error", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(1)); |
| 195 | ((TH1D*) output.GetHistogram(Form("W_M%i%i_error", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(2)); |
| 196 | ((TH1D*) output.GetHistogram(Form("p0_M%i%i_error", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(3)); |
| 197 | ((TH1D*) output.GetHistogram(Form("p1_M%i%i_error", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(4)); |
| 198 | ((TH1D*) output.GetHistogram(Form("p2_M%i%i_error", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(5)); |
| 199 | } else { |
| 200 | ((TH1D*) output.GetHistogram(Form("A_M%i%i", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(0)); |
| 201 | ((TH1D*) output.GetHistogram(Form("M_M%i%i", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(1)); |
| 202 | ((TH1D*) output.GetHistogram(Form("W_M%i%i", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(2)); |
| 203 | ((TH1D*) output.GetHistogram(Form("p0_M%i%i", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(3)); |
| 204 | ((TH1D*) output.GetHistogram(Form("p1_M%i%i", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(4)); |
| 205 | ((TH1D*) output.GetHistogram(Form("p2_M%i%i", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(5)); |
| 206 | ((TH1D*) output.GetHistogram(Form("A_M%i%i", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(0)); |
| 207 | ((TH1D*) output.GetHistogram(Form("M_M%i%i", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(1)); |
| 208 | ((TH1D*) output.GetHistogram(Form("W_M%i%i", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(2)); |
| 209 | ((TH1D*) output.GetHistogram(Form("p0_M%i%i", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(3)); |
| 210 | ((TH1D*) output.GetHistogram(Form("p1_M%i%i", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(4)); |
| 211 | ((TH1D*) output.GetHistogram(Form("p2_M%i%i", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(5)); |
| 212 | } |
| 213 | } |
| 214 | } |
| 215 | } |
| 216 | } |
| 217 | void MakePi0FitTCP(Input& input, const TriCenPidBin& outBin, Output& output) { |
| 47747ea2 |
218 | MakePtBins(); |
| 20fb7b55 |
219 | Printf("\nMakePi0Fit(%s)", outBin.Key().Data()); |
| 6aa36480 |
220 | output.SetDir(outBin); |
| a9e9e511 |
221 | |
| 47747ea2 |
222 | TH1F * hTotSelEvents = (TH1F*) input.GetHistogram("hTotSelEvents"); |
| 223 | TH2F * hCentrality = (TH2F*) input.GetHistogram("hCenPHOSCells"); |
| 224 | TH1D * hCentralityX = hCentrality->ProjectionX(); |
| 20fb7b55 |
225 | TH2F *hPi0 = (TH2F*)GetHistogram_cent(input, Form("hPi0%s", outBin.PID().Data()), outBin.Centrality()); |
| 226 | TH2F *hPi0Mix = (TH2F*)GetHistogram_cent(input, Form("hMiPi0%s", outBin.PID().Data()), outBin.Centrality()); |
| 47747ea2 |
227 | |
| 228 | printf("TotSelEvents (4): %.0f \n", hTotSelEvents->GetBinContent(4)) ; |
| 229 | printf("Centrality: %.0f \n", hCentralityX->Integral()) ; |
| a9e9e511 |
230 | |
| 47747ea2 |
231 | if( !hPi0 || !hPi0Mix ) { |
| 232 | Printf(Form("no histogram(0x%p, 0x%p), returning", hPi0, hPi0Mix)); |
| 233 | return; |
| 234 | } |
| 235 | |
| 236 | // for temp convas for drawing/monitoring |
| 6aa36480 |
237 | TCanvas* canvas = new TCanvas("cMakePi0Fit", Form("MakePi0Fit Canvas, %s", outBin.Key().Data()),10,10,1200,800); |
| a9e9e511 |
238 | |
| 239 | |
| 47747ea2 |
240 | // Peak Parameterization |
| 241 | // 1. Linear Bg |
| 242 | TF1 * funcRatioFit1 = new TF1("funcRatioFit1",CGausPol1,0.,1.,5) ; |
| 243 | funcRatioFit1->SetParNames("A", "m_{0}", "#sigma", "a_{0}", "a_{1}"); |
| 244 | funcRatioFit1->SetLineWidth(2) ; |
| 245 | funcRatioFit1->SetLineColor(2) ; |
| 246 | // 2. Quadratic Bg |
| 247 | TF1 * funcRatioFit2 = new TF1("funcRatioFit2",CGausPol2,0.,1.,6) ; |
| 248 | funcRatioFit2->SetParNames("A", "m_{0}", "#sigma", "a_{0}", "a_{1}", "a_{2}"); |
| 249 | funcRatioFit2->SetLineWidth(2) ; |
| 250 | funcRatioFit2->SetLineColor(4) ; |
| 251 | funcRatioFit2->SetLineStyle(2) ; |
| 252 | // Other |
| 253 | TF1 * fgs = new TF1("gs",CGausPol0,0.,1.,4) ; |
| 254 | fgs->SetParNames("A", "m_{0}", "#sigma", "B") ; |
| 255 | fgs->SetLineColor(2) ; |
| 256 | fgs->SetLineWidth(1) ; |
| 257 | TF1 * fbg1 = new TF1("bg1",CPol1,0.,1.,2) ; |
| 258 | TF1 * fbg2 = new TF1("bg2",CPol2,0.,1.,3) ; |
| 259 | |
| 260 | // Adding Histograms: |
| 261 | // 1. Linear Bg |
| 262 | TStringToken names("mr1;mr1r;sr1;sr1r;ar1;br1;yr1;yr1int", ";"); |
| 263 | TStringToken titles("Mass;Mass, Ratio Fit;Width;Width, Ratio Fit;a;b;Raw Yield; Raw Yield, integrated", ";"); |
| 6aa36480 |
264 | while( names.NextToken() && titles.NextToken() ) { |
| 265 | new TH1D(names.Data(), titles.Data(), nPtBins,ptBinEdges); |
| 266 | new TH1D(Form("%s_error", names.Data()), titles.Data(), nPtBins,ptBinEdges); |
| 267 | } |
| 47747ea2 |
268 | // 2. Quadratic Bg |
| 269 | TStringToken names2("mr2;mr2r;sr2;sr2r;ar2;br2;cr2;yr2;yr2int", ";"); |
| 270 | TStringToken titles2("Mass;Mass, Ratio Fit;Width;Width, Ratio Fit;a;b;c;Raw Yield; Raw Yield, integrated", ";"); |
| 6aa36480 |
271 | while( names2.NextToken() && titles2.NextToken() ) { |
| 272 | new TH1D(names2.Data(), titles2.Data(), nPtBins,ptBinEdges); |
| 273 | new TH1D(Form("%s_error", names2.Data()), titles2.Data(), nPtBins,ptBinEdges); |
| 274 | } |
| a9e9e511 |
275 | |
| 6aa36480 |
276 | TH1D* hMixRawRatio =new TH1D("hMixRawRatio", "ratio of statistics in RAW and Mixed", nPtBins, ptBinEdges); |
| a9e9e511 |
277 | |
| 47747ea2 |
278 | // Pt slice loop |
| 279 | for(Int_t ptBin=1; ptBin<=nPtBins; ptBin++){ |
| 280 | //gSystem->Sleep(2000); |
| 281 | canvas->Clear(); |
| 282 | canvas->Divide(2,2); |
| 283 | canvas->cd(1); |
| 284 | printf("pt bin %i, %3.1f to %3.1f, ", ptBin, ptBinEdges[ptBin-1], ptBinEdges[ptBin]); |
| a9e9e511 |
285 | |
| 47747ea2 |
286 | //TList* ptBinOutputList = outputList; |
| 287 | TAxis * ptAxis=hPi0->GetYaxis() ; |
| 288 | Int_t imin=ptAxis->FindBin(ptBinEdges[ptBin-1]+0.0001); |
| 289 | Int_t imax=ptAxis->FindBin(ptBinEdges[ptBin]+0.0001) -1; |
| 290 | if( TMath::Abs( ptAxis->GetBinLowEdge(imin) - ptBinEdges[ptBin -1] ) >0.0001 ) { |
| 291 | Printf("\nError: hPi0 lower bin edge (%f) different from ptBinEdges lower (%f)", ptAxis->GetBinLowEdge(imin), ptBinEdges[ptBin -1]); |
| 292 | continue; |
| 293 | } |
| 294 | if( TMath::Abs( ptAxis->GetBinLowEdge(imax+1) - ptBinEdges[ptBin] ) >0.0001 ) { |
| 295 | Printf("\nError: hPi0 upper bin edge (%f) different from ptBinEdges upper (%f)", ptAxis->GetBinLowEdge(imax+1), ptBinEdges[ptBin]); |
| 296 | continue; |
| 297 | } |
| a9e9e511 |
298 | |
| 47747ea2 |
299 | Double_t pt=(ptBinEdges[ptBin]+ptBinEdges[ptBin-1])/2. ; |
| 300 | Double_t dpt= ptBinEdges[ptBin] - ptBinEdges[ptBin-1]; |
| 301 | |
| 302 | TH1D * hPi0Proj = hPi0->ProjectionX(Form("pt%03i_hPi0",ptBin), imin, imax); |
| 303 | hPi0Proj->Sumw2(); |
| 304 | TH1D * hPi0ProjMix= hPi0Mix->ProjectionX(Form("pt%03i_hPi0Mix",ptBin), imin, imax) ; |
| 305 | hPi0ProjMix->Sumw2(); |
| a9e9e511 |
306 | |
| 47747ea2 |
307 | const Int_t pi0Entries = hPi0Proj->Integral(hPi0Proj->FindBin(lowerMass), hPi0Proj->FindBin(upperMass)); |
| 308 | const Int_t mixEntries = hPi0ProjMix->Integral(hPi0Proj->FindBin(lowerMass), hPi0Proj->FindBin(upperMass)); |
| 6aa36480 |
309 | if( pi0Entries >0) { |
| 310 | hMixRawRatio->SetBinContent(ptBin, mixEntries/pi0Entries); |
| 311 | hMixRawRatio->SetBinError(ptBin, TMath::Sqrt(mixEntries/pi0Entries/pi0Entries + mixEntries*mixEntries/pi0Entries/pi0Entries/pi0Entries)); |
| 312 | } |
| 47747ea2 |
313 | printf("statistics in bin is %i, mixed %i, ", pi0Entries, mixEntries); |
| 314 | if( pi0Entries < 10 ) { |
| 315 | Printf("to few entries"); |
| 316 | continue; |
| 317 | } |
| a9e9e511 |
318 | |
| 47747ea2 |
319 | const bool rebin = pi0Entries < 1000; |
| 320 | if( rebin && pi0Entries >= 500 ) { |
| 321 | printf("rebin by factor 2, "); |
| 322 | hPi0Proj->Rebin(2); |
| 323 | hPi0ProjMix->Rebin(2); |
| 324 | } else if( rebin && pi0Entries >= 200) { |
| 325 | printf("rebin by factor 4, "); |
| 326 | hPi0Proj->Rebin(4); |
| 327 | hPi0ProjMix->Rebin(4); |
| 328 | } else if( rebin && pi0Entries >= 10) { |
| 329 | printf("rebin by factor 5, "); |
| 330 | hPi0Proj->Rebin(5); |
| 331 | hPi0ProjMix->Rebin(5); |
| 332 | } |
| 333 | std::cout << std::endl; |
| a9e9e511 |
334 | |
| 335 | hPi0Proj->SetTitle(Form("M_{#gamma#gamma}, %.1f<p_{T}<%.1f, PID=%s, %s", ptBinEdges[ptBin-1], ptBinEdges[ptBin], outBin.PID().Data(), outBin.Trigger().Data())); |
| 47747ea2 |
336 | hPi0Proj->SetXTitle("M_{#gamma#gamma} (GeV/c^{2})"); |
| a9e9e511 |
337 | hPi0ProjMix->SetTitle(Form("M_{#gamma#gamma}^{Mix}, %.1f<p_{T}<%.1f, PID=%s, %s",ptBinEdges[ptBin-1],ptBinEdges[ptBin],outBin.PID().Data(), outBin.Trigger().Data())); |
| 47747ea2 |
338 | hPi0ProjMix->SetXTitle("M_{#gamma#gamma}^{Mix} (GeV/c^{2})"); |
| 339 | |
| 340 | |
| 341 | // Error Fix |
| 342 | for(Int_t ib=1; ib<=hPi0Proj->GetNbinsX();ib++){if(hPi0Proj ->GetBinContent(ib)==0)hPi0Proj ->SetBinError(ib,1.);} |
| 343 | for(Int_t ib=1; ib<=hPi0Proj->GetNbinsX();ib++){if(hPi0ProjMix->GetBinContent(ib)==0)hPi0ProjMix->SetBinError(ib,1.);} |
| a9e9e511 |
344 | |
| 47747ea2 |
345 | // Signal-Mix Ratio |
| 346 | canvas->cd(2); |
| 347 | TH1D * hPi0Ratio = (TH1D*)hPi0Proj->Clone( Form("pt%03i_hPi0Ratio",ptBin) ) ; |
| 348 | hPi0Ratio->SetTitle(Form("#frac{M_{#gamma#gamma}}{M_{#gamma#gamma}^{Mix}}, %.1f<p_{T}<%.1f GeV/c", ptBinEdges[ptBin-1], ptBinEdges[ptBin])); |
| 349 | hPi0Ratio->Divide(hPi0ProjMix) ; |
| 350 | hPi0Ratio->SetMarkerStyle(20) ; |
| 351 | hPi0Ratio->SetMarkerSize(0.7) ; |
| 352 | |
| 915c7da0 |
353 | if( ptBinEdges[ptBin] < 2.1 ) { |
| 354 | rangeMin = 0.1; |
| 355 | rangeMax = 0.2; |
| 356 | } else if( ptBinEdges[ptBin] < 3.1 ) { |
| 357 | rangeMin = 0.08; |
| 358 | rangeMax = 0.22; |
| 359 | } else if( ptBinEdges[ptBin] < 5.1 ) { |
| 360 | rangeMin = 0.08; |
| 361 | rangeMax = 0.25; |
| 362 | } else if( ptBinEdges[ptBin] < 10.1 ) { |
| 363 | rangeMin = 0.07; |
| 364 | rangeMax = 0.3; |
| 365 | } else { |
| 366 | rangeMin = 0.05; |
| 367 | rangeMax = 0.3; |
| 368 | } |
| a9e9e511 |
369 | |
| 47747ea2 |
370 | // ================================================ |
| 371 | // Fit Pol1 ratio |
| 372 | // ================================================ |
| 373 | printf("Pol1 ratio Fit, "); |
| 374 | canvas->cd(2); |
| a9e9e511 |
375 | |
| 9583e896 |
376 | funcRatioFit1->SetParameters(GetADiffB(hPi0Ratio), 0.136, 0.0055, GetP0(hPi0Ratio), GetP1(hPi0Ratio)) ; |
| 47747ea2 |
377 | funcRatioFit1->SetParLimits(0,0.000,1.000) ; |
| 378 | funcRatioFit1->SetParLimits(1,lowerMass,upperMass) ; |
| 379 | funcRatioFit1->SetParLimits(2,lowerWidth,upperWidth) ; |
| 380 | |
| a9e9e511 |
381 | |
| d1f44d99 |
382 | TFitResultPtr ratioFitResultPtr1 = hPi0Ratio->Fit(funcRatioFit1,"MSQ" ,"",rangeMin,rangeMax) ; |
| 47747ea2 |
383 | if( int(ratioFitResultPtr1) % 4000 ) // "More" error is acceptable |
| d1f44d99 |
384 | ratioFitResultPtr1 = hPi0Ratio->Fit(funcRatioFit1,"MSQ" ,"",rangeMin,rangeMax) ; |
| 47747ea2 |
385 | |
| 386 | Int_t ratioFitError1 = ratioFitResultPtr1; |
| 387 | ratioFitError1 = ratioFitError1 % 4000; // "More" error is acceptable |
| 6aa36480 |
388 | ratioFitError1 = ratioFitError1 || TMath::Abs( funcRatioFit1->GetParameter(1) - lowerMass) < 0.0001 || TMath::Abs( funcRatioFit1->GetParameter(1) - upperMass) < 0.0001;// "center" mass converged to limit |
| 389 | ratioFitError1 = ratioFitError1 || funcRatioFit1->GetParError(1) > (upperMass - lowerMass)/2; // to large of an error |
| 390 | ratioFitError1 = ratioFitError1 || TMath::Abs( funcRatioFit1->GetParameter(2) - lowerWidth) < 0.0001 || TMath::Abs( funcRatioFit1->GetParameter(2) - upperWidth) < 0.0001;// st. error converged to limit |
| 391 | ratioFitError1 = ratioFitError1 || funcRatioFit1->GetParError(2) > (upperWidth - lowerWidth)/2; // to large of an error |
| a9e9e511 |
392 | |
| 47747ea2 |
393 | if( ratioFitError1) { |
| 394 | Printf("in ERROR, %i", ratioFitError1); |
| 6aa36480 |
395 | ((TH1D*) output.GetHistogram("mr1r_error", outBin))->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
| 396 | ((TH1D*) output.GetHistogram("mr1r_error", outBin))->SetBinError (ptBin,fgs->GetParError(1) ) ; |
| 397 | ((TH1D*) output.GetHistogram("sr1r_error", outBin))->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; |
| 398 | ((TH1D*) output.GetHistogram("sr1r_error", outBin))->SetBinError (ptBin,fgs->GetParError(2) ) ; |
| 399 | ((TH1D*) output.GetHistogram("ar1_error", outBin))->SetBinContent(ptBin,funcRatioFit1->GetParameter(3)) ; |
| 400 | ((TH1D*) output.GetHistogram("ar1_error", outBin))->SetBinError (ptBin,funcRatioFit1->GetParError(3)) ; |
| 401 | ((TH1D*) output.GetHistogram("br1_error", outBin))->SetBinContent(ptBin,funcRatioFit1->GetParameter(4)) ; |
| 402 | ((TH1D*) output.GetHistogram("br1_error", outBin))->SetBinError (ptBin,funcRatioFit1->GetParError(4)) ; |
| a9e9e511 |
403 | } |
| 6aa36480 |
404 | if( !ratioFitError1 || ignoreErrors ) { |
| 47747ea2 |
405 | Printf("converged, status:%i, covMatrixStatus: %i", ratioFitResultPtr1->Status(), ratioFitResultPtr1->CovMatrixStatus()); |
| 6aa36480 |
406 | ((TH1D*) output.GetHistogram("mr1r", outBin))->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
| 407 | ((TH1D*) output.GetHistogram("mr1r", outBin))->SetBinError (ptBin,fgs->GetParError(1) ) ; |
| 408 | ((TH1D*) output.GetHistogram("sr1r", outBin))->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; |
| 409 | ((TH1D*) output.GetHistogram("sr1r", outBin))->SetBinError (ptBin,fgs->GetParError(2) ) ; |
| 20fb7b55 |
410 | ((TH1D*) output.GetHistogram("ar1", outBin))->SetBinContent(ptBin,funcRatioFit1->GetParameter(3)) ; |
| 411 | ((TH1D*) output.GetHistogram("ar1", outBin))->SetBinError (ptBin,funcRatioFit1->GetParError(3)) ; |
| 412 | ((TH1D*) output.GetHistogram("br1", outBin))->SetBinContent(ptBin,funcRatioFit1->GetParameter(4)) ; |
| 413 | ((TH1D*) output.GetHistogram("br1", outBin))->SetBinError (ptBin,funcRatioFit1->GetParError(4)) ; |
| 47747ea2 |
414 | } |
| 415 | |
| a9e9e511 |
416 | |
| 417 | |
| 47747ea2 |
418 | // ================================================ |
| 419 | // Fit Pol2 ratio |
| 420 | // ================================================ |
| 421 | printf("Pol2 ratio Fit, "); |
| 422 | if( ratioFitError1 ) { |
| 9583e896 |
423 | funcRatioFit2->SetParameters(GetADiffB(hPi0Ratio), 0.136, 0.0055, GetP0(hPi0Ratio),GetP1(hPi0Ratio), 0.) ; |
| 47747ea2 |
424 | funcRatioFit2->SetParLimits(0,0.000,1.000) ; |
| 425 | funcRatioFit2->SetParLimits(1,lowerMass,upperMass) ; |
| 426 | funcRatioFit2->SetParLimits(2,lowerWidth,upperWidth) ; |
| 427 | } else { |
| 428 | funcRatioFit2->SetParameters(funcRatioFit1->GetParameters()) ; |
| 429 | funcRatioFit2->SetParameter(5, 0); |
| 430 | funcRatioFit2->SetParLimits(0,0.000,1.000) ; |
| 431 | funcRatioFit2->SetParLimits(1,lowerMass,upperMass) ; |
| 432 | funcRatioFit2->SetParLimits(2,lowerWidth,upperWidth) ; |
| 433 | } |
| d1f44d99 |
434 | TFitResultPtr ratioFitResultPtr2 = hPi0Ratio->Fit(funcRatioFit2,"+MSQ" ,"",rangeMin,rangeMax) ; |
| 47747ea2 |
435 | if( int(ratioFitResultPtr2) != 4000 ) // if error, "More" error is acceptable |
| d1f44d99 |
436 | ratioFitResultPtr2 = hPi0Ratio->Fit(funcRatioFit2,"MSQ" ,"",rangeMin,rangeMax) ; |
| 47747ea2 |
437 | |
| 438 | Int_t ratioFitError2 = ratioFitResultPtr2; |
| 439 | ratioFitError2 = ratioFitError2 % 4000; // "More" error is acceptable |
| 6aa36480 |
440 | ratioFitError2 = ratioFitError2 || TMath::Abs( funcRatioFit2->GetParameter(1) - lowerMass) < 0.0001 || TMath::Abs( funcRatioFit2->GetParameter(1) - upperMass) < 0.0001; // "center" mass converged to limit |
| 441 | //ratioFitError2 = ratioFitError2 || funcRatioFit2->GetParError(1) > (upperMass - lowerMass)/2; // to large of an error |
| 442 | ratioFitError2 = ratioFitError2 || TMath::Abs( funcRatioFit2->GetParameter(2) - lowerWidth) < 0.0001 || TMath::Abs( funcRatioFit2->GetParameter(2) - upperWidth) < 0.0001; // st. error converged to limit |
| 443 | //ratioFitError2 = ratioFitError2 || funcRatioFit2->GetParError(2) > (upperWidth - lowerWidth)/2; // to large of an error |
| 47747ea2 |
444 | |
| 445 | if( ratioFitError2) { |
| 446 | Printf("in ERROR, %i", ratioFitError2); |
| 6aa36480 |
447 | ((TH1D*) output.GetHistogram("mr2r_error", outBin))->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
| 448 | ((TH1D*) output.GetHistogram("mr2r_error", outBin))->SetBinError (ptBin,fgs->GetParError(1) ) ; |
| 449 | ((TH1D*) output.GetHistogram("sr2r_error", outBin))->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; |
| 450 | ((TH1D*) output.GetHistogram("sr2r_error", outBin))->SetBinError (ptBin,fgs->GetParError(2) ) ; |
| 451 | ((TH1D*) output.GetHistogram("ar2_error", outBin))->SetBinContent(ptBin,funcRatioFit2->GetParameter(3)) ; |
| 452 | ((TH1D*) output.GetHistogram("ar2_error", outBin))->SetBinError (ptBin,funcRatioFit2->GetParError(3)) ; |
| 453 | ((TH1D*) output.GetHistogram("br2_error", outBin))->SetBinContent(ptBin,funcRatioFit2->GetParameter(4)) ; |
| 454 | ((TH1D*) output.GetHistogram("br2_error", outBin))->SetBinError (ptBin,funcRatioFit2->GetParError(4)) ; |
| 455 | ((TH1D*) output.GetHistogram("cr2_error", outBin))->SetBinContent(ptBin,funcRatioFit2->GetParameter(5)) ; |
| 456 | ((TH1D*) output.GetHistogram("cr2_error", outBin))->SetBinError (ptBin,funcRatioFit2->GetParError(5)) ; |
| a9e9e511 |
457 | } |
| 6aa36480 |
458 | if( !ratioFitError2 || ignoreErrors ) { |
| 47747ea2 |
459 | Printf("converged, status:%i, covMatrixStatus: %i", ratioFitResultPtr2->Status(), ratioFitResultPtr2->CovMatrixStatus()); |
| 6aa36480 |
460 | ((TH1D*) output.GetHistogram("mr2r", outBin))->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
| 461 | ((TH1D*) output.GetHistogram("mr2r", outBin))->SetBinError (ptBin,fgs->GetParError(1) ) ; |
| 462 | ((TH1D*) output.GetHistogram("sr2r", outBin))->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; |
| 463 | ((TH1D*) output.GetHistogram("sr2r", outBin))->SetBinError (ptBin,fgs->GetParError(2) ) ; |
| 20fb7b55 |
464 | ((TH1D*) output.GetHistogram("ar2", outBin))->SetBinContent(ptBin,funcRatioFit2->GetParameter(3)) ; |
| 465 | ((TH1D*) output.GetHistogram("ar2", outBin))->SetBinError (ptBin,funcRatioFit2->GetParError(3)) ; |
| 466 | ((TH1D*) output.GetHistogram("br2", outBin))->SetBinContent(ptBin,funcRatioFit2->GetParameter(4)) ; |
| 467 | ((TH1D*) output.GetHistogram("br2", outBin))->SetBinError (ptBin,funcRatioFit2->GetParError(4)) ; |
| 468 | ((TH1D*) output.GetHistogram("cr2", outBin))->SetBinContent(ptBin,funcRatioFit2->GetParameter(5)) ; |
| 469 | ((TH1D*) output.GetHistogram("cr2", outBin))->SetBinError (ptBin,funcRatioFit2->GetParError(5)) ; |
| 47747ea2 |
470 | } |
| 47747ea2 |
471 | |
| a9e9e511 |
472 | |
| 473 | |
| 47747ea2 |
474 | // ================================================ |
| 475 | // Plot Ratio Fits |
| 476 | // ================================================ |
| 477 | hPi0Ratio->GetXaxis()->SetRangeUser(rangeMin, rangeMax); |
| 478 | hPi0Ratio->Draw(); |
| 479 | canvas->Update(); |
| 480 | |
| 481 | |
| a9e9e511 |
482 | |
| 47747ea2 |
483 | // ================================================ |
| 484 | // Pol1 Scaled Background Subtraction |
| 485 | // ================================================ |
| 486 | canvas->cd(3); |
| 487 | Double_t intRangeMin = PeakPosition(pt)-3.*PeakWidth(pt) ; |
| 488 | Double_t intRangeMax = PeakPosition(pt)+3.*PeakWidth(pt) ; |
| 489 | Int_t intBinMin = hPi0Proj->GetXaxis()->FindBin(intRangeMin) ; |
| 490 | Int_t intBinMax = hPi0Proj->GetXaxis()->FindBin(intRangeMax) ; |
| 491 | Double_t mixInt = hPi0ProjMix->Integral(intBinMin,intBinMax); |
| 492 | |
| 6aa36480 |
493 | if( ! ratioFitError1 || true) { |
| 47747ea2 |
494 | printf("Pol1 BS Fit, "); |
| 495 | TH1D * hPi0MixScaledPol1 = (TH1D*)hPi0ProjMix->Clone(Form("pt%03i_hPi0MixScaledPol1", ptBin)) ; |
| 496 | TH1D * hPi0BSPol1 = (TH1D*)hPi0Proj->Clone(Form("pt%03i_hPi0BSPol1", ptBin)) ; |
| 497 | |
| 498 | // Scale Mix by linear part of ratio, yielding approx background |
| 499 | fbg1->SetParameters(funcRatioFit1->GetParameter(3), funcRatioFit1->GetParameter(4)); |
| 500 | hPi0MixScaledPol1 ->Multiply(fbg1) ; |
| 501 | hPi0BSPol1->Add(hPi0MixScaledPol1 ,-1.) ; |
| 502 | |
| 503 | Int_t binPi0 = hPi0BSPol1->FindBin(funcRatioFit1->GetParameter(1)); |
| 504 | Int_t nWidPi0 = 2 * (Int_t) (funcRatioFit1->GetParameter(2)/hPi0BSPol1->GetBinWidth(1)); |
| 505 | Int_t integral = TMath::Abs(hPi0BSPol1->Integral(binPi0-nWidPi0,binPi0+nWidPi0)); |
| 9583e896 |
506 | fgs->SetParameters(GetADiffB(hPi0BSPol1), funcRatioFit1->GetParameter(1), funcRatioFit1->GetParameter(2)) ; |
| 47747ea2 |
507 | fgs->SetParLimits(0,0.,pi0Entries) ; |
| 508 | fgs->SetParLimits(1,lowerMass,upperMass) ; |
| 509 | fgs->SetParLimits(2,lowerWidth,upperWidth) ; |
| 510 | |
| 511 | // Fit |
| d1f44d99 |
512 | TFitResultPtr bs1FitResultPtr = hPi0BSPol1->Fit(fgs,"MSQ","",rangeMin,rangeMax) ; |
| 47747ea2 |
513 | if( int(bs1FitResultPtr) != 4000 ) // if error, "More" error is acceptable |
| d1f44d99 |
514 | bs1FitResultPtr = hPi0BSPol1->Fit(fgs,"MSQ","",rangeMin,rangeMax) ; |
| a9e9e511 |
515 | |
| 47747ea2 |
516 | Int_t bs1FitError = bs1FitResultPtr; |
| 517 | bs1FitError = bs1FitError % 4000; // "More" error is acceptable |
| 6aa36480 |
518 | bs1FitError = bs1FitError || TMath::Abs( fgs->GetParameter(1) - lowerMass) < 0.0001 || TMath::Abs( fgs->GetParameter(1) - upperMass) < 0.0001; // "center" mass converged to limit |
| 519 | //bs1FitError = bs1FitError || fgs->GetParError(1) > (upperMass - lowerMass)/2; // to large of an error |
| 520 | bs1FitError = bs1FitError || TMath::Abs( fgs->GetParameter(2) - lowerWidth) < 0.0001 || TMath::Abs( fgs->GetParameter(2) - upperWidth) < 0.0001; // st. error converged to limit |
| 521 | //bs1FitError = bs1FitError || fgs->GetParError(2) > (upperWidth - lowerWidth)/2; // to large of an error |
| 47747ea2 |
522 | |
| 6aa36480 |
523 | Double_t y=fgs->GetParameter(0)/hPi0BSPol1->GetXaxis()->GetBinWidth(1) ; |
| 524 | Double_t ey=fgs->GetParError(0)/hPi0BSPol1->GetXaxis()->GetBinWidth(1) ; |
| 525 | Double_t npiInt = hPi0BSPol1->Integral(intBinMin,intBinMax) ; |
| 526 | Double_t norm = fbg1->GetParameter(0) ; |
| 527 | Double_t normErr= fbg1->GetParError(0) ; |
| 47747ea2 |
528 | if( bs1FitError) { |
| 529 | Printf("in ERROR, %i", bs1FitError); |
| 6aa36480 |
530 | ((TH1D*) output.GetHistogram("mr1_error", outBin))->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
| 531 | ((TH1D*) output.GetHistogram("mr1_error", outBin))->SetBinError (ptBin,fgs->GetParError(1) ) ; |
| 532 | ((TH1D*) output.GetHistogram("sr1_error", outBin))->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; |
| 533 | ((TH1D*) output.GetHistogram("sr1_error", outBin))->SetBinError (ptBin,fgs->GetParError(2) ) ; |
| 534 | |
| 535 | ((TH1D*) output.GetHistogram("yr1_error", outBin))->SetBinContent(ptBin,y/dpt) ; |
| 536 | ((TH1D*) output.GetHistogram("yr1_error", outBin))->SetBinError(ptBin,ey/dpt) ; |
| 537 | if(npiInt>0.){ |
| 538 | ((TH1D*) output.GetHistogram("yr1int_error", outBin))->SetBinContent(ptBin,npiInt/dpt) ; |
| 539 | ((TH1D*) output.GetHistogram("yr1int_error", outBin))->SetBinError(ptBin,TMath::Sqrt(npiInt + norm*mixInt + normErr*normErr*mixInt*mixInt + norm*norm*mixInt)/dpt) ; |
| 540 | } |
| a9e9e511 |
541 | } |
| 6aa36480 |
542 | if( !bs1FitError || ignoreErrors ) { |
| 47747ea2 |
543 | Printf("converged, status:%i, covMatrixStatus: %i", bs1FitResultPtr->Status(), bs1FitResultPtr->CovMatrixStatus()); |
| 20fb7b55 |
544 | ((TH1D*) output.GetHistogram("mr1", outBin))->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
| 545 | ((TH1D*) output.GetHistogram("mr1", outBin))->SetBinError (ptBin,fgs->GetParError(1) ) ; |
| 546 | ((TH1D*) output.GetHistogram("sr1", outBin))->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; |
| 547 | ((TH1D*) output.GetHistogram("sr1", outBin))->SetBinError (ptBin,fgs->GetParError(2) ) ; |
| 47747ea2 |
548 | |
| 20fb7b55 |
549 | ((TH1D*) output.GetHistogram("yr1", outBin))->SetBinContent(ptBin,y/dpt) ; |
| 20fb7b55 |
550 | ((TH1D*) output.GetHistogram("yr1", outBin))->SetBinError(ptBin,ey/dpt) ; |
| 47747ea2 |
551 | if(npiInt>0.){ |
| 20fb7b55 |
552 | ((TH1D*) output.GetHistogram("yr1int", outBin))->SetBinContent(ptBin,npiInt/dpt) ; |
| 553 | ((TH1D*) output.GetHistogram("yr1int", outBin))->SetBinError(ptBin,TMath::Sqrt(npiInt + norm*mixInt + normErr*normErr*mixInt*mixInt + norm*norm*mixInt)/dpt) ; |
| 47747ea2 |
554 | // maybe we should use TH1::IntegralAndError |
| 555 | } |
| 556 | } |
| 557 | // Ploting |
| 558 | hPi0BSPol1->SetAxisRange(rangeMin, rangeMax); |
| 559 | hPi0BSPol1->SetMaximum(hPi0BSPol1->GetMaximum()*1.5) ; |
| 560 | hPi0BSPol1->SetMinimum(hPi0BSPol1->GetMinimum()*0.9) ; |
| 561 | hPi0BSPol1->SetMarkerStyle(20) ; |
| 562 | hPi0BSPol1->SetMarkerSize(0.7) ; |
| 563 | hPi0BSPol1->Draw(); |
| 564 | canvas->Update(); |
| 565 | //printf("end pt"); |
| 566 | } |
| a9e9e511 |
567 | |
| 568 | |
| 47747ea2 |
569 | // ================================================ |
| a9e9e511 |
570 | // Pol2 Scaled Background Subtraction |
| 47747ea2 |
571 | // ================================================ |
| 572 | canvas->cd(4); |
| 573 | fbg2->SetParameters(funcRatioFit2->GetParameter(3),funcRatioFit2->GetParameter(4),funcRatioFit2->GetParameter(5)); |
| 6aa36480 |
574 | if( ! ratioFitError2 || true) { |
| 47747ea2 |
575 | printf("Pol1 Scaled Background Subtraction, "); |
| 576 | TH1D * hPi0MixScaledPol2 = (TH1D*)hPi0ProjMix->Clone(Form("pt%03i_hPi0MixScaledPol2", ptBin)) ; |
| 577 | TH1D * hPi0BSPol2 = (TH1D*)hPi0Proj ->Clone(Form("pt%03i_hPi0BSPol2", ptBin)) ; |
| a9e9e511 |
578 | |
| 47747ea2 |
579 | hPi0MixScaledPol2->Multiply(fbg2) ; |
| 580 | hPi0BSPol2 ->Add(hPi0MixScaledPol2,-1.) ; |
| 6aa36480 |
581 | hPi0BSPol2->SetOption(); |
| 47747ea2 |
582 | |
| 583 | Int_t binPi0 = hPi0BSPol2->FindBin(funcRatioFit2->GetParameter(1)); |
| 584 | Int_t nWidPi0 = 2 * (Int_t) (funcRatioFit2->GetParameter(2)/hPi0BSPol2->GetBinWidth(1)); |
| 585 | Int_t integral = TMath::Abs(hPi0BSPol2->Integral(binPi0-nWidPi0,binPi0+nWidPi0)); |
| 9583e896 |
586 | fgs->SetParameters(GetADiffB(hPi0BSPol2), funcRatioFit2->GetParameter(1), funcRatioFit2->GetParameter(2)) ; |
| 47747ea2 |
587 | fgs->SetParLimits(0,0.,pi0Entries) ; |
| 588 | fgs->SetParLimits(1,lowerMass,upperMass) ; |
| 589 | fgs->SetParLimits(2,lowerWidth,upperWidth) ; |
| 590 | |
| 591 | // Fit |
| d1f44d99 |
592 | TFitResultPtr bs2FitResultPtr = hPi0BSPol2->Fit(fgs,"MSQ","",rangeMin,rangeMax) ; |
| 47747ea2 |
593 | if( int(bs2FitResultPtr) != 4000 ) // if error, "More" error is acceptable |
| d1f44d99 |
594 | bs2FitResultPtr = hPi0BSPol2->Fit(fgs,"MSQ","",rangeMin,rangeMax) ; |
| a9e9e511 |
595 | |
| 47747ea2 |
596 | Int_t bs2FitError = bs2FitResultPtr; |
| 597 | bs2FitError = bs2FitError % 4000; // "More" error is acceptable |
| 6aa36480 |
598 | bs2FitError = bs2FitError || TMath::Abs( fgs->GetParameter(1) - lowerMass) < 0.0001 || TMath::Abs( fgs->GetParameter(1) - upperMass) < 0.0001; // "center" mass converged to limit |
| 47747ea2 |
599 | // bs2FitError = bs2FitError || fgs->GetParError(1) > (upperMass - lowerMass)/2; // to large of an error |
| 6aa36480 |
600 | bs2FitError = bs2FitError || TMath::Abs( fgs->GetParameter(2) - lowerWidth) < 0.0001 || TMath::Abs( fgs->GetParameter(2) - upperWidth) < 0.0001; // st. error converged to limit |
| 47747ea2 |
601 | // bs2FitError = bs2FitError || fgs->GetParError(2) > (upperWidth - lowerWidth)/2; // to large of an error |
| 602 | |
| 6aa36480 |
603 | Double_t y=fgs->GetParameter(0)/hPi0BSPol2->GetXaxis()->GetBinWidth(1) ; |
| 604 | Double_t ey=fgs->GetParError(0)/hPi0BSPol2->GetXaxis()->GetBinWidth(1) ; |
| 605 | Double_t npiInt = hPi0BSPol2->Integral(intBinMin,intBinMax) ; |
| 606 | Double_t norm = fbg2->GetParameter(0) ; |
| 607 | Double_t normErr= fbg2->GetParError(0) ; |
| 47747ea2 |
608 | if( bs2FitError ) { |
| 609 | Printf("in ERROR, %i", bs2FitError); |
| 6aa36480 |
610 | ((TH1D*) output.GetHistogram("mr2_error", outBin))->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
| 611 | ((TH1D*) output.GetHistogram("mr2_error", outBin))->SetBinError (ptBin,fgs->GetParError(1) ) ; |
| 612 | ((TH1D*) output.GetHistogram("sr2_error", outBin))->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; |
| 613 | ((TH1D*) output.GetHistogram("sr2_error", outBin))->SetBinError (ptBin,fgs->GetParError(2) ) ; |
| 614 | |
| 615 | ((TH1D*) output.GetHistogram("yr2_error", outBin))->SetBinContent(ptBin,y/dpt) ; |
| 616 | ((TH1D*) output.GetHistogram("yr2_error", outBin))->SetBinError(ptBin,ey/dpt) ; |
| 617 | if(npiInt>0.){ |
| 618 | ((TH1D*) output.GetHistogram("yr2int_error", outBin))->SetBinContent(ptBin,npiInt/dpt) ; |
| 619 | ((TH1D*) output.GetHistogram("yr2int_error", outBin))->SetBinError(ptBin,TMath::Sqrt(npiInt + norm*mixInt + normErr*normErr*mixInt*mixInt + norm*norm*mixInt)/dpt) ; |
| 620 | // maybe we should use TH1::IntegralAndError |
| 621 | } |
| a9e9e511 |
622 | } |
| 6aa36480 |
623 | if( !bs2FitError || ignoreErrors ) { |
| 47747ea2 |
624 | Printf("converged, status:%i, covMatrixStatus: %i", bs2FitResultPtr->Status(), bs2FitResultPtr->CovMatrixStatus()); |
| 20fb7b55 |
625 | ((TH1D*) output.GetHistogram("mr2", outBin))->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
| 626 | ((TH1D*) output.GetHistogram("mr2", outBin))->SetBinError (ptBin,fgs->GetParError(1) ) ; |
| 627 | ((TH1D*) output.GetHistogram("sr2", outBin))->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; |
| 628 | ((TH1D*) output.GetHistogram("sr2", outBin))->SetBinError (ptBin,fgs->GetParError(2) ) ; |
| 47747ea2 |
629 | |
| 20fb7b55 |
630 | ((TH1D*) output.GetHistogram("yr2", outBin))->SetBinContent(ptBin,y/dpt) ; |
| 20fb7b55 |
631 | ((TH1D*) output.GetHistogram("yr2", outBin))->SetBinError(ptBin,ey/dpt) ; |
| 47747ea2 |
632 | if(npiInt>0.){ |
| 20fb7b55 |
633 | ((TH1D*) output.GetHistogram("yr2int", outBin))->SetBinContent(ptBin,npiInt/dpt) ; |
| 634 | ((TH1D*) output.GetHistogram("yr2int", outBin))->SetBinError(ptBin,TMath::Sqrt(npiInt + norm*mixInt + normErr*normErr*mixInt*mixInt + norm*norm*mixInt)/dpt) ; |
| 47747ea2 |
635 | // maybe we should use TH1::IntegralAndError |
| 636 | } |
| 637 | } |
| 638 | // Plotting |
| 639 | hPi0BSPol2->SetAxisRange(rangeMin, rangeMax); |
| 640 | hPi0BSPol2->SetMaximum(hPi0BSPol2->GetMaximum()*1.5) ; |
| 641 | hPi0BSPol2->SetMinimum(hPi0BSPol2->GetMinimum()*0.9) ; |
| 642 | hPi0BSPol2->SetMarkerStyle(20) ; |
| 643 | hPi0BSPol2->SetMarkerSize(0.7) ; |
| 644 | hPi0BSPol2->Draw(); |
| 645 | canvas->Update(); |
| 646 | |
| 647 | } |
| a9e9e511 |
648 | |
| 47747ea2 |
649 | canvas->cd(1); |
| 650 | TH1D * hPi0MixScaled = (TH1D*)hPi0ProjMix ->Clone(Form("%sScaled", hPi0Proj->GetName())) ; |
| 04df14a2 |
651 | int fromBin = hPi0Proj->FindBin(scalarEMBSFrom); |
| 652 | int toBin = TMath::Min(hPi0Proj->FindBin(scalarEMBSTo), hPi0Proj->GetNbinsX()); |
| 653 | Double_t scalar = hPi0Proj->Integral(fromBin, toBin) / hPi0ProjMix->Integral(fromBin, toBin); |
| 654 | hPi0MixScaled->Scale(scalar); |
| 47747ea2 |
655 | hPi0MixScaled->SetLineColor(kRed); |
| 656 | hPi0MixScaled->SetTitle(Form("%s, Scaled", hPi0Proj->GetName())); |
| 657 | hPi0Proj->SetAxisRange(rangeMin, 1.); |
| 47747ea2 |
658 | hPi0Proj->Draw(); |
| 659 | hPi0MixScaled->Draw("same"); |
| a9e9e511 |
660 | |
| 47747ea2 |
661 | canvas->Update(); |
| a9e9e511 |
662 | |
| 20fb7b55 |
663 | canvas->Print(Form("imgs/%s_ptBin:%03i.pdf", outBin.Key().Data(), ptBin)); |
| 664 | canvas->Print(Form("imgs/%s_ptBin:%03i.png", outBin.Key().Data(), ptBin)); |
| a9e9e511 |
665 | |
| 47747ea2 |
666 | std::cout << std::endl; |
| 667 | }// end of Pt slice loop |
| 668 | |
| 669 | |
| 670 | //Normalize by the number of events |
| 671 | Int_t cMin=0, cMax=0; |
| d1f44d99 |
672 | if( input.Bin().Trigger().EqualTo("kCentral") ) |
| a9e9e511 |
673 | switch(outBin.Centrality()) { |
| 26d3ece9 |
674 | // case 0: cMin = 1; cMax = 5; break; |
| 675 | // case 1: cMin = 6; cMax = 10; break; |
| 6aa36480 |
676 | case -1: cMin = 1; cMax = 10; break; |
| 26d3ece9 |
677 | default: Printf("ERROR: cent bin %d not defined for trigger kCentral", outBin.Centrality()); |
| 47747ea2 |
678 | } |
| d1f44d99 |
679 | else if( input.Bin().Trigger().EqualTo("kSemiCentral") ) |
| a9e9e511 |
680 | switch(outBin.Centrality()) { |
| 26d3ece9 |
681 | // case 0: cMin = 11; cMax = 20; break; |
| 682 | // case 1: cMin = 21; cMax = 30; break; |
| 683 | // case 2: cMin = 31; cMax = 40; break; |
| 684 | // case 3: cMin = 41; cMax = 50; break; |
| 6aa36480 |
685 | case -2: cMin = 11; cMax = 20; break; |
| 686 | case -3: cMin = 21; cMax = 30; break; |
| 687 | case -4: cMin = 31; cMax = 40; break; |
| 688 | case -5: cMin = 41; cMax = 50; break; |
| 26d3ece9 |
689 | case -11: cMin = 11; cMax = 50; break; |
| 690 | default: Printf("ERROR: cent bin %d not defined for trigger kSemiCentral", outBin.Centrality()); |
| 47747ea2 |
691 | } |
| d1f44d99 |
692 | else if ( input.Bin().Trigger().EqualTo("kMB") || input.Bin().Trigger().EqualTo("kPHOSPb") ) |
| 20fb7b55 |
693 | switch(outBin.Centrality()) { |
| 26d3ece9 |
694 | // case 0: cMin = 1; cMax = 5; break; |
| 695 | // case 1: cMin = 6; cMax = 10; break; |
| 696 | // case 2: cMin = 11; cMax = 20; break; |
| 697 | // case 3: cMin = 21; cMax = 30; break; |
| 698 | // case 4: cMin = 31; cMax = 40; break; |
| 699 | // case 5: cMin = 41; cMax = 50; break; |
| 700 | // case 6: cMin = 51; cMax = 80; break; |
| 47747ea2 |
701 | case -10: cMin=1; cMax = 80; break; |
| 26d3ece9 |
702 | case -11: cMin = 11; cMax = 50; break; |
| 6aa36480 |
703 | case -1: cMin = 1; cMax = 10; break; |
| 704 | case -2: cMin = 11; cMax = 20; break; |
| 705 | case -3: cMin = 21; cMax = 30; break; |
| 706 | case -4: cMin = 31; cMax = 40; break; |
| 707 | case -5: cMin = 41; cMax = 50; break; |
| 708 | case -6: cMin = 51; cMax = 80; break; |
| 26d3ece9 |
709 | case -7: cMin = 51; cMax = 60; break; |
| 710 | case -8: cMin = 61; cMax = 70; break; |
| 711 | case -9: cMin = 71; cMax = 80; break; |
| 712 | default: Printf("ERROR: cent bin %d not defined for trigger %s", outBin.Centrality(), input.Bin().Trigger().Data()); |
| 47747ea2 |
713 | } |
| 714 | else |
| d1f44d99 |
715 | Printf("ERROR: cent bins not defined for trigger, %s", input.Bin().Trigger().Data()); |
| 47747ea2 |
716 | |
| 717 | Double_t nevents = hCentralityX->Integral(cMin,cMax); |
| 718 | if ( nevents > 0.9 ) { |
| 20fb7b55 |
719 | ((TH1D*) output.GetHistogram("yr1", outBin)) ->Scale(1./nevents) ; |
| 720 | ((TH1D*) output.GetHistogram("yr1int", outBin)) ->Scale(1./nevents) ; |
| 721 | ((TH1D*) output.GetHistogram("yr2", outBin)) ->Scale(1./nevents) ; |
| 722 | ((TH1D*) output.GetHistogram("yr2int", outBin)) ->Scale(1./nevents) ; |
| a9e9e511 |
723 | |
| 6aa36480 |
724 | ((TH1D*) output.GetHistogram("yr1_error", outBin)) ->Scale(1./nevents) ; |
| 725 | ((TH1D*) output.GetHistogram("yr1int_error", outBin)) ->Scale(1./nevents) ; |
| 726 | ((TH1D*) output.GetHistogram("yr2_error", outBin)) ->Scale(1./nevents) ; |
| 727 | ((TH1D*) output.GetHistogram("yr2int_error", outBin)) ->Scale(1./nevents) ; |
| 47747ea2 |
728 | } else { |
| 729 | Printf("WARNING: non positive nEvents in centrality range, cMin:%d, cMax:%d, nEvents:%f", cMin, cMax, nevents ); |
| 730 | |
| 731 | } |
| 732 | |
| a9e9e511 |
733 | |
| 47747ea2 |
734 | // store to file |
| 735 | // TFile* outputFile = TFile::Open(saveToFileName.Data(), "UPDATE"); |
| 736 | //outputList->Write(input.KeySuggestion(), TObject::kSingleKey); |
| 737 | //outputList->Write(); |
| 20fb7b55 |
738 | // outputFile->Write(); |
| 739 | // outputFile->Close(); |
| 740 | // delete outputFile; |
| 47747ea2 |
741 | // delete list and content from memory |
| 742 | //outputList->SetOwner(kTRUE); |
| 743 | //outputList->Delete("slow"); |
| 744 | //delete outputList; |
| 20fb7b55 |
745 | //output.Write(); |
| 6aa36480 |
746 | delete canvas; |
| 47747ea2 |
747 | } |
| 748 | |
| 9583e896 |
749 | double GetA(TH1* hist, double from, double to) |
| 750 | { |
| 751 | int fbin = hist->FindBin(from); |
| 752 | if( from > (hist->GetBinLowEdge(fbin)+hist->GetBinLowEdge(fbin+1)) ) |
| 753 | fbin=fbin+1; |
| 754 | int tbin = hist->FindBin(to); |
| 755 | if( to < (hist->GetBinLowEdge(tbin)+hist->GetBinLowEdge(tbin+1)) ) |
| 756 | tbin=tbin-1; |
| 757 | double integer = hist->Integral(fbin, tbin); |
| 758 | return integer/(tbin-fbin+1); |
| 759 | } |
| 760 | |
| 761 | double GetADiffB(TH1* hist, TH1* mixHist) |
| 762 | { |
| 763 | double rawSig = GetA(hist); |
| 764 | // If we have SCALED Event Mixing Histogram |
| 765 | if(mixHist) { |
| 766 | double mixSig = GetA(mixHist); |
| 767 | double A = TMath::Abs(rawSig - mixSig ); |
| 768 | return A; |
| 769 | } |
| 770 | // else, interpolate between B1 and B2. |
| 771 | double background = GetP0(hist)*(ATo-AFrom) + GetP1(hist)*(ATo*ATo-AFrom*AFrom)/2; |
| 772 | double A = GetA(hist); |
| 773 | return TMath::Abs(A-background); |
| 774 | } |
| 775 | |
| 776 | double GetP0(TH1* hist) |
| 777 | { |
| 778 | double b1=GetA(hist,B1From,B1To)*(B1To-B1From); |
| 779 | double b2=GetA(hist,B2From,B2To)*(B2To-B2From); |
| 780 | double a11 = B1To - B1From; |
| 781 | double a12 = (B1To*B1To - B1From*B1From)/2; |
| 782 | double a21 = B2To - B2From; |
| 783 | double a22 = (B2To*B2To - B2From*B2From)/2; |
| 784 | return (a22*b1-a12*b2)/(a11*a22-a12*a21); |
| 785 | } |
| 786 | double GetP1(TH1* hist) |
| 787 | { |
| 788 | double b1=GetA(hist,B1From,B1To)*(B1To-B1From); |
| 789 | double b2=GetA(hist,B2From,B2To)*(B2To-B2From); |
| 790 | double a11 = B1To - B1From; |
| 791 | double a12 = (B1To*B1To - B1From*B1From)/2; |
| 792 | double a21 = B2To - B2From; |
| 793 | double a22 = (B2To*B2To - B2From*B2From)/2; |
| 794 | return (-a21*b1+a11*b2)/(a11*a22-a12*a21); |
| 795 | } |
| 796 | |
| 47747ea2 |
797 | |
| 798 | //----------------------------------------------------------------------------- |
| 799 | double GetPtBin(int bin){ |
| 800 | // recusive function used by MakePtBins |
| 801 | |
| 802 | if( bin==0 ) |
| 803 | return 1.; |
| 804 | |
| 805 | // return GetPtBin(bin-1) * 1.1; |
| 806 | |
| 807 | // if ( bin % 2 ) |
| 808 | // return GetPtBin(bin-1) + 0.4; |
| 809 | // else |
| 810 | // return GetPtBin(bin-1) + 0.2; |
| 811 | |
| 812 | double previousBin = GetPtBin(bin-1); |
| 813 | double linInc = 1.; |
| 814 | double threshold = 5.; |
| 815 | double logFact = 1 + linInc/threshold; |
| 816 | if ( previousBin < threshold ) // linear |
| 817 | return double(int((previousBin + linInc)*10))/10; |
| 818 | else { // logarithmic |
| 819 | return double(int((previousBin * logFact)*10))/10; |
| 820 | } |
| 821 | } |
| 822 | |
| 823 | //----------------------------------------------------------------------------- |
| 824 | void MakePtBins() { |
| 825 | // function for setting Pt bins |
| 826 | |
| 827 | int bin = -1; |
| 828 | do { |
| 829 | ++bin; |
| 830 | ptBinEdges[bin] = GetPtBin(bin); |
| 831 | } while(ptBinEdges[bin] < 40.); |
| 832 | nPtBins = bin -2; |
| 833 | |
| 834 | printf("Making Pt Bins:\n"); |
| 835 | for(int b=0; b < nPtBins+1; ++b) |
| 836 | printf("%.1f, ", ptBinEdges[b]); |
| 837 | printf("\n N. Bins: %d\n", nPtBins); |
| 838 | |
| 839 | |
| 840 | // for(int bin = 0; bin <= nPtBins; ++bin){ |
| 841 | // ptBinEdges[bin] = GetPtBin(bin); |
| 842 | // printf("%.1f, ", ptBinEdges[bin]); |
| 843 | // } |
| 844 | // printf("\n"); |
| 845 | } |
| 846 | |
| 847 | //----------------------------------------------------------------------------- |
| 848 | Double_t CGausPol1(Double_t * x, Double_t * par){ |
| 849 | //Parameterization of Real/Mixed ratio |
| 850 | Double_t m=par[1] ; |
| 851 | Double_t s=par[2] ; |
| 852 | Double_t dx=(x[0]-m)/s ; |
| 853 | return par[0]*exp(-dx*dx/2.)+par[3]+par[4]*(x[0]-kMean) ; |
| 854 | } |
| 855 | |
| 856 | //----------------------------------------------------------------------------- |
| 857 | Double_t CGausPol2(Double_t * x, Double_t * par){ |
| 858 | //Another parameterization of Real/Mixed ratio7TeV/README |
| 859 | Double_t m=par[1] ; |
| 860 | Double_t s=par[2] ; |
| 861 | Double_t dx=(x[0]-m)/s ; |
| 862 | return par[0]*exp(-dx*dx/2.)+par[3]+par[4]*(x[0]-kMean)+par[5]*(x[0]-kMean)*(x[0]-kMean) ; |
| 863 | } |
| 864 | //----------------------------------------------------------------------------- |
| 865 | Double_t CGausPol0(Double_t * x, Double_t * par){ |
| 866 | //Parameterizatin of signal |
| 867 | Double_t m=par[1] ; |
| 868 | Double_t s=par[2] ; |
| 869 | Double_t dx=(x[0]-m)/s ; |
| 870 | return par[0]*exp(-dx*dx/2.)/TMath::Sqrt(TMath::TwoPi())/s+par[3] ; |
| 871 | } |
| 872 | //----------------------------------------------------------------------------- |
| 873 | Double_t CPol1(Double_t * x, Double_t * par){ |
| 874 | //Normalizatino of Mixed |
| 875 | return par[0]+par[1]*(x[0]-kMean) ; |
| 876 | } |
| 877 | //----------------------------------------------------------------------------- |
| 878 | Double_t CPol2(Double_t * x, Double_t * par){ |
| 879 | //Another normalization of Mixed |
| 880 | return par[0]+par[1]*(x[0]-kMean)+par[2]*(x[0]-kMean)*(x[0]-kMean) ; |
| 881 | } |
| 882 | |
| a9e9e511 |
883 | |
| 47747ea2 |
884 | // Input Definitions |
| 885 | TFile* Input::fFile = 0; |
| a9e9e511 |
886 | Input::Input(const TString& fileName, const RawProduction::TriggerBin& inputBin, TString listPath) |
| d1f44d99 |
887 | : fList(0x0), fBin(inputBin.Trigger()) |
| 47747ea2 |
888 | { |
| 889 | // File |
| 890 | if(fFile && !fileName.EqualTo(fFile->GetName())){ |
| 891 | fFile->Close(); |
| 892 | delete fFile; |
| 893 | fFile = 0x0; |
| 8f90c6d5 |
894 | } |
| 895 | |
| 896 | if(! fFile) { |
| 47747ea2 |
897 | Printf("Opening %s", fileName.Data()); |
| a9e9e511 |
898 | fFile = TFile::Open(fileName.Data(), "READ"); |
| 47747ea2 |
899 | } |
| a9e9e511 |
900 | |
| 47747ea2 |
901 | if( listPath.EqualTo("") ) { |
| 902 | char cstr[256] = ""; |
| d1f44d99 |
903 | sprintf(cstr, "PHOSPi0Flow_%s/PHOSPi0Flow_%sCoutput1", fBin.Trigger().Data(), fBin.Trigger().Data()); |
| 47747ea2 |
904 | listPath = cstr; |
| 905 | } |
| a9e9e511 |
906 | |
| 47747ea2 |
907 | Printf("Getting list, %s", listPath.Data()); |
| 908 | fFile->GetObject(listPath.Data(), fList); |
| a9e9e511 |
909 | if( !fList ) |
| 47747ea2 |
910 | Printf("ERROR: list not found"); |
| 911 | } |
| a9e9e511 |
912 | |
| 47747ea2 |
913 | TH1* Input::GetHistogram(const char* name){ |
| 914 | TObject* obj = fList->FindObject(name); |
| 915 | TH1* hist = dynamic_cast<TH1*> (obj); |
| 916 | if( ! hist) |
| 917 | Printf("MakePi0FitInput::GetHistogram: Error, could not find object of name: %s or cast to hist", name); |
| 918 | return hist; |
| 919 | } |
| 920 | |
| d1f44d99 |
921 | //OutputBin Definitions |
| a53ea5d3 |
922 | TriggerBin::TriggerBin(const TString& trigger) |
| 3374aa5b |
923 | : fTrigger(trigger), fDir(trigger), fKey(trigger) |
| d1f44d99 |
924 | { } |
| 925 | |
| a53ea5d3 |
926 | TriggerBin::TriggerBin(const char* trigger) |
| 3374aa5b |
927 | : fTrigger(trigger), fDir(trigger), fKey(trigger) |
| d1f44d99 |
928 | { } |
| 929 | |
| a53ea5d3 |
930 | TriCenPidBin::TriCenPidBin(Int_t centrality, const TString& pid, const TString& trigger) |
| 931 | : TriggerBin(trigger), fCentrality(centrality), fPID(pid) |
| d1f44d99 |
932 | { |
| 3374aa5b |
933 | fDir.Form("%s/c%03i/%s", trigger.Data(), centrality, pid.Data()); |
| 934 | fKey.Form("%s_c%03i_%s", trigger.Data(), centrality, pid.Data()); |
| d1f44d99 |
935 | } |
| 936 | |
| 937 | Output::Output(const TString& fileName, const char* options) |
| 6aa36480 |
938 | : fFile(0x0) |
| d1f44d99 |
939 | { |
| 940 | fFile = TFile::Open(fileName.Data(), options); |
| d1f44d99 |
941 | } |
| 8375919f |
942 | Output::~Output() |
| 943 | { |
| 4a78c315 |
944 | if( fFile ) |
| 945 | fFile->DeleteAll(); |
| 8375919f |
946 | delete fFile; |
| 947 | } |
| 948 | |
| a9e9e511 |
949 | |
| a53ea5d3 |
950 | void Output::SetDir(const TriggerBin& inBin) |
| d1f44d99 |
951 | { |
| 3374aa5b |
952 | TDirectory* dir = (TDirectoryFile*) fFile; |
| 953 | TStringToken dirs(inBin.Dir(), "/"); |
| 954 | while( dirs.NextToken() ) { |
| 955 | Printf("%s", dirs.Data()); |
| 956 | TDirectory* ndir = dir->GetDirectory(dirs.Data()); |
| 957 | if(!ndir) |
| 958 | dir = dir->mkdir(dirs.Data()); |
| 959 | else |
| 960 | dir = ndir; |
| d1f44d99 |
961 | } |
| 3374aa5b |
962 | dir->cd(); |
| d1f44d99 |
963 | } |
| a9e9e511 |
964 | |
| a53ea5d3 |
965 | TH1* Output::GetHistogram(const TString& name, const RawProduction::TriggerBin& inBin) |
| d1f44d99 |
966 | { |
| 3374aa5b |
967 | TDirectory* dir = fFile->GetDirectory(inBin.Dir().Data(), true); |
| a9e9e511 |
968 | TH1* hist = dynamic_cast<TH1*>( dir->Get(name.Data()) ); |
| d1f44d99 |
969 | if( hist ) |
| 970 | return hist; |
| 971 | else { |
| a53ea5d3 |
972 | Printf("ERROR: Output::GetHistogram: %s could not be found", name.Data()); |
| d1f44d99 |
973 | return 0x0; |
| 974 | } |
| 975 | } |
| baacba22 |
976 | |
| 977 | TH1* Output::GetHistogram(const TString& name) { |
| 978 | TH1* hist = dynamic_cast<TH1*>( fFile->Get(name.Data()) ); |
| 979 | if( hist ) |
| 980 | return hist; |
| 981 | else { |
| 982 | Printf("ERROR: Output::GetHistogram: %s could not be found", name.Data()); |
| 983 | return 0x0; |
| 984 | } |
| 985 | } |
| 986 | |
| 6aa36480 |
987 | |
| 988 | |
| d1f44d99 |
989 | |
| 990 | void Output::Write() |
| 991 | { |
| 6aa36480 |
992 | fFile->Write(); |
| d1f44d99 |
993 | } |
| 994 | |
| 995 | |
| a9e9e511 |
996 | |
| 997 | |
| 47747ea2 |
998 | TH1* GetHistogram_cent(Input& input, const TString& name, int centrality) |
| 999 | { |
| 1000 | // Getter (and merger) for histograms following the naming patern of %s_cen%i |
| 1001 | // |
| 1002 | // For certain negeative values, function is defined to merge centrality bins. |
| 1003 | // -1 0-10% |
| 1004 | // -2 10-20% |
| 1005 | // -3 20-30% |
| 1006 | // -4 30-40% |
| 1007 | // -5 40-50% |
| 04f85d20 |
1008 | // -6 50-80% |
| 1009 | // -7 50-60% |
| 1010 | // -8 60-70% |
| 1011 | // -9 70-80% |
| a9e9e511 |
1012 | // -10 0-80% |
| 04f85d20 |
1013 | // -11 10-50% |
| 47747ea2 |
1014 | if( centrality >= 0 ) { |
| 1015 | char cname[256] = ""; |
| 1016 | sprintf(cname, "%s_cen%i", name.Data(), centrality); |
| 1017 | input.GetHistogram(cname); |
| 1018 | } |
| a9e9e511 |
1019 | |
| 47747ea2 |
1020 | TH1* hist = 0x0; |
| 04f85d20 |
1021 | if( 1 == centBinVersion ) { |
| 1022 | if( input.Bin().Trigger().EqualTo("kMB") || input.Bin().Trigger().EqualTo("kPHOSPb") ) { |
| 1023 | switch(centrality) { |
| 1024 | case -10: hist = MergeHistogram_cent(input, name, centrality, 0, 7); break; |
| 1025 | case -1: hist = MergeHistogram_cent(input, name, centrality, 0, 2); break; |
| 1026 | case -2: hist = MergeHistogram_cent(input, name, centrality, 2, 3); break; |
| 1027 | case -3: hist = MergeHistogram_cent(input, name, centrality, 3, 4); break; |
| 1028 | case -4: hist = MergeHistogram_cent(input, name, centrality, 4, 5); break; |
| 1029 | case -5: hist = MergeHistogram_cent(input, name, centrality, 5, 6); break; |
| 1030 | case -6: hist = MergeHistogram_cent(input, name, centrality, 6, 7); break; |
| 1031 | } |
| 1032 | } else if ( input.Bin().Trigger().EqualTo("kCentral") ) { |
| 1033 | switch( centrality ) { |
| 1034 | case -1: return MergeHistogram_cent(input, name, centrality, 0, 2); break; |
| 1035 | } |
| 1036 | } else if ( input.Bin().Trigger().EqualTo("kSemiCentral") ) { |
| 1037 | switch( centrality ) { |
| 1038 | case -2: return MergeHistogram_cent(input, name, centrality, 0, 1); break; |
| 1039 | case -3: return MergeHistogram_cent(input, name, centrality, 1, 2); break; |
| 1040 | case -4: return MergeHistogram_cent(input, name, centrality, 2, 3); break; |
| 1041 | case -5: return MergeHistogram_cent(input, name, centrality, 3, 4); break; |
| 1042 | } |
| 47747ea2 |
1043 | } |
| 04f85d20 |
1044 | } |
| 1045 | else if ( 2 == centBinVersion ) { |
| 1046 | if( input.Bin().Trigger().EqualTo("kMB") || input.Bin().Trigger().EqualTo("kPHOSPb") ) { |
| 1047 | switch(centrality) { |
| 1048 | case -10: hist = MergeHistogram_cent(input, name, centrality, 0, 8); break; |
| 8a59c037 |
1049 | case -11: hist = MergeHistogram_cent(input, name, centrality, 1, 5); break; |
| 04f85d20 |
1050 | case -1: hist = MergeHistogram_cent(input, name, centrality, 0, 1); break; |
| 1051 | case -2: hist = MergeHistogram_cent(input, name, centrality, 1, 2); break; |
| 1052 | case -3: hist = MergeHistogram_cent(input, name, centrality, 2, 3); break; |
| 1053 | case -4: hist = MergeHistogram_cent(input, name, centrality, 3, 4); break; |
| 1054 | case -5: hist = MergeHistogram_cent(input, name, centrality, 4, 5); break; |
| 1055 | case -6: hist = MergeHistogram_cent(input, name, centrality, 5, 8); break; |
| 1056 | case -7: hist = MergeHistogram_cent(input, name, centrality, 5, 6); break; |
| 1057 | case -8: hist = MergeHistogram_cent(input, name, centrality, 6, 7); break; |
| 1058 | case -9: hist = MergeHistogram_cent(input, name, centrality, 7, 8); break; |
| 1059 | } |
| 1060 | } else if ( input.Bin().Trigger().EqualTo("kCentral") ) { |
| 1061 | switch( centrality ) { |
| 8a59c037 |
1062 | case -1: hist = MergeHistogram_cent(input, name, centrality, 0, 4); break; |
| 04f85d20 |
1063 | } |
| 1064 | } else if ( input.Bin().Trigger().EqualTo("kSemiCentral") ) { |
| 1065 | switch( centrality ) { |
| 8a59c037 |
1066 | case -11:hist = MergeHistogram_cent(input, name, centrality, 0, 8); break; |
| 1067 | case -2: hist = MergeHistogram_cent(input, name, centrality, 0, 5); break; |
| 1068 | case -3: hist = MergeHistogram_cent(input, name, centrality, 5, 6); break; |
| 1069 | case -4: hist = MergeHistogram_cent(input, name, centrality, 6, 7); break; |
| 1070 | case -5: hist = MergeHistogram_cent(input, name, centrality, 7, 8); break; |
| 04f85d20 |
1071 | } |
| 47747ea2 |
1072 | } |
| 04f85d20 |
1073 | |
| 47747ea2 |
1074 | } |
| 1075 | // in case not defined above |
| 1076 | if( ! hist ) { |
| d1f44d99 |
1077 | Printf("ERROR:GetHistogram_cent: %i not possible for %s trigger", centrality, input.Bin().Trigger().Data()); |
| 47747ea2 |
1078 | return 0x0; |
| 1079 | } |
| a9e9e511 |
1080 | |
| 499f0073 |
1081 | hist->SetTitle(Form("%s, %s cent.", hist->GetTitle(), GetCentString(centrality))); |
| 47747ea2 |
1082 | switch(centrality) { |
| 47747ea2 |
1083 | } |
| 1084 | return hist; |
| 1085 | } |
| a9e9e511 |
1086 | |
| 499f0073 |
1087 | const char* GetCentString(int centrality) |
| 1088 | { |
| 1089 | switch(centrality) { |
| 1090 | case -10: return "0-80%"; |
| 1091 | case -11: return "10-50%"; |
| 1092 | case -1: return "0-10%"; |
| 1093 | case -2: return "10-20%"; |
| 1094 | case -3: return "20-30%"; |
| 1095 | case -4: return "30-40%"; |
| 1096 | case -5: return "40-50%"; |
| 1097 | case -6: return "50-80%"; |
| 1098 | case -7: return "50-60%"; |
| 1099 | case -8: return "60-70%"; |
| 1100 | case -9: return "70-80%"; |
| 1101 | default: |
| 3610f0d0 |
1102 | static char cstr[64] =""; |
| 499f0073 |
1103 | sprintf(cstr, "centBin:%i", centrality); |
| 1104 | return cstr; |
| 1105 | } |
| 1106 | } |
| 1107 | |
| 47747ea2 |
1108 | TH1* MergeHistogram_cent(Input& input, const TString& name, int newCentIndex, int fromCentIndex, int toCentIndex) |
| 1109 | { |
| 1110 | // Merger (All cent) for histograms following the naming patern of %s_cen%i, from including to excluding |
| 1111 | // |
| 1112 | // Merges centralites bins into one histogram, from including to excluding, and names the histogram using the patern above. |
| 1113 | // If an histogram with that name Allready exist in the current directory (gDirectory), then no merge |
| 1114 | // occurs and this hist. is simply returned. |
| a9e9e511 |
1115 | |
| 47747ea2 |
1116 | char mergeHistName[256] = ""; |
| 1117 | sprintf(mergeHistName, "%s_cen%i", name.Data(), newCentIndex); |
| a9e9e511 |
1118 | |
| 47747ea2 |
1119 | // Check if histogram allready exists in current directory. |
| 1120 | TH1* mergeHist = dynamic_cast<TH1*>( gDirectory->FindObject(mergeHistName) ); |
| 1121 | if ( mergeHist ) |
| 1122 | return mergeHist; |
| a9e9e511 |
1123 | |
| 47747ea2 |
1124 | // If not so; get the first hist, clone, and add the others |
| 1125 | char cname[256] = ""; |
| 1126 | sprintf(cname, "%s_cen%i", name.Data(), fromCentIndex); |
| 1127 | TH1* hist0 = input.GetHistogram(cname); |
| 1128 | sprintf(cname, "%s_cen%i", name.Data(), newCentIndex); |
| 1129 | TH1 * histMerged = (TH1*) hist0->Clone(cname); |
| a9e9e511 |
1130 | |
| 47747ea2 |
1131 | for(int cent=fromCentIndex+1; cent < toCentIndex; ++cent) { |
| 1132 | sprintf(cname, "%s_cen%i", name.Data(), cent); |
| 1133 | TH1* nextHist = input.GetHistogram(cname); |
| 1134 | if( ! nextHist ) {Printf("ERROR: Merge of histograms failed"); delete histMerged; return 0x0; } |
| 1135 | histMerged->Add(nextHist); |
| 1136 | } |
| 1137 | return histMerged; |
| 1138 | } |
| 1139 | |
| 1140 | |
| a9e9e511 |
1141 | |
| 47747ea2 |
1142 | //----------------------------------------------------------------------------- |
| 1143 | void PPRstyle() |
| 1144 | { |
| 1145 | |
| 1146 | ////////////////////////////////////////////////////////////////////// |
| 1147 | // |
| 1148 | // ROOT style macro for the TRD TDR |
| 1149 | // |
| 1150 | ////////////////////////////////////////////////////////////////////// |
| 1151 | |
| 1152 | gStyle->SetPalette(1); |
| 1153 | gStyle->SetCanvasBorderMode(-1); |
| 1154 | gStyle->SetCanvasBorderSize(1); |
| 1155 | gStyle->SetCanvasColor(10); |
| 1156 | |
| 1157 | gStyle->SetFrameFillColor(10); |
| 1158 | gStyle->SetFrameBorderSize(1); |
| 1159 | gStyle->SetFrameBorderMode(-1); |
| 1160 | gStyle->SetFrameLineWidth(1.2); |
| 1161 | gStyle->SetFrameLineColor(1); |
| 1162 | |
| 1163 | gStyle->SetHistFillColor(0); |
| 1164 | gStyle->SetHistLineWidth(1); |
| 1165 | gStyle->SetHistLineColor(1); |
| 1166 | |
| 1167 | gStyle->SetPadColor(10); |
| 1168 | gStyle->SetPadBorderSize(1); |
| 1169 | gStyle->SetPadBorderMode(-1); |
| 1170 | |
| 1171 | gStyle->SetStatColor(10); |
| 1172 | gStyle->SetTitleColor(kBlack,"X"); |
| 1173 | gStyle->SetTitleColor(kBlack,"Y"); |
| 1174 | |
| 1175 | gStyle->SetLabelSize(0.04,"X"); |
| 1176 | gStyle->SetLabelSize(0.04,"Y"); |
| 1177 | gStyle->SetLabelSize(0.04,"Z"); |
| 1178 | gStyle->SetTitleSize(0.04,"X"); |
| 1179 | gStyle->SetTitleSize(0.04,"Y"); |
| 1180 | gStyle->SetTitleSize(0.04,"Z"); |
| 1181 | gStyle->SetTitleFont(42,"X"); |
| 1182 | gStyle->SetTitleFont(42,"Y"); |
| 1183 | gStyle->SetTitleFont(42,"X"); |
| 1184 | gStyle->SetLabelFont(42,"X"); |
| 1185 | gStyle->SetLabelFont(42,"Y"); |
| 1186 | gStyle->SetLabelFont(42,"Z"); |
| 1187 | gStyle->SetStatFont(42); |
| 1188 | |
| 1189 | gStyle->SetTitleOffset(1.0,"X"); |
| 1190 | gStyle->SetTitleOffset(1.4,"Y"); |
| 1191 | |
| 1192 | gStyle->SetFillColor(kWhite); |
| 1193 | gStyle->SetTitleFillColor(kWhite); |
| 1194 | |
| 1195 | gStyle->SetOptDate(0); |
| 1196 | gStyle->SetOptTitle(1); |
| 1197 | gStyle->SetOptStat(0); |
| 1198 | gStyle->SetOptFit(111); |
| 1199 | |
| 1200 | } |
| 1201 | } |
| 1202 | |
| 1203 | |
| 1204 | void MakeRawProduction() |
| 1205 | { |
| 6aa36480 |
1206 | RawProduction::Output output; |
| 1207 | |
| 1208 | |
| 1209 | //TStringToken triggers("kMB kCentral kSemiCentral kPHOSPb", " "); |
| 1210 | TStringToken triggers("kMB kPHOSPb", " "); |
| 1211 | while(triggers.NextToken()) { |
| a53ea5d3 |
1212 | RawProduction::TriggerBin inBin(triggers); |
| 6aa36480 |
1213 | RawProduction::Input input("AnalysisResults.root", inBin); |
| 1214 | TStringToken pids("All Allcore Allwou Disp Disp2 Dispcore Dispwou CPV CPVcore CPV2 Both Bothcore", " "); |
| 1215 | //TStringToken pids("Bothcore", " "); |
| 1216 | while(pids.NextToken()) { |
| a53ea5d3 |
1217 | RawProduction::TriCenPidBin outBin(-10, pids, inBin.Trigger()); |
| 6aa36480 |
1218 | RawProduction::MakePi0Fit(input, outBin, output); |
| 1219 | } |
| 1220 | } |
| 1221 | output.Write(); |
| 1222 | |
| 1223 | } |
| 1224 | |
| 1225 | void MakeRawProductionAll() |
| 1226 | { |
| 1227 | //gStyle->SetOptStat(0); |
| 1228 | gStyle->SetOptFit(1); |
| 6aa36480 |
1229 | |
| a9e9e511 |
1230 | RawProduction::Output output; |
| 6aa36480 |
1231 | |
| 26d3ece9 |
1232 | TStringToken triggers("kMB kCentral kSemiCentral kPHOSPb", " "); |
| 1233 | //TStringToken triggers("kCentral", " "); |
| 6aa36480 |
1234 | while(triggers.NextToken()) { |
| a9e9e511 |
1235 | RawProduction::TriggerBin triggerBin(triggers); |
| 1236 | RawProduction::Input input("AnalysisResults.root", triggerBin); |
| 1237 | |
| 1238 | RawProduction::MakePi0Fit(input, triggerBin, output); |
| 1239 | |
| 6aa36480 |
1240 | //TStringToken pids("All Allcore Allwou Disp Disp2 Dispcore Dispwou CPV CPVcore CPV2 Both Bothcore", " "); |
| 1241 | //TStringToken pids("All Allcore Allwou Disp Disp2 Dispcore Dispwou", " "); |
| b06cd7b3 |
1242 | //TStringToken pids("All", " "); |
| 1243 | TStringToken pids("All Allcore Allwou Disp Disp2 Dispcore Disp2core Dispwou CPV CPVcore CPV2 CPV2core Both Bothcore Both2 Both2core", " "); |
| 6aa36480 |
1244 | while(pids.NextToken()) { |
| b06cd7b3 |
1245 | for(int cent = -11; cent < 0; ++cent) { |
| a9e9e511 |
1246 | RawProduction::TriCenPidBin tcpBin(cent, pids, triggerBin.Trigger()); |
| b06cd7b3 |
1247 | if(triggers.EqualTo("kMB") || triggers.EqualTo("kPHOSPb")) { |
| 1248 | if( -1 == cent || -11 == cent || -10 == cent || -6 == cent ) |
| 1249 | RawProduction::MakePi0FitTCP(input, tcpBin, output); |
| 1250 | } |
| 1251 | if(triggers.EqualTo("kCentral") ) { |
| 1252 | if( -1 == cent ) |
| 1253 | RawProduction::MakePi0FitTCP(input, tcpBin, output); |
| 1254 | } |
| 1255 | if(triggers.EqualTo("kSemiCentral") ) { |
| 1256 | if( -11 == cent ) |
| 1257 | RawProduction::MakePi0FitTCP(input, tcpBin, output); |
| 1258 | } |
| 1259 | } // cent |
| 1260 | } // pid |
| 1261 | } // trigger |
| 20fb7b55 |
1262 | output.Write(); |
| 47747ea2 |
1263 | } |
| 8375919f |
1264 | |
| 1265 | |
| 1266 | void MakeRawProductionRanges() |
| 1267 | { |
| 1268 | //gStyle->SetOptStat(0); |
| 1269 | gStyle->SetOptFit(1); |
| 1270 | |
| 1271 | float fromRanges[4] = {0.04, 0.05, 0.07, 0.1}; |
| 1272 | float toRanges[4] = {0.2, 0.25, 0.3, 0.4}; |
| 1273 | |
| 4a78c315 |
1274 | TStringToken triggers("kMB kCentral kSemiCentral kPHOSPb", " "); |
| 1275 | while(triggers.NextToken()) { |
| 1276 | RawProduction::TriggerBin triggerBin(triggers); |
| 1277 | RawProduction::Input input("AnalysisResults.root", triggerBin); |
| 1278 | for(int fidx=0; fidx<4; fidx++) { |
| 1279 | for(int tidx=0; tidx<4; tidx++) { |
| 1280 | RawProduction::rangeMin = fromRanges[fidx]; |
| 1281 | RawProduction::rangeMax = toRanges[tidx]; |
| 1282 | |
| 1283 | Printf(" RawProduction_%.2f_%.2f.root", RawProduction::rangeMin, RawProduction::rangeMax); |
| 1284 | RawProduction::Output output(Form("RawProduction_%.2f_%.2f.root", RawProduction::rangeMin, RawProduction::rangeMax)); |
| 8375919f |
1285 | |
| 8375919f |
1286 | |
| 1287 | RawProduction::MakePi0Fit(input, triggerBin, output); |
| 1288 | |
| 1289 | TStringToken pids("All Allcore Allwou Disp Disp2 Dispcore Disp2core Dispwou CPV CPVcore CPV2 CPV2core Both Bothcore Both2 Both2core", " "); |
| 1290 | while(pids.NextToken()) { |
| 1291 | for(int cent = -11; cent < 0; ++cent) { |
| 1292 | RawProduction::TriCenPidBin tcpBin(cent, pids, triggerBin.Trigger()); |
| 1293 | if(triggers.EqualTo("kMB") || triggers.EqualTo("kPHOSPb")) { |
| 1294 | if( -1 == cent || -11 == cent || -10 == cent || -6 == cent ) |
| 1295 | RawProduction::MakePi0FitTCP(input, tcpBin, output); |
| 1296 | } |
| 1297 | if(triggers.EqualTo("kCentral") ) { |
| 1298 | if( -1 == cent ) |
| 1299 | RawProduction::MakePi0FitTCP(input, tcpBin, output); |
| 1300 | } |
| 1301 | if(triggers.EqualTo("kSemiCentral") ) { |
| 1302 | if( -11 == cent ) |
| 1303 | RawProduction::MakePi0FitTCP(input, tcpBin, output); |
| 1304 | } |
| 1305 | } // cent |
| 1306 | } // pid |
| 4a78c315 |
1307 | output.Write(); |
| 1308 | } |
| 8375919f |
1309 | } |
| 4a78c315 |
1310 | }// trigger |
| 8375919f |
1311 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame