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