22 #include "GaudiKernel/SmartDataPtr.h"
23 #include "GaudiKernel/DataSvc.h"
25 #include "GaudiKernel/ITHistSvc.h"
49 float pt_binvals[
n_ptbin+1]={0.0,0.25,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,6.0,8.0,12.0,16.0,20.0,30.0,40.0};
50 float eta_bin_max = 4.0;
53 SmartIF<ITHistSvc> rootHistSvc{Gaudi::svcLocator()->service(
"THistSvc")};
56 return StatusCode::FAILURE;
59 std::string histPath =
"/FlowOutPut/";
61 for (
int ihar=0;ihar<6;ihar++){
64 sprintf(
name,
"hist_Psi_%d_true",ihar+1);
65 sprintf(
name1,
"Truth Psi_{%d} distribution;%dPsi_{%d} Truth;events",ihar+1,ihar+1,ihar+1);
73 sprintf(
name,
"hist_Psi_%d_reco",ihar+1);
74 sprintf(
name1,
"Reconstructed Psi_{%d} distribution;%dPsi_{%d} Reco;events",ihar+1,ihar+1,ihar+1);
83 for (
int ihar2=0;ihar2<6;ihar2++)
85 int ihar_i=ihar*6+ihar2;
87 sprintf(
name,
"hist_Psi_corr_true_%d_%d",ihar+1,ihar2+1);
88 sprintf(
name1,
"true Psi_{%d} -Psi_{%d};%dPsi_{%d} -%dPsi_{%d} ;events",ihar+1,ihar2+1,ihar+1,ihar+1,ihar2+1,ihar2+1);
97 sprintf(
name,
"hist_Psi_corr_reco_%d_%d",ihar+1,ihar2+1);
98 sprintf(
name1,
"reco Psi_{%d} -Psi_{%d};%dPsi_{%d} -%dPsi_{%d} ;events",ihar+1,ihar2+1,ihar+1,ihar+1,ihar2+1,ihar2+1);
110 sprintf(
name,
"hist_v%d_ebe",ihar+1);
111 sprintf(
name1,
"v%d;v%d;events",ihar+1,ihar+1);
119 sprintf(
name ,
"hist_Psi%d_ebe",ihar+1);
120 sprintf(
name1,
"%d#Delta#Psi;%d(#Psi_{reco}-#Psi_{Truth});events",ihar+1,ihar+1);
128 sprintf(
name ,
"hist_Psi%d_ebe_pt",ihar+1);
129 sprintf(
name1,
"%d#Delta#Psi (pT weighted);%d(#Psi_{reco}-#Psi_{Truth});events",ihar+1,ihar+1);
142 for(
int ieta=0;ieta<
n_etabin;ieta++){
143 sprintf(
name ,
"profile_pt_dep_%d_eta%d" ,ihar+1,ieta);
144 sprintf(
name1,
"v%d vs pT (eta%d);pT;v%d",ihar+1,ieta,ihar+1);
153 for(
int ipt=0;ipt<
n_ptbin;ipt++){
154 sprintf(
name ,
"profile_eta_dep_%d_pt%d",ihar+1,ipt);
155 sprintf(
name1,
"v%d vs #eta; (ipt%d)#eta;v%d",ihar+1,ipt,ihar+1);
165 for(
int ieta=0;ieta<
n_etabin;ieta++){
166 sprintf(
name ,
"profile_pt_dep_reco_%d_eta%d",ihar+1,ieta);
167 sprintf(
name1,
"v%d vs pT (eta%d);pT;v%d",ihar+1,ieta,ihar+1);
177 for(
int ipt=0;ipt<
n_ptbin;ipt++){
178 sprintf(
name ,
"profile_eta_dep_reco_%d_pt%d",ihar+1,ipt);
179 sprintf(
name1,
"v%d vs #eta (pt%d);#eta;v%d",ihar+1,ipt,ihar+1);
197 return StatusCode::SUCCESS;
205 float pt_binvals[
n_ptbin+1]={0.0,0.25,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,6.0,8.0,12.0,16.0,20.0,30.0,40.0};
206 float eta_bin_max = 4.0;
209 float b = hijing_pars->get_b();
210 float Psi_n[6],Psi_n_reco[6];
211 float Psi_n_reco_pos[6],Psi_n_reco_neg[6];
212 for(
int ihar=0;ihar<6;ihar++){Psi_n[ihar]=hijing_pars->get_psi(ihar+1);}
213 msg(
MSG::INFO)<<
"SOUMYA "<<hijing_pars->get_psi(1)<<
" "<<hijing_pars->get_psi(2)<<
" "<<hijing_pars->get_psi(3)
214 <<hijing_pars->get_psi(4)<<
" "<<hijing_pars->get_psi(5)<<
" "<<hijing_pars->get_psi(6)<<
" "<<
b <<
endmsg;
218 if(b<m_bcut_min || b>
m_bcut_max)
return StatusCode::SUCCESS;
221 double ngenerated_pos = 0,ngenerated_pt_pos=0;
222 double cos_n_pos[6],sin_n_pos[6],cos_n_pt_pos[6],sin_n_pt_pos[6];
223 double ngenerated_neg = 0,ngenerated_pt_neg=0;
224 double cos_n_neg[6],sin_n_neg[6],cos_n_pt_neg[6],sin_n_pt_neg[6];
225 for(
int ihar=0;ihar<6;ihar++){cos_n_pos[ihar]=0;sin_n_pos[ihar]=0; cos_n_pt_pos[ihar]=0;sin_n_pt_pos[ihar]=0;
226 cos_n_neg[ihar]=0;sin_n_neg[ihar]=0; cos_n_pt_neg[ihar]=0;sin_n_pt_neg[ihar]=0;}
229 std::vector<HepMC::ConstGenParticlePtr>
particles;
231 if (
stat.isFailure()) {
236 int pid = pitr->pdg_id();
237 double pt = pitr->momentum().perp();
238 double rapid = pitr->momentum().pseudoRapidity();
239 double phi = pitr->momentum().phi();
241 <<
" PID = " <<
pid <<
" Status = " << pitr->status()
242 <<
" Eta = " << rapid <<
" Phi = " <<
phi<<
endmsg;
247 for(
int ihar=0;ihar<6;ihar++){
248 float temp=(ihar+1)*(
phi-Psi_n[ihar]);
250 int ieta= (
int)(std::abs(rapid)*
n_etabin/eta_bin_max);
254 float temp_pt=
pt/1000;
255 for(
int ipt=0;ipt<
n_ptbin;ipt++){
256 if(temp_pt<pt_binvals[ipt+1]){
262 if( rapid >3.2 && rapid< 4.9){
269 ngenerated_pt_pos +=
pt;
271 if( rapid <-3.2 && rapid >-4.9){
278 ngenerated_pt_neg +=
pt;
288 float cos_n[6],sin_n[6],cos_n_pt[6],sin_n_pt[6];
289 for(
int ihar=0;ihar<6;ihar++){
290 cos_n[ihar] = ( cos_n_pos[ihar]+ cos_n_neg[ihar] ) / (ngenerated_pos+ngenerated_neg);
291 sin_n[ihar] = ( sin_n_pos[ihar]+ sin_n_neg[ihar] ) / (ngenerated_pos+ngenerated_neg);
293 float psi_reco=std::atan2(sin_n[ihar],cos_n[ihar])/(ihar+1);
295 m_hist_vn_ebe [ihar]->Fill(std::sqrt(cos_n[ihar]*cos_n[ihar] +sin_n[ihar]*sin_n[ihar] ));
297 Psi_n_reco_pos[ihar]=std::atan2(sin_n_pos[ihar],cos_n_pos[ihar])/ (ihar+1);
298 Psi_n_reco_neg[ihar]=std::atan2(sin_n_neg[ihar],cos_n_neg[ihar])/ (ihar+1);
299 Psi_n_reco [ihar]=psi_reco;
302 cos_n_pt[ihar] = ( cos_n_pt_pos[ihar]+ cos_n_pt_neg[ihar] ) / (ngenerated_pt_pos+ngenerated_pt_neg);
303 sin_n_pt[ihar] = ( sin_n_pt_pos[ihar]+ sin_n_pt_neg[ihar] ) / (ngenerated_pt_pos+ngenerated_pt_neg);
305 psi_reco=std::atan2(sin_n_pt[ihar],cos_n_pt[ihar])/(ihar+1);
311 for(
int ihar=0;ihar<6;ihar++){
316 for(
int ihar2=0;ihar2<6;ihar2++){
317 psi1=(ihar+1)*Psi_n[ihar];psi2=(ihar2+1)*Psi_n[ihar2];
320 psi1=(ihar+1)*Psi_n_reco[ihar];psi2=(ihar2+1)*Psi_n_reco[ihar2];
330 for(
int ihar=0;ihar<6;ihar++){
334 double pt = pitr->momentum().perp();
335 double rapid = pitr->momentum().pseudoRapidity();
336 double phi = pitr->momentum().phi();
340 for(
int ihar=0;ihar<6;ihar++){
341 float temp=(ihar+1)*(
phi-Psi_n_reco_pos[ihar]);
342 if(rapid>0) temp=(ihar+1)*(
phi-Psi_n_reco_neg[ihar]);
345 int ieta= (
int)(std::abs(rapid)*
n_etabin/eta_bin_max);
348 float temp_pt=
pt/1000;
349 for(
int ipt=0;ipt<
n_ptbin;ipt++){
350 if(temp_pt<pt_binvals[ipt+1]){
363 return StatusCode::SUCCESS;
369 return StatusCode::SUCCESS;