ATLAS Offline Software
CheckFlow.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // File: Generators/FlowAfterburnber/CheckFlow.h
6 // Description:
7 // This is a simple algorithm to histogram particle properties
8 // for diagnosing of flow generation
9 //
10 // AuthorList:
11 // Andrzej Olszewski: Initial Code February 2006
12 // Andrzej Olszewski: Converted to ROOT histograms July 2007
13 #include "FlowAfterburner/CheckFlow.h"
14 #include "GeneratorObjects/McEventCollection.h"
15 
16 #include <TH1F.h>
17 #include <TH2F.h>
18 #include <TH3F.h>
19 
20 #include "AtlasHepMC/GenEvent.h"
21 #include "AtlasHepMC/GenParticle.h"
22 #include "AtlasHepMC/GenVertex.h"
23 
24 #include "GeneratorObjects/HijingEventParams.h"
25 
26 
27 typedef std::vector<HepMC::ConstGenParticlePtr> MCparticleCollection ;
28 
29 CheckFlow::CheckFlow(const std::string& name, ISvcLocator* pSvcLocator) :
30  AthAlgorithm(name, pSvcLocator),
31  m_hgenerated ( NULL ),
32  m_b ( NULL ),
33  m_phi ( NULL ),
34  m_phiR ( NULL ),
35  m_phi_vs_phiR ( NULL ),
36  m_phiv1reco_vs_phiR ( NULL ),
37  m_phiv2reco_vs_phiR ( NULL ),
38  m_phi_vs_phiR_etap ( NULL ),
39  m_phi_vs_phiR_etan ( NULL ),
40  m_v2betapth ( NULL ),
41  m_ebetapth ( NULL ),
42  m_sgSvc ( NULL ),
43  m_tesIO ( NULL )
44 {
45  //Declare the algorithm's properties
46  declareProperty("McEventKey", m_key="FLOW_EVENT");
47  declareProperty("HistogramFlag", m_produceHistogram = true );
48  declareProperty("ImpactCutMin", m_bcut_min = 0 );
49  declareProperty("ImpactCutMax", m_bcut_max = 99 );
50  declareProperty("PtCutMin", m_ptcut_min = 0 );
51  declareProperty("PtCutMax", m_ptcut_max = 999999 );
52  declareProperty("RapidityCutMin", m_rapcut_min = 0 );
53  declareProperty("RapidityCutMax", m_rapcut_max = 5.5 );
54 }
55 
56 StatusCode CheckFlow::initialize(){
57  StatusCode result = StatusCode::SUCCESS;
58  msg(MSG::INFO) << ">>> CheckFlow from Initialize" << endmsg;
59 
60  StatusCode sc = service("StoreGateSvc", m_sgSvc);
61  if (sc.isFailure()) {
62  msg(MSG::ERROR) << "Could not find StoreGateSvc" << endmsg;
63  return sc;
64  }
65 
66  m_hgenerated = new TH1F("ngen","Generated",100,0,100000);
67  m_b = new TH1F("b","Impact parameter",35,0.,35.0);
68  m_phi = new TH1F("phi","Phi",100,-M_PI,M_PI);
69  m_phiR = new TH1F("phiR","PhiR",100,0,2*M_PI);
70  m_phi_vs_phiR = new TH1F("phi_vs_phiR","Phi - PhiR",100,-M_PI,M_PI);
71  m_phiv1reco_vs_phiR = new TH2F("phiv1reco_vs_phiR",
72  "PhiV1Reco vs PhiR", 30,-M_PI,M_PI,30,-M_PI,M_PI);
73  m_phiv2reco_vs_phiR = new TH2F("phiv2reco_vs_phiR",
74  "PhiV2Reco vs PhiR", 30,-M_PI/2,M_PI/2,30,-M_PI/2,M_PI/2);
75  m_phi_vs_phiR_etap = new TH1F("phi_vs_phiR_etap",
76  "Phi - PhiR positive eta",
77  100,-M_PI,M_PI);
78  m_phi_vs_phiR_etan = new TH1F("phi_vs_phiR_etan",
79  "Phi - PhiR negative eta",
80  100,-M_PI,M_PI);
81  m_v2betapth = new TH3F("v2betapth",
82  "V2 vs b, eta, pt",
83  20,0,20, 30,-7.5,7.5, 25,0,5000);
84  m_ebetapth = new TH3F("ebetapth",
85  "Tracks vs b, eta, pt",
86  20,0,20, 30,-7.5,7.5, 25,0,5000);
87 
88  ITHistSvc *rootHistSvc;
89  if (!service("THistSvc", rootHistSvc, true).isSuccess()) {
90  msg(MSG::ERROR) << "Unable to locate THistSvc" << endmsg;
91  return StatusCode::FAILURE;
92  }
93 
94  std::string StreamAndPath="/FlowOutPut/";
95  std::string histPath = StreamAndPath;
96  if ( rootHistSvc->regHist(histPath+m_hgenerated->GetName(),
97  m_hgenerated).isFailure() )
98  msg(MSG::WARNING) << "Can't book "
99  << histPath+m_hgenerated->GetName() << endmsg;
100 
101  if ( rootHistSvc->regHist(histPath+m_b->GetName(),
102  m_b).isFailure() )
103  msg(MSG::WARNING) << "Can't book "
104  << histPath+m_b->GetName() << endmsg;
105 
106  if ( rootHistSvc->regHist(histPath+m_phi->GetName(),
107  m_phi).isFailure() )
108  msg(MSG::WARNING) << "Can't book "
109  << histPath+m_phi->GetName() << endmsg;
110 
111  if ( rootHistSvc->regHist(histPath+m_phiR->GetName(),
112  m_phiR).isFailure() )
113  msg(MSG::WARNING) << "Can't book "
114  << histPath+m_phiR->GetName() << endmsg;
115 
116  if ( rootHistSvc->regHist(histPath+m_phi_vs_phiR->GetName(),
117  m_phi_vs_phiR).isFailure() )
118  msg(MSG::WARNING) << "Can't book "
119  << histPath+m_phi_vs_phiR->GetName() << endmsg;
120 
121  if ( rootHistSvc->regHist(histPath+m_phiv1reco_vs_phiR->GetName(),
122  m_phiv1reco_vs_phiR).isFailure() )
123  msg(MSG::WARNING) << "Can't book "
124  << histPath+m_phiv1reco_vs_phiR->GetName() << endmsg;
125 
126  if ( rootHistSvc->regHist(histPath+m_phiv2reco_vs_phiR->GetName(),
127  m_phiv2reco_vs_phiR).isFailure() )
128  msg(MSG::WARNING) << "Can't book "
129  << histPath+m_phiv2reco_vs_phiR->GetName() << endmsg;
130 
131  if ( rootHistSvc->regHist(histPath+m_phi_vs_phiR_etap->GetName(),
132  m_phi_vs_phiR_etap).isFailure() )
133  msg(MSG::WARNING) << "Can't book "
134  << histPath+m_phi_vs_phiR_etap->GetName() << endmsg;
135 
136  if ( rootHistSvc->regHist(histPath+m_phi_vs_phiR_etan->GetName(),
137  m_phi_vs_phiR_etan).isFailure() )
138  msg(MSG::WARNING) << "Can't book "
139  << histPath+m_phi_vs_phiR_etan->GetName() << endmsg;
140 
141  if ( rootHistSvc->regHist(histPath+m_v2betapth->GetName(),
142  m_v2betapth).isFailure() )
143  msg(MSG::WARNING) << "Can't book "
144  << histPath+m_v2betapth->GetName() << endmsg;
145 
146  if ( rootHistSvc->regHist(histPath+m_ebetapth->GetName(),
147  m_ebetapth).isFailure() )
148  msg(MSG::WARNING) << "Can't book "
149  << histPath+m_ebetapth->GetName() << endmsg;
150 
151  msg(MSG::DEBUG) << "Histograms have been booked " << endmsg;
152 
153  m_tesIO = new TruthHelper::GenAccessIO();
154 
155  return result;
156 }
157 
158 StatusCode CheckFlow::execute() {
159  msg(MSG::INFO) << ">>> CheckFlow from execute" << endmsg;
160 
161  //
162  // Event parameters
163  //
164 
165 
166 
167 //---------------------------------------------------------------------------------------------------
168  const HijingEventParams *hijing_pars;
169  //HijingEventParams *hijing_pars;
170 //---------------------------------------------------------------------------------------------------
171 
172 
173 
174 
175 
176 
177  if ( m_sgSvc->retrieve(hijing_pars, "Hijing_event_params").isFailure() ) {
178  msg(MSG::ERROR) << "Could not retrieve Hijing_event_params"
179  << endmsg;
180  return StatusCode::FAILURE;
181  }
182  float b = hijing_pars->get_b();
183  float phiR = hijing_pars->get_bphi();
184 
185 
186  // Check cut on impact parameter b
187  if(b<m_bcut_min || b>m_bcut_max)
188  return StatusCode::SUCCESS;
189 
190  m_b->Fill(b, 1.);
191  m_phiR->Fill(phiR, 1.);
192 
193  float ngenerated = 0;
194  double phiv1_recon, phiv1_recop;
195  double phiv1_reco = 0, phiv2_reco = 0;
196  double phi_reco_sin1phip = 0, phi_reco_cos1phip = 0;
197  double phi_reco_sin1phin = 0, phi_reco_cos1phin = 0;
198  double phi_reco_sin2phi = 0, phi_reco_cos2phi = 0;
199 
200  std::vector<HepMC::ConstGenParticlePtr> particles;
201  StatusCode stat = m_tesIO->getMC(particles, false, m_key);
202  if (stat.isFailure()) {
203  msg(MSG::ERROR) << "Could not find " << m_key << endmsg;
204  return stat;
205  }
206 
207  for (auto pitr: particles) {
208  int pid = pitr->pdg_id();
209  int p_stat = pitr->status();
210  double pt = pitr->momentum().perp();
211  double rapid = pitr->momentum().pseudoRapidity();
212  double phi = pitr->momentum().phi();
213  msg(MSG::DEBUG)
214  << " PID = " << pid << " Status = " << p_stat
215  << " Eta = " << rapid << " Phi = " << phi
216  << " PhiR = " << phiR << endmsg;
217 
218  if( (std::abs(rapid) >= m_rapcut_min) && (std::abs(rapid) <= m_rapcut_max) &&
219  (std::abs(pt) >= m_ptcut_min) && (std::abs(pt) <= m_ptcut_max) ) {
220  ngenerated++;
221  m_phi->Fill(phi, 1.);
222  double phi_corr = phi - phiR;
223  // v2 vs b,eta,pt histograms
224  m_v2betapth->Fill(b,rapid,pt,cos(2*phi_corr));
225  m_ebetapth->Fill(b,rapid,pt,1);
226  // -----------------
227  // convert to (-pi,pi) range
228  int kpi = (int)(phi_corr/(2*M_PI));
229  phi_corr -= kpi * 2*M_PI;
230  if (phi_corr < -M_PI) phi_corr += 2*M_PI;
231  if (phi_corr > M_PI) phi_corr -= 2*M_PI;
232  // --------------------------------------
233  m_phi_vs_phiR->Fill(phi_corr, 1.);
234  if(rapid>=0) m_phi_vs_phiR_etap->Fill(phi_corr, 1.);
235  else m_phi_vs_phiR_etan->Fill(phi_corr, 1.);
236  // -------------------------------------------------
237  if( rapid >= 0 ) {
238  phi_reco_sin1phip += std::sin(1*phi);
239  phi_reco_cos1phip += std::cos(1*phi);
240  } else {
241  phi_reco_sin1phin += std::sin(1*phi);
242  phi_reco_cos1phin += std::cos(1*phi);
243  }
244  phi_reco_sin2phi += std::sin(2*phi);
245  phi_reco_cos2phi += std::cos(2*phi);
246  }
247  }
248  m_hgenerated->Fill(ngenerated, 1.);
249 
250  // calculate event plane position
251  phiv1_recop = std::atan2( phi_reco_sin1phip,phi_reco_cos1phip );
252  phiv1_recon = std::atan2( phi_reco_sin1phin,phi_reco_cos1phin ) + M_PI;
253  if( phiv1_recon > M_PI ) phiv1_recon -= 2*M_PI;
254  // averaged v1 plane position (in pos and neg eta ranges)
255  phiv1_reco = (phiv1_recop + phiv1_recon)/2;
256  phiv2_reco = 0.5 * std::atan2( phi_reco_sin2phi,phi_reco_cos2phi );
257  msg(MSG::INFO)
258  << " PhiR = " << phiR
259  << " PhiV1Reco = " << phiv1_reco
260  << " PhiV2Reco = " << phiv2_reco << endmsg;
261 
262  // convert phiR (0, 2*pi) to (-pi,pi) range
263  double phiR_v1corr = phiR;
264  if( phiR > M_PI ) phiR_v1corr = phiR_v1corr - 2*M_PI;
265  // convert phiR (0, 2*pi) to (-pi/2,pi/2) range
266  double phiR_v2corr = phiR;
267  if (phiR > M_PI/2) phiR_v2corr -= M_PI;
268  if (phiR > 3*M_PI/2) phiR_v2corr -= 2*M_PI;
269  m_phiv1reco_vs_phiR->Fill(phiR_v1corr, phiv1_reco);
270  m_phiv2reco_vs_phiR->Fill(phiR_v2corr, phiv2_reco);
271 
272  // End of execution for each event
273  return StatusCode::SUCCESS;
274 }
275 
276 StatusCode CheckFlow::finalize() {
277  msg(MSG::INFO) << ">>> CheckFlow from finalize" << endmsg;
278 
279  return StatusCode::SUCCESS;
280 }
281 
CheckFlow.h
CheckFlow::execute
StatusCode execute()
Definition: CheckFlow.cxx:158
MCparticleCollection
std::vector< HepMC::ConstGenParticlePtr > MCparticleCollection
Definition: CheckFlow.cxx:27
GenEvent.h
get_generator_info.result
result
Definition: get_generator_info.py:21
CheckFlow::m_phi
TH1F * m_phi
Definition: CheckFlow.h:50
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:64
HijingEventParams.h
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
CheckFlow::m_bcut_max
double m_bcut_max
Definition: CheckFlow.h:41
GenVertex.h
test_pyathena.pt
pt
Definition: test_pyathena.py:11
M_PI
#define M_PI
Definition: ActiveFraction.h:11
CheckFlow::m_v2betapth
TH3F * m_v2betapth
Definition: CheckFlow.h:57
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
python.TrigEgammaMonitorHelper.TH2F
def TH2F(name, title, nxbins, bins_par2, bins_par3, bins_par4, bins_par5=None, bins_par6=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:45
CheckFlow::m_phi_vs_phiR_etap
TH1F * m_phi_vs_phiR_etap
Definition: CheckFlow.h:55
GenParticle.h
CheckFlow::m_phiv1reco_vs_phiR
TH2F * m_phiv1reco_vs_phiR
Definition: CheckFlow.h:53
CheckFlow::m_ebetapth
TH3F * m_ebetapth
Definition: CheckFlow.h:58
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
StoreGateSvc::retrieve
StatusCode retrieve(const T *&ptr) const
Retrieve the default object into a const T*.
CheckFlow::finalize
StatusCode finalize()
Definition: CheckFlow.cxx:276
CheckFlow::m_ptcut_max
double m_ptcut_max
Definition: CheckFlow.h:43
CheckFlow::m_rapcut_max
double m_rapcut_max
Definition: CheckFlow.h:45
McEventCollection.h
CheckFlow::m_rapcut_min
double m_rapcut_min
Definition: CheckFlow.h:44
CheckFlow::m_ptcut_min
double m_ptcut_min
Definition: CheckFlow.h:42
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
CheckFlow::initialize
StatusCode initialize()
Definition: CheckFlow.cxx:56
HijingEventParams
Definition: HijingEventParams.h:23
CheckFlow::m_phi_vs_phiR
TH1F * m_phi_vs_phiR
Definition: CheckFlow.h:52
CheckFlow.CheckFlow
CheckFlow
Definition: CheckFlow.py:50
ParticleGun_EoverP_Config.pid
pid
Definition: ParticleGun_EoverP_Config.py:62
CheckFlow::m_phiR
TH1F * m_phiR
Definition: CheckFlow.h:51
CheckFlow::m_hgenerated
TH1F * m_hgenerated
Definition: CheckFlow.h:48
TruthHelper::GenAccessIO::getMC
StatusCode getMC(MCParticleCollection &mcParticles, const bool ifgen=false, const std::string &key="GEN_EVENT") const
Definition: Generators/FlowAfterburner/FlowAfterburner/GenAccessIO.h:32
CheckFlow::m_phi_vs_phiR_etan
TH1F * m_phi_vs_phiR_etan
Definition: CheckFlow.h:56
HijingEventParams::get_bphi
float get_bphi() const
Definition: HijingEventParams.h:92
beamspotman.stat
stat
Definition: beamspotman.py:266
AthAlgorithm
Definition: AthAlgorithm.h:47
CheckFlow::m_tesIO
TruthHelper::GenAccessIO * m_tesIO
Definition: CheckFlow.h:61
CheckFlow::m_key
std::string m_key
Definition: CheckFlow.h:37
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:195
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
CheckFlow::m_phiv2reco_vs_phiR
TH2F * m_phiv2reco_vs_phiR
Definition: CheckFlow.h:54
TruthHelper::GenAccessIO
Definition: Generators/FlowAfterburner/FlowAfterburner/GenAccessIO.h:24
TH3F
Definition: rootspy.cxx:495
HijingEventParams::get_b
float get_b() const
Definition: HijingEventParams.h:91
CheckFlow::m_sgSvc
StoreGateSvc * m_sgSvc
Definition: CheckFlow.h:60
DEBUG
#define DEBUG
Definition: page_access.h:11
AthCommonMsg< Algorithm >::msg
MsgStream & msg() const
Definition: AthCommonMsg.h:24
LArG4FSStartPointFilter.particles
list particles
Definition: LArG4FSStartPointFilter.py:84
CheckFlow::m_b
TH1F * m_b
Definition: CheckFlow.h:49
declareProperty
#define declareProperty(n, p, h)
Definition: BaseFakeBkgTool.cxx:15
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
Generated on Fri Jul 5 2024 21:12:22 for ATLAS Offline Software by   doxygen 1.8.18