C++ TNtuple::GetEntries方法代码示例

void PixelMergeSmallFiles(){

    TFile * oFilebDist = new TFile("bDistr2.root", "RECREATE");
    TNtuple* SimEventsGlobal = new TNtuple("SimEventGlob", "SimEventGlob", "EvN:b");
    
    TFile * oFileLinks = new TFile("Links2.root", "RECREATE");
    TNtuple* LinksGlobal = new TNtuple("LinksGlob", "LinksGlob", "EvN:fedid:linkn:nHits");
    
    char FileInNumber[5];
    string FileInPath = "/net/pstore01/d00/scratch/icali/CMSSW_2_1_11/PixelAnalysis/PixelNTuple_hydjet_x2_mb_oldPL_d20081106/";
    string FileInNameRoot= "hydjet_x2_mb_oldPL_d20081106_r0";
    
    Float_t j = 0;    
    for(int FileN = 901; FileN <1801; ++FileN){
    
         sprintf(FileInNumber, "%.5d", FileN);
         string FileInName = FileInPath+FileInNameRoot+FileInNumber+".root";
         TFile *iFile = new TFile((const char*)FileInName.c_str());
         if(!iFile->IsZombie()){
           PixelAnalyzer->cd();
    
           TNtuple *FEDLinks = (TNtuple*)iFile->FindObjectAny("Links");
           TNtuple *SimEvents = (TNtuple*)iFile->FindObjectAny("SimEvent");
    
           Float_t FEDEvN, fedid, linkn, nHits;
           FEDLinks->SetBranchAddress("EventN",&FEDEvN);
           FEDLinks->SetBranchAddress("fedid",&fedid);
           FEDLinks->SetBranchAddress("linkn",&linkn);
           FEDLinks->SetBranchAddress("nHits",&nHits);
           Long64_t FEDLinkLenght =FEDLinks->GetEntries(); 
    
           Float_t SimEvN, b;
           SimEvents->SetBranchAddress("EventN",&SimEvN);
           SimEvents->SetBranchAddress("mult",&b);
           Long64_t SimEventLenght =SimEvents->GetEntries(); 
           
           Long64_t i;
           for(i=0; i < SimEventLenght; ++i){
             SimEvents->GetEntry(i);
             SimEventsGlobal->Fill(i, b);
           } 
           
          Float_t OldEvN= -1;        
          for(i=0; i < FEDLinkLenght; ++i){
             FEDLinks->GetEntry(i);
             if(OldEvN != FEDEvN){
                OldEvN= FEDEvN;
                ++j;
             }
             LinksGlobal->Fill(j, fedid, linkn, nHits);
             
          }
      }
    }
    oFilebDist->Write();
    oFileLinks->Write();
}

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

  gROOT->SetStyle ("Plain") ;

  // Open a stream connected to an event file:
  if (argc < 3) exit (1) ;

  TNtuple ggFntuple ("ggFntuple", "ggFntuple", "Minv:Deta:DR:Dphi") ;
  fillHistos (argv[1], ggFntuple) ;
  TNtuple vbFntuple ("vbFntuple", "vbFntuple", "Minv:Deta:DR:Dphi") ;
  fillHistos (argv[2], vbFntuple) ;

  int bins = 50 ;
  TH2F Minv_vs_deta_cont ("Minv_vs_deta_cont", "contamination in VBF selections", 
                          bins, 0, 1500, bins, 0, 5) ;
  Minv_vs_deta_cont.SetStats (0) ;                          
  Minv_vs_deta_cont.GetXaxis ()->SetTitle ("M_{inv}") ;                          
  Minv_vs_deta_cont.GetYaxis ()->SetTitle ("#Delta#eta") ;                          

  double num_ggF = ggFntuple.GetEntries () ;
  double num_vbF = vbFntuple.GetEntries () ;
  double ggF_Xsec = 14.761 ; // pb NNLO
  double vbF_Xsec = 1.6864 ; // pb NLO
 
  for (double iMinv = 0.5 * (1500./bins) ; iMinv < 1500. ; iMinv += 1500./bins)
    for (double iDeta = 0.5 * (5./bins) ; iDeta < 5. ; iDeta += 5./bins)
      {
        char cut[20] ;
        sprintf (cut, "Minv>%f&&Deta>%f", iMinv, iDeta) ;
//        std::cout << cut << std::endl ;
        double num_ggF_sel = ggFntuple.GetEntries (cut) ;
        double num_vbF_sel = vbFntuple.GetEntries (cut) ;
        double val = (ggF_Xsec * num_ggF_sel / num_ggF) /
          (ggF_Xsec * num_ggF_sel / num_ggF + vbF_Xsec * num_vbF_sel / num_vbF) ;
        Minv_vs_deta_cont.Fill (iMinv, iDeta, val) ;
      }

  std::cout << "out of loop" << std::endl ;

  TCanvas c1 ;
  Minv_vs_deta_cont.Draw ("COLZ") ;
  c1.Print ("Minv_vs_deta_cont.gif","gif") ;

  TFile output ("ouptut.root","recreate") ;
  output.cd () ;
  Minv_vs_deta_cont.Write () ;
  output.Close () ;
    
  return 0 ;
}

void read_ntuple_from_file(){

  int i,j,k,n;

  TFile *in = new TFile("ntupleoutputsample.root");
  TNtuple *data = (TNtuple*) in->GetObjectChecked("data","TNtuple");

  double pot,cur,temp,pres;
  float *row_content; //Must necessarily be float and not a double... WHY?

  TH1D *histo = new TH1D("histo","HISTO",100,0,10);

  cout << "Potential\tCurrent\tTemperature\tPressure" << endl;
  for(i=0;i<data->GetEntries();++i){
      data->GetEntry(i);
      row_content = data->GetArgs();
      pot = row_content[0];
      cur = row_content[1];
      temp = row_content[2];
      pres = row_content[3];
      cout << pot << "\t" << cur << "\t" << temp << "\t" << pres << endl;
      histo->Fill(pot);
  }

  histo->Draw();
}

void ROCOccupancy(Int_t StartRun, Int_t EndRun){

    using namespace std;
       
   
   //string FileRoot = "/home/l_tester/log/bt05r";
   string FileInName = "ProcessedData.root";
   char RunNumber[6];
    
    
    Float_t j = 0;    
    for(int RunN = StartRun; RunN <=EndRun; ++RunN){
      string FileRoot = "/d00/icali/PixelHwStudies/PSIStuff/log/bt05r";        
      sprintf(RunNumber, "%.6d", RunN);
      FileRoot = FileRoot + RunNumber + "/";
      FileInName = FileRoot + FileInName;
        
      TFile *iFile = new TFile((const char*)FileInName.c_str());
      if(!iFile->IsZombie()){
        TNtuple *Digis = (TNtuple*)iFile->FindObjectAny("events");
        Int_t EvN, roc, ROCHits, MeanPh, NCols;
        Digis->SetBranchAddress("ROCHits", &ROCHits);
        Digis->SetBranchAddress("MeanPh", &MeanPh);
        Digis->SetBranchAddress("roc", &roc);
        Digis->SetBranchAddress("NCols", &NCols);
        Digis->SetBranchAddress("EvN", &EvN);
        Int_t iNTlenght = Digis->GetEntries();
        
        Int_t MaxROCHits =0;
        for(int i =0 ; i < iNTlenght; ++i){
          Digis->GetEntry(i);
          if(ROCHits > MaxROCHits) MaxROCHits = ROCHits;
        }
        cout << "Run: " << RunN << " MaxROCHits: " << MaxROCHits << endl;
      }
   }
}

void fitBjetJES(int ppPbPb=1, int cbinlo=12, int cbinhi=40){

  if(!ppPbPb){
    cbinlo=0;
    cbinhi=40;
  }

  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);

  TFile *fL;
  
  if(!ppPbPb)fL=new TFile("histos/ppMC_hiReco_jetTrig_highPurity_JEC.root");
  else     fL=new TFile("histos/PbPbQCDMC_pt30by3_ipHICalibCentWeight.root");

  // these are dummy files for pp
  TFile *fB=new TFile("histos/PbPbBMC_pt30by3_ipHICalibCentWeight.root");
  TFile *fC=new TFile("histos/PbPbCMC_pt30by3_ipHICalibCentWeight.root");


  TNtuple *tL = (TNtuple*) fL->Get("nt");
  TNtuple *tB = (TNtuple*) fB->Get("nt");
  TNtuple *tC = (TNtuple*) fC->Get("nt");
  
  float jtptL, refptL, jtetaL, weightL, refparton_flavorForBL, binL;

  tL->SetBranchAddress("jtpt",&jtptL);
  tL->SetBranchAddress("jteta",&jtetaL);
  tL->SetBranchAddress("refpt",&refptL);
  tL->SetBranchAddress("weight",&weightL);
  if(ppPbPb)tL->SetBranchAddress("bin",&binL);
  tL->SetBranchAddress("refparton_flavorForB",&refparton_flavorForBL);

  float jtptB, refptB, jtetaB, weightB, refparton_flavorForBB, binB;

  tB->SetBranchAddress("jtpt",&jtptB);
  tB->SetBranchAddress("jteta",&jtetaB);
  tB->SetBranchAddress("refpt",&refptB);
  tB->SetBranchAddress("weight",&weightB);
  if(ppPbPb)tB->SetBranchAddress("bin",&binB);
  tB->SetBranchAddress("refparton_flavorForB",&refparton_flavorForBB);

  float jtptC, refptC, jtetaC, weightC, refparton_flavorForBC, binC;

  tC->SetBranchAddress("jtpt",&jtptC);
  tC->SetBranchAddress("jteta",&jtetaC);
  tC->SetBranchAddress("refpt",&refptC);
  tC->SetBranchAddress("weight",&weightC);
  if(ppPbPb)tC->SetBranchAddress("bin",&binC);
  tC->SetBranchAddress("refparton_flavorForB",&refparton_flavorForBC);

  TProfile  *hL = new TProfile("hL","hL",250,50,300,0,10);
  TProfile  *hB = new TProfile("hB","hB",250,50,300,0,10);
  TProfile  *hC = new TProfile("hC","hC",250,50,300,0,10);
  hL->Sumw2(),hB->Sumw2(),hC->Sumw2();


  for(int i=0;i<tL->GetEntries();i++){
    tL->GetEntry(i);
    if(!ppPbPb) binL=39;

    if(fabs(jtetaL)<2 && binL>=cbinlo && binL<cbinhi)
      hL->Fill(refptL,jtptL/refptL,weightL); 


    if(!ppPbPb){
      if(fabs(jtetaL)<2 && binL>=cbinlo && binL<cbinhi && abs(refparton_flavorForBL)==5)
	hB->Fill(refptL,jtptL/refptL,weightL);
      
      if(fabs(jtetaL)<2 && binL>=cbinlo && binL<cbinhi && abs(refparton_flavorForBL)==4)
	hC->Fill(refptL,jtptL/refptL,weightL);      
    }
  }

  if(ppPbPb){
    for(int i=0;i<tB->GetEntries();i++){
      tB->GetEntry(i);
      if(fabs(jtetaB)<2 && binB>=cbinlo && binB<cbinhi && abs(refparton_flavorForBB)==5)
	hB->Fill(refptB,jtptB/refptB,weightB);
    }
    for(int i=0;i<tC->GetEntries();i++){
      tC->GetEntry(i);
      if(fabs(jtetaC)<2 && binC>=cbinlo && binC<cbinhi && abs(refparton_flavorForBC)==4)
	hC->Fill(refptC,jtptC/refptC,weightC);
    }
  }

  
 
  hL->SetMinimum(0.);
  
  hL->SetLineColor(kBlue);
  hB->SetLineColor(kRed);
  hC->SetLineColor(kGreen);

  hL->SetMarkerColor(kBlue);
  hB->SetMarkerColor(kRed);
  hC->SetMarkerColor(kGreen);
  
  //hL->SetMarkerStyle(4);
//.........这里部分代码省略.........

int main(int argc, char **argv)
{
	TFile *f;
	TFile *newf;
	TNtuple *ntuple;
	TNtuple *newntuple;
	TH1F *h1;
	TF1 *func;
	TString Metal;
	TString dataLoc;
	Int_t dataSL;
	
	if(argc < 5 || argc > 6){
		std::cout << "The number of arguments is incorrect" << std::endl;
		return 0;
	}
	
	dataLoc = argv[1];
	PHI_MIN = (Double_t) std::stod(argv[2]);
	PHI_MAX = (Double_t) std::stod(argv[3]);
	N_PHI = (Int_t) std::stoi(argv[4]);
	if(argc >= 6) Metal = argv[5];
		
	dataSL = dataLoc.Length();
	if(dataLoc(dataSL-1, 1) != "/"){
		dataLoc = dataLoc + "/";
	}	
	
	std::cout << "The following settings are being used" << std::endl;
	std::cout << "Data Directory = " << dataLoc << std::endl;
	std::cout << PHI_MIN << " < PHI < " << PHI_MAX << ", N_PHI = " << N_PHI << std::endl;
	if(argc == 6) std::cout << "Running only for Metal "<< Metal << std::endl;
	else std::cout << "Running all Metals" << std::endl;	
	
	Float_t Q2, Xb, Zh, Pt, Phi, Val, Err;
	Float_t A, Ac, Acc;
	Float_t AErr, AcErr, AccErr;
	Float_t ChiSQ;

	Int_t nentries, empty;

	if(Metal == "") N_METAL = 6;
	else N_METAL = 1;
	
	for(Int_t met = 0; met < N_METAL; met++){
		if(Metal == ""){
			if(met == 0) Metal = "C";
			else if(met == 1) Metal = "Fe";
			else if(met == 2) Metal = "Pb";
			else if(met == 3) Metal = "D_C";
			else if(met == 4) Metal = "D_Fe";
			else if(met == 5) Metal = "D_Pb";
		}
		if(dataLoc == "")
			f = new TFile("../Gen5DimData/" + Metal + "_5_dim_dist.root", "READ");
		else
			f = new TFile(dataLoc + Metal + "_5_dim_dist.root", "READ");
		ntuple = (TNtuple*) f->Get("fit_data");

		nentries = ntuple->GetEntries();
		ntuple->SetBranchAddress("Xb", &Xb);
		ntuple->SetBranchAddress("Q2", &Q2);
		ntuple->SetBranchAddress("Xb", &Xb);
		ntuple->SetBranchAddress("Zh", &Zh);
		ntuple->SetBranchAddress("Pt", &Pt);
		ntuple->SetBranchAddress("Phi", &Phi);
		ntuple->SetBranchAddress("Val", &Val);
		ntuple->SetBranchAddress("Err", &Err);

		newntuple = new TNtuple("AAcAcc_data", "AAcAcc_data", "Q2:Xb:Zh:Pt:A:AErr:Ac:AcErr:Acc:AccErr:ChiSQ");
		func = new TF1("fit", "[0]+TMath::Cos(x*TMath::Pi()/180)*[1]+TMath::Cos(2*x*TMath::Pi()/180)*[2]");

		newf = new TFile(Metal + "newphihist.root", "RECREATE");
		newf->cd();

		for(Int_t i = 0; i < nentries; i = i + N_PHI){
			ntuple->GetEntry(i);
			h1 = new TH1F((const char*) Form("PhiDist Q2=%.3f Xb=%.3f Zh=%.3f Pt=%.3f", Q2, Xb, Zh, Pt), 
							(const char*) Form("PhiDist Q2=%.3f Xb=%.3f Zh=%.3f Pt=%.3f", Q2, Xb, Zh, Pt), N_PHI, PHI_MIN, PHI_MAX);
			empty = 0;
			for(Int_t j = 1; j <= N_PHI; j++){
				ntuple->GetEntry(i+j-1);
				h1->SetBinContent(j, Val);
				if(h1->GetBinContent(j) == 0)
					empty++;
				h1->SetBinError(j, Err*1.04);
			}
			if(empty <= (N_PHI - MIN_BIN)){
				h1->Fit(func, "q");
				h1->Write();
				if(func->GetNDF() != 0){
					ChiSQ = func->GetChisquare();
					A = func->GetParameter(0);
					AErr = func->GetParError(0);
					Ac = func->GetParameter(1);
					AcErr = func->GetParError(1);
					Acc = func->GetParameter(2);
					AccErr = func->GetParError(2);
					newntuple->Fill(Q2, Xb, Zh, Pt, A, AErr, Ac, AcErr, Acc, AccErr, ChiSQ);
				}
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	TString Fit_Cut = "(abs(" + Fit_Status + "-3.0)<"+Double_Tolerance+")";


	//	param1_gen == notgen	&&	param1_err == 0
	TString Param_1_Cut = "(abs(" + param1_gen + "-" + notgen +")<" + Double_Tolerance + ")&&(abs("+ param1_err + "-0.0)<"+ Double_Tolerance + ")";
	//	param2_gen == notgen	&&	param2_err == 0
	TString Param_2_Cut = "(abs(" + param2_gen + "-" + notgen + ")<"+ Double_Tolerance + ")&&(abs(" +param2_err +"-0.0)<" + Double_Tolerance + ")";

	//	Combine the individual Cuts
	TString Fit_Cut_String = Param_1_Cut + "&&" + Param_2_Cut + "&&" + Fit_Cut;



	//	Toys have a defined generation Value, Fits to Data DO NOT

	//	param1_gen != notgen
	TString p1isatoy = "(abs(" + param1_gen + "-" + notgen + ")>" + Double_Tolerance + ")";
	//	param2_gen != notgen
	TString p2isatoy = "(abs(" + param2_gen + "-" + notgen + ")>" + Double_Tolerance + ")";

	//	Combine the individual Cuts
	TString Toy_Cut_String = p1isatoy + "&&" + p2isatoy;




	//	Check for Toys in the file and wether I should run the FC code

	allresults->Draw( NLL, Toy_Cut_String, "goff" );
	bool Has_Toys = allresults->GetSelectedRows() > 0;

	cout << endl << "NUMBER OF TOYS IN FILE:\t" << allresults->GetSelectedRows() << endl;

	if( int(allresults->GetEntries() - allresults->GetSelectedRows()) == 0 )
	{
		cerr << "SERIOUS ERROR:\tSOMETHING HAS REALLY GOTTEN SCREWED UP!" << endl;
		exit(-3498);
	}


	//	Fit values for the global fit are now stored in:
	//
	//	Global_Best_NLL,	best_fit_values,	x_point,	y_point

	//	Tell the user
	cout << "GLOBAL DATA BEST FIT NLL:\t" << setprecision(10) << Global_Best_NLL << "\tAT:\tX:" << setprecision(10) << x_point << "\tY:\t" <<setprecision(10)<< y_point << endl;

	//	Check wether the minima as defined from the Global fit is the true minima within the phase-space

	//Check_Minima( allresults, Fit_Cut_String, &Global_Best_NLL, NLL, Double_Tolerance, param1_val, param2_val );


	TString NLL_Min;		//	Of course ROOT doesn't have USEFUL constructors!
	NLL_Min+=Global_Best_NLL;

	//	Plot the distribution of successfully fitted grid points for the PLL scan
	//	NB:	For FC this will likely saturate due to multiple layers of fits
	cout << endl << "FOUND UNIQUE GRID POINTS, PLOTTING" << endl;

	LL2D_Grid( allresults, Fit_Cut_String, param1string, param2string, rand_gen, "LL", outputdir );

	//	Construct a plot string for the NLL plot and plot it

	//	PLL part of Draw_String
	TString PLL = "(" + NLL + "-" + NLL_Min + ")";
	//	Gridding part of Draw_String

//iteration code
void iterate(TrkSettings s,int iter, int stepType, bool doCondor, bool testErrors = false)
{
  float pt, eta, phi, weight, centPU, rmin, maxJetPt,trkStatus,pNRec,mpt,mtrkQual,nEv; 
  
  TFile * histFile;
  std::string ifPP = "";
  if(s.nPb==0) ifPP = "pp_";
  if(iter==0) histFile = TFile::Open(Form("%scorrHists_job%d.root",ifPP.c_str(),s.job),"recreate");
  else        histFile = TFile::Open(Form("%scorrHists_job%d.root",ifPP.c_str(),s.job),"update");

  //make needed gen skim if it hasn't been done yet
  if(iter==0)
  {
    TH1D *genPre[20], *mrecoPre[20];
    TH2D *genPre2[20], *mrecoPre2[20];
    std::cout << "Denominator info not yet calculated; calculating and saving it..." << std::endl;
    for(int i = 0; i<8; i++)
    {
      if(i==6) continue;
      if(i != 1 && i!=7)
      {
        genPre[i] = makeTH1(s,i,"gen");
        mrecoPre[i] = makeTH1(s,i,"mreco");
      }
      else
      {
        genPre2[i] = makeTH2(s,i,"gen");
        mrecoPre2[i]= makeTH2(s,i,"mreco");
      }
    }
   
    TFile * skim;
    if(doCondor)
    {
      if(s.reuseSkim) skim = TFile::Open(Form("/mnt/hadoop/cms/store/user/abaty/tracking_Efficiencies/ntuples/%strackSkim_job%d.root",ifPP.c_str(),s.job),"read");
      else skim = TFile::Open(Form("%strackSkim_job%d.root",ifPP.c_str(),s.job),"read");
    }
    else skim = TFile::Open(Form("/export/d00/scratch/abaty/trackingEff/ntuples/%strackSkim_job%d.root",ifPP.c_str(),s.job),"read");
    //for efficiency
    std::cout << "Doing Efficiency denominator" << std::endl;   
    TNtuple * gen = (TNtuple*)  skim->Get("Gen");
    gen->SetBranchAddress("genPt",&pt);
    gen->SetBranchAddress("genEta",&eta); 
    gen->SetBranchAddress("genPhi",&phi);
    gen->SetBranchAddress("weight",&weight);
    gen->SetBranchAddress("centPU",&centPU);
    gen->SetBranchAddress("rmin",&rmin);
    gen->SetBranchAddress("jtpt",&maxJetPt);
    gen->SetBranchAddress("pNRec",&pNRec); 
    gen->SetBranchAddress("nEv",&nEv); 
	   
    for(int i = 0; i<gen->GetEntries(); i++)
    {
      gen->GetEntry(i);
      if(s.doSplit && nEv==1) continue;
      genPre[0]->Fill(pt,weight);
      genPre2[1]->Fill(eta,phi,weight); 
      genPre[2]->Fill(centPU,weight);
      genPre[3]->Fill(maxJetPt,weight);
      genPre[4]->Fill(eta,weight); 
      genPre[5]->Fill(rmin,weight);
      genPre2[7]->Fill(eta,pt,weight);
    }

    //for fake
    std::cout << "Doing Fake denominator" << std::endl;   
    TNtuple * reco = (TNtuple*)  skim->Get("Reco"); 
    reco->SetBranchAddress("trkPt",&pt);
    reco->SetBranchAddress("trkEta",&eta);
    reco->SetBranchAddress("trkPhi",&phi);
    reco->SetBranchAddress("weight",&weight);
    reco->SetBranchAddress("centPU",&centPU);
    reco->SetBranchAddress("rmin",&rmin);
    reco->SetBranchAddress("jtpt",&maxJetPt);
    reco->SetBranchAddress("trkStatus",&trkStatus);
    reco->SetBranchAddress("nEv",&nEv); 
    for(int i = 0; i<reco->GetEntries(); i++)
    {
      reco->GetEntry(i);
      if(s.doSplit && nEv==1) continue;
      if(trkStatus<-100) continue;
      mrecoPre[0]->Fill(pt,weight);
      mrecoPre2[1]->Fill(eta,phi,weight); 
      mrecoPre[2]->Fill(centPU,weight);
      mrecoPre[3]->Fill(maxJetPt,weight);
      mrecoPre[4]->Fill(eta,weight); 
      mrecoPre[5]->Fill(rmin,weight);
      mrecoPre2[7]->Fill(eta,pt,weight);
    }
  
    //Secondary calculation (no iterations)
    std::cout << "Quickly calculating the Secondary Rate from the reco tree (No further iterations needed)" << std::endl;
    TH2D * Secondary_Matched = new TH2D("Secondary_Matched",";pt;",s.multiRecoBins.at(s.job/s.nPtBinCoarse),s.ptMin,s.ptMax,24,-2.4,2.4); 
    TH2D * Secondary_Secondaries = new TH2D("Secondary_Secondaries",";pt;",s.multiRecoBins.at(s.job/s.nPtBinCoarse),s.ptMin,s.ptMax,24,-2.4,2.4); 
    for(int i = 0; i<reco->GetEntries(); i++)
    {
      reco->GetEntry(i);
      if(s.doSplit && nEv==1) continue;
      if(trkStatus>-100)
//.........这里部分代码省略.........

float n3HeSim(int nev_=0, float rfsf_phi_=300, int sev=0)
{
  // commandline option '-var=val' to change any global
  for (int i=1; i<gApplication->Argc(); ++i) 
  {
      TString opt(gApplication->Argv(i));
      if ((opt[0]=='-' || opt[0]=='+') && opt.Contains("="))
	  gROOT->ProcessLine(opt=opt.Append(";").Strip(TString::kLeading,'-'));
  }

  if (nev_) 
      nev = nev_;
  if (rfsf_phi_!=0) 
      rfsf_phi=rfsf_phi_/360.;
  TFile* f = new TFile("/home/siplu/GIT/n3He_Soft/Github/Simulation/data/mcstas.root");
  TNtuple* ntp = (TNtuple*)f->Get("mcstas");

  // set up ntuple
  BRANCH(l);
  BRANCH(pol);
  BRANCH(p5);
  BRANCH(p6);
  BRANCH(t6);
  BRANCH(x6);
  BRANCH(y6); 
  BRANCH(vx6);
  BRANCH(vy6);
  BRANCH(vz6);
  
  float chop[4]; 
  TBranch* b_chop[4];
  for (int j=0;j<nch;++j) 
  {
      b_chop[j]=ntp->GetBranch(vch[j]);
      ntp->SetBranchAddress(vch[j],&chop[j]);
  }

  // precalculate track projections
  float dwc = lwc/nz;			// WC parameters
  float da  = 2./na, a0=-1+da/2;	// polar [cos(theta)] integration
  float db  = 2*PI/nb, b0=-PI+db/2;	// azimuthal [phi] integration
  float cxyz[3], &cx=cxyz[0],&cy=cxyz[1],&cz=cxyz[2];	// projection cosines
  float ccz[NDIV], ccx[NDIV][NDIV], ccy[NDIV][NDIV], as; // cached values
  float aa=a0; 
  for (int ia=0; ia<na; ++ia, aa+=da) 
  {
      ccz[ia]=aa; 
      if (ccz[ia]==0) 
	  ccz[ia]=1e-11; 
      as=sqrt(1-aa*aa);
      float bb=b0; 
      for (int ib=0; ib<nb; ++ib, bb+=db) 
      {
	  ccx[ia][ib]=as*cos(bb); 
	  if (ccx[ia][ib]==0) 
	      ccx[ia][ib]=1e-11;
	  ccy[ia][ib]=as*sin(bb); 
	  if (ccy[ia][ib]==0) 
	      ccy[ia][ib]=1e-11;
      }
  }

  // cumulative beta distribution
  for (int i=0;i<nbet;++i) 
      beta_int[i+1]=beta_int[i]+beta_pt[i];

  // initialize WC response matrices
  double  n5 = 0;			// pre-polarizer neutrons
  double   N[NTOF]={0};			// chopped neutrons at end of SMpol
  double P2N[NTOF]={0};			//  " * P^2
  double  WN[NTOF]={0};			// captures in 3He target
  double  BN[NTOF][nwc]={{0}};		// N_i = N beta_i    = \sum{n_i} ions
  double GBN[NTOF][nwc]={{0}};		// M_i = N_i gamma_i = \sum{n_i gamma}
  double D2N[NTOF][nwc][nwc]={{{0}}};	// delta^2 N_ij = covariance
  double  GW[NTOF][nwc]={{0}};		// weight of each element

  // loop over events, fill "n","p2n" histograms, wc response matrices
  if (ntp->GetEntries()<nev) 
      nev=ntp->GetEntries();
  for (int iev=sev;iev<sev+nev;++iev) 
  {
      if (verbose && cnt) 
      {		// :........:....
	  if      (iev%(50*cnt)==0) 
	      cout<<endl<<" :"<<flush;
	  else if (iev%(10*cnt)==0) 
	      cout<<":"<<flush;
	  else if (iev%cnt==0) cout<<"."<<flush;
      }

      // read variables
      ntp->GetEntry(iev);
      hl5.Fill(l,p5*1.4/2*60); 
      n5 += p5;
      hl6.Fill(l,p6*1.4/2*60);

      // polarizer
      float tx = 1; //Transmission
      b_pol->GetEntry(iev);
      if (pol!=pol) 
//.........这里部分代码省略.........

//.........这里部分代码省略.........
for(int k = 1; k <= 6; k++){
		TString Znum = Form("%i", k);
	hist1->SetTitle("Z=" + Znum);


	//ALL UNITS ARE cm!
	Double_t detectorSideLength = 4; //40mm, as e.haettner H1 detector
	Double_t scatteringDistance = detectorDistance - phantomCenterDistance; //temporarily hard-coded, should be distance from target-center to detector
	Double_t degrees; //< actually radians
	
	Double_t r, rMin, rMax, deltaOmega, normFloat;
	TString rMinString, rMaxString, normString;
	TString same = "";
	TString histName;
	TCanvas *c3 = new TCanvas("histograms", "Distribution (at different angles)");
	int i = 0; //so that the degree steps can be varied to unevenly spaced values separate counter is used

	std::cout << "The following numbers also make it possible to make number of fragments comparison to the graph in A1 of E.Haettner\n";
	for(Double_t j = 0.0; j <= 8.0; j=j+1.0){
		i++;
		degrees = j * TMath::DegToRad();
		//std::cout << "plotting for Z = " << Znum << " at " << j << " degrees\n";
		//Distance from straight beam at the requested angle
		r = scatteringDistance * TMath::Tan(degrees);
		//now the "detector is rotated around all possible perpendicularlynangle values to beamline".
		//This forms an annulus with rMin and RMax as otuer and inner radiuses
		//Notice this will give a bit of approximation at small angles where at 0 degrees this gives a round sensor.
		Double_t deltaPhi = TMath::ATan((TMath::Cos(degrees)*detectorSideLength)/(2*scatteringDistance));
		rMin = TMath::Max(0.0,r - (detectorSideLength/(2*TMath::Cos(degrees))));
		rMax = rMin + ((detectorSideLength*TMath::Sin(degrees))/TMath::Tan((TMath::Pi()/2) - degrees - deltaPhi)) + (detectorSideLength*TMath::Cos(degrees));
		rMinString = Form("%f", rMin);
		rMaxString = Form("%f", rMax);
		//normalization of the bins.
		deltaPhi = degrees - TMath::ATan(TMath::Tan(degrees) - detectorSideLength/(2*scatteringDistance)); // this should be around arctan(detectorsidelength/sd)
		if(j != 0.0){
		deltaOmega = 2*TMath::Pi()*(TMath::Cos(TMath::Max(0.0,degrees-deltaPhi)) - TMath::Cos(degrees+deltaPhi));
		}else{
		deltaOmega = 4 * TMath::ASin(pow(detectorSideLength,2.0) / (4*pow(scatteringDistance,2) + pow(detectorSideLength,2)) );
		}
		normFloat = deltaOmega * events * binWidth;
		normString = Form("/%f", normFloat);

	// The following is veryvery ugly relies on a bunch of hardcoded histograms	because other solutions did not work
		histName = Form("hist%i", i);
		if(j != 0.0){
		fragments->Project(histName,"energy", "(Z == " + Znum + " && energy > 0 && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")" + normString);
		}else{
		fragments->Project(histName,"energy", "(Z == " + Znum + " && energy > 0 && posZ < " + rMaxString + "&& posY < " + rMaxString + " && posY > 0 && posZ > 0)" + normString);
		}
		int numEntries = fragments->GetEntries("(Z == " + Znum + " && energy > 0 && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")");
		std::cout << "\nj: "<< numEntries/(deltaOmega * events) << " entries for " << j;
		
		}

//the ugly hardcoded histograms being plotted
		//0 degrees
		hist1->SetLineColor(kBlue);
		hist1->Draw();
		//1 degree
		hist2->SetLineColor(kGreen);
		hist2->Draw("same"); //add "same" when also plotting 0 degrees
		//2 degrees
		hist3->SetLineColor(kRed);
		hist3->Draw("same");
		//3 degrees
		//hist4->SetLineColor(kGreen + 5);
		//hist4->Draw("same");
		//4 degrees
		hist5->SetLineColor(kGreen + 3); //gives a darker shade of green
		hist5->Draw("same");
		//5 degrees
		//hist6->SetLineColor(kRed);
		//hist6->Draw("same");
		//6 degrees
		hist7->SetLineColor(kRed);
		hist7->Draw("same");
		//7 degrees
		//hist8->SetLineColor(kRed);
		//hist8->Draw("same");
		//8 degrees
		//hist9->SetLineColor(kRed);
		//hist9->Draw("same");

   // Legends for the data
   leg = new TLegend(0.9,0.7,1,1);  //coordinates are fractions
   leg->SetHeader("Angles");
   leg->AddEntry(hist1,"0","l");
   leg->AddEntry(hist2,"1","l");
   leg->AddEntry(hist3,"2","l");
   leg->AddEntry(hist5,"4","l");
   leg->AddEntry(hist7,"6","l");
   leg->Draw();

   c3->SaveAs("AEDistrib" + Znum + ".png");
}

   in.close();

   f->Write();
}

void DoTDeriMax1090Correction(TString SpectrumFileInput = "/lustre/miifs05/scratch/him-specf/hyp/steinen/COSYBeamtestAna/COSYnewMogon/June2014/COSYJune2014Dataset11_200,100,0,5339_SR1.root", TString FitFileInput = "/lustre/miifs05/scratch/him-specf/hyp/steinen/COSYBeamtestAna/COSYnewMogon/Fit/FitCOSYJune2014Dataset11_200,100,0,5339_SR1.root")
{
	TH2D *hSpectrumTDeriMax1090_EnergyChannel;
	TH2D *hSpectrumTDeriMax1090Rel_EnergyChannel;
	
	TH2D *hSpectrumT1090_EnergyChannelCorr1;
	TNtuple *DataNTuple;
	TFile *SpectrumInput = new TFile(SpectrumFileInput.Data());
	
	hSpectrumTDeriMax1090_EnergyChannel = (TH2D*) SpectrumInput->Get("Histograms/Energy_DeriMaxT90/Energy_DeriMaxT90_01");
	hSpectrumTDeriMax1090_EnergyChannel->SetDirectory(0);
	hSpectrumTDeriMax1090Rel_EnergyChannel = (TH2D*) SpectrumInput->Get("Histograms/Energy_DeriMaxT90Rel/Energy_DeriMaxT90Rel_01");
	hSpectrumTDeriMax1090Rel_EnergyChannel->SetDirectory(0);
	hSpectrumT1090_EnergyChannelCorr1 = (TH2D*) SpectrumInput->Get("Histograms/EnergyRt1090/EnergyRt1090CorrectionRt_01");
	hSpectrumT1090_EnergyChannelCorr1->SetDirectory(0);
	SpectrumInput->Close();
	//hSpectrumTDeriMax1090_EnergyChannel->Draw("colz");
	
	
	TFile *FitInput = new TFile(FitFileInput.Data(),"Update");
	DataNTuple = (TNtuple*)FitInput->Get("DataNTuple");
	DataNTuple->Scan();
	
	Int_t entries = (Int_t)DataNTuple->GetEntries();
	cout<<"Number of Entries: "<<entries<< endl;
	const int entriesArrayValue =entries;
	TF1 *FitFunc[entriesArrayValue];
	float Energy=0;
	float ChannelPeakPos=0;
	float ChannelRangeMin=0;
	float ChannelRangeMax=0;
	DataNTuple->SetBranchAddress("Energy",&Energy);
	DataNTuple->SetBranchAddress("ChannelPeakPos",&ChannelPeakPos);
	DataNTuple->SetBranchAddress("ChannelRangeMin",&ChannelRangeMin);
	DataNTuple->SetBranchAddress("ChannelRangeMax",&ChannelRangeMax);
	
	TCanvas* can=new TCanvas();
	gPad->SetLogz();
	hSpectrumTDeriMax1090_EnergyChannel->GetXaxis()->SetRangeUser(-30,250);
	hSpectrumTDeriMax1090Rel_EnergyChannel->GetXaxis()->SetRangeUser(-0.1,0.95);
	
	for (Int_t ki=0;ki<entries;ki++)
	{
		DataNTuple->GetEntry(ki);
		//if (int(Energy) == 1332)
		//if (int(Energy) == 510)
		//if (ki == entries-1)
		{
			cout << ChannelRangeMin << " " << ChannelRangeMax << endl;
			
			
			//first correction via TDeriMaxT90Rel
			hSpectrumTDeriMax1090_EnergyChannel->GetYaxis()->SetRangeUser(ChannelRangeMin,ChannelRangeMax);
			hSpectrumTDeriMax1090Rel_EnergyChannel->GetYaxis()->SetRangeUser(ChannelRangeMin,ChannelRangeMax);
			//TF1* FitFuncSlices = new TF1("FitFuncSlices","gaus(0)+pol0(3)",ChannelRangeMin,ChannelRangeMax);
			//FitFuncSlices->SetParameters(1000,ChannelPeakPos-10,4,10);
			//FitFuncSlices->SetParLimits(1,ChannelRangeMin,ChannelRangeMax);
			//FitFuncSlices->SetParLimits(2,0,10);
			//FitFuncSlices->SetParLimits(3,0,100);
	
			//gDirectory->ls();
			
			TH1D *hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually=new TH1D("hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually","",hSpectrumTDeriMax1090Rel_EnergyChannel->GetNbinsX(),-0.3,1.3);
			//cout <<hSpectrumTDeriMax1090_EnergyChannel_MaxPos->GetEntries()<< endl;
			for(int binX = hSpectrumTDeriMax1090Rel_EnergyChannel->GetXaxis()->FindBin(-0.1);binX <= hSpectrumTDeriMax1090Rel_EnergyChannel->GetXaxis()->FindBin(0.90);binX++)
			{
				cout << "binx " << binX << endl;
				TH1D *hProfileY =hSpectrumTDeriMax1090Rel_EnergyChannel->ProjectionY("_py",binX,binX);	
				double MaxValue=hProfileY->GetBinCenter(hProfileY->GetMaximumBin());
				TF1* FitFuncSlices = new TF1("FitFuncSlices","gaus(0)+[3]",MaxValue-20,MaxValue+20);
				TF1* FitFuncGausSlices = new TF1("FitFuncGausSlices","gaus(0)",MaxValue-20,MaxValue+20);
				FitFuncGausSlices->SetParameters(hProfileY->GetBinContent(hProfileY->GetMaximumBin()),MaxValue,4);
				
				hProfileY->Fit(FitFuncGausSlices,"RNQ");
				
				FitFuncSlices->SetParameters(FitFuncGausSlices->GetParameter(0),FitFuncGausSlices->GetParameter(1),FitFuncGausSlices->GetParameter(2),10);
				
				FitFuncSlices->SetParLimits(0,0,10000);
				FitFuncSlices->SetParLimits(1,MaxValue-10,MaxValue+10);
				FitFuncSlices->SetParLimits(2,0,10);
				FitFuncSlices->SetParLimits(3,0,100);
				hProfileY->Fit(FitFuncSlices,"RNQ");
				
				cout <<MaxValue<<"  " << FitFuncSlices->GetParameter(1) << "   " << FitFuncSlices->GetParError(1) <<endl;
				hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->SetBinContent(binX, FitFuncSlices->GetParameter(1));
				hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->SetBinError(binX, FitFuncSlices->GetParError(1));
				
				
			}
			hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->GetYaxis()->SetRangeUser(ChannelPeakPos-100,ChannelPeakPos+50);
			hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->GetXaxis()->SetRangeUser(-0.05,0.9);
			TF1 * funcCorrConst = new TF1("funcCorrConst","pol0",-0.03,0.03);
			funcCorrConst->SetLineColor(kRed);
			TF1 * funcCorrPol = new TF1("funcCorrPol",PolPiecewise,-0.05,0.9,6);
			
			funcCorrPol->SetLineColor(kBlue);
			funcCorrPol->SetParameter(0,0.04);
			funcCorrPol->SetParameter(1,ChannelPeakPos);
			hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->Fit(funcCorrConst,"R");
			hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->Fit(funcCorrPol,"R+");
//.........这里部分代码省略.........

void TakeDataAnalyzer(Int_t RunN){
 
  //Int_t RunN = 4591;   
  //string FileRoot = "/d00/icali/PixelHwStudies/PSIStuff/log/bt05r";
 string FileRoot = "/home/l_tester/log/bt05r";
 string FileInName = "thistos.root";
 string FileOutName = "ProcessedData.root";
 char RunNumber[6];
 sprintf(RunNumber, "%.6d", RunN);
 FileRoot = FileRoot + RunNumber + "/";
 FileInName = FileRoot + FileInName;
 FileOutName = FileRoot + FileOutName;
 
 TFile *iFile = new TFile((const char*)FileInName.c_str());
 TNtuple *Digis = (TNtuple*)iFile->FindObjectAny("events");
 Int_t EvN, row, col, roc, ph;
 Float_t vcal;
 Digis->SetBranchAddress("row", &row);
 Digis->SetBranchAddress("col", &col);
 Digis->SetBranchAddress("roc", &roc);
 Digis->SetBranchAddress("ph", &ph);
 Digis->SetBranchAddress("vcal", &vcal);
 Digis->SetBranchAddress("eventNr", &EvN);
 Int_t iNTlenght = Digis->GetEntries();
 
 TFile *oFile = new TFile((const char*)FileOutName.c_str(), "RECREATE");
 TNtuple *Multiplicity = new TNtuple("ModMult", "ModMult", "EvN:roc:NModHits:NROCHits:NRows:NCols:MeanPh:MeanPhPuls",10000000);
 
 //******************Calculating Number of Events**************
 Digis->GetEntry(iNTlenght -1);
 std::cout << " MaxEvent " << EvN << std::endl;
  
 //*****************Creating Multiplicity NTuple***************
 
 ModuleMultiplicity  ModMult;
 Int_t OldEvN = 1;
 for(Int_t i =0; i< iNTlenght; ++i){
    Digis->GetEntry(i);
    if(OldEvN != EvN ||( i == iNTlenght -1) ){
      for(Int_t j=0; j < 16; ++j){ 
        Float_t NModHits = ModMult.ModMultiplicity();
        Float_t NROCHits = ModMult.ROCMultiplicity(j);
        Float_t NRows = ModMult.NumberOfRows(j);
        Float_t NCols = ModMult.NumberOfColumns(j);
        Float_t MeanPh = ModMult.MeanPh(j);  
        Float_t MeanPhPuls = ModMult.MeanPhPuls(j); 
        if(NROCHits >0) Multiplicity->Fill(EvN-1, j,NModHits ,NROCHits ,NRows ,NCols ,MeanPh, MeanPhPuls); 
      }
      ModMult.Clean();
      OldEvN = EvN;
    }
    
    bool SelPixel = 0;
    if(row ==5 && col ==5)SelPixel = 1;
    if(roc>-1&&col>-1)ModMult.Add(roc, row, col, ph, SelPixel);
 }
 
 oFile->Write();
 oFile->Close();          
 
}

void gammaJetDifferentialxjg()
{
  TH1::SetDefaultSumw2();

  TFile *inFile = TFile::Open("gammaJets_inclusive_dphi7pi8_pPbData_v2.root");
  TNtuple *inTuple = (TNtuple*)inFile->Get("gammaJets");

  //book histos
  TH1D *gammaJetEta[nhfBins];
  TH1D *x_jg[nhfBins];
  TH1D *hfEnergy[nhfBins];
  Int_t numEvts[nhfBins];
  Int_t numEvtsEtaPlus[nhfBins];
  
  TCanvas *c[4];
  c[0] = new TCanvas();
  c[1] = new TCanvas();
  TLegend *leg = new TLegend(0.65,0.65,0.9,0.9);
  leg->SetFillColor(0);
  
  for(Int_t i = 0; i < nhfBins; ++i)
  {
    TString s_gammaJetEta = "gammaJetEta";
    s_gammaJetEta += i;
    gammaJetEta[i] = new TH1D(s_gammaJetEta,";Avg #eta_{#gamma jet};Event Fraction",14,-3,3);

    TString s_hfEnergy = "hfEnergy";
    s_hfEnergy += i;
    hfEnergy[i] = new TH1D(s_hfEnergy,"",1000,0,1000);

    TString s_x_jg = "x_jg";
    s_x_jg += i;
    x_jg[i] = new TH1D(s_x_jg,";p_{T}^{jet}/p_{T}^{#gamma};Event Fraction",20,0,1.5);

    TCut isolationCut = "(cc4 + cr4 + ct4PtCut20 < 1.0)";
    //TCut etaCut = "(abs(gEta) > 1.479)";
    //TCut showerShapeCut = "(sigmaIetaIeta < 0.035)";
    TCut etaShowerCut = "((abs(gEta) < 1.479) && (sigmaIetaIeta < 0.01) || (abs(gEta) > 1.479) && (sigmaIetaIeta < 0.035))";
    TCut ptCut = "(jPt > 30) && (gPt > 50)";
    TString hfCut = Form("((HFplusEta4+HFminusEta4) > %f) && ((HFplusEta4+HFminusEta4) > %f)", hfBins[i], hfBins[i+1]);
    //TCut cut = isolationCut && etaCut && showerShapeCut && ptCut && hfCut;
    TCut cut = isolationCut && etaShowerCut && ptCut && hfCut;
    //printf("cut: %s\n",(const char*)cut);
    TString runFlip = "((run > 211257)*(-1) + (run < 211257))";
    
    numEvts[i] = inTuple->Project(s_gammaJetEta,"avgEta*"+runFlip,cut);
    numEvtsEtaPlus[i] = inTuple->GetEntries(
      cut &&
      "( avgEta*(run < 211257) >0 ) || ( avgEta*(run > 211257) <0 )"
      );
    inTuple->Project(s_hfEnergy,"(HFplusEta4+HFminusEta4)",cut);
    inTuple->Project(s_x_jg,"jPt/gPt",cut);

    gammaJetEta[i]->Scale(1./numEvts[i]);
    gammaJetEta[i]->SetMarkerColor(i+1);
    gammaJetEta[i]->SetLineColor(i+1);

    x_jg[i]->Scale(1./numEvts[i]);
    x_jg[i]->SetMarkerColor(i+1);
    x_jg[i]->SetLineColor(i+1);

    TString label;
    label += hfBins[i];
    label += " < E_{T}^{HF[|#eta|>4]} < ";
    label += hfBins[i+1];
    
    leg->AddEntry(gammaJetEta[i], label, "lp");
    
    if(i==0)
    {
      //TLatex *lnorm = new TLatex(0.2,0.85, "Monte Carlo");
      //lnorm->SetNDC(1);
      //lnorm->SetTextSize(0.05);
     
      gammaJetEta[i]->GetXaxis()->CenterTitle();
      gammaJetEta[i]->GetYaxis()->CenterTitle();
      gammaJetEta[i]->SetMaximum(gammaJetEta[i]->GetMaximum()*1.7);
      x_jg[i]->GetXaxis()->CenterTitle();
      x_jg[i]->GetYaxis()->CenterTitle();
      x_jg[i]->SetMaximum(x_jg[i]->GetMaximum()*1.7);

      c[0]->cd();
      gammaJetEta[i]->Draw("");
      //lnorm->Draw("same");
      c[1]->cd();
      x_jg[i]->Draw("");
      //lnorm->Draw("same");
    }
    else
    {
      c[0]->cd();
      gammaJetEta[i]->Draw("same");
      c[1]->cd();
      x_jg[i]->Draw("same");
    }

    c[0]->cd();
    gammaJetEta[i]->Draw("same Lhist");
    leg->Draw();
    c[1]->cd();
//.........这里部分代码省略.........

void jettrigger()
{
  TFile *f = new TFile("jettrig.root");
  TNtuple *nt = (TNtuple *)f->Get("nt");

//Declaration of leaves types
   Float_t         pt;
   Float_t         eta;
   Float_t         phi;
   Float_t         mass;
   Float_t         jt40;
   Float_t         jt60;
   Float_t         jt80;
   Float_t         jt100;
   Float_t         pscl40;
   Float_t         pscl60;
   Float_t         pscl80;
   Float_t         pscl100;
  

   // Set branch addresses.
   nt->SetBranchAddress("pt",&pt);
   nt->SetBranchAddress("eta",&eta);
   nt->SetBranchAddress("phi",&phi);
   nt->SetBranchAddress("mass",&mass);
   nt->SetBranchAddress("jt40",&jt40);
   nt->SetBranchAddress("jt60",&jt60);
   nt->SetBranchAddress("jt80",&jt80);
   nt->SetBranchAddress("jt100",&jt100);
   nt->SetBranchAddress("pscl40",&pscl40);
   nt->SetBranchAddress("pscl60",&pscl60);
   nt->SetBranchAddress("pscl80",&pscl80);
   nt->SetBranchAddress("pscl100",&pscl100);
  
   TFile *fout = new TFile("jettrig_withweight.root","recreate");
   TNtuple *ntout = new TNtuple("ntweight","ntweight","pt:eta:phi:mass:jt40:jt60:jt80:jt100:pscl40:pscl60:pscl80:weightJet:weight12003");
   //TTree *ntout=new TTree("ntweight","ntweight");
   ntout->SetDirectory(fout);

   /*   ntout->Branch("pt",&pt);
   ntout->Branch("eta",&eta);
   ntout->Branch("phi",&phi);
   ntout->Branch("mass",&mass);
   ntout->Branch("jt40",&jt40);
   ntout->Branch("jt60",&jt60);
   ntout->Branch("jt80",&jt80);
   ntout->Branch("jt100",&jt100);
   ntout->Branch("pscl40",&pscl40);
   ntout->Branch("pscl60",&pscl60);
   ntout->Branch("pscl80",&pscl80);
   ntout->Branch("pscl100",&pscl100);
   ntout->Branch("weight",&weight);
   */
   Long64_t nentries = nt->GetEntries();
   cout<<nentries<<endl;

   int oneperc = nentries/100;
   Long64_t nbytes = 0;

   //#ifndef __CINT__
   //#pragma omp parallel for ordered schedule(dynamic)
   //#endif
   for (Long64_t i=0; i<nentries;i++) {
           nbytes += nt->GetEntry(i);

      if (i % oneperc == 0) cout<<"\r"<<i/oneperc<<"%   "<<flush;

      float weightJet = 0, weight12003 = 0;

      if (jt40 && !jt60 && !jt80 && !jt100) weight12003 = 1/pscl40;
      if (jt60 && !jt80 && !jt100) weight12003 = 1;
      if (jt80 && !jt100) weight12003 = 1;
      if (jt100) weight12003 = 1;

      if (jt40 && pt>40 && pt<60) weightJet = pscl40;
      if (jt60 && pt>60 && pt<80) weightJet = pscl60;
      if (jt80 && pt>80 && pt<100) weightJet = pscl80;
      if (jt100 && pt>100) weightJet = pscl100;

      //#ifndef __CINT__
      //#pragma omp ordered
      //#endif
	ntout->Fill(pt,eta,phi,mass,jt40,jt60,jt80,jt100,pscl40,pscl60,pscl80,weightJet,weight12003);

   }
   cout<<endl;
   cout<<ntout->GetEntries()<<endl;
   ntout->Write();
   fout->Close();
   f->Close();

   f = new TFile("jettrig_withweight.root");
   nt = (TNtuple *)f->Get("ntweight");
   cout<<nt->GetEntries()<<endl;
   f->Close();

}

void fragmentYieldsPlot() {
   gStyle->SetOptStat(0000000000); //remove the for this graphs totally redundant statbox
   int ZNumGiven;
   cout << "Enter fragment Z-number (eg. 1): ";
   cin >> ZNumGiven;

   TCanvas *c1 = new TCanvas("fragmentYieldsPlot", "Total yield of fragments zero to ten degrees as function of depth");
   
   TString fragmentNameChoices[6] = {"H","He","Li","Be","B","C"};
   
   TString fragmentName = fragmentNameChoices[ZNumGiven-1];  
   
   std::cout << fragmentName << endl;
   
   TH1F* dummyHisto = new TH1F("dummyHisto", fragmentName + " yields 0-10 degrees" ,100, 0.0,40); //Dummyhisto fix for missing TNtuple methods.
   dummyHisto->SetXTitle("Depth (cm)");
   dummyHisto->SetYTitle("N/N0");

   ifstream in;
   TString experimentalDataPath = "experimentalData/iaeaBenchmark/yields/TDK" + fragmentName + ".dat";
   
   ifstream in;

   //Pull in ascii/exfor-style data
   in.open(experimentalDataPath);

   Float_t f1,f2;
   Int_t nlines = 0;
   TFile *f = new TFile("fragmentAngularDistribution.root","RECREATE");
   TNtuple *ntuple = new TNtuple("ntuple","Data from ascii file","x:y");
	  
   Char_t DATAFLAG[4];
   Int_t NDATA;
   Char_t n1[15], n2[15];
   in >> DATAFLAG >> NDATA ; // Read EXFOR line: 'DATA 6'
   in >> n1 >> n2; // Read  column titles: 'Energy He B [...]'

   cout <<n1<<"   "<<n2<<"\n";
   while (1) {
      in >> f1 >> f2;
      if (!in.good()) break;
      if (nlines < 500 ) printf("%f %f\n",f1,f2);
      ntuple->Fill(f1,f2);
      nlines++;
   }
   std::cout << "Imported " << nlines << " lines from data-file" << endl;




   TNtuple *simData = new TNtuple("ntuple","Data from ascii file","depth:H:He:Li:Be:B:C");
   
//   gROOT->SetStyle("clearRetro");
 //this will be used as base for pulling the experimental data
   TString dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
   dir.ReplaceAll("fragmentAngularDistribution.C","");
   dir.ReplaceAll("/./","/");
   ifstream in;
   simData->Fill(0.0,0.0,0.0,0.0,0.0,0.0,1.0);
for(int j = 1; j <= 40;j=j=j+1){
   TString pDepth, fragment, Znum, normToOneAtZeroAngle;
   pDepth = Form("%i",j);
/*
   cout << "Enter phantom depth (eg. 27.9, see experimentalData directory for choices): ";
   cin >> pDepth;
*/
   TString simulationDataPath = "IAEA_" + pDepth + ".root";

   //Let's pull in the simulation-data
   //TFile *MCData = TFile::Open("IAEA_" + pDepth + ".root");
   TFile *MCData = TFile::Open(simulationDataPath);
   TNtuple *fragments = (TNtuple*) MCData->Get("fragmentNtuple");

   //Block bellow pulls out the simulation's metadata from the metadata ntuple.
   TNtuple *metadata = (TNtuple*) MCData->Get("metaData");
   Float_t events, detectorDistance,waterThickness,beamEnergy,energyError,phantomCenterDistance;
   metadata->SetBranchAddress("events",&events);
   metadata->SetBranchAddress("waterThickness",&waterThickness);
   metadata->SetBranchAddress("detectorDistance",&detectorDistance);
   metadata->SetBranchAddress("beamEnergy",&beamEnergy);
   metadata->SetBranchAddress("energyError",&energyError);
   metadata->SetBranchAddress("phantomCenterDistance",&phantomCenterDistance);
   metadata->GetEntry(0); //there is just one row to consider.
	//ALL UNITS ARE cm!
	Double_t scatteringDistance = detectorDistance - phantomCenterDistance; //temporarily hard-coded, should be distance from target-center to detector

	Double_t degrees = 10.0;
	Double_t r, rMin, rMax, graphMaximum = 0.0;
	Double_t norming = events*.999;
	TString rMinString;
	TString rMaxString;

	rMinString = "0.00";
	rMaxString = Form("%f", scatteringDistance*TMath::ATan(degrees*TMath::DegToRad()));

		Double_t H = ((Double_t*) fragments->GetEntries("(Z == " + TString::Format("%i",1) + "  && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")")) / norming;
		Double_t He = ((Double_t*) fragments->GetEntries("(Z == " + TString::Format("%i",2) + "  && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")")) / norming;
		Double_t Li = ((Double_t*) fragments->GetEntries("(Z == " + TString::Format("%i",3) + "  && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")")) / norming;
		Double_t Be = ((Double_t*) fragments->GetEntries("(Z == " + TString::Format("%i",4) + "  && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")")) / norming;
		Double_t B = ((Double_t*) fragments->GetEntries("(Z == " + TString::Format("%i",5) + "  && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")")) / norming;
//.........这里部分代码省略.........