tau_mode.c

Go to the documentation of this file.

#include <stdlib.h>
#ifndef __CINT__
#include <TSQLServer.h>
#include <TSQLResult.h>
#include <TSQLRow.h>
#endif
 
//
// get  the online beam energy from the database for a specific run
// within ROOT
//          <[email protected]>
//
 
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 =event;
 
    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
    Int_t tot_time = data->GetMaximum("etsT1");
    Double_t difft;
    if ( tot_time > 0 ) 
    {
        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>108000) // 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;
    }
 
    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.491&&(e2/E_cms)>0.417&&(e2/E_cms)<0.491
                           &&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;
    
        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[11],lum[10];
            for(Int_t i=0;i<11;i++)
            {
                xtime[i]=timemin+(difft/10.)*i;
            }
    
    
            //Get 10 histograms of dltphi distribution between time nodes
            TH1D  *hdltphi[10];
                for(Int_t i=0;i<10;i++)
            {
                        hdltphi[i] = new TH1D("","dltphi distribution",50,-50.,50.);
                }
                TH1 *htime = new TH1D("htime","time", 80, -10., 10.);
 
            for(Int_t i=0;i<10;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[10];
            Double_t x[10];
            for(Int_t i=0;i<10;i++)
            {
                nevt[i] = hdltphi[i]->GetEntries();
                x[i]=xtime[i+1]-xtime[0];
                lum[i]=nevt[i];
            }
 
            //Set draw options and fit the luminosity curve
            TCanvas *c2 = new TCanvas("c2"," ",700,500);
            c2->SetFillColor(kWhite);
            c2->SetGrid();
    
            Int_t imin;
            for(Int_t i=0;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<10;i++)
            {
                lumsum += lum[i];
            }
            Double_t lumsum2 = lumsum*lumsum;
 
            const Int_t n = 10;
            TF1 *g1    = new TF1("g1","[0]*exp(-x/[1])",xmin,10000.);
            g1->SetParameters(lummax,1e+4.);
            g1->SetLineColor(2);
            TGraph* gr = new TGraph(n,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;
}