AddTaskFemto for train update
[u/mrichter/AliRoot.git] / ITS / GetGainModuleLevel.C
CommitLineData
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
17class TList;\r
18class TH2F;\r
19class TH1F;\r
20class TH1D;\r
21class TProfile;\r
22class TFile;\r
23class TCanvas;\r
24class TStyle;\r
25class TF1;\r
26class TGrid;\r
27\r
28\r
29\r
30\r
31\r
32Double_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
90void 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
105 listin=(TList*)file_data->Get("PWG1dEdxSSDQA/output");\r
106 if(!listin) \r
107 listin=(TList*)file_data->Get("PWG1dEdxSSDQA/SSDdEdxQA");\r
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