Generating Madgraph Events

Information on how to generate madgraph events over a farm (and restart the failed jobs!) can be found here

So the analysis

Ok so I've chained together my 5 .ntpl files (output of CMKIN), into a file called analysis.root (which is what the h101 object uses) and I now want to analyise them, here's how.

Open root:

> root 

Instantiate a h101 object and run its loop method:

// This loads the class files from the directory above

root [3] .L ../h101.C

// This makes our h101 object

root [4] h101 c

// This runs the loop method of the h101 object, which makes the graphs etc.

root [5] c.Loop()

Analysis Code

The code here failed to work initially, I'm assuming it's faily old, and probably a bit out of date. Below is my updated version, it works for me, and produces four graphs.

#define h101_cxx

#include "h101.h"

#include 

#include 

#include 

                                                                            

void h101::Loop()

{

//   In a ROOT session, you can do:

//      Root > .L h101.C

//      Root > h101 t

//      Root > t.GetEntry(12); // Fill t data members with entry number 12

//      Root > t.Show();       // Show values of entry 12

//      Root > t.Show(16);     // Read and show values of entry 16

//      Root > t.Loop();       // Loop on all entries

//

                                                                            

//     This is the loop skeleton where:

//    jentry is the global entry number in the chain

//    ientry is the entry number in the current Tree

//  Note that the argument to GetEntry must be:

//    jentry for TChain::GetEntry

//    ientry for TTree::GetEntry and TBranch::GetEntry

//

//       To read only selected branches, Insert statements like:

// METHOD1:

//    fChain->SetBranchStatus("*",0);  // disable all branches

//    fChain->SetBranchStatus("branchname",1);  // activate branchname

// METHOD2: replace line

//    fChain->GetEntry(jentry);       //read all branches

//by  b_branchname->GetEntry(ientry); //read only this branch

                                                                            

Int_t num_events = h101->GetEntries();

Float_t eta, theta, p, pt;

Int_t num_Z = 0;

Int_t num_g = 0;

// Make an empty histogram named "eta"

// (declared on stack memory so it can be used after the fuction is over)

                                                                            

cout << "\t Number of Events = " << num_events << endl;

cout << "\t Nhep = " << this->Nhep << endl;

TH1F *myZHisto1 = new TH1F("Z eta","Z in CMS eta (H0 150)",25, -5, 5);

TH1F *myZHisto2 = new TH1F("Z theta","Z in CMS theta (H0 150)", 10, -3, 3);

TH1F *mygHisto1 = new TH1F("g eta","gluons in CMS eta (H0 150)",25, -5, 5);

TH1F *mygHisto2 = new TH1F("g theta","glouns in CMS theta (H0 150)", 10, -3, 3);

                                                                                                      

   for (Int_t i=0; i< num_events;i++)

   {

      GetEntry(i);  // This sets pointer to a particular Event #

                                                                                                      

      for (Int_t j=0; j < this->Nhep; j++)

      {

                                                                                                      

          if(abs(this->Idhep[j]) == 23) //Sets Particle type (Idhep)

         {

                //cout << "\t Event (" << i << ", " << j << ") has ID = " << this->Idhep[j] << endl;

                //cout << "Px=" << Phep[0][j]<< " Py=" << Phep[1][j]<< " Pz=" << Phep[2][j]<                p = sqrt(this->Phep[0][j]*this->Phep[0][j] + this->Phep[1][j]*this->Phep[1][j] + this->Phep[2][j]*this->Phep[2][j]);

                pt = sqrt(this->Phep[0][j]*this->Phep[0][j] + this->Phep[1][j]*this->Phep[1][j]);

                //cout << "p=" << p << endl << "pt=" << pt << endl;

                theta = asin(pt/p);

                eta = -log(tan(theta/2));

                if(this->Phep[2][j] < 0) eta = -1*eta;  // put sign on theta

                if(this->Phep[2][j] < 0) theta = -1*theta;

                //cout << "====>  eta for " << this->Idhep[j] << "\t is " << eta << endl;

                myZHisto1->Fill(eta);

                myZHisto2->Fill(theta);

                num_Z++;

          } // end if particle type

          if(abs(this->Idhep[j]) == 21) //Sets Particle type (Idhep)

         {

                if (Phep[0][j] != 0){

                        //cout << "\t Event (" << i << ", " << j << ") has ID = " << this->Idhep[j] << endl;

                        //cout << "Px=" << Phep[0][j]<< " Py=" << Phep[1][j]<< " Pz=" << Phep[2][j]<                        p = sqrt(this->Phep[0][j]*this->Phep[0][j] + this->Phep[1][j]*this->Phep[1][j] + this->Phep[2][j]*this->Phep[2][j]);

                        pt = sqrt(this->Phep[0][j]*this->Phep[0][j] + this->Phep[1][j]*this->Phep[1][j]);

                        //cout << "p=" << p << endl << "pt=" << pt << endl;

                        theta = asin(pt/p);

                        eta = -log(tan(theta/2));

                        if(this->Phep[2][j] < 0) eta = -1*eta;  // put sign on theta

                        if(this->Phep[2][j] < 0) theta = -1*theta;

                        //cout << "====>  eta for " << this->Idhep[j] << "\t is " << eta << endl;

                        mygHisto1->Fill(eta);

                        mygHisto2->Fill(theta);

                        num_g++;

                }

          } // end if particle type

                                                                                                      

      } // end loop on particles(j)

   } // end loop on events (i)

                                                                                                      

  TCanvas *can1 = new TCanvas("can1", "Z Eta", 1);

  TCanvas *can2 = new TCanvas("can2", "Z Theta", 1);

  TCanvas *can3 = new TCanvas("can3", "g Eta", 1);

  TCanvas *can4 = new TCanvas("can4", "g Theta", 1);

                                                                                                      

//Plot Z's

  can1->cd();

  myZHisto1->Draw();

  can2->cd();

  myZHisto2->Draw();

  can3->cd();

  mygHisto1->Draw();

  can4->cd();

  mygHisto2->Draw();

                                                                                                      

  cout << "number of Z's = " << num_Z << endl << "number of g's = " << num_g << endl;

                                                                                                      

} // END myclass::Loop

Well started running root analysis yesterday. Quite a learning curve. Take output from madgraph, run through CMKIN to get ntpls then h2root to get root files. Once in root chain the files together:

root [0] TChain *myChain = new TChain("h101","");
root [1] myChain->Add("*.root");
root [2] myChain->Merge("bigfile2.root")

now open the file up to look at quantities etc.

first put it in a tree:
TTree *T = _file0->Get("h101;1")

then following from this analyise some quantities

(I posted this before but it failed to work :( )

Toot toot, time to use root

Simon At Home: yay - it's created simon.root
Jim: now you need to fire up Root and read the file
Jim: find out what's in it, and try plotting some stuff
Simon At Home: :shudder:
Jim: for example, generator jet et distribution
Simon At Home: ...time to learn root :)
Jim: yep, and to work out the meaning of all the jetmet quantities ;-)

Now to load my root file and plot some graphs.....

Trying again

With advice from Jim:
"- create fresh orca area and check out JetMetAnalysis
- build in domain/test
- extract .orcarc contents from run_db.csh
- run JetMet"

So deleted the Orca directory, redid the check out etc.
- cd JetMetAnalysis/domain/test
- eval `scram runtime -csh`
- scram build JetMet
- rehash

Set up .orcarc (link)
- JetMet

Result: It looks to have worked! It's created simon.root at least.... (log file)

Nope it failed

It really looked like it was doing something but get this:

---

------ Information on L1 jets :
*** isolated em objects : 0
*** central jets : 1
*** forward jets : 3
*** tau objects : 1
*** global objects : 1
------ MEMORY use statistics -------
CobraSignal: SEGV received.
in thread 14331 2051
0x42ad0b8d BackTrace::trace(std::ostream&) const + 0x2d [/afs/cern.ch/cms/Releases/COBRA/COBRA_7_8_6/lib/Linux__2.4/libGenUtil.so]
0x42a34156 CobraSignalHandler::TermSignalHandler(int) + 0x246 [/afs/cern.ch/cms/Releases/COBRA/COBRA_7_8_6/lib/Linux__2.4/libBaseMonitor.so]
0x43bb0f05 [/lib/i686/libpthread.so.0]
0x43be5188 [/lib/i686/libc.so.6]
0x43c36d24 __libc_free + 0xb4 [/lib/i686/libc.so.6]
0x43b58953 operator delete(void*) + 0x23 [/usr/local/gcc-alt-3.2.3/lib/libstdc++.so.5]
0x4204cdca ConcreteJet::~ConcreteJet() + 0x66 [/afs/cern.ch/cms/Releases/ORCA/ORCA_8_3_0/lib/Linux__2.4/libJetTypes.so]
0x42356fee TReconstructor::~TReconstructor() + 0x4e [/afs/cern.ch/cms/Releases/COBRA/COBRA_7_8_6/lib/Linux__2.4/libReco.so]
0x4235542e TRecEvent::~TRecEvent() + 0x8e [/afs/cern.ch/cms/Releases/COBRA/COBRA_7_8_6/lib/Linux__2.4/libReco.so]
0x422f8326 TRecEventWP::~TRecEventWP() + 0xa6 [/afs/cern.ch/cms/Releases/COBRA/COBRA_7_8_6/lib/Linux__2.4/libEvent.so]
0x4226ffc8 G3EventProxy::reset() + 0x418 [/afs/cern.ch/cms/Releases/COBRA/COBRA_7_8_6/lib/Linux__2.4/libG3Event.so]
0x4226f85e G3EventProxy::~G3EventProxy() + 0x3e [/afs/cern.ch/cms/Releases/COBRA/COBRA_7_8_6/lib/Linux__2.4/libG3Event.so]
0x4219f64a SimpleRecEventProcessor::clearEvent() + 0xaa [/afs/cern.ch/cms/Releases/COBRA/COBRA_7_8_6/lib/Linux__2.4/libRecReader.so]
0x42197698 TEventReader::nextEvent(pool::Ref const&) + 0x48 [/afs/cern.ch/cms/Releases/COBRA/COBRA_7_8_6/lib/Linux__2.4/libRecReader.so]
0x421938a5 RecEventProcessor::consume(std::pair > const + 0xa5 [/afs/cern.ch/cms/Releases/COBRA/COBRA_7_8_6/lib/Linux__2.4/libRecReader.so]
0x42171eaf CollectionProcessingThread::run() + 0xff [/afs/cern.ch/cms/Releases/COBRA/COBRA_7_8_6/lib/Linux__2.4/libBatchRecReader.so]
0x42afd9e3 Thread::operator()() + 0x313 [/afs/cern.ch/cms/Releases/COBRA/COBRA_7_8_6/lib/Linux__2.4/libGenUtil.so]
0x42afe3b1 boost::detail::function::void_function_obj_invoker0::invoke(boost::detail::functio + 0x21 [/afs/cern.ch/cms/Releases/COBRA/COBRA_7_8_6/lib/Linux__2.4/libGenUtil.so]
0x43a41b8a [/afs/cern.ch/cms/external/lcg/external/Boost/1.30.2/rh73_gcc32/lib/libboost_thread.so]
0x43badfaf [/lib/i686/libpthread.so.0]
0x43ca392a __clone + 0x3a [/lib/i686/libc.so.6]

Forced KILL,

Killed

---

Balls!

Google is your friend

It seems that all the orca documentation is spread out over the net and not really connected or anything, so google is your friend (cheers Jim for the pointer).
So I was getting: "CobraSignal: SEGV received." googled it, and it turns out its caused by not setting the OutputDataSet in .orcarc. Just set that and it's now running - taking ages mind! - but it is doing stuff.....

Couple of useful links

USCMS JetMet group
The Algorithms

Starting from scratch

I want to do some JetMet analysis. To this end I'm going to check out from CVS the JetMet code and try and run their example jobs on lxplus.
Lets see how it goes (sort of following this):

1. scram project ORCA ORCA_8_3_0
   
2. cd ORCA_8_3_0 /src
   
3. eval `scram runtime -csh`
   
4. cmscvsroot ORCA
   
5. cvs login
   
6. cvs co -r ORCA_8_3_0 JetMetAnalysis

That sets up the workspace and gets all the relevant files....

Now to try and build some executables, and maybe even a root tree!

1. cd JetMetAnalysis/JetMetTest/test
   
2. scram b shared > JetMetTest_scram_b_shared.log & (to store it for posterity)
   
3. scram b bin > JetMetTest_scram_b_bin.log &
   
4. eval `scram runtime -csh`
   
5. rehash

Now to run an exe:

1. JetMetTest (to run the exe - this will probably fail as I have no .orcarc set up)
   
2. Failed with a seg fault - no data to run on time to set the orcarc
   
3. vi .orcarc and enter the following:
   
   PoolCatalogFile = @{xmlcatalog_http://cmsdoc.cern.ch/orca/catalog/PoolFileCatalog_8_3_0.xml}@
   InputCollections=/System/StW830DST2x1033/jets50100/jets50100
   MaxEvents=10
   
4. JetMetTest > JetMetTest.log &
   
5. Takes a while....
   
6. ...but it is doing something...
   
7. ...and then it fails:

Try with a different exe JetMetSmallTest:

1. JetMetSmallTest > JetMetSmallTest.log
2. runs and doesn't crash, but produces no output

First Post

Well I'm struggling to learn ORCA which is bordering on a total nightmare! So going to use this blog as an online work book and keep track of the steps I do, you never know it might come in handy for thesis or for some other poor soul learning the CMS analysis framework.