@@ -29,6 +29,8 @@ ExodusDecayer::ExodusDecayer():
2929 fEPMassPhiDalitz(0 ),
3030 fEPMassPhiDalitz_toPi0(0 ),
3131 fEPMassJPsi(0 ),
32+ fEPMassPsi2S(0 ),
33+ fEPMassUpsilon(0 ),
3234 fPol(new TF1(" dsigdcostheta" ," 1.+[0]*x*x" ,-1 .,1 .)), /* Polarization Function for resonances */
3335 fInit(0 ),
3436 fDecayToDimuon(0 )
@@ -51,6 +53,8 @@ ExodusDecayer::~ExodusDecayer()
5153 delete fEPMassPhiDalitz ;
5254 delete fEPMassPhiDalitz_toPi0 ;
5355 delete fEPMassJPsi ;
56+ delete fEPMassPsi2S ;
57+ delete fEPMassUpsilon ;
5458 delete fPol ;
5559}
5660
@@ -74,16 +78,16 @@ void ExodusDecayer::Init()
7478
7579 Float_t binwidth;
7680 Float_t mass_bin, mass_min, mass_max;
77- Double_t vmass_eta, vmass_rho, vmass_omega, vmass_phi, vmass_jpsi, vwidth_eta, vwidth_rho, vwidth_omega, vwidth_phi, vwidth_jpsi;
78- Double_t weight_eta, weight_rho, weight_omega, weight_phi, weight_jpsi;
81+ Double_t vmass_eta, vmass_rho, vmass_omega, vmass_phi, vmass_jpsi, vmass_psi2s, vmass_upsilon, vwidth_eta, vwidth_rho, vwidth_omega, vwidth_phi, vwidth_jpsi, vwidth_psi2s, vwidth_upsilon ;
82+ Double_t weight_eta, weight_rho, weight_omega, weight_phi, weight_jpsi, weight_psi2s, weight_upsilon ;
7983
8084// ================================================================================//
8185// Create electron pair mass histograms from dalitz decays //
8286// ================================================================================//
8387
8488 // Get the particle masses
8589 // parent
86- nbins = 1000 ;
90+ nbins = 2000 ;
8791 mass_min = 0 .;
8892 mass_max = 0 .;
8993 pionmass = (TDatabasePDG::Instance ()->GetParticle (111 ))->Mass ();
@@ -315,29 +319,41 @@ void ExodusDecayer::Init()
315319 vmass_omega = (TDatabasePDG::Instance ()->GetParticle (223 ))->Mass ();
316320 vmass_phi = (TDatabasePDG::Instance ()->GetParticle (333 ))->Mass ();
317321 vmass_jpsi = (TDatabasePDG::Instance ()->GetParticle (443 ))->Mass ();
322+ vmass_psi2s = (TDatabasePDG::Instance ()->GetParticle (100443 ))->Mass ();
323+ vmass_upsilon = (TDatabasePDG::Instance ()->GetParticle (553 ))->Mass ();
318324 // Get the particle widths
319325 // parent
320326 vwidth_eta = (TDatabasePDG::Instance ()->GetParticle (221 ))->Width ();
321327 vwidth_rho = (TDatabasePDG::Instance ()->GetParticle (113 ))->Width ();
322328 vwidth_omega = (TDatabasePDG::Instance ()->GetParticle (223 ))->Width ();
323329 vwidth_phi = (TDatabasePDG::Instance ()->GetParticle (333 ))->Width ();
324330 vwidth_jpsi = (TDatabasePDG::Instance ()->GetParticle (443 ))->Width ();
331+ vwidth_psi2s = (TDatabasePDG::Instance ()->GetParticle (100443 ))->Width ();
332+ vwidth_upsilon = (TDatabasePDG::Instance ()->GetParticle (553 ))->Width ();
325333
326334 if ( vwidth_jpsi < 1e-6 ){
327335 printf (Form (" ExodusDecayer: Warning: Width JPsi = %f (Mass = %f) < 1e-6, set to 1e-6\n " ,vwidth_jpsi,vmass_jpsi));
328336 vwidth_jpsi = 1e-6 ;
329337 }
338+ if ( vwidth_psi2s < 1e-6 ){
339+ printf (Form (" ExodusDecayer: Warning: Width Psi2S = %f (Mass = %f) < 1e-6, set to 1e-6\n " ,vwidth_psi2s,vmass_psi2s));
340+ vwidth_psi2s = 1e-6 ;
341+ }
342+ if ( vwidth_upsilon < 1e-6 ){
343+ printf (Form (" ExodusDecayer: Warning: Width Upsilon = %f (Mass = %f) < 1e-6, set to 1e-6\n " ,vwidth_upsilon,vmass_upsilon));
344+ vwidth_upsilon = 1e-6 ;
345+ }
330346
331347
332348 if ( mass_min == 0 . && mass_max == 0 . )
333349 {
334350 mass_min = 2 .*pimass;
335- mass_max = 5 .;
351+ mass_max = 10 .;
336352 }
337353
338354 binwidth = (mass_max-mass_min)/(Double_t)nbins;
339355
340- // create pair mass histograms for resonances of rho, omega, phi and jpsi
356+ // create pair mass histograms for resonances of rho, omega, phi, jpsi, psi2s and upsilon
341357 fEPMassEta = new TH1F (" fEPMassEta" ," mass eta" ,nbins,mass_min,mass_max);
342358 fEPMassEta ->SetDirectory (0 );
343359 fEPMassRho = new TH1F (" fEPMassRho" ," mass rho" ,nbins,mass_min,mass_max);
@@ -348,6 +364,10 @@ void ExodusDecayer::Init()
348364 fEPMassPhi ->SetDirectory (0 );
349365 fEPMassJPsi = new TH1F (" fEPMassJPsi" ," mass jpsi" ,nbins,mass_min,mass_max);
350366 fEPMassJPsi ->SetDirectory (0 );
367+ fEPMassPsi2S = new TH1F (" fEPMassPsi2S" ," mass psi2s" ,nbins,mass_min,mass_max);
368+ fEPMassPsi2S ->SetDirectory (0 );
369+ fEPMassUpsilon = new TH1F (" fEPMassUpsilon" ," mass upsilon" ,nbins,mass_min,mass_max);
370+ fEPMassUpsilon ->SetDirectory (0 );
351371
352372 for (ibin=1 ; ibin<=nbins; ibin++ )
353373 {
@@ -360,13 +380,17 @@ void ExodusDecayer::Init()
360380 weight_omega = (Float_t)GounarisSakurai (mass_bin,vmass_omega,vwidth_omega,emass);
361381 weight_phi = (Float_t)GounarisSakurai (mass_bin,vmass_phi,vwidth_phi,emass);
362382 weight_jpsi = (Float_t)Lorentz (mass_bin,vmass_jpsi,vwidth_jpsi);
383+ weight_psi2s = (Float_t)Lorentz (mass_bin,vmass_psi2s,vwidth_psi2s);
384+ weight_upsilon = (Float_t)Lorentz (mass_bin,vmass_upsilon,vwidth_upsilon);
363385
364386 // Fill histograms of electron pair masses from resonance decays
365387 fEPMassEta ->AddBinContent (ibin,weight_eta);
366388 fEPMassRho ->AddBinContent (ibin,weight_rho);
367389 fEPMassOmega ->AddBinContent (ibin,weight_omega);
368390 fEPMassPhi ->AddBinContent (ibin,weight_phi);
369391 fEPMassJPsi ->AddBinContent (ibin,weight_jpsi);
392+ fEPMassPsi2S ->AddBinContent (ibin,weight_psi2s);
393+ fEPMassUpsilon ->AddBinContent (ibin,weight_upsilon);
370394 }
371395 fInit =1 ;
372396}
@@ -478,7 +502,11 @@ void ExodusDecayer::Decay(Int_t idpart, TLorentzVector* pparent)
478502
479503 // -------- get id of the partner from idpart ------- //
480504 Int_t idpartner=idpart/1000 ;
481- idpart=idpart%1000 ;
505+ if (idpart!=100443 ){ // protect psi(2S) pdg code (100443)
506+ idpart=idpart%1000 ;
507+ } else {
508+ idpartner=0 ;
509+ }
482510
483511 // local variables for dalitz/2-body decay:
484512 Double_t pmass, epmass, realp_mass, e1 , p1, e3 , p3;
@@ -492,6 +520,8 @@ void ExodusDecayer::Decay(Int_t idpart, TLorentzVector* pparent)
492520 Int_t idPi0=111 ;
493521 Int_t idEta=221 ;
494522 Int_t idEtaPrime=331 ;
523+ Int_t idPsi2S=100443 ;
524+ Int_t idUpsilon=553 ;
495525
496526 // Get the particle masses of daughters
497527 Double_t emass, proton_mass, omass_pion, omass_eta, omass_gamma, omass_omega;
@@ -651,11 +681,9 @@ void ExodusDecayer::Decay(Int_t idpart, TLorentzVector* pparent)
651681
652682
653683// -----------------------------------------------------------------------------//
654- // Generate 2-body resonance decays: Rho/Omega/Phi/JPsi //
684+ // Generate 2-body resonance decays: Rho/Omega/Phi/JPsi/Psi2S/Upsilon //
655685// -----------------------------------------------------------------------------//
656-
657- if (((idpart==idEta&&fDecayToDimuon ==1 )||idpart==idRho||idpart==idOmega||idpart==idPhi||idpart==idJPsi)&&(idpartner==0 )){
658-
686+ if (((idpart==idEta&&fDecayToDimuon ==1 )||idpart==idRho||idpart==idOmega||idpart==idPhi||idpart==idJPsi||idpart==idPsi2S||idpart==idUpsilon)&&(idpartner==0 )){
659687
660688 // get the parent mass
661689 mp_res = pparent->M ();
@@ -666,7 +694,6 @@ void ExodusDecayer::Decay(Int_t idpart, TLorentzVector* pparent)
666694 printf (" ExodusDecayer: res into ee Decay kinematically impossible! \n " );
667695 return ;
668696 }
669-
670697 // Sample the electron pair mass from a histogram and set Polarization
671698 for ( ;; ) {
672699 if (idpart==idEta&&fDecayToDimuon ==1 ){
@@ -687,6 +714,12 @@ void ExodusDecayer::Decay(Int_t idpart, TLorentzVector* pparent)
687714 }else if (idpart==idJPsi){
688715 epmass_res = fEPMassJPsi ->GetRandom ();
689716 PolPar=0 .;
717+ }else if (idpart==idPsi2S){
718+ epmass_res = fEPMassPsi2S ->GetRandom ();
719+ PolPar=0 .;
720+ } else if (idpart==idUpsilon){
721+ epmass_res = fEPMassUpsilon ->GetRandom ();
722+ PolPar=0 .;
690723 }else { printf (" ExodusDecayer: ERROR: Resonance mass G-S parametrization not found \n " );
691724 return ;
692725 }
@@ -757,6 +790,12 @@ void ExodusDecayer::Decay(Int_t idpart, TLorentzVector* pparent)
757790 }else if (idpart==idJPsi){
758791 fProducts_jpsi [0 ]=fProducts_res [0 ];
759792 fProducts_jpsi [1 ]=fProducts_res [1 ];
793+ }else if (idpart==idPsi2S){
794+ fProducts_psi2s [0 ]=fProducts_res [0 ];
795+ fProducts_psi2s [1 ]=fProducts_res [1 ];
796+ }else if (idpart==idUpsilon){
797+ fProducts_upsilon [0 ]=fProducts_res [0 ];
798+ fProducts_upsilon [1 ]=fProducts_res [1 ];
760799 }
761800
762801 }
0 commit comments