tau_mode.c

Go to the documentation of this file.

#include <stdlib.h>
#ifndef __CINT__
#  include <TSQLResult.h>
#  include <TSQLRow.h>
#  include <TSQLServer.h>
#endif
 
//
// get  the online beam energy from the database for a specific run
// within ROOT
//          <jixb@ihep.ac.cn>
//
 
const Int_t MaxRunTime = 108000;
Double_t    BeamEnergyOnline( Int_t );
 
void tau_mode( int flag = 0 ) {
  // Reset root
  gROOT->Reset();
 
  // Get the root file and E_cms
  TFile* f1 = new TFile( "YYYY/m_root_dir/LumTau_XXXX.root" );
  TTree* data;
  f1->GetObject( "event", data );
 
  if ( flag != 0 )
  {
    // Define the canvas
    TCanvas* myCanvas = new TCanvas();
    myCanvas->Divide( 1, 1 );
    TPad* c1_1 = new TPad( "c1_1", "c1_1", 0.01, 0.01, 0.99, 0.99 );
    c1_1->Draw();
    c1_1->cd();
    c1_1->Range( 0.425458, -114.842, 0.674993, 802.951 );
    c1_1->SetBorderSize( 2 );
    c1_1->SetBottomMargin( 0.125129 );
    c1_1->SetFrameFillColor( 0 );
  }
 
  // Initialize the variables and set the branches' address
  Bool_t   topup = false, is_topup = false;
  Int_t    run  = 0;
  Double_t time = 0.0, etsT1 = 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;
  Int_t nentries = 0;
 
  data->SetBranchAddress( "topup", &topup );
  data->SetBranchAddress( "run", &run );
  data->SetBranchAddress( "time", &time );
  data->SetBranchAddress( "etsT1", &etsT1 );
  data->SetBranchAddress( "dltphi", &dltphi );
  data->SetBranchAddress( "costht1", &costht1 );
  data->SetBranchAddress( "costht2", &costht2 );
  data->SetBranchAddress( "e1", &e1 );
  data->SetBranchAddress( "e2", &e2 );
  data->SetBranchAddress( "etot", &etot );
  data->SetBranchAddress( "phi1", &phi1 );
  data->SetBranchAddress( "phi2", &phi2 );
 
  nentries = (Int_t)data->GetEntries();
 
  Int_t    runno   = 0;
  Double_t timemin = 0.0, timemax = 0.0;
  Int_t    init_time = 1234567890; // corresponding to 2009-2-14 7:31:30
  Int_t    n_first   = 0;          // the number of the good first entry
  do {
    data->GetEntry( n_first );
    // timemin = time;
    n_first += 1;
  } while ( time < init_time && n_first < nentries );
  runno          = run;
  Double_t E_cms = 2 * BeamEnergyOnline( run );
  cout << "E_cms = " << E_cms << endl;
  is_topup = topup;
  timemin  = time;
  // get the total time of the run
  // etsT1 has some problem before round 13 data, so do not use it
  Int_t tot_time = 0;
  if ( run >= 63075 )
    tot_time = data->GetMaximum( "etsT1" ); // 63075 is the 1st run of round 13 data
  Double_t difft;
  if ( tot_time > 0 && tot_time < MaxRunTime ) { difft = tot_time; }
  else
  {                              // etsT1 can not be used, use another method
    Int_t n_last = nentries - 1; // last event in the ntuple
    // Deal with bad events with wrong time value.
    // Check the contents of the event can judge whether it is a good event
    do {
      data->GetEntry( n_last );
      n_last -= 1;
    } while ( etot <= 0 && n_last >= n_first );
    timemax = time;
    if ( timemin < init_time || timemax < init_time )
    { // protection
      cout << "\n the time of the  event is wrong!";
      cout << "time range : " << timemin << " - " << timemax << endl;
      return;
    }
    difft = timemax - timemin;
  }
 
  Double_t dlt = 0.0, mean = 0.0, width = 0.0;
  Double_t luminitial = 0.;
  Double_t tau        = 5000.;
 
  if ( difft <= 0. ||
       difft > MaxRunTime ) // protection. 30h. The time of a run shold < 10 hours
  {
    cout << "\nThe runtime is abnormal!!!" << endl;
    cout << "It's value is " << difft << endl;
    cout << "Please check it carefully!" << endl;
    return;
  }
 
  // Set draw options and fit the luminosity curve
  TCanvas* c2 = new TCanvas( "c2", " ", 700, 500 );
  c2->SetFillColor( kWhite );
  c2->SetGrid();
 
  if ( difft < 300. )
  { // the run is too short
    cout << "\nRuntime less than 300 " << endl;
    cout << "Using default tau value " << tau << endl;
  }
  else if ( is_topup )
  { // topup mode
    cout << "\n topup mode." << endl;
    luminitial = 2000;
    tau        = 99999999;
  }
  else
  {
    // Fill the histogram of dltphi distribution
    TH1F* dltphi1 = new TH1F( "dltphi1", "dltphi", 150, -30, 30 );
    for ( Int_t i = 0; i < nentries; i++ )
    {
      data->GetEntry( i );
      bool cut = ( e1 / E_cms ) > 0.417 && ( e1 / E_cms ) < 0.55 && ( e2 / E_cms ) > 0.417 &&
                 ( e2 / E_cms ) < 0.55 && fabs( costht1 ) < 0.8 && fabs( costht2 ) < 0.8 &&
                 etot < E_cms + 0.5;
      if ( cut ) { dltphi1->Fill( dltphi ); }
    }
 
    // Fit the Di-photon peak and get the value of dlt
    TF1* h1 = new TF1( "h1", "gaus", -6, 4 );
    dltphi1->Fit( h1, "R+" );
    dlt = h1->GetParameter( 1 );
 
    // Fit the Bhabha peak and get the values of mean and width
    TF1* h2 = new TF1( "h2", "gaus", 4, 30 );
    dltphi1->Fit( h2, "R+" );
    mean           = h2->GetParameter( 1 );
    Double_t sigma = h2->GetParameter( 2 );
    width          = 3 * sigma;
 
    const Int_t nbin = 10; // number of bin
    if ( dlt < -6 || dlt > 4 || mean < 4 || mean > 30 || ( mean - dlt ) - 0.5 * width < 0 )
      cout << "\nSomething wrong with cut Bhabha entries. " << endl;
    else
    {
      // Get 11 time nodes
      Double_t xtime[nbin + 1], lum[nbin];
      for ( Int_t i = 0; i < nbin + 1; i++ ) { xtime[i] = timemin + ( difft / 10. ) * i; }
 
      // Get 10 histograms of dltphi distribution between time nodes
      TH1D* hdltphi[nbin];
      for ( Int_t i = 0; i < nbin; i++ )
      { hdltphi[i] = new TH1D( "", "dltphi distribution", 50, -50., 50. ); }
      TH1D* htime = new TH1D( "htime", "time", 80, -10., 10. );
 
      for ( Int_t i = 0; i < nbin; i++ )
      {
        for ( Int_t j = 0; j < nentries; j++ )
        {
          data->GetEntry( j );
          bool lumcut = ( e1 / E_cms ) > 0.417 && ( e2 / E_cms ) > 0.417 &&
                        fabs( costht1 ) < 0.8 && fabs( costht2 ) < 0.8 &&
                        fabs( dltphi - dlt ) > ( mean - dlt ) - width &&
                        fabs( dltphi - dlt ) < ( mean - dlt ) + width && etot < E_cms + 0.5;
          if ( lumcut && time > xtime[i] && time < xtime[i + 1] )
          { hdltphi[i]->Fill( dltphi ); }
        }
      }
 
      // Get the average luminosity between time nodes
      Int_t    nevt[nbin];
      Double_t x[nbin];
      for ( Int_t i = 0; i < nbin; i++ )
      {
        nevt[i] = hdltphi[i]->GetEntries();
        x[i]    = xtime[i + 1] - xtime[0];
        lum[i]  = nevt[i];
      }
 
      Int_t imin = 0;
      for ( Int_t i = 0; i < nbin - 1 && lum[i] < lum[i + 1]; i++ ) { imin = i + 1; }
      Double_t lummax = lum[imin];
      Double_t xmin   = x[imin];
      Double_t lumsum;
      for ( Int_t i = imin; i < nbin; i++ ) { lumsum += lum[i]; }
      Double_t lumsum2 = lumsum * lumsum;
 
      TF1* g1 = new TF1( "g1", "[0]*exp(-x/[1])", xmin, 10000. );
      g1->SetParameters( lummax, 1.e+4 );
      g1->SetLineColor( 2 );
      TGraph* gr = new TGraph( nbin, x, lum );
      gr->SetLineColor( 2 );
      gr->SetLineWidth( 2 );
      gr->SetMarkerColor( 4 );
      gr->SetMarkerStyle( 21 );
      gr->SetTitle( "BbLum" );
      gr->GetXaxis()->SetTitle( "Time" );
      gr->GetYaxis()->SetTitle( "Luminosity" );
      gr->Fit( "g1", "R" );
      gr->Draw( "ap" );
 
      // pick the fit result
      luminitial = g1->GetParameter( 0 );
      tau        = g1->GetParameter( 1 );
    }
  }
 
  if ( !is_topup && tau > 99999. )
  {
    cout << "\nThe result of fit is terrible. " << endl << endl;
    cout << "set the tau value = 99999" << endl;
    tau = 99999;
  }
 
  // Export the fit result
  std::ofstream m_outputFile;
  m_outputFile.open( "YYYY/m_txt_dir/LumTau_XXXX.txt", ios_base::app );
  m_outputFile << runno << "        " << difft << "     " << luminitial
               << "     " << ( luminitial * exp( -1.0 * difft / ( tau ) ) )
               << "     " << tau << endl;
  m_outputFile.close();
 
  if ( flag != 0 )
  {
    // Get the value of tau and print it
    TString result = Form( "tau = %6.2f", tau );
    TText*  text01 = new TText( 0.7, 0.8, result );
    text01->SetNDC();
    text01->Draw( "same" );
    c2->Update();
    c2->GetFrame()->SetFillColor( kWhite );
    c2->GetFrame()->SetBorderSize( 1 );
    c2->Modified();
    c2->Print( "lum_0000XXXX.ps" );
  }
 
  if ( flag != 0 )
  {
    // Print some parameters
    cout << "E_cms = " << E_cms << endl;
    cout << "\nruntime = " << difft << endl;
    cout << "\ndlt = " << dlt << endl;
    cout << "mean = " << mean << endl;
    cout << "width = " << width << endl;
    cout << "\ntau = " << tau << endl << endl;
  }
}
 
Double_t BeamEnergyOnline( Int_t runno ) {
  Double_t dbe = 0; // initialization
 
  // connect to the BESIII database
  TSQLServer* db = TSQLServer::Connect( "mysql://bes3db2.ihep.ac.cn", "guest", "guestpass" );
  //  printf("Server info: %s\n", db->ServerInfo());
 
  // Int_t runno = 33743;
  // read beam energy from the online database 'run'
  const char* ins = "select BER_PRB from run.RunParams where run_number = %d";
  char        sql[4096];
  sprintf( sql, ins, runno );
  TSQLResult* res = db->Query( sql );
  if ( res->GetRowCount() == 0 )
  { // not found
    printf( "Cannot find the beam energy of run %d from the database.\n", runno );
    delete res;
    delete db;
    return dbe;
  }
  TSQLRow* row = res->Next();
  TString  sbe = row->GetField( 0 ); // beam energy of electron
  dbe          = atof( sbe );        // convert string to float
  //  printf("beam energy is : %5.3f\n", dbe);
  delete row;
  delete res;
  delete db;
 
  return dbe;
}