Tracking/ImpactPointEstimator/pca_global_impactpoint.C

0001 #include "DD4hep/Detector.h"
0002 #include "DD4hep/DetElement.h"
0003
0004 #include "Acts/Plugins/DD4hep/ConvertDD4hepDetector.hpp"
0005 #include "Acts/Plugins/DD4hep/DD4hepFieldAdapter.hpp"
0006 #include "Acts/Geometry/TrackingGeometry.hpp"
0007 #include "Acts/Geometry/GeometryContext.hpp"
0008 #include "Acts/Plugins/DD4hep/ConvertDD4hepDetector.hpp"
0009 #include "Acts/Utilities/Logger.hpp"
0010 #include "Acts/Utilities/BinningType.hpp"
0011 #include "Acts/Definitions/Units.hpp"
0012 #include "Acts/Surfaces/PerigeeSurface.hpp"
0013 #include "Acts/EventData/TrackParameters.hpp"
0014 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
0015 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0016 #include "Acts/Propagator/EigenStepper.hpp"
0017 #include "Acts/Propagator/Propagator.hpp"
0018 #include "Acts/Vertexing/ImpactPointEstimator.hpp"
0019
0020 void pca_global_impactpoint(){
0021
0022     // Load DD4Hep geometry

0023     dd4hep::Detector& detector = dd4hep::Detector::getInstance();
0024     detector.fromCompact("/opt/detector/epic-main/share/epic/epic_craterlake.xml");
0025     dd4hep::DetElement geometry = detector.world();
0026
0027     // Convert DD4Hep geometry to tracking geometry

0028     Acts::GeometryContext trackingGeoCtx;
0029     auto logger = Acts::getDefaultLogger("DD4hepConversion", Acts::Logging::Level::INFO);
0030     Acts::BinningType bTypePhi = Acts::equidistant;
0031     Acts::BinningType bTypeR = Acts::equidistant;
0032     Acts::BinningType bTypeZ = Acts::equidistant;
0033     double layerEnvelopeR = Acts::UnitConstants::mm;
0034     double layerEnvelopeZ = Acts::UnitConstants::mm;
0035     double defaultLayerThickness = Acts::UnitConstants::fm;
0036     using Acts::sortDetElementsByID;
0037
0038     std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry{nullptr};
0039     trackingGeometry = Acts::convertDD4hepDetector(geometry,*logger,bTypePhi,bTypeR,bTypeZ,layerEnvelopeR,layerEnvelopeZ,defaultLayerThickness,sortDetElementsByID,trackingGeoCtx);
0040
0041     // Define Perigee surface at which reconstructed track parameters are set

0042     auto perigee = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3(0,0,0));
0043
0044     // Get Magnetic field context

0045     Acts::MagneticFieldContext fieldctx;
0046     std::shared_ptr<const Acts::DD4hepFieldAdapter> field_provider = std::make_shared<const Acts::DD4hepFieldAdapter>(detector.field());
0047     Acts::MagneticFieldProvider::Cache field_cache = field_provider->makeCache(fieldctx);
0048
0049     // Stepper and Propagator

0050     using Stepper    = Acts::EigenStepper<>;
0051     using Propagator = Acts::Propagator<Stepper>;
0052
0053     Stepper stepper(field_provider);
0054     Propagator propagator(stepper);
0055
0056     // Create Impact Point Estimator

0057     Acts::ImpactPointEstimator::Config ImPoEs_cfg(field_provider,std::make_shared<Propagator>(propagator));
0058
0059     Acts::ImpactPointEstimator::State ImPoEs_state;
0060     ImPoEs_state.fieldCache = field_cache;
0061
0062     Acts::ImpactPointEstimator ImPoEs(ImPoEs_cfg);
0063
0064     // Create 'vertex' at particle's creation point -- which is (x,y,z) = (1,0,0) mm

0065     Acts::Vector3 vtx_pos(1.0 * Acts::UnitConstants::mm, 0, 0);
0066
0067     // Create another 'vertex' at (2,0,0) mm and check the distance at the DCA to this point

0068     Acts::Vector3 vtx_pos2(2.0 * Acts::UnitConstants::mm, 0, 0);
0069
0070     // Define Style

0071     gStyle->SetOptStat(0);
0072     gStyle->SetPadBorderMode(0);
0073     gStyle->SetFrameBorderMode(0);
0074     gStyle->SetFrameLineWidth(2);
0075     gStyle->SetLabelSize(0.035,"X");
0076     gStyle->SetLabelSize(0.035,"Y");
0077     //gStyle->SetLabelOffset(0.01,"X");

0078     //gStyle->SetLabelOffset(0.01,"Y");

0079     gStyle->SetTitleXSize(0.04);
0080     gStyle->SetTitleXOffset(0.9);
0081     gStyle->SetTitleYSize(0.04);
0082     gStyle->SetTitleYOffset(0.9);
0083
0084     // Define histograms

0085     TH2 *h1 = new TH2D("h1","Single particle generated at (x,y,z) = (1,0,0) mm",100,-1.5,1.5,100,-1.5,1.5);
0086     h1->GetXaxis()->SetTitle("Track global x at beamline (z-axis) POCA [mm]");h1->GetXaxis()->CenterTitle();
0087     h1->GetYaxis()->SetTitle("Track global y at beamline (z-axis) POCA [mm]");h1->GetYaxis()->CenterTitle();
0088
0089     TH1 *h2 = new TH1D("h2","Single particle generated at (x,y,z) = (1,0,0) mm",100,-0.02,1);
0090     h2->GetXaxis()->SetTitle("Track distance to generation point at 3D DCA [mm]");h2->GetXaxis()->CenterTitle();
0091     h2->SetLineWidth(2);h2->SetLineColor(kBlue);
0092
0093     TH1 *h3 = new TH1D("h3","Single particle generated at (x,y,z) = (1,0,0) mm",100,0,2);
0094     h3->GetXaxis()->SetTitle("Track distance to (x,y,z) = (2,0,0) mm at 3D DCA [mm]");h3->GetXaxis()->CenterTitle();
0095     h3->SetLineWidth(2);h3->SetLineColor(kBlue);
0096
0097     TH1 *h3a = new TH1D("h3a","Single particle generated at (x,y,z) = (1,0,0) mm",100,0,2);
0098     h3a->GetXaxis()->SetTitle("Track xy distance to (x,y,z) = (2,0,0) mm at 3D DCA [mm]");h3a->GetXaxis()->CenterTitle();
0099     h3a->SetLineWidth(2);h3a->SetLineColor(kBlue);
0100
0101     TH1 *h3b = new TH1D("h3b","Single particle generated at (x,y,z) = (1,0,0) mm",100,-1,1);
0102     h3b->GetXaxis()->SetTitle("Track z signed distance to (x,y,z) = (2,0,0) mm at 3D DCA [mm]");h3b->GetXaxis()->CenterTitle();
0103     h3b->SetLineWidth(2);h3b->SetLineColor(kBlue);
0104
0105     TH1 *h4a = new TH2D("h4a","Single particle generated at (x,y,z) = (1,0,0) mm",100,-3.2,3.2,100,0,2);
0106     h4a->GetXaxis()->SetTitle("Generated particle #phi [Rad]");h4a->GetXaxis()->CenterTitle();
0107     h4a->GetYaxis()->SetTitle("Track xy distance to (x,y,z) = (2,0,0) mm at 3D DCA [mm]");h4a->GetYaxis()->CenterTitle();
0108
0109     TH1 *h4b = new TH2D("h4b","Single particle generated at (x,y,z) = (1,0,0) mm",100,0,3.2,100,-1,1);
0110     h4b->GetXaxis()->SetTitle("Generated particle #theta [Rad]");h4b->GetXaxis()->CenterTitle();
0111     h4b->GetYaxis()->SetTitle("Track z signed distance to (x,y,z) = (2,0,0) mm at 3D DCA [mm]");h4b->GetYaxis()->CenterTitle();
0112
0113     // Load input ROOT file

0114     TString run_name;
0115     TString path = "./input/";
0116     run_name = "eicrecon_out_1_0_0.root"; //Single negative muons generated at (1,0,0) mm

0117
0118     // Define variables

0119     int pid_code = 13;
0120     std::string coll = "CentralCKFTrackParameters"; //Real-seeded track parameters

0121     TLorentzVector gen_vec; //Generated particle four momentum

0122
0123     TString input = path + run_name;
0124     TFile *f = new TFile(input.Data());
0125     TTree *tree = (TTree*) f->Get("events");
0126
0127     //Create Array Reader

0128     TTreeReader tr(tree);
0129
0130     TTreeReaderArray<int>   gen_status(tr, "MCParticles.generatorStatus");
0131     TTreeReaderArray<int>   gen_pid(tr, "MCParticles.PDG");
0132     TTreeReaderArray<float> gen_px(tr, "MCParticles.momentum.x");
0133     TTreeReaderArray<float> gen_py(tr, "MCParticles.momentum.y");
0134     TTreeReaderArray<float> gen_pz(tr, "MCParticles.momentum.z");
0135     TTreeReaderArray<double> gen_mass(tr, "MCParticles.mass");
0136     TTreeReaderArray<float> gen_charge(tr, "MCParticles.charge");
0137     TTreeReaderArray<double> gen_vx(tr, "MCParticles.vertex.x");
0138     TTreeReaderArray<double> gen_vy(tr, "MCParticles.vertex.y");
0139     TTreeReaderArray<double> gen_vz(tr, "MCParticles.vertex.z");
0140
0141     TTreeReaderArray<float> track_qoverp(tr, Form("%s.qOverP",coll.c_str()));
0142     TTreeReaderArray<float> track_theta(tr, Form("%s.theta",coll.c_str()));
0143     TTreeReaderArray<float> track_phi(tr, Form("%s.phi",coll.c_str()));
0144     TTreeReaderArray<float> track_loca(tr, Form("%s.loc.a",coll.c_str()));
0145     TTreeReaderArray<float> track_locb(tr, Form("%s.loc.b",coll.c_str()));
0146     TTreeReaderArray<float> track_time(tr, Form("%s.time",coll.c_str()));
0147
0148     //Loop over events

0149     int counter(0);
0150     while (tr.Next()) {
0151
0152         if(counter%100==0) cout<<"Analyzing event "<<counter<<endl;
0153         counter++;
0154
0155         //Loop over generated particles, select primary particle (assuming single particle)

0156         for(size_t igen=0;igen<gen_status.GetSize();igen++){
0157
0158             if(gen_status[igen]==1 && gen_pid[igen]==pid_code){ //PID code requirement not really needed...

0159                 gen_vec.SetXYZM(gen_px[igen],gen_py[igen],gen_pz[igen],gen_mass[igen]);
0160                 break;
0161             }
0162         } //End loop over generated particles

0163
0164         //Loop over tracks

0165         size_t track_mult = track_qoverp.GetSize();
0166         for(size_t itrack=0;itrack<track_mult;itrack++){
0167
0168             auto loc_a = track_loca[itrack];
0169             auto loc_b = track_locb[itrack];
0170             auto phi = track_phi[itrack];
0171             auto theta = track_theta[itrack];
0172             auto qoverP = track_qoverp[itrack];
0173             auto time = track_time[itrack];
0174
0175             // Create BoundTrackParamters

0176             Acts::BoundVector params;
0177
0178         params(Acts::eBoundLoc0)   = loc_a;
0179         params(Acts::eBoundLoc1)   = loc_b;
0180         params(Acts::eBoundPhi)    = phi;
0181         params(Acts::eBoundTheta)  = theta;
0182         params(Acts::eBoundQOverP) = qoverP;
0183         params(Acts::eBoundTime)   = time;
0184
0185             //FIXME: Set covariance matrix based on input ROOT file information

0186         Acts::BoundSquareMatrix cov = Acts::BoundSquareMatrix::Zero();
0187
0188         Acts::BoundTrackParameters track_parameters(perigee,params,cov,Acts::ParticleHypothesis::pion());
0189
0190             //---- Part 1: Convert from local coordinates to global coordinates at beamline (z-axis) POCA ----

0191             Acts::Vector2 localpos( loc_a, loc_b );
0192             Acts::Vector3 direction(sin(theta)*cos(phi),
0193                                     sin(theta)*sin(phi),
0194                                     cos(theta));
0195
0196             // Convert to global coordinates for PCA at perigee surface

0197             auto global_perigee = perigee->localToGlobal(trackingGeoCtx,localpos,direction);
0198
0199             if( !global_perigee.isZero() ){
0200                 h1->Fill(global_perigee.x(),global_perigee.y());
0201             }
0202
0203             //---- Part 2: Get track parameters at 3D DCA to creation point at (x,y,z) = (1,0,0) mm ----

0204             auto result = ImPoEs.estimate3DImpactParameters(trackingGeoCtx,fieldctx,track_parameters,vtx_pos,ImPoEs_state);
0205
0206         if(result.ok()){
0207         Acts::BoundTrackParameters trk_boundpar_vtx = result.value();
0208         const auto& trk_vtx_params  = trk_boundpar_vtx.parameters();
0209
0210         h2->Fill(trk_vtx_params[Acts::eBoundLoc0]);
0211             }
0212
0213             //---- Part 2a: Get track parameters at 3D DCA to (x,y,z) = (2,0,0) mm ----

0214             auto result2 = ImPoEs.estimate3DImpactParameters(trackingGeoCtx,fieldctx,track_parameters,vtx_pos2,ImPoEs_state);
0215
0216         if(result2.ok()){
0217         Acts::BoundTrackParameters trk_boundpar_vtx2 = result2.value();
0218         const auto& trk_vtx_params2  = trk_boundpar_vtx2.parameters();
0219
0220                 h3->Fill(trk_vtx_params2[Acts::eBoundLoc0]);
0221
0222                 //Get global position at 3D DCA

0223                 auto trk_vtx2_gbl_pos = trk_boundpar_vtx2.position(trackingGeoCtx);
0224
0225         auto delta_x = trk_vtx2_gbl_pos.x() - 2.;
0226                 auto delta_y = trk_vtx2_gbl_pos.y() - 0.;
0227                 auto delta_z = trk_vtx2_gbl_pos.z() - 0.;
0228
0229                 auto xy_distance = std::hypot(delta_x,delta_y);
0230
0231                 h3a->Fill(xy_distance);
0232                 h3b->Fill(delta_z);
0233                 h4a->Fill(gen_vec.Phi(),xy_distance);
0234                 h4b->Fill(gen_vec.Theta(),delta_z);
0235             }
0236
0237             //---- Part 3: Propagate track to Perigee surface at (x,y,z) = (2,0,0) mm ----

0238
0239             // Define Perigee surface to which to propagate track

0240             auto perigee2 = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3(2,0,0));
0241
0242             // Create propagator options

0243             Acts::PropagatorOptions<> pOptions(trackingGeoCtx, fieldctx);
0244             auto intersection = perigee2->intersect(trackingGeoCtx, track_parameters.position(trackingGeoCtx),
0245                                 track_parameters.direction(), Acts::BoundaryCheck(false)).closest();
0246
0247             pOptions.direction = Acts::Direction::fromScalarZeroAsPositive(intersection.pathLength());
0248
0249             // Do the propagation to linPoint

0250             auto result_perigee2 = propagator.propagateToSurface(track_parameters, *perigee2, pOptions);
0251
0252             if(result_perigee2.ok()){
0253                 Acts::BoundTrackParameters trk_boundpar_perigee2 = result2.value();
0254                 const auto& trk_perigee2  = trk_boundpar_perigee2.parameters();
0255              }
0256
0257         } // End loop over tracks

0258
0259     } // End loop over events

0260
0261     // Make plots

0262     TCanvas *c1 = new TCanvas("c1");
0263     h1->Draw();
0264
0265     TCanvas *c2 = new TCanvas("c2");
0266     h2->Draw();
0267
0268     TCanvas *c3 = new TCanvas("c3");
0269     h3->Draw();
0270
0271     TCanvas *c3a = new TCanvas("c3a");
0272     h3a->Draw();
0273
0274     TCanvas *c3b = new TCanvas("c3b");
0275     h3b->Draw();
0276
0277     TCanvas *c4a = new TCanvas("c4a");
0278     h4a->Draw();
0279
0280     TCanvas *c4b = new TCanvas("c4b");
0281     h4b->Draw();
0282
0283     //Print plots to file

0284     c1->Print("plots/pca_global_impactpoint.pdf[");
0285     c1->Print("plots/pca_global_impactpoint.pdf");
0286     c2->Print("plots/pca_global_impactpoint.pdf");
0287     c3->Print("plots/pca_global_impactpoint.pdf");
0288     c3a->Print("plots/pca_global_impactpoint.pdf");
0289     c3b->Print("plots/pca_global_impactpoint.pdf");
0290     c4a->Print("plots/pca_global_impactpoint.pdf");
0291     c4b->Print("plots/pca_global_impactpoint.pdf");
0292     c3b->Print("plots/pca_global_impactpoint.pdf]");
0293 }