MBDLaser/macros/plot_deltT_vs_Timestamp.C

0001 #include <iostream>
0002 #include <fstream>
0003 #include <string>
0004 #include <TSQLServer.h>
0005 #include <TSQLResult.h>
0006 #include <TSQLRow.h>
0007 #include <TSystem.h>
0008 #include <TDatime.h>
0009 #include <fstream>
0010 #include <TCanvas.h>
0011 #include <TFile.h>
0012 #include <TGraph.h>
0013 #include <TH2F.h>
0014 #include <TDatime.h>
0015
0016 #include <fun4all/Fun4AllUtils.h>
0017 R__LOAD_LIBRARY(libfun4all.so)
0018
0019   // varaibles needed
0020   const int MAXRUNS = 20000;
0021   TFile *savefile;
0022   TFile *tfile[MAXRUNS];
0023   TTree *t;
0024   Double_t runindex[MAXRUNS];
0025   TCanvas *ac[100];
0026   TFile *refFile_all;
0027   TString dir;
0028
0029   //declare TGraph
0030   TGraph *graph;
0031   TGraph *gt_chBych;
0032   TGraph *gq_chBych;
0033
0034   //Time in all runs
0035   TH2 *h2_mbdtt[MAXRUNS];
0036   Double_t T_current;
0037   TFile *refFile;
0038   TH2 *h2_ref;
0039   Double_t T_ref;
0040   Double_t deltaT;
0041
0042   //charge in each ch
0043   TH1 *h_mbdq[MAXRUNS][128];
0044   Double_t Q_current[MAXRUNS][128];
0045   TH1 *h_ref_Q[128];
0046   Double_t Q_ref[128];
0047   Double_t Q_ratio[MAXRUNS][128];
0048
0049   //Time in each ch
0050   TH1 *h_mbdt[MAXRUNS][128];
0051   //Double_t btmean[128][MAXRUNS];      // bt mean in each ch
0052   //Double_t btmeanerr[128][MAXRUNS];
0053   Double_t T_current_all[MAXRUNS][128];
0054   TH1 *h_ref_all[128];
0055   Double_t T_ref_all[128];
0056   Double_t deltaT_all[MAXRUNS][128];
0057
0058
0059   //set of variable to get the time from database
0060   UInt_t runSeconds; //Double_t runSeconds; this wrong insted I should use UInt_t
0061   Int_t startSeconds;
0062
0063   //Canvas
0064   TCanvas *g_Canvas;
0065   TCanvas *gt_Canvas;
0066   TCanvas *gq_Canvas;
0067
0068 unsigned int get_time( Double_t run_number)
0069 {
0070   // connect to the database
0071   TString servername = "pgsql://sphnxdaqdbreplica.sdcc.bnl.gov:5432/daq";
0072   TString username = "phnxrc";
0073   TString password = "";
0074   TSQLServer* dbserver = TSQLServer::Connect(servername, username, password);
0075
0076   if (!dbserver)
0077   {
0078     cout << "Failed to connect " << endl;
0079     return -1;
0080   }
0081
0082   // check the connection status
0083   if (dbserver->IsConnected())
0084   {
0085     cout << "Connected to the db" <<  endl;
0086   }
0087   else
0088   {
0089     cout << "Failed to connect to db" <<  endl;
0090     return -1;
0091   }
0092
0093   TDatime dt;
0094
0095   // query the database for the timestamp using the run_number number; constructs a SQL query to retriveve the ctime column from the filelist table in db, basically the query searches for rows where the filename column contains the extracted sRUN_number then stored the query result in res variable
0096   TString query = "SELECT ctime FROM filelist WHERE filename LIKE '%";
0097   query += run_number;
0098   query += "%';";
0099
0100
0101   TSQLResult* res = dbserver->Query(query);
0102
0103   //checking that the query is not null (res !=nullptr), then if there is a result it retieves the next row using res->Next() function and then if a row exisit it retrieves the value of the first filed in the row which is index 0 and stores it in the timestampe variable
0104   if (res==nullptr)
0105   {
0106     cout << "Failed to query the database for runnumber " << run_number <<  endl;
0107     return 0;
0108   }
0109
0110   TSQLRow* row = res->Next();
0111   if (row==nullptr)
0112   {
0113     cout << "No entry found in the database for runnumber " << run_number <<  endl;
0114     delete res;
0115     return 0;
0116   }
0117
0118   dt.Set( row->GetField(0) );
0119
0120   // checking
0121   //cout << "Timestamp for run " << run_number << ": " << dt.AsString() << "\t"<< " field = " <<  row->GetField(0) << endl;
0122
0123   delete row;
0124   delete res;
0125
0126   // disconnect from the database as Mickey asked
0127   dbserver->Close();
0128   delete dbserver;
0129
0130   return dt.Convert();
0131
0132 }
0133
0134 void plot_deltT_vs_time(const TDatime& dt, const char *flistname = "analyzed_Laser_data2023_To_current.list")
0135 {
0136   ifstream flist(flistname);
0137   TString dstfname;
0138   int nruns = 0;
0139   TString name;
0140   TString title;
0141   int icv = 0;
0142   savefile = new TFile(dir + "Tlaser.root", "RECREATE");
0143
0144   //initialize TGraph
0145   graph = new TGraph();
0146   gt_chBych = new TGraph();
0147   gq_chBych = new TGraph();
0148
0149
0150   // start time (Tzero) is the timestamp from the test_ctime function
0151   TDatime startTime(dt);
0152
0153   // calculate the time in seconds since the start time based on the timestamp from the test_ctime function
0154   TDatime runTime;
0155   startTime.Convert();
0156
0157   // looping over all the root file in the list and doing the analysis
0158   while (flist >> dstfname)
0159   {
0160     tfile[nruns] = new TFile(dstfname, "READ");
0161     title = dstfname;
0162     // cout << "Title = " << dstfname << endl;
0163
0164     h2_mbdtt[nruns] = (TH2 *)tfile[nruns]->Get("h2_mbdtt");
0165
0166     if (h2_mbdtt[nruns] != nullptr)
0167     {
0168       // get the time of the laser in the current run
0169       T_current = h2_mbdtt[nruns]->GetMean(1);
0170
0171       // get the time of the laser in the reference run
0172       refFile = new TFile("LASER_UNCALMBD-00040759-0000.root", "READ");
0173       h2_ref = (TH2 *)refFile->Get("h2_mbdtt");
0174       T_ref = h2_ref->GetMean(1);
0175
0176       // calculate delta-T to feed the TGraph Y-axis
0177       deltaT = T_current - T_ref;
0178
0179       //seting the runTime
0180       runTime.Set(dt.Convert() + nruns * 86400); // add one day (86400 seconds) for each run
0181
0182
0183       //seting the run_number using the GetRunSegment() function in Fun4AllUtils class which returns a pair of integers and assigned the first elements of this pair to run_number
0184       pair<int, int> runseg = Fun4AllUtils::GetRunSegment(dstfname.Data());
0185       int run_number = runseg.first;
0186
0187
0188
0189       //calling the get_time function
0190       runSeconds = get_time( run_number); // this function return time and pass run number
0191
0192
0193       //cout<< " runsSconds = " << runSeconds << "\t" << "dt.convert = " << dt.Convert() << "\t" << " runTime.Convert =" << runTime.Convert()<< endl;
0194
0195       // skip if the time can't be found
0196       if ( runSeconds==0 )
0197       {
0198         cout << "ERROR, no timestamp for run " << run_number << endl;
0199         continue;
0200       }
0201
0202       //fill
0203       graph->SetPoint(nruns,(Double_t)runSeconds, deltaT);
0204
0205       delete refFile;
0206
0207       // to get the meean and the mean error from  North MBD/MBD time && South MBD/MBD time
0208       for (int ipmt=0; ipmt<128; ipmt++)
0209       {
0210         // appending the value of ipmt to the string h_mbdq and then assigns it to the varaible name
0211         name = "h_mbdt"; name += ipmt; // h_mbdt0, h_mbdt1, ... so on
0212
0213         h_mbdt[nruns][ipmt] = (TH1*)tfile[nruns]->Get(name);
0214
0215         // btmean[ipmt][nruns] = h_mbdt[nruns][ipmt]->GetMean();
0216         // btmeanerr[ipmt][nruns] = h_mbdt[nruns][ipmt]->GetMeanError();
0217
0218         name = "h_mbdq"; name += ipmt;
0219         h_mbdq[nruns][ipmt] = (TH1*)tfile[nruns]->Get(name);
0220
0221
0222         //cout<< name<<" \t"<< btmean[ipmt][nruns]<< "\t"<< btmeanerr[ipmt][nruns]<<endl;
0223
0224         if (h_mbdt[nruns][ipmt] != nullptr)
0225         {
0226
0227           // get the time of the laser in the current run for each channls (we should have 128 values )
0228           T_current_all[nruns][ipmt] = h_mbdt[nruns][ipmt]->GetMean(1);  //cout<< T_current_all << "\t"<< name<< " \t" << ipmt<<endl;
0229
0230           Q_current[nruns][ipmt] = h_mbdq[nruns][ipmt]->GetMean(1);
0231
0232           // get the time of the laser in the reference run for each ch
0233           refFile_all = new TFile("LASER_UNCALMBD-00041145-0000.root", "READ");
0234
0235           name = "h_mbdt"; name += ipmt;
0236           h_ref_all[ipmt] = (TH1*)refFile_all->Get(name);
0237           T_ref_all[ipmt] = h_ref_all[ipmt]->GetMean(1);
0238
0239
0240           name = "h_mbdq"; name += ipmt;
0241           h_ref_Q[ipmt] = (TH1*)refFile_all->Get(name);
0242           Q_ref[ipmt] = h_ref_Q[ipmt]->GetMean(1);
0243
0244
0245
0246           //calculation
0247           deltaT_all[nruns][ipmt] = T_current_all[nruns][ipmt] - T_ref_all[ipmt];
0248
0249           Q_ratio[nruns][ipmt] = Q_current[nruns][ipmt] / Q_ref[ipmt];
0250
0251
0252           //cout<< "deltaT = " << deltaT_all[nruns][ipmt]<< " \t" <<"T_current in each ch = " << T_current_all[nruns][ipmt] << "\t" << "T_ref in each ch = "<<T_ref_all[ipmt] << " \t" << name <<"\t"<< nruns<< endl;
0253
0254          // cout<< "Q_ratio = " << Q_ratio[nruns][ipmt]<< " \t" <<"Q_current in each ch = " << Q_current[nruns][ipmt] << "\t" << "Q_ref in each ch = "<<Q_ref[ipmt] << " \t" << name <<"\t"<< nruns<< endl;
0255
0256           //Fill
0257           gt_chBych->SetPoint(nruns*128+ipmt,(Double_t)runSeconds, deltaT_all[nruns][ipmt]);
0258
0259           gq_chBych->SetPoint(nruns*128+ipmt,(Double_t)runSeconds, Q_ratio[nruns][ipmt] );
0260
0261           refFile_all->Close(); // I had to add it avoiding error like too many files open
0262           delete refFile_all;
0263
0264
0265         }
0266
0267       }
0268
0269
0270     }
0271     else
0272     {
0273       cout << "Histogram 'h2_mbdtt' not found in file: " << dstfname << endl;
0274     }
0275
0276     nruns++;
0277     icv++;
0278
0279   }
0280   if(tfile[nruns] != nullptr) // to avoid such error like opening too many files
0281   {
0282     tfile[nruns]->Close();
0283     delete tfile[nruns];
0284   }
0285
0286
0287   // Create the TGraph
0288   g_Canvas = new TCanvas("deltaT_graph", "Delta-T vs Time", 1000, 500);
0289   g_Canvas->cd();
0290   g_Canvas->Update();
0291   graph->SetTitle("Delta-T vs Date");
0292   graph->GetXaxis()->SetTimeDisplay(1);
0293   //graph->GetXaxis()->SetLabelSize(0.02);
0294  // graph->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%Y}");
0295   graph->GetXaxis()->SetTimeFormat("%m/%d/%Y ");
0296  // graph->GetXaxis()->SetLabelOffset(0.15);
0297   graph->GetXaxis()->SetLabelSize(0.02);
0298   graph->GetXaxis()->SetTimeOffset(0,"gmt");
0299   graph->GetXaxis()->SetTitle("Date");
0300   graph->GetYaxis()->SetTitle("Delta-T[ns]");
0301   graph->SetMarkerSize(0.6);
0302   graph->SetMarkerStyle(20);
0303   graph->SetMarkerColor(2.);
0304   graph->Draw("AP");
0305   g_Canvas->SaveAs("Delta-T_vs_Time.png");
0306
0307
0308
0309   gt_Canvas= new TCanvas("deltaT_graph_128ch", "Delta-T vs Date for 128 ch ", 1000, 500);
0310   gt_chBych->SetTitle("Delta-T vs Date");
0311   gt_chBych->GetXaxis()->SetTitle("Date");
0312   gt_chBych->GetYaxis()->SetTitle("Delta-T[ns]");
0313   gt_chBych->GetXaxis()->SetTimeDisplay(1);
0314   gt_chBych->GetXaxis()->SetTimeFormat("%m/%d/%Y ");
0315   gt_chBych->GetXaxis()->SetLabelSize(0.02);
0316   gt_chBych->SetFillColor(kBlue);
0317   // add space between the x-axis labels to move it vartcal
0318   //gt_graph->GetXaxis()->SetLabelOffset(0.08);
0319   gt_chBych->GetXaxis()->SetTimeOffset(0,"gmt");
0320   gt_chBych->SetMarkerSize(0.2);
0321   gt_chBych->SetMarkerStyle(5);
0322   gt_chBych->SetMarkerColor(2.);
0323   gt_chBych->Draw("AP");
0324   gt_Canvas->SaveAs("Delta-T_vs_Time128ch.png");
0325
0326
0327
0328
0329   gq_Canvas= new TCanvas("Q_ratio_graph_128ch", "Q_ratio vs Date for 128 ", 1000, 500);
0330   gq_chBych->SetTitle("Q_ratio vs Date");
0331   gq_chBych->GetXaxis()->SetTitle("Date");
0332   gq_chBych->GetYaxis()->SetTitle("Q_ratio[ADC]");
0333   gq_chBych->GetXaxis()->SetTimeDisplay(1);
0334   gq_chBych->GetXaxis()->SetTimeFormat("%Y/%m/%d ");
0335   gq_chBych->GetXaxis()->SetLabelSize(0.02);
0336   gq_chBych->GetXaxis()->SetTimeOffset(0,"gmt");
0337   gq_chBych->SetMarkerSize(0.2);
0338   gq_chBych->SetMarkerStyle(5);
0339   gq_chBych->SetMarkerColor(kBlue);
0340   gq_chBych->Draw("AP");
0341   gq_Canvas->SaveAs("Q_ration_vs_Time128ch.png");
0342
0343
0344
0345   gPad->Modified();
0346   savefile->Write();
0347   savefile->Close();
0348
0349   cout << "Processed " << nruns << " runs." << endl;
0350 }
0351
0352 int plot_deltT_vs_Timestamp()
0353 {  Double_t run_number;
0354   TDatime dt;
0355   get_time(run_number);
0356   plot_deltT_vs_time(dt);
0357   return 0;
0358 }