[u/mrichter/AliRoot.git] / ITS / GetGainModuleLevel.C

#include "TMath.h"
#include "TList.h"
#include "TH2F.h"
#include "TH1F.h"
#include "TF1.h"
#include "TH1D.h"
#include "TProfile.h"
#include "TFile.h"
#include "TCanvas.h"
#include <Riostream.h>
#include <TROOT.h>
#include <TString.h>
#include "TStyle.h"
#include "TGrid.h"
 
 
class TList;
class TH2F;
class TH1F;
class TH1D;
class TProfile;
class TFile;
class TCanvas;
class TStyle;
class TF1;
class TGrid;


Double_t convolution(Double_t* x,Double_t *par)
{
	  //Fit parameters:
   //par[0]=Width (scale) parameter of Landau density
   //par[1]=Most Probable (MP, location) parameter of Landau density
   //par[2]=Total area (integral -inf to inf, normalization constant)
   //par[3]=Width (sigma) of convoluted Gaussian function
   //
   //In the Landau distribution (represented by the CERNLIB approximation),
   //the maximum is located at x=-0.22278298 with the location parameter=0.
   //This shift is corrected within this function, so that the actual
   //maximum is identical to the MP parameter.

      // Numeric constants
      Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
      Double_t mpshift  = -0.22278298;       // Landau maximum location

      // Control constants
      Double_t np = 200.0;      // number of convolution steps
      Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas

      // Variables
      Double_t xx;
      Double_t mpc;
      Double_t fland;
      Double_t sum = 0.0;
      Double_t xlow,xupp;
      Double_t step;
      Double_t i;


      // MP shift correction
      mpc = par[1] - mpshift * par[0];

      // Range of convolution integral
      xlow = x[0] - sc * par[3];
      xupp = x[0] + sc * par[3];

      step = (xupp-xlow) / np;

      // Convolution integral of Landau and Gaussian by sum
      for(i=1.0; i<=np/2; i++) 
      {
        	 xx = xlow + (i-.5) * step;
         	fland = TMath::Landau(xx,mpc,par[0]) / par[0];
         	sum += fland * TMath::Gaus(x[0],xx,par[3]);

         	xx = xupp - (i-.5) * step;
         	fland = TMath::Landau(xx,mpc,par[0]) / par[0];
         	sum += fland * TMath::Gaus(x[0],xx,par[3]);
      }

      return (par[2] * step * sum * invsq2pi / par[3]);
}


void GetGainModuleLevel(TString filename,Bool_t normal=1,Int_t ntofit=500, Bool_t grid=0)
{
	gROOT->SetStyle("Plain");
	gStyle->SetPalette(1,0);
	gStyle->SetOptStat(111111);
	gStyle->SetOptFit(1);	
	if(grid)
		TGrid::Connect("alien://");
	
	TFile* file_data=TFile::Open(filename);
	if(!file_data)
		return;
	TList *listin=0x0;
	listin=(TList*)file_data->Get("output");
	if(!listin)
		listin=(TList*)file_data->Get("PWGPPdEdxSSDQA/output");
	if(!listin)	
		listin=(TList*)file_data->Get("PWGPPdEdxSSDQA/SSDdEdxQA");
	if(!listin)	
		listin=(TList*)file_data->Get("SSDdEdxQA");
	if(!listin)	
		return;
	TH2F* fHistQ=0x0;
	if (normal)		
		fHistQ=(TH2F*)listin ->FindObject("QACharge");
	else
		fHistQ=(TH2F*)listin ->FindObject("QAChargeCorrected");
	if(!fHistQ)
		return;
	TH2F* fHistCR=(TH2F*)listin ->FindObject("QAChargeRatio");
	if(!fHistCR)
		return;
	TList *listout1=new TList();
	
	TList *listout2=new TList();
	
	TH1F* fHistMPVs=new TH1F("HistMPVS","HistMPVs;MPV;N",75,70,95);
	fHistMPVs->SetDirectory(0);
	listout2->Add(fHistMPVs);
	
	TH1F* fHistSL=new TH1F("HistSL","#sigma_{Landau};#sigma_{Landau};N",40,0,16);
	fHistSL->SetDirectory(0);
	listout2->Add(fHistSL);
	
	TH1F* fHistSG=new TH1F("HistSG","#sigma_{Gaus};#sigma_{Gaus};N",40,0,16);
	fHistSG->SetDirectory(0);
	listout2->Add(fHistSG);
	
	TH1F* fHistCRmean=new TH1F("HistCRmean","HistCRmean;CRmean;N",200,-1,1);
	fHistCRmean->SetDirectory(0);
	listout2->Add(fHistCRmean);
	
	TH1F* fHistCRRMS=new TH1F("HistCRRMS","HistCRRMS;CRRMS;N",100,0,1);
	fHistCRRMS->SetDirectory(0);
	listout2->Add(fHistCRRMS);
	
	TH1F* fHistGainP=new TH1F("HistGainP","HistGainP;CRGainPcorr;N",120,0.5,2.0);
	fHistGainP->SetDirectory(0);
	listout2->Add(fHistGainP);
	
	TH1F* fHistGainN=new TH1F("HistGainN","HistGainN;CRGainNcorr;N",120,0.5,2.0);
	fHistGainN->SetDirectory(0);
	listout2->Add(fHistGainN);
	
	TH1F *fMPVGraph = new TH1F("MPVgraph","MPVgraph;Module number;MPV",1698,-0.5,1697.5);
	fMPVGraph->SetMarkerColor(kRed);
	fMPVGraph->SetMarkerStyle(22);
	listout2->Add(fMPVGraph);
	
	TH1F *fCRmeanGraph = new TH1F("CRmeangraph","CRmeangraph;Module number;MPV",1698,-0.5,1697.5);
	fCRmeanGraph->SetMarkerColor(kBlue);
	fCRmeanGraph->SetMarkerStyle(23);
	listout2->Add(fCRmeanGraph);
	
	Float_t gainP[1698];
	Float_t gainN[1698];
	Float_t mpv[1698];
	Int_t flag[1698];	
	
	ofstream outfiletxt;
	outfiletxt.open("gain.txt");
	 outfiletxt.width(10) ;
	outfiletxt.setf(outfiletxt.left);
	outfiletxt<<"MODULE"<<"\t";
	outfiletxt.width(10);
	outfiletxt.setf(outfiletxt.left);
	outfiletxt<<"FLAG"<<"\t";
	outfiletxt.width(10) ;
	outfiletxt.setf(outfiletxt.left);
	outfiletxt<<"GainPcorr"<<"\t";
	outfiletxt.width(10) ;
	outfiletxt.setf(outfiletxt.left);
	outfiletxt<<"GainNcorr"<<"\t";
	outfiletxt.width(10) ;
	outfiletxt.setf(outfiletxt.left);
	outfiletxt<<"MPV"<<endl;
	
	
	ofstream outfiletxtbad;
	outfiletxtbad.open("badModules.txt");
	
	
	for (int i =0;i<1698;i++)
	{
		cout<<i<<endl;
		TString tmpQ("Q");
		tmpQ+=i;
		TString tmpCR("CR");
		tmpCR+=i;
		TH1D* fHist1DCR= fHistCR->ProjectionY(tmpCR,i+1,i+1);
		Double_t mean=fHist1DCR->GetMean();
		if(!(TMath::Abs(mean)<1.0)||fHist1DCR->GetEntries()<10)
		{		
			flag[i]=-2;
			gainN[i]=1.0;
			gainP[i]=1.0;
			mpv[i]=1.0;
			continue;
		}
		fHistCRmean->Fill(mean);
		fHistCRRMS->Fill(fHist1DCR->GetRMS());
		gainN[i]=1.0/(1.0+mean);
		gainP[i]=1.0/(1.0-mean);
		fHistGainP->Fill(gainP[i]);
		fHistGainN->Fill(gainN[i]);
		fCRmeanGraph->SetBinContent(i+1,mean);
		fCRmeanGraph->SetBinError(i+1,fHist1DCR->GetRMS());
		
		
		TH1D* fHist1DQ=fHistQ->ProjectionY(tmpQ,i+1,i+1);
		 fHist1DQ->SetDirectory(0);
		 listout1->Add(fHist1DQ);
		if(fHist1DQ->GetEntries()<ntofit)
		{
			flag[i]=-1;
			mpv[i]=1.0;
			 outfiletxtbad<<"Low statistic \t module= "<<i<<" netries="<<fHist1DQ->GetEntries()<<endl;
			continue;
		}
		else
		{
			tmpQ+="fit";
			Float_t range=fHist1DQ->GetBinCenter(fHist1DQ->GetMaximumBin());
			TF1 *f1 = new TF1(tmpQ,convolution,range*0.45,range*3.0,4);
			f1->SetParameters(7.0,range,1.0,5.5);
			Float_t normalization=fHist1DQ->GetEntries()*fHist1DQ->GetXaxis()->GetBinWidth(2)/f1->Integral(range*0.45,range*3.0);
			f1->SetParameters(7.0,range,normalization,5.5);
			//f1->SetParameters(7.0,range,fHist1DQ->GetMaximum(),5.5);
			f1->SetParNames("sigma Landau","MPV","N","sigma Gaus");
			f1->SetParLimits(0,2.0,100.0);
			f1->SetParLimits(3,0.0,100.0);
			if(fHist1DQ->Fit(tmpQ,"BRQ")==0)
			{
				mpv[i]=f1->GetParameter(1);
				fHistMPVs->Fill(mpv[i]);
				fHistSL->Fill(f1->GetParameter(0));
				fHistSG->Fill(f1->GetParameter(3));
				flag[i]=1;		
				fMPVGraph->SetBinContent(i+1,f1->GetParameter(1));
				fMPVGraph->SetBinError(i+1,f1->GetParError(1));
				if(mpv[i]<75.0)
				{
					outfiletxtbad<<"MPV lower than 75 \t module="<<i<<endl;
					flag[i]=0;
				}	
				if(mpv[i]>100.0)
				{
					outfiletxtbad<<"MPV higher than 100 \t module="<<i<<endl;
					flag[i]=0;
					
				}
				if(f1->GetParError(1)>1.0)
				{
					outfiletxtbad<<"MPV high error on MPV  \t module="<<i<<endl;				
					flag[i]=0;
				}
			}
			else
			{
				mpv[i]=1;
				flag[i]=0;
				 outfiletxtbad<<"BAD FIT \t module="<<i<<endl;
				//continue;
			}	
		}	
	}	
	
	for (int i=0;i<1698;i++)
	{	
		outfiletxt.setf(outfiletxt.scientific);
		outfiletxt.precision(2);	
		 outfiletxt.width(10) ;
		outfiletxt.setf(outfiletxt.left);
		outfiletxt<<i<<"\t";
		 outfiletxt.width(10) ;
		outfiletxt.setf(outfiletxt.left);
		outfiletxt<<flag[i]<<"\t";
		 outfiletxt.width(10) ;
		outfiletxt.setf(outfiletxt.left);
		outfiletxt<<gainP[i]<<"\t";
		 outfiletxt.width(10) ;
		outfiletxt.setf(outfiletxt.left);
		outfiletxt<<gainN[i]<<"\t";
		 outfiletxt.width(10) ;
		outfiletxt.setf(outfiletxt.left);
		outfiletxt<<mpv[i]<<endl;
	}
		 
	TCanvas *c1 = new TCanvas("1","1",1200,800);
	c1->Divide(2,1);
	c1->cd(1);
	fHistQ->Draw("colz");
	c1->cd(2);
	fHistCR->Draw("colz");
	
	
	TFile* fout1=TFile::Open("gain_all_fits.root","recreate");
	listout1->Write("output",TObject::kSingleKey);	
	fout1->Close();
	
	TFile* fout2=TFile::Open("gain_control_plots.root","recreate");
	listout2->Write("output",TObject::kSingleKey);	
	fout2->Close();
	
	
}
Commit	Line	Data
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
	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,range0.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(range0.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	}