]>
Commit | Line | Data |
---|---|---|
a65a7e70 | 1 | #include "TMath.h" |
2 | #include "TList.h" | |
3 | #include "TH2F.h" | |
4 | #include "TH1F.h" | |
5 | #include "TF1.h" | |
6 | #include "TH1D.h" | |
7 | #include "TProfile.h" | |
8 | #include "TFile.h" | |
9 | #include "TCanvas.h" | |
10 | #include <Riostream.h> | |
11 | #include <TROOT.h> | |
12 | #include <TString.h> | |
13 | #include "TStyle.h" | |
14 | #include "TGrid.h" | |
15 | ||
16 | ||
17 | class TList; | |
18 | class TH2F; | |
19 | class TH1F; | |
20 | class TH1D; | |
21 | class TProfile; | |
22 | class TFile; | |
23 | class TCanvas; | |
24 | class TStyle; | |
25 | class TF1; | |
26 | class TGrid; | |
27 | ||
28 | ||
29 | ||
30 | ||
31 | ||
32 | Double_t convolution(Double_t* x,Double_t *par) | |
33 | { | |
34 | //Fit parameters: | |
35 | //par[0]=Width (scale) parameter of Landau density | |
36 | //par[1]=Most Probable (MP, location) parameter of Landau density | |
37 | //par[2]=Total area (integral -inf to inf, normalization constant) | |
38 | //par[3]=Width (sigma) of convoluted Gaussian function | |
39 | // | |
40 | //In the Landau distribution (represented by the CERNLIB approximation), | |
41 | //the maximum is located at x=-0.22278298 with the location parameter=0. | |
42 | //This shift is corrected within this function, so that the actual | |
43 | //maximum is identical to the MP parameter. | |
44 | ||
45 | // Numeric constants | |
46 | Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2) | |
47 | Double_t mpshift = -0.22278298; // Landau maximum location | |
48 | ||
49 | // Control constants | |
50 | Double_t np = 200.0; // number of convolution steps | |
51 | Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas | |
52 | ||
53 | // Variables | |
54 | Double_t xx; | |
55 | Double_t mpc; | |
56 | Double_t fland; | |
57 | Double_t sum = 0.0; | |
58 | Double_t xlow,xupp; | |
59 | Double_t step; | |
60 | Double_t i; | |
61 | ||
62 | ||
63 | // MP shift correction | |
64 | mpc = par[1] - mpshift * par[0]; | |
65 | ||
66 | // Range of convolution integral | |
67 | xlow = x[0] - sc * par[3]; | |
68 | xupp = x[0] + sc * par[3]; | |
69 | ||
70 | step = (xupp-xlow) / np; | |
71 | ||
72 | // Convolution integral of Landau and Gaussian by sum | |
73 | for(i=1.0; i<=np/2; i++) | |
74 | { | |
75 | xx = xlow + (i-.5) * step; | |
76 | fland = TMath::Landau(xx,mpc,par[0]) / par[0]; | |
77 | sum += fland * TMath::Gaus(x[0],xx,par[3]); | |
78 | ||
79 | xx = xupp - (i-.5) * step; | |
80 | fland = TMath::Landau(xx,mpc,par[0]) / par[0]; | |
81 | sum += fland * TMath::Gaus(x[0],xx,par[3]); | |
82 | } | |
83 | ||
84 | return (par[2] * step * sum * invsq2pi / par[3]); | |
85 | } | |
86 | ||
87 | ||
88 | ||
89 | ||
90 | void GetGainModuleLevel(TString filename,Bool_t normal=1,Int_t ntofit=500, Bool_t grid=0) | |
91 | { | |
92 | gROOT->SetStyle("Plain"); | |
93 | gStyle->SetPalette(1,0); | |
94 | gStyle->SetOptStat(111111); | |
95 | gStyle->SetOptFit(1); | |
96 | if(grid) | |
97 | TGrid::Connect("alien://"); | |
98 | ||
99 | TFile* file_data=TFile::Open(filename); | |
100 | if(!file_data) | |
101 | return; | |
102 | TList *listin=0x0; | |
103 | listin=(TList*)file_data->Get("output"); | |
104 | if(!listin) | |
105 | listin=(TList*)file_data->Get("PWGPPdEdxSSDQA/output"); | |
106 | if(!listin) | |
107 | listin=(TList*)file_data->Get("PWGPPdEdxSSDQA/SSDdEdxQA"); | |
108 | if(!listin) | |
109 | listin=(TList*)file_data->Get("SSDdEdxQA"); | |
110 | if(!listin) | |
111 | return; | |
112 | TH2F* fHistQ=0x0; | |
113 | if (normal) | |
114 | fHistQ=(TH2F*)listin ->FindObject("QACharge"); | |
115 | else | |
116 | fHistQ=(TH2F*)listin ->FindObject("QAChargeCorrected"); | |
117 | if(!fHistQ) | |
118 | return; | |
119 | TH2F* fHistCR=(TH2F*)listin ->FindObject("QAChargeRatio"); | |
120 | if(!fHistCR) | |
121 | return; | |
122 | TList *listout1=new TList(); | |
123 | ||
124 | TList *listout2=new TList(); | |
125 | ||
126 | TH1F* fHistMPVs=new TH1F("HistMPVS","HistMPVs;MPV;N",75,70,95); | |
127 | fHistMPVs->SetDirectory(0); | |
128 | listout2->Add(fHistMPVs); | |
129 | ||
130 | TH1F* fHistSL=new TH1F("HistSL","#sigma_{Landau};#sigma_{Landau};N",40,0,16); | |
131 | fHistSL->SetDirectory(0); | |
132 | listout2->Add(fHistSL); | |
133 | ||
134 | TH1F* fHistSG=new TH1F("HistSG","#sigma_{Gaus};#sigma_{Gaus};N",40,0,16); | |
135 | fHistSG->SetDirectory(0); | |
136 | listout2->Add(fHistSG); | |
137 | ||
138 | TH1F* fHistCRmean=new TH1F("HistCRmean","HistCRmean;CRmean;N",200,-1,1); | |
139 | fHistCRmean->SetDirectory(0); | |
140 | listout2->Add(fHistCRmean); | |
141 | ||
142 | TH1F* fHistCRRMS=new TH1F("HistCRRMS","HistCRRMS;CRRMS;N",100,0,1); | |
143 | fHistCRRMS->SetDirectory(0); | |
144 | listout2->Add(fHistCRRMS); | |
145 | ||
146 | TH1F* fHistGainP=new TH1F("HistGainP","HistGainP;CRGainPcorr;N",120,0.5,2.0); | |
147 | fHistGainP->SetDirectory(0); | |
148 | listout2->Add(fHistGainP); | |
149 | ||
150 | TH1F* fHistGainN=new TH1F("HistGainN","HistGainN;CRGainNcorr;N",120,0.5,2.0); | |
151 | fHistGainN->SetDirectory(0); | |
152 | listout2->Add(fHistGainN); | |
153 | ||
154 | TH1F *fMPVGraph = new TH1F("MPVgraph","MPVgraph;Module number;MPV",1698,-0.5,1697.5); | |
155 | fMPVGraph->SetMarkerColor(kRed); | |
156 | fMPVGraph->SetMarkerStyle(22); | |
157 | listout2->Add(fMPVGraph); | |
158 | ||
159 | TH1F *fCRmeanGraph = new TH1F("CRmeangraph","CRmeangraph;Module number;MPV",1698,-0.5,1697.5); | |
160 | fCRmeanGraph->SetMarkerColor(kBlue); | |
161 | fCRmeanGraph->SetMarkerStyle(23); | |
162 | listout2->Add(fCRmeanGraph); | |
163 | ||
164 | Float_t gainP[1698]; | |
165 | Float_t gainN[1698]; | |
166 | Float_t mpv[1698]; | |
167 | Int_t flag[1698]; | |
168 | ||
169 | ofstream outfiletxt; | |
170 | outfiletxt.open("gain.txt"); | |
171 | outfiletxt.width(10) ; | |
172 | outfiletxt.setf(outfiletxt.left); | |
173 | outfiletxt<<"MODULE"<<"\t"; | |
174 | outfiletxt.width(10); | |
175 | outfiletxt.setf(outfiletxt.left); | |
176 | outfiletxt<<"FLAG"<<"\t"; | |
177 | outfiletxt.width(10) ; | |
178 | outfiletxt.setf(outfiletxt.left); | |
179 | outfiletxt<<"GainPcorr"<<"\t"; | |
180 | outfiletxt.width(10) ; | |
181 | outfiletxt.setf(outfiletxt.left); | |
182 | outfiletxt<<"GainNcorr"<<"\t"; | |
183 | outfiletxt.width(10) ; | |
184 | outfiletxt.setf(outfiletxt.left); | |
185 | outfiletxt<<"MPV"<<endl; | |
186 | ||
187 | ||
188 | ofstream outfiletxtbad; | |
189 | outfiletxtbad.open("badModules.txt"); | |
190 | ||
191 | ||
192 | ||
193 | ||
194 | for (int i =0;i<1698;i++) | |
195 | { | |
196 | cout<<i<<endl; | |
197 | TString tmpQ("Q"); | |
198 | tmpQ+=i; | |
199 | TString tmpCR("CR"); | |
200 | tmpCR+=i; | |
201 | TH1D* fHist1DCR= fHistCR->ProjectionY(tmpCR,i+1,i+1); | |
202 | Double_t mean=fHist1DCR->GetMean(); | |
203 | if(!(TMath::Abs(mean)<1.0)||fHist1DCR->GetEntries()<10) | |
204 | { | |
205 | flag[i]=-2; | |
206 | gainN[i]=1.0; | |
207 | gainP[i]=1.0; | |
208 | mpv[i]=1.0; | |
209 | continue; | |
210 | } | |
211 | fHistCRmean->Fill(mean); | |
212 | fHistCRRMS->Fill(fHist1DCR->GetRMS()); | |
213 | gainN[i]=1.0/(1.0+mean); | |
214 | gainP[i]=1.0/(1.0-mean); | |
215 | fHistGainP->Fill(gainP[i]); | |
216 | fHistGainN->Fill(gainN[i]); | |
217 | fCRmeanGraph->SetBinContent(i+1,mean); | |
218 | fCRmeanGraph->SetBinError(i+1,fHist1DCR->GetRMS()); | |
219 | ||
220 | ||
221 | TH1D* fHist1DQ=fHistQ->ProjectionY(tmpQ,i+1,i+1); | |
222 | fHist1DQ->SetDirectory(0); | |
223 | listout1->Add(fHist1DQ); | |
224 | if(fHist1DQ->GetEntries()<ntofit) | |
225 | { | |
226 | flag[i]=-1; | |
227 | mpv[i]=1.0; | |
228 | outfiletxtbad<<"Low statistic \t module= "<<i<<" netries="<<fHist1DQ->GetEntries()<<endl; | |
229 | continue; | |
230 | } | |
231 | else | |
232 | { | |
233 | tmpQ+="fit"; | |
234 | Float_t range=fHist1DQ->GetBinCenter(fHist1DQ->GetMaximumBin()); | |
235 | TF1 *f1 = new TF1(tmpQ,convolution,range*0.45,range*3.0,4); | |
236 | f1->SetParameters(7.0,range,1.0,5.5); | |
237 | Float_t normalization=fHist1DQ->GetEntries()*fHist1DQ->GetXaxis()->GetBinWidth(2)/f1->Integral(range*0.45,range*3.0); | |
238 | f1->SetParameters(7.0,range,normalization,5.5); | |
239 | //f1->SetParameters(7.0,range,fHist1DQ->GetMaximum(),5.5); | |
240 | f1->SetParNames("sigma Landau","MPV","N","sigma Gaus"); | |
241 | f1->SetParLimits(0,2.0,100.0); | |
242 | f1->SetParLimits(3,0.0,100.0); | |
243 | if(fHist1DQ->Fit(tmpQ,"BRQ")==0) | |
244 | { | |
245 | mpv[i]=f1->GetParameter(1); | |
246 | fHistMPVs->Fill(mpv[i]); | |
247 | fHistSL->Fill(f1->GetParameter(0)); | |
248 | fHistSG->Fill(f1->GetParameter(3)); | |
249 | flag[i]=1; | |
250 | fMPVGraph->SetBinContent(i+1,f1->GetParameter(1)); | |
251 | fMPVGraph->SetBinError(i+1,f1->GetParError(1)); | |
252 | if(mpv[i]<75.0) | |
253 | { | |
254 | outfiletxtbad<<"MPV lower than 75 \t module="<<i<<endl; | |
255 | flag[i]=0; | |
256 | } | |
257 | if(mpv[i]>100.0) | |
258 | { | |
259 | outfiletxtbad<<"MPV higher than 100 \t module="<<i<<endl; | |
260 | flag[i]=0; | |
261 | ||
262 | } | |
263 | if(f1->GetParError(1)>1.0) | |
264 | { | |
265 | outfiletxtbad<<"MPV high error on MPV \t module="<<i<<endl; | |
266 | flag[i]=0; | |
267 | } | |
268 | } | |
269 | else | |
270 | { | |
271 | mpv[i]=1; | |
272 | flag[i]=0; | |
273 | outfiletxtbad<<"BAD FIT \t module="<<i<<endl; | |
274 | //continue; | |
275 | } | |
276 | } | |
277 | } | |
278 | ||
279 | for (int i=0;i<1698;i++) | |
280 | { | |
281 | outfiletxt.setf(outfiletxt.scientific); | |
282 | outfiletxt.precision(2); | |
283 | outfiletxt.width(10) ; | |
284 | outfiletxt.setf(outfiletxt.left); | |
285 | outfiletxt<<i<<"\t"; | |
286 | outfiletxt.width(10) ; | |
287 | outfiletxt.setf(outfiletxt.left); | |
288 | outfiletxt<<flag[i]<<"\t"; | |
289 | outfiletxt.width(10) ; | |
290 | outfiletxt.setf(outfiletxt.left); | |
291 | outfiletxt<<gainP[i]<<"\t"; | |
292 | outfiletxt.width(10) ; | |
293 | outfiletxt.setf(outfiletxt.left); | |
294 | outfiletxt<<gainN[i]<<"\t"; | |
295 | outfiletxt.width(10) ; | |
296 | outfiletxt.setf(outfiletxt.left); | |
297 | outfiletxt<<mpv[i]<<endl; | |
298 | } | |
299 | ||
300 | TCanvas *c1 = new TCanvas("1","1",1200,800); | |
301 | c1->Divide(2,1); | |
302 | c1->cd(1); | |
303 | fHistQ->Draw("colz"); | |
304 | c1->cd(2); | |
305 | fHistCR->Draw("colz"); | |
306 | ||
307 | ||
308 | ||
309 | TFile* fout1=TFile::Open("gain_all_fits.root","recreate"); | |
310 | listout1->Write("output",TObject::kSingleKey); | |
311 | fout1->Close(); | |
312 | ||
313 | TFile* fout2=TFile::Open("gain_control_plots.root","recreate"); | |
314 | listout2->Write("output",TObject::kSingleKey); | |
315 | fout2->Close(); | |
316 | ||
317 | ||
318 | ||
319 | } |