/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/LumTauAlg/LumTauAlg-00-00-08/share/tau_mode.c File Reference

#include <stdlib.h>
#include <TSQLServer.h>
#include <TSQLResult.h>
#include <TSQLRow.h>

Go to the source code of this file.

Functions

Double_t BeamEnergyOnline (Int_t)
void tau_mode (int flag=0)


Function Documentation

Double_t BeamEnergyOnline ( Int_t   ) 

Definition at line 236 of file tau_mode.c.

Referenced by tau_mode().

00237 {
00238         Double_t dbe = 0;             // initialization
00239 
00240         // connect to the BESIII database
00241         TSQLServer *db = TSQLServer::Connect("mysql://bes3db2.ihep.ac.cn", "guest", "guestpass");
00242         //  printf("Server info: %s\n", db->ServerInfo());
00243         
00244         //Int_t runno = 33743;
00245         //read beam energy from the online database 'run'
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) { // not found
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); // beam energy of electron
00258         dbe = atof(sbe);              // convert string to float
00259         //  printf("beam energy is : %5.3f\n", dbe);
00260         delete row;
00261         delete res;
00262         delete db;
00263 
00264         return dbe;
00265 }

void tau_mode ( int  flag = 0  ) 

Definition at line 16 of file tau_mode.c.

References BeamEnergyOnline(), cut, exp(), f1, genRecEmupikp::i, ganga-rec::j, phi1, phi2, check_raw_filter::runno, and x.

00017 {
00018         //Reset root
00019         gROOT->Reset();
00020 
00021         //Get the root file and E_cms
00022         TFile *f1 = new TFile("YYYY/m_root_dir/LumTau_XXXX.root");
00023         TTree *data =event;
00024 
00025         if (flag!=0)
00026         {
00027                 //Define the canvas
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         //Initialize the variables and set the branches' address
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                 //Fill the histogram of dltphi distribution
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                 //Fit the Di-photon peak and get the value of dlt
00092                 TF1* h1=new TF1("h1","gaus",-6,4);
00093                 dltphi1->Fit(h1,"R+");
00094                 dlt = h1->GetParameter(1);
00095 
00096                 //Fit the Bhabha peak and get the values of mean and width
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                         //Get 11 time nodes
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                         //Get 10 histograms of dltphi distribution between time nodes
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                         //Get the average luminosity between time nodes
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                         //Set draw options and fit the luminosity curve
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                         //pick the fit result
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         //Export the fit result
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                 //Get the value of tau and print it
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                 //Print some parameters
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 }


Generated on Tue Nov 29 23:16:42 2016 for BOSS_7.0.2 by  doxygen 1.4.7