]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/macros/PlotImpParResVsPt.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix.
[u/mrichter/AliRoot.git] / PWGPP / macros / PlotImpParResVsPt.C
CommitLineData
db287dca 1TF1 *CreateFuncGausTail(TH1F *h,TString funcname);\r
2TF1 *CreateFuncGaussTail(TH1F *h,TString funcname,Double_t parone,Double_t partwo);\r
3\r
4void PlotImpParResVsPt()\r
5{\r
6 gStyle->SetOptStat(0);\r
7\r
8\r
9 TFile *file1=TFile::Open("QAresults_LHC11h_p2_170593.root","READ");\r
10 TDirectoryFile *dir1=(TDirectoryFile*)file1->GetDirectory("ImpParRes_Performance");\r
11 \r
12 //TFile *d0dist1=TFile::Open("/home/luojb/2011/ImpactParameter/OutputFile/PP/LHC10b/ImpParRes.Performance_LHC10bcde_pass2_ESDdiamond_2011_01_26.root","READ");\r
13 \r
14 //if (!d0dist1->IsOpen()) return;\r
15 \r
16 \r
17\r
18 TList *d0allpointRec1 = (TList*)dir1->Get("coutputd0allPointRec_0_1000000");\r
19 \r
20 TList *d0allpointSkip1 = (TList*)dir1->Get("coutputd0allPointSkip_0_1000000");\r
21 \r
22 TList *d0particlePID1 = (TList*)dir1->Get("coutputd0particlePID_0_1000000");\r
23 \r
24 TList *d0pt = (TList*)dir1->Get("coutputd0Pt_0_1000000");\r
25 //TFile *d0ppr=TFile::Open("../ResolutionsAnalysis_def.root","READ"); \r
26 //if(!d0ppr->IsOpen()) return;\r
27 \r
28 //define the histgram\r
29 \r
30 TH1F **d0AllpointrphiSkip1_=new TH1F*[26]; \r
31 \r
32 TH1F **d0AllpointrphiSkipTail1_=new TH1F*[26]; \r
33 TH1F **d0AllpointrphiSkipGaus1_=new TH1F*[26]; \r
34 \r
35 TH1F **d0Pt_=new TH1F*[26];\r
36 Int_t nbins = 51; \r
37 Double_t xbins[52]={0.22,0.23,0.26,0.27,0.35,0.36,0.45,0.46,0.55,0.56,0.65,0.66,0.75,0.76,0.85,0.86,1.05,1.06,1.25,1.27,1.45,1.47,1.65,1.67,1.85,1.87,2.15,2.17,2.45,2.47,2.65,2.67,2.85,2.87,3.25,3.27,3.75,3.77,4.15,4.20,4.95,5.15,5.35,5.55,6.0,6.8,8.5,10.5,12,19.,21.,32};\r
38 \r
39 //TGraph *Resrphi_allpoint =new TGraph("Resrphi_allpot","d_0 impact parameter resolution_rphi [#mu m]; p_{t} [GeV/c]",nbins,xbins); \r
40 TGraphErrors *Resrphi_allpointRec =new TGraphErrors(26);\r
41 //Resrphi_allpointRec->SetMaximum(400); \r
42 //Resrphi_allpointRec->SetMinimum(0); \r
43 //Resrphi_allpoint->SetXTitle("p_{t} [GeV/c]");\r
44 //Resrphi_Rec->SetYTitle("d^{0} impact parameter resolution_rphi [#mu m]");\r
45 Resrphi_allpointRec->SetLineWidth(1);\r
46 Resrphi_allpointRec->SetMarkerStyle(20); \r
47 \r
48 TGraphErrors *Resrphi_allpointSkip1 =(TGraphErrors*)Resrphi_allpointRec->Clone("Resrphi_allpointSkip1");\r
49 \r
50 TGraphErrors *Resrphi_allpointSkipmean1 =(TGraphErrors*)Resrphi_allpointRec->Clone("Resrphi_allpointSkipmean1");\r
51 \r
52 \r
53 // fit gaussian\r
54 /*\r
55 TF1 *allpointRec1 = new TF1("allpointRec1","gaus",-10000.,10000.);\r
56 TF1 *allpointRec2 = new TF1("allpointRec2","gaus",-10000.,10000.);\r
57 TF1 *allpointRec3 = new TF1("allpointRec3","gaus",-10000.,10000.);\r
58 TF1 *allpointRec4 = new TF1("allpointRec4","gaus",-10000.,10000.);\r
59 TF1 *allpointRec5 = new TF1("allpointRec5","gaus",-10000.,10000.);\r
60 TF1 *allpointRec6 = new TF1("allpointRec6","gaus",-10000.,10000.);\r
61 */\r
62\r
63 \r
64 TF1 *f1 = new TF1("f1","gaus",-10000.,10000.);\r
65 TF1 *wn1 = new TF1("wn1","gaus",-10000.,10000.);\r
66 //fitting histogram\r
67 for (Int_t i=0;i<26;i++){\r
68 TString getstr;\r
69 Int_t jj=i+1;\r
70 Double_t j =3.; \r
71 getstr="d0allpointrphiSkip_"; \r
72 getstr += jj;\r
73 d0AllpointrphiSkip1_[i] = (TH1F*)d0allpointSkip1->FindObject(getstr);//->Clone(getstr+"d0allrphiRec");\r
74\r
75 getstr="d0pt_";\r
76 getstr += jj;\r
77 d0Pt_[i] = (TH1F*)d0pt->FindObject(getstr);\r
78\r
79 \r
80 //Set fitting function\r
81 //normalization\r
82 TF1 *h = new TF1("h","gaus",-10000.,10000.);\r
83 d0AllpointrphiSkip1_[i]->Fit("h","NR,Q");\r
84 Double_t d0rphirange_allpointSkip1 = h->GetParameter(2); \r
85 Double_t Entries = d0AllpointrphiSkip1_[i]->GetEntries();\r
86 Double_t Integral = d0AllpointrphiSkip1_[i]->Integral();\r
87 d0AllpointrphiSkip1_[i]->Scale(1./Integral);\r
88\r
89 //cp histogram\r
90 d0AllpointrphiSkipTail1_[i] = new TH1F(*d0AllpointrphiSkip1_[i]); // = new TH1F(*(TH1F*)d0AllpointrphiSkip1_[i]); \r
91 d0AllpointrphiSkipGaus1_[i] = new TH1F(*d0AllpointrphiSkip1_[i]);\r
92 Double_t wholerangeint=d0AllpointrphiSkip1_[i]->Integral();\r
93 Double_t cutleft1= -j*(d0rphirange_allpointSkip1);\r
94 Double_t cutright1 =j*d0rphirange_allpointSkip1;\r
95 d0AllpointrphiSkipTail1_[i]->Reset(0);\r
96 d0AllpointrphiSkipGaus1_[i]->Reset(0);\r
97 for (Int_t bin=1;bin<d0AllpointrphiSkip1_[i]->GetNbinsX();bin++){\r
98 Int_t bincontent = d0AllpointrphiSkip1_[i]->GetBinCenter(bin);\r
99 if(bincontent<cutleft1 || bincontent>cutright1) {\r
100 d0AllpointrphiSkipTail1_[i]->SetBinContent(bin,d0AllpointrphiSkip1_[i]->GetBinContent(bin));\r
101 d0AllpointrphiSkipTail1_[i]->SetBinError(bin,d0AllpointrphiSkip1_[i]->GetBinError(bin));\r
102 d0AllpointrphiSkipGaus1_[i]->SetBinContent(bin,0.); \r
103 //This sentence is very important,otherwise we will get the information the data is empty when we fit it .\r
104 }\r
105 else if(bincontent>=cutleft1 && bincontent<=cutright1){\r
106 d0AllpointrphiSkipTail1_[i]->SetBinContent(bin,0);\r
107 d0AllpointrphiSkipGaus1_[i]->SetBinContent(bin,d0AllpointrphiSkip1_[i]->GetBinContent(bin));\r
108 d0AllpointrphiSkipGaus1_[i]->SetBinError(bin,d0AllpointrphiSkip1_[i]->GetBinError(bin)); \r
109 }\r
110 //d0AllpointrphiSkipGaus_[i]->SetBinContent(bin,d0AllpointrphiSkip1_[i]->GetBinContent(bin));\r
111 //d0AllpointrphiSkipGaus_[i]->SetBinError(bin,d0AllpointrphiSkip1_[i]->GetBinError(bin)); \r
112 }\r
113 //return;\r
114 TF1 *hh;\r
115 hh =CreateFuncTail(d0AllpointrphiSkipTail1_[i],"hh");\r
116 hh->SetLineColor(4);\r
117 d0AllpointrphiSkipTail1_[i]->SetLineColor(5);\r
118 d0AllpointrphiSkipTail1_[i]->Fit("hh","NR,Q","",-2500,2500); \r
119 Double_t Sigmatail_allpointSkip1 = hh->GetParameter(2);\r
120 d0AllpointrphiSkipGaus1_[i]->Fit("h","NR,Q","",-2500,2500);\r
121 Double_t Sigmagaus_allpointSkip1 = h->GetParameter(2);\r
122 TF1 *allpointSkip1;\r
123 allpointSkip1=CreateFuncGaussTail(d0AllpointrphiSkip1_[i],"allpointSkip1",Sigmagaus_allpointSkip1,Sigmatail_allpointSkip1);\r
124 \r
125 \r
126\r
127 Double_t pt =d0Pt_[i]->GetMean();\r
128 Double_t RMS=d0Pt_[i]->GetRMS();\r
129 cout<<"i:"<<i<<endl; \r
130 cout<<"pt:"<<pt<<endl;\r
131\r
132 //fill the histogram \r
133 //Double_t ptpion=ptt*ptt/TMath::Sqrt(ptt*ptt+mass*mass); \r
134 //allpointSkip1->SetRange(-j*d0rphirange_allpointSkip1,j*d0rphirange_allpointSkip1);\r
135 d0AllpointrphiSkip1_[i]->Fit("allpointSkip1","NR,Q");\r
136 Resrphi_allpointSkip1->SetPoint(jj,pt,allpointSkip1->GetParameter(2));\r
137 Resrphi_allpointSkip1->SetPointError(jj,RMS,allpointSkip1->GetParError(2));\r
138 Resrphi_allpointSkipmean1->SetPoint(jj,pt,allpointSkip1->GetParameter(1));\r
139 Resrphi_allpointSkipmean1->SetPointError(jj,RMS,allpointSkip1->GetParError(1));\r
140 \r
141 }\r
142 \r
143 //********************************************************************************************** \r
144 //************************************ Plot figures ******************************************\r
145 //**********************************************************************************************\r
146\r
147 TCanvas *d0ResrphiSkip1 = new TCanvas("d0ResrphiSkip1","d0ResrphiSkip1",0,0,540,420); \r
148 d0ResrphiSkip1->cd();\r
149 gPad->SetLogx(); \r
150 TH2F *hFrame111 = new TH2F("hFrame","; p_{t} [GeV/c]; d0_resolution [#mum]",2,0.15,38,2,0,350); \r
151 hFrame111->Draw();\r
152 \r
153 Resrphi_allpointSkip1->SetMarkerStyle(20);\r
154 Resrphi_allpointSkip1->SetMarkerColor(2);\r
155 Resrphi_allpointSkip1->Draw("p"); \r
156 /*\r
157 TLegend *legend = new TLegend(0.28,0.75,0.65,0.88);\r
158 legend-> SetTextSize(0.025);\r
159 legend->SetBorderSize(0);\r
160 legend->SetFillColor(0);\r
161 legend->AddEntry(Resrphi_allpointSkip1,"d0ResrphiSkip_LHCbcd","pl");\r
162 //legend->AddEntry(Resrphi_allpointRec4,"d0ResrphiSkip_LHC11a10a","pl");\r
163 //legend->AddEntry(Resrphi_allpointRec2,"d0ResrphiSkip_LHC11a10a_m","pl");\r
164 legend->Draw();\r
165 */\r
166 d0ResrphiSkip1->Update();\r
167\r
168 TCanvas *d0MeanrphiSkip1 = new TCanvas("d0MeanrphiSkip1","d0MeanrphiSkip1",0,0,540,420); \r
169 d0MeanrphiSkip1->cd();\r
170 gPad->SetLogx(); \r
171 TH2F *hFrame1 = new TH2F("hFrame1","; p_{t} [GeV/c]; d0_mean [#mum]",2,0.15,38,2,-35,30);\r
172 hFrame1->Draw(); \r
173 Resrphi_allpointSkipmean1->SetMarkerStyle(20);\r
174 Resrphi_allpointSkipmean1->SetMarkerColor(2);\r
175 Resrphi_allpointSkipmean1->Draw("p"); \r
176 //Resrphi_allpointSkipmean2->SetMarkerStyle(20);\r
177 //Resrphi_allpointSkipmean2->SetMarkerColor(3);\r
178 //Resrphi_allpointSkipmean2->Draw("p,sames"); \r
179 \r
180 /*\r
181 TLegend *legend = new TLegend(0.2,0.7,0.7,0.88);\r
182 legend-> SetTextSize(0.025);\r
183 legend->SetBorderSize(0);\r
184 legend->SetFillColor(0);\r
185 legend->AddEntry(Resrphi_allpointSkipmean1,"d0MeanrphiSkip_LHC11a_partpointSkip_pass2_with_SDD","pl");\r
186 //legend->AddEntry(Resrphi_allpointSkipmean2,"d0MeanrphiSkip_LHC11a_partpointSkip_pass3_with_SDD","pl");\r
187 //legend->AddEntry(Resrphi_allpointSkipmean3,"d0MeanrphiSkip_LHC10d_partpointSkip_diamond","pl");\r
188 legend->Draw();\r
189 */\r
190 d0MeanrphiSkip1->Update();\r
191\r
192 return;\r
193\r
194}\r
195//---------------------------------------------------------------------------\r
196TF1 *CreateFuncTail(TH1F *hh,TString funcname,Double_t wholeRangeInt=-1.)\r
197{\r
198 /* TF1 *f3 = new TF1("f3","gaus",300,10000.);\r
199 cout<<"here:"<<2222222<<endl;\r
200 hh->Fit("f3","NR,Q");\r
201 Double_t value =f3->GetParameter(2);\r
202 cout<<"value:"<<value<<endl;*/\r
203 TF1 *tail=new TF1(funcname.Data(),"[0]*(1./(2.*[2])*TMath::Exp(-TMath::Abs(x-[1])/[2]))",-2500.,2500.);\r
204 //tail->SetParLimits(0,0.,250.);//Set the first parameter [0] range\r
205 tail->SetParLimits(1,-50,50.);//Set the second parameter [1] range\r
206 tail->SetParLimits(2,0,15000.);\r
207 Double_t binwidth=hh->GetBinWidth(10);\r
208 Double_t integral=hh->Integral();\r
209 if(wholeRangeInt!=-1.)tail->SetParLimits(0,(0.2)*wholeRangeInt*binwidth,(0.5)*wholeRangeInt*binwidth);\r
210 Double_t RMS1=TMath::Abs(hh->GetRMS());\r
211 Double_t firstvalue1=binwidth*integral;\r
212 //cout<<"binwidth:"<<binwidth<<endl;\r
213 //cout<<"integral:"<<integral<<endl;\r
214 //cout<<"RMS:"<<RMS<<endl;\r
215 //cout<<"firstvalue1:"<<firstvalue1<<endl;\r
216 tail->SetParameters(1.,0,100.);//Set the initial value of parameter\r
217 // tail->SetParameters(1,0,RMS1);\r
218 \r
219 return tail;\r
220\r
221} \r
222//----------------------------------------------------------------------------\r
223TF1 *CreateFuncGausTail(TH1F *h,TString funcname)\r
224{\r
225 TF1 *gaustail=new TF1(funcname.Data(),"[0]*([4]/(TMath::Sqrt(2.*TMath::Pi())*[2])*TMath::Exp(-1.*(x-[1])*(x-[1])/(2.*[2]*[2]))+(1.-[4])/(2.*[3])*TMath::Exp(-TMath::Abs(x-[1])/[3]))",-2500.,2500.);\r
226 gaustail->SetParLimits(0,0.,250000000.);//Set the first parameter [0] range\r
227 gaustail->SetParLimits(1,-20,20.);//Set the second parameter [1] range\r
228 gaustail->SetParLimits(2,0.,500.);\r
229 gaustail->SetParLimits(3,0.,500.);\r
230 gaustail->SetParLimits(4,0.,1.);\r
231 Double_t binwidth=h->GetBinWidth(10);\r
232 Double_t integral=h->Integral();\r
233 Double_t RMS1=TMath::Abs(h->GetRMS());\r
234 Double_t firstvalue=binwidth*integral;\r
235 //cout<<"binwidth:"<<binwidth<<endl;\r
236 //cout<<"integral:"<<integral<<endl;\r
237 //cout<<"RMS:"<<RMS<<endl;\r
238 cout<<"firstvalue:"<<firstvalue<<endl;\r
239 gaustail->SetParameters(firstvalue,0,RMS1,RMS1/20,0.82);//Set the initial value of parameter\r
240 //gaustail->FixParameter(0,firstvalue);\r
241 return gaustail;\r
242\r
243} \r
244//----------------------------------------------------------------------------\r
245TF1 *CreateFuncGaussTail(TH1F *h,TString funcname,Double_t parone,Double_t partwo)\r
246{\r
247 TF1 *f = new TF1("f","gaus",-10000.,10000.);\r
248 h->Fit("f","NR,Q");\r
249 Double_t value =f->GetParameter(2);\r
250 Double_t par1=1.1*parone;\r
251 Double_t par2=partwo*1.2;//par2*(1.-0.2),par2*(1.+0.2));\r
252 TF1 *gaustail=new TF1(funcname.Data(),"[0]*([4]/(TMath::Sqrt(2.*TMath::Pi())*[2])*TMath::Exp(-1.*(x-[1])*(x-[1])/(2.*[2]*[2]))+(1.-[4])/(2.*[3])*TMath::Exp(-TMath::Abs(x-[1])/[3]))",-2500.,2500.);\r
253 //gaustail->SetParLimits(0,0.,10.);//Set the first parameter [0] range\r
254 gaustail->SetParLimits(1,-50.,50.);//Set the second parameter [1] range\r
255 //gaustail->SetParLimits(2,0.,par1);\r
256 gaustail->SetParLimits(3,0,par2);\r
257 gaustail->SetParLimits(4,0.,1.);\r
258 Double_t binwidth=h->GetBinWidth(10);\r
259 Double_t integral=h->Integral();\r
260 Double_t RMS1=TMath::Abs(h->GetRMS());\r
261 Double_t firstvalue=binwidth*integral;\r
262 //cout<<"binwidth:"<<binwidth<<endl;\r
263 //cout<<"integral:"<<integral<<endl;\r
264 //cout<<"RMS:"<<RMS<<endl;\r
265 //cout<<"firstvalue:"<<firstvalue<<endl;\r
266 gaustail->SetParameters(1,0,value,partwo,0.82);//Set the initial value of parameter\r
267 \r
268 delete f;\r
269 return gaustail;\r
270\r
271} \r