Back to home page

EIC code displayed by LXR
 
 

    

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

0001 #ifdef __CINT__
0002 
0003 #pragma link off all globals;
0004 #pragma link off all classes;
0005 #pragma link off all functions;
0006 
0007 #pragma link C++ class PlotFile;
0008 #endif
0009 
0010 #ifndef __CINT__
0011 #include <stdio.h>
0012 #include <stdlib.h>
0013 #include <fstream>
0014 #include <iostream>
0015 #include <iomanip>
0016 #include <string>
0017 #include <sys/types.h>
0018 #include <sys/stat.h>
0019 #include <dirent.h>
0020 #include "math.h"
0021 #include "string.h"
0022 
0023 #include "TROOT.h"
0024 #include "TFile.h"
0025 #include "TChain.h"
0026 #include "TH1D.h"
0027 #include "TH2D.h"
0028 #include "TH3D.h"
0029 #include "THnSparse.h"
0030 #include "TStyle.h"
0031 #include "TCanvas.h"
0032 #include "TProfile.h"
0033 #include "TTree.h"
0034 #include "TNtuple.h"
0035 #include "TRandom3.h"
0036 #include "TMath.h"
0037 #include "TSystem.h"
0038 #include "TUnixSystem.h"
0039 #include "TVector2.h"
0040 #include "TVector3.h"
0041 #include "TLorentzVector.h"
0042 #include "TTreeReader.h"
0043 #include "TTreeReaderValue.h"
0044 #include "TTreeReaderArray.h"
0045 #include "TLatex.h"
0046 #include "TMinuit.h"
0047 #include "Math/Functor.h"
0048 #include "Fit/Fitter.h"
0049 #include "Math/Minimizer.h"
0050 #endif
0051 
0052 
0053 #include "StPhysicalHelix.h"
0054 #include "SystemOfUnits.h"
0055 #include "PhysicalConstants.h"
0056 
0057 using namespace std;
0058 
0059 StPhysicalHelix* gHelix1 = nullptr;
0060 StPhysicalHelix* gHelix2 = nullptr;
0061 
0062 const double gPionMass = 0.13957;
0063 const double gKaonMass = 0.493677;
0064 
0065 const double twoPi = 2.*3.1415927;
0066 const double eMass = 0.000511;
0067 
0068 const double bField = -1.7; // Tesla
0069 
0070 TVector3 getDcaToVtx(const int index, TVector3 vtx);
0071 void fcnVertexFit(int& npar, double* grad, double& fval, double* par, int iflag);
0072 void getDecayVertex_Chi2fit(const int index1, const int index2, double &s1, double &s2, TVector3 &vertex, double &chi2_ndf, double * parFitErr);
0073 
0074 TLorentzVector getPairParent(const int index1, const int index2, TVector3 vtx,
0075                  float &dcaDaughters, float &cosTheta, float &cosTheta_xy, float &decayLength, float &V0DcaToVtx, float &sigma_vtx, TVector3 &decayVertex,TVector3 &decayVertex_ana, double &chi2_ndf,  double * parFitErr);
0076 
0077 TTreeReaderArray<float> *rcMomPx2;
0078 TTreeReaderArray<float> *rcMomPy2;
0079 TTreeReaderArray<float> *rcMomPz2;
0080 TTreeReaderArray<float> *rcCharge2;
0081 
0082 TTreeReaderArray<float> *rcTrkLoca2;
0083 TTreeReaderArray<float> *rcTrkLocb2;
0084 TTreeReaderArray<float> *rcTrkTheta2;
0085 TTreeReaderArray<float> *rcTrkPhi2;
0086 TTreeReaderArray<std::array<float, 21>> *rcTrkCov;
0087 
0088 
0089 
0090 // Define function for Chi2 minimization between two helices    
0091 struct Chi2Minimization {    
0092      StPhysicalHelix fhelix1, fhelix2;
0093      std::array<float, 21> fcov1, fcov2; // full covariance matrix
0094      
0095      Chi2Minimization(StPhysicalHelix helix1, StPhysicalHelix helix2, std::array<float, 21> cov1, std::array<float, 21> cov2) : fhelix1(helix1),fhelix2(helix2), fcov1(cov1), fcov2(cov2) {}
0096     // Implementation of the function to be minimized
0097     double operator() (const double *par) {
0098     double x = par[0];
0099     double y = par[1];
0100     double z = par[2];
0101     double s1 = par[3];
0102     double s2 = par[4];   
0103     double f = 0;
0104     TVector3 vertex(x, y, z);
0105     TVector3 p1 = fhelix1.at(s1);
0106     TVector3 p2 = fhelix2.at(s2);
0107     TVector3 mom1 = fhelix1.momentumAt(s1,  bField * tesla);
0108     TVector3 mom2 = fhelix2.momentumAt(s2,  bField * tesla);
0109     
0110    // x= −l0 sinϕ , y=l0 cosϕ , z=l
0111    // Recalculate l0 at PCA for error propagation
0112      float l0_track1 = p1.Pt(); float l1_track1 = p1.Z(); double phi_track1 = mom1.Phi();
0113      float l0_track2 = p2.Pt(); float l1_track2 = p2.Z(); double phi_track2 = mom2.Phi();
0114      // Track1: σx^2​=sin^2ϕ⋅σℓ0​^2​+ℓ0^2​cos^2ϕ⋅σϕ^2​+2⋅ℓ0​sinϕcosϕ⋅Cov(ℓ0​,ϕ)
0115      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]; 
0116      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];
0117   
0118      // σy^2​=cos^2ϕ⋅σℓ0^​2​+ℓ0^2​sin^2ϕ⋅σϕ^2​−2⋅ℓ0​sinϕcosϕ⋅Cov(ℓ0​,ϕ)
0119      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]; 
0120      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]; 
0121      
0122      // σz^2​
0123      float sigz1_2 = fcov1[2];
0124      float sigz2_2 = fcov2[2];  
0125      // convert to mm
0126     double d1_x = 10.*(vertex - p1).X(); double d2_x = 10.*(vertex - p2).X(); 
0127     double d1_y = 10.*(vertex - p1).Y(); double d2_y = 10.*(vertex - p2).Y();
0128     double d1_z = 10.*(vertex - p1).Z(); double d2_z = 10.*(vertex - p2).Z();
0129     
0130      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
0131       return f;
0132     }
0133     };
0134 
0135 int main(int argc, char **argv)
0136 {
0137     
0138   if(argc!=3 && argc!=1) return 0;
0139 
0140   TString listname;
0141   TString outname;
0142 
0143   if(argc==1)
0144     {
0145       listname  = "/lustrehome/skumar/eic/my_analysis/D0_Sample_Chi2/test.list";
0146       outname = "test.root";
0147     }
0148 
0149   if(argc==3)
0150     {
0151       listname = argv[1];
0152       outname = argv[2];
0153     }  
0154 
0155   TChain *chain = new TChain("events");
0156 
0157   int nfiles = 0;
0158   char filename[512];
0159   ifstream *inputstream = new ifstream;
0160   inputstream->open(listname.Data());
0161   if (!inputstream->is_open())
0162   {
0163   printf("[e] Cannot open file list: %s\n", listname.Data());
0164   return 0; // or handle as needed
0165   }
0166 
0167   while (inputstream->good())
0168   {
0169   inputstream->getline(filename, 512);
0170   if (inputstream->good())
0171   {
0172     TFile *ftmp = TFile::Open(filename, "READ");
0173     if (!ftmp || !ftmp->IsOpen() || !ftmp->GetNkeys())
0174     {
0175       printf("[e] Skipping bad file: %s\n", filename);
0176       if (ftmp) { ftmp->Close(); delete ftmp; }
0177       continue; 
0178     }
0179     cout << "[i] Add " << nfiles << "th file: " << filename << endl;
0180     chain->Add(filename);
0181     nfiles++;
0182 
0183     ftmp->Close(); // cleanup
0184     delete ftmp;
0185   }
0186 }
0187 
0188   inputstream->close();
0189 
0190   printf("[i] Read in %d files with %lld events in total\n", nfiles, chain->GetEntries());
0191 
0192   TH1F *hEventStat = new TH1F("hEventStat", "Event statistics", 7, 0, 7);
0193   hEventStat->GetXaxis()->SetBinLabel(1, "MC events");
0194   hEventStat->GetXaxis()->SetBinLabel(2, "D0");
0195   hEventStat->GetXaxis()->SetBinLabel(3, "D0 -> pi+K");
0196   hEventStat->GetXaxis()->SetBinLabel(4, "Reco D0");
0197   hEventStat->GetXaxis()->SetBinLabel(5, "Reco Signal D0");
0198   hEventStat->GetXaxis()->SetBinLabel(6, "Reco Signal D0bar");
0199   hEventStat->GetXaxis()->SetBinLabel(7, "Reco Bkg D0");
0200 
0201   TH1F *hMcMult = new TH1F("hMcMult", "MC multiplicity (|#eta| < 3.5);N_{MC}", 50, 0, 50);
0202   
0203   // Secondary vertex with chi2 fit method
0204   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); 
0205   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); 
0206   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);   
0207  
0208   // Secondary vertex with analytical method
0209   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); 
0210   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); 
0211   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);   
0212  
0213   // Secondary vertex  (XY) fit method
0214   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); 
0215 
0216   // Secondary vertex  (XY) analytical method
0217   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); 
0218     
0219   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); 
0220   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); 
0221   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); 
0222   
0223   TH1F *hchi2_vtx = new TH1F("hchi2_vtx", "Helix Calculation: Chi2/ndf; #chi^{2}/ndf; Entries (a.u.)", 1000, 0.0, 50.0); 
0224   TH1F *hchi2_vtx_sig = new TH1F("hchi2_vtx_sig", "Helix Calculation: Chi2/ndf; #chi^{2}/ndf; Entries (a.u.)", 1000, 0.0, 50.0); 
0225   TH1F *hchi2_vtx_bkg = new TH1F("hchi2_vtx_bkg", "Helix Calculation: Chi2/ndf; #chi^{2}/ndf; Entries (a.u.)", 1000, 0.0, 50.0);   
0226 
0227   TH1F *hMcVtxX = new TH1F("hMcVtxX", "x position of MC vertex;x (mm)", 100, -5.05, 4.95);
0228   TH1F *hMcVtxY = new TH1F("hMcVtxY", "y position of MC vertex;y (mm)", 500, -5.01, 4.99);
0229   TH1F *hMcVtxZ = new TH1F("hMcVtxZ", "z position of MC vertex;z (mm)", 400, -200, 200);
0230   
0231   TH1F *hPullVtxX = new TH1F("hPullVtxX", "Pull x position of MC vertex;(Vx_{rec}-Vx_{mc})/#sigma_{vx}", 100, -5.05, 4.95);
0232   TH1F *hPullVtxY = new TH1F("hPullVtxY", "Pull y position of MC vertex;(Vy_{rec}-Vy_{mc})/#sigma_{vy}", 500, -5.01, 4.99);
0233   TH1F *hPullVtxZ = new TH1F("hPullVtxZ", "Pull z position of MC vertex;(Vz_{rec}-Vz_{mc})/#sigma_{vz}", 400, -200, 200);
0234 
0235   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);
0236   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);
0237 
0238   TH2F *hMCD0PtRap = new TH2F("hMCD0PtRap", "MC D^{0};y;p_{T} (GeV/c)", 20, -5, 5, 100, 0, 10);
0239 
0240   TH2F *hMcPiPtEta = new TH2F("hMcPiPtEta", "MC #pi from D^{0} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0241   TH2F *hMcPiPtEtaReco = new TH2F("hMcPiPtEtaReco", "RC #pi from D^{0} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0242 
0243   TH2F *hMcKPtEta = new TH2F("hMcKPtEta", "MC K from D^{0} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0244   TH2F *hMcKPtEtaReco = new TH2F("hMcKPtEtaReco", "RC K from D^{0} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0245 
0246   TH1F *hNRecoVtx = new TH1F("hNRecoVtx", "Number of reconstructed vertices;N", 10, 0, 10);
0247 
0248   const char* part_name[3] = {"Pi", "K", "P"};
0249   const char* part_title[3] = {"#pi", "K", "P"};
0250   TH3F *hRcSecPartLocaToRCVtx[2];
0251   TH3F *hRcSecPartLocbToRCVtx[2];
0252   TH3F *hRcPrimPartLocaToRCVtx[2];
0253   TH3F *hRcPrimPartLocbToRCVtx[2];
0254   for(int i=0; i<2; i++)
0255     {
0256       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);
0257       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);
0258       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);
0259       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);
0260     }
0261 
0262   const char* axis_name[3] = {"x", "y", "z"};
0263   const int nDimDca = 4;
0264   const int nBinsDca[nDimDca] = {50, 20, 500, 50};
0265   const double minBinDca[nDimDca] = {0, -5, -1+0.002, 0};
0266   const double maxBinDca[nDimDca] = {5, 5, 1+0.002, 50};
0267   THnSparseF *hPrimTrkDcaToRCVtx[3][3];
0268   for(int i=0; i<3; i++)
0269     {
0270       for(int j=0; j<3; j++)
0271     {
0272       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);
0273     }
0274     }
0275 
0276   TH3F *h3PairDca12[2];
0277   TH3F *h3PairCosTheta[2];
0278   TH3F *h3PairDca[2];
0279   TH3F *h3PairDecayLength[2];
0280   const char* pair_name[2] = {"signal", "bkg"};
0281   const char* pair_title[2] = {"Signal", "Background"};
0282   for(int i=0; i<2; i++)
0283     {
0284       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);
0285 
0286       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);
0287 
0288       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);
0289 
0290       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);
0291     }
0292 
0293   // Invariant mass
0294   const char* cut_name[2] = {"all", "DCA"};
0295   TH3F *h3InvMass[2][2];
0296   for(int i=0; i<2; i++)
0297     {
0298       for(int j=0; j<2; j++)
0299     {
0300       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);
0301     }
0302     }
0303 
0304   TTreeReader treereader(chain);
0305   // MC  
0306   TTreeReaderArray<int> mcPartGenStatus = {treereader, "MCParticles.generatorStatus"};
0307   TTreeReaderArray<int> mcPartPdg = {treereader, "MCParticles.PDG"};
0308   TTreeReaderArray<float> mcPartCharge = {treereader, "MCParticles.charge"};
0309   TTreeReaderArray<unsigned int> mcPartParent_begin = {treereader, "MCParticles.parents_begin"};
0310   TTreeReaderArray<unsigned int> mcPartParent_end = {treereader, "MCParticles.parents_end"};
0311   TTreeReaderArray<int> mcPartParent_index = {treereader, "_MCParticles_parents.index"};
0312   TTreeReaderArray<unsigned int> mcPartDaughter_begin = {treereader, "MCParticles.daughters_begin"};
0313   TTreeReaderArray<unsigned int> mcPartDaughter_end = {treereader, "MCParticles.daughters_end"};
0314   TTreeReaderArray<int> mcPartDaughter_index = {treereader, "_MCParticles_daughters.index"};
0315   TTreeReaderArray<double> mcPartMass = {treereader, "MCParticles.mass"};
0316   TTreeReaderArray<double> mcPartVx = {treereader, "MCParticles.vertex.x"};
0317   TTreeReaderArray<double> mcPartVy = {treereader, "MCParticles.vertex.y"};
0318   TTreeReaderArray<double> mcPartVz = {treereader, "MCParticles.vertex.z"};
0319   TTreeReaderArray<float> mcMomPx = {treereader, "MCParticles.momentum.x"};
0320   TTreeReaderArray<float> mcMomPy = {treereader, "MCParticles.momentum.y"};
0321   TTreeReaderArray<float> mcMomPz = {treereader, "MCParticles.momentum.z"};
0322   TTreeReaderArray<double> mcEndPointX = {treereader, "MCParticles.endpoint.x"};
0323   TTreeReaderArray<double> mcEndPointY = {treereader, "MCParticles.endpoint.y"};
0324   TTreeReaderArray<double> mcEndPointZ = {treereader, "MCParticles.endpoint.z"};
0325 
0326   TTreeReaderArray<unsigned int> assocChSimID = {treereader, "ReconstructedChargedParticleAssociations.simID"};
0327   TTreeReaderArray<unsigned int> assocChRecID = {treereader, "ReconstructedChargedParticleAssociations.recID"};
0328   TTreeReaderArray<float> assocWeight = {treereader, "ReconstructedChargedParticleAssociations.weight"};
0329  
0330   TTreeReaderArray<float> rcMomPx = {treereader, "ReconstructedChargedParticles.momentum.x"};
0331   TTreeReaderArray<float> rcMomPy = {treereader, "ReconstructedChargedParticles.momentum.y"};
0332   TTreeReaderArray<float> rcMomPz = {treereader, "ReconstructedChargedParticles.momentum.z"};
0333   TTreeReaderArray<float> rcPosx = {treereader, "ReconstructedChargedParticles.referencePoint.x"};
0334   TTreeReaderArray<float> rcPosy = {treereader, "ReconstructedChargedParticles.referencePoint.y"};
0335   TTreeReaderArray<float> rcPosz = {treereader, "ReconstructedChargedParticles.referencePoint.z"};
0336   TTreeReaderArray<float> rcCharge = {treereader, "ReconstructedChargedParticles.charge"};
0337   TTreeReaderArray<int>   rcPdg = {treereader, "ReconstructedChargedParticles.PDG"};
0338 
0339   TTreeReaderArray<float> rcTrkLoca = {treereader, "CentralCKFTrackParameters.loc.a"};
0340   TTreeReaderArray<float> rcTrkLocb = {treereader, "CentralCKFTrackParameters.loc.b"};
0341   TTreeReaderArray<float> rcTrkqOverP = {treereader, "CentralCKFTrackParameters.qOverP"};
0342   TTreeReaderArray<float> rcTrkTheta = {treereader, "CentralCKFTrackParameters.theta"};
0343   TTreeReaderArray<float> rcTrkPhi = {treereader, "CentralCKFTrackParameters.phi"};
0344 
0345   rcMomPx2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.momentum.x"};
0346   rcMomPy2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.momentum.y"};
0347   rcMomPz2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.momentum.z"};
0348   rcCharge2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.charge"};
0349 
0350   rcTrkLoca2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.loc.a"};
0351   rcTrkLocb2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.loc.b"};
0352   rcTrkTheta2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.theta"};
0353   rcTrkPhi2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.phi"};
0354   rcTrkCov = new TTreeReaderArray<std::array<float, 21>>{treereader, "CentralCKFTrackParameters.covariance.covariance[21]"};
0355 
0356   TTreeReaderArray<float> CTVx = {treereader, "CentralTrackVertices.position.x"};
0357   TTreeReaderArray<float> CTVy = {treereader, "CentralTrackVertices.position.y"};
0358   TTreeReaderArray<float> CTVz = {treereader, "CentralTrackVertices.position.z"};
0359   TTreeReaderArray<int> CTVndf = {treereader, "CentralTrackVertices.ndf"};
0360   TTreeReaderArray<float> CTVchi2 = {treereader, "CentralTrackVertices.chi2"};
0361   TTreeReaderArray<float> CTVerr_xx = {treereader, "CentralTrackVertices.positionError.xx"};
0362   TTreeReaderArray<float> CTVerr_yy = {treereader, "CentralTrackVertices.positionError.yy"};
0363   TTreeReaderArray<float> CTVerr_zz = {treereader, "CentralTrackVertices.positionError.zz"};
0364 
0365   TTreeReaderArray<int> prim_vtx_index = {treereader, "PrimaryVertices_objIdx.index"};
0366 
0367   TTreeReaderArray<unsigned int> vtxAssocPart_begin = {treereader, "CentralTrackVertices.associatedParticles_begin"};
0368   TTreeReaderArray<unsigned int> vtxAssocPart_end = {treereader, "CentralTrackVertices.associatedParticles_end"};
0369   TTreeReaderArray<int> vtxAssocPart_index = {treereader, "_CentralTrackVertices_associatedParticles.index"};
0370 
0371   
0372     // Create a ROOT file to store the Ntuple
0373    TFile *file_signal = new TFile("SignalD0.root", "RECREATE");
0374    TTree *tree_sig = new TTree("treeMLSig", "treeMLSig"); 
0375 
0376   // Define variables to store in the Ntuple
0377   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;
0378   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;
0379   
0380     // Link the variables to the TTree branches
0381   tree_sig->Branch("d0_pi", &d0_pi_sig, "d0_pi/F");
0382   tree_sig->Branch("d0_k", &d0_k_sig, "d0_k/F"); 
0383   tree_sig->Branch("d0xy_pi", &d0xy_pi_sig, "d0xy_pi/F");
0384   tree_sig->Branch("d0xy_k", &d0xy_k_sig, "d0xy_k/F");
0385   tree_sig->Branch("sum_d0xy", &sum_d0xy_sig, "sum_d0xy/F");        
0386   tree_sig->Branch("dca_12", &dca_12_sig, "dca_12/F");
0387   tree_sig->Branch("dca_D0", &dca_D0_sig, "dca_D0/F");
0388   tree_sig->Branch("pt_D0", &pt_D0_sig, "pt_D0/F");  
0389   tree_sig->Branch("y_D0", &y_D0_sig, "y_D0/F");  
0390   tree_sig->Branch("mass_D0", &mass_D0_sig, "mass_D0/F");              
0391   tree_sig->Branch("decay_length", &decay_length_sig, "decay_length/F");   
0392   tree_sig->Branch("costheta", &costheta_sig, "costheta/F"); 
0393   tree_sig->Branch("costheta_xy", &costhetaxy_sig, "costheta_xy/F"); 
0394   tree_sig->Branch("sigma_vtx", &sigma_vtx_sig, "sigma_vtx/F"); 
0395   tree_sig->Branch("mult", &mult_sig, "mult/F"); 
0396   tree_sig->Branch("signif_d0xy_pi", &signif_d0xy_pi_sig, "signif_d0xy_pi/F");               
0397   tree_sig->Branch("signif_d0xy_k", &signif_d0xy_k_sig, "signif_d0xy_k/F"); 
0398   
0399   TFile *file_bkg = new TFile("BkgD0.root", "RECREATE");
0400   TTree *tree_bkg = new TTree("treeMLBkg", "treeMLBkg"); 
0401   
0402   // Define variables to store in the Ntuple
0403   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;
0404   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;
0405   // Link the variables to the TTree branches
0406   tree_bkg->Branch("d0_pi", &d0_pi_bkg, "d0_pi/F");
0407   tree_bkg->Branch("d0_k", &d0_k_bkg, "d0_k/F");
0408   tree_bkg->Branch("d0xy_pi", &d0xy_pi_bkg, "d0xy_pi/F");
0409   tree_bkg->Branch("d0xy_k", &d0xy_k_bkg, "d0xy_k/F");
0410   tree_bkg->Branch("sum_d0xy", &sum_d0xy_bkg, "sum_d0xy/F");            
0411   tree_bkg->Branch("dca_12", &dca_12_bkg, "dca_12/F");
0412   tree_bkg->Branch("dca_D0", &dca_D0_bkg, "dca_D0/F");
0413   tree_bkg->Branch("pt_D0", &pt_D0_bkg, "pt_D0/F");  
0414   tree_bkg->Branch("y_D0", &y_D0_bkg, "y_D0/F");  
0415   tree_bkg->Branch("mass_D0", &mass_D0_bkg, "mass_D0/F");              
0416   tree_bkg->Branch("decay_length", &decay_length_bkg, "decay_length/F");   
0417   tree_bkg->Branch("costheta", &costheta_bkg, "costheta/F");  
0418   tree_bkg->Branch("costheta_xy", &costhetaxy_bkg, "costheta_xy/F"); 
0419   tree_bkg->Branch("sigma_vtx", &sigma_vtx_bkg, "sigma_vtx/F"); 
0420   tree_bkg->Branch("mult", &mult_bkg, "mult/F");
0421   tree_bkg->Branch("signif_d0xy_pi", &signif_d0xy_pi_bkg, "signif_d0xy_pi/F");  
0422   tree_bkg->Branch("signif_d0xy_k", &signif_d0xy_k_bkg, "signif_d0xy_k/F");  
0423   
0424   int nevents = 0;
0425   int mult_charged = 0;
0426   while(treereader.Next())
0427     {
0428       if(nevents%1000==0) printf("\n[i] New event %d\n",nevents);
0429       //if(nevents==20) break;
0430 
0431       // find MC primary vertex
0432       int nMCPart = mcPartMass.GetSize();
0433 
0434       TVector3 vertex_mc(-999., -999., -999.);
0435       for(int imc=0; imc<nMCPart; imc++)
0436     {
0437       if(mcPartGenStatus[imc] == 4 && mcPartPdg[imc] == 11)
0438         {
0439           vertex_mc.SetXYZ(mcEndPointX[imc], mcEndPointY[imc], mcEndPointZ[imc]);
0440           break;
0441         }
0442     }
0443       hEventStat->Fill(0.5);
0444       hMcVtxX->Fill(vertex_mc.x());
0445       hMcVtxY->Fill(vertex_mc.y());
0446       hMcVtxZ->Fill(vertex_mc.z());
0447 
0448       // get RC primary vertex
0449       TVector3 vertex_rc(-999., -999., -999.);
0450       TVector3 err_vertex_rc(-999., -999., -999.);
0451       if(prim_vtx_index.GetSize()>0)
0452     {
0453       int rc_vtx_index = prim_vtx_index[0];
0454       vertex_rc.SetXYZ(CTVx[rc_vtx_index], CTVy[rc_vtx_index], CTVz[rc_vtx_index]);
0455       err_vertex_rc.SetXYZ(sqrt(CTVerr_xx[rc_vtx_index]), sqrt(CTVerr_yy[rc_vtx_index]), sqrt(CTVerr_zz[rc_vtx_index]));
0456     }
0457     hPullVtxX->Fill((vertex_rc.x()-vertex_mc.x())/err_vertex_rc.x()); 
0458     hPullVtxY->Fill((vertex_rc.y()-vertex_mc.y())/err_vertex_rc.y()); 
0459     hPullVtxZ->Fill((vertex_rc.z()-vertex_mc.z())/err_vertex_rc.z());   
0460 
0461       // map MC and RC particles
0462       int nAssoc = assocChRecID.GetSize();
0463       map<int, int> assoc_map_to_rc;
0464       map<int, int> assoc_map_to_mc;
0465       
0466       for(unsigned int rc_index=0; rc_index<rcMomPx.GetSize(); rc_index++)
0467     {
0468       // loop over the association to find the matched MC particle
0469       // with largest weight
0470       double max_weight = 0;
0471       int matched_mc_index = -1;
0472       for(int j=0; j<nAssoc; j++)
0473         {
0474           if(assocChRecID[j] != rc_index) continue;
0475           if(assocWeight[j] > max_weight)
0476         {
0477           max_weight = assocWeight[j];
0478           matched_mc_index = assocChSimID[j];
0479         }
0480         }
0481 
0482       // build the map
0483       assoc_map_to_rc[matched_mc_index] = rc_index;
0484       assoc_map_to_mc[rc_index] = matched_mc_index;
0485     }
0486 
0487       // Loop over primary particles
0488       int nMcPart = 0;
0489       for(int imc=0; imc<nMCPart; imc++)
0490     {
0491       if(mcPartGenStatus[imc] == 1 && mcPartCharge[imc] != 0)
0492         {
0493           double dist = sqrt( pow(mcPartVx[imc]-vertex_mc.x(),2) + pow(mcPartVy[imc]-vertex_mc.y(),2) + pow(mcPartVz[imc]-vertex_mc.z(),2));      
0494           if(dist < 1e-4)
0495         {
0496           // count charged particles within |eta| < 3.5
0497           TVector3 mc_mom(mcMomPx[imc], mcMomPy[imc], mcMomPz[imc]);
0498           double mcEta = mc_mom.PseudoRapidity();
0499           if(fabs(mcEta) < 3.5) nMcPart++;
0500         }
0501         }
0502     }
0503       mult_charged = nMcPart;
0504       hMcMult->Fill(nMcPart);
0505       
0506       for(int imc=0; imc<nMCPart; imc++)
0507     {
0508       if(mcPartGenStatus[imc] == 1 && mcPartCharge[imc] != 0)
0509         {
0510           double dist = sqrt( pow(mcPartVx[imc]-vertex_mc.x(),2) + pow(mcPartVy[imc]-vertex_mc.y(),2) + pow(mcPartVz[imc]-vertex_mc.z(),2));      
0511           if(dist < 1e-4)
0512         {         
0513           // check if the MC particle is reconstructed
0514           int rc_index = -1;
0515           if(assoc_map_to_rc.find(imc) != assoc_map_to_rc.end()) rc_index = assoc_map_to_rc[imc];
0516 
0517           if(rc_index>=0)
0518             {
0519               TVector3 dcaToVtx = getDcaToVtx(rc_index, vertex_rc);
0520               
0521               int ip = -1;
0522               if(fabs(mcPartPdg[imc]) == 211) ip = 0;
0523               if(fabs(mcPartPdg[imc]) == 321) ip = 1;
0524               if(fabs(mcPartPdg[imc]) == 2212) ip = 2;
0525               if(ip>=0)
0526             {
0527               TVector3 mom(rcMomPx[rc_index], rcMomPy[rc_index], rcMomPz[rc_index]);
0528               if(ip<2)
0529                 {
0530                   hRcPrimPartLocaToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.Perp());
0531                   hRcPrimPartLocbToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.z());
0532                 }
0533 
0534               double fill1[] = {mom.Pt(), mom.Eta(), dcaToVtx.x(), nMcPart*1.};
0535               double fill2[] = {mom.Pt(), mom.Eta(), dcaToVtx.y(), nMcPart*1.};
0536               double fill3[] = {mom.Pt(), mom.Eta(), dcaToVtx.z(), nMcPart*1.};
0537               hPrimTrkDcaToRCVtx[ip][0]->Fill(fill1);
0538               hPrimTrkDcaToRCVtx[ip][1]->Fill(fill2);
0539               hPrimTrkDcaToRCVtx[ip][2]->Fill(fill3);
0540             }
0541             }
0542         }
0543             }
0544     }
0545 
0546       // look for D0
0547       bool hasD0 = false;
0548       vector<int> mc_index_D0_pi;
0549       vector<int> mc_index_D0_k;
0550       mc_index_D0_pi.clear();
0551       mc_index_D0_k.clear();
0552       
0553       for(int imc=0; imc<nMCPart; imc++)
0554     {
0555       if(fabs(mcPartPdg[imc]) != 421) continue;
0556       hEventStat->Fill(1.5);
0557       
0558       int nDuaghters = mcPartDaughter_end[imc]-mcPartDaughter_begin[imc];
0559       if(nDuaghters!=2) continue;
0560 
0561       // find D0 that decay into pi+K
0562       bool is_pik_decay = false;      
0563       int daug_index_1 = mcPartDaughter_index[mcPartDaughter_begin[imc]];
0564       int daug_index_2 = mcPartDaughter_index[mcPartDaughter_begin[imc]+1];
0565       int daug_pdg_1 = mcPartPdg[daug_index_1];
0566       int daug_pdg_2 = mcPartPdg[daug_index_2];
0567       if( (fabs(daug_pdg_1)==321 && fabs(daug_pdg_2)==211) || (fabs(daug_pdg_1)==211 && fabs(daug_pdg_2)==321) )
0568         {
0569           is_pik_decay = true;
0570         }
0571       if(!is_pik_decay) continue;
0572       if(fabs(daug_pdg_1)==211)
0573         {
0574           mc_index_D0_pi.push_back(daug_index_1);
0575           mc_index_D0_k.push_back(daug_index_2);
0576         }
0577       else
0578         {
0579           mc_index_D0_pi.push_back(daug_index_2);
0580           mc_index_D0_k.push_back(daug_index_1);
0581         }
0582       hasD0 = true;
0583       hEventStat->Fill(2.5);
0584 
0585       // D0 kinematics
0586       TLorentzVector mc_mom_vec;
0587       mc_mom_vec.SetXYZM(mcMomPx[imc], mcMomPy[imc], mcMomPz[imc], mcPartMass[imc]);
0588       double mcRap = mc_mom_vec.Rapidity();
0589       double mcPt = mc_mom_vec.Pt();
0590       hMCD0PtRap->Fill(mcRap, mcPt);
0591 
0592       // decay dauther kinematics
0593       for(int ip = 0; ip<2; ip++)
0594         {
0595           int mc_part_index;
0596           if(ip==0) mc_part_index = mc_index_D0_pi[mc_index_D0_pi.size()-1];
0597           if(ip==1) mc_part_index = mc_index_D0_k[mc_index_D0_k.size()-1];
0598           
0599           TLorentzVector mc_part_vec;
0600           mc_part_vec.SetXYZM(mcMomPx[mc_part_index], mcMomPy[mc_part_index], mcMomPz[mc_part_index], mcPartMass[mc_part_index]);
0601           if(ip==0) hMcPiPtEta->Fill(mc_part_vec.Eta(), mc_part_vec.Pt());
0602           if(ip==1) hMcKPtEta->Fill(mc_part_vec.Eta(), mc_part_vec.Pt());
0603           
0604           int rc_part_index = -1;
0605           if(assoc_map_to_rc.find(mc_part_index) != assoc_map_to_rc.end()) rc_part_index = assoc_map_to_rc[mc_part_index];
0606           if(rc_part_index>=0)
0607         {
0608           TVector3 dcaToVtx = getDcaToVtx(rc_part_index, vertex_rc);
0609           
0610           TVector3 mom(rcMomPx[rc_part_index], rcMomPy[rc_part_index], rcMomPz[rc_part_index]);
0611           hRcSecPartLocaToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.Pt());
0612           hRcSecPartLocbToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.z());
0613           
0614           //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]);
0615         }
0616         }
0617     }
0618 
0619       // Get reconstructed pions and kaons
0620       hNRecoVtx->Fill(CTVx.GetSize());
0621       const int pid_mode = 1; // 0 - truth; 1 - realistic
0622       vector<unsigned int> pi_index;
0623       vector<unsigned int> k_index;
0624       pi_index.clear();
0625       k_index.clear();
0626       for(unsigned int rc_index=0; rc_index<rcMomPx.GetSize(); rc_index++)
0627     {     
0628       if(pid_mode==0)
0629         {
0630           int iSimPartID = -1;
0631           if(assoc_map_to_mc.find(rc_index) != assoc_map_to_mc.end()) iSimPartID = assoc_map_to_mc[rc_index];
0632           if(iSimPartID>=0)
0633         {
0634           if(fabs(mcPartPdg[iSimPartID]) == 211) pi_index.push_back(rc_index);
0635           if(fabs(mcPartPdg[iSimPartID]) == 321) k_index.push_back(rc_index);
0636         }
0637         }
0638       else if(pid_mode==1)
0639         {
0640           if(fabs(rcPdg[rc_index]) == 211) pi_index.push_back(rc_index);
0641           if(fabs(rcPdg[rc_index]) == 321) k_index.push_back(rc_index);
0642         }
0643     }
0644 
0645       // pair pion and kaon
0646       for(unsigned int i=0; i<pi_index.size(); i++)
0647     {
0648       TVector3 dcaToVtx = getDcaToVtx(pi_index[i], vertex_rc);
0649        std::array<float, 21>& cov_pion = rcTrkCov->At(pi_index[i]);
0650      int q_pion = rcCharge[pi_index[i]];
0651       for(unsigned int j=0; j<k_index.size(); j++)
0652         {
0653           TVector3 dcaToVtx2 = getDcaToVtx(k_index[j], vertex_rc);
0654           std::array<float, 21>& cov_kaon = rcTrkCov->At(k_index[j]); 
0655         int q_kaon = rcCharge[k_index[j]];
0656           if(rcCharge[pi_index[i]]*rcCharge[k_index[j]]<0)
0657         {
0658           //printf("[i] Check pair (%d, %d)\n", pi_index[i], k_index[j]);
0659           // -- only look at unlike-sign pi+k pair
0660           bool is_D0_pik = false;
0661           int mc_index_pi = -1, mc_index_k = -1;
0662           if(assoc_map_to_mc.find(pi_index[i]) != assoc_map_to_mc.end()) mc_index_pi = assoc_map_to_mc[pi_index[i]];
0663           if(assoc_map_to_mc.find(k_index[j])  != assoc_map_to_mc.end()) mc_index_k  = assoc_map_to_mc[k_index[j]];
0664 
0665           for(unsigned int k=0; k<mc_index_D0_pi.size(); k++)
0666             {
0667               if(mc_index_pi==mc_index_D0_pi[k] && mc_index_k==mc_index_D0_k[k])
0668             {
0669               is_D0_pik = true;
0670               break;
0671             }
0672             }
0673 
0674           float dcaDaughters, cosTheta, decayLength, V0DcaToVtx, cosTheta_xy, sigma_vtx;
0675           TVector3 decayVertex, decayVertex_ana; 
0676           double chi2_ndf;
0677           double err_Par[5];  // or whatever size is appropriate
0678       //double* err_Par = errParArray;
0679           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);
0680           hchi2_vtx->Fill(chi2_ndf);
0681                   
0682           if(is_D0_pik)
0683             {
0684            TVector3 MCVertex_Kaon(mcPartVx[mc_index_k], mcPartVy[mc_index_k], mcPartVz[mc_index_k]);
0685            TVector3 MCVertex_Pion(mcPartVx[mc_index_pi], mcPartVy[mc_index_pi], mcPartVz[mc_index_pi]);
0686            
0687           // printf("Signal MC Vertex Kaon = (%f, %f, %f)\n",MCVertex_Kaon.X(), MCVertex_Kaon.Y(), MCVertex_Kaon.Z());      
0688           // printf("Signal MC Vertex Pion = (%f, %f, %f)\n",MCVertex_Pion.X(), MCVertex_Pion.Y(), MCVertex_Pion.Z());
0689               hRes_SVx_Helixfit->Fill(parent.Pt(), parent.Rapidity(), (decayVertex.X()-MCVertex_Kaon.X()*0.1)*10);
0690               hRes_SVy_Helixfit->Fill(parent.Pt(), parent.Rapidity(), (decayVertex.Y()-MCVertex_Kaon.Y()*0.1)*10);
0691               hRes_SVz_Helixfit->Fill(parent.Pt(), parent.Rapidity(), (decayVertex.Z()-MCVertex_Kaon.Z()*0.1)*10);  
0692               
0693               hRes_SVx_Helixana->Fill(parent.Pt(), parent.Rapidity(), (decayVertex_ana.X()-MCVertex_Kaon.X()*0.1)*10);
0694               hRes_SVy_Helixana->Fill(parent.Pt(), parent.Rapidity(), (decayVertex_ana.Y()-MCVertex_Kaon.Y()*0.1)*10);
0695               hRes_SVz_Helixana->Fill(parent.Pt(), parent.Rapidity(), (decayVertex_ana.Z()-MCVertex_Kaon.Z()*0.1)*10);
0696 
0697                 // Signed Δv_xy
0698               TVector3 dV = (decayVertex-MCVertex_Kaon*0.1)*10; dV.SetZ(0.);   // remove z-component
0699               TVector3 pTD0V = parent.Vect(); pTD0V.SetZ(0.);
0700               double dvxy_sign = dV.Dot(pTD0V.Unit()); 
0701               int sign = (dvxy_sign >= 0 ? +1 : -1);    
0702               hRes_SVxy_Helixfit->Fill(parent.Pt(), parent.Rapidity(), sign*dV.Mag());
0703               hRes_SVxy_Helixana->Fill(parent.Pt(), parent.Rapidity(), sign*dV.Mag());
0704               
0705               hRes_SVx_Helixfit_pull->Fill(parent.Pt(), parent.Rapidity(), ((decayVertex.X()-MCVertex_Kaon.X()*0.1))/(err_Par[0]));
0706               hRes_SVy_Helixfit_pull->Fill(parent.Pt(), parent.Rapidity(), ((decayVertex.Y()-MCVertex_Kaon.Y()*0.1))/(err_Par[1]));
0707               hRes_SVz_Helixfit_pull->Fill(parent.Pt(), parent.Rapidity(), ((decayVertex.Z()-MCVertex_Kaon.Z()*0.1))/(err_Par[2])); 
0708               
0709               hchi2_vtx_sig->Fill(chi2_ndf);              
0710                                
0711               hEventStat->Fill(3.5);
0712               if (q_kaon == -1 && q_pion == 1)
0713               hEventStat->Fill(4.5);   // D0
0714               else if (q_kaon == 1 && q_pion == -1)
0715               hEventStat->Fill(5.5);   // D0bar
0716               h3PairDca12[0]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters);
0717               h3PairCosTheta[0]->Fill(parent.Pt(), parent.Rapidity(), cosTheta);
0718               h3PairDca[0]->Fill(parent.Pt(), parent.Rapidity(), V0DcaToVtx);
0719               h3PairDecayLength[0]->Fill(parent.Pt(), parent.Rapidity(), decayLength);
0720               //printf("Signal: dca12 = %2.4f, cosTheta = %2.4f, D0dca = %2.4f, decay = %2.4f\n", dcaDaughters, cosTheta, V0DcaToVtx, decayLength);
0721               h3InvMass[0][0]->Fill(parent.Pt(), parent.Rapidity(), parent.M());
0722                
0723            // Toplogical Variables for Signal
0724               d0_pi_sig = dcaToVtx.Mag();
0725               d0_k_sig = dcaToVtx2.Mag();
0726               d0xy_pi_sig = dcaToVtx.Perp();
0727               signif_d0xy_pi_sig = d0xy_pi_sig/sqrt(cov_pion[0]);
0728               d0xy_k_sig = dcaToVtx2.Perp();
0729               signif_d0xy_k_sig = d0xy_k_sig/sqrt(cov_kaon[0]);
0730               sum_d0xy_sig = sqrt(d0xy_pi_sig*d0xy_pi_sig+d0xy_k_sig*d0xy_k_sig);             
0731               dca_12_sig = dcaDaughters;
0732                 dca_D0_sig = V0DcaToVtx;
0733                 decay_length_sig = decayLength;
0734                 costheta_sig = cosTheta;
0735                 costhetaxy_sig = cosTheta_xy;
0736                 pt_D0_sig = parent.Pt();
0737                 y_D0_sig = parent.Rapidity();
0738                 mass_D0_sig = parent.M();
0739                 sigma_vtx_sig = sigma_vtx;
0740                 mult_sig = mult_charged;
0741                 tree_sig->Fill();
0742             }
0743           else
0744             {
0745           hEventStat->Fill(6.5);   // Bkg
0746           hchi2_vtx_bkg->Fill(chi2_ndf);                    
0747               h3PairDca12[1]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters);
0748               h3PairCosTheta[1]->Fill(parent.Pt(), parent.Rapidity(), cosTheta);
0749               h3PairDca[1]->Fill(parent.Pt(), parent.Rapidity(), V0DcaToVtx);
0750               h3PairDecayLength[1]->Fill(parent.Pt(), parent.Rapidity(), decayLength);
0751               
0752               //printf("Bkg: dca12 = %2.4f, cosTheta = %2.4f, D0dca = %2.4f, decay = %2.4f\n", dcaDaughters, cosTheta, V0DcaToVtx, decayLength);
0753               h3InvMass[1][0]->Fill(parent.Pt(), parent.Rapidity(), parent.M());
0754 
0755               // Toplogical Variables for Bkg
0756               d0_pi_bkg = dcaToVtx.Mag();
0757               d0_k_bkg = dcaToVtx2.Mag();
0758               d0xy_pi_bkg = dcaToVtx.Perp();
0759               signif_d0xy_pi_bkg = d0xy_pi_bkg/sqrt(cov_pion[0]);
0760               d0xy_k_bkg = dcaToVtx2.Perp();
0761               signif_d0xy_k_bkg = d0xy_k_bkg/sqrt(cov_kaon[0]);
0762               sum_d0xy_bkg = sqrt(d0xy_pi_bkg*d0xy_pi_bkg+d0xy_k_bkg*d0xy_k_bkg);             
0763               dca_12_bkg = dcaDaughters;
0764                 dca_D0_bkg = V0DcaToVtx;
0765                 decay_length_bkg = decayLength;
0766                 costheta_bkg = cosTheta;
0767                 costhetaxy_bkg = cosTheta_xy;
0768                 pt_D0_bkg = parent.Pt();
0769                 y_D0_bkg = parent.Rapidity();
0770                 mass_D0_bkg = parent.M();
0771                 sigma_vtx_bkg = sigma_vtx;
0772                 mult_bkg = mult_charged;
0773                 tree_bkg->Fill();
0774             }
0775 
0776           if(dcaToVtx.Perp() >= 0.02 && dcaToVtx2.Perp() >= 0.02 &&
0777              dcaDaughters < 0.07 && cosTheta > 0.95 && decayLength > 0.05 && V0DcaToVtx < 0.1)
0778             {
0779               if(is_D0_pik)
0780             {
0781               h3InvMass[0][1]->Fill(parent.Pt(), parent.Rapidity(), parent.M());
0782             }
0783               else
0784             {
0785               h3InvMass[1][1]->Fill(parent.Pt(), parent.Rapidity(), parent.M());
0786             }
0787             }
0788         }
0789         }
0790     }
0791           
0792       nevents++;
0793     }
0794 
0795   file_signal->cd();  
0796   tree_sig->Write();
0797   file_signal->Close();  
0798   
0799   file_bkg->cd();  
0800   tree_bkg->Write();
0801   file_bkg->Close();
0802 
0803   TFile *outfile = new TFile(outname.Data(), "recreate");
0804 
0805   hEventStat->Write();
0806   hMcMult->Write();
0807   hMcVtxX->Write();
0808   hMcVtxY->Write();
0809   hMcVtxZ->Write();
0810   
0811   hPullVtxX->Write();
0812   hPullVtxY->Write();  
0813   hPullVtxZ->Write();
0814   
0815   hRes_SVx_Helixfit->Write();
0816   hRes_SVy_Helixfit->Write();
0817   hRes_SVz_Helixfit->Write(); 
0818   hRes_SVxy_Helixfit->Write();  
0819  
0820   hRes_SVx_Helixana->Write();
0821   hRes_SVy_Helixana->Write();
0822   hRes_SVz_Helixana->Write();
0823   hRes_SVxy_Helixana->Write();    
0824       
0825   hchi2_vtx->Write();
0826   hchi2_vtx_sig->Write();
0827   hchi2_vtx_bkg->Write();  
0828   hRes_SVx_Helixfit_pull->Write();
0829   hRes_SVy_Helixfit_pull->Write(); 
0830   hRes_SVz_Helixfit_pull->Write(); 
0831     
0832   hD0DecayVxVy->Write();
0833   hD0DecayVrVz->Write();
0834   
0835   hMCD0PtRap->Write();
0836 
0837   hMcPiPtEta->Write();
0838   hMcPiPtEtaReco->Write();
0839   hMcKPtEta->Write();
0840   hMcKPtEtaReco->Write();
0841   
0842   hNRecoVtx->Write();
0843 
0844   for(int ip=0; ip<2; ip++)
0845     {
0846       hRcSecPartLocaToRCVtx[ip]->Write();
0847       hRcSecPartLocbToRCVtx[ip]->Write();
0848       hRcPrimPartLocaToRCVtx[ip]->Write();
0849       hRcPrimPartLocbToRCVtx[ip]->Write();
0850     }
0851 
0852   for(int i=0; i<3; i++)
0853     {
0854       for(int j=0; j<3; j++)
0855     {
0856       hPrimTrkDcaToRCVtx[i][j]->Write();
0857     }
0858     }
0859 
0860   for(int i=0; i<2; i++)
0861     {
0862       h3PairDca12[i]->Write();
0863       h3PairCosTheta[i]->Write();
0864       h3PairDca[i]->Write();
0865       h3PairDecayLength[i]->Write();
0866     }
0867 
0868   for(int i=0; i<2; i++)
0869     {
0870       for(int j=0; j<2; j++)
0871     {
0872       h3InvMass[i][j]->Write();
0873     }
0874     }
0875   
0876   
0877   outfile->Close();
0878 }
0879 
0880 //======================================
0881 TVector3 getDcaToVtx(const int index, TVector3 vtx)
0882 {
0883   TVector3 pos(rcTrkLoca2->At(index) * sin(rcTrkPhi2->At(index)) * -1 * millimeter, rcTrkLoca2->At(index) * cos(rcTrkPhi2->At(index)) * millimeter, rcTrkLocb2->At(index) * millimeter);
0884   TVector3 mom(rcMomPx2->At(index), rcMomPy2->At(index), rcMomPz2->At(index));
0885    
0886   StPhysicalHelix pHelix(mom, pos, bField * tesla, rcCharge2->At(index));
0887 
0888   TVector3 vtx_tmp;
0889   vtx_tmp.SetXYZ(vtx.x()*millimeter, vtx.y()*millimeter, vtx.z()*millimeter);
0890   
0891   pHelix.moveOrigin(pHelix.pathLength(vtx_tmp));
0892   TVector3 dcaToVtx = pHelix.origin() - vtx_tmp;
0893 
0894   dcaToVtx.SetXYZ(dcaToVtx.x()/millimeter, dcaToVtx.y()/millimeter, dcaToVtx.z()/millimeter);
0895   
0896   return dcaToVtx;
0897 }
0898 
0899 //======================================
0900 TLorentzVector getPairParent(const int index1, const int index2, TVector3 vtx,
0901                  float &dcaDaughters, float &cosTheta, float &cosTheta_xy, float &decayLength, float &V0DcaToVtx, float &sigma_vtx, TVector3 &decayVertex, TVector3 &decayVertex_ana, double &chi2_ndf, double * parFitErr)
0902 {
0903   // -- get helix
0904   TVector3 pos1(rcTrkLoca2->At(index1) * sin(rcTrkPhi2->At(index1)) * -1 * millimeter, rcTrkLoca2->At(index1) * cos(rcTrkPhi2->At(index1)) * millimeter, rcTrkLocb2->At(index1) * millimeter);
0905   TVector3 mom1(rcMomPx2->At(index1), rcMomPy2->At(index1), rcMomPz2->At(index1));
0906 
0907   TVector3 pos2(rcTrkLoca2->At(index2) * sin(rcTrkPhi2->At(index2)) * -1 * millimeter, rcTrkLoca2->At(index2) * cos(rcTrkPhi2->At(index2)) * millimeter, rcTrkLocb2->At(index2) * millimeter);
0908   TVector3 mom2(rcMomPx2->At(index2), rcMomPy2->At(index2), rcMomPz2->At(index2));
0909 
0910   float charge1 = rcCharge2->At(index1);
0911   float charge2 = rcCharge2->At(index2);
0912   
0913   StPhysicalHelix p1Helix(mom1, pos1, bField * tesla, charge1);
0914   StPhysicalHelix p2Helix(mom2, pos2, bField * tesla, charge2);
0915 
0916   TVector3 vtx_tmp;
0917   vtx_tmp.SetXYZ(vtx.x()*millimeter, vtx.y()*millimeter, vtx.z()*millimeter);
0918   
0919   double s1, s2;
0920   getDecayVertex_Chi2fit(index1, index2, s1, s2, decayVertex, chi2_ndf, parFitErr);
0921  
0922   TVector3 const p1AtDcaToP2 = p1Helix.at(s1);
0923   TVector3 const p2AtDcaToP1 = p2Helix.at(s2);
0924   // printf("p1AtDcaToP2 origin = (%2.4f, %2.4f, %2.4f)\n", p1AtDcaToP2.x(), p1AtDcaToP2.y(), p1AtDcaToP2.z());
0925   // printf("p2AtDcaToP1 origin = (%2.4f, %2.4f, %2.4f)\n", p2AtDcaToP1.x(), p2AtDcaToP1.y(), p2AtDcaToP1.z());
0926   
0927   // -- calculate DCA of particle1 to particle2 at their DCA
0928   dcaDaughters = (p1AtDcaToP2 - p2AtDcaToP1).Mag()/millimeter;
0929     
0930   // -- calculate Lorentz vector of particle1-particle2 pair
0931   TVector3 const p1MomAtDca = p1Helix.momentumAt(s1,  bField * tesla);
0932   TVector3 const p2MomAtDca = p2Helix.momentumAt(s2, bField * tesla);
0933   
0934   TLorentzVector p1FourMom(p1MomAtDca, sqrt(p1MomAtDca.Mag2()+gPionMass*gPionMass));
0935   TLorentzVector p2FourMom(p2MomAtDca, sqrt(p2MomAtDca.Mag2()+gKaonMass*gKaonMass));
0936   
0937   TLorentzVector parent = p1FourMom + p2FourMom;
0938 
0939   // -- calculate decay vertex (secondary or tertiary)
0940   decayVertex_ana = (p1AtDcaToP2 + p2AtDcaToP1) * 0.5 ;
0941   sigma_vtx = sqrt((p1AtDcaToP2-decayVertex).Mag2()+(p2AtDcaToP1-decayVertex).Mag2())/millimeter;
0942     
0943   // -- calculate pointing angle and decay length with respect to primary vertex
0944   //    if decay vertex is a tertiary vertex
0945   //    -> only rough estimate -> needs to be updated after secondary vertex is found
0946   TVector3 vtxToV0 = decayVertex - vtx_tmp;
0947   TVector3 vtxToV0_xy(vtxToV0.x(), vtxToV0.y(), 0.);
0948   TVector3 parent_xy(parent.Vect().x(),parent.Vect().y(),0.);
0949   float pointingAngle = vtxToV0.Angle(parent.Vect());
0950   float pointingAngle_xy = vtxToV0_xy.Angle(parent_xy);
0951   cosTheta = std::cos(pointingAngle);
0952   cosTheta_xy = std::cos(pointingAngle_xy);  
0953   decayLength = vtxToV0.Mag()/millimeter;
0954 
0955   // -- calculate V0 DCA to primary vertex
0956   V0DcaToVtx = decayLength * std::sin(pointingAngle);
0957     
0958   //TVector3 dcaToVtx = getDcaToVtx(parent.Vect(), decayVertex, 0, vtx);
0959   //V0DcaToVtx = dcaToVtx.Mag();
0960   return parent;
0961 }
0962 
0963 
0964 
0965 void getDecayVertex_Chi2fit(const int index1, const int index2, double &s1, double &s2, TVector3 &vertex, double &chi2_ndf, double *parFitErr)
0966 {
0967     TVector3 pos1(rcTrkLoca2->At(index1) * sin(rcTrkPhi2->At(index1)) * -1 * millimeter,
0968                   rcTrkLoca2->At(index1) * cos(rcTrkPhi2->At(index1)) * millimeter,
0969                   rcTrkLocb2->At(index1) * millimeter);
0970 
0971     TVector3 mom1(rcMomPx2->At(index1), rcMomPy2->At(index1), rcMomPz2->At(index1));
0972 
0973     TVector3 pos2(rcTrkLoca2->At(index2) * sin(rcTrkPhi2->At(index2)) * -1 * millimeter,
0974                   rcTrkLoca2->At(index2) * cos(rcTrkPhi2->At(index2)) * millimeter,
0975                   rcTrkLocb2->At(index2) * millimeter);
0976 
0977     TVector3 mom2(rcMomPx2->At(index2), rcMomPy2->At(index2), rcMomPz2->At(index2));
0978 
0979 
0980     float charge1 = rcCharge2->At(index1);
0981     float charge2 = rcCharge2->At(index2);
0982 
0983 
0984     StPhysicalHelix helix1(mom1, pos1, bField * tesla, charge1);
0985     StPhysicalHelix helix2(mom2, pos2, bField * tesla, charge2);
0986    
0987     std::array<float, 21>& cov_track1 = rcTrkCov->At(index1);
0988     std::array<float, 21>& cov_track2 = rcTrkCov->At(index2);  
0989     
0990     pair<double, double> const ss = helix1.pathLengths(helix2);
0991     TVector3 const p1_init = helix1.at(ss.first);
0992     TVector3 const p2_init = helix2.at(ss.second); 
0993     TVector3 const mid_point = 0.5*(p1_init+p2_init); 
0994 
0995       // Perform Minimization
0996     const Int_t nPar = 5;
0997     Chi2Minimization d2Function(helix1,helix2,cov_track1,cov_track2);
0998     ROOT::Math::Functor fcn(d2Function,nPar); // 5 parameters
0999     ROOT::Fit::Fitter fitter;
1000 
1001     double pStart[nPar] = {mid_point.X(),mid_point.Y(),mid_point.Z(),ss.first,ss.second};
1002     fitter.SetFCN(fcn, pStart,nPar,1);
1003     
1004     fitter.Config().ParSettings(0).SetName("x0");
1005     fitter.Config().ParSettings(0).SetStepSize(0.01);
1006    // fitter.Config().ParSettings(0).SetLimits(-1., 1.);    
1007     // No limits for x, y, z
1008 
1009     fitter.Config().ParSettings(1).SetName("y0");
1010     fitter.Config().ParSettings(1).SetStepSize(0.01);
1011     //fitter.Config().ParSettings(1).SetLimits(-1., 1.);
1012 
1013     fitter.Config().ParSettings(2).SetName("z0");
1014     fitter.Config().ParSettings(2).SetStepSize(0.01);
1015    // fitter.Config().ParSettings(2).SetLimits(-10., 10.);    
1016 
1017     fitter.Config().ParSettings(3).SetName("s1");
1018     fitter.Config().ParSettings(3).SetValue(0.0);
1019     fitter.Config().ParSettings(3).SetStepSize(0.01);
1020     //fitter.Config().ParSettings(3).SetLimits(-1., 1.);
1021 
1022     fitter.Config().ParSettings(4).SetName("s2");
1023     fitter.Config().ParSettings(4).SetValue(0.0);
1024     fitter.Config().ParSettings(4).SetStepSize(0.01);
1025     //fitter.Config().ParSettings(4).SetLimits(-1., 1.);      
1026     
1027     fitter.Config().MinimizerOptions().SetMaxIterations(10000); 
1028     // do the fit 
1029 
1030     Bool_t ok = fitter.FitFCN();
1031     if (!ok) Error("Fitting","Fitting failed");
1032     const ROOT::Fit::FitResult & result = fitter.Result();
1033    // double chi2 = fitter.Result().Chi2();
1034     double ndf = 2*3-nPar;
1035     chi2_ndf = fitter.Result().MinFcnValue()/ndf;  // Minimum value of your function
1036     
1037     int status = fitter.Result().Status();
1038   //  if (status>0 ) {printf("Fit Failed!!!!\n");}
1039    // if (status>0 || chi2_ndf>10. ) return;
1040    // cout <<"\033[1;31m Fit Result Chi2:\033[0m"<<chi2_ndf<<endl;
1041    // result.Print(std::cout);
1042    
1043    // Get the covariance matrix
1044  //  TMatrixDSym covMatrix(5);
1045   // result.GetCovarianceMatrix(covMatrix); // Matrix for the parameter errors
1046   // covMatrix.Print();
1047    
1048    const double * parFit = result.GetParams();
1049    const double *FitErr = result.GetErrors();
1050     
1051    for (int i = 0; i < nPar; ++i) parFitErr[i] = FitErr[i];
1052    vertex.SetXYZ(parFit[0], parFit[1], parFit[2]);
1053    s1 = parFit[3]; s2 = parFit[4];
1054 }
1055 
1056
[ 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!