LBNEQuickPiToNu.hh

Go to the documentation of this file.

//
// LBNEQuickPitoNu.hh
//

#ifndef LBNEQuickPitoNu_H
#define LBNEQuickPitoNu_H 1
#include <fstream>
#include <iostream>
#include "globals.hh"
#include "G4Event.hh"
#include "G4EventManager.hh"
#include "G4RunManager.hh"
#include "G4ExceptionSeverity.hh"

class G4EventManager;
class LBNEEventAction;
class LBNEQuickPitoNu;
class LBNERunManager;
//
// These actions are for either physics studies, or debugging. They are expected to run either in " Proton on target mode", 
// or on neutrinos.  This a singleton class. Means that an implicit global pointer, tough 
// Hopefull nobody will complain. Takes very little memory if not filled. 
//


class  LBNEQuickPiToNu
{
// see below for use and reasons.. 
  public:
   LBNEQuickPiToNu();
   int evtNum;
   int trNum;
   int sign;
   double piEnergy;
   double nuEnergy;
   std::vector<double> piMomentumTarget;
   std::vector<double> piPosAtHangerRing;
   std::vector<double> piMomentumAtHangerRing;
   std::vector<double> nuMomentum;
};

//
// Investigating why fixing the Ring Hanger diameter gave no change at 3GeV. 
// We will try to correlate the pion that crosses the Horn 2 Hanger Ring (Wrong diameter!)
// With the direction of the corresponding neutrino. 
// This version is probably just a protoype: Investigating how much flux is lost by an 
// intervening volume is likely to be an on-going problem.. 
//
class  LBNEQuickPiToNuVect {
  private:
    LBNEQuickPiToNuVect();
    LBNEQuickPiToNuVect(LBNEQuickPiToNuVect const &); // not implemented 
    LBNEQuickPiToNuVect& operator=(LBNEQuickPiToNuVect const&); // not implemented 
    static LBNEQuickPiToNuVect* m_pInstance;
    bool doItForReal;
    std::vector<LBNEQuickPiToNu> data;
    std::ofstream fOut; 
   
  public:
     static LBNEQuickPiToNuVect* Instance();
     inline void YesDoThis() { doItForReal = true;}
     inline bool doIt() const {return doItForReal;}
     inline void open(const std::string &prefix) {
       if (fOut.is_open()) {
          std::cerr << " From LBNEQuickPiToNuVect, file has been already opened, close it " << std::endl;
          fOut.close();
       }
       std::string fName(prefix); fName += std::string("QuickPitoNu_v1.txt");
       fOut.open(fName.c_str());
       fOut << " evtNum trNum sign piEnergy nuEnergy piPxTgt piPyTgt pixH2 piyH2 pizH2 piPxH2 piPyH2 ";
       fOut << " nuPx nuPy " << std::endl; // P stands for slope ,Px/Pz... 
     }
     inline void close() {fOut.close();} 
     inline void addPion(int trNum1, int sign1, double e, G4ThreeVector p) {
       LBNERunManager* pRunManager=(LBNERunManager*)LBNERunManager::GetRunManager();
       int evtno = pRunManager->GetCurrentEvent()->GetEventID();       
       LBNEQuickPiToNu a; a.evtNum = evtno; a.sign = sign1;
       a.trNum = trNum1; a.piEnergy=e;
        for (size_t k=0; k!=3; k++) a.piMomentumTarget[k]= p[k];
        data.push_back(a);
     }
     inline void RingHangerPion(int trNum1, int sign, G4ThreeVector p,  G4ThreeVector xPos) {
       bool foundPion = false;
       for (std::vector<LBNEQuickPiToNu>::iterator it=data.begin(); it!=data.end(); it++) {
         if ((it->trNum == -9999) || (it->trNum == -9998)) continue;
         if (it->trNum != trNum1) continue;
         for (size_t k=0; k!= 3; k++) {
           it->piMomentumAtHangerRing[k] = p[k];
           it->piPosAtHangerRing[k] = xPos[k];
         }
         foundPion = true; return;
       }
       if (!foundPion) { 
       const double e = std::sqrt((139.58*139.58) + p[0]*p[0] +  p[1]*p[1] + p[2]*p[2]);
        std::cerr << " Found new secondary pion, declare it with a higher track number " << std::endl;
        addPion(10000+trNum1, sign, e, p);
       }
     }
     
     inline void blessPion(int trNum1, double e, G4ThreeVector p) {
      bool foundPion = false;
      for (std::vector<LBNEQuickPiToNu>::iterator it=data.begin(); it!=data.end(); it++) {
        if ((it->trNum == -9999) || (it->trNum == -9998)) continue;
        if ((it->trNum != trNum1) && ((it->trNum-10000) != trNum1)) continue;
        for (size_t k=0; k!= 3; k++) {
           it->nuMomentum[k] = p[k];
         }
         it->nuEnergy = e;
         foundPion = true; return;
        if (!foundPion) { 
            std::ostringstream mStrStr;
            mStrStr << "Found Unknown pion.. Should not happen.. Stop here. Fatal " << std::endl;
            G4String mStr(mStrStr.str());
            G4Exception("LBNERunAction::BeginOfRunAction", " ", RunMustBeAborted, mStr.c_str()); 
        }
      }
    }
    inline void reset () { data.clear();}
    
    inline void evtOut() {
      for (std::vector<LBNEQuickPiToNu>::const_iterator it=data.begin(); it!=data.end(); it++) {
        if ((it->trNum == -9999) || (it->trNum == -9998)) continue;
        if (std::abs(it->nuEnergy) < 0.001) continue; // gave no neutrino..
        fOut << " " << it->evtNum << " " << it->trNum << " " << it->sign << " " << it->piEnergy;
        fOut << " " << it->nuEnergy;
        fOut << " " << it->piMomentumTarget[0]/it->piMomentumTarget[2] << " " 
                    << it->piMomentumTarget[1]/it->piMomentumTarget[2];
        for (size_t k=0; k != 3; k++) fOut << " " << it->piPosAtHangerRing[k];
        const double rSq = it->piPosAtHangerRing[0]*it->piPosAtHangerRing[0] + 
                           it->piPosAtHangerRing[1]*it->piPosAtHangerRing[1];
        if (rSq > 1.0e-6)                  
           fOut << " " << it->piMomentumAtHangerRing[0]/it->piMomentumAtHangerRing[2] << " "
                       << it->piMomentumAtHangerRing[1]/it->piMomentumAtHangerRing[2];
        else fOut << " 0. 0. ";        
        fOut << " " << it->nuMomentum[0]/it->nuMomentum[2] << " " << it->nuMomentum[1]/it->nuMomentum[2];
        fOut << std::endl;       
      }
      this->reset();
    }
};
#endif