PWGDQ/dielectron/TakuAlberica/pair/get_spectrum_dir2.C

   1 void get_spectrum_dir2(int icent=1){
   2   gStyle->SetPalette(1);
   3   gStyle->SetOptStat(0);
   4   gStyle->SetOptFit(111111);
   5   //gStyle->SetOptFit(0);
   6   gStyle->SetOptTitle(0);
   7   //gStyle->SetFillColor(10);
   8   gStyle->SetCanvasColor(10);
   9   gStyle->SetFrameBorderMode(0);
  10   gStyle->SetFrameFillColor(0);
  11   gStyle->SetCanvasColor(0);
  12   gStyle->SetPadBorderSize(0);
  13   gStyle->SetCanvasBorderSize(0);
  14   //gStyle->SetPadLeftMargin(0.15);
  15   gStyle->SetPadLeftMargin(0.125);
  16   gStyle->SetPadBottomMargin(0.125);
  17   gStyle->SetPadTopMargin(0.1);
  18   gStyle->SetTitleYOffset(1.3);
  19   //gStyle->SetPadLeftMargin(0.1);
  20   gStyle->SetTitleW(0.7);
  21   gStyle->SetTitleH(0.1);
  22   cout<<"physics is always fun! "<<endl;
  23
  24   const int ncentbin=9;
  25   //  int cent_low[ncentbin] ={  0,  0, 20,  40, 40,  0,  0,  0,  30};
  26   //  int cent_high[ncentbin]={100, 10, 40, 100, 80, 80, 40, 30, 100};
  27
  28   //  int cent_low_bin[ncentbin]={0, 0, 2, 4, 4, 0, 0, 0, 3 };
  29   //  int cent_high_bin[ncentbin]={10, 1, 4, 10, 8, 8, 4, 3, 10};
  30
  31   int cent_low[ncentbin] ={ 0, 10, 20, 30};
  32   int cent_high[ncentbin]={10, 20, 30, 50};
  33
  34   int cent_low_bin[ncentbin] ={0, 1, 2, 3};
  35   int cent_high_bin[ncentbin]={1, 2, 3, 5};
  36
  37
  38   TH2D *hmasspt_inc[11][7];
  39   TH1D *hmasspt[4][5];
  40   TH1D *hmasspt_ls[4][5];
  41   TH1D *hdiv[5];
  42
  43   char outname[100];
  44   //sprintf(outname,"spec_histo_loose_eid_pt04_2_cent%d_no_gst.root",icent);
  45   sprintf(outname,"spec_histo_pt04_2_cent%d_no_gst.root",icent);
  46   //sprintf(outname,"output/spec_histo_loose_eid_2_pt04_cent%d_6.root",icent);
  47   //sprintf(outname,"spec_histo_loose_eid_pt04_cent%d_conv2.root",icent);
  48   char name[100];
  49   //  TFile *fin = new TFile("result_histo_nmix40_pt04.root");
  50   // TFile *fin = new TFile("result_histo_loose_eid_nmix10_conv2_pt04.root");
  51   //TFile *fin = new TFile("output/result_histo_loose_eid_nmix40_rev1_6.root");
  52   //TFile *fin = new TFile("result_histo_loose_eid_nmix40_pt04_no_gst.root");
  53   TFile *fin = new TFile("result_trig1_cut1.root");
  54   //TFile *fin = new TFile("result_histo_nmix20_pt04.root");
  55   for(int i=0;i<11;i++){
  56     for(int j=0;j<7;j++){
  57       sprintf(name,"hmasspt_weight_cent%d_pair%d",i,j);
  58       //sprintf(name,"hmasspt_cent%d_pair%d",i,j);
  59       hmasspt_inc[i][j] = (TH2D*)fin->Get(name);
  60       //cout<<i<<" "<<j<<" "<<hmasspt_inc[i][j]->GetEntries()<<endl;
  61     }
  62   }
  63
  64   TH1D *hpt = (TH1D*)hmasspt_inc[0][0]->ProjectionY("hpt");
  65
  66   TH1F *hstat = (TH1F*)fin->Get("hCentrality");
  67   float nevents = hstat->Integral(cent_low[icent], cent_high[icent]-1);
  68
  69   cout<<nevents<<endl;
  70
  71   const int nptbin=4;
  72   float pt_low[nptbin] =  {0.5, 1.0, 1.5, 2.0};
  73   float pt_high[nptbin] = {1.0, 1.5, 2.0, 5.0};
  74   double mass_r_ptbin[nptbin+2]={0, 0.5, 1.0, 1.5, 2.0, 4.0};
  75
  76   for(int i=0;i<4;i++){
  77     for(int j=0;j<nptbin;j++){
  78       sprintf(name,"hmass2_w_pair%d_pt%d", i, j);
  79       //sprintf(name,"hmass2_pair%d_pt%d", i, j);
  80       hmasspt[i][j] = new TH1D(name, name, 500, 0, 5);
  81       hmasspt[i][j]->Sumw2();
  82
  83       sprintf(name,"hmass2_ls_w_pair%d_pt%d", i, j);
  84       //sprintf(name,"hmass2_pair%d_pt%d", i, j);
  85       hmasspt_ls[i][j] = new TH1D(name, name, 500, 0, 5);
  86       hmasspt_ls[i][j]->Sumw2();
  87     }
  88   }
  89
  90   /// first []
  91   /// 0-->unlike
  92   /// 1-->like
  93   /// 2-->mixed unlike
  94   /// 3-->mixed like
  95
  96   /// second [] --> pt
  97   /// third [] --> ee or pp
  98   float npairs[4][nptbin][2];
  99
 100   for(int j=0;j<nptbin;j++){
 101     int binl = hpt->FindBin(pt_low[j]);
 102     int binh = hpt->FindBin(pt_high[j]);
 103
 104     npairs[0][j][0]=0; npairs[0][j][1]=0;
 105     npairs[1][j][0]=0; npairs[1][j][1]=0;
 106     npairs[2][j][0]=0; npairs[2][j][1]=0;
 107     npairs[3][j][0]=0; npairs[3][j][1]=0;
 108
 109
 110
 111     for(int i=cent_low_bin[icent]; i<cent_high_bin[icent]; i++){
 112
 113       ///////////////////////////////////////////
 114       ///// pair0
 115       ////////////////////////////////////////////
 116       sprintf(name,"hmass_%d_%d_0", i, j);
 117       hmasspt[0][j]->Add((TH1D*)hmasspt_inc[i][0]->ProjectionX(name,binl, binh-1));
 118       //      npairs[0][j][0] += hmasspt_inc[i][0]->ProjectionX(name,binl, binh-1)->Integral();
 119       npairs[0][j][0] += hmasspt_inc[i][0]->Integral();
 120
 121       /////////////////////////////////////////
 122       //// pair2 (like sign)
 123       ////////////////////////////////////////
 124       sprintf(name,"hmass_%d_%d_1", i, j);
 125       hmasspt[2][j]->Add((TH1D*)hmasspt_inc[i][1]->ProjectionX(name,binl, binh-1));
 126       hmasspt_ls[0][j]->Add((TH1D*)hmasspt_inc[i][1]->ProjectionX(name,binl, binh-1));
 127       //npairs[2][j][0] += hmasspt_inc[i][1]->ProjectionX(name,binl, binh-1)->Integral();
 128       npairs[2][j][0] += hmasspt_inc[i][1]->Integral();
 129
 130       sprintf(name,"hmass_%d_%d_2", i, j);
 131       hmasspt[2][j]->Add((TH1D*)hmasspt_inc[i][2]->ProjectionX(name,binl, binh-1));
 132       //npairs[2][j][1] += hmasspt_inc[i][2]->ProjectionX(name,binl, binh-1)->Integral();
 133       hmasspt_ls[1][j]->Add((TH1D*)hmasspt_inc[i][2]->ProjectionX(name,binl, binh-1));
 134       npairs[2][j][1] += hmasspt_inc[i][2]->Integral();
 135
 136       /////////////////////////////////////////
 137       //// pair1 (mixed unlike sign)
 138       /////////////////////////////////////////
 139       sprintf(name,"hmass_%d_%d_1", i, j);
 140       hmasspt[1][j]->Add((TH1D*)hmasspt_inc[i][3]->ProjectionX(name,binl, binh-1));
 141       //npairs[1][j][0] += hmasspt_inc[i][3]->ProjectionX(name,binl, binh-1)->Integral();
 142       npairs[1][j][0] += hmasspt_inc[i][3]->Integral();
 143
 144       sprintf(name,"hmass_%d_%d_2", i, j);
 145       hmasspt[1][j]->Add((TH1D*)hmasspt_inc[i][4]->ProjectionX(name,binl, binh-1));
 146       //      npairs[1][j][0] += hmasspt_inc[i][4]->ProjectionX(name,binl, binh-1)->Integral();
 147       npairs[1][j][0] += hmasspt_inc[i][4]->Integral();
 148
 149       /////////////////////////////////////////
 150       //// pair3 (mixed like sign)
 151       ////////////////////////////////////////
 152       sprintf(name,"hmass_%d_%d_1", i, j);
 153       //      hmasspt[3][j]->Add((TH1D*)hmasspt_inc[i][5]->ProjectionX(name,binl, binh-1));
 154       hmasspt_ls[2][j]->Add((TH1D*)hmasspt_inc[i][5]->ProjectionX(name,binl, binh-1));
 155       //npairs[3][j][0] += hmasspt_inc[i][5]->ProjectionX(name,binl, binh-1)->Integral();
 156       npairs[3][j][0] += hmasspt_inc[i][5]->Integral();
 157
 158       sprintf(name,"hmass_%d_%d_2", i, j);
 159       //hmasspt[3][j]->Add((TH1D*)hmasspt_inc[i][6]->ProjectionX(name,binl, binh-1));
 160       hmasspt_ls[3][j]->Add((TH1D*)hmasspt_inc[i][6]->ProjectionX(name,binl, binh-1));
 161       //npairs[3][j][1] += hmasspt_inc[i][6]->ProjectionX(name,binl, binh-1)->Integral();
 162       npairs[3][j][1] += hmasspt_inc[i][6]->Integral();
 163
 164     }
 165   }
 166
 167   /// for mixed like sign pairs
 168   /// entry should be same as like sign pairs
 169   /// re-add the hist for mixed like sign pairs using weights
 170   for(int j=0;j<nptbin;j++){
 171     float frac1 = (npairs[2][j][0]*(npairs[3][j][0]+npairs[3][j][1]))/(npairs[3][j][0]*(npairs[2][j][0]+npairs[2][j][1]));
 172     float frac2 = (npairs[2][j][1]*(npairs[3][j][0]+npairs[3][j][1]))/(npairs[3][j][1]*(npairs[2][j][0]+npairs[2][j][1]));
 173     hmasspt[3][j]->Add( hmasspt_ls[2][j], frac1); ///mixed like sign
 174     hmasspt[3][j]->Add( hmasspt_ls[3][j], frac2); // mixed like sign
 175     //hmasspt[3][j]->Add( hmasspt_ls[2][j], 1);
 176     //hmasspt[3][j]->Add( hmasspt_ls[3][j], 1);
 177   }
 178
 179
 180   for(int i=0;i<4;i++){
 181     for(int j=0;j<nptbin;j++){
 182       hmasspt[i][j]->Scale(1.0/nevents);
 183       hmasspt_ls[i][j]->Scale(1.0/nevents);
 184       cout<<hmasspt[i][j]->Integral()<<endl;
 185       cout<<npairs[i][j][0]<<" "<<npairs[i][j][1]<<endl;
 186     }
 187   }
 188
 189   TCanvas *c0 = new TCanvas("c0","c0",1500,500);
 190   c0->Divide(3,1);
 191   for(int iptbin=1;iptbin<nptbin;iptbin++){
 192     c0->cd(iptbin);
 193     gPad->SetGrid(1);
 194     gPad->SetLogy();
 195     hmasspt[0][iptbin]->SetXTitle("M_{ee}[GeV/c^{2}]");
 196     hmasspt[0][iptbin]->SetYTitle("N_{ee}/N_{evt}");
 197     hmasspt[0][iptbin]->Draw();
 198   }
 199
 200
 201   TCanvas *c1 = new TCanvas("c1","c1",1500,1000);
 202   c1->Divide(3,3);
 203
 204   TF1 *f1 = new TF1("f1","pol0",0.05,2);
 205
 206   double r[5];
 207
 208   for(int iptbin=1;iptbin<nptbin;iptbin++){
 209     sprintf(name,"hdiv_pt%d",iptbin);
 210     hdiv[iptbin] = new TH1D(name, name, 500,0,5);
 211     hdiv[iptbin]->SetXTitle("M_{ee}[GeV/c^{2}]");
 212     hdiv[iptbin]->SetYTitle("mixed unlike/mixed like");
 213     hdiv[iptbin]->Sumw2();
 214     //hmasspt[3][iptbin]->Scale(npairs[1][iptbin][0]/(npairs[3][iptbin][0]+npairs[3][iptbin][1]));
 215     hdiv[iptbin]->Divide(hmasspt[1][iptbin], hmasspt[3][iptbin]); ///mixed unlike/mixed like
 216     /*
 217     for(int ii=0;ii<hdiv[iptbin]->GetNbinsX(); ii++){
 218       double vent1 = hmasspt[1][iptbin]->GetBinContent(ii+1);
 219       double e_vent1 = hmasspt[1][iptbin]->GetBinError(ii+1);
 220
 221
 222       double vent2 = hmasspt_ls[2][iptbin]->GetBinContent(ii+1);
 223       double e_vent2 = hmasspt_ls[2][iptbin]->GetBinError(ii+1);
 224       double vent3 = hmasspt_ls[3][iptbin]->GetBinContent(ii+1);
 225       double e_vent3 = hmasspt_ls[3][iptbin]->GetBinError(ii+1);
 226
 227
 228
 229       double val1 = vent1;
 230       double val2 = 2*sqrt(vent2*vent3);
 231
 232       if(val2>0){
 233         hdiv[iptbin]->SetBinContent(ii+1, val1/val2);
 234         double e_val1 = e_vent1;
 235         double e_val2 = val2/2*sqrt(pow(e_vent2/vent2,2)+pow(e_vent3/vent3,2));
 236         hdiv[iptbin]->SetBinError(ii+1, val1/val2*sqrt(pow(e_val1/val1,2)+pow(e_val2/val2,2)));
 237       }
 238     }
 239     */
 240
 241     c1->cd(iptbin);
 242     gPad->SetGrid(1);
 243     hmasspt[1][iptbin]->SetXTitle("M_{ee}[GeV/c^{2}]");
 244     hmasspt[0][iptbin]->SetYTitle("N_{ee}/N_{evt}");
 245     hmasspt[1][iptbin]->Draw();
 246     hmasspt[3][iptbin]->SetLineColor(2);
 247     hmasspt[3][iptbin]->Draw("same");
 248
 249     c1->cd(iptbin+3);
 250     gPad->SetGrid(1);
 251     //hdiv[iptbin]->Rebin(4); hdiv[iptbin]->Scale(1.0/4);
 252
 253     hdiv[iptbin]->SetAxisRange(0,3);
 254     hdiv[iptbin]->SetMinimum(0);
 255     hdiv[iptbin]->SetMaximum(2);
 256     hdiv[iptbin]->Draw();
 257     ///hdiv[iptbin]->Fit(f1,"R","");
 258
 259     //r[iptbin] = f1->GetParameter(0);
 260     //int bin1 =  hmasspt[2][iptbin]->FindBin(0.1);
 261     //int bin2 =  hmasspt[2][iptbin]->FindBin(1.0);
 262     //r[iptbin] = hmasspt[2][iptbin]->Integral(bin1, bin2)/hmasspt[3][iptbin]->Integral(bin1, bin2);
 263     //cout<<" normalization  factor "<<r[iptbin]<<endl;
 264
 265     /*
 266     for(int k=0;k<hmasspt[1][iptbin]->GetNbinsX();k++){
 267       float w = f1->Eval(hmasspt[1][iptbin]->GetBinCenter(k+1));
 268       float rr = hmasspt[1][iptbin]->GetBinContent(k+1);
 269       float err = hmasspt[1][iptbin]->GetBinError(k+1);
 270       hmasspt[1][iptbin]->SetBinContent(k+1, rr*w);
 271       hmasspt[1][iptbin]->SetBinError(k+1, err*w);
 272     }
 273     */
 274     //hmasspt[1][iptbin]->Scale(r[iptbin]);
 275     hmasspt[2][iptbin]->Multiply(hdiv[iptbin]);
 276     //c/out<<npairs[1][iptbin][0]/(2*sqrt(npairs[3][iptbin][0]*npairs[3][iptbin][1]))<<endl;
 277     //cout<<" aaa "<<(2*sqrt(npairs[2][iptbin][0]*npairs[2][iptbin][1]))<<" "<<npairs[2][iptbin][0]+npairs[2][iptbin][1]<<" "<<hmasspt[2][iptbin]->Integral()*nevents<<endl;
 278     //hmasspt[2][iptbin]->Scale((2*sqrt(npairs[2][iptbin][0]*npairs[2][iptbin][1])/(npairs[2][iptbin][0]+npairs[2][iptbin][1])));
 279     //hmasspt[2][iptbin]->Scale(npairs[0][iptbin][0]/(2*sqrt(npairs[2][iptbin][0]*npairs[2][iptbin][1])));
 280
 281
 282     for(int ii=0;ii<hdiv[iptbin]->GetNbinsX(); ii++){
 283       double d_val1 = 2*sqrt(hmasspt_ls[0][iptbin]->GetBinContent(ii+1)*hmasspt_ls[1][iptbin]->GetBinContent(ii+1));
 284       if(d_val1>0){
 285         double d_eval1 = d_val1/2*sqrt(pow(hmasspt_ls[0][iptbin]->GetBinError(ii+1)/hmasspt_ls[0][iptbin]->GetBinContent(ii+1),2)+
 286                                         pow(hmasspt_ls[1][iptbin]->GetBinError(ii+1)/hmasspt_ls[1][iptbin]->GetBinContent(ii+1),2));
 287         //      double accep = hdiv[iptbin]->GetBinContent(ii+1)-hdiv[iptbin]->GetBinError(ii+1);
 288         double accep = hdiv[iptbin]->GetBinContent(ii+1);
 289         hmasspt[2][iptbin]->SetBinContent(ii+1, d_val1*accep);
 290         hmasspt[2][iptbin]->SetBinError(ii+1, d_eval1*accep);
 291       }
 292     }
 293
 294     /*
 295     double norm1=0;
 296     double norm2=0;
 297     for(int ii=0;ii<hmasspt[2][iptbin]->GetNbinsX(); ii++){
 298       double d_bg1 = hmasspt[1][iptbin]->GetBinContent(ii+1);
 299       double d_bg_pp = hmasspt_ls[2][iptbin]->GetBinContent(ii+1);
 300       double d_bg_mm = hmasspt_ls[3][iptbin]->GetBinContent(ii+1);
 301       double d_fg_pp = hmasspt_ls[0][iptbin]->GetBinContent(ii+1);
 302       double d_fg_mm = hmasspt_ls[1][iptbin]->GetBinContent(ii+1);
 303       double d_e_fg_pp = hmasspt_ls[0][iptbin]->GetBinError(ii+1);
 304       double d_e_fg_mm = hmasspt_ls[1][iptbin]->GetBinError(ii+1);
 305
 306       if(d_fg_pp>0 && d_fg_mm>0 && d_bg_pp>0 && d_bg_mm>0){
 307         //double d_val = 2*sqrt(d_fg_pp*d_fg_mm)*d_bg1/(2*sqrt(d_bg_pp*d_bg_mm));
 308         //double d_e_val = sqrt(d_fg_pp*d_fg_mm)*sqrt(pow(d_e_fg_pp/d_fg_pp,2)+pow(d_e_fg_mm/d_fg_mm,2))*d_bg1/(2*sqrt(d_bg_pp*d_bg_mm));
 309         double d_val = 2*(d_fg_pp*(0.5*d_bg1/d_bg_pp)*sqrt(npairs[2][iptbin][0])+d_fg_mm*(0.5*d_bg1/d_bg_mm)*sqrt(npairs[2][iptbin][1]))/(sqrt(npairs[2][iptbin][0])+sqrt(npairs[2][iptbin][1]));
 310         double d_e_val = d_e_fg_pp*(0.5*d_bg1/d_bg_pp)+d_e_fg_mm*(0.5*d_bg1/d_bg_mm);
 311
 312         hmasspt[2][iptbin]->SetBinContent(ii+1, d_val);
 313         hmasspt[2][iptbin]->SetBinError(ii+1, d_e_val);
 314
 315       }
 316     }
 317     */
 318     //hmasspt[2][iptbin]->Scale(hmasspt[0][iptbin]->Integral()/(2*sqrt(norm1*norm2)));
 319     cout<<" entry ="<<hmasspt[0][iptbin]->Integral()<<" "<<hmasspt[2][iptbin]->Integral()<<" "<<hmasspt_ls[0][iptbin]->Integral()<<" "<<hmasspt_ls[1][iptbin]->Integral()<<endl;
 320     c1->cd(iptbin+6);
 321     hmasspt[0][iptbin]->Draw();
 322     hmasspt[2][iptbin]->SetLineColor(2);
 323     hmasspt[2][iptbin]->Draw("same");
 324   }
 325
 326
 327   const int nmassbin=8;
 328   float mass_low[nmassbin] ={0,    0.05, 0.09, 0.14, 0.20, 0.30, 0.1, 0.09};
 329   float mass_high[nmassbin]={0.03, 0.09, 0.14, 0.20, 0.30, 0.40, 0.3, 0.30};
 330   float nsig[5][nmassbin];
 331   float nbg[5][nmassbin];
 332   float ensig[5][nmassbin];
 333   float enbg[5][nmassbin];
 334
 335   for(int ii=0;ii<nptbin; ii++){
 336     for(int j=0;j<nmassbin; j++){
 337       int binlow = hmasspt[0][ii]->FindBin(mass_low[j]);
 338       int binhigh = hmasspt[0][ii]->FindBin(mass_high[j]);
 339       nsig[ii][j] = hmasspt[0][ii]->Integral(binlow, binhigh-1);
 340       nbg[ii][j] = hmasspt[2][ii]->Integral(binlow, binhigh-1);
 341       ensig[ii][j]=0;
 342       enbg[ii][j]=0;
 343       for(int ib=binlow; ib<binhigh;ib++){
 344         ensig[ii][j] += pow(hmasspt[0][ii]->GetBinError(ib),2);
 345         enbg[ii][j] += pow(hmasspt[2][ii]->GetBinError(ib),2);
 346       }
 347
 348       ensig[ii][j] = sqrt(ensig[ii][j]);
 349       enbg[ii][j] = sqrt(enbg[ii][j]);
 350
 351
 352       float r1 = nsig[ii][j]-nbg[ii][j];
 353       float r0 =  nsig[ii][0]-nbg[ii][0];
 354       float er1 = sqrt(pow(ensig[ii][j],2)+pow(enbg[ii][j],2));
 355       float er0 = sqrt(pow(ensig[ii][0],2)+pow(enbg[ii][0],2));
 356
 357
 358       if(r0>0 && r1>0){
 359         cout<<" pt : "<<ii<<" mass "<<mass_low[j]<<"-"<<mass_high[j]<<" --> Nsig = "<<nsig[ii][j]<<", Nbg = "<<nbg[ii][j]
 360             <<" : Nsig-Nbg "<<nsig[ii][j]-nbg[ii][j]<<" +- "<<sqrt(pow(ensig[ii][j],2)+pow(enbg[ii][j],2))<<" ("<<ensig[ii][j]<<" "<<enbg[ii][j]<<")"
 361             <<" : Ratio = "<<r1/r0<<" +- "<<r1/r0*sqrt(pow(er1/r1,2)+pow(er0/r0,2))
 362             <<endl;
 363       }else{
 364         cout<<" pt : "<<ii<<" mass "<<mass_low[j]<<"-"<<mass_high[j]<<" --> Nsig = "<<nsig[ii][j]<<", Nbg = "<<nbg[ii][j]
 365             <<" : Nsig-Nbg "<<nsig[ii][j]-nbg[ii][j]<<" +- "<<sqrt(pow(ensig[ii][j],2)+pow(enbg[ii][j],2))<<" ("<<ensig[ii][j]<<" "<<enbg[ii][j]<<")"
 366             <<endl;
 367       }
 368     }
 369   }
 370
 371
 372   /// for histogram
 373   TH1D *hsub[5];
 374   const int nbins = 55;
 375   double x[nbins+1]={0, 0.01, 0.02, 0.03, 0.05, 0.09, 0.14, 0.20, 0.30, 0.40,
 376                      0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4,
 377                      1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4,
 378                      2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4,
 379                      3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4,
 380                      4.5, 4.6, 4.7, 4.8, 4.9, 5.0};
 381
 382 //  const int nbins = 53;
 383 //  double x[nbins+1]={0, 0.01, 0.02, 0.03, 0.05, 0.1, 0.20, 0.40,
 384 //                   0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4,
 385 //                   1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4,
 386 //                   2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4,
 387 //                   3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4,
 388 //                   4.5, 4.6, 4.7, 4.8, 4.9, 5.0};
 389
 390   int narrows[nptbin];
 391   TArrow *Arrow[nptbin][10];
 392
 393   for(int i=0; i<nptbin; i++){
 394     narrows[i]=0;
 395
 396
 397     sprintf(name,"hsub_mass_pt%d", i);
 398     hsub[i] = new TH1D(name, name, nbins, x);
 399
 400     for(int j=0; j<nbins; j++){
 401       int binl = hmasspt[0][i]->FindBin(x[j]);
 402       int binh = hmasspt[0][i]->FindBin(x[j+1]);
 403       float dw = x[j+1]-x[j];
 404
 405       float ent = 0;
 406       float e_ent = 0;
 407       float sig = 0;
 408       float bg = 0;
 409       for(int k=binl; k<binh; k++){
 410         sig += hmasspt[0][i]->GetBinContent(k);
 411         bg += hmasspt[2][i]->GetBinContent(k);
 412         ent += (hmasspt[0][i]->GetBinContent(k)-hmasspt[2][i]->GetBinContent(k));
 413         e_ent += (pow(hmasspt[0][i]->GetBinError(k),2)+
 414                   pow(hmasspt[2][i]->GetBinError(k),2));
 415       }
 416       e_ent = sqrt(e_ent);
 417       hsub[i]->SetBinContent(j+1, ent/dw);
 418       hsub[i]->SetBinError(j+1, e_ent/dw);
 419
 420       if(sig<bg && sig>0 && x[j]<0.5){ // calculate 90% confidence level
 421         double conf90=0;
 422         sig = sig*nevents;
 423         bg = bg*nevents;
 424         TH1F *htmp = new TH1F("htmp","htmp",1000,-2*bg, 2*bg);
 425         for(int kk=0;kk<100000;kk++){
 426           double a1 = gRandom->PoissonD(sig);
 427           double a2 = gRandom->PoissonD(bg);
 428           htmp->Fill(a1-a2);
 429         }
 430         int bin0 = htmp->FindBin(0);
 431         float ent0 = htmp->Integral(bin0, 1000);
 432         htmp->Scale(1.0/ent0);
 433         //htmp->Scale(1.0/100000);
 434         int bincent=0;
 435         for(int ib=bin0;ib<1000;ib++){
 436           if(htmp->Integral(bin0, ib)>0.9){
 437             bincent = ib;
 438             break;
 439           }
 440         }
 441         conf90 = htmp->GetBinCenter(bincent);
 442         cout<<i<<" "<<hsub[i]->GetBinCenter(j+1)<<" : 90% conf: "<<conf90<<" from sig :"<<sig<<" bg:"<<bg<<endl;
 443         Arrow[i][narrows[i]] = new TArrow(hsub[i]->GetBinCenter(j+1), 1.0e-07, hsub[i]->GetBinCenter(j+1), conf90/(dw*nevents),0.01,"<|");
 444         narrows[i]++;
 445         htmp->Clear();
 446       }// calculate 90% confidence level
 447
 448     }
 449   }
 450
 451   TH1D *hsub_clone[5];
 452
 453   TCanvas *c3 = new TCanvas("c3","c3",1200,600);
 454   c3->Divide(4,1);
 455
 456   c3->cd(1);
 457   gPad->SetGrid(1);  gPad->SetLogy();
 458   TH2F *hw = new TH2F("hw","hw",10,0,0.5,10,1.0e-07, 10);
 459   hw->Draw();
 460   hw->SetXTitle("M_{ee} [GeV]");
 461   hw->SetYTitle("dN_{ee}/dM_{ee}");
 462   c3->cd(2);
 463   gPad->SetGrid(1);  gPad->SetLogy();
 464   hw->Draw();
 465   c3->cd(3);
 466   gPad->SetGrid(1);  gPad->SetLogy();
 467   hw->Draw();
 468   c3->cd(4);
 469   gPad->SetGrid(1);  gPad->SetLogy();
 470   hw->Draw();
 471
 472   for(int i=0; i<nptbin; i++){
 473     hsub_clone[i] = (TH1D*)hsub[i]->Clone();
 474     sprintf(name,"%s_cln",hsub[i]->GetName());
 475     hsub_clone[i]->SetName(name);
 476     hsub_clone[i]->SetLineColor(i+1);
 477     hsub_clone[i]->SetMarkerColor(i+1);
 478     hsub_clone[i]->SetMarkerStyle(20);
 479   }
 480
 481
 482   //  int bin = hsub_clone[1]->FindBin(0.03);
 483   //  float ent = hsub_clone[1]->Integral(0,bin);
 484   for(int i=0; i<nptbin; i++){
 485     //float ent1 = hsub_clone[i]->Integral(0,bin);
 486     //hsub_clone[i]->Scale(ent/ent1);
 487     c3->cd(i+1);
 488     hsub_clone[i]->Draw("same");
 489     for(int k=0;k<narrows[i];k++){
 490       Arrow[i][k]->SetLineColor(i+1);
 491       Arrow[i][k]->Draw();
 492     }
 493   }
 494
 495   TCanvas *c5 = new TCanvas("c5","c5",1200,600);
 496   c5->Divide(3,1);
 497
 498   c5->cd(1);
 499   gPad->SetGrid(1);  gPad->SetLogy();
 500   TH2F *hw2 = new TH2F("hw2","hw2",10,0,0.5,10,1.0e-07, 10);
 501   hw2->Draw();
 502   hw2->SetXTitle("M_{ee} [GeV]");
 503   hw2->SetYTitle("1/N_{evt}dN_{ee}/dM_{ee}");
 504   c5->cd(2);
 505   gPad->SetGrid(1);  gPad->SetLogy();
 506   hw2->Draw();
 507   c5->cd(3);
 508   gPad->SetGrid(1);  gPad->SetLogy();
 509   hw2->Draw();
 510
 511   TH1D *hmass_bin[2][nptbin];
 512   for(int i=1; i<nptbin; i++){
 513     c5->cd(i);
 514     hmass_bin[0][i] = (TH1D*)hmasspt[0][i]->Clone();
 515     hmass_bin[1][i] = (TH1D*)hmasspt[2][i]->Clone();
 516
 517
 518     float binw = hmass_bin[0][i]->GetBinWidth(10);
 519     hmass_bin[0][i]->Scale(1/binw);
 520     hmass_bin[1][i]->Scale(1/binw);
 521
 522     hmass_bin[1][i]->SetLineColor(2);
 523     hmass_bin[0][i]->Draw("same");
 524     hmass_bin[1][i]->Draw("same");
 525     hsub_clone[i]->Draw("same");
 526     for(int k=0;k<narrows[i];k++){
 527       Arrow[i][k]->SetLineColor(i+1);
 528       Arrow[i][k]->Draw();
 529     }
 530   }
 531
 532
 533
 534
 535
 536
 537   //  TLegend *leg = new TLegend(0.5,0.7,0.875,0.875);
 538   //  leg->SetFillColor(10);
 539   //  if(nptbin==4){
 540   //    leg->AddEntry(hsub_clone[0],"0.5<p_{T}<1.0 GeV","pl");
 541   //    leg->AddEntry(hsub_clone[1],"1.0<p_{T}<1.5 GeV","pl");
 542   //    leg->AddEntry(hsub_clone[2],"1.5<p_{T}<2.0 GeV","pl");
 543   //    leg->AddEntry(hsub_clone[3],"2.0<p_{T}","pl");
 544   //  }else{
 545   //  }
 546   //  leg->Draw("same");
 547
 548   /////draw the results
 549   TH1D *hmass_cocktail[4];
 550   TH1D *hmass_cocktail_comp[4][6];
 551   TFile *fin1 = new TFile("../../ana2/cocktail/cocktail_allpt_pt04_no_smearing.root");
 552   fin1->cd();
 553   hmass_cocktail[0] = (TH1D*)fin1->Get("hmass_all");
 554   hmass_cocktail[0]->SetName("hmass_all_pt0");
 555   hmass_cocktail_comp[0][0] =  (TH1D*)fin1->Get("hmass_3");
 556   hmass_cocktail_comp[0][0]->SetName("hmass_pt0_eta_prime");
 557   hmass_cocktail_comp[0][1] =  (TH1D*)fin1->Get("hmass_4");
 558   hmass_cocktail_comp[0][1]->SetName("hmass_pt0_rho");
 559   hmass_cocktail_comp[0][2] =  (TH1D*)fin1->Get("hmass_7");
 560   hmass_cocktail_comp[0][2]->SetName("hmass_pt0_pi0");
 561   hmass_cocktail_comp[0][3] =  (TH1D*)fin1->Get("hmass_8");
 562   hmass_cocktail_comp[0][3]->SetName("hmass_pt0_eta");
 563   hmass_cocktail_comp[0][4] =  (TH1D*)fin1->Get("hmass_9");
 564   hmass_cocktail_comp[0][4]->SetName("hmass_pt0_omega");
 565   hmass_cocktail_comp[0][5] =  (TH1D*)fin1->Get("hmass_10");
 566   hmass_cocktail_comp[0][5]->SetName("hmass_pt0_phi");
 567
 568
 569   //TFile *fin2 = new TFile("./cocktail//cocktail_pt10_15.root");
 570   TFile *fin2 = new TFile("../../ana2/cocktail//cocktail_pt04_pairpt_10_15_no_smearing.root");
 571   fin2->cd();
 572   hmass_cocktail[1] = (TH1D*)fin2->Get("hmass_all");
 573   hmass_cocktail[1]->SetName("hmass_all_pt1");
 574   hmass_cocktail_comp[1][0] =  (TH1D*)fin1->Get("hmass_3");
 575   hmass_cocktail_comp[1][0]->SetName("hmass_pt1_eta_prime");
 576   hmass_cocktail_comp[1][1] =  (TH1D*)fin1->Get("hmass_4");
 577   hmass_cocktail_comp[1][1]->SetName("hmass_pt1_rho");
 578   hmass_cocktail_comp[1][2] =  (TH1D*)fin1->Get("hmass_7");
 579   hmass_cocktail_comp[1][2]->SetName("hmass_pt1_pi0");
 580   hmass_cocktail_comp[1][3] =  (TH1D*)fin1->Get("hmass_8");
 581   hmass_cocktail_comp[1][3]->SetName("hmass_pt1_eta");
 582   hmass_cocktail_comp[1][4] =  (TH1D*)fin1->Get("hmass_9");
 583   hmass_cocktail_comp[1][4]->SetName("hmass_pt1_omega");
 584   hmass_cocktail_comp[1][5] =  (TH1D*)fin1->Get("hmass_10");
 585   hmass_cocktail_comp[1][5]->SetName("hmass_pt1_phi");
 586
 587
 588
 589   //TFile *fin3 = new TFile("./cocktail//cocktail_pt15_20.root");
 590   //TFile *fin3 = new TFile("./cocktail//cocktail_pt15_20_no_smearing.root");
 591   TFile *fin3 = new TFile("../../ana2/cocktail//cocktail_pt04_pairpt_15_20_no_smearing.root");
 592   fin3->cd();
 593   hmass_cocktail[2] = (TH1D*)fin3->Get("hmass_all");
 594   hmass_cocktail[2]->SetName("hmass_all_pt2");
 595   hmass_cocktail_comp[2][0] =  (TH1D*)fin1->Get("hmass_3");
 596   hmass_cocktail_comp[2][0]->SetName("hmass_pt2_eta_prime");
 597   hmass_cocktail_comp[2][1] =  (TH1D*)fin1->Get("hmass_4");
 598   hmass_cocktail_comp[2][1]->SetName("hmass_pt2_rho");
 599   hmass_cocktail_comp[2][2] =  (TH1D*)fin1->Get("hmass_7");
 600   hmass_cocktail_comp[2][2]->SetName("hmass_pt2_pi0");
 601   hmass_cocktail_comp[2][3] =  (TH1D*)fin1->Get("hmass_8");
 602   hmass_cocktail_comp[2][3]->SetName("hmass_pt2_eta");
 603   hmass_cocktail_comp[2][4] =  (TH1D*)fin1->Get("hmass_9");
 604   hmass_cocktail_comp[2][4]->SetName("hmass_pt2_omega");
 605   hmass_cocktail_comp[2][5] =  (TH1D*)fin1->Get("hmass_10");
 606   hmass_cocktail_comp[2][5]->SetName("hmass_pt2_phi");
 607
 608
 609
 610   //TFile *fin4 = new TFile("./cocktail//cocktail_pt20_no_smearing.root");
 611   //TFile *fin4 = new TFile("./cocktail//cocktail_pt20.root");
 612 TFile *fin4 = new TFile("../../ana2/cocktail//cocktail_pt04_pairpt_20_no_smearing.root");
 613   fin4->cd();
 614   hmass_cocktail[3] = (TH1D*)fin4->Get("hmass_all");
 615   hmass_cocktail[3]->SetName("hmass_all_pt3");
 616   hmass_cocktail_comp[3][0] =  (TH1D*)fin1->Get("hmass_3");
 617   hmass_cocktail_comp[3][0]->SetName("hmass_pt3_eta_prime");
 618   hmass_cocktail_comp[3][1] =  (TH1D*)fin1->Get("hmass_4");
 619   hmass_cocktail_comp[3][1]->SetName("hmass_pt3_rho");
 620   hmass_cocktail_comp[3][2] =  (TH1D*)fin1->Get("hmass_7");
 621   hmass_cocktail_comp[3][2]->SetName("hmass_pt3_pi0");
 622   hmass_cocktail_comp[3][3] =  (TH1D*)fin1->Get("hmass_8");
 623   hmass_cocktail_comp[3][3]->SetName("hmass_pt3_eta");
 624   hmass_cocktail_comp[3][4] =  (TH1D*)fin1->Get("hmass_9");
 625   hmass_cocktail_comp[3][4]->SetName("hmass_pt3_omega");
 626   hmass_cocktail_comp[3][5] =  (TH1D*)fin1->Get("hmass_10");
 627   hmass_cocktail_comp[3][5]->SetName("hmass_pt3_phi");
 628
 629
 630
 631   for(int i=0;i<4;i++){
 632     hmass_cocktail[i]->Rebin(4);
 633     float ent = 0;
 634     for(int j=0;j<hsub_clone[i]->GetNbinsX();j++){
 635       if(hsub_clone[i]->GetBinCenter(j+1)<0.03){
 636         ent += hsub_clone[i]->GetBinWidth(j+1)*hsub_clone[i]->GetBinContent(j+1);
 637       }
 638     }
 639     float ent2=0;
 640     for(int j=0;j<hmass_cocktail[i]->GetNbinsX();j++){
 641       if(hmass_cocktail[i]->GetBinCenter(j+1)<0.03){
 642         ent2 += hmass_cocktail[i]->GetBinWidth(j+1)*hmass_cocktail[i]->GetBinContent(j+1);
 643       }
 644     }
 645     hmass_cocktail[i]->Scale(ent/ent2);
 646     hmass_cocktail[i]->SetLineColor(6);
 647     hmass_cocktail[i]->SetLineWidth(2);
 648     c3->cd(i+1);
 649     hmass_cocktail[i]->Draw("same,l");
 650     if(i>=1){
 651       c5->cd(i);
 652       hmass_cocktail[i]->Draw("same,l");
 653     }
 654
 655     for(int j=0;j<6;j++){
 656       hmass_cocktail_comp[i][j]->Rebin(4);
 657       hmass_cocktail_comp[i][j]->Scale(ent/ent2*0.2);
 658       hmass_cocktail_comp[i][j]->SetLineColor(1);
 659       //hmass_cocktail_comp[i][j]->Draw("same");
 660     }
 661   }
 662
 663   float r_cocktail[nptbin];
 664   cout<<" ***************** "<<endl;
 665   cout<<" mass ratio ("<<mass_low[nmassbin-1]<<"-"<<mass_high[nmassbin-1]<<"/0-30) in cocktail"<<endl;
 666   TH1F *hmass_r_ck = new TH1F("hmass_r_cocktail","mass ratio cocktail", nptbin+1, mass_r_ptbin);
 667   for(int i=1;i<nptbin;i++){
 668     int bin1=hmass_cocktail[i]->FindBin(0);
 669     int bin2=hmass_cocktail[i]->FindBin(0.03);
 670     float ent1 = hmass_cocktail[i]->Integral(bin1, bin2-1);
 671     bin1=hmass_cocktail[i]->FindBin(mass_low[nmassbin-1]);
 672     bin2=hmass_cocktail[i]->FindBin(mass_high[nmassbin-1]);
 673     float ent2 = hmass_cocktail[i]->Integral(bin1, bin2-1);
 674     cout<<" pt : "<<i<<" "
 675         <<ent2/ent1
 676         <<endl;
 677     int bin = hmass_r_ck->FindBin((0.5*(pt_low[i]+pt_high[i])));
 678     hmass_r_ck->SetBinContent(bin, ent2/ent1);
 679     hmass_r_ck->SetBinError(bin,  ent2/ent1*0.05);
 680     r_cocktail[i] = ent2/ent1;
 681   }
 682
 683   float r_dir = 0.37; // this is from 90-300MeV
 684
 685   TH1F *hmass_r = new TH1F("hmass_r","mass ratio", nptbin+1, mass_r_ptbin);
 686   TH1F *hrgamma = new TH1F("hrgamma","r gmamma", nptbin+1, mass_r_ptbin);
 687   cout<<" ***************** "<<endl;
 688   cout<<" mass ratio (100-300/0-30)"<<endl;
 689   for(int i=1; i<nptbin; i++){
 690     cout<<" pt : "<<i<<" "
 691         <<nsig[i][nmassbin-1]-nbg[i][nmassbin-1]<<" +- "<<sqrt(pow(ensig[i][nmassbin-1],2)+pow(enbg[i][nmassbin-1],2))<<" / "
 692         <<(nsig[i][0]-nbg[i][0])<<" ratio = "
 693         <<(nsig[i][nmassbin-1]-nbg[i][nmassbin-1])/(nsig[i][0]-nbg[i][0])
 694         <<" +- "<<sqrt(pow(ensig[i][nmassbin-1],2)+pow(enbg[i][nmassbin-1],2))/(nsig[i][0]-nbg[i][0])
 695         <<endl;
 696
 697     int bin = hmass_r->FindBin((0.5*(pt_low[i]+pt_high[i])));
 698     float r_data = (nsig[i][nmassbin-1]-nbg[i][nmassbin-1])/(nsig[i][0]-nbg[i][0]);
 699     float er_data = sqrt(pow(ensig[i][nmassbin-1],2)+pow(enbg[i][nmassbin-1],2))/(nsig[i][0]-nbg[i][0]);
 700
 701     hmass_r->SetBinContent(bin, r_data);
 702     hmass_r->SetBinError(bin, er_data);
 703
 704     float r_gamma = (r_data-r_cocktail[i])/(r_dir - r_cocktail[i]);
 705     float e_r_gamma = (er_data)/(r_dir - r_cocktail[i]);
 706     hrgamma->SetBinContent(bin, r_gamma);
 707     hrgamma->SetBinError(bin, e_r_gamma);
 708
 709   }
 710
 711   hrgamma->SetXTitle("p_{T} [GeV]");
 712   hrgamma->SetYTitle("R_{#gamma}^{0-30}");
 713   hrgamma->SetMarkerColor(icent+1);
 714   hrgamma->SetLineColor(icent+1);
 715   hrgamma->SetMarkerStyle(20);
 716   hrgamma->SetMinimum(0);
 717   hrgamma->SetMaximum(1);
 718
 719   TCanvas *c4 = new TCanvas("c4","c4");
 720   gPad->SetGrid(1);
 721   hmass_r->SetXTitle("p_{T} [GeV]");
 722   hmass_r->SetYTitle("R_{data}^{90-300}");
 723   hmass_r->SetMarkerColor(icent+1);
 724   hmass_r->SetLineColor(icent+1);
 725   hmass_r->SetMarkerStyle(20);
 726   hmass_r->SetMinimum(0);
 727   hmass_r->SetMaximum(0.5);
 728
 729   hmass_r_ck->SetMarkerStyle(24);
 730   hmass_r->Draw(); hmass_r_ck->Draw("same");
 731
 732
 733
 734
 735
 736
 737   TFile *fout = new TFile(outname,"recreate");
 738   hw->Write();
 739   for(int i=0; i<nptbin; i++){
 740     hsub[i]->Write();
 741     hsub_clone[i]->Write();
 742   }
 743   for(int j=0;j<4;j++){
 744     hmass_cocktail[j]->Write();
 745     for(int k=0;k<6;k++){
 746       hmass_cocktail_comp[j][k]->Write();
 747     }
 748   }
 749   hmass_r->Write();
 750   hmass_r_ck->Write();
 751   hrgamma->Write();
 752
 753 }
 754
 755
 756
 757
 758