Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*- 
00002 // Author: A. Knutsson
00003 // Version: 1.0, 10/11-2010
00004 //              
00005 //
00006 // Note: You currently need to link ROOT when compiling. The histograms are
00007 //       standard ROOT histograms which are stored in a ROOT file. AIDA is not
00008 //       used.
00009 //
00010 #include "Rivet/Analysis.hh"
00011 #include "Rivet/RivetAIDA.hh"
00012 #include "Rivet/Projections/ChargedFinalState.hh"
00013 #include "Rivet/Projections/Beam.hh"
00014 #include "TFile.h"
00015 #include "TH1F.h"
00016 #include "TH2F.h"
00017 
00018 
00019 #include <TH1D.h>
00020 #include <TH2D.h>
00021 #include <TNtuple.h>
00022 #include <TROOT.h>
00023 #include <TSystem.h>
00024 #include <TString.h>
00025 #include <TCanvas.h>
00026 #include <TVector3.h>
00027 #include <TRandom.h>
00028 
00029 namespace Rivet {
00030 
00031   class CMS_2010_S8808686 : public Analysis {
00032   public:
00033 
00034     CMS_2010_S8808686() : Analysis("CMS_2010_S8808686") {
00035        setBeams(PROTON, PROTON);
00036        setNeedsCrossSection(false);
00037     }
00038 
00039 //AK =====================================================INIT
00040     void init() {
00041       ChargedFinalState cfs(-2.4, 2.4, 0.1*GeV); //Note the eta selection for charged particles is here
00042       addProjection(cfs, "CFS");
00043       addProjection(Beam(), "Beam");
00044         
00045         _N110events = 0;
00046         _Nevt_after_cuts = 0;
00047         
00048          file = new TFile("CMS_2010_S8808686.root","recreate");
00049         
00050         for (int ibin = 0; ibin < 16; ibin++) {
00051            
00052           Char_t hname[100];
00053         
00054           sprintf(hname,"S_DeltaPhi_%i",ibin+1);
00055           _h_S_DeltaPhi_Nptbins[ibin] = new TH1F(hname, hname,17,0.0,PI);
00056         
00057           sprintf(hname,"B_DeltaPhi_%i",ibin+1);
00058           _h_B_DeltaPhi_Nptbins[ibin] = new TH1F(hname, hname,17,0.0,PI);
00059         
00060           sprintf(hname,"R_DeltaPhi_%i",ibin+1);
00061           _h_R_DeltaPhi_Nptbins[ibin] = new TH1F(hname, hname,17,0.0,PI);
00062         
00063           sprintf(hname,"S_3D_Nptbin_%i",ibin+1);
00064           _h_S_3D_Nptbins[ibin] = new TH2F(hname,hname,24,-4.8,4.8,30,-0.5*PI,3./2.*PI);
00065         
00066           sprintf(hname,"B_3D_Nptbin_%i",ibin+1);
00067           _h_B_3D_Nptbins[ibin] = new TH2F(hname,hname,24,-4.8,4.8,30,-0.5*PI,3./2.*PI);
00068                 
00069           sprintf(hname,"R_3D_Nptbin_%i",ibin+1);
00070           _h_R_3D_Nptbins[ibin] = new TH2F(hname,hname,24,-4.8,4.8,30,-0.5*PI,3./2.*PI);
00071 
00072           sprintf(hname,"mult_%i",ibin+1);
00073           _hMult_Nptbins[ibin] = new  TH1F(hname,"N",250,0,250); 
00074                 
00075         }
00076 
00077         _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);
00078         _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);
00079         _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);
00080         _hMult_mb_01pt = new  TH1F("mult_mb_01pt","N",250,0,250); 
00081         
00082         _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);
00083         _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);
00084         _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);
00085         _hMult_N110_01pt = new  TH1F("mult_N110_01pt","N",250,0,250); 
00086        
00087         _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);
00088         _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);
00089         _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);
00090         _hMult_mb_1pt3 = new  TH1F("mult_mb_1pt3","N",250,0,250); 
00091 
00092         _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);
00093         _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);
00094         _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);
00095         _hMult_N110_1pt3 = new  TH1F("mult_N110_1pt3","N",250,0,250); 
00096         
00097         _hPhi_01pt = new TH1F("phi_pt01","phi 0.1<pt",100,-2*PI,2*PI);
00098     }
00099 
00100 //AK =====================================================ANALYZE
00101     void analyze(const Event& event) {
00102         
00103       const double _ptbinslimits[5] = {0.1,1.0,2.0,3.0,4.0};
00104       const unsigned int _Nbinslimits[5] = {1, 35, 90, 110, 9999};
00105 
00106       //const double weight = event.weight();
00107       const ChargedFinalState& charged = applyProjection<ChargedFinalState>(event, "CFS");
00108       const ParticleVector ChrgParticles = charged.particles();
00109         
00110       if (ChrgParticles.size() <= 1) vetoEvent;
00111 
00112       _Nevt_after_cuts++;
00113  
00114       //count particles for the signal event
00115       unsigned int Nparts_01pt=0;
00116       unsigned int Nparts_04pt=0;
00117       unsigned int Nparts_1pt3=0;
00118       unsigned int Nparts_ptbin[4]={0};
00119       foreach (const Particle& p, charged.particles()) {
00120         double pT = p.momentum().pT();
00121         if (pT>0.1){Nparts_01pt++;}
00122         if (pT>0.4){Nparts_04pt++;}
00123         for (int iPtbin = 0; iPtbin < 4; iPtbin++) {
00124           if (pT > _ptbinslimits[iPtbin] && pT < _ptbinslimits[iPtbin+1]){                               
00125              Nparts_ptbin[iPtbin]++;
00126            }
00127          } //end - pt-bin loop   
00128       }
00129       Nparts_1pt3 =  Nparts_ptbin[1] + Nparts_ptbin[2];
00130             
00131       //determine the multiplcity bin and fill particle multiplcity in pt bins
00132       int Nbin=-99;
00133       for (int iNbin = 0; iNbin < 4; iNbin++) {
00134         if(Nparts_04pt > _Nbinslimits[iNbin] && Nparts_04pt < _Nbinslimits[iNbin+1]) {
00135             Nbin=iNbin; //Nbin now starts at 0
00136             for (int iPtbin = 0; iPtbin < 4; iPtbin++) {
00137               _hMult_Nptbins[iPtbin + Nbin*4]->Fill(Nparts_ptbin[iPtbin]);
00138             }       
00139          }    
00140        }
00141        
00142        
00143        _hMult_mb_01pt->Fill(Nparts_01pt);
00144        _hMult_mb_1pt3->Fill(Nparts_1pt3);
00145        if (Nbin == 3) {
00146           _hMult_N110_01pt->Fill(Nparts_01pt);
00147           _hMult_N110_1pt3->Fill(Nparts_1pt3);
00148        } 
00149        
00150        
00151        //particle count - MB background event
00152       unsigned int oldNpartsMB_01pt=0;
00153       unsigned int oldNpartsMB_1pt3=0;
00154       foreach (const Particle& p, _oldpartvecMB) {
00155         double pT = p.momentum().pT();
00156         if (pT>0.1){oldNpartsMB_01pt++;}
00157         if (pT>1.0 && pT<3.0){oldNpartsMB_1pt3++;}
00158       }         
00159  
00160       //Get the background for the N classified event   
00161       ParticleVector oldpartvecNBin;
00162       if (Nbin == 0){ oldpartvecNBin = _oldpartvec1N35;}
00163       if (Nbin == 1){ oldpartvecNBin = _oldpartvec35N90;}
00164       if (Nbin == 2){ oldpartvecNBin = _oldpartvec90N110;}
00165       if (Nbin == 3){ oldpartvecNBin = _oldpartvec110N;} 
00166       
00167       //particle count for the N classified background event    
00168       unsigned int oldNparts_01pt=0;
00169       unsigned int oldNparts_1pt3=0;
00170       unsigned int oldNparts_ptbin[4]={0};
00171       foreach (const Particle& p, oldpartvecNBin)  {
00172          double pT = p.momentum().pT();
00173          if(pT>0.1) oldNparts_01pt++;
00174          for (int iPtbin = 0; iPtbin < 4; iPtbin++) {
00175            if (pT > _ptbinslimits[iPtbin] && pT <= _ptbinslimits[iPtbin+1]){                             
00176             oldNparts_ptbin[iPtbin]++;
00177            }
00178          } //end - pt-bin loop   
00179        }
00180        oldNparts_1pt3 =  oldNparts_ptbin[1] + oldNparts_ptbin[2];
00181         
00182 
00183       if(oldpartvecNBin.size() > _Nbinslimits[Nbin] && _oldpartvecMB.size() > 1 ) {  //only carry on with filling plots if we already have a one background event
00184       
00185       
00186       
00187       for (unsigned int ip1 = 0; ip1 < ChrgParticles.size(); ip1++) {
00188         const Particle& p1 = ChrgParticles[ip1];
00189 
00190         double pT1 = p1.momentum().pT();        
00191         double eta1 = p1.momentum().eta();
00192         double phi1 = p1.momentum().phi();
00193         
00194         _hPhi_01pt->Fill(phi1, 1.0); //just test histo    
00195 
00196         
00197         // Loop same event for S-distributions
00198        for (unsigned int ip2 = ip1+1; ip2 < ChrgParticles.size(); ip2++) {
00199         const Particle& p2 = ChrgParticles[ip2];
00200           
00201           double pT2 = p2.momentum().pT();              
00202           double eta2 = p2.momentum().eta();
00203           double phi2 = p2.momentum().phi();
00204         
00205           double deta = fabs(eta1-eta2);
00206           double dphi = phi1-phi2;
00207                           
00208         if(dphi>PI) dphi=dphi-2*PI;
00209         if(dphi<-PI) dphi=dphi+2*PI;
00210           
00211           //1D (R vs DeltaPhi) - Signal
00212             for (int iPtbin = 0; iPtbin < 4; iPtbin++) {
00213               if (pT1 > _ptbinslimits[iPtbin] && pT1 < _ptbinslimits[iPtbin+1] &&  pT2 >  _ptbinslimits[iPtbin] && pT2 <  _ptbinslimits[iPtbin+1]){                              
00214                 int ibin = iPtbin + Nbin*4;  //which histo to fill, 4x4
00215                 _h_S_3D_Nptbins[ibin]->Fill(fabs(deta),fabs(dphi),1.0/4.0/Nparts_ptbin[iPtbin]);
00216                 _h_S_3D_Nptbins[ibin]->Fill(-fabs(deta),fabs(dphi),1.0/4.0/Nparts_ptbin[iPtbin]);
00217                 _h_S_3D_Nptbins[ibin]->Fill(fabs(deta),-fabs(dphi),1.0/4.0/Nparts_ptbin[iPtbin]);
00218                 _h_S_3D_Nptbins[ibin]->Fill(-fabs(deta),-fabs(dphi),1.0/4.0/Nparts_ptbin[iPtbin]);
00219                 _h_S_3D_Nptbins[ibin]->Fill(fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_ptbin[iPtbin]);
00220                 _h_S_3D_Nptbins[ibin]->Fill(-fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_ptbin[iPtbin]);
00221                 
00222                 if(deta > 2.0 && deta < 4.8){
00223                   _h_S_DeltaPhi_Nptbins[ibin]->Fill(fabs(dphi),1.0/Nparts_ptbin[iPtbin]);
00224                 }
00225                 
00226               }
00227             } //end - pt-bin loop 
00228           
00229           
00230           //3D plots (R vs DeltaPhi vs DeltaEta) - Signal
00231           if (pT1 >= 0.1 && pT2 >= 0.1){                  
00232                 _h2_S_dphivsdeta_mb_01pt->Fill(fabs(deta),fabs(dphi),1.0/4.0/Nparts_01pt);
00233                 _h2_S_dphivsdeta_mb_01pt->Fill(-fabs(deta),fabs(dphi),1.0/4.0/Nparts_01pt);
00234                 _h2_S_dphivsdeta_mb_01pt->Fill(fabs(deta),-fabs(dphi),1.0/4.0/Nparts_01pt);
00235                 _h2_S_dphivsdeta_mb_01pt->Fill(-fabs(deta),-fabs(dphi),1.0/4.0/Nparts_01pt);
00236                 _h2_S_dphivsdeta_mb_01pt->Fill(fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_01pt);
00237                 _h2_S_dphivsdeta_mb_01pt->Fill(-fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_01pt);
00238 
00239               //N>110 (3D - S - low pt)
00240              if(Nparts_04pt>=110){
00241                 _h2_S_dphivsdeta_N110_01pt->Fill(fabs(deta),fabs(dphi),1.0/4.0/Nparts_01pt);
00242                 _h2_S_dphivsdeta_N110_01pt->Fill(-fabs(deta),fabs(dphi),1.0/4.0/Nparts_01pt);
00243                 _h2_S_dphivsdeta_N110_01pt->Fill(fabs(deta),-fabs(dphi),1.0/4.0/Nparts_01pt);
00244                 _h2_S_dphivsdeta_N110_01pt->Fill(-fabs(deta),-fabs(dphi),1.0/4.0/Nparts_01pt);
00245                 _h2_S_dphivsdeta_N110_01pt->Fill(fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_01pt);
00246                 _h2_S_dphivsdeta_N110_01pt->Fill(-fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_01pt);
00247              }  //end - N>110        
00248           } //end - pt cuts
00249         
00250                   
00251           if (pT1 >= 1 && pT1 <= 3 && pT2 >= 1 && pT2 <= 3){ 
00252                 _h2_S_dphivsdeta_mb_1pt3->Fill(fabs(deta),fabs(dphi),1.0/4.0/Nparts_1pt3);
00253                 _h2_S_dphivsdeta_mb_1pt3->Fill(-fabs(deta),fabs(dphi),1.0/4.0/Nparts_1pt3);
00254                 _h2_S_dphivsdeta_mb_1pt3->Fill(fabs(deta),-fabs(dphi),1.0/4.0/Nparts_1pt3);
00255                 _h2_S_dphivsdeta_mb_1pt3->Fill(-fabs(deta),-fabs(dphi),1.0/4.0/Nparts_1pt3);
00256                 _h2_S_dphivsdeta_mb_1pt3->Fill(fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_1pt3);
00257                 _h2_S_dphivsdeta_mb_1pt3->Fill(-fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_1pt3);
00258               //N>110 (3D - S  1<pt<3)
00259              if(Nparts_04pt>=110){
00260                 _h2_S_dphivsdeta_N110_1pt3->Fill(fabs(deta),fabs(dphi),1.0/4.0/Nparts_1pt3);
00261                 _h2_S_dphivsdeta_N110_1pt3->Fill(-fabs(deta),fabs(dphi),1.0/4.0/Nparts_1pt3);
00262                 _h2_S_dphivsdeta_N110_1pt3->Fill(fabs(deta),-fabs(dphi),1.0/4.0/Nparts_1pt3);
00263                 _h2_S_dphivsdeta_N110_1pt3->Fill(-fabs(deta),-fabs(dphi),1.0/4.0/Nparts_1pt3);
00264                 _h2_S_dphivsdeta_N110_1pt3->Fill(fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_1pt3);
00265                 _h2_S_dphivsdeta_N110_1pt3->Fill(-fabs(deta),2*PI-fabs(dphi),1.0/4.0/Nparts_1pt3);
00266              } //end - N>110                 
00267           }//end - pt cuts
00268         
00269                 
00270       } //end - 2nd particle for the current event
00271 
00272 
00276 
00277       //Background bussiness MB
00278      
00279       if (_oldpartvecMB.size() > 1){ //only if background is already filled
00280        
00281       // Loop old MB event for B-distributions
00282       for (unsigned int ip2 = 0; ip2 < _oldpartvecMB.size(); ip2++) {
00283           const Particle& p2 = _oldpartvecMB[ip2];
00284                   
00285           double pT2 = p2.momentum().pT();              
00286           double eta2 = p2.momentum().eta();
00287           double phi2 = p2.momentum().phi();
00288         
00289           double deta = fabs(eta1-eta2);
00290           double dphi = phi1-phi2;
00291                                                           
00292         if(dphi>PI) dphi=dphi-2*PI;
00293         if(dphi<-PI) dphi=dphi+2*PI;
00294           
00295                 
00296           //MB (3D - B - low pt)
00297           if (pT1 >= 0.1 && pT2 >= 0.1){                  
00298               double pweight=1.0/(Nparts_01pt*oldNpartsMB_01pt);
00299              _h2_B_dphivsdeta_mb_01pt->Fill(deta,fabs(dphi),pweight);
00300              _h2_B_dphivsdeta_mb_01pt->Fill(-deta,fabs(dphi),pweight);
00301              _h2_B_dphivsdeta_mb_01pt->Fill(deta,-fabs(dphi),pweight);
00302              _h2_B_dphivsdeta_mb_01pt->Fill(-deta,-fabs(dphi),pweight);
00303              _h2_B_dphivsdeta_mb_01pt->Fill(deta,2*PI-fabs(dphi),pweight);
00304              _h2_B_dphivsdeta_mb_01pt->Fill(-deta,2*PI-fabs(dphi),pweight);
00305           } //end - pt cuts
00306 
00307           //MB (3D - B - 1<pt<3)
00308           if (pT1 >= 1 && pT1 <= 3 && pT2 >= 1 && pT2 <= 3){ 
00309               double pweight=1.0/(Nparts_1pt3*oldNpartsMB_1pt3);
00310              _h2_B_dphivsdeta_mb_1pt3->Fill(deta,fabs(dphi),pweight);
00311              _h2_B_dphivsdeta_mb_1pt3->Fill(-deta,fabs(dphi),pweight);
00312              _h2_B_dphivsdeta_mb_1pt3->Fill(deta,-fabs(dphi),pweight);
00313              _h2_B_dphivsdeta_mb_1pt3->Fill(-deta,-fabs(dphi),pweight);
00314              _h2_B_dphivsdeta_mb_1pt3->Fill(deta,2*PI-fabs(dphi),pweight);
00315              _h2_B_dphivsdeta_mb_1pt3->Fill(-deta,2*PI-fabs(dphi),pweight);
00316           }//end - pt cuts
00317 
00318 
00319         } //end - particle loop over saved MB particle vector
00320 
00321        } //end - need atleast 1 background event 
00322       
00323       
00324         //Background bussiness N>110
00325         
00326       if (oldpartvecNBin.size() > 100 &&  Nbin==3){ //only if it is already filled
00327       
00328       // Loop old N110 event for B-distributions
00329       for (unsigned int ip2 = 0; ip2 < oldpartvecNBin.size(); ip2++) {
00330           const Particle& p2 = oldpartvecNBin[ip2];
00331          
00332           double pT2 = p2.momentum().pT();              
00333           double eta2 = p2.momentum().eta();
00334           double phi2 = p2.momentum().phi();
00335         
00336           double deta = fabs(eta1-eta2);
00337           double dphi = phi1-phi2;
00338                                                           
00339         if(dphi>PI) dphi=dphi-2*PI;
00340         if(dphi<-PI) dphi=dphi+2*PI;
00341 //        if(dphi>-PI && dphi<-PI/2.) dphi=dphi+2*PI;
00342           
00343         
00344           if (pT1 >= 0.1 && pT2 >= 0.1){                  
00345               //Fill 
00346               double pweight=1.0/(Nparts_01pt*oldNparts_01pt);
00347              _h2_B_dphivsdeta_N110_01pt->Fill(deta,fabs(dphi),pweight);
00348              _h2_B_dphivsdeta_N110_01pt->Fill(-deta,fabs(dphi),pweight);
00349              _h2_B_dphivsdeta_N110_01pt->Fill(deta,-fabs(dphi),pweight);
00350              _h2_B_dphivsdeta_N110_01pt->Fill(-deta,-fabs(dphi),pweight);
00351              _h2_B_dphivsdeta_N110_01pt->Fill(deta,2*PI-fabs(dphi),pweight);
00352              _h2_B_dphivsdeta_N110_01pt->Fill(-deta,2*PI-fabs(dphi),pweight);
00353           } //end - pt cuts
00354 
00355           if (pT1 >= 1 && pT1 <= 3 && pT2 >= 1 && pT2 <= 3){ 
00356               //Fill 
00357               double pweight=1.0/(Nparts_1pt3*oldNparts_1pt3);
00358              _h2_B_dphivsdeta_N110_1pt3->Fill(deta,fabs(dphi),pweight);
00359              _h2_B_dphivsdeta_N110_1pt3->Fill(-deta,fabs(dphi),pweight);
00360              _h2_B_dphivsdeta_N110_1pt3->Fill(deta,-fabs(dphi),pweight);
00361              _h2_B_dphivsdeta_N110_1pt3->Fill(-deta,-fabs(dphi),pweight);
00362              _h2_B_dphivsdeta_N110_1pt3->Fill(deta,2*PI-fabs(dphi),pweight);
00363              _h2_B_dphivsdeta_N110_1pt3->Fill(-deta,2*PI-fabs(dphi),pweight);
00364           }//end - pt cuts
00365 
00366 
00367          } //end - particle loop old N110
00368         
00369         } //end - need atleast 1 event already 
00370 
00371 
00372 
00376       if (oldpartvecNBin.size() > _Nbinslimits[Nbin]){ //only if we already have BG particles                   
00377      
00378       // Loop old event for B-distributions
00379         for (unsigned int ip2 = 0; ip2 < oldpartvecNBin.size(); ip2++) {
00380           const Particle& p2 = oldpartvecNBin[ip2];
00381         
00382           double pT2 = p2.momentum().pT();              
00383           double eta2 = p2.momentum().eta();
00384           double phi2 = p2.momentum().phi();
00385         
00386           double deta = fabs(eta1-eta2);
00387           double dphi = phi1-phi2;
00388                                                           
00389         if(dphi>PI) dphi=dphi-2*PI;
00390         if(dphi<-PI) dphi=dphi+2*PI;
00391          
00392           //loop the pt bins for the DeltaPhi - 1D - Background
00393             for (int iPtbin = 0; iPtbin < 4; iPtbin++) {
00394               if (pT1 > _ptbinslimits[iPtbin] && pT1 < _ptbinslimits[iPtbin+1] &&  pT2 >  _ptbinslimits[iPtbin] && pT2 <  _ptbinslimits[iPtbin+1]){               
00395                    int ibin = iPtbin + Nbin*4;  //which histo to fill, 4x4 matix
00396                    double pweight=1.0/(Nparts_ptbin[iPtbin]*oldNparts_ptbin[iPtbin]);
00397                 _h_B_3D_Nptbins[ibin]->Fill(deta,fabs(dphi),pweight);
00398                 _h_B_3D_Nptbins[ibin]->Fill(-deta,fabs(dphi),pweight);
00399                 _h_B_3D_Nptbins[ibin]->Fill(deta,-fabs(dphi),pweight);
00400                 _h_B_3D_Nptbins[ibin]->Fill(-deta,-fabs(dphi),pweight);
00401                 _h_B_3D_Nptbins[ibin]->Fill(deta,2*PI-fabs(dphi),pweight);
00402                 _h_B_3D_Nptbins[ibin]->Fill(-deta,2*PI-fabs(dphi),pweight);
00403                 if(deta > 2.0 && deta < 4.8){
00404                   _h_B_DeltaPhi_Nptbins[ibin]->Fill(fabs(dphi),pweight);
00405                 }
00406               } //end check pt-bin
00407             } //end - pt-bin loop 
00408           
00409         
00410          } //end - particle loop old N110       
00411         } //end - need atleast 1 event already 
00412       
00413 
00414       } //end - main particle loop for current event    
00415     } //end - if background events found
00416         
00417         
00418 //save the old particle vector
00419     if (Nbin == 0){ _oldpartvec1N35 = ChrgParticles;}  
00420     if (Nbin == 1){ _oldpartvec35N90 = ChrgParticles;}
00421     if (Nbin == 2){ _oldpartvec90N110 = ChrgParticles;}
00422     if (Nbin == 3){ _oldpartvec110N = ChrgParticles; _N110events++;}
00423     if (Nparts_01pt >= 1){ _oldpartvecMB = ChrgParticles;}
00424     
00425     } //end - analyze()
00426     
00427 //AK =====================================================FINALIZE
00428     void finalize() {
00429 
00430         getLog() << Log::INFO << "Number of events after event selection: " << _Nevt_after_cuts << endl;        
00431         getLog() << Log::INFO << "Number of events with N>110:" << _N110events << endl;
00432                         
00433         int nEvent = -99.0;
00434         
00435         nEvent = _hMult_mb_01pt->Integral(3,10000);
00436            _h2_S_dphivsdeta_mb_01pt->Scale(1.0/nEvent);
00437            _h2_B_dphivsdeta_mb_01pt->Scale(_h2_S_dphivsdeta_mb_01pt->Integral()/_h2_B_dphivsdeta_mb_01pt->Integral());
00438            _h2_R_dphivsdeta_mb_01pt->Add((TH2F*)_h2_S_dphivsdeta_mb_01pt);
00439            _h2_R_dphivsdeta_mb_01pt->Add(_h2_B_dphivsdeta_mb_01pt,-1);
00440            _h2_R_dphivsdeta_mb_01pt->Divide(_h2_B_dphivsdeta_mb_01pt);
00441            _h2_R_dphivsdeta_mb_01pt->Scale(_h2_S_dphivsdeta_mb_01pt->Integral());
00442 
00443          nEvent = _hMult_N110_01pt->Integral(3,10000);
00444            _h2_S_dphivsdeta_N110_01pt->Scale(1.0/nEvent);
00445            _h2_B_dphivsdeta_N110_01pt->Scale(_h2_S_dphivsdeta_N110_01pt->Integral()/_h2_B_dphivsdeta_N110_01pt->Integral());
00446            _h2_R_dphivsdeta_N110_01pt->Add(_h2_S_dphivsdeta_N110_01pt);
00447            _h2_R_dphivsdeta_N110_01pt->Add(_h2_B_dphivsdeta_N110_01pt,-1);
00448            _h2_R_dphivsdeta_N110_01pt->Divide(_h2_B_dphivsdeta_N110_01pt);
00449            _h2_R_dphivsdeta_N110_01pt->Scale(_h2_S_dphivsdeta_N110_01pt->Integral());
00450 
00451          nEvent = _hMult_mb_1pt3->Integral(3,10000);
00452            _h2_S_dphivsdeta_mb_1pt3->Scale(1.0/nEvent);
00453            _h2_B_dphivsdeta_mb_1pt3->Scale(_h2_S_dphivsdeta_mb_1pt3->Integral()/_h2_B_dphivsdeta_mb_1pt3->Integral());
00454            _h2_R_dphivsdeta_mb_1pt3->Add(_h2_S_dphivsdeta_mb_1pt3);
00455            _h2_R_dphivsdeta_mb_1pt3->Add(_h2_B_dphivsdeta_mb_1pt3,-1);
00456            _h2_R_dphivsdeta_mb_1pt3->Divide(_h2_B_dphivsdeta_mb_1pt3);
00457            _h2_R_dphivsdeta_mb_1pt3->Scale(_h2_S_dphivsdeta_mb_1pt3->Integral());
00458 
00459          nEvent = _hMult_N110_1pt3->Integral(3,10000);
00460            _h2_S_dphivsdeta_N110_1pt3->Scale(1.0/nEvent);
00461            _h2_B_dphivsdeta_N110_1pt3->Scale(_h2_S_dphivsdeta_N110_1pt3->Integral()/_h2_B_dphivsdeta_N110_1pt3->Integral());
00462            _h2_R_dphivsdeta_N110_1pt3->Add(_h2_S_dphivsdeta_N110_1pt3);
00463            _h2_R_dphivsdeta_N110_1pt3->Add(_h2_B_dphivsdeta_N110_1pt3,-1);
00464            _h2_R_dphivsdeta_N110_1pt3->Divide(_h2_B_dphivsdeta_N110_1pt3);
00465            _h2_R_dphivsdeta_N110_1pt3->Scale(_h2_S_dphivsdeta_N110_1pt3->Integral());
00466 
00467         for (int ibin = 0; ibin < 16; ibin++) {
00468            
00469            int nEvent = _hMult_Nptbins[ibin]->Integral(3,10000);
00470 
00471            _h_S_3D_Nptbins[ibin]->Scale(1.0/nEvent);
00472            _h_B_3D_Nptbins[ibin]->Scale(_h_S_3D_Nptbins[ibin]->Integral()/_h_B_3D_Nptbins[ibin]->Integral());
00473            _h_R_3D_Nptbins[ibin]->Add(_h_S_3D_Nptbins[ibin]);
00474            _h_R_3D_Nptbins[ibin]->Add(_h_B_3D_Nptbins[ibin],-1);
00475            _h_R_3D_Nptbins[ibin]->Divide(_h_B_3D_Nptbins[ibin]);
00476            _h_R_3D_Nptbins[ibin]->Scale(_h_S_3D_Nptbins[ibin]->Integral());
00477 
00478 
00479             _h_S_DeltaPhi_Nptbins[ibin]->Scale(1.0/nEvent);         
00480             // Wei's correction - to get the average multiplicity correct --
00481             _hMult_Nptbins[ibin]->SetAxisRange(2,10000,"X");
00482             double rescale = (_hMult_Nptbins[ibin]->GetMean()-1)/_h_S_DeltaPhi_Nptbins[ibin]->Integral();
00483             _h_S_DeltaPhi_Nptbins[ibin]->Scale(rescale);
00484            //---------------------------------------------------------------
00485            _h_B_DeltaPhi_Nptbins[ibin]->Scale(_h_S_DeltaPhi_Nptbins[ibin]->Integral()/_h_B_DeltaPhi_Nptbins[ibin]->Integral());
00486            _h_R_DeltaPhi_Nptbins[ibin]->Add(_h_S_DeltaPhi_Nptbins[ibin]);
00487            _h_R_DeltaPhi_Nptbins[ibin]->Add(_h_B_DeltaPhi_Nptbins[ibin],-1);
00488            _h_R_DeltaPhi_Nptbins[ibin]->Divide(_h_B_DeltaPhi_Nptbins[ibin]);
00489            _h_R_DeltaPhi_Nptbins[ibin]->Scale(_h_S_DeltaPhi_Nptbins[ibin]->Integral());
00490  
00491 
00492         }
00493 
00494         file -> Write();
00495         
00496     }
00497 
00498 
00499 //AK =====================================================DECLARATIONS
00500   private:
00501 
00502     double detamin;
00503     double detamax;     
00504     int Ndetabins;
00505     
00506     double _N110events;
00507     double _Nevt_after_cuts;
00508 
00509 
00510         TH1F *_hMult;
00511         TH1F *_hPhi_01pt;
00512         TH1F *_h_S_DeltaPhi_Nptbins[16];
00513         TH1F *_h_B_DeltaPhi_Nptbins[16];
00514         TH1F *_h_R_DeltaPhi_Nptbins[16];
00515 
00516         TH2F *_h_S_3D_Nptbins[16];
00517         TH2F *_h_B_3D_Nptbins[16];
00518         TH2F *_h_R_3D_Nptbins[16];
00519         TH1F *_hMult_Nptbins[16];
00520         
00521         TH2F *_h2_S_dphivsdeta_mb_01pt;
00522         TH2F *_h2_B_dphivsdeta_mb_01pt;
00523         TH2F *_h2_R_dphivsdeta_mb_01pt;
00524         TH1F *_hMult_mb_01pt;
00525                 
00526         TH2F *_h2_S_dphivsdeta_mb_1pt3;
00527         TH2F *_h2_B_dphivsdeta_mb_1pt3;
00528         TH2F *_h2_R_dphivsdeta_mb_1pt3;
00529         TH1F *_hMult_mb_1pt3;
00530         
00531         TH2F *_h2_S_dphivsdeta_N110_01pt;
00532         TH2F *_h2_B_dphivsdeta_N110_01pt;
00533         TH2F *_h2_R_dphivsdeta_N110_01pt;
00534         TH1F *_hMult_N110_01pt;
00535         
00536         TH2F *_h2_S_dphivsdeta_N110_1pt3;
00537         TH2F *_h2_B_dphivsdeta_N110_1pt3;       
00538         TH2F *_h2_R_dphivsdeta_N110_1pt3;
00539         TH1F *_hMult_N110_1pt3;
00540         
00541         ParticleVector _oldpartvec1N35;
00542         ParticleVector _oldpartvec35N90;
00543         ParticleVector _oldpartvec90N110; //TODO: this is a bit ugly -> arrays or vec<vec<>>
00544         ParticleVector _oldpartvec110N;
00545         ParticleVector _oldpartvecMB;
00546         
00547         TFile *file;
00548    };
00549 
00550 
00551 
00552   // This global object acts as a hook for the plugin system
00553   AnalysisBuilder<CMS_2010_S8808686> plugin_CMS_2010_S8808686;
00554 
00555 }
00556
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/GeneratorInterface/RivetInterface/src/CMS_2010_S8808686.cc
