[u/mrichter/AliRoot.git] / MUON / MUONQuarkoniaPtDistri.C

//
// Macro for analysis of 2d InvMass v Pt distribution for the generation of:
// - Pt distribution of quarkonia, 
// - Mass and sigma of the quarkonia resonance as a function of pT
// - and invariant mass distributions ( all, background and subtracted)
//   as a function of pT
// Usual utilisation:
// QuarkoniaPtDistri("ProdMuonCocktailPer1/MUONmassPlot.root",4,2.9,3.3)
//
//  MUONmassPlotFile  root file form MUONmassPlot_ESD.C
//  NRes=1 bining of the pT distribution 
//    ==1, equal to MUONmassPlotFile 2D histo, 
//    ==2, bin twice larger
//  InvMass1, InvMass2 Range of the quarkonia resonance
//
// Gines MARTINEZ, Subatech, May 04
//

void MUONQuarkoniaPtDistri(char * MUONmassPlotFile="MUONmassPlot.root", Int_t NRes=4, Float_t InvMass1=2.9, Float_t InvMass2=3.3)
{
  TFile MUONmassPlot(MUONmassPlotFile);

  TH2F * hInvMassAll_vs_Pt = (TH2F*) MUONmassPlot.Get("hInvMassAll_vs_Pt");
  TH2F * hInvMassBgk_vs_Pt = (TH2F*) MUONmassPlot.Get("hInvMassBgk_vs_Pt");
  hInvMassAll_vs_Pt->Sumw2();
  hInvMassBgk_vs_Pt->Sumw2();

  gROOT->cd();

  Int_t   Nbins = hInvMassAll_vs_Pt->GetYaxis()->GetNbins();
  Float_t ptMin = hInvMassAll_vs_Pt->GetYaxis()->GetXmin();
  Float_t ptMax = hInvMassAll_vs_Pt->GetYaxis()->GetXmax();

  const Int_t NumberOfPoints = Nbins/NRes;
  Float_t *pT = new Float_t[NumberOfPoints];
  Float_t *e_pT = new Float_t[NumberOfPoints];
  Float_t *dNdpT = new Float_t[NumberOfPoints];
  Float_t *e_dNdpT = new Float_t[NumberOfPoints];

  Float_t *pT1 = new Float_t[NumberOfPoints];
  Float_t *e_pT1 = new Float_t[NumberOfPoints];
  Float_t *Cent = new Float_t[NumberOfPoints];
  Float_t *e_Cent = new Float_t[NumberOfPoints];

  Float_t *pT2 = new Float_t[NumberOfPoints];
  Float_t *e_pT2 = new Float_t[NumberOfPoints];
  Float_t *Sigma = new Float_t[NumberOfPoints];
  Float_t *e_Sigma = new Float_t[NumberOfPoints];

  TH1D * hInvMassAll[NumberOfPoints+1];
  TH1D * hInvMassBgk[NumberOfPoints+1];
  TH1D * hInvMassSub[NumberOfPoints+1];

  char histoname[60];

  Int_t ibin, ibin1,ibin2,  ipoint;
  Int_t integral1, integral2;
  for(ipoint=0; ipoint<=NumberOfPoints; ipoint++) {

    if (ipoint!=NumberOfPoints) {
      ibin1 = ipoint * NRes;
      ibin2 = ibin1+  NRes;
      pT[ipoint]=0.;
      e_pT[ipoint]=0.;
      dNdpT[ipoint]=0.;
      e_dNdpT[ipoint]=0.;
      for(ibin=ibin1; ibin<ibin2; ibin++) {
	pT[ipoint]+=hInvMassAll_vs_Pt->GetYaxis()->GetBinCenter(ibin);
      }
      pT[ipoint]/=NRes;
      pT1[ipoint]=pT[ipoint];
      pT2[ipoint]=pT[ipoint];
      
      e_pT[ipoint] = (hInvMassAll_vs_Pt->GetYaxis()->GetBinCenter(ibin2)-hInvMassAll_vs_Pt->GetYaxis()->GetBinCenter(ibin1))/2.;
      e_pT1[ipoint]=e_pT[ipoint];
      e_pT2[ipoint]=e_pT[ipoint];
    } 
    else {
      ibin1 = 0;
      ibin2 = Nbins;
    }


    if (ipoint!=NumberOfPoints) sprintf(histoname,"hInvMassAll_%2.0f",100*pT[ipoint]);
    else sprintf(histoname,"hInvMassAll_full");
    TH1D * nume = hInvMassAll_vs_Pt->ProjectionX("nume",ibin1,ibin2-1);
    hInvMassAll[ipoint] =  new TH1D(*nume);
    hInvMassAll[ipoint]->SetName(histoname); hInvMassAll[ipoint]->SetTitle(histoname);
    hInvMassAll[ipoint]->Sumw2();

    if (ipoint!=NumberOfPoints)  sprintf(histoname,"hInvMassSub_%2.0f",100*pT[ipoint]);
    else sprintf(histoname,"hInvMassSub_full");
    hInvMassSub[ipoint] =  new TH1D(*nume);
    hInvMassSub[ipoint]->SetName(histoname); hInvMassSub[ipoint]->SetTitle(histoname);
    hInvMassSub[ipoint]->Sumw2();

    if (ipoint!=NumberOfPoints)  sprintf(histoname,"hInvMassBgk_%2.0f",100*pT[ipoint]);
    else sprintf(histoname,"hInvMassBgk_full");
    TH1D * deno = hInvMassBgk_vs_Pt->ProjectionX("deno",ibin1,ibin2-1);
    hInvMassBgk[ipoint]=  new TH1D(*deno);
    hInvMassBgk[ipoint]->SetName(histoname);hInvMassBgk[ipoint]->SetTitle(histoname);
    hInvMassBgk[ipoint]->Sumw2();

    hInvMassSub[ipoint]->Add(hInvMassBgk[ipoint],-1.);
    hInvMassSub[ipoint]->Fit("gaus","","",InvMass1,InvMass2);

    integral1 = hInvMassSub[ipoint]->FindBin( (hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(1)-2.*hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(2) ) );
    integral2 = hInvMassSub[ipoint]->FindBin( (hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(1)+2.*hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(2) ) );

    // Renormalisation to take into account the resonance statistics
    if ( hInvMassBgk[ipoint]->Integral()>0.) {
      Float_t renormalisation = (hInvMassBgk[ipoint]->Integral() - hInvMassSub[ipoint]->Integral(integral1,integral2))/hInvMassBgk[ipoint]->Integral();
      hInvMassBgk[ipoint]->Scale(renormalisation);
      hInvMassSub[ipoint]->Reset("ICE");
      hInvMassSub[ipoint]->Add(hInvMassAll[ipoint],+1.);
      hInvMassSub[ipoint]->Add(hInvMassBgk[ipoint],-1.);
    }
    if (ipoint!=NumberOfPoints) { 
      Cent[ipoint] = hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(1);
      e_Cent[ipoint] = hInvMassSub[ipoint]->GetFunction("gaus")->GetParError(1);
      Sigma[ipoint] = hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(2);
      e_Sigma[ipoint] = hInvMassSub[ipoint]->GetFunction("gaus")->GetParError(2);
      dNdpT[ipoint] = hInvMassSub[ipoint]->Integral(integral1,integral2);
      for(ibin=integral1; ibin<=integral2; ibin++) {
	e_dNdpT[ipoint]+=hInvMassSub[ipoint]->GetBinError(ibin)*hInvMassSub[ipoint]->GetBinError(ibin);
      }
      e_dNdpT[ipoint] = TMath::Sqrt( e_dNdpT[ipoint]);
    }
  }

  MUONmassPlot.Close();
    
  TFile out("QuarkoniaPtDistri.root","RECREATE");
  TGraphErrors * gePtDistri = new TGraphErrors(NumberOfPoints,pT,dNdpT,e_pT, e_dNdpT);
  TGraphErrors * geCenter   = new TGraphErrors(NumberOfPoints,pT1,Cent,e_pT1, e_Cent);
  TGraphErrors * geSigma   = new TGraphErrors(NumberOfPoints,pT2,Sigma,e_pT2, e_Sigma);
  gePtDistri->Write();
  geCenter->Write();
  geSigma->Write();
  for(ipoint=0;ipoint<NumberOfPoints; ipoint++) {
    hInvMassAll[ipoint]->Write();
    hInvMassBgk[ipoint]->Write();
    hInvMassSub[ipoint]->Write();
  }
  out.Write();
  out.Close();

  TCanvas *c1 =new TCanvas("c1","c1");
  TH2D * hframe1 = new TH2D("hframe1","hframe1",10,0,12,10,InvMass1,InvMass2);
  hframe1->SetXTitle("Transverse Momentum (GeV/c)");
  hframe1->SetYTitle("Resonance Centroide (GeV/C^{2})");
  hframe1->Draw();
  geCenter->SetMarkerStyle(27);
  geCenter->Draw("p");
  
  TCanvas *c2 =new TCanvas("c2","c2");
  TH2D * hframe2 = new TH2D("hframe2","hframe2",10,0,12,10,0.03,0.200); 
  hframe2->SetXTitle("Transverse Momentum (GeV/c)");
  hframe2->SetYTitle("Resonance Sigma (GeV/C^2)");
  hframe2->Draw();
  geSigma->SetMarkerStyle(28);
  geSigma->Draw("p");
  
  TCanvas *c3 =new TCanvas("c3","c3");
  TH2D * hframe3 = new TH2D("hframe3","hframe3",10,0,12,10,0.1,10000);
  hframe3->SetXTitle("Transverse Momentum (GeV/c)");
  hframe3->SetYTitle("dN/dpT (GeV^{-1})");
  hframe3->Draw();
  c3->SetLogy();
  hframe3->Draw();
  gePtDistri->SetMarkerStyle(29);
  gePtDistri->Draw("p");

  TCanvas *c4 =new TCanvas("c4","c4");
  c4->SetLogy();
  c4->Divide(3,3);
  Int_t i; 
  Float_t izoom1;
  Float_t izoom2; 
  for (i=1;i<10;i++) {
    izoom1 = hInvMassSub[i]->FindBin( (hInvMassSub[i]->GetFunction("gaus")->GetParameter(1)-15.*hInvMassSub[i]->GetFunction("gaus")->GetParameter(2) ) );
    izoom2 = hInvMassSub[i]->FindBin( (hInvMassSub[i]->GetFunction("gaus")->GetParameter(1)+15.*hInvMassSub[i]->GetFunction("gaus")->GetParameter(2) ) );
    c4->cd(i);
    hInvMassAll[i]->SetXTitle("Invariant Mass (GeV/C^{2})");
    hInvMassAll[i]->SetYTitle("Counts (1/25 MeV)");
    hInvMassAll[i]->GetXaxis()->SetRange(izoom1,izoom2);
    hInvMassAll[i]->Draw();
    hInvMassBgk[i]->SetLineColor(4);
    hInvMassBgk[i]->Draw("histo,same");
  }
 
  TCanvas *c5 =new TCanvas("c5","c5"); 
  c5->Divide(3,3);
  Int_t i;
  for (i=1;i<10;i++) {
    c5->cd(i);
    izoom1 = hInvMassSub[i]->FindBin( (hInvMassSub[i]->GetFunction("gaus")->GetParameter(1)-15.*hInvMassSub[i]->GetFunction("gaus")->GetParameter(2) ) );
    izoom2 = hInvMassSub[i]->FindBin( (hInvMassSub[i]->GetFunction("gaus")->GetParameter(1)+15.*hInvMassSub[i]->GetFunction("gaus")->GetParameter(2) ) );
    hInvMassSub[i]->GetXaxis()->SetRange(izoom1,izoom2);
    hInvMassSub[i]->SetXTitle("Invariant Mass (GeV/C^{2})");
    hInvMassSub[i]->SetYTitle("Counts (1/25 MeV)");
    hInvMassSub[i]->Draw();
  }

  TCanvas *c6 =new TCanvas("c6","c6"); 
  c6->Divide(2,1);
  c6->cd(1);
  izoom1 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)-15.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
  izoom2 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)+15.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
  c4->cd(i);
  hInvMassAll[NumberOfPoints]->SetXTitle("Invariant Mass (GeV/C^{2})");
  hInvMassAll[NumberOfPoints]->SetYTitle("Counts (1/25 MeV)");
  hInvMassAll[NumberOfPoints]->GetXaxis()->SetRange(izoom1,izoom2);
  hInvMassAll[NumberOfPoints]->Draw();
  hInvMassBgk[NumberOfPoints]->SetLineColor(4);
  hInvMassBgk[NumberOfPoints]->Draw("histo,same");
  c6->cd(2);
  izoom1 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)-15.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
  izoom2 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)+15.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
  hInvMassSub[NumberOfPoints]->GetXaxis()->SetRange(izoom1,izoom2);
  hInvMassSub[NumberOfPoints]->SetXTitle("Invariant Mass (GeV/C^{2})");
  hInvMassSub[NumberOfPoints]->SetYTitle("Counts (1/25 MeV)");
  hInvMassSub[NumberOfPoints]->Draw(); 

  //Number of resonances
  
    integral1 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)-2.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
    integral2 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)+2.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
    
    Int_t   NumberOfResonances = hInvMassSub[NumberOfPoints]->Integral(integral1,integral2);
    Int_t e_NumberOfResonances = 0;
    for(ibin=integral1; ibin<=integral2; ibin++) {
      e_NumberOfResonances+=hInvMassSub[NumberOfPoints]->GetBinError(ibin)*hInvMassSub[NumberOfPoints]->GetBinError(ibin);
    }
    e_NumberOfResonances= TMath::Sqrt( e_NumberOfResonances);
    printf(">>> Number of resonances is %d+-%d in the inv mass range (2sigma) %f, %f \n",NumberOfResonances,e_NumberOfResonances,
	   hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)-2.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ,
	   hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)+2.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) );;

    printf(">>> Centroide is %f +-%f GeV\n", hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1), hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParError(1));
    printf(">>> Sigma is %f +-%f MeV \n", 1000.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2), 1000.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParError(2));

}
Commit	Line	Data
4692d494	1	//
	2	// Macro for analysis of 2d InvMass v Pt distribution for the generation of:
	3	// - Pt distribution of quarkonia,
	4	// - Mass and sigma of the quarkonia resonance as a function of pT
	5	// - and invariant mass distributions ( all, background and subtracted)
	6	// as a function of pT
	7	// Usual utilisation:
	8	// QuarkoniaPtDistri("ProdMuonCocktailPer1/MUONmassPlot.root",4,2.9,3.3)
	9	//
	10	// MUONmassPlotFile root file form MUONmassPlot_ESD.C
	11	// NRes=1 bining of the pT distribution
	12	// ==1, equal to MUONmassPlotFile 2D histo,
	13	// ==2, bin twice larger
	14	// InvMass1, InvMass2 Range of the quarkonia resonance
	15	//
	16	// Gines MARTINEZ, Subatech, May 04
	17	//
	18
	19	void MUONQuarkoniaPtDistri(char * MUONmassPlotFile="MUONmassPlot.root", Int_t NRes=4, Float_t InvMass1=2.9, Float_t InvMass2=3.3)
	20	{
	21	TFile MUONmassPlot(MUONmassPlotFile);
	22
	23	TH2F * hInvMassAll_vs_Pt = (TH2F*) MUONmassPlot.Get("hInvMassAll_vs_Pt");
	24	TH2F * hInvMassBgk_vs_Pt = (TH2F*) MUONmassPlot.Get("hInvMassBgk_vs_Pt");
	25	hInvMassAll_vs_Pt->Sumw2();
	26	hInvMassBgk_vs_Pt->Sumw2();
	27
	28	gROOT->cd();
	29
	30	Int_t Nbins = hInvMassAll_vs_Pt->GetYaxis()->GetNbins();
	31	Float_t ptMin = hInvMassAll_vs_Pt->GetYaxis()->GetXmin();
	32	Float_t ptMax = hInvMassAll_vs_Pt->GetYaxis()->GetXmax();
	33
	34	const Int_t NumberOfPoints = Nbins/NRes;
	35	Float_t *pT = new Float_t[NumberOfPoints];
	36	Float_t *e_pT = new Float_t[NumberOfPoints];
	37	Float_t *dNdpT = new Float_t[NumberOfPoints];
	38	Float_t *e_dNdpT = new Float_t[NumberOfPoints];
	39
	40	Float_t *pT1 = new Float_t[NumberOfPoints];
	41	Float_t *e_pT1 = new Float_t[NumberOfPoints];
	42	Float_t *Cent = new Float_t[NumberOfPoints];
	43	Float_t *e_Cent = new Float_t[NumberOfPoints];
	44
	45	Float_t *pT2 = new Float_t[NumberOfPoints];
	46	Float_t *e_pT2 = new Float_t[NumberOfPoints];
	47	Float_t *Sigma = new Float_t[NumberOfPoints];
	48	Float_t *e_Sigma = new Float_t[NumberOfPoints];
	49
642ca698	50	TH1D * hInvMassAll[NumberOfPoints+1];
	51	TH1D * hInvMassBgk[NumberOfPoints+1];
	52	TH1D * hInvMassSub[NumberOfPoints+1];
4692d494	53
	54	char histoname[60];
	55
	56	Int_t ibin, ibin1,ibin2, ipoint;
	57	Int_t integral1, integral2;
642ca698	58	for(ipoint=0; ipoint<=NumberOfPoints; ipoint++) {
	59
	60	if (ipoint!=NumberOfPoints) {
	61	ibin1 = ipoint * NRes;
	62	ibin2 = ibin1+ NRes;
	63	pT[ipoint]=0.;
	64	e_pT[ipoint]=0.;
	65	dNdpT[ipoint]=0.;
	66	e_dNdpT[ipoint]=0.;
	67	for(ibin=ibin1; ibin<ibin2; ibin++) {
	68	pT[ipoint]+=hInvMassAll_vs_Pt->GetYaxis()->GetBinCenter(ibin);
	69	}
	70	pT[ipoint]/=NRes;
	71	pT1[ipoint]=pT[ipoint];
	72	pT2[ipoint]=pT[ipoint];
	73
	74	e_pT[ipoint] = (hInvMassAll_vs_Pt->GetYaxis()->GetBinCenter(ibin2)-hInvMassAll_vs_Pt->GetYaxis()->GetBinCenter(ibin1))/2.;
	75	e_pT1[ipoint]=e_pT[ipoint];
	76	e_pT2[ipoint]=e_pT[ipoint];
	77	}
	78	else {
	79	ibin1 = 0;
	80	ibin2 = Nbins;
4692d494	81	}
4692d494	82
4692d494	83
642ca698	84	if (ipoint!=NumberOfPoints) sprintf(histoname,"hInvMassAll_%2.0f",100*pT[ipoint]);
642ca698	85	else sprintf(histoname,"hInvMassAll_full");
4692d494	86	TH1D * nume = hInvMassAll_vs_Pt->ProjectionX("nume",ibin1,ibin2-1);
	87	hInvMassAll[ipoint] = new TH1D(*nume);
	88	hInvMassAll[ipoint]->SetName(histoname); hInvMassAll[ipoint]->SetTitle(histoname);
	89	hInvMassAll[ipoint]->Sumw2();
	90
642ca698	91	if (ipoint!=NumberOfPoints) sprintf(histoname,"hInvMassSub_%2.0f",100*pT[ipoint]);
642ca698	92	else sprintf(histoname,"hInvMassSub_full");
4692d494	93	hInvMassSub[ipoint] = new TH1D(*nume);
	94	hInvMassSub[ipoint]->SetName(histoname); hInvMassSub[ipoint]->SetTitle(histoname);
	95	hInvMassSub[ipoint]->Sumw2();
	96
642ca698	97	if (ipoint!=NumberOfPoints) sprintf(histoname,"hInvMassBgk_%2.0f",100*pT[ipoint]);
642ca698	98	else sprintf(histoname,"hInvMassBgk_full");
4692d494	99	TH1D * deno = hInvMassBgk_vs_Pt->ProjectionX("deno",ibin1,ibin2-1);
	100	hInvMassBgk[ipoint]= new TH1D(*deno);
	101	hInvMassBgk[ipoint]->SetName(histoname);hInvMassBgk[ipoint]->SetTitle(histoname);
	102	hInvMassBgk[ipoint]->Sumw2();
	103
	104	hInvMassSub[ipoint]->Add(hInvMassBgk[ipoint],-1.);
	105	hInvMassSub[ipoint]->Fit("gaus","","",InvMass1,InvMass2);
642ca698	106
4692d494	107	integral1 = hInvMassSub[ipoint]->FindBin( (hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(1)-2.*hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(2) ) );
4692d494	108	integral2 = hInvMassSub[ipoint]->FindBin( (hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(1)+2.*hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(2) ) );
642ca698	109
	110	// Renormalisation to take into account the resonance statistics
	111	if ( hInvMassBgk[ipoint]->Integral()>0.) {
	112	Float_t renormalisation = (hInvMassBgk[ipoint]->Integral() - hInvMassSub[ipoint]->Integral(integral1,integral2))/hInvMassBgk[ipoint]->Integral();
	113	hInvMassBgk[ipoint]->Scale(renormalisation);
	114	hInvMassSub[ipoint]->Reset("ICE");
	115	hInvMassSub[ipoint]->Add(hInvMassAll[ipoint],+1.);
	116	hInvMassSub[ipoint]->Add(hInvMassBgk[ipoint],-1.);
	117	}
	118	if (ipoint!=NumberOfPoints) {
	119	Cent[ipoint] = hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(1);
	120	e_Cent[ipoint] = hInvMassSub[ipoint]->GetFunction("gaus")->GetParError(1);
	121	Sigma[ipoint] = hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(2);
	122	e_Sigma[ipoint] = hInvMassSub[ipoint]->GetFunction("gaus")->GetParError(2);
	123	dNdpT[ipoint] = hInvMassSub[ipoint]->Integral(integral1,integral2);
	124	for(ibin=integral1; ibin<=integral2; ibin++) {
a4103cd4	125	e_dNdpT[ipoint]+=hInvMassSub[ipoint]->GetBinError(ibin)*hInvMassSub[ipoint]->GetBinError(ibin);
642ca698	126	}
642ca698	127	e_dNdpT[ipoint] = TMath::Sqrt( e_dNdpT[ipoint]);
4692d494	128	}
4692d494	129	}
	130
	131	MUONmassPlot.Close();
	132
	133	TFile out("QuarkoniaPtDistri.root","RECREATE");
	134	TGraphErrors * gePtDistri = new TGraphErrors(NumberOfPoints,pT,dNdpT,e_pT, e_dNdpT);
	135	TGraphErrors * geCenter = new TGraphErrors(NumberOfPoints,pT1,Cent,e_pT1, e_Cent);
	136	TGraphErrors * geSigma = new TGraphErrors(NumberOfPoints,pT2,Sigma,e_pT2, e_Sigma);
	137	gePtDistri->Write();
	138	geCenter->Write();
	139	geSigma->Write();
	140	for(ipoint=0;ipoint<NumberOfPoints; ipoint++) {
	141	hInvMassAll[ipoint]->Write();
	142	hInvMassBgk[ipoint]->Write();
	143	hInvMassSub[ipoint]->Write();
	144	}
	145	out.Write();
	146	out.Close();
	147
	148	TCanvas *c1 =new TCanvas("c1","c1");
	149	TH2D * hframe1 = new TH2D("hframe1","hframe1",10,0,12,10,InvMass1,InvMass2);
	150	hframe1->SetXTitle("Transverse Momentum (GeV/c)");
	151	hframe1->SetYTitle("Resonance Centroide (GeV/C^{2})");
	152	hframe1->Draw();
	153	geCenter->SetMarkerStyle(27);
	154	geCenter->Draw("p");
	155
	156	TCanvas *c2 =new TCanvas("c2","c2");
	157	TH2D * hframe2 = new TH2D("hframe2","hframe2",10,0,12,10,0.03,0.200);
	158	hframe2->SetXTitle("Transverse Momentum (GeV/c)");
	159	hframe2->SetYTitle("Resonance Sigma (GeV/C^2)");
	160	hframe2->Draw();
	161	geSigma->SetMarkerStyle(28);
	162	geSigma->Draw("p");
	163
	164	TCanvas *c3 =new TCanvas("c3","c3");
	165	TH2D * hframe3 = new TH2D("hframe3","hframe3",10,0,12,10,0.1,10000);
	166	hframe3->SetXTitle("Transverse Momentum (GeV/c)");
	167	hframe3->SetYTitle("dN/dpT (GeV^{-1})");
	168	hframe3->Draw();
	169	c3->SetLogy();
	170	hframe3->Draw();
	171	gePtDistri->SetMarkerStyle(29);
	172	gePtDistri->Draw("p");
	173
	174	TCanvas *c4 =new TCanvas("c4","c4");
	175	c4->SetLogy();
	176	c4->Divide(3,3);
	177	Int_t i;
	178	Float_t izoom1;
	179	Float_t izoom2;
	180	for (i=1;i<10;i++) {
	181	izoom1 = hInvMassSub[i]->FindBin( (hInvMassSub[i]->GetFunction("gaus")->GetParameter(1)-15.*hInvMassSub[i]->GetFunction("gaus")->GetParameter(2) ) );
	182	izoom2 = hInvMassSub[i]->FindBin( (hInvMassSub[i]->GetFunction("gaus")->GetParameter(1)+15.*hInvMassSub[i]->GetFunction("gaus")->GetParameter(2) ) );
	183	c4->cd(i);
	184	hInvMassAll[i]->SetXTitle("Invariant Mass (GeV/C^{2})");
	185	hInvMassAll[i]->SetYTitle("Counts (1/25 MeV)");
	186	hInvMassAll[i]->GetXaxis()->SetRange(izoom1,izoom2);
	187	hInvMassAll[i]->Draw();
	188	hInvMassBgk[i]->SetLineColor(4);
	189	hInvMassBgk[i]->Draw("histo,same");
	190	}
	191
	192	TCanvas *c5 =new TCanvas("c5","c5");
193	c5->Divide(3,3);
194	Int_t i;
195	for (i=1;i<10;i++) {
196	c5->cd(i);
197	izoom1 = hInvMassSub[i]->FindBin( (hInvMassSub[i]->GetFunction("gaus")->GetParameter(1)-15.*hInvMassSub[i]->GetFunction("gaus")->GetParameter(2) ) );
198	izoom2 = hInvMassSub[i]->FindBin( (hInvMassSub[i]->GetFunction("gaus")->GetParameter(1)+15.*hInvMassSub[i]->GetFunction("gaus")->GetParameter(2) ) );
199	hInvMassSub[i]->GetXaxis()->SetRange(izoom1,izoom2);
200	hInvMassSub[i]->SetXTitle("Invariant Mass (GeV/C^{2})");
201	hInvMassSub[i]->SetYTitle("Counts (1/25 MeV)");
202	hInvMassSub[i]->Draw();
203	}
204
642ca698	205	TCanvas *c6 =new TCanvas("c6","c6");
	206	c6->Divide(2,1);
	207	c6->cd(1);
	208	izoom1 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)-15.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
	209	izoom2 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)+15.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
	210	c4->cd(i);
	211	hInvMassAll[NumberOfPoints]->SetXTitle("Invariant Mass (GeV/C^{2})");
	212	hInvMassAll[NumberOfPoints]->SetYTitle("Counts (1/25 MeV)");
	213	hInvMassAll[NumberOfPoints]->GetXaxis()->SetRange(izoom1,izoom2);
	214	hInvMassAll[NumberOfPoints]->Draw();
	215	hInvMassBgk[NumberOfPoints]->SetLineColor(4);
	216	hInvMassBgk[NumberOfPoints]->Draw("histo,same");
	217	c6->cd(2);
	218	izoom1 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)-15.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
	219	izoom2 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)+15.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
	220	hInvMassSub[NumberOfPoints]->GetXaxis()->SetRange(izoom1,izoom2);
	221	hInvMassSub[NumberOfPoints]->SetXTitle("Invariant Mass (GeV/C^{2})");
	222	hInvMassSub[NumberOfPoints]->SetYTitle("Counts (1/25 MeV)");
	223	hInvMassSub[NumberOfPoints]->Draw();
	224
	225	//Number of resonances
	226
	227	integral1 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)-2.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
	228	integral2 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)+2.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
	229
	230	Int_t NumberOfResonances = hInvMassSub[NumberOfPoints]->Integral(integral1,integral2);
	231	Int_t e_NumberOfResonances = 0;
	232	for(ibin=integral1; ibin<=integral2; ibin++) {
a4103cd4	233	e_NumberOfResonances+=hInvMassSub[NumberOfPoints]->GetBinError(ibin)*hInvMassSub[NumberOfPoints]->GetBinError(ibin);
642ca698	234	}
	235	e_NumberOfResonances= TMath::Sqrt( e_NumberOfResonances);
	236	printf(">>> Number of resonances is %d+-%d in the inv mass range (2sigma) %f, %f \n",NumberOfResonances,e_NumberOfResonances,
	237	hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)-2.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ,
	238	hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)+2.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) );;
	239
	240	printf(">>> Centroide is %f +-%f GeV\n", hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1), hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParError(1));
	241	printf(">>> Sigma is %f +-%f MeV \n", 1000.hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2), 1000.hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParError(2));
	242
4692d494	243	}