Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 10:37:38 

0001 // Λc⁺ → p K⁻ π⁺ reconstruction with (6.28 ± 0.32)%
0002 // Shyam Kumar; INFN Bari, Italy
0003 // shyam.kumar@ba.infn.it; shyam055119@gmail.com 
0004 
0005 #ifdef __CINT__
0006 
0007 #pragma link off all globals;
0008 #pragma link off all classes;
0009 #pragma link off all functions;
0010 
0011 #pragma link C++ class PlotFile;
0012 #endif
0013 
0014 #ifndef __CINT__
0015 #include <stdio.h>
0016 #include <stdlib.h>
0017 #include <fstream>
0018 #include <iostream>
0019 #include <iomanip>
0020 #include <string>
0021 #include <sys/types.h>
0022 #include <sys/stat.h>
0023 #include <dirent.h>
0024 #include "math.h"
0025 #include "string.h"
0026 
0027 #include "TROOT.h"
0028 #include "TFile.h"
0029 #include "TChain.h"
0030 #include "TH1D.h"
0031 #include "TH2D.h"
0032 #include "TH3D.h"
0033 #include "THnSparse.h"
0034 #include "TStyle.h"
0035 #include "TCanvas.h"
0036 #include "TProfile.h"
0037 #include "TTree.h"
0038 #include "TNtuple.h"
0039 #include "TRandom3.h"
0040 #include "TMath.h"
0041 #include "TSystem.h"
0042 #include "TUnixSystem.h"
0043 #include "TVector2.h"
0044 #include "TVector3.h"
0045 #include "TLorentzVector.h"
0046 #include "TTreeReader.h"
0047 #include "TTreeReaderValue.h"
0048 #include "TTreeReaderArray.h"
0049 #include "TLatex.h"
0050 #include "TMinuit.h"
0051 #include "Math/Functor.h"
0052 #include "Fit/Fitter.h"
0053 #include "Math/Minimizer.h"
0054 #endif
0055 #include "StPhysicalHelix.h"
0056 #include "SystemOfUnits.h"
0057 #include "PhysicalConstants.h"
0058 
0059 using namespace std;
0060 
0061 const double gPionMass = 0.13957;
0062 const double gKaonMass = 0.493677;
0063 const double gProtonMass = 0.938272;
0064 
0065 const double twoPi = 2.*3.1415927;
0066 const double eMass = 0.000511;
0067 const int pdg_Lcp = 4122, pdg_p = 2212, pdg_k = 321, pdg_pi = 211;
0068 bool debug = false;
0069 
0070 
0071 const double bField = -1.7; // Tesla
0072 
0073 void getDecayVertex_chi2fit(const int index1, const int index2, const int index3, double &s1, double &s2, double &s3, TVector3 &vertex, float &chi2, double *parFitErr);
0074 TVector3 GetDCAToPrimaryVertex(const int index, TVector3 vtx);
0075 
0076 TLorentzVector GetDCADaughters(const int index1, const int index2, const int index3, TVector3 vtx,
0077   float *dcaDaughters, float &cosTheta, float &cosTheta_xy, float &decayLength, float &V0DcaToVtx, float &sigma_vtx, float & chi2_ndf, TVector3 &decayVertex, double *parFitErr);
0078 
0079 TTreeReaderArray<float> *rcMomPx2;
0080 TTreeReaderArray<float> *rcMomPy2;
0081 TTreeReaderArray<float> *rcMomPz2;
0082 TTreeReaderArray<float> *rcCharge2;
0083 
0084 TTreeReaderArray<float> *rcTrkLoca2;
0085 TTreeReaderArray<float> *rcTrkLocb2;
0086 TTreeReaderArray<float> *rcTrkTheta2;
0087 TTreeReaderArray<float> *rcTrkPhi2;
0088 TTreeReaderArray<std::array<float, 21>> *rcTrkCov;
0089 
0090 // Define function for Chi2 minimization between two helices    
0091 struct Chi2Minimization {    
0092      StPhysicalHelix fhelix1, fhelix2, fhelix3;
0093      std::array<float, 21> fcov1, fcov2, fcov3; // full covariance matrix
0094      
0095      Chi2Minimization(StPhysicalHelix helix1, StPhysicalHelix helix2, StPhysicalHelix helix3, std::array<float, 21> cov1, std::array<float, 21> cov2, std::array<float, 21> cov3) : fhelix1(helix1),fhelix2(helix2), fhelix3(helix3), fcov1(cov1), fcov2(cov2), fcov3(cov3) {}
0096     // Implementation of the function to be minimized
0097     double operator() (const double *par) {
0098     double x = par[0];
0099     double y = par[1];
0100     double z = par[2];
0101     double s1 = par[3];
0102     double s2 = par[4];
0103     double s3 = par[5];    
0104     double f = 0;
0105     TVector3 vertex(x, y, z);
0106     TVector3 p1 = fhelix1.at(s1);
0107     TVector3 p2 = fhelix2.at(s2);
0108     TVector3 p3 = fhelix3.at(s3);
0109     TVector3 mom1 = fhelix1.momentumAt(s1,  bField * tesla);
0110     TVector3 mom2 = fhelix2.momentumAt(s2,  bField * tesla);
0111     TVector3 mom3 = fhelix3.momentumAt(s3,  bField * tesla);
0112      // x= −l0 sinϕ , y=l0 cosϕ , z=l
0113    // Recalculate l0 at PCA for error propagation
0114   float l0_track1 = p1.Pt(); float l1_track1 = p1.Z(); double phi_track1 = mom1.Phi();
0115   float l0_track2 = p2.Pt(); float l1_track2 = p2.Z(); double phi_track2 = mom2.Phi();
0116   float l0_track3 = p3.Pt(); float l1_track3 = p3.Z(); double phi_track3 = mom3.Phi();  
0117      // Track1: σx^2​=sin^2ϕ⋅σℓ0​^2​+ℓ0^2​cos^2ϕ⋅σϕ^2​+2⋅ℓ0​sinϕcosϕ⋅Cov(ℓ0​,ϕ)
0118   float sigx1_2 = sin(phi_track1)*sin(phi_track1)*fcov1[0] + l0_track1*l0_track1*cos(phi_track1)*cos(phi_track1)*fcov1[5]+ 2.0*l0_track1*sin(phi_track1)*cos(phi_track1)*fcov1[3]; 
0119   float sigx2_2 = sin(phi_track2)*sin(phi_track2)*fcov2[0] + l0_track2*l0_track2*cos(phi_track2)*cos(phi_track2)*fcov2[5]+ 2.0*l0_track2*sin(phi_track2)*cos(phi_track2)*fcov2[3];
0120   float sigx3_2 = sin(phi_track3)*sin(phi_track3)*fcov3[0] + l0_track3*l0_track3*cos(phi_track3)*cos(phi_track3)*fcov3[5]+ 2.0*l0_track3*sin(phi_track3)*cos(phi_track3)*fcov3[3];
0121   
0122      // σy^2​=cos^2ϕ⋅σℓ0^​2​+ℓ0^2​sin^2ϕ⋅σϕ^2​−2⋅ℓ0​sinϕcosϕ⋅Cov(ℓ0​,ϕ)
0123   float sigy1_2 = cos(phi_track1)*cos(phi_track1)*fcov1[0] + l0_track1*l0_track1*sin(phi_track1)*sin(phi_track1)*fcov1[5]-2.0*l0_track1*sin(phi_track1)*cos(phi_track1)*fcov1[3]; 
0124   float sigy2_2 = cos(phi_track2)*cos(phi_track2)*fcov2[0] + l0_track2*l0_track2*sin(phi_track2)*sin(phi_track2)*fcov2[5]-2.0*l0_track2*sin(phi_track2)*cos(phi_track2)*fcov2[3]; 
0125   float sigy3_2 = cos(phi_track3)*cos(phi_track3)*fcov3[0] + l0_track3*l0_track3*sin(phi_track3)*sin(phi_track3)*fcov3[5]-2.0*l0_track3*sin(phi_track3)*cos(phi_track3)*fcov3[3]; 
0126      
0127      // σz^2​
0128   float sigz1_2 = fcov1[2];
0129   float sigz2_2 = fcov2[2]; 
0130   float sigz3_2 = fcov3[2];   
0131   double d1_x = 10.*(vertex - p1).X(); double d2_x = 10.*(vertex - p2).X(); double d3_x = 10.*(vertex - p3).X();
0132   double d1_y = 10.*(vertex - p1).Y(); double d2_y = 10.*(vertex - p2).Y(); double d3_y = 10.*(vertex - p3).Y();
0133   double d1_z = 10.*(vertex - p1).Z(); double d2_z = 10.*(vertex - p2).Z(); double d3_z = 10.*(vertex - p3).Z();
0134     
0135   f = d1_x*d1_x/sigx1_2 + d2_x*d2_x/sigx2_2 + d3_x*d3_x/sigx3_2 + d1_y*d1_y/sigy1_2 + d2_y*d2_y/sigy2_2 + d3_y*d3_y/sigy3_2 + d1_z*d1_z/sigz1_2+ d2_z*d2_z/sigz2_2 + d3_z*d3_z/sigz3_2; // chi2 
0136 
0137       return f;
0138     }
0139     };
0140 
0141 
0142 // Define function for Distance2 minimization between three helices 
0143 struct Distance2Minimization {    
0144      StPhysicalHelix fhelix1, fhelix2, fhelix3;
0145      Distance2Minimization(StPhysicalHelix helix1, StPhysicalHelix helix2, StPhysicalHelix helix3) : fhelix1(helix1),fhelix2(helix2), fhelix3(helix3) {}
0146     // Implementation of the function to be minimized
0147     double operator() (const double *par) {
0148     double x = par[0];
0149     double y = par[1];
0150     double z = par[2];
0151     double s1 = par[3];
0152     double s2 = par[4]; 
0153     double s3 = par[5];   
0154     double f = 0;
0155     
0156     TVector3 vertex(x, y, z);
0157     TVector3 p1 = fhelix1.at(s1);
0158     TVector3 p2 = fhelix2.at(s2);
0159     TVector3 p3 = fhelix3.at(s3);
0160 
0161     double d1 = (vertex - p1).Mag2();
0162     double d2 = (vertex - p2).Mag2();
0163     double d3 = (vertex - p3).Mag2();
0164     f = d1 + d2 + d3; // minimize total squared distance
0165     return f;   
0166     }
0167     };
0168 
0169 int main(int argc, char **argv)
0170 {
0171   if(argc!=3 && argc!=1) return 0;
0172 
0173   TString listname;
0174   TString outname;
0175 
0176   if(argc==1)
0177   {
0178     listname  = "test.list";
0179     outname = "test.root";
0180   }
0181 
0182   if(argc==3)
0183   {
0184     listname = argv[1];
0185     outname = argv[2];
0186   }  
0187 
0188   TChain *chain = new TChain("events");
0189 
0190   int nfiles = 0;
0191   char filename[512];
0192   ifstream *inputstream = new ifstream;
0193   inputstream->open(listname.Data());
0194   if(!inputstream)
0195   {
0196     printf("[e] Cannot open file list: %s\n", listname.Data());
0197   }
0198   while(inputstream->good())
0199   {
0200     inputstream->getline(filename, 512);
0201     if(inputstream->good())
0202     {
0203      TFile *ftmp = TFile::Open(filename, "read");
0204      if(!ftmp||!(ftmp->IsOpen())||!(ftmp->GetNkeys())) 
0205      {
0206        printf("[e] Could you open file: %s\n", filename);
0207      } 
0208      else
0209      {
0210        cout<<"[i] Add "<<nfiles<<"th file: "<<filename<<endl;
0211        chain->Add(filename);
0212        nfiles++;
0213      }
0214    }
0215  }
0216  inputstream->close();
0217  printf("[i] Read in %d files with %lld events in total\n", nfiles, chain->GetEntries());
0218 
0219  TH1F *hEventStat = new TH1F("hEventStat", "Event statistics", 7, 0, 7);
0220  hEventStat->GetXaxis()->SetBinLabel(1, "MC events");
0221  hEventStat->GetXaxis()->SetBinLabel(2, "#Lambda_{c}^{+}");
0222  hEventStat->GetXaxis()->SetBinLabel(3, "#Lambda_{c}^{+} -> p + K^{-} + #pi^{+}");
0223  hEventStat->GetXaxis()->SetBinLabel(4, "Reco Signal #Lambda_{c}^{+}/#Lambda_{c}^{-}");
0224  hEventStat->GetXaxis()->SetBinLabel(5, "Reco Signal #Lambda_{c}^{+}");
0225  hEventStat->GetXaxis()->SetBinLabel(6, "Reco Signal #Lambda_{c}^{-}");
0226  hEventStat->GetXaxis()->SetBinLabel(7, "Reco Bkg #Lambda_{c}^{+}");
0227 
0228 
0229  TH1F *hMcMult = new TH1F("hMcMult", "MC multiplicity (|#eta| < 3.5);N_{MC}", 50, 0, 50);
0230 
0231  TH1F *hMcVtxX = new TH1F("hMcVtxX", "x position of MC vertex;x (mm)", 500, -5.0, 5.0);
0232  TH1F *hMcVtxY = new TH1F("hMcVtxY", "y position of MC vertex;y (mm)", 500, -5.0, 5.0);
0233  TH1F *hMcVtxZ = new TH1F("hMcVtxZ", "z position of MC vertex;z (mm)", 800, -200, 200);
0234 
0235  TH1F *hRecVtxX = new TH1F("hRecVtxX", "x position of Rec vertex;x (mm)", 500, -5.0, 5.0);
0236  TH1F *hRecVtxY = new TH1F("hRecVtxY", "y position of Rec vertex;y (mm)", 500, -5.0, 5.0);
0237  TH1F *hRecVtxZ = new TH1F("hRecVtxZ", "z position of Rec vertex;z (mm)", 800, -200, 200);
0238 
0239  TH1F *hPullVtxX = new TH1F("hPullVtxX", "Pull x position of MC vertex;(Vx_{rec}-Vx_{mc})/#sigma_{vx}; Entries (a.u.)", 1000, -10., 10.);
0240  TH1F *hPullVtxY = new TH1F("hPullVtxY", "Pull y position of MC vertex;(Vy_{rec}-Vy_{mc})/#sigma_{vy}; Entries (a.u.)", 1000, -10., 10.);
0241  TH1F *hPullVtxZ = new TH1F("hPullVtxZ", "Pull z position of MC vertex;(Vz_{rec}-Vz_{mc})/#sigma_{vz}; Entries (a.u.)", 2000, -20, 20);
0242  
0243   TH1F *hRes_SVx_Helixfit = new TH1F("hRes_SVx_Helixfit", "Fit method: Residual of SVx; SVx_{rec}-SVx_{mc} (mm); Entries (a.u.)", 200, -1.0, 1.0);
0244   TH1F *hRes_SVy_Helixfit= new TH1F("hRes_SVy_Helixfit", "Fit method: Residual of SVy; SVy_{rec}-SVy_{mc} (mm); Entries (a.u.)", 200, -1.0, 1.0);
0245   TH1F *hRes_SVz_Helixfit = new TH1F("hRes_SVz_Helixfit", "Fit method: Residual of SVz; SVz_{rec}-SVz_{mc} (mm); Entries (a.u.)", 1000, -5.0, 5.0);
0246   
0247   TH1F *hRes_SVx_Helixfit_pull = new TH1F("hRes_SVx_Helixfit_pull", "Fit method: Pull of SVx; SVx_{rec}-SVx_{mc}/#sigma; Entries (a.u.)", 200, -5.0, 5.0);
0248   TH1F *hRes_SVy_Helixfit_pull = new TH1F("hRes_SVy_Helixfit_pull", "Fit method: Pull of SVy; SVy_{rec}-SVy_{mc}/#sigma; Entries (a.u.)", 200, -5.0, 5.0);
0249   TH1F *hRes_SVz_Helixfit_pull = new TH1F("hRes_SVz_Helixfit_pull", "Fit method: Pull of SVz; SVz_{rec}-SVz_{mc}/#sigma; Entries (a.u.)", 200, -5.0, 5.0);
0250  
0251  
0252  TH1F *hchi2_vtx = new TH1F("hchi2_vtx", "Helix Calculation: Chi2/ndf; #chi^{2}/ndf; Entries (a.u.)", 1000, 0.0, 50.0); 
0253  TH1F *hchi2_vtx_sig = new TH1F("hchi2_vtx_sig", "Helix Calculation: Chi2/ndf; #chi^{2}/ndf; Entries (a.u.)", 1000, 0.0, 50.0); 
0254  TH1F *hchi2_vtx_bkg = new TH1F("hchi2_vtx_bkg", "Helix Calculation: Chi2/ndf; #chi^{2}/ndf; Entries (a.u.)", 1000, 0.0, 50.0); 
0255 
0256  TH2F *hLcpDecayVxVy = new TH2F("hLcpDecayVxVy", "#Lambda_{c}^{+} decay vertex to primary vertex;#Deltav_{x} (mm);#Deltav_{y} (mm)", 400, -1-0.0025, 1-0.0025, 400, -1-0.0025, 1-0.0025);
0257  TH2F *hLcpDecayVrVz = new TH2F("hLcpDecayVrVz", "#Lambda_{c}^{+} decay vertex to primary vertex;#Deltav_{z} (mm);#Deltav_{r} (mm)", 100, -2, 2, 100, -0.2, 1.8);
0258 
0259  TH2F *hMCLcpPtRap = new TH2F("hMCLcpPtRap", "MC #Lambda_{c}^{+};y;p_{T} (GeV/c)", 20, -5, 5, 100, 0, 10);
0260 
0261  TH2F *hMcPPtEta = new TH2F("hMcPPtEta", "MC P from #Lambda_{c}^{+} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0262  TH2F *hMcPPtEtaReco = new TH2F("hMcPPtEtaReco", "RC P from #Lambda_{c}^{+} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0263 
0264  TH2F *hMcKPtEta = new TH2F("hMcKPtEta", "MC K from #Lambda_{c}^{+} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0265  TH2F *hMcKPtEtaReco = new TH2F("hMcKPtEtaReco", "RC K from #Lambda_{c}^{+} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0266 
0267  TH2F *hMcPiPtEta = new TH2F("hMcPiPtEta", "MC #pi from #Lambda_{c}^{+} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0268  TH2F *hMcPiPtEtaReco = new TH2F("hMcPiPtEtaReco", "RC #pi from #Lambda_{c}^{+} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0269 
0270  TH1F *hNRecoVtx = new TH1F("hNRecoVtx", "Number of reconstructed vertices;N", 10, 0, 10);
0271 
0272  const char* part_name[3] = {"P", "K", "Pi"};
0273  const char* part_title[3] = {"P", "K", "#pi"};
0274  TH3F *hRcSecPartLocaToRCVtx[3];
0275  TH3F *hRcSecPartLocbToRCVtx[3];
0276  TH3F *hRcPrimPartLocaToRCVtx[3];
0277  TH3F *hRcPrimPartLocbToRCVtx[3];
0278  for(int i=0; i<3; i++)
0279  {
0280   hRcSecPartLocaToRCVtx[i] = new TH3F(Form("hRcSec%sLocaToRCVtx",part_name[i]), Form( "DCA_{xy} distribution for #Lambda_{c}^{+} decayed %s;p_{T} (GeV/c);#eta;DCA_{xy} (mm)", part_title[i]), 100, 0, 10, 20, -5, 5, 100, 0, 1);
0281   hRcSecPartLocbToRCVtx[i] = new TH3F(Form("hRcSec%sLocbToRCVtx",part_name[i]), Form( "DCA_{z} distribution for #Lambda_{c}^{+} decayed %s;p_{T} (GeV/c);#eta;DCA_{z} (mm)", part_title[i]), 100, 0, 10, 20, -5, 5, 100, -0.5, 0.5);
0282   hRcPrimPartLocaToRCVtx[i] = new TH3F(Form("hRcPrim%sLocaToRCVtx",part_name[i]), Form( "DCA_{xy} distribution for primary %s;p_{T} (GeV/c);#eta;DCA_{xy} (mm)", part_title[i]), 100, 0, 10, 20, -5, 5, 100, 0, 1);
0283   hRcPrimPartLocbToRCVtx[i] = new TH3F(Form("hRcPrim%sLocbToRCVtx",part_name[i]), Form( "DCA_{z} distribution for primary %s;p_{T} (GeV/c);#eta;DCA_{z} (mm)", part_title[i]), 100, 0, 10, 20, -5, 5, 100, -0.5, 0.5);
0284 }
0285 
0286 const char* axis_name[3] = {"x", "y", "z"};
0287 const int nDimDca = 4;
0288 const int nBinsDca[nDimDca] = {50, 20, 500, 50};
0289 const double minBinDca[nDimDca] = {0, -5, -1+0.002, 0};
0290 const double maxBinDca[nDimDca] = {5, 5, 1+0.002, 50};
0291 THnSparseF *hPrimTrkDcaToRCVtx[3][3];
0292 for(int i=0; i<3; i++)
0293 {
0294   for(int j=0; j<3; j++)
0295   {
0296    hPrimTrkDcaToRCVtx[i][j] = new THnSparseF(Form("hPrim%sDca%sToRCVtx",part_name[i],axis_name[j]), Form("DCA_{%s} distribution for primary %s;p_{T} (GeV/c);#eta;DCA_{%s} (mm);N_{MC}",axis_name[j],part_title[i],axis_name[j]), nDimDca, nBinsDca, minBinDca, maxBinDca);
0297  }
0298 }
0299 
0300 TH3F *h3PairDca12[2], *h3PairDca23[2], *h3PairDca13[2];
0301 TH3F *h3PairCosTheta[2];
0302 TH3F *h3PairDca[2];
0303 TH3F *h3PairDecayLength[2];
0304 const char* pair_name[2] = {"signal", "bkg"};
0305 const char* pair_title[2] = {"Signal", "Background"};
0306 for(int i=0; i<2; i++)
0307 {
0308   h3PairDca12[i] = new TH3F(Form("h3PairDca12_%s", pair_name[i]), Form("%s pair DCA_{12};p_{T} (GeV/c);#eta;DCA_{12} (mm)", pair_title[i]), 100, 0, 10, 20, -5, 5, 100, 0, 1);
0309   h3PairDca23[i] = new TH3F(Form("h3PairDca23_%s", pair_name[i]), Form("%s pair DCA_{23};p_{T} (GeV/c);#eta;DCA_{23} (mm)", pair_title[i]), 100, 0, 10, 20, -5, 5, 100, 0, 1);
0310   h3PairDca13[i] = new TH3F(Form("h3PairDca13_%s", pair_name[i]), Form("%s pair DCA_{13};p_{T} (GeV/c);#eta;DCA_{13} (mm)", pair_title[i]), 100, 0, 10, 20, -5, 5, 100, 0, 1);
0311   h3PairCosTheta[i] = new TH3F(Form("h3PairCosTheta_%s", pair_name[i]), Form("%s pair cos(#theta);p_{T} (GeV/c);#eta;cos(#theta)", pair_title[i]), 100, 0, 10, 20, -5, 5, 100, -1, 1);
0312 
0313   h3PairDca[i] = new TH3F(Form("h3PairDca_%s", pair_name[i]), Form("%s pair DCA;p_{T} (GeV/c);#eta;DCA_{pair} (mm)", pair_title[i]), 100, 0, 10, 20, -5, 5, 100, 0, 1);
0314 
0315   h3PairDecayLength[i] = new TH3F(Form("h3PairDecayLength_%s", pair_name[i]), Form("%s pair decay length;p_{T} (GeV/c);#eta;L (mm)", pair_title[i]), 100, 0, 10, 20, -5, 5, 100, 0, 1);
0316 }
0317 
0318   // Invariant mass
0319 const char* cut_name[2] = {"all", "DCA"};
0320 TH3F *h3InvMass[2][2];
0321 for(int i=0; i<2; i++)
0322 {
0323   for(int j=0; j<2; j++)
0324   {
0325    h3InvMass[i][j] = new TH3F(Form("h3InvMass_%s_%s", pair_name[i], cut_name[j]), "Invariant mass of unlike-sign #piK pairs;p_{T} (GeV/c);y;M_{#piK} (GeV/c^{2})", 100, 0, 10, 20, -5, 5, 100, 2.0, 3.0);
0326  }
0327 }
0328 
0329 TTreeReader treereader(chain);
0330   // MC  
0331 TTreeReaderArray<int> mcPartGenStatus = {treereader, "MCParticles.generatorStatus"};
0332 TTreeReaderArray<int> mcPartPdg = {treereader, "MCParticles.PDG"};
0333 TTreeReaderArray<float> mcPartCharge = {treereader, "MCParticles.charge"};
0334 TTreeReaderArray<unsigned int> mcPartParent_begin = {treereader, "MCParticles.parents_begin"};
0335 TTreeReaderArray<unsigned int> mcPartParent_end = {treereader, "MCParticles.parents_end"};
0336 TTreeReaderArray<int> mcPartParent_index = {treereader, "_MCParticles_parents.index"};
0337 TTreeReaderArray<unsigned int> mcPartDaughter_begin = {treereader, "MCParticles.daughters_begin"};
0338 TTreeReaderArray<unsigned int> mcPartDaughter_end = {treereader, "MCParticles.daughters_end"};
0339 TTreeReaderArray<int> mcPartDaughter_index = {treereader, "_MCParticles_daughters.index"};
0340 TTreeReaderArray<double> mcPartMass = {treereader, "MCParticles.mass"};
0341 TTreeReaderArray<double> mcPartVx = {treereader, "MCParticles.vertex.x"};
0342 TTreeReaderArray<double> mcPartVy = {treereader, "MCParticles.vertex.y"};
0343 TTreeReaderArray<double> mcPartVz = {treereader, "MCParticles.vertex.z"};
0344 TTreeReaderArray<float> mcMomPx = {treereader, "MCParticles.momentum.x"};
0345 TTreeReaderArray<float> mcMomPy = {treereader, "MCParticles.momentum.y"};
0346 TTreeReaderArray<float> mcMomPz = {treereader, "MCParticles.momentum.z"};
0347 TTreeReaderArray<double> mcEndPointX = {treereader, "MCParticles.endpoint.x"};
0348 TTreeReaderArray<double> mcEndPointY = {treereader, "MCParticles.endpoint.y"};
0349 TTreeReaderArray<double> mcEndPointZ = {treereader, "MCParticles.endpoint.z"};
0350 
0351 TTreeReaderArray<unsigned int> assocChSimID = {treereader, "ReconstructedChargedParticleAssociations.simID"};
0352 TTreeReaderArray<unsigned int> assocChRecID = {treereader, "ReconstructedChargedParticleAssociations.recID"};
0353 TTreeReaderArray<float> assocWeight = {treereader, "ReconstructedChargedParticleAssociations.weight"};
0354 
0355 TTreeReaderArray<float> rcMomPx = {treereader, "ReconstructedChargedParticles.momentum.x"};
0356 TTreeReaderArray<float> rcMomPy = {treereader, "ReconstructedChargedParticles.momentum.y"};
0357 TTreeReaderArray<float> rcMomPz = {treereader, "ReconstructedChargedParticles.momentum.z"};
0358 TTreeReaderArray<float> rcPosx = {treereader, "ReconstructedChargedParticles.referencePoint.x"};
0359 TTreeReaderArray<float> rcPosy = {treereader, "ReconstructedChargedParticles.referencePoint.y"};
0360 TTreeReaderArray<float> rcPosz = {treereader, "ReconstructedChargedParticles.referencePoint.z"};
0361 TTreeReaderArray<float> rcCharge = {treereader, "ReconstructedChargedParticles.charge"};
0362 TTreeReaderArray<int>   rcPdg = {treereader, "ReconstructedChargedParticles.PDG"};
0363 
0364 TTreeReaderArray<float> rcTrkLoca = {treereader, "CentralCKFTrackParameters.loc.a"};
0365 TTreeReaderArray<float> rcTrkLocb = {treereader, "CentralCKFTrackParameters.loc.b"};
0366 TTreeReaderArray<float> rcTrkqOverP = {treereader, "CentralCKFTrackParameters.qOverP"};
0367 TTreeReaderArray<float> rcTrkTheta = {treereader, "CentralCKFTrackParameters.theta"};
0368 TTreeReaderArray<float> rcTrkPhi = {treereader, "CentralCKFTrackParameters.phi"};
0369 
0370 rcMomPx2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.momentum.x"};
0371 rcMomPy2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.momentum.y"};
0372 rcMomPz2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.momentum.z"};
0373 rcCharge2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.charge"};
0374 
0375 rcTrkLoca2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.loc.a"};
0376 rcTrkLocb2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.loc.b"};
0377 rcTrkTheta2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.theta"};
0378 rcTrkPhi2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.phi"};
0379 rcTrkCov = new TTreeReaderArray<std::array<float, 21>>{treereader, "CentralCKFTrackParameters.covariance.covariance[21]"};
0380 
0381 TTreeReaderArray<float> CTVx = {treereader, "CentralTrackVertices.position.x"};
0382 TTreeReaderArray<float> CTVy = {treereader, "CentralTrackVertices.position.y"};
0383 TTreeReaderArray<float> CTVz = {treereader, "CentralTrackVertices.position.z"};
0384 TTreeReaderArray<int> CTVndf = {treereader, "CentralTrackVertices.ndf"};
0385 TTreeReaderArray<float> CTVchi2 = {treereader, "CentralTrackVertices.chi2"};
0386 TTreeReaderArray<float> CTVerr_xx = {treereader, "CentralTrackVertices.positionError.xx"};
0387 TTreeReaderArray<float> CTVerr_yy = {treereader, "CentralTrackVertices.positionError.yy"};
0388 TTreeReaderArray<float> CTVerr_zz = {treereader, "CentralTrackVertices.positionError.zz"};
0389 
0390 TTreeReaderArray<int> prim_vtx_index = {treereader, "PrimaryVertices_objIdx.index"};
0391 
0392 TTreeReaderArray<unsigned int> vtxAssocPart_begin = {treereader, "CentralTrackVertices.associatedParticles_begin"};
0393 TTreeReaderArray<unsigned int> vtxAssocPart_end = {treereader, "CentralTrackVertices.associatedParticles_end"};
0394 TTreeReaderArray<int> vtxAssocPart_index = {treereader, "_CentralTrackVertices_associatedParticles.index"};
0395 
0396   // File to create efficiency of D0 meson
0397 TFile *file_gen = new TFile("SignalLcpGen.root", "RECREATE");
0398 TTree *tree_gen = new TTree("treeMLSigGen", "treeMLSigGen"); 
0399 
0400    // Define variables to store in the Ntuple
0401 float pt_Lcp_gen, y_Lcp_gen;
0402 tree_gen->Branch("pt_Lcp_gen", &pt_Lcp_gen, "pt_Lcp_gen/F");  
0403 tree_gen->Branch("y_Lcp_gen", &y_Lcp_gen, "y_Lcp_gen/F"); 
0404 
0405 
0406     // Create a ROOT file to store the Ntuple
0407 TFile *file_signal = new TFile("SignalLcp.root", "RECREATE");
0408 TTree *tree_sig = new TTree("treeMLSig", "treeMLSig"); 
0409 
0410   // Define variables to store in the Ntuple
0411 float d0_p_sig, d0_k_sig, d0_pi_sig, d0xy_p_sig, d0xy_k_sig, d0xy_pi_sig, sum_d0xy_sig, dca_12_sig, dca_Lcp_sig, decay_length_sig;
0412 float costheta_sig, costhetaxy_sig, pt_Lcp_sig, y_Lcp_sig, mass_Lcp_sig, sigma_vtx_sig, mult_sig, chi2_sig;
0413 
0414     // Link the variables to the TTree branches
0415 tree_sig->Branch("d0_p", &d0_p_sig, "d0_p/F");
0416 tree_sig->Branch("d0_k", &d0_k_sig, "d0_k/F");
0417 tree_sig->Branch("d0_pi", &d0_pi_sig, "d0_pi/F"); 
0418 tree_sig->Branch("d0xy_p", &d0xy_p_sig, "d0xy_p/F");
0419 tree_sig->Branch("d0xy_k", &d0xy_k_sig, "d0xy_k/F");
0420 tree_sig->Branch("d0xy_pi", &d0xy_pi_sig, "d0xy_pi/F");  
0421 tree_sig->Branch("sum_d0xy", &sum_d0xy_sig, "sum_d0xy/F");        
0422 tree_sig->Branch("dca_12", &dca_12_sig, "dca_12/F");
0423 tree_sig->Branch("dca_Lcp", &dca_Lcp_sig, "dca_Lcp/F");
0424 tree_sig->Branch("pt_Lcp", &pt_Lcp_sig, "pt_Lcp/F");  
0425 tree_sig->Branch("y_Lcp", &y_Lcp_sig, "y_Lcp/F");  
0426 tree_sig->Branch("mass_Lcp", &mass_Lcp_sig, "mass_Lcp/F");              
0427 tree_sig->Branch("decay_length", &decay_length_sig, "decay_length/F");   
0428 tree_sig->Branch("costheta", &costheta_sig, "costheta/F"); 
0429 tree_sig->Branch("costheta_xy", &costhetaxy_sig, "costheta_xy/F"); 
0430 tree_sig->Branch("sigma_vtx", &sigma_vtx_sig, "sigma_vtx/F"); 
0431 tree_sig->Branch("mult", &mult_sig, "mult/F");
0432 tree_sig->Branch("chi2", &chi2_sig, "chi2/F");                          
0433 
0434 TFile *file_bkg = new TFile("BkgLcp.root", "RECREATE");
0435 TTree *tree_bkg = new TTree("treeMLBkg", "treeMLBkg"); 
0436 
0437   // Define variables to store in the Ntuple
0438 float d0_p_bkg, d0_k_bkg, d0_pi_bkg, d0xy_p_bkg, d0xy_k_bkg, d0xy_pi_bkg, sum_d0xy_bkg, dca_12_bkg, dca_Lcp_bkg, decay_length_bkg;
0439 float costheta_bkg, costhetaxy_bkg, pt_Lcp_bkg, y_Lcp_bkg, mass_Lcp_bkg, sigma_vtx_bkg, mult_bkg, chi2_bkg;
0440   // Link the variables to the TTree branches
0441 tree_bkg->Branch("d0_p", &d0_p_bkg, "d0_p/F");
0442 tree_bkg->Branch("d0_k", &d0_k_bkg, "d0_k/F");
0443 tree_bkg->Branch("d0_pi", &d0_pi_bkg, "d0_pi/F");
0444 tree_bkg->Branch("d0xy_p", &d0xy_p_bkg, "d0xy_p/F");
0445 tree_bkg->Branch("d0xy_k", &d0xy_k_bkg, "d0xy_k/F");
0446 tree_bkg->Branch("d0xy_pi", &d0xy_pi_bkg, "d0xy_pi/F");  
0447 tree_bkg->Branch("sum_d0xy", &sum_d0xy_bkg, "sum_d0xy/F");            
0448 tree_bkg->Branch("dca_12", &dca_12_bkg, "dca_12/F");
0449 tree_bkg->Branch("dca_Lcp", &dca_Lcp_bkg, "dca_Lcp/F");
0450 tree_bkg->Branch("pt_Lcp", &pt_Lcp_bkg, "pt_Lcp/F");  
0451 tree_bkg->Branch("y_Lcp", &y_Lcp_bkg, "y_Lcp/F");  
0452 tree_bkg->Branch("mass_Lcp", &mass_Lcp_bkg, "mass_Lcp/F");              
0453 tree_bkg->Branch("decay_length", &decay_length_bkg, "decay_length/F");   
0454 tree_bkg->Branch("costheta", &costheta_bkg, "costheta/F");  
0455 tree_bkg->Branch("costheta_xy", &costhetaxy_bkg, "costheta_xy/F"); 
0456 tree_bkg->Branch("sigma_vtx", &sigma_vtx_bkg, "sigma_vtx/F"); 
0457 tree_bkg->Branch("mult", &mult_bkg, "mult/F");
0458 tree_bkg->Branch("chi2", &chi2_bkg, "chi2/F");                                        
0459 
0460 
0461 int nevents = 0;
0462 int count = 0;
0463 
0464 while(treereader.Next())
0465 {
0466   if(nevents%1000==0) printf("\nEvent No.-----> %d\n",nevents);
0467   // find MC primary vertex
0468   int nMCPart = mcPartMass.GetSize();
0469 
0470   TVector3 vertex_mc(-999., -999., -999.);
0471   for(int imc=0; imc<nMCPart; imc++)
0472   {
0473    // Status 4 describes the incoming electron/proton and the end point of electron or proton is MC Vertex  
0474    if(mcPartGenStatus[imc] == 4 && mcPartPdg[imc] == 11)
0475    {
0476      vertex_mc.SetXYZ(mcEndPointX[imc], mcEndPointY[imc], mcEndPointZ[imc]);
0477      break;
0478    }
0479  }
0480  hEventStat->Fill(0.5);
0481  hMcVtxX->Fill(vertex_mc.x());
0482  hMcVtxY->Fill(vertex_mc.y());
0483  hMcVtxZ->Fill(vertex_mc.z());
0484 
0485 
0486       // get RC primary vertex
0487  TVector3 vertex_rc(-999., -999., -999.);
0488  TVector3 err_vertex_rc(-999., -999., -999.);
0489 
0490  if(prim_vtx_index.GetSize()>0)
0491  {
0492    int rc_vtx_index = prim_vtx_index[0];
0493    vertex_rc.SetXYZ(CTVx[rc_vtx_index], CTVy[rc_vtx_index], CTVz[rc_vtx_index]);
0494    err_vertex_rc.SetXYZ(sqrt(CTVerr_xx[rc_vtx_index]), sqrt(CTVerr_yy[rc_vtx_index]), sqrt(CTVerr_zz[rc_vtx_index]));
0495 
0496  }
0497  hRecVtxX->Fill(vertex_rc.x());
0498  hRecVtxY->Fill(vertex_rc.y());
0499  hRecVtxZ->Fill(vertex_rc.z());
0500 
0501   hPullVtxX->Fill((vertex_rc.x()-vertex_mc.x())/err_vertex_rc.x()); 
0502   hPullVtxY->Fill((vertex_rc.y()-vertex_mc.y())/err_vertex_rc.y()); 
0503   hPullVtxZ->Fill((vertex_rc.z()-vertex_mc.z())/err_vertex_rc.z());
0504 
0505 
0506       // map MC and RC particles
0507       int nAssoc = assocChRecID.GetSize();
0508       map<int, int> assoc_map_mc_to_rc;  // MC --> RC Associations
0509       map<int, int> assoc_map_rc_to_mc;  // RC --> MC Associations
0510 
0511       for(unsigned int rc_index=0; rc_index<rcMomPx.GetSize(); rc_index++)
0512       {
0513       // loop over the association to find the matched MC particle with largest weight
0514        double max_weight = 0;
0515        int matched_mc_index = -1;
0516        for(int j=0; j<nAssoc; j++)
0517        {
0518          if(assocChRecID[j] != rc_index) continue;
0519          if(assocWeight[j] > max_weight)
0520          {
0521           max_weight = assocWeight[j];
0522           matched_mc_index = assocChSimID[j];
0523         }
0524       }
0525 
0526       // build the map
0527       assoc_map_mc_to_rc[matched_mc_index] = rc_index;
0528       assoc_map_rc_to_mc[rc_index] = matched_mc_index;
0529     }
0530         // MC --> RC associations (kv; key,value for map)
0531    // std::cout << "MC to RC map:\n";
0532    // for (const auto& kv : assoc_map_mc_to_rc) {
0533    //     std::cout << "MC index " << kv.first << " ---> RC index " << kv.second << "\n";
0534    // }
0535 
0536     // RC --> MC associations
0537    // std::cout << "\nRC to MC map:\n";
0538    // for (const auto& kv : assoc_map_rc_to_mc) {
0539     //    std::cout << "RC index " << kv.first << " --->  MC index " << kv.second << "\n";
0540    // }
0541 
0542 
0543       // Loop over primary particles
0544     int nMcPart = 0;
0545     for(int imc=0; imc<nMCPart; imc++)
0546     {
0547      if(mcPartGenStatus[imc] == 1 && mcPartCharge[imc] != 0)
0548      {
0549        double dist = sqrt( pow(mcPartVx[imc]-vertex_mc.x(),2) + pow(mcPartVy[imc]-vertex_mc.y(),2) + pow(mcPartVz[imc]-vertex_mc.z(),2));      
0550        if(dist < 1e-4)
0551        {
0552           // count primary charged particles within |eta| < 3.5
0553         TVector3 mc_mom(mcMomPx[imc], mcMomPy[imc], mcMomPz[imc]);
0554         double mcEta = mc_mom.PseudoRapidity();
0555         if(fabs(mcEta) < 3.5) nMcPart++;
0556       }
0557     }
0558   }
0559   // Primary charged multiplicity at generated level with acceptance cut
0560   hMcMult->Fill(nMcPart);
0561 
0562 
0563   for(int imc=0; imc<nMCPart; imc++)
0564   {
0565    if(mcPartGenStatus[imc] == 1 && mcPartCharge[imc] != 0)
0566    {
0567      double dist = sqrt( pow(mcPartVx[imc]-vertex_mc.x(),2) + pow(mcPartVy[imc]-vertex_mc.y(),2) + pow(mcPartVz[imc]-vertex_mc.z(),2));      
0568      if(dist < 1e-4)
0569      {        
0570           // check if the MC particle is reconstructed
0571       int rc_index = -1;
0572       if(assoc_map_mc_to_rc.find(imc) != assoc_map_mc_to_rc.end()) rc_index = assoc_map_mc_to_rc[imc];
0573 
0574       if(rc_index>=0)
0575       {
0576         TVector3 dcaToVtx = GetDCAToPrimaryVertex(rc_index, vertex_rc);
0577 
0578         int ip = -1;
0579         if(fabs(mcPartPdg[imc]) == pdg_p) ip = 0;
0580         if(fabs(mcPartPdg[imc]) == pdg_k) ip = 1;
0581         if(fabs(mcPartPdg[imc]) == pdg_pi) ip = 2;
0582         if(ip>=0)
0583         {
0584          TVector3 mom(rcMomPx[rc_index], rcMomPy[rc_index], rcMomPz[rc_index]);
0585          if(ip<3)
0586          {
0587            hRcPrimPartLocaToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.Perp());
0588            hRcPrimPartLocbToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.z());
0589          }
0590 
0591          double fill1[] = {mom.Pt(), mom.Eta(), dcaToVtx.x(), nMcPart*1.};
0592          double fill2[] = {mom.Pt(), mom.Eta(), dcaToVtx.y(), nMcPart*1.};
0593          double fill3[] = {mom.Pt(), mom.Eta(), dcaToVtx.z(), nMcPart*1.};
0594          hPrimTrkDcaToRCVtx[ip][0]->Fill(fill1);
0595          hPrimTrkDcaToRCVtx[ip][1]->Fill(fill2);
0596          hPrimTrkDcaToRCVtx[ip][2]->Fill(fill3);
0597        }
0598      }
0599    }
0600  }
0601 }
0602 
0603  // Look for Λc⁺ → p K⁻ π⁺
0604 bool hasLc = false;
0605 vector<int> mc_index_Lcp_p;
0606 vector<int> mc_index_Lcp_k;
0607 vector<int> mc_index_Lcp_pi;
0608 mc_index_Lcp_p.clear();      
0609 mc_index_Lcp_k.clear();
0610 mc_index_Lcp_pi.clear();
0611 
0612 for(int imc=0; imc<nMCPart; imc++)
0613 {
0614  if(fabs(mcPartPdg[imc]) != pdg_Lcp) continue;
0615  hEventStat->Fill(1.5);
0616 
0617   // printf("MC Particle (Momentum) = (%f, %f, %f) \n",mcMomPx[imc], mcMomPy[imc], mcMomPz[imc]);           
0618  int nDaughters = mcPartDaughter_end[imc]-mcPartDaughter_begin[imc];
0619  if(nDaughters!=3) continue;
0620 
0621       // find Lc+ that decay into pKPi
0622  bool is_pkpi_decay = false;      
0623  int daug_index_1 = mcPartDaughter_index[mcPartDaughter_begin[imc]];
0624  int daug_index_2 = mcPartDaughter_index[mcPartDaughter_begin[imc]+1];
0625  int daug_index_3 = mcPartDaughter_index[mcPartDaughter_begin[imc]+2];    
0626  int daug_pdg_1 = mcPartPdg[daug_index_1];
0627  int daug_pdg_2 = mcPartPdg[daug_index_2];
0628  int daug_pdg_3 = mcPartPdg[daug_index_3];
0629 
0630  if( (fabs(daug_pdg_1)==pdg_p && fabs(daug_pdg_2)==pdg_k && fabs(daug_pdg_3)==pdg_pi) || (fabs(daug_pdg_1)==pdg_p && fabs(daug_pdg_2)==pdg_pi && fabs(daug_pdg_3)==pdg_k)  || 
0631    (fabs(daug_pdg_1)==pdg_k && fabs(daug_pdg_2)==pdg_p && fabs(daug_pdg_3)==pdg_pi) || (fabs(daug_pdg_1)==pdg_k && fabs(daug_pdg_2)==pdg_pi && fabs(daug_pdg_3)==pdg_p) ||
0632    (fabs(daug_pdg_1)==pdg_pi && fabs(daug_pdg_2)==pdg_k && fabs(daug_pdg_3)==pdg_p) || (fabs(daug_pdg_1)==pdg_pi && fabs(daug_pdg_2)==pdg_p && fabs(daug_pdg_3)==pdg_k) )
0633  {
0634    is_pkpi_decay = true;
0635  }
0636  if(!is_pkpi_decay) continue;
0637    // Efficiency for Lc+
0638  TLorentzVector mc_part_gen;
0639  mc_part_gen.SetXYZM(mcMomPx[imc], mcMomPy[imc], mcMomPz[imc], mcPartMass[imc]);
0640  // Fill the momentum at Gen level Aug3, 2025
0641  float mcRap_gen = mc_part_gen.Rapidity();
0642  float mcPt_gen = mc_part_gen.Pt();  
0643  pt_Lcp_gen = mcPt_gen;
0644  y_Lcp_gen = mcRap_gen;
0645  tree_gen->Fill();
0646 
0647  hEventStat->Fill(2.5);
0648 
0649  if((fabs(daug_pdg_1)==pdg_p && fabs(daug_pdg_2)==pdg_k && fabs(daug_pdg_3)==pdg_pi))
0650  {
0651    mc_index_Lcp_p.push_back(daug_index_1);
0652    mc_index_Lcp_k.push_back(daug_index_2);
0653    mc_index_Lcp_pi.push_back(daug_index_3);
0654  }
0655  else if((fabs(daug_pdg_1)==pdg_p && fabs(daug_pdg_2)==pdg_pi && fabs(daug_pdg_3)==pdg_k))
0656  {
0657    mc_index_Lcp_p.push_back(daug_index_1);
0658    mc_index_Lcp_k.push_back(daug_index_3);
0659    mc_index_Lcp_pi.push_back(daug_index_2);
0660  }
0661  else if((fabs(daug_pdg_1)==pdg_k && fabs(daug_pdg_2)==pdg_p && fabs(daug_pdg_3)==pdg_pi))
0662  {
0663    mc_index_Lcp_p.push_back(daug_index_2);
0664    mc_index_Lcp_k.push_back(daug_index_1);
0665    mc_index_Lcp_pi.push_back(daug_index_3);
0666  }
0667  else if((fabs(daug_pdg_1)==pdg_k && fabs(daug_pdg_2)==pdg_pi && fabs(daug_pdg_3)==pdg_p))
0668  {
0669    mc_index_Lcp_p.push_back(daug_index_3);
0670    mc_index_Lcp_k.push_back(daug_index_1);
0671    mc_index_Lcp_pi.push_back(daug_index_2);
0672  }
0673  else if((fabs(daug_pdg_1)==pdg_pi && fabs(daug_pdg_2)==pdg_k && fabs(daug_pdg_3)==pdg_p))
0674  {
0675    mc_index_Lcp_p.push_back(daug_index_3);
0676    mc_index_Lcp_k.push_back(daug_index_2);
0677    mc_index_Lcp_pi.push_back(daug_index_1);
0678  }                              
0679  else 
0680    {
0681      mc_index_Lcp_p.push_back(daug_index_2);
0682      mc_index_Lcp_k.push_back(daug_index_3);
0683      mc_index_Lcp_pi.push_back(daug_index_1);
0684    }
0685    hasLc = true;
0686 
0687       // Lcp kinematics
0688    TLorentzVector mc_mom_vec;
0689    mc_mom_vec.SetXYZM(mcMomPx[imc], mcMomPy[imc], mcMomPz[imc], mcPartMass[imc]);
0690 //printf("Mass  = %f \n",mcPartMass[imc]);
0691 
0692    double mcRap = mc_mom_vec.Rapidity();
0693    double mcPt = mc_mom_vec.Pt();
0694    hMCLcpPtRap->Fill(mcRap, mcPt); 
0695 
0696       // decay daughter kinematics
0697    for(int ip = 0; ip<3; ip++)
0698    {
0699      int mc_part_index;
0700      if(ip==0) mc_part_index = mc_index_Lcp_p[mc_index_Lcp_p.size()-1];
0701      if(ip==1) mc_part_index = mc_index_Lcp_k[mc_index_Lcp_k.size()-1];
0702      if(ip==2) mc_part_index = mc_index_Lcp_pi[mc_index_Lcp_pi.size()-1];
0703 
0704      TLorentzVector mc_part_vec;
0705      mc_part_vec.SetXYZM(mcMomPx[mc_part_index], mcMomPy[mc_part_index], mcMomPz[mc_part_index], mcPartMass[mc_part_index]);
0706      if(ip==0) hMcPPtEta->Fill(mc_part_vec.Eta(), mc_part_vec.Pt());
0707      if(ip==1) hMcKPtEta->Fill(mc_part_vec.Eta(), mc_part_vec.Pt());
0708      if(ip==2) hMcPiPtEta->Fill(mc_part_vec.Eta(), mc_part_vec.Pt());         
0709 
0710      int rc_part_index = -1;
0711      if(assoc_map_mc_to_rc.find(mc_part_index) != assoc_map_mc_to_rc.end()) rc_part_index = assoc_map_mc_to_rc[mc_part_index];
0712      if (debug) printf("Rec: ProngNo., MC Index, Reco Index = (%d, %d, %d) \n",ip,mc_part_index,rc_part_index);
0713      if(rc_part_index>=0)
0714      {
0715       TVector3 dcaToVtx = GetDCAToPrimaryVertex(rc_part_index, vertex_rc);
0716       TVector3 mom(rcMomPx[rc_part_index], rcMomPy[rc_part_index], rcMomPz[rc_part_index]);
0717       hRcSecPartLocaToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.Pt());
0718       hRcSecPartLocbToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.z());
0719 
0720           //printf("Sec %d: (%2.4f, %2.4f, %2.4f), mcStartPoint = (%2.4f, %2.4f, %2.4f)\n", rc_part_index, pos.x(), pos.y(), pos.z(), mcPartVx[mc_part_index], mcPartVy[mc_part_index], mcPartVz[mc_part_index]);
0721     }
0722   }
0723 }
0724 
0725       // Get reconstructed protons, kaons and pions
0726       hNRecoVtx->Fill(CTVx.GetSize());
0727       const int pid_mode = 0; // 0 - truth; 1 - realistic
0728       vector<unsigned int> p_index;
0729       vector<unsigned int> k_index;      
0730       vector<unsigned int> pi_index;
0731       p_index.clear();
0732       k_index.clear();      
0733       pi_index.clear();
0734 
0735       for(unsigned int rc_index=0; rc_index<rcMomPx.GetSize(); rc_index++)
0736       {   
0737        if(pid_mode==0)
0738        {
0739          int iSimPartID = -1;
0740          if(assoc_map_rc_to_mc.find(rc_index) != assoc_map_rc_to_mc.end()) iSimPartID = assoc_map_rc_to_mc[rc_index];
0741          if(iSimPartID>=0)
0742          {
0743           if(fabs(mcPartPdg[iSimPartID]) == pdg_p) p_index.push_back(rc_index); 
0744           if(fabs(mcPartPdg[iSimPartID]) == pdg_k) k_index.push_back(rc_index);
0745           if(fabs(mcPartPdg[iSimPartID]) == pdg_pi) pi_index.push_back(rc_index);
0746         }
0747       }
0748       else if(pid_mode==1)
0749       {
0750        if(fabs(rcPdg[rc_index]) == pdg_p) p_index.push_back(rc_index);
0751        if(fabs(rcPdg[rc_index]) == pdg_k) k_index.push_back(rc_index);
0752        if(fabs(rcPdg[rc_index]) == pdg_pi) pi_index.push_back(rc_index);          
0753      }
0754    }
0755 
0756   // printf("Proton, Kaon, and Pion vector sizes: = (%d, %d, %d) \n",p_index.size(), k_index.size(),pi_index.size()); 
0757       // Combinatorics proton, kaon, and pion 
0758    for(unsigned int i=0; i<p_index.size(); i++)
0759    {
0760      TVector3 dcaToVtx_p = GetDCAToPrimaryVertex(p_index[i], vertex_rc);
0761      int q_proton = rcCharge[p_index[i]];
0762 
0763      for(unsigned int j=0; j<k_index.size(); j++)
0764      {
0765       TVector3 dcaToVtx_k = GetDCAToPrimaryVertex(k_index[j], vertex_rc);
0766       int q_kaon = rcCharge[k_index[j]];
0767 
0768       for(unsigned int k=0; k<pi_index.size(); k++)
0769       {
0770        TVector3 dcaToVtx_pi = GetDCAToPrimaryVertex(pi_index[k], vertex_rc);
0771        int q_pion = rcCharge[pi_index[k]];
0772 
0773        if ((q_proton == +1 && q_kaon == -1 && q_pion == +1) || (q_proton == -1 && q_kaon == +1 && q_pion == -1))       
0774        {
0775         bool is_Lcp_pkpi = false;
0776 
0777         int  mc_index_p = -1, mc_index_k = -1, mc_index_pi = -1;
0778         if(assoc_map_rc_to_mc.find(p_index[i]) != assoc_map_rc_to_mc.end()) mc_index_p = assoc_map_rc_to_mc[p_index[i]];
0779         if(assoc_map_rc_to_mc.find(k_index[j])  != assoc_map_rc_to_mc.end()) mc_index_k  = assoc_map_rc_to_mc[k_index[j]];
0780         if(assoc_map_rc_to_mc.find(pi_index[k])  != assoc_map_rc_to_mc.end()) mc_index_pi  = assoc_map_rc_to_mc[pi_index[k]];
0781 
0782         for(unsigned int idaugh=0; idaugh<mc_index_Lcp_pi.size(); idaugh++)
0783         {
0784           if(mc_index_p==mc_index_Lcp_p[idaugh] && mc_index_k==mc_index_Lcp_k[idaugh] && mc_index_pi==mc_index_Lcp_pi[idaugh])
0785           {
0786            is_Lcp_pkpi = true;
0787            break;
0788          }
0789        }
0790 
0791        float dcaDaughters[3], cosTheta, decayLength, V0DcaToVtx, cosTheta_xy, sigma_vtx;
0792        float chi2_ndf = 0.;
0793        TVector3 decayVertex;
0794        double err_Par[6];  // or whatever size is appropriate
0795        TLorentzVector parent = GetDCADaughters(p_index[i], k_index[j], pi_index[k], vertex_rc, dcaDaughters, cosTheta, cosTheta_xy, decayLength, V0DcaToVtx, sigma_vtx, chi2_ndf, decayVertex,err_Par);
0796            hchi2_vtx->Fill(chi2_ndf);
0797 
0798        if(is_Lcp_pkpi)
0799        {
0800         hEventStat->Fill(3.5);
0801         hchi2_vtx_sig->Fill(chi2_ndf);
0802            TVector3 MCVertex_Kaon(mcPartVx[mc_index_k], mcPartVy[mc_index_k], mcPartVz[mc_index_k]);
0803            TVector3 MCVertex_Pion(mcPartVx[mc_index_pi], mcPartVy[mc_index_pi], mcPartVz[mc_index_pi]);
0804            
0805          //  printf("Signal MC Vertex Kaon = (%f, %f, %f)\n",MCVertex_Kaon.X(), MCVertex_Kaon.Y(), MCVertex_Kaon.Z());      
0806           // printf("Signal MC Vertex Pion = (%f, %f, %f)\n",MCVertex_Pion.X(), MCVertex_Pion.Y(), MCVertex_Pion.Z());
0807            
0808           // cout<<"Signal MC Vertex Kaon (cm): "<<MCVertex_Kaon.X()*0.1<<"\t"<<MCVertex_Kaon.Y()*0.1<<"\t"<<MCVertex_Kaon.Z()*0.1<<endl;
0809           // cout<<"Signal MC Vertex Pion (cm): "<<MCVertex_Pion.X()*0.1<<"\t"<<MCVertex_Pion.Y()*0.1<<"\t"<<MCVertex_Pion.Z()*0.1<<endl; 
0810            
0811              hRes_SVx_Helixfit->Fill((decayVertex.X()-MCVertex_Kaon.X()*0.1)*10);
0812               hRes_SVy_Helixfit->Fill((decayVertex.Y()-MCVertex_Kaon.Y()*0.1)*10);
0813               hRes_SVz_Helixfit->Fill((decayVertex.Z()-MCVertex_Kaon.Z()*0.1)*10);  
0814               
0815               hRes_SVx_Helixfit_pull->Fill(((decayVertex.X()-MCVertex_Kaon.X()*0.1))/err_Par[0]);
0816               hRes_SVy_Helixfit_pull->Fill(((decayVertex.Y()-MCVertex_Kaon.Y()*0.1))/err_Par[1]);
0817               hRes_SVz_Helixfit_pull->Fill(((decayVertex.Z()-MCVertex_Kaon.Z()*0.1))/err_Par[2]);
0818                       
0819         if (q_proton == 1 && q_kaon == -1 && q_pion == 1)
0820         hEventStat->Fill(4.5);   // Λc⁺
0821       else if (q_proton == -1 && q_kaon == 1 && q_pion == -1)
0822         hEventStat->Fill(5.5);   // Λc⁻
0823 
0824       h3PairDca12[0]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters[0]);
0825       h3PairDca23[0]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters[1]); 
0826       h3PairDca13[0]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters[2]);              
0827       h3PairCosTheta[0]->Fill(parent.Pt(), parent.Rapidity(), cosTheta);
0828       h3PairDca[0]->Fill(parent.Pt(), parent.Rapidity(), V0DcaToVtx);
0829       h3PairDecayLength[0]->Fill(parent.Pt(), parent.Rapidity(), decayLength);
0830               //printf("Signal: dca12 = %2.4f, cosTheta = %2.4f, D0dca = %2.4f, decay = %2.4f\n", dcaDaughters, cosTheta, V0DcaToVtx, decayLength);
0831       h3InvMass[0][0]->Fill(parent.Pt(), parent.Rapidity(), parent.M());
0832 
0833                 // Toplogical Variables for Signal
0834       d0_p_sig = dcaToVtx_p.Mag();
0835       d0_k_sig = dcaToVtx_k.Mag();
0836       d0_pi_sig = dcaToVtx_pi.Mag();              
0837       d0xy_p_sig = dcaToVtx_p.Perp();
0838       d0xy_k_sig = dcaToVtx_k.Perp();
0839       d0xy_pi_sig = dcaToVtx_pi.Perp();           
0840       sum_d0xy_sig = sqrt(d0xy_p_sig*d0xy_p_sig+d0xy_k_sig*d0xy_k_sig+d0xy_pi_sig*d0xy_pi_sig);           
0841       dca_12_sig = *min_element(dcaDaughters, dcaDaughters + 3);
0842       dca_Lcp_sig = V0DcaToVtx;
0843       decay_length_sig = decayLength;
0844       costheta_sig = cosTheta;
0845       costhetaxy_sig = cosTheta_xy;
0846       pt_Lcp_sig = parent.Pt();
0847       y_Lcp_sig = parent.Rapidity();
0848       mass_Lcp_sig = parent.M();
0849       sigma_vtx_sig = sigma_vtx;
0850       mult_sig = nMcPart;
0851       chi2_sig = chi2_ndf;
0852       tree_sig->Fill();  
0853 
0854     }
0855     else
0856     {
0857       hEventStat->Fill(6.5);
0858       hchi2_vtx_bkg->Fill(chi2_ndf);              
0859       h3PairDca12[1]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters[0]);
0860       h3PairDca23[1]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters[1]); 
0861       h3PairDca13[1]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters[2]);                  
0862       h3PairCosTheta[1]->Fill(parent.Pt(), parent.Rapidity(), cosTheta);
0863       h3PairDca[1]->Fill(parent.Pt(), parent.Rapidity(), V0DcaToVtx);
0864       h3PairDecayLength[1]->Fill(parent.Pt(), parent.Rapidity(), decayLength);
0865 
0866               //printf("Bkg: dca12 = %2.4f, cosTheta = %2.4f, D0dca = %2.4f, decay = %2.4f\n", dcaDaughters, cosTheta, V0DcaToVtx, decayLength);
0867       h3InvMass[1][0]->Fill(parent.Pt(), parent.Rapidity(), parent.M());
0868 
0869               // Toplogical Variables for Bkg
0870       d0_p_bkg = dcaToVtx_p.Mag();
0871       d0_k_bkg = dcaToVtx_k.Mag();
0872       d0_pi_bkg = dcaToVtx_pi.Mag();              
0873       d0xy_p_bkg = dcaToVtx_p.Perp();
0874       d0xy_k_bkg = dcaToVtx_k.Perp();
0875       d0xy_pi_bkg = dcaToVtx_pi.Perp();           
0876       sum_d0xy_bkg = sqrt(d0xy_p_bkg*d0xy_p_bkg+d0xy_k_bkg*d0xy_k_bkg+d0xy_pi_bkg*d0xy_pi_bkg);       
0877       dca_12_bkg = *min_element(dcaDaughters, dcaDaughters + 3);
0878       dca_Lcp_bkg = V0DcaToVtx;
0879       decay_length_bkg = decayLength;
0880       costheta_bkg = cosTheta;
0881       costhetaxy_bkg = cosTheta_xy;
0882       pt_Lcp_bkg = parent.Pt();
0883       y_Lcp_bkg = parent.Rapidity();
0884       mass_Lcp_bkg = parent.M();
0885       sigma_vtx_bkg = sigma_vtx;
0886       mult_bkg = nMcPart;
0887       chi2_bkg = chi2_ndf;
0888       tree_bkg->Fill();
0889     } 
0890 
0891   }
0892 }
0893 }
0894 }
0895 
0896 nevents++;
0897 }
0898 
0899 
0900 file_gen->cd();  
0901 tree_gen->Write();
0902 file_gen->Close();
0903 
0904 file_signal->cd();  
0905 tree_sig->Write();
0906 file_signal->Close();  
0907 
0908 file_bkg->cd();  
0909 tree_bkg->Write();
0910 file_bkg->Close(); 
0911 
0912 TFile *outfile = new TFile(outname.Data(), "recreate");
0913 hEventStat->SetMarkerSize(2);
0914 hEventStat->Write("");
0915 hMcMult->Write();
0916 hMcVtxX->Write();
0917 hMcVtxY->Write();
0918 hMcVtxZ->Write();
0919 hRecVtxX->Write();
0920 hRecVtxY->Write();
0921 hRecVtxZ->Write();
0922 hPullVtxX->Write();
0923 hPullVtxY->Write();
0924 hPullVtxZ->Write();
0925 hchi2_vtx->Write();
0926 hchi2_vtx_sig->Write();           
0927 hchi2_vtx_bkg->Write();
0928 hLcpDecayVxVy->Write();
0929 hLcpDecayVrVz->Write();
0930 hRes_SVx_Helixfit->Write();
0931 hRes_SVy_Helixfit->Write();
0932 hRes_SVz_Helixfit->Write();
0933 hRes_SVx_Helixfit_pull->Write();
0934 hRes_SVy_Helixfit_pull->Write(); 
0935 hRes_SVz_Helixfit_pull->Write(); 
0936 hMCLcpPtRap->Write();
0937 hMcPPtEta->Write();
0938 hMcPPtEtaReco->Write();
0939 hMcPiPtEta->Write();
0940 hMcPiPtEtaReco->Write();
0941 hMcKPtEta->Write();
0942 hMcKPtEtaReco->Write();
0943 
0944 hNRecoVtx->Write();
0945 
0946 for(int ip=0; ip<3; ip++)
0947 {
0948   hRcSecPartLocaToRCVtx[ip]->Write();
0949   hRcSecPartLocbToRCVtx[ip]->Write();
0950   hRcPrimPartLocaToRCVtx[ip]->Write();
0951   hRcPrimPartLocbToRCVtx[ip]->Write();
0952 }
0953 
0954 for(int i=0; i<3; i++)
0955 {
0956   for(int j=0; j<3; j++)
0957   {
0958    hPrimTrkDcaToRCVtx[i][j]->Write();
0959  }
0960 }
0961 
0962 for(int i=0; i<2; i++)
0963 {
0964   h3PairDca12[i]->Write();
0965   h3PairDca23[i]->Write();
0966   h3PairDca13[i]->Write();            
0967   h3PairCosTheta[i]->Write();
0968   h3PairDca[i]->Write();
0969   h3PairDecayLength[i]->Write();
0970 }
0971 
0972 for(int i=0; i<2; i++)
0973 {
0974   for(int j=0; j<2; j++)
0975   {
0976    h3InvMass[i][j]->Write();
0977  }
0978 }
0979 
0980 
0981 outfile->Close();
0982 
0983 }
0984 
0985 //======================================
0986 TVector3 GetDCAToPrimaryVertex(const int index, TVector3 vtx)
0987 {
0988   //printf("check %d: (%2.4f, %2.4f, %2.4f) =? (%2.4f, %2.4f, %2.4f)\n", index, mom.x(), mom.y(), mom.z(), rcMomPx2->At(index), rcMomPy2->At(index), rcMomPz2->At(index));
0989 
0990   TVector3 pos(rcTrkLoca2->At(index) * sin(rcTrkPhi2->At(index)) * -1 * millimeter, rcTrkLoca2->At(index) * cos(rcTrkPhi2->At(index)) * millimeter, rcTrkLocb2->At(index) * millimeter);
0991   TVector3 mom(rcMomPx2->At(index), rcMomPy2->At(index), rcMomPz2->At(index));
0992 
0993   StPhysicalHelix pHelix(mom, pos, bField * tesla, rcCharge2->At(index));
0994 
0995   TVector3 vtx_tmp;
0996   vtx_tmp.SetXYZ(vtx.x()*millimeter, vtx.y()*millimeter, vtx.z()*millimeter);
0997   
0998   pHelix.moveOrigin(pHelix.pathLength(vtx_tmp));
0999   TVector3 dcaToVtx = pHelix.origin() - vtx_tmp;
1000 
1001   dcaToVtx.SetXYZ(dcaToVtx.x()/millimeter, dcaToVtx.y()/millimeter, dcaToVtx.z()/millimeter);
1002   
1003   return dcaToVtx;
1004 }
1005 
1006 //======================================
1007 // Local to global conversion 
1008 // x = - l0 Sin (phi), y = l0 Cos (phi), z = l1 
1009 TLorentzVector GetDCADaughters(const int index1, const int index2, const int index3, TVector3 vtx,
1010   float *dcaDaughters, float &cosTheta, float &cosTheta_xy, float &decayLength, float &V0DcaToVtx, float &sigma_vtx, float & chi2_ndf, TVector3 &decayVertex, double *parFitErr)
1011 {
1012   // -- get helix
1013   TVector3 pos1(rcTrkLoca2->At(index1) * sin(rcTrkPhi2->At(index1)) * -1 * millimeter, rcTrkLoca2->At(index1) * cos(rcTrkPhi2->At(index1)) * millimeter, rcTrkLocb2->At(index1) * millimeter);
1014   TVector3 pos2(rcTrkLoca2->At(index2) * sin(rcTrkPhi2->At(index2)) * -1 * millimeter, rcTrkLoca2->At(index2) * cos(rcTrkPhi2->At(index2)) * millimeter, rcTrkLocb2->At(index2) * millimeter);
1015   TVector3 pos3(rcTrkLoca2->At(index3) * sin(rcTrkPhi2->At(index3)) * -1 * millimeter, rcTrkLoca2->At(index3) * cos(rcTrkPhi2->At(index3)) * millimeter, rcTrkLocb2->At(index3) * millimeter);
1016   
1017   TVector3 mom1(rcMomPx2->At(index1), rcMomPy2->At(index1), rcMomPz2->At(index1));
1018   TVector3 mom2(rcMomPx2->At(index2), rcMomPy2->At(index2), rcMomPz2->At(index2));
1019   TVector3 mom3(rcMomPx2->At(index3), rcMomPy2->At(index3), rcMomPz2->At(index3));
1020 
1021   float charge1 = rcCharge2->At(index1);
1022   float charge2 = rcCharge2->At(index2);
1023   float charge3 = rcCharge2->At(index3);  
1024   
1025   StPhysicalHelix p1Helix(mom1, pos1, bField * tesla, charge1);
1026   StPhysicalHelix p2Helix(mom2, pos2, bField * tesla, charge2);
1027   StPhysicalHelix p3Helix(mom3, pos3, bField * tesla, charge3);
1028 
1029   TVector3 vtx_tmp;
1030   vtx_tmp.SetXYZ(vtx.x()*millimeter, vtx.y()*millimeter, vtx.z()*millimeter);
1031   
1032   // Decay vertex with distance minimization
1033   double s1, s2, s3;
1034   
1035   getDecayVertex_chi2fit(index1,index2,index3,s1,s2,s3,decayVertex,chi2_ndf, parFitErr);
1036   
1037   // Points of closest approach among three helices
1038   TVector3 const p1 = p1Helix.at(s1);
1039   TVector3 const p2 = p2Helix.at(s2);
1040   TVector3 const p3 = p3Helix.at(s3);
1041 
1042   printf("Chi2/ndf = %f \n",chi2_ndf);
1043  // Three dcaDaughters; dcaDaughters_12, dcaDaughters_23, dcaDaughters_13
1044   dcaDaughters[0] = (p1 - p2).Mag()/millimeter;
1045   dcaDaughters[1] = (p2 - p3).Mag()/millimeter;
1046   dcaDaughters[2] = (p3 - p1).Mag()/millimeter; 
1047 
1048   // -- calculate Momentum of each daughters at the DCA point
1049   TVector3 const p1MomAtDca = p1Helix.momentumAt(s1,  bField * tesla);  
1050   TVector3 const p2MomAtDca = p2Helix.momentumAt(s2,  bField * tesla); 
1051   TVector3 const p3MomAtDca = p3Helix.momentumAt(s3,  bField * tesla);
1052   
1053   TLorentzVector p1FourMom(p1MomAtDca, sqrt(p1MomAtDca.Mag2()+gProtonMass*gProtonMass));
1054   TLorentzVector p2FourMom(p2MomAtDca, sqrt(p2MomAtDca.Mag2()+gKaonMass*gKaonMass));
1055   TLorentzVector p3FourMom(p3MomAtDca, sqrt(p3MomAtDca.Mag2()+gPionMass*gPionMass));  
1056   
1057   TLorentzVector parent = p1FourMom + p2FourMom + p3FourMom;
1058 
1059   // Confirm by evaluating the distance with vertices
1060 
1061   sigma_vtx = sqrt((p1-decayVertex).Mag2()+(p2-decayVertex).Mag2()+(p3-decayVertex).Mag2())/millimeter;
1062 
1063   // -- calculate pointing angle and decay length with respect to primary vertex
1064   //    if decay vertex is a tertiary vertex
1065   //    -> only rough estimate -> needs to be updated after secondary vertex is found
1066   TVector3 vtxToV0 = decayVertex - vtx_tmp;
1067   TVector3 vtxToV0_xy(vtxToV0.x(), vtxToV0.y(), 0.);
1068   TVector3 parent_xy(parent.Vect().x(),parent.Vect().y(),0.);
1069   float pointingAngle = vtxToV0.Angle(parent.Vect());
1070   float pointingAngle_xy = vtxToV0_xy.Angle(parent_xy);
1071   cosTheta = std::cos(pointingAngle);
1072   cosTheta_xy = std::cos(pointingAngle_xy);  
1073   decayLength = vtxToV0.Mag()/millimeter;
1074 
1075   // -- calculate V0 DCA to primary vertex
1076   V0DcaToVtx = decayLength * std::sin(pointingAngle);
1077 
1078   //TVector3 dcaToVtx = GetDCAToPrimaryVertex(parent.Vect(), decayVertex, 0, vtx);
1079   //V0DcaToVtx = dcaToVtx.Mag();
1080 
1081   return parent;
1082 }
1083 
1084 void getDecayVertex_chi2fit(const int index1, const int index2, const int index3,  double &s1, double &s2, double &s3, TVector3 &vertex, float &chi2, double *parFitErr)
1085 {
1086     TVector3 pos1(rcTrkLoca2->At(index1) * sin(rcTrkPhi2->At(index1)) * -1 * millimeter,
1087                   rcTrkLoca2->At(index1) * cos(rcTrkPhi2->At(index1)) * millimeter,
1088                   rcTrkLocb2->At(index1) * millimeter);
1089 
1090     TVector3 mom1(rcMomPx2->At(index1), rcMomPy2->At(index1), rcMomPz2->At(index1));
1091 
1092     TVector3 pos2(rcTrkLoca2->At(index2) * sin(rcTrkPhi2->At(index2)) * -1 * millimeter,
1093                   rcTrkLoca2->At(index2) * cos(rcTrkPhi2->At(index2)) * millimeter,
1094                   rcTrkLocb2->At(index2) * millimeter);
1095 
1096     TVector3 mom2(rcMomPx2->At(index2), rcMomPy2->At(index2), rcMomPz2->At(index2));
1097     
1098     TVector3 pos3(rcTrkLoca2->At(index3) * sin(rcTrkPhi2->At(index3)) * -1 * millimeter,
1099                   rcTrkLoca2->At(index3) * cos(rcTrkPhi2->At(index3)) * millimeter,
1100                   rcTrkLocb2->At(index3) * millimeter);
1101 
1102     TVector3 mom3(rcMomPx2->At(index3), rcMomPy2->At(index3), rcMomPz2->At(index3));
1103 
1104 
1105     float charge1 = rcCharge2->At(index1);
1106     float charge2 = rcCharge2->At(index2);
1107     float charge3 = rcCharge2->At(index3);
1108 
1109 
1110     StPhysicalHelix helix1(mom1, pos1, bField * tesla, charge1);
1111     StPhysicalHelix helix2(mom2, pos2, bField * tesla, charge2);
1112     StPhysicalHelix helix3(mom3, pos3, bField * tesla, charge3); 
1113     
1114     pair<double, double> const ss12 = helix1.pathLengths(helix2);
1115     pair<double, double> const ss13 = helix1.pathLengths(helix3);
1116     TVector3 const p1_init = helix1.at(ss12.first);
1117     TVector3 const p2_init = helix2.at(ss12.second); 
1118     TVector3 const p3_init = helix3.at(ss13.second);
1119     TVector3 const centroid = 1./3*(p1_init+p2_init+p3_init); 
1120     
1121     std::array<float, 21>& fcov1 = rcTrkCov->At(index1);
1122     std::array<float, 21>& fcov2 = rcTrkCov->At(index2); 
1123     std::array<float, 21>& fcov3 = rcTrkCov->At(index3);  
1124 
1125       // Perform Minimization
1126     const Int_t nPar = 6;
1127     Chi2Minimization d2Function(helix1,helix2,helix3,fcov1,fcov2,fcov3);
1128     ROOT::Math::Functor fcn(d2Function,nPar); // 5 parameters
1129     ROOT::Fit::Fitter fitter;
1130 
1131     double pStart[nPar] = {centroid.X(),centroid.Y(),centroid.Z(),ss12.first,ss12.second,ss13.second};
1132     fitter.SetFCN(fcn, pStart,nPar,1);
1133     fitter.Config().ParSettings(0).SetName("x0");
1134     fitter.Config().ParSettings(0).SetStepSize(0.01);
1135    // fitter.Config().ParSettings(0).SetLimits(-1., 1.);    
1136     // No limits for x, y, z
1137 
1138     fitter.Config().ParSettings(1).SetName("y0");
1139     fitter.Config().ParSettings(1).SetStepSize(0.01);
1140    // fitter.Config().ParSettings(1).SetLimits(-1., 1.);
1141 
1142     fitter.Config().ParSettings(2).SetName("z0");
1143     fitter.Config().ParSettings(2).SetStepSize(0.01);
1144    // fitter.Config().ParSettings(2).SetLimits(-10., 10.);    
1145 
1146     fitter.Config().ParSettings(3).SetName("s1");
1147     fitter.Config().ParSettings(3).SetValue(0.0);
1148     fitter.Config().ParSettings(3).SetStepSize(0.01);
1149    // fitter.Config().ParSettings(3).SetLimits(-1., 1.);
1150 
1151     fitter.Config().ParSettings(4).SetName("s2");
1152     fitter.Config().ParSettings(4).SetValue(0.0);
1153     fitter.Config().ParSettings(4).SetStepSize(0.01);
1154    // fitter.Config().ParSettings(4).SetLimits(-1., 1.);  
1155 
1156     fitter.Config().ParSettings(5).SetName("s3");
1157     fitter.Config().ParSettings(5).SetValue(0.0);
1158     fitter.Config().ParSettings(5).SetStepSize(0.01);
1159    // fitter.Config().ParSettings(5).SetLimits(-1., 1.);          
1160     
1161     fitter.Config().MinimizerOptions().SetMaxIterations(10000); 
1162     // do the fit 
1163 
1164     Bool_t ok = fitter.FitFCN();
1165     if (!ok) Error("Fitting","Fitting failed");
1166     const ROOT::Fit::FitResult & result = fitter.Result();
1167     //chi2 = fitter.Result().Chi2()/3.0;
1168     chi2 = fitter.Result().MinFcnValue()/3.0;  // Minimum value of your function
1169     int status = fitter.Result().Status();
1170     if (status>0 ) {printf("Fit Failed!!!!\n");}
1171    // if (status>0 || chi2_ndf>10. ) return;
1172    // cout <<"\033[1;31m Fit Result Chi2:\033[0m"<<chi2<<endl;
1173     //result.Print(std::cout);
1174    
1175    // Get the covariance matrix
1176  //  TMatrixDSym covMatrix(5);
1177   // result.GetCovarianceMatrix(covMatrix); // Matrix for the parameter errors
1178   // covMatrix.Print();
1179    
1180    const double * parFit = result.GetParams();
1181    const double *FitErr = result.GetErrors();
1182    for (int i = 0; i < nPar; ++i) parFitErr[i] = FitErr[i];
1183 
1184     vertex.SetXYZ(parFit[0], parFit[1], parFit[2]);
1185     s1 = parFit[3]; s2 = parFit[4]; s3 = parFit[5];
1186 
1187 
1188 }
1189
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!