StRoot: StSpinPool/StWalgoB2009/macros/plWJacob.C Source File

00001 int blueCol=4;
00002 int yellCol=kYellow;
00003 int stage=3;
00004 float yMax=175; 
00005 float etMax=70;
00006 
00007 void plWJacob(int charge=0) {
00008   gStyle->SetPalette(1,0);
00009   gStyle->SetOptStat(1000000);
00010 
00011   TString iPath="/star/data05/scratch/stevens4/wAnalysisOut/apsXsec/data/";
00012   TString core="run9setABCD";
00013 
00014   TFile *fd = TFile::Open(iPath+core+".wana.hist.root"); 
00015  
00016   assert(fd->IsOpen());
00017   
00018   if(charge==0){
00019     //charge summed histos 
00020     hwE=(TH1F*)fd->Get("muclustPtBal"); assert(hwE);
00021     hnE=(TH1F*)fd->Get("muclustPtBalnoE"); assert(hnE);
00022     hQCD=(TH1F*)fd->Get("muclustPtBal_bckgrd"); assert(hQCD);
00023   }
00024   else if(charge==1){
00025     //charge summed histos 
00026     hwE=(TH1F*)fd->Get("pubclustPtBalP"); assert(hwE);
00027     hnE=(TH1F*)fd->Get("pubclustPtBalnoEP"); assert(hnE);
00028     hQCD=(TH1F*)fd->Get("pubclustPtBal_bckgrdP"); assert(hQCD);
00029   }
00030   else if(charge==-1){
00031     //charge summed histos 
00032     hwE=(TH1F*)fd->Get("pubclustPtBalN"); assert(hwE);
00033     hnE=(TH1F*)fd->Get("pubclustPtBalnoEN"); assert(hnE);
00034     hQCD=(TH1F*)fd->Get("pubclustPtBal_bckgrdN"); assert(hQCD);
00035   }
00036 
00037   // plot raw spectra ..........  
00038   c=new TCanvas();
00039   c->Divide(2,3);
00040   c->cd(2);
00041   hwE->Draw();  
00042   c->cd(4);
00043   hnE->Draw();  // hnE->SetAxisRange(0.,etMax);
00044 
00045   //.. compute Intermediate results
00046   TAxis* axX=hwE->GetXaxis();
00047   float x1=axX->GetXmin()+1.0;
00048   float x2=axX->GetXmax()+1.0;
00049   int   nbX=axX->GetNbins();
00050 
00051   // assumption: output histos will have 100 bins [1,101] GeV
00052 
00053   // those 3 have errors =sqrt(N)
00054   hA=new TH1F("hA","W candidates (input)",nbX,x1,x2);
00055   hAA=new TH1F("hAA","W candidates- 2nd Endcap",nbX,x1,x2);
00056 
00057   hB=new TH1F("hB","rejected by including Endcap",nbX,x1,x2);
00058   hC=new TH1F("hC","rejected by inverting  cut",nbX,x1,x2);
00059 
00060  // .......copy & subtract input spectra, subtract 2nd endcap
00061   float sumA=0,sumB=0,sumC=0; // over ET<20
00062   for(int i=1;i<=nbX-1;i++) {
00063     float y1=hwE->GetBinContent(i+1);
00064     float y2=hnE->GetBinContent(i+1);
00065     float y3=hQCD->GetBinContent(i+1);
00066     //  printf("%d , %.1f,  %.0f, %.0f, %.0f \n",i,hA->GetBinCenter(i),y1,y2,y3);
00067     // compute yields & comulative variance 
00068     float yA=y1;
00069     float yB=y2-y1; // Endcap background, no adding errors in quadrature
00070     float yC=y3; // rejected by away ET
00071     hA->SetBinContent(i,yA); // W candidates
00072     hB->SetBinContent(i,yB);  // endcapBCG
00073     hC->SetBinContent(i,yC);  // away BCG
00074   }
00075   
00076   hC->Sumw2();   hAA->Sumw2();
00077   // subtract 2nd endcap
00078   hAA->Add(hA,1.); hAA->Add(hB,-1.);
00079 
00080   c->cd(1); hA->Draw();//   hB->SetFillColor(8);
00081 
00082   c->cd(3); hB->Draw();   hB->SetFillColor(8);
00083 
00084 
00085   c->cd(5); hC->Draw("he");
00086   hC->SetLineColor(kRed);
00087   hC->SetLineWidth(2);
00088 
00089   hA->Rebin();   hAA->Rebin();  hB->Rebin();  hC->Rebin(); 
00090   nbX/=2;
00091   hA->SetLineWidth(2.);  hA->SetAxisRange(0.,etMax); 
00092   hB->SetAxisRange(0.,etMax);
00093   hC->SetAxisRange(0.,etMax);
00094 
00095   float sumA=0, sumC=0; // over ET =[15,21] GeV, delET=6 GeV
00096   for(int i=1;i<=nbX;i++) {// decides on backg norm. range
00097     float y1=hAA->GetBinContent(i);
00098     float y2=hC->GetBinContent(i);
00099     //    printf("%d %f %f    \n",i,y1,y2);
00100     if(i>=8 && i<=9) { sumA+=y1; sumC+=y2;}
00101   }
00102   float fact=sumA/sumC;
00103   printf("sumA=%.1f sumC=%.1f fac=%.3f\n",sumA, sumC, fact);
00104   hC->Scale(fact);
00105 
00106   //..........................................
00107   //..........................................
00108   // ..... "Preliminary" Jacobian peak
00109   //..........................................
00110   //..........................................
00111   TString tt="W2009-prelimB"; tt+=stage;
00112 
00113   c=new TCanvas(tt,tt,600,600);
00114   gStyle->SetOptStat(0); 
00115   c->SetFillColor(kWhite);
00116   c->SetLeftMargin(0.15);
00117   c->SetFrameBorderMode(0);
00118 
00119 
00120 
00121   // input W spectrum
00122   hA3=(TH1F*) hA->Clone();
00123   hA3->SetTitle("PT Balance Cut; EMC cluster #font[72]{E_{T}} (GeV)       ;Counts"); 
00124   hA3->GetYaxis()->SetTitleOffset(1.5);
00125   hA3->SetAxisRange(7,69); hA3->SetMinimum(-39); hA3->SetMaximum(yMax);
00126   hA3->Draw("h"); hA3->SetLineWidth(2);
00127   ln=new TLine(8,0,70,0); ln->Draw();
00128   printf("raw sum=%.0f\n",hA3->Integral());
00129   hA3->SetTitle(iPath+core);
00130 
00131   // pure Ws
00132   hDD=(TH1F*) hAA->Clone(); 
00133   hDD->Add(hC,-1.);  // subtracted scaled away ET 
00134   hDD->SetMarkerStyle(8); 
00135   hDD2=(TH1F*) hDD->Clone(); 
00136   hDD2->SetMarkerStyle(8); hDD2->SetMarkerSize(1);
00137   hDD->SetFillColor(yellCol);  //hDD2->SetLineWidth(2);
00138   
00139   if( stage>2) { hDD->Draw("same h");  hDD2->Draw("same h"); }
00140 
00141   // away backg+ Endcap backg
00142   hCB=(TH1F*) hC->Clone();
00143   hCB->Add(hB,1.);
00144   hCB->SetMarkerStyle(8); hCB->SetMarkerSize(1); hCB->SetMarkerColor(blueCol);
00145   hCB->SetLineColor(blueCol);  hCB->SetLineWidth(2);
00146   
00147   hCB->SetAxisRange(14.5,65);
00148   if(stage>1) hCB->Draw("p e h same");
00149 
00150   // draw inut again w/ dashed style
00151   hA3b=(TH1F*) hA3->Clone();
00152   hA3b->SetLineStyle(2);
00153   hA3b->Draw("same");   hA3b->SetAxisRange(15,18);
00154 
00155 
00156   lg=new TLegend(0.32,0.64,0.66,0.80); // top
00157   TDatime dt;lg->SetHeader(dt.AsString());
00158   //lg->SetHeader("     cluster |#eta|_{ }<_{ }1");
00159   if(charge==0) 
00160     lg->AddEntry(hA,"#font[52]{W} candidates (charge summed)","l");
00161   else if(charge==1)
00162     lg->AddEntry(hA,"#font[52]{W+} candidates ","l");
00163   else 
00164     lg->AddEntry(hA,"#font[52]{W-} candidates ","l");
00165   lg->AddEntry(hCB,"QCD backg. est.","lfep");
00166   lg->AddEntry(hDD,"Backg. subtr. #font[52]{W}s","lfep");
00167   lg->SetFillColor(kWhite);
00168   lg->SetLineColor(kWhite);
00169   lg->Draw();
00170 
00171   TLatex *lat1 = new TLatex(0.32,0.84,"STAR   #font[72]{#vec{p}_{ }+_{ }#vec{p}}   #sqrt{#font[72]{s}}_{ }=_{ }500 GeV");
00172   lat1->SetNDC();  lat1->Draw("same");
00173 
00174 
00175   tx=new TText(14,20,"software threshold"); tx->SetTextAngle(90); 
00176   tx->SetTextSize(0.03);//tx->SetTextColor(kMagenta);
00177   tx->Draw();
00178 
00179   //....... Joe's additional variables
00180   // Calcate the total relative error squared for the cross section
00181   // which is roughly (S+2B)/S^2 
00182   // and for the AL which is roughly (S+B?)/S^2
00183 
00184   float tot_xsec_err_num = 0.;
00185   float tot_xsec_err_den = 0.;
00186   float tot_al_err_num = 0.;
00187   float tot_al_err_den = 0.;
00188   for (int i=11; i<=hA->GetNbinsX(); i++) {// defines ET range for x-section
00189     tot_xsec_err_num += hA->GetBinError(i)*hA->GetBinError(i);
00190     tot_xsec_err_num += hB->GetBinError(i)*hB->GetBinError(i);
00191     tot_xsec_err_num += hC->GetBinError(i)*hC->GetBinError(i);
00192     tot_xsec_err_den += hDD->GetBinContent(i);
00193     // cout << i << " " << hA->GetBinContent(i) << endl;
00194   }
00195 
00196   for (int i=13; i<=hA->GetNbinsX(); i++) { // defines ET range for AL
00197     tot_al_err_num += hA->GetBinError(i)*hA->GetBinError(i);
00198     tot_al_err_den += hDD->GetBinContent(i); 
00199   }
00200   
00201   float tot_xsec_rel_err = sqrt(tot_xsec_err_num)/tot_xsec_err_den;
00202   float tot_al_err = sqrt(tot_al_err_num)/tot_al_err_den;
00203   printf(" Joe measures: xsec_rel_err=%.3f  AL_err=%.3f\n", tot_xsec_rel_err, tot_al_err);
00204 
00205   char text[100];
00206   sprintf(text,"Error estim:  xSec_rel=%.4f    AL=%.3f", tot_xsec_rel_err, tot_al_err);
00207   tx=new TText(21,-18,text); tx->Draw(); tx->SetTextSize(0.03);
00208 
00209   //...... end Joe' calculation
00210 
00211   return; // not print
00212   for(int i=1;i<=nbX;i++) {
00213     printf("%d , %.1f,",i,hA->GetBinCenter(i));
00214     printf("   %.1f, %.1f, ", hA->GetBinContent(i), hA->GetBinError(i));
00215     printf("   %.1f, %.1f, ", hCB->GetBinContent(i), hCB->GetBinError(i));
00216     printf("   %.1f, %.1f\n ", hDD->GetBinContent(i), hDD->GetBinError(i));
00217   }
00218 
00219 }