C++ TVector3::Phi方法代码示例

int get_circ_wg_efield_vert(TVector3 &pos, TVector3 &efield)
{
  //returns the TE_11 efield inside an infinite circ. waveguide 
  //efield normalized the jackson way with 1/cm units
  //waveguide extends in x-direction
  // wavelength of 27 GHz radiation is 1.1 cm
  double p11 = 1.841;            //1st zero of the derivate of bessel function
  //k11 is angular wavenumber for cutoff frequency for TE11
  double k11 = p11 / tl_data.rO;
  //convert position to cylindrical
  pos.SetZ(0);
  double radius = pos.Mag();
  double phi = pos.Phi();//azimuthal position
  //double phi = acos(pos.Y()/radius);//azimuthal position, see definition of phase
  double J1 = TMath::BesselJ1(k11 * radius);//this term cancels in dot product w/ vel
  double Jp = TMath::BesselJ0(k11 * radius) - TMath::BesselJ1(k11 * radius) / k11 / radius;
  double e_amp = 1.63303/tl_data.rO;//cm
  int status = 0;

  if (radius > tl_data.rO) {
    cout << "Problem!!! Electron hit a wall! At radius " << radius << endl;
    status = 1;
  }
  efield.SetX(e_amp * (J1 / k11 / radius * cos(phi) * sin(phi) - Jp * sin(phi) * cos(phi)));
  efield.SetY(e_amp * (J1 / k11 / radius * sin(phi) * sin(phi) + Jp * cos(phi) * cos(phi)));
  efield.SetZ(0);                //only true for TE mode
  if (radius == 0) {
    Jp = 1.0/2; 
    phi = 0;
    efield.SetX(e_amp * (1. / 2 * cos(phi) * sin(phi) - Jp * sin(phi) * cos(phi)));
    efield.SetY(e_amp * (1. / 2 * sin(phi) * sin(phi) + Jp * cos(phi) * cos(phi)));
  }
  return status;
}

  static Bool_t is_in_phi_fiducial(const TVector3 &v) {
    // converted from emid::is_in_phi_fiducial
    
    const double PHI_CRACK_WIDTH=0.02; 
    
    double phi = v.Phi()>=0 ? v.Phi() : (v.Phi()+TMath::TwoPi()); 
    
    // CC phi cracks
    Bool_t outside_phi_crack=true;
    if (TMath::Abs(v.Z()) < 150.) {
      double phimod=fmod(phi+0.1,TMath::Pi()/16.);
      outside_phi_crack=(phimod<0.1-PHI_CRACK_WIDTH) || 
	(phimod>0.1+PHI_CRACK_WIDTH);
    }
    return outside_phi_crack;
  }

void KVFAZIALNS2016::BuildFAZIA()
{
   //Build geometry of FAZIASYM
   //All telescopes are : Si(300µm)-Si(500µm)-CsI(10cm)
   //No attempt has been made to implement real thicknesses
   //
   Info("BuildFAZIA", "Compact geometry, %f cm from target",
        fFDist);

   TGeoVolume* top = gGeoManager->GetTopVolume();

   Double_t distance_block_cible = fFDist * KVUnits::cm;
   Double_t thick_si1 = 300 * KVUnits::um;
   TGeoTranslation trans;
   trans.SetDz(distance_block_cible + thick_si1 / 2.);

   KVFAZIABlock* block = new KVFAZIABlock;

   TGeoRotation rot1, rot2;
   TGeoHMatrix h;
   TGeoHMatrix* ph = 0;
   Double_t theta = 0;
   Double_t phi = 0;

   Double_t theta_min = fFThetaMin;//smallest lab polar angle in degrees
   Double_t centre_hole = 2.*tan(theta_min * TMath::DegToRad()) * distance_block_cible;
   Double_t dx = (block->GetTotalSideWithBlindage()) / 2.;

   TVector3 centre;
   for (Int_t bb = 0; bb < fNblocks; bb += 1) {

      if (bb == 1)        centre.SetXYZ(-1 * (dx - centre_hole / 2), -dx - centre_hole / 2, distance_block_cible);
      else if (bb == 2)   centre.SetXYZ(-1 * (dx + centre_hole / 2), dx - centre_hole / 2, distance_block_cible);
      else if (bb == 3)   centre.SetXYZ(-1 * (-dx + centre_hole / 2), dx + centre_hole / 2, distance_block_cible);
      else if (bb == 0)   centre.SetXYZ(-1 * (-dx - centre_hole / 2), -dx + centre_hole / 2, distance_block_cible);
      else if (bb == 4)   centre.SetXYZ(-1 * (-dx - centre_hole / 2), -3 * dx + centre_hole / 2, distance_block_cible); //centre.SetXYZ(-1 * (dx - centre_hole / 2), -3 * dx - centre_hole / 2, distance_block_cible);
      else {
         Warning("BuildFAZIA", "Block position definition is done only for %d blocks", fNblocks);
      }
      theta = centre.Theta() * TMath::RadToDeg();
      phi = centre.Phi() * TMath::RadToDeg();
      printf("BLK #%d => theta=%1.2lf - phi=%1.2lf\n", bb, theta, phi);

      rot2.SetAngles(phi + 90., theta, 0.);
      rot1.SetAngles(-1.*phi, 0., 0.);
      h = rot2 * trans * rot1;
      ph = new TGeoHMatrix(h);
      top->AddNode(block, bb, ph);
   }

   // add telescope for elastic scattering monitoring
//   RutherfordTelescope();
   // Change default geometry import angular range for rutherford telescope
   SetGeometryImportParameters(.25, 1., 1.84);
}

//_________________________________________________________________
Double_t KVTenseur3::GetPhiPlan(void)
{
   // Obtention du Phi du plan de reaction (en degrès).
   // Cet angle suit les conventions de KaliVeda i.e. il est entre 0 et 360 degrees
   if (is_diago == 0)
      Diago();

   TVector3 vp = GetVep(3);

   if (vp.Z() < 0)
      vp = -vp;
   Double_t phi = vp.Phi() * TMath::RadToDeg();
   return (phi < 0 ? 360. + phi : phi);
}

//.........这里部分代码省略.........
	  float smearing = overSmear.randOverSmearing(phoSCEta[i],phoR9[i],isInGAP_EB(i),overSmearSyst);
	  phoRegrE[i] *= (1 + smearing); 
	  phoE[i]     *= (1 + smearing); 
	  /// from MassFactorized in gglobe:   energyCorrectedError[ipho] *=(l.pho_isEB[ipho]) ? 1.07 : 1.045 ;
	  float smearFactor = 1;
	  if( _config->setup() == "ReReco2011" ) smearFactor = fabs(phoSCEta[i]) < 1.45 ? 1.07: 1.045;
	  phoRegrErr[i] *= smearFactor;
	}  
      
      
      // 2. reweighting of photon id variables (R9...)
      for (int i=0; i<nPho; ++i) ReweightMC_phoIdVar(i);
      //// ************************************************************* ////
    }
    
    //// ********** Apply regression energy ************* ////
    float phoStdE[500];
    for( int i=0; i<nPho; ++i)
      if( fabs(phoSCEta[i]) <= 2.5 ) {
	if( isData ){
	  float enCorrSkim = 1;//enScaleSkimEOS.energyScale( phoR9[i], phoSCEta[i], run);
	  float phoEnScale = enScale.energyScale( phoR9[i], phoSCEta[i], run)/enCorrSkim;
	  phoRegrE[i]  *= phoEnScale;
	  phoE[i]      *= phoEnScale;
	}
	phoStdE[i] = phoE[i];
	phoE[i]   = phoRegrE[i];
	
	/// transform calo position abd etaVtx, phiVtx with SC position
	for( int x = 0 ; x < 3; x++ )  phoCaloPos[i][x] = phoSCPos[i][x];
	for( int ivtx = 0 ; ivtx < nVtxBS; ivtx++ ) {
	  TVector3 xxi = getCorPhotonTVector3(i,ivtx);
	  phoEtaVtx[i][ivtx] = xxi.Eta();
	  phoPhiVtx[i][ivtx] = xxi.Phi();
	}
	
	/// additionnal smearing to go to data energy resolution
	phoRegrSmear[i] = phoE[i]*overSmearHCP.meanOverSmearing(phoSCEta[i],phoR9[i],isInGAP_EB(i),0);
	phoEt[i] = phoE[i] / cosh(phoEta[i]);  
      }    
    //// ************************************************* ////
     

    /// lepton selection
    int iElecVtx(-1), iMuonVtx(-1);
    vector<int> elecIndex = selectElectronsHCP2012( wei, iElecVtx );
    vector<int> muonIndex = selectMuonsHCP2012(     wei, iMuonVtx );

    vector<int>             event_vtx; 
    vector<int>             event_ilead ;
    vector<int>             event_itrail;
    vector<TLorentzVector>  event_plead ;
    vector<TLorentzVector>  event_ptrail;
    TLorentzVector leptag;
    
    int lepCat = -1;
    bool exitLoop = false;
    for( int ii = 0     ; ii < nPho ; ++ii ) {
      for( int jj = (ii+1); jj < nPho ; ++jj ) {
	// Preselection 2nd leg
	if (DoPreselection && !preselectPhoton(ii,phoEt[ii])) continue;
	if (DoPreselection && !preselectPhoton(jj,phoEt[jj])) continue;
	
	/// define i, j locally, so when they are inverted this does not mess up the loop
	int i = ii; 
	int j = jj;

StatusCode CorrectECalBarrelSliWinCluster::execute() {
  // Get the input collection with clusters
  const fcc::CaloClusterCollection* inClusters = m_inClusters.get();
  fcc::CaloClusterCollection* correctedClusters = m_correctedClusters.createAndPut();

  // for single particle events compare with truth particles
  TVector3 momentum;
  double phiVertex = 0;
  double etaVertex = 0;
  double thetaVertex = 0;
  double zVertex = 0;
  const auto particle = m_particle.get();
  const auto vertex = m_vertex.get();
  if (particle->size() == 1 && vertex->size() == 1) {
    for(const auto& part : *particle) {
      momentum = TVector3(part.core().p4.px, part.core().p4.py, part.core().p4.pz);
      etaVertex = momentum.Eta();
      phiVertex = momentum.Phi();
      zVertex =  vertex->begin()->position().z;
      thetaVertex = 2 * atan( exp( - etaVertex ) );
    verbose() << " vertex eta " << etaVertex << "   phi = " << phiVertex << " theta = " << thetaVertex << " z = " << zVertex << endmsg;
   }
  }

  // TODO change that so all systems can be used
  uint systemId = m_systemId[0];
  const dd4hep::DDSegmentation::FCCSWGridPhiEta* segmentation = nullptr;
  if (m_segmentationPhiEta[systemId] != nullptr) {
    segmentation = m_segmentationPhiEta[systemId];
  }

  std::vector<TLorentzVector> clustersMassInv;
  std::vector<TLorentzVector> clustersMassInvScaled;
  for (const auto& cluster : *inClusters) {
    double oldEnergy = 0;
    TVector3 pos(cluster.core().position.x, cluster.core().position.y, cluster.core().position.z);
    double oldEta = pos.Eta();
    double oldPhi = pos.Phi();
    for (auto cell = cluster.hits_begin(); cell != cluster.hits_end(); cell++) {
      oldEnergy += cell->core().energy;
    }
    verbose() << " OLD ENERGY = " << oldEnergy << " from " << cluster.hits_size() << " cells" << endmsg;
    verbose() << " OLD CLUSTER ENERGY = " << cluster.core().energy << endmsg;

    // Do everything only using the first defined calorimeter (default: Ecal barrel)
    double oldEtaId = -1;
    double oldPhiId = -1;
    if (m_segmentationPhiEta[systemId] != nullptr) {
      oldEtaId = int(floor((oldEta + 0.5 * segmentation->gridSizeEta() - segmentation->offsetEta()) / segmentation->gridSizeEta()));
      oldPhiId = int(floor((oldPhi + 0.5 * segmentation->gridSizePhi() - segmentation->offsetPhi()) / segmentation->gridSizePhi()));
    }

    // 0. Create new cluster, copy information from input
    fcc::CaloCluster newCluster = correctedClusters->create();
    double energy = 0;
    newCluster.core().position.x = cluster.core().position.x;
    newCluster.core().position.y = cluster.core().position.y;
    newCluster.core().position.z = cluster.core().position.z;
    for (auto cell = cluster.hits_begin(); cell != cluster.hits_end(); cell++) {
      if (m_segmentationMulti[systemId] != nullptr) {
        segmentation = dynamic_cast<const dd4hep::DDSegmentation::FCCSWGridPhiEta*>(&m_segmentationMulti[systemId]->subsegmentation(cell->core().cellId));
        oldEtaId = int(floor((oldEta + 0.5 * segmentation->gridSizeEta() - segmentation->offsetEta()) / segmentation->gridSizeEta()));
        oldPhiId = int(floor((oldPhi + 0.5 * segmentation->gridSizePhi() - segmentation->offsetPhi()) / segmentation->gridSizePhi()));
      }
      if (m_decoder[systemId]->get(cell->core().cellId, "system") == systemId) {
        uint layerId = m_decoder[systemId]->get(cell->core().cellId, "layer");
        if(m_nPhiFinal[layerId] > 0 && m_nEtaFinal[layerId] > 0) {
          uint etaId = m_decoder[systemId]->get(cell->core().cellId, "eta");
          uint phiId = m_decoder[systemId]->get(cell->core().cellId, "phi");
          if ( etaId >= (oldEtaId - m_halfEtaFin[layerId]) &&  etaId <= (oldEtaId + m_halfEtaFin[layerId]) &&
               phiId >= phiNeighbour((oldPhiId - m_halfPhiFin[layerId]), segmentation->phiBins()) &&  phiId <= phiNeighbour((oldPhiId + m_halfPhiFin[layerId]), segmentation->phiBins()) ) {
            if (m_ellipseFinalCluster) {
              if ( pow( (etaId - oldEtaId) / (m_nEtaFinal[layerId] / 2.), 2) + pow( (phiId - oldPhiId) / (m_nPhiFinal[layerId] / 2.), 2) < 1) {
                newCluster.addhits(*cell);
                energy += cell->core().energy;
              }
            } else {
              newCluster.addhits(*cell);
              energy += cell->core().energy;
            }
          }
        }
      }
    }
    newCluster.core().energy = energy;

    // 1. Correct eta position with log-weighting
    double sumEnFirstLayer = 0;
    // get current pseudorapidity
    std::vector<double> sumEnLayer;
    std::vector<double> sumEnLayerSorted;
    std::vector<double> sumEtaLayer;
    std::vector<double> sumWeightLayer;
    sumEnLayer.assign(m_numLayers, 0);
    sumEnLayerSorted.assign(m_numLayers, 0);
    sumEtaLayer.assign(m_numLayers, 0);
    sumWeightLayer.assign(m_numLayers, 0);
    // first check the energy deposited in each layer
    for (auto cell = newCluster.hits_begin(); cell != newCluster.hits_end(); cell++) {
      dd4hep::DDSegmentation::CellID cID = cell->core().cellId;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      tr = new TGeoTranslation(coefx[nt - 1]*shift, coefy[nt - 1]*shift, thick_si2 / 2. + distance_si2_si1);
      quartet->AddNode(si, ndet++, tr);
      ((TGeoNodeMatrix*)quartet->GetNodes()->Last())->SetName(Form("DET_SI2-T%d", nt));

      shift = side_si / 2 + inter_si / 2 + adjust_csi;

      csi = gGeoManager->MakeTrd2(Form("DET_CSI-T%d", nt), CesiumIodide, side_csi_front / 2, side_csi_back / 2, side_csi_front / 2, side_csi_back / 2, thick_csi / 2.);
      tr = new TGeoTranslation(coefx[nt - 1]*shift, coefy[nt - 1]*shift, thick_csi / 2. + distance_csi_si2);
      quartet->AddNode(csi, ndet++, tr);
      ((TGeoNodeMatrix*)quartet->GetNodes()->Last())->SetName(Form("DET_CSI-T%d", nt));

   }

   Int_t nbl = 1;
   TGeoVolume* blindage = 0;
//Double_t thick_bld = thick_si1+distance_si2_si1+thick_si2;
   /* l'epaisseur du si1 est compris dans la distance_si2_si1 */
   Double_t thick_bld = distance_si2_si1 + thick_si2;
   Double_t shift_bld = (side_si + inter_si) / 2.;
///Croix inter quartet
//
// Separation des 4 télescopes
//
//
   blindage = gGeoManager->MakeBox("DEADZONE_BLINDAGE_1", Plomb, inter_si / 2, (side_si + inter_si / 2), thick_bld / 2.);
   //printf("%s\n", blindage->GetMaterial()->GetTitle());
   tr = new TGeoTranslation(0, 0, thick_bld / 2.);
   quartet->AddNode(blindage, nbl++, tr);

   blindage = gGeoManager->MakeBox("DEADZONE_BLINDAGE_2", Plomb, (side_si / 2), inter_si / 2, thick_bld / 2.);
   tr = new TGeoTranslation(-1 * shift_bld, 0, thick_bld / 2.);
   quartet->AddNode(blindage, nbl++, tr);
   tr = new TGeoTranslation(+1 * shift_bld, 0, thick_bld / 2.);
   quartet->AddNode(blindage, nbl++, tr);

///Contour de l ensemble du quartet
//
//Délimiation des bords exterieurs
//
//
   shift_bld = (side_si + inter_si);

   blindage = gGeoManager->MakeBox("DEADZONE_BLINDAGE_3", Plomb, (side_si + inter_si / 2), inter_si / 2, thick_bld / 2.);
   tr = new TGeoTranslation(0, shift_bld, thick_bld / 2.);
   quartet->AddNode(blindage, nbl++, tr);
   tr = new TGeoTranslation(0, -1 * shift_bld, thick_bld / 2.);
   quartet->AddNode(blindage, nbl++, tr);
///
   blindage = gGeoManager->MakeBox("DEADZONE_BLINDAGE_4", Plomb, inter_si / 2, (side_si + inter_si * 1.5), thick_bld / 2.);
   tr = new TGeoTranslation(shift_bld, 0, thick_bld / 2.);
   quartet->AddNode(blindage, nbl++, tr);
   tr = new TGeoTranslation(-1 * shift_bld, 0, thick_bld / 2.);
   quartet->AddNode(blindage, nbl++, tr);

   fTotSidWBlind = 4 * side_si + 5 * inter_si;
//Coordonnées extraite des côtes données par Yvan M.
//vecteur pointant le milieu d un quartet
//X=-2.231625
//Y=-2.230525
//Z=99.950350
// Mag=100.000139
// Theta=1.808104
// Phi = -135.014124
   TVector3* placement = new TVector3(-2.231625, -2.230525, 99.950350);
   TVector3* Centre = new TVector3();

   TGeoRotation rot1, rot2;
   TGeoTranslation trans;
   TGeoTranslation invZtrans(0, 0, -100);

   TGeoHMatrix h;
   TGeoHMatrix* ph = 0;

//Boucle sur les 4 quartets d un block
   Double_t tx[4] = {1, -1, -1, 1};
   Double_t ty[4] = {1, 1, -1, -1};

   Double_t theta = 0;
   Double_t phi = 0;
   Double_t trans_z = 0;

   for (Int_t nq = 1; nq <= 4; nq += 1) {

      Centre->SetXYZ(placement->X()*tx[nq - 1], placement->Y()*ty[nq - 1], placement->Z());

      theta = Centre->Theta() * TMath::RadToDeg();
      phi = Centre->Phi() * TMath::RadToDeg();

      trans_z = Centre->Mag() + thick_si1 / 2.;

      rot2.SetAngles(phi + 90., theta, 0.);
      rot1.SetAngles(-1.*phi, 0., 0.);
      trans.SetDz(trans_z);

      h = invZtrans * rot2 * trans * rot1;
      ph = new TGeoHMatrix(h);

      AddNode(quartet, nq, ph);
   }
}

void HwwDYBkgEstimate(const string inputFilename, Double_t mHiggs,
                 const string Label) 
{  
  gBenchmark->Start("WWTemplate");

  
  //--------------------------------------------------------------------------------------------------------------
  // Settings 
  //==============================================================================================================
  
  Double_t lumi = 35.5;              // luminosity (pb^-1)
  Int_t ChargeSelection = 0;
//   ChargeSelection = 1;
    
  Double_t fPtMaxLowerCut;
  Double_t fPtMinLowerCut;
  Double_t fDileptonMassUpperCut;
  Double_t fDeltaPhiCut;

  Double_t fHiggsMass = mHiggs;

  //Configure Higgs Mass Dependant Cuts
  if (fHiggsMass == 120) { fPtMaxLowerCut = 20.0;        fPtMinLowerCut = 20.0; 
                           fDileptonMassUpperCut = 40.0; fDeltaPhiCut = 60.0;
  }
  if (fHiggsMass == 130) { fPtMaxLowerCut = 25.0;        fPtMinLowerCut = 20.0;
                           fDileptonMassUpperCut = 45.0; fDeltaPhiCut = 60.0;
  }
  if (fHiggsMass == 140) { fPtMaxLowerCut = 25.0;        fPtMinLowerCut = 20.0;
                           fDileptonMassUpperCut = 45.0; fDeltaPhiCut = 60.0;
  }
  if (fHiggsMass == 150) { fPtMaxLowerCut = 27.0;        fPtMinLowerCut = 25.0;
                           fDileptonMassUpperCut = 50.0; fDeltaPhiCut = 60.0;
  }
  if (fHiggsMass == 160) { fPtMaxLowerCut = 30.0;        fPtMinLowerCut = 25.0;
                           fDileptonMassUpperCut = 50.0; fDeltaPhiCut = 60.0;
  }
  if (fHiggsMass == 170) { fPtMaxLowerCut = 34.0;        fPtMinLowerCut = 25.0;
                           fDileptonMassUpperCut = 50.0; fDeltaPhiCut = 60.0;
  }
  if (fHiggsMass == 180) { fPtMaxLowerCut = 36.0;        fPtMinLowerCut = 25.0;
                           fDileptonMassUpperCut = 60.0; fDeltaPhiCut = 70.0;
  }
  if (fHiggsMass == 190) { fPtMaxLowerCut = 38.0;        fPtMinLowerCut = 25.0;
                           fDileptonMassUpperCut = 80.0; fDeltaPhiCut = 90.0;
  }
  if (fHiggsMass == 200) { fPtMaxLowerCut = 40.0;        fPtMinLowerCut = 25.0;
                           fDileptonMassUpperCut = 90.0; fDeltaPhiCut = 100.0;
  }
  if (fHiggsMass == 210) { fPtMaxLowerCut = 44.0;         fPtMinLowerCut = 25.0;
                           fDileptonMassUpperCut = 110.0; fDeltaPhiCut = 110.0;
  }
  if (fHiggsMass == 220) { fPtMaxLowerCut = 48.0;         fPtMinLowerCut = 25.0;
                           fDileptonMassUpperCut = 120.0; fDeltaPhiCut = 120.0;
  }
  if (fHiggsMass == 230) { fPtMaxLowerCut = 52.0;         fPtMinLowerCut = 25.0;
                           fDileptonMassUpperCut = 130.0; fDeltaPhiCut = 130.0;
  }
  if (fHiggsMass == 250) { fPtMaxLowerCut = 55.0;         fPtMinLowerCut = 25.0;
                           fDileptonMassUpperCut = 150.0; fDeltaPhiCut = 140.0;
  }
  if (fHiggsMass == 300) { fPtMaxLowerCut = 70.0;         fPtMinLowerCut = 25.0;
                           fDileptonMassUpperCut = 200.0; fDeltaPhiCut = 175.0;
  }
  if (fHiggsMass == 350) { fPtMaxLowerCut = 80.0;         fPtMinLowerCut = 25.0;
                           fDileptonMassUpperCut = 250.0; fDeltaPhiCut = 175.0;
  }
  if (fHiggsMass == 400) { fPtMaxLowerCut = 90.0;         fPtMinLowerCut = 25.0;
                           fDileptonMassUpperCut = 300.0; fDeltaPhiCut = 175.0;
  }
  if (fHiggsMass == 450) { fPtMaxLowerCut = 110.0;        fPtMinLowerCut = 25.0;
                           fDileptonMassUpperCut = 350.0; fDeltaPhiCut = 175.0;
  }
  if (fHiggsMass == 500) { fPtMaxLowerCut = 120.0;        fPtMinLowerCut = 25.0;
                           fDileptonMassUpperCut = 400.0; fDeltaPhiCut = 175.0;
  }
  if (fHiggsMass == 550) { fPtMaxLowerCut = 130.0;        fPtMinLowerCut = 25.0;
                           fDileptonMassUpperCut = 450.0; fDeltaPhiCut = 175.0;
  }
  if (fHiggsMass == 600) { fPtMaxLowerCut = 140.0;        fPtMinLowerCut = 25.0;
                           fDileptonMassUpperCut = 500.0; fDeltaPhiCut = 175.0;
  }


  //--------------------------------------------------------------------------------------------------------------
  // Histograms
  //==============================================================================================================  
  TH1D *DileptonMass_allOtherCuts_ee = new TH1D("DileptonMass_allOtherCuts_ee", ";Mass_{ll};Number of Events",150,0.,300.);
  TH1D *DileptonMass_allOtherCuts_emu = new TH1D("DileptonMass_allOtherCuts_emu", ";Mass_{ll};Number of Events",150,0.,300.);
  TH1D *DileptonMass_allOtherCuts_mumu = new TH1D("DileptonMass_allOtherCuts_mumu", ";Mass_{ll};Number of Events",150,0.,300.);
  TH1D *DileptonMass_allOtherCutsExceptMetCut_ee = new TH1D("DileptonMass_allOtherCutsExceptMetCut_ee", ";Mass_{ll};Number of Events",150,0.,300.);
  TH1D *DileptonMass_allOtherCutsExceptMetCut_emu = new TH1D("DileptonMass_allOtherCutsExceptMetCut_emu", ";Mass_{ll};Number of Events",150,0.,300.);
  TH1D *DileptonMass_allOtherCutsExceptMetCut_mumu = new TH1D("DileptonMass_allOtherCutsExceptMetCut_mumu", ";Mass_{ll};Number of Events",150,0.,300.);

  Double_t NEventsIn_BeforeMetCut_ee = 0;
  Double_t NEventsIn_BeforeMetCut_em = 0;
  Double_t NEventsIn_BeforeMetCut_mm = 0;
  Double_t NEventsOut_BeforeMetCut_ee = 0;
  Double_t NEventsOut_BeforeMetCut_em = 0;
  Double_t NEventsOut_BeforeMetCut_mm = 0;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
          map<float, PHG4Hit *> time_sort;
          map<float, PHG4Hit *> layer_sort;
          for (auto & hit : g4hits)
            {
              if (hit)
                {
                  time_sort[hit->get_avg_t()] = hit;
                }
            }

          for (auto & hit_pair : time_sort)
            {

              if (hit_pair.second->get_layer() != UINT_MAX
                  and layer_sort.find(hit_pair.second->get_layer())
                      == layer_sort.end())
                {
                  layer_sort[hit_pair.second->get_layer()] = hit_pair.second;
                }
            }

          if (layer_sort.size() > 5 and mom.Pt() > _pT_threshold) // minimal track length cut
            {

              stringstream spts;

              TVector3 last_pos(0, 0, 0);

              bool first = true;
              for (auto & hit_pair : layer_sort)
                {
                  TVector3 pos(hit_pair.second->get_avg_x(),
                      hit_pair.second->get_avg_y(),
                      hit_pair.second->get_avg_z());

                  // hit step cuts
                  if ((pos - last_pos).Mag() < _min_track_hit_dist
                      and hit_pair.first != (layer_sort.rbegin()->first)
                      and hit_pair.first != (layer_sort.begin()->first))
                    continue;

                  last_pos = pos;

                  if (first)
                    {
                      first = false;
                    }
                  else
                    spts << ",";

                  spts << "[";
                  spts << pos.x();
                  spts << ",";
                  spts << pos.y();
                  spts << ",";
                  spts << pos.z();
                  spts << "]";
                }

              const int abs_pid = abs(g4particle->get_pid());
              int t = 5;
              if (abs_pid
                  == TDatabasePDG::Instance()->GetParticle("pi+")->PdgCode())
                {
                  t = 1;
                }
              else if (abs_pid
                  == TDatabasePDG::Instance()->GetParticle("proton")->PdgCode())
                {
                  t = 2;
                }
              else if (abs_pid
                  == TDatabasePDG::Instance()->GetParticle("K+")->PdgCode())
                {
                  t = 3;
                }
              else if (abs_pid
                  == TDatabasePDG::Instance()->GetParticle("e-")->PdgCode())
                {
                  t = 3;
                }

              const TParticlePDG * pdg_part =
                  TDatabasePDG::Instance()->GetParticle(11);
              const int c =
                  (pdg_part != nullptr) ? (copysign(1, pdg_part->Charge())) : 0;

              fdata
                  << (boost::format(
                      "{ \"pt\": %1%, \"t\": %2%, \"e\": %3%, \"p\": %4%, \"c\": %5%, \"pts\":[ %6% ]}")
                      % mom.Pt() % t % mom.PseudoRapidity() % mom.Phi() % c
                      % spts.str()) << endl;

            }
        }
    }

  fdata.close();
  return 0;
}