1 /********************************************************************
   2 
   3 NAME:     PdfTreeMaker.cxx
   4 PACKAGE:  offline/PhysicsAnalysis/TopPhys/SingleTopDPDMaker/
   5 
   6 AUTHOR:   Clemens Lange <clemens.lange@SPAMNOTdesy.de>
   7 CREATED:  October 2009
   8 
   9 EDITED :  --
  10 
  11 PURPOSE:  SingleTopDPDMaker Tool that inserts PdfInfo using two 
  12           different methods for mc08 and mc09.
  13           mc09: uses genEvt->pdf_info()
  14           mc08: as pdf_info not available, fall-back method fills
  15                 first two partons (usually barcode 3 & 4) from 
  16                 TruthParticleContainer
  17 
  18 ********************************************************************/
  19 
  20 #include "SingleTopDPDMaker/PdfTreeMaker.h"
  21 
  22 PdfTreeMaker::PdfTreeMaker(const std::string& name, ISvcLocator* pSvcLocator) : Algorithm(name, pSvcLocator)
  23 {
  24   
  25   declareProperty("TreeName", m_TreeName = "PdfTree");
  26   
  27   declareProperty("McEventContainerName",       m_McEventContainerName = "GEN_AOD"      );
  28   declareProperty("TruthParticleContainerName", m_TruthParticleContainerName = "SpclMC" );
  29   declareProperty("PdfMethod",                  m_PdfMethod = "mc09"                    );
  30   //mc09: get genEvt->pdf_info, mc08: try to obtain pdf information from TruthParticleContainer
  31   declareProperty("Mc08_barcode",               m_Mc08_barcode                          );
  32   //initial partons usually found at barcodes 3 and 4 (1 and 2 are incoming protons) 
  33   declareProperty("Mc08_pdgId",                 m_Mc08_pdgId                            );
  34   // all quarks (but top) and gluons
  35   
  36 }
  37 
  38 PdfTreeMaker::~PdfTreeMaker(){}
  39 
  40 StatusCode PdfTreeMaker::initialize()
  41 {
  42 
  43   MsgStream MsgLog(msgSvc(), name());
  44   
  45   MsgLog << MSG::INFO << "Initializing " << name() << endreq;
  46     
  47   MsgLog << MSG::INFO << "Tree Name = " << m_TreeName << endreq;
  48   
  49   MsgLog << MSG::INFO << "MC Event Container Name = " << m_McEventContainerName << endreq;
  50   MsgLog << MSG::INFO << "TruthParticle Container Name = " << m_TruthParticleContainerName << endreq;
  51   if (m_PdfMethod.find( "mc08", 0 ) != std::string::npos) {
  52     if (m_Mc08_barcode.empty()) {
  53       MsgLog << MSG::INFO << "Mc08_barcode not set. Using standard: 1,2,3,4 " << endreq;
  54       m_set_Mc08_barcode.insert(1); m_set_Mc08_barcode.insert(2); m_set_Mc08_barcode.insert(3); m_set_Mc08_barcode.insert(4);
  55     }
  56     else m_set_Mc08_barcode.insert(m_Mc08_barcode.begin(), m_Mc08_barcode.end());
  57     if (m_Mc08_pdgId.empty()) {
  58       MsgLog << MSG::INFO << "Mc08_pdgId not set. Using standard: 1,2,3,4,5,21 " << endreq;
  59       m_set_Mc08_pdgId.insert(1); m_set_Mc08_pdgId.insert(2); m_set_Mc08_pdgId.insert(3); m_set_Mc08_pdgId.insert(4); m_set_Mc08_pdgId.insert(5); m_set_Mc08_pdgId.insert(21);
  60     }
  61     else m_set_Mc08_pdgId.insert(m_Mc08_pdgId.begin(), m_Mc08_pdgId.end());
  62   }
  63 
  64   StatusCode status;
  65   
  66   status = service("StoreGateSvc", m_StoreGateSvc);
  67 
  68   if(status.isFailure()) 
  69   {
  70     MsgLog << MSG::FATAL << "StoreGateSvc Unavailable" << endreq;
  71     return StatusCode::FAILURE;
  72   }
  73   
  74   status = service("THistSvc", m_THistSvc);
  75 
  76   if(status.isFailure()) 
  77   {
  78     MsgLog << MSG::FATAL << "THistSvc Unavailable" << endreq;
  79     return StatusCode::FAILURE;
  80   }
  81   
  82   std::string TreeTitle = m_TreeName;
  83   m_Tree = new TTree(m_TreeName.c_str(), TreeTitle.c_str());
  84   
  85   std::string FullTreeName = "/SingleTop/"+m_TreeName;
  86   status = m_THistSvc->regTree(FullTreeName, m_Tree);
  87 
  88   if(status.isFailure()) 
  89   {
  90     MsgLog << MSG::FATAL << "regTree Unavailable" << endreq;
  91     return StatusCode::FAILURE;
  92   }
  93             
  94   m_InputEventNumber = 0;
  95   
  96   this->SetBranches();
  97     
  98   return StatusCode::SUCCESS;
  99   
 100 }
 101 
 102 StatusCode PdfTreeMaker::finalize()
 103 {
 104   
 105   MsgStream MsgLog(msgSvc(), name());
 106   
 107   MsgLog << MSG::INFO << "Finalizing " << name() << endreq;
 108   
 109   MsgLog << MSG::INFO << "Event Number = " << m_InputEventNumber << endreq;
 110   
 111   this->DeleteBranches();
 112       
 113   return StatusCode::SUCCESS;
 114   
 115 }
 116 
 117 StatusCode PdfTreeMaker::execute()
 118 {
 119 
 120   MsgStream MsgLog(msgSvc(), name());
 121   
 122   StatusCode status;
 123 
 124   m_InputEventNumber++;
 125   
 126   status = m_StoreGateSvc->retrieve(m_McEventContainer, m_McEventContainerName);
 127   
 128   if(status.isFailure())
 129   {
 130     m_McEventContainer = 0;
 131     if(m_InputEventNumber == 1) MsgLog << MSG::INFO << "No MC Event Container" << endreq;
 132   }
 133   
 134   status = m_StoreGateSvc->retrieve(m_TruthParticleContainer, m_TruthParticleContainerName);
 135   
 136   if(status.isFailure())
 137   {
 138     m_TruthParticleContainer = 0;
 139     if(m_InputEventNumber == 1) MsgLog << MSG::FATAL << "No Truth Particle Container" << endreq;
 140   }
 141                           
 142   this->ClearBranches();
 143   this->FillBranches();
 144           
 145   return StatusCode::SUCCESS;
 146   
 147 }
 148 
 149 void PdfTreeMaker::SetBranches()
 150 {
 151 
 152   m_Tree->Branch("id1",      &m_id1,      "id1/I"     );
 153   m_Tree->Branch("id2",      &m_id2,      "id2/I"     ); 
 154   m_Tree->Branch("x1",       &m_x1,       "x1/D"      ); 
 155   m_Tree->Branch("x2",       &m_x2,       "x2/D"      ); 
 156   m_Tree->Branch("scalePDF", &m_scalePDF, "scalePDF/D"); 
 157 
 158   m_Tree->Branch("InitialParton_px",   &m_InitialParton_px);
 159   m_Tree->Branch("InitialParton_py",   &m_InitialParton_py);
 160   m_Tree->Branch("InitialParton_pz",   &m_InitialParton_pz);
 161   m_Tree->Branch("InitialParton_pt",   &m_InitialParton_pt);
 162   m_Tree->Branch("InitialParton_e",   &m_InitialParton_e);
 163   m_Tree->Branch("InitialParton_pdgId",   &m_InitialParton_pdgId);
 164   m_Tree->Branch("InitialParton_barcode",   &m_InitialParton_barcode);
 165   m_Tree->Branch("InitialParton_N",   &m_InitialParton_N,   "InitialParton_N/I" );
 166   
 167   
 168 }
 169 
 170 void PdfTreeMaker::ClearBranches()
 171 {
 172 
 173   m_id1      = 0;
 174   m_id2      = 0;
 175   m_x1       = 0.0;
 176   m_x2       = 0.0;
 177   m_scalePDF = 0.0;
 178   
 179   m_InitialParton_px->clear();
 180   m_InitialParton_py->clear();
 181   m_InitialParton_pz->clear();
 182   m_InitialParton_pt->clear();
 183   m_InitialParton_e->clear();
 184   m_InitialParton_pdgId->clear();
 185   m_InitialParton_barcode->clear();
 186   m_InitialParton_N = 0;
 187 
 188 }
 189 
 190 void PdfTreeMaker::FillBranches()
 191 {
 192 
 193   MsgStream MsgLog(msgSvc(), name());
 194   
 195   // mc09 method
 196   if (m_PdfMethod.find( "mc09", 0 ) != std::string::npos) {
 197     if(m_McEventContainer)
 198     {
 199       if(m_InputEventNumber == 1) MsgLog << MSG::INFO << "Using mc09 Pdf Info method" << endreq;
 200       const HepMC::GenEvent* genEvt = *(m_McEventContainer->begin());
 201       const HepMC::PdfInfo *pdfInformation = genEvt->pdf_info();
 202       if(pdfInformation) {
 203         m_id1 = pdfInformation->id1();
 204         m_id2 = pdfInformation->id2();
 205         m_x1  = pdfInformation->x1();
 206         m_x2  = pdfInformation->x2();
 207         m_scalePDF = pdfInformation->scalePDF();
 208       }
 209       else if(m_InputEventNumber == 1) MsgLog << MSG::WARNING << "No mc09 PdfInfo available, cannot write mc09 PdfInfo!" << endreq;
 210     }
 211     else if(m_InputEventNumber == 1) MsgLog << MSG::WARNING << "No McEventContainer available, cannot write mc09 PdfInfo!" << endreq;
 212   }
 213     
 214   // mc08 method
 215   if (m_PdfMethod.find( "mc08", 0 ) != std::string::npos) {
 216     if(m_InputEventNumber == 1) MsgLog << MSG::INFO << "Using mc08 Pdf TruthParticle method" << endreq;
 217     if (m_TruthParticleContainer) {
 218       TruthParticleContainer::const_iterator Itr;
 219       for(Itr=m_TruthParticleContainer->begin(); Itr!=m_TruthParticleContainer->end(); ++Itr) {
 220         const TruthParticle * aParticle = *Itr;
 221         if ((m_set_Mc08_barcode.find(aParticle->barcode()) != m_set_Mc08_barcode.end()) && (m_set_Mc08_pdgId.find(fabs(aParticle->pdgId())) != m_set_Mc08_pdgId.end())) {
 222           m_InitialParton_px->push_back(aParticle->px());
 223           m_InitialParton_py->push_back(aParticle->py());
 224           m_InitialParton_pz->push_back(aParticle->pz());
 225           m_InitialParton_pt->push_back(aParticle->pt());
 226           m_InitialParton_e->push_back(aParticle->e());
 227           m_InitialParton_pdgId->push_back(aParticle->pdgId());
 228           m_InitialParton_barcode->push_back(aParticle->barcode());
 229           ++m_InitialParton_N;
 230         }
 231       }
 232     }
 233     else if(m_InputEventNumber == 1) MsgLog << MSG::WARNING << "No TruthParticleContainer info available, cannot write mc08 PdfInfo" << endreq;
 234   }
 235 
 236   m_Tree->Fill();
 237 
 238 }
 239 
 240 void PdfTreeMaker::DeleteBranches()
 241 {
 242 
 243   delete m_InitialParton_px;
 244   delete m_InitialParton_py;
 245   delete m_InitialParton_pz;
 246   delete m_InitialParton_pt;
 247   delete m_InitialParton_e;
 248   delete m_InitialParton_pdgId;
 249   delete m_InitialParton_barcode;
 250 
 251 }