StRoot: StEmcPool/StPhotonAnalysis/MyDecay.cxx Source File

00001 #include <Riostream.h>
00002 #include <TCanvas.h>
00003 #include <TRandom.h>
00004 #include <TGenPhaseSpace.h>
00005 #include "TDatabasePDG.h"
00006 #include <TParticle.h>
00007 #include <TDecayChannel.h>
00008 #include <assert.h>
00009 
00010 #include "MyDecay.h"
00011 
00012 ClassImp(MyDecay)
00013 
00014 MyDecay::MyDecay(const char *file)
00015 {
00016   mFileName=TString(file);
00017   isPP05=kFALSE;
00018   isDAU=kFALSE;
00019   cout<<"construct"<<endl;
00020 }
00021 MyDecay::~MyDecay()
00022 {
00023   cout<<"destruct"<<endl;
00024 }
00025 
00026 TObjArray *MyDecay::decay(TParticle *part) 
00027 {
00028   Float_t frac=gRandom->Rndm();
00029   TParticlePDG *part_pdg=part->GetPDG();
00030   Int_t i_dec=0;
00031   Double_t sum_br=0;
00032   while(i_dec<part_pdg->NDecayChannels()) 
00033     {
00034       sum_br+=part_pdg->DecayChannel(i_dec)->BranchingRatio();
00035       if(frac<sum_br) break;
00036       i_dec++;
00037     }
00038   if(i_dec==part_pdg->NDecayChannels()) return 0;
00039 
00040   TDecayChannel *decay=part_pdg->DecayChannel(i_dec);
00041   Int_t n_part=decay->NDaughters();
00042 
00043   Double_t *mass_arr=new Double_t[n_part];
00044   for (Int_t i=0; i<n_part; i++) 
00045     {
00046       mass_arr[i]=TDatabasePDG::Instance()->GetParticle(decay->DaughterPdgCode(i))->Mass();
00047     }
00048 
00049   TGenPhaseSpace *gs=new TGenPhaseSpace();//del
00050   Double_t mass=part_pdg->Mass();
00051   TLorentzVector v;
00052   part->Momentum(v);
00053   v.SetE(sqrt(mass*mass+v.P()*v.P()));
00054   gs->SetDecay(v,n_part,mass_arr);
00055   delete [] mass_arr;
00056   Float_t max_wt=gs->GetWtMax();
00057   while(gRandom->Rndm()*max_wt>gs->Generate());
00058 
00059   TObjArray *dec_part=new TObjArray(n_part);
00060   for(Int_t i_part=0;i_part<n_part;i_part++)
00061     {
00062       TParticle *dght=new TParticle(decay->DaughterPdgCode(i_part),0,0,0,0,0,0,0,0,0,0,0,0,0);
00063       dght->SetMomentum(*gs->GetDecay(i_part));
00064       dec_part->Add(dght);      
00065     }
00066   gs->Delete();
00067   
00068   return dec_part;
00069 }
00070 
00071 void MyDecay::doDecay(Int_t ndecays)
00072 {
00073   if(!isPP05&&!isDAU){
00074     cout<<"label not set"<<endl;
00075     return;
00076   }
00077   //limits:
00078   Double_t ptmin=.0;
00079   Double_t ptmax=30.;
00080   Double_t ptminbin=.0;
00081   Double_t ptmaxbin=15.;
00082   Int_t ptbins=60;
00083   Double_t phimin=-TMath::Pi();
00084   Double_t phimax=TMath::Pi();
00085   Double_t etamin=-1.;
00086   Double_t etamax=2.;
00087 
00088   Double_t rapmin=0.0;
00089   Double_t rapmax=1.0;
00090 
00091   //create pi0:
00092   TParticle *part=new TParticle(111,0,0,0,0,0,0,0,0,0,0,0,0,0);
00093   TLorentzVector *v=new TLorentzVector();
00094   //create eta:
00095   TParticle *part2=new TParticle(221,0,0,0,0,0,0,0,0,0,0,0,0,0);
00096   TLorentzVector *v2=new TLorentzVector();
00097   //create omega:
00098   TParticle *part3=new TParticle(223,0,0,0,0,0,0,0,0,0,0,0,0,0);
00099   TLorentzVector *v3=new TLorentzVector();
00100 
00101   //output hists:
00102   TH1F *pion=new TH1F("pion","# versus pT",ptbins,ptminbin,ptmaxbin);
00103   TH1F *etameson=new TH1F("etameson","# versus pT",ptbins,ptminbin,ptmaxbin);
00104   TH1F *omega=new TH1F("omega","# versus pT",ptbins,ptminbin,ptmaxbin);
00105   TH1F *gamma=new TH1F("gamma","# versus pT",ptbins,ptminbin,ptmaxbin);
00106   TH1F *gamma2=new TH1F("gamma2","# versus pT",ptbins,ptminbin,ptmaxbin);
00107   TH1F *gamma3=new TH1F("gamma3","# versus pT",ptbins,ptminbin,ptmaxbin);
00108 
00109   for(Int_t i=0;i<ndecays;i++) 
00110     {
00111       if(i%10000==0) cout<<i<<endl;
00112 
00113       //set part:
00114       Float_t pT=ptmin+gRandom->Rndm()*(ptmax-ptmin);
00115       Float_t phi=phimin+gRandom->Rndm()*(phimax-phimin);
00116       Float_t eta=etamin+gRandom->Rndm()*(etamax-etamin);
00117       Float_t mass=part->GetMass();
00118       Float_t mass2=part2->GetMass();
00119       Float_t mass3=part3->GetMass();
00120       v->SetPtEtaPhiM(pT,eta,phi,mass);
00121       v2->SetPtEtaPhiM(pT,eta,phi,mass2);
00122       v3->SetPtEtaPhiM(pT,eta,phi,mass3);
00123       part->SetMomentum(*v);
00124       part2->SetMomentum(*v2);
00125       part3->SetMomentum(*v3);
00126 
00127       Float_t WEIGHT=ptweight(pT)*phiweight(phi)*etaweight(eta);
00128       Float_t WEIGHT2=ptweightETA(pT)*phiweight(phi)*etaweight(eta);      
00129       Float_t WEIGHT3=ptweightOMEGA(pT)*phiweight(phi)*etaweight(eta);
00130 
00131       if(v->PseudoRapidity()>rapmin && v->PseudoRapidity()<rapmax) pion->Fill(v->Pt(),WEIGHT);
00132       if(v2->PseudoRapidity()>rapmin && v2->PseudoRapidity()<rapmax) etameson->Fill(v2->Pt(),WEIGHT2);
00133       if(v3->PseudoRapidity()>rapmin && v3->PseudoRapidity()<rapmax) omega->Fill(v3->Pt(),WEIGHT3);
00134 
00135       //gen. decay:
00136       TObjArray *dec_dght=decay(part);
00137       TObjArray *dec_dght2=decay(part2);
00138       TObjArray *dec_dght3=decay(part3);
00139       
00140       //get pion daughters:
00141       Int_t nD=dec_dght->GetEntries();
00142       for(Int_t j=0;j<nD;j++)
00143         {
00144           if(nD==2)
00145             {
00146               TParticle *d=(TParticle*)dec_dght->At(j);
00147               if(d->GetPdgCode()==22 && d->Eta()>rapmin &&  d->Eta()<rapmax) 
00148                 {
00149                   gamma->Fill(d->Pt(),WEIGHT);
00150                 }
00151             }
00152           else if(nD==3)
00153             {
00154               TParticle *d=(TParticle*)dec_dght->At(j);
00155               if(d->GetPdgCode()==22 && d->Eta()>rapmin &&  d->Eta()<rapmax) 
00156                 {
00157                   gamma->Fill(d->Pt(),WEIGHT);
00158                 }
00159             }
00160         }      
00161       //delete array:
00162       dec_dght->Delete();
00163       delete dec_dght;
00164 
00165       //get eta daughters:
00166       Int_t nD2=dec_dght2->GetEntries();
00167       for(Int_t j=0;j<nD2;j++)
00168         {
00169           if(nD2==2)
00170             {
00171               TParticle *d=(TParticle*)dec_dght2->At(j);
00172               if(d->GetPdgCode()==22 && d->Eta()>rapmin &&  d->Eta()<rapmax)
00173                 {
00174                   gamma2->Fill(d->Pt(),WEIGHT2);
00175                 }
00176             }
00177           else if(nD2==3)
00178             {
00179               TParticle *d=(TParticle*)dec_dght2->At(j);
00180               if(d->GetPdgCode()==22 && d->Eta()>rapmin &&  d->Eta()<rapmax)
00181                 {
00182                   gamma2->Fill(d->Pt(),WEIGHT2);
00183                 }
00184             }
00185           //there is e+e-pi+pi- for eta!!
00186         }
00187       //delete array:
00188       dec_dght2->Delete();
00189       delete dec_dght2;
00190 
00191       //get omega daughters:
00192       Int_t nD3=dec_dght3->GetEntries();
00193       for(Int_t j=0;j<nD3;j++)
00194         {
00195           if(nD3==2)
00196             {
00197               TParticle *d=(TParticle*)dec_dght3->At(j);
00198               if(d->GetPdgCode()==22 && d->Eta()>rapmin &&  d->Eta()<rapmax)
00199                 {
00200                   gamma3->Fill(d->Pt(),WEIGHT3);
00201                 }
00202             }
00203           else if(nD3==3)
00204             {
00205               TParticle *d=(TParticle*)dec_dght3->At(j);
00206               if(d->GetPdgCode()==22 && d->Eta()>rapmin &&  d->Eta()<rapmax)
00207                 {
00208                   gamma3->Fill(d->Pt(),WEIGHT3);
00209                 }
00210             }
00211         }
00212       //delete array:
00213       dec_dght3->Delete();
00214       delete dec_dght3;
00215 
00216     }
00217 
00218 
00219   TCanvas *c=new TCanvas("c","c",600,900);
00220   c->Divide(2,3);
00221   c->cd(1);
00222   gPad->SetLogy();
00223   gamma->SetTitle("pion(daughters) vs pT");
00224   gamma->Draw();
00225   gamma->SetLineStyle(2);
00226   pion->Draw("same");
00227   c->cd(2);
00228   gPad->SetLogy();
00229   gamma2->SetTitle("eta(daughters) vs pT");
00230   gamma2->Draw();
00231   gamma2->SetLineStyle(2);
00232   etameson->Draw("same");
00233   c->cd(3);
00234   gPad->SetLogy();
00235   gamma3->SetTitle("omega782(daughters) vs pT");
00236   gamma3->Draw();
00237   gamma3->SetLineStyle(2);
00238   omega->Draw("same");
00239   c->cd(4);
00240   gPad->SetLogy();
00241   pion->SetTitle("pions/etas/omegas vs pT");
00242   pion->Draw();
00243   pion->SetLineColor(2);
00244   etameson->Draw("same");
00245   etameson->SetLineColor(4);
00246   omega->Draw("same");
00247   omega->SetLineColor(3);
00248   c->cd(5);
00249   gPad->SetLogy();
00250   TH1F *f=new TH1F(*gamma);
00251   f->SetLineStyle(1);
00252   TH1F *f2=new TH1F(*gamma2);
00253   f2->SetLineStyle(1);
00254   TH1F *f3=new TH1F(*gamma3);
00255   f3->SetLineStyle(1);
00256   TH1F *fsum=new TH1F(*f);
00257   fsum->Add(f2);
00258   fsum->Add(f3);
00259   TH1F *g=new TH1F(*pion);
00260   f->Divide(g);
00261   f2->Divide(g);
00262   f3->Divide(g);
00263   fsum->Divide(g);
00264   fsum->Draw();
00265   fsum->SetTitle("bg photons per pion");
00266   f->Draw("same");
00267   f->SetLineColor(2);
00268   f2->Draw("same");
00269   f2->SetLineColor(4);
00270   f3->Draw("same");
00271   f3->SetLineColor(3);
00272   fsum->SetMinimum(.001);
00273   c->cd(6);
00274   TH1F *etatopi=new TH1F(*etameson);
00275   TH1F *omegatopi=new TH1F(*omega);
00276   etatopi->Divide(pion);
00277   omegatopi->Divide(pion);
00278   omegatopi->SetTitle("ratio of #eta/#omega to #pi");
00279   omegatopi->Draw();
00280   etatopi->Draw("same");
00281   c->cd(0);
00282   c->SaveAs((mFileName+".ps").Data());
00283   c->SaveAs((mFileName+".root").Data());
00284 
00285   TFile *dec_out=new TFile((mFileName+"Sum.root").Data(),"recreate");
00286   fsum->Write();
00287   f->Write("gamma_pion");
00288   f2->Write("gamma_eta");
00289   f3->Write("gamma_omega");
00290   dec_out->Close();
00291 }
00292 
00293 Double_t MyDecay::ptweight(Double_t x)
00294 {
00295   float weight=1.;
00296   if(x<0.2) x=0.2;
00297   if(isPP05){
00298     weight=300.*x*pow(1.+x,-9.0);//error: +/- 1.0%
00299   }
00300   if(isDAU){
00301     weight=8.24601e-02*x*pow(1.+x,-9.01897e+00);//error: +/- 1.0%
00302   }
00303   return weight;
00304 }
00305 Double_t MyDecay::ptweightETA(Double_t x)
00306 {
00307   Float_t ME=0.548;//eta meson mass
00308   Float_t MP=0.135;//pi0 mass
00309   Float_t weight=fETAMTSCALE * x/sqrt(x*x + ME*ME - MP*MP) * ptweight(sqrt(x*x + ME*ME - MP*MP));
00310   return weight;
00311 }
00312 Double_t MyDecay::ptweightOMEGA(Double_t x)
00313 {
00314   Float_t MO=0.782;//omega mass
00315   Float_t MP=0.135;//pi0 mass
00316   Float_t weight=1. * x/sqrt(x*x + MO*MO - MP*MP) * ptweight(sqrt(x*x + MO*MO - MP*MP));
00317   return weight;
00318 }
00319 Double_t MyDecay::etaweight(Double_t /*x*/)
00320 {
00321   return 1.;
00322 }  
00323 Double_t MyDecay::phiweight(Double_t /*x*/)
00324 {
00325   return 1.;
00326 }