Migration of FlowAfterburner to HepMC3

Generators/FlowAfterburner/src/AddFlowByShifting.cxx

 /*
-  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 // File:  Generators/FlowAfterburnber/AddFlowByShifting.cxx
@@ -46,7 +46,6 @@ double AddFlowByShifting::vn_func(double x, void *params)
 double AddFlowByShifting::vn_func_derivative(double x, void *params)
 {
    float *par_float = (float*) params;
-   //float phi_0  = par_float[0];
    float *vn    = par_float+1;
    float *psi_n = vn+6;
    double val=1   +2*(   vn[0]*cos(1*(x-psi_n[0]))/1.0 + vn[1]*cos(2*(x-psi_n[1]))/2.0 +
@@ -239,22 +238,20 @@ StatusCode AddFlowByShifting::execute() {
   for (itr = mcFlowCollptr->begin(); itr!=mcFlowCollptr->end(); ++itr) {
     ATH_MSG_DEBUG("Next event in the bag ...");
 
-    //int g_id = (*itr)->signal_process_id();
-    //GeneratorName_print(g_id);
-    //std::cout << std::endl;
 
 
     auto mainvtx=HepMC::barcode_to_vertex(*itr,-1);
     if(m_flow_fluctuations) Set_EbE_Fluctuation_Multipliers(mainvtx,hijing_pars->get_b());
 
+#ifdef HEPMC3
+    int particles_in_event = (*itr)->particles().size();
+#else
     int particles_in_event = (*itr)->particles_size();
+#endif
     m_particles_processed = 0;
-    for ( HepMC::GenVertex::particle_iterator partit =
-	    (*mainvtx).particles_begin(HepMC::children);
-	  partit != (*mainvtx).particles_end(HepMC::children); partit++ ) {
+    for ( auto parent: *mainvtx) {
 
       // Process particles from main vertex
-      auto parent = (*partit);
       CLHEP::HepLorentzVector momentum(parent->momentum().px(),
 				parent->momentum().py(),
 				parent->momentum().pz(),
@@ -371,11 +368,7 @@ void AddFlowByShifting::MoveDescendantsToParent
       }
 
       // now rotate their associated particles
-      for ( HepMC::GenVertex::particle_iterator descpartit
-	    = descvtx->particles_begin(HepMC::children);
-	  descpartit != descvtx->particles_end(HepMC::children);
-	  ++descpartit ) {
-	auto descpart = (*descpartit);
+      for (auto descpart: *descvtx){
         CLHEP::HepLorentzVector momentum(descpart->momentum().px(),
 				  descpart->momentum().py(),
 				  descpart->momentum().pz(),
@@ -492,11 +485,6 @@ double AddFlowByShifting::AddFlowToParent (HepMC::GenParticlePtr parent, const H
 
     if (iter>=1000) return 0;
 
- /*
-    float val=phi   +2*( v1*sin(1*(phi-m_psi_n[0]))/1.0 + v2*sin(2*(phi-m_psi_n[1]))/2.0 +
-                         v3*sin(3*(phi-m_psi_n[2]))/3.0 + v4*sin(4*(phi-m_psi_n[3]))/4.0 +
-                         v5*sin(5*(phi-m_psi_n[4]))/5.0 + v6*sin(6*(phi-m_psi_n  [5]))/6.0 );
-    std::cout<<phi_0<<"  "<<val<<"  "<<phi_0-val<<std::endl;  // */
     phishift = phi-phi_0;
   }
 
@@ -658,10 +646,7 @@ void AddFlowByShifting::Set_EbE_Fluctuation_Multipliers(HepMC::GenVertexPtr main
     double EbE_Vn[6];
     for(int ihar=0;ihar<6;ihar++){m_EbE_Multiplier_vn[ihar]=1.0;EbE_Vn[ihar]=0.0;}
   
-    for(auto  partit  =(*mainvtx).particles_begin(HepMC::children);
-	      partit !=(*mainvtx).particles_end  (HepMC::children); 
-              partit++ ) {
-       auto parent = (*partit);
+    for(auto  parent: *mainvtx) {
        float eta= parent->momentum().pseudoRapidity();
        float pT = parent->momentum().perp();
 

Generators/FlowAfterburner/src/CheckFlow.cxx

@@ -156,7 +156,6 @@ StatusCode CheckFlow::initialize(){
 
   m_tesIO = new GenAccessIO();
 
-  // return StatusCode::SUCCESS;
   return result;
 }
 
@@ -180,7 +179,6 @@ StatusCode CheckFlow::execute() {
 
 
   if ( m_sgSvc->retrieve(hijing_pars, "Hijing_event_params").isFailure() ) {
-//  if ( evtStore()->retrieve(hijing_pars, "Hijing_event_params").isFailure() ) {
     msg(MSG::ERROR) << "Could not retrieve Hijing_event_params"
 	   << endmsg;
     return StatusCode::FAILURE;
@@ -188,13 +186,6 @@ StatusCode CheckFlow::execute() {
   float b = hijing_pars->get_b();
   float phiR = hijing_pars->get_bphi();
 
-  //msglog << MSG::DEBUG<<"SOUMYA  "
-  // <<hijing_pars->get_psi(1)
-  // <<"   "<<hijing_pars->get_psi(2)
-  // <<"  "<<hijing_pars->get_psi(3)
-  // <<hijing_pars->get_psi(4)
-  // <<"   "<<hijing_pars->get_psi(5)
-  // <<"  "<<hijing_pars->get_psi(6)  << endmsg;
 
   // Check cut on impact parameter b
   if(b<m_bcut_min || b>m_bcut_max) 

Generators/FlowAfterburner/src/CheckFlow_New.cxx

 /*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 // File:  Generators/FlowAfterburnber/CheckFlow_New.h
@@ -254,7 +254,6 @@ StatusCode CheckFlow_New::execute() {
 
 
   if ( m_sgSvc->retrieve(hijing_pars, "Hijing_event_params").isFailure() ) {
-//  if ( evtStore()->retrieve(hijing_pars, "Hijing_event_params").isFailure() ) {
     msg(MSG::ERROR) << "Could not retrieve Hijing_event_params"<< endmsg;
     return StatusCode::FAILURE;
   }
@@ -389,16 +388,16 @@ StatusCode CheckFlow_New::execute() {
     double pt     = pitr->momentum().perp();
     double rapid  = pitr->momentum().pseudoRapidity();
     double phi    = pitr->momentum().phi();
-    if( (fabs(rapid) >= m_rapcut_min) && (fabs(rapid) <= m_rapcut_max) &&
-	(fabs(pt) >= m_ptcut_min) && (fabs(pt) <= m_ptcut_max) ) {
+    if( (std::abs(rapid) >= m_rapcut_min) && (std::abs(rapid) <= m_rapcut_max) &&
+	(std::abs(pt) >= m_ptcut_min) && (std::abs(pt) <= m_ptcut_max) ) {
 
       for(int ihar=0;ihar<6;ihar++){
         float       temp=(ihar+1)*(phi-Psi_n_reco_pos[ihar]);
         if(rapid>0) temp=(ihar+1)*(phi-Psi_n_reco_neg[ihar]);
 
 
-        int ieta= (int)(fabs(rapid)*n_etabin/eta_bin_max);
-        if(ieta>=0 && ieta<n_etabin) m_profile_pt_dep_reco [ihar][ieta]->Fill(pt/1000,cos(temp));
+        int ieta= (int)(std::abs(rapid)*n_etabin/eta_bin_max);
+        if(ieta>=0 && ieta<n_etabin) m_profile_pt_dep_reco [ihar][ieta]->Fill(pt/1000,std::cos(temp));
 
         float temp_pt=pt/1000;
         for(int ipt=0;ipt<n_ptbin;ipt++){
@@ -421,19 +420,6 @@ StatusCode CheckFlow_New::execute() {
 StatusCode CheckFlow_New::finalize() {
   msg(MSG::INFO) << ">>> CheckFlow_New from finalize" << endmsg;
 
-/*
-  for(int ihar=0;ihar<6;ihar++){
-    double reso=m_profile_resolution->GetBinContent(ihar+1);
-    if (reso >=0) reso= sqrt( reso);
-    else          reso=-sqrt(-reso);
-    for(int ieta=0;ieta<n_etabin;ieta++){
-      m_profile_pt_dep_reco [ihar][ieta]->Scale(1.0/reso);
-    }
-    for(int ipt=0;ipt<n_ptbin;ipt++){
-     m_profile_eta_dep_reco[ihar][ipt]->Scale(1.0/reso);
-    }
-  }
-*/
   return StatusCode::SUCCESS;
 }
 

Generators/FlowAfterburner/src/CheckFlow_New_Minbias.cxx

 /*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 // File:  Generators/FlowAfterburnber/CheckFlow_New_Minbias.h
@@ -314,8 +314,7 @@ StatusCode CheckFlow_New_Minbias::execute() {
     }
 
     //EbE vn for ID acceptance end pt>0.5GeV
-    //if(fabs(rapid)<=2.5 &&fabs(pt)>=500){
-    if(fabs(pt)>=500){
+    if(std::abs(pt)>=500){
       tot_ID1++;
       for(int ihar=0;ihar<6;ihar++){
         cos_ID1[ihar]+=cos((ihar+1)*phi);