tpc/TpcDVCalib/drift.C

0001 /*
0002  * This macro shows a minimum working example of running the tracking
0003  * hit unpackers with some basic seeding algorithms to try to put together
0004  * tracks. There are some analysis modules run at the end which package
0005  * hits, clusters, and clusters on tracks into trees for analysis.
0006  */
0007
0008 #include <GlobalVariables.C>
0009
0010 #include <Trkr_TpcReadoutInit.C>
0011 #include <algorithm>
0012
0013 #include <cdbobjects/CDBTTree.h>
0014
0015 #include <ffamodules/CDBInterface.h>
0016
0017 #include <fun4all/Fun4AllUtils.h>
0018
0019 #include <phool/recoConsts.h>
0020
0021 #include <TCanvas.h>
0022 #include <TFile.h>
0023 #include <TGraph.h>
0024 #include <TH1.h>
0025 #include <TLatex.h>
0026 #include <TLine.h>
0027 #include <TROOT.h>
0028 #include <TStyle.h>
0029 #include <TSystem.h>
0030
0031 R__LOAD_LIBRARY(libtpc.so)
0032 R__LOAD_LIBRARY(libphool.so)
0033 R__LOAD_LIBRARY(libcdbobjects.so)
0034
0035 double driftGET(int runnumber = 51279);
0036 bool CheckHV(int irun = 51279);
0037
0038 //------------------------------------------------------------------------
0039 //  get drift velocity from cdb...  wjl aug 2024
0040 //
0041 void drift()
0042 {
0043   //    ~70 um/ns -> 0.007 in sphenix units
0044
0045   TGraph *gdriftrun = new TGraph();
0046   gdriftrun->SetMarkerStyle(20);
0047   gdriftrun->SetMarkerSize(0.9);
0048   gdriftrun->SetMarkerColor(kBlue);
0049   gdriftrun->SetLineColor(kBlue);
0050
0051   double drifts[2000] = {0};
0052   const int NPHASE = 9;
0053   int iphase = 0;
0054   const std::string phasename[NPHASE] = {"6x6", "28x28", "56x56", "111x111 8/9", "111x111 8/11", "", "BP#downarrow", "", ""};
0055   int lastruns[NPHASE] = {49722, 50025, 50458, 50936, 51109, 51191, 51301, 51495, 51617};
0056   double phaseN[NPHASE] = {0};
0057   double phaseA[NPHASE] = {0};
0058   //    double phaseNhv[NPHASE] = {0};
0059   //    double phaseAhv[NPHASE] = {0};
0060   double lowest[NPHASE] = {0};
0061   for (double &iph : lowest)
0062   {
0063     iph = 99.;
0064   }
0065
0066   //---- loop over run numbers...
0067   //
0068   for (int i = 0; i < 2000; i++)
0069   {
0070     //
0071     int irun = 49700 + i;
0072     double val = driftGET(irun);
0073     //
0074     if (irun <= 49722)
0075     {
0076       iphase = 0;
0077     }
0078     else if (irun <= 50025)
0079     {
0080       iphase = 1;
0081     }
0082     else if (irun <= 50458)
0083     {
0084       iphase = 2;
0085     }
0086     else if (irun <= 50936)
0087     {
0088       iphase = 3;
0089     }
0090     else if (irun <= 51109)
0091     {
0092       iphase = 4;
0093     }
0094     else if (irun <= 51191)
0095     {
0096       iphase = 5;
0097     }
0098     else if (irun <= 51301)
0099     {
0100       iphase = 6;
0101     }
0102     else if (irun <= 51495)
0103     {
0104       iphase = 7;
0105     }
0106     else
0107     {
0108       iphase = 8;
0109     }
0110     //
0111     //        bool CorrectHV    = CheckHV(irun);
0112     //
0113     drifts[i] = val;
0114     if (val > 0.)
0115     {
0116       gdriftrun->SetPoint(gdriftrun->GetN(), irun, val);
0117       phaseA[iphase] += val;
0118       phaseN[iphase] += 1.0;
0119       //            if (CorrectHV){
0120       //                phaseAhv[iphase]    += val;
0121       //                phaseNhv[iphase]    += 1.0;
0122       //            }
0123       lowest[iphase] = std::min(val, lowest[iphase]);
0124     }
0125   }
0126
0127   int nrunstot = 0;
0128   //    int nrunstothv  = 0;
0129   for (int iph = 0; iph < NPHASE; iph++)
0130   {
0131     if (phaseN[iph] > 0)
0132     {
0133       nrunstot += phaseN[iph];
0134       //            nrunstothv      += phaseNhv[iph];
0135       phaseA[iph] /= phaseN[iph];
0136       //            phaseAhv[iph]   /= phaseNhv[iph];
0137       //            std::cout<<"RunNum<="<<lastruns[iph]<<"\t Nruns="<<phaseN[iph]<<" "<<phaseNhv[iph]<<"\t <v>="<<phaseA[iph]<<std::endl;
0138     }
0139   }
0140   std::cout << "Nruns with custom <v> found = " << nrunstot << std::endl;
0141
0142   //---- paint setup
0143   int ican = -1;
0144   int itext = -1;
0145   //    int iline=-1;
0146   //    int ivline=-1;
0147   //    int iframe=-1;
0148   TCanvas *ccan[100];
0149   //    TH1F    *frame[100];
0150   TLatex *text[1000];
0151   TLine *line[100];
0152   TLine *vline[100];
0153   for (int i = 0; i < 1000; i++)
0154   {
0155     text[i] = new TLatex();
0156     text[i]->SetNDC();
0157     text[i]->SetTextFont(42);
0158     text[i]->SetTextSize(0.05);
0159     text[i]->SetTextAlign(33);
0160     line[i] = new TLine();
0161     line[i]->SetLineColor(1);
0162     line[i]->SetLineStyle(2);
0163     line[i]->SetLineWidth(2);
0164     vline[i] = new TLine(0, -1, 0, 1);
0165     vline[i]->SetLineColor(17);
0166     vline[i]->SetLineWidth(2);
0167   }
0168   gROOT->SetStyle("Modern");
0169   gStyle->SetOptStat(0);
0170   gStyle->SetTitleAlign(13);
0171   gStyle->SetTitleFontSize(0.08);
0172   gStyle->SetTitleX(0.20);
0173   gStyle->SetTitleY(0.98);
0174   gStyle->SetPadRightMargin(0.01);
0175   gStyle->SetPadTopMargin(0.01);
0176   gStyle->SetPadBottomMargin(0.08);
0177   gStyle->SetPadLeftMargin(0.12);
0178   gStyle->SetTitleSize(0.03, "xyzt");
0179   gStyle->SetLabelSize(0.03, "xyzt");
0180   //
0181   //---- end paint setup..
0182
0183   //---- paint...
0184   //
0185   ++ican;
0186   std::string hname = "ccan" + std::to_string(ican);
0187   ccan[ican] = new TCanvas(hname.c_str(), hname.c_str(), 20 + ican * 30, 30 + ican * 30, 0.7 * 1100, 0.7 * 850);
0188   ccan[ican]->cd();
0189   ccan[ican]->Divide(1, 1, 0.0001, 0.0001);
0190   ccan[ican]->cd(1);
0191   gdriftrun->Draw("AP");
0192   for (int iph = 0; iph < NPHASE; iph++)
0193   {
0194     ++itext;
0195     text[itext]->SetNDC(kFALSE);
0196     text[itext]->SetTextAlign(23);
0197     text[itext]->SetTextSize(0.03);
0198
0199     text[itext]->DrawLatex(lastruns[iph] - 20, 0.996 * lowest[iph], Form("%.5f", phaseA[iph]));  // NOLINT(hicpp-vararg)
0200     //++itext; text[itext]->SetNDC(kFALSE); text[itext]->SetTextAlign(23); text[itext]->SetTextSize(0.03); text[itext]->SetTextColor(kGreen+1);
0201     //  text[itext]->DrawLatex(lastruns[iph],0.996*lowest[iph]-0.00005,Form("%.5f",phaseAhv[iph]));
0202     //
0203     ++itext;
0204     text[itext]->SetNDC(kFALSE);
0205     text[itext]->SetTextAlign(21);
0206     text[itext]->SetTextSize(0.03);
0207     text[itext]->SetTextAngle(90);
0208     text[itext]->DrawLatex(lastruns[iph], 0.0067, phasename[iph].c_str());
0209   }
0210   ccan[ican]->cd();
0211   ccan[ican]->Update();
0212   ccan[ican]->Print("drift.ps(");
0213
0214   //---- close-out...
0215   //
0216   ccan[ican]->Print("drift.ps]");
0217   int iSuccess = gSystem->Exec("ps2pdf drift.ps drift.pdf");
0218   std::cout << "isuccess=" << iSuccess << std::endl;
0219   if (iSuccess == -1)
0220   {
0221     std::cout << "drift.pdf produced, deleting drift.ps file..." << std::endl;
0222     TString sexec = TString("/bin/rm drift.ps");
0223     std::cout << sexec.Data() << std::endl;
0224     gSystem->Exec(sexec.Data());
0225   }
0226   std::cout << "Done..." << std::endl;
0227
0228   std::cout << "Writing drift.root..." << std::endl;
0229   TFile *fout = new TFile("drift.root", "RECREATE");
0230   fout->cd();
0231   gdriftrun->Write("gdriftrun");
0232   fout->Close();
0233 }
0234
0235 //------------------------------------------------------------
0236 //
0237 //  return true if HV=3.3 kV
0238 //
0239 bool CheckHV(const int irun)
0240 {
0241   if (irun >= 49706 && irun <= 49713)
0242   {
0243     return true;
0244   }
0245   if (irun >= 50006 && irun <= 50014)
0246   {
0247     return true;
0248   }
0249   if (irun >= 50021 && irun <= 50025)
0250   {
0251     return true;
0252   }
0253   if (irun >= 50438 && irun <= 50440)
0254   {
0255     return true;
0256   }
0257   if (irun >= 50444 && irun <= 50449)
0258   {
0259     return true;
0260   }
0261   if (irun >= 50916 && irun <= 50921)
0262   {
0263     return true;
0264   }
0265   if (irun >= 50933 && irun <= 50936)
0266   {
0267     return true;
0268   }
0269   if (irun >= 51099 && irun <= 51107)
0270   {
0271     return true;
0272   }
0273   //    if (irun>=&&irun<=){ return true; } else
0274   //    if (irun>=&&irun<=){ return true; } else
0275   //    if (irun>=&&irun<=){ return true; } else
0276   //    if (irun>=&&irun<=){ return true; } else
0277   return false;
0278 }
0279
0280 //-----------------------------------------------------------
0281 //
0282 double driftGET(const int runnumber)
0283 {
0284   //
0285   TpcReadoutInit(runnumber);
0286   double driftOLD = G4TPC::tpc_drift_velocity_reco;
0287   //
0288   auto *rc = recoConsts::instance();
0289   rc->set_IntFlag("RUNNUMBER", runnumber);
0290   //
0291   Enable::CDB = true;
0292   rc->set_StringFlag("CDB_GLOBALTAG", "ProdA_2024");
0293   rc->set_uint64Flag("TIMESTAMP", 6);
0294   //
0295   // G4TPC::tpc_drift_velocity_reco = (8.0 / 1000) * 107.0 / 105.0;
0296   //
0297   //---- get drift velocity
0298   double driftCDB = 0;
0299   rc->set_uint64Flag("TIMESTAMP", runnumber);
0300   CDBInterface *cdb = CDBInterface::instance();
0301   std::string tpc_dv_calib_dir = cdb->getUrl("TPC_DRIFT_VELOCITY");
0302   if (!tpc_dv_calib_dir.empty())
0303   {  // path to driftvelocity calib file was found...
0304     // std::cout << "DRIFT: Path of tpc_drift_velocity_calib_file = "<< tpc_dv_calib_dir << std::endl;
0305     CDBTTree *cdbttree = new CDBTTree(tpc_dv_calib_dir);
0306     cdbttree->LoadCalibrations();
0307     driftCDB = cdbttree->GetSingleFloatValue("tpc_drift_velocity");
0308     std::cout << "DRIFT: RunNum=" << runnumber << "\t driftOLD=" << driftOLD << "\t driftCDB=" << driftCDB << std::endl;
0309     // if (driftCDB>0.&&driftCDB<0.01){
0310     //  G4TPC::tpc_drift_velocity_reco  = driftCDB;
0311     // }
0312   }
0313   else
0314   {
0315     // std::cout << "DRIFT: Path to tpc_drift_velocity_calib_file NOT FOUND...."<< std::endl;
0316   }
0317   //
0318   return driftCDB;
0319 }