StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MyDecay.cxx
1 #include <Riostream.h>
2 #include <TCanvas.h>
3 #include <TRandom.h>
4 #include <TGenPhaseSpace.h>
5 #include "TDatabasePDG.h"
6 #include <TParticle.h>
7 #include <TDecayChannel.h>
8 #include <assert.h>
9 
10 #include "MyDecay.h"
11 using namespace std;
12 
13 ClassImp(MyDecay)
14 
15 MyDecay::MyDecay(const char *file)
16 {
17  mFileName=TString(file);
18  isPP05=kFALSE;
19  isDAU=kFALSE;
20  cout<<"construct"<<endl;
21 }
22 MyDecay::~MyDecay()
23 {
24  cout<<"destruct"<<endl;
25 }
26 
27 TObjArray *MyDecay::decay(TParticle *part)
28 {
29  Float_t frac=gRandom->Rndm();
30  TParticlePDG *part_pdg=part->GetPDG();
31  Int_t i_dec=0;
32  Double_t sum_br=0;
33  while(i_dec<part_pdg->NDecayChannels())
34  {
35  sum_br+=part_pdg->DecayChannel(i_dec)->BranchingRatio();
36  if(frac<sum_br) break;
37  i_dec++;
38  }
39  if(i_dec==part_pdg->NDecayChannels()) return 0;
40 
41  TDecayChannel *decay=part_pdg->DecayChannel(i_dec);
42  Int_t n_part=decay->NDaughters();
43 
44  Double_t *mass_arr=new Double_t[n_part];
45  for (Int_t i=0; i<n_part; i++)
46  {
47  mass_arr[i]=TDatabasePDG::Instance()->GetParticle(decay->DaughterPdgCode(i))->Mass();
48  }
49 
50  TGenPhaseSpace *gs=new TGenPhaseSpace();//del
51  Double_t mass=part_pdg->Mass();
52  TLorentzVector v;
53  part->Momentum(v);
54  v.SetE(sqrt(mass*mass+v.P()*v.P()));
55  gs->SetDecay(v,n_part,mass_arr);
56  delete [] mass_arr;
57  Float_t max_wt=gs->GetWtMax();
58  while(gRandom->Rndm()*max_wt>gs->Generate());
59 
60  TObjArray *dec_part=new TObjArray(n_part);
61  for(Int_t i_part=0;i_part<n_part;i_part++)
62  {
63  TParticle *dght=new TParticle(decay->DaughterPdgCode(i_part),0,0,0,0,0,0,0,0,0,0,0,0,0);
64  dght->SetMomentum(*gs->GetDecay(i_part));
65  dec_part->Add(dght);
66  }
67  gs->Delete();
68 
69  return dec_part;
70 }
71 
72 void MyDecay::doDecay(Int_t ndecays)
73 {
74  if(!isPP05&&!isDAU){
75  cout<<"label not set"<<endl;
76  return;
77  }
78  //limits:
79  Double_t ptmin=.0;
80  Double_t ptmax=30.;
81  Double_t ptminbin=.0;
82  Double_t ptmaxbin=15.;
83  Int_t ptbins=60;
84  Double_t phimin=-TMath::Pi();
85  Double_t phimax=TMath::Pi();
86  Double_t etamin=-1.;
87  Double_t etamax=2.;
88 
89  Double_t rapmin=0.0;
90  Double_t rapmax=1.0;
91 
92  //create pi0:
93  TParticle *part=new TParticle(111,0,0,0,0,0,0,0,0,0,0,0,0,0);
94  TLorentzVector *v=new TLorentzVector();
95  //create eta:
96  TParticle *part2=new TParticle(221,0,0,0,0,0,0,0,0,0,0,0,0,0);
97  TLorentzVector *v2=new TLorentzVector();
98  //create omega:
99  TParticle *part3=new TParticle(223,0,0,0,0,0,0,0,0,0,0,0,0,0);
100  TLorentzVector *v3=new TLorentzVector();
101 
102  //output hists:
103  TH1F *pion=new TH1F("pion","# versus pT",ptbins,ptminbin,ptmaxbin);
104  TH1F *etameson=new TH1F("etameson","# versus pT",ptbins,ptminbin,ptmaxbin);
105  TH1F *omega=new TH1F("omega","# versus pT",ptbins,ptminbin,ptmaxbin);
106  TH1F *gamma=new TH1F("gamma","# versus pT",ptbins,ptminbin,ptmaxbin);
107  TH1F *gamma2=new TH1F("gamma2","# versus pT",ptbins,ptminbin,ptmaxbin);
108  TH1F *gamma3=new TH1F("gamma3","# versus pT",ptbins,ptminbin,ptmaxbin);
109 
110  for(Int_t i=0;i<ndecays;i++)
111  {
112  if(i%10000==0) cout<<i<<endl;
113 
114  //set part:
115  Float_t pT=ptmin+gRandom->Rndm()*(ptmax-ptmin);
116  Float_t phi=phimin+gRandom->Rndm()*(phimax-phimin);
117  Float_t eta=etamin+gRandom->Rndm()*(etamax-etamin);
118  Float_t mass=part->GetMass();
119  Float_t mass2=part2->GetMass();
120  Float_t mass3=part3->GetMass();
121  v->SetPtEtaPhiM(pT,eta,phi,mass);
122  v2->SetPtEtaPhiM(pT,eta,phi,mass2);
123  v3->SetPtEtaPhiM(pT,eta,phi,mass3);
124  part->SetMomentum(*v);
125  part2->SetMomentum(*v2);
126  part3->SetMomentum(*v3);
127 
128  Float_t WEIGHT=ptweight(pT)*phiweight(phi)*etaweight(eta);
129  Float_t WEIGHT2=ptweightETA(pT)*phiweight(phi)*etaweight(eta);
130  Float_t WEIGHT3=ptweightOMEGA(pT)*phiweight(phi)*etaweight(eta);
131 
132  if(v->PseudoRapidity()>rapmin && v->PseudoRapidity()<rapmax) pion->Fill(v->Pt(),WEIGHT);
133  if(v2->PseudoRapidity()>rapmin && v2->PseudoRapidity()<rapmax) etameson->Fill(v2->Pt(),WEIGHT2);
134  if(v3->PseudoRapidity()>rapmin && v3->PseudoRapidity()<rapmax) omega->Fill(v3->Pt(),WEIGHT3);
135 
136  //gen. decay:
137  TObjArray *dec_dght=decay(part);
138  TObjArray *dec_dght2=decay(part2);
139  TObjArray *dec_dght3=decay(part3);
140 
141  //get pion daughters:
142  Int_t nD=dec_dght->GetEntries();
143  for(Int_t j=0;j<nD;j++)
144  {
145  if(nD==2)
146  {
147  TParticle *d=(TParticle*)dec_dght->At(j);
148  if(d->GetPdgCode()==22 && d->Eta()>rapmin && d->Eta()<rapmax)
149  {
150  gamma->Fill(d->Pt(),WEIGHT);
151  }
152  }
153  else if(nD==3)
154  {
155  TParticle *d=(TParticle*)dec_dght->At(j);
156  if(d->GetPdgCode()==22 && d->Eta()>rapmin && d->Eta()<rapmax)
157  {
158  gamma->Fill(d->Pt(),WEIGHT);
159  }
160  }
161  }
162  //delete array:
163  dec_dght->Delete();
164  delete dec_dght;
165 
166  //get eta daughters:
167  Int_t nD2=dec_dght2->GetEntries();
168  for(Int_t j=0;j<nD2;j++)
169  {
170  if(nD2==2)
171  {
172  TParticle *d=(TParticle*)dec_dght2->At(j);
173  if(d->GetPdgCode()==22 && d->Eta()>rapmin && d->Eta()<rapmax)
174  {
175  gamma2->Fill(d->Pt(),WEIGHT2);
176  }
177  }
178  else if(nD2==3)
179  {
180  TParticle *d=(TParticle*)dec_dght2->At(j);
181  if(d->GetPdgCode()==22 && d->Eta()>rapmin && d->Eta()<rapmax)
182  {
183  gamma2->Fill(d->Pt(),WEIGHT2);
184  }
185  }
186  //there is e+e-pi+pi- for eta!!
187  }
188  //delete array:
189  dec_dght2->Delete();
190  delete dec_dght2;
191 
192  //get omega daughters:
193  Int_t nD3=dec_dght3->GetEntries();
194  for(Int_t j=0;j<nD3;j++)
195  {
196  if(nD3==2)
197  {
198  TParticle *d=(TParticle*)dec_dght3->At(j);
199  if(d->GetPdgCode()==22 && d->Eta()>rapmin && d->Eta()<rapmax)
200  {
201  gamma3->Fill(d->Pt(),WEIGHT3);
202  }
203  }
204  else if(nD3==3)
205  {
206  TParticle *d=(TParticle*)dec_dght3->At(j);
207  if(d->GetPdgCode()==22 && d->Eta()>rapmin && d->Eta()<rapmax)
208  {
209  gamma3->Fill(d->Pt(),WEIGHT3);
210  }
211  }
212  }
213  //delete array:
214  dec_dght3->Delete();
215  delete dec_dght3;
216 
217  }
218 
219 
220  TCanvas *c=new TCanvas("c","c",600,900);
221  c->Divide(2,3);
222  c->cd(1);
223  gPad->SetLogy();
224  gamma->SetTitle("pion(daughters) vs pT");
225  gamma->Draw();
226  gamma->SetLineStyle(2);
227  pion->Draw("same");
228  c->cd(2);
229  gPad->SetLogy();
230  gamma2->SetTitle("eta(daughters) vs pT");
231  gamma2->Draw();
232  gamma2->SetLineStyle(2);
233  etameson->Draw("same");
234  c->cd(3);
235  gPad->SetLogy();
236  gamma3->SetTitle("omega782(daughters) vs pT");
237  gamma3->Draw();
238  gamma3->SetLineStyle(2);
239  omega->Draw("same");
240  c->cd(4);
241  gPad->SetLogy();
242  pion->SetTitle("pions/etas/omegas vs pT");
243  pion->Draw();
244  pion->SetLineColor(2);
245  etameson->Draw("same");
246  etameson->SetLineColor(4);
247  omega->Draw("same");
248  omega->SetLineColor(3);
249  c->cd(5);
250  gPad->SetLogy();
251  TH1F *f=new TH1F(*gamma);
252  f->SetLineStyle(1);
253  TH1F *f2=new TH1F(*gamma2);
254  f2->SetLineStyle(1);
255  TH1F *f3=new TH1F(*gamma3);
256  f3->SetLineStyle(1);
257  TH1F *fsum=new TH1F(*f);
258  fsum->Add(f2);
259  fsum->Add(f3);
260  TH1F *g=new TH1F(*pion);
261  f->Divide(g);
262  f2->Divide(g);
263  f3->Divide(g);
264  fsum->Divide(g);
265  fsum->Draw();
266  fsum->SetTitle("bg photons per pion");
267  f->Draw("same");
268  f->SetLineColor(2);
269  f2->Draw("same");
270  f2->SetLineColor(4);
271  f3->Draw("same");
272  f3->SetLineColor(3);
273  fsum->SetMinimum(.001);
274  c->cd(6);
275  TH1F *etatopi=new TH1F(*etameson);
276  TH1F *omegatopi=new TH1F(*omega);
277  etatopi->Divide(pion);
278  omegatopi->Divide(pion);
279  omegatopi->SetTitle("ratio of #eta/#omega to #pi");
280  omegatopi->Draw();
281  etatopi->Draw("same");
282  c->cd(0);
283  c->SaveAs((mFileName+".ps").Data());
284  c->SaveAs((mFileName+".root").Data());
285 
286  TFile *dec_out=new TFile((mFileName+"Sum.root").Data(),"recreate");
287  fsum->Write();
288  f->Write("gamma_pion");
289  f2->Write("gamma_eta");
290  f3->Write("gamma_omega");
291  dec_out->Close();
292 }
293 
294 Double_t MyDecay::ptweight(Double_t x)
295 {
296  float weight=1.;
297  if(x<0.2) x=0.2;
298  if(isPP05){
299  weight=300.*x*pow(1.+x,-9.0);//error: +/- 1.0%
300  }
301  if(isDAU){
302  weight=8.24601e-02*x*pow(1.+x,-9.01897e+00);//error: +/- 1.0%
303  }
304  return weight;
305 }
306 Double_t MyDecay::ptweightETA(Double_t x)
307 {
308  Float_t ME=0.548;//eta meson mass
309  Float_t MP=0.135;//pi0 mass
310  Float_t weight=fETAMTSCALE * x/sqrt(x*x + ME*ME - MP*MP) * ptweight(sqrt(x*x + ME*ME - MP*MP));
311  return weight;
312 }
313 Double_t MyDecay::ptweightOMEGA(Double_t x)
314 {
315  Float_t MO=0.782;//omega mass
316  Float_t MP=0.135;//pi0 mass
317  Float_t weight=1. * x/sqrt(x*x + MO*MO - MP*MP) * ptweight(sqrt(x*x + MO*MO - MP*MP));
318  return weight;
319 }
320 Double_t MyDecay::etaweight(Double_t /*x*/)
321 {
322  return 1.;
323 }
324 Double_t MyDecay::phiweight(Double_t /*x*/)
325 {
326  return 1.;
327 }
Generated by   doxygen 1.8.5