e506910949caac23a9860a13decbb6b30cffa385
[u/mrichter/AliRoot.git] / MUON / MUONQuarkoniaPtDistri.C
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
50   TH1D * hInvMassAll[NumberOfPoints];
51   TH1D * hInvMassBgk[NumberOfPoints];
52   TH1D * hInvMassSub[NumberOfPoints];
53
54   char histoname[60];
55
56   Int_t ibin, ibin1,ibin2,  ipoint;
57   Int_t integral1, integral2;
58   for(ipoint=0; ipoint<NumberOfPoints; ipoint++) {
59     ibin1 = ipoint * NRes;
60     ibin2 = ibin1+  NRes;
61     pT[ipoint]=0.;
62     e_pT[ipoint]=0.;
63     dNdpT[ipoint]=0.;
64     e_dNdpT[ipoint]=0.;
65     for(ibin=ibin1; ibin<ibin2; ibin++) {
66       pT[ipoint]+=hInvMassAll_vs_Pt->GetYaxis()->GetBinCenter(ibin);
67     }
68     pT[ipoint]/=NRes;
69     pT1[ipoint]=pT[ipoint];
70     pT2[ipoint]=pT[ipoint];
71
72     e_pT[ipoint] = (hInvMassAll_vs_Pt->GetYaxis()->GetBinCenter(ibin2)-hInvMassAll_vs_Pt->GetYaxis()->GetBinCenter(ibin1))/2.;
73     e_pT1[ipoint]=e_pT[ipoint];
74     e_pT2[ipoint]=e_pT[ipoint];
75
76     sprintf(histoname,"hInvMassAll_%2.0f",100*pT[ipoint]);
77     TH1D * nume = hInvMassAll_vs_Pt->ProjectionX("nume",ibin1,ibin2-1);
78     hInvMassAll[ipoint] =  new TH1D(*nume);
79     hInvMassAll[ipoint]->SetName(histoname); hInvMassAll[ipoint]->SetTitle(histoname);
80     hInvMassAll[ipoint]->Sumw2();
81
82     sprintf(histoname,"hInvMassSub_%2.0f",100*pT[ipoint]);
83     hInvMassSub[ipoint] =  new TH1D(*nume);
84     hInvMassSub[ipoint]->SetName(histoname); hInvMassSub[ipoint]->SetTitle(histoname);
85     hInvMassSub[ipoint]->Sumw2();
86
87     sprintf(histoname,"hInvMassBgk_%2.0f",100*pT[ipoint]);
88     TH1D * deno = hInvMassBgk_vs_Pt->ProjectionX("deno",ibin1,ibin2-1);
89     hInvMassBgk[ipoint]=  new TH1D(*deno);
90     hInvMassBgk[ipoint]->SetName(histoname);hInvMassBgk[ipoint]->SetTitle(histoname);
91     hInvMassBgk[ipoint]->Sumw2();
92
93     hInvMassSub[ipoint]->Add(hInvMassBgk[ipoint],-1.);
94     hInvMassSub[ipoint]->Fit("gaus","","",InvMass1,InvMass2);
95     Cent[ipoint] = hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(1);
96     e_Cent[ipoint] = hInvMassSub[ipoint]->GetFunction("gaus")->GetParError(1);
97     Sigma[ipoint] = hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(2);
98     e_Sigma[ipoint] = hInvMassSub[ipoint]->GetFunction("gaus")->GetParError(2);
99     integral1 = hInvMassSub[ipoint]->FindBin( (hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(1)-2.*hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(2) ) );
100     integral2 = hInvMassSub[ipoint]->FindBin( (hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(1)+2.*hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(2) ) );
101     dNdpT[ipoint] = hInvMassSub[ipoint]->Integral(integral1,integral2);
102     for(ibin=integral1; ibin<=integral2; ibin++) {
103       e_dNdpT[ipoint]+=hInvMassSub[ipoint]->GetBinContent(ibin)*hInvMassSub[ipoint]->GetBinContent(ibin);
104     }
105     e_dNdpT[ipoint] = TMath::Sqrt( e_dNdpT[ipoint]);
106   }
107
108   MUONmassPlot.Close();
109     
110   TFile out("QuarkoniaPtDistri.root","RECREATE");
111   TGraphErrors * gePtDistri = new TGraphErrors(NumberOfPoints,pT,dNdpT,e_pT, e_dNdpT);
112   TGraphErrors * geCenter   = new TGraphErrors(NumberOfPoints,pT1,Cent,e_pT1, e_Cent);
113   TGraphErrors * geSigma   = new TGraphErrors(NumberOfPoints,pT2,Sigma,e_pT2, e_Sigma);
114   gePtDistri->Write();
115   geCenter->Write();
116   geSigma->Write();
117   for(ipoint=0;ipoint<NumberOfPoints; ipoint++) {
118     hInvMassAll[ipoint]->Write();
119     hInvMassBgk[ipoint]->Write();
120     hInvMassSub[ipoint]->Write();
121   }
122   out.Write();
123   out.Close();
124
125   TCanvas *c1 =new TCanvas("c1","c1");
126   TH2D * hframe1 = new TH2D("hframe1","hframe1",10,0,12,10,InvMass1,InvMass2);
127   hframe1->SetXTitle("Transverse Momentum (GeV/c)");
128   hframe1->SetYTitle("Resonance Centroide (GeV/C^{2})");
129   hframe1->Draw();
130   geCenter->SetMarkerStyle(27);
131   geCenter->Draw("p");
132   
133   TCanvas *c2 =new TCanvas("c2","c2");
134   TH2D * hframe2 = new TH2D("hframe2","hframe2",10,0,12,10,0.03,0.200); 
135   hframe2->SetXTitle("Transverse Momentum (GeV/c)");
136   hframe2->SetYTitle("Resonance Sigma (GeV/C^2)");
137   hframe2->Draw();
138   geSigma->SetMarkerStyle(28);
139   geSigma->Draw("p");
140   
141   TCanvas *c3 =new TCanvas("c3","c3");
142   TH2D * hframe3 = new TH2D("hframe3","hframe3",10,0,12,10,0.1,10000);
143   hframe3->SetXTitle("Transverse Momentum (GeV/c)");
144   hframe3->SetYTitle("dN/dpT (GeV^{-1})");
145   hframe3->Draw();
146   c3->SetLogy();
147   hframe3->Draw();
148   gePtDistri->SetMarkerStyle(29);
149   gePtDistri->Draw("p");
150
151   TCanvas *c4 =new TCanvas("c4","c4");
152   c4->SetLogy();
153   c4->Divide(3,3);
154   Int_t i; 
155   Float_t izoom1;
156   Float_t izoom2; 
157   for (i=1;i<10;i++) {
158     izoom1 = hInvMassSub[i]->FindBin( (hInvMassSub[i]->GetFunction("gaus")->GetParameter(1)-15.*hInvMassSub[i]->GetFunction("gaus")->GetParameter(2) ) );
159     izoom2 = hInvMassSub[i]->FindBin( (hInvMassSub[i]->GetFunction("gaus")->GetParameter(1)+15.*hInvMassSub[i]->GetFunction("gaus")->GetParameter(2) ) );
160     c4->cd(i);
161     hInvMassAll[i]->SetXTitle("Invariant Mass (GeV/C^{2})");
162     hInvMassAll[i]->SetYTitle("Counts (1/25 MeV)");
163     hInvMassAll[i]->GetXaxis()->SetRange(izoom1,izoom2);
164     hInvMassAll[i]->Draw();
165     hInvMassBgk[i]->SetLineColor(4);
166     hInvMassBgk[i]->Draw("histo,same");
167   }
168  
169   TCanvas *c5 =new TCanvas("c5","c5"); 
170   c5->Divide(3,3);
171   Int_t i;
172   for (i=1;i<10;i++) {
173     c5->cd(i);
174     izoom1 = hInvMassSub[i]->FindBin( (hInvMassSub[i]->GetFunction("gaus")->GetParameter(1)-15.*hInvMassSub[i]->GetFunction("gaus")->GetParameter(2) ) );
175     izoom2 = hInvMassSub[i]->FindBin( (hInvMassSub[i]->GetFunction("gaus")->GetParameter(1)+15.*hInvMassSub[i]->GetFunction("gaus")->GetParameter(2) ) );
176     hInvMassSub[i]->GetXaxis()->SetRange(izoom1,izoom2);
177     hInvMassSub[i]->SetXTitle("Invariant Mass (GeV/C^{2})");
178     hInvMassSub[i]->SetYTitle("Counts (1/25 MeV)");
179     hInvMassSub[i]->Draw();
180   }
181
182 }