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