00001 #include <stdlib.h>
00002 #ifndef __CINT__
00003 #include <TSQLServer.h>
00004 #include <TSQLResult.h>
00005 #include <TSQLRow.h>
00006 #endif
00007
00008
00009
00010
00011
00012
00013
00014 Double_t BeamEnergyOnline(Int_t);
00015
00016 void tau_mode(int flag=0)
00017 {
00018
00019 gROOT->Reset();
00020
00021
00022 TFile *f1 = new TFile("YYYY/m_root_dir/LumTau_XXXX.root");
00023 TTree *data =event;
00024
00025 if (flag!=0)
00026 {
00027
00028 TCanvas *myCanvas = new TCanvas();
00029 myCanvas->Divide(1,1);
00030 TPad *c1_1 = new TPad("c1_1", "c1_1",0.01,0.01,0.99,0.99);
00031 c1_1->Draw();
00032 c1_1->cd();
00033 c1_1->Range(0.425458,-114.842,0.674993,802.951);
00034 c1_1->SetBorderSize(2);
00035 c1_1->SetBottomMargin(0.125129);
00036 c1_1->SetFrameFillColor(0);
00037 }
00038
00039
00040 Int_t run = 0;
00041 Double_t time = 0.0,dltphi = 0.0,costht1 = 0.0,costht2 = 0.0,e1 = 0.0,e2 = 0.0,etot = 0.0,phi1 = 0.0,phi2 = 0.0;
00042 Int_t nentries = 0;
00043
00044 data->SetBranchAddress("run",&run);
00045 data->SetBranchAddress("time",&time);
00046 data->SetBranchAddress("dltphi",&dltphi);
00047 data->SetBranchAddress("costht1",&costht1);
00048 data->SetBranchAddress("costht2",&costht2);
00049 data->SetBranchAddress("e1",&e1);
00050 data->SetBranchAddress("e2",&e2);
00051 data->SetBranchAddress("etot",&etot);
00052 data->SetBranchAddress("phi1",&phi1);
00053 data->SetBranchAddress("phi2",&phi2);
00054
00055 data->GetEntry(0);
00056 Int_t runno = run;
00057 Double_t E_cms = 2*BeamEnergyOnline(run);
00058 cout << "E_cms = " << E_cms << endl;
00059
00060 nentries = (Int_t)data->GetEntries();
00061
00062 Double_t timemin = 0.0,timemax = 0.0;
00063 data->GetEntry(0);
00064 timemin = time;
00065 data->GetEntry((nentries - 1));
00066 timemax = time;
00067 Double_t difft = timemax-timemin;
00068
00069 Double_t dlt = 0.0,mean = 0.0,width = 0.0;
00070 Double_t luminitial = 0.;
00071 Double_t tau = 5000.;
00072
00073 if (difft<300.)
00074 cout << "\nRuntime less than 300 " << endl;
00075 else
00076 {
00077
00078 TH1F *dltphi1 = new TH1F("dltphi1","dltphi",150,-30,30);
00079 for (Int_t i = 0;i<nentries;i++)
00080 {
00081 data->GetEntry(i);
00082 bool cut = (e1/E_cms)>0.417&&(e1/E_cms)<0.491&&(e2/E_cms)>0.417&&(e2/E_cms)<0.491
00083 &&fabs(costht1)<0.8&&fabs(costht2)<0.8
00084 &&etot<E_cms+0.5;
00085 if (cut)
00086 {
00087 dltphi1->Fill(dltphi);
00088 }
00089 }
00090
00091
00092 TF1* h1=new TF1("h1","gaus",-6,4);
00093 dltphi1->Fit(h1,"R+");
00094 dlt = h1->GetParameter(1);
00095
00096
00097 TF1* h2=new TF1("h2","gaus",4,30);
00098 dltphi1->Fit(h2,"R+");
00099 mean = h2->GetParameter(1);
00100 Double_t sigma = h2->GetParameter(2);
00101 width = 3*sigma;
00102
00103 if (dlt<-6||dlt>4||mean<4||mean>30||(mean-dlt)-0.5*width<0)
00104 cout << "\nSomething wrong with cut Bhabha entries. " << endl;
00105 else
00106 {
00107
00108 Double_t xtime[11],lum[10];
00109 for(Int_t i=0;i<11;i++)
00110 {
00111 xtime[i]=timemin+(difft/10.)*i;
00112 }
00113
00114
00115
00116 TH1D *hdltphi[10];
00117 for(Int_t i=0;i<10;i++)
00118 {
00119 hdltphi[i] = new TH1D("","dltphi distribution",50,-50.,50.);
00120 }
00121 TH1 *htime = new TH1D("htime","time", 80, -10., 10.);
00122
00123 for(Int_t i=0;i<10;i++)
00124 {
00125 for(Int_t j=0;j<nentries;j++)
00126 {
00127 data->GetEntry(j);
00128 bool lumcut = (e1/E_cms)>0.417&&(e2/E_cms)>0.417&&fabs(costht1)<0.8&&
00129 fabs(costht2)<0.8&&fabs(dltphi-dlt)>(mean-dlt)-width&&
00130 fabs(dltphi-dlt)<(mean-dlt)+width&&etot<E_cms+0.5;
00131 if(lumcut&&time>xtime[i]&&time<xtime[i+1])
00132 {
00133 hdltphi[i]->Fill(dltphi);
00134 }
00135 }
00136 }
00137
00138
00139 Int_t nevt[10];
00140 Double_t x[10];
00141 for(Int_t i=0;i<10;i++)
00142 {
00143 nevt[i] = hdltphi[i]->GetEntries();
00144 x[i]=xtime[i+1]-xtime[0];
00145 lum[i]=nevt[i];
00146 }
00147
00148
00149 TCanvas *c2 = new TCanvas("c2"," ",700,500);
00150 c2->SetFillColor(kWhite);
00151 c2->SetGrid();
00152
00153 Int_t imin;
00154 for(Int_t i=0;lum[i]<lum[i+1];i++)
00155 {
00156 imin = i+1;
00157 }
00158 Double_t lummax = lum[imin];
00159 Double_t xmin = x[imin];
00160 Double_t lumsum;
00161 for(Int_t i=imin;i<10;i++)
00162 {
00163 lumsum += lum[i];
00164 }
00165 Double_t lumsum2 = lumsum*lumsum;
00166
00167 const Int_t n = 10;
00168 TF1 *g1 = new TF1("g1","[0]*exp(-x/[1])",xmin,10000.);
00169 g1->SetParameters(lummax,1e+4.);
00170 g1->SetLineColor(2);
00171 TGraph* gr = new TGraph(n,x,lum);
00172 gr->SetLineColor(2);
00173 gr->SetLineWidth(2);
00174 gr->SetMarkerColor(4);
00175 gr->SetMarkerStyle(21);
00176 gr->SetTitle("BbLum");
00177 gr->GetXaxis()->SetTitle("Time");
00178 gr->GetYaxis()->SetTitle("Luminosity");
00179 gr->Fit("g1","R");
00180 gr->Draw("ap");
00181
00182
00183 luminitial = g1->GetParameter(0);
00184 tau = g1->GetParameter(1);
00185 }
00186 }
00187
00188
00189 if (tau>99999.)
00190 {
00191 cout << "\nThe result of fit is terrible. " << endl << endl;
00192 cout << "set the tau value = 99999" << endl;
00193 tau = 99999;
00194 }
00195
00196 if (difft<0.)
00197 {
00198 cout << "\nThe runtime is minus." << endl << endl;
00199 difft = fabs(difft);
00200 }
00201
00202
00203 std::ofstream m_outputFile;
00204 m_outputFile.open("YYYY/m_txt_dir/LumTau_XXXX.txt", ios_base::app);
00205 m_outputFile<< runno <<" "<<difft<<" "<<luminitial<<" "<<(luminitial * exp(- 1.0 * difft / (tau)))<<" "<< tau <<endl;
00206 m_outputFile.close();
00207
00208 if (flag!=0)
00209 {
00210
00211 TString result = Form("tau = %6.2f",tau);
00212 TText* text01 = new TText(0.7,0.8,result);
00213 text01->SetNDC();
00214 text01->Draw("same");
00215 c2->Update();
00216 c2->GetFrame()->SetFillColor(kWhite);
00217 c2->GetFrame()->SetBorderSize(1);
00218 c2->Modified();
00219 c2->Print("lum_0000XXXX.ps");
00220 }
00221
00222 if (flag!=0)
00223 {
00224
00225 cout << "E_cms = " << E_cms << endl;
00226 cout << "\nruntime = " << difft << endl;
00227 cout << "\ndlt = " << dlt << endl;
00228 cout << "mean = " << mean << endl;
00229 cout << "width = " << width << endl;
00230 cout << "\ntau = " << tau << endl << endl;
00231 }
00232 }
00233
00234
00235
00236 Double_t BeamEnergyOnline(Int_t runno)
00237 {
00238 Double_t dbe = 0;
00239
00240
00241 TSQLServer *db = TSQLServer::Connect("mysql://bes3db2.ihep.ac.cn", "guest", "guestpass");
00242
00243
00244
00245
00246 const char *ins = "select BER_PRB from run.RunParams where run_number = %d";
00247 char sql[4096];
00248 sprintf(sql, ins, runno);
00249 TSQLResult *res = db->Query(sql);
00250 if (res->GetRowCount() == 0) {
00251 printf("Cannot find the beam energy of run %d from the database.\n", runno);
00252 delete res;
00253 delete db;
00254 return dbe;
00255 }
00256 TSQLRow *row = res->Next();
00257 TString sbe = row->GetField(0);
00258 dbe = atof(sbe);
00259
00260 delete row;
00261 delete res;
00262 delete db;
00263
00264 return dbe;
00265 }
00266
00267