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

/*Given the impact point on the front face (vin) and the incoming particle LorentzVector (chi for invisible decay, A' for visible),
 *  determine the interaction point within the fiducial volume and save it in vhit.
 Use a random distribution along the chi flight path, with uniform probability
 This function returns the length (in m) of the trajectory within the fiducial volume.
 Displacement is the lateral displacement (in m) of the detector, along x
 */
double KinUtils::findInteractionPoint(const TLorentzVector &chi, const TVector3 &fiducialV, const TVector3 &vin, TVector3 &vout, TVector3 &vhit) {

	double tz, tx, ty, tout, L;
	int sigPx, sigPy;

	sigPx = (chi.Px() > 0 ? 1 : -1);
	sigPy = (chi.Py() > 0 ? 1 : -1);

	tz = fiducialV.Z() / chi.Pz();
	tx = (sigPx * fiducialV.X() / 2 - vin.X()) / chi.Px();
	ty = (sigPy * fiducialV.Y() / 2 - vin.Y()) / chi.Py();
	tout = 0;

	if ((tz < tx) && (tz < ty)) {
		tout = tz;
	} else if ((tx < ty) && (tx < tz)) {
		tout = tx;
	} else if ((ty < tx) && (ty < tz)) {
		tout = ty;
	}
	vout.SetXYZ(tout * chi.Px() + vin.X(), tout * chi.Py() + vin.Y(), tout * chi.Pz() + vin.Z());
	vhit.SetXYZ(Rand.Uniform(vin.X(), vout.X()), Rand.Uniform(vin.Y(), vout.Y()), Rand.Uniform(vin.Z(), vout.Z()));
	L = (vout - vin).Mag();

	return L;
}

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);
}

TVector3 KVParticle::GetVelocity() const
{
   //returns velocity vector in cm/ns units
   TVector3 beta;
   if (E()) {
      beta = GetMomentum() * (1. / E());
   } else {
      beta.SetXYZ(0, 0, 0);
   }
   return (kSpeedOfLight * beta);
}

//This is the function to determine the interaction point in case of a cylinder (case1), for BDXmini.
//The cylinder is with the axis along y(vertical), center at x=0,y=0,z=ldet, radius is R,height is h
double KinUtils::findInteractionPointCylinder1(const TLorentzVector &chi,double ldet,double h,double R,TVector3 &vin,TVector3 &vout,TVector3 &vhit){

	double tIN,tOUT,t2,t3,t,L;

	double px=chi.Px();
	double py=chi.Py();
	double pz=chi.Pz();

	double delta=pz*pz*ldet*ldet-(pz*pz+px*px)*(ldet*ldet-R*R);
	if (delta<0){
		cout<<"KinUtils::findInteractionPointCylinder1 error, the delta is: "<<delta<<endl;
		return 0;
	}

	//entry point
	tIN = (pz*ldet - sqrt(delta))/(px*px+pz*pz);
	t2 = (pz*ldet + sqrt(delta))/(px*px+pz*pz);

	t3 = (h/2)/py;


	if ((t3>0)&&(t3<t2)&&(t3>tIN)){
		tOUT=t3;
	}else{
		tOUT=t2;
	}

	t=Rand.Uniform(tIN,tOUT);

	vin.SetXYZ(tIN*px,tIN*py,tIN*pz);
	vout.SetXYZ(tOUT*px,tOUT*py,tOUT*pz);

	vhit.SetXYZ(t*px,t*py,t*pz);


	L = (vout - vin).Mag();



	return L;
}

// This is the pt corrected MR - here, we assume the pt
// of the Higgs is (MET+Pt+Qt);
// L1 and L2 are the 4-vectors for the 2 hemispheres, or in you case,
// the two leptons - setting mass to 0 should be fine
// MET is the MET 3 vector (don't forget to set the z-component of
// MET to 0)
// Also, 2 times this variable should give you the Higgs mass
double HWWKinematics::CalcMRNEW(){
  TVector3 vI = MET+L1.Vect()+L2.Vect();
  vI.SetZ(0.0);
  double L1pL2 = CalcMR(); //Note - this calls the old MR function
  double vptx = (L1+L2).Px();
  double vpty = (L1+L2).Py();
  TVector3 vpt;
  vpt.SetXYZ(vptx,vpty,0.0);
  
  float MR2 = 0.5*(L1pL2*L1pL2-vpt.Dot(vI)+L1pL2*sqrt(L1pL2*L1pL2+vI.Dot(vI)-2.*vI.Dot(vpt)));
  
  return sqrt(MR2);
  
}

void tv3Write() {
   //creates the Tree
   TVector3 *v = new TVector3();
   TVector3::Class()->IgnoreTObjectStreamer();
   TFile *f = new TFile("v3.root","recreate");
   TTree *T = new TTree("T","v3 Tree");
   T->Branch("v3","TVector3",&v,32000,1);
   TRandom r;
   for (Int_t i=0;i<10000;i++) {
      v->SetXYZ(r.Gaus(0,1),r.Landau(0,1),r.Gaus(100,10));
      T->Fill();
   }
   T->Write();
   T->Print();
   delete f;
}

//.........这里部分代码省略.........


	if (fDoRoot){
		cout<<"Geant4 recon: start filling histograms"<<endl;
		cout<<"There are: "<<NeventsMC<<" MonteCarlo events"<<endl;
		for (int ii=0;ii<NeventsMC;ii++){
			cMC->GetEntry(ii);

			digi=event->getCollection(OpNoviceDigi::Class(),"DetDigiMC");

			Nhits=digi->GetEntries();
			for (int jj=0;jj<Nhits;jj++){
				faceNumber=((OpNoviceDigi*)digi->At(jj))->GetFaceNumber();
				detNumber=((OpNoviceDigi*)digi->At(jj))->GetDetectorNumber();
				pixelNumber=((OpNoviceDigi*)digi->At(jj))->GetPixelNumber();
				hSimChargeMC[faceNumber][detNumber].at(pixelNumber)->Fill(((OpNoviceDigi*)digi->At(jj))->GetPheCount());
				chargeMC[faceNumber][detNumber].at(pixelNumber)+=((OpNoviceDigi*)digi->At(jj))->GetPheCount();
			}
		}
		cout<<"Geant4 recon: mean"<<endl;
		for (int ii=0;ii<6;ii++){
			for (int jj=0;jj<m_detector->getNdet(ii);jj++){
				if (m_detector->isDetPresent(ii,jj)){
					for (int kk=0;kk<m_detector->getNPixels(ii,jj);kk++){
						qMeanMC=chargeMC[ii][jj].at(kk)/NeventsMC;
						fMC[ii][jj][kk]=qMeanMC;
					}
				}
			}
		}
	}


	vin.SetXYZ(xSource,3.,ySource);
	for (int ii=0;ii<6;ii++){
		for (int jj=0;jj<m_detector->getNdet(ii);jj++){
			if (m_detector->isDetPresent(ii,jj)){
				for (int kk=0;kk<m_detector->getNPixels(ii,jj);kk++){
					F=m_utils->SinglePixelAverageCharge(vin,ii,jj,kk);
					//Corr=CorrectionPixel(vin,ii,jj,kk,m_detector);
				//	cout<<"Setting: "<<F<<" "<<ii<<" "<<jj<<" "<<kk<<endl;
				//	cin.get();
					fTeo[ii][jj][kk]=F;
					if (!fDoRoot) {
						fMC[ii][jj][kk]=0;
						for (int qq=0;qq<1000;qq++){
							hSimChargeMC[ii][jj].at(kk)->Fill(gRandom->Poisson(F));
						}
					}
				}
			}
		}
	}



	cout<<"Comparison, fixing normalization"<<endl;
	/*Now compute normalization and fill histograms*/
	fNormTeo=fNormExp=fNormMC=0;
	for (int iface=0;iface<6;iface++){
		for (int idetector=0;idetector<m_detector->getNdet(iface);idetector++){
			if (m_detector->isDetPresent(iface,idetector)){
				for (int ipixel=0;ipixel<m_detector->getNPixels(iface,idetector);ipixel++){
					fNormTeo+=fTeo[iface][idetector][ipixel];
					fNormMC+=fMC[iface][idetector][ipixel];
					fNormExp+=fExp[iface][idetector][ipixel];

//.........这里部分代码省略.........
  eventTree->SetBranchAddress("Info",       &info);      infoBr       = eventTree->GetBranch("Info");
  eventTree->SetBranchAddress("Electron", &electronArr); electronBr = eventTree->GetBranch("Electron");
  eventTree->SetBranchAddress("Muon", &muonArr);         muonBr = eventTree->GetBranch("Muon");
  eventTree->SetBranchAddress("PFJet", &jetArr);         jetBr = eventTree->GetBranch("PFJet");
  
  for(UInt_t ientry=0; ientry<eventTree->GetEntries(); ientry++) {       	
    infoBr->GetEntry(ientry);
    if (ientry % 100000 == 0) cout << "Event " << ientry << endl;
		
    mithep::RunLumiRangeMap::RunLumiPairType rl(info->runNum, info->lumiSec);      
    if(hasJSON && !rlrm.HasRunLumi(rl)) continue;  // not certified run? Skip to next event...
     
    //for the skimmed input, I already required the HLT bits.
    //    if (!passHLT(info->triggerBits, info->runNum, kTRUE)) continue;

    //********************************************************
    // Load the branches
    //********************************************************
    electronArr->Clear(); 
    muonArr->Clear(); 
    jetArr->Clear(); 
    electronBr->GetEntry(ientry);
    muonBr->GetEntry(ientry);
    jetBr->GetEntry(ientry);

    //event weight
    Double_t eventweight = info->eventweight * lumi;

    //********************************************************
    // TcMet
    //********************************************************
    TVector3 met;        
    if(info->tcMEx!=0 || info->tcMEy!=0) {       
      met.SetXYZ(info->tcMEx, info->tcMEy, 0);
    }

    //********************************************************
    // TcMet
    //********************************************************

    Int_t NSoftMuons = 0;
    Int_t NLeptons = 0;
    vector<Int_t> leptonType;
    vector<Int_t> leptonIndex;
    vector<Double_t> leptonPt;
    vector<Double_t> leptonEta;
    vector<Double_t> leptonPhi;
    vector<Int_t> leptonCharge;

    Int_t NJets = 0;
    const mithep::TJet *leadingJet = 0;
    
    for(Int_t i=0; i<muonArr->GetEntries(); i++) {
      const mithep::TMuon *mu = (mithep::TMuon*)((*muonArr)[i]);      
      if ( (0==0)
           &&
           passMuonCuts(mu)
           &&
           fabs(mu->eta) < 2.4
           && 
           mu->pt > 10.0
        ) {
        leptonPt.push_back(mu->pt);
        leptonEta.push_back(mu->eta);
        leptonPhi.push_back(mu->phi);
        leptonType.push_back(13);

int main(int argc, char* argv[])
{
  TApplication theApp(srcName.Data(), &argc, argv);
//=============================================================================

  if (argc<5) return -1;
  TString sPath = argv[1]; if (sPath.IsNull()) return -1;
  TString sFile = argv[2]; if (sFile.IsNull()) return -1;
  TString sJetR = argv[3]; if (sJetR.IsNull()) return -1;
  TString sSjeR = argv[4]; if (sSjeR.IsNull()) return -1;
//=============================================================================

  sPath.ReplaceAll("#", "/");
//=============================================================================

  double dJetR = -1.;
  if (sJetR=="JetR02") dJetR = 0.2;
  if (sJetR=="JetR03") dJetR = 0.3;
  if (sJetR=="JetR04") dJetR = 0.4;
  if (sJetR=="JetR05") dJetR = 0.5;

  if (dJetR<0.) return -1;
  cout << "Jet R = " << dJetR << endl;
//=============================================================================

  double dSjeR = -1.;
  if (sSjeR=="SjeR01") dSjeR = 0.1;
  if (sSjeR=="SjeR02") dSjeR = 0.2;
  if (sSjeR=="SjeR03") dSjeR = 0.3;
  if (sSjeR=="SjeR04") dSjeR = 0.4;

  if (dSjeR<0.) return -1;
  cout << "Sub-jet R = " << dSjeR << endl;
//=============================================================================

  const double dJetsPtMin  = 0.001;
  const double dCutEtaMax  = 1.6;
  const double dJetEtaMax  = 1.;
  const double dJetAreaRef = TMath::Pi() * dJetR * dJetR;

  fastjet::GhostedAreaSpec areaSpc(dCutEtaMax);
  fastjet::JetDefinition   jetsDef(fastjet::antikt_algorithm, dJetR, fastjet::BIpt_scheme, fastjet::Best);

//fastjet::AreaDefinition  areaDef(fastjet::active_area,areaSpc);
  fastjet::AreaDefinition  areaDef(fastjet::active_area_explicit_ghosts,areaSpc);

//fastjet::JetDefinition   bkgsDef(fastjet::kt_algorithm, 0.2, fastjet::BIpt_scheme, fastjet::Best);
//fastjet::AreaDefinition  aBkgDef(fastjet::active_area_explicit_ghosts, areaSpc);

  fastjet::Selector selectJet = fastjet::SelectorAbsEtaMax(dJetEtaMax);
//fastjet::Selector selectRho = fastjet::SelectorAbsEtaMax(dCutEtaMax-0.2);
//fastjet::Selector selecHard = fastjet::SelectorNHardest(2);
//fastjet::Selector selectBkg = selectRho * (!(selecHard));
//fastjet::JetMedianBackgroundEstimator bkgsEstimator(selectBkg, bkgsDef, aBkgDef);
//fastjet::Subtractor                   bkgSubtractor(&bkgsEstimator);

  fastjet::JetDefinition subjDef(fastjet::antikt_algorithm, dSjeR, fastjet::BIpt_scheme, fastjet::Best);
//=============================================================================

  std::vector<fastjet::PseudoJet> fjInput;
//=============================================================================

  TList *list = new TList();
  TH1D *hWeightSum = new TH1D("hWeightSum", "", 1, 0., 1.); list->Add(hWeightSum);

  TH1D *hJet = new TH1D("hJet", "", 1000, 0., 1000.); hJet->Sumw2(); list->Add(hJet);
  TH2D *hJetNsj = new TH2D("hJetNsj", "", 1000, 0., 1000., 101, -0.5, 100.5); hJetNsj->Sumw2(); list->Add(hJetNsj);

  TH2D *hJetIsj = new TH2D("hJetIsj", "", 1000, 0., 1000., 1000, 0., 1000.); hJetIsj->Sumw2(); list->Add(hJetIsj);
  TH2D *hJet1sj = new TH2D("hJet1sj", "", 1000, 0., 1000., 1000, 0., 1000.); hJet1sj->Sumw2(); list->Add(hJet1sj);
  TH2D *hJet2sj = new TH2D("hJet2sj", "", 1000, 0., 1000., 1000, 0., 1000.); hJet2sj->Sumw2(); list->Add(hJet2sj);
  TH2D *hJetDsj = new TH2D("hJetDsj", "", 1000, 0., 1000., 1000, 0., 1000.); hJetDsj->Sumw2(); list->Add(hJetDsj);

  TH2D *hJetIsz = new TH2D("hJetIsz", "", 1000, 0., 1000., 120, 0., 1.2); hJetIsz->Sumw2(); list->Add(hJetIsz);
  TH2D *hJet1sz = new TH2D("hJet1sz", "", 1000, 0., 1000., 120, 0., 1.2); hJet1sz->Sumw2(); list->Add(hJet1sz);
  TH2D *hJet2sz = new TH2D("hJet2sz", "", 1000, 0., 1000., 120, 0., 1.2); hJet2sz->Sumw2(); list->Add(hJet2sz);
  TH2D *hJetDsz = new TH2D("hJetDsz", "", 1000, 0., 1000., 120, 0., 1.2); hJetDsz->Sumw2(); list->Add(hJetDsz);
//=============================================================================

  HepMC::IO_GenEvent ascii_in(Form("%s/%s.hepmc",sPath.Data(),sFile.Data()), std::ios::in);
  HepMC::GenEvent *evt = ascii_in.read_next_event();

  while (evt) {
    fjInput.resize(0);
    double dXsect  = evt->cross_section()->cross_section() / 1e9;
    double dWeight = evt->weights().back();
    double dNorm = dWeight * dXsect;
    hWeightSum->Fill(0.5, dWeight);

    TVector3 vPar;
    for (HepMC::GenEvent::particle_const_iterator p=evt->particles_begin(); p!=evt->particles_end(); ++p) if ((*p)->status()==1) {
      vPar.SetXYZ((*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz());

      if ((TMath::Abs(vPar.Eta())<dCutEtaMax)) {
        fjInput.push_back(fastjet::PseudoJet(vPar.Px(), vPar.Py(), vPar.Pz(), vPar.Mag()));
      }
    }
//=============================================================================

    fastjet::ClusterSequenceArea clustSeq(fjInput, jetsDef, areaDef);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      pfcandidateArr->Clear(); 
      electronBr->GetEntry(ientry);
      muonBr->GetEntry(ientry);
      photonBr->GetEntry(ientry);
      jetBr->GetEntry(ientry);
      pfcandidateBr->GetEntry(ientry);


      //********************************************************
      // Pileup Energy Density
      //********************************************************
      Double_t rhoEleIso = 0;
      UInt_t EleEAEra = 0;

      if (DataEra == kDataEra_2011_MC) {     
        if (!(isnan(info->RhoKt6PFJetsForIso25) || 
              isinf(info->RhoKt6PFJetsForIso25))) {
          rhoEleIso = info->RhoKt6PFJetsForIso25;
        }
        EleEAEra = kDataEra_2011_Data;
      } else if (DataEra == kDataEra_2012_MC) {
        if (!(isnan(info->RhoKt6PFJets) || 
              isinf(info->RhoKt6PFJets))) {
          rhoEleIso = info->RhoKt6PFJets;
        }
        EleEAEra = kDataEra_2012_Data;
      }

      //********************************************************
      // TcMet
      //********************************************************
      TVector3 pfMet;        
      if(info->pfMEx!=0 || info->pfMEy!=0) {       
        pfMet.SetXYZ(info->pfMEx, info->pfMEy, 0);
      }
      Double_t met = pfMet.Pt();

      Int_t NElectrons = electronArr->GetEntries();
 

      //********************************************************
      // Event Selection Cuts
      //********************************************************
      //veto events with more than 1 reco electron
      if (NElectrons > 1) continue;
      //met cut removed W events
      if (met > 20) continue;


      //******************************************************************************
      //loop over electrons 
      //******************************************************************************
      for(Int_t i=0; i<electronArr->GetEntries(); i++) {
        const higgsana::TElectron *ele = (higgsana::TElectron*)((*electronArr)[i]);
 

        //make cut on dz
        if (fabs(ele->dz) > 0.1) continue;

        //protect against pathologies
        if (TMath::IsNaN(ele->sigiPhiiPhi)) {
          cout << "Pathological SigmaIPhiIPhi : " 
               << info->runNum << " " << info->lumiSec << " " << info->evtNum << endl;
          continue;
        }
        

//.........这里部分代码省略.........
      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 testRotation(){
	

	//LOAD LIBS
	cout << "\n";
	gROOT->Macro("StRoot/LoadLibs.C");
	gSystem->Load("pionPair");
	cout << " loading of pionPair library done" << endl;
	
	
	
	
	TFile* infile = new TFile("/star/u/klandry/ucladisk/2012IFF/schedOutputAll/all_0.root");

	
	//SET UP TREE TO RECEIVE INPUT
	pionPair* pair1 = new pionPair();
	TTree* pairTree = infile->Get("pionPairTree");
	pairTree->SetBranchAddress("pionPair", &pair1);

	
	for (int iPair = 0; iPair < pairTree->GetEntries(); iPair++)
	{
		if (iPair%10000 == 0) {cout << "processing pair number " << iPair << endl;}
		//cout << "processing pair number " << iPair << endl;
		
		
		//if (iPair == 661){continue;}
		
		pairTree->GetEntry(iPair);
		
		TVector3 spinVec;
		
		
		if (pair1->withinRadius(0.05, 0.3))
		{
			
			
			int spinB = pair1->spinBit();
			
			if (spinB == 5 || spinB == 9) //yelow down
			{
				spinVec.SetXYZ(0, -1, 0);
			}
			
			if (spinB == 6 || spinB == 10) //yellow up
			{
				spinVec.SetXYZ(0, 1, 0);
			}
			
			if (spinB == 5 || spinB == 6 || spinB == 9 || spinB == 10)
			{
				
				
				
				TVector3 Ph = pair1->piPlusLV().Vect() + pair1->piMinusLV().Vect();
				TVector3 Rh  = pair1->piPlusLV().Vect() - pair1->piMinusLV().Vect();
				
				
				TVector3 Pa;
				Pa.SetXYZ(0, 0, 1);   //blue is unpolarized beam
				
				TVector3 Pb;
				Pb.SetXYZ(0, 0, -1);  //yellow is polarized beam
				
				
				//ROTATE EVERYTHING BY PI AROUND Y AXIS
				
				Ph.RotateY(TMath::Pi());
				Rh.RotateY(TMath::Pi());
				Pa.RotateY(TMath::Pi());
				Pb.RotateY(TMath::Pi());
				
				
			//	cout << Ph << endl;
			//	cout << Rh << endl;
			//	cout << Pa << endl;
			//	cout << Pb << endl;
				
				
				//cout << "\n";
				
				//cout << Ph.Unit().Cross(Pa) << endl;
				//cout << Ph.Unit().Cross(Rh) << endl;
				
				double cosPhi_S = Pb.Unit().Cross(Ph).Unit() * Pb.Unit().Cross(spinVec).Unit();
				
				double cosPhi_R = Ph.Unit().Cross(Pa).Unit() * Ph.Unit().Cross(Rh).Unit();
				
				double sinPhi_S = Ph.Cross(spinVec) * Pb.Unit() / (Pb.Unit().Cross(Ph).Mag() * Pb.Unit().Cross(spinVec).Mag());
				
				double sinPhi_R = Pa.Cross(Rh) * Ph.Unit() / (Ph.Unit().Cross(Pa).Mag() * Ph.Unit().Cross(Rh).Mag());
				
				
				
				double sinPhi_S_R = sinPhi_S*cosPhi_R - cosPhi_S*sinPhi_R;
				
				double cosPhi_S_R = cosPhi_S*cosPhi_R + sinPhi_S*sinPhi_R;
				
				
//.........这里部分代码省略.........

// Main macro function
//--------------------------------------------------------------------------------------------------
void SkimMultipleFiles_fakes(string inputFilename, string outputFilename) 
{
  gBenchmark->Start("SkimNtuples");
    
  TTree::SetMaxTreeSize(kMaxLong64);
  
  // Don't write TObject part of the objects
  mithep::TEventInfo::Class()->IgnoreTObjectStreamer();
  mithep::TMuon::Class()->IgnoreTObjectStreamer();
  mithep::TElectron::Class()->IgnoreTObjectStreamer();
  mithep::TJet::Class()->IgnoreTObjectStreamer();
  mithep::TPhoton::Class()->IgnoreTObjectStreamer();

  // Data structures to store info from TTrees
  mithep::TEventInfo *info  = new mithep::TEventInfo();
  TClonesArray *muonArr     = new TClonesArray("mithep::TMuon");
  TClonesArray *electronArr     = new TClonesArray("mithep::TElectron");
  TClonesArray *jetArr    = new TClonesArray("mithep::TJet");
  TClonesArray *photonArr     = new TClonesArray("mithep::TPhoton");
    
  UInt_t nInputEvts = 0;
  UInt_t nPassEvts  = 0;
  UInt_t nEventsTotal = 0;

  TFile* outfile = new TFile(outputFilename.c_str(), "RECREATE");
  
  //
  // Initialize data trees and structs
  // 
  TTree *outEventTree = new TTree("Events","Events"); 
  outEventTree->Branch("Info",     &info);
  outEventTree->Branch("Muon",     &muonArr);
  outEventTree->Branch("Electron", &electronArr);
  outEventTree->Branch("PFJet",    &jetArr);
  outEventTree->Branch("Photon",   &photonArr);


  // list input ntuple files to be skimmed
  vector<string> infilenames;  
  ifstream ifs;
  ifs.open(inputFilename.c_str()); 
  assert(ifs.is_open());
  string line;
  while(getline(ifs,line)) { infilenames.push_back(line); }
  ifs.close();

  for(UInt_t ifile=0; ifile<infilenames.size(); ifile++) {
    cout << "Skimming " << infilenames[ifile] << "..." << endl;
    TTree *eventTree = getTreeFromFile(infilenames[ifile].c_str(),"Events"); 
    nEventsTotal += getNEvents(infilenames[ifile].c_str()); 
    assert(eventTree);
    
    // Set branch address to structures that will store the info  
    eventTree->SetBranchAddress("Info",     &info);        TBranch *infoBr     = eventTree->GetBranch("Info");
    eventTree->SetBranchAddress("Muon",     &muonArr);     TBranch *muonBr     = eventTree->GetBranch("Muon");
    eventTree->SetBranchAddress("Electron", &electronArr); TBranch *electronBr = eventTree->GetBranch("Electron");
    eventTree->SetBranchAddress("PFJet",    &jetArr);      TBranch *jetBr      = eventTree->GetBranch("PFJet");
    eventTree->SetBranchAddress("Photon",    &photonArr);  TBranch *photonBr   = eventTree->GetBranch("Photon");
     
    for(UInt_t ientry=0; ientry<eventTree->GetEntries(); ientry++) { 
      infoBr->GetEntry(ientry);
      muonArr->Clear();     muonBr->GetEntry(ientry);
      electronArr->Clear(); electronBr->GetEntry(ientry);
      jetArr->Clear();      jetBr->GetEntry(ientry);
      
      if (ientry % 100000 == 0) cout << "Events: " << ientry << endl;

      nInputEvts++;
      
      //********************************************************
      // TcMet
      //********************************************************
      TVector3 met;        
      if(info->tcMEx!=0 || info->tcMEy!=0) {       
        met.SetXYZ(info->tcMEx, info->tcMEy, 0);
      }
      
      Double_t tempLeadingJetPt = 0;
      for(Int_t i=0; i<jetArr->GetEntries(); i++) {
        const mithep::TJet *jet = (mithep::TJet*)((*jetArr)[i]);
        if( jet->pt > tempLeadingJetPt) tempLeadingJetPt = jet->pt;
      }


      Bool_t keep = kTRUE;
    
      Int_t NRecoMuon = 0;
      Int_t NDenominatorMuon = 0;
      Int_t NMuons = 0;
      for(Int_t i=0; i<muonArr->GetEntries(); i++) {
        const mithep::TMuon *mu = (mithep::TMuon*)((*muonArr)[i]);      
        if (fabs(mu->eta) < 2.4 && mu->pt > 10.0) {
          NRecoMuon++;
          if (passMuonDenominatorCuts(mu)) {
            NDenominatorMuon++;
          }
        }
      }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                    randomSpinBit = 10;
                }


                //CHECK CUTS
                if (sumB.Pt() > ptCutLo && sumB.Pt() < ptCutHi && sumB.M() > mCutLo && sumB.M() < mCutHi && sumB.Eta() > etaCutLo && sumB.Eta() < etaCutHi && useBlueBeam == true)
                {

                    //BLUE BEAM SPIN UP: spin bin 9 and 10
                    if (pair1->spinBit() == 9 || pair1->spinBit() == 10)
                    {
                        bin = hNumberUp->FindBin(pair1->phiSR('b'));
                        //hNumberUp->Fill(pair1->phiSR('b'));



                        polOfBinSumUp[bin] += polarizationOfFill_B[blueFillNo];
                        pErrorOfBinUp[bin] += polErrOfFill_B[blueFillNo];
                    }

                    //BLUE BEAM SPIN DOWN: spin bin 5 and 6
                    if (pair1->spinBit() == 5 || pair1->spinBit() == 6)
                    {
                        bin = hNumberDown->FindBin(pair1->phiSR('b'));
                        hNumberDown->Fill(pair1->phiSR('b'));

                        polOfBinSumDown[bin] += polarizationOfFill_B[blueFillNo];
                        pErrorOfBinDown[bin] += polErrOfFill_B[blueFillNo];
                    }

                }//end blue cuts

                TVector3 Pa;
                Pa.SetXYZ(0, 0, 1);   //blue is unpolarized beam

                TVector3 Pb;
                Pb.SetXYZ(0, 0, -1);  //yellow is polarized beam


                if (sumY.Pt()>ptCutLo && sumY.Pt() < ptCutHi && sumY.M() > mCutLo && sumY.M() < mCutHi && sumY.Eta() > etaCutLo && sumY.Eta() < etaCutHi && useYellowBeam == true)
                {

                    //YELLOW BEAM SPIN UP: spin bin 6 and 10
                    //if (pair1->spinBit() == 6 || pair1->spinBit() == 10)
                    if (randomSpinBit == 6 || randomSpinBit == 10)

                    {

                        totalPairsFinal++;


                        spinVec.SetXYZ(0, 1, 0);

                        TVector3 Ph = pair1->piPlusLV().Vect() + pair1->piMinusLV().Vect();
                        TVector3 Rh  = pair1->piPlusLV().Vect() - pair1->piMinusLV().Vect();


                        double cosPhi_S = -Pb.Unit().Cross(Ph).Unit() * Pb.Unit().Cross(spinVec).Unit();

                        double cosPhi_R = Ph.Unit().Cross(Pb).Unit() * Ph.Unit().Cross(Rh).Unit();

                        double sinPhi_S = Ph.Cross(spinVec) * Pb.Unit() / (Pb.Unit().Cross(Ph).Mag() * Pb.Unit().Cross(spinVec).Mag());

                        double sinPhi_R = Pb.Cross(Rh) * Ph.Unit() / (Ph.Unit().Cross(Pb).Mag() * Ph.Unit().Cross(Rh).Mag());

//.........这里部分代码省略.........
	KIdx=D0d1Idx[D0Idx];
	Pi0Idx=D0d3Idx[D0Idx];
	GammaIdx=Pi0d1Idx[Pi0Idx];

	PitrkIdx=PiTrkIdx[PiIdx];
	SlowPitrkIdx=PiTrkIdx[SlowPiIdx];
	KsPi1trkIdx=PiTrkIdx[KsPi1Idx];
	KsPi2trkIdx=PiTrkIdx[KsPi2Idx];


	if(abs(Dstard1Lund[DstarIdx]-D0Lund[D0Idx])>0)HLundCheck.Fill(1);
	else if(abs(Dstard2Lund[DstarIdx]-PiLund[SlowPiIdx])>0)HLundCheck.Fill(2);
	else if(abs(D0d1Lund[D0Idx]-KLund[KIdx])>0)HLundCheck.Fill(3);
	else if(abs(D0d2Lund[D0Idx]-PiLund[PiIdx])>0)HLundCheck.Fill(4);
	else if(abs(TRKLund[SlowPitrkIdx]-PiLund[SlowPiIdx])>0)HLundCheck.Fill(5);
	else if(abs(TRKLund[PitrkIdx]-PiLund[PiIdx])>0)HLundCheck.Fill(6);
	else if(abs(Ds1pd1Lund[Ds1pIdx]-DstarLund[DstarIdx])>0)HLundCheck.Fill(7);
	else if(abs(Ds1pd2Lund[Ds1pIdx]-KsLund[KsIdx])>0)HLundCheck.Fill(8);
	else if(abs(TRKLund[KsPi1trkIdx]-PiLund[KsPi1Idx])>0)HLundCheck.Fill(9);
	else if(abs(TRKLund[KsPi2trkIdx]-PiLund[KsPi2Idx])>0)HLundCheck.Fill(10);
	else if(abs(Ksd1Lund[KsIdx]-PiLund[KsPi1Idx])>0)HLundCheck.Fill(11);
	else if(abs(Ksd2Lund[KsIdx]-PiLund[KsPi2Idx])>0)HLundCheck.Fill(12);
	else if(abs(D0d3Lund[D0Idx]-Pi0Lund[Pi0Idx])>0)HLundCheck.Fill(13);
	else if(abs(Pi0d1Lund[Pi0Idx]-GammaLund[GammaIdx])>0)HLundCheck.Fill(14);	
	else HLundCheck.Fill(0);

      }


      ///compute some quantities
      D0Probab=TMath::Prob(D0Chi2[D0Idx],D0nDof[D0Idx]);     
      KsProbab=TMath::Prob(KsChi2[D0Idx],KsnDof[D0Idx]);     
         
      Ksp3Direction.SetXYZ(sin(acos(Kscosth[KsIdx]))*cos(Ksphi[KIdx]),
			   sin(acos(Kscosth[KsIdx]))*sin(Ksphi[KIdx]),
			   Kscosth[KsIdx]); 
      KsFlightDirection.SetXYZ(KsVtxx[KsIdx]-Ds1pVtxx[Ds1pIdx],
			       KsVtxy[KsIdx]-Ds1pVtxy[Ds1pIdx],
			       KsVtxz[KsIdx]-Ds1pVtxz[Ds1pIdx]);
      KsCosine=Ksp3Direction*KsFlightDirection/KsFlightDirection.Mag();

      ///Fill Distributions  
      if(Mode=="D0ToKPiPi0"){        
	HPi0MassPreCuts.Fill(Pi0Mass[Pi0Idx]);  
      }

      HD0p3CM.Fill(D0p3CM[D0Idx]);
      HD0MassPreCuts.Fill(D0Mass[D0Idx]);
      HD0Probab.Fill(log(D0Probab));

      HDstarMassPreCuts.Fill(DstarMass[DstarIdx]);   

      HKsMassPreCuts.Fill(KsMass[KsIdx]);
      HKsProbab.Fill(log(KsProbab));
      HKsCosine.Fill(KsCosine);

      HDs1pMassPreCuts.Fill(Ds1pMass[Ds1pIdx]);  


      //Apply Cuts       
      //     if((PiSelectorsMap[KsPi1trkIdx] & (1<<4) ) != 0 && (PiSelectorsMap[KsPi2trkIdx] & (1<<4) ) != 0){
      //	HDs1pMassPID.Fill(Ds1pMass[Ds1pIdx]);

      if(D0p3CM[D0Idx] > D0p3CMCut){
	if(D0Probab > D0ProbabCut){
	  if(KsProbab > KsProbabCut && KsCosine>KsCosineCut){