Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-08-28 09:14:09 

0001 // Λc⁺ → p K⁻ π⁺ reconstruction with (6.28 ± 0.32)% 
0002 // Mass:  2.28646 ± 0.00014 GeV/c²
0003 // Shyam Kumar; INFN Bari, Italy
0004 // shyam.kumar@ba.infn.it; shyam055119@gmail.com 
0005 
0006 #ifdef __CINT__
0007 
0008 #pragma link off all globals;
0009 #pragma link off all classes;
0010 #pragma link off all functions;
0011 
0012 #pragma link C++ class PlotFile;
0013 #endif
0014 
0015 #ifndef __CINT__
0016 #include <stdio.h>
0017 #include <stdlib.h>
0018 #include <fstream>
0019 #include <iostream>
0020 #include <iomanip>
0021 #include <string>
0022 #include <sys/types.h>
0023 #include <sys/stat.h>
0024 #include <dirent.h>
0025 #include "math.h"
0026 #include "string.h"
0027 
0028 #include "TROOT.h"
0029 #include "TFile.h"
0030 #include "TChain.h"
0031 #include "TH1D.h"
0032 #include "TH2D.h"
0033 #include "TH3D.h"
0034 #include "THnSparse.h"
0035 #include "TStyle.h"
0036 #include "TCanvas.h"
0037 #include "TProfile.h"
0038 #include "TTree.h"
0039 #include "TNtuple.h"
0040 #include "TRandom3.h"
0041 #include "TMath.h"
0042 #include "TSystem.h"
0043 #include "TUnixSystem.h"
0044 #include "TVector2.h"
0045 #include "TVector3.h"
0046 #include "TLorentzVector.h"
0047 #include "TTreeReader.h"
0048 #include "TTreeReaderValue.h"
0049 #include "TTreeReaderArray.h"
0050 #include "TLatex.h"
0051 #endif
0052 #include "StPhysicalHelix.h"
0053 #include "SystemOfUnits.h"
0054 #include "PhysicalConstants.h"
0055 
0056 using namespace std;
0057 
0058 const double gPionMass = 0.13957;
0059 const double gKaonMass = 0.493677;
0060 const double gProtonMass = 0.938272;
0061 
0062 const double twoPi = 2.*3.1415927;
0063 const double eMass = 0.000511;
0064 const int pdg_Lcp = 4122, pdg_p = 2212, pdg_k = 321, pdg_pi = 211;
0065 bool debug = false;
0066 
0067 
0068 const double bField = -1.7; // Tesla
0069 
0070 TVector3 GetDCAToPrimaryVertex(const int index, TVector3 vtx);
0071 
0072 TLorentzVector GetDCADaughters(const int index1, const int index2, const int index3, TVector3 vtx,
0073   float *dcaDaughters, float &cosTheta, float &cosTheta_xy, float &decayLength, float &V0DcaToVtx, float &sigma_vtx);
0074 
0075 TTreeReaderArray<float> *rcMomPx2;
0076 TTreeReaderArray<float> *rcMomPy2;
0077 TTreeReaderArray<float> *rcMomPz2;
0078 TTreeReaderArray<float> *rcCharge2;
0079 
0080 TTreeReaderArray<float> *rcTrkLoca2;
0081 TTreeReaderArray<float> *rcTrkLocb2;
0082 TTreeReaderArray<float> *rcTrkTheta2;
0083 TTreeReaderArray<float> *rcTrkPhi2;
0084 
0085 int main(int argc, char **argv)
0086 {
0087   if(argc!=3 && argc!=1) return 0;
0088 
0089   TString listname;
0090   TString outname;
0091 
0092   if(argc==1)
0093   {
0094     listname  = "test.list";
0095     outname = "test.root";
0096   }
0097 
0098   if(argc==3)
0099   {
0100     listname = argv[1];
0101     outname = argv[2];
0102   }  
0103 
0104   TChain *chain = new TChain("events");
0105 
0106   int nfiles = 0;
0107   char filename[512];
0108   ifstream *inputstream = new ifstream;
0109   inputstream->open(listname.Data());
0110   if(!inputstream)
0111   {
0112     printf("[e] Cannot open file list: %s\n", listname.Data());
0113   }
0114   while(inputstream->good())
0115   {
0116     inputstream->getline(filename, 512);
0117     if(inputstream->good())
0118     {
0119      TFile *ftmp = TFile::Open(filename, "read");
0120      if(!ftmp||!(ftmp->IsOpen())||!(ftmp->GetNkeys())) 
0121      {
0122        printf("[e] Could you open file: %s\n", filename);
0123      } 
0124      else
0125      {
0126        cout<<"[i] Add "<<nfiles<<"th file: "<<filename<<endl;
0127        chain->Add(filename);
0128        nfiles++;
0129      }
0130    }
0131  }
0132  inputstream->close();
0133  printf("[i] Read in %d files with %lld events in total\n", nfiles, chain->GetEntries());
0134 
0135  TH1F *hEventStat = new TH1F("hEventStat", "Event statistics", 7, 0, 7);
0136  hEventStat->GetXaxis()->SetBinLabel(1, "MC events");
0137  hEventStat->GetXaxis()->SetBinLabel(2, "#Lambda_{c}^{+}");
0138  hEventStat->GetXaxis()->SetBinLabel(3, "#Lambda_{c}^{+} -> p + K^{-} + #pi^{+}");
0139  hEventStat->GetXaxis()->SetBinLabel(4, "Reco Signal #Lambda_{c}^{+}/#Lambda_{c}^{-}");
0140  hEventStat->GetXaxis()->SetBinLabel(5, "Reco Signal #Lambda_{c}^{+}");
0141  hEventStat->GetXaxis()->SetBinLabel(6, "Reco Signal #Lambda_{c}^{-}");
0142  hEventStat->GetXaxis()->SetBinLabel(7, "Reco Bkg #Lambda_{c}^{+}");
0143 
0144 
0145  TH1F *hMcMult = new TH1F("hMcMult", "MC multiplicity (|#eta| < 3.5);N_{MC}", 50, 0, 50);
0146 
0147  TH1F *hMcVtxX = new TH1F("hMcVtxX", "x position of MC vertex;x (mm)", 500, -5.0, 5.0);
0148  TH1F *hMcVtxY = new TH1F("hMcVtxY", "y position of MC vertex;y (mm)", 500, -5.0, 5.0);
0149  TH1F *hMcVtxZ = new TH1F("hMcVtxZ", "z position of MC vertex;z (mm)", 800, -200, 200);
0150 
0151  TH1F *hRecVtxX = new TH1F("hRecVtxX", "x position of Rec vertex;x (mm)", 500, -5.0, 5.0);
0152  TH1F *hRecVtxY = new TH1F("hRecVtxY", "y position of Rec vertex;y (mm)", 500, -5.0, 5.0);
0153  TH1F *hRecVtxZ = new TH1F("hRecVtxZ", "z position of Rec vertex;z (mm)", 800, -200, 200);
0154 
0155  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);
0156  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);
0157 
0158  TH2F *hMCLcpPtRap = new TH2F("hMCLcpPtRap", "MC #Lambda_{c}^{+};y;p_{T} (GeV/c)", 20, -5, 5, 100, 0, 10);
0159 
0160  TH2F *hMcPPtEta = new TH2F("hMcPPtEta", "MC P from #Lambda_{c}^{+} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0161  TH2F *hMcPPtEtaReco = new TH2F("hMcPPtEtaReco", "RC P from #Lambda_{c}^{+} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0162 
0163  TH2F *hMcKPtEta = new TH2F("hMcKPtEta", "MC K from #Lambda_{c}^{+} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0164  TH2F *hMcKPtEtaReco = new TH2F("hMcKPtEtaReco", "RC K from #Lambda_{c}^{+} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0165 
0166  TH2F *hMcPiPtEta = new TH2F("hMcPiPtEta", "MC #pi from #Lambda_{c}^{+} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0167  TH2F *hMcPiPtEtaReco = new TH2F("hMcPiPtEtaReco", "RC #pi from #Lambda_{c}^{+} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0168 
0169  TH1F *hNRecoVtx = new TH1F("hNRecoVtx", "Number of reconstructed vertices;N", 10, 0, 10);
0170 
0171  const char* part_name[3] = {"P", "K", "Pi"};
0172  const char* part_title[3] = {"P", "K", "#pi"};
0173  TH3F *hRcSecPartLocaToRCVtx[3];
0174  TH3F *hRcSecPartLocbToRCVtx[3];
0175  TH3F *hRcPrimPartLocaToRCVtx[3];
0176  TH3F *hRcPrimPartLocbToRCVtx[3];
0177  for(int i=0; i<3; i++)
0178  {
0179   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);
0180   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);
0181   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);
0182   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);
0183 }
0184 
0185 const char* axis_name[3] = {"x", "y", "z"};
0186 const int nDimDca = 4;
0187 const int nBinsDca[nDimDca] = {50, 20, 500, 50};
0188 const double minBinDca[nDimDca] = {0, -5, -1+0.002, 0};
0189 const double maxBinDca[nDimDca] = {5, 5, 1+0.002, 50};
0190 THnSparseF *hPrimTrkDcaToRCVtx[3][3];
0191 for(int i=0; i<3; i++)
0192 {
0193   for(int j=0; j<3; j++)
0194   {
0195    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);
0196  }
0197 }
0198 
0199 TH3F *h3PairDca12[2], *h3PairDca23[2], *h3PairDca13[2];
0200 TH3F *h3PairCosTheta[2];
0201 TH3F *h3PairDca[2];
0202 TH3F *h3PairDecayLength[2];
0203 const char* pair_name[2] = {"signal", "bkg"};
0204 const char* pair_title[2] = {"Signal", "Background"};
0205 for(int i=0; i<2; i++)
0206 {
0207   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);
0208   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);
0209   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);
0210   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);
0211 
0212   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);
0213 
0214   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);
0215 }
0216 
0217   // Invariant mass
0218 const char* cut_name[2] = {"all", "DCA"};
0219 TH3F *h3InvMass[2][2];
0220 for(int i=0; i<2; i++)
0221 {
0222   for(int j=0; j<2; j++)
0223   {
0224    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);
0225  }
0226 }
0227 
0228 TTreeReader treereader(chain);
0229   // MC  
0230 TTreeReaderArray<int> mcPartGenStatus = {treereader, "MCParticles.generatorStatus"};
0231 TTreeReaderArray<int> mcPartPdg = {treereader, "MCParticles.PDG"};
0232 TTreeReaderArray<float> mcPartCharge = {treereader, "MCParticles.charge"};
0233 TTreeReaderArray<unsigned int> mcPartParent_begin = {treereader, "MCParticles.parents_begin"};
0234 TTreeReaderArray<unsigned int> mcPartParent_end = {treereader, "MCParticles.parents_end"};
0235 TTreeReaderArray<int> mcPartParent_index = {treereader, "_MCParticles_parents.index"};
0236 TTreeReaderArray<unsigned int> mcPartDaughter_begin = {treereader, "MCParticles.daughters_begin"};
0237 TTreeReaderArray<unsigned int> mcPartDaughter_end = {treereader, "MCParticles.daughters_end"};
0238 TTreeReaderArray<int> mcPartDaughter_index = {treereader, "_MCParticles_daughters.index"};
0239 TTreeReaderArray<double> mcPartMass = {treereader, "MCParticles.mass"};
0240 TTreeReaderArray<double> mcPartVx = {treereader, "MCParticles.vertex.x"};
0241 TTreeReaderArray<double> mcPartVy = {treereader, "MCParticles.vertex.y"};
0242 TTreeReaderArray<double> mcPartVz = {treereader, "MCParticles.vertex.z"};
0243 TTreeReaderArray<float> mcMomPx = {treereader, "MCParticles.momentum.x"};
0244 TTreeReaderArray<float> mcMomPy = {treereader, "MCParticles.momentum.y"};
0245 TTreeReaderArray<float> mcMomPz = {treereader, "MCParticles.momentum.z"};
0246 TTreeReaderArray<double> mcEndPointX = {treereader, "MCParticles.endpoint.x"};
0247 TTreeReaderArray<double> mcEndPointY = {treereader, "MCParticles.endpoint.y"};
0248 TTreeReaderArray<double> mcEndPointZ = {treereader, "MCParticles.endpoint.z"};
0249 
0250 TTreeReaderArray<unsigned int> assocChSimID = {treereader, "ReconstructedChargedParticleAssociations.simID"};
0251 TTreeReaderArray<unsigned int> assocChRecID = {treereader, "ReconstructedChargedParticleAssociations.recID"};
0252 TTreeReaderArray<float> assocWeight = {treereader, "ReconstructedChargedParticleAssociations.weight"};
0253 
0254 TTreeReaderArray<float> rcMomPx = {treereader, "ReconstructedChargedParticles.momentum.x"};
0255 TTreeReaderArray<float> rcMomPy = {treereader, "ReconstructedChargedParticles.momentum.y"};
0256 TTreeReaderArray<float> rcMomPz = {treereader, "ReconstructedChargedParticles.momentum.z"};
0257 TTreeReaderArray<float> rcPosx = {treereader, "ReconstructedChargedParticles.referencePoint.x"};
0258 TTreeReaderArray<float> rcPosy = {treereader, "ReconstructedChargedParticles.referencePoint.y"};
0259 TTreeReaderArray<float> rcPosz = {treereader, "ReconstructedChargedParticles.referencePoint.z"};
0260 TTreeReaderArray<float> rcCharge = {treereader, "ReconstructedChargedParticles.charge"};
0261 TTreeReaderArray<int>   rcPdg = {treereader, "ReconstructedChargedParticles.PDG"};
0262 
0263 TTreeReaderArray<float> rcTrkLoca = {treereader, "CentralCKFTrackParameters.loc.a"};
0264 TTreeReaderArray<float> rcTrkLocb = {treereader, "CentralCKFTrackParameters.loc.b"};
0265 TTreeReaderArray<float> rcTrkqOverP = {treereader, "CentralCKFTrackParameters.qOverP"};
0266 TTreeReaderArray<float> rcTrkTheta = {treereader, "CentralCKFTrackParameters.theta"};
0267 TTreeReaderArray<float> rcTrkPhi = {treereader, "CentralCKFTrackParameters.phi"};
0268 
0269 rcMomPx2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.momentum.x"};
0270 rcMomPy2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.momentum.y"};
0271 rcMomPz2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.momentum.z"};
0272 rcCharge2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.charge"};
0273 
0274 rcTrkLoca2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.loc.a"};
0275 rcTrkLocb2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.loc.b"};
0276 rcTrkTheta2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.theta"};
0277 rcTrkPhi2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.phi"};
0278 
0279 TTreeReaderArray<float> CTVx = {treereader, "CentralTrackVertices.position.x"};
0280 TTreeReaderArray<float> CTVy = {treereader, "CentralTrackVertices.position.y"};
0281 TTreeReaderArray<float> CTVz = {treereader, "CentralTrackVertices.position.z"};
0282 TTreeReaderArray<int> CTVndf = {treereader, "CentralTrackVertices.ndf"};
0283 TTreeReaderArray<float> CTVchi2 = {treereader, "CentralTrackVertices.chi2"};
0284 TTreeReaderArray<float> CTVerr_xx = {treereader, "CentralTrackVertices.positionError.xx"};
0285 TTreeReaderArray<float> CTVerr_yy = {treereader, "CentralTrackVertices.positionError.yy"};
0286 TTreeReaderArray<float> CTVerr_zz = {treereader, "CentralTrackVertices.positionError.zz"};
0287 
0288 TTreeReaderArray<int> prim_vtx_index = {treereader, "PrimaryVertices_objIdx.index"};
0289 
0290 TTreeReaderArray<unsigned int> vtxAssocPart_begin = {treereader, "CentralTrackVertices.associatedParticles_begin"};
0291 TTreeReaderArray<unsigned int> vtxAssocPart_end = {treereader, "CentralTrackVertices.associatedParticles_end"};
0292 TTreeReaderArray<int> vtxAssocPart_index = {treereader, "_CentralTrackVertices_associatedParticles.index"};
0293 
0294   // File to create efficiency of D0 meson
0295 TFile *file_gen = new TFile("SignalLcpGen.root", "RECREATE");
0296 TTree *tree_gen = new TTree("treeMLSigGen", "treeMLSigGen"); 
0297 
0298    // Define variables to store in the Ntuple
0299 float pt_Lcp_gen, y_Lcp_gen;
0300 tree_gen->Branch("pt_Lcp_gen", &pt_Lcp_gen, "pt_Lcp_gen/F");  
0301 tree_gen->Branch("y_Lcp_gen", &y_Lcp_gen, "y_Lcp_gen/F"); 
0302 
0303 
0304     // Create a ROOT file to store the Ntuple
0305 TFile *file_signal = new TFile("SignalLcp.root", "RECREATE");
0306 TTree *tree_sig = new TTree("treeMLSig", "treeMLSig"); 
0307 
0308   // Define variables to store in the Ntuple
0309 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;
0310 float costheta_sig, costhetaxy_sig, pt_Lcp_sig, y_Lcp_sig, mass_Lcp_sig, sigma_vtx_sig, mult_sig;
0311 
0312     // Link the variables to the TTree branches
0313 tree_sig->Branch("d0_p", &d0_p_sig, "d0_p/F");
0314 tree_sig->Branch("d0_k", &d0_k_sig, "d0_k/F");
0315 tree_sig->Branch("d0_pi", &d0_pi_sig, "d0_pi/F"); 
0316 tree_sig->Branch("d0xy_p", &d0xy_p_sig, "d0xy_p/F");
0317 tree_sig->Branch("d0xy_k", &d0xy_k_sig, "d0xy_k/F");
0318 tree_sig->Branch("d0xy_pi", &d0xy_pi_sig, "d0xy_pi/F");  
0319 tree_sig->Branch("sum_d0xy", &sum_d0xy_sig, "sum_d0xy/F");        
0320 tree_sig->Branch("dca_12", &dca_12_sig, "dca_12/F");
0321 tree_sig->Branch("dca_Lcp", &dca_Lcp_sig, "dca_Lcp/F");
0322 tree_sig->Branch("pt_Lcp", &pt_Lcp_sig, "pt_Lcp/F");  
0323 tree_sig->Branch("y_Lcp", &y_Lcp_sig, "y_Lcp/F");  
0324 tree_sig->Branch("mass_Lcp", &mass_Lcp_sig, "mass_Lcp/F");              
0325 tree_sig->Branch("decay_length", &decay_length_sig, "decay_length/F");   
0326 tree_sig->Branch("costheta", &costheta_sig, "costheta/F"); 
0327 tree_sig->Branch("costheta_xy", &costhetaxy_sig, "costheta_xy/F"); 
0328 tree_sig->Branch("sigma_vtx", &sigma_vtx_sig, "sigma_vtx/F"); 
0329 tree_sig->Branch("mult", &mult_sig, "mult/F");             
0330 
0331 TFile *file_bkg = new TFile("BkgLcp.root", "RECREATE");
0332 TTree *tree_bkg = new TTree("treeMLBkg", "treeMLBkg"); 
0333 
0334   // Define variables to store in the Ntuple
0335 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;
0336 float costheta_bkg, costhetaxy_bkg, pt_Lcp_bkg, y_Lcp_bkg, mass_Lcp_bkg, sigma_vtx_bkg, mult_bkg;
0337   // Link the variables to the TTree branches
0338 tree_bkg->Branch("d0_p", &d0_p_bkg, "d0_p/F");
0339 tree_bkg->Branch("d0_k", &d0_k_bkg, "d0_k/F");
0340 tree_bkg->Branch("d0_pi", &d0_pi_bkg, "d0_pi/F");
0341 tree_bkg->Branch("d0xy_p", &d0xy_p_bkg, "d0xy_p/F");
0342 tree_bkg->Branch("d0xy_k", &d0xy_k_bkg, "d0xy_k/F");
0343 tree_bkg->Branch("d0xy_pi", &d0xy_pi_bkg, "d0xy_pi/F");  
0344 tree_bkg->Branch("sum_d0xy", &sum_d0xy_bkg, "sum_d0xy/F");            
0345 tree_bkg->Branch("dca_12", &dca_12_bkg, "dca_12/F");
0346 tree_bkg->Branch("dca_Lcp", &dca_Lcp_bkg, "dca_Lcp/F");
0347 tree_bkg->Branch("pt_Lcp", &pt_Lcp_bkg, "pt_Lcp/F");  
0348 tree_bkg->Branch("y_Lcp", &y_Lcp_bkg, "y_Lcp/F");  
0349 tree_bkg->Branch("mass_Lcp", &mass_Lcp_bkg, "mass_Lcp/F");              
0350 tree_bkg->Branch("decay_length", &decay_length_bkg, "decay_length/F");   
0351 tree_bkg->Branch("costheta", &costheta_bkg, "costheta/F");  
0352 tree_bkg->Branch("costheta_xy", &costhetaxy_bkg, "costheta_xy/F"); 
0353 tree_bkg->Branch("sigma_vtx", &sigma_vtx_bkg, "sigma_vtx/F"); 
0354 tree_bkg->Branch("mult", &mult_bkg, "mult/F");                    
0355 
0356 
0357 int nevents = 0;
0358 int count = 0;
0359 
0360 while(treereader.Next())
0361 {
0362   if(nevents%1000==0) printf("\nEvent No.-----> %d\n",nevents);
0363   // find MC primary vertex
0364   int nMCPart = mcPartMass.GetSize();
0365 
0366   TVector3 vertex_mc(-999., -999., -999.);
0367   for(int imc=0; imc<nMCPart; imc++)
0368   {
0369    // Status 4 describes the incoming electron/proton and the end point of electron or proton is MC Vertex  
0370    if(mcPartGenStatus[imc] == 4 && mcPartPdg[imc] == 11)
0371    {
0372      vertex_mc.SetXYZ(mcEndPointX[imc], mcEndPointY[imc], mcEndPointZ[imc]);
0373      break;
0374    }
0375  }
0376  hEventStat->Fill(0.5);
0377  hMcVtxX->Fill(vertex_mc.x());
0378  hMcVtxY->Fill(vertex_mc.y());
0379  hMcVtxZ->Fill(vertex_mc.z());
0380 
0381 
0382       // get RC primary vertex
0383  TVector3 vertex_rc(-999., -999., -999.);
0384  if(prim_vtx_index.GetSize()>0)
0385  {
0386    int rc_vtx_index = prim_vtx_index[0];
0387    vertex_rc.SetXYZ(CTVx[rc_vtx_index], CTVy[rc_vtx_index], CTVz[rc_vtx_index]);
0388  }
0389 
0390  hRecVtxX->Fill(vertex_rc.x());
0391  hRecVtxY->Fill(vertex_rc.y());
0392  hRecVtxZ->Fill(vertex_rc.z());
0393 
0394       // map MC and RC particles
0395  int nAssoc = assocChRecID.GetSize();
0396       map<int, int> assoc_map_mc_to_rc;  // MC --> RC Associations
0397       map<int, int> assoc_map_rc_to_mc;  // RC --> MC Associations
0398 
0399       for(unsigned int rc_index=0; rc_index<rcMomPx.GetSize(); rc_index++)
0400       {
0401       // loop over the association to find the matched MC particle
0402       // with largest weight
0403        double max_weight = 0;
0404        int matched_mc_index = -1;
0405        for(int j=0; j<nAssoc; j++)
0406        {
0407          if(assocChRecID[j] != rc_index) continue;
0408          if(assocWeight[j] > max_weight)
0409          {
0410           max_weight = assocWeight[j];
0411           matched_mc_index = assocChSimID[j];
0412         }
0413       }
0414 
0415       // build the map
0416       assoc_map_mc_to_rc[matched_mc_index] = rc_index;
0417       assoc_map_rc_to_mc[rc_index] = matched_mc_index;
0418     }
0419         // MC --> RC associations (kv; key,value for map)
0420    // std::cout << "MC to RC map:\n";
0421    // for (const auto& kv : assoc_map_mc_to_rc) {
0422    //     std::cout << "MC index " << kv.first << " ---> RC index " << kv.second << "\n";
0423    // }
0424 
0425     // RC --> MC associations
0426    // std::cout << "\nRC to MC map:\n";
0427    // for (const auto& kv : assoc_map_rc_to_mc) {
0428     //    std::cout << "RC index " << kv.first << " --->  MC index " << kv.second << "\n";
0429    // }
0430 
0431 
0432       // Loop over primary particles
0433     int nMcPart = 0;
0434     for(int imc=0; imc<nMCPart; imc++)
0435     {
0436      if(mcPartGenStatus[imc] == 1 && mcPartCharge[imc] != 0)
0437      {
0438        double dist = sqrt( pow(mcPartVx[imc]-vertex_mc.x(),2) + pow(mcPartVy[imc]-vertex_mc.y(),2) + pow(mcPartVz[imc]-vertex_mc.z(),2));      
0439        if(dist < 1e-4)
0440        {
0441           // count primary charged particles within |eta| < 3.5
0442         TVector3 mc_mom(mcMomPx[imc], mcMomPy[imc], mcMomPz[imc]);
0443         double mcEta = mc_mom.PseudoRapidity();
0444         if(fabs(mcEta) < 3.5) nMcPart++;
0445       }
0446     }
0447   }
0448   // Primary charged multiplicity at generated level with acceptance cut
0449   hMcMult->Fill(nMcPart);
0450 
0451 
0452   for(int imc=0; imc<nMCPart; imc++)
0453   {
0454    if(mcPartGenStatus[imc] == 1 && mcPartCharge[imc] != 0)
0455    {
0456      double dist = sqrt( pow(mcPartVx[imc]-vertex_mc.x(),2) + pow(mcPartVy[imc]-vertex_mc.y(),2) + pow(mcPartVz[imc]-vertex_mc.z(),2));      
0457      if(dist < 1e-4)
0458      {        
0459           // check if the MC particle is reconstructed
0460       int rc_index = -1;
0461       if(assoc_map_mc_to_rc.find(imc) != assoc_map_mc_to_rc.end()) rc_index = assoc_map_mc_to_rc[imc];
0462 
0463       if(rc_index>=0)
0464       {
0465         TVector3 dcaToVtx = GetDCAToPrimaryVertex(rc_index, vertex_rc);
0466 
0467         int ip = -1;
0468         if(fabs(mcPartPdg[imc]) == pdg_p) ip = 0;
0469         if(fabs(mcPartPdg[imc]) == pdg_k) ip = 1;
0470         if(fabs(mcPartPdg[imc]) == pdg_pi) ip = 2;
0471         if(ip>=0)
0472         {
0473          TVector3 mom(rcMomPx[rc_index], rcMomPy[rc_index], rcMomPz[rc_index]);
0474          if(ip<3)
0475          {
0476            hRcPrimPartLocaToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.Perp());
0477            hRcPrimPartLocbToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.z());
0478          }
0479 
0480          double fill1[] = {mom.Pt(), mom.Eta(), dcaToVtx.x(), nMcPart*1.};
0481          double fill2[] = {mom.Pt(), mom.Eta(), dcaToVtx.y(), nMcPart*1.};
0482          double fill3[] = {mom.Pt(), mom.Eta(), dcaToVtx.z(), nMcPart*1.};
0483          hPrimTrkDcaToRCVtx[ip][0]->Fill(fill1);
0484          hPrimTrkDcaToRCVtx[ip][1]->Fill(fill2);
0485          hPrimTrkDcaToRCVtx[ip][2]->Fill(fill3);
0486        }
0487      }
0488    }
0489  }
0490 }
0491 
0492  // Look for Λc⁺ → p K⁻ π⁺
0493 bool hasLc = false;
0494 vector<int> mc_index_Lcp_p;
0495 vector<int> mc_index_Lcp_k;
0496 vector<int> mc_index_Lcp_pi;
0497 mc_index_Lcp_p.clear();      
0498 mc_index_Lcp_k.clear();
0499 mc_index_Lcp_pi.clear();
0500 
0501 for(int imc=0; imc<nMCPart; imc++)
0502 {
0503  if(fabs(mcPartPdg[imc]) != pdg_Lcp) continue;
0504  hEventStat->Fill(1.5);
0505 
0506   // printf("MC Particle (Momentum) = (%f, %f, %f) \n",mcMomPx[imc], mcMomPy[imc], mcMomPz[imc]);           
0507  int nDuaghters = mcPartDaughter_end[imc]-mcPartDaughter_begin[imc];
0508  if(nDuaghters!=3) continue;
0509 
0510       // find Lc+ that decay into pKPi
0511  bool is_pkpi_decay = false;      
0512  int daug_index_1 = mcPartDaughter_index[mcPartDaughter_begin[imc]];
0513  int daug_index_2 = mcPartDaughter_index[mcPartDaughter_begin[imc]+1];
0514  int daug_index_3 = mcPartDaughter_index[mcPartDaughter_begin[imc]+2];    
0515  int daug_pdg_1 = mcPartPdg[daug_index_1];
0516  int daug_pdg_2 = mcPartPdg[daug_index_2];
0517  int daug_pdg_3 = mcPartPdg[daug_index_3];
0518 
0519  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)  || 
0520    (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) ||
0521    (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) )
0522  {
0523    is_pkpi_decay = true;
0524  }
0525  if(!is_pkpi_decay) continue;
0526    // Efficiency for Lc+
0527  TLorentzVector mc_part_gen;
0528  mc_part_gen.SetXYZM(mcMomPx[imc], mcMomPy[imc], mcMomPz[imc], mcPartMass[imc]);
0529  // Fill the momentum at Gen level Aug3, 2025
0530  float mcRap_gen = mc_part_gen.Rapidity();
0531  float mcPt_gen = mc_part_gen.Pt();  
0532  pt_Lcp_gen = mcPt_gen;
0533  y_Lcp_gen = mcRap_gen;
0534  tree_gen->Fill();
0535 
0536  hEventStat->Fill(2.5);
0537 
0538  if((fabs(daug_pdg_1)==pdg_p && fabs(daug_pdg_2)==pdg_k && fabs(daug_pdg_3)==pdg_pi))
0539  {
0540    mc_index_Lcp_p.push_back(daug_index_1);
0541    mc_index_Lcp_k.push_back(daug_index_2);
0542    mc_index_Lcp_pi.push_back(daug_index_3);
0543  }
0544  else if((fabs(daug_pdg_1)==pdg_p && fabs(daug_pdg_2)==pdg_pi && fabs(daug_pdg_3)==pdg_k))
0545  {
0546    mc_index_Lcp_p.push_back(daug_index_1);
0547    mc_index_Lcp_k.push_back(daug_index_3);
0548    mc_index_Lcp_pi.push_back(daug_index_2);
0549  }
0550  else if((fabs(daug_pdg_1)==pdg_k && fabs(daug_pdg_2)==pdg_p && fabs(daug_pdg_3)==pdg_pi))
0551  {
0552    mc_index_Lcp_p.push_back(daug_index_2);
0553    mc_index_Lcp_k.push_back(daug_index_1);
0554    mc_index_Lcp_pi.push_back(daug_index_3);
0555  }
0556  else if((fabs(daug_pdg_1)==pdg_k && fabs(daug_pdg_2)==pdg_pi && fabs(daug_pdg_3)==pdg_p))
0557  {
0558    mc_index_Lcp_p.push_back(daug_index_3);
0559    mc_index_Lcp_k.push_back(daug_index_1);
0560    mc_index_Lcp_pi.push_back(daug_index_2);
0561  }
0562  else if((fabs(daug_pdg_1)==pdg_pi && fabs(daug_pdg_2)==pdg_k && fabs(daug_pdg_3)==pdg_p))
0563  {
0564    mc_index_Lcp_p.push_back(daug_index_3);
0565    mc_index_Lcp_k.push_back(daug_index_2);
0566    mc_index_Lcp_pi.push_back(daug_index_1);
0567  }                              
0568    else 
0569    {
0570      mc_index_Lcp_p.push_back(daug_index_2);
0571      mc_index_Lcp_k.push_back(daug_index_3);
0572      mc_index_Lcp_pi.push_back(daug_index_1);
0573    }
0574    hasLc = true;
0575 
0576       // Lcp kinematics
0577    TLorentzVector mc_mom_vec;
0578    mc_mom_vec.SetXYZM(mcMomPx[imc], mcMomPy[imc], mcMomPz[imc], mcPartMass[imc]);
0579 //printf("Mass  = %f \n",mcPartMass[imc]);
0580 
0581    double mcRap = mc_mom_vec.Rapidity();
0582    double mcPt = mc_mom_vec.Pt();
0583    hMCLcpPtRap->Fill(mcRap, mcPt); 
0584 
0585       // decay daughter kinematics
0586    for(int ip = 0; ip<3; ip++)
0587    {
0588      int mc_part_index;
0589      if(ip==0) mc_part_index = mc_index_Lcp_p[mc_index_Lcp_p.size()-1];
0590      if(ip==1) mc_part_index = mc_index_Lcp_k[mc_index_Lcp_k.size()-1];
0591      if(ip==2) mc_part_index = mc_index_Lcp_pi[mc_index_Lcp_pi.size()-1];
0592 
0593      TLorentzVector mc_part_vec;
0594      mc_part_vec.SetXYZM(mcMomPx[mc_part_index], mcMomPy[mc_part_index], mcMomPz[mc_part_index], mcPartMass[mc_part_index]);
0595      if(ip==0) hMcPPtEta->Fill(mc_part_vec.Eta(), mc_part_vec.Pt());
0596      if(ip==1) hMcKPtEta->Fill(mc_part_vec.Eta(), mc_part_vec.Pt());
0597      if(ip==2) hMcPiPtEta->Fill(mc_part_vec.Eta(), mc_part_vec.Pt());         
0598 
0599      int rc_part_index = -1;
0600      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];
0601      if (debug) printf("Rec: ProngNo., MC Index, Reco Index = (%d, %d, %d) \n",ip,mc_part_index,rc_part_index);
0602      if(rc_part_index>=0)
0603      {
0604       TVector3 dcaToVtx = GetDCAToPrimaryVertex(rc_part_index, vertex_rc);
0605       TVector3 mom(rcMomPx[rc_part_index], rcMomPy[rc_part_index], rcMomPz[rc_part_index]);
0606       hRcSecPartLocaToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.Pt());
0607       hRcSecPartLocbToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.z());
0608 
0609           //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]);
0610     }
0611   }
0612 }
0613 
0614       // Get reconstructed protons, kaons and pions
0615       hNRecoVtx->Fill(CTVx.GetSize());
0616       const int pid_mode = 0; // 0 - truth; 1 - realistic
0617       vector<unsigned int> p_index;
0618       vector<unsigned int> k_index;      
0619       vector<unsigned int> pi_index;
0620       p_index.clear();
0621       k_index.clear();      
0622       pi_index.clear();
0623 
0624       for(unsigned int rc_index=0; rc_index<rcMomPx.GetSize(); rc_index++)
0625       {   
0626        if(pid_mode==0)
0627        {
0628          int iSimPartID = -1;
0629          if(assoc_map_rc_to_mc.find(rc_index) != assoc_map_rc_to_mc.end()) iSimPartID = assoc_map_rc_to_mc[rc_index];
0630          if(iSimPartID>=0)
0631          {
0632           if(fabs(mcPartPdg[iSimPartID]) == pdg_p) p_index.push_back(rc_index); 
0633           if(fabs(mcPartPdg[iSimPartID]) == pdg_k) k_index.push_back(rc_index);
0634           if(fabs(mcPartPdg[iSimPartID]) == pdg_pi) pi_index.push_back(rc_index);
0635         }
0636       }
0637       else if(pid_mode==1)
0638       {
0639        if(fabs(rcPdg[rc_index]) == pdg_p) p_index.push_back(rc_index);
0640        if(fabs(rcPdg[rc_index]) == pdg_k) k_index.push_back(rc_index);
0641        if(fabs(rcPdg[rc_index]) == pdg_pi) pi_index.push_back(rc_index);          
0642      }
0643    }
0644 
0645   // printf("Proton, Kaon, and Pion vector sizes: = (%d, %d, %d) \n",p_index.size(), k_index.size(),pi_index.size()); 
0646       // Combinatorics proton, kaon, and pion 
0647    for(unsigned int i=0; i<p_index.size(); i++)
0648    {
0649      TVector3 dcaToVtx_p = GetDCAToPrimaryVertex(p_index[i], vertex_rc);
0650      int q_proton = rcCharge[p_index[i]];
0651 
0652      for(unsigned int j=0; j<k_index.size(); j++)
0653      {
0654       TVector3 dcaToVtx_k = GetDCAToPrimaryVertex(k_index[j], vertex_rc);
0655       int q_kaon = rcCharge[k_index[j]];
0656 
0657       for(unsigned int k=0; k<pi_index.size(); k++)
0658       {
0659        TVector3 dcaToVtx_pi = GetDCAToPrimaryVertex(pi_index[k], vertex_rc);
0660        int q_pion = rcCharge[pi_index[k]];
0661 
0662        if ((q_proton == +1 && q_kaon == -1 && q_pion == +1) || (q_proton == -1 && q_kaon == +1 && q_pion == -1))       
0663        {
0664         bool is_Lcp_pkpi = false;
0665 
0666         int  mc_index_p = -1, mc_index_k = -1, mc_index_pi = -1;
0667         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]];
0668         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]];
0669         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]];
0670 
0671         for(unsigned int idaugh=0; idaugh<mc_index_Lcp_pi.size(); idaugh++)
0672         {
0673           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])
0674           {
0675            is_Lcp_pkpi = true;
0676            break;
0677          }
0678        }
0679 
0680        float dcaDaughters[3], cosTheta, decayLength, V0DcaToVtx, cosTheta_xy, sigma_vtx;
0681        TLorentzVector parent = GetDCADaughters(p_index[i], k_index[j], pi_index[k], vertex_rc, dcaDaughters, cosTheta, cosTheta_xy, decayLength, V0DcaToVtx, sigma_vtx);
0682 
0683 
0684        if(is_Lcp_pkpi)
0685        {
0686         hEventStat->Fill(3.5);
0687         if (q_proton == 1 && q_kaon == -1 && q_pion == 1)
0688         hEventStat->Fill(4.5);   // Λc⁺
0689       else if (q_proton == -1 && q_kaon == 1 && q_pion == -1)
0690         hEventStat->Fill(5.5);   // Λc⁻
0691 
0692       h3PairDca12[0]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters[0]);
0693       h3PairDca23[0]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters[1]); 
0694       h3PairDca13[0]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters[2]);              
0695       h3PairCosTheta[0]->Fill(parent.Pt(), parent.Rapidity(), cosTheta);
0696       h3PairDca[0]->Fill(parent.Pt(), parent.Rapidity(), V0DcaToVtx);
0697       h3PairDecayLength[0]->Fill(parent.Pt(), parent.Rapidity(), decayLength);
0698       h3InvMass[0][0]->Fill(parent.Pt(), parent.Rapidity(), parent.M());
0699 
0700       // Toplogical Variables for Signal
0701       d0_p_sig = dcaToVtx_p.Mag();
0702       d0_k_sig = dcaToVtx_k.Mag();
0703       d0_pi_sig = dcaToVtx_pi.Mag();              
0704       d0xy_p_sig = dcaToVtx_p.Perp();
0705       d0xy_k_sig = dcaToVtx_k.Perp();
0706       d0xy_pi_sig = dcaToVtx_pi.Perp();           
0707       sum_d0xy_sig = sqrt(d0xy_p_sig*d0xy_p_sig+d0xy_k_sig*d0xy_k_sig+d0xy_pi_sig*d0xy_pi_sig);           
0708       dca_12_sig = *min_element(dcaDaughters, dcaDaughters + 3);
0709       dca_Lcp_sig = V0DcaToVtx;
0710       decay_length_sig = decayLength;
0711       costheta_sig = cosTheta;
0712       costhetaxy_sig = cosTheta_xy;
0713       pt_Lcp_sig = parent.Pt();
0714       y_Lcp_sig = parent.Rapidity();
0715       mass_Lcp_sig = parent.M();
0716       sigma_vtx_sig = sigma_vtx;
0717       mult_sig = nMcPart;
0718       tree_sig->Fill();  
0719 
0720     }
0721     else
0722     {
0723       hEventStat->Fill(6.5);
0724       h3PairDca12[1]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters[0]);
0725       h3PairDca23[1]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters[1]); 
0726       h3PairDca13[1]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters[2]);                  
0727       h3PairCosTheta[1]->Fill(parent.Pt(), parent.Rapidity(), cosTheta);
0728       h3PairDca[1]->Fill(parent.Pt(), parent.Rapidity(), V0DcaToVtx);
0729       h3PairDecayLength[1]->Fill(parent.Pt(), parent.Rapidity(), decayLength);
0730       h3InvMass[1][0]->Fill(parent.Pt(), parent.Rapidity(), parent.M());
0731 
0732       // Toplogical Variables for Bkg
0733       d0_p_bkg = dcaToVtx_p.Mag();
0734       d0_k_bkg = dcaToVtx_k.Mag();
0735       d0_pi_bkg = dcaToVtx_pi.Mag();              
0736       d0xy_p_bkg = dcaToVtx_p.Perp();
0737       d0xy_k_bkg = dcaToVtx_k.Perp();
0738       d0xy_pi_bkg = dcaToVtx_pi.Perp();           
0739       sum_d0xy_bkg = sqrt(d0xy_p_bkg*d0xy_p_bkg+d0xy_k_bkg*d0xy_k_bkg+d0xy_pi_bkg*d0xy_pi_bkg);       
0740       dca_12_bkg = *min_element(dcaDaughters, dcaDaughters + 3);
0741       dca_Lcp_bkg = V0DcaToVtx;
0742       decay_length_bkg = decayLength;
0743       costheta_bkg = cosTheta;
0744       costhetaxy_bkg = cosTheta_xy;
0745       pt_Lcp_bkg = parent.Pt();
0746       y_Lcp_bkg = parent.Rapidity();
0747       mass_Lcp_bkg = parent.M();
0748       sigma_vtx_bkg = sigma_vtx;
0749       mult_bkg = nMcPart;
0750       tree_bkg->Fill();
0751     } 
0752 
0753   }
0754 }
0755 }
0756 }
0757 
0758 nevents++;
0759 }
0760 
0761 
0762 file_gen->cd();  
0763 tree_gen->Write();
0764 file_gen->Close();
0765 
0766 file_signal->cd();  
0767 tree_sig->Write();
0768 file_signal->Close();  
0769 
0770 file_bkg->cd();  
0771 tree_bkg->Write();
0772 file_bkg->Close(); 
0773 
0774 TFile *outfile = new TFile(outname.Data(), "recreate");
0775 hEventStat->SetMarkerSize(2);
0776 hEventStat->Write("");
0777 hMcMult->Write();
0778 hMcVtxX->Write();
0779 hMcVtxY->Write();
0780 hMcVtxZ->Write();
0781 hRecVtxX->Write();
0782 hRecVtxY->Write();
0783 hRecVtxZ->Write();
0784 
0785 hLcpDecayVxVy->Write();
0786 hLcpDecayVrVz->Write();
0787 
0788 hMCLcpPtRap->Write();
0789 hMcPPtEta->Write();
0790 hMcPPtEtaReco->Write();
0791 hMcPiPtEta->Write();
0792 hMcPiPtEtaReco->Write();
0793 hMcKPtEta->Write();
0794 hMcKPtEtaReco->Write();
0795 
0796 hNRecoVtx->Write();
0797 
0798 for(int ip=0; ip<3; ip++)
0799 {
0800   hRcSecPartLocaToRCVtx[ip]->Write();
0801   hRcSecPartLocbToRCVtx[ip]->Write();
0802   hRcPrimPartLocaToRCVtx[ip]->Write();
0803   hRcPrimPartLocbToRCVtx[ip]->Write();
0804 }
0805 
0806 for(int i=0; i<3; i++)
0807 {
0808   for(int j=0; j<3; j++)
0809   {
0810    hPrimTrkDcaToRCVtx[i][j]->Write();
0811  }
0812 }
0813 
0814 for(int i=0; i<2; i++)
0815 {
0816   h3PairDca12[i]->Write();
0817   h3PairDca23[i]->Write();
0818   h3PairDca13[i]->Write();            
0819   h3PairCosTheta[i]->Write();
0820   h3PairDca[i]->Write();
0821   h3PairDecayLength[i]->Write();
0822 }
0823 
0824 for(int i=0; i<2; i++)
0825 {
0826   for(int j=0; j<2; j++)
0827   {
0828    h3InvMass[i][j]->Write();
0829  }
0830 }
0831 
0832 
0833 outfile->Close();
0834 
0835 }
0836 
0837 //======================================
0838 TVector3 GetDCAToPrimaryVertex(const int index, TVector3 vtx)
0839 {
0840   //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));
0841 
0842   TVector3 pos(rcTrkLoca2->At(index) * sin(rcTrkPhi2->At(index)) * -1 * millimeter, rcTrkLoca2->At(index) * cos(rcTrkPhi2->At(index)) * millimeter, rcTrkLocb2->At(index) * millimeter);
0843   TVector3 mom(rcMomPx2->At(index), rcMomPy2->At(index), rcMomPz2->At(index));
0844 
0845   StPhysicalHelix pHelix(mom, pos, bField * tesla, rcCharge2->At(index));
0846 
0847   TVector3 vtx_tmp;
0848   vtx_tmp.SetXYZ(vtx.x()*millimeter, vtx.y()*millimeter, vtx.z()*millimeter);
0849   
0850   pHelix.moveOrigin(pHelix.pathLength(vtx_tmp));
0851   TVector3 dcaToVtx = pHelix.origin() - vtx_tmp;
0852 
0853   dcaToVtx.SetXYZ(dcaToVtx.x()/millimeter, dcaToVtx.y()/millimeter, dcaToVtx.z()/millimeter);
0854   
0855   return dcaToVtx;
0856 }
0857 
0858 //======================================
0859 // Local to global conversion 
0860 // x = - l0 Sin (phi), y = l0 Cos (phi), z = l1 
0861 TLorentzVector GetDCADaughters(const int index1, const int index2, const int index3, TVector3 vtx,
0862   float *dcaDaughters, float &cosTheta, float &cosTheta_xy, float &decayLength, float &V0DcaToVtx, float &sigma_vtx)
0863 {
0864   // -- get helix
0865   TVector3 pos1(rcTrkLoca2->At(index1) * sin(rcTrkPhi2->At(index1)) * -1 * millimeter, rcTrkLoca2->At(index1) * cos(rcTrkPhi2->At(index1)) * millimeter, rcTrkLocb2->At(index1) * millimeter);
0866   TVector3 pos2(rcTrkLoca2->At(index2) * sin(rcTrkPhi2->At(index2)) * -1 * millimeter, rcTrkLoca2->At(index2) * cos(rcTrkPhi2->At(index2)) * millimeter, rcTrkLocb2->At(index2) * millimeter);
0867   TVector3 pos3(rcTrkLoca2->At(index3) * sin(rcTrkPhi2->At(index3)) * -1 * millimeter, rcTrkLoca2->At(index3) * cos(rcTrkPhi2->At(index3)) * millimeter, rcTrkLocb2->At(index3) * millimeter);
0868   
0869   TVector3 mom1(rcMomPx2->At(index1), rcMomPy2->At(index1), rcMomPz2->At(index1));
0870   TVector3 mom2(rcMomPx2->At(index2), rcMomPy2->At(index2), rcMomPz2->At(index2));
0871   TVector3 mom3(rcMomPx2->At(index3), rcMomPy2->At(index3), rcMomPz2->At(index3));
0872 
0873   float charge1 = rcCharge2->At(index1);
0874   float charge2 = rcCharge2->At(index2);
0875   float charge3 = rcCharge2->At(index3);  
0876   
0877   StPhysicalHelix p1Helix(mom1, pos1, bField * tesla, charge1);
0878   StPhysicalHelix p2Helix(mom2, pos2, bField * tesla, charge2);
0879   StPhysicalHelix p3Helix(mom3, pos3, bField * tesla, charge3);
0880 
0881   TVector3 vtx_tmp;
0882   vtx_tmp.SetXYZ(vtx.x()*millimeter, vtx.y()*millimeter, vtx.z()*millimeter);
0883 
0884   // -- Full calculation for pathlengths
0885   // Point1
0886   pair<double, double> const ss_12 = p1Helix.pathLengths(p2Helix); 
0887   pair<double, double> const ss_13 = p1Helix.pathLengths(p3Helix);
0888   // Point2 
0889   pair<double, double> const ss_23 = p2Helix.pathLengths(p3Helix); 
0890   pair<double, double> const ss_21 = p2Helix.pathLengths(p1Helix); 
0891   // Point3
0892   pair<double, double> const ss_31 = p3Helix.pathLengths(p1Helix); 
0893   pair<double, double> const ss_32 = p3Helix.pathLengths(p2Helix);  
0894   
0895   // Point1 w.r.t. helices 2 and 3, later average will be taken 
0896   TVector3 const P1P2 = p1Helix.at(ss_12.first); TVector3 const P1P3 = p1Helix.at(ss_13.first); 
0897   // Point2 w.r.t. helices 1 and 3, later average will be taken 
0898   TVector3 const P2P1 = p2Helix.at(ss_21.first); TVector3 const P2P3 = p2Helix.at(ss_23.first); 
0899   // Point3 w.r.t. helices 1 and 2, later average will be taken 
0900   TVector3 const P3P1 = p3Helix.at(ss_31.first); TVector3 const P3P2 = p3Helix.at(ss_32.first); 
0901 
0902  // printf("P1P2 = (%2.4f, %2.4f, %2.4f) \t P1P3 = (%2.4f, %2.4f, %2.4f)\n", P1P2.x(), P1P2.y(), P1P2.z(), P1P3.x(), P1P3.y(), P1P3.z());
0903  // printf("P2P1 = (%2.4f, %2.4f, %2.4f) \t P2P3 = (%2.4f, %2.4f, %2.4f)\n", P2P1.x(), P2P1.y(), P2P1.z(), P2P3.x(), P2P3.y(), P2P3.z());
0904  // printf("P3P1 = (%2.4f, %2.4f, %2.4f) \t P3P2 = (%2.4f, %2.4f, %2.4f)\n", P3P1.x(), P3P1.y(), P3P1.z(), P3P2.x(), P3P2.y(), P3P2.z());
0905   // Evaluate the DCA Points
0906   // Between track1 and track2
0907   TVector3 const p1 = 0.5*(P1P2+P1P3);
0908   TVector3 const p2 = 0.5*(P2P1+P2P3);
0909   TVector3 const p3 = 0.5*(P3P1+P3P2);
0910 
0911 
0912   // cout << ss.first << "  " << ss.second << endl;
0913   // printf("p1AtDcaToP2 origin = (%2.4f, %2.4f, %2.4f)\n", p1AtDcaToP2.x(), p1AtDcaToP2.y(), p1AtDcaToP2.z());
0914   // printf("p2AtDcaToP1 origin = (%2.4f, %2.4f, %2.4f)\n", p2AtDcaToP1.x(), p2AtDcaToP1.y(), p2AtDcaToP1.z());
0915   
0916   // -- calculate DCA between daughters 
0917   // Three dcaDaughters; dcaDaughters_12, dcaDaughters_23, dcaDaughters_13
0918   dcaDaughters[0] = (p1 - p2).Mag()/millimeter;
0919   dcaDaughters[1] = (p2 - p3).Mag()/millimeter;
0920   dcaDaughters[2] = (p3 - p1).Mag()/millimeter; 
0921 
0922   // -- calculate Momentum of each daughters at the DCA point
0923   TVector3 const p1MomAtDcap2 = p1Helix.momentumAt(ss_12.first,  bField * tesla);
0924   TVector3 const p1MomAtDcap3 = p1Helix.momentumAt(ss_13.first,  bField * tesla);
0925   TVector3 const p2MomAtDca3 = p2Helix.momentumAt(ss_23.first, bField * tesla);
0926   TVector3 const p2MomAtDca1 = p2Helix.momentumAt(ss_21.first, bField * tesla);
0927   TVector3 const p3MomAtDcap1 = p3Helix.momentumAt(ss_31.first, bField * tesla); 
0928   TVector3 const p3MomAtDcap2 = p3Helix.momentumAt(ss_32.first, bField * tesla); 
0929 
0930   // printf("p1MomAtDcap2 = (%2.4f, %2.4f, %2.4f) \t p1MomAtDcap3 = (%2.4f, %2.4f, %2.4f)\n", p1MomAtDcap2.x(), p1MomAtDcap2.y(), p1MomAtDcap2.z(), p1MomAtDcap3.x(), p1MomAtDcap3.y(), p1MomAtDcap3.z());
0931   // printf("p2MomAtDca3 = (%2.4f, %2.4f, %2.4f) \t p2MomAtDca1 = (%2.4f, %2.4f, %2.4f)\n", p2MomAtDca3.x(), p2MomAtDca3.y(), p2MomAtDca3.z(), p2MomAtDca1.x(), p2MomAtDca1.y(), p2MomAtDca1.z());
0932   // printf("p3MomAtDcap1 = (%2.4f, %2.4f, %2.4f) \t p3MomAtDcap2 = (%2.4f, %2.4f, %2.4f)\n", p3MomAtDcap1.x(), p3MomAtDcap1.y(), p3MomAtDcap1.z(), p3MomAtDcap2.x(), p3MomAtDcap2.y(), p3MomAtDcap2.z());
0933 
0934   
0935   TVector3 const p1MomAtDca = 0.5*(p1MomAtDcap2+p1MomAtDcap3);  
0936   TVector3 const p2MomAtDca = 0.5*(p2MomAtDca3+p2MomAtDca1); 
0937   TVector3 const p3MomAtDca = 0.5*(p3MomAtDcap1+p3MomAtDcap2); 
0938   
0939   TLorentzVector p1FourMom(p1MomAtDca, sqrt(p1MomAtDca.Mag2()+gProtonMass*gProtonMass));
0940   TLorentzVector p2FourMom(p2MomAtDca, sqrt(p2MomAtDca.Mag2()+gKaonMass*gKaonMass));
0941   TLorentzVector p3FourMom(p3MomAtDca, sqrt(p3MomAtDca.Mag2()+gPionMass*gPionMass));  
0942   
0943   TLorentzVector parent = p1FourMom + p2FourMom + p3FourMom;
0944 
0945   // -- calculate cosThetaStar
0946   // TLorentzVector pairFourMomReverse(-parent.Px(), -parent.Py(), -parent.Pz(), parent.E());
0947   // TLorentzVector p1FourMomStar = p1FourMom;
0948   // p1FourMomStar.Boost(pairFourMomReverse.Vect());
0949   // float cosThetaStar = std::cos(p1FourMomStar.Vect().Angle(parent.Vect()));
0950 
0951   //Approach1: Assign the decay vertex the centroid coordinate of the triangle
0952    TVector3 decayVertex = 1./3*(p1 + p2+ p3);
0953  
0954   //Approach2: Assign the decay vertex as the circumcenter of the triangle (p1, p2, p3)
0955   // Compute side vectors of triangle; a = p1, b = p2, c = p3
0956     // Side vectors
0957   TVector3 ab = p2 - p1; 
0958   TVector3 ac = p3 - p1;
0959 
0960     // Cross product of ab and ac
0961   TVector3 abXac = ab.Cross(ac);
0962 
0963   double ab2 = ab.Mag2();
0964   double ac2 = ac.Mag2();
0965   double abXac2 = abXac.Mag2();
0966 
0967     // Vector from point a to the circhumspere center
0968   TVector3 a_to_CircumsphereCenter = (abXac.Cross(ab) * ac2 + ac.Cross(abXac) * ab2) * (1.0 / (2.0 * abXac2));
0969 
0970     // Magnitude of the above vector is radius
0971   double radius = a_to_CircumsphereCenter.Mag();
0972 
0973     // Coordinate of the Circumsphere center is the decayvertex in 3D (A point at equal distance of three PCA)
0974  
0975   //TVector3 decayVertex = p1 + a_to_CircumsphereCenter; // uncomment it if you want to define as circumsphere center
0976 
0977   // printf("dist_p1_to_decayvertex = (%1.4f) \t dist_p2_to_decayvertex = (%1.4f) \t dist_p3_to_decayvertex = (%1.4f) \n", (p1-decayVertex).Mag(), (p2-decayVertex).Mag(), (p3-decayVertex).Mag());
0978 
0979     // Confirm by evaluating the distance with vertices
0980 
0981   sigma_vtx = sqrt((p1-decayVertex).Mag2()+(p2-decayVertex).Mag2()+(p3-decayVertex).Mag2())/millimeter;
0982 
0983   // -- calculate pointing angle and decay length with respect to primary vertex
0984   //    if decay vertex is a tertiary vertex
0985   //    -> only rough estimate -> needs to be updated after secondary vertex is found
0986   TVector3 vtxToV0 = decayVertex - vtx_tmp;
0987   TVector3 vtxToV0_xy(vtxToV0.x(), vtxToV0.y(), 0.);
0988   TVector3 parent_xy(parent.Vect().x(),parent.Vect().y(),0.);
0989   float pointingAngle = vtxToV0.Angle(parent.Vect());
0990   float pointingAngle_xy = vtxToV0_xy.Angle(parent_xy);
0991   cosTheta = std::cos(pointingAngle);
0992   cosTheta_xy = std::cos(pointingAngle_xy);  
0993   decayLength = vtxToV0.Mag()/millimeter;
0994 
0995   // -- calculate V0 DCA to primary vertex
0996   V0DcaToVtx = decayLength * std::sin(pointingAngle);
0997 
0998   //TVector3 dcaToVtx = GetDCAToPrimaryVertex(parent.Vect(), decayVertex, 0, vtx);
0999   //V0DcaToVtx = dcaToVtx.Mag();
1000 
1001   return parent;
1002 }
1003
[ 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!