JetsAndHF/SecondaryVertex_Chi2/pca_global_impactpoint.C

0001 // Example code provide by Barak Schmookler
0002 // for calculating DCA to a specific point
0003
0004 #include "DD4hep/Detector.h"
0005 #include "DD4hep/DetElement.h"
0006
0007 #include "Acts/Plugins/DD4hep/ConvertDD4hepDetector.hpp"
0008 #include "Acts/Plugins/DD4hep/DD4hepFieldAdapter.hpp"
0009 #include "Acts/Geometry/TrackingGeometry.hpp"
0010 #include "Acts/Geometry/GeometryContext.hpp"
0011 #include "Acts/Plugins/DD4hep/ConvertDD4hepDetector.hpp"
0012 #include "Acts/Utilities/Logger.hpp"
0013 #include "Acts/Utilities/BinningType.hpp"
0014 #include "Acts/Definitions/Units.hpp"
0015 #include "Acts/Surfaces/PerigeeSurface.hpp"
0016 #include "Acts/EventData/TrackParameters.hpp"
0017 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
0018 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0019 #include "Acts/Propagator/EigenStepper.hpp"
0020 #include "Acts/Propagator/Propagator.hpp"
0021 #include "Acts/Vertexing/ImpactPointEstimator.hpp"
0022
0023 void pca_global_impactpoint(){
0024
0025     // Load DD4Hep geometry
0026     dd4hep::Detector& detector = dd4hep::Detector::getInstance();
0027     detector.fromCompact("/opt/detector/epic-main/share/epic/epic_craterlake.xml");
0028     dd4hep::DetElement geometry = detector.world();
0029
0030     // Convert DD4Hep geometry to tracking geometry
0031     Acts::GeometryContext trackingGeoCtx;
0032     auto logger = Acts::getDefaultLogger("DD4hepConversion", Acts::Logging::Level::INFO);
0033     Acts::BinningType bTypePhi = Acts::equidistant;
0034     Acts::BinningType bTypeR = Acts::equidistant;
0035     Acts::BinningType bTypeZ = Acts::equidistant;
0036     double layerEnvelopeR = Acts::UnitConstants::mm;
0037     double layerEnvelopeZ = Acts::UnitConstants::mm;
0038     double defaultLayerThickness = Acts::UnitConstants::fm;
0039     using Acts::sortDetElementsByID;
0040
0041     std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry{nullptr};
0042     trackingGeometry = Acts::convertDD4hepDetector(geometry,*logger,bTypePhi,bTypeR,bTypeZ,layerEnvelopeR,layerEnvelopeZ,defaultLayerThickness,sortDetElementsByID,trackingGeoCtx);
0043
0044     // Define Perigee surface at which reconstructed track parameters are set
0045     auto perigee = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3(0,0,0));
0046
0047     // Get Magnetic field context
0048     Acts::MagneticFieldContext fieldctx;
0049     std::shared_ptr<const Acts::DD4hepFieldAdapter> field_provider = std::make_shared<const Acts::DD4hepFieldAdapter>(detector.field());
0050     Acts::MagneticFieldProvider::Cache field_cache = field_provider->makeCache(fieldctx);
0051
0052     // Stepper and Propagator
0053     using Stepper    = Acts::EigenStepper<>;
0054     using Propagator = Acts::Propagator<Stepper>;
0055
0056     Stepper stepper(field_provider);
0057     Propagator propagator(stepper);
0058
0059     // Create Impact Point Estimator
0060     Acts::ImpactPointEstimator::Config ImPoEs_cfg(field_provider,std::make_shared<Propagator>(propagator));
0061
0062     Acts::ImpactPointEstimator::State ImPoEs_state;
0063     ImPoEs_state.fieldCache = field_cache;
0064
0065     Acts::ImpactPointEstimator ImPoEs(ImPoEs_cfg);
0066
0067     // Create 'vertex' at particle's creation point -- which is (x,y,z) = (1,0,0) mm
0068     Acts::Vector3 vtx_pos(1.0 * Acts::UnitConstants::mm, 0, 0);
0069
0070     // Create another 'vertex' at (2,0,0) mm and check the distance at the DCA to this point
0071     Acts::Vector3 vtx_pos2(2.0 * Acts::UnitConstants::mm, 0, 0);
0072
0073     // Define Style
0074     gStyle->SetOptStat(0);
0075     gStyle->SetPadBorderMode(0);
0076     gStyle->SetFrameBorderMode(0);
0077     gStyle->SetFrameLineWidth(2);
0078     gStyle->SetLabelSize(0.035,"X");
0079     gStyle->SetLabelSize(0.035,"Y");
0080     //gStyle->SetLabelOffset(0.01,"X");
0081     //gStyle->SetLabelOffset(0.01,"Y");
0082     gStyle->SetTitleXSize(0.04);
0083     gStyle->SetTitleXOffset(0.9);
0084     gStyle->SetTitleYSize(0.04);
0085     gStyle->SetTitleYOffset(0.9);
0086
0087     // Define histograms
0088     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);
0089     h1->GetXaxis()->SetTitle("Track global x at beamline (z-axis) POCA [mm]");h1->GetXaxis()->CenterTitle();
0090     h1->GetYaxis()->SetTitle("Track global y at beamline (z-axis) POCA [mm]");h1->GetYaxis()->CenterTitle();
0091
0092     TH1 *h2 = new TH1D("h2","Single particle generated at (x,y,z) = (1,0,0) mm",100,-0.02,1);
0093     h2->GetXaxis()->SetTitle("Track distance to generation point at 3D DCA [mm]");h2->GetXaxis()->CenterTitle();
0094     h2->SetLineWidth(2);h2->SetLineColor(kBlue);
0095
0096     TH1 *h3 = new TH1D("h3","Single particle generated at (x,y,z) = (1,0,0) mm",100,0,2);
0097     h3->GetXaxis()->SetTitle("Track distance to (x,y,z) = (2,0,0) mm at 3D DCA [mm]");h3->GetXaxis()->CenterTitle();
0098     h3->SetLineWidth(2);h3->SetLineColor(kBlue);
0099
0100     TH1 *h3a = new TH1D("h3a","Single particle generated at (x,y,z) = (1,0,0) mm",100,0,2);
0101     h3a->GetXaxis()->SetTitle("Track xy distance to (x,y,z) = (2,0,0) mm at 3D DCA [mm]");h3a->GetXaxis()->CenterTitle();
0102     h3a->SetLineWidth(2);h3a->SetLineColor(kBlue);
0103
0104     TH1 *h3b = new TH1D("h3b","Single particle generated at (x,y,z) = (1,0,0) mm",100,-1,1);
0105     h3b->GetXaxis()->SetTitle("Track z signed distance to (x,y,z) = (2,0,0) mm at 3D DCA [mm]");h3b->GetXaxis()->CenterTitle();
0106     h3b->SetLineWidth(2);h3b->SetLineColor(kBlue);
0107
0108     TH1 *h4a = new TH2D("h4a","Single particle generated at (x,y,z) = (1,0,0) mm",100,-3.2,3.2,100,0,2);
0109     h4a->GetXaxis()->SetTitle("Generated particle #phi [Rad]");h4a->GetXaxis()->CenterTitle();
0110     h4a->GetYaxis()->SetTitle("Track xy distance to (x,y,z) = (2,0,0) mm at 3D DCA [mm]");h4a->GetYaxis()->CenterTitle();
0111
0112     TH1 *h4b = new TH2D("h4b","Single particle generated at (x,y,z) = (1,0,0) mm",100,0,3.2,100,-1,1);
0113     h4b->GetXaxis()->SetTitle("Generated particle #theta [Rad]");h4b->GetXaxis()->CenterTitle();
0114     h4b->GetYaxis()->SetTitle("Track z signed distance to (x,y,z) = (2,0,0) mm at 3D DCA [mm]");h4b->GetYaxis()->CenterTitle();
0115
0116     // Load input ROOT file
0117     TString run_name;
0118     TString path = "./input/";
0119     run_name = "eicrecon_out_1_0_0.root"; //Single negative muons generated at (1,0,0) mm
0120
0121     // Define variables
0122     int pid_code = 13;
0123     std::string coll = "CentralCKFTrackParameters"; //Real-seeded track parameters
0124     TLorentzVector gen_vec; //Generated particle four momentum
0125
0126     TString input = path + run_name;
0127     TFile *f = new TFile(input.Data());
0128     TTree *tree = (TTree*) f->Get("events");
0129
0130     //Create Array Reader
0131     TTreeReader tr(tree);
0132
0133     TTreeReaderArray<int>   gen_status(tr, "MCParticles.generatorStatus");
0134     TTreeReaderArray<int>   gen_pid(tr, "MCParticles.PDG");
0135     TTreeReaderArray<float> gen_px(tr, "MCParticles.momentum.x");
0136     TTreeReaderArray<float> gen_py(tr, "MCParticles.momentum.y");
0137     TTreeReaderArray<float> gen_pz(tr, "MCParticles.momentum.z");
0138     TTreeReaderArray<double> gen_mass(tr, "MCParticles.mass");
0139     TTreeReaderArray<float> gen_charge(tr, "MCParticles.charge");
0140     TTreeReaderArray<double> gen_vx(tr, "MCParticles.vertex.x");
0141     TTreeReaderArray<double> gen_vy(tr, "MCParticles.vertex.y");
0142     TTreeReaderArray<double> gen_vz(tr, "MCParticles.vertex.z");
0143
0144     TTreeReaderArray<float> track_qoverp(tr, Form("%s.qOverP",coll.c_str()));
0145     TTreeReaderArray<float> track_theta(tr, Form("%s.theta",coll.c_str()));
0146     TTreeReaderArray<float> track_phi(tr, Form("%s.phi",coll.c_str()));
0147     TTreeReaderArray<float> track_loca(tr, Form("%s.loc.a",coll.c_str()));
0148     TTreeReaderArray<float> track_locb(tr, Form("%s.loc.b",coll.c_str()));
0149
0150     //Loop over events
0151     int counter(0);
0152     while (tr.Next()) {
0153
0154         if(counter%100==0) cout<<"Analyzing event "<<counter<<endl;
0155         counter++;
0156
0157         //Loop over generated particles, select primary particle (assuming single particle)
0158         for(size_t igen=0;igen<gen_status.GetSize();igen++){
0159
0160             if(gen_status[igen]==1 && gen_pid[igen]==pid_code){ //PID code requirement not really needed...
0161                 gen_vec.SetXYZM(gen_px[igen],gen_py[igen],gen_pz[igen],gen_mass[igen]);
0162                 break;
0163             }
0164         } //End loop over generated particles
0165
0166         //Loop over tracks
0167         size_t track_mult = track_qoverp.GetSize();
0168         for(size_t itrack=0;itrack<track_mult;itrack++){
0169
0170             auto loc_a = track_loca[itrack];
0171             auto loc_b = track_locb[itrack];
0172             auto phi = track_phi[itrack];
0173             auto theta = track_theta[itrack];
0174             auto qoverP = track_qoverp[itrack];
0175
0176             // Create BoundTrackParamters
0177             Acts::BoundVector params;
0178
0179             params(Acts::eBoundLoc0)   = loc_a;
0180             params(Acts::eBoundLoc1)   = loc_b;
0181             params(Acts::eBoundPhi)    = phi;
0182             params(Acts::eBoundTheta)  = theta;
0183             params(Acts::eBoundQOverP) = qoverP;
0184             params(Acts::eBoundTime)   = 0;
0185
0186             //FIXME: Set covariance matrix based on input ROOT file information
0187             Acts::BoundSquareMatrix cov = Acts::BoundSquareMatrix::Zero();
0188
0189             Acts::BoundTrackParameters track_parameters(perigee,params,cov,Acts::ParticleHypothesis::pion());
0190
0191             //---- Part 1: Convert from local coordinates to global coordinates at beamline (z-axis) POCA ----
0192             Acts::Vector2 localpos( loc_a, loc_b );
0193             Acts::Vector3 direction(sin(theta)*cos(phi),
0194                                     sin(theta)*sin(phi),
0195                                     cos(theta));
0196
0197             // Convert to global coordinates for PCA at perigee surface
0198             auto global_perigee = perigee->localToGlobal(trackingGeoCtx,localpos,direction);
0199
0200             if( !global_perigee.isZero() ){
0201                 h1->Fill(global_perigee.x(),global_perigee.y());
0202             }
0203
0204             //---- Part 2: Get track parameters at 3D DCA to creation point at (x,y,z) = (1,0,0) mm ----
0205             auto result = ImPoEs.estimate3DImpactParameters(trackingGeoCtx,fieldctx,track_parameters,vtx_pos,ImPoEs_state);
0206
0207             if(result.ok()){
0208                 Acts::BoundTrackParameters trk_boundpar_vtx = result.value();
0209                 const auto& trk_vtx_params  = trk_boundpar_vtx.parameters();
0210
0211                 h2->Fill(trk_vtx_params[Acts::eBoundLoc0]);
0212             }
0213
0214             //---- Part 2a: Get track parameters at 3D DCA to (x,y,z) = (2,0,0) mm ----
0215             auto result2 = ImPoEs.estimate3DImpactParameters(trackingGeoCtx,fieldctx,track_parameters,vtx_pos2,ImPoEs_state);
0216
0217             if(result2.ok()){
0218                 Acts::BoundTrackParameters trk_boundpar_vtx2 = result2.value();
0219                 const auto& trk_vtx_params2  = trk_boundpar_vtx2.parameters();
0220
0221                 h3->Fill(trk_vtx_params2[Acts::eBoundLoc0]);
0222
0223                 //Get global position at 3D DCA
0224                 auto trk_vtx2_gbl_pos = trk_boundpar_vtx2.position(trackingGeoCtx);
0225
0226                 auto delta_x = trk_vtx2_gbl_pos.x() - 2.;
0227                 auto delta_y = trk_vtx2_gbl_pos.y() - 0.;
0228                 auto delta_z = trk_vtx2_gbl_pos.z() - 0.;
0229
0230                 auto xy_distance = std::hypot(delta_x,delta_y);
0231
0232                 h3a->Fill(xy_distance);
0233                 h3b->Fill(delta_z);
0234                 h4a->Fill(gen_vec.Phi(),xy_distance);
0235                 h4b->Fill(gen_vec.Theta(),delta_z);
0236             }
0237
0238         } // End loop over tracks
0239
0240     } // End loop over events
0241
0242     // Make plots
0243     TCanvas *c1 = new TCanvas("c1");
0244     h1->Draw();
0245
0246     TCanvas *c2 = new TCanvas("c2");
0247     h2->Draw();
0248
0249     TCanvas *c3 = new TCanvas("c3");
0250     h3->Draw();
0251
0252     TCanvas *c3a = new TCanvas("c3a");
0253     h3a->Draw();
0254
0255     TCanvas *c3b = new TCanvas("c3b");
0256     h3b->Draw();
0257
0258     TCanvas *c4a = new TCanvas("c4a");
0259     h4a->Draw();
0260
0261     TCanvas *c4b = new TCanvas("c4b");
0262     h4b->Draw();
0263
0264     //Print plots to file
0265     c1->Print("plots/pca_global_impactpoint.pdf[");
0266     c1->Print("plots/pca_global_impactpoint.pdf");
0267     c2->Print("plots/pca_global_impactpoint.pdf");
0268     c3->Print("plots/pca_global_impactpoint.pdf");
0269     c3a->Print("plots/pca_global_impactpoint.pdf");
0270     c3b->Print("plots/pca_global_impactpoint.pdf");
0271     c4a->Print("plots/pca_global_impactpoint.pdf");
0272     c4b->Print("plots/pca_global_impactpoint.pdf");
0273     c3b->Print("plots/pca_global_impactpoint.pdf]");
0274 }