本文整理汇总了C++中TVector3::Phi方法的典型用法代码示例。如果您正苦于以下问题:C++ TVector3::Phi方法的具体用法?C++ TVector3::Phi怎么用?C++ TVector3::Phi使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类TVector3
的用法示例。
在下文中一共展示了TVector3::Phi方法的9个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: get_circ_wg_efield_vert
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;
}
示例2: is_in_phi_fiducial
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;
}
示例3: if
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);
}
示例4: Diago
//_________________________________________________________________
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);
}
示例5: HggTreeWriteLoop
//.........这里部分代码省略.........
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;
示例6: execute
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;
//.........这里部分代码省略.........
示例7: mat_si
//.........这里部分代码省略.........
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);
}
}
示例8: HwwDYBkgEstimate
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;
//.........这里部分代码省略.........
示例9: fdata
//.........这里部分代码省略.........
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;
}