Analysis.cxx

/*****************************************************************************
 * Copyright (C) 2009-2014    this file is part of the NPTool Project        *
 *                                                                           *
 * For the licensing terms see $NPTOOL/Licence/NPTool_Licence                *
 * For the list of contributors see $NPTOOL/Licence/Contributors             *
 *****************************************************************************/

/*****************************************************************************
 * Original Author: Adrien MATTA  contact address: a.matta@surrey.ac.uk      *
 *                                                                           *
 * Creation Date  : july  2020                                               *
 * Last update    :                                                          *
 *---------------------------------------------------------------------------*
 * Decription:                                                               *
 * Class describing the property of an Analysis object                       *
 *                                                                           *
 *---------------------------------------------------------------------------*
 * Comment:                                                                  *
 *                                                                           *
 *                                                                           *
 *****************************************************************************/
#include<iostream>
using namespace std;
#include"Analysis.h"
#include"NPAnalysisFactory.h"
#include"NPDetectorManager.h"
#include"NPOptionManager.h"
#include"NPFunction.h"
#include"NPTrackingUtility.h"
#include"NPPhysicalConstants.h"
////////////////////////////////////////////////////////////////////////////////
Analysis::Analysis(){
}
////////////////////////////////////////////////////////////////////////////////
Analysis::~Analysis(){
}

////////////////////////////////////////////////////////////////////////////////
void Analysis::Init(){
  IC= new TInitialConditions;
  DC= new TInteractionCoordinates;
  RC= new TReactionConditions;
  
  InitOutputBranch();
  InitInputBranch();
  
  Strasse = (TStrassePhysics*)  m_DetectorManager -> GetDetector("Strasse");
  Catana = (TCatanaPhysics*)  m_DetectorManager -> GetDetector("Catana");
  // reaction properties
  m_QFS = new NPL::QFS();
  m_QFS->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
  // reaction properties
  myBeam = new NPL::Beam();
  myBeam->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
  InitialBeamEnergy = myBeam->GetEnergy()/* * myBeam->GetA()*/;
  // target thickness
  TargetThickness = m_DetectorManager->GetTargetThickness();
  string TargetMaterial = m_DetectorManager->GetTargetMaterial();
  // EnergyLoss Tables
  string BeamName = NPL::ChangeNameToG4Standard(myBeam->GetName());
  BeamTarget = NPL::EnergyLoss(BeamName+"_"+TargetMaterial+".G4table","G4Table",100);
  protonTarget = NPL::EnergyLoss("proton_"+TargetMaterial+".G4table","G4Table",100);
  protonAl = NPL::EnergyLoss("proton_Al.G4table","G4Table",100);
  protonSi = NPL::EnergyLoss("proton_Si.G4table","G4Table",100);
  LV_T.SetVectM(TVector3(0,0,0),NPUNITS::proton_mass_c2);
} 

////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent(){
  // Reinitiate calculated variable
  ReInitValue();
  unsigned int size = Strasse->GetEventMultiplicity();
  if(size==2){ // 2 proton detected
    // Proton 1
    TVector3 InnerPos1 = Strasse->GetInnerPositionOfInteraction(0);
    TVector3 OuterPos1 = Strasse->GetOuterPositionOfInteraction(0);
    TVector3 Proton1 = OuterPos1-InnerPos1;
    Proton1=Proton1.Unit();
    // Proton 2
    TVector3 InnerPos2 = Strasse->GetInnerPositionOfInteraction(1);
    TVector3 OuterPos2 = Strasse->GetOuterPositionOfInteraction(1);
    TVector3 Proton2 = OuterPos2-InnerPos2;
    Proton2=Proton2.Unit();

    deltaPhi = abs(Proton1.Phi()/deg-Proton2.Phi()/deg);
    sumTheta = Proton1.Theta()/deg+Proton2.Theta()/deg;
    Theta12  = Proton1.Angle(Proton2)/deg;

    // reject event that make no physical sense
    /*if(deltaPhi<170 && sumTheta<80){
      return;
      }
    */
    // computing minimum distance of the two lines
    TVector3 Vertex;
    TVector3 delta;
    Distance = NPL::MinimumDistanceTwoLines(InnerPos1,OuterPos1,InnerPos2,OuterPos2,Vertex,delta);
    VertexX=Vertex.X();
    VertexY=Vertex.Y();
    VertexZ=Vertex.Z();
    deltaX=delta.X();
    deltaY=delta.Y();
    deltaZ=delta.Z();
    
    // Look for associated Catana event
    double d1,d2;
    unsigned int i1,i2;
    i1 = Catana->FindClosestHitToLine(InnerPos1,OuterPos1,d1);
    i2 = Catana->FindClosestHitToLine(InnerPos2,OuterPos2,d2);
    //if(i1!=i2){
    if(true){
      double TA = BeamTarget.Slow(InitialBeamEnergy,abs(VertexZ-75),RC->GetBeamDirection().Angle(TVector3(0,0,1)));
       
      ////////////////////////////////////
      // From Reaction Conditions
      double E1s = RC->GetKineticEnergy(0); 
      double E2s = RC->GetKineticEnergy(1); 
      TVector3 Proton1s=RC->GetParticleMomentum(0).Unit();
      TVector3 Proton2s=RC->GetParticleMomentum(1).Unit();
      // Matching the right energy with the right proton
      if((Proton1s-Proton1).Mag()<(Proton1s-Proton2).Mag()){
        E1=E1s;
        E2=E2s;
        alpha=Proton1s.Angle(Proton1)/deg;
        Theta1=Proton1.Theta();
        Phi1=Proton1.Phi();
        Theta2=Proton2.Theta();
        Phi2=Proton2.Phi();
        Theta1s=Proton1s.Theta();
        Phi1s=Proton1s.Phi();
        Theta2s=Proton2s.Theta();
        Phi2s=Proton2s.Phi();
        
        }
      else{
        E2=E1s;
        E1=E2s;
        alpha=Proton2s.Angle(Proton1)/deg;
        Theta1=Proton1.Theta()/deg;
        Phi1=Proton1.Phi()/deg;
        Theta2=Proton2.Theta()/deg;
        Phi2=Proton2.Phi()/deg;
        Theta1s=Proton2s.Theta()/deg;
        Phi1s=Proton2s.Phi()/deg;
        Theta2s=Proton1s.Theta()/deg;
        Phi2s=Proton1s.Phi()/deg;
        }
      // From detectors
      E1 = ReconstructProtonEnergy(Vertex,Proton1,Catana->Energy[i1]); 
      E2 = ReconstructProtonEnergy(Vertex,Proton2,Catana->Energy[i2]);

     // Vertex = RC->GetVertexPosition();  
      //TA = RC->GetBeamEnergy();
      //////////////////////////////////// 
      
     // setting up Lorentz Vector from measured trajectories and energies
     // TVector3 PA(0,0,sqrt(TA*(TA+2*m_QFS->GetParticleA()->Mass()))); // for like there is no BDC
     double beam_mom=sqrt(TA*(TA+2*m_QFS->GetParticleA()->Mass()));
    // TVector3 PA(0,0,beam_mom); // for like there is no BDC
       TVector3 PA=beam_mom*RC->GetBeamDirection().Unit();

      LV_A.SetVectM(PA,m_QFS->GetParticleA()->Mass());
      double P1= sqrt(E1*(E1+2*NPUNITS::proton_mass_c2));
      double P2= sqrt(E2*(E2+2*NPUNITS::proton_mass_c2));

      LV_p1.SetVectM(Proton1.Unit()*P1,NPUNITS::proton_mass_c2); 
      LV_p2.SetVectM(Proton2.Unit()*P2,NPUNITS::proton_mass_c2); 
      // computing Ex from Missing Mass
      LV_B = LV_A + LV_T - LV_p1 - LV_p2;
      //LV_B = RC->GetParticleMomentum(2);
      Ex = LV_B.M() - m_QFS->GetParticleB()->Mass();
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
double Analysis::ReconstructProtonEnergy(const TVector3& x0, const TVector3& dir,const double& Ecatana){
    TVector3 Normal = TVector3(0,1,0);
    Normal.SetPhi(dir.Phi());
    double Theta = dir.Angle(Normal);  
    // Catana Al housing 
    double E = protonAl.EvaluateInitialEnergy(Ecatana,0.5*mm,Theta);
    // Strasse Chamber
    E = protonAl.EvaluateInitialEnergy(E,3*mm,Theta);
    // Outer Barrel
    E = protonSi.EvaluateInitialEnergy(E,300*micrometer,Theta);
    // Inner Barrel
    E = protonSi.EvaluateInitialEnergy(E,200*micrometer,Theta);
    // LH2 target
    static TVector3 x1;
    x1= x0+dir;
    TVector3 T1(0,15,0);
    TVector3 T2(0,15,1);
    T1.SetPhi(dir.Phi());
    T2.SetPhi(dir.Phi());
    double d = NPL::MinimumDistancePointLine(T1,T2,x0);
    E = protonTarget.EvaluateInitialEnergy(E,d,Theta);
  }


////////////////////////////////////////////////////////////////////////////////
TVector3 Analysis::InterpolateInPlaneZ(TVector3 V0, TVector3 V1, double Zproj){
    TVector3 Vproj(-999,-999,-999);
    double t = (Zproj - V1.Z()) / (V1.Z()-V0.Z());
    double Xproj= V1.X() + (V1.X()-V0.X()) * t;
    double Yproj= V1.Y() + (V1.Y()-V0.Y()) * t; 
    Vproj.SetXYZ(Xproj,Yproj,Zproj);
    return Vproj;
}

////////////////////////////////////////////////////////////////////////////////
void Analysis::End(){
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch() {
  RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
  RootOutput::getInstance()->GetTree()->Branch("E1",&E1,"E1/D");
  RootOutput::getInstance()->GetTree()->Branch("E2",&E2,"E2/D");
  RootOutput::getInstance()->GetTree()->Branch("Theta12",&Theta12,"Theta12/D");
  RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
  RootOutput::getInstance()->GetTree()->Branch("VertexX",&VertexX,"VertexX/D");
  RootOutput::getInstance()->GetTree()->Branch("VertexY",&VertexY,"VertexY/D");
  RootOutput::getInstance()->GetTree()->Branch("VertexZ",&VertexZ,"VertexZ/D");
  RootOutput::getInstance()->GetTree()->Branch("deltaX",&deltaX,"deltaX/D");
  RootOutput::getInstance()->GetTree()->Branch("deltaY",&deltaY,"deltaY/D");
  RootOutput::getInstance()->GetTree()->Branch("deltaZ",&deltaZ,"deltaZ/D");
  RootOutput::getInstance()->GetTree()->Branch("deltaPhi",&deltaPhi,"deltaPhi/D");
  RootOutput::getInstance()->GetTree()->Branch("sumTheta",&sumTheta,"sumTheta/D");

  RootOutput::getInstance()->GetTree()->Branch("alpha",&alpha,"alpha/D");

  RootOutput::getInstance()->GetTree()->Branch("Theta1",&Theta1,"Theta1/D");
  RootOutput::getInstance()->GetTree()->Branch("Phi1",&Phi1,"Phi1/D");
  RootOutput::getInstance()->GetTree()->Branch("Theta2",&Theta2,"Theta2/D");
  RootOutput::getInstance()->GetTree()->Branch("Phi2",&Phi2,"Phi2/D");
  RootOutput::getInstance()->GetTree()->Branch("Theta1s",&Theta1s,"Theta1s/D");
  RootOutput::getInstance()->GetTree()->Branch("Phi1s",&Phi1s,"Phi1s/D");
  RootOutput::getInstance()->GetTree()->Branch("Theta2s",&Theta2s,"Theta2s/D");
  RootOutput::getInstance()->GetTree()->Branch("Phi2s",&Phi2s,"Phi2s/D");


  RootOutput::getInstance()->GetTree()->Branch("Distance",&Distance,"Distance/D");
  RootOutput::getInstance()->GetTree()->Branch("InteractionCoordinates","TInteractionCoordinates",&DC);
  RootOutput::getInstance()->GetTree()->Branch("ReactionConditions","TReactionConditions",&RC);
}

////////////////////////////////////////////////////////////////////////////////
void Analysis::InitInputBranch(){
    RootInput:: getInstance()->GetChain()->SetBranchAddress("InteractionCoordinates",&DC);
    RootInput:: getInstance()->GetChain()->SetBranchAddress("ReactionConditions",&RC);
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInitValue(){
  Ex = -1000 ;
  E1= -1000;
  E2 = -1000;
  Theta12 = -1000;
  ThetaCM = -1000;
  VertexX=-1000;
  VertexY=-1000;
  VertexZ=-1000;
  deltaX=-1000;
  deltaY=-1000;
  deltaZ=-1000;
  Distance=-1000;
  sumTheta=-1000;
  deltaPhi=-1000;
  alpha=-1000;
  Theta1=-1000;
  Phi1=-1000;
  Theta2=-1000;
  Phi2=-1000;
  Theta1s=-1000;
  Phi1s=-1000;
  Theta2s=-1000;
  Phi2s=-1000;
}


////////////////////////////////////////////////////////////////////////////////
//            Construct Method to be pass to the AnalysisFactory              //
////////////////////////////////////////////////////////////////////////////////
NPL::VAnalysis* Analysis::Construct(){
  return (NPL::VAnalysis*) new Analysis();
}

////////////////////////////////////////////////////////////////////////////////
//            Registering the construct method to the factory                 //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
class proxy{
  public:
    proxy(){
      NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
    }
};

proxy p;
}