CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CMS_2010_S8808686.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // Author: A. Knutsson
3 // Version: 1.0, 10/11-2010
4 //
5 //
6 // Note: You currently need to link ROOT when compiling. The histograms are
7 // standard ROOT histograms which are stored in a ROOT file. AIDA is not
8 // used.
9 //
10 #include "Rivet/Analysis.hh"
11 #include "Rivet/RivetAIDA.hh"
12 #include "Rivet/Projections/ChargedFinalState.hh"
13 #include "Rivet/Projections/Beam.hh"
14 #include "TFile.h"
15 #include "TH1F.h"
16 #include "TH2F.h"
17 
18 
19 #include <TH1D.h>
20 #include <TH2D.h>
21 #include <TNtuple.h>
22 #include <TROOT.h>
23 #include <TSystem.h>
24 #include <TString.h>
25 #include <TCanvas.h>
26 #include <TVector3.h>
27 #include <TRandom.h>
28 
29 namespace Rivet {
30 
31  class CMS_2010_S8808686 : public Analysis {
32  public:
33 
34  CMS_2010_S8808686() : Analysis("CMS_2010_S8808686") {
35  setBeams(PROTON, PROTON);
36  setNeedsCrossSection(false);
37  }
38 
39 //AK =====================================================INIT
40  void init() override {
41  ChargedFinalState cfs(-2.4, 2.4, 0.1*GeV); //Note the eta selection for charged particles is here
42  addProjection(cfs, "CFS");
43  addProjection(Beam(), "Beam");
44 
45  _N110events = 0;
46  _Nevt_after_cuts = 0;
47 
48  file = new TFile("CMS_2010_S8808686.root","recreate");
49 
50  for (int ibin = 0; ibin < 16; ibin++) {
51 
52  Char_t hname[100];
53 
54  sprintf(hname,"S_DeltaPhi_%i",ibin+1);
55  _h_S_DeltaPhi_Nptbins[ibin] = new TH1F(hname, hname,17,0.0,PI);
56 
57  sprintf(hname,"B_DeltaPhi_%i",ibin+1);
58  _h_B_DeltaPhi_Nptbins[ibin] = new TH1F(hname, hname,17,0.0,PI);
59 
60  sprintf(hname,"R_DeltaPhi_%i",ibin+1);
61  _h_R_DeltaPhi_Nptbins[ibin] = new TH1F(hname, hname,17,0.0,PI);
62 
63  sprintf(hname,"S_3D_Nptbin_%i",ibin+1);
64  _h_S_3D_Nptbins[ibin] = new TH2F(hname,hname,24,-4.8,4.8,30,-0.5*PI,3./2.*PI);
65 
66  sprintf(hname,"B_3D_Nptbin_%i",ibin+1);
67  _h_B_3D_Nptbins[ibin] = new TH2F(hname,hname,24,-4.8,4.8,30,-0.5*PI,3./2.*PI);
68 
69  sprintf(hname,"R_3D_Nptbin_%i",ibin+1);
70  _h_R_3D_Nptbins[ibin] = new TH2F(hname,hname,24,-4.8,4.8,30,-0.5*PI,3./2.*PI);
71 
72  sprintf(hname,"mult_%i",ibin+1);
73  _hMult_Nptbins[ibin] = new TH1F(hname,"N",250,0,250);
74 
75  }
76 
77  _h2_S_dphivsdeta_mb_01pt = new TH2F("S_MB_01pt","S MB 0.1<pt",33,-4.8,4.8,30,-0.5*PI,3./2.*PI);
78  _h2_B_dphivsdeta_mb_01pt = new TH2F("B_MB_01pt","B MB 0.1<pt",33,-4.8,4.8,30,-0.5*PI,3./2.*PI);
79  _h2_R_dphivsdeta_mb_01pt = new TH2F("R_MB_01pt","R MB 0.1<pt",33,-4.8,4.8,30,-0.5*PI,3./2.*PI);
80  _hMult_mb_01pt = new TH1F("mult_mb_01pt","N",250,0,250);
81 
82  _h2_S_dphivsdeta_N110_01pt = new TH2F("S_N110_01pt","S N>110 0.1<pt",33,-4.8,4.8,30,-0.5*PI,3./2.*PI);
83  _h2_B_dphivsdeta_N110_01pt = new TH2F("B_N110_01pt","B N>110 0.1<pt",33,-4.8,4.8,30,-0.5*PI,3./2.*PI);
84  _h2_R_dphivsdeta_N110_01pt = new TH2F("R_N110_01pt","R N>110 0.1<pt",33,-4.8,4.8,30,-0.5*PI,3./2.*PI);
85  _hMult_N110_01pt = new TH1F("mult_N110_01pt","N",250,0,250);
86 
87  _h2_S_dphivsdeta_mb_1pt3 = new TH2F("S_MB_1pt3","S MB 1<pt<3",33,-4.8,4.8,30,-0.5*PI,3./2.*PI);
88  _h2_B_dphivsdeta_mb_1pt3 = new TH2F("B_MB_1pt3","B MB 1<pt<3",33,-4.8,4.8,30,-0.5*PI,3./2.*PI);
89  _h2_R_dphivsdeta_mb_1pt3 = new TH2F("R_MB_1pt3","R MB 1<pt<3",33,-4.8,4.8,30,-0.5*PI,3./2.*PI);
90  _hMult_mb_1pt3 = new TH1F("mult_mb_1pt3","N",250,0,250);
91 
92  _h2_S_dphivsdeta_N110_1pt3 = new TH2F("S_N110_1pt3","S N>110 1<pt<3",33,-4.8,4.8,30,-0.5*PI,3./2.*PI);
93  _h2_B_dphivsdeta_N110_1pt3 = new TH2F("B_N110_1pt3","B N>110 1<pt<3",33,-4.8,4.8,30,-0.5*PI,3./2.*PI);
94  _h2_R_dphivsdeta_N110_1pt3 = new TH2F("R_N110_1pt3","R N>110 1<pt<3",33,-4.8,4.8,30,-0.5*PI,3./2.*PI);
95  _hMult_N110_1pt3 = new TH1F("mult_N110_1pt3","N",250,0,250);
96 
97  _hPhi_01pt = new TH1F("phi_pt01","phi 0.1<pt",100,-2*PI,2*PI);
98  }
99 
100 //AK =====================================================ANALYZE
101  void analyze(const Event& event) override {
102 
103  const double _ptbinslimits[5] = {0.1,1.0,2.0,3.0,4.0};
104  const unsigned int _Nbinslimits[5] = {1, 35, 90, 110, 9999};
105 
106  //const double weight = event.weight();
107  const ChargedFinalState& charged = applyProjection<ChargedFinalState>(event, "CFS");
108  const ParticleVector ChrgParticles = charged.particles();
109 
110  if (ChrgParticles.size() <= 1) vetoEvent;
111 
112  _Nevt_after_cuts++;
113 
114  //count particles for the signal event
115  unsigned int Nparts_01pt=0;
116  unsigned int Nparts_04pt=0;
117  unsigned int Nparts_1pt3=0;
118  unsigned int Nparts_ptbin[4]={0};
119  foreach (const Particle& p, charged.particles()) {
120  double pT = p.momentum().pT();
121  if (pT>0.1){Nparts_01pt++;}
122  if (pT>0.4){Nparts_04pt++;}
123  for (int iPtbin = 0; iPtbin < 4; iPtbin++) {
124  if (pT > _ptbinslimits[iPtbin] && pT < _ptbinslimits[iPtbin+1]){
125  Nparts_ptbin[iPtbin]++;
126  }
127  } //end - pt-bin loop
128  }
129  Nparts_1pt3 = Nparts_ptbin[1] + Nparts_ptbin[2];
130 
131  //determine the multiplcity bin and fill particle multiplcity in pt bins
132  int Nbin=-99;
133  for (int iNbin = 0; iNbin < 4; iNbin++) {
134  if(Nparts_04pt > _Nbinslimits[iNbin] && Nparts_04pt < _Nbinslimits[iNbin+1]) {
135  Nbin=iNbin; //Nbin now starts at 0
136  for (int iPtbin = 0; iPtbin < 4; iPtbin++) {
137  _hMult_Nptbins[iPtbin + Nbin*4]->Fill(Nparts_ptbin[iPtbin]);
138  }
139  }
140  }
141 
142 
143  _hMult_mb_01pt->Fill(Nparts_01pt);
144  _hMult_mb_1pt3->Fill(Nparts_1pt3);
145  if (Nbin == 3) {
146  _hMult_N110_01pt->Fill(Nparts_01pt);
147  _hMult_N110_1pt3->Fill(Nparts_1pt3);
148  }
149 
150 
151  //particle count - MB background event
152  unsigned int oldNpartsMB_01pt=0;
153  unsigned int oldNpartsMB_1pt3=0;
154  foreach (const Particle& p, _oldpartvecMB) {
155  double pT = p.momentum().pT();
156  if (pT>0.1){oldNpartsMB_01pt++;}
157  if (pT>1.0 && pT<3.0){oldNpartsMB_1pt3++;}
158  }
159 
160  //Get the background for the N classified event
161  ParticleVector oldpartvecNBin;
162  if (Nbin == 0){ oldpartvecNBin = _oldpartvec1N35;}
163  if (Nbin == 1){ oldpartvecNBin = _oldpartvec35N90;}
164  if (Nbin == 2){ oldpartvecNBin = _oldpartvec90N110;}
165  if (Nbin == 3){ oldpartvecNBin = _oldpartvec110N;}
166 
167  //particle count for the N classified background event
168  unsigned int oldNparts_01pt=0;
169  unsigned int oldNparts_1pt3=0;
170  unsigned int oldNparts_ptbin[4]={0};
171  foreach (const Particle& p, oldpartvecNBin) {
172  double pT = p.momentum().pT();
173  if(pT>0.1) oldNparts_01pt++;
174  for (int iPtbin = 0; iPtbin < 4; iPtbin++) {
175  if (pT > _ptbinslimits[iPtbin] && pT <= _ptbinslimits[iPtbin+1]){
176  oldNparts_ptbin[iPtbin]++;
177  }
178  } //end - pt-bin loop
179  }
180  oldNparts_1pt3 = oldNparts_ptbin[1] + oldNparts_ptbin[2];
181 
182 
183  if(oldpartvecNBin.size() > _Nbinslimits[Nbin] && _oldpartvecMB.size() > 1 ) { //only carry on with filling plots if we already have a one background event
184 
185 
186 
187  for (unsigned int ip1 = 0; ip1 < ChrgParticles.size(); ip1++) {
188  const Particle& p1 = ChrgParticles[ip1];
189 
190  double pT1 = p1.momentum().pT();
191  double eta1 = p1.momentum().eta();
192  double phi1 = p1.momentum().phi();
193 
194  _hPhi_01pt->Fill(phi1, 1.0); //just test histo
195 
196 
197  // Loop same event for S-distributions
198  for (unsigned int ip2 = ip1+1; ip2 < ChrgParticles.size(); ip2++) {
199  const Particle& p2 = ChrgParticles[ip2];
200 
201  double pT2 = p2.momentum().pT();
202  double eta2 = p2.momentum().eta();
203  double phi2 = p2.momentum().phi();
204 
205  double deta = fabs(eta1-eta2);
206  double dphi = phi1-phi2;
207 
208  if(dphi>PI) dphi=dphi-2*PI;
209  if(dphi<-PI) dphi=dphi+2*PI;
210 
211  //1D (R vs DeltaPhi) - Signal
212  for (int iPtbin = 0; iPtbin < 4; iPtbin++) {
213  if (pT1 > _ptbinslimits[iPtbin] && pT1 < _ptbinslimits[iPtbin+1] && pT2 > _ptbinslimits[iPtbin] && pT2 < _ptbinslimits[iPtbin+1]){
214  int ibin = iPtbin + Nbin*4; //which histo to fill, 4x4
215  _h_S_3D_Nptbins[ibin]->Fill(fabs(deta),fabs(dphi),1.0/4.0/Nparts_ptbin[iPtbin]);
216  _h_S_3D_Nptbins[ibin]->Fill(-fabs(deta),fabs(dphi),1.0/4.0/Nparts_ptbin[iPtbin]);
217  _h_S_3D_Nptbins[ibin]->Fill(fabs(deta),-fabs(dphi),1.0/4.0/Nparts_ptbin[iPtbin]);
218  _h_S_3D_Nptbins[ibin]->Fill(-fabs(deta),-fabs(dphi),1.0/4.0/Nparts_ptbin[iPtbin]);
219  _h_S_3D_Nptbins[ibin]->Fill(fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_ptbin[iPtbin]);
220  _h_S_3D_Nptbins[ibin]->Fill(-fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_ptbin[iPtbin]);
221 
222  if(deta > 2.0 && deta < 4.8){
223  _h_S_DeltaPhi_Nptbins[ibin]->Fill(fabs(dphi),1.0/Nparts_ptbin[iPtbin]);
224  }
225 
226  }
227  } //end - pt-bin loop
228 
229 
230  //3D plots (R vs DeltaPhi vs DeltaEta) - Signal
231  if (pT1 >= 0.1 && pT2 >= 0.1){
232  _h2_S_dphivsdeta_mb_01pt->Fill(fabs(deta),fabs(dphi),1.0/4.0/Nparts_01pt);
233  _h2_S_dphivsdeta_mb_01pt->Fill(-fabs(deta),fabs(dphi),1.0/4.0/Nparts_01pt);
234  _h2_S_dphivsdeta_mb_01pt->Fill(fabs(deta),-fabs(dphi),1.0/4.0/Nparts_01pt);
235  _h2_S_dphivsdeta_mb_01pt->Fill(-fabs(deta),-fabs(dphi),1.0/4.0/Nparts_01pt);
236  _h2_S_dphivsdeta_mb_01pt->Fill(fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_01pt);
237  _h2_S_dphivsdeta_mb_01pt->Fill(-fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_01pt);
238 
239  //N>110 (3D - S - low pt)
240  if(Nparts_04pt>=110){
241  _h2_S_dphivsdeta_N110_01pt->Fill(fabs(deta),fabs(dphi),1.0/4.0/Nparts_01pt);
242  _h2_S_dphivsdeta_N110_01pt->Fill(-fabs(deta),fabs(dphi),1.0/4.0/Nparts_01pt);
243  _h2_S_dphivsdeta_N110_01pt->Fill(fabs(deta),-fabs(dphi),1.0/4.0/Nparts_01pt);
244  _h2_S_dphivsdeta_N110_01pt->Fill(-fabs(deta),-fabs(dphi),1.0/4.0/Nparts_01pt);
245  _h2_S_dphivsdeta_N110_01pt->Fill(fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_01pt);
246  _h2_S_dphivsdeta_N110_01pt->Fill(-fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_01pt);
247  } //end - N>110
248  } //end - pt cuts
249 
250 
251  if (pT1 >= 1 && pT1 <= 3 && pT2 >= 1 && pT2 <= 3){
252  _h2_S_dphivsdeta_mb_1pt3->Fill(fabs(deta),fabs(dphi),1.0/4.0/Nparts_1pt3);
253  _h2_S_dphivsdeta_mb_1pt3->Fill(-fabs(deta),fabs(dphi),1.0/4.0/Nparts_1pt3);
254  _h2_S_dphivsdeta_mb_1pt3->Fill(fabs(deta),-fabs(dphi),1.0/4.0/Nparts_1pt3);
255  _h2_S_dphivsdeta_mb_1pt3->Fill(-fabs(deta),-fabs(dphi),1.0/4.0/Nparts_1pt3);
256  _h2_S_dphivsdeta_mb_1pt3->Fill(fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_1pt3);
257  _h2_S_dphivsdeta_mb_1pt3->Fill(-fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_1pt3);
258  //N>110 (3D - S 1<pt<3)
259  if(Nparts_04pt>=110){
260  _h2_S_dphivsdeta_N110_1pt3->Fill(fabs(deta),fabs(dphi),1.0/4.0/Nparts_1pt3);
261  _h2_S_dphivsdeta_N110_1pt3->Fill(-fabs(deta),fabs(dphi),1.0/4.0/Nparts_1pt3);
262  _h2_S_dphivsdeta_N110_1pt3->Fill(fabs(deta),-fabs(dphi),1.0/4.0/Nparts_1pt3);
263  _h2_S_dphivsdeta_N110_1pt3->Fill(-fabs(deta),-fabs(dphi),1.0/4.0/Nparts_1pt3);
264  _h2_S_dphivsdeta_N110_1pt3->Fill(fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_1pt3);
265  _h2_S_dphivsdeta_N110_1pt3->Fill(-fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_1pt3);
266  } //end - N>110
267  }//end - pt cuts
268 
269 
270  } //end - 2nd particle for the current event
271 
272 
276 
277  //Background bussiness MB
278 
279  if (_oldpartvecMB.size() > 1){ //only if background is already filled
280 
281  // Loop old MB event for B-distributions
282  for (unsigned int ip2 = 0; ip2 < _oldpartvecMB.size(); ip2++) {
283  const Particle& p2 = _oldpartvecMB[ip2];
284 
285  double pT2 = p2.momentum().pT();
286  double eta2 = p2.momentum().eta();
287  double phi2 = p2.momentum().phi();
288 
289  double deta = fabs(eta1-eta2);
290  double dphi = phi1-phi2;
291 
292  if(dphi>PI) dphi=dphi-2*PI;
293  if(dphi<-PI) dphi=dphi+2*PI;
294 
295 
296  //MB (3D - B - low pt)
297  if (pT1 >= 0.1 && pT2 >= 0.1){
298  double pweight=1.0/(Nparts_01pt*oldNpartsMB_01pt);
299  _h2_B_dphivsdeta_mb_01pt->Fill(deta,fabs(dphi),pweight);
300  _h2_B_dphivsdeta_mb_01pt->Fill(-deta,fabs(dphi),pweight);
301  _h2_B_dphivsdeta_mb_01pt->Fill(deta,-fabs(dphi),pweight);
302  _h2_B_dphivsdeta_mb_01pt->Fill(-deta,-fabs(dphi),pweight);
303  _h2_B_dphivsdeta_mb_01pt->Fill(deta,2*PI-fabs(dphi),pweight);
304  _h2_B_dphivsdeta_mb_01pt->Fill(-deta,2*PI-fabs(dphi),pweight);
305  } //end - pt cuts
306 
307  //MB (3D - B - 1<pt<3)
308  if (pT1 >= 1 && pT1 <= 3 && pT2 >= 1 && pT2 <= 3){
309  double pweight=1.0/(Nparts_1pt3*oldNpartsMB_1pt3);
310  _h2_B_dphivsdeta_mb_1pt3->Fill(deta,fabs(dphi),pweight);
311  _h2_B_dphivsdeta_mb_1pt3->Fill(-deta,fabs(dphi),pweight);
312  _h2_B_dphivsdeta_mb_1pt3->Fill(deta,-fabs(dphi),pweight);
313  _h2_B_dphivsdeta_mb_1pt3->Fill(-deta,-fabs(dphi),pweight);
314  _h2_B_dphivsdeta_mb_1pt3->Fill(deta,2*PI-fabs(dphi),pweight);
315  _h2_B_dphivsdeta_mb_1pt3->Fill(-deta,2*PI-fabs(dphi),pweight);
316  }//end - pt cuts
317 
318 
319  } //end - particle loop over saved MB particle vector
320 
321  } //end - need atleast 1 background event
322 
323 
324  //Background bussiness N>110
325 
326  if (oldpartvecNBin.size() > 100 && Nbin==3){ //only if it is already filled
327 
328  // Loop old N110 event for B-distributions
329  for (unsigned int ip2 = 0; ip2 < oldpartvecNBin.size(); ip2++) {
330  const Particle& p2 = oldpartvecNBin[ip2];
331 
332  double pT2 = p2.momentum().pT();
333  double eta2 = p2.momentum().eta();
334  double phi2 = p2.momentum().phi();
335 
336  double deta = fabs(eta1-eta2);
337  double dphi = phi1-phi2;
338 
339  if(dphi>PI) dphi=dphi-2*PI;
340  if(dphi<-PI) dphi=dphi+2*PI;
341 // if(dphi>-PI && dphi<-PI/2.) dphi=dphi+2*PI;
342 
343 
344  if (pT1 >= 0.1 && pT2 >= 0.1){
345  //Fill
346  double pweight=1.0/(Nparts_01pt*oldNparts_01pt);
347  _h2_B_dphivsdeta_N110_01pt->Fill(deta,fabs(dphi),pweight);
348  _h2_B_dphivsdeta_N110_01pt->Fill(-deta,fabs(dphi),pweight);
349  _h2_B_dphivsdeta_N110_01pt->Fill(deta,-fabs(dphi),pweight);
350  _h2_B_dphivsdeta_N110_01pt->Fill(-deta,-fabs(dphi),pweight);
351  _h2_B_dphivsdeta_N110_01pt->Fill(deta,2*PI-fabs(dphi),pweight);
352  _h2_B_dphivsdeta_N110_01pt->Fill(-deta,2*PI-fabs(dphi),pweight);
353  } //end - pt cuts
354 
355  if (pT1 >= 1 && pT1 <= 3 && pT2 >= 1 && pT2 <= 3){
356  //Fill
357  double pweight=1.0/(Nparts_1pt3*oldNparts_1pt3);
358  _h2_B_dphivsdeta_N110_1pt3->Fill(deta,fabs(dphi),pweight);
359  _h2_B_dphivsdeta_N110_1pt3->Fill(-deta,fabs(dphi),pweight);
360  _h2_B_dphivsdeta_N110_1pt3->Fill(deta,-fabs(dphi),pweight);
361  _h2_B_dphivsdeta_N110_1pt3->Fill(-deta,-fabs(dphi),pweight);
362  _h2_B_dphivsdeta_N110_1pt3->Fill(deta,2*PI-fabs(dphi),pweight);
363  _h2_B_dphivsdeta_N110_1pt3->Fill(-deta,2*PI-fabs(dphi),pweight);
364  }//end - pt cuts
365 
366 
367  } //end - particle loop old N110
368 
369  } //end - need atleast 1 event already
370 
371 
372 
376  if (oldpartvecNBin.size() > _Nbinslimits[Nbin]){ //only if we already have BG particles
377 
378  // Loop old event for B-distributions
379  for (unsigned int ip2 = 0; ip2 < oldpartvecNBin.size(); ip2++) {
380  const Particle& p2 = oldpartvecNBin[ip2];
381 
382  double pT2 = p2.momentum().pT();
383  double eta2 = p2.momentum().eta();
384  double phi2 = p2.momentum().phi();
385 
386  double deta = fabs(eta1-eta2);
387  double dphi = phi1-phi2;
388 
389  if(dphi>PI) dphi=dphi-2*PI;
390  if(dphi<-PI) dphi=dphi+2*PI;
391 
392  //loop the pt bins for the DeltaPhi - 1D - Background
393  for (int iPtbin = 0; iPtbin < 4; iPtbin++) {
394  if (pT1 > _ptbinslimits[iPtbin] && pT1 < _ptbinslimits[iPtbin+1] && pT2 > _ptbinslimits[iPtbin] && pT2 < _ptbinslimits[iPtbin+1]){
395  int ibin = iPtbin + Nbin*4; //which histo to fill, 4x4 matix
396  double pweight=1.0/(Nparts_ptbin[iPtbin]*oldNparts_ptbin[iPtbin]);
397  _h_B_3D_Nptbins[ibin]->Fill(deta,fabs(dphi),pweight);
398  _h_B_3D_Nptbins[ibin]->Fill(-deta,fabs(dphi),pweight);
399  _h_B_3D_Nptbins[ibin]->Fill(deta,-fabs(dphi),pweight);
400  _h_B_3D_Nptbins[ibin]->Fill(-deta,-fabs(dphi),pweight);
401  _h_B_3D_Nptbins[ibin]->Fill(deta,2*PI-fabs(dphi),pweight);
402  _h_B_3D_Nptbins[ibin]->Fill(-deta,2*PI-fabs(dphi),pweight);
403  if(deta > 2.0 && deta < 4.8){
404  _h_B_DeltaPhi_Nptbins[ibin]->Fill(fabs(dphi),pweight);
405  }
406  } //end check pt-bin
407  } //end - pt-bin loop
408 
409 
410  } //end - particle loop old N110
411  } //end - need atleast 1 event already
412 
413 
414  } //end - main particle loop for current event
415  } //end - if background events found
416 
417 
418 //save the old particle vector
419  if (Nbin == 0){ _oldpartvec1N35 = ChrgParticles;}
420  if (Nbin == 1){ _oldpartvec35N90 = ChrgParticles;}
421  if (Nbin == 2){ _oldpartvec90N110 = ChrgParticles;}
422  if (Nbin == 3){ _oldpartvec110N = ChrgParticles; _N110events++;}
423  if (Nparts_01pt >= 1){ _oldpartvecMB = ChrgParticles;}
424 
425  } //end - analyze()
426 
427 //AK =====================================================FINALIZE
428  void finalize() override {
429 
430  getLog() << Log::INFO << "Number of events after event selection: " << _Nevt_after_cuts << endl;
431  getLog() << Log::INFO << "Number of events with N>110:" << _N110events << endl;
432 
433  int nEvent = -99.0;
434 
435  nEvent = _hMult_mb_01pt->Integral(3,10000);
436  _h2_S_dphivsdeta_mb_01pt->Scale(1.0/nEvent);
437  _h2_B_dphivsdeta_mb_01pt->Scale(_h2_S_dphivsdeta_mb_01pt->Integral()/_h2_B_dphivsdeta_mb_01pt->Integral());
438  _h2_R_dphivsdeta_mb_01pt->Add((TH2F*)_h2_S_dphivsdeta_mb_01pt);
439  _h2_R_dphivsdeta_mb_01pt->Add(_h2_B_dphivsdeta_mb_01pt,-1);
440  _h2_R_dphivsdeta_mb_01pt->Divide(_h2_B_dphivsdeta_mb_01pt);
441  _h2_R_dphivsdeta_mb_01pt->Scale(_h2_S_dphivsdeta_mb_01pt->Integral());
442 
443  nEvent = _hMult_N110_01pt->Integral(3,10000);
444  _h2_S_dphivsdeta_N110_01pt->Scale(1.0/nEvent);
445  _h2_B_dphivsdeta_N110_01pt->Scale(_h2_S_dphivsdeta_N110_01pt->Integral()/_h2_B_dphivsdeta_N110_01pt->Integral());
446  _h2_R_dphivsdeta_N110_01pt->Add(_h2_S_dphivsdeta_N110_01pt);
447  _h2_R_dphivsdeta_N110_01pt->Add(_h2_B_dphivsdeta_N110_01pt,-1);
448  _h2_R_dphivsdeta_N110_01pt->Divide(_h2_B_dphivsdeta_N110_01pt);
449  _h2_R_dphivsdeta_N110_01pt->Scale(_h2_S_dphivsdeta_N110_01pt->Integral());
450 
451  nEvent = _hMult_mb_1pt3->Integral(3,10000);
452  _h2_S_dphivsdeta_mb_1pt3->Scale(1.0/nEvent);
453  _h2_B_dphivsdeta_mb_1pt3->Scale(_h2_S_dphivsdeta_mb_1pt3->Integral()/_h2_B_dphivsdeta_mb_1pt3->Integral());
454  _h2_R_dphivsdeta_mb_1pt3->Add(_h2_S_dphivsdeta_mb_1pt3);
455  _h2_R_dphivsdeta_mb_1pt3->Add(_h2_B_dphivsdeta_mb_1pt3,-1);
456  _h2_R_dphivsdeta_mb_1pt3->Divide(_h2_B_dphivsdeta_mb_1pt3);
457  _h2_R_dphivsdeta_mb_1pt3->Scale(_h2_S_dphivsdeta_mb_1pt3->Integral());
458 
459  nEvent = _hMult_N110_1pt3->Integral(3,10000);
460  _h2_S_dphivsdeta_N110_1pt3->Scale(1.0/nEvent);
461  _h2_B_dphivsdeta_N110_1pt3->Scale(_h2_S_dphivsdeta_N110_1pt3->Integral()/_h2_B_dphivsdeta_N110_1pt3->Integral());
462  _h2_R_dphivsdeta_N110_1pt3->Add(_h2_S_dphivsdeta_N110_1pt3);
463  _h2_R_dphivsdeta_N110_1pt3->Add(_h2_B_dphivsdeta_N110_1pt3,-1);
464  _h2_R_dphivsdeta_N110_1pt3->Divide(_h2_B_dphivsdeta_N110_1pt3);
465  _h2_R_dphivsdeta_N110_1pt3->Scale(_h2_S_dphivsdeta_N110_1pt3->Integral());
466 
467  for (int ibin = 0; ibin < 16; ibin++) {
468 
469  int nEvent = _hMult_Nptbins[ibin]->Integral(3,10000);
470 
471  _h_S_3D_Nptbins[ibin]->Scale(1.0/nEvent);
472  _h_B_3D_Nptbins[ibin]->Scale(_h_S_3D_Nptbins[ibin]->Integral()/_h_B_3D_Nptbins[ibin]->Integral());
473  _h_R_3D_Nptbins[ibin]->Add(_h_S_3D_Nptbins[ibin]);
474  _h_R_3D_Nptbins[ibin]->Add(_h_B_3D_Nptbins[ibin],-1);
475  _h_R_3D_Nptbins[ibin]->Divide(_h_B_3D_Nptbins[ibin]);
476  _h_R_3D_Nptbins[ibin]->Scale(_h_S_3D_Nptbins[ibin]->Integral());
477 
478 
479  _h_S_DeltaPhi_Nptbins[ibin]->Scale(1.0/nEvent);
480  // Wei's correction - to get the average multiplicity correct --
481  _hMult_Nptbins[ibin]->SetAxisRange(2,10000,"X");
482  double rescale = (_hMult_Nptbins[ibin]->GetMean()-1)/_h_S_DeltaPhi_Nptbins[ibin]->Integral();
483  _h_S_DeltaPhi_Nptbins[ibin]->Scale(rescale);
484  //---------------------------------------------------------------
485  _h_B_DeltaPhi_Nptbins[ibin]->Scale(_h_S_DeltaPhi_Nptbins[ibin]->Integral()/_h_B_DeltaPhi_Nptbins[ibin]->Integral());
486  _h_R_DeltaPhi_Nptbins[ibin]->Add(_h_S_DeltaPhi_Nptbins[ibin]);
487  _h_R_DeltaPhi_Nptbins[ibin]->Add(_h_B_DeltaPhi_Nptbins[ibin],-1);
488  _h_R_DeltaPhi_Nptbins[ibin]->Divide(_h_B_DeltaPhi_Nptbins[ibin]);
489  _h_R_DeltaPhi_Nptbins[ibin]->Scale(_h_S_DeltaPhi_Nptbins[ibin]->Integral());
490 
491 
492  }
493 
494  file -> Write();
495 
496  }
497 
498 
499 //AK =====================================================DECLARATIONS
500  private:
501 
502  double detamin;
503  double detamax;
504  int Ndetabins;
505 
506  double _N110events;
507  double _Nevt_after_cuts;
508 
509 
510  TH1F *_hMult;
511  TH1F *_hPhi_01pt;
512  TH1F *_h_S_DeltaPhi_Nptbins[16];
513  TH1F *_h_B_DeltaPhi_Nptbins[16];
514  TH1F *_h_R_DeltaPhi_Nptbins[16];
515 
516  TH2F *_h_S_3D_Nptbins[16];
517  TH2F *_h_B_3D_Nptbins[16];
518  TH2F *_h_R_3D_Nptbins[16];
519  TH1F *_hMult_Nptbins[16];
520 
521  TH2F *_h2_S_dphivsdeta_mb_01pt;
522  TH2F *_h2_B_dphivsdeta_mb_01pt;
523  TH2F *_h2_R_dphivsdeta_mb_01pt;
524  TH1F *_hMult_mb_01pt;
525 
526  TH2F *_h2_S_dphivsdeta_mb_1pt3;
527  TH2F *_h2_B_dphivsdeta_mb_1pt3;
528  TH2F *_h2_R_dphivsdeta_mb_1pt3;
529  TH1F *_hMult_mb_1pt3;
530 
531  TH2F *_h2_S_dphivsdeta_N110_01pt;
532  TH2F *_h2_B_dphivsdeta_N110_01pt;
533  TH2F *_h2_R_dphivsdeta_N110_01pt;
534  TH1F *_hMult_N110_01pt;
535 
536  TH2F *_h2_S_dphivsdeta_N110_1pt3;
537  TH2F *_h2_B_dphivsdeta_N110_1pt3;
538  TH2F *_h2_R_dphivsdeta_N110_1pt3;
539  TH1F *_hMult_N110_1pt3;
540 
541  ParticleVector _oldpartvec1N35;
542  ParticleVector _oldpartvec35N90;
543  ParticleVector _oldpartvec90N110; //TODO: this is a bit ugly -> arrays or vec<vec<>>
544  ParticleVector _oldpartvec110N;
545  ParticleVector _oldpartvecMB;
546 
547  TFile *file;
548  };
549 
550 
551 
552  // This global object acts as a hook for the plugin system
553  AnalysisBuilder<CMS_2010_S8808686> plugin_CMS_2010_S8808686;
554 
555 }
556 
Rivet::CMS_2010_S8808686::analyze
void analyze(const Event &event) override
Definition: CMS_2010_S8808686.cc:101
Rivet::CMS_2010_S8808686::_hMult_mb_01pt
TH1F * _hMult_mb_01pt
Definition: CMS_2010_S8808686.cc:524
Rivet::CMS_2010_S8808686::_oldpartvec110N
ParticleVector _oldpartvec110N
Definition: CMS_2010_S8808686.cc:544
Rivet::CMS_2010_S8808686::_h2_S_dphivsdeta_N110_1pt3
TH2F * _h2_S_dphivsdeta_N110_1pt3
Definition: CMS_2010_S8808686.cc:536
GeV
const double GeV
Definition: MathUtil.h:16
Rivet::CMS_2010_S8808686::_h2_B_dphivsdeta_mb_01pt
TH2F * _h2_B_dphivsdeta_mb_01pt
Definition: CMS_2010_S8808686.cc:522
Rivet::CMS_2010_S8808686::_h_R_DeltaPhi_Nptbins
TH1F * _h_R_DeltaPhi_Nptbins[16]
Definition: CMS_2010_S8808686.cc:514
dtDQMClient_cfg.INFO
tuple INFO
Definition: dtDQMClient_cfg.py:19
PI
#define PI
Definition: HcalBeamMonitor.cc:14
Rivet::CMS_2010_S8808686::_hMult_N110_01pt
TH1F * _hMult_N110_01pt
Definition: CMS_2010_S8808686.cc:534
Rivet::CMS_2010_S8808686::_hPhi_01pt
TH1F * _hPhi_01pt
Definition: CMS_2010_S8808686.cc:511
Rivet::CMS_2010_S8808686::_h_S_DeltaPhi_Nptbins
TH1F * _h_S_DeltaPhi_Nptbins[16]
Definition: CMS_2010_S8808686.cc:512
Rivet::CMS_2010_S8808686::_oldpartvec35N90
ParticleVector _oldpartvec35N90
Definition: CMS_2010_S8808686.cc:542
nEvent
int nEvent
Definition: myFastSimVal.cc:49
Rivet::plugin_CMS_2010_S8808686
AnalysisBuilder< CMS_2010_S8808686 > plugin_CMS_2010_S8808686
Definition: CMS_2010_S8808686.cc:553
Rivet::CMS_2010_S8808686::_hMult_N110_1pt3
TH1F * _hMult_N110_1pt3
Definition: CMS_2010_S8808686.cc:539
Rivet::CMS_2010_S8808686
Definition: CMS_2010_S8808686.cc:31
Rivet::CMS_2010_S8808686::_hMult_Nptbins
TH1F * _hMult_Nptbins[16]
Definition: CMS_2010_S8808686.cc:519
Rivet::CMS_2010_S8808686::file
TFile * file
Definition: CMS_2010_S8808686.cc:547
Rivet::CMS_2010_S8808686::_h_S_3D_Nptbins
TH2F * _h_S_3D_Nptbins[16]
Definition: CMS_2010_S8808686.cc:516
Rivet::CMS_2010_S8808686::_h_B_3D_Nptbins
TH2F * _h_B_3D_Nptbins[16]
Definition: CMS_2010_S8808686.cc:517
Rivet::CMS_2010_S8808686::Ndetabins
int Ndetabins
Definition: CMS_2010_S8808686.cc:504
Rivet::CMS_2010_S8808686::_h2_R_dphivsdeta_mb_01pt
TH2F * _h2_R_dphivsdeta_mb_01pt
Definition: CMS_2010_S8808686.cc:523
event
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
Rivet::CMS_2010_S8808686::_h2_R_dphivsdeta_N110_1pt3
TH2F * _h2_R_dphivsdeta_N110_1pt3
Definition: CMS_2010_S8808686.cc:538
Rivet::CMS_2010_S8808686::init
void init() override
Definition: CMS_2010_S8808686.cc:40
p2
double p2[4]
Definition: TauolaWrapper.h:90
Rivet::CMS_2010_S8808686::_h_R_3D_Nptbins
TH2F * _h_R_3D_Nptbins[16]
Definition: CMS_2010_S8808686.cc:518
Rivet::CMS_2010_S8808686::_Nevt_after_cuts
double _Nevt_after_cuts
Definition: CMS_2010_S8808686.cc:507
Rivet::CMS_2010_S8808686::_hMult_mb_1pt3
TH1F * _hMult_mb_1pt3
Definition: CMS_2010_S8808686.cc:529
Rivet::CMS_2010_S8808686::_h2_B_dphivsdeta_N110_1pt3
TH2F * _h2_B_dphivsdeta_N110_1pt3
Definition: CMS_2010_S8808686.cc:537
Rivet::CMS_2010_S8808686::_oldpartvecMB
ParticleVector _oldpartvecMB
Definition: CMS_2010_S8808686.cc:545
Rivet::CMS_2010_S8808686::_h2_R_dphivsdeta_mb_1pt3
TH2F * _h2_R_dphivsdeta_mb_1pt3
Definition: CMS_2010_S8808686.cc:528
Rivet::CMS_2010_S8808686::_hMult
TH1F * _hMult
Definition: CMS_2010_S8808686.cc:510
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
p1
double p1[4]
Definition: TauolaWrapper.h:89
Rivet::CMS_2010_S8808686::_h2_S_dphivsdeta_mb_1pt3
TH2F * _h2_S_dphivsdeta_mb_1pt3
Definition: CMS_2010_S8808686.cc:526
Rivet::CMS_2010_S8808686::CMS_2010_S8808686
CMS_2010_S8808686()
Definition: CMS_2010_S8808686.cc:34
Rivet::CMS_2010_S8808686::detamax
double detamax
Definition: CMS_2010_S8808686.cc:503
Rivet::CMS_2010_S8808686::_oldpartvec90N110
ParticleVector _oldpartvec90N110
Definition: CMS_2010_S8808686.cc:543
Rivet::CMS_2010_S8808686::detamin
double detamin
Definition: CMS_2010_S8808686.cc:502
Rivet::CMS_2010_S8808686::_h2_S_dphivsdeta_mb_01pt
TH2F * _h2_S_dphivsdeta_mb_01pt
Definition: CMS_2010_S8808686.cc:521
Rivet::CMS_2010_S8808686::_oldpartvec1N35
ParticleVector _oldpartvec1N35
Definition: CMS_2010_S8808686.cc:541
Rivet::CMS_2010_S8808686::_h_B_DeltaPhi_Nptbins
TH1F * _h_B_DeltaPhi_Nptbins[16]
Definition: CMS_2010_S8808686.cc:513
Rivet::CMS_2010_S8808686::_h2_B_dphivsdeta_mb_1pt3
TH2F * _h2_B_dphivsdeta_mb_1pt3
Definition: CMS_2010_S8808686.cc:527
Rivet::CMS_2010_S8808686::_h2_R_dphivsdeta_N110_01pt
TH2F * _h2_R_dphivsdeta_N110_01pt
Definition: CMS_2010_S8808686.cc:533
Rivet::CMS_2010_S8808686::_N110events
double _N110events
Definition: CMS_2010_S8808686.cc:506
Rivet::CMS_2010_S8808686::_h2_B_dphivsdeta_N110_01pt
TH2F * _h2_B_dphivsdeta_N110_01pt
Definition: CMS_2010_S8808686.cc:532
Rivet::CMS_2010_S8808686::finalize
void finalize() override
Definition: CMS_2010_S8808686.cc:428
Rivet::CMS_2010_S8808686::_h2_S_dphivsdeta_N110_01pt
TH2F * _h2_S_dphivsdeta_N110_01pt
Definition: CMS_2010_S8808686.cc:531
Generated for CMSSW Reference Manual by doxygen 1.8.5