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

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 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();
}

void fragmentEnergyDistributionDifferentAngles() {

gStyle->SetOptStat(0000000000); //remove the for this graphs totally redundant statbox

//   gROOT->SetStyle("clearRetro");

   TString pDepth, fragment, Znum, normToOneAtZeroAngle;
   cout << "Enter phantom depth (eg. 27.9, see experimentalData directory for choices): ";
   cin >> pDepth;  
   TString simulationDataPath = "IAEA_" + pDepth + ".root";
   TString dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
   dir.ReplaceAll("basic.C","");
   dir.ReplaceAll("/./","/");
   ifstream in;
   in.open(Form("experimentalData/iaeaBenchmark/fragmentEnergySpctra279mmWater0deg.dat",dir.Data()));
   Float_t f1,f2,f3, f4,f5,f6;
   Int_t nlines = 0;
   TFile *f = new TFile("fragmentEnergyWithAngularDistribution.root","RECREATE");
   TNtuple *ntuple = new TNtuple("ntuple","Data from ascii file","Energy:He:B:H:Li:Be");
	  
   Char_t DATAFLAG[4];
   Int_t NDATA;
   Char_t n1[6], n2[2], n3[2], n4[2], n5[2], n6[2];
   in >> DATAFLAG >> NDATA ; // Read EXFOR line: 'DATA 6'
   in >> n1 >> n2 >> n3 >> n4 >> n5 >> n6; // Read  column titles: 'Energy He B [...]'

   cout <<n1<<" "<<n2<<" "<<n3<<"    "<<n4<<"    "<<n5<<"   "<<n6<<"\n";
   while (1) {
      in >> f1 >> f2 >> f3 >>f4 >> f5 >> f6;
      if (!in.good()) break;
      if (nlines < 500 ) printf("%f  %0.2f %0.2f %0.2f %0.2f %0.2f \n",f1,f2,f3,f4,f5,f6);
      ntuple->Fill(f1,f2,f3,f4,f5,f6);
      nlines++;
   }
   

   //Let's pull in the simulation-data
   TFile *MCData = TFile::Open("IAEA_200000.root");
   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.
   
	//good to keep for ref. G4 might give weird units due to change.
	metadata->Scan();
	std::cout << "Recieved metadata-row: " << events << " " <<  detectorDistance << " " << waterThickness << " " << beamEnergy << " " << energyError << " " << phantomCenterDistance;


//A lot of hardcoded histograms, ugly
   Double_t binAmount = 50.0; //casting from int failed somehow, so in float temporarily, fixme
   Double_t maxEnergy = 450.0;
   Double_t binWidth = maxEnergy / binAmount;
   TH1F *hist1 = new TH1F("hist1", "", binAmount, 0.0, maxEnergy);
   TH1F *hist2 = new TH1F("hist2", "", binAmount, 0.0, maxEnergy);   
   TH1F *hist3 = new TH1F("hist3", "", binAmount, 0.0, maxEnergy);
   TH1F *hist4 = new TH1F("hist4", "", binAmount, 0.0, maxEnergy);
   TH1F *hist5 = new TH1F("hist5", "", binAmount, 0.0, maxEnergy);
   TH1F *hist6 = new TH1F("hist6", "", binAmount, 0.0, maxEnergy);
   TH1F *hist7 = new TH1F("hist7", "", binAmount, 0.0, maxEnergy);
   TH1F *hist8 = new TH1F("hist8", "", binAmount, 0.0, maxEnergy);
   TH1F *hist9 = new TH1F("hist9", "", binAmount, 0.0, maxEnergy);
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);
//.........这里部分代码省略.........

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+");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
   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");
}