PWGCF/Correlations/macros/minijet/analyse_pA.C

   1 /* $Id:  $ */
   2 //--------------------------------------------------
   3 //
   4 // macro to do the final analysis step
   5 // uses input of analysis class AliAnalysisTaskPhiCorrelation
   6 //
   7 //  Author : Emilia Leogrande (University of Utrecht)
   8 //
   9 //-------------------------------------------------
  10
  11 #include <TChain.h>
  12 #include <TList.h>
  13 #include <TTree.h>
  14 #include <TH1F.h>
  15 #include <TH2F.h>
  16 #include <TH3F.h>
  17 #include <THnSparse.h>
  18 #include <TProfile.h>
  19 #include <TCanvas.h>
  20 #include "TRandom.h"
  21 #include "TGraphErrors.h"
  22 #include "TFile.h"
  23 #include "TF1.h"
  24 #include "TMath.h"
  25 #include "TDirectory.h"
  26 #include "TStyle.h"
  27 #include "TROOT.h"
  28 #include "TColor.h"
  29
  30 #include <iostream>
  31 using namespace std;
  32
  33 void analyseEmy2(Bool_t zyam = kTRUE);  // if zyam = kFALSE, fit is used
  34 Double_t fitFunction(Double_t *x ,Double_t *par); // fit function using constant + 3 gaussians
  35 Double_t fitFunction2Gaus(Double_t *x ,Double_t *par); // fit function using constant + 2 gaussians
  36
  37 //input file and mixed event removed file
  38 TFile *fileData=0x0;
  39 TFile *fileDataEMremoved = 0x0;
  40
  41 const int multclass = 20;
  42
  43 TH1D *fDeltaPhiNch[multclass];
  44 TH1D *fDeltaEtaNch[multclass];
  45 TH1D *fSignalDPhi[multclass];
  46 TH1D *fSignalNSDPhi[multclass];
  47 TH1D *fSignalASDPhi[multclass];
  48 TH1D *fRidge1DPhi[multclass];
  49 TH1D *fRidge2DPhi[multclass];
  50 TH1D *fRidgeDPhi[multclass];
  51 TH1D *fSymmRidgeNotScaled[multclass];
  52 TH1D *fSymmRidge[multclass];
  53 TH1D *fFinal1DPhi[multclass];
  54 TH1D *fFinalDPhi[multclass];
  55
  56 TString flag = "R";
  57 TF1 *fTotal2Gaus[multclass];       // fit with 2 gaussians + const
  58 TF1 *fTotal[multclass];            // fit with 3 gaussians + const
  59
  60 //properties of histogram
  61 const int bins = 72; //
  62 Double_t binWidth=2*TMath::Pi()/bins;
  63
  64 const int binsDeta = 48;
  65
  66
  67 Double_t max_bin_for_etagap =   1.2;
  68 Double_t min_bin_for_etagap =   -1.2;
  69 Double_t max_eta = 1.8;
  70 Double_t min_eta = -1.8;
  71
  72 //________________________________________________________________________________________________________________
  73 //
  74 Double_t fitFunction(Double_t *x ,Double_t *par)
  75 {
  76   // fit function for 3 gaus + constant
  77
  78   // parameters for Gaussian
  79   Double_t A1     = par[0];
  80   Double_t sigma1 = par[1];
  81   Double_t A2     = par[2];
  82   Double_t sigma2 = par[3];
  83   Double_t A3     = par[4];
  84   Double_t sigma3 = par[5];
  85   Double_t integral = par[6];
  86
  87   Double_t constante = (integral-
  88                         TMath::Sqrt(TMath::Pi()*2)/ binWidth*
  89                         (A1 * sigma1 + A2 * sigma2 + A3*sigma3))/bins;
  90   Double_t q  = x[0];
  91
  92   //fit value
  93   Double_t fitval = constante +
  94     (q>-0.5*TMath::Pi()&&q<0.5*TMath::Pi())*(
  95                                              A1 * exp(- q * q / (2 * sigma1 *sigma1)) +
  96                                              A1 * exp(-((q - TMath::TwoPi())) * ((q - TMath::TwoPi())) / ( 2 * sigma1 * sigma1))
  97                                              )
  98     +
  99     (q>-0.2*TMath::Pi()&&q<0.2*TMath::Pi())*(
 100                                              A2 * exp(- q * q / (2 * sigma2 *sigma2)) +
 101                                              A2 * exp(-((q - TMath::TwoPi())) * ((q - TMath::TwoPi())) / ( 2 * sigma2 * sigma2))
 102                                              )
 103     +
 104     (q>0.5*TMath::Pi()&&q<1.5*TMath::Pi())*(
 105                                             A3 * exp(-((q - TMath::Pi())) * ((q - TMath::Pi())) / ( 2 * sigma3 * sigma3)) +
 106                                             A3 * exp(-((q + TMath::Pi())) * ((q + TMath::Pi())) / (2 * sigma3 * sigma3))
 107                                             );
 108   return fitval;
 109 }
 110
 111 //________________________________________________________________________________________________________________
 112 //
 113 Double_t fitFunction2Gaus(Double_t *x ,Double_t *par)
 114 {
 115   // fit function for 2 gaus + constant
 116
 117   // parameters for Gaussian
 118   Double_t A1     = par[0];
 119   Double_t sigma1 = par[1];
 120   Double_t A3     = par[2];
 121   Double_t sigma3 = par[3];
 122   Double_t integral = par[4];
 123
 124   Double_t constante = (integral -
 125                         TMath::Sqrt(TMath::Pi()*2)/ binWidth*
 126                         (A1 * sigma1 + A3*sigma3))/bins;
 127   Double_t q  = x[0];
 128
 129   //fit value
 130   Double_t fitval = constante +
 131     (q>-0.5*TMath::Pi()&&q<0.5*TMath::Pi())*(
 132                                              A1 * exp(- q * q / (2 * sigma1 *sigma1)) +
 133                                              A1 * exp(-((q - TMath::TwoPi())) * ((q - TMath::TwoPi())) / ( 2 * sigma1 * sigma1))
 134                                              )
 135     +
 136     (q>0.5*TMath::Pi()&&q<1.5*TMath::Pi())*(
 137                                             A3 * exp(-((q - TMath::Pi())) * ((q - TMath::Pi())) / ( 2 * sigma3 * sigma3)) +
 138                                             A3 * exp(-((q + TMath::Pi())) * ((q + TMath::Pi())) / (2 * sigma3 * sigma3))
 139                                             );
 140   return fitval;
 141 }
 142
 143 //_______________________________________________________________________________________________________________
 144 //
 145 Double_t fline(Double_t *x, Double_t *par){
 146
 147     if(x[0]>-1.8 && x[0]<=0){
 148         return par[0]+par[1]*x[0];
 149     }
 150     else if(x[0]>0 && x[0]<1.8){
 151         return par[2]+par[3]*x[0];
 152     }
 153     else
 154         return 0;
 155 }
 156
 157
 158 //________________________________________________________________________________________________________________
 159 //
 160 void analyseEmy2(Bool_t zyam){
 161
 162
 163   // plot style
 164   gStyle->SetOptStat(0);
 165   const Int_t NRGBs = 5;
 166   const Int_t NCont = 500;
 167   Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
 168   Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
 169   Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
 170   Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
 171   TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
 172   gStyle->SetNumberContours(NCont);
 173
 174   //style
 175   gROOT->SetStyle("Plain");
 176   gStyle->SetOptStat(0);
 177   gStyle->SetPalette(1);
 178
 179   //-------------- TRIGGERS AND EVENTS
 180
 181   TH2D *dphideta[multclass];
 182   TH1D * trigger = 0x0;
 183   TH1D * event = 0x0;
 184
 185   fileData = TFile::Open("dphi_corr.root");
 186   trigger = (TH1D*)fileData->Get("triggers_0");
 187   event = (TH1D*)fileData->Get("events");
 188
 189   // get average trigger particles per event
 190   TProfile *p0 = (TProfile*)trigger->Clone();
 191   TProfile *p1 = (TProfile*)event->Clone();
 192   p0->Sumw2();
 193   p1->Sumw2();
 194   p0->Divide(p0,p1,1,1,"B");
 195
 196   // copy triggers and events in the new dphi_corr with the Mixed Event removed
 197   TH1D *triggerCopy = 0x0;
 198   TH1D *eventCopy = 0x0;
 199
 200   triggerCopy = (TH1D*)trigger->Clone();
 201   eventCopy = (TH1D*)event->Clone();
 202
 203   fileDataEMremoved = TFile::Open("dphi_corr_MEremoved.root","RECREATE");
 204   triggerCopy->SetName("triggers_0");
 205   triggerCopy->Write();
 206   eventCopy->SetName("events");
 207   eventCopy->Write();
 208   fileDataEMremoved->Close();
 209
 210
 211   //-------------- MIXED EVENT REMOVAL: restores the right number of particles in the detector acceptance but keeps the detector azimuthal unefficiencies corrections and cures the dip in (0,0) from two-trak cuts
 212   // Removing the event mixing: S/M (from dphi_corr) * M (from the triangle)
 213
 214     Double_t triangle_factor[binsDeta]={0};
 215
 216     TH2D *s_over_m[multclass];
 217     TH1D *s_m_deta[multclass];
 218     TH2D *s_over_m_x_m[multclass];
 219
 220     for(Int_t i=0;i<multclass;i++){
 221         s_over_m[i] = (TH2D*)fileData->Get(Form("dphi_0_0_%d",i));
 222         s_m_deta[i] = (TH1D*)s_over_m[i]->ProjectionY()->Clone();
 223         s_over_m_x_m[i] = (TH2D*)s_over_m[i]->Clone();
 224         s_over_m_x_m[i]->Reset();
 225     }
 226
 227
 228     TF1 *f2 = new TF1("f2",fline,min_eta,max_eta,4);
 229
 230     f2->FixParameter(0,1);
 231     f2->FixParameter(1,1/max_eta);
 232     f2->FixParameter(2,1);
 233     f2->FixParameter(3,-1/max_eta);
 234
 235     for(Int_t i=0;i<binsDeta;i++){
 236
 237         triangle_factor[i] = f2->Eval(s_m_deta[0]->GetBinCenter(i+1));
 238
 239     }
 240
 241
 242
 243     //--scale each deta bin of the old TH2 with the triangle_factor[deta]
 244
 245     for(Int_t i=0;i<multclass;i++){
 246         for(Int_t j=0;j<binsDeta;j++){
 247             for(Int_t k=0;k<bins;k++){
 248                     s_over_m_x_m[i] -> SetBinContent(k+1,j+1,(s_over_m[i]->GetBinContent(k+1,j+1))*triangle_factor[j]);
 249                     s_over_m_x_m[i]->SetBinError(k+1,j+1,(s_over_m[i]->GetBinError(k+1,j+1))*triangle_factor[j]);
 250             }
 251         }
 252     }
 253
 254     fileDataEMremoved = TFile::Open("dphi_corr_MEremoved.root","UPDATE");
 255
 256     for(Int_t i=0;i<multclass;i++){
 257
 258         s_over_m_x_m[i]->SetName(Form("dphiNoMixed_%d",i));
 259         s_over_m_x_m[i]->Write();
 260
 261     }
 262
 263
 264
 265     //-------------- DOUBLE RIDGE SUBTRACTION: gets rid of no-jet related components (v3 is still kept => effect added to the systematics)
 266
 267     // the ridge, estimated via an etagap, has to be scaled since it sits on the triangle
 268     Double_t scale_for_ridge_NS = 0, scale_for_ridge_AS = 0;
 269
 270
 271     scale_for_ridge_NS = f2->Integral(min_bin_for_etagap,max_bin_for_etagap)/(f2->Integral(min_eta,min_bin_for_etagap)+f2->Integral(max_bin_for_etagap,max_eta)); //there is etagap in the NS
 272     cout<<"scaling NS:"<<scale_for_ridge_NS<<endl;
 273
 274     scale_for_ridge_AS = f2->Integral(min_eta,max_eta)/(f2->Integral(min_eta,min_bin_for_etagap)+f2->Integral(max_bin_for_etagap,max_eta)); // there is no etagap in the AS
 275     cout<<"scaling AS:"<<scale_for_ridge_AS<<endl;
 276
 277   // Double ridge subtraction
 278
 279   TCanvas *c = new TCanvas();
 280   c->Divide(5,4);
 281
 282   for(Int_t i=0;i<multclass;i++){
 283   c->cd(i+1);
 284
 285
 286       dphideta[i] = (TH2D*)fileDataEMremoved->Get(Form("dphiNoMixed_%d",i));
 287
 288
 289       // phi and eta projections
 290       fDeltaPhiNch[i] = (TH1D*)dphideta[i]->ProjectionX()->Clone();
 291       if(!zyam)
 292           fDeltaPhiNch[i]->Scale(binWidth);    //gaussians include the binwidth, so when using the fit, the histograms must be scaled first
 293       fDeltaPhiNch[i]->Draw();
 294
 295       fDeltaEtaNch[i] = (TH1D*)dphideta[i]->ProjectionY()->Clone();
 296
 297       // signal NS: |DEta|<max_bin_for_etagap; signal AS: |DEta|<max_eta
 298       fSignalNSDPhi[i] = (TH1D*)dphideta[i]->ProjectionX(Form("|DEta|<%f",max_bin_for_etagap),fDeltaEtaNch[i]->FindBin(min_bin_for_etagap+0.0001),fDeltaEtaNch[i]->FindBin(max_bin_for_etagap-0.0001))->Clone();
 299       fSignalASDPhi[i] = (TH1D*)dphideta[i]->ProjectionX(Form("|DEta|<%f",max_eta))->Clone();
 300
 301       fSignalDPhi[i] = (TH1D*)fSignalASDPhi[i]->Clone();
 302       fSignalDPhi[i]->Reset();
 303       fSignalDPhi[i]->Sumw2();
 304
 305       for(Int_t k=0;k<bins/2;k++){
 306           fSignalDPhi[i]->SetBinContent(k+1,fSignalNSDPhi[i]->GetBinContent(k+1));
 307           fSignalDPhi[i]->SetBinError(k+1, fSignalNSDPhi[i]->GetBinError(k+1));
 308       }
 309       for(Int_t k=bins/2;k<bins;k++){
 310           fSignalDPhi[i]->SetBinContent(k+1,fSignalASDPhi[i]->GetBinContent(k+1));
 311           fSignalDPhi[i]->SetBinError(k+1, fSignalASDPhi[i]->GetBinError(k+1));
 312       }
 313       if(!zyam)
 314           fSignalDPhi[i]->Scale(binWidth);
 315
 316       // ridge1 DEta<min_bin_for_etagap
 317       fRidge1DPhi[i] = (TH1D*)dphideta[i]->ProjectionX(Form("DEta<%f",min_bin_for_etagap),1,fDeltaEtaNch[i]->FindBin(min_bin_for_etagap-0.0001))->Clone();
 318       if(!zyam)
 319           fRidge1DPhi[i]->Scale(binWidth);
 320       fRidge1DPhi[i]->SetMarkerColor(kRed);
 321
 322       // ridge2 DEta>max_bin_for_etagap
 323       fRidge2DPhi[i] = (TH1D*)dphideta[i]->ProjectionX(Form("DEta>%f",max_bin_for_etagap),fDeltaEtaNch[i]->FindBin(max_bin_for_etagap+0.0001),fDeltaEtaNch[i]->GetNbinsX())->Clone();
 324       if(!zyam)
 325           fRidge2DPhi[i]->Scale(binWidth);
 326       fRidge2DPhi[i]->SetMarkerColor(kBlue);
 327
 328       // ridge = ridge1 + ridge2
 329       fRidgeDPhi[i] = (TH1D*)fRidge1DPhi[i]->Clone("fRidge");
 330       fRidgeDPhi[i]->Reset();
 331       fRidgeDPhi[i]->Sumw2();
 332       fRidgeDPhi[i]->Add(fRidge1DPhi[i],fRidge2DPhi[i],1,1);
 333       //fRidgeDPhi[i]->Scale(scale_for_ridge);
 334
 335       // symmetrize NS ridge in the AS
 336       fSymmRidgeNotScaled[i] = (TH1D*)fRidgeDPhi[i]->Clone("fSymmRidgeNotScaled");
 337
 338       for(Int_t k=fSymmRidgeNotScaled[i]->GetNbinsX()/2+1;k<=fSymmRidgeNotScaled[i]->GetNbinsX();k++){
 339
 340           fSymmRidgeNotScaled[i]->SetBinContent(k,fSymmRidgeNotScaled[i]->GetBinContent(fSymmRidgeNotScaled[i]->GetNbinsX()+1-k));
 341
 342       }
 343
 344       // scale the symmetrized ridge according to NS or AS
 345       fSymmRidge[i] = (TH1D*)fSymmRidgeNotScaled[i]->Clone("fSymmRidge");
 346
 347       for(Int_t k=0;k<bins/2;k++){
 348           fSymmRidge[i]->SetBinContent(k+1,(fSymmRidgeNotScaled[i]->GetBinContent(k+1))*scale_for_ridge_NS);
 349       }
 350       for(Int_t k=bins/2;k<bins;k++){
 351           fSymmRidge[i]->SetBinContent(k+1,(fSymmRidgeNotScaled[i]->GetBinContent(k+1))*scale_for_ridge_AS);
 352       }
 353
 354
 355       // signal - symmetric ridge
 356
 357       if(zyam){
 358           fFinal1DPhi[i] = new TH1D(Form("fFinal1DPhi[%d]",i),Form("fFinal1DPhi[%d]",i),bins,-0.5*TMath::Pi(),1.5*TMath::Pi());
 359           fFinal1DPhi[i]->Add(fSignalDPhi[i],fSymmRidge[i],1,-1);
 360           fFinal1DPhi[i]->Sumw2();
 361           fFinalDPhi[i] = (TH1D*)fFinal1DPhi[i]->Clone("fFinal"); // zyam: average between the two min values => sum first half of NS in the second half and second half of AS in the first half, so zyam = min/2
 362           fFinalDPhi[i]->Reset();
 363           fFinalDPhi[i]->Sumw2();
 364
 365           for(Int_t k=1;k<=bins/4;k++){
 366               fFinalDPhi[i]->SetBinContent(k,0.);
 367               fFinalDPhi[i]->SetBinContent(k+bins/4,fFinal1DPhi[i]->GetBinContent(k+bins/4)+fFinal1DPhi[i]->GetBinContent(bins/4+1-k));
 368               fFinalDPhi[i]->SetBinError(k+bins/4,TMath::Sqrt(pow(fFinal1DPhi[i]->GetBinError(k+bins/4),2)+pow(fFinal1DPhi[i]->GetBinError(bins/4+1-k),2)));
 369               fFinalDPhi[i]->SetBinContent(k+bins/2,fFinal1DPhi[i]->GetBinContent(k+bins/2)+fFinal1DPhi[i]->GetBinContent(bins+1-k));
 370               fFinalDPhi[i]->SetBinError(k+bins/2,TMath::Sqrt(pow(fFinal1DPhi[i]->GetBinError(k+bins/2),2)+pow(fFinal1DPhi[i]->GetBinError(bins+1-k),2)));
 371               fFinalDPhi[i]->SetBinContent(k+bins/4*3,0.);
 372
 373           }
 374       }
 375
 376       else{
 377
 378           fFinalDPhi[i] = (TH1D*)fSignalDPhi[i]->Clone();
 379           fFinalDPhi[i]->Reset();
 380           fFinalDPhi[i]->Sumw2();
 381           fFinalDPhi[i]->Add(fSignalDPhi[i],fSymmRidge[i],1,-1);
 382       }
 383
 384   }
 385
 386   // store the pair yields in a file (the yields are *not* normalized to the Ntriggers)
 387
 388   TFile* file_yields = 0x0;
 389   if(zyam)
 390       file_yields = TFile::Open("PairYields_zyam.root","RECREATE");
 391   else
 392       file_yields = TFile::Open("PairYields_fit.root","RECREATE");
 393
 394
 395   for(Int_t i=0;i<multclass;i++){
 396       fDeltaEtaNch[i]->SetName(Form("DeltaEta_0_0_%d",i));
 397       fDeltaEtaNch[i]->Write();
 398       fDeltaPhiNch[i]->SetName(Form("Correlation bin %d in dphi",i));
 399       fDeltaPhiNch[i]->Write();
 400       fSignalDPhi[i]->SetName(Form("Signal_0_0_%d",i));
 401       fSignalDPhi[i]->Write();
 402       fRidgeDPhi[i]->SetName(Form("Ridge_0_0_%d",i));
 403       fRidgeDPhi[i]->Write();
 404       fSymmRidgeNotScaled[i]->SetName(Form("Symmetric_Ridge_NotScaled_0_0_%d",i));
 405       fSymmRidgeNotScaled[i]->Write();
 406       fSymmRidge[i]->SetName(Form("Symmetric_Ridge_0_0_%d",i));
 407       fSymmRidge[i]->Write();
 408       fFinalDPhi[i]->SetName(Form("Pure_Signal_0_0_%d",i));
 409       fFinalDPhi[i]->Write();
 410   }
 411   file_yields->Close();
 412
 413   //-------------- CORRELATION OBSERVABLES: per-trigger yields, triggers and uncorrelated seeds
 414
 415   Float_t baseline[multclass]={0};
 416
 417   TGraphErrors *fNearSideIntegral = new TGraphErrors();
 418   fNearSideIntegral->SetName("fNearSideIntegral");
 419   fNearSideIntegral->SetMarkerColor(kGreen+2);
 420   fNearSideIntegral->SetLineColor(kGreen+2);
 421   fNearSideIntegral->SetLineWidth(1);
 422   fNearSideIntegral->SetMarkerStyle(4);
 423
 424   TGraphErrors *fAwaySideIntegral = new TGraphErrors();
 425   fAwaySideIntegral->SetName("fAwaySideIntegral");
 426   fAwaySideIntegral->SetMarkerColor(kBlue);
 427   fAwaySideIntegral->SetLineColor(kBlue);
 428   fAwaySideIntegral->SetLineWidth(1);
 429   fAwaySideIntegral->SetMarkerStyle(4);
 430
 431   TGraphErrors *fBothSideIntegral = new TGraphErrors();
 432   fBothSideIntegral->SetName("fBothSideIntegral");
 433   fBothSideIntegral->SetMarkerColor(kMagenta);
 434   fBothSideIntegral->SetLineColor(kMagenta);
 435   fBothSideIntegral->SetLineWidth(1);
 436   fBothSideIntegral->SetMarkerStyle(4);
 437
 438
 439   TGraphErrors *fNjets = new TGraphErrors();
 440   fNjets->SetName("fNjets");
 441   fNjets->SetMarkerColor(kCyan+2);
 442   fNjets->SetLineColor(kCyan+2);
 443   fNjets->SetLineWidth(1);
 444   fNjets->SetMarkerStyle(4);
 445
 446   TGraphErrors *fTriggerAverage = new TGraphErrors();
 447   fTriggerAverage->SetName("fTriggerAverage");
 448   fTriggerAverage->SetMarkerColor(kBlack);
 449   fTriggerAverage->SetLineColor(kBlack);
 450   fTriggerAverage->SetLineWidth(1);
 451   fTriggerAverage->SetMarkerStyle(4);
 452
 453   Int_t points=0;
 454   Double_t minbin[multclass] = {0};
 455
 456   //  extract information out of dphi histograms
 457   TCanvas * cYields= new TCanvas("cYields", "cYields", 150, 150, 820, 620);
 458   cYields->Divide(5,4);
 459
 460   for(Int_t i=0;i<multclass;i++){
 461   cYields->cd(i+1);
 462
 463
 464   if(zyam) {
 465
 466       if(fFinalDPhi[i]->Integral()>0){
 467           fFinalDPhi[i]->GetXaxis()->SetRange(bins/4+1,bins/4*3);
 468           baseline[i]=fFinalDPhi[i]->GetMinimum()/2;
 469           minbin[i] = fFinalDPhi[i]->GetMinimumBin();
 470           fFinalDPhi[i]->GetXaxis()->UnZoom();
 471
 472           for(Int_t k=0;k<bins;k++){
 473               if(fFinalDPhi[i]->GetBinContent(k+1)!=0)
 474                   fFinalDPhi[i]->SetBinContent(k+1,fFinalDPhi[i]->GetBinContent(k+1)-baseline[i]);
 475               else
 476                   fFinalDPhi[i]->SetBinContent(k+1,0.);
 477           }
 478
 479           fFinalDPhi[i]->DrawClone("");
 480
 481           fFinalDPhi[i]->SetTitle(Form("0.7<p_{T,trig}<5.0 - 0.7<p_{T,assoc}<5.0 - %d-%d %",i*5,(i+1)*5));
 482           fFinalDPhi[i]->SetTitle("1/N_{trig} dN_{assoc}/d#Delta#varphi (rad^{-1})");
 483           //-
 484           Double_t errorNS = 0;
 485           Double_t nearSideResult = (fFinalDPhi[i]->IntegralAndError(0,minbin[i],errorNS,"width"))/trigger->GetBinContent(i+1);
 486           Double_t nearSideError = errorNS/trigger->GetBinContent(i+1);
 487           fNearSideIntegral->SetPoint(points,i, nearSideResult);
 488           fNearSideIntegral->SetPointError(points,0.5,errorNS/trigger->GetBinContent(i+1));
 489           //-
 490
 491           //--
 492           Double_t errorAS = 0;
 493           Double_t awaySideResult = (fFinalDPhi[i]->IntegralAndError(minbin[i],bins,errorAS,"width"))/trigger->GetBinContent(i+1);
 494           Double_t awaySideError = errorAS/trigger->GetBinContent(i+1);
 495           fAwaySideIntegral->SetPoint(points,i, awaySideResult );
 496           fAwaySideIntegral->SetPointError(points,0.5, errorAS/trigger->GetBinContent(i+1));
 497           //--
 498
 499           //---
 500           Double_t bothSideResult = nearSideResult + awaySideResult;
 501           Double_t bothSideError = bothSideResult * TMath::Sqrt(pow(errorNS,2)+pow(errorAS,2))/trigger->GetBinContent(i+1);
 502           fBothSideIntegral->SetPoint(points,i, bothSideResult );
 503           fBothSideIntegral->SetPointError(points,0.5, bothSideError );
 504           //---
 505
 506
 507
 508       }
 509       else{
 510           fNearSideIntegral->SetPoint(points,i, 0);
 511           fAwaySideIntegral->SetPoint(points,i, 0);
 512           fBothSideIntegral->SetPoint(points,i,0);
 513       }
 514       Double_t p0BinContent=p0->GetBinContent(i+1);
 515       Double_t p0BinError=p0->GetBinError(i+1);
 516
 517       //--------
 518       Double_t njets =  p0BinContent/(1+bothSideResult);
 519       Double_t njetsError = njets*TMath::Sqrt(bothSideError*bothSideError/(1+bothSideResult)/(1+bothSideResult)+p0BinError*p0BinError/p0BinContent/p0BinContent);
 520       fNjets->SetPoint(points,i, njets );
 521       fNjets->SetPointError(points,0.5,njetsError );
 522
 523       //-------
 524
 525       fTriggerAverage->SetPoint(points,i, p0BinContent);
 526       fTriggerAverage->SetPointError(points,0.5, p0BinError);
 527
 528   }
 529
 530   else if (!zyam){
 531
 532       if(fFinalDPhi[i]->Integral()>0){
 533
 534           //first fit function: 2 gauss + const
 535           fTotal2Gaus[i] = new TF1(Form("gaus3and2_%d",i), fitFunction2Gaus , -0.5*TMath::Pi(), 1.5*TMath::Pi(), 5);
 536           fTotal2Gaus[i]->SetName(Form("gaus3_%d",i));
 537           fTotal2Gaus[i]->SetParNames ("A1","sigma1","A3", "sigma3");
 538           fTotal2Gaus[i]->SetLineColor(kRed);
 539           fTotal2Gaus[i]->SetLineWidth(2);
 540
 541           baseline[i]=fFinalDPhi[i]->GetMinimum();
 542           Double_t integr_for_const_2 = fFinalDPhi[i]->Integral();
 543
 544           fTotal2Gaus[i]->FixParameter(4,integr_for_const_2);
 545           fTotal2Gaus[i]->SetParameters( fFinalDPhi[i]->GetBinContent(fFinalDPhi[i]->GetXaxis()->FindFixBin(0)) - baseline[i] , 0.6 , fFinalDPhi[i]->GetBinContent(fFinalDPhi[i]->GetXaxis()->FindFixBin(TMath::Pi()))-baseline[i] , 0.6);
 546
 547           fTotal2Gaus[i]->SetParLimits(0, 0, (fFinalDPhi[i]->GetBinContent(fFinalDPhi[i]->GetXaxis()->FindFixBin(0))-baseline[i])*2);
 548           fTotal2Gaus[i]->SetParLimits(1, 0.01, 10);
 549           fTotal2Gaus[i]->SetParLimits(2, 0, (fFinalDPhi[i]->GetBinContent(fFinalDPhi[i]->GetXaxis()->FindFixBin(TMath::Pi()))-baseline[i])*2);
 550           fTotal2Gaus[i]->SetParLimits(3, 0.01, 10);
 551
 552           fTotal2Gaus[i]->SetLineColor(kRed);
 553           fTotal2Gaus[i]->SetLineWidth(2);
 554
 555           fFinalDPhi[i]->Fit(fTotal2Gaus[i],flag);
 556           fFinalDPhi[i]->SetMinimum(0);
 557           fFinalDPhi[i]->DrawClone("");
 558           fTotal2Gaus[i] ->DrawClone("same");
 559
 560           Double_t A11     = fTotal2Gaus[i]->GetParameter(0);
 561           Double_t sigma11 = fTotal2Gaus[i]->GetParameter(1);
 562           Double_t A31     = fTotal2Gaus[i]->GetParameter(2);
 563           Double_t sigma31 = fTotal2Gaus[i]->GetParameter(3);
 564
 565           Double_t a1e1 = fTotal2Gaus[i]->GetParError(0);
 566           Double_t s1e1 = fTotal2Gaus[i]->GetParError(1);
 567           Double_t a3e1 = fTotal2Gaus[i]->GetParError(2);
 568           Double_t s3e1 = fTotal2Gaus[i]->GetParError(3);
 569
 570
 571           Double_t T11 = A11*sigma11;
 572           Double_t T31 = A31*sigma31;
 573           Double_t t11 = T11*TMath::Sqrt(a1e1*a1e1/A11/A11 + s1e1*s1e1/sigma11/sigma11);
 574           Double_t t31 = T31*TMath::Sqrt(a3e1*a3e1/A31/A31 + s3e1*s3e1/sigma31/sigma31);
 575
 576
 577           //second fit: 3 gauss + const
 578           fTotal[i] = new TF1(Form("gaus3_%d",i), fitFunction , -0.5*TMath::Pi(), 1.5*TMath::Pi(), 7);
 579           fTotal[i]->SetName(Form("gaus3_%d",i));
 580           fTotal[i]->SetParNames ("A1","sigma1","A2","sigma2", "A3", "sigma3","integral");
 581           fTotal[i]->SetLineColor(kRed);
 582           fTotal[i]->SetLineWidth(2);
 583
 584           Double_t integr_for_const = fFinalDPhi[i]->Integral();
 585
 586
 587           fTotal[i]->FixParameter(0,A11);
 588           fTotal[i]->FixParameter(1,sigma11*1.2);
 589           fTotal[i]->FixParameter(2,A11);
 590           fTotal[i]->FixParameter(3,sigma11*0.7);
 591           fTotal[i]->FixParameter(4,A31);
 592           fTotal[i]->FixParameter(5,sigma31);
 593           fTotal[i]->FixParameter(6,integr_for_const);
 594
 595           fTotal[i]->SetParLimits(0, 0, (fFinalDPhi[i]->GetBinContent(fFinalDPhi[i]->GetXaxis()->FindFixBin(0))-baseline[i])*2);
 596           fTotal[i]->SetParLimits(1, 0.3, 10);
 597           fTotal[i]->SetParLimits(2, 0, (fFinalDPhi[i]->GetBinContent(fFinalDPhi[i]->GetXaxis()->FindFixBin(0))-baseline[i])*2);
 598           fTotal[i]->SetParLimits(3, 0.12, 0.4);
 599           fTotal[i]->SetParLimits(4, 0, (fFinalDPhi[i]->GetBinContent(fFinalDPhi[i]->
 600                                                               GetXaxis()->FindFixBin(TMath::Pi()))-baseline[i])*2);
 601           fTotal[i]->SetParLimits(5, 0.01, 10);
 602
 603           fTotal[i]->SetLineColor(kRed);
 604           fTotal[i]->SetLineWidth(2);
 605
 606
 607           fFinalDPhi[i]->Fit(fTotal[i],flag);
 608           fFinalDPhi[i]->SetMinimum(0);
 609           fFinalDPhi[i]->DrawClone("");
 610           fFinalDPhi[i]->SetTitle(Form("0.7<p_{T,trig}<5.0 - 0.7<p_{T,assoc}<5.0 - %d-%d %",i*5,(i+1)*5));
 611           fFinalDPhi[i]->SetTitle("1/N_{trig} dN_{assoc}/d#Delta#varphi (rad^{-1})");
 612           fTotal[i]->DrawClone("same");
 613
 614           Double_t A1     = fTotal[i]->GetParameter(0);
 615           Double_t sigma1 = fTotal[i]->GetParameter(1);
 616           Double_t A2     = fTotal[i]->GetParameter(2);
 617           Double_t sigma2 = fTotal[i]->GetParameter(3);
 618           Double_t A3     = fTotal[i]->GetParameter(4);
 619           Double_t sigma3 = fTotal[i]->GetParameter(5);
 620
 621
 622           //define each gaussian and constant to be drawn with different colors on top of each other
 623
 624           TF1 * fConstant = new TF1("konst", "pol0(0)",-0.5*TMath::Pi(), 1.5*TMath::Pi());
 625           fConstant->SetParameter(0,(integr_for_const - TMath::Sqrt(TMath::Pi()*2)/binWidth*(A1*sigma1+A2*sigma2+A3*sigma3))/bins);
 626           fConstant->SetLineColor(kBlue);
 627           fConstant->Draw("same");
 628
 629           //gaus 1 NS
 630           TF1 * fGaussian1 = new TF1("fGaussian1", "[0]*exp(-x*x/(2*[1]*[1])) +[0] * exp(-(x-TMath::TwoPi())*(x-TMath::TwoPi())/(2*[1]*[1]))",-0.5*TMath::Pi(), 1.5*TMath::Pi());
 631           fGaussian1->SetParameters(fTotal[i]->GetParameter(0),fTotal[i]->GetParameter(1));
 632           fGaussian1->SetLineColor(kMagenta);
 633           fGaussian1->SetLineStyle(1);
 634           fGaussian1->Draw("same");
 635
 636           //gaus 2 NS
 637           TF1 * fGaussian2 = new TF1("fGaussian2", "[0]*exp(-x*x/(2*[1]*[1])) +[0] * exp(-(x-TMath::TwoPi())*(x-TMath::TwoPi())/(2*[1]*[1]))",-0.5*TMath::Pi(), 1.5*TMath::Pi());
 638           fGaussian2->SetLineColor(kGreen+2);
 639           fGaussian2->SetParameters(fTotal[i]->GetParameter(2),fTotal[i]->GetParameter(3));
 640           fGaussian2->Draw("same");
 641
 642           //gaus 3 AS
 643           TF1 * fGaussian3 = new TF1("fGaussian3", "[0] * exp(-((x-TMath::Pi()))*((x-TMath::Pi()))/(2*[1]*[1]))+[0] * exp(-((x+TMath::Pi()))*((x+TMath::Pi()))/(2*[1]*[1]))",-0.5*TMath::Pi(), 1.5*TMath::Pi());
 644           fGaussian3->SetLineColor(kCyan);
 645           fGaussian3->SetParameters(fTotal[i]->GetParameter(4), fTotal[i]->GetParameter(5));
 646           fGaussian3->Draw("same");
 647
 648
 649           Double_t a1e = fTotal[i]->GetParError(0);
 650           Double_t s1e = fTotal[i]->GetParError(1);
 651           Double_t a2e = fTotal[i]->GetParError(2);
 652           Double_t s2e = fTotal[i]->GetParError(3);
 653           Double_t a3e = fTotal[i]->GetParError(4);
 654           Double_t s3e = fTotal[i]->GetParError(5);
 655
 656           Double_t T1 = A1*sigma1;
 657           Double_t T2 = A2*sigma2;
 658           Double_t T3 = A3*sigma3;
 659           Double_t t1 = T1*TMath::Sqrt(a1e*a1e/A1/A1 + s1e*s1e/sigma1/sigma1);
 660           Double_t t2 = T2*TMath::Sqrt(a2e*a2e/A2/A2 + s2e*s2e/sigma2/sigma2);
 661           Double_t t3 = T3*TMath::Sqrt(a3e*a3e/A3/A3 + s3e*s3e/sigma3/sigma3);
 662
 663           //-
 664           Double_t nearSideResult = TMath::Sqrt(TMath::Pi()*2)/ binWidth* (A1 * sigma1 + A2 * sigma2)/trigger->GetBinContent(i+1);
 665           Double_t nearSideError = nearSideResult * TMath::Sqrt((t1*t1 + t2*t2)/(T1+T2)/(T1+T2)+ 1./trigger->GetBinContent(i+1));
 666           fNearSideIntegral->SetPoint(points,i, nearSideResult);
 667           fNearSideIntegral->SetPointError(points,0.5,nearSideError);
 668
 669           //-
 670
 671           //--
 672           Double_t awaySideResult = TMath::Sqrt(TMath::Pi()*2)/ binWidth*
 673         (A3 * sigma3)/trigger->GetBinContent(i+1);
 674           Double_t awaySideError = awaySideResult*TMath::Sqrt(a3e*a3e/A3/A3 + s3e*s3e/sigma3/sigma3 + 1/trigger->GetBinContent(i+1));
 675           fAwaySideIntegral->SetPoint(points,i, awaySideResult );
 676           fAwaySideIntegral->SetPointError(points,0.5, awaySideError );
 677           //--
 678
 679           //---
 680           bothSideResult = TMath::Sqrt(TMath::Pi()*2)/ binWidth* (A1 * sigma1 + A2 * sigma2 + A3 * sigma3 )/trigger->GetBinContent(i+1);
 681           bothSideError = nearSideResult *  TMath::Sqrt((t1*t1 + t2*t2 + t3*t3)/(T1+T2+T3)/(T1+T2+T3)+ 1./trigger->GetBinContent(i+1));
 682           fBothSideIntegral->SetPoint(points,i, bothSideResult );
 683           fBothSideIntegral->SetPointError(points,0.5, bothSideError );
 684           //---
 685
 686     }
 687     else{
 688
 689         fNearSideIntegral->SetPoint(points,i, 0);
 690         fAwaySideIntegral->SetPoint(points,i, 0);
 691         fBothSideIntegral->SetPoint(points,i,0);
 692
 693     }
 694     Double_t p0BinContent=p0->GetBinContent(i+1);
 695     Double_t p0BinError=p0->GetBinError(i+1);
 696
 697     //--------
 698     Double_t njets =  p0BinContent/(1+bothSideResult);
 699     Double_t njetsError = njets*TMath::Sqrt(bothSideError*bothSideError/(1+bothSideResult)/(1+bothSideResult) + p0BinError*p0BinError/p0BinContent/p0BinContent);
 700     fNjets->SetPoint(points,i, njets );
 701     fNjets->SetPointError(points,0.5,njetsError );
 702     //-------
 703
 704     fTriggerAverage->SetPoint(points,i, p0BinContent);
 705     fTriggerAverage->SetPointError(points,0.5, p0BinError);
 706
 707
 708   }
 709       points++;
 710   }
 711
 712
 713   TFile* file = 0x0;
 714   if(zyam)
 715       file = TFile::Open("njet_zyam.root","RECREATE");
 716   else
 717       file = TFile::Open("njet_fit.root","RECREATE");
 718
 719   fNearSideIntegral->Write();
 720   fAwaySideIntegral->Write();
 721   fBothSideIntegral->Write();
 722   fNjets->Write();
 723   fTriggerAverage->Write();
 724
 725   file->Close();
 726
 727
 728
 729 }
 730
 731