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

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 jackknife( I_XMLConfigReader * xmlFile, MinimiserConfiguration * theMinimiser, 
		FitFunctionConfiguration * theFunction, ParameterSet* argumentParameterSet, vector<string> CommandLineParam, int start, int stop )
	{
		cout << "Starting JackKnife" << endl;
		PDFWithData  * pdfAndData  = xmlFile->GetPDFsAndData()[0];
		ParameterSet * parSet      = xmlFile->GetFitParameters( CommandLineParam );
		PhaseSpaceBoundary * phase = xmlFile->GetPhaseSpaceBoundaries()[0];
		pdfAndData->SetPhysicsParameters( parSet );
			
		MemoryDataSet * dataset = (MemoryDataSet*) pdfAndData->GetDataSet();
		vector<IDataSet*> data;
		data.push_back(dataset);
	
		IPDF * pdf = pdfAndData->GetPDF();
		vector<IPDF*> vectorPDFs;
		vectorPDFs.push_back(pdf);
	
		// Repeat nominal fit to get the nominal value of the physics parameter	
		FitResult * nominal = FitAssembler::DoFit( theMinimiser, theFunction, argumentParameterSet, vectorPDFs, data, xmlFile->GetConstraints() );
		double nominal_value = nominal->GetResultParameterSet()->GetResultParameter("tau")->GetValue();
		cout << nominal_value << endl;
		double jackknifed_value = 0.;

		MemoryDataSet * subset = new MemoryDataSet( phase );
		int nData = dataset->Yield();	
                string fileName = ResultFormatter::GetOutputFolder();
                fileName.append("/jackknife.root");         
		TFile * outputFile = new TFile("jackknife.root", "RECREATE");	
                TNtuple * jack = new TNtuple("jack", "jackknifed - nominal", "diff_jackknifed_nominal:reco_time:true_time");

		double reco_time = 0.;
		double true_time = 0.;

		for ( int i = start; i < stop; i++ )
		{
			for ( int j = 0; j < nData; j++ )
			{
				if ( j == i ) {
					reco_time = dataset->GetDataPoint( j )->GetObservable( "time" )->GetValue();
					true_time = dataset->GetDataPoint( j )->GetObservable( "truetime" )->GetValue();
					continue;
				}
				subset->AddDataPoint( dataset->GetDataPoint( j ) );
			}
			cout << "Creating subset " << i << " containing " << subset->Yield() << " candidates" << endl;
			vector<IDataSet*> vectorData;
			vectorData.push_back(subset);
			FitResult * result = FitAssembler::DoFit( theMinimiser, theFunction, argumentParameterSet, vectorPDFs, vectorData, xmlFile->GetConstraints() );
			if ( result->GetFitStatus() == 3 ) {
				jackknifed_value = result->GetResultParameterSet()->GetResultParameter("tau")->GetValue();
				jack->Fill( (Float_t)(jackknifed_value - nominal_value), (Float_t)reco_time, (Float_t)true_time);
			}
			subset->Clear();
			delete result;
		}
		outputFile->Write();
		outputFile->Close();
	}

int execQuery() 
{
   TNtuple *ed = new TNtuple("ed","ed","pt:val:stat");
   ed->SetScanField(0);
   ed->Fill(1.0, 5.0, 4.0);
   TSQLRow *row = ed->Query("pt:val:stat")->Next();
   fprintf(stdout,"1: %s 2: %s 3: %s\n",row->GetField(0),row->GetField(1),row->GetField(2));
   return 0;
}

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

  float E, p, px, py, pz;
  float phi,cos_theta;

  const float PI=3.14159;
  const float mpi=0.13957018;
  const float mrho=0.77549; 
  TLorentzVector P1, P2;

  TFile *f = new TFile("decayang.root","RECREATE");

  TNtuple *cmspi = new TNtuple("cmspi","Pi decay to two photons in CMS","px1:py1:pz1:E1:px2:py2:pz2:E2");
  TNtuple *cmsrho = new TNtuple("cmsrho","Rho decay to two pions in CMS","px1:py1:pz1:E1:px2:py2:pz2:E2");
  TNtuple *boostpi = new TNtuple("boostpi","Pi decay to two photons in lab frames","b1px1:b1py1:b1pz1:b1e1:b1px2:b1py2:b1pz2:b1e2");


  for (int i=0; i<10000; i++) {
    phi = 2.0 * PI *(gRandom->Rndm());                // generate phi
    cos_theta = 2.0 * (gRandom->Rndm())-1;            //generate cos(theta)

     //Start with the pi to photon decays
    E=mpi/2.0;
    p=E;
    px=p*cos(phi)*sqrt(1.0-pow(cos_theta,2));
    py=p*sin(phi)*sqrt(1.0-pow(cos_theta,2));
    pz=p*cos_theta;
    cmspi->Fill(px, py, pz, E, -px, -py, -pz, E);
    
    

     //Now do the rho to pi decays
    E=mrho/2.0;
    p=sqrt(E*E-mpi*mpi);
    px = p * cos(phi) * sqrt( 1.0 - pow(cos_theta,2));
    py = p * sin(phi) * sqrt( 1.0 - pow(cos_theta,2));
    pz = p * cos_theta;
    cmsrho->Fill(px, py, pz, E, -px, -py, -pz, E);
  }

  f->Write(); 
}

		void visit(const SmoothParticle &p) {
			size_t overlaps;
			//float directRho = dmf->getRho(particles[i].position, overlaps);
			//float r = particles[i].position.distanceTo(ComaPositionKpc);
#if GADGET_ROOT_ENABLED
				ntuple->Fill(r, particles[i].rho, directRho, overlaps);
#endif
			//out << r << " " << particles[i].rho << " " << directRho << " "
			//		<< overlaps << std::endl;
		}

	void writeVolumeWeightedProfile() {
		std::cout << "** Write Volume Weighted Profile" << std::endl;

		const size_t steps = 50;
		const float rMin = 5;
		const float rMax = ComaRadiusKpc;
		const size_t nPoints = 10000;

#if GADGET_ROOT_ENABLED
		TFile *file = new TFile("coma_profile_volume_weighted.root", "RECREATE",
				"Rho Weighted");
		if (file->IsZombie())
		throw std::runtime_error(
				"Root output file cannot be properly opened");
		TNtuple *ntuple = new TNtuple("events", "Rho", "r:B");
#endif
		std::ofstream out("coma_profile_volume_weighted.csv");
		out << "r B" << std::endl;

		float stepLog = (log10(rMax) - log10(rMin)) / (steps - 1);
		for (size_t i = 0; i < steps; i++) {
			std::cout << ".";
			std::cout.flush();

			float r = pow(10, log10(rMin) + stepLog * i);

			float avgB = 0;
			for (size_t i = 0; i < nPoints; i++) {
				Vector3f p = randomUnitVector();
				p *= r;
				p += ComaPositionKpc;
				Vector3f b;
				bool isGood = dmf->getField(p, b);
				avgB += b.length(); // / nPoints;
			}
			avgB /= nPoints;

#if GADGET_ROOT_ENABLED
			ntuple->Fill(r, avgB);
#endif
			out << r << " " << avgB << std::endl;

		}

#if GADGET_ROOT_ENABLED
		file->Write();
		file->Close();
#endif
		std::cout << std::endl;
	}

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 CheckDALIDet(char * ridffile="test.ridf"){

  gSystem->Load("libanacore.so");

  TArtEventStore *estore = new TArtEventStore();
  estore->Open(ridffile);
  TArtRawEventObject *rawevent = estore->GetRawEventObject();

  TFile *fout = new TFile("dataid.root","RECREATE");
  TNtuple *ntp = new TNtuple("ntp","ntp","geo:ch:val");

  int neve = 0;
  while(estore->GetNextEvent() && neve<100){

    int eve_number=rawevent->GetEventNumber();
    for(int i=0;i<rawevent->GetNumSeg();i++){
      TArtRawSegmentObject *seg = rawevent->GetSegment(i);
      int device = seg->GetDevice();
      int fp = seg->GetFP();
      int detector = seg->GetDetector();
      int module = seg->GetModule();
      if(DALI==device&&DALIA==detector){
	cout << "    seg:"<< i <<" dev:"<< device 
	     << " fp:"<<fp<< " det:"<<detector<< " mod:"<<module
	     << " #data=" << seg->GetNumData() << endl;
	for(int j=0;j<seg->GetNumData();j++){
	  TArtRawDataObject *d = seg->GetData(j);
	  int geo = d->GetGeo();
	  int ch = d->GetCh();
	  unsigned int val = d->GetVal();
	  cout << "       geo:" << geo 
	       << " ch:" << ch << " val:" << val << endl;
	  ntp->Fill((float)geo,(float)ch,(float)val);
	}
      }
    }

    rawevent->Clear();
    neve ++;
  }
  fout->Write();
  fout->Close();

}

void RIDF2DataId(){

  gSystem->Load("libanacore.so");

  TArtEventStore *estore = new TArtEventStore();
  estore->Open("test.ridf");
  TArtRawEventObject *rawevent = estore->GetRawEventObject();

  TFile *fout = new TFile("dataid.root","RECREATE");
  TNtuple *ntp = new TNtuple("ntp","ntp","sys:cat:det:dtype:val");

  int neve = 0;
  while(estore->GetNextEvent() && neve < 10000){

    for(int i=0;i<rawevent->GetNumSeg();i++){
      TArtRawSegmentObject *seg = rawevent->GetSegment(i);
      int device = seg->GetDevice();
      int fp = seg->GetFP();
      int detector = seg->GetDetector();
      int module = seg->GetModule();

      for(int j=0;j<seg->GetNumData();j++){
	TArtRawDataObject *d = seg->GetData(j);
	int geo = d->GetGeo(); 
	int ch = d->GetCh();
	int val = d->GetVal(); 
	int cat = d->GetCategoryID(); 
	int det = d->GetDetectorID(); 
	int id = d->GetDatatypeID(); 

	ntp->Fill((float)detector,(float)cat,(float)det,(float)id,(float)val);
      }
    }

    estore->ClearData();
    neve ++;
  }

  fout->Write();
  fout->Close();

}

void ReadAsciiCreateRootFile(){

  //  Read data from an ascii file and create a root file with an histogram and an ntuple.
//   see a variant of this macro in basic2.C
  //Author: Rene Brun


  // read file basic.data
  //this file has 3 columns of float data

  TString dir = gSystem->UnixPathName(__FILE__);
  dir.ReplaceAll("ReadAsciiCreateRootFile.C","");
  dir.ReplaceAll("/./","/");
  ifstream in;
  in.open(Form("%sbasic.dat", dir.Data()));

  Float_t x,y,z;
  Int_t nlines = 0;
  TFile *f = new TFile("basic.root","RECREATE");
  TH1F *h1= new TH1F("h1","x distribution", 100, -4,4);
  TNtuple *ntuple = new TNtuple("ntuple","data from ascii file","x:y:z");

  while(1){

    in >> x >> y >> z; //format of the data
    if (!in.good()) break;
    if(nlines<5) printf("x=%8f, y=%8f,z=%8f \n", x,y,z);

    h1->Fill(x);
    ntuple->Fill(x,y,z  );
    nlines++;

  }

  printf("found %d points \n", nlines);

  in.close();
  f->Write();


}

TNtuple* LikelihoodSpace::get_contour(float delta) {
  TNtuple* contour = (TNtuple*) this->samples->Clone("lscontour");
  contour->Reset();

  // Get list of branch names
  std::vector<std::string> names;
  for (int i=0; i<this->samples->GetListOfBranches()->GetEntries(); i++) {
    std::string name = this->samples->GetListOfBranches()->At(i)->GetName();
    if (name == "likelihood") {
      continue;
    }
    names.push_back(name);
  }

  float* params_branch = new float[names.size()];
  for (size_t i=0; i<names.size(); i++) {
    this->samples->SetBranchAddress(names[i].c_str(), &params_branch[i]);
  }

  float ml_branch;
  this->samples->SetBranchAddress("likelihood", &ml_branch);

  // Build a new TNtuple with samples inside the contour
  float* v = new float[names.size() + 1];
  for (int i=0; i<this->samples->GetEntries(); i++) {
    this->samples->GetEntry(i);
    if (ml_branch < this->ml + delta) {
      for (size_t j=0; j<names.size(); j++) {
        v[j] = params_branch[j];
      }
      v[names.size()] = ml_branch;
      contour->Fill(v);
    }
  }

  this->samples->ResetBranchAddresses();
  delete[] v;

  return contour;
}

void read18() {
    ifstream in;
    in.open("/Users/Yuichi/root/macros/mytext18.txt");
    
    Float_t a,b,c,d,e,k,g,h,i,j;
    Int_t nlines = 0;
    TFile *f = new TFile("read18.root","RECREATE");
    TNtuple *ntuple = new TNtuple("ntuple","data from reading18 file","a:b:c:d:e:k:g:h:i:j");
    
    while (1) {
        in >> a >> b >> c >> d >> e >> k >> g >> h >> i >> j;
        if (!in.good()) break;
        if (nlines < 100) printf("a=%8f, b=%8f, c=%8f, d=%8f, e=%8f, k=%8f, g=%8f, h=%8f, i=%8f, j=%8f\n",a,b,c,d,e,k,g,h,i,j);
        ntuple->Fill(a,b,c,d,e,k,g,h,i,j);
        nlines++;
    }
    printf(" found %d points\n",nlines);
    ntuple->Scan("a:b:c:d:e:k:g:h:i:j");
    
    in.close();
    
    f->Write();
}

void testPicoD0EventRead(TString filename)
{
	gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
	loadSharedLibraries();

	gSystem->Load("StPicoDstMaker");
	gSystem->Load("StPicoD0Maker");

	TFile* f = new TFile(filename.Data());
	TTree* T = (TTree*)f->Get("T");
	StPicoD0Event* event = new StPicoD0Event();
	T->SetBranchAddress("dEvent",&event);

	TFile ff("read_test.root","RECREATE");
  TNtuple* nt = new TNtuple("nt","","m:pt:eta:phi:theta:"
                                     "decayL:kDca:pDca:dca12:cosThetaStar");

	StKaonPion* kp = 0;

	for(Int_t i=0;i<100000;++i)
	{
		T->GetEntry(i);

		TClonesArray* arrKPi = event->kaonPionArray();

		for(int idx=0;idx<event->nKaonPion();++idx)
		{
			kp = (StKaonPion*)arrKPi->At(idx);

      nt->Fill(kp->m(),kp->pt(),kp->eta(),kp->phi(),kp->pointingAngle(),
              kp->decayLength(),kp->kaonDca(),kp->pionDca(),kp->dcaDaughters(),kp->cosThetaStar());
		}
	}

  nt->Write();
	ff.Close();
}

void read21n() {
    ifstream in;
    in.open("/Users/Yuichi/root/macros/mytext21.txt");    //Check
    
    Float_t hm,zm;      //Values
    Int_t nlines = 0;
    TFile *f = new TFile("read21.root","RECREATE");    //Output name
    TNtuple *ntuple = new TNtuple("ntuple","data from read20_3 file","hm:zm");    //Values
    
    while (1) {
        in >> hm >> zm;   //input
        if (!in.good()) break;
        if (nlines < 100) printf("hm=%8f, zm=%8f\n",hm,zm);  //Values
        ntuple->Fill(hm,zm);       //Fill values
        nlines++;
    }
    printf(" found %d points\n",nlines);
    ntuple->Scan("hm:zm");  //Scan Values
    
    in.close();
    
    f->Write();
    f->Close();
}

void read25cn() {
    ifstream in;
    in.open("/Users/Yuichi/root/macros/mytext25c.txt");    //Check
    
    Float_t dRm,dRmbar,dRb,dRbbar;      //Values
    Int_t nlines = 0;
    TFile *f = new TFile("read25c.root","RECREATE");    //Output name
    TNtuple *ntuple = new TNtuple("ntuple","data from read25c file","dRm:dRmbar:dRb:dRbbar");    //Values
    
    while (1) {
        in >> dRm >> dRmbar >> dRb >> dRbbar;   //input
        if (!in.good()) break;
        if (nlines < 100) printf("dRm=%8f, dRmbar=%8f, dRb=%8f, dRbbar=%8f\n",dRm,dRmbar,dRb,dRbbar);  //Values
        ntuple->Fill(dRm,dRmbar,dRb,dRbbar);       //Fill values
        nlines++;
    }
    printf(" found %d points\n",nlines);
    ntuple->Scan("dRm:dRmbar:dRb:dRbbar");  //Scan Values
    
    in.close();
    
    f->Write();
    f->Close();
}