StRoot: StEEmcPool/StEEmcIUPi0/macros/Calall.C Source File

00001 //this macro is designed to calculate the single spin asymmetries from yellow and blue beams, and double spin asymmetry from pi0 study, based on the generated/fitted/normalized pi0 yield from FitPt2Mass.C in our fitting philosopy. 
00002 //author: Weihong He
00003 
00004 #include <iostream>
00005 #include <string>
00006 #include <fstream>
00007   using namespace std;
00008 //X:yellow Y:blue
00009 //background_LL is calculated from BackgroundLL.C because we want to use another method: side-band analysis to background study, although you can also calculate background all here by simply replacing the all array to back array.
00010 double epsilon_ll[7],ep_error[7],ptmean[7],epsilon_real[7],epsilon_back[7],ds[7],dr[7],dw[7],db[7],epsilon_X[7],epsilon_Y[7],dx[7],dy[7],ep_error_rb[7];
00011 Calall(){
00012   //x*[*] save epsilon from yellow beam, y*[*] save epsilon from blue beam, all*[*] save all values.
00013   double x1[7],x2[7],x3[7],y1[7],y2[7],y3[7],all1[7],all2[7],all3[7];
00014   //corresponding errors
00015   double ex1[7],ex2[7],ex3[7],ey1[7],ey2[7],ey3[7],eall1[7],eall2[7],eall3[7];
00016   double pt1[7],pt2[7],pt3[7];
00017   gStyle->SetPadGridX(0);
00018   gStyle->SetPadGridY(0);
00019   gStyle->SetCanvasColor(0);
00020   gStyle->SetOptStat(0);
00021   c0=new TCanvas("epsilony","epsilony",700,500);
00022   c0->Divide(1,1);
00023   c1=new TCanvas("epsilonx","epsilonx",700,500);
00024   c1->Divide(1,1);
00025   c2=new TCanvas("all","all",700,500);
00026   c2->Divide(1,1);
00027   TH1F* h1=new TH1F("epsilon_LB","#epsilon_{LB} vs pT",60,4.,14.);
00028   TH1F* h2=new TH1F("epsilon_LY","#epsilon_{LY} vs pT",60,4.,14.);
00029   //theory curve histograms
00030   TH1F* h5=new TH1F("A_LL","#vec{p}+#vec{p}->#pi^{0}X, #sqrt{s}=200Gev, 1.0#leq#eta#leq2.0;pT[Gev/c];A_{LL}",22,1.0,12.0);
00031   TH1F* h6=new TH1F("A_LL","#vec{p}+#vec{p}->#pi^{0}X, #sqrt{s}=200Gev, 1.089#leq#eta#leq2.0,P=60%;[Gev/c];A_{LL}",22,1.0,12.0);
00032   TH1F* h7=new TH1F("A_LL","#vec{p}+#vec{p}->#pi^{0}X, #sqrt{s}=200Gev, 1.089#leq#eta#leq2.0,P=60%;[Gev/c];A_{LL}",22,1.0,12.0);
00033   TH1F* h8=new TH1F("A_LL","#vec{p}+#vec{p}->#pi^{0}X, #sqrt{s}=200Gev, 1.089#leq#eta#leq2.0,P=60%;[Gev/c];A_{LL}",22,1.0,12.0);
00034 
00035 
00036   gStyle->SetOptStat(0); 
00037 
00038   //dataset 1
00039   ifstream yield1("yield3335557.txt");  
00040   doCalall(yield1);
00041   for(int j=0;j<7;j++){
00042     x1[j]=epsilon_X[j];
00043     ex1[j]=dx[j];
00044     y1[j]=epsilon_Y[j];
00045     ey1[j]=dy[j];
00046     all1[j]=epsilon_ll[j];
00047     eall1[j]=ep_error_rb[j];
00048     pt1[j]=ptmean[j];
00049 
00050   }
00051 
00052   //dataset 2
00053   ifstream yield2("yield4444446.txt");  
00054   doCalall(yield2);
00055   for(int j=0;j<7;j++){
00056     x2[j]=epsilon_X[j];
00057     ex2[j]=dx[j];
00058     y2[j]=epsilon_Y[j];
00059     ey2[j]=dy[j];
00060     all2[j]=epsilon_ll[j];
00061     eall2[j]=ep_error_rb[j];
00062     pt2[j]=ptmean[j];
00063 
00064   }
00065 
00066   //dataset3
00067   ifstream yield3("yield5556668.txt");  
00068   doCalall(yield3);
00069   for(int j=0;j<7;j++){
00070     x3[j]=epsilon_X[j];
00071     ex3[j]=dx[j];
00072     y3[j]=epsilon_Y[j];
00073     ey3[j]=dy[j];
00074     all3[j]=epsilon_ll[j];
00075     eall3[j]=ep_error_rb[j];
00076     pt3[j]=ptmean[j];
00077 
00078   }
00079   double big[7],small[7];
00080   for(int k=0;k<7;k++){
00081     if(all1[k]>all2[k]) big[k]=all1[k];
00082     else big[k]=all2[k];
00083     if(all3[k]>big[k]) big[k]=all3[k];
00084 
00085     if(all1[k]<all2[k]) small[k]=all1[k];
00086     else small[k]=all2[k];
00087     if(all3[k]<small[k]) small[k]=all3[k];
00088   }
00089   double sys_error[7];
00090   for(int m=0;m<7;m++){
00091     sys_error[m]=big[m]-small[m];
00092     //cout<<"syserror="<<sys_error[m]<<endl;
00093   }
00094   //calculate All
00095 
00096   const int nBins = 7;
00097   float bins[nBins+1] = {4,5,6,7,8,9,10,12};
00098 #if 1
00099   TH1F* h4 = new TH1F("A_LL","#vec{p}+#vec{p}->#pi^{0}X, #sqrt{s}=200Gev, 1.0#leq#eta#leq2.0;pT[Gev/c];A_{LL};",nBins,bins);
00100   h4->SetMarkerStyle(8);
00101 #endif
00102   TH1F* hSys = new TH1F("A_LL","#vec{p}+#vec{p}->#pi^{0}X, #sqrt{s}=200Gev, 1.0#leq#eta#leq2.0;pT[Gev/c];A_{LL};",nBins,bins);
00103   TH1F* hB = new TH1F("A_LL","#vec{p}+#vec{p}->#pi^{0}X, #sqrt{s}=200Gev, 1.0#leq#eta#leq2.0;pT[Gev/c];A_{LL};",nBins,bins);
00104   float baseLine=-0.17;
00105   float sysErr[7] = {0.00239,0.00521,0.00485,0.00405,0.012,0.0108,0.0104};
00106 
00107   //point 1
00108   h4->Fill(pt3[0],all3[0]);
00109   int nbin0=h4->FindBin(pt3[0]);
00110   h4->SetBinError(nbin0,eall3[0]);//+sys_error[0]);//+0.0005315);
00111   h1->Fill(pt3[0],y3[0]);
00112   int nbbin0=h1->FindBin(pt3[0]);
00113   h1->SetBinError(nbbin0,ey3[0]);
00114   h2->Fill(pt3[0],x3[0]);
00115   int nbinn0=h2->FindBin(pt3[0]);
00116   h2->SetBinError(nbinn0,ex3[0]);
00117   cout<<"point1 error="<<eall3[0]<<" all="<<all3[0]<<endl;
00118   hSys->SetBinContent(nbin0,baseLine + sysErr[0]);
00119   hB->SetBinContent(nbin0,baseLine);
00120   //hSys->SetBinError(nbin0,sysErr[0]/2.);
00121   //point 2
00122   h4->Fill(pt2[1],all2[1]);
00123   int nbin1=h4->FindBin(pt2[1]);
00124   h4->SetBinError(nbin1,eall2[1]);//+sys_error[1]);//+0.0002291);
00125   h1->Fill(pt2[1],y2[1]);
00126   int nbbin1=h1->FindBin(pt2[1]);
00127   h1->SetBinError(nbbin1,ey2[1]);
00128   h2->Fill(pt2[1],x2[1]);
00129   int nbinn1=h2->FindBin(pt2[1]);
00130   h2->SetBinError(nbinn1,ex2[1]);
00131   cout<<"point2 error="<<eall2[1]<<" all="<<all2[1]<<endl;
00132   hSys->SetBinContent(nbin1,baseLine + sysErr[1]);
00133   //hSys->SetBinError(nbin1,sysErr[1]/2.);
00134   hB->SetBinContent(nbin1,baseLine);
00135   //point 3
00136   h4->Fill(pt1[2],all1[2]);
00137   int nbin2=h4->FindBin(pt1[2]);
00138   h4->SetBinError(nbin2,eall1[2]);//+sys_error[2]);//+0.0005046);
00139   h1->Fill(pt1[2],y1[2]);
00140   int nbbin2=h1->FindBin(pt1[2]);
00141   h1->SetBinError(nbbin2,ey1[2]);
00142   h2->Fill(pt1[2],x1[2]);
00143   int nbinn2=h2->FindBin(pt1[2]);
00144   h2->SetBinError(nbinn2,ex1[2]);
00145   cout<<"point3 error="<<eall1[2]<<" all="<<all1[2]<<endl;
00146   hSys->SetBinContent(nbin2,baseLine + sysErr[2]);
00147   //hSys->SetBinError(nbin2,sysErr[2]/2.);
00148   hB->SetBinContent(nbin2,baseLine);
00149   //point 4
00150   h4->Fill(pt2[3],all2[3]);
00151   int nbin3=h4->FindBin(pt2[3]);
00152   h4->SetBinError(nbin3,eall2[3]);//+sys_error[3]);//+0.0004908);
00153   h1->Fill(pt2[3],y2[3]);
00154   int nbbin3=h1->FindBin(pt2[3]);
00155   h1->SetBinError(nbbin3,ey2[3]);
00156   h2->Fill(pt2[3],x2[3]);
00157   int nbinn3=h2->FindBin(pt2[3]);
00158   h2->SetBinError(nbinn3,ex2[3]);
00159   cout<<"point4 error="<<eall2[3]<<" all="<<all2[3]<<endl;
00160   hSys->SetBinContent(nbin3,baseLine + sysErr[3]);
00161   //hSys->SetBinError(nbin3,sysErr[3]/2.);
00162   hB->SetBinContent(nbin3,baseLine);
00163   //point 5
00164   h4->Fill(pt2[4],all2[4]);
00165   int nbin4=h4->FindBin(pt2[4]);
00166   h4->SetBinError(nbin4,eall2[4]);//+sys_error[4]);//+0.0002051);
00167   h1->Fill(pt2[4],y2[4]);
00168   int nbbin4=h1->FindBin(pt2[4]);
00169   h1->SetBinError(nbbin4,ey2[4]);
00170   h2->Fill(pt2[4],x2[4]);
00171   int nbinn4=h2->FindBin(pt2[4]);
00172   h2->SetBinError(nbinn4,ex2[4]);
00173   cout<<"point5 error="<<eall2[4]<<" all="<<all2[4]<<endl;
00174   hSys->SetBinContent(nbin4,baseLine + sysErr[4]);
00175   //hSys->SetBinError(nbin4,sysErr[4]/2.);
00176   hB->SetBinContent(nbin4,baseLine);
00177   //point 6
00178   h4->Fill(pt1[5],all1[5]);
00179   int nbin5=h4->FindBin(pt1[5]);
00180   h4->SetBinError(nbin5,eall1[5]);//+sys_error[5]);//+0.00055);
00181   h1->Fill(pt1[5],y1[5]);
00182   int nbbin5=h1->FindBin(pt1[5]);
00183   h1->SetBinError(nbbin5,ey1[5]);
00184   h2->Fill(pt1[5],x1[5]);
00185   int nbinn5=h2->FindBin(pt1[5]);
00186   h2->SetBinError(nbbin5,ex1[5]);
00187   cout<<"point6 error="<<eall1[5]<<" all="<<all1[5]<<endl;
00188   hSys->SetBinContent(nbin5,baseLine + sysErr[5]);
00189   //hSys->SetBinError(nbin5,sysErr[5]/2.);
00190   hB->SetBinContent(nbin5,baseLine);
00191   //point 7
00192   h4->Fill(pt2[6],all2[6]);
00193   int nbin6=h4->FindBin(pt2[6]);
00194   h4->SetBinError(nbin6,eall2[6]);//+sys_error[6]);//+0.001953);
00195   h1->Fill(pt2[6],y2[6]);
00196   int nbbin6=h1->FindBin(pt2[6]);
00197   h1->SetBinError(nbbin6,ey2[6]);
00198   h2->Fill(pt2[6],x2[6]);
00199   int nbinn6=h2->FindBin(pt2[6]);
00200   h2->SetBinError(nbinn6,ex2[6]);
00201   cout<<"point7 error="<<eall2[6]<<" all="<<all2[6]<<endl;
00202   hSys->SetBinContent(nbin6,baseLine + sysErr[6]);
00203   //hSys->SetBinError(nbin6,sysErr[6]/2.);
00204   hB->SetBinContent(nbin6,baseLine);
00205   //open theory curve
00206   ifstream unp("theorycurve/pion-unp-cteq6-rap1to2.dat");
00207   ifstream g0("theorycurve/pion-pol-g0-rap1to2.dat");
00208   ifstream plusg("theorycurve/pion-pol-max-rap1to2.dat");
00209   ifstream minusg("theorycurve/pion-pol-maxminus-rap1to2.dat");
00210   ifstream stdg("theorycurve/pion-pol-std-rap1to2.dat");
00211   double tx1=0.0,tx2=0.0,tx3=0.0,tx4=0.0,tx5=0.0;
00212   double td1=0.0,td2=0.0,td3=0.0,td4=0.0,td5=0.0;
00213   double tb1=0.0,tb2=0.0,tb3=0.0,tb4=0.0,tb5=0.0;
00214   double tc1=0.0,tc2=0.0,tc3=0.0,tc4=0.0,tc5=0.0;
00215   double ty1=0.0,ty2=0.0,ty3=0.0,ty4=0.0,ty5=0.0;
00216   int ccount=0;
00217   double theoryd[18],theoryb[18],theoryc[18],theoryg[18];;
00218   while(1){
00219     unp>>ty1>>ty2>>ty3>>ty4>>ty5;
00220     g0>>tx1>>tx2>>tx3>>tx4>>tx5;
00221     plusg>>td1>>td2>>td3>>td4>>td5;
00222     minusg>>tb1>>tb2>>tb3>>tb4>>tb5;
00223     stdg>>tc1>>tc2>>tc3>>tc4>>tc5;
00224     theoryd[ccount]=td5/ty5;
00225     theoryb[ccount]=tb5/ty5;
00226     theoryc[ccount]=tc5/ty5;
00227     theoryg[ccount]=tx5/ty5;
00228     //cout<<"x3="<<d3<<" theoryd="<<theoryg[ccount]<<endl;
00229     h5->Fill(td3,theoryd[ccount]);
00230     h6->Fill(tb3,theoryb[ccount]);
00231     h7->Fill(tc3,theoryc[ccount]);
00232     h8->Fill(tx3,theoryg[ccount]);
00233     ccount++;
00234     if(ccount>17) break;
00235   }
00236   h5->Fill(8.0,0.03);
00237   h5->SetAxisRange(0,25,"X");
00238   h5->Fill(9.0,0.0305);
00239   h5->Fill(10.0,0.031);
00240   h5->Fill(10.5,0.0312);
00241   h5->Fill(11.5,0.031);
00242   h6->Fill(8.0,-0.001);
00243   h6->Fill(9.0,-0.0035);
00244   h6->Fill(10.0,-0.0065);
00245   h6->Fill(10.5,-0.0075);
00246   h6->Fill(11.5,-0.01);
00247   h7->Fill(8.0,0.0072);
00248   h7->Fill(9.0,0.0079);
00249   h7->Fill(10.0,0.0087);
00250   h7->Fill(10.5,0.0091);
00251   h7->Fill(11.5,0.0099);
00252   h8->Fill(8.0,0.00152);
00253   h8->Fill(9.0,0.00183);
00254   h8->Fill(10.0,0.00204);
00255   h8->Fill(10.5,0.0022);
00256   h8->Fill(11.5,0.0026);
00257   //plotting epsilon_ll, h1 for blue beam, h2 for yellow beam
00258 #if 1
00259   c0->cd(1);
00260   h1->Draw();
00261   h1->Fit("pol0","R","",4.0,12.0);
00262   h1->SetMinimum(-0.05);
00263   h1->SetMaximum(0.05);
00264   TLegend* leg1 = new TLegend(0.15,0.65,0.5,0.85);
00265   leg1->SetHeader("STAR 2006 Preliminary #pi^{0}");
00266   leg1->SetBorderSize(0.0);
00267   leg1->SetFillColor(0.0);
00268   leg1->Draw();
00269   c0->Print("epsilony.gif");
00270   c1->cd(1);
00271   h2->Draw();
00272   h2->Fit("pol0","R","",4.0,12.0);
00273   h2->SetMinimum(-0.05);
00274   h2->SetMaximum(0.05);
00275   TLegend* leg2 = new TLegend(0.15,0.65,0.5,0.85);
00276   leg2->SetHeader("STAR 2006 Preliminary #pi^{0}");
00277   leg2->SetBorderSize(0.0);
00278   leg2->SetFillColor(0.0);
00279   leg2->Draw();
00280   c1->Print("epsilonx.gif");
00281 #endif
00282 
00283 #if 1
00284   //A_LL plotting
00285   c2->cd(1);
00286   h5->SetLineColor(kRed);
00287   h5->Draw("C");
00288   h5->SetMinimum(-0.18);
00289   h5->SetMaximum(0.15);
00290   hB->SetFillColor(11);
00291   hB->Draw("same");
00292   hSys->SetFillColor(10);
00293   hSys->Draw("same");
00294   h6->SetLineColor(kGreen);
00295   h6->Draw("Csame");
00296   h7->Draw("Csame");
00297   h8->SetLineColor(kBlue);
00298   h8->Draw("Csame");
00299 
00300   h4->Draw("E1same,p");
00301 
00302   TLegend* leg = new TLegend(0.15,0.6,0.4,0.8);
00303   leg->AddEntry(h5,"#Delta{G}=G","L");
00304   leg->AddEntry(h6,"#Delta{G}=-G","L");
00305   leg->AddEntry(h7,"#Delta{G}=std","L");
00306   leg->AddEntry(h8,"#Delta{G}=0","L");
00307   leg->SetHeader("Vogelsang Prediction(GRSV)");
00308   leg->SetBorderSize(0.0);
00309   leg->SetFillColor(0.0);
00310   leg->Draw();
00311   TLegend* leg3 = new TLegend(0.6,0.75,0.9,0.85);
00312   leg3->SetHeader("STAR 2006 Preliminary #pi^{0}");
00313   leg3->SetBorderSize(0.0);
00314   leg3->SetFillColor(0.0);
00315   leg3->Draw();
00316   c2->Print("all.gif");
00317 #endif
00318   
00319 }
00320 
00321 void doCalall(ifstream& yield){
00322   double epsilon_raw[7],weight[7],eps_lx_raw[7],eps_ly_raw[7],eps_lx_back[7],eps_ly_back[7];
00323   string s1,s2,s3,s4,s5;
00324   yield>>s1>>s2>>s3>>s4>>s5;
00325   double totsum=0,realsum=0,backsum=0;
00326   double psquare=0.3;
00327   //cout<<s1<<endl;
00328   for(int i=0;i<7;i++){
00329     int flag=0;
00330     int spin=0;
00331     double totpi=0,realpi=0,backpi=0;
00332     double pt=0.0;
00333     double Nuu_sig=0,Nuu_real=0,Nuu_back=0,Nud_sig=0,Nud_real=0,Nud_back=0,Ndu_sig=0,Ndu_real=0,Ndu_back=0,Ndd_sig=0,Ndd_real=0,Ndd_back=0;
00334     double sumtot=0,sumreal=0,sumback=0;
00335     double sumpt=0.0;
00336     double xx1=0.0,xx2=0.0,yy1=0.0,yy2=0.0;
00337     double des=0.0,deb=0.0,der=0.0,dew=0.0,eepsilon_raw=0.0,eepsilon_real=0.0,eeps_lx_raw=0.0,eeps_ly_raw=0.0,eeps_lx_back=0.0,eeps_ly_back=0.0;
00338     while(1)
00339       {
00340         yield>>spin>>pt>>totpi>>realpi>>backpi;
00341         //cout<<"spin="<<spin<<" pt="<<pt<<" totpi="<<totpi<<" realpi="<<realpi<<" backpi="<<backpi<<endl;
00342         sumtot+=totpi;
00343         sumreal+=realpi;
00344         sumback+=backpi;
00345         sumpt+=pt*totpi;
00346         totsum+=totpi;
00347         realsum+=realpi;
00348         backsum+=backpi;
00349         //cout<<"sumpt="<<sumpt<<endl;
00350         if(spin==0)
00351           {
00352             Nuu_sig+=totpi;
00353             Nuu_real+=realpi;
00354             Nuu_back+=backpi;
00355           }
00356         if(spin==1)
00357           {
00358             Nud_sig+=totpi;
00359             Nud_real+=realpi;
00360             Nud_back+=backpi;
00361           }
00362         if(spin==2)
00363           {
00364             Ndu_sig+=totpi;
00365             Ndu_real+=realpi;
00366             Ndu_back+=backpi;
00367           }
00368         if(spin==3)
00369           {
00370             Ndd_sig+=totpi;
00371             Ndd_real+=realpi;
00372             Ndd_back+=backpi;
00373           }
00374         flag++;
00375         if(flag==4) break;
00376         //if(flag==8) break;
00377       }//end of while
00378     xx1=Ndd_real+Ndu_real;
00379     xx2=Nud_real+Nuu_real;
00380     yy1=Ndd_real+Nud_real;
00381     yy2=Ndu_real+Nuu_real;
00382     epsilon_X[i]=-(xx1-xx2)/(xx1+xx2);
00383     epsilon_Y[i]=-(yy1-yy2)/(yy1+yy2);
00384     ptmean[i]=sumpt/sumtot;
00385 
00386     eepsilon_raw=(Nuu_sig-Nud_sig-Ndu_sig+Ndd_sig)/sumtot;
00387     eepsilon_real=(Nuu_real-Nud_real-Ndu_real+Ndd_real)/sumreal;
00388     eeps_lx_raw=(Ndd_sig-Nud_sig-Nuu_sig+Ndu_sig)/sumback;
00389     eeps_ly_raw=(Ndd_sig+Nud_sig-Nuu_sig-Ndu_sig)/sumback;
00390     eeps_lx_back=(Ndd_back-Nud_back-Nuu_back+Ndu_back)/sumback;
00391     eeps_ly_back=(Ndd_back+Nud_back-Nuu_back-Ndu_back)/sumback; 
00392     epsilon_raw[i]=1.0/psquare*eepsilon_raw;
00393     //all calculation from pure pi0 yield after background subtraction.
00394     epsilon_real[i]=1.0/psquare*eepsilon_real;
00395     eps_lx_raw[i]=1.0/psquare*eeps_lx_raw;
00396     eps_ly_raw[i]=1.0/psquare*eeps_ly_raw;
00397     eps_lx_back[i]=1.0/psquare*eeps_lx_back;
00398     eps_ly_back[i]=1.0/psquare*eeps_ly_back;
00399     weight[i]=sumback/sumtot;
00400 
00401     epsilon_back[i]=1.0/psquare*(Nuu_back-Nud_back-Ndu_back+Ndd_back)/sumback;
00402     //all calculation
00403     epsilon_ll[i]=(epsilon_raw[i]-weight[i]*epsilon_back[i])/(1.0-weight[i]);
00404     //stat. error calculation
00405     des=sqrt(sumtot)/sumtot;
00406     der=sqrt(sumreal)/sumreal;
00407     deb=sqrt(sumback)/sumback;
00408     dew=sqrt(sumback)/sumtot; 
00409     ds[i]=1.0/psquare*des;
00410     dr[i]=1.0/psquare*der;
00411     db[i]=1.0/psquare*deb;
00412     dw[i]=1.0/psquare*dew;
00413 
00414 
00415     ep_error[i]=sqrt(pow(ds[i]/(1-weight[i]),2)+pow((weight[i]*db[i])/(1-weight[i]),2)+pow(dw[i]*(epsilon_raw[i]-epsilon_back[i])/pow(1-weight[i],2),2));
00416     dx[i]=sqrt(pow(des/(1-weight[i]),2)+pow((weight[i]*deb)/(1-weight[i]),2)+pow(dew*(eeps_lx_raw-eeps_lx_back)/pow(1-weight[i],2),2));
00417     dy[i]=sqrt(pow(des/(1-weight[i]),2)+pow((weight[i]*deb)/(1-weight[i]),2)+pow(dew*(eeps_ly_raw-eeps_ly_back)/pow(1-weight[i],2),2));
00418     //the second way to calculate error of all
00419     ep_error_rb[i]=sqrt(pow(dr[i],2)+pow(db[i],2));
00420   
00421 
00422   }//end of for looping pt bins
00423  
00424 
00425 }