loon: NuBeam.cxx Source File

00001 
00002 
00003 // Coded by Jeff Hartnell Mar/2006       
00004 //                                                                    
00005 // Contact: j.j.hartnell@rl.ac.uk
00007 
00008 #include <cmath>
00009 
00010 #include "TH1F.h"
00011 #include "TMath.h"
00012 #include "TROOT.h"
00013 
00014 #include "BeamDataUtil/BeamMonSpill.h"
00015 //#include "MCReweight/BeamType.h"
00016 #include "Conventions/BeamType.h"//used here from 25/Feb/07
00017 #include "BeamDataUtil/BMSpillAna.h"
00018 #include "DataUtil/MCInfo.h"
00019 #include "MessageService/MsgService.h"
00020 #include "Conventions/Munits.h"
00021 #include "BeamDataNtuple/NtpBDLiteRecord.h"
00022 #include "StandardNtuple/NtpStRecord.h"
00023 #include "Conventions/ReleaseType.h"//used here from 25/Feb/07
00024 
00025 #include "NtupleUtils/NuEvent.h"
00026 #include "NtupleUtils/NuUtilities.h"
00027 
00028 #include "NuMuBar/NuConfig.h"
00029 #include "NuMuBar/NuBeam.h"
00030 
00031 CVSID("$Id: NuBeam.cxx,v 1.11 2009/12/20 17:49:40 ahimmel Exp $");
00032 
00033 //......................................................................
00034 
00035 NuBeam::NuBeam()
00036 {
00037   MSG("NuBeam",Msg::kDebug)
00038     <<"Running NuBeam Constructor..."<<endl;
00039 
00040 
00041   MSG("NuBeam",Msg::kDebug)
00042     <<"Finished NuBeam Constructor"<<endl;
00043 }
00044 
00045 //......................................................................
00046 
00047 NuBeam::~NuBeam()
00048 {
00049   MSG("NuBeam",Msg::kDebug)
00050     <<"Running NuBeam Destructor..."<<endl;
00051   
00052 
00053   MSG("NuBeam",Msg::kDebug)
00054     <<"Finished NuBeam Destructor"<<endl;
00055 }
00056 
00057 //......................................................................
00058 
00059 Bool_t NuBeam::IsGoodSpillAndFillPot(const NtpBDLiteRecord* recBD,
00060                                      const NuConfig& config,
00061                                      NuEvent& nu) const
00062 {
00067 
00068   static TH1F* hPottor101=0;
00069   static TH1F* hPottr101d=0;
00070   static TH1F* hPottortgt=0;
00071   static TH1F* hPottrtgtd=0;
00072 
00073   //good pot
00074   static TH1F* hPotGoodtor101=0;
00075   static TH1F* hPotGoodtr101d=0;
00076   static TH1F* hPotGoodtortgt=0;
00077   static TH1F* hPotGoodtrtgtd=0;
00078 
00079   //bad pot
00080   static TH1F* hPotBadtor101=0;
00081   static TH1F* hPotBadtr101d=0;
00082   static TH1F* hPotBadtortgt=0;
00083   static TH1F* hPotBadtrtgtd=0;
00084   
00085   //horncur, etc
00086   static TH1F* hHorncur=0;
00087   static TH1F* hHorncurBad=0;
00088   static TH1F* hBwidx=0;
00089   static TH1F* hBwidxBad=0;
00090   static TH1F* hBwidy=0;
00091   static TH1F* hBwidyBad=0;
00092   static TH1F* hBposx=0;
00093   static TH1F* hBposxBad=0;
00094   static TH1F* hBposy=0;
00095   static TH1F* hBposyBad=0;
00096   static TH1F* hBpmint=0;
00097   static TH1F* hBpmintBad=0;
00098 
00099   static TH1F* hAvBposx=0;
00100   static TH1F* hAvBposxBad=0;
00101   static TH1F* hAvBposy=0;
00102   static TH1F* hAvBposyBad=0;
00103 
00104   static TH1F* hDeltaBposx=0;
00105   static TH1F* hDeltaBposxBad=0;
00106   static TH1F* hDeltaBposy=0;
00107   static TH1F* hDeltaBposyBad=0;
00108 
00109   static TH1F* hBDTimeDiff=0;
00110   static TH1F* hBDTimeDiffBad=0;
00111   static TH1F* hBDTimeDiffBadBad=0;
00112 
00113   static TH1F* hBeamType=0;
00114   static TH1F* hBeamTypeBad=0;
00115   static TH1F* hTargetIn=0;
00116   static TH1F* hTargetInBad=0;
00117 
00118   static TH1F* hSpillsPerFile=0;
00119   
00120   if (!hPottor101){
00121     MAXMSG("NuBeam",Msg::kDebug,1)
00122       <<"Creating POT plots..."<<endl;
00123     
00124     //these 1st four POT histograms also take into account the coil
00125     //the PotGood histos below are for all spills where the beam
00126     //was good (regardless of the coil)
00127     hPottor101=new TH1F("hPottor101","hPottor101",100000,-10,90);
00128     hPottor101->GetXaxis()->SetTitle("POT (x10^{12})");
00129     hPottor101->GetXaxis()->CenterTitle();
00130     hPottor101->GetYaxis()->SetTitle("");
00131     hPottor101->GetYaxis()->CenterTitle();
00132     hPottor101->SetFillColor(0);
00133     hPottor101->SetLineColor(1);
00134     //hPottor101->SetBit(TH1::kCanRebin);
00135 
00136     hPottr101d=new TH1F("hPottr101d","hPottr101d",100000,-10,90);
00137     hPottr101d->GetXaxis()->SetTitle("POT (x10^{12})");
00138     hPottr101d->GetXaxis()->CenterTitle();
00139     hPottr101d->GetYaxis()->SetTitle("");
00140     hPottr101d->GetYaxis()->CenterTitle();
00141     hPottr101d->SetFillColor(0);
00142     hPottr101d->SetLineColor(1);
00143     //hPottr101d->SetBit(TH1::kCanRebin);
00144 
00145     hPottortgt=new TH1F("hPottortgt","hPottortgt",100000,-10,90);
00146     hPottortgt->GetXaxis()->SetTitle("POT (x10^{12})");
00147     hPottortgt->GetXaxis()->CenterTitle();
00148     hPottortgt->GetYaxis()->SetTitle("");
00149     hPottortgt->GetYaxis()->CenterTitle();
00150     hPottortgt->SetFillColor(0);
00151     hPottortgt->SetLineColor(1);
00152     //hPottortgt->SetBit(TH1::kCanRebin);
00153 
00154     hPottrtgtd=new TH1F("hPottrtgtd","hPottrtgtd",100000,-10,90);
00155     hPottrtgtd->GetXaxis()->SetTitle("POT (x10^{12})");
00156     hPottrtgtd->GetXaxis()->CenterTitle();
00157     hPottrtgtd->GetYaxis()->SetTitle("");
00158     hPottrtgtd->GetYaxis()->CenterTitle();
00159     hPottrtgtd->SetFillColor(0);
00160     hPottrtgtd->SetLineColor(1);
00161     //hPottrtgtd->SetBit(TH1::kCanRebin);
00162 
00163     
00164     //Good POT
00165     hPotGoodtor101=new TH1F("hPotGoodtor101","hPotGoodtor101",
00166                            100000,-10,90);
00167     hPotGoodtor101->GetXaxis()->SetTitle("POT (x10^{12})");
00168     hPotGoodtor101->GetXaxis()->CenterTitle();
00169     hPotGoodtor101->GetYaxis()->SetTitle("");
00170     hPotGoodtor101->GetYaxis()->CenterTitle();
00171     hPotGoodtor101->SetFillColor(0);
00172     hPotGoodtor101->SetLineColor(1);
00173     //hPotGoodtor101->SetBit(TH1::kCanRebin);
00174 
00175     hPotGoodtr101d=new TH1F("hPotGoodtr101d","hPotGoodtr101d",
00176                            100000,-10,90);
00177     hPotGoodtr101d->GetXaxis()->SetTitle("POT (x10^{12})");
00178     hPotGoodtr101d->GetXaxis()->CenterTitle();
00179     hPotGoodtr101d->GetYaxis()->SetTitle("");
00180     hPotGoodtr101d->GetYaxis()->CenterTitle();
00181     hPotGoodtr101d->SetFillColor(0);
00182     hPotGoodtr101d->SetLineColor(1);
00183     //hPotGoodtr101d->SetBit(TH1::kCanRebin);
00184 
00185     hPotGoodtortgt=new TH1F("hPotGoodtortgt","hPotGoodtortgt",
00186                            100000,-10,90);
00187     hPotGoodtortgt->GetXaxis()->SetTitle("POT (x10^{12})");
00188     hPotGoodtortgt->GetXaxis()->CenterTitle();
00189     hPotGoodtortgt->GetYaxis()->SetTitle("");
00190     hPotGoodtortgt->GetYaxis()->CenterTitle();
00191     hPotGoodtortgt->SetFillColor(0);
00192     hPotGoodtortgt->SetLineColor(1);
00193     //hPotGoodtortgt->SetBit(TH1::kCanRebin);
00194 
00195     hPotGoodtrtgtd=new TH1F("hPotGoodtrtgtd","hPotGoodtrtgtd",
00196                            100000,-10,90);
00197     hPotGoodtrtgtd->GetXaxis()->SetTitle("POT (x10^{12})");
00198     hPotGoodtrtgtd->GetXaxis()->CenterTitle();
00199     hPotGoodtrtgtd->GetYaxis()->SetTitle("");
00200     hPotGoodtrtgtd->GetYaxis()->CenterTitle();
00201     hPotGoodtrtgtd->SetFillColor(0);
00202     hPotGoodtrtgtd->SetLineColor(1);
00203     //hPotGoodtrtgtd->SetBit(TH1::kCanRebin);
00204 
00205 
00206     //bad POT
00207     hPotBadtor101=new TH1F("hPotBadtor101","hPotBadtor101",
00208                            100000,-10,90);
00209     hPotBadtor101->GetXaxis()->SetTitle("POT (x10^{12})");
00210     hPotBadtor101->GetXaxis()->CenterTitle();
00211     hPotBadtor101->GetYaxis()->SetTitle("");
00212     hPotBadtor101->GetYaxis()->CenterTitle();
00213     hPotBadtor101->SetFillColor(0);
00214     hPotBadtor101->SetLineColor(1);
00215     //hPotBadtor101->SetBit(TH1::kCanRebin);
00216 
00217     hPotBadtr101d=new TH1F("hPotBadtr101d","hPotBadtr101d",
00218                            100000,-10,90);
00219     hPotBadtr101d->GetXaxis()->SetTitle("POT (x10^{12})");
00220     hPotBadtr101d->GetXaxis()->CenterTitle();
00221     hPotBadtr101d->GetYaxis()->SetTitle("");
00222     hPotBadtr101d->GetYaxis()->CenterTitle();
00223     hPotBadtr101d->SetFillColor(0);
00224     hPotBadtr101d->SetLineColor(1);
00225     //hPotBadtr101d->SetBit(TH1::kCanRebin);
00226 
00227     hPotBadtortgt=new TH1F("hPotBadtortgt","hPotBadtortgt",
00228                            100000,-10,90);
00229     hPotBadtortgt->GetXaxis()->SetTitle("POT (x10^{12})");
00230     hPotBadtortgt->GetXaxis()->CenterTitle();
00231     hPotBadtortgt->GetYaxis()->SetTitle("");
00232     hPotBadtortgt->GetYaxis()->CenterTitle();
00233     hPotBadtortgt->SetFillColor(0);
00234     hPotBadtortgt->SetLineColor(1);
00235     //hPotBadtortgt->SetBit(TH1::kCanRebin);
00236 
00237     hPotBadtrtgtd=new TH1F("hPotBadtrtgtd","hPotBadtrtgtd",
00238                            100000,-10,90);
00239     hPotBadtrtgtd->GetXaxis()->SetTitle("POT (x10^{12})");
00240     hPotBadtrtgtd->GetXaxis()->CenterTitle();
00241     hPotBadtrtgtd->GetYaxis()->SetTitle("");
00242     hPotBadtrtgtd->GetYaxis()->CenterTitle();
00243     hPotBadtrtgtd->SetFillColor(0);
00244     hPotBadtrtgtd->SetLineColor(1);
00245     //hPotBadtrtgtd->SetBit(TH1::kCanRebin);
00246 
00247 
00248     //horn current
00249     hHorncur=new TH1F("hHorncur","hHorncur",
00250                       1000,-250,250);
00251     hHorncur->GetXaxis()->SetTitle("Horn Current");
00252     hHorncur->GetXaxis()->CenterTitle();
00253     hHorncur->GetYaxis()->SetTitle("");
00254     hHorncur->GetYaxis()->CenterTitle();
00255     hHorncur->SetFillColor(0);
00256     hHorncur->SetLineColor(1);
00257     //hHorncur->SetBit(TH1::kCanRebin);
00258 
00259     hHorncurBad=new TH1F("hHorncurBad","hHorncurBad",
00260                          1000,-250,250);
00261     hHorncurBad->GetXaxis()->SetTitle("Horn Current");
00262     hHorncurBad->GetXaxis()->CenterTitle();
00263     hHorncurBad->GetYaxis()->SetTitle("");
00264     hHorncurBad->GetYaxis()->CenterTitle();
00265     hHorncurBad->SetFillColor(0);
00266     hHorncurBad->SetLineColor(1);
00267     //hHorncurBad->SetBit(TH1::kCanRebin);
00268 
00269 
00270 
00271     hBwidx=new TH1F("hBwidx","hBwidx",
00272                     1000,-0.002,0.01);
00273     hBwidx->GetXaxis()->SetTitle("Beam Width x");
00274     hBwidx->GetXaxis()->CenterTitle();
00275     hBwidx->GetYaxis()->SetTitle("");
00276     hBwidx->GetYaxis()->CenterTitle();
00277     hBwidx->SetFillColor(0);
00278     hBwidx->SetLineColor(1);
00279     //hBwidx->SetBit(TH1::kCanRebin);
00280 
00281     hBwidxBad=new TH1F("hBwidxBad","hBwidxBad",
00282                        1000,-0.002,0.01);
00283     hBwidxBad->GetXaxis()->SetTitle("Beam Width x");
00284     hBwidxBad->GetXaxis()->CenterTitle();
00285     hBwidxBad->GetYaxis()->SetTitle("");
00286     hBwidxBad->GetYaxis()->CenterTitle();
00287     hBwidxBad->SetFillColor(0);
00288     hBwidxBad->SetLineColor(1);
00289     //hBwidxBad->SetBit(TH1::kCanRebin);
00290 
00291     hBwidy=new TH1F("hBwidy","hBwidy",
00292                     1000,-0.002,0.01);
00293     hBwidy->GetXaxis()->SetTitle("Beam Width y");
00294     hBwidy->GetXaxis()->CenterTitle();
00295     hBwidy->GetYaxis()->SetTitle("");
00296     hBwidy->GetYaxis()->CenterTitle();
00297     hBwidy->SetFillColor(0);
00298     hBwidy->SetLineColor(1);
00299     //hBwidy->SetBit(TH1::kCanRebin);
00300 
00301     hBwidyBad=new TH1F("hBwidyBad","hBwidyBad",
00302                        1000,-0.002,0.01);
00303     hBwidyBad->GetXaxis()->SetTitle("Beam Width y");
00304     hBwidyBad->GetXaxis()->CenterTitle();
00305     hBwidyBad->GetYaxis()->SetTitle("");
00306     hBwidyBad->GetYaxis()->CenterTitle();
00307     hBwidyBad->SetFillColor(0);
00308     hBwidyBad->SetLineColor(1);
00309     //hBwidyBad->SetBit(TH1::kCanRebin);
00310 
00311 
00312     //beam intensity per batch
00313     hBpmint=new TH1F("hBpmint","hBpmint",
00314                      2010,-10,2000);
00315     hBpmint->GetXaxis()->SetTitle("Beam intensity");
00316     hBpmint->GetXaxis()->CenterTitle();
00317     hBpmint->GetYaxis()->SetTitle("");
00318     hBpmint->GetYaxis()->CenterTitle();
00319     hBpmint->SetFillColor(0);
00320     hBpmint->SetLineColor(1);
00321     //hBpmint->SetBit(TH1::kCanRebin);
00322     
00323     hBpmintBad=new TH1F("hBpmintBad","hBpmintBad",
00324                         2010,-10,2000);
00325     hBpmintBad->GetXaxis()->SetTitle("Beam intensity");
00326     hBpmintBad->GetXaxis()->CenterTitle();
00327     hBpmintBad->GetYaxis()->SetTitle("");
00328     hBpmintBad->GetYaxis()->CenterTitle();
00329     hBpmintBad->SetFillColor(0);
00330     hBpmintBad->SetLineColor(1);
00331     //hBpmintBad->SetBit(TH1::kCanRebin);
00332 
00333 
00334     hBposx=new TH1F("hBposx","hBposx",
00335                     5000,-0.01,0.01);
00336     hBposx->GetXaxis()->SetTitle("Beam Position x");
00337     hBposx->GetXaxis()->CenterTitle();
00338     hBposx->GetYaxis()->SetTitle("");
00339     hBposx->GetYaxis()->CenterTitle();
00340     hBposx->SetFillColor(0);
00341     hBposx->SetLineColor(1);
00342     //hBposx->SetBit(TH1::kCanRebin);
00343 
00344     hBposxBad=new TH1F("hBposxBad","hBposxBad",
00345                        5000,-0.01,0.01);
00346     hBposxBad->GetXaxis()->SetTitle("Beam Position x");
00347     hBposxBad->GetXaxis()->CenterTitle();
00348     hBposxBad->GetYaxis()->SetTitle("");
00349     hBposxBad->GetYaxis()->CenterTitle();
00350     hBposxBad->SetFillColor(0);
00351     hBposxBad->SetLineColor(1);
00352     //hBposxBad->SetBit(TH1::kCanRebin);
00353 
00354     hBposy=new TH1F("hBposy","hBposy",
00355                     5000,-0.01,0.01);//
00356     hBposy->GetXaxis()->SetTitle("Beam Position y");
00357     hBposy->GetXaxis()->CenterTitle();
00358     hBposy->GetYaxis()->SetTitle("");
00359     hBposy->GetYaxis()->CenterTitle();
00360     hBposy->SetFillColor(0);
00361     hBposy->SetLineColor(1);
00362     //hBposy->SetBit(TH1::kCanRebin);
00363 
00364     hBposyBad=new TH1F("hBposyBad","hBposyBad",
00365                        5000,-0.01,0.01);
00366     hBposyBad->GetXaxis()->SetTitle("Beam Position y");
00367     hBposyBad->GetXaxis()->CenterTitle();
00368     hBposyBad->GetYaxis()->SetTitle("");
00369     hBposyBad->GetYaxis()->CenterTitle();
00370     hBposyBad->SetFillColor(0);
00371     hBposyBad->SetLineColor(1);
00372     //hBposyBad->SetBit(TH1::kCanRebin);
00373 
00374 
00375 
00376     hAvBposx=new TH1F("hAvBposx","hAvBposx",
00377                     5000,-0.01,0.01);
00378     hAvBposx->GetXaxis()->SetTitle("Av. Beam Position x");
00379     hAvBposx->GetXaxis()->CenterTitle();
00380     hAvBposx->GetYaxis()->SetTitle("");
00381     hAvBposx->GetYaxis()->CenterTitle();
00382     hAvBposx->SetFillColor(0);
00383     hAvBposx->SetLineColor(1);
00384     //hAvBposx->SetBit(TH1::kCanRebin);
00385 
00386     hAvBposxBad=new TH1F("hAvBposxBad","hAvBposxBad",
00387                        5000,-0.01,0.01);
00388     hAvBposxBad->GetXaxis()->SetTitle("Av. Beam Position x");
00389     hAvBposxBad->GetXaxis()->CenterTitle();
00390     hAvBposxBad->GetYaxis()->SetTitle("");
00391     hAvBposxBad->GetYaxis()->CenterTitle();
00392     hAvBposxBad->SetFillColor(0);
00393     hAvBposxBad->SetLineColor(1);
00394     //hAvBposxBad->SetBit(TH1::kCanRebin);
00395 
00396     hAvBposy=new TH1F("hAvBposy","hAvBposy",
00397                     5000,-0.01,0.01);//
00398     hAvBposy->GetXaxis()->SetTitle("Av. Beam Position y");
00399     hAvBposy->GetXaxis()->CenterTitle();
00400     hAvBposy->GetYaxis()->SetTitle("");
00401     hAvBposy->GetYaxis()->CenterTitle();
00402     hAvBposy->SetFillColor(0);
00403     hAvBposy->SetLineColor(1);
00404     //hAvBposy->SetBit(TH1::kCanRebin);
00405 
00406     hAvBposyBad=new TH1F("hAvBposyBad","hAvBposyBad",
00407                        5000,-0.01,0.01);
00408     hAvBposyBad->GetXaxis()->SetTitle("Av. Beam Position y");
00409     hAvBposyBad->GetXaxis()->CenterTitle();
00410     hAvBposyBad->GetYaxis()->SetTitle("");
00411     hAvBposyBad->GetYaxis()->CenterTitle();
00412     hAvBposyBad->SetFillColor(0);
00413     hAvBposyBad->SetLineColor(1);
00414     //hAvBposyBad->SetBit(TH1::kCanRebin);
00415 
00416 
00417 
00418     hDeltaBposx=new TH1F("hDeltaBposx","hDeltaBposx",
00419                          5000,-0.01,0.01);//
00420     hDeltaBposx->GetXaxis()->SetTitle("Delta Beam Position x");
00421     hDeltaBposx->GetXaxis()->CenterTitle();
00422     hDeltaBposx->GetYaxis()->SetTitle("");
00423     hDeltaBposx->GetYaxis()->CenterTitle();
00424     hDeltaBposx->SetFillColor(0);
00425     hDeltaBposx->SetLineColor(1);
00426     //hDeltaBposx->SetBit(TH1::kCanRebin);
00427 
00428     hDeltaBposxBad=new TH1F("hDeltaBposxBad","hDeltaBposxBad",
00429                             5000,-0.01,0.01);
00430     hDeltaBposxBad->GetXaxis()->SetTitle("Delta Beam Position x");
00431     hDeltaBposxBad->GetXaxis()->CenterTitle();
00432     hDeltaBposxBad->GetYaxis()->SetTitle("");
00433     hDeltaBposxBad->GetYaxis()->CenterTitle();
00434     hDeltaBposxBad->SetFillColor(0);
00435     hDeltaBposxBad->SetLineColor(1);
00436     //hDeltaBposxBad->SetBit(TH1::kCanRebin);
00437 
00438     hDeltaBposy=new TH1F("hDeltaBposy","hDeltaBposy",
00439                          5000,-0.01,0.01);//
00440     hDeltaBposy->GetXaxis()->SetTitle("Delta Beam Position y");
00441     hDeltaBposy->GetXaxis()->CenterTitle();
00442     hDeltaBposy->GetYaxis()->SetTitle("");
00443     hDeltaBposy->GetYaxis()->CenterTitle();
00444     hDeltaBposy->SetFillColor(0);
00445     hDeltaBposy->SetLineColor(1);
00446     //hDeltaBposy->SetBit(TH1::kCanRebin);
00447 
00448     hDeltaBposyBad=new TH1F("hDeltaBposyBad","hDeltaBposyBad",
00449                             5000,-0.01,0.01);
00450     hDeltaBposyBad->GetXaxis()->SetTitle("Delta Beam Position y");
00451     hDeltaBposyBad->GetXaxis()->CenterTitle();
00452     hDeltaBposyBad->GetYaxis()->SetTitle("");
00453     hDeltaBposyBad->GetYaxis()->CenterTitle();
00454     hDeltaBposyBad->SetFillColor(0);
00455     hDeltaBposyBad->SetLineColor(1);
00456     //hDeltaBposyBad->SetBit(TH1::kCanRebin);
00457 
00458 
00459 
00460     //timediff
00461     hBDTimeDiff=new TH1F("hBDTimeDiff","hBDTimeDiff",
00462                          5000,-2.5,2.5);
00463     hBDTimeDiff->GetXaxis()->SetTitle("Beam/ND Time Difference (s)");
00464     hBDTimeDiff->GetXaxis()->CenterTitle();
00465     hBDTimeDiff->GetYaxis()->SetTitle("");
00466     hBDTimeDiff->GetYaxis()->CenterTitle();
00467     hBDTimeDiff->SetFillColor(0);
00468     hBDTimeDiff->SetLineColor(1);
00469     //hBDTimeDiff->SetBit(TH1::kCanRebin);
00470 
00471     hBDTimeDiffBad=new TH1F("hBDTimeDiffBad","hBDTimeDiffBad",
00472                             5000,-2.5,2.5);
00473     hBDTimeDiffBad->GetXaxis()->SetTitle("Beam/ND Time Difference (s)");
00474     hBDTimeDiffBad->GetXaxis()->CenterTitle();
00475     hBDTimeDiffBad->GetYaxis()->SetTitle("");
00476     hBDTimeDiffBad->GetYaxis()->CenterTitle();
00477     hBDTimeDiffBad->SetFillColor(0);
00478     hBDTimeDiffBad->SetLineColor(1);
00479     //hBDTimeDiffBad->SetBit(TH1::kCanRebin);
00480 
00481     hBDTimeDiffBadBad=new TH1F("hBDTimeDiffBadBad","hBDTimeDiffBadBad",
00482                                5000,-2.5,1000);
00483     hBDTimeDiffBadBad->GetXaxis()->
00484       SetTitle("Beam/ND Time Difference (s)");
00485     hBDTimeDiffBadBad->GetXaxis()->CenterTitle();
00486     hBDTimeDiffBadBad->GetYaxis()->SetTitle("");
00487     hBDTimeDiffBadBad->GetYaxis()->CenterTitle();
00488     hBDTimeDiffBadBad->SetFillColor(0);
00489     hBDTimeDiffBadBad->SetLineColor(1);
00490     //hBDTimeDiffBadBad->SetBit(TH1::kCanRebin);
00491 
00492 
00493     
00494     hBeamType=new TH1F("hBeamType","hBeamType",
00495                        100,-10,90);//
00496     hBeamType->GetXaxis()->SetTitle("Beam Type");
00497     hBeamType->GetXaxis()->CenterTitle();
00498     hBeamType->GetYaxis()->SetTitle("");
00499     hBeamType->GetYaxis()->CenterTitle();
00500     hBeamType->SetFillColor(0);
00501     hBeamType->SetLineColor(1);
00502     //hBeamType->SetBit(TH1::kCanRebin);
00503 
00504     hBeamTypeBad=new TH1F("hBeamTypeBad","hBeamTypeBad",
00505                           100,-10,90);
00506     hBeamTypeBad->GetXaxis()->SetTitle("Beam Type");
00507     hBeamTypeBad->GetXaxis()->CenterTitle();
00508     hBeamTypeBad->GetYaxis()->SetTitle("");
00509     hBeamTypeBad->GetYaxis()->CenterTitle();
00510     hBeamTypeBad->SetFillColor(0);
00511     hBeamTypeBad->SetLineColor(1);
00512     //hBeamTypeBad->SetBit(TH1::kCanRebin);
00513 
00514 
00515     hTargetIn=new TH1F("hTargetIn","hTargetIn",
00516                        100,-10,90);//
00517     hTargetIn->GetXaxis()->SetTitle("Target In");
00518     hTargetIn->GetXaxis()->CenterTitle();
00519     hTargetIn->GetYaxis()->SetTitle("");
00520     hTargetIn->GetYaxis()->CenterTitle();
00521     hTargetIn->SetFillColor(0);
00522     hTargetIn->SetLineColor(1);
00523     //hTargetIn->SetBit(TH1::kCanRebin);
00524 
00525     hTargetInBad=new TH1F("hTargetInBad","hTargetInBad",
00526                           100,-10,90);
00527     hTargetInBad->GetXaxis()->SetTitle("Target In");
00528     hTargetInBad->GetXaxis()->CenterTitle();
00529     hTargetInBad->GetYaxis()->SetTitle("");
00530     hTargetInBad->GetYaxis()->CenterTitle();
00531     hTargetInBad->SetFillColor(0);
00532     hTargetInBad->SetLineColor(1);
00533     //hTargetInBad->SetBit(TH1::kCanRebin);
00534 
00535     hSpillsPerFile=new TH1F("hSpillsPerFile","hSpillsPerFile",
00536                             10000,0,10000);
00537     hSpillsPerFile->GetXaxis()->SetTitle("hSpillsPerFile");
00538     hSpillsPerFile->GetXaxis()->CenterTitle();
00539     hSpillsPerFile->GetYaxis()->SetTitle("");
00540     hSpillsPerFile->GetYaxis()->CenterTitle();
00541     hSpillsPerFile->SetFillColor(0);
00542     hSpillsPerFile->SetLineColor(1);
00543     //hSpillsPerFile->SetBit(TH1::kCanRebin);
00544   }
00545   
00546   if (!recBD) { //this is the MC case
00547     //this is the pot for carrot MC at LE-10
00548     Float_t potPerSpill=24.2;//x1e12 but ntp assumes e12
00549     Float_t spillsPerFile=400;
00550     
00551     //This function, for e.g., returns the value 24.2, since PoT/snarl 
00552     //for Daikon L010185 Near MC is 2.42*1013.
00553 
00554     //Convert to canonical BeamType in the case of intensity
00555     //weighted MC
00556     NuUtilities utils;
00557     BeamType::BeamType_t convBeamType =
00558       utils.ReturnConventionsBeamType(config);
00559     MAXMSG("NuBeam",Msg::kInfo,1)
00560       << "Beam type to be used for PoT counting: "
00561       << BeamType::AsString(convBeamType)
00562       << endl;
00563 
00564     potPerSpill=MCInfo::GetMCPoT
00565       (static_cast<Detector::Detector_t>(config.detector),
00566        static_cast<BeamType::BeamType_t>(convBeamType),
00567        static_cast<ReleaseType::Release_t>(config.mcVersion));
00568     spillsPerFile=MCInfo::GetNoSnarlPerFile
00569       (static_cast<Detector::Detector_t>(config.detector),
00570        static_cast<BeamType::BeamType_t>(convBeamType),
00571        static_cast<ReleaseType::Release_t>(config.mcVersion),
00572        config.runPeriod);
00573     
00574     //remove the FD factor of 1e8 so that it can be stored in the histo
00575     //also just put a 1 for spills per file, since it's pot/file really
00576     if (config.detector==Detector::kFar) {
00577       if (potPerSpill < 0) {
00578         MAXMSG("NuBeam",Msg::kWarning,50) << "*** Warning: potPerSpill undefined for beam type: " 
00579         << BeamType::AsString(convBeamType) << " (" << convBeamType << ") in the far detector." << endl
00580         << "ASSUMING 6.5e20 per file, which is standard for Daikon.  " << endl
00581         << "This is dangerous - make sure the above beam type is added to DataUtil/MCInfo::GetMCPoT() ASAP." << endl;
00582         potPerSpill = 6.5e8;
00583       }
00584       potPerSpill/=1e8;
00585       if (spillsPerFile<=0) {
00586         MAXMSG("NuBeam",Msg::kInfo,1)
00587           <<"Overriding spills per file (="<<spillsPerFile
00588           <<"), now set to 1"<<endl;
00589         spillsPerFile=1;
00590       }
00591     }
00592 
00593     //sanity check: histogram overflow
00594     Int_t n=hSpillsPerFile->GetNbinsX();
00595     if (hSpillsPerFile->GetBinContent(n+1)>0){
00596       MSG("NuBeam",Msg::kError)
00597         <<"n bins="<<n
00598         <<", overflow="<<hSpillsPerFile->GetBinContent(n+1)<<endl; 
00599     }
00600 
00601     MAXMSG("NuBeam",Msg::kInfo,1)
00602       <<"No NtpBDLiteRecord, so using potPerSpill="<<potPerSpill
00603       <<"x10^12 ND (10^20 FD)"<<" and spillsPerFile="<<spillsPerFile
00604       <<" for beamType="<<convBeamType
00605       <<", hornCurrent="<<config.hornCurrent
00606       <<", mcVersion="<<config.mcVersion<<endl;
00607     
00608     //fill histos
00609     hPottor101->Fill(potPerSpill);
00610     hPottr101d->Fill(potPerSpill);
00611     hPottortgt->Fill(potPerSpill);
00612     hPottrtgtd->Fill(potPerSpill);
00613     hSpillsPerFile->Fill(spillsPerFile);
00614 
00615     //write out the potPerSpill
00616     nu.pot=potPerSpill;
00617     nu.goodBeamSntp=true;//always true for MC
00618     
00619     //don't fill the bad histos ever in this case
00620   }
00622   //this is the data case
00624   else {
00625     MAXMSG("NuBeam",Msg::kInfo,1)
00626       <<"Using NtpBDLiteRecord for POT"<<endl;
00627     
00628     const NtpBDLiteRecord& ntpBD=*recBD;
00629     
00630     static BeamMonSpill spill;
00631     static BMSpillAna bmsa;
00632         
00633     //tell it to use the database defined values of the cuts
00634     bmsa.UseDatabaseCuts();
00635     
00636     Bool_t changeBeamCuts=false;
00637     //change default values (only once though)
00638     static Registry* pr=0;
00639     if (!pr && changeBeamCuts){
00640       MAXMSG("NuBeam",Msg::kInfo,1)
00641         <<"Changing default values of BMSpillAna"<<endl;
00642       pr=new Registry();
00643       Registry& r=*pr;//get a reference
00644       r.UnLockValues();
00645       //r.Set("TargetIn",1);
00646       //r.Set("BeamType",-1);//don't care if ME,LE,etc
00647       //r.Set("PosTgtXMin",-2.0*Munits::mm);
00648       //have to make sure that zero is not included
00649       r.Set("PosTgtXMax",-0.0000001*Munits::mm);
00650       r.Set("PosTgtYMin", 0.0000001*Munits::mm);
00651       //r.Set("PosTgtYMax", 2.0*Munits::mm);
00652       r.LockValues();
00653       bmsa.Config(r);
00654       
00655       //print the registry
00656       r.PrettyPrint(cout);
00657     }
00658     
00659     //give it the spill
00660     bmsa.SetSpill(ntpBD,spill);
00661 
00662     Double_t xmean=0;
00663     Double_t ymean=0;
00664     Double_t xrms=0;
00665     Double_t yrms=0;
00666     bmsa.GetSpill().BpmAtTarget(xmean,ymean,xrms,yrms);
00667     
00668     Int_t beamType=spill.GetStatusBits().beam_type;
00669     Int_t targetIn=spill.GetStatusBits().target_in;
00670     
00671     //write to the object whether the beam is good
00672     nu.goodBeamSntp=bmsa.SelectSpill();
00673     //set the pot value
00674     nu.pot=ntpBD.tortgt;
00675     if (nu.pot<0) nu.pot=0;//bad spills can be <0 but don't want to sum
00676     Float_t tor101=ntpBD.tor101;
00677     if (tor101<0) tor101=0;
00678     Float_t tr101d=ntpBD.tr101d;
00679     if (tr101d<0) tr101d=0;
00680     Float_t tortgt=ntpBD.tortgt;
00681     if (tortgt<0) tortgt=0;
00682     Float_t trtgtd=ntpBD.trtgtd;
00683     if (trtgtd<0) trtgtd=0;    
00684 
00685     //check if a good spill
00686     if (!nu.goodBeamSntp) {
00687 
00688       MAXMSG("NuBeam",Msg::kWarning, 20)
00689       <<"NuBeam: BAD SPILL:" << endl
00690       << "  > Maximum time diffrence (s):     " << ntpBD.GetHeader().GetTimeDiffStreamSpill() << endl
00691       << "  > Spill intensity (1e12 pot):     " << trtgtd << endl
00692       << "  > Horn Current (kA):              " << ntpBD.horncur << endl
00693       << "  > fBeamType:                      " << BeamType::AsString(static_cast<BeamType::BeamType_t>(spill.BeamType())) <<" ("<<spill.BeamType()<<")"<<endl
00694       << "  > Horizontal beam position (mm):  " << xmean << endl
00695       << "  > Vertical beam position (mm):    " << ymean << endl
00696       << "  > Horizontal beam width a (mm):   " << ntpBD.bwidx << endl
00697       << "  > Vertical beam width a (mm):     " << ntpBD.bwidy << endl;
00698 
00699       //for the Bad plots to have any meaning you have to make sure
00700       //that the spill was in time, otherwise you just plot
00701       //data from spills ages ago.
00702       if (ntpBD.GetHeader().GetTimeDiffStreamSpill()>1){
00703         hBDTimeDiffBadBad->Fill(ntpBD.GetHeader().
00704                                 GetTimeDiffStreamSpill());
00705       }
00706       else {
00707         hPotBadtor101->Fill(tor101);
00708         hPotBadtr101d->Fill(tr101d);
00709         hPotBadtortgt->Fill(tortgt);
00710         hPotBadtrtgtd->Fill(trtgtd);
00711         hHorncurBad->Fill(ntpBD.horncur);
00712         hBwidxBad->Fill(ntpBD.bwidx);
00713         hBwidyBad->Fill(ntpBD.bwidy);
00714         hBDTimeDiffBad->Fill(ntpBD.GetHeader().
00715                              GetTimeDiffStreamSpill());
00716         
00717         hAvBposxBad->Fill(xmean);
00718         hAvBposyBad->Fill(ymean);
00719         
00720         for (Int_t batch=0;batch<6;batch++){
00721           hBpmintBad->Fill(ntpBD.bpmint[batch]);
00722           //check if batch existed in spill
00723           if (ntpBD.bpmint[batch]>1) {
00724             hBposxBad->Fill(ntpBD.bposx[batch]);
00725             hBposyBad->Fill(ntpBD.bposy[batch]);
00726             
00727             //calc delta from mean position
00728             //should probably use some sort of intensity weighting
00729             hDeltaBposxBad->Fill(ntpBD.bposx[batch]-xmean);
00730             hDeltaBposyBad->Fill(ntpBD.bposy[batch]-ymean);
00731           }
00732         }
00733 
00734         hTargetInBad->Fill(targetIn);
00735         hBeamTypeBad->Fill(beamType);
00736       }
00737       return nu.goodBeamSntp;
00738     }
00739     
00740     MAXMSG("NuBeam",Msg::kInfo,1)
00741       <<"Found first good spill"<<", trtgtd="<<trtgtd
00742       <<", timeDiff="<<ntpBD.GetHeader().GetTimeDiffStreamSpill() 
00743       <<", beamType="<<beamType<<", tgtIn="<<targetIn<<endl;
00744   
00745     static Bool_t firstTime=true;
00746     if (firstTime) {
00747       bmsa.Print(); 
00748       cout<<endl;
00749       firstTime=false;
00750     }
00751 
00752     //perform sanity check that hornCurrent is 
00753     //within +/- 4% of the config one
00754     if (!(fabs(ntpBD.horncur)>fabs(config.hornCurrent)*0.96 &&
00755           fabs(ntpBD.horncur)<fabs(config.hornCurrent)*1.04)){
00756       MSG("NuBeam",Msg::kError)
00757         <<"HornCurrent problem for BeamType (sntp)="<<beamType
00758         <<", event hornCur="<<ntpBD.horncur
00759         <<" config.hornCurrent="<<config.hornCurrent
00760         <<endl;
00761       TH1F* hError=dynamic_cast<TH1F*>(gROOT->FindObject("hError"));
00762       if (hError) hError->Fill(1);//set an error flag
00763     }
00764 
00765     //these histograms are the ones used for POT counting
00766     //have to check the coil and data quality as well
00767     if (nu.goodBeamSntp && nu.coilIsOk && 
00768         (nu.isGoodDataQuality || !nu.useDBForDataQuality)) {
00769       hPottor101->Fill(tor101);
00770       hPottr101d->Fill(tr101d);
00771       hPottortgt->Fill(tortgt);
00772       hPottrtgtd->Fill(trtgtd);
00773     }
00774     
00775     //these are just for when the beam is good (regardless of coil)
00776     hPotGoodtor101->Fill(tor101);
00777     hPotGoodtr101d->Fill(tr101d);
00778     hPotGoodtortgt->Fill(tortgt);
00779     hPotGoodtrtgtd->Fill(trtgtd);
00780     hHorncur->Fill(ntpBD.horncur);
00781     hBwidx->Fill(ntpBD.bwidx);
00782     hBwidy->Fill(ntpBD.bwidy);
00783     hBDTimeDiff->Fill(ntpBD.GetHeader().GetTimeDiffStreamSpill());
00784 
00785     hAvBposx->Fill(xmean);
00786     hAvBposy->Fill(ymean);
00787     
00788     for (Int_t batch=0;batch<6;batch++){
00789       hBpmint->Fill(ntpBD.bpmint[batch]);
00790       //check if batch existed in spill
00791       if (ntpBD.bpmint[batch]>1) {
00792         hBposx->Fill(ntpBD.bposx[batch]);
00793         hBposy->Fill(ntpBD.bposy[batch]);
00794 
00795             
00796         //calc delta from mean position
00797         //should probably use some sort of intensity weighting
00798         hDeltaBposx->Fill(ntpBD.bposx[batch]-xmean);
00799         hDeltaBposy->Fill(ntpBD.bposy[batch]-ymean);
00800       }
00801     }
00802 
00803     hTargetIn->Fill(targetIn);
00804     hBeamType->Fill(beamType);
00805   }
00806 
00807   return nu.goodBeamSntp;
00808 }
00809 
00810 //......................................................................
00811 
00812 
00813 /*
00814 void NuBeam::GetMCPoTAndSnarlsPerFile(const NuConfig& config) const
00815 {
00817   assert(true);
00818   
00819   
00820   //bug fix to pot for MC: then there was another bug fix...
00821   //this is wrong now I believe...
00822   //L010170  2.493e13      400 snarls per file
00823   //L010185  2.486e13      400
00824   //L010200  2.497e13      400
00825   //L100200  1.034e13      500
00826   //L250200  1.141e13      130
00827 
00828   //FD MC is 440 files, 1.28e+023 pot total
00829   //1.28/440=2.909090909e20
00830   //Robert H says: 2.90619e20 pot per file
00831 
00832   //20070404 JJH:
00833   //from http://www-numi.fnal.gov/workgrps/mc/checkMCenstore/summary.html
00840 
00848 
00849 
00850   if (config.detector==Detector::kNear) {
00851     if (config.beamType==BeamType::kLE || 
00852         config.beamType==BeamType::k010){
00853       if (config.hornCurrent==-170){
00854         if (config.mcVersion==ReleaseType::kCarrot) {
00855           potPerSpill=25.05;//was 24.93;
00856           MAXMSG("NuBeam",Msg::kInfo,1)
00857             <<"Using potPerSpill="<<potPerSpill<<" for Carrot"<<endl;
00858           spillsPerFile=400;
00859         }
00860         else if (config.mcVersion==ReleaseType::kDaikon) {
00861           potPerSpill=24.2;
00862           MAXMSG("NuBeam",Msg::kInfo,1)
00863             <<"Using potPerSpill="<<potPerSpill<<" for Daikon"<<endl;
00864           spillsPerFile=800;
00865         }
00866         else {
00867           //default to Carrot value
00868           potPerSpill=25.05;
00869           MAXMSG("NuBeam",Msg::kWarning,1)
00870             <<"Defaulting to potPerSpill="<<potPerSpill
00871             <<" for unknown MC"<<endl;
00872           spillsPerFile=400;
00873         }
00874       }
00875       else if (config.hornCurrent==-185){
00876         if (config.mcVersion==ReleaseType::kCarrot) {
00877           potPerSpill=25.06;//was 24.86;
00878           MAXMSG("NuBeam",Msg::kInfo,1)
00879             <<"Using potPerSpill="<<potPerSpill<<" for Carrot"<<endl;
00880           spillsPerFile=400;
00881         }
00882         else if (config.mcVersion==ReleaseType::kDaikon) {
00883           potPerSpill=24.2;
00884           MAXMSG("NuBeam",Msg::kInfo,1)
00885             <<"Using potPerSpill="<<potPerSpill<<" for Daikon"<<endl;
00886           spillsPerFile=800;
00887         }
00888         else {
00889           //default to Carrot value
00890           potPerSpill=25.06;
00891           MAXMSG("NuBeam",Msg::kWarning,1)
00892             <<"Defaulting to potPerSpill="<<potPerSpill
00893             <<" for unknown MC"<<endl;
00894           spillsPerFile=400;
00895         }
00896       }
00897       else if (config.hornCurrent==-200){
00898         if (config.mcVersion==ReleaseType::kCarrot) {
00899           potPerSpill=25.23;//was 24.86;
00900           MAXMSG("NuBeam",Msg::kInfo,1)
00901             <<"Using potPerSpill="<<potPerSpill<<" for Carrot"<<endl;
00902           spillsPerFile=400;
00903         }
00904         else if (config.mcVersion==ReleaseType::kDaikon) {
00905           potPerSpill=24.2;
00906           MAXMSG("NuBeam",Msg::kInfo,1)
00907             <<"Using potPerSpill="<<potPerSpill<<" for Daikon"<<endl;
00908           spillsPerFile=800;
00909         }
00910         else {
00911           //default to Carrot value
00912           potPerSpill=25.23;
00913           MAXMSG("NuBeam",Msg::kWarning,1)
00914             <<"Defaulting to potPerSpill="<<potPerSpill
00915             <<" for unknown MC"<<endl;
00916           spillsPerFile=400;
00917         }
00918       }
00919       else {
00920         MAXMSG("NuBeam",Msg::kError,10)
00921           <<"Ahhh, don't recognise hornCurrent="
00922           <<config.hornCurrent<<endl;
00923         assert(false);
00924       }
00925     }
00926     else if (config.beamType==BeamType::kME){
00927       if (config.mcVersion==ReleaseType::kCarrot) {
00928         potPerSpill=10.17;//was 10.34;
00929         MAXMSG("NuBeam",Msg::kInfo,1)
00930           <<"Using potPerSpill="<<potPerSpill<<" for Carrot"<<endl;
00931         spillsPerFile=500;
00932       }
00933       else if (config.mcVersion==ReleaseType::kDaikon) {
00934         potPerSpill=11.40;//was 10.34
00935         MAXMSG("NuBeam",Msg::kInfo,1)
00936           <<"Using potPerSpill="<<potPerSpill<<" for Daikon"<<endl;
00937         spillsPerFile=1000;
00938       }
00939       else {
00940         //default to Carrot value
00941         potPerSpill=10.17;
00942         MAXMSG("NuBeam",Msg::kWarning,1)
00943           <<"Defaulting to potPerSpill="<<potPerSpill
00944           <<" for unknown MC"<<endl;
00945         spillsPerFile=500;
00946       }
00947     }
00948     else if (config.beamType==BeamType::kHE){
00949       if (config.mcVersion==ReleaseType::kCarrot) {
00950         potPerSpill=11.14;//was 11.41;
00951         MAXMSG("NuBeam",Msg::kInfo,1)
00952           <<"Using potPerSpill="<<potPerSpill<<" for Carrot"<<endl;
00953         spillsPerFile=130;
00954       }
00955       else if (config.mcVersion==ReleaseType::kDaikon) {
00956         potPerSpill=11.40;
00957         MAXMSG("NuBeam",Msg::kInfo,1)
00958           <<"Using potPerSpill="<<potPerSpill<<" for Daikon"<<endl;
00959         spillsPerFile=1000;
00960       }
00961       else {
00962         //default to Carrot value
00963         potPerSpill=11.14;
00964         MAXMSG("NuBeam",Msg::kWarning,1)
00965           <<"Defaulting to potPerSpill="<<potPerSpill
00966           <<" for unknown MC"<<endl;
00967         spillsPerFile=130;
00968       }
00969     }
00970     else {
00971       MAXMSG("NuBeam",Msg::kError,10)
00972         <<"Ahhh, don't know beamType="<<config.beamType<<endl;
00973       assert(false);
00974     }
00975   } 
00976   else if (config.detector==Detector::kFar){
00977     //can't fill hSpillsPerFile properly because the snarls with no 
00978     //detector activity are not written out...
00979     potPerSpill=2.90619;//e20: this is actually the number per file
00980     spillsPerFile=1;
00981   }
00982   else {
00983     MAXMSG("NuBeam",Msg::kError,10)
00984       <<"Ahhhh, don't know detector="<<config.detector<<endl;
00985     assert(false);
00986   }
00987 }
00988 */
00989 
00990 //......................................................................