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+1];
  51   TH1D * hInvMassBgk[NumberOfPoints+1];
  52   TH1D * hInvMassSub[NumberOfPoints+1];
  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
  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;
  81     }
  82
  83
  84     if (ipoint!=NumberOfPoints) sprintf(histoname,"hInvMassAll_%2.0f",100*pT[ipoint]);
  85     else sprintf(histoname,"hInvMassAll_full");
  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
  91     if (ipoint!=NumberOfPoints)  sprintf(histoname,"hInvMassSub_%2.0f",100*pT[ipoint]);
  92     else sprintf(histoname,"hInvMassSub_full");
  93     hInvMassSub[ipoint] =  new TH1D(*nume);
  94     hInvMassSub[ipoint]->SetName(histoname); hInvMassSub[ipoint]->SetTitle(histoname);
  95     hInvMassSub[ipoint]->Sumw2();
  96
  97     if (ipoint!=NumberOfPoints)  sprintf(histoname,"hInvMassBgk_%2.0f",100*pT[ipoint]);
  98     else sprintf(histoname,"hInvMassBgk_full");
  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);
 106
 107     integral1 = hInvMassSub[ipoint]->FindBin( (hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(1)-2.*hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(2) ) );
 108     integral2 = hInvMassSub[ipoint]->FindBin( (hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(1)+2.*hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(2) ) );
 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++) {
 125         e_dNdpT[ipoint]+=hInvMassSub[ipoint]->GetBinError(ibin)*hInvMassSub[ipoint]->GetBinError(ibin);
 126       }
 127       e_dNdpT[ipoint] = TMath::Sqrt( e_dNdpT[ipoint]);
 128     }
 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
 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++) {
 233       e_NumberOfResonances+=hInvMassSub[NumberOfPoints]->GetBinError(ibin)*hInvMassSub[NumberOfPoints]->GetBinError(ibin);
 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
 243 }