|
| 1 | +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. |
| 2 | +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. |
| 3 | +// All rights not expressly granted are reserved. |
| 4 | +// |
| 5 | +// This software is distributed under the terms of the GNU General Public |
| 6 | +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". |
| 7 | +// |
| 8 | +// In applying this license CERN does not waive the privileges and immunities |
| 9 | +// granted to it by virtue of its status as an Intergovernmental Organization |
| 10 | +// or submit itself to any jurisdiction. |
| 11 | + |
| 12 | +#include "FastTracker.h" |
| 13 | +#include "TrackUtilities.h" |
| 14 | + |
| 15 | +#include <DataFormatsParameters/GRPMagField.h> |
| 16 | +#include <DetectorsBase/Propagator.h> |
| 17 | +#include <ReconstructionDataFormats/Vertex.h> |
| 18 | + |
| 19 | +#include <TAxis.h> |
| 20 | +#include <TCanvas.h> |
| 21 | +#include <TDatabasePDG.h> |
| 22 | +#include <TEfficiency.h> |
| 23 | +#include <TGraph.h> |
| 24 | +#include <TH1F.h> |
| 25 | +#include <TLatex.h> |
| 26 | +#include <TLorentzVector.h> |
| 27 | + |
| 28 | +void drawFastTracker(float magneticField = 5.f, // in units of kGauss |
| 29 | + const int nch = 100, // number of charged particles per unit rapidity |
| 30 | + const int pdg = 211) // PDG code of the particle to track |
| 31 | +{ |
| 32 | + TDatabasePDG* db = TDatabasePDG::Instance(); |
| 33 | + TParticlePDG* p = 0; |
| 34 | + p = db->GetParticle(pdg); |
| 35 | + if (!p) { |
| 36 | + LOG(fatal) << "Particle with PDG code " << pdg << " not found in TDatabasePDG"; |
| 37 | + return; |
| 38 | + } |
| 39 | + const float mass = p->Mass(); // particle mass in GeV/c^2 |
| 40 | + const float q = p->Charge() / 3.0; // charge in e |
| 41 | + |
| 42 | + o2::parameters::GRPMagField grpmag; |
| 43 | + grpmag.setFieldUniformity(true); |
| 44 | + grpmag.setL3Current(30000.f * (magneticField / 5.0f)); |
| 45 | + auto field = grpmag.getNominalL3Field(); |
| 46 | + o2::base::Propagator::initFieldFromGRP(&grpmag); |
| 47 | + |
| 48 | + fair::Logger::SetVerbosity(fair::Verbosity::verylow); |
| 49 | + o2::fastsim::FastTracker fastTracker = o2::fastsim::FastTracker(); |
| 50 | + if (0) { |
| 51 | + fastTracker.SetApplyEffCorrection(false); |
| 52 | + Double_t x0IB = 0.001; |
| 53 | + Double_t x0OB = 0.01; |
| 54 | + Double_t xrhoIB = 2.3292e-02; // 100 mum Si |
| 55 | + Double_t xrhoOB = 2.3292e-01; // 1000 mum Si |
| 56 | + Double_t resRPhiIB = 0.00025; |
| 57 | + Double_t resZIB = 0.00025; |
| 58 | + Double_t resRPhiOB = 0.00100; |
| 59 | + Double_t resZOB = 0.00100; |
| 60 | + Double_t eff = 0.98; |
| 61 | + // fastTracker.AddLayer("vertex", 0.0, 250, 0, 0); // dummy vertex for matrix calculation |
| 62 | + fastTracker.AddLayer("bpipe0", 0.48, 250, 0.00042, 2.772e-02); // 150 mum Be |
| 63 | + fastTracker.AddLayer("B00", 0.50, 250, x0IB, xrhoIB, resRPhiIB, resZIB, eff, 1); |
| 64 | + fastTracker.AddLayer("B01", 1.20, 250, x0IB, xrhoIB, resRPhiIB, resZIB, eff, 1); |
| 65 | + fastTracker.AddLayer("B02", 2.50, 250, x0IB, xrhoIB, resRPhiIB, resZIB, eff, 1); |
| 66 | + fastTracker.AddLayer("bpipe1", 3.7, 250, 0.0014, 9.24e-02); // 500 mum Be |
| 67 | + fastTracker.AddLayer("B03", 3.75, 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1); |
| 68 | + fastTracker.AddLayer("B04", 7.00, 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1); |
| 69 | + fastTracker.AddLayer("B05", 12.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1); |
| 70 | + fastTracker.AddLayer("B06", 20.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1); |
| 71 | + fastTracker.AddLayer("B07", 30.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1); |
| 72 | + fastTracker.AddLayer("B08", 45.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1); |
| 73 | + fastTracker.AddLayer("B09", 60.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1); |
| 74 | + fastTracker.AddLayer("B10", 80.0, 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1); |
| 75 | + fastTracker.AddLayer("B11", 100., 250, x0OB, xrhoOB, resRPhiOB, resZOB, eff, 1); |
| 76 | + } else { |
| 77 | + std::vector<float> pixelRes{0.025, 0.025, 0.01, 0.01}; |
| 78 | + fastTracker.AddSiliconALICE3v4(pixelRes); |
| 79 | + } |
| 80 | + |
| 81 | + fastTracker.Print(); |
| 82 | + fastTracker.SetMagneticField(magneticField); |
| 83 | + |
| 84 | + TAxis ptBinning(1000, 0., 10); |
| 85 | + TGraph* gPt = new TGraph(); |
| 86 | + gPt->GetXaxis()->SetTitle("#it{p}_{T} (GeV/c)"); |
| 87 | + gPt->GetYaxis()->SetTitle("Efficiency"); |
| 88 | + TEfficiency* hEfficiency = new TEfficiency("hEfficiency", ";#it{p}_{T} (GeV/c);Efficiency", ptBinning.GetNbins(), ptBinning.GetBinLowEdge(1), ptBinning.GetBinUpEdge(ptBinning.GetNbins())); |
| 89 | + TH1F* hFastTrackerQA = new TH1F("hFastTrackerQA", ";#it{p}_{T} (GeV/c;Tracking code", ptBinning.GetNbins(), ptBinning.GetBinLowEdge(1), ptBinning.GetBinUpEdge(ptBinning.GetNbins())); |
| 90 | + |
| 91 | + TLorentzVector tlv; |
| 92 | + |
| 93 | + // Now we compute the efficiency for a range of pt values and plot the results |
| 94 | + o2::track::TrackParCov trkIn; |
| 95 | + o2::track::TrackParCov trkOut; |
| 96 | + for (int i = 1; i <= ptBinning.GetNbins(); ++i) { |
| 97 | + const float pt = ptBinning.GetBinCenter(i); |
| 98 | + tlv.SetPtEtaPhiM(pt, 0., 0., mass); |
| 99 | + o2::upgrade::convertTLorentzVectorToO2Track(q, tlv, {0.0, 0., 0.}, trkIn); |
| 100 | + |
| 101 | + bool prop = o2::base::Propagator::Instance()->propagateToR(trkIn, 0.4f, true, o2::base::Propagator::MAX_SIN_PHI, o2::base::Propagator::MAX_STEP, o2::base::Propagator::MatCorrType::USEMatCorrNONE); |
| 102 | + if (!prop) { |
| 103 | + LOG(info) << "Initial propagation to R=0 failed for pt = " << pt; |
| 104 | + continue; |
| 105 | + } |
| 106 | + |
| 107 | + o2::dataformats::VertexBase vertex{{0.f, 0.f, 0.f}, {0.f, 0.f, 0.f, 0.f, 0.f, 0.f}}; |
| 108 | + // Now we propagate to the DCA of the vertex |
| 109 | + if (!trkIn.propagateToDCA(vertex, magneticField)) { |
| 110 | + LOG(info) << "Propagation to DCA failed for pt = " << pt; |
| 111 | + continue; |
| 112 | + } |
| 113 | + |
| 114 | + for (int trial = 0; trial < 200; trial++) { |
| 115 | + hEfficiency->Fill(fastTracker.FastTrack(trkIn, trkOut, nch) > 0, pt); |
| 116 | + } |
| 117 | + int status = 4; |
| 118 | + status = fastTracker.FastTrack(trkIn, trkOut, nch); |
| 119 | + hFastTrackerQA->Fill(pt, status); |
| 120 | + if (status < 0) { |
| 121 | + LOG(debug) << " --- fatSolve: FastTrack failed with status " << status << " --- "; |
| 122 | + // tlv.Print(); |
| 123 | + gPt->AddPoint(pt, 0); |
| 124 | + continue; |
| 125 | + } |
| 126 | + // define the efficiency |
| 127 | + float eff = 1.; |
| 128 | + for (size_t l = 1; l < fastTracker.GetNLayers(); ++l) { |
| 129 | + if (fastTracker.IsLayerInert(l)) { |
| 130 | + continue; // skip inert layers |
| 131 | + } |
| 132 | + float igoodhit = 0.f; |
| 133 | + igoodhit = fastTracker.GetGoodHitProb(l); |
| 134 | + if (igoodhit <= 0.) { |
| 135 | + continue; |
| 136 | + } |
| 137 | + eff *= igoodhit; |
| 138 | + } |
| 139 | + gPt->AddPoint(pt, eff); |
| 140 | + } |
| 141 | + |
| 142 | + TLatex latex; |
| 143 | + latex.SetTextFont(42); |
| 144 | + latex.SetNDC(); |
| 145 | + latex.SetTextSize(0.037); |
| 146 | + |
| 147 | + new TCanvas("gPt"); |
| 148 | + gPt->Draw("ALP"); |
| 149 | + hEfficiency->Draw("SAME"); |
| 150 | + latex.DrawLatex(0.1, 0.93, Form("#LT #frac{dN_{ch}}{d#eta} #GT = %i", nch)); |
| 151 | + |
| 152 | + new TCanvas("hEfficiency"); |
| 153 | + hEfficiency->Draw(""); |
| 154 | + latex.DrawLatex(0.1, 0.93, Form("#LT #frac{dN_{ch}}{d#eta} #GT = %i", nch)); |
| 155 | + new TCanvas("hFastTrackerQA"); |
| 156 | + hFastTrackerQA->Print(); |
| 157 | + hFastTrackerQA->Draw(""); |
| 158 | + latex.DrawLatex(0.1, 0.93, Form("#LT #frac{dN_{ch}}{d#eta} #GT = %i", nch)); |
| 159 | +} |
0 commit comments