]>
Commit | Line | Data |
---|---|---|
db1bbbe1 | 1 | #if !defined(__CINT__) || defined(__MAKECINT__) |
2 | #include "TList.h" | |
3 | #include "TFile.h" | |
4 | #include "TStyle.h" | |
5 | #include "TH1F.h" | |
6 | #include "TH2F.h" | |
7 | #include "THnSparse.h" | |
8 | #include "TLegend.h" | |
9 | #include "TSystem.h" | |
10 | #include "TMath.h" | |
11 | #include "TCanvas.h" | |
12 | #include "TLegend.h" | |
13 | #include "TLatex.h" | |
14 | #include "TF1.h" | |
15 | #include "TLine.h" | |
16 | #include "TPaletteAxis.h" | |
17 | #include "TArrayD.h" | |
18 | #include "TGraphErrors.h" | |
19 | // | |
20 | // | |
21 | #endif | |
22 | ||
23 | ||
24 | TH2 *hMCAccTr,*hMCAccClsS,*hMCAccCls,*hMCGen,*hzvMC2; | |
25 | TH2 *hDtAccTr,*hDtAccClsS,*hDtAccCls,*hzvDt2; | |
26 | TH1 *hzvMC,*hzvDt,*hMCTrue; | |
27 | // | |
28 | TH2 *hTrAcceptance,*hTrDtAccNrm,*hTrDtCorr2D; | |
29 | TH2 *hClAcceptance,*hClDtAccNrm,*hClDtCorr2D; | |
30 | TH1 *hTrDtCorr1D,*hClDtCorr1D; | |
31 | // | |
32 | // for test of MC on MC only | |
33 | TH2 *hMCGenData = 0; | |
34 | TH1 *hMCTrueData = 0; | |
35 | ||
36 | TH1* GetHisto(TList* lst, const char* name, const char* pref=""); | |
37 | void CorrectMult(TH2* hMCAcc, TH2* hMCGn, | |
38 | TH2* hDtAcc, TH1* hzDt, | |
39 | const char* pref, | |
40 | // output | |
41 | TH2 *&hDtAccNrm, // raw data normalized per Z bin | |
42 | TH2 *&hAcceptance, // 2D acceptance | |
43 | TH2 *&hDtCorr2D, // 2D corrected dN/dEta per Z bin (norm up to bin width) | |
44 | TH1 *&hDtCorr1D // 1D corrected dN/dEta | |
45 | ); | |
46 | ||
47 | ||
48 | const double minZStat = 50; | |
49 | ||
50 | void ppcor(const char* flMC, const char* flDt) | |
51 | { | |
52 | TFile *fileMC = TFile::Open(flMC); | |
53 | TList* clistMC = (TList*)fileMC->Get("clist"); | |
54 | TFile *fileDt = TFile::Open(flDt); | |
55 | TList* clistDt = (TList*)fileDt->Get("clist"); | |
56 | // | |
57 | hMCAccTr = (TH2*)GetHisto(clistMC,"b0_TrData_ZvEtaCutT","_mc"); | |
58 | hMCAccClsS = (TH2*)GetHisto(clistMC,"b0_TrData_ZvEtaSPD1","_mc"); | |
59 | hMCAccCls = (TH2*)hMCAccClsS->Clone( Form("MCClusSPD_%s","_mc") ); | |
60 | hMCAccCls->Add(hMCAccTr); | |
61 | hMCGen = (TH2*)GetHisto(clistMC,"b0_zvEtaPrimMC","_mc"); | |
62 | hzvMC2 = (TH2*)GetHisto(clistMC,"zv","_mc"); | |
63 | hzvMC = hzvMC2->ProjectionX("zvgenMC"); | |
64 | // | |
65 | hMCTrue = hMCGen->ProjectionX("MCTrue"); | |
66 | if (hzvMC2->Integral()>0) hMCTrue->Scale(1./hzvMC2->Integral()); | |
67 | hMCTrue->Scale(1./hMCTrue->GetBinWidth(1)); | |
68 | // | |
69 | hDtAccTr = (TH2*)GetHisto(clistDt,"b0_TrData_ZvEtaCutT","_dt"); | |
70 | hDtAccClsS = (TH2*)GetHisto(clistDt,"b0_TrData_ZvEtaSPD1","_dt"); | |
71 | hDtAccCls = (TH2*)hDtAccClsS->Clone( Form("DtClusSPD_%s","_dt") ); | |
72 | hDtAccCls->Add(hDtAccTr); | |
73 | hzvDt2 = (TH2*)GetHisto(clistDt,"zv","_dt"); | |
74 | hzvDt = hzvDt2->ProjectionX("zvDt"); | |
75 | // | |
76 | hMCGenData = (TH2*)GetHisto(clistDt,"b0_zvEtaPrimMC","_data"); | |
77 | if (hMCGenData) { | |
78 | hMCTrueData = hMCGenData->ProjectionX("MCTrueData"); | |
79 | if (hzvDt2->Integral()>0) hMCTrueData->Scale(1./hzvDt2->Integral()); | |
80 | hMCTrueData->Scale(1./hMCTrueData->GetBinWidth(1)); | |
81 | } | |
82 | // | |
83 | CorrectMult(hMCAccTr,hMCGen, | |
84 | hDtAccTr,hzvDt, "trc", | |
85 | hTrDtAccNrm,hTrAcceptance,hTrDtCorr2D,hTrDtCorr1D); | |
86 | // | |
87 | CorrectMult(hMCAccCls,hMCGen, | |
88 | hDtAccCls,hzvDt, "cls", | |
89 | hClDtAccNrm,hClAcceptance,hClDtCorr2D,hClDtCorr1D); | |
90 | // | |
91 | // | |
92 | hTrDtCorr1D->Draw(); // dndeta from tracklets | |
93 | hClDtCorr1D->Draw(); // dndeta from clusters | |
94 | } | |
95 | ||
96 | ||
97 | TH1* GetHisto(TList* lst, const char* name, const char* pref) | |
98 | { | |
99 | TH1* hst = (TH1*)lst->FindObject(name); | |
100 | if (!hst) return 0; | |
101 | TH1* hcl = (TH1*) hst->Clone( Form("%s%ss",hst->GetName(),pref) ); | |
102 | return hcl; | |
103 | } | |
104 | ||
105 | void CorrectMult(TH2* hMCAcc, TH2* hMCGn, | |
106 | TH2* hDtAcc, TH1* hzDt, | |
107 | const char* pref, | |
108 | // output | |
109 | TH2 *&hDtAccNrm, // raw data normalized per Z bin | |
110 | TH2 *&hAcceptance, // 2D acceptance | |
111 | TH2 *&hDtCorr2D, // 2D corrected dN/dEta per Z bin (norm up to bin width) | |
112 | TH1 *&hDtCorr1D // 1D corrected dN/dEta | |
113 | ) | |
114 | { | |
115 | // | |
116 | int neta = hDtAcc->GetXaxis()->GetNbins(); | |
117 | int nz = hDtAcc->GetYaxis()->GetNbins(); | |
118 | // | |
119 | hDtAccNrm = (TH2*) hDtAcc->Clone( Form("rawDataNormZ_%s",pref) ); | |
120 | hDtAccNrm->Reset(); | |
121 | for (int iz=1;iz<=nz;iz++) { // normalize data histo per Z bin | |
122 | double zv = hDtAcc->GetYaxis()->GetBinCenter(iz); | |
123 | int iz0 = hzDt->FindBin(zv); | |
124 | double zst = hzDt->GetBinContent(iz0); | |
125 | double zstE = hzDt->GetBinError(iz0); | |
126 | if (zst<minZStat) continue; | |
127 | for (int ix=1;ix<=neta;ix++) { | |
128 | double vl = hDtAcc->GetBinContent(ix,iz); | |
129 | double er = hDtAcc->GetBinError(ix,iz); | |
130 | if (vl<1e-9 || er<1e-6) continue; | |
131 | double rat = vl/zst; | |
132 | double ratE = rat*TMath::Sqrt( er*er/(vl*vl) + (zstE*zstE)/(zst*zst) ); | |
133 | //printf("%2d %2d %+e(%+e) -> %+e(%+e)\n",ix,iz,vl,er,rat,ratE); | |
134 | hDtAccNrm->SetBinContent(ix,iz,rat); | |
135 | hDtAccNrm->SetBinError(ix,iz,ratE); | |
136 | } | |
137 | } | |
138 | // | |
139 | hAcceptance = (TH2*)hMCAcc->Clone( Form("Acceptance_%s",pref) ); | |
140 | hAcceptance->Divide(hMCGn); | |
141 | // | |
142 | hDtCorr2D = (TH2*)hDtAccNrm->Clone( Form("dNdEtaCor2D_%s",pref) ); | |
143 | hDtCorr2D->Divide(hAcceptance); | |
144 | // | |
145 | hDtCorr1D = hDtCorr2D->ProjectionX( Form("dNdEtaCor1D_%s",pref) ); | |
146 | hDtCorr1D->Reset(); | |
147 | // | |
148 | for (int ix=1;ix<=neta;ix++) { // get weighted average | |
149 | double sm=0,sme=0; | |
150 | double bw = hDtCorr1D->GetBinWidth(ix); | |
151 | for (int iz=1;iz<=nz;iz++) { | |
152 | double vl = hDtCorr2D->GetBinContent(ix,iz); | |
153 | double er = hDtCorr2D->GetBinError(ix,iz); | |
154 | if (vl<1e-6 || er<1e-12) continue; | |
155 | sm += vl/(er*er); | |
156 | sme += 1./(er*er); | |
157 | } | |
158 | if (sme<1e-6) continue; | |
159 | sm /= sme; | |
160 | sme = 1./TMath::Sqrt(sme); | |
161 | hDtCorr1D->SetBinContent(ix,sm/bw); | |
162 | hDtCorr1D->SetBinError(ix,sme/bw); | |
163 | } | |
164 | // | |
165 | } |