c_cpp 代码用muons浏览Pythia事件并使用一些玩具MC测试对背景mehtods

Posted 2021-05-14

//////////////////////////////////////////////////////////
// This class has been automatically generated on
// Tue Jan 31 19:24:36 2017 by ROOT version 6.06/04
// from TTree pythia/Pythia pp 200 GeV
// found on file: out_100.root
//////////////////////////////////////////////////////////

#ifndef pythia_h
#define pythia_h

#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>

// Header file for the classes stored in the TTree if any.
#include "TClonesArray.h"
#include "TObject.h"

class pythia {
public :
   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
   Int_t           fCurrent; //!current Tree number in a TChain

// Fixed size dimensions of array or collections stored in the TTree if any.
   static const Int_t kMaxTracks = 500;

   // Declaration of leaf types
   Int_t           Tracks_;
   UInt_t          Tracks_fUniqueID[kMaxTracks];   //[Tracks_]
   UInt_t          Tracks_fBits[kMaxTracks];   //[Tracks_]
   Int_t           Tracks_mId[kMaxTracks];   //[Tracks_]
   Float_t         Tracks_mEnergy[kMaxTracks];   //[Tracks_]
   Int_t           Tracks_mKF[kMaxTracks];   //[Tracks_]
   Int_t           Tracks_mKS[kMaxTracks];   //[Tracks_]
   Float_t         Tracks_mLifetime[kMaxTracks];   //[Tracks_]
   Float_t         Tracks_mMass[kMaxTracks];   //[Tracks_]
   Int_t           Tracks_mParent[kMaxTracks];   //[Tracks_]
   Float_t         Tracks_mPx[kMaxTracks];   //[Tracks_]
   Float_t         Tracks_mPy[kMaxTracks];   //[Tracks_]
   Float_t         Tracks_mPz[kMaxTracks];   //[Tracks_]
   Float_t         Tracks_mTime[kMaxTracks];   //[Tracks_]
   Float_t         Tracks_mVx[kMaxTracks];   //[Tracks_]
   Float_t         Tracks_mVy[kMaxTracks];   //[Tracks_]
   Float_t         Tracks_mVz[kMaxTracks];   //[Tracks_]

   // List of branches
   TBranch        *b_Tracks_;   //!
   TBranch        *b_Tracks_fUniqueID;   //!
   TBranch        *b_Tracks_fBits;   //!
   TBranch        *b_Tracks_mId;   //!
   TBranch        *b_Tracks_mEnergy;   //!
   TBranch        *b_Tracks_mKF;   //!
   TBranch        *b_Tracks_mKS;   //!
   TBranch        *b_Tracks_mLifetime;   //!
   TBranch        *b_Tracks_mMass;   //!
   TBranch        *b_Tracks_mParent;   //!
   TBranch        *b_Tracks_mPx;   //!
   TBranch        *b_Tracks_mPy;   //!
   TBranch        *b_Tracks_mPz;   //!
   TBranch        *b_Tracks_mTime;   //!
   TBranch        *b_Tracks_mVx;   //!
   TBranch        *b_Tracks_mVy;   //!
   TBranch        *b_Tracks_mVz;   //!

   pythia(TTree *tree=0);
   virtual ~pythia();
   virtual Int_t    Cut(Long64_t entry);
   virtual Int_t    GetEntry(Long64_t entry);
   virtual Long64_t LoadTree(Long64_t entry);
   virtual void     Init(TTree *tree);
   virtual void     Loop();
   virtual Bool_t   Notify();
   virtual void     Show(Long64_t entry = -1);
};

#endif

#ifdef pythia_cxx
pythia::pythia(TTree *tree) : fChain(0) 
{
// if parameter tree is not specified (or zero), connect the file
// used to generate this class and read the Tree.
   if (tree == 0) {
      TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("out_100.root");
      if (!f || !f->IsOpen()) {
         f = new TFile("out_100.root");
      }
      f->GetObject("pythia",tree);

   }
   Init(tree);
}

pythia::~pythia()
{
   if (!fChain) return;
   delete fChain->GetCurrentFile();
}

Int_t pythia::GetEntry(Long64_t entry)
{
// Read contents of entry.
   if (!fChain) return 0;
   return fChain->GetEntry(entry);
}
Long64_t pythia::LoadTree(Long64_t entry)
{
// Set the environment to read one entry
   if (!fChain) return -5;
   Long64_t centry = fChain->LoadTree(entry);
   if (centry < 0) return centry;
   if (fChain->GetTreeNumber() != fCurrent) {
      fCurrent = fChain->GetTreeNumber();
      Notify();
   }
   return centry;
}

void pythia::Init(TTree *tree)
{
   // The Init() function is called when the selector needs to initialize
   // a new tree or chain. Typically here the branch addresses and branch
   // pointers of the tree will be set.
   // It is normally not necessary to make changes to the generated
   // code, but the routine can be extended by the user if needed.
   // Init() will be called many times when running on PROOF
   // (once per file to be processed).

   // Set branch addresses and branch pointers
   if (!tree) return;
   fChain = tree;
   fCurrent = -1;
   fChain->SetMakeClass(1);

   fChain->SetBranchAddress("Tracks", &Tracks_, &b_Tracks_);
   fChain->SetBranchAddress("Tracks.fUniqueID", Tracks_fUniqueID, &b_Tracks_fUniqueID);
   fChain->SetBranchAddress("Tracks.fBits", Tracks_fBits, &b_Tracks_fBits);
   fChain->SetBranchAddress("Tracks.mId", Tracks_mId, &b_Tracks_mId);
   fChain->SetBranchAddress("Tracks.mEnergy", Tracks_mEnergy, &b_Tracks_mEnergy);
   fChain->SetBranchAddress("Tracks.mKF", Tracks_mKF, &b_Tracks_mKF);
   fChain->SetBranchAddress("Tracks.mKS", Tracks_mKS, &b_Tracks_mKS);
   fChain->SetBranchAddress("Tracks.mLifetime", Tracks_mLifetime, &b_Tracks_mLifetime);
   fChain->SetBranchAddress("Tracks.mMass", Tracks_mMass, &b_Tracks_mMass);
   fChain->SetBranchAddress("Tracks.mParent", Tracks_mParent, &b_Tracks_mParent);
   fChain->SetBranchAddress("Tracks.mPx", Tracks_mPx, &b_Tracks_mPx);
   fChain->SetBranchAddress("Tracks.mPy", Tracks_mPy, &b_Tracks_mPy);
   fChain->SetBranchAddress("Tracks.mPz", Tracks_mPz, &b_Tracks_mPz);
   fChain->SetBranchAddress("Tracks.mTime", Tracks_mTime, &b_Tracks_mTime);
   fChain->SetBranchAddress("Tracks.mVx", Tracks_mVx, &b_Tracks_mVx);
   fChain->SetBranchAddress("Tracks.mVy", Tracks_mVy, &b_Tracks_mVy);
   fChain->SetBranchAddress("Tracks.mVz", Tracks_mVz, &b_Tracks_mVz);
   Notify();
}

Bool_t pythia::Notify()
{
   // The Notify() function is called when a new file is opened. This
   // can be either for a new TTree in a TChain or when when a new TTree
   // is started when using PROOF. It is normally not necessary to make changes
   // to the generated code, but the routine can be extended by the
   // user if needed. The return value is currently not used.

   return kTRUE;
}

void pythia::Show(Long64_t entry)
{
// Print contents of entry.
// If entry is not specified, print current entry
   if (!fChain) return;
   fChain->Show(entry);
}
Int_t pythia::Cut(Long64_t entry)
{
// This function may be called from Loop.
// returns  1 if entry is accepted.
// returns -1 otherwise.
   return 1;
}
#endif // #ifdef pythia_cxx

#define pythia_cxx
#include "pythia.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>

void pythia::Loop()
{
//   In a ROOT session, you can do:
//      root> .L pythia.C
//      root> pythia t
//      root> t.GetEntry(12); // Fill t data members with entry number 12
//      root> t.Show();       // Show values of entry 12
//      root> t.Show(16);     // Read and show values of entry 16
//      root> t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntriesFast();

   int totalPip = 0;
   int totalPim = 0;
   int totalKp = 0;
   int totalKm = 0;

   int totalPiMup = 0;
   int totalPiMum = 0;

   int totalPiULS = 0;
   int totalPiLSP  = 0;
   int totalPiLSN  = 0;


   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;


      // cout << "n tracks : " << Tracks_ << endl;

      vector<int> mup;
      vector<int> mum;

      int Npip = 0;
      int Npim = 0;

      int Nkp = 0;
      int Nkm = 0;

      for ( int i = 0; i < Tracks_; i ++ ){

      	TLorentzVector lv;
      	lv.SetXYZM( Tracks_mPx[i], Tracks_mPy[i], Tracks_mPz[i], Tracks_mMass[i] );

      	// cout << "Eta = " << lv.Eta() << endl;
      	if ( lv.Pt() < 0.000001 ) continue;

      	if ( abs( lv.Eta() ) < 0.5  ){
      		if ( Tracks_mKF[i] == 211 ){
      			Npip++;
      			totalPip ++;
      		}
      		if ( Tracks_mKF[i] == -211 ){
      			Npim++;
      			totalPim ++;
      		}

      		if ( Tracks_mKF[i] == 321 ){
      			Nkp ++;
      			totalKp ++;
      		}
      		if ( Tracks_mKF[i] == -321 ){
      			Nkm ++;
      			totalKm ++;
      		}

      		if ( Tracks_mKF[i] == -13 && Tracks_mKF[ Tracks_mParent[ i ] ] == 211){
      			mup.push_back( i );
      			totalPiMup ++;
      		}

      		if ( Tracks_mKF[i] == 13 && Tracks_mKF[ Tracks_mParent[ i ] ] == -211){
      			mum.push_back( i );
      			totalPiMum ++;
      		}
      	}
      }



      TRandom3 rnd;
      rnd.SetSeed(0);

      int Np = rnd.Binomial( Npip, 0.99 );
      int Nm = rnd.Binomial( Npim, 0.99 );

      Np += rnd.Binomial( Nkm*1.2, 0.63 );
      Nm += rnd.Binomial( Nkm, 0.63 );

      // Np += rnd.Binomial( Npip, 0.1 );
      // Nm += rnd.Binomial( Npim, 0.1 );

      Np = rnd.Binomial( Np, 0.05 );
      Nm = rnd.Binomial( Nm, 0.05 );




      if ( Np + Nm == 2 ){
	      totalPiULS += ( Nm * Np );
	      totalPiLSN += ( Nm * (Nm-1)/2.0 ) ;
	      totalPiLSP += ( Np * (Np-1)/2.0 );
  		}


   }


   


   cout << "pi+ / event = " << totalPip / (float)nentries << endl;
   cout << "pi- / event = " << totalPim / (float)nentries << endl;

   cout << "K+ / event = " << totalKp / (float)nentries << endl;
   cout << "K+ / event = " << totalKm / (float)nentries << endl;

   cout << "mu+ (pi+) / event = " << totalPiMup / (float)nentries << endl;
   cout << "mu+ (pi+) / event = " << totalPiMum / (float)nentries << endl;

   cout << "totalPiULS = " << totalPiULS << endl;
   cout << "totalPiLSP = " << totalPiLSP << endl; 
   cout << "totalPiLSN = " << totalPiLSN << endl; 

   cout << "ULS / gmLS = " << totalPiULS / (2*sqrt( totalPiLSN * totalPiLSP )) << endl;

}