Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-09-01 09:34:09 

0001 // D0-tagged Jets
0002 
0003 #ifdef __CINT__
0004 
0005 #pragma link off all globals;
0006 #pragma link off all classes;
0007 #pragma link off all functions;
0008 
0009 #pragma link C++ class PlotFile;
0010 #endif
0011 
0012 #ifndef __CINT__
0013 #include <stdio.h>
0014 #include <stdlib.h>
0015 #include <fstream>
0016 #include <iostream>
0017 #include <iomanip>
0018 #include <string>
0019 #include <sys/types.h>
0020 #include <sys/stat.h>
0021 #include <dirent.h>
0022 #include "math.h"
0023 #include "string.h"
0024 
0025 #include "TROOT.h"
0026 #include "TFile.h"
0027 #include "TChain.h"
0028 #include "TH1D.h"
0029 #include "TH2D.h"
0030 #include "TH3D.h"
0031 #include "THnSparse.h"
0032 #include "TStyle.h"
0033 #include "TCanvas.h"
0034 #include "TProfile.h"
0035 #include "TTree.h"
0036 #include "TNtuple.h"
0037 #include "TRandom3.h"
0038 #include "TMath.h"
0039 #include "TSystem.h"
0040 #include "TUnixSystem.h"
0041 #include "TVector2.h"
0042 #include "TVector3.h"
0043 #include "TLorentzVector.h"
0044 #include "TTreeReader.h"
0045 #include "TTreeReaderValue.h"
0046 #include "TTreeReaderArray.h"
0047 #include "TLatex.h"
0048 #include "TMinuit.h"
0049 #include "Math/Functor.h"
0050 #include "Fit/Fitter.h"
0051 #include "Math/Minimizer.h"
0052 #endif
0053 
0054 
0055 #include "StPhysicalHelix.h"
0056 #include "SystemOfUnits.h"
0057 #include "PhysicalConstants.h"
0058 
0059 #include "fastjet/ClusterSequence.hh"
0060 #include "fastjet/ClusterSequenceArea.hh"
0061 #include "fastjet/GhostedAreaSpec.hh"
0062 #include "fastjet/PseudoJet.hh"
0063 #include "fastjet/AreaDefinition.hh"
0064 
0065 using namespace std;
0066 
0067 StPhysicalHelix* gHelix1 = nullptr;
0068 StPhysicalHelix* gHelix2 = nullptr;
0069 
0070 const double gPionMass = 0.13957;
0071 const double gKaonMass = 0.493677;
0072 
0073 const double twoPi = 2.*3.1415927;
0074 const double eMass = 0.000511;
0075 
0076 const double bField = -1.7; // Tesla
0077 
0078 //JetClustering Parameters:
0079 double minCstPt            = 0.2 ;               // minimum pT of objects
0080 double maxCstPt            = 100.;               // maximum pT of objects
0081 double minJetPt            = 1.0 ;               // minimum jet pT
0082 double ghostMaxRap         = 3.5;                // maximum rapidity of ghosts
0083 double ghostArea           = 0.01;               // area per ghost
0084 int numGhostRepeat         = 1;                              // reuse count of ghosts
0085 int removeelectrons        = 1;
0086 int nhitcut                = 0;
0087 
0088 TVector3 getDcaToVtx(const int index, TVector3 vtx);
0089 void fcnVertexFit(int& npar, double* grad, double& fval, double* par, int iflag);
0090 void getDecayVertex_Chi2fit(const int index1, const int index2, double &s1, double &s2, TVector3 &vertex, double &chi2_ndf, double * parFitErr);
0091 
0092 TLorentzVector getPairParent(const int index1, const int index2, TVector3 vtx,
0093                  float &dcaDaughters, float &cosTheta, float &cosTheta_xy, float &decayLength, float &V0DcaToVtx, float &sigma_vtx, TVector3 &decayVertex,TVector3 &decayVertex_ana, double &chi2_ndf,  double * parFitErr);
0094 
0095 TTreeReaderArray<float> *rcMomPx2;
0096 TTreeReaderArray<float> *rcMomPy2;
0097 TTreeReaderArray<float> *rcMomPz2;
0098 TTreeReaderArray<float> *rcCharge2;
0099 
0100 TTreeReaderArray<float> *rcTrkLoca2;
0101 TTreeReaderArray<float> *rcTrkLocb2;
0102 TTreeReaderArray<float> *rcTrkTheta2;
0103 TTreeReaderArray<float> *rcTrkPhi2;
0104 TTreeReaderArray<std::array<float, 21>> *rcTrkCov;
0105 
0106 
0107 
0108 // Define function for Chi2 minimization between two helices    
0109 struct Chi2Minimization {    
0110      StPhysicalHelix fhelix1, fhelix2;
0111      std::array<float, 21> fcov1, fcov2; // full covariance matrix
0112      
0113      Chi2Minimization(StPhysicalHelix helix1, StPhysicalHelix helix2, std::array<float, 21> cov1, std::array<float, 21> cov2) : fhelix1(helix1),fhelix2(helix2), fcov1(cov1), fcov2(cov2) {}
0114     // Implementation of the function to be minimized
0115     double operator() (const double *par) {
0116     double x = par[0];
0117     double y = par[1];
0118     double z = par[2];
0119     double s1 = par[3];
0120     double s2 = par[4];   
0121     double f = 0;
0122     TVector3 vertex(x, y, z);
0123     TVector3 p1 = fhelix1.at(s1);
0124     TVector3 p2 = fhelix2.at(s2);
0125     TVector3 mom1 = fhelix1.momentumAt(s1,  bField * tesla);
0126     TVector3 mom2 = fhelix2.momentumAt(s2,  bField * tesla);
0127     
0128    // x= −l0 sinϕ , y=l0 cosϕ , z=l
0129    // Recalculate l0 at PCA for error propagation
0130      float l0_track1 = p1.Pt(); float l1_track1 = p1.Z(); double phi_track1 = mom1.Phi();
0131      float l0_track2 = p2.Pt(); float l1_track2 = p2.Z(); double phi_track2 = mom2.Phi();
0132      // Track1: σx^2​=sin^2ϕ⋅σℓ0​^2​+ℓ0^2​cos^2ϕ⋅σϕ^2​+2⋅ℓ0​sinϕcosϕ⋅Cov(ℓ0​,ϕ)
0133      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]; 
0134      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];
0135   
0136      // σy^2​=cos^2ϕ⋅σℓ0^​2​+ℓ0^2​sin^2ϕ⋅σϕ^2​−2⋅ℓ0​sinϕcosϕ⋅Cov(ℓ0​,ϕ)
0137      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]; 
0138      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]; 
0139      
0140      // σz^2​
0141      float sigz1_2 = fcov1[2];
0142      float sigz2_2 = fcov2[2];  
0143      // convert to mm
0144     double d1_x = 10.*(vertex - p1).X(); double d2_x = 10.*(vertex - p2).X(); 
0145     double d1_y = 10.*(vertex - p1).Y(); double d2_y = 10.*(vertex - p2).Y();
0146     double d1_z = 10.*(vertex - p1).Z(); double d2_z = 10.*(vertex - p2).Z();
0147     
0148      f = d1_x*d1_x/sigx1_2 + d2_x*d2_x/sigx2_2 + d1_y*d1_y/sigy1_2 + d2_y*d2_y/sigy2_2 + d1_z*d1_z/sigz1_2+ d2_z*d2_z/sigz2_2; // chi2
0149       return f;
0150     }
0151     };
0152 
0153 int main(int argc, char **argv)
0154 {
0155     
0156   TString listname;
0157   TString outname;
0158   TString collname;
0159   float R_value;
0160   TString signal, bkg, gen;
0161   TString signalmc; 
0162 
0163   if(argc==1)
0164     {
0165       listname  = "test.list";
0166       outname = "test.root";
0167       collname = "ep";   
0168       R_value = 1.0;
0169       signal = "signal.root";
0170       bkg = "bkg.root";
0171       gen = "gen.root";
0172       signalmc = "signalmc.root";
0173     }
0174 
0175   if(argc>=8)
0176     {
0177       listname = argv[1];
0178       outname = argv[2];
0179       collname = argv[3];      
0180       R_value = std::atof(argv[4]);       
0181       signal = argv[5]; 
0182       bkg = argv[6];
0183       gen = argv[7];
0184       signalmc = "signalmc.root";
0185     }
0186 
0187   TChain *chain = new TChain("events");
0188 
0189   int nfiles = 0;
0190   char filename[512];
0191   ifstream *inputstream = new ifstream;
0192   inputstream->open(listname.Data());
0193   if (!inputstream->is_open())
0194   {
0195   printf("[e] Cannot open file list: %s\n", listname.Data());
0196   return 0; // or handle as needed
0197   }
0198 
0199   while (inputstream->good())
0200   {
0201   inputstream->getline(filename, 512);
0202   if (inputstream->good())
0203   {
0204     TFile *ftmp = TFile::Open(filename, "READ");
0205     if (!ftmp || !ftmp->IsOpen() || !ftmp->GetNkeys())
0206     {
0207       printf("[e] Skipping bad file: %s\n", filename);
0208       if (ftmp) { ftmp->Close(); delete ftmp; }
0209       continue; 
0210     }
0211     cout << "[i] Add " << nfiles << "th file: " << filename << endl;
0212     chain->Add(filename);
0213     nfiles++;
0214 
0215     ftmp->Close(); // cleanup
0216     delete ftmp;
0217   }
0218 }
0219 
0220   inputstream->close();
0221 
0222   printf("[i] Read in %d files with %lld events in total\n", nfiles, chain->GetEntries());
0223 
0224   while (inputstream->good())
0225   {
0226   inputstream->getline(filename, 512);
0227   if (inputstream->good())
0228   {
0229     TFile *ftmp = TFile::Open(filename, "READ");
0230     if (!ftmp || !ftmp->IsOpen() || !ftmp->GetNkeys())
0231     {
0232       printf("[e] Skipping bad file: %s\n", filename);
0233       if (ftmp) { ftmp->Close(); delete ftmp; }
0234       continue; 
0235     }
0236     cout << "[i] Add " << nfiles << "th file: " << filename << endl;
0237     chain->Add(filename);
0238     nfiles++;
0239 
0240     ftmp->Close(); // cleanup
0241     delete ftmp;
0242   }
0243 }
0244 
0245   inputstream->close();
0246 
0247   printf("[i] Read in %d files with %lld events in total\n", nfiles, chain->GetEntries());
0248 
0249   TH1F *hEventStat = new TH1F("hEventStat", "Event statistics", 7, 0, 7);
0250   hEventStat->GetXaxis()->SetBinLabel(1, "MC events");
0251   hEventStat->GetXaxis()->SetBinLabel(2, "D0");
0252   hEventStat->GetXaxis()->SetBinLabel(3, "D0 -> pi+K");
0253   hEventStat->GetXaxis()->SetBinLabel(4, "Reco D0");
0254   hEventStat->GetXaxis()->SetBinLabel(5, "Reco Signal D0");
0255   hEventStat->GetXaxis()->SetBinLabel(6, "Reco Signal D0bar");
0256   hEventStat->GetXaxis()->SetBinLabel(7, "Reco Bkg D0");
0257 
0258   TH1F *hMcMult = new TH1F("hMcMult", "MC multiplicity (|#eta| < 3.5);N_{MC}", 50, 0, 50);
0259 
0260   // Secondary vertex with chi2 fit method
0261   TH3F *hRes_SVx_Helixfit = new TH3F("hRes_SVx_Helixfit", "Fit method: Residual of SV (X); p_{T} (GeV/c); y ; SVx_{rec}-SVx_{mc} (mm)", 500, 0.0, 50.0, 80, -4.0, 4.0, 1000, -5.0, 5.0); 
0262   TH3F *hRes_SVy_Helixfit = new TH3F("hRes_SVy_Helixfit", "Fit method: Residual of SV (Y); p_{T} (GeV/c); y ; SVy_{rec}-SVy_{mc} (mm)", 500, 0.0, 50.0, 80, -4.0, 4.0, 1000, -5.0, 5.0); 
0263   TH3F *hRes_SVz_Helixfit = new TH3F("hRes_SVz_Helixfit", "Fit method: Residual of SV (Z); p_{T} (GeV/c); y ; SVz_{rec}-SVz_{mc} (mm)", 500, 0.0, 50.0, 80, -4.0, 4.0, 1000, -5.0, 5.0);   
0264  
0265   // Secondary vertex with analytical method
0266   TH3F *hRes_SVx_Helixana = new TH3F("hRes_SVx_Helixana", "Analytical method: Residual of SV (X); p_{T} (GeV/c); y ; SVx_{rec}-SVx_{mc} (mm)", 500, 0.0, 50.0, 80, -4.0, 4.0, 1000, -5.0, 5.0); 
0267   TH3F *hRes_SVy_Helixana = new TH3F("hRes_SVy_Helixana", "Analytical method: Residual of SV (Y); p_{T} (GeV/c); y ; SVy_{rec}-SVy_{mc} (mm)", 500, 0.0, 50.0, 80, -4.0, 4.0, 1000, -5.0, 5.0); 
0268   TH3F *hRes_SVz_Helixana = new TH3F("hRes_SVz_Helixana", "Analytical method: Residual of SV (Z); p_{T} (GeV/c); y ; SVz_{rec}-SVz_{mc} (mm)", 500, 0.0, 50.0, 80, -4.0, 4.0, 1000, -5.0, 5.0);   
0269  
0270   // Secondary vertex  (XY) fit method
0271   TH3F *hRes_SVxy_Helixfit = new TH3F("hRes_SVxy_Helixfit", "Fit method: Residual of SV (XY); p_{T} (GeV/c); y ; SVxy_{rec}-SVxy_{mc} (mm)", 500, 0.0, 50.0, 80, -4.0, 4.0, 1000, -5.0, 5.0); 
0272 
0273   // Secondary vertex  (XY) analytical method
0274   TH3F *hRes_SVxy_Helixana = new TH3F("hRes_SVxy_Helixana", "Analytical method: Residual of SV (XY); p_{T} (GeV/c); y ; SVxy_{rec}-SVxy_{mc} (mm)", 500, 0.0, 50.0, 80, -4.0, 4.0, 1000, -5.0, 5.0); 
0275     
0276   TH3F *hRes_SVx_Helixfit_pull = new TH3F("hRes_SVx_Helixfit_pull", "Fit method: Pull of SV (X); p_{T} (GeV/c); y ; SVx_{rec}-SVx_{mc} (mm)", 500, 0.0, 50.0, 80, -4.0, 4.0, 1000, -5.0, 5.0); 
0277   TH3F *hRes_SVy_Helixfit_pull = new TH3F("hRes_SVy_Helixfit_pull", "Fit method: Pull of SV (Y); p_{T} (GeV/c); y ; SVy_{rec}-SVy_{mc} (mm)", 500, 0.0, 50.0, 80, -4.0, 4.0, 1000, -5.0, 5.0); 
0278   TH3F *hRes_SVz_Helixfit_pull = new TH3F("hRes_SVz_Helixfit_pull", "Fit method: Pull of SV (Z); p_{T} (GeV/c); y ; SVz_{rec}-SVz_{mc} (mm)", 500, 0.0, 50.0, 80, -4.0, 4.0, 1000, -5.0, 5.0); 
0279   
0280   TH1F *hchi2_vtx = new TH1F("hchi2_vtx", "Helix Calculation: Chi2/ndf; #chi^{2}/ndf; Entries (a.u.)", 1000, 0.0, 50.0); 
0281   TH1F *hchi2_vtx_sig = new TH1F("hchi2_vtx_sig", "Helix Calculation: Chi2/ndf; #chi^{2}/ndf; Entries (a.u.)", 1000, 0.0, 50.0); 
0282   TH1F *hchi2_vtx_bkg = new TH1F("hchi2_vtx_bkg", "Helix Calculation: Chi2/ndf; #chi^{2}/ndf; Entries (a.u.)", 1000, 0.0, 50.0);   
0283 
0284   TH1F *hMcVtxX = new TH1F("hMcVtxX", "x position of MC vertex;x (mm)", 100, -5.05, 4.95);
0285   TH1F *hMcVtxY = new TH1F("hMcVtxY", "y position of MC vertex;y (mm)", 500, -5.01, 4.99);
0286   TH1F *hMcVtxZ = new TH1F("hMcVtxZ", "z position of MC vertex;z (mm)", 400, -200, 200);
0287   
0288   TH1F *hPullVtxX = new TH1F("hPullVtxX", "Pull x position of MC vertex;(Vx_{rec}-Vx_{mc})/#sigma_{vx}", 100, -5.05, 4.95);
0289   TH1F *hPullVtxY = new TH1F("hPullVtxY", "Pull y position of MC vertex;(Vy_{rec}-Vy_{mc})/#sigma_{vy}", 500, -5.01, 4.99);
0290   TH1F *hPullVtxZ = new TH1F("hPullVtxZ", "Pull z position of MC vertex;(Vz_{rec}-Vz_{mc})/#sigma_{vz}", 400, -200, 200);
0291 
0292   TH2F *hD0DecayVxVy = new TH2F("hD0DecayVxVy", "D^{0} 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);
0293   TH2F *hD0DecayVrVz = new TH2F("hD0DecayVrVz", "D^{0} decay vertex to primary vertex;#Deltav_{z} (mm);#Deltav_{r} (mm)", 100, -2, 2, 100, -0.2, 1.8);
0294 
0295   TH2F *hMCD0PtRap = new TH2F("hMCD0PtRap", "MC D^{0};y;p_{T} (GeV/c)", 20, -5, 5, 100, 0, 10);
0296 
0297   TH2F *hMcPiPtEta = new TH2F("hMcPiPtEta", "MC #pi from D^{0} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0298   TH2F *hMcPiPtEtaReco = new TH2F("hMcPiPtEtaReco", "RC #pi from D^{0} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0299 
0300   TH2F *hMcKPtEta = new TH2F("hMcKPtEta", "MC K from D^{0} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0301   TH2F *hMcKPtEtaReco = new TH2F("hMcKPtEtaReco", "RC K from D^{0} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0302 
0303   TH1F *hNRecoVtx = new TH1F("hNRecoVtx", "Number of reconstructed vertices;N", 10, 0, 10);
0304 
0305   const char* part_name[3] = {"Pi", "K", "P"};
0306   const char* part_title[3] = {"#pi", "K", "P"};
0307   TH3F *hRcSecPartLocaToRCVtx[2];
0308   TH3F *hRcSecPartLocbToRCVtx[2];
0309   TH3F *hRcPrimPartLocaToRCVtx[2];
0310   TH3F *hRcPrimPartLocbToRCVtx[2];
0311   for(int i=0; i<2; i++)
0312     {
0313       hRcSecPartLocaToRCVtx[i] = new TH3F(Form("hRcSec%sLocaToRCVtx",part_name[i]), Form( "DCA_{xy} distribution for D^{0} decayed %s;p_{T} (GeV/c);#eta;DCA_{xy} (mm)", part_title[i]), 100, 0, 10, 20, -5, 5, 100, 0, 1);
0314       hRcSecPartLocbToRCVtx[i] = new TH3F(Form("hRcSec%sLocbToRCVtx",part_name[i]), Form( "DCA_{z} distribution for D^{0} decayed %s;p_{T} (GeV/c);#eta;DCA_{z} (mm)", part_title[i]), 100, 0, 10, 20, -5, 5, 100, -0.5, 0.5);
0315       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);
0316       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);
0317     }
0318 
0319   const char* axis_name[3] = {"x", "y", "z"};
0320   const int nDimDca = 4;
0321   const int nBinsDca[nDimDca] = {50, 20, 500, 50};
0322   const double minBinDca[nDimDca] = {0, -5, -1+0.002, 0};
0323   const double maxBinDca[nDimDca] = {5, 5, 1+0.002, 50};
0324   THnSparseF *hPrimTrkDcaToRCVtx[3][3];
0325   for(int i=0; i<3; i++)
0326     {
0327       for(int j=0; j<3; j++)
0328     {
0329       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);
0330     }
0331     }
0332 
0333   TH3F *h3PairDca12[2];
0334   TH3F *h3PairCosTheta[2];
0335   TH3F *h3PairDca[2];
0336   TH3F *h3PairDecayLength[2];
0337   const char* pair_name[2] = {"signal", "bkg"};
0338   const char* pair_title[2] = {"Signal", "Background"};
0339   for(int i=0; i<2; i++)
0340     {
0341       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);
0342 
0343       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);
0344 
0345       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);
0346 
0347       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);
0348     }
0349 
0350   // Invariant mass
0351   const char* cut_name[2] = {"all", "DCA"};
0352   TH3F *h3InvMass[2][2];
0353   for(int i=0; i<2; i++)
0354     {
0355       for(int j=0; j<2; j++)
0356     {
0357       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, 1.6, 2.0);
0358     }
0359     }
0360   
0361   // Fragmentaion variable  
0362   TH3F *h3sig_z = new TH3F("h3sig_z","Signal_D0;Z;y;M_D0(GeV/c^{2})",1000, -10, 10, 20, -5, 5, 100, 1.6, 2.0);    
0363   TH3F *h3bkg_z = new TH3F("h3bkg_z","Bkg_D0;Z;y;M_D0(GeV/c^{2})",1000, -10, 10, 20, -5, 5, 100, 1.6, 2.0);
0364   TH1F *hreco_eta = new TH1F("hreco_eta", "Eta distribution of reco particles ;eta", 500, -5,+5);
0365 
0366   TH1F *hmc_eta   = new TH1F("hmc_eta",   "MC Particle #eta;#eta;Number of Particles", 500, -5, 5);
0367   TH1F *hmc_eta_e = new TH1F("hmc_eta_e", "Electron #eta;#eta;Number of Particles",       500, -5, 5);
0368   TH1F *hmc_eta_pi = new TH1F("hmc_eta_pi", "Pion #eta;#eta;Number of Particles",         500, -5, 5);
0369   TH1F *hmc_eta_k  = new TH1F("hmc_eta_k",  "Kaon #eta;#eta;Number of Particles",         500, -5, 5);
0370   TH1F *hmc_eta_p  = new TH1F("hmc_eta_p",  "Proton #eta;#eta;Number of Particles",       500, -5, 5);
0371 
0372   TH1F *hreco_eta_e  = new TH1F("hreco_eta_e",  "Electron #eta;#eta;Number of Particles", 500, -5, 5);
0373   TH1F *hreco_eta_pi = new TH1F("hreco_eta_pi", "Pion #eta;#eta;Number of Particles",     500, -5, 5);
0374   TH1F *hreco_eta_k  = new TH1F("hreco_eta_k",  "Kaon #eta;#eta;Number of Particles",     500, -5, 5);
0375   TH1F *hreco_eta_p  = new TH1F("hreco_eta_p",  "Proton #eta;#eta;Number of Particles",   500, -5, 5);
0376   
0377   
0378   TTreeReader treereader(chain);
0379   // MC Particles 
0380   TTreeReaderArray<int> mcPartGenStatus = {treereader, "MCParticles.generatorStatus"};
0381   TTreeReaderArray<int> mcPartPdg = {treereader, "MCParticles.PDG"};
0382   TTreeReaderArray<float> mcPartCharge = {treereader, "MCParticles.charge"};
0383   TTreeReaderArray<unsigned int> mcPartParent_begin = {treereader, "MCParticles.parents_begin"};
0384   TTreeReaderArray<unsigned int> mcPartParent_end = {treereader, "MCParticles.parents_end"};
0385   TTreeReaderArray<int> mcPartParent_index = {treereader, "_MCParticles_parents.index"};
0386   TTreeReaderArray<unsigned int> mcPartDaughter_begin = {treereader, "MCParticles.daughters_begin"};
0387   TTreeReaderArray<unsigned int> mcPartDaughter_end = {treereader, "MCParticles.daughters_end"};
0388   TTreeReaderArray<int> mcPartDaughter_index = {treereader, "_MCParticles_daughters.index"};
0389   TTreeReaderArray<double> mcPartMass = {treereader, "MCParticles.mass"};
0390   TTreeReaderArray<double> mcPartVx = {treereader, "MCParticles.vertex.x"};
0391   TTreeReaderArray<double> mcPartVy = {treereader, "MCParticles.vertex.y"};
0392   TTreeReaderArray<double> mcPartVz = {treereader, "MCParticles.vertex.z"};
0393   TTreeReaderArray<double> mcMomPx = {treereader, "MCParticles.momentum.x"};
0394   TTreeReaderArray<double> mcMomPy = {treereader, "MCParticles.momentum.y"};
0395   TTreeReaderArray<double> mcMomPz = {treereader, "MCParticles.momentum.z"};
0396   TTreeReaderArray<double> mcEndPointX = {treereader, "MCParticles.endpoint.x"};
0397   TTreeReaderArray<double> mcEndPointY = {treereader, "MCParticles.endpoint.y"};
0398   TTreeReaderArray<double> mcEndPointZ = {treereader, "MCParticles.endpoint.z"};
0399   // Reconstructed Particles
0400   TTreeReaderArray<int> assocChSimID = {treereader, "_ReconstructedChargedParticleAssociations_sim.index"};
0401   TTreeReaderArray<int> assocChRecID = {treereader, "_ReconstructedChargedParticleAssociations_rec.index"};
0402   TTreeReaderArray<float> assocWeight = {treereader, "ReconstructedChargedParticleAssociations.weight"};
0403  
0404   TTreeReaderArray<float> rcMomPx = {treereader, "ReconstructedChargedParticles.momentum.x"};
0405   TTreeReaderArray<float> rcMomPy = {treereader, "ReconstructedChargedParticles.momentum.y"};
0406   TTreeReaderArray<float> rcMomPz = {treereader, "ReconstructedChargedParticles.momentum.z"};
0407   TTreeReaderArray<float> rcPosx = {treereader, "ReconstructedChargedParticles.referencePoint.x"};
0408   TTreeReaderArray<float> rcPosy = {treereader, "ReconstructedChargedParticles.referencePoint.y"};
0409   TTreeReaderArray<float> rcPosz = {treereader, "ReconstructedChargedParticles.referencePoint.z"};
0410   TTreeReaderArray<float> rcCharge = {treereader, "ReconstructedChargedParticles.charge"};
0411   TTreeReaderArray<int>   rcPdg = {treereader, "ReconstructedChargedParticles.PDG"};
0412   TTreeReaderArray<float> TrkRecoE = {treereader, "ReconstructedChargedParticles.energy"};  
0413   TTreeReaderArray<float> TrkRecoM = {treereader, "ReconstructedChargedParticles.mass"};
0414   
0415   TTreeReaderArray<float> rcTrkLoca = {treereader, "CentralCKFTrackParameters.loc.a"};
0416   TTreeReaderArray<float> rcTrkLocb = {treereader, "CentralCKFTrackParameters.loc.b"};
0417   TTreeReaderArray<float> rcTrkqOverP = {treereader, "CentralCKFTrackParameters.qOverP"};
0418   TTreeReaderArray<float> rcTrkTheta = {treereader, "CentralCKFTrackParameters.theta"};
0419   TTreeReaderArray<float> rcTrkPhi = {treereader, "CentralCKFTrackParameters.phi"};
0420 
0421   rcMomPx2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.momentum.x"};
0422   rcMomPy2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.momentum.y"};
0423   rcMomPz2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.momentum.z"};
0424   rcCharge2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.charge"};
0425  
0426 
0427   rcTrkLoca2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.loc.a"};
0428   rcTrkLocb2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.loc.b"};
0429   rcTrkTheta2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.theta"};
0430   rcTrkPhi2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.phi"};
0431   rcTrkCov = new TTreeReaderArray<std::array<float, 21>>{treereader, "CentralCKFTrackParameters.covariance.covariance[21]"};
0432 
0433   TTreeReaderArray<float> CTVx = {treereader, "CentralTrackVertices.position.x"};
0434   TTreeReaderArray<float> CTVy = {treereader, "CentralTrackVertices.position.y"};
0435   TTreeReaderArray<float> CTVz = {treereader, "CentralTrackVertices.position.z"};
0436   TTreeReaderArray<int> CTVndf = {treereader, "CentralTrackVertices.ndf"};
0437   TTreeReaderArray<float> CTVchi2 = {treereader, "CentralTrackVertices.chi2"};
0438   TTreeReaderArray<float> CTVerr_xx = {treereader, "CentralTrackVertices.positionError.xx"};
0439   TTreeReaderArray<float> CTVerr_yy = {treereader, "CentralTrackVertices.positionError.yy"};
0440   TTreeReaderArray<float> CTVerr_zz = {treereader, "CentralTrackVertices.positionError.zz"};
0441   TTreeReaderArray<float> incQ2 = {treereader, "InclusiveKinematicsTruth.Q2"};
0442   TTreeReaderArray<float> incxB = {treereader, "InclusiveKinematicsTruth.x"};
0443 
0444   TTreeReaderArray<int> prim_vtx_index = {treereader, "PrimaryVertices_objIdx.index"};
0445 
0446   TTreeReaderArray<unsigned int> vtxAssocPart_begin = {treereader, "CentralTrackVertices.associatedParticles_begin"};
0447   TTreeReaderArray<unsigned int> vtxAssocPart_end = {treereader, "CentralTrackVertices.associatedParticles_end"};
0448   TTreeReaderArray<int> vtxAssocPart_index = {treereader, "_CentralTrackVertices_associatedParticles.index"};
0449 
0450 
0451   TTreeReaderArray<float> EvtQ2 ={treereader, "InclusiveKinematicsElectron.Q2"};
0452   TTreeReaderArray<float> Evtx =  {treereader, "InclusiveKinematicsElectron.x"};
0453   TTreeReaderArray<float> EvtQ2Gen = {treereader, "InclusiveKinematicsTruth.Q2"};
0454   TTreeReaderArray<float> EvtxGen = {treereader, "InclusiveKinematicsTruth.x"};
0455 
0456   TTreeReaderArray<int>  ScatElecRecoId = {treereader, "_InclusiveKinematicsElectron_scat.index"};
0457   TTreeReaderArray<int> ScatElecGenId = {treereader, "MCScatteredElectrons_objIdx.index"};
0458    
0459   TTreeReaderArray<int> TrkPartAssocRec = {treereader, "_ReconstructedChargedParticleAssociations_rec.index"};
0460   TTreeReaderArray<int> TrkPartAssocSim = {treereader, "_ReconstructedChargedParticleAssociations_sim.index"};
0461   TTreeReaderArray<unsigned int> TrkRecoNhits = {treereader, "CentralCKFTrajectories.nMeasurements"};
0462 
0463   // Gen particles
0464   TTreeReaderArray<float>  TrkGenE = {treereader, "GeneratedParticles.energy"};
0465   TTreeReaderArray<float> TrkGenPx = {treereader, "GeneratedParticles.momentum.x"};
0466   TTreeReaderArray<float> TrkGenPy = {treereader, "GeneratedParticles.momentum.y"};
0467   TTreeReaderArray<float> TrkGenPz = {treereader, "GeneratedParticles.momentum.z"};
0468   TTreeReaderArray<float> TrkGenM = {treereader, "GeneratedParticles.mass"};
0469   TTreeReaderArray<int> TrkGenPDG = {treereader, "GeneratedParticles.PDG"};
0470   TTreeReaderArray<float> TrkGenCharge = {treereader, "GeneratedParticles.charge"};
0471    
0472 
0473   
0474     // Create a ROOT file to store the Ntuple
0475    TFile *file_signal = new TFile(signal.Data(), "RECREATE");
0476    TTree *tree_sig = new TTree("treeMLSig", "treeMLSig"); 
0477 
0478   // Define variables to store in the Ntuple
0479   float d0_pi_sig, d0_k_sig, d0xy_pi_sig, d0xy_k_sig, sum_d0xy_sig, dca_12_sig, dca_D0_sig, decay_length_sig, xB_sig, Q2_sig;
0480   float costheta_sig, costhetaxy_sig, pt_D0_sig, y_D0_sig, mass_D0_sig, sigma_vtx_sig, mult_sig, signif_d0xy_pi_sig, signif_d0xy_k_sig, chi2_dca_sig,z_sig,dr_sig,angle_sig,etajet_sig,pTjet_sig, D0eta_sig, pt_parent_sig;
0481   
0482     // Link the variables to the TTree branches
0483   tree_sig->Branch("d0_pi", &d0_pi_sig, "d0_pi/F");
0484   tree_sig->Branch("d0_k", &d0_k_sig, "d0_k/F"); 
0485   tree_sig->Branch("d0xy_pi", &d0xy_pi_sig, "d0xy_pi/F");
0486   tree_sig->Branch("d0xy_k", &d0xy_k_sig, "d0xy_k/F");
0487   tree_sig->Branch("sum_d0xy", &sum_d0xy_sig, "sum_d0xy/F");        
0488   tree_sig->Branch("dca_12", &dca_12_sig, "dca_12/F");
0489   tree_sig->Branch("dca_D0", &dca_D0_sig, "dca_D0/F");
0490   tree_sig->Branch("pt_D0", &pt_D0_sig, "pt_D0/F");  
0491   tree_sig->Branch("y_D0", &y_D0_sig, "y_D0/F");  
0492   tree_sig->Branch("mass_D0", &mass_D0_sig, "mass_D0/F");              
0493   tree_sig->Branch("decay_length", &decay_length_sig, "decay_length/F");   
0494   tree_sig->Branch("costheta", &costheta_sig, "costheta/F"); 
0495   tree_sig->Branch("costheta_xy", &costhetaxy_sig, "costheta_xy/F"); 
0496   tree_sig->Branch("sigma_vtx", &sigma_vtx_sig, "sigma_vtx/F"); 
0497   tree_sig->Branch("mult", &mult_sig, "mult/F"); 
0498   tree_sig->Branch("signif_d0xy_pi", &signif_d0xy_pi_sig, "signif_d0xy_pi/F");               
0499   tree_sig->Branch("signif_d0xy_k", &signif_d0xy_k_sig, "signif_d0xy_k/F"); 
0500   tree_sig->Branch("chi2_dca", &chi2_dca_sig, "chi2_dca/F");
0501   tree_sig->Branch("xB", &xB_sig, "xB/F");
0502   tree_sig->Branch("Q2", &Q2_sig, "Q2/F");
0503   tree_sig->Branch("z", &z_sig, "z/F");
0504   tree_sig->Branch("dr", &dr_sig, "dr/F");
0505   tree_sig->Branch("angle", &angle_sig, "angle/F");
0506   tree_sig->Branch("etajet",&etajet_sig,"etajet/F");
0507   tree_sig->Branch("pTjet",&pTjet_sig,"pTjet/F");
0508   tree_sig->Branch("etaD0",&D0eta_sig,"etaD0/F");
0509   tree_sig->Branch("pt_parent",&pt_parent_sig,"pt_parent/F");
0510 
0511   
0512   TFile *file_bkg = new TFile(bkg.Data(), "RECREATE");
0513   TTree *tree_bkg = new TTree("treeMLBkg", "treeMLBkg"); 
0514   
0515   // Define variables to store in the Ntuple
0516   float d0_pi_bkg, d0_k_bkg, d0xy_pi_bkg, d0xy_k_bkg, sum_d0xy_bkg, dca_12_bkg, dca_D0_bkg, decay_length_bkg, xB_bkg, Q2_bkg; 
0517   float costheta_bkg, costhetaxy_bkg, pt_D0_bkg, y_D0_bkg, mass_D0_bkg, sigma_vtx_bkg, mult_bkg, signif_d0xy_pi_bkg, signif_d0xy_k_bkg, chi2_dca_bkg,z_bkg,dr_bkg,angle_bkg, etajet_bkg, pTjet_bkg, D0eta_bkg,pt_parent_bkg;
0518   
0519   // Link the variables to the TTree branches
0520   tree_bkg->Branch("d0_pi", &d0_pi_bkg, "d0_pi/F");
0521   tree_bkg->Branch("d0_k", &d0_k_bkg, "d0_k/F");
0522   tree_bkg->Branch("d0xy_pi", &d0xy_pi_bkg, "d0xy_pi/F");
0523   tree_bkg->Branch("d0xy_k", &d0xy_k_bkg, "d0xy_k/F");
0524   tree_bkg->Branch("sum_d0xy", &sum_d0xy_bkg, "sum_d0xy/F");            
0525   tree_bkg->Branch("dca_12", &dca_12_bkg, "dca_12/F");
0526   tree_bkg->Branch("dca_D0", &dca_D0_bkg, "dca_D0/F");
0527   tree_bkg->Branch("pt_D0", &pt_D0_bkg, "pt_D0/F");  
0528   tree_bkg->Branch("y_D0", &y_D0_bkg, "y_D0/F");  
0529   tree_bkg->Branch("mass_D0", &mass_D0_bkg, "mass_D0/F");              
0530   tree_bkg->Branch("decay_length", &decay_length_bkg, "decay_length/F");   
0531   tree_bkg->Branch("costheta", &costheta_bkg, "costheta/F");  
0532   tree_bkg->Branch("costheta_xy", &costhetaxy_bkg, "costheta_xy/F"); 
0533   tree_bkg->Branch("sigma_vtx", &sigma_vtx_bkg, "sigma_vtx/F"); 
0534   tree_bkg->Branch("mult", &mult_bkg, "mult/F");
0535   tree_bkg->Branch("signif_d0xy_pi", &signif_d0xy_pi_bkg, "signif_d0xy_pi/F");  
0536   tree_bkg->Branch("signif_d0xy_k", &signif_d0xy_k_bkg, "signif_d0xy_k/F");  
0537   tree_bkg->Branch("chi2_dca", &chi2_dca_bkg, "chi2_dca/F");   
0538   tree_bkg->Branch("xB", &xB_bkg, "xB/F"); 
0539   tree_bkg->Branch("Q2", &Q2_bkg, "Q2/F"); 
0540   tree_bkg->Branch("z", &z_bkg, "z/F");
0541   tree_bkg->Branch("dr", &dr_bkg, "dr/F");
0542   tree_bkg->Branch("angle", &angle_bkg, "angle/F");
0543   tree_bkg->Branch("etajet",&etajet_bkg,"etajet/F");
0544   tree_bkg->Branch("pTjet",&pTjet_bkg,"pTjet/F");
0545   tree_bkg->Branch("etaD0",&D0eta_bkg,"etaD0/F");
0546   tree_bkg->Branch("pt_parent",&pt_parent_bkg,"pt_parent/F");
0547 
0548 
0549   // Generate the tree with true D0 Properties
0550   TFile *fout_mcgen = new TFile(signalmc.Data(),"RECREATE");
0551   TTree *tree_D0 = new TTree("D0Tree","D0 Meson Kinematics");
0552 
0553   float d0_px, d0_py, d0_pz;
0554   float d0_pt, d0_eta, d0_y, d0_mass;
0555 
0556   tree_D0->Branch("px",&d0_px,"px/F");
0557   tree_D0->Branch("py",&d0_py,"py/F");
0558   tree_D0->Branch("pz",&d0_pz,"pz/F");
0559   tree_D0->Branch("pt",&d0_pt,"pt/F");
0560   tree_D0->Branch("eta",&d0_eta,"eta/F");
0561   tree_D0->Branch("rapidity",&d0_y,"rapidity/F");
0562   tree_D0->Branch("mass",&d0_mass,"mass/F");
0563 
0564   //ROOT file for Gen Jets 
0565   TFile *file_gen = new TFile(gen.Data(), "RECREATE");
0566   TTree *tree_gen_sig = new TTree("GenTree_sig","GenTree_sig");//Gen Jet Tree sig
0567   TTree *tree_gen_bkg = new TTree("GenTree_bkg","GenTree_bkg");
0568   float pt_gjet_sig,eta_gjet_sig,y_gjet_sig,genZ_sig, gdr_eta_sig, gdr_y_sig,gD0jetAngle_sig;
0569   float pt_gjet_bkg,eta_gjet_bkg,y_gjet_bkg,genZ_bkg, gdr_eta_bkg, gdr_y_bkg,gD0jetAngle_bkg;
0570     
0571   tree_gen_sig->Branch("pt_gjet",&pt_gjet_sig,"pt_gjet/F");
0572   tree_gen_sig->Branch("eta_gjet",&eta_gjet_sig,"eta_gjet/F");
0573   tree_gen_sig->Branch("y_gjet",&y_gjet_sig,"y_gjet/F");
0574   tree_gen_sig->Branch("Z", &genZ_sig,"Z/F");
0575   tree_gen_sig->Branch("DeltaR_inEtaPhi", &gdr_eta_sig,"DeltaR_inEtaPhi/F");
0576   tree_gen_sig->Branch("DeltaR_inYPhi", &gdr_y_sig,"DeltaR_inYPhi/F");
0577   tree_gen_sig->Branch("Angle_bw_D0andJet", &gD0jetAngle_sig,"Angle_bw_D0andJet/F");
0578   
0579   tree_gen_bkg->Branch("pt_gjet",&pt_gjet_bkg,"pt_gjet/F");
0580   tree_gen_bkg->Branch("eta_gjet",&eta_gjet_bkg,"eta_gjet/F");
0581   tree_gen_bkg->Branch("y_gjet",&y_gjet_bkg,"y_gjet/F");
0582   tree_gen_bkg->Branch("Z", &genZ_bkg,"Z/F");
0583   tree_gen_bkg->Branch("DeltaR_inEtaPhi", &gdr_eta_bkg,"DeltaR_inEtaPhi/F");
0584   tree_gen_bkg->Branch("DeltaR_inYPhi", &gdr_y_bkg,"DeltaR_inYPhi/F");
0585   tree_gen_bkg->Branch("Angle_bw_D0andJet", &gD0jetAngle_bkg,"Angle_bw_D0andJet/F");
0586   
0587 //Variables for Jet Clustering  
0588       int NEVENTS = 0;
0589       int EVETMULTRECO = 0; 
0590       int EVETMULTGEN = 0;
0591       int ScatteredERecId = 0;
0592       int ScatteredEGenId = 0;
0593       float EventQ2 = -999;
0594       float Eventx = -999;
0595       float EventQ2Gen = -999;
0596       float EventxGen = -999;
0597 
0598 //vectors to store jet clustering outputs
0599 // Vertex
0600     std::vector<float> Vertex_x;
0601     std::vector<float> Vertex_y;
0602     std::vector<float> Vertex_z;
0603     std::vector<int>   Vertex_ndf;
0604     std::vector<float> Vertex_chi2;
0605     std::vector<float> VertexErr_xx;
0606     std::vector<float> VertexErr_yy;
0607     std::vector<float> VertexErr_zz;
0608 
0609     // Reco Jets (Variable-length vectors for multiple jets per event)
0610     std::vector<float> RecoJet_pt;
0611     std::vector<float> RecoJet_eta;
0612     std::vector<float> RecoJet_rapidity;
0613     std::vector<float> RecoJet_phi;
0614     std::vector<float> RecoJet_E;
0615     std::vector<float> RecoJet_M;
0616     std::vector<bool>  RecoJet_hasElectron;
0617     std::vector<float> RecoJet_maxPtPart_pt;
0618         std::vector<float> RecoJet_rapi;
0619     // Reco jet constituents
0620     std::vector<std::vector<float>> RecoJet_constituent_pt;
0621     std::vector<std::vector<float>> RecoJet_constituent_eta;
0622     std::vector<std::vector<float>> RecoJet_constituent_phi;
0623     std::vector<std::vector<int>> RecoJet_constituent_nhits;
0624     std::vector<std::vector<int>> RecoJet_constituent_pdgid;
0625     std::vector<std::vector<int>> RecoJet_constituent_pdgidTruth;
0626     std::vector<std::vector<int>> RecoJet_constituent_idx;
0627         std::vector<std::vector<float>> RecoJet_constituent_energy;
0628         std::vector<std::vector<float>> RecoJet_constituent_rapi;
0629     
0630     // Gen Jets (Variable-length vectors for multiple jets per event)
0631     std::vector<float> GenJet_pt;
0632     std::vector<float> GenJet_eta;
0633     std::vector<float> GenJet_phi;
0634     std::vector<float> GenJet_E;
0635     std::vector<float> GenJet_M;
0636     std::vector<bool> GenJet_hasElectron;
0637     std::vector<bool> GenJet_hasNeutral;
0638     std::vector<float> GenJet_maxPtPart_pt;
0639         std::vector<float> GenJet_rapi;
0640     
0641     // Gen jet constituents
0642     std::vector<std::vector<float>> GenJet_constituent_pt;
0643     std::vector<std::vector<float>> GenJet_constituent_eta;
0644     std::vector<std::vector<float>> GenJet_constituent_phi;
0645     std::vector<std::vector<int>> GenJet_constituent_pdgid;
0646     std::vector<std::vector<float>> GenJet_constituent_rapi;
0647     std::vector<std::vector<int>> GenJet_constituent_idx;
0648     std::vector<std::vector<float>> GenJet_constituent_energy;
0649     
0650   int nevents = 0;
0651   int mult_charged = 0;
0652   while(treereader.Next())
0653     {
0654       if(nevents%1000==0) printf("\n[i] New event %d\n",nevents);
0655       float Q2_mc = -1000.;
0656       float xB_mc = -1000.;
0657       if (incQ2.GetSize()>0){
0658       Q2_mc = incQ2[0];
0659       xB_mc = incxB[0];
0660       }
0661       // find MC primary vertex
0662       int nMCPart = mcPartMass.GetSize();
0663 
0664       TVector3 vertex_mc(-999., -999., -999.);
0665       for(int imc=0; imc<nMCPart; imc++)
0666     {
0667       if(mcPartGenStatus[imc] == 4 && mcPartPdg[imc] == 11)
0668         {
0669           vertex_mc.SetXYZ(mcEndPointX[imc], mcEndPointY[imc], mcEndPointZ[imc]);
0670           break;
0671         }
0672     }
0673       hEventStat->Fill(0.5);
0674       hMcVtxX->Fill(vertex_mc.x());
0675       hMcVtxY->Fill(vertex_mc.y());
0676       hMcVtxZ->Fill(vertex_mc.z());
0677 
0678       // get RC primary vertex
0679       TVector3 vertex_rc(-999., -999., -999.);
0680       TVector3 err_vertex_rc(-999., -999., -999.);
0681       if(prim_vtx_index.GetSize()>0)
0682     {
0683       int rc_vtx_index = prim_vtx_index[0];
0684       vertex_rc.SetXYZ(CTVx[rc_vtx_index], CTVy[rc_vtx_index], CTVz[rc_vtx_index]);
0685       err_vertex_rc.SetXYZ(sqrt(CTVerr_xx[rc_vtx_index]), sqrt(CTVerr_yy[rc_vtx_index]), sqrt(CTVerr_zz[rc_vtx_index]));
0686     }
0687     hPullVtxX->Fill((vertex_rc.x()-vertex_mc.x())/err_vertex_rc.x()); 
0688     hPullVtxY->Fill((vertex_rc.y()-vertex_mc.y())/err_vertex_rc.y()); 
0689     hPullVtxZ->Fill((vertex_rc.z()-vertex_mc.z())/err_vertex_rc.z());   
0690 
0691       // map MC and RC particles
0692       int nAssoc = assocChRecID.GetSize();
0693       map<int, int> assoc_map_to_rc;
0694       map<int, int> assoc_map_to_mc;
0695       
0696       for(unsigned int rc_index=0; rc_index<rcMomPx.GetSize(); rc_index++)
0697     {
0698       // reco level eta using Real PID
0699       TVector3 mom(rcMomPx[rc_index], rcMomPy[rc_index], rcMomPz[rc_index]);
0700       hreco_eta->Fill(mom.Eta());
0701 
0702       if(abs(rcPdg[rc_index]) == 11) hreco_eta_e->Fill(mom.Eta());
0703       if(abs(rcPdg[rc_index]) == 211) hreco_eta_pi->Fill(mom.Eta());
0704       if(abs(rcPdg[rc_index]) == 321 ) hreco_eta_k->Fill(mom.Eta());
0705       if(abs(rcPdg[rc_index]) == 2212) hreco_eta_p->Fill(mom.Eta());
0706       
0707       // loop over the association to find the matched MC particle
0708       // with largest weight
0709       double max_weight = 0;
0710       int matched_mc_index = -1;
0711       for(int j=0; j<nAssoc; j++)
0712         {
0713           if(assocChRecID[j] != rc_index) continue;
0714           if(assocWeight[j] > max_weight)
0715         {
0716           max_weight = assocWeight[j];
0717           matched_mc_index = assocChSimID[j];
0718         }
0719         }
0720 
0721       // build the map
0722       assoc_map_to_rc[matched_mc_index] = rc_index;
0723       assoc_map_to_mc[rc_index] = matched_mc_index;
0724     }
0725 
0726       // Loop over primary particles
0727       int nMcPart = 0;
0728       for(int imc=0; imc<nMCPart; imc++)
0729     {
0730       // mc level eta using Truth PID
0731       TVector3 mom(mcMomPx[imc], mcMomPy[imc], mcMomPz[imc]);
0732           hmc_eta->Fill(mom.Eta());
0733 
0734           if(abs(mcPartPdg[imc]) == 11) hmc_eta_e->Fill(mom.Eta());
0735           if(abs(mcPartPdg[imc]) == 211) hmc_eta_pi->Fill(mom.Eta());
0736           if(abs(mcPartPdg[imc]) == 321 ) hmc_eta_k->Fill(mom.Eta());
0737           if(abs(mcPartPdg[imc]) == 2212) hmc_eta_p->Fill(mom.Eta());
0738 
0739       
0740       if(mcPartGenStatus[imc] == 1 && mcPartCharge[imc] != 0)
0741         {
0742           double dist = sqrt( pow(mcPartVx[imc]-vertex_mc.x(),2) + pow(mcPartVy[imc]-vertex_mc.y(),2) + pow(mcPartVz[imc]-vertex_mc.z(),2));      
0743           if(dist < 1e-4)
0744         {
0745           // count charged particles within |eta| < 3.5
0746           TVector3 mc_mom(mcMomPx[imc], mcMomPy[imc], mcMomPz[imc]);
0747           double mcEta = mc_mom.PseudoRapidity();
0748           if(fabs(mcEta) < 3.5) nMcPart++;
0749         }
0750         }
0751     }
0752       mult_charged = nMcPart;
0753       hMcMult->Fill(nMcPart);
0754       
0755       for(int imc=0; imc<nMCPart; imc++)
0756     {
0757       if(mcPartGenStatus[imc] == 1 && mcPartCharge[imc] != 0)
0758         {
0759           double dist = sqrt( pow(mcPartVx[imc]-vertex_mc.x(),2) + pow(mcPartVy[imc]-vertex_mc.y(),2) + pow(mcPartVz[imc]-vertex_mc.z(),2));      
0760           if(dist < 1e-4)
0761         {         
0762           // check if the MC particle is reconstructed
0763           int rc_index = -1;
0764           if(assoc_map_to_rc.find(imc) != assoc_map_to_rc.end()) rc_index = assoc_map_to_rc[imc];
0765 
0766           if(rc_index>=0)
0767             {
0768               TVector3 dcaToVtx = getDcaToVtx(rc_index, vertex_rc);
0769               
0770               int ip = -1;
0771               if(fabs(mcPartPdg[imc]) == 211) ip = 0;
0772               if(fabs(mcPartPdg[imc]) == 321) ip = 1;
0773               if(fabs(mcPartPdg[imc]) == 2212) ip = 2;
0774               if(ip>=0)
0775             {
0776               TVector3 mom(rcMomPx[rc_index], rcMomPy[rc_index], rcMomPz[rc_index]);
0777               if(ip<2)
0778                 {
0779                   hRcPrimPartLocaToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.Perp());
0780                   hRcPrimPartLocbToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.z());
0781                 }
0782 
0783               double fill1[] = {mom.Pt(), mom.Eta(), dcaToVtx.x(), nMcPart*1.};
0784               double fill2[] = {mom.Pt(), mom.Eta(), dcaToVtx.y(), nMcPart*1.};
0785               double fill3[] = {mom.Pt(), mom.Eta(), dcaToVtx.z(), nMcPart*1.};
0786               hPrimTrkDcaToRCVtx[ip][0]->Fill(fill1);
0787               hPrimTrkDcaToRCVtx[ip][1]->Fill(fill2);
0788               hPrimTrkDcaToRCVtx[ip][2]->Fill(fill3);
0789             }
0790             }
0791         }
0792             }
0793     }
0794       
0795       // Filling some values: 
0796       // For Vertex
0797       Vertex_x.clear();
0798       Vertex_y.clear();
0799       Vertex_z.clear();
0800       Vertex_ndf.clear();
0801       Vertex_chi2.clear();
0802       VertexErr_xx.clear();
0803       VertexErr_yy.clear();
0804       VertexErr_zz.clear();
0805       
0806       for (unsigned int ivtx = 0; ivtx < CTVx.GetSize(); ++ivtx) {
0807     Vertex_x.push_back(CTVx[ivtx]);
0808     Vertex_y.push_back(CTVy[ivtx]);
0809     Vertex_z.push_back(CTVz[ivtx]);
0810     Vertex_ndf.push_back(CTVndf[ivtx]);
0811     Vertex_chi2.push_back(CTVchi2[ivtx]);
0812     VertexErr_xx.push_back(CTVerr_xx[ivtx]);
0813     VertexErr_yy.push_back(CTVerr_yy[ivtx]);
0814     VertexErr_zz.push_back(CTVerr_zz[ivtx]);        
0815         }
0816         
0817       //  For Scattered electron
0818       
0819       if(ScatElecRecoId.GetSize() > 0){
0820     int iscatR = ScatElecRecoId[0];
0821     ScatteredERecId = iscatR;    
0822       }else{ScatteredERecId = -999;}
0823       
0824       
0825       if(ScatElecGenId.GetSize() > 0){
0826     int iscatS = ScatElecGenId[0];
0827     ScatteredEGenId = iscatS;    
0828       }else{ScatteredEGenId = -999;}
0829         
0830       // other reco and gen quantities
0831       if(EvtQ2.GetSize() > 0){
0832     float q2valueR = EvtQ2[0];
0833     EventQ2 = q2valueR;    
0834       }else{EventQ2 = -999;}
0835       
0836       if(Evtx.GetSize() > 0){
0837     float xvalueR = Evtx[0];
0838     Eventx = xvalueR;    
0839       }else{Eventx = -999;}
0840       
0841       if(EvtQ2Gen.GetSize() > 0){
0842     float q2valueG = EvtQ2Gen[0];
0843     EventQ2Gen = q2valueG;    
0844       }else{EventQ2Gen = -999;}
0845         
0846       if(EvtxGen.GetSize() > 0){
0847     float xvalueG = EvtxGen[0];
0848     EventxGen = xvalueG;    
0849       }else{EventxGen = -999;}
0850       
0851       
0852       // look for D0
0853       bool hasD0 = false;
0854       vector<int> mc_index_D0_pi;
0855       vector<int> mc_index_D0_k;
0856       mc_index_D0_pi.clear();
0857       mc_index_D0_k.clear();
0858       
0859       for(int imc=0; imc<nMCPart; imc++)
0860     {
0861       if(fabs(mcPartPdg[imc]) != 421) continue;
0862       hEventStat->Fill(1.5);
0863    //----------------------------------- Fill D0 gen tree with properties-----------------------------------------
0864     TLorentzVector d0_vec;
0865 
0866     d0_px = mcMomPx[imc];
0867     d0_py = mcMomPy[imc];
0868     d0_pz = mcMomPz[imc];
0869     d0_mass = mcPartMass[imc];
0870 
0871     d0_vec.SetXYZM(d0_px, d0_py, d0_pz, d0_mass);
0872 
0873     d0_pt  = d0_vec.Pt();
0874     d0_eta = d0_vec.Eta();
0875     d0_y   = d0_vec.Rapidity();
0876     tree_D0->Fill();
0877     //--------------------------------------------------------------------------------------------------------------------
0878       
0879       int nDuaghters = mcPartDaughter_end[imc]-mcPartDaughter_begin[imc];
0880       if(nDuaghters!=2) continue;
0881 
0882       // find D0 that decay into pi+K
0883       bool is_pik_decay = false;      
0884       int daug_index_1 = mcPartDaughter_index[mcPartDaughter_begin[imc]];
0885       int daug_index_2 = mcPartDaughter_index[mcPartDaughter_begin[imc]+1];
0886       int daug_pdg_1 = mcPartPdg[daug_index_1];
0887       int daug_pdg_2 = mcPartPdg[daug_index_2];
0888       if( (fabs(daug_pdg_1)==321 && fabs(daug_pdg_2)==211) || (fabs(daug_pdg_1)==211 && fabs(daug_pdg_2)==321) )
0889         {
0890           is_pik_decay = true;
0891           
0892         }
0893       if(!is_pik_decay) continue;
0894       if(fabs(daug_pdg_1)==211)
0895         {
0896           mc_index_D0_pi.push_back(daug_index_1);
0897           mc_index_D0_k.push_back(daug_index_2);
0898         }
0899       else
0900         {
0901           mc_index_D0_pi.push_back(daug_index_2);
0902           mc_index_D0_k.push_back(daug_index_1);
0903         }
0904       hasD0 = true;
0905       hEventStat->Fill(2.5);
0906 
0907       // D0 kinematics
0908       TLorentzVector mc_mom_vec;
0909       mc_mom_vec.SetXYZM(mcMomPx[imc], mcMomPy[imc], mcMomPz[imc], mcPartMass[imc]);
0910       double mcRap = mc_mom_vec.Rapidity();
0911       double mcPt = mc_mom_vec.Pt();
0912       hMCD0PtRap->Fill(mcRap, mcPt);
0913 
0914       // decay dauther kinematics
0915       for(int ip = 0; ip<2; ip++)
0916         {
0917           int mc_part_index;
0918           if(ip==0) mc_part_index = mc_index_D0_pi[mc_index_D0_pi.size()-1];
0919           if(ip==1) mc_part_index = mc_index_D0_k[mc_index_D0_k.size()-1];
0920           
0921           TLorentzVector mc_part_vec;
0922           mc_part_vec.SetXYZM(mcMomPx[mc_part_index], mcMomPy[mc_part_index], mcMomPz[mc_part_index], mcPartMass[mc_part_index]);
0923           if(ip==0) hMcPiPtEta->Fill(mc_part_vec.Eta(), mc_part_vec.Pt());
0924           if(ip==1) hMcKPtEta->Fill(mc_part_vec.Eta(), mc_part_vec.Pt());
0925           
0926           int rc_part_index = -1;
0927           if(assoc_map_to_rc.find(mc_part_index) != assoc_map_to_rc.end()) rc_part_index = assoc_map_to_rc[mc_part_index];
0928           if(rc_part_index>=0)
0929         {
0930           TVector3 dcaToVtx = getDcaToVtx(rc_part_index, vertex_rc);
0931           
0932           TVector3 mom(rcMomPx[rc_part_index], rcMomPy[rc_part_index], rcMomPz[rc_part_index]);
0933           hRcSecPartLocaToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.Pt());
0934           hRcSecPartLocbToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.z());
0935           
0936           //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]);
0937         }
0938         }
0939     }
0940 
0941       // Get reconstructed pions and kaons
0942       hNRecoVtx->Fill(CTVx.GetSize());
0943       const int pid_mode = 1; // 0 - truth; 1 - realistic
0944       vector<unsigned int> pi_index;
0945       vector<unsigned int> k_index;
0946       pi_index.clear();
0947       k_index.clear();
0948       for(unsigned int rc_index=0; rc_index<rcMomPx.GetSize(); rc_index++)
0949     {     
0950       if(pid_mode==0)
0951         {
0952           int iSimPartID = -1;
0953           if(assoc_map_to_mc.find(rc_index) != assoc_map_to_mc.end()) iSimPartID = assoc_map_to_mc[rc_index];
0954           if(iSimPartID>=0)
0955         {
0956           if(fabs(mcPartPdg[iSimPartID]) == 211) pi_index.push_back(rc_index);
0957           if(fabs(mcPartPdg[iSimPartID]) == 321) k_index.push_back(rc_index);
0958         }
0959         }
0960       else if(pid_mode==1)
0961         {
0962           if(fabs(rcPdg[rc_index]) == 211) pi_index.push_back(rc_index);
0963           if(fabs(rcPdg[rc_index]) == 321) k_index.push_back(rc_index);
0964         }
0965     }
0966 
0967       // pair pion and kaon
0968       for(unsigned int i=0; i<pi_index.size(); i++)
0969     {
0970 
0971       //Check for pion using Truth level info
0972       /*int mcIdx_pii = -1;
0973       if(assoc_map_to_mc.find(pi_index[i]) != assoc_map_to_mc.end()) mcIdx_pii = assoc_map_to_mc[pi_index[i]];
0974       
0975       if(mcIdx_pii >=0 ){ 
0976         if(abs(mcPartPdg[mcIdx_pii]) != 211 ) continue ;
0977       }
0978       */
0979       bool not_pi = false;
0980       TVector3 dcaToVtx = getDcaToVtx(pi_index[i], vertex_rc);
0981       std::array<float, 21>& cov_pion = rcTrkCov->At(pi_index[i]);
0982       int q_pion = rcCharge[pi_index[i]];
0983        
0984       for(unsigned int j=0; j<k_index.size(); j++)
0985         {
0986           bool not_k = false;
0987           //Check for kaon using Truth level info
0988           /*int mcIdx_kk = -1;
0989           if(assoc_map_to_mc.find(k_index[i]) != assoc_map_to_mc.end()) mcIdx_kk = assoc_map_to_mc[k_index[i]];
0990 
0991           if(mcIdx_kk>=0){
0992         if(abs(mcPartPdg[mcIdx_kk]) != 321) continue ;
0993           }
0994           */
0995           TVector3 dcaToVtx2 = getDcaToVtx(k_index[j], vertex_rc);
0996           std::array<float, 21>& cov_kaon = rcTrkCov->At(k_index[j]); 
0997           int q_kaon = rcCharge[k_index[j]];
0998           
0999           if(rcCharge[pi_index[i]]*rcCharge[k_index[j]]<0)
1000         {
1001           //printf("[i] Check pair (%d, %d)\n", pi_index[i], k_index[j]);
1002           // -- only look at unlike-sign pi+k pair
1003           bool is_D0_pik = false;
1004           int mc_index_pi = -1, mc_index_k = -1;
1005           if(assoc_map_to_mc.find(pi_index[i]) != assoc_map_to_mc.end()) mc_index_pi = assoc_map_to_mc[pi_index[i]];
1006           if(assoc_map_to_mc.find(k_index[j])  != assoc_map_to_mc.end()) mc_index_k  = assoc_map_to_mc[k_index[j]];
1007 
1008           for(unsigned int k=0; k<mc_index_D0_pi.size(); k++)
1009             {
1010               if(mc_index_pi==mc_index_D0_pi[k] && mc_index_k==mc_index_D0_k[k])
1011             {
1012               is_D0_pik = true;
1013               break;
1014             }
1015                 }
1016 
1017           float dcaDaughters, cosTheta, decayLength, V0DcaToVtx, cosTheta_xy, sigma_vtx;
1018           TVector3 decayVertex, decayVertex_ana; 
1019           double chi2_ndf;
1020           double err_Par[5];  // or whatever size is appropriate
1021           //double* err_Par = errParArray;
1022           TLorentzVector parent = getPairParent(pi_index[i], k_index[j], vertex_rc, dcaDaughters, cosTheta, cosTheta_xy, decayLength, V0DcaToVtx, sigma_vtx,decayVertex,decayVertex_ana,chi2_ndf, err_Par);
1023           hchi2_vtx->Fill(chi2_ndf);
1024 
1025           // Storing rc index of Pion and Kaon
1026               int pii_index = pi_index[i]; // rc index of pion
1027           int kk_index =  k_index[j];  // rc index of kaon
1028 
1029               //Signal D0
1030           if(is_D0_pik) 
1031             {
1032            TVector3 MCVertex_Kaon(mcPartVx[mc_index_k], mcPartVy[mc_index_k], mcPartVz[mc_index_k]);
1033            TVector3 MCVertex_Pion(mcPartVx[mc_index_pi], mcPartVy[mc_index_pi], mcPartVz[mc_index_pi]);
1034            
1035           // printf("Signal MC Vertex Kaon = (%f, %f, %f)\n",MCVertex_Kaon.X(), MCVertex_Kaon.Y(), MCVertex_Kaon.Z());      
1036           // printf("Signal MC Vertex Pion = (%f, %f, %f)\n",MCVertex_Pion.X(), MCVertex_Pion.Y(), MCVertex_Pion.Z());
1037               hRes_SVx_Helixfit->Fill(parent.Pt(), parent.Rapidity(), (decayVertex.X()-MCVertex_Kaon.X()*0.1)*10);
1038               hRes_SVy_Helixfit->Fill(parent.Pt(), parent.Rapidity(), (decayVertex.Y()-MCVertex_Kaon.Y()*0.1)*10);
1039               hRes_SVz_Helixfit->Fill(parent.Pt(), parent.Rapidity(), (decayVertex.Z()-MCVertex_Kaon.Z()*0.1)*10);  
1040               
1041               hRes_SVx_Helixana->Fill(parent.Pt(), parent.Rapidity(), (decayVertex_ana.X()-MCVertex_Kaon.X()*0.1)*10);
1042               hRes_SVy_Helixana->Fill(parent.Pt(), parent.Rapidity(), (decayVertex_ana.Y()-MCVertex_Kaon.Y()*0.1)*10);
1043               hRes_SVz_Helixana->Fill(parent.Pt(), parent.Rapidity(), (decayVertex_ana.Z()-MCVertex_Kaon.Z()*0.1)*10);
1044               
1045               hRes_SVxy_Helixfit->Fill(parent.Pt(), parent.Rapidity(), (decayVertex.Perp()-MCVertex_Kaon.Perp()*0.1)*10);
1046               hRes_SVxy_Helixana->Fill(parent.Pt(), parent.Rapidity(), (decayVertex_ana.Perp()-MCVertex_Kaon.Perp()*0.1)*10);
1047               
1048               hRes_SVx_Helixfit_pull->Fill(parent.Pt(), parent.Rapidity(), ((decayVertex.X()-MCVertex_Kaon.X()*0.1))/(err_Par[0]));
1049               hRes_SVy_Helixfit_pull->Fill(parent.Pt(), parent.Rapidity(), ((decayVertex.Y()-MCVertex_Kaon.Y()*0.1))/(err_Par[1]));
1050               hRes_SVz_Helixfit_pull->Fill(parent.Pt(), parent.Rapidity(), ((decayVertex.Z()-MCVertex_Kaon.Z()*0.1))/(err_Par[2])); 
1051               
1052               hchi2_vtx_sig->Fill(chi2_ndf);              
1053                                
1054               hEventStat->Fill(3.5);
1055               if (q_kaon == -1 && q_pion == 1)
1056               hEventStat->Fill(4.5);   // D0
1057               else if (q_kaon == 1 && q_pion == -1)
1058               hEventStat->Fill(5.5);   // D0bar
1059               h3PairDca12[0]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters);
1060               h3PairCosTheta[0]->Fill(parent.Pt(), parent.Rapidity(), cosTheta);
1061               h3PairDca[0]->Fill(parent.Pt(), parent.Rapidity(), V0DcaToVtx);
1062               h3PairDecayLength[0]->Fill(parent.Pt(), parent.Rapidity(), decayLength);
1063               //printf("Signal: dca12 = %2.4f, cosTheta = %2.4f, D0dca = %2.4f, decay = %2.4f\n", dcaDaughters, cosTheta, V0DcaToVtx, decayLength);
1064               h3InvMass[0][0]->Fill(parent.Pt(), parent.Rapidity(), parent.M());
1065                
1066            // Toplogical Variables for Signal
1067               d0_pi_sig = dcaToVtx.Mag();
1068               d0_k_sig = dcaToVtx2.Mag();
1069               d0xy_pi_sig = dcaToVtx.Perp();
1070               signif_d0xy_pi_sig = d0xy_pi_sig/sqrt(cov_pion[0]);
1071               d0xy_k_sig = dcaToVtx2.Perp();
1072               signif_d0xy_k_sig = d0xy_k_sig/sqrt(cov_kaon[0]);
1073               sum_d0xy_sig = sqrt(d0xy_pi_sig*d0xy_pi_sig+d0xy_k_sig*d0xy_k_sig);             
1074               dca_12_sig = dcaDaughters;
1075               dca_D0_sig = V0DcaToVtx;
1076               decay_length_sig = decayLength;
1077               costheta_sig = cosTheta;
1078               costhetaxy_sig = cosTheta_xy;
1079               pt_D0_sig = parent.Pt();
1080               y_D0_sig = parent.Rapidity();
1081               mass_D0_sig = parent.M();
1082               sigma_vtx_sig = sigma_vtx;
1083               mult_sig = mult_charged;
1084               chi2_dca_sig = chi2_ndf;
1085               xB_sig = xB_mc;
1086               Q2_sig = Q2_mc;        
1087                 
1088 //Jet Clustering for the signal D0 (signal pik pair)   
1089 // Build pseudojets (using reco-particles)
1090 
1091               std::vector<fastjet::PseudoJet> particles_reco; // Reco PseudoJet type Vector
1092               int c1 =0;                                                           
1093               
1094               for (unsigned int i = 0; i < rcMomPx.GetSize(); ++i) { 
1095             
1096             bool is_PiorK = false;
1097             if(i == pii_index || i == kk_index){  //looking for pi and k daughters of signal D0  
1098               is_PiorK = true;
1099             }
1100  
1101             if(is_PiorK) continue;
1102  
1103             TVector3 mom(rcMomPx[i], rcMomPy[i],rcMomPz[i]);
1104             if ( mom.Pt() < minCstPt || mom.Pt() > maxCstPt ) continue;
1105             if ( TrkRecoNhits[i] < nhitcut ) continue;
1106             
1107             if ( removeelectrons == 1 ){
1108               // Find electron
1109               int chargePartIndex = i; 
1110               int elecIndex = -1;
1111               float elecIndexWeight = -1.0;
1112               for(unsigned int itrkass = 0; itrkass < assocWeight.GetSize(); itrkass++){ // Loop Over All ReconstructedChargedParticleAssociations
1113                 
1114                 if( assocWeight.GetSize() > 0 ){
1115                   if(TrkPartAssocRec[itrkass] == chargePartIndex){ // Select Entry Matching the ReconstructedChargedParticle Index
1116                 if(assocWeight[itrkass] > elecIndexWeight){ // Find Particle with Greatest Weight = Contributed Most Hits to Track
1117                   elecIndex = TrkPartAssocSim[itrkass]; // Get Index of MCParticle Associated with ReconstructedChargedParticle
1118                   elecIndexWeight = assocWeight[itrkass];
1119                 }
1120                   }
1121                 }
1122               }
1123 
1124               if(mcPartPdg[elecIndex] == 11) continue;
1125             }
1126             
1127             if ( removeelectrons == 2 && i == ScatteredERecId ) continue; 
1128 
1129             fastjet::PseudoJet p(rcMomPx[i], rcMomPy[i], rcMomPz[i], TrkRecoE[i]);
1130             p.set_user_index(i);
1131             particles_reco.push_back(p);
1132               }
1133 
1134               //----------Introducing D0 parent into recotype PseudoJets---------------------   
1135               float dcaDaughters, cosTheta, decayLength, V0DcaToVtx, cosTheta_xy, sigma_vtx;
1136               TVector3 decayVertex, decayVertex_ana; 
1137               double chi2_ndf;
1138               double err_Par[5];  // or whatever size is appropriate
1139               //double* err_Par = errParArray;
1140               TLorentzVector parent = getPairParent(pii_index,kk_index, vertex_rc, dcaDaughters, cosTheta, cosTheta_xy, decayLength, V0DcaToVtx, sigma_vtx,decayVertex,decayVertex_ana,chi2_ndf, err_Par);
1141               pt_parent_sig = parent.Pt();
1142   
1143               if ( parent.Pt() < minCstPt || parent.Pt() > maxCstPt ){ // if D0 parent don't fit in Jet Constituent Pt limits => introduce pi-k daughters  to pseudojets                                                //we should not do this, it would give jets with no D0!!!
1144     
1145             for (int l=0;l<2;l++){ 
1146               int idx;
1147               if (l==0) idx = pii_index;
1148               if(l==1)  idx = kk_index;
1149             
1150             
1151               TVector3 mom(rcMomPx[idx], rcMomPy[idx],rcMomPz[idx]);
1152               if ( mom.Pt() < minCstPt || mom.Pt() > maxCstPt ) continue;
1153               if ( TrkRecoNhits[i] < nhitcut ) continue;
1154               
1155               fastjet::PseudoJet p(rcMomPx[idx], rcMomPy[idx], rcMomPz[idx], TrkRecoE[idx]);
1156               p.set_user_index(idx);
1157               particles_reco.push_back(p);
1158               
1159             } 
1160               }
1161               
1162               if ( parent.Pt() > minCstPt && parent.Pt() < maxCstPt ){
1163             fastjet::PseudoJet p(parent.Px(), parent.Py(),parent.Pz(), parent.E());  
1164             p.set_user_index(-99999);                                 // setting D0 particle index = negative integer (-99999) 
1165             particles_reco.push_back(p);
1166               }
1167               
1168       // Define algorithm
1169     fastjet::JetAlgorithm algo = fastjet::antikt_algorithm;
1170     fastjet::RecombinationScheme scheme = fastjet::E_scheme;
1171       // Jet definition
1172     fastjet::JetDefinition jet_def(algo, R_value, scheme);
1173     fastjet::GhostedAreaSpec ghost_spec(ghostMaxRap, numGhostRepeat, ghostArea);
1174     fastjet::AreaType atype = fastjet::active_area;
1175     fastjet::AreaDefinition area_def(atype, ghost_spec);
1176                         
1177       // Clear all vectors for the Reco Jet Clustering for new Signal D0 
1178       RecoJet_pt.clear();
1179       RecoJet_eta.clear();
1180       RecoJet_rapidity.clear();
1181       RecoJet_phi.clear();
1182       RecoJet_E.clear();
1183       RecoJet_M.clear();
1184       RecoJet_rapi.clear();
1185       RecoJet_hasElectron.clear();
1186       RecoJet_maxPtPart_pt.clear();
1187       RecoJet_constituent_pt.clear(); 
1188       RecoJet_constituent_eta.clear();
1189       RecoJet_constituent_phi.clear(); 
1190       RecoJet_constituent_nhits.clear();
1191       RecoJet_constituent_pdgid.clear();
1192       RecoJet_constituent_pdgidTruth.clear();
1193       RecoJet_constituent_idx.clear();
1194           RecoJet_constituent_energy.clear();
1195       RecoJet_constituent_rapi.clear();
1196       
1197       EVETMULTRECO = (int)rcMomPx.GetSize();        
1198 
1199       //Reco Jet clustering(Sig)
1200       fastjet::ClusterSequenceArea cs_reco(particles_reco, jet_def, area_def);
1201       std::vector<fastjet::PseudoJet> jets_reco = fastjet::sorted_by_pt(cs_reco.inclusive_jets(minJetPt));
1202       for (auto &jet : jets_reco) {
1203         RecoJet_pt.push_back(jet.pt());
1204         RecoJet_eta.push_back(jet.eta());
1205         RecoJet_rapidity.push_back(jet.rapidity());
1206         RecoJet_phi.push_back(jet.phi_std());
1207         RecoJet_E.push_back(jet.e());
1208         RecoJet_M.push_back(jet.m());
1209         RecoJet_rapi.push_back(jet.rap());
1210         
1211 
1212         bool hasElectron = false;
1213         float maxPtReco = -1.0;
1214         std::vector<float> cpt, ceta, cphi, cenergy,crapi;
1215         std::vector<int> chits, cpdgid, cpdgidtruth,cindex;
1216         cpt.clear(); ceta.clear(); cphi.clear(); chits.clear(); cpdgid.clear(); cpdgidtruth.clear(), cindex.clear(), cenergy.clear(),crapi.clear();
1217         
1218         fastjet::PseudoJet D0_parent;
1219 
1220         for (auto &c : jet.constituents()) {
1221           int idx = c.user_index();
1222           
1223         //Look for the Signal D0 constituent of the Jet 
1224           bool is_D0_parent = false;
1225           
1226           if(idx == -99999){ //check for the signal D0
1227           
1228         for(const auto& p : particles_reco ){ // loop over all the jet particles
1229         
1230           if (p.user_index() == idx) {
1231             D0_parent = p;
1232             is_D0_parent = true;
1233             break;
1234           }
1235         }
1236                       
1237         cpt.push_back(D0_parent.pt());
1238         ceta.push_back(D0_parent.eta());
1239         cphi.push_back(D0_parent.phi());
1240         cenergy.push_back(D0_parent.E());
1241         chits.push_back(0);           // 0 hits for neutral D0 
1242         cpdgid.push_back(421);        // PID = 421 for D0
1243         cpdgidtruth.push_back(421);
1244         cindex.push_back(idx);
1245         crapi.push_back(D0_parent.rap());
1246                                                                   // Here, hasElectron = false
1247         if (D0_parent.pt() > maxPtReco) maxPtReco = D0_parent.pt();   
1248           }
1249           if (is_D0_parent) continue; // skip the rest of loop for signal D0 particle jet
1250             
1251           TVector3 v3(rcMomPx[idx], rcMomPy[idx], rcMomPz[idx]);
1252           float rcEnergy = TrkRecoE[idx];
1253           TLorentzVector lv(rcMomPx[idx],  rcMomPy[idx], rcMomPz[idx], rcEnergy);
1254           
1255           cpt.push_back(v3.Pt());
1256           ceta.push_back(v3.Eta());
1257           cphi.push_back(v3.Phi());
1258           chits.push_back(TrkRecoNhits[idx]);
1259           cindex.push_back(idx);
1260           cenergy.push_back(rcEnergy);
1261           crapi.push_back(lv.Rapidity());
1262           if (v3.Pt() > maxPtReco) maxPtReco = v3.Pt();
1263           
1264           // Find electron
1265           int chargePartIndex = idx; 
1266           int elecIndex = -1;
1267           float elecIndexWeight = -1.0;
1268           for(unsigned int itrkass = 0; itrkass < assocWeight.GetSize(); itrkass++){ // Loop Over All ReconstructedChargedParticleAssociations
1269 
1270         if( assocWeight.GetSize() > 0 ){
1271           if(TrkPartAssocRec[itrkass] == chargePartIndex){ // Select Entry Matching the ReconstructedChargedParticle Index
1272             if(assocWeight[itrkass] > elecIndexWeight){ // Find Particle with Greatest Weight = Contributed Most Hits to Track
1273               elecIndex = TrkPartAssocSim[itrkass]; // Get Index of MCParticle Associated with ReconstructedChargedParticle
1274               elecIndexWeight = assocWeight[itrkass];
1275             }
1276           }
1277         }
1278           }
1279 
1280           if(mcPartPdg[elecIndex] == 11){
1281         hasElectron = true;
1282           }
1283         
1284           cpdgid.push_back(rcPdg[idx]);
1285           cpdgidtruth.push_back(mcPartPdg[elecIndex]);
1286         }
1287             
1288         RecoJet_constituent_pt.push_back(cpt);  
1289         RecoJet_constituent_eta.push_back(ceta);
1290         RecoJet_constituent_phi.push_back(cphi);
1291         RecoJet_constituent_nhits.push_back(chits);
1292         RecoJet_constituent_pdgid.push_back(cpdgid);
1293         RecoJet_constituent_pdgidTruth.push_back(cpdgidtruth);
1294         RecoJet_hasElectron.push_back(hasElectron);
1295         RecoJet_maxPtPart_pt.push_back(maxPtReco);
1296         RecoJet_constituent_idx.push_back(cindex);
1297             RecoJet_constituent_energy.push_back(cenergy);
1298             RecoJet_constituent_rapi.push_back(crapi);
1299         
1300       } 
1301 
1302 // Reading Reco Jet(Signal):
1303 
1304     for(unsigned int i=0; i<RecoJet_pt.size(); i++){  //loop over reco jets  
1305       
1306       int  jetParticles = RecoJet_constituent_idx[i].size();  // Total no. of constituents of ith jet 
1307      
1308       float jetPt  = RecoJet_pt[i];
1309       float jetEta = RecoJet_eta[i];
1310       float jetRapidity = RecoJet_rapi[i];
1311       float jetPhi = RecoJet_phi[i];
1312       float jetE = RecoJet_E[i];
1313      
1314       TLorentzVector LvJet;
1315       LvJet.SetPtEtaPhiE(jetPt,jetEta, jetPhi, jetE); // TLorentzVector of Jet 
1316       TVector3 jetMom = LvJet.Vect(); // Jet Momentum Vector
1317           bool c2= false;
1318      
1319       for(unsigned int j=0; j<RecoJet_constituent_idx[i].size(); j++){ //Loop over all constituent particles of the ith Jet
1320 
1321         int particle_idx = RecoJet_constituent_idx[i][j]; // index of jth particle of ith jet
1322         bool is_D0 = false;
1323             
1324         //check if constituent particle is D0
1325         if(particle_idx == -99999){        
1326           is_D0 = true;
1327           c2 = true;
1328         }
1329      
1330         if(is_D0){  // if Jet constituent is the signal D0
1331 
1332           float D0Pt  = RecoJet_constituent_pt[i][j];
1333           float D0Eta = RecoJet_constituent_eta[i][j];
1334           float D0Phi = RecoJet_constituent_phi[i][j];
1335           float D0E = RecoJet_constituent_energy[i][j];
1336           
1337           TLorentzVector lvD0; 
1338           lvD0.SetPtEtaPhiE(D0Pt, D0Eta, D0Phi, D0E);
1339           
1340           TVector3 D0_Mom;
1341           D0_Mom.SetPtEtaPhi(D0Pt, D0Eta, D0Phi);
1342               float D0Mass = lvD0.M();
1343           float D0Rapidity = lvD0.Rapidity();
1344     
1345           float sig_z = float(jetMom.Dot(D0_Mom))/float(jetMom.Dot(jetMom)) ;  // Fragmentation Variable z for the sinal D0
1346               z_sig =sig_z;
1347               h3sig_z->Fill(sig_z,D0Rapidity,D0Mass);
1348               
1349                      
1350           float dPhi = jetPhi - D0Phi;
1351               dPhi = TVector2::Phi_mpi_pi(dPhi);
1352               float dEta = jetEta - D0Eta;
1353               float deltaR = TMath::Sqrt(dEta*dEta + dPhi*dPhi); // sig_D0 distance from JetAxis in eta-Phi space
1354               
1355               float angleRad = jetMom.Angle(D0_Mom); // angle bw sig_D0 momentum vector and jetMomentum vector
1356               float angleDeg = angleRad * 180.0 / TMath::Pi(); 
1357                         
1358               angle_sig = angleDeg;
1359               dr_sig = deltaR;
1360               etajet_sig = jetEta;
1361           pTjet_sig = jetPt;
1362           D0eta_sig = D0Eta;
1363           
1364           break;//leaving the JetConstituent loop because we found the only signal D0 present
1365        }     
1366     }
1367       if(c2) break; // leaving the Jet loop because we found the only signal D0 present
1368       }
1369                   
1370       tree_sig->Fill();
1371        
1372             }
1373 
1374           //Background:   
1375           else
1376             {
1377 
1378               hEventStat->Fill(6.5);   
1379               hchi2_vtx_bkg->Fill(chi2_ndf);                    
1380               h3PairDca12[1]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters);
1381               h3PairCosTheta[1]->Fill(parent.Pt(), parent.Rapidity(), cosTheta);
1382               h3PairDca[1]->Fill(parent.Pt(), parent.Rapidity(), V0DcaToVtx);
1383               h3PairDecayLength[1]->Fill(parent.Pt(), parent.Rapidity(), decayLength);
1384               
1385               //printf("Bkg: dca12 = %2.4f, cosTheta = %2.4f, D0dca = %2.4f, decay = %2.4f\n", dcaDaughters, cosTheta, V0DcaToVtx, decayLength);
1386               h3InvMass[1][0]->Fill(parent.Pt(), parent.Rapidity(), parent.M());
1387 
1388               // Toplogical Variables for Bkg
1389               d0_pi_bkg = dcaToVtx.Mag();
1390               d0_k_bkg = dcaToVtx2.Mag();
1391               d0xy_pi_bkg = dcaToVtx.Perp();
1392               signif_d0xy_pi_bkg = d0xy_pi_bkg/sqrt(cov_pion[0]);
1393               d0xy_k_bkg = dcaToVtx2.Perp();
1394               signif_d0xy_k_bkg = d0xy_k_bkg/sqrt(cov_kaon[0]);
1395               sum_d0xy_bkg = sqrt(d0xy_pi_bkg*d0xy_pi_bkg+d0xy_k_bkg*d0xy_k_bkg);             
1396               dca_12_bkg = dcaDaughters;
1397               dca_D0_bkg = V0DcaToVtx;
1398               decay_length_bkg = decayLength;
1399               costheta_bkg = cosTheta;
1400               costhetaxy_bkg = cosTheta_xy;
1401               pt_D0_bkg = parent.Pt();
1402               y_D0_bkg = parent.Rapidity();
1403               mass_D0_bkg = parent.M();
1404               sigma_vtx_bkg = sigma_vtx;
1405               mult_bkg = mult_charged;
1406                       chi2_dca_bkg = chi2_ndf;
1407                   xB_bkg = xB_mc;          
1408                   Q2_bkg = Q2_mc; 
1409               
1410 //Jet Clustering for the Bkg D0 (Bkg pi-k pair)      
1411 //Build pseudojets (recotype)
1412 
1413 
1414     std::vector<fastjet::PseudoJet> particles_reco;                                                      
1415         
1416         for (unsigned int i = 0; i < rcMomPx.GetSize(); ++i) { 
1417                                                                              
1418       
1419     TVector3 mom(rcMomPx[i], rcMomPy[i],rcMomPz[i]);
1420     if ( mom.Pt() < minCstPt || mom.Pt() > maxCstPt ) continue;
1421     if ( TrkRecoNhits[i] < nhitcut ) continue;
1422     
1423     if ( removeelectrons == 1 ){ // ******
1424       // Find electron
1425       int chargePartIndex = i; 
1426       int elecIndex = -1;
1427       float elecIndexWeight = -1.0;
1428       for(unsigned int itrkass = 0; itrkass < assocWeight.GetSize(); itrkass++){ // Loop Over All ReconstructedChargedParticleAssociations
1429 
1430         if( assocWeight.GetSize() > 0 ){
1431         if(TrkPartAssocRec[itrkass] == chargePartIndex){ // Select Entry Matching the ReconstructedChargedParticle Index
1432         if(assocWeight[itrkass] > elecIndexWeight){ // Find Particle with Greatest Weight = Contributed Most Hits to Track
1433           elecIndex = TrkPartAssocSim[itrkass]; // Get Index of MCParticle Associated with ReconstructedChargedParticle
1434           elecIndexWeight = assocWeight[itrkass];
1435         }
1436           }
1437         }
1438       }
1439       if(mcPartPdg[elecIndex] == 11){
1440 
1441         if( i == pii_index || i == kk_index ){
1442           if(i == pii_index) not_pi = true;
1443           if(i == kk_index) not_k = true;
1444           break;
1445         }
1446         continue;
1447       }
1448     }
1449     
1450     if ( removeelectrons == 2 && i == ScatteredERecId ) continue; // *****
1451 
1452     fastjet::PseudoJet p(rcMomPx[i], rcMomPy[i], rcMomPz[i], TrkRecoE[i]);
1453     p.set_user_index(i);
1454     particles_reco.push_back(p);
1455         }
1456 
1457     if(not_k) continue;  // skip to next rc k
1458     if(not_pi) goto next_pi; // get out of k loop and move to next pi
1459        
1460     //Introducing Signal D0 parent into recotype PseudoJets 
1461     float dcaDaughters, cosTheta, decayLength, V0DcaToVtx, cosTheta_xy, sigma_vtx;
1462     TVector3 decayVertex, decayVertex_ana; 
1463     double chi2_ndf;
1464     double err_Par[5];  // or whatever size is appropriate
1465         //double* err_Par = errParArray;
1466         TLorentzVector parent = getPairParent(pii_index,kk_index, vertex_rc, dcaDaughters, cosTheta, cosTheta_xy, decayLength, V0DcaToVtx, sigma_vtx,decayVertex,decayVertex_ana,chi2_ndf, err_Par);
1467     pt_parent_bkg = parent.Pt();
1468 
1469 //if bkg D0 Pt < minCstPt or > maxCstPt, then instead of D0 parent, introducing its pi-k daughters to recotype pseudoJet    
1470     if ( parent.Pt() < minCstPt || parent.Pt() > maxCstPt ){ //we should not do this as this will may give jet with no D0
1471         
1472       for (int l=0;l<2;l++){ 
1473         int idx;
1474             
1475             if (l==0){
1476               idx = pii_index;
1477             }
1478       
1479       if(l==1){
1480         idx = kk_index;
1481       }
1482       
1483         TVector3 mom(rcMomPx[idx], rcMomPy[idx],rcMomPz[idx]);
1484         if ( mom.Pt() < minCstPt || mom.Pt() > maxCstPt ) continue;
1485         if ( TrkRecoNhits[i] < nhitcut ) continue;
1486 
1487         fastjet::PseudoJet p(rcMomPx[idx], rcMomPy[idx], rcMomPz[idx], TrkRecoE[idx]);
1488         p.set_user_index(idx);
1489         particles_reco.push_back(p);
1490     
1491     } 
1492     }
1493 
1494 //If bkg D0 Pt is b/w minCstPt and maxCstPt
1495     if ( parent.Pt() > minCstPt && parent.Pt() < maxCstPt ){ 
1496         fastjet::PseudoJet p(parent.Px(), parent.Py(),parent.Pz(), parent.E());  
1497         p.set_user_index(-99999);// setting signal D0 index = negative integer (-99999) 
1498         particles_reco.push_back(p);
1499         }
1500               
1501       // Define algorithm
1502     fastjet::JetAlgorithm algo = fastjet::antikt_algorithm;
1503     fastjet::RecombinationScheme scheme = fastjet::E_scheme;
1504       // Jet definition
1505     fastjet::JetDefinition jet_def(algo, R_value, scheme);
1506     fastjet::GhostedAreaSpec ghost_spec(ghostMaxRap, numGhostRepeat, ghostArea);
1507     fastjet::AreaType atype = fastjet::active_area;
1508     fastjet::AreaDefinition area_def(atype, ghost_spec);
1509                         
1510       // Clear all vectors for the new bkg D0
1511       RecoJet_pt.clear();
1512       RecoJet_eta.clear();
1513       RecoJet_rapidity.clear();
1514       RecoJet_phi.clear();
1515       RecoJet_E.clear();
1516       RecoJet_M.clear();
1517       RecoJet_rapi.clear();
1518       RecoJet_hasElectron.clear();
1519       RecoJet_maxPtPart_pt.clear();
1520       RecoJet_constituent_pt.clear(); 
1521       RecoJet_constituent_eta.clear();
1522       RecoJet_constituent_phi.clear(); 
1523       RecoJet_constituent_nhits.clear();
1524       RecoJet_constituent_pdgid.clear();
1525       RecoJet_constituent_pdgidTruth.clear();
1526       RecoJet_constituent_idx.clear();
1527           RecoJet_constituent_energy.clear();
1528       RecoJet_constituent_rapi.clear();
1529       
1530       EVETMULTRECO = (int)rcMomPx.GetSize();  
1531             
1532       //Reco Jet clustering 
1533       fastjet::ClusterSequenceArea cs_reco(particles_reco, jet_def, area_def);
1534       std::vector<fastjet::PseudoJet> jets_reco = fastjet::sorted_by_pt(cs_reco.inclusive_jets(minJetPt));
1535       for (auto &jet : jets_reco) {
1536         RecoJet_pt.push_back(jet.pt());
1537         RecoJet_eta.push_back(jet.eta());
1538         RecoJet_rapidity.push_back(jet.rapidity());
1539         RecoJet_phi.push_back(jet.phi_std());
1540         RecoJet_E.push_back(jet.e());
1541         RecoJet_M.push_back(jet.m());
1542         RecoJet_rapi.push_back(jet.rap());  
1543 
1544         bool hasElectron = false;
1545         float maxPtReco = -1.0;
1546         std::vector<float> cpt, ceta, cphi, cenergy,crapi;
1547         std::vector<int> chits, cpdgid, cpdgidtruth,cindex;
1548         cpt.clear(); ceta.clear(); cphi.clear(); chits.clear(); cpdgid.clear(); cpdgidtruth.clear(), cindex.clear(), cenergy.clear(),crapi.clear();
1549 
1550         fastjet::PseudoJet D0_parent;
1551 
1552         for (auto &c : jet.constituents()) {
1553           int idx = c.user_index();
1554           
1555         //Look for the bkg D0 particle of the Jet 
1556           bool is_D0_parent = false;
1557           if(idx == -99999){ // check for bkg D0 particle
1558         for(const auto& p : particles_reco ){ // loop over all the particles of the jet
1559           if (p.user_index() == idx) {
1560             
1561             D0_parent = p;
1562             is_D0_parent = true;
1563             break;
1564           }
1565         }
1566                       
1567         cpt.push_back(D0_parent.pt());
1568         ceta.push_back(D0_parent.eta());
1569         cphi.push_back(D0_parent.phi());
1570         cenergy.push_back(D0_parent.E());
1571         chits.push_back(0);           // 0 hits for neutral D0 
1572         cpdgid.push_back(421);        // PID = 421 for D0
1573         cpdgidtruth.push_back(421);
1574         cindex.push_back(idx);
1575         crapi.push_back(D0_parent.rap());
1576                                                                   // Here, hasElectron = false
1577         if (D0_parent.pt() > maxPtReco) maxPtReco = D0_parent.pt();   
1578           }
1579           if (is_D0_parent) continue;// skip the rest of the loop for bkg D0
1580             
1581           TVector3 v3(rcMomPx[idx], rcMomPy[idx], rcMomPz[idx]);
1582           float rcEnergy = TrkRecoE[idx];
1583           TLorentzVector lv(rcMomPx[idx],  rcMomPy[idx], rcMomPz[idx], rcEnergy);
1584 
1585           cpt.push_back(v3.Pt());
1586           ceta.push_back(v3.Eta());
1587           cphi.push_back(v3.Phi());
1588           chits.push_back(TrkRecoNhits[idx]);
1589           cindex.push_back(idx);
1590           cenergy.push_back(rcEnergy);
1591           crapi.push_back(lv.Rapidity());
1592           if (v3.Pt() > maxPtReco) maxPtReco = v3.Pt();
1593           
1594           // Find electron
1595           int chargePartIndex = idx; 
1596           int elecIndex = -1;
1597           float elecIndexWeight = -1.0;
1598           for(unsigned int itrkass = 0; itrkass < assocWeight.GetSize(); itrkass++){ // Loop Over All ReconstructedChargedParticleAssociations
1599 
1600         if( assocWeight.GetSize() > 0 ){
1601           if(TrkPartAssocRec[itrkass] == chargePartIndex){ // Select Entry Matching the ReconstructedChargedParticle Index
1602             if(assocWeight[itrkass] > elecIndexWeight){ // Find Particle with Greatest Weight = Contributed Most Hits to Track
1603               elecIndex = TrkPartAssocSim[itrkass]; // Get Index of MCParticle Associated with ReconstructedChargedParticle
1604               elecIndexWeight = assocWeight[itrkass];
1605             }
1606           }
1607         }
1608           }
1609 
1610           if(mcPartPdg[elecIndex] == 11){
1611                 hasElectron = true;
1612               }
1613 
1614           cpdgid.push_back(rcPdg[idx]);
1615           cpdgidtruth.push_back(mcPartPdg[elecIndex]);
1616         }
1617 
1618         RecoJet_constituent_pt.push_back(cpt);  
1619         RecoJet_constituent_eta.push_back(ceta);
1620         RecoJet_constituent_phi.push_back(cphi);
1621         RecoJet_constituent_nhits.push_back(chits);
1622         RecoJet_constituent_pdgid.push_back(cpdgid);
1623         RecoJet_constituent_pdgidTruth.push_back(cpdgidtruth);
1624         RecoJet_hasElectron.push_back(hasElectron);
1625         RecoJet_maxPtPart_pt.push_back(maxPtReco);
1626         RecoJet_constituent_idx.push_back(cindex);
1627             RecoJet_constituent_energy.push_back(cenergy);
1628         RecoJet_constituent_rapi.push_back(crapi);
1629 
1630       } 
1631  
1632 
1633 //Reading Reco Jet(BKG):
1634     for(unsigned int i=0; i<RecoJet_pt.size(); i++){  //loop over reco Jets  
1635       
1636       int  jetParticles = RecoJet_constituent_idx[i].size();  // Total no. of particles of the ith Jet
1637      
1638       float jetPt  = RecoJet_pt[i];
1639       float jetEta = RecoJet_eta[i];
1640       float jetRapidity = RecoJet_rapi[i];
1641       float jetPhi = RecoJet_phi[i];
1642       float jetE = RecoJet_E[i];
1643      
1644       TLorentzVector LvJet;
1645       LvJet.SetPtEtaPhiE(jetPt,jetEta, jetPhi, jetE); // Jet TLorentzVector
1646       TVector3 jetMom = LvJet.Vect(); // Jet Momentum Vector
1647           bool c2= false;
1648      
1649       for(unsigned int j=0; j<RecoJet_constituent_idx[i].size(); j++){   //Loop over all particles of the ith Jet
1650 
1651         int particle_idx = RecoJet_constituent_idx[i][j]; // index of jth particle of ith jet
1652         bool is_D0 = false;
1653             
1654         //check if particle is the bkg D0 
1655         if(particle_idx == -99999){        
1656           is_D0 = true;
1657           c2 = true;
1658         }
1659      
1660         if(is_D0){  // if Jet constituent is the bkg D0
1661 
1662           float D0Pt  = RecoJet_constituent_pt[i][j];
1663           float D0Eta = RecoJet_constituent_eta[i][j];
1664           float D0Phi = RecoJet_constituent_phi[i][j];
1665           float D0E = RecoJet_constituent_energy[i][j];
1666           
1667           TLorentzVector lvD0;
1668           lvD0.SetPtEtaPhiE(D0Pt, D0Eta, D0Phi, D0E);
1669           
1670           TVector3 D0_Mom;
1671           D0_Mom.SetPtEtaPhi(D0Pt, D0Eta, D0Phi);
1672               float D0Mass = lvD0.M();
1673           float D0Rapidity = lvD0.Rapidity();
1674     
1675           float bkg_z = float(jetMom.Dot(D0_Mom))/float(jetMom.Dot(jetMom)) ;  // Fragmentation Variable z for bkg D0
1676 
1677               z_bkg =bkg_z;
1678               h3bkg_z->Fill(bkg_z,D0Rapidity,D0Mass);
1679                             
1680               float dPhi = jetPhi - D0Phi;
1681               dPhi = TVector2::Phi_mpi_pi(dPhi);
1682               float dEta = jetEta - D0Eta;
1683               float deltaR = TMath::Sqrt(dEta*dEta + dPhi*dPhi); // bkg_D0 distance from JetAxis in eta-Phi space
1684               
1685               float angleRad = jetMom.Angle(D0_Mom); // angle bw bkg_D0 Momentum vector and jetMomentum vector
1686               float angleDeg = angleRad * 180.0 / TMath::Pi();
1687         
1688               angle_bkg = angleDeg;
1689               dr_bkg = deltaR;
1690               pTjet_bkg = jetPt;
1691               etajet_bkg= jetEta;
1692           D0eta_bkg = D0Eta;
1693 
1694           break; // leave the jet constituent loop once we find the only bkg D0 present
1695         }    
1696       }
1697       if(c2) break; // leave the jet constituent loop once we find the only bkg D0 present 
1698     }
1699     
1700        tree_bkg->Fill();
1701             }
1702 
1703           if(dcaToVtx.Perp() >= 0.02 && dcaToVtx2.Perp() >= 0.02 &&
1704              dcaDaughters < 0.07 && cosTheta > 0.95 && decayLength > 0.05 && V0DcaToVtx < 0.1)
1705             {
1706               if(is_D0_pik)
1707             {
1708               h3InvMass[0][1]->Fill(parent.Pt(), parent.Rapidity(), parent.M());
1709             }
1710               else
1711             {
1712               h3InvMass[1][1]->Fill(parent.Pt(), parent.Rapidity(), parent.M());
1713             }
1714             }
1715         }// if(unlike charge rc pi and k) 
1716         }//for(rc k)
1717     next_pi:
1718       continue;
1719     }//for(rc pi)
1720 
1721       //--------------------------------------------------GenJetClustering--------------------------------------------------------------------------
1722 
1723 
1724       vector<unsigned int> mc_pi_index;
1725       vector<unsigned int> mc_k_index;
1726       mc_pi_index.clear();
1727       mc_k_index.clear();
1728       
1729       // get all MC pi and K
1730       for(unsigned int mc_index=0; mc_index<mcPartMass.GetSize(); mc_index++)
1731     {     
1732       if(abs(mcPartPdg[mc_index]) == 211) mc_pi_index.push_back(mc_index);
1733       if(abs(mcPartPdg[mc_index]) == 321) mc_k_index.push_back(mc_index);
1734       
1735     }
1736           
1737       // pair pion and kaon
1738       for(unsigned int i=0; i<mc_pi_index.size(); i++)//loop over all mc pi
1739     {
1740       for(unsigned int j=0; j<mc_k_index.size(); j++)// loop over all mc k
1741         {
1742           if(mcPartCharge[mc_pi_index[i]]*mcPartCharge[mc_k_index[j]]<0)// only unlike charged mc pi and mc k
1743         {
1744           bool is_D0_pik = false;
1745           for(unsigned int k=0; k<mc_index_D0_pi.size(); k++)
1746             {
1747               if(mc_pi_index[i]==mc_index_D0_pi[k] && mc_k_index[j]==mc_index_D0_k[k])// checking for if mc pi-k pair coming from a D0 decay
1748             {
1749               is_D0_pik = true;
1750               break;
1751             }
1752                 }
1753 
1754           // Storing mc index of pi and k of the pair
1755           int pii_index = mc_pi_index[i]; // mc index of pion
1756           int kk_index =  mc_k_index[j];  // mc index of kaon
1757 
1758               //Signal mc pi-k pair
1759           if(is_D0_pik) 
1760             {  
1761               
1762               // PseudoJet Building
1763               std::vector<fastjet::PseudoJet> particles_gen;       
1764               for (unsigned int i = 0; i < mcPartMass.GetSize(); ++i) {
1765             
1766             TVector3 mom(mcMomPx[i], mcMomPy[i], mcMomPz[i]);
1767             if ( mom.Pt() < minCstPt || mom.Pt() > maxCstPt ) continue;
1768             if ( mcPartCharge[i] == 0 ) continue;
1769             if ( removeelectrons == 1 && mcPartPdg[i] == 11 ) continue;
1770             if ( removeelectrons == 2 && i == ScatteredEGenId ) continue;
1771             
1772             if(i == pii_index || i == kk_index) continue; //skipping the pi-k pair in pseudojet making                                      
1773             
1774             float E =sqrt(mcMomPx[i]*mcMomPx[i] + mcMomPy[i]*mcMomPy[i] + mcMomPz[i]*mcMomPz[i] + mcPartMass[i]*mcPartMass[i]);//energyof mc
1775             fastjet::PseudoJet p(mcMomPx[i], mcMomPy[i], mcMomPz[i], E);   
1776             p.set_user_index(i);
1777             particles_gen.push_back(p);
1778               }
1779           
1780               TLorentzVector lv_pii, lv_kk; //making TLorentzVector of the pair pion and kaon
1781               lv_pii.SetXYZM(mcMomPx[pii_index], mcMomPy[pii_index], mcMomPz[pii_index], mcPartMass[pii_index]); 
1782               lv_kk.SetXYZM(mcMomPx[kk_index], mcMomPy[kk_index], mcMomPz[kk_index], mcPartMass[kk_index]);
1783 
1784               // calculating lorentz vector of the parent D0 of the pi-k pair
1785               TLorentzVector lv_D0 = lv_pii + lv_kk;    
1786     
1787               if ( lv_D0.Pt() < minCstPt || lv_D0.Pt() > maxCstPt ){// if Pt of parent D0 is not in require Pt limits---> Introducing pi and k of the pair to pseudoJet making
1788             
1789             for (int l=0;l<2;l++){
1790               int idx;
1791               
1792               if (l==0){
1793                 idx = pii_index;
1794               }
1795               
1796               if(l==1){
1797                 idx = kk_index;
1798               }
1799 
1800               TVector3 mom(mcMomPx[idx], mcMomPy[idx], mcMomPz[idx]);
1801               if ( mom.Pt() < minCstPt || mom.Pt() > maxCstPt ) continue;
1802 
1803               float E = sqrt(mcMomPx[idx]*mcMomPx[idx] + mcMomPy[idx]*mcMomPy[idx] + mcMomPz[idx]*mcMomPz[idx] + mcPartMass[idx]*mcPartMass[idx]);
1804        
1805               fastjet::PseudoJet p(mcMomPx[idx], mcMomPy[idx], mcMomPz[idx], E);     
1806               p.set_user_index(idx);
1807               particles_gen.push_back(p);
1808               
1809             }
1810               }
1811         
1812                                                                                                     
1813               if ( lv_D0.Pt() > minCstPt && lv_D0.Pt() < maxCstPt ){// introducing parent D0 to pseudojets
1814             
1815             fastjet::PseudoJet p(lv_D0.Px(), lv_D0.Py(),lv_D0.Pz(), lv_D0.E());
1816             p.set_user_index(-99999);// setting signal parent D0 index = negative integer (-99999)                                          
1817             particles_gen.push_back(p);
1818               }
1819         
1820                             
1821               // Define algorithm
1822               fastjet::JetAlgorithm algo = fastjet::antikt_algorithm;
1823               fastjet::RecombinationScheme scheme = fastjet::E_scheme;
1824               // Jet definition
1825               fastjet::JetDefinition jet_def(algo, R_value, scheme);
1826               fastjet::GhostedAreaSpec ghost_spec(ghostMaxRap, numGhostRepeat, ghostArea);
1827               fastjet::AreaType atype = fastjet::active_area;
1828               fastjet::AreaDefinition area_def(atype, ghost_spec);      
1829 
1830               //clearing the vectors
1831               GenJet_pt.clear();
1832               GenJet_eta.clear();
1833               GenJet_phi.clear();
1834               GenJet_E.clear();
1835               GenJet_M.clear();
1836               GenJet_rapi.clear();
1837               GenJet_hasElectron.clear();
1838               GenJet_hasNeutral.clear();
1839               GenJet_maxPtPart_pt.clear();
1840               GenJet_constituent_pt.clear(); 
1841               GenJet_constituent_eta.clear();
1842               GenJet_constituent_phi.clear(); 
1843               GenJet_constituent_pdgid.clear();
1844               GenJet_constituent_rapi.clear();
1845               GenJet_constituent_idx.clear();
1846               GenJet_constituent_energy.clear();
1847               
1848               EVETMULTGEN = (int)TrkGenPx.GetSize();
1849               
1850               //Gen clustering 
1851               fastjet:: ClusterSequenceArea cs_gen(particles_gen, jet_def, area_def);
1852               std::vector<fastjet::PseudoJet> jets_gen = fastjet::sorted_by_pt(cs_gen.inclusive_jets(minJetPt));
1853               for (auto &jet : jets_gen) { //loop over genJets
1854             GenJet_pt.push_back(jet.pt());
1855             GenJet_eta.push_back(jet.eta());
1856             GenJet_phi.push_back(jet.phi_std());
1857             GenJet_E.push_back(jet.e());
1858             GenJet_M.push_back(jet.m());
1859             GenJet_rapi.push_back(jet.rap());
1860 
1861     
1862             bool hasGenElectron = false;
1863             bool hasGenNeutral = false;
1864             float maxPtGen = -1.0;
1865             std::vector<float> gpt, geta, gphi,grapi,genergy;
1866             std::vector<int> gpdgid,gidx;
1867             gpt.clear(); geta.clear(); gphi.clear(); gpdgid.clear();gidx.clear();grapi.clear();genergy.clear();
1868             
1869             fastjet::PseudoJet D0_parent;
1870             
1871             for (auto &c : jet.constituents()) {//loop over genJetConstituents
1872               int idx = c.user_index();
1873               
1874               bool is_D0_parent = false;
1875               if(idx == -99999){     //looking for D0 parent                                                                                        
1876                 for(const auto& p : particles_gen ){// loop over all the pseudojets                                                         
1877 
1878                   if (p.user_index() == idx) {
1879             
1880                 D0_parent = p;
1881                 is_D0_parent = true;
1882                 break;
1883                   }
1884                 }
1885         
1886                 gpt.push_back(D0_parent.pt());
1887                 geta.push_back(D0_parent.eta());
1888                 gphi.push_back(D0_parent.phi());
1889                 genergy.push_back(D0_parent.E());                                             
1890                 gpdgid.push_back(421);                                                                                   
1891                 gidx.push_back(idx);
1892                 grapi.push_back(D0_parent.rap());
1893                 
1894                 if (D0_parent.pt() > maxPtGen) maxPtGen = D0_parent.pt();
1895                 hasGenElectron = false;
1896                 hasGenNeutral = true;
1897               }
1898               
1899               if (is_D0_parent) continue;// skip the rest of the loop for the parent D0
1900       
1901       
1902               TVector3 gv(mcMomPx[idx], mcMomPy[idx], mcMomPz[idx]);
1903               TLorentzVector glv;
1904               glv.SetXYZM(mcMomPx[idx], mcMomPy[idx], mcMomPz[idx], mcPartMass[idx]);   
1905 
1906               gpt.push_back(gv.Pt());
1907               geta.push_back(gv.Eta());
1908               gphi.push_back(gv.Phi());
1909               gpdgid.push_back(mcPartPdg[idx]);
1910               gidx.push_back(idx);
1911               grapi.push_back(glv.Rapidity());
1912               genergy.push_back(glv.E());
1913               
1914               if (gv.Pt() > maxPtGen) maxPtGen = gv.Pt();
1915               
1916               if (mcPartPdg[idx] == 11) hasGenElectron = true;
1917               if (mcPartCharge[idx] == 0) hasGenNeutral = true;
1918             }
1919             
1920             GenJet_constituent_pt.push_back(gpt);
1921             GenJet_constituent_eta.push_back(geta);
1922             GenJet_constituent_phi.push_back(gphi);
1923             GenJet_hasElectron.push_back(hasGenElectron);
1924             GenJet_hasNeutral.push_back(hasGenNeutral);
1925             GenJet_maxPtPart_pt.push_back(maxPtGen);
1926             GenJet_constituent_pdgid.push_back(gpdgid);
1927             GenJet_constituent_rapi.push_back(grapi);
1928             GenJet_constituent_idx.push_back(gidx);
1929             GenJet_constituent_energy.push_back(genergy);
1930               }
1931               
1932               //Gen Jet Reading
1933 
1934               for(int i=0;i<GenJet_eta.size();i++){//loop over all the genJets
1935 
1936             int  jetmulti = GenJet_constituent_idx[i].size();  // Total no. of particles of the ith genJet
1937      
1938             float jetPt  = GenJet_pt[i];
1939             float jetEta = GenJet_eta[i];
1940             float jetRapi = GenJet_rapi[i];
1941             float jetPhi = GenJet_phi[i];
1942             float jetE = GenJet_E[i];
1943             
1944             TLorentzVector lvJet;
1945             lvJet.SetPtEtaPhiE(jetPt,jetEta, jetPhi, jetE); // Jet TLorentzVector     
1946             TVector3 jetMom = lvJet.Vect(); // Jet Momentum Vector
1947             bool c2 = false;
1948             
1949             for(unsigned int j=0; j<GenJet_constituent_idx[i].size(); j++){   //Loop over all particles of the ith Jet
1950 
1951               int genIdx = GenJet_constituent_idx[i][j]; // index of jth particle of ith jet
1952               bool isD0 = false;
1953              
1954               
1955               //check if jet constituent is the D0 parent
1956               if( genIdx == -99999 ){        
1957                 isD0 = true;
1958                 c2 = true;
1959               }
1960      
1961               if(isD0){  // if Jet constituent is the D0 parent
1962 
1963                 float D0Pt  = GenJet_constituent_pt[i][j];
1964                 float D0Eta = GenJet_constituent_eta[i][j];
1965                 float D0Phi = GenJet_constituent_phi[i][j];
1966                 float D0E = GenJet_constituent_energy[i][j];
1967                 float D0Rapi = GenJet_constituent_rapi[i][j];
1968           
1969                 TLorentzVector lvD0;
1970                 lvD0.SetPtEtaPhiE(D0Pt, D0Eta, D0Phi, D0E);
1971                 
1972                 TVector3 D0_Mom;
1973                 D0_Mom.SetPtEtaPhi(D0Pt, D0Eta, D0Phi);
1974                 
1975                 float gz = float(jetMom.Dot(D0_Mom))/float(jetMom.Dot(jetMom)) ;  // Fragmentation Variable z 
1976                 cout<<"genZ = "<<gz<<endl;
1977                             
1978                 float dPhi = jetPhi - D0Phi;
1979                 dPhi = TVector2::Phi_mpi_pi(dPhi);
1980                 
1981                 float dEta = jetEta - D0Eta;
1982                 float drapi = jetRapi - D0Rapi;
1983           
1984                 float deltaR_eta = TMath::Sqrt(dEta*dEta + dPhi*dPhi); // D0 distance from JetAxis in eta-Phi space
1985                 float deltaR_y = TMath::Sqrt(drapi*drapi + dPhi*dPhi); // D0 distance from JetAxis in y-Phi space 
1986           
1987                 float angleRad = jetMom.Angle(D0_Mom); // angle bw parent D0 Momentum vector and jetMomentum vector
1988                 float angleDeg = angleRad * 180.0 / TMath::Pi();
1989 
1990                 //storing D0 parameters to branches of gen signal tree
1991                 genZ_sig = gz;
1992                 gdr_eta_sig = deltaR_eta;
1993                 gdr_y_sig = deltaR_y;
1994                 gD0jetAngle_sig = angleDeg;
1995                 pt_gjet_sig = jetPt;
1996                 eta_gjet_sig = jetEta;
1997                 y_gjet_sig = jetRapi;
1998           
1999                 break;// break the jetConstituent loop immediately the D0 parent is found
2000               }
2001             }
2002             if(c2) break;// break the jet loop immediately when the D0 parent is found
2003               }
2004               tree_gen_sig->Fill(); // tree filled only fill once, in the that jet case which have D0 in it
2005             }
2006 
2007 
2008           //gen BKG D0
2009           if(!is_D0_pik) 
2010             {  
2011 
2012                           
2013               // PseudoJet Building
2014               std::vector<fastjet::PseudoJet> particles_gen;       
2015               for (unsigned int i = 0; i < mcPartMass.GetSize(); ++i) {
2016             
2017             TVector3 mom(mcMomPx[i], mcMomPy[i], mcMomPz[i]);
2018             if ( mom.Pt() < minCstPt || mom.Pt() > maxCstPt ) continue;
2019             if ( mcPartCharge[i] == 0 ) continue;
2020             if ( removeelectrons == 1 && mcPartPdg[i] == 11 ) continue;
2021             if ( removeelectrons == 2 && i == ScatteredEGenId ) continue;
2022             
2023             if(i == pii_index || i == kk_index) continue; //skipping the pi-k pair in pseudojet making                                      
2024             
2025             float E =sqrt(mcMomPx[i]*mcMomPx[i] + mcMomPy[i]*mcMomPy[i] + mcMomPz[i]*mcMomPz[i] + mcPartMass[i]*mcPartMass[i]);//energyof mc
2026             fastjet::PseudoJet p(mcMomPx[i], mcMomPy[i], mcMomPz[i], E);   
2027             p.set_user_index(i);
2028             particles_gen.push_back(p);
2029               }
2030           
2031               TLorentzVector lv_pii, lv_kk; //making TLorentzVector of the pair pion and kaon
2032               lv_pii.SetXYZM(mcMomPx[pii_index], mcMomPy[pii_index], mcMomPz[pii_index], mcPartMass[pii_index]); 
2033               lv_kk.SetXYZM(mcMomPx[kk_index], mcMomPy[kk_index], mcMomPz[kk_index], mcPartMass[kk_index]);
2034 
2035               // calculating lorentz vector of the parent D0 of the pi-k pair
2036               TLorentzVector lv_D0 = lv_pii + lv_kk;    
2037     
2038               if ( lv_D0.Pt() < minCstPt || lv_D0.Pt() > maxCstPt ){// if Pt of parent D0 is not in require Pt limits---> Introducing pi and k of the pair to pseudoJet making
2039             
2040             for (int l=0;l<2;l++){
2041               int idx;
2042               
2043               if (l==0){
2044                 idx = pii_index;
2045               }
2046               
2047               if(l==1){
2048                 idx = kk_index;
2049               }
2050 
2051               TVector3 mom(mcMomPx[idx], mcMomPy[idx], mcMomPz[idx]);
2052               if ( mom.Pt() < minCstPt || mom.Pt() > maxCstPt ) continue;
2053 
2054               float E = sqrt(mcMomPx[idx]*mcMomPx[idx] + mcMomPy[idx]*mcMomPy[idx] + mcMomPz[idx]*mcMomPz[idx] + mcPartMass[idx]*mcPartMass[idx]);
2055        
2056               fastjet::PseudoJet p(mcMomPx[idx], mcMomPy[idx], mcMomPz[idx], E);     
2057               p.set_user_index(idx);
2058               particles_gen.push_back(p);
2059               
2060             }
2061               }
2062         
2063                                                                                                     
2064               if ( lv_D0.Pt() > minCstPt && lv_D0.Pt() < maxCstPt ){// introducing parent D0 to pseudojets
2065             
2066             fastjet::PseudoJet p(lv_D0.Px(), lv_D0.Py(),lv_D0.Pz(), lv_D0.E());
2067             p.set_user_index(-99999);// setting bkg  parent D0 index = negative integer (-99999)                                          
2068             particles_gen.push_back(p);
2069               }
2070         
2071                             
2072               // Define algorithm
2073               fastjet::JetAlgorithm algo = fastjet::antikt_algorithm;
2074               fastjet::RecombinationScheme scheme = fastjet::E_scheme;
2075               // Jet definition
2076               fastjet::JetDefinition jet_def(algo, R_value, scheme);
2077               fastjet::GhostedAreaSpec ghost_spec(ghostMaxRap, numGhostRepeat, ghostArea);
2078               fastjet::AreaType atype = fastjet::active_area;
2079               fastjet::AreaDefinition area_def(atype, ghost_spec);      
2080 
2081               //clearing the vectors
2082               GenJet_pt.clear();
2083               GenJet_eta.clear();
2084               GenJet_phi.clear();
2085               GenJet_E.clear();
2086               GenJet_M.clear();
2087               GenJet_rapi.clear();
2088               GenJet_hasElectron.clear();
2089               GenJet_hasNeutral.clear();
2090               GenJet_maxPtPart_pt.clear();
2091               GenJet_constituent_pt.clear(); 
2092               GenJet_constituent_eta.clear();
2093               GenJet_constituent_phi.clear(); 
2094               GenJet_constituent_pdgid.clear();
2095               GenJet_constituent_rapi.clear();
2096               GenJet_constituent_idx.clear();
2097               GenJet_constituent_energy.clear();
2098               
2099               EVETMULTGEN = (int)TrkGenPx.GetSize();
2100               
2101               //Gen clustering 
2102               fastjet:: ClusterSequenceArea cs_gen(particles_gen, jet_def, area_def);
2103               std::vector<fastjet::PseudoJet> jets_gen = fastjet::sorted_by_pt(cs_gen.inclusive_jets(minJetPt));
2104               for (auto &jet : jets_gen) { //loop over genJets
2105             GenJet_pt.push_back(jet.pt());
2106             GenJet_eta.push_back(jet.eta());
2107             GenJet_phi.push_back(jet.phi_std());
2108             GenJet_E.push_back(jet.e());
2109             GenJet_M.push_back(jet.m());
2110             GenJet_rapi.push_back(jet.rap());
2111 
2112     
2113             bool hasGenElectron = false;
2114             bool hasGenNeutral = false;
2115             float maxPtGen = -1.0;
2116             std::vector<float> gpt, geta, gphi,grapi,genergy;
2117             std::vector<int> gpdgid,gidx;
2118             gpt.clear(); geta.clear(); gphi.clear(); gpdgid.clear();gidx.clear();grapi.clear();genergy.clear();
2119             
2120             fastjet::PseudoJet D0_parent;
2121             
2122             for (auto &c : jet.constituents()) {//loop over genJetConstituents
2123               int idx = c.user_index();
2124               
2125               bool is_D0_parent = false;
2126               if(idx == -99999){     //looking for D0 parent                                                                                        
2127                 for(const auto& p : particles_gen ){// loop over all the pseudojets                                                         
2128 
2129                   if (p.user_index() == idx) {
2130             
2131                 D0_parent = p;
2132                 is_D0_parent = true;
2133                 break;
2134                   }
2135                 }
2136         
2137                 gpt.push_back(D0_parent.pt());
2138                 geta.push_back(D0_parent.eta());
2139                 gphi.push_back(D0_parent.phi());
2140                 genergy.push_back(D0_parent.E());                                             
2141                 gpdgid.push_back(421);                                                                                   
2142                 gidx.push_back(idx);
2143                 grapi.push_back(D0_parent.rap());
2144                 
2145                 if (D0_parent.pt() > maxPtGen) maxPtGen = D0_parent.pt();
2146                 hasGenElectron = false;
2147                 hasGenNeutral = true;
2148               }
2149               
2150               if (is_D0_parent) continue;// skip the rest of the loop for the parent D0
2151       
2152       
2153               TVector3 gv(mcMomPx[idx], mcMomPy[idx], mcMomPz[idx]);
2154               TLorentzVector glv;
2155               glv.SetXYZM(mcMomPx[idx], mcMomPy[idx], mcMomPz[idx], mcPartMass[idx]);   
2156 
2157               gpt.push_back(gv.Pt());
2158               geta.push_back(gv.Eta());
2159               gphi.push_back(gv.Phi());
2160               gpdgid.push_back(mcPartPdg[idx]);
2161               gidx.push_back(idx);
2162               grapi.push_back(glv.Rapidity());
2163               genergy.push_back(glv.E());
2164               
2165               if (gv.Pt() > maxPtGen) maxPtGen = gv.Pt();
2166               
2167               if (mcPartPdg[idx] == 11) hasGenElectron = true;
2168               if (mcPartCharge[idx] == 0) hasGenNeutral = true;
2169             }
2170             
2171             GenJet_constituent_pt.push_back(gpt);
2172             GenJet_constituent_eta.push_back(geta);
2173             GenJet_constituent_phi.push_back(gphi);
2174             GenJet_hasElectron.push_back(hasGenElectron);
2175             GenJet_hasNeutral.push_back(hasGenNeutral);
2176             GenJet_maxPtPart_pt.push_back(maxPtGen);
2177             GenJet_constituent_pdgid.push_back(gpdgid);
2178             GenJet_constituent_rapi.push_back(grapi);
2179             GenJet_constituent_idx.push_back(gidx);
2180             GenJet_constituent_energy.push_back(genergy);
2181               }
2182               
2183               //Gen Bkg Jet Reading
2184 
2185               for(int i=0;i<GenJet_eta.size();i++){//loop over all the genJets
2186 
2187             int  jetmulti = GenJet_constituent_idx[i].size();  // Total no. of particles of the ith genJet
2188      
2189             float jetPt  = GenJet_pt[i];
2190             float jetEta = GenJet_eta[i];
2191             float jetRapi = GenJet_rapi[i];
2192             float jetPhi = GenJet_phi[i];
2193             float jetE = GenJet_E[i];
2194             
2195             TLorentzVector lvJet;
2196             lvJet.SetPtEtaPhiE(jetPt,jetEta, jetPhi, jetE); // Jet TLorentzVector     
2197             TVector3 jetMom = lvJet.Vect(); // Jet Momentum Vector
2198             bool c2 = false;
2199             
2200             for(unsigned int j=0; j<GenJet_constituent_idx[i].size(); j++){   //Loop over all particles of the ith Jet
2201 
2202               int genIdx = GenJet_constituent_idx[i][j]; // index of jth particle of ith jet
2203               bool isD0 = false;
2204              
2205               
2206               //check if jet constituent is the D0 parent
2207               if( genIdx == -99999 ){        
2208                 isD0 = true;
2209                 c2 = true;
2210               }
2211      
2212               if(isD0){  // if Jet constituent is the D0 parent
2213 
2214                 float D0Pt  = GenJet_constituent_pt[i][j];
2215                 float D0Eta = GenJet_constituent_eta[i][j];
2216                 float D0Phi = GenJet_constituent_phi[i][j];
2217                 float D0E = GenJet_constituent_energy[i][j];
2218                 float D0Rapi = GenJet_constituent_rapi[i][j];
2219           
2220                 TLorentzVector lvD0;
2221                 lvD0.SetPtEtaPhiE(D0Pt, D0Eta, D0Phi, D0E);
2222                 
2223                 TVector3 D0_Mom;
2224                 D0_Mom.SetPtEtaPhi(D0Pt, D0Eta, D0Phi);
2225                 
2226                 float gz = float(jetMom.Dot(D0_Mom))/float(jetMom.Dot(jetMom)) ;  // Fragmentation Variable z 
2227                 cout<<"genZ = "<<gz<<endl;
2228                             
2229                 float dPhi = jetPhi - D0Phi;
2230                 dPhi = TVector2::Phi_mpi_pi(dPhi);
2231                 
2232                 float dEta = jetEta - D0Eta;
2233                 float drapi = jetRapi - D0Rapi;
2234           
2235                 float deltaR_eta = TMath::Sqrt(dEta*dEta + dPhi*dPhi); // D0 distance from JetAxis in eta-Phi space
2236                 float deltaR_y = TMath::Sqrt(drapi*drapi + dPhi*dPhi); // D0 distance from JetAxis in y-Phi space 
2237           
2238                 float angleRad = jetMom.Angle(D0_Mom); // angle bw parent D0 Momentum vector and jetMomentum vector
2239                 float angleDeg = angleRad * 180.0 / TMath::Pi();
2240 
2241                 //storing D0 parameters to branches of gen signal tree
2242                 genZ_bkg = gz;
2243                 gdr_eta_bkg = deltaR_eta;
2244                 gdr_y_bkg = deltaR_y;
2245                 gD0jetAngle_bkg = angleDeg;
2246                 pt_gjet_bkg = jetPt;
2247                 eta_gjet_bkg = jetEta;
2248                 y_gjet_bkg = jetRapi;
2249           
2250                 break;// break the jetConstituent loop immediately the D0 parent is found
2251               }
2252             }
2253             if(c2) break;// break the jet loop immediately when the D0 parent is found
2254               }
2255               tree_gen_bkg->Fill(); // tree filled only fill once, in the that jet case which have D0 in it
2256             }
2257 
2258         }
2259         }
2260     }
2261 
2262       
2263       //--------------------------------------------------GenJetClustering-----------------------------------------------------------------------------
2264       cout<<nevents<<" Completed !!"<<endl;
2265       nevents++;
2266     } 
2267 
2268   file_signal->cd();  
2269   tree_sig->Write();
2270   file_signal->Close();  
2271   
2272   file_bkg->cd();  
2273   tree_bkg->Write();
2274   file_bkg->Close();
2275 
2276   file_gen->cd();
2277   tree_gen_sig->Write();
2278   tree_gen_bkg->Write();
2279   file_gen->Close();
2280   
2281   fout_mcgen->cd();
2282   tree_D0->Write();
2283   fout_mcgen->Close();
2284   
2285   TFile *outfile = new TFile(outname.Data(), "recreate");
2286 
2287   hEventStat->Write();
2288   hMcMult->Write();
2289   hMcVtxX->Write();
2290   hMcVtxY->Write();
2291   hMcVtxZ->Write();
2292   
2293   hPullVtxX->Write();
2294   hPullVtxY->Write();  
2295   hPullVtxZ->Write();
2296   
2297   hRes_SVx_Helixfit->Write();
2298   hRes_SVy_Helixfit->Write();
2299   hRes_SVz_Helixfit->Write(); 
2300   hRes_SVxy_Helixfit->Write();  
2301  
2302   hRes_SVx_Helixana->Write();
2303   hRes_SVy_Helixana->Write();
2304   hRes_SVz_Helixana->Write();
2305   hRes_SVxy_Helixana->Write();    
2306       
2307   hchi2_vtx->Write();
2308   hchi2_vtx_sig->Write();
2309   hchi2_vtx_bkg->Write();  
2310   hRes_SVx_Helixfit_pull->Write();
2311   hRes_SVy_Helixfit_pull->Write(); 
2312   hRes_SVz_Helixfit_pull->Write(); 
2313     
2314   hD0DecayVxVy->Write();
2315   hD0DecayVrVz->Write();
2316   
2317   hMCD0PtRap->Write();
2318 
2319   hMcPiPtEta->Write();
2320   hMcPiPtEtaReco->Write();
2321   hMcKPtEta->Write();
2322   hMcKPtEtaReco->Write();
2323   
2324   hNRecoVtx->Write();
2325   h3sig_z->Write();
2326   h3bkg_z->Write();
2327 
2328   hmc_eta->Write();
2329   hmc_eta_e->Write();
2330   hmc_eta_pi->Write();
2331   hmc_eta_k->Write();
2332   hmc_eta_p->Write();
2333 
2334   hreco_eta->Write();
2335   hreco_eta_e->Write();
2336   hreco_eta_pi->Write();
2337   hreco_eta_k->Write();
2338   hreco_eta_p->Write();
2339 
2340   
2341 
2342   for(int ip=0; ip<2; ip++)
2343     {
2344       hRcSecPartLocaToRCVtx[ip]->Write();
2345       hRcSecPartLocbToRCVtx[ip]->Write();
2346       hRcPrimPartLocaToRCVtx[ip]->Write();
2347       hRcPrimPartLocbToRCVtx[ip]->Write();
2348     }
2349 
2350   for(int i=0; i<3; i++)
2351     {
2352       for(int j=0; j<3; j++)
2353     {
2354       hPrimTrkDcaToRCVtx[i][j]->Write();
2355     }
2356     }
2357 
2358   for(int i=0; i<2; i++)
2359     {
2360       h3PairDca12[i]->Write();
2361       h3PairCosTheta[i]->Write();
2362       h3PairDca[i]->Write();
2363       h3PairDecayLength[i]->Write();
2364     }
2365 
2366   for(int i=0; i<2; i++)
2367     {
2368       for(int j=0; j<2; j++)
2369     {
2370       h3InvMass[i][j]->Write();
2371     }
2372     }
2373   
2374   
2375   outfile->Close();
2376 }
2377 
2378 //======================================
2379 TVector3 getDcaToVtx(const int index, TVector3 vtx)
2380 {
2381   TVector3 pos(rcTrkLoca2->At(index) * sin(rcTrkPhi2->At(index)) * -1 * millimeter, rcTrkLoca2->At(index) * cos(rcTrkPhi2->At(index)) * millimeter, rcTrkLocb2->At(index) * millimeter);
2382   TVector3 mom(rcMomPx2->At(index), rcMomPy2->At(index), rcMomPz2->At(index));
2383    
2384   StPhysicalHelix pHelix(mom, pos, bField * tesla, rcCharge2->At(index));
2385 
2386   TVector3 vtx_tmp;
2387   vtx_tmp.SetXYZ(vtx.x()*millimeter, vtx.y()*millimeter, vtx.z()*millimeter);
2388   
2389   pHelix.moveOrigin(pHelix.pathLength(vtx_tmp));
2390   TVector3 dcaToVtx = pHelix.origin() - vtx_tmp;
2391 
2392   dcaToVtx.SetXYZ(dcaToVtx.x()/millimeter, dcaToVtx.y()/millimeter, dcaToVtx.z()/millimeter);
2393   
2394   return dcaToVtx;
2395 }
2396 
2397 //======================================
2398 TLorentzVector getPairParent(const int index1, const int index2, TVector3 vtx,
2399                  float &dcaDaughters, float &cosTheta, float &cosTheta_xy, float &decayLength, float &V0DcaToVtx, float &sigma_vtx, TVector3 &decayVertex, TVector3 &decayVertex_ana, double &chi2_ndf, double * parFitErr)
2400 {
2401   // -- get helix
2402   TVector3 pos1(rcTrkLoca2->At(index1) * sin(rcTrkPhi2->At(index1)) * -1 * millimeter, rcTrkLoca2->At(index1) * cos(rcTrkPhi2->At(index1)) * millimeter, rcTrkLocb2->At(index1) * millimeter);
2403   TVector3 mom1(rcMomPx2->At(index1), rcMomPy2->At(index1), rcMomPz2->At(index1));
2404 
2405   TVector3 pos2(rcTrkLoca2->At(index2) * sin(rcTrkPhi2->At(index2)) * -1 * millimeter, rcTrkLoca2->At(index2) * cos(rcTrkPhi2->At(index2)) * millimeter, rcTrkLocb2->At(index2) * millimeter);
2406   TVector3 mom2(rcMomPx2->At(index2), rcMomPy2->At(index2), rcMomPz2->At(index2));
2407 
2408   float charge1 = rcCharge2->At(index1);
2409   float charge2 = rcCharge2->At(index2);
2410   
2411   StPhysicalHelix p1Helix(mom1, pos1, bField * tesla, charge1);
2412   StPhysicalHelix p2Helix(mom2, pos2, bField * tesla, charge2);
2413 
2414   TVector3 vtx_tmp;
2415   vtx_tmp.SetXYZ(vtx.x()*millimeter, vtx.y()*millimeter, vtx.z()*millimeter);
2416   
2417   double s1, s2;
2418   getDecayVertex_Chi2fit(index1, index2, s1, s2, decayVertex, chi2_ndf, parFitErr);
2419  
2420   TVector3 const p1AtDcaToP2 = p1Helix.at(s1);
2421   TVector3 const p2AtDcaToP1 = p2Helix.at(s2);
2422   // printf("p1AtDcaToP2 origin = (%2.4f, %2.4f, %2.4f)\n", p1AtDcaToP2.x(), p1AtDcaToP2.y(), p1AtDcaToP2.z());
2423   // printf("p2AtDcaToP1 origin = (%2.4f, %2.4f, %2.4f)\n", p2AtDcaToP1.x(), p2AtDcaToP1.y(), p2AtDcaToP1.z());
2424   
2425   // -- calculate DCA of particle1 to particle2 at their DCA
2426   dcaDaughters = (p1AtDcaToP2 - p2AtDcaToP1).Mag()/millimeter;
2427     
2428   // -- calculate Lorentz vector of particle1-particle2 pair
2429   TVector3 const p1MomAtDca = p1Helix.momentumAt(s1,  bField * tesla);
2430   TVector3 const p2MomAtDca = p2Helix.momentumAt(s2, bField * tesla);
2431   
2432   TLorentzVector p1FourMom(p1MomAtDca, sqrt(p1MomAtDca.Mag2()+gPionMass*gPionMass));
2433   TLorentzVector p2FourMom(p2MomAtDca, sqrt(p2MomAtDca.Mag2()+gKaonMass*gKaonMass));
2434   
2435   TLorentzVector parent = p1FourMom + p2FourMom;
2436 
2437   // -- calculate decay vertex (secondary or tertiary)
2438   decayVertex_ana = (p1AtDcaToP2 + p2AtDcaToP1) * 0.5 ;
2439   sigma_vtx = sqrt((p1AtDcaToP2-decayVertex).Mag2()+(p2AtDcaToP1-decayVertex).Mag2())/millimeter;
2440     
2441   // -- calculate pointing angle and decay length with respect to primary vertex
2442   //    if decay vertex is a tertiary vertex
2443   //    -> only rough estimate -> needs to be updated after secondary vertex is found
2444   TVector3 vtxToV0 = decayVertex - vtx_tmp;
2445   TVector3 vtxToV0_xy(vtxToV0.x(), vtxToV0.y(), 0.);
2446   TVector3 parent_xy(parent.Vect().x(),parent.Vect().y(),0.);
2447   float pointingAngle = vtxToV0.Angle(parent.Vect());
2448   float pointingAngle_xy = vtxToV0_xy.Angle(parent_xy);
2449   cosTheta = std::cos(pointingAngle);
2450   cosTheta_xy = std::cos(pointingAngle_xy);  
2451   decayLength = vtxToV0.Mag()/millimeter;
2452 
2453   // -- calculate V0 DCA to primary vertex
2454   V0DcaToVtx = decayLength * std::sin(pointingAngle);
2455     
2456   //TVector3 dcaToVtx = getDcaToVtx(parent.Vect(), decayVertex, 0, vtx);
2457   //V0DcaToVtx = dcaToVtx.Mag();
2458   return parent;
2459 }
2460 
2461 
2462 
2463 void getDecayVertex_Chi2fit(const int index1, const int index2, double &s1, double &s2, TVector3 &vertex, double &chi2_ndf, double *parFitErr)
2464 {
2465     TVector3 pos1(rcTrkLoca2->At(index1) * sin(rcTrkPhi2->At(index1)) * -1 * millimeter,
2466                   rcTrkLoca2->At(index1) * cos(rcTrkPhi2->At(index1)) * millimeter,
2467                   rcTrkLocb2->At(index1) * millimeter);
2468 
2469     TVector3 mom1(rcMomPx2->At(index1), rcMomPy2->At(index1), rcMomPz2->At(index1));
2470 
2471     TVector3 pos2(rcTrkLoca2->At(index2) * sin(rcTrkPhi2->At(index2)) * -1 * millimeter,
2472                   rcTrkLoca2->At(index2) * cos(rcTrkPhi2->At(index2)) * millimeter,
2473                   rcTrkLocb2->At(index2) * millimeter);
2474 
2475     TVector3 mom2(rcMomPx2->At(index2), rcMomPy2->At(index2), rcMomPz2->At(index2));
2476 
2477 
2478     float charge1 = rcCharge2->At(index1);
2479     float charge2 = rcCharge2->At(index2);
2480 
2481 
2482     StPhysicalHelix helix1(mom1, pos1, bField * tesla, charge1);
2483     StPhysicalHelix helix2(mom2, pos2, bField * tesla, charge2);
2484    
2485     std::array<float, 21>& cov_track1 = rcTrkCov->At(index1);
2486     std::array<float, 21>& cov_track2 = rcTrkCov->At(index2);  
2487     
2488     pair<double, double> const ss = helix1.pathLengths(helix2);
2489     TVector3 const p1_init = helix1.at(ss.first);
2490     TVector3 const p2_init = helix2.at(ss.second); 
2491     TVector3 const mid_point = 0.5*(p1_init+p2_init); 
2492 
2493       // Perform Minimization
2494     const Int_t nPar = 5;
2495     Chi2Minimization d2Function(helix1,helix2,cov_track1,cov_track2);
2496     ROOT::Math::Functor fcn(d2Function,nPar); // 5 parameters
2497     ROOT::Fit::Fitter fitter;
2498 
2499     double pStart[nPar] = {mid_point.X(),mid_point.Y(),mid_point.Z(),ss.first,ss.second};
2500     fitter.SetFCN(fcn, pStart,nPar,1);
2501     
2502     fitter.Config().ParSettings(0).SetName("x0");
2503     fitter.Config().ParSettings(0).SetStepSize(0.01);
2504    // fitter.Config().ParSettings(0).SetLimits(-1., 1.);    
2505     // No limits for x, y, z
2506 
2507     fitter.Config().ParSettings(1).SetName("y0");
2508     fitter.Config().ParSettings(1).SetStepSize(0.01);
2509     //fitter.Config().ParSettings(1).SetLimits(-1., 1.);
2510 
2511     fitter.Config().ParSettings(2).SetName("z0");
2512     fitter.Config().ParSettings(2).SetStepSize(0.01);
2513    // fitter.Config().ParSettings(2).SetLimits(-10., 10.);    
2514 
2515     fitter.Config().ParSettings(3).SetName("s1");
2516     fitter.Config().ParSettings(3).SetValue(0.0);
2517     fitter.Config().ParSettings(3).SetStepSize(0.01);
2518     //fitter.Config().ParSettings(3).SetLimits(-1., 1.);
2519 
2520     fitter.Config().ParSettings(4).SetName("s2");
2521     fitter.Config().ParSettings(4).SetValue(0.0);
2522     fitter.Config().ParSettings(4).SetStepSize(0.01);
2523     //fitter.Config().ParSettings(4).SetLimits(-1., 1.);      
2524     
2525     fitter.Config().MinimizerOptions().SetMaxIterations(10000); 
2526     // do the fit 
2527 
2528     Bool_t ok = fitter.FitFCN();
2529     if (!ok) Error("Fitting","Fitting failed");
2530     const ROOT::Fit::FitResult & result = fitter.Result();
2531    // double chi2 = fitter.Result().Chi2();
2532     double ndf = 2*3-nPar;
2533     chi2_ndf = fitter.Result().MinFcnValue()/ndf;  // Minimum value of your function
2534     
2535     int status = fitter.Result().Status();
2536   //  if (status>0 ) {printf("Fit Failed!!!!\n");}
2537    // if (status>0 || chi2_ndf>10. ) return;
2538    // cout <<"\033[1;31m Fit Result Chi2:\033[0m"<<chi2_ndf<<endl;
2539    // result.Print(std::cout);
2540    
2541    // Get the covariance matrix
2542  //  TMatrixDSym covMatrix(5);
2543   // result.GetCovarianceMatrix(covMatrix); // Matrix for the parameter errors
2544   // covMatrix.Print();
2545    
2546    const double * parFit = result.GetParams();
2547    const double *FitErr = result.GetErrors();
2548     
2549    for (int i = 0; i < nPar; ++i) parFitErr[i] = FitErr[i];
2550    vertex.SetXYZ(parFit[0], parFit[1], parFit[2]);
2551    s1 = parFit[3]; s2 = parFit[4];
2552 }
2553
[ 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!