[u/mrichter/AliRoot.git] / PWGPP / macros / PlotImpParResVsPt.C

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