Another histos for lumi
[u/mrichter/AliRoot.git] / ITS / GetGainModuleLevel.C
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
105                 listin=(TList*)file_data->Get("PWGPPdEdxSSDQA/output");\r
106         if(!listin)     \r
107                 listin=(TList*)file_data->Get("PWGPPdEdxSSDQA/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