loon: DataStabilityAnalysisClass.cxx Source File

00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 #include "TH1.h"
00009 #include "TH1D.h"
00010 #include "TH2D.h"
00011 #include "TProfile2D.h"
00012 #include "TMath.h"
00013 #include "TDirectory.h"
00014 #include "TROOT.h"
00015 #include <iomanip>
00016 #include "TChain.h"
00017 #include "TChainElement.h"
00018 #include "TLeaf.h"
00019 #include <math.h>
00020 #include <fstream>
00021 #include <sstream>
00022 #include <string>
00023 #include <iostream>
00024 #include <iomanip>
00025 
00026 #include "Conventions/Detector.h"
00027 #include "Conventions/SimFlag.h"
00028 #include "MessageService/MsgService.h"
00029 #include "NtupleUtils/NuCutter.h"
00030 #include "NtupleUtils/NuEvent.h"
00031 #include "NtupleUtils/NuInputEvents.h"
00032 #include "NtupleUtils/NuLibrary.h"
00033 #include "NtupleUtils/NuSystematic.h"
00034 #include "NuMuBar/NuTransition.h"
00035 #include "NtupleUtils/NuXMLConfig.h"
00036 #include "NuMuBar/DataStabilityAnalysisClass.h"
00037 #include "Validity/VldTimeStamp.h"
00038 #include "BeamDataUtil/BDSpillAccessor.h"
00039 #include "BeamDataUtil/BeamMonSpill.h"
00040 #include "Util/MsgUtil.h"
00041 
00042 using std::vector;
00043 
00044 CVSID("$Id: DataStabilityAnalysisClass.cxx,v 1.6 2015/09/30 15:02:09 ashley90 Exp $");
00045 
00047 DataStabilityAnalysisClass::DataStabilityAnalysisClass(NuXMLConfig* xmlConfig)
00048 {
00049   if (!xmlConfig){
00050     MSG("DataStabilityAnalysisClass",Msg::kWarning)
00051       << "Analysis object configured with a non-existent xml object."
00052       << endl
00053       << "Asserting now."
00054       << endl;
00055     assert(false);
00056   }
00057   fxmlConfig = xmlConfig;
00058 
00059   // Get the anaVersion and anaVersionNC from the XML config
00060   TString anaVer = fxmlConfig->AnaVersionString();
00061 
00062   fCutter = new NuCutter(anaVer);
00063   timeTree = new TTree("timeTree","Tree of times");
00064 
00065   // Get NuAnaVersion_t from the cutter objects
00066   fAnaVersion = static_cast<NuCuts::NuAnaVersion_t>(fCutter->AnaVersion());
00067   
00068   //create NuSystematic object:
00069   fsyst = new NuSystematic(*xmlConfig);
00070   fsyst->SetNuBarSelector(fCutter->GetCut());
00071   fsyst->SetCCSelector(fCutter->GetCut());
00072   fsyst->SetNCSelector(fCutter->GetCut());
00073   
00074   fLib = &NuLibrary::Instance();
00075   
00076   // Intialisating a few things
00077   prevSnarl = 0.0;
00078   tmpSubRun = 0;
00079   tmprun = 0;
00080   newSubRun = false;
00081 
00082   // Fixed event time
00083   // i cant use nu->timeSeconds 
00084   // because i only know the POT to the entry(0)
00085   FixedTime = 0;
00086 
00087   // MINOS+ start run time (in unix)
00088   //StartTime = 1378080000; // this is 2nd Sep 2013
00089   //EndTime = 1693612800; // this is 1st Sep 2023
00090 
00091   // MINOS start time (in unix)
00092   StartTime = 1062547200; // this is 3rd Sep 2003
00093   EndTime = 1378080000; // this is 2nd Sep 2013
00094 
00095   fHistogramsCreated = false;
00096   fTreeCreated = false;
00097 
00098   DSTpotCounter = 0.0;
00099   DSTspillcounter = 0;
00100   delta = 0;
00101 
00102   // initalise tree variables
00103   tselectionevent = -999;
00104   ttimeRelToSpill = -999;
00105   ttimeSec = -999;
00106   ttimeNanoSec = -999;
00107   tnearestSpillSec = -999;
00108   tnearestSpillNanosec = -999;
00109   tnshw = -999;
00110   trun = -999;
00111   tsubRun = -999;
00112   tsnarl = -999;
00113   ttrkEn = -999.;
00114   ttrkEnRange = -999.;
00115   ttrkEnCurv = -999.;
00116   tshwEn = -999.;
00117   tshwEnLinCCCor = -999.;
00118   tshwEnkNN = -999.;
00119   txEvtVtx = -999.;
00120   tyEvtVtx = -999.;
00121   tzEvtVtx = -999.;
00122   fTor101 = -999.;
00123   fTr101d = -999.;
00124   fTortgt = -999.;
00125   fTrtgtd = -999.;
00126   fHornCur = -999.;
00127   troID = -999.;
00128   fTargProfX = -999.;
00129   fTargProfY = -999.;
00130   fProfWidX = -999.;
00131   fProfWidY = -999.;
00132   fHadInt = 0.;
00133   fMuInt1 = 0.;
00134   fMuInt2 = 0.;
00135   fMuInt3 = 0.;
00136   totBpm =  0.;
00137   whichBatch = -1;
00138   numBatches = 0;
00139   //batchPot[6] = {0.,0.,0.,0.,0.,0.};
00140   beamTimeSec = 0.;
00141   beamTimeNSec = 0.;
00142 
00143   batchBin1summer = 0.0;
00144 }
00145 
00147 DataStabilityAnalysisClass::~DataStabilityAnalysisClass()
00148 {
00149   if (fCutter) {
00150     delete fCutter;
00151     fCutter = 0;
00152   }
00153   fLib = 0; //I didn't create this, so don't have the right to delete it.
00154   this->DeleteHistograms();
00155 }
00157 NuInputEvents& DataStabilityAnalysisClass::GetEventListForLooping()
00158 {
00159   NuInputEvents& input = this->DoIO(0, "", "null");
00160   fxmlConfig->Write();
00161   input.ResetNuEventLoopPositionToStart();
00162   return input;
00163 }
00164 
00166 void DataStabilityAnalysisClass::MakePlots()
00167 {
00168   NuInputEvents& input = this->GetEventListForLooping();
00169 
00170   BDSpillAccessor& bsAccess = BDSpillAccessor::Get();
00171   BeamMonSpill *beamSpill = 0x0;
00172 
00173   for (Int_t i = 0; i < input.GetEntriesNuEvent(); ++i){
00174 
00175     NuUtilities::ProgressBar(i, input.GetEntriesNuEvent(), 5);
00176     NuEvent &nu = const_cast<NuEvent&>(input.GetNextNuEvent(Msg::kDebug));
00177 
00178     if (!fHistogramsCreated){ this->CreateHistograms(nu); }
00179     if (!fTreeCreated){ this->CreateTree(nu); }
00180 
00181     this->PrepareEventForPlotting(nu);
00182 
00183     // I want to fix the time of all neutrinos to the nu->timeSeconds of the first 
00184     // event in each subrun, so that it matches up with the PotPerTime histogram
00185     //if(nu.subRun != tmpSubRun || nu.run != tmprun){ // uncomment if doing on subrun basis
00186       if(nu.run != tmprun){ // umcomment it doing only rn by run basis
00187       FixedTime = nu.timeSeconds;
00188     }
00189 
00190     // comment this function out if only doing DSTs run by run not subrun bu subrun
00191     //if(nu.subRun != tmpSubRun){
00192     //tmpSubRun = nu.subRun;
00193     //newSubRun = true;
00194     //}
00195 
00196     if(nu.run != tmprun){
00197       tmprun = nu.run;
00198     }
00199 
00200     // Now lets make some cuts
00201     fCutter->MakeCuts(nu);
00202 
00203     // the tree is created because it matches events with data base infomation                                                                                              
00204     // this is useful for pot per batch studies 
00205     Bool_t passedSelectionCut = false; 
00206 
00207     // if true then its a selection event
00208     if( fCutter->AnaVersion() == NuCuts::kNCPRLNC    ||
00209         fCutter->AnaVersion() == NuCuts::kCC1070Std  ||
00210         fCutter->AnaVersion() == NuCuts::kCC0720Std ){
00211       
00212       if(fCutter->Passed() == true){
00213         passedSelectionCut = true;
00214       }
00215       else{
00216         passedSelectionCut = false;
00217       }
00218     }
00219     if( fCutter->AnaVersion() == NuCuts:: kNMB0720Bravo ){
00220       if(fCutter->Passed() == true && nu.charge == +1){
00221         passedSelectionCut = true;
00222       } 
00223       else{
00224         passedSelectionCut = false;
00225       }
00226     }
00227 
00228     // Fill the tree for all events.. well i will cut on only good events inside.
00229     //this->FillTTree(nu, bsAccess, beamSpill, passedSelectionCut);
00230 
00231     // If it passed
00232     if ( fCutter->Passed() ){ 
00233       
00234       // Check if its CC, NC or Preselection
00235       if( fCutter->AnaVersion() == NuCuts::kNCPRLNC    || 
00236           fCutter->AnaVersion() == NuCuts::kCC1070Std  || 
00237           fCutter->AnaVersion() == NuCuts::kCC0720Std ){
00238 
00239         // the tree is created because it matches events with data base infomation
00240         // this is useful for pot per batch studies
00241 
00242         // this just fills a set of histograms
00243         this->FillPlots(nu);
00244       }
00245 
00246       // Check if its NuMubar
00247       if( fCutter->AnaVersion() == NuCuts:: kNMB0720Bravo ){
00248 
00249         // NuMuBar Selector doesnt cut on charge so i will.
00250           if(nu.charge == +1){
00251             this->FillPlots(nu);
00252 
00253           }
00254       }
00255     }
00256     if (fCutter->Failed()){
00257     }
00258   }
00259 
00260   this->WritePlots();
00261   fCutter->PrintSummary();
00262   fCutter->PrintNMinusOneSummary();
00263   fOutFile->Close();
00264 }
00265 
00267 void DataStabilityAnalysisClass::PrepareEventForPlotting(NuEvent& event) const
00268 {
00269     event.anaVersion = fAnaVersion;
00270     fxmlConfig->ConfigureWeights(event);
00271 
00272     //RE-RECONSTRUCT the neutrino energy
00273     //fLib->reco.GetEvtEnergy(event, false);
00274     //fLib->reco.ApplyReweights(event);
00275     //fLib->reco.CalcResolution(event);
00276    
00277     //apply systematics & oscillations
00278     if (event.simFlag!=SimFlag::kData) {
00279       fsyst->Shift(event);
00280       //Apply oscillations here if desired
00281     }
00282     
00283   //Used for NC selection comparisons 
00284     if(NuCuts::kNCPRLNC == fAnaVersion){
00285       //Do nothing for now
00286     }
00287     
00288 }
00290 void DataStabilityAnalysisClass::CreateTree(const NuEvent& event)
00291 {
00292   std::cout << "Creating TTree." << std::endl;
00293   
00294   // This suppresses most output when accessing database, expect fatal errors
00295   Util::MsgLevel("Dbi","Fatal");
00296   Util::MsgLevel("BDU","Fatal");
00297 
00298   if (event.simFlag == SimFlag::kMC){
00299     // do something with MC
00300   }
00301 
00302   // event stuff
00303   timeTree->Branch("nshw", &tnshw, "nshw/I");
00304   timeTree->Branch("selectionevent", &tselectionevent, "selectionevent/I");
00305   timeTree->Branch("timeRelToSpill",&ttimeRelToSpill,"timeRelToSpill/D");
00306   timeTree->Branch("timeSec",&ttimeSec,"timeSec/I");
00307   timeTree->Branch("timeNanoSec",&ttimeNanoSec,"timeNanoSec/I");
00308   timeTree->Branch("nearestSpillSec",&tnearestSpillSec,"nearestSpillSec/I");
00309   timeTree->Branch("nearestSpillNanosec",&tnearestSpillNanosec,"nearestSpillNanosec/I");;
00310   timeTree->Branch("run",&trun,"run/I");
00311   timeTree->Branch("subRun",&tsubRun,"subRun/I");
00312   timeTree->Branch("snarl",&tsnarl,"snarl/I");
00313   timeTree->Branch("trkEn",&ttrkEn,"trkEn/F");
00314   timeTree->Branch("trkEnRange",&ttrkEnRange,"trkEnRange/F");
00315   timeTree->Branch("trkEnCurv",&ttrkEnCurv,"trkEnCurv/F");
00316   timeTree->Branch("shwEn",&tshwEn,"shwEn/F");
00317   timeTree->Branch("shwEnLinCCCor",&tshwEnLinCCCor,"shwEnLinCCCor/F");
00318   timeTree->Branch("shwEnkNN",&tshwEnkNN,"shwEnkNN/F");
00319   timeTree->Branch("roID", &troID, "roID/F");
00320   timeTree->Branch("xEvtVtx",&txEvtVtx,"xEvtVtx/F");
00321   timeTree->Branch("yEvtVtx",&tyEvtVtx,"yEvtVtx/F");
00322   timeTree->Branch("zEvtVtx",&tzEvtVtx,"zEvtVtx/F");
00323 
00324   // beam stuff
00325   timeTree->Branch("batchPot", batchPot,"batchPot[6]/F");
00326   timeTree->Branch("numBatches",&numBatches,"numBatches/I");
00327   timeTree->Branch("whichBatch",&whichBatch,"whichBatch/I");
00328   timeTree->Branch("spillTimeSec",&beamTimeSec,"spillTimeSec/l");
00329   timeTree->Branch("spillTimeNSec",&beamTimeNSec,"spillTimeNSec/l");
00330   timeTree->Branch("fTor101",&fTor101,"fTor101/F");
00331   timeTree->Branch("fTr101d",&fTr101d,"fTr101d/F");
00332   timeTree->Branch("fTrtgtd",&fTrtgtd,"fTrtgtd/F");
00333   timeTree->Branch("fTortgt",&fTortgt,"fTortgt/F");
00334   timeTree->Branch("fHornCur",&fHornCur,"fHornCur/F");
00335   timeTree->Branch("fTargBpmX",fTargBpmX,"fTargBpmX[6]/F");
00336   timeTree->Branch("fTargBpmY",fTargBpmY,"fTargBpmY[6]/F");
00337   timeTree->Branch("fBpmInt",fBpmInt,"fBpmInt[6]/F");
00338   timeTree->Branch("fTargProfX",&fTargProfX,"fTargProfX/F");
00339   timeTree->Branch("fTargProfY",&fTargProfY,"fTargProfY/F");
00340   timeTree->Branch("fProfWidX",&fProfWidX,"fProfWidX/F");
00341   timeTree->Branch("fProfWidY",&fProfWidY,"fProfWidY/F");
00342   timeTree->Branch("fHadInt",&fHadInt,"fHadInt/F");
00343   timeTree->Branch("fMuInt1",&fMuInt1,"fMuInt1/F");
00344   timeTree->Branch("fMuInt2",&fMuInt2,"fMuInt2/F");
00345   timeTree->Branch("fMuInt3",&fMuInt3,"fMuInt3/F");
00346 
00347   fTreeCreated = true;
00348 }
00350 void DataStabilityAnalysisClass::CreateHistograms(const NuEvent& event)
00351 {
00352   std::cout << "Creating histograms" << std::endl;
00353 
00354   // Define binning scheme
00355   const vector<Double_t> recoBins = NuUtilities::RecoBins(NuBinningScheme::kSterile);
00356   const Int_t numRecoBins = recoBins.size() - 1;
00357   Double_t* recoBinsArray = new Double_t[numRecoBins + 1];
00358   Int_t i = 0;
00359 
00360   // Not too sure why its being converted into an array... but i will go with the flow.
00361   for (vector<Double_t>::const_iterator itBin = recoBins.begin();
00362        itBin != recoBins.end();
00363        ++itBin, ++i){
00364 
00365     recoBinsArray[i] = *itBin;
00366   }
00367 
00368   hRecoEnergyVsTime = new TH2D("hRecoEnergyVsTime","",3652, StartTime, EndTime, 400, 0, 200);
00369 
00370   // Neutrino Energy histograms
00371   //hRecoEnergyAll = new TH1D("hRecoEnergyAll","", 400, 0, 200);
00372   hRecoEnergyAll = new TH1D("hRecoEnergyAll","",numRecoBins, recoBinsArray);
00373   hRecoEnergyNQ = new TH1D("hRecoEnergyNQ","", 400, 0, 200);
00374   hRecoEnergyPQ = new TH1D("hRecoEnergyPQ","", 400, 0, 200);
00375 
00376   hEventsPerTimePerPOTLT3GeV = new TH1D("hEventsPerTimePerPOTLT3GeV", "", 3652, StartTime, EndTime);
00377   hEventsPerTimePerPOTGT3GeV = new TH1D("hEventsPerTimePerPOTGT3GeV", "", 3652, StartTime, EndTime);
00378   hEventsPerTimePerPOTLT6GeV = new TH1D("hEventsPerTimePerPOTLT6GeV", "", 3652, StartTime, EndTime);
00379   hEventsPerTimePerPOTGT6GeV = new TH1D("hEventsPerTimePerPOTGT6GeV", "", 3652, StartTime, EndTime);
00380   hEventsPerTimePerPOTLT8GeV = new TH1D("hEventsPerTimePerPOTLT8GeV", "", 3652, StartTime, EndTime);
00381   hEventsPerTimePerPOTGT8GeV = new TH1D("hEventsPerTimePerPOTGT8GeV", "", 3652, StartTime, EndTime);
00382   hEventsPerTimePerPOTLT14GeV = new TH1D("hEventsPerTimePerPOTLT14GeV", "", 3652, StartTime, EndTime);
00383   hEventsPerTimePerPOTGT14GeV = new TH1D("hEventsPerTimePerPOTGT14GeV", "", 3652, StartTime, EndTime);
00384 
00385   // Muon Energy Histograms
00386   hRecoEnergyAllmuon = new TH1D("hRecoEnergyAllmuon","",400,0,200);
00387   hEventsPerTimePerPOTLT3GeVmuon = new TH1D("hEventsPerTimePerPOTLT3GeVmuon", "", 3652, StartTime, EndTime);
00388   hEventsPerTimePerPOTGT3GeVmuon = new TH1D("hEventsPerTimePerPOTGT3GeVmuon", "", 3652, StartTime, EndTime);
00389   hEventsPerTimePerPOTLT6GeVmuon = new TH1D("hEventsPerTimePerPOTLT6GeVmuon", "", 3652, StartTime, EndTime);
00390   hEventsPerTimePerPOTGT6GeVmuon = new TH1D("hEventsPerTimePerPOTGT6GeVmuon", "", 3652, StartTime, EndTime);
00391   hEventsPerTimePerPOTLT8GeVmuon = new TH1D("hEventsPerTimePerPOTLT8GeVmuon", "", 3652, StartTime, EndTime);
00392   hEventsPerTimePerPOTGT8GeVmuon = new TH1D("hEventsPerTimePerPOTGT8GeVmuon", "", 3652, StartTime, EndTime);
00393   hEventsPerTimePerPOTLT14GeVmuon = new TH1D("hEventsPerTimePerPOTLT14GeVmuon", "", 3652, StartTime, EndTime);
00394   hEventsPerTimePerPOTGT14GeVmuon = new TH1D("hEventsPerTimePerPOTGT14GeVmuon", "", 3652, StartTime, EndTime);
00395 
00396   // Batch pot                                                                                                                                                                   
00397   hbatchPOT_1batchmode = new TH2D( "batchPOT_1batchmode", "batchPOT_1batchmode; Batch Index; Batch Intensity (x10^{12} POT)", 7, 0, 7, 200, 0, 20);
00398   hbatchPOT_2batchmode = new TH2D( "batchPOT_2batchmode", "batchPOT_2batchmode; Batch Index; Batch Intensity (x10^{12} POT)", 7, 0, 7, 200, 0, 20);
00399   hbatchPOT_3batchmode = new TH2D( "batchPOT_3batchmode", "batchPOT_3batchmode; Batch Index; Batch Intensity (x10^{12} POT)", 7, 0, 7, 200, 0, 20);
00400   hbatchPOT_4batchmode = new TH2D( "batchPOT_4batchmode", "batchPOT_4batchmode; Batch Index; Batch Intensity (x10^{12} POT)", 7, 0, 7, 200, 0, 20);
00401   hbatchPOT_5batchmode = new TH2D( "batchPOT_5batchmode", "batchPOT_5batchmode; Batch Index; Batch Intensity (x10^{12} POT)", 7, 0, 7, 200, 0, 20);
00402   hbatchPOT_6batchmode = new TH2D( "batchPOT_6batchmode", "batchPOT_6batchmode; Batch Index; Batch Intensity (x10^{12} POT)", 7, 0, 7, 200, 0, 20);
00403 
00404   hBatchPOTSum_1batchmode = new TH1D( "hBatchPOTSum_1batchmode","hBatchPOTSum_1batchmode; Batch Index; Batch POT (x10^{12} POT) ", 7, 0, 7);
00405   hBatchPOTSum_2batchmode = new TH1D( "hBatchPOTSum_2batchmode","hBatchPOTSum_2batchmode; Batch Index; Batch POT (x10^{12} POT) ", 7, 0, 7);
00406   hBatchPOTSum_3batchmode = new TH1D( "hBatchPOTSum_3batchmode","hBatchPOTSum_3batchmode; Batch Index; Batch POT (x10^{12} POT) ", 7, 0, 7);
00407   hBatchPOTSum_4batchmode = new TH1D( "hBatchPOTSum_4batchmode","hBatchPOTSum_4batchmode; Batch Index; Batch POT (x10^{12} POT) ", 7, 0, 7);
00408   hBatchPOTSum_5batchmode = new TH1D( "hBatchPOTSum_5batchmode","hBatchPOTSum_5batchmode; Batch Index; Batch POT (x10^{12} POT) ", 7, 0, 7);
00409   hBatchPOTSum_6batchmode = new TH1D( "hBatchPOTSum_6batchmode","hBatchPOTSum_6batchmode; Batch Index; Batch POT (x10^{12} POT) ", 7, 0, 7);
00410 
00411   // selection event batches histogram
00412   hSelectionBatchesAll         = new TH1D("hSelectionBatchesAll"        , "", 4000, 0, 12e-6); // in units of ns
00413   hSelectionBatches_BatchMode6 = new TH1D("hSelectionBatches_BatchMode6", "", 4000, 0, 12e-6);
00414   hSelectionBatches_BatchMode5 = new TH1D("hSelectionBatches_BatchMode5", "", 4000, 0, 12e-6);
00415   hSelectionBatches_BatchMode4 = new TH1D("hSelectionBatches_BatchMode4", "", 4000, 0, 12e-6);
00416   hSelectionBatches_BatchMode3 = new TH1D("hSelectionBatches_BatchMode3", "", 4000, 0, 12e-6);
00417   hSelectionBatches_BatchMode2 = new TH1D("hSelectionBatches_BatchMode2", "", 4000, 0, 12e-6);
00418   hSelectionBatches_BatchMode1 = new TH1D("hSelectionBatches_BatchMode1", "", 4000, 0, 12e-6);
00419 
00420   // Calibrate POT histogram
00421   hCalibrateBatchPOT = new TH2D("hCalibrateBatchPOT","CalibrateBatchPOT; fTrtgtd (x10^{12} POT); #Sigma fBpmInt", 250, 0, 40, 500, 0, 4000);
00422 
00423   // Special histogram - look at z vertex as function of energy
00424   hVertexZvsEnergy = new TH2D("hVertexZvsEnergy","",400,0,200,100,0,8);
00425 
00426   // For overall spill
00427   hPotPerIntensity = new TH1D("hPotPerIntensity", "hPotPerIntensity; Intensity per Spill (x10^{12} POT); Summed POT x10^{12}", 450, 0, 50);// pot per intensity
00428  
00429   // Now just for batches
00430   hPotPerIntensity_5BatchMode_Batch[0] = new TH1D("hPotPerIntensity_5BatchMode_Batch1", 
00431                                                   "hPotPerIntensity_5BatchMode_Batch1; Intensity per Batch (x10^{12} POT); Summed POT x10^{12}", 400, 0, 20);
00432   hPotPerIntensity_5BatchMode_Batch[1] = new TH1D("hPotPerIntensity_5BatchMode_Batch2",
00433                                                   "hPotPerIntensity_5BatchMode_Batch2; Intensity per Batch (x10^{12} POT); Summed POT x10^{12}", 400, 0, 20);
00434   hPotPerIntensity_5BatchMode_Batch[2] = new TH1D("hPotPerIntensity_5BatchMode_Batch3",
00435                                                   "hPotPerIntensity_5BatchMode_Batch3; Intensity per Batch (x10^{12} POT); Summed POT x10^{12}", 400, 0, 20);
00436   hPotPerIntensity_5BatchMode_Batch[3] = new TH1D("hPotPerIntensity_5BatchMode_Batch4",
00437                                                   "hPotPerIntensity_5BatchMode_Batch4; Intensity per Batch (x10^{12} POT); Summed POT x10^{12}", 400, 0, 20);
00438   hPotPerIntensity_5BatchMode_Batch[4] = new TH1D("hPotPerIntensity_5BatchMode_Batch5",
00439                                                   "hPotPerIntensity_5BatchMode_Batch5; Intensity per Batch (x10^{12} POT); Summed POT x10^{12}", 400, 0, 20);
00440 
00441   hPotPerIntensity_6BatchMode_Batch[0] = new TH1D("hPotPerIntensity_6BatchMode_Batch1",
00442                                                   "hPotPerIntensity_6BatchMode_Batch1; Intensity per Batch (x10^{12} POT); Summed POT x10^{12}", 400, 0, 20);
00443   hPotPerIntensity_6BatchMode_Batch[1] = new TH1D("hPotPerIntensity_6BatchMode_Batch2",
00444                                                   "hPotPerIntensity_6BatchMode_Batch2; Intensity per Batch (x10^{12} POT); Summed POT x10^{12}", 400, 0, 20);
00445   hPotPerIntensity_6BatchMode_Batch[2] = new TH1D("hPotPerIntensity_6BatchMode_Batch3",
00446                                                   "hPotPerIntensity_6BatchMode_Batch3; Intensity per Batch (x10^{12} POT); Summed POT x10^{12}", 400, 0, 20);
00447   hPotPerIntensity_6BatchMode_Batch[3] = new TH1D("hPotPerIntensity_6BatchMode_Batch4",
00448                                                   "hPotPerIntensity_6BatchMode_Batch4; Intensity per Batch (x10^{12} POT); Summed POT x10^{12}", 400, 0, 20);
00449   hPotPerIntensity_6BatchMode_Batch[4] = new TH1D("hPotPerIntensity_6BatchMode_Batch5",
00450                                                   "hPotPerIntensity_6BatchMode_Batch5; Intensity per Batch (x10^{12} POT); Summed POT x10^{12}", 400, 0, 20);
00451   hPotPerIntensity_6BatchMode_Batch[5] = new TH1D("hPotPerIntensity_6BatchMode_Batch6",
00452                                                   "hPotPerIntensity_6BatchMode_Batch6; Intensity per Batch (x10^{12} POT); Summed POT x10^{12}", 400, 0, 20);
00453 
00454 
00455   // number of spills with that intensity per day 
00456   IntensityPerTime = new TH2D("IntensityPerTime", "IntensityPerTime; Time; Spill Intensity (x10^{12} ", 3652, StartTime, EndTime, 450, 0, 50); 
00457   hNumberofGoodSpillsPerIntensity = new TH1D("NumberofGoodSpillsPerIntensity", "NumberofGoodSpillsPerIntensity; Intensity per Spill (x10^{12} POT); # Spills", 400, 0, 45);  // # of spills
00458 
00459   IntensityPerTime_5BatchMode_Batch[0] = new TH2D("IntensityPerTime_5BatchMode_Batch1",
00460                                                   "IntensityPerTime_5BatchMode_Batch1; Time; Spill Intensity (x10^{12} ", 3652, StartTime, EndTime, 400, 0, 20);
00461   IntensityPerTime_5BatchMode_Batch[1] = new TH2D("IntensityPerTime_5BatchMode_Batch2",
00462                                                   "IntensityPerTime_5BatchMode_Batch2; Time; Spill Intensity (x10^{12} ", 3652, StartTime, EndTime, 400, 0, 20);
00463   IntensityPerTime_5BatchMode_Batch[2] = new TH2D("IntensityPerTime_5BatchMode_Batch3",
00464                                                   "IntensityPerTime_5BatchMode_Batch3; Time; Spill Intensity (x10^{12} ", 3652, StartTime, EndTime, 400, 0, 20);
00465   IntensityPerTime_5BatchMode_Batch[3] = new TH2D("IntensityPerTime_5BatchMode_Batch4",
00466                                                   "IntensityPerTime_5BatchMode_Batch4; Time; Spill Intensity (x10^{12} ", 3652, StartTime, EndTime, 400, 0, 20);
00467   IntensityPerTime_5BatchMode_Batch[4] = new TH2D("IntensityPerTime_5BatchMode_Batch5",
00468                                                   "IntensityPerTime_5BatchMode_Batch5; Time; Spill Intensity (x10^{12} ", 3652, StartTime, EndTime, 400, 0, 20);
00469 
00470   IntensityPerTime_6BatchMode_Batch[0] = new TH2D("IntensityPerTime_6BatchMode_Batch1",
00471                                                   "IntensityPerTime_6BatchMode_Batch1; Time; Spill Intensity (x10^{12} ", 3652, StartTime, EndTime, 400, 0, 20);
00472   IntensityPerTime_6BatchMode_Batch[1] = new TH2D("IntensityPerTime_6BatchMode_Batch2",
00473                                                   "IntensityPerTime_6BatchMode_Batch2; Time; Spill Intensity (x10^{12} ", 3652, StartTime, EndTime, 400, 0, 20);
00474   IntensityPerTime_6BatchMode_Batch[2] = new TH2D("IntensityPerTime_6BatchMode_Batch3",
00475                                                   "IntensityPerTime_6BatchMode_Batch3; Time; Spill Intensity (x10^{12} ", 3652, StartTime, EndTime, 400, 0, 20);
00476   IntensityPerTime_6BatchMode_Batch[3] = new TH2D("IntensityPerTime_6BatchMode_Batch4",
00477                                                   "IntensityPerTime_6BatchMode_Batch4; Time; Spill Intensity (x10^{12} ", 3652, StartTime, EndTime, 400, 0, 20);
00478   IntensityPerTime_6BatchMode_Batch[4] = new TH2D("IntensityPerTime_6BatchMode_Batch5",
00479                                                   "IntensityPerTime_6BatchMode_Batch5; Time; Spill Intensity (x10^{12} ", 3652, StartTime, EndTime, 400, 0, 20);
00480   IntensityPerTime_6BatchMode_Batch[5] = new TH2D("IntensityPerTime_6BatchMode_Batch6",
00481                                                   "hPotPerIntensity_6BatchMode_Batch6; Time; Spill Intensity (x10^{12} ", 3652, StartTime, EndTime, 400, 0, 20);
00482 
00483   // NC Cleaning Histograms
00484   hedgeActivityPH     = new TH1D("hedgeActivityPH","",30,0,10000);
00485   hedgeActivityStrips = new TH1D("hedgeActivityStrips","",20,0,20);
00486   hminTimeSeparation  = new TH1D("hminTimeSeparation","",100,-1000,1000);
00487   hoppEdgePH          = new TH1D("hoppEdgePH","",30,0,5000);
00488   hoppEdgeStrips      = new TH1D("hoppEdgeStrips","",20,0,16);
00489   hcloseTimeDeltaZ    = new TH1D("hcloseTimeDeltaZ","",100,-10,10);
00490   hnstripEvt          = new TH1D("hnstripEvt","",100,0,300);
00491   hstripsOverSqPlanes = new TH1D("hstripsOverSqPlanes","",80,0,0.4);
00492 
00493   // Event Plotting Histograms
00494   hydistribution   = new TH1D("hydistribution","",200,0,1);
00495   hplanesEvent     = new TH1D("hplanesEvent","",200,0,200);
00496   hPhEvent         = new TH1D("hPhEvent","",300,0,6e5);
00497   hPhPerPlane      = new TH1D("hPhPerPlane","",200,0,1e4);
00498   hPhPerStrip      = new TH1D("hPhPerStrip","",200,0,1e4);
00499   hStripsPerPlane  = new TH1D("hStripsPerPlane","",15,0,15);
00500   hTracksInEvent   = new TH1D("hTracksInEvent","",8,0,8);  
00501   hShowersInEvent  = new TH1D("hShowersInEvent","",10,0,10);
00502   hEventXVertex    = new TH1D("hEventXVertex","",100,-1,3);
00503   hEventYVertex    = new TH1D("hEventYVertex","",100,-2,2);
00504   hEventZVertex    = new TH1D("hEventZVertex","",100,0,8);
00505   hEventRVertex    = new TH1D("hEventRVertex","",100,0,3);
00506 
00507   // Track Plotting Histograms
00508   hTrkLength          = new TH1D("hTrkLength","",280,0,280);
00509   hPlanesCrossedInTrk = new TH1D("hPlanesCrossedInTrk","",300,0,300);
00510   hTrkCosTheta        = new TH1D("hTrkCosTheta","",100,-1,1);
00511   hTrkdzds            = new TH1D("hTrkdzds","",100,0,1);
00512   hTrkPh              = new TH1D("hTrkPh","",50,0,20e4);
00513   hTrkMomentumRange   = new TH1D("hTrkMomentumRange","",400,0,100);
00514   hTrkPhPerPlane      = new TH1D("hTrkPhPerPlane","",20,0,2);
00515   hTrkqpSigma         = new TH1D("hTrkqpSigma","",400,0,100);
00516   hTrkPassedFit       = new TH1D("hTrkPassedFit","",2,0,2);
00517   hTrkChi2PerNdf      = new TH1D("hTrkChi2PerNdf","",100,0,6);
00518   hTrkXVertex         = new TH1D("hTrkXVertex","",50,-1,4);
00519   hTrkYVertex         = new TH1D("hTrkYVertex","",50,-2,3);
00520   hTrkZVertex         = new TH1D("hTrkZVertex","",200,0,8);
00521 
00522   // Shower Plotting Histograms
00523   hPlanesCrossedInShower = new TH1D("hPlanesCrossedInShower","",100,0,100);
00524   hShowerEnergy          = new TH1D("hShowerEnergy","",100,0,50);
00525   hXShwVtx               = new TH1D("hXShwVtx","",50,-1,4);
00526   hYShwVtx               = new TH1D("hYShwVtx","",50,-2,3);
00527   hZShwVtx               = new TH1D("hZShwVtx","",50,0,10);
00528   
00529   if (event.simFlag == SimFlag::kMC){
00530     //Create MC-only plots. You can similarly create ND-only plot, etc
00531     // Good for truth information that real data doesnt have.
00532   }
00533 
00534   fvHistos.push_back(hRecoEnergyVsTime);
00535   
00536   fvHistos.push_back(hRecoEnergyAll);
00537   fvHistos.push_back(hRecoEnergyNQ);
00538   fvHistos.push_back(hRecoEnergyPQ);
00539 
00540   // Batch pot                                                                                                                                    
00541   fvHistos.push_back(hbatchPOT_1batchmode);
00542   fvHistos.push_back(hbatchPOT_2batchmode);
00543   fvHistos.push_back(hbatchPOT_3batchmode);
00544   fvHistos.push_back(hbatchPOT_4batchmode);
00545   fvHistos.push_back(hbatchPOT_5batchmode);
00546   fvHistos.push_back(hbatchPOT_6batchmode);
00547 
00548   fvHistos.push_back(hCalibrateBatchPOT);
00549 
00550   fvHistos.push_back(hBatchPOTSum_1batchmode);
00551   fvHistos.push_back(hBatchPOTSum_2batchmode);
00552   fvHistos.push_back(hBatchPOTSum_3batchmode);
00553   fvHistos.push_back(hBatchPOTSum_4batchmode);
00554   fvHistos.push_back(hBatchPOTSum_5batchmode);
00555   fvHistos.push_back(hBatchPOTSum_6batchmode);
00556 
00557   fvHistos.push_back(hEventsPerTimePerPOTLT3GeV);
00558   fvHistos.push_back(hEventsPerTimePerPOTGT3GeV);
00559   fvHistos.push_back(hEventsPerTimePerPOTLT6GeV);
00560   fvHistos.push_back(hEventsPerTimePerPOTGT6GeV);
00561   fvHistos.push_back(hEventsPerTimePerPOTLT8GeV);
00562   fvHistos.push_back(hEventsPerTimePerPOTGT8GeV);
00563   fvHistos.push_back(hEventsPerTimePerPOTLT14GeV);
00564   fvHistos.push_back(hEventsPerTimePerPOTGT14GeV);
00565 
00566   // Muon Energy Histograms
00567   fvHistos.push_back(hRecoEnergyAllmuon);
00568   fvHistos.push_back(hEventsPerTimePerPOTLT3GeVmuon);
00569   fvHistos.push_back(hEventsPerTimePerPOTGT3GeVmuon);
00570   fvHistos.push_back(hEventsPerTimePerPOTLT6GeVmuon);
00571   fvHistos.push_back(hEventsPerTimePerPOTGT6GeVmuon);
00572   fvHistos.push_back(hEventsPerTimePerPOTLT8GeVmuon);
00573   fvHistos.push_back(hEventsPerTimePerPOTGT8GeVmuon);
00574   fvHistos.push_back(hEventsPerTimePerPOTLT14GeVmuon);
00575   fvHistos.push_back(hEventsPerTimePerPOTGT14GeVmuon);
00576 
00577   fvHistos.push_back(hSelectionBatchesAll);
00578   fvHistos.push_back(hSelectionBatches_BatchMode6);
00579   fvHistos.push_back(hSelectionBatches_BatchMode5);
00580   fvHistos.push_back(hSelectionBatches_BatchMode4);
00581   fvHistos.push_back(hSelectionBatches_BatchMode3);
00582   fvHistos.push_back(hSelectionBatches_BatchMode2);
00583   fvHistos.push_back(hSelectionBatches_BatchMode1);
00584   fvHistos.push_back(hVertexZvsEnergy);
00585 
00586   fvHistos.push_back(hPotPerIntensity);
00587 
00588   fvHistos.push_back(hPotPerIntensity_5BatchMode_Batch[0]);
00589   fvHistos.push_back(hPotPerIntensity_5BatchMode_Batch[1]);
00590   fvHistos.push_back(hPotPerIntensity_5BatchMode_Batch[2]);
00591   fvHistos.push_back(hPotPerIntensity_5BatchMode_Batch[3]);
00592   fvHistos.push_back(hPotPerIntensity_5BatchMode_Batch[4]);
00593 
00594   fvHistos.push_back(hPotPerIntensity_6BatchMode_Batch[0]);
00595   fvHistos.push_back(hPotPerIntensity_6BatchMode_Batch[1]);
00596   fvHistos.push_back(hPotPerIntensity_6BatchMode_Batch[2]);
00597   fvHistos.push_back(hPotPerIntensity_6BatchMode_Batch[3]);
00598   fvHistos.push_back(hPotPerIntensity_6BatchMode_Batch[4]);
00599   fvHistos.push_back(hPotPerIntensity_6BatchMode_Batch[5]);
00600 
00601   fvHistos.push_back(IntensityPerTime);
00602   
00603   fvHistos.push_back(IntensityPerTime_5BatchMode_Batch[0]);
00604   fvHistos.push_back(IntensityPerTime_5BatchMode_Batch[1]);
00605   fvHistos.push_back(IntensityPerTime_5BatchMode_Batch[2]);
00606   fvHistos.push_back(IntensityPerTime_5BatchMode_Batch[3]);
00607   fvHistos.push_back(IntensityPerTime_5BatchMode_Batch[4]);
00608 
00609   fvHistos.push_back(IntensityPerTime_6BatchMode_Batch[0]);
00610   fvHistos.push_back(IntensityPerTime_6BatchMode_Batch[1]);
00611   fvHistos.push_back(IntensityPerTime_6BatchMode_Batch[2]);
00612   fvHistos.push_back(IntensityPerTime_6BatchMode_Batch[3]);
00613   fvHistos.push_back(IntensityPerTime_6BatchMode_Batch[4]);
00614   fvHistos.push_back(IntensityPerTime_6BatchMode_Batch[5]);
00615 
00616   fvHistos.push_back(hNumberofGoodSpillsPerIntensity);
00617 
00618   fvHistos.push_back(hedgeActivityPH);
00619   fvHistos.push_back(hedgeActivityStrips);
00620   fvHistos.push_back(hminTimeSeparation);
00621   fvHistos.push_back(hoppEdgePH);
00622   fvHistos.push_back(hoppEdgeStrips);
00623   fvHistos.push_back(hcloseTimeDeltaZ);
00624   fvHistos.push_back(hnstripEvt);
00625   fvHistos.push_back(hstripsOverSqPlanes);
00626 
00627   fvHistos.push_back(hydistribution);
00628   fvHistos.push_back(hplanesEvent);
00629   fvHistos.push_back(hPhEvent);
00630   fvHistos.push_back(hPhPerPlane);
00631   fvHistos.push_back(hPhPerStrip);
00632   fvHistos.push_back(hStripsPerPlane);
00633   fvHistos.push_back(hTracksInEvent);
00634   fvHistos.push_back(hShowersInEvent);
00635   fvHistos.push_back(hEventXVertex);
00636   fvHistos.push_back(hEventYVertex);
00637   fvHistos.push_back(hEventZVertex);
00638   fvHistos.push_back(hEventRVertex);
00639 
00640   fvHistos.push_back(hTrkLength);
00641   fvHistos.push_back(hPlanesCrossedInTrk);
00642   fvHistos.push_back(hTrkCosTheta);
00643   fvHistos.push_back(hTrkdzds);
00644   fvHistos.push_back(hTrkPh);
00645   fvHistos.push_back(hTrkMomentumRange);
00646   fvHistos.push_back(hTrkPhPerPlane);
00647   fvHistos.push_back(hTrkqpSigma);
00648   fvHistos.push_back(hTrkPassedFit);
00649   fvHistos.push_back(hTrkChi2PerNdf);
00650   fvHistos.push_back(hTrkXVertex);
00651   fvHistos.push_back(hTrkYVertex);
00652   fvHistos.push_back(hTrkZVertex);
00653 
00654   fvHistos.push_back(hPlanesCrossedInShower);
00655   fvHistos.push_back(hShowerEnergy);
00656   fvHistos.push_back(hXShwVtx);
00657   fvHistos.push_back(hYShwVtx);
00658   fvHistos.push_back(hZShwVtx);
00659 
00660   fHistogramsCreated = true;
00661 }
00662 
00664 void DataStabilityAnalysisClass::DeleteHistograms()
00665 {
00666   for (vector<TH1*>::iterator it = fvHistos.begin();
00667        it != fvHistos.end();
00668        ++it){
00669     this->Delete(*it);
00670   }
00671 }
00672 
00674 void DataStabilityAnalysisClass::Delete(TH1* p)
00675 {
00676   if (p){delete p; p=0;}
00677 }
00678 
00680 void DataStabilityAnalysisClass::FillTTree(NuEvent& event, BDSpillAccessor& bsAccess, BeamMonSpill* beamSpill, Bool_t passedSelectionCut)
00681 {
00682 
00683   // only fill this tree with good events
00684   if( event.goodBeamSntp && event.coilIsOk && ( event.isGoodDataQuality || !event.useDBForDataQuality ) ){
00685 
00686     // Filling Tree
00687     ttimeRelToSpill      = (event.timeEvtMin + event.crateT0);
00688     ttimeSec             = event.timeSec;
00689     ttimeNanoSec         = event.timeNanoSec;
00690     tnearestSpillSec     = event.nearestSpillSec;
00691     tnearestSpillNanosec = event.nearestSpillNanosec;
00692     trun                 = event.run;
00693     tnshw                = event.nshw;
00694     tsubRun              = event.subRun;                
00695     tsnarl               = event.snarl;
00696     ttrkEn               = event.trkEn;
00697     ttrkEnRange          = event.trkEnRange;
00698     ttrkEnCurv           = event.trkEnCurv;
00699     tshwEn               = event.shwEn;
00700     tshwEnLinCCCor       = event.shwEnCor;
00701     troID                = event.roID;
00702 
00703     // if its a selection event... give it a one so i can cut on it later
00704     if(passedSelectionCut == true){
00705       tselectionevent = 1;  
00706     }
00707     else{
00708       tselectionevent = 0;
00709     }
00710     
00711     // get a time stamp from event
00712     VldTimeStamp timeStamp(event.timeSec,event.timeNanoSec);
00713     
00714     if( event.snarl != prevSnarl ) {
00715       
00716       DSTspillcounter++;
00717 
00718       beamSpill = const_cast<BeamMonSpill*>( bsAccess.LoadSpill(timeStamp) );
00719       VldTimeStamp vldBeamTime = beamSpill->SpillTime();
00720       
00721       beamTimeSec  = (ULong64_t)vldBeamTime.GetSec(); // Cast time_t to ULong64_t
00722       beamTimeNSec = (ULong64_t)vldBeamTime.GetNanoSec(); // nano_sec is Int_t by default 
00723       
00724       fTor101  = beamSpill->fTor101;
00725 
00726       fTr101d  = beamSpill->fTr101d;
00727       fTortgt  = beamSpill->fTortgt;
00728       fTrtgtd  = beamSpill->fTrtgtd;
00729 
00730       DSTpotCounter += fTrtgtd;
00731 
00732       fHornCur = beamSpill->fHornCur;
00733       
00734       //Per batch X,Y positions projected to the target via BPMs 121 and TGT.
00735       // memcpy - copies a whole array in one go, you have delcare size of course
00736       memcpy( fTargBpmX, beamSpill->fTargBpmX, 6*sizeof(float) );
00737       memcpy( fTargBpmY, beamSpill->fTargBpmY, 6*sizeof(float) );
00738       memcpy( fBpmInt,   beamSpill->fBpmInt,   6*sizeof(float) );
00739       
00740       fTargProfX = beamSpill->fTargProfX;
00741       fTargProfY = beamSpill->fTargProfY;
00742       fProfWidX  = beamSpill->fProfWidX;
00743       // Widths via profile monitors at TGT. No projection.                                                                                                     
00744       fProfWidY  = beamSpill->fProfWidY;
00745       fHadInt    = beamSpill->fHadInt;
00746       //  hadron/muon mon. intensities
00747       fMuInt1    = beamSpill->fMuInt1;
00748       fMuInt2    = beamSpill->fMuInt2;
00749       fMuInt3    = beamSpill->fMuInt3;
00750       totBpm     = fBpmInt[0] + fBpmInt[1] + fBpmInt[2] + fBpmInt[3] + fBpmInt[4] + fBpmInt[5];
00751       
00752       numBatches = 0;
00753 
00754       for( int batch = 0; batch < 6; batch++ ){
00755         //if( fBpmInt[batch] > 1e-20 ){
00756         if( fBpmInt[batch] > 0 ){
00757           numBatches++;
00758         }
00759         else{
00760           // Do nothing
00761         }
00762       }
00763 
00764       hCalibrateBatchPOT->Fill(fTrtgtd, totBpm);
00765 
00766       batchPot[0] = 0;
00767       batchPot[1] = 0;
00768       batchPot[2] = 0;
00769       batchPot[3] = 0;
00770       batchPot[4] = 0;
00771       batchPot[5] = 0;
00772     
00773       for(int batchnum = 0; batchnum < 6; ++batchnum){
00774 
00775         batchPot[batchnum] = fBpmInt[batchnum] * fTrtgtd / totBpm;
00776 
00777         if(numBatches == 1){
00778           hbatchPOT_1batchmode->Fill(batchnum, batchPot[batchnum]);
00779           hBatchPOTSum_1batchmode->Fill(batchnum, batchPot[batchnum]);
00780         }
00781         else if(numBatches == 2){
00782           hbatchPOT_2batchmode->Fill(batchnum, batchPot[batchnum]);
00783           hBatchPOTSum_2batchmode->Fill(batchnum, batchPot[batchnum]);
00784         }
00785         else if(numBatches == 3){
00786           hbatchPOT_3batchmode->Fill(batchnum, batchPot[batchnum]);
00787           hBatchPOTSum_3batchmode->Fill(batchnum, batchPot[batchnum]);
00788         }
00789         else if(numBatches == 4){
00790           hbatchPOT_4batchmode->Fill(batchnum, batchPot[batchnum]);
00791           hBatchPOTSum_4batchmode->Fill(batchnum, batchPot[batchnum]);
00792         }
00793         else if(numBatches == 5){
00794           hbatchPOT_5batchmode->Fill(batchnum, batchPot[batchnum]);
00795           hBatchPOTSum_5batchmode->Fill(batchnum, batchPot[batchnum]);
00796         }
00797         else if(numBatches == 6){
00798           hbatchPOT_6batchmode->Fill(batchnum, batchPot[batchnum]);
00799           hBatchPOTSum_6batchmode->Fill(batchnum, batchPot[batchnum]);
00800           if(batchnum == 0) {batchBin1summer += batchPot[batchnum];}
00801         }
00802         else{
00803           assert(false);
00804           //cout<< "numbatches = " << numBatches << ", this is wrong! asserting..." << endl; assert(false);
00805         }
00806       }
00807 
00808       cout << "snarl " << event.snarl << "Run number = " << event.run << ", POT 6 batch Bin1 = " << batchBin1summer << endl;
00809 
00810       // Now i get intensity profiles from this new snarl
00811       this->IntensityProfile(event, batchPot, fTrtgtd, numBatches);
00812 
00813 
00814       // cout some stuff
00815       //cout << " " << endl;
00816       //cout << "Snarl  = " << event.snarl << endl;
00817       //cout << "subRun = " << event.subRun << endl;
00818       //cout << "Run    = " << event.run << endl;
00819       //cout << "totBpm = " << totBpm << endl;
00820       //cout << "fTrtgtd= " << fTrtgtd << endl;
00821 
00822     } // end of new snarl loop
00823   
00824     // Fill some histograms
00825     if(numBatches == 6)       hSelectionBatches_BatchMode6->Fill((event.timeEvtMin + event.crateT0), event.rw);
00826     else if (numBatches == 5) hSelectionBatches_BatchMode5->Fill((event.timeEvtMin + event.crateT0), event.rw);
00827     else if (numBatches == 4) hSelectionBatches_BatchMode4->Fill((event.timeEvtMin + event.crateT0), event.rw);
00828     else if (numBatches == 3) hSelectionBatches_BatchMode3->Fill((event.timeEvtMin + event.crateT0), event.rw);
00829     else if (numBatches == 2) hSelectionBatches_BatchMode2->Fill((event.timeEvtMin + event.crateT0), event.rw);
00830     else if (numBatches == 1) hSelectionBatches_BatchMode1->Fill((event.timeEvtMin + event.crateT0), event.rw);
00831     else{
00832       // Do Nothing
00833     }
00834 
00835     // This is for every event
00836     whichBatch=-1;
00837     
00838     //These are the batches for MINOS+ RUN XI ad RUN XII
00839     if(ttimeRelToSpill > 1.28e-6   && ttimeRelToSpill <= 2.9e-6  ) whichBatch = 0;
00840     else if(ttimeRelToSpill > 2.9e-6  && ttimeRelToSpill <= 4.49e-6 ) whichBatch = 1;
00841     else if(ttimeRelToSpill > 4.49e-6 && ttimeRelToSpill <= 6.06e-6 ) whichBatch = 2;
00842     else if(ttimeRelToSpill > 6.06e-6 && ttimeRelToSpill <= 7.65e-6 ) whichBatch = 3;
00843     else if(ttimeRelToSpill > 7.65e-6 && ttimeRelToSpill <= 9.23e-6 ) whichBatch = 4;
00844     else if(ttimeRelToSpill > 9.23e-6 && ttimeRelToSpill <= 10.82e-6) whichBatch = 5;
00845 
00846     timeTree->Fill();
00847   }
00848   
00849   prevSnarl = event.snarl;
00850   
00851 }   
00853 void DataStabilityAnalysisClass::FillPlots(NuEvent& event)
00854 {
00855   this->PrepareEventForPlotting(event);
00856 
00857   Double_t energy = 0.0;
00858 
00859   // event.energy defaults to energyCC, i want to change this for the NC selector
00860 
00861   if(fCutter->AnaVersion() == NuCuts::kCC1070Std){
00862     energy = event.energy; 
00863   }
00864   if(fCutter->AnaVersion() == NuCuts::kNCPRLNC){
00865     energy = event.shwEnLinNCCor;
00866   }
00867 
00868   // Energy Graphs
00869   hRecoEnergyVsTime->Fill(FixedTime, energy);
00870   //hRecoEnergyAll->Fill( energy, event.rw);
00871   hRecoEnergyAll->Fill( energy, 1);
00872 
00873   if ( -1 == event.charge ){ hRecoEnergyNQ->Fill(energy, event.rw); }
00874   if (  1 == event.charge ){ hRecoEnergyPQ->Fill(energy, event.rw); }
00875 
00876   // These are the different energy regions
00877   if ( energy <  8.0 ){ hEventsPerTimePerPOTLT8GeV->Fill( FixedTime ); }
00878   if ( energy >= 8.0 ){ hEventsPerTimePerPOTGT8GeV->Fill( FixedTime ); }
00879   if ( energy <  6.0 ){ hEventsPerTimePerPOTLT6GeV->Fill( FixedTime ); }
00880   if ( energy >= 6.0 ){ hEventsPerTimePerPOTGT6GeV->Fill( FixedTime ); }
00881   if ( energy <  3.0 ){ hEventsPerTimePerPOTLT3GeV->Fill( FixedTime ); }
00882   if ( energy >= 3.0 ){ hEventsPerTimePerPOTGT3GeV->Fill( FixedTime ); }
00883   if ( energy <  14.0){ hEventsPerTimePerPOTLT14GeV->Fill( FixedTime); }
00884   if ( energy >= 14.0){ hEventsPerTimePerPOTGT14GeV->Fill( FixedTime); }
00885 
00886   // if CC then make muon graphs
00887   if(fCutter->AnaVersion()==29){
00888 
00889     if ( event.trkEn <  8.0){ hEventsPerTimePerPOTLT8GeVmuon->Fill( FixedTime ); }
00890     if ( event.trkEn >= 8.0){ hEventsPerTimePerPOTGT8GeVmuon->Fill( FixedTime ); }
00891     if ( event.trkEn <  6.0){ hEventsPerTimePerPOTLT6GeVmuon->Fill( FixedTime ); }
00892     if ( event.trkEn >= 6.0){ hEventsPerTimePerPOTGT6GeVmuon->Fill( FixedTime ); }
00893     if ( event.trkEn <  3.0){ hEventsPerTimePerPOTLT3GeVmuon->Fill( FixedTime ); }
00894     if ( event.trkEn >= 3.0){ hEventsPerTimePerPOTGT3GeVmuon->Fill( FixedTime ); }
00895     if ( event.trkEn <  14.0){hEventsPerTimePerPOTLT14GeVmuon->Fill( FixedTime); }
00896     if ( event.trkEn >= 14.0){hEventsPerTimePerPOTGT14GeVmuon->Fill( FixedTime); }
00897 
00898     hRecoEnergyAllmuon->Fill(event.trkEn, event.rw);
00899   }
00900 
00901   // Filling The Weekly Plot Histograms
00902   // NC Cleaning Variables
00903   hedgeActivityPH->Fill(    event.edgeActivityPH);
00904   hedgeActivityStrips->Fill(event.edgeActivityStrips);
00905   hminTimeSeparation->Fill( event.minTimeSeparation*1e9); // So It's in nanoseconds
00906   hoppEdgePH->Fill(         event.oppEdgePH);
00907   hoppEdgeStrips->Fill(     event.oppEdgeStrips);
00908   hcloseTimeDeltaZ->Fill(   event.closeTimeDeltaZ);
00909   hnstripEvt->Fill(         event.nstripEvt);
00910 
00911   Double_t stripsPerSqPlanes = 0;
00912   stripsPerSqPlanes = ((Double_t)event.nstripEvt)/(event.planeEvtN*event.planeEvtN);
00913   hstripsOverSqPlanes->Fill( stripsPerSqPlanes);
00914 
00915   // Eventhistograms
00916   hydistribution->Fill(event.y, event.rw);
00917   hplanesEvent->Fill(   event.planeEvtN);
00918   hPhEvent->Fill(       event.rawPhEvt);
00919   hPhPerPlane->Fill( ((Double_t)event.rawPhEvt/(Double_t)event.planeEvtN));
00920   hPhPerStrip->Fill( (Double_t)(event.rawPhEvt/(Double_t)event.nstripEvt));
00921   hStripsPerPlane->Fill((Double_t)(event.nstripEvt/event.planeEvtN));
00922   hTracksInEvent->Fill( event.ntrk);
00923   hShowersInEvent->Fill(event.nshw);
00924   hEventXVertex->Fill(  event.xEvtVtx);
00925   hEventYVertex->Fill(  event.yEvtVtx);
00926   hEventZVertex->Fill(  event.zEvtVtx);
00927 
00928   hSelectionBatchesAll->Fill((event.timeEvtMin + event.crateT0), event.rw);
00929   hVertexZvsEnergy    ->Fill(energy, event.zEvtVtx, event.rw);
00930 
00931   // Simply R = (x^2 + y^2)^0.5 as cylindrial coords
00932   Double_t RVertex = TMath::Sqrt((event.xEvtVtx * event.xEvtVtx) + 
00933                                  (event.yEvtVtx * event.yEvtVtx)); 
00934   hEventRVertex->Fill(RVertex);
00935 
00936   // Track Histograms1
00937   hTrkLength->Fill(         event.trkLength);
00938   hPlanesCrossedInTrk->Fill(event.trknplane);
00939   hTrkCosTheta->Fill(       event.dirCosNu);
00940   //hTrkCosTheta->Fill(     event.cosPrTrkVtx);
00941   hTrkdzds->Fill(           event.trkvtxdcosz);
00942   hTrkPh->Fill(             event.trkphsigcor);
00943   hTrkMomentumRange->Fill(  event.trkMomentumRange);
00944   hTrkPhPerPlane->Fill(     event.knn10TrkMeanPh);
00945   hTrkqpSigma->Fill(        event.qp_sigqp);
00946   hTrkPassedFit->Fill(      event.trkfitpass);
00947   hTrkChi2PerNdf->Fill(     event.chi2PerNdof);
00948   hTrkXVertex->Fill(        event.xTrkVtx);
00949   hTrkYVertex->Fill(        event.yTrkVtx);
00950   hTrkZVertex->Fill(        event.zTrkVtx);
00951 
00952   // Shower Histograms
00953   hPlanesCrossedInShower->Fill(event.nplaneShw);
00954   hShowerEnergy->Fill(event.shwEn);
00955   hXShwVtx->Fill(event.xShwVtx);
00956   hYShwVtx->Fill(event.yShwVtx);
00957   hZShwVtx->Fill(event.zShwVtx);
00958 
00959 }
00960 
00962 void DataStabilityAnalysisClass::WritePlots()
00963 {
00964   for (vector<TH1*>::iterator it = fvHistos.begin();
00965        it != fvHistos.end();
00966        ++ it){
00967     (*it)->Write();
00968   }
00969 
00970   // Write the tree
00971   timeTree->Write("timeTree");
00972 
00973   // Its nice to cout this at the end.
00974   std::cout << " " << std::endl;
00975   std::cout << "Spill POT     = " << DSTpotCounter << std::endl;
00976   std::cout << "Spill counter = " << DSTspillcounter  << std::endl;
00977   std::cout << " " << std::endl;
00978 }
00980 void DataStabilityAnalysisClass::IntensityProfile(NuEvent& nu, Double_t* batchPot, Double_t intensity, Int_t numBatches)
00981 {
00982   // This function is only called once per spill
00983 
00984   // if i want to resolve stuff to within a day - FixedTime is (on purpose) only good for resolving days
00985   //Double_t AccurateTime = nu.timeSeconds;
00986 
00987   // Recording the intensities for batches and overall spill
00988   hPotPerIntensity               ->Fill( intensity, intensity );
00989   hNumberofGoodSpillsPerIntensity->Fill( intensity );
00990   IntensityPerTime               ->Fill( FixedTime, intensity, 1); // where 1 is number of spills
00991 
00992   // Fill batch plots depending on batch mode 5 or 6 
00993   if(numBatches == 5){
00994     for(int i = 0; i < 5; ++i){
00995       hPotPerIntensity_5BatchMode_Batch[i]->Fill(batchPot[i]);
00996       IntensityPerTime_5BatchMode_Batch[i]->Fill(FixedTime, batchPot[i], 1);
00997     }
00998   }
00999   else if(numBatches == 6){
01000     for(int i = 0; i < 6; ++i){
01001       hPotPerIntensity_6BatchMode_Batch[i]->Fill(batchPot[i]);
01002       IntensityPerTime_6BatchMode_Batch[i]->Fill(FixedTime, batchPot[i], 1);
01003     }
01004   }
01005   else{
01006     // Do nohting
01007   }
01008 
01009 }
01011 void DataStabilityAnalysisClass::ObtainPotPerRun()
01012 {
01013   std::cout << "\n----------------------Obtaining Pot Per Run-----------------" << std::endl;
01014 
01015   string inputFileName = this->GetInputFileName();
01016 
01017   std::cout << "The InputFilename is: " << this->GetInputFileName() << endl;
01018 
01019   std::cout << "Creating the TChain" << std::endl;
01020 
01021   // connect them via a common tree, 's'
01022   TChain tChain("s");
01023   tChain.Add(inputFileName.c_str());
01024   TObjArray* fileList = tChain.GetListOfFiles();
01025 
01026   TH1D* hRunPot = new TH1D("hRunPot", "hRunPot; Run; Pot", 20000, 7000, 27000);
01027   TH1D* hGoodRunSpills = new TH1D("hGoodRunSpills", "hGoodRunSpills; Run; # Good Spills", 20000, 7000, 27000);
01028   TH1D* hTotalRunSpills = new TH1D("hTotalRunSpills", "hTotalRunSpills; Run; # Total Spills", 20000, 7000, 27000);
01029   TH1D* hPotPerTime = new TH1D("hPotPerTime", "hPotPerTime", 3652, StartTime, EndTime);
01030 
01031   TH1D* hTotalSnarlCount = new TH1D("TotalSnarlCount","TotalSnarlCount",2,0,2);
01032 
01033   TH1D* hTotalBwidx = new TH1D("hTotalBwidx", "hTotalBwidx; beam width x; # of spills", 1000,  -2e-3, 10e-3);
01034   TH1D* hTotalBwidy = new TH1D("hTotalBwidy", "hTotalBwidy; beam width y; # of spills", 1000,  -2e-3, 10e-3);
01035   TH1D* hTotalBposx = new TH1D("hTotalBposx", "hTotalBposx; beam pos x; # of spills", 5000, -10e-3, 10e-3);
01036   TH1D* hTotalBposy = new TH1D("hTotalBposy", "hTotalBposy; beam pos y; # of spills", 5000, -10e-3, 10e-3);
01037 
01038   TH2D* TimeVswidx = new TH2D("TimeVswidx","TimeVswidx; beam width x; time (s)", 500, -0.002, 0.002, 3652, StartTime, EndTime);
01039   TH2D* TimeVswidy = new TH2D("TimeVswidy","TimeVswidy; beam width y; time (s)", 500, -0.002, 0.002, 3652, StartTime, EndTime);
01040   TH2D* TimeVsposx = new TH2D("TimeVsposx","TimeVsposx; beam pos x; time (s)", 500, 0.0005, 0.0013, 3652, StartTime, EndTime);
01041   TH2D* TimeVsposy = new TH2D("TimeVsposy","TimeVsposy; beam pos y; time (s)", 500, 0.0014, 0.0022, 3652, StartTime, EndTime);
01042   
01043   NuEvent *event(0);
01044 
01045   // Total Snarl Count
01046   Int_t TotalSnarlCount = -999;
01047   
01048   // Create this for converting times
01049   TTimeStamp* s = new TTimeStamp();
01050 
01051   std::cout << "There are " << fileList->GetEntries() << " root files to loop over." << std::endl;
01052 
01053   for ( Int_t counter=0; counter<fileList->GetEntries(); ++counter ){
01054 
01055     TChainElement* chainElement = dynamic_cast<TChainElement*> (fileList->At(counter));
01056 
01057     if ( chainElement ){
01058       
01059       string fileName = chainElement->GetTitle();
01060       TFile f(fileName.c_str(),"READ");
01061       
01062       cout << "Filename = " << fileName.c_str() << endl;
01063       
01064       TTree *tree= (TTree*)f.Get("s");
01065       tree->SetBranchAddress("s",&event);
01066       tree->GetEntry(0);
01067 
01068       // I get these from the tfiles (sub run granularity)
01069       TH1D* hTotalPot = (TH1D*)f.Get("hTotalPot"); assert(hTotalPot);
01070       TH1D* hSnarlCount = (TH1D*)f.Get("hSnarlCount"); assert(hSnarlCount);
01071       TH1D* hBwidx    = (TH1D*)f.Get("hBwidx"); assert(hBwidx);
01072       TH1D* hBwidy    = (TH1D*)f.Get("hBwidy"); assert(hBwidy);
01073       TH1D* hBposx    = (TH1D*)f.Get("hBposx"); assert(hBposx);
01074       TH1D* hBposy    = (TH1D*)f.Get("hBposy"); assert(hBposy);
01075 
01076       std::cout << "The Pot for file " << fileName.c_str() << " is " << hTotalPot->Integral() << std::endl;
01077       hPotPerTime->Fill( event->timeSeconds, hTotalPot->Integral() );
01078 
01079       TotalSnarlCount += hSnarlCount->Integral();
01080 
01081       // Time to convert time, only so i can cout it for sanity checks
01082       s->SetSec(event->timeSeconds);
01083 
01084       cout << std::setprecision(10) << "the time = " << event->timeSeconds << endl;
01085 
01086       // Print time in sensible units.
01087       s->Print();
01088 
01089       hRunPot->Fill( event->run, hTotalPot->Integral() );
01090 
01091       // the integral of any of the hBid histograms is basically the number of good spills, so i use this
01092       hGoodRunSpills->Fill( event->run, hBwidx->GetEntries() );
01093       hTotalRunSpills->Fill( event->run, hSnarlCount->Integral() );
01094 
01095       hTotalBwidx->Add(hBwidx); 
01096       hTotalBwidy->Add(hBwidy);
01097       hTotalBposx->Add(hBposx);
01098       hTotalBposy->Add(hBposy);
01099 
01100       for(Int_t bin=1; bin<=hBwidx->GetNbinsX(); ++bin){
01101           TimeVswidx->Fill(hBwidx->GetXaxis()->GetBinCenter(bin), event->timeSeconds, hBwidx->GetBinContent(bin));
01102       }
01103       for(Int_t bin=1; bin<=hBwidy->GetNbinsX(); ++bin){
01104           TimeVswidy->Fill(hBwidy->GetXaxis()->GetBinCenter(bin), event->timeSeconds, hBwidy->GetBinContent(bin));
01105       }
01106       for(Int_t bin=1; bin<=hBposx->GetNbinsX(); ++bin){
01107           TimeVsposx->Fill(hBposx->GetXaxis()->GetBinCenter(bin), event->timeSeconds, hBposx->GetBinContent(bin));
01108       }
01109       for(Int_t bin=1; bin<=hBposy->GetNbinsX(); ++bin){
01110           TimeVsposy->Fill(hBposy->GetXaxis()->GetBinCenter(bin), event->timeSeconds, hBposy->GetBinContent(bin));
01111       }
01112 
01113       f.Close();
01114     }
01115   }
01116 
01117   std::cout << " " << std::endl;
01118   std::cout << "TotalSnarlCount for this set of files = " << TotalSnarlCount << std::endl;
01119   std::cout << " " << std::endl;
01120 
01121   hTotalSnarlCount->Fill(1, TotalSnarlCount);
01122 
01123   fOutFile->cd();
01124 
01125   hRunPot->Write();
01126   hTotalSnarlCount->Write();
01127   hTotalRunSpills->Write();
01128   hGoodRunSpills->Write();
01129   hPotPerTime->Write();
01130   hTotalBwidx->Write();
01131   hTotalBwidy->Write();
01132   hTotalBposx->Write();
01133   hTotalBposy->Write();
01134   TimeVswidx->Write();
01135   TimeVswidy->Write();
01136   TimeVsposx->Write();
01137   TimeVsposy->Write();
01138   
01139   std::cout << " " << std::endl;
01140 
01141 }
01142 
01143 
01144 
01145