Macro for generation of quarkonia pT distribution (Gines)

[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];
  TH1D * hInvMassBgk[NumberOfPoints];
  TH1D * hInvMassSub[NumberOfPoints];

  char histoname[60];

  Int_t ibin, ibin1,ibin2,  ipoint;
  Int_t integral1, integral2;
  for(ipoint=0; ipoint<NumberOfPoints; ipoint++) {
    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];

    sprintf(histoname,"hInvMassAll_%2.0f",100*pT[ipoint]);
    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();

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

    sprintf(histoname,"hInvMassBgk_%2.0f",100*pT[ipoint]);
    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);
    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);
    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) ) );
    dNdpT[ipoint] = hInvMassSub[ipoint]->Integral(integral1,integral2);
    for(ibin=integral1; ibin<=integral2; ibin++) {
      e_dNdpT[ipoint]+=hInvMassSub[ipoint]->GetBinContent(ibin)*hInvMassSub[ipoint]->GetBinContent(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();
  }

}