What we want to do :

    Generate a bunch of particles via Monte Carlo in PYTHIA and do some sort of simple (but meaningful analysis).  Lets worry about the detector later ...
 

To generate a bunch of particles :

     Go to Han's cms production web page and click on the "Creating CMSIM Jobs" link. Read it.  Hans has made a way to generate PYTHIA scripts in a very easy way.  It will also allow you to easily feed these particles into the simulated CMS detector (CMSIM) later on.  Using Han's scripts also ensures that everything you do is comparable to everyone else in CMS.
 

     To setup PYTHIA and CMSIM (compliments of Han's cms_production.html web page):

        cd /home/leer/cmsim
        rm -rf cms_production
        source ~/cms_setup.csh
        cmsim cms120
        cvs -d :pserver:anonymous@cdcvs.fnal.gov:/cvs/cd_read_only login
        ENTER PASSWORD--> anoncvs
        cvs -d :pserver:anonymous@cdcvs.fnal.gov:/cvs/cd_read_only checkout cms_production
        cd cms_production
        cd src
        cmsbuild120
        kine_make_ntpl_pyt_120.com cms120
        cd cmsim_single/
        cmsbuild120
        ls ../../bin
        cd ../../
        cp ./scripts/setup.csh .

     Now you should make up a partial script to run.  The script Hans wrote will place all the necessary information at the front of the card and initialize PYTHIA to CMS conditions for you, but it is still necessary to specify exactly what process or processes you want to simulate.  To do that make a "title card" in either the btau, egamma, jetsmet, level1, minbias_no_ckin, or muon directory - which are located in your cms_production/Template directory.  You should carefully read "Creating CMSIM Jobs" on the cms production web page and remember that your file should be named <something>_gen.tit.

     If this is your first time running PYTHIA and you don't know much about the data card particle and/or decay codes, download the PYTHIA miniguide from the PYTHIA website and look at the section called "Practical Information".  You can also find a sample datacard here.  If you already know what your doing and want to edit or have a look at what else will be automatically put on your datacard, you can look at the cms_production/scripts/makeJob.pl file.

The format of the command to generate a script is :  setup.csh [dir] [script] [jobs] [#events] stuff...

So, to run a job (10 events) with only CMGEN (ie no CMSIM):

 setup.csh muon higgs_muon 1 10 800001 800000 800000 1 cms120 kine_make_ntpl_pyt_cms120.exe  none
 cd /home/leer/cmsim/cms_production/Projects/higgs_muon/pending
 mv ../scripts/* .
 higgs_muon_800001.run > & ../output/higgs_muon_800001.out &
 cd ../output

Note:  If you really want to make a huge number of events, you should use a more than one job with a smaller number of events.  For example, rather than run 10000 events, run 10 jobs with 1000 events each - there is appearently a max number of events you can run and it will stop running if you exceed that maximum.  You can always merge the ntuples or root files later anyway.

The 'none' entry at the end of the script indicates that we do not want to run cmsim with this job.  See "How To Simulate CMS Detector Response to Particles in CMSIM" if you are interested in turning on the detector simulation as well.

You should see your ntuple in the directory once the job is finished ... it is a standard HEPEVT ntuple.  You can find information about whats inside the ntuple and how it is organized here.

Now we want to look at what is inside the ntuple ...
 
 

To Analyze the ntuple in ROOT

To need to get your ntuple into root format :

     h2root yourfile.ntpl yourfile.root

Now we can look at the ntuple data in at least 3 ways.  I recommend the first way since it is the jumping off point for more complex things, the other two ways are things I stumbled upon which may or may not be helpful to you.

1) The "Real" way (full blown analysis with unlimited options)  -- Recommended
 

-- Start a root session

   Let ROOT make a class for you with (it will generate .C and .h files for you):

 TFile *f = new TFile("yourfile.root")
 h100->MakeClass("yourclass")

-- To look at your data, you can start ROOT again anytime and type

.L myclass.C             // Load the class you created earlier
myclass c                  // make an object based on the class ROOT created

c.Loop()                       // Execute Loop method for analysis
                                     // Inside here is where you can do analysis and sorting
 

-- If you were to look inside the myclass.C file you created, you will see one method, myclass::Loop
     You should edit this code to loop over all particles in myfile.root and do an analysis
     For example, I replaced the default myclass::Loop method with:

void myclass::Loop()
{

// Purpose: Create a histogram of all muon (initial) eta

Int_t num_events = h100->GetEntries();
Float_t eta, theta, p, pt;

// Make an empty histogram named "eta"
// (declared on stack memory so it can be used after the fuction is over)

TH1F *myHisto = new TH1F("eta","Muons in CMS eta (H0 150)", 100, -2.4, 2.4);

   for (Int_t i=0; i< num_events;i++)
   {
      GetEntry(i);  // This sets pointer to a particular Event #

      for (Int_t j=0; j< c.Nhep;j++)
      {
          if(abs(c.Idhep[j]) == 13 && c.Isthep[j] == 1)
         {
            p = sqrt(c.Phep[j][0]*c.Phep[j][0] + c.Phep[j][1]*c.Phep[j][1] + c.Phep[j][2]*c.Phep[j][2]);
            pt = sqrt(c.Phep[j][0]*c.Phep[j][0] + c.Phep[j][1]*c.Phep[j][1]);
            theta = asin(pt/p);
            eta = -log(tan(theta/2));
            if(c.Phep[j][2] < 0) eta = -1*eta;  // put sign on theta

            cout << "\t Event (" << i << ") ";
            cout << "====>  eta for " << c.Idhep[j] << "\t is " << eta << endl;

            myHisto->Fill(eta);
          } // end if muon
      } // end loop on particles(j)
   } // end loop on events (i)

 eta->Draw();

} // END myclass::Loop
 

You can also look at individual things like

c.GetEntry(0)           // Look into Event1
c.Idhep[10]                // Particle Type for entry 10 in 1st Event
c.Phep[0][3]              // Energy for 0th particle in 1st Event
c.GetEntry(1)           // Look into Event2
c.Phep[10][2]            // Pz for 10th particle in 2nd Event

2)  To make quick plots and see some basic info with only a couple cuts inside a Root session

Start a Root session:

TFile f("yourfile.root");
TTree *t100 = (TTree*)f->Get("h100")

To draw the E of 13th particle in 1st event
t100->Draw("Phep[12][3]","Ievt == 1")

To Draw the pz of the 12 particle in 1st event
t100->Draw("Phep[11][2]","Ievt == 1")

To Draw the pz of all particles in the 1st event
t100->Draw("Phep[][2]","Ievt == 1")

To Draw the E of all undecayed muons in 1st event
t100->Draw("Phep[][3]", "Ievt == 1 && Isthep == 1 && Idhep == -13")

To Draw the E of all Higgs immediately before their decay in 1st event
t100->Draw("Phep[][3]", "Ievt == 1 && Isthep == 1 && Idhep == -13")

Show pz entry of 0th particle (for all events)
t100->Scan("Phep[0][2]")

Show pz and E entries of 0th particle (for events)
t100->Scan("Phep[0][2]:Phep[0][3]")

Show pz and E entries of 0th particle for 1st event
t100->Scan("Phep[0][2]:Phep[0][3]", "Ievt == 1")
 
 

3)  To quicky get a particular variable of your data into an array in a Root session.
     (This trick only works for 3 variables at one time, 4th overwrites 1st)

All Energy's from 1st event
t100->Draw("Phep[][3]", "Ievt == 1")
b= t100->GetV1()

Test With (array index is particle number):
for(int i = 0; i < 10; i++) cout << b[i] << " "
cout << endl                                                        // to get the outout to screen