JetsAndHF/SecondaryVertex_Chi2/Check_decaylength.C

0001 // Code to check the decaylength of D0 at the simulation level.
0002 // Shyam Kumar; shyam.kumar@ba.infn.it; shyam055119@gmail.com
0003
0004 void Check_decaylength(TString filelist = "test.list"){
0005
0006     // Read all files
0007     int nfiles = 0;
0008      char filename[512];
0009     TChain *chain = new TChain("events");
0010
0011     ifstream *inputstream = new ifstream;
0012      inputstream->open(filelist.Data());
0013      if(!inputstream) printf("[Error] Cannot open file list: %s\n", filelist.Data());
0014
0015     while (inputstream->good())
0016       {
0017       inputstream->getline(filename, 512);
0018       if (inputstream->good())
0019       {
0020         TFile *ftmp = TFile::Open(filename, "READ");
0021         if (!ftmp || !ftmp->IsOpen() || !ftmp->GetNkeys())
0022         {
0023          printf("[e] Skipping bad file: %s\n", filename);
0024          if (ftmp) { ftmp->Close(); delete ftmp; }
0025          continue;
0026         }
0027         cout << "[i] Add " << nfiles << "th file: " << filename << endl;
0028         chain->Add(filename);
0029         nfiles++;
0030
0031         ftmp->Close();
0032         delete ftmp;
0033       }
0034     }
0035
0036     inputstream->close();
0037     printf("[i] Read in %d files with %lld events in total\n", nfiles, chain->GetEntries());
0038     TTreeReader myReader(chain);
0039
0040      // MC Track Properties
0041      TTreeReaderArray<int> mcPartGenStatus = {myReader, "MCParticles.generatorStatus"};
0042      TTreeReaderArray<int> mcPartPdg = {myReader, "MCParticles.PDG"};
0043      TTreeReaderArray<double> mcPartMass = {myReader, "MCParticles.mass"};
0044      TTreeReaderArray<double> mcEndPointX = {myReader, "MCParticles.endpoint.x"};
0045      TTreeReaderArray<double> mcEndPointY = {myReader, "MCParticles.endpoint.y"};
0046      TTreeReaderArray<double> mcEndPointZ = {myReader, "MCParticles.endpoint.z"};
0047      TTreeReaderArray<double> mcMomPx = {myReader, "MCParticles.momentum.x"};
0048      TTreeReaderArray<double> mcMomPy = {myReader, "MCParticles.momentum.y"};
0049      TTreeReaderArray<double> mcMomPz = {myReader, "MCParticles.momentum.z"};
0050      TTreeReaderArray<double> mcPartVx = {myReader, "MCParticles.vertex.x"};
0051      TTreeReaderArray<double> mcPartVy = {myReader, "MCParticles.vertex.y"};
0052      TTreeReaderArray<double> mcPartVz = {myReader, "MCParticles.vertex.z"};
0053      TTreeReaderArray<int> mcPartParent_index = {myReader, "_MCParticles_parents.index"};
0054      TTreeReaderArray<unsigned int> mcPartDaughter_begin = {myReader, "MCParticles.daughters_begin"};
0055      TTreeReaderArray<unsigned int> mcPartParent_begin = {myReader, "MCParticles.parents_begin"};
0056      TTreeReaderArray<unsigned int> mcPartParent_end = {myReader, "MCParticles.parents_end"};
0057      TTreeReaderArray<unsigned int> mcPartDaughter_end = {myReader, "MCParticles.daughters_end"};
0058      TTreeReaderArray<int> mcPartDaughter_index = {myReader, "_MCParticles_daughters.index"};
0059
0060   TFile *fout = new TFile("output_file.root","recreate");
0061   TH2D *hptD0vsdl = new TH2D("hptD0vsdl","hptD0vsdl;pD0 (GeV/c); dl (mm)",100.,0.,20.,100,0.,5.0);
0062   int count = -1;
0063   TDatabasePDG *pdgDB = TDatabasePDG::Instance();
0064   while (myReader.Next())
0065   {
0066    count++;
0067  //  if (count!=nevent) continue;
0068    printf("Event No. = %d\n",count);
0069     // find MC primary vertex
0070     int nMCPart = mcPartMass.GetSize();
0071
0072     TVector3 vertex_mc(-999., -999., -999.);
0073     TVector3 SV_mc(-999., -999., -999.);
0074     TVector3 pD0(0.,0.,0.);
0075    for(int imc=0; imc<nMCPart; imc++)
0076       {
0077       if(mcPartGenStatus[imc] == 4 && mcPartPdg[imc] == 11)
0078         {
0079           vertex_mc.SetXYZ(mcEndPointX[imc], mcEndPointY[imc], mcEndPointZ[imc]);
0080           break;
0081         }
0082        }
0083
0084     for(int imc=0; imc<nMCPart; imc++)
0085     {
0086       if(fabs(mcPartPdg[imc]) != 421) continue;
0087       pD0.SetXYZ(mcMomPx[imc],mcMomPy[imc],mcMomPz[imc]);
0088       int nDaughters = mcPartDaughter_end[imc]-mcPartDaughter_begin[imc];
0089       if(nDaughters!=2) continue;
0090       // A "status 21" particle might represent the before-hadronization D⁰.
0091        //A "status 1" particle is the final, stable D⁰ after decays or transport.
0092       int parent_index_begin = mcPartParent_index[mcPartParent_begin[imc]];
0093       int parent_index_end = mcPartParent_index[mcPartParent_end[imc]];
0094       TString particle1 = pdgDB->GetParticle(mcPartPdg[parent_index_begin])->GetName();
0095       TString particle2 = pdgDB->GetParticle(mcPartPdg[parent_index_end])->GetName();
0096     printf("Parent Name  = (%s, %s) \n",particle1.Data(),particle2.Data());
0097       // find D0 that decay into pi+K
0098       bool is_pik_decay = false;
0099       int daug_index_1 = mcPartDaughter_index[mcPartDaughter_begin[imc]];
0100       int daug_index_2 = mcPartDaughter_index[mcPartDaughter_begin[imc]+1];
0101       int daug_pdg_1 = mcPartPdg[daug_index_1];
0102       int daug_pdg_2 = mcPartPdg[daug_index_2];
0103       if( (fabs(daug_pdg_1)==321 && fabs(daug_pdg_2)==211) || (fabs(daug_pdg_1)==211 && fabs(daug_pdg_2)==321) )
0104       {
0105           is_pik_decay = true;
0106           SV_mc.SetXYZ(mcPartVx[daug_index_1],mcPartVy[daug_index_1],mcPartVz[daug_index_1]);
0107           SV_mc.Print();
0108           hptD0vsdl->Fill(pD0.Mag(),(SV_mc-vertex_mc).Mag());
0109       }
0110
0111     }
0112
0113      }
0114      fout->cd();
0115    hptD0vsdl->Write();
0116    fout->Close();
0117
0118 }