pythia_display.C

Go to the documentation of this file.

/// \file
/// \ingroup tutorial_eve
/// Demo showing H -> ZZ -> 4 mu generated by Pythia.
/// Requires libPythia6.
///
/// \macro_code
///
/// \author Matevz Tadel

#ifndef __RUN_PYTHIA_DISPLAY__

void pythia_display()
{
   TString dir = gSystem->UnixPathName(__FILE__);
   dir.ReplaceAll("pythia_display.C","");
   dir.ReplaceAll("/./","/");
   gROOT->LoadMacro(dir +"MultiView.C+");

#ifndef G__WIN32 // libPythia6 is a static library on Windoze
   if (gSystem->Load("libPythia6") < 0)
   {
      Error("pythia_display()",
            "Could not load 'libPythia6', make sure it is available!");
      return;
   }
#endif
   gSystem->Load("libEGPythia6");

   gROOT->ProcessLine("#define __RUN_PYTHIA_DISPLAY__ 1");
   gROOT->ProcessLine("#include \"pythia_display.C\"");
   gROOT->ProcessLine("run_pythia_display()");
   gROOT->ProcessLine("#undef __RUN_PYTHIA_DISPLAY__");
}

#else

//==============================================================================
// Constants.
//------------------------------------------------------------------------------

const Double_t kR_min = 240;
const Double_t kR_max = 250;
const Double_t kZ_d   = 300;

// Solenoid field along z, in Tesla.
const Double_t kMagField = 4;

// Color for Higgs, Zs and muons
const Color_t  kColors[3] = { kRed, kGreen, kYellow };

//==============================================================================
// Global variables.
//------------------------------------------------------------------------------

#include "TPythia6.h"
#include "TGeoTube.h"
#include "TMCParticle.h"

TPythia6   *g_pythia = 0;

// Implemented in MultiView.C
class MultiView;
MultiView* gMultiView = 0;

TEveTrackList *gTrackList = 0;

//==============================================================================
// Forward decalarations of CINT functions.
//------------------------------------------------------------------------------

void pythia_next_event();
void pythia_make_gui();

//==============================================================================
// Main - pythia_display()
//------------------------------------------------------------------------------

void run_pythia_display()
{
   if (g_pythia != 0)
   {
      Warning("pythia_display()", "Already initialized.");
      return;
   }

   //========================================================================
   //========================================================================

   // Create an instance of the Pythia event generator ...
   g_pythia = new TPythia6;
   TPythia6& P = * g_pythia;

   P.SetMSEL(0);           // full user controll;
   P.SetMSUB(102, 1);      // g + g -> H0
   //P.SetMSUB(123, 1);    // f + f' -> f + f' + H0
   //P.SetMSUB(124, 1);    // f + f' -> f" + f"' + H0

   P.SetPMAS(6,  1, 175);  // mass of TOP
   P.SetPMAS(25, 1, 180);  // mass of Higgs


   P.SetCKIN(1, 170.0);    // range of allowed mass
   P.SetCKIN(2, 190.0);

   P.SetMSTP(61, 0);   // switch off ISR
   P.SetMSTP(71, 0);   // switch off FSR
   P.SetMSTP(81, 0);   // switch off multiple interactions

   P.SetMSTP(111, 0);  // Switch off fragmentation

   // Force h0 -> ZZ
   for (Int_t i = 210; i <= 288; ++i)
      P.SetMDME(i, 1, 0);
   P.SetMDME(225, 1, 1);

   // Force Z -> mumu
   for (Int_t i = 174; i <= 189; ++i)
      P.SetMDME(i, 1, 0);
   P.SetMDME(184, 1, 1);


   P.Initialize("cms", "p", "p", 14000);

   //========================================================================
   // Create views and containers.
   //========================================================================

   TEveManager::Create();

   TEveElementList *fake_geom = new TEveElementList("Geometry");

   TEveGeoShape *b;

   b = new TEveGeoShape("Barell 1");
   b->SetShape(new TGeoTube(kR_min, kR_max, kZ_d));
   b->SetMainColor(kCyan);
   b->SetMainTransparency(80);
   fake_geom->AddElement(b);

   b = new TEveGeoShape("Barell 2");
   b->SetShape(new TGeoTube(2*kR_min, 2*kR_max, 2*kZ_d));
   b->SetMainColor(kPink-3);
   b->SetMainTransparency(80);
   fake_geom->AddElement(b);

   gEve->AddGlobalElement(fake_geom);


   gMultiView = new MultiView;

   gMultiView->ImportGeomRPhi(fake_geom);
   gMultiView->ImportGeomRhoZ(fake_geom);

   gEve->GetBrowser()->GetTabRight()->SetTab(1);

   gTrackList = new TEveTrackList("Pythia Tracks");
   gTrackList->SetMainColor(kYellow);
   gTrackList->SetMarkerColor(kRed);
   gTrackList->SetMarkerStyle(4);
   gTrackList->SetMarkerSize(0.5);
   gEve->AddElement(gTrackList);

   TEveTrackPropagator* trkProp = gTrackList->GetPropagator();
   trkProp->SetMagField(kMagField);
   trkProp->SetMaxR(2*kR_max);
   trkProp->SetMaxZ(2*kZ_d);

   //========================================================================
   //========================================================================

   pythia_make_gui();
   pythia_next_event();

   gEve->Redraw3D(kTRUE);
}


//==============================================================================
// Next event
//------------------------------------------------------------------------------

void pythia_next_event()
{
   gTrackList->DestroyElements();

   TPythia6& P = * g_pythia;

   P.GenerateEvent();

   int nh = P.GetMSTU(72);

   // printf("N = %d, Nhard = %d :: NumSec = %d, separators (%d,%d,%d,%d)\n",
   //    P.GetN(), nh, P.GetMSTU(70), P.GetMSTU(71), P.GetMSTU(72), P.GetMSTU(73), P.GetMSTU(74));
   //                          2->2                 hard                  postfrag              final

   TEveTrackPropagator *trkProp = gTrackList->GetPropagator();
   TClonesArray        &MC      = * (TClonesArray*) P.GetListOfParticles();
   for (Int_t i = 0; i < 7; ++i)
   {
      TMCParticle& p = (TMCParticle&)*MC[nh+i];
      TParticle pb(p.GetKF(), p.GetKS(), 0, 0,
                   p.GetFirstChild()-nh-1, p.GetLastChild()-nh-1,
                   p.GetPx(), p.GetPy(), p.GetPz(), p.GetEnergy(),
                   p.GetVx(), p.GetVy(), p.GetVz(), p.GetTime());

      TEveTrack* track = new TEveTrack(&pb, i, trkProp);
      track->SetName(Form("%s [%d]", pb.GetName(), i));
      track->SetStdTitle();
      track->SetAttLineAttMarker(gTrackList);
      if (i == 0)
         track->SetLineColor(kColors[0]);
      else if (i <= 2)
         track->SetLineColor(kColors[1]);

      gTrackList->AddElement(track);

      /*
        printf("%d - %d %d %d %d %d %d\n", i,
        p.GetKF(), p.GetKS(), 0, 0,
        p.GetFirstChild()-nh-1, p.GetLastChild()-nh-1);
        printf("%d - %f %f %f %f\n", i,
        p.GetPx(), p.GetPy(), p.GetPz(), p.GetEnergy(),
        printf("%d - %f %f %f %f\n", i,
        p.GetVx(), p.GetVy(), p.GetVz(), p.GetTime());
      */
   }

   gTrackList->MakeTracks();


   TEveElement* top = gEve->GetCurrentEvent();

   gMultiView->DestroyEventRPhi();
   gMultiView->ImportEventRPhi(top);

   gMultiView->DestroyEventRhoZ();
   gMultiView->ImportEventRhoZ(top);

   gEve->Redraw3D();
}


//==============================================================================
// GUI stuff
//------------------------------------------------------------------------------
class EvNavHandler
{
public:
   void Fwd()
   {
      pythia_next_event();
   }
   void Bck()
   {}
};

//______________________________________________________________________________
void pythia_make_gui()
{
   // Create minimal GUI for event navigation.

   TEveBrowser* browser = gEve->GetBrowser();
   browser->StartEmbedding(TRootBrowser::kLeft);

   TGMainFrame* frmMain = new TGMainFrame(gClient->GetRoot(), 1000, 600);
   frmMain->SetWindowName("XX GUI");
   frmMain->SetCleanup(kDeepCleanup);

   TGHorizontalFrame* hf = new TGHorizontalFrame(frmMain);
   {

      TString icondir( Form("%s/icons/", gSystem->Getenv("ROOTSYS")) );
      TGPictureButton* b = 0;
      EvNavHandler    *fh = new EvNavHandler;

      b = new TGPictureButton(hf, gClient->GetPicture(icondir+"GoBack.gif"));
      b->SetEnabled(kFALSE);
      b->SetToolTipText("Go to previous event - not supported.");
      hf->AddFrame(b);
      b->Connect("Clicked()", "EvNavHandler", fh, "Bck()");

      b = new TGPictureButton(hf, gClient->GetPicture(icondir+"GoForward.gif"));
      b->SetToolTipText("Generate new event.");
      hf->AddFrame(b);
      b->Connect("Clicked()", "EvNavHandler", fh, "Fwd()");
   }
   frmMain->AddFrame(hf);

   frmMain->MapSubwindows();
   frmMain->Resize();
   frmMain->MapWindow();

   browser->StopEmbedding();
   browser->SetTabTitle("Event Control", 0);
}

#endif