avoid histos tied to directory to avoid their automatic deletion
[u/mrichter/AliRoot.git] / ITS / GetGainModuleLevel.C
CommitLineData
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
17class TList;
18class TH2F;
19class TH1F;
20class TH1D;
21class TProfile;
22class TFile;
23class TCanvas;
24class TStyle;
25class TF1;
26class TGrid;
27
28
29
30
31
32Double_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
90void 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}