2424#include " CommonConstants/MathConstants.h"
2525#include " TDatabasePDG.h"
2626
27+ #include < TF1.h>
28+ #include < TRandom2.h>
29+
2730using namespace o2 ;
2831using namespace o2 ::framework;
2932using namespace o2 ::framework::expressions;
@@ -37,6 +40,11 @@ struct PseudorapidityDensity {
3740 {VARIABLE_WIDTH , 0 ., 0.01 , 0.1 , 0.5 , 1 , 5 , 10 , 15 , 20 , 30 , 40 , 50 , 70 , 100 },
3841 " Centrality/multiplicity percentile binning" };
3942
43+ Configurable<std::vector<float >> parRatio{" parRatio" , {0.773047 , 0.00131658 , 0.00232244 , 1.53395e-05 , -2.36494e-05 , -3.71754e-08 , 4.4547e-08 }, " parameters of pol6 of mes /reco tracks versus zvtx" };
44+
45+ TF1 * fRatio = 0 ;
46+ TRandom2* rand = 0 ;
47+
4048 HistogramRegistry registry{
4149 " registry" ,
4250 {
@@ -76,6 +84,12 @@ struct PseudorapidityDensity {
7684 x->SetBinLabel (3 , " Reconstructed" );
7785 x->SetBinLabel (4 , " Selected" );
7886 x->SetBinLabel (5 , " Selected INEL>0" );
87+
88+ fRatio = new TF1 (" fRatio" , " pol6" , -15 , 15 );
89+ auto param = (std::vector<float >)parRatio;
90+ fRatio ->SetParameters (param[0 ], param[1 ], param[2 ], param[3 ], param[4 ], param[5 ], param[6 ]);
91+
92+ rand = new TRandom2 ();
7993 }
8094
8195 if (doprocessBinned) {
@@ -100,23 +114,54 @@ struct PseudorapidityDensity {
100114 void process (soa::Join<aod::Collisions, aod::EvSels>::iterator const & collision, Trks const & tracks)
101115 {
102116 registry.fill (HIST (" EventSelection" ), 1 .);
103- if (!useEvSel || (useEvSel && collision.sel7 ())) {
104- registry.fill (HIST (" EventSelection" ), 2 .);
105- auto z = collision.posZ ();
106- auto perCollisionSample = sample->sliceByCached (aod::track::collisionId, collision.globalIndex ());
107- if (perCollisionSample.size () > 0 ) {
108- registry.fill (HIST (" EventSelection" ), 3 .);
109- }
110- registry.fill (HIST (" EventsNtrkZvtx" ), perCollisionSample.size (), z);
111- for (auto & track : tracks) {
112- registry.fill (HIST (" TracksEtaZvtx" ), track.eta (), z);
113- registry.fill (HIST (" TracksPhiEta" ), track.phi (), track.eta ());
117+ if (doprocessGen) // we are sudying the MC, that needs correction
118+ {
119+ if (!useEvSel || (useEvSel && collision.sel7 ())) {
120+
121+ auto z = collision.posZ ();
122+ auto perCollisionSample = sample->sliceByCached (aod::track::collisionId, collision.globalIndex ());
123+ auto v = fRatio ->Eval (z);
124+
125+ double r = rand->Rndm (); // random number between 0 and 1
126+ if ((1 - v) > r) // we discard the event
127+ {
128+ registry.fill (HIST (" EventSelection" ), 4 .);
129+ return ;
130+ }
131+ registry.fill (HIST (" EventSelection" ), 2 .);
114132 if (perCollisionSample.size () > 0 ) {
115- registry.fill (HIST (" TracksEtaZvtx_gt0 " ), track. eta (), z );
133+ registry.fill (HIST (" EventSelection " ), 3 . );
116134 }
135+ registry.fill (HIST (" EventsNtrkZvtx" ), perCollisionSample.size (), z);
136+ for (auto & track : tracks) {
137+ registry.fill (HIST (" TracksEtaZvtx" ), track.eta (), z);
138+ registry.fill (HIST (" TracksPhiEta" ), track.phi (), track.eta ());
139+ if (perCollisionSample.size () > 0 ) {
140+ registry.fill (HIST (" TracksEtaZvtx_gt0" ), track.eta (), z);
141+ }
142+ }
143+ } else {
144+ registry.fill (HIST (" EventSelection" ), 4 .);
117145 }
118146 } else {
119- registry.fill (HIST (" EventSelection" ), 4 .);
147+ if (!useEvSel || (useEvSel && collision.sel7 ())) {
148+ registry.fill (HIST (" EventSelection" ), 2 .);
149+ auto z = collision.posZ ();
150+ auto perCollisionSample = sample->sliceByCached (aod::track::collisionId, collision.globalIndex ());
151+ if (perCollisionSample.size () > 0 ) {
152+ registry.fill (HIST (" EventSelection" ), 3 .);
153+ }
154+ registry.fill (HIST (" EventsNtrkZvtx" ), perCollisionSample.size (), z);
155+ for (auto & track : tracks) {
156+ registry.fill (HIST (" TracksEtaZvtx" ), track.eta (), z);
157+ registry.fill (HIST (" TracksPhiEta" ), track.phi (), track.eta ());
158+ if (perCollisionSample.size () > 0 ) {
159+ registry.fill (HIST (" TracksEtaZvtx_gt0" ), track.eta (), z);
160+ }
161+ }
162+ } else {
163+ registry.fill (HIST (" EventSelection" ), 4 .);
164+ }
120165 }
121166 }
122167
0 commit comments