diff --git a/Detector/DetCRD/compact/CRD_common_v01/DC_Simple_v01_03.xml b/Detector/DetCRD/compact/CRD_common_v01/DC_Simple_v01_03.xml new file mode 100644 index 000000000..e3c372fdd --- /dev/null +++ b/Detector/DetCRD/compact/CRD_common_v01/DC_Simple_v01_03.xml @@ -0,0 +1,79 @@ + + + + + Test with Drift Chamber + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:5,layer:7:9,chamber:8,cellID:32:16 + + + + diff --git a/Detector/DetCRD/compact/CRD_common_v01/DC_Simple_v01_05.xml b/Detector/DetCRD/compact/CRD_common_v01/DC_Simple_v01_05.xml new file mode 100644 index 000000000..383019ecb --- /dev/null +++ b/Detector/DetCRD/compact/CRD_common_v01/DC_Simple_v01_05.xml @@ -0,0 +1,79 @@ + + + + + Test with Drift Chamber + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:5,layer:7:9,chamber:8,cellID:32:16 + + + + diff --git a/Detector/DetCRD/compact/CRD_common_v01/DC_Stero_v01_01.xml b/Detector/DetCRD/compact/CRD_common_v01/DC_Stero_v01_01.xml new file mode 100644 index 000000000..fe303a521 --- /dev/null +++ b/Detector/DetCRD/compact/CRD_common_v01/DC_Stero_v01_01.xml @@ -0,0 +1,89 @@ + + + + + Test with Drift Chamber + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:5,layer:7:9,chamber:8,cellID:32:16 + + + + diff --git a/Detector/DetCRD/compact/CRD_common_v01/DC_Simple_v01_01.xml b/Detector/DetCRD/compact/CRD_common_v01/DC_Straight_v01_01.xml similarity index 86% rename from Detector/DetCRD/compact/CRD_common_v01/DC_Simple_v01_01.xml rename to Detector/DetCRD/compact/CRD_common_v01/DC_Straight_v01_01.xml index a1da46679..73594f243 100644 --- a/Detector/DetCRD/compact/CRD_common_v01/DC_Simple_v01_01.xml +++ b/Detector/DetCRD/compact/CRD_common_v01/DC_Straight_v01_01.xml @@ -28,7 +28,7 @@ - + @@ -42,7 +42,7 @@ - + @@ -53,8 +53,10 @@ - - + + + + @@ -78,7 +80,7 @@ - + system:5,layer:7:9,chamber:8,cellID:32:16 diff --git a/Detector/DetCRD/compact/CRD_common_v01/DC_Simple_v01_02.xml b/Detector/DetCRD/compact/CRD_common_v01/DC_Straight_v01_02.xml similarity index 88% rename from Detector/DetCRD/compact/CRD_common_v01/DC_Simple_v01_02.xml rename to Detector/DetCRD/compact/CRD_common_v01/DC_Straight_v01_02.xml index ae4ccbe02..0b6e8f67e 100644 --- a/Detector/DetCRD/compact/CRD_common_v01/DC_Simple_v01_02.xml +++ b/Detector/DetCRD/compact/CRD_common_v01/DC_Straight_v01_02.xml @@ -28,7 +28,7 @@ - + @@ -54,7 +54,9 @@ - + + + @@ -79,7 +81,7 @@ - + system:5,layer:7:9,chamber:8,cellID:32:16 diff --git a/Detector/DetCRD/compact/CRD_o1_v01/CRD_Dimensions_v01_01.xml b/Detector/DetCRD/compact/CRD_o1_v01/CRD_Dimensions_v01_01.xml index cd1e9b2e5..9a79f6fb8 100644 --- a/Detector/DetCRD/compact/CRD_o1_v01/CRD_Dimensions_v01_01.xml +++ b/Detector/DetCRD/compact/CRD_o1_v01/CRD_Dimensions_v01_01.xml @@ -84,12 +84,16 @@ + + + + - - - + + + @@ -99,13 +103,6 @@ - - - - - - - diff --git a/Detector/DetCRD/compact/CRD_o1_v01/CRD_o1_v01-onlyTracker.xml b/Detector/DetCRD/compact/CRD_o1_v01/CRD_o1_v01-onlyTracker.xml index b35e32655..d80301fb2 100644 --- a/Detector/DetCRD/compact/CRD_o1_v01/CRD_o1_v01-onlyTracker.xml +++ b/Detector/DetCRD/compact/CRD_o1_v01/CRD_o1_v01-onlyTracker.xml @@ -31,7 +31,7 @@ - + diff --git a/Detector/DetCRD/compact/CRD_o1_v02/CRD_Dimensions_v01_02.xml b/Detector/DetCRD/compact/CRD_o1_v02/CRD_Dimensions_v01_02.xml index cd1e9b2e5..9a79f6fb8 100644 --- a/Detector/DetCRD/compact/CRD_o1_v02/CRD_Dimensions_v01_02.xml +++ b/Detector/DetCRD/compact/CRD_o1_v02/CRD_Dimensions_v01_02.xml @@ -84,12 +84,16 @@ + + + + - - - + + + @@ -99,13 +103,6 @@ - - - - - - - diff --git a/Detector/DetDriftChamber/compact/det.xml b/Detector/DetDriftChamber/compact/det.xml index fb4bffdee..22548148c 100644 --- a/Detector/DetDriftChamber/compact/det.xml +++ b/Detector/DetDriftChamber/compact/det.xml @@ -1,15 +1,13 @@ - + - Test with Drift Chamber + version="v2"> + Drift Chamber @@ -17,51 +15,49 @@ - + - + - - + + + + + - - + + - - + - - - - + + + + + + - - - + + - - - + + - - - - + + - - @@ -76,7 +72,7 @@ - + @@ -86,11 +82,13 @@ - - + + + + - + @@ -112,7 +110,7 @@ - + system:5,layer:7:9,chamber:8,cellID:32:16 diff --git a/Detector/DetDriftChamber/compact/det_stero.xml b/Detector/DetDriftChamber/compact/det_stero.xml new file mode 100644 index 000000000..44e919087 --- /dev/null +++ b/Detector/DetDriftChamber/compact/det_stero.xml @@ -0,0 +1,129 @@ + + + + + Test with Drift Chamber + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:5,layer:7:9,chamber:8,cellID:32:16 + + + + + + + + + + + + diff --git a/Detector/DetDriftChamber/src/driftchamber/DriftChamber.cpp b/Detector/DetDriftChamber/src/driftchamber/DriftChamber.cpp index be98c1c14..1ffdd929e 100644 --- a/Detector/DetDriftChamber/src/driftchamber/DriftChamber.cpp +++ b/Detector/DetDriftChamber/src/driftchamber/DriftChamber.cpp @@ -16,6 +16,9 @@ #include "DDSegmentation/Segmentation.h" #include "DetSegmentation/GridDriftChamber.h" +#include +#include + using namespace dd4hep; using namespace dd4hep::detail; using namespace dd4hep::rec ; @@ -27,40 +30,57 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector, xml_h e, dd4hep::SensitiveDetector sens) { // ------- Lambda functions ---- // - auto delta_a_func = [](auto x, auto y) { return 0.5 * ( x + y ); }; + auto delta_b_func = [](auto x, auto y) { return 2 * std::sqrt((x + y) * (x - y)); }; + auto epsilon_func = [](auto delta_a, auto L) { return std::atan(delta_a / L); }; // ======================================================================= // Parameter Definition // ======================================================================= + auto start = std::chrono::steady_clock::now(); + xml_det_t x_det = e; xml_coll_t c(x_det,_U(chamber)); xml_comp_t x_chamber = c; + xml_coll_t cc(x_det,_U(side)); + xml_comp_t x_side = cc; + std::string det_name = x_det.nameStr(); std::string det_type = x_det.typeStr(); dd4hep::SensitiveDetector sd = sens; // - global - double chamber_half_length = theDetector.constant("DC_half_length"); + double chamber_half_length = theDetector.constant("Gas_half_length"); + double chamber_length = chamber_half_length*2; - // - chamber - double chamber_radius_min = theDetector.constant("SDT_chamber_radius_min"); - double chamber_radius_max = theDetector.constant("SDT_chamber_radius_max"); - double SDT_half_length = theDetector.constant("SDT_chamber_half_length"); + double alpha = theDetector.constant("Alpha"); + double Pi = 0; + if(0!=alpha) Pi = 0.5*M_PI; + + double safe_distance = theDetector.constant("DC_safe_distance"); + + double DC_inner_wall_thickness = theDetector.constant("DC_inner_wall_thickness"); + double DC_outer_wall_thickness = theDetector.constant("DC_outer_wall_thickness"); + + // - chamber gas + double DC_rend = theDetector.constant("DC_rend"); + + double chamber_radius_min = theDetector.constant("Gas_radius_min"); + double chamber_radius_max = DC_rend-safe_distance-DC_outer_wall_thickness; + double layer_radius_max = chamber_radius_max*std::cos(alpha); int chamberID = x_chamber.id(); // - layer - double chamber_layer_width = theDetector.constant("SDT_chamber_layer_width"); - double chamber_cell_width = theDetector.constant("SDT_chamber_cell_width"); - double chamber_layer_rbegin = theDetector.constant("DC_chamber_layer_rbegin"); - double chamber_layer_rend = theDetector.constant("DC_chamber_layer_rend"); - int chamber_layer_number = floor((chamber_layer_rend-chamber_layer_rbegin)/chamber_layer_width); + int DC_layer_number = theDetector.constant("DC_layer_number"); - double safe_distance = theDetector.constant("DC_safe_distance"); - double epsilon = theDetector.constant("Epsilon"); + double layer_width = (layer_radius_max-chamber_radius_min)/DC_layer_number; + + double cell_width = theDetector.constant("DC_cell_width"); + + int DC_construct_wire = theDetector.constant("DC_construct_wire"); // ======================================================================= // Detector Construction @@ -74,106 +94,112 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector, if( theDetector.buildType() == BUILD_ENVELOPE ) return sdet ; + dd4hep::Material det_mat(theDetector.material(x_det.materialStr())); + dd4hep::Material chamber_mat = theDetector.material(x_chamber.materialStr()); - dd4hep::Material det_mat(theDetector.material("Air")); - dd4hep::Material chamber_mat(theDetector.material("GasHe_90Isob_10")); - - // - global + /// - global Assembly det_vol( det_name+"_assembly" ) ; - // - chamber volume + /// - chamber volume dd4hep::Tube det_chamber_solid(chamber_radius_min, chamber_radius_max, chamber_half_length); dd4hep::Volume det_chamber_vol(det_name+"_chamber_vol", det_chamber_solid, chamber_mat); - // - wall - double chamber_inner_wall_rmin = theDetector.constant("SDT_chamber_inner_wall_radius_min"); - double chamber_inner_wall_rmax = theDetector.constant("SDT_chamber_inner_wall_radius_max"); - double chamber_outer_wall_rmin = theDetector.constant("SDT_chamber_outer_wall_radius_min"); - double chamber_outer_wall_rmax = theDetector.constant("SDT_chamber_outer_wall_radius_max"); + /// - wall + double chamber_inner_wall_rmin = theDetector.constant("DC_inner_wall_radius_min"); + double chamber_inner_wall_rmax = theDetector.constant("DC_inner_wall_radius_max"); + double chamber_outer_wall_rmin = chamber_radius_max+safe_distance; + double chamber_outer_wall_rmax = chamber_outer_wall_rmin+DC_outer_wall_thickness; - dd4hep::Material wall_mat(theDetector.material("CarbonFiber")); + dd4hep::Material wall_mat(theDetector.material(x_side.materialStr())); double wall_rmin[2] = {chamber_inner_wall_rmin, chamber_outer_wall_rmin}; double wall_rmax[2] = {chamber_inner_wall_rmax, chamber_outer_wall_rmax}; - // - wire - dd4hep::Volume module_vol; - dd4hep::Volume Module_vol; - for(xml_coll_t c(x_det,_U(module)); c; ++c) { - xml_comp_t x_module = c; - double module_rmin = x_module.rmin(); - double module_rmax = x_module.rmax(); - std::string module_name = x_module.nameStr(); - dd4hep::Tube module_solid(module_rmin,module_rmax,chamber_half_length); - if(x_module.id()==0) { - module_vol = dd4hep::Volume(module_name,module_solid,det_mat); - module_vol.setVisAttributes(theDetector.visAttributes(x_module.visStr())); - } else { - Module_vol = dd4hep::Volume(module_name,module_solid,det_mat); - Module_vol.setVisAttributes(theDetector.visAttributes(x_module.visStr())); - } - - for(xml_coll_t l(x_module,_U(tubs)); l; ++l) { - xml_comp_t x_tube =l; - double tube_rmin = x_tube.rmin(); - double tube_rmax = x_tube.rmax(); - std::string tube_name = x_tube.nameStr(); - std::string wire_name= module_name + tube_name; - dd4hep::Material tube_mat = theDetector.material(x_tube.materialStr()); - dd4hep::Tube wire_solid(tube_rmin,tube_rmax,chamber_half_length); - dd4hep::Volume wire_vol(wire_name,wire_solid,tube_mat); - dd4hep::Transform3D transform_wire(dd4hep::Rotation3D(),dd4hep::Position(0.,0.,0.)); - dd4hep::PlacedVolume wire_phy; - if(x_module.id()==0) { - wire_phy = module_vol.placeVolume(wire_vol,transform_wire); - } else { - wire_phy = Module_vol.placeVolume(wire_vol,transform_wire); - } - } - } - - // End cap - double Endcap_rmin = theDetector.constant("DC_Endcap_rmin"); - double Endcap_rmax = theDetector.constant("DC_Endcap_rmax"); - double Endcap_z = theDetector.constant("DC_Endcap_dz"); - dd4hep::Tube det_Endcap_solid(Endcap_rmin,Endcap_rmax,0.5*Endcap_z); + /// - construct wires + dd4hep::Volume signalWireVolume; + dd4hep::Volume fieldWireVolume; + for(xml_coll_t dcModule(x_det,_U(module)); dcModule; ++dcModule) { + xml_comp_t x_module = dcModule; + std::string module_name = x_module.nameStr(); + dd4hep::Tube module_solid(x_module.rmin(),x_module.rmax(),chamber_half_length); + if(0==x_module.id()) { + signalWireVolume = dd4hep::Volume(module_name,module_solid,det_mat); + signalWireVolume.setVisAttributes(theDetector.visAttributes(x_module.visStr())); + } else { + fieldWireVolume = dd4hep::Volume(module_name,module_solid,det_mat); + fieldWireVolume.setVisAttributes(theDetector.visAttributes(x_module.visStr())); + } + + /// construct wire tubes + for(xml_coll_t l(x_module,_U(tubs)); l; ++l) { + xml_comp_t x_tube =l; + std::string wire_name= module_name + x_tube.nameStr(); + dd4hep::Material tube_mat = theDetector.material(x_tube.materialStr()); + dd4hep::Tube wire_solid(x_tube.rmin(),x_tube.rmax(),chamber_half_length); + dd4hep::Volume wire_vol(wire_name,wire_solid,tube_mat); + dd4hep::Transform3D transform_wire(dd4hep::Rotation3D(),dd4hep::Position(0.,0.,0.)); + if(0==x_module.id()) { + signalWireVolume.placeVolume(wire_vol,transform_wire); + } else { + fieldWireVolume.placeVolume(wire_vol,transform_wire); + } + }//end of construct tubes + }//end of construct wire + + /// construct End cap + double endcap_rmin = theDetector.constant("DC_Endcap_rmin"); + double endcap_rmax = theDetector.constant("DC_Endcap_rmax"); + double endcap_z = theDetector.constant("DC_Endcap_dz"); + dd4hep::Tube det_Endcap_solid(endcap_rmin,endcap_rmax,0.5*endcap_z); dd4hep::Volume det_Endcap_vol(det_name+"Endcap",det_Endcap_solid,det_mat); det_Endcap_vol.setVisAttributes(theDetector,"YellowVis"); - //Initialize the segmentation - dd4hep::Readout readout = sd.readout(); - dd4hep::Segmentation geomseg = readout.segmentation(); - dd4hep::Segmentation* _geoSeg = &geomseg; + ///Initialize the segmentation + auto segmentationDC = + dynamic_cast + (sd.readout().segmentation().segmentation()); - auto DCHseg = dynamic_cast(_geoSeg->segmentation()); - - // - layer + /// - construct layers int chamber_id = 0; - int layerIndex = -1; - for(int layer_id = 0; layer_id < chamber_layer_number; layer_id++) { - double rmin,rmax,offset=0; + for(int layer_id = 0; layer_id < DC_layer_number; layer_id++) { dd4hep::Volume* current_vol_ptr = nullptr; - current_vol_ptr = &det_chamber_vol; - rmin = chamber_layer_rbegin+(layer_id*chamber_layer_width); - rmax = rmin+chamber_layer_width; - layerIndex = layer_id; + double layer_rmin = chamber_radius_min+(layer_id*layer_width); + double layer_rmax = layer_rmin+layer_width; + double rmid_zZero = (layer_rmin+layer_rmax)/2.; // z=0 + double rmid_zEnd = rmid_zZero/std::cos(alpha/2); // z=endcap + int nCell = floor((2. * M_PI * rmid_zZero) / layer_width); + int nWire = nCell; + if(!DC_construct_wire) nWire =0; + double cell_phi = 2*M_PI / nCell; + double offset=0;//phi offset of first cell in each layer + double sign_eps = 1;// setero angle sign + if(0==(layer_id%2)) { + offset = 0.; + sign_eps = -1; + } else { + offset = 0.5 * cell_phi; + } + double epsilon = sign_eps*std::atan(rmid_zZero * std::tan(alpha / 2.0) / chamber_half_length); + + segmentationDC->setGeomParams(chamber_id, layer_id, cell_phi, rmid_zEnd , epsilon, offset); + segmentationDC->setWiresInLayer(chamber_id, layer_id, nCell); - //Construction of drift chamber layers - double rmid = delta_a_func(rmin,rmax); - double ilayer_cir = 2 * M_PI * rmid; - double ncell = ilayer_cir / chamber_layer_width; - int ncell_layer = floor(ncell); - int numWire = ncell_layer; - double layer_Phi = 2*M_PI / ncell_layer; + double r_in_test = rmid_zZero*std::cos(alpha / 2.0); - if(layer_id %2 ==0){ offset = 0.; } - else { offset = 0.5 * layer_Phi; } + double r_in0 = rmid_zZero - layer_width / 2.0;// + double r_in = r_in0 / std::cos(alpha / 2.0);//layer min at z=half length - DCHseg->setGeomParams(chamber_id, layerIndex, layer_Phi, rmid, epsilon, offset); - DCHseg->setWiresInLayer(chamber_id, layerIndex, numWire); + double r_out0 = rmid_zZero + layer_width / 2.0; + double r_out = r_out0 / std::cos(alpha / 2.0); + double delta_a_in = delta_b_func(r_in, r_in0); + double delta_a_out = delta_b_func(r_out, r_out0); - dd4hep::Tube layer_vol_solid(rmin,rmax,chamber_half_length); + double eps_in = epsilon_func(delta_a_in, chamber_length ); + double eps_out = epsilon_func(delta_a_out, chamber_length ); + + /// create hyper layer volume + dd4hep::Hyperboloid layer_vol_solid(r_in0, eps_in, r_out0, eps_out, chamber_half_length); dd4hep::Volume layer_vol(det_name+"_layer_vol",layer_vol_solid,det_mat); current_vol_ptr = &layer_vol; @@ -184,46 +210,56 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector, sd.setType("tracker"); } - // - wire vol - //phi <-------------------> -phi - // | F8 F7 F6 F5| Only on the outermost layer. - // | | - // | S F4| - // | | - // | F0 F1 F2 F3| - // ----------------------- - for(int icell=0; icell< numWire; icell++) { - double wire_phi = (icell+0.5)*layer_Phi + offset; + // - create wire volume + //phi <---------------> -phi + // | F4 F3| Only on the outermost layer. + // | | + // | S F2| + // | | + // | F0 F1| + //-------------------- + // loop over cells + for(int icell=0; icell elapsed_seconds = end-start; + std::cout << "elapsed time: " << elapsed_seconds.count() << "s\n"; return sdet; } diff --git a/Detector/DetSegmentation/include/DetSegmentation/GridDriftChamber.h b/Detector/DetSegmentation/include/DetSegmentation/GridDriftChamber.h index 0233ccf82..77e9bfa56 100644 --- a/Detector/DetSegmentation/include/DetSegmentation/GridDriftChamber.h +++ b/Detector/DetSegmentation/include/DetSegmentation/GridDriftChamber.h @@ -34,7 +34,6 @@ typedef struct CID { int chamberID; int layerID; -// CID(){} CID(int i, int j): chamberID(i),layerID(j){} // the operator < defines the operation used in map friend bool operator < (const CID &c1, const CID &c2); @@ -60,23 +59,25 @@ class GridDriftChamber : public Segmentation { virtual CellID cellID(const Vector3D& aLocalPosition, const Vector3D& aGlobalPosition, const VolumeID& aVolumeID) const; virtual double distanceTrackWire(const CellID& cID, const TVector3& hit_start, const TVector3& hit_end) const; + virtual double distanceTrackWire2(const CellID& cID, const TVector3& hit_pos) const; virtual void cellposition(const CellID& cID, TVector3& Wstart, TVector3& Wend) const; + virtual void cellposition2(int chamber, int layer, int cell, TVector3& Wstart, TVector3& Wend) const; TVector3 LineLineIntersect(TVector3& p1, TVector3& p2, TVector3& p3, TVector3& p4) const; - virtual TVector3 distanceClosestApproach(const CellID& cID, const TVector3& hitPos) const; + virtual TVector3 distanceClosestApproach(const CellID& cID, const TVector3& hitPos, TVector3& PCA) const; virtual TVector3 Line_TrackWire(const CellID& cID, const TVector3& hit_start, const TVector3& hit_end) const; virtual TVector3 IntersectionTrackWire(const CellID& cID, const TVector3& hit_start, const TVector3& hit_end) const; virtual TVector3 wirePos_vs_z(const CellID& cID, const double& zpos) const; + virtual double Distance(const CellID& cID, const TVector3& pointIn, const TVector3& pointOut, TVector3& hitPosition, TVector3& PCA) const; + virtual TVector3 returnPhi0(int chamber,int layer, double z) const; + // double phi(const CellID& cID) const; inline double cell_Size() const { return m_cellSize; } - inline double epsilon0() const { return m_epsilon0; } + inline double epsilon() const { return m_epsilon; } inline double detectorLength() const { return m_detectorLength; } - inline double safe_distance() const { return m_safe_distance; } inline double layer_width() const { return m_layer_width; } inline double DC_rbegin() const { return m_DC_rbegin; } inline double DC_rend() const { return m_DC_rend; } - inline double DC_rmax() const { return m_DC_rmax; } - inline double DC_rmin() const { return m_DC_rmin; } inline const std::string& fieldNamePhi() const { return m_phiID; } inline const std::string& Layerid() const { return layer_id; } @@ -88,7 +89,13 @@ class GridDriftChamber : public Segmentation { return hit_phi; } - inline void setGeomParams(int chamberID, int layerID, double layerphi, double R, double eps, double offset) { + inline double phiFromXY2(const TVector3& aposition) const { + double hit_phi = std::atan2(aposition.Y(), aposition.X()) ; + if( hit_phi < 0 ) { hit_phi += 2 * M_PI; } + return hit_phi; + } + + inline void setGeomParams(int chamberID, int layerID, double layerphi, double R, double eps, double offset) { layer_params.insert(std::pair(vID(chamberID,layerID),LAYER(layerphi,R,eps,offset))); @@ -135,7 +142,18 @@ class GridDriftChamber : public Segmentation { _currentRadius = Radius; m_epsilon = Eps; m_offset = Offset; - } + } + + inline Vector3D returnPosWire0(double z) const { + double alpha = returnAlpha(); + double t = 0.5 * (1 - 2.0 * z / m_detectorLength); + double x = _currentRadius * (1 + t * (std::cos(alpha) - 1)); + double y = _currentRadius * t * std::sin(alpha); + + Vector3D vec(x, y, z); + return vec; + } + protected: @@ -152,19 +170,15 @@ class GridDriftChamber : public Segmentation { } inline double returnAlpha() const { - double alpha = 2 * std::asin(m_detectorLength * std::tan(m_epsilon0)/(2 * _currentRadius)); + double alpha = 2 * std::asin(m_detectorLength * std::tan(m_epsilon)/(2 * _currentRadius)); return alpha; } double m_cellSize; - double m_epsilon0; double m_detectorLength; double m_layer_width; - double m_safe_distance; double m_DC_rbegin; double m_DC_rend; - double m_DC_rmax; - double m_DC_rmin; std::string m_phiID; std::string layer_id; diff --git a/Detector/DetSegmentation/src/GridDriftChamber.cpp b/Detector/DetSegmentation/src/GridDriftChamber.cpp index a4302ebf5..5c59f2006 100644 --- a/Detector/DetSegmentation/src/GridDriftChamber.cpp +++ b/Detector/DetSegmentation/src/GridDriftChamber.cpp @@ -22,16 +22,12 @@ GridDriftChamber::GridDriftChamber(const BitFieldCoder* decoder) : Segmentation( _description = "Drift chamber segmentation in the global coordinates"; registerParameter("cell_size", "cell size", m_cellSize, 1., SegmentationParameter::LengthUnit); - registerParameter("epsilon0", "epsilon", m_epsilon0, 0., SegmentationParameter::AngleUnit, true); registerParameter("detector_length", "Length of the wire", m_detectorLength, 1., SegmentationParameter::LengthUnit); registerIdentifier("identifier_phi", "Cell ID identifier for phi", m_phiID, "cellID"); registerIdentifier("layerID", "layer id", layer_id, "layer"); - registerParameter("safe_distance", "safe_distance", m_safe_distance, 0., SegmentationParameter::LengthUnit); registerParameter("layer_width", "layer_width", m_layer_width, 0., SegmentationParameter::LengthUnit); registerParameter("DC_rbegin", "DC_rbegin", m_DC_rbegin, 0., SegmentationParameter::LengthUnit); registerParameter("DC_rend", "DC_rend", m_DC_rend, 0., SegmentationParameter::LengthUnit); - registerParameter("DC_rmin", "DC_rmin", m_DC_rmin, 0., SegmentationParameter::LengthUnit); - registerParameter("DC_rmax", "DC_rmax", m_DC_rmax, 0., SegmentationParameter::LengthUnit); } Vector3D GridDriftChamber::position(const CellID& /*cID*/) const { @@ -45,214 +41,293 @@ CellID GridDriftChamber::cellID(const Vector3D& /*localPosition*/, const Vector3 CellID cID = vID; int chamberID = _decoder->get(cID, "chamber"); + int layerid = _decoder->get(cID, "layer"); double posx = globalPosition.X; double posy = globalPosition.Y; - double radius = sqrt(posx*posx+posy*posy); - - int m_DC_layer_number = floor((m_DC_rend-m_DC_rbegin)/m_layer_width); - double DC_layerdelta = m_layer_width; - - int layerid; - if( radius<= m_DC_rend && radius>= m_DC_rbegin) { - layerid = floor((radius - m_DC_rbegin)/DC_layerdelta); - } else if ( radius>= (m_DC_rmin-m_safe_distance) && radius < m_DC_rbegin) { - layerid = 0; - } else if ( radius> m_DC_rend && radius <= (m_DC_rmax+m_safe_distance)) { - layerid = m_DC_layer_number-1; - } + double posz = globalPosition.Z; updateParams(chamberID,layerid); + TVector3 Phi0 = returnPhi0(chamberID,layerid,posz); + double phi0 = phiFromXY2(Phi0); + double phi_hit = phiFromXY(globalPosition); double offsetphi= m_offset; - int _lphi; + double _lphi; - if(phi_hit >= offsetphi) { - _lphi = (int) ((phi_hit - offsetphi)/ _currentLayerphi); - } - else { - _lphi = (int) ((phi_hit - offsetphi + 2 * M_PI)/ _currentLayerphi); + _lphi = phi_hit - phi0 + _currentLayerphi/2.; + if(_lphi<0.){ + _lphi+=2*M_PI; + }else if(_lphi>2*M_PI){ + _lphi=fmod(_lphi,2*M_PI); } + int cellID=floor(_lphi/_currentLayerphi); - int lphi = _lphi; - _decoder->set(cID, layer_id, layerid); - _decoder->set(cID, m_phiID, lphi); + _decoder->set(cID, m_phiID, cellID); return cID; } double GridDriftChamber::phi(const CellID& cID) const { - CellID phiValue = _decoder->get(cID, m_phiID); - return binToPosition(phiValue, _currentLayerphi, m_offset); + CellID phiValue = _decoder->get(cID, m_phiID); + return binToPosition(phiValue, _currentLayerphi, m_offset); } void GridDriftChamber::cellposition(const CellID& cID, TVector3& Wstart, - TVector3& Wend) const { + TVector3& Wend) const { + + auto chamberIndex = _decoder->get(cID, "chamber"); + auto layerIndex = _decoder->get(cID, "layer"); + updateParams(chamberIndex,layerIndex); + + double phi_start = phi(cID); + double phi_end = phi_start + returnAlpha(); + + Wstart = returnWirePosition(phi_start, -1); + Wend = returnWirePosition(phi_end, 1); +} + +TVector3 GridDriftChamber::returnPhi0(int chamber,int layer, double z) const +{ + updateParams(chamber,layer); + + double phi_wire_start = binToPosition(0 , _currentLayerphi, m_offset); + double phi_wire_end = phi_wire_start + returnAlpha(); + + TVector3 wire_start = returnWirePosition(phi_wire_start, -1); + TVector3 wire_end = returnWirePosition(phi_wire_end, 1); + + double ratio = fabs(z - wire_start.Z())/fabs(wire_end.Z() - wire_start.Z()); + double x_pos = ratio * (wire_end.X() - wire_start.X()) + wire_start.X(); + double y_pos = ratio * (wire_end.Y() - wire_start.Y()) + wire_start.Y(); + + return TVector3(x_pos,y_pos,z); +} - auto chamberIndex = _decoder->get(cID, "chamber"); - auto layerIndex = _decoder->get(cID, "layer"); - updateParams(chamberIndex,layerIndex); - double phi_start = phi(cID); - double phi_mid = phi_start + _currentLayerphi/2.; - double phi_end = phi_mid + returnAlpha(); +void GridDriftChamber::cellposition2(int chamber,int layer, int cell, + TVector3& Wstart, TVector3& Wend) const { + updateParams(chamber,layer); + double phi_start = binToPosition(cell, _currentLayerphi, m_offset); + double phi_end = phi_start + returnAlpha(); - Wstart = returnWirePosition(phi_mid, -1); - Wend = returnWirePosition(phi_end, 1); + Wstart = returnWirePosition(phi_start, -1); + Wend = returnWirePosition(phi_end, 1); } TVector3 GridDriftChamber::LineLineIntersect(TVector3& p1, TVector3& p2, TVector3& p3, TVector3& p4) const { - TVector3 p13, p43, p21; - double d1343, d4321, d1321, d4343, d2121; - double numer, denom; - double mua, mub; - TVector3 pa, pb; - - p13.SetX(p1.X() - p3.X()); - p13.SetY(p1.Y() - p3.Y()); - p13.SetZ(p1.Z() - p3.Z()); - p43.SetX(p4.X() - p3.X()); - p43.SetY(p4.Y() - p3.Y()); - p43.SetZ(p4.Z() - p3.Z()); - /* if (ABS(p43.X()) < EPS && ABS(p43.Y()) < EPS && ABS(p43.Z()) < EPS) */ - /* return(FALSE); */ - p21.SetX(p2.X() - p1.X()); - p21.SetY(p2.Y() - p1.Y()); - p21.SetZ(p2.Z() - p1.Z()); - /* if (ABS(p21.X()) < EPS && ABS(p21.Y()) < EPS && ABS(p21.Z()) < EPS) */ - /* return(FALSE); */ - - d1343 = p13.X() * p43.X() + p13.Y() * p43.Y() + p13.Z() * p43.Z(); - d4321 = p43.X() * p21.X() + p43.Y() * p21.Y() + p43.Z() * p21.Z(); - d1321 = p13.X() * p21.X() + p13.Y() * p21.Y() + p13.Z() * p21.Z(); - d4343 = p43.X() * p43.X() + p43.Y() * p43.Y() + p43.Z() * p43.Z(); - d2121 = p21.X() * p21.X() + p21.Y() * p21.Y() + p21.Z() * p21.Z(); - - denom = d2121 * d4343 - d4321 * d4321; - /* if (ABS(denom) < EPS) */ - /* return(FALSE); */ - numer = d1343 * d4321 - d1321 * d4343; - - mua = numer / denom; - mub = (d1343 + d4321 * (mua)) / d4343; - - pa.SetX(p1.X() + mua * p21.X()); - pa.SetY(p1.Y() + mua * p21.Y()); - pa.SetZ(p1.Z() + mua * p21.Z()); - pb.SetX(p3.X() + mub * p43.X()); - pb.SetY(p3.Y() + mub * p43.Y()); - pb.SetZ(p3.Z() + mub * p43.Z()); - - return pb - pa; + TVector3 p13, p43, p21; + double d1343, d4321, d1321, d4343, d2121; + double numer, denom; + double mua, mub; + TVector3 pa, pb; + + p13.SetX(p1.X() - p3.X()); + p13.SetY(p1.Y() - p3.Y()); + p13.SetZ(p1.Z() - p3.Z()); + p43.SetX(p4.X() - p3.X()); + p43.SetY(p4.Y() - p3.Y()); + p43.SetZ(p4.Z() - p3.Z()); + /* if (ABS(p43.X()) < EPS && ABS(p43.Y()) < EPS && ABS(p43.Z()) < EPS) */ + /* return(FALSE); */ + p21.SetX(p2.X() - p1.X()); + p21.SetY(p2.Y() - p1.Y()); + p21.SetZ(p2.Z() - p1.Z()); + /* if (ABS(p21.X()) < EPS && ABS(p21.Y()) < EPS && ABS(p21.Z()) < EPS) */ + /* return(FALSE); */ + + d1343 = p13.X() * p43.X() + p13.Y() * p43.Y() + p13.Z() * p43.Z(); + d4321 = p43.X() * p21.X() + p43.Y() * p21.Y() + p43.Z() * p21.Z(); + d1321 = p13.X() * p21.X() + p13.Y() * p21.Y() + p13.Z() * p21.Z(); + d4343 = p43.X() * p43.X() + p43.Y() * p43.Y() + p43.Z() * p43.Z(); + d2121 = p21.X() * p21.X() + p21.Y() * p21.Y() + p21.Z() * p21.Z(); + + denom = d2121 * d4343 - d4321 * d4321; + /* if (ABS(denom) < EPS) */ + /* return(FALSE); */ + numer = d1343 * d4321 - d1321 * d4343; + + mua = numer / denom; + mub = (d1343 + d4321 * (mua)) / d4343; + + pa.SetX(p1.X() + mua * p21.X()); + pa.SetY(p1.Y() + mua * p21.Y()); + pa.SetZ(p1.Z() + mua * p21.Z()); + pb.SetX(p3.X() + mub * p43.X()); + pb.SetY(p3.Y() + mub * p43.Y()); + pb.SetZ(p3.Z() + mub * p43.Z()); + + return pb - pa; } double GridDriftChamber::distanceTrackWire(const CellID& cID, const TVector3& hit_start, - const TVector3& hit_end) const { + const TVector3& hit_end) const { - TVector3 Wstart = {0,0,0}; - TVector3 Wend = {0,0,0}; - cellposition(cID,Wstart,Wend); + TVector3 Wstart = {0,0,0}; + TVector3 Wend = {0,0,0}; + cellposition(cID,Wstart,Wend); - TVector3 a = (hit_end - hit_start).Unit(); - TVector3 b = (Wend - Wstart).Unit(); - TVector3 c = Wstart - hit_start; + TVector3 a = (hit_end - hit_start).Unit(); + TVector3 b = (Wend - Wstart).Unit(); + TVector3 c = Wstart - hit_start; - double num = std::abs(c.Dot(a.Cross(b))); - double denum = (a.Cross(b)).Mag(); + double num = std::abs(c.Dot(a.Cross(b))); + double denum = (a.Cross(b)).Mag(); - double DCA = 0; + double DCA = 0; - if (denum) { - DCA = num / denum; - } + if (denum) { + DCA = num / denum; + } - return DCA; + return DCA; } -TVector3 GridDriftChamber::distanceClosestApproach(const CellID& cID, const TVector3& hitPos) const { - // Distance of the closest approach between a single hit (point) and the closest wire +double GridDriftChamber::distanceTrackWire2(const CellID& cID, const TVector3& hit_pos) const { - TVector3 Wstart = {0,0,0}; - TVector3 Wend = {0,0,0}; - cellposition(cID,Wstart,Wend); + TVector3 Wstart = {0,0,0}; + TVector3 Wend = {0,0,0}; + cellposition(cID,Wstart,Wend); - TVector3 temp = (Wend + Wstart); - TVector3 Wmid(temp.X() / 2.0, temp.Y() / 2.0, temp.Z() / 2.0); + TVector3 denominator = Wend - Wstart; + TVector3 numerator = denominator.Cross(Wstart-hit_pos); - double hitPhi = hitPos.Phi(); - if (hitPhi < 0) { - hitPhi = hitPhi + 2 * M_PI; - } + double DCA = numerator.Mag()/denominator.Mag() ; + + return DCA; +} - TVector3 PCA = Wstart + ((Wend - Wstart).Unit()).Dot((hitPos - Wstart)) * ((Wend - Wstart).Unit()); - TVector3 dca = hitPos - PCA; +TVector3 GridDriftChamber::distanceClosestApproach(const CellID& cID, const TVector3& hitPos, TVector3& PCA) const { + // Distance of the closest approach between a single hit (point) and the closest wire - return dca; + TVector3 Wstart = {0,0,0}; + TVector3 Wend = {0,0,0}; + cellposition(cID,Wstart,Wend); + + TVector3 temp = (Wend + Wstart); + TVector3 Wmid(temp.X() / 2.0, temp.Y() / 2.0, temp.Z() / 2.0); + + double hitPhi = hitPos.Phi(); + if (hitPhi < 0) { + hitPhi = hitPhi + 2 * M_PI; + } + + PCA = Wstart + ((Wend - Wstart).Unit()).Dot((hitPos - Wstart)) * ((Wend - Wstart).Unit()); + TVector3 dca = hitPos - PCA; + + return dca; } TVector3 GridDriftChamber::Line_TrackWire(const CellID& cID, const TVector3& hit_start, const TVector3& hit_end) const { - // The line connecting a particle track to the closest wire - // Returns the vector connecting the both - TVector3 Wstart = {0,0,0}; - TVector3 Wend = {0,0,0}; - cellposition(cID,Wstart,Wend); - - TVector3 P1 = hit_start; - TVector3 P2 = hit_end; - TVector3 P3 = Wstart; - TVector3 P4 = Wend; - - TVector3 intersect = LineLineIntersect(P1, P2, P3, P4); - return intersect; + // The line connecting a particle track to the closest wire + // Returns the vector connecting the both + TVector3 Wstart = {0,0,0}; + TVector3 Wend = {0,0,0}; + cellposition(cID,Wstart,Wend); + + TVector3 P1 = hit_start; + TVector3 P2 = hit_end; + TVector3 P3 = Wstart; + TVector3 P4 = Wend; + + TVector3 intersect = LineLineIntersect(P1, P2, P3, P4); + return intersect; } // Get the wire position for a z TVector3 GridDriftChamber::wirePos_vs_z(const CellID& cID, const double& zpos) const { - TVector3 Wstart = {0,0,0}; - TVector3 Wend = {0,0,0}; - cellposition(cID,Wstart,Wend); + TVector3 Wstart = {0,0,0}; + TVector3 Wend = {0,0,0}; + cellposition(cID,Wstart,Wend); - double t = (zpos - Wstart.Z())/(Wend.Z()-Wstart.Z()); - double x = Wstart.X()+t*(Wend.X()-Wstart.X()); - double y = Wstart.Y()+t*(Wend.Y()-Wstart.Y()); + double t = (zpos - Wstart.Z())/(Wend.Z()-Wstart.Z()); + double x = Wstart.X()+t*(Wend.X()-Wstart.X()); + double y = Wstart.Y()+t*(Wend.Y()-Wstart.Y()); - TVector3 wireCoord(x, y, zpos); - return wireCoord; + TVector3 wireCoord(x, y, zpos); + return wireCoord; } TVector3 GridDriftChamber::IntersectionTrackWire(const CellID& cID, const TVector3& hit_start, const TVector3& hit_end) const { - // Intersection between the particle track and the wire assuming that the track between hit_start and hit_end is linear + // Intersection between the particle track and the wire assuming that the track between hit_start and hit_end is linear - TVector3 Wstart = {0,0,0}; - TVector3 Wend = {0,0,0}; - cellposition(cID,Wstart,Wend); + TVector3 Wstart = {0,0,0}; + TVector3 Wend = {0,0,0}; + cellposition(cID,Wstart,Wend); - TVector3 P1 = hit_start; - TVector3 V1 = hit_end-hit_start; + TVector3 P1 = hit_start; + TVector3 V1 = hit_end-hit_start; - TVector3 P2 = Wstart; - TVector3 V2 = Wend - Wstart; + TVector3 P2 = Wstart; + TVector3 V2 = Wend - Wstart; - TVector3 denom = V1.Cross(V2); - double mag_denom = denom.Mag(); + TVector3 denom = V1.Cross(V2); + double mag_denom = denom.Mag(); - TVector3 intersect(0, 0, 0); + TVector3 intersect(0, 0, 0); - if (mag_denom !=0) + if (mag_denom !=0) { - TVector3 num = ((P2-P1)).Cross(V2); - double mag_num = num.Mag(); - double a = mag_num / mag_denom; + TVector3 num = ((P2-P1)).Cross(V2); + double mag_num = num.Mag(); + double a = mag_num / mag_denom; - intersect = P1 + a * V1; + intersect = P1 + a * V1; } - return intersect; + return intersect; } +double GridDriftChamber::Distance(const CellID& cID, const TVector3& pointIn, const TVector3& pointOut, TVector3& hitPosition, TVector3& PCA) const { + + //For two lines r=r1+t1.v1 & r=r2+t2.v2 + //the closest approach is d=|(r2-r1).(v1 x v2)|/|v1 x v2| + //the point where closest approach are + //t1=(v1 x v2).[(r2-r1) x v2]/[(v1 x v2).(v1 x v2)] + //t2=(v1 x v2).[(r2-r1) x v1]/[(v1 x v2).(v1 x v2)] + //if v1 x v2=0 means two lines are parallel + //d=|(r2-r1) x v1|/|v1| + + double t1,distance,dInOut,dHitIn,dHitOut; + //Get wirepoint @ endplate + TVector3 west = {0,0,0}; + TVector3 east = {0,0,0}; + cellposition(cID,west,east); + TVector3 wireLine=east - west; + TVector3 hitLine=pointOut - pointIn; + + TVector3 hitXwire=hitLine.Cross(wireLine); + TVector3 wire2hit=east-pointOut; + //Hitposition is the position on hit line where closest approach + //of two lines, but it may out the area from pointIn to pointOut + if(hitXwire.Mag()==0){ + distance=wireLine.Cross(wire2hit).Mag()/wireLine.Mag(); + hitPosition=pointIn; + }else{ + t1=hitXwire.Dot(wire2hit.Cross(wireLine))/hitXwire.Mag2(); + hitPosition=pointOut+t1*hitLine; + + dInOut=(pointOut-pointIn).Mag(); + dHitIn=(hitPosition-pointIn).Mag(); + dHitOut=(hitPosition-pointOut).Mag(); + if(dHitIn<=dInOut && dHitOut<=dInOut){ //Between point in & out + distance=fabs(wire2hit.Dot(hitXwire)/hitXwire.Mag()); + }else if(dHitOut>dHitIn){ // out pointIn + distance=wireLine.Cross(pointIn-east).Mag()/wireLine.Mag(); + hitPosition=pointIn; + }else{ // out pointOut + distance=wireLine.Cross(pointOut-east).Mag()/wireLine.Mag(); + hitPosition=pointOut; + } + } + + PCA = west + ((east - west).Unit()).Dot((hitPosition - west)) * ((east - west).Unit()); + + return distance; +} } } diff --git a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp index 8f5a8cefa..21a840411 100644 --- a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp +++ b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp @@ -28,15 +28,15 @@ DCHDigiAlg::DCHDigiAlg(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc), _nEvt(0) { - + // Input collections declareProperty("SimDCHitCollection", r_SimDCHCol, "Handle of the Input SimHit collection"); - + // Output collections declareProperty("DigiDCHitCollection", w_DigiDCHCol, "Handle of Digi DCHit collection"); - + declareProperty("AssociationCollection", w_AssociationCol, "Handle of Association collection"); - + } StatusCode DCHDigiAlg::initialize() @@ -53,7 +53,10 @@ StatusCode DCHDigiAlg::initialize() return StatusCode::FAILURE; } + fRandom.SetSeed(105105);//FIXME: set by users + if(m_WriteAna){ + NTuplePtr nt( ntupleSvc(), "MyTuples/DCH_digi_evt" ); if ( nt ) m_tuple = nt; else { @@ -66,17 +69,28 @@ StatusCode DCHDigiAlg::initialize() m_tuple->addItem( "simhit_x", m_n_sim, m_simhit_x).ignore(); m_tuple->addItem( "simhit_y", m_n_sim, m_simhit_y).ignore(); m_tuple->addItem( "simhit_z", m_n_sim, m_simhit_z).ignore(); + m_tuple->addItem( "Simdca" , m_n_sim, m_Simdca ).ignore(); + m_tuple->addItem( "simhitT", m_n_sim, m_simhitT ).ignore(); + m_tuple->addItem( "simhitmom", m_n_sim, m_simhitmom ).ignore(); + m_tuple->addItem( "simPDG", m_n_sim, m_simPDG ).ignore(); m_tuple->addItem( "chamber" , m_n_digi, m_chamber ).ignore(); m_tuple->addItem( "layer" , m_n_digi, m_layer ).ignore(); m_tuple->addItem( "cell" , m_n_digi, m_cell ).ignore(); m_tuple->addItem( "cell_x" , m_n_digi, m_cell_x ).ignore(); m_tuple->addItem( "cell_y" , m_n_digi, m_cell_y ).ignore(); + m_tuple->addItem( "cell1_x" , m_n_digi, m_cell1_x ).ignore(); + m_tuple->addItem( "cell1_y" , m_n_digi, m_cell1_y ).ignore(); m_tuple->addItem( "hit_x" , m_n_digi,m_hit_x ).ignore(); m_tuple->addItem( "hit_y" , m_n_digi,m_hit_y ).ignore(); m_tuple->addItem( "hit_z" , m_n_digi,m_hit_z ).ignore(); - m_tuple->addItem( "dca" , m_n_digi,m_dca ).ignore(); + m_tuple->addItem( "mom_x" , m_n_digi,m_mom_x ).ignore(); + m_tuple->addItem( "mom_y" , m_n_digi,m_mom_y ).ignore(); + m_tuple->addItem( "dca" , m_n_digi, m_dca ).ignore(); + m_tuple->addItem( "poca_x" , m_n_digi, m_poca_x ).ignore(); + m_tuple->addItem( "poca_y" , m_n_digi, m_poca_y ).ignore(); m_tuple->addItem( "hit_dE" , m_n_digi,m_hit_dE ).ignore(); m_tuple->addItem( "hit_dE_dx", m_n_digi,m_hit_dE_dx ).ignore(); + m_tuple->addItem( "truthlength", m_n_digi,m_truthlength).ignore(); } else { // did not manage to book the N tuple.... info() << " Cannot book N-tuple:" << long( m_tuple ) << endmsg; } @@ -92,7 +106,7 @@ StatusCode DCHDigiAlg::execute() m_start = clock(); info() << "Processing " << _nEvt << " events " << endmsg; - m_evt = _nEvt; + if(m_WriteAna) m_evt = _nEvt; edm4hep::TrackerHitCollection* Vec = w_DigiDCHCol.createAndPut(); edm4hep::MCRecoTrackerAssociationCollection* AssoVec = w_AssociationCol.createAndPut(); const edm4hep::SimTrackerHitCollection* SimHitCol = r_SimDCHCol.get(); @@ -102,8 +116,8 @@ StatusCode DCHDigiAlg::execute() debug()<<"input sim hit size="<< SimHitCol->size() <at(0); - std::map > id_hits_map; - //std::map > id_hits_map; + std::map> id_hits_map; + for( int i = 0; i < SimHitCol->size(); i++ ) { @@ -111,7 +125,12 @@ StatusCode DCHDigiAlg::execute() unsigned long long id = SimHit.getCellID(); float sim_hit_mom = sqrt( SimHit.getMomentum()[0]*SimHit.getMomentum()[0] + SimHit.getMomentum()[1]*SimHit.getMomentum()[1] + SimHit.getMomentum()[2]*SimHit.getMomentum()[2] );//GeV if(sim_hit_mom < m_mom_threshold) continue; - if(SimHit.getEDep() <= 0) continue; + if(sim_hit_mom > m_mom_threshold_high) continue; + if(SimHit.getEDep() <= m_edep_threshold) continue; + + //Wire efficiency + double hitProb = fRandom.Uniform(1.); + if(hitProb > m_wireEff) continue; if ( id_hits_map.find(id) != id_hits_map.end()) id_hits_map[id].push_back(SimHit); else @@ -127,120 +146,128 @@ StatusCode DCHDigiAlg::execute() } for(auto iter = id_hits_map.begin(); iter != id_hits_map.end(); iter++) { - unsigned long long wcellid = iter->first; - auto trkHit = Vec->create(); - trkHit.setCellID(wcellid); - double tot_edep = 0 ; - double tot_length = 0 ; - double pos_x = 0 ; - double pos_y = 0 ; - double pos_z = 0 ; - int simhit_size = iter->second.size(); - for(unsigned int i=0; i< simhit_size; i++) - { - tot_edep += iter->second.at(i).getEDep();//GeV - } - int chamber = m_decoder->get(wcellid, "chamber"); - int layer = m_decoder->get(wcellid, "layer" ); - int cellID = m_decoder->get(wcellid, "cellID" ); - TVector3 Wstart(0,0,0); - TVector3 Wend (0,0,0); - m_segmentation->cellposition(wcellid, Wstart, Wend); - float dd4hep_mm = dd4hep::mm; - //std::cout<<"dd4hep_mm="<get(wcellid, "layer" ); + int cellID = m_decoder->get(wcellid, "cellID" ); + TVector3 Wstart(0,0,0); + TVector3 Wend (0,0,0); + m_segmentation->cellposition(wcellid, Wstart, Wend); + float dd4hep_mm = dd4hep::mm; + if(m_debug) std::cout<<"DCHDigi wcellid ="<second.at(i)<second.at(i).getTime(); + m_simhitmom[m_n_sim] = sim_hit_mom; + m_simPDG[m_n_sim] = iter->second.at(i).getMCParticle().getPDG(); + m_n_sim ++ ; + } + } + trkHit.setTime(min_distance*1e3/m_velocity);//m_velocity is um/ns, drift time in ns + trkHit.setEDep(tot_edep);// GeV + //trkHit.setEdx (tot_edep/tot_length); // GeV/mm + trkHit.setPosition (edm4hep::Vector3d(pos_x, pos_y, pos_z));//position of closest sim hit + trkHit.setCovMatrix(std::array{m_res_x, 0, m_res_y, 0, 0, m_res_z});//cov(x,x) , cov(y,x) , cov(y,y) , cov(z,x) , cov(z,y) , cov(z,z) in mm + + if(m_WriteAna && (nullptr!=m_tuple)) { + m_chamber [m_n_digi] = chamber; + m_layer [m_n_digi] = layer ; + m_cell [m_n_digi] = cellID; + m_cell_x [m_n_digi] = Wstart.x(); + m_cell_y [m_n_digi] = Wstart.y(); + m_cell1_x [m_n_digi] = Wend.x(); + m_cell1_y [m_n_digi] = Wend.y(); + m_hit_x [m_n_digi] = pos_x; + m_hit_y [m_n_digi] = pos_y; + m_hit_z [m_n_digi] = pos_z; + m_mom_x [m_n_digi] = momx ; + m_mom_y [m_n_digi] = momy ; + m_dca [m_n_digi] = min_distance; + m_poca_x [m_n_digi] = PCA.x(); + m_poca_y [m_n_digi] = PCA.y(); + m_hit_dE [m_n_digi] = trkHit.getEDep(); + //m_hit_dE_dx[m_n_digi] = trkHit.getEdx() ; + m_truthlength[m_n_digi] = tot_length ; + m_n_digi ++ ; + } + } debug()<<"output digi DCHhit size="<< Vec->size() <write(); if ( status.isFailure() ) { - error() << " Cannot fill N-tuple:" << long( m_tuple ) << endmsg; - return StatusCode::FAILURE; + error() << " Cannot fill N-tuple:" << long( m_tuple ) << endmsg; + return StatusCode::FAILURE; } } m_end = clock(); - m_time = (m_end - m_start); + if(m_WriteAna){ + m_time = (m_end - m_start); + } return StatusCode::SUCCESS; } StatusCode DCHDigiAlg::finalize() { - info() << "Processed " << _nEvt << " events " << endmsg; - return GaudiAlgorithm::finalize(); + info() << "Processed " << _nEvt << " events " << endmsg; + return GaudiAlgorithm::finalize(); } diff --git a/Digitisers/DCHDigi/src/DCHDigiAlg.h b/Digitisers/DCHDigi/src/DCHDigiAlg.h index d6a82961f..946363c1c 100644 --- a/Digitisers/DCHDigi/src/DCHDigiAlg.h +++ b/Digitisers/DCHDigi/src/DCHDigiAlg.h @@ -7,6 +7,7 @@ #include "edm4hep/SimTrackerHitCollection.h" #include "edm4hep/TrackerHitCollection.h" #include "edm4hep/MCRecoTrackerAssociationCollection.h" +#include "edm4hep/MCParticle.h" #include #include @@ -15,8 +16,7 @@ #include "DetSegmentation/GridDriftChamber.h" #include "TVector3.h" - - +#include "TRandom3.h" class DCHDigiAlg : public GaudiAlgorithm { @@ -44,6 +44,8 @@ class DCHDigiAlg : public GaudiAlgorithm typedef std::vector FloatVec; int _nEvt ; + TRandom3 fRandom; + NTuple::Tuple* m_tuple = nullptr ; NTuple::Item m_evt; NTuple::Item m_n_sim; @@ -52,40 +54,56 @@ class DCHDigiAlg : public GaudiAlgorithm NTuple::Array m_chamber ; NTuple::Array m_layer ; NTuple::Array m_cell ; + //The position of wire at -z NTuple::Array m_cell_x ; NTuple::Array m_cell_y ; + //The position of wire at +z + NTuple::Array m_cell1_x ; + NTuple::Array m_cell1_y ; + NTuple::Array m_simhit_x ; NTuple::Array m_simhit_y ; NTuple::Array m_simhit_z ; + NTuple::Array m_simhitT ; + NTuple::Array m_simhitmom ; + NTuple::Array m_simPDG ; NTuple::Array m_hit_x ; NTuple::Array m_hit_y ; NTuple::Array m_hit_z ; + NTuple::Array m_mom_x ; + NTuple::Array m_mom_y ; NTuple::Array m_dca ; + NTuple::Array m_Simdca ; + NTuple::Array m_poca_x ; + NTuple::Array m_poca_y ; NTuple::Array m_hit_dE ; NTuple::Array m_hit_dE_dx ; + NTuple::Array m_truthlength ; clock_t m_start,m_end; dd4hep::rec::CellIDPositionConverter* m_cellIDConverter; dd4hep::DDSegmentation::GridDriftChamber* m_segmentation; dd4hep::DDSegmentation::BitFieldCoder* m_decoder; - + Gaudi::Property m_readout_name{ this, "readout", "DriftChamberHitsCollection"};//readout for getting segmentation - + Gaudi::Property m_res_x { this, "res_x", 0.11};//mm Gaudi::Property m_res_y { this, "res_y", 0.11};//mm Gaudi::Property m_res_z { this, "res_z", 1 };//mm Gaudi::Property m_velocity { this, "drift_velocity", 40};// um/ns Gaudi::Property m_mom_threshold { this, "mom_threshold", 0};// GeV + Gaudi::Property m_mom_threshold_high { this, "mom_threshold_high", 1e9};// GeV + Gaudi::Property m_edep_threshold{ this, "edep_threshold", 0};// GeV Gaudi::Property m_WriteAna { this, "WriteAna", false}; - Gaudi::Property m_Doca { this, "Doca", false};//1:line dca 0:point dca - + Gaudi::Property m_debug{ this, "debug", false}; + Gaudi::Property m_wireEff{ this, "wireEff", 1.0}; // Input collections DataHandle r_SimDCHCol{"DriftChamberHitsCollection", Gaudi::DataHandle::Reader, this}; // Output collections DataHandle w_DigiDCHCol{"DigiDCHitCollection", Gaudi::DataHandle::Writer, this}; DataHandle w_AssociationCol{"DCHitAssociationCollection", Gaudi::DataHandle::Writer, this}; - }; + #endif diff --git a/Reconstruction/PFA/Arbor/src/MarlinArbor.cc b/Reconstruction/PFA/Arbor/src/MarlinArbor.cc old mode 100755 new mode 100644 diff --git a/Reconstruction/RecGenfitAlg/CMakeLists.txt b/Reconstruction/RecGenfitAlg/CMakeLists.txt index ee9b04300..f7c5b6961 100644 --- a/Reconstruction/RecGenfitAlg/CMakeLists.txt +++ b/Reconstruction/RecGenfitAlg/CMakeLists.txt @@ -2,11 +2,15 @@ if (GenFit_FOUND) gaudi_add_module(RecGenfitAlg - SOURCES src/RecGenfitAlgDC.cpp + SOURCES src/RecGenfitAlgSDT.cpp src/GenfitTrack.cpp + src/GenfitHit.cpp src/GenfitField.cpp src/GenfitFitter.cpp src/GenfitMaterialInterface.cpp + src/LSFitting.cpp + src/WireMeasurementDC.cpp + src/PlanarMeasurementSDT.cpp LINK GearSvc Gaudi::GaudiAlgLib Gaudi::GaudiKernel diff --git a/Reconstruction/RecGenfitAlg/src/GenfitField.cpp b/Reconstruction/RecGenfitAlg/src/GenfitField.cpp index 744084254..e0b8e3548 100644 --- a/Reconstruction/RecGenfitAlg/src/GenfitField.cpp +++ b/Reconstruction/RecGenfitAlg/src/GenfitField.cpp @@ -1,4 +1,5 @@ #include "GenfitField.h" +#include "GenfitUnit.h" //External #include "DD4hep/DD4hepUnits.h" @@ -39,16 +40,21 @@ GenfitField::get(const double& posX, const double& posY, const double& posZ, { /// get field from dd4hepField const dd4hep::Direction& field=m_dd4hepField.magneticField( - dd4hep::Position(posX/dd4hep::cm,posY/dd4hep::cm,posZ/dd4hep::cm)); - //m_dd4hepField.magneticField(pos,B); - - Bx=field.X()/dd4hep::kilogauss; - By=field.Y()/dd4hep::kilogauss; - Bz=field.Z()/dd4hep::kilogauss; - -// std::cout<<"GenfitField " -// < #include +#include +//#define GENFIT_MY_DEBUG 1 GenfitFitter::~GenfitFitter(){ delete m_absKalman; } -GenfitFitter::GenfitFitter(const char* type, const char* name): +GenfitFitter::GenfitFitter(const char* type,int debug,const char* name): m_absKalman(nullptr) ,m_genfitField(nullptr) ,m_geoMaterial(nullptr) @@ -49,7 +58,7 @@ GenfitFitter::GenfitFitter(const char* type, const char* name): ,m_maxIterations(10) ,m_deltaPval(1e-3) ,m_relChi2Change(0.2) - ,m_blowUpFactor(1e3) + ,m_blowUpFactor(500) ,m_resetOffDiagonals(true) ,m_blowUpMaxVal(1.e6) ,m_multipleMeasurementHandling(genfit::unweightedClosestToPredictionWire) @@ -66,9 +75,10 @@ GenfitFitter::GenfitFitter(const char* type, const char* name): ,m_noiseBrems(false) ,m_ignoreBoundariesBetweenEqualMaterials(true) ,m_mscModelName("GEANE") - ,m_debug(0) + //,m_debug(0) ,m_hist(0) { + m_debug=debug; /// Initialize genfit fitter init(); } @@ -78,15 +88,20 @@ void GenfitFitter::setField(const GenfitField* field) if(nullptr==m_genfitField) m_genfitField=field; } + /// Set geometry for material, use geometry from IOADatabase -void GenfitFitter::setGeoMaterial(const dd4hep::Detector* dd4hepGeo) +void GenfitFitter::setGeoMaterial(const dd4hep::Detector* dd4hepGeo, + double extDistCut, bool skipWireMaterial) { if(nullptr==m_geoMaterial){ m_geoMaterial=GenfitMaterialInterface::getInstance(dd4hepGeo); } + m_geoMaterial->setMinSafetyDistanceCut(extDistCut); + m_geoMaterial->setSkipWireMaterial(skipWireMaterial); + m_geoMaterial->setDebugLvl(m_debug); } -/// initialize genfit fitter +/// initialize genfit fitter, old fitter will be deleted int GenfitFitter::init(bool deleteOldFitter) { if(deleteOldFitter && m_absKalman) delete m_absKalman; @@ -95,17 +110,21 @@ int GenfitFitter::init(bool deleteOldFitter) <=2)std::cout<<" m_fitterType==DAFRef "<=2)std::cout<<" m_fitterType==DAF"<=2)std::cout<<" m_fitterType==KalmanFitter"<=2)std::cout<<" m_fitterType==KalmanFitterRefTrack"<=2)std::cout<<"Fitter type is "<setDebugLvl(m_debug); +#ifdef GENFIT_MY_DEBUG + //m_absKalman->setDebugLvlLocal(m_debugLocal); +#endif return 0; } @@ -125,6 +147,10 @@ int GenfitFitter::processTrackWithRep(GenfitTrack* track,int repID,bool resort) { if(m_debug>=2)std::cout<< "In ProcessTrackWithRep rep "<2) print(""); + if(track->getNumPoints()<=0){ + if(m_debug>=2)std::cout<<"skip track w.o. hit"<0){ if(m_debug>=2)std::cout<<"Print track seed "<=2)std::cout<<"In ProcessTrack"<getNumPoints()<=0){ + if(m_debug>=2)std::cout<<"skip track w.o. hit"<2) print(""); /// Do the fitting @@ -170,142 +200,6 @@ void GenfitFitter::setFitterType(const char* val) init(oldFitterType==val); } -/// Extrapolate track to the cloest point of approach(POCA) to the wire of hit, -/// return StateOnPlane of this POCA -/// inputs -/// pos,mom ... position & momentum at starting point, unit= [mm]&[GeV/c] -/// (converted to [cm]&[GeV/c] inside this function) -/// hit ... destination -/// outputs poca [mm] (converted from [cm] in this function) ,global -double GenfitFitter::extrapolateToHit( TVector3& poca, TVector3& pocaDir, - TVector3& pocaOnWire, double& doca, const GenfitTrack* track, - TVector3 pos, TVector3 mom, - TVector3 end0,//one end of the hit wire - TVector3 end1,//the orhter end of the hit wire - int debug, - int repID, - bool stopAtBoundary, - bool calcJacobianNoise) -{ - - return track->extrapolateToHit(poca,pocaDir,pocaOnWire,doca,pos,mom,end0, - end1,debug,repID,stopAtBoundary,calcJacobianNoise); -}//end of ExtrapolateToHit - - -/// Extrapolate the track to the cyliner at fixed raidus -/// position & momentum as starting point -/// position and momentum at global coordinate in dd4hepUnit -/// return trackLength in dd4hepUnit - double -GenfitFitter::extrapolateToCylinder(TVector3& pos, TVector3& mom, - GenfitTrack* track, double radius, const TVector3 linePoint, - const TVector3 lineDirection, int hitID, int repID, - bool stopAtBoundary, bool calcJacobianNoise) -{ - double trackLength(1e9*dd4hep::cm); - if(!track->getFitStatus(repID)->isFitted()) return trackLength; - try{ - // get track rep - genfit::AbsTrackRep* rep = track->getRep(repID); - if(nullptr == rep) { - if(m_debug>=2)std::cout<<"In ExtrapolateToCylinder rep is null" - <getTrack()->getPointWithFitterInfo(hitID,rep); - if(nullptr == tp) { - if(m_debug>=2)std::cout<< - "In ExtrapolateToCylinder TrackPoint is null"<( - tp->getFitterInfo(rep))->getBackwardUpdate(); - - if(nullptr == state){ - if(m_debug>=2)std::cout<<"In ExtrapolateToCylinder " - <<"no KalmanFittedStateOnPlane in backwardUpdate"<setPosMom(*state, pos*(1/dd4hep::cm), mom*(1/dd4hep::GeV)); - - /// extrapolate - trackLength = rep->extrapolateToCylinder(*state, - radius/dd4hep::cm, linePoint*(1/dd4hep::cm), lineDirection, - stopAtBoundary, calcJacobianNoise); - // get pos&mom at extrapolated point on the cylinder - rep->getPosMom(*state,pos,mom);//FIXME exception exist - pos = pos*dd4hep::cm; - mom = mom*dd4hep::GeV; - } catch(genfit::Exception& e){ - if(m_debug>=3)std::cout - <<"Exception in GenfitFitter::extrapolateToCylinder " - << e.what()<getFitStatus(repID)->isFitted()) return trackLength; - try{ - // get track rep - genfit::AbsTrackRep* rep = track->getRep(repID); - if(nullptr == rep) { - if(m_debug>=2)std::cout<<"In ExtrapolateToPoint rep " - <getTrack()->getFittedState(0,rep)))); - - // get track point with fitter info - genfit::TrackPoint* tp = - track->getTrack()->getPointWithFitterInfo(0,rep); - if(nullptr == tp) { - if(m_debug>=2)std::cout<< - "In ExtrapolateToPoint TrackPoint is null"<( - tp->getFitterInfo(rep))->getBackwardUpdate(); - - if(nullptr == state) { - if(m_debug>=2)std::cout<< - "In ExtrapolateToPoint KalmanFittedStateOnPlane is null"<extrapolateToPoint(*state, - point*(1/dd4hep::cm),stopAtBoundary, calcJacobianNoise); - rep->getPosMom(*state,pos,mom);//FIXME exception exist - pos = pos*dd4hep::cm; - mom = mom*dd4hep::GeV; - } catch(genfit::Exception& e){ - if(m_debug>=3)std::cout - <<"Exception in GenfitFitter::extrapolateToPoint" - << e.what()<getKalman()->setBlowUpFactor(m_blowUpFactor); + } } void GenfitFitter::setResetOffDiagonals(bool val) { m_absKalman->setResetOffDiagonals(val); m_resetOffDiagonals = val; + if (m_fitterType=="DAFRef" || m_fitterType=="DAF") { + getDAF()->getKalman()->setResetOffDiagonals(m_resetOffDiagonals); + } } void GenfitFitter::setBlowUpMaxVal(double val) { m_absKalman->setBlowUpMaxVal(val); m_blowUpMaxVal = val; + if (m_fitterType=="DAFRef" || m_fitterType=="DAF") { + getDAF()->getKalman()->setBlowUpMaxVal(m_blowUpMaxVal); + } } void GenfitFitter::setMultipleMeasurementHandling( @@ -558,16 +461,33 @@ void GenfitFitter::setMaterialDebugLvl(unsigned int val) genfit::MaterialEffects::getInstance()->setDebugLvl(val); } -///Set GenfitFitter parameters void GenfitFitter::setDebug(unsigned int val) { - if(m_debug>=2)std::cout<<"set fitter debugLvl "<setDebugLvl(val); - if(nullptr!=getDAF()) getDAF()->setDebugLvl(val); - if(val>10) genfit::MaterialEffects::getInstance()->setDebugLvl(val); m_debug = val; } +void GenfitFitter::setDebugGenfit(unsigned int val) +{ + //std::cout<<"set fitter debugGenfit "<setDebugLvl(val); +} + +void GenfitFitter::setDebugLocal(unsigned int val) +{ + //std::cout<<"set fitter debugLvlLocal "<setDebugLvlLocal(val); + //if(0==strncmp(m_fitterType.c_str(),"DAF",3)){ + //std::cout<<" GenfitFitter::setDebugLvlLocal DAF "<setDebugLvlLocal(val+1); + //} + //getDAF()->setDebugLvlLocal(); +#endif +} + void GenfitFitter::setHist(unsigned int val) {m_hist = val;} genfit::DAF* GenfitFitter::getDAF() @@ -591,21 +511,33 @@ genfit::KalmanFitterRefTrack* GenfitFitter::getKalRef() }catch(...){ if(m_debug>=3)std::cout << "dynamic_cast m_rom AbsFitter to KalmanFitterRefTrack m_ailed!" - <=2)std::cout<<"GenfitFitter::initHist "<initHist(name); + if(m_debug)std::cout<<"GenfitFitter::initHist "<initHist(name); +#endif } void GenfitFitter::writeHist() { - if(m_debug>=2)std::cout<<"GenfitFitter::writeHist "<writeHist(); + if(m_debug)std::cout<<"GenfitFitter::writeHist "<initialized()){ + genfit::GenfitHist::instance()->writeHist(); + } +#endif } +void GenfitFitter::SetRunEvent(int event){ +#ifdef GENFIT_MY_DEBUG + m_absKalman->SetRunEvent(event); +#endif +} + diff --git a/Reconstruction/RecGenfitAlg/src/GenfitFitter.h b/Reconstruction/RecGenfitAlg/src/GenfitFitter.h index 9ca466445..7b3363e6a 100644 --- a/Reconstruction/RecGenfitAlg/src/GenfitFitter.h +++ b/Reconstruction/RecGenfitAlg/src/GenfitFitter.h @@ -18,6 +18,7 @@ #include "AbsKalmanFitter.h" #include +#include "MaterialEffects.h" class GenfitField; class GenfitMaterialInterface; @@ -29,6 +30,7 @@ namespace genfit{ class AbsKalmanFitter; class KalmanFitterRefTrack; class DAF; + class MaterialEffects; } namespace dd4hep{ class OverlayedField; @@ -40,49 +42,22 @@ namespace dd4hep{ class GenfitFitter{ public: /// Type of fitters are :DAFRef,DAF,KalmanFitter,KalmanFitterRefTrack - GenfitFitter(const char* type="DAFRef",const char* name="GenfitFitter"); + GenfitFitter(const char* type="DAFRef",int debug=0,const char* name="GenfitFitter"); virtual ~GenfitFitter(); /// Magnetic field and geometry for material effect in genfit /// please SET before use !!!! void setField(const GenfitField* field); /// please SET before use !!!! - void setGeoMaterial(const dd4hep::Detector* dd4hepGeo); + void setGeoMaterial(const dd4hep::Detector* dd4hepGeo, + double extDistCut=1e-4, bool skipWireMaterial=false); /// Main fitting function - int processTrack(GenfitTrack* track, bool resort=false); + int processTrack(GenfitTrack* track, bool resort=true); /// fitting with rep int processTrackWithRep(GenfitTrack* track,int repID=0, - bool resort=false); - - /// Extrapolate the track to the CDC hit - /// Output: poca pos and dir and poca distance to the hit wire - /// Input: genfit track, pos and mom, two ends of a wire - /// pos, and mom are position & momentum at starting point - double extrapolateToHit(TVector3& poca, TVector3& pocaDir, - TVector3& pocaOnWire, double& doca, const GenfitTrack* track, - TVector3 pos, TVector3 mom, TVector3 end0, TVector3 end1, - int debug=0, int repID=0, bool stopAtBoundary=false, - bool calcJacobianNoise=false); - - /// Extrapolate the track to the cyliner at fixed raidus - /// Output: pos and mom at fixed radius - /// Input: genfitTrack, radius of cylinder at center of the origin, - /// repID, stopAtBoundary and calcAverageState - double extrapolateToCylinder(TVector3& pos, TVector3& mom, - GenfitTrack* track, double radius, const TVector3 linePoint, - const TVector3 lineDirection, int hitID =0, int repID=0, - bool stopAtBoundary=false, bool calcJacobianNoise=false); - - /// Extrapolate the track to the point - /// Output: pos and mom of POCA point to point - /// Input: genfitTrack,point,repID,stopAtBoundary and calcAverageState - /// repID same with pidType - double extrapolateToPoint(TVector3& pos, TVector3& mom, - const GenfitTrack* genfitTrack, const TVector3& point, - int repID=0, bool stopAtBoundary = false, - bool calcJacobianNoise = false) const; + bool resort=true); /// setters of fitter properties void setFitterType(const char* val); @@ -110,10 +85,16 @@ class GenfitFitter{ void setMscModelName(std::string val); void setMaterialDebugLvl(unsigned int val); void setDebug(unsigned int val); + void setDebugGenfit(unsigned int val); + void setDebugLocal(unsigned int val); void setHist(unsigned int val); /// getters of fitter properties std::string getFitterType() const {return m_fitterType;} + GenfitMaterialInterface* getGeoMaterial() const {return m_geoMaterial;} + genfit::MaterialEffects* getMaterialEffects() const { + return genfit::MaterialEffects::getInstance(); + } unsigned int getMinIterations() const { return m_minIterations; } unsigned int getMaxIterations() const { return m_maxIterations; } double getDeltaPval() const { return m_deltaPval; } @@ -137,7 +118,10 @@ class GenfitFitter{ {return m_ignoreBoundariesBetweenEqualMaterials;} std::string getMscModelName(){return m_mscModelName;} int getDebug() const {return m_debug;} + int getDebugGenfit() const {return m_debugGenfit;} + int getDebugLocal() const {return m_debugLocal;} int getHist() const {return m_hist;} + void SetRunEvent(int event); /// Printer void print(const char* name=""); @@ -194,7 +178,9 @@ class GenfitFitter{ std::string m_mscModelName; int m_debug; /// debug - int m_hist; /// hist + int m_debugGenfit; /// debug + int m_debugLocal; /// debug + int m_hist; /// hist }; #endif diff --git a/Reconstruction/RecGenfitAlg/src/GenfitHit.cpp b/Reconstruction/RecGenfitAlg/src/GenfitHit.cpp new file mode 100644 index 000000000..85af21fec --- /dev/null +++ b/Reconstruction/RecGenfitAlg/src/GenfitHit.cpp @@ -0,0 +1,94 @@ +#include "GenfitHit.h" +#include "DataHelper/TrackerHitHelper.h" +#include "DetSegmentation/GridDriftChamber.h" +#include "GenfitUnit.h" + +#include "DD4hep/Detector.h" +#include "DD4hep/DetElement.h" +#include "DD4hep/Segmentations.h" +#include "DD4hep/DD4hepUnits.h" +#include "edm4hep/TrackerHit.h" +#include "edm4hep/SimTrackerHit.h" +#include "TRandom.h" + +#include + +GenfitHit::GenfitHit(const edm4hep::TrackerHit* trackerHit, + const edm4hep::SimTrackerHit* simTrackerHit, + const dd4hep::DDSegmentation::BitFieldCoder* decoder, + const dd4hep::DDSegmentation::GridDriftChamber* gridDriftChamber, + double driftVelocity,double driftDistanceErr){ + m_trackerHit=trackerHit; + m_simTrackerHit=simTrackerHit; + m_decoder=decoder; + m_gridDriftChamber=gridDriftChamber; + m_driftVelocity=driftVelocity; + m_driftDistanceErr=driftDistanceErr*GenfitUnit::cm; + //driftVelocity um/ns + m_driftDistanceTruth=m_trackerHit->getTime()*driftVelocity*GenfitUnit::um; + m_driftDistance=m_driftDistanceTruth*GenfitUnit::cm; +// if(driftDistanceErr>0) m_driftDistance+=gRandom->Gaus(0,fabs(driftDistanceErr*GenfitUnit::cm));//FIXME +} + +unsigned long long GenfitHit::getCellID() const { + return m_trackerHit->getCellID(); +} + +int GenfitHit::getLayer() const { + return m_decoder->get(getCellID(),"layer"); +} + +int GenfitHit::getCell() const { + return m_decoder->get(getCellID(),"cellID"); +} + +int GenfitHit::getLeftRightAmbig() const { + TVector3 momTruth=getTruthMom(); + TVector3 pocaOnTrack=getTruthPos();//FIXME, not poca on track + TVector3 trackDir=momTruth.Unit(); + TVector3 wireDir=(getEnd1()-getEnd0()).Unit(); + TVector3 pocaOnWire= + m_gridDriftChamber->wirePos_vs_z(getCellID(),pocaOnTrack.Z()); + TVector3 pocaDir=(pocaOnWire-pocaOnTrack).Unit(); + //TVector3 a=pocaDir.Cross(trackDir); + int lrAmbig=(pocaDir.Cross(trackDir))*wireDir; + return fabs(lrAmbig)/lrAmbig; +} + +TVector3 GenfitHit::getEnd0() const { + TVector3 end0; + TVector3 end1; + m_gridDriftChamber->cellposition(m_trackerHit->getCellID(),end0,end1);//dd4hep unit + end0*=(1./dd4hep::cm)*GenfitUnit::cm; + return end0; +} + +TVector3 GenfitHit::getEnd1() const { + TVector3 end0; + TVector3 end1; + m_gridDriftChamber->cellposition(m_trackerHit->getCellID(),end0,end1);//dd4hep unit + end1*=(1./dd4hep::cm)*GenfitUnit::cm; + return end1; +} + +TVector3 GenfitHit::getTruthPos()const{ + edm4hep::Vector3d pos=m_simTrackerHit->getPosition(); + return TVector3(pos.x/dd4hep::cm*GenfitUnit::cm, + pos.y/dd4hep::cm*GenfitUnit::cm, + pos.z/dd4hep::cm*GenfitUnit::cm); +} + +TVector3 GenfitHit::getTruthMom()const{ + edm4hep::Vector3f mom=m_simTrackerHit->getMomentum(); + return TVector3(mom.x/dd4hep::GeV*GenfitUnit::GeV, + mom.y/dd4hep::GeV*GenfitUnit::GeV, + mom.z/dd4hep::GeV*GenfitUnit::GeV); +} + +void GenfitHit::print()const{ + //TODO + std::cout<<"driftDistanceTruth cm "< #include @@ -39,6 +39,8 @@ GenfitMaterialInterface::GenfitMaterialInterface( assert(nullptr!=dd4hepGeo); m_geoManager=&(dd4hepGeo->manager()); assert(nullptr!=m_geoManager); + m_skipWireMaterial=false; + debugLvl_=0; } GenfitMaterialInterface* GenfitMaterialInterface::getInstance( @@ -64,7 +66,6 @@ TGeoManager* GenfitMaterialInterface::getGeoManager() GenfitMaterialInterface::initTrack(double posX, double posY, double posZ, double dirX, double dirY, double dirZ) { - //debugLvl_ = 1; #ifdef DEBUG_GENFITGEO std::cout << "GenfitMaterialInterface::initTrack. \n"; std::cout << "Pos "; TVector3(posX, posY, posZ).Print(); @@ -77,11 +78,13 @@ GenfitMaterialInterface::initTrack(double posX, double posY, // Set the intended direction. setCurrentDirection(dirX, dirY, dirZ); +#ifdef DEBUG_GENFITGEO if (debugLvl_ > 0) { std::cout << " GenfitMaterialInterface::initTrack at \n"; std::cout << " position: "; TVector3(posX, posY, posZ).Print(); std::cout << " direction: "; TVector3(dirX, dirY, dirZ).Print(); } +#endif return result; } @@ -91,6 +94,9 @@ genfit::Material GenfitMaterialInterface::getMaterialParameters() { TGeoMaterial* mat = getGeoManager()->GetCurrentVolume()->GetMedium()->GetMaterial(); +#ifdef DEBUG_GENFITGEO + getGeoManager()->GetCurrentVolume()->Print(); +#endif //Scalar density; /// Density in g / cm^3 //Scalar Z; /// Atomic number //Scalar A; /// Mass number in g / mol @@ -98,6 +104,17 @@ genfit::Material GenfitMaterialInterface::getMaterialParameters() //Scalar mEE; /// Mean excitaiton energy in eV //Material from CEPCSW is NOT follow the dd4hep?FIXME + //getGeoManager()->GetCurrentVolume()->Print(); + if(m_skipWireMaterial){//FIXME + if((strncmp(getGeoManager()->GetCurrentVolume()->GetName(),"FieldWire",9) == 0) || + (strncmp(getGeoManager()->GetCurrentVolume()->GetName(),"SenseWire",9) == 0)){ + //Air: den 0.0012 radlen 30528.8 Z 7.366 A 14.7844 + std::cout<<" skipWireMateria "<GetCurrentVolume()->GetName()<Print(); return genfit::Material(mat->GetDensity(), mat->GetZ(), mat->GetA(), @@ -110,7 +127,6 @@ GenfitMaterialInterface::findNextBoundary(const genfit::RKTrackRep* rep, double sMax, // signed bool varField) { - //debugLvl_ = 1; // cm, distance limit beneath which straight-line steps are taken. const double delta(1.E-2); const double epsilon(1.E-1); // cm, allowed upper bound on arch @@ -131,6 +147,8 @@ GenfitMaterialInterface::findNextBoundary(const genfit::RKTrackRep* rep, findNextBoundary(fabs(sMax) - s); double safety = getSafeDistance(); // >= 0 double slDist = getStep(); + if (debugLvl_ > 0) + std::cout << " slDist=getStep()= " << slDist<< "; safety=getSafeDistance()=" << safety<< "\n"; // this should not happen, but it happens sometimes. // The reason are probably overlaps in the geometry. @@ -175,6 +193,11 @@ GenfitMaterialInterface::findNextBoundary(const genfit::RKTrackRep* rep, state7 = stateOrig; rep->RKPropagate(state7, nullptr, SA, stepSign*(s + step), varField); + if(debugLvl_){ + std::cout<<" RKTrackRep at state "<GetSafeDistance(); + return getGeoManager()->GetSafeDistance();//yzhang debug FIXME } double GenfitMaterialInterface::getStep() { - return getGeoManager()->GetSafeDistance(); + return getGeoManager()->GetStep();//yzhang debug FIXME } TGeoNode* GenfitMaterialInterface::findNextBoundary(double stepmax, @@ -340,8 +363,16 @@ bool GenfitMaterialInterface::isSameLocation(double posX, double posY, { //std::cout<<" MatInt "<<__LINE__<<" posXYZ*dd4hep::cm "<IsSameLocation(posX,posY, - posZ,change); +#ifdef DEBUG_GENFITGEO + std::cout<<" MatInt "<<" posXYZ "<IsSameLocation( + posX/GenfitUnit::cm*rootUnitCM, + posY/GenfitUnit::cm*rootUnitCM, + posZ/GenfitUnit::cm*rootUnitCM,change); } void GenfitMaterialInterface::setCurrentDirection(double nx, double ny, diff --git a/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.h b/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.h index 83da9fc1e..dbf602610 100644 --- a/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.h +++ b/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.h @@ -11,6 +11,7 @@ #define RECGENFITALG_GENFITMATERIALINTERFACE_H #include "AbsMaterialInterface.h" +//#include "MaterialProperties.h" class RKTrackRep; class TGeoManager; @@ -33,6 +34,9 @@ class GenfitMaterialInterface : public genfit::AbsMaterialInterface{ //Set Detector pointer of dd4hep void setDetector(dd4hep::Detector*); + //void getMaterialParameters(double& density,double& Z,double& A, + // double& radiationLength, double& mEE) {return;} + //void getMaterialParameters(genfit::MaterialProperties& parameters) {return;} /** @brief Initialize the navigator at given position and with given direction. Returns true if the volume changed. */ @@ -54,6 +58,11 @@ class GenfitMaterialInterface : public genfit::AbsMaterialInterface{ bool varField = true) override; // ClassDef(GenfitMaterialInterface, 1); + void setMinSafetyDistanceCut(double safeDistCut=1e-7) + {m_safeDistCut=safeDistCut;} + void setSkipWireMaterial(bool skipWireMateria) + {m_skipWireMaterial=skipWireMateria;} + virtual void setDebugLvl(unsigned int lvl = 1) {debugLvl_ = lvl;} private: static GenfitMaterialInterface* m_instance; @@ -66,6 +75,8 @@ class GenfitMaterialInterface : public genfit::AbsMaterialInterface{ bool isSameLocation(double posX, double posY, double posZ, bool change=false); void setCurrentDirection(double nx, double ny, double nz); + double m_safeDistCut; + bool m_skipWireMaterial; }; diff --git a/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp b/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp index 9d1da9c30..19bb23c79 100644 --- a/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp +++ b/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp @@ -1,121 +1,135 @@ #include "GenfitTrack.h" +#include "GenfitHit.h" #include "GenfitField.h" //CEPCSW #include "DataHelper/HelixClass.h" +#include "DataHelper/TrackHelper.h" +#include "DataHelper/TrackerHitHelper.h" +#include "DetInterface/IGeomSvc.h" +#include "GenfitUnit.h" #include "UTIL/ILDConf.h" +#include "WireMeasurementDC.h" //Externals #include "DD4hep/DD4hepUnits.h" -#include "edm4hep/MCParticle.h" +#include "DD4hep/Detector.h" +#include "DD4hep/DetElement.h" +#include "DD4hep/Segmentations.h" +#include "DDRec/ISurface.h" +#include "DDRec/SurfaceManager.h" +#include "DDRec/Vector3D.h" +#include "DetSegmentation/GridDriftChamber.h" +#include "edm4hep/MutableReconstructedParticle.h" +#include "edm4hep/SimTrackerHit.h" #include "edm4hep/Track.h" #include "edm4hep/MutableTrack.h" #include "edm4hep/TrackerHit.h" -#include "edm4hep/SimTrackerHit.h" -#include "edm4hep/ReconstructedParticle.h" -#include "edm4hep/MutableReconstructedParticle.h" #include "edm4hep/TrackerHitCollection.h" +#include "edm4hep/MCParticle.h" #include "edm4hep/MCRecoTrackerAssociationCollection.h" #include "edm4hep/Vector3d.h" -#include "DetSegmentation/GridDriftChamber.h" +#include "GaudiKernel/SmartIF.h" //genfit -#include "Track.h" -#include "MeasuredStateOnPlane.h" -#include "RKTrackRep.h" -#include "TrackPoint.h" -#include "StateOnPlane.h" -#include "KalmanFitterInfo.h" -#include "KalmanFittedStateOnPlane.h" +#include "AbsMeasurement.h" #include "AbsTrackRep.h" #include "FitStatus.h" +#include "KalmanFitterInfo.h" +#include "KalmanFittedStateOnPlane.h" +#include "RKTrackRep.h" +//#include "PlanarMeasurement.h" +#include "PlanarMeasurementSDT.h" #include "SpacepointMeasurement.h" +#include "StateOnPlane.h" +#include "RectangularFinitePlane.h" +#include "Track.h" +#include "TrackPoint.h" +#include "MeasuredStateOnPlane.h" #include "WireMeasurementNew.h" +#include "AbsMeasurement.h" +#include "TrackPoint.h" //ROOT -#include "TRandom.h" -#include "TVector3.h" #include "TLorentzVector.h" #include "TMatrixDSym.h" +#include "TRandom.h" +#include "TVector3.h" //cpp #include #include +#undef GENFIT_MY_DEBUG +//#define GENFIT_MY_DEBUG 1 + const int GenfitTrack::s_PDG[2][5] ={{-11,-13,211,321,2212},{11,13,-211,-321,-2212}}; bool -sortDCHit(edm4hep::SimTrackerHit hit1,edm4hep::SimTrackerHit hit2) +sortDCSimHit(edm4hep::SimTrackerHit hit1,edm4hep::SimTrackerHit hit2) { - //std::cout<<"hit1"< hitPair1,std::pair hitPair2) { + bool isEarly=hitPair1.first geom, const char* name) +:m_name(name),m_track(nullptr),m_debug(0),m_debugLocal(0),m_geomSvc(geom), + m_genfitField(genfitField),m_gridDriftChamber(seg), + m_decoderDC(geom->getDecoder("DriftChamberHitsCollection")) +{ } GenfitTrack::~GenfitTrack() { + clearGenfitHitVec(); ///Note: track reps and points will be deleted in the destructor of track ///implemented in genfit::Track::Clear() delete m_track; } -/// create a Genfit track from track state, without trackRep -/// Initialize track with seed states -/// NO unit conversion here +void GenfitTrack::clearGenfitHitVec(){ + for(auto h:m_genfitHitVec){delete h;} + m_genfitHitVec.clear(); + std::vector().swap(m_genfitHitVec); +} + +/// create a Genfit track from track state +/// Initialize track with seed state and cov +/// unit conversion here!!! bool GenfitTrack::createGenfitTrack(int pdgType,int charge, - TLorentzVector posInit, TVector3 momInit, TMatrixDSym covMInit) + TVectorD trackParam, TMatrixDSym covMInit_6) { - TVectorD seedState(6); - TMatrixDSym seedCov(6); - - //for(int i = 0; i < 6; ++i) { - // for(int j = 0; j < 6; ++j) { - // seedCov(i,j)=covMInit(i,j); - // } - //} - //yzhang FIXME - //seed position - for(int i = 0; i < 3; ++i) { - seedState(i)=posInit[i]; - //yzhang TODO from covMInit to seedCov - double resolution = 0.1;//*dd4hep::mm/dd4hep::cm; - seedCov(i,i)=resolution*resolution; - if(i==2) seedCov(i,i)=0.5*0.5; - } - //seed momentum - for(int i = 3; i < 6; ++i){ - //seedState(i)=momInit[i-3]*(dd4hep::GeV); - seedState(i)=momInit[i-3]; - //yzhang TODO from covMInit to seedCov - seedCov(i,i)=0.01;//pow(resolution / sqrt(3),2); - } - - if(nullptr==m_track) m_track=new genfit::Track(); - m_track->setStateSeed(seedState); - m_track->setCovSeed(seedCov); - - /// new a track representation and add to the track - int chargeId=0; - charge>0 ? chargeId=0 : chargeId=1; + if(m_debug){ + std::cout<<"createGenfitTrack pdgType "<setStateSeed(trackParam); + m_track->setCovSeed(covMInit_6); - if(m_debug>=2)std::cout<=2)std::cout<<"seedCov "<setDebugLvlLocal(m_debugLocal); +#endif - addTrackRep(s_PDG[chargeId][pdgType]); + ///new a track representation and add to the track + addTrackRep(pdgType,charge); - if(m_debug>0) seedCov.Print(); + if(m_debug>0) printSeed(); return true; } @@ -123,228 +137,321 @@ bool GenfitTrack::createGenfitTrack(int pdgType,int charge, bool GenfitTrack::createGenfitTrackFromMCParticle(int pidType, const edm4hep::MCParticle& mcParticle, double eventStartTime) { - ///get track parameters from McParticle - const auto& mcPocaPos = mcParticle.getVertex();//mm - const auto& mcPocaMom = mcParticle.getMomentum();//GeV - if(m_debug>=2)std::cout<<"seedPos poca "<< mcPocaPos.x - <<" "<=2)std::cout<<"seedMom poca "<< mcPocaMom.x - <<" "<>; - MomentumT firstLayerMom{1e9,1e9,1e9}; - pivotToFirstLayer(mcPocaPos,mcPocaMom,firstLayerPos,firstLayerMom); - - //TODO convert unit - ///Get seed position and momentum - TLorentzVector seedPos(firstLayerPos.x,firstLayerPos.y,firstLayerPos.z, - eventStartTime); - TVector3 seedMom(firstLayerMom.x,firstLayerMom.y,firstLayerMom.z); - if(m_debug>=2)std::cout<<"seedPos "<< firstLayerPos.x - <<" "<=2)std::cout<<"seedMom "<< firstLayerMom.x - <<" "<=2)std::cout<<"createGenfitTrack " ; - if(!GenfitTrack::createGenfitTrack(pidType,mcParticle.getCharge(), - seedPos,seedMom,covM)){ - if(m_debug>=2)std::cout<<"GenfitTrack" - <<" Error in createGenfitTrack" <=2)std::cout<<"GenfitTrack " - <<"createGenfitTrackFromMCParticle track created" <getBz({0.,0.,0.})/dd4hep::kilogauss<getBz({0.,0.,0.})/dd4hep::tesla<getBz({0.,0.,0.})<< dd4hep::kilogauss <<" "<getBz({0.,0.,0.})*dd4hep::kilogauss/dd4hep::tesla); - TLorentzVector posInit(helixClass.getReferencePoint()[0], - helixClass.getReferencePoint()[1], - helixClass.getReferencePoint()[2],eventStartTime); - posInit.SetX(posInit.X()*dd4hep::mm); - posInit.SetY(posInit.Y()*dd4hep::mm); - posInit.SetZ(posInit.Z()*dd4hep::mm); - TVector3 momInit(helixClass.getMomentum()[0], - helixClass.getMomentum()[1],helixClass.getMomentum()[2]); - momInit.SetX(momInit.x()*dd4hep::GeV); - momInit.SetY(momInit.y()*dd4hep::GeV); - momInit.SetZ(momInit.z()*dd4hep::GeV); - TMatrixDSym covMInit; - if(!createGenfitTrack(pidType, - int(trackStat.omega/fabs(trackStat.omega)),//charge,FIXME? - posInit,momInit,covMInit)){ + + ///Skip track w.o. hit + if(track.trackerHits_size()<=0) { + if(m_debug){ + std::cout<<"createGenfitTrackFromEDM4HepTrack skip track w/o hit"< sigmaUVec,std::vector sigmaVVec,int cellID,int hitID) { - edm4hep::Vector3d pos=hit.getPosition(); - TVectorD p(3); - p[0]=pos.x*dd4hep::mm; - p[1]=pos.y*dd4hep::mm; - p[2]=pos.z*dd4hep::mm; + TVectorD pos_smeared(3); + for(int i=0;i<3;i++) pos_smeared[i]=pos(i)*GenfitUnit::mm; - if(m_debug>=2)std::cout<=2)std::cout<<"cov "<insertPoint(trackPoint); + /// space point error matrix and smear, unit cm + TMatrixDSym hitCov(3); + hitCov.Zero(); + + //get sigmas + float sigmaU,sigmaV; + getSigmas(cellID,sigmaUVec,sigmaVVec,sigmaU,sigmaV); + float sigmaX=sigmaU;//sigmaU*cos(atan2(pos_smeared[1],pos_smeared[0])); + float sigmaY=sigmaU;//*sin(atan2(pos_smeared[1],pos_smeared[0])); + float sigmaZ=sigmaV; + + //smear 3d track position + bool smear=(sigmaUVec[0]>0); + if(smear){ + pos_smeared[0]+=gRandom->Gaus(0,sigmaX); + pos_smeared[1]+=gRandom->Gaus(0,sigmaX); + pos_smeared[2]+=gRandom->Gaus(0,sigmaX); + } + hitCov(0,0)=sigmaX*sigmaX;//FIXME + hitCov(1,1)=sigmaX*sigmaX; + hitCov(2,2)=sigmaX*sigmaX; + + if(m_debug>=2){ + std::cout<<" addSpacePointMeasurement pos "<=2)std::cout<<"end of addSpacePointTrakerHit"<=2){ + std::cout<Print(); + } + m_track->insertPoint(trackPoint); return true; +}//end of addSpacePointMeasurement + +/// Return isurface of a silicon hit +const dd4hep::rec::ISurface* +GenfitTrack::getISurface(edm4hep::TrackerHit* hit){ +//GenfitTrack::getISurface(edm4hep::ConstTrackerHit* hit){ + dd4hep::rec::SurfaceManager surfaceManager(*m_geomSvc->lcdd()); + + std::string detectorName; + unsigned long long cellID=hit->getCellID(); + int detTypeID=getDetTypeID(hit->getCellID()); + if(detTypeID==lcio::ILDDetID::VXD){ + detectorName="VXD"; + if(hit->getPosition().z>0){ + cellID+=0x100000000; + }else{ + cellID+=0x300000000; + } + }else if(detTypeID==lcio::ILDDetID::SIT){ + detectorName="SIT"; + }else if(detTypeID==lcio::ILDDetID::SET){ + detectorName="SET"; + }else if(detTypeID==lcio::ILDDetID::FTD){ + detectorName="FTD"; + }else{ + std::cout << "ERROR:getISurface iSurface = NULL!" << std::endl; + return nullptr; + } + if(m_debug>2) std::cout<getPosition()<<" hit.cellID "<getCellID() + <<" cellId "<find(cellID); + dd4hep::rec::ISurface* iSurface=nullptr; + if(iter!=surfaceMap->end()){iSurface=(*iter).second;} + + return iSurface; } -/// Add a 3d SpacepointMeasurement with MC truth position smeared by sigma -bool GenfitTrack::addSpacePointMeasurement(const TVectorD& pos, - double sigma, int detID, int hitID, bool smear) +/// Add a 1d strip or 2d pixel smeared by sigma + bool +GenfitTrack::addSiliconMeasurement(edm4hep::TrackerHit* hit, + float sigmaU,float sigmaV,int cellID,int hitID) { - double sigma_t=sigma*dd4hep::mm; - /// Convert from CEPCSW unit to genfit unit, cm - TVectorD pos_t(3); - pos_t(0)=pos(0)*dd4hep::mm; - pos_t(1)=pos(1)*dd4hep::mm; - pos_t(2)=pos(2)*dd4hep::mm; + if(m_debug>0) std::cout<<"addSiliconMeasurement "<<*hit<Gaus(0,sigma_t/TMath::Sqrt(3.)); + ///get surface by cellID + const dd4hep::rec::ISurface* iSurface = getISurface(hit); + if(nullptr==iSurface){ + std::cout<=2)std::cout<insertPoint(trackPoint); + ///Get detector plane parameter + //VXD + //SIT + //SET + + ///Get measurement and cov + TVectorD hitCoords(2); + TVector3 p; + TMatrixDSym hitCov(2); + getMeasurementAndCov(hit,p,hitCov); + hitCoords(0)=(p-o).Dot(u); + hitCoords(1)=(p-o).Dot(v); + if(m_debug) std::cout<<"yzhang debug hitCoords cm "<setFinitePlane(pixelOrStripPlane); + genfit::PlanarMeasurementSDT* planarMeasurement=new genfit::PlanarMeasurementSDT( + hitCoords,hitCov,cellID,hitID,nullptr); + planarMeasurement->setTrackerHit(hit); + planarMeasurement->setPlane(plane); + m_track->insertPoint(new genfit::TrackPoint(planarMeasurement,m_track)); + + if(m_debug>2){ + std::cout<<"hitID "< sigmaUVec,std::vector sigmaVVec) +{ + ///Get TrackerHit on Track + int nHitAdd=0; + for(unsigned int iHit=0;iHitgetCellID(); + float sigmaU,sigmaV; + int sigmaUID=getSigmas(cellID,sigmaUVec,sigmaVVec,sigmaU,sigmaV); + if(0==sigmaUID) continue;//skip DC track + addSiliconMeasurement(trackerHit,sigmaU,sigmaV,cellID,nHitAdd++); + GenfitHit* genfitHit= + new GenfitHit(trackerHit,nullptr,nullptr, + nullptr,0.,0.); + if(nullptr==genfitHit) continue; + m_genfitHitVec.push_back(genfitHit); + } + return nHitAdd; +} -/// Add a WireMeasurement, no Unit conversion here -void GenfitTrack::addWireMeasurement(double driftDistance, - double driftDistanceError, const TVector3& endPoint1, - const TVector3& endPoint2, int lrAmbig, int detID, int hitID) +//Add wire measurement on wire, unit conversion here +int GenfitTrack::addWireMeasurementsFromListTrF(const edm4hep::TrackerHitCollection* trkHits, + float sigma,int sortMethod) { - try{ - /// New a WireMeasurement - genfit::WireMeasurementNew* wireMeas = new genfit::WireMeasurementNew( - driftDistance, driftDistanceError, endPoint1, endPoint2, detID, - hitID, nullptr); - wireMeas->setMaxDistance(0.5);//0.5 cm FIXME - wireMeas->setLeftRightResolution(lrAmbig); - - if(m_debug>=2)std::cout< hits; + for(auto trkHit:*trkHits){ + hits.push_back(&trkHit); + } + + std::vector sortedTrackerHits; + getSortedTrackerHitsTrF(hits,sortedTrackerHits,sortMethod); + + if(m_debug>0){ + std::cout<<"n sortedTrackerHits "<getCellID(),end0,end1); ///New a TrackPoint,create connection between meas. and trackPoint - genfit::TrackPoint* trackPoint=new genfit::TrackPoint(wireMeas,m_track); - wireMeas->setTrackPoint(trackPoint); + WireMeasurementDC* wireMeasurementDC= + new WireMeasurementDC(1e-3*driftVelocity*trackerHit->getTime(),sigma, + end0,end1,getDetTypeID(trackerHit->getCellID()),nHitAdd,nullptr); - m_track->insertPoint(trackPoint); + int layer = m_decoderDC->get(trackerHit->getCellID(),"layer"); + int cellID = m_decoderDC->get(trackerHit->getCellID(),"cellID"); + const edm4hep::TrackerHit trackerHit_; + wireMeasurementDC->setTrackerHit(trackerHit_,layer,cellID); - }catch(genfit::Exception& e){ - if(m_debug>=2)std::cout<insertPoint(new genfit::TrackPoint(wireMeasurementDC,m_track)); + if(m_debug>=2){ std::cout<Print(); } + nHitAdd++; + }//end of loop over sotred hits + return nHitAdd; +}//end of addWireMeasurementsFromList + + +//Add wire measurement on wire, unit conversion here +//int GenfitTrack::addWireMeasurementsFromList(podio::RelationRange hits,float sigma, +int GenfitTrack::addWireMeasurementsFromList(std::vector& hits,float sigma, + const edm4hep::MCRecoTrackerAssociationCollection* assoHits, + int sortMethod, bool truthAmbig,float skipCorner,float skipNear) +{ + if(m_debug>0){ std::cout<<"addWireMeasurementsFromList"< sortedTrackerHits; + getSortedTrackerHits(hits,assoHits,sortedTrackerHits,sortMethod); + + if(m_debug>0){ + std::cout<<"n sortedTrackerHits "<insertPoint(new genfit::TrackPoint(wireMeasurementDC,m_track)); + if(m_debug>=2){ std::cout<Print(); } + }//end of loop over sotred hits + return nHitAdd; +}//end of addWireMeasurementsFromList + //Add wire measurement on wire, unit conversion here -bool GenfitTrack::addWireMeasurementOnTrack(edm4hep::Track track,double sigma) +int GenfitTrack::addWireMeasurementsOnTrack(edm4hep::Track& track,float sigma, + const edm4hep::MCRecoTrackerAssociationCollection* assoHits, + int sortMethod, bool truthAmbig,float skipCorner,float skipNear) { - for(unsigned int iHit=0;iHit0){ std::cout<<"addWireMeasurementsOnTrack"<cellposition(hit.getCellID(),endPointStart, - endPointEnd); - int lrAmbig=0; - if(m_debug>=2)std::cout<=2)std::cout<0){ + std::cout<<"n trackerHits size"<0) std::cout<<"No hit on track"< hits_t=track.getTrackerHits(); + std::vector hits; + //hits.reserve(1000); + for(auto &h:hits_t){ + edm4hep::TrackerHit* hit = const_cast(&h); + hits.push_back(hit); + } + nHitAdd=addWireMeasurementsFromList(hits,sigma,assoHits,sortMethod,truthAmbig,skipCorner,skipNear); + + return nHitAdd; +}//end of addWireMeasurements /// Get MOP -bool GenfitTrack::getMOP(int hitID, - genfit::MeasuredStateOnPlane& mop, genfit::AbsTrackRep* trackRep) const +bool GenfitTrack::getMOP(int hitID,genfit::MeasuredStateOnPlane& mop, + genfit::AbsTrackRep* trackRep) const { + if(nullptr == trackRep) trackRep = getRep(); try{ mop = m_track->getFittedState(hitID,trackRep); @@ -356,30 +463,22 @@ bool GenfitTrack::getMOP(int hitID, } /// New and add a track representation to track -genfit::RKTrackRep* GenfitTrack::addTrackRep(int pdg) +genfit::RKTrackRep* GenfitTrack::addTrackRep(int pdgType,int charge) { - /// create a new track representation - genfit::RKTrackRep* rep = new genfit::RKTrackRep(pdg); - m_reps.push_back(rep); + int chargeId=0; + charge>0 ? chargeId=0 : chargeId=1;//s_PDG[0]: positive particle + + genfit::RKTrackRep* rep=new genfit::RKTrackRep(s_PDG[chargeId][pdgType]); m_track->addTrackRep(rep); - //try{ - // genfit::MeasuredStateOnPlane stateInit(rep); - // rep->setPosMomCov(stateInit, pos, mom, covM); - //}catch(genfit::Exception e){ - // if(m_debug>=2)std::cout<getStateSeed(); - TVector3 p(seedStat[0],seedStat[1],seedStat[2]); + TVectorD seedState(6); + seedState = m_track->getStateSeed(); + TVector3 p(seedState[0],seedState[1],seedState[2]); p = p*dd4hep::cm; TLorentzVector pos(p[0],p[1],p[2],9999);//FIXME return pos; @@ -388,18 +487,18 @@ const TLorentzVector GenfitTrack::getSeedStatePos()const /// Get the momentum from genfit::Track::getStateSeed const TVector3 GenfitTrack::getSeedStateMom() const { - TVectorD seedStat(6); seedStat = m_track->getStateSeed(); - TVector3 mom(seedStat[3],seedStat[4],seedStat[5]); + TVectorD seedState(6); seedState = m_track->getStateSeed(); + TVector3 mom(seedState[3],seedState[4],seedState[5]); return mom*dd4hep::GeV; } /// Get the seed states of momentum and position void GenfitTrack::getSeedStateMom(TLorentzVector& pos, TVector3& mom) const { - TVectorD seedStat(6); seedStat = m_track->getStateSeed(); - mom = TVector3(seedStat[3],seedStat[4],seedStat[5])*dd4hep::GeV; - seedStat = m_track->getStateSeed(); - TVector3 p = TVector3(seedStat[0],seedStat[1],seedStat[2])*dd4hep::cm; + TVectorD seedState(6); seedState = m_track->getStateSeed(); + mom = TVector3(seedState[3],seedState[4],seedState[5])*dd4hep::GeV; + seedState = m_track->getStateSeed(); + TVector3 p = TVector3(seedState[0],seedState[1],seedState[2])*dd4hep::cm; pos.SetXYZT(p[0],p[1],p[2],9999);//FIXME time } @@ -408,6 +507,23 @@ unsigned int GenfitTrack::getNumPoints() const return m_track->getNumPoints(); } +GenfitHit* GenfitTrack::GetHit(long unsigned int i) const { + if(i>=m_genfitHitVec.size())return nullptr; + return m_genfitHitVec[i]; +} + +unsigned int GenfitTrack::getNumPointsDet(int detTypeID) const +{ + unsigned int nHit=0; + const std::vector tps=m_track->getPoints(); + for(auto tp:tps){ + const genfit::AbsMeasurement* m=nullptr; + if(tp->hasRawMeasurements()) m=tp->getRawMeasurement(); + if(nullptr!=m && detTypeID==getDetTypeID(m->getDetId())) nHit++; + } + return nHit; +} + /// Test the fit result FIXME bool GenfitTrack::fitSuccess(int repID) const { @@ -419,8 +535,8 @@ bool GenfitTrack::fitSuccess(int repID) const ||fitStatus->isFitConvergedFully()) { if(m_debug>=2)std::cout<isFitted()<<" , isFitConverged " - <isFitConverged()<<", isFitConvergedFully " - <isFitConvergedFully()<isFitConverged()<<", isFitConvergedFully " + <isFitConvergedFully()<setDebugLvl(debug); } + if(m_track){ + for(unsigned int i=0;igetNumReps();i++){ + m_track->getTrackRep(i)->setDebugLvl(debug); + } + } +} +void GenfitTrack::setDebugLocal(int debug){ +#ifdef GENFIT_MY_DEBUG + if(m_track){ + for(unsigned int i=0;igetNumReps();i++){ + } + } + m_debugLocal=debug; +#endif + if(m_debug)std::cout<<"GenfitTrack:setDebugLvlLocal "<getCovSeed(); + std::cout << " covSeed = " << std::endl; + covSeed.Print(); } void GenfitTrack::printFitted(int repID) const @@ -481,16 +616,15 @@ void GenfitTrack::print( TLorentzVector pos, TVector3 mom, if(m_debug>=2)std::cout<1){ - for(unsigned int i=0;iPrint(); } + for(unsigned int i=0;igetNumReps();i++){ + m_track->getTrackRep(i)->Print(); + } } - //for(unsigned int i=0; igetNumPoints(); i++){ - // m_track->getPoint(i)->print(); - //} } /// Get position, momentum, cov on plane of hitID-th hit bool GenfitTrack::getPosMomCovMOP(int hitID, TLorentzVector& pos, - TVector3& mom, TMatrixDSym& cov, int repID) const + TVector3& mom, TMatrixDSym& cov,int repID) const { TVector3 p; genfit::MeasuredStateOnPlane mop; @@ -507,7 +641,9 @@ bool GenfitTrack::getPosMomCovMOP(int hitID, TLorentzVector& pos, int GenfitTrack::getNumPointsWithFittedInfo(int repID) const { int nHitWithFittedInfo = 0; + //check number of hit with fitter info int nHit = m_track->getNumPointsWithMeasurement(); + //check number of hit with fitter info for(int i=0; igetPointWithFitterInfo(i,getRep(repID))){ nHitWithFittedInfo++; @@ -517,7 +653,7 @@ int GenfitTrack::getNumPointsWithFittedInfo(int repID) const } int GenfitTrack::getFittedState(TLorentzVector& pos, TVector3& mom, - TMatrixDSym& cov, int repID, bool biased) const + TMatrixDSym& cov, int trackPointId, int repID, bool biased) const { //check number of hit with fitter info if(getNumPointsWithFittedInfo(repID)<=2) return 1; @@ -529,14 +665,13 @@ int GenfitTrack::getFittedState(TLorentzVector& pos, TVector3& mom, //get first or last measured state on plane genfit::MeasuredStateOnPlane mop; try{ - mop = m_track->getFittedState(biased); + mop = m_track->getFittedState(trackPointId,rep,biased); }catch(genfit::Exception& e){ - std::cout<<" getNumPointsWithFittedInfo " + std::cout<<" getNumPointsWithFittedInfo=" <=2)std::cout<getPDG(); + return m_track->getTrackRep(id)->getPDG(); } int GenfitTrack::getPDGCharge(int id) const { - return m_reps[id]->getPDGCharge(); + return m_track->getTrackRep(id)->getPDGCharge(); } const genfit::FitStatus* @@ -574,258 +709,916 @@ GenfitTrack::getFitStatus(int repID) const } /// Extrapolate track to the cloest point of approach(POCA) to the wire of hit, -/// return StateOnPlane of this POCA -/// inputs -/// pos,mom ... position & momentum at starting point, unit= [mm]&[GeV/c] -/// (converted to [cm]&[GeV/c] inside this function) -/// hit ... destination -/// outputs poca [mm] (converted from [cm] in this function) ,global -double GenfitTrack::extrapolateToHit( TVector3& poca, TVector3& pocaDir, - TVector3& pocaOnWire, double& doca, TVector3 pos, TVector3 mom, - TVector3 end0,//one end of the hit wire - TVector3 end1,//the orhter end of the hit wire - int debug, - int repID, - bool stopAtBoundary, - bool calcJacobianNoise)const +/// return doca, poca on wire and track +/// outputs poca [cm] +/// unit conversion here +double GenfitTrack::extrapolateToHit(TVector3& poca, TVector3& pocaDir, + TVector3& pocaOnWire, double& doca,edm4hep::MCParticle mcParticle, + int cellID, int repID, bool stopAtBoundary, bool calcJacobianNoise)const { - - //genfit::MeasuredStateOnPlane state = getMOP(iPoint); // better? - //genfit::MeasuredStateOnPlane state = getMOP(0); // better? - //To do the extrapolation without IHitSelection,above 2 lines are not used. - pos = pos*dd4hep::cm;//FIXME - mom = mom*dd4hep::GeV; - - //std::cout<<__LINE__<<" extrapolate to Hit pos and mom"<getPDG()); - rep->setDebugLvl(debug); - genfit::MeasuredStateOnPlane state(rep); + ///Get MCParticle and convert unit + TVector3 pos; + TVector3 mom; + getPosMomFromMCPartile(mcParticle,pos,mom); + + TVector3 end0(0,0,0); + TVector3 end1(0,0,0); + getEndPointsOfWire(cellID,end0,end1); + + int chargeId; + mcParticle.getCharge() >0 ? chargeId=0 : chargeId=1;//s_PDG[0]: positive particle + genfit::RKTrackRep* rep = new genfit::RKTrackRep(s_PDG[chargeId][repID]); + genfit::StateOnPlane state(rep); rep->setPosMom(state, pos, mom); - - //genfit::MeasuredStateOnPlane state(m_track->getRep(repID)); - //m_track->getRep(repID)->setPosMom(state, pos, mom); - - //m_track->PrintSeed(); - double extrapoLen(0); - //std::cout<<" wire1 "< (m_track->getRep(repID)); - //extrapoLen = rkRep->extrapolateToLine(state, end0, end1, poca, - // pocaDir, pocaOnWire, stopAtBoundary, calcJacobianNoise); - extrapoLen = rep->extrapolateToLine(state, end0, end1, poca, - pocaDir, pocaOnWire, stopAtBoundary, calcJacobianNoise); + extrapoLen = rep->extrapolateToLine(state,end0,end0-end1, + poca,pocaDir,pocaOnWire,stopAtBoundary,calcJacobianNoise); }catch(genfit::Exception& e) { if(m_debug>=3)std::cout<< - "Exception in GenfigFitter::ExtrapolateToHit"<=2)std::cout<< " poca "<=2)std::cout<< " pocaOnWire "<=2){ + std::cout<< " poca "< sigmaU,std::vector sigmaV) +{ + ///Get TrackerHit on Track + int nHitAdd=0; + for(unsigned int iHit=0;iHitlcdd()->constant("DetID_DC")==detTypeID) continue; + edm4hep::TrackerHit hit=track.getTrackerHits(iHit); + edm4hep::Vector3d pos=hit.getPosition(); + + TVector3 p(pos.x,pos.y,pos.z); + + p.Print(); + unsigned long long cellID = hit.getCellID(); + if(addSpacePointMeasurement(p,sigmaU,sigmaV,cellID,nHitAdd)){ + if(m_debug>=2)std::cout<<"add silicon space point"<=2)std::cout<<"addSpacePointMeasurement" + <=2) std::cout< sortedDCTrackHitCol; for(unsigned int iHit=0;iHit=2)std::cout<=2)std::cout<<"add slicon space point"<=2)std::cout<<"silicon addSpacePointTrakerHit" - <=2)std::cout<<"add DC space point"<=2)std::cout<<"addSpacePointMeasurement" - // <size();iSimHit++){ - if(assoHits->at(iSimHit).getRec()==hit && - assoHits->at(iSimHit).getSim().getTime()at(iSimHit).getSim(); - minTime=assoHits->at(iSimHit).getSim().getTime(); - } + int detTypeID=getDetTypeID(track.getTrackerHits(iHit).getCellID()); + if(m_geomSvc->lcdd()->constant("DetID_DC")!=detTypeID) continue; + + ///Select the SimTrakerHit with least time + float minTime=FLT_MAX; + edm4hep::SimTrackerHit minTimeSimHit; + for(int iSimHit=0;iSimHit<(int) assoHits->size();iSimHit++){ + if(assoHits->at(iSimHit).getRec()==hit && + assoHits->at(iSimHit).getSim().getTime()at(iSimHit).getSim(); + minTime=assoHits->at(iSimHit).getSim().getTime(); } - //std::cout<<"minTimeSimHit "<=2)std::cout<<" Skip add this hit!"<=2)std::cout<<" addSimTrakerHits trackerHits_size=" - <=2)std::cout<<"add DC space point"<=2)std::cout<<"addSpacePointMeasurement" - <=2) std::cout<getTrack()->getFittedState(0,rep)))); + + // get track point with fitter info + genfit::TrackPoint* tp = getTrack()->getPointWithFitterInfo(0,rep); + if(nullptr == tp) { + if(m_debug>=2)std::cout<< + "In extrapolateToPoint TrackPoint is null"<( + tp->getFitterInfo(rep))->getBackwardUpdate(); + genfit::StateOnPlane orignalState(*state); + if(m_debug>3){ + tp->Print(); + std::cout<<" original state before extrapolate "<Print(); + } + + if(nullptr == state) { + if(m_debug>=2)std::cout<< + "In extrapolateToPoint KalmanFittedStateOnPlane is null"<extrapolateToPoint(*state, + point*(1/dd4hep::cm),stopAtBoundary, calcJacobianNoise); + + rep->getPosMomCov(*state,pos,mom,cov);//FIXME exception exist + + pos = pos*dd4hep::cm;//FIXME unit + mom = mom*dd4hep::GeV;//FIXME unit + if(m_debug>3){ + std::cout<<" original state before extrapolate "<Print(); + } + + } catch(genfit::Exception& e){ + if(m_debug>=3)std::cout + <<"Exception in GenfitTrack::extrapolateToPoint" + << e.what()<isFitted()) return trackLength; + try{ + // get track rep + genfit::AbsTrackRep* rep = getRep(repID); + if(nullptr == rep) { + if(m_debug>=2)std::cout<<"In extrapolateToCylinder rep is null" + <getPointWithFitterInfo(hitID,rep); + if(nullptr == tp) { + if(m_debug>=2)std::cout<< + "In extrapolateToCylinder TrackPoint is null"<=2)std::cout<<"GenfitTrack nHitAdded "<( + tp->getFitterInfo(rep))->getBackwardUpdate(); + + if(nullptr == state){ + if(m_debug>=2)std::cout<<"In extrapolateToCylinder " + <<"no KalmanFittedStateOnPlane in backwardUpdate"<setPosMom(*state, pos*(1/dd4hep::cm), mom*(1/dd4hep::GeV));//FIXME unit + + /// extrapolate + trackLength = rep->extrapolateToCylinder(*state, + radius/dd4hep::cm, linePoint*(1/dd4hep::cm), lineDirection, + stopAtBoundary, calcJacobianNoise); + // get pos&mom at extrapolated point on the cylinder + rep->getPosMom(*state,pos,mom);//FIXME exception exist + pos = pos*dd4hep::cm; + mom = mom*dd4hep::GeV; + } catch(genfit::Exception& e){ + if(m_debug>=3)std::cout + <<"Exception in GenfitTrack::extrapolateToCylinder " + << e.what()<& smearDistance, + std::vector& truthDistance,double driftVelocity) +{ + + int DCHit = 0; + int SDTHit = 0; + unsigned int nPoints = m_track->getNumPoints(); + for(unsigned int i = 0; igetPoint(i); + genfit::AbsMeasurement* absMea = point->getRawMeasurement(); + genfit::PlanarMeasurementSDT* sdtMea = + dynamic_cast(absMea); + if(sdtMea){ + const edm4hep::TrackerHit* TrackerHit_ = sdtMea->getTrackerHit(); + SDTHit++; + }else{ + WireMeasurementDC* dcMea = + dynamic_cast(absMea); + if(dcMea){ + const edm4hep::TrackerHit* TrackerHit_ = dcMea->getTrackerHit(); + smearDistance.push_back(1e-3*driftVelocity*TrackerHit_->getTime()); + DCHit++; + for(auto dcDigi: *dCDigiCol){ + if(dcDigi.getCellID() == TrackerHit_->getCellID()) + { + truthDistance.push_back(1e-3*driftVelocity*(dcDigi.getTime())); + } + } + } + } + } + + DCHit = nFittedDC; + SDTHit = nFittedSDT; + ngenfitHit = nFittedSDT+nFittedDC; + + return true; +} + +/// Get doca on plane of hitID-th hit +bool GenfitTrack::GetDocaRes(int hitID, double& DriftDis,double& fittedDoca, + double& res, int repID, bool biased) const +{ + genfit::TrackPoint* trackPoint = m_track->getPointWithFitterInfo(hitID, getRep(repID)); + if(!trackPoint)return false; + genfit::AbsMeasurement* thismea = trackPoint->getRawMeasurement(0); + const TVectorD measereddoca = thismea->getRawHitCoords(); + DriftDis = abs(measereddoca[0]); + genfit::MeasuredStateOnPlane mop; + if(!getMOP(hitID,mop,getRep(repID))) return false; + TVectorD state = mop.getState(); + if(state.GetNrows() != 5)return false; + fittedDoca = abs(state[3]); + res = 10*(fittedDoca - DriftDis); + return true; } bool GenfitTrack::storeTrack(edm4hep::MutableReconstructedParticle& recParticle, - int pidType, int ndfCut, double chi2Cut) + edm4hep::MutableTrack& track,TVector3& pocaToOrigin_pos, + TVector3& pocaToOrigin_mom,TMatrixDSym& pocaToOrigin_cov, + // edm4hep::TrackState& trackState, + int pidType, int ndfCut, double chi2Cut, + int& nFittedDC, int& nFittedSDT, int& ngenfitHit, + std::vector& trackL, std::vector& hitMom, + std::vector& truthMomEdep, + const edm4hep::MCRecoTrackerAssociationCollection* assoHits, + std::vector& driftDis, + std::vector& FittedDoca, + std::vector& truthDoca, + std::vector& Res, + std::vector& truthRes) { + int id = 0; + int fitid = 0; + int DCHit = 0; + int dcFit = 0; + int sdtFit = 0; + int SDTHit = 0; + + // getNumRawMeasurements + unsigned int nPoints = m_track->getNumPoints(); + unsigned int nPointsWithMea = m_track->getNumPointsWithMeasurement(); + + std::vector hitMomMag; + + while(1){ + genfit::TrackPoint* point = m_track->getPointWithFitterInfo(id); + if(!point)break; + id++; + genfit::AbsMeasurement* absMea = point->getRawMeasurement(); + + // getPoint FitterInfo + genfit::AbsFitterInfo* FitterInfo = point->getFitterInfo(); + genfit::KalmanFitterInfo *KalmanfitterInfo = + dynamic_cast(FitterInfo); + + unsigned int nNumMea = KalmanfitterInfo->getNumMeasurements(); + + bool flag = false; + for(int i=0;igetMeasurementOnPlane(i); + + double weight = MeaOnPlane->getWeight(); + if(weight>0.8) flag = true; + } + if(flag) fitid++; + + genfit::MeasurementOnPlane * Mea = + KalmanfitterInfo->getMeasurementOnPlane(); + + genfit::PlanarMeasurementSDT* sdtMea = + dynamic_cast(absMea); + if(sdtMea){ + SDTHit++; + if(flag) sdtFit++; + const edm4hep::TrackerHit* TrackerHit_ = sdtMea->getTrackerHit(); + track.addToTrackerHits(*TrackerHit_); + }else{ + WireMeasurementDC* dcMea = + dynamic_cast(absMea); + if(dcMea){ + const edm4hep::TrackerHit* TrackerHit_ = dcMea->getTrackerHit(); + if(nullptr==TrackerHit_) std::cout << __LINE__ << std::endl; + track.addToTrackerHits(*TrackerHit_); + + + DCHit++; + if(flag) + { + + truthMomEdep.push_back(TrackerHit_->getEDep()); + double DriftDis,fittedDoca,res = 0; + GetDocaRes((id-1),DriftDis,fittedDoca,res); + driftDis.push_back(10*DriftDis); + FittedDoca.push_back(10*fittedDoca); + truthDoca.push_back(40*1e-3*TrackerHit_->getTime()); + truthRes.push_back((40*1e-3*TrackerHit_->getTime())-10*fittedDoca); + Res.push_back(res); + + + dcFit++; + } + } + } + } + + nFittedDC = dcFit; + nFittedSDT = SDTHit; + + if(fitid<1e-9) return false; + + if(m_debug>0)std::cout<getFitStatus(); + const genfit::FitStatus* fitState = getFitStatus(); double ndf = fitState->getNdf(); - if(ndf>ndfCut) return false; + if(ndf>ndfCut){ + if(m_debug>0){ + std::cout<"<getChi2(); - if(chi2>chi2Cut) return false; + if(chi2>chi2Cut){ + if(m_debug>0){ + std::cout<"<getCharge(); int isFitted = fitState->isFitted(); int isConverged = fitState->isFitConverged(); int isConvergedFully = fitState->isFitConvergedFully(); - TMatrixDSym fittedCov; + + TMatrixDSym fittedCov(6);//cm, GeV TLorentzVector fittedPos; TVector3 fittedMom; int fittedState=getFittedState(fittedPos,fittedMom,fittedCov); - if(m_debug>=2)std::cout<<"fit result: get status OK? pidType " - <=2)std::cout<<" fitting failed"<0)std::cout<ndfCut)){ + if(m_debug>0)std::cout<=2)std::cout<<"fit result: Pos("<< - fittedPos.X()<<" "<< - fittedPos.Y()<<" "<< - fittedPos.Z()<<") mom("<< - fittedMom.X()<<" "<< - fittedMom.Y()<<" "<< - fittedMom.Z()<<") p_tot "<< - fittedMom.Mag()<<" pt "<< - fittedMom.Perp()<< - " chi2 "<getBz(fittedPos.Vect())); - - - /////track status at POCA to origin - //TVector3 origin(pivot); - //TVector3 pocaToOrigin_pos(1e9*dd4hep::cm,1e9*dd4hep::cm,1e9*dd4hep::cm); - //TVector3 pocaToOrigin_mom(1e9*dd4hep::GeV,1e9*dd4hep::GeV,1e9*dd4hep::GeV); - - //if(ExtrapolateToPoint(pocaToOrigin_pos,pocaToOrigin_mom, - // m_track,origin) > 1e6*dd4hep::cm){ - // log<<"extrapolate to origin failed"; - // return false; - //} - - - //new TrackState - edm4hep::TrackState* trackState = new edm4hep::TrackState(); - trackState->location=0;//edm4hep::AtIP; - trackState->D0=helix.getD0(); - trackState->phi=helix.getPhi0(); - trackState->omega=helix.getOmega(); - trackState->Z0=helix.getZ0(); - trackState->tanLambda=helix.getTanLambda(); - trackState->referencePoint=helix.getReferencePoint(); - - // std::array covMatrix; - // int k=0; - // for(int i=0;i<5;i++){ - // for(int j=0;j<5;j++){ - // if(i<=j) covMatrix[k]=;//FIXME - // } - // } - // trackState.covMatrix= - - //new Track - edm4hep::MutableTrack* track = new edm4hep::MutableTrack(); - //track->setType(); - track->setChi2(fitState->getChi2()); - track->setNdf(fitState->getNdf()); - //track->setDEdx(); - //track->setRadiusOfInnermostHit();//FIXME - //track->addToTrackerHits(); - - //new ReconstructedParticle - //recParticle->setType(); - //dcRecParticle->setEnergy(); - - edm4hep::Vector3f momVec3(helix.getMomentum()[0], - helix.getMomentum()[1],helix.getMomentum()[2]); - recParticle.setMomentum(momVec3); - //recParticle->setReferencePoint(referencePoint); - recParticle.setCharge(helix.getCharge()); - // recParticle->setMass(); - // recParticle->setCovMatrix(); - // rcRecParticle->setStartVertex(); - //recParticle->addToTracks(track); + if(m_debug>0){ + const TLorentzVector seedPos=getSeedStatePos(); + const TVector3 seedMom=getSeedStateMom(); + std::cout< 1e6*dd4hep::cm){ + if(m_debug>0)std::cout<mm + pocaToOrigin_pos.SetXYZ(pocaToOrigin_pos.X()/dd4hep::mm, + pocaToOrigin_pos.Y()/dd4hep::mm, + pocaToOrigin_pos.Z()/dd4hep::mm); + pocaToOrigin_mom.SetXYZ(pocaToOrigin_mom.X(), + pocaToOrigin_mom.Y(), + pocaToOrigin_mom.Z()); + + //unit conversion of error matrix + TMatrixDSym covMatrix_6=fittedCov; + for(int i=0;i<5;i++){ + covMatrix_6[0][i]=fittedCov[0][i]/dd4hep::mm;//d0 column + covMatrix_6[2][i]=fittedCov[1][i]/dd4hep::mm;//omega column + covMatrix_6[3][i]=fittedCov[2][i]/dd4hep::mm;//z0 column + covMatrix_6[i][0]=fittedCov[i][0]/dd4hep::mm;//d0 row + covMatrix_6[i][2]=fittedCov[i][1]/dd4hep::mm;//omega row + covMatrix_6[i][3]=fittedCov[i][2]/dd4hep::mm;//z0 row + } + + double Bz=m_genfitField->getBz(referencePoint)/GenfitUnit::tesla; + + if(m_debug>0){ + std::cout<2){ + std::cout<2){ + std::cout<getTrackRep(id); +} + +void GenfitTrack::getSeedCov(TMatrixDSym& cov){ + cov.Zero(); + for(int i=0;i<3;i++) { + double posResolusion=1.; + cov(i,i)=posResolusion*posResolusion; //seed position + double momResolusion=5.; + cov(i+3,i+3)=momResolusion*momResolusion; //seed momentum + } +} + +//unit conversion here +void GenfitTrack::getEndPointsOfWire(int cellID,TVector3& end0,TVector3& end1) const{ + m_gridDriftChamber->cellposition(cellID,end0,end1);//dd4hep unit + end0*=(1./dd4hep::cm*GenfitUnit::cm); + end1*=(1./dd4hep::cm*GenfitUnit::cm); +} +//unit conversion here +void GenfitTrack::getTrackFromMCPartile(const edm4hep::MCParticle mcParticle, + TVectorD& trackParam, TMatrixDSym& cov) const{ + TVector3 pos; + TVector3 mom; + getPosMomFromMCPartile(mcParticle,pos,mom); + trackParam[0]=pos[0]; + trackParam[1]=pos[1]; + trackParam[2]=pos[2]; + trackParam[3]=mom[0]; + trackParam[4]=mom[1]; + trackParam[5]=mom[2]; + //cov TODO +} +//unit conversion here +void GenfitTrack::getPosMomFromMCPartile(const edm4hep::MCParticle mcParticle, + TVector3& pos,TVector3& mom) const{ + const edm4hep::Vector3d mcParticleVertex=mcParticle.getVertex();//mm + const edm4hep::Vector3f mcParticleMom=mcParticle.getMomentum();//GeV + pos[0]=mcParticleVertex.x*GenfitUnit::mm; + pos[1]=mcParticleVertex.y*GenfitUnit::mm; + pos[2]=mcParticleVertex.z*GenfitUnit::mm; + mom[0]=mcParticleMom.x*GenfitUnit::GeV; + mom[1]=mcParticleMom.y*GenfitUnit::GeV; + mom[2]=mcParticleMom.z*GenfitUnit::GeV; + if(m_debug>2){ + std::cout<<"getPosMomFromTrackState pos("<getBz(TVector3{0.,0.,0.})/GenfitUnit::tesla; + double charge_double; + CEPC::getPosMomFromTrackState(edm4HepTrack.getTrackStates(0),Bz,pos,mom,charge_double,cov); + + charge=(int) charge_double; + trackParam[0]=pos[0]*GenfitUnit::mm; + trackParam[1]=pos[1]*GenfitUnit::mm; + trackParam[2]=pos[2]*GenfitUnit::mm; + trackParam[3]=mom[0]*GenfitUnit::GeV; + trackParam[4]=mom[1]*GenfitUnit::GeV; + trackParam[5]=mom[2]*GenfitUnit::GeV; + + //cov unit conversion + for(int i=0;i<6;i++){ + for(int j=0;j<6;j++){ + cov(i,j) = cov(i,j)*GenfitUnit::mm; + } + } +} + +void GenfitTrack::getTrackFromEDMTrackFinding(const edm4hep::Track& edm4HepTrack, + double& charge, TVectorD& trackParam, TMatrixDSym& cov,TVector3& pos, + TVector3& mom){ + + // FIXME + double Bz=3*GenfitUnit::tesla; + double charge_double; + CEPC::getPosMomFromTrackState(edm4HepTrack.getTrackStates(1),Bz,pos,mom,charge_double,cov); + + charge=(int) charge_double; + trackParam[0]=pos[0]*GenfitUnit::mm; + trackParam[1]=pos[1]*GenfitUnit::mm; + trackParam[2]=pos[2]*GenfitUnit::mm; + trackParam[3]=mom[0]*GenfitUnit::GeV; + trackParam[4]=mom[1]*GenfitUnit::GeV; + trackParam[5]=mom[2]*GenfitUnit::GeV; + + //cov unit conversion + for(int i=0;i<6;i++){ + for(int j=0;j<6;j++){ + cov(i,j) = cov(i,j)*GenfitUnit::mm; + } + } +} + +//to genfit unit +void GenfitTrack::getISurfaceOUV(const dd4hep::rec::ISurface* iSurface,TVector3& o, + TVector3& u,TVector3& v, double& lengthU,double& lengthV){ + o.SetXYZ(iSurface->origin().x()/dd4hep::mm*GenfitUnit::mm, + iSurface->origin().y()/dd4hep::mm*GenfitUnit::mm, + iSurface->origin().z()/dd4hep::mm*GenfitUnit::mm); + u.SetXYZ(iSurface->u().x()/dd4hep::mm*GenfitUnit::mm, + iSurface->u().y()/dd4hep::mm*GenfitUnit::mm, + iSurface->u().z()/dd4hep::mm*GenfitUnit::mm); + v.SetXYZ(iSurface->v().x()/dd4hep::mm*GenfitUnit::mm, + iSurface->v().y()/dd4hep::mm*GenfitUnit::mm, + iSurface->v().z()/dd4hep::mm*GenfitUnit::mm); + + lengthU=iSurface->length_along_u(); + lengthV=iSurface->length_along_v(); +} + +void GenfitTrack::getMeasurementAndCov(edm4hep::TrackerHit* hit,TVector3& pos,TMatrixDSym& cov){ + + pos.SetXYZ(hit->getPosition().x*GenfitUnit::mm, + hit->getPosition().y*GenfitUnit::mm, + hit->getPosition().z*GenfitUnit::mm); +} + + +int GenfitTrack::getSigmas(int cellID,std::vector sigmaUVec, + std::vector sigmaVVec,float& sigmaU,float& sigmaV)const{ + int detTypeID=getDetTypeID(cellID); + int sigmaUID=0; + int sigmaVID=0; + int layer=0; + dd4hep::DDSegmentation::BitFieldCoder* decoder; + if(m_geomSvc->lcdd()->constant("DetID_DC")==detTypeID){ + sigmaUID=0; + sigmaVID=0; + }else if(detTypeID==lcio::ILDDetID::VXD){ + decoder=m_geomSvc->getDecoder("VXDCollection");//FIXME + layer=decoder->get(cellID,"layer");//FIXME + sigmaUID=layer+1; + sigmaVID=layer+1; + }else if(detTypeID==lcio::ILDDetID::SIT){ + sigmaUID=7; + sigmaVID=7; + }else if(detTypeID==lcio::ILDDetID::SET){ + sigmaUID=8; + sigmaVID=8; + }else if(detTypeID==lcio::ILDDetID::FTD){ + decoder=m_geomSvc->getDecoder("FTDCollection");//FIXME + layer=decoder->get(cellID,"layer");//FIXME + sigmaUID=layer+9; + sigmaVID=layer+9; + }else{ + if(m_debug>=0) std::cout<lcdd()->constant("DetID_DC")== + getDetTypeID(hit->getCellID()); +} + +// TrackFinding +void GenfitTrack::getSortedTrackerHitsTrF( + std::vector trackerHits, + std::vector& sortedDCTrackerHits, + int sortMethod){ + + std::vector > sortedDCTrackerHitPair; + for(auto trackerHit : trackerHits){ + + double time=trackerHit->getTime(); + if(0==sortMethod){ + //by time + sortedDCTrackerHitPair.push_back(std::make_pair(time,trackerHit)); + if(m_debug>0){ std::cout<<"sorted DC digi by time"<get(trackerHit->getCellID(),"layer"),trackerHit)); + if(m_debug>0){ std::cout<<"sorted DC digi by layer"<0){ + std::cout<<"trackerHits on track after sort\n"; + for(auto trackerHit:sortedDCTrackerHits){ + std::cout<getCellID()<<"("<get(trackerHit->getCellID(),"layer") + <<","<get(trackerHit->getCellID(),"cellID")<<")\n"; + } + std::cout<<"\n"; std::cout.unsetf(std::ios_base::floatfield); + } + return; +}//end of getSortedTrackerHits + +void GenfitTrack::getSortedTrackerHits( + std::vector& trackerHits, + const edm4hep::MCRecoTrackerAssociationCollection* assoHits, + std::vector& sortedDCTrackerHits, + int sortMethod){ + + std::vector > sortedDCTrackerHitPair; + for(auto trackerHit : trackerHits){ + edm4hep::TrackerHit* thisHit = trackerHit; + if(!isCDCHit(thisHit))continue;//skip non-DC trackerHit + + edm4hep::SimTrackerHit simTrackerHitAsso; + getAssoSimTrackerHit(assoHits,trackerHit,simTrackerHitAsso); + double time=simTrackerHitAsso.getTime(); + if(0==sortMethod){ + //by time + sortedDCTrackerHitPair.push_back(std::make_pair(time,trackerHit)); + if(m_debug>0){ std::cout<<"sorted DC digi by time"<get(trackerHit->getCellID(),"layer"),trackerHit)); + if(m_debug>0){ std::cout<<"sorted DC digi by layer"<0){ + std::cout<<"trackerHits on track after sort\n"; + for(auto trackerHit:sortedDCTrackerHits){ + std::cout<getCellID()<<"("<get(trackerHit->getCellID(),"layer") + <<","<get(trackerHit->getCellID(),"cellID")<<")\n"; + } + std::cout<<"\n"; std::cout.unsetf(std::ios_base::floatfield); + } + return; +}//end of getSortedTrackerHits + +//make a genfit hit and do hits selection +GenfitHit* GenfitTrack::makeAGenfitHit(edm4hep::TrackerHit* trackerHit, + edm4hep::SimTrackerHit* simTrackerHitAsso, + double sigma,bool truthAmbig,double skipCorner,double skipNear){ + + //TODO truthAmbig + double driftVelocity=40.;//FIXME + GenfitHit* genfitHit=new GenfitHit(trackerHit,simTrackerHitAsso,m_decoderDC, + m_gridDriftChamber,driftVelocity,sigma*GenfitUnit::mm); + //skip corner hit + if(fabs(genfitHit->getDriftDistance())>skipCorner){ + if(m_debug) std::cout<<" skip hit dd > skipCorner"<getDriftDistance() class TLorentzVector; +class IGeomSvc; +class WireMeasurementDC; namespace genfit{ class Track; @@ -44,106 +50,133 @@ namespace edm4hep{ class MutableReconstructedParticle; class MCRecoTrackerAssociationCollection; class Track; - class Track; + class MutableTrack; + class TrackCollection; class TrackerHit; + class SimTrackerHit; class Vector3d; class Vector3f; + class TrackerHitCollection; } namespace dd4hep { namespace DDSegmentation{ class GridDriftChamber; + class BitFieldCoder; + } + namespace rec{ + class ISurface; } } class GenfitTrack { friend int GenfitFitter::processTrack( GenfitTrack* track, bool resort); - + friend int GenfitFitter::processTrackWithRep( GenfitTrack* track, int repID, bool resort); - friend double GenfitFitter::extrapolateToHit(TVector3& poca, - TVector3& pocaDir, - TVector3& pocaOnWire, double& doca, const GenfitTrack* track, - TVector3 pos, TVector3 mom, TVector3 end0, TVector3 end1, int debug, - int repID, bool stopAtBoundary, bool calcJacobianNoise); - - friend double GenfitFitter::extrapolateToCylinder(TVector3& pos, - TVector3& mom, - GenfitTrack* track, double radius, const TVector3 linePoint, - const TVector3 lineDirection, int hitID, int repID, - bool stopAtBoundary, bool calcJacobianNoise); - - friend double GenfitFitter::extrapolateToPoint(TVector3& pos, TVector3& mom, - const GenfitTrack* genfitTrack, const TVector3& point, int repID, - bool stopAtBoundary, bool calcJacobianNoise) const; - public: GenfitTrack(const GenfitField* field, const dd4hep::DDSegmentation::GridDriftChamber* seg, + SmartIF geom, const char* name="GenfitTrack"); virtual ~GenfitTrack(); - /// Add a Genfit track - virtual bool createGenfitTrack(int pdgType,int charge,TLorentzVector pos, TVector3 mom, - TMatrixDSym covM); - //virtual bool createGenfitTrack(TLorentzVector posInit,TVector3 momInit, - //TMatrixDSym covMInit); - ///Create genfit track from MCParticle bool createGenfitTrackFromMCParticle(int pidTyep,const edm4hep::MCParticle& mcParticle, double eventStartTime=0.); - bool createGenfitTrackFromEDM4HepTrack(int pidType, edm4hep::Track track, - double eventStartTime); - - // /// Prepare a hit list, return number of hits on track - // int PrepareHits();//TODO - - /// Add a space point measurement, return number of hits on track - bool addSpacePointTrakerHit(edm4hep::TrackerHit hit, int hitID); - - /// Add a space point measurement, return number of hits on track - virtual bool addSpacePointMeasurement(const TVectorD&, double, - int detID=-1, int hitID=-1, bool smear=false); - - /// Add a WireMeasurement with MC truth position smeared by sigma - virtual void addWireMeasurement(double driftDistance, - double driftDistanceError, const TVector3& endPoint1, - const TVector3& endPoint2, int lrAmbig, int detID, int hitID); - - /// Add a WireMeasurement with DC digi - virtual bool addWireMeasurementOnTrack(edm4hep::Track track, double sigma); - - ///Add space point from truth to track - int addSimTrackerHits( edm4hep::Track track, - const edm4hep::MCRecoTrackerAssociationCollection* assoHits, - float sigma,bool smear=false);// float nSigmaSelection + bool createGenfitTrackFromEDM4HepTrack(int pidType,const edm4hep::Track& track, + double eventStartTime,bool isUseCovTrack); + + /// ---------Add measurements--------- + ///Add one space point measurement, return number of hits on track + virtual bool addSpacePointMeasurement(const TVector3&,std::vector + sigmaU,std::vector sigmaV,int cellID,int hitID); + + ///Add silicon space points from edm4hep::track + int addSpacePointsSi(const edm4hep::Track& track, + std::vector sigmaU,std::vector sigmaV); + + ///Add drift chamber space points from edm4hep::track + int addSpacePointsDC(const edm4hep::Track& track, + const edm4hep::MCRecoTrackerAssociationCollection* assoHits, + std::vector sigmaU,std::vector sigmaV); + + ///Add WireMeasurements of hits on track + virtual int addWireMeasurementsOnTrack(edm4hep::Track& track,float sigma, + const edm4hep::MCRecoTrackerAssociationCollection* assoHits, + int sortMethod,bool truthAmbig,float skipCorner, float skipNear); + + ///Add WireMeasurements of hits on track from hit selection + virtual int addWireMeasurementsFromList(std::vector& hits,float sigma, + const edm4hep::MCRecoTrackerAssociationCollection* assoHits, + int sortMethod,bool truthAmbig,float skipCorner, float skipNear); + + virtual int addWireMeasurementsFromListTrF(const edm4hep::TrackerHitCollection* trkHits, + float sigma,int sortMethod); + + ///Add one silicon hits + bool addSiliconMeasurement(edm4hep::TrackerHit* hit, + float sigmaU,float sigmaV,int cellID,int hitID); + + ///Add silicon measurements, return number of hits on track + int addSiliconMeasurements(edm4hep::Track& track, + std::vector sigmaU,std::vector sigmaV); + + bool debugDistance(const edm4hep::TrackerHitCollection* dCDigiCol, + int& nFittedDC,int& nFittedSDT,int& ngenfitHit, + std::vector& smearDistance, + std::vector& truthDistance,double driftVelocity); + bool GetDocaRes(int hitID, double& DriftDis,double& fittedDoca, + double& res,int repID=0, bool biased=true) const; ///Store track to ReconstructedParticle - bool storeTrack(edm4hep::MutableReconstructedParticle& dcRecParticle,int pidType, - int ndfCut=1e9, double chi2Cut=1.e9); + bool storeTrack(edm4hep::MutableReconstructedParticle& recParticle, + edm4hep::MutableTrack& track, + TVector3& pocaToOrigin_pos, + TVector3& pocaToOrigin_mom, + TMatrixDSym& pocaToOrigin_cov, + int pidType, int ndfCut, double chi2Cut, + int& nFittedDC, int& nFittedSDT,int& ngenfitHit, + std::vector& trackL, std::vector& hitMom, + std::vector& truthMomEdep, + const edm4hep::MCRecoTrackerAssociationCollection* assoHits, + std::vector& driftDis, + std::vector& FittedDoca, + std::vector& truthDoca, + std::vector& Res, + std::vector& truthRes); ///A tool to convert track to the first layer of DC - template - void pivotToFirstLayer(const edm4hep::Vector3d& pos, const MomT& mom, - edm4hep::Vector3d& firstPos, MomT& firstMom) { - //FIXME, TODO - firstPos=pos; - firstMom=mom; - } + + void pivotToFirstLayer(const edm4hep::Vector3d& pos, + const edm4hep::Vector3f& mom, edm4hep::Vector3d& firstPos, + edm4hep::Vector3f& firstMom); /// Copy a track to event //void CopyATrack()const; - ///Extrapolate to Hit - double extrapolateToHit( TVector3& poca, TVector3& pocaDir, - TVector3& pocaOnWire, double& doca, TVector3 pos, TVector3 mom, - TVector3 end0,//one end of the hit wire - TVector3 end1,//the orhter end of the hit wire - int debug, - int repID, - bool stopAtBoundary, - bool calcJacobianNoise) const; + //return dd4hep unit + double extrapolateToHit(TVector3& poca, TVector3& pocaDir, + TVector3& pocaOnWire, double& doca,edm4hep::MCParticle mcParticle, + int cellID, int repID, bool stopAtBoundary, bool calcJacobianNoise)const; + /// Extrapolate the track to the point + /// Output: pos and mom of POCA point to point + /// Input: genfitTrack,point,repID,stopAtBoundary and calcAverageState + /// repID same with pidType + double extrapolateToPoint(TVector3& pos, TVector3& mom, TMatrixDSym& cov, + const TVector3& point, int repID=0, + bool stopAtBoundary = false, bool calcJacobianNoise = true) const; + + /// Extrapolate the track to the cyliner at fixed raidus + /// Output: pos and mom at fixed radius + /// Input: genfitTrack, radius of cylinder at center of the origin, + /// repID, stopAtBoundary and calcAverageState + double extrapolateToCylinder(TVector3& pos, TVector3& mom, + double radius, const TVector3 linePoint, + const TVector3 lineDirection, int hitID =0, int repID=0, + bool stopAtBoundary=false, bool calcJacobianNoise=true); + bool getPosMomCovMOP(int hitID, TLorentzVector& pos, TVector3& mom, TMatrixDSym& cov, int repID=0) const; @@ -153,11 +186,12 @@ class GenfitTrack { const TVector3 getSeedStateMom() const; void getSeedStateMom(TLorentzVector& pos, TVector3& mom) const; unsigned int getNumPoints() const; + unsigned int getNumPointsDet(int cellID) const; /// get the fitted track status const genfit::FitStatus* getFitStatus(int repID=0) const; int getFittedState(TLorentzVector& pos, TVector3& mom, TMatrixDSym& cov, - int repID=0, bool biased=true) const; + int trackPointId=0, int repID=0, bool biased=true) const; int getNumPointsWithFittedInfo(int repID=0) const; bool getFirstPointWithFittedInfo(int repID=0) const; bool fitSuccess(int repID) const; @@ -174,7 +208,9 @@ class GenfitTrack { void print(TLorentzVector pos, TVector3 mom, const char* comment="") const; /// set and get debug level - void setDebug(int debug); + void setDebug(int debug){m_debug=debug;} + void setDebugGenfit(int debug); + void setDebugLocal(int debug); int getDebug(void) const { return m_debug;} /// get name of this object @@ -182,30 +218,74 @@ class GenfitTrack { genfit::Track* getTrack() const{return m_track;} /// Add a track representation - genfit::RKTrackRep* addTrackRep(int pdg); - + genfit::RKTrackRep* addTrackRep(int pdgType,int charge); + + /// Get a hit according to index + GenfitHit* GetHit(long unsigned int i) const; + + static void getTrackFromEDMTrackFinding(const edm4hep::Track& edm4HepTrack, + double& charge, TVectorD& trackParam, TMatrixDSym& cov,TVector3& pos, + TVector3& mom); protected: //genfit::Track* getTrack() {return m_track;} private: + /// ---------Add a Genfit track------- + bool createGenfitTrack(int pdgType,int charge, + TVectorD trackParam, TMatrixDSym covMInit_6); + + int getDetTypeID(unsigned long long cellID) const; const char* m_name; ///Note! private functions are using genfit unit, cm and MeV - genfit::AbsTrackRep* getRep(int id=0) const {return m_reps[id];} + genfit::AbsTrackRep* getRep(int id=0) const; bool getMOP(int hitID, genfit::MeasuredStateOnPlane& mop, genfit::AbsTrackRep* trackRep=nullptr) const; + const dd4hep::rec::ISurface* getISurface(edm4hep::TrackerHit* hit); + void getSeedCov(TMatrixDSym& cov); + void getAssoSimTrackerHit( + const edm4hep::MCRecoTrackerAssociationCollection*& assoHits, + edm4hep::TrackerHit* trackerHit, + edm4hep::SimTrackerHit& simTrackerHit) const; + void getEndPointsOfWire(int cellID,TVector3& end0,TVector3& end1)const; + void getTrackFromEDMTrack(const edm4hep::Track& edm4HepTrack, + double& charge, TVectorD& trackParam, TMatrixDSym& cov) const; + void getTrackFromMCPartile(const edm4hep::MCParticle mcParticle, + TVectorD& trackParam, TMatrixDSym& cov) const; + void getPosMomFromMCPartile(const edm4hep::MCParticle mcParticle, + TVector3& pos,TVector3& mom) const; + void clearGenfitHitVec(); + void getISurfaceOUV(const dd4hep::rec::ISurface* iSurface,TVector3& o, + TVector3& u,TVector3& v,double& lengthU,double& lengthV); + void getMeasurementAndCov(edm4hep::TrackerHit* hit,TVector3& pos,TMatrixDSym& cov); + int getSigmas(int cellID,std::vector sigmaUVec, + std::vector sigmaVVec,float& sigmaU,float& sigmaV)const; + bool isCDCHit(edm4hep::TrackerHit* hit); + GenfitHit* makeAGenfitHit(edm4hep::TrackerHit* trackerHit, + edm4hep::SimTrackerHit* simTrackerHitAsso, + double sigma,bool truthAmbig,double skipCorner,double skipNear); + void getSortedTrackerHits(std::vector& hits, + const edm4hep::MCRecoTrackerAssociationCollection* assoHits, + std::vector& sortedDCTrackerHits, + int sortMethod); + void getSortedTrackerHitsTrF(std::vector trackerHits, + std::vector& sortedDCTrackerHits, + int sortMethod=1); genfit::Track* m_track;/// track - std::vector m_reps;/// track repesentations int m_debug;/// debug level + int m_debugLocal;/// debug level local + SmartIF m_geomSvc; const GenfitField* m_genfitField;//pointer to genfit field const dd4hep::DDSegmentation::GridDriftChamber* m_gridDriftChamber; + const dd4hep::DDSegmentation::BitFieldCoder* m_decoderDC; static const int s_PDG[2][5]; -}; + std::vector m_genfitHitVec; +}; #endif diff --git a/Reconstruction/RecGenfitAlg/src/GenfitUnit.h b/Reconstruction/RecGenfitAlg/src/GenfitUnit.h new file mode 100644 index 000000000..5601fcc28 --- /dev/null +++ b/Reconstruction/RecGenfitAlg/src/GenfitUnit.h @@ -0,0 +1,13 @@ +#ifndef GENFITUNIT_H +#define GENFITUNIT_H + +namespace GenfitUnit{ + static constexpr double cm=1.; + static constexpr double mm=0.1; + static constexpr double um=mm*1e-3; + static constexpr double GeV=1.; + static constexpr double kilogauss=1.; + static constexpr double tesla = 10.*kilogauss; +} + +#endif diff --git a/Reconstruction/RecGenfitAlg/src/LSFitting.cpp b/Reconstruction/RecGenfitAlg/src/LSFitting.cpp new file mode 100644 index 000000000..6763e4233 --- /dev/null +++ b/Reconstruction/RecGenfitAlg/src/LSFitting.cpp @@ -0,0 +1,92 @@ +#include "LSFitting.h" +#include +#include + + +std::vector LSFitting::m_wireX; +std::vector LSFitting::m_wireY; +std::vector LSFitting::m_driftDist; + +void LSFitting::Fitting(double& fitXc, double& fitYc, double& fitRadius) +{ + TMinuit* tMinuit = new TMinuit(3); //initialize TMinuit with a maximum of 3 params + tMinuit->SetFCN(FCN); + double arglist[2]; + arglist[0] = 0; + int ierflg = 0; + tMinuit->SetPrintLevel(-1); // no print + tMinuit->mnexcm("SET NOW", arglist,1,ierflg); // no warning + arglist[0] = 1; + tMinuit->mnexcm("SET ERR", arglist ,1,ierflg); + double xcErr= 1; + double ycErr= 1; + double radiusErr= 5; + //fitXc=0; + //fitYc=0; + //fitRadius=1000; + //std::cout<<" test Fitting---------------"<mnparm(0,"xc",fitXc,xcErr/1.e4,fitXc-xcErr,fitXc+xcErr,ierflg); + tMinuit->mnparm(1,"yc",fitYc,ycErr/1.e4,fitYc-ycErr,fitYc+ycErr,ierflg); + tMinuit->mnparm(2,"r",fitRadius,radiusErr/1.e4,fitRadius-radiusErr,fitRadius+radiusErr,ierflg); + arglist[0] = 0.0; + arglist[1] = 1.0; + tMinuit->mnexcm("MIGRAD", arglist ,2,ierflg); + double temp; + tMinuit->GetParameter(0, fitXc,temp); + tMinuit->GetParameter(1, fitYc,temp); + tMinuit->GetParameter(2, fitRadius,temp); + //std::cout<<" test Fitting---------------"< tmp1; + m_wireX.swap(tmp1); + std::vector tmp2; + m_wireY.swap(tmp2); + std::vector tmp3; + m_driftDist.swap(tmp3); +} + +void LSFitting::setWire(double x,double y){ + m_wireX.push_back(x); + m_wireY.push_back(y); +} +void LSFitting::setDrift(double driftDist){ + m_driftDist.push_back(driftDist); +} + +void LSFitting::print(){ + std::cout<<" nHit "< +#include + +class LSFitting{ +public: + void Fitting(double& fitXc, double& fitYc, double& fitRadisu); + static void FCN(int &npar, double *gin, double &f, double *par, int iflag); + void setWire(double x,double y); + void setDrift(double driftDist); + void print(); + void clear(); + static std::vector m_wireX; + static std::vector m_wireY; + static std::vector m_driftDist; +}; +#endif diff --git a/Reconstruction/RecGenfitAlg/src/PlanarMeasurementSDT.cpp b/Reconstruction/RecGenfitAlg/src/PlanarMeasurementSDT.cpp new file mode 100644 index 000000000..cfa03ce9d --- /dev/null +++ b/Reconstruction/RecGenfitAlg/src/PlanarMeasurementSDT.cpp @@ -0,0 +1,127 @@ +/* Copyright 2008-2010, Technische Universitaet Muenchen, + Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch + + This file is part of GENFIT. + + GENFIT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published + by the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + GENFIT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with GENFIT. If not, see . +*/ + +#include "PlanarMeasurementSDT.h" + +#include +#include +#include +#include +#include + +#include +#include + +namespace genfit { + +PlanarMeasurementSDT::PlanarMeasurementSDT(int nDim) + : AbsMeasurement(nDim), physicalPlane_(), planeId_(-1), stripV_(false) +{ + assert(nDim >= 1); +} + +PlanarMeasurementSDT::PlanarMeasurementSDT(const TVectorD& rawHitCoords, const TMatrixDSym& rawHitCov, int detId, int hitId, TrackPoint* trackPoint) + : AbsMeasurement(rawHitCoords, rawHitCov, detId, hitId, trackPoint), physicalPlane_(), planeId_(-1), stripV_(false) +{ + assert(rawHitCoords_.GetNrows() >= 1); +} + + +SharedPlanePtr PlanarMeasurementSDT::constructPlane(const StateOnPlane&) const { + if (!physicalPlane_) { + Exception exc("PlanarMeasurementSDT::constructPlane(): No plane has been set!", __LINE__,__FILE__); + throw exc; + } + return physicalPlane_; +} + + +std::vector PlanarMeasurementSDT::constructMeasurementsOnPlane(const StateOnPlane& state) const { + + MeasurementOnPlane* mop = new MeasurementOnPlane(rawHitCoords_, + rawHitCov_, + state.getPlane(), state.getRep(), constructHMatrix(state.getRep())); + + std::vector retVal; + retVal.push_back(mop); + return retVal; +} + + +const AbsHMatrix* PlanarMeasurementSDT::constructHMatrix(const AbsTrackRep* rep) const { + + if (dynamic_cast(rep) == nullptr) { + Exception exc("SpacepointMeasurement default implementation can only handle state vectors of type RKTrackRep!", __LINE__,__FILE__); + throw exc; + } + + switch(rawHitCoords_.GetNrows()) { + case 1: + if (stripV_) + return new HMatrixV(); + return new HMatrixU(); + + case 2: + return new HMatrixUV(); + + default: + Exception exc("PlanarMeasurementSDT default implementation can only handle 1D (strip) or 2D (pixel) measurements!", __LINE__,__FILE__); + throw exc; + } + +} + +void PlanarMeasurementSDT::Streamer(TBuffer &R__b) +{ + // Stream an object of class genfit::PlanarMeasurementSDT. + + //This works around a msvc bug and should be harmless on other platforms + typedef ::genfit::PlanarMeasurementSDT thisClass; + UInt_t R__s, R__c; + if (R__b.IsReading()) { + Version_t R__v = R__b.ReadVersion(&R__s, &R__c); if (R__v) { } + //This works around a msvc bug and should be harmless on other platforms + typedef genfit::AbsMeasurement baseClass0; + baseClass0::Streamer(R__b); + char flag; + R__b >> flag; + physicalPlane_.reset(); + if (flag) { + physicalPlane_.reset(new DetPlane()); + physicalPlane_->Streamer(R__b); + } + R__b >> planeId_; + R__b.CheckByteCount(R__s, R__c, thisClass::IsA()); + } else { + R__c = R__b.WriteVersion(thisClass::IsA(), kTRUE); + //This works around a msvc bug and should be harmless on other platforms + typedef genfit::AbsMeasurement baseClass0; + baseClass0::Streamer(R__b); + if (physicalPlane_) { + R__b << (char)1; + physicalPlane_->Streamer(R__b); + } else { + R__b << (char)0; + } + R__b << planeId_; + R__b.SetByteCount(R__c, kTRUE); + } +} + +} /* End of namespace genfit */ diff --git a/Reconstruction/RecGenfitAlg/src/PlanarMeasurementSDT.h b/Reconstruction/RecGenfitAlg/src/PlanarMeasurementSDT.h new file mode 100644 index 000000000..27867a9fa --- /dev/null +++ b/Reconstruction/RecGenfitAlg/src/PlanarMeasurementSDT.h @@ -0,0 +1,94 @@ +/* Copyright 2008-2010, Technische Universitaet Muenchen, + Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch + + This file is part of GENFIT. + + GENFIT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published + by the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + GENFIT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with GENFIT. If not, see . +*/ +/** @addtogroup genfit + * @{ + */ + +#ifndef genfit_PlanarMeasurementSDT_h +#define genfit_PlanarMeasurementSDT_h + +#include "AbsMeasurement.h" +#include "AbsHMatrix.h" +#include "MeasurementOnPlane.h" + +namespace edm4hep{ + class TrackerHit; + class SimTrackerHit; +} + +namespace genfit { + +class AbsTrackRep; + +/** @brief Measurement class implementing a planar hit geometry (1 or 2D). + * + * @author Christian Höppner (Technische Universität München, original author) + * @author Sebastian Neubert (Technische Universität München, original author) + * @author Johannes Rauch (Technische Universität München, original author) + * + * The main feature of this type of hit is, that the detector plane + * is defined by the detector hardware. + */ +class PlanarMeasurementSDT : public AbsMeasurement { + + public: + PlanarMeasurementSDT(int nDim = 1); + PlanarMeasurementSDT(const TVectorD& rawHitCoords, const TMatrixDSym& rawHitCov, int detId, int hitId, TrackPoint* trackPoint); + + virtual ~PlanarMeasurementSDT() {;} + + virtual AbsMeasurement* clone() const override {return new PlanarMeasurementSDT(*this);} + + int getPlaneId() const {return planeId_;} + + virtual SharedPlanePtr constructPlane(const StateOnPlane& state) const override; + + virtual std::vector constructMeasurementsOnPlane(const StateOnPlane& state) const override; + + virtual const AbsHMatrix* constructHMatrix(const AbsTrackRep*) const override; + + virtual void setPlane(const SharedPlanePtr& physicalPlane, int planeId = -1) {physicalPlane_ = physicalPlane; planeId_ = planeId;} + + /** @brief Use if the coordinate for 1D hits measured in V direction. + * + * Per default for 1D planar hits, the coordinate is measured in U direction. + * With this function you can set it to be measured in V direction. + * This affects the outcoe of constructHMatrix(). + */ + void setStripV(bool v = true) {stripV_ = v;} + + void setTrackerHit(const edm4hep::TrackerHit* hit){trackerHit_=hit;} + const edm4hep::TrackerHit* getTrackerHit(){return trackerHit_;} + + protected: + SharedPlanePtr physicalPlane_; //! This is persistent, but '!' makes ROOT shut up. + int planeId_; // planeId id is -1 per default + bool stripV_; + const edm4hep::SimTrackerHit* simTrackerHit_; + const edm4hep::TrackerHit* trackerHit_; + + public: + + ClassDefOverride(PlanarMeasurementSDT,1) +}; + +} /* End of namespace genfit */ +/** @} */ + +#endif // genfit_PlanarMeasurementSDT_h diff --git a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.cpp b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.cpp new file mode 100644 index 000000000..c15dee5e0 --- /dev/null +++ b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.cpp @@ -0,0 +1,972 @@ +#include "RecGenfitAlgSDT.h" +#include "GenfitTrack.h" +#include "GenfitFitter.h" +#include "GenfitField.h" +#include "GenfitUnit.h" + +//genfit +#include "EventDisplay.h" +#include "ReferenceStateOnPlane.h" + +//cepcsw +#include "DetInterface/IGeomSvc.h" +#include "DataHelper/HelixClass.h" +#include "DataHelper/TrackHelper.h" +#include "DataHelper/TrackerHitHelper.h" +#include "DetSegmentation/GridDriftChamber.h" +#include "UTIL/ILDConf.h" + +//externals +#include "edm4hep/EventHeaderCollection.h" +#include "edm4hep/MCParticle.h" +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/SimTrackerHitCollection.h" +#include "edm4hep/TrackerHitCollection.h" +#include "edm4hep/TrackCollection.h" +#include "edm4hep/MCRecoTrackerAssociationCollection.h" +#include "edm4hep/ReconstructedParticle.h" +#include "edm4hep/ReconstructedParticleCollection.h" +#include "edm4hep/Track.h" +#include "edm4hep/TrackCollection.h" +#include "DD4hep/DD4hepUnits.h" +#include "DD4hep/Detector.h" +#include "UTIL/BitField64.h" +#include "DDSegmentation/Segmentation.h" +#include "TRandom.h" +#include "TLorentzVector.h" + +//stl +#include +#include "time.h" +#include +#include +#include + +#include + +#include +#include + +#include "map" + +DECLARE_COMPONENT( RecGenfitAlgSDT ) + + ///////////////////////////////////////////////////////////////////// + RecGenfitAlgSDT::RecGenfitAlgSDT(const std::string& name, + ISvcLocator* pSvcLocator): + GaudiAlgorithm(name, pSvcLocator),m_nPDG(5),m_dd4hepDetector(nullptr), + m_gridDriftChamber(nullptr),m_decoder(nullptr) +{ + declareProperty("EventHeaderCollection", m_headerCol); + declareProperty("MCParticleCollection", m_mcParticleCol, + "Handle of the input MCParticle collection"); + declareProperty("DigiDCHitCollection", m_DCDigiCol, + "Handle of DC digi(TrakerHit) collection"); + declareProperty("DCHitAssociationCollection", m_DCHitAssociationCol, + "Handle of DCsimTrackerHit and DCTrackerHit association collection"); + declareProperty("SDTTrackCollection", m_SDTTrackCol, + "Handle of input silicon track collection"); + declareProperty("SDTRecTrackCollection",m_SDTRecTrackCol, + "Handle of input silicon rec. track collection"); + declareProperty("DCTrackCollection", m_dcTrackCol, + "Handle of DC track collection"); + declareProperty("SDTRecParticleCollection", m_SDTRecParticleCol, + "Handle of silicon+drift chamber rec. particle collection"); + + declareProperty("SimTrackerHitCollection",m_simVXDHitCol, + "Handle of the VXDsimTrackerHit collection"); + declareProperty("SimTrackerHitCollection",m_simSETHitCol, + "Handle of the SETsimTrackerHit collection"); + declareProperty("SimTrackerHitCollection",m_simSITHitCol, + "Handle of the SITsimTrackerHit collection"); + declareProperty("SimTrackerHitCollection",m_simFTDHitCol, + "Handle of the FTDsimTrackerHit collection"); + declareProperty("SimTrackerHitCollection",m_simDCHitCol, + "Handle of the DCsimTrackerHit collection"); + declareProperty("SimTrackerHitCollection",m_simVXDHitCol, + "Handle of the VXDsimTrackerHit collection"); + declareProperty("NoiseDCHitAssociationCollection",r_NoiseAssociationCol, + "Handle of the DCSimTrackerHits and Noise TrackerHit collection"); + + declareProperty("SmearDCHitAssociationCollection", r_SmearAssociationCol, + "Handle of output smear simulationhit and TrackerHit collection"); + + declareProperty("DCTrackFindingHitCollection", m_DCTrackFindingCol, + "Handle of output TrackFinding TrackerHit collection"); + +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * + +StatusCode RecGenfitAlgSDT::initialize() +{ + MsgStream log(msgSvc(), name()); + info()<<" RecGenfitAlgSDT initialize()"<service("GeomSvc"); + if (nullptr==m_geomSvc) { + std::cout<<"Failed to find GeomSvc"<lcdd(); + + ///Get Field + m_dd4hepField=m_geomSvc->lcdd()->field(); + + /// New a genfit fitter + m_genfitFitter=new GenfitFitter(m_fitterType.toString().c_str()); + m_genfitField=new GenfitField(m_dd4hepField); + m_genfitFitter->setField(m_genfitField); + m_genfitFitter->setGeoMaterial(m_geomSvc->lcdd(),m_extMinDistCut, + m_skipWireMaterial); + m_genfitFitter->setEnergyLossBrems(m_correctBremsstrahlung); + m_genfitFitter->setNoiseBrems(m_correctBremsstrahlung); + //genfit::eMultipleMeasurementHandling(m_multipleMeasurementHandling.value())); + if(m_debug>10) m_genfitFitter->setDebug(m_debug-10); + if(m_noMaterialEffects) m_genfitFitter->setNoEffects(true); + if(-1==m_debugPid) m_genfitFitter->setNoEffects(true); + if(-1==m_debugPid) m_debugPid=0;//charged geantino with electron pid + if(m_fitterType=="DAF"||m_fitterType=="DafRef"){ + m_genfitFitter->setMaxIterationsBetas(m_bStart,m_bFinal,m_maxIteration); + } else { + m_genfitFitter->setMaxIterations(m_maxIteration); + } + //print genfit parameters + if(m_debug) m_genfitFitter->print(); + if(""!=m_genfitHistRootName) m_genfitFitter->initHist(m_genfitHistRootName); + + //initialize member vairables + for(int i=0;i<5;i++) m_fitSuccess[i]=0; + m_nRecTrack=0; + ///Get Readout + dd4hep::Readout readout=m_dd4hepDetector->readout(m_readout_name); + ///Get Segmentation + m_gridDriftChamber=dynamic_cast + (readout.segmentation().segmentation()); + if(nullptr==m_gridDriftChamber){ + error() << "Failed to get the GridDriftChamber" << endmsg; + return StatusCode::FAILURE; + } + + m_cell_width = 0.5*(m_gridDriftChamber->cell_Size()); //dd4hep::cm + debug() << " m_cell_width = " << m_cell_width << " dd4hep::cm"<< endmsg; + m_skipCorner = m_cell_width/dd4hep::cm; + + ///Get Decoder + m_decoder = m_geomSvc->getDecoder(m_readout_name); + if (nullptr==m_decoder) { + error() << "Failed to get the decoder" << endmsg; + return StatusCode::FAILURE; + } + + + ///book tuple + NTuplePtr nt(ntupleSvc(), "RecGenfitAlgSDT/recGenfitAlgSDT"); + if(nt){ + m_tuple=nt; + }else{ + m_tuple=ntupleSvc()->book("RecGenfitAlgSDT/recGenfitAlgSDT", + CLID_ColumnWiseTuple,"RecGenfitAlgSDT"); + if(m_tuple){ + StatusCode sc; + sc=m_tuple->addItem("run",m_run); + sc=m_tuple->addItem("evt",m_evt); + sc=m_tuple->addItem("tkId",m_tkId); + sc=m_tuple->addItem("nStdTrack",m_nSdtTrack); + sc=m_tuple->addItem("nStdTrackHit",m_nSdtTrackHit,0,1000); + + sc=m_tuple->addItem("nSdtRecTrack",m_nSdtRecTrack); + + + sc=m_tuple->addItem("mcIndex",m_mcIndex,0,100);//max. 100 particles + sc=m_tuple->addItem("seedMomP",m_mcIndex,m_seedMomP);//for some track debug + sc=m_tuple->addItem("seedMomPt",m_mcIndex,m_seedMomPt); + sc=m_tuple->addItem("seedMomQ",m_mcIndex,m_seedMomQ); + sc=m_tuple->addItem("seedPos",m_mcIndex,m_seedPos,3); + sc=m_tuple->addItem("seedMom",m_mcIndex,m_seedMom,3); + sc=m_tuple->addItem("truthPocaMc",m_mcIndex,m_truthPocaMc,3); + sc=m_tuple->addItem("pocaPosMc",m_mcIndex,m_pocaPosMc,3); + sc=m_tuple->addItem("pocaMomMc",m_mcIndex,m_pocaMomMc,3); + sc=m_tuple->addItem("pocaMomMcP",m_mcIndex,m_pocaMomMcP); + sc=m_tuple->addItem("pocaMomMcPt",m_mcIndex,m_pocaMomMcPt); + sc=m_tuple->addItem("pocaPosMdc",m_mcIndex,m_pocaPosMdc,3); + sc=m_tuple->addItem("pocaMomMdc",m_mcIndex,m_pocaMomMdc,3); + sc=m_tuple->addItem("index",m_pidIndex, 0, 5); + //sc=m_tuple->addItem("firstPosKalP",5,3,m_firstPosKal); + //sc=m_tuple->addItem("firstMomKalP",5,m_firstMomKalP); + //sc=m_tuple->addItem("firstMomKalPt",5,m_firstMomKalPt); + + sc=m_tuple->addItem("ErrorcovMatrix6",m_mcIndex,m_ErrorcovMatrix6,6); + sc=m_tuple->addItem("McErrCov",m_mcIndex,m_McErrCov,6); + sc=m_tuple->addItem("posx",m_mcIndex,m_posx); + sc=m_tuple->addItem("posy",m_mcIndex,m_posy); + sc=m_tuple->addItem("posz",m_mcIndex,m_posz); + + sc=m_tuple->addItem("momx",m_mcIndex,m_momx); + sc=m_tuple->addItem("momy",m_mcIndex,m_momy); + sc=m_tuple->addItem("momz",m_mcIndex,m_momz); + + sc=m_tuple->addItem("PosMcX",m_mcIndex,m_PosMcX); + sc=m_tuple->addItem("PosMcY",m_mcIndex,m_PosMcY); + sc=m_tuple->addItem("PosMcZ",m_mcIndex,m_PosMcZ); + + sc=m_tuple->addItem("MomMcX",m_mcIndex,m_MomMcX); + sc=m_tuple->addItem("MomMcY",m_mcIndex,m_MomMcY); + sc=m_tuple->addItem("MomMcZ",m_mcIndex,m_MomMcZ); + + sc=m_tuple->addItem("PocaPosX",m_mcIndex,m_PocaPosX); + sc=m_tuple->addItem("PocaPosY",m_mcIndex,m_PocaPosY); + sc=m_tuple->addItem("PocaPosZ",m_mcIndex,m_PocaPosZ); + + sc=m_tuple->addItem("PocaMomX",m_mcIndex,m_PocaMomX); + sc=m_tuple->addItem("PocaMomY",m_mcIndex,m_PocaMomY); + sc=m_tuple->addItem("PocaMomZ",m_mcIndex,m_PocaMomZ); + + sc=m_tuple->addItem("PocaErrCov",m_mcIndex,m_PocaErrCov,6); + + sc=m_tuple->addItem("ErrorcovMatrix",m_mcIndex,m_ErrorcovMatrix,15); + sc=m_tuple->addItem("D0",m_mcIndex,m_D0); + sc=m_tuple->addItem("phi",m_mcIndex,m_phi); + sc=m_tuple->addItem("omega",m_mcIndex,m_omega); + sc=m_tuple->addItem("Z0",m_mcIndex,m_Z0); + sc=m_tuple->addItem("tanLambda",m_mcIndex,m_tanLambda); + + sc=m_tuple->addItem("ErrorcovMatrix_Origin",m_mcIndex,m_ErrorcovMatrix_Origin,15); + sc=m_tuple->addItem("D0_Origin",m_mcIndex,m_D0_Origin); + sc=m_tuple->addItem("phi_Origin",m_mcIndex,m_phi_Origin); + sc=m_tuple->addItem("omega_Origin",m_mcIndex,m_omega_Origin); + sc=m_tuple->addItem("Z0_Origin",m_mcIndex,m_Z0_Origin); + sc=m_tuple->addItem("tanLambda_Origin",m_mcIndex,m_tanLambda_Origin); + + sc=m_tuple->addItem("mcP_D0",m_mcIndex,mcP_D0); + sc=m_tuple->addItem("mcP_phi",m_mcIndex,mcP_phi); + sc=m_tuple->addItem("mcP_omega",m_mcIndex,mcP_omega); + sc=m_tuple->addItem("mcP_Z0",m_mcIndex,mcP_Z0); + sc=m_tuple->addItem("mcP_tanLambda",m_mcIndex,mcP_tanLambda); + + sc=m_tuple->addItem("pocaPosKal",5,3,m_pocaPosKal); + sc=m_tuple->addItem("pocaMomKal",5,3,m_pocaMomKal); + sc=m_tuple->addItem("pocaMomKalP",m_mcIndex,m_pocaMomKalP,5); + sc=m_tuple->addItem("pocaMomKalPt",m_mcIndex,m_pocaMomKalPt,5); + sc=m_tuple->addItem("chargeKal",m_mcIndex,m_chargeKal,5); + sc=m_tuple->addItem("nDofKal",m_mcIndex,m_nDofKal,5); + sc=m_tuple->addItem("chi2Kal",m_mcIndex,m_chi2Kal,5); + sc=m_tuple->addItem("isFitted",m_mcIndex,m_isFitted,5); + sc=m_tuple->addItem("isFitConverged",m_mcIndex,m_isFitConverged,5); + sc=m_tuple->addItem("isFitConvergedFully",m_mcIndex, + m_isFitConvergedFully,5); + sc=m_tuple->addItem("fittedState",m_mcIndex,m_fittedState,5); + sc=m_tuple->addItem("nHitFailedKal",m_mcIndex,m_nHitFailedKal,5); + sc=m_tuple->addItem("nHitFitted",m_mcIndex,m_nHitFitted,5); + sc=m_tuple->addItem("nDCDigi",m_nDCDigi,0,50000); + sc=m_tuple->addItem("nTruthDCDigi",m_nTruthDCDigi,0,50000); + sc=m_tuple->addItem("nHitKalInput",m_nHitKalInput,0,300000); + //10 is greater than # of tracking detectors + sc=m_tuple->addItem("hitDetID",10,m_nHitDetType); + sc=m_tuple->addItem("nHitWithFitInfo",m_mcIndex,m_nHitWithFitInfo,5); + sc=m_tuple->addItem("nSimDCHit",m_nSimDCHit,0,10000000); + sc=m_tuple->addItem("mdcHitDriftT",m_nSimDCHit,m_mdcHitDriftT); + sc=m_tuple->addItem("mdcHitDriftDl",m_nSimDCHit,m_mdcHitDriftDl); + sc=m_tuple->addItem("mdcHitDriftDr",m_nSimDCHit,m_mdcHitDriftDr); + sc=m_tuple->addItem("mdcHitLr",m_nSimDCHit,m_mdcHitLr); + sc=m_tuple->addItem("mdcHitLayer",m_nSimDCHit,m_mdcHitLayer); + sc=m_tuple->addItem("mdcHitWire",m_nSimDCHit,m_mdcHitWire); + sc=m_tuple->addItem("mdcHitExpDoca",m_nSimDCHit,m_mdcHitExpDoca); + sc=m_tuple->addItem("mdcHitExpMcDoca",m_nSimDCHit,m_mdcHitExpMcDoca); + sc=m_tuple->addItem("mdcHitErr",m_nSimDCHit,m_mdcHitErr); + sc=m_tuple->addItem("exeTime",m_exeTime); + sc=m_tuple->addItem("mdcHitMcTkId",m_nSimDCHit,m_mdcHitMcTkId); + sc=m_tuple->addItem("mdcHitMcLr",m_nSimDCHit,m_mdcHitMcLr); + sc=m_tuple->addItem("mdcHitMcDrift",m_nSimDCHit,m_mdcHitMcDrift); + sc=m_tuple->addItem("mdcHitMcX",m_nSimDCHit,m_mdcHitMcX); + sc=m_tuple->addItem("mdcHitMcY",m_nSimDCHit,m_mdcHitMcY); + sc=m_tuple->addItem("mdcHitMcZ",m_nSimDCHit,m_mdcHitMcZ); + sc=m_tuple->addItem("mcPocaX",m_nSimDCHit,m_mdcHitExpMcPocaX); + sc=m_tuple->addItem("mcPocaY",m_nSimDCHit,m_mdcHitExpMcPocaY); + sc=m_tuple->addItem("mcPocaZ",m_nSimDCHit,m_mdcHitExpMcPocaZ); + sc=m_tuple->addItem("mcPocaWireY",m_nSimDCHit,m_mdcHitExpMcPocaWireY); + sc=m_tuple->addItem("mcPocaWireZ",m_nSimDCHit,m_mdcHitExpMcPocaWireZ); + + sc=m_tuple->addItem("dcDigiChamber",m_nDCDigi,m_dcDigiChamber); + sc=m_tuple->addItem("dcDigiLayer",m_nDCDigi,m_dcDigiLayer); + sc=m_tuple->addItem("dcDigiCell",m_nDCDigi,m_dcDigiCell); + sc=m_tuple->addItem("dcDigiTime",m_nDCDigi,m_dcDigiTime); + sc=m_tuple->addItem("dcDigiDrift",m_nDCDigi,m_dcDigiDrift); + sc=m_tuple->addItem("dcDigiDocaMC",m_nDCDigi,m_dcDigiDocaMC); + sc=m_tuple->addItem("dcDigiPocaOnWireMCX",m_nDCDigi,m_dcDigiPocaOnWireMCX); + sc=m_tuple->addItem("dcDigiPocaOnWireMCY",m_nDCDigi,m_dcDigiPocaOnWireMCY); + sc=m_tuple->addItem("dcDigiWireStartX",m_nDCDigi,m_dcDigiWireStartX); + sc=m_tuple->addItem("dcDigiWireStartY",m_nDCDigi,m_dcDigiWireStartY); + sc=m_tuple->addItem("dcDigiWireStartZ",m_nDCDigi,m_dcDigiWireStartZ); + sc=m_tuple->addItem("dcDigiWireEndX",m_nDCDigi,m_dcDigiWireEndX); + sc=m_tuple->addItem("dcDigiWireEndY",m_nDCDigi,m_dcDigiWireEndY); + sc=m_tuple->addItem("dcDigiWireEndZ",m_nDCDigi,m_dcDigiWireEndZ); + sc=m_tuple->addItem("dcDigiMcMomX",m_nDCDigi,m_dcDigiMcMomX); + sc=m_tuple->addItem("dcDigiMcMomY",m_nDCDigi,m_dcDigiMcMomY); + sc=m_tuple->addItem("dcDigiMcMomZ",m_nDCDigi,m_dcDigiMcMomZ); + sc=m_tuple->addItem("dcDigiMcPosX",m_nDCDigi,m_dcDigiMcPosX); + sc=m_tuple->addItem("dcDigiMcPosY",m_nDCDigi,m_dcDigiMcPosY); + sc=m_tuple->addItem("dcDigiMcPosZ",m_nDCDigi,m_dcDigiMcPosZ); + sc=m_tuple->addItem("firstMomMc",m_firstMomMc); + + sc=m_tuple->addItem("dcDigiDocaExt",m_nDCDigi,m_dcDigiDocaExt); + sc=m_tuple->addItem("dcDigiPocaExtX",m_nDCDigi,m_dcDigiPocaExtX); + sc=m_tuple->addItem("dcDigiPocaExtY",m_nDCDigi,m_dcDigiPocaExtY); + sc=m_tuple->addItem("dcDigiPocaExtZ",m_nDCDigi,m_dcDigiPocaExtZ); + + sc=m_tuple->addItem("nTrackerHitDC",m_nTrackerHitDC,0,1000); + sc=m_tuple->addItem("trackLength",m_nTrackerHitDC,m_trackLength); + sc=m_tuple->addItem("hitMomEdep",m_nTrackerHitDC,m_hitMomEdep); + sc=m_tuple->addItem("truthMomedep",m_nTrackerHitDC,m_truthMomEdep); + sc=m_tuple->addItem("driftDis",m_nTrackerHitDC,m_driftDis); + sc=m_tuple->addItem("FittedDoca",m_nTrackerHitDC,m_FittedDoca); + sc=m_tuple->addItem("truthDoca",m_nTrackerHitDC,m_truthDoca); + sc=m_tuple->addItem("Res",m_nTrackerHitDC,m_Res); + sc=m_tuple->addItem("truthRes",m_nTrackerHitDC,m_truthRes); + sc=m_tuple->addItem("nTrackerHitSDT",m_nTrackerHitSDT); + sc=m_tuple->addItem("nGenFitTrackerHit",m_nGenFitTrackerHit); + debug()<< "Book tuple RecGenfitAlgSDT/recGenfitAlgSDT" << endmsg; + }else{ + warning()<<"Tuple RecGenfitAlgSDT/recGenfitAlgSDT not booked"< start; + if(m_tuple) start=std::chrono::high_resolution_clock::now(); + + ///retrieve silicon Track and TrackHits + const edm4hep::TrackCollection* sdtTrackCol=nullptr; + if(m_SDTTrackCol.exist())sdtTrackCol=m_SDTTrackCol.get(); + if(nullptr==sdtTrackCol || sdtTrackCol->size()<=0) { + debug()<<"TrackCollection not found or sdtTrackCol size=0"<setDebug(m_debug); + + if(!genfitTrack->createGenfitTrackFromEDM4HepTrack(pidType, + sdtTrack, eventStartTime,m_isUseCovTrack)){ + debug()<<"createGenfitTrackFromEDM4HepTrack from SDT track failed!"<addSpacePointsSi(sdtTrack, + m_sigmaHitU,m_sigmaHitV); + }else if(1==m_measurementTypeSi.value()){ + nHitAdded+=genfitTrack->addSiliconMeasurements(sdtTrack, + m_sigmaHitU,m_sigmaHitV); + } + + //add DC hits + if(0==m_measurementTypeDC.value()){ + nHitAdded+=genfitTrack->addSpacePointsDC(sdtTrack, + assoDCHitsCol,m_sigmaHitU,m_sigmaHitV); + }else if(1==m_measurementTypeDC.value()){ + if(m_selectDCHit){ + std::vector selectedHits; + selectHits(sdtTrack,selectedHits); + nHitAdded+=genfitTrack->addWireMeasurementsFromList(selectedHits, + m_sigmaHitU[0],assoDCHitsCol,m_sortMethod,m_truthAmbig, + m_skipCorner,m_skipNear);//mm + std::vector tmp; + selectedHits.swap(tmp); + }else{ + if(m_useNoiseDCHit){ + nHitAdded+=genfitTrack->addWireMeasurementsOnTrack(sdtTrack, + m_sigmaHitU[0],r_NoiseAssociationCol.get(),m_sortMethod,m_truthAmbig, + m_skipCorner,m_skipNear);//mm + } else { + nHitAdded+=genfitTrack->addWireMeasurementsOnTrack(sdtTrack, + m_sigmaHitU[0],r_SmearAssociationCol.get(),m_sortMethod,m_truthAmbig, + m_skipCorner,m_skipNear);//mm + } + } + } + + + // skip events w.o hits + if(0==nHitAdded){ + debug()<printSeed(); + + ///----------------------------------- + ///call genfit fitting procedure + ///----------------------------------- + m_genfitFitter->setDebug(m_debug); + m_genfitFitter->setDebugGenfit(m_debugGenfit); + m_genfitFitter->processTrack(genfitTrack,m_resortHits.value()); + + ///----------------------------------- + ///Store track + ///----------------------------------- + auto dcRecParticle=sdtRecParticleCol->create(); + auto dcRecTrack=sdtRecTrackCol->create(); + + TVector3 pocaToOrigin_pos,pocaToOrigin_mom; + TMatrixDSym pocaToOrigin_cov; + edm4hep::TrackState pocaToOrigin_trackState; + if(!genfitTrack->storeTrack(dcRecParticle,dcRecTrack, + pocaToOrigin_pos,pocaToOrigin_mom,pocaToOrigin_cov, + pidType,m_ndfCut,m_chi2Cut,nFittedDC,nFittedSDT, + ngenfitHit,trackL,hitMom,truthMomEdep,assoDCHitsCol, + driftDis,FittedDoca,truthDoca,Res,truthRes)){ + debug()<<"Fitting failed!"<addEvent(genfitTrack->getTrack()); + m_genfitDisplay->open(); + + using namespace std::chrono_literals; + std::this_thread::sleep_for(1000000000ms); + system("pause"); + }else{ + delete genfitTrack; + } + }//end loop over particle type + ++iSdtTrack; + }//end loop over a track + m_nRecTrack++; + + if(m_tuple) { + m_nTrackerHitDC = nFittedDC; + m_nTrackerHitSDT = nFittedSDT; + m_nGenFitTrackerHit = ngenfitHit; + for(int i=0;i elapsed = finish - start; + debug() << "Elapsed time: " << elapsed.count() << " s"<write(); + + return StatusCode::SUCCESS; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * + +StatusCode RecGenfitAlgSDT::finalize() +{ + MsgStream log(msgSvc(), name()); + info()<< " RecGenfitAlgSDT in finalize()" << endmsg; + + m_genfitFitter->writeHist(); + delete m_genfitFitter; + info()<<"RecGenfitAlgSDT nRecTrack="<0){ + std::cout<<"RecGenfitAlgSDT Success rate = "<getFitStatus(); + int charge= fitState->getCharge(); + + if(m_firstTuple){ + m_nHitKalInput=genfitTrack->getNumPoints(); + debug()<<"m_nHitKalInput "<getNumPointsDet(detIDs[i]); + if(m_debug){ + debug()<<" "<getNumPointsWithFittedInfo(); + m_chi2Kal[iStrack][pidType]=fitState->getChi2(); + m_nDofKal[iStrack][pidType]=fitState->getNdf(); + m_isFitted[iStrack][pidType]=(int)fitState->isFitted(); + m_isFitConverged[iStrack][pidType]=(int) fitState->isFitConverged(); + m_isFitConvergedFully[iStrack][pidType]=(int) fitState->isFitConvergedFully(); + + ///get fitted state of track + TMatrixDSym fittedCov; + TLorentzVector fittedPos; + TVector3 fittedMom; + int fittedState=genfitTrack->getFittedState(fittedPos,fittedMom,fittedCov); + const TLorentzVector seedPos=genfitTrack->getSeedStatePos(); + const TVector3 seedMom=genfitTrack->getSeedStateMom(); + m_fittedState[iStrack][pidType]=fittedState; + HelixClass helix;//mm and GeV + HelixClass helix_origin;//mm and GeV + + + double pos[3]={(fittedPos.X()/dd4hep::mm),(fittedPos.Y()/dd4hep::mm), + (fittedPos.Z()/dd4hep::mm)}; + double mom[3]={(fittedMom.X()),(fittedMom.Y()),(fittedMom.Z())}; + + m_posx[iStrack] = fittedPos.X(); + m_posy[iStrack] = fittedPos.Y(); + m_posz[iStrack] = fittedPos.Z(); + + m_momx[iStrack] = fittedMom.X(); + m_momy[iStrack] = fittedMom.Y(); + m_momz[iStrack] = fittedMom.Z(); + + m_PosMcX[iStrack] = seedPos.X(); + m_PosMcY[iStrack] = seedPos.Y(); + m_PosMcZ[iStrack] = seedPos.Z(); + + m_MomMcX[iStrack] = seedMom.X(); + m_MomMcY[iStrack] = seedMom.Y(); + m_MomMcZ[iStrack] = seedMom.Z(); + + m_PocaPosX[iStrack] = pocaToOrigin_pos.X()*dd4hep::mm; + m_PocaPosY[iStrack] = pocaToOrigin_pos.Y()*dd4hep::mm; + m_PocaPosZ[iStrack] = pocaToOrigin_pos.Z()*dd4hep::mm; + + m_PocaMomX[iStrack] = pocaToOrigin_mom.X(); + m_PocaMomY[iStrack] = pocaToOrigin_mom.Y(); + m_PocaMomZ[iStrack] = pocaToOrigin_mom.Z(); + + for(int i=0;i<6;i++) + { + m_ErrorcovMatrix6[iStrack][i] = fittedCov(i,i); + m_PocaErrCov[iStrack][i] = pocaToOrigin_cov(i,i); + } + + double pocaToOrigin_Pos[3] = {pocaToOrigin_pos.X(),pocaToOrigin_pos.Y(),pocaToOrigin_pos.Z()}; + double pocaToOrigin_Mom[3] = {pocaToOrigin_mom.X(),pocaToOrigin_mom.Y(),pocaToOrigin_mom.Z()}; + TLorentzVector pocaToOrigin_POS; + pocaToOrigin_POS.SetXYZT(pocaToOrigin_pos.X()*dd4hep::mm,pocaToOrigin_pos.Y()*dd4hep::mm, + pocaToOrigin_pos.Z()*dd4hep::mm,999); + helix.Initialize_VP(pos,mom,charge,m_genfitField->getBz(fittedPos.Vect())/GenfitUnit::tesla); + helix_origin.Initialize_VP(pocaToOrigin_Pos,pocaToOrigin_Mom,charge,m_genfitField->getBz(pocaToOrigin_POS.Vect())/GenfitUnit::tesla); + m_pocaMomKalP[iStrack][pidType]=fittedMom.Mag(); + + TMatrixDSym covMatrix_6=pocaToOrigin_cov; + for(int i=0;i<5;i++){ + covMatrix_6[0][i]=pocaToOrigin_cov[0][i]/dd4hep::mm;//d0 column + covMatrix_6[1][i]=pocaToOrigin_cov[1][i]/dd4hep::mm;//omega column + covMatrix_6[2][i]=pocaToOrigin_cov[2][i]/dd4hep::mm;//z0 column + covMatrix_6[i][0]=pocaToOrigin_cov[i][0]/dd4hep::mm;//d0 row + covMatrix_6[i][1]=pocaToOrigin_cov[i][1]/dd4hep::mm;//omega row + covMatrix_6[i][2]=pocaToOrigin_cov[i][2]/dd4hep::mm;//z0 row + } + edm4hep::TrackState trackState_Origin; + CEPC::getTrackStateFromPosMom(trackState_Origin,m_genfitField->getBz(pocaToOrigin_POS.Vect())/GenfitUnit::tesla,pocaToOrigin_pos, + pocaToOrigin_mom,charge,covMatrix_6); + std::array errorCov_Origin; + errorCov_Origin = trackState_Origin.covMatrix; + for(int j=0; j<15; j++) { + m_ErrorcovMatrix_Origin[iStrack][j] = errorCov_Origin[j]; + } + m_D0_Origin[iStrack] = helix_origin.getD0(); + m_phi_Origin[iStrack] = helix_origin.getPhi0(); + m_omega_Origin[iStrack] = helix_origin.getOmega(); + m_Z0_Origin[iStrack] = helix_origin.getZ0(); + m_tanLambda_Origin[iStrack] = helix_origin.getTanLambda(); + + m_evt=m_eventNo; + /// Get fit status + if((0!=fittedState)||(!m_isFitted[pidType])||(m_nDofKal[iStrack][pidType]>m_ndfCut)){ + debug()<<"evt "<size(); + for(auto sdtTrack: *sdtTrackCol){ + m_nSdtTrackHit = sdtTrack.trackerHits_size(); + if(sdtTrack.trackerHits_size()<1e-9) continue; + for(int ihit=0;ihit0) break;//TODO debug for some track only + edm4hep::TrackState trackStat=sdtTrack.getTrackStates(0);//FIXME? + HelixClass helixClass; + helixClass.Initialize_Canonical(trackStat.phi,trackStat.D0, + trackStat.Z0,trackStat.omega,trackStat.tanLambda, + m_genfitField->getBz({0.,0.,0.})/GenfitUnit::tesla); + + TLorentzVector posInit(helixClass.getReferencePoint()[0], + helixClass.getReferencePoint()[1], + helixClass.getReferencePoint()[2],eventStartTime); + m_seedPos[iSdtTrack][0]=posInit.X(); + m_seedPos[iSdtTrack][1]=posInit.Y(); + m_seedPos[iSdtTrack][2]=posInit.Z(); + TVector3 momInit(helixClass.getMomentum()[0], + helixClass.getMomentum()[1],helixClass.getMomentum()[2]); + m_seedMomP[iSdtTrack]=momInit.Mag(); + m_seedMomPt[iSdtTrack]=momInit.Perp(); + m_seedMom[iSdtTrack][0]=momInit.X(); + m_seedMom[iSdtTrack][1]=momInit.Y(); + m_seedMom[iSdtTrack][2]=momInit.Z(); + TVector3 pos,mom; + TMatrixDSym cov(6); + double charge; + CEPC::getPosMomFromTrackState(trackStat, + m_genfitField->getBz({0.,0.,0.})/GenfitUnit::tesla,pos,mom,charge,cov); + for(int i =0;i<6;i++) + { + m_McErrCov[iSdtTrack][i] = cov(i,i); + } + m_seedMomQ[iSdtTrack]=charge; + iSdtTrack++; + } + + const edm4hep::MCParticleCollection* mcParticleCol = nullptr; + const edm4hep::SimTrackerHitCollection* simHitCol=nullptr; + + m_pidIndex=5; + + if((m_mcParticleCol.get())!=nullptr){ + mcParticleCol=m_mcParticleCol.get(); + int iMcParticle=0; + HelixClass helix_mcP; + for(auto mcParticle : *mcParticleCol){ + edm4hep::Vector3f mcPocaMom = mcParticle.getMomentum();//GeV + edm4hep::Vector3d mcPocaPos = mcParticle.getVertex(); + + double mcPos[3]={(mcPocaPos.x),(mcPocaPos.y),(mcPocaPos.z)}; + double mcMom[3]={(mcPocaMom.x),(mcPocaMom.y),(mcPocaMom.z)}; + for(int i=0;i<3;i++){debug()<<"mcPos "<getBz(mcPos)/GenfitUnit::tesla); + + mcP_D0[iMcParticle] = helix_mcP.getD0(); + mcP_phi[iMcParticle] = helix_mcP.getPhi0(); + mcP_omega[iMcParticle] = helix_mcP.getOmega(); + mcP_Z0[iMcParticle] = helix_mcP.getZ0(); + mcP_tanLambda[iMcParticle] = helix_mcP.getTanLambda(); + + debug()<< "debugEvent Bz " << m_genfitField->getBz(mcPos)/GenfitUnit::tesla + << "Tesla mc d0= " << mcP_D0 + << " phi0= " << mcP_phi + << " omega= " << mcP_omega + << " Z0= " << mcP_Z0 + << " tanLambda= " << mcP_tanLambda << endmsg; + + float px=mcPocaMom.x; + float py=mcPocaMom.y; + float pz=mcPocaMom.z; + debug()<<"mc pos("<size(); + + m_nSdtRecTrack=sdtRecTrackCol->size(); + int isdttrack=0; + for(auto sdtRecTrack: *sdtRecTrackCol){ + for(int iHit=0;iHit errorCov; + errorCov = trackStat.covMatrix; + for(int j=0; j<15; j++) { + m_ErrorcovMatrix[isdttrack][j] = errorCov[j]; + if(m_debug)debug()<<"errorCov "<size(); } + int iDCDigi=0; + for(auto dcDigi: *dCDigiCol){ + m_dcDigiChamber[iDCDigi]=m_decoder->get(dcDigi.getCellID(),"chamber"); + m_dcDigiLayer[iDCDigi]=m_decoder->get(dcDigi.getCellID(),"layer"); + m_dcDigiCell[iDCDigi]=m_decoder->get(dcDigi.getCellID(),"cellID"); + m_dcDigiTime[iDCDigi]=dcDigi.getTime(); + m_dcDigiDrift[iDCDigi]=dcDigi.getTime()*m_driftVelocity.value()/10000.; //cm + TVector3 endPointStart(0,0,0); + TVector3 endPointEnd(0,0,0); + m_gridDriftChamber->cellposition(dcDigi.getCellID(),endPointStart, + endPointEnd); + m_dcDigiWireStartX[iDCDigi]=endPointStart.X(); + m_dcDigiWireStartY[iDCDigi]=endPointStart.Y(); + m_dcDigiWireStartZ[iDCDigi]=endPointStart.Z(); + m_dcDigiWireEndX[iDCDigi]=endPointEnd.X(); + m_dcDigiWireEndY[iDCDigi]=endPointEnd.Y(); + m_dcDigiWireEndZ[iDCDigi]=endPointEnd.Z(); + + + + //get information from associated simTrackerHit + //edm4hep::SimTrackerHit dcSimTrackerHit; + auto dcSimTrackerHit=CEPC::getAssoClosestSimTrackerHit(m_DCHitAssociationCol.get(),dcDigi,m_gridDriftChamber,0); + //const edm4hep::MCRecoTrackerAssociationCollection* assoHits=m_DCHitAssociationCol.get(); + m_dcDigiMcMomX[iDCDigi]=dcSimTrackerHit.getMomentum().x*GenfitUnit::GeV; + m_dcDigiMcMomY[iDCDigi]=dcSimTrackerHit.getMomentum().y*GenfitUnit::GeV; + m_dcDigiMcMomZ[iDCDigi]=dcSimTrackerHit.getMomentum().z*GenfitUnit::GeV; + m_dcDigiMcPosX[iDCDigi]=dcSimTrackerHit.getPosition().x*GenfitUnit::mm; + m_dcDigiMcPosY[iDCDigi]=dcSimTrackerHit.getPosition().y*GenfitUnit::mm; + m_dcDigiMcPosZ[iDCDigi]=dcSimTrackerHit.getPosition().z*GenfitUnit::mm; + m_dcDigiDocaMC[iDCDigi]=dcDigi.getTime()*m_driftVelocity.value()/10000.;//cm + TVector3 pocaOnWire=m_gridDriftChamber->wirePos_vs_z(dcDigi.getCellID(), + dcSimTrackerHit.getPosition().z*dd4hep::mm); + m_dcDigiPocaOnWireMCX[iDCDigi]=pocaOnWire.X(); + m_dcDigiPocaOnWireMCY[iDCDigi]=pocaOnWire.Y(); + double firstMom=sqrt(dcSimTrackerHit.getMomentum().x* + dcSimTrackerHit.getMomentum().x+dcSimTrackerHit.getMomentum().y + *dcSimTrackerHit.getMomentum().y+dcSimTrackerHit.getMomentum().z + *dcSimTrackerHit.getMomentum().z); + if(0==m_decoder->get(dcDigi.getCellID(),"layer")){ + m_firstMomMc=firstMom; + if(m_debug.value()>0){ + std::cout<<" firstMomMc "<get(dcDigi.getCellID(),"layer")<<"," + <get(dcDigi.getCellID(),"cellID")<<") truth mom GeV" + <& dcDigiSelected) +{ + //for single track only, FIXME + double eventStartTime=0; + unsigned int pidType=1;//mu + GenfitTrack* genfitTrack=new GenfitTrack(m_genfitField, + m_gridDriftChamber,m_geomSvc); + genfitTrack->setDebug(m_debug); + const edm4hep::MCParticleCollection* mcParticleCol=nullptr; + mcParticleCol=m_mcParticleCol.get();//FIXME get error when call exist() + if(genfitTrack->createGenfitTrackFromMCParticle( + pidType,*(mcParticleCol->begin()), + eventStartTime)){ + const edm4hep::TrackerHitCollection* dCDigiCol=nullptr; + dCDigiCol=m_DCDigiCol.get(); + int iDCDigi=0; + for(auto dcDigi:*dCDigiCol){ + TVector3 poca,pocaDir,pocaOnWire; + + double docaExt=1e9; + bool stopAtBoundary=false; + bool calcJacobianNoise=true; + edm4hep::MCParticle mcParticle=*(mcParticleCol->begin());//FIXME single track only + + genfitTrack->extrapolateToHit(poca,pocaDir,pocaOnWire,docaExt, + mcParticle,dcDigi.getCellID(),pidType,stopAtBoundary,calcJacobianNoise); + m_dcDigiDocaExt[iDCDigi]=docaExt; + m_dcDigiPocaExtX[iDCDigi]=poca.X(); + m_dcDigiPocaExtY[iDCDigi]=poca.Y(); + m_dcDigiPocaExtZ[iDCDigi]=poca.Z(); + double docaMC=dcDigi.getTime()*m_driftVelocity.value()/10000.;//cm + auto dcSimTrackerHit=CEPC::getAssoClosestSimTrackerHit( + m_DCHitAssociationCol.get(),dcDigi,m_gridDriftChamber,0); + + debug()<<" CellID "<get(dcDigi.getCellID(),"layer") + <<","<get(dcDigi.getCellID(),"cellID") + <<") by extdoca:"<m_docaCut*GenfitUnit::mm){ + debug()<<"Skip hit delta doca "<0){ + std::cout<<"selectHits "<size() + <<" after "< +#include +#include "AbsKalmanFitter.h" + +class GenfitFitter; +class GenfitField; +class GenfitTrack; +class IGeomSvc; +class time; +namespace genfit{ + class EventDisplay; +} +namespace dd4hep { + class Detector; + //class rec::CellIDPositionConverter; + namespace DDSegmentation{ + class GridDriftChamber; + class BitFieldCoder; + } +} +namespace edm4hep{ + class EventHeaderCollection; + class MCParticleCollection; + class MutableReconstructedParticle; + class SimTrackerHitCollection; + class TrackCollection; + class TrackerHitCollection; + class MCRecoTrackerAssociationCollection; + class ReconstructedParticle; + class ReconstructedParticleCollection; + class MutableReconstructedParticleCollection; + class TrackerHit; + class Track; +} + +///////////////////////////////////////////////////////////////////////////// + +class RecGenfitAlgSDT:public GaudiAlgorithm { + public: + RecGenfitAlgSDT (const std::string& name, ISvcLocator* pSvcLocator); + StatusCode initialize() override; StatusCode execute() override; StatusCode finalize() override; + + private: + GenfitFitter* m_genfitFitter;//The pointer to a GenfitFitter + const GenfitField* m_genfitField;//The pointer to a GenfitField + + void debugTrack(int iStrack,int pidType,const GenfitTrack* genfitTrack, + TVector3 pocaToOrigin_pos,TVector3 pocaToOrigin_mom, + TMatrixDSym pocaToOrigin_cov); + void debugEvent(const edm4hep::TrackCollection* sdtTrackCol, + const edm4hep::TrackCollection* sdtRecTrackCol, + double eventStartTime,int nFittedSDT); + + void selectHits(const edm4hep::Track&, std::vector& dcDigiSelected); + + DataHandle m_headerCol{ + "EventHeaderCol", Gaudi::DataHandle::Reader, this}; + DataHandle m_DCDigiCol{ + "DigiDCHitCollection", Gaudi::DataHandle::Reader, this}; + //Mc truth + DataHandle m_simVXDHitCol{ + "VXDCollection", Gaudi::DataHandle::Reader, this}; + DataHandle m_simSETHitCol{ + "SETCollection", Gaudi::DataHandle::Reader, this}; + DataHandle m_simSITHitCol{ + "SITCollection", Gaudi::DataHandle::Reader, this}; + DataHandle m_simFTDHitCol{ + "FTDCollection", Gaudi::DataHandle::Reader, this}; + DataHandle m_mcParticleCol{ + "MCParticle", Gaudi::DataHandle::Reader, this}; + DataHandle m_simDCHitCol{ + "DriftChamberHitsCollection" , Gaudi::DataHandle::Reader, this}; + DataHandle + m_DCHitAssociationCol{"DCHitAssociationCollection", + Gaudi::DataHandle::Reader, this}; + DataHandle m_dcTrackCol{ + "DCTrackCollection", Gaudi::DataHandle::Reader, this}; + DataHandle + m_VXDHitAssociationCol{"VXDTrackerHitAssociation", + Gaudi::DataHandle::Reader, this}; + DataHandle + m_SITHitAssociationCol{"SITTrackerHitAssociation", + Gaudi::DataHandle::Reader, this}; + DataHandle + m_SETHitAssociationCol{"SETTrackerHitAssociation", + Gaudi::DataHandle::Reader, this}; + DataHandle + m_FTDHitAssociationCol{"FTDTrackerHitAssociation", + Gaudi::DataHandle::Reader, this}; + + //Track from silicon detectors + DataHandle m_SDTTrackCol{"SDTTrackCollection", + Gaudi::DataHandle::Writer, this}; + DataHandle m_SDTRecTrackCol{"SDTRecTrackCollection", + Gaudi::DataHandle::Writer, this}; + DataHandle m_SDTDebugRecTrackCol{"SDTDebugRecTrackCollection", + Gaudi::DataHandle::Writer, this}; + + //Smear + DataHandle r_SmearAssociationCol{ + "SmearDCHitAssociationCollection", Gaudi::DataHandle::Reader, this}; + + + //Noise + DataHandle r_NoiseAssociationCol{ + "NoiseDCHitAssociationCollection", Gaudi::DataHandle::Reader, this}; + + //Track Finding + DataHandle m_DCTrackFindingCol{ + "DCTrackFindingHitCollection",Gaudi::DataHandle::Reader, this}; + + //Output hits and particles + DataHandle m_SDTRecParticleCol{ + "SDTRecParticleCollection", Gaudi::DataHandle::Writer, this}; + + const unsigned int m_nPDG;//5:e,mu,pi,K,proton + int m_eventNo; + SmartIF m_geomSvc; + dd4hep::OverlayedField m_dd4hepField; + dd4hep::Detector* m_dd4hepDetector; + double m_cell_width; + dd4hep::DDSegmentation::GridDriftChamber* m_gridDriftChamber; + dd4hep::DDSegmentation::BitFieldCoder* m_decoder; + Gaudi::Property m_readout_name{this, + "readout", "DriftChamberHitsCollection"}; + Gaudi::Property m_debug{this,"debug",0}; + Gaudi::Property m_debugGenfit{this,"debugGenfit",0}; + Gaudi::Property m_debugPid{this,"debugPid",-99}; + Gaudi::Property m_eventNoSelection{this,"eventNoSelection",1e9}; + Gaudi::Property > m_sigmaHitU{this, + "sigmaHitU",{0.11, // DC z mm + 0.0028,0.006,0.004,0.004,0.004,0.004, //VXD U + 0.0072, //SIT U 4 layers same resolusion + 0.0072, //SET U + 0.003,0.003,0.0072,0.0072,0.0072,0.0072,0.0072}};//FTD V + //mm, 0:DC, 1~7:VXD, 8:SIT, 9:SET, FTD:10~16 + Gaudi::Property > m_sigmaHitV{this, + "sigmaHitV",{0.11, // DC z mm + 0.0028,0.006,0.004,0.004,0.004,0.004, //VXD V + 0.086, //SIT V + 0.086, //SET V + 0.003,0.003,0.0072,0.0072,0.0072,0.0072,0.0072}};//FTD V + Gaudi::Property m_measurementTypeSi{this,"measurementTypeSi",0}; + //-1: not use, 0, space point, 1, pixel or planer measurement + Gaudi::Property m_measurementTypeDC{this,"measurementTypeDC",0}; + //-1: not use, 0, space point, 1, wire measurement + //Gaudi::Property m_smearHit{this,"smearHit",true}; + //Gaudi::Property m_nSigmaHit{this,"nSigmaHit",5}; + Gaudi::Property m_sortMethod{this,"sortMethod",0}; + Gaudi::Property m_truthAmbig{this,"truthAmbig",true}; + Gaudi::Property m_skipCorner{this,"skipCorner",999.}; + Gaudi::Property m_skipNear{this,"skipNear",0.}; + //Gaudi::Property m_initCovResPos{this,"initCovResPos",1}; + //Gaudi::Property m_initCovResMom{this,"initCovResMom",0.1}; + Gaudi::Property m_isUseCovTrack{this,"isUseCovTrack",false}; + //Fitter type default is DAFRef. + //Candidates are DAF,DAFRef,KalmanFitter and KalmanFitterRefTrack. + Gaudi::Property m_fitterType{this,"fitterType","DAF"}; + Gaudi::Property m_correctBremsstrahlung{this, + "correctBremsstrahlung",false}; + Gaudi::Property m_noMaterialEffects{this, + "noMaterialEffects",false}; + Gaudi::Property m_skipWireMaterial{this, + "skipWireMaterial",true}; + Gaudi::Property m_maxIteration{this,"maxIteration",20}; + Gaudi::Property m_resortHits{this,"resortHits",true}; + Gaudi::Property m_bStart{this,"bStart",100}; + Gaudi::Property m_bFinal{this,"bFinal",0.01}; + Gaudi::Property m_DCCornerCuts{this,"dcCornerCuts",-999}; + Gaudi::Property m_ndfCut{this,"ndfCut",1e9}; + Gaudi::Property m_chi2Cut{this,"chi2Cut",1e9}; + //-1,chargedGeantino;0,1,2,3,4:e,mu,pi,K,proton + Gaudi::Property m_useTruthTrack{this,"useTruthTrack",false}; + Gaudi::Property m_genfitHistRootName{this, + "genfitHistRootName",""}; + Gaudi::Property m_showDisplay{this,"showDisplay",false}; + Gaudi::Property m_fitSiliconOnly{this,"fitSiliconOnly",false}; + Gaudi::Property m_isUseFixedSiHitError{this,"isUseFixedSiHitError",false}; + Gaudi::Property > m_hitError{this,"hitError", + {0.007,0.007,0.03}}; + Gaudi::Property m_extMinDistCut{this,"extMinDistCut",1e-4}; + Gaudi::Property m_multipleMeasurementHandling{this, + "multipleMeasurementHandling", + (int) genfit::eMultipleMeasurementHandling::unweightedClosestToPredictionWire}; + Gaudi::Property m_driftVelocity{this,"drift_velocity",40};//um/ns + Gaudi::Property m_selectDCHit{this,"selectDCHit",false}; + Gaudi::Property m_useNoiseDCHit{this,"useNoiseDCHit",false}; + Gaudi::Property m_docaCut{this,"docaCut",3.3};//mm + int m_fitSuccess[5]; + int m_nRecTrack; + bool m_firstTuple; + //bool m_useRecLRAmbig; + + genfit::EventDisplay* m_genfitDisplay; + + /// tuples + NTuple::Tuple* m_tuple; + NTuple::Item m_run; + NTuple::Item m_evt; + NTuple::Item m_tkId; + NTuple::Item m_mcIndex;//number of navigated mcParicle + NTuple::Matrix m_truthPocaMc;//2 dim matched particle and 3 pos. + NTuple::Array m_seedMomP;//for some track + NTuple::Array m_seedMomPt; + NTuple::Array m_seedMomQ; + NTuple::Matrix m_seedMom; + NTuple::Matrix m_seedPos; + NTuple::Matrix m_pocaPosMc;//2 dim matched particle and 3 pos. + NTuple::Matrix m_pocaMomMc;//2 dim matched particle and 3 mom. + NTuple::Array m_pocaMomMcP;//2 dim matched particle and p + NTuple::Array m_pocaMomMcPt;//2 dim matched particle and pt + NTuple::Matrix m_pocaPosMdc;//pos 0:x,1:y,2:z + NTuple::Matrix m_pocaMomMdc;//mom. 0:px,1:py,2:pz + NTuple::Item m_pidIndex; + NTuple::Matrix m_firstPosKal;//5 hyposis and pos. at first + NTuple::Array m_firstMomKalP;//5 hyposis and mom. at first + NTuple::Array m_firstMomKalPt;//5 hyposis and mom. at first + + NTuple::Matrix m_ErrorcovMatrix6; + NTuple::Array m_posx; + NTuple::Array m_posy; + NTuple::Array m_posz; + + NTuple::Array m_momx; + NTuple::Array m_momy; + NTuple::Array m_momz; + + NTuple::Array m_PosMcX; + NTuple::Array m_PosMcY; + NTuple::Array m_PosMcZ; + + NTuple::Array m_MomMcX; + NTuple::Array m_MomMcY; + NTuple::Array m_MomMcZ; + + + NTuple::Array m_PocaPosX; + NTuple::Array m_PocaPosY; + NTuple::Array m_PocaPosZ; + + NTuple::Array m_PocaMomX; + NTuple::Array m_PocaMomY; + NTuple::Array m_PocaMomZ; + + NTuple::Matrix m_McErrCov; + NTuple::Matrix m_PocaErrCov; + + NTuple::Matrix m_ErrorcovMatrix; + NTuple::Array m_D0; + NTuple::Array m_phi; + NTuple::Array m_omega; + NTuple::Array m_Z0; + NTuple::Array m_tanLambda; + + NTuple::Matrix m_ErrorcovMatrix_Origin; + NTuple::Array m_D0_Origin; + NTuple::Array m_phi_Origin; + NTuple::Array m_omega_Origin; + NTuple::Array m_Z0_Origin; + NTuple::Array m_tanLambda_Origin; + + NTuple::Array mcP_D0; + NTuple::Array mcP_phi; + NTuple::Array mcP_omega; + NTuple::Array mcP_Z0; + NTuple::Array mcP_tanLambda; + + NTuple::Matrix m_pocaPosKal;//5 hyposis and 3 mom. + NTuple::Matrix m_pocaMomKal;//5 hyposis and 3 mom. + NTuple::Matrix m_pocaMomKalP;//5 hyposis and p + NTuple::Matrix m_pocaMomKalPt;//5 hyposis and pt + NTuple::Matrix m_chargeKal; + NTuple::Matrix m_chi2Kal; + NTuple::Matrix m_nDofKal; + NTuple::Matrix m_isFitConverged; + NTuple::Matrix m_isFitConvergedFully; + NTuple::Matrix m_isFitted; + NTuple::Matrix m_fittedState; + NTuple::Item m_nDCDigi; + NTuple::Item m_nTruthDCDigi; + NTuple::Item m_nHitMc; + NTuple::Item m_nSdtTrack; + NTuple::Item m_nSdtTrackHit; + + NTuple::Item m_nSdtRecTrack; + + NTuple::Item m_nSimDCHit; + NTuple::Matrix m_nHitWithFitInfo; + NTuple::Item m_nHitKalInput; + NTuple::Array m_nHitDetType; + NTuple::Array m_mdcHitDriftT; + NTuple::Array m_mdcHitDriftDl; + NTuple::Array m_mdcHitDriftDr; + NTuple::Array m_mdcHitLr; + NTuple::Array m_mdcHitLayer; + NTuple::Array m_mdcHitWire; + NTuple::Array m_mdcHitExpDoca; + NTuple::Array m_mdcHitExpMcDoca; + NTuple::Array m_mdcHitErr; + NTuple::Matrix m_nHitFailedKal; + NTuple::Matrix m_nHitFitted; + NTuple::Item m_exeTime; + //truth + NTuple::Array m_mdcHitMcLr; + NTuple::Array m_mdcHitMcTkId; + NTuple::Array m_mdcHitMcDrift; + NTuple::Array m_mdcHitMcX; + NTuple::Array m_mdcHitMcY; + NTuple::Array m_mdcHitMcZ; + NTuple::Array m_mdcHitExpMcPocaX; + NTuple::Array m_mdcHitExpMcPocaY; + NTuple::Array m_mdcHitExpMcPocaZ; + NTuple::Array m_mdcHitExpMcPocaWireX; + NTuple::Array m_mdcHitExpMcPocaWireY; + NTuple::Array m_mdcHitExpMcPocaWireZ; + + NTuple::Array m_dcDigiChamber; + NTuple::Array m_dcDigiLayer; + NTuple::Array m_dcDigiCell; + NTuple::Array m_dcDigiTime; + NTuple::Array m_dcDigiDrift; + NTuple::Array m_dcDigiDocaMC; + NTuple::Array m_dcDigiPocaOnWireMCX; + NTuple::Array m_dcDigiPocaOnWireMCY; + NTuple::Array m_dcDigiWireStartX; + NTuple::Array m_dcDigiWireStartY; + NTuple::Array m_dcDigiWireStartZ; + NTuple::Array m_dcDigiWireEndX; + NTuple::Array m_dcDigiWireEndY; + NTuple::Array m_dcDigiWireEndZ; + NTuple::Array m_dcDigiMcPosX; + NTuple::Array m_dcDigiMcPosY; + NTuple::Array m_dcDigiMcPosZ; + NTuple::Array m_dcDigiMcMomX; + NTuple::Array m_dcDigiMcMomY; + NTuple::Array m_dcDigiMcMomZ; + NTuple::Array m_dcDigiDocaExt; + NTuple::Array m_dcDigiPocaExtX; + NTuple::Array m_dcDigiPocaExtY; + NTuple::Array m_dcDigiPocaExtZ; + NTuple::Item m_firstMomMc; + + NTuple::Item m_nTrackerHitSDT; + NTuple::Item m_nTrackerHitDC; + NTuple::Item m_nGenFitTrackerHit; + + NTuple::Array m_trackLength; + NTuple::Array m_hitMomEdep; + NTuple::Array m_truthMomEdep; + NTuple::Array m_FittedDoca; + NTuple::Array m_truthDoca; + NTuple::Array m_driftDis; + NTuple::Array m_Res; + NTuple::Array m_truthRes; + +}; +#endif diff --git a/Reconstruction/RecGenfitAlg/src/WireMeasurementDC.cpp b/Reconstruction/RecGenfitAlg/src/WireMeasurementDC.cpp new file mode 100644 index 000000000..0070d7ff0 --- /dev/null +++ b/Reconstruction/RecGenfitAlg/src/WireMeasurementDC.cpp @@ -0,0 +1,194 @@ +/* Copyright 2008-2010, Technische Universitaet Muenchen, + 2014, Ludwig-Maximilians-Universität München + Authors: Tobias Schlüter + + This file is part of GENFIT. + + GENFIT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published + by the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + GENFIT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with GENFIT. If not, see . +*/ +/** @addtogroup genfit + * @{ + */ + +#include "WireMeasurementDC.h" +#include "edm4hep/SimTrackerHit.h" +#include "edm4hep/TrackerHit.h" + +#include +#include + +#include +#include +#include + +#include + + + + +WireMeasurementDC::WireMeasurementDC() + : AbsMeasurement(1), maxDistance_(2), leftRight_(0) +{ + memset(wireEndPoint1_, 0, 3*sizeof(double)); + memset(wireEndPoint2_, 0, 3*sizeof(double)); +} + +WireMeasurementDC::WireMeasurementDC(double driftDistance, double driftDistanceError, const TVector3& endPoint1, const TVector3& endPoint2, int detId, int hitId, TrackPoint* trackPoint) + : AbsMeasurement(1), maxDistance_(2), leftRight_(0) +{ + TVectorD coords(1); + coords[0] = driftDistance; + this->setRawHitCoords(coords); + + TMatrixDSym cov(1); + cov(0,0) = driftDistanceError*driftDistanceError; + this->setRawHitCov(cov); + + this->setWireEndPoints(endPoint1, endPoint2); + + this->setDetId(detId); + this->setHitId(hitId); + this->setTrackPoint(trackPoint); +} + +SharedPlanePtr WireMeasurementDC::constructPlane(const StateOnPlane& state) const { + + // copy state. Neglect covariance. + StateOnPlane st(state); + + TVector3 wire1(wireEndPoint1_); + TVector3 wire2(wireEndPoint2_); + + // unit vector along the wire (V) + TVector3 wireDirection = wire2 - wire1; + wireDirection.SetMag(1.); + + // point of closest approach + const AbsTrackRep* rep = state.getRep(); + rep->extrapolateToLine(st, wire1, wireDirection); + const TVector3& poca = rep->getPos(st); + TVector3 dirInPoca = rep->getMom(st); + dirInPoca.SetMag(1.); + const TVector3& pocaOnWire = wire1 + wireDirection.Dot(poca - wire1)*wireDirection; + + // check if direction is parallel to wire + if (fabs(wireDirection.Angle(dirInPoca)) < 0.01){ + Exception exc("WireMeasurementDC::detPlane(): Cannot construct detector plane, direction is parallel to wire", __LINE__,__FILE__); + throw exc; + } + + // construct orthogonal vector + TVector3 U = wireDirection.Cross(dirInPoca); + // U.SetMag(1.); automatically assured + + return SharedPlanePtr(new DetPlane(pocaOnWire, U, wireDirection)); +} + + +std::vector WireMeasurementDC::constructMeasurementsOnPlane(const StateOnPlane& state) const +{ + double mR = getRawHitCoords()(0); + double mL = -mR; + + MeasurementOnPlane* mopL = new MeasurementOnPlane(TVectorD(1, &mL), + getRawHitCov(), + state.getPlane(), state.getRep(), constructHMatrix(state.getRep())); + + MeasurementOnPlane* mopR = new MeasurementOnPlane(TVectorD(1, &mR), + getRawHitCov(), + state.getPlane(), state.getRep(), constructHMatrix(state.getRep())); + + // set left/right weights + if (leftRight_ < 0) { + mopL->setWeight(1); + mopR->setWeight(0); + } + else if (leftRight_ > 0) { + mopL->setWeight(0); + mopR->setWeight(1); + } + else { + double val = 0.5 * pow(std::max(0., 1 - mR/maxDistance_), 2.); + mopL->setWeight(val); + mopR->setWeight(val); + } + + std::vector retVal; + retVal.push_back(mopL); + retVal.push_back(mopR); + return retVal; +} + +const AbsHMatrix* WireMeasurementDC::constructHMatrix(const AbsTrackRep* rep) const { + if (dynamic_cast(rep) == nullptr) { + Exception exc("WireMeasurementDC default implementation can only handle state vectors of type RKTrackRep!", __LINE__,__FILE__); + throw exc; + } + + return new HMatrixU(); +} + +void WireMeasurementDC::setWireEndPoints(const TVector3& endPoint1, const TVector3& endPoint2) +{ + wireEndPoint1_[0] = endPoint1.X(); + wireEndPoint1_[1] = endPoint1.Y(); + wireEndPoint1_[2] = endPoint1.Z(); + + wireEndPoint2_[0] = endPoint2.X(); + wireEndPoint2_[1] = endPoint2.Y(); + wireEndPoint2_[2] = endPoint2.Z(); +} + +void WireMeasurementDC::setLeftRightResolution(int lr){ + if (lr==0) leftRight_ = 0; + else if (lr<0) leftRight_ = -1; + else leftRight_ = 1; +} + +void WireMeasurementDC::print(){ + std::cout<<"("<getTime() + <<"ns time "<getTime() + <<"ns dd "<getMaxDistance()), + leftRight_(genfitHit->getLeftRightAmbig()), + trackerHit_(genfitHit->getTrackerHit()), + simTrackerHit_(genfitHit->getSimTrackerHit()), + layer_(genfitHit->getLayer()), + cell_(genfitHit->getCell()) +{ + /// New a WireMeasurement + try{ + new (this)WireMeasurementDC( + genfitHit->getDriftDistance(), + genfitHit->getDriftDistanceErr(), + genfitHit->getEnd0(),genfitHit->getEnd1(), + genfitHit->getCellID(),iHit,nullptr); + }catch(genfit::Exception& e){ + std::cout<<"New WireMeasurementDC exception"<. +*/ +/** @addtogroup genfit + * @{ + */ + +#ifndef WireMeasurementDC_h +#define WireMeasurementDC_h + +#include "AbsMeasurement.h" +#include "TrackPoint.h" +#include "SharedPlanePtr.h" +#include "MeasurementOnPlane.h" +#include "StateOnPlane.h" +#include "AbsHMatrix.h" +#include "GenfitHit.h" + +namespace edm4hep{ + class TrackerHit; + class SimTrackerHit; +} + + + +/** @brief Class for measurements in wire detectors (Straw tubes and drift chambers) + * which do not measure the coordinate along the wire. + * + * @author Tobias Schlüter + * + * This is similar to WireMeasurement, but since WireMeasurement + * stores a 7x7 covariance matrix for what is a one-dimensional + * measurement, this class is preferable. Protected inheritance of + * rawHitCoords_ and rawHitCov_ makes it impossible to rewrite + * WireMeasurement, as subclasses will access these members. + * + * This hit class is not valid for arbitrary choices of plane + * orientation: to use it you MUST choose a plane described by u + * and v axes with v coincident with the wire (and u orthogonal + * to it, obviously). + * The hit will be described by 7 coordinates: + * w_x1, w_y1, w_z1, w_x2, w_y2, w_z2, rdrift + * where w_ji (with j = x, y, z and i = 1, 2) are the wire + * extremities coordinates; rdrift = distance from the wire (u + * coordinate in the plane) + * + */ +using namespace genfit; +class WireMeasurementDC : public genfit::AbsMeasurement{ + + public: + WireMeasurementDC(); + WireMeasurementDC(double driftDistance, double driftDistanceError, const TVector3& endPoint1, const TVector3& endPoint2, int detId, int hitId, TrackPoint* trackPoint); + WireMeasurementDC(const GenfitHit* genfitHit,int iHit); + + virtual ~WireMeasurementDC() {;} + + virtual WireMeasurementDC* clone() const override {return new WireMeasurementDC(*this);} + + virtual SharedPlanePtr constructPlane(const StateOnPlane& state) const override; + + /** Hits with a small drift distance get a higher weight, whereas hits with + * big drift distances become weighted down. + * When these initial weights are used by the DAF, the smoothed track will be closer to the real + * trajectory than if both sides are weighted with 0.5 regardless of the drift distance. + * This helps a lot when resolving l/r ambiguities with the DAF. + * The idea is that for the first iteration of the DAF, the wire positions are taken. + * For small drift radii, the wire position does not bend the fit away from the + * trajectory, whereas the wire position for hits with large drift radii is further away + * from the trajectory and will therefore bias the fit if not weighted down. + */ + virtual std::vector constructMeasurementsOnPlane(const StateOnPlane& state) const override; + + virtual const AbsHMatrix* constructHMatrix(const AbsTrackRep*) const override; + + /** Reset the wire end points. + */ + void setWireEndPoints(const TVector3& endPoint1, const TVector3& endPoint2); + + /** Set maximum drift distance. This is used to calculate the start weights of the two + * measurementsOnPlane. + */ + void setMaxDistance(double d){maxDistance_ = d;} + /** + * select how to resolve the left/right ambiguity: + * -1: negative (left) side on vector (wire direction) x (track direction) + * 0: mirrors enter with same weight, DAF will decide. + * 1: positive (right) side on vector (wire direction) x (track direction) + * where the wire direction is pointing from endPoint1 to endPoint2 + */ + void setLeftRightResolution(int lr); + + virtual bool isLeftRigthMeasurement() const {return true;} + double getMaxDistance(){return maxDistance_;} + int getLeftRightResolution() const override {return leftRight_;} + void setTrackerHit(const edm4hep::TrackerHit& trackerHit,int l,int c){ + trackerHit_=&trackerHit; + layer_=l; + cell_=c; + } + void setSimTrackerHit(const edm4hep::SimTrackerHit& simTrackerHit){ + simTrackerHit_=&simTrackerHit; + } + void print(); + + const edm4hep::TrackerHit* getTrackerHit(){return trackerHit_;} + + protected: + + double wireEndPoint1_[3]; //! Wire end point 1 (X, Y, Z) + double wireEndPoint2_[3]; //! Wire end point 2 (X, Y, Z) + double maxDistance_; + double leftRight_; + + int layer_; + int cell_; + const edm4hep::SimTrackerHit* simTrackerHit_; + const edm4hep::TrackerHit* trackerHit_; + + public: + + //ClassDefOverride(WireMeasurementDC, 1) + +}; + +/** @} */ + +#endif // WireMeasurementDC_h diff --git a/Reconstruction/Tracking/include/Tracking/TrackingHelper.h b/Reconstruction/Tracking/include/Tracking/TrackingHelper.h index 4724fb02c..6235b3fb2 100644 --- a/Reconstruction/Tracking/include/Tracking/TrackingHelper.h +++ b/Reconstruction/Tracking/include/Tracking/TrackingHelper.h @@ -62,5 +62,4 @@ inline int getLayer(const edm4hep::TrackerHit hit) { return layer; } - #endif diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp old mode 100755 new mode 100644 diff --git a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp index 6c2876b8b..611b7789d 100644 --- a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp +++ b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp @@ -7,6 +7,7 @@ #include "DD4hep/IDDescriptor.h" #include "DD4hep/Plugins.h" #include "DD4hep/DD4hepUnits.h" +#include "UTIL/ILDConf.h" //external #include "CLHEP/Random/RandGauss.h" @@ -25,37 +26,71 @@ TruthTrackerAlg::TruthTrackerAlg(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc),m_dd4hep(nullptr),m_gridDriftChamber(nullptr), m_decoder(nullptr) { + declareProperty("NoiseSimHitsCollection", w_SimNoiseHCol, + "Handle of the output simulation DC Noise hits collection"); + declareProperty("NoiseDCHitsCollection", w_NoiseHitCol, + "Handle of the input DC Noise hits collection"); + declareProperty("NoiseDCHitAssociationCollection", w_NoiseAssociationCol, + "Handle of the simNoiseDCHit and noiseDC Noise hits collection"); declareProperty("MCParticle", m_mcParticleCol, "Handle of the input MCParticle collection"); + declareProperty("DriftChamberHitsCollection", m_DCSimTrackerHitCol, + "Handle of DC SimTrackerHit collection"); declareProperty("DigiDCHitCollection", m_DCDigiCol, "Handle of DC digi(TrackerHit) collection"); declareProperty("DCHitAssociationCollection", m_DCHitAssociationCol, "Handle of association collection"); declareProperty("DCTrackCollection", m_DCTrackCol, "Handle of DC track collection"); - declareProperty("SiSubsetTrackCollection", m_siSubsetTrackCol, + declareProperty("SubsetTracks", m_siSubsetTrackCol, "Handle of silicon subset track collection"); declareProperty("SDTTrackCollection", m_SDTTrackCol, "Handle of SDT track collection"); - declareProperty("DCRecParticleCollection", m_DCRecParticleCol, - "Handle of drift chamber reconstructed particle collection"); - declareProperty("DCRecParticleAssociationCollection", - m_DCRecParticleAssociationCol, - "Handle of drift chamber reconstructed particle collection"); + declareProperty("VXDTrackerHits", m_VXDTrackerHits, + "Handle of input VXD tracker hit collection"); + declareProperty("SITTrackerHits", m_SITTrackerHits, + "Handle of input SIT tracker hit collection"); + declareProperty("SETTrackerHits", m_SETTrackerHits, + "Handle of input SET tracker hit collection"); + declareProperty("FTDTrackerHits", m_FTDTrackerHits, + "Handle of input FTD tracker hit collection"); + declareProperty("SITSpacePoints", m_SITSpacePointCol, + "Handle of input SIT hit collection"); +// declareProperty("SETSpacePoints", m_SETSpacePointCol, +// "Handle of input SET hit collection"); + declareProperty("FTDSpacePoints", m_FTDSpacePointCol, + "Handle of input FTD hit collection"); + declareProperty("VXDCollection", m_VXDCollection, + "Handle of input VXD hit collection"); + declareProperty("SITCollection", m_SITCollection, + "Handle of input SIT hit collection"); + declareProperty("SETCollection", m_SETCollection, + "Handle of input SET hit collection"); + declareProperty("FTDCollection", m_FTDCollection, + "Handle of input FTD hit collection"); + declareProperty("TruthTrackerHitCollection", m_truthTrackerHitCol, + "Handle of output truth TrackerHit collection"); + declareProperty("SmearTrackerHitCollection", m_SmeartruthTrackerHitCol, + "Handle of output smear truth TrackerHit collection"); + declareProperty("SmearSimHitsCollection", w_SimSmearHCol, + "Handle of output smear SimTrackerHit collection"); + declareProperty("SmearDCHitAssociationCollection", w_SmearAssociationCol, + "Handle of output smear simulationhit and TrackerHit collection"); } + StatusCode TruthTrackerAlg::initialize() { ///Get geometry m_geomSvc=service("GeomSvc"); if (!m_geomSvc) { - error() << "Failed to get GeomSvc." << endmsg; + error()<<"Failed to get GeomSvc."<lcdd(); if (nullptr==m_dd4hep) { - error() << "Failed to get dd4hep::Detector." << endmsg; + error()<<"Failed to get dd4hep::Detector."< (readout.segmentation().segmentation()); if(nullptr==m_gridDriftChamber){ - error() << "Failed to get the GridDriftChamber" << endmsg; + error()<<"Failed to get the GridDriftChamber"<getDecoder(m_readout_name); if (nullptr==m_decoder) { - error() << "Failed to get the decoder" << endmsg; + error()<<"Failed to get the decoder"<book("TruthTrackerAlg/truthTrackerAlg", + CLID_ColumnWiseTuple,"TruthTrackerAlg"); + if(m_tuple){ + StatusCode sc; + sc=m_tuple->addItem("run",m_run); + sc=m_tuple->addItem("evt",m_evt); + sc=m_tuple->addItem("siMom",3,m_siMom); + sc=m_tuple->addItem("siPos",3,m_siPos); + sc=m_tuple->addItem("mcMom",3,m_mcMom); + + sc=m_tuple->addItem("mcPos",3,m_mcPos); + sc=m_tuple->addItem("nDCTrackHit",m_nDCTrackHit,0,100000); + sc=m_tuple->addItem("DriftDistance",m_nDCTrackHit,m_DriftDistance); + sc=m_tuple->addItem("nSmearDCTrackHit",m_nSmearDCTrackHit,0,100000); + sc=m_tuple->addItem("SmearDriftDistance",m_nSmearDCTrackHit,m_SmearDriftDistance); + sc=m_tuple->addItem("nNoiseDCTrackHit",m_nNoiseDCTrackHit,0,100000); + sc=m_tuple->addItem("NoiseDriftDistance",m_nNoiseDCTrackHit,m_NoiseDriftDistance); + + sc=m_tuple->addItem("nSimTrackerHitVXD",m_nSimTrackerHitVXD); + sc=m_tuple->addItem("nSimTrackerHitSIT",m_nSimTrackerHitSIT); + sc=m_tuple->addItem("nSimTrackerHitSET",m_nSimTrackerHitSET); + sc=m_tuple->addItem("nSimTrackerHitFTD",m_nSimTrackerHitFTD); + sc=m_tuple->addItem("nSimTrackerHitDC",m_nSimTrackerHitDC); + sc=m_tuple->addItem("nTrackerHitVXD",m_nTrackerHitVXD); + sc=m_tuple->addItem("nTrackerHitSIT",m_nTrackerHitSIT); + sc=m_tuple->addItem("nTrackerHitSET",m_nTrackerHitSET); + sc=m_tuple->addItem("nTrackerHitFTD",m_nTrackerHitFTD); + sc=m_tuple->addItem("nTrackerHitDC",m_nTrackerHitDC); + sc=m_tuple->addItem("nSpacePointSIT",m_nSpacePointSIT); + sc=m_tuple->addItem("nSpacePointSET",m_nSpacePointSET); + sc=m_tuple->addItem("nSpacePointFTD",m_nSpacePointFTD); + sc=m_tuple->addItem("nHitOnSiTkXVD",m_nHitOnSiTkVXD); + sc=m_tuple->addItem("nHitOnSiTkSIT",m_nHitOnSiTkSIT); + sc=m_tuple->addItem("nHitOnSiTkSET",m_nHitOnSiTkSET); + sc=m_tuple->addItem("nHitOnSiTkFTD",m_nHitOnSiTkFTD); + sc=m_tuple->addItem("nHitOnSdtTkVXD",m_nHitOnSdtTkVXD); + sc=m_tuple->addItem("nHitOnSdtTkSIT",m_nHitOnSdtTkSIT); + sc=m_tuple->addItem("nHitOnSdtTkSET",m_nHitOnSdtTkSET); + sc=m_tuple->addItem("nHitOnSdtTkFTD",m_nHitOnSdtTkFTD); + sc=m_tuple->addItem("nHitOnSdtTkDC",m_nHitOnSdtTkDC); + sc=m_tuple->addItem("nHitSdt",m_nHitOnSdtTk); + sc=m_tuple->addItem("nNoiseOnSdt",m_nNoiseOnSdtTk); + } + } + } return GaudiAlgorithm::initialize(); } StatusCode TruthTrackerAlg::execute() { + info()<<"In execute()"<size()<size(); + } + digiDCHitsCol=m_DCDigiCol.get(); + m_nDCTrackHit = digiDCHitsCol->size(); + if(nullptr==digiDCHitsCol){ + debug()<<"TrackerHitCollection not found"<size()<size()>m_maxDCDigiCut){ + debug()<<"Track cut by m_maxDCDigiCut "<size()>m_maxDCDigiCut) return StatusCode::SUCCESS; - ///Output Track collection - edm4hep::TrackCollection* dcTrackCol=m_DCTrackCol.createAndPut(); - edm4hep::TrackCollection* sdtTrackCol=m_SDTTrackCol.createAndPut(); - edm4hep::ReconstructedParticleCollection* dcRecParticleCol(nullptr); - if(m_writeRecParticle){ - dcRecParticleCol=m_DCRecParticleCol.createAndPut(); - } + // make noise + edm4hep::SimTrackerHitCollection* SimVec = w_SimNoiseHCol.createAndPut(); + edm4hep::MCRecoTrackerAssociationCollection* AssoVec = + w_NoiseAssociationCol.createAndPut(); + edm4hep::TrackerHitCollection* Vec = w_NoiseHitCol.createAndPut(); + ////TODO //Output MCRecoTrackerAssociationCollection collection //const edm4hep::MCRecoTrackerAssociationCollection* @@ -115,50 +230,184 @@ StatusCode TruthTrackerAlg::execute() //} //mcRecoTrackerAssociationCol=m_mcRecoParticleAssociation.get(); - ///Retrieve silicon Track - const edm4hep::TrackCollection* siTrackCol=nullptr; - if(m_siSubsetTrackCol.exist()) siTrackCol=m_siSubsetTrackCol.get(); - if(nullptr!=siTrackCol) { - ///New SDT track - for(auto siTrack:*siTrackCol){ - auto sdtTrack=sdtTrackCol->create(); - edm4hep::TrackState sdtTrackState; - edm4hep::TrackState siTrackStat=siTrack.getTrackStates(0);//FIXME? - sdtTrackState.location=siTrackStat.location; - sdtTrackState.D0=siTrackStat.D0; - sdtTrackState.phi=siTrackStat.phi; - sdtTrackState.omega=siTrackStat.omega; - sdtTrackState.Z0=siTrackStat.Z0; - sdtTrackState.tanLambda=siTrackStat.tanLambda; - sdtTrackState.referencePoint=siTrackStat.referencePoint; - for(int k=0;k<15;k++){ - sdtTrackState.covMatrix[k]=siTrackStat.covMatrix[k]; - } - sdtTrack.addToTrackStates(sdtTrackState); - sdtTrack.setType(siTrack.getType()); - sdtTrack.setChi2(siTrack.getChi2()); - sdtTrack.setNdf(siTrack.getNdf()); - sdtTrack.setDEdx(siTrack.getDEdx()); - sdtTrack.setDEdxError(siTrack.getDEdxError()); - sdtTrack.setRadiusOfInnermostHit(siTrack.getRadiusOfInnermostHit()); - debug()<<"siTrack trackerHits_size="<create(); + + int nVXDHit=0; + int nSITHit=0; + int nSETHit=0; + int nFTDHit=0; + int nDCHitDCTk=0; + int nDCHitSDTTk=0; + + ///Create track with mcParticle + edm4hep::TrackState trackStateMc; + getTrackStateFromMcParticle(mcParticleCol,trackStateMc); + if(m_useTruthTrack.value()||!m_useSi){sdtTk.addToTrackStates(trackStateMc);} + + if(m_useSi){ + ///Retrieve silicon Track + const edm4hep::TrackCollection* siTrackCol=nullptr; + //if(m_siSubsetTrackCol.exist()){ + siTrackCol=m_siSubsetTrackCol.get(); + if(nullptr==siTrackCol){ + debug()<<"SDTTrackCollection is empty"<size()<size()){ + return StatusCode::SUCCESS; } - //TODO tracks - for(auto digiDC:*digiDCHitsCol){ - //if(Sim->MCParti!=current) continue;//TODO - sdtTrack.addToTrackerHits(digiDC); + } + //} + + for(auto siTk:*siTrackCol){ + if(!m_useTruthTrack.value()){ + debug()<<"siTk: "<create(); + + //Create TrackState + edm4hep::TrackState trackStateFirstDCHit; + float charge=trackStateMc.omega/fabs(trackStateMc.omega); + if(m_useFirstHitForDC&&getTrackStateFirstHit(m_DCSimTrackerHitCol, + charge,trackStateFirstDCHit)){ + dcTrack.addToTrackStates(trackStateFirstDCHit); + dcTrack.addToTrackStates(trackStateMc); + }else{ + dcTrack.addToTrackStates(trackStateMc); + dcTrack.addToTrackStates(trackStateFirstDCHit); + } + //sdtTk.addToTrackStates(trackStateFirstDCHit); + + ///Add other track properties + dcTrack.setNdf(dcTrack.trackerHits_size()-5); + + ///Add DC hits to tracks after track state set + if(m_useIdealHit){ + nDCHitDCTk=addIdealHitsToTk(m_DCDigiCol,truthTrackerHitCol,dcTrack, + "DC digi",nDCHitDCTk); + }else{ + nDCHitDCTk=addHitsToTk(m_DCDigiCol,dcTrack,"DC digi",nDCHitDCTk); + } +// if(m_useSi) nDCHitSDTTk=addHitsToTk(m_DCDigiCol,sdtTk,"DC digi",nDCHitSDTTk); + // if(m_useSi) + // { + if(m_smearHits){ + m_nSmearDCTrackHit = smearDCTkhit(m_DCDigiCol, + m_SmeartruthTrackerHitCol,m_DCSimTrackerHitCol, + w_SimSmearHCol,m_DCHitAssociationCol, + w_SmearAssociationCol,m_resX,m_resY,m_resZ); } + + if(!m_useNoiseHits) + { + //nDCHitSDTTk=addHitsToTk(m_DCDigiCol,sdtTk,"DC digi",nDCHitSDTTk); + nDCHitSDTTk=addHitsToTk(m_SmeartruthTrackerHitCol,sdtTk,"DC digi",nDCHitSDTTk); + } else { + m_nNoiseOnSdtTk = makeNoiseHit(SimVec,Vec,AssoVec, + m_SmeartruthTrackerHitCol.get(),w_SmearAssociationCol.get()); + if(debugNoiseHitsCol(Vec))std::cout << " Make noise hits successfully!" << std::endl; + nDCHitSDTTk=addHitsToTk(w_NoiseHitCol,sdtTk, + "DC digi",nDCHitSDTTk); + } + //} + + //track.setType();//TODO + //track.setChi2(gauss(digiDCHitsCol->size-5(),1));//FIXME + //track.setDEdx();//TODO + + debug()<<"dcTrack nHit "<write(); } - ///Convert MCParticle to Track and ReconstructedParticle + return StatusCode::SUCCESS; +} + +StatusCode TruthTrackerAlg::finalize() +{ + return GaudiAlgorithm::finalize(); +} + +void TruthTrackerAlg::getTrackStateFromMcParticle( + const edm4hep::MCParticleCollection* mcParticleCol, + edm4hep::TrackState& trackState) +{ + ///Convert MCParticle to DC Track and ReconstructedParticle debug()<<"MCParticleCol size="<size()<create(); - edm4hep::TrackState trackState; trackState.D0=helix.getD0(); trackState.phi=helix.getPhi0(); trackState.omega=helix.getOmega(); trackState.Z0=helix.getZ0(); trackState.tanLambda=helix.getTanLambda(); trackState.referencePoint=helix.getReferencePoint(); + decltype(trackState.covMatrix) covMatrix; for(int i=0;isize-5(),1));//FIXME - dcTrack.setNdf(digiDCHitsCol->size()-5); - //dcTrack.setDEdx();//TODO - //set hits - double radiusOfInnermostHit=1e9; - debug()<<"digiDCHitsCol size"<size()<MCParti!=current) continue;//TODO - edm4hep::Vector3d digiPos=digiDC.getPosition(); - double r=sqrt(digiPos.x*digiPos.x+digiPos.y*digiPos.y); - if(rcreate(); - ///new ReconstructedParticle - //dcRecParticle.setType();//TODO - double mass=mcParticle.getMass(); - double p=sqrt(mcParticleMomSmeared.x*mcParticleMomSmeared.x+ - mcParticleMomSmeared.y*mcParticleMomSmeared.y+ - mcParticleMomSmeared.z*mcParticleMomSmeared.z); - dcRecParticle.setEnergy(sqrt(mass*mass+p*p)); - dcRecParticle.setMomentum(mcParticleMomSmeared); - dcRecParticle.setReferencePoint(mcParticleVertexSmeared); - dcRecParticle.setCharge(mcParticle.getCharge()); - dcRecParticle.setMass(mass); - //dcRecParticle.setGoodnessOfPID();//TODO - //std::array,10> covMat=?;//TODO - //dcRecParticle.setCovMatrix(covMat);//TODO - //dcRecParticle.setStartVertex();//TODO - edm4hep::ParticleID particleID(0,mcParticle.getPDG(),0,1);//FIXME - dcRecParticle.setParticleIDUsed(particleID); - dcRecParticle.addToTracks(dcTrack); - }//end of write RecParticle + getCircleFromPosMom(pos,mom,B[2]/dd4hep::tesla,mcParticle.getCharge(),m_helixRadius,m_helixXC,m_helixYC); + + debug()<<"dd4hep::mm "<size()<size()){ + edm4hep::SimTrackerHit firstHit; + for(auto dcSimTrackerHit:*col){ + const edm4hep::Vector3f mom=dcSimTrackerHit.getMomentum(); + if(abs(sqrt(mom[0]*mom[0]+mom[1]*mom[1])) covMatrix; + for(int i=0;i<21;i++){covMatrix[i]=100.;}//FIXME + trackState.covMatrix=covMatrix; + debug()<<"first hit trackState "<& + colHandle, edm4hep::TrackerHitCollection*& truthTrackerHitCol, + edm4hep::MutableTrack& track, const char* msg,int nHitAdded) { - return GaudiAlgorithm::finalize(); + if(nHitAdded>0) return nHitAdded; + int nHit=0; + const edm4hep::TrackerHitCollection* col=colHandle.get(); + debug()<<"add "<size()<<" trackerHit"<cellposition(hit.getCellID(),endPointStart, + endPointEnd);//cm + + //calc. doca of helix to wire + TVector3 wire(endPointStart.X()/dd4hep::mm,endPointStart.Y()/dd4hep::mm,0);//to mm + TVector3 center(m_helixXC,m_helixYC,0);//mm + double docaIdeal=(center-wire).Mag()-m_helixRadius;//mm + TVector3 centerFirst(m_helixXCFirst,m_helixYCFirst,0);//mm + double docaIdealFirst=(centerFirst-wire).Mag()-m_helixRadiusFirst;//mm + + //add modified hit + auto tmpHit = truthTrackerHitCol->create(); + tmpHit=hit.clone(); + tmpHit.setTime(fabs(docaIdeal)*1e3/40.);//40#um/ns, drift time in ns + track.addToTrackerHits(tmpHit); + + long long int detID=hit.getCellID(); + debug()<<" addIdealHitsToTk "<get(detID,"layer")<<"," + <get(detID,"cellID")<<") "<size(); +std::cout << " m_nNoiseDCTrackHit = " << m_nNoiseDCTrackHit << std::endl; + int ihit=0; + for(auto Vechit: *Vec) + { + m_NoiseDriftDistance[ihit] = (Vechit.getTime())*m_driftVelocity*1e-3; + ihit++; + } + + return true; +} + +int TruthTrackerAlg::makeNoiseHit(edm4hep::SimTrackerHitCollection* SimVec, + edm4hep::TrackerHitCollection* Vec, + edm4hep::MCRecoTrackerAssociationCollection* AssoVec, + const edm4hep::TrackerHitCollection* digiDCHitsCol, + const edm4hep::MCRecoTrackerAssociationCollection* assoHits) +{ + int nNoise = 0; + //loop all hits + for(unsigned int i = 0; i<(digiDCHitsCol->size()); i++) + { + edm4hep::TrackerHit mcHit = digiDCHitsCol->at(i); + unsigned long long wcellid = mcHit.getCellID(); + + float pocaTime = mcHit.getTime(); + + // store trackHit + auto trkHit = Vec->create(); + + trkHit.setCellID(wcellid); + trkHit.setTime(pocaTime); + trkHit.setEDep(mcHit.getEDep()); + //trkHit.setEdx(mcHit.getEdx()); + trkHit.setPosition(mcHit.getPosition()); + trkHit.setCovMatrix(mcHit.getCovMatrix()); + for(int iAsso=0;iAsso<(int) assoHits->size();iAsso++) + { + + if(assoHits->at(iAsso).getRec().getCellID()==wcellid ) + { + auto SimtrkHit = SimVec->create(); + + SimtrkHit.setCellID(assoHits->at(iAsso).getSim().getCellID()); + SimtrkHit.setEDep(assoHits->at(iAsso).getSim().getEDep()); + SimtrkHit.setTime(assoHits->at(iAsso).getSim().getTime()); + SimtrkHit.setPathLength(assoHits->at(iAsso).getSim().getPathLength()); + SimtrkHit.setQuality(assoHits->at(iAsso).getSim().getQuality()); + SimtrkHit.setPosition(assoHits->at(iAsso).getSim().getPosition()); + SimtrkHit.setMomentum(assoHits->at(iAsso).getSim().getMomentum()); + SimtrkHit.setMCParticle(assoHits->at(iAsso).getSim().getMCParticle()); + + auto asso = AssoVec->create(); + asso.setRec(trkHit); + asso.setSim(SimtrkHit); + asso.setWeight(SimtrkHit.getEDep()/trkHit.getEDep()); + } + } + + double prob1 = fRandom.Uniform(1.); + if(prob1 > m_fHitPurity) continue; + + float noiseTime = m_pocaTime*fRandom.Uniform(1.); + if(noiseTime>pocaTime) continue; + nNoise++; + trkHit.setTime(noiseTime); + + for(int iAsso=0;iAsso<(int) assoHits->size();iAsso++) + { + + if(assoHits->at(iAsso).getRec().getCellID()==wcellid ) + { + + AssoVec->at(iAsso).setRec(trkHit); + + } + + } + + } + +std::cout << " Truth Noise number = " << nNoise << std::endl; + return nNoise; +} + +int TruthTrackerAlg::smearDCTkhit(DataHandle& + colHandle,DataHandle& smearCol, + DataHandle& SimDCHitCol, + DataHandle& SimSmearDCHitCol, + DataHandle& AssoDCHitCol, + DataHandle& AssoSmearDCHitCol, + double resX, double resY, double resZ) +{ + int nHit=0; + const edm4hep::TrackerHitCollection* col=colHandle.get(); + auto smearTrackerHitCol = smearCol.createAndPut(); + + auto simDCCol = SimDCHitCol.get(); + auto SimVec = SimSmearDCHitCol.createAndPut(); + + auto assoHits = AssoDCHitCol.get(); + auto AssoVec = AssoSmearDCHitCol.createAndPut(); + + //sort,FIXME + for(auto hit:*col){ + + auto smearHit=smearTrackerHitCol->create(); + + float truthD = (hit.getTime())*m_driftVelocity*1e-3; + m_DriftDistance[nHit] = truthD; + + truthD+=gRandom->Gaus(0,fabs(m_resX)); + m_SmearDriftDistance[nHit] = truthD; + float smearT = (truthD*1e3)/m_driftVelocity; + + smearHit.setTime(smearT); + smearHit.setCellID(hit.getCellID()); + smearHit.setType(hit.getCellID()); + smearHit.setQuality(hit.getQuality()); + smearHit.setEDep(hit.getEDep()); + smearHit.setEDepError(hit.getEDepError()); + //smearHit.setEdx(hit.getEdx()); + smearHit.setPosition(hit.getPosition()); + smearHit.setCovMatrix(hit.getCovMatrix()); + smearHit.addToRawHits(hit.getObjectID()); + + ++nHit; + for(int iAsso=0;iAsso<(int) assoHits->size();iAsso++) + { + + if(assoHits->at(iAsso).getRec().getCellID()== smearHit.getCellID()) + { + auto SimtrkHit = SimVec->create(); + + SimtrkHit.setCellID(assoHits->at(iAsso).getSim().getCellID()); + SimtrkHit.setEDep(assoHits->at(iAsso).getSim().getEDep()); + SimtrkHit.setTime(assoHits->at(iAsso).getSim().getTime()); + SimtrkHit.setPathLength(assoHits->at(iAsso).getSim().getPathLength()); + SimtrkHit.setQuality(assoHits->at(iAsso).getSim().getQuality()); + SimtrkHit.setPosition(assoHits->at(iAsso).getSim().getPosition()); + SimtrkHit.setMomentum(assoHits->at(iAsso).getSim().getMomentum()); + SimtrkHit.setMCParticle(assoHits->at(iAsso).getSim().getMCParticle()); + + auto asso = AssoVec->create(); + asso.setRec(smearHit); + asso.setSim(SimtrkHit); + asso.setWeight(SimtrkHit.getEDep()/smearHit.getEDep()); + } + } + } + std::cout << " SimSmear size = " << SimVec->size() << " Asso size = " << AssoVec->size() << std::endl; + std::cout << " Sim size = " << simDCCol->size() << " Asso size = " << assoHits->size() << std::endl; + + return nHit; +} + +int TruthTrackerAlg::addHitsToTk(DataHandle& + colHandle, edm4hep::MutableTrack& track, const char* msg,int nHitAdded) +{ + if(nHitAdded>0) return nHitAdded; + int nHit=0; + const edm4hep::TrackerHitCollection* col=colHandle.get(); + debug()<<"add "<size()<<" trackerHit"<& colHandle, + edm4hep::TrackerHitCollection*& truthTrackerHitCol, + edm4hep::MutableTrack& track, const char* msg,int nHitAdded) +{ + if(nHitAdded>0) return nHitAdded; + int nHit=0; + const edm4hep::SimTrackerHitCollection* col=colHandle.get(); + for(auto simTrackerHit:*col){ + auto trackerHit=truthTrackerHitCol->create(); + if(m_skipSecondaryHit&&simTrackerHit.isProducedBySecondary()) { + debug()<<"skip secondary simTrackerHit "< cov; + cov[0]=resolution[0]*resolution[0]; + cov[1]=0.; + cov[2]=resolution[1]*resolution[1]; + cov[3]=0.; + cov[4]=0.; + cov[5]=resolution[2]*resolution[2]; + trackerHit.setCovMatrix(cov); + debug()<<"add simTrackerHit "<0) return nHitAdded; + int nHit=0; + for(unsigned int iHit=0;iHit& col) +{ + const edm4hep::TrackerHitCollection* c=col.get(); + if(nullptr!=c) return c->size(); + return 0; +} + +int TruthTrackerAlg::simTrackerHitColSize(DataHandle& col) +{ + const edm4hep::SimTrackerHitCollection* c=col.get(); + if(nullptr!=c) return c->size(); + return 0; +} +//unit length is mm +void TruthTrackerAlg::getCircleFromPosMom(double pos[3],double mom[3], + double Bz,double q,double& helixRadius,double& helixXC,double& helixYC) +{ + double FCT = 2.99792458E-4;//mm + double pxy = sqrt(mom[0]*mom[0]+mom[1]*mom[1]); + helixRadius = pxy / (FCT*Bz); + double phiMomRefPoint = atan2(mom[1],mom[0]); + helixXC= pos[0] + helixRadius*cos(phiMomRefPoint-M_PI*0.5*q); + helixYC= pos[1] + helixRadius*sin(phiMomRefPoint-M_PI*0.5*q); } diff --git a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h index 30efd0a8a..563d98543 100644 --- a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h +++ b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h @@ -4,6 +4,9 @@ #include "GaudiAlg/GaudiAlgorithm.h" #include "k4FWCore/DataHandle.h" #include "DD4hep/Fields.h" +#include "GaudiKernel/NTuple.h" + +#include "TRandom3.h" class IGeomSvc; namespace dd4hep { @@ -15,10 +18,14 @@ namespace dd4hep { } namespace edm4hep { class MCParticleCollection; + class SimTrackerHitCollection; class TrackerHitCollection; class TrackCollection; - class MCRecoTrackerAssociationCollection; + class Track; + class MutableTrack; + class TrackState; class ReconstructedParticleCollection; + class MCRecoTrackerAssociationCollection; class MCRecoParticleAssociationCollection; } @@ -32,47 +39,199 @@ class TruthTrackerAlg: public GaudiAlgorithm virtual StatusCode finalize() override; private: + void getTrackStateFromMcParticle(const edm4hep::MCParticleCollection* + mcParticleCol, edm4hep::TrackState& stat); + int addSimHitsToTk(DataHandle& + colHandle, edm4hep::TrackerHitCollection*& truthTrackerHitCol, + edm4hep::MutableTrack& track, const char* msg,int nHitAdded); + int smearDCTkhit(DataHandle& + colHandle,DataHandle& smearCol, + DataHandle& SimDCHitCol, + DataHandle& SimSmearDCHitCol, + DataHandle& AssoDCHitCol, + DataHandle& AssoSmearDCHitCol, + double resX, double resY, double resZ); + int addHitsToTk(DataHandle& + colHandle, edm4hep::MutableTrack& track, const char* msg,int nHitAdded); + int addIdealHitsToTk(DataHandle& + colHandle, edm4hep::TrackerHitCollection*& truthTrackerHitCol, + edm4hep::MutableTrack& track, const char* msg,int nHitAdded); + + int addHotsToTk(edm4hep::Track& sourceTrack,edm4hep::MutableTrack& + targetTrack, int hitType,const char* msg,int nHitAdded); + int nHotsOnTrack(edm4hep::Track& track, int hitType); + int trackerHitColSize(DataHandle& hitCol); + int simTrackerHitColSize(DataHandle& hitCol); + bool getTrackStateFirstHit( + DataHandle& dcSimTrackerHitCol, + float charge,edm4hep::TrackState& trackState); SmartIF m_geomSvc; dd4hep::Detector* m_dd4hep; dd4hep::OverlayedField m_dd4hepField; dd4hep::DDSegmentation::GridDriftChamber* m_gridDriftChamber; dd4hep::DDSegmentation::BitFieldCoder* m_decoder; + void debugEvent(); + //unit length is mm + void getCircleFromPosMom(double pos[3],double mom[3], + double Bz,double q,double& helixRadius,double& helixXC,double& helixYC); + int makeNoiseHit(edm4hep::SimTrackerHitCollection* SimVec, + edm4hep::TrackerHitCollection* Vec, + edm4hep::MCRecoTrackerAssociationCollection* AssoVec, + const edm4hep::TrackerHitCollection* digiDCHitsCol, + const edm4hep::MCRecoTrackerAssociationCollection* assoHits); + bool debugNoiseHitsCol(edm4hep::TrackerHitCollection* Vec); - //reader + //reader + // DataHandle m_NoiseHitCol{ + // "NoiseDCHitsCollection", Gaudi::DataHandle::Reader, this}; DataHandle m_mcParticleCol{ "MCParticle", Gaudi::DataHandle::Reader, this}; + DataHandle m_DCSimTrackerHitCol{ + "DriftChamberHitsCollection", Gaudi::DataHandle::Reader, this}; DataHandle m_DCDigiCol{ "DigiDCHitCollection", Gaudi::DataHandle::Reader, this}; DataHandle m_DCHitAssociationCol{ "DCHitAssociationCollection", Gaudi::DataHandle::Reader, this}; DataHandle - m_siSubsetTrackCol{ "SiSubsetTrackCollection", + m_siSubsetTrackCol{ "SubsetTracks", Gaudi::DataHandle::Reader, this}; + DataHandle m_SITSpacePointCol{ + "SITSpacePoints" , Gaudi::DataHandle::Reader, this}; + // DataHandle m_SETSpacePointCol{ + // "SETSpacePoints" , Gaudi::DataHandle::Reader, this}; + DataHandle m_FTDSpacePointCol{ + "FTDSpacePoints" , Gaudi::DataHandle::Reader, this}; + DataHandle m_VXDTrackerHits{ + "VXDTrackerHits" , Gaudi::DataHandle::Reader, this}; + DataHandle m_SETTrackerHits{ + "SETTrackerHits" , Gaudi::DataHandle::Reader, this}; + DataHandle m_SITTrackerHits{ + "SITTrackerHits" , Gaudi::DataHandle::Reader, this}; + DataHandle m_FTDTrackerHits{ + "FTDTrackerHits" , Gaudi::DataHandle::Reader, this}; + DataHandle m_VXDCollection{ + "VXDCollection" , Gaudi::DataHandle::Reader, this}; + DataHandle m_SETCollection{ + "SETCollection" , Gaudi::DataHandle::Reader, this}; + DataHandle m_SITCollection{ + "SITCollection" , Gaudi::DataHandle::Reader, this}; + DataHandle m_FTDCollection{ + "FTDCollection" , Gaudi::DataHandle::Reader, this}; //writer DataHandle m_DCTrackCol{ "DCTrackCollection", Gaudi::DataHandle::Writer, this}; DataHandle m_SDTTrackCol{ "SDTTrackCollection", Gaudi::DataHandle::Writer, this}; - DataHandle m_DCRecParticleCol{ - "DCRecParticleCollection", Gaudi::DataHandle::Writer, this}; - DataHandle - m_DCRecParticleAssociationCol{ - "DCRecMCRecoParticleAssociationCollection", - Gaudi::DataHandle::Writer, this}; + DataHandle m_truthTrackerHitCol{ + "TruthTrackerHitCollection", Gaudi::DataHandle::Writer, this}; + + // Smear hit + DataHandle w_SimSmearHCol{ + "SmearSimHitsCollection", Gaudi::DataHandle::Writer, this}; + DataHandle m_SmeartruthTrackerHitCol{ + "SmearTrackerHitCollection", Gaudi::DataHandle::Writer, this}; + DataHandle w_SmearAssociationCol{ + "SmearDCHitAssociationCollection", Gaudi::DataHandle::Writer, this}; + // make noise hit + DataHandle w_SimNoiseHCol{ + "NoiseSimHitsCollection", Gaudi::DataHandle::Writer, this}; + DataHandle w_NoiseHitCol{ + "NoiseDCHitsCollection", Gaudi::DataHandle::Writer, this}; + DataHandle w_NoiseAssociationCol{ + "NoiseDCHitAssociationCollection", Gaudi::DataHandle::Writer, this}; + //readout for getting segmentation Gaudi::Property m_readout_name{this, "readout", "DriftChamberHitsCollection"}; - Gaudi::Property m_writeRecParticle{this,"writeRecParticle",false}; + Gaudi::Property m_hist{this,"hist",false}; + Gaudi::Property m_useDC{this,"useDC",true}; + Gaudi::Property m_useSi{this,"useSi",true}; + Gaudi::Property m_useTruthTrack{this,"useTruthTrack",false}; + Gaudi::Property m_useSiTruthHit{this,"useSiTruthHit",false}; + Gaudi::Property m_skipSecondaryHit{this,"skipSecondaryHit",true}; + Gaudi::Property m_useFirstHitForDC{this,"useFirstHitForDC",false}; + Gaudi::Property m_useSiSpacePoint{this,"useSiSpacePoint",false}; + Gaudi::Property m_useIdealHit{this,"useIdealHit",false}; + + Gaudi::Property m_useNoiseHits{this,"useNoiseHits",false}; + Gaudi::Property m_smearHits{this,"smearHits",false}; + + Gaudi::Property m_momentumCut{this,"momentumCut",0.1};//momentum cut for the first hit Gaudi::Property m_resPT{this,"resPT",0};//ratio Gaudi::Property m_resPz{this,"resPz",0};//ratio + Gaudi::Property m_resX{this,"resX",0.11};//mm + Gaudi::Property m_resY{this,"resY",0.11};//mm + Gaudi::Property m_resZ{this,"resZ",0.11};//mm + Gaudi::Property m_driftVelocity{this,"driftVelocity",40};// um/us Gaudi::Property m_resMomPhi{this,"resMomPhi",0};//radian Gaudi::Property m_resMomTheta{this,"resMomTheta",0};//radian Gaudi::Property m_resVertexX{this,"resVertexX",0.003};//3um Gaudi::Property m_resVertexY{this,"resVertexY",0.003};//3um Gaudi::Property m_resVertexZ{this,"resVertexZ",0.003};//3um Gaudi::Property m_maxDCDigiCut{this,"maxDigiCut",1e6}; + Gaudi::Property > m_resVXD{this,"resVXD",{0.003,0.003,0.003}};//mm + Gaudi::Property > m_resSIT{this,"resSIT",{0.003,0.003,0.003}};//mm + Gaudi::Property > m_resSET{this,"resSET",{0.003,0.003,0.003}};//mm + Gaudi::Property > m_resFTDPixel{this,"resFTDPixel",{0.003,0.003,0.003}};//mm + Gaudi::Property > m_resFTDStrip{this,"resFTDStrip",{0.003,0.003,0.003}};//mm + + Gaudi::Property m_fHitPurity{this,"fHitPurity",0.1}; + Gaudi::Property m_pocaTime { this, "pocaTime", 225};// ns + + double m_helixRadius,m_helixXC,m_helixYC; + double m_helixRadiusFirst,m_helixXCFirst,m_helixYCFirst; + + NTuple::Tuple* m_tuple; + NTuple::Item m_run; + NTuple::Item m_evt; + NTuple::Array m_siMom; + NTuple::Array m_siPos; + NTuple::Array m_mcMom; + NTuple::Array m_mcPos; + + NTuple::Item m_nDCTrackHit; + NTuple::Item m_nSmearDCTrackHit; + NTuple::Item m_nNoiseDCTrackHit; + NTuple::Array m_DriftDistance; + NTuple::Array m_SmearDriftDistance; + NTuple::Array m_NoiseDriftDistance; + + NTuple::Item m_nSimTrackerHitVXD; + NTuple::Item m_nSimTrackerHitSIT; + NTuple::Item m_nSimTrackerHitSET; + NTuple::Item m_nSimTrackerHitFTD; + NTuple::Item m_nSimTrackerHitDC; + NTuple::Item m_nTrackerHitVXD; + NTuple::Item m_nTrackerHitSIT; + NTuple::Item m_nTrackerHitSET; + NTuple::Item m_nTrackerHitFTD; + NTuple::Item m_nTrackerHitDC; + NTuple::Item m_nTrackerHitErrVXD; + NTuple::Item m_nTrackerHitErrSIT; + NTuple::Item m_nTrackerHitErrSET; + NTuple::Item m_nTrackerHitErrFTD; + NTuple::Item m_nSpacePointSIT; + NTuple::Item m_nSpacePointSET; + NTuple::Item m_nSpacePointFTD; + NTuple::Item m_nSpacePointErrVXD; + NTuple::Item m_nSpacePointErrSIT; + NTuple::Item m_nSpacePointErrSET; + NTuple::Item m_nSpacePointErrFTD; + NTuple::Item m_nHitOnSiTkVXD; + NTuple::Item m_nHitOnSiTkSIT; + NTuple::Item m_nHitOnSiTkSET; + NTuple::Item m_nHitOnSiTkFTD; + NTuple::Item m_nHitOnSdtTkVXD; + NTuple::Item m_nHitOnSdtTkSIT; + NTuple::Item m_nHitOnSdtTkSET; + NTuple::Item m_nHitOnSdtTkFTD; + NTuple::Item m_nHitOnSdtTkDC; + NTuple::Item m_nHitOnSdtTk; + NTuple::Item m_nNoiseOnSdtTk; + + TRandom3 fRandom; }; #endif diff --git a/Utilities/DataHelper/CMakeLists.txt b/Utilities/DataHelper/CMakeLists.txt index 8398a6557..74ec620fb 100644 --- a/Utilities/DataHelper/CMakeLists.txt +++ b/Utilities/DataHelper/CMakeLists.txt @@ -15,9 +15,11 @@ gaudi_add_library(DataHelperLib src/TrackerHitExtended.cc src/TrackExtended.cc src/TrackHitPair.cc - src/TrackerHitHelper.cpp + src/TrackerHitHelper.cpp + src/TrackHelper.cc LINK EDM4HEP::edm4hep EDM4HEP::edm4hepDict + DetSegmentation GSL::gsl ${CLHEP_LIBRARIES} Identifier diff --git a/Utilities/DataHelper/include/DataHelper/HelixClass.h b/Utilities/DataHelper/include/DataHelper/HelixClass.h index 794f81586..f573ee9d5 100644 --- a/Utilities/DataHelper/include/DataHelper/HelixClass.h +++ b/Utilities/DataHelper/include/DataHelper/HelixClass.h @@ -69,6 +69,7 @@ class HelixClass { * - particle charge : q;
* - magnetic field (in Tesla) : B
*/ + void Initialize_VP(double * pos, double * mom, double q, double B); void Initialize_VP(float * pos, float * mom, float q, float B); /** @@ -266,35 +267,35 @@ class HelixClass { float getCharge(); private: - float _momentum[3]; // momentum @ ref point + float _momentum[3]; // momentum @ ref point float _referencePoint[3]; // coordinates @ ref point - float _phi0; // phi0 in canonical parameterization - float _d0; // d0 in canonical parameterisation - float _z0; // z0 in canonical parameterisation - float _omega; // signed curvuture in canonical parameterisation - float _tanLambda; // TanLambda - float _pxy; // Transverse momentum - float _charge; // Particle Charge - float _bField; // Magnetic field (assumed to point to Z>0) - float _radius; // radius of circle in XY plane - float _xCentre; // X of circle centre - float _yCentre; // Y of circle centre - float _phiRefPoint; // Phi w.r.t. (X0,Y0) of circle @ ref point - float _phiAtPCA; // Phi w.r.t. (X0,Y0) of circle @ PCA - float _xAtPCA; // X @ PCA - float _yAtPCA; // Y @ PCA - float _pxAtPCA; // PX @ PCA - float _pyAtPCA; // PY @ PCA - float _phiMomRefPoint; // Phi of Momentum vector @ ref point - float _const_pi; // PI - float _const_2pi; // 2*PI - float _const_pi2; // PI/2 - float _FCT; // 2.99792458E-4 + double _phi0; // phi0 in canonical parameterization + double _d0; // d0 in canonical parameterisation + double _z0; // z0 in canonical parameterisation + double _omega; // signed curvuture in canonical parameterisation + double _tanLambda; // TanLambda + double _pxy; // Transverse momentum + double _charge; // Particle Charge + double _bField; // Magnetic field (assumed to point to Z>0) + double _radius; // radius of circle in XY plane + double _xCentre; // X of circle centre + double _yCentre; // Y of circle centre + double _phiRefPoint; // Phi w.r.t. (X0,Y0) of circle @ ref point + double _phiAtPCA; // Phi w.r.t. (X0,Y0) of circle @ PCA + double _xAtPCA; // X @ PCA + double _yAtPCA; // Y @ PCA + double _pxAtPCA; // PX @ PCA + double _pyAtPCA; // PY @ PCA + double _phiMomRefPoint; // Phi of Momentum vector @ ref point + double _const_pi; // PI + double _const_2pi; // 2*PI + double _const_pi2; // PI/2 + double _FCT; // 2.99792458E-4 float _xStart[3]; // Starting point of track segment float _xEnd[3]; // Ending point of track segment - float _bZ; - float _phiZ; + double _bZ; + double _phiZ; }; diff --git a/Utilities/DataHelper/include/DataHelper/TrackHelper.h b/Utilities/DataHelper/include/DataHelper/TrackHelper.h new file mode 100644 index 000000000..cb76b6bea --- /dev/null +++ b/Utilities/DataHelper/include/DataHelper/TrackHelper.h @@ -0,0 +1,19 @@ +#ifndef TRACKHELPER_H +#define TRACKHELPER_H +#include "edm4hep/TrackState.h" +#include "TMatrixDSym.h" +#include "TVector3.h" + +namespace CEPC{ + //get track position and momentum from TrackState + void getPosMomFromTrackState(const edm4hep::TrackState& trackState, + double Bz, TVector3& pos,TVector3& mom,double& charge, + TMatrixDSym& covMatrix_6); + + //Set track state from position, momentum and charge + void getTrackStateFromPosMom(edm4hep::TrackState& trackState,double Bz, + TVector3 pos,TVector3 mom,double charge,TMatrixDSym covMatrix_6); + +} + +#endif diff --git a/Utilities/DataHelper/include/DataHelper/TrackerHitHelper.h b/Utilities/DataHelper/include/DataHelper/TrackerHitHelper.h index fcd9dd146..319bbebf6 100644 --- a/Utilities/DataHelper/include/DataHelper/TrackerHitHelper.h +++ b/Utilities/DataHelper/include/DataHelper/TrackerHitHelper.h @@ -2,13 +2,29 @@ #define TrackerHitHelper_H #include "edm4hep/TrackerHit.h" +#include "edm4hep/SimTrackerHit.h" +#include "edm4hep/MCRecoTrackerAssociationCollection.h" +#include "DDSegmentation/Segmentation.h" +#include "DetSegmentation/GridDriftChamber.h" #include +//namespace dd4hep { +// class Detector; +// namespace DDSegmentation{ +// class GridDriftChamber; +// } +//} + namespace CEPC{ std::array GetCovMatrix(edm4hep::TrackerHit& hit, bool useSpacePointerBuilderMethod = false); float GetResolutionRPhi(edm4hep::TrackerHit& hit); float GetResolutionZ(edm4hep::TrackerHit& hit); std::array ConvertToCovXYZ(float dU, float thetaU, float phiU, float dV, float thetaV, float phiV, bool useSpacePointBuilderMethod = false); + const edm4hep::SimTrackerHit getAssoClosestSimTrackerHit( + const edm4hep::MCRecoTrackerAssociationCollection* assoHits, + const edm4hep::TrackerHit trackerHit, + const dd4hep::DDSegmentation::GridDriftChamber* segmentation, + int docaMehtod); } #endif diff --git a/Utilities/DataHelper/src/HelixClass.cc b/Utilities/DataHelper/src/HelixClass.cc index 591acc514..a9ef6ca70 100644 --- a/Utilities/DataHelper/src/HelixClass.cc +++ b/Utilities/DataHelper/src/HelixClass.cc @@ -12,6 +12,75 @@ HelixClass::HelixClass() { HelixClass::~HelixClass() {} +void HelixClass::Initialize_VP(double * pos, double * mom, double q, double B) { + _referencePoint[0] = (float) pos[0]; + _referencePoint[1] = (float) pos[1]; + _referencePoint[2] = (float) pos[2]; + _momentum[0] = (float) mom[0]; + _momentum[1] = (float) mom[1]; + _momentum[2] = (float) mom[2]; + _charge = q; + _bField = B; + _pxy = sqrt(mom[0]*mom[0]+mom[1]*mom[1]); + _radius = _pxy / (_FCT*B); + _omega = q/_radius; + _tanLambda = mom[2]/_pxy; + _phiMomRefPoint = atan2(mom[1],mom[0]); + _xCentre = pos[0] + _radius*cos(_phiMomRefPoint-_const_pi2*q); + _yCentre = pos[1] + _radius*sin(_phiMomRefPoint-_const_pi2*q); + _phiRefPoint = atan2(pos[1]-_yCentre,pos[0]-_xCentre); + _phiAtPCA = atan2(-_yCentre,-_xCentre); + _phi0 = -_const_pi2*q + _phiAtPCA; + while (_phi0<0) _phi0+=_const_2pi; + while (_phi0>=_const_2pi) _phi0-=_const_2pi; + _xAtPCA = _xCentre + _radius*cos(_phiAtPCA); + _yAtPCA = _yCentre + _radius*sin(_phiAtPCA); + // _d0 = -_xAtPCA*sin(_phi0) + _yAtPCA*cos(_phi0); + + if (q>0) { + _d0 = double(q)*_radius - double(sqrt(_xCentre*_xCentre+_yCentre*_yCentre)); + } + else { + _d0 = double(q)*_radius + double(sqrt(_xCentre*_xCentre+_yCentre*_yCentre)); + } + +// if (fabs(_d0)>0.001 ) { +// std::cout << "New helix : " << std::endl; +// std::cout << " Position : " << pos[0] +// << " " << pos[1] +// << " " << pos[2] << std::endl; +// std::cout << " Radius = " << _radius << std::endl; +// std::cout << " RC = " << sqrt(_xCentre*_xCentre+_yCentre*_yCentre) << std::endl; +// std::cout << " D0 = " << _d0 << std::endl; +// } + + _pxAtPCA = _pxy*cos(_phi0); + _pyAtPCA = _pxy*sin(_phi0); + float deltaPhi = _phiRefPoint - _phiAtPCA; + float xCircles = -pos[2]*q/(_radius*_tanLambda) - deltaPhi; + xCircles = xCircles/_const_2pi; + int nCircles; + int n1,n2; + + if (xCircles >= 0.) { + n1 = int(xCircles); + n2 = n1 + 1; + } + else { + n1 = int(xCircles) - 1; + n2 = n1 + 1; + } + + if (fabs(n1-xCircles) < fabs(n2-xCircles)) { + nCircles = n1; + } + else { + nCircles = n2; + } + _z0 = pos[2] + _radius*_tanLambda*q*(deltaPhi + _const_2pi*nCircles); + +} + void HelixClass::Initialize_VP(float * pos, float * mom, float q, float B) { _referencePoint[0] = pos[0]; _referencePoint[1] = pos[1]; diff --git a/Utilities/DataHelper/src/TrackHelper.cc b/Utilities/DataHelper/src/TrackHelper.cc new file mode 100644 index 000000000..9e1f437fb --- /dev/null +++ b/Utilities/DataHelper/src/TrackHelper.cc @@ -0,0 +1,134 @@ +#include "DataHelper/TrackHelper.h" +#include "DD4hep/DD4hepUnits.h" +#include "DataHelper/HelixClass.h" +#include "TVector3.h" +#include "TLorentzVector.h" +#include + +//get track position and momentum from TrackState +void CEPC::getPosMomFromTrackState(const edm4hep::TrackState& trackState, + double Bz, TVector3& pos,TVector3& mom,double& charge, + TMatrixDSym& covMatrix_6){ + double D0=trackState.D0; + double phi=trackState.phi; + double omega=trackState.omega; + double Z0=trackState.Z0; + double tanLambda=trackState.tanLambda; + + ::edm4hep::Vector3f referencePoint=trackState.referencePoint; + charge=omega/fabs(omega); + const double FCT = 2.99792458E-4; + double radius = 1./fabs(omega); + double pxy = FCT*Bz*radius; + for(int i=0;i<3;i++){ + pos[i]=referencePoint[i]; + } + mom[0]=pxy*cos(phi); + mom[1]=pxy*sin(phi); + mom[2]=pxy*tanLambda; + TMatrixDSym covMatrix_5(5); + ///< lower triangular covariance matrix of the track parameters. + /// the order of parameters is d0, phi, omega, z0, tan(lambda). + std::array covMatrix=trackState.covMatrix; + int k=0; + for(int i=0;i<5;i++){ + for(int j=0;j<5;j++){ + if(i>=j) { covMatrix_5(i,j)=covMatrix[k++]; } + } + } + + ///Error propagation + //V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , helix.covariance() = V(X) + TMatrix mS(covMatrix_6.GetNrows(),covMatrix_5.GetNrows()); + mS.Zero(); + mS[0][0]=-sin(phi); + mS[0][1]=-1*D0*cos(phi); + mS[1][0]=cos(phi); + mS[1][1]=-1*D0*sin(phi); + mS[2][3]=1; + mS[3][1]=FCT*Bz*(1/omega)*sin(phi); + mS[3][2]=charge*FCT*Bz*(1/(omega*omega))*cos(phi); + mS[4][1]=-1*FCT*Bz*(1/omega)*cos(phi); + mS[4][2]=charge*FCT*Bz*(1/(omega*omega))*sin(phi); + mS[5][2]=charge*tanLambda*Bz*(1/(omega*omega)); + mS[5][4]=-FCT*Bz/omega*(1+tanLambda*tanLambda); + + covMatrix_6= covMatrix_5.Similarity(mS); + + const char* m_name="TrackHelper"; + int m_debug=0; + if(m_debug>=2){ + //std::cout<=2){ + std::cout<<__FILE__<<" pos mom"< covMatrix; + int k=0; + int k1; + for(int i=0;i<6;i++){ + for(int j=0;j<6;j++){ + if(i>=j) { + k1=k++; + //covMatrix[k++]=covMatrix_5(i,j); + if(5==i){ + covMatrix[k1]=-999; + continue; + } + covMatrix[k1]=covMatrix_5(i,j); + } + } + } + + trackState.covMatrix = covMatrix; +} diff --git a/Utilities/DataHelper/src/TrackerHitHelper.cpp b/Utilities/DataHelper/src/TrackerHitHelper.cpp index e12eb886c..f12b16873 100644 --- a/Utilities/DataHelper/src/TrackerHitHelper.cpp +++ b/Utilities/DataHelper/src/TrackerHitHelper.cpp @@ -7,6 +7,8 @@ #include "CLHEP/Vector/ThreeVector.h" #include "CLHEP/Vector/Rotation.h" #include +#include "TVector3.h" +#include "DD4hep/DD4hepUnits.h" std::array CEPC::GetCovMatrix(edm4hep::TrackerHit& hit, bool useSpacePointBuilderMethod){ if(hit.isAvailable()){ @@ -115,3 +117,44 @@ std::array CEPC::ConvertToCovXYZ(float dU, float thetaU, float phiU, f } return cov; } + +const edm4hep::SimTrackerHit CEPC::getAssoClosestSimTrackerHit( + const edm4hep::MCRecoTrackerAssociationCollection* assoHits, + const edm4hep::TrackerHit trackerHit, + const dd4hep::DDSegmentation::GridDriftChamber* segmentation, + int docaMehtod) +{ + std::vector hits; + for(auto assoHit: *assoHits){ + if(assoHit.getRec()==trackerHit){ + hits.push_back(assoHit.getSim()); + } + } + edm4hep::SimTrackerHit minSimTrackerHit; + double min_distance = 999 ; + double tmp_distance =0.; + for(auto hit:hits){ + unsigned long long wcellid=hit.getCellID(); + TVector3 pos(hit.getPosition()[0]*dd4hep::mm, hit.getPosition()[1]*dd4hep::mm, hit.getPosition()[2]*dd4hep::mm); + + TVector3 sim_mon(hit.getMomentum()[0],hit.getMomentum()[1],hit.getMomentum()[2]); + float Steplength = hit.getPathLength(); + TVector3 pos_start=pos - 0.5 * Steplength * sim_mon.Unit(); + TVector3 pos_end=pos + 0.5 * Steplength * sim_mon.Unit(); + TVector3 hitPosition; + TVector3 PCA; + tmp_distance = segmentation->Distance(wcellid,pos_start,pos_end,hitPosition,PCA); + tmp_distance = tmp_distance/dd4hep::mm; //mm + + if(tmp_distance < min_distance){ + min_distance = tmp_distance; + //pos_x = hitPosition.x(); //pos.x(); + //pos_y = hitPosition.y(); //pos.y(); + //pos_z = pos.z(); + //momx = hit.getMomentum()[0]; + //momy = hit.getMomentum()[1]; + minSimTrackerHit=hit; + } + } + return minSimTrackerHit; +} diff --git a/run.sh b/run.sh index 634faa1d6..461f88885 100755 --- a/run.sh +++ b/run.sh @@ -69,6 +69,7 @@ function run-job() { lcg_platform=x86_64-centos7-gcc11-opt lcg_version=103.0.2 + bldtool=${CEPCSW_BLDTOOL} # make, ninja # set in env var check-cepcsw-envvar || exit -1