Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-09-17 09:23:23 

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, 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  = "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   TH1F *hRes_SVx_Helixfit = new TH1F("hRes_SVx_Helixfit", "Fit method: Residual of SVx; SVx_{rec}-SVx_{mc} (mm); Entries (a.u.)", 200, -1.0, 1.0);
0203   TH1F *hRes_SVy_Helixfit= new TH1F("hRes_SVy_Helixfit", "Fit method: Residual of SVy; SVy_{rec}-SVy_{mc} (mm); Entries (a.u.)", 200, -1.0, 1.0);
0204   TH1F *hRes_SVz_Helixfit = new TH1F("hRes_SVz_Helixfit", "Fit method: Residual of SVz; SVz_{rec}-SVz_{mc} (mm); Entries (a.u.)", 1000, -5.0, 5.0);
0205   
0206   TH1F *hRes_SVx_Helixfit_pull = new TH1F("hRes_SVx_Helixfit_pull", "Fit method: Pull of SVx; SVx_{rec}-SVx_{mc}/#sigma; Entries (a.u.)", 200, -5.0, 5.0);
0207   TH1F *hRes_SVy_Helixfit_pull = new TH1F("hRes_SVy_Helixfit_pull", "Fit method: Pull of SVy; SVy_{rec}-SVy_{mc}/#sigma; Entries (a.u.)", 200, -5.0, 5.0);
0208   TH1F *hRes_SVz_Helixfit_pull = new TH1F("hRes_SVz_Helixfit_pull", "Fit method: Pull of SVz; SVz_{rec}-SVz_{mc}/#sigma; Entries (a.u.)", 200, -5.0, 5.0);
0209   
0210   TH1F *hchi2_vtx = new TH1F("hchi2_vtx", "Helix Calculation: Chi2/ndf; #chi^{2}/ndf; Entries (a.u.)", 1000, 0.0, 50.0); 
0211   TH1F *hchi2_vtx_sig = new TH1F("hchi2_vtx_sig", "Helix Calculation: Chi2/ndf; #chi^{2}/ndf; Entries (a.u.)", 1000, 0.0, 50.0); 
0212   TH1F *hchi2_vtx_bkg = new TH1F("hchi2_vtx_bkg", "Helix Calculation: Chi2/ndf; #chi^{2}/ndf; Entries (a.u.)", 1000, 0.0, 50.0);   
0213 
0214   TH1F *hMcVtxX = new TH1F("hMcVtxX", "x position of MC vertex;x (mm)", 100, -5.05, 4.95);
0215   TH1F *hMcVtxY = new TH1F("hMcVtxY", "y position of MC vertex;y (mm)", 500, -5.01, 4.99);
0216   TH1F *hMcVtxZ = new TH1F("hMcVtxZ", "z position of MC vertex;z (mm)", 400, -200, 200);
0217   
0218   TH1F *hPullVtxX = new TH1F("hPullVtxX", "Pull x position of MC vertex;(Vx_{rec}-Vx_{mc})/#sigma_{vx}", 100, -5.05, 4.95);
0219   TH1F *hPullVtxY = new TH1F("hPullVtxY", "Pull y position of MC vertex;(Vy_{rec}-Vy_{mc})/#sigma_{vy}", 500, -5.01, 4.99);
0220   TH1F *hPullVtxZ = new TH1F("hPullVtxZ", "Pull z position of MC vertex;(Vz_{rec}-Vz_{mc})/#sigma_{vz}", 400, -200, 200);
0221 
0222   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);
0223   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);
0224 
0225   TH2F *hMCD0PtRap = new TH2F("hMCD0PtRap", "MC D^{0};y;p_{T} (GeV/c)", 20, -5, 5, 100, 0, 10);
0226 
0227   TH2F *hMcPiPtEta = new TH2F("hMcPiPtEta", "MC #pi from D^{0} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0228   TH2F *hMcPiPtEtaReco = new TH2F("hMcPiPtEtaReco", "RC #pi from D^{0} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0229 
0230   TH2F *hMcKPtEta = new TH2F("hMcKPtEta", "MC K from D^{0} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0231   TH2F *hMcKPtEtaReco = new TH2F("hMcKPtEtaReco", "RC K from D^{0} decay;#eta^{MC};p_{T}^{MC} (GeV/c)", 20, -5, 5, 100, 0, 10);
0232 
0233   TH1F *hNRecoVtx = new TH1F("hNRecoVtx", "Number of reconstructed vertices;N", 10, 0, 10);
0234 
0235   const char* part_name[3] = {"Pi", "K", "P"};
0236   const char* part_title[3] = {"#pi", "K", "P"};
0237   TH3F *hRcSecPartLocaToRCVtx[2];
0238   TH3F *hRcSecPartLocbToRCVtx[2];
0239   TH3F *hRcPrimPartLocaToRCVtx[2];
0240   TH3F *hRcPrimPartLocbToRCVtx[2];
0241   for(int i=0; i<2; i++)
0242     {
0243       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);
0244       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);
0245       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);
0246       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);
0247     }
0248 
0249   const char* axis_name[3] = {"x", "y", "z"};
0250   const int nDimDca = 4;
0251   const int nBinsDca[nDimDca] = {50, 20, 500, 50};
0252   const double minBinDca[nDimDca] = {0, -5, -1+0.002, 0};
0253   const double maxBinDca[nDimDca] = {5, 5, 1+0.002, 50};
0254   THnSparseF *hPrimTrkDcaToRCVtx[3][3];
0255   for(int i=0; i<3; i++)
0256     {
0257       for(int j=0; j<3; j++)
0258     {
0259       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);
0260     }
0261     }
0262 
0263   TH3F *h3PairDca12[2];
0264   TH3F *h3PairCosTheta[2];
0265   TH3F *h3PairDca[2];
0266   TH3F *h3PairDecayLength[2];
0267   const char* pair_name[2] = {"signal", "bkg"};
0268   const char* pair_title[2] = {"Signal", "Background"};
0269   for(int i=0; i<2; i++)
0270     {
0271       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);
0272 
0273       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);
0274 
0275       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);
0276 
0277       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);
0278     }
0279 
0280   // Invariant mass
0281   const char* cut_name[2] = {"all", "DCA"};
0282   TH3F *h3InvMass[2][2];
0283   for(int i=0; i<2; i++)
0284     {
0285       for(int j=0; j<2; j++)
0286     {
0287       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);
0288     }
0289     }
0290 
0291   TTreeReader treereader(chain);
0292   // MC  
0293   TTreeReaderArray<int> mcPartGenStatus = {treereader, "MCParticles.generatorStatus"};
0294   TTreeReaderArray<int> mcPartPdg = {treereader, "MCParticles.PDG"};
0295   TTreeReaderArray<float> mcPartCharge = {treereader, "MCParticles.charge"};
0296   TTreeReaderArray<unsigned int> mcPartParent_begin = {treereader, "MCParticles.parents_begin"};
0297   TTreeReaderArray<unsigned int> mcPartParent_end = {treereader, "MCParticles.parents_end"};
0298   TTreeReaderArray<int> mcPartParent_index = {treereader, "_MCParticles_parents.index"};
0299   TTreeReaderArray<unsigned int> mcPartDaughter_begin = {treereader, "MCParticles.daughters_begin"};
0300   TTreeReaderArray<unsigned int> mcPartDaughter_end = {treereader, "MCParticles.daughters_end"};
0301   TTreeReaderArray<int> mcPartDaughter_index = {treereader, "_MCParticles_daughters.index"};
0302   TTreeReaderArray<double> mcPartMass = {treereader, "MCParticles.mass"};
0303   TTreeReaderArray<double> mcPartVx = {treereader, "MCParticles.vertex.x"};
0304   TTreeReaderArray<double> mcPartVy = {treereader, "MCParticles.vertex.y"};
0305   TTreeReaderArray<double> mcPartVz = {treereader, "MCParticles.vertex.z"};
0306   TTreeReaderArray<double> mcMomPx = {treereader, "MCParticles.momentum.x"};
0307   TTreeReaderArray<double> mcMomPy = {treereader, "MCParticles.momentum.y"};
0308   TTreeReaderArray<double> mcMomPz = {treereader, "MCParticles.momentum.z"};
0309   TTreeReaderArray<double> mcEndPointX = {treereader, "MCParticles.endpoint.x"};
0310   TTreeReaderArray<double> mcEndPointY = {treereader, "MCParticles.endpoint.y"};
0311   TTreeReaderArray<double> mcEndPointZ = {treereader, "MCParticles.endpoint.z"};
0312 
0313   TTreeReaderArray<unsigned int> assocChSimID = {treereader, "ReconstructedChargedParticleAssociations.simID"};
0314   TTreeReaderArray<unsigned int> assocChRecID = {treereader, "ReconstructedChargedParticleAssociations.recID"};
0315   TTreeReaderArray<float> assocWeight = {treereader, "ReconstructedChargedParticleAssociations.weight"};
0316  
0317   TTreeReaderArray<float> rcMomPx = {treereader, "ReconstructedChargedParticles.momentum.x"};
0318   TTreeReaderArray<float> rcMomPy = {treereader, "ReconstructedChargedParticles.momentum.y"};
0319   TTreeReaderArray<float> rcMomPz = {treereader, "ReconstructedChargedParticles.momentum.z"};
0320   TTreeReaderArray<float> rcPosx = {treereader, "ReconstructedChargedParticles.referencePoint.x"};
0321   TTreeReaderArray<float> rcPosy = {treereader, "ReconstructedChargedParticles.referencePoint.y"};
0322   TTreeReaderArray<float> rcPosz = {treereader, "ReconstructedChargedParticles.referencePoint.z"};
0323   TTreeReaderArray<float> rcCharge = {treereader, "ReconstructedChargedParticles.charge"};
0324   TTreeReaderArray<int>   rcPdg = {treereader, "ReconstructedChargedParticles.PDG"};
0325 
0326   TTreeReaderArray<float> rcTrkLoca = {treereader, "CentralCKFTrackParameters.loc.a"};
0327   TTreeReaderArray<float> rcTrkLocb = {treereader, "CentralCKFTrackParameters.loc.b"};
0328   TTreeReaderArray<float> rcTrkqOverP = {treereader, "CentralCKFTrackParameters.qOverP"};
0329   TTreeReaderArray<float> rcTrkTheta = {treereader, "CentralCKFTrackParameters.theta"};
0330   TTreeReaderArray<float> rcTrkPhi = {treereader, "CentralCKFTrackParameters.phi"};
0331 
0332   rcMomPx2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.momentum.x"};
0333   rcMomPy2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.momentum.y"};
0334   rcMomPz2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.momentum.z"};
0335   rcCharge2 = new TTreeReaderArray<float>{treereader, "ReconstructedChargedParticles.charge"};
0336 
0337   rcTrkLoca2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.loc.a"};
0338   rcTrkLocb2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.loc.b"};
0339   rcTrkTheta2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.theta"};
0340   rcTrkPhi2 = new TTreeReaderArray<float>{treereader, "CentralCKFTrackParameters.phi"};
0341   rcTrkCov = new TTreeReaderArray<std::array<float, 21>>{treereader, "CentralCKFTrackParameters.covariance.covariance[21]"};
0342 
0343   TTreeReaderArray<float> CTVx = {treereader, "CentralTrackVertices.position.x"};
0344   TTreeReaderArray<float> CTVy = {treereader, "CentralTrackVertices.position.y"};
0345   TTreeReaderArray<float> CTVz = {treereader, "CentralTrackVertices.position.z"};
0346   TTreeReaderArray<int> CTVndf = {treereader, "CentralTrackVertices.ndf"};
0347   TTreeReaderArray<float> CTVchi2 = {treereader, "CentralTrackVertices.chi2"};
0348   TTreeReaderArray<float> CTVerr_xx = {treereader, "CentralTrackVertices.positionError.xx"};
0349   TTreeReaderArray<float> CTVerr_yy = {treereader, "CentralTrackVertices.positionError.yy"};
0350   TTreeReaderArray<float> CTVerr_zz = {treereader, "CentralTrackVertices.positionError.zz"};
0351 
0352   TTreeReaderArray<int> prim_vtx_index = {treereader, "PrimaryVertices_objIdx.index"};
0353 
0354   TTreeReaderArray<unsigned int> vtxAssocPart_begin = {treereader, "CentralTrackVertices.associatedParticles_begin"};
0355   TTreeReaderArray<unsigned int> vtxAssocPart_end = {treereader, "CentralTrackVertices.associatedParticles_end"};
0356   TTreeReaderArray<int> vtxAssocPart_index = {treereader, "_CentralTrackVertices_associatedParticles.index"};
0357 
0358   
0359     // Create a ROOT file to store the Ntuple
0360    TFile *file_signal = new TFile("SignalD0.root", "RECREATE");
0361    TTree *tree_sig = new TTree("treeMLSig", "treeMLSig"); 
0362 
0363   // Define variables to store in the Ntuple
0364   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;
0365   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;
0366   
0367     // Link the variables to the TTree branches
0368   tree_sig->Branch("d0_pi", &d0_pi_sig, "d0_pi/F");
0369   tree_sig->Branch("d0_k", &d0_k_sig, "d0_k/F"); 
0370   tree_sig->Branch("d0xy_pi", &d0xy_pi_sig, "d0xy_pi/F");
0371   tree_sig->Branch("d0xy_k", &d0xy_k_sig, "d0xy_k/F");
0372   tree_sig->Branch("sum_d0xy", &sum_d0xy_sig, "sum_d0xy/F");        
0373   tree_sig->Branch("dca_12", &dca_12_sig, "dca_12/F");
0374   tree_sig->Branch("dca_D0", &dca_D0_sig, "dca_D0/F");
0375   tree_sig->Branch("pt_D0", &pt_D0_sig, "pt_D0/F");  
0376   tree_sig->Branch("y_D0", &y_D0_sig, "y_D0/F");  
0377   tree_sig->Branch("mass_D0", &mass_D0_sig, "mass_D0/F");              
0378   tree_sig->Branch("decay_length", &decay_length_sig, "decay_length/F");   
0379   tree_sig->Branch("costheta", &costheta_sig, "costheta/F"); 
0380   tree_sig->Branch("costheta_xy", &costhetaxy_sig, "costheta_xy/F"); 
0381   tree_sig->Branch("sigma_vtx", &sigma_vtx_sig, "sigma_vtx/F"); 
0382   tree_sig->Branch("mult", &mult_sig, "mult/F"); 
0383   tree_sig->Branch("signif_d0xy_pi", &signif_d0xy_pi_sig, "signif_d0xy_pi/F");               
0384   tree_sig->Branch("signif_d0xy_k", &signif_d0xy_k_sig, "signif_d0xy_k/F"); 
0385   tree_sig->Branch("chi2_dca", &chi2_dca_sig, "chi2_dca/F"); 
0386   
0387   TFile *file_bkg = new TFile("BkgD0.root", "RECREATE");
0388   TTree *tree_bkg = new TTree("treeMLBkg", "treeMLBkg"); 
0389   
0390   // Define variables to store in the Ntuple
0391   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;
0392   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;
0393   // Link the variables to the TTree branches
0394   tree_bkg->Branch("d0_pi", &d0_pi_bkg, "d0_pi/F");
0395   tree_bkg->Branch("d0_k", &d0_k_bkg, "d0_k/F");
0396   tree_bkg->Branch("d0xy_pi", &d0xy_pi_bkg, "d0xy_pi/F");
0397   tree_bkg->Branch("d0xy_k", &d0xy_k_bkg, "d0xy_k/F");
0398   tree_bkg->Branch("sum_d0xy", &sum_d0xy_bkg, "sum_d0xy/F");            
0399   tree_bkg->Branch("dca_12", &dca_12_bkg, "dca_12/F");
0400   tree_bkg->Branch("dca_D0", &dca_D0_bkg, "dca_D0/F");
0401   tree_bkg->Branch("pt_D0", &pt_D0_bkg, "pt_D0/F");  
0402   tree_bkg->Branch("y_D0", &y_D0_bkg, "y_D0/F");  
0403   tree_bkg->Branch("mass_D0", &mass_D0_bkg, "mass_D0/F");              
0404   tree_bkg->Branch("decay_length", &decay_length_bkg, "decay_length/F");   
0405   tree_bkg->Branch("costheta", &costheta_bkg, "costheta/F");  
0406   tree_bkg->Branch("costheta_xy", &costhetaxy_bkg, "costheta_xy/F"); 
0407   tree_bkg->Branch("sigma_vtx", &sigma_vtx_bkg, "sigma_vtx/F"); 
0408   tree_bkg->Branch("mult", &mult_bkg, "mult/F");
0409   tree_bkg->Branch("signif_d0xy_pi", &signif_d0xy_pi_bkg, "signif_d0xy_pi/F");  
0410   tree_bkg->Branch("signif_d0xy_k", &signif_d0xy_k_bkg, "signif_d0xy_k/F");  
0411   tree_bkg->Branch("chi2_dca", &chi2_dca_bkg, "chi2_dca/F");    
0412   
0413   int nevents = 0;
0414   int mult_charged = 0;
0415   while(treereader.Next())
0416     {
0417       if(nevents%1000==0) printf("\n[i] New event %d\n",nevents);
0418       //if(nevents==20) break;
0419 
0420       // find MC primary vertex
0421       int nMCPart = mcPartMass.GetSize();
0422 
0423       TVector3 vertex_mc(-999., -999., -999.);
0424       for(int imc=0; imc<nMCPart; imc++)
0425     {
0426       if(mcPartGenStatus[imc] == 4 && mcPartPdg[imc] == 11)
0427         {
0428           vertex_mc.SetXYZ(mcEndPointX[imc], mcEndPointY[imc], mcEndPointZ[imc]);
0429           break;
0430         }
0431     }
0432       hEventStat->Fill(0.5);
0433       hMcVtxX->Fill(vertex_mc.x());
0434       hMcVtxY->Fill(vertex_mc.y());
0435       hMcVtxZ->Fill(vertex_mc.z());
0436 
0437       // get RC primary vertex
0438       TVector3 vertex_rc(-999., -999., -999.);
0439       TVector3 err_vertex_rc(-999., -999., -999.);
0440       if(prim_vtx_index.GetSize()>0)
0441     {
0442       int rc_vtx_index = prim_vtx_index[0];
0443       vertex_rc.SetXYZ(CTVx[rc_vtx_index], CTVy[rc_vtx_index], CTVz[rc_vtx_index]);
0444       err_vertex_rc.SetXYZ(sqrt(CTVerr_xx[rc_vtx_index]), sqrt(CTVerr_yy[rc_vtx_index]), sqrt(CTVerr_zz[rc_vtx_index]));
0445     }
0446     hPullVtxX->Fill((vertex_rc.x()-vertex_mc.x())/err_vertex_rc.x()); 
0447     hPullVtxY->Fill((vertex_rc.y()-vertex_mc.y())/err_vertex_rc.y()); 
0448     hPullVtxZ->Fill((vertex_rc.z()-vertex_mc.z())/err_vertex_rc.z());   
0449 
0450       // map MC and RC particles
0451       int nAssoc = assocChRecID.GetSize();
0452       map<int, int> assoc_map_to_rc;
0453       map<int, int> assoc_map_to_mc;
0454       
0455       for(unsigned int rc_index=0; rc_index<rcMomPx.GetSize(); rc_index++)
0456     {
0457       // loop over the association to find the matched MC particle
0458       // with largest weight
0459       double max_weight = 0;
0460       int matched_mc_index = -1;
0461       for(int j=0; j<nAssoc; j++)
0462         {
0463           if(assocChRecID[j] != rc_index) continue;
0464           if(assocWeight[j] > max_weight)
0465         {
0466           max_weight = assocWeight[j];
0467           matched_mc_index = assocChSimID[j];
0468         }
0469         }
0470 
0471       // build the map
0472       assoc_map_to_rc[matched_mc_index] = rc_index;
0473       assoc_map_to_mc[rc_index] = matched_mc_index;
0474     }
0475 
0476       // Loop over primary particles
0477       int nMcPart = 0;
0478       for(int imc=0; imc<nMCPart; imc++)
0479     {
0480       if(mcPartGenStatus[imc] == 1 && mcPartCharge[imc] != 0)
0481         {
0482           double dist = sqrt( pow(mcPartVx[imc]-vertex_mc.x(),2) + pow(mcPartVy[imc]-vertex_mc.y(),2) + pow(mcPartVz[imc]-vertex_mc.z(),2));      
0483           if(dist < 1e-4)
0484         {
0485           // count charged particles within |eta| < 3.5
0486           TVector3 mc_mom(mcMomPx[imc], mcMomPy[imc], mcMomPz[imc]);
0487           double mcEta = mc_mom.PseudoRapidity();
0488           if(fabs(mcEta) < 3.5) nMcPart++;
0489         }
0490         }
0491     }
0492       mult_charged = nMcPart;
0493       hMcMult->Fill(nMcPart);
0494       
0495       for(int imc=0; imc<nMCPart; imc++)
0496     {
0497       if(mcPartGenStatus[imc] == 1 && mcPartCharge[imc] != 0)
0498         {
0499           double dist = sqrt( pow(mcPartVx[imc]-vertex_mc.x(),2) + pow(mcPartVy[imc]-vertex_mc.y(),2) + pow(mcPartVz[imc]-vertex_mc.z(),2));      
0500           if(dist < 1e-4)
0501         {         
0502           // check if the MC particle is reconstructed
0503           int rc_index = -1;
0504           if(assoc_map_to_rc.find(imc) != assoc_map_to_rc.end()) rc_index = assoc_map_to_rc[imc];
0505 
0506           if(rc_index>=0)
0507             {
0508               TVector3 dcaToVtx = getDcaToVtx(rc_index, vertex_rc);
0509               
0510               int ip = -1;
0511               if(fabs(mcPartPdg[imc]) == 211) ip = 0;
0512               if(fabs(mcPartPdg[imc]) == 321) ip = 1;
0513               if(fabs(mcPartPdg[imc]) == 2212) ip = 2;
0514               if(ip>=0)
0515             {
0516               TVector3 mom(rcMomPx[rc_index], rcMomPy[rc_index], rcMomPz[rc_index]);
0517               if(ip<2)
0518                 {
0519                   hRcPrimPartLocaToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.Perp());
0520                   hRcPrimPartLocbToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.z());
0521                 }
0522 
0523               double fill1[] = {mom.Pt(), mom.Eta(), dcaToVtx.x(), nMcPart*1.};
0524               double fill2[] = {mom.Pt(), mom.Eta(), dcaToVtx.y(), nMcPart*1.};
0525               double fill3[] = {mom.Pt(), mom.Eta(), dcaToVtx.z(), nMcPart*1.};
0526               hPrimTrkDcaToRCVtx[ip][0]->Fill(fill1);
0527               hPrimTrkDcaToRCVtx[ip][1]->Fill(fill2);
0528               hPrimTrkDcaToRCVtx[ip][2]->Fill(fill3);
0529             }
0530             }
0531         }
0532             }
0533     }
0534 
0535       // look for D0
0536       bool hasD0 = false;
0537       vector<int> mc_index_D0_pi;
0538       vector<int> mc_index_D0_k;
0539       mc_index_D0_pi.clear();
0540       mc_index_D0_k.clear();
0541       
0542       for(int imc=0; imc<nMCPart; imc++)
0543     {
0544       if(fabs(mcPartPdg[imc]) != 421) continue;
0545       hEventStat->Fill(1.5);
0546       
0547       int nDuaghters = mcPartDaughter_end[imc]-mcPartDaughter_begin[imc];
0548       if(nDuaghters!=2) continue;
0549 
0550       // find D0 that decay into pi+K
0551       bool is_pik_decay = false;      
0552       int daug_index_1 = mcPartDaughter_index[mcPartDaughter_begin[imc]];
0553       int daug_index_2 = mcPartDaughter_index[mcPartDaughter_begin[imc]+1];
0554       int daug_pdg_1 = mcPartPdg[daug_index_1];
0555       int daug_pdg_2 = mcPartPdg[daug_index_2];
0556       if( (fabs(daug_pdg_1)==321 && fabs(daug_pdg_2)==211) || (fabs(daug_pdg_1)==211 && fabs(daug_pdg_2)==321) )
0557         {
0558           is_pik_decay = true;
0559         }
0560       if(!is_pik_decay) continue;
0561       if(fabs(daug_pdg_1)==211)
0562         {
0563           mc_index_D0_pi.push_back(daug_index_1);
0564           mc_index_D0_k.push_back(daug_index_2);
0565         }
0566       else
0567         {
0568           mc_index_D0_pi.push_back(daug_index_2);
0569           mc_index_D0_k.push_back(daug_index_1);
0570         }
0571       hasD0 = true;
0572       hEventStat->Fill(2.5);
0573 
0574       // D0 kinematics
0575       TLorentzVector mc_mom_vec;
0576       mc_mom_vec.SetXYZM(mcMomPx[imc], mcMomPy[imc], mcMomPz[imc], mcPartMass[imc]);
0577       double mcRap = mc_mom_vec.Rapidity();
0578       double mcPt = mc_mom_vec.Pt();
0579       hMCD0PtRap->Fill(mcRap, mcPt);
0580 
0581       // decay dauther kinematics
0582       for(int ip = 0; ip<2; ip++)
0583         {
0584           int mc_part_index;
0585           if(ip==0) mc_part_index = mc_index_D0_pi[mc_index_D0_pi.size()-1];
0586           if(ip==1) mc_part_index = mc_index_D0_k[mc_index_D0_k.size()-1];
0587           
0588           TLorentzVector mc_part_vec;
0589           mc_part_vec.SetXYZM(mcMomPx[mc_part_index], mcMomPy[mc_part_index], mcMomPz[mc_part_index], mcPartMass[mc_part_index]);
0590           if(ip==0) hMcPiPtEta->Fill(mc_part_vec.Eta(), mc_part_vec.Pt());
0591           if(ip==1) hMcKPtEta->Fill(mc_part_vec.Eta(), mc_part_vec.Pt());
0592           
0593           int rc_part_index = -1;
0594           if(assoc_map_to_rc.find(mc_part_index) != assoc_map_to_rc.end()) rc_part_index = assoc_map_to_rc[mc_part_index];
0595           if(rc_part_index>=0)
0596         {
0597           TVector3 dcaToVtx = getDcaToVtx(rc_part_index, vertex_rc);
0598           
0599           TVector3 mom(rcMomPx[rc_part_index], rcMomPy[rc_part_index], rcMomPz[rc_part_index]);
0600           hRcSecPartLocaToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.Pt());
0601           hRcSecPartLocbToRCVtx[ip]->Fill(mom.Pt(), mom.Eta(), dcaToVtx.z());
0602           
0603           //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]);
0604         }
0605         }
0606     }
0607 
0608       // Get reconstructed pions and kaons
0609       hNRecoVtx->Fill(CTVx.GetSize());
0610       const int pid_mode = 1; // 0 - truth; 1 - realistic
0611       vector<unsigned int> pi_index;
0612       vector<unsigned int> k_index;
0613       pi_index.clear();
0614       k_index.clear();
0615       for(unsigned int rc_index=0; rc_index<rcMomPx.GetSize(); rc_index++)
0616     {     
0617       if(pid_mode==0)
0618         {
0619           int iSimPartID = -1;
0620           if(assoc_map_to_mc.find(rc_index) != assoc_map_to_mc.end()) iSimPartID = assoc_map_to_mc[rc_index];
0621           if(iSimPartID>=0)
0622         {
0623           if(fabs(mcPartPdg[iSimPartID]) == 211) pi_index.push_back(rc_index);
0624           if(fabs(mcPartPdg[iSimPartID]) == 321) k_index.push_back(rc_index);
0625         }
0626         }
0627       else if(pid_mode==1)
0628         {
0629           if(fabs(rcPdg[rc_index]) == 211) pi_index.push_back(rc_index);
0630           if(fabs(rcPdg[rc_index]) == 321) k_index.push_back(rc_index);
0631         }
0632     }
0633 
0634       // pair pion and kaon
0635       for(unsigned int i=0; i<pi_index.size(); i++)
0636     {
0637       TVector3 dcaToVtx = getDcaToVtx(pi_index[i], vertex_rc);
0638        std::array<float, 21>& cov_pion = rcTrkCov->At(pi_index[i]);
0639      int q_pion = rcCharge[pi_index[i]];
0640       for(unsigned int j=0; j<k_index.size(); j++)
0641         {
0642           TVector3 dcaToVtx2 = getDcaToVtx(k_index[j], vertex_rc);
0643           std::array<float, 21>& cov_kaon = rcTrkCov->At(k_index[j]); 
0644         int q_kaon = rcCharge[k_index[j]];
0645           if(rcCharge[pi_index[i]]*rcCharge[k_index[j]]<0)
0646         {
0647           //printf("[i] Check pair (%d, %d)\n", pi_index[i], k_index[j]);
0648           // -- only look at unlike-sign pi+k pair
0649           bool is_D0_pik = false;
0650           int mc_index_pi = -1, mc_index_k = -1;
0651           if(assoc_map_to_mc.find(pi_index[i]) != assoc_map_to_mc.end()) mc_index_pi = assoc_map_to_mc[pi_index[i]];
0652           if(assoc_map_to_mc.find(k_index[j])  != assoc_map_to_mc.end()) mc_index_k  = assoc_map_to_mc[k_index[j]];
0653 
0654           for(unsigned int k=0; k<mc_index_D0_pi.size(); k++)
0655             {
0656               if(mc_index_pi==mc_index_D0_pi[k] && mc_index_k==mc_index_D0_k[k])
0657             {
0658               is_D0_pik = true;
0659               break;
0660             }
0661             }
0662 
0663           float dcaDaughters, cosTheta, decayLength, V0DcaToVtx, cosTheta_xy, sigma_vtx;
0664           TVector3 decayVertex; 
0665           double chi2_ndf;
0666           double err_Par[5];  // or whatever size is appropriate
0667       //double* err_Par = errParArray;
0668           TLorentzVector parent = getPairParent(pi_index[i], k_index[j], vertex_rc, dcaDaughters, cosTheta, cosTheta_xy, decayLength, V0DcaToVtx, sigma_vtx,decayVertex,chi2_ndf, err_Par);
0669           hchi2_vtx->Fill(chi2_ndf);
0670                   
0671           if(is_D0_pik)
0672             {
0673            TVector3 MCVertex_Kaon(mcPartVx[mc_index_k], mcPartVy[mc_index_k], mcPartVz[mc_index_k]);
0674            TVector3 MCVertex_Pion(mcPartVx[mc_index_pi], mcPartVy[mc_index_pi], mcPartVz[mc_index_pi]);
0675            
0676           // printf("Signal MC Vertex Kaon = (%f, %f, %f)\n",MCVertex_Kaon.X(), MCVertex_Kaon.Y(), MCVertex_Kaon.Z());      
0677           // printf("Signal MC Vertex Pion = (%f, %f, %f)\n",MCVertex_Pion.X(), MCVertex_Pion.Y(), MCVertex_Pion.Z());
0678               hRes_SVx_Helixfit->Fill((decayVertex.X()-MCVertex_Kaon.X()*0.1)*10);
0679               hRes_SVy_Helixfit->Fill((decayVertex.Y()-MCVertex_Kaon.Y()*0.1)*10);
0680               hRes_SVz_Helixfit->Fill((decayVertex.Z()-MCVertex_Kaon.Z()*0.1)*10);  
0681               
0682               hRes_SVx_Helixfit_pull->Fill(((decayVertex.X()-MCVertex_Kaon.X()*0.1))/err_Par[0]);
0683               hRes_SVy_Helixfit_pull->Fill(((decayVertex.Y()-MCVertex_Kaon.Y()*0.1))/err_Par[1]);
0684               hRes_SVz_Helixfit_pull->Fill(((decayVertex.Z()-MCVertex_Kaon.Z()*0.1))/err_Par[2]);   
0685               
0686               hchi2_vtx_sig->Fill(chi2_ndf);              
0687                                
0688               hEventStat->Fill(3.5);
0689               if (q_kaon == -1 && q_pion == 1)
0690               hEventStat->Fill(4.5);   // D0
0691               else if (q_kaon == 1 && q_pion == -1)
0692               hEventStat->Fill(5.5);   // D0bar
0693               h3PairDca12[0]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters);
0694               h3PairCosTheta[0]->Fill(parent.Pt(), parent.Rapidity(), cosTheta);
0695               h3PairDca[0]->Fill(parent.Pt(), parent.Rapidity(), V0DcaToVtx);
0696               h3PairDecayLength[0]->Fill(parent.Pt(), parent.Rapidity(), decayLength);
0697               //printf("Signal: dca12 = %2.4f, cosTheta = %2.4f, D0dca = %2.4f, decay = %2.4f\n", dcaDaughters, cosTheta, V0DcaToVtx, decayLength);
0698               h3InvMass[0][0]->Fill(parent.Pt(), parent.Rapidity(), parent.M());
0699                
0700            // Toplogical Variables for Signal
0701               d0_pi_sig = dcaToVtx.Mag();
0702               d0_k_sig = dcaToVtx2.Mag();
0703               d0xy_pi_sig = dcaToVtx.Perp();
0704               signif_d0xy_pi_sig = d0xy_pi_sig/sqrt(cov_pion[0]);
0705               d0xy_k_sig = dcaToVtx2.Perp();
0706               signif_d0xy_k_sig = d0xy_k_sig/sqrt(cov_kaon[0]);
0707               sum_d0xy_sig = sqrt(d0xy_pi_sig*d0xy_pi_sig+d0xy_k_sig*d0xy_k_sig);             
0708               dca_12_sig = dcaDaughters;
0709                 dca_D0_sig = V0DcaToVtx;
0710                 decay_length_sig = decayLength;
0711                 costheta_sig = cosTheta;
0712                 costhetaxy_sig = cosTheta_xy;
0713                 pt_D0_sig = parent.Pt();
0714                 y_D0_sig = parent.Rapidity();
0715                 mass_D0_sig = parent.M();
0716                 sigma_vtx_sig = sigma_vtx;
0717                 mult_sig = mult_charged;
0718                 chi2_dca_sig = chi2_ndf;
0719                 tree_sig->Fill();
0720             }
0721           else
0722             {
0723           hEventStat->Fill(6.5);   // Bkg
0724           hchi2_vtx_bkg->Fill(chi2_ndf);                    
0725               h3PairDca12[1]->Fill(parent.Pt(), parent.Rapidity(), dcaDaughters);
0726               h3PairCosTheta[1]->Fill(parent.Pt(), parent.Rapidity(), cosTheta);
0727               h3PairDca[1]->Fill(parent.Pt(), parent.Rapidity(), V0DcaToVtx);
0728               h3PairDecayLength[1]->Fill(parent.Pt(), parent.Rapidity(), decayLength);
0729               
0730               //printf("Bkg: dca12 = %2.4f, cosTheta = %2.4f, D0dca = %2.4f, decay = %2.4f\n", dcaDaughters, cosTheta, V0DcaToVtx, decayLength);
0731               h3InvMass[1][0]->Fill(parent.Pt(), parent.Rapidity(), parent.M());
0732 
0733               // Toplogical Variables for Bkg
0734               d0_pi_bkg = dcaToVtx.Mag();
0735               d0_k_bkg = dcaToVtx2.Mag();
0736               d0xy_pi_bkg = dcaToVtx.Perp();
0737               signif_d0xy_pi_bkg = d0xy_pi_bkg/sqrt(cov_pion[0]);
0738               d0xy_k_bkg = dcaToVtx2.Perp();
0739               signif_d0xy_k_bkg = d0xy_k_bkg/sqrt(cov_kaon[0]);
0740               sum_d0xy_bkg = sqrt(d0xy_pi_bkg*d0xy_pi_bkg+d0xy_k_bkg*d0xy_k_bkg);             
0741               dca_12_bkg = dcaDaughters;
0742                 dca_D0_bkg = V0DcaToVtx;
0743                 decay_length_bkg = decayLength;
0744                 costheta_bkg = cosTheta;
0745                 costhetaxy_bkg = cosTheta_xy;
0746                 pt_D0_bkg = parent.Pt();
0747                 y_D0_bkg = parent.Rapidity();
0748                 mass_D0_bkg = parent.M();
0749                 sigma_vtx_bkg = sigma_vtx;
0750                 mult_bkg = mult_charged;
0751                 chi2_dca_bkg = chi2_ndf;                
0752                 tree_bkg->Fill();
0753             }
0754 
0755           if(dcaToVtx.Perp() >= 0.02 && dcaToVtx2.Perp() >= 0.02 &&
0756              dcaDaughters < 0.07 && cosTheta > 0.95 && decayLength > 0.05 && V0DcaToVtx < 0.1)
0757             {
0758               if(is_D0_pik)
0759             {
0760               h3InvMass[0][1]->Fill(parent.Pt(), parent.Rapidity(), parent.M());
0761             }
0762               else
0763             {
0764               h3InvMass[1][1]->Fill(parent.Pt(), parent.Rapidity(), parent.M());
0765             }
0766             }
0767         }
0768         }
0769     }
0770           
0771       nevents++;
0772     }
0773 
0774   file_signal->cd();  
0775   tree_sig->Write();
0776   file_signal->Close();  
0777   
0778   file_bkg->cd();  
0779   tree_bkg->Write();
0780   file_bkg->Close();
0781 
0782   TFile *outfile = new TFile(outname.Data(), "recreate");
0783 
0784   hEventStat->Write();
0785   hMcMult->Write();
0786   hMcVtxX->Write();
0787   hMcVtxY->Write();
0788   hMcVtxZ->Write();
0789   
0790   hPullVtxX->Write();
0791   hPullVtxY->Write();  
0792   hPullVtxZ->Write();
0793   
0794   hRes_SVx_Helixfit->Write();
0795   hRes_SVy_Helixfit->Write();
0796   hRes_SVz_Helixfit->Write(); 
0797   
0798   hchi2_vtx->Write();
0799   hchi2_vtx_sig->Write();
0800   hchi2_vtx_bkg->Write();  
0801   hRes_SVx_Helixfit_pull->Write();
0802   hRes_SVy_Helixfit_pull->Write(); 
0803   hRes_SVz_Helixfit_pull->Write(); 
0804     
0805   hD0DecayVxVy->Write();
0806   hD0DecayVrVz->Write();
0807   
0808   hMCD0PtRap->Write();
0809 
0810   hMcPiPtEta->Write();
0811   hMcPiPtEtaReco->Write();
0812   hMcKPtEta->Write();
0813   hMcKPtEtaReco->Write();
0814   
0815   hNRecoVtx->Write();
0816 
0817   for(int ip=0; ip<2; ip++)
0818     {
0819       hRcSecPartLocaToRCVtx[ip]->Write();
0820       hRcSecPartLocbToRCVtx[ip]->Write();
0821       hRcPrimPartLocaToRCVtx[ip]->Write();
0822       hRcPrimPartLocbToRCVtx[ip]->Write();
0823     }
0824 
0825   for(int i=0; i<3; i++)
0826     {
0827       for(int j=0; j<3; j++)
0828     {
0829       hPrimTrkDcaToRCVtx[i][j]->Write();
0830     }
0831     }
0832 
0833   for(int i=0; i<2; i++)
0834     {
0835       h3PairDca12[i]->Write();
0836       h3PairCosTheta[i]->Write();
0837       h3PairDca[i]->Write();
0838       h3PairDecayLength[i]->Write();
0839     }
0840 
0841   for(int i=0; i<2; i++)
0842     {
0843       for(int j=0; j<2; j++)
0844     {
0845       h3InvMass[i][j]->Write();
0846     }
0847     }
0848   
0849   
0850   outfile->Close();
0851 }
0852 
0853 //======================================
0854 TVector3 getDcaToVtx(const int index, TVector3 vtx)
0855 {
0856   TVector3 pos(rcTrkLoca2->At(index) * sin(rcTrkPhi2->At(index)) * -1 * millimeter, rcTrkLoca2->At(index) * cos(rcTrkPhi2->At(index)) * millimeter, rcTrkLocb2->At(index) * millimeter);
0857   TVector3 mom(rcMomPx2->At(index), rcMomPy2->At(index), rcMomPz2->At(index));
0858    
0859   StPhysicalHelix pHelix(mom, pos, bField * tesla, rcCharge2->At(index));
0860 
0861   TVector3 vtx_tmp;
0862   vtx_tmp.SetXYZ(vtx.x()*millimeter, vtx.y()*millimeter, vtx.z()*millimeter);
0863   
0864   pHelix.moveOrigin(pHelix.pathLength(vtx_tmp));
0865   TVector3 dcaToVtx = pHelix.origin() - vtx_tmp;
0866 
0867   dcaToVtx.SetXYZ(dcaToVtx.x()/millimeter, dcaToVtx.y()/millimeter, dcaToVtx.z()/millimeter);
0868   
0869   return dcaToVtx;
0870 }
0871 
0872 //======================================
0873 TLorentzVector getPairParent(const int index1, const int index2, TVector3 vtx,
0874                  float &dcaDaughters, float &cosTheta, float &cosTheta_xy, float &decayLength, float &V0DcaToVtx, float &sigma_vtx, TVector3 &decayVertex, double &chi2_ndf, double * parFitErr)
0875 {
0876   // -- get helix
0877   TVector3 pos1(rcTrkLoca2->At(index1) * sin(rcTrkPhi2->At(index1)) * -1 * millimeter, rcTrkLoca2->At(index1) * cos(rcTrkPhi2->At(index1)) * millimeter, rcTrkLocb2->At(index1) * millimeter);
0878   TVector3 mom1(rcMomPx2->At(index1), rcMomPy2->At(index1), rcMomPz2->At(index1));
0879 
0880   TVector3 pos2(rcTrkLoca2->At(index2) * sin(rcTrkPhi2->At(index2)) * -1 * millimeter, rcTrkLoca2->At(index2) * cos(rcTrkPhi2->At(index2)) * millimeter, rcTrkLocb2->At(index2) * millimeter);
0881   TVector3 mom2(rcMomPx2->At(index2), rcMomPy2->At(index2), rcMomPz2->At(index2));
0882 
0883   float charge1 = rcCharge2->At(index1);
0884   float charge2 = rcCharge2->At(index2);
0885   
0886   StPhysicalHelix p1Helix(mom1, pos1, bField * tesla, charge1);
0887   StPhysicalHelix p2Helix(mom2, pos2, bField * tesla, charge2);
0888 
0889   TVector3 vtx_tmp;
0890   vtx_tmp.SetXYZ(vtx.x()*millimeter, vtx.y()*millimeter, vtx.z()*millimeter);
0891   
0892   double s1, s2;
0893   getDecayVertex_Chi2fit(index1, index2, s1, s2, decayVertex, chi2_ndf, parFitErr);
0894  
0895   TVector3 const p1AtDcaToP2 = p1Helix.at(s1);
0896   TVector3 const p2AtDcaToP1 = p2Helix.at(s2);
0897   // printf("p1AtDcaToP2 origin = (%2.4f, %2.4f, %2.4f)\n", p1AtDcaToP2.x(), p1AtDcaToP2.y(), p1AtDcaToP2.z());
0898   // printf("p2AtDcaToP1 origin = (%2.4f, %2.4f, %2.4f)\n", p2AtDcaToP1.x(), p2AtDcaToP1.y(), p2AtDcaToP1.z());
0899   
0900   // -- calculate DCA of particle1 to particle2 at their DCA
0901   dcaDaughters = (p1AtDcaToP2 - p2AtDcaToP1).Mag()/millimeter;
0902     
0903   // -- calculate Lorentz vector of particle1-particle2 pair
0904   TVector3 const p1MomAtDca = p1Helix.momentumAt(s1,  bField * tesla);
0905   TVector3 const p2MomAtDca = p2Helix.momentumAt(s2, bField * tesla);
0906   
0907   TLorentzVector p1FourMom(p1MomAtDca, sqrt(p1MomAtDca.Mag2()+gPionMass*gPionMass));
0908   TLorentzVector p2FourMom(p2MomAtDca, sqrt(p2MomAtDca.Mag2()+gKaonMass*gKaonMass));
0909   
0910   TLorentzVector parent = p1FourMom + p2FourMom;
0911 
0912   // -- calculate decay vertex (secondary or tertiary)
0913  // decayVertex = (p1AtDcaToP2 + p2AtDcaToP1) * 0.5 ;
0914   sigma_vtx = sqrt((p1AtDcaToP2-decayVertex).Mag2()+(p2AtDcaToP1-decayVertex).Mag2())/millimeter;
0915     
0916   // -- calculate pointing angle and decay length with respect to primary vertex
0917   //    if decay vertex is a tertiary vertex
0918   //    -> only rough estimate -> needs to be updated after secondary vertex is found
0919   TVector3 vtxToV0 = decayVertex - vtx_tmp;
0920   TVector3 vtxToV0_xy(vtxToV0.x(), vtxToV0.y(), 0.);
0921   TVector3 parent_xy(parent.Vect().x(),parent.Vect().y(),0.);
0922   float pointingAngle = vtxToV0.Angle(parent.Vect());
0923   float pointingAngle_xy = vtxToV0_xy.Angle(parent_xy);
0924   cosTheta = std::cos(pointingAngle);
0925   cosTheta_xy = std::cos(pointingAngle_xy);  
0926   decayLength = vtxToV0.Mag()/millimeter;
0927 
0928   // -- calculate V0 DCA to primary vertex
0929   V0DcaToVtx = decayLength * std::sin(pointingAngle);
0930     
0931   //TVector3 dcaToVtx = getDcaToVtx(parent.Vect(), decayVertex, 0, vtx);
0932   //V0DcaToVtx = dcaToVtx.Mag();
0933   return parent;
0934 }
0935 
0936 
0937 
0938 void getDecayVertex_Chi2fit(const int index1, const int index2, double &s1, double &s2, TVector3 &vertex, double &chi2_ndf, double *parFitErr)
0939 {
0940     TVector3 pos1(rcTrkLoca2->At(index1) * sin(rcTrkPhi2->At(index1)) * -1 * millimeter,
0941                   rcTrkLoca2->At(index1) * cos(rcTrkPhi2->At(index1)) * millimeter,
0942                   rcTrkLocb2->At(index1) * millimeter);
0943 
0944     TVector3 mom1(rcMomPx2->At(index1), rcMomPy2->At(index1), rcMomPz2->At(index1));
0945 
0946     TVector3 pos2(rcTrkLoca2->At(index2) * sin(rcTrkPhi2->At(index2)) * -1 * millimeter,
0947                   rcTrkLoca2->At(index2) * cos(rcTrkPhi2->At(index2)) * millimeter,
0948                   rcTrkLocb2->At(index2) * millimeter);
0949 
0950     TVector3 mom2(rcMomPx2->At(index2), rcMomPy2->At(index2), rcMomPz2->At(index2));
0951 
0952 
0953     float charge1 = rcCharge2->At(index1);
0954     float charge2 = rcCharge2->At(index2);
0955 
0956 
0957     StPhysicalHelix helix1(mom1, pos1, bField * tesla, charge1);
0958     StPhysicalHelix helix2(mom2, pos2, bField * tesla, charge2);
0959    
0960     std::array<float, 21>& cov_track1 = rcTrkCov->At(index1);
0961     std::array<float, 21>& cov_track2 = rcTrkCov->At(index2);  
0962     
0963     pair<double, double> const ss = helix1.pathLengths(helix2);
0964     TVector3 const p1_init = helix1.at(ss.first);
0965     TVector3 const p2_init = helix2.at(ss.second); 
0966     TVector3 const mid_point = 0.5*(p1_init+p2_init); 
0967 
0968       // Perform Minimization
0969     const Int_t nPar = 5;
0970     Chi2Minimization d2Function(helix1,helix2,cov_track1,cov_track2);
0971     ROOT::Math::Functor fcn(d2Function,nPar); // 5 parameters
0972     ROOT::Fit::Fitter fitter;
0973 
0974     double pStart[nPar] = {mid_point.X(),mid_point.Y(),mid_point.Z(),ss.first,ss.second};
0975     fitter.SetFCN(fcn, pStart,nPar,1);
0976     
0977     fitter.Config().ParSettings(0).SetName("x0");
0978     fitter.Config().ParSettings(0).SetStepSize(0.01);
0979    // fitter.Config().ParSettings(0).SetLimits(-1., 1.);    
0980     // No limits for x, y, z
0981 
0982     fitter.Config().ParSettings(1).SetName("y0");
0983     fitter.Config().ParSettings(1).SetStepSize(0.01);
0984     //fitter.Config().ParSettings(1).SetLimits(-1., 1.);
0985 
0986     fitter.Config().ParSettings(2).SetName("z0");
0987     fitter.Config().ParSettings(2).SetStepSize(0.01);
0988    // fitter.Config().ParSettings(2).SetLimits(-10., 10.);    
0989 
0990     fitter.Config().ParSettings(3).SetName("s1");
0991     fitter.Config().ParSettings(3).SetValue(0.0);
0992     fitter.Config().ParSettings(3).SetStepSize(0.01);
0993     //fitter.Config().ParSettings(3).SetLimits(-1., 1.);
0994 
0995     fitter.Config().ParSettings(4).SetName("s2");
0996     fitter.Config().ParSettings(4).SetValue(0.0);
0997     fitter.Config().ParSettings(4).SetStepSize(0.01);
0998     //fitter.Config().ParSettings(4).SetLimits(-1., 1.);      
0999     
1000     fitter.Config().MinimizerOptions().SetMaxIterations(10000); 
1001     // do the fit 
1002 
1003     Bool_t ok = fitter.FitFCN();
1004     if (!ok) Error("Fitting","Fitting failed");
1005     const ROOT::Fit::FitResult & result = fitter.Result();
1006    // double chi2 = fitter.Result().Chi2();
1007     double ndf = 2*3-nPar;
1008     chi2_ndf = fitter.Result().MinFcnValue()/ndf;  // Minimum value of your function
1009     
1010     int status = fitter.Result().Status();
1011     if (status>0 ) {printf("Fit Failed!!!!\n");}
1012    // if (status>0 || chi2_ndf>10. ) return;
1013    // cout <<"\033[1;31m Fit Result Chi2:\033[0m"<<chi2_ndf<<endl;
1014    // result.Print(std::cout);
1015    
1016    // Get the covariance matrix
1017  //  TMatrixDSym covMatrix(5);
1018   // result.GetCovarianceMatrix(covMatrix); // Matrix for the parameter errors
1019   // covMatrix.Print();
1020    
1021    const double * parFit = result.GetParams();
1022    const double *FitErr = result.GetErrors();
1023     
1024    for (int i = 0; i < nPar; ++i) parFitErr[i] = FitErr[i];
1025    vertex.SetXYZ(parFit[0], parFit[1], parFit[2]);
1026    s1 = parFit[3]; s2 = parFit[4];  
1027 }
1028 
1029
[ 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!