Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages
MiniBooNEAna

CountPot.cxx

Go to the documentation of this file.
00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 #include "TStyle.h"
00009 #include "TFile.h"
00010 #include "TMath.h"
00011 #include "MiniBooNEAna/CountPot.h"
00012 #include "MiniBooNEAna/MBSpill.h"
00013 #include "MessageService/MsgService.h"
00014 #include "MinosObjectMap/MomNavigator.h"
00015 #include "RawData/RawRecord.h"
00016 #include "RawData/RawHeader.h"
00017 #include "RawData/RawDaqHeaderBlock.h"
00018 #include "RawData/RawDaqSnarlHeader.h"
00019 #include "Record/RecCandHeader.h"
00020 #include "CandData/CandHeader.h"
00021 #include "DatabaseInterface/DbiWriter.h"
00022 #include "DatabaseInterface/DbiResultPtr.h"
00023 #include "DatabaseInterface/DbiTableProxy.h"
00024 #include "DatabaseInterface/DbiDBProxy.h"
00025 #include "DatabaseInterface/DbiCache.h"
00026 #include "DatabaseInterface/DbiSqlContext.h"
00027 #include "Validity/VldContext.h"
00028 #include "DataUtil/GetCandHeader.h"
00029 #include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
00030 
00031 JOBMODULE(CountPot, "CountPot",
00032           "a module to count MiniBooNE pots based on Minos ND physics runs");
00033 CVSID("$Id: CountPot.cxx,v 1.7 2007/11/11 07:01:52 rhatcher Exp $");
00034 //......................................................................
00035 
00036 CountPot::CountPot()  
00037   : fTimeGap(0.0)
00038   , fTimeBuffer(0.0)
00039   , fPOTcut(0.0)
00040   , fhIlo(0.0)
00041   , fhIhi(0.0)
00042   , year(0)
00043   , month(0)
00044   , runtype(-1)
00045   , runtypepre(-1)
00046   , totpotmb(0.0)
00047   , totpot(0.0)
00048   , totpotphysics(0.0)
00049   , totpotphysicst(0.0)
00050   , totpotphysicsm(0.0)
00051   , totpotphysicstm(0.0)
00052   , totpotunknown(0.0)
00053   , pot(0.0)
00054   , potphys(0)
00055   , potphys_t(0)
00056   , potmb(0)
00057   , hImb(0)
00058   , meanhI(0)
00059   , potsum(0)
00060   , fFile(0)
00061 {
00062 }
00063 //......................................................................
00064 
00065 CountPot::~CountPot()
00066 {
00067 }
00068 
00069 //......................................................................
00070 
00071 void CountPot::BeginJob()
00072 {
00073   //gStyle->SetTimeOffset(0);
00074   timestart = VldTimeStamp(2020,1,1,0,0,0);
00075   timeend = timestart;
00076   timestart_month = timestart;
00077   timeend_month = VldTimeStamp(1990,1,1,0,0,0);
00078   yesterday = timestart;
00079   meanhI = 0;
00080   potsum = 0;
00081 }
00082 
00083 //......................................................................
00084 
00085 void CountPot::EndJob()
00086 {
00087   if (timestart!=timeend){//ignore the single pulse
00088     if (timestart<timestart_month) timestart_month = timestart;
00089     if (timeend>timeend_month) timeend_month = timeend;
00090     pot = 0.;
00091     DbiResultPtr<MBSpill> fBegin, fEnd;
00092     VldContext vcbeg(Detector::kNear,SimFlag::kData,timestart);
00093     VldContext vcend(Detector::kNear,SimFlag::kData,timeend);
00094     fBegin.NewQuery(vcbeg,0,true);
00095     fEnd.NewQuery(vcend,0,true);
00096     const DbiValidityRec* vrbeg = fBegin.GetValidityRec();
00097     const DbiValidityRec* vrend = fEnd.GetValidityRec();
00098     
00099     VldTimeStamp tsStart = vrbeg->GetVldRange().GetTimeStart();
00100     VldTimeStamp tsEnd = vrend->GetVldRange().GetTimeEnd();
00101     //    tsStart=tsStart-VldTimeStamp(fTimeBuffer);
00102     //    tsEnd.Add(fTimeBuffer);
00103     DbiSqlContext extContext(DbiSqlContext::kStarts,
00104                              tsStart,tsEnd,
00105                              Detector::kNear,SimFlag::kData);     
00106     DbiResultPtr<MBSpill> rsMBSpill("MBSPILL",extContext,Dbi::kAnyTask);  
00107     const int numrows = rsMBSpill.GetNumRows();
00108     for (int ir = 0; ir<numrows; ir++){
00109       const MBSpill *rs = rsMBSpill.GetRow(ir);
00110       if (double(rs->GetTimeStamp()-timestart)>-fTimeBuffer
00111           &&double(rs->GetTimeStamp()-timeend)<fTimeBuffer){
00112         if(rs->GetPOT()>fPOTcut
00113            &&rs->GetHorn_current()>fhIlo
00114            &&rs->GetHorn_current()<fhIhi){
00115           pot += rs->GetPOT();
00116           if (runtypepre==0||runtypepre==2){//test run
00117             potphys_t->Fill(double(rs->GetTimeStamp()),rs->GetPOT());
00118           }
00119           else if (runtypepre==1||runtypepre==3){//physics run
00120             potphys->Fill(double(rs->GetTimeStamp()),rs->GetPOT());
00121           }
00122         }
00123         if (rs->GetPOT()>fPOTcut){
00124           if (yesterday>rs->GetTimeStamp()){//start counting
00125             yesterday = rs->GetTimeStamp();
00126             meanhI = rs->GetHorn_current()*rs->GetPOT();
00127             potsum = rs->GetPOT();
00128           }
00129           else {
00130             UInt_t y1,m1,d1;
00131             UInt_t y2,m2,d2;
00132             yesterday.GetDate(true,0,&y1,&m1,&d1);
00133             rs->GetTimeStamp().GetDate(true,0,&y2,&m2,&d2);
00134             if (d1!=d2){//tomorrow is another day
00135               if (potsum>0){
00136                 meanhI/=potsum;
00137                 hImb->Fill(double(yesterday),meanhI);
00138               }
00139               yesterday = rs->GetTimeStamp();
00140               meanhI = rs->GetHorn_current()*rs->GetPOT();
00141               potsum = rs->GetPOT();                                
00142             }
00143             else {
00144               meanhI += rs->GetHorn_current()*rs->GetPOT();
00145               potsum += rs->GetPOT();
00146             }
00147           }
00148         }
00149       }
00150     }
00151     totpot += pot;
00152     switch(runtypepre){
00153     case 0:
00154       totpotphysicst += pot;
00155       break;
00156     case 1:
00157       totpotphysicsm += pot;
00158       break;
00159     case 2:
00160       totpotphysicstm += pot;
00161       break;
00162     case 3:
00163       totpotphysics += pot;
00164       break;
00165     case 4:
00166       totpotunknown += pot;
00167       break;
00168     }
00169     cout<<"POTs between "<<timestart<<"("<<begrun<<"/"<<begsubrun<<"/"<<begsnarl<<")"<<" and "<<timeend<<"("<<endrun<<"/"<<endsubrun<<"/"<<endsnarl<<")"<<"(runtype:"<<runtypepre<<")"<<":"<<pot<<" Tot:"<<totpot<<endl;
00170   }
00171   if (potsum>0){
00172     meanhI/=potsum;
00173     hImb->Fill(double(yesterday),meanhI);
00174   }
00175 
00176   FillMBPOT();
00177 
00178   char ch[1000];
00179   sprintf(ch,"POTs vs Time(%d%d/%d%d) Total: %.3e; ND on: %.3e; test run: %.3e",month/10,month%10,year%100/10,year%100%10,totpotmb*1e12,(totpotphysics+totpotphysicsm)*1e12,(totpotphysicst+totpotphysicstm)*1e12);
00180   potphys->SetTitle(ch);
00181 
00182   cout<<"*******************************************************"<<endl;
00183   cout<<"physics: "<<totpotphysics<<endl;
00184   cout<<"physics;t: "<<totpotphysicst<<endl;
00185   cout<<"physics;m: "<<totpotphysicsm<<endl;
00186   cout<<"physics;tm: "<<totpotphysicstm<<endl;
00187   cout<<"unknown: "<<totpotunknown<<endl;
00188   cout<<"Tot: "<<totpot<<endl;
00189   cout<<"*******************************************************"<<endl;
00190 
00191   if (fFile){
00192     fFile->Write();
00193     fFile->Close();
00194   }
00195 }
00196 //......................................................................
00197 
00198 JobCResult CountPot::Ana(const MomNavigator* mom)
00199 {
00200   int run, subrun, snarl;
00201 
00202   TIter next(mom->GetFragmentArray());
00203   TObject* frag;
00204   while((frag = next())){
00205     if(frag->InheritsFrom("RawRecord")) {
00206       
00207       RawRecord* rawrec = dynamic_cast<RawRecord*>(frag);
00208       
00209       if (!rawrec) { 
00210         MSG("CountPot", Msg::kWarning) << "No RawRecord in MOM." << endl; 
00211         return JobCResult::kFailed;
00212       }
00213 
00214       const RawSnarlHeaderBlock* rshb = 
00215         dynamic_cast<const RawSnarlHeaderBlock*>
00216         (rawrec->FindRawBlock("RawSnarlHeaderBlock"));
00217       if (!rshb) {
00218         MSG("CountPot", Msg::kWarning) << "No RawSnarlHeaderBlock." << endl; 
00219         return JobCResult::kFailed;
00220       }
00221       if ((rshb->GetRunType()&0xfff)!=0x301){//not a physics run(0x301)
00222         MSG("CountPot", Msg::kWarning) << "Not a physics run." <<endl; 
00223         return JobCResult::kPassed;     
00224       }
00225       if ((rshb->GetRunType()&0xf000&0x8000)&&!(rshb->GetRunType()&0xf000&0x4000)){//physics;t
00226         runtype = 0;
00227       }
00228       else if (!(rshb->GetRunType()&0xf000&0x8000)&&(rshb->GetRunType()&0xf000&0x4000)){//physics;m
00229         runtype = 1;
00230       }
00231       else if ((rshb->GetRunType()&0xf000&0x8000)&&(rshb->GetRunType()&0xf000&0x4000)){//physics;tm
00232         runtype = 2;
00233       }
00234       else if (!(rshb->GetRunType()&0xf000&0x8000)&&!(rshb->GetRunType()&0xf000&0x4000)){//physics
00235         runtype = 3;
00236       }
00237       else runtype = 4;
00238 
00239       if (rshb->GetTriggerSource()&4==0&&rshb->GetTriggerSource()&16==0){
00240         MSG("CountPot", Msg::kDebug) << "Not passed plane trigger or activity trigger" << endl; 
00241         return JobCResult::kPassed;
00242       }
00243 
00244       run = rshb->GetRun();
00245       subrun = rshb->GetSubRun();
00246       snarl = rshb->GetSnarl();
00247 
00248       static bool first = true;
00249       VldTimeStamp snarltime = rshb->GetTriggerTime();
00250       if (first) {
00251         timestart = snarltime;
00252         timeend = snarltime;
00253         begrun = run;
00254         endrun = run;
00255         begsubrun = subrun;
00256         endsubrun = subrun;
00257         begsnarl = snarl;
00258         endsnarl = snarl;
00259         runtypepre = runtype;
00260         //find out which year and which month we are looking at
00261         UInt_t day;
00262         snarltime.GetDate(true,0,&year,&month,&day);
00263         if (day>20) {//It's very likely we are looking at the last day of the previous month
00264           if (month==12){//happy new year
00265             year++;
00266             month=1;
00267           }
00268           else month++;
00269         }
00270         this->bookhist();
00271         first = false;
00272       }
00273       else if (double(snarltime-timeend)>fTimeGap){//find a big time gap
00274         if (timestart==timeend){//ignore the single pulse
00275           timestart = snarltime;
00276           timeend = snarltime;
00277           begrun = run;
00278           endrun = run;
00279           begsubrun = subrun;
00280           endsubrun = subrun;
00281           begsnarl = snarl;
00282           endsnarl = snarl;
00283         }
00284         else {//calculate the POTs
00285           if (timestart<timestart_month) timestart_month = timestart;
00286           if (timeend>timeend_month) timeend_month = timeend;
00287           pot = 0.;
00288           DbiResultPtr<MBSpill> fBegin, fEnd;
00289           VldContext vcbeg(Detector::kNear,SimFlag::kData,timestart);
00290           VldContext vcend(Detector::kNear,SimFlag::kData,timeend);
00291           fBegin.NewQuery(vcbeg,0,true);
00292           fEnd.NewQuery(vcend,0,true);
00293           const DbiValidityRec* vrbeg = fBegin.GetValidityRec();
00294           const DbiValidityRec* vrend = fEnd.GetValidityRec();
00295 
00296           VldTimeStamp tsStart = vrbeg->GetVldRange().GetTimeStart();
00297           VldTimeStamp tsEnd = vrend->GetVldRange().GetTimeEnd();
00298 //        tsStart=tsStart-VldTimeStamp(fTimeBuffer);
00299 //        tsEnd.Add(fTimeBuffer);
00300           DbiSqlContext extContext(DbiSqlContext::kStarts,
00301                                    tsStart,tsEnd,
00302                                    Detector::kNear,SimFlag::kData);       
00303           DbiResultPtr<MBSpill> rsMBSpill("MBSPILL",extContext,Dbi::kAnyTask);  
00304           const int numrows = rsMBSpill.GetNumRows();
00305           for (int ir = 0; ir<numrows; ir++){
00306             const MBSpill *rs = rsMBSpill.GetRow(ir);
00307             if (double(rs->GetTimeStamp()-timestart)>-fTimeBuffer
00308                 &&double(rs->GetTimeStamp()-timeend)<fTimeBuffer){
00309               if(rs->GetPOT()>fPOTcut
00310                  &&TMath::Abs(rs->GetHorn_current())>fhIlo
00311                  &&TMath::Abs(rs->GetHorn_current())<fhIhi){
00312                 pot += rs->GetPOT();
00313                 if (runtypepre==0||runtypepre==2){//test run
00314                   potphys_t->Fill(double(rs->GetTimeStamp()),rs->GetPOT());
00315                 }
00316                 else if (runtypepre==1||runtypepre==3){//physics run
00317                   potphys->Fill(double(rs->GetTimeStamp()),rs->GetPOT());
00318                 }
00319               }
00320               if (rs->GetPOT()>fPOTcut){
00321                 if (yesterday>rs->GetTimeStamp()){//start counting
00322                   yesterday = rs->GetTimeStamp();
00323                   meanhI = rs->GetHorn_current()*rs->GetPOT();
00324                   potsum = rs->GetPOT();
00325                 }
00326                 else {
00327                   UInt_t y1,m1,d1;
00328                   UInt_t y2,m2,d2;
00329                   yesterday.GetDate(true,0,&y1,&m1,&d1);
00330                   rs->GetTimeStamp().GetDate(true,0,&y2,&m2,&d2);
00331                   if (d1!=d2){//tomorrow is another day
00332                     if (potsum>0){
00333                       meanhI/=potsum;
00334                       hImb->Fill(double(yesterday),meanhI);
00335                     }
00336                     yesterday = rs->GetTimeStamp();
00337                     meanhI = rs->GetHorn_current()*rs->GetPOT();
00338                     potsum = rs->GetPOT();                                  
00339                   }
00340                   else {
00341                     meanhI += rs->GetHorn_current()*rs->GetPOT();
00342                     potsum += rs->GetPOT();
00343                   }
00344                 }
00345               }
00346             }
00347           }
00348           totpot += pot;
00349           switch(runtypepre){
00350           case 0:
00351             totpotphysicst += pot;
00352             break;
00353           case 1:
00354             totpotphysicsm += pot;
00355             break;
00356           case 2:
00357             totpotphysicstm += pot;
00358             break;
00359           case 3:
00360             totpotphysics += pot;
00361             break;
00362           case 4:
00363             totpotunknown += pot;
00364             break;
00365           }
00366           cout<<"POTs between "<<timestart<<"("<<begrun<<"/"<<begsubrun<<"/"<<begsnarl<<")"<<" and "<<timeend<<"("<<endrun<<"/"<<endsubrun<<"/"<<endsnarl<<")"<<"(runtype:"<<runtypepre<<")"<<":"<<pot<<" Tot:"<<totpot<<endl;
00367           timestart = snarltime;
00368           timeend = snarltime;
00369           begrun = run;
00370           endrun = run;
00371           begsubrun = subrun;
00372           endsubrun = subrun;
00373           begsnarl = snarl;
00374           endsnarl = snarl;
00375           runtypepre = runtype;
00376         }
00377       }//find a big time gap
00378       else {
00379         timeend = snarltime;
00380         endrun = run;
00381         endsubrun = subrun;
00382         endsnarl = snarl;
00383       } 
00384     }
00385   }
00386 
00387   return JobCResult::kPassed; // kNoDecision, kFailed, etc.
00388 }
00389 
00390 //......................................................................
00391 
00392 void CountPot::bookhist(){
00393 
00394   VldTimeStamp t1(year,month,1,0,0,0);
00395   VldTimeStamp t2;
00396   if (month==12){
00397     t2 = VldTimeStamp(year+1,1,1,0,0,0);
00398   }
00399   else t2 = VldTimeStamp(year,month+1,1,0,0,0);
00400   fFile = new TFile(Form("beaminfo%d%d%d%d.root",month/10,month%10,year%100/10,year%100%10),"recreate");
00401   potphys = new TH1D("potphys","Number of POTs vs Time",int(double(t2-t1)/86400),(double)t1,(double)t2);
00402   potphys_t = new TH1D("potphys_t","Number of POTs vs Time",int(double(t2-t1)/86400),(double)t1,(double)t2);
00403   potmb = new TH1D("potmb","Number of POTs vs Time",int(double(t2-t1)/86400),(double)t1,(double)t2);
00404   hImb = new TH1D("hImb","horn current vs Time",int(double(t2-t1)/86400),(double)t1,(double)t2);
00405 
00406 }
00407 //......................................................................
00408 
00409 void CountPot::FillMBPOT(){
00410   if (timeend_month < timestart_month) {
00411     MSG("CountPot", Msg::kWarning) << "timeend<timestart, somthing is wrong" << endl; 
00412     return;
00413   }
00414   
00415   VldTimeStamp t1 = timestart_month;
00416   VldTimeStamp t2 = t1;
00417   t2.Add(VldTimeStamp(86400,0));
00418   if (t2>timeend_month) t2 = timeend_month;
00419   cout<<t1<<endl;
00420   cout<<t2<<endl;
00421   cout<<timestart_month<<endl;
00422   cout<<timeend_month<<endl;
00423   while (t1<timeend_month){
00424     DbiResultPtr<MBSpill> fBegin, fEnd;
00425     VldContext vcbeg(Detector::kNear,SimFlag::kData,t1);
00426     VldContext vcend(Detector::kNear,SimFlag::kData,t2);
00427     fBegin.NewQuery(vcbeg,0,true);
00428     fEnd.NewQuery(vcend,0,true);
00429     const DbiValidityRec* vrbeg = fBegin.GetValidityRec();
00430     const DbiValidityRec* vrend = fEnd.GetValidityRec();
00431     
00432     VldTimeStamp tsStart = vrbeg->GetVldRange().GetTimeStart();
00433     VldTimeStamp tsEnd = vrend->GetVldRange().GetTimeEnd();
00434     DbiSqlContext extContext(DbiSqlContext::kStarts,
00435                              tsStart,tsEnd,
00436                              Detector::kNear,SimFlag::kData);     
00437     DbiResultPtr<MBSpill> rsMBSpill("MBSPILL",extContext,Dbi::kAnyTask);  
00438     const int numrows = rsMBSpill.GetNumRows();
00439     for (int ir = 0; ir<numrows; ir++){
00440       const MBSpill *rs = rsMBSpill.GetRow(ir);
00441       if (double(rs->GetTimeStamp()-t1)>=0
00442           &&double(rs->GetTimeStamp()-t2)<=0){
00443         if(rs->GetPOT()>fPOTcut){
00444           potmb->Fill(double(rs->GetTimeStamp()),rs->GetPOT());
00445           totpotmb += rs->GetPOT();
00446         }
00447       }
00448     }
00449     t1=t2;
00450     t2.Add(VldTimeStamp(86400,0));
00451     if (t2>timeend_month) t2 = timeend_month;
00452   }
00453   cout<<"Total MiniBooNE POTs "<<totpotmb<<endl;
00454 }
00455 //......................................................................
00456 const Registry& CountPot::DefaultConfig() const
00457 {
00461   static Registry r; // Default configuration for module
00462 
00463   // Set name of config
00464   std::string name = this->GetName();
00465   name += ".config.default";
00466   r.SetName(name.c_str());
00467 
00468   // Set values in configuration
00469   r.UnLockValues();
00470   r.Set("TimeGap", 30.0);
00471   r.Set("TimeBuffer",0.1);
00472   r.Set("POTcut",0.1);
00473   r.Set("hIlo",173);
00474   r.Set("hIhi",175);
00475   r.LockValues();
00476 
00477   return r;
00478 }
00479 
00480 //......................................................................
00481 
00482 void CountPot::Config(const Registry& r)
00483 {
00487   fTimeGap = r.GetDouble("TimeGap");
00488   fTimeBuffer = r.GetDouble("TimeBuffer");
00489   fPOTcut = r.GetDouble("POTcut");
00490   fhIlo = r.GetDouble("hIlo");
00491   fhIhi = r.GetDouble("hIhi");
00492 }
00493
Generated on Mon Feb 15 11:06:33 2010 for loon by  doxygen 1.3.9.1