JetsAndHF/SecondaryVertex_Chi2/Check_Origin_D0meson.C

0001 // Code to check the orgin of D0 (prompt vs non-prompt) at the simulation level.
0002 // Shyam Kumar; shyam.kumar@ba.infn.it; shyam055119@gmail.com
0003
0004 int GetD0OriginFlavor(int imc, const TTreeReaderArray<int> &mcPartPdg, const TTreeReaderArray<unsigned int> &mcPartParent_begin, const TTreeReaderArray<unsigned int> &mcPartParent_end, const TTreeReaderArray<int> &mcPartParent_index);
0005 void Check_Origin_D0meson(TString filelist = "test.list"){
0006      gStyle->SetOptStat(0);
0007     // Read all files
0008     int nfiles = 0;
0009      char filename[512];
0010     TChain *chain = new TChain("events");
0011
0012     ifstream *inputstream = new ifstream;
0013      inputstream->open(filelist.Data());
0014      if(!inputstream) printf("[Error] Cannot open file list: %s\n", filelist.Data());
0015
0016     while (inputstream->good())
0017       {
0018       inputstream->getline(filename, 512);
0019       if (inputstream->good())
0020       {
0021         TFile *ftmp = TFile::Open(filename, "READ");
0022         if (!ftmp || !ftmp->IsOpen() || !ftmp->GetNkeys())
0023         {
0024          printf("[e] Skipping bad file: %s\n", filename);
0025          if (ftmp) { ftmp->Close(); delete ftmp; }
0026          continue;
0027         }
0028         cout << "[i] Add " << nfiles << "th file: " << filename << endl;
0029         chain->Add(filename);
0030         nfiles++;
0031
0032         ftmp->Close();
0033         delete ftmp;
0034       }
0035     }
0036
0037     inputstream->close();
0038     printf("[i] Read in %d files with %lld events in total\n", nfiles, chain->GetEntries());
0039     TTreeReader myReader(chain);
0040
0041      // MC Track Properties
0042      TTreeReaderArray<int> mcPartGenStatus = {myReader, "MCParticles.generatorStatus"};
0043      TTreeReaderArray<int> mcPartPdg = {myReader, "MCParticles.PDG"};
0044      TTreeReaderArray<double> mcPartMass = {myReader, "MCParticles.mass"};
0045     TTreeReaderArray<double> mcEndPointX = {myReader, "MCParticles.endpoint.x"};
0046     TTreeReaderArray<double> mcEndPointY = {myReader, "MCParticles.endpoint.y"};
0047     TTreeReaderArray<double> mcEndPointZ = {myReader, "MCParticles.endpoint.z"};
0048     TTreeReaderArray<double> mcMomPx = {myReader, "MCParticles.momentum.x"};
0049      TTreeReaderArray<double> mcMomPy = {myReader, "MCParticles.momentum.y"};
0050      TTreeReaderArray<double> mcMomPz = {myReader, "MCParticles.momentum.z"};
0051      TTreeReaderArray<double> mcPartVx = {myReader, "MCParticles.vertex.x"};
0052      TTreeReaderArray<double> mcPartVy = {myReader, "MCParticles.vertex.y"};
0053      TTreeReaderArray<double> mcPartVz = {myReader, "MCParticles.vertex.z"};
0054      TTreeReaderArray<int> mcPartParent_index = {myReader, "_MCParticles_parents.index"};
0055      TTreeReaderArray<unsigned int> mcPartDaughter_begin = {myReader, "MCParticles.daughters_begin"};
0056      TTreeReaderArray<unsigned int> mcPartParent_begin = {myReader, "MCParticles.parents_begin"};
0057      TTreeReaderArray<unsigned int> mcPartParent_end = {myReader, "MCParticles.parents_end"};
0058      TTreeReaderArray<unsigned int> mcPartDaughter_end = {myReader, "MCParticles.daughters_end"};
0059      TTreeReaderArray<int> mcPartDaughter_index = {myReader, "_MCParticles_daughters.index"};
0060
0061      TCanvas *c = new TCanvas("c","c",2000,1200);
0062      c->SetMargin(0.10, 0.05 ,0.12,0.07);
0063      TH1F *hEventStat = new TH1F("hEventStat", ";;Entries (a.u.)", 6, 0, 6);
0064      hEventStat->GetXaxis()->SetBinLabel(1, "MC events");
0065      hEventStat->GetXaxis()->SetBinLabel(2, "D^{0}");
0066      hEventStat->GetXaxis()->SetBinLabel(3, "NonPrompt D^{0}");
0067      hEventStat->GetXaxis()->SetBinLabel(4, "Prompt D^{0}");
0068      hEventStat->GetXaxis()->SetBinLabel(5, "originHF==#pm 5 (Beauty)");
0069      hEventStat->GetXaxis()->SetBinLabel(6, "originHF==#pm 4 (Charm)");
0070      hEventStat->SetLineWidth(2);
0071      hEventStat->SetMarkerSize(2);
0072      hEventStat->SetTitleSize(0.04,"XY");
0073      hEventStat->SetLabelSize(0.04,"XY");
0074
0075   TFile *fout = new TFile("output_file.root","recreate");
0076   TH2D *hptD0vsdl = new TH2D("hptD0vsdl","hptD0vsdl;pD0 (GeV/c); dl (mm)",100.,0.,20.,100,0.,5.0);
0077   int count = -1;
0078   const double eps = 1e-6; // 1 nm in mm
0079   TDatabasePDG *pdgDB = TDatabasePDG::Instance();
0080   while (myReader.Next())
0081   {
0082    count++;
0083  //  if (count!=nevent) continue;
0084  //  printf("Event No. = %d\n",count);
0085    hEventStat->Fill(0.5);
0086     // find MC primary vertex
0087     int nMCPart = mcPartMass.GetSize();
0088
0089     TVector3 vertex_mc(-999., -999., -999.);
0090     TVector3 SV_mc(-999., -999., -999.);
0091     TVector3 pD0(0.,0.,0.);
0092    for(int imc=0; imc<nMCPart; imc++)
0093       {
0094       if(mcPartGenStatus[imc] == 4 && mcPartPdg[imc] == 11)
0095         {
0096           vertex_mc.SetXYZ(mcEndPointX[imc], mcEndPointY[imc], mcEndPointZ[imc]); // collision vertex
0097           break;
0098         }
0099        }
0100
0101     for(int imc=0; imc<nMCPart; imc++)
0102     {
0103       if(fabs(mcPartPdg[imc]) != 421) continue;
0104        int parent_index_begin = mcPartParent_index[mcPartParent_begin[imc]];
0105        int parent_index_end   = mcPartParent_index[mcPartParent_end[imc]];
0106       hEventStat->Fill(1.5);
0107       TVector3 prodVD0 (mcPartVx[imc],mcPartVy[imc],mcPartVz[imc]);
0108       double dist = (prodVD0-vertex_mc).Mag();
0109       if (dist > eps) hEventStat->Fill(2.5);
0110       else hEventStat->Fill(3.5);
0111       //printf("MCParticle = %d, PDG = %d, Dist = %f \n",imc,mcPartPdg[imc],dist);
0112       // classify origin: 4 = charm, 5 = beauty
0113        int originHF = GetD0OriginFlavor(imc,
0114                                      mcPartPdg,
0115                                      mcPartParent_begin,
0116                                      mcPartParent_end,
0117                                      mcPartParent_index);
0118          // printf("MC D0 (index %d): PDG=%d → originHF = %d\n",imc, mcPartPdg[imc], originHF);
0119        if(fabs(originHF)==5) hEventStat->Fill(4.5);
0120        else if(fabs(originHF)==4) hEventStat->Fill(5.5);
0121       }
0122      }
0123       c->cd();
0124       hEventStat->Draw("hist-text");
0125       c->SaveAs("hist_origin_D0meson.png");
0126      fout->cd();
0127      hEventStat->Write();
0128       fout->Close();
0129 }
0130 // return 5 if any beauty ancestor is found, 4 if (at least) a charm ancestor but no beauty,
0131 int GetD0OriginFlavor(int imc,
0132                       const TTreeReaderArray<int> &mcPartPdg,
0133                       const TTreeReaderArray<unsigned int> &mcPartParent_begin,
0134                       const TTreeReaderArray<unsigned int> &mcPartParent_end,
0135                       const TTreeReaderArray<int> &mcPartParent_index)
0136 {
0137   int origin = 0;        // 0 = unknown, 4 = charm, 5 = beauty
0138   int current = imc;
0139   const int maxDepth = 20;  // safety
0140
0141   for (int depth = 0; depth < maxDepth; ++depth) {
0142
0143     int pdg = TMath::Abs(mcPartPdg[current]);
0144
0145     // charm quark or charm hadron
0146     if (pdg == 4 || (pdg/100 == 4) || (pdg/1000 == 4)) {
0147       if (origin == 0) origin = 4;
0148     }
0149
0150     // beauty quark or beauty hadron
0151     if (pdg == 5 || (pdg/100 == 5) || (pdg/1000 == 5)) {
0152       origin = 5;
0153       break;
0154     }
0155
0156     int firstParentIdx = mcPartParent_begin[current];
0157     int lastParentIdx  = mcPartParent_end[current];
0158
0159     if (firstParentIdx >= lastParentIdx) {
0160       break;
0161     }
0162
0163     int parentParticleIndex = mcPartParent_index[firstParentIdx];
0164     if (parentParticleIndex < 0) break;
0165
0166     current = parentParticleIndex;
0167   }
0168
0169   return origin;
0170 }