CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h

Go to the documentation of this file.
00001 #ifndef __HiEvtPlaneFlatten__
00002 #define __HiEvtPlaneFlatten__
00003 // -*- C++ -*-
00004 //
00005 // Package:    HiEvtPlaneFlatten
00006 // Class:      HiEvtPlaneFlatten
00007 // 
00008 
00009 //
00010 //
00011 // Original Author:  Stephen Sanders
00012 //         Created:  Mon Jun  7 14:40:12 EDT 2010
00013 // $Id: HiEvtPlaneFlatten.h,v 1.5 2012/02/15 10:33:36 eulisse Exp $
00014 //
00015 //
00016 
00017 // system include files
00018 #include <memory>
00019 #include <iostream>
00020 #include <string>
00021 
00022 // user include files
00023 #include "FWCore/Framework/interface/Frameworkfwd.h"
00024 #include "FWCore/Framework/interface/EDProducer.h"
00025 
00026 #include "FWCore/Framework/interface/Event.h"
00027 #include "FWCore/Framework/interface/MakerMacros.h"
00028 
00029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00030 
00031 #include "DataFormats/Common/interface/Handle.h"
00032 #include "FWCore/Framework/interface/ESHandle.h"
00033 
00034 #include "DataFormats/HeavyIonEvent/interface/EvtPlane.h"
00035 
00036 #include "TMath.h"
00037 #include <vector>
00038 
00039 #define MAXCUT 5000
00040 
00041 //
00042 // class declaration
00043 //
00044 
00045 class HiEvtPlaneFlatten {
00046 public:
00047 
00048   explicit HiEvtPlaneFlatten()
00049   {
00050     pi = TMath::Pi();
00051     hcentbins = 1;
00052     hOrder = 20;
00053     vorder = 2;    //sets order of event plane
00054     minvtx = -25;
00055     delvtx = 5;
00056     nvtxbins = 10;
00057   }
00058   void Init(int order, int ncentbins,const int centbinCompression, std::string tag, int vord)
00059   {
00060     hOrder = order;  //order of flattening
00061     vorder = vord;   //1(v1), 2(v2), 3(v3), 4(v4)       
00062     hcentbins = ncentbins;
00063     centbinComp = centbinCompression;
00064     if(hcentbins<=0) hcentbins = 1;
00065     hbins = hcentbins*nvtxbins*hOrder;
00066     if(hbins>MAXCUT) {
00067       std::cout<<"Too many cuts for flattening calculation.  RESET to deaults"<<std::endl;
00068       hcentbins = 1;
00069       hOrder = 21;
00070     }
00071     for(int i = 0; i<hbins; i++) {
00072       flatX[i]=0;
00073       flatY[i]=0;
00074       flatCnt[i]=0;
00075     } 
00076   }
00077 
00078   int GetCutIndx(int centbin, double vtx, int iord)
00079   {
00080     int cut;
00081     int icent = centbin/centbinComp;
00082     if(icent < 0 || icent > hcentbins) return -1;
00083     int ivtx = (vtx-minvtx)/delvtx;
00084     if(ivtx < 0 || ivtx > nvtxbins) return -1;
00085     cut = hOrder*nvtxbins*icent + hOrder*ivtx + iord;
00086     if(cut<0 || cut>hbins) return -1;
00087     return cut;
00088   }
00089   
00090   void Fill(double psi, double vtx, int centbin)
00091   {
00092     if(fabs(psi)>4 ) return;
00093     for(int k = 0; k<hOrder; k++) {
00094       double fsin = sin(vorder*(k+1)*psi);
00095       double fcos = cos(vorder*(k+1)*psi);
00096       int indx = GetCutIndx(centbin,vtx,k);
00097       if(indx>=0) {
00098         flatX[indx]+=fcos;
00099         flatY[indx]+=fsin;
00100         ++flatCnt[indx];
00101       }
00102     }
00103   }
00104 
00105   double GetFlatPsi(double psi, double vtx, double cent)
00106   {
00107     double correction = 0;
00108     for(int k = 0; k<hOrder; k++) {
00109       int indx = GetCutIndx(cent,vtx,k);
00110       correction+=(2./(double)((k+1)*vorder))*(flatXDB[indx]*sin(vorder*(k+1)*psi)-flatYDB[indx]*cos(vorder*(k+1)*psi));
00111     }
00112     psi+=correction;
00113     psi=bounds(psi);
00114     psi=bounds2(psi);
00115     return psi;
00116   }
00117   
00118   ~HiEvtPlaneFlatten(){}
00119   int GetHBins(){return hbins;}
00120   double GetX(int bin){return flatX[bin];}
00121   double GetY(int bin){return flatY[bin];}
00122   double GetCnt(int bin) {return flatCnt[bin];}
00123   void SetXDB(int indx, double val) {flatXDB[indx]=val;}
00124   void SetYDB(int indx, double val) {flatYDB[indx]=val;}
00125   Double_t bounds(Double_t ang) {
00126     if(ang<-pi) ang+=2.*pi;
00127     if(ang>pi)  ang-=2.*pi;
00128     return ang;
00129   }
00130   Double_t bounds2(Double_t ang) {
00131     double range = pi/(double) vorder;
00132     if(ang<-range) ang+=2*range;
00133     if(ang>range)  ang-=2*range;
00134     return ang;
00135   }
00136 private:
00137   double flatX[MAXCUT];
00138   double flatY[MAXCUT];
00139   double flatXDB[MAXCUT];
00140   double flatYDB[MAXCUT];
00141   double flatCnt[MAXCUT];
00142   int hOrder;    //flattening order
00143   int hcentbins; //# of centrality bins
00144   int centbinComp;
00145   int hbins;
00146   int vorder; //order of flattened event plane
00147   double pi;
00148 
00149   int nvtxbins;
00150   double minvtx;
00151   double delvtx;
00152 
00153 };
00154 
00155 
00156 
00157 #endif