Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 #include "Rivet/Analysis.hh"
00003 #include "Rivet/RivetAIDA.hh"
00004 #include "Rivet/Projections/ChargedFinalState.hh"
00005 #include "Rivet/Projections/HadronicFinalState.hh"
00006 #include "Rivet/Tools/ParticleIdUtils.hh"
00007 #include "Rivet/Projections/Beam.hh"
00008 using namespace std;
00009 
00010 namespace Rivet {
00011 
00012 
00013   class CMS_2011_S8884919 : public Analysis {
00014   public:
00015 
00016     CMS_2011_S8884919() : Analysis("CMS_2011_S8884919") {
00017        setBeams(PROTON, PROTON);
00018        setNeedsCrossSection(false);
00019     }
00020     
00021     //Container for the moments
00022    /* struct IMoments {
00023       int n;
00024       double mean;
00025       vector<double> cq;
00026       
00027       //default constructor
00028       IMoments():n(6),mean(0),cq(vector<double>(n,0)){}
00029     };
00030     
00031     IMoments _moments_eta05;
00032     IMoments _moments_eta24; */
00033     
00034     
00035 //AK =====================================================DECLARATIONS
00036   private:
00037 
00038 
00039     vector<AIDA::IHistogram1D*> _h_dNch_dn;
00040     
00041     //vector<double>      _v_dNch_dn_binning1_eta05;
00042     //vector<double>      _v_dNch_dn_binning1_eta24;
00043     
00044     AIDA::IHistogram1D* _h_dNch_dn_pt500_eta24;
00045 
00046     //AIDA::IHistogram1D* _h_cq_eta05;
00047     //AIDA::IHistogram1D* _h_cq_eta24;
00048 
00049     AIDA::IProfile1D*   _h_dmpt_dNch_eta24;
00050 
00051     AIDA::IHistogram1D* _h_KNO_eta05;
00052     AIDA::IHistogram1D* _h_KNO_eta24;
00053         
00054     double _Nevt_after_cuts;
00055     vector<double> _etabins;
00056     vector<int> _nch_in_Evt;
00057     vector<int> _nch_in_Evt_pt500;
00058     double _sumpt;
00059     int nch_max;
00060     
00061     //mandatory functions
00062     void init();
00063     void analyze(const Event&);
00064     void finalize();
00065     
00066     
00067     /*void moments_add(IMoments& , double , double = 1);
00068     void makeMoments(AIDA::IHistogram1D* , IMoments&);
00069     void makeMoments(vector<double>& , IMoments&);
00070     
00071     void makeKNO(AIDA::IHistogram1D* nch , AIDA::IHistogram1D* kno){
00072       
00073     }*/
00074 
00075 
00076   };
00077   
00078 
00079 
00080 //----------------------------------------------------------------------------------------------
00081 
00082 
00083   /*void CMS_2011_S8884919::moments_add(IMoments& moments , double value , double weight){
00084     moments.mean += value*weight;
00085 
00086     for(int m = 0 ; m < moments.n ; ++m)
00087       (moments.cq[m]) += pow(value,m) * weight;
00088   }
00089 
00090 
00091   void CMS_2011_S8884919::makeMoments(AIDA::IHistogram1D* h , IMoments& moments){
00092     int sumHeight = 0;
00093   
00094     for(int i=0 ; i < h->axis().bins() ; ++i){
00095       moments_add(moments , (h->axis().binLowerEdge(i) + h->axis().binUpperEdge(i)) / 2. , h->binHeight(i));
00096       sumHeight += h->binHeight(i);
00097     }
00098       
00099     
00100     //finishing moments
00101     if(sumHeight != 0){
00102       moments.mean /= sumHeight;
00103       for(int m = 0 ; m < moments.n ; ++m)
00104         if( moments.mean != 0) moments.cq[m] = (moments.cq[m] / sumHeight) / pow(moments.mean , m) ;
00105     }
00106   }
00107   
00108   void CMS_2011_S8884919::makeMoments(vector<double>& v , IMoments& moments){
00109     int sumHeight = 0;
00110   
00111     for(unsigned i=0 ; i < v.size() ; ++i){
00112       moments_add(moments , i , v[i]);
00113       sumHeight += v[i];
00114     }
00115       
00116     
00117     //finishing moments
00118     if(sumHeight != 0){
00119       moments.mean /= sumHeight;
00120       for(int m = 0 ; m < moments.n ; ++m)
00121         if( moments.mean != 0) moments.cq[m] = (moments.cq[m] / sumHeight) / pow(moments.mean , m) ;
00122     }
00123   }
00124   */
00125 
00126 
00127 
00128 
00129 //AK =====================================================INIT
00130   void CMS_2011_S8884919::init() {
00131     ChargedFinalState cfs(-2.4, 2.4, 0.0*GeV);
00132     addProjection(cfs, "CFS");
00133     addProjection(Beam(), "Beam");
00134 
00135     nch_max = 400;
00136     _Nevt_after_cuts = 0;
00137 
00138     //eta bins      
00139     _etabins.push_back(0.5) ; _etabins.push_back(1.0) ; _etabins.push_back(1.5) ; _etabins.push_back(2.0) ; _etabins.push_back(2.4) ;
00140     _h_dNch_dn.reserve( _etabins.size() );
00141     ostringstream t("");
00142     
00143     if(fuzzyEquals(sqrtS(), 900, 1E-3)){
00144       
00145       for (unsigned ietabin=0; ietabin < _etabins.size(); ietabin++){
00146         t.str("") ; t << "$|\\eta| <$ " << _etabins[ietabin] << " , $\\sqrt(s)$ = 0.9 TeV" ;
00147         _h_dNch_dn.push_back( bookHistogram1D( 2 + ietabin, 1, 1 , t.str() , "n" , "$P_{n}$") );
00148       }
00149         
00150       _h_dNch_dn_pt500_eta24 = bookHistogram1D( 20 , 1, 1 , "$p_{T} >$ 500 GeV/c , $|\\eta| <$ 2.4 , $\\sqrt(s)$ = 0.9 TeV" , "n" , "$P_{n}$");
00151 
00152       //_h_cq_eta05 = bookHistogram1D( 17 , 1, 1 , "$|\\eta| <$ 0.5 , $\\sqrt(s)$ = 0.9 TeV" , "q" , "$C_{q}$");
00153       //_h_cq_eta24 = bookHistogram1D( 17 , 1, 2 , "$|\\eta| <$ 2.4 , $\\sqrt(s)$ = 0.9 TeV" , "q" , "$C_{q}$");
00154 
00155       _h_dmpt_dNch_eta24 = bookProfile1D( 23 , 1, 1 , "$|\\eta| <$ 2.4 , $\\sqrt(s)$ = 0.9 TeV" , "n" , "$< p_{T}> [ GeV/c ]$");
00156 
00157     }
00158     
00159     if(fuzzyEquals(sqrtS(), 2360, 1E-3)){
00160       for (unsigned ietabin=0; ietabin < _etabins.size(); ietabin++){
00161         t.str("") ; t << "$|\\eta| <$ " << _etabins[ietabin] << " , $\\sqrt(s)$ = 2.36 TeV" ;
00162         _h_dNch_dn.push_back( bookHistogram1D( 7 + ietabin, 1, 1 , t.str() , "n" , "$P_{n}$") );
00163       }
00164         
00165       _h_dNch_dn_pt500_eta24 = bookHistogram1D( 21 , 1, 1 , "$p_{T} >$ 500 GeV/c , $|\\eta| <$ 2.4 , $\\sqrt(s)$ = 2.36 TeV" , "n" , "$P_{n}$");
00166 
00167       //_h_cq_eta05 = bookHistogram1D( 18 , 1, 1 , "$|\\eta| <$ 0.5 , $\\sqrt(s)$ = 2.36 TeV" , "q" , "$C_{q}$");
00168       //_h_cq_eta24 = bookHistogram1D( 18 , 1, 2 , "$|\\eta| <$ 2.4 , $\\sqrt(s)$ = 2.36 TeV" , "q" , "$C_{q}$");
00169 
00170       _h_dmpt_dNch_eta24 = bookProfile1D( 24 , 1, 1 , "$|\\eta| <$ 2.4 , $\\sqrt(s)$ = 2.36 TeV" , "n" , "$< p_{T}> [ GeV/c ]$");
00171     }
00172     
00173     if(fuzzyEquals(sqrtS(), 7000, 1E-3)){
00174       for (unsigned ietabin=0; ietabin < _etabins.size(); ietabin++){
00175         t.str("") ; t << "$|\\eta| <$ " << _etabins[ietabin] << " , $\\sqrt(s)$ = 7 TeV" ;
00176         _h_dNch_dn.push_back( bookHistogram1D( 12 + ietabin, 1, 1 , t.str() , "n" , "$P_{n}$") );
00177       }
00178         
00179       _h_dNch_dn_pt500_eta24 = bookHistogram1D( 22 , 1, 1 , "$p_{T} >$ 500 GeV/c , $|\\eta| <$ 2.4 , $\\sqrt(s)$ = 7 TeV" , "n" , "$P_{n}$");
00180 
00181       //_h_cq_eta05 = bookHistogram1D( 19 , 1, 1 , "$|\\eta| <$ 0.5 , $\\sqrt(s)$ = 7 TeV" , "q" , "$C_{q}$");
00182       //_h_cq_eta24 = bookHistogram1D( 19 , 1, 2 , "$|\\eta| <$ 2.4 , $\\sqrt(s)$ = 7 TeV" , "q" , "$C_{q}$");
00183 
00184       _h_dmpt_dNch_eta24 = bookProfile1D( 25 , 1, 1 , "$|\\eta| <$2.4 , $\\sqrt(s)$ = 7 TeV" , "n" , "$< p_{T}> [ GeV/c ]$");
00185 
00186     }
00187     
00188     //_v_dNch_dn_binning1_eta05.assign(nch_max+1,0);
00189     //_v_dNch_dn_binning1_eta24.assign(nch_max+1,0);    
00190     
00191   }
00192 
00193 //AK =====================================================ANALYZE
00194   void CMS_2011_S8884919::analyze(const Event& event) {
00195     const double weight = event.weight();
00196 
00197     //charge particles
00198     const ChargedFinalState& charged = applyProjection<ChargedFinalState>(event, "CFS");
00199     
00200     //This cut is not needed
00201     /*if (charged.particles().size()<1) {
00202       vetoEvent;
00203     }*/ 
00204     
00205     
00206     _Nevt_after_cuts += weight;
00207     
00208     //resetting the multiplicity for the event to 0;
00209     _nch_in_Evt.assign(_etabins.size() , 0);
00210     _nch_in_Evt_pt500.assign(_etabins.size() , 0);
00211     _sumpt = 0;
00212     
00213     //std::cout << charged.size() << std::endl;
00214 
00215     //Loop over particles in event
00216     foreach (const Particle& p, charged.particles()) {
00217       
00218       //selecting only charged hadrons
00219       if(! PID::isHadron(p.pdgId())) continue;
00220 
00221       double pT = p.momentum().pT();          
00222       double eta = p.momentum().eta();
00223 
00224       _sumpt+=pT;
00225 
00226       //cout << "eta : " << eta << "   pT : " << pT << endl;
00227 
00228       for (int ietabin=_etabins.size()-1; ietabin >= 0 ; --ietabin){   
00229         //cout << "  etabin : " << _etabins[ietabin] << "  res : " << (fabs(eta) <= _etabins[ietabin]) << endl;
00230         if (fabs(eta) <= _etabins[ietabin]){
00231           ++(_nch_in_Evt[ietabin]);
00232           
00233           if(pT>0.5)
00234             ++(_nch_in_Evt_pt500[ietabin]);  
00235         }
00236         //break loop to go faster
00237         else
00238           break;
00239         
00240       }
00241 
00242     }
00243     
00244     //filling mutliplicity dependent histogramms
00245     for (unsigned ietabin=0; ietabin < _etabins.size(); ietabin++){
00246       _h_dNch_dn[ietabin]->fill(_nch_in_Evt[ietabin] , weight);
00247     }
00248     
00249     //Do only if eta bins are the needed ones
00250     if(_etabins[4] == 2.4 && _etabins[0] == 0.5){
00251       if(_nch_in_Evt[4] != 0) _h_dmpt_dNch_eta24->fill(_nch_in_Evt[4] , _sumpt / _nch_in_Evt[4] , weight);
00252       
00253       _h_dNch_dn_pt500_eta24->fill(_nch_in_Evt_pt500[4] , weight);
00254       
00255       /*if( _nch_in_Evt[4] < nch_max ){
00256         (_v_dNch_dn_binning1_eta05[_nch_in_Evt[0]])+=weight;
00257         (_v_dNch_dn_binning1_eta24[_nch_in_Evt[0]])+=weight;
00258       }*/
00259       
00260     }
00261     else
00262       getLog() << Log::WARNING << "You changed the number of eta bins, but forgot to propagate it everywhere !! " << endl;   
00263         
00264   }
00265   
00266 //AK =====================================================FINALIZE
00267   void CMS_2011_S8884919::finalize() {
00268     getLog() << Log::INFO << "Number of events after event selection: " << _Nevt_after_cuts << endl;        
00269     
00270     /*makeMoments(_v_dNch_dn_binning1_eta05 , _moments_eta05);
00271     makeMoments(_v_dNch_dn_binning1_eta24 , _moments_eta24);
00272     
00273     for(int m = 0 ; m < _moments_eta05.n ; ++m){
00274       _h_cq_eta05->fill(m , _moments_eta05.cq[m]);
00275       _h_cq_eta24->fill(m , _moments_eta24.cq[m]);
00276     }*/
00277     
00278     for (unsigned ietabin=0; ietabin < _etabins.size(); ietabin++){
00279       normalize(_h_dNch_dn[ietabin]);
00280     }
00281     normalize(_h_dNch_dn_pt500_eta24);
00282   }
00283 
00284 
00285   // This global object acts as a hook for the plugin system
00286   AnalysisBuilder<CMS_2011_S8884919> plugin_CMS_2011_S8884919;
00287 
00288 }
00289
/data/refman/pasoursint/CMSSW_5_3_0/src/GeneratorInterface/RivetInterface/src/CMS_2011_S8884919.cc
