EWCosmics_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       EWCosmics
// Module Type: analyzer
// File:        EWCosmics_module.cc
////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "canvas/Utilities/InputTag.h"
#include "canvas/Persistency/Common/FindMany.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "Geometry/Geometry.h"
#include "RecoBase/Track.h"
#include "RecoBase/FitSum.h"

#include <iostream>
#include <cassert>
#include <vector>
#include <fstream>
#include <string>
#include <map>
#include "TMath.h"
#include "TF1.h"
#include "TH2D.h"
#include "TGraph.h"
#include "TTree.h"
#include "TGeoManager.h"

#include "RecoBase/RecoHit.h"
#include "RecoBase/Cluster.h"
#include "MCCheater/BackTracker.h"
#include "Simulation/FLSHit.h"
#include "Simulation/FLSHitList.h"
#include "nusimdata/SimulationBase/MCTruth.h"

#include "RecoBase/FilterList.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Vertex.h" 
#include "RecoBase/Prong.h"
#include "Utilities/AssociationUtil.h"
#include "GeometryObjects/GeometryBase.h"

#include "CLHEP/Random/RandFlat.h"
#include "EWCosmicAnalyzer/OverburdenCalculator.h"

namespace red {
        class MatterProbe;
        class EWCosmics;
        class Energystimatter;
        class EnergystimatterConst;
        class EnergystimatterVar;
}

class red::Energystimatter{
        public:
                Energystimatter(){};
                virtual double EforX(double)=0;
};

class red::EnergystimatterConst : virtual public Energystimatter{
        public:
                EnergystimatterConst(){};
                void SetS(double);
                double EforX(double);
        private:
                double S;
};

void red::EnergystimatterConst::SetS(double alpha){
        S=alpha;
}

double red::EnergystimatterConst::EforX(double x){
        return S*x;
}

class red::EnergystimatterVar : virtual public Energystimatter{
        public:
                EnergystimatterVar(){};
        private:
};

class red::MatterProbe {
        public:
                MatterProbe(){};
                double mattersum(TVector3 startp, TVector3 startd);
                void SetGeoManager(TGeoManager* gm){fGeoManager=gm;}
                void SetMap();
                void ECalc(TVector3 startp, TVector3 startd);
                double GetDelE();
        private:
                TGeoManager* fGeoManager = nullptr;
                std::map<std::string,TGraph> tab,tabT;
                double deltaE;
};

double red::MatterProbe::mattersum(TVector3 startp, TVector3 startd){
        fGeoManager->InitTrack(startp.X(),startp.Y(),startp.Z(),-startd.X(),-startd.Y(),-startd.Z());
        double step,dens,sum=0;
        TGeoNode *nextnode;
        nextnode = fGeoManager->FindNextBoundaryAndStep();
        while(nextnode){
                step = fGeoManager->GetStep();
                if (step>1E10) break;
                
        //check that we are not in the detector volume anymore
        TString path=fGeoManager->GetPath();
        if(!path.Contains("/vDetEnclosure_")){
            TGeoVolume *vol;
                    vol = fGeoManager->GetCurrentVolume();
                    TGeoMaterial *mat;
                    mat = vol->GetMaterial();
                    dens = mat->GetDensity();
                    sum += step*dens;
                }
        nextnode = fGeoManager->FindNextBoundaryAndStep();
        }
        return sum;
}

void red::MatterProbe::SetMap(){
        const char* temp="%*lg %lg %*lg %*lg %*lg %*lg %*lg %*lg %lg %*lg %*lg";
        tab["Vacuum"]      =TGraph("EWCosmicAnalyzer/materials/muE_air_dry_1_atm.txt",temp);
        tab["Air"]         =TGraph("EWCosmicAnalyzer/materials/muE_air_dry_1_atm.txt",temp);
        tab["Water"]       =TGraph("EWCosmicAnalyzer/materials/muE_water_liquid.txt",temp);
        tab["ShotRock"]    =TGraph("EWCosmicAnalyzer/materials/muE_standard_rock.txt",temp);
        tab["Granite"]     =TGraph("EWCosmicAnalyzer/materials/muE_standard_rock.txt",temp);
        tab["Dirt"]        =TGraph("EWCosmicAnalyzer/materials/muE_standard_rock.txt",temp);
        tab["Concrete"]    =TGraph("EWCosmicAnalyzer/materials/muE_shielding_concrete.txt",temp);
        tab["BariteRock"]  =TGraph("EWCosmicAnalyzer/materials/muE_barium_sulfate.txt",temp);
        tab["Scintillator"]=TGraph("EWCosmicAnalyzer/materials/muE_paraffin.txt",temp);
        tab["PVC"]         =TGraph("EWCosmicAnalyzer/materials/muE_polyvinylchloride_PVC.txt",temp);
        tab["WLSFiber"]    =TGraph("EWCosmicAnalyzer/materials/muE_polyvinylchloride_PVC.txt",temp);
        tab["Glue"]        =TGraph("EWCosmicAnalyzer/materials/muE_carbon_amorphous_C.txt",temp);
        tab["Steel"]       =TGraph("EWCosmicAnalyzer/materials/muE_iron_Fe.txt",temp);
        tab["PivoterSteel"]=TGraph("EWCosmicAnalyzer/materials/muE_iron_Fe.txt",temp);
        for(auto it : tab){
                tabT[it.first]=TGraph(it.second.GetN(),it.second.GetY(),it.second.GetX());
        }

}

void red::MatterProbe::ECalc(TVector3 startp, TVector3 startd){
        TGraph *tg,*tgT;
        art::ServiceHandle<art::TFileService> tfs2;
        tg =tfs2->make<TGraph>();
        tgT=tfs2->make<TGraph>();

        fGeoManager->InitTrack(startp.X(),startp.Y(),startp.Z(),-startd.X(),-startd.Y(),-startd.Z());
        double step,dens,cX;
        TGeoNode *nextnode;
        nextnode = fGeoManager->FindNextBoundaryAndStep();
        double dE=0;
        while(nextnode){
                step = fGeoManager->GetStep();
                if (step>1E10) break;
        
        //check that we are not in the detector volume anymore
        TString path=fGeoManager->GetPath();
        if(!path.Contains("/vDetEnclosure_")){

            TGeoVolume *vol;
            vol = fGeoManager->GetCurrentVolume();
                TGeoMaterial *mat;
                    mat = vol->GetMaterial();
                    dens= mat->GetDensity();
                    cX  = step*dens;

                std::string matname=mat->GetName();
                auto it=tab.find(matname);
                if(it!=tab.end()){
                        *tg=TGraph(it->second);
                        cX+=tg->Eval(dE);
                }
                    auto itT=tabT.find(matname);
                    if(itT!=tabT.end()){
                            *tgT=TGraph(itT->second);
                            dE=tgT->Eval(cX);
                }
        }
                nextnode = fGeoManager->FindNextBoundaryAndStep();
        }
        deltaE=dE/1000.;
}

double red::MatterProbe::GetDelE(){
        return deltaE;
}

class red::EWCosmics : public art::EDAnalyzer {
        public:
                explicit EWCosmics(fhicl::ParameterSet const & p);
                EWCosmics(EWCosmics const &) = delete;
                EWCosmics(EWCosmics &&) = delete;
                EWCosmics & operator = (EWCosmics const &) = delete;
                EWCosmics & operator = (EWCosmics &&) = delete;

                void analyze(art::Event const & e) override;

                void beginJob() override;
                void endJob() override;

                double distToBorder(const TVector3& p);

                void FirstWay(TVector3, double&, double&);
                void SecondWay(const TVector3&, double&, double&);
        private:
                red::MatterProbe fMatterProber;
                red::EnergystimatterConst fEMC;

                std::string fSliceLabel;  
        std::string fTrackLabel;
        std::string fFitSumLabel;
        std::string fVertexLabel;
        std::string fClusterLabel;
        std::string fClusterInstance;
        std::string tracksLabel;
                double theta,phi,matter,Esurfv,Esurfc,th_d,ph_d,Ecal,EBPF,Erec;
                unsigned long Ne=0; // Event number
                TVector3 point,start,stop,startdir;
                TH2D *hEarth;
                TTree *tMu;

                const double phitoN=0;//26;
                const double FDmaxZ=5962.;//5959.08;//geo->DetLength();
                const double FDminZ=0.;
                const double FDmaxX=782.95;//geo->DetHalfWidth();
                const double FDminX=-FDmaxX;
                const double FDmaxY=774.7;//geo->DetHalfHeight();
                const double FDminY=-FDmaxY;
};

red::EWCosmics::EWCosmics(fhicl::ParameterSet const & p):
        EDAnalyzer(p)  // ,
        // More initializers here.
{
        std::cout<<"EWCosmicslyzer starting"<<std::endl;

    fSliceLabel=p.get<std::string>("SliceLabel");
    fTrackLabel=p.get<std::string>("TrackLabel");
    fFitSumLabel=p.get<std::string>("FitSumLabel");
    fVertexLabel=p.get<std::string>("VertexLabel");
    fClusterLabel=p.get<std::string>("ClusterLabel");
    fClusterInstance=p.get<std::string>("ClusterInstance");   
}

void red::EWCosmics::beginJob(){
        art::ServiceHandle<art::TFileService> tfs;
        //      hEarth=tfs->make<TH2D>("hEarth","Cosmic tracks; #phi;#theta",360,0,360,90,0,90);
        tMu=tfs->make<TTree>("tMu","Stopped cosmic mu tracks");
        //tMu->Branch("theta",&theta,"theta/D");
        //tMu->Branch("phi",&phi,"phi/D");
        tMu->Branch("th_d",&th_d,"th_d/D");
        tMu->Branch("ph_d",&ph_d,"ph_d/D");
        tMu->Branch("start",&start);
        tMu->Branch("stop",&stop);
        tMu->Branch("Ne",&Ne,"Ne/l");
        tMu->Branch("Esurfv",&Esurfv, "Esurfv/D");
        //tMu->Branch("Esurfc",&Esurfc, "Esurfc/D");
        //tMu->Branch("Ecal",&Ecal, "Ecal/D");
    tMu->Branch("EBPF",&EBPF, "EBPF/D");
//    tMu->Branch("Erec",&Erec, "Erec/D");
}

void red::EWCosmics::endJob(){ 
        std::cout<<"bye-bye!"<<std::endl;
}

double red::EWCosmics::distToBorder(const TVector3& p){
        double dx=std::max(std::min(p.X()-FDminX,FDmaxX-p.X()),0.);
        double dy=std::max(std::min(p.Y()-FDminY,FDmaxY-p.Y()),0.);
        double dz=std::max(std::min(p.Z()-FDminZ,FDmaxZ-p.Z()),0.);
        return std::min(dz,std::min(dx,dy));
}

void red::EWCosmics::FirstWay(TVector3 v, double& T, double& P){
        //Okay kid, this is where it gets complicated
        //We have Z pointing along the detector, X pointing left, and Y pointing up
        //We want X pointing along the detector, and Z pointing down
        //Let's rotate the vector
        v.RotateX(-TMath::PiOver2()); //now for vector new Y is old Z, Z is opposite old Y
        v.RotateZ(-TMath::PiOver2()); //now for vector the newest Y is opposite the first X. Phi is going clockwise
        double theta=v.Theta()*TMath::RadToDeg();
        //theta=180-theta; //because Z pointing down, but
        //we are interested in angles where CR comes FROM, not where TO, so
        //theta=180-theta; //again. Just ignore both operations
        double phi=v.Phi()*TMath::RadToDeg();
        //We are interested in angles where CR comes FROM, not where TO
        phi=-phi;
        //Phi is going clockwise. 0 is correspond particles coming along the detector from the North, 90 more or less from East
        phi-=phitoN; //now 0-particles are strict from the North, 90 from East, 270 from West
        //But azimuth should be counted from the South
        phi+=180; //now 180 is correspond particles coming from the North, 270 from East, 90(+360) from West. Oh, shut up
        while(phi<   0){phi+=360;}
        while(phi>=360){phi-=360;}
        T=theta;
        P=phi;
}

void red::EWCosmics::SecondWay(const TVector3& v, double& T, double& P){
        //Les just substitute axes in formulas
        double theta=TMath::ACos(v.Y()/TMath::Sqrt(v.X()*v.X()+v.Y()*v.Y()+v.Z()*v.Z()))*TMath::RadToDeg();
        double phi=TMath::ATan2(-v.X(),v.Z())*TMath::RadToDeg();
        //We are interested in angles where CR comes FROM, not where TO
        theta=180-theta;
        phi=-phi;
        //Phi is going clockwise. 0 is correspond particles coming along the detector from the North, 90 more or less from East
        phi-=phitoN; //now 0-particles are strict from the North, 90 from East, 270 from West
        //But azimuth should be counted from the South
        phi+=180; //now 180 is correspond particles coming from the North, 270 from East, 90(+360) from West. Oh, shut up
        while(phi<   0){phi+=360;}
        while(phi>=360){phi-=360;}
        T=theta;
        P=phi;
}

void red::EWCosmics::analyze(art::Event const & e){
    Ne = e.id().event();//event number

    art::ServiceHandle<geo::Geometry> geo;
        auto geomngr = geo->ROOTGeoManager();

        fMatterProber.SetGeoManager(geomngr);
        //fEMC.SetS(0.002);

    fMatterProber.SetMap();

    art::Handle<std::vector<rb::Cluster>> slice;
    e.getByLabel(fSliceLabel, slice);
    if((!slice.failedToGet())&&(!slice->empty())){
    art::FindManyP<rb::Track> trackh(slice, e, fTrackLabel);
    art::ServiceHandle<red::OverburdenCalculator> overcalc;

    for(size_t i = 0; i < trackh.size(); ++i) {
        // Find the breakpoint tracks associated to this slice
        std::vector< art::Ptr<rb::Track> > tracks = trackh.at(i);
        // Only look at slices with one track
        if (tracks.size()!=1)continue;

        const class rb::Track t=*tracks.at(0);

        start   =t.Start();
        stop    =t.Stop();
        startdir=t.Dir();
        
        if((red::EWCosmics::distToBorder(stop)<100)||(red::EWCosmics::distToBorder(start)>50)||(t.NXCell()<5)||(t.NYCell()<5)){continue;}
        red::EWCosmics::FirstWay(startdir,theta,phi);
        if(theta>90)continue;

        th_d=startdir.Theta()*TMath::RadToDeg();
        ph_d=startdir.Phi()*TMath::RadToDeg();
        while(ph_d<   0){ph_d+=360;}
        while(ph_d>=360){ph_d-=360;}
        
        fMatterProber.ECalc(stop,startdir);
        Esurfv=fMatterProber.GetDelE();
        matter=overcalc->traceBack(t);
        std::cout<<"Matter = "<<matter<<std::endl;
        //Esurfc=fEMC.EforX(matter);

        art::FindManyP<rb::FitSum> kininf(tracks, e, fFitSumLabel);
        std::vector<art::Ptr<rb::FitSum>> kinin = kininf.at(0);
        art::Ptr<rb::FitSum> kin = kinin.at(0); //track kinetic information
        TLorentzVector mom = kin->FourMom();
        EBPF = mom.E();
        
//        Erec=Esurfv+EBPF;
        tMu->Fill();}
    }
}

DEFINE_ART_MODULE(red::EWCosmics)