ITS: extend macros for pull distributions

Signed-off-by: Felix Schlepper <[email protected]>

Author: f3sch

Detectors/ITSMFT/ITS/macros/test/CheckTracksCA.C

@@ -31,6 +31,7 @@
 #include "SimulationDataFormat/MCEventHeader.h"
 #include "DetectorsBase/Propagator.h"
 #include "SimulationDataFormat/TrackReference.h"
+#include "SimulationDataFormat/O2DatabasePDG.h"
 #include "SimulationDataFormat/MCTrack.h"
 #include "SimulationDataFormat/MCCompLabel.h"
 #include "SimulationDataFormat/MCTruthContainer.h"
@@ -60,11 +61,15 @@ struct ParticleInfo {
   bool isPrimary = 0u;
   unsigned char storedStatus = 2; /// not stored = 2, fake = 1, good = 0
   o2::its::TrackITS track;
+  o2::MCTrack mcTrack;
 };
 
 #pragma link C++ class ParticleInfo + ;
 
-void CheckTracksCA(bool doFakeClStud = false,
+void CheckTracksCA(bool doEffStud = false,
+                   bool doFakeClStud = false,
+                   bool doPullStud = true,
+                   bool createOutput = false,
                    std::string tracfile = "o2trac_its.root",
                    std::string magfile = "o2sim_grp.root",
                    std::string clusfile = "o2clus_its.root",
@@ -124,7 +129,10 @@ void CheckTracksCA(bool doFakeClStud = false,
     info[n].resize(mcArr->size());
     hZvertex->Fill(mcEvent->GetZ());
     for (unsigned int mcI{0}; mcI < mcArr->size(); ++mcI) {
-      auto part = mcArr->at(mcI);
+      const auto part = mcArr->at(mcI);
+      if (!o2::O2DatabasePDG::Instance()->GetParticle(part.GetPdgCode())) {
+        continue;
+      }
       info[n][mcI].event = n;
       info[n][mcI].pdg = part.GetPdgCode();
       info[n][mcI].pvx = mcEvent->GetX();
@@ -134,6 +142,7 @@ void CheckTracksCA(bool doFakeClStud = false,
       info[n][mcI].phi = part.GetPhi();
       info[n][mcI].eta = part.GetEta();
       info[n][mcI].isPrimary = part.isPrimary();
+      info[n][mcI].mcTrack = part;
     }
   }
   std::cout << "done." << std::endl;
@@ -214,117 +223,108 @@ void CheckTracksCA(bool doFakeClStud = false,
   std::cout << "\t- Total number of fakes: " << fakes << " (" << fakes * 100. / total << "%)" << std::endl;
   std::cout << "\t- Total number of good: " << good << " (" << good * 100. / total << "%)" << std::endl;
 
-  int nb = 100;
-  double xbins[nb + 1], ptcutl = 0.01, ptcuth = 10.;
-  double a = std::log(ptcuth / ptcutl) / nb;
-  for (int i = 0; i <= nb; i++)
-    xbins[i] = ptcutl * std::exp(i * a);
-  TH1D* num = new TH1D("num", ";#it{p}_{T} (GeV/#it{c});Efficiency (fake-track rate)", nb, xbins);
-  num->Sumw2();
-  TH1D* numEta = new TH1D("numEta", ";#eta;Number of tracks", 60, -3, 3);
-  numEta->Sumw2();
-  TH1D* numChi2 = new TH1D("numChi2", ";#it{p}_{T} (GeV/#it{c});Efficiency (fake-track rate)", 200, 0, 100);
-
-  TH1D* fak = new TH1D("fak", ";#it{p}_{T} (GeV/#it{c});Fak", nb, xbins);
-  fak->Sumw2();
-  TH1D* multiFak = new TH1D("multiFak", ";#it{p}_{T} (GeV/#it{c});Fak", nb, xbins);
-  multiFak->Sumw2();
-  TH1D* fakChi2 = new TH1D("fakChi2", ";#it{p}_{T} (GeV/#it{c});Fak", 200, 0, 100);
-
-  TH1D* clone = new TH1D("clone", ";#it{p}_{T} (GeV/#it{c});Clone", nb, xbins);
-  clone->Sumw2();
-
-  TH1D* den = new TH1D("den", ";#it{p}_{T} (GeV/#it{c});Den", nb, xbins);
-  den->Sumw2();
-
-  for (auto& evInfo : info) {
-    for (auto& part : evInfo) {
-      if ((part.clusters & 0x7f) != 0x7f) {
-        // part.clusters != 0x3f && part.clusters != 0x3f << 1 &&
-        // part.clusters != 0x1f && part.clusters != 0x1f << 1 && part.clusters != 0x1f << 2 &&
-        // part.clusters != 0x0f && part.clusters != 0x0f << 1 && part.clusters != 0x0f << 2 && part.clusters != 0x0f << 3) {
-        continue;
-      }
-      if (!part.isPrimary) {
-        continue;
-      }
-      den->Fill(part.pt);
-      if (part.isReco) {
-        num->Fill(part.pt);
-        numEta->Fill(part.eta);
-        if (part.isReco > 1) {
-          for (int _i{0}; _i < part.isReco - 1; ++_i) {
-            clone->Fill(part.pt);
+  TFile* file{nullptr};
+  if (createOutput) {
+    file = TFile::Open("CheckTracksCA.root", "recreate");
+  }
+
+  if (doEffStud) {
+    std::cout << "Calculating efficiencies... ";
+    const int nb = 100;
+    double xbins[nb + 1], ptcutl = 0.01, ptcuth = 10.;
+    double a = std::log(ptcuth / ptcutl) / nb;
+    for (int i = 0; i <= nb; i++)
+      xbins[i] = ptcutl * std::exp(i * a);
+    TH1D* num = new TH1D("num", ";#it{p}_{T} (GeV/#it{c});Efficiency (fake-track rate)", nb, xbins);
+    num->Sumw2();
+    TH1D* numEta = new TH1D("numEta", ";#eta;Number of tracks", 60, -3, 3);
+    numEta->Sumw2();
+    TH1D* numChi2 = new TH1D("numChi2", ";#it{p}_{T} (GeV/#it{c});Efficiency (fake-track rate)", 200, 0, 100);
+
+    TH1D* fak = new TH1D("fak", ";#it{p}_{T} (GeV/#it{c});Fak", nb, xbins);
+    fak->Sumw2();
+    TH1D* multiFak = new TH1D("multiFak", ";#it{p}_{T} (GeV/#it{c});Fak", nb, xbins);
+    multiFak->Sumw2();
+    TH1D* fakChi2 = new TH1D("fakChi2", ";#it{p}_{T} (GeV/#it{c});Fak", 200, 0, 100);
+
+    TH1D* clone = new TH1D("clone", ";#it{p}_{T} (GeV/#it{c});Clone", nb, xbins);
+    clone->Sumw2();
+
+    TH1D* den = new TH1D("den", ";#it{p}_{T} (GeV/#it{c});Den", nb, xbins);
+    den->Sumw2();
+
+    for (auto& evInfo : info) {
+      for (auto& part : evInfo) {
+        if ((part.clusters & 0x7f) != 0x7f) {
+          // part.clusters != 0x3f && part.clusters != 0x3f << 1 &&
+          // part.clusters != 0x1f && part.clusters != 0x1f << 1 && part.clusters != 0x1f << 2 &&
+          // part.clusters != 0x0f && part.clusters != 0x0f << 1 && part.clusters != 0x0f << 2 && part.clusters != 0x0f << 3) {
+          continue;
+        }
+        if (!part.isPrimary) {
+          continue;
+        }
+        den->Fill(part.pt);
+        if (part.isReco) {
+          num->Fill(part.pt);
+          numEta->Fill(part.eta);
+          if (part.isReco > 1) {
+            for (int _i{0}; _i < part.isReco - 1; ++_i) {
+              clone->Fill(part.pt);
+            }
           }
         }
-      }
-      if (part.isFake) {
-        fak->Fill(part.pt);
-        if (part.isFake > 1) {
-          for (int _i{0}; _i < part.isFake - 1; ++_i) {
-            multiFak->Fill(part.pt);
+        if (part.isFake) {
+          fak->Fill(part.pt);
+          if (part.isFake > 1) {
+            for (int _i{0}; _i < part.isFake - 1; ++_i) {
+              multiFak->Fill(part.pt);
+            }
           }
         }
       }
     }
-  }
 
-  TCanvas* c1 = new TCanvas;
-  c1->SetLogx();
-  c1->SetGridx();
-  c1->SetGridy();
-  TH1* sum = (TH1*)num->Clone("sum");
-  sum->Add(fak);
-  sum->Divide(sum, den, 1, 1);
-  sum->SetLineColor(kBlack);
-  sum->Draw("hist");
-  num->Divide(num, den, 1, 1, "b");
-  num->Draw("histesame");
-  fak->Divide(fak, den, 1, 1, "b");
-  fak->SetLineColor(2);
-  fak->Draw("histesame");
-  multiFak->Divide(multiFak, den, 1, 1, "b");
-  multiFak->SetLineColor(kRed + 1);
-  multiFak->Draw("histsame");
-  clone->Divide(clone, den, 1, 1, "b");
-  clone->SetLineColor(3);
-  clone->Draw("histesame");
-  TCanvas* c2 = new TCanvas;
-  c2->SetGridx();
-  c2->SetGridy();
-  hZvertex->DrawClone();
-
-  std::cout << "** Streaming output TTree to file ... " << std::flush;
-  TFile file("CheckTracksCA.root", "recreate");
-  TTree tree("ParticleInfo", "ParticleInfo");
-  ParticleInfo pInfo;
-  tree.Branch("particle", &pInfo);
-  for (auto& event : info) {
-    for (auto& part : event) {
-      int nCl{0};
-      for (unsigned int bit{0}; bit < sizeof(pInfo.clusters) * 8; ++bit) {
-        nCl += bool(part.clusters & (1 << bit));
-      }
-      if (nCl < 3) {
-        continue;
-      }
-      pInfo = part;
-      tree.Fill();
+    TCanvas* c1 = new TCanvas;
+    c1->SetLogx();
+    c1->SetGridx();
+    c1->SetGridy();
+    TH1* sum = (TH1*)num->Clone("sum");
+    sum->Add(fak);
+    sum->Divide(sum, den, 1, 1);
+    sum->SetLineColor(kBlack);
+    sum->Draw("hist");
+    num->Divide(num, den, 1, 1, "b");
+    num->Draw("histesame");
+    fak->Divide(fak, den, 1, 1, "b");
+    fak->SetLineColor(2);
+    fak->Draw("histesame");
+    multiFak->Divide(multiFak, den, 1, 1, "b");
+    multiFak->SetLineColor(kRed + 1);
+    multiFak->Draw("histsame");
+    clone->Divide(clone, den, 1, 1, "b");
+    clone->SetLineColor(3);
+    clone->Draw("histesame");
+    TCanvas* c2 = new TCanvas;
+    c2->SetGridx();
+    c2->SetGridy();
+    hZvertex->DrawClone();
+
+    if (createOutput) {
+      sum->Write("total");
+      fak->Write("singleFake");
+      num->Write("efficiency");
+      numEta->Write("etaDist");
+      multiFak->Write("multiFake");
+      clone->Write("clones");
     }
+    std::cout << " done." << std::endl;
   }
-  tree.Write();
-  sum->Write("total");
-  fak->Write("singleFake");
-  num->Write("efficiency");
-  numEta->Write("etaDist");
-  multiFak->Write("multiFake");
-  clone->Write("clones");
-  file.Close();
-  std::cout << " done." << std::endl;
 
   //////////////////////
   // Fake clusters study
   if (doFakeClStud) {
+    std::cout << "Creating fake cluster study... ";
     std::vector<TH1I*> histLength, histLength1Fake, histLengthNoCl, histLength1FakeNoCl;
     std::vector<THStack*> stackLength, stackLength1Fake;
     std::vector<TLegend*> legends, legends1Fake;
@@ -364,10 +364,10 @@ void CheckTracksCA(bool doFakeClStud = false,
       stackLength1Fake[iH - 4]->Add(histLength1FakeNoCl[iH - 4]);
     }
 
-    for (auto& event : info) {
-      for (auto& part : event) {
+    for (const auto& event : info) {
+      for (const auto& part : event) {
         int nCl{0};
-        for (unsigned int bit{0}; bit < sizeof(pInfo.clusters) * 8; ++bit) {
+        for (unsigned int bit{0}; bit < sizeof(part.clusters) * 8; ++bit) {
           nCl += bool(part.clusters & (1 << bit));
         }
         if (nCl < 3) {
@@ -409,5 +409,93 @@ void CheckTracksCA(bool doFakeClStud = false,
       gPad->BuildLegend();
     }
     canvas->SaveAs("fakeClusters.png", "recreate");
+    std::cout << " done\n";
+  }
+
+  if (doPullStud) {
+    std::cout << "Creating pull study... ";
+    const int nBins{30};
+    const float xWidth{10};
+    // Pulls
+    auto hYPull = new TH1F("hYPull", "Pull Y", nBins, -xWidth, xWidth);
+    auto hZPull = new TH1F("hZPull", "Pull Z", nBins, -xWidth, xWidth);
+    auto hSPhiPull = new TH1F("hSPhiPull", "Pull Sin(Phi)", nBins, -xWidth, xWidth);
+    auto hTglPull = new TH1F("hTglPull", "Pull Tg(Lambda)", nBins, -xWidth, xWidth);
+    auto hQoPtPull = new TH1F("hQoPtPull", "Pull Q/Pt", nBins, -xWidth, xWidth);
+    // Correlation
+    auto hCorYZ = new TH2F("hCorYZ", ";#sigma_{Z}^{2};#sigma_{Y}^{2}", nBins, 0, 5e-5, nBins, 0, 5e-5);
+    auto hCorYSPhi = new TH2F("hCorYSPhi", ";#sigma_{snp}^{2};#sigma_{Y}^{2}", nBins, 0, 5e-5, nBins, 0, 5e-5);
+
+    for (const auto& event : info) {
+      for (const auto& part : event) {
+        if (((part.clusters & 0x7f) != 0x7f) && !part.isPrimary) {
+          continue;
+        }
+
+        // prepare mc truth parameters
+        std::array<float, 3> xyz{(float)part.mcTrack.GetStartVertexCoordinatesX(), (float)part.mcTrack.GetStartVertexCoordinatesY(), (float)part.mcTrack.GetStartVertexCoordinatesZ()};
+        std::array<float, 3> pxyz{(float)part.mcTrack.GetStartVertexMomentumX(), (float)part.mcTrack.GetStartVertexMomentumY(), (float)part.mcTrack.GetStartVertexMomentumZ()};
+        o2::track::TrackPar mcTrack(xyz, pxyz, TMath::Nint(o2::O2DatabasePDG::Instance()->GetParticle(part.mcTrack.GetPdgCode())->Charge() / 3), false);
+        if (!mcTrack.rotate(part.track.getAlpha()) &&
+            o2::base::Propagator::Instance()->propagateTo(mcTrack, part.track.getX())) {
+          continue;
+        }
+
+        hYPull->Fill((part.track.getY() - mcTrack.getY()) / std::sqrt(part.track.getSigmaY2()));
+        hZPull->Fill((part.track.getZ() - mcTrack.getZ()) / std::sqrt(part.track.getSigmaZ2()));
+        hSPhiPull->Fill((part.track.getSnp() - mcTrack.getSnp()) / std::sqrt(part.track.getSigmaSnp2()));
+        hTglPull->Fill((part.track.getTgl() - mcTrack.getTgl()) / std::sqrt(part.track.getSigmaTgl2()));
+        hQoPtPull->Fill((part.track.getQ2Pt() - mcTrack.getQ2Pt()) / std::sqrt(part.track.getSigma1Pt2()));
+
+        hCorYZ->Fill(part.track.getSigmaZ2(), part.track.getSigmaY2());
+        hCorYSPhi->Fill(part.track.getSigmaZ2(), part.track.getSigmaSnp2());
+      }
+    }
+
+    // normalise, set axis, fit and draw
+    auto doPullCalc = [](TH1F* h) {
+      h->Scale(1. / h->Integral("width"));
+      h->GetYaxis()->SetRangeUser(1e-5, 1.);
+      h->Fit("gaus", "", "2QMR");
+      gPad->SetLogy();
+      h->Draw("hist");
+    };
+
+    auto c = new TCanvas("cPull", "", 1000, 1000);
+    c->Divide(5, 5);
+    c->cd(1);
+    doPullCalc(hYPull);
+    c->cd(2);
+    hCorYZ->Draw("colz");
+    c->cd(3);
+    hCorYSPhi->Draw("colz");
+    c->cd(1 + 6);
+    doPullCalc(hZPull);
+    c->cd(1 + 12);
+    doPullCalc(hSPhiPull);
+    c->cd(1 + 18);
+    doPullCalc(hTglPull);
+    c->cd(1 + 24);
+    doPullCalc(hQoPtPull);
+    c->Draw();
+    std::cout << " done\n";
+  }
+
+  if (createOutput) {
+    std::cout << "** Streaming output TTree to file ... " << std::flush;
+    TTree tree("ParticleInfo", "ParticleInfo");
+    ParticleInfo pInfo;
+    tree.Branch("particle", &pInfo);
+    for (const auto& event : info) {
+      for (const auto& part : event) {
+        if (((part.clusters & 0x7f) != 0x7f) && !part.isPrimary) {
+          continue;
+        }
+        pInfo = part;
+        tree.Fill();
+      }
+    }
+    tree.Write();
+    file->Close();
   }
 }