PWGPP/macros/PlotImpParResVsPt.C

   1 TF1 *CreateFuncGausTail(TH1F *h,TString funcname);
   2 TF1 *CreateFuncGaussTail(TH1F *h,TString funcname,Double_t parone,Double_t partwo);
   3
   4 void 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 //---------------------------------------------------------------------------
 203 TF1 *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 //----------------------------------------------------------------------------
 230 TF1 *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 //----------------------------------------------------------------------------
 252 TF1 *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 }