sPhenixAj/src/PHAJMaker.cc

0001 //----------------------------------------------------------
0002 // Module for flow jet reconstruction in sPHENIX
0003 // J. Orjuela-Koop
0004 // May 5 2015
0005 //----------------------------------------------------------
0006
0007 #include"PHAJMaker.h"
0008 #include"PseudoJetPlus.h"
0009
0010 #include<calobase/RawCluster.h>
0011 #include<calobase/RawClusterContainer.h>
0012
0013 #include <g4hough/SvtxTrackMap.h>
0014 #include <g4hough/SvtxTrack.h>
0015
0016 #include<phool/PHCompositeNode.h>
0017 #include<phool/PHNodeIterator.h>
0018 #include<phool/PHNodeReset.h>
0019 #include<fun4all/getClass.h>
0020
0021 #include <fastjet/JetDefinition.hh>
0022 #include <fastjet/PseudoJet.hh>
0023 #include <fastjet/ClusterSequence.hh>
0024 #include <fastjet/SISConePlugin.hh>
0025
0026 #include <TF1.h>
0027 #include <TFile.h>
0028 #include <TH1.h>
0029 #include <TH2.h>
0030 #include <TLorentzVector.h>
0031 #include <TNtuple.h>
0032
0033 #include <cmath>
0034
0035 using namespace std;
0036
0037 typedef std::map<int,TLorentzVector*> tlvmap;
0038
0039 /*
0040  * Constructor
0041  */
0042
0043 // Standard ctor
0044 PHAJMaker::PHAJMaker( const std::string &name,
0045               double R,
0046               double jet_ptmin, double jet_ptmax,
0047               double LeadPtMin, double SubLeadPtMin,
0048               double max_track_rap, double PtConsLo, double PtConsHi,
0049               double dPhiCut
0050               )
0051   :   SubsysReco(name),
0052       R(R),
0053       jet_ptmin(jet_ptmin), jet_ptmax(jet_ptmax),
0054       LeadPtMin(LeadPtMin), SubLeadPtMin(SubLeadPtMin),
0055       max_track_rap (max_track_rap), PtConsLo (PtConsLo), PtConsHi (PtConsHi),
0056       dPhiCut (dPhiCut),
0057       pJAhi (0), pJAlo(0), pOtherJAlo(0)
0058 {
0059   // derived rapidity cuts
0060   // ---------------------
0061   max_rap      = max_track_rap-R;
0062   ghost_maxrap = max_rap + 2.0 * R;
0063
0064   // Constituent selectors
0065   // ---------------------
0066   select_track_rap = fastjet::SelectorAbsRapMax(max_track_rap);
0067   select_lopt      = fastjet::SelectorPtMin( PtConsLo );
0068   select_hipt      = fastjet::SelectorPtMin( PtConsHi );
0069
0070   // Provide but turn off charge selector
0071   // fastjet::Selector select_track_charge= SelectorChargeRange( -3, 3);
0072   fastjet::Selector select_track_charge= fastjet::SelectorIdentity();
0073   slo     = select_track_rap * select_lopt * select_track_charge;
0074   shi     = select_track_rap * select_hipt * select_track_charge;
0075
0076   // Jet candidate selectors
0077   // -----------------------
0078   select_jet_rap     = fastjet::SelectorAbsRapMax(max_rap);
0079   select_jet_pt_min  = fastjet::SelectorPtMin( jet_ptmin );
0080   select_jet_pt_max  = fastjet::SelectorPtMax( jet_ptmax );
0081   sjet = select_jet_rap && select_jet_pt_min && select_jet_pt_max;
0082
0083   // Choose a jet and area definition
0084   // --------------------------------
0085   jet_def = fastjet::JetDefinition(fastjet::antikt_algorithm, R);
0086
0087   // create an area definition for the clustering
0088   //----------------------------------------------------------
0089   // ghosts should go up to the acceptance of the detector or
0090   // (with infinite acceptance) at least 2R beyond the region
0091   // where you plan to investigate jets.
0092   area_spec = fastjet::GhostedAreaSpec( ghost_maxrap, AjParameters::ghost_repeat, AjParameters::ghost_area );
0093   // // DEBUG
0094   // // area_spec.set_grid_scatter(0);
0095   // // area_spec.set_pt_scatter(0);
0096   area_def = fastjet::AreaDefinition(fastjet::active_area_explicit_ghosts, area_spec);
0097
0098   // DEBUG
0099   // area_spec = fastjet::VoronoiAreaSpec( 0.9 );
0100   // area_def = fastjet::AreaDefinition( fastjet::VoronoiAreaSpec( 0.9 ) );
0101
0102   std::cout << " ################################################### " << std::endl;
0103   std::cout << "Leading jet above " << LeadPtMin << std::endl;
0104   std::cout << "Sub-Leading jet above " << SubLeadPtMin << std::endl;
0105   std::cout << "Clustered with " << jet_def.description() << std::endl;
0106   std::cout << "Area Spec " << area_spec.description() << std::endl;
0107   std::cout << "Area Def  " << area_def.description() << std::endl;
0108   std::cout << " ################################################### " << std::endl;
0109
0110   // std::cout << slo.applies_jet_by_jet() << std::endl;
0111   // std::cout << shi.applies_jet_by_jet() << std::endl;
0112   // std::cout << sjet.applies_jet_by_jet() << std::endl;
0113   // std::cout << " ################################################### " << std::endl;
0114
0115 }
0116
0117 /*
0118  * Empty destructor
0119  */
0120 PHAJMaker::~PHAJMaker()
0121 {
0122 }
0123
0124 /*
0125  * Initialize module
0126  */
0127 int PHAJMaker::Init(PHCompositeNode* topNode)
0128 {
0129   cout << "------PHAJMaker::Init(PHCompositeNode*)------" << endl;
0130
0131   // Histo Manager
0132   // -------------
0133   MyHistos = new Fun4AllHistoManager("MyHistos");
0134   MyHistos->setOutfileName( "AjHistos.root");
0135   Fun4AllServer *se = Fun4AllServer::instance();
0136   se->registerHistoManager(MyHistos);
0137
0138   // Histos
0139   // ------
0140   TH1::SetDefaultSumw2(true);
0141   TH2::SetDefaultSumw2(true);
0142   TH3::SetDefaultSumw2(true);
0143
0144   UnmatchedAJ_hi = new TH2D( "UnmatchedAJ_hi","Unmatched A_{J} for hard constituent jets;A_{J};Refmult;fraction", 40, -0.3, 0.9, 800, -0.5, 799.5 );
0145   AJ_hi = new TH2D( "AJ_hi","A_{J} for hard constituent jets;A_{J};Refmult;fraction", 40, -0.3, 0.9, 800, -0.5, 799.5 );
0146   AJ_lo = new TH2D( "AJ_lo","A_{J} for soft constituent jets;A_{J};Refmult;fraction", 40, -0.3, 0.9, 800, -0.5, 799.5 );
0147
0148   MyHistos->registerHisto(UnmatchedAJ_hi);
0149   MyHistos->registerHisto(AJ_hi);
0150   MyHistos->registerHisto(AJ_lo);
0151
0152   return 0;
0153 }
0154
0155 /*
0156  * Process event
0157  */
0158 int PHAJMaker::process_event(PHCompositeNode* topNode)
0159 {
0160   cout << "------PHAJMaker::process_event(PHCompositeNode*)------" << endl;
0161
0162   //-------------------------------------------------
0163   // Get Information from Node Tree
0164   //-------------------------------------------------
0165   PJContainer* particlesCont = findNode::getClass<PJContainer> (topNode,"PJ");
0166   std::vector<fastjet::PseudoJet>& particles = particlesCont->data;
0167
0168   // Classifier, such as Centrality.
0169   double EventClassifier=0;
0170
0171   // We want to hold onto the jetanalyzer objects, so they're created dynamically
0172   // Need to delete them by hand
0173   if (pJAhi){    delete pJAhi;    pJAhi=0;  }
0174   if (pJAlo){    delete pJAlo;    pJAlo=0;  }
0175   if ( pOtherJAlo ){    delete pOtherJAlo; pOtherJAlo=0;  }
0176
0177
0178   DiJetsHi.clear();
0179   DiJetsLo.clear();
0180   OtherDiJetsLo.clear();
0181
0182   pLo.clear();
0183   pHi.clear();
0184
0185   Has10Gev=false;
0186
0187   // Select particles to perform analysis on
0188   // ---------------------------------------
0189   pLo = slo( particles );
0190   pHi = shi( particles );
0191
0192   // Background selector
0193   // -------------------
0194   // It is unclear to me why, but it leads to segfaults if only a once-initialized member :-/
0195   fastjet::Selector selector_bkgd = fastjet::SelectorAbsRapMax( max_rap ) * (!fastjet::SelectorNHardest(2));
0196   // selector_bkgd=fastjet::SelectorAbsRapMax( max_rap );
0197
0198   // find high constituent pT jets
0199   // -----------------------------
0200   pJAhi = new JetAnalyzer( pHi, jet_def ); // NO background subtraction
0201   JetAnalyzer& JAhi = *pJAhi;
0202   // std::cout << 0 << std::endl;
0203   JAhiResult = fastjet::sorted_by_pt( sjet ( JAhi.inclusive_jets() ) );
0204   // std::cout << 1 << std::endl;
0205
0206
0207   // DEBUG
0208   // -----
0209   cout << " -------------------------------------------------------------------- " << endl;
0210   cout << "particles.size() = " << particles.size() << endl;
0211   cout << "pLo.size() = " << pLo.size() << endl;
0212   cout << "pHi.size() = " << pHi.size() << endl;
0213   cout << "JAhiResult.size() = " << JAhiResult.size() << endl;
0214   for ( unsigned int i=0; i<JAhiResult.size() ; ++i ){
0215     cout << "JAhiResult.at(" << i<< ").pT() = " << JAhiResult.at(i).pt() << endl;
0216   }
0217   cout << " -------------------------------------------------------------------- " << endl;
0218
0219
0220   if ( JAhiResult.size() < 1 )                 {     return 0; }
0221   if ( JAhiResult.at(0).pt() > 10 )            { Has10Gev=true; }
0222
0223   if ( JAhiResult.size() < 2 )                 {     return 0; }
0224   if ( JAhiResult.at(0).pt() < LeadPtMin )     {     return 0; }
0225   if ( JAhiResult.at(1).pt() < SubLeadPtMin )  {     return 0; }
0226
0227   // back to back? Answer this question with a selector
0228   // ---------------------------------------------------
0229   DiJetsHi = SelectorDijets( dPhiCut ) ( JAhiResult );
0230   if ( DiJetsHi.size() == 0 ) {
0231     // std::cout << " NO dijet found" << std::endl;
0232     return 0;
0233   }
0234   assert ( DiJetsHi.size() == 2 && "SelectorDijets returned impossible number of Dijets." );
0235
0236   // // FOR EMBEDDING: At least one matched to the reference jet?
0237   // // ---------------------------------------------------------
0238   // if ( ToMatch ){
0239   //   if ( !IsMatched( DiJetsHi, *ToMatch, R ) ) return 0;
0240   // }
0241
0242   // Calculate Aj and fill histos -- for unmatched high constituent pT di-jets
0243   // -------------------------------------------------------------------------
0244   UnmatchedAJ_hi->Fill ( CalcAj( DiJetsHi ), EventClassifier );
0245
0246   // find corresponding jets with soft constituents
0247   // ----------------------------------------------
0248   pJAlo = new JetAnalyzer( pLo, jet_def, area_def, selector_bkgd ); // WITH background subtraction
0249   JetAnalyzer& JAlo = *pJAlo;
0250   fastjet::Subtractor* BackgroundSubtractor =  JAlo.GetBackgroundSubtractor();
0251   JAloResult = fastjet::sorted_by_pt( (*BackgroundSubtractor)( JAlo.inclusive_jets() ) );
0252
0253   // Using selectors mostly because I can :)
0254   fastjet::Selector SelectClose = fastjet::SelectorCircle( R );
0255
0256   // Leading:
0257   SelectClose.set_reference( DiJetsHi.at(0) );
0258   // Pulling apart this one-liner,   DiJetsLo.push_back ( ((SelectClose && fastjet::SelectorNHardest(1)) ( JAloResult )).at(0));
0259   // More readable and we can catch abnormalities
0260   std::vector<fastjet::PseudoJet> MatchedToLead = sorted_by_pt(SelectClose( JAloResult ));
0261   if ( MatchedToLead.size() == 0 ) {
0262     std::cerr << "PROBLEM: SelectorClose returned no match to leading jet." << std::endl;
0263     // return 2;
0264     return 0;
0265   }
0266   DiJetsLo.push_back ( MatchedToLead.at(0) );
0267
0268   SelectClose.set_reference( DiJetsHi.at(1) );
0269   std::vector<fastjet::PseudoJet> MatchedToSubLead = sorted_by_pt(SelectClose( JAloResult ));
0270   if ( MatchedToSubLead.size() == 0 ) {
0271     std::cerr << "PROBLEM: SelectorClose returned no match to sub-leading jet." << std::endl;
0272     //     return 2;
0273     return 0;
0274   }
0275   DiJetsLo.push_back ( MatchedToSubLead.at(0) );
0276
0277   if ( fabs(DiJetsLo.at(0).eta())>max_rap ) std:: cerr << "Uh-oh... Lead jet eta = " << DiJetsLo.at(0).eta() << std::endl;
0278   if ( fabs(DiJetsLo.at(1).eta())>max_rap ) std:: cerr << "Uh-oh... SubLead jet eta = " << DiJetsLo.at(1).eta() << std::endl;
0279
0280   // And we're done! Calculate A_J
0281   // -----------------------------
0282   AJ_hi->Fill( CalcAj( DiJetsHi ), EventClassifier );
0283   AJ_lo->Fill( CalcAj( DiJetsLo ), EventClassifier );
0284
0285   // DEBUG
0286   // -----
0287   cout << " -------------------------------------------------------------------- " << endl;
0288   cout << "DiJetsHi.at(0).pT() = " << DiJetsHi.at(0).pt() << endl;
0289   cout << "DiJetsHi.at(1).pT() = " << DiJetsHi.at(1).pt() << endl;
0290   cout << "AJ, hi: " << CalcAj( DiJetsHi ) << endl;
0291   cout << "DiJetsLo.at(0).pT() = " << DiJetsLo.at(0).pt() << endl;
0292   cout << "DiJetsLo.at(1).pT() = " << DiJetsLo.at(1).pt() << endl;
0293   cout << "AJ, lo: " << CalcAj( DiJetsLo ) << endl;
0294   cout << " -------------------------------------------------------------------- " << endl;
0295
0296
0297
0298   return 0;
0299 }
0300
0301 /*
0302  * Write jets to node tree
0303  */
0304 /*
0305 todo - we are working on a jet object
0306 void PHAJMaker::save_jets_to_nodetree()
0307 {
0308   flow_jet_container = new PHPyJetContainerV2();
0309 }
0310 */