loon: LossModule.cxx Source File

00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 #include "LossModule.h"
00013 
00014 #include "StripHist.h"
00015 
00016 #include <JobControl/JobCModuleRegistry.h>
00017 #include <JobControl/JobCResult.h>
00018 
00019 #include <RawData/RawRecord.h>
00020 #include <RawData/RawBeamMonBlock.h>
00021 #include <RawData/RawBeamData.h>
00022 
00023 #include <DataUtil/GetRecords.h>
00024 #include <BeamDataUtil/BDDevices.h>
00025 
00026 #include <MessageService/MsgService.h>
00027 
00028 #include <Conventions/Munits.h>
00029 
00030 #include <TGraph.h>
00031 #include <TCanvas.h>
00032 #include <TH2F.h>
00033 #include <TLatex.h>
00034 #include <TStyle.h>
00035 #include <TGaxis.h>
00036 
00037 #include <vector>
00038 
00039 //ClassImp(LossModule)
00040 
00041 #define minimum(a,b) ((a) > (b) ? (b) : (a)) 
00042 #define maximum(a,b) ((a) > (b) ? (a) : (b)) 
00043 
00044 using namespace std;
00045 
00046 CVSID("$Id: LossModule.cxx,v 1.9 2005/06/09 15:26:43 thosieck Exp $");
00047 JOBMODULE(LossModule,"MonLoss","Proton loss related for Monitoring");
00048 
00049 LossModule::LossModule()
00050 {
00051  
00052 
00053   for (int bafn = 1; bafn <= 2; ++bafn) {
00054     
00055     StripHist* sh = new StripHist(Form("Baffle Scrape %d",bafn),
00056                                   Form("Estimate of beam loss on baffle, thermocouple %d",bafn),
00057                                   100,-3000,100);
00058     sh->SetStripRange(1*Munits::day);
00059     sh->GetHist().SetXTitle("(Baffle Temp - Air Temp)/Toroid (degC/1E12 protons)");
00060     fStripHist[Form("Baffle Scrape %d",bafn)] = sh;
00061   }
00062 
00063   TGraph *sh1 = new TGraph("BPM Location along Beam Line H=Red V=Blue","BPM Location along Beam Line H=Red V=Blue");
00064   fTGraph["Beam Position H"] = sh1;
00065   TGraph *sh2 = new TGraph("BPM Location along Beam Line H=Red V=Blue","BPM Location along Beam Line H=Red V=Blue");
00066   fTGraph["Beam Position V"] = sh2;
00067   
00068   sh1->GetXaxis()->SetTitle("Distance from Q806 (ft)");
00069   sh2->GetXaxis()->SetTitle("Distance from Q806 (ft)");
00070 
00071   TH1F *s1 = new TH1F("Beam Loss","Beam Loss along tunnel",100,-40,1200);
00072 
00073   s1->SetXTitle("Distance from Q806 (ft)");
00074   s1->SetYTitle("Beam Loss (Rad / s)");
00075   s1->GetXaxis()->CenterTitle();
00076   s1->GetYaxis()->CenterTitle();
00077   s1->GetYaxis()->SetTitleOffset(1.3);
00078   //TGaxis::SetMaxDigits(1);
00079 
00080   fTH1F["Beam Loss"] = s1;
00081 }
00082 
00083 LossModule::~LossModule()
00084 {
00085 
00086 }
00087 
00088 typedef std::map<std::string,StripHist*>  StripHistMap;
00089 typedef std::map<std::string,TH1F*> TH1FMap;
00090 
00091 void LossModule::BeginJob()
00092 {
00093   HistMan hm = this->GetHistMan();
00094 
00095   StripHist *sht1 = fStripHist["Baffle Scrape 1"];
00096   TCanvas* canvas = new TCanvas("Baffle Scrape 1",sht1->GetHist().GetTitle(),500,400);
00097   sht1->Draw("AL");
00098   hm.Adopt("Loss",canvas);
00099 
00100   StripHist *sht2 = fStripHist["Baffle Scrape 2"];
00101   TCanvas* canvas2 = new TCanvas("Baffle Scrape 2",sht2->GetHist().GetTitle(),500,400);
00102   sht2->Draw("AL");
00103   hm.Adopt("Loss",canvas2);
00104 
00105   TGraph *sht3 = fTGraph["Beam Position H"];
00106   TGraph *sht4 = fTGraph["Beam Position V"];
00107 
00108   sht3->SetLineColor(2);
00109   sht4->SetLineColor(4);
00110 
00111   TCanvas *canvas3 = new TCanvas("Beam Position","BPM Location along Beam Line",500,400);
00112   sht3->Draw("AL");
00113   sht4->Draw("LP");
00114   hm.Adopt("Loss",canvas3);
00115 
00116   TH1FMap::iterator mit2, done2=fTH1F.end();
00117   for(mit2=fTH1F.begin(); mit2 != done2; ++mit2) {
00118     TH1F *sh2 = mit2->second;
00119     TCanvas *canvas3 = new TCanvas(sh2->GetTitle(),sh2->GetTitle(),500,400);
00120     TVirtualPad* pad = gPad, *oldpad = gPad;    pad->cd();
00121     sh2->Draw();
00122     oldpad->cd();
00123     hm.Adopt("Loss",canvas3);
00124   }
00125 
00126   TH1F *th1 = fTH1F["Beam Loss"];
00127 
00128   TCanvas *canvas4 = new TCanvas("Proton Beam Transport Line Summary Canvas",
00129                                  "Proton Beam Transport Line Summary Canvas",800,800);
00130   canvas4->Divide(1,3);
00131   canvas4->cd(1);
00132   sht3->Draw("AL");
00133   sht4->Draw("LP");
00134   canvas4->cd(2);
00135   th1->Draw();
00136   canvas4->cd(3);
00137   
00138   hm.Adopt("Summary Canvases",canvas4);
00139 
00140 }
00141 
00142 void LossModule::Fill(const RawBeamMonHeaderBlock& /*head*/, const RawBeamMonBlock& block)
00143 {
00144   const RawBeamData* tpcast = block["E:TPCAST"];
00145   const RawBeamData* baft1  = block["E:BAFT1"];
00146   const RawBeamData* baft2  = block["E:BAFT2"];
00147   const RawBeamData* tortgt = block["E:TORTGT"];
00148   const RawBeamData *Etrtgtd = block["E:TRTGTD"];
00149 
00150   if(Etrtgtd && Etrtgtd->GetDataLength()) {
00151     double P2tgt  = Etrtgtd->GetData()[0];
00152     if(P2tgt < 0.1) {
00153       //If Beam intensity low, get out.
00154       return;
00155     }
00156   }
00157   
00158   if (!(tpcast && baft1 && baft2 && tortgt)) return;
00159   if (!(tpcast->GetDataLength() &&
00160         baft1->GetDataLength() &&
00161         baft2->GetDataLength() &&
00162         tortgt->GetDataLength()))
00163     return;
00164   
00165   StripHist* sh1 = fStripHist["Baffle Scrape 1"];
00166   StripHist* sh2 = fStripHist["Baffle Scrape 2"];
00167   
00168   double t = tpcast->GetData()[0];
00169   double pot = tortgt->GetData()[0];
00170   double y1 = Munits::ToCelcius(Munits::FromFahrenheit(baft1->GetData()[0] - t))/pot;
00171   double y2 = Munits::ToCelcius(Munits::FromFahrenheit(baft2->GetData()[0] - t))/pot;
00172   double dae = tortgt->GetSeconds() + 1.0e-6*tortgt->GetMsecs();
00173   
00174   sh1->Fill(dae,y1);
00175   sh2->Fill(dae,y2);
00176 
00177   
00178   float feet[40] = {-14.76, -2.02, 
00179                     27.55,   38.55,   49.57,   61.80,  72.80, 
00180                     83.80,   96.61,   108.84,  133.79, 
00181                     148.05,  204.60,  261.15,  317.78, 330.13, 
00182                     351.21,  374.22,  386.34,  407.42, 430.31, 
00183                     437.55,  458.63,  487.31,  502.01, 519.6, 
00184                     737.51,  752.20,  814.51,  843.55, 923.94, 
00185                     938.60,  945.05,  976.62,  997.70, 1020.59, 
00186                     1027.83, 1048.91, 1070.12, 1081.53};
00187 
00188   std::string lossdevice[40] = {"E:LMV100","E:LMQ101",
00189                                 "E:LM1011","E:LM1012","E:LMQ102","E:LM1013","E:LM1014",
00190                                 "E:LM1015","E:LMQ103","E:LM1016","E:LMQ104",
00191                                 "E:LMQ105","E:LMQ106","E:LMQ107","E:LMQ108","E:LM1081",
00192                                 "E:LM1082","E:LMQ109","E:LM1083","E:LM1084","E:LMQ110",
00193                                 "E:LM1085","E:LM1086","E:LMQ111","E:LMQ112","E:LMBBBP",
00194                                 "E:LMQ113","E:LMQ114","E:LMQ115","E:LMQ116","E:LMQ117",
00195                                 "E:LMH117","E:LMQ118","E:LM1181","E:LM1182","E:LMQ119",
00196                                 "E:LM1183","E:LM1184","E:LMQ120","E:LMQ121"};
00197   
00198   TH1F *sh3 = fTH1F["Beam Loss"];
00199   sh3->Reset();
00200   for(int i=0; i<40; i++) {
00201     const RawBeamData* tmp = block[lossdevice[i]];
00202     double blah = tmp->GetData()[0];
00203     sh3->Fill(feet[i], blah );
00204   }
00205   
00206   float feet2[23] = {16.79, 17.70, 60.27, 107.31, 132.58, 145.21, 203.39,
00207                      258.31, 314.94, 385.04, 436.18,  498.16, 512.88, 736.29,
00208                      762.90, 811.78, 842.33, 921.10, 943.83, 1095.81, 1096.89,
00209                      1134.89, 1135.80};
00210   
00211   TGraph *sh10 = fTGraph["Beam Position H"];
00212   TGraph *sh11 = fTGraph["Beam Position V"];
00213   //MSG("BD",Msg::kDebug) << "A) sh10 has " << sh10->GetN() << "entries" << endl;
00214   //MSG("BD",Msg::kDebug) << "A) sh11 has " << sh11->GetN() << "entries" << endl;
00215 
00216   for(int i = 0; i < 13; i++) { //sh10->GetN(); i++) {
00217     sh10->RemovePoint(0);
00218   }
00219   for(int j = 0; j < 10; j++) {//sh11->GetN(); j++) {
00220     sh11->RemovePoint(0);
00221   }
00222   
00223   //MSG("BD",Msg::kDebug) << "B) sh10 has " << sh10->GetN() << "entries" << endl;
00224   //MSG("BD",Msg::kDebug) << "B) sh11 has " << sh11->GetN() << "entries" << endl;
00225 
00226 
00228     std::string bpm_positions[23] = {
00229       "E:VP101",    "E:HP101",    "E:HP102",    "E:VP103",    "E:HP104",    
00230       "E:HP105",    "E:VP106",    "E:HP107",    "E:VP108",    "E:HP109",    
00231       "E:VP110",    "E:VP111",    "E:HP112",    "E:VP113",    "E:HP114", 
00232       "E:HP115",    "E:VP116",    "E:HP117",    "E:VP118",    "E:HP119",  
00233       "E:HP121",    "E:HPTGT",    "E:VPTGT"};
00234     
00235     for(int i=0; i<23; i++) {
00236       const RawBeamData* datum = block[bpm_positions[i]];
00237       if(datum && datum->GetDataLength()) {
00238         double datum2 = datum->GetData()[0];
00239         if(datum2 < 990) {
00240           /*
00241           Explanation for the above cut
00242           
00243           FYI, according to Peter Lucas:
00244           
00245           "We now understand the rationale for zero intensity leading to large
00246           values of position. This is intentional and is connected with the BPM
00247           averaging mask. when there is no beam in the masked-in batches, 1
00248           meter is returned as the position. - P. Lucas"
00249           */
00250           std::string tempor = bpm_positions[i];
00251           char ch = tempor[2];
00252           //MSG("BD",Msg::kDebug) << "ch=[" << ch << "]" << endl; 
00253           if(ch == 'V') {
00254             sh11->SetPoint(sh11->GetN(), feet2[i], datum2);
00255           } else {
00256             sh10->SetPoint(sh10->GetN(), feet2[i], datum2);
00257           }
00258         }
00259       }
00260     }
00261 
00262     //Since using TGraph's, have to do this every update.  blah.
00263     TH1F* frame = sh10->GetHistogram();
00264     TAxis* axis_x = frame->GetXaxis();
00265     axis_x->SetTitle("Distance from Q806 (ft)");
00266     axis_x->CenterTitle();
00267     TAxis *axis_y = frame->GetYaxis();
00268     axis_y->SetTitle("Beam Position (mm)");
00269     axis_y->CenterTitle();
00270     
00271     TH1F* frame2 = sh11->GetHistogram();
00272     TAxis* axis_x2 = frame2->GetXaxis();
00273     axis_x2->SetTitle("Distance from Q806 (ft)");
00274     axis_x2->CenterTitle();
00275     TAxis *axis_y2 = frame2->GetYaxis();
00276     axis_y2->SetTitle("Beam Position (mm)");
00277     axis_y2->CenterTitle();
00278 
00279     double min1 = 0.0, max1 = 0.0;
00280     RangeFinder(sh10, sh11, min1, max1);
00281     sh10->GetYaxis()->SetRangeUser(min1, max1);
00282     sh11->GetYaxis()->SetRangeUser(min1, max1);
00283 
00284     const RawBeamData* datum3 = block["I:BEAM"];
00285     const RawBeamData* datum4 = block["E:TR101D"];
00286     const RawBeamData* datum5 = block["E:TRTGTD"];
00287 
00288     double datum6 = -1.0, datum7 = -1.0, datum8 = -1.0;
00289 
00290     if(datum3 && datum3->GetDataLength()) {
00291       datum6 = datum3->GetData()[0];
00292     }
00293     if(datum4 && datum4->GetDataLength()) {
00294       datum7 = datum4->GetData()[0];
00295     }
00296     if(datum5 && datum5->GetDataLength()) {
00297       datum8 = datum5->GetData()[0];
00298     }
00299     
00300     HistMan hm = this->GetHistMan();
00301     TCanvas *can = (TCanvas*) hm.GetObject("Summary Canvases/Proton Beam Transport Line Summary Canvas");
00302     TVirtualPad *pad = can->cd(3);
00303     pad->Clear();
00304     TLatex l;
00305     l.SetTextSize(0.1);
00306     l.DrawLatex(0.1,0.7,Form("Beam From Main Injector = %f x10^{12} (I:BEAM)",datum6));
00307     l.DrawLatex(0.1,0.4,Form("Beam into NuMI Line = %f x10^{12} (E:TR101D)",datum7));
00308     l.DrawLatex(0.1,0.1,Form("Beam at NuMI Target = %f x10^{12} (E:TRTGTD)",datum8));
00309 }
00310 
00311 void LossModule::RangeFinder(TGraph *sh1, TGraph *sh2, double &min, double &max)
00312 {
00313   double max1 = GetMax(sh1), max2 = GetMax(sh2);
00314   double min1 = GetMin(sh1), min2 = GetMin(sh2);
00315 
00316   min = minimum(min1, min2);
00317   max = maximum(max1, max2);
00318 
00319   if(min != 0.0) {
00320     if(min < 0) {
00321       min *= 1.5;
00322     } else {
00323       min *= 0.5;
00324     }
00325   } else {
00326     min = -0.5;
00327   }
00328 
00329   if(max != 0.0) {
00330     if(max < 0) {
00331       max *= 0.5;
00332     } else {
00333       max *= 1.5;
00334     }
00335   } else {
00336     max = 0.5;
00337   }
00338 
00339 }
00340 
00341 double LossModule::GetMax(TGraph *gr)
00342 {
00343     if(gr->GetN() == 0) return 0.0;
00344     double max = 0.0, x = 0.0, y = 0.0;
00345     for(int i=0; i<gr->GetN(); i++) {
00346         gr->GetPoint(i, x, y);
00347         if(y > max) max = y;
00348     }
00349     return max;
00350 }
00351  
00352 double LossModule::GetMin(TGraph *gr)
00353 {
00354     if(gr->GetN() == 0) return 0.0;
00355     double max = GetMax(gr);
00356     double min = max, x = 0.0, y = 0.0;
00357     for(int i=0; i<gr->GetN(); i++) {
00358         gr->GetPoint(i, x, y);
00359         if(y < min) min = y;
00360     }
00361     return min;
00362 }