end-of-line normalization
[u/mrichter/AliRoot.git] / ITS / GetGainModuleLevel.C
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 }