PWGPP/macros/PlotImpParResVsPt.C

   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
   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
 196 TF1 *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
 223 TF1 *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
 245 TF1 *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