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 R__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 "TEveTrack.h"
#include "TEveTrackPropagator.h"
#include "TEveElement.h"
#include "TEveGeoShape.h"
 
#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 = static_cast<TEveElement *>(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