JetsAndHF/Lcp_Reconstruction/Check_daughters.C

0001 // Macro to check daughters of Lc+ baryon
0002 // Shyam Kumar; INFN Bari, shyam.kumar@ba.infn.it
0003 void Check_daughters(TString filelist = "test_Lc.list")
0004 {
0005
0006
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 treereader(chain);
0040
0041   // MC
0042   TTreeReaderArray<int> mcPartGenStatus = {treereader, "MCParticles.generatorStatus"};
0043   TTreeReaderArray<int> mcPartPdg = {treereader, "MCParticles.PDG"};
0044   TTreeReaderArray<float> mcPartCharge = {treereader, "MCParticles.charge"};
0045   TTreeReaderArray<unsigned int> mcPartParent_begin = {treereader, "MCParticles.parents_begin"};
0046   TTreeReaderArray<unsigned int> mcPartParent_end = {treereader, "MCParticles.parents_end"};
0047   TTreeReaderArray<int> mcPartParent_index = {treereader, "_MCParticles_parents.index"};
0048   TTreeReaderArray<unsigned int> mcPartDaughter_begin = {treereader, "MCParticles.daughters_begin"};
0049   TTreeReaderArray<unsigned int> mcPartDaughter_end = {treereader, "MCParticles.daughters_end"};
0050   TTreeReaderArray<int> mcPartDaughter_index = {treereader, "_MCParticles_daughters.index"};
0051   TTreeReaderArray<double> mcPartMass = {treereader, "MCParticles.mass"};
0052   TTreeReaderArray<double> mcPartVx = {treereader, "MCParticles.vertex.x"};
0053   TTreeReaderArray<double> mcPartVy = {treereader, "MCParticles.vertex.y"};
0054   TTreeReaderArray<double> mcPartVz = {treereader, "MCParticles.vertex.z"};
0055   TTreeReaderArray<double> mcMomPx = {treereader, "MCParticles.momentum.x"};
0056   TTreeReaderArray<double> mcMomPy = {treereader, "MCParticles.momentum.y"};
0057   TTreeReaderArray<double> mcMomPz = {treereader, "MCParticles.momentum.z"};
0058   TTreeReaderArray<double> mcEndPointX = {treereader, "MCParticles.endpoint.x"};
0059   TTreeReaderArray<double> mcEndPointY = {treereader, "MCParticles.endpoint.y"};
0060   TTreeReaderArray<double> mcEndPointZ = {treereader, "MCParticles.endpoint.z"};
0061
0062     TDatabasePDG *pdgDB = TDatabasePDG::Instance();
0063     Int_t iEvent = 0;
0064     while(treereader.Next())
0065     {
0066    //   printf("Event No. = %d \n", iEvent);
0067         int nMCPart = mcPartMass.GetSize();
0068         for(int imc=0; imc<nMCPart; imc++)
0069     {
0070       if(fabs(mcPartPdg[imc]) != 4122) continue;
0071
0072       TString particle_name = pdgDB->GetParticle(mcPartPdg[imc])->GetName();
0073
0074       int nDaughters = mcPartDaughter_end[imc]-mcPartDaughter_begin[imc];
0075       if(nDaughters!=3) continue;
0076       int daug_index_1 = mcPartDaughter_index[mcPartDaughter_begin[imc]];
0077       int daug_index_2 = mcPartDaughter_index[mcPartDaughter_begin[imc]+1];
0078       int daug_index_3 = mcPartDaughter_index[mcPartDaughter_begin[imc]+2];
0079      // int daug_index_4 = mcPartDaughter_index[mcPartDaughter_begin[imc]+3];
0080       int daug_pdg_1 = mcPartPdg[daug_index_1];
0081       int daug_pdg_2 = mcPartPdg[daug_index_2];
0082       int daug_pdg_3 = mcPartPdg[daug_index_3];
0083      // int daug_pdg_4 = mcPartPdg[daug_index_4];
0084       if (fabs(daug_pdg_1)!=13 && fabs(daug_pdg_2)!=13 && fabs(daug_pdg_3)!=13) continue;
0085       printf("Event No. = %d \n", iEvent);
0086       printf("Parent Particle = %s \n", particle_name.Data());
0087       printf("PDG name = %d \t %d \t %d \t \n",daug_pdg_1, daug_pdg_2, daug_pdg_3);
0088       printf("Daughters name = %s \t %s \t %s \t \n", pdgDB->GetParticle(daug_pdg_1)->GetName(),pdgDB->GetParticle(daug_pdg_2)->GetName(), pdgDB->GetParticle(daug_pdg_3)->GetName());
0089
0090
0091
0092
0093     }
0094
0095     iEvent++;
0096
0097     }
0098     delete inputstream;
0099     delete chain;
0100 }