---

debugMCinfo.cpp

Go to the documentation of this file.

#include <XANADOO/domain/interface/XANAEsdEvent.h>
#include <XANADOO/XANAClusters/interface/XANACluster.h>
#include <XANADOO/XANAClusters/interface/XANACaloRecHit.h>
#include <XANADOO/XANAClusters/interface/XANAEmCluster.h>
#include <XANADOO/XANAClusters/interface/XANASuperCluster.h>
#include <XANADOO/XANAClusters/interface/XANAHadCluster.h>
#include <XANADOO/XANATracks/interface/XANATrack.h>
#include <XANADOO/XANAMcInfo/interface/XANAGeneratorEvent.h>
#include <XANADOO/XANAMcInfo/interface/XANAGeantEvent.h>
#include <XANADOO/XANAMcInfo/interface/XANAGeneratorParticle.h>
#include <XANADOO/XANAElectronCandidate/interface/XANAElectronCandidate.h>
#include <XANADOO/XANAElectronCandidate/interface/XANAElectronTrack.h>
#include <XANADOO/XANAElectronCandidate/interface/XANAElectronSeed.h>
#include <XANADOO/XANAMuons/interface/XANAMuonCandidate.h>

#include <XANADOO/XANAMcInfo/interface/XANAGeantTrack.h>
#include <XANADOO/XANAMcInfo/interface/XANAGeantVertex.h>

#include "TTree.h"
#include "TFile.h"
#include "TBranch.h"
#include "TClonesArray.h"
#include "TObjArray.h"
#include "TH1.h"
#include "TH2.h"
#include "TSystem.h"
#include "TROOT.h"

#include <string>


using namespace std;

int main(int argc, char** argv)
{
   //   Connect to file generated
   TFile *f;
   if (argc==1) f = new TFile("h300eemm.root");
   else f = new TFile(*++argv);

   //   Read Tree named "T" in memory. Tree pointer is assigned the same name
   TTree *T = (TTree*)f->Get("ESD");


   // Reconstructed event
   XANAEsdEvent *event = 0;  
   T->SetBranchAddress("Event",&event);

   // MC information (generator and GEANT event)
   XANAGeneratorEvent *generatorEvent = new XANAGeneratorEvent();
   TBranch *branchGeneratorEvent = T->GetBranch("GeneratorEvent");
   branchGeneratorEvent->SetAddress(&generatorEvent);

   XANAGeantEvent *geantEvent = new XANAGeantEvent();
   TBranch *branchGeantEvent = T->GetBranch("GeantEvent");
   branchGeantEvent->SetAddress(&geantEvent);

   Short_t mc_ngt = 0, mc_ngv = 0;  

   TFile *fout = new TFile("debugMC.root","RECREATE");
   TH2F *h_xy = new TH2F("x-y view","x-y view",400,-200.,200.,400,-200.,200.);
   TH2F *h_rz = new TH2F("r-z view","r-z view",800,-400.,400.,200,0.,200.);
   TH1F *h_r = new TH1F("r","r",200,0.,200.);
   
   //   Start main loop on all events
   Int_t nb = 0; 
   Int_t nevent = (Int_t)T->GetEntries();
   cout << "Going to read " << nevent << " events" << endl;
   for (Int_t i=0; i<nevent; i++) {
      nb += T->GetEntry(i); 

      cout << "================ Event: " << i << " ================" << endl;    


         // generator event
         cout << "--> new generator event" << endl;
         cout << endl;

         // GEANT event
         mc_ngt = geantEvent->getNumberOfGeantTracks();
         mc_ngv = geantEvent->getNumberOfGeantVertices();
         cout << "--> new GEANT event" << endl;
         cout << "      num of Geant tracks        = " << mc_ngt << endl;
         cout << "      num of Geant vertices      = " << mc_ngv << endl;
         cout << "    tracks ( piD --> momentum, mother index ):" << endl;
         
         XANAGeantTrack *geantTrack = 0;
         for (Int_t ig=0; ig<mc_ngt; ig++) {
            geantTrack = (XANAGeantTrack*)geantEvent->getGeantTracks()->At(ig);
            Int_t pID = geantTrack->getPartId();
            Int_t mother = geantTrack->getMotherIndex();
            cout << pID << "  -->  " << geantTrack->getMomentum() << ", " << mother<<  endl;
            XANAGeneratorParticle *genPart = geantTrack->getGeneratorParticle();
            if (genPart) {
              cout << " Parent generator particle " << genPart << endl;
              genPart->print();
            }
         }

         gROOT->SetStyle("Plain");
         // impact point histos
         for (Int_t ig=0; ig<mc_ngv; ig++) {
            XANAGeantVertex *geantVertex = (XANAGeantVertex*)geantEvent->getGeantVertices()->At(ig);
            Int_t mc_not = geantVertex->getNumberOfTracks();
            for (Int_t ig=0; ig<mc_not; ig++) {
               XANAGeantTrack *geantTrack = (XANAGeantTrack*)geantVertex->getTracks()->At(ig);
//             cout << geantTrack->getPartId() << endl;
               if( geantTrack->getPartId() == -103) {
                  Float_t x = geantVertex->getPosition().x();
                  Float_t y = geantVertex->getPosition().y();
                  Float_t z = geantVertex->getPosition().z();
                  if (z<275) {
                     h_xy->Fill(x,y);
                     h_rz->Fill(z,sqrt(x*x+y*y));
                     h_r->Fill(sqrt(x*x+y*y));
                  }
               }
            }
         }


         event->clear(); 
         generatorEvent->clear();
         geantEvent->clear();
      }

   cout << endl;
   cout << nevent << " events and " << nb << " bytes read " << endl;  
   cout << endl;

   fout->Write();
   fout->Close();

   f->Close();  
}















         /*
         cout << "    vertices:" << endl;
         for (Int_t ig=0; ig<mc_ngv; ig++) {
            XANAGeantVertex *geantVertex = (XANAGeantVertex*)geantEvent->getGeantVertices()->At(ig);
            cout << "      vertex position  = " << geantVertex->getPosition() << endl;
            Int_t mc_not = geantVertex->getNumberOfTracks();
            cout << "      vertex tracks: " << mc_not << endl;
            for (Int_t ig=0; ig<mc_not; ig++) {
               XANAGeantTrack *geantTrack = (XANAGeantTrack*)geantVertex->getTracks()->At(ig);
//             cout << "        particle Id                = " << geantTrack->getPartId() << "   (" << ig << ")" << endl;
//             cout << "        momentum                   = " << geantTrack->getMomentum() << "   (" << ig << ")" << endl;
               XANAGeneratorParticle *genPart = geantTrack->getGeneratorParticle();
               cout << geantTrack << " " << genPart << endl;
               cout << " genPart -- " << genPart->getPartId() << "  " << genPart->getMo1() << endl;
            }
         }
         */

---