void FindAbundances(){


//Create MC histograms

Double_t NetOut_MC,type_MC;

Int_t Np_MC=0;
Int_t NFe_MC=0;
Int_t Nev;


TH1F *protonMC= new TH1F("protonMC","",100,-0.5,1.5);
TH1F *ironMC= new TH1F("ironMC","",100,-0.5,1.5);


TTree *DataMC= new TTree("Data","");
DataMC->ReadFile("MCSample.txt","NetOut_MC/D:type_MC/D"); 
Nev= DataMC->GetEntries();

DataMC->SetBranchAddress("NetOut_MC",&NetOut_MC); 
DataMC->SetBranchAddress("type_MC",&type_MC); 

for(Int_t i=0;i<Nev;i++){
	
    DataMC->GetEntry(i);   
	if(type_MC==0) protonMC->Fill(NetOut_MC);
	if(type_MC==1) ironMC->Fill(NetOut_MC);

}    

//Normalize MC histograms
protonMC_norm= protonMC->DrawNormalized();
ironMC_norm= ironMC->DrawNormalized("same");


//Create Data Histogram
Double_t NetOut_data,type_data;
Float_t Np_data=0;
Float_t NFe_data=0;
Int_t Ntot;

TH1F *DataSample= new TH1F("DataSample","",100,-0.5,1.5);

TTree *Data= new TTree("Data","");
Data->ReadFile("DataSample.txt","NetOut_data/D:type_data/D"); 
Ntot= Data->GetEntries();

Data->SetBranchAddress("NetOut_data",&NetOut_data); 
Data->SetBranchAddress("type_data",&type_data); 

for(Int_t i=0;i<Ntot;i++){
	Data->GetEntry(i);   
	DataSample->Fill(NetOut_data);
	if(type_data==0) Np_data++;
	if(type_data==1) NFe_data++;
}   
	
DataSample_norm= DataSample->DrawNormalized("same");
		
//Get the TRUE fractions
Double_t TrueFract_proton= Np_data/Ntot;
Double_t TrueFract_iron= NFe_data/Ntot;


//Fit the DataSample_norm with protonMC_norm and ironMC_norm
Double_t RecFract_proton, RecFractErr_proton;
Double_t RecFract_iron, RecFractErr_iron;

TObjArray *mc = new TObjArray(2);
mc->Add(protonMC_norm);
mc->Add(ironMC_norm);

TFractionFitter* fit = new TFractionFitter(DataSample_norm, mc); // initialise
fit->Constrain(1,0.,1.);               // constrain fraction 1 to be between 0 and 1
fit->Constrain(2,0.,1.); 
//fit->SetRangeX(1,15);                    // use only the first 15 bins in the fit
Int_t status = fit->Fit();               // perform the fit
cout << "fit status: " << status << endl;

if (status == 0) {                       // check on fit status
	TH1F* result = (TH1F*) fit->GetPlot();
	DataSample_norm->Draw();
	result->Draw("same");
	fit->GetResult(0,RecFract_proton,RecFractErr_proton);
    fit->GetResult(1,RecFract_iron,RecFractErr_iron);
	cout<<""<<endl;
	cout<<"***TRUE FRACTIONS***"<<endl;
	cout<<"True proton fraction "<<TrueFract_proton<<endl;
	cout<<"True iron fraction "<<TrueFract_iron<<endl;
	cout<<"***RECONSTRUCTED FRACTIONS***"<<endl;
	cout<<"Rec proton fraction "<<RecFract_proton<<" +- "<<RecFractErr_proton<<endl;
	cout<<"Rec iron fraction "<<RecFract_iron<<" +- "<<RecFractErr_iron<<endl;

}




}