PWGCF/Correlations/macros/jcorran/Util.C

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //Functions used by the other macros
  17 //Author: Jason Glyndwr Ulery, ulery@uni-frankfurt.de
  18
  19 void MakeProjections(Float_t &xTPt1, Float_t &xTPt2, Float_t &xAPt1, Float_t &xAPt2, Int_t xCent, TList *xList, TH1F *xOutHistPhi, TH1F *xOutHistPhiMix, TH1F *xOutHistEtaN, TH1F *xOutHistEtaNMix, TH1F *xOutHistEtaA, TH1F *xOutHistEtaAMix, TH2F *xOutHistPhiEta, TH2F *xOutHistPhiEtaMix, Float_t xMPt[4], Int_t xptweighted=0, Int_t xMC=0, Int_t xNotDelta=0, Int_t xSign=0){//(Trigger Pt Low, Trigger Pt High, Assoc Pt Low, Assoc Pt High, Centrality Bin, List, Hists, mean pt array, Pt weighted?, Monty Carlo?, gives phi from triggered and from all for delta phi and delta phi mixed respectively)
  20   // TFile *xFile=new TFile(xFileName);
  21   //TList *xList=xFile->Get("julery_DiHadron");
  22   //This needs changed if Input binning is changed
  23   //cout << xptweighted << " " << xMC << " " << xNotDelta << " " << xSign << endl;
  24   const int xnpt=12;
  25   //float xptbins[(xnpt+1)]={2.5,3,4,6,8,10,15,20,30,40,50,75,100};
  26   //IF you change also change in EffFit
  27   float xptbins[(xnpt+1)]={2,2.5,3,4,5,6,8,10,15,20,30,40,50};
  28   const int xapt=30;
  29   float xaptbins[(xapt+1)]={0.25,0.5,0.75,1,1.5,2,2.5,3,3.5,4,4.5,5,6,7,8,9,10,12,15,20,25,30,35,40,45,50,60,70,80,90,100};
  30   float xPi=3.1415926535898;
  31
  32   //After this no changes should be needed
  33   float xsum1,xsum2,xerr1,xerr2;
  34   int xtpt1=-1,xtpt2=-1;
  35   TH2F *xtemphist1;
  36   TH2F *xtemphist2;
  37   TH3F *xtemphist3;
  38   TH3F *xtemphist4;
  39   TH2F *xtemphistEN;
  40   TH2F *xtemphistENM;
  41   TH2F *xtemphistEA;
  42   TH2F *xtemphistEAM;
  43   int xapt1,xapt2,xapt3,xapt4;
  44
  45   for(int xi=0;xi<=xnpt;xi++){
  46     if(fabs(xTPt1-xptbins[xi])<0.2)xtpt1=xi;
  47     if(fabs(xTPt2-xptbins[xi])<0.2)xtpt2=xi;
  48   }
  49   //cout << xtpt1 << " " << xtpt2 << endl;
  50   if(xtpt1<0)cout << "Invalid Trigger Pt1" << endl;
  51   if(xtpt2<0)cout << "Invalid Trigger Pt2" << endl;
  52   if(xtpt1<0||xtpt2<0)break;
  53
  54
  55   if(xAPt2>xTPt1){
  56     cout << "Associateds must be less then trigger (Automatically forced)" << endl;
  57     xAPt2=xTPt1;
  58   }
  59   for(int xi=0;xi<=xapt;xi++){
  60     if(fabs(xAPt1-xaptbins[xi])<0.1)xapt1=xi+1;
  61     if(fabs(xAPt2-xaptbins[xi])<0.1)xapt2=xi;
  62   }
  63   // cout << xapt1 << " " << xapt2 << endl;
  64   //since there are some tolerences reset the values to what are actually used
  65   // xAPt1=(xapt1-1)*xasswidth;
  66   // xAPt2=(xapt2-1)*xasswidth;
  67   xAPt1=xaptbins[xapt1-1];
  68   xAPt2=xaptbins[xapt2];
  69   xTPt1=xptbins[xtpt1];
  70   xTPt2=xptbins[xtpt2];
  71   cout << "xTPt: "<< xTPt1 <<  " - " << xTPt2 << endl;
  72   cout << "xAPt: " << xAPt1 << " - " << xAPt2 << endl;
  73   float xntrig=0,xnmix=0,xntrigpt=0;
  74   float exntrig=0,exntrigpt=0;
  75   char xname[150];
  76   char *cpt1[2]={"","Pt"};
  77   char *cdelta2[2]={"DeltaPhiEta","PhiEtaTrig"};
  78   char *cdelta1[2]={"DeltaPhi","PhiTrig"};
  79   char *cmc1[2]={"","_MC"};
  80   char *sign1[3]={"","_LS","_ULS"};
  81   char *sign2[3]={""," LikeSign"," UnlikeSign"};
  82   char *xEtaN[2]={"DeltaEtaN","EtaTrig"};
  83   char *xEtaA[2]={"DeltaEtaA","EtaTrig"};
  84
  85
  86   sprintf(xname,"fHistNTrigger_C%d%s",xCent,cmc1[xMC]);
  87   TH1F *xNTrig=(TH1F*)xList->FindObject(xname);
  88   sprintf(xname,"fHistNMix_C%d%s",xCent,cmc1[xMC]);
  89   TH1F *xNMix=(TH1F*)xList->FindObject(xname);
  90   sprintf(xname,"fHistNTriggerPt_C%d%s",xCent,cmc1[xMC]);
  91   TH1F *xNTrigPt=(TH1F*)xList->FindObject(xname);
  92   sprintf(xname,"fHistMult%s",cmc1[xMC]);
  93   TH1F *xMult=(TH1F*)xList->FindObject(xname);
  94   int xnevent=xMult->GetEntries();
  95
  96   for(int xi=xtpt1;xi<xtpt2;xi++){//Do for all trigger pt ranges we are adding
  97     sprintf(xname,"fHist%s%s_P%d_C%d%s%s",cdelta1[xNotDelta],cpt1[xptweighted],xi,xCent,sign1[xSign],cmc1[xMC]);
  98     xtemphist1=(TH2F*)xList->FindObject(xname);
  99     sprintf(xname,"fHistDeltaPhiMix%s_P%d_C%d%s%s",cpt1[xptweighted],xi,xCent,sign1[xSign],cmc1[xMC]);
 100     if(xNotDelta) sprintf(xname,"fHistPhi%s_C%d%s",cpt1[xptweighted],xCent,cmc1[xMC]);
 101     xtemphist2=(TH2F*)xList->FindObject(xname);
 102
 103     sprintf(xname,"fHist%s%s_P%d_C%d%s%s",xEtaN[xNotDelta],cpt1[xptweighted],xi,xCent,sign1[xSign],cmc1[xMC]);
 104     xtemphistEN=(TH2F*)xList->FindObject(xname);
 105
 106     sprintf(xname,"fHistDeltaEtaNMix%s_P%d_C%d%s%s",cpt1[xptweighted],xi,xCent,sign1[xSign],cmc1[xMC]);
 107     if(xNotDelta) sprintf(xname,"fHistEta%s_C%d%s",cpt1[xptweighted],xCent,cmc1[xMC]);
 108     xtemphistENM=(TH2F*)xList->FindObject(xname);
 109
 110
 111     sprintf(xname,"fHist%s%s_P%d_C%d%s%s",xEtaA[xNotDelta],cpt1[xptweighted],xi,xCent,sign1[xSign],cmc1[xMC]);
 112     xtemphistEA=(TH2F*)xList->FindObject(xname);
 113     sprintf(xname,"fHistDeltaEtaAMix%s_P%d_C%d%s%s",cpt1[xptweighted],xi,xCent,sign1[xSign],cmc1[xMC]);
 114     if(xNotDelta) sprintf(xname,"fHistEta%s_C%d%s",cpt1[xptweighted],xCent,cmc1[xMC]);
 115     xtemphistEAM=(TH2F*)xList->FindObject(xname);
 116
 117     sprintf(xname,"fHist%s%s_P%d_C%d%s",cdelta2[xNotDelta],cpt1[xptweighted],xi,xCent,cmc1[xMC]);
 118     xtemphist3=(TH3F*)xList->FindObject(xname);
 119     sprintf(xname,"fHistDeltaPhiEtaMix%s_P%d_C%d%s",cpt1[xptweighted],xi,xCent,cmc1[xMC]);
 120     if(xNotDelta) sprintf(xname,"fHistPhiEta%s_C%d%s",cpt1[xptweighted],xCent,cmc1[xMC]);
 121     xtemphist4=(TH3F*)xList->FindObject(xname);
 122     //c12355=new TCanvas("c12355","c12355");
 123     // xtemphist4->Draw();
 124
 125     xntrig+=xNTrig->GetBinContent(xi+1);
 126     xnmix+=xNMix->GetBinContent(xi+1);
 127     exntrig+=pow(xNTrig->GetBinError(xi+1),2);
 128     xntrigpt+=xNTrigPt->GetBinContent(xi+1);
 129     exntrigpt+=pow(xNTrigPt->GetBinError(xi+1),2);
 130     if(xi==xtpt1){
 131       TH2F *xDeltaPhi=(TH2F*)xtemphist1->Clone();
 132       xDeltaPhi->SetName("xDeltaPhi");
 133       TH2F *xDeltaPhiMix=(TH2F*)xtemphist2->Clone();
 134       xDeltaPhiMix->SetName("xDeltaPhiMix");
 135       int xnbins=xDeltaPhi->GetNbinsX();
 136       float xmin=xDeltaPhi->GetBinCenter(1)-(xDeltaPhi->GetBinWidth(1))/2.;
 137       float xmax=xDeltaPhi->GetBinCenter(xnbins)+(xDeltaPhi->GetBinWidth(xnbins))/2.;
 138
 139       TH2F *xDeltaEtaN=(TH2F*)xtemphistEN->Clone();
 140       xDeltaEtaN->SetName("xDeltaEtaN");
 141
 142       TH2F *xDeltaEtaNMix=(TH2F*)xtemphistENM->Clone();
 143       xDeltaEtaNMix->SetName("xDeltaEtaNMix");
 144       int xnEbins=xDeltaEtaN->GetNbinsX();
 145       //cout << xnEbins << "EtaBins" << endl;
 146       float xEmin=xDeltaEtaN->GetBinCenter(1)-(xDeltaEtaN->GetBinWidth(1))/2.;
 147       float xEmax=xDeltaEtaN->GetBinCenter(xnEbins)+(xDeltaEtaN->GetBinWidth(xnEbins))/2.;
 148
 149       TH2F *xDeltaEtaA=(TH2F*)xtemphistEA->Clone();
 150       xDeltaEtaA->SetName("xDeltaEtaA");
 151       TH2F *xDeltaEtaAMix=(TH2F*)xtemphistEAM->Clone();
 152       xDeltaEtaAMix->SetName("xDeltaEtaAMix");
 153
 154       TH3F *xDeltaPhiEta=(TH3F*)xtemphist3->Clone();
 155       xDeltaPhiEta->SetName("xDeltaPhiEta");
 156       TH3F *xDeltaPhiEtaMix=(TH3F*)xtemphist4->Clone();
 157       xDeltaPhiEtaMix->SetName("xDeltaPhiEtaMix");
 158       TAxis *xAxis=xDeltaPhiEta->GetXaxis();
 159       TAxis *xAxis2=xDeltaPhiEta->GetXaxis();
 160       TAxis *yAxis2=xDeltaPhiEta->GetYaxis();
 161       int xnbins2=xAxis2->GetNbins();
 162       int ynbins2=yAxis2->GetNbins();
 163       float xmin2=xAxis2->GetBinCenter(1)-(xAxis2->GetBinWidth(1))/2.;
 164       float xmax2=xAxis2->GetBinCenter(xnbins2)+(xAxis2->GetBinWidth(xnbins2))/2.;
 165       float ymin2=yAxis2->GetBinCenter(1)-(yAxis2->GetBinWidth(1))/2.;
 166       float ymax2=yAxis2->GetBinCenter(ynbins2)+(yAxis2->GetBinWidth(ynbins2))/2.;
 167       //cout << "Sum Phi " << xDeltaPhi->GetSum() << endl;
 168     }
 169     else{
 170       for(int xj=1;xj<=xnbins;xj++){
 171         for(int xk=xapt1;xk<=xapt2;xk++){
 172           xDeltaPhi->SetBinContent(xj,xk,(xDeltaPhi->GetBinContent(xj,xk)+xtemphist1->GetBinContent(xj,xk)));
 173           xDeltaPhi->SetBinError(xj,xk,sqrt(pow(xDeltaPhi->GetBinError(xj,xk),2)+pow(xtemphist1->GetBinError(xj,xk),2)));
 174           if(!xNotDelta)xDeltaPhiMix->SetBinContent(xj,xk,(xDeltaPhiMix->GetBinContent(xj,xk)+xtemphist2->GetBinContent(xj,xk)));
 175           if(!xNotDelta)xDeltaPhiMix->SetBinError(xj,xk,sqrt(pow(xDeltaPhiMix->GetBinError(xj,xk),2)+pow(xtemphist2->GetBinError(xj,xk),2)));
 176         }
 177       }//end phi loop
 178
 179       for(int xj=1;xj<=xnEbins;xj++){
 180         for(int xk=xapt1;xk<=xapt2;xk++){
 181           xDeltaEtaN->SetBinContent(xj,xk,(xDeltaEtaN->GetBinContent(xj,xk)+xtemphistEN->GetBinContent(xj,xk)));
 182           xDeltaEtaN->SetBinError(xj,xk,sqrt(pow(xDeltaEtaN->GetBinError(xj,xk),2)+pow(xtemphistEN->GetBinError(xj,xk),2)));
 183           if(!xNotDelta)xDeltaEtaNMix->SetBinContent(xj,xk,(xDeltaEtaNMix->GetBinContent(xj,xk)+xtemphistENM->GetBinContent(xj,xk)));
 184           if(!xNotDelta)xDeltaEtaNMix->SetBinError(xj,xk,sqrt(pow(xDeltaEtaNMix->GetBinError(xj,xk),2)+pow(xtemphistENM->GetBinError(xj,xk),2)));
 185
 186           xDeltaEtaA->SetBinContent(xj,xk,(xDeltaEtaA->GetBinContent(xj,xk)+xtemphistEA->GetBinContent(xj,xk)));
 187           xDeltaEtaA->SetBinError(xj,xk,sqrt(pow(xDeltaEtaA->GetBinError(xj,xk),2)+pow(xtemphistEA->GetBinError(xj,xk),2)));
 188           if(!xNotDelta)xDeltaEtaAMix->SetBinContent(xj,xk,(xDeltaEtaAMix->GetBinContent(xj,xk)+xtemphistEAM->GetBinContent(xj,xk)));
 189           if(!xNotDelta)xDeltaEtaAMix->SetBinError(xj,xk,sqrt(pow(xDeltaEtaAMix->GetBinError(xj,xk),2)+pow(xtemphistEAM->GetBinError(xj,xk),2)));
 190         }
 191       }//end eta loop
 192
 193       for(int xj=1;xj<=xnbins2;xj++){
 194         for(int xk=1;xk<=ynbins2;xk++){
 195           for(int xl=xapt1;xl<=xapt2;xl++){
 196             xDeltaPhiEta->SetBinContent(xj,xk,xl,(xDeltaPhiEta->GetBinContent(xj,xk,xl)+xtemphist3->GetBinContent(xj,xk,xl)));
 197             xDeltaPhiEta->SetBinError(xj,xk,xl,sqrt(pow(xDeltaPhiEta->GetBinError(xj,xk,xl),2)+pow(xtemphist3->GetBinError(xj,xk,xl),2)));
 198             if(!xNotDelta)xDeltaPhiEtaMix->SetBinContent(xj,xk,xl,(xDeltaPhiEtaMix->GetBinContent(xj,xk,xl)+xtemphist4->GetBinContent(xj,xk,xl)));
 199             if(!xNotDelta)xDeltaPhiEtaMix->SetBinError(xj,xk,xl,sqrt(pow(xDeltaPhiEtaMix->GetBinError(xj,xk,xl),2)+pow(xtemphist4->GetBinError(xj,xk,xl),2)));
 200           }
 201         }
 202       }//end eta-phi loop
 203     }//else
 204   }//xi
 205
 206   /*
 207     c124=new TCanvas("c124");
 208     xtemphist4->Draw();
 209     cout << xtemphist4->GetBinContent(10,10,2) << " " << xDeltaPhiEtaMix->GetBinContent(10,10,2) << endl;
 210     cout << xtemphist4->GetBinError(10,10,2) << " " << xDeltaPhiEtaMix->GetBinError(10,10,2) << endl;
 211     cout << xtemphist4->GetBinContent(10,1,2) << " " << xDeltaPhiEtaMix->GetBinContent(10,1,2) << endl;
 212   */
 213
 214
 215   cout << "Number of Triggers: " << xntrig << endl;
 216   cout << "Number of Mixed Events: " << xnmix << endl;
 217   sprintf(xname,"fHistNEvents_C%d",xCent);
 218   TH1F *xNEvents=(TH1F*)xList->FindObject(xname);
 219   cout << "Number of Events: " << xNEvents->GetBinContent(1) << endl;
 220   if(xptweighted)cout << "<P_T Trigger>: " << xntrigpt/xntrig << endl;
 221   //cout << xntrigpt << " " << exntrigpt << " " << xntrig << " " << exntrig << endl;
 222   if(xptweighted)xMPt[0]=xntrigpt/xntrig;
 223   exntrigpt=sqrt(exntrigpt);
 224   exntrig=sqrt(exntrig);
 225   if(xptweighted)xMPt[1]=xMPt[0]*sqrt(pow(exntrigpt/xntrigpt,2)+pow(exntrig/xntrig,2));
 226   xMPt[2]=xntrig;
 227   xMPt[3]=xnmix;
 228   xDeltaPhi->Scale(xnbins/(2.*xPi)/xntrig);
 229   xDeltaEtaN->Scale(xnEbins/(xEmax-xEmin)/xntrig);
 230   xDeltaEtaA->Scale(xnEbins/(xEmax-xEmin)/xntrig);
 231   xDeltaPhiEta->Scale(xnbins2/(2.*xPi)*ynbins2/(ymax2-ymin2)/xntrig);
 232   if(!xNotDelta){
 233     xDeltaPhiMix->Scale(xnbins/(2.*xPi)/xnmix);
 234     xDeltaEtaNMix->Scale(xnEbins/(xEmax-xEmin)/xnmix);
 235     xDeltaEtaAMix->Scale(xnEbins/(xEmax-xEmin)/xnmix);
 236     xDeltaPhiEtaMix->Scale(xnbins2/(2.*xPi)*ynbins2/(ymax2-ymin2)/xnmix);
 237   }
 238   else{
 239     xDeltaPhiMix->Scale(xnbins/(2.*xPi)/xnevent);
 240     xDeltaEtaNMix->Scale(xnEbins/(xEmax-xEmin)/xnevent);
 241     xDeltaEtaAMix->Scale(xnEbins/(xEmax-xEmin)/xnevent);
 242     xDeltaPhiEtaMix->Scale(xnbins2/(2.*xPi)*ynbins2/(ymax2-ymin2)/xnevent);
 243   }
 244
 245   char *tit1="#Delta#phi";
 246   if(xNotDelta)tit1="#phi";
 247   char *tit11="";
 248   if(xptweighted)tit11="p_{T} Weighted";
 249   char *tit12="";
 250   if(xNotDelta)tit12="Triggered";
 251   char *tit2="1";
 252   if(xptweighted)tit2="p_{T}";
 253   char *tit3="Mixed #Delta#phi";
 254   if(xNotDelta)tit3="#phi";
 255
 256   sprintf(xname,"%s %s %s Distribution for %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f%s",tit12,tit11,tit1,xTPt1,xTPt2,xAPt1,xAPt2,sign2[xSign]);
 257   *xOutHistPhi=new TH1F("dPhi",xname,xnbins,xmin,xmax);
 258   xOutHistPhi->Sumw2();
 259   sprintf(xname,"%s (radians)     ",tit1);
 260   xOutHistPhi->GetXaxis()->SetTitle(xname);
 261   sprintf(xname,"#frac{%s}{N_{Trig}}#frac{dN}{d%s}   ",tit2,tit1);
 262   xOutHistPhi->GetYaxis()->SetTitle(xname);
 263   SetTitles1D(xOutHistPhi);
 264   xOutHistPhi->SetMarkerStyle(20);
 265   xOutHistPhi->SetMarkerColor(2);
 266   xOutHistPhi->SetLineColor(2);
 267
 268   sprintf(xname,"%s %s Distribution for %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f%s",tit11,tit3, xTPt1,xTPt2,xAPt1,xAPt2,sign2[xSign]);
 269   *xOutHistPhiMix=new TH1F("dMix",xname,xnbins,xmin,xmax);
 270   xOutHistPhiMix->Sumw2();
 271   sprintf(xname,"%s (radians)     ",tit1);
 272   xOutHistPhiMix->GetXaxis()->SetTitle(xname);
 273   sprintf(xname,"#frac{%s}{N_{Trig}}#frac{dN}{d%s}   ",tit2,tit1);
 274   xOutHistPhiMix->GetYaxis()->SetTitle(xname);
 275   SetTitles1D(xOutHistPhiMix);
 276   xOutHistPhiMix->SetMarkerStyle(21);
 277   xOutHistPhiMix->SetMarkerColor(1);
 278   xOutHistPhiMix->SetLineColor(1);
 279
 280   char *tit14="#Delta#eta";
 281   if(xNotDelta) tit14="#eta";
 282   char *tit15="Mixed #Delta#eta";
 283   if(xNotDelta) tit15="#eta";
 284   char *titNear="Near-Side";
 285   if(xNotDelta)titNear="";
 286
 287   sprintf(xname,"%s %s %s %s Distribution for %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f%s",tit12,titNear,tit11,tit14,xTPt1,xTPt2,xAPt1,xAPt2,sign2[xSign]);
 288   *xOutHistEtaN=new TH1F("dEtaN",xname,xnEbins,xEmin,xEmax);
 289   xOutHistEtaN->Sumw2();
 290   sprintf(xname,"%s     ",tit14);
 291   xOutHistEtaN->GetXaxis()->SetTitle(xname);
 292   sprintf(xname,"#frac{%s}{N_{Trig}}#frac{dN}{d%s}   ",tit2,tit14);
 293   xOutHistEtaN->GetYaxis()->SetTitle(xname);
 294   SetTitles1D(xOutHistEtaN);
 295   xOutHistEtaN->SetMarkerStyle(20);
 296   xOutHistEtaN->SetMarkerColor(2);
 297   xOutHistEtaN->SetLineColor(2);
 298
 299   sprintf(xname,"%s %s %s Distribution for %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f%s",tit11,titNear,tit15,xTPt1,xTPt2,xAPt1,xAPt2,sign2[xSign]);
 300   *xOutHistEtaNMix=new TH1F("dEtaMixN",xname,xnEbins,xEmin,xEmax);
 301   xOutHistEtaNMix->Sumw2();
 302   sprintf(xname,"%s (radians)     ",tit1);
 303   xOutHistEtaNMix->GetXaxis()->SetTitle(xname);
 304   sprintf(xname,"#frac{%s}{N_{Trig}}#frac{dN}{d%s}   ",tit2,tit14);
 305   xOutHistEtaNMix->GetYaxis()->SetTitle(xname);
 306   SetTitles1D(xOutHistEtaNMix);
 307   xOutHistEtaNMix->SetMarkerStyle(21);
 308   xOutHistEtaNMix->SetMarkerColor(1);
 309   xOutHistEtaNMix->SetLineColor(1);
 310
 311   char *titAway="Away-Side";
 312   if(xNotDelta)titAway="";
 313
 314   sprintf(xname,"%s %s %s %s Distribution for %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f%s",tit12,titAway,tit11,tit14,xTPt1,xTPt2,xAPt1,xAPt2,sign2[xSign]);
 315   *xOutHistEtaA=new TH1F("dEtaA",xname,xnEbins,xEmin,xEmax);
 316   xOutHistEtaA->Sumw2();
 317   sprintf(xname,"%s     ",tit14);
 318   xOutHistEtaA->GetXaxis()->SetTitle(xname);
 319   sprintf(xname,"#frac{%s}{N_{Trig}}#frac{dN}{d%s}   ",tit2,tit14);
 320   xOutHistEtaA->GetYaxis()->SetTitle(xname);
 321   SetTitles1D(xOutHistEtaA);
 322   xOutHistEtaA->SetMarkerStyle(20);
 323   xOutHistEtaA->SetMarkerColor(2);
 324   xOutHistEtaA->SetLineColor(2);
 325
 326   sprintf(xname,"%s %s %s Distribution for %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f%s",tit11,titAway,tit15,xTPt1,xTPt2,xAPt1,xAPt2,sign2[xSign]);
 327   *xOutHistEtaAMix=new TH1F("dEtaMixA",xname,xnEbins,xEmin,xEmax);
 328   xOutHistEtaAMix->Sumw2();
 329   sprintf(xname,"%s (radians)     ",tit1);
 330   xOutHistEtaAMix->GetXaxis()->SetTitle(xname);
 331   sprintf(xname,"#frac{%s}{N_{Trig}}#frac{dN}{d%s}   ",tit2,tit14);
 332   xOutHistEtaAMix->GetYaxis()->SetTitle(xname);
 333   SetTitles1D(xOutHistEtaAMix);
 334   xOutHistEtaAMix->SetMarkerStyle(21);
 335   xOutHistEtaAMix->SetMarkerColor(1);
 336   xOutHistEtaAMix->SetLineColor(1);
 337
 338   sprintf(xname,"%s %s %s-%s Distribution for %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f%s",tit12,tit11,tit1,tit14,xTPt1,xTPt2,xAPt1,xAPt2,sign2[xSign]);
 339   *xOutHistPhiEta=new TH2F("dPhiEta",xname,xnbins2,xmin2,xmax2,ynbins2,ymin2,ymax2);
 340   xOutHistPhiEta->Sumw2();
 341   sprintf(xname,"%s (radians)     ",tit1);
 342   xOutHistPhiEta->GetXaxis()->SetTitle(xname);
 343   sprintf(xname,"%s     ",tit14);
 344   xOutHistPhiEta->GetYaxis()->SetTitle(xname);
 345   sprintf(xname,"#frac{%s}{N_{Trig}}#frac{dN}{d%sd%s}   ",tit2,tit1,tit14);
 346   xOutHistPhiEta->GetZaxis()->SetTitle(xname);
 347   SetTitles2D(xOutHistPhiEta);
 348
 349   sprintf(xname,"%s %s-%s Distribution for %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f%s",tit11,tit3,tit14,xTPt1,xTPt2,xAPt1,xAPt2,sign2[xSign]);
 350   *xOutHistPhiEtaMix=new TH2F("dPhiEtaMix",xname,xnbins2,xmin2,xmax2,ynbins2,ymin2,ymax2);
 351   xOutHistPhiEtaMix->Sumw2();
 352   sprintf(xname,"%s (radians)     ",tit1);
 353   xOutHistPhiEtaMix->GetXaxis()->SetTitle(xname);
 354   sprintf(xname,"%s     ",tit14);
 355   xOutHistPhiEtaMix->GetYaxis()->SetTitle(xname);
 356   sprintf(xname,"#frac{%s}{N_{Trig}}#frac{dN}{d%sd%s}   ",tit2,tit1,tit14);
 357   xOutHistPhiEtaMix->GetZaxis()->SetTitle(xname);
 358   SetTitles2D(xOutHistPhiEtaMix);
 359
 360   float xsum1,xsum2,xerr1,xerr2,xsum3,xsum4,xerr3,xerr4;
 361   for(int xi=1;xi<=xnbins;xi++){
 362     xsum1=0; xerr1=0; xsum2=0; xerr2=0;
 363     for(int xj=xapt1; xj<=xapt2;xj++){
 364       xsum1+=xDeltaPhi->GetBinContent(xi,xj);
 365       xerr1+=pow(xDeltaPhi->GetBinError(xi,xj),2);
 366       xsum2+=xDeltaPhiMix->GetBinContent(xi,xj);
 367       xerr2+=pow(xDeltaPhiMix->GetBinError(xi,xj),2);
 368     }
 369
 370     xOutHistPhi->SetBinContent(xi,xsum1);
 371     xOutHistPhi->SetBinError(xi,sqrt(xerr1));
 372     xOutHistPhiMix->SetBinContent(xi,xsum2);
 373     xOutHistPhiMix->SetBinError(xi,sqrt(xerr2));
 374   }//philoop
 375   for(int xi=1;xi<=xnEbins;xi++){
 376     xsum1=0; xerr1=0; xsum2=0; xerr2=0;
 377     xsum3=0; xerr3=0; xsum4=0; xerr4=0;
 378     //cout << "apt" << xapt1 << " " << xapt2 << endl;
 379     for(int xj=xapt1; xj<=xapt2;xj++){
 380       xsum1+=xDeltaEtaN->GetBinContent(xi,xj);
 381       xerr1+=pow(xDeltaEtaN->GetBinError(xi,xj),2);
 382       xsum2+=xDeltaEtaNMix->GetBinContent(xi,xj);
 383       xerr2+=pow(xDeltaEtaNMix->GetBinError(xi,xj),2);
 384       xsum3+=xDeltaEtaA->GetBinContent(xi,xj);
 385       xerr3+=pow(xDeltaEtaA->GetBinError(xi,xj),2);
 386       xsum4+=xDeltaEtaAMix->GetBinContent(xi,xj);
 387       xerr4+=pow(xDeltaEtaAMix->GetBinError(xi,xj),2);
 388     }
 389     xOutHistEtaN->SetBinContent(xi,xsum1);
 390     xOutHistEtaN->SetBinError(xi,sqrt(xerr1));
 391     xOutHistEtaNMix->SetBinContent(xi,xsum2);
 392     xOutHistEtaNMix->SetBinError(xi,sqrt(xerr2));
 393     xOutHistEtaA->SetBinContent(xi,xsum3);
 394     xOutHistEtaA->SetBinError(xi,sqrt(xerr3));
 395     xOutHistEtaAMix->SetBinContent(xi,xsum4);
 396     xOutHistEtaAMix->SetBinError(xi,sqrt(xerr4));
 397   }//etaloop
 398
 399   for(int xi=1;xi<=xnbins2;xi++){
 400     for(int xj=1;xj<=ynbins2;xj++){
 401       xsum1=0; xerr1=0; xsum2=0; xerr2=0;
 402       for(int xk=xapt1; xk<=xapt2;xk++){
 403         xsum1+=xDeltaPhiEta->GetBinContent(xi,xj,xk);
 404         xerr1+=pow(xDeltaPhiEta->GetBinError(xi,xj,xk),2);
 405         xsum2+=xDeltaPhiEtaMix->GetBinContent(xi,xj,xk);
 406         xerr2+=pow(xDeltaPhiEtaMix->GetBinError(xi,xj,xk),2);
 407       }
 408       //cout << xDeltaPhiEtaMix->GetYaxis()->GetBinCenter(xj) << " " << xOutHistPhiEtaMix->GetYaxis()->GetBinCenter(xj) << " " << xsum2 << endl;
 409       xOutHistPhiEta->SetBinContent(xi,xj,xsum1);
 410       xOutHistPhiEta->SetBinError(xi,xj,sqrt(xerr1));
 411       xOutHistPhiEtaMix->SetBinContent(xi,xj,xsum2);
 412       xOutHistPhiEtaMix->SetBinError(xi,xj,sqrt(xsum2));
 413     }
 414   }//phi-eta loop
 415
 416   //c10=new TCanvas("c10","",800,600);
 417   //xOutHistPhi->Draw();
 418   //delete xDeltaPhi;
 419   // delete xDeltaMix;
 420   //delete xtemphist1;
 421   //delete xtemphist2;
 422   //xFile->Close();
 423
 424 }
 425
 426 //--------------------------------------------------------------------------------------------------------
 427
 428 void ZYA1(TH1F *yHistSig, TH1F *yHistBg, float scalea[3], float ZYAMCent=1.,float ZYAMWidth=0.2){
 429   //yerr=0 no error
 430   //yerr=1 no mixed event
 431   //yerr=2 no signal
 432   //yerr=3 no signal and no mixed
 433
 434   //float ZYAMCent=1.3;
 435   //float ZYAMWidth=0.2;//window is cent +/- width so its really a 1/2 width of the window
 436   float a,yea,yerr;
 437   //static float scalea[3];
 438   float ysumsig=0,ysumbg=0,yerrsig=0,yerrbg=0;
 439   float ycent;
 440   int ybins=yHistSig->GetNbinsX();
 441   float yPi=3.1415926535898;
 442   //cout << ybins << endl;
 443   //cout << yHistSig->GetSum() << endl;
 444   //c20=new TCanvas("c20","",800,600);
 445   //yHistSig->Draw();
 446   cout << "ZYAMRegion: " << (ZYAMCent-ZYAMWidth) << "-" << (ZYAMCent+ZYAMWidth) << endl;
 447
 448   for(int yi=1;yi<=ybins;yi++){
 449     ycent=yHistSig->GetBinCenter(yi);
 450     if((fabs(fabs(ycent)-ZYAMCent)<ZYAMWidth)||(fabs(fabs(ycent)-(2*yPi-ZYAMCent))<ZYAMWidth)){
 451       ysumsig+=yHistSig->GetBinContent(yi);
 452       yerrsig+=pow(yHistSig->GetBinError(yi),2);
 453       ysumbg+=yHistBg->GetBinContent(yi);
 454       yerrbg+=pow(yHistBg->GetBinError(yi),2);
 455     }
 456   }
 457   cout << ysumsig << " " << yerrsig << " " << ysumbg << " " << yerrbg << endl;
 458   if(ysumbg==0&&ysumsig==0){
 459     a=0;
 460     yea=0;
 461     yerr=3;
 462   }
 463   else if(ysumbg==0){
 464     a=0;
 465     yea=0;
 466     yerr=1;
 467
 468   }
 469   else if(ysumsig==0){
 470     a=0;
 471     yea=0;
 472     yerr=2;
 473   }
 474   else{
 475     a=ysumsig/ysumbg;
 476     yea=a*sqrt(yerrsig/ysumsig/ysumsig+yerrbg/ysumbg/ysumbg);
 477     yerr=0;
 478   }
 479   cout << a << endl;
 480   scalea[0]=a;
 481   scalea[1]=yea;
 482   scalea[2]=yerr;
 483   if(yerr) cout << "ZYA1 Error:" << yerr << endl;
 484   //return scalea;
 485 }
 486 //----------------------------------------------------------------------------------------------------
 487 void ZYAM2D(TH2F *yHistSig, TH2F *yHistBg, Float_t scalea[3], Int_t nMin=10, Int_t nIter=10){
 488   //for pp nIter=1 might be fine, for heavy-ions where the background has shape it shouldn't
 489
 490   const Int_t NMin=100;
 491
 492   float sumSig, sumBg;
 493   float esumSig, esumBg;
 494   float sumMin, esumMin;
 495   scalea[0]=1;
 496   scalea[1]=0;
 497   float MinArray[NMin][2];
 498   float MinArraySig[NMin][2];
 499   float MinArrayBg[NMin][2];
 500
 501   float maxBinVal;
 502   int maxBin;
 503   float bin, err;
 504   float s,b,es,eb;
 505   for(int i=0;i<nIter;i++){
 506     for(int q=0;q<nMin;q++){
 507       MinArray[q][0]=10000;
 508     }
 509     for(int x=1;x<=yHistSig->GetNbinsX();x++){
 510       for(int y=2;y<=(yHistSig->GetNbinsY()-1);y++){
 511         maxBinVal=-10000;
 512         for(int j=0;j<nMin;j++){//Find highest bin
 513           if(MinArray[j][0]>maxBinVal){
 514             maxBinVal=MinArray[j][0];
 515             maxBin=j;
 516           }
 517         }
 518         s=yHistSig->GetBinContent(x,y);
 519         b=yHistBg->GetBinContent(x,y);
 520         es=yHistSig->GetBinError(x,y);
 521         eb=yHistBg->GetBinError(x,y);
 522         bin=s-scalea[0]*b;
 523         //float a11=scalea[0];
 524         //bin2=s-a11*b;
 525         //cout << "s" << s << " " << "b"  << b << " " << bin << endl;
 526         //cout << bin << " " << maxBinVal << " " << maxBin << endl;
 527         if(bin<maxBinVal){//replce highest bin with current if lower
 528           MinArray[maxBin][0]=bin;
 529           // cout << maxBin << " " << maxBinVal << " " << bin << " " << MinArray[maxBin][0] << endl;
 530           //divide by 0 protections
 531           if(scalea[0]==0)scalea[0]=1E-10;
 532           if(scalea[1]==0)scalea[1]=1E-10;
 533           if(es==0)es=1E-10;
 534           if(eb==0)eb=1E-10;
 535           if(s==0)s=1E-10;
 536           if(b==0)b=1E-10;
 537           MinArray[maxBin][1]=pow(es,2)+pow(scalea[0]*b,2)*(pow(scalea[1]/scalea[0],2)+pow(eb/b,2));
 538           MinArraySig[maxBin][0]=s;
 539           MinArraySig[maxBin][1]=es*es;
 540           MinArrayBg[maxBin][0]=b;
 541           MinArrayBg[maxBin][1]=eb*eb;
 542           //cout << s << " " << b << " " << bin  << " " << maxBin << endl;
 543         }
 544       }
 545     }
 546     sumSig=0, sumBg=0, esumSig=0, esumBg=0, sumMin=0, esumMin=0;
 547     for(int j=0;j<nMin;j++){
 548       //cout << MinArray[j][0] << endl;
 549       if(MinArray[j][0]!=10000&&MinArraySig[j][0]>1E-9&&MinArrayBg[j][0]>1E-9){
 550         //      cout << MinArraySig[j][0] << endl;
 551       sumMin+=MinArray[j][0];
 552       esumMin+=MinArray[j][1];
 553       sumSig+=MinArraySig[j][0];
 554       esumSig+=MinArraySig[j][1];
 555       sumBg+=MinArrayBg[j][0];
 556       esumBg+=MinArrayBg[j][1];
 557       }
 558     }
 559     if(sumSig==0){sumSig=1E-10; cout << "Warning: Zero sumSig" << endl;}
 560     if(sumBg==0){sumBg=1E-10; cout << "Warning: Zero sumBg" << endl;}
 561     scalea[0]=sumSig/sumBg;
 562     //cout << sumSig << endl;
 563     if(sumSig==1E-10)scalea[0]=0;
 564     scalea[1]=scalea[0]*pow(esumSig/sumSig/sumSig+esumBg/sumBg/sumBg,0.5);
 565     esumMin=pow(esumMin,0.5);
 566     esumSig=pow(esumSig,0.5);
 567     esumBg=pow(esumBg,0.5);
 568     cout << "Iter:  " << i << "   Min=" << sumMin << " +/- " << esumMin << "  a=" << scalea[0] << endl;
 569   }
 570 }
 571 //--------------------------------------------------------------------------------------------------
 572 void ZYAM2D2(TH2F *yHistSig, TH2F *yHistBg, float scalea[3], float ZYAMCent=1.,float ZYAMWidth=0.2){
 573   //yerr=0 no error
 574   //yerr=1 no mixed event
 575   //yerr=2 no signal
 576   //yerr=3 no signal and no mixed
 577
 578   //float ZYAMCent=1.3;
 579   //float ZYAMWidth=0.2;//window is cent +/- width so its really a 1/2 width of the window
 580   float a,yea,yerr;
 581   //static float scalea[3];
 582   float ysumsig=0,ysumbg=0,yerrsig=0,yerrbg=0;
 583   float xcent,ycent;
 584   int xbins=yHistSig->GetNbinsX();
 585   int ybins=yHistSig->GetNbinsY();
 586   float yPi=3.1415926535898;
 587   //cout << ybins << endl;
 588   //cout << yHistSig->GetSum() << endl;
 589   //c20=new TCanvas("c20","",800,600);
 590   //yHistSig->Draw();
 591   cout << "ZYAMRegion: " << (ZYAMCent-ZYAMWidth) << "-" << (ZYAMCent+ZYAMWidth) << endl;
 592   for(int xi=1;xi<=xbins;xi++){
 593    xcent=yHistSig->GetBinCenter(xi);
 594     if((fabs(fabs(xcent)-ZYAMCent)<ZYAMWidth)||(fabs(fabs(xcent)-(2*yPi-ZYAMCent))<ZYAMWidth)){
 595     for(int yi=1;yi<=ybins;yi++){
 596       ysumsig+=yHistSig->GetBinContent(xi,yi);
 597       yerrsig+=pow(yHistSig->GetBinError(xi,yi),2);
 598       ysumbg+=yHistBg->GetBinContent(xi,yi);
 599       yerrbg+=pow(yHistBg->GetBinError(xi,yi),2);
 600     }
 601     }
 602   }
 603   cout << ysumsig << " " << yerrsig << " " << ysumbg << " " << yerrbg << endl;
 604   if(ysumbg==0&&ysumsig==0){
 605     a=0;
 606     yea=0;
 607     yerr=3;
 608   }
 609   else if(ysumbg==0){
 610     a=0;
 611     yea=0;
 612     yerr=1;
 613
 614   }
 615   else if(ysumsig==0){
 616     a=0;
 617     yea=0;
 618     yerr=2;
 619   }
 620   else{
 621     a=ysumsig/ysumbg;
 622     yea=a*sqrt(yerrsig/ysumsig/ysumsig+yerrbg/ysumbg/ysumbg);
 623     yerr=0;
 624   }
 625   cout << a << endl;
 626   scalea[0]=a;
 627   scalea[1]=yea;
 628   scalea[2]=yerr;
 629   if(yerr) cout << "ZYA1 Error:" << yerr << endl;
 630   //return scalea;
 631 }
 632 //-------------------------------------------------------------------------------------------------------
 633 void EffCorr(Float_t &yTPt1, Float_t &yTPt2, Float_t &yAPt1, Float_t &yAPt2, Int_t  yCent, char *yFileName, char *yEffFileName, TH1F *yPhiRaw, TH1F *yPhiCorr, TH1F *yPhiEff,  TH1F *yPhiMixRaw, TH1F *yPhiMixCorr, TH1F *yPhiMixEff, TH2F *yPhiEtaRaw, TH2F *yPhiEtaCorr, TH2F *yPhiEtaEff,  TH2F *yPhiEtaMixRaw, TH2F *yPhiEtaMixCorr, TH2F *yPhiEtaMixEff, Float_t yMPt[3], Int_t yptweighted=0, Int_t yNotDelta=0, Int_t yMethod=1){
 634   char name[120], name2[105];
 635
 636   float yMPt2[2];
 637
 638   MakeProjections(yTPt1,yTPt2,yAPt1,yAPt2,yCent,yFileName,yPhiRaw,yPhiMixRaw,yPhiEtaRaw,yPhiEtaMixRaw,yMPt,yptweighted,0,yNotDelta);
 639
 640   sprintf(name,"yPhiRaw_%2.2fPT%2.2f_%2.2fpt%2.2f_%d",yTPt1,yTPt2,yAPt1,yAPt2,yCent);
 641   yPhiRaw->SetName(name);
 642
 643   sprintf(name,"yPhiMixRaw_%2.2fPT%2.2f_%2.2fpt%2.2f_%d",yTPt1,yTPt2,yAPt1,yAPt2,yCent);
 644   yPhiMixRaw->SetName(name);
 645
 646   sprintf(name,"yPhiEtaRaw_%2.2fPT%2.2f_%2.2fpt%2.2f_%d",yTPt1,yTPt2,yAPt1,yAPt2,yCent);
 647   yPhiEtaRaw->SetName(name);
 648
 649   sprintf(name,"yPhiEtaMixRaw_%2.2fPT%2.2f_%2.2fpt%2.2f_%d",yTPt1,yTPt2,yAPt1,yAPt2,yCent);
 650   yPhiEtaMixRaw->SetName(name);
 651
 652
 653   //Projections from MC (from eff file)
 654   TH1F *yPhiMC=new TH1F("yPhiMC","",1,0,1);
 655   TH1F *yPhiMixMC=new TH1F("yPhiMixMC","",1,0,1);
 656   TH2F *yPhiEtaMC=new TH2F("yPhiEtaMC","",1,0,1,1,0,1);
 657   TH2F *yPhiEtaMixMC=new TH2F("yPhiEtaMixMC","",1,0,1,1,0,1);
 658   MakeProjections(yTPt1,yTPt2,yAPt1,yAPt2,yCent,yEffFileName,yPhiMC,yPhiMixMC,yPhiEtaMC,yPhiEtaMixMC,yMPt2,yptweighted,1,yNotDelta);
 659   sprintf(name,"yPhiMC_%2.2fPT%2.2f_%2.2fpt%2.2f_%d",yTPt1,yTPt2,yAPt1,yAPt2,yCent);
 660   yPhiMC->SetName(name);
 661   sprintf(name,"yPhiMixMC_%2.2fPT%2.2f_%2.2fpt%2.2f_%d",yTPt1,yTPt2,yAPt1,yAPt2,yCent);
 662   yPhiMixMC->SetName(name);
 663   sprintf(name,"yPhiEtaMC_%2.2fPT%2.2f_%2.2fpt%2.2f_%d",yTPt1,yTPt2,yAPt1,yAPt2,yCent);
 664   yPhiEtaMC->SetName(name);
 665   sprintf(name,"yPhiEtaMixMC_%2.2fPT%2.2f_%2.2fpt%2.2f_%d",yTPt1,yTPt2,yAPt1,yAPt2,yCent);
 666   yPhiEtaMixMC->SetName(name);
 667
 668   //Efficiency Histograms
 669   MakeProjections(yTPt1,yTPt2,yAPt1,yAPt2,yCent,yEffFileName,yPhiEff,yPhiMixEff,yPhiEtaEff,yPhiEtaMixEff,yMPt2,yptweighted,0,yNotDelta);
 670   sprintf(name,"yPhiEff_%2.2fPT%2.2f_%2.2fpt%2.2f_%d",yTPt1,yTPt2,yAPt1,yAPt2,yCent);
 671   yPhiEff->SetName(name);
 672   sprintf(name2,"Efficiency %s",yPhiEff->GetTitle());
 673   yPhiEff->SetTitle(name2);
 674
 675   sprintf(name,"yPhiMixEff_%2.2fPT%2.2f_%2.2fpt%2.2f_%d",yTPt1,yTPt2,yAPt1,yAPt2,yCent);
 676   yPhiMixEff->SetName(name);
 677   sprintf(name2,"Efficiency %s",yPhiMixEff->GetTitle());
 678   yPhiMixEff->SetTitle(name2);
 679
 680   sprintf(name,"yPhiEtaEff_%2.2fPT%2.2f_%2.2fpt%2.2f_%d",yTPt1,yTPt2,yAPt1,yAPt2,yCent);
 681   yPhiEtaEff->SetName(name);
 682   sprintf(name2,"Efficiency %s",yPhiEtaEff->GetTitle());
 683   yPhiEtaEff->SetTitle(name2);
 684
 685   sprintf(name,"yPhiEtaMixEff_%2.2fPT%2.2f_%2.2fpt%2.2f_%d",yTPt1,yTPt2,yAPt1,yAPt2,yCent);
 686   yPhiEtaMixEff->SetName(name);
 687   sprintf(name2,"Efficiency %s",yPhiEtaMixEff->GetTitle());
 688   yPhiEtaMixEff->SetTitle(name2);
 689
 690
 691   //Calculate Eff
 692   yPhiEff->Divide(yPhiMC);
 693   yPhiMixEff->Divide(yPhiMixMC);
 694   yPhiEtaEff->Divide(yPhiEtaMC);
 695   yPhiEtaMixEff->Divide(yPhiEtaMixMC);
 696
 697   // if(yMethod==0){
 698   // yPhiEff
 699
 700   //Corrected Plots
 701   *yPhiCorr=(TH1F*)yPhiRaw->Clone();
 702   sprintf(name,"yPhiCorr_%2.2fPT%2.2f_%2.2fpt%2.2f_%d",yTPt1,yTPt2,yAPt1,yAPt2,yCent);
 703   yPhiCorr->SetName(name);
 704   yPhiCorr->Divide(yPhiEff);
 705
 706   *yPhiMixCorr=(TH1F*)yPhiMixRaw->Clone();
 707   sprintf(name,"yPhiMixCorr_%2.2fPT%2.2f_%2.2fpt%2.2f_%d",yTPt1,yTPt2,yAPt1,yAPt2,yCent);
 708   yPhiMixCorr->SetName(name);
 709   yPhiMixCorr->Divide(yPhiMixEff);
 710
 711   *yPhiEtaCorr=(TH2F*)yPhiEtaRaw->Clone();
 712   sprintf(name,"yPhiEtaCorr_%2.2fPT%2.2f_%2.2fpt%2.2f_%d",yTPt1,yTPt2,yAPt1,yAPt2,yCent);
 713   yPhiEtaCorr->SetName(name);
 714   yPhiEtaCorr->Divide(yPhiEtaEff);
 715
 716   *yPhiEtaMixCorr=(TH2F*)yPhiEtaMixRaw->Clone();
 717   sprintf(name,"yPhiMixCorr_%2.2fPT%2.2f_%2.2fpt%2.2f_%d",yTPt1,yTPt2,yAPt1,yAPt2,yCent);
 718   yPhiEtaMixCorr->SetName(name);
 719   yPhiEtaMixCorr->Divide(yPhiEtaMixEff);
 720
 721   //Change some titles
 722   sprintf(name2,"Uncorrected %s",yPhiRaw->GetTitle());
 723   yPhiRaw->SetTitle(name2);
 724   sprintf(name2,"Uncorrected %s",yPhiMixRaw->GetTitle());
 725   yPhiMixRaw->SetTitle(name2);
 726   sprintf(name2,"Uncorrected %s",yPhiEtaRaw->GetTitle());
 727   yPhiEtaRaw->SetTitle(name2);
 728   sprintf(name2,"Uncorrected %s",yPhiEtaMixRaw->GetTitle());
 729   yPhiEtaMixRaw->SetTitle(name2);
 730
 731   //Geometry Correction for the eta-phi plots
 732   if(!yNotDelta){
 733     int xbins=yPhiEtaCorr->GetNbinsX();
 734     int ybins=yPhiEtaCorr->GetNbinsY();
 735     float corr;
 736     float center;
 737     float scale=0;
 738     float sumy;
 739     //uses the mc for the geomery correction
 740     for(int yx=1;yx<=xbins;yx++){
 741       scale+=yPhiEtaMixMC->GetBinContent(yx,ybins/2);
 742       scale+=yPhiEtaMixMC->GetBinContent(yx,ybins/2+1);
 743     }
 744     scale/=2;
 745     for(int yy=1;yy<=ybins;yy++){
 746       sumy=0;
 747       for(int yx=1;yx<=xbins;yx++){
 748         sumy+=yPhiEtaMixMC->GetBinContent(yx,yy);
 749       }
 750       if(sumy==0)sumy=0.00001;
 751       corr=scale/sumy;
 752       for(int yx=1;yx<=xbins;yx++){
 753         yPhiEtaCorr->SetBinContent(yx,yy,corr*yPhiEtaCorr->GetBinContent(yx,yy));
 754         yPhiEtaCorr->SetBinError(yx,yy,corr*yPhiEtaCorr->GetBinError(yx,yy));
 755         yPhiEtaMixCorr->SetBinContent(yx,yy,corr*yPhiEtaMixCorr->GetBinContent(yx,yy));
 756         yPhiEtaMixCorr->SetBinError(yx,yy,corr*yPhiEtaMixCorr->GetBinError(yx,yy));
 757       }
 758     }
 759   }
 760   /*
 761   //uses equation for the geometry correct, must be changed if trigger or associtated eta cuts change
 762   TAxis *yAxis=yPhiEtaCorr->GetYaxis();
 763   for(int yx=1;yx<=xbins;yx++){
 764   for(int yy=1;yy<=xbins;yy++){
 765   center=yAxis->GetBinCenter(yy);
 766   if(center<0)corr=1/(1+center/1.8);
 767   else corr=1/(1-center/1.8);
 768   yPhiEtaCorr->SetBinContent(yx,yy,corr*yPhiEtaCorr->GetBinContent(yx,yy));
 769   yPhiEtaCorr->SetBinError(yx,yy,corr*yPhiEtaCorr->GetBinError(yx,yy));
 770   yPhiEtaMixCorr->SetBinContent(yx,yy,corr*yPhiEtaMixCorr->GetBinContent(yx,yy));
 771   yPhiEtaMixCorr->SetBinError(yx,yy,corr*yPhiEtaMixCorr->GetBinError(yx,yy));
 772   }
 773   }
 774   */
 775 }
 776 //----------------------------------------------------------------------------------------------------------
 777 void EffCorr2(Float_t yTPt1, Float_t yTPt2, Float_t yAPt1, Float_t yAPt2, Int_t yCent, TH1F *yPhiEff, TH1F *yPhiMC, TH1F *yPhiMixEff, TH1F *yPhiMixMC, TH1F *yEtaNEff, TH1F *yEtaNMC, TH1F *yEtaNMixEff, TH1F *yEtaNMixMC, TH1F *yEtaAEff, TH1F *yEtaAMC, TH1F *yEtaAMixEff, TH1F *yEtaAMixMC, TH2F *yPhiEtaEff,  TH2F *yPhiEtaMC, TH2F *yPhiEtaMixEff, TH2F *yPhiEtaMixMC, Int_t yMethod=1){//MovingMax is special for yMethod==4
 778   char effname[120], effname2[105];
 779   if(yMethod==0){
 780     for(int i=1;i<=yPhiEff->GetNbinsX();i++){
 781       yPhiEff->SetBinContent(i,1);
 782       yPhiEff->SetBinError(i,0);
 783       yPhiMixEff->SetBinContent(i,1);
 784       yPhiMixEff->SetBinError(i,0);
 785     }
 786     for(int i=1;i<=yEtaNEff->GetNbinsX();i++){
 787       yEtaNEff->SetBinContent(i,1);
 788       yEtaNEff->SetBinError(i,0);
 789       yEtaAEff->SetBinContent(i,1);
 790       yEtaAEff->SetBinError(i,0);
 791       yEtaNMixEff->SetBinContent(i,1);
 792       yEtaNMixEff->SetBinError(i,0);
 793       yEtaAMixEff->SetBinContent(i,1);
 794       yEtaAMixEff->SetBinError(i,0);
 795     }
 796     for(int i=1;i<=yPhiEtaEff->GetXaxis()->GetNbins();i++){
 797       for(int j=1;j<=yPhiEtaEff->GetYaxis()->GetNbins();j++){
 798         yPhiEtaEff->SetBinContent(i,j,1);
 799         yPhiEtaEff->SetBinError(i,j,0);
 800         yPhiEtaMixEff->SetBinContent(i,j,1);
 801         yPhiEtaMixEff->SetBinError(i,j,0);
 802       }
 803     }
 804   }//yMethod==0
 805   else{
 806     //Calculate Eff (Method 1)
 807     yPhiEff->Divide(yPhiMC);
 808     yPhiMixEff->Divide(yPhiMixMC);
 809     yEtaNEff->Divide(yEtaNMC);
 810     yEtaNMixEff->Divide(yEtaNMixMC);
 811     yEtaAEff->Divide(yEtaAMC);
 812     yEtaAMixEff->Divide(yEtaAMixMC);
 813     yPhiEtaEff->Divide(yPhiEtaMC);
 814     yPhiEtaMixEff->Divide(yPhiEtaMixMC);
 815
 816     if(yMethod==2){//Use mixed event efficiencys for triggered events
 817       /*
 818         yPhiEff=(TH1F*)yPhiMixEff->Clone();
 819         yPhiEff->SetName("yPhiEff");
 820       yEtaNEff=(TH1F*)yEtaNMixEff->Clone();
 821       yEtaNEff->SetName("yEtaNEff");
 822       yEtaAEff=(TH1F*)yEtaAMixEff->Clone();
 823       yEtaAEff->SetName("yEtaAEff");
 824       yPhiEtaEff=(TH1F*)yPhiEtaMixEff->Clone();
 825       yPhiEtaEff->SetName("yPhiEtaEff");
 826       */
 827       for(int i=1;i<=yPhiEff->GetNbinsX();i++){
 828         yPhiEff->SetBinContent(i,yPhiMixEff->GetBinContent(i));
 829         yPhiEff->SetBinError(i,yPhiMixEff->GetBinError(i));
 830       }
 831       for(int i=1;i<=yEtaNEff->GetNbinsX();i++){
 832         yEtaNEff->SetBinContent(i,yEtaNMixEff->GetBinContent(i));
 833         yEtaNEff->SetBinError(i,yEtaNMixEff->GetBinError(i));
 834         yEtaAEff->SetBinContent(i,yEtaAMixEff->GetBinContent(i));
 835         yEtaAEff->SetBinError(i,yEtaAMixEff->GetBinError(i));
 836       }
 837       for(int i=1;i<=yPhiEtaEff->GetXaxis()->GetNbins();i++){
 838         for(int j=1;j<=yPhiEtaEff->GetYaxis()->GetNbins();j++){
 839           yPhiEtaEff->SetBinContent(i,j,yPhiEtaMixEff->GetBinContent(i,j));
 840           yPhiEtaEff->SetBinError(i,j,yPhiEtaMixEff->GetBinError(i,j));
 841         }
 842       }
 843     }
 844     if(yMethod==3){//Sum one distribution and use it for all
 845       float eff=0, effE=0, effM=0, effME=0;
 846       for(int i=1;i<=yPhiEff->GetNbinsX();i++){
 847         eff+=yPhiEff->GetBinContent(i);
 848         effE+=pow(yPhiEff->GetBinError(i),2);
 849         effM+=yPhiMixEff->GetBinContent(i);
 850         effME+=pow(yPhiMixEff->GetBinError(i),2);
 851       }
 852       eff/=yPhiEff->GetNbinsX();
 853       effM/=yPhiEff->GetNbinsX();
 854       effE=sqrt(effE)/yPhiEff->GetNbinsX();
 855       effME=sqrt(effME)/yPhiEff->GetNbinsX();
 856       cout << "Triggered Event Eff: " << eff << " +/- " << effE << endl;
 857       cout << "Mixed Event Efficie: " << effM << " +/- " << effME << endl;
 858
 859       for(int i=1;i<=yPhiEff->GetNbinsX();i++){
 860         yPhiEff->SetBinContent(i,eff);
 861         yPhiEff->SetBinError(i,effE);
 862         yPhiMixEff->SetBinContent(i,effM);
 863         yPhiMixEff->SetBinError(i,effME);
 864       }
 865       for(int i=1;i<=yEtaNEff->GetNbinsX();i++){
 866         yEtaNEff->SetBinContent(i,eff);
 867         yEtaNEff->SetBinError(i,effE);
 868         yEtaNMixEff->SetBinContent(i,effM);
 869         yEtaNMixEff->SetBinError(i,effME);
 870         yEtaAEff->SetBinContent(i,eff);
 871         yEtaAEff->SetBinError(i,effE);
 872         yEtaAMixEff->SetBinContent(i,effM);
 873         yEtaAMixEff->SetBinError(i,effME);
 874       }
 875       for(int i=1;i<=yPhiEtaEff->GetXaxis()->GetNbins();i++){
 876         for(int j=1;j<=yPhiEtaEff->GetYaxis()->GetNbins();j++){
 877           yPhiEtaEff->SetBinContent(i,j,eff);
 878           yPhiEtaEff->SetBinError(i,j,effE);
 879           yPhiEtaMixEff->SetBinContent(i,j,effM);
 880           yPhiEtaMixEff->SetBinError(i,j,effME);
 881         }
 882       }
 883     }
 884   }
 885
 886
 887   yPhiEff->GetYaxis()->SetTitle("Efficiency+Contamination");
 888   yPhiMixEff->GetYaxis()->SetTitle("Efficiency+Contamination");
 889   yEtaNEff->GetYaxis()->SetTitle("Efficiency+Contamination");
 890   yEtaNMixEff->GetYaxis()->SetTitle("Efficiency+Contamination");
 891   yEtaAEff->GetYaxis()->SetTitle("Efficiency+Contamination");
 892   yEtaAMixEff->GetYaxis()->SetTitle("Efficiency+Contamination");
 893   yPhiEtaEff->GetZaxis()->SetTitle("Efficiency+Contamination");
 894   yPhiEtaMixEff->GetZaxis()->SetTitle("Efficiency+Contamination");
 895   sprintf(effname,"Efficiency Triggered %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f",yTPt1,yTPt2,yAPt1,yAPt2);
 896   yPhiEff->SetTitle(effname);
 897   sprintf(effname,"Efficiency Mixed %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f", yTPt1,yTPt2,yAPt1,yAPt2);
 898   yPhiMixEff->SetTitle(effname);
 899
 900   sprintf(effname,"Near-Side Efficiency Triggered %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f",yTPt1,yTPt2,yAPt1,yAPt2);
 901   yEtaNEff->SetTitle(effname);
 902   sprintf(effname,"Near-Side Efficiency Mixed %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f", yTPt1,yTPt2,yAPt1,yAPt2);
 903   yEtaNMixEff->SetTitle(effname);
 904
 905   sprintf(effname,"Away-Side Efficiency Triggered %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f",yTPt1,yTPt2,yAPt1,yAPt2);
 906   yEtaAEff->SetTitle(effname);
 907   sprintf(effname,"Away-Side Efficiency Mixed %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f", yTPt1,yTPt2,yAPt1,yAPt2);
 908   yEtaAMixEff->SetTitle(effname);
 909
 910   sprintf(effname,"Efficiency Triggered %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f",yTPt1,yTPt2,yAPt1,yAPt2);
 911   yPhiEtaEff->SetTitle(effname);
 912   sprintf(effname,"Efficiency Mixed %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f",yTPt1,yTPt2,yAPt1,yAPt2);
 913   yPhiEtaMixEff->SetTitle(effname);
 914
 915
 916   //Corrected Plots
 917   //if still problems try shifting this to the other code
 918   // yPhiCorr->Divide(yPhiEff);
 919   /*
 920   //Change some titles
 921   sprintf(effname2,"Uncorrected %s",yPhiRaw->GetTitle());
 922   yPhiRaw->SetTitle(effname2);
 923   sprintf(effname2,"Uncorrected %s",yPhiMixRaw->GetTitle());
 924   yPhiMixRaw->SetTitle(effname2);
 925   sprintf(effname2,"Uncorrected %s",yPhiEtaRaw->GetTitle());
 926   yPhiEtaRaw->SetTitle(effname2);
 927   sprintf(effname2,"Uncorrected %s",yPhiEtaMixRaw->GetTitle());
 928   yPhiEtaMixRaw->SetTitle(effname2);
 929   */
 930 }
 931
 932 //Efficiency Fits-----------------------------------------------------------------------
 933 void EffFit(Float_t yAPt1,Float_t yAPt2,Int_t Cent,TList *effList, TH1F *yPhiEff, TH1F *yPhiMixEff, TH1F *yEtaNEff, TH1F *yEtaNMixEff, TH1F* yEtaAEff, TH1F *yEtaAMixEff, TH2F *yPhiEtaEff, TH2F *yPhiEtaMixEff, Int_t LSign=0, Int_t VariablePt=0){
 934
 935   const int xnpt=12;
 936   //float xptbins[(xnpt+1)]={2.5,3,4,6,8,10,15,20,30,40,50,75,100};
 937   float xptbins[(xnpt+1)]={2,2.5,3,4,5,6,8,10,15,20,30,40,50};
 938
 939   int maxPtTrig=10;//we don't want to go beyond our resolution/statistics
 940   //if(yAPt1>0.5)maxPtTrig=15;
 941   if(yAPt1>1)maxPtTrig=8;
 942   if(yAPt1>1.5)maxPtTrig=6;
 943
 944   int minPtTrig=yAPt1;
 945   //set up trigger array
 946   int size1=0;
 947   for(int i=0;i<xnpt;i++)if(xptbins[i]>(yAPt1+0.0001)&&xptbins[i]<(maxPtTrig-0.0001))size1++;
 948   const int NTPt=size1;
 949   float yPi=3.1415962;
 950   Float_t X[NTPt];
 951   Float_t T0[2][NTPt];
 952   Float_t T1[2][NTPt];
 953   Float_t T2[2][NTPt];
 954   Float_t T3[2][NTPt];
 955   Float_t T4[2][NTPt];
 956   Float_t M0[2][NTPt];
 957   Float_t NT0[2][NTPt];
 958   Float_t NT1[2][NTPt];
 959   Float_t NT2[2][NTPt];
 960   Float_t NT3[2][NTPt];
 961   Float_t NT4[2][NTPt];
 962   Float_t NM0[2][NTPt];
 963   Float_t AT0[2][NTPt];
 964   Float_t AM0[2][NTPt];
 965
 966   TH1F *hPhiEff=new TH1F("hPhiEff","",1,0,1);
 967   TH1F *hPhiMC=new TH1F("hPhiMC","",1,0,1);
 968   TH1F *hPhiMixEff=new TH1F("hPhiMixEff","",1,0,1);
 969   TH1F *hPhiMixMC=new TH1F("hPhiMixMC","",1,0,1);
 970
 971   TH1F *hEtaNEff=new TH1F("hEtaNEff","",1,0,1);
 972   TH1F *hEtaNMC=new TH1F("hEtaNMC","",1,0,1);
 973   TH1F *hEtaNMixEff=new TH1F("hEtaNMixEff","",1,0,1);
 974   TH1F *hEtaNMixMC=new TH1F("hEtaNMixMC","",1,0,1);
 975
 976   TH1F *hEtaAEff=new TH1F("hEtaAEff","",1,0,1);
 977   TH1F *hEtaAMC=new TH1F("hEtaAMC","",1,0,1);
 978   TH1F *hEtaAMixEff=new TH1F("hEtaAMixEff","",1,0,1);
 979   TH1F *hEtaAMixMC=new TH1F("hEtaAMixMC","",1,0,1);
 980
 981   TH2F *hPhiEtaEff=new TH2F("hPhiEtaEff","",1,0,1,1,0,1);
 982   TH2F *hPhiEtaMC=new TH2F("hPhiEtaMC","",1,0,1,1,0,1);
 983   TH2F *hPhiEtaMixEff=new TH2F("hPhiEtaMixEff","",1,0,1,1,0,1);
 984   TH2F *hPhiEtaMixMC=new TH2F("hPhiEtaMixMC","",1,0,1,1,0,1);
 985
 986   TF1 *fit1=new TF1("fit1","[0]+[1]/[2]*exp(-0.5*pow(x/[2],2))+[3]/[4]*exp(-0.5*pow((x-3.1415962)/[4],2))",-1.5,6);
 987   TF1 *fit2=new TF1("fit2","[0]",-2,6);
 988  TF1 *fit3=new TF1("fit3","[0]+[1]/[2]*exp(-0.5*pow(x/[2],2))",-2,6);
 989   fit1->SetParameters(1,0.1,0.1,0.1,0.1);
 990   fit2->SetParameter(0,1);
 991   fit3->SetParameters(1,0.1,0.1);
 992   fit1->SetParLimits(0,0,1000);
 993   fit1->SetParLimits(1,0,10);
 994   fit1->SetParLimits(2,0,10);
 995   fit1->SetParLimits(3,0,10);
 996   fit1->SetParLimits(4,0,10);
 997   fit3->SetParLimits(0,0,1000);
 998   fit3->SetParLimits(1,0,10);
 999   fit3->SetParLimits(2,0,10);
1000
1001   float MPt[4];
1002   float TPt1,TPt2,APt1,APt2;
1003   for(int i=0;i<NTPt;i++){
1004     TPt1=xptbins[i];
1005     TPt2=xptbins[i+1];
1006     APt1=yAPt1;
1007     APt2=yAPt2;
1008     if(VariablePt)APt2=maxPtTrig;
1009
1010     MakeProjections(TPt1,TPt2,APt1,APt2,Cent,effList,hPhiEff,hPhiMixEff,hEtaNEff,hEtaNMixEff,hEtaAEff,hEtaAMixEff,hPhiEtaEff,hPhiEtaMixEff,MPt,0,0,0,LSign);
1011     MakeProjections(TPt1,TPt2,APt1,APt2,Cent,effList,hPhiMC,hPhiMixMC,hEtaNMC,hEtaNMixMC,hEtaAMC,hEtaAMixMC,hPhiEtaMC,hPhiEtaMixMC,MPt,0,1,0,LSign);
1012     // EffCorr2(TPt1,TPt2,APt1,APt2,Cent,hPhiEff,hPhiMC,hPhiMixEff,hPhiMixMC,hEtaNEff,hEtaNMC,hEtaNMixEff,hEtaNMixMC,hEtaAEff,hEtaAMC,hEtaAMixEff,hEtaAMixMC,hPhiEtaEff,hPhiEtaMC,hPhiEtaMixEff,hPhiEtaMixMC,1);
1013     EffCorr2(TPt1,TPt2,APt1,APt2,Cent,hPhiEff,hPhiMC,hPhiMixEff,hPhiMixMC,hEtaNEff,hEtaNMC,hEtaNMixEff,hEtaNMixMC,hEtaAEff,hEtaAMC,hEtaAMixEff,hEtaAMixMC,hPhiEtaEff,hPhiEtaMC,hPhiEtaMixEff,hPhiEtaMixMC,1);
1014     X[i]=(TPt1+TPt2)*0.5;
1015
1016     hPhiEff->Fit("fit1");
1017     T0[0][i]=fit1->GetParameter(0);
1018     T0[1][i]=fit1->GetParError(0);
1019     T1[0][i]=fit1->GetParameter(1);
1020     T1[1][i]=fit1->GetParError(1);
1021     T2[0][i]=fit1->GetParameter(2);
1022     T2[1][i]=fit1->GetParError(2);
1023     T3[0][i]=fit1->GetParameter(3);
1024     T3[1][i]=fit1->GetParError(3);
1025     T4[0][i]=fit1->GetParameter(4);
1026     T4[1][i]=fit1->GetParError(4);
1027
1028     hPhiMixEff->Fit("fit2");
1029     M0[0][i]=fit2->GetParameter(0);
1030     M0[1][i]=fit2->GetParError(0);
1031
1032     hEtaNEff->Fit("fit3");
1033     NT0[0][i]=fit3->GetParameter(0);
1034     NT0[1][i]=fit3->GetParError(0);
1035     NT1[0][i]=fit3->GetParameter(1);
1036     NT1[1][i]=fit3->GetParError(1);
1037     NT2[0][i]=fit3->GetParameter(2);
1038     NT2[1][i]=fit3->GetParError(2);
1039
1040
1041     hEtaNMixEff->Fit("fit2");
1042     NM0[0][i]=fit2->GetParameter(0);
1043     NM0[1][i]=fit2->GetParError(0);
1044
1045     hEtaAEff->Fit("fit2");
1046     AT0[0][i]=fit2->GetParameter(0);
1047     AT0[1][i]=fit2->GetParError(0);
1048
1049     hEtaNMixEff->Fit("fit2");
1050     AM0[0][i]=fit2->GetParameter(0);
1051     AM0[1][i]=fit2->GetParError(0);
1052   }
1053
1054   gT0=new TGraphErrors(NTPt,X,T0[0],0,T0[1]);
1055   gT1=new TGraphErrors(NTPt,X,T1[0],0,T1[1]);
1056   gT2=new TGraphErrors(NTPt,X,T2[0],0,T2[1]);
1057   gT3=new TGraphErrors(NTPt,X,T3[0],0,T3[1]);
1058   gT4=new TGraphErrors(NTPt,X,T4[0],0,T4[1]);
1059
1060   gM0=new TGraphErrors(NTPt,X,M0[0],0,M0[1]);
1061
1062   gNT0=new TGraphErrors(NTPt,X,NT0[0],0,NT0[1]);
1063   gNT1=new TGraphErrors(NTPt,X,NT1[0],0,NT1[1]);
1064   gNT2=new TGraphErrors(NTPt,X,NT2[0],0,NT2[1]);
1065
1066   gNM0=new TGraphErrors(NTPt,X,NM0[0],0,NM0[1]);
1067
1068   gAT0=new TGraphErrors(NTPt,X,AT0[0],0,AT0[1]);
1069   gAM0=new TGraphErrors(NTPt,X,AM0[0],0,AM0[1]);
1070
1071   float pT0[2],pT1[2],pT2[2],pT3[2],pT4[2];
1072   float pM0[2];
1073   float pNT0[2],pNT1[2],pNT2[2],pNT3[2],pNT4[2];
1074   float pNM0[2];
1075   float pAT0[2];
1076   float pAM0[2];
1077
1078   gT0->Fit("fit2");
1079   pT0[0]=fit2->GetParameter(0);
1080   pT0[1]=fit2->GetParError(0);
1081   gT1->Fit("fit2");
1082   pT1[0]=fit2->GetParameter(0);
1083   pT1[1]=fit2->GetParError(0);
1084   gT2->Fit("fit2");
1085   pT2[0]=fit2->GetParameter(0);
1086   pT2[1]=fit2->GetParError(0);
1087   gT3->Fit("fit2");
1088   pT3[0]=fit2->GetParameter(0);
1089   pT3[1]=fit2->GetParError(0);
1090   gT4->Fit("fit2");
1091   pT4[0]=fit2->GetParameter(0);
1092   pT4[1]=fit2->GetParError(0);
1093
1094   gNT0->Fit("fit2");
1095   pNT0[0]=fit2->GetParameter(0);
1096   pNT0[1]=fit2->GetParError(0);
1097   gNT1->Fit("fit2");
1098   pNT1[0]=fit2->GetParameter(0);
1099   pNT1[1]=fit2->GetParError(0);
1100   gNT2->Fit("fit2");
1101   pNT2[0]=fit2->GetParameter(0);
1102   pNT2[1]=fit2->GetParError(0);
1103
1104   gM0->Fit("fit2");
1105   pM0[0]=fit2->GetParameter(0);
1106   pM0[1]=fit2->GetParError(0);
1107
1108   gNM0->Fit("fit2");
1109   pNM0[0]=fit2->GetParameter(0);
1110   pNM0[1]=fit2->GetParError(0);
1111
1112   gAM0->Fit("fit2");
1113   pAM0[0]=fit2->GetParameter(0);
1114   pAM0[1]=fit2->GetParError(0);
1115
1116   gAT0->Fit("fit2");
1117   pAT0[0]=fit2->GetParameter(0);
1118   pAT0[1]=fit2->GetParError(0);
1119
1120   float eff,gaus,err,cent;
1121   float gausA,errA;
1122   for(int x=1;x<=yPhiEff->GetNbinsX();x++){
1123     cent=yPhiEff->GetBinCenter(x);
1124     if(pT2[0]==0)cout << "Error pT2[0]" << endl;
1125     gaus=pT1[0]/pT2[0]*exp(-0.5*pow(cent/pT2[0],2));
1126     gausA=pT3[0]/pT4[0]*exp(-0.5*pow((cent-yPi)/pT4[0],2));
1127     eff=gaus+gausA+pT0[0];
1128     if(pT1[0]==0) cout << "Error pT1[0]" << endl;
1129     if(gaus==0)gaus+=1E-10;//divide by 0 protection
1130     if(gausA==0)gausA+=1E-10;
1131     err=pow(pow(pT0[1],2)+pow(pT1[1]*gaus/pT1[0],2)+pow(pT2[1]*gaus/pT2[0]*(pow(cent/pT2[0],2)-1),2)+pow(pT3[1]*gausA/pT3[0],2)+pow(pT4[1]*gausA/pT4[0]*(pow((cent-yPi)/pT4[0],2)-1),2),0.5);
1132     yPhiEff->SetBinContent(x,eff);
1133     yPhiEff->SetBinError(x,err);
1134
1135     yPhiMixEff->SetBinContent(x,pM0[0]);
1136     yPhiMixEff->SetBinError(x,pM0[1]);
1137   }
1138   for(int x=1;x<=yEtaNEff->GetNbinsX();x++){
1139     cent=yEtaNEff->GetBinCenter(x);
1140     if(pNT2[0]==0)cout << "Error pNT2[0]" << endl;
1141     gaus=pNT1[0]/pNT2[0]*exp(-0.5*pow(cent/pNT2[0],2));
1142     eff=gaus+pNT0[0];
1143     if(pNT1[0]==0) cout << "Error pNT1[0]" << endl;
1144     err=pow(pow(pNT0[1],2)+pow(pNT1[1]*gaus/pNT1[0],2)+pow(pNT2[1]*gaus/pNT2[0]*(pow(cent/pNT2[0],2)-1),2),0.5);
1145     yEtaNEff->SetBinContent(x,eff);
1146     yEtaNEff->SetBinError(x,err);
1147
1148     yEtaNMixEff->SetBinContent(x,pNM0[0]);
1149     yEtaNMixEff->SetBinError(x,pNM0[1]);
1150
1151     yEtaAEff->SetBinContent(x,pAT0[0]);
1152     yEtaAEff->SetBinError(x,pAT0[1]);
1153
1154     yEtaAMixEff->SetBinContent(x,pAM0[0]);
1155     yEtaAMixEff->SetBinError(x,pAM0[1]);
1156   }
1157
1158   float gaus2,eff2,err2,cent2;
1159   float eff3,err3;
1160   for(int x=1;x<=yPhiEtaEff->GetXaxis()->GetNbins();x++){
1161     cent=yPhiEtaEff->GetXaxis()->GetBinCenter(x);
1162     for(int y=1;y<=yPhiEtaEff->GetYaxis()->GetNbins();y++){
1163       if(fabs(cent)<1.5||fabs(cent+2*yPi)<1.5||fabs(cent-2*yPi)<1.5){
1164       gaus=pT1[0]/pT2[0]*exp(-0.5*pow(cent/pT2[0],2));
1165       eff=gaus;
1166       if(eff==0)eff+=1E-10;//divide by 0 protection
1167       err=pow(pow(pT1[1]*gaus/pT1[0],2)+pow(pT2[1]*gaus/pT2[0]*(pow(cent/pT2[0],2)-1),2),0.5);
1168       cent2=yPhiEtaEff->GetYaxis()->GetBinCenter(y);
1169       gaus2=pNT1[0]/pNT2[0]*exp(-0.5*pow(cent2/pNT2[0],2));
1170       eff2=gaus2;
1171       if(eff2==0)eff2+=1E-10;//divide by 0 protection
1172       err2=pow(pow(pNT1[1]*gaus2/pNT1[0],2)+pow(pNT2[1]*gaus2/pNT2[0]*(pow(cent2/pNT2[0],2)-1),2),0.5);
1173       // eff3=pow(eff*eff2,0.5);
1174       eff3=pow(gaus*gaus2,0.5)+pow((pT0[0])*pNT0[0],0.5);
1175       // eff=pow(pT0[0]*pNT0[0],0.5);
1176       err3=pow(gaus*gaus2*(pow(err/eff,2)+pow(err2/eff2,2))+pT0[0]*pNT0[0]*(pow(pT0[1]/pT0[0],2)+pow(pNT0[1]/pNT0[0],2)),0.5);
1177
1178       }
1179       else{
1180         gausA=pT3[0]/pT4[0]*exp(-0.5*pow((cent-yPi)/pT4[0],2));
1181         if(gausA==0)gausA+=1E-10;
1182         eff3=pow((pT0[0]+gausA)*pAT0[0],0.5);
1183         //err3=pow(pow(pT0[1]/pT0[0],2)+pow(pAT0[1]/pAT0[0],2),0.5);
1184         err3=eff3*pow((pow(pT0[1],2)+pow(pT3[1]*gausA/pT3[0],2)+pow(pT4[1]*gausA/pT4[0]*(pow((cent-yPi)/pT4[0],2)-1),2))/pow(pT0[0]+gausA,2)+pow(pAT0[1]/pAT0[1],2),0.5);
1185       }
1186        yPhiEtaEff->SetBinContent(x,y,eff3);
1187       yPhiEtaEff->SetBinError(x,y,err3);
1188
1189       if(fabs(cent)<1.5){
1190         eff3=pow(pM0[0]*pNM0[0],0.5);
1191         err3=pow(pow(pM0[1]/pM0[0],2)+pow(pNM0[1]/pNM0[0],2),0.5);
1192       }
1193       else{
1194         eff3=pow(pM0[0]*pAM0[0],0.5);
1195         err3=pow(pow(pM0[1]/pM0[0],2)+pow(pAM0[1]/pAM0[0],2),0.5);
1196       }
1197       yPhiEtaMixEff->SetBinContent(x,y,eff3);
1198       yPhiEtaMixEff->SetBinError(x,y,err3);
1199     }
1200   }
1201   }
1202
1203  //Geometry Correction for the eta-phi plots
1204  //\-------------------------------------------------------------------
1205 void GeoCorr(TH1F *yEtaNCorr, TH1F *yEtaNMixCorr, TH1F *yEtaNMixMC, TH1F *yEtaACorr, TH1F *yEtaAMixCorr, TH1F *yEtaAMixMC, TH2F *yPhiEtaCorr, TH2F *yPhiEtaMixCorr, TH2F *yPhiEtaMixMC){
1206
1207   int xbins=yEtaNCorr->GetNbinsX();
1208
1209   float corr;
1210   float center;
1211   float scale=0,scaleN,scaleA;
1212   float sumy;
1213   float escaleN,escaleA, ecorr;
1214   float con,err;
1215   scaleN=(yEtaNMix->GetBinContent(xbins/2)+yEtaNMix->GetBinContent(xbins/2+1))/2.;
1216   escaleN=pow(yEtaNMix->GetBinError(xbins/2),2)+pow(yEtaNMix->GetBinError(xbins/2+1),2)/4.;
1217   scaleA=(yEtaAMix->GetBinContent(xbins/2)+yEtaAMix->GetBinContent(xbins/2+1))/2.;
1218   escaleA=pow(yEtaAMix->GetBinError(xbins/2),2)+pow(yEtaAMix->GetBinError(xbins/2+1),2)/4.;
1219   if(scaleN==0)scaleN=1E-5;
1220   if(scaleA==0)scaleA=1E-5;
1221   for(int yx=1;yx<=xbins;yx++){
1222     con=yEtaNMixMC->GetBinContent(yx);
1223     if(con==0){con=1E-5;corr=1E-5;}//divide by 0 protection
1224     else corr=scaleN/con;
1225     ecorr=corr*corr*(escaleN/scaleN/scaleN+pow(yEtaNMixMC->GetBinError(yx),2)/con/con);
1226     con=yEtaNCorr->GetBinContent(yx);
1227     err=pow(yEtaNCorr->GetBinError(yx),2);
1228     yEtaNCorr->SetBinContent(yx,corr*con);
1229     if(con==0)con=1E-5;//Divide by 0 protection
1230     if(corr==0)corr=1E-5;
1231     yEtaNCorr->SetBinError(yx,corr*con*pow(ecorr/corr/corr+err/con/con,0.5));
1232     con=yEtaNMixCorr->GetBinContent(yx);
1233     err=pow(yEtaNMixCorr->GetBinError(yx),2);
1234     yEtaNMixCorr->SetBinContent(yx,corr*con);
1235     if(con==0)con=1E-5;
1236     if(corr==0)corr=1E-5;
1237     yEtaNMixCorr->SetBinError(yx,corr*con*pow(ecorr/corr/corr+err/con/con,0.5));
1238
1239     con=yEtaAMixMC->GetBinContent(yx);
1240     if(con==0){con=1E-5;corr=1E-5;}
1241     else corr=scaleA/con;
1242     ecorr=corr*corr*(escaleA/scaleA/scaleA+pow(yEtaAMixMC->GetBinError(yx),2)/con/con);
1243     con=yEtaACorr->GetBinContent(yx);
1244     err=pow(yEtaACorr->GetBinError(yx),2);
1245     yEtaACorr->SetBinContent(yx,corr*con);
1246     if(con==0)con=1E-5;
1247     if(corr==0)corr=1E-5;
1248     yEtaACorr->SetBinError(yx,corr*con*pow(ecorr/corr/corr+err/con/con,0.5));
1249     con=yEtaAMixCorr->GetBinContent(yx);
1250     err=pow(yEtaAMixCorr->GetBinError(yx),2);
1251     yEtaAMixCorr->SetBinContent(yx,corr*con);
1252     if(con==0)con=1E-5;
1253     if(corr==0)corr=1E-5;
1254     yEtaAMixCorr->SetBinError(yx,corr*con*pow(ecorr/corr/corr+err/con/con,0.5));
1255   }//end eta loop
1256
1257   xbins=yPhiEtaCorr->GetNbinsX();
1258   int ybins=yPhiEtaCorr->GetNbinsY();
1259   scale=0;
1260   //uses the mc for the geomery correction
1261   for(int yx=1;yx<=xbins;yx++){
1262     scale+=yPhiEtaMixMC->GetBinContent(yx,ybins/2);
1263     scale+=yPhiEtaMixMC->GetBinContent(yx,ybins/2+1);
1264   }
1265   scale/=2.;
1266   for(int yy=1;yy<=ybins;yy++){
1267     sumy=0;
1268     for(int yx=1;yx<=xbins;yx++){
1269       sumy+=yPhiEtaMixMC->GetBinContent(yx,yy);
1270     }
1271     if(sumy==0)sumy=0.00001;
1272     corr=scale/sumy;
1273     for(int yx=1;yx<=xbins;yx++){
1274       yPhiEtaCorr->SetBinContent(yx,yy,corr*yPhiEtaCorr->GetBinContent(yx,yy));
1275       yPhiEtaCorr->SetBinError(yx,yy,corr*yPhiEtaCorr->GetBinError(yx,yy));
1276       yPhiEtaMixCorr->SetBinContent(yx,yy,corr*yPhiEtaMixCorr->GetBinContent(yx,yy));
1277       yPhiEtaMixCorr->SetBinError(yx,yy,corr*yPhiEtaMixCorr->GetBinError(yx,yy));
1278     }
1279   }
1280 }
1281 //-----------------------------------------------------------------------------------------------
1282 void GeoCorr2(TH1F *yEtaNCorr, TH1F *yEtaNMixCorr, TH1F *yEtaACorr, TH1F *yEtaAMixCorr, TH2F *yPhiEtaCorr, TH2F *yPhiEtaMixCorr){
1283   int nbins=yEtaNCorr->GetNbinsX();
1284   int nbinsPhi=yPhiEtaCorr->GetNbinsX();
1285   int nbinsEta=yPhiEtaCorr->GetNbinsY();
1286   float EtaCut=yEtaNCorr->GetBinCenter(nbins)+yEtaNCorr->GetBinWidth(nbins)/2.;
1287   float corr;
1288   float cent;
1289   for(int i=1;i<=nbins;i++){
1290     cent=yEtaNCorr->GetBinCenter(i);
1291     corr=1-fabs(cent)/EtaCut;
1292     yEtaNCorr->SetBinContent(i,yEtaNCorr->GetBinContent(i)/corr);
1293     yEtaNCorr->SetBinError(i,yEtaNCorr->GetBinError(i)/corr);
1294     yEtaNMixCorr->SetBinContent(i,yEtaNMixCorr->GetBinContent(i)/corr);
1295     yEtaNMixCorr->SetBinError(i,yEtaNMixCorr->GetBinError(i)/corr);
1296     yEtaACorr->SetBinContent(i,yEtaACorr->GetBinContent(i)/corr);
1297     yEtaACorr->SetBinError(i,yEtaACorr->GetBinError(i)/corr);
1298     yEtaAMixCorr->SetBinContent(i,yEtaAMixCorr->GetBinContent(i)/corr);
1299     yEtaAMixCorr->SetBinError(i,yEtaAMixCorr->GetBinError(i)/corr);
1300   }
1301   for(int i=1;i<=nbinsEta;i++){
1302     cent=yPhiEtaCorr->GetYaxis()->GetBinCenter(i);
1303     corr=1-fabs(cent)/EtaCut;
1304     for(int j=1;j<=nbinsPhi;j++){
1305       yPhiEtaCorr->SetBinContent(j,i,yPhiEtaCorr->GetBinContent(j,i)/corr);
1306       yPhiEtaCorr->SetBinError(j,i,yPhiEtaCorr->GetBinError(j,i)/corr);
1307       yPhiEtaMixCorr->SetBinContent(j,i,yPhiEtaMixCorr->GetBinContent(j,i)/corr);
1308       yPhiEtaMixCorr->SetBinError(j,i,yPhiEtaMixCorr->GetBinError(j,i)/corr);
1309     }
1310   }
1311 }
1312 //---------------------------------------------------------------------------------
1313 void Load3Particle(Float_t &xTPt1, Float_t &xTPt2, Float_t &xAPt1, Float_t &xAPt2,Int_t xCent, TList *xList, TH2F *xPhiPhiRaw,TH2F *xPhiPhiSS,TH2F *xPhiPhiMix, TH2F *xEtaEtaRaw, TH2F *xEtaEtaSS, TH2F *xEtaEtaMix,Int_t xMC=0, Int_t xSign=0){
1314   const int xnAPt=7;
1315   float xAPt31[xnAPt]={0.5,1.0,1.5,2.0,3,4,1};
1316   float xAPt32[xnAPt]={1.0,1.5,2.0,3.0,4,5,2};
1317   const int xnTPt=12;
1318   float xTPt[(xnTPt+1)]={2,2.5,3,4,5,6,8,10,15,20,30,40,50};//any element in this array should also been in the one 2 lines down
1319  float xPi=3.1415962;
1320   int APtBin=-1;
1321   int TPtBin1=-1;
1322   int TPtBin2=-1;
1323   for(int i=0;i<xnAPt;i++){
1324     if(fabs(xAPt1-xAPt31[i])<0.1&&fabs(xAPt2-xAPt32[i])<0.1){
1325       xAPt1=xAPt31[i];
1326       xAPt2=xAPt32[i];
1327       APtBin=i;
1328     }
1329   }
1330   if(APtBin==-1){
1331     cout << "Invalid Associated Pt: " << xAPt1 << " " << xAPt2 << endl;
1332     cout << "Valid values are:" << endl;
1333     for(int i=0;i<xnAPt;i++){
1334       cout << "APt1: " << xAPt31[i] << "  APt2: " << xAPt32[i] << endl;
1335     }
1336     break;
1337   }
1338   for(int i=0;i<xnTPt;i++){
1339     if(fabs(xTPt1-xTPt[i])<0.1){
1340      xTPt1=xTPt[i];
1341      TPtBin1=i;
1342     }
1343     if(fabs(xTPt2-xTPt[i+1])<0.1){
1344       xTPt2=xTPt[i+1];
1345       TPtBin2=i;
1346     }
1347   }
1348   if(TPtBin1==-1||TPtBin2==-1){
1349     cout << "Invalid Trigger Pt: " << xTPt1 << " " << xTPt2 << end;
1350      cout << "Valid values are:" << endl;
1351     for(int i=0;i<=xnTPt;i++){
1352       cout <<  xTPt[i] << "  ";
1353     }
1354     cout << endl;
1355     break;
1356   }
1357   char xname[100];
1358   char xtit[200];
1359  char *cmc1[2]={"","_MC"};
1360   sprintf(xname,"fHistNTrigger_C%d%s",xCent,cmc1[xMC]);
1361   TH1F *xNTrig=(TH1F*)xList->FindObject(xname);
1362   sprintf(xname,"fHistNMix_C%d%s",xCent,cmc1[xMC]);
1363   TH1F *xNMix=(TH1F*)xList->FindObject(xname);
1364   int ntriggers=xNTrig->GetBinContent(1);
1365   int nmix=xNMix->GetBinContent(1);
1366   // c100=new TCanvas("c100");
1367   // xNTrig->Draw();
1368   TH2F *xtemphistP;
1369   TH2F *xtemphistPSS;
1370   TH2F *xtemphistPMix;
1371   TH2F *xtemphistE;
1372   TH2F *xtemphistESS;
1373   TH2F *xtemphistEMix;
1374   //char *cmc1[2]={"","_MC"};
1375   char *sign1[3]={"","_LS","_ULS"};
1376   char *sign2[3]={""," LikeSign"," UnlikeSign"};
1377   int nbins;
1378   float min,max;
1379   float conP, conPSS, conPMix, conE, conESS, conEMix;
1380   float errP, errPSS, errPMix, errE, errESS, errEMix;
1381   float xntrig=0;
1382   float xnmix=0;
1383
1384
1385   for(int xi=TPtBin1;xi<=TPtBin2;xi++){//Do for all trigger pt ranges we are adding
1386     xntrig+=xNTrig->GetBinContent(xi+1);
1387     //cout << xi << " " << xNTrig->GetBinContent(xi+1) << endl;
1388     xnmix+=xNMix->GetBinContent(xi+1);
1389     sprintf(xname,"fHistDeltaPhiPhi_P%dp%d_C%d%s%s",xi,APtBin,xCent,cmc1[xMC],sign1[xSign]);
1390     xtemphistP=(TH2F*)xList->FindObject(xname);
1391
1392     sprintf(xname,"fHistDeltaPhiPhiSS_P%dp%d_C%d%s%s",xi,APtBin,xCent,cmc1[xMC],sign1[xSign]);
1393     xtemphistPSS=(TH2F*)xList->FindObject(xname);
1394
1395     sprintf(xname,"fHistDeltaPhiPhiMix_P%dp%d_C%d%s%s",xi,APtBin,xCent,cmc1[xMC],sign1[xSign]);
1396     xtemphistPMix=(TH2F*)xList->FindObject(xname);
1397
1398     sprintf(xname,"fHistDeltaEtaEta_P%dp%d_C%d%s%s",xi,APtBin,xCent,cmc1[xMC],sign1[xSign]);
1399     xtemphistE=(TH2F*)xList->FindObject(xname);
1400
1401     sprintf(xname,"fHistDeltaEtaEtaSS_P%dp%d_C%d%s%s",xi,APtBin,xCent,cmc1[xMC],sign1[xSign]);
1402     xtemphistESS=(TH2F*)xList->FindObject(xname);
1403
1404     sprintf(xname,"fHistDeltaEtaEtaMix_P%dp%d_C%d%s%s",xi,APtBin,xCent,cmc1[xMC],sign1[xSign]);
1405     xtemphistEMix=(TH2F*)xList->FindObject(xname);
1406
1407     if(xi==TPtBin1){
1408       nbins=xtemphistP->GetNbinsX();
1409       min=xtemphistP->GetBinCenter(1)-xtemphistP->GetBinWidth(1)/2.;
1410       max=xtemphistP->GetBinCenter(nbins)+xtemphistP->GetBinWidth(nbins)/2.;
1411
1412       sprintf(xtit,"#Delta#phi-#Delta#phi Distribution for  %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f%s",xTPt1,xTPt2,xAPt1,xAPt2,sign2[xSign]);
1413       *xPhiPhiRaw=new TH2F("dPhiPhiRaw",xtit,nbins,min,max,nbins,min,max);
1414       xPhiPhiRaw->Sumw2();
1415       xPhiPhiRaw->GetXaxis()->SetTitle("#Delta#phi");
1416       xPhiPhiRaw->GetYaxis()->SetTitle("#Delta#phi");
1417       xPhiPhiRaw->GetZaxis()->SetTitle("#frac{1}{N_{Trig}}#frac{d^{2}N_{pairs}}{d#Delta#phid#Delta#phi}");
1418       SetTitles2D(xPhiPhiRaw);
1419
1420       sprintf(xtit,"Soft-Soft #Delta#phi-#Delta#phi Distribution for  %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f%s",xTPt1,xTPt2,xAPt1,xAPt2,sign2[xSign]);
1421       *xPhiPhiSS=new TH2F("dPhiPhiSS",xtit,nbins,min,max,nbins,min,max);
1422       xPhiPhiSS->Sumw2();
1423       xPhiPhiSS->GetXaxis()->SetTitle("#Delta#phi");
1424       xPhiPhiSS->GetYaxis()->SetTitle("#Delta#phi");
1425       xPhiPhiSS->GetZaxis()->SetTitle("#frac{1}{N_{Trig}}#frac{d^{2}N_{pairs}}{d#Delta#phid#Delta#phi}");
1426       SetTitles2D(xPhiPhiSS);
1427
1428        sprintf(xtit,"Mixed #Delta#phi-#Delta#phi Distribution for  %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f%s",xTPt1,xTPt2,xAPt1,xAPt2,sign2[xSign]);
1429       *xPhiPhiMix=new TH2F("dPhiPhiMix",xtit,nbins,min,max,nbins,min,max);
1430       xPhiPhiMix->Sumw2();
1431       xPhiPhiMix->GetXaxis()->SetTitle("#Delta#phi");
1432       xPhiPhiMix->GetYaxis()->SetTitle("#Delta#phi");
1433       xPhiPhiMix->GetZaxis()->SetTitle("#frac{1}{N_{Trig}}#frac{d^{2}N_{pairs}}{d#Delta#phid#Delta#phi}");
1434       SetTitles2D(xPhiPhiMix);
1435       nbins=xtemphistE->GetNbinsX();
1436       min=xtemphistE->GetBinCenter(1)-xtemphistE->GetBinWidth(1)/2.;
1437       max=xtemphistE->GetBinCenter(nbins)+xtemphistE->GetBinWidth(nbins)/2.;
1438       sprintf(xtit,"#Delta#eta-#Delta#eta Distribution for  %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f%s",xTPt1,xTPt2,xAPt1,xAPt2,sign2[xSign]);
1439       *xEtaEtaRaw=new TH2F("dEtaEtaRaw",xtit,nbins,min,max,nbins,min,max);
1440       xEtaEtaRaw->Sumw2();
1441       xEtaEtaRaw->GetXaxis()->SetTitle("#Delta#eta");
1442       xEtaEtaRaw->GetYaxis()->SetTitle("#Delta#eta");
1443       xEtaEtaRaw->GetZaxis()->SetTitle("#frac{1}{N_{Trig}}#frac{d^{2}N_{pairs}}{d#Delta#etad#Delta#eta}");
1444       SetTitles2D(xEtaEtaRaw);
1445
1446        sprintf(xtit,"Soft-Soft #Delta#eta-#Delta#eta Distribution for  %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f%s",xTPt1,xTPt2,xAPt1,xAPt2,sign2[xSign]);
1447        *xEtaEtaSS=new TH2F("dEtaEtaSS",xtit,nbins,min,max,nbins,min,max);
1448       xEtaEtaSS->Sumw2();
1449       xEtaEtaSS->GetXaxis()->SetTitle("#Delta#eta");
1450       xEtaEtaSS->GetYaxis()->SetTitle("#Delta#eta");
1451       xEtaEtaSS->GetZaxis()->SetTitle("#frac{1}{N_{Trig}}#frac{d^{2}N_{pairs}}{d#Delta#etad#Delta#eta}");
1452       SetTitles2D(xEtaEtaSS);
1453
1454        sprintf(xtit,"Mixed #Delta#eta-#Delta#eta Distribution for  %3.1f<p_{T}^{Trig}<%3.1f %2.2f<p_{T}^{Assoc}<%2.2f%s",xTPt1,xTPt2,xAPt1,xAPt2,sign2[xSign]);
1455       *xEtaEtaMix=new TH2F("dEtaEtaMix",xtit,nbins,min,max,nbins,min,max);
1456       xEtaEtaMix->Sumw2();
1457       xEtaEtaMix->GetXaxis()->SetTitle("#Delta#eta");
1458       xEtaEtaMix->GetYaxis()->SetTitle("#Delta#eta");
1459       xEtaEtaMix->GetZaxis()->SetTitle("#frac{1}{N_{Trig}}#frac{d^{2}N_{pairs}}{d#Delta#etad#Delta#eta}");
1460       SetTitles2D(xEtaEtaMix);
1461
1462     }
1463     for(int x=1;x<=xtemphistP->GetNbinsX();x++){
1464       for(int y=1;y<=xtemphistP->GetNbinsY();y++){
1465         if(xi==TPtBin1){
1466           conP=0;
1467           errP=0;
1468           conPSS=0;
1469           errPSS=0;
1470           conPMix=0;
1471           errPMix=0;
1472         }
1473         else{
1474           conP=xPhiPhiRaw->GetBinContent(x,y);
1475           errP=xPhiPhiRaw->GetBinError(x,y);
1476           conPSS=xPhiPhiSS->GetBinContent(x,y);
1477           errPSS=xPhiPhiSS->GetBinError(x,y);
1478           conPMix=xPhiPhiMix->GetBinContent(x,y);
1479           errPMix=xPhiPhiMix->GetBinError(x,y);
1480         }
1481         xPhiPhiRaw->SetBinContent(x,y,conP+xtemphistP->GetBinContent(x,y));
1482         xPhiPhiRaw->SetBinError(x,y,pow(errP*errP+pow(xtemphistP->GetBinError(x,y),2),0.5));
1483         xPhiPhiSS->SetBinContent(x,y,conPSS+xtemphistPSS->GetBinContent(x,y));
1484         xPhiPhiSS->SetBinError(x,y,pow(errPSS*errPSS+pow(xtemphistPSS->GetBinError(x,y),2),0.5));
1485         xPhiPhiMix->SetBinContent(x,y,conPMix+xtemphistPMix->GetBinContent(x,y));
1486         xPhiPhiMix->SetBinError(x,y,pow(errPMix*errPMix+pow(xtemphistPMix->GetBinError(x,y),2),0.5));
1487       }
1488     }
1489
1490     for(int x=1;x<=xtemphistE->GetNbinsX();x++){
1491       for(int y=1;y<=xtemphistE->GetNbinsY();y++){
1492         if(xi==TPtBin1){
1493           conE=0;
1494           errE=0;
1495           conESS=0;
1496           errESS=0;
1497           conEMix=0;
1498           errEMix=0;
1499         }
1500         else{
1501           conE=xEtaEtaRaw->GetBinContent(x,y);
1502           errE=xEtaEtaRaw->GetBinError(x,y);
1503           conESS=xEtaEtaSS->GetBinContent(x,y);
1504           errESS=xEtaEtaSS->GetBinError(x,y);
1505           conEMix=xEtaEtaMix->GetBinContent(x,y);
1506           errEMix=xEtaEtaMix->GetBinError(x,y);
1507         }
1508         xEtaEtaRaw->SetBinContent(x,y,conE+xtemphistE->GetBinContent(x,y));
1509         xEtaEtaRaw->SetBinError(x,y,pow(errE*errE+pow(xtemphistE->GetBinError(x,y),2),0.5));
1510         xEtaEtaSS->SetBinContent(x,y,conESS+xtemphistESS->GetBinContent(x,y));
1511         xEtaEtaSS->SetBinError(x,y,pow(errESS*errESS+pow(xtemphistESS->GetBinError(x,y),2),0.5));
1512         xEtaEtaMix->SetBinContent(x,y,conEMix+xtemphistEMix->GetBinContent(x,y));
1513         xEtaEtaMix->SetBinError(x,y,pow(errEMix*errEMix+pow(xtemphistEMix->GetBinError(x,y),2),0.5));
1514
1515       }
1516     }
1517   }
1518
1519 //Normalize the histograms
1520   //cout << xPi << " " << max << " " << xntrig << " " << xnmix << endl;
1521   float phibins=xPhiPhiRaw->GetNbinsX();
1522   float etabins=xEtaEtaRaw->GetNbinsX();
1523   float scalephi=pow(phibins/(2*xPi),2);
1524   float scaleeta=pow(etabins/(max-min),2);//max still set for eta from making histograms
1525   xPhiPhiRaw->Scale(scalephi/xntrig);
1526   xPhiPhiSS->Scale(scalephi/xnmix);
1527   xPhiPhiMix->Scale(scalephi/(2.*xnmix));
1528   xEtaEtaRaw->Scale(scaleeta/xntrig);
1529   xEtaEtaSS->Scale(scaleeta/xnmix);
1530   xEtaEtaMix->Scale(scaleeta/(2.*xnmix));
1531 }
1532
1533 //------------------------------------------------------------------------------------
1534 void EffCorr3Part(Float_t TPt1, Float_t TPt2, Float_t APt1, Float_t APt2, Int_t Cent, TH2F *yPhiPhiEff, TH2F *yPhiPhiMC, TH2F *yPhiPhiSSEff, TH2F *yPhiPhiSSMC, TH2F *yPhiPhiMixEff, TH2F *yPhiPhiMixMC, TH2F *yEtaEtaEff, TH2F *yEtaEtaMC, TH2F *yEtaEtaSSEff, TH2F *yEtaEtaSSMC, TH2F *yEtaEtaMixEff, TH2F *yEtaEtaMixMC, Int_t yMethod){
1535   //Right now only putting in for case where efficiency is taken care of before hand.  I will add others if used
1536   if(yMethod==0){
1537     int nBinPhi=yPhiPhiEff->GetNbinsX();
1538     int nBinEta=yEtaEtaEff->GetNbinsX();
1539     for(int x=1;x<=nBinPhi;x++){
1540       for(int y=1;y<=nBinPhi;y++){
1541         yPhiPhiEff->SetBinContent(x,y,1);
1542         yPhiPhiEff->SetBinError(x,y,0);
1543         yPhiPhiSSEff->SetBinContent(x,y,1);
1544         yPhiPhiSSEff->SetBinError(x,y,0);
1545         yPhiPhiMixEff->SetBinContent(x,y,1);
1546         yPhiPhiMixEff->SetBinError(x,y,0);
1547       }
1548     }
1549   for(int x=1;x<=nBinEta;x++){
1550       for(int y=1;y<=nBinEta;y++){
1551         yEtaEtaEff->SetBinContent(x,y,1);
1552         yEtaEtaEff->SetBinError(x,y,0);
1553         yEtaEtaSSEff->SetBinContent(x,y,1);
1554         yEtaEtaSSEff->SetBinError(x,y,0);
1555         yEtaEtaMixEff->SetBinContent(x,y,1);
1556         yEtaEtaMixEff->SetBinError(x,y,0);
1557       }
1558     }
1559   }
1560   else{
1561     Printf("Method needs to be inserted into EffCorr3Part");
1562     break;
1563   }
1564
1565
1566 }
1567
1568 //-------------------------------------------------------------------------------------------------------
1569 void GeoCorr3Part2(TH2F *yEtaEtaCorr, TH2F *yEtaEtaSSCorr, TH2F *yEtaEtaMixCorr, TH2F *yEtaEtaHSCorr){//using equation to correct
1570   int bin=yEtaEtaCorr->GetNbinsX();
1571   float centx, centy;
1572   float corr,corrx,corry,corrz;
1573   float EtaCut=yEtaEtaCorr->GetBinCenter(bin)+yEtaEtaCorr->GetBinWidth(bin)/2.;
1574   cout << "EtaCut " << EtaCut << endl;
1575   for(int x=1;x<=bin;x++){
1576     centx=yEtaEtaCorr->GetXaxis()->GetBinCenter(x);
1577     corrx=1-fabs(centx)/EtaCut;
1578     for(int y=1;y<=bin;y++){
1579       centy=yEtaEtaCorr->GetYaxis()->GetBinCenter(y);
1580       corry=1-fabs(centy)/EtaCut;
1581       corrz=1-fabs(centx-centy)/EtaCut;
1582       if(fabs(centx-centy)>EtaCut)corrz=0;
1583       corr=corrz*pow(corrx*corry,1/2.);
1584       if(corr==0)corr=1;
1585       corr=1;
1586       yEtaEtaCorr->SetBinContent(x,y,yEtaEtaCorr->GetBinContent(x,y)/corr);
1587       yEtaEtaCorr->SetBinError(x,y,yEtaEtaCorr->GetBinError(x,y)/corr);
1588       yEtaEtaMixCorr->SetBinContent(x,y,yEtaEtaMixCorr->GetBinContent(x,y)/corr);
1589       yEtaEtaMixCorr->SetBinError(x,y,yEtaEtaMixCorr->GetBinError(x,y)/corr);
1590       yEtaEtaSSCorr->SetBinContent(x,y,yEtaEtaSSCorr->GetBinContent(x,y)/corr);
1591       yEtaEtaSSCorr->SetBinError(x,y,yEtaEtaSSCorr->GetBinError(x,y)/corr);
1592       yEtaEtaHSCorr->SetBinContent(x,y,yEtaEtaHSCorr->GetBinContent(x,y)/corr);
1593       yEtaEtaHSCorr->SetBinError(x,y,yEtaEtaHSCorr->GetBinError(x,y)/corr);
1594     }
1595   }
1596 }
1597 //--------------------------------------------------------------------------------------------------------------------
1598 void ZYAM3(TH1F *yPhi, TH1F *yPhiMix, TH2F *yPhiPhi, TH2F *yPhiPhiSS, TH2F *yPhiPhiMix, Float_t a[3], Float_t b, Int_t nMin=10, Int_t nIter=10){//xxx scale the mixed event to the SS before it comes in
1599   const Int_t NMin=100;
1600
1601   // float sumSig2, sumMix2, sumSig3, sumSS3, sumMix3;
1602   //float esumSig2, esumMix2, esumSig3, esumSS3, esumMix3;
1603   //float sumMin, esumMin;
1604   float itersize=0.1;
1605   float oldmin=1000;
1606   a[0]=1;
1607   a[1]=0;
1608   float MinArray[NMin][2];
1609   float Phi[NMin][2];
1610   float SS[NMin][2];
1611   float Mix[NMin][2];
1612   float PhiX[NMin][2];
1613   float PhiY[NMin][2];
1614    float MixX[NMin][2];
1615   float MixY[NMin][2];
1616   int nBins2=yPhi->GetNbinsX();
1617   int nBins3=yPhiPhi->GetNbinsX();
1618   int ReBin=nBins2/nBins3;
1619     if(fabs(1.0*nBins2/nBins3-ReBin)>0.01)Printf("Please use a divisor of the 2-particle binning for the 3-particle");
1620
1621   float maxBinVal;
1622   int maxBin;
1623   float bin, err;
1624   float phi, ss, mix, phiX, mixX, phiY, mixY;
1625   float ephi, ess, emix, ephiX, emixX, ephiY, emixY;
1626   //float b1=yPhiPhi->GetSum()/pow(yPhi->GetSum(),2);
1627   //float b2=yPhiPhiSS->GetSum()/pow(yPhiMix->GetSum(),2);
1628   // b=b1/b2;
1629   for(int i=0;i<nIter;i++){
1630     for(int q=0;q<nMin;q++){
1631       MinArray[q][0]=10000;
1632     }
1633     for(int x=1;x<=yPhiPhi->GetNbinsX();x++){
1634       phiX=0, ephiX=0, mixX=0, emixX=0;
1635       for(int j=(ReBin*(x-1)+1);j<ReBin*x;j++){
1636         phiX+=yPhi->GetBinContent(j);
1637         ephiX+=pow(yPhi->GetBinError(j),2);
1638         mixX+=yPhiMix->GetBinContent(j);
1639         emixX+=pow(yPhiMix->GetBinError(j),2);
1640       }
1641       for(int y=1;y<=yPhiPhi->GetNbinsY();y++){
1642         phiY=0, ephiY=0, mixY=0, emixY=0;
1643         for(int j=(ReBin*(y-1)+1);j<ReBin*y;j++){
1644         phiY+=yPhi->GetBinContent(j);
1645         ephiY+=pow(yPhi->GetBinError(j),2);
1646         mixY+=yPhiMix->GetBinContent(j);
1647         emixY+=pow(yPhiMix->GetBinError(j),2);
1648       }
1649         maxBinVal=-10000;
1650         for(int j=0;j<nMin;j++){//Find highest bin
1651           if(MinArray[j][0]>maxBinVal){
1652             maxBinVal=MinArray[j][0];
1653             maxBin=j;
1654           }
1655         }
1656         phi=yPhiPhi->GetBinContent(x,y);
1657         ephi=yPhiPhi->GetBinError(x,y);
1658         ss=yPhiPhiSS->GetBinContent(x,y);
1659         ess=yPhiPhiSS->GetBinError(x,y);
1660         mix=yPhiPhiMix->GetBinContent(x,y);
1661         emix=yPhiPhiMix->GetBinError(x,y);
1662         //xxx
1663         bin=phi-a[0]*mixY*(phiX-a[0]*mixX)-a[0]*mixX*(phiY-a[0]*mixY)-a[0]*a[0]*b*ss;
1664         //cout << bin << " " << maxBinVal << " " << maxBin << endl;
1665         if(bin<maxBinVal){//replce highest bin with current if lower
1666           MinArray[maxBin][0]=bin;
1667           //divide by 0 protections
1668           if(a[0]==0)a[0]=1E-10;
1669           if(a[1]==0)a[1]=1E-10;
1670           if(ephi==0)ephi=1E-10;
1671           if(ess==0)ess=1E-10;
1672           if(emix==0)emix=1E-10;
1673           if(phi==0)phi=1E-10;
1674           if(ss==0)ss=1E-10;
1675           if(mix==0)mix=1E-10;
1676           MinArray[maxBin][1]=ephi*ephi+a[0]*mixY*phiX*(pow(emixY/mixY,2)+pow(ephiX/phiX,2))+2*a[0]*a[0]*mixY*mixX*(pow(emixY/mixY,2)+pow(emixX/mixX,2))+a[0]*mixX*phiY*(pow(emixX/mixX,2)+pow(ephiY/phiY,2))+pow(a[0]*a[0]*b*ess,2);
1677           Phi[maxBin][0]=phi;
1678           SS[maxBin][0]=ss;
1679           Mix[maxBin][0]=mix;
1680           PhiX[maxBin][0]=phiX;
1681           PhiY[maxBin][0]=phiY;
1682           MixX[maxBin][0]=MixX;
1683           MixY[maxBin][0]=MixY;
1684
1685           //cout << s << " " << b << " " << bin  << " " << maxBin << endl;
1686         }
1687       }
1688     }
1689     //sumSig=0, sumBg=0, esumSig=0, esumBg=0, sumMin=0, esumMin=0;
1690     float sumPhi=0,sumHS=0,sumSS=0,sumMix=0,sumMixXY=0, sumMin=0, esumMin=0;
1691     for(int j=0;j<nMin;j++){
1692       sumMin+=MinArray[j][0];
1693       esumMin+=MinArray[j][1];
1694       sumPhi+=Phi[j][0];
1695       sumHS+=(PhiX[j][0]*MixY[j][0]+PhiY[j][0]*MixX[j][0]);
1696       sumMix+=mix;
1697       sumSS+=b*ss;
1698       sumMixXY+=(2*MixY[j][0]*MixX[j][0]);
1699     }
1700     // if(sumSig==0){sumSig=1E-10; cout << "Warning: Zero sumSig" << endl;}
1701     // if(sumBg==0){sumBg=1E-10; cout << "Warning: Zero sumBg" << endl;}
1702     float a1,a2;
1703     // cout << pow(sumHS*sumHS-4*sumPhi*(sumMixXY-sumSS),0.5) << endl;
1704     a1=(sumHS+pow(sumHS*sumHS-4*sumPhi*(sumMixXY-sumSS),0.5))/(2*sumPhi*(sumMixXY-sumSS));
1705     a2=(sumHS-pow(sumHS*sumHS-4*sumPhi*(sumMixXY-sumSS),0.5))/(2*sumPhi*(sumMixXY-sumSS));
1706     if(fabs(oldmin)<fabs(sumMin))itersize/=2;
1707     oldmin=sumMin;
1708     if(a1>0&&a2>0){
1709       if(fabs(a1-1)<fabs(a2-1)){
1710         a[0]=a1;
1711       }
1712       else a[0]=a2;
1713     }
1714     else if(a1<0)a[0]=a2;
1715     else if(a2<0)a[0]=a1;
1716     else if(sumMin<0) a[0]-=itersize;
1717     else a[0]+=itersize;
1718
1719
1720
1721     cout << "a1: " << a1 << "    a2: " << a2 << endl;
1722     //scalea[1]=scalea[0]*pow(esumSig/sumSig/sumSig+esumBg/sumBg/sumBg,0.5);
1723     esumMin=pow(esumMin,0.5);
1724     //esumSig=pow(esumSig,0.5);
1725     //esumBg=pow(esumBg,0.5);
1726     cout << "Iter:  " << i << "   Min=" << sumMin/nMin << " +/- " << esumMin/nMin << "  a=" << a[0] << endl;
1727   }
1728 }
1729
1730
1731 //----------------------------------------------------------------------------------------------
1732 void Add1D(TH1F *Histo, TH1F *HistAdd, float xscale){
1733   int nbinsx=Histo->GetNbinsX();
1734   for(int xi=1;xi<=nbinsx;xi++){
1735     if(fabs(Histo->GetBinContent(xi)>1E-10)&&fabs(HistAdd->GetBinContent(xi)>1E-10)){//only add if contents in both
1736       Histo->SetBinContent(xi,Histo->GetBinContent(xi)+xscale*HistAdd->GetBinContent(xi));
1737       Histo->SetBinError(xi,sqrt(pow(Histo->GetBinError(xi),2)+pow(xscale*HistAdd->GetBinError(xi),2)));
1738     }
1739   }
1740 }
1741
1742 //----------------------------------------------------------------------------------
1743 void Add2D(TH2F *Histo, TH2F *HistAdd,float xscale){
1744   int nbinsx=Histo->GetNbinsX();
1745   int nbinsy=Histo->GetNbinsY();
1746   //histOut=(TH2F*)Histo->Clone();
1747   // histOut->SetName("Added Histo");
1748
1749   for(int xi=1;xi<=nbinsx;xi++){
1750     for(int yi=1;yi<=nbinsy;yi++){
1751       //cout << Histo->GetBinContent(xi,yi) << endl;
1752       if(fabs(Histo->GetBinContent(xi,yi)>1E-10)&&fabs(HistAdd->GetBinContent(xi,yi))>1E-10){
1753       Histo->SetBinContent(xi,yi,Histo->GetBinContent(xi,yi)+xscale*HistAdd->GetBinContent(xi,yi));
1754       Histo->SetBinError(xi,yi,sqrt(pow(Histo->GetBinError(xi,yi),2)+pow(xscale*HistAdd->GetBinError(xi,yi),2)));
1755       }
1756     }
1757   }
1758 }
1759 //----------------------------------------------------------------------------------
1760 void DivideMean1D(TH1F *Hist, TH1F *Hist2,float a1, float a2,float ntrig){
1761     int nbinsx=Hist->GetNbinsX();
1762     float bin,err,bin2,err2;
1763     float sigma=(a2-a1)/2;
1764     for(int xi=1;xi<=nbinsx;xi++){
1765       bin=Hist->GetBinContent(xi);
1766       err=Hist->GetBinError(xi);
1767       bin2=Hist2->GetBinContent(xi);
1768       err2=Hist2->GetBinError(xi);
1769       if(bin2==0){
1770         bin2=1;
1771         bin=0;
1772         err2=1;
1773       }
1774       Hist->SetBinContent(xi,bin/bin2);
1775       Hist->SetBinError(xi,sigma/(err2*ntrig*Hist->GetBinWidth(1)));
1776     }
1777
1778 }
1779 //---------------------------------------------------------------------------------
1780 //Correction factor screws with the statistics
1781 void DivideMean1DCorr(TH1F *Hist, TH1F *Hist2, TH1F *Hist3, float a1, float a2,float ntrig){
1782     int nbinsx=Hist->GetNbinsX();
1783     float bin,err,bin2,err2, bin3, err3;
1784     float sigma=(a2-a1)/2;
1785     for(int xi=1;xi<=nbinsx;xi++){
1786       bin=Hist->GetBinContent(xi);
1787       err=Hist->GetBinError(xi);
1788       bin2=Hist2->GetBinContent(xi);
1789       err2=Hist2->GetBinError(xi);
1790       bin3=Hist3->GetBinContent(xi);
1791       err3=Hist3->GetBinError(xi);
1792       if(bin2==0){
1793         bin2=1;
1794         bin=0;
1795         err2=1;
1796         err3=1;
1797       }
1798       Hist->SetBinContent(xi,bin/bin2);
1799       Hist->SetBinError(xi,sigma/(err3*ntrig*Hist->GetBinWidth(1)));//error should go with uncorrected counts
1800     }
1801 }
1802
1803 //-----------------------------
1804 //Background subtraction screws with them more
1805 void DivideMean1DBgSub(TH1F *Hist, TH1F *Hist2,TH1F *Hist3,TH1F *Hist4,float a1, float a2,float ntrig, float nmix){
1806     int nbinsx=Hist->GetNbinsX();
1807     float bin,err,bin2,err2,bin3,err3,bin4,err4;
1808     float sigma=(a2-a1);
1809     float e1,e2;
1810     for(int xi=1;xi<=nbinsx;xi++){
1811       bin=Hist->GetBinContent(xi);
1812       err=Hist->GetBinError(xi);
1813       bin2=Hist2->GetBinContent(xi);
1814       err2=Hist2->GetBinError(xi);
1815       bin3=Hist3->GetBinContent(xi);
1816       err3=Hist3->GetBinError(xi);
1817       bin4=Hist4->GetBinContent(xi);
1818       err4=Hist4->GetBinError(xi);
1819       if(bin2==0){
1820         bin2=1;
1821         bin=0;
1822         err2=1;
1823         err3=1;
1824         err4=1;
1825       }
1826       Hist->SetBinContent(xi,bin/bin2);
1827       cout << sigma << endl;
1828       e1=pow(sigma/(err3*ntrig*Hist->GetBinWidth(xi)),2);
1829       e2=pow(sigma/(err4*nmix*Hist->GetBinWidth(xi)),2);
1830       cout << e1 << endl;
1831
1832       Hist->SetBinError(xi,sqrt(e1));
1833     }
1834 }
1835 //----------------------------------------------------------------------------------
1836 void DivideMean2D(TH2F *Hista, TH2F *Hista2, float a1, float a2,float ntrig){
1837   int nbinsx=Hista->GetXaxis()->GetNbins();
1838   int nbinsy=Hista->GetXaxis()->GetNbins();
1839     float bin,err,bin2,err2;
1840     float sigma=(a2-a1)/2;
1841     for(int xi=1;xi<=nbinsx;xi++){
1842       for(int yi=1;yi<=nbinsy;yi++){
1843         bin=Hista->GetBinContent(xi,yi);
1844         err=Hista->GetBinError(xi,yi);
1845         bin2=Hista2->GetBinContent(xi,yi);
1846         err2=Hista2->GetBinError(xi,yi);
1847         if(bin2==0){
1848           bin2=1;
1849           bin=0;
1850           err2=1;
1851         }
1852
1853         Hista->SetBinContent(xi,yi,bin/bin2);
1854         Hista->SetBinError(xi,yi,sigma/(err2*ntrig*HistGetXaxis()->GetBinWidth(1)*HistGetYaxis()->GetBinWidth(1)));
1855       }
1856     }
1857 }
1858 //----------------------------------------------------------------------------------------------------------------
1859 void SetMargins1D(TCanvas *xc){
1860   xc->SetRightMargin(0.02);
1861   xc->SetLeftMargin(0.16);
1862   xc->SetBottomMargin(0.15);
1863 }
1864
1865 void SetMargins1DCyl(TCanvas *xc){
1866   xc->SetRightMargin(0.1);
1867   xc->SetLeftMargin(0.1);
1868   xc->SetTopMargin(0.05);
1869   xc->SetBottomMargin(0.05);
1870   xc->SetTheta(90);
1871   xc->SetPhi(209);//151,209
1872 }
1873
1874 void SetMargins2DSurf(TCanvas *xc){
1875   xc->SetBottomMargin(0.15);
1876   xc->SetLeftMargin(0.21);
1877   xc->SetRightMargin(0.05);
1878   xc->SetTopMargin(0.05);
1879 }
1880 void SetMargins2D(TCanvas *xc){
1881   xc->SetBottomMargin(0.15);
1882   xc->SetLeftMargin(0.12);
1883   xc->SetRightMargin(0.25);
1884 }
1885 void SetTitles1D(TH1 *xhisto){
1886   xhisto->GetXaxis()->SetTitleColor(1);
1887   xhisto->GetXaxis()->SetTitleSize(0.06);
1888   xhisto->GetYaxis()->SetTitleSize(0.06);
1889   xhisto->GetXaxis()->SetLabelSize(0.05);
1890   xhisto->GetYaxis()->SetLabelSize(0.05);
1891   xhisto->GetYaxis()->SetTitleOffset(1.2);
1892 }
1893 void SetTitles2D(TH1 *xhisto){
1894   xhisto->GetXaxis()->SetTitleColor(1);
1895   xhisto->GetXaxis()->SetTitleSize(0.06);
1896   xhisto->GetYaxis()->SetTitleSize(0.06);
1897   xhisto->GetZaxis()->SetTitleSize(0.06);
1898   xhisto->GetXaxis()->SetLabelSize(0.05);
1899   xhisto->GetYaxis()->SetLabelSize(0.05);
1900   xhisto->GetZaxis()->SetLabelSize(0.05);
1901   xhisto->GetZaxis()->SetTitleOffset(1.43);
1902 }
1903
1904 void Style(int i)
1905 {
1906   switch (i) {
1907   case 1:
1908     gStyle->SetStatColor(kWhite);
1909     gStyle->SetTitleColor(kWhite);
1910     gStyle->SetPadColor(kWhite);
1911     gStyle->SetCanvasColor(kWhite);
1912     gStyle->SetPadBorderMode(0);
1913     gStyle->SetPadBorderSize(0);
1914     gStyle->SetCanvasBorderMode(0);
1915     gStyle->SetCanvasBorderSize(0);
1916     gStyle->SetFrameBorderMode(0);
1917     gStyle->SetFrameBorderSize(1);
1918     gStyle->SetFrameFillStyle(0);
1919     gStyle->SetOptStat(0);
1920     gStyle->SetOptFit(0);
1921     gStyle->SetPalette(1);
1922     break;
1923   case 2:  //paper figure style
1924     gStyle->SetStatColor(kWhite);
1925     gStyle->SetTitleColor(kWhite);
1926     gStyle->SetPadColor(kWhite);
1927     gStyle->SetCanvasColor(kWhite);
1928     gStyle->SetPadBorderMode(0);
1929     gStyle->SetPadBorderSize(0);
1930     gStyle->SetCanvasBorderMode(0);
1931     gStyle->SetCanvasBorderSize(0);
1932     gStyle->SetFrameBorderMode(0);
1933     gStyle->SetFrameBorderSize(1);
1934     gStyle->SetFrameFillStyle(0);
1935     gStyle->SetOptStat(0);
1936     gStyle->SetOptFit(1);
1937     gStyle->SetPalette(1);
1938     break;
1939   default:
1940     break;
1941   }
1942 }
1943 //------------------------------------------------------------------
1944 void keySymbol(Float_t x, Float_t y, const Char_t* tt,
1945                Int_t color=1, Int_t marker=1, Float_t tSize=0.04, Float_t Msize=1.0){
1946
1947   gPad->Update();
1948
1949   TLatex *t = new TLatex();
1950   t->SetNDC();
1951   t->SetTextAlign(12);
1952   t->SetTextSize(tSize);
1953   t->SetTextColor(color);
1954   t->DrawLatex(x+0.025,y,tt);
1955
1956
1957   TMarker *l = new TMarker();
1958   l->SetNDC(kTRUE);
1959   l->SetMarkerStyle(marker);
1960   l->SetMarkerColor(color);
1961   l->SetMarkerSize(Msize);
1962   l->PaintMarkerNDC(x-0.0025,y);
1963   l->SetX(x-0.0025);
1964   l->SetY(y);
1965   l->Draw();
1966 }
1967 //------------------------------------------------
1968 void keyText(Float_t x, Float_t y, const Char_t* tt,
1969              Int_t color=1,  Float_t tSize=0.04, Float_t tangle=0){
1970
1971   gPad->Update();
1972   Float_t x1 = gPad->GetFrame()->GetX1();
1973   Float_t x2 = gPad->GetFrame()->GetX2();
1974   Float_t dx = x2-x1;
1975   x = x1 + dx*x;
1976   Float_t y1 = gPad->GetFrame()->GetY1();
1977   Float_t y2 = gPad->GetFrame()->GetY2();
1978   Float_t dy = y2-y1;
1979   y = y1 + dy*y;
1980 TLatex *t = new TLatex(x,y,tt);
1981   t->SetTextAlign(12);
1982   t->SetTextSize(tSize);
1983   t->SetTextColor(color);
1984   t->SetTextAngle(tangle);
1985   t->Draw();
1986 }
1987 ////////////////////////////////////////////////
1988 void rebin1D(TH1 *&temp1D, int rbins){
1989   int nbinxx1=temp1D->GetNbinsX();
1990   int nbinxx2=nbinxx1/rbins;
1991   float widxx=temp1D->GetBinWidth(1);
1992   float minxx=temp1D->GetBinCenter(1)-widxx/2;
1993   float maxxx=temp1D->GetBinCenter(nbinxx1)+widxx/2;
1994   char *titxx=temp1D->GetTitle();
1995
1996   TH1D *temp1D2=new TH1D("temp1D2",titxx,nbinxx2,minxx,maxxx);
1997   for(int xx1=1;xx1<=nbinxx2;xx1++){
1998     float conx1=0,errx1=0;
1999     for(int xx2=1;xx2<=rbins;xx2++){
2000       conx1+=temp1D->GetBinContent(rbins*(xx1-1)+xx2);
2001
2002       errx1+=pow(temp1D->GetBinError(rbins*(xx1-1)+xx2),2);
2003
2004     }
2005
2006     conx1/=rbins;
2007     errx1=sqrt(errx1)/rbins;
2008     temp1D2->SetBinContent(xx1,conx1);
2009     temp1D2->SetBinError(xx1,errx1);
2010   }
2011   char *namex=temp1D->GetName();
2012   char *xtit=temp1D->GetXaxis()->GetTitle();
2013   char *ytit=temp1D->GetYaxis()->GetTitle();
2014   int markx=temp1D->GetMarkerStyle();
2015   int colx=temp1D->GetMarkerColor();
2016   temp1D=(TH1D*)temp1D2->Clone();
2017   temp1D->SetName(namex);
2018   temp1D->GetXaxis()->SetTitle(xtit);
2019   temp1D->GetYaxis()->SetTitle(ytit);
2020   temp1D->SetMarkerStyle(markx);
2021   temp1D->SetLineColor(colx);
2022   temp1D->SetMarkerColor(colx);
2023   SetTitles1D(temp1D);
2024
2025 }
2026 ////////////////////////////////
2027 //Broken need to fix before I use
2028 void rebin2D(TH2D *temp2D,int rbins){
2029   int nbinxx1=temp2D->GetNbinsX();
2030   int nbinxx2=nbinxx2/rbins;
2031   float widxx=temp2D->GetBinWidth(1);
2032   float minxx=temp2D->GetBinCenter(1)-widxx/2;
2033   float maxxx=temp2D->GetBinCenter(nbinxx1)+widxx/2;
2034   char titxx=temp2D->GetTitle();
2035   TH2D *temp2D2=new TH2D("temp2D2",titxx,nbinxx2,minxx,maxxx,nbinxx2,minxx,maxxx);
2036   for(int xx1=1;xx1<=nbinxx2;xx1++){
2037     for(int yy1=1;yy1<=nbinxx2;yy1++){
2038       float conx1=0,errx1=0;
2039       for(int xx2=1;xx2<=rbins;xx2++){
2040         for(int yy2=1;yy2<=rbins;yy2++){
2041           conx1+=temp2D->GetBinContent(rbins*(xx1-1)+xx2,rebins*(yy1-1)+yy2);
2042           errx1+=pow(temp2D->GetBinError(rbins(xx1-1)+xx2,rebins*(yy1-1)+yy2),2);
2043         }
2044       }
2045       conx1/=pow(rbins,2);
2046       errx1=sqrt(errx1)/pow(rbins,2);
2047       temp2D2->SetBinContent(xx1,conx1);
2048       temp2D2->SetBinError(xx1,errx1);
2049     }
2050   }
2051   return temp2D2;
2052 }
2053 //----------------------------------
2054 void DrawLego1D(TCanvas *&tc, TH1 *&thist){
2055
2056   //cout << "off" << thist->GetXaxis()->GetLabelOffset() << endl;
2057   if(thist->GetMinimum()>0)thist->SetMinimum(0);
2058   thist->GetXaxis()->SetLabelOffset(-0.1);
2059   tc1=new TPad("tc1","ytit",0,0.0,0.1,1);
2060   tc2=new TPad("tc2","main",0.1,0,1,1);
2061   tc1->Draw();
2062   tc2->Draw();
2063   tc2->cd();
2064   tc2->SetTheta(0);
2065   tc2->SetPhi(0);
2066   tc2->SetRightMargin(0.02);
2067   //xc->SetLeftMargin(0.16);
2068   tc2->SetBottomMargin(0.15);
2069   thist->Draw("lego2");
2070   tc1->cd();
2071   keyText(.45,.65,"#frac{1}{N_{Trig}}#frac{dN}{d#Delta#phi}",1,0.4,90);
2072   // thist->GetXaxis()->SetLabelOffset(0);
2073 }
2074 //---------------------------------------
2075 void SaveAsText(TH1 *histo, char *filename){
2076
2077   //I should also make it works with TH2
2078   FILE *fp;
2079   fp=fopen(filename,"w");
2080   fprintf(fp,"Bin\t%s\tWidth\t%s\tError\n",histo->GetXaxis()->GetTitle(),histo->GetYaxis()->GetTitle());
2081   for(int x=1;x<=histo->GetNbinsX();x++){
2082     fprintf(fp,"%d\t%3.4f\t%3.4f\t%3.4f\t%3.4f\n",x,histo->GetBinCenter(x),histo->GetBinWidth(x),histo->GetBinContent(x),histo->GetBinError(x));
2083   }
2084   fclose(fp);
2085 }
2086
2087 //==========================================
2088 void MixedCorrect(TH1F *hPhi, TH1F *hPhiMix, TH1F *hEtaN, TH1F *hEtaNMix, TH1F *hEtaA, TH1F *hEtaAMix, TH2F *hPhiEta, TH2F *hPhiEtaMix){
2089
2090 //Do not use GeoCorr if using this
2091   // float avePhiMix=hPhiMix->GetSum()/hPhiMix->GetNbinsX();
2092   float avePhiMix=0, eavePhiMix=0;
2093   for(int i=0;i<=hPhiMix->GetNbinsX();i++){
2094     avePhiMix+=hPhiMix->GetBinContent(i);
2095     eavePhiMix+=pow(hPhiMix->GetBinError(i),2);
2096   }
2097   avePhiMix/=hPhiMix->GetNbinsX();
2098   eavePhiMix=pow(eavePhiMix,0.5)/hPhiMix->GetNbinsX();
2099   hPhiMix->Scale(1/avePhiMix);
2100   hPhi->Divide(hPhiMix);
2101
2102
2103   float maxEtaNMix=hEtaNMix->GetBinContent(hEtaNMix->GetNbinsX()/2);
2104   // float maxEtaNMix=0, emaxEtaNMix=0;
2105   // for(int i=0;i<=hEtaNMix->GetNbinsX();i++){
2106   //   aveEtaMix+=hEtaMix->GetBinContent(i);
2107   //   eaveEtaMix+=pow(hEtaMix->GetBinError(i),2);
2108   //  }
2109   //aveEtaMix/=hEtaMix->GetNbinsX();
2110   // eaveEtaMix=pow(eaveEtaMix,0.5)/hEtaMix->GetNbinsX();
2111   hEtaNMix->Scale(1/maxEtaNMix);
2112   hEtaN->Divide(hEtaNMix);
2113
2114   float maxEtaAMix=hEtaAMix->GetBinContent(hEtaAMix->GetNbinsX()/2);
2115   hEtaAMix->Scale(1/maxEtaAMix);
2116   hEtaA->Divide(hEtaAMix);
2117
2118   float maxavePhiEta=0,emaxavePhiEta=0;
2119   float nbin=0;
2120   for(int i=1;i<=hPhiEtaMix->GetNbinsX();i++){
2121     maxavePhiEta+=hPhiEtaMix->GetBinContent(i,hPhiEtaMix->GetNbinsY()/2);
2122     emaxavePhiEta+=pow(hPhiEtaMix->GetBinError(i,hPhiEtaMix->GetNbinsY()/2),2);
2123   }
2124   maxavePhiEta/=hPhiEtaMix->GetNbinsX();
2125   emaxavePhiEta=pow(emaxavePhiEta,0.5)/hPhiEtaMix->GetNbinsX();
2126   hPhiEtaMix->Scale(1/maxavePhiEta);
2127   hPhiEta->Divide(hPhiEtaMix);
2128
2129   for(int i=1;i<=hPhiMix->GetNbinsX();i++){
2130     hPhiMix->SetBinContent(i,avePhiMix);
2131     hPhiMix->SetBinError(i,0);
2132     }
2133
2134   for(int i=1;i<=hEtaNMix->GetNbinsX();i++){
2135     hEtaNMix->SetBinContent(i,maxEtaNMix);
2136     hEtaAMix->SetBinContent(i,maxEtaAMix);
2137     hEtaNMix->SetBinError(i,0);
2138     hEtaAMix->SetBinError(i,0);
2139   }
2140
2141   for(int i=1;i<=hPhiEta->GetNbinsX();i++){
2142     for(int j=1;j<=hPhiEta->GetNbinsY();j++){
2143       hPhiEtaMix->SetBinContent(i,j,maxavePhiEta);
2144       hPhiEtaMix->SetBinError(i,j,0);
2145     }
2146   }
2147 }
2148
2149 void MixedCorrect3(TH2F *hPhiPhi, TH2F *hPhiPhiSS, TH2F *hPhiPhiMix,TH2F *hEtaEta, TH2F *hEtaEtaSS, TH2F *hEtaEtaMix){
2150   float avePhiPhiMix=hPhiPhiMix->GetSum()/pow(hPhiPhiMix->GetNbinsX(),2);
2151   cout << "hPhiPhiMix: " << hPhiPhiMix->GetBinContent(1,1) << " " << hPhiPhiMix->GetBinError(1,1) << endl;
2152   hPhiPhiMix->Scale(1/avePhiPhiMix);
2153   cout << "hPhiPhiMix: " << hPhiPhiMix->GetBinContent(1,1) << " " << hPhiPhiMix->GetBinError(1,1) << endl;
2154   hPhiPhi->Divide(hPhiPhiMix);
2155   hPhiPhiSS->Divide(hPhiPhiMix);
2156
2157   float maxEtaEtaMix=hEtaEtaMix->GetBinContent(hEtaEtaMix->GetNbinsX()/2,hEtaEtaMix->GetNbinsX()/2);
2158   hEtaEtaMix->Scale(1/maxEtaEtaMix);
2159   hEtaEta->Divide(hEtaEtaMix);
2160   hEtaEtaSS->Divide(hEtaEtaMix);
2161
2162   float conX, conY;
2163   int nbin=hEtaEta->GetNbinsX();
2164   float EtaCut=hEtaEta->GetBinCenter(nbin)+hEtaEta->GetBinWidth(1)/2.;
2165
2166   for(int i=1;i<=hPhiPhi->GetNbinsX();i++){
2167
2168     for(int j=0;j<=hPhiPhi->GetNbinsX();j++){
2169       hPhiPhiMix->SetBinContent(i,j,avePhiPhiMix);
2170       hPhiPhiMix->SetBinError(i,j,0);
2171     }
2172   }
2173   for(int i=1;i<=hEtaEta->GetNbinsX();i++){
2174     conX=hEtaEta->GetBinCenter(i);
2175     for(int j=0;j<=hEtaEta->GetNbinsX();j++){
2176       conY=hEtaEta->GetBinCenter(j);
2177       if(hEtaEtaMix->GetBinContent(i,j)){
2178         hEtaEtaMix->SetBinContent(i,j,maxEtaEtaMix);
2179         hEtaEtaMix->SetBinError(i,j,0);
2180       }
2181     }
2182   }
2183 }
2184 //=================================================================
2185
2186 void ProjPhiPhi(TH2F *shist, TH1F *Non1D, TH1F *Noff1D, TH1F *Aon1D, TH1F *Aoff1D , Float_t NearCut=1.5, float DiagWidth=0.35,int nrand=1000){
2187   delete gRandom;
2188   float pi=3.14159265;
2189   int binxy=shist->GetNbinsX();
2190   float widthxy=shist->GetXaxis()->GetBinWidth(1);
2191   gRandom=new TRandom3();
2192   gRandom->SetSeed(1);
2193   float s,d;
2194   int Nbinxy=2.8*NearCut/widthxy;
2195   float AwayCut=pi-NearCut;
2196   int Abinxy=2.8*AwayCut/widthxy;
2197   *Non1D=new TH1F("Non1D","Near On-Diag",Nbinxy,-NearCut,NearCut);
2198   *Noff1D=new TH1F("Noff1D","Near Off-Diag",Nbinxy,-NearCut,NearCut);
2199   *Aon1D=new TH1F("Aon1D","Away On-Diag",Abinxy,-AwayCut,AwayCut);
2200   *Aoff1D=new TH1F("Aoff1D","Away Off-Diag",Abinxy,-AwayCut,AwayCut);
2201
2202   Non1D->GetXaxis()->SetTitle("#Sigma=(#Delta#phi_{1}+#Delta#phi_{2})/2  or #Delta=(#Delta#phi_{1}-#Delta#phi_{2})/2     ");
2203   Noff1D->GetXaxis()->SetTitle("#Sigma=(#Delta#phi_{1}+#Delta#phi_{2})/2 or #Delta=(#Delta#phi_{1}-#Delta#phi_{2})/2     ");
2204   Aon1D->GetXaxis()->SetTitle("#Sigma=(#Delta#phi_{1}+#Delta#phi_{2})/2-#pi or #Delta=(#Delta#phi_{1}-#Delta#phi_{2})/2     ");
2205   Aoff1D->GetXaxis()->SetTitle("#Sigma=(#Delta#phi_{1}+#Delta#phi_{2})/2-#pi or #Delta=(#Delta#phi_{1}-#Delta#phi_{2})/2     ");
2206   Non1D->GetYaxis()->SetTitle("1/N_{Trig} dN/d(#Sigma or #Delta)   ");
2207   Noff1D->GetYaxis()->SetTitle("1/N_{Trig} dN/d(#Sigma or #Delta)   ");
2208   Aon1D->GetYaxis()->SetTitle("1/N_{Trig} dN/d(#Sigma or #Delta)   ");
2209   Aoff1D->GetYaxis()->SetTitle("1/N_{Trig} dN/d(#Sigma or #Delta)   ");
2210   SetTitles1D(Non1D);
2211   SetTitles1D(Noff1D);
2212   SetTitles1D(Aon1D);
2213   SetTitles1D(Aoff1D);
2214
2215   TH1F *eNon1D=(TH1F*)Non1D->Clone();
2216   eNon1D->SetName("eNon1D");
2217   TH1F *eNoff1D=(TH1F*)Noff1D->Clone();
2218   eNoff1D->SetName("eNoff1D");
2219   TH1F *eAon1D=(TH1F*)Aon1D->Clone();
2220   eAon1D->SetName("eAon1D");
2221   TH1F *eAoff1D=(TH1F*)Aoff1D->Clone();
2222   eAoff1D->SetName("eAoff1D");
2223
2224   //float binsw=2*pi/binxy*2*pi/(2*pi-2);
2225   float binsw;
2226   float NbinsW=Non1D->GetBinWidth(1);
2227   float AbinsW=Aon1D->GetBinWidth(1);
2228   float x10,y10,x20,y20,conp1,errp1;
2229   int nears=10;
2230   for(int xx1=1;xx1<=binxy;xx1++){
2231     x10=shist->GetBinCenter(xx1);
2232     for(int yy1=1;yy1<=binxy;yy1++){
2233       y10=shist->GetBinCenter(yy1);
2234       if((x10<NearCut&&y10<NearCut)){binsw=NbinsW/widthxy; nears=1;}
2235       else if((x10>NearCut&&y10>NearCut)){binsw=AbinsW/widthxy;nears=0;}
2236       else continue;
2237       conp1=shist->GetBinContent(xx1,yy1)*binsw/nrand;
2238       errp1=pow(shist->GetBinError(xx1,yy1)*binsw,2)/nrand;
2239       //     cout << conp1 << " " << errp1 << " " << endl;
2240       // cout << shist->GetBinContent(xx1,yy1) << " " << shist->GetBinError(xx1,yy1) << endl;
2241       for(int r1=0;r1<nrand;r1++){
2242         x20=x10+binsw*(gRandom->Rndm()-0.5);
2243         y20=y10+binsw*(gRandom->Rndm()-0.5);
2244         s=(x20+y20)/2;
2245         d=(x20-y20)/2;
2246         if(fabs(s-pi)<DiagWidth&&nears==0){
2247           Aoff1D->Fill(d,conp1);
2248           eAoff1D->Fill(d,errp1);
2249
2250         }
2251
2252         if(fabs(s)<DiagWidth&&nears==1){
2253           Noff1D->Fill(d,conp1);
2254           eNoff1D->Fill(d,errp1);
2255         }
2256         if(d<DiagWidth&&d>0){
2257           if(nears==0){
2258             Aon1D->Fill((s-pi),2*conp1);//conp1*2
2259             eAon1D->Fill((s-pi),4*errp1);//errp1*4
2260           }
2261           else if(nears==1){
2262             Non1D->Fill(s,2*conp1);//conp1*2
2263             eNon1D->Fill(s,4*errp1);//errp1*4
2264           }
2265         }
2266       }
2267     }
2268   }
2269   for(int xx1=1;xx1<=Nbinxy;xx1++){
2270     Noff1D->SetBinError(xx1,sqrt(eNoff1D->GetBinContent(xx1)));
2271     Non1D->SetBinError(xx1,sqrt(eNon1D->GetBinContent(xx1)));
2272   }
2273   for(int xx1=1;xx1<=Abinxy;xx1++){
2274     Aoff1D->SetBinError(xx1,sqrt(eAoff1D->GetBinContent(xx1)));
2275     Aon1D->SetBinError(xx1,sqrt(eAon1D->GetBinContent(xx1)));
2276   }
2277   Non1D->SetMarkerStyle(25);
2278   Non1D->SetMarkerColor(4);
2279   Non1D->SetLineColor(4);
2280   Noff1D->SetMarkerStyle(20);
2281   Noff1D->SetMarkerColor(2);
2282   Noff1D->SetLineColor(2);
2283  Aon1D->SetMarkerStyle(25);
2284  Aon1D->SetMarkerColor(4);
2285  Aon1D->SetLineColor(4);
2286  Aoff1D->SetMarkerStyle(20);
2287   Aoff1D->SetMarkerColor(2);
2288   Aoff1D->SetLineColor(2);
2289 }
2290 //========================================
2291 ProjEtaEta(TH2F *shist, TH1F *Non1D, TH1F *Noff1D, float DiagWidth=0.35,int nrand=1000){
2292   float pi=3.14159265;
2293   int binxy=shist->GetNbinsX();
2294   float widthxy=shist->GetXaxis()->GetBinWidth(1);
2295   gRandom=new TRandom3();
2296   gRandom->SetSeed(1);
2297   float NearCut=shist->GetBinCenter(binxy)+shist->GetBinWidth(binxy)/2.;
2298   float s,d;
2299   int Nbinxy=2.8*NearCut/widthxy;
2300
2301   *Non1D=new TH1F("NEon1D","Near On-Diag",Nbinxy,-NearCut,NearCut);
2302   *Noff1D=new TH1F("NEoff1D","Near Off-Diag",Nbinxy,-NearCut,NearCut);
2303
2304   Non1D->GetXaxis()->SetTitle("#Sigma=(#Delta#phi_{1}+#Delta#phi_{2})/2  or #Delta=(#Delta#phi_{1}-#Delta#phi_{2})/2     ");
2305   Noff1D->GetXaxis()->SetTitle("#Sigma=(#Delta#phi_{1}+#Delta#phi_{2})/2 or #Delta=(#Delta#phi_{1}-#Delta#phi_{2})/2     ");
2306   Non1D->GetYaxis()->SetTitle("1/N_{Trig} dN/d(#Sigma or #Delta)   ");
2307   Noff1D->GetYaxis()->SetTitle("1/N_{Trig} dN/d(#Sigma or #Delta)   ");
2308   SetTitles1D(Non1D);
2309   SetTitles1D(Noff1D);
2310
2311  TH1F *eNon1D=(TH1F*)Non1D->Clone();
2312   eNon1D->SetName("eNEon1D");
2313   TH1F *eNoff1D=(TH1F*)Noff1D->Clone();
2314   eNoff1D->SetName("eNEoff1D");
2315
2316   float binsw=Non1D->GetBinWidth(1);
2317   float x10,y10,x20,y20,conp1,errp1;
2318
2319   for(int xx1=1;xx1<=binxy;xx1++){
2320     x10=shist->GetBinCenter(xx1);
2321     for(int yy1=1;yy1<=binxy;yy1++){
2322       y10=shist->GetBinCenter(yy1);
2323       conp1=shist->GetBinContent(xx1,yy1)*binsw/nrand;
2324       errp1=pow(shist->GetBinError(xx1,yy1)*binsw,2)/nrand;
2325       for(int r1=0;r1<nrand;r1++){
2326         x20=x10+binsw*(gRandom->Rndm()-0.5);
2327         y20=y10+binsw*(gRandom->Rndm()-0.5);
2328         s=(x20+y20)/2;
2329         d=(x20-y20)/2;
2330
2331         if(fabs(s)<DiagWidth){
2332           Noff1D->Fill(d,conp1);
2333           eNoff1D->Fill(d,errp1);
2334         }
2335         if(d<DiagWidth&&d>0){
2336           Non1D->Fill(s,2*conp1);//conp1*2
2337           eNon1D->Fill(s,4*errp1);//errp1*4
2338         }
2339       }
2340     }
2341   }
2342  for(int xx1=1;xx1<=Nbinxy;xx1++){
2343     Noff1D->SetBinError(xx1,sqrt(eNoff1D->GetBinContent(xx1)));
2344     Non1D->SetBinError(xx1,sqrt(eNon1D->GetBinContent(xx1)));
2345   }
2346 Non1D->SetMarkerStyle(25);
2347   Non1D->SetMarkerColor(4);
2348   Non1D->SetLineColor(4);
2349   Noff1D->SetMarkerStyle(20);
2350   Noff1D->SetMarkerColor(2);
2351   Noff1D->SetLineColor(2);
2352 }