C++ TTree::Print方法代码示例

int main(int argc, char *argv[]){

  TFile* fout = new TFile("cascade.result.root","RECREATE");

  HepMC::IO_GenEvent fin("../CASCADE_HEPMC/bin/example_test_out.dat",std::ios::in);
  
  HepMC::GenEvent* evt = new HepMC::GenEvent;

  TTree *tree = new TTree("tree","cascade result");
  tree->Branch("event",&evt);

  Int_t nevt = 0;

  while (evt = fin.read_next_event() ) {

    nevt++;

    if(debug){ 
      cout<<"nevt: "<<nevt<<endl;
      getchar();
      evt->print();
    }

    tree->Fill();
  }

  tree->Print();

  fout->cd();
  tree->Write();
  fout->Close();
  
  delete fout;

  return(0);  
}

void CreateParentTree() {
   // create a simple TTree with 5 branches
   // Two branches ("Run" and "Event") will be used to index the Tree
   TFile *f = new TFile("treeparent.root","recreate");
   TTree *T = new TTree("T","test friend trees");
   T->Branch("Run",&Run,"Run/I");
   T->Branch("Event",&Event,"Event/I");
   T->Branch("x",&x,"x/F");
   T->Branch("y",&y,"y/F");
   T->Branch("z",&z,"z/F");
   TRandom r;
   for (Int_t i=0;i<10000;i++) {
      if (i < 5000) Run = 1;
      else          Run = 2;
      Event = i;
      x = r.Gaus(10,1);
      y = r.Gaus(20,2);
      z = r.Landau(2,1);
      T->Fill();
   }
   T->Print();
   T->Write();
   delete f;
}

// nEvents is how many events we want.
int makeEventSample(Int_t nEvents)
{
  // Load needed libraries
  loadLibraries();

  // Create an instance of the Pythia event generator ...
  TPythia6* pythia = new TPythia6;

  // ... and initialise it to run p+p at sqrt(200) GeV in CMS
  pythia->Initialize("cms", "p", "p", 200);

  // Open an output file
  TFile* file = TFile::Open(FILENAME, "RECREATE");
  if (!file || !file->IsOpen()) {
    Error("makeEventSample", "Couldn;t open file %s", FILENAME);
    return 1;
  }

  // Make a tree in that file ...
  TTree* tree = new TTree(TREENAME, "Pythia 6 tree");

  // ... and register a the cache of pythia on a branch (It's a
  // TClonesArray of TMCParticle objects. )
  TClonesArray* particles = (TClonesArray*)pythia->GetListOfParticles();
  tree->Branch(BRANCHNAME, &particles);

  // Now we make some events
  for (Int_t i = 0; i < nEvents; i++) {
    // Show how far we got every 100'th event.
    if (i % 100 == 0)
      cout << "Event # " << i << endl;

    // Make one event.
    pythia->GenerateEvent();

    // Maybe you want to have another branch with global event
    // information.  In that case, you should process that here.
    // You can also filter out particles here if you want.

    // Now we're ready to fill the tree, and the event is over.
    tree->Fill();
  }

  // Show tree structure
  tree->Print();

  // After the run is over, we may want to do some summary plots:
  TH1D* hist = new TH1D(HISTNAME, "p_{#perp}  spectrum for  #pi^{+}",
                        100, 0, 3);
  hist->SetXTitle("p_{#perp}");
  hist->SetYTitle("dN/dp_{#perp}");
  char expression[64];
  sprintf(expression,"sqrt(pow(%s.fPx,2)+pow(%s.fPy,2))>>%s",
          BRANCHNAME, BRANCHNAME, HISTNAME);
  char selection[64];
  sprintf(selection,"%s.fKF==%d", BRANCHNAME, PDGNUMBER);
  tree->Draw(expression,selection);

  // Normalise to the number of events, and the bin sizes.
  hist->Sumw2();
  hist->Scale(3 / 100. / hist->Integral());
  hist->Fit("expo", "QO+", "", .25, 1.75);
  TF1* func = hist->GetFunction("expo");
  func->SetParNames("A", "- 1 / T");
  // and now we flush and close the file
  file->Write();
  file->Close();

  return 0;
}

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

  if      (whichSkims == 0) path += "merged/";
  else if (whichSkims == 1) path += "minimal/";
  else if (whichSkims == 2) path += "complete/";
  else {
    std::cout << "NO SKIMS TYPE DEFINED!\n";
    return;
  }



  if      ( dataset == 0 ) path += "DYJetsToLL_";
  else if ( dataset == 1 ) path += "DYToMuMu_";
  else if ( dataset == 2 ) path += "TTJets_";
  else if ( dataset == 3 ) path += "DYToTauTau_";
  else if ( dataset == 4 ) path += "WJetsToLNu_";
  else if ( dataset == 5 ) path += "WW_";
  else if ( dataset == 6 ) path += "WZ_";
  else if ( dataset == 7 ) path += "ZZ_";
  else if ( dataset == 8 ) path += "QCD_Pt_20_MuEnrichedPt_15_";
  else if ( dataset == 9 ) path += "DY2JetsToLL_";


  else if ( dataset == 9990 ) path += "DYJetsToLL_PtZ-50To70_";
  else if ( dataset == 9990 ) path += "DYJetsToLL_PtZ-70To100_";
  else if ( dataset == 9992 ) path += "DYJetsToLL_PtZ-100_";


  else if ( dataset == 10 ) path += "SingleMuRun2012A-13Jul2012-v1_";
  else if ( dataset == 11 ) path += "SingleMuRun2012A-recover-06Aug2012-v1_";
  else if ( dataset == 12 ) path += "SingleMuRun2012B-13Jul2012-v1_";
  else if ( dataset == 13 ) path += "SingleMuRun2012C-24Aug2012-v1_";
  else if ( dataset == 14 ) path += "SingleMuRun2012C-PromptReco-v2_";
  else if ( dataset == 15 ) path += "SingleMuRun2012D-PromptReco-v1_";


  else if ( dataset == 20 ) path += "DoubleMuRun2012A-13Jul2012-v1_";
  else if ( dataset == 21 ) path += "DoubleMuRun2012B-13Jul2012-v4_";
  else if ( dataset == 22 ) path += "DoubleMuRun2012C-24Aug2012-v1_";
  else if ( dataset == 23 ) path += "DoubleMuRun2012C-PromptReco-v2_";

  else {
    std::cout << "NO DATASET DEFINED!\n";
    return;
  }


  if      (whichSkims == 0) newfile = new TFile(path+"merged.root","recreate");
  else if (whichSkims == 1) newfile = new TFile(path+"minimal.root","recreate");
  else if (whichSkims == 2) newfile = new TFile(path+"complete.root","recreate");
  else {
    std::cout << "NO SKIMS TYPE DEFINED!\n";
    return;
  }


  std::cout << "saving in file " << newfile -> GetName() << std::endl;
  TTree *newtree = tree->CloneTree(0);


  // ======== GET THE HANDLES FOR SKIMMING ========
  float recoCandMass;
  tree->SetBranchAddress("recoCandMass", &recoCandMass);

  float vertexNormChiSquare, angleDiMuons;

  tree->SetBranchAddress("vertexNormChiSquare", &vertexNormChiSquare);
  tree->SetBranchAddress("angleDiMuons"       , &angleDiMuons );
  
  _MuonInfo reco1, reco2;
  
  tree->SetBranchAddress("reco1", &reco1);
  tree->SetBranchAddress("reco2", &reco2);


  // ======== PERFORM THE SKIMMING ========
  cout<<"Loop over the " << tree->GetEntries() << " entries ...\n";
  for (int iEvt=0; iEvt < tree->GetEntries(); iEvt++) {
    
    if ( (iEvt % 500000)==0 ) cout << "event " << iEvt << endl;
    tree -> GetEntry(iEvt);

    // additional selection cuts
    if (recoCandMass <  60 && whichSkims > 0) continue;
    //if (recoCandMass > 160 && whichSkims > 0) continue;
    if (vertexNormChiSquare > 10 && whichSkims > 0) continue; 
    if (angleDiMuons > TMath::Pi()-0.02 && whichSkims > 0) continue;

    if (!isKinTight_2012_noAcc(reco1) && whichSkims > 1) continue;
    if (!isKinTight_2012_noAcc(reco2) && whichSkims > 1) continue;

    newtree->Fill(); 
  }

  newtree->Print();
  newtree->AutoSave();
  
  std::cout << "new tree has " << newtree -> GetEntries() << std::endl;
  delete newfile;
}

//.........这里部分代码省略.........
       float value = histoRecoDelta_yVsTTPt->GetBinContent(bin_index1,bin_index2)/histoGenDelta_yVsTTPt->GetBinContent(bin_index1,bin_index2);
       histoDelta_yVsTTPtEfficiency->SetBinContent(bin_index1, bin_index2, value);
       float value_fb = histoFineBinning_RecoDelta_yVsTTPt->GetBinContent(bin_index1,bin_index2)/histoFineBinning_GenDelta_yVsTTPt->GetBinContent(bin_index1,bin_index2);
       histoFineBinning_Delta_yVsTTPtEfficiency->SetBinContent(bin_index1, bin_index2, value_fb);

     }
    }


    for (unsigned int bin_index1 = 1; bin_index1 < 11; bin_index1++) {
     for (unsigned int bin_index2 = 1; bin_index2 < 6; bin_index2++) {
       float value = histoRecoDelta_yVsTTMass->GetBinContent(bin_index1,bin_index2)/histoGenDelta_yVsTTMass->GetBinContent(bin_index1,bin_index2);
       histoDelta_yVsTTMassEfficiency->SetBinContent(bin_index1, bin_index2, value);
       float value_fb = histoFineBinning_RecoDelta_yVsTTPt->GetBinContent(bin_index1,bin_index2)/histoFineBinning_GenDelta_yVsTTPt->GetBinContent(bin_index1,bin_index2);
       histoFineBinning_Delta_yVsTTPtEfficiency->SetBinContent(bin_index1, bin_index2, value_fb);

     }
    }






  }				// end of loop over the datasets


//  if(verbosity>0) {
    cout<<"#########################"<<endl;
    cout<<"    End of the program   "<<endl;
    cout<<"#########################"<<endl;
//  }

    tree->Print();

    TupleFile->Write();
    TupleFile->Close();


    TFile * fileCA = new TFile("ChargeAsymmetryResults_Skimmed_TTbarSemileptonicPowheg.root","RECREATE");
    //  fileCA->cd();

    // GEN
    histoGenTopRapidity1->Write();
    histoGenTopRapidity2->Write();
    histoGenTopMass1->Write();
    histoGenTopMass2->Write();
    histoGenTopMass1VsTopMass2->Write();
    histoGenDelta_y->Write();
    histoGenTTRapidity->Write();
    histoGenTTPt->Write();
    histoGenTTMass->Write();
    histoGenDelta_yVsTTRapidity->Write();
    histoGenDelta_yVsTTPt->Write();
    histoGenDelta_yVsTTMass->Write();
    histoGenN_plusTTRapidity->Write();
    histoGenN_plusTTPt->Write();
    histoGenN_plusTTMass->Write();
    histoGenN_minusTTRapidity->Write();
    histoGenN_minusTTPt->Write();
    histoGenN_minusTTMass->Write();
    histoGenN_plusMinusN_minusTTRapidity->Write();
    histoGenN_plusPlusN_minusTTRapidity->Write();
    histoGenN_plusMinusN_minusTTPt->Write();
    histoGenN_plusPlusN_minusTTPt->Write();
    histoGenN_plusMinusN_minusTTMass->Write();

//.........这里部分代码省略.........
          oldtree->SetBranchStatus("*", 0);
          oldtree->SetBranchStatus("runnumber", 1);
          oldtree->SetBranchStatus("eventnumber", 1);
          oldtree->SetBranchStatus("luminosityBlock", 1);
          oldtree->SetBranchStatus("HLT_Mu22_Photon22_CaloIdL", 1);
          oldtree->SetBranchStatus("HLT_Mu8_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL", 1);
          oldtree->SetBranchStatus("HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL", 1);
          oldtree->SetBranchStatus("prescale_HLT_Mu22_Photon22_CaloIdL", 1);
          oldtree->SetBranchStatus("prescale_HLT_Mu8_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL", 1);
          oldtree->SetBranchStatus("prescale_HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL", 1);
          oldtree->SetBranchStatus("rho", 1);
          oldtree->SetBranchStatus("pfmet", 1);
          oldtree->SetBranchStatus("pvsize", 1);
          oldtree->SetBranchStatus("JetColl_size", 1);
          oldtree->SetBranchStatus("Jet_pt", 1);
          oldtree->SetBranchStatus("muon_size", 1);
          oldtree->SetBranchStatus("muon_pt", 1);
          oldtree->SetBranchStatus("muon_ptError", 1);
          oldtree->SetBranchStatus("muon_eta", 1);
          oldtree->SetBranchStatus("muon_phi", 1);
          oldtree->SetBranchStatus("muon_charge", 1);
          oldtree->SetBranchStatus("muon_nhitspixel", 1);
          oldtree->SetBranchStatus("muon_nhitstrack", 1);
          oldtree->SetBranchStatus("muon_nhitsmuons", 1);
          oldtree->SetBranchStatus("muon_nlayerswithhits", 1);
          oldtree->SetBranchStatus("muon_nSegmentMatch", 1);
          oldtree->SetBranchStatus("muon_isTrackerMuon", 1);
          oldtree->SetBranchStatus("muon_normChi2", 1);
          oldtree->SetBranchStatus("muon_dz_beamSpot", 1);
          oldtree->SetBranchStatus("muon_dz_firstPVtx", 1);
          oldtree->SetBranchStatus("muon_dxy_cmsCenter", 1);
          oldtree->SetBranchStatus("muon_dxy_beamSpot", 1);
          oldtree->SetBranchStatus("muon_dxy_firstPVtx", 1);
          oldtree->SetBranchStatus("muon_trackIso03", 1);
          oldtree->SetBranchStatus("muon_emIso03", 1);
          oldtree->SetBranchStatus("muon_hadIso03", 1);
          oldtree->SetBranchStatus("gsf_size", 1);
          oldtree->SetBranchStatus("gsf_eta", 1);
          oldtree->SetBranchStatus("gsf_phi", 1);
          oldtree->SetBranchStatus("gsf_theta", 1);
          oldtree->SetBranchStatus("gsf_charge", 1);
          oldtree->SetBranchStatus("gsf_sigmaetaeta", 1);
          oldtree->SetBranchStatus("gsf_sigmaIetaIeta", 1);
          oldtree->SetBranchStatus("gsf_dxy_firstPVtx", 1);
          oldtree->SetBranchStatus("gsf_dz_beamSpot", 1);
          oldtree->SetBranchStatus("gsf_dz_firstPVtx", 1);
          oldtree->SetBranchStatus("gsf_nLostInnerHits", 1);
          oldtree->SetBranchStatus("gsf_deltaeta", 1);
          oldtree->SetBranchStatus("gsf_deltaphi", 1);
          oldtree->SetBranchStatus("gsf_hovere", 1);
          oldtree->SetBranchStatus("gsf_trackiso", 1);
          oldtree->SetBranchStatus("gsf_ecaliso", 1);
          oldtree->SetBranchStatus("gsf_hcaliso1", 1);
          oldtree->SetBranchStatus("gsf_hcaliso2", 1);
          oldtree->SetBranchStatus("gsf_isecaldriven", 1);
          oldtree->SetBranchStatus("gsfsc_e", 1);
          oldtree->SetBranchStatus("gsfsc_eta", 1);
          oldtree->SetBranchStatus("gsfsc_phi", 1);
          oldtree->SetBranchStatus("gsf_e2x5overe5x5", 1);
          oldtree->SetBranchStatus("gsf_e1x5overe5x5", 1);
          oldtree->SetBranchStatus("gsf_gsfet", 1);
          oldtree->SetBranchStatus("genPair_mass", 1);
          oldtree->SetBranchStatus("trueNVtx", 1);

  
	  //Create a new file + a clone of old tree header. Do not copy events
	  TFile *newfile = new TFile(outputfilepath,"recreate");
	  newfile->mkdir("gsfcheckerjob");
	  newfile->cd("gsfcheckerjob");
	  TTree *newtree = oldtree->CloneTree(0);
	  
	  
	  for (Long64_t i = 0; i < nentries; ++i) {
	    if (i %100000 == 0) cout << "entry nb: " <<i << endl;
	    oldtree->GetEntry(i);
	    if(gsf_size < 1) continue;
	    if(muon_size < 1) continue;
	    bool passEleSelection = false;
	    bool passMuSelection = false;
	    for (int it = 0; it < gsf_size; ++it) {
	      if (gsf_gsfet[it] < 35.) continue;
		passEleSelection = true;
	    }
	    for (int muIt = 0; muIt < muon_size; ++muIt) {
	      if (muon_pt[muIt] < 35.) continue;
		passMuSelection = true;
	    }
	    if (passEleSelection && passMuSelection) newtree->Fill();
     
	  }


	  newtree->Print();
	  newfile->Write();
	  delete oldfile;
	  delete newfile;
	}
      myfile.close();
    } 
}

int stdhep2root(int nEvent=100) {

  gROOT->Reset();

  // In order to refer to shareable libraries in this simple form,
  // include their paths in a .rootrc file.
  gSystem->Load("libLCDEvent");
  gSystem->Load("libLCDRootAppsUtil");
  gSystem->Load("libLCDFastMC");
  gSystem->Load("libLCDStdHEPUtil");
    gSystem->Load("libEG");
  gSystem->Load("libEGPythia6");
  gSystem->Load("libLCDEvent");
  gSystem->Load("libLCDRootAppsUtil");
  gSystem->Load("libLCDFastMC");
  gSystem->Load("libLCDPhUtil");
  gSystem->Load("libLCDGenUtil");


  TString parfile_dir;
  parfile_dir += gSystem->Getenv("LCDROOT")   ;parfile_dir += "/";
  parfile_dir += "ParFiles/";

  LCDEvent *event= new LCDEvent(); //Create Physics event container first.

  // Input stdhep data
  Char_t* stdFile = "test.hep";
  // Char_t* stdFile = "../../../GenData/eetobbZ0.dat";
  LCDreadStdFile source(stdFile,event);

  //
  // Open root file.
  //
  Int_t comp=2;
  Char_t* rootFile = "eetobbZ0.root";
  TFile* hfile = new TFile(rootFile,"RECREATE","LCD Event ROOT file");
  hfile->SetCompressionLevel(comp);

  //
  // Create a tree with one superbranch.
  //
  TTree* tree = new TTree( "T", "LCD ROOT tree" );
  tree->SetAutoSave( 1000000000 );
  tree->Branch( "LCDEvent", "LCDEvent", &event);
  
  // Event loop
  Int_t iEvent;
  for (iEvent = 0; iEvent < nEvent; iEvent++) {
    if (!source.GetEvent()) break; // Read an event from the stdhep file.
    tree->Fill();
  }

  //
  // Clean up the ROOT file.
  //
  hfile->Write();
  tree->Print();
  hfile->Close();

  return iEvent;
}

void CutTree(){

	int NMin=300000;

/*	char outputfilename[200];
	sprintf(outputfilename,"Results/ProblemBin.txt");
	printf("output filename is: %s\n", outputfilename);
	FILE *outputFile = fopen(outputfilename,"w");
*/
	TFile* fInData = new TFile(fileNameData);

	TTree *treeData = (TTree*)fInData->Get("data");

	treeData->Print();


	TTree *fChain;
	fChain=treeData;

	Long64_t nentries = fChain->GetEntries();
	cout<<nentries<<endl;
	Long64_t nb = 0;

	Double_t        JpsiMass;
	Double_t        JpsiPt;
	Double_t        JpsiRap;
	Double_t        Jpsict;
	Int_t        JpsiType_idx;
	Int_t        MCType_idx;
	Double_t        costh_CS;
	Double_t        phi_CS;
	Double_t        costh_HX;
	Double_t        phi_HX;

	TBranch        *b_JpsiMass;   //!
	TBranch        *b_JpsiPt;   //!
	TBranch        *b_JpsiRap;   //!
	TBranch        *b_Jpsict;   //!
	TBranch        *b_JpsiType_idx;   //!
	TBranch        *b_MCType_idx;   //!
	TBranch        *b_costh_CS;   //!
	TBranch        *b_phi_CS;   //!
	TBranch        *b_costh_HX;   //!
	TBranch        *b_phi_HX;   //!

	fChain->SetBranchAddress("JpsiMass", &JpsiMass, &b_JpsiMass);
	fChain->SetBranchAddress("JpsiPt", &JpsiPt, &b_JpsiPt);
	fChain->SetBranchAddress("JpsiRap", &JpsiRap, &b_JpsiRap);
	fChain->SetBranchAddress("Jpsict", &Jpsict, &b_Jpsict);
	fChain->SetBranchAddress("JpsiType_idx", &JpsiType_idx, &b_JpsiType_idx);
	fChain->SetBranchAddress("MCType_idx", &MCType_idx, &b_MCType_idx);
	fChain->SetBranchAddress("costh_CS", &costh_CS, &b_costh_CS);
	fChain->SetBranchAddress("phi_CS", &phi_CS, &b_phi_CS);
	fChain->SetBranchAddress("costh_HX", &costh_HX, &b_costh_HX);
	fChain->SetBranchAddress("phi_HX", &phi_HX, &b_phi_HX);


//	Long64_t jentry=NMin;

    TNtuple* pseudodata = new TNtuple("data","data","JpsiMass:JpsiPt:JpsiRap:Jpsict:costh_CS:phi_CS:costh_HX:phi_HX:JpsiType_idx:MCType_idx");
    TNtuple* data = new TNtuple("data","data","JpsiMass:JpsiPt:JpsiRap:Jpsict:costh_CS:phi_CS:costh_HX:phi_HX:JpsiType_idx:MCType_idx");

    for (Long64_t jentry=0; jentry<NMin;jentry++) {
       		    	 if(jentry % 1000 == 0) printf("event %d\n", (Int_t) jentry);

       	   nb = fChain->GetEntry(jentry);   //nbytes += nb;

       	pseudodata->Fill(JpsiMass,JpsiPt,JpsiRap,Jpsict,costh_CS,phi_CS,costh_HX,phi_HX,JpsiType_idx,MCType_idx);

        }

    for (Long64_t jentry=0; jentry<nentries;jentry++) {
   		    	 if(jentry % 1000 == 0) printf("event %d\n", (Int_t) jentry);

   	   nb = fChain->GetEntry(jentry);   //nbytes += nb;

// 	  fprintf(outputFile, "%f\n",Jpsict);


   	data->Fill(JpsiMass,JpsiPt,JpsiRap,Jpsict,costh_CS,phi_CS,costh_HX,phi_HX,JpsiType_idx,MCType_idx);

    }

    pseudodata->Print();

    pseudodata->SaveAs("/scratch/knuenz/Polarization/RootInput/TNtuple_red_PR_pseudo.root");

	data->Print();

	data->SaveAs("/scratch/knuenz/Polarization/RootInput/TNtuple_red_PR.root");

//	RooDataSetData.Print();


//	fclose(outputFile);


//	treeData.Print();

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

TFile *cernbuild(Int_t getFile=0, Int_t print=1) {

   Int_t           Category;
   UInt_t          Flag;
   Int_t           Age;
   Int_t           Service;
   Int_t           Children;
   Int_t           Grade;
   Int_t           Step;
   Int_t           Hrweek;
   Int_t           Cost;
   Char_t          Division[4];
   Char_t          Nation[3];

   //The input file cern.dat is a copy of the CERN staff data base
   //from 1988
   TString filename = "cernstaff.root";
   TString dir = gROOT->GetTutorialDir();
   dir.Append("/tree/");
   dir.ReplaceAll("/./","/");
   FILE *fp = fopen(Form("%scernstaff.dat",dir.Data()),"r");

   TFile *hfile = 0;
   if (getFile) {
      // if the argument getFile =1 return the file "cernstaff.root"
      // if the file does not exist, it is created
      if (!gSystem->AccessPathName(dir+"cernstaff.root",kFileExists)) {
         hfile = TFile::Open(dir+"cernstaff.root"); //in $ROOTSYS/tutorials/tree
         if (hfile) return hfile;
      }
      //otherwise try $PWD/cernstaff.root
      if (!gSystem->AccessPathName("cernstaff.root",kFileExists)) {
         hfile = TFile::Open("cernstaff.root"); //in current dir
         if (hfile) return hfile;
      }
   }
   //no cernstaff.root file found. Must generate it !
   //generate cernstaff.root in $ROOTSYS/tutorials/tree if we have write access
   if (gSystem->AccessPathName(".",kWritePermission)) {
      printf("you must run the script in a directory with write access\n");
      return 0;
   }
   hfile = TFile::Open(filename,"RECREATE");
   TTree *tree = new TTree("T","CERN 1988 staff data");
   tree->Branch("Category",&Category,"Category/I");
   tree->Branch("Flag",&Flag,"Flag/i");
   tree->Branch("Age",&Age,"Age/I");
   tree->Branch("Service",&Service,"Service/I");
   tree->Branch("Children",&Children,"Children/I");
   tree->Branch("Grade",&Grade,"Grade/I");
   tree->Branch("Step",&Step,"Step/I");
   tree->Branch("Hrweek",&Hrweek,"Hrweek/I");
   tree->Branch("Cost",&Cost,"Cost/I");
   tree->Branch("Division",Division,"Division/C");
   tree->Branch("Nation",Nation,"Nation/C");
   char line[80];
   while (fgets(line,80,fp)) {
      sscanf(&line[0],"%d %d %d %d %d %d %d  %d %d %s %s",
      &Category,&Flag,&Age,&Service,&Children,&Grade,&Step,&Hrweek,&Cost,Division,Nation);
      tree->Fill();
   }
   if (print) tree->Print();
   tree->Write();

   fclose(fp);
   delete hfile;
   if (getFile) {
      //we come here when the script is executed outside $ROOTSYS/tutorials/tree
      hfile = TFile::Open(filename);
      return hfile;
   }
   return 0;
}

void create_tree_for_toyMC()()
{
  gROOT->Reset();

  Double_t Y;
  Double_t Y_true;
  Double_t dY;
  Double_t Tk;
  Double_t Tk_true;
  Double_t dTk;
  Double_t Sa;
  Double_t Sa_true;
  Double_t dSa;
  Double_t Iso;
  Double_t Iso_true;
  Double_t dIso;
  Double_t Hlt;
  Double_t Hlt_true;
  Double_t dHlt;
  Double_t chi2;

  TFile *f;
  TTree *tree;
  
  f = new TFile("fitResult.root","RECREATE");
  tree = new TTree("tree"," C data from ASCII file");

  tree->Branch("Y",&Y,"Y/D");
  tree->Branch("Y_true",&Y_true,"Y/D");
  tree->Branch("dY",&dY,"dY/D");
  tree->Branch("Tk",&Tk," Tk/D");
  tree->Branch("Tk_true",&Tk_true," Tk_true/D");
  tree->Branch("dTk",&dTk," dTk/D");
  tree->Branch("Sa",&Sa," Sa/D");
  tree->Branch("Sa_true",&Sa_true," Sa_true/D");
  tree->Branch("dSa",&dSa," dSa/D");
  tree->Branch("Iso",&Iso," Iso/D");
  tree->Branch("Iso_true",&Iso_true," Iso_true/D");
  tree->Branch("dIso",&dIso," dIso/D");
  tree->Branch("Hlt",&Hlt," Hlt/D");
  tree->Branch("Hlt_true",&Hlt_true," Hlt_true/D");
  tree->Branch("dHlt",&dHlt," dHlt/D");
  tree->Branch("chi2",&chi2," chi2/D");

  ifstream fin;
  fin.open("fitResult.txt");

  char line[1024];

  fin.getline(line, 1024);
  cout << line << endl;
  fin >> Y_true >> Tk_true >> Sa_true >> Iso_true >> Hlt_true;
  cout << "Yield = " << Y_true;
  cout << " eff_trk = " << Tk_true;
  cout << " eff_sa = " << Sa_true;
  cout << " eff_iso = " << Iso_true;
  cout << " eff_hlt = " << Hlt_true << endl;
  while(!(fin.eof())){
    Y = 0;
    fin >> Y >> dY >> Tk >>  dTk >>  Sa >>
      dSa >>  Iso >> dIso >>  Hlt >>  dHlt >>chi2;
    if(Y > 0)
      tree->Fill();
  }
  
  tree->Print();
  f->Write();
  f->Close();


}

int main(){
	
	float data[1024], x[1024], min, area, energy, timing;
	int start, cut, count;
	TFile *fout = TFile::Open("/projects/physics/dmice/Fermilab/Co60_Temp.root", "recreate");
		TTree *tree = new TTree("tree", "NaI Detector Data");
		tree->Branch("Waveforms", data, "data[1024]/F");
		tree->Branch("Height", &min, "min/F");
		tree->Branch("Area", &area, "area/F");
		tree->Branch("Saturated", &cut, "cut/I");
		tree->Branch("Energy", &energy, "energy/F");
		tree->Branch("Timing" , &timing, "timing/F");
	
	TChain t("tree");
		t.Add("/projects/physics/dmice/Fermilab/sourcedata/Co60raw/Waveforms_0.root");
		//t.Add("/projects/physics/dmice/Fermilab/1V_runs/Fermilab_1V_raw.root");
		t.SetBranchAddress("Waveforms", data);
	
	for(int e = 0; t.GetEntry(e) > 0; e++){  //&& e < 10
		Float_t weight=0, avg = 0, y[1024], offavg=0, bin_sum = 0;
		int minLoc, n0;
		min = 0;
		for(int i = 0; i < 200; i++)
			offavg += data[i];
		offavg = offavg/200;
		for(int i = 0; i < 1024; i++){
			y[i] = data[i] - offavg;
			x[i] = i;
			if(i < 1021 && cut == 0)
				if(data[i] == data[i+1] && data[i] == data[i+2] && data[i] < -.2)
					cut = 1;
			if(y[i] < min){
				min = y[i];
				minLoc = i;
			}

		}
		int zero_counter = 0, start, end;
		bool foundstart = false, foundend = false;
		for(int i = 2; i < 997; i++){
			float local_avg = y[i-2] + y[i-1] + y[i] + y[i+1] + y[i+2];
			if(local_avg < -.05 && !foundstart){
				start = i-2;
				foundstart = true;
			}
			if(local_avg >= 0 && !foundend && foundstart){
				zero_counter++;
				if(zero_counter == 15){
					end = i;
					foundend = true;
				}
			}
		
		}
		//cout <<"Pulse Start: " << start << endl;
		//cout <<"Pulse End: " << end << endl;
		

	//	for(inti i = start; i <= end; i++){
	//		bin_sum += y[i];
	//	}

		TGraph gr(1024, x, y);
		area = gr.Integral(start, end);
	//	area = -bin_sum;

		energy = (area*69.674)+1.8148;
		
		for(int i = minLoc; i > 0; i--){
			if(y[i] < min/2)
				n0 = i;
		}		
		for(int i = n0; i < (n0+100); i++)
		{
			weight += (i-n0)*y[i];
			avg += y[i];
			if(i == 1023)
			{
				weight = 0;
				avg = 0;
				break;
			}
		}
		if(avg != 0)
			timing = (weight/avg)*4;
		tree->Fill();
		//cout << "min: " << min << endl;
		//cout << "Area: "<< area << endl;
		//cout << "Saturated: " << cut << endl;
		//cout << "Meantime: " << timing << endl;
		//cout << endl;

		count++;
	}
	tree->Write();
	tree->Print();
	fout ->Close("R");
	delete fout;
}

//.........这里部分代码省略.........
   newtree->Branch("nJets",&nJets,"nJets/I");
   newtree->Branch("nJets30",&nJets30,"nJets30/I");
   newtree->Branch("nHiggs",&nHiggs,"nHiggs/I");
   newtree->Branch("partID[nPart]",partID,"partID[nPart]/I");
   newtree->Branch("Motherindex[nPart]",Motherindex,"Motherindex[nPart]/I");
   newtree->Branch("b_mothers[bPartons]",b_mothers,"b_mothers[bPartons]/I");
   newtree->Branch("Higgs_mother[nHiggs]",Higgs_mother,"Higgs_mother[nHiggs]/I");
   newtree->Branch("MCpt[nPart]",MCpt,"MCpt[nPart]/F");
   newtree->Branch("b_pt[bPartons]",b_pt,"b_pt[bPartons]/F");
   newtree->Branch("Higgs_pt[nHiggs]",Higgs_pt,"Higgs_pt[nHiggs]/F");
   newtree->Branch("MCeta[nPart]",MCeta,"MCeta[nPart]/F");
   newtree->Branch("b_eta[bPartons]",b_eta,"b_eta[bPartons]/F");
   newtree->Branch("Higgs_eta[nHiggs]",Higgs_eta,"Higgs_eta[nHiggs]/F");
   newtree->Branch("MCphi[nPart]",MCphi,"MCphi[nPart]/F");
   newtree->Branch("b_phi[bPartons]",b_phi,"b_phi[bPartons]/F");
   newtree->Branch("Higgs_phi[nHiggs]",Higgs_phi,"Higgs_phi[nHiggs]/F");
   newtree->Branch("jet_pt[nJets]",jet_pt,"jet_pt[nJets]/F");
   newtree->Branch("jet_eta[nJets]",jet_eta,"jet_eta[nJets]/F");
   newtree->Branch("jet_phi[nJets]",jet_phi,"jet_phi[nJets]/F");
   newtree->Branch("jet_bdisc[nJets]",jet_bdisc,"jet_bdisc[nJets]/F");
   newtree->Branch("Ht",&Ht,"Ht/F");

// Address old tree's branches
   oldtree->SetBranchAddress("nGenPart", &nGenPart);
   oldtree->SetBranchAddress("nPFJets", &nPFJets);
   oldtree->SetBranchAddress("MCpx[nGenPart]",MCpx);
   oldtree->SetBranchAddress("MCpy[nGenPart]",MCpy);
   oldtree->SetBranchAddress("MCpz[nGenPart]",MCpz);
   oldtree->SetBranchAddress("pdgID[nGenPart]",pdgID);
   oldtree->SetBranchAddress("MCmotherind[nGenPart]",MCmotherind);
   oldtree->SetBranchAddress("jet_PF_pt[nPFJets]",jet_PF_pt);
   oldtree->SetBranchAddress("jet_PF_eta[nPFJets]",jet_PF_eta);
   oldtree->SetBranchAddress("jet_PF_phi[nPFJets]",jet_PF_phi);
   oldtree->SetBranchAddress("jet_PF_px[nPFJets]",jet_PF_px);
   oldtree->SetBranchAddress("jet_PF_py[nPFJets]",jet_PF_py);
   oldtree->SetBranchAddress("bdiscCSV_PF[nPFJets]",bdiscCSV_PF);

// Loop over each event
   for(int i = 0; oldtree->GetEntry(i) > 0; ++i)
   {	oldtree->GetEntry(i);
	bPartons = 0;
	bPartons30 = 0;
	nHiggs = 0;
	Ht = 0;	
	nPart = nGenPart;
	nJets = nPFJets;
	nJets30 = 0;

//     	Loop over each particle	
	for(int j = 0; j < nGenPart; ++j)
	{	partID[j] = pdgID[j];					//copies particle ID array
		
		MCpt[j] = sqrt(pow(MCpx[j],2)+pow(MCpy[j],2));		//calculate pT
	  float p_mag = sqrt(pow(MCpx[j],2)+pow(MCpy[j],2)+pow(MCpz[j],2));	//magnitude of 3p vector
		MCeta[j] = .5*log((p_mag + MCpz[j])/(p_mag - MCpz[j]));			//calculate eta
		MCphi[j] = atan2(MCpy[j],MCpx[j]);			//calculate phi
		
		Motherindex[j] = MCmotherind[j];			//copies mother index array

		if(fabs(pdgID[j]) == 5)	
		{	b_pt[bPartons] = MCpt[j];		//pt, eta and phi of b partons
			b_eta[bPartons] = MCeta[j];
			b_phi[bPartons] = MCphi[j];
			b_mothers[bPartons] = pdgID[MCmotherind[j]];	//stores ID of mother, rather than index
			++bPartons;		//count b partons
			if((MCpt[j]>30)&&(fabs(MCeta[j])<2.4)) ++bPartons30;
		}
		if(fabs(pdgID[j]) == 25) 
		{	Higgs_pt[nHiggs] = MCpt[j];		//pt, eta and phi of Higgs
			Higgs_eta[nHiggs] = MCeta[j];
			Higgs_phi[nHiggs] = MCphi[j];
			Higgs_mother[nHiggs] = pdgID[MCmotherind[j]];	//stores ID of mother, rather than index
			++nHiggs;		//count Higgs
		}
	}

//	Loop over each jet
	for(int j = 0; j < nPFJets; ++j)
	{	jet_pt[j] = jet_PF_pt[j]; 	//Copy values from old tree
		jet_eta[j] = jet_PF_eta[j];
		jet_phi[j] = jet_PF_phi[j];
		jet_bdisc[j] = bdiscCSV_PF[j];	

		if((jet_PF_pt[j] > 30.0)&&(fabs(jet_PF_eta[j])<2.4))	
		{	Ht += jet_PF_pt[j];	//Ht is sum of pT of jets with pT > 30 GeV
			++nJets30;
		}
	}
	
//	Fill new tree
	newtree->Fill();
   }

// Write new tree to output file
   newtree->Print();
   newtree->Write();

// Write output file
   ofile.Write();
}

int
TrackParametrization( TString csvfile="fitslices_out.csv" )
{

  /* Read data from input file */
  TTree *tres = new TTree();
  tres->ReadFile( csvfile, "ptrue:etatrue:psig:psig_err:pmean:pmean_err:norm", ',' );

  /* Print read-in tree */
  tres->Print();

  /* colors array */
  unsigned colors[8] = {1,2,3,4,6,7,14,16};

  /* Create vector of theta values to include for visualization*/
  vector< double > etas_vis;
  etas_vis.push_back(-2.75);
  etas_vis.push_back(-2.25);
  etas_vis.push_back(-1.75);
  etas_vis.push_back(-0.25);
  etas_vis.push_back( 0.25);
  etas_vis.push_back( 1.75);
  etas_vis.push_back( 2.25);

//  etas_vis.push_back(-3.25);
//  etas_vis.push_back(-2.25);
//  etas_vis.push_back(-1.25);
//  etas_vis.push_back(-0.25);
//  etas_vis.push_back( 0.25);
//  etas_vis.push_back( 1.25);
//  etas_vis.push_back( 2.25);
//  etas_vis.push_back( 3.25);

  /* Create vector of theta values to include for fitting*/
  vector< double > etas_fit;
  for ( double eta = -4.45; eta < 4.5; eta += 0.1 )
    etas_fit.push_back( eta );

  /* Create fit function */
  TF1* f_momres = new TF1("f_momres", "sqrt( [0]*[0] + [1]*[1]*x*x )" );

  cout << "\nFit function: " << f_momres->GetTitle() << "\n" << endl;

  /* Create scratch canvas */
  TCanvas *cscratch = new TCanvas("cscratch");

  /* Create framehistogram */
  TH1F* hframe = new TH1F("hframe","",100,0,40);
  hframe->GetYaxis()->SetRangeUser(0,0.15);
  hframe->GetYaxis()->SetNdivisions(505);
  hframe->GetXaxis()->SetTitle("Momentum (GeV/c)");
  hframe->GetYaxis()->SetTitle("#sigma_{p}/p");

  /* create combined canvas plot */
  TCanvas *c1 = new TCanvas();
  hframe->Draw();

  /* Create legend */
  TLegend* leg_eta = new TLegend( 0.2, 0.6, 0.5, 0.9);
  leg_eta->SetNColumns(2);

  /* Create ofstream to write fit parameter results */
  ofstream ofsfit("track_momres_new.csv");
  ofsfit<<"eta,par1,par1err,par2,par2err"<<endl;

  /* Create resolution-vs-momentum plot with fits for each selected theta value */
  for ( int i = 0; i < etas_fit.size(); i++ )
    {
      /* Switch to scratch canvas */
      cscratch->cd();

      double eta = etas_fit.at(i);

      /* No tracking outside -4 < eta < 4 */
      if ( eta < -4 || eta > 4 )
	continue;

      cout << "\n***Eta = " << eta << endl;

      /* Define range of theta because float comparison with fixed value doesn't work
	 too well for cuts in ROOT trees */
      double eta_min = eta * 0.999;
      double eta_max = eta * 1.001;

      /* Cut for tree */
      TCut cutx( Form("ptrue > 1 && ( (etatrue > 0 && (etatrue > %f && etatrue < %f)) || (etatrue < 0 && (etatrue < %f && etatrue > %f)) )", eta_min, eta_max, eta_min, eta_max) );

      /* "Draw" tree on scratch canvas to fill V1...V4 arrays */
      tres->Draw("psig:ptrue:psig_err:0", cutx );

      /* Create TGraphErrors with selected data from tree */
      TGraphErrors *gres = new TGraphErrors( tres->GetEntries(cutx),
					     &(tres->GetV2())[0],
					     &(tres->GetV1())[0],
					     &(tres->GetV4())[0],
					     &(tres->GetV3())[0] );

      /* reset function parameters before fit */
      f_momres->SetParameter(0,0.1);
      f_momres->SetParameter(1,0.1);
//.........这里部分代码省略.........

int main(int argc, char ** argv)
{
  // select fed
  int month = atoi(argv[1]);
  string month_name [12] = {"Gen","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"};

  if (month < 0 || month > 12) {
   cout << "Give as argument a number between 1 and 12 corresponding to the chosen month" << endl;
   return 0;
  }

  if (month > 10) {
   cout << "No data corresponding to the chosen month..yet!" << endl;
   return 0;
  }


  //Get old tree 
  TChain *tx = new TChain("x");
  tx->Add(argv[2]);
 
  init_ttree(tx, &x);

  TTree *oldtree = (TTree*)tx;

  Long64_t nentries = oldtree->GetEntries();
  cout<< "Number of entries in the tree : " << nentries << endl;


  //Create a new file + a clone of old tree in new file
  char fname[1000];
  sprintf(fname,argv[3],(month_name[month - 1]).c_str());
  TFile *newfile = new TFile(fname,"recreate");
  TTree *newtree = oldtree->CloneTree(0);

  int lowerRun, upperRun = 0;

  if (month == 1) {
    lowerRun = 153941;
    upperRun = 156054;
  }

  if (month == 2) {
    lowerRun = 156225;
    upperRun = 159130;
  }

  if (month == 3) {
    lowerRun = 159248;
    upperRun = 161732;
  }

  if (month == 4) {
    lowerRun = 161846;
    upperRun = 163726;
  }

  if (month == 5) {
    lowerRun = 163762;
    upperRun = 166139;
  }

  if (month == 6) {
    lowerRun = 166243;
    upperRun = 168148;
  }

  if (month == 7) {
    lowerRun = 168266;
    upperRun = 172266;
  }

  if (month == 8) {
    lowerRun = 172319;
    upperRun = 174912;
  }

  if (month == 9) {
    lowerRun = 175118;
    upperRun = 177519;
  }

  if (month == 10) {
    lowerRun = 177624;
    upperRun = 179998;
  }


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

    if (i%10000000==0) cout << "Analyzing entry " << i << endl;
    oldtree->GetEntry(i);

    if ( x.run > lowerRun && x.run < upperRun ) newtree->Fill();
  }
  
  newtree->Print();
  newtree->AutoSave();
  
  delete newfile;
//.........这里部分代码省略.........

int makePlot_LoI_Momentum(
                          const TString infile = "data/pythia.ep.100.test.root"
                          )
{
  gStyle->SetOptStat(kFALSE);

  /*--------------Get Input File--------------*/

  TFile *f_pyth = new TFile( infile, "OPEN");
  TTree* T = (TTree*)f_pyth->Get("tree");
  T->Print();

  /*------------------------------------------*/
  /*--------------Define Cuts-----------------*/

  //TCut Q2_cut = "Q2 > 1";

  TCut electron_cut = "p.fKF == 11";
  TCut hadron_cut = "abs(p.fKF) > 100";
  TCut proton_cut = "p.fKF == 2212";
  TCut neutron_cut = "p.fKF == 2112";
  TCut Kaon_cut = "abs(p.fKF) == 321 || p.fKF == 311";
  TCut Pion_charged_cut = "abs(p.fKF) == 211";
  //TCut Pion_cut = "abs(p.fKF) == 211 || p.fKF == 111";
  TCut photon_cut = "p.fKF == 22";

  TCut mother_cut = "p.fParent == 0";
  TCut status_cut = "p.fKS < 10";

  TCut eta_cut_n3n2 = "TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2)) > -3 && TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2)) < -2";
  TCut eta_cut_n2n1 = "TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2)) > -2 && TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2)) < -1";
  TCut eta_cut_n1z0 = "TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2)) > -1 && TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2)) < -0";

  /*------------------------------------------*/
  /*-------Momentum vs. Pseudorapidity--------*/
  /*------------------------------------------*/

  /*------------------------------------------*/
  /*---------Electrons (LoI Fig. 2-1)---------*/
  /*------------------------------------------*/
  TH2F *h_peta_e = new TH2F("h_peta_e", "", 100, -4, 3, 100, 0, 50 ); //250x015 10M

  TCanvas *c_peta_e = new TCanvas( "c_peta_e" );
  c_peta_e->SetLogz();

  T->Draw("sqrt((p.fPx)**2 + (p.fPy)**2 + (p.fPz)**2 ):TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2))>>h_peta_e", electron_cut && mother_cut && "Q2 > 1", "colz");

  h_peta_e->GetXaxis()->SetTitle("Pseudorapidity #eta");
  h_peta_e->GetYaxis()->SetTitle("Electron Momentum p_{e-} [GeV]");

  //  c_p_eta_e->Print("Plots/Pythia_peta_e_10M_250x010.eps");
  //  c_p_eta_e->Print("Plots/Pythia_peta_e_10M_250x010.png");

  /*------------------------------------------*/
  /*----------Hadrons (LoI Fig. 2-4)----------*/
  /*------------------------------------------*/
  TH2F *h_p_eta_h = new TH2F("h_p_eta_h", "", 60,-6,6, 200,0,225); //250x015 10M

  TCanvas *c_p_eta_h = new TCanvas( "c_p_eta_h" );
  // T->Draw("sqrt((p.fPx)**2 + (p.fPy)**2 + (p.fPz)**2 ):TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2))>>h_p_eta_h",hadron_cut && "Q2 > 1 && 0.01 < y < 0.80 && W2 > 10", "colz");
  T->Draw("sqrt((p.fPx)**2 + (p.fPy)**2 + (p.fPz)**2 ):TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2))>>h_p_eta_h",hadron_cut && "Q2 > 1 && 0.01 < y < 0.80", "colz");

  h_p_eta_h->GetXaxis()->SetTitle("Pseudorapidity #eta");
  h_p_eta_h->GetYaxis()->SetTitle("Hadron Momentum p_{Hadron} [GeV]");
  c_p_eta_h->SetLogz();

  //  c_p_eta_h->Print("Plots/Pythia_peta_h_10M_250x010.eps");
  //  c_p_eta_h->Print("Plots/Pythia_peta_h_10M_250x010.png");


  /*------------Momentum Spectra (Figure 2.2)--------------*/

  /*-------- -3 < Eta < -2 ---------*/

  TH1F* hp_e_n3n2 = new TH1F("hp_e_n3n2", "dN/dp vs. p", 60, 0, 30);
  TH1F* hp_p_n3n2 = new TH1F("hp_p_n3n2", "dN/dp vs. p", 60, 0, 30);
  TH1F* hp_y_n3n2 = new TH1F("hp_y_n3n2", "dN/dp vs. p", 60, 0, 30);

  TCanvas *cp_e_n3n2 = new TCanvas("cp_e_n3n2");
  T->Draw("sqrt(p.fPx**2 + p.fPy**2 + p.fPz**2)>>hp_e_n3n2", electron_cut && eta_cut_n3n2 && "Q2 > 0.01", "goff");
  T->Draw("sqrt(p.fPx**2 + p.fPy**2 + p.fPz**2)>>hp_p_n3n2", Pion_charged_cut && eta_cut_n3n2 && "Q2 > 0.01", "goff");
  T->Draw("sqrt(p.fPx**2 + p.fPy**2 + p.fPz**2)>>hp_y_n3n2", photon_cut && eta_cut_n3n2 && status_cut && "Q2 > 0.01", "goff");


  TH1F* htmp_n3n2 = hp_e_n3n2->Clone();
  htmp_n3n2->SetTitle("");
  htmp_n3n2->GetXaxis()->SetTitle("p [GeV]");
  htmp_n3n2->GetYaxis()->SetTitle("dN/dp");
  htmp_n3n2->SetMaximum( 0.99e7);
  htmp_n3n2->Draw();

  hp_e_n3n2->SetLineColor(2);
  hp_e_n3n2->Draw("same");
  hp_p_n3n2->SetLineColor(1);
  hp_p_n3n2->Draw("same");
  hp_y_n3n2->SetLineColor(4);
  hp_y_n3n2->Draw("same");
  cp_e_n3n2->SetLogy();

  TLegend* leg_n3n2 = new TLegend(0.53,0.67,0.73,0.90);
//.........这里部分代码省略.........