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

void makeTTree(){

	//	Declaring 
	Double_t num;
	Double_t xNumBins, yNumBins;

	//	Setting up Canvas to put plot the 2D histogram (from Tamii) and the 2 Bracnh Tree (created by this program)
	TCanvas *canvas = new TCanvas("canvas","Canvas");
	canvas->Divide(2,1);
	canvas->cd(1);


	//	Opening up histogram file from tamii's analyzer
	//	Remember to do h2root before this program runs.
	TFile *f1 = new TFile("run6106_test44.root");
	TH2F *hist2D = (TH2F*)f1->Get("h159");
	//	Plotting the 2D histogram on the first part of the Canvas
	hist2D->Draw("COLZ");
	
	//	Going through 2D hist and getting the total number of x and y bins.
	xNumBins = hist2D->GetNbinsX();
	yNumBins = hist2D->GetNbinsY();
	
	//	Creating a new file to save the 2 branch tree to.
	TFile *t1 = new TFile("~/e329/root/13C/runs/run3014.root","recreate");

	//	This is the actual tree being created.
	//	!!! For you YingYing, instead of x vs theta, you would have theta(fp or target) vs Yfp !!!
	TNtuple *DATA = new TNtuple("DATA","Xpos vs Theta","Xpos:Theta");
	

	//	This is the main nested for loop that goes through each bin in the 2D hist 
	//	takes the number of counts in that (x,y)bin and creates that many events
	//	in the 2 branch tree  
	for (int binx = 0; binx < xNumBins+2; binx++){

		for (int biny = 0; biny < yNumBins+2; biny++){

			//getting the number of events in a particular (x,y)bin
			num = hist2D->GetBinContent(binx,biny);
			double x = hist2D->GetXaxis()->GetBinCenter(binx); 
			double y = hist2D->GetYaxis()->GetBinCenter(biny); 
			//cout << num << " at x= " << x << ", y= " << y << endl;

			//Filling the Tree 'num' times with xbin and ybin events
			for (int i = 0; i < num; i++){
				DATA->Fill(x,y);
			}
		}
	}
	
	//Saving the Tree that was just created and filled
	t1->Write();

	//Ploting the Tree in the second part of the canvas.
	canvas->cd(2);
	DATA->Draw("Theta:Xpos>>(2000,-600,600,1000,-3,3)","","COLZ");
}

void parallelcoord() {

   TNtuple *nt = NULL;

   Double_t s1x, s1y, s1z;
   Double_t s2x, s2y, s2z;
   Double_t s3x, s3y, s3z;
   r = new TRandom();;

   new TCanvas("c1", "c1",0,0,800,700);
   gStyle->SetPalette(1);

   nt = new TNtuple("nt","Demo ntuple","x:y:z:u:v:w");

   for (Int_t i=0; i<20000; i++) {
      r->Sphere(s1x, s1y, s1z, 0.1);
      r->Sphere(s2x, s2y, s2z, 0.2);
      r->Sphere(s3x, s3y, s3z, 0.05);

      generate_random(i);
      nt->Fill(r1, r2, r3, r4, r5, r6);

      generate_random(i);
      nt->Fill(s1x, s1y, s1z, s2x, s2y, s2z);

      generate_random(i);
      nt->Fill(r1, r2, r3, r4, r5, r6);

      generate_random(i);
      nt->Fill(s2x-1, s2y-1, s2z, s1x+.5, s1y+.5, s1z+.5);

      generate_random(i);
      nt->Fill(r1, r2, r3, r4, r5, r6);

      generate_random(i);
      nt->Fill(s1x+1, s1y+1, s1z+1, s3x-2, s3y-2, s3z-2);

      generate_random(i);
      nt->Fill(r1, r2, r3, r4, r5, r6);
   }
   nt->Draw("x:y:z:u:v:w","","para",5000);
   TParallelCoord* para = (TParallelCoord*)gPad->GetListOfPrimitives()->FindObject("ParaCoord");
   para->SetDotsSpacing(5);
   TParallelCoordVar* firstaxis = (TParallelCoordVar*)para->GetVarList()->FindObject("x");
   firstaxis->AddRange(new TParallelCoordRange(firstaxis,0.846018,1.158469));
   para->AddSelection("violet");
   para->GetCurrentSelection()->SetLineColor(kViolet);
   firstaxis->AddRange(new TParallelCoordRange(firstaxis,-0.169447,0.169042));
   para->AddSelection("Orange");
   para->GetCurrentSelection()->SetLineColor(kOrange+9);
   firstaxis->AddRange(new TParallelCoordRange(firstaxis,-1.263024,-0.755292));
}

void Example2(){
    gROOT->Reset();
    gROOT->Clear();
    gROOT->SetStyle("Plain");

    gStyle->SetTextSize(0.01908148);
    gStyle->SetTitleFontSize(0.07);
    gStyle->SetOptTitle(1);
    gStyle->SetOptStat(1110);
    gStyle->SetOptFit(1111);
    gStyle->SetTitleXOffset(1.1);
    gStyle->SetTitleYOffset(1.55);
    gStyle->SetPadTopMargin(0.15);
    gStyle->SetPadBottomMargin(0.15);
    gStyle->SetPadLeftMargin(0.15);

    // select the one of the versus : Hf energy, Centrality, # Charged , # Tracks
    int Vselect = 1;

    // number of Trigger you will use , 
    const int Ntr = 2; 
    //if you want to use more than two trigger please specify the Ntr2 for dividivg the canvas
    const int Ntr2 = 1;
    TFile* inf = new TFile("openhlt2.root");

    string triggers[6] = {"L1Tech_BSC_minBias_threshold1.v0","L1Tech_BSC_minBias_threshold2.v0","L1_SingleJet30","L1_TripleJet30","L1_QuadJet15","L1_DoubleJet70"};
    string vers[4] = {"hiHF","hiBin","Ncharged","hiNtracks"};

    double bins [4] = {170,40,240,140};
    double limits[4] = {170000,40,24000,1400};

    TCanvas* c1 = new TCanvas("c1","c1",800,400);
    c1->Divide(Ntr,Ntr2);

    TNtuple* nt = (TNtuple*)inf->Get("HltTree");

    TProfile* p1[Ntr];

    for(unsigned int i =0; i<Ntr ; i++){
        c1->cd(i+1);

        p1[i] = new TProfile("p1",Form(";%s;%s",vers[Vselect].data(),triggers[i].data()),bins[1],0,limits[1]);
        nt->SetAlias("trigger",triggers[i].data());

        nt->SetAlias("versus",vers[Vselect].data());
        nt->Draw("trigger:versus>>p1","","prof");

    }

    c1->Print(Form("%s_vsCent.gif",triggers[Vselect].data()));

}

double getUncertaintyValue(double m0,double m12, int channel, string valuename, TFile * input)
{
  bool DEBUGuncertaintyValue=false;
/*
  string buffer = doubletostr(m0)+doubletostr(m12)+inttostr(channel)+valuename;
  map <std::string, double>::iterator iter;
  iter=UncertaintyValueBuffer.find(buffer);
  if (bufferValue && iter!=UncertaintyValueBuffer.end()) return iter->second;
*/
  double returnvalue2 ;
  double min=-10;
  double max=+10;

  //cout << "v" ;
  string leaveName="";
  if (valuename=="K") leaveName="K";
  if (valuename=="crossSection") {leaveName="crossSection"; min=0; max=100000;};
  if (valuename=="PDFUncertainty") leaveName="relUncPDF";
  if (valuename=="ScaleUncertainty2Q") leaveName="relScaleUnc2Q";
  if (valuename=="ScaleUncertainty12Q") leaveName="relScaleUncHalfQ";
  if (leaveName=="") throw "no leave found for that value";

  string commandstring=leaveName+">>h1";
  string conditionstring= "finalState>" + doubletostr(channel-0.1)  +"&& finalState<"+ doubletostr(channel+0.1) +"&& m0>" + doubletostr(m0-0.5) + " && m0<"+doubletostr(m0+0.5) +" && m12>" + doubletostr(m12-0.5) + " && m12<"+doubletostr(m12+0.5);

  if (DEBUGuncertaintyValue){
  	  cout << "	 commandstring " << commandstring << endl;
  	  cout << "	 conditionstring " << conditionstring << endl;
  	  cout << "	 leaveName " << leaveName << endl;
  	  cout << "	 valuename :" << valuename << " channel: " << channel << " m0: " << m0 << " m12: " << m12   << endl;
  	  }

  TNtuple *ntuple = (TNtuple*)file_robin->Get("SignalUncertainties");
  if (ntuple==0) throw "Ntuple not found";

  if (DEBUGuncertaintyValue) cout << "(min,max) "<<min << " "<<max<<endl;
  TH1F * h1 = new TH1F("h1","h1",100000,min,max);
  ntuple->Draw(commandstring.c_str(),conditionstring.c_str());

  if (DEBUGuncertaintyValue) cout << "mean " << h1->GetMean()<<endl;
  
  returnvalue2 = fabs (h1->GetMean());
//  UncertaintyValueBuffer[buffer]=returnvalue2;
  
  delete h1;

  return returnvalue2;
}

void rootFit(){

  TCanvas *c1 = new TCanvas("c1","",600,600);
  c1->Divide(1,2);
  TNtuple *T = new TNtuple("T","","x:y:z");
  T->ReadFile("histogramData.dat");
  T->Draw("y:x");

  Double_t *X = T->GetV1();
  Double_t *Y = T->GetV2();

  for(int loop = 0; loop < 120; loop++){
    cout << X[loop] << "  " << Y[loop] << endl;
  }

  c1->Clear();
  TGraph *graphT = new TGraph(120,Y,X);
  TGraph *graphB = new TGraph(120,Y,X);
  c1->cd(1);
  graphT->Draw("APL");
  graphT->Fit("gaus+pol2","RME");

}

void unfoldPt(int mode=0)
{

   // Matched Tracklets
   TFile *inf = new TFile("match-10TeV-12.root");
   TNtuple *nt = (TNtuple*)inf->FindObjectAny("nt");

   // Test sample
   TFile *infTest = new TFile("./TrackletTree-Run123596.root");
   TNtuple *ntTest = (TNtuple*)infTest->FindObjectAny("TrackletTree12");
   
   TFile *pdfFile;
   if (mode==0) pdfFile = new TFile("pdf.root","recreate");
           else pdfFile = new TFile("pdf.root");

   double nPtBin=15;
   double minPt=log(0.05);
   double maxPt=log(10);
   double nDphiBin=600;
   double maxDphi=0.4;
   
   char* mycut = Form("abs(eta)<2&&log(pt)>%f&&log(pt)<%f",minPt,maxPt);
   char* mycut1=Form("abs(eta)<2&&log(pt)>%f&&log(pt)<%f&&abs(eta-eta1)<0.01&&abs(deta)<0.01",minPt,maxPt);

   TH2F *h;
   TH1F *hdphi = new TH1F("hdphi","",nDphiBin,0,maxDphi);
   TH1F *hdphi2;
   TH1F *hpt;
   TH1F *hptH = new TH1F("hptH","",nPtBin,minPt,maxPt);
   
   TH1F *hptUnfold = new TH1F("hptUnfold","",nPtBin,minPt,maxPt);
   TH1F *hptMC = new TH1F("hptMC","",nPtBin,minPt,maxPt);
   TH1F *hptTemp = new TH1F("hptTemp","",nPtBin,minPt,maxPt);

   // Delta phi as a function of matched genparticle transverse momentum
   TCanvas *c = new TCanvas("c","",600,600);
   
   if (mode == 0) {
      h = new TH2F("h","",nPtBin,minPt,maxPt,nDphiBin,0,maxDphi);
      hdphi2 = new TH1F("hdphiMC","",nDphiBin,0,maxDphi);
      hpt = new TH1F("hpt","",nPtBin,minPt,maxPt);      
      h->SetXTitle("ln(P_{T}) GeV/c");
      h->SetYTitle("|#Delta#phi|");
      nt->Draw("abs(dphi):log(pt)>>h",mycut1,"col");
      // used to generate pdf
      nt->Draw("abs(dphi)>>hdphiMC",mycut,"");
      nt->Draw("log(pt)>>hpt",mycut,"");
      h->Write();      
      hpt->Write();      
      hdphi2->Write();      
   } else {
      h = (TH2F*) pdfFile->FindObjectAny("h");
      hdphi2 = (TH1F*) pdfFile->FindObjectAny("hdphiMC");
      hpt = (TH1F*) pdfFile->FindObjectAny("hpt");
   }
   // Delta phi fit
   TCanvas *c2 = new TCanvas("c2","",600,600);
   c2->SetLogy();
   c2->SetLogx();

 
   // dphi for unfolding and MC truth:  
   ntTest->Draw("abs(dphi)>>hdphi","abs(eta1)<2&&abs(deta)<0.1","",200000);
   ntTest->Draw("log(pt)>>hptH",mycut,"",200000);
   
   histFunction2D *myfun = new histFunction2D(h);
   
   TF1 *test = new TF1("histFun",myfun,&histFunction2D::evaluate,0,maxDphi,nPtBin+1);
   TF1 *test2 = new TF1("histFunMC",myfun,&histFunction2D::evaluate,0,maxDphi,nPtBin+1);

   for (int i=0;i<nPtBin+1;i++)
   {  
      test->SetParameter(i,1);   
   }


   hdphi2->SetXTitle("|#Delta#phi|");
   hdphi2->SetYTitle("Arbitrary Normalization");
   hdphi2->Fit("histFunMC","M");

   hdphi->SetXTitle("|#Delta#phi|");
   hdphi->SetYTitle("Arbitrary Normalization");
   hdphi->Fit("histFun","M");
   hdphi->SetStats(0);
   hdphi->Draw();

   
   for (int i=0;i<nPtBin+1;i++) {
      TF1 *testPlot = new TF1(Form("histFun%d",i),myfun,&histFunction2D::evaluate,0,maxDphi,nPtBin+1);

      testPlot->SetParameter(i,test->GetParameter(i));
      testPlot->SetLineColor(i+2);
      testPlot->Draw("same");
   }
   
   int total=0,totalMC=0;


   for (int i=0;i<nPtBin;i++){
      if (test->GetParameter(i)==0) continue;
//.........这里部分代码省略.........

void canvas_write()
{
   //just in case this script is executed multiple times
   delete gROOT->GetListOfFiles()->FindObject("hsimple.root");
   delete gROOT->GetListOfCanvases()->FindObject("c1");

   gBenchmark->Start("ntuple1");
   //
   // Connect ROOT histogram/ntuple demonstration file
   // generated by example hsimple.C.
   TFile *f1 = new TFile("hsimple.root");
   //
   // Create a canvas, with 4 pads
   //
   TCanvas *c1 = new TCanvas("c1","The Ntuple canvas",200,10,700,780);
   TPad *pad1 = new TPad("pad1","This is pad1",0.02,0.52,0.48,0.98,21);
   TPad *pad2 = new TPad("pad2","This is pad2",0.52,0.52,0.98,0.98,21);
   TPad *pad3 = new TPad("pad3","This is pad3",0.02,0.02,0.48,0.48,21);
   TPad *pad4 = new TPad("pad4","This is pad4",0.52,0.02,0.98,0.48,1);
   pad1->Draw();
   pad2->Draw();
   pad3->Draw();
   pad4->Draw();
   //
   // Change default style for the statistics box
   gStyle->SetStatW(0.30);
   gStyle->SetStatH(0.20);
   gStyle->SetStatColor(42);
   //
   // Display a function of one ntuple column imposing a condition
   // on another column.
   pad1->cd();
   pad1->SetGrid();
   pad1->SetLogy();
   pad1->GetFrame()->SetFillColor(15);
   TNtuple *ntuple = (TNtuple*)f1->Get("ntuple");
   ntuple->SetLineColor(1);
   ntuple->SetFillStyle(1001);
   ntuple->SetFillColor(45);
   ntuple->Draw("3*px+2","px**2+py**2>1");
   ntuple->SetFillColor(38);
   ntuple->Draw("2*px+2","pz>2","same");
   ntuple->SetFillColor(5);
   ntuple->Draw("1.3*px+2","(px^2+py^2>4) && py>0","same");
   pad1->RedrawAxis();
   //
   // Display the profile of two columns
   // The profile histogram produced is saved in the current directory with
   // the name hprofs
   pad2->cd();
   pad2->SetGrid();
   pad2->GetFrame()->SetFillColor(32);
   ntuple->Draw("pz:px>>hprofs","","goffprofs");
   TProfile *hprofs = (TProfile*)gDirectory->Get("hprofs");
   hprofs->SetMarkerColor(5);
   hprofs->SetMarkerSize(0.7);
   hprofs->SetMarkerStyle(21);
   hprofs->Fit("pol2");
   // Get pointer to fitted function and modify its attributes
   TF1 *fpol2 = hprofs->GetFunction("pol2");
   fpol2->SetLineWidth(4);
   fpol2->SetLineColor(2);
   //
   // Display a scatter plot of two columns with a selection.
   // Superimpose the result of another cut with a different marker color
   pad3->cd();
   pad3->GetFrame()->SetFillColor(38);
   pad3->GetFrame()->SetBorderSize(8);
   ntuple->SetMarkerColor(1);
   ntuple->Draw("py:px","pz>1");
   ntuple->SetMarkerColor(2);
   ntuple->Draw("py:px","pz<1","same");
   //
   // Display a 3-D scatter plot of 3 columns. Superimpose a different selection.
   pad4->cd();
   ntuple->Draw("pz:py:px","(pz<10 && pz>6)+(pz<4 && pz>3)");
   ntuple->SetMarkerColor(4);
   ntuple->Draw("pz:py:px","pz<6 && pz>4","same");
   ntuple->SetMarkerColor(5);
   ntuple->Draw("pz:py:px","pz<4 && pz>3","same");
   TPaveText *l4 = new TPaveText(-0.9,0.5,0.9,0.95);
   l4->SetFillColor(42);
   l4->SetTextAlign(12);
   l4->AddText("You can interactively rotate this view in 2 ways:");
   l4->AddText("  - With the RotateCube in clicking in this pad");
   l4->AddText("  - Selecting View with x3d in the View menu");
   l4->Draw();
   //
   c1->cd();
   c1->Update();
   gStyle->SetStatColor(19);
   gBenchmark->Show("ntuple1");

   TSQLFile* fsql1 = new TSQLFile(dbname, "recreate", username, userpass);
   if (fsql1->IsZombie()) { delete fsql1; return; }

//  changing TSQLFile configuration, you may improve speed
//  of reading or writing object to/from sql database

//   fsql1->SetUseSuffixes(kFALSE);
//.........这里部分代码省略.........

void TestEnergy() {
  TFile *file = TFile::Open("PPCollisionUps1sMilNoTune.root");
  TFile *file2 = TFile::Open("PPCollisionMilNoTune.root");
  TNtuple* PPTupleUps = (TNtuple*)file->Get("PPTuple");
  TNtuple* PPTupleND = (TNtuple*)file2->Get("PPTuple");
  TH1D* NCharUps = new TH1D("NCharUps", "Sum of Energy at mid rapdity", 30000, 0, 300);
  PPTupleUps->Draw("(ETmid-mu1pT-mu2pT)>>NCharUps", "NChar>2 && abs(deta1) < 1.93 && abs(deta2) < 1.93  && mu1pT > 4 && mu2pT > 4");
//  NCharUps->Scale(3.35e-6);
  double sum, percent = 0, centAr[21];
  int c, centa[21], i;
  sum = 0;
  centAr[0] = 180;
  centAr[20] = 0;
  c = 1;
	//Find Activity Ranges
  double min = 0, average[20], bins;
  i = 0;
  //Activity Range Calculations for Track Count
  TH1D* Tracks = new TH1D("Tracks", "Number of Tracks", 10000, 0, 300);
  PPTupleND->Draw("ETmid>>Tracks", "NChar > 0");
  double MinBiasEvents;
  MinBiasEvents = Tracks->GetEntries();
  NCharUps->Draw();
  Tracks->Draw();
  Tracks->Scale(1./Tracks->GetEntries());
  int TrackCent[21];
  percent = .02;
  TrackCent[0] = 300;
  TrackCent[20] = 0;
  c = 1;
  double Ranges[12];
  int counter = 1;
  Ranges[0] = 0;
	//Loop To Get Activity Range Values
  for (i=0; percent <= .02; percent += .02) {
        sum = 0;
  	cout << c << "		";
        for(i=0; sum < percent && i != 10000; i++) {
                sum += Tracks->GetBinContent(10000-i); }
       	 cout << 100*percent - 2 << "-" << 100*percent << "%" << "      " << MinBiasEvents*sum  << "      " << i << " " << Tracks->GetBinCenter(10000-i) << "\n"; // Use to check Activity Range Values
	Ranges[counter] = sum;
	counter++;		
        TrackCent[c] = Tracks->GetBinCenter(10000-i);
        c++;}
  percent = .08;
  for (i=0; percent <= .08; percent += .02) {
        sum = 0;
        cout << c << "		";
        for(i=0; sum < percent || percent == 1; i++) {
                sum += Tracks->GetBinContent(10000-i); }
                         cout << 100*percent - 2 << "-" << 100*percent << "%" << "      " << MinBiasEvents*sum  << "      " << i << " " << Tracks->GetBinCenter(10000-i) << "\n"; // Use to check Activity Range Values
	Ranges[counter] = sum;
        counter++;
        TrackCent[c] = Tracks->GetBinCenter(10000-i);
        c++;}
  percent = .20;
  for (i=0; percent <= .20; percent += .04) {
        sum = 0;
        cout << c << "		";
        for(i=0; sum < percent || percent == 1; i++) {
                sum += Tracks->GetBinContent(10000-i); }
                         cout << 100*percent - 4 << "-" << 100*percent << "%" << "      " << MinBiasEvents*sum  << "      " << i << " " << Tracks->GetBinCenter(10000-i) << "\n"; // Use to check Activity Range Values
	Ranges[counter] = sum;
        counter++;
        TrackCent[c] = Tracks->GetBinCenter(10000-i);
        c++;}
	percent = .30;
//  for (i=0; percent <= .30; percent += .05) {
  //      sum = 0;
//	cout << c << "\n";
//        for(i=0; sum < percent || percent == 1; i++) {
 //               sum += Tracks->GetBinContent(10000-i); }
  //                       cout << 100*percent - 5 << "-" << 100*percent << "%" << "      " << sum  << "      " << i << " " << Tracks->GetBinCenter(10000-i) << "\n"; // Use to check Activity Range Values
 //       TrackCent[c] = Tracks->GetBinCenter(10000-i) + .5;
   //     c++; }
//	percent = .40;
  for (i=0; percent <= 1; percent += .10) {
        sum = 0;
	cout << c << "		";
        for(i=0; sum <= percent && i != 10000; i++) {
                sum += Tracks->GetBinContent(10000-i); }
                         cout << 100*percent - 10 << "-" << 100*percent << "%" << "      " << MinBiasEvents*sum  << "      " << i << " " << Tracks->GetBinCenter(10000-i) << "\n"; // Use to check Activity Range Values
	Ranges[counter] = sum;
        counter++;
        TrackCent[c] = Tracks->GetBinCenter(10000-i);
        c++; }
  cout << "\n" << "Activity Range Values from Array" << "\n";
  for(c = 0; c < 11; c++) {
    	cout << c << ")   "  << TrackCent[c] << "\n"; }
  min = 0;
  double TrackAverage[20];

  cout << "\n" << NCharUps->GetBinCenter(30000) << "\n";

  i = 0;
  double TotalUps = 0;
  //cout << "\n" << "Activity Bin Averages" << "\n";
	//Loop To Get Average In Each Activity Region
  for(c = 0; NCharUps->GetBinCenter(30000-i) > TrackCent[c+1] && NCharUps->GetBinCenter(30000-i) <= TrackCent[c]; c++) {
        sum = 0;
//.........这里部分代码省略.........

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


	//	General Cuts to be applied for various plots

	//	Fit_Status == 3
	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 );

void PFcorr(bool isMC=false) {
  
  gStyle->SetOptStat(0);
  gStyle->SetTitleOffset(1.8,"Y");

  TString filename = "../test/pftuple";
  if(isMC) filename += "MC";
  filename += ".root";

  TFile *f = new TFile(filename.Data());
  TNtuple *nt = (TNtuple*) f->Get("pfCandidateAnalyzer/nt");
  
  TCanvas *c1 = new TCanvas("c1","c1",800,800);
  c1->Divide(2,2);

  TString hadron = "tkptmax>10";
  //hadron += " && type==1";
  nt->SetAlias("trksel","(algo<=7 && nhits>=5 && relpterr<=0.05 && abs(nd0)<3 && abs(ndz)<3)");
  TString seltrk = " && trksel";
  TString fakehadron = hadron + " && fake";

  Float_t axisrange=500.;
  if(!isMC) axisrange=350.;

  c1->cd(1);
  nt->Draw(Form("tkptsum:eetsum+hetsum>>h1(50,0,%f,50,0,%f)",axisrange,axisrange),hadron.Data(),"goff");
  h1->SetTitle("PF candidate with p_{T}>10 GeV/c track element; PF calo cluster E_{T} sum (ECAL+HCAL); Track p_{T} sum");
  h1->Draw("colz");
  if(isMC) {
    nt->Draw(Form("tkptsum:eetsum+hetsum>>h1f(50,0,%f,50,0,%f)",axisrange,axisrange),fakehadron.Data(),"goff");
    TH2F* h1 = h1;
    h1f->Divide(h1); h1f->Scale(100.);
    h1f->Draw("boxsame");
  }
  gPad->SetLogz();

  c1->cd(2);
  nt->Draw(Form("tkptsum:eetsum+hetsum>>h2(50,0,%f,50,0,%f)",axisrange,axisrange),(hadron+seltrk).Data(),"goff");
  h2->SetTitle("PF candidate with p_{T}>10 GeV/c track element; PF calo cluster E_{T} sum (ECAL+HCAL); SELECTED Track p_{T} sum");
  h2->Draw("colz");
  if(isMC) {
    nt->Draw(Form("tkptsum:eetsum+hetsum>>h2f(50,0,%f,50,0,%f)",axisrange,axisrange),(fakehadron+seltrk).Data(),"goff");
    h2f->Draw("boxsame");
  }
  gPad->SetLogz();

  c1->cd(3);
  nt->Draw(Form("tkptmax:eetsum+hetsum>>h3(50,0,%f,50,0,%f)",axisrange,axisrange),hadron.Data(),"goff");
  h3->SetTitle("PF candidate with p_{T}>10 GeV/c track element; PF calo cluster E_{T} sum (ECAL+HCAL); Track with MAX p_{T}");
  h3->Draw("colz");
  if(isMC) {
    nt->Draw(Form("tkptmax:eetsum+hetsum>>h3f(50,0,%f,50,0,%f)",axisrange,axisrange),fakehadron.Data(),"goff");
    h3f->Draw("boxsame");
  }
  gPad->SetLogz();

  c1->cd(4);
  nt->Draw(Form("tkptmax:eetsum+hetsum>>h4(50,0,%f,50,0,%f)",axisrange,axisrange),(hadron+seltrk).Data(),"goff");
  h4->SetTitle("PF candidate with p_{T}>10 GeV/c track element; PF calo cluster E_{T} sum (ECAL+HCAL); SELECTED Track with MAX p_{T}");
  h4->Draw("colz");
  if(isMC) {
    nt->Draw(Form("tkptmax:eetsum+hetsum>>h4f(50,0,%f,50,0,%f)",axisrange,axisrange),(fakehadron+seltrk).Data(),"goff");
    h4f->Draw("boxsame");
  }
  gPad->SetLogz();

}

void parallelcoordtrans() {


   Double_t x,y,z,u,v,w,a,b,c;
   Double_t s1x, s1y, s1z;
   Double_t s2x, s2y, s2z;
   Double_t s3x, s3y, s3z;
   r = new TRandom();;

   TCanvas *c1 = new TCanvas("c1", "c1",0,0,900,1000);
   c1->Divide(1,2);

   if (gVirtualX && !gVirtualX->InheritsFrom("TGCocoa")) {
      std::cout<<"This macro works only on MacOS X with --enable-cocoa\n";
      delete c1;
      delete r;
      return;
   }

   TNtuple *nt = new TNtuple("nt","Demo ntuple","x:y:z:u:v:w:a:b:c");

   int n=0;
   for (Int_t i=0; i<1500; i++) {
      r->Sphere(s1x, s1y, s1z, 0.1);
      r->Sphere(s2x, s2y, s2z, 0.2);
      r->Sphere(s3x, s3y, s3z, 0.05);

      generate_random(i);
      nt->Fill(r1, r2, r3, r4, r5, r6, r7, r8, r9);
      n++;

      generate_random(i);
      nt->Fill(s1x, s1y, s1z, s2x, s2y, s2z, r7, r8, r9);
      n++;

      generate_random(i);
      nt->Fill(r1, r2, r3, r4, r5, r6, r7, s3y, r9);
      n++;

      generate_random(i);
      nt->Fill(s2x-1, s2y-1, s2z, s1x+.5, s1y+.5, s1z+.5, r7, r8, r9);
      n++;

      generate_random(i);
      nt->Fill(r1, r2, r3, r4, r5, r6, r7, r8, r9);
      n++;

      generate_random(i);
      nt->Fill(s1x+1, s1y+1, s1z+1, s3x-2, s3y-2, s3z-2, r7, r8, r9);
      n++;

      generate_random(i);
      nt->Fill(r1, r2, r3, r4, r5, r6, s3x, r8, s3z );
      n++;
   }

   TParallelCoordVar* pcv;

   c1->cd(1);

   // ||-Coord plot without transparency
   nt->Draw("x:y:z:u:v:w:a:b:c","","para");
   TParallelCoord* para1 = (TParallelCoord*)gPad->GetListOfPrimitives()->FindObject("ParaCoord");
   para1->SetLineColor(25);
   pcv = (TParallelCoordVar*)para1->GetVarList()->FindObject("x"); pcv->SetHistogramHeight(0.);
   pcv = (TParallelCoordVar*)para1->GetVarList()->FindObject("y"); pcv->SetHistogramHeight(0.);
   pcv = (TParallelCoordVar*)para1->GetVarList()->FindObject("z"); pcv->SetHistogramHeight(0.);
   pcv = (TParallelCoordVar*)para1->GetVarList()->FindObject("a"); pcv->SetHistogramHeight(0.);
   pcv = (TParallelCoordVar*)para1->GetVarList()->FindObject("b"); pcv->SetHistogramHeight(0.);
   pcv = (TParallelCoordVar*)para1->GetVarList()->FindObject("c"); pcv->SetHistogramHeight(0.);
   pcv = (TParallelCoordVar*)para1->GetVarList()->FindObject("u"); pcv->SetHistogramHeight(0.);
   pcv = (TParallelCoordVar*)para1->GetVarList()->FindObject("v"); pcv->SetHistogramHeight(0.);
   pcv = (TParallelCoordVar*)para1->GetVarList()->FindObject("w"); pcv->SetHistogramHeight(0.);


   // ||-Coord plot with transparency
   TColor *col26 = gROOT->GetColor(26); col26->SetAlpha(0.01);
   c1->cd(2);
   nt->Draw("x:y:z:u:v:w:a:b:c","","para");
   TParallelCoord* para2 = (TParallelCoord*)gPad->GetListOfPrimitives()->FindObject("ParaCoord");
   para2->SetLineColor(26);
   pcv = (TParallelCoordVar*)para2->GetVarList()->FindObject("x"); pcv->SetHistogramHeight(0.);
   pcv = (TParallelCoordVar*)para2->GetVarList()->FindObject("y"); pcv->SetHistogramHeight(0.);
   pcv = (TParallelCoordVar*)para2->GetVarList()->FindObject("z"); pcv->SetHistogramHeight(0.);
   pcv = (TParallelCoordVar*)para2->GetVarList()->FindObject("a"); pcv->SetHistogramHeight(0.);
   pcv = (TParallelCoordVar*)para2->GetVarList()->FindObject("b"); pcv->SetHistogramHeight(0.);
   pcv = (TParallelCoordVar*)para2->GetVarList()->FindObject("c"); pcv->SetHistogramHeight(0.);
   pcv = (TParallelCoordVar*)para2->GetVarList()->FindObject("u"); pcv->SetHistogramHeight(0.);
   pcv = (TParallelCoordVar*)para2->GetVarList()->FindObject("v"); pcv->SetHistogramHeight(0.);
   pcv = (TParallelCoordVar*)para2->GetVarList()->FindObject("w"); pcv->SetHistogramHeight(0.);
}

void ntuple1() {
   //Small tree analysis script
   // To see the output of this macro, click begin_html <a href="gif/ntuple1.gif">here</a> end_html
   //Author:: Rene Brun
   
   //just in case this script is executed multiple times
   delete gROOT->GetListOfFiles()->FindObject("hsimple.root");
   delete gROOT->GetListOfCanvases()->FindObject("c1");

   gBenchmark->Start("ntuple1");
   //
   // Connect ROOT histogram/ntuple demonstration file
   // generated by example $ROOTSYS/tutorials/hsimple.C.
   TFile *f1 = TFile::Open("hsimple.root");
   if (!f1) return;
   //
   // Create a canvas, with 4 pads
   //
   TCanvas *c1 = new TCanvas("c1","The Ntuple canvas",200,10,700,780);
   TPad *pad1 = new TPad("pad1","This is pad1",0.02,0.52,0.48,0.98,21);
   TPad *pad2 = new TPad("pad2","This is pad2",0.52,0.52,0.98,0.98,21);
   TPad *pad3 = new TPad("pad3","This is pad3",0.02,0.02,0.48,0.48,21);
   TPad *pad4 = new TPad("pad4","This is pad4",0.52,0.02,0.98,0.48,1);
   pad1->Draw();
   pad2->Draw();
   pad3->Draw();
   pad4->Draw();
   //
   // Change default style for the statistics box
   gStyle->SetStatW(0.30);
   gStyle->SetStatH(0.20);
   gStyle->SetStatColor(42);
   //
   // Display a function of one ntuple column imposing a condition
   // on another column.
   pad1->cd();
   pad1->SetGrid();
   pad1->SetLogy();
   pad1->GetFrame()->SetFillColor(15);
   TNtuple *ntuple = (TNtuple*)f1->Get("ntuple");
   ntuple->SetLineColor(1);
   ntuple->SetFillStyle(1001);
   ntuple->SetFillColor(45);
   ntuple->Draw("3*px+2","px**2+py**2>1");
   ntuple->SetFillColor(38);
   ntuple->Draw("2*px+2","pz>2","same");
   ntuple->SetFillColor(5);
   ntuple->Draw("1.3*px+2","(px^2+py^2>4) && py>0","same");
   pad1->RedrawAxis();
   //
   // Display the profile of two columns
   // The profile histogram produced is saved in the current directory with
   // the name hprofs
   pad2->cd();
   pad2->SetGrid();
   pad2->GetFrame()->SetFillColor(32);
   ntuple->Draw("pz:px>>hprofs","","goffprofs");
   TProfile *hprofs = (TProfile*)gDirectory->Get("hprofs");
   hprofs->SetMarkerColor(5);
   hprofs->SetMarkerSize(0.7);
   hprofs->SetMarkerStyle(21);
   hprofs->Fit("pol2");
   // Get pointer to fitted function and modify its attributes
   TF1 *fpol2 = hprofs->GetFunction("pol2");
   fpol2->SetLineWidth(4);
   fpol2->SetLineColor(2);
   //
   // Display a scatter plot of two columns with a selection.
   // Superimpose the result of another cut with a different marker color
   pad3->cd();
   pad3->GetFrame()->SetFillColor(38);
   pad3->GetFrame()->SetBorderSize(8);
   ntuple->SetMarkerColor(1);
   ntuple->Draw("py:px","pz>1");
   ntuple->SetMarkerColor(2);
   ntuple->Draw("py:px","pz<1","same");
   //
   // Display a 3-D scatter plot of 3 columns. Superimpose a different selection.
   pad4->cd();
   ntuple->Draw("pz:py:px","(pz<10 && pz>6)+(pz<4 && pz>3)");
   ntuple->SetMarkerColor(4);
   ntuple->Draw("pz:py:px","pz<6 && pz>4","same");
   ntuple->SetMarkerColor(5);
   ntuple->Draw("pz:py:px","pz<4 && pz>3","same");
   TPaveText *l4 = new TPaveText(-0.9,0.5,0.9,0.95);
   l4->SetFillColor(42);
   l4->SetTextAlign(12);
   l4->AddText("You can interactively rotate this view in 2 ways:");
   l4->AddText("  - With the RotateCube in clicking in this pad");
   l4->AddText("  - Selecting View with x3d in the View menu");
   l4->Draw();
   //
   c1->cd();
   c1->Update();
   gStyle->SetStatColor(19);
   gBenchmark->Show("ntuple1");
}

//.........这里部分代码省略.........
  bool draw95CLlimitModel = true;
  bool drawcan = false;
  RooPlot** frame = new RooPlot*[20];//define an array of pointer;frame[i] is a pointer
  TCanvas* thisCan = new TCanvas("thisCan","MyCanvas",1200,(numCat+1)*300); 
  TCanvas* extraCan = new TCanvas("extraCan","extraCan",800,600); 
  thisCan->Divide(3,numCat+1);
  gStyle->SetMarkerSize(0.5);
  //---------count option------------------
  bool enableCount = true;
  //============category boundaries===============
  std::stringstream mvaCut;
  std::stringstream mvaCutLabel;
  TString SIGstring;
  TString BGstring;
  TString ObsDatastring;
  for(int i=0; i<numCat+1; ++i) {
    if(i!= numCat) {
      mvaCut.str("");
      mvaCut << " && mva > "<<catVals[i]<<" && mva < "<<catVals[i+1];//cut values for cat i
      std::cout<<mvaCut.str().c_str()<<std::endl;  
      mvaCutLabel.str("");
      mvaCutLabel << " mva > "<<catVals[i]<<" && mva < "<<catVals[i+1];//cut values for cat i 
    }
    if(i== numCat) {
     mvaCutLabel.str("");
     mvaCutLabel << " Dijet-TAG ";
    }
    //-----obsData for count-----
    //read in the obsdata mass hist for category i
    std::stringstream pSS; 
    pSS.str("");
    pSS << "ObsDataH" <<i;    
    ObsDataH[i] = new TH1D(pSS.str().c_str(),"",320,100.,180.);//0.25 GeV per bin
    TString theDrawString0 = TString("mass>>")+TString(pSS.str().c_str());
    if(i!=numCat){
      ObsDatastring = baseObsData+TString(mvaCut.str().c_str())+TString(" && vbftag==0");
    }
    if(i==numCat){
      ObsDatastring = baseObsData+TString(" && vbftag==1 && mva>0.05");
    }
    extraCan->cd();
    tup->Draw(theDrawString0.Data(),ObsDatastring.Data());
    numObsData[i] = ObsDataH[i]->GetEntries();//number of events in the signal mass hist
     
    //-----signal model-----
    //read in the signal mass hist for category i
    pSS.str("");
    pSS << "sigH" <<i;    
    sigH[i] = new TH1D(pSS.str().c_str(),"",320,100.,180.);//0.25 GeV per bin
    TString theDrawString = TString("mass>>")+TString(pSS.str().c_str());
    //TString SIGstring = TString("weight*(")+baseSIG+TString(mvaCut.str().c_str())+TString(")");//signal cut along with cat cut
    if(i!=numCat){
      SIGstring = TString("weight*(")+baseSIG+TString(mvaCut.str().c_str())+TString(" && vbftag==0")+TString(")");
    }
    if(i==numCat){
      SIGstring = TString("weight*(")+baseSIG+TString(" && vbftag==1")+TString(" && mva>0.05")+TString(")");
    }
    extraCan->cd();
    tup->Draw(theDrawString.Data(),SIGstring.Data());
    //sigH[i] = (TH1D*)gPad->GetPrimitive(pSS.str().c_str());
    nominalSignal[i] = sigH[i]->GetSumOfWeights();//number of events in the signal mass hist
    //obtain the signal RooDataHist to be fit for category i
    pSS.str("");
    pSS << "sigdata" << i;    
    sigNames[i]    = TString(pSS.str().c_str());
    sigCat[i]      = new RooDataHist( sigNames[i].Data(),"",*mass,sigH[i]);//hist used to extract signal model for category i

//.........这里部分代码省略.........
   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;
		Double_t C = ((Double_t*) fragments->GetEntries("(Z == " + TString::Format("%i",6) + "  && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")")) / norming;
		simData->Fill(waterThickness,H,He,Li,Be,B,C);

	}
	simData->Scan();
	simData->SetMarkerStyle(2); //filled dot
	simData->SetMarkerColor(kBlue);
	graphMaximum = TMath::Max(graphMaximum, simData->GetMaximum(fragmentName));
	graphMaximum = TMath::Max(graphMaximum, ntuple->GetMaximum("y"));
	dummyHisto->SetMaximum(graphMaximum + .05*graphMaximum);
	dummyHisto->Draw();
	
	simData->Draw(fragmentName + ":depth","","p,same");

	ntuple->SetMarkerStyle(22); //triangle
    ntuple->SetMarkerColor(kRed);
	ntuple->Draw("y:x","","p,same");
	c1->SaveAs("fragmentYieldsFor" + fragmentName + ".png");
}