#include<TH1F>
#include<TProfile>
#include<iostream>
#include<iomanip>
#include"TFile.h"
#include"TTree.h"
#include"TLeaf.h"
#include"TChain.h"
//Functions in this macro///
void Initialize();
////////////////////////////


//Files and chains
TChain* chain;//= new TChain("CaloTowerTree");
TChain* chain2;//= new TChain("hiGoodTightMergedTracksTree");

//Calo Tower Tree
//chain->Add("/hadoop/store/user/jgomez2/ForwardTrees/2010/PanicTime/Forward*");
//Tracks Tree
//chain2->Add("/hadoop/store/user/jgomez2/ForwardTrees/2010/PanicTime/Forward*");


/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
////////////           GLOBAL VARIABLES            //////////////////
/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
Float_t pi=TMath::Pi();
Int_t vterm=1;//Set which order harmonic that this code is meant to measure
Int_t jMax=14;////Set out to which order correction we would like to apply
Int_t NumberOfEvents=0;
//    NumberOfEvents=1;
//NumberOfEvents=2;
NumberOfEvents=10;
//NumberOfEvents=100;
//  NumberOfEvents=50000;
//NumberOfEvents=300000;
//NumberOfEvents=3000000;
//  NumberOfEvents = chain->GetEntries();

const Int_t nCent=5;//Number of Centrality classes

///Looping Variables
Int_t Centrality=0; //This will be the centrality variable later
Int_t NumberOfHits=0;//This will be for both tracks and Hits
Float_t pT=0.;
Float_t phi=0.;
Float_t eta=0.;
Float_t Energy=0.;



Float_t centlo[nCent]={0,10,10,30,30};
Float_t centhi[nCent]={30,20,60,40,60};


//Create the output ROOT file
TFile *myFile;// = new TFile("blah.root","RECREATE");
//TTree *myTree;

//PT Bin Centers
TProfile *PTCenters[nCent];

//EP Resolution Plots

//For Resolution of Standard EP Method
TProfile *HFPMinusHFM[nCent];
TProfile *HFMMiusTrack[nCent];
TProfile *HFPMinusTrack[nCent];


//For Second Order EP
TProfile *SecondHFPMinusHFM[nCent];
TProfile *SecondHFMMiusTrack[nCent];
TProfile *SecondHFPMinusTrack[nCent];

//For Resolution of V1 Even
TProfile *TRPMinusTRM[nCent];
TProfile *TRMMinusTRC[nCent];
TProfile *TRPMinusTRC[nCent];


//Make Subdirectories for what will follow
TDirectory *myPlots;

//Comparing Second order event planes
TH1F *EricWay[nCent];
TH1F *SteveWay[nCent];

//Flattened Vs Unflattened First order EP
TH1F *PsiOne;
TH1F *FlatPsiOne;


//TProfiles to save <pT> and <pT^2> info ....All this is for Ollitrault weights
TProfile *PtStatsWhole[nCent];//Both sides of the tracker
TProfile *PtStatsPos[nCent];//Positive eta tracker
TProfile *PtStatsNeg[nCent];//Negative eta tracker
TProfile *PtStatsMid[nCent];//Middle eta tracker
Float_t ptavwhole[nCent]={0.},pt2avwhole[nCent]={0.};
Float_t ptavpos[nCent]={0.},pt2avpos[nCent]={0.};
Float_t ptavneg[nCent]={0.},pt2avneg[nCent]={0.};
Float_t ptavmid[nCent]={0.},pt2avmid[nCent]={0.};

/////////////////////////
//The following variables are for the FlowVectors Function
/////////////////////////
TProfile *WholeHF[nCent];
TProfile *PosHF[nCent];
TProfile *NegHF[nCent];
TProfile *WholeHFEven[nCent];
//Second Order
TProfile *SecondTracker[nCent];
TProfile *SecondHF[nCent];
TProfile *SecondHFP[nCent];
TProfile *SecondHFM[nCent];
//v1even vectors
TProfile *WholeTracker[nCent];
TProfile *PosTracker[nCent];
TProfile *NegTracker[nCent];
TProfile *MidTracker[nCent];

//Looping Variables
Float_t X_wholehf=0.,Y_wholehf=0.,
  X_wholehfeven=0.,Y_wholehfeven=0.,
  X_poshf=0.,Y_poshf=0.,
  X_neghf=0.,Y_neghf=0.,
  X_2tracker=0.,Y_2tracker=0.,
  X_2hf=0.,Y_2hf=0.,
  X_2hfp=0.,Y_2hfp=0.,
  X_2hfm=0.,Y_2hfm=0.;
//v1 even
Float_t X_wholetracker=0.,Y_wholetracker=0.,
  X_postracker=0.,Y_postracker=0.,
  X_negtracker=0.,Y_negtracker=0.,
  X_midtracker=0.,Y_midtracker=0.;

//Permanent Variables
Float_t Xav_wholehf[nCent]={0.},Yav_wholehf[nCent]={0.},
  Xav_wholehfeven[nCent]={0.},Yav_wholehfeven[nCent]={0.},
  Xav_poshf[nCent]={0.},Yav_poshf[nCent]={0.},
  Xav_neghf[nCent]={0.},Yav_neghf[nCent]={0.},
  Xav_2tracker[nCent]={0.},Yav_2tracker[nCent]={0.},
  Xav_2hf[nCent]={0.},Yav_2hf[nCent]={0.},
  Xav_2hfp[nCent]={0.},Yav_2hfp[nCent]={0.},
  Xav_2hfm[nCent]={0.},Yav_2hfm[nCent]={0.};

//v1 even
Float_t Xav_wholetracker[nCent]={0.},Yav_wholetracker[nCent]={0.},
  Xav_postracker[nCent]={0.},Yav_postracker[nCent]={0.},
  Xav_negtracker[nCent]={0.},Yav_negtracker[nCent]={0.},
  Xav_midtracker[nCent]={0.},Yav_midtracker[nCent]={0.};

//Standard Deviations

Float_t Xstdev_wholehf[nCent]={0.},Ystdev_wholehf[nCent]={0.},
  Xstdev_wholehfeven[nCent]={0.},Ystdev_wholehfeven[nCent]={0.},
  Xstdev_poshf[nCent]={0.},Ystdev_poshf[nCent]={0.},
  Xstdev_neghf[nCent]={0.},Ystdev_neghf[nCent]={0.},
  Xstdev_2tracker[nCent]={0.},Ystdev_2tracker[nCent]={0.},
  Xstdev_2hf[nCent]={0.},Ystdev_2hf[nCent]={0.},
  Xstdev_2hfp[nCent]={0.},Ystdev_2hfp[nCent]={0.},
  Xstdev_2hfm[nCent]={0.},Ystdev_2hfm[nCent]={0.};

//v1 even
Float_t Xstdev_wholetracker[nCent]={0.},Ystdev_wholetracker[nCent]={0.},
  Xstdev_postracker[nCent]={0.},Ystdev_postracker[nCent]={0.},
  Xstdev_negtracker[nCent]={0.},Ystdev_negtracker[nCent]={0.},
  Xstdev_midtracker[nCent]={0.},Ystdev_midtracker[nCent]={0.};



//////////////////////////////////////
// The following variables and plots
// are for the AngularCorrections
// function
///////////////////////////////////////


///Looping Variables
Float_t Xcorr_wholehf=0.,Ycorr_wholehf=0., EPwholehf=0.,
  Xcorr_wholehfeven=0.,Ycorr_wholehfeven=0., EPwholehfeven=0.,
  Xcorr_poshf=0.,Ycorr_wholehf=0., EPposhf=0.,
  Xcorr_neghf=0.,Ycorr_wholehf=0., EPneghf=0.,
  Xcorr_2hf=0.,Ycorr_2hf=0., EP2hf=0.,
  Xcorr_2hfp=0.,Ycorr_2hfp=0., EP2hfp=0.,
  Xcorr_2hfm=0.,Ycorr_2hfm=0., EP2hfm=0.,
  Xcorr_2tracker=0.,Ycorr_2tracker=0., EP2tracker=0.;

//V1 even stuff
Float_t Xcorr_wholetracker=0.,Ycorr_wholetracker=0.,EPwholetracker=0.,
  Xcorr_postracker=0.,Ycorr_postracker=0.,EPpostracker=0.,
  Xcorr_negtracker=0.,Ycorr_negtracker=0.,EPnegtracker=0.,
  Xcorr_midtracker=0.,Ycorr_midtracker=0.,EPmidtracker=0.;

//These Will store the angular correction factors

//Whole HF
TProfile *Coswholehf[nCent];
TProfile *Sinwholehf[nCent];

//Whole HF
TProfile *Coswholehfeven[nCent];
TProfile *Sinwholehfeven[nCent];


//Pos HF
TProfile *Cosposhf[nCent];
TProfile *Sinposhf[nCent];

//Neg HF
TProfile *Cosneghf[nCent];
TProfile *Sinneghf[nCent];

//Second Order HF
TProfile *Cos2hf[nCent];
TProfile *Sin2hf[nCent];

//Second Order HF
TProfile *Cos2hfp[nCent];
TProfile *Sin2hfp[nCent];

//Second Order HF
TProfile *Cos2hfm[nCent];
TProfile *Sin2hfm[nCent];

//Second Order Tracker
TProfile *Cos2tracker[nCent];
TProfile *Sin2tracker[nCent];

//v1 even
//Whole Tracker
TProfile *Coswholetracker[nCent];
TProfile *Sinwholetracker[nCent];

//Pos Tracker
TProfile *Cospostracker[nCent];
TProfile *Sinpostracker[nCent];

//Neg Tracker
TProfile *Cosnegtracker[nCent];
TProfile *Sinnegtracker[nCent];

//Mid Tracker
TProfile *Cosmidtracker[nCent];
TProfile *Sinmidtracker[nCent];


/////////////////////////////////////////
/// Variables that are used in the //////
// Flow Analysis function////////////////
/////////////////////////////////////////

Float_t AngularCorrectionWholeHF=0.,EPfinalwholehf=0.,
  AngularCorrectionWholeHFeven=0.,EPfinalwholehfeven=0.,
  AngularCorrectionPosHF=0.,EPfinalposhf=0.,
  AngularCorrectionSecondHF=0.,EPfinalsecondhf=0.,
  AngularCorrectionNegHF=0.,EPfinalneghf=0.,
  AngularCorrectionSecondHFP=0.,EPfinalsecondhfp=0.,
  AngularCorrectionSecondHFM=0.,EPfinalsecondhfm=0.,

  AngularCorrectionSecondTracker=0.,EPfinalsecondtracker=0.;

//v1 even stuff
Float_t AngularCorrectionWholeTracker=0.,EPfinalwholetracker=0.,
  AngularCorrectionPosTracker=0.,EPfinalpostracker=0.,
  AngularCorrectionNegTracker=0.,EPfinalnegtracker=0.,
  AngularCorrectionMidTracker=0.,EPfinalmidtracker=0.,EPfinalsteveway=0.;


/////////////////
///FLOW PLOTS////
////////////////
TProfile *V1EtaOdd[nCent];
TProfile *V1EtaEven[nCent];
TProfile *V1PtEven[nCent];
TProfile *V1PtOdd[nCent];
TProfile *V1EtaMH[nCent];
TProfile *V1PtMH[nCent];


///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////
/////////////////// END OF GLOBAL VARIABLES ///////////////////////
///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////


//Running the Macro
Int_t EPFinalKolja(){//put functions in here
  Initialize();
  //PTStats();
  //FlowVectors();
  //AngularCorrections();
  //  FlowAnalysis();
  //PrettyPlotting();
    Analyze();
  return 0;
}

void Initialize(){

  //  float eta_bin[11]={-2.5,-2.0,-1.5,-1.0,-0.5,0.0,0.5,1.0,1.5,2.0,2.5};
  float eta_bin_small[7]={-1.5,-1.0,-0.5,0.0,0.5,1.0,1.5};//I will make two eta distributions just in case the results differ somehow

  double pt_bin[17]={0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.4,2.8,3.2,3.6,4.5,6.5,9.5,12};

  chain= new TChain("CaloTowerTree");
  chain2= new TChain("hiGoodTightMergedTracksTree");

  //Calo Tower Tree
  //chain->Add("/hadoop/store/user/jgomez2/ForwardTrees/2010/PanicTime/Forward*");
  //Tracks Tree
  //chain2->Add("/hadoop/store/user/jgomez2/ForwardTrees/2010/PanicTime/Forward*");


  //Create the output ROOT file
  myFile = new TFile("EventPlaneAnalysis.root","recreate");

  //Make Subdirectories for what will follow
  myPlots = myFile->mkdir("Plots");
  myPlots->cd();


  Float_t centlo[nCent]={0,10,10,30,30};
  Float_t centhi[nCent]={30,20,60,40,60};


  char ptcentname[1024];
  char ptcenttitle[1024];

  char res1name[1024],res1title[1024];
  char res2name[1024],res2title[1024];
  char res3name[1024],res3title[1024];
  char res4name[1024],res4title[1024];
  char res5name[1024],res5title[1024];
  char res6name[1024],res6title[1024];
  char res7name[1024],res7title[1024];
  char res8name[1024],res8title[1024];
  char res9name[1024],res9title[1024];

  char restrackname[1024];
  char restracktitle[1024];

  char ptwholename[1024];
  char ptwholetitle[1024];

  char ptposname[1024];
  char ptpostitle[1024];

  char ptnegname[1024];
  char ptnegtitle[1024];

  char ptmidname[1024];
  char ptmidtitle[1024];

  char stevename[1024];
  char stevetitle[1024];

  char ericname[1024];
  char erictitle[1024];

  //Vector Plots
  char wholetrackername[1024],wholetrackertitle[1024];
  char postrackername[1024],postrackertitle[1024];
  char negtrackername[1024],negtrackertitle[1024];
  char midtrackername[1024],midtrackertitle[1024];
  char wholehfname[1024],wholehftitle[1024];
  char wholehfevenname[1024],wholehfeventitle[1024];
  char poshfname[1024],poshftitle[1024];
  char neghfname[1024],neghftitle[1024];
  char secondtrackername[1024],secondtrackertitle[1024];
  char secondhfname[1024],secondhftitle[1024];
  char secondhfpname[1024],secondhfptitle[1024];
  char secondhfmname[1024],secondhfmtitle[1024];

  // <Cos> <Sin> plots
  char coswholetrackername[1024],coswholetrackertitle[1024];
  char cospostrackername[1024],cospostrackertitle[1024];
  char cosnegtrackername[1024],cosnegtrackertitle[1024];
  char cosmidtrackername[1024],cosmidtrackertitle[1024];
  char coswholehfname[1024],coswholehftitle[1024];
  char coswholehfevenname[1024],coswholehfeventitle[1024];
  char cosposhfname[1024],cosposhftitle[1024];
  char cosneghfname[1024],cosneghftitle[1024];
  char cossecondtrackername[1024],cossecondtrackertitle[1024];
  char cossecondhfname[1024],cossecondhftitle[1024];
  char cossecondhfpname[1024],cossecondhfptitle[1024];
  char cossecondhfmname[1024],cossecondhfmtitle[1024];

  char sinwholetrackername[1024],sinwholetrackertitle[1024];
  char sinpostrackername[1024],sinpostrackertitle[1024];
  char sinnegtrackername[1024],sinnegtrackertitle[1024];
  char sinmidtrackername[1024],sinmidtrackertitle[1024];
  char sinwholehfname[1024],sinwholehftitle[1024];
  char sinwholehfevenname[1024],sinwholehfeventitle[1024];
  char sinposhfname[1024],sinposhftitle[1024];
  char sinneghfname[1024],sinneghftitle[1024];
  char sinsecondtrackername[1024],sinsecondtrackertitle[1024];
  char sinsecondhfname[1024],sinsecondhftitle[1024];
  char sinsecondhfpname[1024],sinsecondhfptitle[1024];
  char sinsecondhfmname[1024],sinsecondhfmtitle[1024];

  //V1 Plots
  char v1etaoddname[1024],v1etaoddtitle[1024];
  char v1etaevenname[1024],v1etaeventitle[1024];
  char v1ptoddname[1024],v1ptoddtitle[1024];
  char v1ptevenname[1024],v1pteventitle[1024];
  char v1etamhname[1024],v1etamhtitle[1024];
  char v1ptmhname[1024],v1ptmhtitle[1024];

  for (int i=0;i<nCent;i++)
    {
      //V1 odd eta
      sprintf(v1etaoddname,"V1EtaOdd_%ito%i",centlo[i],centhi[i]);
      sprintf(v1etaoddtitle,"v_{1}^{odd}(#eta) for %i-%i %%",centlo[i],centhi[i]);
      V1EtaOdd[i]= new TProfile(v1etaoddname,v1etaoddtitle,6,eta_bin_small);

      //v1 even eta
      sprintf(v1etaevenname,"V1EtaEven_%ito%i",centlo[i],centhi[i]);
      sprintf(v1etaeventitle,"v_{1}^{even}(#eta) for %i-%i %%",centlo[i],centhi[i]);
      V1EtaEven[i]= new TProfile(v1etaevenname,v1etaeventitle,6,eta_bin_small);

      //v1 pt even
      sprintf(v1ptevenname,"V1PtEven_%ito%i",centlo[i],centhi[i]);
      sprintf(v1pteventitle,"v_{1}^{even}(p_{T}) for %i-%i %%",centlo[i],centhi[i]);
      V1PtEven[i]= new TProfile(v1ptevenname,v1pteventitle,16,pt_bin);


      //v1 pt odd
      sprintf(v1ptoddname,"V1PtOdd_%ito%i",centlo[i],centhi[i]);
      sprintf(v1ptoddtitle,"v_{1}^{odd}(p_{T}) for %i-%i %%",centlo[i],centhi[i]);
      V1PtOdd[i]= new TProfile(v1ptoddname,v1ptoddtitle,16,pt_bin);

      //v1 eta mh
      sprintf(v1etamhname,"V1EtaMH_%ito%i",centlo[i],centhi[i]);
      sprintf(v1etamhtitle,"v_{1}^{MH}(#eta) for %i-%i %%",centlo[i],centhi[i]);
      V1EtaMH[i]= new TProfile(v1etamhname,v1etamhtitle,6,eta_bin_small);

      //v1 pt mh
      sprintf(v1ptmhname,"V1PtMH_%ito%i",centlo[i],centhi[i]);
      sprintf(v1ptmhtitle,"v_{1}^{MH}(p_{T}) for %i-%i %%",centlo[i],centhi[i]);
      V1PtMH[i]= new TProfile(v1ptmhname,v1ptmhtitle,16,pt_bin);


      sprintf(ptcentname,"pTcenter_%ito%i",centlo[i],centhi[i]);
      sprintf(ptcenttitle,"Bin Center for %i-%i %%",centlo[i],centhi[i]);
      PTCenters[i]= new TProfile(ptcentname,ptcenttitle,16,pt_bin); //or can make a TH1f and fill a specific range

      sprintf(res1name,"HFPMINUSHFM_EPResolution_%ito%i",centlo[i],centhi[i]);
      sprintf(res1title,"First Order EP resolution HFPMINUSHFM %i-%i %%",centlo[i],centhi[i]);
      HFPMinusHFM[i]= new TProfile(res1name,res1title,1,0,1);

      sprintf(res2name,"HFMMINUSTrack_EPResolution_%ito%i",centlo[i],centhi[i]);
      sprintf(res2title,"First Order EP resolution HFMMINUSTRACK %i-%i %%",centlo[i],centhi[i]);
      HFMMinusTrack[i]= new TProfile(res2name,res2title,1,0,1);

      sprintf(res3name,"HFPMINUSTrack_EPResolution_%ito%i",centlo[i],centhi[i]);
      sprintf(res3title,"First Order EP resolution HFPMINUSTRACK %i-%i %%",centlo[i],centhi[i]);
      HFPMinusTrack[i]= new TProfile(res3name,res3title,1,0,1);

      //For Second Order EP
      sprintf(res7name,"SecondOrderHFPMINUSHFM_EPResolution_%ito%i",centlo[i],centhi[i]);
      sprintf(res7title,"Second Order EP resolution HFPMINUSHFM %i-%i %%",centlo[i],centhi[i]);
      SecondHFPMinusHFM[i]= new TProfile(res7name,res7title,1,0,1);

      sprintf(res8name,"SecondOrderHFMMINUSTrack_EPResolution_%ito%i",centlo[i],centhi[i]);
      sprintf(res8title,"Second Order EP resolution HFMMINUSTRACK %i-%i %%",centlo[i],centhi[i]);
      SecondHFMMinusTrack[i]= new TProfile(res8name,res8title,1,0,1);

      sprintf(res9name,"2OHFPMINUSTrack_Res_%ito%i",i,i+1);
      sprintf(res9title,"SO resolution HFPMINUSTRACK %i-%i %%",centlo[i],centhi[i]);
      SecondHFPMinusTrack[i]= new TProfile(res9name,res9title,1,0,1);
      
      //For Resolution of V1 Even
      sprintf(res4name,"TRPMinusTRM_EPResolution_%ito%i",centlo[i],centhi[i]);
      sprintf(res4title,"First Order EP resolution TRPMinusTRM %i-%i %%",centlo[i],centhi[i]);
      TRPMinusTRM[i]= new TProfile(res4name,res4title,1,0,1);
      
      sprintf(res5name,"TRMMinusTRC_EPResolution_%ito%i",centlo[i],centhi[i]);
      sprintf(res5title,"First Order EP resolution TRMMinusTRC %i-%i %%",centlo[i],centhi[i]);
      TRMMinusTRC[i]= new TProfile(res5name,res5title,1,0,1);

      sprintf(res6name,"TRPMinusTRC_EPResolution_%ito%i",centlo[i],centhi[i]);
      sprintf(res6title,"First Order EP resolution TRPMinusTRC %i-%i %%",centlo[i],centhi[i]);
      TRPMinusTRC[i]= new TProfile(res6name,res6title,1,0,1);


      //PT stats plots

      //whole tracker
      sprintf(ptwholename,"PtStatsWhole_%ito%i",centlo[i],centhi[i]);
      sprintf(ptwholetitle,"p_{T} stats for whole tracker %i-%i %%",centlo[i],centhi[i]);
      PtStatsWhole[i]= new TProfile(ptwholename,ptwholetitle,2,0,2);
      PtStatsWhole[i]->GetXaxis()->SetBinLabel(1,"<p_{T}>");
      PtStatsWhole[i]->GetXaxis()->SetBinLabel(2,"<p_{T}^{2}>");

      //Pos tracker
      sprintf(ptposname,"PtStatsPos_%ito%i",centlo[i],centhi[i]);
      sprintf(ptpostitle,"p_{T} stats for positive #eta tracker %i-%i %%",centlo[i],centhi[i]);
      PtStatsPos[i]= new TProfile(ptposname,ptpostitle,2,0,2);
      PtStatsPos[i]->GetXaxis()->SetBinLabel(1,"<p_{T}>");
      PtStatsPos[i]->GetXaxis()->SetBinLabel(2,"<p_{T}^{2}>");

      //Neg tracker
      sprintf(ptnegname,"PtStatsNeg_%ito%i",centlo[i],centhi[i]);
      sprintf(ptnegtitle,"p_{T} stats for negative #eta tracker %i-%i %%",centlo[i],centhi[i]);
      PtStatsNeg[i]= new TProfile(ptnegname,ptnegtitle,2,0,2);
      PtStatsNeg[i]->GetXaxis()->SetBinLabel(1,"<p_{T}>");
      PtStatsNeg[i]->GetXaxis()->SetBinLabel(2,"<p_{T}^{2}>");

      //Mid tracker
      sprintf(ptmidname,"PtStatsMid_%ito%i",centlo[i],centhi[i]);
      sprintf(ptmidtitle,"p_{T} stats for mid-rapidity tracker %i-%i %%",centlo[i],centhi[i]);
      PtStatsMid[i]= new TProfile(ptmidname,ptmidtitle,2,0,2);
      PtStatsMid[i]->GetXaxis()->SetBinLabel(1,"<p_{T}>");
      PtStatsMid[i]->GetXaxis()->SetBinLabel(2,"<p_{T}^{2}>");


      //////////////////////
      ////// X,Y Vectors////
      //////////////////////

      //Whole HF
      sprintf(wholehfname,"VecInfo_FirstOrderHF_%ito%i",centlo[i],centhi[i]);
      sprintf(wholehftitle,"VecInfo_FirstOrderHF_%ito%i",centlo[i],centhi[i]);
      WholeHF[i] = new TProfile(wholehfname,wholehftitle,2,0,2);
      WholeHF[i]->GetXaxis()->SetBinLabel(1,"<X>");
      WholeHF[i]->GetXaxis()->SetBinLabel(1,"<Y>");

      //Whole HF
      sprintf(wholehfevenname,"VecInfo_FirstOrderHFEven_%ito%i",centlo[i],centhi[i]);
      sprintf(wholehfeventitle,"VecInfo_FirstOrderHFEven_%ito%i",centlo[i],centhi[i]);
      WholeHFEven[i] = new TProfile(wholehfevenname,wholehfeventitle,2,0,2);
      WholeHFEven[i]->GetXaxis()->SetBinLabel(1,"<X>");
      WholeHFEven[i]->GetXaxis()->SetBinLabel(1,"<Y>");

      //Pos HF
      sprintf(poshfname,"VecInfo_FirstOrderPosHF_%ito%i",centlo[i],centhi[i]);
      sprintf(poshftitle,"VecInfo_FirstOrderPosHF_%ito%i",centlo[i],centhi[i]);
      PosHF[i] = new TProfile(poshfname,poshftitle,2,0,2);
      PosHF[i]->GetXaxis()->SetBinLabel(1,"<X>");
      PosHF[i]->GetXaxis()->SetBinLabel(1,"<Y>");

      //Neg HF
      sprintf(neghfname,"VecInfo_FirstOrderNegHF_%ito%i",centlo[i],centhi[i]);
      sprintf(neghftitle,"VecInfo_FirstOrderNegHF_%ito%i",centlo[i],centhi[i]);
      NegHF[i] = new TProfile(neghfname,neghftitle,2,0,2);
      NegHF[i]->GetXaxis()->SetBinLabel(1,"<X>");
      NegHF[i]->GetXaxis()->SetBinLabel(1,"<Y>");

      //Psi2 HF
      sprintf(secondhfname,"VecInfo_SecondOrderHF_%ito%i",centlo[i],centhi[i]);
      sprintf(secondhftitle,"VecInfo_SecondOrderHF_%ito%i",centlo[i],centhi[i]);
      SecondHF[i] = new TProfile(secondhfname,secondhftitle,2,0,2);
      SecondHF[i]->GetXaxis()->SetBinLabel(1,"<X>");
      SecondHF[i]->GetXaxis()->SetBinLabel(1,"<Y>");

      //Psi2 HFP
      sprintf(secondhfpname,"VecInfo_SecondOrderHFP_%ito%i",centlo[i],centhi[i]);
      sprintf(secondhfptitle,"VecInfo_SecondOrderHFP_%ito%i",centlo[i],centhi[i]);
      SecondHFP[i] = new TProfile(secondhfpname,secondhfptitle,2,0,2);
      SecondHFP[i]->GetXaxis()->SetBinLabel(1,"<X>");
      SecondHFP[i]->GetXaxis()->SetBinLabel(1,"<Y>");

      //Psi2 HFM
      sprintf(secondhfmname,"VecInfo_SecondOrderHFM_%ito%i",centlo[i],centhi[i]);
      sprintf(secondhfmtitle,"VecInfo_SecondOrderHFM_%ito%i",centlo[i],centhi[i]);
      SecondHFM[i] = new TProfile(secondhfmname,secondhfmtitle,2,0,2);
      SecondHFM[i]->GetXaxis()->SetBinLabel(1,"<X>");
      SecondHFM[i]->GetXaxis()->SetBinLabel(1,"<Y>");


      //Psi2 Tracker
      sprintf(secondtrackername,"VecInfo_SecondOrderTracker_%ito%i",centlo[i],centhi[i]);
      sprintf(secondtrackertitle,"VecInfo_SecondOrderTracker_%ito%i",centlo[i],centhi[i]);
      SecondTracker[i] = new TProfile(secondtrackername,secondtrackertitle,2,0,2);
      SecondTracker[i]->GetXaxis()->SetBinLabel(1,"<X>");
      SecondTracker[i]->GetXaxis()->SetBinLabel(1,"<Y>");

      //Whole Tracker
      sprintf(wholetrackername,"VecInfo_FirstOrderWholeTracker_%ito%i",centlo[i],centhi[i]);
      sprintf(wholetrackertitle,"VecInfo_WholeTracker_%ito%i",centlo[i],centhi[i]);
      WholeTracker[i] = new TProfile(wholetrackername,wholetrackertitle,2,0,2);
      WholeTracker[i]->GetXaxis()->SetBinLabel(1,"<X>");
      WholeTracker[i]->GetXaxis()->SetBinLabel(1,"<Y>");

      //Pos Tracker
      sprintf(postrackername,"VecInfo_FirstOrderPosTracker_%ito%i",centlo[i],centhi[i]);
      sprintf(postrackertitle,"VecInfo_FirstOrderPosTracker_%ito%i",centlo[i],centhi[i]);
      PosTracker[i] = new TProfile(postrackername,postrackertitle,2,0,2);
      PosTracker[i]->GetXaxis()->SetBinLabel(1,"<X>");
      PosTracker[i]->GetXaxis()->SetBinLabel(1,"<Y>");

      //Neg Tracker
      sprintf(negtrackername,"VecInfo_FirstOrderNegTracker_%ito%i",centlo[i],centhi[i]);
      sprintf(negtrackertitle,"VecInfo_FirstOrderNegTracker_%ito%i",centlo[i],centhi[i]);
      NegTracker[i] = new TProfile(negtrackername,negtrackertitle,2,0,2);
      NegTracker[i]->GetXaxis()->SetBinLabel(1,"<X>");
      NegTracker[i]->GetXaxis()->SetBinLabel(1,"<Y>");

      //MidTracker
      sprintf(midtrackername,"VecInfo_FirstOrderMidTracker_%ito%i",centlo[i],centhi[i]);
      sprintf(midtrackertitle,"VecInfo_FirstOrderMidTracker_%ito%i",centlo[i],centhi[i]);
      MidTracker[i] = new TProfile(midtrackername,midtrackertitle,2,0,2);
      MidTracker[i]->GetXaxis()->SetBinLabel(1,"<X>");
      MidTracker[i]->GetXaxis()->SetBinLabel(1,"<Y>");

      ///////////////////////////////
      ////////<cos>,<sin> plots//////
      ///////////////////////////////

      //Whole HF
      sprintf(coswholehfname,"CosValues_WholeHF_%ito%i",centlo[i],centhi[i]);
      sprintf(coswholehftitle,"CosValues_WholeHF_%ito%i",centlo[i],centhi[i]);
      Coswholehf[i] = new TProfile(coswholehfname,coswholehftitle,jMax,0,jMax);
      Coswholehf[i]->GetYaxis()->SetTitle("<cos(Xbin*Psi)>");

      sprintf(sinwholehfname,"SinValues_WholeHF_%ito%i",centlo[i],centhi[i]);
      sprintf(sinwholehftitle,"SinValues_WholeHF_%ito%i",centlo[i],centhi[i]);
      Sinwholehf[i] = new TProfile(sinwholehfname,sinwholehftitle,jMax,0,jMax);
      Sinwholehf[i]->GetYaxis()->SetTitle("<sin(Xbin*Psi)>");


      //Whole HF Even
      sprintf(coswholehfevenname,"CosValues_WholeHFEven_%ito%i",centlo[i],centhi[i]);
      sprintf(coswholehfeventitle,"CosValues_WholeHFEven_%ito%i",centlo[i],centhi[i]);
      Coswholehfeven[i] = new TProfile(coswholehfevenname,coswholehfeventitle,jMax,0,jMax);
      Coswholehfeven[i]->GetYaxis()->SetTitle("<cos(Xbin*Psi)>");

      sprintf(sinwholehfevenname,"SinValues_WholeHFEven_%ito%i",centlo[i],centhi[i]);
      sprintf(sinwholehfeventitle,"SinValues_WholeHFEven_%ito%i",centlo[i],centhi[i]);
      Sinwholehfeven[i] = new TProfile(sinwholehfevenname,sinwholehfeventitle,jMax,0,jMax);
      Sinwholehfeven[i]->GetYaxis()->SetTitle("<sin(Xbin*Psi)>");


      //Pos HF
      sprintf(cosposhfname,"CosValues_PosHF_%ito%i",centlo[i],centhi[i]);
      sprintf(cosposhftitle,"CosValues_PosHF_%ito%i",centlo[i],centhi[i]);
      Cosposhf[i] = new TProfile(cosposhfname,cosposhftitle,jMax,0,jMax);
      Cosposhf[i]->GetYaxis()->SetTitle("<cos(Xbin*Psi)>");


      sprintf(sinposhfname,"SinValues_PosHF_%ito%i",centlo[i],centhi[i]);
      sprintf(sinposhftitle,"SinValues_PosHF_%ito%i",centlo[i],centhi[i]);
      Sinposhf[i] = new TProfile(sinposhfname,sinposhftitle,jMax,0,jMax);
      Sinposhf[i]->GetYaxis()->SetTitle("<sin(Xbin*Psi)>");

      //Neg HF
      sprintf(cosneghfname,"CosValues_NegHF_%ito%i",centlo[i],centhi[i]);
      sprintf(cosneghftitle,"CosValues_NegHF_%ito%i",centlo[i],centhi[i]);
      Cosneghf[i] = new TProfile(cosneghfname,cosneghftitle,jMax,0,jMax);
      Cosneghf[i]->GetYaxis()->SetTitle("<cos(Xbin*Psi)>");


      sprintf(sinneghfname,"SinValues_NegHF_%ito%i",centlo[i],centhi[i]);
      sprintf(sinneghftitle,"SinValues_NegHF_%ito%i",centlo[i],centhi[i]);
      Sinneghf[i] = new TProfile(sinneghfname,sinneghftitle,jMax,0,jMax);
      Sinneghf[i]->GetYaxis()->SetTitle("<sin(Xbin*Psi)>");

      //Second Order HF
      sprintf(cossecondhfname,"CosValues_SecondOrderHF_%ito%i",centlo[i],centhi[i]);
      sprintf(cossecondhftitle,"CosValues_SecondOrderHF_%ito%i",centlo[i],centhi[i]);
      Cos2hf[i] = new TProfile(cossecondhfname,cossecondhftitle,jMax,0,jMax);
      Cos2hf[i]->GetYaxis()->SetTitle("<cos(Xbin*Psi)>");


      sprintf(sinsecondhfname,"SinValues_SecondOrderHF_%ito%i",centlo[i],centhi[i]);
      sprintf(sinsecondhftitle,"SinValues_SecondOrderHF_%ito%i",centlo[i],centhi[i]);
      Sin2hf[i] = new TProfile(sinsecondhfname,sinsecondhftitle,jMax,0,jMax);
      Sin2hf[i]->GetYaxis()->SetTitle("<sin(Xbin*Psi)>");

      //Second Order HFP
      sprintf(cossecondhfpname,"CosValues_SecondOrderHFP_%ito%i",centlo[i],centhi[i]);
      sprintf(cossecondhfptitle,"CosValues_SecondOrderHFP_%ito%i",centlo[i],centhi[i]);
      Cos2hfp[i] = new TProfile(cossecondhfpname,cossecondhfptitle,jMax,0,jMax);
      Cos2hfp[i]->GetYaxis()->SetTitle("<cos(Xbin*Psi)>");


      sprintf(sinsecondhfpname,"SinValues_SecondOrderHFP_%ito%i",centlo[i],centhi[i]);
      sprintf(sinsecondhfptitle,"SinValues_SecondOrderHFP_%ito%i",centlo[i],centhi[i]);
      Sin2hfp[i] = new TProfile(sinsecondhfpname,sinsecondhfptitle,jMax,0,jMax);
      Sin2hfp[i]->GetYaxis()->SetTitle("<sin(Xbin*Psi)>");

      //Second Order HFM
      sprintf(cossecondhfmname,"CosValues_SecondOrderHFM_%ito%i",centlo[i],centhi[i]);
      sprintf(cossecondhfmtitle,"CosValues_SecondOrderHFM_%ito%i",centlo[i],centhi[i]);
      Cos2hfm[i] = new TProfile(cossecondhfmname,cossecondhfmtitle,jMax,0,jMax);
      Cos2hfm[i]->GetYaxis()->SetTitle("<cos(Xbin*Psi)>");


      sprintf(sinsecondhfmname,"SinValues_SecondOrderHFM_%ito%i",centlo[i],centhi[i]);
      sprintf(sinsecondhfmtitle,"SinValues_SecondOrderHFM_%ito%i",centlo[i],centhi[i]);
      Sin2hfm[i] = new TProfile(sinsecondhfmname,sinsecondhfmtitle,jMax,0,jMax);
      Sin2hfm[i]->GetYaxis()->SetTitle("<sin(Xbin*Psi)>");

      //Second Order Tracker
      sprintf(cossecondtrackername,"CosValues_SecondOrderTracker_%ito%i",centlo[i],centhi[i]);
      sprintf(cossecondtrackertitle,"CosValues_SecondOrderTracker_%ito%i",centlo[i],centhi[i]);
      Cos2tracker[i] = new TProfile(cossecondtrackername,cossecondtrackertitle,jMax,0,jMax);
      Cos2tracker[i]->GetYaxis()->SetTitle("<cos(Xbin*Psi)>");


      sprintf(sinsecondtrackername,"SinValues_SecondOrderTracker_%ito%i",centlo[i],centhi[i]);
      sprintf(sinsecondtrackertitle,"SinValues_SecondOrderTracker_%ito%i",centlo[i],centhi[i]);
      Sin2tracker[i] = new TProfile(sinsecondtrackername,sinsecondtrackertitle,jMax,0,jMax);
      Sin2tracker[i]->GetYaxis()->SetTitle("<sin(Xbin*Psi)>");

      //v1 even

      //Whole tracker
      sprintf(coswholetrackername,"CosValues_WholeTracker_%ito%i",centlo[i],centhi[i]);
      sprintf(coswholetrackertitle,"CosValues_WholeTracker_%ito%i",centlo[i],centhi[i]);
      Coswholetracker[i] = new TProfile(coswholetrackername,coswholetrackertitle,jMax,0,jMax);
      Coswholetracker[i]->GetYaxis()->SetTitle("<cos(Xbin*Psi)>");


      sprintf(sinwholetrackername,"SinValues_WholeTracker_%ito%i",centlo[i],centhi[i]);
      sprintf(sinwholetrackertitle,"SinValues_WholeTracker_%ito%i",centlo[i],centhi[i]);
      Sinwholetracker[i] = new TProfile(sinwholetrackername,sinwholetrackertitle,jMax,0,jMax);
      Sinwholetracker[i]->GetYaxis()->SetTitle("<sin(Xbin*Psi)>");

      //Pos Tracker
      sprintf(cospostrackername,"CosValues_PosTracker_%ito%i",centlo[i],centhi[i]);
      sprintf(cospostrackertitle,"CosValues_PosTracker_%ito%i",centlo[i],centhi[i]);
      Cospostracker[i] = new TProfile(cospostrackername,cospostrackertitle,jMax,0,jMax);
      Cospostracker[i]->GetYaxis()->SetTitle("<cos(Xbin*Psi)>");


      sprintf(sinpostrackername,"SinValues_PosTracker_%ito%i",centlo[i],centhi[i]);
      sprintf(sinpostrackertitle,"SinValues_PosTracker_%ito%i",centlo[i],centhi[i]);
      Sinpostracker[i] = new TProfile(sinpostrackername,sinpostrackertitle,jMax,0,jMax);
      Sinpostracker[i]->GetYaxis()->SetTitle("<sin(Xbin*Psi)>");

      //Neg Tracker
      sprintf(cosnegtrackername,"CosValues_NegTracker_%ito%i",centlo[i],centhi[i]);
      sprintf(cosnegtrackertitle,"CosValues_NegTracker_%ito%i",centlo[i],centhi[i]);
      Cosnegtracker[i] = new TProfile(cosnegtrackername,cosnegtrackertitle,jMax,0,jMax);
      Cosnegtracker[i]->GetYaxis()->SetTitle("<cos(Xbin*Psi)>");


      sprintf(sinnegtrackername,"SinValues_NegTracker_%ito%i",centlo[i],centhi[i]);
      sprintf(sinnegtrackertitle,"SinValues_NegTracker_%ito%i",centlo[i],centhi[i]);
      Sinnegtracker[i] = new TProfile(sinnegtrackername,sinnegtrackertitle,jMax,0,jMax);
      Sinnegtracker[i]->GetYaxis()->SetTitle("<sin(Xbin*Psi)>");

      //Mid Tracker
      sprintf(cosmidtrackername,"CosValues_MidTracker_%ito%i",centlo[i],centhi[i]);
      sprintf(cosmidtrackertitle,"CosValues_MidTracker_%ito%i",centlo[i],centhi[i]);
      Cosmidtracker[i] = new TProfile(cosmidtrackername,cosmidtrackertitle,jMax,0,jMax);
      Cosmidtracker[i]->GetYaxis()->SetTitle("<cos(Xbin*Psi)>");


      sprintf(sinmidtrackername,"SinValues_MidTracker_%ito%i",centlo[i],centhi[i]);
      sprintf(sinmidtrackertitle,"SinValues_MidTracker_%ito%i",centlo[i],centhi[i]);
      Sinmidtracker[i] = new TProfile(sinmidtrackername,sinmidtrackertitle,jMax,0,jMax);
      Sinmidtracker[i]->GetYaxis()->SetTitle("<sin(Xbin*Psi)>");

      sprintf(stevename,"StevePsi2_%ito%i",centlo[i],centhi[i]);
      sprintf(stevetitle,"StevePsi2_%ito%i",centlo[i],centhi[i]);
      SteveWay[i] = new TH1F(stevename,stevetitle,168,-TMath::Pi(),TMath::Pi());
      SteveWay[i]->GetXaxis()->SetTitle("#Psi_{2}");

      sprintf(ericname,"EricPsi2_%ito%i",centlo[i],centhi[i]);
      sprintf(erictitle,"EricPsi2_%ito%i",centlo[i],centhi[i]);
      EricWay[i] = new TH1F(ericname,erictitle,168,-TMath::Pi(),TMath::Pi());
      EricWay[i]->GetXaxis()->SetTitle("#Psi_{2}");
    }//end of centrality loop

  //Flattened Vs Unflattened First order EP
  PsiOne = new TH1F("FirstOrderNoFlat","#Psi_{1} unflattened dist",180,-TMath::Pi(),TMath::Pi());

  FlatPsiOne = new TH1F("FirstOrderFlat","#Psi_{1} flattened dist",180,-TMath::Pi(),TMath::Pi());

}//end of initialize function


void Analyze(){
  /*  for (int i=0;i<NumberOfEvents;i++)
      {
      // some event cuts
      //get Centrality value


      // Noe do for ALL qualifying centralities
      for (Int_t c=0 ; c<nCent ; ++i ){
      if ( centvalue > centhi[c] ) continue;
      if ( centvalue < centlo[c] ) continue;

      // NOW do stuff and fill resulthisto[c]
      }

      }//end of loop over events
  */
  //  myTree->Fill();
  myFile->Write();
}//end of analyze function